Index: /trunk/SophyaProg/Tests/tsttminv.cc
===================================================================
--- /trunk/SophyaProg/Tests/tsttminv.cc	(revision 2554)
+++ /trunk/SophyaProg/Tests/tsttminv.cc	(revision 2555)
@@ -1,12 +1,39 @@
-// Test de l'inversion de matrice (cmv 14/04/00)
-// Pb numeriques: dapax22 et daplxa115
-//                tsttminv -v 55,100,130,1000000000 -l 0 -a 1 
-
-// if defined COMPLEX , if not REAL
+// Test de l'inversion de matrices et valeurs propres (avec Lapack) (cmv 21/07/04)
+//       cmvtminv -a 1234 -l 0 -s 1 -b 25,10000 -n 50 -S
+
+///////////////////////////////////////////////////
+///////////////////////////////////////////////////
+// PARTIE POUVANT ETRE CHANGEE PAR L'UTILISATEUR //
+///////////////////////////////////////////////////
+///////////////////////////////////////////////////
+
+// --- Choix de travailler avec des matrices complexes ?
 //#define COMPLEX
-// if defined 32 bits precision, if not 64 bits
-// #define PRECIS32
-#define ALSO_LAPACK
-
+
+//////////////////////////////////////////////////
+// --- Choix de travailler en simple precision ?
+//#define PRECIS32
+
+//////////////////////////////////////////////////
+// --- Choix GausPiv + Lapack ?
+#define USE_GAUSPIV
+#define USE_LAPACK
+
+// --- Choix de ce que doit faire Lapack
+#ifdef USE_LAPACK
+#define ALSO_LAPACK_INV
+#define ALSO_LAPACK_INV_SYM
+#define ALSO_LAPACK_INV_LSS
+#define ALSO_LAPACK_EV
+#define ALSO_LAPACK_EV_SYM
+#endif
+
+//////////////////////////////////////////////////
+//////////////////////////////////////////////////
+//           NE RIEN CHANGER CI-APRES           //
+//////////////////////////////////////////////////
+//////////////////////////////////////////////////
+
+//////////////////////////////////////////////////
 #include "machdefs.h"
 #include <iostream>
@@ -21,28 +48,54 @@
 #include "array.h"
 #include "srandgen.h"
-
-#ifdef ALSO_LAPACK
+#if defined(USE_LAPACK)
 #include "intflapack.h"
 #endif
 
-
+//////////////////////////////////////////////////
 #if defined(COMPLEX)
-  #define ABS_VAL(_x_) sqrt((double)((_x_).real()*(_x_).real() + (_x_).imag()*(_x_).imag()))
   #if defined(PRECIS32)
     #define TYPE complex<r_4>
+    #define TYPER r_4
   #else
     #define TYPE complex<r_8>
+    #define TYPER r_8
   #endif
+  #define REAL_PART(_x_) (TYPE((_x_).real(),0.))
+  #define CONJ_VAL(_x_) (TYPE((_x_).real(),-(_x_).imag()))
+  #define ABS_VAL(_x_) sqrt((double)((_x_).real()*(_x_).real() + (_x_).imag()*(_x_).imag()))
 #else
-  #define ABS_VAL(_x_) fabs((double)_x_)
   #if defined(PRECIS32)
     #define TYPE r_4
+    #define TYPER r_4
   #else
     #define TYPE r_8
+    #define TYPER r_8
   #endif
