Golang FloatMatrix.SetAt方法代碼示例

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

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

// Solve secular function arising in symmetric eigenproblems. On exit 'Y' contains new
// eigenvalues and 'V' the rank-one update vector corresponding new eigenvalues.
// The matrix Qd holds for each eigenvalue then computed deltas as row vectors.
// On entry 'D' holds original eigenvalues and 'Z' is the rank-one update vector.
func TRDSecularSolveAll(y, v, Qd, d, z *cmat.FloatMatrix, rho float64, confs ...*gomas.Config) (err *gomas.Error) {
	var delta cmat.FloatMatrix
	var lmbda float64
	var e, ei int

	ei = 0
	err = nil
	if y.Len() != d.Len() || z.Len() != d.Len() || m(Qd) != n(Qd) || m(Qd) != d.Len() {
		err = gomas.NewError(gomas.ESIZE, "TRDSecularSolveAll")
		return
	}
	for i := 0; i < d.Len(); i++ {
		delta.Row(Qd, i)
		lmbda, e = trdsecRoot(d, z, &delta, i, rho)
		if e < 0 && ei == 0 {
			ei = -(i + 1)
		}
		y.SetAt(i, lmbda)
	}
	if ei == 0 {
		trdsecUpdateVecDelta(v, Qd, d, rho)
	} else {
		err = gomas.NewError(gomas.ECONVERGE, "TRDSecularSolveAll", ei)
	}
	return
}

// 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 eigenvector corresponding precomputed deltas
func trdsecEigenVecDelta(qi, delta, z *cmat.FloatMatrix) {
	var dk, zk float64

	for k := 0; k < delta.Len(); k++ {
		zk = z.GetAt(k)
		dk = delta.GetAt(k)
		qi.SetAt(k, zk/dk)
	}
	s := blasd.Nrm2(qi)
	blasd.InvScale(qi, s)
}

func sortVec(D *cmat.FloatMatrix, updown int) {
	var cval, tmpval float64

	j := 0
	for k := 0; k < D.Len(); k++ {
		cval = D.GetAt(k)
		for j = k; j > 0; j-- {
			tmpval = D.GetAt(j - 1)
			if updown > 0 && tmpval >= cval {
				break
			}
			if updown < 0 && tmpval <= cval {
				break
			}
			D.SetAt(j, tmpval)
		}
		D.SetAt(j, cval)
	}
}

// Solve secular function arising in symmetric eigenproblems. On exit 'Y' contains new
// eigenvalues On entry 'D' holds original eigenvalues and 'Z' is the rank-one update vector.
// Parameter 'delta' is workspace needed for computation.
func TRDSecularSolve(Y, D, Z, delta *cmat.FloatMatrix, rho float64, confs ...*gomas.Config) (err *gomas.Error) {
	var lmbda float64
	var e, ei int
	ei = 0
	err = nil
	if delta.Len() != D.Len() || Y.Len() != D.Len() || Z.Len() != D.Len() {
		err = gomas.NewError(gomas.ESIZE, "TRDSecularSolve")
		return
	}
	for i := 0; i < D.Len(); i++ {
		lmbda, e = trdsecRoot(D, Z, delta, i, rho)
		if e < 0 && ei == 0 {
			ei = -(i + 1)
		}
		Y.SetAt(i, lmbda)
	}
	if ei != 0 {
		err = gomas.NewError(gomas.ECONVERGE, "TRDSecularSolve", ei)
	}
	return
}

/*
 * Generic rank update of diagonal matrix.
 *   diag(D) = diag(D) + alpha * x * y.T
 *
 * Arguments:
 *   D     N element column or row vector or N-by-N matrix
 *
 *   x, y  N element vectors
 *
 *   alpha scalar
 */
func MVUpdateDiag(D, x, y *cmat.FloatMatrix, alpha float64, confs ...*gomas.Config) *gomas.Error {
	var d *cmat.FloatMatrix
	var d0 cmat.FloatMatrix

	if !x.IsVector() || !y.IsVector() {
		return gomas.NewError(gomas.ENEED_VECTOR, "MvUpdateDiag")
	}

	d = D
	if !D.IsVector() {
		d0.Diag(D)
		d = &d0
	}

	for k := 0; k < d.Len(); k++ {
		val := d.GetAt(k)
		val += x.GetAt(k) * y.GetAt(k) * alpha
		d.SetAt(k, val)
	}
	return nil
}

