Golang la.PrintMat函数代码示例

func (o *Rjoint) debug_print_init() {
	sldNn := o.Sld.Cell.Shp.Nverts
	rodNn := o.Rod.Cell.Shp.Nverts
	rodNp := len(o.Rod.IpsElem)
	io.Pf("Nmat =\n")
	for i := 0; i < sldNn; i++ {
		for j := 0; j < rodNn; j++ {
			io.Pf("%g ", o.Nmat[i][j])
		}
		io.Pf("\n")
	}
	io.Pf("\nPmat =\n")
	for i := 0; i < sldNn; i++ {
		for j := 0; j < rodNp; j++ {
			io.Pf("%g ", o.Pmat[i][j])
		}
		io.Pf("\n")
	}
	io.Pf("\n")
	la.PrintMat("e0", o.e0, "%20.13f", false)
	io.Pf("\n")
	la.PrintMat("e1", o.e1, "%20.13f", false)
	io.Pf("\n")
	la.PrintMat("e2", o.e2, "%20.13f", false)
}

func Test_stat02(tst *testing.T) {

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

	x := [][]float64{
		{100, 100, 102, 98, 77, 99, 70, 105, 98},
		{80, 101, 12, 58, 47, 80, 20, 111, 89},
		{50, 130, 72, 38, 71, 15, 10, 12, 55},
	}

	y, z := StatTable(x, true, true)
	la.PrintMat("x", x, "%5g", false)
	la.PrintMat("y", y, "%13.6f", false)
	la.PrintMat("z", z, "%13.6f", false)
	io.Pforan("\nmin\n")
	chk.Scalar(tst, "y00=min(x[0,:])", 1e-17, y[0][0], 70)
	chk.Scalar(tst, "y01=min(x[1,:])", 1e-17, y[0][1], 12)
	chk.Scalar(tst, "y02=min(x[2,:])", 1e-17, y[0][2], 10)
	io.Pforan("\nave\n")
	chk.Scalar(tst, "y10=ave(x[0,:])", 1e-17, y[1][0], 849.0/9.0)
	chk.Scalar(tst, "y11=ave(x[1,:])", 1e-17, y[1][1], 598.0/9.0)
	chk.Scalar(tst, "y12=ave(x[2,:])", 1e-17, y[1][2], 453.0/9.0)
	io.Pforan("\nmax\n")
	chk.Scalar(tst, "y20=max(x[0,:])", 1e-17, y[2][0], 105)
	chk.Scalar(tst, "y21=max(x[1,:])", 1e-17, y[2][1], 111)
	chk.Scalar(tst, "y22=max(x[2,:])", 1e-17, y[2][2], 130)
	io.Pforan("\ndev\n")
	chk.Scalar(tst, "y30=dev(x[0,:])", 1e-17, y[3][0], 12.134661099511597)
	chk.Scalar(tst, "y31=dev(x[1,:])", 1e-17, y[3][1], 34.688294535444918)
	chk.Scalar(tst, "y32=dev(x[2,:])", 1e-17, y[3][2], 38.343839140075687)
	io.Pfyel("\nmin\n")
	chk.Scalar(tst, "z00=min(y[0,:])=min(min)", 1e-17, z[0][0], 10)
	chk.Scalar(tst, "z01=min(y[1,:])=min(ave)", 1e-17, z[0][1], 453.0/9.0)
	chk.Scalar(tst, "z02=min(y[2,:])=min(max)", 1e-17, z[0][2], 105)
	chk.Scalar(tst, "z03=min(y[3,:])=min(dev)", 1e-17, z[0][3], 12.134661099511597)
	io.Pfyel("\nave\n")
	chk.Scalar(tst, "z10=ave(y[0,:])=ave(min)", 1e-17, z[1][0], 92.0/3.0)
	chk.Scalar(tst, "z11=ave(y[1,:])=ave(ave)", 1e-17, z[1][1], ((849.0+598.0+453.0)/9.0)/3.0)
	chk.Scalar(tst, "z12=ave(y[2,:])=ave(max)", 1e-17, z[1][2], 346.0/3.0)
	chk.Scalar(tst, "z13=ave(y[3,:])=ave(dev)", 1e-17, z[1][3], (12.134661099511597+34.688294535444918+38.343839140075687)/3.0)
	io.Pfyel("\nmax\n")
	chk.Scalar(tst, "z20=max(y[0,:])=max(min)", 1e-17, z[2][0], 70)
	chk.Scalar(tst, "z21=max(y[1,:])=max(ave)", 1e-17, z[2][1], 849.0/9.0)
	chk.Scalar(tst, "z22=max(y[2,:])=max(max)", 1e-17, z[2][2], 130)
	chk.Scalar(tst, "z23=max(y[3,:])=max(dev)", 1e-17, z[2][3], 38.343839140075687)
	io.Pfyel("\ndev\n")
	chk.Scalar(tst, "z30=dev(y[0,:])=dev(min)", 1e-17, z[3][0], 34.078341117685483)
	chk.Scalar(tst, "z31=dev(y[1,:])=dev(ave)", 1e-17, z[3][1], 22.261169573539771)
	chk.Scalar(tst, "z32=dev(y[2,:])=dev(max)", 1e-17, z[3][2], 13.051181300301263)
	chk.Scalar(tst, "z33=dev(y[3,:])=dev(dev)", 1e-17, z[3][3], 14.194778389023206)
}

