svdrep.cpp

Go to the documentation of this file.
00001 /* Kernel Independent Fast Multipole Method
00002    Copyright (C) 2004 Lexing Ying, New York University
00003 
00004 This program is free software; you can redistribute it and/or modify
00005 it under the terms of the GNU General Public License as published by
00006 the Free Software Foundation; either version 2, or (at your option)
00007 any later version.
00008 
00009 This program is distributed in the hope that it will be useful, but WITHOUT
00010 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00011 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00012 for more details.
00013 
00014 You should have received a copy of the GNU General Public License
00015 along with this program; see the file COPYING.  If not, write to the Free
00016 Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
00017 02111-1307, USA.  */
00019 #include "blas.h"
00020 #include "lapack.h"
00021 #include "svdrep.hpp"
00022 
00023 using std::min;
00024 using std::max;
00025 using std::abs;
00026 
00027 //int    SVDRep::_wssize = 4194304;
00028 //double SVDRep::_wsbuf[4194304];
00029 
00037 /* ********************************************************************** */
00038 int SVDRep::construct(double epsilon, const DblNumMat& K)
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 }
00081 
00082 /* ********************************************************************** */
00083 int SVDRep::dgemv(double alpha, const DblNumVec& X, double beta, DblNumVec& Y, double tol)
00084 {
00085   iA(Y.m() == _matU.m());
00086   iA(X.m() == _matVT.n());
00087   iC( dgemv(alpha, X.data(), beta, Y.data(), tol) );
00088   
00089   return 0;
00090 }
00091 
00092 /* ********************************************************************** */
00093 int SVDRep::dgemv(double alpha, double* X, double beta, double* Y, double tol)
00094 {
00095   int K = 1; //prevent matrix of zero size
00096   while(K<_matS.m() && _matS(K)>=tol) 
00097          K++;
00098   //buf = VT(1:K,:) * X
00099   //double* buf = _wsbuf;
00100   double* buf = new double[K];
00101   {
00102          char TRANS = 'N';
00103          int M = _matVT.m();
00104          int N = _matVT.n();
00105          double ALPHA = 1.0;
00106          double BETA = 0.0;
00107          int INC = 1;
00108          DGEMV(&TRANS, &K, &N, &ALPHA, _matVT.data(), &M, X, &INC, &BETA, buf, &INC); //first K rows
00109   }
00110   // buf = S(1:K) .* buf;
00111   for(int i=0; i<K; i++)
00112          buf[i] = buf[i] * _matS(i);
00113   // y = U(:,1:K) * buf
00114   {
00115          char TRANS = 'N';
00116          int M = _matU.m(); //int N = _matU.n();
00117          double ALPHA = alpha;
00118          double BETA = beta;
00119          int INC = 1;
00120          DGEMV(&TRANS, &M, &K, &ALPHA, _matU.data(), &M, buf, &INC, &BETA, Y, &INC);    
00121   }
00122   delete [] buf;
00123   
00124   return 0;
00125 }

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