// MatInvG returns the matrix inverse of 'a' in 'ai'. 'a' can be of any size,
// even non-square; in this case, the pseudo-inverse is returned
func MatInvG(ai, a [][]float64, tol float64) (err error) {
	if len(a) < 1 {
		return chk.Err(_matinvg_err1)
	}
	m, n := len(a), len(a[0])
	if m == n && m < 4 { // call simple function
		_, err = MatInv(ai, a, tol)
		return
	}
	am := MatToColMaj(a)
	ami := make([]float64, n*m) // column-major inverse matrix
	if m == n {                 // general matrix inverse
		status := C.lapack_square_inverse((*C.double)(unsafe.Pointer(&ami[0])), (C.long)(m), (*C.double)(unsafe.Pointer(&am[0])))
		if status != 0 {
			return chk.Err(_matinvg_err2, m, n, "lapack_square_inverse", status)
		}
	} else { // pseudo inverse
		status := C.lapack_pseudo_inverse((*C.double)(unsafe.Pointer(&ami[0])), (C.long)(m), (C.long)(n), (*C.double)(unsafe.Pointer(&am[0])), (C.double)(tol))
		if status != 0 {
			return chk.Err(_matinvg_err2, m, n, "lapack_pseudo_inverse", status)
		}
	}
	ColMajToMat(ai, ami)
	return
}

// Init initialises the function
func (o *Halo) Init(prms Prms) (err error) {
	ndim := 2
	for _, p := range prms {
		if p.N == "zc" {
			ndim = 3
			break
		}
	}
	o.xc = make([]float64, ndim)
	e := prms.Connect(&o.r, "r", "halo function")
	e += prms.Connect(&o.xc[0], "xc", "halo function")
	e += prms.Connect(&o.xc[1], "yc", "halo function")
	if ndim == 3 {
		e += prms.Connect(&o.xc[2], "zc", "halo function")
	}
	if e != "" {
		err = chk.Err("%v\n", e)
		return
	}
	rtol := 1e-10
	if o.r < rtol {
		return chk.Err("halo: radius must be greater than %g", rtol)
	}
	return
}

// SaveSol saves solution (o.Sol) to a file which name is set with tidx (time output index)
func (o Domain) SaveSol(tidx int, verbose bool) (err error) {

	// skip if root
	if o.Proc != 0 {
		return
	}

	// buffer and encoder
	var buf bytes.Buffer
	enc := GetEncoder(&buf, o.Sim.EncType)

	// encode Sol
	err = enc.Encode(o.Sol.T)
	if err != nil {
		return chk.Err("cannot encode Domain.Sol.T\n%v", err)
	}
	err = enc.Encode(o.Sol.Y)
	if err != nil {
		return chk.Err("cannot encode Domain.Sol.Y\n%v", err)
	}
	err = enc.Encode(o.Sol.Dydt)
	if err != nil {
		return chk.Err("cannot encode Domain.Sol.Dydt\n%v", err)
	}
	err = enc.Encode(o.Sol.D2ydt2)
	if err != nil {
		return chk.Err("cannot encode Domain.Sol.D2ydt2\n%v", err)
	}

	// save file
	fn := out_nod_path(o.Sim.DirOut, o.Sim.Key, o.Sim.EncType, tidx, o.Proc)
	return save_file(fn, &buf, verbose)
}

// DeleteMany deletes users
func (o *Control) DeleteMany(w http.ResponseWriter, r *http.Request, datUser interface{}) (res *web.Results, err error) {

	// parse data
	if datUser == nil {
		err = chk.Err("user.Control.DeleteMany: invalid datUser = %v", datUser)
		return
	}
	user := datUser.(*User)

	// find users
	users, err := o.find(user, true)
	if err != nil {
		err = chk.Err("user.Control.DeleteMany: cannot find users with datUser = %v", datUser)
		return
	}

	// delete
	for _, c := range users {
		if c == nil {
			continue
		}
		err = o.collection.RemoveId(c.Id)
		if err != nil {
			chk.Err("user.Control.DeleteMany: cannot remove user = %v; err = %v", c, err)
			return
		}
	}
	return
}

