/*
 Solution of a triangular and banded set of equations.

 Tbsv(A, X, uplo=PLower, trans=PNoTrans, diag=PNonDiag, n=A.Cols,
 k=max(0,A.Rows-1), ldA=A.size[0], incx=1, offsetA=0, offsetx=0)

PURPOSE
  X := A^{-1}*X, if trans is PNoTrans
  X := A^{-T}*X, if trans is PTrans
  X := A^{-H}*X, if trans is PConjTrans

 A is banded triangular of order n and with bandwidth k.

 ARGUMENTS
  A         float or complex m*k matrix.
  X         float or complex k*1 matrix. Must have the same type as A.

 OPTIONS
  uplo      PLower   or PUpper
  trans     PNoTrans, PTrans or PConjTrans
  diag      PNoNUnit or PUnit
  n         nonnegative integer.  If negative, the default value is used.
  k         nonnegative integer.  If negative, the default value is used.
  ldA       nonnegative integer.  ldA >= 1+k.
            If zero the default value is used.
  incx      nonzero integer
  offsetA   nonnegative integer
  offsetx   nonnegative integer;
*/
func Tbsv(A, X matrix.Matrix, opts ...linalg.Option) (err error) {

	var params *linalg.Parameters
	if !matrix.EqualTypes(A, X) {
		err = errors.New("Parameters not of same type")
		return
	}
	params, err = linalg.GetParameters(opts...)
	if err != nil {
		return
	}
	ind := linalg.GetIndexOpts(opts...)
	err = check_level2_func(ind, ftbsv, X, nil, A, params)
	if err != nil {
		return
	}
	if ind.N == 0 {
		return
	}
	switch X.(type) {
	case *matrix.FloatMatrix:
		Xa := X.FloatArray()
		Aa := A.FloatArray()
		uplo := linalg.ParamString(params.Uplo)
		trans := linalg.ParamString(params.Trans)
		diag := linalg.ParamString(params.Diag)
		dtbsv(uplo, trans, diag, ind.N, ind.K,
			Aa[ind.OffsetA:], ind.LDa, Xa[ind.OffsetX:], ind.IncX)
	case *matrix.ComplexMatrix:
		return errors.New("Not implemented yet for complx.Matrix")
	default:
		return errors.New("Unknown type, not implemented")
	}
	return
}

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
}

/*
 Symmetric rank-2 update.
 syr2(x, y, A, uplo='L', alpha=1.0, n=A.size[0], incx=1, incy=1,
     ldA=max(1,A.size[0]), offsetx=0, offsety=0, offsetA=0)
 PURPOSE
 Computes A := A + alpha*(x*y^T + y*x^T) with A real symmetric matrix of order n.
 ARGUMENTS
 x         float matrix
 y         float matrix
 A         float matrix
 alpha     real number (int or float)

 OPTIONS
 uplo      'L' or 'U'
 n         integer.  If negative, the default value is used.
 incx      nonzero integer
 incy      nonzero integer
 ldA       nonnegative integer.  ldA >= max(1,n).
           If zero the default value is used.
 offsetx   nonnegative integer
 offsety   nonnegative integer
 offsetA   nonnegative integer;
*/
func Syr2(X, Y, A matrix.Matrix, alpha matrix.Scalar, opts ...linalg.Option) (err error) {

	var params *linalg.Parameters
	params, err = linalg.GetParameters(opts...)
	if err != nil {
		return
	}
	ind := linalg.GetIndexOpts(opts...)
	err = check_level2_func(ind, fsyr2, X, Y, A, params)
	if err != nil {
		return
	}
	if !matrix.EqualTypes(A, X, Y) {
		return errors.New("Parameters not of same type")
	}
	switch X.(type) {
	case *matrix.FloatMatrix:
		Xa := X.FloatArray()
		Ya := X.FloatArray()
		Aa := A.FloatArray()
		aval := alpha.Float()
		if math.IsNaN(aval) {
			return errors.New("alpha not a number")
		}
		uplo := linalg.ParamString(params.Uplo)
		dsyr2(uplo, ind.N, aval, Xa[ind.OffsetX:], ind.IncX,
			Ya[ind.OffsetY:], ind.IncY,
			Aa[ind.OffsetA:], ind.LDa)
	case *matrix.ComplexMatrix:
		return errors.New("Not implemented yet for complx.Matrix")
	default:
		return errors.New("Unknown type, not implemented")
	}
	return
}

