Changeset 975 in Sophya for trunk/SophyaProg

Timestamp:

Apr 27, 2000, 7:52:13 PM (25 years ago)

Author:

ansari

Message:

cm 27/4/00

File:

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

trunk/SophyaProg/Tests/tsttminv.cc

@@ -18 +18 @@
 
 #if defined(COMPLEX)
-#define ABS_VAL(_x_) sqrt((double)(_x_.real()*_x_.real() + _x_.imag()*_x_.imag()))
-#if defined(PRECIS32)
-#define TYPE complex<r_4>
+  #define ABS_VAL(_x_) sqrt((double)((_x_).real()*(_x_).real() + (_x_).imag()*(_x_).imag()))
+  #if defined(PRECIS32)
+    #define TYPE complex<r_4>
+  #else
+    #define TYPE complex<r_8>
+  #endif
 #else
-#define TYPE complex<r_8>
-#endif
-#else
-#define ABS_VAL(_x_) fabs((double)_x_)
-#if defined(PRECIS32)
-#define TYPE r_4
-#else
-#define TYPE r_8
-#endif
+  #define ABS_VAL(_x_) fabs((double)_x_)
+  #if defined(PRECIS32)
+    #define TYPE r_4
+  #else
+    #define TYPE r_8
+  #endif
 #endif
 
@@ -37 +37 @@
 //--------------------------------------------------------
 // number of lines/columns
-int_4 N = 5;
+uint_4 N = 5;
 // scale of the value (if =1 values between -1 and 1)
 r_8 scale = 1.;
 // number of values change by +/- vbig
-int_4 nbig = N;
+uint_4 nbig = N;
 r_8 vbig = 1000.;
-// Nombre de ligne de matrice a imprimer
- int_4 nprline = 10;
+// Nombre de lignes de matrice a imprimer
+uint_4 nprline = 10;
+// Initialisation du pauvre de l'aleatoire
+ uint_4 nalea = 0;
 //--------------------------------------------------------
 
 //-- Decodage arguments
 char c;