/*
 * Bidiagonal bottom to top implicit QL sweep
 */
func bdQLsweep(D, E, Cr, Sr, Cl, Sl *cmat.FloatMatrix, f0, g0 float64, saves bool) int {
	var d1, e1, d2, e2, f, g, cosr, sinr, cosl, sinl, r float64
	N := D.Len()

	d1 = D.GetAtUnsafe(N - 1)
	e1 = E.GetAtUnsafe(N - 2)
	f = f0
	g = g0
	for k := N - 1; k > 0; k-- {
		d2 = D.GetAtUnsafe(k - 1)
		cosr, sinr, r = ComputeGivens(f, g)
		if k < N-1 {
			E.SetAt(k, r)
		}
		f, e1 = RotateGivens(d1, e1, cosr, sinr)
		g, d2 = RotateGivens(0.0, d2, cosr, sinr)
		cosl, sinl, r = ComputeGivens(f, g)
		d1 = r
		f, d2 = RotateGivens(e1, d2, cosl, sinl)
		if k > 1 {
			e2 = E.GetAtUnsafe(k - 2)
			g, e2 = RotateGivens(0.0, e2, cosl, sinl)
			E.SetAtUnsafe(k-2, e2)
		}
		D.SetAtUnsafe(k, d1)
		//D.SetAtUnsafe(k-1, d2)
		d1 = d2
		e1 = e2
		if saves {
			Cr.SetAtUnsafe(k-1, cosr)
			Sr.SetAtUnsafe(k-1, -sinr)
			Cl.SetAtUnsafe(k-1, cosl)
			Sl.SetAtUnsafe(k-1, -sinl)
		}
	}
	D.SetAtUnsafe(0, d1)
	E.SetAt(0, f)
	return N - 1
}

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

func setDiagonals(A *cmat.FloatMatrix, offdiag, kind int) string {
	var sD, sE cmat.FloatMatrix
	var desc string

	switch kind {
	case 2:
		desc = "middle"
	case 1:
		desc = "top   "
	default:
		desc = "bottom"
	}

	sD.Diag(A, 0)
	sE.Diag(A, offdiag)
	N := sD.Len()
	for k := 0; k < N; k++ {
		if k < N-1 {
			sE.SetAt(k, 1.0)
		}
		switch kind {
		case 2: // midheavy
			if k < N/2 {
				sD.SetAt(k, float64(k+1))
			} else {
				sD.SetAt(k, float64(N-k))
			}
			break
		case 1: // top heavy
			sD.SetAt(N-1-k, float64(k+1))
			break
		default: // bottom heavy
			sD.SetAt(k, float64(k+1))
			break
		}
	}
	return desc
}

func updateDelta(delta *cmat.FloatMatrix, eta float64) {
	for i := 0; i < delta.Len(); i++ {
		delta.SetAt(i, delta.GetAt(i)-eta)
	}
}

func computeDelta(delta, D *cmat.FloatMatrix, index int, tau float64) {
	d0 := D.GetAt(index)
	for i := 0; i < delta.Len(); i++ {
		delta.SetAt(i, D.GetAt(i)-d0-tau)
	}
}