-#endif
-
+  #define REAL_PART(_x_) (_x_)
+  #define CONJ_VAL(_x_) (_x_)
+  #define ABS_VAL(_x_) fabs((double)_x_)
+#endif
+
+//////////////////////////////////////////////////
+void Symetrize(TMatrix< TYPE >& A);
+void Hermitian(TMatrix< TYPE >& A);
+r_8  Check_Mat_Ident(TMatrix< TYPE >& A);
+r_8  Check_Mat_VecCol_0(TMatrix< TYPE >& A);
+void  Check_Mat_VecCol_2(TMatrix< complex<TYPER> >& A);
+#if defined(USE_LAPACK)
+/*
+-- Pour faire ce test il faut passer la methode ilaenv_en_C()
+   de LapackServer en methode "public" (dans intflapack.h)
+   et recompiler la librairie externe sophya
+*/
+// void TestIlaEnv(int_4 n);
+#endif
+
+
+//////////////////////////////////////////////////
 int main(int narg,char *arg[])
 {
+//--------------------------------------------------------
+//-- Initialisation
 //--------------------------------------------------------
 // number of lines/columns
@@ -52,7 +105,7 @@
 // number of values change by +/- vbig
 uint_4 nbig = N;
-r_8 vbig = 1000.;
+r_8 vbig = 100.;
 // Nombre de lignes de matrice a imprimer
-uint_4 nprline = 10;
+uint_4 nprline = N;
 // Initialisation du pauvre de l'aleatoire
 uint_4 nalea = 0;
@@ -60,13 +113,32 @@
 int tscal = 1;
 bool detok=false;
-bool timok = false;
-//--------------------------------------------------------
-
+// Please symetrize the input matrice
+bool symetok=false;
+// Please symetrize the input matrice
+bool gaussok=false;
+
+//--------------------------------------------------------
 //-- Decodage arguments
+//--------------------------------------------------------
 char c;
-while((c = getopt(narg,arg,"hdtv:l:a:n:")) != -1) {
+while((c = getopt(narg,arg,"Sdgn:s:b:l:a:t:h")) != -1) {
   switch (c) {
-  case 'v' :
-    sscanf(optarg,"%d,%lf,%d,%lf",&N,&scale,&nbig,&vbig);
+  case 'S' :
+    symetok = true;
+    break;
+  case 'd' :
+    detok = true;
+    break;
+  case 'g' :
+    gaussok = true;
+    break;
+  case 'n' :
+    sscanf(optarg,"%d",&N);
+    break;
+  case 's' :
+    sscanf(optarg,"%lf",&scale);
+    break;
+  case 'b' :
+    sscanf(optarg,"%d,%lf",&nbig,&vbig);
     break;
   case 'l' :
@@ -76,151 +148,408 @@
     sscanf(optarg,"%d",&nalea);
     break;
-  case 'n' :
+  case 't' :
     sscanf(optarg,"%d",&tscal);
     break;
-  case 'd' :
-    detok = true;
-    break;
-  case 't' :
-    timok = true;
-    break;
   case 'h' :
-    cout<<"tsttminv [-h] [-v N,scale,nbig,vbig] [-l nprline] [-a nalea]"<<endl
-        <<"         [-n typs] [-d] [-t]"<<endl
-        <<"matrix filled with : {[-1,1] +/- vbig(nbig time)}*scale"<<endl
-	<<"-n 0/1/2 : data scaling 0=no, 1=global (def), 2=row-by-row"<<endl
-	<<"-d : also compute determinant"<<endl
-	<<"-t : do timing"<<endl;
-    exit(-1);
+    cout<<"tsttminv [-h] [-n N] [-S] [-s scale] [-b nbig,vbig] [-g]"<<endl
+        <<"         [-l nprline] [-a nalea] [-t tscal] [-d]"<<endl
+        <<"-- matrix A(N,N) filled with {[-1,1] +/- vbig(nbig time)}*scale --"<<endl
+	<<"-g : instead of flat [-1,1] use normal gaussian distribution for A(i,j)"<<endl
+	<<"-S : symetrize the input matrix"<<endl
+	<<"-l : print nprline of input and test matrices"<<endl
+	<<"-a : for random (pseudo) changing"<<endl
+        <<"-- Only GausPiv --"<<endl
+	<<"-t 0/1/2 : data scaling 0=no, 1=global (def), 2=row-by-row"<<endl
+	<<"-d : also compute determinant"<<endl;
+    return(-1);
   }
 }
 if(N<=1) N = 1;
 cout<<"Taille matrice NxN, N = "<<N<<endl;
-cout<<"Elements entre +/- scale = "<<scale<<endl;
+if(gaussok) cout<<"Elements gaussian normal * scale = "<<scale<<endl;
+  else      cout<<"Elements entre +/- 1 * scale = "<<scale<<endl;
 cout<<"Nombre de valeurs hors standard nbig = "<<nbig<<endl;
-cout<<"Valeurs hors standard (+/- vbig = "<<vbig<<" ) * scale"<<endl;
+cout<<"Valeurs hors standard (+/- vbig = "<<vbig<<" ) * scale = "<<vbig*scale<<endl;
 cout<<"Nombre de lignes de matrice a imprimer "<<nprline<<endl;
 cout<<"Initialisation de l aleatoire par "<<nalea<<" tirages"<<endl;
- cout<<"Data scaling "<<tscal<<" determinant="<<detok<<" timing="<<timok<<endl;
+cout<<"Data scaling "<<tscal<<" determinant="<<detok<<endl;
+if(symetok) cout<<"Input matrix has been symetrized "<<symetok<<endl;
 cout<<endl;
 
+//--------------------------------------------------------
+// TestIlaEnv(N); return -41;
+//--------------------------------------------------------
+
+//--------------------------------------------------------
 //-- Initialization arrays
+//--------------------------------------------------------
 SophyaInit();
-if(timok) InitTim();
-#ifdef ALSO_LAPACK
+InitTim();
+#if defined(USE_LAPACK)
   BaseArray::SetDefaultMemoryMapping(BaseArray::FortranMemoryMapping);
 #endif
-
-//-- Definition arrays
-TMatrix< TYPE > A(N,N);
-TMatrix< TYPE > InvA(N,N);
-TMatrix< TYPE > AiA(N,N);
-TMatrix< TYPE > B(N,N);
-TMatrix< TYPE > Adum(N,N);
-TVector< int >  Vind(N*N);
-A.Show();
-
-SimpleMatrixOperation< TYPE >::SetGausPivScal(tscal);
-
-//-- Mise a zero
-A    = (TYPE) 0;
-InvA = (TYPE) 0;
-AiA  = (TYPE) 0;
-B    = (TYPE) 0;
-Adum = (TYPE) 0;
+if(nalea>0) for(int i=0;i<nalea;i++) drand01();
 BaseArray::SetMaxPrint(nprline*N,0);
 
-//-- Fill matrices
-double vmax=-1.;
-uint_8 k;
-uint_4 i,j; double s;
-uint_4 nbr=0, nbi=0;
-Vind = 0;
-if(nalea>0) for(i=0;i<nalea;i++) drand01();
-A = RandomSequence(RandomSequence::Flat,0.,1.);
-if(nbig>0) for(i=0;i<nbig;i++) {
-  k=(uint_8)(drand01()*N*N); if(k>=N*N) k--;
-  s=(drand01()>0.5)?1.:-1.;
-  if(Vind[k]==0) {A[k] += (TYPE) s*vbig; Vind[k]=1; nbr++;}
-}
-A *= (TYPE) scale;
+//--------------------------------------------------------
+//-- Definition global arrays
+//--------------------------------------------------------
+TMatrix< TYPE > Ainput(N,N); Ainput = (TYPE) 0;
+TMatrix< TYPE > A(N,N); A = (TYPE) 0;
+Ainput.Show();
+
+//--------------------------------------------------------
+//-- Fill matrices with flat random
+//--------------------------------------------------------
+if(gaussok) Ainput = RandomSequence(RandomSequence::Gaussian,0.,1.);
+  else      Ainput = RandomSequence(RandomSequence::Flat,0.,1.);
 #if defined(COMPLEX)
-Vind = 0;
-B = RandomSequence(RandomSequence::Flat,0.,1.);
-if(nbig>0) for(i=0;i<nbig;i++) {
-  k=(uint_8)(drand01()*N*N); if(k>=N*N) k--;
-  s=(drand01()>0.5)?1.:-1.;
-  if(Vind[k]==0) {B[k] += (TYPE) s*vbig; Vind[k]=1; nbi++;}
-}
-B *= (TYPE) scale;
-A += TYPE(0.,1.)*B;
-#endif
-cout<<"Nombre de valeurs BIG reelles = "<<nbr<<", imaginaires = "<<nbi<<endl;
-Adum = A;
-
-//-- Print matrice A
-cout<<"------------ TMatrix A :"<<endl;
-if(nprline>0) {cout<<A; cout<<endl<<endl;}
-
-//-- Inversion
-cout<<"------------ Inversion"<<endl;
-//InvA = Inverse(A);
-InvA = IdentityMatrix(1.,N);
-if(timok) PrtTim("--- End of filling ---");
-TYPE det = GausPiv(A,InvA,detok);
-if(timok) PrtTim("--- End of GausPiv ---");
-cout<<"Det = "<<det<<endl;
-cout<<"------------ TMatrix InvA = A^(-1):"<<endl;
-if(nprline>0) {cout<<InvA; cout<<endl<<endl;}
-
-//-- AiA = A * InvA
-cout<<"AiA = A * InvA"<<endl;
-AiA = A * InvA;
-cout<<"------------ TMatrix AiA = A * InvA:"<<endl;
-if(nprline>0) {cout<<AiA; cout<<endl;}
-
-//-- Check
-vmax=-1.;
-for(i=0;i<N;i++) {
-  double absv = ABS_VAL( 1. - AiA(i,i) );
-  if( absv > vmax ) vmax = absv;
-}
-cout<<"Ecart maximum par rapport a 1 sur diagonale = "<<vmax<<endl;
-vmax = -1.;
-for(i=0;i<N;i++) for(j=0;j<N;j++) {
-  if( i == j ) continue;
-  double absv = ABS_VAL( AiA(i,j) );
-  if( absv > vmax ) vmax = absv;
-}
-cout<<"Ecart maximum par rapport a 0 hors diagonale = "<<vmax<<endl;
-
+if(gaussok) A = RandomSequence(RandomSequence::Gaussian,0.,1.);
+  else      A = RandomSequence(RandomSequence::Flat,0.,1.);
+Ainput += TYPE(0.,1.)*A;
+#endif
+
+//--------------------------------------------------------
+//-- Fill matrices with big values
+//--------------------------------------------------------
+if(nbig>0) {
+#if defined(COMPLEX)
+  nbig = (nbig+1)/2;
+#endif
+  TMatrix< uint_2 > Vind(N,N); Vind = 0;
+  // for real part
+  uint_4 nbr=0;
+  for(int k=0;k<nbig;k++) {
+    int i = (int) (drand01()*N); int j = (int) (drand01()*N);
+    double s=(drand01()>0.5)?1.:-1.;
+    if(Vind(i,j)==0) {Ainput(i,j) += (TYPER) s*vbig; Vind(i,j)+=1; nbr++;}
+  }
+  cout<<"Nombre de valeurs BIG reelles = "<<nbr<<endl;
+#if defined(COMPLEX)
+  // for imaginary part
+  uint_4 nbi=0;
+  for(int k=0;k<nbig;k++) {
+    int i = (int) (drand01()*N); int j = (int) (drand01()*N);
+    double s=(drand01()>0.5)?1.:-1.;
+    if(Vind(i,j)<=1) {Ainput(i,j) += TYPE(0.,(TYPER)s*vbig); Vind(i,j)+=2; nbi++;}
+  }
+  cout<<"Nombre de valeurs BIG imaginaires = "<<nbi<<endl;
+  cout<<"Nombre de valeurs BIG = "<<nbr+nbi<<endl;
+#endif
+}
+
+//--------------------------------------------------------
+//-- Scale matrix for machine precision tests
+//--------------------------------------------------------
+Ainput *= (TYPE) scale;
+
+//--------------------------------------------------------
+//-- Create symetric matrix for all A if requested
+//--------------------------------------------------------
+if(symetok) Symetrize(Ainput);
+
+//--------------------------------------------------------
+//-- Print matrice Ainput
+//--------------------------------------------------------
+cout<<"------------ TMatrix Ainput :"<<endl;
+if(nprline>0) {cout<<Ainput; cout<<endl;}
+PrtTim("--- End of Matrix filling ---");
+
+
+#ifdef ALSO_LAPACK_INV
 ////////////////////////////////////
 ///////// Test avec Lapack /////////
 ////////////////////////////////////
-#ifdef ALSO_LAPACK
-
-cout<<endl<<"LAPACK Cross-Check"<<endl;
-InvA = IdentityMatrix(1.,N);
-if(timok) PrtTim("--- End of check ---");
-LapackLinSolve(Adum,InvA);
-if(timok) PrtTim("--- End of LapackLinSolve ---");
-AiA = A * InvA;
-vmax=-1.;
-for(i=0;i<N;i++) {
-  double absv = ABS_VAL( 1. - AiA(i,i) );
-  if( absv > vmax ) vmax = absv;
-}
-cout<<"Ecart maximum par rapport a 1 sur diagonale = "<<vmax<<endl;
-vmax = -1.;
-for(i=0;i<N;i++) for(j=0;j<N;j++) {
-  if( i == j ) continue;
-  double absv = ABS_VAL( AiA(i,j) );
-  if( absv > vmax ) vmax = absv;
-}
-cout<<"Ecart maximum par rapport a 0 hors diagonale = "<<vmax<<endl;
-
-#endif
-
-if(timok) PrtTim("--- End of Job ---");
+{
+cout<<"\n=========================================="<<endl;
+cout<<"------------ Inversion LAPACK"<<endl;
+A = Ainput;
+//-- Inversion
+TMatrix< TYPE > InvA(N,N); InvA = IdentityMatrix(1.,N);
+int_4 info = LapackLinSolve(A,InvA);
+cout<<"info="<<info<<endl;
+PrtTim("--- End of LapackLinSolve Inversion ---");
+//-- AiA = A * InvA
+cout<<"Compute AiA = A * InvA"<<endl;
+TMatrix< TYPE > AiA(N,N); AiA = Ainput * InvA;
+cout<<"------------ TMatrix AiA = A * InvA:"<<endl;
+if(nprline>0) {cout<<AiA; cout<<endl;}
+//-- Check
+Check_Mat_Ident(AiA);
+PrtTim("--- End of LapackLinSolve Test ---");
+}
+#endif
+
+
+#ifdef ALSO_LAPACK_INV_SYM
+////////////////////////////////////////
+///////// Test avec Lapack sym /////////
+////////////////////////////////////////
+{
+cout<<"\n=========================================="<<endl;
+cout<<"------------ Inversion LAPACK sym"<<endl;
+TMatrix< TYPE > Asym(N,N); Asym=Ainput; Symetrize(Asym); A=Asym;
+//-- Inversion
+TMatrix< TYPE > InvA(N,N); InvA = IdentityMatrix(1.,N);
+int_4 info = LapackLinSolveSym(A,InvA);
+cout<<"info="<<info<<endl;
+PrtTim("--- End of LapackLinSolveSym Inversion ---");
+//-- AiA = A * InvA
+cout<<"Compute AiA = A * InvA"<<endl;
+TMatrix< TYPE > AiA(N,N); AiA = Asym * InvA;
+cout<<"------------ TMatrix AiA = A * InvA:"<<endl;
+if(nprline>0) {cout<<AiA; cout<<endl;} 
+//-- Check
+Check_Mat_Ident(AiA);
+PrtTim("--- End of LapackLinSolveSym Test ---");
+}
+#endif
+
+
+#ifdef ALSO_LAPACK_INV_LSS
+////////////////////////////////////////////////
+///////// Test avec Lapack LeastSquare /////////
+////////////////////////////////////////////////
+{
+cout<<"\n=========================================="<<endl;
+cout<<"------------ Inversion LAPACK LeastSquare"<<endl;
+A = Ainput;
+//-- Inversion
+TMatrix< TYPE > InvA(N,N); InvA = IdentityMatrix(1.,N);
+int_4 info = LapackLeastSquareSolve(A,InvA);
+cout<<"info="<<info<<endl;
+PrtTim("--- End of LapackLeastSquareSolve Inversion ---");
+//-- AiA = A * InvA
+cout<<"Compute AiA = A * InvA"<<endl;
+TMatrix< TYPE > AiA(N,N); AiA = Ainput * InvA;
+cout<<"------------ TMatrix AiA = A * InvA:"<<endl;
+if(nprline>0) {cout<<AiA; cout<<endl;} 
+//-- Check
+Check_Mat_Ident(AiA);
+PrtTim("--- End of LapackLeastSquareSolve Test ---");
+}
+#endif
+
+
+#ifdef ALSO_LAPACK_EV
+////////////////////////////////////////////////
+///////// Test avec Lapack pour EV /////////
+////////////////////////////////////////////////
+{
+cout<<"\n=========================================="<<endl;
+cout<<"------------ Eigen decompositon LapackEigen "<<endl;
+A=Ainput;
+TMatrix< TYPE > Evec(N,N); Evec = (TYPE) 0;
+TVector< complex<r_8> > Eval(N); Eval = complex<r_8>(0,0);
+//-- Decompositon
+int_4 info = LapackEigen(A,Eval,Evec,true);
+cout<<"info="<<info<<endl;
+PrtTim("--- End of LapackEigen decompositon ---");
+if(nprline>0) {cout<<Eval; cout<<endl; cout<<Evec; cout<<endl;}
+#ifndef COMPLEX
+//-- Find the complex conjugate pairs
+TVector< uint_2 > Evalconj(N); Evalconj = 0;
+int_4 nconj=0;
+for(int i=0;i<N-1;i++) {
+  if(Evalconj(i)!=0) continue;   // deja traite
+  if(Eval(i).imag()==0.) continue; // real eigenvalue
+  if(fabs(Eval(i).imag()+Eval(i+1).imag())>1e-150) continue; // les 2 consecutives ne sont pas conjuguees
+  if(Eval(i).imag()<0.) continue; // first conjugate have positive imaginary part
+  if(Eval(i+1).imag()>0.) continue;
+  Evalconj(i) = 1; Evalconj(i+1) = 2; nconj++;
+}
+//cout<<Evalconj<<endl;
+cout<<"Number of conjugate eigen values: "<<nconj<<" *2 = "<<2*nconj<<" / "<<N<<endl;
+#endif
+//-- Azmlz = A*z(l) - l*z(l)
+cout<<"Compute Azmlz(l) = A*z(l) - l*z(l)"<<endl;
+TMatrix< complex<TYPER> > Azmlz(N,N); Azmlz = (complex<TYPER>) 0;
+for(int l=0;l<N;l++) { // eigen value
+  complex<TYPER> Eval_l = complex<TYPER>(Eval(l).real(),Eval(l).imag());
+  for(int i=0;i<N;i++) {
+    complex<TYPER> Evec_il;
+    #ifdef COMPLEX
+    Evec_il = Evec(i,l);
+    #else
+    Evec_il = complex<TYPER>(Evec(i,l),0.);
+    if(Evalconj(l)==1)        Evec_il = complex<TYPER>(Evec(i,l),Evec(i,l+1));
+      else if(Evalconj(l)==2) Evec_il = complex<TYPER>(Evec(i,l-1),-Evec(i,l));
+    #endif
+    for(int j=0;j<N;j++) {
+      complex<TYPER> Evec_jl;
+      #ifdef COMPLEX
+      Evec_jl = Evec(j,l);
+      #else
+      Evec_jl = complex<TYPER>(Evec(j,l),0.);
+      if(Evalconj(l)==1)        Evec_jl = complex<TYPER>(Evec(j,l),Evec(j,l+1));
+        else if(Evalconj(l)==2) Evec_jl = complex<TYPER>(Evec(j,l-1),-Evec(j,l));
+      #endif
+      Azmlz(i,l) += Ainput(i,j) * Evec_jl;
+    }
+    Azmlz(i,l) -= Eval_l*Evec_il;
+  }
+}
+if(nprline>0) {cout<<Azmlz; cout<<endl;} 
+//-- Check
+Check_Mat_VecCol_2(Azmlz);
+PrtTim("--- End of LapackEigen Test ---");
+}
+#endif
+
+
+#ifdef ALSO_LAPACK_EV_SYM
+////////////////////////////////////////////////
+///////// Test avec Lapack sym pour EV /////////
+////////////////////////////////////////////////
+{
+cout<<"\n=========================================="<<endl;
+cout<<"------------ Eigen decompositon LapackEigenSym "<<endl;
+TMatrix< TYPE > Aher(N,N); Aher=Ainput; Hermitian(Aher); A=Aher;
+TVector<r_8> Eval;
+//-- Decompositon
+int_4 info = LapackEigenSym(A,Eval,true);
+cout<<"info="<<info<<endl;
+PrtTim("--- End of LapackEigenSym decompositon ---");
+if(nprline>0) {cout<<Eval; cout<<endl; cout<<A; cout<<endl;} 
+//-- Azmlz = A*z(l) - l*z(l)
+// le vecteur propre z pour la l-ieme valeur propre est dans A(.,l):
+//    z_i = A(i,l)  ou "l" est la l-ieme valeur propre
+cout<<"Compute Azmlz(l) = A*z(l) - l*z(l)"<<endl;
+TMatrix< TYPE > Azmlz(N,N); Azmlz = (TYPE) 0;
+for(int l=0;l<N;l++)  // eigen value
+  for(int i=0;i<N;i++)
+    {for(int j=0;j<N;j++) Azmlz(i,l)+=Aher(i,j)*A(j,l); Azmlz(i,l)-=(TYPER)Eval(l)*A(i,l);}
+if(nprline>0) {cout<<Azmlz; cout<<endl;} 
+//-- Check
+Check_Mat_VecCol_0(Azmlz);
+PrtTim("--- End of LapackEigenSym Test ---");
+}
+#endif
+
+
+#ifdef USE_GAUSPIV
+////////////////////////////////////
+///////// Test avec GausPiv /////////
+////////////////////////////////////
+{
+cout<<"\n==========================================\n"
+    <<"------------ Inversion GausPiv"<<endl;
+SimpleMatrixOperation< TYPE >::SetGausPivScal(tscal);
+A = Ainput;
+//-- Inversion
+TMatrix< TYPE > InvA(N,N); InvA = IdentityMatrix(1.,N);
+TYPE det = GausPiv(A,InvA,detok);
+PrtTim("--- End of GausPiv Inversion ---");
+cout<<"Det = "<<det<<endl;
+cout<<"------------ TMatrix InvA = A^(-1):"<<endl;
+//-- AiA = A * InvA
+cout<<"Compute AiA = A * InvA"<<endl;
+TMatrix< TYPE > AiA(N,N); AiA = Ainput * InvA;
+cout<<"------------ TMatrix AiA = A * InvA:"<<endl;
+if(nprline>0) {cout<<AiA; cout<<endl;}
+//-- Check
+Check_Mat_Ident(AiA);
+PrtTim("--- End of GausPiv Test ---");
+}
+#endif
+
+
+PrtTim("--- End of Job ---");
 exit(0);
 }
+
+
+
+////////////////////////////////////////////////////////////
+////-------------------------------------------------------
+void Symetrize(TMatrix< TYPE >& A)
+// Symetrize A
+{
+ int_4 N = A.NRows();
+ for(int i=0;i<N-1;i++) for(int j=i+1;j<N;j++) A(j,i) = A(i,j);
+}
+
+////-------------------------------------------------------
+void Hermitian(TMatrix< TYPE >& A)
+// Put A hermitian
+{
+ int_4 N = A.NRows();
+ for(int i=0;i<N-1;i++) for(int j=i+1;j<N;j++) A(j,i) = CONJ_VAL(A(i,j));
+ for(int i=0;i<N;i++) A(i,i) = REAL_PART(A(i,i));
+}
+
+////-------------------------------------------------------
+r_8  Check_Mat_Ident(TMatrix< TYPE >& A)
+// Compute the biggest difference element by element of A / Identity
+{
+ int_4 N = A.NRows();
+ r_8 vmaxd=-1.;
+ for(int i=0;i<N;i++)
+   if( ABS_VAL((TYPER)1.-A(i,i)) > vmaxd ) vmaxd = ABS_VAL((TYPER)1.-A(i,i));
+ cout<<"Ecart maximum par rapport a 1 sur diagonale = "<<vmaxd<<endl;
+ r_8 vmaxh = -1.;
+ for(int i=0;i<N;i++) for(int j=0;j<N;j++) {
+   if(i==j) continue;
+   if( ABS_VAL(A(i,j)) > vmaxh ) vmaxh = ABS_VAL(A(i,j));
+ }
+ cout<<"Ecart maximum par rapport a 0 hors diagonale = "<<vmaxh<<endl;
+ return (vmaxd>vmaxh)? vmaxd: vmaxh;
+}
+
+////-------------------------------------------------------
+r_8  Check_Mat_VecCol_0(TMatrix< TYPE >& A)
+// Return the biggest norm of the N vectors column of matrix
+{
+ int_4 N = A.NRows();
+ r_8 vmax=-1.;
+ for(int l=0;l<N;l++) {
+   r_8 absv = 0.;
+   for(int i=0;i<N;i++) absv += ABS_VAL(A(i,l)) * ABS_VAL(A(i,l));
+   if( absv > vmax ) vmax = absv;
+ }
+ vmax = sqrt(vmax);
+ cout<<"Longueur max de ||A*z-l*z|| pour tous l = "<<vmax<<endl;
+ return vmax;
+}
+
+////-------------------------------------------------------
+void  Check_Mat_VecCol_2(TMatrix< complex<TYPER> >& A)
+// Return the biggest norm of :
+//   - the real part of the N vectors column of matrix
+//   - the imaginary part of the N vectors column of matrix
+//   - the module of the N vectors column of matrix
+{
+ int_4 N = A.NRows();
+ r_8 vmaxr=-1., vmaxi=-1., vmaxn=-1.;
+ for(int l=0;l<N;l++) {
+   double absvr = 0., absvi = 0., absvn = 0.;
+   for(int i=0;i<N;i++) {
+     absvr += A(i,l).real()*A(i,l).real();
+     absvi += A(i,l).imag()*A(i,l).imag();
+     absvn += A(i,l).real()*A(i,l).real() + A(i,l).imag()*A(i,l).imag();
+   }
+   if( absvr > vmaxr ) vmaxr = absvr;
+   if( absvi > vmaxi ) vmaxi = absvi;
+   if( absvn > vmaxn ) vmaxn = absvn;
+ }
+ vmaxr=sqrt(vmaxr); vmaxi=sqrt(vmaxi); vmaxn=sqrt(vmaxn);
+ cout<<"Longueur max de ||A*z-l*z|| pour tous l, reel = "<<vmaxr
+     <<", imag = "<<vmaxi<<", module = "<<vmaxn<<endl;
+}
+
+
+/*
+void TestIlaEnv(int_4 n)
+{
+LapackServer<TYPE> lps;
+cout<<"TestIlaEnv n="<<n<<endl;
+cout<<lps.ilaenv_en_C(1,"SSYTRF","U",n,-1,-1,-1)<<endl;
+cout<<lps.ilaenv_en_C(1,"DSYTRF","U",n,-1,-1,-1)<<endl;
+cout<<lps.ilaenv_en_C(1,"CSYTRF","U",n,-1,-1,-1)<<endl;
+cout<<lps.ilaenv_en_C(1,"ZSYTRF","U",n,-1,-1,-1)<<endl;
+cout<<lps.ilaenv_en_C(1,"SSYTRF","L",n,-1,-1,-1)<<endl;
+cout<<lps.ilaenv_en_C(1,"DSYTRF","L",n,-1,-1,-1)<<endl;
+cout<<lps.ilaenv_en_C(1,"CSYTRF","L",n,-1,-1,-1)<<endl;
+cout<<lps.ilaenv_en_C(1,"ZSYTRF","L",n,-1,-1,-1)<<endl;
+}
+*/
