本文整理汇总了Golang中github.com/hrautila/gomas/util.Partition2x2函数的典型用法代码示例。如果您正苦于以下问题:Golang Partition2x2函数的具体用法?Golang Partition2x2怎么用?Golang Partition2x2使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了Partition2x2函数的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的Golang代码示例。
示例1: unblkBlockReflectorRQ
/*
* 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)
}
}
示例2: unblkQLBlockReflector
/*
* 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)
}
}
示例3: unblkBlockReflectorLQ
/*
* 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)
}
}
示例4: unblockedQRT
/*
* 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
}
示例5: buildQRTReflector
/*
* 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
}
示例6: blockedLQ
/*
* 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)
}
}
示例7: blockedLUnoPiv
// 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
}
示例8: blockedQL
/*
* 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)
}
}
示例9: unblockedLUnoPiv
// 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
}
示例10: unblkReduceTridiagUpper
/*
* 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)
}
}
示例11: blkReduceTridiagUpper
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)
}
}
示例12: blockedRQ
/*
* 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)
}
}
示例13: unblkReduceTridiagLower
/*
* 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)
}
}
示例14: blockedCHOL
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
}
示例15: unblkBuildQL
/*
* 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)
}
}