Golang cvx.SolverOptions類代碼示例

func main() {
	flag.Parse()
	if len(spPath) > 0 {
		checkpnt.Reset(spPath)
		checkpnt.Activate()
		checkpnt.Verbose(spVerbose)
		checkpnt.Format("%.17f")
	}

	adata := [][]float64{
		[]float64{0.3, -0.4, -0.2, -0.4, 1.3},
		[]float64{0.6, 1.2, -1.7, 0.3, -0.3},
		[]float64{-0.3, 0.0, 0.6, -1.2, -2.0}}

	A := matrix.FloatMatrixFromTable(adata, matrix.ColumnOrder)
	b := matrix.FloatVector([]float64{1.5, 0.0, -1.2, -0.7, 0.0})

	_, n := A.Size()
	N := n + 1 + n

	h := matrix.FloatZeros(N, 1)
	h.SetIndex(n, 1.0)

	I0 := matrix.FloatDiagonal(n, -1.0)
	I1 := matrix.FloatIdentity(n)
	G, _ := matrix.FloatMatrixStacked(matrix.StackDown, I0, matrix.FloatZeros(1, n), I1)

	At := A.Transpose()
	P := At.Times(A)
	q := At.Times(b).Scale(-1.0)

	dims := sets.NewDimensionSet("l", "q", "s")
	dims.Set("l", []int{n})
	dims.Set("q", []int{n + 1})

	var solopts cvx.SolverOptions
	solopts.MaxIter = 20
	solopts.ShowProgress = true
	if maxIter > 0 {
		solopts.MaxIter = maxIter
	}
	if len(solver) > 0 {
		solopts.KKTSolverName = solver
	}
	sol, err := cvx.ConeQp(P, q, G, h, nil, nil, dims, &solopts, nil)
	if err == nil {
		x := sol.Result.At("x")[0]
		s := sol.Result.At("s")[0]
		z := sol.Result.At("z")[0]
		fmt.Printf("Optimal\n")
		fmt.Printf("x=\n%v\n", x.ToString("%.9f"))
		fmt.Printf("s=\n%v\n", s.ToString("%.9f"))
		fmt.Printf("z=\n%v\n", z.ToString("%.9f"))
		check(x, s, z)
	}

}

func main() {

	flag.Parse()
	if len(spPath) > 0 {
		checkpnt.Reset(spPath)
		checkpnt.Activate()
		checkpnt.Verbose(spVerbose)
		checkpnt.Format("%.17f")
	}

	gdata := [][]float64{
		[]float64{16., 7., 24., -8., 8., -1., 0., -1., 0., 0., 7.,
			-5., 1., -5., 1., -7., 1., -7., -4.},
		[]float64{-14., 2., 7., -13., -18., 3., 0., 0., -1., 0., 3.,
			13., -6., 13., 12., -10., -6., -10., -28.},
		[]float64{5., 0., -15., 12., -6., 17., 0., 0., 0., -1., 9.,
			6., -6., 6., -7., -7., -6., -7., -11.}}

	hdata := []float64{-3., 5., 12., -2., -14., -13., 10., 0., 0., 0., 68.,
		-30., -19., -30., 99., 23., -19., 23., 10.}

	c := matrix.FloatVector([]float64{-6., -4., -5.})
	G := matrix.FloatMatrixFromTable(gdata)
	h := matrix.FloatVector(hdata)

	dims := sets.NewDimensionSet("l", "q", "s")
	dims.Set("l", []int{2})
	dims.Set("q", []int{4, 4})
	dims.Set("s", []int{3})

	var solopts cvx.SolverOptions
	solopts.MaxIter = 30
	solopts.ShowProgress = true
	if maxIter > 0 {
		solopts.MaxIter = maxIter
	}
	if len(solver) > 0 {
		solopts.KKTSolverName = solver
	}
	sol, err := cvx.ConeLp(c, G, h, nil, nil, dims, &solopts, nil, nil)
	if err == nil {
		x := sol.Result.At("x")[0]
		s := sol.Result.At("s")[0]
		z := sol.Result.At("z")[0]
		fmt.Printf("Optimal\n")
		fmt.Printf("x=\n%v\n", x.ToString("%.9f"))
		fmt.Printf("s=\n%v\n", s.ToString("%.9f"))
		fmt.Printf("z=\n%v\n", z.ToString("%.9f"))
		check(x, s, z)
	} else {
		fmt.Printf("status: %s\n", err)
	}
}

