Golang FloatMatrix.SetFrom方法代碼示例

/*
 * Solve a system of linear equations A*X = B with general M-by-N
 * matrix A using the QR factorization computed by DecomposeQRT().
 *
 * If flags&gomas.TRANS != 0:
 *   find the minimum norm solution of an overdetermined system A.T * X = B.
 *   i.e min ||X|| s.t A.T*X = B
 *
 * Otherwise:
 *   find the least squares solution of an overdetermined system, i.e.,
 *   solve the least squares problem: min || B - A*X ||.
 *
 * Arguments:
 *  B     On entry, the right hand side N-by-P matrix B. On exit, the solution matrix X.
 *
 *  A     The elements on and above the diagonal contain the min(M,N)-by-N upper
 *        trapezoidal matrix R. The elements below the diagonal with the matrix 'T',
 *        represent the ortogonal matrix Q as product of elementary reflectors.
 *        Matrix A and T are as returned by DecomposeQRT()
 *
 *  T     The block reflector computed from elementary reflectors as returned by
 *        DecomposeQRT() or computed from elementary reflectors and scalar coefficients
 *        by BuildT()
 *
 *  W     Workspace, size as returned by WorkspaceMultQT()
 *
 *  flags Indicator flag
 *
 *  conf  Blocking configuration
 *
 * Compatible with lapack.GELS (the m >= n part)
 */
func QRTSolve(B, A, T, W *cmat.FloatMatrix, flags int, confs ...*gomas.Config) *gomas.Error {
	var err *gomas.Error = nil
	var R, BT cmat.FloatMatrix
	conf := gomas.CurrentConf(confs...)

	if flags&gomas.TRANS != 0 {
		// Solve overdetermined system A.T*X = B

		// B' = R.-1*B
		R.SubMatrix(A, 0, 0, n(A), n(A))
		BT.SubMatrix(B, 0, 0, n(A), n(B))
		err = blasd.SolveTrm(&BT, &R, 1.0, gomas.LEFT|gomas.UPPER|gomas.TRANSA, conf)

		// Clear bottom part of B
		BT.SubMatrix(B, n(A), 0)
		BT.SetFrom(cmat.NewFloatConstSource(0.0))

		// X = Q*B'
		err = QRTMult(B, A, T, W, gomas.LEFT, conf)
	} else {
		// solve least square problem min ||A*X - B||

		// B' = Q.T*B
		err = QRTMult(B, A, T, W, gomas.LEFT|gomas.TRANS, conf)
		if err != nil {
			return err
		}

		// X = R.-1*B'
		R.SubMatrix(A, 0, 0, n(A), n(A))
		BT.SubMatrix(B, 0, 0, n(A), n(B))
		err = blasd.SolveTrm(&BT, &R, 1.0, gomas.LEFT|gomas.UPPER, conf)
	}
	return err
}

func TestDTrmmUnitUpperRight(t *testing.T) {
	var d cmat.FloatMatrix
	N := 563
	K := 171

	A := cmat.NewMatrix(N, N)
	B := cmat.NewMatrix(K, N)
	B0 := cmat.NewMatrix(K, N)
	C := cmat.NewMatrix(K, N)

	zeros := cmat.NewFloatConstSource(0.0)
	ones := cmat.NewFloatConstSource(1.0)
	zeromean := cmat.NewFloatNormSource()

	A.SetFrom(zeromean, cmat.UPPER|cmat.UNIT)
	B.SetFrom(ones)
	B0.SetFrom(ones)
	// B = B*A
	blasd.MultTrm(B, A, 1.0, gomas.UPPER|gomas.RIGHT|gomas.UNIT)
	d.Diag(A).SetFrom(ones)
	blasd.Mult(C, B0, A, 1.0, 0.0, gomas.NONE)
	ok := C.AllClose(B)
	t.Logf("trmm(B, A, R|U|N|U) == gemm(C, B, TriUU(A))   : %v\n", ok)

	B.SetFrom(ones)
	// B = B*A.T
	d.SetFrom(zeros)
	blasd.MultTrm(B, A, 1.0, gomas.UPPER|gomas.RIGHT|gomas.TRANSA|gomas.UNIT)
	d.SetFrom(ones)
	blasd.Mult(C, B0, A, 1.0, 0.0, gomas.TRANSB)
	ok = C.AllClose(B)
	t.Logf("trmm(B, A, R|U|T|U) == gemm(C, B, TriUU(A).T) : %v\n", ok)
}

