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

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

/*************** Pour memoire  (Christophe) ***************
Les dispositions memoires (FORTRAN) pour les vecteurs et matrices LAPACK:

1./ --- REAL X(N):
    if an array X of dimension (N) holds a vector x,
    then X(i) holds "x_i" for i=1,...,N

2./ --- REAL A(LDA,N):
    if a two-dimensional array A of dimension (LDA,N) holds an m-by-n matrix A,
    then A(i,j) holds "a_ij" for i=1,...,m et j=1,...,n  (LDA must be at least m).
    Note that array arguments are usually declared in the software as assumed-size
    arrays (last dimension *), for example: REAL A(LDA,*)
    --- Rangement en memoire:
                            | 11 12 13 14 |
    Ex: Real A(4,4):    A = | 21 22 23 24 |
                            | 31 32 33 34 |
                            | 41 42 43 44 |
        memoire: {11 21 31 41} {12 22 32 42} {13 23 33 43} {14 24 34 44}
                 First indice (line) "i" varies then the second (column):
                 (put all the first column, then put all the second column,
                  ..., then put all the last column)
***********************************************************/

/*!
   \defgroup LinAlg LinAlg module
   This module contains classes and functions for complex linear
   algebra on arrays. This module is intended mainly to have 
   classes implementing C++ interfaces between Sophya objects
   and external linear algebra libraries, such as LAPACK.  
*/

/*! 
  \class SOPHYA::LapackServer
  \ingroup LinAlg
  This class implements an interface to LAPACK library driver routines.
  The LAPACK (Linear Algebra PACKage) is a collection high performance 
  routines to solve common problems in numerical linear algebra. 
  its is available from http://www.netlib.org.
  
  The present version of our LapackServer (Feb 2001) provides only 
  interfaces for the linear system solver and singular value 
  decomposition (SVD). Only arrays with BaseArray::FortranMemoryMapping 
  can be handled by LapackServer. LapackServer can be instanciated 
  for simple and double precision real or complex array types.

  The example below shows solving a linear system A*X = B

  \code
  #include "intflapack.h"
  // ...
  // Use FortranMemoryMapping as default
  BaseArray::SetDefaultMemoryMapping(BaseArray::FortranMemoryMapping);
  // Create an fill the arrays A and B
  int n = 20;
  Matrix A(n, n);
  A = RandomSequence();
  Vector X(n),B(n);
  X = RandomSequence();
  B = A*X;
  // Solve the linear system A*X = B
  LapackServer<r_8> lps;
  lps.LinSolve(A,B);
  // We get the result in B, which should be equal to X ...
  // Compute the difference B-X ;
  Vector diff = B-X;
  \endcode

*/

////////////////////////////////////////////////////////////////////////////////////
extern "C" {
// Le calculateur de workingspace
  int_4 ilaenv_(int_4 *ispec,char *name,char *opts,int_4 *n1,int_4 *n2,int_4 *n3,int_4 *n4,
                int_4 nc1,int_4 nc2);

// Drivers pour resolution de systemes lineaires
  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);

// Drivers pour resolution de systemes lineaires symetriques
  void ssysv_(char* uplo, int_4* n, int_4* nrhs, r_4* a, int_4* lda, 
              int_4* ipiv, r_4* b, int_4* ldb,
              r_4* work, int_4* lwork, int_4* info);
  void dsysv_(char* uplo, int_4* n, int_4* nrhs, r_8* a, int_4* lda, 
              int_4* ipiv, r_8* b, int_4* ldb,
              r_8* work, int_4* lwork, int_4* info);
  void csysv_(char* uplo, int_4* n, int_4* nrhs, complex<r_4>* a, int_4* lda, 
              int_4* ipiv, complex<r_4>* b, int_4* ldb,
              complex<r_4>* work, int_4* lwork, int_4* info);
  void zsysv_(char* uplo, int_4* n, int_4* nrhs, complex<r_8>* a, int_4* lda, 
              int_4* ipiv, complex<r_8>* b, int_4* ldb,
              complex<r_8>* work, int_4* lwork, int_4* info);

// Driver pour resolution de systemes au sens de Xi2
  void sgels_(char * trans, int_4* m, int_4* n, int_4* nrhs, r_4* a, int_4* lda, 
              r_4* b, int_4* ldb, r_4* work, int_4* lwork, int_4* info);
  void dgels_(char * trans, int_4* m, int_4* n, int_4* nrhs, r_8* a, int_4* lda, 
              r_8* b, int_4* ldb, r_8* work, int_4* lwork, int_4* info);
  void cgels_(char * trans, int_4* m, int_4* n, int_4* nrhs, complex<r_4>* a, int_4* lda, 
              complex<r_4>* b, int_4* ldb, complex<r_4>* work, int_4* lwork, int_4* info);
  void zgels_(char * trans, int_4* m, int_4* n, int_4* nrhs, complex<r_8>* a, int_4* lda, 
              complex<r_8>* b, int_4* ldb, complex<r_8>* work, int_4* lwork, int_4* info);