func main() {

	// input matrix in Triplet format
	// including repeated positions. e.g. (0,0)
	var A la.Triplet
	A.Init(5, 5, 13)
	A.Put(0, 0, 1.0) // << repeated
	A.Put(0, 0, 1.0) // << repeated
	A.Put(1, 0, 3.0)
	A.Put(0, 1, 3.0)
	A.Put(2, 1, -1.0)
	A.Put(4, 1, 4.0)
	A.Put(1, 2, 4.0)
	A.Put(2, 2, -3.0)
	A.Put(3, 2, 1.0)
	A.Put(4, 2, 2.0)
	A.Put(2, 3, 2.0)
	A.Put(1, 4, 6.0)
	A.Put(4, 4, 1.0)

	// right-hand-side
	b := []float64{8.0, 45.0, -3.0, 3.0, 19.0}

	// solve
	x, err := la.SolveRealLinSys(&A, b)
	if err != nil {
		io.Pfred("solver failed:\n%v", err)
		return
	}

	// output
	la.PrintMat("a", A.ToMatrix(nil).ToDense(), "%5g", false)
	la.PrintVec("b", b, "%v ", false)
	la.PrintVec("x", x, "%v ", false)
}

func (o *Rjoint) debug_print_K() {
	sldNn := o.Sld.Cell.Shp.Nverts
	rodNn := o.Rod.Cell.Shp.Nverts
	K := la.MatAlloc(o.Ny, o.Ny)
	start := o.Sld.Nu
	for i := 0; i < o.Ndim; i++ {
		for m := 0; m < sldNn; m++ {
			r := i + m*o.Ndim
			for j := 0; j < o.Ndim; j++ {
				for n := 0; n < sldNn; n++ {
					c := j + n*o.Ndim
					K[r][c] = o.Kss[r][c]
				}
				for n := 0; n < rodNn; n++ {
					c := j + n*o.Ndim
					K[r][start+c] = o.Ksr[r][c]
					K[start+c][r] = o.Krs[c][r]
				}
			}
		}
	}
	for i := 0; i < o.Ndim; i++ {
		for m := 0; m < rodNn; m++ {
			r := i + m*o.Ndim
			for j := 0; j < o.Ndim; j++ {
				for n := 0; n < rodNn; n++ {
					c := j + n*o.Ndim
					K[start+r][start+c] = o.Krr[r][c]
				}
			}
		}
	}
	la.PrintMat("K", K, "%20.10f", false)
}

