Changeset 2558 in Sophya for trunk/SophyaProg/Tests/tsttminv.cc

Timestamp:

Jul 22, 2004, 3:38:05 PM (21 years ago)

Author:

cmv

Message:

add test of svd decomposition (cmv 22/07/04)

File:

• trunk/SophyaProg/Tests/tsttminv.cc (modified) (6 diffs)

trunk/SophyaProg/Tests/tsttminv.cc

@@ -26 +26 @@
 #define ALSO_LAPACK_EV
 #define ALSO_LAPACK_EV_SYM
+#define ALSO_LAPACK_SVD
 
 //////////////////////////////////////////////////
@@ -79 +80 @@
 void Hermitian(TMatrix< TYPE >& A);
 r_8  Check_Mat_Ident(TMatrix< TYPE >& A);
+r_8  Check_Mat_Zero(TMatrix< TYPE >& A);
 r_8  Check_Mat_VecCol_0(TMatrix< TYPE >& A);
 void  Check_Mat_VecCol_2(TMatrix< complex<TYPER> >& A);
@@ -271 +273 @@
 cout<<"Compute AiA = A * InvA"<<endl;
 TMatrix< TYPE > AiA(N,N); AiA = Ainput * InvA;
-cout<<"------------ TMatrix AiA = A * InvA:"<<endl;
 if(nprline>0) {cout<<AiA; cout<<endl;}
 //-- Check
@@ -321 +322 @@
 cout<<"Compute AiA = A * InvA"<<endl;
 TMatrix< TYPE > AiA(N,N); AiA = Ainput * InvA;
-cout<<"------------ TMatrix AiA = A * InvA:"<<endl;
 if(nprline>0) {cout<<AiA; cout<<endl;} 
 //-- Check
@@ -426 +426 @@
 
 
+////////////////////////////////////////////////
+///////////// Test avec Lapack SVD /////////////
+////////////////////////////////////////////////
+#if defined(USE_LAPACK) && defined(ALSO_LAPACK_SVD)
+{
+cout<<"\n=========================================="<<endl;
+cout<<"------------ SVD decompositon LapackSVD "<<endl;
+A=Ainput;
+TVector< TYPE > S(N,N); TMatrix< TYPE > U(N,N); TMatrix< TYPE > VT(N,N); 
+//-- Decompositon
+int_4 info = LapackSVD(A,S,U,VT);
+cout<<"info="<<info<<endl;
+PrtTim("--- End of LapackSVD decompositon ---");
+if(nprline>0) {cout<<S; cout<<endl;}
+double smax = ABS_VAL(S(0)), smin = ABS_VAL(S(N-1));
+cout<<" Smin = |"<<S(N-1)<<"| = "<<smin<<", "
+    <<" Smax = |"<<S(0)<<"| = "<<smax<<", "
+    <<" --> Smin/Smax = "<<smin/smax<<endl;
+//-- A = U*S*VT
+cout<<"Compute A = U*S*VT"<<endl;
+TMatrix< TYPE > AmUSVt(N,N); AmUSVt = U;
+for(int i=0;i<N;i++) for(int j=0;j<N;j++) AmUSVt(i,j) *= S(j);
+AmUSVt *= VT; AmUSVt -= Ainput;
+if(nprline>0) {cout<<AmUSVt; cout<<endl;} 
+//-- Check
+Check_Mat_Zero(AmUSVt);
+PrtTim("--- End of LapackSVD Test ---");
+}
+#endif
+
+
 ///////////////////////////////////////////////
 ///////// Test Inversion avec GausPiv /////////
@@ -493 +524 @@
  cout<<"Ecart maximum par rapport a 0 hors diagonale = "<<vmaxh<<endl;
  return (vmaxd>vmaxh)? vmaxd: vmaxh;
+}
+
+////-------------------------------------------------------
+r_8  Check_Mat_Zero(TMatrix< TYPE >& A)
+// Compute the biggest difference element by element of A / Zero matrix
+{
+ int_4 N = A.NRows();
+ r_8 vmax = -1.;
+ for(int i=0;i<N;i++) for(int j=0;j<N;j++)
+   if( ABS_VAL(A(i,j)) > vmax ) vmax = ABS_VAL(A(i,j));
+ cout<<"Ecart maximum par rapport a zero = "<<vmax<<endl;
+ return vmax;
 }