Golang FloatMatrix.SetIndexes方法代码示例

/*
 * 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&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 N-by-N block reflector which, together with trilu(A) represent
 *        the ortogonal matrix Q as Q = I - Y*T*Y.T where Y = trilu(A).
 *
 *  W     Workspace, P-by-nb matrix used for work space in blocked invocations.
 *
 *  flags Indicator flag
 *
 *  nb    The block size used in blocked invocations. If nb is zero default
 *        value N is used.
 *
 * Compatible with lapack.GELS (the m >= n part)
 */
func SolveQRT(B, A, T, W *matrix.FloatMatrix, flags Flags, nb int) error {
	var err error = nil
	var R, BT matrix.FloatMatrix
	if flags&TRANS != 0 {
		// Solve overdetermined system A.T*X = B

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

		// Clear bottom part of B
		B.SubMatrix(&BT, A.Cols(), 0)
		BT.SetIndexes(0.0)

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

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

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

func blockedBuildQ(A, tau, W *matrix.FloatMatrix, nb int) error {
	var err error = nil
	var ATL, ATR, ABL, ABR, AL matrix.FloatMatrix
	var A00, A01, A02, A10, A11, A12, A20, A21, A22 matrix.FloatMatrix
	var tT, tB matrix.FloatMatrix
	var t0, tau1, t2, Tw, Wrk matrix.FloatMatrix
	var mb int

	mb = A.Rows() - A.Cols()
	Twork := matrix.FloatZeros(nb, nb)

	partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, mb, 0, pBOTTOMRIGHT)
	partition2x1(
		&tT,
		&tB, tau, 0, pBOTTOM)

	// clearing of the columns of the right and setting ABR to unit diagonal
	// (only if not applying all reflectors, kb > 0)

	for ATL.Rows() > 0 && ATL.Cols() > 0 {
		repartition2x2to3x3(&ATL,
			&A00, &A01, &A02,
			&A10, &A11, &A12,
			&A20, &A21, &A22, A, nb, pTOPLEFT)
		repartition2x1to3x1(&tT,
			&t0,
			&tau1,
			&t2, tau, nb, pTOP)

		// --------------------------------------------------------

		// build block reflector from current block
		merge2x1(&AL, &A11, &A21)
		Twork.SubMatrix(&Tw, 0, 0, A11.Cols(), A11.Cols())
		unblkQRBlockReflector(&Tw, &AL, &tau1)

		// update with current block reflector (I - Y*T*Y.T)*Atrailing
		W.SubMatrix(&Wrk, 0, 0, A12.Cols(), A11.Cols())
		updateWithQT(&A12, &A22, &A11, &A21, &Tw, &Wrk, nb, false)

		// use unblocked version to compute current block
		W.SubMatrix(&Wrk, 0, 0, 1, A11.Cols())
		unblockedBuildQ(&AL, &tau1, &Wrk, 0)

		// zero upper part
		A01.SetIndexes(0.0)

		// --------------------------------------------------------
		continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &A11, &A22, A, pTOPLEFT)
		continue3x1to2x1(
			&tT,
			&tB, &t0, &tau1, tau, pTOP)
	}
	return err
}

func blockedBuildQT(A, T, W *matrix.FloatMatrix, nb int) error {
	var err error = nil
	var ATL, ATR, ABL, ABR, AL matrix.FloatMatrix
	var A00, A01, A11, A12, A21, A22 matrix.FloatMatrix
	var TTL, TTR, TBL, TBR matrix.FloatMatrix
	var T00, T01, T02, T11, T12, T22 matrix.FloatMatrix
	var tau1, Wrk matrix.FloatMatrix
	var mb int

	mb = A.Rows() - A.Cols()

	partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, mb, 0, pBOTTOMRIGHT)
	partition2x2(
		&TTL, &TTR,
		&TBL, &TBR, T, 0, 0, pBOTTOMRIGHT)

	// clearing of the columns of the right and setting ABR to unit diagonal
	// (only if not applying all reflectors, kb > 0)

	for ATL.Rows() > 0 && ATL.Cols() > 0 {
		repartition2x2to3x3(&ATL,
			&A00, &A01, nil,
			nil, &A11, &A12,
			nil, &A21, &A22, A, nb, pTOPLEFT)
		repartition2x2to3x3(&TTL,
			&T00, &T01, &T02,
			nil, &T11, &T12,
			nil, nil, &T22, T, nb, pTOPLEFT)

		// --------------------------------------------------------

		// update with current block reflector (I - Y*T*Y.T)*Atrailing
		W.SubMatrix(&Wrk, 0, 0, A12.Cols(), A11.Cols())
		updateWithQT(&A12, &A22, &A11, &A21, &T11, &Wrk, nb, false)

		// use unblocked version to compute current block
		W.SubMatrix(&Wrk, 0, 0, 1, A11.Cols())
		// elementary scalar coefficients on the diagonal, column vector
		T11.Diag(&tau1)
		merge2x1(&AL, &A11, &A21)
		// do an unblocked update to current block
		unblockedBuildQ(&AL, &tau1, &Wrk, 0)

		// zero upper part
		A01.SetIndexes(0.0)
		// --------------------------------------------------------
		continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &A11, &A22, A, pTOPLEFT)
		continue3x3to2x2(
			&TTL, &TTR,
			&TBL, &TBR, &T00, &T11, &T22, T, pTOPLEFT)
	}
	return err
}

