Golang FloatMatrix.Data方法代码示例

/*
 * Compute LDL^T factorization of real symmetric matrix.
 *
 * Computes of a real symmetric matrix A using Bunch-Kauffman pivoting method.
 * The form of factorization is
 *
 *    A = L*D*L.T  or A = U*D*U.T
 *
 * where L (or U) is product of permutation and unit lower (or upper) triangular matrix
 * and D is block diagonal symmetric matrix with 1x1 and 2x2 blocks.
 *
 * Arguments
 *  A     On entry, the N-by-N symmetric matrix A. If flags bit LOWER (or UPPER) is set then
 *        lower (or upper) triangular matrix and strictly upper (or lower) part is not
 *        accessed. On exit, the block diagonal matrix D and lower (or upper) triangular
 *        product matrix L (or U).
 *
 *  W     Workspace, size as returned by WorksizeBK().
 *
 *  ipiv  Pivot vector. On exit details of interchanges and the block structure of D. If
 *        ipiv[k] > 0 then D[k,k] is 1x1 and rows and columns k and ipiv[k]-1 were changed.
 *        If ipiv[k] == ipiv[k+1] < 0 then D[k,k] is 2x2. If A is lower then rows and
 *        columns k+1 and ipiv[k]-1  were changed. And if A is upper then rows and columns
 *        k and ipvk[k]-1 were changed.
 *
 *  flags Indicator bits, LOWER or UPPER.
 *
 *  confs Optional blocking configuration. If not provided then default blocking
 *        as returned by DefaultConf() is used.
 *
 *  Unblocked algorithm is used if blocking configuration LB is zero or if N < LB.
 *
 *  Compatible with lapack.SYTRF.
 */
func BKFactor(A, W *cmat.FloatMatrix, ipiv Pivots, flags int, confs ...*gomas.Config) *gomas.Error {
	var err *gomas.Error = nil
	conf := gomas.CurrentConf(confs...)

	for k, _ := range ipiv {
		ipiv[k] = 0
	}
	wsz := BKFactorWork(A, conf)
	if W.Len() < wsz {
		return gomas.NewError(gomas.EWORK, "DecomposeBK", wsz)
	}

	var Wrk cmat.FloatMatrix
	if n(A) < conf.LB || conf.LB == 0 {
		// make workspace rows(A)*2 matrix
		Wrk.SetBuf(m(A), 2, m(A), W.Data())
		if flags&gomas.LOWER != 0 {
			err, _ = unblkDecompBKLower(A, &Wrk, ipiv, conf)
		} else if flags&gomas.UPPER != 0 {
			err, _ = unblkDecompBKUpper(A, &Wrk, ipiv, conf)
		}
	} else {
		// make workspace rows(A)*(LB+1) matrix
		Wrk.SetBuf(m(A), conf.LB+1, m(A), W.Data())
		if flags&gomas.LOWER != 0 {
			err = blkDecompBKLower(A, &Wrk, &ipiv, conf)
		} else if flags&gomas.UPPER != 0 {
			err = blkDecompBKUpper(A, &Wrk, &ipiv, conf)
		}
	}
	return err
}

/*
 * Compute RQ factorization of a M-by-N matrix A: A = R*Q
 *
 * Arguments:
 *  A    On entry, the M-by-N matrix A, M <= N. On exit, upper triangular matrix R
 *       and the orthogonal matrix Q as product of elementary reflectors.
 *
 * tau  On exit, the scalar factors of the elementary reflectors.
 *
 * W    Workspace, M-by-nb matrix used for work space in blocked invocations.
 *
 * conf The blocking configuration. If nil then default blocking configuration
 *      is used. Member conf.LB defines blocking size of blocked algorithms.
 *      If it is zero then unblocked algorithm is used.
 *
 * Returns:
 *      Error indicator.
 *
 * Additional information
 *
 *  Ortogonal matrix Q is product of elementary reflectors H(k)
 *
 *    Q = H(0)H(1),...,H(K-1), where K = min(M,N)
 *
 *  Elementary reflector H(k) is stored on row k of A right of the diagonal with
 *  implicit unit value on diagonal entry. The vector TAU holds scalar factors of
 *  the elementary reflectors.
 *
 *  Contents of matrix A after factorization is as follow:
 *
 *    ( v0 v0 r  r  r  r )  M=4, N=6
 *    ( v1 v1 v1 r  r  r )
 *    ( v2 v2 v2 v2 r  r )
 *    ( v3 v3 v3 v3 v3 r )
 *
 *  where l is element of L, vk is element of H(k).
 *
 *  RQFactor is compatible with lapack.DGERQF
 */
