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 Test_spo751re(tst *testing.T) {

	//verbose()
	chk.PrintTitle("spo751re. Richardson extrapolation")

	// run simulation
	if !Start("data/spo751re.sim", true, chk.Verbose) {
		io.Pfred("start failed\n")
		return
	}

	// make sure to flush log
	defer End()

	// run simulation
	if !Run() {
		io.Pfred("run failed\n")
	}

	// plot
	//if true {
	if false {
		plot_spo751("spo751re")
	}
}

// StatF computes statistical information corresponding to objective function idxF
func StatF(o *Optimiser, idxF int, verbose bool) (fmin, fave, fmax, fdev float64, F []float64) {
	nsamples := len(o.BestOvas[idxF])
	if nsamples == 0 {
		if verbose {
			io.Pfred("there are no samples for statistical analysis\n")
		}
		return
	}
	F = make([]float64, nsamples)
	if nsamples == 1 {
		F[0] = o.BestOvas[idxF][0]
		fmin, fave, fmax = F[0], F[0], F[0]
		return
	}
	for i, f := range o.BestOvas[idxF] {
		F[i] = f
	}
	fmin, fave, fmax, fdev = rnd.StatBasic(F, true)
	if verbose {
		str := "\n"
		if len(o.RptFref) == o.Nova {
			str = io.Sf(" (%g)\n", o.RptFref[idxF])
		}
		io.Pf("fmin = %g\n", fmin)
		io.Pf("fave = %g"+str, fave)
		io.Pf("fmax = %g\n", fmax)
		io.Pf("fdev = %g\n", fdev)
		o.fix_formatting_data()
		io.Pf(rnd.BuildTextHist(nice(fmin, o.HistNdig)-o.HistDelFmin, nice(fmax, o.HistNdig)+o.HistDelFmax,
			o.HistNsta, F, o.HistFmt, o.HistLen))
	}
	return
}

// StatMulti prints statistical analysis for multi-objective problems
//  emin, eave, emax, edev -- errors on f1(f0)
//  key -- "IGD" if IGD values are available. In this case e{...} are IGD values
func StatMulti(o *Optimiser, verbose bool) (key string, emin, eave, emax, edev float64, E []float64) {
	if len(o.Multi_err) < 2 && len(o.Multi_IGD) < 2 {
		io.Pfred("there are no samples for statistical analysis\n")
		return
	}
	o.fix_formatting_data()
	n := len(o.Multi_err)
	key = "E"
	if n < 2 {
		n = len(o.Multi_IGD)
		key = "IGD"
	}
	E = make([]float64, n)
	if key == "E" {
		copy(E, o.Multi_err)
	} else {
		copy(E, o.Multi_IGD)
	}
	emin, eave, emax, edev = rnd.StatBasic(E, true)
	if verbose {
		io.Pf("\nerror on Pareto front (multi)\n")
		io.Pf("%smin = %g\n", key, emin)
		io.Pf("%save = %g\n", key, eave)
		io.Pf("%smax = %g\n", key, emax)
		io.Pf("%sdev = %g\n", key, edev)
		io.Pf(rnd.BuildTextHist(nice(emin, o.HistNdig)-o.HistDelEmin, nice(emax, o.HistNdig)+o.HistDelEmax,
			o.HistNsta, E, o.HistFmt, o.HistLen))
	}
	return
}

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

// Init initialises model
func (o *DruckerPrager) Init(ndim int, pstress bool, prms fun.Prms) (err error) {

	// parse parameters
	err = o.SmallElasticity.Init(ndim, pstress, prms)
	if err != nil {
		return
	}
	for _, p := range prms {
		switch p.N {
		case "M":
			o.M = p.V
		case "Mb":
			o.Mb = p.V
		case "qy0":
			o.qy0 = p.V
		case "H":
			o.H = p.V
		case "E", "nu", "l", "G", "K", "rho":
		case "c", "phi", "typ":
			io.Pfred("dp: warning: handling of 'c', 'phi' and 'typ' parameters is not implemented yet\n")
		default:
			return chk.Err("dp: parameter named %q is incorrect\n", p.N)
		}
	}

	// auxiliary structures
	o.ten = make([]float64, o.Nsig)
	return
}

