Golang FloatMatrix.Cols方法代码示例

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

//    Solves a quadratic program
//
//        minimize    (1/2)*x'*P*x + q'*x
//        subject to  G*x <= h
//                    A*x = b.
//
//
//    Input arguments.
//
//        P is a n x n float matrix with the lower triangular part of P stored
//        in the lower triangle.  Must be positive semidefinite.
//
//        q is an n x 1 matrix.
//
//        G is an m x n matrix or nil.
//
//        h is an m x 1 matrix or nil.
//
//        A is a p x n matrix or nil.
//
//        b is a p x 1 matrix or nil.
//
//        The default values for G, h, A and b are empty matrices with zero rows.
//
//
func Qp(P, q, G, h, A, b *matrix.FloatMatrix, solopts *SolverOptions, initvals *FloatMatrixSet) (sol *Solution, err error) {

	sol = nil
	if P == nil || P.Rows() != P.Cols() {
		err = errors.New("'P' must a non-nil square matrix")
		return
	}
	if q == nil {
		err = errors.New("'q' must a non-nil matrix")
		return
	}
	if q.Rows() != P.Rows() || q.Cols() > 1 {
		err = errors.New(fmt.Sprintf("'q' must be matrix of size (%d,1)", P.Rows()))
		return
	}
	if G == nil {
		G = matrix.FloatZeros(0, P.Rows())
	}
	if G.Cols() != P.Rows() {
		err = errors.New(fmt.Sprintf("'G' must be matrix of %d columns", P.Rows()))
		return
	}
	if h == nil {
		h = matrix.FloatZeros(G.Rows(), 1)
	}
	if h.Rows() != G.Rows() || h.Cols() > 1 {
		err = errors.New(fmt.Sprintf("'h' must be matrix of size (%d,1)", G.Rows()))
		return
	}
	if A == nil {
		A = matrix.FloatZeros(0, P.Rows())
	}
	if A.Cols() != P.Rows() {
		err = errors.New(fmt.Sprintf("'A' must be matrix of %d columns", P.Rows()))
		return
	}
	if b == nil {
		b = matrix.FloatZeros(A.Rows(), 1)
	}
	if b.Rows() != A.Rows() {
		err = errors.New(fmt.Sprintf("'b' must be matrix of size (%d,1)", A.Rows()))
		return
	}
	return ConeQp(P, q, G, h, A, b, nil, solopts, initvals)
}

//    Solves a pair of primal and dual LPs
//
//        minimize    c'*x
//        subject to  G*x + s = h
//                    A*x = b
//                    s >= 0
//
//        maximize    -h'*z - b'*y
//        subject to  G'*z + A'*y + c = 0
//                    z >= 0.
//
func Lp(c, G, h, A, b *matrix.FloatMatrix, solopts *SolverOptions, primalstart, dualstart *FloatMatrixSet) (sol *Solution, err error) {

	if c == nil {
		err = errors.New("'c' must a column matrix")
		return
	}
	n := c.Rows()
	if n < 1 {
		err = errors.New("Number of variables must be at least 1")
		return
	}
	if G == nil || G.Cols() != n {
		err = errors.New(fmt.Sprintf("'G' must be matrix with %d columns", n))
		return
	}
	m := G.Rows()
	if h == nil || !h.SizeMatch(m, 1) {
		err = errors.New(fmt.Sprintf("'h' must be matrix of size (%d,1)", m))
		return
	}
	if A == nil {
		A = matrix.FloatZeros(0, n)
	}
	if A.Cols() != n {
		err = errors.New(fmt.Sprintf("'A' must be matrix with %d columns", n))
		return
	}
	p := A.Rows()
	if b == nil {
		b = matrix.FloatZeros(0, 1)
	}
	if !b.SizeMatch(p, 1) {
		err = errors.New(fmt.Sprintf("'b' must be matrix of size (%d,1)", p))
		return
	}
	dims := DSetNew("l", "q", "s")
	dims.Set("l", []int{m})

	return ConeLp(c, G, h, A, b, dims, solopts, primalstart, dualstart)
}

