Golang util.Repartition2x2to3x3函数代码示例

/*
 * like LAPACK/dlafrt.f
 *
 * Build block reflector T from HH reflector stored in TriLU(A) and coefficients
 * in tau.
 *
 * Q = I - Y*T*Y.T; Householder H = I - tau*v*v.T
 *
 * T = | T  0 |   z = -tau*T*Y.T*v
 *     | z  c |   c = tau
 *
 * Q = H(1)H(2)...H(k) building forward here.
 */
func unblkBlockReflectorRQ(T, A, tau *cmat.FloatMatrix) {
	var ATL, ABR cmat.FloatMatrix
	var A00, a10, A20, a11, a21, A22 cmat.FloatMatrix
	var TTL, TBR cmat.FloatMatrix
	var T00, t11, t21, T22 cmat.FloatMatrix
	var tT, tB cmat.FloatMatrix
	var t0, tau1, t2 cmat.FloatMatrix

	util.Partition2x2(
		&ATL, nil,
		nil, &ABR /**/, A, 0, 0, util.PBOTTOMRIGHT)
	util.Partition2x2(
		&TTL, nil,
		nil, &TBR /**/, T, 0, 0, util.PBOTTOMRIGHT)
	util.Partition2x1(
		&tT,
		&tB /**/, tau, 0, util.PBOTTOM)

	for m(&ATL) > 0 && n(&ATL) > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, nil, nil,
			&a10, &a11, nil,
			&A20, &a21, &A22 /**/, A, 1, util.PTOPLEFT)
		util.Repartition2x2to3x3(&TTL,
			&T00, nil, nil,
			nil, &t11, nil,
			nil, &t21, &T22 /**/, T, 1, util.PTOPLEFT)
		util.Repartition2x1to3x1(&tT,
			&t0,
			&tau1,
			&t2 /**/, tau, 1, util.PTOP)
		// --------------------------------------------------

		// t11 := tau
		tauval := tau1.Get(0, 0)
		if tauval != 0.0 {
			t11.Set(0, 0, tauval)

			// t21 := -tauval*(a21 + A20*a10)
			blasd.Axpby(&t21, &a21, 1.0, 0.0)
			blasd.MVMult(&t21, &A20, &a10, -tauval, -tauval, gomas.NONE)
			// t21 := T22*t21
			blasd.MVMultTrm(&t21, &T22, 1.0, gomas.LOWER)
		}

		// --------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, nil,
			nil, &ABR /**/, &A00, &a11, &A22, A, util.PTOPLEFT)
		util.Continue3x3to2x2(
			&TTL, nil,
			nil, &TBR /**/, &T00, &t11, &T22, T, util.PTOPLEFT)
		util.Continue3x1to2x1(
			&tT,
			&tB /**/, &t0, &tau1, tau, util.PTOP)
	}
}

/*
 * like LAPACK/dlafrt.f
 *
 * Build block reflector T from HH reflector stored in TriLU(A) and coefficients
 * in tau.
 *
 * Q = I - Y*T*Y.T; Householder H = I - tau*v*v.T
 *
 * T = | T  z |   z = -tau*T*Y.T*v
 *     | 0  c |   c = tau
 *
 * Q = H(1)H(2)...H(k) building forward here.
 */
func unblkQLBlockReflector(T, A, tau *cmat.FloatMatrix) {
	var ATL, ABR cmat.FloatMatrix
	var A00, a01, a11, A02, a12, A22 cmat.FloatMatrix
	var TTL, TBR cmat.FloatMatrix
	var T00, t11, t21, T22 cmat.FloatMatrix
	var tT, tB cmat.FloatMatrix
	var t0, tau1, t2 cmat.FloatMatrix

	util.Partition2x2(
		&ATL, nil,
		nil, &ABR, A, 0, 0, util.PBOTTOMRIGHT)
	util.Partition2x2(
		&TTL, nil,
		nil, &TBR, T, 0, 0, util.PBOTTOMRIGHT)
	util.Partition2x1(
		&tT,
		&tB, tau, 0, util.PBOTTOM)

	for m(&ATL) > 0 && n(&ATL) > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, &a01, &A02,
			nil, &a11, &a12,
			nil, nil, &A22, A, 1, util.PTOPLEFT)
		util.Repartition2x2to3x3(&TTL,
			&T00, nil, nil,
			nil, &t11, nil,
			nil, &t21, &T22, T, 1, util.PTOPLEFT)
		util.Repartition2x1to3x1(&tT,
			&t0,
			&tau1,
			&t2, tau, 1, util.PTOP)
		// --------------------------------------------------

		// t11 := tau
		tauval := tau1.Get(0, 0)
		if tauval != 0.0 {
			t11.Set(0, 0, tauval)

			// t21 := -tauval*(a12.T + &A02.T*a12)
			blasd.Axpby(&t21, &a12, 1.0, 0.0)
			blasd.MVMult(&t21, &A02, &a01, -tauval, -tauval, gomas.TRANSA)
			// t21 := T22*t01
			blasd.MVMultTrm(&t21, &T22, 1.0, gomas.LOWER)
		}

		// --------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, nil,
			nil, &ABR, &A00, &a11, &A22, A, util.PTOPLEFT)
		util.Continue3x3to2x2(
			&TTL, nil,
			nil, &TBR, &T00, &t11, &T22, T, util.PTOPLEFT)
		util.Continue3x1to2x1(
			&tT,
			&tB, &t0, &tau1, tau, util.PTOP)
	}
}

