Golang DimensionSet.Sum方法代碼示例

/*
   Matrix-vector multiplication.

   A is a matrix or spmatrix of size (m, n) where

       N = dims['l'] + sum(dims['q']) + sum( k**2 for k in dims['s'] )

   representing a mapping from R^n to S.

   If trans is 'N':

       y := alpha*A*x + beta * y   (trans = 'N').

   x is a vector of length n.  y is a vector of length N.

   If trans is 'T':

       y := alpha*A'*x + beta * y  (trans = 'T').

   x is a vector of length N.  y is a vector of length n.

   The 's' components in S are stored in unpacked 'L' storage.
*/
func sgemv(A, x, y *matrix.FloatMatrix, alpha, beta float64, dims *sets.DimensionSet, opts ...la_.Option) error {

	m := dims.Sum("l", "q") + dims.SumSquared("s")
	n := la_.GetIntOpt("n", -1, opts...)
	if n == -1 {
		n = A.Cols()
	}
	trans := la_.GetIntOpt("trans", int(la_.PNoTrans), opts...)
	offsetX := la_.GetIntOpt("offsetx", 0, opts...)
	offsetY := la_.GetIntOpt("offsety", 0, opts...)
	offsetA := la_.GetIntOpt("offseta", 0, opts...)

	if trans == int(la_.PTrans) && alpha != 0.0 {
		trisc(x, dims, offsetX)
		//fmt.Printf("trisc x=\n%v\n", x.ConvertToString())
	}
	//fmt.Printf("alpha=%.4f beta=%.4f m=%d n=%d\n", alpha, beta, m, n)
	//fmt.Printf("A=\n%v\nx=\n%v\ny=\n%v\n", A, x.ConvertToString(), y.ConvertToString())
	err := blas.GemvFloat(A, x, y, alpha, beta, &la_.IOpt{"trans", trans},
		&la_.IOpt{"n", n}, &la_.IOpt{"m", m}, &la_.IOpt{"offseta", offsetA},
		&la_.IOpt{"offsetx", offsetX}, &la_.IOpt{"offsety", offsetY})
	//fmt.Printf("gemv y=\n%v\n", y.ConvertToString())

	if trans == int(la_.PTrans) && alpha != 0.0 {
		triusc(x, dims, offsetX)
	}
	return err
}

// In-place version of pack(), which also accepts matrix arguments x.
// The columns of x are elements of S, with the 's' components stored
// in unpacked storage.  On return, the 's' components are stored in
// packed storage and the off-diagonal entries are scaled by sqrt(2).
//
func pack2(x *matrix.FloatMatrix, dims *sets.DimensionSet, mnl int) (err error) {
	if len(dims.At("s")) == 0 {
		return nil
	}

	const sqrt2 = 1.41421356237309504880

	iu := mnl + dims.Sum("l", "q")
	ip := iu
	row := matrix.FloatZeros(1, x.Cols())
	//fmt.Printf("x.size = %d %d\n", x.Rows(), x.Cols())
	for _, n := range dims.At("s") {
		for k := 0; k < n; k++ {
			cnt := n - k
			row = x.GetRow(iu+(n+1)*k, row)
			//fmt.Printf("%02d: %v\n", iu+(n+1)*k, x.FloatArray())
			x.SetRow(ip, row)
			for i := 1; i < n-k; i++ {
				row = x.GetRow(iu+(n+1)*k+i, row)
				//fmt.Printf("%02d: %v\n", iu+(n+1)*k+i, x.FloatArray())
				x.SetRow(ip+i, row.Scale(sqrt2))
			}
			ip += cnt
		}
		iu += n * n
	}
	return nil
}

