source: Sophya/trunk/SophyaExt/LinAlg/intflapack.cc@ 825

#include <iostream.h>
#include "intflapack.h"
#include <typeinfo>

extern "C" {
void sgesv_(int_4* n, int_4* nrhs, r_4* a, int_4* lda, 
            int_4* ipiv, r_4* b, int_4* ldb, int_4* info);
void dgesv_(int_4* n, int_4* nrhs, r_8* a, int_4* lda, 
            int_4* ipiv, r_8* b, int_4* ldb, int_4* info);
void cgesv_(int_4* n, int_4* nrhs, complex<r_4>* a, int_4* lda, 
            int_4* ipiv, complex<r_4>* b, int_4* ldb, int_4* info);
void zgesv_(int_4* n, int_4* nrhs, complex<r_8>* a, int_4* lda, 
            int_4* ipiv, complex<r_8>* b, int_4* ldb, int_4* info);
}

template <class T>
TArray<T>& LapackServer<T>::LinSolve(TArray<T>& a, TArray<T> & b)
{
  if ( ( a.NbDimensions() != 2 ) || ( b.NbDimensions() != 2 ) )
    throw(SzMismatchError("LapackServer::LinSolve(a,b) a Or b NbDimensions() != 2"));
  
  uint_4 rowa = a.RowsKA();
  uint_4 cola = a.ColsKA();
  uint_4 rowb = b.RowsKA();
  uint_4 colb = b.ColsKA();
  if ( a.Size(rowa) !=  a.Size(cola)) 
    throw(SzMismatchError("LapackServer::LinSolve(a,b) a Not a square Array"));
  if ( a.Size(rowa) !=  b.Size(rowa)) 
    throw(SzMismatchError("LapackServer::LinSolve(a,b) RowSize(a <> b) "));

  if (!a.IsPacked(rowa) || !b.IsPacked(rowb))
     throw(SzMismatchError("LapackServer::LinSolve(a,b) a Or b Not Packed Columns "));

  int_4 n = a.Size(rowa);
  int_4 nrhs = b.Size(colb);
  int_4 lda = a.Step(cola);
  int_4 ldb = b.Step(colb);
  int_4 info;
  int_4* ipiv = new int_4[n];

  if (typeid(T) == typeid(r_4) )
    sgesv_(&n, &nrhs, (r_4 *)a.Data(), &lda, ipiv, (r_4 *)b.Data(), &ldb, &info);
  else if (typeid(T) == typeid(r_8) ) 
    dgesv_(&n, &nrhs, (r_8 *)a.Data(), &lda, ipiv, (r_8 *)b.Data(), &ldb, &info);
  else if (typeid(T) == typeid(complex<r_4>) ) 
    cgesv_(&n, &nrhs, (complex<r_4> *)a.Data(), &lda, ipiv, 
           (complex<r_4> *)b.Data(), &ldb, &info);
  else if (typeid(T) == typeid(complex<r_8>) ) 
    zgesv_(&n, &nrhs, (complex<r_8> *)a.Data(), &lda, ipiv, 
           (complex<r_8> *)b.Data(), &ldb, &info);
  else {
    delete[] ipiv;
    string tn = typeid(T).name();
    cerr << " LapackServer::LinSolve(a,b) - Unsupported DataType T = " << tn << endl;
    throw TypeMismatchExc("LapackServer::LinSolve(a,b) - Unsupported DataType (T)");
  }
  delete[] ipiv;
  return(b);
}

void rztest_lapack(TArray<r_4>& aa, TArray<r_4>& bb)
{
  if ( aa.NbDimensions() != 2 ) throw(SzMismatchError("rztest_lapack(TMatrix<r_4> A Not a Matrix"));
  if ( aa.SizeX() !=  aa.SizeY()) throw(SzMismatchError("rztest_lapack(TMatrix<r_4> A Not a square Matrix"));
  if ( bb.NbDimensions() != 2 ) throw(SzMismatchError("rztest_lapack(TMatrix<r_4> A Not a Matrix"));
  if ( bb.SizeX() !=  aa.SizeX() ) throw(SzMismatchError("rztest_lapack(TMatrix<r_4> A <> B "));
  if ( !bb.IsPacked() || !bb.IsPacked() ) 
    throw(SzMismatchError("rztest_lapack(TMatrix<r_4> Not packed A or B "));
  
  int_4 n = aa.SizeX();
  int_4 nrhs = bb.SizeY();
  int_4 lda = n;
  int_4 ldb = bb.SizeX();
  int_4 info;
  int_4* ipiv = new int_4[n];
  sgesv_(&n, &nrhs, aa.Data(), &lda, ipiv, bb.Data(), &ldb, &info);
  delete[] ipiv;
  cout << "rztest_lapack/Info= " << info << endl;
  cout << aa << "\n" << bb << endl;
  return;
}

///////////////////////////////////////////////////////////////
#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template LapackServer<r_4>
#pragma define_template LapackServer<r_8>
#pragma define_template LapackServer< complex<r_4> >
#pragma define_template LapackServer< complex<r_8> >
#endif

#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
template class LapackServer<r_4>;
template class LapackServer<r_8>;
template class LapackServer< complex<r_4> >;
template class LapackServer< complex<r_8> >;
#endif

#if defined(OS_LINUX)
//  Pour le link avec f2c sous Linux
extern "C" {
    void MAIN__();
}

void MAIN__()
{
  cerr << "MAIN__() function for linking with libf2c.a " << endl;
  cerr << "  This function should never be called !!! " << endl;
  throw PError("MAIN__() should not be called - see intflapack.cc");
}
#endif