Golang lapackd.NormP函数代码示例

func TestBidiagReduceUnblocked(t *testing.T) {

	N := 217
	M := 269
	conf := gomas.NewConf()
	conf.LB = 0

	zeromean := cmat.NewFloatNormSource()
	A := cmat.NewMatrix(M, N)
	A.SetFrom(zeromean)
	tauq := cmat.NewMatrix(N, 1)
	taup := cmat.NewMatrix(N, 1)

	At := cmat.NewMatrix(N, M)
	blasd.Transpose(At, A)
	tauqt := cmat.NewMatrix(N, 1)
	taupt := cmat.NewMatrix(N, 1)

	W := lapackd.Workspace(M + N)

	lapackd.BDReduce(A, tauq, taup, W, conf)
	lapackd.BDReduce(At, tauqt, taupt, W, conf)

	// BiRed(A) == BiRed(A.T).T
	blasd.Plus(At, A, 1.0, -1.0, gomas.TRANSB)
	blasd.Axpy(tauqt, taup, -1.0)
	blasd.Axpy(taupt, tauq, -1.0)

	nrm := lapackd.NormP(At, lapackd.NORM_ONE)
	t.Logf("M=%d, N=%d || BiRed(A) - BiRed(A.T).T||_1 : %e\n", M, N, nrm)
	nrm = lapackd.NormP(taupt, lapackd.NORM_ONE)
	t.Logf("  || BiRed(A).tauq - BiRed(A.T).taup||_1 : %e\n", nrm)
	nrm = lapackd.NormP(tauqt, lapackd.NORM_ONE)
	t.Logf("  || BiRed(A).taup - BiRed(A.T).tauq||_1 : %e\n", nrm)
}

func test_bdsvd(N, flags, kind int, verbose bool, t *testing.T) {
	var At, sD, sE, tmp cmat.FloatMatrix

	uplo := "upper"
	offdiag := 1
	if flags&gomas.LOWER != 0 {
		offdiag = -1
		uplo = "lower"
	}
	A0 := cmat.NewMatrix(N, N)
	desc := setDiagonals(A0, offdiag, kind)
	At.SubMatrix(A0, 0, 0, N, N)
	sD.Diag(A0, 0)
	sE.Diag(A0, offdiag)
	D := cmat.NewCopy(&sD)
	E := cmat.NewCopy(&sE)

	// unit singular vectors
	U := cmat.NewMatrix(N, N)
	sD.Diag(U, 0)
	sD.Add(1.0)

	V := cmat.NewMatrix(N, N)
	sD.Diag(V, 0)
	sD.Add(1.0)

	W := cmat.NewMatrix(4*N, 1)
	C := cmat.NewMatrix(N, N)

	lapackd.BDSvd(D, E, U, V, W, flags|gomas.WANTU|gomas.WANTV)

	blasd.Mult(C, U, U, 1.0, 0.0, gomas.TRANSA)
	sD.Diag(C)
	sD.Add(-1.0)
	nrmu := lapackd.NormP(C, lapackd.NORM_ONE)

	blasd.Mult(C, V, V, 1.0, 0.0, gomas.TRANSB)
	sD.Add(-1.0)
	nrmv := lapackd.NormP(C, lapackd.NORM_ONE)

	blasd.Mult(C, U, A0, 1.0, 0.0, gomas.TRANSA)
	blasd.Mult(&At, C, V, 1.0, 0.0, gomas.TRANSB)
	if verbose && N < 10 {
		t.Logf("D:\n%v\n", asRow(&tmp, D))
		t.Logf("U:\n%v\n", U)
		t.Logf("V:\n%v\n", V)
		t.Logf("U.T*A*V\n%v\n", &At)
	}
	sD.Diag(&At)
	blasd.Axpy(&sD, D, -1.0)
	nrma := lapackd.NormP(&At, lapackd.NORM_ONE)

	t.Logf("N=%d [%s,%s] ||U.T*A*V - bdsvd(A)||_1: %e\n", N, uplo, desc, nrma)
	t.Logf("  ||I - U.T*U||_1: %e\n", nrmu)
	t.Logf("  ||I - V.T*V||_1: %e\n", nrmv)
}

