Golang FloatMatrix.Scale方法代码示例

func _TestViewUpdate(t *testing.T) {
	Adata2 := [][]float64{
		[]float64{4.0, 2.0, 2.0},
		[]float64{6.0, 4.0, 2.0},
		[]float64{4.0, 6.0, 1.0},
	}

	A := matrix.FloatMatrixFromTable(Adata2, matrix.RowOrder)
	N := A.Rows()

	// simple LU decomposition without pivoting
	var A11, a10, a01, a00 matrix.FloatMatrix
	for k := 1; k < N; k++ {
		a00.SubMatrixOf(A, k-1, k-1, 1, 1)
		a01.SubMatrixOf(A, k-1, k, 1, A.Cols()-k)
		a10.SubMatrixOf(A, k, k-1, A.Rows()-k, 1)
		A11.SubMatrixOf(A, k, k)
		//t.Logf("A11: %v  a01: %v\n", A11, a01)
		a10.Scale(1.0 / a00.Float())
		MVRankUpdate(&A11, &a10, &a01, -1.0)
	}

	Ld := TriLU(A.Copy())
	Ud := TriU(A)
	t.Logf("Ld:\n%v\nUd:\n%v\n", Ld, Ud)
	An := matrix.FloatZeros(N, N)
	Mult(An, Ld, Ud, 1.0, 1.0, NOTRANS)
	t.Logf("A == Ld*Ud: %v\n", An.AllClose(An))
}

func (g *matrixFs) Gf(x, y *matrix.FloatMatrix, alpha, beta float64, trans linalg.Option) error {
	//
	minor := 0
	if !checkpnt.MinorEmpty() {
		minor = checkpnt.MinorTop()
	} else {
		loopg += 1
		minor = loopg
	}
	checkpnt.Check("00-Gfunc", minor)

	m, n := g.A.Size()
	y.Scale(beta)

	// x_n = x[:n]
	//x_n := matrix.FloatVector(x.FloatArray()[:n])
	x_n := x.SubMatrix(0, 0, n, 1).Copy()

	// x_n_2n = x[n:2*n]
	//x_n_2n := matrix.FloatVector(x.FloatArray()[n : 2*n])
	x_n_2n := x.SubMatrix(n, 0, n, 1).Copy()

	if linalg.Equal(trans, linalg.OptNoTrans) {
		// y += alpha * G * x

		// y[:n] += alpha * (x[:n] - x[n:2*n])
		y_n := matrix.Minus(x_n, x_n_2n).Scale(alpha)
		y.SubMatrix(0, 0, n, 1).Plus(y_n)
		//y.AddIndexes(matrix.Indexes(n), y_n.FloatArray())

		// y[n:2*n] += alpha * (-x[:n] - x[n:2*n]) = -alpha * (x[:n]+x[n:2*n])
		y_n = matrix.Plus(x_n, x_n_2n).Scale(-alpha)
		y.SubMatrix(n, 0, n, 1).Plus(y_n)
		//y.AddIndexes(matrix.Indexes(n, 2*n), y_n.FloatArray())

		// y[2*n+1:] += -alpha * A * x[:n]
		y_2n := matrix.Times(g.A, x_n).Scale(-alpha)
		//y.AddIndexes(matrix.Indexes(2*n+1, y.NumElements()), y_2n.FloatArray())
		y.SubMatrix(2*n+1, 0, y_2n.NumElements(), 1).Plus(y_2n)
	} else {
		// x_m = x[-m:]
		//x_m := matrix.FloatVector(x.FloatArray()[x.NumElements()-m:])
		x_m := x.SubMatrix(x.NumElements()-m, 0)

		// x_tmp = (x[:n] - x[n:2*n] - A.T * x[-m:])
		x_tmp := matrix.Minus(x_n, x_n_2n, matrix.Times(g.A.Transpose(), x_m))

		// y[:n] += alpha * (x[:n] - x[n:2*n] - A.T * x[-m:])
		//y.AddIndexes(matrix.Indexes(n), x_tmp.Scale(alpha).FloatArray())
		y.SubMatrix(0, 0, n, 1).Plus(x_tmp.Scale(alpha))

		x_tmp = matrix.Plus(x_n, x_n_2n).Scale(-alpha)
		//y.AddIndexes(matrix.Indexes(n, y.NumElements()), x_tmp.FloatArray())
		y.SubMatrix(n, 0).Plus(x_tmp)
	}
	checkpnt.Check("10-Gfunc", minor)
	return nil
}

