// See function Symm.
func SymmFloat(A, B, C *matrix.FloatMatrix, alpha, beta float64, 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, fsymm, A, B, C, params)
	if err != nil {
		return
	}
	if ind.M == 0 || ind.N == 0 {
		return
	}
	Aa := A.FloatArray()
	Ba := B.FloatArray()
	Ca := C.FloatArray()
	uplo := linalg.ParamString(params.Uplo)
	side := linalg.ParamString(params.Side)
	dsymm(side, uplo, ind.M, ind.N, alpha, Aa[ind.OffsetA:], ind.LDa,
		Ba[ind.OffsetB:], ind.LDb, beta, Ca[ind.OffsetC:], ind.LDc)

	return
}

// See function Syrk2.
func Syr2kFloat(A, B, C *matrix.FloatMatrix, alpha, beta float64, 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, fsyr2k, A, B, C, params)
	if err != nil {
		return
	}
	if ind.N == 0 {
		return
	}
	Aa := A.FloatArray()
	Ba := B.FloatArray()
	Ca := C.FloatArray()
	uplo := linalg.ParamString(params.Uplo)
	trans := linalg.ParamString(params.Trans)
	//diag := linalg.ParamString(params.Diag)
	dsyr2k(uplo, trans, ind.N, ind.K, alpha, Aa[ind.OffsetA:], ind.LDa,
		Ba[ind.OffsetB:], ind.LDb, beta, Ca[ind.OffsetC:], ind.LDc)
	return
}

// See function Trsm.
func TrsmFloat(A, B *matrix.FloatMatrix, alpha float64, 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 ind.N == 0 || ind.M == 0 {
		return
	}
	Aa := A.FloatArray()
	Ba := B.FloatArray()
	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, alpha,
		Aa[ind.OffsetA:], ind.LDa, Ba[ind.OffsetB:], ind.LDb)
	return
}

// See function Gemm.
func GemmFloat(A, B, C *matrix.FloatMatrix, alpha, beta float64, 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
	}
	Aa := A.FloatArray()
	Ba := B.FloatArray()
	Ca := C.FloatArray()
	transB := linalg.ParamString(params.TransB)
	transA := linalg.ParamString(params.TransA)
	//diag := linalg.ParamString(params.Diag)
	dgemm(transA, transB, ind.M, ind.N, ind.K, alpha,
		Aa[ind.OffsetA:], ind.LDa, Ba[ind.OffsetB:], ind.LDb, beta,
		Ca[ind.OffsetC:], ind.LDc)
	return
}

func GesvdFloat(A, S, U, Vt *matrix.FloatMatrix, opts ...linalg.Option) error {
	pars, err := linalg.GetParameters(opts...)
	if err != nil {
		return err
	}
	ind := linalg.GetIndexOpts(opts...)
	err = checkGesvd(ind, pars, A, S, U, Vt)
	if err != nil {
		return err
	}
	if ind.M == 0 || ind.N == 0 {
		return nil
	}
	Aa := A.FloatArray()
	Sa := S.FloatArray()
	var Ua, Va []float64
	Ua = nil
	Va = nil
	if U != nil {
		Ua = U.FloatArray()[ind.OffsetU:]
	}
	if Vt != nil {
		Va = Vt.FloatArray()[ind.OffsetVt:]
	}
	info := dgesvd(linalg.ParamString(pars.Jobu), linalg.ParamString(pars.Jobvt),
		ind.M, ind.N, Aa[ind.OffsetA:], ind.LDa, Sa[ind.OffsetS:], Ua, ind.LDu, Va, ind.LDvt)
	if info != 0 {
		return onError(fmt.Sprintf("GesvdFloat lapack error: %d", info))
	}
	return nil
}

