Golang FloatMatrix.Add方法代码示例

func _TestPartition1Dh(t *testing.T) {
	var AL, AR, A0, a1, A2 matrix.FloatMatrix
	A := matrix.FloatZeros(1, 6)
	partition1x2(&AL, &AR, A, 0, pRIGHT)
	t.Logf("m(AL)=%d, m(AR)=%d\n", AL.Cols(), AR.Cols())
	for AL.Cols() < A.Cols() {
		AR.Add(1.0)
		t.Logf("m(AR)=%d; %v\n", AR.Cols(), AR)
		repartition1x2to1x3(&AL, &A0, &a1, &A2, A, 1, pRIGHT)
		t.Logf("m(A0)=%d, m(A2)=%d, a1=%.1f\n", A0.Cols(), A2.Cols(), a1.Float())
		continue1x3to1x2(&AL, &AR, &A0, &a1, A, pRIGHT)
	}
}

func (F *cpProg) F1(xa MatrixVariable) (f MatrixVariable, Df MatrixVarDf, err error) {
	f = nil
	Df = nil
	err = nil
	x, x_ok := xa.(*epigraph)
	if !x_ok {
		err = errors.New("'x' argument not an epigraph")
	}
	var f0, Df0 *matrix.FloatMatrix
	f0, Df0, err = F.convexF.F1(x.m())
	if err != nil {
		fmt.Printf("'cpProg.F1' error: %v\n%s\n", err, x)
		return
	}
	f0.Add(-x.t(), 0)
	f = &matrixVar{f0}
	Df = &epigraphDf{Df0}
	return
}

func (F *cpProg) F2(xa MatrixVariable, z *matrix.FloatMatrix) (f MatrixVariable, Df MatrixVarDf, H MatrixVarH, err error) {
	f = nil
	Df = nil
	H = nil
	err = nil
	x, x_ok := xa.(*epigraph)
	if !x_ok {
		err = errors.New("'x' argument not an epigraph")
		return
	}
	var f0, Df0, H0 *matrix.FloatMatrix
	f0, Df0, H0, err = F.convexF.F2(x.m(), z)
	if err != nil {
		return
	}
	f0.Add(-x.t(), 0)
	f = &matrixVar{f0}
	Df = &epigraphDf{Df0}
	H = &epigraphH{H0}
	return
}

func _TestPartition2D(t *testing.T) {
	var ATL, ATR, ABL, ABR, As matrix.FloatMatrix
	var A00, a01, A02, a10, a11, a12, A20, a21, A22 matrix.FloatMatrix

	A := matrix.FloatZeros(6, 6)
	As.SubMatrixOf(A, 1, 1, 4, 4)
	As.SetIndexes(1.0)
	partition2x2(&ATL, &ATR, &ABL, &ABR, &As, 0)
	t.Logf("ATL:\n%v\n", &ATL)

	for ATL.Rows() < As.Rows() {
		repartition2x2to3x3(&ATL,
			&A00, &a01, &A02,
			&a10, &a11, &a12,
			&A20, &a21, &A22, &As, 1)
		t.Logf("m(a12)=%d [%d], m(a11)=%d\n", a12.Cols(), a12.NumElements(), a11.NumElements())
		a11.Add(1.0)
		a21.Add(-2.0)

		continue3x3to2x2(&ATL, &ATR, &ABL, &ABR, &A00, &a11, &A22, &As)
	}
	t.Logf("A:\n%v\n", A)
}

