func Test_invs07(tst *testing.T) {

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

	smp_a, smp_b, smp_β, smp_ϵ := -1.0, 0.0, 1.0, 1e-3

	σ := []float64{-1, -1, 0, 0}

	pcam, qcam, _ := M_pqw(σ)
	poct, qoct := pcam*SQ3, qcam*SQ2by3

	N := make([]float64, 3)
	n := make([]float64, 3)
	m := SmpDirector(N, σ, smp_a, smp_b, smp_β, smp_ϵ)
	SmpUnitDirector(n, m, N)
	psmp1, qsmp1, err := GenInvs(σ, n, smp_a)
	if err != nil {
		chk.Panic("M_GenInvs failed:\n%v", err)
	}

	psmp2, qsmp2, err := M_pq_smp(σ, smp_a, smp_b, smp_β, smp_ϵ)
	if err != nil {
		chk.Panic("M_pq_smp failed:\n%v", err)
	}
	io.Pforan("pcam,  qcam  = %v, %v\n", pcam, qcam)
	io.Pforan("poct,  qoct  = %v, %v\n", poct, qoct)
	io.Pforan("psmp1, qsmp1 = %v, %v\n", psmp1, qsmp1)
	io.Pforan("psmp2, qsmp2 = %v, %v\n", psmp2, qsmp2)
	chk.Scalar(tst, "p", 1e-15, psmp1, psmp2)
	chk.Scalar(tst, "q", 1e-15, qsmp1, qsmp2)
}

// ReadMsh reads .msh file
func ReadMsh(fnk string) (nurbss []*Nurbs) {

	// read file
	fn := fnk + ".msh"
	buf, err := io.ReadFile(fn)
	if err != nil {
		chk.Panic("ReadMsh cannot read file = '%s', %v'", fn, err)
	}

	// decode
	var dat Data
	err = json.Unmarshal(buf, &dat)
	if err != nil {
		chk.Panic("ReadMsh cannot unmarshal file = '%s', %v'", fn, err)
	}

	// list of vertices
	verts := make([][]float64, len(dat.Verts))
	for i, _ := range dat.Verts {
		verts[i] = make([]float64, 4)
		for j := 0; j < 4; j++ {
			verts[i][j] = dat.Verts[i].C[j]
		}
	}

	// allocate NURBSs
	nurbss = make([]*Nurbs, len(dat.Nurbss))
	for i, v := range dat.Nurbss {
		nurbss[i] = new(Nurbs)
		nurbss[i].Init(v.Gnd, v.Ords, v.Knots)
		nurbss[i].SetControl(verts, v.Ctrls)
	}
	return
}

// HashPoint returns a unique id of a point
func HashPoint(x, xmin, xdel []float64, tol float64) int {
	if tol < 1e-15 {
		chk.Panic("HashPoint: minimum tolerance must be 1e-15. %v is invalid", tol)
	}
	coefs := []float64{11, 101, 1001}
	n := utl.Imin(len(x), 3)
	var hash, xbar float64
	for i := 0; i < n; i++ {
		if x[i] < xmin[i] {
			chk.Panic("HashPoint: coordinate is outside range: %v < %v", x[i], xmin[i])
		}
		if x[i] > xmin[i]+xdel[i] {
			chk.Panic("HashPoint: coordinate is outside range: %v > %v", x[i], xmin[i]+xdel[i])
		}
		if xdel[i] > 0 {
			xbar = (x[i] - xmin[i]) / xdel[i]
			if xbar < 0 {
				xbar = 0
			}
			if xbar > 1 {
				xbar = 1
			}
			hash += (xbar / tol) * coefs[i]
		}
	}
	return int(hash)
}

