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 {
		io.PfYel("\nTest MPI 03\n")
	}
	if mpi.Size() != 3 {
		chk.Panic("this test needs 3 processors")
	}
	x := []int{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}
	n := len(x)
	id, sz := mpi.Rank(), mpi.Size()
	start, endp1 := (id*n)/sz, ((id+1)*n)/sz
	for i := start; i < endp1; i++ {
		x[i] = i
	}

	//io.Pforan("x = %v\n", x)

	// IntAllReduceMax
	w := make([]int, n)
	mpi.IntAllReduceMax(x, w)
	var tst testing.T
	chk.Ints(&tst, fmt.Sprintf("IntAllReduceMax: x @ proc # %d", id), x, []int{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10})

	//io.Pfred("x = %v\n", x)
}

// CheckFront0 returns front0 and number of failed/success
func CheckFront0(opt *Optimiser, verbose bool) (nfailed int, front0 []*Solution) {
	front0 = make([]*Solution, 0)
	var nsuccess int
	for _, sol := range opt.Solutions {
		var failed bool
		for _, oor := range sol.Oor {
			if oor > 0 {
				failed = true
				break
			}
		}
		if failed {
			nfailed++
		} else {
			nsuccess++
			if sol.FrontId == 0 {
				front0 = append(front0, sol)
			}
		}
	}
	if verbose {
		if nfailed > 0 {
			io.PfRed("N failed = %d out of %d\n", nfailed, opt.Nsol)
		} else {
			io.PfGreen("N success = %d out of %d\n", nsuccess, opt.Nsol)
		}
		io.PfYel("N front 0 = %d\n", len(front0))
	}
	return
}

// skip skips test based on it and/or t
func (o testKb) skip() bool {
	if o.itmin >= 0 {
		if o.it < o.itmin {
			return true // skip
		}
	}
	if o.itmax >= 0 {
		if o.it > o.itmax {
			return true // skip
		}
	}
	if o.tmin >= 0 {
		if o.t < o.tmin {
			return true // skip
		}
	}
	if o.tmax >= 0 {
		if o.t > o.tmax {
			return true // skip
		}
	}
	if o.verb {
		io.PfYel("\nit=%2d t=%v\n", o.it, o.t)
	}
	return false
}

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)
	}
}

// PrintMemStat prints memory statistics
func PrintMemStat(msg string) {
	var mem runtime.MemStats
	runtime.ReadMemStats(&mem)
	io.PfYel("%s\n", msg)
	io.Pfyel("Alloc      = %v [KB]  %v [MB]  %v [GB]\n", mem.Alloc/KBSIZE, mem.Alloc/MBSIZE, mem.Alloc/GBSIZE)
	io.Pfyel("HeapAlloc  = %v [KB]  %v [MB]  %v [GB]\n", mem.HeapAlloc/KBSIZE, mem.HeapAlloc/MBSIZE, mem.HeapAlloc/GBSIZE)
	io.Pfyel("Sys        = %v [KB]  %v [MB]  %v [GB]\n", mem.Sys/KBSIZE, mem.Sys/MBSIZE, mem.Sys/GBSIZE)
	io.Pfyel("HeapSys    = %v [KB]  %v [MB]  %v [GB]\n", mem.HeapSys/KBSIZE, mem.HeapSys/MBSIZE, mem.HeapSys/GBSIZE)
	io.Pfyel("TotalAlloc = %v [KB]  %v [MB]  %v [GB]\n", mem.TotalAlloc/KBSIZE, mem.TotalAlloc/MBSIZE, mem.TotalAlloc/GBSIZE)
	io.Pfyel("Mallocs    = %v\n", mem.Mallocs)
	io.Pfyel("Frees      = %v\n", mem.Frees)
}