-while((c = getopt(narg,arg,"hv:l:")) != -1) {
+while((c = getopt(narg,arg,"hv:l:a:")) != -1) {
   switch (c) {
   case 'v' :
@@ -57 +59 @@
     sscanf(optarg,"%d",&nprline);
     break;
+  case 'a' :
+    sscanf(optarg,"%d",&nalea);
+    break;
   case 'h' :
-    cout<<"tsttminv [-h] [-v N,scale,nbig,vbig] [-l nprline]"<<endl;
+    cout<<"tsttminv [-h] [-v N,scale,nbig,vbig] [-l nprline] [-a nalea]"<<endl;
     cout<<"matrix filled with : {[-1,1] +/- vbig(nbig time)}*scale"<<endl;
     exit(-1);
@@ -64 +69 @@
 }
 if(N<=1) N = 1;
-cout<<"Taille matrice N = "<<N<<endl;
-cout<<"Elements entre +/- "<<scale<<endl;
-cout<<"Nombre de valeurs hors standard "<<nbig<<endl;
-cout<<"Valeurs hors standard v = v*(+/- "<<vbig<<" )"<<endl;
+cout<<"Taille matrice NxN, N = "<<N<<endl;
+cout<<"Elements entre +/- scale = "<<scale<<endl;
+cout<<"Nombre de valeurs hors standard nbig = "<<nbig<<endl;
+cout<<"Valeurs hors standard (+/- vbig = "<<vbig<<" ) * scale"<<endl;
 cout<<"Nombre de lignes de matrice a imprimer "<<nprline<<endl;
+cout<<"Initialisation de l aleatoire par "<<nalea<<" tirages"<<endl;
 cout<<endl;
 
@@ -80 +86 @@
 TMatrix< TYPE > B(N,N);
 TMatrix< TYPE > C(N,N);
+TVector< int >  Vind(N*N);
 A.Show();
 
 //-- Mise a zero
-A = (TYPE) 0;
+A    = (TYPE) 0;
 InvA = (TYPE) 0;
-AiA = (TYPE) 0;
-B = (TYPE) 0;
-C = (TYPE) 0;
+AiA  = (TYPE) 0;
+B    = (TYPE) 0;
+C    = (TYPE) 0;
 BaseArray::SetMaxPrint(nprline*N,0);
 
 //-- Fill matrices
 uint_8 k;
-int_4 i,j; double s;
+uint_4 i,j; double s;
+uint_4 nbr=0, nbi=0;
+Vind = 0;
+if(nalea>0) for(i=0;i<nalea;i++) drand01();
 A = Sequence(RandomSequence(RandomSequence::Flat,0.,1.));
-if(nbig>0) for(i=0;i<nbig;i++)
-  {k=(uint_8)(drand01()*N*N); if(k>=N*N) k--;
-   s=(drand01()>0.5)?1.:-1.; A[k] += (TYPE) s*vbig;}
+if(nbig>0) for(i=0;i<nbig;i++) {
+  k=(uint_8)(drand01()*N*N); if(k>=N*N) k--;
+  s=(drand01()>0.5)?1.:-1.;
+  if(Vind[k]==0) {A[k] += (TYPE) s*vbig; Vind[k]=1; nbr++;}
+}
 A *= (TYPE) scale;
 #if defined(COMPLEX)
+Vind = 0;
 B = Sequence(RandomSequence(RandomSequence::Flat,0.,1.));
-if(nbig>0) for(i=0;i<nbig;i++)
-  {k=(uint_8)(drand01()*N*N); if(k>=N*N) k--;
-   s=(drand01()>0.5)?1.:-1.; B[k] += (TYPE) s*vbig;}
+if(nbig>0) for(i=0;i<nbig;i++) {
+  k=(uint_8)(drand01()*N*N); if(k>=N*N) k--;
+  s=(drand01()>0.5)?1.:-1.;
+  if(Vind[k]==0) {B[k] += (TYPE) s*vbig; Vind[k]=1; nbi++;}
+}
 B *= (TYPE) scale;
 A += TYPE(0.,1.)*B;
 #endif
+cout<<"Nombre de valeurs BIG reelles = "<<nbr<<", imaginaires = "<<nbi<<endl;
 
 //-- Print matrice A
@@ -112 +128 @@
 
 //-- Inversion
-cout<<"Inversion"<<endl;
+cout<<"------------ Inversion"<<endl;
 InvA = Inverse(A);
 cout<<"------------ TMatrix InvA = A^(-1):"<<endl;
@@ -124 +140 @@
 
 //-- Check
-double vmin,vmax;
-k = 0;
-for(int i=0;i<N;i++) {
-  double absv = ABS_VAL( AiA(i,i) );
-  if(k==0) {vmin = vmax = absv; k++; continue;}
+double vmax=-1.;
+for(i=0;i<N;i++) {
+  double absv = ABS_VAL( 1. - AiA(i,i) );
   if( absv > vmax ) vmax = absv;
-  if( absv < vmin ) vmin = absv;
-  k++;
 }
-cout<<"Limites sur la diagonale  "
-    <<vmin<<" , "<<vmax<<"  n="<<k<<endl;
-k = 0;
-for(int i=0;i<N;i++) for(int j=0;j<N;j++) {
+cout<<"Ecart maximum par rapport a 1 sur diagonale = "<<vmax<<endl;
+vmax = -1.;
+for(i=0;i<N;i++) for(j=0;j<N;j++) {
   if( i == j ) continue;
   double absv = ABS_VAL( AiA(i,j) );
-  if(k==0) {vmin = vmax = absv; k++; continue;}
   if( absv > vmax ) vmax = absv;
-  if( absv < vmin ) vmin = absv;
-  k++;
 }
-cout<<"Limites hors diagonale  "
-    <<vmin<<" , "<<vmax<<"  n="<<k<<endl;
+cout<<"Ecart maximum par rapport a 0 hors diagonale = "<<vmax<<endl;
 
 exit(0);