// Init initialises the model
func (o *RefDecSp1) Init(prms Prms) (err error) {

	// parameters
	e := prms.Connect(&o.β, "bet", "ref-dec-sp1 function")
	e += prms.Connect(&o.λ1, "lam1", "ref-dec-sp1 function")
	e += prms.Connect(&o.ya, "ya", "ref-dec-sp1 function")
	e += prms.Connect(&o.yb, "yb", "ref-dec-sp1 function")
	if e != "" {
		err = chk.Err("%v\n", e)
		return
	}

	// check
	if o.yb >= o.ya {
		return chk.Err("yb(%g) must be smaller than ya(%g)", o.yb, o.ya)
	}

	// constants
	o.c1 = o.β * o.λ1
	o.c2 = math.Exp(-o.β * o.ya)
	o.c3 = math.Exp(-o.β*o.yb) - o.c2
	o.c1timestmax = 400

	// check
	if math.IsInf(o.c2, 0) || math.IsInf(o.c3, 0) {
		return chk.Err("β*ya or β*yb is too large:\n β=%v, ya=%v, yb=%v\n c1=%v, c2=%v, c3=%v", o.β, o.ya, o.yb, o.c1, o.c2, o.c3)
	}
	return
}

// CalcD computes D = dσ_new/dε_new (consistent)
func (o SmallElasticity) CalcD(D [][]float64, s *State) (err error) {
	if o.Pse {
		if o.Nsig != 4 {
			return chk.Err("for plane-stress analyses, D must be 4x4. nsig = %d is incorrect.\n", o.Nsig)
		}
		if o.Kgc != nil {
			return chk.Err("plane-stress analysis does not work with nonlinear K and G\n")
		}
		c := o.E / (1.0 - o.Nu*o.Nu)
		la.MatFill(D, 0)
		D[0][0] = c
		D[0][1] = c * o.Nu
		D[1][0] = c * o.Nu
		D[1][1] = c
		D[3][3] = c * (1.0 - o.Nu)
		return
	}
	if o.Kgc != nil {
		o.K, o.G = o.Kgc.Calc(s)
	}
	for i := 0; i < o.Nsig; i++ {
		for j := 0; j < o.Nsig; j++ {
			D[i][j] = o.K*tsr.Im[i]*tsr.Im[j] + 2*o.G*tsr.Psd[i][j]
		}
	}
	return
}

// SolveC solves the linear Complex system A.x = b
//  NOTES:
//    1) sum_b_to_root is a flag for MUMPS; it tells Solve to sum the values in 'b' arrays to the root processor
func (o *LinSolMumps) SolveC(xR, xC, bR, bC []float64, sum_b_to_root bool) (err error) {

	// check
	if !o.cmplx {
		return chk.Err(_linsol_mumps_err11)
	}

	// start time
	if o.ton {
		o.tini = time.Now()
	}

	// message
	if o.verb {
		io.Pfgreen("\n . . . . . . . . . . . . . . LinSolMumps.SolveC . . . . . . . . . . . . . . . \n\n")
	}

	// MUMPS: set RHS in processor # 0
	if sum_b_to_root {
		mpi.SumToRoot(xR, bR)
		mpi.SumToRoot(xC, bC)
		// join complex values
		if mpi.Rank() == 0 {
			for i := 0; i < len(xR); i++ {
				o.xRC[i*2], o.xRC[i*2+1] = xR[i], xC[i]
			}
		}
	} else {
		// join complex values
		if mpi.Rank() == 0 {
			for i := 0; i < len(xR); i++ {
				o.xRC[i*2], o.xRC[i*2+1] = bR[i], bC[i]
			}
		}
	}

	// MUMPS: solve
	o.mz.job = 3      // solution code
	C.zmumps_c(&o.mz) // solve
	if o.mz.info[1-1] < 0 {
		return chk.Err(_linsol_mumps_err12, mumps_error(o.mz.info[1-1], o.mz.info[2-1]))
	}

	// MUMPS: split complex values
	if mpi.Rank() == 0 {
		for i := 0; i < len(xR); i++ {
			xR[i], xC[i] = o.xRC[i*2], o.xRC[i*2+1]
		}
	}

	// MUMPS: broadcast from root
	mpi.BcastFromRoot(xR)
	mpi.BcastFromRoot(xC)

	// duration
	if o.ton {
		io.Pfcyan("%s: Time spent in LinSolMumps.Solve = %v\n", o.name, time.Now().Sub(o.tini))
	}
	return
}