func Test_cxint01(tst *testing.T) {

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

	var ops OpsData
	ops.SetDefault()
	ops.Pc = 1
	ops.Ncuts = 1

	A := []int{1, 2}
	B := []int{-1, -2}
	a := make([]int, len(A))
	b := make([]int, len(A))
	IntCrossover(a, b, A, B, 0, &ops)
	io.Pfred("A = %2d\n", A)
	io.PfRed("B = %2d\n", B)
	io.Pfcyan("a = %2d\n", a)
	io.Pfblue2("b = %2d\n", b)
	chk.Ints(tst, "a", a, []int{1, -2})
	chk.Ints(tst, "b", b, []int{-1, 2})
	io.Pf("\n")

	A = []int{1, 2, 3, 4, 5, 6, 7, 8}
	B = []int{-1, -2, -3, -4, -5, -6, -7, -8}
	a = make([]int, len(A))
	b = make([]int, len(A))
	ops.Cuts = []int{1, 3}
	IntCrossover(a, b, A, B, 0, &ops)
	io.Pfred("A = %2v\n", A)
	io.PfRed("B = %2v\n", B)
	io.Pfcyan("a = %2v\n", a)
	io.Pfblue2("b = %2v\n", b)
	chk.Ints(tst, "a", a, []int{1, -2, -3, 4, 5, 6, 7, 8})
	chk.Ints(tst, "b", b, []int{-1, 2, 3, -4, -5, -6, -7, -8})

	ops.Cuts = []int{5, 7}
	IntCrossover(a, b, A, B, 0, &ops)
	io.Pfred("A = %2v\n", A)
	io.PfRed("B = %2v\n", B)
	io.Pfcyan("a = %2v\n", a)
	io.Pfblue2("b = %2v\n", b)
	chk.Ints(tst, "a", a, []int{1, 2, 3, 4, 5, -6, -7, 8})
	chk.Ints(tst, "b", b, []int{-1, -2, -3, -4, -5, 6, 7, -8})
}

// rjoint_DebugKb defines a global function to debug Kb for rjoint-elements
func rjoint_DebugKb(fem *FEM, o *testKb) {
	fem.DebugKb = func(d *Domain, it int) {

		elem := d.Elems[o.eid]
		if e, ok := elem.(*Rjoint); ok {

			// skip?
			o.it = it
			o.t = d.Sol.T
			if o.skip() {
				return
			}

			// copy states and solution
			nip := len(e.Rod.IpsElem)
			states := make([]*msolid.OnedState, nip)
			statesBkp := make([]*msolid.OnedState, nip)
			for i := 0; i < nip; i++ {
				states[i] = e.States[i].GetCopy()
				statesBkp[i] = e.StatesBkp[i].GetCopy()
			}
			o.aux_arrays(d)

			// make sure to restore states and solution
			defer func() {
				for i := 0; i < nip; i++ {
					e.States[i].Set(states[i])
					e.StatesBkp[i].Set(statesBkp[i])
				}
				copy(d.Sol.ΔY, o.ΔYbkp)
			}()

			// define restore function
			restore := func() {
				if it == 0 {
					for k := 0; k < nip; k++ {
						e.States[k].Set(states[k])
					}
					return
				}
				for k := 0; k < nip; k++ {
					e.States[k].Set(statesBkp[k])
				}
			}

			// check
			o.check("Krr", d, e, e.Rod.Umap, e.Rod.Umap, e.Krr, restore)
			o.check("Krs", d, e, e.Rod.Umap, e.Sld.Umap, e.Krs, restore)
			o.check("Ksr", d, e, e.Sld.Umap, e.Rod.Umap, e.Ksr, restore)
			o.check("Kss", d, e, e.Sld.Umap, e.Sld.Umap, e.Kss, restore)
		} else {
			io.Pfred("warning: eid=%d does not correspond to Rjoint element\n", o.eid)
		}
	}
	return
}

// DistPointLine computes the distance from p to line passing through a -> b
func DistPointLine(p, a, b *Point, tol float64, verbose bool) float64 {
	ns := NewSegment(a, b)
	vs := NewSegment(p, a)
	nn := ns.Len()
	if nn < tol { // point-point distance
		if verbose {
			io.Pfred("basicgeom.go: DistPointLine: __WARNING__ point-point distance too small:\n p=%v a=%v b=%v\n", p, a, b)
		}
		return vs.Len()
	}
	n := ns.Vector(1.0 / nn)
	v := vs.Vector(1.0)
	s := VecDot(v, n)
	l := VecNewAdd(1, v, -s, n) // l := v - dot(v,n) * n
	return VecNorm(l)
}