// Driver pour decomposition SVD 
  void sgesvd_(char* jobu, char* jobvt, int_4* m, int_4* n, r_4* a, int_4* lda,
               r_4* s, r_4* u, int_4* ldu, r_4* vt, int_4* ldvt, 
               r_4* work, int_4* lwork, int_4* info);
  void dgesvd_(char* jobu, char* jobvt, int_4* m, int_4* n, r_8* a, int_4* lda,
               r_8* s, r_8* u, int_4* ldu, r_8* vt, int_4* ldvt, 
               r_8* work, int_4* lwork, int_4* info);
  void cgesvd_(char* jobu, char* jobvt, int_4* m, int_4* n, complex<r_4>* a, int_4* lda,
               r_4* s, complex<r_4>* u, int_4* ldu, complex<r_4>* vt, int_4* ldvt, 
               complex<r_4>* work, int_4* lwork, r_4* rwork, int_4* info);
  void zgesvd_(char* jobu, char* jobvt, int_4* m, int_4* n, complex<r_8>* a, int_4* lda,
               r_8* s, complex<r_8>* u, int_4* ldu, complex<r_8>* vt, int_4* ldvt, 
               complex<r_8>* work, int_4* lwork, r_8* rwork, int_4* info);

// Driver pour decomposition SVD Divide and Conquer
  void sgesdd_(char* jobz, int_4* m, int_4* n, r_4* a, int_4* lda,
               r_4* s, r_4* u, int_4* ldu, r_4* vt, int_4* ldvt, 
               r_4* work, int_4* lwork, int_4* iwork, int_4* info);
  void dgesdd_(char* jobz, int_4* m, int_4* n, r_8* a, int_4* lda,
               r_8* s, r_8* u, int_4* ldu, r_8* vt, int_4* ldvt, 
               r_8* work, int_4* lwork, int_4* iwork, int_4* info);
  void cgesdd_(char* jobz, int_4* m, int_4* n, complex<r_4>* a, int_4* lda,
               r_4* s, complex<r_4>* u, int_4* ldu, complex<r_4>* vt, int_4* ldvt, 
               complex<r_4>* work, int_4* lwork, r_4* rwork, int_4* iwork, int_4* info);
  void zgesdd_(char* jobz, int_4* m, int_4* n, complex<r_8>* a, int_4* lda,
               r_8* s, complex<r_8>* u, int_4* ldu, complex<r_8>* vt, int_4* ldvt, 
               complex<r_8>* work, int_4* lwork, r_8* rwork, int_4* iwork, int_4* info);

// Driver pour eigen decomposition for symetric/hermitian matrices
  void ssyev_(char* jobz, char* uplo, int_4* n, r_4* a, int_4* lda, r_4* w,
              r_4* work, int_4 *lwork, int_4* info);
  void dsyev_(char* jobz, char* uplo, int_4* n, r_8* a, int_4* lda, r_8* w,
              r_8* work, int_4 *lwork, int_4* info);
  void cheev_(char* jobz, char* uplo, int_4* n, complex<r_4>* a, int_4* lda, r_4* w,
              complex<r_4>* work, int_4 *lwork, r_4* rwork, int_4* info);
  void zheev_(char* jobz, char* uplo, int_4* n, complex<r_8>* a, int_4* lda, r_8* w,
              complex<r_8>* work, int_4 *lwork, r_8* rwork, int_4* info);

// Driver pour eigen decomposition for general squared matrices
  void sgeev_(char* jobl, char* jobvr, int_4* n, r_4* a, int_4* lda, r_4* wr, r_4* wi,
              r_4* vl, int_4* ldvl, r_4* vr, int_4* ldvr, 
              r_4* work, int_4 *lwork, int_4* info);
  void dgeev_(char* jobl, char* jobvr, int_4* n, r_8* a, int_4* lda, r_8* wr, r_8* wi,
              r_8* vl, int_4* ldvl, r_8* vr, int_4* ldvr, 
              r_8* work, int_4 *lwork, int_4* info);
  void cgeev_(char* jobl, char* jobvr, int_4* n,  complex<r_4>* a, int_4* lda, complex<r_4>* w,
              complex<r_4>* vl, int_4* ldvl, complex<r_4>* vr, int_4* ldvr, 
              complex<r_4>* work, int_4 *lwork, r_4* rwork, int_4* info);
  void zgeev_(char* jobl, char* jobvr, int_4* n,  complex<r_8>* a, int_4* lda, complex<r_8>* w,
              complex<r_8>* vl, int_4* ldvl, complex<r_8>* vr, int_4* ldvr, 
              complex<r_8>* work, int_4 *lwork, r_8* rwork, int_4* info);

}

//   -------------- Classe LapackServer<T> --------------

////////////////////////////////////////////////////////////////////////////////////
template <class T>
LapackServer<T>::LapackServer()
{
  SetWorkSpaceSizeFactor();
}

template <class T>
LapackServer<T>::~LapackServer()
{
}