// SolveC solves the linear Complex system A.x = b
func (o *LinSolUmfpack) SolveC(xR, xC, bR, bC []float64, dummy bool) (err error) {

	// check
	if !o.cmplx {
		return chk.Err(_linsol_umfpack_err12)
	}

	// start time
	if o.ton {
		o.tini = time.Now()
	}

	// message
	if o.verb {
		io.Pfgreen("\n . . . . . . . . . . . . . . LinSolUmfpack.SolveC . . . . . . . . . . . . . . . \n\n")
	}

	// UMFPACK: pointers
	pxR := (*C.double)(unsafe.Pointer(&xR[0]))
	pxC := (*C.double)(unsafe.Pointer(&xC[0]))
	pbR := (*C.double)(unsafe.Pointer(&bR[0]))
	pbC := (*C.double)(unsafe.Pointer(&bC[0]))

	// UMFPACK: solve
	st := C.umfpack_zl_solve(C.UMFPACK_A, o.ap, o.ai, o.ax, o.az, pxR, pxC, pbR, pbC, o.unum, o.uctrl, nil)
	if st != C.UMFPACK_OK {
		chk.Err(_linsol_umfpack_err13, Uerr2Text[int(st)])
	}

	// duration
	if o.ton {
		io.Pfcyan("%s: Time spent in LinSolUmfpack.Solve = %v\n", o.name, time.Now().Sub(o.tini))
	}
	return
}

// SetInitial sets the initial state
func (o *Domain) SetInitial(stg *inp.Stage) (err error) {

	// check
	if len(stg.Initial.Fcns) != len(stg.Initial.Dofs) {
		return chk.Err("number of functions (fcns) must be equal to number of dofs for setting initial values. %d != %d", len(stg.Initial.Fcns), len(stg.Initial.Dofs))
	}

	// loop over functions
	for i, fname := range stg.Initial.Fcns {

		// get function
		fcn := o.Sim.Functions.Get(fname)
		if fcn == nil {
			return chk.Err("cannot get function named %q", fname)
		}

		// set nodes
		key := stg.Initial.Dofs[i]
		for _, nod := range o.Nodes {
			eq := nod.GetEq(key)
			if eq < 0 {
				return chk.Err("dof=%q cannot be found in node=%d for setting initial values", key, nod.Vert.Id)
			}
			o.Sol.Y[eq] = fcn.F(0, nod.Vert.C)
		}
	}
	return
}

// Delete deletes one user
func (o *Control) Delete(w http.ResponseWriter, r *http.Request, datUser interface{}) (res *web.Results, err error) {

	// parse data
	if datUser == nil {
		err = chk.Err("user.Control.Delete: invalid datUser = %v", datUser)
		return
	}
	user := datUser.(*User)

	// find users
	users, err := o.find(user, false)
	if err != nil {
		err = chk.Err("user.Control.Delete: cannot find user with datUser = %v. err = %v", datUser, err)
		return
	}

	// check
	if users[0] == nil {
		chk.Err("user.Control.Delete: cannot find user to delete. datUser = %v", datUser)
		return
	}

	// delete
	err = o.collection.RemoveId(users[0].Id)
	if err != nil {
		chk.Err("user.Control.Delete: cannot remove user = %v; err = %v", users[0], err)
		return
	}
	return
}