func main() {
	flag.Parse()

	gdata0 := [][]float64{
		[]float64{12., 13., 12.},
		[]float64{6., -3., -12.},
		[]float64{-5., -5., 6.}}

	gdata1 := [][]float64{
		[]float64{3., 3., -1., 1.},
		[]float64{-6., -6., -9., 19.},
		[]float64{10., -2., -2., -3.}}

	c := matrix.FloatVector([]float64{-2.0, 1.0, 5.0})
	g0 := matrix.FloatMatrixFromTable(gdata0, matrix.ColumnOrder)
	g1 := matrix.FloatMatrixFromTable(gdata1, matrix.ColumnOrder)
	Ghq := sets.FloatSetNew("Gq", "hq")
	Ghq.Append("Gq", g0, g1)

	h0 := matrix.FloatVector([]float64{-12.0, -3.0, -2.0})
	h1 := matrix.FloatVector([]float64{27.0, 0.0, 3.0, -42.0})
	Ghq.Append("hq", h0, h1)

	var Gl, hl, A, b *matrix.FloatMatrix = nil, nil, nil, nil
	var solopts cvx.SolverOptions
	solopts.MaxIter = 30
	solopts.ShowProgress = true
	if maxIter > -1 {
		solopts.MaxIter = maxIter
	}
	if len(solver) > 0 {
		solopts.KKTSolverName = solver
	}

	sol, err := cvx.Socp(c, Gl, hl, A, b, Ghq, &solopts, nil, nil)
	fmt.Printf("status: %v\n", err)
	if sol != nil && sol.Status == cvx.Optimal {
		x := sol.Result.At("x")[0]
		fmt.Printf("x=\n%v\n", x.ToString("%.9f"))
		for i, m := range sol.Result.At("sq") {
			fmt.Printf("sq[%d]=\n%v\n", i, m.ToString("%.9f"))
		}
		for i, m := range sol.Result.At("zq") {
			fmt.Printf("zq[%d]=\n%v\n", i, m.ToString("%.9f"))
		}
		sq0 := sol.Result.At("sq")[0]
		sq1 := sol.Result.At("sq")[1]
		zq0 := sol.Result.At("zq")[0]
		zq1 := sol.Result.At("zq")[1]
		check(x, sq0, sq1, zq0, zq1)
	}
}

func main() {
	flag.Parse()
	if len(spPath) > 0 {
		checkpnt.Reset(spPath)
		checkpnt.Activate()
		checkpnt.Verbose(spVerbose)
		checkpnt.Format("%.17f")
	}

	A := matrix.FloatNew(2, 3, []float64{1.0, -1.0, 0.0, 1.0, 0.0, 1.0})
	b := matrix.FloatNew(2, 1, []float64{1.0, 0.0})
	c := matrix.FloatNew(3, 1, []float64{0.0, 1.0, 0.0})
	G := matrix.FloatNew(1, 3, []float64{0.0, -1.0, 1.0})
	h := matrix.FloatNew(1, 1, []float64{0.0})
	//dims := sets.NewDimensionSet("l", "q", "s")
	//dims.Set("l", []int{1})

	fmt.Printf("A=\n%v\n", A)
	fmt.Printf("b=\n%v\n", b)
	fmt.Printf("G=\n%v\n", G)
	fmt.Printf("h=\n%v\n", h)
	fmt.Printf("c=\n%v\n", c)

	var solopts cvx.SolverOptions
	solopts.MaxIter = 30
	solopts.ShowProgress = true
	if maxIter > -1 {
		solopts.MaxIter = maxIter
	}
	if len(solver) > 0 {
		solopts.KKTSolverName = solver
	}
	sol, err := cvx.Lp(c, G, h, A, b, &solopts, nil, nil)
	if sol != nil && sol.Status == cvx.Optimal {
		x := sol.Result.At("x")[0]
		s := sol.Result.At("s")[0]
		z := sol.Result.At("z")[0]
		fmt.Printf("x=\n%v\n", x.ToString("%.9f"))
		fmt.Printf("s=\n%v\n", s.ToString("%.9f"))
		fmt.Printf("z=\n%v\n", z.ToString("%.9f"))
		check(x, s, z)
	} else {
		fmt.Printf("status: %v\n", err)
	}
}