/*
 * like LAPACK/dlafrt.f
 *
 * Build block reflector T from HH reflector stored in TriLU(A) and coefficients
 * in tau.
 *
 * Q = I - Y*T*Y.T; Householder H = I - tau*v*v.T
 *
 * T = | T  z |   z = -tau*T*Y.T*v
 *     | 0  c |   c = tau
 *
 * Q = H(1)H(2)...H(k) building forward here.
 */
func unblkBlockReflectorLQ(T, A, tau *cmat.FloatMatrix) {
	var ATL, ATR, ABL, ABR cmat.FloatMatrix
	var A00, a01, A02, a11, a12, A22 cmat.FloatMatrix
	var TTL, TTR, TBL, TBR cmat.FloatMatrix
	var T00, t01, T02, t11, t12, T22 cmat.FloatMatrix
	var tT, tB cmat.FloatMatrix
	var t0, tau1, t2 cmat.FloatMatrix

	util.Partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, 0, 0, util.PTOPLEFT)
	util.Partition2x2(
		&TTL, &TTR,
		&TBL, &TBR, T, 0, 0, util.PTOPLEFT)
	util.Partition2x1(
		&tT,
		&tB, tau, 0, util.PTOP)

	for m(&ABR) > 0 && n(&ABR) > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, &a01, &A02,
			nil, &a11, &a12,
			nil, nil, &A22, A, 1, util.PBOTTOMRIGHT)
		util.Repartition2x2to3x3(&TTL,
			&T00, &t01, &T02,
			nil, &t11, &t12,
			nil, nil, &T22, T, 1, util.PBOTTOMRIGHT)
		util.Repartition2x1to3x1(&tT,
			&t0,
			&tau1,
			&t2, tau, 1, util.PBOTTOM)
		// --------------------------------------------------

		// t11 := tau
		tauval := tau1.Get(0, 0)
		if tauval != 0.0 {
			t11.Set(0, 0, tauval)

			// t01 := -tauval*(a01 + A02*a12)
			blasd.Axpby(&t01, &a01, 1.0, 0.0)
			blasd.MVMult(&t01, &A02, &a12, -tauval, -tauval, gomas.NONE)
			// t01 := T00*t01
			blasd.MVMultTrm(&t01, &T00, 1.0, gomas.UPPER)
		}

		// --------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &a11, &A22, A, util.PBOTTOMRIGHT)
		util.Continue3x3to2x2(
			&TTL, &TTR,
			&TBL, &TBR, &T00, &t11, &T22, T, util.PBOTTOMRIGHT)
		util.Continue3x1to2x1(
			&tT,
			&tB, &t0, &tau1, tau, util.PBOTTOM)
	}
}

/*
 * Unblocked QR decomposition with block reflector T.
 */
