Golang FloatMatrix.Column方法代码示例

/*
 * 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 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 mNorm1(A *cmat.FloatMatrix) float64 {
	var amax float64 = 0.0
	var col cmat.FloatMatrix
	_, acols := A.Size()
	for k := 0; k < acols; k++ {
		col.Column(A, k)
		cmax := blasd.ASum(&col)
		if cmax > amax {
			amax = cmax
		}
	}
	return amax
}

/*
 * Generates the real orthogonal matrix Q which is defined as the product of K elementary
 * reflectors of order N embedded in matrix A as returned by TRDReduce().
 *
 *   A     On entry tridiagonal reduction as returned by TRDReduce().
 *         On exit the orthogonal matrix Q.
 *
 *  tau    Scalar coefficients of elementary reflectors.
 *
 *  W      Workspace
 *
 *  K      Number of reflectors , 0 < K < N
 *
 *  flags  LOWER or UPPER
 *
 *  confs  Optional blocking configuration
 *
 * If flags has UPPER set then
 *    Q = H(K)...H(1)H(0) where 0 < K < N-1
 *
 * If flags has LOWR set then
 *    Q = H(0)H(1)...H(K) where 0 < K < N-1
 */
func TRDBuild(A, tau, W *cmat.FloatMatrix, K, flags int, confs ...*gomas.Config) *gomas.Error {
	var err *gomas.Error = nil
	var Qh, tauh cmat.FloatMatrix
	var s, d cmat.FloatMatrix

	if K > m(A)-1 {
		K = m(A) - 1
	}

	switch flags & (gomas.LOWER | gomas.UPPER) {
	case gomas.LOWER:
		// Shift Q matrix embedded in A right and fill first column
		// unit column vector
		for j := m(A) - 1; j > 0; j-- {
			s.SubMatrix(A, j, j-1, m(A)-j, 1)
			d.SubMatrix(A, j, j, m(A)-j, 1)
			blasd.Copy(&d, &s)
			A.Set(0, j, 0.0)
		}
		// zero first column and set first entry to one
		d.Column(A, 0)
		blasd.Scale(&d, 0.0)
		d.Set(0, 0, 1.0)

		Qh.SubMatrix(A, 1, 1, m(A)-1, m(A)-1)
		tauh.SubMatrix(tau, 0, 0, m(A)-1, 1)
		err = QRBuild(&Qh, &tauh, W, K, confs...)

	case gomas.UPPER:
		// Shift Q matrix embedded in A left and fill last column
		// unit column vector
		for j := 1; j < m(A); j++ {
			s.SubMatrix(A, 0, j, j, 1)
			d.SubMatrix(A, 0, j-1, j, 1)
			blasd.Copy(&d, &s)
			A.Set(-1, j-1, 0.0)
		}
		// zero last column and set last entry to one
		d.Column(A, m(A)-1)
		blasd.Scale(&d, 0.0)
		d.Set(-1, 0, 1.0)

		Qh.SubMatrix(A, 0, 0, m(A)-1, m(A)-1)
		tauh.SubMatrix(tau, 0, 0, m(A)-1, 1)
		err = QLBuild(&Qh, &tauh, W, K, confs...)
	}
	if err != nil {
		err.Update("TRDBuild")
	}
	return err
}

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

func trdsecEigenBuildInplace(Q, z *cmat.FloatMatrix) {
	var QTL, QBR, Q00, q11, q12, q21, Q22, qi cmat.FloatMatrix
	var zk0, zk1, dk0, dk1 float64

	util.Partition2x2(
		&QTL, nil,
		nil, &QBR /**/, Q, 0, 0, util.PTOPLEFT)

	for m(&QBR) > 0 {
		util.Repartition2x2to3x3(&QTL,
			&Q00, nil, nil,
			nil, &q11, &q12,
			nil, &q21, &Q22 /**/, Q, 1, util.PBOTTOMRIGHT)
		//---------------------------------------------------------------
		k := m(&Q00)
		zk0 = z.GetAt(k)
		dk0 = q11.Get(0, 0)
		q11.Set(0, 0, zk0/dk0)

		for i := 0; i < q12.Len(); i++ {
			zk1 = z.GetAt(k + i + 1)
			dk0 = q12.GetAt(i)
			dk1 = q21.GetAt(i)
			q12.SetAt(i, zk0/dk1)
			q21.SetAt(i, zk1/dk0)
		}
		//---------------------------------------------------------------
		util.Continue3x3to2x2(
			&QTL, nil,
			nil, &QBR /**/, &Q00, &q11, &Q22 /**/, Q, util.PBOTTOMRIGHT)
	}
	// scale column eigenvectors
	for k := 0; k < z.Len(); k++ {
		qi.Column(Q, k)
		t := blasd.Nrm2(&qi)
		blasd.InvScale(&qi, t)
	}
}

/*
 * Generate one of the orthogonal matrices Q or P.T determined by BDReduce() when
 * reducing a real matrix A to bidiagonal form. Q and P.T are defined as products
 * elementary reflectors H(i) or G(i) respectively.
 *
 * Orthogonal matrix Q is generated if flag WANTQ is set. And matrix P respectively
 * if flag WANTP is set.
 */