// Make A tridiagonal, lower, non-unit matrix by clearing the strictly upper part
// of the matrix.
func TriL(A *matrix.FloatMatrix) *matrix.FloatMatrix {
	var Ac matrix.FloatMatrix
	mlen := imin(A.Rows(), A.Cols())
	for k := 1; k < mlen; k++ {
		Ac.SubMatrixOf(A, 0, k, k, 1)
		Ac.SetIndexes(0.0)
	}
	if A.Cols() > A.Rows() {
		Ac.SubMatrixOf(A, 0, A.Rows())
		Ac.SetIndexes(0.0)
	}
	return A
}

// Make A tridiagonal, upper, non-unit matrix by clearing the strictly lower part
// of the matrix.
func TriU(A *matrix.FloatMatrix) *matrix.FloatMatrix {
	var Ac matrix.FloatMatrix
	var k int
	mlen := imin(A.Rows(), A.Cols())
	for k = 0; k < mlen; k++ {
		Ac.SubMatrixOf(A, k+1, k, A.Rows()-k-1, 1)
		Ac.SetIndexes(0.0)
	}
	if A.Cols() < A.Rows() {
		Ac.SubMatrixOf(A, A.Cols(), 0)
		Ac.SetIndexes(0.0)
	}
	return A
}

// Make A tridiagonal, upper, unit matrix by clearing the strictly lower part
// of the matrix and setting diagonal elements to one.
func TriUU(A *matrix.FloatMatrix) *matrix.FloatMatrix {
	var Ac matrix.FloatMatrix
	var k int
	mlen := imin(A.Rows(), A.Cols())
	for k = 0; k < mlen; k++ {
		Ac.SubMatrixOf(A, k+1, k, A.Rows()-k-1, 1)
		Ac.SetIndexes(0.0)
		A.SetAt(k, k, 1.0)
	}
	// last element on diagonal
	A.SetAt(k, k, 1.0)
	if A.Cols() < A.Rows() {
		Ac.SubMatrixOf(A, A.Cols(), 0)
		Ac.SetIndexes(0.0)
	}
	return A
}

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

	A := matrix.FloatZeros(6, 6)
	As.SubMatrixOf(A, 1, 1, 4, 4)
	As.SetIndexes(1.0)
	partition2x2(&ATL, &ATR, &ABL, &ABR, &As, 0)
	t.Logf("ATL:\n%v\n", &ATL)

	for ATL.Rows() < As.Rows() {
		repartition2x2to3x3(&ATL,
			&A00, &a01, &A02,
			&a10, &a11, &a12,
			&A20, &a21, &A22, &As, 1)
		t.Logf("m(a12)=%d [%d], m(a11)=%d\n", a12.Cols(), a12.NumElements(), a11.NumElements())
		a11.Add(1.0)
		a21.Add(-2.0)

		continue3x3to2x2(&ATL, &ATR, &ABL, &ABR, &A00, &a11, &A22, &As)
	}
	t.Logf("A:\n%v\n", A)
}

// Build Q in place by applying elementary reflectors in reverse order to
// an implied identity matrix.  This forms Q = H(1)H(2) ... H(k)
//
// this is compatibe with lapack.DORG2R
func unblockedBuildQ(A, tau, w *matrix.FloatMatrix, kb int) error {
	var err error = nil
	var ATL, ATR, ABL, ABR matrix.FloatMatrix
	var A00, a01, A02, a10t, a11, a12t, A20, a21, A22 matrix.FloatMatrix
	var tT, tB matrix.FloatMatrix
	var t0, tau1, t2, w1 matrix.FloatMatrix
	var mb int
	var rowvec bool

	mb = A.Rows() - A.Cols()
	rowvec = tau.Rows() == 1

	partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, mb, 0, pBOTTOMRIGHT)

	if rowvec {
		partition1x2(
			&tT, &tB, tau, 0, pRIGHT)
	} else {
		partition2x1(
			&tT,
			&tB, tau, 0, pBOTTOM)
	}

	// clearing of the columns of the right and setting ABR to unit diagonal
	// (only if not applying all reflectors, kb > 0)

	for ATL.Rows() > 0 && ATL.Cols() > 0 {
		repartition2x2to3x3(&ATL,
			&A00, &a01, &A02,
			&a10t, &a11, &a12t,
			&A20, &a21, &A22, A, 1, pTOPLEFT)
		if rowvec {
			repartition1x2to1x3(&tT,
				&t0, &tau1, &t2, tau, 1, pLEFT)
		} else {
			repartition2x1to3x1(&tT,
				&t0,
				&tau1,
				&t2, tau, 1, pTOP)
		}

		// --------------------------------------------------------

		// adjust workspace to correct size
		w.SubMatrix(&w1, 0, 0, 1, a12t.Cols())
		// apply Householder reflection from left
		applyHHTo2x1(&tau1, &a21, &a12t, &A22, &w1, LEFT)

		// apply (in-place) current elementary reflector to unit vector
		a21.Scale(-tau1.Float())
		a11.SetAt(0, 0, 1.0-tau1.Float())

		// zero the upper part
		a01.SetIndexes(0.0)

		// --------------------------------------------------------
		continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &a11, &A22, A, pTOPLEFT)
		if rowvec {
			continue1x3to1x2(
				&tT, &tB, &t0, &tau1, tau, pLEFT)
		} else {
			continue3x1to2x1(
				&tT,
				&tB, &t0, &tau1, tau, pTOP)
		}
	}
	return err
}