func _TestMultMV(t *testing.T) {
	bM := 100 * M
	bN := 100 * N
	A := matrix.FloatNormal(bM, bN)
	X := matrix.FloatNormal(bN, 1)
	Y1 := matrix.FloatZeros(bM, 1)
	Y0 := matrix.FloatZeros(bM, 1)

	Ar := A.FloatArray()
	Xr := X.FloatArray()
	Y1r := Y1.FloatArray()

	blas.GemvFloat(A, X, Y0, 1.0, 1.0)

	DMultMV(Y1r, Ar, Xr, 1.0, 1.0, NOTRANS, 1, A.LeadingIndex(), 1, 0, bN, 0, bM, 32, 32)
	t.Logf("Y0 == Y1: %v\n", Y0.AllClose(Y1))
	/*
	   if ! Y0.AllClose(Y1) {
	       y0 := Y0.SubMatrix(0, 0, 2, 1)
	       y1 := Y1.SubMatrix(0, 0, 2, 1)
	       t.Logf("y0=\n%v\n", y0)
	       t.Logf("y1=\n%v\n", y1)
	   }
	*/
}

func runTest(A *matrix.FloatMatrix, ntest, LB int) time.Duration {
	var W *matrix.FloatMatrix = nil
	var mintime time.Duration

	N := A.Cols()
	tau := matrix.FloatZeros(N, 1)
	if LB > 0 {
		W = matrix.FloatZeros(A.Rows(), LB)
	}
	fnc := func() {
		_, ERRmatops = matops.DecomposeQR(A, tau, W, LB)
	}

	A0 := A.Copy()
	for n := 0; n < ntest; n++ {
		if n > 0 {
			// restore original A
			A0.CopyTo(A)
			tau.Scale(0.0)
		}
		mperf.FlushCache()
		time0 := mperf.Timeit(fnc)
		if n == 0 || time0 < mintime {
			mintime = time0
		}
		if verbose {
			fmt.Printf("%.4f ms\n", time0.Seconds()*1000.0)
		}
	}
	return mintime
}

func TestMultQT(t *testing.T) {
	M := 60
	N := 40
	K := 30
	nb := 12
	A := matrix.FloatUniform(M, N)
	B := matrix.FloatUniform(M, K)
	W := matrix.FloatZeros(N, nb)
	T := matrix.FloatZeros(N, N)

	// QR = Q*R
	QR, err := DecomposeQRT(A.Copy(), T, W, nb)
	if err != nil {
		t.Logf("decompose-err: %v\n", err)
	}

	// compute: B - Q*Q.T*B = 0

	// X = Q*Q.T*B
	X := B.Copy()
	MultQT(X, QR, T, W, LEFT|TRANS, nb)
	MultQT(X, QR, T, W, LEFT, nb)
	B.Minus(X)

	// ||B - Q*Q.T*B||_1
	nrm := NormP(B, NORM_ONE)
	t.Logf("||B - Q*Q.T*B||_1: %e\n", nrm)
}

func _TestMultMVTransA(t *testing.T) {
	bM := 1000 * M
	bN := 1000 * N
	A := matrix.FloatNormal(bN, bM)
	X := matrix.FloatWithValue(bN, 1, 1.0)
	Y1 := matrix.FloatZeros(bM, 1)
	Y0 := matrix.FloatZeros(bM, 1)

	Ar := A.FloatArray()
	Xr := X.FloatArray()
	Y1r := Y1.FloatArray()

	blas.GemvFloat(A, X, Y0, 1.0, 1.0, linalg.OptTrans)

	DMultMV(Y1r, Ar, Xr, 1.0, 1.0, TRANSA, 1, A.LeadingIndex(), 1, 0, bN, 0, bM, 4, 4)
	ok := Y0.AllClose(Y1)
	t.Logf("Y0 == Y1: %v\n", ok)
	if !ok {
		var y1, y0 matrix.FloatMatrix
		Y1.SubMatrix(&y1, 0, 0, 5, 1)
		t.Logf("Y1[0:5]:\n%v\n", y1)
		Y0.SubMatrix(&y0, 0, 0, 5, 1)
		t.Logf("Y0[0:5]:\n%v\n", y0)
	}
}