func TestPartition2D(t *testing.T) {
	var ATL, ATR, ABL, ABR, As cmat.FloatMatrix
	var A00, a01, A02, a10, a11, a12, A20, a21, A22 cmat.FloatMatrix

	csource := cmat.NewFloatConstSource(1.0)
	A := cmat.NewMatrix(6, 6)
	As.SubMatrix(A, 1, 1, 4, 4)
	As.SetFrom(csource)

	Partition2x2(&ATL, &ATR, &ABL, &ABR, &As, 0, 0, PTOPLEFT)
	t.Logf("ATL:\n%v\n", &ATL)

	t.Logf("n(ATL)=%d, n(As)=%d\n", n(&ATL), n(&As))
	k := 0
	for n(&ATL) < n(&As) && k < n(&As) {
		Repartition2x2to3x3(&ATL,
			&A00, &a01, &A02,
			&a10, &a11, &a12,
			&A20, &a21, &A22, &As, 1, PBOTTOMRIGHT)
		t.Logf("n(A00)=%d, n(a01)=%d, n(A02)=%d\n", n(&A00), n(&a01), n(&A02))
		t.Logf("n(a10)=%d, n(a11)=%d, n(a12)=%d\n", n(&a10), n(&a11), n(&a12))
		t.Logf("n(A20)=%d, n(a21)=%d, n(A22)=%d\n", n(&A20), n(&a21), n(&A22))
		//t.Logf("n(a12)=%d [%d], n(a11)=%d\n", n(&a12), a12.Len(), a11.Len())
		a11.Set(0, 0, a11.Get(0, 0)+1.0)
		addConst(&a21, -2.0)

		Continue3x3to2x2(&ATL, &ATR, &ABL, &ABR, &A00, &a11, &A22, &As, PBOTTOMRIGHT)
		t.Logf("n(ATL)=%d, n(As)=%d\n", n(&ATL), n(&As))
		k += 1
	}
	t.Logf("A:\n%v\n", A)
}

// Simple and slow LQ decomposition with Givens rotations
func TestGivensLQ(t *testing.T) {
	var d cmat.FloatMatrix
	M := 149
	N := 167
	A := cmat.NewMatrix(M, N)
	A1 := cmat.NewCopy(A)

	ones := cmat.NewFloatConstSource(1.0)
	src := cmat.NewFloatNormSource()
	A.SetFrom(src)
	A0 := cmat.NewCopy(A)

	Qt := cmat.NewMatrix(N, N)
	d.Diag(Qt)
	d.SetFrom(ones)

	// R = G(n)...G(2)G(1)*A; Q = G(1).T*G(2).T...G(n).T ;  Q.T = G(n)...G(2)G(1)
	for i := 0; i < M; i++ {
		// zero elements right of diagonal
		for j := N - 2; j >= i; j-- {
			c, s, r := lapackd.ComputeGivens(A.Get(i, j), A.Get(i, j+1))
			A.Set(i, j, r)
			A.Set(i, j+1, 0.0)
			// apply rotation to this column starting from row i+1
			lapackd.ApplyGivensRight(A, j, j+1, i+1, M-i-1, c, s)
			// update Qt = G(k)*Qt
			lapackd.ApplyGivensRight(Qt, j, j+1, 0, N, c, s)
		}
	}
	// A = L*Q
	blasd.Mult(A1, A, Qt, 1.0, 0.0, gomas.TRANSB)
	blasd.Plus(A0, A1, 1.0, -1.0, gomas.NONE)
	nrm := lapackd.NormP(A0, lapackd.NORM_ONE)
	t.Logf("M=%d, N=%d ||A - L*G(1)..G(n)||_1: %e\n", M, N, nrm)
}

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

