Golang Matrix.Cols方法代码示例

/*
 Inverse of a real or complex matrix.

 Getri(A, ipiv, n=A.Rows, ldA = max(1,A.Rows), offsetA=0)

 PURPOSE

 Computes the inverse of real or complex matrix of order n.  On
 entry, A and ipiv contain the LU factorization, as returned by
 gesv() or getrf().  On exit A is replaced by the inverse.

 ARGUMENTS
  A         float or complex matrix
  ipiv      int vector

 OPTIONS
  n         nonnegative integer.  If negative, the default value is used.
  ldA       positive integer.  ldA >= max(1,n).  If zero, the default
            value is used.
  offsetA   nonnegative integer;
*/
func Getri(A matrix.Matrix, ipiv []int32, opts ...linalg.Option) error {
	ind := linalg.GetIndexOpts(opts...)
	if ind.N < 0 {
		ind.N = A.Cols()
	}
	if ind.N == 0 {
		return nil
	}
	if ind.LDa == 0 {
		ind.LDa = max(1, A.Rows())
	}
	if ind.OffsetA < 0 {
		return errors.New("lda")
	}
	sizeA := A.NumElements()
	if sizeA < ind.OffsetA+(ind.N-1)*ind.LDa+ind.N {
		return errors.New("sizeA")
	}
	if ipiv != nil && len(ipiv) < ind.N {
		return errors.New("size ipiv")
	}
	info := -1
	switch A.(type) {
	case *matrix.FloatMatrix:
		Aa := A.FloatArray()
		info = dgetri(ind.N, Aa[ind.OffsetA:], ind.LDa, ipiv)
	case *matrix.ComplexMatrix:
	}
	if info != 0 {
		return errors.New("Getri call error")
	}
	return nil
}

func checkSyevd(ind *linalg.IndexOpts, A, W matrix.Matrix) error {
	if ind.N < 0 {
		ind.N = A.Rows()
		if ind.N != A.Cols() {
			return errors.New("Syevd: A not square")
		}
	}
	if ind.N == 0 {
		return nil
	}
	if ind.LDa == 0 {
		ind.LDa = max(1, A.Rows())
	}
	if ind.LDa < max(1, ind.N) {
		return errors.New("Syevd: lda")
	}
	if ind.OffsetA < 0 {
		return errors.New("Syevd: offsetA")
	}
	sizeA := A.NumElements()
	if sizeA < ind.OffsetA+(ind.N-1)*ind.LDa+ind.N {
		return errors.New("Syevd: sizeA")
	}
	if ind.OffsetW < 0 {
		return errors.New("Syevd: offsetW")
	}
	sizeW := W.NumElements()
	if sizeW < ind.OffsetW+ind.N {
		return errors.New("Syevd: sizeW")
	}
	return nil
}

func checkSytrf(ind *linalg.IndexOpts, A matrix.Matrix, ipiv []int32) error {
	if ind.N < 0 {
		ind.N = A.Rows()
		if ind.N != A.Cols() {
			return errors.New("A not square")
		}
	}
	if ind.N == 0 {
		return nil
	}
	if ind.LDa == 0 {
		ind.LDa = max(1, A.Rows())
	}
	if ind.LDa < max(1, ind.N) {
		return errors.New("Sytrf: lda")
	}
	if ind.OffsetA < 0 {
		return errors.New("Sytrf: offsetA")
	}
	sizeA := A.NumElements()
	if sizeA < ind.OffsetA+(ind.N-1)*ind.LDa+ind.N {
		return errors.New("Sytrf: sizeA")
	}
	if ipiv != nil && len(ipiv) < ind.N {
		return errors.New("Sytrf: size ipiv")
	}
	return nil
}