func Test_bspline02(tst *testing.T) {

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

	//               0 1 2 3 4 5 6 7 8 9 10
	T := []float64{0, 0, 0, 1, 2, 3, 4, 4, 5, 5, 5}
	sol := []int{2, 2, 3, 3, 4, 4, 5, 5, 7, 7, 7}
	var s Bspline
	s.Init(T, 2)
	s.SetControl([][]float64{{0, 0}, {0.5, 1}, {1, 0}, {1.5, 0}, {2, 1}, {2.5, 1}, {3, 0.5}, {3.5, 0}})

	tt := utl.LinSpace(0, 5, 11)
	for k, t := range tt {
		span := s.find_span(t)
		io.Pforan("t=%.4f  =>  span=%v\n", t, span)
		if span != sol[k] {
			chk.Panic("find_span failed: t=%v  span => %d != %d", t, span, sol[k])
		}
	}

	tol := 1e-14
	np := len(tt)
	xx, yy := make([]float64, np), make([]float64, np)
	for k, t := range tt {
		t0 := time.Now()
		pa := s.Point(t, 0) // 0 => CalcBasis
		io.Pf("Point(rec): dtime = %v\n", time.Now().Sub(t0))
		t0 = time.Now()
		pb := s.Point(t, 1) // 1 => RecursiveBasis
		io.Pf("Point:      dtime = %v\n", time.Now().Sub(t0))
		xx[k], yy[k] = pb[0], pb[1]
		io.Pfred("pa - pb = %v, %v\n", pa[0]-pb[0], pa[1]-pb[1])
		chk.Vector(tst, "Point", tol, pa, pb)
	}

	if chk.Verbose {
		npts := 201
		plt.SetForPng(0.75, 300, 150)
		str0 := ",lw=2"
		str1 := ",ls='none',marker='+',color='cyan',markevery=10"
		s.Draw2d(str0, "", npts, 0) // 0 => CalcBasis
		s.Draw2d(str1, "", npts, 1) // 1 => RecursiveBasis
		plt.Plot(xx, yy, "'bo', clip_on=0")
		plt.SaveD("/tmp/gosl/gm", "bspline02.png")
	}
}

func Test_bh14c(tst *testing.T) {

	//verbose()
	chk.PrintTitle("bh14c. using go-routines")

	// channels
	nch := 2
	done := make(chan int, nch)

	// allocate structures and set stage
	analyses := make([]*FEM, nch)
	for i := 0; i < nch; i++ {
		analyses[i] = NewFEM("data/bh14.sim", io.Sf("ch%d", i), true, true, false, false, false, i)
		err := analyses[i].SetStage(0)
		if err != nil {
			tst.Errorf("SetStage failed:\n%v", err)
			return
		}
	}

	// run all analyses
	for i := 0; i < nch; i++ {
		go func(analysis *FEM) {
			err := analysis.SolveOneStage(0, true)
			if err != nil {
				io.Pfred("SolveOneStage failed:\n%v", err)
			}
			done <- 1
		}(analyses[i])
	}

	// wait
	for i := 0; i < nch; i++ {
		<-done
	}

	// check
	skipK := false
	tolK := 1e-10
	tolu := 1e-15
	tols := 1e-13
	for i := 0; i < nch; i++ {
		TestingCompareResultsU(tst, "data/bh14.sim", "cmp/bh14.cmp", io.Sf("ch%d", i), tolK, tolu, tols, skipK, chk.Verbose)
	}
}

// StatF1F0 prints statistical analysis for two-objective problems
//  emin, eave, emax, edev -- errors on f1(f0)
//  lmin, lave, lmax, ldev -- arc-lengths along f1(f0) curve
func StatF1F0(o *Optimiser, verbose bool) (emin, eave, emax, edev float64, E []float64, lmin, lave, lmax, ldev float64, L []float64) {
	if len(o.F1F0_err) == 0 && len(o.F1F0_arcLen) == 0 {
		io.Pfred("there are no samples for statistical analysis\n")
		return
	}
	o.fix_formatting_data()
	if len(o.F1F0_err) > 2 {
		E = make([]float64, len(o.F1F0_err))
		copy(E, o.F1F0_err)
		emin, eave, emax, edev = rnd.StatBasic(E, true)
		if verbose {
			io.Pf("\nerror on Pareto front\n")
			io.Pf("emin = %g\n", emin)
			io.Pf("eave = %g\n", eave)
			io.Pf("emax = %g\n", emax)
			io.Pf("edev = %g\n", edev)
			io.Pf(rnd.BuildTextHist(nice(emin, o.HistNdig)-o.HistDelEmin, nice(emax, o.HistNdig)+o.HistDelEmax,
				o.HistNsta, E, o.HistFmt, o.HistLen))
		}
	}
	if len(o.F1F0_arcLen) > 2 {
		den := 1.0
		if o.F1F0_arcLenRef > 0 {
			den = o.F1F0_arcLenRef
		}
		L := make([]float64, len(o.F1F0_arcLen))
		for i, l := range o.F1F0_arcLen {
			L[i] = l / den
		}
		lmin, lave, lmax, ldev = rnd.StatBasic(L, true)
		if verbose {
			io.Pf("\nnormalised arc length along Pareto front (ref = %g)\n", o.F1F0_arcLenRef)
			io.Pf("lmin = %g\n", lmin)
			io.Pf("lave = %g\n", lave)
			io.Pf("lmax = %g\n", lmax)
			io.Pf("ldev = %g\n", ldev)
			io.Pf(rnd.BuildTextHist(nice(lmin, o.HistNdig)-o.HistDelEmin, nice(lmax, o.HistNdig)+o.HistDelEmax,
				o.HistNsta, L, o.HistFmt, o.HistLen))
		}
	}
	return
}