func Test_up01b(tst *testing.T) {

	// capture errors and flush log
	defer End()

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

	// start simulation
	if !Start("data/up01.sim", true, chk.Verbose) {
		chk.Panic("cannot start simulation")
	}

	// for debugging Kb
	if true {
		defer up_DebugKb(&testKb{
			tst: tst, eid: 3, tol: 1e-8, verb: chk.Verbose,
			ni: 1, nj: 1, itmin: 1, itmax: -1, tmin: 800, tmax: 1000,
		})()
	}

	// run simulation
	if !Run() {
		chk.Panic("cannot run simulation\n")
	}
}

func Test_brent03(tst *testing.T) {

	//verbose()
	chk.PrintTitle("brent03. minimum finding")

	ffcn := func(x float64) (res float64, err error) {
		return x*x*x - 2.0*x - 5.0, nil
	}

	var o Brent
	o.Init(ffcn)
	xa, xb := 0.0, 1.0
	x, err := o.Min(xa, xb, false)
	if err != nil {
		chk.Panic("%v", err)
	}
	y, err := ffcn(x)
	if err != nil {
		chk.Panic("%v", err)
	}
	xcor := math.Sqrt(2.0 / 3.0)
	io.Pforan("x      = %v (correct=%g)\n", x, xcor)
	io.Pforan("f(x)   = %v\n", y)
	io.Pforan("nfeval = %v\n", o.NFeval)
	io.Pforan("nit    = %v\n", o.It)

	//save := true
	save := false
	PlotYxe(ffcn, "results", "brent03.png", x, -1, 3, 101, "Brent", "'b-'", save, false, nil)
	chk.Scalar(tst, "xcorrect", 1e-8, x, xcor)
}