// get_porous_parameters extracts parameters based on region data
func get_porous_parameters(mdb *inp.MatDb, reg *inp.Region, ctag int) (RhoS0, nf0 float64, err error) {
	edat := reg.Etag2data(ctag)
	mat := mdb.Get(edat.Mat)
	if mat.Model != "group" {
		err = chk.Err("geost: material type describing layer must be 'group' with porous data")
		return
	}
	if matname, found := io.Keycode(mat.Extra, "p"); found {
		m := mdb.Get(matname)
		for _, p := range m.Prms {
			switch p.N {
			case "RhoS0":
				RhoS0 = p.V
			case "nf0":
				nf0 = p.V
			}
		}
	}
	if RhoS0 < 1e-7 {
		err = chk.Err("geost: initial density of solids RhoS0=%g is incorrect", RhoS0)
		return
	}
	if nf0 < 1e-7 {
		err = chk.Err("geost: initial porosity nf0=%g is incorrect", nf0)
	}
	return
}

// Init initialises the function
func (o *Cdist) Init(prms Prms) (err error) {
	var xc, yc, zc float64
	is3d := false
	for _, p := range prms {
		switch p.N {
		case "xc":
			xc = p.V
		case "yc":
			yc = p.V
		case "zc":
			zc = p.V
			is3d = true
		case "r":
			o.r = p.V
		default:
			return chk.Err("cdist: parameter named %q is invalid", p.N)
		}
	}
	rtol := 1e-10
	if o.r < rtol {
		return chk.Err("cdist: radius must be greater than %g", rtol)
	}
	if is3d {
		o.xc = []float64{xc, yc, zc}
	} else {
		o.xc = []float64{xc, yc}
	}
	return
}

// ReadSol reads Solution from a file which name is set with tidx (time output index)
func (o *Domain) ReadSol(dir, fnkey, enctype string, tidx int) (err error) {

	// open file
	fn := out_nod_path(dir, fnkey, enctype, tidx, 0) // 0 => reading always from proc # 0
	fil, err := os.Open(fn)
	if err != nil {
		return
	}
	defer func() { err = fil.Close() }()

	// get decoder
	dec := GetDecoder(fil, enctype)

	// decode Sol
	err = dec.Decode(&o.Sol.T)
	if err != nil {
		return chk.Err("cannot decode Domain.Sol.T\n%v", err)
	}
	err = dec.Decode(&o.Sol.Y)
	if err != nil {
		return chk.Err("cannot decode Domain.Sol.Y\n%v", err)
	}
	err = dec.Decode(&o.Sol.Dydt)
	if err != nil {
		return chk.Err("cannot decode Domain.Sol.Dydt\n%v", err)
	}
	err = dec.Decode(&o.Sol.D2ydt2)
	if err != nil {
		return chk.Err("cannot decode Domain.Sol.D2ydt2\n%v", err)
	}
	return
}

// ReadIvs reads elements's internal values from a file which name is set with tidx (time output index)
func (o *Domain) ReadIvs(dir, fnkey, enctype string, tidx, proc int) (err error) {

	// open file
	fn := out_ele_path(dir, fnkey, enctype, tidx, proc)
	fil, err := os.Open(fn)
	if err != nil {
		return
	}
	defer func() { err = fil.Close() }()

	// decoder
	dec := GetDecoder(fil, enctype)

	// elements that are in file
	err = dec.Decode(&o.MyCids)
	if err != nil {
		return chk.Err("cannot decode elements ids:\n%v", err)
	}

	// decode internal variables
	for _, cid := range o.MyCids {
		elem := o.Cid2elem[cid]
		if elem == nil {
			return chk.Err("cannot find element with cid=%d", cid)
		}
		err = elem.Decode(dec)
		if err != nil {
			return chk.Err("cannot decode element:\n%v", err)
		}
	}
	return
}