func RQFactor(A, tau, W *cmat.FloatMatrix, confs ...*gomas.Config) *gomas.Error {
	var err *gomas.Error = nil
	conf := gomas.CurrentConf(confs...)

	// must have: M <= N
	if m(A) > n(A) {
		return gomas.NewError(gomas.ESIZE, "RQFactor")
	}

	wsmin := wsLQ(A, 0)
	if W == nil || W.Len() < wsmin {
		return gomas.NewError(gomas.EWORK, "RQFactor", wsmin)
	}
	lb := estimateLB(A, W.Len(), wsRQ)
	lb = imin(lb, conf.LB)
	if lb == 0 || m(A) <= lb {
		unblockedRQ(A, tau, W)
	} else {
		var Twork, Wrk cmat.FloatMatrix
		// block reflector T in first LB*LB elements in workspace
		// the rest, m(A)-LB*LB, is workspace for intermediate matrix operands
		Twork.SetBuf(lb, lb, lb, W.Data())
		Wrk.SetBuf(m(A)-lb, lb, m(A)-lb, W.Data()[Twork.Len():])
		blockedRQ(A, tau, &Twork, &Wrk, lb, conf)
	}
	return err
}

/*
 * Generate the M by N matrix Q with orthogonal rows which
 * are defined as the first M rows of the product of K first elementary
 * reflectors.
 *
 * Arguments
 *   A     On entry, the elementary reflectors as returned by LQFactor().
 *         stored right of diagonal of the M by N matrix A.
 *         On exit, the orthogonal matrix Q
 *
 *   tau   Scalar coefficents of elementary reflectors
 *
 *   W     Workspace
 *
 *   K     The number of elementary reflector whose product define the matrix Q
 *
 *   conf  Optional blocking configuration.
 *
 * Compatible with lapackd.ORGLQ.
 */
func LQBuild(A, tau, W *cmat.FloatMatrix, K int, confs ...*gomas.Config) *gomas.Error {
	var err *gomas.Error = nil
	conf := gomas.CurrentConf(confs...)
	if K <= 0 || K > n(A) {
		return gomas.NewError(gomas.EVALUE, "LQBuild", K)
	}
	wsz := wsBuildLQ(A, 0)
	if W == nil || W.Len() < wsz {
		return gomas.NewError(gomas.EWORK, "LQBuild", wsz)
	}

	// adjust blocking factor for workspace size
	lb := estimateLB(A, W.Len(), wsBuildLQ)
	//lb = imin(lb, conf.LB)
	lb = conf.LB
	if lb == 0 || m(A) <= lb {
		unblkBuildLQ(A, tau, W, m(A)-K, n(A)-K, true)
	} else {
		var Twork, Wrk cmat.FloatMatrix
		Twork.SetBuf(lb, lb, lb, W.Data())
		Wrk.SetBuf(m(A)-lb, lb, m(A)-lb, W.Data()[Twork.Len():])
		blkBuildLQ(A, tau, &Twork, &Wrk, K, lb, conf)
	}
	return err
}

/*
 * Compute QL factorization of a M-by-N matrix A: A = Q * L.
 *
 * Arguments:
 *  A    On entry, the M-by-N matrix A, M >= N. On exit, lower triangular matrix L
 *       and the orthogonal matrix Q as product of elementary reflectors.
 *
 * tau  On exit, the scalar factors of the elemenentary reflectors.
 *
 * W    Workspace, N-by-nb matrix used for work space in blocked invocations.
 *
 * conf The blocking configuration. If nil then default blocking configuration
 *      is used. Member conf.LB defines blocking size of blocked algorithms.
 *      If it is zero then unblocked algorithm is used.
 *
 * Returns:
 *      Error indicator.
 *
 * Additional information
 *
 *  Ortogonal matrix Q is product of elementary reflectors H(k)
 *
 *    Q = H(K-1)...H(1)H(0), where K = min(M,N)
 *
 *  Elementary reflector H(k) is stored on column k of A above the diagonal with
 *  implicit unit value on diagonal entry. The vector TAU holds scalar factors
 *  of the elementary reflectors.
 *
 *  Contents of matrix A after factorization is as follow:
 *
 *    ( v0 v1 v2 v3 )   for M=6, N=4
 *    ( v0 v1 v2 v3 )
 *    ( l  v1 v2 v3 )
 *    ( l  l  v2 v3 )
 *    ( l  l  l  v3 )
 *    ( l  l  l  l  )
 *
 *  where l is element of L, vk is element of H(k).
 *
 * DecomposeQL is compatible with lapack.DGEQLF
 */