func (o *RichardsonExtrap) divergence_control(d *Domain, name string) (docontinue bool) {
	if o.diverging {
		if Global.Verbose {
			io.Pfred(". . . %s: iterations diverging (%2d) . . .\n", name, o.ndiverg+1)
		}
		d.restore()
		o.Δtcpy = o.Δt
		o.Δt *= 0.5
		o.ndiverg += 1
		o.nreject += 1
		o.laststep = false
		o.prevdiv = true
		return true
	}
	if o.prevdiv && false {
		o.Δt = o.Δtcpy
		o.ndiverg = 0
		o.prevdiv = false
	}
	return false
}

// Update updates state
func (o *PrincStrainsUp) Update(s *State, ε, Δε []float64, eid, ipid int, time float64) (err error) {

	// debugging
	if o.DbgOn {
		o.dbg_init(s, ε, Δε, eid, ipid)
	}

	// trial strains
	for i := 0; i < o.Nsig; i++ {
		s.EpsE[i] += Δε[i]
	}
	copy(s.EpsTr, s.EpsE)

	// eigenvalues/projectors of trial elastic strain
	err = tsr.M_EigenValsProjsNum(o.P, o.Lεetr, s.EpsTr)
	if err != nil {
		return
	}

	// trial stresses
	o.Mdl.E_CalcSig(o.Lσ, o.Lεetr)

	// debugging
	if o.DbgOn {
		defer o.dbg_end(s, ε, eid, ipid)() // TODO: fix this
	}

	// check loading condition => elastic update?
	ftr := o.Mdl.L_YieldFunc(o.Lσ, s.Alp)
	if ftr <= o.Fzero {
		s.Dgam = 0
		s.Loading = false
		for i := 0; i < o.Nsig; i++ {
			s.Sig[i] = o.Lσ[0]*o.P[0][i] + o.Lσ[1]*o.P[1][i] + o.Lσ[2]*o.P[2][i]
		}
		return
	}

	// initial values
	for i := 0; i < 3; i++ {
		o.x[i] = o.Lεetr[i]
	}
	for i := 0; i < o.Nalp; i++ {
		o.αn[i] = s.Alp[i]
		o.x[3+i] = s.Alp[i]
	}
	o.x[3+o.Nalp] = 0 // Δγ

	// check Jacobian
	if o.ChkJac {
		var cnd float64
		cnd, err = o.nls.CheckJ(o.x, o.ChkJacTol, true, o.ChkSilent)
		io.Pfred("before: cnd(J) = %v\n", cnd)
	}

	// modify b
	bsmp := o.Mdl.Get_bsmp()
	if bsmp > 0 && o.Nbsmp > 1 {
		o.Mdl.Set_bsmp(0)
		defer func() { o.Mdl.Set_bsmp(bsmp) }()
		δb := bsmp / float64(o.Nbsmp-1)
		for i := 0; i < o.Nbsmp; i++ {
			b := float64(i) * δb
			err = o.do_solve(b, eid, ipid, time)
			if err != nil {
				return
			}
		}
	} else {
		err = o.do_solve(bsmp, eid, ipid, time)
		if err != nil {
			return
		}
	}

	// check Jacobian again
	if o.ChkJac {
		var cnd float64
		cnd, err = o.nls.CheckJ(o.x, o.ChkJacTol, true, o.ChkSilent)
		io.Pfred("after: cnd(J) = %v\n", cnd)
		if err != nil {
			return
		}
	}

	// set new state
	εe, α, Δγ := o.x[:3], o.x[3:3+o.Nalp], o.x[3+o.Nalp]
	o.Mdl.E_CalcSig(o.Lσ, εe)
	for i := 0; i < o.Nsig; i++ {
		s.Sig[i] = o.Lσ[0]*o.P[0][i] + o.Lσ[1]*o.P[1][i] + o.Lσ[2]*o.P[2][i]
		s.EpsE[i] = εe[0]*o.P[0][i] + εe[1]*o.P[1][i] + εe[2]*o.P[2][i]
	}
	copy(s.Alp, α)
	s.Dgam = Δγ
	s.Loading = true
	return
}