// Init initialises the model
func (o *RefDecSp1) Init(prms Prms) (err error) {

	// parameters
	for _, p := range prms {
		switch p.N {
		case "bet":
			o.β = p.V
		case "lam1":
			o.λ1 = p.V
		case "ya":
			o.ya = p.V
		case "yb":
			o.yb = p.V
		default:
			return chk.Err("ref-dec-sp1: parameter named %q is invalid", p.N)
		}
	}

	// check
	if o.yb >= o.ya {
		return chk.Err("yb(%g) must be smaller than ya(%g)", o.yb, o.ya)
	}

	// constants
	o.c1 = o.β * o.λ1
	o.c2 = math.Exp(-o.β * o.ya)
	o.c3 = math.Exp(-o.β*o.yb) - o.c2
	o.c1timestmax = 400

	// check
	if math.IsInf(o.c2, 0) || math.IsInf(o.c3, 0) {
		return chk.Err("β*ya or β*yb is too large:\n β=%v, ya=%v, yb=%v\n c1=%v, c2=%v, c3=%v", o.β, o.ya, o.yb, o.c1, o.c2, o.c3)
	}
	return
}

// Append adds a new entry {x, id} to the bins structure
func (o *Bins) Append(x []float64, id int) (err error) {
	idx := o.CalcIdx(x)
	if idx < 0 {
		return chk.Err("point %v is out of range", x)
	}
	bin := o.FindBinByIndex(idx)
	if bin == nil {
		return chk.Err("bin index %v is out of range", idx)
	}
	entry := BinEntry{id, x}
	bin.Entries = append(bin.Entries, &entry)
	return
}

// checkprm checks if parameter value is correct
func checkprm(name string, val, minval, maxval float64, usemin, usemax bool) (err error) {
	if usemin {
		if val < minval {
			return chk.Err("%q parameter: wrong value: %g < %g", name, val, minval)
		}
	}
	if usemax {
		if val > maxval {
			return chk.Err("%q parameter: wrong value: %g > %g", name, val, maxval)
		}
	}
	return
}

// CalcΔεElast calculates Δε corresponding to an elastic loading with Δp and Δq
func CalcΔεElast(Δε []float64, K, G float64, Δp, Δq float64, axsym bool) (Δεv, Δεd float64, err error) {
	Δεv = -Δp / K
	Δεd = Δq / (3.0 * G)
	var Δεx, Δεy, Δεz float64
	if axsym { // axisymmetric
		compression := true
		if compression {
			Δεx = Δεv/3.0 + Δεd/2.0
			Δεy = Δεx
			Δεz = Δεv/3.0 - Δεd
		} else {
			Δεx = Δεv/3.0 - Δεd/2.0
			Δεy = Δεv/3.0 + Δεd
			Δεz = Δεx
		}
	} else { // plane-strain with Δεy = Δεx / 2
		c := 9.0 * Δεd * Δεd / (4.0 * Δεv * Δεv)
		α := 0.0
		if math.Abs(c-1.0) > 1e-15 {
			d := 3.0 * (4.0*c - 1.0)
			if d < 0.0 {
				return 0, 0, chk.Err("discriminant < 0:  c=%v  d=%v", c, d)
			}
			α1 := (1.0 + 2.0*c + math.Sqrt(d)) / (2.0 - 2.0*c)
			α2 := (1.0 + 2.0*c - math.Sqrt(d)) / (2.0 - 2.0*c)
			α = α1
			io.Pfyel("d, α1, α2 = %v, %v, %v\n", d, α1, α2)
		}
		io.Pfyel("c, α = %v, %v\n", c, α)
		Δεy = Δεv / (1.0 + α)
		Δεx = α * Δεy
		Δεz = 0
	}
	//io.Pfpink("Δp=%v, Δq=%v => Δεv=%v, Δεd=%v => Δεx=%v, Δεy=%v, Δεz=%v\n", Δp, Δq, Δεv, Δεd, Δεx, Δεy, Δεz)
	Δε[0] = Δεx
	Δε[1] = Δεy
	Δε[2] = Δεz
	Δε[3] = 0
	Δεv_ := tsr.M_εv(Δε)
	Δεd_ := tsr.M_εd(Δε)
	if math.Abs(Δεv-Δεv_) > 1e-15 {
		return 0, 0, chk.Err(_path_err09, Δεv, Δεv_)
	}
	if Δεd < 0 {
		Δεd_ = -Δεd_ // allow negative values
	}
	if math.Abs(Δεd-Δεd_) > 1e-15 {
		return 0, 0, chk.Err(_path_err10, Δεd, Δεd_)
	}
	return
}

