Golang FloatMatrix.Len方法代码示例

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

/*
 * Compute
 *   B = B*diag(D).-1      flags & RIGHT == true
 *   B = diag(D).-1*B      flags & LEFT  == true
 *
 * If flags is LEFT (RIGHT) then element-wise divides columns (rows) of B with vector D.
 *
 * Arguments:
 *   B     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 SolveDiag(B, D *cmat.FloatMatrix, flags int, confs ...*gomas.Config) *gomas.Error {
	var c, d0 cmat.FloatMatrix
	var d *cmat.FloatMatrix

	conf := gomas.CurrentConf(confs...)
	d = D
	if !D.IsVector() {
		d0.Diag(D)
		d = &d0
	}
	dn := d0.Len()
	br, bc := B.Size()
	switch flags & (gomas.LEFT | gomas.RIGHT) {
	case gomas.LEFT:
		if br != dn {
			return gomas.NewError(gomas.ESIZE, "SolveDiag")
		}
		// scale rows;
		for k := 0; k < dn; k++ {
			c.Row(B, k)
			blasd.InvScale(&c, d.GetAt(k), conf)
		}
	case gomas.RIGHT:
		if bc != dn {
			return gomas.NewError(gomas.ESIZE, "SolveDiag")
		}
		// scale columns
		for k := 0; k < dn; k++ {
			c.Column(B, k)
			blasd.InvScale(&c, d.GetAt(k), conf)
		}
	}
	return nil
}

// Compute i'th updated element of rank-1 update vector
func trdsecUpdateElemDelta(d, delta *cmat.FloatMatrix, index int, rho float64) float64 {
	var n0, n1, dn, dk, val, p0, p1 float64
	var k, N int

	N = d.Len()
	dk = d.GetAt(index)
	dn = delta.GetAt(N - 1)

	// compute; prod j; (lambda_j - d_k)/(d_j - d_k), j = 0 .. index-1
	p0 = 1.0
	for k = 0; k < index; k++ {
		n0 = delta.GetAt(k)
		n1 = d.GetAt(k) - dk
		p0 = p0 * (n0 / n1)
	}

	p1 = 1.0
	for k = index; k < N-1; k++ {
		n0 = delta.GetAt(k)
		n1 = d.GetAt(k+1) - dk
		p1 = p1 * (n0 / n1)
	}
	val = p0 * p1 * (dn / rho)
	return math.Sqrt(math.Abs(val))
}

func sortEigenVec(D, U, V, C *cmat.FloatMatrix, updown int) {
	var sD, m0, m1 cmat.FloatMatrix

	N := D.Len()
	for k := 0; k < N-1; k++ {
		sD.SubVector(D, k, N-k)
		pk := vecMinMax(&sD, -updown)
		if pk != 0 {
			t0 := D.GetAt(k)
			D.SetAt(k, D.GetAt(pk+k))
			D.SetAt(k+pk, t0)
			if U != nil {
				m0.Column(U, k)
				m1.Column(U, k+pk)
				blasd.Swap(&m1, &m0)
			}
			if V != nil {
				m0.Row(V, k)
				m1.Row(V, k+pk)
				blasd.Swap(&m1, &m0)
			}
			if C != nil {
				m0.Column(C, k)
				m1.Column(C, k+pk)
				blasd.Swap(&m1, &m0)
			}
		}
	}
}

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

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

/*
 * Reduce a general M-by-N matrix A to upper or lower bidiagonal form B
 * by an ortogonal transformation A = Q*B*P.T,  B = Q.T*A*P
 *
 *
 * Arguments
 *   A     On entry, the real M-by-N matrix. On exit the upper/lower
 *         bidiagonal matrix and ortogonal matrices Q and P.
 *
 *   tauq  Scalar factors for elementary reflector forming the
 *         ortogonal matrix Q.
 *
 *   taup  Scalar factors for elementary reflector forming the
 *         ortogonal matrix P.
 *
 *   W     Workspace needed for reduction.
 *
 *   conf  Current blocking configuration. Optional.
 *
 *
 * Details
 *
 * Matrices Q and P are products of elementary reflectors H(k) and G(k)
 *
 * If M > N:
 *     Q = H(1)*H(2)*...*H(N)   and P = G(1)*G(2)*...*G(N-1)
 *
 * where H(k) = 1 - tauq*u*u.T and G(k) = 1 - taup*v*v.T
 *
 * Elementary reflector H(k) are stored on columns of A below the diagonal with
 * implicit unit value on diagonal entry. Vector TAUQ holds corresponding scalar
 * factors. Reflector G(k) are stored on rows of A right of first superdiagonal
 * with implicit unit value on superdiagonal. Corresponding scalar factors are
 * stored on vector TAUP.
 *
 * If M < N:
 *   Q = H(1)*H(2)*...*H(N-1)   and P = G(1)*G(2)*...*G(N)
 *
 * where H(k) = 1 - tauq*u*u.T and G(k) = 1 - taup*v*v.T
 *
 * Elementary reflector H(k) are stored on columns of A below the first sub diagonal
 * with implicit unit value on sub diagonal entry. Vector TAUQ holds corresponding
 * scalar factors. Reflector G(k) are sotre on rows of A right of diagonal with
 * implicit unit value on superdiagonal. Corresponding scalar factors are stored
 * on vector TAUP.
 *
 * Contents of matrix A after reductions are as follows.
 *
 *    M = 6 and N = 5:                  M = 5 and N = 6:
 *
 *    (  d   e   v1  v1  v1 )           (  d   v1  v1  v1  v1  v1 )
 *    (  u1  d   e   v2  v2 )           (  e   d   v2  v2  v2  v2 )
 *    (  u1  u2  d   e   v3 )           (  u1  e   d   v3  v3  v3 )
 *    (  u1  u2  u3  d   e  )           (  u1  u2  e   d   v4  v4 )
 *    (  u1  u2  u3  u4  d  )           (  u1  u2  u3  e   d   v5 )
 *    (  u1  u2  u3  u4  u5 )
 */