// Run runs optimisations
func (o *SimpleFltProb) Run(verbose bool) {

	// benchmark
	if verbose {
		time0 := time.Now()
		defer func() {
			io.Pfblue2("\ncpu time = %v\n", time.Now().Sub(time0))
		}()
	}

	// run all trials
	for itrial := 0; itrial < o.C.Ntrials; itrial++ {

		// reset populations
		if itrial > 0 {
			for id, isl := range o.Evo.Islands {
				isl.Pop = o.C.PopFltGen(id, o.C)
				isl.CalcOvs(isl.Pop, 0)
				isl.CalcDemeritsAndSort(isl.Pop)
			}
		}

		// run evolution
		o.Evo.Run()

		// results
		xbest := o.Evo.Best.GetFloats()
		o.Fcn(o.ff[0], o.gg[0], o.hh[0], xbest)

		// check if best is unfeasible
		unfeasible := false
		for _, g := range o.gg[0] {
			if g < 0 {
				unfeasible = true
				break
			}
		}
		for _, h := range o.hh[0] {
			if math.Abs(h) > o.C.Eps1 {
				unfeasible = true
				break
			}
		}

		// feasible results
		if !unfeasible {
			for i, x := range xbest {
				o.Xbest[o.Nfeasible][i] = x
			}
			o.Nfeasible++
		}

		// message
		if verbose {
			io.Pfyel("%3d x*="+o.NumfmtX+" f="+o.NumfmtF, itrial, xbest, o.ff[0])
			if unfeasible {
				io.Pfred(" unfeasible\n")
			} else {
				io.Pfgreen(" ok\n")
			}
		}

		// best populations
		if o.C.DoPlot {
			if o.Nfeasible == 1 {
				o.PopsBest = o.Evo.GetPopulations()
			} else {
				fcur := utl.DblCopy(o.ff[0])
				o.Fcn(o.ff[0], o.gg[0], o.hh[0], o.Xbest[o.Nfeasible-1])
				cur_dom, _ := utl.DblsParetoMin(fcur, o.ff[0])
				if cur_dom {
					o.PopsBest = o.Evo.GetPopulations()
				}
			}
		}
	}
}

