Golang chk.Matrix函数代码示例

func Test_invs04(tst *testing.T) {

	//verbose()
	chk.PrintTitle("invs04")

	a := []float64{-10.0, -20.0, -30.0, 4.0 * SQ2, 5.0 * SQ2, 6.0 * SQ2}
	at := Alloc2()
	Man2Ten(at, a)
	io.Pf("a = %v\n", a)
	chk.Matrix(tst, "Man2Ten", 1e-17, at, [][]float64{
		{-10, 4, 6},
		{4, -20, 5},
		{6, 5, -30},
	})

	b := []float64{-88, -77, -55, -3 * SQ2}
	bt := Alloc2()
	Man2Ten(bt, b)
	io.Pf("b = %v\n", b)
	chk.Matrix(tst, "Man2Ten", 1e-17, bt, [][]float64{
		{-88, -3, 0},
		{-3, -77, 0},
		{0, 0, -55},
	})

	ver := chk.Verbose
	run_invs_tests(tst, a, ver)
	run_invs_tests(tst, b, ver)
}

func Test_conv06(tst *testing.T) {

	//verbose()
	chk.PrintTitle("conv06")

	a := [][]float64{
		{1, 2, 3},
		{4, 5, 6},
		{7, 8, 9},
		{10, 11, 12},
	}
	a_ := MatAlloc(4, 3)
	am := MatToColMaj(a)
	aa := ColMajToMatNew(am, 4, 3)
	ColMajToMat(a_, am)
	io.Pforan("a  = %v\n", a)
	io.Pforan("am = %v\n", am)
	io.Pforan("aa = %v\n", aa)
	chk.Vector(tst, "a => am", 1e-17, am, []float64{1, 4, 7, 10, 2, 5, 8, 11, 3, 6, 9, 12})
	chk.Matrix(tst, "am => a", 1e-17, aa, a)
	chk.Matrix(tst, "am => a", 1e-17, a_, a)

	b := [][]float64{
		{1, 2, 3, 4},
		{5, 6, 7, 8},
		{9, 0, -1, -2},
	}
	bm := MatToColMaj(b)
	bb := ColMajToMatNew(bm, 3, 4)
	io.Pforan("b  = %v\n", b)
	io.Pforan("bm = %v\n", bm)
	io.Pforan("bb = %v\n", bb)
	chk.Vector(tst, "b => bm", 1e-15, bm, []float64{1, 5, 9, 2, 6, 0, 3, 7, -1, 4, 8, -2})
	chk.Matrix(tst, "bm => b", 1e-15, bb, b)
}

func Test_matinvSmall01(tst *testing.T) {

	//verbose()
	chk.PrintTitle("matinvSmall01")

	noise := 0.0
	tol := 1.0e-16

	// 1 x 1 matrix
	res := MatAlloc(1, 1)
	det, err := MatInv(res, [][]float64{{2.0}}, tol)
	if err != nil {
		chk.Panic("%v", err.Error())
	}
	chk.Scalar(tst, "1 x 1 matrix: det ", tol, det, 2.0)
	chk.Matrix(tst, "1 x 1 matrix: ", tol, res, [][]float64{{0.5}})

	// matrix: inverse
	A := [][]float64{{1.0, 2.0}, {3.0, 2.0}}
	Aicorr := [][]float64{{-0.5, 0.5}, {0.75, -0.25 + noise}}
	Ai := MatAlloc(2, 2)
	detA, err := MatInv(Ai, A, tol)
	if err != nil {
		chk.Panic("%v", err.Error())
	}
	chk.Scalar(tst, "matrix: inv (det) ", tol, detA, -4.0+noise)
	chk.Matrix(tst, "matrix: inv (A)   ", tol, Ai, Aicorr)

	// using MatInvG
	Ai_ := MatAlloc(2, 2)
	err = MatInvG(Ai_, A, tol)
	if err != nil {
		chk.Panic("%v", err.Error())
	}
	chk.Matrix(tst, "matrix: inv with MatInvG", tol, Ai_, Aicorr)

	// another test
	B := [][]float64{{9.0, 3.0, 5.0}, {-6.0, -9.0, 7.0}, {-1.0, -8.0, 1.0}}
	Bicorr := [][]float64{
		{7.642276422764227e-02, -6.991869918699187e-02, 1.073170731707317e-01},
		{-1.626016260162601e-03, 2.276422764227642e-02, -1.512195121951219e-01},
		{6.341463414634146e-02, 1.121951219512195e-01, -1.024390243902439e-01 + noise},
	}
	Bi := MatAlloc(3, 3)
	detB, err := MatInv(Bi, B, tol)
	if err != nil {
		chk.Panic("%v", err.Error())
	}
	chk.Scalar(tst, "matrix: det ", tol, detB, 615.0+noise)
	chk.Matrix(tst, "matrix: inv ", tol, Bi, Bicorr)

	// using MatInvG
	Bi_ := MatAlloc(3, 3)
	err = MatInvG(Bi_, B, tol)
	if err != nil {
		chk.Panic("%v", err.Error())
	}
	chk.Matrix(tst, "matrix: inv with MatInvG", tol, Bi_, Bicorr)
}