func main() {
	flag.Parse()
	if len(spPath) > 0 {
		checkpnt.Reset(spPath)
		checkpnt.Activate()
		checkpnt.Verbose(spVerbose)
		checkpnt.Format("%.17f")
	}

	gdata := [][]float64{
		[]float64{2.0, 1.0, -1.0, 0.0},
		[]float64{1.0, 2.0, 0.0, -1.0}}

	c := matrix.FloatVector([]float64{-4.0, -5.0})
	G := matrix.FloatMatrixFromTable(gdata, matrix.ColumnOrder)
	h := matrix.FloatVector([]float64{3.0, 3.0, 0.0, 0.0})

	var solopts cvx.SolverOptions
	solopts.MaxIter = 30
	solopts.ShowProgress = true
	if maxIter > -1 {
		solopts.MaxIter = maxIter
	}
	if len(solver) > 0 {
		solopts.KKTSolverName = solver
	}
	sol, err := cvx.Lp(c, G, h, nil, nil, &solopts, nil, nil)
	if sol != nil && sol.Status == cvx.Optimal {
		x := sol.Result.At("x")[0]
		s := sol.Result.At("s")[0]
		z := sol.Result.At("z")[0]
		fmt.Printf("x=\n%v\n", x.ToString("%.9f"))
		fmt.Printf("s=\n%v\n", s.ToString("%.9f"))
		fmt.Printf("z=\n%v\n", z.ToString("%.9f"))
		check(x, s, z)
	} else {
		fmt.Printf("status: %v\n", err)
	}
}

func acenter() *matrix.FloatMatrix {

	F := &acenterProg{3, 1}

	gdata := [][]float64{
		[]float64{0., -1., 0., 0., -21., -11., 0., -11., 10., 8., 0., 8., 5.},
		[]float64{0., 0., -1., 0., 0., 10., 16., 10., -10., -10., 16., -10., 3.},
		[]float64{0., 0., 0., -1., -5., 2., -17., 2., -6., 8., -17., -7., 6.}}

	G := matrix.FloatMatrixFromTable(gdata)
	h := matrix.FloatVector(
		[]float64{1.0, 0.0, 0.0, 0.0, 20., 10., 40., 10., 80., 10., 40., 10., 15.})

	var solopts cvx.SolverOptions
	solopts.MaxIter = 40
	solopts.ShowProgress = true
	if maxIter > -1 {
		solopts.MaxIter = maxIter
	}
	if len(solver) > 0 {
		solopts.KKTSolverName = solver
	}

	dims := sets.NewDimensionSet("l", "q", "s")
	dims.Set("l", []int{0})
	dims.Set("q", []int{4})
	dims.Set("s", []int{3})

	var err error
	var sol *cvx.Solution

	sol, err = cvx.Cp(F, G, h, nil, nil, dims, &solopts)
	if err == nil && sol.Status == cvx.Optimal {
		return sol.Result.At("x")[0]
	} else {
		fmt.Printf("result: %v\n", err)
	}
	return nil
}

func main() {

	Sdata := [][]float64{
		[]float64{4e-2, 6e-3, -4e-3, 0.0},
		[]float64{6e-3, 1e-2, 0.0, 0.0},
		[]float64{-4e-3, 0.0, 2.5e-3, 0.0},
		[]float64{0.0, 0.0, 0.0, 0.0}}

	pbar := matrix.FloatVector([]float64{.12, .10, .07, .03})
	S := matrix.FloatMatrixFromTable(Sdata)
	n := pbar.Rows()
	G := matrix.FloatDiagonal(n, -1.0)
	h := matrix.FloatZeros(n, 1)
	A := matrix.FloatWithValue(1, n, 1.0)
	b := matrix.FloatNew(1, 1, []float64{1.0})

	var solopts cvx.SolverOptions
	solopts.MaxIter = 30
	solopts.ShowProgress = true

	mu := 1.0
	Smu := matrix.Scale(S, mu)
	pbarNeg := matrix.Scale(pbar, -1.0)
	fmt.Printf("Smu=\n%v\n", Smu.String())
	fmt.Printf("-pbar=\n%v\n", pbarNeg.String())

	sol, err := cvx.Qp(Smu, pbarNeg, G, h, A, b, &solopts, nil)

	fmt.Printf("status: %v\n", err)
	if sol != nil && sol.Status == cvx.Optimal {
		x := sol.Result.At("x")[0]
		ret := blas.DotFloat(x, pbar)
		risk := math.Sqrt(blas.DotFloat(x, S.Times(x)))
		fmt.Printf("ret=%.3f, risk=%.3f\n", ret, risk)
		fmt.Printf("x=\n%v\n", x)
	}
}