// test: C = C*Q.T
func TestQLMultRightTrans(t *testing.T) {
	var d, di0, di1 cmat.FloatMatrix
	M := 891
	N := 853
	lb := 36
	conf := gomas.NewConf()

	A := cmat.NewMatrix(M, N)
	src := cmat.NewFloatNormSource()
	A.SetFrom(src)

	C0 := cmat.NewMatrix(N, M)
	d.Diag(C0, M-N)
	ones := cmat.NewFloatConstSource(1.0)
	d.SetFrom(ones)
	C1 := cmat.NewCopy(C0)

	I0 := cmat.NewMatrix(N, N)
	I1 := cmat.NewCopy(I0)
	di0.Diag(I0)
	di1.Diag(I1)

	tau := cmat.NewMatrix(N, 1)
	W := cmat.NewMatrix(lb*(M+N), 1)

	conf.LB = lb
	lapackd.QLFactor(A, tau, W, conf)

	conf.LB = 0
	lapackd.QLMult(C0, A, tau, W, gomas.RIGHT|gomas.TRANS, conf)
	// I = Q*Q.T - I
	blasd.Mult(I0, C0, C0, 1.0, 0.0, gomas.TRANSB, conf)
	blasd.Add(&di0, -1.0)
	n0 := lapackd.NormP(I0, lapackd.NORM_ONE)

	conf.LB = lb
	lapackd.QLMult(C1, A, tau, W, gomas.RIGHT|gomas.TRANS, conf)
	// I = Q*Q.T - I
	blasd.Mult(I1, C1, C1, 1.0, 0.0, gomas.TRANSB, conf)
	blasd.Add(&di1, -1.0)
	n1 := lapackd.NormP(I1, lapackd.NORM_ONE)

	if N < 10 {
		t.Logf("unblk C0*Q:\n%v\n", C0)
		t.Logf("blk. C2*Q:\n%v\n", C1)
	}
	blasd.Plus(C0, C1, 1.0, -1.0, gomas.NONE)
	n2 := lapackd.NormP(C0, lapackd.NORM_ONE)

	t.Logf("M=%d, N=%d ||unblk.QLMult(C) - blk.QLMult(C)||_1: %e\n", M, N, n2)
	t.Logf("unblk M=%d, N=%d ||I - Q*Q.T||_1: %e\n", M, N, n0)
	t.Logf("blk   M=%d, N=%d ||I - Q*Q.T||_1: %e\n", M, N, n1)
}

func TestUpperCHOL(t *testing.T) {
	N := 311
	K := 43
	nb := 0

	conf := gomas.NewConf()
	conf.LB = nb

	Z := cmat.NewMatrix(N, N)
	A := cmat.NewMatrix(N, N)
	A0 := cmat.NewMatrix(N, N)
	B := cmat.NewMatrix(N, K)
	X := cmat.NewMatrix(N, K)

	unitrand := cmat.NewFloatUniformSource()
	Z.SetFrom(unitrand)

	blasd.Mult(A, Z, Z, 1.0, 0.0, gomas.TRANSB)
	A0.Copy(A)

	B.SetFrom(unitrand)
	X.Copy(B)

	// A = chol(A) = U.T*U
	t.Logf("Unblocked version: nb=%d\n", conf.LB)
	lapackd.CHOLFactor(A, gomas.UPPER, conf)
	// X = A.-1*B = U.-1*(U.-T*B)
	lapackd.CHOLSolve(X, A, gomas.UPPER)
	// B = B - A*X
	blasd.Mult(B, A0, X, -1.0, 1.0, gomas.NONE)
	// ||B - A*X||_1
	nrm := lapackd.NormP(B, lapackd.NORM_ONE)
	t.Logf("N=%d:  ||B - A*X||_1: %e\n", N, nrm)

	// A = chol(A) = U.T*U
	A.Copy(A0)
	B.SetFrom(unitrand)
	X.Copy(B)
	conf.LB = 16
	t.Logf("Blocked version: nb=%d\n", conf.LB)
	lapackd.CHOLFactor(A, gomas.UPPER, conf)
	// X = A.-1*B = U.-1*(U.-T*B)
	lapackd.CHOLSolve(X, A, gomas.UPPER)
	// B = B - A*X
	blasd.Mult(B, A0, X, -1.0, 1.0, gomas.NONE)
	// ||B - A*X||_1
	nrm = lapackd.NormP(B, lapackd.NORM_ONE)
	t.Logf("N=%d:  ||B - A*X||_1: %e\n", N, nrm)

}

