Changeset 2647 in Sophya for trunk/SophyaProg/Tests/lpk.cc


Ignore:
Timestamp:
Feb 7, 2005, 5:43:49 PM (21 years ago)
Author:
ansari
Message:

Ajout programme test inversion matrice par SimpleMatrixOperation et ajout test inversion par lapack ds lpk.cc - Reza 7 Fev 2005

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaProg/Tests/lpk.cc

    r2615 r2647  
    1212
    1313
    14 int lpk_rzt_arr(int n);
    15 int lpk_linsolve_mtx(int n);
    16 int lpk_svd_mtx(int l, int c, int wsf, bool covu=true);
    17 int lpk_leastsquare_mtx(int l, int c);
     14int lpk_linsolve_mtx(int n);                             // ErrorCode = 8
     15int lpk_svd_mtx(int l, int c, int wsf, bool covu=true);  // ErrorCode = 16
     16int lpk_leastsquare_mtx(int l, int c);                   // ErrorCode = 32   
     17int lpk_inverse_mtx(int l);                              // ErrorCode = 64
    1818
    1919static double TOLERANCE = 1.e-6;
     20static int nprt = 100;
     21static int prtlev = 1;
    2022
    2123int main(int narg, char* arg[])
    … …  
    2931  if (narg < 2) {
    3032    cout << " lpk - LinAlg/LapackServer test \n"
    31          << " Usage: lpk linsolve/svd/rzt [sizeL,C=5,5] [prtlev=0] \n"
     33         << " Usage: lpk select [sizeL,C=5,5] [prtlev=1] \n"
    3234         << "             [nprtmax=100] [WorkSpaceSizeFactor=2] \n"
     35         << "   select= linsolve svd svds lss inverse / all= \n"
    3336         << "   linsolve:  lpk_linsolve_mtx() LapackServer::LinSolve with TMatrix<r_8> \n"
    3437         << "   svd:  lpk_svd_mtx() LapackServer::SVD(a,s,u,vt) with TMatrix<r_8> \n"
    3538         << "   svds:  lpk_svd_mtx() LapackServer::SVD(a,s); with TMatrix<r_8> \n"
    3639         << "   lss:  lpk_leastsquare_mtx() LapackServer::LeastSquareSolve() r_8 \n"
    37          << "   rzt:  lpk_rzt_tarr() rztest_lapack with TArray<r_4> \n" << endl;
     40         << "   inverse:  LapackServer::ComputeInverse() r_8 \n" << endl;
    3841    exit(0);
    3942  }
    … …  
    4346  string opt = arg[1];
    4447  if (narg > 2) sscanf(arg[2], "%d,%d", &l, &c);
    45   int nprt = 100;
    46   int prtlev = 1;
    4748  if (narg > 3) prtlev = atoi(arg[3]);
    4849  if (narg > 4) nprt = atoi(arg[4]);
    4950  if (narg > 5) wsf = atoi(arg[5]);
    5051
     52  cout << " lpk - LinAlg/LapackServer_Test sizeL,C=" << l << "," << c
     53       << " opt= " << opt << endl;
    5154  int rc = 0;
    5255  BaseArray::SetMaxPrint(nprt, prtlev);
    5356  try {
    54     if (opt == "linsolve") rc = lpk_linsolve_mtx(l);
     57    if (opt == "all") {
     58      rc += lpk_linsolve_mtx(l);
     59      rc += lpk_svd_mtx(l,c,wsf,true);
     60      rc += lpk_leastsquare_mtx(l,c);
     61      rc += lpk_inverse_mtx(l);
     62    }
     63    else if (opt == "linsolve") rc = lpk_linsolve_mtx(l);
    5564    else if (opt == "svd") rc = lpk_svd_mtx(l,c,wsf,true);
    5665    else if (opt == "svds") rc = lpk_svd_mtx(l,c,wsf,false);
    5766    else if (opt == "lss") rc = lpk_leastsquare_mtx(l,c);
    58     else if (opt == "rzt") rc = lpk_rzt_arr(l);
     67    else if (opt == "inverse") rc = lpk_inverse_mtx(l);
    5968    else { cout << " Unknown option " << opt << " ! " << endl; rc = 66; }
    6069  }
    … …  
    8796  Vector x(n), b;
    8897  //  Matrix  x(n,1), b;
    89   cout << " ------------ Vector X = \n " << x << "\n" << endl;
     98  if (prtlev > 0)
     99    cout << " ------------ Vector X = \n " << x << "\n" << endl;
    90100  for(i=0; i<n; i++) x(i) = GauRnd(2., 1.5);
    91101  b = a*x;
    92102
    93   cout << " ---- lpk_tmtx() LapackServer::LinSolve Test Using TMatrix<r_8> ----- " << endl;
    94   cout << " ------------ Matrix A = \n " << a << "\n" << endl;
    95   cout << " ------------ Matrix X = \n " << x << "\n" << endl;
    96   cout << " ------------ Matrix B = \n " << b << "\n" << endl;
    97 
     103  if (prtlev > 0) {
     104    cout << " ---- lpk_tmtx() LapackServer::LinSolve Test Using TMatrix<r_8> ----- " << endl;
     105    cout << " ------------ Matrix A = \n " << a << "\n" << endl;
     106    cout << " ------------ Matrix X = \n " << x << "\n" << endl;
     107    cout << " ------------ Matrix B = \n " << b << "\n" << endl;
     108  }
    98109  cout << "\n   Calling LapackLinSolve(a,b) .... " << endl;
    99110  PrtTim(" Calling LapackLinSolve(a,b) ");
    … …  
    101112  PrtTim(" End LapackLinSolve(a,b) ");
    102113
    103   cout << " ------------ Result B(=X ?) = \n " << b << "\n" << endl;
     114  if (prtlev > 0)
     115    cout << " ------------ Result B(=X ?) = \n " << b << "\n" << endl;
    104116  Vector diff = b-x;
    105117  PrtTim(" End of Compute(diff)");
    106   cout << " ------------ Vector diff B-X = \n " << diff << "\n" << endl;
     118  if (prtlev > 0)
     119    cout << " ------------ Vector diff B-X = \n " << diff << "\n" << endl;
    107120  double min,max;
    108121  diff.MinMax(min, max);
    … …  
    112125    cout << " !!! Difference exceeding tolerance (=" << TOLERANCE << ") !!!"
    113126         << endl;
    114     return(99);
     127    return(8);
    115128  }
    116129  return(0);
    … …  
    131144  Matrix u, vt;
    132145  cout << " ---- lpk_svd_tmtx() LapackServer::SVD Test Using TMatrix<r_8> ---- " << endl;
    133   cout << " ------------ Matrix A = \n " << a << "\n" << endl;
     146  if (prtlev > 0)
     147    cout << " ------------ Matrix A = \n " << a << "\n" << endl;
    134148
    135149  cout << "\n   Calling LapackSVD(a,s,u,vt) .... " << endl;
    … …  
    141155  PrtTim(" End LapackSVD() ");
    142156
    143   cout << " ------------ Result S  = \n " << s << "\n" << endl;
     157  if (prtlev > 0)
     158    cout << " ------------ Result S  = \n " << s << "\n" << endl;
    144159  if (!covu) return(0);
    145   cout << " ------------ Result U  = \n " << u << "\n" << endl;
    146   cout << " ------------ Result VT = \n " << vt << "\n" << endl;
     160  if (prtlev > 0) {
     161    cout << " ------------ Result U  = \n " << u << "\n" << endl;
     162    cout << " ------------ Result VT = \n " << vt << "\n" << endl;
     163  }
    147164  Matrix sm(m,n);
    148165  int minmn = (m<n) ? m : n ;
    … …  
    150167  Matrix diff = u*(sm*vt) - a;
    151168  PrtTim(" End of Compute(diff)");
    152   cout << " ------------ Matrix diff U*S*Vt - A = \n " << diff << "\n" << endl;
     169  if (prtlev > 0)
     170    cout << " ------------ Matrix diff U*S*Vt - A = \n " << diff << "\n" << endl;
    153171  double min,max;
    154172  diff.MinMax(min, max);
    … …  
    158176    cout << " !!! Difference exceeding tolerance (=" << TOLERANCE << ") !!!"
    159177         << endl;
    160     return(99);
     178    return(16);
    161179  }
    162180  return(0);
    … …  
    182200 
    183201  cout << " ---- lpk_leastsquare_tmtx() LapackLeastSquareSolve<r_8> ---- " << endl;
    184   cout << " ------------ Matrix A = \n " << a << "\n" << endl;
    185   cout << " ------------ Vector B = \n " << bx << "\n" << endl;
    186   cout << " ------------ Vector X = \n " << x << "\n" << endl;
    187 
     202  if (prtlev > 0) {
     203    cout << " ------------ Matrix A = \n " << a << "\n" << endl;
     204    cout << " ------------ Vector B = \n " << bx << "\n" << endl;
     205    cout << " ------------ Vector X = \n " << x << "\n" << endl;
     206  }
    188207  cout << "\n   Calling  LapackLeastSquareSolve(a,bx) .... " << endl;
    189208  PrtTim(" Calling LapackLeastSquareSolve() ");
    … …  
    191210  PrtTim(" End LapackLeastSquareSolve() ");
    192211  bx.Share(bx(Range(0,0,n)));
    193   cout << " ------------ Result X  = \n " << bx << "\n" << endl;
    194   cout << " ------------ X-X_Real  = \n " << bx-x << "\n" << endl;
     212  if (prtlev > 0) {
     213    cout << " ------------ Result X  = \n " << bx << "\n" << endl;
     214    cout << " ------------ X-X_Real  = \n " << bx-x << "\n" << endl;
     215  }
    195216  Vector diff = b-a*bx;
    196217  PrtTim(" End of Compute(diff)");
    197   cout << " ------------ Matrix diff b-a*x = \n " << diff << "\n" << endl;
     218  if (prtlev > 0)
     219    cout << " ------------ Matrix diff b-a*x = \n " << diff << "\n" << endl;
    198220  double mean,sigma;
    199221  MeanSigma(diff, mean, sigma);
    … …  
    203225    cout << " !!! Difference exceeding tolerance (=" << TOLERANCE << ") !!!"
    204226         << endl;
    205     return(99);
    206   }
    207   return(0);
    208 
    209 }
    210 
    211 
    212 void rztest_lapack(TArray<r_4>& aa, TArray<r_4>& bb);
    213 
    214 int  lpk_rzt_arr(int n)
    215 {
    216   int i,j;
    217   TArray<r_4> a(n,n);
    218   for(i=0; i<n; i++)
    219     for(j=0; j<n; j++)  a(i,j,0) = GauRnd(0., 1.);
    220  
    221   TArray<r_4> x(n,1), b(n,1);
    222   r_4 sum ;
    223   for(i=0; i<n; i++) x(i,0,0) = GauRnd(2., 1.5);
    224   for(i=0; i<n; i++) {
    225     sum = 0.;
    226     for(j=0; j<n; j++) sum +=  a(i,j,0)*x(j,0,0);
    227     b(i,0,0) = sum;
    228   }
    229 
    230   //  cout << ":::::::: rztest_lapack - Size=" << n << "  ::::::::: " << endl;
    231   cout << " ------- lpk_tarr() LapaackTest Using TArray<r_4> -------- " << endl;
    232   cout << " ------------ Array A = \n " << a << "\n" << endl;
    233   cout << " ------------ Array X = \n " << x << "\n" << endl;
    234   cout << " ------------ Array B = \n " << b << "\n" << endl;
    235 
    236   cout << "\n     Calling rztest_lapack ... " << endl;
    237 
    238   rztest_lapack(a, b);
    239 
    240   cout << " ------------ Result B(=X ?) = \n " << b << "\n" << endl;
    241   return(0);
    242 }
     227    return 32;
     228  }
     229  return 0;
     230
     231}
     232// -----------------------------------------------------------------------
     233/* Nouvelle-Fonction */
     234int lpk_inverse_mtx(int m)
     235{
     236  BaseArray::SetDefaultMemoryMapping(BaseArray::FortranMemoryMapping);
     237  cout << " lpk_inverse_mtx() - Test of LapackComputeInverse<r_8> " << endl;
     238 
     239  Matrix a(m , m), aa, ainv, idmx(m,m), mxprod, diff;
     240  a = RandomSequence(RandomSequence::Gaussian, 0., 4.);
     241  aa = a;  // We make a copy of a, as it is modified by LapackServer
     242  PrtTim(" End of a_inint");
     243
     244  cout << "\n   Calling  LapackInverse() .... " << endl;
     245  ainv = LapackInverse(aa);
     246  PrtTim(" End of LapackInverse(a)");
     247  mxprod = a*ainv;
     248  PrtTim(" End of mxprod = a*ainv");
     249  idmx = IdentityMatrix();
     250  PrtTim(" End of idmx = IdentityMatrix()");
     251  diff = mxprod-idmx;
     252  PrtTim(" End of Compute(diff)");
     253
     254  if (prtlev > 0) {
     255    cout << " ------------ Matrix A = \n " << a << "\n" << endl;
     256    cout << " ------------ Matrix Inverse(A) = \n " << ainv << "\n" << endl;
     257    cout << " ------------ Matrix mxprod = A*Inverse(A) = \n " << mxprod << "\n" << endl;
     258  }
     259  double min,max;
     260  diff.MinMax(min, max);
     261  PrtTim(" End of diffCheck");
     262  cout << " Min/Max difference Matrix (?=0) , Min= " << min
     263       << " Max= " << max << endl;
     264  if ((fabs(min) > TOLERANCE) || (fabs(max) > TOLERANCE)) {
     265    cout << " !!! Difference exceeding tolerance (=" << TOLERANCE << ") !!!"
     266         << endl;
     267    return 64;
     268  }
     269  return 0;
     270
     271}
     272
Note: See TracChangeset for help on using the changeset viewer.