// test: A - Q*Hess(A)*Q.T  == 0
func TestMultHess(t *testing.T) {
	N := 377
	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)

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

	var Hlow cmat.FloatMatrix
	H := cmat.NewCopy(A)

	// set triangular part below first subdiagonal to zeros
	zeros := cmat.NewFloatConstSource(0.0)
	Hlow.SubMatrix(H, 1, 0, N-1, N-1)
	Hlow.SetFrom(zeros, cmat.LOWER|cmat.UNIT)
	H1 := cmat.NewCopy(H)

	// H := Q*H*Q.T
	conf.LB = nb
	lapackd.HessMult(H, A, tau, W, gomas.LEFT, conf)
	lapackd.HessMult(H, A, tau, W, gomas.RIGHT|gomas.TRANS, conf)

	// H := Q*H*Q.T
	conf.LB = 0
	lapackd.HessMult(H1, A, tau, W, gomas.LEFT, conf)
	lapackd.HessMult(H1, A, tau, W, gomas.RIGHT|gomas.TRANS, conf)

	// compute ||Q*Hess(A)*Q.T - A||_1
	blasd.Plus(H, A0, 1.0, -1.0, gomas.NONE)
	nrm := lapackd.NormP(H, lapackd.NORM_ONE)
	t.Logf("  blk.|| Q*Hess(A)*Q.T - A ||_1 : %e\n", nrm)

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

func testEigen(N int, bits int, t *testing.T) {
	var A, A0, W, D, V *cmat.FloatMatrix
	var sD cmat.FloatMatrix
	var s string = "lower"

	if bits&gomas.UPPER != 0 {
		s = "upper"
	}

	wsize := N * N
	if wsize < 100 {
		wsize = 100
	}

	D = cmat.NewMatrix(N, 1)
	A = cmat.NewMatrix(N, N)
	V = cmat.NewMatrix(N, N)

	src := cmat.NewFloatNormSource()
	A.SetFrom(src, cmat.SYMM)
	A0 = cmat.NewCopy(A)
	W = cmat.NewMatrix(wsize, 1)

	if err := lapackd.EigenSym(D, A, W, bits|gomas.WANTV); err != nil {
		t.Errorf("EigenSym error: %v\n", err)
		return
	}

	// ||I - V.T*V||
	sD.Diag(V)
	blasd.Mult(V, A, A, 1.0, 0.0, gomas.TRANSA)
	blasd.Add(&sD, -1.0)
	nrm1 := lapackd.NormP(V, lapackd.NORM_ONE)

	// left vectors are M-by-N
	V.Copy(A)
	lapackd.MultDiag(V, D, gomas.RIGHT)
	blasd.Mult(A0, V, A, -1.0, 1.0, gomas.TRANSB)
	nrm2 := lapackd.NormP(A0, lapackd.NORM_ONE)

	t.Logf("N=%d, [%s] ||A - V*D*V.T||_1 :%e\n", N, s, nrm2)
	t.Logf("  ||I - V.T*V||_1 : %e\n", nrm1)
}

/*
 * Solve a system of linear equations A.T*X = B with general M-by-N
 * matrix A using the QR factorization computed by LQFactor().
 *
 * If flags&TRANS != 0:
 *   find the minimum norm solution of an overdetermined system A.T * X = B.
 *   i.e min ||X|| s.t A.T*X = B
 *
 * Otherwise:
 *   find the least squares solution of an overdetermined system, i.e.,
 *   solve the least squares problem: min || B - A*X ||.
 *
 * Arguments:
 *  B     On entry, the right hand side N-by-P matrix B. On exit, the solution matrix X.
 *
 *  A     The elements on and below the diagonal contain the M-by-min(M,N) lower
 *        trapezoidal matrix L. The elements right of the diagonal with the vector 'tau',
 *        represent the ortogonal matrix Q as product of elementary reflectors.
 *        Matrix A is as returned by LQFactor()
 *
 *  tau   The vector of N scalar coefficients that together with trilu(A) define
 *        the ortogonal matrix Q as Q = H(N)H(N-1)...H(1)
 *
 *  W     Workspace, size required returned WorksizeMultLQ().
 *
 *  flags Indicator flags
 *
 *  conf  Optinal blocking configuration. If not given default will be used. Unblocked
 *        invocation is indicated with conf.LB == 0.
 *
 * Compatible with lapack.GELS (the m < n part)
 */