func TestTrdMultUpper(t *testing.T) {
	var dt, et, da, ea cmat.FloatMatrix
	N := 843
	nb := 48
	conf := gomas.NewConf()
	conf.LB = nb

	A := cmat.NewMatrix(N, N)
	tau := cmat.NewMatrix(N, 1)
	src := cmat.NewFloatNormSource()
	// create symmetric matrix
	A.SetFrom(src, cmat.SYMM)
	A0 := cmat.NewCopy(A)

	W := lapackd.Workspace(lapackd.TRDReduceWork(A, conf))
	lapackd.TRDReduce(A, tau, W, gomas.UPPER, conf)

	// make tridiagonal matrix T
	T0 := cmat.NewMatrix(N, N)
	dt.Diag(T0)
	da.Diag(A)
	blasd.Copy(&dt, &da)

	ea.Diag(A, 1)
	et.Diag(T0, 1)
	blasd.Copy(&et, &ea)
	et.Diag(T0, -1)
	blasd.Copy(&et, &ea)
	T1 := cmat.NewCopy(T0)

	// compute Q*T*Q.T (unblocked)
	conf.LB = 0
	lapackd.TRDMult(T0, A, tau, W, gomas.LEFT|gomas.UPPER, conf)
	lapackd.TRDMult(T0, A, tau, W, gomas.RIGHT|gomas.TRANS|gomas.UPPER, conf)

	blasd.Plus(T0, A0, 1.0, -1.0, gomas.NONE)
	nrm := lapackd.NormP(T0, lapackd.NORM_ONE)
	t.Logf("N=%d, unblk.||A - Q*T*Q.T||_1: %e\n", N, nrm)

	// compute Q*T*Q.T (blocked)
	conf.LB = nb
	W = lapackd.Workspace(lapackd.TRDMultWork(A, gomas.LEFT|gomas.UPPER, conf))
	lapackd.TRDMult(T1, A, tau, W, gomas.LEFT|gomas.UPPER, conf)
	lapackd.TRDMult(T1, A, tau, W, gomas.RIGHT|gomas.TRANS|gomas.UPPER, conf)

	blasd.Plus(T1, A0, 1.0, -1.0, gomas.NONE)
	nrm = lapackd.NormP(T1, lapackd.NORM_ONE)
	t.Logf("N=%d,   blk.||A - Q*T*Q.T||_1: %e\n", N, nrm)
}

