Golang FloatMatrix.SubMatrix方法代码示例

/*
 * 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 TestSubMatrixGob(t *testing.T) {
	var B, As cmat.FloatMatrix
	var network bytes.Buffer
	N := 32
	A := cmat.NewMatrix(N, N)
	zeromean := cmat.NewFloatNormSource()
	A.SetFrom(zeromean)
	As.SubMatrix(A, 3, 3, N-6, N-6)

	enc := gob.NewEncoder(&network)
	dec := gob.NewDecoder(&network)

	// encode to network
	err := enc.Encode(&As)
	if err != nil {
		t.Logf("encode error: %v\n", err)
		t.FailNow()
	}

	// decode from network
	err = dec.Decode(&B)
	if err != nil {
		t.Logf("decode error: %v\n", err)
		t.FailNow()
	}

	ar, ac := As.Size()
	br, bc := B.Size()
	t.Logf("As[%d,%d] == B[%d,%d]: %v\n", ar, ac, br, bc, B.AllClose(&As))
}

/*
 * Compute QL factorization of a M-by-N matrix A: A = Q * L.
 *
 * Arguments:
 *  A    On entry, the M-by-N matrix A, M >= N. On exit, lower triangular matrix L
 *       and the orthogonal matrix Q as product of elementary reflectors.
 *
 * tau  On exit, the scalar factors of the elemenentary reflectors.
 *
 * W    Workspace, N-by-nb matrix used for work space in blocked invocations.
 *
 * conf The blocking configuration. If nil then default blocking configuration
 *      is used. Member conf.LB defines blocking size of blocked algorithms.
 *      If it is zero then unblocked algorithm is used.
 *
 * Returns:
 *      Error indicator.
 *
 * Additional information
 *
 *  Ortogonal matrix Q is product of elementary reflectors H(k)
 *
 *    Q = H(K-1)...H(1)H(0), where K = min(M,N)
 *
 *  Elementary reflector H(k) is stored on column k of A above the diagonal with
 *  implicit unit value on diagonal entry. The vector TAU holds scalar factors
 *  of the elementary reflectors.
 *
 *  Contents of matrix A after factorization is as follow:
 *
 *    ( v0 v1 v2 v3 )   for M=6, N=4
 *    ( v0 v1 v2 v3 )
 *    ( l  v1 v2 v3 )
 *    ( l  l  v2 v3 )
 *    ( l  l  l  v3 )
 *    ( l  l  l  l  )
 *
 *  where l is element of L, vk is element of H(k).
 *
 * DecomposeQL is compatible with lapack.DGEQLF
 */
func QLFactor(A, tau, W *cmat.FloatMatrix, confs ...*gomas.Config) *gomas.Error {
	var err *gomas.Error = nil
	var tauh cmat.FloatMatrix
	conf := gomas.CurrentConf(confs...)

	if m(A) < n(A) {
		return gomas.NewError(gomas.ESIZE, "QLFactor")
	}
	wsmin := wsQL(A, 0)
	if W == nil || W.Len() < wsmin {
		return gomas.NewError(gomas.EWORK, "QLFactor", wsmin)
	}
	if tau.Len() < n(A) {
		return gomas.NewError(gomas.ESIZE, "QLFactor")
	}
	tauh.SubMatrix(tau, 0, 0, n(A), 1)
	lb := estimateLB(A, W.Len(), wsQL)
	lb = imin(lb, conf.LB)

	if lb == 0 || n(A) <= lb {
		unblockedQL(A, &tauh, W)
	} else {
		var Twork, Wrk cmat.FloatMatrix
		// block reflector T in first LB*LB elements in workspace
		// the rest, n(A)-LB*LB, is workspace for intermediate matrix operands
		Twork.SetBuf(conf.LB, conf.LB, -1, W.Data())
		Wrk.SetBuf(n(A)-conf.LB, conf.LB, -1, W.Data()[Twork.Len():])
		blockedQL(A, &tauh, &Twork, &Wrk, lb, conf)
	}
	return err
}

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

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