func unblockedQRT(A, T, W *cmat.FloatMatrix) *gomas.Error {
	var err *gomas.Error = nil
	var ATL, ATR, ABL, ABR cmat.FloatMatrix
	var A00, a10, a11, a12, A20, a21, A22 cmat.FloatMatrix
	var TTL, TTR, TBL, TBR cmat.FloatMatrix
	var T00, t01, T02, t11, t12, T22, w12 cmat.FloatMatrix

	util.Partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, 0, 0, util.PTOPLEFT)
	util.Partition2x2(
		&TTL, &TTR,
		&TBL, &TBR, T, 0, 0, util.PTOPLEFT)

	for m(&ABR) > 0 && n(&ABR) > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, nil, nil,
			&a10, &a11, &a12,
			&A20, &a21, &A22, A, 1, util.PBOTTOMRIGHT)
		util.Repartition2x2to3x3(&TTL,
			&T00, &t01, &T02,
			nil, &t11, &t12,
			nil, nil, &T22, T, 1, util.PBOTTOMRIGHT)

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

		computeHouseholder(&a11, &a21, &t11)

		// H*[a12 A22].T
		w12.SubMatrix(W, 0, 0, a12.Len(), 1)
		applyHouseholder2x1(&t11, &a21, &a12, &A22, &w12, gomas.LEFT)

		// update T
		tauval := t11.Get(0, 0)
		if tauval != 0.0 {
			// t01 := -tauval*(a10.T + &A20.T*a21)
			//a10.CopyTo(&t01)
			blasd.Axpby(&t01, &a10, 1.0, 0.0)
			blasd.MVMult(&t01, &A20, &a21, -tauval, -tauval, gomas.TRANSA)
			// t01 := T00*t01
			blasd.MVMultTrm(&t01, &T00, 1.0, gomas.UPPER)
		}

		// ------------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &a11, &A22, A, util.PBOTTOMRIGHT)
		util.Continue3x3to2x2(
			&TTL, &TTR,
			&TBL, &TBR, &T00, &t11, &T22, T, util.PBOTTOMRIGHT)
	}
	return err
}

/*
 * Build full block reflect T for nc columns from sequence of reflector stored in S.
 * Reflectors in S are the diagonal of T, off-diagonal values of reflector are computed
 * from elementary reflector store in lower triangular part of A.
 */
func buildQRTReflector(T, A, S *cmat.FloatMatrix, nc int, conf *gomas.Config) *gomas.Error {
	var ATL, ATR, ABL, ABR cmat.FloatMatrix
	var A00, A10, A11, A20, A21, A22 cmat.FloatMatrix
	var TTL, TTR, TBL, TBR cmat.FloatMatrix
	var T00, T01, T02, T11, T12, T22 cmat.FloatMatrix
	var SL, SR cmat.FloatMatrix
	var S00, S01, S02 cmat.FloatMatrix

	util.Partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, 0, 0, util.PTOPLEFT)
	util.Partition2x2(
		&TTL, &TTR,
		&TBL, &TBR, T, 0, 0, util.PTOPLEFT)
	util.Partition1x2(
		&SL, &SR, S, 0, util.PLEFT)

	nb := conf.LB
	for m(&ABR)-nb > 0 && n(&ABR)-nb > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, nil, nil,
			&A10, &A11, nil,
			&A20, &A21, &A22, A, nb, util.PBOTTOMRIGHT)
		util.Repartition2x2to3x3(&TTL,
			&T00, &T01, &T02,
			nil, &T11, &T12,
			nil, nil, &T22, T, nb, util.PBOTTOMRIGHT)
		util.Repartition1x2to1x3(&SL,
			&S00, &S01, &S02, S, nb, util.PRIGHT)
		// --------------------------------------------------------
		// update T01: T01 = -T00*Y1.T*Y2*T11
		//  Y1 = /A10\   Y2 = /A11\
		//       \A20/        \A21/
		//
		T11.Copy(&S01)
		updateQRTReflector(&T01, &A10, &A20, &A11, &A21, &T00, &S01, conf)
		// --------------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &A11, &A22, A, util.PBOTTOMRIGHT)
		util.Continue3x3to2x2(
			&TTL, &TTR,
			&TBL, &TBR, &T00, &T11, &T22, T, util.PBOTTOMRIGHT)
		util.Continue1x3to1x2(
			&SL, &SR, &S00, &S01, S, util.PRIGHT)
	}
	if m(&ABR) > 0 && n(&ABR) > 0 {
	}
	return nil
}

/*
 * Blocked LQ decomposition with compact WY transform. As implemented
 * in lapack.DGELQF subroutine.
 */
