本文整理汇总了Golang中github.com/hrautila/matrix.FloatMatrix.Float方法的典型用法代码示例。如果您正苦于以下问题:Golang FloatMatrix.Float方法的具体用法?Golang FloatMatrix.Float怎么用?Golang FloatMatrix.Float使用的例子?那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类github.com/hrautila/matrix.FloatMatrix
的用法示例。
在下文中一共展示了FloatMatrix.Float方法的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的Golang代码示例。
示例1: unblockedLUnoPiv
// unblocked LU decomposition w/o pivots, FLAME LU nopivots variant 5
func unblockedLUnoPiv(A *matrix.FloatMatrix) (err error) {
var ATL, ATR, ABL, ABR matrix.FloatMatrix
var A00, a01, A02, a10, a11, a12, A20, a21, A22 matrix.FloatMatrix
err = nil
partition2x2(
&ATL, &ATR,
&ABL, &ABR, A, 0, 0, pTOPLEFT)
for ATL.Rows() < A.Rows() {
repartition2x2to3x3(&ATL,
&A00, &a01, &A02,
&a10, &a11, &a12,
&A20, &a21, &A22, A, 1, pBOTTOMRIGHT)
// a21 = a21/a11
//a21.Scale(1.0/a11.Float())
InvScale(&a21, a11.Float())
// A22 = A22 - a21*a12
err = MVRankUpdate(&A22, &a21, &a12, -1.0)
continue3x3to2x2(
&ATL, &ATR,
&ABL, &ABR, &A00, &a11, &A22, A, pBOTTOMRIGHT)
}
return
}
示例2: unblockedInverseLower
// Inverse NON-UNIT diagonal tridiagonal matrix
func unblockedInverseLower(A *matrix.FloatMatrix) (err error) {
var ATL, ATR, ABL, ABR matrix.FloatMatrix
var A00, a10t, a11, A20, a21, A22 matrix.FloatMatrix
err = nil
partition2x2(
&ATL, &ATR,
&ABL, &ABR, A, 0, 0, pTOPLEFT)
for ATL.Rows() < A.Rows() {
repartition2x2to3x3(&ATL,
&A00, nil, nil,
&a10t, &a11, nil,
&A20, &a21, &A22, A, 1, pBOTTOMRIGHT)
// -------------------------------------------------
aval := a11.Float()
// a21 = -a21/a11
InvScale(&a21, -aval)
// A20 = A20 + a21*a10.t
MVRankUpdate(&A20, &a21, &a10t, 1.0)
// a10 = a10/a11
InvScale(&a10t, aval)
// a11 = 1.0/a11
a11.SetAt(0, 0, 1.0/aval)
// -------------------------------------------------
continue3x3to2x2(
&ATL, &ATR,
&ABL, &ABR, &A00, &a11, &A22, A, pBOTTOMRIGHT)
}
return
}
示例3: unblockedInverseUpper
// Inverse NON-UNIT diagonal tridiagonal matrix
func unblockedInverseUpper(A *matrix.FloatMatrix) (err error) {
var ATL, ATR, ABL, ABR matrix.FloatMatrix
var A00, a01, A02, a11, a12t, A22 matrix.FloatMatrix
err = nil
partition2x2(
&ATL, &ATR,
&ABL, &ABR, A, 0, 0, pTOPLEFT)
for ATL.Rows() < A.Rows() {
repartition2x2to3x3(&ATL,
&A00, &a01, &A02,
nil, &a11, &a12t,
nil, nil, &A22, A, 1, pBOTTOMRIGHT)
// -------------------------------------------------
aval := a11.Float()
// a12 = -a12/a11
InvScale(&a12t, -aval)
// A02 = A02 + a01*a12
MVRankUpdate(&A02, &a01, &a12t, 1.0)
// a01 = a01/a11
InvScale(&a01, aval)
// a11 = 1.0/a11
a11.SetAt(0, 0, 1.0/aval)
// -------------------------------------------------
continue3x3to2x2(
&ATL, &ATR,
&ABL, &ABR, &A00, &a11, &A22, A, pBOTTOMRIGHT)
}
return
}
示例4: _TestViewUpdate
func _TestViewUpdate(t *testing.T) {
Adata2 := [][]float64{
[]float64{4.0, 2.0, 2.0},
[]float64{6.0, 4.0, 2.0},
[]float64{4.0, 6.0, 1.0},
}
A := matrix.FloatMatrixFromTable(Adata2, matrix.RowOrder)
N := A.Rows()
// simple LU decomposition without pivoting
var A11, a10, a01, a00 matrix.FloatMatrix
for k := 1; k < N; k++ {
a00.SubMatrixOf(A, k-1, k-1, 1, 1)
a01.SubMatrixOf(A, k-1, k, 1, A.Cols()-k)
a10.SubMatrixOf(A, k, k-1, A.Rows()-k, 1)
A11.SubMatrixOf(A, k, k)
//t.Logf("A11: %v a01: %v\n", A11, a01)
a10.Scale(1.0 / a00.Float())
MVRankUpdate(&A11, &a10, &a01, -1.0)
}
Ld := TriLU(A.Copy())
Ud := TriU(A)
t.Logf("Ld:\n%v\nUd:\n%v\n", Ld, Ud)
An := matrix.FloatZeros(N, N)
Mult(An, Ld, Ud, 1.0, 1.0, NOTRANS)
t.Logf("A == Ld*Ud: %v\n", An.AllClose(An))
}
示例5: unblkUpperLDLnoPiv
/*
* ( A11 a12 ) ( U11 u12 )( D1 0 )( U11.t 0 )
* ( a21 a22 ) ( 0 1 )( 0 d2 )( u12.t 1 )
*
* a22 = d2
* a01 = u12*d2 => u12 = a12/d2
* A11 = u12*d2*u12.t + U11*D1*U11.t => U11 = A11 - u12*d2*u12.t
*/
func unblkUpperLDLnoPiv(A *matrix.FloatMatrix) (err error) {
var ATL, ATR, ABL, ABR matrix.FloatMatrix
var A00, a01, A02, a11, a12, A22 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, 1, pTOPLEFT)
// --------------------------------------------------------
// A00 = A00 - u01*d11*u01.T = A00 - a01*a01.T/a11; triangular update
err = MVUpdateTrm(&A00, &a01, &a01, -1.0/a11.Float(), UPPER)
// u01 = a01/a11
InvScale(&a01, a11.Float())
// ---------------------------------------------------------
continue3x3to2x2(
&ATL, &ATR,
&ABL, &ABR, &A00, &a11, &A22, A, pTOPLEFT)
}
return
}
示例6: unblkLowerLDLnoPiv
/*
* ( a11 a12 ) ( 1 0 )( d1 0 )( l l21.t )
* ( a21 A22 ) ( l21 L22 )( 0 A22 )( 0 L22.t )
*
* a11 = d1
* a21 = l21*d1 => l21 = a21/d1
* A22 = l21*d1*l21.t + L22*D2*L22.t => L22 = A22 - l21*d1*l21t
*/
func unblkLowerLDLnoPiv(A *matrix.FloatMatrix) (err error) {
var ATL, ATR, ABL, ABR matrix.FloatMatrix
var A00, a10, a11, A20, a21, A22 matrix.FloatMatrix
err = nil
partition2x2(
&ATL, &ATR,
&ABL, &ABR, A, 0, 0, pTOPLEFT)
for ATL.Rows() < A.Rows() {
repartition2x2to3x3(&ATL,
&A00, nil, nil,
&a10, &a11, nil,
&A20, &a21, &A22, A, 1, pBOTTOMRIGHT)
// --------------------------------------------------------
// d11 = a11; no-op
// A22 = A22 - l21*d11*l21.T = A22 - a21*a21.T/a11; triangular update
err = MVUpdateTrm(&A22, &a21, &a21, -1.0/a11.Float(), LOWER)
// l21 = a21/a11
InvScale(&a21, a11.Float())
// ---------------------------------------------------------
continue3x3to2x2(
&ATL, &ATR,
&ABL, &ABR, &A00, &a11, &A22, A, pBOTTOMRIGHT)
}
return
}
示例7: _TestPartition1Dh
func _TestPartition1Dh(t *testing.T) {
var AL, AR, A0, a1, A2 matrix.FloatMatrix
A := matrix.FloatZeros(1, 6)
partition1x2(&AL, &AR, A, 0, pRIGHT)
t.Logf("m(AL)=%d, m(AR)=%d\n", AL.Cols(), AR.Cols())
for AL.Cols() < A.Cols() {
AR.Add(1.0)
t.Logf("m(AR)=%d; %v\n", AR.Cols(), AR)
repartition1x2to1x3(&AL, &A0, &a1, &A2, A, 1, pRIGHT)
t.Logf("m(A0)=%d, m(A2)=%d, a1=%.1f\n", A0.Cols(), A2.Cols(), a1.Float())
continue1x3to1x2(&AL, &AR, &A0, &a1, A, pRIGHT)
}
}
示例8: unblkUpperLDL
/*
* ( A11 a12 ) ( U11 u12 )( D1 0 )( U11.t 0 )
* ( a21 a22 ) ( 0 1 )( 0 d2 )( u12.t 1 )
*
* a22 = d2
* a01 = u12*d2 => u12 = a12/d2
* A11 = u12*d2*u12.t + U11*D1*U11.t => U11 = A11 - u12*d2*u12.t
*/
func unblkUpperLDL(A *matrix.FloatMatrix, p *pPivots) (err error) {
var ATL, ATR, ABL, ABR matrix.FloatMatrix
var A00, a01, A02, a11, a12, A22 matrix.FloatMatrix
var AL, AR, acol 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, 1, pTOPLEFT)
repartPivot2x1to3x1(&pT,
&p0, &p1, &p2 /**/, p, 1, pTOP)
// --------------------------------------------------------
// search diagonal; diag(A00;a11)
ATL.Diag(&acol)
//merge2x1(&acol, &a01, &a11)
imax := IAMax(&acol)
if imax < ATL.Rows()-1 {
merge1x2(&AL, &ATL, &ATR)
merge1x2(&AR, &a11, &a12)
// pivot diagonal in symmetric matrix; will swap a11 and [imax,imax]
applyPivotSym(&AL, &AR, imax, UPPER)
p1.pivots[0] = imax + 1
} else {
p1.pivots[0] = 0
}
if a11.Float() == 0.0 {
err = onError("zero on diagonal.")
return
}
// A00 = A00 - u01*d11*u01.T = A00 - a01*a01.T/a11; triangular update
err = MVUpdateTrm(&A00, &a01, &a01, -1.0/a11.Float(), UPPER)
// u01 = a01/a11
InvScale(&a01, a11.Float())
// ---------------------------------------------------------
continue3x3to2x2(
&ATL, &ATR,
&ABL, &ABR, &A00, &a11, &A22, A, pTOPLEFT)
contPivot3x1to2x1(
&pT,
&pB, &p0, &p1, p, pTOP)
}
return
}
示例9: unblockedCHOL
func unblockedCHOL(A *matrix.FloatMatrix, flags Flags, nr int) (err error) {
var ATL, ATR, ABL, ABR matrix.FloatMatrix
var A00, a01, A02, a10, a11, a12, A20, a21, A22 matrix.FloatMatrix
err = nil
partition2x2(
&ATL, &ATR,
&ABL, &ABR, A, 0, 0, pTOPLEFT)
for ATL.Rows() < A.Rows() {
repartition2x2to3x3(&ATL,
&A00, &a01, &A02,
&a10, &a11, &a12,
&A20, &a21, &A22, A, 1, pBOTTOMRIGHT)
// a11 = sqrt(a11)
aval := math.Sqrt(a11.Float())
if math.IsNaN(aval) {
panic(fmt.Sprintf("illegal value at %d: %e", nr+ATL.Rows(), a11.Float()))
}
a11.SetAt(0, 0, aval)
if flags&LOWER != 0 {
// a21 = a21/a11
InvScale(&a21, a11.Float())
// A22 = A22 - a21*a21' (SYR)
err = MVRankUpdateSym(&A22, &a21, -1.0, flags)
} else {
// a21 = a12/a11
InvScale(&a12, a11.Float())
// A22 = A22 - a12'*a12 (SYR)
err = MVRankUpdateSym(&A22, &a12, -1.0, flags)
}
continue3x3to2x2(
&ATL, &ATR,
&ABL, &ABR, &A00, &a11, &A22, A, pBOTTOMRIGHT)
}
return
}
示例10: 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
//.........这里部分代码省略.........
示例11: 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
}
示例12: 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,
//.........这里部分代码省略.........
示例13: unblkBoundedBKLower
/*
* Unblocked, bounded Bunch-Kauffman LDL factorization for at most ncol columns.
* At most ncol columns are factorized and trailing matrix updates are restricted
* to ncol columns. Also original columns are accumulated to working matrix, which
* is used by calling blocked algorithm to update the trailing matrix with BLAS3
* update.
*
* Corresponds lapack.DLASYF
*/
func unblkBoundedBKLower(A, wrk *matrix.FloatMatrix, p *pPivots, ncol int) (error, int) {
var err error
var ATL, ATR, ABL, ABR matrix.FloatMatrix
var A00, a10t, a11, A20, a21, A22, a11inv matrix.FloatMatrix
var w00, w10, w11 matrix.FloatMatrix
var cwrk matrix.FloatMatrix
//var s, d 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, 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 && nc < ncol {
partition2x2(
&w00, nil,
&w10, &w11, wrk, nc, nc, pTOPLEFT)
//fmt.Printf("ABR:\n%v\n", &ABR)
r, np := findAndBuildBKPivotLower(&ABL, &ABR, &w10, &w11, nc)
//fmt.Printf("after find: r=%d, np=%d, ncol=%d, nc=%d\n", r, np, ncol, nc)
if np > ncol-nc {
// next pivot does not fit into ncol columns, restore last column,
// return with number of factorized columns
//fmt.Printf("np > ncol-nc: %d > %d\n", np, ncol-nc)
return err, nc
//goto undo
}
if r != 0 && r != np-1 {
// pivoting needed; do swaping here
applyBKPivotSym(&ABR, np-1, r, LOWER)
// swap left hand rows to get correct updates
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))
//.........这里部分代码省略.........
示例14: 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]
//.........这里部分代码省略.........
示例15: unblkSolveBKUpper
func unblkSolveBKUpper(B, A *matrix.FloatMatrix, p *pPivots, phase int) error {
var err error
var ATL, ATR, ABL, ABR matrix.FloatMatrix
var A00, a01, A02, a11, a12t, A22 matrix.FloatMatrix
var Aref *matrix.FloatMatrix
var BT, BB, B0, b1, B2, Bx matrix.FloatMatrix
var pT, pB, p0, p1, p2 pPivots
var aStart, aDir, bStart, bDir pDirection
var nc int
err = nil
np := 0
if phase == 2 {
aStart = pTOPLEFT
aDir = pBOTTOMRIGHT
bStart = pTOP
bDir = pBOTTOM
nc = 1
Aref = &ABR
} else {
aStart = pBOTTOMRIGHT
aDir = pTOPLEFT
bStart = pBOTTOM
bDir = pTOP
nc = A.Rows()
Aref = &ATL
}
partition2x2(
&ATL, &ATR,
&ABL, &ABR, A, 0, 0, aStart)
partition2x1(
&BT,
&BB, B, 0, bStart)
partitionPivot2x1(
&pT,
&pB, p, 0, bStart)
// ABR.Cols() == 0 is end of matrix,
for Aref.Cols() > 0 {
// see if next diagonal block is 1x1 or 2x2
np = 1
if p.pivots[nc-1] < 0 {
np = 2
}
fmt.Printf("nc=%d, np=%d, m(ABR)=%d\n", nc, np, m(&ABR))
// repartition according the pivot size
repartition2x2to3x3(&ATL,
&A00, &a01, &A02,
nil, &a11, &a12t,
nil, nil, &A22 /**/, A, np, aDir)
repartition2x1to3x1(&BT,
&B0,
&b1,
&B2 /**/, B, np, bDir)
repartPivot2x1to3x1(&pT,
&p0,
&p1,
&p2 /**/, p, np, bDir)
// ------------------------------------------------------------
switch phase {
case 1:
// computes D.-1*(L.-1*B)
if np == 1 {
if p1.pivots[0] != nc {
// swap rows in top part of B
//fmt.Printf("1x1 pivot top with %d [%d]\n", p1.pivots[0], p1.pivots[0]-BT.Rows())
swapRows(&BT, BT.Rows()-1, p1.pivots[0]-1)
}
// B2 = B2 - a21*b1
MVRankUpdate(&B2, &a01, &b1, -1.0)
// b1 = b1/d1
InvScale(&b1, a11.Float())
nc += 1
} else if np == 2 {
if p1.pivots[0] != -nc {
// swap rows on bottom part of B
//fmt.Printf("2x2 pivot %d with %d [%d]\n", nc+1, -p1.pivots[0])
//fmt.Printf("pre :\n%v\n", B)
swapRows(&BT, BT.Rows()-2, -p1.pivots[0]-1)
//fmt.Printf("post:\n%v\n", B)
}
b := a11.GetAt(0, 1)
apb := a11.GetAt(0, 0) / b
dpb := a11.GetAt(1, 1) / b
// (a/b)*(d/b)-1.0 == (a*d - b^2)/b^2
scale := apb*dpb - 1.0
scale *= b
// B2 = B2 - a21*b1
Mult(&B2, &a01, &b1, -1.0, 1.0, NOTRANS)
// b1 = a11.-1*b1.T
//(2x2 block, no subroutine for doing this in-place)
for k := 0; k < b1.Cols(); k++ {
s0 := b1.GetAt(0, k)
s1 := b1.GetAt(1, k)
b1.SetAt(0, k, (dpb*s0-s1)/scale)
//.........这里部分代码省略.........