func main() {

	// input matrix in Triplet format
	// including repeated positions. e.g. (0,0)
	var A la.Triplet
	A.Init(5, 5, 13)
	A.Put(0, 0, 1.0) // << repeated
	A.Put(0, 0, 1.0) // << repeated
	A.Put(1, 0, 3.0)
	A.Put(0, 1, 3.0)
	A.Put(2, 1, -1.0)
	A.Put(4, 1, 4.0)
	A.Put(1, 2, 4.0)
	A.Put(2, 2, -3.0)
	A.Put(3, 2, 1.0)
	A.Put(4, 2, 2.0)
	A.Put(2, 3, 2.0)
	A.Put(1, 4, 6.0)
	A.Put(4, 4, 1.0)

	// right-hand-side
	b := []float64{8.0, 45.0, -3.0, 3.0, 19.0}

	// allocate solver
	lis := la.GetSolver("umfpack")
	defer lis.Clean()

	// info
	symmetric := false
	verbose := false
	timing := false

	// initialise solver (R)eal
	err := lis.InitR(&A, symmetric, verbose, timing)
	if err != nil {
		io.Pfred("solver failed:\n%v", err)
		return
	}

	// factorise
	err = lis.Fact()
	if err != nil {
		io.Pfred("solver failed:\n%v", err)
		return
	}

	// solve (R)eal
	var dummy bool
	x := make([]float64, len(b))
	err = lis.SolveR(x, b, dummy) // x := inv(a) * b
	if err != nil {
		io.Pfred("solver failed:\n%v", err)
		return
	}

	// output
	la.PrintMat("a", A.ToMatrix(nil).ToDense(), "%5g", false)
	la.PrintVec("b", b, "%v ", false)
	la.PrintVec("x", x, "%v ", false)
}

func (o ElemUP) debug_print_K() {
	la.PrintMat("Kpp", o.P.Kpp, "%20.10f", false)
	la.PrintMat("Kpf", o.P.Kpf, "%20.10f", false)
	la.PrintMat("Kfp", o.P.Kfp, "%20.10f", false)
	la.PrintMat("Kff", o.P.Kff, "%20.10f", false)
	la.PrintMat("Kpu", o.Kpu, "%20.10f", false)
	la.PrintMat("Kup", o.Kup, "%20.10f", false)
	la.PrintMat("Kuu", o.U.K, "%20.10f", false)
}

func Test_nurbs02(tst *testing.T) {

	//verbose()
	chk.PrintTitle("nurbs02. square with initial stress. run")

	// fem
	analysis := NewFEM("data/nurbs02.sim", "", true, false, false, false, chk.Verbose, 0)

	// run simulation
	err := analysis.Run()
	if err != nil {
		tst.Errorf("Run failed\n%v", err)
		return
	}

	// domain
	dom := analysis.Domains[0]

	e := dom.Elems[0].(*ElemU)
	io.PfYel("fex = %v\n", e.fex)
	io.PfYel("fey = %v\n", e.fey)
	la.PrintMat("K", e.K, "%10.2f", false)

	// solution
	var sol ana.CteStressPstrain
	sol.Init(fun.Prms{
		&fun.Prm{N: "qnH0", V: -20},
		&fun.Prm{N: "qnV0", V: -20},
		&fun.Prm{N: "qnH", V: -50},
		&fun.Prm{N: "qnV", V: -100},
	})

	// check displacements
	t := dom.Sol.T
	tolu := 1e-16
	for _, n := range dom.Nodes {
		eqx := n.GetEq("ux")
		eqy := n.GetEq("uy")
		u := []float64{dom.Sol.Y[eqx], dom.Sol.Y[eqy]}
		io.Pfyel("u = %v\n", u)
		sol.CheckDispl(tst, t, u, n.Vert.C, tolu)
	}

	// check stresses
	tols := 1e-13
	for idx, ip := range e.IpsElem {
		x := e.Cell.Shp.IpRealCoords(e.X, ip)
		σ := e.States[idx].Sig
		io.Pforan("σ = %v\n", σ)
		sol.CheckStress(tst, t, σ, x, tols)
	}
}

