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

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

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