Commit 1e362fda authored by Chris Cantwell's avatar Chris Cantwell

Merge branch 'feature/IncNSOpt' into 'master'

Improve performance of incompressible solver

This MR contains improvements to the performance of the incompressible solver (some of the changes might also help in other cases).

I worked on these changes based on a profiling of the `Cyl.xml` example (2D, N=5, using both tris and quads). For this particular case (in serial), I got a reduction of computational cost close to 30%, although I suspect simulations with higher order or in 3D will get more modest improvements.

The main changes are:

- Use a virtual function to make UnsteadySystem skip the FwdTrans it performs after each time step, since this is completely redundant for the incompressible solver
- Introduce a parameter to choose the frequency at which the solver checks for NaNs
- Tweak the default value of the mdswitch parameter for the multilevelstaticcond, and also allow to specify it in the session file
- Include PhysDeriv in a single direction in the Collections, since we need it to compute div(u) for the pressure forcing
- Use symmetric matrices for the interior part of the static cond (whenever possible)
- Several small changes to reduce overhead costs, especially in matrix vector multiplications.


See merge request !645
parents c976e139 07a2eddc
......@@ -75,6 +75,15 @@ class BwdTrans_StdMat : public Operator
0.0, output.get(), m_stdExp->GetTotPoints());
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
DNekMatSharedPtr m_mat;
......@@ -155,6 +164,15 @@ class BwdTrans_IterPerExp : public Operator
}
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
private:
BwdTrans_IterPerExp(
vector<StdRegions::StdExpansionSharedPtr> pCollExp,
......@@ -229,6 +247,15 @@ class BwdTrans_NoCollection : public Operator
}
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
vector<StdRegions::StdExpansionSharedPtr> m_expList;
......@@ -310,6 +337,15 @@ class BwdTrans_SumFac_Seg : public Operator
}
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
const int m_nquad0;
......@@ -406,6 +442,15 @@ class BwdTrans_SumFac_Quad : public Operator
}
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
const int m_nquad0;
const int m_nquad1;
......@@ -495,6 +540,15 @@ class BwdTrans_SumFac_Tri : public Operator
&output[0], m_nquad0);
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
const int m_nquad0;
......@@ -596,6 +650,15 @@ class BwdTrans_SumFac_Hex : public Operator
}
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
const int m_nquad0;
const int m_nquad1;
......@@ -755,6 +818,15 @@ class BwdTrans_SumFac_Tet : public Operator
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
const int m_nquad0;
const int m_nquad1;
......@@ -887,6 +959,15 @@ class BwdTrans_SumFac_Prism : public Operator
0.0, output.get(), m_nquad0);
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
const int m_nquad0;
......
......@@ -79,6 +79,12 @@ class Collection
Array<OneD, NekDouble> &output1,
Array<OneD, NekDouble> &output2);
inline void ApplyOperator(
const OperatorType &op,
int dir,
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &output);
inline bool HasOperator(const OperatorType &op);
protected:
......@@ -135,6 +141,19 @@ inline void Collection::ApplyOperator(
(*m_ops[op])(inarray, output0, output1, output2, wsp);
}
/**
*
*/
inline void Collection::ApplyOperator(
const OperatorType &op,
int dir,
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &output)
{
Array<OneD, NekDouble> wsp(m_ops[op]->GetWspSize());
(*m_ops[op])(dir, inarray, output, wsp);
}
inline bool Collection::HasOperator(const OperatorType &op)
{
return (m_ops.find(op) != m_ops.end());
......
......@@ -81,6 +81,15 @@ class IProductWRTBase_StdMat : public Operator
0.0, output.get(), m_stdExp->GetNcoeffs());
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
DNekMatSharedPtr m_mat;
Array<OneD, const NekDouble> m_jac;
......@@ -170,6 +179,15 @@ class IProductWRTBase_IterPerExp : public Operator
}
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
Array<OneD, NekDouble> m_jac;
......@@ -270,6 +288,15 @@ class IProductWRTBase_NoCollection : public Operator
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
vector<StdRegions::StdExpansionSharedPtr> m_expList;
......@@ -364,6 +391,15 @@ class IProductWRTBase_SumFac_Seg : public Operator
}
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
const int m_nquad0;
const int m_nmodes0;
......@@ -422,6 +458,15 @@ class IProductWRTBase_SumFac_Quad : public Operator
m_jac, input, output, wsp);
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
const int m_nquad0;
const int m_nquad1;
......@@ -487,6 +532,15 @@ class IProductWRTBase_SumFac_Tri : public Operator
output,wsp);
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
const int m_nquad0;
const int m_nquad1;
......@@ -561,6 +615,15 @@ class IProductWRTBase_SumFac_Hex : public Operator
m_jac,input,output,wsp);
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
const int m_nquad0;
const int m_nquad1;
......@@ -639,6 +702,15 @@ class IProductWRTBase_SumFac_Tet : public Operator
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
const int m_nquad0;
const int m_nquad1;
......@@ -722,6 +794,15 @@ class IProductWRTBase_SumFac_Prism : public Operator
m_jac,input,output,wsp);
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
const int m_nquad0;
const int m_nquad1;
......
......@@ -124,6 +124,15 @@ class IProductWRTDerivBase_StdMat : public Operator
}
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
Array<OneD, DNekMatSharedPtr> m_iProdWRTStdDBase;
Array<TwoD, const NekDouble> m_derivFac;
......@@ -284,6 +293,15 @@ class IProductWRTDerivBase_IterPerExp : public Operator
}
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
Array<TwoD, const NekDouble> m_derivFac;
Array<OneD, const NekDouble> m_jac;
......@@ -403,6 +421,15 @@ class IProductWRTDerivBase_NoCollection : public Operator
}
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
int m_dim;
int m_coordim;
......@@ -498,6 +525,15 @@ class IProductWRTDerivBase_SumFac_Seg : public Operator
&output[0], m_nmodes0);
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
const int m_nquad0;
const int m_nmodes0;
......@@ -591,6 +627,15 @@ class IProductWRTDerivBase_SumFac_Quad : public Operator
Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
const int m_nquad0;
const int m_nquad1;
......@@ -730,6 +775,15 @@ class IProductWRTDerivBase_SumFac_Tri : public Operator
Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
const int m_nquad0;
const int m_nquad1;
......@@ -890,6 +944,15 @@ class IProductWRTDerivBase_SumFac_Hex : public Operator
Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
const int m_nquad0;
const int m_nquad1;
......@@ -1077,6 +1140,15 @@ class IProductWRTDerivBase_SumFac_Tet : public Operator
Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
const int m_nquad0;
const int m_nquad1;
......@@ -1284,6 +1356,15 @@ class IProductWRTDerivBase_SumFac_Prism : public Operator
Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
const int m_nquad0;
......
......@@ -127,6 +127,13 @@ class Operator
Array<OneD, NekDouble> &wsp
= NullNekDouble1DArray) = 0;
COLLECTIONS_EXPORT virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp
= NullNekDouble1DArray) = 0;
COLLECTIONS_EXPORT virtual ~Operator();
/// Get the size of the required workspace
......
This diff is collapsed.
......@@ -79,7 +79,7 @@ namespace Blas
void F77NAME(dtpmv) (const char& uplo, const char& trans, const char& diag,
const int& n, const double* ap, double* x, const int& incx);
void F77NAME(dspmv) (const char& trans, const int& n, const double& alpha,
void F77NAME(dspmv) (const char& uplo, const int& n, const double& alpha,
const double* a, const double* x, const int& incx,
const double& beta, double* y, const int& incy);
......@@ -192,11 +192,11 @@ namespace Blas
/// \brief BLAS level 2: Matrix vector multiply y = A \e x where A
/// is symmetric packed
static inline void Dspmv (const char& trans, const int& n, const double& alpha,
static inline void Dspmv (const char& uplo, const int& n, const double& alpha,
const double* a, const double* x, const int& incx,
const double& beta, double* y, const int& incy)
{
F77NAME(dspmv) (trans,n,alpha,a,x,incx,beta,y,incy);
F77NAME(dspmv) (uplo,n,alpha,a,x,incx,beta,y,incy);
}
static inline void Dsbmv (const char& uplo, const int& m, const int& k,
......
......@@ -40,14 +40,13 @@ namespace Nektar
{
template<typename DataType, typename InnerMatrixType>
NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::NekMatrix(MatrixStorage type) :
BaseType(0,0),
BaseType(0,0,type),
m_data(),
m_rowSizes(),
m_columnSizes(),
m_storageSize(),
m_numberOfBlockRows(0),
m_numberOfBlockColumns(0),
m_storageType(type)
m_numberOfBlockColumns(0)
{
}
......@@ -55,14 +54,13 @@ namespace Nektar
NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::NekMatrix(unsigned int numberOfBlockRows, unsigned int numberOfBlockColumns,
unsigned int rowsPerBlock, unsigned int columnsPerBlock,
MatrixStorage type) :
BaseType(numberOfBlockRows*rowsPerBlock, numberOfBlockColumns*columnsPerBlock),
BaseType(numberOfBlockRows*rowsPerBlock, numberOfBlockColumns*columnsPerBlock,type),
m_data(),
m_rowSizes(numberOfBlockRows),
m_columnSizes(numberOfBlockColumns),
m_storageSize(0),
m_numberOfBlockRows(numberOfBlockRows),
m_numberOfBlockColumns(numberOfBlockColumns),
m_storageType(type)
m_numberOfBlockColumns(numberOfBlockColumns)
{
m_storageSize = GetRequiredStorageSize();
m_data = Array<OneD, boost::shared_ptr<InnerType> >(m_storageSize, boost::shared_ptr<InnerType>());
......@@ -82,14 +80,14 @@ namespace Nektar
const unsigned int* rowsPerBlock, const unsigned int* columnsPerBlock,
MatrixStorage type) :
BaseType(std::accumulate(rowsPerBlock, rowsPerBlock + numberOfBlockRows, 0),
std::accumulate(columnsPerBlock, columnsPerBlock + numberOfBlockColumns, 0)),
std::accumulate(columnsPerBlock, columnsPerBlock + numberOfBlockColumns, 0),
type),
m_data(),
m_rowSizes(numberOfBlockRows),
m_columnSizes(numberOfBlockColumns),
m_storageSize(0),
m_numberOfBlockRows(numberOfBlockRows),
m_numberOfBlockColumns(numberOfBlockColumns),
m_storageType(type)
m_numberOfBlockColumns(numberOfBlockColumns)
{
m_storageSize = GetRequiredStorageSize();
m_data = Array<OneD, boost::shared_ptr<InnerType> >(m_storageSize, boost::shared_ptr<InnerType>());
......@@ -101,14 +99,14 @@ namespace Nektar
const Array<OneD, const unsigned int>& rowsPerBlock, const Array<OneD, const unsigned int>& columnsPerBlock,
MatrixStorage type) :
BaseType(std::accumulate(rowsPerBlock.data(), rowsPerBlock.data() + numberOfBlockRows, 0),
std::accumulate(columnsPerBlock.data(), columnsPerBlock.data() + numberOfBlockColumns, 0)),
std::accumulate(columnsPerBlock.data(), columnsPerBlock.data() + numberOfBlockColumns, 0),
type),
m_data(),
m_rowSizes(numberOfBlockRows),
m_columnSizes(numberOfBlockColumns),
m_storageSize(0),
m_numberOfBlockRows(numberOfBlockRows),
m_numberOfBlockColumns(numberOfBlockColumns),
m_storageType(type)
m_numberOfBlockColumns(numberOfBlockColumns)
{
m_storageSize = GetRequiredStorageSize();
m_data = Array<OneD, boost::shared_ptr<InnerType> >(m_storageSize, boost::shared_ptr<InnerType>());
......@@ -120,14 +118,14 @@ namespace Nektar
const Array<OneD, const unsigned int>& columnsPerBlock,
MatrixStorage type) :
BaseType(std::accumulate(rowsPerBlock.begin(), rowsPerBlock.end(), 0),
std::accumulate(columnsPerBlock.begin(), columnsPerBlock.end(), 0)),
std::accumulate(columnsPerBlock.begin(), columnsPerBlock.end(), 0),
type),
m_data(),
m_rowSizes(rowsPerBlock.num_elements()),
m_columnSizes(columnsPerBlock.num_elements()),
m_storageSize(0),
m_numberOfBlockRows(rowsPerBlock.num_elements()),
m_numberOfBlockColumns(columnsPerBlock.num_elements()),
m_storageType(type)
m_numberOfBlockColumns(columnsPerBlock.num_elements())
{
m_storageSize = GetRequiredStorageSize();
m_data = Array<OneD, boost::shared_ptr<InnerType> >(m_storageSize, boost::shared_ptr<InnerType>());
......@@ -142,8 +140,7 @@ namespace Nektar
m_columnSizes(rhs.m_columnSizes),
m_storageSize(rhs.m_storageSize),
m_numberOfBlockRows(rhs.m_numberOfBlockRows),
m_numberOfBlockColumns(rhs.m_numberOfBlockColumns),
m_storageType(rhs.m_storageType)
m_numberOfBlockColumns(rhs.m_numberOfBlockColumns)
{
for(unsigned int i = 0; i < rhs.m_data.num_elements(); ++i)
{
......@@ -161,16 +158,8 @@ namespace Nektar
template<typename DataType, typename InnerMatrixType>
unsigned int NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::CalculateBlockIndex(unsigned int row, unsigned int column) const
{
unsigned int blockRows = GetNumberOfBlockRows();
unsigned int blockCols = GetNumberOfBlockColumns();
if( this->GetTransposeFlag() == 'T' )
{
std::swap(blockRows, blockCols);
}
return BaseType::CalculateIndex(this->GetStorageType(),
row, column, blockRows, blockCols, this->GetTransposeFlag());
row, column, m_numberOfBlockRows, m_numberOfBlockColumns, this->GetTransposeFlag());
}
template<typename DataType, typename InnerMatrixType>
......@@ -303,12 +292,6 @@ namespace Nektar
return m_storageSize;
}
template<typename DataType, typename InnerMatrixType>
MatrixStorage NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::GetType() const
{
return m_storageType;
}
template<typename DataType, typename InnerMatrixType>
unsigned int NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::GetNumberOfBlockRows() const
{
......@@ -361,6 +344,23 @@ namespace Nektar
}
}
template<typename DataType, typename InnerMatrixType>
void NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::GetBlockSizes(
Array<OneD, unsigned int>& rowSizes,
Array<OneD, unsigned int>& colSizes) const
{
if( this->GetTransposeFlag() == 'T' )
{
rowSizes = m_columnSizes;
colSizes = m_rowSizes;
}
else
{
rowSizes = m_rowSizes;
colSizes = m_columnSizes;
}
}
template<typename DataType, typename InnerMatrixType>
typename NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::iterator NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::begin() { return iterator(*this, 0, 0); }
......@@ -438,12 +438,6 @@ namespace Nektar
return this->GetStorageSize();
}
template<typename DataType, typename InnerMatrixType>
MatrixStorage NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::v_GetStorageType() const
{
return this->GetType();
}
template<typename DataType, typename InnerMatrixType>
void NekMatrix<NekMatrix<DataType, InnerMatrixType>, BlockMatrixTag>::v_Transpose()
{
......
......@@ -167,8 +167,6 @@ namespace Nektar
LIB_UTILITIES_EXPORT unsigned int GetStorageSize() const;
LIB_UTILITIES_EXPORT MatrixStorage GetType() const;
LIB_UTILITIES_EXPORT unsigned int GetNumberOfBlockRows() const;
LIB_UTILITIES_EXPORT unsigned int GetNumberOfBlockColumns() const;
......@@ -177,6 +175,9 @@ namespace Nektar
LIB_UTILITIES_EXPORT unsigned int GetNumberOfColumnsInBlockColumn(unsigned int blockCol) const;
LIB_UTILITIES_EXPORT void GetBlockSizes(Array<OneD, unsigned int>& rowSizes,
Array<OneD, unsigned int>& colSizes) const;
LIB_UTILITIES_EXPORT iterator begin();
LIB_UTILITIES_EXPORT iterator end();
LIB_UTILITIES_EXPORT const_iterator begin() const;
......@@ -197,8 +198,6 @@ namespace Nektar
LIB_UTILITIES_EXPORT virtual unsigned int v_GetStorageSize() const;
LIB_UTILITIES_EXPORT virtual MatrixStorage v_GetStorageType() const;
LIB_UTILITIES_EXPORT virtual void v_Transpose();
Array<OneD, boost::shared_ptr<InnerType> > m_data;
......@@ -208,7 +207,6 @@ namespace Nektar
unsigned int m_storageSize;
unsigned int m_numberOfBlockRows;
unsigned int m_numberOfBlockColumns;
MatrixStorage m_storageType;
static NumberType m_zeroElement;
};
......
......@@ -48,6 +48,10 @@ namespace Lapack
const int& nrhs, const double* ap,
const int *ipiv, double* b,
const int& ldb, int& info);
void F77NAME(dsptri) (const char& uplo, const int& n,
const double* ap,
const int *ipiv, double* work,
int& info);
void F77NAME(dtrtrs) (const char& uplo, const char& trans, const char& diag,