func blockedLQ(A, Tvec, Twork, W *cmat.FloatMatrix, lb int, conf *gomas.Config) {
	var ATL, ATR, ABL, ABR, AR cmat.FloatMatrix
	var A00, A11, A12, A21, A22 cmat.FloatMatrix
	var TT, TB cmat.FloatMatrix
	var t0, tau, t2 cmat.FloatMatrix
	var Wrk, w1 cmat.FloatMatrix

	util.Partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, 0, 0, util.PTOPLEFT)
	util.Partition2x1(
		&TT,
		&TB, Tvec, 0, util.PTOP)

	//nb := conf.LB
	for m(&ABR)-lb > 0 && n(&ABR)-lb > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, nil, nil,
			nil, &A11, &A12,
			nil, &A21, &A22, A, lb, util.PBOTTOMRIGHT)
		util.Repartition2x1to3x1(&TT,
			&t0,
			&tau,
			&t2, Tvec, lb, util.PBOTTOM)

		// current block size
		cb, rb := A11.Size()
		if rb < cb {
			cb = rb
		}
		// --------------------------------------------------------
		// decompose left side AL == /A11\
		//                           \A21/
		w1.SubMatrix(W, 0, 0, cb, 1)
		util.Merge1x2(&AR, &A11, &A12)
		unblockedLQ(&AR, &tau, &w1)

		// build block reflector
		unblkBlockReflectorLQ(Twork, &AR, &tau)

		// update A'tail i.e. A21 and A22 with A'*(I - Y*T*Y.T).T
		// compute: C - Y*(C.T*Y*T).T
		ar, ac := A21.Size()
		Wrk.SubMatrix(W, 0, 0, ar, ac)
		updateRightLQ(&A21, &A22, &A11, &A12, Twork, &Wrk, true, conf)
		// --------------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &A11, &A22, A, util.PBOTTOMRIGHT)
		util.Continue3x1to2x1(
			&TT,
			&TB, &t0, &tau, Tvec, util.PBOTTOM)
	}

	// last block with unblocked
	if m(&ABR) > 0 && n(&ABR) > 0 {
		w1.SubMatrix(W, 0, 0, m(&ABR), 1)
		unblockedLQ(&ABR, &t2, &w1)
	}
}

// blocked LU decomposition w/o pivots, FLAME LU nopivots variant 5
func blockedLUnoPiv(A *cmat.FloatMatrix, nb int, conf *gomas.Config) *gomas.Error {
	var err *gomas.Error = nil
	var ATL, ATR, ABL, ABR cmat.FloatMatrix
	var A00, A01, A02, A10, A11, A12, A20, A21, A22 cmat.FloatMatrix

	util.Partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, 0, 0, util.PTOPLEFT)

	for m(&ATL) < m(A)-nb {
		util.Repartition2x2to3x3(&ATL,
			&A00, &A01, &A02,
			&A10, &A11, &A12,
			&A20, &A21, &A22, A, nb, util.PBOTTOMRIGHT)

		// A00 = LU(A00)
		unblockedLUnoPiv(&A11, conf)
		// A12 = trilu(A00)*A12.-1  (TRSM)
		blasd.SolveTrm(&A12, &A11, 1.0, gomas.LEFT|gomas.LOWER|gomas.UNIT)
		// A21 = A21.-1*triu(A00) (TRSM)
		blasd.SolveTrm(&A21, &A11, 1.0, gomas.RIGHT|gomas.UPPER)
		// A22 = A22 - A21*A12
		blasd.Mult(&A22, &A21, &A12, -1.0, 1.0, gomas.NONE)

		util.Continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &A11, &A22, A, util.PBOTTOMRIGHT)
	}
	// last block
	if m(&ATL) < m(A) {
		unblockedLUnoPiv(&ABR, conf)
	}
	return err
}

/*
 * Blocked QR decomposition with compact WY transform.
 *
 * Compatible with lapack.DGEQRF.
 */