func BDReduce(A, tauq, taup, W *cmat.FloatMatrix, confs ...*gomas.Config) *gomas.Error {
	var err *gomas.Error = nil
	conf := gomas.CurrentConf(confs...)
	_ = conf

	wmin := wsBired(A, 0)
	wsz := W.Len()
	if wsz < wmin {
		return gomas.NewError(gomas.EWORK, "ReduceBidiag", wmin)
	}
	lb := conf.LB
	wneed := wsBired(A, lb)
	if wneed > wsz {
		lb = estimateLB(A, wsz, wsBired)
	}
	if m(A) >= n(A) {
		if lb > 0 && n(A) > lb {
			blkBidiagLeft(A, tauq, taup, W, lb, conf)
		} else {
			unblkReduceBidiagLeft(A, tauq, taup, W)
		}
	} else {
		if lb > 0 && m(A) > lb {
			blkBidiagRight(A, tauq, taup, W, lb, conf)
		} else {
			unblkReduceBidiagRight(A, tauq, taup, W)
		}
	}
	return err
}

func MVMult(Y, A, X *cmat.FloatMatrix, alpha, beta float64, bits int, confs ...*gomas.Config) *gomas.Error {
	ok := true
	yr, yc := Y.Size()
	ar, ac := A.Size()
	xr, xc := X.Size()

	if ar*ac == 0 {
		return nil
	}
	if yr != 1 && yc != 1 {
		return gomas.NewError(gomas.ENEED_VECTOR, "MVMult")
	}
	if xr != 1 && xc != 1 {
		return gomas.NewError(gomas.ENEED_VECTOR, "MVMult")
	}
	nx := X.Len()
	ny := Y.Len()

	if bits&gomas.TRANSA != 0 {
		bits |= gomas.TRANS
	}
	if bits&gomas.TRANS != 0 {
		ok = ny == ac && nx == ar
	} else {
		ok = ny == ar && nx == ac
	}
	if !ok {
		return gomas.NewError(gomas.ESIZE, "MVMult")
	}
	if beta != 1.0 {
		vscal(Y, beta, ny)
	}
	gemv(Y, A, X, alpha, beta, bits, 0, nx, 0, ny)
	return nil
}

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

