Golang FloatMatrix.GetAt方法代碼示例

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

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

/*
 * 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 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 absVecMinMax(D *cmat.FloatMatrix, minmax int) int {
	var cval, tmpval float64
	cval = math.Abs(D.GetAt(0))
	ix := 0
	for k := 1; k < D.Len(); k++ {
		tmpval = math.Abs(D.GetAt(k))
		if minmax > 0 && tmpval > cval {
			cval = tmpval
			ix = k
		} else if minmax < 0 && tmpval < cval {
			cval = tmpval
			ix = k
		}
	}
	return ix
}

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

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

// estimat singular values
func estimateSval(D, E *cmat.FloatMatrix) (smin, smax float64) {
	var ssmin, ssmax, mu, e0, e1, d1 float64

	N := D.Len()

	d1 = math.Abs(D.GetAt(0))
	e1 = math.Abs(E.GetAt(0))
	e0 = e1
	ssmax = d1
	if e1 > ssmax {
		ssmax = e1
	}
	mu = d1
	ssmin = d1
	for k := 1; k < N; k++ {
		d1 = math.Abs(D.GetAt(k))
		if k < N-1 {
			e1 = math.Abs(E.GetAt(k))
		}
		if d1 > ssmax {
			ssmax = d1
		}
		if e1 > ssmax {
			ssmax = e1
		}
		if ssmin != 0.0 {
			mu = d1 * (mu / (mu + e0))
			if mu < ssmin {
				ssmin = mu
			}
		}
		e0 = e1
	}
	smax = ssmax
	smin = ssmin / math.Sqrt(float64(N))
	return
}

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

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

func UpdateGivens(A *cmat.FloatMatrix, start int, C, S *cmat.FloatMatrix, nrot, flags int) int {
	var k, l int
	var cos, sin float64
	end := start + nrot

	if flags&gomas.BACKWARD != 0 {
		if flags&gomas.LEFT != 0 {
			end = imin(m(A), end)
			k = end
			l = nrot
			for l > 0 && k > start {
				cos = C.GetAt(l - 1)
				sin = S.GetAt(l - 1)
				if cos != 1.0 || sin != 0.0 {
					ApplyGivensLeft(A, k-1, k, 0, n(A), cos, sin)
				}
				l--
				k--
			}
		} else {
			end = imin(n(A), end)
			k = end
			l = nrot
			for l > 0 && k > start {
				cos = C.GetAt(l - 1)
				sin = S.GetAt(l - 1)
				if cos != 1.0 || sin != 0.0 {
					ApplyGivensRight(A, k-1, k, 0, m(A), cos, sin)
				}
				l--
				k--
			}
		}
	} else {
		if flags&gomas.LEFT != 0 {
			end = imin(m(A), end)
			k = start
			l = 0
			for l < nrot && k < end {
				cos = C.GetAt(l)
				sin = S.GetAt(l)
				if cos != 1.0 || sin != 0.0 {
					ApplyGivensLeft(A, k, k+1, 0, n(A), cos, sin)
				}
				l++
				k++
			}
		} else {
			end = imin(n(A), end)
			k = start
			l = 0
			for l < nrot && k < end {
				cos = C.GetAt(l)
				sin = S.GetAt(l)
				if cos != 1.0 || sin != 0.0 {
					ApplyGivensRight(A, k, k+1, 0, m(A), cos, sin)
				}
				l++
				k++
			}
		}
	}
	return nrot
}

// Initial guess for iteration
func trdsecInitialGuess(tau, taulow, tauhigh *float64,
	D, Z, delta *cmat.FloatMatrix, index int, rho float64) int {

	var d_k, d_k1, z_k, z_k1, diff, mpnt, A, B, C, F, G0, G1, Hx, dd float64
	var iN, iK, N int
	N = D.Len()
	last := false
	if index == N-1 {
		iN = N - 2
		last = true
	} else {
		iN = index
	}

	d_k = D.GetAt(iN)
	d_k1 = D.GetAt(iN + 1)
	diff = d_k1 - d_k
	mpnt = diff / 2.0
	if last {
		mpnt = rho / 2.0
	}

	// compute delta = D[i] - D[index] - mpnt
	computeDelta(delta, D, index, mpnt)
	G0, _ = rationalForward(Z, delta, 0, iN+1)
	G1, _ = rationalBackward(Z, delta, iN+1, N)

	d_k = delta.GetAt(iN)
	d_k1 = delta.GetAt(iN + 1)
	z_k = Z.GetAt(iN)
	z_k1 = Z.GetAt(iN + 1)
	// F is f(x) at initial point, 1/rho + g(y) + h(y)
	F = 1.0/rho + G0 + G1
	// Hx is h(y)
	Hx = z_k*(z_k/d_k) + z_k1*(z_k1/d_k1)
	// C is g(y) at initial point
	C = F - Hx
	if last {
		goto LastEntry
	}

	if F > 0 {
		iK = index
		A = z_k * z_k * diff
		B = C*diff + z_k*z_k + z_k1*z_k1
		*taulow = 0.0
		*tauhigh = mpnt
	} else {
		iK = index + 1
		A = -z_k1 * z_k1 * diff
		B = -C*diff + z_k*z_k + z_k1*z_k1
		*taulow = -mpnt
		*tauhigh = 0.0
	}
	B = B / 2.0
	dd = discriminant(A, B, C)
	if B > 0.0 {
		*tau = A / (B + math.Sqrt(dd))
	} else {
		*tau = (B - math.Sqrt(dd)) / C
	}
	computeDelta(delta, D, iK, *tau)
	return iK

LastEntry:
	A = -z_k1 * z_k1 * diff
	B = (-C*diff + z_k*z_k + z_k1*z_k1) / 2.0
	Hx = z_k*z_k/(diff+rho) + z_k1*z_k1/rho
	if F <= 0.0 && C <= Hx {
		*tau = rho
	} else {
		dd = discriminant(A, B, C)
		if B < 0.0 {
			*tau = A / (B - math.Sqrt(dd))
		} else {
			*tau = (B + math.Sqrt(dd)) / C
		}
	}
	if F < 0.0 {
		*taulow = mpnt
		*tauhigh = rho
	} else {
		*taulow = 0.0
		*tauhigh = mpnt
	}
	computeDelta(delta, D, N-1, *tau)
	return index - 1
}

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

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