Golang FloatMatrix.Copy方法代碼示例

/*
 * Build full block reflect T for nc columns from sequence of reflector stored in S.
 * Reflectors in S are the diagonal of T, off-diagonal values of reflector are computed
 * from elementary reflector store in lower triangular part of A.
 */
func buildQRTReflector(T, A, S *cmat.FloatMatrix, nc int, conf *gomas.Config) *gomas.Error {
	var ATL, ATR, ABL, ABR cmat.FloatMatrix
	var A00, A10, A11, A20, A21, A22 cmat.FloatMatrix
	var TTL, TTR, TBL, TBR cmat.FloatMatrix
	var T00, T01, T02, T11, T12, T22 cmat.FloatMatrix
	var SL, SR cmat.FloatMatrix
	var S00, S01, S02 cmat.FloatMatrix

	util.Partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, 0, 0, util.PTOPLEFT)
	util.Partition2x2(
		&TTL, &TTR,
		&TBL, &TBR, T, 0, 0, util.PTOPLEFT)
	util.Partition1x2(
		&SL, &SR, S, 0, util.PLEFT)

	nb := conf.LB
	for m(&ABR)-nb > 0 && n(&ABR)-nb > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, nil, nil,
			&A10, &A11, nil,
			&A20, &A21, &A22, A, nb, util.PBOTTOMRIGHT)
		util.Repartition2x2to3x3(&TTL,
			&T00, &T01, &T02,
			nil, &T11, &T12,
			nil, nil, &T22, T, nb, util.PBOTTOMRIGHT)
		util.Repartition1x2to1x3(&SL,
			&S00, &S01, &S02, S, nb, util.PRIGHT)
		// --------------------------------------------------------
		// update T01: T01 = -T00*Y1.T*Y2*T11
		//  Y1 = /A10\   Y2 = /A11\
		//       \A20/        \A21/
		//
		T11.Copy(&S01)
		updateQRTReflector(&T01, &A10, &A20, &A11, &A21, &T00, &S01, conf)
		// --------------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &A11, &A22, A, util.PBOTTOMRIGHT)
		util.Continue3x3to2x2(
			&TTL, &TTR,
			&TBL, &TBR, &T00, &T11, &T22, T, util.PBOTTOMRIGHT)
		util.Continue1x3to1x2(
			&SL, &SR, &S00, &S01, S, util.PRIGHT)
	}
	if m(&ABR) > 0 && n(&ABR) > 0 {
	}
	return nil
}

func testEigen(N int, bits int, t *testing.T) {
	var A, A0, W, D, V *cmat.FloatMatrix
	var sD cmat.FloatMatrix
	var s string = "lower"

	if bits&gomas.UPPER != 0 {
		s = "upper"
	}

	wsize := N * N
	if wsize < 100 {
		wsize = 100
	}

	D = cmat.NewMatrix(N, 1)
	A = cmat.NewMatrix(N, N)
	V = cmat.NewMatrix(N, N)

	src := cmat.NewFloatNormSource()
	A.SetFrom(src, cmat.SYMM)
	A0 = cmat.NewCopy(A)
	W = cmat.NewMatrix(wsize, 1)

	if err := lapackd.EigenSym(D, A, W, bits|gomas.WANTV); err != nil {
		t.Errorf("EigenSym error: %v\n", err)
		return
	}

	// ||I - V.T*V||
	sD.Diag(V)
	blasd.Mult(V, A, A, 1.0, 0.0, gomas.TRANSA)
	blasd.Add(&sD, -1.0)
	nrm1 := lapackd.NormP(V, lapackd.NORM_ONE)

	// left vectors are M-by-N
	V.Copy(A)
	lapackd.MultDiag(V, D, gomas.RIGHT)
	blasd.Mult(A0, V, A, -1.0, 1.0, gomas.TRANSB)
	nrm2 := lapackd.NormP(A0, lapackd.NORM_ONE)

	t.Logf("N=%d, [%s] ||A - V*D*V.T||_1 :%e\n", N, s, nrm2)
	t.Logf("  ||I - V.T*V||_1 : %e\n", nrm1)
}

