Golang matrix.FloatZeros函数代码示例

func _TestMultMVTransA(t *testing.T) {
	bM := 1000 * M
	bN := 1000 * N
	A := matrix.FloatNormal(bN, bM)
	X := matrix.FloatWithValue(bN, 1, 1.0)
	Y1 := matrix.FloatZeros(bM, 1)
	Y0 := matrix.FloatZeros(bM, 1)

	Ar := A.FloatArray()
	Xr := X.FloatArray()
	Y1r := Y1.FloatArray()

	blas.GemvFloat(A, X, Y0, 1.0, 1.0, linalg.OptTrans)

	DMultMV(Y1r, Ar, Xr, 1.0, 1.0, TRANSA, 1, A.LeadingIndex(), 1, 0, bN, 0, bM, 4, 4)
	ok := Y0.AllClose(Y1)
	t.Logf("Y0 == Y1: %v\n", ok)
	if !ok {
		var y1, y0 matrix.FloatMatrix
		Y1.SubMatrix(&y1, 0, 0, 5, 1)
		t.Logf("Y1[0:5]:\n%v\n", y1)
		Y0.SubMatrix(&y0, 0, 0, 5, 1)
		t.Logf("Y0[0:5]:\n%v\n", y0)
	}
}

func (gp *gpConvexProg) F1(x *matrix.FloatMatrix) (f, Df *matrix.FloatMatrix, err error) {
    f = nil
    Df = nil
    err = nil
    f = matrix.FloatZeros(gp.mnl+1, 1)
    Df = matrix.FloatZeros(gp.mnl+1, gp.n)
    y := gp.g.Copy()
    blas.GemvFloat(gp.F, x, y, 1.0, 1.0)

    for i, s := range gp.ind {
        start := s[0]
        stop := s[1]
        // yi := exp(yi) = exp(Fi*x+gi)
        ymax := maxvec(y.FloatArray()[start:stop])
        // ynew = exp(y[start:stop] - ymax)
        ynew := matrix.Exp(matrix.FloatVector(y.FloatArray()[start:stop]).Add(-ymax))
        y.SetIndexesFromArray(ynew.FloatArray(), matrix.Indexes(start, stop)...)

        // fi = log sum yi = log sum exp(Fi*x+gi)
        ysum := blas.AsumFloat(y, &la.IOpt{"n", stop - start}, &la.IOpt{"offset", start})
        f.SetIndex(i, ymax+math.Log(ysum))

        blas.ScalFloat(y, 1.0/ysum, &la.IOpt{"n", stop - start}, &la.IOpt{"offset", start})
        blas.GemvFloat(gp.F, y, Df, 1.0, 0.0, la.OptTrans, &la.IOpt{"m", stop - start},
            &la.IOpt{"incy", gp.mnl + 1}, &la.IOpt{"offseta", start},
            &la.IOpt{"offsetx", start}, &la.IOpt{"offsety", i})
    }
    return
}