// Scales a vector by a constant (X := alpha*X).
//
// ARGUMENTS
//  X         float or complex matrix
//  alpha     number (float or complex singleton matrix).  Complex alpha is only
//            allowed if X is complex.
// 
// OPTIONS
//  n         integer.  If n<0, the default value of n is used.
//            The default value is equal to 1+(len(x)-offset-1)/inc or 0
//            if len(x) > offset+1.
//  inc       positive integer, default = 1
//  offset    nonnegative integer, default = 0
//
func Scal(X matrix.Matrix, alpha matrix.Scalar, opts ...linalg.Option) (err error) {
	ind := linalg.GetIndexOpts(opts...)
	err = check_level1_func(ind, fscal, X, nil)
	if err != nil {
		return
	}
	if ind.Nx == 0 {
		return
	}
	switch X.(type) {
	case *matrix.ComplexMatrix:
		Xa := X.ComplexArray()
		cval := alpha.Complex()
		zscal(ind.Nx, cval, Xa[ind.OffsetX:], ind.IncX)
	case *matrix.FloatMatrix:
		Xa := X.FloatArray()
		rval := alpha.Float()
		if math.IsNaN(rval) {
			return errors.New("alpha not float value")
		}
		dscal(ind.Nx, rval, Xa[ind.OffsetX:], ind.IncX)
	default:
		err = errors.New("not implemented for parameter types", )
	}
	return
}

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

// Returns Y = X^H*Y for real or complex X, Y.
//
// ARGUMENTS
//  X         float or complex matrix
//  Y         float or complex matrix.  Must have the same type as X.
//
// OPTIONS
//  n         integer.  If n<0, the default value of n is used.
//            The default value is equal to nx = 1+(len(x)-offsetx-1)/incx or 0 if
//            len(x) > offsetx+1.  If the default value is used, it must be equal to
//            ny = 1+(len(y)-offsetx-1)/|incy| or 0 if len(y) > offsety+1
//  incx      nonzero integer [default=1]
//  incy      nonzero integer [default=1]
//  offsetx   nonnegative integer [default=0]
//  offsety   nonnegative integer [default=0]
//
func Dot(X, Y matrix.Matrix, opts ...linalg.Option) (v matrix.Scalar) {
	v = matrix.FScalar(math.NaN())
	//cv = cmplx.NaN()
	ind := linalg.GetIndexOpts(opts...)
	err := check_level1_func(ind, fdot, X, Y)
	if err != nil {
		return
	}
	if ind.Nx == 0 {
		return matrix.FScalar(0.0)
	}
	sameType := matrix.EqualTypes(X, Y)
	if ! sameType {
		err = errors.New("arrays not of same type")
		return
	}
	switch X.(type) {
	case *matrix.ComplexMatrix:
		Xa := X.ComplexArray()
		Ya := Y.ComplexArray()
		v = matrix.CScalar(zdotc(ind.Nx, Xa[ind.OffsetX:], ind.IncX, Ya[ind.OffsetY:], ind.IncY))
	case *matrix.FloatMatrix:
		Xa := X.FloatArray()
		Ya := Y.FloatArray()
		v = matrix.FScalar(ddot(ind.Nx, Xa[ind.OffsetX:], ind.IncX, Ya[ind.OffsetY:], ind.IncY))
	//default:
	//	err = errors.New("not implemented for parameter types", )
	}
	return
}