func Test_mylab01(tst *testing.T) {

	//verbose()
	chk.PrintTitle("mylab01")

	I := make([]int, 5)
	IntFill(I, 666)
	J := IntVals(5, 666)
	Js := StrVals(5, "666")
	M := IntsAlloc(3, 4)
	N := DblsAlloc(3, 4)
	S := StrsAlloc(2, 3)
	A := IntRange(-1)
	a := IntRange2(0, 0)
	b := IntRange2(0, 1)
	c := IntRange2(0, 5)
	C := IntRange3(0, -5, -1)
	d := IntRange2(2, 5)
	D := IntRange2(-2, 5)
	e := IntAddScalar(D, 2)
	f := DblOnes(5)
	ff := DblVals(5, 666)
	g := []int{1, 2, 3, 4, 3, 4, 2, 1, 1, 2, 3, 4, 4, 2, 3, 7, 8, 3, 8, 3, 9, 0, 11, 23, 1, 2, 32, 12, 4, 32, 4, 11, 37}
	h := IntUnique(g)
	G := []int{1, 2, 3, 38, 3, 5, 3, 1, 2, 15, 38, 1, 11}
	H := IntUnique(D, C, G, []int{16, 39})
	X, Y := MeshGrid2D(3, 6, 10, 20, 4, 3)
	P := [][]int{
		{1, 2, 3, 4, 5},
		{-1, -2, -3, -4, -5},
		{6, 7, 8, 9, 10},
	}
	Pc := IntsClone(P)
	chk.Ints(tst, "I", I, []int{666, 666, 666, 666, 666})
	chk.Ints(tst, "J", J, []int{666, 666, 666, 666, 666})
	chk.Strings(tst, "Js", Js, []string{"666", "666", "666", "666", "666"})
	chk.Ints(tst, "A", A, []int{})
	chk.Ints(tst, "a", a, []int{})
	chk.Ints(tst, "b", b, []int{0})
	chk.Ints(tst, "c", c, []int{0, 1, 2, 3, 4})
	chk.Ints(tst, "C", C, []int{0, -1, -2, -3, -4})
	chk.Ints(tst, "d", d, []int{2, 3, 4})
	chk.Ints(tst, "D", D, []int{-2, -1, 0, 1, 2, 3, 4})
	chk.Ints(tst, "e", e, []int{0, 1, 2, 3, 4, 5, 6})
	chk.Vector(tst, "f", 1e-16, f, []float64{1, 1, 1, 1, 1})
	chk.Vector(tst, "ff", 1e-16, ff, []float64{666, 666, 666, 666, 666})
	chk.Ints(tst, "h", h, []int{0, 1, 2, 3, 4, 7, 8, 9, 11, 12, 23, 32, 37})
	chk.Ints(tst, "H", H, []int{-4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 11, 15, 16, 38, 39})
	chk.IntMat(tst, "M", M, [][]int{{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}})
	chk.Matrix(tst, "N", 1e-17, N, [][]float64{{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}})
	chk.Matrix(tst, "X", 1e-17, X, [][]float64{{3, 4, 5, 6}, {3, 4, 5, 6}, {3, 4, 5, 6}})
	chk.Matrix(tst, "Y", 1e-17, Y, [][]float64{{10, 10, 10, 10}, {15, 15, 15, 15}, {20, 20, 20, 20}})
	chk.StrMat(tst, "S", S, [][]string{{"", "", ""}, {"", "", ""}})
	chk.IntMat(tst, "Pc", Pc, P)
}

