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

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

func main() {

	// catch errors
	defer func() {
		if err := recover(); err != nil {
			if mpi.Rank() == 0 {
				chk.Verbose = true
				for i := 8; i > 3; i-- {
					chk.CallerInfo(i)
				}
				io.PfRed("ERROR: %v\n", err)
			}
		}
		mpi.Stop(false)
	}()
	mpi.Start(false)

	// default input parameters

	// read input parameters
	fnamepath, _ := io.ArgToFilename(0, "", ".sim", true)
	verbose := io.ArgToBool(1, true)
	erasePrev := io.ArgToBool(2, true)
	saveSummary := io.ArgToBool(3, true)
	allowParallel := io.ArgToBool(4, true)
	alias := io.ArgToString(5, "")

	// message
	if mpi.Rank() == 0 && verbose {
		io.PfWhite("\nGofem v3 -- Go Finite Element Method\n\n")
		io.Pf("Copyright 2015 Dorival Pedroso and Raul Durand. All rights reserved.\n")
		io.Pf("Use of this source code is governed by a BSD-style\n")
		io.Pf("license that can be found in the LICENSE file.\n\n")

		io.Pf("\n%v\n", io.ArgsTable(
			"filename path", "fnamepath", fnamepath,
			"show messages", "verbose", verbose,
			"erase previous results", "erasePrev", erasePrev,
			"save summary", "saveSummary", saveSummary,
			"allow parallel run", "allowParallel", allowParallel,
			"word to add to results", "alias", alias,
		))
	}

	// profiling?
	defer utl.DoProf(false)()

	// analysis data
	readSummary := false
	analysis := fem.NewFEM(fnamepath, alias, erasePrev, saveSummary, readSummary, allowParallel, verbose, 0)

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

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

// SolveR solves the linear Real 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) SolveR(xR, bR []float64, sum_b_to_root bool) (err error) {

	// check
	if !o.is_initialised {
		return chk.Err("linear solver must be initialised first\n")
	}
	if o.cmplx {
		return chk.Err(_linsol_mumps_err09)
	}

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

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

	// MUMPS: set RHS in processor # 0
	if sum_b_to_root {
		mpi.SumToRoot(xR, bR)
	} else {
		if mpi.Rank() == 0 {
			copy(xR, bR) // x := b
		}
	}

	// only proc # 0 needs the RHS
	if mpi.Rank() == 0 {
		o.m.rhs = (*C.double)(unsafe.Pointer(&xR[0]))
	}

	// MUMPS: solve
	o.m.job = 3      // solution code
	C.dmumps_c(&o.m) // solve
	if o.m.info[1-1] < 0 {
		return chk.Err(_linsol_mumps_err10, mumps_error(o.m.info[1-1], o.m.info[2-1]))
	}
	mpi.BcastFromRoot(xR) // broadcast from root

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

func main() {

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

	if mpi.Rank() == 0 {
		chk.PrintTitle("TestJacobian 02b (MPI)")
	}
	if mpi.Size() > 6 {
		io.Pf("this tests works with 6 or less MPI processors\n")
		return
	}

	ffcn := func(fx, x []float64) error {
		fx[0] = 2.0*x[0] - x[1] + sin(x[2]) - cos(x[3]) - x[5]*x[5] - 1.0      // 0
		fx[1] = -x[0] + 2.0*x[1] + cos(x[2]) - sin(x[3]) + x[5] - 1.0          // 1
		fx[2] = x[0] + 3.0*x[1] + sin(x[3]) - cos(x[4]) - x[5]*x[5] - 1.0      // 2
		fx[3] = 2.0*x[0] + 4.0*x[1] + cos(x[3]) - cos(x[4]) + x[5] - 1.0       // 3
		fx[4] = x[0] + 5.0*x[1] - sin(x[2]) + sin(x[4]) - x[5]*x[5]*x[5] - 1.0 // 4
		fx[5] = x[0] + 6.0*x[1] - cos(x[2]) + cos(x[4]) + x[5] - 1.0           // 5
		return nil
	}
	Jfcn := func(dfdx *la.Triplet, x []float64) error {
		dfdx.Start()
		J := [][]float64{
			{2.0, -1.0, cos(x[2]), sin(x[3]), 0.0, -2.0 * x[5]},
			{-1.0, 2.0, -sin(x[2]), -cos(x[3]), 0.0, 1.0},
			{1.0, 3.0, 0.0, cos(x[3]), sin(x[4]), -2.0 * x[5]},
			{2.0, 4.0, 0.0, -sin(x[3]), sin(x[4]), 1.0},
			{1.0, 5.0, -cos(x[2]), 0.0, cos(x[4]), -3.0 * x[5] * x[5]},
			{1.0, 6.0, sin(x[2]), 0.0, -sin(x[4]), 1.0},
		}
		id, sz, ndim := mpi.Rank(), mpi.Size(), 6
		start, endp1 := (id*ndim)/sz, ((id+1)*ndim)/sz
		for col := 0; col < 6; col++ {
			for row := start; row < endp1; row++ {
				dfdx.Put(row, col, J[row][col])
			}
		}
		//la.PrintMat(fmt.Sprintf("J @ %d",mpi.Rank()), dfdx.ToMatrix(nil).ToDense(), "%12.6f", false)
		return nil
	}
	x := []float64{5.0, 5.0, pi, pi, pi, 5.0}
	var tst testing.T
	num.CompareJac(&tst, ffcn, Jfcn, x, 1e-6, true)
}

func main() {

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

	myrank := mpi.Rank()
	if myrank == 0 {
		chk.PrintTitle("Test MUMPS Sol 05")
	}

	ndim := 10
	id, sz := mpi.Rank(), mpi.Size()
	start, endp1 := (id*ndim)/sz, ((id+1)*ndim)/sz

	if mpi.Size() > ndim {
		chk.Panic("the number of processors must be smaller than or equal to %d", ndim)
	}

	n := 10
	b := make([]complex128, n)
	x_correct := make([]complex128, n)

	// Let exact solution = 1 + 0.5i
	for i := 0; i < ndim; i++ {
		x_correct[i] = complex(float64(i+1), float64(i+1)/10.0)
	}

	var t la.TripletC
	t.Init(ndim, ndim, ndim, true)

	// assemble a and b
	for i := start; i < endp1; i++ {

		// Some very fake diagonals. Should take exactly 20 GMRES steps
		ar := 10.0 + float64(i)/(float64(ndim)/10.0)
		ac := 10.0 - float64(i)/(float64(ndim)/10.0)
		t.Put(i, i, ar, ac)

		// Generate RHS to match exact solution
		b[i] = complex(ar*real(x_correct[i])-ac*imag(x_correct[i]),
			ar*imag(x_correct[i])+ac*real(x_correct[i]))
	}

	sum_b_to_root := true
	la.RunMumpsTestC(&t, 1e-14, b, x_correct, sum_b_to_root)
}

func main() {

	// catch errors
	var tst testing.T
	defer func() {
		if mpi.Rank() == 0 {
			if err := recover(); err != nil {
				io.PfRed("ERROR: %v\n", err)
			}
			if tst.Failed() {
				io.PfRed("test failed\n")
			}
		}
		mpi.Stop(false)
	}()
	mpi.Start(false)

	// start global variables and log
	analysis := fem.NewFEM("data/bh16.sim", "", true, true, false, true, true, 0)

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

	// check
	skipK := true
	tolK := 1e-12
	tolu := 1e-15
	tols := 1e-12
	fem.TestingCompareResultsU(&tst, "data/bh16.sim", "cmp/bh16.cmp", "", tolK, tolu, tols, skipK, true)
}

func main() {

	// catch errors
	var tst testing.T
	defer func() {
		if mpi.Rank() == 0 {
			if err := recover(); err != nil {
				io.PfRed("ERROR: %v\n", err)
			}
			if tst.Failed() {
				io.PfRed("test failed\n")
			}
		}
		mpi.Stop(false)
	}()
	mpi.Start(false)

	// start global variables and log
	analysis := fem.NewFEM("data/p01.sim", "", true, true, false, true, true, 0)

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

func main() {

	// catch errors
	var tst testing.T
	defer func() {
		if mpi.Rank() == 0 {
			if err := recover(); err != nil {
				io.PfRed("ERROR: %v\n", err)
			}
			if tst.Failed() {
				io.PfRed("test failed\n")
			}
		}
		mpi.Stop(false)
	}()
	mpi.Start(false)

	// start global variables and log
	if !fem.Start("data/p01.sim", true, true) {
		tst.Error("Start failed\n")
		return
	}

	// make sure to flush log
	defer fem.End()

	// run simulation
	if !fem.Run() {
		tst.Error("Run failed\n")
		return
	}
}

func (o *Solver) init_mpi() {
	if mpi.IsOn() {
		o.root = (mpi.Rank() == 0)
		if mpi.Size() > 1 {
			o.Distr = true
		}
	}
}

func main() {

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

	if mpi.Rank() == 0 {
		chk.PrintTitle("TestJacobian 01b (MPI)")
	}
	if mpi.Size() != 2 {
		io.Pf("this tests needs MPI 2 processors\n")
		return
	}

	ffcn := func(fx, x []float64) error {
		fx[0] = math.Pow(x[0], 3.0) + x[1] - 1.0
		fx[1] = -x[0] + math.Pow(x[1], 3.0) + 1.0
		return nil
	}
	Jfcn := func(dfdx *la.Triplet, x []float64) error {
		dfdx.Start()
		if false {
			if mpi.Rank() == 0 {
				dfdx.Put(0, 0, 3.0*x[0]*x[0])
				dfdx.Put(1, 0, -1.0)
			} else {
				dfdx.Put(0, 1, 1.0)
				dfdx.Put(1, 1, 3.0*x[1]*x[1])
			}
		} else {
			if mpi.Rank() == 0 {
				dfdx.Put(0, 0, 3.0*x[0]*x[0])
				dfdx.Put(0, 1, 1.0)
			} else {
				dfdx.Put(1, 0, -1.0)
				dfdx.Put(1, 1, 3.0*x[1]*x[1])
			}
		}
		return nil
	}
	x := []float64{0.5, 0.5}
	var tst testing.T
	num.CompareJacMpi(&tst, ffcn, Jfcn, x, 1e-8, true)
}

func main() {

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

	myrank := mpi.Rank()
	if myrank == 0 {
		chk.PrintTitle("Test MUMPS Sol 01a")
	}

	var t la.Triplet
	switch mpi.Size() {
	case 1:
		t.Init(5, 5, 13)
		t.Put(0, 0, 1.0)
		t.Put(0, 0, 1.0)
		t.Put(1, 0, 3.0)
		t.Put(0, 1, 3.0)
		t.Put(2, 1, -1.0)
		t.Put(4, 1, 4.0)
		t.Put(1, 2, 4.0)
		t.Put(2, 2, -3.0)
		t.Put(3, 2, 1.0)
		t.Put(4, 2, 2.0)
		t.Put(2, 3, 2.0)
		t.Put(1, 4, 6.0)
		t.Put(4, 4, 1.0)
	case 2:
		if myrank == 0 {
			t.Init(5, 5, 6)
			t.Put(0, 0, 1.0)
			t.Put(0, 0, 1.0)
			t.Put(1, 0, 3.0)
			t.Put(0, 1, 3.0)
			t.Put(2, 1, -1.0)
			t.Put(4, 1, 4.0)
		} else {
			t.Init(5, 5, 7)
			t.Put(1, 2, 4.0)
			t.Put(2, 2, -3.0)
			t.Put(3, 2, 1.0)
			t.Put(4, 2, 2.0)
			t.Put(2, 3, 2.0)
			t.Put(1, 4, 6.0)
			t.Put(4, 4, 1.0)
		}
	default:
		chk.Panic("this test needs 1 or 2 procs")
	}

	b := []float64{8.0, 45.0, -3.0, 3.0, 19.0}
	x_correct := []float64{1, 2, 3, 4, 5}
	sum_b_to_root := false
	la.RunMumpsTestR(&t, 1e-14, b, x_correct, sum_b_to_root)
}

func setslice(x []float64) {
	switch mpi.Rank() {
	case 0:
		copy(x, []float64{0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3})
	case 1:
		copy(x, []float64{10, 10, 10, 20, 20, 20, 30, 30, 30, 40, 40})
	case 2:
		copy(x, []float64{100, 100, 100, 1000, 1000, 1000, 2000, 2000, 2000, 3000, 3000})
	}
}

func main() {

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

	myrank := mpi.Rank()
	if myrank == 0 {
		chk.PrintTitle("Test MUMPS Sol 04")
	}

	ndim := 10
	id, sz := mpi.Rank(), mpi.Size()
	start, endp1 := (id*ndim)/sz, ((id+1)*ndim)/sz

	if mpi.Size() > ndim {
		chk.Panic("the number of processors must be smaller than or equal to %d", ndim)
	}

	b := make([]complex128, ndim)
	var t la.TripletC
	t.Init(ndim, ndim, ndim*ndim, true)

	for i := start; i < endp1; i++ {
		j := i
		if i > 0 {
			j = i - 1
		}
		for ; j < 10; j++ {
			val := 10.0 - float64(j)
			if i > j {
				val -= 1.0
			}
			t.Put(i, j, val, 0)
		}
		b[i] = complex(float64(i+1), 0.0)
	}

	x_correct := []complex128{-1, 8, -65, 454, -2725, 13624, -54497, 163490, -326981, 326991}
	sum_b_to_root := true
	la.RunMumpsTestC(&t, 1e-4, b, x_correct, sum_b_to_root)
}

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

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

// Start initialises 'global' and starts logging
func Start(simfilepath string, erasefiles, verbose bool) (startisok bool) {

	// multiprocessing data
	Global.Rank = 0
	Global.Nproc = 1
	Global.Root = true
	Global.Distr = false
	if mpi.IsOn() {
		Global.Rank = mpi.Rank()
		Global.Nproc = mpi.Size()
		Global.Root = Global.Rank == 0
		Global.Distr = Global.Nproc > 1
	}
	Global.Verbose = verbose
	if !Global.Root {
		Global.Verbose = false
	}
	Global.WspcStop = make([]int, Global.Nproc)
	Global.WspcInum = make([]int, Global.Nproc)

	// simulation and convenience variables
	dir := filepath.Dir(simfilepath)
	fn := filepath.Base(simfilepath)
	Global.Sim = inp.ReadSim(dir, fn, Global.LogPrefix, erasefiles)
	LogErrCond(Global.Sim == nil, "ReadSim failed\n")
	if Stop() {
		return
	}
	Global.Ndim = Global.Sim.Ndim
	Global.Dirout = Global.Sim.Data.DirOut
	Global.Fnkey = Global.Sim.Data.FnameKey
	Global.Enc = Global.Sim.Data.Encoder
	Global.Stat = Global.Sim.Data.Stat
	Global.LogBcs = Global.Sim.Data.LogBcs
	Global.Debug = Global.Sim.Data.Debug

	// fix show residual flag
	if !Global.Root {
		Global.Sim.Data.ShowR = false
	}

	// auxiliar structures
	Global.DynCoefs = new(DynCoefs)
	if !Global.DynCoefs.Init(&Global.Sim.Solver) {
		return
	}
	Global.HydroSt = new(HydroStatic)
	Global.HydroSt.Init()

	// success
	return true
}

/*  Jacobian
    ========
        Calculates (with N=n-1):
            df0dx0, df0dx1, df0dx2, ... df0dxN
            df1dx0, df1dx1, df1dx2, ... df1dxN
                 . . . . . . . . . . . . .
            dfNdx0, dfNdx1, dfNdx2, ... dfNdxN
    INPUT:
        ffcn : f(x) function
        x    : station where dfdx has to be calculated
        fx   : f @ x
        w    : workspace with size == n == len(x)
    RETURNS:
        J : dfdx @ x [must be pre-allocated]        */
func Jacobian(J *la.Triplet, ffcn Cb_f, x, fx, w []float64, distr bool) (err error) {
	ndim := len(x)
	start, endp1 := 0, ndim
	if distr {
		id, sz := mpi.Rank(), mpi.Size()
		start, endp1 = (id*ndim)/sz, ((id+1)*ndim)/sz
		if J.Max() == 0 {
			J.Init(ndim, ndim, (endp1-start)*ndim)
		}
	} else {
		if J.Max() == 0 {
			J.Init(ndim, ndim, ndim*ndim)
		}
	}
	J.Start()
	// NOTE: cannot split calculation by columns unless the f function is
	//       independently calculated by each MPI processor.
	//       Otherwise, the AllReduce in f calculation would
	//       join pieces of f from different processors calculated for
	//       different x values (δx[col] from different columns).
	/*
	   for col := start; col < endp1; col++ {
	       xsafe := x[col]
	       delta := math.Sqrt(EPS * max(CTE1, math.Abs(xsafe)))
	       x[col] = xsafe + delta
	       ffcn(w, x) // fnew
	       io.Pforan("x = %v, f = %v\n", x, w)
	       for row := 0; row < ndim; row++ {
	           J.Put(row, col, (w[row]-fx[row])/delta)
	       }
	       x[col] = xsafe
	   }
	*/
	var df float64
	for col := 0; col < ndim; col++ {
		xsafe := x[col]
		delta := math.Sqrt(EPS * max(CTE1, math.Abs(xsafe)))
		x[col] = xsafe + delta
		err = ffcn(w, x) // w := f(x+δx[col])
		if err != nil {
			return
		}
		for row := start; row < endp1; row++ {
			df = w[row] - fx[row]
			//if math.Abs(df) > EPS {
			J.Put(row, col, df/delta)
			//}
		}
		x[col] = xsafe
	}
	return
}

// InitLogFile initialises logger
func InitLogFile(dirout, fnamekey string) (err error) {

	// create log file
	var rank int
	if mpi.IsOn() {
		rank = mpi.Rank()
	}
	LogFile, err = os.Create(io.Sf("%s/%s_p%d.log", dirout, fnamekey, rank))
	if err != nil {
		return
	}

	// connect logger to output file
	log.SetOutput(LogFile)
	return
}