Changeset 1567 in Sophya

Timestamp:

Jul 5, 2001, 5:46:30 PM (24 years ago)

Author:

ansari

Message:

Ajout test de LapackServer<T>::LeastSquareSolve() ds lpk.cc - Reza 5/7/2001

File:

• trunk/SophyaProg/Tests/lpk.cc (modified) (4 diffs)

trunk/SophyaProg/Tests/lpk.cc

@@ -14 +14 @@
 int lpk_linsolve_mtx(int n);
 int lpk_svd_mtx(int l, int c, int wsf, bool covu=true);
+int lpk_leastsquare_mtx(int l, int c);
 
 static double TOLERANCE = 1.e-6;
@@ -31 +32 @@
          << "   linsolve:  lpk_linsolve_mtx() LapackServer::LinSolve with TMatrix<r_8> \n" 
          << "   svd:  lpk_svd_mtx() LapackServer::SVD(a,s,u,vt) with TMatrix<r_8> \n" 
-         << "   svds:  lpk_svd_mtx() LapackServer::SVD(a,s); with TMatrix<r_8> \n" 
+         << "   svds:  lpk_svd_mtx() LapackServer::SVD(a,s); with TMatrix<r_8> \n"
+         << "   lss:  lpk_leastsquare_mtx() LapackServer::LeastSquareSolve() r_8 \n"
          << "   rzt:  lpk_rzt_tarr() rztest_lapack with TArray<r_4> \n" << endl;
     exit(0);
@@ -52 +54 @@
     else if (opt == "svd") rc = lpk_svd_mtx(l,c,wsf,true);
     else if (opt == "svds") rc = lpk_svd_mtx(l,c,wsf,false);
+    else if (opt == "lss") rc = lpk_leastsquare_mtx(l,c);
     else if (opt == "rzt") rc = lpk_rzt_arr(l);
     else { cout << " Unknown option " << opt << " ! " << endl; rc = 66; }
@@ -160 +163 @@
 }
 
+// -----------------------------------------------------------------------
+/* Nouvelle-Fonction */
+int lpk_leastsquare_mtx(int m, int n)
+{
+  BaseArray::SetDefaultMemoryMapping(BaseArray::FortranMemoryMapping);
+  cout << " lpk_leastsquare_mtx() - Test of LapackLeastSquareSolve " << endl;
+  
+  Matrix a(m , n), aa;
+  a = RandomSequence(RandomSequence::Gaussian, 0., 4.);
+  aa = a;
+  Vector x(n),noise(m),b,bx;
+  x = RandomSequence(RandomSequence::Gaussian, 2., 1.);
+  double signoise = 0.1;
+  noise = RandomSequence(RandomSequence::Gaussian, 0., 0.1);
+  b = a*x;
+  bx = b+noise;
+ 
+  cout << " ---- lpk_leastsquare_tmtx() LapackLeastSquareSolve<r_8> ---- " << endl; 
+  cout << " ------------ Matrix A = \n " << a << "\n" << endl; 
+  cout << " ------------ Vector B = \n " << bx << "\n" << endl; 
+  cout << " ------------ Vector X = \n " << x << "\n" << endl; 
+
+  cout << "\n   Calling  LapackLeastSquareSolve(a,bx) .... " << endl;
+  PrtTim(" Calling LapackLeastSquareSolve() ");
+  LapackLeastSquareSolve(aa, bx);
+  PrtTim(" End LapackLeastSquareSolve() ");
+  bx.Share(bx(Range(0,0,n)));
+  cout << " ------------ Result X  = \n " << bx << "\n" << endl; 
+  cout << " ------------ X-X_Real  = \n " << bx-x << "\n" << endl; 
+  Vector diff = b-a*bx;
+  PrtTim(" End of Compute(diff)");
+  cout << " ------------ Matrix diff b-a*x = \n " << diff << "\n" << endl; 
+  double mean,sigma;
+  MeanSigma(diff, mean, sigma);
+  cout << " MeanSigma(diff) , Mean= " << mean << " Sigma=" << sigma
+       << " SigNoise= " << signoise << endl;
+  if ((fabs(mean) > signoise) || (fabs(sigma-signoise) > signoise)) {
+    cout << " !!! Difference exceeding tolerance (=" << TOLERANCE << ") !!!" 
+         << endl;
+    return(99);
+  } 
+  return(0);
+
+}
+
+
 void rztest_lapack(TArray<r_4>& aa, TArray<r_4>& bb);