func ltqnorm(p float64) float64 {
	LOW := 0.02425
	HIGH := 0.97575
	var q, r float64
	if p < 0 || p > 1 {
		chk.Panic("input value (%g) is outside allowed range [0, 1]", p)
	} else if p == 0 {
		chk.Panic("p==0 => -infinity")
	} else if p == 1 {
		chk.Panic("p==1 => infinity")
	} else if p < LOW {
		/* Rational approximation for lower region */
		q = math.Sqrt(-2 * math.Log(p))
		return (((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q + c[5]) /
			((((d[0]*q+d[1])*q+d[2])*q+d[3])*q + 1)
	} else if p > HIGH {
		/* Rational approximation for upper region */
		q = math.Sqrt(-2 * math.Log(1-p))
		return -(((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q + c[5]) /
			((((d[0]*q+d[1])*q+d[2])*q+d[3])*q + 1)
	} else {
		/* Rational approximation for central region */
		q = p - 0.5
		r = q * q
		return (((((a[0]*r+a[1])*r+a[2])*r+a[3])*r+a[4])*r + a[5]) * q /
			(((((b[0]*r+b[1])*r+b[2])*r+b[3])*r+b[4])*r + 1)
	}
	return 0
}

// Init initialises this structure
func (o *DynCoefs) Init(dat *inp.SolverData) {

	// hmin
	o.hmin = dat.DtMin

	// HHT
	o.HHT = dat.HHT

	// θ-method
	o.θ = dat.Theta
	if o.θ < 1e-5 || o.θ > 1.0 {
		chk.Panic("θ-method requires 1e-5 <= θ <= 1.0 (θ = %v is incorrect)", o.θ)
	}

	// HHT method
	if dat.HHT {
		o.α = dat.HHTalp
		if o.α < -1.0/3.0 || o.α > 0.0 {
			chk.Panic("HHT method requires: -1/3 <= α <= 0 (α = %v is incorrect)", o.α)
		}
		o.θ1 = (1.0 - 2.0*o.α) / 2.0
		o.θ2 = (1.0 - o.α) * (1.0 - o.α) / 2.0

		// Newmark's method
	} else {
		o.θ1, o.θ2 = dat.Theta1, dat.Theta2
		if o.θ1 < 0.0001 || o.θ1 > 1.0 {
			chk.Panic("θ1 must be between 0.0001 and 1.0 (θ1 = %v is incorrect)", o.θ1)
		}
		if o.θ2 < 0.0001 || o.θ2 > 1.0 {
			chk.Panic("θ2 must be between 0.0001 and 1.0 (θ2 = %v is incorrect)", o.θ2)
		}
	}
}

func Test_basicgeom05(tst *testing.T) {

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

	zero := 1e-1
	told := 1e-9
	tolin := 1e-3

	o := &Point{0, 0, 0}
	a := &Point{0.5, 0.5, 0.5}
	b := &Point{0.1, 0.5, 0.8}
	c := &Point{1, 1, 1}
	cmin, cmax := PointsLims([]*Point{o, a, b, c})
	chk.Vector(tst, "cmin", 1e-17, cmin, []float64{0, 0, 0})
	chk.Vector(tst, "cmax", 1e-17, cmax, []float64{1, 1, 1})
	if !IsPointIn(a, cmin, cmax, zero) {
		chk.Panic("a=%v must be in box")
	}
	if !IsPointIn(b, cmin, cmax, zero) {
		chk.Panic("b=%v must be in box")
	}
	if !IsPointIn(c, cmin, cmax, zero) {
		chk.Panic("c=%v must be in box")
	}

	p := &Point{0.5, 0.5, 0.5}
	q := &Point{1.5, 1.5, 1.5}
	if !IsPointInLine(p, o, c, zero, told, tolin) {
		chk.Panic("p=%v must be in line")
	}
	if IsPointInLine(q, o, c, zero, told, tolin) {
		chk.Panic("q=%v must not be in line")
	}
}

// ReadLines reads lines from a file and calls ReadLinesCallback to process each line being read
func ReadLines(fn string, cb ReadLinesCallback) {
	fil, err := os.Open(os.ExpandEnv(fn))
	if err != nil {
		chk.Panic(_fileio_err05, fn)
	}
	defer fil.Close()
	r := bufio.NewReader(fil)
	idx := 0
	for {
		lin, prefix, errl := r.ReadLine()
		if prefix {
			chk.Panic(_fileio_err06, fn)
		}
		if errl == io.EOF {
			break
		}
		if errl != nil {
			chk.Panic(_fileio_err07, fn)
		}
		stop := cb(idx, string(lin))
		if stop {
			break
		}
		idx++
	}
}

func Test_smpinvs01(tst *testing.T) {

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

	a, b, β, ϵ := -1.0, 0.5, 1e-3, 1e-3

	L := []float64{-8.0, -8.0, -8.0}
	N := make([]float64, 3)
	n := make([]float64, 3)
	m := SmpDirector(N, L, a, b, β, ϵ)
	SmpUnitDirector(n, m, N)
	io.Pforan("L = %v\n", L)
	io.Pforan("N = %v\n", N)
	io.Pforan("m = %v\n", m)
	io.Pforan("n = %v\n", n)
	chk.Vector(tst, "n", 1e-15, n, []float64{a / SQ3, a / SQ3, a / SQ3})

	p, q, err := GenInvs(L, n, a)
	if err != nil {
		chk.Panic("GenInvs failed:\n%v", err.Error())
	}
	io.Pforan("p = %v\n", p)
	io.Pforan("q = %v\n", q)
	if q < 0.0 || q > GENINVSQEPS {
		chk.Panic("q=%g is incorrect", q)
	}
	if math.Abs(p-a*L[0]) > 1e-14 {
		chk.Panic("p=%g is incorrect. err = %g", p, math.Abs(p-a*L[0]))
	}
}

// IntegOnPlane integrates the results of nodes located on a plane perpedicular to either {x,y,z}
// The mesh on the plane must be structured and with rectangle elements in order to build
// an {u,v} coordinates system
// Input:
//   key   -- node variable; e.g. "uz", "ex_sz"
//   alias -- alias of plane; e.g. "surf"
// Output:
//   res -- the result of the integration res = ∫_u ∫_v f(u,v) du dv for all output times
func IntegOnPlane(key, alias string) (res []float64) {

	// get plane
	plane, ok := Planes[alias]
	if !ok {
		chk.Panic("cannot get plane with alias=%q", alias)
	}

	// compute integral
	res = make([]float64, len(TimeInds))
	u := make([]float64, 2)
	for idxI, k := range TimeInds {
		for j := 0; j < plane.Nu[1]; j++ {
			u[1] = plane.Umin[1] + float64(j)*plane.Du[1]
			for i := 0; i < plane.Nu[0]; i++ {
				u[0] = plane.Umin[0] + float64(i)*plane.Du[0]
				id := plane.Ubins.Find(u)
				if id < 0 {
					chk.Panic("IntegOnPlane with key=%q and alias=%q\nError: cannot find {u,v} coordinate in any bin. u=%v vid/ipid=%d", key, alias, u, id)
				}
				vals, ok := plane.id2pt[id].Vals[key]
				if !ok {
					chk.Panic("cannot find results with key=%q for plane with alias=%q", key, alias)
				}
				plane.F[i][j] = vals[idxI]
			}
		}
		res[k] = num.Simps2D(plane.Du[0], plane.Du[1], plane.F)
	}
	return
}

// Define defines aliases
//  alias -- an alias to a group of points, an individual point, or to a set of points.
//           Example: "A", "left-column" or "a b c". If the number of points found is different
//           than the number of aliases, a group is created.
//  Note:
//    To use spaces in aliases, prefix the alias with an exclamation mark; e.g "!right column"
func Define(alias string, loc Locator) {

	// check
	if len(alias) < 1 {
		chk.Panic("alias must have at least one character. %q is invalid", alias)
	}

	// set planes map
	if plane, ok := loc.(*PlaneData); ok {
		Planes[alias] = plane
		defer func() {
			plane.ConnectResults(alias)
		}()
	}

	// locate points
	pts := loc.Locate()
	if len(pts) < 1 {
		chk.Panic("cannot define entities with alias=%q and locator=%v", alias, loc)
	}

	// set results map
	if alias[0] == '!' {
		Results[alias[1:]] = pts
		return
	}
	lbls := strings.Fields(alias)
	if len(lbls) == len(pts) {
		for i, l := range lbls {
			Results[l] = []*Point{pts[i]}
		}
		return
	}
	Results[alias] = pts
}

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

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

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

	// start simulation
	if !Start("data/frees01.sim", true, chk.Verbose) {
		chk.Panic("cannot start FE simulation")
	}

	// domain
	distr := false
	dom := NewDomain(Global.Sim.Regions[0], distr)
	if dom == nil {
		chk.Panic("cannot run FE simulation")
	}

	// set stage
	if !dom.SetStage(0, Global.Sim.Stages[0], distr) {
		chk.Panic("cannot set stage\n")
	}

	// nodes and elements
	chk.IntAssert(len(dom.Nodes), 62)
	chk.IntAssert(len(dom.Elems), 15)

	// vertices with "fl"
	seepverts := map[int]bool{3: true, 45: true, 7: true, 49: true, 11: true, 53: true, 15: true, 57: true, 19: true, 61: true, 23: true}

	// check dofs
	var seepeqs []int
	for _, nod := range dom.Nodes {
		if seepverts[nod.Vert.Id] {
			chk.IntAssert(len(nod.Dofs), 2)
			seepeqs = append(seepeqs, nod.Dofs[1].Eq)
		} else {
			chk.IntAssert(len(nod.Dofs), 1)
		}
	}
	sort.Ints(seepeqs)
	io.Pforan("seepeqs = %v\n", seepeqs)
	chk.Ints(tst, "seepeqs", seepeqs, []int{14, 16, 19, 30, 32, 43, 45, 56, 58, 69, 71})

	// check Fmap
	e2 := dom.Elems[2].(*ElemP)
	chk.Ints(tst, "e2.Fmap", e2.Fmap, []int{14, 16, 19})
	e5 := dom.Elems[5].(*ElemP)
	chk.Ints(tst, "e5.Fmap", e5.Fmap, []int{16, 30, 32})
	e8 := dom.Elems[8].(*ElemP)
	chk.Ints(tst, "e8.Fmap", e8.Fmap, []int{30, 43, 45})
	e11 := dom.Elems[11].(*ElemP)
	chk.Ints(tst, "e11.Fmap", e11.Fmap, []int{43, 56, 58})
	e14 := dom.Elems[14].(*ElemP)
	chk.Ints(tst, "e14.Fmap", e14.Fmap, []int{56, 69, 71})
}

// SetDefault sets default options for report
func NewTexReport(opts []*Optimiser) (o *TexReport) {

	// new struct
	o = new(TexReport)

	// options
	o.DirOut = "/tmp/goga"
	o.Fnkey = "texreport"
	o.Title = "Goga Report"
	o.RefLabel = ""
	o.TextSize = `\scriptsize  \setlength{\tabcolsep}{0.5em}`
	o.MiniPageSz = "4.1cm"
	o.HistTextSize = `\fontsize{5pt}{6pt}`
	o.Type = 4
	o.NRowPerTab = -1
	o.UseGeom = true
	o.RunPDF = true
	o.ShowNsol = false
	o.ShowNcpu = false
	o.ShowTmax = false
	o.ShowDtExc = false
	o.ShowDEC = false
	o.ShowX01 = true

	// constants
	o.DroundCte = 0.001e9 // 0.0001e9

	// input
	o.Opts = opts

	// buffers
	o.buf = new(bytes.Buffer)
	o.binp = new(bytes.Buffer)
	o.bxres = new(bytes.Buffer)

	// check
	if len(o.Opts) < 1 {
		chk.Panic("slice Opts must be set with at least one item")
	}
	if o.Type < 1 || o.Type > 4 {
		chk.Panic("type of table must be in [1,4]")
	}

	// number of samples
	o.nsamples = o.Opts[0].Nsamples

	// symbol for f
	o.symbF = o.Opts[0].RptWordF
	if o.symbF == "" {
		o.symbF = "f"
	}

	// flag
	o.singleObj = o.Opts[0].Nova == 1
	return
}

func Test_out02(tst *testing.T) {

	// finalise analysis process and catch errors
	defer func() {
		if err := recover(); err != nil {
			tst.Fail()
			io.PfRed("ERROR: %v\n", err)
		} else {
			fem.End()
		}
	}()

	// test title
	//verbose()
	chk.PrintTitle("out02")

	// run FE simulation
	if !fem.Start("data/twoqua4.sim", true, chk.Verbose) {
		chk.Panic("cannot start FE simulation")
	}
	if !fem.Run() {
		chk.Panic("cannot run FE simulation")
	}

	// start analysis process
	Start("data/twoqua4.sim", 0, 0)

	// get second ip coordinates
	xip := Ipoints[1].X
	io.Pfcyan("xip = %v\n", xip)

	// define points
	Define("A", N{-1})
	Define("ips", Along{{xip[0], 0}, {xip[0], 1}})

	// load results
	LoadResults(nil)

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

	// check displacements
	tolu := 1e-15
	x := GetCoords("A")
	ux := GetRes("ux", "A", 0)
	uy := GetRes("uy", "A", 0)
	io.Pforan("ux=%v uy=%v\n", ux, uy)
	for j, t := range Times {
		io.Pfyel("t=%g\n", t)
		sol.CheckDispl(tst, t, []float64{ux[j], uy[j]}, x, tolu)
	}
}

// WriteBytesToFile writes bytes to a new file
func WriteBytesToFile(fn string, b []byte) {
	fil, err := os.Create(os.ExpandEnv(fn))
	if err != nil {
		chk.Panic(_fileio_err03, fn)
	}
	defer fil.Close()
	if _, err = fil.Write(b); err != nil {
		chk.Panic(_fileio_err04, err.Error())
	}
}

func Test_bins02(tst *testing.T) {

	//verbose()
	chk.PrintTitle("bins02. find along line (2D)")

	// bins
	var bins Bins
	bins.Init([]float64{-0.2, -0.2}, []float64{0.8, 1.8}, 5)

	// fill bins structure
	maxit := 5 // number of entries
	ID := make([]int, maxit)
	for k := 0; k < maxit; k++ {
		x := float64(k) / float64(maxit)
		ID[k] = k
		err := bins.Append([]float64{x, 2*x + 0.2}, ID[k])
		if err != nil {
			chk.Panic(err.Error())
		}
	}

	// add more points to bins
	for i := 0; i < 5; i++ {
		err := bins.Append([]float64{float64(i) * 0.1, 1.8}, 100+i)
		if err != nil {
			chk.Panic(err.Error())
		}
	}

	// message
	for _, bin := range bins.All {
		if bin != nil {
			io.Pf("%v\n", bin)
		}
	}

	// find points along diagonal
	ids := bins.FindAlongSegment([]float64{0.0, 0.2}, []float64{0.8, 1.8}, 1e-8)
	io.Pforan("ids = %v\n", ids)
	chk.Ints(tst, "ids", ids, ID)

	// find additional points
	ids = bins.FindAlongSegment([]float64{-0.2, 1.8}, []float64{0.8, 1.8}, 1e-8)
	io.Pfcyan("ids = %v\n", ids)
	chk.Ints(tst, "ids", ids, []int{100, 101, 102, 103, 104, 4})

	// draw
	if chk.Verbose {
		plt.SetForPng(1, 500, 150)
		bins.Draw2d(true, true, true, true, map[int]bool{8: true, 9: true, 10: true})
		plt.SetXnticks(15)
		plt.SetYnticks(15)
		plt.SaveD("/tmp/gosl/gm", "test_bins02.png")
	}
}

func RunMumpsTestC(t *TripletC, tol_cmp float64, b, x_correct []complex128, sum_b_to_root bool) {

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

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

	// initialise solver
	err := lis.InitC(t, symmetric, verbose, timing)
	if err != nil {
		chk.Panic("%v", err.Error())
	}

	// factorise
	err = lis.Fact()
	if err != nil {
		chk.Panic("%v", err.Error())
	}

	// solve
	bR, bC := ComplexToRC(b)
	xR := make([]float64, len(b))
	xC := make([]float64, len(b))
	err = lis.SolveC(xR, xC, bR, bC, sum_b_to_root) // x := inv(A) * b
	if err != nil {
		chk.Panic("%v", err.Error())
	}
	x := RCtoComplex(xR, xC)

	if mpi.Rank() == 0 {
		// output
		A := t.ToMatrix(nil)
		io.Pforan("A.x = b\n")
		PrintMatC("A", A.ToDense(), "(%g+", "%gi) ", false)
		PrintVecC("x", x, "(%g+", "%gi) ", false)
		PrintVecC("b", b, "(%g+", "%gi) ", false)

		// check
		xR_correct, xC_correct := ComplexToRC(x_correct)
		errR := VecMaxDiff(xR, xR_correct)
		if errR > tol_cmp {
			chk.Panic("test failed: errR = %g", errR)
		}
		errC := VecMaxDiff(xC, xC_correct)
		if errC > tol_cmp {
			chk.Panic("test failed: errC = %g", errC)
		}
		io.Pf("err(xR) = %g [1;32mOK[0m\n", errR)
		io.Pf("err(xC) = %g [1;32mOK[0m\n", errC)
	}
}

// Read reads configuration parameters from JSON file
func (o *Parameters) Read(filenamepath string) {
	o.Default()
	b, err := io.ReadFile(filenamepath)
	if err != nil {
		chk.Panic("cannot read parameters file %q", filenamepath)
	}
	err = json.Unmarshal(b, o)
	if err != nil {
		chk.Panic("cannot unmarshal parameters file %q", filenamepath)
	}
	return
}