Golang FloatMatrix.GetAt方法代碼示例

/* From LAPACK/dlarfg.f
 *
 * DLARFG generates a real elementary reflector H of order n, such
 * that
 *
 *       H * ( alpha ) = ( beta ),   H**T * H = I.
 *           (   x   )   (   0  )
 *
 * where alpha and beta are scalars, and x is an (n-1)-element real
 * vector. 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 (n-1)-element
 * vector.
 *
 * If the elements of x are all zero, then tau = 0 and H is taken to be
 * the unit matrix.
 *
 * Otherwise  1 <= tau <= 2.
 */
func computeHouseholder(a11, x, tau *matrix.FloatMatrix, flags Flags) {

	// norm_x2 = ||x||_2
	norm_x2 := Norm2(x)
	if norm_x2 == 0.0 {
		//a11.SetAt(0, 0, -a11.GetAt(0, 0))
		tau.SetAt(0, 0, 0.0)
		return
	}

	alpha := a11.GetAt(0, 0)
	sign := 1.0
	if math.Signbit(alpha) {
		sign = -1.0
	}
	// beta = -(alpha / |alpha|) * ||alpha x||
	//      = -sign(alpha) * sqrt(alpha**2, norm_x2**2)
	beta := -sign * sqrtX2Y2(alpha, norm_x2)

	// x = x /(a11 - beta)
	InvScale(x, alpha-beta)

	tau.SetAt(0, 0, (beta-alpha)/beta)
	a11.SetAt(0, 0, beta)
}

/*
 * 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)
	}
}

/*
 * 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)
	}
}

// Find largest absolute value on column
func pivotIndex(A *matrix.FloatMatrix, p *pPivots) {
	max := math.Abs(A.GetAt(0, 0))
	for k := 1; k < A.Rows(); k++ {
		v := math.Abs(A.GetAt(k, 0))
		if v > max {
			p.pivots[0] = k
			max = v
		}
	}
}

/*
 * 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)
	}
}

func applyHHTo1x1(tau, v, A2, w1 *matrix.FloatMatrix, flags Flags) {

	tval := tau.GetAt(0, 0)
	if tval == 0.0 {
		return
	}
	if flags&LEFT != 0 {
		// w1 = A2.T*v
		MVMult(w1, A2, v, 1.0, 0.0, TRANSA)
	} else {
		// w1 = A2*v
		MVMult(w1, A2, v, 1.0, 0.0, NOTRANS)
	}

	// A2 = A2 - tau*v*w1
	MVRankUpdate(A2, v, w1, -tval)
}

/*
 * Apply diagonal pivot (row and column swapped) to symmetric matrix blocks.
 * AR[0,0] is on diagonal and AL is block to the left of diagonal and AR the
 * triangular diagonal block. Need to swap row and column.
 *
 * LOWER triangular; moving from top-left to bottom-right
 *
 *    d
 *    x  d
 *    x  x  d  |
 *    --------------------------
 *    S1 S1 S1 | P1 x  x  x  P2     -- current row
 *    x  x  x  | S2 d  x  x  x
 *    x  x  x  | S2 x  d  x  x
 *    x  x  x  | S2 x  x  d  x
 *    D1 D1 D1 | P2 D2 D2 D2 P3     -- swap with row 'index'
 *    x  x  x  | S3 x  x  x  D3 d
 *    x  x  x  | S3 x  x  x  D3 x d
 *       (ABL)          (ABR)
 *
 * UPPER triangular; moving from bottom-right to top-left
 *
 *         (ATL)             (ATR)
 *    d  x  x  D3 x  x  x | S3 x  x
 *       d  x  D3 x  x  x | S3 x  x
 *          d  D3 x  x  x | S3 x  x
 *             P3 D2 D2 D2| P2 D1 D1
 *                d  x  x | S2 x  x
 *                   d  x | S2 x  x
 *                      d | S2 x  x
 *    -----------------------------
 *                        | P1 S1 S1
 *                        |    d  x
 *                        |       d
 *                           (ABR)
 */