// Copies a vector X to a vector Y (Y := X).
//
// ARGUMENTS
//  X         float or complex matrix
//  Y         float or complex matrix.  Must have the same type as X.
//
// OPTIONS
//  n         integer.  If n<0, the default value of n is used.
//            The default value is given by 1+(len(x)-offsetx-1)/incx or 0
//            if len(x) > offsetx+1
//  incx      nonzero integer
//  incy      nonzero integer
//  offsetx   nonnegative integer
//  offsety   nonnegative integer;
//
func Copy(X, Y matrix.Matrix, opts ...linalg.Option) (err error) {
	ind := linalg.GetIndexOpts(opts...)
	err = check_level1_func(ind, fcopy, X, Y)
	if err != nil {
		return
	}
	if ind.Nx == 0 {
		return
	}
	sameType := matrix.EqualTypes(X, Y)
	if ! sameType {
		err = errors.New("arrays not same type")
		return
	}
	switch X.(type) {
	case *matrix.ComplexMatrix:
		Xa := X.ComplexArray()
		Ya := Y.ComplexArray()
		zcopy(ind.Nx, Xa[ind.OffsetX:], ind.IncX, Ya[ind.OffsetY:], ind.IncY)
	case *matrix.FloatMatrix:
		Xa := X.FloatArray()
		Ya := Y.FloatArray()
		dcopy(ind.Nx, Xa[ind.OffsetX:], ind.IncX, Ya[ind.OffsetY:], ind.IncY)
	default:
		err = errors.New("not implemented for parameter types", )
	}
	return
}

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
}

// Returns ||Re x||_1 + ||Im x||_1.
//
// ARGUMENTS
//  X       float or complex matrix
// 
// OPTIONS
//  n       integer.  If n<0, the default value of n is used.
//          The default value is equal to n = 1+(len(x)-offset-1)/inc or 0 if
//          len(x) > offset+1
//  inc     positive integer
//  offset  nonnegative integer
//
func Asum(X matrix.Matrix, opts ...linalg.Option) (v matrix.Scalar) {
	v = matrix.FScalar(math.NaN())
	ind := linalg.GetIndexOpts(opts...)
	err := check_level1_func(ind, fasum, X, nil)
	if err != nil {
		return
	}
	if ind.Nx == 0 {
		return
	}
	switch X.(type) {
	case *matrix.ComplexMatrix:
		Xa := X.ComplexArray()
		v = matrix.FScalar(dzasum(ind.Nx, Xa[ind.OffsetX:], ind.IncX))
	case *matrix.FloatMatrix:
		Xa := X.FloatArray()
		v =  matrix.FScalar(dasum(ind.Nx, Xa[ind.OffsetX:], ind.IncX))
	//default:
	//	err = errors.New("not implemented for parameter types", )
	}
	return
}