// --- ATTENTION BUG POSSIBLE dans l'avenir (CMV) --- REZA A LIRE S.T.P.
// -> Cette connerie de Fortran/C interface
// Dans les routines fortran de lapack:
// Appel depuis le C avec:
//   int_4 lwork = -1;
// SUBROUTINE SSYSV( UPLO,N,NRHS,A,LDA,IPIV,B,LDB,WORK,LWORK,INFO)
//   INTEGER            INFO, LDA, LDB, LWORK, N, NRHS
//   LOGICAL            LQUERY
//   LQUERY = ( LWORK.EQ.-1 )
//   ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
//    ==> le test est bien interprete sous Linux mais pas sous OSF
//    ==> Sous OSF "LWORK.EQ.-1" est FALSE quand on passe lwork=-1 par argument
//    ==> POUR REZA: confusion entier 4 / 8 bits ??? (bizarre on l'aurait vu avant?)
////////////////////////////////////////////////////////////////////////////////////
template <class T>
int_4 LapackServer<T>::ilaenv_en_C(int_4 ispec,char *name,char *opts,int_4 n1,int_4 n2,int_4 n3,int_4 n4)
{
 int_4 nc1 = strlen(name);
 int_4 nc2 = strlen(opts);
 int_4 rc=0;
 rc = ilaenv_(&ispec,name,opts,&n1,&n2,&n3,&n4,nc1,nc2);
 //cout<<"ilaenv_en_C("<<ispec<<","<<name<<"("<<nc1<<"),"<<opts<<"("<<nc2<<"),"
 //    <<n1<<","<<n2<<","<<n3<<","<<n4<<") = "<<rc<<endl;
 return rc;
}

////////////////////////////////////////////////////////////////////////////////////
//! Interface to Lapack linear system solver driver s/d/c/zgesvd().
/*! Solve the linear system a * x = b. Input arrays
  should have FortranMemory mapping (column packed).
  \param a : input matrix, overwritten on output
  \param b : input-output, input vector b, contains x on exit
  \return : return code from lapack driver _gesv()
 */
template <class T>
int 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"));
  
  int_4 rowa = a.RowsKA();
  int_4 cola = a.ColsKA();
  int_4 rowb = b.RowsKA();
  int_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(rowb)) 
    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 Column Packed"));

  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(info);
}

////////////////////////////////////////////////////////////////////////////////////
//! Interface to Lapack linear system solver driver s/d/c/zsysvd().
/*! Solve the linear system a * x = b with a symetric. Input arrays
  should have FortranMemory mapping (column packed).
  \param a : input matrix symetric , overwritten on output
  \param b : input-output, input vector b, contains x on exit
  \return : return code from lapack driver
 */
template <class T>
int LapackServer<T>::LinSolveSym(TArray<T>& a, TArray<T> & b)
// --- REMARQUES DE CMV ---
// 1./ contrairement a ce qui est dit dans la doc, il s'agit
//     de matrices SYMETRIQUES complexes et non HERMITIENNES !!!
// 2./ pourquoi les routines de LinSolve pour des matrices symetriques
//     sont plus de deux fois plus lentes que les LinSolve generales sur OSF
//     et sensiblement plus lentes sous Linux ???
{
  if ( ( a.NbDimensions() != 2 ) || ( b.NbDimensions() != 2 ) )
    throw(SzMismatchError("LapackServer::LinSolveSym(a,b) a Or b NbDimensions() != 2"));
  int_4 rowa = a.RowsKA();
  int_4 cola = a.ColsKA();
  int_4 rowb = b.RowsKA();
  int_4 colb = b.ColsKA();
  if ( a.Size(rowa) !=  a.Size(cola)) 
    throw(SzMismatchError("LapackServer::LinSolveSym(a,b) a Not a square Array"));
  if ( a.Size(rowa) !=  b.Size(rowb)) 
    throw(SzMismatchError("LapackServer::LinSolveSym(a,b) RowSize(a <> b) "));

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

  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 = 0;
  int_4* ipiv = new int_4[n];
  int_4 lwork = -1;
  T * work = NULL;

  char uplo = 'U'; // char uplo = 'L';
  char struplo[5]; struplo[0] = uplo; struplo[1] = '\0';

  if (typeid(T) == typeid(r_4) ) {
    lwork = ilaenv_en_C(1,"SSYTRF",struplo,n,-1,-1,-1) * n  +5;
    work = new T[lwork];
    ssysv_(&uplo, &n, &nrhs, (r_4 *)a.Data(), &lda, ipiv, (r_4 *)b.Data(), &ldb,
          (r_4 *)work, &lwork, &info);
  } else if (typeid(T) == typeid(r_8) )  {
    lwork = ilaenv_en_C(1,"DSYTRF",struplo,n,-1,-1,-1) * n  +5;
    work = new T[lwork];
    dsysv_(&uplo, &n, &nrhs, (r_8 *)a.Data(), &lda, ipiv, (r_8 *)b.Data(), &ldb,
          (r_8 *)work, &lwork, &info);
  } else if (typeid(T) == typeid(complex<r_4>) )  {
    lwork = ilaenv_en_C(1,"CSYTRF",struplo,n,-1,-1,-1) * n  +5;
    work = new T[lwork];
    csysv_(&uplo, &n, &nrhs, (complex<r_4> *)a.Data(), &lda, ipiv,
          (complex<r_4> *)b.Data(), &ldb,
          (complex<r_4> *)work, &lwork, &info);
  } else if (typeid(T) == typeid(complex<r_8>) )  {
    lwork = ilaenv_en_C(1,"ZSYTRF",struplo,n,-1,-1,-1) * n  +5;
    work = new T[lwork];
    zsysv_(&uplo, &n, &nrhs, (complex<r_8> *)a.Data(), &lda, ipiv,
          (complex<r_8> *)b.Data(), &ldb,
          (complex<r_8> *)work, &lwork, &info);
  } else {
    if(work) delete[] work;
    delete[] ipiv;
    string tn = typeid(T).name();
    cerr << " LapackServer::LinSolveSym(a,b) - Unsupported DataType T = " << tn << endl;
    throw TypeMismatchExc("LapackServer::LinSolveSym(a,b) - Unsupported DataType (T)");
  }
  if(work) delete[] work;
  delete[] ipiv;
  return(info);
}