func blockedQL(A, Tvec, Twork, W *cmat.FloatMatrix, lb int, conf *gomas.Config) {
	var ATL, ATR, ABL, ABR, AL cmat.FloatMatrix
	var A00, A01, A10, A11, A22 cmat.FloatMatrix
	var TT, TB cmat.FloatMatrix
	var t0, tau, t2 cmat.FloatMatrix
	var Wrk, w1 cmat.FloatMatrix

	util.Partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, 0, 0, util.PBOTTOMRIGHT)
	util.Partition2x1(
		&TT,
		&TB, Tvec, 0, util.PBOTTOM)

	nb := lb
	for m(&ATL)-nb > 0 && n(&ATL)-nb > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, &A01, nil,
			&A10, &A11, nil,
			nil, nil, &A22, A, nb, util.PTOPLEFT)
		util.Repartition2x1to3x1(&TT,
			&t0,
			&tau,
			&t2, Tvec, nb, util.PTOP)

		// current block size
		cb, rb := A11.Size()
		if rb < cb {
			cb = rb
		}
		// --------------------------------------------------------
		// decompose righ side AL == /A01\
		//                           \A11/
		w1.SubMatrix(W, 0, 0, cb, 1)
		util.Merge2x1(&AL, &A01, &A11)
		unblockedQL(&AL, &tau, &w1)

		// build block reflector
		unblkQLBlockReflector(Twork, &AL, &tau)

		// update A'tail i.e. A10 and A00 with (I - Y*T*Y.T).T * A'tail
		// compute: C - Y*(C.T*Y*T).T
		ar, ac := A10.Size()
		Wrk.SubMatrix(W, 0, 0, ac, ar)
		updateQLLeft(&A10, &A00, &A11, &A01, Twork, &Wrk, true, conf)
		// --------------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &A11, &A22, A, util.PTOPLEFT)
		util.Continue3x1to2x1(
			&TT,
			&TB, &t0, &tau, Tvec, util.PTOP)
	}

	// last block with unblocked
	if m(&ATL) > 0 && n(&ATL) > 0 {
		w1.SubMatrix(W, 0, 0, n(&ATL), 1)
		unblockedQL(&ATL, &t0, &w1)
	}
}

// unblocked LU decomposition w/o pivots, FLAME LU nopivots variant 5
func unblockedLUnoPiv(A *cmat.FloatMatrix, conf *gomas.Config) *gomas.Error {
	var ATL, ATR, ABL, ABR cmat.FloatMatrix
	var A00, a01, A02, a10, a11, a12, A20, a21, A22 cmat.FloatMatrix
	var err *gomas.Error = nil

	util.Partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, 0, 0, util.PTOPLEFT)

	for m(&ATL) < m(A) {
		util.Repartition2x2to3x3(&ATL,
			&A00, &a01, &A02,
			&a10, &a11, &a12,
			&A20, &a21, &A22, A, 1, util.PBOTTOMRIGHT)

		// a21 = a21/a11
		blasd.InvScale(&a21, a11.Get(0, 0))
		// A22 = A22 - a21*a12
		blasd.MVUpdate(&A22, &a21, &a12, -1.0)

		util.Continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &a11, &A22, A, util.PBOTTOMRIGHT)
	}
	return err
}

/*
 * Reduce upper triangular matrix to tridiagonal.
 *
 * Elementary reflectors Q = H(n-1)...H(2)H(1) are stored on upper
 * triangular part of A. Reflector H(n-1) saved at column A(n) and
 * scalar multiplier to tau[n-1]. If parameter `tail` is true then
 * this function is used to reduce tail part of partially reduced
 * matrix and tau-vector partitioning is starting from last position.
 */
func unblkReduceTridiagUpper(A, tauq, W *cmat.FloatMatrix, tail bool) {
	var ATL, ABR cmat.FloatMatrix
	var A00, a01, a11, A22 cmat.FloatMatrix
	var tqT, tqB, tq0, tauq1, tq2 cmat.FloatMatrix
	var y21 cmat.FloatMatrix
	var v0 float64

	toff := 1
	if tail {
		toff = 0
	}
	util.Partition2x2(
		&ATL, nil,
		nil, &ABR, A, 0, 0, util.PBOTTOMRIGHT)
	util.Partition2x1(
		&tqT,
		&tqB, tauq, toff, util.PBOTTOM)

	for n(&ATL) > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, &a01, nil,
			nil, &a11, nil,
			nil, nil, &A22, A, 1, util.PTOPLEFT)
		util.Repartition2x1to3x1(&tqT,
			&tq0,
			&tauq1,
			&tq2, tauq, 1, util.PTOP)
		// set temp vectors for this round
		y21.SetBuf(n(&A00), 1, n(&A00), W.Data())
		// ------------------------------------------------------

		// Compute householder to zero super-diagonal entries
		computeHouseholderRev(&a01, &tauq1)
		tauqv := tauq1.Get(0, 0)

		// set superdiagonal to unit
		v0 = a01.Get(-1, 0)
		a01.Set(-1, 0, 1.0)

		// y21 := A22*a12t
		blasd.MVMultSym(&y21, &A00, &a01, tauqv, 0.0, gomas.UPPER)
		// beta := tauq*a12t*y21
		beta := tauqv * blasd.Dot(&a01, &y21)
		// y21  := y21 - 0.5*beta*a125
		blasd.Axpy(&y21, &a01, -0.5*beta)
		// A22 := A22 - a12t*y21.T - y21*a12.T
		blasd.MVUpdate2Sym(&A00, &a01, &y21, -1.0, gomas.UPPER)

		// restore superdiagonal value
		a01.Set(-1, 0, v0)
		// ------------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, nil,
			nil, &ABR, &A00, &a11, &A22, A, util.PTOPLEFT)
		util.Continue3x1to2x1(
			&tqT,
			&tqB, &tq0, &tauq1, tauq, util.PTOP)
	}
}