func TestQLBuildwithK(t *testing.T) {
	var dc cmat.FloatMatrix
	M := 711
	N := 707
	K := 691
	lb := 36
	conf := gomas.NewConf()

	A := cmat.NewMatrix(M, N)
	src := cmat.NewFloatNormSource()
	A.SetFrom(src)
	tau := cmat.NewMatrix(N, 1)

	W := cmat.NewMatrix(M+N, 1)
	C := cmat.NewMatrix(N, N)

	conf.LB = lb
	lapackd.QLFactor(A, tau, W, conf)
	A1 := cmat.NewCopy(A)

	conf.LB = 0
	lapackd.QLBuild(A, tau, W, K, conf)
	blasd.Mult(C, A, A, 1.0, 0.0, gomas.TRANSA, conf)
	dc.Diag(C)
	blasd.Add(&dc, -1.0)
	if N < 10 {
		t.Logf("unblk.QLBuild Q:\n%v\n", A)
		t.Logf("unblk.QLBuild Q.T*Q:\n%v\n", C)
	}
	n0 := lapackd.NormP(C, lapackd.NORM_ONE)

	conf.LB = lb
	W1 := lapackd.Workspace(lapackd.QLBuildWork(A1, conf))
	lapackd.QLBuild(A1, tau, W1, K, conf)
	if N < 10 {
		t.Logf("blk.QLBuild Q:\n%v\n", A1)
	}
	// compute: I - Q.T*Q
	blasd.Mult(C, A1, A1, 1.0, 0.0, gomas.TRANSA, conf)
	blasd.Add(&dc, -1.0)
	n1 := lapackd.NormP(C, lapackd.NORM_ONE)

	blasd.Plus(A, A1, 1.0, -1.0, gomas.NONE)
	n2 := lapackd.NormP(A, lapackd.NORM_ONE)

	t.Logf("M=%d, N=%d, K=%d ||unblk.QLBuild(A) - blk.QLBuild(A)||_1 :%e\n", M, N, K, n2)
	t.Logf("unblk M=%d, N=%d, K=%d ||Q.T*Q - I||_1 : %e\n", M, N, K, n0)
	t.Logf("blk   M=%d, N=%d, K=%d ||Q.T*Q - I||_1 : %e\n", M, N, K, n1)
}

func test_trdevd(N, flags, kind int, verbose bool, t *testing.T) {
	var At, sD, sE, tmp cmat.FloatMatrix

	A0 := cmat.NewMatrix(N, N)
	desc := setTrdDiagonals(A0, kind)
	At.SubMatrix(A0, 0, 0, N, N)
	sD.Diag(A0, 0)
	sE.Diag(A0, 1)
	D := cmat.NewCopy(&sD)
	E := cmat.NewCopy(&sE)

	V := cmat.NewMatrix(N, N)
	sD.Diag(V, 0)
	sD.Add(1.0)

	W := cmat.NewMatrix(4*N, 1)
	C := cmat.NewMatrix(N, N)

	if verbose && N < 10 {
		t.Logf("A0:\n%v\n", A0.ToString("%6.3f"))
		t.Logf("V.pre:\n%v\n", V.ToString("%6.3f"))
	}
	lapackd.TRDEigen(D, E, V, W, flags|gomas.WANTV)
	for k := 0; k < N-1; k++ {
		if E.GetAt(k) != 0.0 {
			t.Logf("E[%d] != 0.0 (%e)\n", k, E.GetAt(k))
		}
	}

	blasd.Mult(C, V, V, 1.0, 0.0, gomas.TRANSB)
	sD.Diag(C)
	sD.Add(-1.0)
	nrmv := lapackd.NormP(C, lapackd.NORM_ONE)

	blasd.Mult(C, V, A0, 1.0, 0.0, gomas.TRANSA)
	blasd.Mult(&At, C, V, 1.0, 0.0, gomas.NONE)
	if verbose && N < 10 {
		t.Logf("D:\n%v\n", asRow(&tmp, D).ToString("%6.3f"))
		t.Logf("V:\n%v\n", V.ToString("%6.3f"))
		t.Logf("V.T*A*V\n%v\n", At.ToString("%6.3f"))
	}
	sD.Diag(&At)
	blasd.Axpy(&sD, D, -1.0)
	nrma := lapackd.NormP(&At, lapackd.NORM_ONE)

	t.Logf("N=%d [%s] ||V.T*A*V - eigen(A)||_1: %e\n", N, desc, nrma)
	t.Logf("  ||I - V.T*V||_1: %e\n", nrmv)
}