/*
   Copy x to y using packed storage.

   The vector x is an element of S, with the 's' components stored in
   unpacked storage.  On return, x is copied to y with the 's' components
   stored in packed storage and the off-diagonal entries scaled by
   sqrt(2).
*/
func pack(x, y *matrix.FloatMatrix, dims *sets.DimensionSet, opts ...la_.Option) (err error) {
	/*DEBUGGED*/
	err = nil
	mnl := la_.GetIntOpt("mnl", 0, opts...)
	offsetx := la_.GetIntOpt("offsetx", 0, opts...)
	offsety := la_.GetIntOpt("offsety", 0, opts...)

	nlq := mnl + dims.At("l")[0] + dims.Sum("q")
	blas.Copy(x, y, &la_.IOpt{"n", nlq}, &la_.IOpt{"offsetx", offsetx},
		&la_.IOpt{"offsety", offsety})

	iu, ip := offsetx+nlq, offsety+nlq
	for _, n := range dims.At("s") {
		for k := 0; k < n; k++ {
			blas.Copy(x, y, &la_.IOpt{"n", n - k}, &la_.IOpt{"offsetx", iu + k*(n+1)},
				&la_.IOpt{"offsety", ip})
			y.SetIndex(ip, (y.GetIndex(ip) / math.Sqrt(2.0)))
			ip += n - k
		}
		iu += n * n
	}
	np := dims.SumPacked("s")
	blas.ScalFloat(y, math.Sqrt(2.0), &la_.IOpt{"n", np}, &la_.IOpt{"offset", offsety + nlq})
	return
}

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

/*
   Scales the strictly lower triangular part of the 's' components of x
   by 0.5.

*/
func triusc(x *matrix.FloatMatrix, dims *sets.DimensionSet, offset int) error {

	//m := dims.Sum("l", "q") + dims.SumSquared("s")
	ind := offset + dims.Sum("l", "q")

	for _, mk := range dims.At("s") {
		for j := 1; j < mk; j++ {
			blas.ScalFloat(x, 0.5, &la_.IOpt{"n", mk - j}, &la_.IOpt{"offset", ind + mk*(j-1) + j})
		}
		ind += mk * mk
	}
	return nil
}

// Inner product of two vectors in S.
func sdot(x, y *matrix.FloatMatrix, dims *sets.DimensionSet, mnl int) float64 {
	/*DEBUGGED*/
	ind := mnl + dims.At("l")[0] + dims.Sum("q")
	a := blas.DotFloat(x, y, &la_.IOpt{"n", ind})
	for _, m := range dims.At("s") {
		a += blas.DotFloat(x, y, &la_.IOpt{"offsetx", ind}, &la_.IOpt{"offsety", ind},
			&la_.IOpt{"incx", m + 1}, &la_.IOpt{"incy", m + 1}, &la_.IOpt{"n", m})
		for j := 1; j < m; j++ {
			a += 2.0 * blas.DotFloat(x, y, &la_.IOpt{"offsetx", ind + j}, &la_.IOpt{"offsety", ind + j},
				&la_.IOpt{"incx", m + 1}, &la_.IOpt{"incy", m + 1}, &la_.IOpt{"n", m - j})
		}
		ind += m * m
	}
	return a
}

/*
   Sets upper triangular part of the 's' components of x equal to zero
   and scales the strictly lower triangular part by 2.0.
*/
func trisc(x *matrix.FloatMatrix, dims *sets.DimensionSet, offset int) error {

	//m := dims.Sum("l", "q") + dims.SumSquared("s")
	ind := offset + dims.Sum("l", "q")

	for _, mk := range dims.At("s") {
		for j := 1; j < mk; j++ {
			blas.ScalFloat(x, 0.0, la_.IntOpt("n", mk-j), la_.IntOpt("inc", mk),
				la_.IntOpt("offset", ind+j*(mk+1)-1))
			blas.ScalFloat(x, 2.0, la_.IntOpt("n", mk-j), la_.IntOpt("offset", ind+mk*(j-1)+j))
		}
		ind += mk * mk
	}
	return nil
}

// The product x := y o y.   The 's' components of y are diagonal and
// only the diagonals of x and y are stored.
func ssqr(x, y *matrix.FloatMatrix, dims *sets.DimensionSet, mnl int) (err error) {
	/*DEBUGGED*/
	blas.Copy(y, x)
	ind := mnl + dims.At("l")[0]
	err = blas.Tbmv(y, x, &la_.IOpt{"n", ind}, &la_.IOpt{"k", 0}, &la_.IOpt{"lda", 1})
	if err != nil {
		return
	}

	for _, m := range dims.At("q") {
		v := blas.Nrm2Float(y, &la_.IOpt{"n", m}, &la_.IOpt{"offset", ind})
		x.SetIndex(ind, v*v)
		blas.ScalFloat(x, 2.0*y.GetIndex(ind), &la_.IOpt{"n", m - 1}, &la_.IOpt{"offset", ind + 1})
		ind += m
	}
	err = blas.Tbmv(y, x, &la_.IOpt{"n", dims.Sum("s")}, &la_.IOpt{"k", 0},
		&la_.IOpt{"lda", 1}, &la_.IOpt{"offseta", ind}, &la_.IOpt{"offsetx", ind})
	return
}