func main() {

	// finalise analysis process and catch errors
	defer out.End()

	// input data
	simfn := "d2-simple-flux"
	flag.Parse()
	if len(flag.Args()) > 0 {
		simfn = flag.Arg(0)
	}
	if io.FnExt(simfn) == "" {
		simfn += ".sim"
	}

	// start analysis process
	out.Extrap = []string{"nwlx", "nwly"}
	out.Start(simfn, 0, 0)

	// define entities
	out.Define("top-middle", out.At{5, 3})
	out.Define("section-A", out.N{-1})
	out.Define("section-B", out.Along{{0, 0}, {10, 0}})

	// load results
	out.LoadResults(nil)

	// compute water discharge along section-A
	nwlx_TM := out.GetRes("ex_nwlx", "top-middle", -1)
	Q := out.Integrate("ex_nwlx", "section-A", "y", -1)
	io.PfYel("Q = %g m³/s [answer: 0.0003]\n", Q)

	// plot
	kt := len(out.Times) - 1
	out.Splot("")
	out.Plot("pl", "y", "section-A", plt.Fmt{L: "t=0"}, 0)
	out.Plot("pl", "y", "section-A", plt.Fmt{L: io.Sf("t=%g", out.Times[kt])}, kt)
	out.Splot("")
	out.Plot("x", "pl", "section-B", plt.Fmt{L: "t=0"}, 0)
	out.Plot("x", "pl", "section-B", plt.Fmt{L: io.Sf("t=%g", out.Times[kt])}, kt)
	out.Splot("")
	out.Plot("t", nwlx_TM, "top-middle", plt.Fmt{}, -1)
	out.Csplot.Ylbl = "$n_{\\ell}\\cdot w_{\\ell x}$"

	// show
	if true {
		out.Draw("", "", true, func(i, j, nplots int) {
			if i == 2 && j == 1 {
				plt.Plot([]float64{0, 10}, []float64{10, 9}, "'k--'")
			}
		})
	}
}

func (o *PrincStrainsUp) dbg_plot(eid, ipid int, time float64) {
	if eid == o.DbgEid && ipid == o.DbgIpId {

		// skip if not selected time
		if o.DbgTime > -1 {
			if math.Abs(time-o.DbgTime) > 1e-7 {
				return
			}
		}

		// check Jacobian
		var cnd float64
		io.PfYel("\nchecking Jacobian: eid=%d ipid=%d\n", eid, ipid)
		cnd, err := o.nls.CheckJ(o.x, o.ChkJacTol, true, o.ChkSilent)
		if err != nil {
			io.PfRed("CheckJ failed:\n%v\n", err)
		}
		io.PfYel("after: cnd(J) = %v\n\n", cnd)
		io.Pfcyan("\nN  = %v\n", o.N)
		io.Pfcyan("Nb = %v\n", o.Nb)
		io.Pfcyan("A  = %v\n", o.A)
		io.Pfcyan("h  = %v\n", o.h)
		io.Pfcyan("Mb = %v\n", o.Mb)
		io.Pfcyan("a  = %v\n", o.a)
		io.Pfcyan("b  = %v\n", o.b)
		io.Pfcyan("c  = %v\n", o.c)
		io.Pfcyan("Lσ = %v\n", o.Lσ)
		io.Pforan("Ne = %v\n", o.Ne)
		io.Pforan("Mbe = %v\n", o.Mbe)
		io.Pforan("Fcoef = %v\n\n", o.Fcoef)

		// plot
		var plr Plotter
		plr.SetFig(false, false, 1.5, 400, "/tmp", io.Sf("fig_stress_eid%d_ipid%d", eid, ipid))
		plr.SetModel(o.Mdl)
		plr.PreCor = o.DbgPco
		plr.Plot([]string{"p,q,ys", "oct,ys"}, o.DbgRes, o.DbgSts, true, true)
	}
}