func main() {

	mpi.Start(false)
	defer func() {
		mpi.Stop(false)
	}()

	if mpi.Rank() == 0 {
		chk.PrintTitle("Test SumToRoot 01")
	}

	M := [][]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},
	}

	id, sz, m := mpi.Rank(), mpi.Size(), len(M)
	start, endp1 := (id*m)/sz, ((id+1)*m)/sz

	if sz > 6 {
		chk.Panic("this test works with at most 6 processors")
	}

	var J la.Triplet
	J.Init(m, m, m*m)
	for i := start; i < endp1; i++ {
		for j := 0; j < m; j++ {
			J.Put(i, j, M[i][j])
		}
	}
	la.PrintMat(fmt.Sprintf("J @ proc # %d", id), J.ToMatrix(nil).ToDense(), "%10.1f", false)

	la.SpTriSumToRoot(&J)
	var tst testing.T
	if mpi.Rank() == 0 {
		chk.Matrix(&tst, "J @ proc 0", 1.0e-17, J.ToMatrix(nil).ToDense(), [][]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},
		})
	}
}

func Test_hyperelast02(tst *testing.T) {

	//verbose()
	chk.PrintTitle("hyperelast02 (linear)")

	E, ν := 1500.0, 0.25
	K := Calc_K_from_Enu(E, ν)
	G := Calc_G_from_Enu(E, ν)
	io.Pforan("K = %v\n", K)
	io.Pforan("G = %v\n", G)

	var m HyperElast1
	m.Init(2, false, []*fun.Prm{
		&fun.Prm{N: "K0", V: K},
		&fun.Prm{N: "G0", V: G},
		&fun.Prm{N: "le", V: 1},
	})
	io.Pforan("m = %+v\n", m)

	ε := []float64{-0.001, -0.002, -0.003}
	σ := make([]float64, 3)
	m.L_update(σ, ε)
	io.Pfblue2("ε = %v\n", ε)
	io.Pfcyan("σ = %v\n", σ)

	D := la.MatAlloc(3, 3)
	m.L_CalcD(D, ε)
	la.PrintMat("D", D, "%14.6f", false)

	tol := 1e-11
	verb := io.Verbose
	var tmp float64
	for i := 0; i < 3; i++ {
		for j := 0; j < 3; j++ {
			dnum := num.DerivCen(func(x float64, args ...interface{}) (res float64) {
				tmp, ε[j] = ε[j], x
				m.L_update(σ, ε)
				res = σ[i]
				ε[j] = tmp
				return
			}, ε[j])
			chk.AnaNum(tst, io.Sf("D%d%d", i, j), tol, D[i][j], dnum, verb)
		}
	}
}

func Test_nurbs03(tst *testing.T) {

	//verbose()
	chk.PrintTitle("nurbs03. ini stress free square")

	// fem
	analysis := NewFEM("data/nurbs03.sim", "", true, false, false, false, chk.Verbose, 0)

	// run simulation
	err := analysis.Run()
	if err != nil {
		tst.Errorf("Run failed\n%v", err)
		return
	}

	// domain
	dom := analysis.Domains[0]

	// element
	e := dom.Elems[0].(*ElemU)
	io.PfYel("fex = %v\n", e.fex)
	io.PfYel("fey = %v\n", e.fey)
	la.PrintMat("K", e.K, "%10.2f", false)

	// solution
	var sol ana.CteStressPstrain
	sol.Init(fun.Prms{
		&fun.Prm{N: "qnH", V: -50},
		&fun.Prm{N: "qnV", V: -100},
	})

	// check displacements
	t := dom.Sol.T
	tolu := 1e-16
	for _, n := range dom.Nodes {
		eqx := n.GetEq("ux")
		eqy := n.GetEq("uy")
		u := []float64{dom.Sol.Y[eqx], dom.Sol.Y[eqy]}
		io.Pfyel("u = %v\n", u)
		sol.CheckDispl(tst, t, u, n.Vert.C, tolu)
	}
}