func TestLQBuild(t *testing.T) {
	var dc cmat.FloatMatrix

	M := 877
	N := 913
	K := 831
	lb := 48
	conf := gomas.NewConf()
	_ = lb

	A := cmat.NewMatrix(M, N)
	src := cmat.NewFloatNormSource()
	A.SetFrom(src)
	tau := cmat.NewMatrix(M, 1)
	W := cmat.NewMatrix(M, 1)
	C := cmat.NewMatrix(M, M)
	dc.Diag(C)

	conf.LB = lb
	lapackd.LQFactor(A, tau, W, conf)
	A1 := cmat.NewCopy(A)

	conf.LB = 0
	lapackd.LQBuild(A, tau, W, K, conf)
	if N < 10 {
		t.Logf("unblk.LQBuild Q:\n%v\n", A)
	}
	blasd.Mult(C, A, A, 1.0, 0.0, gomas.TRANSB, conf)
	blasd.Add(&dc, -1.0)
	n0 := lapackd.NormP(C, lapackd.NORM_ONE)

	conf.LB = lb
	W2 := lapackd.Workspace(lapackd.LQBuildWork(A, conf))
	lapackd.LQBuild(A1, tau, W2, K, conf)
	if N < 10 {
		t.Logf("blk.LQBuild Q:\n%v\n", A1)
	}
	blasd.Mult(C, A1, A1, 1.0, 0.0, gomas.TRANSB, conf)
	blasd.Add(&dc, -1.0)
	n1 := lapackd.NormP(C, lapackd.NORM_ONE)

	blasd.Plus(A, A1, 1.0, -1.0, gomas.NONE)
	n2 := lapackd.NormP(A, lapackd.NORM_ONE)

	t.Logf("M=%d, N=%d, K=%d ||unblk.LQBuild(A) - blk.LQBuild(A)||_1 :%e\n", M, N, K, n2)
	t.Logf("unblk M=%d, N=%d, K=%d ||I - Q*Q.T||_1 : %e\n", M, N, K, n0)
	t.Logf("  blk M=%d, N=%d, K=%d ||I - Q*Q.T||_1 : %e\n", M, N, K, n1)
}

func TestSolveBKLowerBig(t *testing.T) {
	N := 427
	normsrc := cmat.NewFloatNormSource(5.0, 10.0)
	A := cmat.NewMatrix(N, N)
	A.SetFrom(normsrc, cmat.LOWER)

	X := cmat.NewMatrix(N, 2)
	X.SetFrom(normsrc)
	B := cmat.NewCopy(X)
	blasd.MultSym(B, A, X, 1.0, 0.0, gomas.LOWER|gomas.LEFT)

	ipiv := lapackd.NewPivots(N)

	conf := gomas.NewConf()
	conf.LB = 16
	W := lapackd.Workspace(lapackd.BKFactorWork(A, conf))
	lapackd.BKFactor(A, W, ipiv, gomas.LOWER, conf)

	lapackd.BKSolve(B, A, ipiv, gomas.LOWER, conf)
	ok := B.AllClose(X)
	t.Logf("N=%d unblk.BK(X) == A.-1*B : %v\n", N, ok)
	blasd.Plus(B, X, 1.0, -1.0, gomas.NONE)
	nrm := lapackd.NormP(B, lapackd.NORM_ONE)
	t.Logf("  ||X - A.-1*B||_1: %.4e\n", nrm)

}