// Stat prints statistical analysis
func (o *SimpleFltProb) Stat(idxF, hlen int, Fref float64) {
	if o.C.Ntrials < 2 || o.Nfeasible < 2 {
		return
	}
	F := make([]float64, o.Nfeasible)
	for i := 0; i < o.Nfeasible; i++ {
		o.Fcn(o.ff[0], o.gg[0], o.hh[0], o.Xbest[i])
		F[i] = o.ff[0][idxF]
	}
	fmin, fave, fmax, fdev := rnd.StatBasic(F, true)
	io.Pf("fmin = %v\n", fmin)
	io.PfYel("fave = %v (%v)\n", fave, Fref)
	io.Pf("fmax = %v\n", fmax)
	io.Pf("fdev = %v\n\n", fdev)
	io.Pf(rnd.BuildTextHist(nice_num(fmin-0.05), nice_num(fmax+0.05), 11, F, "%.2f", hlen))
}

func main() {

	// filename
	filename, fnkey := io.ArgToFilename(0, "d2-simple-flux", ".sim", true)

	// start analysis process
	out.Extrap = []string{"nwlx", "nwly"}
	out.Start(filename, 0, 0)

	// define entities
	out.Define("top-middle", out.At{5, 3})
	out.Define("section-A", out.N{-1})
	out.Define("section-B", out.Along{{0, 0}, {10, 0}})

	// load results
	out.LoadResults(nil)

	// compute water discharge along section-A
	nwlx_TM := out.GetRes("ex_nwlx", "top-middle", -1)
	Q := out.Integrate("ex_nwlx", "section-A", "y", -1)
	io.PfYel("Q = %g m³/s [answer: 0.0003]\n", Q)

	// plot
	kt := len(out.Times) - 1
	out.Splot("")
	out.Plot("pl", "y", "section-A", plt.Fmt{L: "t=0"}, 0)
	out.Plot("pl", "y", "section-A", plt.Fmt{L: io.Sf("t=%g", out.Times[kt])}, kt)
	out.Splot("")
	out.Plot("x", "pl", "section-B", plt.Fmt{L: "t=0"}, 0)
	out.Plot("x", "pl", "section-B", plt.Fmt{L: io.Sf("t=%g", out.Times[kt])}, kt)
	out.Splot("")
	out.Plot("t", nwlx_TM, "top-middle", plt.Fmt{}, -1)
	out.Csplot.Ylbl = "$n_{\\ell}\\cdot w_{\\ell x}$"

	// save
	plt.SetForPng(1.5, 400, 200)
	out.Draw("/tmp", "seep_simple_flux_"+fnkey+".png", false, func(i, j, nplots int) {
		if i == 2 && j == 1 {
			plt.Plot([]float64{0, 10}, []float64{10, 9}, "'k--'")
		}
	})
}

func main() {

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

	if mpi.Rank() == 0 {
		io.PfYel("\nTest MPI 04\n")
	}

	for i := 0; i < 60; i++ {
		time.Sleep(1e9)
		io.Pf("hello from %v\n", mpi.Rank())
		if mpi.Rank() == 2 && i == 3 {
			io.PfGreen("rank = 3 wants to abort (the following error is OK)\n")
			mpi.Abort()
		}
	}
}

// do_solve solve nonlinear problem
func (o *PrincStrainsUp) do_solve(bsmp float64, eid, ipid int, time float64) (err error) {
	silent := true
	if o.DbgOn {
		silent = o.dbg_silent(eid, ipid)
		if !silent {
			io.PfYel("\n\neid=%d ipid=%d time=%g: running with bsmp=%g\n", eid, ipid, time, bsmp)
		}
	}
	err = o.nls.Solve(o.x, silent)
	if err != nil {
		if o.DbgOn {
			o.dbg_plot(eid, ipid, time)
		}
		return
	}
	if o.DbgOn && o.DbgPlot {
		o.dbg_plot(eid, ipid, time)
	}
	return
}