func LQSolve(B, A, tau, W *cmat.FloatMatrix, flags int, confs ...*gomas.Config) *gomas.Error {
	var err *gomas.Error = nil
	var L, BL cmat.FloatMatrix

	conf := gomas.CurrentConf(confs...)

	wsmin := wsMultLQLeft(B, 0)
	if W.Len() < wsmin {
		return gomas.NewError(gomas.EWORK, "SolveLQ", wsmin)
	}

	if flags&gomas.TRANS != 0 {
		// solve: MIN ||A.T*X - B||

		// B' = Q.T*B
		err = LQMult(B, A, tau, W, gomas.LEFT, conf)
		if err != nil {
			return err
		}

		// X = L.-1*B'
		L.SubMatrix(A, 0, 0, m(A), m(A))
		BL.SubMatrix(B, 0, 0, m(A), n(B))
		err = blasd.SolveTrm(&BL, &L, 1.0, gomas.LEFT|gomas.LOWER|gomas.TRANSA, conf)

	} else {
		// Solve underdetermined system A*X = B

		// B' = L.-1*B
		L.SubMatrix(A, 0, 0, m(A), m(A))
		BL.SubMatrix(B, 0, 0, m(A), n(B))
		err = blasd.SolveTrm(&BL, &L, 1.0, gomas.LEFT|gomas.LOWER, conf)

		// Clear bottom part of B
		BL.SubMatrix(B, m(A), 0)
		BL.SetFrom(cmat.NewFloatConstSource(0.0))

		// X = Q.T*B'
		err = LQMult(B, A, tau, W, gomas.LEFT|gomas.TRANS, conf)

	}
	return err
}

// Simple and slow QR decomposition with Givens rotations
func TestGivensQR(t *testing.T) {
	var d cmat.FloatMatrix
	M := 181
	N := 159
	A := cmat.NewMatrix(M, N)
	A1 := cmat.NewCopy(A)

	ones := cmat.NewFloatConstSource(1.0)
	src := cmat.NewFloatNormSource()
	A.SetFrom(src)
	A0 := cmat.NewCopy(A)

	Qt := cmat.NewMatrix(M, M)
	d.Diag(Qt)
	d.SetFrom(ones)

	// R = G(n)...G(2)G(1)*A; Q = G(1).T*G(2).T...G(n).T ;  Q.T = G(n)...G(2)G(1)

	// for all columns ...
	for j := 0; j < N; j++ {
		// ... zero elements below diagonal, starting from bottom
		for i := M - 2; i >= j; i-- {
			c, s, r := lapackd.ComputeGivens(A.Get(i, j), A.Get(i+1, j))
			A.Set(i, j, r)
			A.Set(i+1, j, 0.0)
			// apply rotations on this row starting from column j, N-j column
			lapackd.ApplyGivensLeft(A, i, i+1, j+1, N-j-1, c, s)
			// update Qt = G(k)*Qt
			lapackd.ApplyGivensLeft(Qt, i, i+1, 0, M, c, s)
		}
	}
	// check: A = Q*R
	blasd.Mult(A1, Qt, A, 1.0, 0.0, gomas.TRANSA)
	blasd.Plus(A0, A1, 1.0, -1.0, gomas.NONE)
	nrm := lapackd.NormP(A0, lapackd.NORM_ONE)
	t.Logf("M=%d, N=%d ||A - G(n)..G(1)*R||_1: %e\n", M, N, nrm)
}