func Test_tsr02(tst *testing.T) {

	//verbose()
	chk.PrintTitle("tsr02")

	F := [][]float64{
		{2, 8.0 / 3.0, 0},
		{0, 2, 0},
		{0, 0.0, 1},
	}
	Fi := Alloc2()
	C := Alloc2()
	b := Alloc2()
	J, err := Inv(Fi, F)
	if err != nil {
		chk.Panic("%v", err)
	}
	RightCauchyGreenDef(C, F)
	LeftCauchyGreenDef(b, F)
	io.Pf("F = %v\n", F)
	io.Pf("C = %v\n", C)
	io.Pf("b = %v\n", b)
	chk.Scalar(tst, "J", 1.0e-17, J, 4.0)
	chk.Matrix(tst, "C", 1.0e-17, C, [][]float64{{36.0 / 9.0, 48.0 / 9.0, 0}, {48.0 / 9.0, 100.0 / 9.0, 0}, {0, 0, 1}})
	chk.Matrix(tst, "b", 1.0e-17, b, [][]float64{{100.0 / 9.0, 48.0 / 9.0, 0}, {48.0 / 9.0, 36.0 / 9.0, 0}, {0, 0, 1}})

	λ, μ := 2.0, 3.0
	σ := Alloc2()
	for i := 0; i < 3; i++ {
		for j := 0; j < 3; j++ {
			σ[i][j] = (λ/J)*math.Log(J)*It[i][j] + (μ/J)*(b[i][j]-It[i][j])
		}
	}

	P := Alloc2()
	S := Alloc2()
	σfromP := Alloc2()
	σfromS := Alloc2()
	CauchyToPK1(P, σ, F, Fi, J)
	CauchyToPK2(S, σ, F, Fi, J)
	PK1ToCauchy(σfromP, P, F, Fi, J)
	PK2ToCauchy(σfromS, S, F, Fi, J)

	io.Pf("σ = %v\n", σ)
	io.Pf("P = %v\n", P)
	io.Pf("S = %v\n", S)
	chk.Matrix(tst, "σfromP", 1.0e-17, σfromP, σ)
	chk.Matrix(tst, "σfromS", 1.0e-14, σfromS, σ)
}

func Test_conv03(tst *testing.T) {

	//verbose()
	chk.PrintTitle("conv03")

	var K, L, A Triplet
	K.Init(6, 6, 36+2*6) // 2*6 == number of nonzeros in A
	L.Init(6, 6, 36+2*6) // 2*6 == number of nonzeros in A
	for i := 0; i < 6; i++ {
		for j := 0; j < 6; j++ {
			K.Put(i, j, 1000)
			L.Put(i, j, 1000)
		}
	}
	A.Init(2, 3, 6)
	for i := 0; i < 2; i++ {
		for j := 0; j < 3; j++ {
			A.Put(i, j, float64(10*(i+1)+j+1))
		}
	}
	Am := A.ToMatrix(nil)
	Ad := Am.ToDense()
	if chk.Verbose {
		Kd := K.ToMatrix(nil).ToDense()
		Ld := L.ToMatrix(nil).ToDense()
		PrintMat("K", Kd, "%8g", false)
		PrintMat("L", Ld, "%8g", false)
		PrintMat("A", Ad, "%8g", false)
	}
	K.PutMatAndMatT(&A)
	L.PutCCMatAndMatT(Am)
	Kaug := K.ToMatrix(nil).ToDense()
	Laug := L.ToMatrix(nil).ToDense()
	if chk.Verbose {
		PrintMat("K augmented", Kaug, "%8g", false)
		PrintMat("L augmented", Laug, "%8g", false)
	}
	Cor := [][]float64{
		{1000, 1000, 1000, 1011, 1021, 1000},
		{1000, 1000, 1000, 1012, 1022, 1000},
		{1000, 1000, 1000, 1013, 1023, 1000},
		{1011, 1012, 1013, 1000, 1000, 1000},
		{1021, 1022, 1023, 1000, 1000, 1000},
		{1000, 1000, 1000, 1000, 1000, 1000},
	}
	chk.Matrix(tst, "Kaug", 1.0e-17, Kaug, Cor)
	chk.Matrix(tst, "Laug", 1.0e-17, Laug, Cor)
}