func blkReduceTridiagUpper(A, tauq, Y, W *cmat.FloatMatrix, lb int, conf *gomas.Config) {
	var ATL, ABR cmat.FloatMatrix
	var A00, A01, A11, A22 cmat.FloatMatrix
	var YT, YB, Y0, Y1, Y2 cmat.FloatMatrix
	var tqT, tqB, tq0, tauq1, tq2 cmat.FloatMatrix
	var v0 float64

	util.Partition2x2(
		&ATL, nil,
		nil, &ABR, A, 0, 0, util.PBOTTOMRIGHT)
	util.Partition2x1(
		&YT,
		&YB, Y, 0, util.PBOTTOM)
	util.Partition2x1(
		&tqT,
		&tqB, tauq, 1, util.PBOTTOM)

	for m(&ATL)-lb > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, &A01, nil,
			nil, &A11, nil,
			nil, nil, &A22, A, lb, util.PTOPLEFT)
		util.Repartition2x1to3x1(&YT,
			&Y0,
			&Y1,
			&Y2, Y, lb, util.PTOP)
		util.Repartition2x1to3x1(&tqT,
			&tq0,
			&tauq1,
			&tq2, tauq, lb, util.PTOP)
		// ------------------------------------------------------
		unblkBuildTridiagUpper(&ATL, &tauq1, &YT, W)

		// set subdiagonal entry to unit
		v0 = A01.Get(-1, 0)
		A01.Set(-1, 0, 1.0)

		// A22 := A22 - A01*Y0.T - Y0*A01.T
		blasd.Update2Sym(&A00, &A01, &Y0, -1.0, 1.0, gomas.UPPER, conf)

		// restore subdiagonal entry
		A01.Set(-1, 0, v0)
		// ------------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, nil,
			nil, &ABR, &A00, &A11, &A22, A, util.PTOPLEFT)
		util.Continue3x1to2x1(
			&YT,
			&YB, &Y0, &Y1, Y, util.PTOP)
		util.Continue3x1to2x1(
			&tqT,
			&tqB, &tq0, &tauq1, tauq, util.PTOP)
	}

	if m(&ATL) > 0 {
		unblkReduceTridiagUpper(&ATL, &tqT, W, true)
	}
}

/*
 * Blocked RQ decomposition with compact WY transform. As implemented
 * in lapack.DGERQF subroutine.
 */
func blockedRQ(A, Tvec, Twork, W *cmat.FloatMatrix, lb int, conf *gomas.Config) {
	var ATL, ABR, AL cmat.FloatMatrix
	var A00, A01, A10, A11, A22 cmat.FloatMatrix
	var TT, TB cmat.FloatMatrix
	var t0, tau, t2 cmat.FloatMatrix
	var Wrk, w1 cmat.FloatMatrix

	util.Partition2x2(
		&ATL, nil,
		nil, &ABR /**/, A, 0, 0, util.PBOTTOMRIGHT)
	util.Partition2x1(
		&TT,
		&TB /**/, Tvec, 0, util.PBOTTOM)

	for m(&ATL)-lb > 0 && n(&ATL)-lb > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, &A01, nil,
			&A10, &A11, nil,
			nil, nil, &A22 /**/, A, lb, util.PTOPLEFT)
		util.Repartition2x1to3x1(&TT,
			&t0,
			&tau,
			&t2 /**/, Tvec, n(&A11), util.PTOP)

		// current block size
		cb, rb := A11.Size()
		if rb < cb {
			cb = rb
		}
		// --------------------------------------------------------
		// decompose left side AL == ( A10 A11 )
		w1.SubMatrix(W, 0, 0, cb, 1)
		util.Merge1x2(&AL, &A10, &A11)
		unblockedRQ(&AL, &tau, &w1)

		// build block reflector
		unblkBlockReflectorRQ(Twork, &AL, &tau)

		// compute: (A00 A01)(I - Y*T*Y.T)
		ar, ac := A01.Size()
		Wrk.SubMatrix(W, 0, 0, ar, ac)
		updateRightRQ(&A01, &A00, &A11, &A10, Twork, &Wrk, false, conf)
		// --------------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, nil,
			nil, &ABR, &A00, &A11, &A22, A, util.PTOPLEFT)
		util.Continue3x1to2x1(
			&TT,
			&TB, &t0, &tau, Tvec, util.PTOP)
	}

	// last block with unblocked
	if m(&ATL) > 0 && n(&ATL) > 0 {
		w1.SubMatrix(W, 0, 0, m(&ATL), 1)
		unblockedRQ(&ATL, &TT, &w1)
	}
}

