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

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

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