Golang Matrix.FloatArray方法代码示例

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

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

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

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

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

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

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

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

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

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

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

/*
 Matrix-vector product with a real symmetric or complex hermitian band matrix.

 Computes with A real symmetric and  banded of order n and with bandwidth k.
  Y := alpha*A*X + beta*Y

 ARGUMENTS
  A         float or complex n*n matrix
  X         float or complex n*1 matrix
  Y         float or complex n*1 matrix
  alpha     number (float or complex singleton matrix)
  beta      number (float or complex singleton matrix)

 OPTIONS
  uplo      PLower or PUpper
  n         integer.  If negative, the default value is used.
  k         integer.  If negative, the default value is used.
            The default value is k = max(0,A.Rows()-1).
  ldA       nonnegative integer.  ldA >= k+1.
            If zero, the default vaule is used.
  incx      nonzero integer
  incy      nonzero integer
  offsetA   nonnegative integer
  offsetx   nonnegative integer
  offsety   nonnegative integer

*/
func Hbmv(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, fsbmv, X, Y, A, params)
	if err != nil {
		return
	}
	if 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")
		}
		uplo := linalg.ParamString(params.Uplo)
		dsbmv(uplo, ind.N, ind.K, aval, Aa[ind.OffsetA:], ind.LDa,
			Xa[ind.OffsetX:], ind.IncX, bval, Ya[ind.OffsetY:], ind.IncY)

	case *matrix.ComplexMatrix:
		Xa := X.ComplexArray()
		Ya := Y.ComplexArray()
		Aa := A.ComplexArray()
		aval := alpha.Complex()
		bval := beta.Complex()
		uplo := linalg.ParamString(params.Uplo)
		zhbmv(uplo, ind.N, ind.K, aval, Aa[ind.OffsetA:], ind.LDa,
			Xa[ind.OffsetX:], ind.IncX, bval, Ya[ind.OffsetY:], ind.IncY)
		//zhbmv(uplo, ind.N, aval, Aa[ind.OffsetA:], ind.LDa,
		//	Xa[ind.OffsetX:], ind.IncX,
		//	bval, Ya[ind.OffsetY:], ind.IncY)
	default:
		return errors.New("Unknown type, not implemented")
	}
	return
}