/*
   The vector x is an element of S, with the 's' components stored
   in unpacked storage and off-diagonal entries scaled by sqrt(2).
   On return, x is copied to y with the 's' components stored in
   unpacked storage.

*/
func unpack(x, y *matrix.FloatMatrix, dims *sets.DimensionSet, opts ...la_.Option) (err error) {
	/*DEBUGGED*/
	err = nil
	mnl := la_.GetIntOpt("mnl", 0, opts...)
	offsetx := la_.GetIntOpt("offsetx", 0, opts...)
	offsety := la_.GetIntOpt("offsety", 0, opts...)

	nlq := mnl + dims.At("l")[0] + dims.Sum("q")
	err = blas.Copy(x, y, &la_.IOpt{"n", nlq}, &la_.IOpt{"offsetx", offsetx},
		&la_.IOpt{"offsety", offsety})
	if err != nil {
		return
	}

	ip, iu := offsetx+nlq, offsety+nlq
	for _, n := range dims.At("s") {
		for k := 0; k < n; k++ {
			err = blas.Copy(x, y, &la_.IOpt{"n", n - k}, &la_.IOpt{"offsetx", ip},
				&la_.IOpt{"offsety", iu + k*(n+1)})
			if err != nil {
				return
			}

			ip += n - k
			blas.ScalFloat(y, 1.0/math.Sqrt(2.0),
				&la_.IOpt{"n", n - k - 1}, &la_.IOpt{"offset", iu + k*(n+1) + 1})
		}
		iu += n * n
	}
	/*
		nu := dims.SumSquared("s")
		fmt.Printf("-- UnPack: nu=%d, offset=%d\n", nu, offsety+nlq)
		err = blas.ScalFloat(y,
			&la_.IOpt{"n", nu}, &la_.IOpt{"offset", offsety+nlq})
	*/
	return
}

// Solves a convex optimization problem with a linear objective
//
//        minimize    c'*x
//        subject to  f(x) <= 0
//                    G*x <= h
//                    A*x = b.
//
// using custom KTT equation solver and custom constraints G and A.
//
func CplCustomMatrix(F ConvexProg, c *matrix.FloatMatrix, G MatrixG, h *matrix.FloatMatrix,
	A MatrixA, b *matrix.FloatMatrix, dims *sets.DimensionSet, kktsolver KKTCpSolver,
	solopts *SolverOptions) (sol *Solution, err error) {

	var mnl int
	var x0 *matrix.FloatMatrix

	mnl, x0, err = F.F0()
	if err != nil {
		return
	}

	if x0.Cols() != 1 {
		err = errors.New("'x0' must be matrix with one column")
		return
	}
	if c == nil {
		err = errors.New("'c' must be non nil matrix")
		return
	}
	if !c.SizeMatch(x0.Size()) {
		err = errors.New(fmt.Sprintf("'c' must be matrix of size (%d,1)", x0.Rows()))
		return
	}

	if h == nil {
		h = matrix.FloatZeros(0, 1)
	}
	if h.Cols() > 1 {
		err = errors.New("'h' must be matrix with 1 column")
		return
	}

	if dims == nil {
		dims = sets.NewDimensionSet("l", "q", "s")
		dims.Set("l", []int{h.Rows()})
	}

	cdim := dims.Sum("l", "q") + dims.SumSquared("s")

	if h.Rows() != cdim {
		err = errors.New(fmt.Sprintf("'h' must be float matrix of size (%d,1)", cdim))
		return
	}

	// Check b and set defaults if it is nil
	if b == nil {
		b = matrix.FloatZeros(0, 1)
	}
	if b.Cols() != 1 {
		estr := fmt.Sprintf("'b' must be a matrix with 1 column")
		err = errors.New(estr)
		return
	}

	mc := matrixVar{c}
	mb := matrixVar{b}
	var mG MatrixVarG
	var mA MatrixVarA

	if G == nil {
		mG = &matrixVarG{matrix.FloatZeros(0, c.Rows()), dims}
	} else {
		mG = &matrixIfG{G}
	}
	if A == nil {
		mA = &matrixVarA{matrix.FloatZeros(0, c.Rows())}
	} else {
		mA = &matrixIfA{A}
	}

	return cpl_problem(F, &mc, mG, h, mA, &mb, dims, kktsolver, solopts, x0, mnl)
}

