Index: trunk/SophyaExt/LinAlg/intflapack.cc
===================================================================
--- trunk/SophyaExt/LinAlg/intflapack.cc	(revision 2554)
+++ trunk/SophyaExt/LinAlg/intflapack.cc	(revision 2556)
@@ -4,4 +4,27 @@
 #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)
+***********************************************************/
 
 /*!
@@ -51,4 +74,5 @@
 */
 
+////////////////////////////////////////////////////////////////////////////////////
 extern "C" {
 // Le calculateur de workingspace
@@ -103,10 +127,34 @@
 	       complex<r_8>* s, complex<r_8>* u, int_4* ldu, complex<r_8>* vt, int_4* ldvt, 
 	       complex<r_8>* work, int_4* lwork, 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()
@@ -120,5 +168,5 @@
 }
 
-// --- ATTENTION BUG PREVISIBLE (CMV) --- REZA A LIRE S.T.P.
+// --- 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:
@@ -133,4 +181,5 @@
 //    ==> 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)
@@ -145,4 +194,5 @@
 }
 
+////////////////////////////////////////////////////////////////////////////////////
 //! Interface to Lapack linear system solver driver s/d/c/zgesvd().
 /*! Solve the linear system a * x = b. Input arrays
@@ -197,4 +247,5 @@
 }
 
+////////////////////////////////////////////////////////////////////////////////////
 //! Interface to Lapack linear system solver driver s/d/c/zsysvd().
 /*! Solve the linear system a * x = b with a symetric. Input arrays
@@ -210,5 +261,6 @@
 //     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 ???
+//     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 ) )
@@ -239,27 +291,27 @@
 
   if (typeid(T) == typeid(r_4) ) {
-    lwork = ilaenv_en_C(1,"SSYTRF",struplo,n,-1,-1,-1) * n;
-    work = new T[lwork+5];
+    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;
-    work = new T[lwork+5];
+    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;
-    work = new T[lwork+5];
+    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;
-    work = new T[lwork+5];
+    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 {
-    delete[] work;
+    if(work) delete[] work;
     delete[] ipiv;
     string tn = typeid(T).name();
@@ -267,9 +319,10 @@
     throw TypeMismatchExc("LapackServer::LinSolveSym(a,b) - Unsupported DataType (T)");
   }
-  delete[] work;
+  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 
@@ -354,4 +407,5 @@
 
 
+////////////////////////////////////////////////////////////////////////////////////
 //! Interface to Lapack SVD driver s/d/c/zgesv().
 /*! Computes the vector of singular values of \b a. Input arrays
@@ -489,4 +543,166 @@
 }
 
+
+////////////////////////////////////////////////////////////////////////////////////
+/*! 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 _gesvd()
+ */
+
+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 struplo[5]; struplo[0]=uplo; struplo[1]='\0';
+  char jobz='N'; if(eigenvector) jobz='V';
+  char strjobz[5]; strjobz[0]=jobz; strjobz[1]='\0';
+
+  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_(strjobz,struplo,&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_(strjobz,struplo,&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_(strjobz,struplo,&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_(strjobz,struplo,&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 _gesvd()
+  \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 strjobvl[5]; strjobvl[0] = jobvl; strjobvl[1] = '\0';
+  char jobvr = 'N'; if(eigenvector) jobvr='V';
+  char strjobvr[5]; strjobvr[0] = jobvr; strjobvr[1] = '\0';
+
+  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_(strjobvl,strjobvr,&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_(strjobvl,strjobvr,&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_(strjobvl,strjobvr,&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_(strjobvl,strjobvr,&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)
 {
