Changeset 2555 in Sophya


Ignore:
Timestamp:
Jul 21, 2004, 5:49:10 PM (21 years ago)
Author:
cmv
Message:
  • Nouveau tests inversion matrices (cas des symetriques avec Lapack)
  • Nouveau tests decompositions eigenvalues (gene+sym+hermit)
  • refonte du programme pour + de convivialite ! (cmv 21/07/04)
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaProg/Tests/tsttminv.cc

    r2322 r2555  
    1 // Test de l'inversion de matrice (cmv 14/04/00)
    2 // Pb numeriques: dapax22 et daplxa115
    3 //                tsttminv -v 55,100,130,1000000000 -l 0 -a 1
    4 
    5 // if defined COMPLEX , if not REAL
     1// Test de l'inversion de matrices et valeurs propres (avec Lapack) (cmv 21/07/04)
     2//       cmvtminv -a 1234 -l 0 -s 1 -b 25,10000 -n 50 -S
     3
     4///////////////////////////////////////////////////
     5///////////////////////////////////////////////////
     6// PARTIE POUVANT ETRE CHANGEE PAR L'UTILISATEUR //
     7///////////////////////////////////////////////////
     8///////////////////////////////////////////////////
     9
     10// --- Choix de travailler avec des matrices complexes ?
    611//#define COMPLEX
    7 // if defined 32 bits precision, if not 64 bits
    8 // #define PRECIS32
    9 #define ALSO_LAPACK
    10 
     12
     13//////////////////////////////////////////////////
     14// --- Choix de travailler en simple precision ?
     15//#define PRECIS32
     16
     17//////////////////////////////////////////////////
     18// --- Choix GausPiv + Lapack ?
     19#define USE_GAUSPIV
     20#define USE_LAPACK
     21
     22// --- Choix de ce que doit faire Lapack
     23#ifdef USE_LAPACK
     24#define ALSO_LAPACK_INV
     25#define ALSO_LAPACK_INV_SYM
     26#define ALSO_LAPACK_INV_LSS
     27#define ALSO_LAPACK_EV
     28#define ALSO_LAPACK_EV_SYM
     29#endif
     30
     31//////////////////////////////////////////////////
     32//////////////////////////////////////////////////
     33//           NE RIEN CHANGER CI-APRES           //
     34//////////////////////////////////////////////////
     35//////////////////////////////////////////////////
     36
     37//////////////////////////////////////////////////
    1138#include "machdefs.h"
    1239#include <iostream>
    … …  
    2148#include "array.h"
    2249#include "srandgen.h"
    23 
    24 #ifdef ALSO_LAPACK
     50#if defined(USE_LAPACK)
    2551#include "intflapack.h"
    2652#endif
    2753
    28 
     54//////////////////////////////////////////////////
    2955#if defined(COMPLEX)
    30   #define ABS_VAL(_x_) sqrt((double)((_x_).real()*(_x_).real() + (_x_).imag()*(_x_).imag()))
    3156  #if defined(PRECIS32)
    3257    #define TYPE complex<r_4>
     58    #define TYPER r_4
    3359  #else
    3460    #define TYPE complex<r_8>
     61    #define TYPER r_8
    3562  #endif
     63  #define REAL_PART(_x_) (TYPE((_x_).real(),0.))
     64  #define CONJ_VAL(_x_) (TYPE((_x_).real(),-(_x_).imag()))
     65  #define ABS_VAL(_x_) sqrt((double)((_x_).real()*(_x_).real() + (_x_).imag()*(_x_).imag()))
    3666#else
    37   #define ABS_VAL(_x_) fabs((double)_x_)
    3867  #if defined(PRECIS32)
    3968    #define TYPE r_4
     69    #define TYPER r_4
    4070  #else
    4171    #define TYPE r_8
     72    #define TYPER r_8
    4273  #endif
    43 #endif
    44 
     74  #define REAL_PART(_x_) (_x_)
     75  #define CONJ_VAL(_x_) (_x_)
     76  #define ABS_VAL(_x_) fabs((double)_x_)
     77#endif
     78
     79//////////////////////////////////////////////////
     80void Symetrize(TMatrix< TYPE >& A);
     81void Hermitian(TMatrix< TYPE >& A);
     82r_8  Check_Mat_Ident(TMatrix< TYPE >& A);
     83r_8  Check_Mat_VecCol_0(TMatrix< TYPE >& A);
     84void  Check_Mat_VecCol_2(TMatrix< complex<TYPER> >& A);
     85#if defined(USE_LAPACK)
     86/*
     87-- Pour faire ce test il faut passer la methode ilaenv_en_C()
     88   de LapackServer en methode "public" (dans intflapack.h)
     89   et recompiler la librairie externe sophya
     90*/
     91// void TestIlaEnv(int_4 n);
     92#endif
     93
     94
     95//////////////////////////////////////////////////
    4596int main(int narg,char *arg[])
    4697{
     98//--------------------------------------------------------
     99//-- Initialisation
    47100//--------------------------------------------------------
    48101// number of lines/columns
    … …  
    52105// number of values change by +/- vbig
    53106uint_4 nbig = N;
    54 r_8 vbig = 1000.;
     107r_8 vbig = 100.;
    55108// Nombre de lignes de matrice a imprimer
    56 uint_4 nprline = 10;
     109uint_4 nprline = N;
    57110// Initialisation du pauvre de l'aleatoire
    58111uint_4 nalea = 0;
    … …  
    60113int tscal = 1;
    61114bool detok=false;
    62 bool timok = false;
    63 //--------------------------------------------------------
    64 
     115// Please symetrize the input matrice
     116bool symetok=false;
     117// Please symetrize the input matrice
     118bool gaussok=false;
     119
     120//--------------------------------------------------------
    65121//-- Decodage arguments
     122//--------------------------------------------------------
    66123char c;
    67 while((c = getopt(narg,arg,"hdtv:l:a:n:")) != -1) {
     124while((c = getopt(narg,arg,"Sdgn:s:b:l:a:t:h")) != -1) {
    68125  switch (c) {
    69   case 'v' :
    70     sscanf(optarg,"%d,%lf,%d,%lf",&N,&scale,&nbig,&vbig);
     126  case 'S' :
     127    symetok = true;
     128    break;
     129  case 'd' :
     130    detok = true;
     131    break;
     132  case 'g' :
     133    gaussok = true;
     134    break;
     135  case 'n' :
     136    sscanf(optarg,"%d",&N);
     137    break;
     138  case 's' :
     139    sscanf(optarg,"%lf",&scale);
     140    break;
     141  case 'b' :
     142    sscanf(optarg,"%d,%lf",&nbig,&vbig);
    71143    break;
    72144  case 'l' :
    … …  
    76148    sscanf(optarg,"%d",&nalea);
    77149    break;
    78   case 'n' :
     150  case 't' :
    79151    sscanf(optarg,"%d",&tscal);
    80152    break;
    81   case 'd' :
    82     detok = true;
    83     break;
    84   case 't' :
    85     timok = true;
    86     break;
    87153  case 'h' :
    88     cout<<"tsttminv [-h] [-v N,scale,nbig,vbig] [-l nprline] [-a nalea]"<<endl
    89         <<"         [-n typs] [-d] [-t]"<<endl
    90         <<"matrix filled with : {[-1,1] +/- vbig(nbig time)}*scale"<<endl
    91         <<"-n 0/1/2 : data scaling 0=no, 1=global (def), 2=row-by-row"<<endl
    92         <<"-d : also compute determinant"<<endl
    93         <<"-t : do timing"<<endl;
    94     exit(-1);
     154    cout<<"tsttminv [-h] [-n N] [-S] [-s scale] [-b nbig,vbig] [-g]"<<endl
     155        <<"         [-l nprline] [-a nalea] [-t tscal] [-d]"<<endl
     156        <<"-- matrix A(N,N) filled with {[-1,1] +/- vbig(nbig time)}*scale --"<<endl
     157        <<"-g : instead of flat [-1,1] use normal gaussian distribution for A(i,j)"<<endl
     158        <<"-S : symetrize the input matrix"<<endl
     159        <<"-l : print nprline of input and test matrices"<<endl
     160        <<"-a : for random (pseudo) changing"<<endl
     161        <<"-- Only GausPiv --"<<endl
     162        <<"-t 0/1/2 : data scaling 0=no, 1=global (def), 2=row-by-row"<<endl
     163        <<"-d : also compute determinant"<<endl;
     164    return(-1);
    95165  }
    96166}
    97167if(N<=1) N = 1;
    98168cout<<"Taille matrice NxN, N = "<<N<<endl;
    99 cout<<"Elements entre +/- scale = "<<scale<<endl;
     169if(gaussok) cout<<"Elements gaussian normal * scale = "<<scale<<endl;
     170  else      cout<<"Elements entre +/- 1 * scale = "<<scale<<endl;
    100171cout<<"Nombre de valeurs hors standard nbig = "<<nbig<<endl;
    101 cout<<"Valeurs hors standard (+/- vbig = "<<vbig<<" ) * scale"<<endl;
     172cout<<"Valeurs hors standard (+/- vbig = "<<vbig<<" ) * scale = "<<vbig*scale<<endl;
    102173cout<<"Nombre de lignes de matrice a imprimer "<<nprline<<endl;
    103174cout<<"Initialisation de l aleatoire par "<<nalea<<" tirages"<<endl;
    104  cout<<"Data scaling "<<tscal<<" determinant="<<detok<<" timing="<<timok<<endl;
     175cout<<"Data scaling "<<tscal<<" determinant="<<detok<<endl;
     176if(symetok) cout<<"Input matrix has been symetrized "<<symetok<<endl;
    105177cout<<endl;
    106178
     179//--------------------------------------------------------
     180// TestIlaEnv(N); return -41;
     181//--------------------------------------------------------
     182
     183//--------------------------------------------------------
    107184//-- Initialization arrays
     185//--------------------------------------------------------
    108186SophyaInit();
    109 if(timok) InitTim();
    110 #ifdef ALSO_LAPACK
     187InitTim();
     188#if defined(USE_LAPACK)
    111189  BaseArray::SetDefaultMemoryMapping(BaseArray::FortranMemoryMapping);
    112190#endif
    113 
    114 //-- Definition arrays
    115 TMatrix< TYPE > A(N,N);
    116 TMatrix< TYPE > InvA(N,N);
    117 TMatrix< TYPE > AiA(N,N);
    118 TMatrix< TYPE > B(N,N);
    119 TMatrix< TYPE > Adum(N,N);
    120 TVector< int >  Vind(N*N);
    121 A.Show();
    122 
    123 SimpleMatrixOperation< TYPE >::SetGausPivScal(tscal);
    124 
    125 //-- Mise a zero
    126 A    = (TYPE) 0;
    127 InvA = (TYPE) 0;
    128 AiA  = (TYPE) 0;
    129 B    = (TYPE) 0;
    130 Adum = (TYPE) 0;
     191if(nalea>0) for(int i=0;i<nalea;i++) drand01();
    131192BaseArray::SetMaxPrint(nprline*N,0);
    132193
    133 //-- Fill matrices
    134 double vmax=-1.;
    135 uint_8 k;
    136 uint_4 i,j; double s;
    137 uint_4 nbr=0, nbi=0;
    138 Vind = 0;
    139 if(nalea>0) for(i=0;i<nalea;i++) drand01();
    140 A = RandomSequence(RandomSequence::Flat,0.,1.);
    141 if(nbig>0) for(i=0;i<nbig;i++) {
    142   k=(uint_8)(drand01()*N*N); if(k>=N*N) k--;
    143   s=(drand01()>0.5)?1.:-1.;
    144   if(Vind[k]==0) {A[k] += (TYPE) s*vbig; Vind[k]=1; nbr++;}
    145 }
    146 A *= (TYPE) scale;
     194//--------------------------------------------------------
     195//-- Definition global arrays
     196//--------------------------------------------------------
     197TMatrix< TYPE > Ainput(N,N); Ainput = (TYPE) 0;
     198TMatrix< TYPE > A(N,N); A = (TYPE) 0;
     199Ainput.Show();
     200
     201//--------------------------------------------------------
     202//-- Fill matrices with flat random
     203//--------------------------------------------------------
     204if(gaussok) Ainput = RandomSequence(RandomSequence::Gaussian,0.,1.);
     205  else      Ainput = RandomSequence(RandomSequence::Flat,0.,1.);
    147206#if defined(COMPLEX)
    148 Vind = 0;
    149 B = RandomSequence(RandomSequence::Flat,0.,1.);
    150 if(nbig>0) for(i=0;i<nbig;i++) {
    151   k=(uint_8)(drand01()*N*N); if(k>=N*N) k--;
    152   s=(drand01()>0.5)?1.:-1.;
    153   if(Vind[k]==0) {B[k] += (TYPE) s*vbig; Vind[k]=1; nbi++;}
    154 }
    155 B *= (TYPE) scale;
    156 A += TYPE(0.,1.)*B;
    157 #endif
    158 cout<<"Nombre de valeurs BIG reelles = "<<nbr<<", imaginaires = "<<nbi<<endl;
    159 Adum = A;
    160 
    161 //-- Print matrice A
    162 cout<<"------------ TMatrix A :"<<endl;
    163 if(nprline>0) {cout<<A; cout<<endl<<endl;}
    164 
    165 //-- Inversion
    166 cout<<"------------ Inversion"<<endl;
    167 //InvA = Inverse(A);
    168 InvA = IdentityMatrix(1.,N);
    169 if(timok) PrtTim("--- End of filling ---");
    170 TYPE det = GausPiv(A,InvA,detok);
    171 if(timok) PrtTim("--- End of GausPiv ---");
    172 cout<<"Det = "<<det<<endl;
    173 cout<<"------------ TMatrix InvA = A^(-1):"<<endl;
    174 if(nprline>0) {cout<<InvA; cout<<endl<<endl;}
    175 
    176 //-- AiA = A * InvA
    177 cout<<"AiA = A * InvA"<<endl;
    178 AiA = A * InvA;
    179 cout<<"------------ TMatrix AiA = A * InvA:"<<endl;
    180 if(nprline>0) {cout<<AiA; cout<<endl;}
    181 
    182 //-- Check
    183 vmax=-1.;
    184 for(i=0;i<N;i++) {
    185   double absv = ABS_VAL( 1. - AiA(i,i) );
    186   if( absv > vmax ) vmax = absv;
    187 }
    188 cout<<"Ecart maximum par rapport a 1 sur diagonale = "<<vmax<<endl;
    189 vmax = -1.;
    190 for(i=0;i<N;i++) for(j=0;j<N;j++) {
    191   if( i == j ) continue;
    192   double absv = ABS_VAL( AiA(i,j) );
    193   if( absv > vmax ) vmax = absv;
    194 }
    195 cout<<"Ecart maximum par rapport a 0 hors diagonale = "<<vmax<<endl;
    196 
     207if(gaussok) A = RandomSequence(RandomSequence::Gaussian,0.,1.);
     208  else      A = RandomSequence(RandomSequence::Flat,0.,1.);
     209Ainput += TYPE(0.,1.)*A;
     210#endif
     211
     212//--------------------------------------------------------
     213//-- Fill matrices with big values
     214//--------------------------------------------------------
     215if(nbig>0) {
     216#if defined(COMPLEX)
     217  nbig = (nbig+1)/2;
     218#endif
     219  TMatrix< uint_2 > Vind(N,N); Vind = 0;
     220  // for real part
     221  uint_4 nbr=0;
     222  for(int k=0;k<nbig;k++) {
     223    int i = (int) (drand01()*N); int j = (int) (drand01()*N);
     224    double s=(drand01()>0.5)?1.:-1.;
     225    if(Vind(i,j)==0) {Ainput(i,j) += (TYPER) s*vbig; Vind(i,j)+=1; nbr++;}
     226  }
     227  cout<<"Nombre de valeurs BIG reelles = "<<nbr<<endl;
     228#if defined(COMPLEX)
     229  // for imaginary part
     230  uint_4 nbi=0;
     231  for(int k=0;k<nbig;k++) {
     232    int i = (int) (drand01()*N); int j = (int) (drand01()*N);
     233    double s=(drand01()>0.5)?1.:-1.;
     234    if(Vind(i,j)<=1) {Ainput(i,j) += TYPE(0.,(TYPER)s*vbig); Vind(i,j)+=2; nbi++;}
     235  }
     236  cout<<"Nombre de valeurs BIG imaginaires = "<<nbi<<endl;
     237  cout<<"Nombre de valeurs BIG = "<<nbr+nbi<<endl;
     238#endif
     239}
     240
     241//--------------------------------------------------------
     242//-- Scale matrix for machine precision tests
     243//--------------------------------------------------------
     244Ainput *= (TYPE) scale;
     245
     246//--------------------------------------------------------
     247//-- Create symetric matrix for all A if requested
     248//--------------------------------------------------------
     249if(symetok) Symetrize(Ainput);
     250
     251//--------------------------------------------------------
     252//-- Print matrice Ainput
     253//--------------------------------------------------------
     254cout<<"------------ TMatrix Ainput :"<<endl;
     255if(nprline>0) {cout<<Ainput; cout<<endl;}
     256PrtTim("--- End of Matrix filling ---");
     257
     258
     259#ifdef ALSO_LAPACK_INV
    197260////////////////////////////////////
    198261///////// Test avec Lapack /////////
    199262////////////////////////////////////
    200 #ifdef ALSO_LAPACK
    201 
    202 cout<<endl<<"LAPACK Cross-Check"<<endl;
    203 InvA = IdentityMatrix(1.,N);
    204 if(timok) PrtTim("--- End of check ---");
    205 LapackLinSolve(Adum,InvA);
    206 if(timok) PrtTim("--- End of LapackLinSolve ---");
    207 AiA = A * InvA;
    208 vmax=-1.;
    209 for(i=0;i<N;i++) {
    210   double absv = ABS_VAL( 1. - AiA(i,i) );
    211   if( absv > vmax ) vmax = absv;
    212 }
    213 cout<<"Ecart maximum par rapport a 1 sur diagonale = "<<vmax<<endl;
    214 vmax = -1.;
    215 for(i=0;i<N;i++) for(j=0;j<N;j++) {
    216   if( i == j ) continue;
    217   double absv = ABS_VAL( AiA(i,j) );
    218   if( absv > vmax ) vmax = absv;
    219 }
    220 cout<<"Ecart maximum par rapport a 0 hors diagonale = "<<vmax<<endl;
    221 
    222 #endif
    223 
    224 if(timok) PrtTim("--- End of Job ---");
     263{
     264cout<<"\n=========================================="<<endl;
     265cout<<"------------ Inversion LAPACK"<<endl;
     266A = Ainput;
     267//-- Inversion
     268TMatrix< TYPE > InvA(N,N); InvA = IdentityMatrix(1.,N);
     269int_4 info = LapackLinSolve(A,InvA);
     270cout<<"info="<<info<<endl;
     271PrtTim("--- End of LapackLinSolve Inversion ---");
     272//-- AiA = A * InvA
     273cout<<"Compute AiA = A * InvA"<<endl;
     274TMatrix< TYPE > AiA(N,N); AiA = Ainput * InvA;
     275cout<<"------------ TMatrix AiA = A * InvA:"<<endl;
     276if(nprline>0) {cout<<AiA; cout<<endl;}
     277//-- Check
     278Check_Mat_Ident(AiA);
     279PrtTim("--- End of LapackLinSolve Test ---");
     280}
     281#endif
     282
     283
     284#ifdef ALSO_LAPACK_INV_SYM
     285////////////////////////////////////////
     286///////// Test avec Lapack sym /////////
     287////////////////////////////////////////
     288{
     289cout<<"\n=========================================="<<endl;
     290cout<<"------------ Inversion LAPACK sym"<<endl;
     291TMatrix< TYPE > Asym(N,N); Asym=Ainput; Symetrize(Asym); A=Asym;
     292//-- Inversion
     293TMatrix< TYPE > InvA(N,N); InvA = IdentityMatrix(1.,N);
     294int_4 info = LapackLinSolveSym(A,InvA);
     295cout<<"info="<<info<<endl;
     296PrtTim("--- End of LapackLinSolveSym Inversion ---");
     297//-- AiA = A * InvA
     298cout<<"Compute AiA = A * InvA"<<endl;
     299TMatrix< TYPE > AiA(N,N); AiA = Asym * InvA;
     300cout<<"------------ TMatrix AiA = A * InvA:"<<endl;
     301if(nprline>0) {cout<<AiA; cout<<endl;}
     302//-- Check
     303Check_Mat_Ident(AiA);
     304PrtTim("--- End of LapackLinSolveSym Test ---");
     305}
     306#endif
     307
     308
     309#ifdef ALSO_LAPACK_INV_LSS
     310////////////////////////////////////////////////
     311///////// Test avec Lapack LeastSquare /////////
     312////////////////////////////////////////////////
     313{
     314cout<<"\n=========================================="<<endl;
     315cout<<"------------ Inversion LAPACK LeastSquare"<<endl;
     316A = Ainput;
     317//-- Inversion
     318TMatrix< TYPE > InvA(N,N); InvA = IdentityMatrix(1.,N);
     319int_4 info = LapackLeastSquareSolve(A,InvA);
     320cout<<"info="<<info<<endl;
     321PrtTim("--- End of LapackLeastSquareSolve Inversion ---");
     322//-- AiA = A * InvA
     323cout<<"Compute AiA = A * InvA"<<endl;
     324TMatrix< TYPE > AiA(N,N); AiA = Ainput * InvA;
     325cout<<"------------ TMatrix AiA = A * InvA:"<<endl;
     326if(nprline>0) {cout<<AiA; cout<<endl;}
     327//-- Check
     328Check_Mat_Ident(AiA);
     329PrtTim("--- End of LapackLeastSquareSolve Test ---");
     330}
     331#endif
     332
     333
     334#ifdef ALSO_LAPACK_EV
     335////////////////////////////////////////////////
     336///////// Test avec Lapack pour EV /////////
     337////////////////////////////////////////////////
     338{
     339cout<<"\n=========================================="<<endl;
     340cout<<"------------ Eigen decompositon LapackEigen "<<endl;
     341A=Ainput;
     342TMatrix< TYPE > Evec(N,N); Evec = (TYPE) 0;
     343TVector< complex<r_8> > Eval(N); Eval = complex<r_8>(0,0);
     344//-- Decompositon
     345int_4 info = LapackEigen(A,Eval,Evec,true);
     346cout<<"info="<<info<<endl;
     347PrtTim("--- End of LapackEigen decompositon ---");
     348if(nprline>0) {cout<<Eval; cout<<endl; cout<<Evec; cout<<endl;}
     349#ifndef COMPLEX
     350//-- Find the complex conjugate pairs
     351TVector< uint_2 > Evalconj(N); Evalconj = 0;
     352int_4 nconj=0;
     353for(int i=0;i<N-1;i++) {
     354  if(Evalconj(i)!=0) continue;   // deja traite
     355  if(Eval(i).imag()==0.) continue; // real eigenvalue
     356  if(fabs(Eval(i).imag()+Eval(i+1).imag())>1e-150) continue; // les 2 consecutives ne sont pas conjuguees
     357  if(Eval(i).imag()<0.) continue; // first conjugate have positive imaginary part
     358  if(Eval(i+1).imag()>0.) continue;
     359  Evalconj(i) = 1; Evalconj(i+1) = 2; nconj++;
     360}
     361//cout<<Evalconj<<endl;
     362cout<<"Number of conjugate eigen values: "<<nconj<<" *2 = "<<2*nconj<<" / "<<N<<endl;
     363#endif
     364//-- Azmlz = A*z(l) - l*z(l)
     365cout<<"Compute Azmlz(l) = A*z(l) - l*z(l)"<<endl;
     366TMatrix< complex<TYPER> > Azmlz(N,N); Azmlz = (complex<TYPER>) 0;
     367for(int l=0;l<N;l++) { // eigen value
     368  complex<TYPER> Eval_l = complex<TYPER>(Eval(l).real(),Eval(l).imag());
     369  for(int i=0;i<N;i++) {
     370    complex<TYPER> Evec_il;
     371    #ifdef COMPLEX
     372    Evec_il = Evec(i,l);
     373    #else
     374    Evec_il = complex<TYPER>(Evec(i,l),0.);
     375    if(Evalconj(l)==1)        Evec_il = complex<TYPER>(Evec(i,l),Evec(i,l+1));
     376      else if(Evalconj(l)==2) Evec_il = complex<TYPER>(Evec(i,l-1),-Evec(i,l));
     377    #endif
     378    for(int j=0;j<N;j++) {
     379      complex<TYPER> Evec_jl;
     380      #ifdef COMPLEX
     381      Evec_jl = Evec(j,l);
     382      #else
     383      Evec_jl = complex<TYPER>(Evec(j,l),0.);
     384      if(Evalconj(l)==1)        Evec_jl = complex<TYPER>(Evec(j,l),Evec(j,l+1));
     385        else if(Evalconj(l)==2) Evec_jl = complex<TYPER>(Evec(j,l-1),-Evec(j,l));
     386      #endif
     387      Azmlz(i,l) += Ainput(i,j) * Evec_jl;
     388    }
     389    Azmlz(i,l) -= Eval_l*Evec_il;
     390  }
     391}
     392if(nprline>0) {cout<<Azmlz; cout<<endl;}
     393//-- Check
     394Check_Mat_VecCol_2(Azmlz);
     395PrtTim("--- End of LapackEigen Test ---");
     396}
     397#endif
     398
     399
     400#ifdef ALSO_LAPACK_EV_SYM
     401////////////////////////////////////////////////
     402///////// Test avec Lapack sym pour EV /////////
     403////////////////////////////////////////////////
     404{
     405cout<<"\n=========================================="<<endl;
     406cout<<"------------ Eigen decompositon LapackEigenSym "<<endl;
     407TMatrix< TYPE > Aher(N,N); Aher=Ainput; Hermitian(Aher); A=Aher;
     408TVector<r_8> Eval;
     409//-- Decompositon
     410int_4 info = LapackEigenSym(A,Eval,true);
     411cout<<"info="<<info<<endl;
     412PrtTim("--- End of LapackEigenSym decompositon ---");
     413if(nprline>0) {cout<<Eval; cout<<endl; cout<<A; cout<<endl;}
     414//-- Azmlz = A*z(l) - l*z(l)
     415// le vecteur propre z pour la l-ieme valeur propre est dans A(.,l):
     416//    z_i = A(i,l)  ou "l" est la l-ieme valeur propre
     417cout<<"Compute Azmlz(l) = A*z(l) - l*z(l)"<<endl;
     418TMatrix< TYPE > Azmlz(N,N); Azmlz = (TYPE) 0;
     419for(int l=0;l<N;l++)  // eigen value
     420  for(int i=0;i<N;i++)
     421    {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);}
     422if(nprline>0) {cout<<Azmlz; cout<<endl;}
     423//-- Check
     424Check_Mat_VecCol_0(Azmlz);
     425PrtTim("--- End of LapackEigenSym Test ---");
     426}
     427#endif
     428
     429
     430#ifdef USE_GAUSPIV
     431////////////////////////////////////
     432///////// Test avec GausPiv /////////
     433////////////////////////////////////
     434{
     435cout<<"\n==========================================\n"
     436    <<"------------ Inversion GausPiv"<<endl;
     437SimpleMatrixOperation< TYPE >::SetGausPivScal(tscal);
     438A = Ainput;
     439//-- Inversion
     440TMatrix< TYPE > InvA(N,N); InvA = IdentityMatrix(1.,N);
     441TYPE det = GausPiv(A,InvA,detok);
     442PrtTim("--- End of GausPiv Inversion ---");
     443cout<<"Det = "<<det<<endl;
     444cout<<"------------ TMatrix InvA = A^(-1):"<<endl;
     445//-- AiA = A * InvA
     446cout<<"Compute AiA = A * InvA"<<endl;
     447TMatrix< TYPE > AiA(N,N); AiA = Ainput * InvA;
     448cout<<"------------ TMatrix AiA = A * InvA:"<<endl;
     449if(nprline>0) {cout<<AiA; cout<<endl;}
     450//-- Check
     451Check_Mat_Ident(AiA);
     452PrtTim("--- End of GausPiv Test ---");
     453}
     454#endif
     455
     456
     457PrtTim("--- End of Job ---");
    225458exit(0);
    226459}
     460
     461
     462
     463////////////////////////////////////////////////////////////
     464////-------------------------------------------------------
     465void Symetrize(TMatrix< TYPE >& A)
     466// Symetrize A
     467{
     468 int_4 N = A.NRows();
     469 for(int i=0;i<N-1;i++) for(int j=i+1;j<N;j++) A(j,i) = A(i,j);
     470}
     471
     472////-------------------------------------------------------
     473void Hermitian(TMatrix< TYPE >& A)
     474// Put A hermitian
     475{
     476 int_4 N = A.NRows();
     477 for(int i=0;i<N-1;i++) for(int j=i+1;j<N;j++) A(j,i) = CONJ_VAL(A(i,j));
     478 for(int i=0;i<N;i++) A(i,i) = REAL_PART(A(i,i));
     479}
     480
     481////-------------------------------------------------------
     482r_8  Check_Mat_Ident(TMatrix< TYPE >& A)
     483// Compute the biggest difference element by element of A / Identity
     484{
     485 int_4 N = A.NRows();
     486 r_8 vmaxd=-1.;
     487 for(int i=0;i<N;i++)
     488   if( ABS_VAL((TYPER)1.-A(i,i)) > vmaxd ) vmaxd = ABS_VAL((TYPER)1.-A(i,i));
     489 cout<<"Ecart maximum par rapport a 1 sur diagonale = "<<vmaxd<<endl;
     490 r_8 vmaxh = -1.;
     491 for(int i=0;i<N;i++) for(int j=0;j<N;j++) {
     492   if(i==j) continue;
     493   if( ABS_VAL(A(i,j)) > vmaxh ) vmaxh = ABS_VAL(A(i,j));
     494 }
     495 cout<<"Ecart maximum par rapport a 0 hors diagonale = "<<vmaxh<<endl;
     496 return (vmaxd>vmaxh)? vmaxd: vmaxh;
     497}
     498
     499////-------------------------------------------------------
     500r_8  Check_Mat_VecCol_0(TMatrix< TYPE >& A)
     501// Return the biggest norm of the N vectors column of matrix
     502{
     503 int_4 N = A.NRows();
     504 r_8 vmax=-1.;
     505 for(int l=0;l<N;l++) {
     506   r_8 absv = 0.;
     507   for(int i=0;i<N;i++) absv += ABS_VAL(A(i,l)) * ABS_VAL(A(i,l));
     508   if( absv > vmax ) vmax = absv;
     509 }
     510 vmax = sqrt(vmax);
     511 cout<<"Longueur max de ||A*z-l*z|| pour tous l = "<<vmax<<endl;
     512 return vmax;
     513}
     514
     515////-------------------------------------------------------
     516void  Check_Mat_VecCol_2(TMatrix< complex<TYPER> >& A)
     517// Return the biggest norm of :
     518//   - the real part of the N vectors column of matrix
     519//   - the imaginary part of the N vectors column of matrix
     520//   - the module of the N vectors column of matrix
     521{
     522 int_4 N = A.NRows();
     523 r_8 vmaxr=-1., vmaxi=-1., vmaxn=-1.;
     524 for(int l=0;l<N;l++) {
     525   double absvr = 0., absvi = 0., absvn = 0.;
     526   for(int i=0;i<N;i++) {
     527     absvr += A(i,l).real()*A(i,l).real();
     528     absvi += A(i,l).imag()*A(i,l).imag();
     529     absvn += A(i,l).real()*A(i,l).real() + A(i,l).imag()*A(i,l).imag();
     530   }
     531   if( absvr > vmaxr ) vmaxr = absvr;
     532   if( absvi > vmaxi ) vmaxi = absvi;
     533   if( absvn > vmaxn ) vmaxn = absvn;
     534 }
     535 vmaxr=sqrt(vmaxr); vmaxi=sqrt(vmaxi); vmaxn=sqrt(vmaxn);
     536 cout<<"Longueur max de ||A*z-l*z|| pour tous l, reel = "<<vmaxr
     537     <<", imag = "<<vmaxi<<", module = "<<vmaxn<<endl;
     538}
     539
     540
     541/*
     542void TestIlaEnv(int_4 n)
     543{
     544LapackServer<TYPE> lps;
     545cout<<"TestIlaEnv n="<<n<<endl;
     546cout<<lps.ilaenv_en_C(1,"SSYTRF","U",n,-1,-1,-1)<<endl;
     547cout<<lps.ilaenv_en_C(1,"DSYTRF","U",n,-1,-1,-1)<<endl;
     548cout<<lps.ilaenv_en_C(1,"CSYTRF","U",n,-1,-1,-1)<<endl;
     549cout<<lps.ilaenv_en_C(1,"ZSYTRF","U",n,-1,-1,-1)<<endl;
     550cout<<lps.ilaenv_en_C(1,"SSYTRF","L",n,-1,-1,-1)<<endl;
     551cout<<lps.ilaenv_en_C(1,"DSYTRF","L",n,-1,-1,-1)<<endl;
     552cout<<lps.ilaenv_en_C(1,"CSYTRF","L",n,-1,-1,-1)<<endl;
     553cout<<lps.ilaenv_en_C(1,"ZSYTRF","L",n,-1,-1,-1)<<endl;
     554}
     555*/
Note: See TracChangeset for help on using the changeset viewer.