func main() {
	flag.Parse()

	aflr := 1000.0
	awall := 100.0
	alpha := 0.5
	beta := 2.0
	gamma := 0.5
	delta := 2.0

	fdata := [][]float64{
		[]float64{-1.0, 1.0, 1.0, 0.0, -1.0, 1.0, 0.0, 0.0},
		[]float64{-1.0, 1.0, 0.0, 1.0, 1.0, -1.0, 1.0, -1.0},
		[]float64{-1.0, 0.0, 1.0, 1.0, 0.0, 0.0, -1.0, 1.0}}

	gdata := []float64{1.0, 2.0 / awall, 2.0 / awall, 1.0 / aflr, alpha, 1.0 / beta, gamma, 1.0 / delta}

	g := matrix.FloatNew(8, 1, gdata).Log()
	F := matrix.FloatMatrixFromTable(fdata)
	K := []int{1, 2, 1, 1, 1, 1, 1}

	var solopts cvx.SolverOptions
	solopts.MaxIter = 40
	if maxIter > 0 {
		solopts.MaxIter = maxIter
	}
	if len(spPath) > 0 {
		checkpnt.Reset(spPath)
		checkpnt.Activate()
		checkpnt.Verbose(spVerbose)
		checkpnt.Format("%.7f")
	}
	solopts.ShowProgress = true
	if maxIter > 0 {
		solopts.MaxIter = maxIter
	}
	if len(solver) > 0 {
		solopts.KKTSolverName = solver
	}
	sol, err := cvx.Gp(K, F, g, nil, nil, nil, nil, &solopts)
	if sol != nil && sol.Status == cvx.Optimal {
		x := sol.Result.At("x")[0]
		r := matrix.Exp(x)
		h := r.GetIndex(0)
		w := r.GetIndex(1)
		d := r.GetIndex(2)
		fmt.Printf("x=\n%v\n", x.ToString("%.9f"))
		fmt.Printf("\n h = %f,  w = %f, d = %f.\n", h, w, d)
		check(x)
	} else {
		fmt.Printf("status: %v\n", err)
	}
}

func main() {
	flag.Parse()
	if len(spPath) > 0 {
		checkpnt.Reset(spPath)
		checkpnt.Activate()
		checkpnt.Verbose(spVerbose)
		checkpnt.Format("%.17f")
	}

	gdata0 := [][]float64{
		[]float64{-7., -11., -11., 3.},
		[]float64{7., -18., -18., 8.},
		[]float64{-2., -8., -8., 1.}}

	gdata1 := [][]float64{
		[]float64{-21., -11., 0., -11., 10., 8., 0., 8., 5.},
		[]float64{0., 10., 16., 10., -10., -10., 16., -10., 3.},
		[]float64{-5., 2., -17., 2., -6., 8., -17., -7., 6.}}

	hdata0 := [][]float64{
		[]float64{33., -9.},
		[]float64{-9., 26.}}

	hdata1 := [][]float64{
		[]float64{14., 9., 40.},
		[]float64{9., 91., 10.},
		[]float64{40., 10., 15.}}

	g0 := matrix.FloatMatrixFromTable(gdata0, matrix.ColumnOrder)
	g1 := matrix.FloatMatrixFromTable(gdata1, matrix.ColumnOrder)
	Ghs := sets.FloatSetNew("Gs", "hs")
	Ghs.Append("Gs", g0, g1)

	h0 := matrix.FloatMatrixFromTable(hdata0, matrix.ColumnOrder)
	h1 := matrix.FloatMatrixFromTable(hdata1, matrix.ColumnOrder)
	Ghs.Append("hs", h0, h1)

	c := matrix.FloatVector([]float64{1.0, -1.0, 1.0})

	var Gs, hs, A, b *matrix.FloatMatrix = nil, nil, nil, nil
	var solopts cvx.SolverOptions
	solopts.MaxIter = 30
	solopts.ShowProgress = true
	if maxIter > -1 {
		solopts.MaxIter = maxIter
	}
	if len(solver) > 0 {
		solopts.KKTSolverName = solver
	}

	sol, err := cvx.Sdp(c, Gs, hs, A, b, Ghs, &solopts, nil, nil)
	if sol != nil && sol.Status == cvx.Optimal {
		x := sol.Result.At("x")[0]
		fmt.Printf("x=\n%v\n", x.ToString("%.9f"))
		for i, m := range sol.Result.At("zs") {
			fmt.Printf("zs[%d]=\n%v\n", i, m.ToString("%.9f"))
		}
		ss0 := sol.Result.At("ss")[0]
		ss1 := sol.Result.At("ss")[1]
		zs0 := sol.Result.At("zs")[0]
		zs1 := sol.Result.At("zs")[1]
		check(x, ss0, ss1, zs0, zs1)
	} else {
		fmt.Printf("status: %v\n", err)
	}
	checkpnt.Report()
}