func TestSubMatrixJSON(t *testing.T) {
	var B, As cmat.FloatMatrix
	var network bytes.Buffer
	N := 26
	A := cmat.NewMatrix(N, N)
	zeromean := cmat.NewFloatNormSource()
	A.SetFrom(zeromean)
	As.SubMatrix(A, 3, 3, N-6, N-6)

	enc := json.NewEncoder(&network)
	dec := json.NewDecoder(&network)

	// encode to network
	err := enc.Encode(&As)
	//t.Logf("bytes: %v\n", string(network.Bytes()))
	if err != nil {
		t.Logf("encode error: %v\n", err)
		t.FailNow()
	}

	// decode from network
	err = dec.Decode(&B)
	if err != nil {
		t.Logf("decode error: %v\n", err)
		t.FailNow()
	}
	t.Logf("As == B: %v\n", B.AllClose(&As))
}

/*
 * Merge 1 by 1 block from 2 by 1 block.
 *
 *          AT
 * Abkl <-- --
 *          AB
 *
 */
func Merge2x1(ABLK, AT, AB *cmat.FloatMatrix) {
	tr, tc := AT.Size()
	br, _ := AB.Size()
	if tr > 0 {
		ABLK.SubMatrix(AT, 0, 0, tr+br, tc)
	} else {
		ABLK.SubMatrix(AB, 0, 0, br, tc)
	}
}

/*
 * Merge 1 by 1 block from 1 by 2 block.
 *
 * ABLK <--  AL | AR
 *
 */
func Merge1x2(ABLK, AL, AR *cmat.FloatMatrix) {
	lr, lc := AL.Size()
	_, rc := AR.Size()
	if lc > 0 {
		ABLK.SubMatrix(AL, 0, 0, lr, lc+rc)
	} else {
		ABLK.SubMatrix(AR, 0, 0, lr, rc)
	}
}

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

func swapCols(A *cmat.FloatMatrix, src, dst int) {
	var c0, c1 cmat.FloatMatrix
	ar, _ := A.Size()
	if src == dst || ar == 0 {
		return
	}
	c0.SubMatrix(A, 0, src, ar, 1)
	c1.SubMatrix(A, 0, dst, ar, 1)
	blasd.Swap(&c0, &c1)
}

func swapRows(A *cmat.FloatMatrix, src, dst int) {
	var r0, r1 cmat.FloatMatrix
	ar, ac := A.Size()
	if src == dst || ar == 0 {
		return
	}
	r0.SubMatrix(A, src, 0, 1, ac)
	r1.SubMatrix(A, dst, 0, 1, ac)
	blasd.Swap(&r0, &r1)
}

// Compute: (I - Y*T*Y.T).T*C = C - Y*(C.T*Y*T)
func updateHessLeftWY(C, Y1, Y2, T, V *cmat.FloatMatrix, conf *gomas.Config) {
	var C1, C2 cmat.FloatMatrix

	if C.Len() == 0 {
		return
	}
	C1.SubMatrix(C, 1, 0, m(Y1), n(C))
	C2.SubMatrix(C, 1+n(Y1), 0, m(Y2), n(C))

	updateWithQTLeft(&C1, &C2, Y1, Y2, T, V, true, conf)
}

// Compute: C*(I - Y*T*Y.T) = C - C*Y*T*Y.T = C - V*T*Y.T
func updateHessRightWY(C, Y1, Y2, T, W *cmat.FloatMatrix, conf *gomas.Config) {
	var C1, C2 cmat.FloatMatrix

	if C.Len() == 0 {
		return
	}
	C1.SubMatrix(C, 0, 1, m(C), n(Y1))
	C2.SubMatrix(C, 0, 1+n(Y1), m(C), m(Y2))

	updateWithQTRight(&C1, &C2, Y1, Y2, T, W, false, conf)
}

/*
 * Build block reflector for QL factorized matrix.
 */
func QLReflector(T, A, tau *cmat.FloatMatrix, confs ...*gomas.Config) *gomas.Error {
	var tauh cmat.FloatMatrix

	if n(T) < n(A) || m(T) < n(A) {
		return gomas.NewError(gomas.ESIZE, "QLReflector")
	}

	tauh.SubMatrix(tau, 0, 0, imin(m(A), n(A)), 1)
	unblkQLBlockReflector(T, A, &tauh)
	return nil
}