func Test_hyperelast03(tst *testing.T) {

	//verbose()
	chk.PrintTitle("hyperelast03 (nonlinear)")

	var m HyperElast1
	m.Init(2, false, []*fun.Prm{
		&fun.Prm{N: "kap", V: 0.05},
		&fun.Prm{N: "kapb", V: 20.0},
		&fun.Prm{N: "G0", V: 1500},
		&fun.Prm{N: "pr", V: 2.2},
		&fun.Prm{N: "pt", V: 11.0},
	})
	io.Pforan("m = %+v\n", m)

	ε := []float64{-0.001, -0.002, -0.003}
	σ := make([]float64, 3)
	m.L_update(σ, ε)
	io.Pfblue2("ε = %v\n", ε)
	io.Pfcyan("σ = %v\n", σ)

	D := la.MatAlloc(3, 3)
	m.L_CalcD(D, ε)
	la.PrintMat("D", D, "%14.6f", false)

	tol := 1e-7
	verb := io.Verbose
	var tmp float64
	for i := 0; i < 3; i++ {
		for j := 0; j < 3; j++ {
			dnum := num.DerivCen(func(x float64, args ...interface{}) (res float64) {
				tmp, ε[j] = ε[j], x
				m.L_update(σ, ε)
				res = σ[i]
				ε[j] = tmp
				return
			}, ε[j])
			chk.AnaNum(tst, io.Sf("D%d%d", i, j), tol, D[i][j], dnum, verb)
		}
	}
}

func Test_invs02(tst *testing.T) {

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

	eps := [][]float64{
		{100 / 200.0, 150 / 200.0, 5 / 200.0},
		{150 / 200.0, 100 / 200.0, 10 / 200.0},
		{5 / 200.0, 10 / 200.0, 100 / 200.0},
	}
	ε := make([]float64, 6)
	e := make([]float64, 6)
	e_ := make([]float64, 6)
	Ten2Man(ε, eps)
	εv := M_εv(ε)
	εd := M_εd(ε)
	eno, εv_, εd_ := M_devε(e, ε)
	Lε := make([]float64, 3)
	err := M_EigenValsNum(Lε, ε)
	if err != nil {
		tst.Errorf("test failed: %v\n", err)
		return
	}
	Lεv, Lεd := L_strains(Lε)
	la.MatVecMul(e_, 1, Psd, ε)
	la.PrintMat("eps", eps, "%8g", false)
	io.Pf("ε   = %v\n", ε)
	io.Pf("e   = %v\n", e)
	io.Pf("e_  = %v\n", e_)
	io.Pf("eno = %v\n", eno)
	io.Pf("εv  = %v  Lεv=%v\n", εv, Lεv)
	io.Pf("εd  = %v  Lεd=%v\n", εd, Lεd)
	io.Pf("εd_ = %v\n", εd_)
	chk.Scalar(tst, "Lεv", 1e-17, Lεv, εv)
	chk.Scalar(tst, "Lεd", 1e-15, Lεd, εd)
	chk.Scalar(tst, "εv", 1e-17, εv, εv_)
	chk.Scalar(tst, "εv", 1e-17, εv, eps[0][0]+eps[1][1]+eps[2][2])
	chk.Scalar(tst, "εd", 1e-13, εd, εd_)
	chk.Vector(tst, "e", 1e-17, e, e_)
}