/*
 Rank-k update of symmetric matrix. (L3)

 Herk(A, C, alpha, beta, uplo=PLower, trans=PNoTrans,  n=-1,
 k=-1, ldA=max(1,A.Rows), ldC=max(1,C.Rows), offsetA=0, offsetB=0)

 Computes
  C := alpha*A*A^T + beta*C, if trans is PNoTrans
  C := alpha*A^T*A + beta*C, if trans is PTrans

 C is symmetric (real or complex) of order n. The inner dimension of the matrix
 product is k.  If k=0 this is interpreted as C := beta*C.

 ARGUMENTS
  A         float or complex matrix.
  C         float or complex matrix.  Must have the same type as A.
  alpha     number (float or complex singleton matrix).  Complex alpha is only
            allowed if A is complex.
  beta      number (float or complex singleton matrix).  Complex beta is only
            allowed if A is complex.

 OPTIONS
  uplo      PLower or PUpper
  trans     PNoTrans or PTrans
  n         integer.  If negative, the default value is used.
            The default value is n = A.Rows or if trans == PNoTrans n = A.Cols.
  k         integer.  If negative, the default value is used.
            The default value is k =  A.Cols, or if trans == PNoTrans k = A.Rows.
  ldA       nonnegative integer.
            ldA >= max(1,n) or if trans != PNoTrans ldA >= max(1,k).
            If zero, the default value is used.
  ldC       nonnegative integer.  ldC >= max(1,n).
            If zero, the default value is used.
  offsetA   nonnegative integer
  offsetC   nonnegative integer;
*/
func Herk(A, C matrix.Matrix, alpha, beta matrix.Scalar, opts ...linalg.Option) (err error) {

	params, e := linalg.GetParameters(opts...)
	if e != nil {
		err = e
		return
	}
	ind := linalg.GetIndexOpts(opts...)
	err = check_level3_func(ind, fsyrk, A, nil, C, params)
	if e != nil || err != nil {
		return
	}
	if !matrix.EqualTypes(A, C) {
		return errors.New("Parameters not of same type")
	}
	switch A.(type) {
	case *matrix.FloatMatrix:
		Aa := A.FloatArray()
		Ca := C.FloatArray()
		aval := alpha.Float()
		bval := beta.Float()
		if math.IsNaN(aval) || math.IsNaN(bval) {
			return errors.New("alpha or beta not a number")
		}
		uplo := linalg.ParamString(params.Uplo)
		trans := linalg.ParamString(params.Trans)
		dsyrk(uplo, trans, ind.N, ind.K, aval, Aa[ind.OffsetA:], ind.LDa, bval,
			Ca[ind.OffsetC:], ind.LDc)
	case *matrix.ComplexMatrix:
		Aa := A.ComplexArray()
		Ca := C.ComplexArray()
		aval := alpha.Complex()
		if cmplx.IsNaN(aval) {
			return errors.New("alpha not a real or complex number")
		}
		bval := beta.Float()
		if math.IsNaN(bval) {
			return errors.New("beta not a real number")
		}
		uplo := linalg.ParamString(params.Uplo)
		trans := linalg.ParamString(params.Trans)
		zherk(uplo, trans, ind.N, ind.K, aval, Aa[ind.OffsetA:], ind.LDa, bval,
			Ca[ind.OffsetC:], ind.LDc)
	default:
		return errors.New("Unknown type, not implemented")
	}

	return
}