func Test_nurbs01(tst *testing.T) {

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

	b := get_nurbs_A()
	elems := b.Elements()
	enodes := b.Enodes()
	io.PfYel("enodes = %v\n", enodes)
	chk.Ints(tst, "elem[0]", elems[0], []int{2, 3, 1, 2})
	chk.Ints(tst, "elem[1]", elems[1], []int{3, 4, 1, 2})
	chk.Ints(tst, "elem[2]", elems[2], []int{4, 5, 1, 2})
	chk.Ints(tst, "enodes[0]", enodes[0], []int{0, 1, 2, 5, 6, 7})
	chk.Ints(tst, "enodes[1]", enodes[1], []int{1, 2, 3, 6, 7, 8})
	chk.Ints(tst, "enodes[2]", enodes[2], []int{2, 3, 4, 7, 8, 9})

	if T_NURBS_SAVE {
		do_plot_nurbs_basis(b, 0, 7)
		plt.SaveD("/tmp/gosl", "t_nurbs01.eps")
	}
}

func Test_eigenp01(tst *testing.T) {

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

	// constants
	tolP := 1e-14        // eigenprojectors
	tolS := 1e-13        // spectral decomposition
	toldP := 1e-9        // derivatives of eigenprojectors
	ver := chk.Verbose   // check P verbose
	verdP := chk.Verbose // check dPda verbose

	// run test
	nd := test_nd
	for idxA := 0; idxA < len(test_nd); idxA++ {
		//for idxA := 10; idxA < 11; idxA++ {
		//for idxA := 11; idxA < 12; idxA++ {
		//for idxA := 12; idxA < 13; idxA++ {

		// tensor and eigenvalues
		A := test_AA[idxA]
		a := M_Alloc2(nd[idxA])
		Ten2Man(a, A)
		io.PfYel("\n\ntst # %d ###################################################################################\n", idxA)
		io.Pfblue2("a = %v\n", a)
		io.Pfblue2("λ = %v\n", test_λ[idxA])

		// check eigenprojectors
		io.Pforan("\neigenprojectors\n")
		λsorted := CheckEigenprojs(a, tolP, tolS, ver)
		io.Pfyel("λsorted = %v\n", λsorted)
		λchk := utl.DblGetSorted(test_λ[idxA])
		chk.Vector(tst, "λchk", 1e-12, λsorted, λchk)

		// check derivatives of eigenprojectors
		io.Pforan("\nderivatives\n")
		CheckEigenprojsDerivs(a, toldP, verdP, EV_ZERO)

	}
}

// main function
func main() {

	// flags
	benchmark := false
	ncpuMax := 16

	// benchmarking
	if benchmark {
		var nsol, tf int
		var et time.Duration
		X := make([]float64, ncpuMax)
		T := make([]float64, ncpuMax)
		S := make([]float64, ncpuMax) // speedup
		S[0] = 1
		for i := 0; i < ncpuMax; i++ {
			io.Pf("\n\n")
			nsol, tf, et = runone(i + 1)
			io.PfYel("elaspsedTime = %v\n", et)
			X[i] = float64(i + 1)
			T[i] = et.Seconds()
			if i > 0 {
				S[i] = T[0] / T[i] // Told / Tnew
			}
		}

		plt.SetForEps(0.75, 250)
		plt.Plot(X, S, io.Sf("'b-',marker='.', label='speedup: $N_{sol}=%d,\\,t_f=%d$', clip_on=0, zorder=100", nsol, tf))
		plt.Plot([]float64{1, 16}, []float64{1, 16}, "'k--',zorder=50")
		plt.Gll("$N_{cpu}:\\;$ number of groups", "speedup", "leg_out=1")
		plt.DoubleYscale("$T_{sys}:\\;$ system time [s]")
		plt.Plot(X, T, "'k-',color='gray', clip_on=0")
		plt.SaveD("/tmp/goga", "topology-speedup.eps")
		return
	}

	// normal run
	runone(-1)
}