// Init initialises this structure
func (o *SmallElasticity) Init(ndim int, pstress bool, prms fun.Prms) (err error) {
	o.Nsig = 2 * ndim
	o.Pse = pstress
	var has_E, has_ν, has_l, has_G, has_K bool
	for _, p := range prms {
		switch p.N {
		case "E":
			o.E, has_E = p.V, true
		case "nu":
			o.Nu, has_ν = p.V, true
		case "l":
			o.L, has_l = p.V, true
		case "G":
			o.G, has_G = p.V, true
		case "K":
			o.K, has_K = p.V, true
		}
		if skgc, found := io.Keycode(p.Extra, "kgc"); found {
			o.Kgc = GetKgc(skgc, prms)
			if o.Kgc == nil {
				return chk.Err("cannot find kgc model named %s", skgc)
			}
			err = o.Kgc.Init(prms)
			if err != nil {
				return
			}
		}
	}
	switch {
	case has_E && has_ν:
		o.L = Calc_l_from_Enu(o.E, o.Nu)
		o.G = Calc_G_from_Enu(o.E, o.Nu)
		o.K = Calc_K_from_Enu(o.E, o.Nu)
	case has_l && has_G:
		o.E = Calc_E_from_lG(o.L, o.G)
		o.Nu = Calc_nu_from_lG(o.L, o.G)
		o.K = Calc_K_from_lG(o.L, o.G)
	case has_K && has_G:
		o.E = Calc_E_from_KG(o.K, o.G)
		o.Nu = Calc_nu_from_KG(o.K, o.G)
		o.L = Calc_l_from_KG(o.K, o.G)
	case has_K && has_ν:
		o.E = Calc_E_from_Knu(o.K, o.Nu)
		o.G = Calc_G_from_Knu(o.K, o.Nu)
		o.L = Calc_l_from_Knu(o.K, o.Nu)
	default:
		return chk.Err("combination of Elastic constants is incorrect. options are {E,nu}, {l,G}, {K,G} and {K,nu}\n")
	}
	return
}

// SetIniStress sets the initial state with initial stresses
func (o *Domain) SetIniStress(stg *inp.Stage) (err error) {

	// set elements with homogeneous stress state
	dat := stg.IniStress
	if dat.Hom {

		// isotropic state
		if dat.Iso {
			for _, e := range o.ElemIntvars {

				// build map with isotropic and homogeneus state
				coords := e.Ipoints()
				nip := len(coords)
				v := utl.DblVals(nip, dat.S0)
				ivs := map[string][]float64{"sx": v, "sy": v, "sz": v}

				// set element's states
				err = e.SetIniIvs(o.Sol, ivs)
				if err != nil {
					return chk.Err("homogeneous/isotropic: element's internal values setting failed:\n%v", err)
				}
			}
			return
		}

		// plane-strain state
		if dat.Psa {
			sz := dat.Nu * (dat.Sh + dat.Sv)
			for _, e := range o.ElemIntvars {

				// build map with plane-strain and homogeneus state
				coords := e.Ipoints()
				nip := len(coords)
				vx := utl.DblVals(nip, dat.Sh)
				vy := utl.DblVals(nip, dat.Sv)
				vz := utl.DblVals(nip, sz)
				ivs := map[string][]float64{"sx": vx, "sy": vy, "sz": vz}

				// set element's states
				err = e.SetIniIvs(o.Sol, ivs)
				if err != nil {
					return chk.Err("homogeneous/plane-strain: element's internal values setting failed:\n%v", err)
				}
			}
			return
		}
	}
	return
}