/*
 Solution of a triangular system of equations with multiple righthand sides. (L3)

 Trsm(A, B, alpha, side=PLeft, uplo=PLower, transA=PNoTrans, diag=PNonUnit,
 m=-1, n=-1, ldA=max(1,A.Rows), ldB=max(1,B.Rows), offsetA=0, offsetB=0)

 Computes
  B := alpha*A^{-1}*B if transA is PNoTrans   and side = PLeft
  B := alpha*B*A^{-1} if transA is PNoTrans   and side = PRight
  B := alpha*A^{-T}*B if transA is PTrans     and side = PLeft
  B := alpha*B*A^{-T} if transA is PTrans     and side = PRight
  B := alpha*A^{-H}*B if transA is PConjTrans and side = PLeft
  B := alpha*B*A^{-H} if transA is PConjTrans and side = PRight

 B is m by n and A is triangular.  The code does not verify whether A is nonsingular.

 ARGUMENTS
  A         float or complex matrix.
  B         float or complex matrix.  Must have the same type as A.
  alpha     number (float or complex).  Complex alpha is only
            allowed if A is complex.

 OPTIONS
  side      PLeft or PRight
  uplo      PLower or PUpper
  transA    PNoTrans or PTrans
  diag      PNonUnit or PUnit
  m         integer.  If negative, the default value is used.
            The default value is m = A.Rows or if side == PRight m = B.Rows
            If the default value is used and side is PLeft, m must be equal to A.Cols.
  n         integer.  If negative, the default value is used.
            The default value is n = B.Cols or if side )= PRight n = A.Rows.
            If the default value is used and side is PRight, n must be equal to A.Cols.
  ldA       nonnegative integer.
            ldA >= max(1,m) of if  side == PRight lda >= max(1,n).
            If zero, the default value is used.
  ldB       nonnegative integer.  ldB >= max(1,m).
            If zero, the default value is used.
  offsetA   nonnegative integer
  offsetB   nonnegative integer
*/
func Trsm(A, B matrix.Matrix, alpha matrix.Scalar, opts ...linalg.Option) (err error) {

	params, e := linalg.GetParameters(opts...)
	if e != nil {
		err = e
		return
	}
	ind := linalg.GetIndexOpts(opts...)
	err = check_level3_func(ind, ftrsm, A, B, nil, params)
	if err != nil {
		return
	}
	if !matrix.EqualTypes(A, B) {
		return errors.New("Parameters not of same type")
	}
	switch A.(type) {
	case *matrix.FloatMatrix:
		Aa := A.FloatArray()
		Ba := B.FloatArray()
		aval := alpha.Float()
		if math.IsNaN(aval) {
			return errors.New("alpha or beta not a number")
		}
		uplo := linalg.ParamString(params.Uplo)
		transA := linalg.ParamString(params.TransA)
		side := linalg.ParamString(params.Side)
		diag := linalg.ParamString(params.Diag)
		dtrsm(side, uplo, transA, diag, ind.M, ind.N, aval,
			Aa[ind.OffsetA:], ind.LDa, Ba[ind.OffsetB:], ind.LDb)
	case *matrix.ComplexMatrix:
		Aa := A.ComplexArray()
		Ba := B.ComplexArray()
		aval := alpha.Complex()
		if cmplx.IsNaN(aval) {
			return errors.New("alpha  not a number")
		}
		uplo := linalg.ParamString(params.Uplo)
		transA := linalg.ParamString(params.TransA)
		side := linalg.ParamString(params.Side)
		diag := linalg.ParamString(params.Diag)
		ztrsm(side, uplo, transA, diag, ind.M, ind.N, aval,
			Aa[ind.OffsetA:], ind.LDa, Ba[ind.OffsetB:], ind.LDb)
	default:
		return errors.New("Unknown type, not implemented")
	}
	return
}

// Constant times a vector plus a vector (Y := alpha*X+Y).
//
// ARGUMENTS
//   X         float or complex matrix
//   Y         float or complex matrix.  Must have the same type as X.
//   alpha     number (float or complex singleton matrix).  Complex alpha is only
//             allowed if x is complex.
//
// OPTIONS
//   n         integer.  If n<0, the default value of n is used.
//             The default value is equal to 1+(len(x)-offsetx-1)/incx 
//             or 0 if  len(x) >= offsetx+1 
//   incx      nonzero integer
//   incy      nonzero integer
//   offsetx   nonnegative integer
//   offsety   nonnegative integer;
//
func Axpy(X, Y matrix.Matrix, alpha matrix.Scalar, opts ...linalg.Option) (err error) {
	ind := linalg.GetIndexOpts(opts...)
	err = check_level1_func(ind, faxpy, X, Y)
	if err != nil {
		return
	}
	if ind.Nx == 0 {
		return
	}
	sameType := matrix.EqualTypes(X, Y)
	if ! sameType {
		err = errors.New("arrays not same type")
		return
	}
	switch X.(type) {
	case *matrix.ComplexMatrix:
		Xa := X.ComplexArray()
		Ya := Y.ComplexArray()
		aval := alpha.Complex()
		if cmplx.IsNaN(aval) {
			return errors.New("alpha not complex value")
		}
		zaxpy(ind.Nx, aval, Xa[ind.OffsetX:],
			ind.IncX, Ya[ind.OffsetY:], ind.IncY)
	case *matrix.FloatMatrix:
		Xa := X.FloatArray()
		Ya := Y.FloatArray()
		aval := alpha.Float()
		if math.IsNaN(aval) {
			return errors.New("alpha not float value")
		}
		daxpy(ind.Nx, aval, Xa[ind.OffsetX:],
			ind.IncX, Ya[ind.OffsetY:], ind.IncY)
	default:
		err = errors.New("not implemented for parameter types", )
	}
	return
}

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
}

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
}