// d = |d| - |s|
func absMinus(d, s *cmat.FloatMatrix) *cmat.FloatMatrix {
	for k := 0; k < d.Len(); k++ {
		tmp := math.Abs(d.GetAt(k))
		d.SetAt(k, math.Abs(s.GetAt(k))-tmp)
	}
	return d
}

/*
 * Symmetric matrix-vector multiplication. Y = beta*Y + alpha*A*X
 */
func MVMultSym(Y, A, X *cmat.FloatMatrix, alpha, beta float64, bits int, confs ...*gomas.Config) *gomas.Error {
	ok := true
	yr, yc := Y.Size()
	ar, ac := A.Size()
	xr, xc := X.Size()

	if ar*ac == 0 {
		return nil
	}
	if yr != 1 && yc != 1 {
		return gomas.NewError(gomas.ENEED_VECTOR, "MVMultSym")
	}
	if xr != 1 && xc != 1 {
		return gomas.NewError(gomas.ENEED_VECTOR, "MVMultSym")
	}
	nx := X.Len()
	ny := Y.Len()

	ok = ny == ar && nx == ac && ac == ar
	if !ok {
		return gomas.NewError(gomas.ESIZE, "MVMultSym")
	}
	if beta != 1.0 {
		vscal(Y, beta, ny)
	}
	symv(Y, A, X, alpha, bits, nx)
	return nil
}

/*
 * Bidiagonal top to bottom implicit zero shift QR sweep
 */
func bdQRzero(D, E, Cr, Sr, Cl, Sl *cmat.FloatMatrix, saves bool) int {
	var d1, e1, d2, cosr, sinr, cosl, sinl, r float64
	N := D.Len()

	d1 = D.GetAtUnsafe(0)
	cosr = 1.0
	cosl = 1.0
	for k := 0; k < N-1; k++ {
		e1 = E.GetAtUnsafe(k)
		d2 = D.GetAtUnsafe(k + 1)
		cosr, sinr, r = ComputeGivens(d1*cosr, e1)
		if k > 0 {
			E.SetAtUnsafe(k-1, sinl*r)
		}
		cosl, sinl, r = ComputeGivens(cosl*r, sinr*d2)
		D.SetAtUnsafe(k, r)
		d1 = d2
		if saves {
			Cr.SetAtUnsafe(k, cosr)
			Sr.SetAtUnsafe(k, sinr)
			Cl.SetAtUnsafe(k, cosl)
			Sl.SetAtUnsafe(k, sinl)
		}
	}
	d2 = cosr * d2
	D.SetAtUnsafe(N-1, d2*cosl)
	E.SetAtUnsafe(N-2, d2*sinl)
	return N - 1
}

// Compute the updated rank-1 update vector with precomputed deltas
func trdsecUpdateVecDelta(z, Q, d *cmat.FloatMatrix, rho float64) {
	var delta cmat.FloatMatrix
	for i := 0; i < d.Len(); i++ {
		delta.Column(Q, i)
		zk := trdsecUpdateElemDelta(d, &delta, i, rho)
		z.SetAt(i, zk)
	}
}

// Compute eigenmatrix Q for updated eigenvalues in 'dl'.
func trdsecEigenBuild(Q, z, Q2 *cmat.FloatMatrix) {
	var qi, delta cmat.FloatMatrix

	for k := 0; k < z.Len(); k++ {
		qi.Column(Q, k)
		delta.Row(Q2, k)
		trdsecEigenVecDelta(&qi, &delta, z)
	}
}

func Sum(X *cmat.FloatMatrix, confs ...*gomas.Config) float64 {
	if X.Len() == 0 {
		return 0.0
	}
	xr, xc := X.Size()
	if xr != 1 && xc != 1 {
		return 0.0
	}
	return sum(X, X.Len())
}