/*
 Solves a real symmetric or complex Hermitian positive definite set
 of linear equations, given the Cholesky factorization computed by
 potrf() or posv().

 Potrs(A, B, uplo=PLower, n=A.Rows, nrhs=B.Cols,
 ldA=max(1,A.Rows), ldB=max(1,B.Rows), offsetA=0, offsetB=0)

 PURPOSE

 Solves
   A*X = B

 where A is n by n, real symmetric or complex Hermitian and positive definite,
 and B is n by nrhs. On entry, A contains the Cholesky factor, as
 returned by Posv() or Potrf().  On exit B is replaced by the solution X.

 ARGUMENTS
  A         float or complex matrix
  B         float or complex matrix.  Must have the same type as A.

 OPTIONS
  uplo      PLower or PUpper
  n         nonnegative integer.  If negative, the default value is used.
  nrhs      nonnegative integer.  If negative, the default value is used.
  ldA       positive integer.  ldA >= max(1,n).  If zero, the default
            value is used.
  ldB       positive integer.  ldB >= max(1,n).  If zero, the default
            value is used.
  offsetA   nonnegative integer
  offsetB   nonnegative integer;

*/
func Potrs(A, B matrix.Matrix, opts ...linalg.Option) error {
	pars, err := linalg.GetParameters(opts...)
	if err != nil {
		return err
	}
	ind := linalg.GetIndexOpts(opts...)
	if ind.N < 0 {
		ind.N = A.Rows()
	}
	if ind.Nrhs < 0 {
		ind.Nrhs = B.Cols()
	}
	if ind.N == 0 || ind.Nrhs == 0 {
		return nil
	}
	if ind.LDa == 0 {
		ind.LDa = max(1, A.Rows())
	}
	if ind.LDa < max(1, ind.N) {
		return errors.New("lda")
	}
	if ind.LDb == 0 {
		ind.LDb = max(1, B.Rows())
	}
	if ind.LDb < max(1, ind.N) {
		return errors.New("ldb")
	}
	if ind.OffsetA < 0 {
		return errors.New("offsetA")
	}
	if A.NumElements() < ind.OffsetA+(ind.N-1)*ind.LDa+ind.N {
		return errors.New("sizeA")
	}
	if ind.OffsetB < 0 {
		return errors.New("offsetB")
	}
	if B.NumElements() < ind.OffsetB+(ind.Nrhs-1)*ind.LDb+ind.N {
		return errors.New("sizeB")
	}
	if !matrix.EqualTypes(A, B) {
		return errors.New("types")
	}
	info := -1
	switch A.(type) {
	case *matrix.FloatMatrix:
		Aa := A.FloatArray()
		Ba := B.FloatArray()
		uplo := linalg.ParamString(pars.Uplo)
		info = dpotrs(uplo, ind.N, ind.Nrhs, Aa[ind.OffsetA:], ind.LDa,
			Ba[ind.OffsetB:], ind.LDb)
	case *matrix.ComplexMatrix:
		return errors.New("ComplexMatrx: not implemented yet")
	}
	if info != 0 {
		return errors.New("Potrs failed")
	}
	return nil
}

func checkGbtrs(ind *linalg.IndexOpts, A, B matrix.Matrix, ipiv []int32) error {
	if ind.Kl < 0 {
		return errors.New("Gbtrs: invalid kl")
	}
	if ind.N < 0 {
		ind.N = A.Rows()
	}
	if ind.Nrhs < 0 {
		ind.Nrhs = A.Cols()
	}
	if ind.N == 0 || ind.Nrhs == 0 {
		return nil
	}
	if ind.Ku < 0 {
		ind.Ku = A.Rows() - 2*ind.Kl - 1
	}
	if ind.Ku < 0 {
		return errors.New("Gbtrs: invalid ku")
	}
	if ind.LDa == 0 {
		ind.LDa = max(1, A.Rows())
	}
	if ind.LDa < 2*ind.Kl+ind.Ku+1 {
		return errors.New("Gbtrs: lda")
	}
	if ind.OffsetA < 0 {
		return errors.New("Gbtrs: offsetA")
	}
	sizeA := A.NumElements()
	if sizeA < ind.OffsetA+(ind.N-1)*ind.LDa+2*ind.Kl+ind.Ku+1 {
		return errors.New("Gbtrs: sizeA")
	}
	if ind.LDb == 0 {
		ind.LDb = max(1, B.Rows())
	}
	if ind.OffsetB < 0 {
		return errors.New("Gbtrs: offsetB")
	}
	sizeB := B.NumElements()
	if sizeB < ind.OffsetB+(ind.Nrhs-1)*ind.LDb+ind.N {
		return errors.New("Gbtrs: sizeB")
	}
	if ipiv != nil && len(ipiv) < ind.N {
		return errors.New("Gbtrs: size ipiv")
	}
	return nil
}