func TestSubMatrixCopy(t *testing.T) {
	var subA, subB cmat.FloatMatrix
	M := 9
	N := 9
	A := cmat.NewMatrix(M, N)
	B := cmat.NewMatrix(M, N)
	twos := cmat.NewFloatConstSource(2.0)
	B.SetFrom(twos)
	subA.SubMatrix(A, 1, 1, M-2, N-2)
	subB.SubMatrix(B, 1, 1, M-2, N-2)
	subA.Copy(&subB)
	ok := subA.AllClose(&subB)
	if !ok {
		t.Logf("copy status: %v\n", ok)
		if N < 9 {
			t.Logf("subA\n%v\n", subA)
		}
	}
	if N < 10 {
		t.Logf("A\n%v\n", A)
	}
}

/*
 * Find diagonal pivot and build incrementaly updated block.
 *
 *  (AL)  (AR)                   (WL)  (WR)
 *  --------------------------   ----------    k'th row in W
 *  x x | c1                     w w | k kp1
 *  x x | c1 d                   w w | k kp1
 *  x x | c1 x  d                w w | k kp1
 *  x x | c1 x  x  d             w w | k kp1
 *  x x | c1 r2 r2 r2 r2         w w | k kp1
 *  x x | c1 x  x  x  r2 d       w w | k kp1
 *  x x | c1 x  x  x  r2 x d     w w | k kp1
 *
 * Matrix AR contains the unfactored part of the matrix and AL the already
 * factored columns. Matrix WL is updated values of factored part ie.
 * w(i) = l(i)d(i). Matrix WR will have updated values for next column.
 * Column WR(k) contains updated AR(c1) and WR(kp1) possible pivot row AR(r2).
 */
func findAndBuildBKPivotLower(AL, AR, WL, WR *cmat.FloatMatrix, k int) (int, int) {
	var r, q int
	var rcol, qrow, src, wk, wkp1, wrow cmat.FloatMatrix

	// Copy AR column 0 to WR column 0 and update with WL[0:]
	src.SubMatrix(AR, 0, 0, m(AR), 1)
	wk.SubMatrix(WR, 0, 0, m(AR), 1)
	wk.Copy(&src)
	if k > 0 {
		wrow.SubMatrix(WL, 0, 0, 1, n(WL))
		blasd.MVMult(&wk, AL, &wrow, -1.0, 1.0, gomas.NONE)
	}
	if m(AR) == 1 {
		return 0, 1
	}
	amax := math.Abs(WR.Get(0, 0))

	// find max off-diagonal on first column.
	rcol.SubMatrix(WR, 1, 0, m(AR)-1, 1)

	// r is row index and rmax is its absolute value
	r = blasd.IAmax(&rcol) + 1
	rmax := math.Abs(rcol.Get(r-1, 0))
	if amax >= bkALPHA*rmax {
		// no pivoting, 1x1 diagonal
		return 0, 1
	}
	// Now we need to copy row r to WR[:,1] and update it
	wkp1.SubMatrix(WR, 0, 1, m(AR), 1)
	qrow.SubMatrix(AR, r, 0, 1, r+1)
	blasd.Copy(&wkp1, &qrow)
	if r < m(AR)-1 {
		var wkr cmat.FloatMatrix
		qrow.SubMatrix(AR, r, r, m(AR)-r, 1)
		wkr.SubMatrix(&wkp1, r, 0, m(&wkp1)-r, 1)
		blasd.Copy(&wkr, &qrow)
	}
	if k > 0 {
		// update wkp1
		wrow.SubMatrix(WL, r, 0, 1, n(WL))
		blasd.MVMult(&wkp1, AL, &wrow, -1.0, 1.0, gomas.NONE)
	}

	// set on-diagonal entry to zero to avoid finding it
	p1 := wkp1.Get(r, 0)
	wkp1.Set(r, 0, 0.0)
	// max off-diagonal on r'th column/row at index q
	q = blasd.IAmax(&wkp1)
	qmax := math.Abs(wkp1.Get(q, 0))
	// restore on-diagonal entry
	wkp1.Set(r, 0, p1)

	if amax >= bkALPHA*rmax*(rmax/qmax) {
		// no pivoting, 1x1 diagonal
		return 0, 1
	}
	// if q == r then qmax is not off-diagonal, qmax == WR[r,1] and
	// we get 1x1 pivot as following is always true
	if math.Abs(WR.Get(r, 1)) >= bkALPHA*qmax {
		// 1x1 pivoting and interchange with k, r
		// pivot row in column WR[:,1] to W[:,0]
		src.SubMatrix(WR, 0, 1, m(AR), 1)
		wkp1.SubMatrix(WR, 0, 0, m(AR), 1)
		blasd.Copy(&wkp1, &src)
		wkp1.Set(0, 0, src.Get(r, 0))
		wkp1.Set(r, 0, src.Get(0, 0))
		return r, 1
	} else {
		// 2x2 pivoting and interchange with k+1, r
		return r, 2
	}
	return 0, 1
}