//.........这里部分代码省略.........
		sin2, cos2 := math.Sin(θ2), math.Cos(θ2)
		fcn = func(f, g, h, x []float64, ξ []int, cpu int) {
			c0 := 1.0
			for i := 1; i < len(x); i++ {
				c0 += x[i]
			}
			f[0] = x[0]
			f[1] = c0 * (1.0 - f[0]/c0)
			if true {
				c1 := cos1*(f[1]-e) - sin1*f[0]
				c2 := sin1*(f[1]-e) + cos1*f[0]
				c3 := math.Sin(b * PI * math.Pow(c2, c))
				g[0] = c1 - a*math.Pow(math.Abs(c3), d)
				d1 := cos2*(f[1]-E) - sin2*f[0]
				d2 := sin2*(f[1]-E) + cos2*f[0]
				d3 := math.Sin(B * PI * math.Pow(d2, C))
				g[1] = d1 - A*math.Pow(math.Abs(d3), D)
			}
		}
		extraplot = func() {
			np := 401
			X, Y := utl.MeshGrid2D(0, 1, 0, 20, np, np)
			Z1 := utl.DblsAlloc(np, np)
			Z2 := utl.DblsAlloc(np, np)
			Z3 := utl.DblsAlloc(np, np)
			for j := 0; j < np; j++ {
				for i := 0; i < np; i++ {
					c1 := cos1*(Y[i][j]-e) - sin1*X[i][j]
					c2 := sin1*(Y[i][j]-e) + cos1*X[i][j]
					c3 := math.Sin(b * PI * math.Pow(c2, c))
					d1 := cos2*(Y[i][j]-E) - sin2*X[i][j]
					d2 := sin2*(Y[i][j]-E) + cos2*X[i][j]
					d3 := math.Sin(B * PI * math.Pow(d2, C))
					Z1[i][j] = c1 - a*math.Pow(math.Abs(c3), d)
					Z2[i][j] = d1 - A*math.Pow(math.Abs(d3), D)
					if Z1[i][j] >= 0 && Z2[i][j] >= 0 {
						Z3[i][j] = 1
					} else {
						Z3[i][j] = -1
					}
				}
			}
			plt.Contour(X, Y, Z3, "colors=['white','gray'],clabels=0,cbar=0,lwd=0.5,fsz=5")
			plt.ContourSimple(X, Y, Z1, false, 7, "linestyles=['--'], linewidths=[0.7], colors=['gray'], levels=[0]")
			plt.ContourSimple(X, Y, Z2, false, 7, "linestyles=['--'], linewidths=[0.7], colors=['gray'], levels=[0]")
		}
		opt.Multi_fcnErr = CTPerror1(θ1, a, b, c, d, e)

	default:
		chk.Panic("problem %d is not available", problem)
	}

	// initialise optimiser
	opt.Init(goga.GenTrialSolutions, nil, fcn, nf, ng, nh)

	// initial solutions
	var sols0 []*goga.Solution
	if false {
		sols0 = opt.GetSolutionsCopy()
	}

	// solve
	opt.RunMany("", "")
	goga.StatMulti(opt, true)
	io.PfYel("Tsys = %v\n", opt.SysTime)

	// check
	goga.CheckFront0(opt, true)

	// plot
	if true {
		feasibleOnly := false
		plt.SetForEps(0.8, 300)
		fmtAll := &plt.Fmt{L: "final solutions", M: ".", C: "orange", Ls: "none", Ms: 3}
		fmtFront := &plt.Fmt{L: "final Pareto front", C: "r", M: "o", Ms: 3, Ls: "none"}
		goga.PlotOvaOvaPareto(opt, sols0, 0, 1, feasibleOnly, fmtAll, fmtFront)
		extraplot()
		//plt.AxisYrange(0, f1max)
		if problem > 0 && problem < 6 {
			plt.Text(0.05, 0.05, "unfeasible", "color='gray', ha='left',va='bottom'")
			plt.Text(0.95, f1max-0.05, "feasible", "color='white', ha='right',va='top'")
		}
		if opt.RptName == "CTP6" {
			plt.Text(0.02, 0.15, "unfeasible", "rotation=-7,color='gray', ha='left',va='bottom'")
			plt.Text(0.02, 6.50, "unfeasible", "rotation=-7,color='gray', ha='left',va='bottom'")
			plt.Text(0.02, 13.0, "unfeasible", "rotation=-7,color='gray', ha='left',va='bottom'")
			plt.Text(0.50, 2.40, "feasible", "rotation=-7,color='white', ha='center',va='bottom'")
			plt.Text(0.50, 8.80, "feasible", "rotation=-7,color='white', ha='center',va='bottom'")
			plt.Text(0.50, 15.30, "feasible", "rotation=-7,color='white', ha='center',va='bottom'")
		}
		if opt.RptName == "TNK" {
			plt.Text(0.05, 0.05, "unfeasible", "color='gray', ha='left',va='bottom'")
			plt.Text(0.80, 0.85, "feasible", "color='white', ha='left',va='top'")
			plt.Equal()
			plt.AxisRange(0, 1.22, 0, 1.22)
		}
		plt.SaveD("/tmp/goga", io.Sf("%s.eps", opt.RptName))
	}
	return
}

