source: Sophya/trunk/SophyaProg/Tests/tsttminv.cc@ 2555

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

//////////////////////////////////////////////////
// --- 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>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <unistd.h>
#include "timing.h"
#include "ntoolsinit.h"
#include "pexceptions.h"
#include "array.h"
#include "srandgen.h"
#if defined(USE_LAPACK)
#include "intflapack.h"
#endif

//////////////////////////////////////////////////
#if defined(COMPLEX)
  #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
  #if defined(PRECIS32)
    #define TYPE r_4
    #define TYPER r_4
  #else
    #define TYPE r_8
    #define TYPER r_8
  #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
uint_4 N = 5;
// scale of the value (if =1 values between -1 and 1)
r_8 scale = 1.;
// number of values change by +/- vbig
uint_4 nbig = N;
r_8 vbig = 100.;
// Nombre de lignes de matrice a imprimer
uint_4 nprline = N;
// Initialisation du pauvre de l'aleatoire
uint_4 nalea = 0;
// data scaling for gauss pivoting and determinant
int tscal = 1;
bool detok=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,"Sdgn:s:b:l:a:t:h")) != -1) {
  switch (c) {
  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' :
    sscanf(optarg,"%d",&nprline);
    break;
  case 'a' :
    sscanf(optarg,"%d",&nalea);
    break;
  case 't' :
    sscanf(optarg,"%d",&tscal);
    break;
  case 'h' :
    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;
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 = "<<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<<endl;
if(symetok) cout<<"Input matrix has been symetrized "<<symetok<<endl;
cout<<endl;

//--------------------------------------------------------
// TestIlaEnv(N); return -41;
//--------------------------------------------------------

//--------------------------------------------------------
//-- Initialization arrays
//--------------------------------------------------------
SophyaInit();
InitTim();
#if defined(USE_LAPACK)
  BaseArray::SetDefaultMemoryMapping(BaseArray::FortranMemoryMapping);
#endif
if(nalea>0) for(int i=0;i<nalea;i++) drand01();
BaseArray::SetMaxPrint(nprline*N,0);

//--------------------------------------------------------
//-- 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)
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 /////////
////////////////////////////////////
{
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;
}
*/