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

Last change on this file since 2598 was 2580, checked in by cmv, 21 years ago

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