func Test_ops03(tst *testing.T) {

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

	nonsymTol := 1e-15

	dtol := 1e-9
	dver := chk.Verbose

	nd := test_nd
	for idxA := 0; idxA < len(test_nd)-3; idxA++ {
		//for idxA := 0; idxA < 1; idxA++ {

		// tensor and eigenvalues
		A := test_AA[idxA]
		a := M_Alloc2(nd[idxA])
		Ten2Man(a, A)
		io.PfYel("\n\ntst # %d ###################################################################################\n", idxA)
		io.Pfblue2("a = %v\n", a)

		// inverse
		Ai := Alloc2()
		ai := M_Alloc2(nd[idxA])
		detA, err := Inv(Ai, A)
		if err != nil {
			chk.Panic("%v", err)
		}
		deta_ := M_Det(a)
		deta, err := M_Inv(ai, a, MINDET)
		if err != nil {
			chk.Panic("%v", err)
		}
		Ai_ := Alloc2()
		Man2Ten(Ai_, ai)
		aia := M_Alloc2(nd[idxA])
		err = M_Dot(aia, ai, a, nonsymTol)
		if err != nil {
			chk.Panic("%v", err)
		}
		chk.Scalar(tst, "detA", 1e-14, detA, deta)
		chk.Scalar(tst, "deta", 1e-14, deta, deta_)
		chk.Matrix(tst, "Ai", 1e-14, Ai, Ai_)
		chk.Vector(tst, "ai*a", 1e-15, aia, Im[:2*nd[idxA]])
		io.Pforan("ai*a = %v\n", aia)

		// derivative of inverse
		dtol_tmp := dtol
		if idxA == 5 {
			dtol = 1e-8
		}
		var tmp float64
		ai_tmp := M_Alloc2(nd[idxA])
		daida := M_Alloc4(nd[idxA])
		M_InvDeriv(daida, ai)
		io.Pforan("ai = %v\n", ai)
		for i := 0; i < len(a); i++ {
			for j := 0; j < len(a); j++ {
				//dnum, _ := num.DerivForward(func(x float64, args ...interface{}) (res float64) {
				dnum, _ := num.DerivCentral(func(x float64, args ...interface{}) (res float64) {
					tmp, a[j] = a[j], x
					_, err := M_Inv(ai_tmp, a, MINDET)
					a[j] = tmp
					if err != nil {
						chk.Panic("daida failed:\n%v", err)
					}
					return ai_tmp[i]
				}, a[j], 1e-6)
				chk.AnaNum(tst, io.Sf("dai/da[%d][%d]", i, j), dtol, daida[i][j], dnum, dver)
			}
		}
		dtol = dtol_tmp
	}
}