// QR decompose A, then compute ||A - (R.T*Q.T).T||_1, should be small
func TestUnblkQRMultRight(t *testing.T) {
	M := 711
	N := 593
	A := cmat.NewMatrix(M, N)
	C := cmat.NewMatrix(N, M)
	tau := cmat.NewMatrix(N, 1)

	zeromean := cmat.NewFloatNormSource()
	A.SetFrom(zeromean)
	A0 := cmat.NewCopy(A)

	conf := gomas.NewConf()
	conf.LB = 0

	// QR = A = Q*R
	W := lapackd.Workspace(lapackd.QRFactorWork(A, conf))
	lapackd.QRFactor(A, tau, W, conf)

	// C = transpose(TriU(QR)) = R.T
	C.Transpose(cmat.TriU(cmat.NewCopy(A), cmat.NONE))

	// C = C*Q.T = R.T*Q.T
	W = lapackd.Workspace(lapackd.QRMultWork(C, gomas.RIGHT, conf))
	err := lapackd.QRMult(C, A, tau, W, gomas.RIGHT|gomas.TRANS, conf)
	if err != nil {
		t.Logf("err: %v\n", err)
	}

	// A = A - QR
	blasd.Plus(A0, C, 1.0, -1.0, gomas.TRANSB)
	// ||A - Q*R||_1
	nrm := lapackd.NormP(A0, lapackd.NORM_ONE)
	t.Logf("M=%d,N=%d  ||A - (R.T*Q.T).T||_1: %e\n", M, N, nrm)
}

// test: unblk.ReduceHess(A) == blk.ReduceHess(A)
func TestReduceHess(t *testing.T) {
	N := 375
	nb := 16

	conf := gomas.NewConf()
	conf.LB = nb

	A := cmat.NewMatrix(N, N)
	tau := cmat.NewMatrix(N, 1)
	zeromean := cmat.NewFloatNormSource()
	A.SetFrom(zeromean)

	A0 := cmat.NewCopy(A)
	tau0 := cmat.NewCopy(tau)

	// blocked reduction
	W := lapackd.Workspace(lapackd.HessReduceWork(A, conf))
	lapackd.HessReduce(A, tau, W, conf)

	// unblocked reduction
	conf.LB = 0
	lapackd.HessReduce(A0, tau0, W, conf)

	ok := A.AllClose(A0)
	t.Logf("blk.ReduceHess(A) == unblk.ReduceHess(A): %v\n", ok)

	ok = tau0.AllClose(tau)
	t.Logf("blk HessQ.tau == unblk HessQ.tau: %v\n", ok)

	// ||A - A0||_1
	blasd.Plus(A, A0, 1.0, -1.0, gomas.NONE)
	nrm := lapackd.NormP(A, lapackd.NORM_ONE)
	t.Logf("||H - H0||_1: %e\n", nrm)
}

// test: min ||X|| s.t A.T*X = B
func TestSolveQR(t *testing.T) {
	M := 799
	N := 711
	K := 241
	nb := 32
	conf := gomas.NewConf()
	conf.LB = nb

	tau := cmat.NewMatrix(N, 1)
	A := cmat.NewMatrix(M, N)
	src := cmat.NewFloatNormSource()
	A.SetFrom(src)
	A0 := cmat.NewCopy(A)
	B0 := cmat.NewMatrix(M, K)
	B0.SetFrom(src)
	B := cmat.NewCopy(B0)

	W := lapackd.Workspace(lapackd.QRFactorWork(A, conf))
	lapackd.QRFactor(A, tau, W, conf)

	lapackd.QRSolve(B, A, tau, W, gomas.TRANS, conf)

	var Bmin cmat.FloatMatrix
	Bmin.SubMatrix(B0, 0, 0, N, K)
	blasd.Mult(&Bmin, A0, B, 1.0, -1.0, gomas.TRANSA, conf)

	nrm := lapackd.NormP(&Bmin, lapackd.NORM_ONE)
	t.Logf("M=%d, N=%d ||B - A.T*X||_1: %e\n", M, N, nrm)
}