/*
 * Compute
 *   C = C*diag(D)      flags & RIGHT == true
 *   C = diag(D)*C      flags & LEFT  == true
 *
 * Arguments
 *   C     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 MultDiag(C, D *matrix.FloatMatrix, flags Flags) {
	var c, d0 matrix.FloatMatrix
	if D.Cols() == 1 {
		// diagonal is column vector
		switch flags & (LEFT | RIGHT) {
		case LEFT:
			// scale rows; for each column element-wise multiply with D-vector
			for k := 0; k < C.Cols(); k++ {
				C.SubMatrix(&c, 0, k, C.Rows(), 1)
				c.Mul(D)
			}
		case RIGHT:
			// scale columns
			for k := 0; k < C.Cols(); k++ {
				C.SubMatrix(&c, 0, k, C.Rows(), 1)
				// scale the column
				c.Scale(D.GetAt(k, 0))
			}
		}
	} else {
		// diagonal is row vector
		var d *matrix.FloatMatrix
		if D.Rows() == 1 {
			d = D
		} else {
			D.SubMatrix(&d0, 0, 0, 1, D.Cols(), D.LeadingIndex()+1)
			d = &d0
		}
		switch flags & (LEFT | RIGHT) {
		case LEFT:
			for k := 0; k < C.Rows(); k++ {
				C.SubMatrix(&c, k, 0, 1, C.Cols())
				// scale the row
				c.Scale(d.GetAt(0, k))
			}
		case RIGHT:
			// scale columns
			for k := 0; k < C.Cols(); k++ {
				C.SubMatrix(&c, 0, k, C.Rows(), 1)
				// scale the column
				c.Scale(d.GetAt(0, k))
			}
		}
	}
}

func blkDecompBKLower(A, W *matrix.FloatMatrix, p *pPivots, nb int) (err error) {
	var ATL, ATR, ABL, ABR matrix.FloatMatrix
	var A00, A10, A11, A20, A21, A22 matrix.FloatMatrix
	var wrk matrix.FloatMatrix
	var pT, pB, p0, p1, p2 pPivots
	var nblk int = 0

	err = nil
	partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, 0, 0, pTOPLEFT)
	partitionPivot2x1(
		&pT,
		&pB, p, 0, pTOP)

	for ABR.Cols() >= nb {
		err, nblk = unblkBoundedBKLower(&ABR, W, &pB, nb)

		// repartition nblk size
		repartition2x2to3x3(&ATL,
			&A00, nil, nil,
			&A10, &A11, nil,
			&A20, &A21, &A22, A, nblk, pBOTTOMRIGHT)
		repartPivot2x1to3x1(&pT,
			&p0, &p1, &p2 /**/, p, nblk, pBOTTOM)

		// --------------------------------------------------------
		// here [A11;A21] has been decomposed by unblkBoundedBKLower()
		// Now we need update A22

		// wrk is original A21
		W.SubMatrix(&wrk, nblk, 0, A21.Rows(), nblk)

		// A22 = A22 - L21*D1*L21.T = A22 - L21*W.T
		UpdateTrm(&A22, &A21, &wrk, -1.0, 1.0, LOWER|TRANSB)

		// partially undo row pivots left of diagonal
		for k := nblk; k > 0; k-- {
			var s, d matrix.FloatMatrix
			r := p1.pivots[k-1]
			rlen := k - 1
			if r < 0 {
				r = -r
				rlen--
			}
			if r == k {
				// no pivot
				continue
			}
			ABR.SubMatrix(&s, k-1, 0, 1, rlen)
			ABR.SubMatrix(&d, r-1, 0, 1, rlen)
			Swap(&d, &s)

			if p1.pivots[k-1] < 0 {
				k-- // skip other entry in 2x2 pivots
			}
		}

		// shift pivot values
		for k, n := range p1.pivots {
			if n > 0 {
				p1.pivots[k] += ATL.Rows()
			} else {
				p1.pivots[k] -= ATL.Rows()
			}
		}

		// zero work for debuging
		W.Scale(0.0)

		// ---------------------------------------------------------

		continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &A11, &A22, A, pBOTTOMRIGHT)
		contPivot3x1to2x1(
			&pT,
			&pB, &p0, &p1, p, pBOTTOM)
	}

	// do the last part with unblocked code
	if ABR.Cols() > 0 {
		unblkDecompBKLower(&ABR, W, &pB)
		// shift pivot values
		for k, n := range pB.pivots {
			if n > 0 {
				pB.pivots[k] += ATL.Rows()
			} else {
				pB.pivots[k] -= ATL.Rows()
			}
		}
	}
	return
}

