int pinv const DblNumMat M,
double  epsilon,
DblNumMat R
 

Definition at line 126 of file vecmatop.cpp.

References iC, NumMat< F >::m(), and NumMat< F >::n().

00127 {
00128   assert(M.m() == R.n());  assert(M.n() == R.m());
00129   SVDRep svd;
00130   iC( svd.construct(epsilon, M) );
00131   //invert Svd    //cerr<<svd.S()(0)<<" "<<svd.S()(svd.S().m()-1)<<" "<<svd.S()(0)/svd.S()(svd.S().m()-1)<<endl;
00132   double cutoff = epsilon * svd.S()(0);  //double cutoff = epsilon;
00133   for(int i=0; i<svd.S().m(); i++) {
00134     if( svd.S()(i) >= cutoff) {
00135       svd.S()(i) = 1.0/(svd.S()(i));
00136          } else {
00137                 assert(0);      //svd.S()(i) = 0.0;
00138          }
00139   }
00140   DblNumMat UT(svd.U().n(),  svd.U().m());
00141   DblNumMat V( svd.VT().n(), svd.VT().m());
00142   iC( tran(svd.U(), UT) );
00143   iC( tran(svd.VT(), V) );
00144   for(int i=0; i<V.m(); i++)
00145     for(int j=0; j<V.n(); j++) {
00146       V(i,j) = V(i,j) * svd.S()(j);
00147          }
00148   
00149   char transa = 'N';
00150   char transb = 'N';
00151   double alpha = 1.0;
00152   double beta = 0.0;
00153   int m = V.m();
00154   int n = UT.n();
00155   int k = V.n();
00156   DGEMM(&transa, &transb, &m, &n, &k, &alpha,
00157                   V.data(), &m, UT.data(), &k, 
00158                   &beta, R.data(), &m);  
00159   
00160   return 0;
00161 }


Generated on Sun Dec 4 18:13:13 2005 for common by  doxygen 1.4.5