/*
 Matrix-vector product with a general banded matrix. (L2)

 Computes
   Y := alpha*A*X + beta*Y,   if trans = PNoTrans
   Y := alpha*A^T*X + beta*Y, if trans = PTrans
   Y := beta*y,               if n=0, m>0, and trans = PNoTrans
   Y := beta*y,               if n>0, m=0, and trans = PTrans

 The matrix A is m by n with upper bandwidth ku and lower bandwidth kl.
 Returns immediately if n=0 and trans is 'Trans', or if m=0 and trans is 'N'.


 ARGUMENTS
   X         float n*1 matrix.
   Y         float m*1 matrix
   A         float m*n matrix.
   alpha     number (float).
   beta      number (float).

 OPTIONS
   trans     NoTrans or Trans
   m         nonnegative integer, default A.Rows()
   kl        nonnegative integer
   n         nonnegative integer.  If negative, the default value is used.
   ku        nonnegative integer.  If negative, the default value is used.
   ldA       positive integer.  ldA >= kl+ku+1. If zero, the default value is used.
   incx      nonzero integer, default =1
   incy      nonzero integer, default =1
   offsetA   nonnegative integer, default =0
   offsetx   nonnegative integer, default =0
   offsety   nonnegative integer, default =0

*/
func Gbmv(A, X, Y matrix.Matrix, alpha, beta matrix.Scalar, opts ...linalg.Option) (err error) {

	var params *linalg.Parameters
	params, err = linalg.GetParameters(opts...)
	if err != nil {
		return
	}
	ind := linalg.GetIndexOpts(opts...)
	err = check_level2_func(ind, fgbmv, X, Y, A, params)
	if err != nil {
		return
	}
	if ind.M == 0 && ind.N == 0 {
		return
	}
	if !matrix.EqualTypes(A, X, Y) {
		return errors.New("Parameters not of same type")
	}
	switch X.(type) {
	case *matrix.FloatMatrix:
		Xa := X.FloatArray()
		Ya := Y.FloatArray()
		Aa := A.FloatArray()
		aval := alpha.Float()
		bval := beta.Float()
		if math.IsNaN(aval) || math.IsNaN(bval) {
			return errors.New("alpha or beta not a number")
		}
		if params.Trans == linalg.PNoTrans && ind.N == 0 {
			dscal(ind.M, bval, Ya[ind.OffsetY:], ind.IncY)
		} else if params.Trans == linalg.PTrans && ind.M == 0 {
			dscal(ind.N, bval, Ya[ind.OffsetY:], ind.IncY)
		} else {
			trans := linalg.ParamString(params.Trans)
			dgbmv(trans, ind.M, ind.N, ind.Kl, ind.Ku,
				aval, Aa[ind.OffsetA:], ind.LDa, Xa[ind.OffsetX:], ind.IncX,
				bval, Ya[ind.OffsetY:], ind.IncY)
		}
	case *matrix.ComplexMatrix:
		return errors.New("Not implemented yet for complx.Matrix")
	default:
		return errors.New("Unknown type, not implemented")
	}
	return
}

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
//.........这里部分代码省略.........