func Test_ops02(tst *testing.T) {

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

	nd := []int{2, 2, 3, 3, 3}
	AA := [][][]float64{
		{
			{1, 2, 0},
			{2, -2, 0},
			{0, 0, -2},
		},
		{
			{-100, 33, 0},
			{33, -200, 0},
			{0, 0, 150},
		},
		{
			{1, 2, 4},
			{2, -2, 3},
			{4, 3, -2},
		},
		{
			{-100, -10, 20},
			{-10, -200, 15},
			{20, 15, -300},
		},
		{
			{-100, 0, -10},
			{0, -200, 0},
			{-10, 0, 100},
		},
	}
	BB := [][][]float64{
		{
			{0.13, 1.2, 0},
			{1.2, -20, 0},
			{0, 0, -28},
		},
		{
			{-10, 3.3, 0},
			{3.3, -2, 0},
			{0, 0, 1.5},
		},
		{
			{0.1, 0.2, 0.8},
			{0.2, -1.3, 0.3},
			{0.8, 0.3, -0.2},
		},
		{
			{-10, -1, 2},
			{-1, -20, 1},
			{2, 1, -30},
		},
		{
			{-10, 3, -1},
			{3, -20, 1},
			{-1, 1, 10},
		},
	}

	nonsymTol := 1e-15

	for m := 0; m < len(nd); m++ {

		// tensors
		A := AA[m]
		B := BB[m]
		a := M_Alloc2(nd[m])
		b := M_Alloc2(nd[m])
		Ten2Man(a, A)
		Ten2Man(b, B)
		io.PfYel("\n\ntst # %d ###################################################################################\n", m)
		io.Pfblue2("a = %v\n", a)
		io.Pfblue2("b = %v\n", b)

		// dyadic
		c := M_Dy(a, b)
		c_ := M_Alloc4(nd[m])
		c__ := M_Alloc4(nd[m])
		M_DyAdd(c_, 1, a, b)
		for i := 0; i < len(a); i++ {
			for j := 0; j < len(a); j++ {
				for k := 0; k < len(a); k++ {
					c__[i][j] = a[i] * b[j]
				}
			}
		}
		chk.Matrix(tst, "a dy b", 1e-12, c, c_)
		chk.Matrix(tst, "a dy b", 1e-12, c, c__)

		// dot product
		d := M_Alloc2(nd[m])
		D := Alloc2()
		for i := 0; i < 3; i++ {
			for j := 0; j < 3; j++ {
				for k := 0; k < 3; k++ {
					D[i][j] += A[i][k] * B[k][j]
				}
			}
//.........这里部分代码省略.........

func Test_isofun01(tst *testing.T) {

	//verbose()
	chk.PrintTitle("isofun01. rounded cone")

	// SMP director parameters
	//  Note:
	//   1) eps and ϵ have an effect on how close to DP/MC SMP will be
	//   2) as eps increases, SMP is closer to DP/MC
	//   3) as ϵ increases, SMP is closer to DP/MC
	//   4) eps also changes the shape of FC surface
	a, b, eps, ϵ := -1.0, 0.5, 1e-3, 1e-3
	shift := 1.0

	// radius
	r := 2.0

	// failure crit parameters and number of stress components
	φ, ncp := 30.0, 6

	// q/p coefficient
	μ := SmpCalcμ(φ, a, b, eps, ϵ)
	io.Pforan("μ = %v\n", μ)

	// isotropic functions
	var o IsoFun
	o.Init(a, b, eps, ϵ, shift, ncp, rounded_cone_ffcn, rounded_cone_gfcn, rounded_cone_hfcn)

	// plot
	if false {
		//if true {
		σcCte := 10.0
		M := Phi2M(φ, "oct")
		rmin, rmax := 0.0, 1.28*M*σcCte
		nr, nα := 31, 81
		//nr,   nα   := 31, 1001
		npolarc := true
		simplec := true
		only0 := true
		grads := false
		showpts := false
		ferr := 10.0
		PlotOct("fig_isofun01.png", σcCte, rmin, rmax, nr, nα, φ, o.Fa, o.Ga,
			npolarc, simplec, only0, grads, showpts, true, true, ferr, r, μ)
	}

	// 3D view
	if false {
		//if true {
		grads := true
		gftol := 5e-2
		o.View(10, nil, grads, gftol, func(e *vtk.IsoSurf) {
			e.Nlevels = 7
		}, r, μ)
	}

	// constants
	ver := chk.Verbose
	tol := 1e-6
	tol2 := 1e-6
	tolq := tol2

	// check gradients
	for idxA := 0; idxA < len(test_nd); idxA++ {
		//for idxA := 0; idxA < 1; idxA++ {
		//for idxA := 2; idxA < 3; idxA++ {
		//for idxA := 10; idxA < 11; idxA++ {
		//for idxA := 11; idxA < 12; idxA++ {
		//for idxA := 12; idxA < 13; idxA++ {

		// tensor
		AA := test_AA[idxA]
		A := M_Alloc2(3)
		Ten2Man(A, AA)
		io.PfYel("\n\ntst # %d ###################################################################################\n", idxA)
		io.Pfblue2("A = %v\n", A)

		// function evaluation and shifted eigenvalues
		fval, err := o.Fa(A, r, μ)
		if err != nil {
			chk.Panic("cannot compute F(A):\n%v", err)
		}
		io.Pfpink("shift = %v\n", shift)
		io.Pforan("p, q  = %v, %v\n", o.p, o.q)
		io.Pforan("f(A)  = %v\n", fval)

		// change tolerances
		tol3 := tol2
		tol2_tmp := tol2
		switch idxA {
		case 7:
			tolq = 1e-5
		case 10:
			tolq = 2508  // TODO
			tol3 = 0.772 // TODO: check why test # 10 fails with d2f/dAdA
		case 11:
			tol2 = 0.0442 // TODO: check this
			tol3 = 440    //TODO: check this
		case 12:
			tol2 = 1e-3
//.........这里部分代码省略.........

func main() {

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

	if mpi.Rank() == 0 {
		io.PfYel("\nTest MPI 01\n")
	}
	if mpi.Size() != 3 {
		chk.Panic("this test needs 3 processors")
	}
	n := 11
	x := make([]float64, n)
	id, sz := mpi.Rank(), mpi.Size()
	start, endp1 := (id*n)/sz, ((id+1)*n)/sz
	for i := start; i < endp1; i++ {
		x[i] = float64(i)
	}

	// Barrier
	mpi.Barrier()

	io.Pfgrey("x @ proc # %d = %v\n", id, x)

	// SumToRoot
	r := make([]float64, n)
	mpi.SumToRoot(r, x)
	var tst testing.T
	if id == 0 {
		chk.Vector(&tst, fmt.Sprintf("SumToRoot:       r @ proc # %d", id), 1e-17, r, []float64{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10})
	} else {
		chk.Vector(&tst, fmt.Sprintf("SumToRoot:       r @ proc # %d", id), 1e-17, r, make([]float64, n))
	}

	// BcastFromRoot
	r[0] = 666
	mpi.BcastFromRoot(r)
	chk.Vector(&tst, fmt.Sprintf("BcastFromRoot:   r @ proc # %d", id), 1e-17, r, []float64{666, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10})

	// AllReduceSum
	setslice(x)
	w := make([]float64, n)
	mpi.AllReduceSum(x, w)
	chk.Vector(&tst, fmt.Sprintf("AllReduceSum:    w @ proc # %d", id), 1e-17, w, []float64{110, 110, 110, 1021, 1021, 1021, 2032, 2032, 2032, 3043, 3043})

	// AllReduceSumAdd
	setslice(x)
	y := []float64{-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000}
	mpi.AllReduceSumAdd(y, x, w)
	chk.Vector(&tst, fmt.Sprintf("AllReduceSumAdd: y @ proc # %d", id), 1e-17, y, []float64{-890, -890, -890, 21, 21, 21, 1032, 1032, 1032, 2043, 2043})

	// AllReduceMin
	setslice(x)
	mpi.AllReduceMin(x, w)
	chk.Vector(&tst, fmt.Sprintf("AllReduceMin:    x @ proc # %d", id), 1e-17, x, []float64{0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3})

	// AllReduceMax
	setslice(x)
	mpi.AllReduceMax(x, w)
	chk.Vector(&tst, fmt.Sprintf("AllReduceMax:    x @ proc # %d", id), 1e-17, x, []float64{100, 100, 100, 1000, 1000, 1000, 2000, 2000, 2000, 3000, 3000})
}