func _TestMultSymmLowerSmall(t *testing.T) {
	//bM := 5
	bN := 7
	bP := 7
	Adata := [][]float64{
		[]float64{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
		[]float64{1.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0},
		[]float64{1.0, 2.0, 3.0, 0.0, 0.0, 0.0, 0.0},
		[]float64{1.0, 2.0, 3.0, 4.0, 0.0, 0.0, 0.0},
		[]float64{1.0, 2.0, 3.0, 4.0, 5.0, 0.0, 0.0},
		[]float64{1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 0.0},
		[]float64{1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0}}

	A := matrix.FloatMatrixFromTable(Adata, matrix.RowOrder)
	B := matrix.FloatNormal(bN, bP)
	C0 := matrix.FloatZeros(bN, bP)
	C1 := matrix.FloatZeros(bN, bP)

	Ar := A.FloatArray()
	Br := B.FloatArray()
	C1r := C1.FloatArray()

	blas.SymmFloat(A, B, C0, 1.0, 1.0, linalg.OptLower, linalg.OptRight)

	DMultSymm(C1r, Ar, Br, 1.0, 1.0, LOWER|RIGHT, bN, A.LeadingIndex(), bN,
		bN, 0, bP, 0, bN, 2, 2, 2)
	ok := C0.AllClose(C1)
	t.Logf("C0 == C1: %v\n", ok)
	if !ok {
		t.Logf("A=\n%v\n", A)
		t.Logf("blas: C=A*B\n%v\n", C0)
		t.Logf("C1: C1 = A*X\n%v\n", C1)
	}
}

func TestQRSmal(t *testing.T) {
	data := [][]float64{
		[]float64{12.0, -51.0, 4.0},
		[]float64{6.0, 167.0, -68.0},
		[]float64{-4.0, 24.0, -41.0}}

	A := matrix.FloatMatrixFromTable(data, matrix.RowOrder)
	T := matrix.FloatZeros(A.Cols(), A.Cols())
	T0 := T.Copy()

	M := A.Rows()
	//N := A.Cols()
	Tau := matrix.FloatZeros(M, 1)
	X, _ := DecomposeQR(A.Copy(), Tau, nil, 0)
	t.Logf("A\n%v\n", A)
	t.Logf("X\n%v\n", X)
	t.Logf("Tau\n%v\n", Tau)

	Tau0 := matrix.FloatZeros(M, 1)
	lapack.Geqrf(A, Tau0)
	t.Logf("lapack X\n%v\n", A)
	t.Logf("lapack Tau\n%v\n", Tau0)

	unblkQRBlockReflector(X, Tau, T)
	t.Logf("T:\n%v\n", T)

	V := TriLU(X.Copy())
	lapack.LarftFloat(V, Tau, T0)
	t.Logf("T0:\n%v\n", T0)

}

func _TestMultMV(t *testing.T) {
	bM := 100 * M
	bN := 100 * N
	A := matrix.FloatNormal(bM, bN)
	X := matrix.FloatNormal(bN, 1)
	Y1 := matrix.FloatZeros(bM, 1)
	Y0 := matrix.FloatZeros(bM, 1)

	Ar := A.FloatArray()
	Xr := X.FloatArray()
	Y1r := Y1.FloatArray()

	blas.GemvFloat(A, X, Y0, 1.0, 1.0)

	DMultMV(Y1r, Ar, Xr, 1.0, 1.0, NOTRANS, 1, A.LeadingIndex(), 1, 0, bN, 0, bM, 32, 32)
	t.Logf("Y0 == Y1: %v\n", Y0.AllClose(Y1))
	/*
	   if ! Y0.AllClose(Y1) {
	       y0 := Y0.SubMatrix(0, 0, 2, 1)
	       y1 := Y1.SubMatrix(0, 0, 2, 1)
	       t.Logf("y0=\n%v\n", y0)
	       t.Logf("y1=\n%v\n", y1)
	   }
	*/
}

func TestMultQT(t *testing.T) {
	M := 60
	N := 40
	K := 30
	nb := 12
	A := matrix.FloatUniform(M, N)
	B := matrix.FloatUniform(M, K)
	W := matrix.FloatZeros(N, nb)
	T := matrix.FloatZeros(N, N)

	// QR = Q*R
	QR, err := DecomposeQRT(A.Copy(), T, W, nb)
	if err != nil {
		t.Logf("decompose-err: %v\n", err)
	}

	// compute: B - Q*Q.T*B = 0

	// X = Q*Q.T*B
	X := B.Copy()
	MultQT(X, QR, T, W, LEFT|TRANS, nb)
	MultQT(X, QR, T, W, LEFT, nb)
	B.Minus(X)

	// ||B - Q*Q.T*B||_1
	nrm := NormP(B, NORM_ONE)
	t.Logf("||B - Q*Q.T*B||_1: %e\n", nrm)
}

func TestUpdateTrmMV(t *testing.T) {
	//bM := 5
	bN := 8
	//bP := 4
	nb := 4
	X := matrix.FloatNormal(bN, 1)
	//B := matrix.FloatNormal(bP, bN)
	Y := X.Copy()
	C0 := matrix.FloatZeros(bN, bN)
	C2 := matrix.FloatZeros(bN, bN)
	C1 := matrix.FloatZeros(bN, bN)

	Xr := X.FloatArray()
	Yr := Y.FloatArray()
	C1r := C1.FloatArray()
	C0r := C0.FloatArray()
	C2r := C2.FloatArray()

	// no transpose
	DRankMV(C1r, Xr, Yr, 1.0, C1.LeadingIndex(), 1, 1,
		0, bN, 0, bN, nb, nb)
	DTrmUpdMV(C0r, Xr, Yr, 1.0, LOWER, C0.LeadingIndex(), 1, 1,
		0, bN, nb)
	DTrmUpdMV(C2r, Xr, Yr, 1.0, UPPER, C2.LeadingIndex(), 1, 1,
		0, bN, nb)

	t.Logf("C1:\n%v\nC0:\n%v\nC2:\n%v\n", C1, C0, C2)
	// C0 == C2.T
	t.Logf("C0 == C2.T: %v\n", C0.AllClose(C2.Transpose()))
	// C1 == C1 - C2 + C0.T
	Cn := matrix.Minus(C1, C2)
	Cn.Plus(C0.Transpose())
	t.Logf("C1 == C1 - C2 + C0.T: %v\n", Cn.AllClose(C1))

}

func TestSolveLeastSquaresQRT(t *testing.T) {
	M := 60
	N := 40
	K := 30
	nb := 12

	A := matrix.FloatUniform(M, N)
	B := matrix.FloatZeros(M, K)
	X0 := matrix.FloatUniform(N, K)

	// B = A*X0
	Mult(B, A, X0, 1.0, 0.0, NOTRANS)
	W := matrix.FloatZeros(N, nb)
	T := matrix.FloatZeros(N, N)

	QR, err := DecomposeQRT(A, T, W, nb)
	if err != nil {
		t.Logf("decompose error: %v\n", err)
	}
	// B' = A.-1*B
	err = SolveQRT(B, QR, T, W, NOTRANS, nb)

	// expect B[0:N, 0:K] == X0, B[N:M, 0:K] == 0.0
	var Xref matrix.FloatMatrix
	Bref := matrix.FloatZeros(M, K)
	Bref.SubMatrix(&Xref, 0, 0, N, K)
	Xref.Plus(X0)
	Bref.Minus(B)
	t.Logf("\nmin ||B - A*X0||\n\twhere B = A*X0\n")
	t.Logf("||B - A*X0||_1 ~~ 0.0: %e\n", NormP(Bref, NORM_ONE))
}

func runTest(A *matrix.FloatMatrix, ntest, LB int) time.Duration {
	var W *matrix.FloatMatrix = nil
	var mintime time.Duration

	N := A.Cols()
	tau := matrix.FloatZeros(N, 1)
	if LB > 0 {
		W = matrix.FloatZeros(A.Rows(), LB)
	}
	fnc := func() {
		_, ERRmatops = matops.DecomposeQR(A, tau, W, LB)
	}

	A0 := A.Copy()
	for n := 0; n < ntest; n++ {
		if n > 0 {
			// restore original A
			A0.CopyTo(A)
			tau.Scale(0.0)
		}
		mperf.FlushCache()
		time0 := mperf.Timeit(fnc)
		if n == 0 || time0 < mintime {
			mintime = time0
		}
		if verbose {
			fmt.Printf("%.4f ms\n", time0.Seconds()*1000.0)
		}
	}
	return mintime
}

//
// Solves a geometric program
//
//   minimize    log sum exp (F0*x+g0)
//   subject to  log sum exp (Fi*x+gi) <= 0,  i=1,...,m
//               G*x <= h      
//               A*x = b
//
func Gp(K []int, F, g, G, h, A, b *matrix.FloatMatrix, solopts *SolverOptions) (sol *Solution, err error) {

    if err = checkArgK(K); err != nil {
        return
    }
    l := sumdim(K)

    if F == nil || F.Rows() != l {
        err = errors.New(fmt.Sprintf("'F' must matrix with %d rows", l))
        return
    }

    if g == nil || !g.SizeMatch(l, 1) {
        err = errors.New(fmt.Sprintf("'g' must matrix with size (%d,1)", l))
        return
    }
    n := F.Cols()

    if G == nil {
        G = matrix.FloatZeros(0, n)
    }
    if h == nil {
        h = matrix.FloatZeros(0, 1)
    }
    if G.Cols() != n {
        err = errors.New(fmt.Sprintf("'G' must matrix with size %d columns", n))
        return
    }
    ml := G.Rows()
    if h == nil || !h.SizeMatch(ml, 1) {
        err = errors.New(fmt.Sprintf("'h' must matrix with size (%d,1)", ml))
        return
    }

    if A == nil {
        A = matrix.FloatZeros(0, n)
    }
    if b == nil {
        b = matrix.FloatZeros(0, 1)
    }
    if A.Cols() != n {
        err = errors.New(fmt.Sprintf("'A' must matrix with size %d columns", n))
        return
    }
    p := A.Rows()
    if b == nil || !b.SizeMatch(p, 1) {
        err = errors.New(fmt.Sprintf("'b' must matrix with size (%d,1)", p))
        return
    }

    dims := sets.NewDimensionSet("l", "q", "s")
    dims.Set("l", []int{ml})
    gpProg := createGpProg(K, F, g)

    return Cp(gpProg, G, h, A, b, dims, solopts)
}

// Returns min {t | x + t*e >= 0}, where e is defined as follows
//
//  - For the nonlinear and 'l' blocks: e is the vector of ones.
//  - For the 'q' blocks: e is the first unit vector.
//  - For the 's' blocks: e is the identity matrix.
//
// When called with the argument sigma, also returns the eigenvalues
// (in sigma) and the eigenvectors (in x) of the 's' components of x.
func maxStep(x *matrix.FloatMatrix, dims *sets.DimensionSet, mnl int, sigma *matrix.FloatMatrix) (rval float64, err error) {
	/*DEBUGGED*/

	rval = 0.0
	err = nil
	t := make([]float64, 0, 10)
	ind := mnl + dims.Sum("l")
	if ind > 0 {
		t = append(t, -minvec(x.FloatArray()[:ind]))
	}
	for _, m := range dims.At("q") {
		if m > 0 {
			v := blas.Nrm2Float(x, &la_.IOpt{"offset", ind + 1}, &la_.IOpt{"n", m - 1})
			v -= x.GetIndex(ind)
			t = append(t, v)
		}
		ind += m
	}

	//var Q *matrix.FloatMatrix
	//var w *matrix.FloatMatrix
	ind2 := 0
	//if sigma == nil && len(dims.At("s")) > 0 {
	//	mx := dims.Max("s")
	//	Q = matrix.FloatZeros(mx, mx)
	//	w = matrix.FloatZeros(mx, 1)
	//}
	for _, m := range dims.At("s") {
		if sigma == nil {
			Q := matrix.FloatZeros(m, m)
			w := matrix.FloatZeros(m, 1)
			blas.Copy(x, Q, &la_.IOpt{"offsetx", ind}, &la_.IOpt{"n", m * m})
			err = lapack.SyevrFloat(Q, w, nil, 0.0, nil, []int{1, 1}, la_.OptRangeInt,
				&la_.IOpt{"n", m}, &la_.IOpt{"lda", m})
			if m > 0 && err == nil {
				t = append(t, -w.GetIndex(0))
			}
		} else {
			err = lapack.SyevdFloat(x, sigma, la_.OptJobZValue, &la_.IOpt{"n", m},
				&la_.IOpt{"lda", m}, &la_.IOpt{"offseta", ind}, &la_.IOpt{"offsetw", ind2})
			if m > 0 {
				t = append(t, -sigma.GetIndex(ind2))
			}
		}
		ind += m * m
		ind2 += m
	}

	if len(t) > 0 {
		rval = maxvec(t)
	}
	return
}

// single invocation for matops and lapack functions
func runCheck(A *matrix.FloatMatrix, LB int) (bool, time.Duration, time.Duration) {

	var W *matrix.FloatMatrix = nil
	N := A.Cols()
	tau := matrix.FloatZeros(N, 1)
	if LB > 0 {
		W = matrix.FloatZeros(A.Rows(), LB)
	}
	fnc := func() {
		_, ERRmatops = matops.DecomposeQR(A, tau, W, LB)
	}

	if verbose && N < 10 {
		fmt.Fprintf(os.Stderr, "A start:\n%v\n", A)
	}
	A0 := A.Copy()
	tau0 := tau.Copy()

	mperf.FlushCache()
	time0 := mperf.Timeit(fnc)
	if verbose && N < 10 {
		fmt.Fprintf(os.Stderr, "A end:\n%v\n", A)
		tau.SetSize(1, N, 1)
		fmt.Fprintf(os.Stderr, "tau: %v\n", tau)
	}

	fn2 := func() {
		ERRlapack = lapack.Geqrf(A0, tau0)
	}

	mperf.FlushCache()
	time2 := mperf.Timeit(fn2)
	if verbose && N < 10 {
		fmt.Fprintf(os.Stderr, "A0 end:\n%v\n", A0)
		tau0.SetSize(1, N, 1) // row vector
		fmt.Fprintf(os.Stderr, "tau0: %v\n", tau0)
	}
	// now A == A0 && tau == tau0

	ok := A.AllClose(A0)
	oktau := tau.AllClose(tau0)
	if !ok || !oktau {
		// save result to globals
		Rlapack = A0
		Rmatops = A
		TAUlapack = tau0
		TAUmatops = tau
	}
	return ok && oktau, time0, time2
}

func _TestBK2(t *testing.T) {
	Bdata := [][]float64{
		[]float64{10.0, 20.0},
		[]float64{10.0, 20.0},
		[]float64{10.0, 20.0},
		[]float64{10.0, 20.0},
		[]float64{10.0, 20.0},
		[]float64{10.0, 20.0},
		[]float64{10.0, 20.0}}

	N := 7

	A0 := matrix.FloatNormal(N, N)
	A := matrix.FloatZeros(N, N)
	// A is symmetric, posivite definite
	Mult(A, A0, A0, 1.0, 1.0, TRANSB)

	X := matrix.FloatMatrixFromTable(Bdata, matrix.RowOrder)
	B := matrix.FloatZeros(N, 2)
	MultSym(B, A, X, 1.0, 0.0, LOWER|LEFT)
	t.Logf("initial B:\n%v\n", B)

	nb := 0
	W := matrix.FloatWithValue(A.Rows(), 5, 1.0)
	A.SetAt(4, 1, A.GetAt(4, 1)+1.0)
	A.SetAt(1, 4, A.GetAt(4, 1))

	ipiv := make([]int, N, N)
	L, _ := DecomposeBK(A.Copy(), W, ipiv, LOWER, nb)
	t.Logf("ipiv: %v\n", ipiv)
	t.Logf("L:\n%v\n", L)

	ipiv0 := make([]int, N, N)
	nb = 4
	L0, _ := DecomposeBK(A.Copy(), W, ipiv0, LOWER, nb)
	t.Logf("ipiv: %v\n", ipiv0)
	t.Logf("L:\n%v\n", L0)
	B0 := B.Copy()
	SolveBK(B0, L0, ipiv0, LOWER)
	t.Logf("B0:\n%v\n", B0)

	ipiv2 := make([]int32, N, N)
	lapack.Sytrf(A, ipiv2, linalg.OptLower)
	t.Logf("ipiv2: %v\n", ipiv2)
	t.Logf("lapack A:\n%v\n", A)
	lapack.Sytrs(A, B, ipiv2, linalg.OptLower)
	t.Logf("lapack B:\n%v\n", B)
	t.Logf("B == B0: %v\n", B.AllClose(B0))
}

func runRefTest(A *matrix.FloatMatrix, ntest, LB int) time.Duration {

	var mintime time.Duration

	N := A.Cols()
	tau := matrix.FloatZeros(N, 1)

	fnc := func() {
		ERRlapack = lapack.Geqrf(A, tau)
	}

	A0 := A.Copy()
	for n := 0; n < ntest; n++ {
		if n > 0 {
			// restore original A
			A0.CopyTo(A)
			tau.Scale(0.0)
		}
		mperf.FlushCache()
		time0 := mperf.Timeit(fnc)
		if n == 0 || time0 < mintime {
			mintime = time0
		}
	}
	return mintime
}

func _TestLDLnoPiv(t *testing.T) {
	N := 42
	nb := 8

	A0 := matrix.FloatUniform(N, N)
	A := matrix.FloatZeros(N, N)
	Mult(A, A0, A0, 1.0, 1.0, TRANSB)

	B := matrix.FloatNormal(A.Rows(), 2)
	w := matrix.FloatWithValue(A.Rows(), 2, 1.0)

	// B0 = A*B
	B0 := B.Copy()

	nb = 2
	L, _ := DecomposeLDLnoPiv(A.Copy(), w, LOWER, nb)
	Mult(B0, A, B, 1.0, 0.0, NOTRANS)
	SolveLDLnoPiv(B0, L, LOWER)
	t.Logf("L*D*L.T: ||B - A*X||_1: %e\n", NormP(B0.Minus(B), NORM_ONE))

	U, _ := DecomposeLDLnoPiv(A.Copy(), w, UPPER, nb)
	Mult(B0, A, B, 1.0, 0.0, NOTRANS)
	SolveLDLnoPiv(B0, U, UPPER)
	t.Logf("U*D*U.T: ||B - A*X||_1: %e\n", NormP(B0.Minus(B), NORM_ONE))

}