////////////////////////////////////////////////////////////////////////////////////
//! Interface to Lapack least squares solver driver s/d/c/zgels().
/*! Solves the linear least squares problem defined by an m-by-n matrix 
  \b a and an m element vector \b b .
  A solution \b x to the overdetermined system of linear equations
  b = a * x is computed, minimizing the norm of b-a*x. 
  Underdetermined systems (m<n) are not yet handled.
  Inout arrays should have FortranMemory mapping (column packed).
  \param a : input matrix, overwritten on output
  \param b : input-output, input vector b, contains x on exit.
  \return : return code from lapack driver _gels()
  \warning : b is not resized.
 */
/*
  $CHECK$ - A faire - cas m<n 
       If the linear system is underdetermined, the minimum norm 
       solution is computed. 
*/

template <class T>
int LapackServer<T>::LeastSquareSolve(TArray<T>& a, TArray<T> & b)
{
  if ( ( a.NbDimensions() != 2 ) || ( b.NbDimensions() != 2 ) )
    throw(SzMismatchError("LapackServer::LeastSquareSolve(a,b) a Or b NbDimensions() != 2"));
  
  int_4 rowa = a.RowsKA();
  int_4 cola = a.ColsKA();
  int_4 rowb = b.RowsKA();
  int_4 colb = b.ColsKA();


  if ( a.Size(rowa) !=  b.Size(rowb)) 
    throw(SzMismatchError("LapackServer::LeastSquareSolve(a,b) RowSize(a <> b) "));

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

  if ( a.Size(rowa) <  a.Size(cola)) {  // $CHECK$ - m<n a changer
    cout << " LapackServer<T>::LeastSquareSolve() - m<n - Not yet implemented for "
         << " underdetermined systems ! " << endl;
    throw(SzMismatchError("LapackServer::LeastSquareSolve(a,b) NRows<NCols - "));
  }
  int_4 m = a.Size(rowa);
  int_4 n = a.Size(cola);
  int_4 nrhs = b.Size(colb);

  int_4 lda = a.Step(cola);
  int_4 ldb = b.Step(colb);
  int_4 info;

  int_4 minmn = (m < n) ? m : n;
  int_4 maxmn = (m > n) ? m : n;
  int_4 maxmnrhs = (nrhs > maxmn) ? nrhs : maxmn;
  if (maxmnrhs < 1) maxmnrhs = 1;
  
  int_4 lwork = minmn+maxmnrhs*5;
  T * work = new T[lwork];

  char trans = 'N';
  
  if (typeid(T) == typeid(r_4) )
    sgels_(&trans, &m, &n, &nrhs, (r_4 *)a.Data(), &lda, 
           (r_4 *)b.Data(), &ldb, (r_4 *)work, &lwork, &info);
  else if (typeid(T) == typeid(r_8) ) 
    dgels_(&trans, &m, &n, &nrhs, (r_8 *)a.Data(), &lda, 
           (r_8 *)b.Data(), &ldb, (r_8 *)work, &lwork, &info);
  else if (typeid(T) == typeid(complex<r_4>) ) 
    cgels_(&trans, &m, &n, &nrhs, (complex<r_4> *)a.Data(), &lda, 
           (complex<r_4> *)b.Data(), &ldb, (complex<r_4> *)work, &lwork, &info);
  else if (typeid(T) == typeid(complex<r_8>) ) 
    zgels_(&trans, &m, &n, &nrhs, (complex<r_8> *)a.Data(), &lda, 
           (complex<r_8> *)b.Data(), &ldb, (complex<r_8> *)work, &lwork, &info);
  else {
    delete[] work;
    string tn = typeid(T).name();
    cerr << " LapackServer::LeastSquareSolve(a,b) - Unsupported DataType T = " << tn << endl;
    throw TypeMismatchExc("LapackServer::LeastSquareSolve(a,b) - Unsupported DataType (T)");
  }
  delete[] work;
  return(info);
}