func svdSmall(S, U, V, A, W *cmat.FloatMatrix, bits int, conf *gomas.Config) (err *gomas.Error) {
	var r cmat.FloatMatrix
	var d0, d1, e0 float64

	err = nil
	K := m(A)
	if n(A) < K {
		K = n(A)
	}
	tau := cmat.NewMatrix(K, 1)
	if m(A) >= n(A) {
		if err = QRFactor(A, tau, W, conf); err != nil {
			return
		}
	} else {
		if err = LQFactor(A, tau, W, conf); err != nil {
			return
		}
	}
	if m(A) == 1 || n(A) == 1 {
		// either tall M-by-1 or wide 1-by-N
		S.SetAt(0, math.Abs(A.Get(0, 0)))

		if bits&gomas.WANTU != 0 {
			if n(A) == 1 {
				if n(U) == n(A) {
					// U is M-by-1
					U.Copy(A)
					QRBuild(U, tau, W, n(A), conf)
				} else {
					// U is M-by-M
					eye := cmat.FloatDiagonalSource{1.0}
					U.SetFrom(&eye, cmat.SYMM)
					if err = QRMult(U, A, tau, W, gomas.RIGHT, conf); err != nil {
						return
					}
				}
			} else {
				U.Set(0, 0, -1.0)
			}
		}
		if bits&gomas.WANTV != 0 {
			if m(A) == 1 {
				if m(V) == m(A) {
					// V is 1-by-N
					V.Copy(A)
					LQBuild(V, tau, W, m(A), conf)
				} else {
					// V is N-by-N
					eye := cmat.FloatDiagonalSource{1.0}
					V.SetFrom(&eye, cmat.SYMM)
					if err = LQMult(V, A, tau, W, gomas.RIGHT, conf); err != nil {
						return
					}
				}
			} else {
				// V is 1-by-1
				V.Set(0, 0, -1.0)
			}
		}
		return
	}

	// Use bdSvd2x2 functions
	d0 = A.Get(0, 0)
	d1 = A.Get(1, 1)
	if m(A) >= n(A) {
		e0 = A.Get(0, 1)
	} else {
		e0 = A.Get(1, 0)
	}

	if bits&(gomas.WANTU|gomas.WANTV) == 0 {
		// no vectors
		smin, smax := bdSvd2x2(d0, e0, d1)
		S.SetAt(0, math.Abs(smax))
		S.SetAt(1, math.Abs(smin))
		return
	}

	// at least left or right eigenvector wanted
	smin, smax, cosl, sinl, cosr, sinr := bdSvd2x2Vec(d0, e0, d1)
	if bits&gomas.WANTU != 0 {
		// left eigenvectors
		if m(A) >= n(A) {
			if n(U) == n(A) {
				U.Copy(A)
				if err = QRBuild(U, tau, W, n(A), conf); err != nil {
					return
				}
			} else {
				// U is M-by-M
				eye := cmat.FloatDiagonalSource{1.0}
				U.SetFrom(&eye, cmat.SYMM)
				if err = QRMult(U, A, tau, W, gomas.RIGHT, conf); err != nil {
					return
				}
			}
			ApplyGivensRight(U, 0, 1, 0, m(A), cosl, sinl)
		} else {
//.........這裏部分代碼省略.........

/*
 * Tridiagonal top to bottom implicit QR sweep
 */
func trdQRsweep(D, E, Cr, Sr *cmat.FloatMatrix, f0, g0 float64, saves bool) int {
	var d0, e0, e1, d1, e0r, e0c, w0, f, g, cosr, sinr, r float64
	N := D.Len()

	d0 = D.GetAt(0)
	e0 = E.GetAt(0)
	f = f0
	g = g0
	for k := 0; k < N-1; k++ {
		d1 = D.GetAt(k + 1)
		cosr, sinr, r = ComputeGivens(f, g)
		if k > 0 {
			E.SetAt(k-1, r)
		}
		d0, e0c = RotateGivens(d0, e0, cosr, sinr)
		e0r, d1 = RotateGivens(e0, d1, cosr, sinr)
		d0, e0r = RotateGivens(d0, e0r, cosr, sinr)
		e0c, d1 = RotateGivens(e0c, d1, cosr, sinr)
		// here: e0c == e0r
		if k < N-2 {
			e1 = E.GetAt(k + 1)
			w0, e1 = RotateGivens(0.0, e1, cosr, sinr)
		}
		D.SetAt(k, d0)
		d0 = d1
		e0 = e1
		f = e0r
		g = w0
		if saves {
			Cr.SetAt(k, cosr)
			Sr.SetAt(k, sinr)
		}
	}
	D.SetAt(N-1, d0)
	E.SetAt(N-2, e0r)
	return N - 1
}