svdrank.cpp

Go to the documentation of this file.
00001 //===========================================================================
00043 //===========================================================================
00044 
00045 #include <cmath>
00046 #include <SharkDefs.h>
00047 #include <LinAlg/LinAlg.h>
00048 
00049 #ifdef _WINDOWS
00050 // disable warning C2055: unreferenced formal parameter
00051 #pragma warning(disable: 4100)
00052 #endif
00053 
00054 //===========================================================================
00095 unsigned svdrank
00096 (
00097     const Array2D< double >& amatA,
00098     Array2D< double >& umatA,
00099     Array2D< double >& vmatA,
00100     Array  < double >& wvecA
00101 )
00102 {
00103     // unsigned m = amatA.rows( );
00104     unsigned n = amatA.cols();
00105     // double** const a = amatA.ptrArr( );
00106     // double **u = umatA.ptrArr( );
00107     // double **v = vmatA.ptrArr( );
00108     // double  *w = wvecA.elemvec();
00109 
00110     unsigned r;
00111     double   s, t;
00112 
00113     /*
00114       determine numerical rank: if the last singular values are < 0, the
00115           numerical error for the values > 0 is at least as large as the
00116       abslolute value of the last eigenvalue, thus values below this
00117       threshold are to be discarded. Additional thresholds are the relative
00118       machine accuracy and the relative error in the singular values    
00119      */
00120     for (r = 0; r < n && wvecA( r ) > 0.; r++);
00121 
00122     t = r < n ? fabs(wvecA( n-1 )) : 0.0;
00123     r = 0;
00124     s = 0.0;
00125     while (r < n     &&
00126             wvecA( r ) > t &&
00127             wvecA( r ) + s > s)
00128         s += wvecA( r++ );
00129 
00130     return r;
00131 }
00132 
00133 
00134 
00135 
00136 
00137