// Solves a convex optimization problem with a linear objective
//
//        minimize    c'*x
//        subject to  f(x) <= 0
//                    G*x <= h
//                    A*x = b.
//
// using custom KTT equation solver.
//
func CplCustomKKT(F ConvexProg, c *matrix.FloatMatrix, G, h, A, b *matrix.FloatMatrix,
	dims *sets.DimensionSet, kktsolver KKTCpSolver,
	solopts *SolverOptions) (sol *Solution, err error) {

	var mnl int
	var x0 *matrix.FloatMatrix

	mnl, x0, err = F.F0()
	if err != nil {
		return
	}

	if x0.Cols() != 1 {
		err = errors.New("'x0' must be matrix with one column")
		return
	}
	if c == nil {
		err = errors.New("'c' must be non nil matrix")
		return
	}
	if !c.SizeMatch(x0.Size()) {
		err = errors.New(fmt.Sprintf("'c' must be matrix of size (%d,1)", x0.Rows()))
		return
	}

	if h == nil {
		h = matrix.FloatZeros(0, 1)
	}
	if h.Cols() > 1 {
		err = errors.New("'h' must be matrix with 1 column")
		return
	}

	if dims == nil {
		dims = sets.NewDimensionSet("l", "q", "s")
		dims.Set("l", []int{h.Rows()})
	}

	cdim := dims.Sum("l", "q") + dims.SumSquared("s")

	if h.Rows() != cdim {
		err = errors.New(fmt.Sprintf("'h' must be float matrix of size (%d,1)", cdim))
		return
	}

	if G == nil {
		G = matrix.FloatZeros(0, c.Rows())
	}
	if !G.SizeMatch(cdim, c.Rows()) {
		estr := fmt.Sprintf("'G' must be of size (%d,%d)", cdim, c.Rows())
		err = errors.New(estr)
		return
	}

	// Check A and set defaults if it is nil
	if A == nil {
		// zeros rows reduces Gemv to vector products
		A = matrix.FloatZeros(0, c.Rows())
	}
	if A.Cols() != c.Rows() {
		estr := fmt.Sprintf("'A' must have %d columns", c.Rows())
		err = errors.New(estr)
		return
	}

	// Check b and set defaults if it is nil
	if b == nil {
		b = matrix.FloatZeros(0, 1)
	}
	if b.Cols() != 1 {
		estr := fmt.Sprintf("'b' must be a matrix with 1 column")
		err = errors.New(estr)
		return
	}
	if b.Rows() != A.Rows() {
		estr := fmt.Sprintf("'b' must have length %d", A.Rows())
		err = errors.New(estr)
		return
	}

	var mc = matrixVar{c}
	var mb = matrixVar{b}
	var mA = matrixVarA{A}
	var mG = matrixVarG{G, dims}

	return cpl_problem(F, &mc, &mG, h, &mA, &mb, dims, kktsolver, solopts, x0, mnl)
}