//    Solves a pair of primal and dual cone programs
//
//        minimize    c'*x
//        subject to  G*x + s = h
//                    A*x = b
//                    s >= 0
//
//        maximize    -h'*z - b'*y
//        subject to  G'*z + A'*y + c = 0
//                    z >= 0.
//
//    The inequalities are with respect to a cone C defined as the Cartesian
//    product of N + M + 1 cones:
//
//        C = C_0 x C_1 x .... x C_N x C_{N+1} x ... x C_{N+M}.
//
//    The first cone C_0 is the nonnegative orthant of dimension ml.
//    The next N cones are second order cones of dimension mq[0], ...,
//    mq[N-1].  The second order cone of dimension m is defined as
//
//        { (u0, u1) in R x R^{m-1} | u0 >= ||u1||_2 }.
//
//    The next M cones are positive semidefinite cones of order ms[0], ...,
//    ms[M-1] >= 0.
//
func ConeLp(c, G, h, A, b *matrix.FloatMatrix, dims *DimensionSet, solopts *SolverOptions, primalstart, dualstart *FloatMatrixSet) (sol *Solution, err error) {

	err = nil
	const EXPON = 3
	const STEP = 0.99

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

	//var primalstart *FloatMatrixSet = nil
	//var dualstart *FloatMatrixSet = nil
	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
	if solopts.FeasTol > 0.0 {
		feasTolerance = solopts.FeasTol
	}
	if solopts.AbsTol > 0.0 {
		absTolerance = solopts.AbsTol
	}
	if solopts.RelTol > 0.0 {
		relTolerance = solopts.RelTol
	}

	solvername := solopts.KKTSolverName
	if len(solvername) == 0 {
		if dims != nil && (len(dims.At("q")) > 0 || len(dims.At("s")) > 0) {
			solvername = "qr"
		} else {
			solvername = "chol2"
		}
	}

	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 dims == nil {
		dims = DSetNew("l", "q", "s")
		dims.Set("l", []int{h.Rows()})
	}
	if err = checkConeLpDimensions(dims); err != nil {
		return
	}

	cdim := dims.Sum("l", "q") + dims.SumSquared("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, 100)
	indq = append(indq, dims.At("l")[0])
	for _, k := range dims.At("q") {
		indq = append(indq, indq[len(indq)-1]+k)
//.........这里部分代码省略.........

/*
   Applies Nesterov-Todd scaling or its inverse.

   Computes

        x := W*x        (trans is false 'N', inverse = false 'N')
        x := W^T*x      (trans is true  'T', inverse = false 'N')
        x := W^{-1}*x   (trans is false 'N', inverse = true  'T')
        x := W^{-T}*x   (trans is true  'T', inverse = true  'T').

   x is a dense float matrix.

   W is a MatrixSet with entries:

   - W['dnl']: positive vector
   - W['dnli']: componentwise inverse of W['dnl']
   - W['d']: positive vector
   - W['di']: componentwise inverse of W['d']
   - W['v']: lists of 2nd order cone vectors with unit hyperbolic norms
   - W['beta']: list of positive numbers
   - W['r']: list of square matrices
   - W['rti']: list of square matrices.  rti[k] is the inverse transpose
     of r[k].

   The 'dnl' and 'dnli' entries are optional, and only present when the
   function is called from the nonlinear solver.
*/
func scale(x *matrix.FloatMatrix, W *FloatMatrixSet, trans, inverse bool) (err error) {
	/*DEBUGGED*/
	var wl []*matrix.FloatMatrix
	var w *matrix.FloatMatrix
	ind := 0
	err = nil

	// Scaling for nonlinear component xk is xk := dnl .* xk; inverse
	// scaling is xk ./ dnl = dnli .* xk, where dnl = W['dnl'],
	// dnli = W['dnli'].

	if wl = W.At("dnl"); wl != nil {
		if inverse {
			w = W.At("dnli")[0]
		} else {
			w = W.At("dnl")[0]
		}
		for k := 0; k < x.Cols(); k++ {
			err = blas.TbmvFloat(w, x, &la_.IOpt{"n", w.Rows()}, &la_.IOpt{"k", 0},
				&la_.IOpt{"lda", 1}, &la_.IOpt{"offsetx", k * x.Rows()})
			if err != nil {
				return
			}
		}
		ind += w.Rows()
	}

	// Scaling for linear 'l' component xk is xk := d .* xk; inverse
	// scaling is xk ./ d = di .* xk, where d = W['d'], di = W['di'].

	if inverse {
		w = W.At("di")[0]
	} else {
		w = W.At("d")[0]
	}

	for k := 0; k < x.Cols(); k++ {
		err = blas.TbmvFloat(w, x, &la_.IOpt{"n", w.Rows()}, &la_.IOpt{"k", 0},
			&la_.IOpt{"lda", 1}, &la_.IOpt{"offsetx", k*x.Rows() + ind})
		if err != nil {
			return
		}
	}
	ind += w.Rows()

	// Scaling for 'q' component is
	//
	//    xk := beta * (2*v*v' - J) * xk
	//        = beta * (2*v*(xk'*v)' - J*xk)
	//
	// where beta = W['beta'][k], v = W['v'][k], J = [1, 0; 0, -I].
	//
	//Inverse scaling is
	//
	//    xk := 1/beta * (2*J*v*v'*J - J) * xk
	//        = 1/beta * (-J) * (2*v*((-J*xk)'*v)' + xk).
	//wf := matrix.FloatZeros(x.Cols(), 1)
	w = matrix.FloatZeros(x.Cols(), 1)
	for k, v := range W.At("v") {
		m := v.Rows()
		if inverse {
			blas.ScalFloat(x, -1.0, &la_.IOpt{"offset", ind}, &la_.IOpt{"inc", x.Rows()})
		}
		err = blas.GemvFloat(x, v, w, 1.0, 0.0, la_.OptTrans, &la_.IOpt{"m", m},
			&la_.IOpt{"n", x.Cols()}, &la_.IOpt{"offsetA", ind},
			&la_.IOpt{"lda", x.Rows()})
		if err != nil {
			return
		}

		err = blas.ScalFloat(x, -1.0, &la_.IOpt{"offset", ind}, &la_.IOpt{"inc", x.Rows()})
		if err != nil {
			return
//.........这里部分代码省略.........

//    Solves a pair of primal and dual SDPs
//
//        minimize    c'*x
//        subject to  Gl*x + sl = hl
//                    mat(Gs[k]*x) + ss[k] = hs[k], k = 0, ..., N-1
//                    A*x = b
//                    sl >= 0,  ss[k] >= 0, k = 0, ..., N-1
//
//        maximize    -hl'*z - sum_k trace(hs[k]*zs[k]) - b'*y
//        subject to  Gl'*zl + sum_k Gs[k]'*vec(zs[k]) + A'*y + c = 0
//                    zl >= 0,  zs[k] >= 0, k = 0, ..., N-1.
//
//    The inequalities sl >= 0 and zl >= 0 are elementwise vector
//    inequalities.  The inequalities ss[k] >= 0, zs[k] >= 0 are matrix
//    inequalities, i.e., the symmetric matrices ss[k] and zs[k] must be
//    positive semidefinite.  mat(Gs[k]*x) is the symmetric matrix X with
//    X[:] = Gs[k]*x.  For a symmetric matrix, zs[k], vec(zs[k]) is the
//    vector zs[k][:].
//
func Sdp(c, Gl, hl, A, b *matrix.FloatMatrix, Ghs *FloatMatrixSet, solopts *SolverOptions, primalstart, dualstart *FloatMatrixSet) (sol *Solution, err error) {
	if c == nil {
		err = errors.New("'c' must a column matrix")
		return
	}
	n := c.Rows()
	if n < 1 {
		err = errors.New("Number of variables must be at least 1")
		return
	}
	if Gl == nil {
		Gl = matrix.FloatZeros(0, n)
	}
	if Gl.Cols() != n {
		err = errors.New(fmt.Sprintf("'G' must be matrix with %d columns", n))
		return
	}
	ml := Gl.Rows()
	if hl == nil {
		hl = matrix.FloatZeros(0, 1)
	}
	if !hl.SizeMatch(ml, 1) {
		err = errors.New(fmt.Sprintf("'hl' must be matrix of size (%d,1)", ml))
		return
	}
	Gsset := Ghs.At("Gs")
	ms := make([]int, 0)
	for i, Gs := range Gsset {
		if Gs.Cols() != n {
			err = errors.New(fmt.Sprintf("'Gs' must be list of matrices with %d columns", n))
			return
		}
		sz := int(math.Sqrt(float64(Gs.Rows())))
		if Gs.Rows() != sz*sz {
			err = errors.New(fmt.Sprintf("the squareroot of the number of rows of 'Gq[%d]' is not an integer", i))
			return
		}
		ms = append(ms, sz)
	}

	hsset := Ghs.At("hs")
	if len(Gsset) != len(hsset) {
		err = errors.New(fmt.Sprintf("'hs' must be a list of %d matrices", len(Gsset)))
		return
	}
	for i, hs := range hsset {
		if !hs.SizeMatch(ms[i], ms[i]) {
			s := fmt.Sprintf("hq[%d] has size (%d,%d). Expected size is (%d,%d)",
				i, hs.Rows(), hs.Cols(), ms[i], ms[i])
			err = errors.New(s)
			return
		}
	}
	if A == nil {
		A = matrix.FloatZeros(0, n)
	}
	if A.Cols() != n {
		err = errors.New(fmt.Sprintf("'A' must be matrix with %d columns", n))
		return
	}
	p := A.Rows()
	if b == nil {
		b = matrix.FloatZeros(0, 1)
	}
	if !b.SizeMatch(p, 1) {
		err = errors.New(fmt.Sprintf("'b' must be matrix of size (%d,1)", p))
		return
	}
	dims := DSetNew("l", "q", "s")
	dims.Set("l", []int{ml})
	dims.Set("s", ms)
	N := dims.Sum("l") + dims.SumSquared("s")

	// Map hs matrices to h vector
	h := matrix.FloatZeros(N, 1)
	h.SetIndexes(matrix.MakeIndexSet(0, ml, 1), hl.FloatArray()[:ml])
	ind := ml
	for k, hs := range hsset {
		h.SetIndexes(matrix.MakeIndexSet(ind, ind+ms[k]*ms[k], 1), hs.FloatArray())
		ind += ms[k] * ms[k]
	}
//.........这里部分代码省略.........

//    Solves a pair of primal and dual SOCPs
//
//        minimize    c'*x
//        subject to  Gl*x + sl = hl
//                    Gq[k]*x + sq[k] = hq[k],  k = 0, ..., N-1
//                    A*x = b
//                    sl >= 0,
//                    sq[k] >= 0, k = 0, ..., N-1
//
//        maximize    -hl'*z - sum_k hq[k]'*zq[k] - b'*y
//        subject to  Gl'*zl + sum_k Gq[k]'*zq[k] + A'*y + c = 0
//                    zl >= 0,  zq[k] >= 0, k = 0, ..., N-1.
//
//    The inequalities sl >= 0 and zl >= 0 are elementwise vector
//    inequalities.  The inequalities sq[k] >= 0, zq[k] >= 0 are second
//    order cone inequalities, i.e., equivalent to
//
//        sq[k][0] >= || sq[k][1:] ||_2,  zq[k][0] >= || zq[k][1:] ||_2.
//
func Socp(c, Gl, hl, A, b *matrix.FloatMatrix, Ghq *FloatMatrixSet, solopts *SolverOptions, primalstart, dualstart *FloatMatrixSet) (sol *Solution, err error) {
	if c == nil {
		err = errors.New("'c' must a column matrix")
		return
	}
	n := c.Rows()
	if n < 1 {
		err = errors.New("Number of variables must be at least 1")
		return
	}
	if Gl == nil {
		Gl = matrix.FloatZeros(0, n)
	}
	if Gl.Cols() != n {
		err = errors.New(fmt.Sprintf("'G' must be matrix with %d columns", n))
		return
	}
	ml := Gl.Rows()
	if hl == nil {
		hl = matrix.FloatZeros(0, 1)
	}
	if !hl.SizeMatch(ml, 1) {
		err = errors.New(fmt.Sprintf("'hl' must be matrix of size (%d,1)", ml))
		return
	}
	Gqset := Ghq.At("Gq")
	mq := make([]int, 0)
	for i, Gq := range Gqset {
		if Gq.Cols() != n {
			err = errors.New(fmt.Sprintf("'Gq' must be list of matrices with %d columns", n))
			return
		}
		if Gq.Rows() == 0 {
			err = errors.New(fmt.Sprintf("the number of rows of 'Gq[%d]' is zero", i))
			return
		}
		mq = append(mq, Gq.Rows())
	}
	hqset := Ghq.At("hq")
	if len(Gqset) != len(hqset) {
		err = errors.New(fmt.Sprintf("'hq' must be a list of %d matrices", len(Gqset)))
		return
	}
	for i, hq := range hqset {
		if !hq.SizeMatch(Gqset[i].Rows(), 1) {
			s := fmt.Sprintf("hq[%d] has size (%d,%d). Expected size is (%d,1)",
				i, hq.Rows(), hq.Cols(), Gqset[i].Rows())
			err = errors.New(s)
			return
		}
	}
	if A == nil {
		A = matrix.FloatZeros(0, n)
	}
	if A.Cols() != n {
		err = errors.New(fmt.Sprintf("'A' must be matrix with %d columns", n))
		return
	}
	p := A.Rows()
	if b == nil {
		b = matrix.FloatZeros(0, 1)
	}
	if !b.SizeMatch(p, 1) {
		err = errors.New(fmt.Sprintf("'b' must be matrix of size (%d,1)", p))
		return
	}
	dims := DSetNew("l", "q", "s")
	dims.Set("l", []int{ml})
	dims.Set("q", mq)
	//N := dims.Sum("l", "q")

	hargs := make([]*matrix.FloatMatrix, 0, len(hqset)+1)
	hargs = append(hargs, hl)
	hargs = append(hargs, hqset...)
	h, indh := matrix.FloatMatrixCombined(matrix.StackDown, hargs...)

	Gargs := make([]*matrix.FloatMatrix, 0, len(Gqset)+1)
	Gargs = append(Gargs, Gl)
	Gargs = append(Gargs, Gqset...)
	G, indg := matrix.FloatMatrixCombined(matrix.StackDown, Gargs...)

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

//    Solves a pair of primal and dual convex quadratic cone programs
//
//        minimize    (1/2)*x'*P*x + q'*x
//        subject to  G*x + s = h
//                    A*x = b
//                    s >= 0
//
//        maximize    -(1/2)*(q + G'*z + A'*y)' * pinv(P) * (q + G'*z + A'*y)
//                    - h'*z - b'*y
//        subject to  q + G'*z + A'*y in range(P)
//                    z >= 0.
//
//    The inequalities are with respect to a cone C defined as the Cartesian
//    product of N + M + 1 cones:
//
//        C = C_0 x C_1 x .... x C_N x C_{N+1} x ... x C_{N+M}.
//
//    The first cone C_0 is the nonnegative orthant of dimension ml.
//    The next N cones are 2nd order cones of dimension mq[0], ..., mq[N-1].
//    The second order cone of dimension m is defined as
//
//        { (u0, u1) in R x R^{m-1} | u0 >= ||u1||_2 }.
//
//    The next M cones are positive semidefinite cones of order ms[0], ...,
//    ms[M-1] >= 0.
//
func ConeQp(P, q, G, h, A, b *matrix.FloatMatrix, dims *DimensionSet, solopts *SolverOptions, initvals *FloatMatrixSet) (sol *Solution, err error) {

	err = nil
	EXPON := 3
	STEP := 0.99

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

	var kktsolver func(*FloatMatrixSet) (kktFunc, error) = nil
	var refinement int
	var correction bool = true

	feasTolerance := FEASTOL
	absTolerance := ABSTOL
	relTolerance := RELTOL
	if solopts.FeasTol > 0.0 {
		feasTolerance = solopts.FeasTol
	}
	if solopts.AbsTol > 0.0 {
		absTolerance = solopts.AbsTol
	}
	if solopts.RelTol > 0.0 {
		relTolerance = solopts.RelTol
	}

	solvername := solopts.KKTSolverName
	if len(solvername) == 0 {
		if dims != nil && (len(dims.At("q")) > 0 || len(dims.At("s")) > 0) {
			solvername = "qr"
			//kktsolver = solvers["qr"]
		} else {
			solvername = "chol2"
			//kktsolver = solvers["chol2"]
		}
	}

	if q == nil || q.Cols() != 1 {
		err = errors.New("'q' must be non-nil matrix with one column")
		return
	}
	if P == nil || P.Rows() != q.Rows() || P.Cols() != q.Rows() {
		err = errors.New(fmt.Sprintf("'P' must be non-nil matrix of size (%d, %d)",
			q.Rows(), q.Rows()))
		return
	}
	fP := func(x, y *matrix.FloatMatrix, alpha, beta float64) error {
		return blas.SymvFloat(P, x, y, alpha, beta)
	}

	if h == nil {
		h = matrix.FloatZeros(0, 1)
	}
	if h.Cols() != 1 {
		err = errors.New("'h' must be non-nil matrix with one column")
		return
	}
	if dims == nil {
		dims = DSetNew("l", "q", "s")
		dims.Set("l", []int{h.Rows()})
	}

	err = checkConeQpDimensions(dims)
	if 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 {
//.........这里部分代码省略.........