/*
 General matrix-matrix product. (L3)

 PURPOSE
 Computes
  C := alpha*A*B + beta*C     if transA = PNoTrans   and transB = PNoTrans.
  C := alpha*A^T*B + beta*C   if transA = PTrans     and transB = PNoTrans.
  C := alpha*A^H*B + beta*C   if transA = PConjTrans and transB = PNoTrans.
  C := alpha*A*B^T + beta*C   if transA = PNoTrans   and transB = PTrans.
  C := alpha*A^T*B^T + beta*C if transA = PTrans     and transB = PTrans.
  C := alpha*A^H*B^T + beta*C if transA = PConjTrans and transB = PTrans.
  C := alpha*A*B^H + beta*C   if transA = PNoTrans   and transB = PConjTrans.
  C := alpha*A^T*B^H + beta*C if transA = PTrans     and transB = PConjTrans.
  C := alpha*A^H*B^H + beta*C if transA = PConjTrans and transB = PConjTrans.

 The number of rows of the matrix product is m.  The number of  columns is n.
 The inner dimension is k.  If k=0, this reduces  to C := beta*C.

 ARGUMENTS
  A         float or complex matrix, m*k
  B         float or complex matrix, k*n
  C         float or complex matrix, m*n
  alpha     number (float or complex singleton matrix)
  beta      number (float or complex singleton matrix)

 OPTIONS
  transA    PNoTrans, PTrans or PConjTrans
  transB    PNoTrans, PTrans or PConjTrans
  m         integer.  If negative, the default value is used. The default value is
            m = A.Rows of if transA != PNoTrans m = A.Cols.
  n         integer.  If negative, the default value is used. The default value is
            n = (transB == PNoTrans) ? B.Cols : B.Rows.
  k         integer.  If negative, the default value is used. The default value is
            k=A.Cols or if transA != PNoTrans) k = A.Rows, transA=PNoTrans.
            If the default value is used it should also be equal to
            (transB == PNoTrans) ? B.Rows : B.Cols.
  ldA       nonnegative integer.  ldA >= max(1,m) of if transA != NoTrans max(1,k).
            If zero, the default value is used.
  ldB       nonnegative integer.  ldB >= max(1,k) or if transB != NoTrans max(1,n).
            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
  offsetC   nonnegative integer;
*/
func Gemm(A, B, C matrix.Matrix, alpha, beta matrix.Scalar, opts ...linalg.Option) (err error) {

	params, e := linalg.GetParameters(opts...)
	if e != nil {
		err = e
		return
	}
	ind := linalg.GetIndexOpts(opts...)
	err = check_level3_func(ind, fgemm, A, B, C, params)
	if err != nil {
		return
	}
	if ind.M == 0 || ind.N == 0 {
		return
	}
	if !matrix.EqualTypes(A, B, C) {
		return errors.New("Parameters not of same type")
	}
	switch A.(type) {
	case *matrix.FloatMatrix:
		Aa := A.FloatArray()
		Ba := B.FloatArray()
		Ca := C.FloatArray()
		aval := alpha.Float()
		bval := beta.Float()
		if math.IsNaN(aval) || math.IsNaN(bval) {
			return errors.New("alpha or beta not a number")
		}
		transB := linalg.ParamString(params.TransB)
		transA := linalg.ParamString(params.TransA)
		dgemm(transA, transB, ind.M, ind.N, ind.K, aval,
			Aa[ind.OffsetA:], ind.LDa, Ba[ind.OffsetB:], ind.LDb, bval,
			Ca[ind.OffsetC:], ind.LDc)

	case *matrix.ComplexMatrix:
		Aa := A.ComplexArray()
		Ba := B.ComplexArray()
		Ca := C.ComplexArray()
		aval := alpha.Complex()
		if cmplx.IsNaN(aval) {
			return errors.New("alpha not a number")
		}
		bval := beta.Complex()
		if cmplx.IsNaN(bval) {
			return errors.New("beta not a number")
		}
		transB := linalg.ParamString(params.TransB)
		transA := linalg.ParamString(params.TransA)
		zgemm(transA, transB, ind.M, ind.N, ind.K, aval,
			Aa[ind.OffsetA:], ind.LDa, Ba[ind.OffsetB:], ind.LDb, bval,
			Ca[ind.OffsetC:], ind.LDc)
	default:
		return errors.New("Unknown type, not implemented")
	}
//.........这里部分代码省略.........

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")
			}
//.........这里部分代码省略.........

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