/*
   Evaluates

       x := H(lambda^{1/2}) * x   (inverse is 'N')
       x := H(lambda^{-1/2}) * x  (inverse is 'I').

   H is the Hessian of the logarithmic barrier.

*/
func scale2(lmbda, x *matrix.FloatMatrix, dims *sets.DimensionSet, mnl int, inverse bool) (err error) {
	err = nil

	//var minor int = 0
	//if ! checkpnt.MinorEmpty() {
	//	minor = checkpnt.MinorTop()
	//}

	//fmt.Printf("\n%d.%04d scale2 x=\n%v\nlmbda=\n%v\n", checkpnt.Major(), minor,
	//	x.ToString("%.17f"), lmbda.ToString("%.17f"))

	//if ! checkpnt.MinorEmpty() {
	//	checkpnt.Check("000scale2", minor)
	//}

	// For the nonlinear and 'l' blocks,
	//
	//     xk := xk ./ l   (inverse is 'N')
	//     xk := xk .* l   (inverse is 'I')
	//
	// where l is lmbda[:mnl+dims['l']].
	ind := mnl + dims.Sum("l")
	if !inverse {
		blas.TbsvFloat(lmbda, x, &la_.IOpt{"n", ind}, &la_.IOpt{"k", 0}, &la_.IOpt{"lda", 1})
	} else {
		blas.TbmvFloat(lmbda, x, &la_.IOpt{"n", ind}, &la_.IOpt{"k", 0}, &la_.IOpt{"lda", 1})
	}

	//if ! checkpnt.MinorEmpty() {
	//	checkpnt.Check("010scale2", minor)
	//}

	// For 'q' blocks, if inverse is 'N',
	//
	//     xk := 1/a * [ l'*J*xk;
	//         xk[1:] - (xk[0] + l'*J*xk) / (l[0] + 1) * l[1:] ].
	//
	// If inverse is 'I',
	//
	//     xk := a * [ l'*xk;
	//         xk[1:] + (xk[0] + l'*xk) / (l[0] + 1) * l[1:] ].
	//
	// a = sqrt(lambda_k' * J * lambda_k), l = lambda_k / a.
	for _, m := range dims.At("q") {
		var lx, a, c, x0 float64
		a = jnrm2(lmbda, m, ind) //&la_.IOpt{"n", m}, &la_.IOpt{"offset", ind})
		if !inverse {
			lx = jdot(lmbda, x, m, ind, ind) //&la_.IOpt{"n", m}, &la_.IOpt{"offsetx", ind},
			//&la_.IOpt{"offsety", ind})
			lx /= a
		} else {
			lx = blas.DotFloat(lmbda, x, &la_.IOpt{"n", m}, &la_.IOpt{"offsetx", ind},
				&la_.IOpt{"offsety", ind})
			lx /= a
		}
		x0 = x.GetIndex(ind)
		x.SetIndex(ind, lx)
		c = (lx + x0) / (lmbda.GetIndex(ind)/a + 1.0) / a
		if !inverse {
			c *= -1.0
		}
		blas.AxpyFloat(lmbda, x, c, &la_.IOpt{"n", m - 1}, &la_.IOpt{"offsetx", ind + 1},
			&la_.IOpt{"offsety", ind + 1})
		if !inverse {
			a = 1.0 / a
		}
		blas.ScalFloat(x, a, &la_.IOpt{"offset", ind}, &la_.IOpt{"n", m})
		ind += m
	}

	//if ! checkpnt.MinorEmpty() {
	//	checkpnt.Check("020scale2", minor)
	//}

	// For the 's' blocks, if inverse is 'N',
	//
	//     xk := vec( diag(l)^{-1/2} * mat(xk) * diag(k)^{-1/2}).
	//
	// If inverse is true,
	//
	//     xk := vec( diag(l)^{1/2} * mat(xk) * diag(k)^{1/2}).
	//
	// where l is kth block of lambda.
	//
	// We scale upper and lower triangular part of mat(xk) because the
	// inverse operation will be applied to nonsymmetric matrices.
	ind2 := ind
	sdims := dims.At("s")
	for k := 0; k < len(sdims); k++ {
		m := sdims[k]
		scaleF := func(v, x float64) float64 {
//.........這裏部分代碼省略.........

// Solves a pair of primal and dual cone programs  using custom KKT solver.
//
func ConeLpCustomKKT(c, G, h, A, b *matrix.FloatMatrix, dims *sets.DimensionSet,
	kktsolver KKTConeSolver, solopts *SolverOptions, primalstart,
	dualstart *sets.FloatMatrixSet) (sol *Solution, err error) {

	if c == nil || c.Cols() > 1 {
		err = errors.New("'c' must be matrix with 1 column")
		return
	}
	if h == nil {
		h = matrix.FloatZeros(0, 1)
	}
	if h.Cols() > 1 {
		err = errors.New("'h' must be matrix with 1 column")
		return
	}

	if dims == nil {
		dims = sets.NewDimensionSet("l", "q", "s")
		dims.Set("l", []int{h.Rows()})
	}
	cdim := dims.Sum("l", "q") + dims.SumSquared("s")
	cdim_pckd := dims.Sum("l", "q") + dims.SumPacked("s")
	//cdim_diag := dims.Sum("l", "q", "s")

	if G == nil {
		G = matrix.FloatZeros(0, c.Rows())
	}
	if !G.SizeMatch(cdim, c.Rows()) {
		estr := fmt.Sprintf("'G' must be of size (%d,%d)", cdim, c.Rows())
		err = errors.New(estr)
		return
	}

	// Check A and set defaults if it is nil
	if A == nil {
		// zeros rows reduces Gemv to vector products
		A = matrix.FloatZeros(0, c.Rows())
	}
	if A.Cols() != c.Rows() {
		estr := fmt.Sprintf("'A' must have %d columns", c.Rows())
		err = errors.New(estr)
		return
	}

	// Check b and set defaults if it is nil
	if b == nil {
		b = matrix.FloatZeros(0, 1)
	}
	if b.Cols() != 1 {
		estr := fmt.Sprintf("'b' must be a matrix with 1 column")
		err = errors.New(estr)
		return
	}
	if b.Rows() != A.Rows() {
		estr := fmt.Sprintf("'b' must have length %d", A.Rows())
		err = errors.New(estr)
		return
	}

	if b.Rows() > c.Rows() || b.Rows()+cdim_pckd < c.Rows() {
		err = errors.New("Rank(A) < p or Rank([G; A]) < n")
		return
	}

	mA := &matrixVarA{A}
	mG := &matrixVarG{G, dims}
	mc := &matrixVar{c}
	mb := &matrixVar{b}

	return conelp_problem(mc, mG, h, mA, mb, dims, kktsolver, solopts, primalstart, dualstart)
}

func conelp_solver(c MatrixVariable, G MatrixVarG, h *matrix.FloatMatrix,
	A MatrixVarA, b MatrixVariable, dims *sets.DimensionSet, kktsolver KKTConeSolverVar,
	solopts *SolverOptions, primalstart, dualstart *sets.FloatMatrixSet) (sol *Solution, err error) {

	err = nil
	const EXPON = 3
	const STEP = 0.99

	sol = &Solution{Unknown,
		nil,
		0.0, 0.0, 0.0, 0.0, 0.0,
		0.0, 0.0, 0.0, 0.0, 0.0, 0}

	var refinement int

	if solopts.Refinement > 0 {
		refinement = solopts.Refinement
	} else {
		refinement = 0
		if len(dims.At("q")) > 0 || len(dims.At("s")) > 0 {
			refinement = 1
		}
	}
	feasTolerance := FEASTOL
	absTolerance := ABSTOL
	relTolerance := RELTOL
	maxIter := MAXITERS
	if solopts.FeasTol > 0.0 {
		feasTolerance = solopts.FeasTol
	}
	if solopts.AbsTol > 0.0 {
		absTolerance = solopts.AbsTol
	}
	if solopts.RelTol > 0.0 {
		relTolerance = solopts.RelTol
	}
	if solopts.MaxIter > 0 {
		maxIter = solopts.MaxIter
	}
	if err = checkConeLpDimensions(dims); err != nil {
		return
	}

	cdim := dims.Sum("l", "q") + dims.SumSquared("s")
	//cdim_pckd := dims.Sum("l", "q") + dims.SumPacked("s")
	cdim_diag := dims.Sum("l", "q", "s")

	if h.Rows() != cdim {
		err = errors.New(fmt.Sprintf("'h' must be float matrix of size (%d,1)", cdim))
		return
	}

	// Data for kth 'q' constraint are found in rows indq[k]:indq[k+1] of G.
	indq := make([]int, 0)
	indq = append(indq, dims.At("l")[0])
	for _, k := range dims.At("q") {
		indq = append(indq, indq[len(indq)-1]+k)
	}

	// Data for kth 's' constraint are found in rows inds[k]:inds[k+1] of G.
	inds := make([]int, 0)
	inds = append(inds, indq[len(indq)-1])
	for _, k := range dims.At("s") {
		inds = append(inds, inds[len(inds)-1]+k*k)
	}

	Gf := func(x, y MatrixVariable, alpha, beta float64, trans la.Option) error {
		return G.Gf(x, y, alpha, beta, trans)
	}

	Af := func(x, y MatrixVariable, alpha, beta float64, trans la.Option) error {
		return A.Af(x, y, alpha, beta, trans)
	}

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

	// res() evaluates residual in 5x5 block KKT system
	//
	//     [ vx   ]    [ 0         ]   [ 0   A'  G'  c ] [ ux        ]
	//     [ vy   ]    [ 0         ]   [-A   0   0   b ] [ uy        ]
	//     [ vz   ] += [ W'*us     ] - [-G   0   0   h ] [ W^{-1}*uz ]
	//     [ vtau ]    [ dg*ukappa ]   [-c' -b' -h'  0 ] [ utau/dg   ]
	//
	//           vs += lmbda o (dz + ds)
	//       vkappa += lmbdg * (dtau + dkappa).
	ws3 := matrix.FloatZeros(cdim, 1)
	wz3 := matrix.FloatZeros(cdim, 1)
	checkpnt.AddMatrixVar("ws3", ws3)
	checkpnt.AddMatrixVar("wz3", wz3)

	//
//.........這裏部分代碼省略.........

// Solves a pair of primal and dual cone programs using custom KKT solver and constraint
// interfaces MatrixG and MatrixA
//
func ConeLpCustomMatrix(c *matrix.FloatMatrix, G MatrixG, h *matrix.FloatMatrix,
	A MatrixA, b *matrix.FloatMatrix, dims *sets.DimensionSet, kktsolver KKTConeSolver,
	solopts *SolverOptions, primalstart, dualstart *sets.FloatMatrixSet) (sol *Solution, err error) {

	err = nil

	if c == nil || c.Cols() > 1 {
		err = errors.New("'c' must be matrix with 1 column")
		return
	}
	if h == nil || h.Cols() > 1 {
		err = errors.New("'h' must be matrix with 1 column")
		return
	}

	if err = checkConeLpDimensions(dims); err != nil {
		return
	}

	cdim := dims.Sum("l", "q") + dims.SumSquared("s")
	cdim_pckd := dims.Sum("l", "q") + dims.SumPacked("s")
	//cdim_diag := dims.Sum("l", "q", "s")

	if h.Rows() != cdim {
		err = errors.New(fmt.Sprintf("'h' must be float matrix of size (%d,1)", cdim))
		return
	}

	// Data for kth 'q' constraint are found in rows indq[k]:indq[k+1] of G.
	indq := make([]int, 0)
	indq = append(indq, dims.At("l")[0])
	for _, k := range dims.At("q") {
		indq = append(indq, indq[len(indq)-1]+k)
	}

	// Data for kth 's' constraint are found in rows inds[k]:inds[k+1] of G.
	inds := make([]int, 0)
	inds = append(inds, indq[len(indq)-1])
	for _, k := range dims.At("s") {
		inds = append(inds, inds[len(inds)-1]+k*k)
	}

	// Check b and set defaults if it is nil
	if b == nil {
		b = matrix.FloatZeros(0, 1)
	}
	if b.Cols() != 1 {
		estr := fmt.Sprintf("'b' must be a matrix with 1 column")
		err = errors.New(estr)
		return
	}
	if b.Rows() > c.Rows() || b.Rows()+cdim_pckd < c.Rows() {
		err = errors.New("Rank(A) < p or Rank([G; A]) < n")
		return
	}

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

	var mA MatrixVarA
	var mG MatrixVarG
	if G == nil {
		mG = &matrixVarG{matrix.FloatZeros(0, c.Rows()), dims}
	} else {
		mG = &matrixIfG{G}
	}
	if A == nil {
		mA = &matrixVarA{matrix.FloatZeros(0, c.Rows())}
	} else {
		mA = &matrixIfA{A}
	}
	var mc = &matrixVar{c}
	var mb = &matrixVar{b}

	return conelp_problem(mc, mG, h, mA, mb, dims, kktsolver, solopts, primalstart, dualstart)
}