func Test_conv01(tst *testing.T) {

	//verbose()
	chk.PrintTitle("conv01")

	var t Triplet
	t.Init(3, 3, 10)
	t.Put(0, 0, 5.0)
	t.Put(0, 0, 5.0)
	t.Put(0, 1, 11.0)
	t.Put(0, 2, 12.0)
	t.Put(1, 0, 20.0)
	t.Put(1, 1, 21.0)
	t.Put(1, 2, 22.0)
	t.Put(2, 0, 30.0)
	t.Put(2, 1, 31.0)
	t.Put(2, 2, 32.0)
	a := t.ToMatrix(nil)
	ad := a.ToDense()
	if chk.Verbose {
		PrintMat("a", ad, "%5g", false)
	}
	chk.Matrix(tst, "a", 1e-17, ad, [][]float64{
		{10, 11, 12},
		{20, 21, 22},
		{30, 31, 32},
	})
}

func Test_conv02(tst *testing.T) {

	//verbose()
	chk.PrintTitle("conv02")

	var t Triplet
	t.Init(4, 3, 4*3+2)
	t.Put(0, 0, 1.0)
	t.Put(0, 1, 2.0)
	t.Put(0, 2, 3.0)
	t.Put(1, 0, 4.0)
	t.Put(1, 1, 5.0)
	t.Put(1, 2, 6.0)
	t.Put(2, 0, 7.0)
	t.Put(2, 1, 8.0)
	t.Put(2, 2, 9.0)
	t.Put(3, 0, 4.0)
	t.Put(3, 1, 11.0)
	t.Put(3, 2, 12.0)
	t.Put(3, 0, 3.0) // repeated
	t.Put(3, 0, 3.0) // repeated
	a := t.ToMatrix(nil)
	ad := a.ToDense()
	if chk.Verbose {
		PrintMat("a", ad, "%5g", false)
	}
	chk.Matrix(tst, "a", 1e-17, ad, [][]float64{
		{1, 2, 3},
		{4, 5, 6},
		{7, 8, 9},
		{10, 11, 12},
	})
}

func Test_mmMul01(tst *testing.T) {

	//verbose()
	chk.PrintTitle("mmMul01. MatrixMatrix multiplication")

	a := [][]float64{
		{1.0, 2.0, 3.0},
		{0.5, 0.75, 1.5},
	}
	b := [][]float64{
		{0.1, 0.5, 0.5, 0.75},
		{0.2, 2.0, 2.0, 2.0},
		{0.3, 0.5, 0.5, 0.5},
	}
	c := MatAlloc(2, 4)
	MatMul(c, 1, a, b) // c := 1*a*b
	io.Pf("a = %v\n", a)
	io.Pf("b = %v\n", b)
	io.Pf("c = %v\n", c)

	ccor := [][]float64{
		{1.4, 6.0, 6.0, 6.25},
		{0.65, 2.5, 2.5, 2.625},
	}
	chk.Matrix(tst, "c", 1.0e-15, c, ccor)
}