func Test_cx01(tst *testing.T) {

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

	var ops OpsData
	ops.SetDefault()
	ops.Pc = 1
	ops.Cuts = []int{2, 4}

	A := []float64{1.1, 2.2, 3.3, 4.4, 5.5, 6.6}
	B := []float64{9.1, 9.2, 9.3, 9.4, 9.5, 9.6}
	a := make([]float64, len(A))
	b := make([]float64, len(A))
	FltCrossover(a, b, A, B, 0, &ops)
	io.Pfred("A = %3v\n", A)
	io.PfRed("B = %3v\n", B)
	io.Pfcyan("a = %3v\n", a)
	io.Pfblue2("b = %3v\n", b)
	chk.Vector(tst, "a", 1e-17, a, []float64{1.1, 2.2, 9.3, 9.4, 5.5, 6.6})
	chk.Vector(tst, "b", 1e-17, b, []float64{9.1, 9.2, 3.3, 4.4, 9.5, 9.6})
	io.Pf("\n")

	C := []string{"A", "B", "C", "D", "E", "F"}
	D := []string{"-", "o", "+", "@", "*", "&"}
	c := make([]string, len(A))
	d := make([]string, len(A))
	ops.Cuts = []int{1, 3}
	StrCrossover(c, d, C, D, 0, &ops)
	io.Pfred("C = %3v\n", C)
	io.PfRed("D = %3v\n", D)
	io.Pfcyan("c = %3v\n", c)
	io.Pfblue2("d = %3v\n", d)
	chk.Strings(tst, "c", c, []string{"A", "o", "+", "D", "E", "F"})
	chk.Strings(tst, "d", d, []string{"-", "B", "C", "@", "*", "&"})
	io.Pf("\n")

	E := [][]byte{[]byte("A"), []byte("B"), []byte("C"), []byte("D"), []byte("E"), []byte("F")}
	F := [][]byte{[]byte("-"), []byte("o"), []byte("+"), []byte("@"), []byte("*"), []byte("&")}
	e := make([][]byte, len(A))
	f := make([][]byte, len(A))
	for i := 0; i < len(A); i++ {
		e[i] = make([]byte, 1)
		f[i] = make([]byte, 1)
	}
	ops.Cuts = []int{1, 3}
	BytCrossover(e, f, E, F, 0, &ops)
	io.Pfred("E = %3s\n", E)
	io.PfRed("F = %3s\n", F)
	io.Pfcyan("e = %3s\n", e)
	io.Pfblue2("f = %3s\n", f)
	e_s := make([]string, len(A))
	f_s := make([]string, len(A))
	for i := 0; i < len(A); i++ {
		e_s[i] = string(e[i])
		f_s[i] = string(f[i])
	}
	chk.Strings(tst, "e_s", e_s, []string{"A", "o", "+", "D", "E", "F"})
	chk.Strings(tst, "f_s", f_s, []string{"-", "B", "C", "@", "*", "&"})
	io.Pf("\n")

	G := []byte("ABCDEF")
	H := []byte("[email protected]*&")
	g := make([]byte, len(A))
	h := make([]byte, len(A))
	ops.Cuts = []int{1, 3}
	KeyCrossover(g, h, G, H, 0, &ops)
	io.Pfred("G = %3v\n", G)
	io.PfRed("H = %3v\n", H)
	io.Pfcyan("g = %3v\n", g)
	io.Pfblue2("h = %3v\n", h)
	g_s := make([]string, len(A))
	h_s := make([]string, len(A))
	for i := 0; i < len(A); i++ {
		g_s[i] = string(g[i])
		h_s[i] = string(h[i])
	}
	chk.Strings(tst, "g_s", g_s, []string{"A", "o", "+", "D", "E", "F"})
	chk.Strings(tst, "h_s", h_s, []string{"-", "B", "C", "@", "*", "&"})
	io.Pf("\n")

	M := []Func_t{func(i *Individual) string { return "A" }, func(i *Individual) string { return "B" }, func(i *Individual) string { return "C" }, func(i *Individual) string { return "D" }, func(i *Individual) string { return "E" }, func(i *Individual) string { return "F" }}
	N := []Func_t{func(i *Individual) string { return "-" }, func(i *Individual) string { return "o" }, func(i *Individual) string { return "+" }, func(i *Individual) string { return "@" }, func(i *Individual) string { return "*" }, func(i *Individual) string { return "&" }}
	m := make([]Func_t, len(A))
	n := make([]Func_t, len(A))
	ops.Cuts = []int{1, 3}
	FunCrossover(m, n, M, N, 0, &ops)
	io.Pfred("M = %3v\n", M)
	io.PfRed("N = %3v\n", N)
	io.Pfcyan("m = %3v\n", m)
	io.Pfblue2("n = %3v\n", n)
	m_s := make([]string, len(A))
	n_s := make([]string, len(A))
	for i := 0; i < len(A); i++ {
		m_s[i] = m[i](nil)
		n_s[i] = n[i](nil)
	}
	chk.Strings(tst, "m_s", m_s, []string{"A", "o", "+", "D", "E", "F"})
	chk.Strings(tst, "n_s", n_s, []string{"-", "B", "C", "@", "*", "&"})
	io.Pf("\n")
//.........这里部分代码省略.........

func main() {

	// given the following matrix of complex numbers:
	//      _                                                  _
	//     |  19.73    12.11-i      5i        0          0      |
	//     |  -0.51i   32.3+7i    23.07       i          0      |
	// A = |    0      -0.51i    70+7.3i     3.95    19+31.83i  |
	//     |    0        0        1+1.1i    50.17      45.51    |
	//     |_   0        0          0      -9.351i       55    _|
	//
	// and the following vector:
	//      _                  _
	//     |    77.38+8.82i     |
	//     |   157.48+19.8i     |
	// b = |  1175.62+20.69i    |
	//     |   912.12-801.75i   |
	//     |_     550-1060.4i  _|
	//
	// solve:
	//         A.x = b
	//
	// the solution is:
	//      _            _
	//     |     3.3-i    |
	//     |    1+0.17i   |
	// x = |      5.5     |
	//     |       9      |
	//     |_  10-17.75i _|

	// flag indicating to store (real,complex) values in monolithic form => 1D array
	xzmono := false

	// input matrix in Complex Triplet format
	var A la.TripletC
	A.Init(5, 5, 16, xzmono) // 5 x 5 matrix with 16 non-zeros

	// first column
	A.Put(0, 0, 19.73, 0) // i=0, j=0, real=19.73, complex=0
	A.Put(1, 0, 0, -0.51) // i=1, j=0, real=0, complex=-0.51

	// second column
	A.Put(0, 1, 12.11, -1) // i=0, j=1, real=12.11, complex=-1
	A.Put(1, 1, 32.3, 7)
	A.Put(2, 1, 0, -0.51)

	// third column
	A.Put(0, 2, 0, 5)
	A.Put(1, 2, 23.07, 0)
	A.Put(2, 2, 70, 7.3)
	A.Put(3, 2, 1, 1.1)

	// fourth column
	A.Put(1, 3, 0, 1)
	A.Put(2, 3, 3.95, 0)
	A.Put(3, 3, 50.17, 0)
	A.Put(4, 3, 0, -9.351)

	// fifth column
	A.Put(2, 4, 19, 31.83)
	A.Put(3, 4, 45.51, 0)
	A.Put(4, 4, 55, 0)

	// right-hand-side
	b := []complex128{
		77.38 + 8.82i,
		157.48 + 19.8i,
		1175.62 + 20.69i,
		912.12 - 801.75i,
		550 - 1060.4i,
	}

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

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

	// initialise solver (C)omplex
	err := lis.InitC(&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
	}

	// auxiliary variables
	bR, bC := la.ComplexToRC(b)   // real and complex components of b
	xR := make([]float64, len(b)) // real compoments of x
	xC := make([]float64, len(b)) // complex compoments of x

	// solve (C)omplex
//.........这里部分代码省略.........

func main() {

	// given the following matrix of complex numbers:
	//      _                                                  _
	//     |  19.73    12.11-i      5i        0          0      |
	//     |  -0.51i   32.3+7i    23.07       i          0      |
	// A = |    0      -0.51i    70+7.3i     3.95    19+31.83i  |
	//     |    0        0        1+1.1i    50.17      45.51    |
	//     |_   0        0          0      -9.351i       55    _|
	//
	// and the following vector:
	//      _                  _
	//     |    77.38+8.82i     |
	//     |   157.48+19.8i     |
	// b = |  1175.62+20.69i    |
	//     |   912.12-801.75i   |
	//     |_     550-1060.4i  _|
	//
	// solve:
	//         A.x = b
	//
	// the solution is:
	//      _            _
	//     |     3.3-i    |
	//     |    1+0.17i   |
	// x = |      5.5     |
	//     |       9      |
	//     |_  10-17.75i _|

	// flag indicating to store (real,complex) values in monolithic form => 1D array
	xzmono := false

	// input matrix in Complex Triplet format
	var A la.TripletC
	A.Init(5, 5, 16, xzmono) // 5 x 5 matrix with 16 non-zeros

	// first column
	A.Put(0, 0, 19.73, 0) // i=0, j=0, real=19.73, complex=0
	A.Put(1, 0, 0, -0.51) // i=1, j=0, real=0, complex=-0.51

	// second column
	A.Put(0, 1, 12.11, -1) // i=0, j=1, real=12.11, complex=-1
	A.Put(1, 1, 32.3, 7)
	A.Put(2, 1, 0, -0.51)

	// third column
	A.Put(0, 2, 0, 5)
	A.Put(1, 2, 23.07, 0)
	A.Put(2, 2, 70, 7.3)
	A.Put(3, 2, 1, 1.1)

	// fourth column
	A.Put(1, 3, 0, 1)
	A.Put(2, 3, 3.95, 0)
	A.Put(3, 3, 50.17, 0)
	A.Put(4, 3, 0, -9.351)

	// fifth column
	A.Put(2, 4, 19, 31.83)
	A.Put(3, 4, 45.51, 0)
	A.Put(4, 4, 55, 0)

	// right-hand-side
	b := []complex128{
		77.38 + 8.82i,
		157.48 + 19.8i,
		1175.62 + 20.69i,
		912.12 - 801.75i,
		550 - 1060.4i,
	}

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

	// output
	a := A.ToMatrix(nil)
	io.Pforan("A.x = b\n")
	la.PrintMatC("A", a.ToDense(), "(%5g", "%+6gi) ", false)
	la.PrintVecC("b", b, "(%g", "%+gi) ", false)
	la.PrintVecC("x", x, "(%.3f", "%+.3fi) ", false)
}

func Test_ccm01(tst *testing.T) {

	defer func() {
		if err := recover(); err != nil {
			io.Pfred("error = %v\n", err)
		}
	}()

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

	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)

	// allocate driver
	pr := 1.0
	ndim, pstress := 2, false
	simfnk, modelname := "test", "ccm"
	var drv Driver
	err := drv.Init(simfnk, modelname, ndim, pstress, []*fun.Prm{
		&fun.Prm{N: "phi", V: 25},
		&fun.Prm{N: "Mfix", V: 1},
		&fun.Prm{N: "c", V: 1},
		&fun.Prm{N: "lam", V: 0.1},
		&fun.Prm{N: "ocr", V: 1},
		&fun.Prm{N: "kap", V: 0.05},
		&fun.Prm{N: "kapb", V: 0.01},
		&fun.Prm{N: "G0", V: G},
		&fun.Prm{N: "pr", V: pr},
		&fun.Prm{N: "le", V: 0},
		&fun.Prm{N: "K0", V: K},
	})
	drv.CheckD = true
	//drv.CheckD = false
	drv.TolD = 1e-4
	drv.VerD = io.Verbose // verbose
	if err != nil {
		tst.Errorf("test failed: %v\n", err)
		return
	}

	// model
	ccm := drv.model.(*CamClayMod)

	// path
	p0 := 0.0
	DP := []float64{10, 1}
	DQ := []float64{0, 3}
	nincs := 1
	niout := 1
	noise := 0.0
	var pth Path
	if false {
		err = pth.SetPQstrain(ndim, nincs, niout, K, G, p0, DP, DQ, noise)
		if err != nil {
			tst.Errorf("test failed: %v\n", err)
			return
		}
	} else {
		pth.Sx = []float64{-1}
		pth.Sy = []float64{-2}
		pth.Sz = []float64{-1}
		pth.Ex = []float64{0, 0, 0.001, -0.004}  //, -0.005}
		pth.Ey = []float64{0, -0.005, 0, -0.002} //, -0.005}
		pth.Ez = []float64{0, 0, 0.001, -0.004}  //, -0.005}
		//pth.Ex = []float64{0, -0.0033333333333333335, -0.0028333333333333335}
		//pth.Ey = []float64{0, -0.0033333333333333335, -0.0028333333333333335}
		//pth.Ez = []float64{0, -0.0033333333333333335, -0.005333333333333334}
		pth.UseS = []int{0, 0}
		pth.UseE = []int{0, 1}
		pth.Init(ndim)
	}

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

	// plot
	//if true {
	if false {
		var plr Plotter
		plr.Pr = pr
		prop := 2.0
		plr.SetFig(false, false, prop, 400, "/tmp", "test_ccm01")
		plr.SetModel(ccm)
		n := len(drv.Res)
		plr.Pmin = -ccm.HE.pt
		plr.Pmax = drv.Res[n-1].Alp[0]
		qmax := ccm.CS.Mcs * (plr.Pmax - plr.Pmin)
		plr.PqLims = []float64{plr.Pmin, plr.Pmax, -qmax, qmax}
		plr.UsePmin = true
		plr.UsePmax = true
		plr.PreCor = drv.PreCor
		//plr.Plot(PlotSet8, drv.Res, drv.Eps, true, true)
//.........这里部分代码省略.........

//.........这里部分代码省略.........
			// check for continued divergence
			if LogErrCond(ndiverg >= Global.Sim.Solver.NdvgMax, "continuous divergence after %d steps reached", ndiverg) {
				return
			}

			// time increment
			Δt = Dt.F(t, nil) * md
			if t+Δt >= tf {
				Δt = tf - t
				lasttimestep = true
			}
			if Δt < Global.Sim.Solver.DtMin {
				if md < 1 {
					LogErrCond(true, "Δt increment is too small: %g < %g", Δt, Global.Sim.Solver.DtMin)
					return false
				}
				return true
			}

			// dynamic coefficients
			if LogErr(Global.DynCoefs.CalcBoth(Δt), "cannot compute dynamic coefficients") {
				return
			}

			// time update
			t += Δt
			for _, d := range domains {
				d.Sol.T = t
			}
			Δtout = DtOut.F(t, nil)

			// message
			if Global.Verbose {
				if !Global.Sim.Data.ShowR && !Global.Debug {
					io.PfWhite("%30.15f\r", t)
				}
			}

			// for all domains
			docontinue := false
			for _, d := range domains {

				// backup solution if divergence control is on
				if Global.Sim.Solver.DvgCtrl {
					d.backup()
				}

				// run iterations
				diverging, ok := run_iterations(t, Δt, d, &sum)
				if !ok {
					return
				}

				// restore solution and reduce time step if divergence control is on
				if Global.Sim.Solver.DvgCtrl {
					if diverging {
						if Global.Verbose {
							io.Pfred(". . . iterations diverging (%2d) . . .\n", ndiverg+1)
						}
						d.restore()
						t -= Δt
						d.Sol.T = t
						md *= 0.5
						ndiverg += 1
						docontinue = true
						break
					}
					ndiverg = 0
					md = 1.0
				}
			}
			if docontinue {
				continue
			}

			// perform output
			if t >= tout || lasttimestep {
				sum.OutTimes = append(sum.OutTimes, t)
				for _, d := range domains {
					//if true {
					if false {
						debug_print_p_results(d)
					}
					if false {
						debug_print_up_results(d)
					}
					if !d.Out(tidx) {
						break
					}
				}
				if Stop() {
					return
				}
				tout += Δtout
				tidx += 1
			}
		}
	}
	return true
}