func svdWide(S, U, V, A, W *cmat.FloatMatrix, bits int, conf *gomas.Config) (err *gomas.Error) {
	var uu, vv *cmat.FloatMatrix
	var tauq, taup, Wred, sD, sE, L, Vm cmat.FloatMatrix

	if (bits & (gomas.WANTU | gomas.WANTV)) != 0 {
		if W.Len() < 4*n(A) {
			err = gomas.NewError(gomas.ESIZE, "SVD")
			return
		}
	}
	tauq.SetBuf(m(A)-1, 1, m(A)-1, W.Data())
	taup.SetBuf(m(A), 1, m(A), W.Data()[tauq.Len():])
	wrl := W.Len() - 2*m(A) - 1
	Wred.SetBuf(wrl, 1, wrl, W.Data()[2*m(A)-1:])

	if svdCrossover(n(A), m(A)) {
		goto do_n_much_bigger
	}

	// reduce to bidiagonal form
	if err = BDReduce(A, &tauq, &taup, &Wred, conf); err != nil {
		return
	}

	sD.Diag(A)
	sE.Diag(A, -1)
	blasd.Copy(S, &sD)

	// leftt vectors
	if bits&gomas.WANTU != 0 {
		L.SubMatrix(A, 0, 0, m(A), m(A))
		U.Copy(&L)
		cmat.TriL(U, 0)
		if err = BDBuild(U, &tauq, &Wred, m(U), gomas.WANTQ|gomas.LOWER, conf); err != nil {
			return
		}
		uu = U
	}
	// right vectors
	if bits&gomas.WANTV != 0 {
		if m(V) == m(A) {
			// V is M-by-N; copy and make upper triangular
			V.Copy(A)
			//cmat.TriU(V, 0)
			if err = BDBuild(V, &taup, &Wred, m(V), gomas.WANTP, conf); err != nil {
				return
			}
		} else {
			// V is N-by-N
			eye := cmat.FloatDiagonalSource{1.0}
			V.SetFrom(&eye, cmat.SYMM)
			err = BDMult(V, A, &taup, &Wred, gomas.MULTP|gomas.LEFT|gomas.TRANS, conf)
			if err != nil {
				return
			}
		}
		vv = V
	}
	err = BDSvd(S, &sE, uu, vv, W, bits|gomas.LOWER)
	return

do_n_much_bigger:
	// here N >> M, use LQ factor first
	if err = LQFactor(A, &taup, &Wred, conf); err != nil {
		return
	}
	if bits&gomas.WANTV != 0 {
		if m(V) == m(A) {
			V.Copy(A)
			if err = LQBuild(V, &taup, &Wred, m(A), conf); err != nil {
				return
			}
		} else {
			// V is N-by-N
			eye := cmat.FloatDiagonalSource{1.0}
			V.SetFrom(&eye, cmat.SYMM)
			if err = LQMult(V, A, &taup, &Wred, gomas.RIGHT, conf); err != nil {
				return
			}
		}
	}
	L.SubMatrix(A, 0, 0, m(A), m(A))
	cmat.TriL(&L, 0)

	// resize tauq/taup for UPPER bidiagonal reduction
	tauq.SetBuf(m(A), 1, m(A), W.Data())
	taup.SetBuf(m(A)-1, 1, m(A)-1, W.Data()[tauq.Len():])

	// bidiagonal reduce
	if err = BDReduce(&L, &tauq, &taup, &Wred, conf); err != nil {
		return
	}

	if bits&gomas.WANTV != 0 {
		Vm.SubMatrix(V, 0, 0, m(A), n(A))
		err = BDMult(&Vm, &L, &taup, &Wred, gomas.MULTP|gomas.LEFT|gomas.TRANS, conf)
		if err != nil {
			return
		}
		vv = V
//.........這裏部分代碼省略.........

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 {
//.........這裏部分代碼省略.........

// Compute SVD when m(A) >= n(A)
func svdTall(S, U, V, A, W *cmat.FloatMatrix, bits int, conf *gomas.Config) (err *gomas.Error) {
	var uu, vv *cmat.FloatMatrix
	var tauq, taup, Wred, sD, sE, R, Un cmat.FloatMatrix

	if (bits & (gomas.WANTU | gomas.WANTV)) != 0 {
		if W.Len() < 4*n(A) {
			err = gomas.NewError(gomas.ESIZE, "SVD")
			return
		}
	}
	tauq.SetBuf(n(A), 1, n(A), W.Data())
	taup.SetBuf(n(A)-1, 1, n(A)-1, W.Data()[tauq.Len():])
	wrl := W.Len() - 2*n(A) - 1
	Wred.SetBuf(wrl, 1, wrl, W.Data()[2*n(A)-1:])

	if svdCrossover(m(A), n(A)) {
		goto do_m_much_bigger
	}

	// reduce to bidiagonal form
	if err = BDReduce(A, &tauq, &taup, &Wred, conf); err != nil {
		return
	}

	sD.Diag(A)
	sE.Diag(A, 1)
	blasd.Copy(S, &sD)

	// left vectors
	if bits&gomas.WANTU != 0 {
		if n(U) == n(A) {
			// U is M-by-N; copy and make lower triangular
			U.Copy(A)
			cmat.TriL(U, 0)
			if err = BDBuild(U, &tauq, &Wred, n(U), gomas.WANTQ, conf); err != nil {
				return
			}
		} else {
			// U is M-by-M
			eye := cmat.FloatDiagonalSource{1.0}
			U.SetFrom(&eye, cmat.SYMM)
			if err = BDMult(U, A, &tauq, &Wred, gomas.MULTQ|gomas.RIGHT, conf); err != nil {
				return
			}
		}
		uu = U
	}
	// right vectors
	if bits&gomas.WANTV != 0 {
		R.SubMatrix(A, 0, 0, n(A), n(A))
		V.Copy(&R)
		cmat.TriU(V, 0)
		if err = BDBuild(V, &taup, &Wred, m(V), gomas.WANTP, conf); err != nil {
			return
		}
		vv = V
	}
	err = BDSvd(S, &sE, uu, vv, W, bits|gomas.UPPER)
	return

do_m_much_bigger:
	// M >> N here; first use QR factorization
	if err = QRFactor(A, &tauq, &Wred, conf); err != nil {
		return
	}
	if bits&gomas.WANTU != 0 {
		if n(U) == n(A) {
			U.Copy(A)
			if err = QRBuild(U, &tauq, &Wred, n(U), 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, &tauq, &Wred, gomas.LEFT, conf); err != nil {
				return
			}
		}
	}
	R.SubMatrix(A, 0, 0, n(A), n(A))
	cmat.TriU(&R, 0)

	// bidiagonal reduce
	if err = BDReduce(&R, &tauq, &taup, &Wred, conf); err != nil {
		return
	}

	if bits&gomas.WANTU != 0 {
		Un.SubMatrix(U, 0, 0, m(A), n(A))
		if err = BDMult(&Un, &R, &tauq, &Wred, gomas.MULTQ|gomas.RIGHT, conf); err != nil {
			return
		}
		uu = U
	}
	if bits&gomas.WANTV != 0 {
		V.Copy(&R)
		if err = BDBuild(V, &taup, &Wred, m(V), gomas.WANTP, conf); err != nil {
			return
//.........這裏部分代碼省略.........