func applyPivotSym(AL, AR *matrix.FloatMatrix, index int, flags Flags) {
	var s, d matrix.FloatMatrix
	if flags&LOWER != 0 {
		// AL is [ABL]; AR is [ABR]; P1 is AR[0,0], P2 is AR[index, 0]
		// S1 -- D1
		AL.SubMatrix(&s, 0, 0, 1, AL.Cols())
		AL.SubMatrix(&d, index, 0, 1, AL.Cols())
		Swap(&s, &d)
		// S2 -- D2
		AR.SubMatrix(&s, 1, 0, index-1, 1)
		AR.SubMatrix(&d, index, 1, 1, index-1)
		Swap(&s, &d)
		// S3 -- D3
		AR.SubMatrix(&s, index+1, 0, AR.Rows()-index-1, 1)
		AR.SubMatrix(&d, index+1, index, AR.Rows()-index-1, 1)
		Swap(&s, &d)
		// swap P1 and P3
		p1 := AR.GetAt(0, 0)
		p3 := AR.GetAt(index, index)
		AR.SetAt(0, 0, p3)
		AR.SetAt(index, index, p1)
		return
	}
	if flags&UPPER != 0 {
		// AL is merged from [ATL, ATR], AR is [ABR]; P1 is AR[0, 0]; P2 is AL[index, -1]
		colno := AL.Cols() - AR.Cols()
		// S1 -- D1; S1 is on the first row of AR
		AR.SubMatrix(&s, 0, 1, 1, AR.Cols()-1)
		AL.SubMatrix(&d, index, colno+1, 1, s.Cols())
		Swap(&s, &d)
		// S2 -- D2
		AL.SubMatrix(&s, index+1, colno, AL.Rows()-index-2, 1)
		AL.SubMatrix(&d, index, index+1, 1, colno-index-1)
		Swap(&s, &d)
		// S3 -- D3
		AL.SubMatrix(&s, 0, index, index, 1)
		AL.SubMatrix(&d, 0, colno, index, 1)
		Swap(&s, &d)
		//fmt.Printf("3, AR=%v\n", AR)
		// swap P1 and P3
		p1 := AR.GetAt(0, 0)
		p3 := AL.GetAt(index, index)
		AR.SetAt(0, 0, p3)
		AL.SetAt(index, index, p1)
		return
	}
}

/*
 * Apply diagonal pivot (row and column swapped) to symmetric matrix blocks.
 * AR[0,0] is on diagonal and AL is block to the left of diagonal and AR the
 * triangular diagonal block. Need to swap row and column.
 *
 * LOWER triangular; moving from top-left to bottom-right
 *
 *    d
 *    x  d |
 *    --------------------------
 *    x  x | d
 *    S1 S1| S1 P1 x  x  x  P2     -- current row/col 'srcix'
 *    x  x | x  S2 d  x  x  x
 *    x  x | x  S2 x  d  x  x
 *    x  x | x  S2 x  x  d  x
 *    D1 D1| D1 P2 D2 D2 D2 P3     -- swap with row/col 'dstix'
 *    x  x | x  S3 x  x  x  D3 d
 *    x  x | x  S3 x  x  x  D3 x d
 *    (ABL)          (ABR)
 *
 * UPPER triangular; moving from bottom-right to top-left
 *
 *         (ATL)                  (ATR)
 *    d  x  x  D3 x  x  x  S3 x | x
 *       d  x  D3 x  x  x  S3 x | x
 *          d  D3 x  x  x  S3 x | x
 *             P3 D2 D2 D2 P2 D1| D1  -- dstinx
 *                d  x  x  S2 x | x
 *                   d  x  S2 x | x
 *                      d  S2 x | x
 *                         P1 S1| S1  -- srcinx
 *                            d | x
 *    -----------------------------
 *                              | d
 *                           (ABR)
 */