/*
 QR factorization.

 Geqrf(A, tau, m=A.Rows, n=A.Cols, ldA=max(1,A.Rows), offsetA=0)

 PURPOSE

 QR factorization of an m by n real or complex matrix A:

 A = Q*R = [Q1 Q2] * [R1; 0] if m >= n
 A = Q*R = Q * [R1 R2]       if m <= n,

 where Q is m by m and orthogonal/unitary and R is m by n with R1
 upper triangular.  On exit, R is stored in the upper triangular
 part of A.  Q is stored as a product of k=min(m,n) elementary
 reflectors.  The parameters of the reflectors are stored in the
 first k entries of tau and in the lower triangular part of the
 first k columns of A.

 ARGUMENTS
  A         float or complex matrix
  tau       float or complex  matrix of length at least min(m,n).  Must
            have the same type as A.
  m         integer.  If negative, the default value is used.
  n         integer.  If negative, the default value is used.
  ldA       nonnegative integer.  ldA >= max(1,m).  If zero, the
            default value is used.
  offsetA   nonnegative integer

*/
func Geqrf(A, tau matrix.Matrix, opts ...linalg.Option) error {
	ind := linalg.GetIndexOpts(opts...)
	if ind.N < 0 {
		ind.N = A.Rows()
	}
	if ind.M < 0 {
		ind.M = A.Cols()
	}
	if ind.N == 0 || ind.M == 0 {
		return nil
	}
	if ind.LDa == 0 {
		ind.LDa = max(1, A.Rows())
	}
	if ind.LDa < max(1, ind.M) {
		return errors.New("lda")
	}
	if ind.OffsetA < 0 {
		return errors.New("offsetA")
	}
	if A.NumElements() < ind.OffsetA+ind.K*ind.LDa {
		return errors.New("sizeA")
	}
	if tau.NumElements() < min(ind.M, ind.N) {
		return errors.New("sizeTau")
	}
	if !matrix.EqualTypes(A, tau) {
		return errors.New("not same type")
	}
	info := -1
	switch A.(type) {
	case *matrix.FloatMatrix:
		Aa := A.FloatArray()
		taua := tau.FloatArray()
		info = dgeqrf(ind.M, ind.N, Aa[ind.OffsetA:], ind.LDa, taua)
	case *matrix.ComplexMatrix:
		return errors.New("ComplexMatrx: not implemented yet")
	}
	if info != 0 {
		return errors.New("Geqrf failed")
	}
	return nil
}

func checkPotrf(ind *linalg.IndexOpts, A matrix.Matrix) error {
	if ind.N < 0 {
		ind.N = A.Rows()
		if ind.N != A.Cols() {
			return errors.New("Potrf: not square")
		}
	}
	if ind.N == 0 {
		return nil
	}
	if ind.LDa == 0 {
		ind.LDa = max(1, A.Rows())
	}
	if ind.LDa < max(1, ind.N) {
		return errors.New("Potrf: lda")
	}
	if ind.OffsetA < 0 {
		return errors.New("Potrf: offsetA")
	}
	if A.NumElements() < ind.OffsetA+(ind.N-1)*ind.LDa+ind.N {
		return errors.New("Potrf: sizeA")
	}
	return nil
}