// test: min || B - A.T*X ||
func TestLeastSquaresLQ(t *testing.T) {
	M := 723
	N := 811
	K := 273
	nb := 32
	conf := gomas.NewConf()
	conf.LB = nb

	tau := cmat.NewMatrix(M, 1)
	A := cmat.NewMatrix(M, N)
	src := cmat.NewFloatNormSource()
	A.SetFrom(src)
	B0 := cmat.NewMatrix(M, K)
	B0.SetFrom(src)
	B := cmat.NewMatrix(N, K)

	// B = A.T*B0
	blasd.Mult(B, A, B0, 1.0, 0.0, gomas.TRANSA, conf)

	W := lapackd.Workspace(lapackd.LQFactorWork(A, conf))
	lapackd.LQFactor(A, tau, W, conf)

	// B' = A.-1*B
	lapackd.LQSolve(B, A, tau, W, gomas.TRANS, conf)

	// expect B[0:M,0:K] == B0[0:M,0:K], B[M:N,0:K] == 0
	var X cmat.FloatMatrix

	X.SubMatrix(B, 0, 0, M, K)
	blasd.Plus(&X, B0, 1.0, -1.0, gomas.NONE)
	nrm := lapackd.NormP(&X, lapackd.NORM_ONE)

	t.Logf("M=%d, N=%d  ||B0 - min( ||A.T*X - B0|| ) ||_1: %e\n", M, N, nrm)
}

// m > n: A[m,n], I[m,m] --> A == I*A == Q*Q.T*A
func TestQRTMultLeftIdent(t *testing.T) {
	M := 411
	N := 399
	nb := 16
	A := cmat.NewMatrix(M, N)
	T := cmat.NewMatrix(nb, N)

	zeromean := cmat.NewFloatNormSource()
	A.SetFrom(zeromean)
	A0 := cmat.NewCopy(A)
	C := cmat.NewCopy(A)
	conf := gomas.NewConf()
	conf.LB = nb
	//t.Logf("A0:\n%v\n", A0)

	// QR = A = Q*R
	W := lapackd.Workspace(lapackd.QRTFactorWork(A, conf))
	lapackd.QRTFactor(A, T, W, conf)
	//t.Logf("T:\n%v\n", T)

	// C = Q.T*A
	W = lapackd.Workspace(lapackd.QRTMultWork(C, T, gomas.LEFT, conf))
	lapackd.QRTMult(C, A, T, W, gomas.LEFT|gomas.TRANS, conf)

	// C = Q*C == Q*Q.T*A
	lapackd.QRTMult(C, A, T, W, gomas.LEFT, conf)
	//t.Logf("A*Q*Q.T:\n%v\n", C)

	// A = A - Q*Q.T*A
	blasd.Plus(A0, C, 1.0, -1.0, gomas.NONE)
	// ||A - Q*Q.T*A||_1
	nrm := lapackd.NormP(A0, lapackd.NORM_ONE)
	t.Logf("M=%d,N=%d  ||A - Q*Q.T*A||_1: %e\n", M, N, nrm)
}

func TestTriRedUpper(t *testing.T) {
	N := 843
	nb := 48
	conf := gomas.NewConf()
	conf.LB = 0

	A := cmat.NewMatrix(N, N)
	tau := cmat.NewMatrix(N, 1)
	src := cmat.NewFloatNormSource()
	A.SetFrom(src, cmat.UPPER)
	A1 := cmat.NewCopy(A)
	tau1 := cmat.NewCopy(tau)
	_ = A1

	W := lapackd.Workspace(N)
	W1 := lapackd.Workspace(N * nb)

	lapackd.TRDReduce(A, tau, W, gomas.UPPER, conf)

	conf.LB = nb
	lapackd.TRDReduce(A1, tau1, W1, gomas.UPPER, conf)

	blasd.Plus(A, A1, -1.0, 1.0, gomas.NONE)
	nrm := lapackd.NormP(A, lapackd.NORM_ONE)
	t.Logf("N=%d, ||unblk.Trired(A) - blk.Trired(A)||_1: %e\n", N, nrm)
	blasd.Axpy(tau, tau1, -1.0)
	nrm = blasd.Nrm2(tau)
	t.Logf("   ||unblk.Trired(tau) - blk.Trired(tau)||_1: %e\n", nrm)
}