func applyPivotSym2(AL, AR *matrix.FloatMatrix, srcix, dstix int, flags Flags) {
	var s, d matrix.FloatMatrix
	if flags&LOWER != 0 {
		// AL is [ABL]; AR is [ABR]; P1 is AR[0,0], P2 is AR[index, 0]
		// S1 -- D1
		AL.SubMatrix(&s, srcix, 0, 1, AL.Cols())
		AL.SubMatrix(&d, dstix, 0, 1, AL.Cols())
		Swap(&s, &d)
		if srcix > 0 {
			AR.SubMatrix(&s, srcix, 0, 1, srcix)
			AR.SubMatrix(&d, dstix, 0, 1, srcix)
			Swap(&s, &d)
		}
		// S2 -- D2
		AR.SubMatrix(&s, srcix+1, srcix, dstix-srcix-1, 1)
		AR.SubMatrix(&d, dstix, srcix+1, 1, dstix-srcix-1)
		Swap(&s, &d)
		// S3 -- D3
		AR.SubMatrix(&s, dstix+1, srcix, AR.Rows()-dstix-1, 1)
		AR.SubMatrix(&d, dstix+1, dstix, AR.Rows()-dstix-1, 1)
		Swap(&s, &d)
		// swap P1 and P3
		p1 := AR.GetAt(srcix, srcix)
		p3 := AR.GetAt(dstix, dstix)
		AR.SetAt(srcix, srcix, p3)
		AR.SetAt(dstix, dstix, p1)
		return
	}
	if flags&UPPER != 0 {
		// AL is ATL, AR is ATR; P1 is AL[srcix, srcix];
		// S1 -- D1;
		AR.SubMatrix(&s, srcix, 0, 1, AR.Cols())
		AR.SubMatrix(&d, dstix, 0, 1, AR.Cols())
		Swap(&s, &d)
		if srcix < AL.Cols()-1 {
			// not the corner element
			AL.SubMatrix(&s, srcix, srcix+1, 1, srcix)
			AL.SubMatrix(&d, dstix, srcix+1, 1, srcix)
			Swap(&s, &d)
		}
		// S2 -- D2
		AL.SubMatrix(&s, dstix+1, srcix, srcix-dstix-1, 1)
		AL.SubMatrix(&d, dstix, dstix+1, 1, srcix-dstix-1)
		Swap(&s, &d)
		// S3 -- D3
		AL.SubMatrix(&s, 0, srcix, dstix, 1)
		AL.SubMatrix(&d, 0, dstix, dstix, 1)
		Swap(&s, &d)
		//fmt.Printf("3, AR=%v\n", AR)
		// swap P1 and P3
		p1 := AR.GetAt(0, 0)
		p3 := AL.GetAt(dstix, dstix)
		AR.SetAt(srcix, srcix, p3)
		AL.SetAt(dstix, dstix, p1)
		return
	}
}

/*
 * Apply diagonal pivot (row and column swapped) to symmetric matrix blocks.
 *
 * LOWER triangular; moving from top-left to bottom-right
 *
 *    -----------------------
 *    | d
 *    | x P1 x  x  x  P2     -- current row/col 'srcix'
 *    | x S2 d  x  x  x
 *    | x S2 x  d  x  x
 *    | x S2 x  x  d  x
 *    | x P2 D2 D2 D2 P3     -- swap with row/col 'dstix'
 *    | x S3 x  x  x  D3 d
 *    | x S3 x  x  x  D3 x d
 *         (AR)
 *
 * UPPER triangular; moving from bottom-right to top-left
 *
 *    d x D3 x  x  x  S3 x |
 *      d D3 x  x  x  S3 x |
 *        P3 D2 D2 D2 P2 x |  -- dstinx
 *           d  x  x  S2 x |
 *              d  x  S2 x |
 *                 d  S2 x |
 *                    P1 x |  -- srcinx
 *                       d |
 *    ----------------------
 *               (ABR)
 */