func _TestMultSymmLowerSmall(t *testing.T) {
	//bM := 5
	bN := 7
	bP := 7
	Adata := [][]float64{
		[]float64{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
		[]float64{1.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0},
		[]float64{1.0, 2.0, 3.0, 0.0, 0.0, 0.0, 0.0},
		[]float64{1.0, 2.0, 3.0, 4.0, 0.0, 0.0, 0.0},
		[]float64{1.0, 2.0, 3.0, 4.0, 5.0, 0.0, 0.0},
		[]float64{1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 0.0},
		[]float64{1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0}}

	A := matrix.FloatMatrixFromTable(Adata, matrix.RowOrder)
	B := matrix.FloatNormal(bN, bP)
	C0 := matrix.FloatZeros(bN, bP)
	C1 := matrix.FloatZeros(bN, bP)

	Ar := A.FloatArray()
	Br := B.FloatArray()
	C1r := C1.FloatArray()

	blas.SymmFloat(A, B, C0, 1.0, 1.0, linalg.OptLower, linalg.OptRight)

	DMultSymm(C1r, Ar, Br, 1.0, 1.0, LOWER|RIGHT, bN, A.LeadingIndex(), bN,
		bN, 0, bP, 0, bN, 2, 2, 2)
	ok := C0.AllClose(C1)
	t.Logf("C0 == C1: %v\n", ok)
	if !ok {
		t.Logf("A=\n%v\n", A)
		t.Logf("blas: C=A*B\n%v\n", C0)
		t.Logf("C1: C1 = A*X\n%v\n", C1)
	}
}

func TestSolveLeastSquaresQRT(t *testing.T) {
	M := 60
	N := 40
	K := 30
	nb := 12

	A := matrix.FloatUniform(M, N)
	B := matrix.FloatZeros(M, K)
	X0 := matrix.FloatUniform(N, K)

	// B = A*X0
	Mult(B, A, X0, 1.0, 0.0, NOTRANS)
	W := matrix.FloatZeros(N, nb)
	T := matrix.FloatZeros(N, N)

	QR, err := DecomposeQRT(A, T, W, nb)
	if err != nil {
		t.Logf("decompose error: %v\n", err)
	}
	// B' = A.-1*B
	err = SolveQRT(B, QR, T, W, NOTRANS, nb)

	// expect B[0:N, 0:K] == X0, B[N:M, 0:K] == 0.0
	var Xref matrix.FloatMatrix
	Bref := matrix.FloatZeros(M, K)
	Bref.SubMatrix(&Xref, 0, 0, N, K)
	Xref.Plus(X0)
	Bref.Minus(B)
	t.Logf("\nmin ||B - A*X0||\n\twhere B = A*X0\n")
	t.Logf("||B - A*X0||_1 ~~ 0.0: %e\n", NormP(Bref, NORM_ONE))
}

func TestUpdateTrmMV(t *testing.T) {
	//bM := 5
	bN := 8
	//bP := 4
	nb := 4
	X := matrix.FloatNormal(bN, 1)
	//B := matrix.FloatNormal(bP, bN)
	Y := X.Copy()
	C0 := matrix.FloatZeros(bN, bN)
	C2 := matrix.FloatZeros(bN, bN)
	C1 := matrix.FloatZeros(bN, bN)

	Xr := X.FloatArray()
	Yr := Y.FloatArray()
	C1r := C1.FloatArray()
	C0r := C0.FloatArray()
	C2r := C2.FloatArray()

	// no transpose
	DRankMV(C1r, Xr, Yr, 1.0, C1.LeadingIndex(), 1, 1,
		0, bN, 0, bN, nb, nb)
	DTrmUpdMV(C0r, Xr, Yr, 1.0, LOWER, C0.LeadingIndex(), 1, 1,
		0, bN, nb)
	DTrmUpdMV(C2r, Xr, Yr, 1.0, UPPER, C2.LeadingIndex(), 1, 1,
		0, bN, nb)

	t.Logf("C1:\n%v\nC0:\n%v\nC2:\n%v\n", C1, C0, C2)
	// C0 == C2.T
	t.Logf("C0 == C2.T: %v\n", C0.AllClose(C2.Transpose()))
	// C1 == C1 - C2 + C0.T
	Cn := matrix.Minus(C1, C2)
	Cn.Plus(C0.Transpose())
	t.Logf("C1 == C1 - C2 + C0.T: %v\n", Cn.AllClose(C1))

}

func TestQRSmal(t *testing.T) {
	data := [][]float64{
		[]float64{12.0, -51.0, 4.0},
		[]float64{6.0, 167.0, -68.0},
		[]float64{-4.0, 24.0, -41.0}}

	A := matrix.FloatMatrixFromTable(data, matrix.RowOrder)
	T := matrix.FloatZeros(A.Cols(), A.Cols())
	T0 := T.Copy()

	M := A.Rows()
	//N := A.Cols()
	Tau := matrix.FloatZeros(M, 1)
	X, _ := DecomposeQR(A.Copy(), Tau, nil, 0)
	t.Logf("A\n%v\n", A)
	t.Logf("X\n%v\n", X)
	t.Logf("Tau\n%v\n", Tau)

	Tau0 := matrix.FloatZeros(M, 1)
	lapack.Geqrf(A, Tau0)
	t.Logf("lapack X\n%v\n", A)
	t.Logf("lapack Tau\n%v\n", Tau0)

	unblkQRBlockReflector(X, Tau, T)
	t.Logf("T:\n%v\n", T)

	V := TriLU(X.Copy())
	lapack.LarftFloat(V, Tau, T0)
	t.Logf("T0:\n%v\n", T0)

}

func (gp *gpConvexProg) F1(x *matrix.FloatMatrix) (f, Df *matrix.FloatMatrix, err error) {
	f = nil
	Df = nil
	err = nil
	f = matrix.FloatZeros(gp.mnl+1, 1)
	Df = matrix.FloatZeros(gp.mnl+1, gp.n)
	y := gp.g.Copy()
	blas.GemvFloat(gp.F, x, y, 1.0, 1.0)

	for i, s := range gp.ind {
		start := s[0]
		stop := s[1]
		// yi := exp(yi) = exp(Fi*x+gi)
		ymax := maxvec(y.FloatArray()[start:stop])
		// ynew = exp(y[start:stop] - ymax)
		ynew := matrix.Exp(matrix.FloatVector(y.FloatArray()[start:stop]).Add(-ymax))
		y.SetIndexesFromArray(ynew.FloatArray(), matrix.Indexes(start, stop)...)

		// fi = log sum yi = log sum exp(Fi*x+gi)
		ysum := blas.AsumFloat(y, &la.IOpt{"n", stop - start}, &la.IOpt{"offset", start})
		f.SetIndex(i, ymax+math.Log(ysum))

		blas.ScalFloat(y, 1.0/ysum, &la.IOpt{"n", stop - start}, &la.IOpt{"offset", start})
		blas.GemvFloat(gp.F, y, Df, 1.0, 0.0, la.OptTrans, &la.IOpt{"m", stop - start},
			&la.IOpt{"incy", gp.mnl + 1}, &la.IOpt{"offseta", start},
			&la.IOpt{"offsetx", start}, &la.IOpt{"offsety", i})
	}
	return
}

func main() {
	flag.Parse()
	if len(spPath) > 0 {
		checkpnt.Reset(spPath)
		checkpnt.Activate()
		checkpnt.Verbose(spVerbose)
		checkpnt.Format("%.17f")
	}

	adata := [][]float64{
		[]float64{0.3, -0.4, -0.2, -0.4, 1.3},
		[]float64{0.6, 1.2, -1.7, 0.3, -0.3},
		[]float64{-0.3, 0.0, 0.6, -1.2, -2.0}}

	A := matrix.FloatMatrixFromTable(adata, matrix.ColumnOrder)
	b := matrix.FloatVector([]float64{1.5, 0.0, -1.2, -0.7, 0.0})

	_, n := A.Size()
	N := n + 1 + n

	h := matrix.FloatZeros(N, 1)
	h.SetIndex(n, 1.0)

	I0 := matrix.FloatDiagonal(n, -1.0)
	I1 := matrix.FloatIdentity(n)
	G, _ := matrix.FloatMatrixStacked(matrix.StackDown, I0, matrix.FloatZeros(1, n), I1)

	At := A.Transpose()
	P := At.Times(A)
	q := At.Times(b).Scale(-1.0)

	dims := sets.NewDimensionSet("l", "q", "s")
	dims.Set("l", []int{n})
	dims.Set("q", []int{n + 1})

	var solopts cvx.SolverOptions
	solopts.MaxIter = 20
	solopts.ShowProgress = true
	if maxIter > 0 {
		solopts.MaxIter = maxIter
	}
	if len(solver) > 0 {
		solopts.KKTSolverName = solver
	}
	sol, err := cvx.ConeQp(P, q, G, h, nil, nil, dims, &solopts, nil)
	if err == nil {
		x := sol.Result.At("x")[0]
		s := sol.Result.At("s")[0]
		z := sol.Result.At("z")[0]
		fmt.Printf("Optimal\n")
		fmt.Printf("x=\n%v\n", x.ToString("%.9f"))
		fmt.Printf("s=\n%v\n", s.ToString("%.9f"))
		fmt.Printf("z=\n%v\n", z.ToString("%.9f"))
		check(x, s, z)
	}

}

//
// Solves a geometric program
//
//   minimize    log sum exp (F0*x+g0)
//   subject to  log sum exp (Fi*x+gi) <= 0,  i=1,...,m
//               G*x <= h
//               A*x = b
//
func Gp(K []int, F, g, G, h, A, b *matrix.FloatMatrix, solopts *SolverOptions) (sol *Solution, err error) {

	if err = checkArgK(K); err != nil {
		return
	}
	l := sumdim(K)

	if F == nil || F.Rows() != l {
		err = errors.New(fmt.Sprintf("'F' must matrix with %d rows", l))
		return
	}

	if g == nil || !g.SizeMatch(l, 1) {
		err = errors.New(fmt.Sprintf("'g' must matrix with size (%d,1)", l))
		return
	}
	n := F.Cols()

	if G == nil {
		G = matrix.FloatZeros(0, n)
	}
	if h == nil {
		h = matrix.FloatZeros(0, 1)
	}
	if G.Cols() != n {
		err = errors.New(fmt.Sprintf("'G' must matrix with size %d columns", n))
		return
	}
	ml := G.Rows()
	if h == nil || !h.SizeMatch(ml, 1) {
		err = errors.New(fmt.Sprintf("'h' must matrix with size (%d,1)", ml))
		return
	}

	if A == nil {
		A = matrix.FloatZeros(0, n)
	}
	if b == nil {
		b = matrix.FloatZeros(0, 1)
	}
	if A.Cols() != n {
		err = errors.New(fmt.Sprintf("'A' must matrix with size %d columns", n))
		return
	}
	p := A.Rows()
	if b == nil || !b.SizeMatch(p, 1) {
		err = errors.New(fmt.Sprintf("'b' must matrix with size (%d,1)", p))
		return
	}

	dims := sets.NewDimensionSet("l", "q", "s")
	dims.Set("l", []int{ml})
	gpProg := createGpProg(K, F, g)

	return Cp(gpProg, G, h, A, b, dims, solopts)
}

// Returns min {t | x + t*e >= 0}, where e is defined as follows
//
//  - For the nonlinear and 'l' blocks: e is the vector of ones.
//  - For the 'q' blocks: e is the first unit vector.
//  - For the 's' blocks: e is the identity matrix.
//
// When called with the argument sigma, also returns the eigenvalues
// (in sigma) and the eigenvectors (in x) of the 's' components of x.
func maxStep(x *matrix.FloatMatrix, dims *sets.DimensionSet, mnl int, sigma *matrix.FloatMatrix) (rval float64, err error) {
	/*DEBUGGED*/

	rval = 0.0
	err = nil
	t := make([]float64, 0, 10)
	ind := mnl + dims.Sum("l")
	if ind > 0 {
		t = append(t, -minvec(x.FloatArray()[:ind]))
	}
	for _, m := range dims.At("q") {
		if m > 0 {
			v := blas.Nrm2Float(x, &la_.IOpt{"offset", ind + 1}, &la_.IOpt{"n", m - 1})
			v -= x.GetIndex(ind)
			t = append(t, v)
		}
		ind += m
	}

	//var Q *matrix.FloatMatrix
	//var w *matrix.FloatMatrix
	ind2 := 0
	//if sigma == nil && len(dims.At("s")) > 0 {
	//	mx := dims.Max("s")
	//	Q = matrix.FloatZeros(mx, mx)
	//	w = matrix.FloatZeros(mx, 1)
	//}
	for _, m := range dims.At("s") {
		if sigma == nil {
			Q := matrix.FloatZeros(m, m)
			w := matrix.FloatZeros(m, 1)
			blas.Copy(x, Q, &la_.IOpt{"offsetx", ind}, &la_.IOpt{"n", m * m})
			err = lapack.SyevrFloat(Q, w, nil, 0.0, nil, []int{1, 1}, la_.OptRangeInt,
				&la_.IOpt{"n", m}, &la_.IOpt{"lda", m})
			if m > 0 && err == nil {
				t = append(t, -w.GetIndex(0))
			}
		} else {
			err = lapack.SyevdFloat(x, sigma, la_.OptJobZValue, &la_.IOpt{"n", m},
				&la_.IOpt{"lda", m}, &la_.IOpt{"offseta", ind}, &la_.IOpt{"offsetw", ind2})
			if m > 0 {
				t = append(t, -sigma.GetIndex(ind2))
			}
		}
		ind += m * m
		ind2 += m
	}

	if len(t) > 0 {
		rval = maxvec(t)
	}
	return
}

// single invocation for matops and lapack functions
func runCheck(A *matrix.FloatMatrix, LB int) (bool, time.Duration, time.Duration) {

	var W *matrix.FloatMatrix = nil
	N := A.Cols()
	tau := matrix.FloatZeros(N, 1)
	if LB > 0 {
		W = matrix.FloatZeros(A.Rows(), LB)
	}
	fnc := func() {
		_, ERRmatops = matops.DecomposeQR(A, tau, W, LB)
	}

	if verbose && N < 10 {
		fmt.Fprintf(os.Stderr, "A start:\n%v\n", A)
	}
	A0 := A.Copy()
	tau0 := tau.Copy()

	mperf.FlushCache()
	time0 := mperf.Timeit(fnc)
	if verbose && N < 10 {
		fmt.Fprintf(os.Stderr, "A end:\n%v\n", A)
		tau.SetSize(1, N, 1)
		fmt.Fprintf(os.Stderr, "tau: %v\n", tau)
	}

	fn2 := func() {
		ERRlapack = lapack.Geqrf(A0, tau0)
	}

	mperf.FlushCache()
	time2 := mperf.Timeit(fn2)
	if verbose && N < 10 {
		fmt.Fprintf(os.Stderr, "A0 end:\n%v\n", A0)
		tau0.SetSize(1, N, 1) // row vector
		fmt.Fprintf(os.Stderr, "tau0: %v\n", tau0)
	}
	// now A == A0 && tau == tau0

	ok := A.AllClose(A0)
	oktau := tau.AllClose(tau0)
	if !ok || !oktau {
		// save result to globals
		Rlapack = A0
		Rmatops = A
		TAUlapack = tau0
		TAUmatops = tau
	}
	return ok && oktau, time0, time2
}

func _TestBK2U(t *testing.T) {
	Bdata := [][]float64{
		[]float64{10.0, 20.0},
		[]float64{10.0, 20.0},
		[]float64{10.0, 20.0},
		[]float64{10.0, 20.0},
		[]float64{10.0, 20.0},
		[]float64{10.0, 20.0},
		[]float64{10.0, 20.0}}

	N := 7

	A0 := matrix.FloatNormal(N, N)
	A := matrix.FloatZeros(N, N)
	// A is symmetric, posivite definite
	Mult(A, A0, A0, 1.0, 1.0, TRANSB)

	X := matrix.FloatMatrixFromTable(Bdata, matrix.RowOrder)
	B := matrix.FloatZeros(N, 2)
	MultSym(B, A, X, 1.0, 0.0, LOWER|LEFT)
	t.Logf("initial B:\n%v\n", B)

	nb := 0
	W := matrix.FloatWithValue(A.Rows(), 5, 1.0)
	A.SetAt(4, 1, A.GetAt(4, 1)+1.0)
	A.SetAt(1, 4, A.GetAt(4, 1))

	ipiv := make([]int, N, N)
	L, _ := DecomposeBK(A.Copy(), W, ipiv, LOWER, nb)
	t.Logf("ipiv: %v\n", ipiv)
	t.Logf("L:\n%v\n", L)

	ipiv0 := make([]int, N, N)
	nb = 4
	L0, _ := DecomposeBK(A.Copy(), W, ipiv0, LOWER, nb)
	t.Logf("ipiv: %v\n", ipiv0)
	t.Logf("L:\n%v\n", L0)
	B0 := B.Copy()
	SolveBK(B0, L0, ipiv0, LOWER)
	t.Logf("B0:\n%v\n", B0)

	ipiv2 := make([]int32, N, N)
	lapack.Sytrf(A, ipiv2, linalg.OptLower)
	t.Logf("ipiv2: %v\n", ipiv2)
	t.Logf("lapack A:\n%v\n", A)
	lapack.Sytrs(A, B, ipiv2, linalg.OptLower)
	t.Logf("lapack B:\n%v\n", B)
	t.Logf("B == B0: %v\n", B.AllClose(B0))
}

// In-place version of pack(), which also accepts matrix arguments x.
// The columns of x are elements of S, with the 's' components stored
// in unpacked storage.  On return, the 's' components are stored in
// packed storage and the off-diagonal entries are scaled by sqrt(2).
//
func pack2(x *matrix.FloatMatrix, dims *sets.DimensionSet, mnl int) (err error) {
	if len(dims.At("s")) == 0 {
		return nil
	}

	const sqrt2 = 1.41421356237309504880

	iu := mnl + dims.Sum("l", "q")
	ip := iu
	row := matrix.FloatZeros(1, x.Cols())
	//fmt.Printf("x.size = %d %d\n", x.Rows(), x.Cols())
	for _, n := range dims.At("s") {
		for k := 0; k < n; k++ {
			cnt := n - k
			row = x.GetRow(iu+(n+1)*k, row)
			//fmt.Printf("%02d: %v\n", iu+(n+1)*k, x.FloatArray())
			x.SetRow(ip, row)
			for i := 1; i < n-k; i++ {
				row = x.GetRow(iu+(n+1)*k+i, row)
				//fmt.Printf("%02d: %v\n", iu+(n+1)*k+i, x.FloatArray())
				x.SetRow(ip+i, row.Scale(sqrt2))
			}
			ip += cnt
		}
		iu += n * n
	}
	return nil
}

func _TestViewUpdate(t *testing.T) {
	Adata2 := [][]float64{
		[]float64{4.0, 2.0, 2.0},
		[]float64{6.0, 4.0, 2.0},
		[]float64{4.0, 6.0, 1.0},
	}

	A := matrix.FloatMatrixFromTable(Adata2, matrix.RowOrder)
	N := A.Rows()

	// simple LU decomposition without pivoting
	var A11, a10, a01, a00 matrix.FloatMatrix
	for k := 1; k < N; k++ {
		a00.SubMatrixOf(A, k-1, k-1, 1, 1)
		a01.SubMatrixOf(A, k-1, k, 1, A.Cols()-k)
		a10.SubMatrixOf(A, k, k-1, A.Rows()-k, 1)
		A11.SubMatrixOf(A, k, k)
		//t.Logf("A11: %v  a01: %v\n", A11, a01)
		a10.Scale(1.0 / a00.Float())
		MVRankUpdate(&A11, &a10, &a01, -1.0)
	}

	Ld := TriLU(A.Copy())
	Ud := TriU(A)
	t.Logf("Ld:\n%v\nUd:\n%v\n", Ld, Ud)
	An := matrix.FloatZeros(N, N)
	Mult(An, Ld, Ud, 1.0, 1.0, NOTRANS)
	t.Logf("A == Ld*Ud: %v\n", An.AllClose(An))
}

func MakeData(M, N, P int, randomData, diagonal bool) (A, B, C *matrix.FloatMatrix) {

	if diagonal && M != N {
		diagonal = false
		fmt.Printf("cannot make diagonal if B.rows != B.cols\n")
	}

	if randomData {
		A = matrix.FloatNormal(M, P)
		if diagonal {
			d := matrix.FloatNormal(P, 1)
			B := matrix.FloatDiagonal(P, 0.0)
			B.SetIndexesFromArray(d.FloatArray(), matrix.DiagonalIndexes(B)...)
		} else {
			B = matrix.FloatNormal(P, N)
		}
	} else {
		A = matrix.FloatWithValue(M, P, 1.0)
		if diagonal {
			B = matrix.FloatDiagonal(P, 1.0)
		} else {
			B = matrix.FloatWithValue(P, N, 1.0)
		}
	}
	C = matrix.FloatZeros(M, N)
	return
}

func runTest(A *matrix.FloatMatrix, ntest, LB int) time.Duration {

	var flags matops.Flags
	var mintime time.Duration

	N := A.Rows()
	ipiv := make([]int, N, N)
	flags = matops.LOWER
	if testUpper {
		flags = matops.UPPER
	}

	W := matrix.FloatZeros(A.Rows(), LB+2)
	fnc := func() {
		_, ERRmatops = matops.DecomposeBK(A, W, ipiv, flags, LB)
	}

	A0 := A.Copy()
	for n := 0; n < ntest; n++ {
		if n > 0 {
			// restore original A
			A0.CopyTo(A)
		}
		mperf.FlushCache()
		time0 := mperf.Timeit(fnc)
		if n == 0 || time0 < mintime {
			mintime = time0
		}
		if verbose {
			fmt.Printf("%.4f ms\n", time0.Seconds()*1000.0)
		}
	}
	return mintime
}

func runRefTest(A *matrix.FloatMatrix, ntest, LB int) time.Duration {

	var mintime time.Duration

	N := A.Cols()
	tau := matrix.FloatZeros(N, 1)

	fnc := func() {
		ERRlapack = lapack.Geqrf(A, tau)
	}

	A0 := A.Copy()
	for n := 0; n < ntest; n++ {
		if n > 0 {
			// restore original A
			A0.CopyTo(A)
			tau.Scale(0.0)
		}
		mperf.FlushCache()
		time0 := mperf.Timeit(fnc)
		if n == 0 || time0 < mintime {
			mintime = time0
		}
	}
	return mintime
}

func runRefTest(A *matrix.FloatMatrix, ntest, LB int) time.Duration {

	var flags matops.Flags
	var mintime time.Duration

	N := A.Rows()
	ipiv := make([]int, N, N)
	flags = matops.LOWER
	if testUpper {
		flags = matops.UPPER
	}

	W := matrix.FloatZeros(A.Rows(), LB+2)
	fnc := func() {
		_, ERRref = matops.DecomposeLDL(A, W, ipiv, flags, 0)
	}

	A0 := A.Copy()
	for n := 0; n < ntest; n++ {
		if n > 0 {
			// restore original A
			A0.CopyTo(A)
		}
		mperf.FlushCache()
		time0 := mperf.Timeit(fnc)
		if n == 0 || time0 < mintime {
			mintime = time0
		}
	}
	return mintime
}