////////////////////////////////////////////////////////////////////////////////////
//! Interface to Lapack SVD driver s/d/c/zgesv().
/*! Computes the vector of singular values of \b a. Input arrays
  should have FortranMemoryMapping (column packed).
  \param a : input m-by-n matrix 
  \param s : Vector of min(m,n) singular values (descending order)
  \return : return code from lapack driver _gesvd()
 */

template <class T>
int LapackServer<T>::SVD(TArray<T>& a, TArray<T> & s)
{
  return (SVDDriver(a, s, NULL, NULL) );
}

//! Interface to Lapack SVD driver s/d/c/zgesv().
/*! Computes the vector of singular values of \b a, as well as
  right and left singular vectors of \b a.
  \f[
  A = U \Sigma V^T , ( A = U \Sigma V^H \ complex)
  \f]
  \f[
  A v_i = \sigma_i u_i \ and A^T u_i = \sigma_i v_i \ (A^H \ complex)
  \f]
  U and V are orthogonal (unitary) matrices.
  \param a : input m-by-n matrix (in FotranMemoryMapping)
  \param s : Vector of min(m,n) singular values (descending order)
  \param u : Matrix of left singular vectors
  \param vt : Transpose of right singular vectors.
  \return : return code from lapack driver _gesvd()
 */
template <class T>
int LapackServer<T>::SVD(TArray<T>& a, TArray<T> & s, TArray<T> & u, TArray<T> & vt)
{
  return (SVDDriver(a, s, &u, &vt) );
}


//! Interface to Lapack SVD driver s/d/c/zgesv().
template <class T>
int LapackServer<T>::SVDDriver(TArray<T>& a, TArray<T> & s, TArray<T>* up, TArray<T>* vtp)
{
  if ( ( a.NbDimensions() != 2 )  )
    throw(SzMismatchError("LapackServer::SVDDriver(a, ...) a.NbDimensions() != 2"));

  int_4 rowa = a.RowsKA();
  int_4 cola = a.ColsKA();

  if ( !a.IsPacked(rowa) )
     throw(SzMismatchError("LapackServer::SVDDriver(a, ...) a Not Column Packed "));

  int_4 m = a.Size(rowa);
  int_4 n = a.Size(cola);
  int_4 maxmn = (m > n) ? m : n;
  int_4 minmn = (m < n) ? m : n;

  char jobu, jobvt;
  jobu = 'N';
  jobvt = 'N';

  sa_size_t sz[2];
  if ( up != NULL) {
    if ( dynamic_cast< TVector<T> * > (vtp) ) 
      throw( TypeMismatchExc("LapackServer::SVDDriver() Wrong type (=TVector<T>) for u !") );
    up->SetMemoryMapping(BaseArray::FortranMemoryMapping);
    sz[0] = sz[1] = m;
    up->ReSize(2, sz );
    jobu = 'A';
  }
  else { 
    up = new TMatrix<T>(1,1);
    jobu = 'N';
  }
  if ( vtp != NULL) {
    if ( dynamic_cast< TVector<T> * > (vtp) ) 
      throw( TypeMismatchExc("LapackServer::SVDDriver() Wrong type (=TVector<T>) for vt !") );
    vtp->SetMemoryMapping(BaseArray::FortranMemoryMapping);
    sz[0] = sz[1] = n;
    vtp->ReSize(2, sz );
    jobvt = 'A';
  }
  else { 
    vtp = new TMatrix<T>(1,1);
    jobvt = 'N';
  }

  TVector<T> *vs = dynamic_cast< TVector<T> * > (&s);
  if (vs) vs->ReSize(minmn);
  else {
    TMatrix<T> *ms = dynamic_cast< TMatrix<T> * > (&s);
    if (ms) ms->ReSize(minmn,1);
    else  { 
      sz[0] = minmn; sz[1] = 1;
      s.ReSize(1, sz);
    }
  }
  
  int_4 lda = a.Step(a.ColsKA());
  int_4 ldu = up->Step(up->ColsKA());
  int_4 ldvt = vtp->Step(vtp->ColsKA());

  int_4 lwork = maxmn*5*wspace_size_factor;
  T * work = new T[lwork];
  int_4 info;

  if (typeid(T) == typeid(r_4) ) {
    sgesvd_(&jobu, &jobvt, &m, &n, (r_4 *)a.Data(), &lda,
            (r_4 *)s.Data(), (r_4 *) up->Data(), &ldu, (r_4 *)vtp->Data(), &ldvt, 
            (r_4 *)work, &lwork, &info);
  } else if (typeid(T) == typeid(r_8) ) {
    dgesvd_(&jobu, &jobvt, &m, &n, (r_8 *)a.Data(), &lda,
            (r_8 *)s.Data(), (r_8 *) up->Data(), &ldu, (r_8 *)vtp->Data(), &ldvt, 
            (r_8 *)work, &lwork, &info);
  } else if (typeid(T) == typeid(complex<r_4>) ) {
    r_4 * rwork = new r_4[5*minmn +5];
    r_4 * sloc  = new r_4[minmn];
    cgesvd_(&jobu, &jobvt, &m, &n, (complex<r_4> *)a.Data(), &lda,
            (r_4 *)sloc, (complex<r_4> *) up->Data(), &ldu, 
            (complex<r_4> *)vtp->Data(), &ldvt, 
            (complex<r_4> *)work, &lwork, (r_4 *)rwork, &info);
    for(int_4 i=0;i<minmn;i++) s[i] = sloc[i];
    delete [] rwork; delete [] sloc;
  } else if (typeid(T) == typeid(complex<r_8>) )  {
    r_8 * rwork = new r_8[5*minmn +5];
    r_8 * sloc  = new r_8[minmn];
    zgesvd_(&jobu, &jobvt, &m, &n, (complex<r_8> *)a.Data(), &lda,
            (r_8 *)sloc, (complex<r_8> *) up->Data(), &ldu, 
            (complex<r_8> *)vtp->Data(), &ldvt, 
            (complex<r_8> *)work, &lwork, (r_8 *)rwork, &info);
    for(int_4 i=0;i<minmn;i++) s[i] = sloc[i];
    delete [] rwork; delete [] sloc;
  } else {
    if (jobu == 'N') delete up;
    if (jobvt == 'N') delete vtp;
    string tn = typeid(T).name();
    cerr << " LapackServer::SVDDriver(...) - Unsupported DataType T = " << tn << endl;
    throw TypeMismatchExc("LapackServer::SVDDriver(a,b) - Unsupported DataType (T)");
  }

  if (jobu == 'N') delete up;
  if (jobvt == 'N') delete vtp;
  return(info);
}