func Test_jacobi01(tst *testing.T) {

	//verbose()
	chk.PrintTitle("jacobi01")

	A := [][]float64{
		{1, 2, 3},
		{2, 3, 2},
		{3, 2, 2},
	}
	Q := MatAlloc(3, 3)
	v := make([]float64, 3)
	nit, err := Jacobi(Q, v, A)
	if err != nil {
		chk.Panic("Jacobi failed:\n%v", err)
	}
	io.Pforan("number of iterations = %v\n", nit)
	PrintMat("A", A, "%13.8f", false)
	PrintMat("Q", Q, "%13.8f", false)
	PrintVec("v", v, "%13.8f", false)

	chk.Matrix(tst, "Q", 1e-17, Q, [][]float64{
		{7.81993314738381295e-01, 5.26633230856907386e-01, 3.33382506832158143e-01},
		{-7.14394870018381645e-02, 6.07084171793832561e-01, -7.91419742017035133e-01},
		{-6.19179178753124115e-01, 5.95068272145819699e-01, 5.12358171676802088e-01},
	})
	chk.Vector(tst, "v", 1e-17, v, []float64{-1.55809924785903786e+00, 6.69537390404459476e+00, 8.62725343814443657e-01})
}

func Test_Ts(tst *testing.T) {

	//verbose()
	chk.PrintTitle("Ts")

	nd := test_nd
	for m := 0; m < len(test_nd)-3; m++ {
		//for m := 0; m < 3; m++ {
		A := test_AA[m]
		a := M_Alloc2(nd[m])
		Ten2Man(a, A)
		s := M_Dev(a)
		Ts := M_Alloc4(nd[m])
		M_Ts(Ts, s)
		s2 := M_Alloc2(nd[m])
		ds2ds := M_Alloc4(nd[m])
		M_Sq(s2, s)
		M_SqDeriv(ds2ds, s)
		Ts_ := M_Alloc4(nd[m])
		for i := 0; i < len(a); i++ {
			for j := 0; j < len(a); j++ {
				for k := 0; k < len(a); k++ {
					for l := 0; l < len(a); l++ {
						Ts_[i][j] += Psd[i][k] * ds2ds[k][l] * Psd[l][j]
					}
				}
			}
		}
		chk.Matrix(tst, "Ts", 1e-13, Ts, Ts_)
	}
}

func Test_tsr03(tst *testing.T) {

	//verbose()
	chk.PrintTitle("tsr03")

	a := [][]float64{
		{4.0, 1.0 / SQ2, 0},
		{1.0 / SQ2, 5.0, 0},
		{0, 0, 6.0},
	}
	am := make([]float64, 4)
	aa := Alloc2()
	Ten2Man(am, a)
	Man2Ten(aa, am)
	io.Pf("a  = %v\n", a)
	io.Pf("am = %v\n", am)
	io.Pf("aa = %v\n", aa)
	chk.Matrix(tst, "aa", 1.0e-15, aa, a)

	for i := 0; i < 3; i++ {
		for j := 0; j < 3; j++ {
			chk.Scalar(tst, fmt.Sprintf("am[%d][%d]", i, j), 1.0e-15, M2T(am, i, j), a[i][j])
		}
	}

	b := [][]float64{
		{4.0, 1.0 / SQ2, 3.0 / SQ2},
		{1.0 / SQ2, 5.0, 2.0 / SQ2},
		{3.0 / SQ2, 2.0 / SQ2, 6.0},
	}
	bm := make([]float64, 6)
	bb := Alloc2()
	Ten2Man(bm, b)
	Man2Ten(bb, bm)
	io.Pf("b  = %v\n", b)
	io.Pf("bm = %v\n", bm)
	io.Pf("bb = %v\n", bb)
	chk.Matrix(tst, "bb", 1.0e-15, bb, b)

	for i := 0; i < 3; i++ {
		for j := 0; j < 3; j++ {
			chk.Scalar(tst, fmt.Sprintf("bm[%d][%d]", i, j), 1.0e-15, M2T(bm, i, j), b[i][j])
		}
	}
}