/*
 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 onError("Parameters not of same type")
	}
	switch A.(type) {
	case *matrix.FloatMatrix:
		Aa := A.(*matrix.FloatMatrix).FloatArray()
		Ba := B.(*matrix.FloatMatrix).FloatArray()
		Ca := C.(*matrix.FloatMatrix).FloatArray()
		aval := alpha.Float()
		bval := beta.Float()
		if math.IsNaN(aval) || math.IsNaN(bval) {
			return onError("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.(*matrix.ComplexMatrix).ComplexArray()
		Ba := B.(*matrix.ComplexMatrix).ComplexArray()
		Ca := C.(*matrix.ComplexMatrix).ComplexArray()
		aval := alpha.Complex()
		if cmplx.IsNaN(aval) {
			return onError("alpha not a number")
		}
		bval := beta.Complex()
		if cmplx.IsNaN(bval) {
			return onError("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 onError("Unknown type, not implemented")
	}
//.........这里部分代码省略.........

/*
 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 onError("Parameters not of same type")
	}
	switch A.(type) {
	case *matrix.FloatMatrix:
		Aa := A.(*matrix.FloatMatrix).FloatArray()
		Ca := C.(*matrix.FloatMatrix).FloatArray()
		aval := alpha.Float()
		bval := beta.Float()
		if math.IsNaN(aval) || math.IsNaN(bval) {
			return onError("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.(*matrix.ComplexMatrix).ComplexArray()
		Ca := C.(*matrix.ComplexMatrix).ComplexArray()
		aval := alpha.Complex()
		if cmplx.IsNaN(aval) {
			return onError("alpha not a real or complex number")
		}
		bval := beta.Float()
		if math.IsNaN(bval) {
			return onError("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 onError("Unknown type, not implemented")
	}

	return
}

func GbtrsFloat(A, B *matrix.FloatMatrix, ipiv []int32, KL int, opts ...linalg.Option) error {
	pars, err := linalg.GetParameters(opts...)
	if err != nil {
		return err
	}
	ind := linalg.GetIndexOpts(opts...)

	ind.Kl = KL
	err = checkGbtrs(ind, A, B, ipiv)
	if err != nil {
		return err
	}
	if ind.N == 0 || ind.Nrhs == 0 {
		return nil
	}
	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)
	if info != 0 {
		return onError(fmt.Sprintf("Gbtrs: lapack error: %d", info))
	}
	return nil
}

// See function Gbmv.
func GbmvFloat(A, X, Y *matrix.FloatMatrix, alpha, beta float64, 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
	}

	Xa := X.FloatArray()
	Ya := Y.FloatArray()
	Aa := A.FloatArray()
	if params.Trans == linalg.PNoTrans && ind.N == 0 {
		dscal(ind.M, beta, Ya[ind.OffsetY:], ind.IncY)
	} else if params.Trans == linalg.PTrans && ind.M == 0 {
		dscal(ind.N, beta, Ya[ind.OffsetY:], ind.IncY)
	} else {
		trans := linalg.ParamString(params.Trans)
		dgbmv(trans, ind.M, ind.N, ind.Kl, ind.Ku,
			alpha, Aa[ind.OffsetA:], ind.LDa, Xa[ind.OffsetX:], ind.IncX,
			beta, Ya[ind.OffsetY:], ind.IncY)
	}
	return
}

/*
 Solves a general real or complex set of linear equations.

 PURPOSE

 Solves A*X=B with A m by n real or complex.

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

 OPTIONS:
  trans
  m         nonnegative integer.  If negative, the default value is used.
  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.
*/
func Gels(A, B matrix.Matrix, opts ...linalg.Option) error {
	pars, _ := linalg.GetParameters(opts...)
	ind := linalg.GetIndexOpts(opts...)
	arows := ind.LDa
	brows := ind.LDb
	if ind.M < 0 {
		ind.M = A.Rows()
	}
	if ind.N < 0 {
		ind.N = A.Cols()
	}
	if ind.Nrhs < 0 {
		ind.Nrhs = B.Cols()
	}
	if ind.M == 0 || ind.N == 0 || ind.Nrhs == 0 {
		return nil
	}
	if ind.LDa == 0 {
		ind.LDa = max(1, A.LeadingIndex())
		arows = max(1, A.Rows())
	}
	if ind.LDa < max(1, ind.M) {
		return onError("Gesv: ldA")
	}
	if ind.LDb == 0 {
		ind.LDb = max(1, B.LeadingIndex())
		brows = max(1, B.Rows())
	}
	if ind.LDb < max(ind.M, ind.N) {
		return onError("Gesv: ldB")
	}
	if !matrix.EqualTypes(A, B) {
		return onError("Gesv: arguments not of same type")
	}
	_, _ = arows, brows // todo!! something
	info := -1
	trans := linalg.ParamString(pars.Trans)
	switch A.(type) {
	case *matrix.FloatMatrix:
		Aa := A.(*matrix.FloatMatrix).FloatArray()
		Ba := B.(*matrix.FloatMatrix).FloatArray()
		info = dgels(trans, ind.M, ind.N, ind.Nrhs, Aa[ind.OffsetA:], ind.LDa,
			Ba[ind.OffsetB:], ind.LDb)
	case *matrix.ComplexMatrix:
		Aa := A.(*matrix.ComplexMatrix).ComplexArray()
		Ba := B.(*matrix.ComplexMatrix).ComplexArray()
		info = zgels(trans, ind.M, ind.N, ind.Nrhs, Aa[ind.OffsetA:], ind.LDa,
			Ba[ind.OffsetB:], ind.LDb)
	}
	if info != 0 {
		return onError(fmt.Sprintf("Gels: lapack error: %d", info))
	}
	return nil
}

