int SVDRep::construct double  epsilon,
const DblNumMat M
 

Construct the SVD representation of a matrix M with cutoff epsilon. The smallest singular values are dropped based on epsilon. The Lapack routine DGESVD from lapack.h is called here as well

Definition at line 38 of file svdrep.cpp.

References _matS, _matU, _matVT, abs(), NumVec< F >::data(), NumMat< F >::data(), DGESVD, iA, k(), NumMat< F >::m(), m(), max(), min(), NumMat< F >::n(), n(), NumVec< F >::resize(), and NumMat< F >::resize().

00039 {
00040   int m = K.m();
00041   int n = K.n();
00042   int k = min(m, n);
00043   
00044   DblNumMat tU(m, k);
00045   DblNumVec tS(k);
00046   DblNumMat tVT(k, n);
00047   
00048   //SVD
00049   int INFO;
00050   char JOBU  = 'S';
00051   char JOBVT = 'S';
00052   
00053   int wssize = max(3*min(m,n)+max(m,n), 5*min(m,n));
00054   double* wsbuf = new double[wssize];
00055   DGESVD(&JOBU, &JOBVT, &m, &n, K.data(), &m, tS.data(), tU.data(), &m, tVT.data(), &k, wsbuf, &wssize, &INFO);  iA(INFO==0);
00056   delete [] wsbuf;
00057   
00058   //cutoff
00059   double cutoff = epsilon*tS(0);
00060   int cnt=0;
00061   while(cnt< k)
00062     if(abs(tS(cnt)) >= cutoff)
00063       cnt++;
00064     else
00065       break;
00066   _matU.resize(m, cnt);
00067   _matS.resize(cnt);    
00068   _matVT.resize(cnt,n);
00069   
00070   for(int i=0; i<m; i++)
00071     for(int j=0; j<cnt; j++)
00072       _matU(i,j) = tU(i,j);
00073   for(int i=0; i<cnt; i++)
00074     _matS(i) = tS(i);
00075   for(int i=0; i<cnt; i++)
00076     for(int j=0; j<n; j++)
00077       _matVT(i,j) = tVT(i,j);
00078   
00079   return 0;
00080 }


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