Golang Matrix.NumElements方法代碼示例

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

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
}

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

/*
 LU factorization of a real or complex tridiagonal matrix.

 Gttrf(dl, d, du, du2, ipiv, n=len(d)-offsetd, offsetdl=0, offsetd=0, offsetdu=0)

 PURPOSE

 Factors an n by n real or complex tridiagonal matrix A as A = P*L*U.

 A is specified by its lower diagonal dl, diagonal d, and upper
 diagonal du.  On exit dl, d, du, du2 and ipiv contain the details
 of the factorization.

 ARGUMENTS.
  DL        float or complex matrix
  D         float or complex matrix.  Must have the same type as DL.
  DU        float or complex matrix.  Must have the same type as DL.
  DU2       float or complex matrix of length at least n-2.  Must have the
            same type as DL.
  ipiv      int vector of length at least n

 OPTIONS
  n         nonnegative integer.  If negative, the default value is used.
  offsetdl  nonnegative integer
  offsetd   nonnegative integer
  offsetdu  nonnegative integer
*/
func Gtrrf(DL, D, DU, DU2 matrix.Matrix, ipiv []int32, opts ...linalg.Option) error {
	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 len(ipiv) < ind.N {
		return errors.New("size ipiv")
	}
	info := -1
	switch DL.(type) {
	case *matrix.FloatMatrix:
		DLa := DL.FloatArray()
		Da := D.FloatArray()
		DUa := DU.FloatArray()
		DU2a := DU2.FloatArray()
		info = dgttrf(ind.N, DLa[ind.OffsetDL:], Da[ind.OffsetD:], DUa[ind.OffsetDU:],
			DU2a, ipiv)
	case *matrix.ComplexMatrix:
	}
	if info != 0 {
		return errors.New("Gttrf call error")
	}
	return nil
}

func checkPotri(ind *linalg.IndexOpts, A matrix.Matrix) error {
	if ind.N < 0 {
		ind.N = A.Rows()
	}
	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("Potri: lda")
	}
	if ind.OffsetA < 0 {
		return errors.New("Potri: offsetA")
	}
	if A.NumElements() < ind.OffsetA+(ind.N-1)*ind.LDa+ind.N {
		return errors.New("Potri: sizeA")
	}
	return nil
}

func checkGbtrf(ind *linalg.IndexOpts, A matrix.Matrix, ipiv []int32) error {
	if ind.M < 0 {
		return errors.New("Gbtrf: illegal m")
	}
	if ind.Kl < 0 {
		return errors.New("GBtrf: illegal kl")
	}
	if ind.N < 0 {
		ind.N = A.Rows()
	}
	if ind.M == 0 || ind.N == 0 {
		return nil
	}
	if ind.Ku < 0 {
		ind.Ku = A.Rows() - 2*ind.Kl - 1
	}
	if ind.Ku < 0 {
		return errors.New("Gbtrf: invalid ku")
	}
	if ind.LDa == 0 {
		ind.LDa = max(1, A.Rows())
	}
	if ind.LDa < 2*ind.Kl+ind.Ku+1 {
		return errors.New("Gbtrf: lda")
	}
	if ind.OffsetA < 0 {
		return errors.New("Gbtrf: offsetA")
	}
	sizeA := A.NumElements()
	if sizeA < ind.OffsetA+(ind.N-1)*ind.LDa+2*ind.Kl+ind.Ku+1 {
		return errors.New("Gbtrf: sizeA")
	}
	if ipiv != nil && len(ipiv) < min(ind.N, ind.M) {
		return errors.New("Gbtrf: size ipiv")
	}
	return nil
}

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

/*
 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)
//.........這裏部分代碼省略.........

func check_level1_func(ind *linalg.IndexOpts, fn funcNum, X, Y matrix.Matrix) error {

	nX, nY := 0, 0
	// this is adapted from cvxopt:blas.c python blas interface
	switch fn {
	case fnrm2, fasum, fiamax, fscal, fset:
		if ind.IncX <= 0 {
			return errors.New("incX illegal, <=0")
		}
		if ind.OffsetX < 0 {
			return errors.New("offsetX illegal, <0")
		}
		sizeX := X.NumElements()
		if sizeX >= ind.OffsetX+1 {
			// calculate default size for N based on X size
			nX = 1 + (sizeX-ind.OffsetX-1)/ind.IncX
		}
		if sizeX < ind.OffsetX+1+(ind.Nx-1)*abs(ind.IncX) {
			return errors.New("X size error")
		}
		if ind.Nx < 0 {
			ind.Nx = nX
		}

	case fdot, fswap, fcopy, faxpy, faxpby:
		// vector X
		if ind.IncX <= 0 {
			return errors.New("incX illegal, <=0")
		}
		if ind.OffsetX < 0 {
			return errors.New("offsetX illegal, <0")
		}
		sizeX := X.NumElements()
		if sizeX >= ind.OffsetX+1 {
			// calculate default size for N based on X size
			nX = 1 + (sizeX-ind.OffsetX-1)/ind.IncX
		}
		if sizeX < ind.OffsetX+1+(ind.Nx-1)*abs(ind.IncX) {
			return errors.New("X size error")
		}
		if ind.Nx < 0 {
			ind.Nx = nX
		}
		// vector Y
		if ind.IncY <= 0 {
			return errors.New("incY illegal, <=0")
		}
		if ind.OffsetY < 0 {
			return errors.New("offsetY illegal, <0")
		}
		sizeY := Y.NumElements()
		if sizeY >= ind.OffsetY+1 {
			// calculate default size for N based on Y size
			nY = 1 + (sizeY-ind.OffsetY-1)/ind.IncY
		}
		if ind.Ny < 0 {
			ind.Ny = nY
		}
		if sizeY < ind.OffsetY+1+(ind.Ny-1)*abs(ind.IncY) {
			fmt.Printf("sizeY=%d, inds: %#v\n", sizeY, ind)
			return errors.New("Y size error")
		}

	case frotg, frotmg, frot, frotm:
	}
	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
//.........這裏部分代碼省略.........