func QLFactor(A, tau, W *cmat.FloatMatrix, confs ...*gomas.Config) *gomas.Error {
	var err *gomas.Error = nil
	var tauh cmat.FloatMatrix
	conf := gomas.CurrentConf(confs...)

	if m(A) < n(A) {
		return gomas.NewError(gomas.ESIZE, "QLFactor")
	}
	wsmin := wsQL(A, 0)
	if W == nil || W.Len() < wsmin {
		return gomas.NewError(gomas.EWORK, "QLFactor", wsmin)
	}
	if tau.Len() < n(A) {
		return gomas.NewError(gomas.ESIZE, "QLFactor")
	}
	tauh.SubMatrix(tau, 0, 0, n(A), 1)
	lb := estimateLB(A, W.Len(), wsQL)
	lb = imin(lb, conf.LB)

	if lb == 0 || n(A) <= lb {
		unblockedQL(A, &tauh, W)
	} else {
		var Twork, Wrk cmat.FloatMatrix
		// block reflector T in first LB*LB elements in workspace
		// the rest, n(A)-LB*LB, is workspace for intermediate matrix operands
		Twork.SetBuf(conf.LB, conf.LB, -1, W.Data())
		Wrk.SetBuf(n(A)-conf.LB, conf.LB, -1, W.Data()[Twork.Len():])
		blockedQL(A, &tauh, &Twork, &Wrk, lb, conf)
	}
	return err
}

/*
 * Compute QR factorization of a M-by-N matrix A using compact WY transformation: A = Q * R,
 * where Q = I - Y*T*Y.T, T is block reflector and Y holds elementary reflectors as lower
 * trapezoidal matrix saved below diagonal elements of the matrix A.
 *
 * Arguments:
 *  A    On entry, the M-by-N matrix A. On exit, the elements on and above
 *       the diagonal contain the min(M,N)-by-N upper trapezoidal matrix R.
 *       The elements below the diagonal with the matrix 'T', represent
 *       the ortogonal matrix Q as product of elementary reflectors.
 *
 * T     On exit, the K block reflectors which, together with trilu(A) represent
 *       the ortogonal matrix Q as Q = I - Y*T*Y.T where Y = trilu(A).
 *       K is ceiling(N/LB) where LB is blocking size from used blocking configuration.
 *       The matrix T is LB*N augmented matrix of K block reflectors,
 *       T = [T(0) T(1) .. T(K-1)].  Block reflector T(n) is LB*LB matrix, expect
 *       reflector T(K-1) that is IB*IB matrix  where IB = min(LB, K % LB)
 *
 * W     Workspace, required size returned by QRTFactorWork().
 *
 * conf  Optional blocking configuration. If not provided then default configuration
 *       is used.
 *
 * Returns:
 *      Error indicator.
 *
 * QRTFactor is compatible with lapack.DGEQRT
 */
func QRTFactor(A, T, W *cmat.FloatMatrix, confs ...*gomas.Config) *gomas.Error {
	var err *gomas.Error = nil
	conf := gomas.CurrentConf(confs...)
	ok := false
	rsize := 0

	if m(A) < n(A) {
		return gomas.NewError(gomas.ESIZE, "QRTFactor")
	}
	wsz := QRTFactorWork(A, conf)
	if W == nil || W.Len() < wsz {
		return gomas.NewError(gomas.EWORK, "QRTFactor", wsz)
	}

	tr, tc := T.Size()
	if conf.LB == 0 || conf.LB > n(A) {
		ok = tr == tc && tr == n(A)
		rsize = n(A) * n(A)
	} else {
		ok = tr == conf.LB && tc == n(A)
		rsize = conf.LB * n(A)
	}
	if !ok {
		return gomas.NewError(gomas.ESMALL, "QRTFactor", rsize)
	}

	if conf.LB == 0 || n(A) <= conf.LB {
		err = unblockedQRT(A, T, W)
	} else {
		Wrk := cmat.MakeMatrix(n(A), conf.LB, W.Data())
		err = blockedQRT(A, T, Wrk, conf)
	}
	return err
}