//! Interface to Lapack SVD driver s/d/c/zgesdd().
/*! Same as SVD but with Divide and Conquer method */
template <class T>
int LapackServer<T>::SVD_DC(TMatrix<T>& a, TVector<T>& s, TMatrix<T>& u, TMatrix<T>& vt)
{

  if ( !a.IsPacked() )
     throw(SzMismatchError("LapackServer::SVD_DC(a, ...) a Not Packed "));

  int_4 m = a.NRows();
  int_4 n = a.NCols();
  int_4 maxmn = (m > n) ? m : n;
  int_4 minmn = (m < n) ? m : n;
  int_4 supermax = 4*minmn*minmn+4*minmn; if(maxmn>supermax) supermax=maxmn;

  char jobz = 'A';

  s.ReSize(minmn);
  u.ReSize(m,m);
  vt.ReSize(n,n);
  
  int_4 lda = n;
  int_4 ldu = u.NCols();
  int_4 ldvt = vt.NCols();
  int_4 info;  

  if(typeid(T) == typeid(r_4) ) {
    int_4 lwork = 3*minmn*minmn + supermax;
    r_4* work = new r_4[lwork +5];
    int_4* iwork = new int_4[8*minmn +5];
    sgesdd_(&jobz,&m,&n,(r_4*)a.Data(),&lda,
           (r_4*)s.Data(),(r_4*)u.Data(),&ldu,(r_4*)vt.Data(),&ldvt, 
           (r_4*)work,&lwork,(int_4*)iwork,&info);
    delete [] work; delete [] iwork;
  } else if(typeid(T) == typeid(r_8) ) {
    int_4 lwork = 3*minmn*minmn + supermax;
    r_8* work = new r_8[lwork +5];
    int_4* iwork = new int_4[8*minmn +5];
    dgesdd_(&jobz,&m,&n,(r_8*)a.Data(),&lda,
           (r_8*)s.Data(),(r_8*)u.Data(),&ldu,(r_8*)vt.Data(),&ldvt, 
           (r_8*)work,&lwork,(int_4*)iwork,&info);
    delete [] work; delete [] iwork;
  } else if(typeid(T) == typeid(complex<r_4>) ) {
    r_4* sloc = new r_4[minmn];
    int_4 lwork = minmn*minmn+2*minmn+maxmn;
    complex<r_4>* work = new complex<r_4>[lwork +5];
    r_4* rwork = new r_4[5*minmn*minmn+5*minmn +5];
    int_4* iwork = new int_4[8*minmn +5];
    cgesdd_(&jobz,&m,&n,(complex<r_4>*)a.Data(),&lda,
           (r_4*)sloc,(complex<r_4>*)u.Data(),&ldu,(complex<r_4>*)vt.Data(),&ldvt, 
           (complex<r_4>*)work,&lwork,(r_4*)rwork,(int_4*)iwork,&info);
    for(int_4 i=0;i<minmn;i++) s[i] = sloc[i];
    delete [] sloc; delete [] work; delete [] rwork; delete [] iwork;
  } else if(typeid(T) == typeid(complex<r_8>) )  {
    r_8* sloc = new r_8[minmn];
    int_4 lwork = minmn*minmn+2*minmn+maxmn;
    complex<r_8>* work = new complex<r_8>[lwork +5];
    r_8* rwork = new r_8[5*minmn*minmn+5*minmn +5];
    int_4* iwork = new int_4[8*minmn +5];
    zgesdd_(&jobz,&m,&n,(complex<r_8>*)a.Data(),&lda,
           (r_8*)sloc,(complex<r_8>*)u.Data(),&ldu,(complex<r_8>*)vt.Data(),&ldvt, 
           (complex<r_8>*)work,&lwork,(r_8*)rwork,(int_4*)iwork,&info);
    for(int_4 i=0;i<minmn;i++) s[i] = sloc[i];
    delete [] sloc; delete [] work; delete [] rwork; delete [] iwork;
  } else {
    string tn = typeid(T).name();
    cerr << " LapackServer::SVD_DC(...) - Unsupported DataType T = " << tn << endl;
    throw TypeMismatchExc("LapackServer::SVD_DC - Unsupported DataType (T)");
  }

  return(info);
}