/*
 * Tridiagonal reduction of LOWER triangular symmetric matrix, zero elements below 1st
 * subdiagonal:
 *
 *   A =  (1 - tau*u*u.t)*A*(1 - tau*u*u.T)
 *     =  (I - tau*( 0   0   )) (a11 a12) (I - tau*( 0  0   ))
 *        (        ( 0  u*u.t)) (a21 A22) (        ( 0 u*u.t))
 *
 *  a11, a12, a21 not affected
 *
 *  from LEFT:
 *    A22 = A22 - tau*u*u.T*A22
 *  from RIGHT:
 *    A22 = A22 - tau*A22*u.u.T
 *
 *  LEFT and RIGHT:
 *    A22   = A22 - tau*u*u.T*A22 - tau*(A22 - tau*u*u.T*A22)*u*u.T
 *          = A22 - tau*u*u.T*A22 - tau*A22*u*u.T + tau*tau*u*u.T*A22*u*u.T
 *    [x    = tau*A22*u (vector)]  (SYMV)
 *    A22   = A22 - u*x.T - x*u.T + tau*u*u.T*x*u.T
 *    [beta = tau*u.T*x (scalar)]  (DOT)
 *          = A22 - u*x.T - x*u.T + beta*u*u.T
 *          = A22 - u*(x - 0.5*beta*u).T - (x - 0.5*beta*u)*u.T
 *    [w    = x - 0.5*beta*u]      (AXPY)
 *          = A22 - u*w.T - w*u.T  (SYR2)
 *
 * Result of reduction for N = 5:
 *    ( d  .  .  . . )
 *    ( e  d  .  . . )
 *    ( v1 e  d  . . )
 *    ( v1 v2 e  d . )
 *    ( v1 v2 v3 e d )
 */
func unblkReduceTridiagLower(A, tauq, W *cmat.FloatMatrix) {
	var ATL, ABR cmat.FloatMatrix
	var A00, a11, a21, A22 cmat.FloatMatrix
	var tqT, tqB, tq0, tauq1, tq2 cmat.FloatMatrix
	var y21 cmat.FloatMatrix
	var v0 float64

	util.Partition2x2(
		&ATL, nil,
		nil, &ABR, A, 0, 0, util.PTOPLEFT)
	util.Partition2x1(
		&tqT,
		&tqB, tauq, 0, util.PTOP)

	for m(&ABR) > 0 && n(&ABR) > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, nil, nil,
			nil, &a11, nil,
			nil, &a21, &A22, A, 1, util.PBOTTOMRIGHT)
		util.Repartition2x1to3x1(&tqT,
			&tq0,
			&tauq1,
			&tq2, tauq, 1, util.PBOTTOM)
		// set temp vectors for this round
		y21.SetBuf(n(&A22), 1, n(&A22), W.Data())
		// ------------------------------------------------------

		// Compute householder to zero subdiagonal entries
		computeHouseholderVec(&a21, &tauq1)
		tauqv := tauq1.Get(0, 0)

		// set subdiagonal to unit
		v0 = a21.Get(0, 0)
		a21.Set(0, 0, 1.0)

		// y21 := tauq*A22*a21
		blasd.MVMultSym(&y21, &A22, &a21, tauqv, 0.0, gomas.LOWER)
		// beta := tauq*a21.T*y21
		beta := tauqv * blasd.Dot(&a21, &y21)
		// y21  := y21 - 0.5*beta*a21
		blasd.Axpy(&y21, &a21, -0.5*beta)
		// A22 := A22 - a21*y21.T - y21*a21.T
		blasd.MVUpdate2Sym(&A22, &a21, &y21, -1.0, gomas.LOWER)

		// restore subdiagonal
		a21.Set(0, 0, v0)
		// ------------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, nil,
			nil, &ABR, &A00, &a11, &A22, A, util.PBOTTOMRIGHT)
		util.Continue3x1to2x1(
			&tqT,
			&tqB, &tq0, &tauq1, tauq, util.PBOTTOM)
	}
}