// test: M < N, U=[m,m] and V=[m,n] or V=[n,n] (square)
func testWide(M, N int, square bool, t *testing.T) {
	var A, A0, W, S, U, Uu, V, Vv *cmat.FloatMatrix
	var sD cmat.FloatMatrix
	var s string

	wsize := M * N
	if wsize < 100 {
		wsize = 100
	}

	S = cmat.NewMatrix(M, 1)
	A = cmat.NewMatrix(M, N)
	U = cmat.NewMatrix(M, M)
	Uu = cmat.NewMatrix(M, M)
	if square {
		V = cmat.NewMatrix(N, N)
		Vv = cmat.NewMatrix(N, N)
	} else {
		V = cmat.NewMatrix(M, N)
		Vv = cmat.NewMatrix(M, M)
	}

	src := cmat.NewFloatNormSource()
	A.SetFrom(src)
	A0 = cmat.NewCopy(A)
	W = cmat.NewMatrix(wsize, 1)

	if err := lapackd.SVD(S, U, V, A, W, gomas.WANTU|gomas.WANTV); err != nil {
		t.Errorf("SVD error: %v\n", err)
		return
	}

	// ||I - U.T*U||
	sD.Diag(Uu)
	blasd.Mult(Uu, U, U, 1.0, 0.0, gomas.TRANSA)
	blasd.Add(&sD, -1.0)
	nrm0 := lapackd.NormP(Uu, lapackd.NORM_ONE)

	// ||I - V*V.T||
	sD.Diag(Vv)
	blasd.Mult(Vv, V, V, 1.0, 0.0, gomas.TRANSB)
	blasd.Add(&sD, -1.0)
	nrm1 := lapackd.NormP(Vv, lapackd.NORM_ONE)

	if square {
		// right vectors are N-by-N
		Sg := cmat.NewMatrix(M, N)
		A1 := cmat.NewMatrix(M, N)
		sD.Diag(Sg)
		blasd.Copy(&sD, S)
		blasd.Mult(A1, Sg, V, 1.0, 0.0, gomas.NONE)
		blasd.Mult(A0, U, A1, -1.0, 1.0, gomas.NONE)
		s = "U=[m,m], V=[n,n]"
	} else {
		// right vectors are M-by-N
		lapackd.MultDiag(V, S, gomas.LEFT)
		blasd.Mult(A0, U, V, -1.0, 1.0, gomas.NONE)
		s = "U=[m,m], V=[m,n]"
	}
	nrm2 := lapackd.NormP(A0, lapackd.NORM_ONE)

	if N < 10 {
		t.Logf("A - U*S*V.T:\n%v\n", A0)
	}
	t.Logf("M=%d, N=%d, %s ||A - U*S*V.T||_1 :%e\n", M, N, s, nrm2)
	t.Logf("  ||I - U.T*U||_1 : %e\n", nrm0)
	t.Logf("  ||I - V*V.T||_1 : %e\n", nrm1)
}