func floorplan(Amin *matrix.FloatMatrix) *matrix.FloatMatrix {
	rho := 1.0
	gamma := 5.0

	c := matrix.FloatZeros(22, 1)
	c.SetAtColumnArray(0, []int{0, 1}, []float64{1.0, 1.0})

	G := matrix.FloatZeros(26, 22)
	h := matrix.FloatZeros(26, 1)

	// -x1 <= 0
	G.SetAt(0, 2, -1.0)

	// -x2 <= 0
	G.SetAt(1, 3, -1.0)

	// -x4 <= 0
	G.SetAt(2, 5, -1.0)

	// x1 - x3 + w1 <= -rho
	G.SetAtRowArray(3, []int{2, 4, 12}, []float64{1.0, -1.0, 1.0})
	h.SetAt(3, 0, -rho)

	// x2 - x3 + w2 <= -rho
	G.SetAtRowArray(4, []int{3, 4, 13}, []float64{1.0, -1.0, 1.0})
	h.SetAt(4, 0, -rho)

	// x3 - x5 + w3 <= -rho
	G.SetAtRowArray(5, []int{4, 6, 14}, []float64{1.0, -1.0, 1.0})
	h.SetAt(5, 0, -rho)

	// x4 - x5 + w4 <= -rho
	G.SetAtRowArray(6, []int{5, 6, 15}, []float64{1.0, -1.0, 1.0})
	h.SetAt(6, 0, -rho)

	// -W + x5 + w5 <= 0
	G.SetAtRowArray(7, []int{0, 6, 16}, []float64{-1.0, 1.0, 1.0})

	// -y2 <= 0
	G.SetAt(8, 8, -1.0)

	// -y3 <= 0
	G.SetAt(9, 9, -1.0)

	// -y5 <= 0
	G.SetAt(10, 11, -1.0)

	// -y1 + y2 + h2 <= -rho
	G.SetAtRowArray(11, []int{7, 8, 18}, []float64{-1.0, 1.0, 1.0})
	h.SetAt(11, 0, -rho)

	// y1 - y4 + h1 <= -rho
	G.SetAtRowArray(12, []int{7, 10, 17}, []float64{1.0, -1.0, 1.0})
	h.SetAt(12, 0, -rho)

	// y3 - y4 + h3 <= -rho
	G.SetAtRowArray(13, []int{9, 10, 19}, []float64{1.0, -1.0, 1.0})
	h.SetAt(13, 0, -rho)

	// -H + y4 + h4 <= 0
	G.SetAtRowArray(14, []int{1, 10, 20}, []float64{-1.0, 1.0, 1.0})

	// -H + y5 + h5 <= 0
	G.SetAtRowArray(15, []int{1, 11, 21}, []float64{-1.0, 1.0, 1.0})

	// -w1 + h1/gamma <= 0
	G.SetAtRowArray(16, []int{12, 17}, []float64{-1.0, 1.0 / gamma})

	// w1 - gamma * h1 <= 0
	G.SetAtRowArray(17, []int{12, 17}, []float64{1.0, -gamma})

	// -w2 + h2/gamma <= 0
	G.SetAtRowArray(18, []int{13, 18}, []float64{-1.0, 1.0 / gamma})

	//  w2 - gamma * h2 <= 0
	G.SetAtRowArray(19, []int{13, 18}, []float64{1.0, -gamma})

	// -w3 + h3/gamma <= 0
	G.SetAtRowArray(20, []int{14, 18}, []float64{-1.0, 1.0 / gamma})

	//  w3 - gamma * h3 <= 0
	G.SetAtRowArray(21, []int{14, 19}, []float64{1.0, -gamma})

	// -w4  + h4/gamma <= 0
	G.SetAtRowArray(22, []int{15, 19}, []float64{-1.0, 1.0 / gamma})

	//  w4 - gamma * h4 <= 0
	G.SetAtRowArray(23, []int{15, 20}, []float64{1.0, -gamma})

	// -w5 + h5/gamma <= 0
	G.SetAtRowArray(24, []int{16, 21}, []float64{-1.0, 1.0 / gamma})

	//  w5 - gamma * h5 <= 0.0
	G.SetAtRowArray(25, []int{16, 21}, []float64{1.0, -gamma})

	F := newFloorPlan(Amin)

	var dims *sets.DimensionSet = nil
	var solopts cvx.SolverOptions
	solopts.MaxIter = 50
//.........這裏部分代碼省略.........

func main() {
	flag.Parse()

	m := len(udata)
	nvars := 2 * m
	u := matrix.FloatVector(udata[:m])
	y := matrix.FloatVector(ydata[:m])

	// minimize    (1/2) * || yhat - y ||_2^2
	// subject to  yhat[j] >= yhat[i] + g[i]' * (u[j] - u[i]), j, i = 0,...,m-1
	//
	// Variables  yhat (m), g (m).

	P := matrix.FloatZeros(nvars, nvars)
	// set m first diagonal indexes to 1.0
	//P.SetIndexes(1.0, matrix.DiagonalIndexes(P)[:m]...)
	P.Diag().SubMatrix(0, 0, 1, m).SetIndexes(1.0)
	q := matrix.FloatZeros(nvars, 1)
	q.SubMatrix(0, 0, y.NumElements(), 1).Plus(matrix.Scale(y, -1.0))

	// m blocks (i = 0,...,m-1) of linear inequalities
	//
	//     yhat[i] + g[i]' * (u[j] - u[i]) <= yhat[j], j = 0,...,m-1.

	G := matrix.FloatZeros(m*m, nvars)
	I := matrix.FloatDiagonal(m, 1.0)

	for i := 0; i < m; i++ {
		// coefficients of yhat[i] (column i)
		//G.Set(1.0, matrix.ColumnIndexes(G, i)[i*m:(i+1)*m]...)
		column(G, i).SetIndexes(1.0)

		// coefficients of gi[i] (column i, rows i*m ... (i+1)*m)
		//rows := matrix.Indexes(i*m, (i+1)*m)
		//G.SetAtColumnArray(m+i, rows, matrix.Add(u, -u.GetIndex(i)).FloatArray())

		// coefficients of gi[i] (column i, rows i*m ... (i+1)*m)
		// from column m+i staring at row i*m select m rows and one column
		G.SubMatrix(i*m, m+i, m, 1).Plus(matrix.Add(u, -u.GetIndex(i)))

		// coeffients of yhat[i]) from rows i*m ... (i+1)*m, cols 0 ... m
		//G.SetSubMatrix(i*m, 0, matrix.Minus(G.GetSubMatrix(i*m, 0, m, m), I))
		G.SubMatrix(i*m, 0, m, m).Minus(I)
	}

	h := matrix.FloatZeros(m*m, 1)
	var A, b *matrix.FloatMatrix = nil, nil
	var solopts cvx.SolverOptions
	solopts.ShowProgress = true
	solopts.KKTSolverName = solver

	sol, err := cvx.Qp(P, q, G, h, A, b, &solopts, nil)
	if err != nil {
		fmt.Printf("error: %v\n", err)
		return
	}
	if sol != nil && sol.Status != cvx.Optimal {
		fmt.Printf("status not optimal\n")
		return
	}
	x := sol.Result.At("x")[0]
	//yhat := matrix.FloatVector(x.FloatArray()[:m])
	//g := matrix.FloatVector(x.FloatArray()[m:])
	yhat := x.SubMatrix(0, 0, m, 1).Copy()
	g := x.SubMatrix(m, 0).Copy()

	rangeFunc := func(n int) []float64 {
		r := make([]float64, 0)
		for i := 0; i < n; i++ {
			r = append(r, float64(i)*2.2/float64(n))
		}
		return r
	}
	ts := rangeFunc(1000)
	fitFunc := func(points []float64) []float64 {
		res := make([]float64, len(points))
		for k, t := range points {
			res[k] = matrix.Plus(yhat, matrix.Mul(g, matrix.Scale(u, -1.0).Add(t))).Max()
		}
		return res
	}
	fs := fitFunc(ts)
	plotData("cvxfit.png", u.FloatArray(), y.FloatArray(), ts, fs)
}

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

			// x[n:] += div(z[:n], d1) + div(z[n:2*n], d2)
			z_res = matrix.Plus(matrix.Div(z_n, d1), matrix.Div(z_2n, d2))
			x.SubMatrix(n, 0, z_res.NumElements(), 1).Plus(z_res)
			checkpnt.Check("15-f", minor)

			// Solve for x[:n]:
			//
			//    S*x[:n] = x[:n] - (W1**2 - W2**2)(W1**2 + W2**2)^-1 * x[n:]

			// w1 = (d1**2 - d2**2), w2 = (d1**2 + d2**2)
			w1 := matrix.Minus(matrix.Mul(d1, d1), matrix.Mul(d2, d2))
			w2 := matrix.Plus(matrix.Mul(d1, d1), matrix.Mul(d2, d2))

			// x[:n] += -mul( div(w1, w2), x[n:])
			x_n := x.SubMatrix(n, 0)
			x_val := matrix.Mul(matrix.Div(w1, w2), x_n).Scale(-1.0)
			x.SubMatrix(0, 0, n, 1).Plus(x_val)
			checkpnt.Check("25-f", minor)

			// Solve for x[n:]:
			//
			//    (d1**-2 + d2**-2) * x[n:] = x[n:] + (d1**-2 - d2**-2)*x[:n]

			err = lapack.Potrs(S, x)
			if err != nil {
				fmt.Printf("Potrs error: %s\n", err)
			}
			checkpnt.Check("30-f", minor)

			// Solve for x[n:]:
			//
			//    (d1**-2 + d2**-2) * x[n:] = x[n:] + (d1**-2 - d2**-2)*x[:n]

			// w1 = (d1**-2 - d2**-2), w2 = (d1**-2 + d2**-2)
			w1 = matrix.Minus(matrix.Mul(d1, d1).Inv(), matrix.Mul(d2, d2).Inv())
			w2 = matrix.Plus(matrix.Mul(d1, d1).Inv(), matrix.Mul(d2, d2).Inv())
			x_n = x.SubMatrix(0, 0, n, 1)

			// x[n:] += mul( d1**-2 - d2**-2, x[:n])
			x_val = matrix.Mul(w1, x_n)
			x.SubMatrix(n, 0, x_val.NumElements(), 1).Plus(x_val)
			checkpnt.Check("35-f", minor)

			// x[n:] = div( x[n:], d1**-2 + d2**-2)
			x_n = x.SubMatrix(n, 0)
			x_val = matrix.Div(x_n, w2)
			x.SubMatrix(n, 0, x_val.NumElements(), 1).Set(x_val)
			checkpnt.Check("40-f", minor)

			// x_n = x[:n], x-2n = x[n:2*n]
			x_n = x.SubMatrix(0, 0, n, 1)
			x_2n := x.SubMatrix(n, 0, n, 1)

			// z := z + W^-T * G*x
			// z[:n] += div( x[:n] - x[n:2*n], d1)
			x_val = matrix.Div(matrix.Minus(x_n, x_2n), d1)
			z.SubMatrix(0, 0, n, 1).Plus(x_val)
			checkpnt.Check("44-f", minor)

			// z[n:2*n] += div( -x[:n] - x[n:2*n], d2)
			x_val = matrix.Div(matrix.Plus(x_n, x_2n).Scale(-1.0), d2)
			z.SubMatrix(n, 0, n, 1).Plus(x_val)
			checkpnt.Check("48-f", minor)

			// z[2*n:] += As*x[:n]
			x_val = matrix.Times(As, x_n)
			z.SubMatrix(2*n, 0, x_val.NumElements(), 1).Plus(x_val)

			checkpnt.Check("50-f", minor)

			return nil
		}
		return
	}

	// matrix(n*[0.0] + n*[1.0])
	c := matrix.FloatZeros(2*n, 1)
	c.SubMatrix(n, 0).SetIndexes(1.0)

	h := matrix.FloatZeros(2*n+m+1, 1)
	h.SetIndexes(1.0, 2*n)
	// h[2*n+1:] = -b
	h.SubMatrix(2*n+1, 0).Plus(b).Scale(-1.0)
	G := &matrixFs{A}

	dims := sets.DSetNew("l", "q", "s")
	dims.Set("l", []int{2 * n})
	dims.Set("q", []int{m + 1})

	var solopts cvx.SolverOptions
	solopts.ShowProgress = true
	if maxIter > 0 {
		solopts.MaxIter = maxIter
	}
	if len(solver) > 0 {
		solopts.KKTSolverName = solver
	}
	return cvx.ConeLpCustomMatrix(c, G, h, nil, nil, dims, Fkkt, &solopts, nil, nil)
}

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

	Fkkt := func(W *sets.FloatMatrixSet) (cvx.KKTFunc, error) {

		//    Solve
		//                  -diag(z)                           = bx
		//        -diag(x) - inv(rti*rti') * z * inv(rti*rti') = bs
		//
		//    On entry, x and z contain bx and bs.
		//    On exit, they contain the solution, with z scaled
		//    (inv(rti)'*z*inv(rti) is returned instead of z).
		//
		//    We first solve
		//
		//        ((rti*rti') .* (rti*rti')) * x = bx - diag(t*bs*t)
		//
		//    and take z  = -rti' * (diag(x) + bs) * rti.

		var err error = nil
		rti := W.At("rti")[0]

		// t = rti*rti' as a nonsymmetric matrix.
		t := matrix.FloatZeros(n, n)
		err = blas.GemmFloat(rti, rti, t, 1.0, 0.0, linalg.OptTransB)
		if err != nil {
			return nil, err
		}

		// Cholesky factorization of tsq = t.*t.
		tsq := matrix.Mul(t, t)
		err = lapack.Potrf(tsq)
		if err != nil {
			return nil, err
		}

		f := func(x, y, z *matrix.FloatMatrix) (err error) {
			// tbst := t * zs * t = t * bs * t
			tbst := z.Copy()
			matrix.Reshape(tbst, n, n)
			cngrnc(t, tbst, 1.0)

			// x := x - diag(tbst) = bx - diag(rti*rti' * bs * rti*rti')
			diag := tbst.Diag().Transpose()
			x.Minus(diag)

			// x := (t.*t)^{-1} * x = (t.*t)^{-1} * (bx - diag(t*bs*t))
			err = lapack.Potrs(tsq, x)
			if err != nil {
				fmt.Printf("Fkkt.f.Potrs: %v\n", err)
			}

			// z := z + diag(x) = bs + diag(x)
			// z, x are really column vectors here
			z.AddIndexes(matrix.MakeIndexSet(0, n*n, n+1), x.FloatArray())

			// z := -rti' * z * rti = -rti' * (diag(x) + bs) * rti
			cngrnc(rti, z, -1.0)
			return nil
		}
		return f, nil
	}

	c := matrix.FloatWithValue(n, 1, 1.0)

	// initial feasible x: x = 1.0 - min(lmbda(w))
	lmbda := matrix.FloatZeros(n, 1)
	wp := w.Copy()
	lapack.Syevx(wp, lmbda, nil, 0.0, nil, []int{1, 1}, linalg.OptRangeInt)
	x0 := matrix.FloatZeros(n, 1).Add(-lmbda.GetAt(0, 0) + 1.0)
	s0 := w.Copy()
	// Diag() return a row vector, x0 is column vector
	s0.Diag().Plus(x0.Transpose())
	matrix.Reshape(s0, n*n, 1)

	// initial feasible z is identity
	z0 := matrix.FloatIdentity(n)
	matrix.Reshape(z0, n*n, 1)

	dims := sets.DSetNew("l", "q", "s")
	dims.Set("s", []int{n})

	primalstart := sets.FloatSetNew("x", "s")
	dualstart := sets.FloatSetNew("z")
	primalstart.Set("x", x0)
	primalstart.Set("s", s0)
	dualstart.Set("z", z0)

	var solopts cvx.SolverOptions
	solopts.ShowProgress = true
	if maxIter > 0 {
		solopts.MaxIter = maxIter
	}
	if len(solver) > 0 {
		solopts.KKTSolverName = solver
	}
	h := w.Copy()
	matrix.Reshape(h, h.NumElements(), 1)
	return cvx.ConeLpCustomMatrix(c, G, h, nil, nil, dims, Fkkt, &solopts, primalstart, dualstart)
}