func Test_invs01(tst *testing.T) {

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

	sig := [][]float64{
		{100, 150, 5},
		{150, 100, 10},
		{5, 10, 100},
	}
	σ := make([]float64, 6)
	s := make([]float64, 6)
	s_ := make([]float64, 6)
	Ten2Man(σ, sig) // σ := sig
	p := M_p(σ)
	q := M_q(σ)
	θ := M_θ(σ)
	sno, p_, q_ := M_devσ(s, σ)
	p1, q1, θ1 := M_pqθ(σ)
	la.MatVecMul(s_, 1, Psd, σ)
	la.PrintMat("sig", sig, "%8g", false)
	io.Pf("σ   = %v\n", σ)
	io.Pf("s   = %v\n", s)
	io.Pf("s_  = %v\n", s_)
	io.Pf("sno = %v\n", sno)
	io.Pf("p   = %v\n", p)
	io.Pf("q   = %v\n", q)
	io.Pf("q_  = %v\n", q_)
	io.Pf("θ   = %v\n", θ)
	chk.Scalar(tst, "p", 1e-17, p, p_)
	chk.Scalar(tst, "p", 1e-17, p, -100)
	chk.Scalar(tst, "q", 1e-17, q, 260.52830940226056)
	chk.Scalar(tst, "q", 1e-13, q, q_)
	chk.Vector(tst, "s", 1e-17, s, s_)
	chk.Scalar(tst, "p1", 1e-17, p, p1)
	chk.Scalar(tst, "q1", 1e-13, q, q1)
	chk.Scalar(tst, "θ1", 1e-17, θ, θ1)
}

func TestJacobian03(tst *testing.T) {

	//verbose()
	chk.PrintTitle("TestJacobian 03")

	// grid
	var g fdm.Grid2D
	//g.Init(1.0, 1.0, 4, 4)
	g.Init(1.0, 1.0, 6, 6)
	//g.Init(1.0, 1.0, 11, 11)

	// equations numbering
	var e fdm.Equations
	peq := utl.IntUnique(g.L, g.R, g.B, g.T)
	e.Init(g.N, peq)

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

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

	// prescribed values
	U2 := make([]float64, e.N2)
	for _, eq := range g.L {
		U2[e.FR2[eq]] = 50.0
	}
	for _, eq := range g.R {
		U2[e.FR2[eq]] = 0.0
	}
	for _, eq := range g.B {
		U2[e.FR2[eq]] = 0.0
	}
	for _, eq := range g.T {
		U2[e.FR2[eq]] = 50.0
	}

	// functions
	k11 := K11.ToMatrix(nil)
	k12 := K12.ToMatrix(nil)
	ffcn := func(fU1, U1 []float64) error { // K11*U1 + K12*U2 - F1
		la.VecCopy(fU1, -1, F1)            // fU1 := (-F1)
		la.SpMatVecMulAdd(fU1, 1, k11, U1) // fU1 += K11*U1
		la.SpMatVecMulAdd(fU1, 1, k12, U2) // fU1 += K12*U2
		return nil
	}
	Jfcn := func(dfU1dU1 *la.Triplet, U1 []float64) error {
		fdm.Assemble(dfU1dU1, &K12, F1, nil, &g, &e)
		return nil
	}
	U1 := make([]float64, e.N1)
	CompareJac(tst, ffcn, Jfcn, U1, 0.0075)

	print_jac := false
	if print_jac {
		W1 := make([]float64, e.N1)
		fU1 := make([]float64, e.N1)
		ffcn(fU1, U1)
		var Jnum la.Triplet
		Jnum.Init(e.N1, e.N1, e.N1*e.N1)
		Jacobian(&Jnum, ffcn, U1, fU1, W1)
		la.PrintMat("K11 ", K11.ToMatrix(nil).ToDense(), "%g ", false)
		la.PrintMat("Jnum", Jnum.ToMatrix(nil).ToDense(), "%g ", false)
	}

	test_ffcn := false
	if test_ffcn {
		Uc := []float64{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 50.0, 25.0, 325.0 / 22.0, 100.0 / 11.0, 50.0 / 11.0,
			0.0, 50.0, 775.0 / 22.0, 25.0, 375.0 / 22.0, 100.0 / 11.0, 0.0, 50.0, 450.0 / 11.0, 725.0 / 22.0,
			25.0, 325.0 / 22.0, 0.0, 50.0, 500.0 / 11.0, 450.0 / 11.0, 775.0 / 22.0, 25.0, 0.0, 50.0, 50.0,
			50.0, 50.0, 50.0, 50.0,
		}
		for i := 0; i < e.N1; i++ {
			U1[i] = Uc[e.RF1[i]]
		}
		fU1 := make([]float64, e.N1)
		min, max := la.VecMinMax(fU1)
		io.Pf("min/max fU1 = %v\n", min, max)
	}
}

