本文整理匯總了Golang中github.com/henrylee2cn/algorithm/matrix.FloatMatrix.CopyTo方法的典型用法代碼示例。如果您正苦於以下問題：Golang FloatMatrix.CopyTo方法的具體用法？Golang FloatMatrix.CopyTo怎麽用？Golang FloatMatrix.CopyTo使用的例子？那麽, 這裏精選的方法代碼示例或許可以為您提供幫助。您也可以進一步了解該方法所在類github.com/henrylee2cn/algorithm/matrix.FloatMatrix的用法示例。
在下文中一共展示了FloatMatrix.CopyTo方法的14個代碼示例，這些例子默認根據受歡迎程度排序。您可以為喜歡或者感覺有用的代碼點讚，您的評價將有助於係統推薦出更棒的Golang代碼示例。
示例1: applyHHTo2x1
/*
* Applies a real elementary reflector H to a real m by n matrix A,
* from either the left or the right. H is represented in the form
*
* H = I - tau * ( 1 ) * ( 1 v.T )
* ( v )
*
* where tau is a real scalar and v is a real vector.
*
* If tau = 0, then H is taken to be the unit matrix.
*
* A is /a1\ a1 := a1 - w1
* \A2/ A2 := A2 - v*w1
* w1 := tau*(a1 + A2.T*v) if side == LEFT
* := tau*(a1 + A2*v) if side == RIGHT
*
* Intermediate work space w1 required as parameter, no allocation.
*/
func applyHHTo2x1(tau, v, a1, A2, w1 *matrix.FloatMatrix, flags Flags) {
tval := tau.GetAt(0, 0)
if tval == 0.0 {
return
}
// maybe with Scale(0.0), Axpy(w1, a1, 1.0)
a1.CopyTo(w1)
if flags&LEFT != 0 {
// w1 = a1 + A2.T*v
MVMult(w1, A2, v, 1.0, 1.0, TRANSA)
} else {
// w1 = a1 + A2*v
MVMult(w1, A2, v, 1.0, 1.0, NOTRANS)
}
// w1 = tau*w1
Scale(w1, tval)
// a1 = a1 - w1
a1.Minus(w1)
// A2 = A2 - v*w1
if flags&LEFT != 0 {
MVRankUpdate(A2, v, w1, -1.0)
} else {
MVRankUpdate(A2, w1, v, -1.0)
}
}示例2: blkUpperLDL
func blkUpperLDL(A, W *matrix.FloatMatrix, p *pPivots, nb int) (err error) {
var ATL, ATR, ABL, ABR matrix.FloatMatrix
var A00, A01, A02, A11, A12, A22 matrix.FloatMatrix
var D1, wrk matrix.FloatMatrix
var pT, pB, p0, p1, p2 pPivots
err = nil
partition2x2(
&ATL, &ATR,
&ABL, &ABR, A, 0, 0, pBOTTOMRIGHT)
partitionPivot2x1(
&pT,
&pB, p, 0, pBOTTOM)
for ATL.Rows() > 0 {
repartition2x2to3x3(&ATL,
&A00, &A01, &A02,
nil, &A11, &A12,
nil, nil, &A22, A, nb, pTOPLEFT)
repartPivot2x1to3x1(&pT,
&p0, &p1, &p2 /**/, p, nb, pTOP)
// --------------------------------------------------------
// A11 = LDL(A11)
err = unblkUpperLDL(&A11, &p1)
if err != nil {
return
}
applyColPivots(&A01, &p1, 0, BACKWARD)
applyRowPivots(&A12, &p1, 0, BACKWARD)
scalePivots(&p1, ATL.Rows()-A11.Rows())
A11.Diag(&D1)
// A01 = A01*A11.-T
SolveTrm(&A01, &A11, 1.0, UPPER|UNIT|RIGHT|TRANSA)
// A01 = A01*D1.-1
SolveDiag(&A01, &D1, RIGHT)
// W = D1*U01.T = U01*D1
W.SubMatrix(&wrk, 0, 0, A01.Rows(), nb)
A01.CopyTo(&wrk)
MultDiag(&wrk, &D1, RIGHT)
// A00 = A00 - U01*D1*U01.T = A22 - U01*W.T
UpdateTrm(&A00, &A01, &wrk, -1.0, 1.0, UPPER|TRANSB)
// ---------------------------------------------------------
continue3x3to2x2(
&ATL, &ATR,
&ABL, &ABR, &A00, &A11, &A22, A, pTOPLEFT)
contPivot3x1to2x1(
&pT,
&pB, &p0, &p1, p, pTOP)
}
return
}示例3: unblkQRBlockReflector
/*
* 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 unblkQRBlockReflector(T, A, tau *matrix.FloatMatrix) {
var ATL, ATR, ABL, ABR matrix.FloatMatrix
var A00, a10, a11, A20, a21, A22 matrix.FloatMatrix
var TTL, TTR, TBL, TBR matrix.FloatMatrix
var T00, t01, T02, t11, t12, T22 matrix.FloatMatrix
var tT, tB matrix.FloatMatrix
var t0, tau1, t2 matrix.FloatMatrix
partition2x2(
&ATL, &ATR,
&ABL, &ABR, A, 0, 0, pTOPLEFT)
partition2x2(
&TTL, &TTR,
&TBL, &TBR, T, 0, 0, pTOPLEFT)
partition2x1(
&tT,
&tB, tau, 0, pTOP)
for ABR.Rows() > 0 && ABR.Cols() > 0 {
repartition2x2to3x3(&ATL,
&A00, nil, nil,
&a10, &a11, nil,
&A20, &a21, &A22, A, 1, pBOTTOMRIGHT)
repartition2x2to3x3(&TTL,
&T00, &t01, &T02,
nil, &t11, &t12,
nil, nil, &T22, T, 1, pBOTTOMRIGHT)
repartition2x1to3x1(&tT,
&t0,
&tau1,
&t2, tau, 1, pBOTTOM)
// --------------------------------------------------
// t11 := tau
tauval := tau1.GetAt(0, 0)
if tauval != 0.0 {
t11.SetAt(0, 0, tauval)
// t01 := a10.T + &A20.T*a21
a10.CopyTo(&t01)
MVMult(&t01, &A20, &a21, -tauval, -tauval, TRANSA)
// t01 := T00*t01
MVMultTrm(&t01, &T00, UPPER)
//t01.Scale(-tauval)
}
// --------------------------------------------------
continue3x3to2x2(
&ATL, &ATR,
&ABL, &ABR, &A00, &a11, &A22, A, pBOTTOMRIGHT)
continue3x3to2x2(
&TTL, &TTR,
&TBL, &TBR, &T00, &t11, &T22, T, pBOTTOMRIGHT)
continue3x1to2x1(
&tT,
&tB, &t0, &tau1, tau, pBOTTOM)
}
}示例4: unblockedQRT
/*
* Unblocked QR decomposition with block reflector T.
*/
func unblockedQRT(A, T *matrix.FloatMatrix) {
var ATL, ATR, ABL, ABR matrix.FloatMatrix
var A00, a10, a11, a12, A20, a21, A22 matrix.FloatMatrix
var TTL, TTR, TBL, TBR matrix.FloatMatrix
var T00, t01, T02, t11, t12, T22 matrix.FloatMatrix
//As.SubMatrixOf(A, 0, 0, mlen, nb)
partition2x2(
&ATL, &ATR,
&ABL, &ABR, A, 0, 0, pTOPLEFT)
partition2x2(
&TTL, &TTR,
&TBL, &TBR, T, 0, 0, pTOPLEFT)
for ABR.Rows() > 0 && ABR.Cols() > 0 {
repartition2x2to3x3(&ATL,
&A00, nil, nil,
&a10, &a11, &a12,
&A20, &a21, &A22, A, 1, pBOTTOMRIGHT)
repartition2x2to3x3(&TTL,
&T00, &t01, &T02,
nil, &t11, &t12,
nil, nil, &T22, T, 1, pBOTTOMRIGHT)
// ------------------------------------------------------
computeHouseholder(&a11, &a21, &t11, LEFT)
// H*[a12 A22].T
applyHouseholder(&t11, &a21, &a12, &A22, LEFT)
// update T
tauval := t11.GetAt(0, 0)
if tauval != 0.0 {
// t01 := -tauval*(a10.T + &A20.T*a21)
a10.CopyTo(&t01)
MVMult(&t01, &A20, &a21, -tauval, -tauval, TRANSA)
// t01 := T00*t01
MVMultTrm(&t01, &T00, UPPER)
}
// ------------------------------------------------------
continue3x3to2x2(
&ATL, &ATR,
&ABL, &ABR, &A00, &a11, &A22, A, pBOTTOMRIGHT)
continue3x3to2x2(
&TTL, &TTR,
&TBL, &TBR, &T00, &t11, &T22, T, pBOTTOMRIGHT)
}
}示例5: blkUpperLDLnoPiv
func blkUpperLDLnoPiv(A, W *matrix.FloatMatrix, nb int) (err error) {
var ATL, ATR, ABL, ABR matrix.FloatMatrix
var A00, A01, A02, A11, A12, A22 matrix.FloatMatrix
var D1, wrk matrix.FloatMatrix
err = nil
partition2x2(
&ATL, &ATR,
&ABL, &ABR, A, 0, 0, pBOTTOMRIGHT)
for ATL.Rows() > 0 {
repartition2x2to3x3(&ATL,
&A00, &A01, &A02,
nil, &A11, &A12,
nil, nil, &A22, A, nb, pTOPLEFT)
// --------------------------------------------------------
// A11 = LDL(A11)
unblkUpperLDLnoPiv(&A11)
A11.Diag(&D1)
// A01 = A01*A11.-T
SolveTrm(&A01, &A11, 1.0, UPPER|UNIT|RIGHT|TRANSA)
// A01 = A01*D1.-1
SolveDiag(&A01, &D1, RIGHT)
// W = D1*U01.T = U01*D1
W.SubMatrix(&wrk, 0, 0, A01.Rows(), nb)
A01.CopyTo(&wrk)
MultDiag(&wrk, &D1, RIGHT)
// A00 = A00 - U01*D1*U01.T = A22 - U01*W.T
UpdateTrm(&A00, &A01, &wrk, -1.0, 1.0, UPPER|TRANSB)
// ---------------------------------------------------------
continue3x3to2x2(
&ATL, &ATR,
&ABL, &ABR, &A00, &A11, &A22, A, pTOPLEFT)
}
return
}示例6: findAndBuildPivot
func findAndBuildPivot(AL, AR, WL, WR *matrix.FloatMatrix, k int) int {
var dg, acol, wcol, wrow matrix.FloatMatrix
// updated diagonal values on last column of workspace
WR.SubMatrix(&dg, 0, WR.Cols()-1, AR.Rows(), 1)
// find on-diagonal maximun value
dmax := IAMax(&dg)
//fmt.Printf("dmax=%d, val=%e\n", dmax, dg.GetAt(dmax, 0))
// copy to first column of WR and update with factorized columns
WR.SubMatrix(&wcol, 0, 0, WR.Rows(), 1)
if dmax == 0 {
AR.SubMatrix(&acol, 0, 0, AR.Rows(), 1)
acol.CopyTo(&wcol)
} else {
AR.SubMatrix(&acol, dmax, 0, 1, dmax+1)
acol.CopyTo(&wcol)
if dmax < AR.Rows()-1 {
var wrst matrix.FloatMatrix
WR.SubMatrix(&wrst, dmax, 0, wcol.Rows()-dmax, 1)
AR.SubMatrix(&acol, dmax, dmax, AR.Rows()-dmax, 1)
acol.CopyTo(&wrst)
}
}
if k > 0 {
WL.SubMatrix(&wrow, dmax, 0, 1, WL.Cols())
//fmt.Printf("update with wrow:%v\n", &wrow)
//fmt.Printf("update wcol\n%v\n", &wcol)
MVMult(&wcol, AL, &wrow, -1.0, 1.0, NOTRANS)
//fmt.Printf("updated wcol:\n%v\n", &wcol)
}
if dmax > 0 {
// pivot column in workspace
t0 := WR.GetAt(0, 0)
WR.SetAt(0, 0, WR.GetAt(dmax, 0))
WR.SetAt(dmax, 0, t0)
// pivot on diagonal
t0 = dg.GetAt(0, 0)
dg.SetAt(0, 0, dg.GetAt(dmax, 0))
dg.SetAt(dmax, 0, t0)
}
return dmax
}示例7: unblkBoundedBKUpper
func unblkBoundedBKUpper(A, wrk *matrix.FloatMatrix, p *pPivots, ncol int) (error, int) {
var err error
var ATL, ATR, ABL, ABR matrix.FloatMatrix
var A00, a01, A02, a11, a12t, A22, a11inv matrix.FloatMatrix
var w00, w01, w11 matrix.FloatMatrix
var cwrk matrix.FloatMatrix
var wx, Ax, wz matrix.FloatMatrix
var pT, pB, p0, p1, p2 pPivots
err = nil
nc := 0
if ncol > A.Cols() {
ncol = A.Cols()
}
partition2x2(
&ATL, &ATR,
&ABL, &ABR, A, 0, 0, pBOTTOMRIGHT)
partitionPivot2x1(
&pT,
&pB, p, 0, pBOTTOM)
// permanent working space for symmetric inverse of a11
wrk.SubMatrix(&a11inv, wrk.Rows()-2, 0, 2, 2)
a11inv.SetAt(0, 1, -1.0)
a11inv.SetAt(1, 0, -1.0)
for ATL.Cols() > 0 && nc < ncol {
partition2x2(
&w00, &w01,
nil, &w11, wrk, nc, nc, pBOTTOMRIGHT)
merge1x2(&wx, &w00, &w01)
merge1x2(&Ax, &ATL, &ATR)
//fmt.Printf("ATL:\n%v\n", &ATL)
r, np := findAndBuildBKPivotUpper(&ATL, &ATR, &w00, &w01, nc)
//fmt.Printf("[w00;w01]:\n%v\n", &wx)
//fmt.Printf("after find: r=%d, np=%d, ncol=%d, nc=%d\n", r, np, ncol, nc)
w00.SubMatrix(&wz, 0, w00.Cols()-2, w00.Rows(), 2)
if np > ncol-nc {
// next pivot does not fit into ncol columns, restore last column,
// return with number of factorized columns
return err, nc
}
if r != -1 {
// pivoting needed; np == 1, last row; np == 2; next to last rows
nrow := ATL.Rows() - np
applyBKPivotSym(&ATL, nrow, r, UPPER)
// swap left hand rows to get correct updates
swapRows(&ATR, nrow, r)
swapRows(&w01, nrow, r)
if np == 2 {
/* pivot block on diagonal; -1,-1
* [r, r] | [r ,-1]
* ---------------- 2-by-2 pivot, swapping [1,0] and [r,0]
* [r,-1] | [-1,-1]
*/
t0 := w00.GetAt(-2, -1)
tr := w00.GetAt(r, -1)
//fmt.Printf("nc=%d, t0=%e, tr=%e\n", nc, t0, tr)
w00.SetAt(-2, -1, tr)
w00.SetAt(r, -1, t0)
// interchange diagonal entries on w11[:,1]
t0 = w00.GetAt(-2, -2)
tr = w00.GetAt(r, -2)
w00.SetAt(-2, -2, tr)
w00.SetAt(r, -2, t0)
//fmt.Printf("wrk:\n%v\n", &wz)
}
//fmt.Printf("pivoted A:\n%v\n", &Ax)
//fmt.Printf("pivoted wrk:\n%v\n", &wx)
}
// repartition according the pivot size
repartition2x2to3x3(&ATL,
&A00, &a01, &A02,
nil, &a11, &a12t,
nil, nil, &A22 /**/, A, np, pTOPLEFT)
repartPivot2x1to3x1(&pT,
&p0,
&p1,
&p2 /**/, p, np, pTOP)
// ------------------------------------------------------------
wlc := w00.Cols() - np
//wlr := w00.Rows() - 1
w00.SubMatrix(&cwrk, 0, wlc, a01.Rows(), np)
if np == 1 {
//fmt.Printf("wz:\n%v\n", &wz)
//fmt.Printf("a11 <-- %e\n", w00.GetAt(a01.Rows(), wlc))
//w00.SubMatrix(&cwrk, 0, wlc-np+1, a01.Rows(), np)
a11.SetAt(0, 0, w00.GetAt(a01.Rows(), wlc))
// a21 = a21/a11
//fmt.Printf("np == 1: pre-update a01\n%v\n", &a01)
cwrk.CopyTo(&a01)
InvScale(&a01, a11.Float())
//fmt.Printf("np == 1: cwrk\n%v\na21\n%v\n", &cwrk, &a21)
// store pivot point relative to original matrix
//.........這裏部分代碼省略.........
示例8: findAndBuildBKPivotUpper
func findAndBuildBKPivotUpper(AL, AR, WL, WR *matrix.FloatMatrix, k int) (int, int) {
var r, q int
var rcol, qrow, src, wk, wkp1, wrow matrix.FloatMatrix
lc := AL.Cols() - 1
wc := WL.Cols() - 1
lr := AL.Rows() - 1
// Copy AR[:,lc] to WR[:,wc] and update with WL[0:]
AL.SubMatrix(&src, 0, lc, AL.Rows(), 1)
WL.SubMatrix(&wk, 0, wc, AL.Rows(), 1)
src.CopyTo(&wk)
if k > 0 {
WR.SubMatrix(&wrow, lr, 0, 1, WR.Cols())
//fmt.Printf("wrow: %v\n", &wrow)
MVMult(&wk, AR, &wrow, -1.0, 1.0, NOTRANS)
//fmt.Printf("wk after update:\n%v\n", &wk)
}
if AL.Rows() == 1 {
return -1, 1
}
amax := math.Abs(WL.GetAt(lr, wc))
// find max off-diagonal on first column.
WL.SubMatrix(&rcol, 0, wc, lr, 1)
//fmt.Printf("rcol:\n%v\n", &rcol)
// r is row index and rmax is its absolute value
r = IAMax(&rcol)
rmax := math.Abs(rcol.GetAt(r, 0))
//fmt.Printf("r=%d, amax=%e, rmax=%e\n", r, amax, rmax)
if amax >= bkALPHA*rmax {
// no pivoting, 1x1 diagonal
return -1, 1
}
// Now we need to copy row r to WR[:,wc-1] and update it
WL.SubMatrix(&wkp1, 0, wc-1, AL.Rows(), 1)
if r > 0 {
// above the diagonal part of AL
AL.SubMatrix(&qrow, 0, r, r, 1)
qrow.CopyTo(&wkp1)
}
//fmt.Printf("m(AR)=%d, r=%d, qrow: %v\n", AL.Rows(), r, &qrow)
var wkr matrix.FloatMatrix
AL.SubMatrix(&qrow, r, r, 1, AL.Rows()-r)
wkp1.SubMatrix(&wkr, r, 0, AL.Rows()-r, 1)
qrow.CopyTo(&wkr)
//fmt.Printf("m(AR)=%d, r=%d, qrow: %v\n", AR.Rows(), r, &qrow)
if k > 0 {
// update wkp1
WR.SubMatrix(&wrow, r, 0, 1, WR.Cols())
//fmt.Printf("initial wpk1:\n%v\n", &wkp1)
MVMult(&wkp1, AR, &wrow, -1.0, 1.0, NOTRANS)
}
//fmt.Printf("updated wpk1:\n%v\n", &wkp1)
// set on-diagonal entry to zero to avoid hitting it.
p1 := wkp1.GetAt(r, 0)
wkp1.SetAt(r, 0, 0.0)
// max off-diagonal on r'th column/row at index q
q = IAMax(&wkp1)
qmax := math.Abs(wkp1.GetAt(q, 0))
wkp1.SetAt(r, 0, p1)
//fmt.Printf("blk: r=%d, q=%d, amax=%e, rmax=%e, qmax=%e\n", r, q, amax, rmax, qmax)
if amax >= bkALPHA*rmax*(rmax/qmax) {
// no pivoting, 1x1 diagonal
return -1, 1
}
// if q == r then qmax is not off-diagonal, qmax == WR[r,1] and
// we get 1x1 pivot as following is always true
if math.Abs(WL.GetAt(r, wc-1)) >= bkALPHA*qmax {
// 1x1 pivoting and interchange with k, r
// pivot row in column WR[:,1] to W[:,0]
//p1 := WL.GetAt(r, wc-1)
WL.SubMatrix(&src, 0, wc-1, AL.Rows(), 1)
WL.SubMatrix(&wkp1, 0, wc, AL.Rows(), 1)
src.CopyTo(&wkp1)
wkp1.SetAt(-1, 0, src.GetAt(r, 0))
wkp1.SetAt(r, 0, src.GetAt(-1, 0))
return r, 1
} else {
// 2x2 pivoting and interchange with k+1, r
return r, 2
}
return -1, 1
}示例9: unblkBoundedBKLower
//.........這裏部分代碼省略.........
swapRows(&ABL, np-1, r)
swapRows(&w10, np-1, r)
//ABL.SubMatrix(&s, np-1, 0, 1, ABL.Cols())
//ABL.SubMatrix(&d, r, 0, 1, ABL.Cols())
//Swap(&s, &d)
//w10.SubMatrix(&s, np-1, 0, 1, w10.Cols())
//w10.SubMatrix(&d, r, 0, 1, w10.Cols())
//Swap(&s, &d)
if np == 2 {
/*
* [0,0] | [r,0]
* a11 == ------------- 2-by-2 pivot, swapping [1,0] and [r,0]
* [r,0] | [r,r]
*/
t0 := w11.GetAt(1, 0)
tr := w11.GetAt(r, 0)
//fmt.Printf("nc=%d, t0=%e, tr=%e\n", nc, t0, tr)
w11.SetAt(1, 0, tr)
w11.SetAt(r, 0, t0)
// interchange diagonal entries on w11[:,1]
t0 = w11.GetAt(1, 1)
tr = w11.GetAt(r, 1)
w11.SetAt(1, 1, tr)
w11.SetAt(r, 1, t0)
}
//fmt.Printf("pivoted A:\n%v\n", A)
//fmt.Printf("pivoted wrk:\n%v\n", wrk)
}
// repartition according the pivot size
repartition2x2to3x3(&ATL,
&A00, nil, nil,
&a10t, &a11, nil,
&A20, &a21, &A22 /**/, A, np, pBOTTOMRIGHT)
repartPivot2x1to3x1(&pT,
&p0,
&p1,
&p2 /**/, p, np, pBOTTOM)
// ------------------------------------------------------------
if np == 1 {
//
w11.SubMatrix(&cwrk, np, 0, a21.Rows(), np)
a11.SetAt(0, 0, w11.GetAt(0, 0))
// a21 = a21/a11
//fmt.Printf("np == 1: pre-update a21\n%v\n", &a21)
cwrk.CopyTo(&a21)
InvScale(&a21, a11.Float())
//fmt.Printf("np == 1: cwrk\n%v\na21\n%v\n", &cwrk, &a21)
// store pivot point relative to original matrix
p1.pivots[0] = r + ATL.Rows() + 1
} else if np == 2 {
/*
* See comments for this block in unblkDecompBKLower().
*/
a := w11.GetAt(0, 0)
b := w11.GetAt(1, 0)
d := w11.GetAt(1, 1)
a11inv.SetAt(0, 0, d/b)
a11inv.SetAt(1, 1, a/b)
// denominator: (a/b)*(d/b)-1.0 == (a*d - b^2)/b^2
scale := 1.0 / ((a/b)*(d/b) - 1.0)
scale /= b
w11.SubMatrix(&cwrk, np, 0, a21.Rows(), np)
// a21 = a21*a11.-1
Mult(&a21, &cwrk, &a11inv, scale, 0.0, NOTRANS)
a11.SetAt(0, 0, a)
a11.SetAt(1, 0, b)
a11.SetAt(1, 1, d)
// store pivot point relative to original matrix
p1.pivots[0] = -(r + ATL.Rows() + 1)
p1.pivots[1] = p1.pivots[0]
}
/*
if m(&ABR) < 5 {
var Ablk, wblk, w5 matrix.FloatMatrix
merge1x2(&Ablk, &ABL, &ABR)
merge1x2(&wblk, &w10, &w11)
wblk.SubMatrix(&w5, 0, 0, Ablk.Rows(), wblk.Cols())
fmt.Printf("blocked EOL: Ablk r=%d, nc=%d. np=%d\n%v\n", r, nc, np, &Ablk)
fmt.Printf("wblk m(wblk)=%d:\n%v\n", m(&w5), &w5)
}
*/
// ------------------------------------------------------------
nc += np
continue3x3to2x2(
&ATL, &ATR,
&ABL, &ABR, &A00, &a11, &A22, A, pBOTTOMRIGHT)
contPivot3x1to2x1(
&pT,
&pB, &p0, &p1, p, pBOTTOM)
}
// undo applied partial row pivots (AL, w00)
//undo:
return err, nc
}示例10: findAndBuildBKPivotLower
/*
* Find diagonal pivot and build incrementaly updated block.
*
* (AL) (AR) (WL) (WR)
* -------------------------- ---------- k'th row in W
* x x | c1 w w | k kp1
* x x | c1 d w w | k kp1
* x x | c1 x d w w | k kp1
* x x | c1 x x d w w | k kp1
* x x | c1 r2 r2 r2 r2 w w | k kp1
* x x | c1 x x x r2 d w w | k kp1
* x x | c1 x x x r2 x d w w | k kp1
*
* Matrix AR contains the unfactored part of the matrix and AL the already
* factored columns. Matrix WL is updated values of factored part ie.
* w(i) = l(i)d(i). Matrix WR will have updated values for next column.
* Column WR(k) contains updated AR(c1) and WR(kp1) possible pivot row AR(r2).
*
*
*/
func findAndBuildBKPivotLower(AL, AR, WL, WR *matrix.FloatMatrix, k int) (int, int) {
var r, q int
var rcol, qrow, src, wk, wkp1, wrow matrix.FloatMatrix
// Copy AR column 0 to WR column 0 and update with WL[0:]
AR.SubMatrix(&src, 0, 0, AR.Rows(), 1)
WR.SubMatrix(&wk, 0, 0, AR.Rows(), 1)
src.CopyTo(&wk)
if k > 0 {
WL.SubMatrix(&wrow, 0, 0, 1, WL.Cols())
MVMult(&wk, AL, &wrow, -1.0, 1.0, NOTRANS)
//fmt.Printf("wk after update:\n%v\n", &wk)
}
if AR.Rows() == 1 {
return 0, 1
}
amax := math.Abs(WR.GetAt(0, 0))
// find max off-diagonal on first column.
WR.SubMatrix(&rcol, 1, 0, AR.Rows()-1, 1)
//fmt.Printf("rcol:\n%v\n", &rcol)
// r is row index and rmax is its absolute value
r = IAMax(&rcol) + 1
rmax := math.Abs(rcol.GetAt(r-1, 0))
//fmt.Printf("r=%d, amax=%e, rmax=%e\n", r, amax, rmax)
if amax >= bkALPHA*rmax {
// no pivoting, 1x1 diagonal
return 0, 1
}
// Now we need to copy row r to WR[:,1] and update it
WR.SubMatrix(&wkp1, 0, 1, AR.Rows(), 1)
AR.SubMatrix(&qrow, r, 0, 1, r+1)
qrow.CopyTo(&wkp1)
//fmt.Printf("m(AR)=%d, r=%d, qrow: %v\n", AR.Rows(), r, &qrow)
if r < AR.Rows()-1 {
var wkr matrix.FloatMatrix
AR.SubMatrix(&qrow, r, r, AR.Rows()-r, 1)
wkp1.SubMatrix(&wkr, r, 0, wkp1.Rows()-r, 1)
qrow.CopyTo(&wkr)
//fmt.Printf("m(AR)=%d, r=%d, qrow: %v\n", AR.Rows(), r, &qrow)
}
if k > 0 {
// update wkp1
WL.SubMatrix(&wrow, r, 0, 1, WL.Cols())
//fmt.Printf("initial wpk1:\n%v\n", &wkp1)
MVMult(&wkp1, AL, &wrow, -1.0, 1.0, NOTRANS)
//fmt.Printf("updated wpk1:\n%v\n", &wkp1)
}
// set on-diagonal entry to zero to avoid finding it
p1 := wkp1.GetAt(r, 0)
wkp1.SetAt(r, 0, 0.0)
// max off-diagonal on r'th column/row at index q
q = IAMax(&wkp1)
qmax := math.Abs(wkp1.GetAt(q, 0))
// restore on-diagonal entry
wkp1.SetAt(r, 0, p1)
//arr := math.Abs(WR.GetAt(r, 1))
//fmt.Printf("blk: r=%d, q=%d, amax=%e, rmax=%e, qmax=%e, Arr=%e\n", r, q, amax, rmax, qmax, arr)
if amax >= bkALPHA*rmax*(rmax/qmax) {
// no pivoting, 1x1 diagonal
return 0, 1
}
// if q == r then qmax is not off-diagonal, qmax == WR[r,1] and
// we get 1x1 pivot as following is always true
if math.Abs(WR.GetAt(r, 1)) >= bkALPHA*qmax {
// 1x1 pivoting and interchange with k, r
// pivot row in column WR[:,1] to W[:,0]
//pr := WR.GetAt(r, 1)
//_ = pr
WR.SubMatrix(&src, 0, 1, AR.Rows(), 1)
WR.SubMatrix(&wkp1, 0, 0, AR.Rows(), 1)
src.CopyTo(&wkp1)
wkp1.SetAt(0, 0, src.GetAt(r, 0))
wkp1.SetAt(r, 0, src.GetAt(0, 0))
return r, 1
} else {
// 2x2 pivoting and interchange with k+1, r
return r, 2
//.........這裏部分代碼省略.........
示例11: unblkDecompBKUpper
/*
* Unblocked Bunch-Kauffman LDL factorization.
*
* Corresponds lapack.DSYTF2
*/
func unblkDecompBKUpper(A, wrk *matrix.FloatMatrix, p *pPivots) (error, int) {
var err error
var ATL, ATR, ABL, ABR matrix.FloatMatrix
var A00, a01, A02, a12t, a11, A22, a11inv matrix.FloatMatrix
var cwrk matrix.FloatMatrix
var pT, pB, p0, p1, p2 pPivots
err = nil
nc := 0
partition2x2(
&ATL, &ATR,
&ABL, &ABR, A, 0, 0, pBOTTOMRIGHT)
partitionPivot2x1(
&pT,
&pB, p, 0, pBOTTOM)
// permanent working space for symmetric inverse of a11
wrk.SubMatrix(&a11inv, 0, wrk.Cols()-2, 2, 2)
a11inv.SetAt(1, 0, -1.0)
a11inv.SetAt(0, 1, -1.0)
for ATL.Cols() > 0 {
nr := ATL.Rows() - 1
r, np := findBKPivot(&ATL, UPPER)
if r != -1 /*&& r != np-1*/ {
// pivoting needed; do swaping here
//fmt.Printf("pre-pivot ATL [%d]:\n%v\n", ATL.Rows()-np, &ATL)
applyBKPivotSym(&ATL, ATL.Rows()-np, r, UPPER)
if np == 2 {
/*
* [r,r] | [r, nr]
* a11 == --------------- 2-by-2 pivot, swapping [nr-1,nr] and [r,nr]
* [r,0] | [nr,nr]
*/
t := ATL.GetAt(nr-1, nr)
ATL.SetAt(nr-1, nr, ATL.GetAt(r, nr))
ATL.SetAt(r, nr, t)
}
//fmt.Printf("unblk: ATL after %d pivot [r=%d]:\n%v\n", np, r, &ATL)
}
// repartition according the pivot size
repartition2x2to3x3(&ATL,
&A00, &a01, &A02,
nil, &a11, &a12t,
nil, nil, &A22 /**/, A, np, pTOPLEFT)
repartPivot2x1to3x1(&pT,
&p0,
&p1,
&p2 /**/, p, np, pTOP)
// ------------------------------------------------------------
if np == 1 {
// A00 = A00 - a01*a01.T/a11
MVUpdateTrm(&A00, &a01, &a01, -1.0/a11.Float(), UPPER)
// a01 = a01/a11
InvScale(&a01, a11.Float())
if r == -1 {
p1.pivots[0] = ATL.Rows()
} else {
p1.pivots[0] = r + 1
}
} else if np == 2 {
/*
* See comments on unblkDecompBKLower().
*/
a := a11.GetAt(0, 0)
b := a11.GetAt(0, 1)
d := a11.GetAt(1, 1)
a11inv.SetAt(0, 0, d/b)
a11inv.SetAt(1, 1, a/b)
// denominator: (a/b)*(d/b)-1.0 == (a*d - b^2)/b^2
scale := 1.0 / ((a/b)*(d/b) - 1.0)
scale /= b
// cwrk = a21
wrk.SubMatrix(&cwrk, 2, 0, a01.Rows(), a01.Cols())
a01.CopyTo(&cwrk)
//fmt.Printf("cwrk:\n%v\n", &cwrk)
//fmt.Printf("a11inv:\n%v\n", &a11inv)
// a01 = a01*a11.-1
Mult(&a01, &cwrk, &a11inv, scale, 0.0, NOTRANS)
// A00 = A00 - a01*a11.-1*a01.T = A00 - a01*cwrk.T
UpdateTrm(&A00, &a01, &cwrk, -1.0, 1.0, UPPER|TRANSB)
p1.pivots[0] = -(r + 1)
p1.pivots[1] = p1.pivots[0]
}
// ------------------------------------------------------------
nc += np
continue3x3to2x2(
&ATL, &ATR,
//.........這裏部分代碼省略.........
示例12: unblkDecompBKLower
/*
* Unblocked Bunch-Kauffman LDL factorization.
*
* Corresponds lapack.DSYTF2
*/
func unblkDecompBKLower(A, wrk *matrix.FloatMatrix, p *pPivots) (error, int) {
var err error
var ATL, ATR, ABL, ABR matrix.FloatMatrix
var A00, a10t, a11, A20, a21, A22, a11inv matrix.FloatMatrix
var cwrk matrix.FloatMatrix
var pT, pB, p0, p1, p2 pPivots
err = nil
nc := 0
partition2x2(
&ATL, &ATR,
&ABL, &ABR, A, 0, 0, pTOPLEFT)
partitionPivot2x1(
&pT,
&pB, p, 0, pTOP)
// permanent working space for symmetric inverse of a11
wrk.SubMatrix(&a11inv, 0, wrk.Cols()-2, 2, 2)
a11inv.SetAt(1, 0, -1.0)
a11inv.SetAt(0, 1, -1.0)
for ABR.Cols() > 0 {
r, np := findBKPivot(&ABR, LOWER)
if r != 0 && r != np-1 {
// pivoting needed; do swaping here
applyBKPivotSym(&ABR, np-1, r, LOWER)
if np == 2 {
/*
* [0,0] | [r,0]
* a11 == ------------- 2-by-2 pivot, swapping [1,0] and [r,0]
* [r,0] | [r,r]
*/
t := ABR.GetAt(1, 0)
ABR.SetAt(1, 0, ABR.GetAt(r, 0))
ABR.SetAt(r, 0, t)
}
//fmt.Printf("unblk: ABR after %d pivot [r=%d]:\n%v\n", np, r, &ABR)
}
// repartition according the pivot size
repartition2x2to3x3(&ATL,
&A00, nil, nil,
&a10t, &a11, nil,
&A20, &a21, &A22 /**/, A, np, pBOTTOMRIGHT)
repartPivot2x1to3x1(&pT,
&p0,
&p1,
&p2 /**/, p, np, pBOTTOM)
// ------------------------------------------------------------
if np == 1 {
// A22 = A22 - a21*a21.T/a11
MVUpdateTrm(&A22, &a21, &a21, -1.0/a11.Float(), LOWER)
// a21 = a21/a11
InvScale(&a21, a11.Float())
// store pivot point relative to original matrix
p1.pivots[0] = r + ATL.Rows() + 1
} else if np == 2 {
/* from Bunch-Kaufmann 1977:
* (E2 C.T) = ( I2 0 )( E 0 )( I[n-2] E.-1*C.T )
* (C B ) ( C*E.-1 I[n-2] )( 0 A[n-2] )( 0 I2 )
*
* A[n-2] = B - C*E.-1*C.T
*
* E.-1 is inverse of a symmetric matrix, cannot use
* triangular solve. We calculate inverse of 2x2 matrix.
* Following is inspired by lapack.SYTF2
*
* a | b 1 d | -b b d/b | -1
* inv ----- = ------ * ------ = ----------- * --------
* b | d (ad-b^2) -b | a (a*d - b^2) -1 | a/b
*
*/
a := a11.GetAt(0, 0)
b := a11.GetAt(1, 0)
d := a11.GetAt(1, 1)
a11inv.SetAt(0, 0, d/b)
a11inv.SetAt(1, 1, a/b)
// denominator: (a/b)*(d/b)-1.0 == (a*d - b^2)/b^2
scale := 1.0 / ((a/b)*(d/b) - 1.0)
scale /= b
// cwrk = a21
wrk.SubMatrix(&cwrk, 2, 0, a21.Rows(), a21.Cols())
a21.CopyTo(&cwrk)
// a21 = a21*a11.-1
Mult(&a21, &cwrk, &a11inv, scale, 0.0, NOTRANS)
// A22 = A22 - a21*a11.-1*a21.T = A22 - a21*cwrk.T
UpdateTrm(&A22, &a21, &cwrk, -1.0, 1.0, LOWER|TRANSB)
// store pivot point relative to original matrix
p1.pivots[0] = -(r + ATL.Rows() + 1)
p1.pivots[1] = p1.pivots[0]
//.........這裏部分代碼省略.........
示例13: blockedQR
/*
* Blocked QR decomposition with compact WY transform. As implemented
* in lapack.xGEQRF subroutine.
*/
func blockedQR(A, Tvec, Twork, W *matrix.FloatMatrix, nb int) {
var ATL, ATR, ABL, ABR, AL, AR matrix.FloatMatrix
var A00, A01, A02, A10, A11, A12, A20, A21, A22 matrix.FloatMatrix
var TT, TB matrix.FloatMatrix
var t0, tau, t2, Tdiag, WT, WB, W0, W1, W2 matrix.FloatMatrix
//var Twork, W *matrix.FloatMatrix
Tdiag.DiagOf(Twork)
partition2x2(
&ATL, &ATR,
&ABL, &ABR, A, 0, 0, pTOPLEFT)
partition2x1(
&TT,
&TB, Tvec, 0, pTOP)
partition2x1(
&WT,
&WB, W, 0, pTOP)
for ABR.Rows() > 0 && ABR.Cols() > 0 {
repartition2x2to3x3(&ATL,
&A00, &A01, &A02,
&A10, &A11, &A12,
&A20, &A21, &A22, A, nb, pBOTTOMRIGHT)
repartition2x1to3x1(&TT,
&t0,
&tau,
&t2, Tvec, nb, pBOTTOM)
repartition2x1to3x1(&WT,
&W0,
&W1,
&W2, W, nb, pBOTTOM)
partition1x2(
&AL, &AR, &ABR, nb, pLEFT)
// current block size
cb, rb := A11.Size()
if rb < cb {
cb = rb
}
// --------------------------------------------------------
// decompose left side AL == /A11\
// \A21/
unblockedQRT(&AL, Twork)
// copy scaling from T diagonal to tau-vector
Tdiag.CopyTo(&tau)
// update A'tail i.e. A12 and A22 with (I - Y*T*Y.T).T * A'tail
// compute: C - Y*(C.T*Y*T).T
updateWithQT(&A12, &A22, &A11, &A21, Twork, &W2, cb, true)
// --------------------------------------------------------
continue3x3to2x2(
&ATL, &ATR,
&ABL, &ABR, &A00, &A11, &A22, A, pBOTTOMRIGHT)
continue3x1to2x1(
&TT,
&TB, &t0, &tau, Tvec, pBOTTOM)
continue3x1to2x1(
&WT,
&WB, &W0, &W1, W, pBOTTOM)
}
}示例14: unblkBoundedLowerLDL
func unblkBoundedLowerLDL(A, W *matrix.FloatMatrix, p *pPivots, ncol int) (error, int) {
var ATL, ATR, ABL, ABR matrix.FloatMatrix
var A00, a10, a11, A20, a21, A22, adiag, wcol matrix.FloatMatrix
var w00, w10, w11 matrix.FloatMatrix
var pT, pB, p0, p1, p2 pPivots
var err error = nil
partition2x2(
&ATL, &ATR,
&ABL, &ABR, A, 0, 0, pTOPLEFT)
partitionPivot2x1(
&pT,
&pB, p, 0, pTOP)
// copy current diagonal to last column of workspace
W.SubMatrix(&wcol, 0, W.Cols()-1, A.Rows(), 1)
A.Diag(&adiag)
adiag.CopyTo(&wcol)
//fmt.Printf("initial diagonal:\n%v\n", &wcol)
nc := 0
for ABR.Cols() > 0 && nc < ncol {
partition2x2(
&w00, nil,
&w10, &w11, W, nc, nc, pTOPLEFT)
dmax := findAndBuildPivot(&ABL, &ABR, &w10, &w11, nc)
//fmt.Printf("dmax=%d\n", dmax)
if dmax > 0 {
// pivot diagonal in symmetric matrix; will swap a11 [0,0] and [imax,imax]
applyPivotSym(&ABL, &ABR, dmax, LOWER)
swapRows(&w10, 0, dmax)
pB.pivots[0] = dmax + ATL.Rows() + 1
} else {
pB.pivots[0] = 0
}
//fmt.Printf("blk pivoted %d, A:\n%v\nW:\n%v\n", dmax, A, W)
repartition2x2to3x3(&ATL,
&A00, nil, nil,
&a10, &a11, nil,
&A20, &a21, &A22, A, 1, pBOTTOMRIGHT)
repartPivot2x1to3x1(&pT,
&p0, &p1, &p2 /**/, p, 1, pBOTTOM)
// --------------------------------------------------------
// Copy updated column from working space
w11.SubMatrix(&wcol, 1, 0, a21.Rows(), 1)
wcol.CopyTo(&a21)
a11.SetAt(0, 0, w11.GetAt(0, 0))
// l21 = a21/a11
InvScale(&a21, a11.Float())
// here: wcol == l21*d11 == a21
if ncol-nc > 1 {
// update diagonal in workspace if not last column of block
w11.SubMatrix(&adiag, 1, w11.Cols()-1, a21.Rows(), 1)
MVUpdateDiag(&adiag, &wcol, &wcol, -1.0/a11.Float())
}
//fmt.Printf("nc=%d, a11=%e\n", nc, a11.Float())
//fmt.Printf("l21\n%v\n", &a21)
//fmt.Printf("a21\n%v\n", &wcol)
//fmt.Printf("diag\n%v\n", &adiag)
//var Ablk, wblk matrix.FloatMatrix
//merge1x2(&Ablk, &ABL, &ABR)
//merge1x2(&wblk, &w10, &w11)
//fmt.Printf("unblk Ablk:\n%v\n", &Ablk)
//fmt.Printf("unblk wblk:\n%v\n", &wblk)
// ---------------------------------------------------------
nc++
continue3x3to2x2(
&ATL, &ATR,
&ABL, &ABR, &A00, &a11, &A22, A, pBOTTOMRIGHT)
contPivot3x1to2x1(
&pT,
&pB, &p0, &p1, p, pBOTTOM)
}
return err, nc
}