func BDBuild(A, tau, W *cmat.FloatMatrix, K, flags int, confs ...*gomas.Config) *gomas.Error {
	var Qh, Ph, tauh, d, s cmat.FloatMatrix
	var err *gomas.Error = nil

	if m(A) == 0 || n(A) == 0 {
		return nil
	}

	if m(A) > n(A) || (m(A) == n(A) && flags&gomas.LOWER == 0) {
		switch flags & (gomas.WANTQ | gomas.WANTP) {
		case gomas.WANTQ:
			tauh.SubMatrix(tau, 0, 0, n(A), 1)
			err = QRBuild(A, &tauh, W, K, confs...)

		case gomas.WANTP:
			// Shift P matrix embedded in A down and fill first column and row
			// to unit vector
			for j := n(A) - 1; j > 0; j-- {
				s.SubMatrix(A, j-1, j, 1, n(A)-j)
				d.SubMatrix(A, j, j, 1, n(A)-j)
				blasd.Copy(&d, &s)
				A.Set(j, 0, 0.0)
			}
			// zero  first row and set first entry to one
			d.Row(A, 0)
			blasd.Scale(&d, 0.0)
			d.Set(0, 0, 1.0)

			Ph.SubMatrix(A, 1, 1, n(A)-1, n(A)-1)
			tauh.SubMatrix(tau, 0, 0, n(A)-1, 1)
			if K > n(A)-1 {
				K = n(A) - 1
			}
			err = LQBuild(&Ph, &tauh, W, K, confs...)
		}
	} else {
		switch flags & (gomas.WANTQ | gomas.WANTP) {
		case gomas.WANTQ:
			// Shift Q matrix embedded in A right and fill first column and row
			// to unit vector
			for j := m(A) - 1; j > 0; j-- {
				s.SubMatrix(A, j, j-1, m(A)-j, 1)
				d.SubMatrix(A, j, j, m(A)-j, 1)
				blasd.Copy(&d, &s)
				A.Set(0, j, 0.0)
			}
			// zero first column and set first entry to one
			d.Column(A, 0)
			blasd.Scale(&d, 0.0)
			d.Set(0, 0, 1.0)

			Qh.SubMatrix(A, 1, 1, m(A)-1, m(A)-1)
			tauh.SubMatrix(tau, 0, 0, m(A)-1, 1)
			if K > m(A)-1 {
				K = m(A) - 1
			}
			err = QRBuild(&Qh, &tauh, W, K, confs...)

		case gomas.WANTP:
			tauh.SubMatrix(tau, 0, 0, m(A), 1)
			err = LQBuild(A, &tauh, W, K, confs...)
		}
	}
	if err != nil {
		err.Update("BDBuild")
	}
	return err
}

// unblocked LU decomposition with pivots: FLAME LU variant 3; Left-looking
func unblockedLUpiv(A *cmat.FloatMatrix, p *Pivots, offset int, conf *gomas.Config) *gomas.Error {
	var err *gomas.Error = nil
	var ATL, ATR, ABL, ABR cmat.FloatMatrix
	var A00, a01, A02, a10, a11, a12, A20, a21, A22 cmat.FloatMatrix
	var AL, AR, A0, a1, A2, aB1, AB0 cmat.FloatMatrix
	var pT, pB, p0, p1, p2 Pivots

	err = nil
	util.Partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, 0, 0, util.PTOPLEFT)
	util.Partition1x2(
		&AL, &AR, A, 0, util.PLEFT)
	partitionPivot2x1(
		&pT,
		&pB, *p, 0, util.PTOP)

	for m(&ATL) < m(A) && n(&ATL) < n(A) {
		util.Repartition2x2to3x3(&ATL,
			&A00, &a01, &A02,
			&a10, &a11, &a12,
			&A20, &a21, &A22 /**/, A, 1, util.PBOTTOMRIGHT)
		util.Repartition1x2to1x3(&AL,
			&A0, &a1, &A2 /**/, A, 1, util.PRIGHT)
		repartPivot2x1to3x1(&pT,
			&p0, &p1, &p2 /**/, *p, 1, util.PBOTTOM)

		// apply previously computed pivots on current column
		applyPivots(&a1, p0)

		// a01 = trilu(A00) \ a01 (TRSV)
		blasd.MVSolveTrm(&a01, &A00, 1.0, gomas.LOWER|gomas.UNIT)
		// a11 = a11 - a10 *a01
		aval := a11.Get(0, 0) - blasd.Dot(&a10, &a01)
		a11.Set(0, 0, aval)
		// a21 = a21 -A20*a01
		blasd.MVMult(&a21, &A20, &a01, -1.0, 1.0, gomas.NONE)

		// pivot index on current column [a11, a21].T
		aB1.Column(&ABR, 0)
		p1[0] = pivotIndex(&aB1)
		// pivots to current column
		applyPivots(&aB1, p1)

		// a21 = a21 / a11
		if aval == 0.0 {
			if err == nil {
				ij := m(&ATL) + p1[0] - 1
				err = gomas.NewError(gomas.ESINGULAR, "DecomposeLU", ij)
			}
		} else {
			blasd.InvScale(&a21, a11.Get(0, 0))
		}

		// apply pivots to previous columns
		AB0.SubMatrix(&ABL, 0, 0)
		applyPivots(&AB0, p1)
		// scale last pivots to origin matrix row numbers
		p1[0] += m(&ATL)

		util.Continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &a11, &A22, A, util.PBOTTOMRIGHT)
		util.Continue1x3to1x2(
			&AL, &AR, &A0, &a1, A, util.PRIGHT)
		contPivot3x1to2x1(
			&pT,
			&pB, p0, p1, *p, util.PBOTTOM)
	}
	if n(&ATL) < n(A) {
		applyPivots(&ATR, *p)
		blasd.SolveTrm(&ATR, &ATL, 1.0, gomas.LEFT|gomas.UNIT|gomas.LOWER, conf)
	}
	return err
}