func blockedCHOL(A *cmat.FloatMatrix, flags int, conf *gomas.Config) *gomas.Error {
	var err, firstErr *gomas.Error
	var ATL, ATR, ABL, ABR cmat.FloatMatrix
	var A00, A01, A02, A10, A11, A12, A20, A21, A22 cmat.FloatMatrix

	nb := conf.LB
	err = nil
	firstErr = nil
	util.Partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, 0, 0, util.PTOPLEFT)

	for m(A)-m(&ATL) > nb {
		util.Repartition2x2to3x3(&ATL,
			&A00, &A01, &A02,
			&A10, &A11, &A12,
			&A20, &A21, &A22, A, nb, util.PBOTTOMRIGHT)

		if flags&gomas.LOWER != 0 {
			// A11 = chol(A11)
			err = unblockedLowerCHOL(&A11, flags, m(&ATL))
			// A21 = A21 * tril(A11).-1
			blasd.SolveTrm(&A21, &A11, 1.0, gomas.RIGHT|gomas.LOWER|gomas.TRANSA, conf)
			// A22 = A22 - A21*A21.T
			blasd.UpdateSym(&A22, &A21, -1.0, 1.0, gomas.LOWER, conf)
		} else {
			// A11 = chol(A11)
			err = unblockedUpperCHOL(&A11, flags, m(&ATL))
			// A12 = triu(A11).-1 * A12
			blasd.SolveTrm(&A12, &A11, 1.0, gomas.UPPER|gomas.TRANSA, conf)
			// A22 = A22 - A12.T*A12
			blasd.UpdateSym(&A22, &A12, -1.0, 1.0, gomas.UPPER|gomas.TRANSA, conf)
		}
		if err != nil && firstErr == nil {
			firstErr = err
		}

		util.Continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &A11, &A22, A, util.PBOTTOMRIGHT)
	}

	if m(&ATL) < m(A) {
		// last block
		if flags&gomas.LOWER != 0 {
			unblockedLowerCHOL(&ABR, flags, 0)
		} else {
			unblockedUpperCHOL(&ABR, flags, 0)
		}
	}
	return firstErr
}

/*
 * Unblocked code for generating M by N matrix Q with orthogonal columns which
 * are defined as the last N columns of the product of K first elementary
 * reflectors.
 *
 * Parameter nk is last nk elementary reflectors that are not used in computing
 * the matrix Q. Parameter mk length of the first unused elementary reflectors
 * First nk columns are zeroed and subdiagonal mk-nk is set to unit.
 *
 * Compatible with lapack.DORG2L subroutine.
 */
func unblkBuildQL(A, Tvec, W *cmat.FloatMatrix, mk, nk int, mayClear bool) {
	var ATL, ATR, ABL, ABR cmat.FloatMatrix
	var A00, a01, a10, a11, a21, A22 cmat.FloatMatrix
	var tT, tB cmat.FloatMatrix
	var t0, tau1, t2, w12, D cmat.FloatMatrix

	// (mk, nk) = (rows, columns) of upper left partition
	util.Partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, mk, nk, util.PTOPLEFT)
	util.Partition2x1(
		&tT,
		&tB, Tvec, nk, util.PTOP)

	// zero the left side
	if nk > 0 && mayClear {
		blasd.Scale(&ABL, 0.0)
		blasd.Scale(&ATL, 0.0)
		D.Diag(&ATL, nk-mk)
		blasd.Add(&D, 1.0)
	}

	for m(&ABR) > 0 && n(&ABR) > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, &a01, nil,
			&a10, &a11, nil,
			nil, &a21, &A22, A, 1, util.PBOTTOMRIGHT)
		util.Repartition2x1to3x1(&tT,
			&t0,
			&tau1,
			&t2, Tvec, 1, util.PBOTTOM)
		// ------------------------------------------------------
		w12.SubMatrix(W, 0, 0, a10.Len(), 1)
		applyHouseholder2x1(&tau1, &a01, &a10, &A00, &w12, gomas.LEFT)

		blasd.Scale(&a01, -tau1.Get(0, 0))
		a11.Set(0, 0, 1.0-tau1.Get(0, 0))

		// zero bottom elements
		blasd.Scale(&a21, 0.0)
		// ------------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &a11, &A22, A, util.PBOTTOMRIGHT)
		util.Continue3x1to2x1(
			&tT,
			&tB, &t0, &tau1, Tvec, util.PBOTTOM)
	}
}