func applyBKPivotSym(AR *matrix.FloatMatrix, srcix, dstix int, flags Flags) {
	var s, d matrix.FloatMatrix
	if flags&LOWER != 0 {
		// S2 -- D2
		AR.SubMatrix(&s, srcix+1, srcix, dstix-srcix-1, 1)
		AR.SubMatrix(&d, dstix, srcix+1, 1, dstix-srcix-1)
		Swap(&s, &d)
		// S3 -- D3
		AR.SubMatrix(&s, dstix+1, srcix, AR.Rows()-dstix-1, 1)
		AR.SubMatrix(&d, dstix+1, dstix, AR.Rows()-dstix-1, 1)
		Swap(&s, &d)
		// swap P1 and P3
		p1 := AR.GetAt(srcix, srcix)
		p3 := AR.GetAt(dstix, dstix)
		AR.SetAt(srcix, srcix, p3)
		AR.SetAt(dstix, dstix, p1)
		return
	}
	if flags&UPPER != 0 {
		// AL is ATL, AR is ATR; P1 is AL[srcix, srcix];
		// S2 -- D2
		AR.SubMatrix(&s, dstix+1, srcix, srcix-dstix-1, 1)
		AR.SubMatrix(&d, dstix, dstix+1, 1, srcix-dstix-1)
		Swap(&s, &d)
		// S3 -- D3
		AR.SubMatrix(&s, 0, srcix, dstix, 1)
		AR.SubMatrix(&d, 0, dstix, dstix, 1)
		Swap(&s, &d)
		//fmt.Printf("3, AR=%v\n", AR)
		// swap P1 and P3
		p1 := AR.GetAt(srcix, srcix)
		p3 := AR.GetAt(dstix, dstix)
		AR.SetAt(srcix, srcix, p3)
		AR.SetAt(dstix, dstix, p1)
		return
	}
}

/* From LAPACK/dlarf.f
 *
 * 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
 *
 * Allocates/frees intermediate work space matrix w1.
 */
func applyHouseholder(tau, v, a1, A2 *matrix.FloatMatrix, flags Flags) {

	tval := tau.GetAt(0, 0)
	if tval == 0.0 {
		return
	}
	w1 := a1.Copy()
	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
	MVRankUpdate(A2, v, w1, -1.0)
}

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
}

/*
 * Compute
 *   C = C*diag(D)      flags & RIGHT == true
 *   C = diag(D)*C      flags & LEFT  == true
 *
 * Arguments
 *   C     M-by-N matrix if flags&RIGHT == true or N-by-M matrix if flags&LEFT == true
 *
 *   D     N element column or row vector or N-by-N matrix
 *
 *   flags Indicator bits, LEFT or RIGHT
 */
func MultDiag(C, D *matrix.FloatMatrix, flags Flags) {
	var c, d0 matrix.FloatMatrix
	if D.Cols() == 1 {
		// diagonal is column vector
		switch flags & (LEFT | RIGHT) {
		case LEFT:
			// scale rows; for each column element-wise multiply with D-vector
			for k := 0; k < C.Cols(); k++ {
				C.SubMatrix(&c, 0, k, C.Rows(), 1)
				c.Mul(D)
			}
		case RIGHT:
			// scale columns
			for k := 0; k < C.Cols(); k++ {
				C.SubMatrix(&c, 0, k, C.Rows(), 1)
				// scale the column
				c.Scale(D.GetAt(k, 0))
			}
		}
	} else {
		// diagonal is row vector
		var d *matrix.FloatMatrix
		if D.Rows() == 1 {
			d = D
		} else {
			D.SubMatrix(&d0, 0, 0, 1, D.Cols(), D.LeadingIndex()+1)
			d = &d0
		}
		switch flags & (LEFT | RIGHT) {
		case LEFT:
			for k := 0; k < C.Rows(); k++ {
				C.SubMatrix(&c, k, 0, 1, C.Cols())
				// scale the row
				c.Scale(d.GetAt(0, k))
			}
		case RIGHT:
			// scale columns
			for k := 0; k < C.Cols(); k++ {
				C.SubMatrix(&c, 0, k, C.Rows(), 1)
				// scale the column
				c.Scale(d.GetAt(0, k))
			}
		}
	}
}

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
//.........這裏部分代碼省略.........

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
}

/*
 * 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))
//.........這裏部分代碼省略.........