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

Last change on this file since 3186 was 2615, checked in by cmv, 21 years ago

using namespace sophya enleve de machdefs.h, nouveau sopnamsp.h cmv 10/09/2004
File size: 20.6 KB
Line 
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 ?
11//#define COMPLEX
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#define ALSO_LAPACK_INV
24#define ALSO_LAPACK_INV_SYM
25#define ALSO_LAPACK_INV_LSS
26#define ALSO_LAPACK_INV_LSS_SVD
27#define ALSO_LAPACK_EV
28#define ALSO_LAPACK_EV_SYM
29#define ALSO_LAPACK_SVD
30#define ALSO_LAPACK_SVD_DC
31
32//////////////////////////////////////////////////
33//////////////////////////////////////////////////
34// NE RIEN CHANGER CI-APRES //
35//////////////////////////////////////////////////
36//////////////////////////////////////////////////
37
38//////////////////////////////////////////////////
39#include "sopnamsp.h"
40#include "machdefs.h"
41#include <iostream>
42#include <stdlib.h>
43#include <stdio.h>
44#include <string.h>
45#include <math.h>
46#include <unistd.h>
47#include "timing.h"
48#include "ntoolsinit.h"
49#include "pexceptions.h"
50#include "array.h"
51#include "srandgen.h"
52#if defined(USE_LAPACK)
53#include "intflapack.h"
54#endif
55
56//////////////////////////////////////////////////
57#if defined(COMPLEX)
58 #if defined(PRECIS32)
59 #define TYPE complex<r_4>
60 #define TYPER r_4
61 #else
62 #define TYPE complex<r_8>
63 #define TYPER r_8
64 #endif
65 #define REAL_PART(_x_) (TYPE((_x_).real(),0.))
66 #define CONJ_VAL(_x_) (TYPE((_x_).real(),-(_x_).imag()))
67 #define ABS_VAL(_x_) sqrt((double)((_x_).real()*(_x_).real() + (_x_).imag()*(_x_).imag()))
68#else
69 #if defined(PRECIS32)
70 #define TYPE r_4
71 #define TYPER r_4
72 #else
73 #define TYPE r_8
74 #define TYPER r_8
75 #endif
76 #define REAL_PART(_x_) (_x_)
77 #define CONJ_VAL(_x_) (_x_)
78 #define ABS_VAL(_x_) fabs((double)_x_)
79#endif
80
81//////////////////////////////////////////////////
82void Symetrize(TMatrix< TYPE >& A);
83void Hermitian(TMatrix< TYPE >& A);
84r_8 Check_Mat_Ident(TMatrix< TYPE >& A);
85r_8 Check_Mat_Zero(TMatrix< TYPE >& A);
86r_8 Check_Mat_VecCol_0(TMatrix< TYPE >& A);
87void Check_Mat_VecCol_2(TMatrix< complex<TYPER> >& A);
88
89
90//////////////////////////////////////////////////
91int main(int narg,char *arg[])
92{
93//--------------------------------------------------------
94//-- Initialisation
95//--------------------------------------------------------
96// number of lines/columns
97uint_4 N = 5;
98// scale of the value (if =1 values between -1 and 1)
99r_8 scale = 1.;
100// number of values change by +/- vbig
101uint_4 nbig = N;
102r_8 vbig = 100.;
103// Nombre de lignes de matrice a imprimer
104uint_4 nprline = N;
105// Initialisation du pauvre de l'aleatoire
106uint_4 nalea = 0;
107// data scaling for gauss pivoting and determinant
108int tscal = 1;
109bool detok=false;
110// Please symetrize the input matrice
111bool symetok=false;
112// Please symetrize the input matrice
113bool gaussok=false;
114
115//--------------------------------------------------------
116//-- Decodage arguments
117//--------------------------------------------------------
118char c;
119while((c = getopt(narg,arg,"Sdgn:s:b:l:a:t:h")) != -1) {
120 switch (c) {
121 case 'S' :
122 symetok = true;
123 break;
124 case 'd' :
125 detok = true;
126 break;
127 case 'g' :
128 gaussok = true;
129 break;
130 case 'n' :
131 sscanf(optarg,"%d",&N);
132 break;
133 case 's' :
134 sscanf(optarg,"%lf",&scale);
135 break;
136 case 'b' :
137 sscanf(optarg,"%d,%lf",&nbig,&vbig);
138 break;
139 case 'l' :
140 sscanf(optarg,"%d",&nprline);
141 break;
142 case 'a' :
143 sscanf(optarg,"%d",&nalea);
144 break;
145 case 't' :
146 sscanf(optarg,"%d",&tscal);
147 break;
148 case 'h' :
149 cout<<"tsttminv [-h] [-n N] [-S] [-s scale] [-b nbig,vbig] [-g]"<<endl
150 <<" [-l nprline] [-a nalea] [-t tscal] [-d]"<<endl
151 <<"-- matrix A(N,N) filled with {[-1,1] +/- vbig(nbig time)}*scale --"<<endl
152 <<"-g : instead of flat [-1,1] use normal gaussian distribution for A(i,j)"<<endl
153 <<"-S : symetrize the input matrix"<<endl
154 <<"-l : print nprline of input and test matrices"<<endl
155 <<"-a : for random (pseudo) changing"<<endl
156 <<"-- Only GausPiv --"<<endl
157 <<"-t 0/1/2 : data scaling 0=no, 1=global (def), 2=row-by-row"<<endl
158 <<"-d : also compute determinant (ne pas utiliser si N est grand)"<<endl;
159 return(-1);
160 }
161}
162if(N<=1) N = 1;
163cout<<"Taille matrice NxN, N = "<<N<<endl;
164if(gaussok) cout<<"Elements gaussian normal * scale = "<<scale<<endl;
165 else cout<<"Elements entre +/- 1 * scale = "<<scale<<endl;
166cout<<"Nombre de valeurs hors standard nbig = "<<nbig<<endl;
167cout<<"Valeurs hors standard (+/- vbig = "<<vbig<<" ) * scale = "<<vbig*scale<<endl;
168cout<<"Nombre de lignes de matrice a imprimer "<<nprline<<endl;
169cout<<"Initialisation de l aleatoire par "<<nalea<<" tirages"<<endl;
170cout<<"Data scaling "<<tscal<<" determinant="<<detok<<endl;
171if(symetok) cout<<"Input matrix has been symetrized "<<symetok<<endl;
172cout<<endl;
173
174//--------------------------------------------------------
175//-- Initialization arrays
176//--------------------------------------------------------
177SophyaInit();
178InitTim();
179#if defined(USE_LAPACK)
180 BaseArray::SetDefaultMemoryMapping(BaseArray::FortranMemoryMapping);
181#endif
182if(nalea>0) for(int i=0;i<nalea;i++) drand01();
183BaseArray::SetMaxPrint(nprline*N,0);
184
185//--------------------------------------------------------
186//-- Definition global arrays
187//--------------------------------------------------------
188TMatrix< TYPE > Ainput(N,N); Ainput = (TYPE) 0;
189TMatrix< TYPE > A(N,N); A = (TYPE) 0;
190Ainput.Show();
191
192//--------------------------------------------------------
193//-- Fill matrices with flat random
194//--------------------------------------------------------
195if(gaussok) Ainput = RandomSequence(RandomSequence::Gaussian,0.,1.);
196 else Ainput = RandomSequence(RandomSequence::Flat,0.,1.);
197#if defined(COMPLEX)
198if(gaussok) A = RandomSequence(RandomSequence::Gaussian,0.,1.);
199 else A = RandomSequence(RandomSequence::Flat,0.,1.);
200Ainput += TYPE(0.,1.)*A;
201#endif
202
203//--------------------------------------------------------
204//-- Fill matrices with big values
205//--------------------------------------------------------
206if(nbig>0) {
207#if defined(COMPLEX)
208 nbig = (nbig+1)/2;
209#endif
210 TMatrix< uint_2 > Vind(N,N); Vind = 0;
211 // for real part
212 uint_4 nbr=0;
213 for(int k=0;k<nbig;k++) {
214 int i = (int) (drand01()*N); int j = (int) (drand01()*N);
215 double s=(drand01()>0.5)?1.:-1.;
216 if(Vind(i,j)==0) {Ainput(i,j) += (TYPER) s*vbig; Vind(i,j)+=1; nbr++;}
217 }
218 cout<<"Nombre de valeurs BIG reelles = "<<nbr<<endl;
219#if defined(COMPLEX)
220 // for imaginary part
221 uint_4 nbi=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)<=1) {Ainput(i,j) += TYPE(0.,(TYPER)s*vbig); Vind(i,j)+=2; nbi++;}
226 }
227 cout<<"Nombre de valeurs BIG imaginaires = "<<nbi<<endl;
228 cout<<"Nombre de valeurs BIG = "<<nbr+nbi<<endl;
229#endif
230}
231
232//--------------------------------------------------------
233//-- Scale matrix for machine precision tests
234//--------------------------------------------------------
235Ainput *= (TYPE) scale;
236
237//--------------------------------------------------------
238//-- Create symetric matrix for all A if requested
239//--------------------------------------------------------
240if(symetok) Symetrize(Ainput);
241
242//--------------------------------------------------------
243//-- Print matrice Ainput
244//--------------------------------------------------------
245cout<<"------------ TMatrix Ainput :"<<endl;
246if(nprline>0) {cout<<Ainput; cout<<endl;}
247PrtTim("--- End of Matrix filling ---");
248
249
250//////////////////////////////////////////////
251///////// Test Inversion avec Lapack /////////
252//////////////////////////////////////////////
253#if defined(USE_LAPACK) && defined(ALSO_LAPACK_INV)
254{
255cout<<"\n=========================================="<<endl;
256cout<<"------------ Inversion LAPACK (LU factorization)"<<endl;
257A = Ainput;
258//-- Inversion
259TMatrix< TYPE > InvA(N,N); InvA = IdentityMatrix(1.,N);
260int_4 info = LapackLinSolve(A,InvA);
261cout<<"info="<<info<<endl;
262PrtTim("--- End of LapackLinSolve Inversion ---");
263//-- AiA = A * InvA
264cout<<"Compute AiA = A * InvA"<<endl;
265TMatrix< TYPE > AiA(N,N); AiA = Ainput * InvA;
266if(nprline>0) {cout<<AiA; cout<<endl;}
267//-- Check
268Check_Mat_Ident(AiA);
269PrtTim("--- End of LapackLinSolve Test ---");
270}
271#endif
272
273
274//////////////////////////////////////////////////
275///////// Test Inversion avec Lapack sym /////////
276//////////////////////////////////////////////////
277#if defined(USE_LAPACK) && defined(ALSO_LAPACK_INV_SYM)
278{
279cout<<"\n=========================================="<<endl;
280cout<<"------------ Inversion LAPACK symetric (LU factorization)"<<endl;
281TMatrix< TYPE > Asym(N,N); Asym=Ainput; Symetrize(Asym); A=Asym;
282//-- Inversion
283TMatrix< TYPE > InvA(N,N); InvA = IdentityMatrix(1.,N);
284int_4 info = LapackLinSolveSym(A,InvA);
285cout<<"info="<<info<<endl;
286PrtTim("--- End of LapackLinSolveSym Inversion ---");
287//-- AiA = A * InvA
288cout<<"Compute AiA = A * InvA"<<endl;
289TMatrix< TYPE > AiA(N,N); AiA = Asym * InvA;
290cout<<"------------ TMatrix AiA = A * InvA:"<<endl;
291if(nprline>0) {cout<<AiA; cout<<endl;}
292//-- Check
293Check_Mat_Ident(AiA);
294PrtTim("--- End of LapackLinSolveSym Test ---");
295}
296#endif
297
298
299////////////////////////////////////////////////
300///////// Test avec Lapack LeastSquare /////////
301////////////////////////////////////////////////
302#if defined(USE_LAPACK) && defined(ALSO_LAPACK_INV_LSS)
303{
304cout<<"\n=========================================="<<endl;
305cout<<"------------ Inversion LAPACK LeastSquare (QR or LQ factorization)"<<endl;
306A = Ainput;
307//-- Inversion
308TMatrix< TYPE > InvA(N,N); InvA = IdentityMatrix(1.,N);
309int_4 info = LapackLeastSquareSolve(A,InvA);
310cout<<"info="<<info<<endl;
311PrtTim("--- End of LapackLeastSquareSolve Inversion ---");
312//-- AiA = A * InvA
313cout<<"Compute AiA = A * InvA"<<endl;
314TMatrix< TYPE > AiA(N,N); AiA = Ainput * InvA;
315if(nprline>0) {cout<<AiA; cout<<endl;}
316//-- Check
317Check_Mat_Ident(AiA);
318PrtTim("--- End of LapackLeastSquareSolve Test ---");
319}
320#endif
321
322
323///////////////////////////////////////////////////////
324///////// Test avec Lapack LeastSquare by SVD /////////
325///////////////////////////////////////////////////////
326#if defined(USE_LAPACK) && defined(ALSO_LAPACK_INV_LSS_SVD)
327{
328cout<<"\n=========================================="<<endl;
329cout<<"------------ Inversion LAPACK LeastSquare (SVD decomposition)"<<endl;
330A = Ainput;
331//-- Inversion
332TMatrix< TYPE > InvA(N,N); InvA = IdentityMatrix(1.,N);
333TVector<r_8> S;
334int_4 rank;
335r_8 rcond = -1.;
336int_4 info = LapackLeastSquareSolveSVD_DC(A,InvA,S,rank,rcond);
337cout<<"info="<<info<<" (rank="<<rank<<")"<<endl;
338PrtTim("--- End of LapackLeastSquareSolveSVD_DC Inversion ---");
339if(nprline>0) {cout<<S; cout<<endl;}
340double smax = fabs(S(0)), smin = fabs(S(S.Size()-1));
341cout<<" Smin = |"<<S(S.Size()-1)<<"| = "<<smin<<", "
342 <<" Smax = |"<<S(0)<<"| = "<<smax<<", "
343 <<" --> Smin/Smax = "<<smin/smax<<endl;
344//-- AiA = A * InvA
345cout<<"Compute AiA = A * InvA"<<endl;
346TMatrix< TYPE > AiA(N,N); AiA = Ainput * InvA;
347if(nprline>0) {cout<<AiA; cout<<endl;}
348//-- Check
349Check_Mat_Ident(AiA);
350PrtTim("--- End of LapackLeastSquareSolveSVD_DC Test ---");
351}
352#endif
353
354
355////////////////////////////////////////////////
356///////// Test avec Lapack pour EV /////////
357////////////////////////////////////////////////
358#if defined(USE_LAPACK) && defined(ALSO_LAPACK_EV)
359{
360cout<<"\n=========================================="<<endl;
361cout<<"------------ Eigen decompositon LapackEigen "<<endl;
362A=Ainput;
363TMatrix< TYPE > Evec;
364TVector< complex<r_8> > Eval;
365//-- Decompositon
366int_4 info = LapackEigen(A,Eval,Evec,true);
367cout<<"info="<<info<<endl;
368PrtTim("--- End of LapackEigen decompositon ---");
369if(nprline>0) {cout<<Eval; cout<<endl; cout<<Evec; cout<<endl;}
370//-- Find the complex conjugate pairs
371#if !defined(COMPLEX)
372TVector< uint_2 > Evalconj(N); Evalconj = 0;
373int_4 nconj=0;
374for(int i=0;i<N-1;i++) {
375 if(Evalconj(i)!=0) continue; // deja traite
376 if(Eval(i).imag()==0.) continue; // real eigenvalue
377 if(fabs(Eval(i).imag()+Eval(i+1).imag())>1e-150) continue; // les 2 consecutives ne sont pas conjuguees
378 if(Eval(i).imag()<0.) continue; // first conjugate have positive imaginary part
379 if(Eval(i+1).imag()>0.) continue;
380 Evalconj(i) = 1; Evalconj(i+1) = 2; nconj++;
381}
382//cout<<Evalconj<<endl;
383cout<<"Number of conjugate eigen values: "<<nconj<<" *2 = "<<2*nconj<<" / "<<N<<endl;
384#endif
385//-- Azmlz = A*z(l) - l*z(l)
386cout<<"Compute Azmlz(l) = A*z(l) - l*z(l)"<<endl;
387TMatrix< complex<TYPER> > Azmlz(N,N); Azmlz = (complex<TYPER>) 0;
388for(int l=0;l<N;l++) { // eigen value
389 complex<TYPER> Eval_l = complex<TYPER>(Eval(l).real(),Eval(l).imag());
390 for(int i=0;i<N;i++) {
391 complex<TYPER> Evec_il;
392 #if defined(COMPLEX)
393 Evec_il = Evec(i,l);
394 #else
395 Evec_il = complex<TYPER>(Evec(i,l),0.);
396 if(Evalconj(l)==1) Evec_il = complex<TYPER>(Evec(i,l),Evec(i,l+1));
397 else if(Evalconj(l)==2) Evec_il = complex<TYPER>(Evec(i,l-1),-Evec(i,l));
398 #endif
399 for(int j=0;j<N;j++) {
400 complex<TYPER> Evec_jl;
401 #if defined(COMPLEX)
402 Evec_jl = Evec(j,l);
403 #else
404 Evec_jl = complex<TYPER>(Evec(j,l),0.);
405 if(Evalconj(l)==1) Evec_jl = complex<TYPER>(Evec(j,l),Evec(j,l+1));
406 else if(Evalconj(l)==2) Evec_jl = complex<TYPER>(Evec(j,l-1),-Evec(j,l));
407 #endif
408 Azmlz(i,l) += Ainput(i,j) * Evec_jl;
409 }
410 Azmlz(i,l) -= Eval_l*Evec_il;
411 }
412}
413if(nprline>0) {cout<<Azmlz; cout<<endl;}
414//-- Check
415Check_Mat_VecCol_2(Azmlz);
416PrtTim("--- End of LapackEigen Test ---");
417}
418#endif
419
420
421////////////////////////////////////////////////
422///////// Test avec Lapack sym pour EV /////////
423////////////////////////////////////////////////
424#if defined(USE_LAPACK) && defined(ALSO_LAPACK_EV_SYM)
425{
426cout<<"\n=========================================="<<endl;
427cout<<"------------ Eigen decompositon LapackEigenSym "<<endl;
428TMatrix< TYPE > Aher(N,N); Aher=Ainput; Hermitian(Aher); A=Aher;
429TVector<r_8> Eval;
430//-- Decompositon
431int_4 info = LapackEigenSym(A,Eval,true);
432cout<<"info="<<info<<endl;
433PrtTim("--- End of LapackEigenSym decompositon ---");
434if(nprline>0) {cout<<Eval; cout<<endl; cout<<A; cout<<endl;}
435//-- Azmlz = A*z(l) - l*z(l)
436// le vecteur propre z pour la l-ieme valeur propre est dans A(.,l):
437// z_i = A(i,l) ou "l" est la l-ieme valeur propre
438cout<<"Compute Azmlz(l) = A*z(l) - l*z(l)"<<endl;
439TMatrix< TYPE > Azmlz(N,N); Azmlz = (TYPE) 0;
440for(int l=0;l<N;l++) // eigen value
441 for(int i=0;i<N;i++)
442 {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);}
443if(nprline>0) {cout<<Azmlz; cout<<endl;}
444//-- Check
445Check_Mat_VecCol_0(Azmlz);
446PrtTim("--- End of LapackEigenSym Test ---");
447}
448#endif
449
450
451////////////////////////////////////////////////
452///////////// Test avec Lapack SVD /////////////
453////////////////////////////////////////////////
454#if defined(USE_LAPACK) && defined(ALSO_LAPACK_SVD)
455{
456cout<<"\n=========================================="<<endl;
457cout<<"------------ SVD decompositon LapackSVD "<<endl;
458A=Ainput;
459TVector< TYPE > S; TMatrix< TYPE > U; TMatrix< TYPE > VT;
460//-- Decompositon
461int_4 info = LapackSVD(A,S,U,VT);
462cout<<"info="<<info<<endl;
463PrtTim("--- End of LapackSVD decompositon ---");
464if(nprline>0) {cout<<S; cout<<endl;}
465double smax = ABS_VAL(S(0)), smin = ABS_VAL(S(N-1));
466cout<<" Smin = |"<<S(N-1)<<"| = "<<smin<<", "
467 <<" Smax = |"<<S(0)<<"| = "<<smax<<", "
468 <<" --> Smin/Smax = "<<smin/smax<<endl;
469//-- A = U*S*VT
470cout<<"Compute A = U*S*VT"<<endl;
471TMatrix< TYPE > AmUSVt(N,N); AmUSVt = U;
472for(int i=0;i<N;i++) for(int j=0;j<N;j++) AmUSVt(i,j) *= S(j);
473AmUSVt = AmUSVt * VT; AmUSVt -= Ainput;
474if(nprline>0) {cout<<AmUSVt; cout<<endl;}
475//-- Check
476Check_Mat_Zero(AmUSVt);
477PrtTim("--- End of LapackSVD Test ---");
478}
479#endif
480
481
482///////////////////////////////////////////////////
483///////////// Test avec Lapack SVD_DC /////////////
484///////////////////////////////////////////////////
485#if defined(USE_LAPACK) && defined(ALSO_LAPACK_SVD_DC)
486{
487cout<<"\n=========================================="<<endl;
488cout<<"------------ SVD decompositon LapackSVD_DC "<<endl;
489A=Ainput;
490TVector< r_8 > S; TMatrix< TYPE > U; TMatrix< TYPE > VT;
491//-- Decompositon
492int_4 info = LapackSVD_DC(A,S,U,VT);
493cout<<"info="<<info<<endl;
494PrtTim("--- End of LapackSVD_DC decompositon ---");
495if(nprline>0) {cout<<S; cout<<endl;}
496double smax = fabs(S(0)), smin = fabs(S(N-1));
497 cout<<"CMV: "<<S(0)<<endl;
498//cout<<" Smin = |"<<S(N-1)<<"| = "<<smin<<", "
499// <<" Smax = |"<<S(0)<<"| = "<<smax<<", "
500// <<" --> Smin/Smax = "<<smin/smax<<endl;
501//-- A = U*S*VT
502cout<<"Compute A = U*S*VT"<<endl;
503TMatrix< TYPE > AmUSVt(N,N); AmUSVt = U;
504for(int i=0;i<N;i++) for(int j=0;j<N;j++) AmUSVt(i,j) *= (TYPE) S(j);
505AmUSVt = AmUSVt * VT; AmUSVt -= Ainput;
506if(nprline>0) {cout<<AmUSVt; cout<<endl;}
507//-- Check
508Check_Mat_Zero(AmUSVt);
509PrtTim("--- End of LapackSVD_DC Test ---");
510}
511#endif
512
513
514///////////////////////////////////////////////
515///////// Test Inversion avec GausPiv /////////
516///////////////////////////////////////////////
517#if defined(USE_GAUSPIV)
518{
519cout<<"\n==========================================\n"
520 <<"------------ Inversion GausPiv"<<endl;
521SimpleMatrixOperation< TYPE >::SetGausPivScal(tscal);
522A = Ainput;
523//-- Inversion
524TMatrix< TYPE > InvA(N,N); InvA = IdentityMatrix(1.,N);
525TYPE det = GausPiv(A,InvA,detok);
526PrtTim("--- End of GausPiv Inversion ---");
527cout<<"Det = "<<det<<endl;
528cout<<"------------ TMatrix InvA = A^(-1):"<<endl;
529//-- AiA = A * InvA
530cout<<"Compute AiA = A * InvA"<<endl;
531TMatrix< TYPE > AiA(N,N); AiA = Ainput * InvA;
532cout<<"------------ TMatrix AiA = A * InvA:"<<endl;
533if(nprline>0) {cout<<AiA; cout<<endl;}
534//-- Check
535Check_Mat_Ident(AiA);
536PrtTim("--- End of GausPiv Test ---");
537}
538#endif
539
540
541PrtTim("--- End of Job ---");
542exit(0);
543}
544
545
546
547////////////////////////////////////////////////////////////
548////-------------------------------------------------------
549void Symetrize(TMatrix< TYPE >& A)
550// Symetrize A
551{
552 int_4 N = A.NRows();
553 for(int i=0;i<N-1;i++) for(int j=i+1;j<N;j++) A(j,i) = A(i,j);
554}
555
556////-------------------------------------------------------
557void Hermitian(TMatrix< TYPE >& A)
558// Put A hermitian
559{
560 int_4 N = A.NRows();
561 for(int i=0;i<N-1;i++) for(int j=i+1;j<N;j++) A(j,i) = CONJ_VAL(A(i,j));
562 for(int i=0;i<N;i++) A(i,i) = REAL_PART(A(i,i));
563}
564
565////-------------------------------------------------------
566r_8 Check_Mat_Ident(TMatrix< TYPE >& A)
567// Compute the biggest difference element by element of A / Identity
568{
569 int_4 N = A.NRows();
570 r_8 vmaxd=-1.;
571 for(int i=0;i<N;i++)
572 if( ABS_VAL((TYPER)1.-A(i,i)) > vmaxd ) vmaxd = ABS_VAL((TYPER)1.-A(i,i));
573 cout<<"Ecart maximum par rapport a 1 sur diagonale = "<<vmaxd<<endl;
574 r_8 vmaxh = -1.;
575 for(int i=0;i<N;i++) for(int j=0;j<N;j++) {
576 if(i==j) continue;
577 if( ABS_VAL(A(i,j)) > vmaxh ) vmaxh = ABS_VAL(A(i,j));
578 }
579 cout<<"Ecart maximum par rapport a 0 hors diagonale = "<<vmaxh<<endl;
580 return (vmaxd>vmaxh)? vmaxd: vmaxh;
581}
582
583////-------------------------------------------------------
584r_8 Check_Mat_Zero(TMatrix< TYPE >& A)
585// Compute the biggest difference element by element of A / Zero matrix
586{
587 int_4 N = A.NRows();
588 r_8 vmax = -1.;
589 for(int i=0;i<N;i++) for(int j=0;j<N;j++)
590 if( ABS_VAL(A(i,j)) > vmax ) vmax = ABS_VAL(A(i,j));
591 cout<<"Ecart maximum par rapport a zero = "<<vmax<<endl;
592 return vmax;
593}
594
595////-------------------------------------------------------
596r_8 Check_Mat_VecCol_0(TMatrix< TYPE >& A)
597// Return the biggest norm of the N vectors column of matrix
598{
599 int_4 N = A.NRows();
600 r_8 vmax=-1.;
601 for(int l=0;l<N;l++) {
602 r_8 absv = 0.;
603 for(int i=0;i<N;i++) absv += ABS_VAL(A(i,l)) * ABS_VAL(A(i,l));
604 if( absv > vmax ) vmax = absv;
605 }
606 vmax = sqrt(vmax);
607 cout<<"Longueur max de ||A*z-l*z|| pour tous l = "<<vmax<<endl;
608 return vmax;
609}
610
611////-------------------------------------------------------
612void Check_Mat_VecCol_2(TMatrix< complex<TYPER> >& A)
613// Return the biggest norm of :
614// - the real part of the N vectors column of matrix
615// - the imaginary part of the N vectors column of matrix
616// - the module of the N vectors column of matrix
617{
618 int_4 N = A.NRows();
619 r_8 vmaxr=-1., vmaxi=-1., vmaxn=-1.;
620 for(int l=0;l<N;l++) {
621 double absvr = 0., absvi = 0., absvn = 0.;
622 for(int i=0;i<N;i++) {
623 absvr += A(i,l).real()*A(i,l).real();
624 absvi += A(i,l).imag()*A(i,l).imag();
625 absvn += A(i,l).real()*A(i,l).real() + A(i,l).imag()*A(i,l).imag();
626 }
627 if( absvr > vmaxr ) vmaxr = absvr;
628 if( absvi > vmaxi ) vmaxi = absvi;
629 if( absvn > vmaxn ) vmaxn = absvn;
630 }
631 vmaxr=sqrt(vmaxr); vmaxi=sqrt(vmaxi); vmaxn=sqrt(vmaxn);
632 cout<<"Longueur max de ||A*z-l*z|| pour tous l, reel = "<<vmaxr
633 <<", imag = "<<vmaxi<<", module = "<<vmaxn<<endl;
634}
Note: See TracBrowser for help on using the repository browser.