////////////////////////////////////////////////////////////////////////////////////
/*! Computes the eigen values and eigen vectors of a symetric (or hermitian) matrix \b a.
  Input arrays should have FortranMemoryMapping (column packed).
  \param a : input symetric (or hermitian) n-by-n matrix 
  \param b : Vector of eigenvalues (descending order)
  \param eigenvector : if true compute eigenvectors, if not only eigenvalues
  \param a : on return array of eigenvectors (same order than eval, one vector = one column)
  \return : return code from lapack driver
 */

template <class T>
int LapackServer<T>::LapackEigenSym(TArray<T>& a, TVector<r_8>& b, bool eigenvector)
{
  if ( a.NbDimensions() != 2 )
    throw(SzMismatchError("LapackServer::LapackEigenSym(a,b) a NbDimensions() != 2"));
  int_4 rowa = a.RowsKA();
  int_4 cola = a.ColsKA();
  if ( a.Size(rowa) !=  a.Size(cola)) 
    throw(SzMismatchError("LapackServer::LapackEigenSym(a,b) a Not a square Array"));
  if (!a.IsPacked(rowa))
    throw(SzMismatchError("LapackServer::LapackEigenSym(a,b) a Not Column Packed"));

  char uplo='U';
  char jobz='N'; if(eigenvector) jobz='V';

  int_4 n = a.Size(rowa);
  int_4 lda = a.Step(cola);
  int_4 info = 0;

  b.ReSize(n); b = 0.;

  if (typeid(T) == typeid(r_4) ) {
    int_4 lwork = 3*n-1 +5; r_4* work = new r_4[lwork];
    r_4* w = new r_4[n];
    ssyev_(&jobz,&uplo,&n,(r_4 *)a.Data(),&lda,(r_4 *)w,(r_4 *)work,&lwork,&info);
    if(info==0) for(int i=0;i<n;i++) b(i) = w[i];
    delete [] work; delete [] w;
  } else if (typeid(T) == typeid(r_8) )  {
    int_4 lwork = 3*n-1 +5; r_8* work = new r_8[lwork];
    r_8* w = new r_8[n];
    dsyev_(&jobz,&uplo,&n,(r_8 *)a.Data(),&lda,(r_8 *)w,(r_8 *)work,&lwork,&info);
    if(info==0) for(int i=0;i<n;i++) b(i) = w[i];
    delete [] work; delete [] w;
  } else if (typeid(T) == typeid(complex<r_4>) )  {
    int_4 lwork = 2*n-1 +5; complex<r_4>* work = new complex<r_4>[lwork];
    r_4* rwork = new r_4[3*n-2  +5]; r_4* w = new r_4[n];
    cheev_(&jobz,&uplo,&n,(complex<r_4> *)a.Data(),&lda,(r_4 *)w
          ,(complex<r_4> *)work,&lwork,(r_4 *)rwork,&info);
    if(info==0) for(int i=0;i<n;i++) b(i) = w[i];
    delete [] work; delete [] rwork; delete [] w;
  } else if (typeid(T) == typeid(complex<r_8>) )  {
    int_4 lwork = 2*n-1 +5; complex<r_8>* work = new complex<r_8>[lwork];
    r_8* rwork = new r_8[3*n-2  +5]; r_8* w = new r_8[n];
    zheev_(&jobz,&uplo,&n,(complex<r_8> *)a.Data(),&lda,(r_8 *)w
          ,(complex<r_8> *)work,&lwork,(r_8 *)rwork,&info);
    if(info==0) for(int i=0;i<n;i++) b(i) = w[i];
    delete [] work; delete [] rwork; delete [] w;
  } else {
    string tn = typeid(T).name();
    cerr << " LapackServer::LapackEigenSym(a,b) - Unsupported DataType T = " << tn << endl;
    throw TypeMismatchExc("LapackServer::LapackEigenSym(a,b) - Unsupported DataType (T)");
  }

  return(info);
}