func svdWide(S, U, V, A, W *cmat.FloatMatrix, bits int, conf *gomas.Config) (err *gomas.Error) {
	var uu, vv *cmat.FloatMatrix
	var tauq, taup, Wred, sD, sE, L, Vm cmat.FloatMatrix

	if (bits & (gomas.WANTU | gomas.WANTV)) != 0 {
		if W.Len() < 4*n(A) {
			err = gomas.NewError(gomas.ESIZE, "SVD")
			return
		}
	}
	tauq.SetBuf(m(A)-1, 1, m(A)-1, W.Data())
	taup.SetBuf(m(A), 1, m(A), W.Data()[tauq.Len():])
	wrl := W.Len() - 2*m(A) - 1
	Wred.SetBuf(wrl, 1, wrl, W.Data()[2*m(A)-1:])

	if svdCrossover(n(A), m(A)) {
		goto do_n_much_bigger
	}

	// reduce to bidiagonal form
	if err = BDReduce(A, &tauq, &taup, &Wred, conf); err != nil {
		return
	}

	sD.Diag(A)
	sE.Diag(A, -1)
	blasd.Copy(S, &sD)

	// leftt vectors
	if bits&gomas.WANTU != 0 {
		L.SubMatrix(A, 0, 0, m(A), m(A))
		U.Copy(&L)
		cmat.TriL(U, 0)
		if err = BDBuild(U, &tauq, &Wred, m(U), gomas.WANTQ|gomas.LOWER, conf); err != nil {
			return
		}
		uu = U
	}
	// right vectors
	if bits&gomas.WANTV != 0 {
		if m(V) == m(A) {
			// V is M-by-N; copy and make upper triangular
			V.Copy(A)
			//cmat.TriU(V, 0)
			if err = BDBuild(V, &taup, &Wred, m(V), gomas.WANTP, conf); err != nil {
				return
			}
		} else {
			// V is N-by-N
			eye := cmat.FloatDiagonalSource{1.0}
			V.SetFrom(&eye, cmat.SYMM)
			err = BDMult(V, A, &taup, &Wred, gomas.MULTP|gomas.LEFT|gomas.TRANS, conf)
			if err != nil {
				return
			}
		}
		vv = V
	}
	err = BDSvd(S, &sE, uu, vv, W, bits|gomas.LOWER)
	return

do_n_much_bigger:
	// here N >> M, use LQ factor first
	if err = LQFactor(A, &taup, &Wred, conf); err != nil {
		return
	}
	if bits&gomas.WANTV != 0 {
		if m(V) == m(A) {
			V.Copy(A)
			if err = LQBuild(V, &taup, &Wred, m(A), conf); err != nil {
				return
			}
		} else {
			// V is N-by-N
			eye := cmat.FloatDiagonalSource{1.0}
			V.SetFrom(&eye, cmat.SYMM)
			if err = LQMult(V, A, &taup, &Wred, gomas.RIGHT, conf); err != nil {
				return
			}
		}
	}
	L.SubMatrix(A, 0, 0, m(A), m(A))
	cmat.TriL(&L, 0)

	// resize tauq/taup for UPPER bidiagonal reduction
	tauq.SetBuf(m(A), 1, m(A), W.Data())
	taup.SetBuf(m(A)-1, 1, m(A)-1, W.Data()[tauq.Len():])

	// bidiagonal reduce
	if err = BDReduce(&L, &tauq, &taup, &Wred, conf); err != nil {
		return
	}

	if bits&gomas.WANTV != 0 {
		Vm.SubMatrix(V, 0, 0, m(A), n(A))
		err = BDMult(&Vm, &L, &taup, &Wred, gomas.MULTP|gomas.LEFT|gomas.TRANS, conf)
		if err != nil {
			return
		}
		vv = V
//.........這裏部分代碼省略.........

func svdSmall(S, U, V, A, W *cmat.FloatMatrix, bits int, conf *gomas.Config) (err *gomas.Error) {
	var r cmat.FloatMatrix
	var d0, d1, e0 float64

	err = nil
	K := m(A)
	if n(A) < K {
		K = n(A)
	}
	tau := cmat.NewMatrix(K, 1)
	if m(A) >= n(A) {
		if err = QRFactor(A, tau, W, conf); err != nil {
			return
		}
	} else {
		if err = LQFactor(A, tau, W, conf); err != nil {
			return
		}
	}
	if m(A) == 1 || n(A) == 1 {
		// either tall M-by-1 or wide 1-by-N
		S.SetAt(0, math.Abs(A.Get(0, 0)))

		if bits&gomas.WANTU != 0 {
			if n(A) == 1 {
				if n(U) == n(A) {
					// U is M-by-1
					U.Copy(A)
					QRBuild(U, tau, W, n(A), conf)
				} else {
					// U is M-by-M
					eye := cmat.FloatDiagonalSource{1.0}
					U.SetFrom(&eye, cmat.SYMM)
					if err = QRMult(U, A, tau, W, gomas.RIGHT, conf); err != nil {
						return
					}
				}
			} else {
				U.Set(0, 0, -1.0)
			}
		}
		if bits&gomas.WANTV != 0 {
			if m(A) == 1 {
				if m(V) == m(A) {
					// V is 1-by-N
					V.Copy(A)
					LQBuild(V, tau, W, m(A), conf)
				} else {
					// V is N-by-N
					eye := cmat.FloatDiagonalSource{1.0}
					V.SetFrom(&eye, cmat.SYMM)
					if err = LQMult(V, A, tau, W, gomas.RIGHT, conf); err != nil {
						return
					}
				}
			} else {
				// V is 1-by-1
				V.Set(0, 0, -1.0)
			}
		}
		return
	}

	// Use bdSvd2x2 functions
	d0 = A.Get(0, 0)
	d1 = A.Get(1, 1)
	if m(A) >= n(A) {
		e0 = A.Get(0, 1)
	} else {
		e0 = A.Get(1, 0)
	}

	if bits&(gomas.WANTU|gomas.WANTV) == 0 {
		// no vectors
		smin, smax := bdSvd2x2(d0, e0, d1)
		S.SetAt(0, math.Abs(smax))
		S.SetAt(1, math.Abs(smin))
		return
	}

	// at least left or right eigenvector wanted
	smin, smax, cosl, sinl, cosr, sinr := bdSvd2x2Vec(d0, e0, d1)
	if bits&gomas.WANTU != 0 {
		// left eigenvectors
		if m(A) >= n(A) {
			if n(U) == n(A) {
				U.Copy(A)
				if err = QRBuild(U, tau, W, n(A), conf); err != nil {
					return
				}
			} else {
				// U is M-by-M
				eye := cmat.FloatDiagonalSource{1.0}
				U.SetFrom(&eye, cmat.SYMM)
				if err = QRMult(U, A, tau, W, gomas.RIGHT, conf); err != nil {
					return
				}
			}
			ApplyGivensRight(U, 0, 1, 0, m(A), cosl, sinl)
		} else {
//.........這裏部分代碼省略.........

// Compute SVD when m(A) >= n(A)
func svdTall(S, U, V, A, W *cmat.FloatMatrix, bits int, conf *gomas.Config) (err *gomas.Error) {
	var uu, vv *cmat.FloatMatrix
	var tauq, taup, Wred, sD, sE, R, Un cmat.FloatMatrix

	if (bits & (gomas.WANTU | gomas.WANTV)) != 0 {
		if W.Len() < 4*n(A) {
			err = gomas.NewError(gomas.ESIZE, "SVD")
			return
		}
	}
	tauq.SetBuf(n(A), 1, n(A), W.Data())
	taup.SetBuf(n(A)-1, 1, n(A)-1, W.Data()[tauq.Len():])
	wrl := W.Len() - 2*n(A) - 1
	Wred.SetBuf(wrl, 1, wrl, W.Data()[2*n(A)-1:])

	if svdCrossover(m(A), n(A)) {
		goto do_m_much_bigger
	}

	// reduce to bidiagonal form
	if err = BDReduce(A, &tauq, &taup, &Wred, conf); err != nil {
		return
	}

	sD.Diag(A)
	sE.Diag(A, 1)
	blasd.Copy(S, &sD)

	// left vectors
	if bits&gomas.WANTU != 0 {
		if n(U) == n(A) {
			// U is M-by-N; copy and make lower triangular
			U.Copy(A)
			cmat.TriL(U, 0)
			if err = BDBuild(U, &tauq, &Wred, n(U), gomas.WANTQ, conf); err != nil {
				return
			}
		} else {
			// U is M-by-M
			eye := cmat.FloatDiagonalSource{1.0}
			U.SetFrom(&eye, cmat.SYMM)
			if err = BDMult(U, A, &tauq, &Wred, gomas.MULTQ|gomas.RIGHT, conf); err != nil {
				return
			}
		}
		uu = U
	}
	// right vectors
	if bits&gomas.WANTV != 0 {
		R.SubMatrix(A, 0, 0, n(A), n(A))
		V.Copy(&R)
		cmat.TriU(V, 0)
		if err = BDBuild(V, &taup, &Wred, m(V), gomas.WANTP, conf); err != nil {
			return
		}
		vv = V
	}
	err = BDSvd(S, &sE, uu, vv, W, bits|gomas.UPPER)
	return

do_m_much_bigger:
	// M >> N here; first use QR factorization
	if err = QRFactor(A, &tauq, &Wred, conf); err != nil {
		return
	}
	if bits&gomas.WANTU != 0 {
		if n(U) == n(A) {
			U.Copy(A)
			if err = QRBuild(U, &tauq, &Wred, n(U), conf); err != nil {
				return
			}
		} else {
			// U is M-by-M
			eye := cmat.FloatDiagonalSource{1.0}
			U.SetFrom(&eye, cmat.SYMM)
			if err = QRMult(U, A, &tauq, &Wred, gomas.LEFT, conf); err != nil {
				return
			}
		}
	}
	R.SubMatrix(A, 0, 0, n(A), n(A))
	cmat.TriU(&R, 0)

	// bidiagonal reduce
	if err = BDReduce(&R, &tauq, &taup, &Wred, conf); err != nil {
		return
	}

	if bits&gomas.WANTU != 0 {
		Un.SubMatrix(U, 0, 0, m(A), n(A))
		if err = BDMult(&Un, &R, &tauq, &Wred, gomas.MULTQ|gomas.RIGHT, conf); err != nil {
			return
		}
		uu = U
	}
	if bits&gomas.WANTV != 0 {
		V.Copy(&R)
		if err = BDBuild(V, &taup, &Wred, m(V), gomas.WANTP, conf); err != nil {
			return
//.........這裏部分代碼省略.........