func Test_imap(tst *testing.T) {

	//utl.Tsilent = false
	chk.PrintTitle("Test imap")

	for name, shape := range factory {
		gndim := shape.Gndim
		if gndim == 1 {
			continue
		}

		io.Pfyel("--------------------------------- %-6s---------------------------------\n", name)

		// check inverse mapping
		tol := 1e-14
		noise := 0.01
		if name == "tri10" {
			tol = 1e-14
		}
		if shape.FaceNvertsMax > 2 {
			noise = 0.0
		}
		nverts := shape.Nverts
		C := la.MatAlloc(gndim, nverts)
		s := []float64{rand.Float64(), rand.Float64(), rand.Float64()} // scale factors
		la.MatCopy(C, 1.0, shape.NatCoords)
		_ = tol
		io.Pf("nverts:%v\n", nverts)
		io.Pf("gndim:%v\n", gndim)
		for i := 0; i < gndim; i++ {
			for j := 0; j < nverts; j++ {
				C[i][j] *= s[i]
				C[i][j] += noise * rand.Float64() // noise
			}
		}

		r := make([]float64, 3)
		x := make([]float64, 3)
		R := la.MatAlloc(gndim, nverts)

		for j := 0; j < nverts; j++ {
			for i := 0; i < gndim; i++ {
				x[i] = C[i][j]
			}
			err := shape.InvMap(r, x, C)
			io.Pf("r:%v\n", r)
			_ = err
			for i := 0; i < gndim; i++ {
				R[i][j] = r[i]
			}
		}

		chk.Matrix(tst, "checking", tol, R, shape.NatCoords)

		io.PfGreen("OK\n")
	}
}

func Test_assemb01(tst *testing.T) {

	chk.PrintTitle("assemb01")

	// grid
	var g Grid2D
	g.Init(1.0, 1.0, 3, 3)

	// equations numbering
	var e Equations
	e.Init(g.N, []int{0, 3, 6})

	// K11 and K12
	var K11, K12 la.Triplet
	InitK11andK12(&K11, &K12, &e)

	// assembly
	F1 := make([]float64, e.N1)
	Assemble(&K11, &K12, F1, nil, &g, &e)

	// check
	K11d := K11.ToMatrix(nil).ToDense()
	K12d := K12.ToMatrix(nil).ToDense()
	K11c := [][]float64{
		{16.0, -4.0, -8.0, 0.0, 0.0, 0.0},
		{-8.0, 16.0, 0.0, -8.0, 0.0, 0.0},
		{-4.0, 0.0, 16.0, -4.0, -4.0, 0.0},
		{0.0, -4.0, -8.0, 16.0, 0.0, -4.0},
		{0.0, 0.0, -8.0, 0.0, 16.0, -4.0},
		{0.0, 0.0, 0.0, -8.0, -8.0, 16.0},
	}
	K12c := [][]float64{
		{-4.0, 0.0, 0.0},
		{0.0, 0.0, 0.0},
		{0.0, -4.0, 0.0},
		{0.0, 0.0, 0.0},
		{0.0, 0.0, -4.0},
		{0.0, 0.0, 0.0},
	}
	chk.Matrix(tst, "K11: ", 1e-16, K11d, K11c)
	chk.Matrix(tst, "K12: ", 1e-16, K12d, K12c)
}

func RunInvCheck(tst *testing.T, key string, M, CorrectInvM [][]float64, checkI bool, Tol, TolI float64) {
	m, n := len(M), len(M[0])
	Mi := MatAlloc(n, m)
	t0 := time.Now()
	err := MatInvG(Mi, M, 1e-13)
	if err != nil {
		chk.Panic("%v", err.Error())
	}
	io.Pfpink("Lapack: time elapsed = %v\n", time.Now().Sub(t0))
	MMi := MatAlloc(m, m)
	MMiM := MatAlloc(m, n)
	MatMul(MMi, 1, M, Mi)   // MMi = M * Mi
	MatMul(MMiM, 1, MMi, M) // MMiM = M * Mi * M == M
	chk.Matrix(tst, io.Sf("%s => Mi - CorrectInvM ", key), Tol, Mi, CorrectInvM)
	chk.Matrix(tst, io.Sf("%s => M*Mi*M = M       ", key), Tol, MMiM, M)
	if checkI {
		I := MatAlloc(m, m)
		MatSetDiag(I, 1)
		chk.Matrix(tst, io.Sf("%s => M*Mi = I         ", key), TolI, MMi, I)
	}
}