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

Timestamp:

Jul 27, 2004, 9:58:18 AM (21 years ago)

Author:

cmv

Message:

add test inversion with LeastSquareSVD_DC cmv 27/7/04

File:

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

trunk/SophyaProg/Tests/tsttminv.cc

@@ -24 +24 @@
 #define ALSO_LAPACK_INV_SYM
 #define ALSO_LAPACK_INV_LSS
+#define ALSO_LAPACK_INV_LSS_SVD
 #define ALSO_LAPACK_EV
 #define ALSO_LAPACK_EV_SYM
@@ -264 +265 @@
 {
 cout<<"\n=========================================="<<endl;
-cout<<"------------ Inversion LAPACK"<<endl;
+cout<<"------------ Inversion LAPACK (LU factorization)"<<endl;
 A = Ainput;
 //-- Inversion
@@ -288 +289 @@
 {
 cout<<"\n=========================================="<<endl;
-cout<<"------------ Inversion LAPACK sym"<<endl;
+cout<<"------------ Inversion LAPACK symetric (LU factorization)"<<endl;
 TMatrix< TYPE > Asym(N,N); Asym=Ainput; Symetrize(Asym); A=Asym;
 //-- Inversion
@@ -313 +314 @@
 {
 cout<<"\n=========================================="<<endl;
-cout<<"------------ Inversion LAPACK LeastSquare"<<endl;
+cout<<"------------ Inversion LAPACK LeastSquare (QR or LQ factorization)"<<endl;
 A = Ainput;
 //-- Inversion
@@ -327 +328 @@
 Check_Mat_Ident(AiA);
 PrtTim("--- End of LapackLeastSquareSolve Test ---");
+}
+#endif
+
+
+////////////////////////////////////////////////
+///////// Test avec Lapack LeastSquare /////////
+////////////////////////////////////////////////
+#if defined(USE_LAPACK) && defined(ALSO_LAPACK_INV_LSS_SVD)
+{
+cout<<"\n=========================================="<<endl;
+cout<<"------------ Inversion LAPACK LeastSquare (SVD decomposition)"<<endl;
+A = Ainput;
+//-- Inversion
+TMatrix< TYPE > InvA(N,N); InvA = IdentityMatrix(1.,N);
+TVector<r_8> S;
+int_4 rank;
+r_8 rcond = -1.;
+int_4 info = LapackLeastSquareSolveSVD_DC(A,InvA,S,rank,rcond);
+cout<<"info="<<info<<"  (rank="<<rank<<")"<<endl;
+PrtTim("--- End of LapackLeastSquareSolveSVD_DC Inversion ---");
+if(nprline>0) {cout<<S; cout<<endl;}
+double smax = fabs(S(0)), smin = fabs(S(S.Size()-1));
+cout<<" Smin = |"<<S(S.Size()-1)<<"| = "<<smin<<", "
+    <<" Smax = |"<<S(0)<<"| = "<<smax<<", "
+    <<" --> Smin/Smax = "<<smin/smax<<endl;
+//-- AiA = A * InvA
+cout<<"Compute AiA = A * InvA"<<endl;
+TMatrix< TYPE > AiA(N,N); AiA = Ainput * InvA;
+if(nprline>0) {cout<<AiA; cout<<endl;} 
+//-- Check
+Check_Mat_Ident(AiA);
+PrtTim("--- End of LapackLeastSquareSolveSVD_DC Test ---");
 }
 #endif
@@ -466 +499 @@
 cout<<"------------ SVD decompositon LapackSVD_DC "<<endl;
 A=Ainput;
-TVector< TYPE > S; TMatrix< TYPE > U; TMatrix< TYPE > VT; 
+TVector< r_8 > S; TMatrix< TYPE > U; TMatrix< TYPE > VT; 
 //-- Decompositon
 int_4 info = LapackSVD_DC(A,S,U,VT);
@@ -472 +505 @@
 PrtTim("--- End of LapackSVD_DC decompositon ---");
 if(nprline>0) {cout<<S; cout<<endl;}
-double smax = ABS_VAL(S(0)), smin = ABS_VAL(S(N-1));
+double smax = fabs(S(0)), smin = fabs(S(N-1));
 cout<<" Smin = |"<<S(N-1)<<"| = "<<smin<<", "
     <<" Smax = |"<<S(0)<<"| = "<<smax<<", "
@@ -479 +512 @@
 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);
+for(int i=0;i<N;i++) for(int j=0;j<N;j++) AmUSVt(i,j) *= (TYPE) S(j);
 AmUSVt *= VT; AmUSVt -= Ainput;
 if(nprline>0) {cout<<AmUSVt; cout<<endl;}