func axpby(Y, X *cmat.FloatMatrix, alpha, beta float64, N int) {
	var x, y C.mvec_t

	xr, _ := X.Size()
	x.md = (*C.double)(unsafe.Pointer(&X.Data()[0]))
	x.inc = C.int(1)
	if xr == 1 {
		x.inc = C.int(X.Stride())
	}
	yr, _ := Y.Size()
	y.md = (*C.double)(unsafe.Pointer(&Y.Data()[0]))
	y.inc = C.int(1)
	if yr == 1 {
		y.inc = C.int(Y.Stride())
	}
	if beta == 1.0 {
		C.__d_vec_axpy(
			(*C.mvec_t)(unsafe.Pointer(&y)),
			(*C.mvec_t)(unsafe.Pointer(&x)),
			C.double(alpha), C.int(N))
	} else {
		C.__d_vec_axpby(
			(*C.mvec_t)(unsafe.Pointer(&y)),
			(*C.mvec_t)(unsafe.Pointer(&x)),
			C.double(alpha), C.double(beta), C.int(N))
	}
	return
}

/*
 * Reduce general matrix A to upper Hessenberg form H by similiarity
 * transformation H = Q.T*A*Q.
 *
 * Arguments:
 *  A    On entry, the general matrix A. On exit, the elements on and
 *       above the first subdiagonal contain the reduced matrix H.
 *       The elements below the first subdiagonal with the vector tau
 *       represent the ortogonal matrix A as product of elementary reflectors.
 *
 *  tau  On exit, the scalar factors of the elementary reflectors.
 *
 *  W    Workspace, as defined by HessReduceWork()
 *
 *  conf The blocking configration.
 *
 * HessReduce is compatible with lapack.DGEHRD.
 */
func HessReduce(A, tau, W *cmat.FloatMatrix, confs ...*gomas.Config) *gomas.Error {
	var err *gomas.Error = nil
	conf := gomas.CurrentConf(confs...)

	wmin := m(A)
	wopt := HessReduceWork(A, conf)
	wsz := W.Len()
	if wsz < wmin {
		return gomas.NewError(gomas.EWORK, "ReduceHess", wmin)
	}
	// use blocked version if workspace big enough for blocksize 4
	lb := conf.LB
	if wsz < wopt {
		lb = estimateLB(A, wsz, wsHess)
	}
	if lb == 0 || n(A) <= lb {
		unblkHessGQvdG(A, tau, W, 0)
	} else {
		// blocked version
		var W0 cmat.FloatMatrix
		// shape workspace for blocked algorithm
		W0.SetBuf(m(A)+lb, lb, m(A)+lb, W.Data())
		blkHessGQvdG(A, tau, &W0, lb, conf)
	}
	return err
}

/*
 * Multiply and replace C with Q*C or Q.T*C where Q is a real orthogonal matrix
 * defined as the product of k elementary reflectors and block reflector T
 *
 *    Q = H(1) H(2) . . . H(k)
 *
 * as returned by DecomposeQRT().
 *
 * Arguments:
 *  C     On entry, the M-by-N matrix C. On exit C is overwritten by Q*C or Q.T*C.
 *
 *  A     QR factorization as returned by QRTFactor() where the lower trapezoidal
 *        part holds the elementary reflectors.
 *
 *  T     The block reflector computed from elementary reflectors as returned by
 *        DecomposeQRT() or computed from elementary reflectors and scalar coefficients
 *        by BuildT()
 *
 *  W     Workspace, size as returned by QRTMultWork()
 *
 *  conf  Blocking configuration
 *
 *  flags Indicators. Valid indicators LEFT, RIGHT, TRANS, NOTRANS
 *
 * Preconditions:
 *   a.   cols(A) == cols(T),
 *          columns A define number of elementary reflector, must match order of block reflector.
 *   b.   if conf.LB == 0, cols(T) == rows(T)
 *          unblocked invocation, block reflector T is upper triangular
 *   c.   if conf.LB != 0, rows(T) == conf.LB
 *          blocked invocation, T is sequence of triangular block reflectors of order LB
 *   d.   if LEFT, rows(C) >= cols(A) && cols(C) >= rows(A)
 *
 *   e.   if RIGHT, cols(C) >= cols(A) && rows(C) >= rows(A)
 *
 * Compatible with lapack.DGEMQRT
 */