// unblocked LU decomposition with pivots: FLAME LU variant 3
func unblockedLUpiv(A *matrix.FloatMatrix, p *pPivots) error {
	var err error
	var ATL, ATR, ABL, ABR matrix.FloatMatrix
	var A00, a01, A02, a10, a11, a12, A20, a21, A22 matrix.FloatMatrix
	var AL, AR, A0, a1, A2, aB1, AB0 matrix.FloatMatrix
	var pT, pB, p0, p1, p2 pPivots

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

	for ATL.Rows() < A.Rows() && ATL.Cols() < A.Cols() {
		repartition2x2to3x3(&ATL,
			&A00, &a01, &A02,
			&a10, &a11, &a12,
			&A20, &a21, &A22 /**/, A, 1, pBOTTOMRIGHT)
		repartition1x2to1x3(&AL,
			&A0, &a1, &A2 /**/, A, 1, pRIGHT)
		repartPivot2x1to3x1(&pT,
			&p0, &p1, &p2 /**/, p, 1, pBOTTOM)

		// apply previously computed pivots
		applyPivots(&a1, &p0)

		// a01 = trilu(A00) \ a01 (TRSV)
		MVSolveTrm(&a01, &A00, 1.0, LOWER|UNIT)
		// a11 = a11 - a10 *a01
		a11.Add(Dot(&a10, &a01, -1.0))
		// a21 = a21 -A20*a01
		MVMult(&a21, &A20, &a01, -1.0, 1.0, NOTRANS)

		// pivot index on current column [a11, a21].T
		aB1.SubMatrixOf(&ABR, 0, 0, ABR.Rows(), 1)
		pivotIndex(&aB1, &p1)

		// pivots to current column
		applyPivots(&aB1, &p1)

		// a21 = a21 / a11
		InvScale(&a21, a11.Float())

		// apply pivots to previous columns
		AB0.SubMatrixOf(&ABL, 0, 0)
		applyPivots(&AB0, &p1)
		// scale last pivots to origin matrix row numbers
		p1.pivots[0] += ATL.Rows()

		continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &a11, &A22, A, pBOTTOMRIGHT)
		continue1x3to1x2(
			&AL, &AR, &A0, &a1, A, pRIGHT)
		contPivot3x1to2x1(
			&pT,
			&pB, &p0, &p1, p, pBOTTOM)
	}
	if ATL.Cols() < A.Cols() {
		applyPivots(&ATR, p)
		SolveTrm(&ATR, &ATL, 1.0, LEFT|UNIT|LOWER)
	}
	return err
}

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

	if G == nil {
		err = errors.New("'G' must be non-nil MatrixG interface.")
		return
	}
	fG := func(x, y MatrixVariable, alpha, beta float64, trans la.Option) error {
		return G.Gf(x, y, alpha, beta, trans)
	}

	// Check A and set defaults if it is nil
	fA := func(x, y MatrixVariable, alpha, beta float64, trans la.Option) error {
		return A.Af(x, y, alpha, beta, trans)
	}

	// Check b and set defaults if it is nil
	if b == nil {
		err = errors.New("'b' must be non-nil MatrixVariable interface.")
		return
	}

	// kktsolver(W) returns a routine for solving 3x3 block KKT system
	//
	//     [ 0   A'  G'*W^{-1} ] [ ux ]   [ bx ]
	//     [ A   0   0         ] [ uy ] = [ by ].
	//     [ G   0   -W'       ] [ uz ]   [ bz ]

	if kktsolver == nil {
		err = errors.New("nil kktsolver not allowed.")
		return
	}

	ws3 := matrix.FloatZeros(cdim, 1)
	wz3 := matrix.FloatZeros(cdim, 1)
	checkpnt.AddMatrixVar("ws3", ws3)
	checkpnt.AddMatrixVar("wz3", wz3)

	//
	res := func(ux, uy MatrixVariable, uz, us *matrix.FloatMatrix, vx, vy MatrixVariable, vz, vs *matrix.FloatMatrix, W *sets.FloatMatrixSet, lmbda *matrix.FloatMatrix) (err error) {
		// Evaluates residual in Newton equations:
		//
		//      [ vx ]    [ vx ]   [ 0     ]   [ P  A'  G' ]   [ ux        ]
		//      [ vy ] := [ vy ] - [ 0     ] - [ A  0   0  ] * [ uy        ]
		//      [ vz ]    [ vz ]   [ W'*us ]   [ G  0   0  ]   [ W^{-1}*uz ]
		//
		//      vs := vs - lmbda o (uz + us).

		// vx := vx - P*ux - A'*uy - G'*W^{-1}*uz
		minor := checkpnt.MinorTop()
		checkpnt.Check("00res", minor)
		fP(ux, vx, -1.0, 1.0)
		fA(uy, vx, -1.0, 1.0, la.OptTrans)
		blas.Copy(uz, wz3)
		scale(wz3, W, true, false)
		fG(&matrixVar{wz3}, vx, -1.0, 1.0, la.OptTrans)
		// vy := vy - A*ux
		fA(ux, vy, -1.0, 1.0, la.OptNoTrans)
		checkpnt.Check("50res", minor)

		// vz := vz - G*ux - W'*us
		fG(ux, &matrixVar{vz}, -1.0, 1.0, la.OptNoTrans)
		blas.Copy(us, ws3)
		scale(ws3, W, true, false)
		blas.AxpyFloat(ws3, vz, -1.0)

		// vs := vs - lmbda o (uz + us)
		blas.Copy(us, ws3)