////////////////////////////////////////////////////////////////////////////////////
/*! Computes the eigen values and eigen vectors of a general squared matrix \b a.
  Input arrays should have FortranMemoryMapping (column packed).
  \param a : input general n-by-n matrix 
  \param eval : Vector of eigenvalues (complex double precision)
  \param evec : Matrix of eigenvector (same order than eval, one vector = one column)
  \param eigenvector : if true compute (right) eigenvectors, if not only eigenvalues
  \param a : on return array of eigenvectors
  \return : return code from lapack driver
  \verbatim
  eval : contains the computed eigenvalues.
         --- For real matrices "a" :
         Complex conjugate pairs of eigenvalues appear consecutively
         with the eigenvalue having the positive imaginary part first.
  evec : the right eigenvectors v(j) are stored one after another
         in the columns of evec, in the same order as their eigenvalues.
         --- For real matrices "a" :
         If the j-th eigenvalue is real, then v(j) = evec(:,j),
           the j-th column of evec.
         If the j-th and (j+1)-st eigenvalues form a complex
           conjugate pair, then v(j) = evec(:,j) + i*evec(:,j+1) and
           v(j+1) = evec(:,j) - i*evec(:,j+1).
  \endverbatim
*/

template <class T>
int LapackServer<T>::LapackEigen(TArray<T>& a, TVector< complex<r_8> >& eval, TMatrix<T>& evec, bool eigenvector)
{
  if ( a.NbDimensions() != 2 )
    throw(SzMismatchError("LapackServer::LapackEigen(a,b) a NbDimensions() != 2"));
  int_4 rowa = a.RowsKA();
  int_4 cola = a.ColsKA();
  if ( a.Size(rowa) !=  a.Size(cola)) 
    throw(SzMismatchError("LapackServer::LapackEigen(a,b) a Not a square Array"));
  if (!a.IsPacked(rowa))
    throw(SzMismatchError("LapackServer::LapackEigen(a,b) a Not Column Packed"));

  char jobvl = 'N';
  char jobvr = 'N'; if(eigenvector) jobvr='V';

  int_4 n = a.Size(rowa);
  int_4 lda = a.Step(cola);
  int_4 info = 0;

  eval.ReSize(n); eval = complex<r_8>(0.,0.);
  if(eigenvector) {evec.ReSize(n,n); evec = (T) 0.;}
  int_4 ldvr = n, ldvl = 1;

  if (typeid(T) == typeid(r_4) ) {
    int_4 lwork = 4*n +5; r_4* work = new r_4[lwork];
    r_4* wr = new r_4[n]; r_4* wi = new r_4[n]; r_4* vl = NULL;
    sgeev_(&jobvl,&jobvr,&n,(r_4 *)a.Data(),&lda,(r_4 *)wr,(r_4 *)wi,
           (r_4 *)vl,&ldvl,(r_4 *)evec.Data(),&ldvr,
           (r_4 *)work,&lwork,&info);
    if(info==0) for(int i=0;i<n;i++) eval(i) = complex<r_8>(wr[i],wi[i]);
    delete [] work; delete [] wr; delete [] wi;
  } else if (typeid(T) == typeid(r_8) )  {
    int_4 lwork = 4*n +5; r_8* work = new r_8[lwork];
    r_8* wr = new r_8[n]; r_8* wi = new r_8[n]; r_8* vl = NULL;
    dgeev_(&jobvl,&jobvr,&n,(r_8 *)a.Data(),&lda,(r_8 *)wr,(r_8 *)wi,
           (r_8 *)vl,&ldvl,(r_8 *)evec.Data(),&ldvr,
           (r_8 *)work,&lwork,&info);
    if(info==0) for(int i=0;i<n;i++) eval(i) = complex<r_8>(wr[i],wi[i]);
    delete [] work; delete [] wr; delete [] wi;
  } else if (typeid(T) == typeid(complex<r_4>) )  {
    int_4 lwork = 2*n +5; complex<r_4>* work = new complex<r_4>[lwork];
    r_4* rwork = new r_4[2*n+5]; r_4* vl = NULL; TVector< complex<r_4> > w(n);
    cgeev_(&jobvl,&jobvr,&n,(complex<r_4> *)a.Data(),&lda,(complex<r_4> *)w.Data(),
           (complex<r_4> *)vl,&ldvl,(complex<r_4> *)evec.Data(),&ldvr,
           (complex<r_4> *)work,&lwork,(r_4 *)rwork,&info);
    if(info==0) for(int i=0;i<n;i++) eval(i) = w(i);
    delete [] work; delete [] rwork;
  } else if (typeid(T) == typeid(complex<r_8>) )  {
    int_4 lwork = 2*n +5; complex<r_8>* work = new complex<r_8>[lwork];
    r_8* rwork = new r_8[2*n+5]; r_8* vl = NULL;
    zgeev_(&jobvl,&jobvr,&n,(complex<r_8> *)a.Data(),&lda,(complex<r_8> *)eval.Data(),
           (complex<r_8> *)vl,&ldvl,(complex<r_8> *)evec.Data(),&ldvr,
           (complex<r_8> *)work,&lwork,(r_8 *)rwork,&info);
    delete [] work; delete [] rwork;
  } else {
    string tn = typeid(T).name();
    cerr << " LapackServer::LapackEigen(a,b) - Unsupported DataType T = " << tn << endl;
    throw TypeMismatchExc("LapackServer::LapackEigen(a,b) - Unsupported DataType (T)");
  }

  return(info);
}





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