func SyevdFloat(A, W *matrix.FloatMatrix, opts ...linalg.Option) error {
	pars, err := linalg.GetParameters(opts...)
	if err != nil {
		return err
	}
	ind := linalg.GetIndexOpts(opts...)
	err = checkSyevd(ind, A, W)
	if err != nil {
		return err
	}
	if ind.N == 0 {
		return nil
	}
	jobz := linalg.ParamString(pars.Jobz)
	uplo := linalg.ParamString(pars.Uplo)
	Aa := A.FloatArray()
	Wa := W.FloatArray()
	info := dsyevd(jobz, uplo, ind.N, Aa[ind.OffsetA:], ind.LDa, Wa[ind.OffsetW:])
	if info != 0 {
		return onError(fmt.Sprintf("Syevd: lapack error %d", info))
	}
	return nil
}

/*
 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 onError("Parameters not of same type")
	}
	switch A.(type) {
	case *matrix.FloatMatrix:
		Aa := A.(*matrix.FloatMatrix).FloatArray()
		Ba := B.(*matrix.FloatMatrix).FloatArray()
		aval := alpha.Float()
		if math.IsNaN(aval) {
			return onError("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.(*matrix.ComplexMatrix).ComplexArray()
		Ba := B.(*matrix.ComplexMatrix).ComplexArray()
		aval := alpha.Complex()
		if cmplx.IsNaN(aval) {
			return onError("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 onError("Unknown type, not implemented")
	}
	return
}

// See function Tbsv.
func TbsvFloat(A, X *matrix.FloatMatrix, 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, ftbsv, X, nil, A, params)
	if err != nil {
		return
	}
	if ind.N == 0 {
		return
	}
	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)
	return
}

func PotrfFloat(A *matrix.FloatMatrix, opts ...linalg.Option) error {
	pars, err := linalg.GetParameters(opts...)
	if err != nil {
		return err
	}
	ind := linalg.GetIndexOpts(opts...)
	err = checkPotrf(ind, A)
	if ind.N == 0 {
		return nil
	}
	Aa := A.FloatArray()
	uplo := linalg.ParamString(pars.Uplo)
	info := dpotrf(uplo, ind.N, Aa[ind.OffsetA:], ind.LDa)
	if info != 0 {
		return onError(fmt.Sprintf("Potrf: lapack error %d", info))
	}
	return nil
}

// See function Syr.
func SyrFloat(X, A *matrix.FloatMatrix, alpha float64, 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, fsyr, X, nil, A, params)
	if err != nil {
		return
	}
	if ind.N == 0 {
		return
	}
	Xa := X.FloatArray()
	Aa := A.FloatArray()
	uplo := linalg.ParamString(params.Uplo)
	dsyr(uplo, ind.N, alpha, Xa[ind.OffsetX:], ind.IncX, Aa[ind.OffsetA:], ind.LDa)
	return
}

func PosvFloat(A, B *matrix.FloatMatrix, opts ...linalg.Option) error {
	pars, err := linalg.GetParameters(opts...)
	if err != nil {
		return err
	}
	ind := linalg.GetIndexOpts(opts...)
	err = checkPosv(ind, A, B)
	if err != nil {
		return err
	}
	if ind.N == 0 || ind.Nrhs == 0 {
		return nil
	}
	Aa := A.FloatArray()
	Ba := B.FloatArray()
	uplo := linalg.ParamString(pars.Uplo)
	info := dposv(uplo, ind.N, ind.Nrhs, Aa[ind.OffsetA:], ind.LDa, Ba[ind.OffsetB:], ind.LDb)
	if info != 0 {
		return onError(fmt.Sprintf("Posv: lapack error %d", info))
	}
	return nil
}

// See function Sbmv.
func SbmvFloat(A, X, Y *matrix.FloatMatrix, alpha, beta float64, 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
	}
	Xa := X.FloatArray()
	Ya := Y.FloatArray()
	Aa := A.FloatArray()
	uplo := linalg.ParamString(params.Uplo)
	dsbmv(uplo, ind.N, ind.K, alpha, Aa[ind.OffsetA:], ind.LDa, Xa[ind.OffsetX:],
		ind.IncX, beta, Ya[ind.OffsetY:], ind.IncY)
	return
}

//.........这里部分代码省略.........
  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.  Must have the same type as dl.
  B         float or complex matrix.  Must have the same type oas dl.
  ipiv      int vector

 OPTIONS
  trans     PNoTrans, PTrans, PConjTrans
  n         nonnegative integer.  If negative, the default value is used.
  nrhs      nonnegative integer.  If negative, the default value is used.
  ldB       positive integer, ldB >= max(1,n). If zero, the default value is used.
  offsetdl  nonnegative integer
  offsetd   nonnegative integer
  offsetdu  nonnegative integer
  offsetB   nonnegative integer

*/
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...)
	brows := ind.LDb
	if ind.OffsetD < 0 {
		return onError("Gttrs: offset D")
	}
	if ind.N < 0 {
		ind.N = D.NumElements() - ind.OffsetD
	}
	if ind.N < 0 {
		return onError("Gttrs: size D")
	}
	if ind.N == 0 {
		return nil
	}
	if ind.OffsetDL < 0 {
		return onError("Gttrs: offset DL")
	}
	sizeDL := DL.NumElements()
	if sizeDL < ind.OffsetDL+ind.N-1 {
		return onError("Gttrs: sizeDL")
	}
	if ind.OffsetDU < 0 {
		return onError("Gttrs: offset DU")
	}
	sizeDU := DU.NumElements()
	if sizeDU < ind.OffsetDU+ind.N-1 {
		return onError("Gttrs: sizeDU")
	}
	sizeDU2 := DU2.NumElements()
	if sizeDU2 < ind.N-2 {
		return onError("Gttrs: sizeDU2")
	}
	if ind.Nrhs < 0 {
		ind.Nrhs = B.Cols()
	}
	if ind.Nrhs == 0 {
		return nil
	}
	if ind.LDb == 0 {
		ind.LDb = max(1, B.LeadingIndex())
		brows = max(1, B.Rows())
	}
	if ind.LDb < max(1, ind.N) {
		return onError("Gttrs: ldB")
	}
	if ind.OffsetB < 0 {
		return onError("Gttrs: offset B")
	}
	sizeB := B.NumElements()
	if sizeB < ind.OffsetB+(ind.Nrhs-1)*brows+ind.N {
		return onError("Gttrs: sizeB")
	}
	if len(ipiv) < ind.N {
		return onError("Gttrs: size ipiv")
	}
	if !matrix.EqualTypes(DL, D, DU, DU2, B) {
		return onError("Gttrs: matrix types")
	}
	var info int = -1
	switch DL.(type) {
	case *matrix.FloatMatrix:
		DLa := DL.(*matrix.FloatMatrix).FloatArray()
		Da := D.(*matrix.FloatMatrix).FloatArray()
		DUa := DU.(*matrix.FloatMatrix).FloatArray()
		DU2a := DU2.(*matrix.FloatMatrix).FloatArray()
		Ba := B.(*matrix.FloatMatrix).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)
	case *matrix.ComplexMatrix:
		return onError("Gttrs: complex valued not yet implemented")
	}
	if info != 0 {
		return onError(fmt.Sprintf("Gttrs lapack error: %d", info))
	}
	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().

 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...)
	ind.Kl = KL
	arows := ind.LDa
	brows := ind.LDb
	if ind.Kl < 0 {
		return onError("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 onError("Gbtrs: invalid ku")
	}
	if ind.LDa == 0 {
		ind.LDa = max(1, A.LeadingIndex())
		arows = max(1, A.Rows())
	}
	if ind.LDa < 2*ind.Kl+ind.Ku+1 {
		return onError("Gbtrs: ldA")
	}
	if ind.OffsetA < 0 {
		return onError("Gbtrs: offsetA")
	}
	sizeA := A.NumElements()
	if sizeA < ind.OffsetA+(ind.N-1)*arows+2*ind.Kl+ind.Ku+1 {
		return onError("Gbtrs: sizeA")
	}
	if ind.LDb == 0 {
		ind.LDb = max(1, B.LeadingIndex())
		brows = max(1, B.Rows())
	}
	if ind.OffsetB < 0 {
		return onError("Gbtrs: offsetB")
	}
	sizeB := B.NumElements()
	if sizeB < ind.OffsetB+(ind.Nrhs-1)*brows+ind.N {
		return onError("Gbtrs: sizeB")
	}
	if ipiv != nil && len(ipiv) < ind.N {
		return onError("Gbtrs: size ipiv")
	}

	if !matrix.EqualTypes(A, B) {
		return onError("Gbtrs: arguments not of same type")
	}
	info := -1
	switch A.(type) {
	case *matrix.FloatMatrix:
		Aa := A.(*matrix.FloatMatrix).FloatArray()
		Ba := B.(*matrix.FloatMatrix).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:
//.........这里部分代码省略.........

/*
 Solves a real or complex symmetric set of linear equations,
 given the LDL^T factorization computed by sytrf() or sysv().

 PURPOSE
 Solves
  A*X = B

 where A is real or complex symmetric and n by n,
 and B is n by nrhs.  On entry, A and ipiv contain the
 factorization of A as returned by Sytrf() or Sysv().  On exit, B is
 replaced by the solution.

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

 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       nonnegative integer.  ldB >= max(1,n).  If zero, the
            default value is used.
  offsetA   nonnegative integer
  offsetB   nonnegative integer;

*/
func Sytrs(A, B matrix.Matrix, ipiv []int32, opts ...linalg.Option) error {
	pars, err := linalg.GetParameters(opts...)
	if err != nil {
		return err
	}
	ind := linalg.GetIndexOpts(opts...)
	arows := ind.LDa
	brows := ind.LDb
	if ind.N < 0 {
		ind.N = A.Rows()
		if ind.N != A.Cols() {
			return onError("Sytrs: 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.LeadingIndex())
		arows = max(1, A.Rows())
	}
	if ind.LDa < max(1, ind.N) {
		return onError("Sytrs: ldA")
	}
	if ind.LDb == 0 {
		ind.LDb = max(1, B.LeadingIndex())
		brows = max(1, B.Rows())
	}
	if ind.LDb < max(1, ind.N) {
		return onError("Sytrs: ldB")
	}
	if ind.OffsetA < 0 {
		return onError("Sytrs: offsetA")
	}
	sizeA := A.NumElements()
	if sizeA < ind.OffsetA+(ind.N-1)*arows+ind.N {
		return onError("Sytrs: sizeA")
	}
	if ind.OffsetB < 0 {
		return onError("Sytrs: offsetB")
	}
	sizeB := B.NumElements()
	if sizeB < ind.OffsetB+(ind.Nrhs-1)*brows+ind.N {
		return onError("Sytrs: sizeB")
	}
	if ipiv != nil && len(ipiv) < ind.N {
		return onError("Sytrs: size ipiv")
	}
	if !matrix.EqualTypes(A, B) {
		return onError("Sytrs: arguments not of same type")
	}
	info := -1
	switch A.(type) {
	case *matrix.FloatMatrix:
		Aa := A.(*matrix.FloatMatrix).FloatArray()
		Ba := B.(*matrix.FloatMatrix).FloatArray()
		uplo := linalg.ParamString(pars.Uplo)
		info = dsytrs(uplo, ind.N, ind.Nrhs, Aa[ind.OffsetA:], ind.LDa, ipiv,
			Ba[ind.OffsetB:], ind.LDb)
	case *matrix.ComplexMatrix:
		return onError("Sytrs: complex not yet implemented")
	}
	if info != 0 {
		return onError(fmt.Sprintf("Sytrs lapack error: %d", info))
	}
	return nil
}