func SyevrFloat(A, W, Z matrix.Matrix, abstol float64, vlimit []float64, ilimit []int, opts ...linalg.Option) error {
	var vl, vu float64
	var il, iu int

	pars, err := linalg.GetParameters(opts...)
	if err != nil {
		return err
	}
	ind := linalg.GetIndexOpts(opts...)
	if ind.N < 0 {
		ind.N = A.Rows()
		if ind.N != A.Cols() {
			return errors.New("Syevr: A not square")
		}
	}
	// Check indexes
	if ind.N == 0 {
		return nil
	}
	if ind.LDa == 0 {
		ind.LDa = max(1, A.Rows())
	}
	if ind.LDa < max(1, A.Rows()) {
		return errors.New("Syevr: lda")
	}
	if pars.Range == linalg.PRangeValue {
		if vlimit == nil {
			return errors.New("Syevr: vlimit is nil")
		}
		vl = vlimit[0]
		vu = vlimit[1]
		if vl >= vu {
			return errors.New("Syevr: must be: vl < vu")
		}
	} else if pars.Range == linalg.PRangeInt {
		if ilimit == nil {
			return errors.New("Syevr: ilimit is nil")
		}
		il = ilimit[0]
		iu = ilimit[1]
		if il < 1 || il > iu || iu > ind.N {
			return errors.New("Syevr: must be:1 <= il <= iu <= N")
		}
	}
	if pars.Jobz == linalg.PJobValue {
		if Z == nil {
			return errors.New("Syevr: Z is nil")
		}
		if ind.LDz == 0 {
			ind.LDz = max(1, Z.Rows())
		}
		if ind.LDz < max(1, ind.N) {
			return errors.New("Syevr: ldz")
		}
	} else {
		if ind.LDz == 0 {
			ind.LDz = 1
		}
		if ind.LDz < 1 {
			return errors.New("Syevr: ldz")
		}
	}
	if ind.OffsetA < 0 {
		return errors.New("Syevr: OffsetA")
	}
	sizeA := A.NumElements()
	if sizeA < ind.OffsetA+(ind.N-1)*ind.LDa+ind.N {
		return errors.New("Syevr: sizeA")
	}
	if ind.OffsetW < 0 {
		return errors.New("Syevr: OffsetW")
	}
	sizeW := W.NumElements()
	if sizeW < ind.OffsetW+ind.N {
		return errors.New("Syevr: sizeW")
	}
	if pars.Jobz == linalg.PJobValue {
		if ind.OffsetZ < 0 {
			return errors.New("Syevr: OffsetW")
		}
		minZ := ind.OffsetZ + (ind.N-1)*ind.LDz + ind.N
		if pars.Range == linalg.PRangeInt {
			minZ = ind.OffsetZ + (iu-il)*ind.LDz + ind.N
		}
		if Z.NumElements() < minZ {
			return errors.New("Syevr: sizeZ")
		}
	}

	Aa := A.FloatArray()
	Wa := W.FloatArray()
	var Za []float64
	if pars.Jobz == linalg.PJobValue {
		Za = Z.FloatArray()
	} else {
		Za = nil
	}
	jobz := linalg.ParamString(pars.Jobz)
	rnge := linalg.ParamString(pars.Range)
	uplo := linalg.ParamString(pars.Uplo)
//.........这里部分代码省略.........

/*
 Solves a real or complex set of linear equations with a banded
 coefficient matrix, given the LU factorization computed by gbtrf()
 or gbsv().

 Gbtrs(A, B, ipiv, kl, trans=PNoTrans, n=A.Cols, ku=A.Rows-2*kl-1,
 nrhs=B.Cols, ldA=max(1,A.Rows), ldB=max(1,B.Rows), offsetA=0, offsetB=0)

 PURPOSE

 Solves linear equations
  A*X = B,   if trans is PNoTrans
  A^T*X = B, if trans is PTrans
  A^H*X = B, if trans is PConjTrans

 On entry, A and ipiv contain the LU factorization of an n by n
 band matrix A as computed by Getrf() or Gbsv().  On exit B is
 replaced by the solution X.

 ARGUMENTS
  A         float or complex matrix
  B         float or complex  matrix.  Must have the same type as A.
  ipiv      int vector
  kl        nonnegative integer

 OPTIONS
  trans     PNoTrans, PTrans or PConjTrans
  n         nonnegative integer.  If negative, the default value is used.
  ku        nonnegative integer.  If negative, the default value is used.
  nrhs      nonnegative integer.  If negative, the default value is used.
  ldA       positive integer, ldA >= 2*kl+ku+1. If zero, the  default value is used.
  ldB       positive integer, ldB >= max(1,n). If zero, the default value is used.
  offsetA   nonnegative integer
  offsetB   nonnegative integer;
*/
func Gbtrs(A, B matrix.Matrix, ipiv []int32, KL int, opts ...linalg.Option) error {
	pars, err := linalg.GetParameters(opts...)
	if err != nil {
		return err
	}
	ind := linalg.GetIndexOpts(opts...)
	//err = lapack_check(ind, fgbtrs, A, B, ipiv, pars)
	ind.Kl = KL
	if ind.Kl < 0 {
		return errors.New("invalid kl")
	}
	if ind.N < 0 {
		ind.N = A.Rows()
	}
	if ind.Nrhs < 0 {
		ind.Nrhs = A.Cols()
	}
	if ind.N == 0 || ind.Nrhs == 0 {
		return nil
	}
	if ind.Ku < 0 {
		ind.Ku = A.Rows() - 2*ind.Kl - 1
	}
	if ind.Ku < 0 {
		return errors.New("invalid ku")
	}
	if ind.LDa == 0 {
		ind.LDa = max(1, A.Rows())
	}
	if ind.LDa < 2*ind.Kl+ind.Ku+1 {
		return errors.New("lda")
	}
	if ind.OffsetA < 0 {
		return errors.New("offsetA")
	}
	sizeA := A.NumElements()
	if sizeA < ind.OffsetA+(ind.N-1)*ind.LDa+2*ind.Kl+ind.Ku+1 {
		return errors.New("sizeA")
	}
	if ind.LDb == 0 {
		ind.LDb = max(1, B.Rows())
	}
	if ind.OffsetB < 0 {
		return errors.New("offsetB")
	}
	sizeB := B.NumElements()
	if sizeB < ind.OffsetB+(ind.Nrhs-1)*ind.LDb+ind.N {
		return errors.New("sizeB")
	}
	if ipiv != nil && len(ipiv) < ind.N {
		return errors.New("size ipiv")
	}

	info := -1
	switch A.(type) {
	case *matrix.FloatMatrix:
		Aa := A.FloatArray()
		Ba := B.FloatArray()
		trans := linalg.ParamString(pars.Trans)
		info = dgbtrs(trans, ind.N, ind.Kl, ind.Ku, ind.Nrhs,
			Aa[ind.OffsetA:], ind.LDa, ipiv, Ba[ind.OffsetB:], ind.LDb)
	case *matrix.ComplexMatrix:
		return errors.New("Gbtrs for complex not yet implemented")
	}
	if info != 0 {
//.........这里部分代码省略.........

func Gtrrs(DL, D, DU, DU2, B matrix.Matrix, ipiv []int32, opts ...linalg.Option) error {
	pars, err := linalg.GetParameters(opts...)
	if err != nil {
		return err
	}
	ind := linalg.GetIndexOpts(opts...)
	if ind.OffsetD < 0 {
		return errors.New("offset D")
	}
	if ind.N < 0 {
		ind.N = D.NumElements() - ind.OffsetD
	}
	if ind.N < 0 {
		return errors.New("size D")
	}
	if ind.N == 0 {
		return nil
	}
	if ind.OffsetDL < 0 {
		return errors.New("offset DL")
	}
	sizeDL := DL.NumElements()
	if sizeDL < ind.OffsetDL+ind.N-1 {
		return errors.New("sizeDL")
	}
	if ind.OffsetDU < 0 {
		return errors.New("offset DU")
	}
	sizeDU := DU.NumElements()
	if sizeDU < ind.OffsetDU+ind.N-1 {
		return errors.New("sizeDU")
	}
	sizeDU2 := DU2.NumElements()
	if sizeDU2 < ind.N-2 {
		return errors.New("sizeDU2")
	}
	if ind.Nrhs < 0 {
		ind.Nrhs = B.Cols()
	}
	if ind.Nrhs == 0 {
		return nil
	}
	if ind.LDb == 0 {
		ind.LDb = max(1, B.Rows())
	}
	if ind.LDb < max(1, ind.N) {
		return errors.New("ldB")
	}
	if ind.OffsetB < 0 {
		return errors.New("offset B")
	}
	sizeB := B.NumElements()
	if sizeB < ind.OffsetB+(ind.Nrhs-1)*ind.LDb+ind.N {
		return errors.New("sizeB")
	}
	if len(ipiv) < ind.N {
		return errors.New("size ipiv")
	}
	DLa := DL.FloatArray()
	Da := D.FloatArray()
	DUa := DU.FloatArray()
	DU2a := DU2.FloatArray()
	Ba := B.FloatArray()
	trans := linalg.ParamString(pars.Trans)
	info := dgttrs(trans, ind.N, ind.Nrhs,
		DLa[ind.OffsetDL:], Da[ind.OffsetD:], DUa[ind.OffsetDU:], DU2a,
		ipiv, Ba[ind.OffsetB:], ind.LDb)
	if info != 0 {
		return errors.New("dgttrs call error")
	}
	return nil
}

/*
 Solution of a triangular set of equations with multiple righthand
 sides.

 Trtrs(A, B, uplo=PLower, trans=PNoTrans, diag=PNonUnit, n=A.Rows,
 nrhs=B.Cols, ldA=max(1,A.Rows), ldB=max(1,B.Rows), offsetA=0, offsetB=0)

 PURPOSE
 Solves set of equations
  A*X = B,   if trans is PNoTrans
  A^T*X = B, if trans is PTrans
  A^H*X = B, if trans is PConjTrans

 B is n by nrhs and A is triangular of order n.

 ARGUMENTS
  A         float or complex matrix
  B         float or complex matrix.  Must have the same type as A.

 OPTIONS
  uplo      PLower or PUpper
  trans     PNoTrans, PTrans, PConjTrans
  diag      PNonUnit, PUnit
  n         nonnegative integer.  If negative, the default value is used.
  nrhs      nonnegative integer.  If negative, the default value is used.
  ldA       positive integer.  ldA >= max(1,n).  If zero, the default
            value is used.
  ldB       positive integer.  ldB >= max(1,n).  If zero, the default
            value is used.
  offsetA   nonnegative integer
  offsetB   nonnegative integer;

*/
func Trtrs(A, B matrix.Matrix, ipiv []int32, opts ...linalg.Option) error {
	pars, err := linalg.GetParameters(opts...)
	if err != nil {
		return err
	}
	ind := linalg.GetIndexOpts(opts...)
	if ind.N < 0 {
		ind.N = A.Rows()
		if ind.N != A.Cols() {
			return errors.New("A not square")
		}
	}
	if ind.Nrhs < 0 {
		ind.Nrhs = B.Cols()
	}
	if ind.N == 0 || ind.Nrhs == 0 {
		return nil
	}
	if ind.LDa == 0 {
		ind.LDa = max(1, A.Rows())
	}
	if ind.LDa < max(1, ind.N) {
		return errors.New("lda")
	}
	if ind.LDb == 0 {
		ind.LDb = max(1, B.Rows())
	}
	if ind.LDb < max(1, ind.N) {
		return errors.New("ldb")
	}
	if ind.OffsetA < 0 {
		return errors.New("offsetA")
	}
	sizeA := A.NumElements()
	if sizeA < ind.OffsetA+(ind.N-1)*ind.LDa+ind.N {
		return errors.New("sizeA")
	}
	if ind.OffsetB < 0 {
		return errors.New("offsetB")
	}
	sizeB := B.NumElements()
	if sizeB < ind.OffsetB+(ind.Nrhs-1)*ind.LDb+ind.N {
		return errors.New("sizeB")
	}
	info := -1
	uplo := linalg.ParamString(pars.Uplo)
	trans := linalg.ParamString(pars.Trans)
	diag := linalg.ParamString(pars.Diag)
	switch A.(type) {
	case *matrix.FloatMatrix:
		Aa := A.FloatArray()
		Ba := B.FloatArray()
		info = dtrtrs(uplo, trans, diag, ind.N, ind.Nrhs, Aa[ind.OffsetA:], ind.LDa,
			Ba[ind.OffsetB:], ind.LDb)
	case *matrix.ComplexMatrix:
	}
	if info != 0 {
		return errors.New("sytrf failed")
	}
	return nil
}

func check_level3_func(ind *linalg.IndexOpts, fn funcNum, A, B, C matrix.Matrix,
	pars *linalg.Parameters) (err error) {

	switch fn {
	case fgemm:
		if ind.M < 0 {
			if pars.TransA == linalg.PNoTrans {
				ind.M = A.Rows()
			} else {
				ind.M = A.Cols()
			}
		}
		if ind.N < 0 {
			if pars.TransB == linalg.PNoTrans {
				ind.N = A.Cols()
			} else {
				ind.N = A.Rows()
			}
		}
		if ind.K < 0 {
			if pars.TransA == linalg.PNoTrans {
				ind.K = A.Cols()
			} else {
				ind.K = A.Rows()
			}
			if pars.TransB == linalg.PNoTrans && ind.K != B.Rows() ||
				pars.TransB != linalg.PNoTrans && ind.K != B.Cols() {
				return errors.New("dimensions of A and B do not match")
			}
		}
		if ind.OffsetB < 0 {
			return errors.New("offsetB illegal, <0")
		}
		if ind.LDa == 0 {
			ind.LDa = max(1, A.Rows())
		}
		if ind.K > 0 {
			if (pars.TransA == linalg.PNoTrans && ind.LDa < max(1, ind.M)) ||
				(pars.TransA != linalg.PNoTrans && ind.LDa < max(1, ind.K)) {
				return errors.New("inconsistent ldA")
			}
			sizeA := A.NumElements()
			if (pars.TransA == linalg.PNoTrans &&
				sizeA < ind.OffsetA+(ind.K-1)*ind.LDa+ind.M) ||
				(pars.TransA != linalg.PNoTrans &&
					sizeA < ind.OffsetA+(ind.M-1)*ind.LDa+ind.K) {
				return errors.New("sizeA")
			}
		}
		// B matrix
		if B != nil {
			if ind.OffsetB < 0 {
				return errors.New("offsetB illegal, <0")
			}
			if ind.LDb == 0 {
				ind.LDb = max(1, B.Rows())
			}
			if ind.K > 0 {
				if (pars.TransB == linalg.PNoTrans && ind.LDb < max(1, ind.K)) ||
					(pars.TransB != linalg.PNoTrans && ind.LDb < max(1, ind.N)) {
					return errors.New("inconsistent ldB")
				}
				sizeB := B.NumElements()
				if (pars.TransB == linalg.PNoTrans &&
					sizeB < ind.OffsetB+(ind.N-1)*ind.LDb+ind.K) ||
					(pars.TransB != linalg.PNoTrans &&
						sizeB < ind.OffsetB+(ind.K-1)*ind.LDb+ind.N) {
					return errors.New("sizeB")
				}
			}
		}
		// C matrix
		if C != nil {
			if ind.OffsetC < 0 {
				return errors.New("offsetC illegal, <0")
			}
			if ind.LDc == 0 {
				ind.LDb = max(1, C.Rows())
			}
			if ind.LDc < max(1, ind.M) {
				return errors.New("inconsistent ldC")
			}
			sizeC := C.NumElements()
			if sizeC < ind.OffsetC+(ind.N-1)*ind.LDc+ind.M {
				return errors.New("sizeC")
			}
		}

	case fsymm, ftrmm, ftrsm:
		if ind.M < 0 {
			ind.M = B.Rows()
			if pars.Side == linalg.PLeft && (ind.M != A.Rows() || ind.M != A.Cols()) {
				return errors.New("dimensions of A and B do not match")
			}
		}
		if ind.N < 0 {
			ind.N = B.Cols()
			if pars.Side == linalg.PRight && (ind.N != A.Rows() || ind.N != A.Cols()) {
				return errors.New("dimensions of A and B do not match")
			}
//.........这里部分代码省略.........

func check_level2_func(ind *linalg.IndexOpts, fn funcNum, X, Y, A matrix.Matrix, pars *linalg.Parameters) error {
	if ind.IncX <= 0 {
		return errors.New("incX")
	}
	if ind.IncY <= 0 {
		return errors.New("incY")
	}

	sizeA := A.NumElements()
	switch fn {
	case fgemv: // general matrix
		if ind.M < 0 {
			ind.M = A.Rows()
		}
		if ind.N < 0 {
			ind.N = A.Cols()
		}
		if ind.LDa == 0 {
			ind.LDa = max(1, A.Rows())
		}
		if ind.OffsetA < 0 {
			return errors.New("offsetA")
		}
		if ind.N > 0 && ind.M > 0 &&
			sizeA < ind.OffsetA+(ind.N-1)*ind.LDa+ind.M {
			return errors.New("sizeA")
		}
		if ind.OffsetX < 0 {
			return errors.New("offsetX")
		}
		if ind.OffsetY < 0 {
			return errors.New("offsetY")
		}
		sizeX := X.NumElements()
		sizeY := Y.NumElements()
		if pars.Trans == linalg.PNoTrans {
			if ind.N > 0 && sizeX < ind.OffsetX+(ind.N-1)*abs(ind.IncX)+1 {
				return errors.New("sizeX")
			}
			if ind.M > 0 && sizeY < ind.OffsetY+(ind.M-1)*abs(ind.IncY)+1 {
				return errors.New("sizeY")
			}
		} else {
			if ind.M > 0 && sizeX < ind.OffsetX+(ind.M-1)*abs(ind.IncX)+1 {
				return errors.New("sizeX")
			}
			if ind.N > 0 && sizeY < ind.OffsetY+(ind.N-1)*abs(ind.IncY)+1 {
				return errors.New("sizeY")
			}
		}
	case fger:
		if ind.M < 0 {
			ind.M = A.Rows()
		}
		if ind.N < 0 {
			ind.N = A.Cols()
		}
		if ind.M == 0 || ind.N == 0 {
			return nil
		}
		if ind.M > 0 && ind.N > 0 {
			if ind.LDa == 0 {
				ind.LDa = max(1, A.Rows())
			}
			if ind.LDa < max(1, ind.M) {
				return errors.New("ldA")
			}
			if ind.OffsetA < 0 {
				return errors.New("offsetA")
			}
			if sizeA < ind.OffsetA+(ind.N-1)*ind.LDa+ind.M {
				return errors.New("sizeA")
			}
			if ind.OffsetX < 0 {
				return errors.New("offsetX")
			}
			if ind.OffsetY < 0 {
				return errors.New("offsetY")
			}
			sizeX := X.NumElements()
			if sizeX < ind.OffsetX+(ind.M-1)*abs(ind.IncX)+1 {
				return errors.New("sizeX")
			}
			sizeY := Y.NumElements()
			if sizeY < ind.OffsetY+(ind.N-1)*abs(ind.IncY)+1 {
				return errors.New("sizeY")
			}
		}
	case fgbmv: // general banded
		if ind.M < 0 {
			ind.M = A.Rows()
		}
		if ind.N < 0 {
			ind.N = A.Cols()
		}
		if ind.Kl < 0 {
			return errors.New("kl")
		}
		if ind.Ku < 0 {
			ind.Ku = A.Rows() - 1 - ind.Kl
//.........这里部分代码省略.........

func checkGesvd(ind *linalg.IndexOpts, pars *linalg.Parameters, A, S, U, Vt matrix.Matrix) error {
	if ind.M < 0 {
		ind.M = A.Rows()
	}
	if ind.N < 0 {
		ind.N = A.Cols()
	}
	if ind.M == 0 || ind.N == 0 {
		return nil
	}
	if pars.Jobu == linalg.PJobO && pars.Jobvt == linalg.PJobO {
		return errors.New("Gesvd: jobu and jobvt cannot both have value PJobO")
	}
	if pars.Jobu == linalg.PJobAll || pars.Jobu == linalg.PJobS {
		if U == nil {
			return errors.New("Gesvd: missing matrix U")
		}
		if ind.LDu == 0 {
			ind.LDu = max(1, U.Rows())
		}
		if ind.LDu < max(1, ind.M) {
			return errors.New("Gesvd: ldU")
		}
	} else {
		if ind.LDu == 0 {
			ind.LDu = 1
		}
		if ind.LDu < 1 {
			return errors.New("Gesvd: ldU")
		}
	}
	if pars.Jobvt == linalg.PJobAll || pars.Jobvt == linalg.PJobS {
		if Vt == nil {
			return errors.New("Gesvd: missing matrix Vt")
		}
		if ind.LDvt == 0 {
			ind.LDvt = max(1, Vt.Rows())
		}
		if pars.Jobvt == linalg.PJobAll && ind.LDvt < max(1, ind.N) {
			return errors.New("Gesvd: ldVt")
		} else if pars.Jobvt != linalg.PJobAll && ind.LDvt < max(1, min(ind.M, ind.N)) {
			return errors.New("Gesvd: ldVt")
		}
	} else {
		if ind.LDvt == 0 {
			ind.LDvt = 1
		}
		if ind.LDvt < 1 {
			return errors.New("Gesvd: ldVt")
		}
	}
	if ind.OffsetA < 0 {
		return errors.New("Gesvd: offsetA")
	}
	sizeA := A.NumElements()
	if sizeA < ind.OffsetA+(ind.N-1)*ind.LDa+ind.M {
		return errors.New("Gesvd: sizeA")
	}

	if ind.OffsetS < 0 {
		return errors.New("Gesvd: offsetS")
	}
	sizeS := S.NumElements()
	if sizeS < ind.OffsetS+min(ind.M, ind.N) {
		return errors.New("Gesvd: sizeA")
	}

	/*
		if U != nil {
			if ind.OffsetU < 0 {
				return errors.New("Gesvd: OffsetU")
			}
			sizeU := U.NumElements()
			if pars.Jobu == linalg.PJobAll && sizeU < ind.LDu*(ind.M-1) {
				return errors.New("Gesvd: sizeU")
			} else if pars.Jobu == linalg.PJobS && sizeU < ind.LDu*(min(ind.M,ind.N)-1) {
				return errors.New("Gesvd: sizeU")
			}
		}

		if Vt != nil {
			if ind.OffsetVt < 0 {
				return errors.New("Gesvd: OffsetVt")
			}
			sizeVt := Vt.NumElements()
			if pars.Jobvt == linalg.PJobAll && sizeVt <  ind.N {
				return errors.New("Gesvd: sizeVt")
			} else if pars.Jobvt == linalg.PJobS && sizeVt < min(ind.M, ind.N) {
				return errors.New("Gesvd: sizeVt")
			}
		}
	*/
	return nil
}

/*
 Product with a real orthogonal matrix.

 Ormqr(A, tau, C, side='L', trans='N', m=C.Rows, n=C.Cols,
 k=len(tau), ldA=max(1,A.Rows), ldC=max(1,C.Rows), offsetA=0, offsetC=0)

 PURPOSE

 Computes
  C := Q*C   if side = PLeft  and trans = PNoTrans
  C := Q^T*C if side = PLeft  and trans = PTrans
  C := C*Q   if side = PRight and trans = PNoTrans
  C := C*Q^T if side = PRight and trans = PTrans

 C is m by n and Q is a square orthogonal matrix computed by geqrf.

 Q is defined as a product of k elementary reflectors, stored as
 the first k columns of A and the first k entries of tau.

 ARGUMENTS
  A         float matrix
  tau       float matrix of length at least k
  C         float matrix

 OPTIONS
  side      PLeft or PRight
  trans     PNoTrans or PTrans
  m         integer.  If negative, the default value is used.
  n         integer.  If negative, the default value is used.
  k         integer.  k <= m if side = PRight and k <= n if side = PLeft.
            If negative, the default value is used.
  ldA       nonnegative integer.  ldA >= max(1,m) if side = PLeft
            and ldA >= max(1,n) if side = PRight.  If zero, the
            default value is used.
  ldC       nonnegative integer.  ldC >= max(1,m).  If zero, the
            default value is used.
  offsetA   nonnegative integer
  offsetB   nonnegative integer

*/
func Ormqf(A, tau, C matrix.Matrix, opts ...linalg.Option) error {
	pars, err := linalg.GetParameters(opts...)
	if err != nil {
		return err
	}
	ind := linalg.GetIndexOpts(opts...)
	if ind.N < 0 {
		ind.N = C.Cols()
	}
	if ind.M < 0 {
		ind.M = C.Rows()
	}
	if ind.K < 0 {
		ind.K = tau.NumElements()
	}
	if ind.N == 0 || ind.M == 0 || ind.K == 0 {
		return nil
	}
	if ind.LDa == 0 {
		ind.LDa = max(1, A.Rows())
	}
	if ind.LDc == 0 {
		ind.LDc = max(1, C.Rows())
	}
	switch pars.Side {
	case linalg.PLeft:
		if ind.K > ind.M {
			errors.New("K")
		}
		if ind.LDa < max(1, ind.M) {
			return errors.New("lda")
		}
	case linalg.PRight:
		if ind.K > ind.N {
			errors.New("K")
		}
		if ind.LDa < max(1, ind.N) {
			return errors.New("lda")
		}
	}
	if ind.OffsetA < 0 {
		return errors.New("offsetA")
	}
	if A.NumElements() < ind.OffsetA+ind.K*ind.LDa {
		return errors.New("sizeA")
	}
	if ind.OffsetC < 0 {
		return errors.New("offsetC")
	}
	if C.NumElements() < ind.OffsetC+(ind.N-1)*ind.LDa+ind.M {
		return errors.New("sizeC")
	}
	if tau.NumElements() < ind.K {
		return errors.New("sizeTau")
	}
	if !matrix.EqualTypes(A, C, tau) {
		return errors.New("not same type")
	}
	info := -1
	side := linalg.ParamString(pars.Side)
//.........这里部分代码省略.........