func QRTMult(C, A, T, W *cmat.FloatMatrix, flags int, confs ...*gomas.Config) *gomas.Error {
	var err *gomas.Error = nil
	conf := gomas.CurrentConf(confs...)

	wsz := QRTMultWork(C, T, flags, conf)
	if W == nil || W.Len() < wsz {
		return gomas.NewError(gomas.EWORK, "QRTMult", wsz)
	}
	ok := false
	switch flags & gomas.RIGHT {
	case gomas.RIGHT:
		ok = n(C) >= m(A)
	default:
		ok = m(C) >= n(A)
	}
	if !ok {
		return gomas.NewError(gomas.ESIZE, "QRTMult")
	}

	var Wrk cmat.FloatMatrix
	if flags&gomas.RIGHT != 0 {
		Wrk.SetBuf(m(C), conf.LB, m(C), W.Data())
		blockedMultQTRight(C, A, T, &Wrk, flags, conf)

	} else {
		Wrk.SetBuf(n(C), conf.LB, n(C), W.Data())
		blockedMultQTLeft(C, A, T, &Wrk, flags, conf)
	}
	return err
}

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

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

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

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

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

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

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

func minvscale(A *cmat.FloatMatrix, alpha float64, M, N int) {
	var a C.mdata_t

	a.md = (*C.double)(unsafe.Pointer(&A.Data()[0]))
	a.step = C.int(A.Stride())
	C.__d_blk_invscale(
		(*C.mdata_t)(unsafe.Pointer(&a)), C.double(alpha), C.int(M), C.int(N))
	return
}

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

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

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

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

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

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

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

func vinvscal(X *cmat.FloatMatrix, alpha float64, N int) {
	var x C.mvec_t

	xr, _ := X.Size()
	x.md = (*C.double)(unsafe.Pointer(&X.Data()[0]))
	x.inc = C.int(1)
	if xr == 1 {
		x.inc = C.int(X.Stride())
	}
	C.__d_vec_invscal(
		(*C.mvec_t)(unsafe.Pointer(&x)), C.double(alpha), C.int(N))
	return
}

func plus(Ac, Bc *cmat.FloatMatrix, alpha, beta float64, bits, S, L, R, E int) {
	var Am, Bm C.mdata_t

	Am.md = (*C.double)(unsafe.Pointer(&Ac.Data()[0]))
	Am.step = C.int(Ac.Stride())
	Bm.md = (*C.double)(unsafe.Pointer(&Bc.Data()[0]))
	Bm.step = C.int(Bc.Stride())

	C.__d_scale_plus(
		(*C.mdata_t)(unsafe.Pointer(&Am)),
		(*C.mdata_t)(unsafe.Pointer(&Bm)),
		C.double(alpha), C.double(beta), C.int(bits),
		C.int(S), C.int(L), C.int(R), C.int(E))
}

func sum(X *cmat.FloatMatrix, N int) float64 {
	var x C.mvec_t
	var dc C.double

	xr, _ := X.Size()
	x.md = (*C.double)(unsafe.Pointer(&X.Data()[0]))
	x.inc = C.int(1)
	if xr == 1 {
		x.inc = C.int(X.Stride())
	}
	dc = C.__d_vec_sum_recursive(
		(*C.mvec_t)(unsafe.Pointer(&x)), C.int(N))
	return float64(dc)
}

func iamax(X *cmat.FloatMatrix, N int) int {
	var x C.mvec_t
	var ix C.int

	xr, _ := X.Size()
	x.md = (*C.double)(unsafe.Pointer(&X.Data()[0]))
	x.inc = C.int(1)
	if xr == 1 {
		x.inc = C.int(X.Stride())
	}
	ix = C.__d_vec_iamax(
		(*C.mvec_t)(unsafe.Pointer(&x)), C.int(N))
	return int(ix)
}