func main() {
	m := 6
	Vdata := [][]float64{
		[]float64{1.0, -1.0, -2.0, -2.0, 0.0, 1.5, 1.0},
		[]float64{1.0, 2.0, 1.0, -1.0, -2.0, -1.0, 1.0}}

	V := matrix.FloatMatrixFromTable(Vdata, matrix.RowOrder)

	// V[1, :m] - V[1,1:]
	a0 := matrix.Minus(V.GetSubMatrix(1, 0, 1, m), V.GetSubMatrix(1, 1, 1))
	// V[0, :m] - V[0,1:]
	a1 := matrix.Minus(V.GetSubMatrix(0, 0, 1, m), V.GetSubMatrix(0, 1, 1))
	A0, _ := matrix.FloatMatrixStacked(matrix.StackDown, a0.Scale(-1.0), a1)
	A0 = A0.Transpose()
	b0 := matrix.Mul(A0, V.GetSubMatrix(0, 0, 2, m).Transpose())
	b0 = matrix.Times(b0, matrix.FloatWithValue(2, 1, 1.0))

	A := make([]*matrix.FloatMatrix, 0)
	b := make([]*matrix.FloatMatrix, 0)
	A = append(A, A0)
	b = append(b, b0)

	// List of symbols
	C := make([]*matrix.FloatMatrix, 0)
	C = append(C, matrix.FloatZeros(2, 1))
	var row *matrix.FloatMatrix = nil
	for k := 0; k < m; k++ {
		row = A0.GetRow(k, row)
		nrm := blas.Nrm2Float(row)
		row.Scale(2.0 * b0.GetIndex(k) / (nrm * nrm))
		C = append(C, row.Transpose())
	}

	// Voronoi set around C[1]
	A1 := matrix.FloatZeros(3, 2)
	A1.SetSubMatrix(0, 0, A0.GetSubMatrix(0, 0, 1).Scale(-1.0))
	A1.SetSubMatrix(1, 0, matrix.Minus(C[m], C[1]).Transpose())
	A1.SetSubMatrix(2, 0, matrix.Minus(C[2], C[1]).Transpose())
	b1 := matrix.FloatZeros(3, 1)
	b1.SetIndex(0, -b0.GetIndex(0))
	v := matrix.Times(A1.GetRow(1, nil), matrix.Plus(C[m], C[1])).Float() * 0.5
	b1.SetIndex(1, v)
	v = matrix.Times(A1.GetRow(2, nil), matrix.Plus(C[2], C[1])).Float() * 0.5
	b1.SetIndex(2, v)
	A = append(A, A1)
	b = append(b, b1)

	// Voronoi set around C[2] ... C[5]
	for k := 2; k < 6; k++ {
		A1 = matrix.FloatZeros(3, 2)
		A1.SetSubMatrix(0, 0, A0.GetSubMatrix(k-1, 0, 1).Scale(-1.0))
		A1.SetSubMatrix(1, 0, matrix.Minus(C[k-1], C[k]).Transpose())
		A1.SetSubMatrix(2, 0, matrix.Minus(C[k+1], C[k]).Transpose())
		b1 = matrix.FloatZeros(3, 1)
		b1.SetIndex(0, -b0.GetIndex(k-1))
		v := matrix.Times(A1.GetRow(1, nil), matrix.Plus(C[k-1], C[k])).Float() * 0.5
		b1.SetIndex(1, v)
		v = matrix.Times(A1.GetRow(2, nil), matrix.Plus(C[k+1], C[k])).Float() * 0.5
		b1.SetIndex(2, v)
		A = append(A, A1)
		b = append(b, b1)
	}

	// Voronoi set around C[6]
	A1 = matrix.FloatZeros(3, 2)
	A1.SetSubMatrix(0, 0, A0.GetSubMatrix(5, 0, 1).Scale(-1.0))
	A1.SetSubMatrix(1, 0, matrix.Minus(C[1], C[6]).Transpose())
	A1.SetSubMatrix(2, 0, matrix.Minus(C[5], C[6]).Transpose())
	b1 = matrix.FloatZeros(3, 1)
	b1.SetIndex(0, -b0.GetIndex(5))
	v = matrix.Times(A1.GetRow(1, nil), matrix.Plus(C[1], C[6])).Float() * 0.5
	b1.SetIndex(1, v)
	v = matrix.Times(A1.GetRow(2, nil), matrix.Plus(C[5], C[6])).Float() * 0.5
	b1.SetIndex(2, v)

	A = append(A, A1)
	b = append(b, b1)

	P := matrix.FloatIdentity(2)
	q := matrix.FloatZeros(2, 1)
	solopts := &cvx.SolverOptions{ShowProgress: false, MaxIter: 30}
	ovals := make([]float64, 0)
	for k := 1; k < 7; k++ {
		sol, err := cvx.Qp(P, q, A[k], b[k], nil, nil, solopts, nil)
		_ = err
		x := sol.Result.At("x")[0]
		ovals = append(ovals, math.Pow(blas.Nrm2Float(x), 2.0))
	}

	optvals := matrix.FloatVector(ovals)
	//fmt.Printf("optvals=\n%v\n", optvals)

	rangeFunc := func(n int) []float64 {
		r := make([]float64, 0)
		for i := 0; i < n; i++ {
			r = append(r, float64(i))
		}
		return r
	}

//.........这里部分代码省略.........

// Build Q in place by applying elementary reflectors in reverse order to
// an implied identity matrix.  This forms Q = H(1)H(2) ... H(k)
//
// this is compatibe with lapack.DORG2R
func unblockedBuildQ(A, tau, w *matrix.FloatMatrix, kb int) error {
	var err error = nil
	var ATL, ATR, ABL, ABR matrix.FloatMatrix
	var A00, a01, A02, a10t, a11, a12t, A20, a21, A22 matrix.FloatMatrix
	var tT, tB matrix.FloatMatrix
	var t0, tau1, t2, w1 matrix.FloatMatrix
	var mb int
	var rowvec bool

	mb = A.Rows() - A.Cols()
	rowvec = tau.Rows() == 1

	partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, mb, 0, pBOTTOMRIGHT)

	if rowvec {
		partition1x2(
			&tT, &tB, tau, 0, pRIGHT)
	} else {
		partition2x1(
			&tT,
			&tB, tau, 0, pBOTTOM)
	}

	// clearing of the columns of the right and setting ABR to unit diagonal
	// (only if not applying all reflectors, kb > 0)

	for ATL.Rows() > 0 && ATL.Cols() > 0 {
		repartition2x2to3x3(&ATL,
			&A00, &a01, &A02,
			&a10t, &a11, &a12t,
			&A20, &a21, &A22, A, 1, pTOPLEFT)
		if rowvec {
			repartition1x2to1x3(&tT,
				&t0, &tau1, &t2, tau, 1, pLEFT)
		} else {
			repartition2x1to3x1(&tT,
				&t0,
				&tau1,
				&t2, tau, 1, pTOP)
		}

		// --------------------------------------------------------

		// adjust workspace to correct size
		w.SubMatrix(&w1, 0, 0, 1, a12t.Cols())
		// apply Householder reflection from left
		applyHHTo2x1(&tau1, &a21, &a12t, &A22, &w1, LEFT)

		// apply (in-place) current elementary reflector to unit vector
		a21.Scale(-tau1.Float())
		a11.SetAt(0, 0, 1.0-tau1.Float())

		// zero the upper part
		a01.SetIndexes(0.0)

		// --------------------------------------------------------
		continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &a11, &A22, A, pTOPLEFT)
		if rowvec {
			continue1x3to1x2(
				&tT, &tB, &t0, &tau1, tau, pLEFT)
		} else {
			continue3x1to2x1(
				&tT,
				&tB, &t0, &tau1, tau, pTOP)
		}
	}
	return err
}