func TestDiffusion1D(tst *testing.T) {

	//verbose()
	chk.PrintTitle("Test Diffusion 1D (cooling)")

	// solution parameters
	silent := false
	fixstp := true
	//fixstp := false
	//method := "FwEuler"
	method := "BwEuler"
	//method := "Dopri5"
	//method := "Radau5"
	//numjac := true
	numjac := false
	rtol := 1e-4
	atol := rtol

	// timestep
	t0, tf, dt := 0.0, 0.2, 0.03

	// problem data
	kx := 1.0 // conductivity
	N := 6    // number of nodes
	//Nb   := N + 2             // augmented system dimension
	xmax := 1.0               // length
	dx := xmax / float64(N-1) // spatial step size
	dxx := dx * dx
	mol := []float64{kx / dxx, -2.0 * kx / dxx, kx / dxx}

	// function
	fcn := func(f []float64, t float64, y []float64, args ...interface{}) error {
		for i := 0; i < N; i++ {
			f[i] = 0
			if i == 0 || i == N-1 {
				continue // skip presc node
			}
			for p, j := range []int{i - 1, i, i + 1} {
				if j < 0 {
					j = i + 1
				} //  left boundary
				if j == N {
					j = i - 1
				} //  right boundary
				f[i] += mol[p] * y[j]
			}
		}
		//io.Pfgrey("y = %v\n", y)
		//io.Pfcyan("f = %v\n", f)
		return nil
	}

	// Jacobian
	jac := func(dfdy *la.Triplet, t float64, y []float64, args ...interface{}) error {
		//chk.Panic("jac is not available")
		if dfdy.Max() == 0 {
			//dfdy.Init(Nb, Nb, 3*N)
			dfdy.Init(N, N, 3*N)
		}
		dfdy.Start()
		for i := 0; i < N; i++ {
			if i == 0 || i == N-1 {
				dfdy.Put(i, i, 0.0)
				continue
			}
			for p, j := range []int{i - 1, i, i + 1} {
				if j < 0 {
					j = i + 1
				} //  left boundary
				if j == N {
					j = i - 1
				} //  right boundary
				dfdy.Put(i, j, mol[p])
			}
		}
		return nil
	}

	// initial values
	x := utl.LinSpace(0.0, xmax, N)
	y := make([]float64, N)
	//y := make([]float64, Nb)
	for i := 0; i < N; i++ {
		y[i] = 4.0*x[i] - 4.0*x[i]*x[i]
	}

	// debug
	f0 := make([]float64, N)
	//f0 := make([]float64, Nb)
	fcn(f0, 0, y)
	if false {
		io.Pforan("y0 = %v\n", y)
		io.Pforan("f0 = %v\n", f0)
		var J la.Triplet
		jac(&J, 0, y)
		la.PrintMat("J", J.ToMatrix(nil).ToDense(), "%8.3f", false)
	}
	//chk.Panic("stop")

	/*
//.........这里部分代码省略.........