core/vnl/algo/vnl_sparse_symmetric_eigensystem.cxx
Go to the documentation of this file.
00001 // This is core/vnl/algo/vnl_sparse_symmetric_eigensystem.cxx
00002 #ifdef VCL_NEEDS_PRAGMA_INTERFACE
00003 #pragma implementation
00004 #endif
00005 //:
00006 // \file
00007 
00008 #include "vnl_sparse_symmetric_eigensystem.h"
00009 #include "vnl_sparse_lu.h"
00010 #include <vnl/vnl_vector_ref.h>
00011 #include <vcl_cassert.h>
00012 #include <vcl_cstring.h>
00013 #include <vcl_cstdlib.h>
00014 #include <vcl_iostream.h>
00015 #include <vcl_vector.h>
00016 
00017 #include <vnl/algo/vnl_netlib.h> // dnlaso_() dseupd_() dsaupd_()
00018 
00019 static vnl_sparse_symmetric_eigensystem * current_system = 0;
00020 
00021 #ifdef VCL_SUNPRO_CC
00022 # define FUNCTION extern "C"
00023 #else
00024 # define FUNCTION static
00025 #endif
00026 
00027 //------------------------------------------------------------
00028 //: Callback for multiplying our matrix by a number of vectors.
00029 //  The input is p, which is an NxM matrix.
00030 //  This function returns q = A p, where A is the current sparse matrix.
00031 FUNCTION
00032 void sse_op_callback(const long* n,
00033                      const long* m,
00034                      const double* p,
00035                      double* q)
00036 {
00037   assert(current_system != 0);
00038 
00039   current_system->CalculateProduct(*n,*m,p,q);
00040 }
00041 
00042 //------------------------------------------------------------
00043 //: Callback for saving the Lanczos vectors as required by dnlaso.
00044 // If k=0, save the m columns of q as the (j-m+1)th through jth
00045 // vectors.  If k=1 then return the (j-m+1)th through jth vectors in
00046 // q.
00047 FUNCTION
00048 void sse_iovect_callback(const long* n,
00049                          const long* m,
00050                          double* q,
00051                          const long* j,
00052                          const long* k)
00053 {
00054   assert(current_system != 0);
00055 
00056   if (*k==0)
00057     current_system->SaveVectors(*n,*m,q,*j-*m);
00058   else if (*k==1)
00059     current_system->RestoreVectors(*n,*m,q,*j-*m);
00060 }
00061 
00062 vnl_sparse_symmetric_eigensystem::vnl_sparse_symmetric_eigensystem()
00063   : nvalues(0), vectors(0), values(0)
00064 {
00065 }
00066 
00067 vnl_sparse_symmetric_eigensystem::~vnl_sparse_symmetric_eigensystem()
00068 {
00069   delete[] vectors; vectors = 0;
00070   delete[] values; values = 0;
00071   for (unsigned i=0; i<temp_store.size(); ++i)
00072     delete temp_store[i];
00073   temp_store.clear();
00074 }
00075 
00076 //------------------------------------------------------------
00077 //: Here is where the fortran converted code gets called.
00078 // The sparse matrix M is assumed to be symmetric.  The n smallest
00079 // eigenvalues and their corresponding eigenvectors are calculated if
00080 // smallest is true (the default).  Otherwise the n largest eigenpairs
00081 // are found.  The accuracy of the eigenvalues is to nfigures decimal
00082 // digits.  Returns 0 if successful, non-zero otherwise.
00083 int vnl_sparse_symmetric_eigensystem::CalculateNPairs(vnl_sparse_matrix<double>& M,
00084                                                       int n,
00085                                                       bool smallest,
00086                                                       long nfigures)
00087 {
00088   mat = &M;
00089 
00090   // Clear current vectors.
00091   if (vectors) {
00092     delete[] vectors; vectors = 0;
00093     delete[] values; values = 0;
00094   }
00095   nvalues = 0;
00096 
00097   current_system = this;
00098 
00099   long dim = mat->columns();
00100   long nvals = (smallest)?-n:n;
00101   long nperm = 0;
00102   long nmval = n;
00103   long nmvec = dim;
00104   vcl_vector<double> temp_vals(n*4);
00105   vcl_vector<double> temp_vecs(n*dim);
00106 
00107   // set nblock = vcl_max(10, dim/6) :
00108   long nblock = (dim<60) ? dim/6 : 10;
00109 
00110   // isn't this rather a lot ? -- fsm
00111   long maxop = dim*10;      // dim*20;
00112 
00113   // set maxj = vcl_max(40, maxop*nblock, 6*nblock+1) :
00114   long maxj = maxop*nblock; // 2*n+1;
00115   long t1 = 6*nblock+1;
00116   if (maxj < t1) maxj = t1;
00117   if (maxj < 40) maxj = 40;
00118 
00119   // Calculate size of workspace needed.  These expressions come from
00120   // the LASO documentation.
00121   int work_size = dim*nblock;
00122   int t2 = maxj*(2*nblock+3) + 2*n + 6 + (2*nblock+2)*(nblock+1);
00123   if (work_size < t2) work_size = t2;
00124   work_size += 2*dim*nblock + maxj*(nblock + n + 2) + 2*nblock*nblock + 3*n;
00125   vcl_vector<double> work(work_size+10);
00126 
00127   // Set starting vectors to zero.
00128   for (int i=0; i<dim*nblock; ++i)
00129     work[i] = 0.0;
00130 
00131   vcl_vector<long> ind(n);
00132 
00133   long ierr = 0;
00134 
00135   v3p_netlib_dnlaso_(sse_op_callback, sse_iovect_callback,
00136                      &dim, &nvals, &nfigures, &nperm,
00137                      &nmval, &temp_vals[0],
00138                      &nmvec, &temp_vecs[0],
00139                      &nblock,
00140                      &maxop,
00141                      &maxj,
00142                      &work[0],
00143                      &ind[0],
00144                      &ierr);
00145   if (ierr > 0)
00146   {
00147     if (ierr & 0x1)
00148       vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem: N < 6*NBLOCK\n";
00149     if (ierr & 0x2)
00150       vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem: NFIG < 0\n";
00151     if (ierr & 0x4)
00152       vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem: NMVEC < N\n";
00153     if (ierr & 0x8)
00154       vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem: NPERM < 0\n";
00155     if (ierr & 0x10)
00156       vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem: MAXJ < 6*NBLOCK\n";
00157     if (ierr & 0x20)
00158       vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem: NVAL < max(1,NPERM)\n";
00159     if (ierr & 0x40)
00160       vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem: NVAL > NMVAL\n";
00161     if (ierr & 0x80)
00162       vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem: NVAL > MAXOP\n";
00163     if (ierr & 0x100)
00164       vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem: NVAL > MAXJ/2\n";
00165     if (ierr & 0x200)
00166       vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem: NBLOCK < 1\n";
00167   }
00168   else if (ierr < 0)
00169   {
00170     if (ierr == -1)
00171       vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem:\n"
00172                << "  poor initial vectors chosen\n";
00173     else if (ierr == -2)
00174       vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem:\n"
00175                << "  reached maximum operations " << maxop
00176                << " without finding all eigenvalues,\n"
00177                << "  found " << nperm << " eigenvalues\n";
00178     else if (ierr == -8)
00179       vcl_cerr << "Error: vnl_sparse_symmetric_eigensystem:\n"
00180                << "  disastrous loss of orthogonality - internal error\n";
00181   }
00182 
00183   // Copy the eigenvalues and vectors.
00184   nvalues = n;
00185   vectors = new vnl_vector<double>[n];
00186   values = new double[n];
00187   for (int i=0; i<n; ++i) {
00188     values[i] = temp_vals[i];
00189 #if 0
00190     vcl_cout << "value " << temp_vals[i]
00191              << " accuracy " << temp_vals[i+n*2] << vcl_endl;
00192 #endif
00193     vnl_vector<double> vec(dim,0.0);
00194     for (int j=0; j<dim; ++j)
00195       vec[j] = temp_vecs[j + dim*i];
00196     vectors[i] = vec;
00197   }
00198 
00199   // Delete temporary space.
00200   for (unsigned i=0; i<temp_store.size(); ++i)
00201     delete [] temp_store[i];
00202   temp_store.clear();
00203 
00204   return ierr;
00205 }
00206 
00207 
00208 //------------------------------------------------------------
00209 //: Here is where the fortran converted code gets called.
00210 // The sparse matrix A is assumed to be symmetric.
00211 // Find n eigenvalue/eigenvectors of the eigenproblem A * x = lambda * B * x.
00212 // !smallest and !magnitude - compute the N largest (algebraic) eigenvalues
00213 //  smallest and !magnitude - compute the N smallest (algebraic) eigenvalues
00214 // !smallest and  magnitude - compute the N largest (magnitude) eigenvalues
00215 //  smallest and  magnitude - compute the nev smallest (magnitude) eigenvalues
00216 // set sigma for shift/invert mode
00217 int vnl_sparse_symmetric_eigensystem::CalculateNPairs(
00218                       vnl_sparse_matrix<double>& A, vnl_sparse_matrix<double>& B, int nEV,
00219                       double tolerance, int numberLanczosVecs,
00220                       bool smallest, bool magnitude,
00221                       int maxIterations,
00222                       double sigma)
00223 {
00224   mat = &A;
00225   Bmat = &B;
00226 
00227   // Clear current vectors.
00228   if (vectors) {
00229     delete[] vectors; vectors = NULL;
00230     delete[] values; values = NULL;
00231   }
00232   nvalues = 0;
00233 
00234   const long whichLength = 2;
00235   char which[whichLength + 1];
00236   which[whichLength] = '\0';
00237   if (smallest)
00238     which[0] = 'S';
00239   else
00240     which[0] = 'L';
00241 
00242   if (magnitude)
00243     which[1] = 'M';
00244   else
00245     which[1] = 'A';
00246 
00247   long  matSize = mat->columns();          // Dimension of the eigenproblem.
00248   long  ido = 0;        // ido == 0 means initialization
00249 
00250   long  nconv = 0L;     // Number of "converged" Ritz values.
00251   long  numberLanczosVecsL = numberLanczosVecs;    // number of vectors to calc
00252   long  nEVL = nEV;    // long number of EVs to calc
00253 
00254   double *resid = new double[matSize];
00255   vcl_memset((void*) resid, 0, sizeof(double)*matSize);
00256 
00257   if (maxIterations <= 0)
00258     maxIterations =  nEVL * 100;
00259 
00260   if (numberLanczosVecsL <= 0)
00261     numberLanczosVecsL = 2 * nEVL + 1;
00262   numberLanczosVecsL = (numberLanczosVecsL > matSize ? matSize : numberLanczosVecsL);
00263   double *V = new double[matSize * numberLanczosVecsL + 1];
00264 
00265 #define DONE 99
00266   const int genEigProblemLength = 1;
00267   char genEigProblem = 'G';
00268   long    info = 0;   // Initialization info (INPUT) and error flag (OUTPUT)
00269 
00270 #define IPARAMSIZE 12
00271   long iParam[IPARAMSIZE];
00272   // for the sake of consistency with parameter indices in FORTRAN,
00273   // start at index 1...
00274   iParam[0] = 0;
00275   iParam[1] = 1;   //  always auto-shift
00276   iParam[2] = 0;   //  no longer referenced
00277   iParam[3] = maxIterations;
00278   iParam[4] = 1;   // NB: blocksize to be used in the recurrence.
00279                    // The code currently works only for NB = 1.
00280 
00281   iParam[5] = 0;   // output - number of converged Ritz values
00282   iParam[6] = 0;   // No longer referenced. Implicit restarting is ALWAYS used
00283 
00284   long mode;
00285 
00286   // if we have a sigma, it's mode 3, otherwise, mode 2
00287   // the mode determines the OP used in the solution
00288   // determine OP
00289   vnl_sparse_matrix<double> OP;
00290   if (sigma != 0.0)
00291   {
00292     // K*x = lambda*M*x, K symmetric, M symmetric positive semi-definite
00293     // OP = (inv[K - sigma*M])*M  and  B = M.
00294     // Shift-and-Invert mode
00295     mode = 3;
00296     // determine OP
00297 
00298     OP = B;
00299     OP *= sigma;
00300     OP = A - OP;
00301 //vsl_print_summary(std::cout, OP);
00302   }
00303   else
00304   {
00305     // A*x = lambda*M*x, A symmetric, M symmetric positive definite
00306     // OP = inv[M]*A  and  B = M.
00307     mode = 2;
00308     OP = B;
00309   }
00310   // iParam[7] is the mode of the solution
00311   iParam[7] = mode;
00312 
00313   // decompose for using in "multiplying" intermediate results
00314   vnl_sparse_lu opLU(OP);
00315 
00316 //std::cout << opLU << std::endl;
00317 
00318   iParam[8] = 0;   //  parameter for user supplied shifts - not used here
00319 
00320   // iParam 9 - 11 are output
00321   iParam[9] = 0;   // total number of OP*x operations
00322   iParam[10] = 0;   // total number of B*x operations if BMAT='G'
00323   iParam[11] = 0;   // total number of steps of re-orthogonalization
00324 
00325   // output vector filled with address information for intermediate data used
00326   // by the solver
00327   // use FORTRAN indexing again...
00328   long iPntr[IPARAMSIZE];
00329   for (int clrIx = 0; clrIx < IPARAMSIZE; clrIx++)
00330     iPntr[clrIx]= 0;
00331 
00332   // Double precision work array of length 3*N.
00333   double *workd = new double[3 * matSize + 1];
00334 
00335   // Double precision work array of length 3*N.
00336   long lworkl = numberLanczosVecsL * (numberLanczosVecsL+9);
00337 
00338   // Double precision work array of length at least NCV**2 + 8*NCV
00339   double *workl = new double[lworkl + 1];
00340 
00341   // start from scratch
00342   bool basisCalculated = false;
00343 
00344   vnl_vector<double> workVector;
00345 
00346   while (!basisCalculated)
00347   {
00348     // Calling arpack routine dsaupd.
00349     v3p_netlib_dsaupd_(
00350       &ido, &genEigProblem, &matSize, which,
00351       &nEVL, &tolerance, resid, &numberLanczosVecsL, &V[1], &matSize,
00352           &iParam[1], &iPntr[1], &workd[1], &workl[1], &lworkl, &info,
00353           genEigProblemLength, whichLength);
00354 
00355     // Checking if aupp is done
00356     if (ido==DONE)
00357     {
00358       nconv = iParam[5];
00359       basisCalculated = true;
00360       break;
00361     }
00362     else
00363     {
00364       switch (info) {
00365         case    -8:  // Could not perform LAPACK eigenvalue calculation
00366         case    -9:  // Starting vector is zero
00367         case -9999:  // Could not build an Arnoldi factorization
00368           return info;
00369           break;
00370         case     0:  // success
00371         case     1:  // hit maxIterations - should be DONE
00372         case     3:  // No shifts could be applied during a cycle of IRAM iteration
00373           break;
00374         default   :  // unknown ARPACK error
00375           return info;
00376       }
00377 
00378       // setting z pointer to ( = Bx) into workd
00379       if (ido == -1)
00380         iPntr[3] = iPntr[2] + matSize;
00381 
00382       vnl_vector_ref<double> x(matSize, &workd[iPntr[1]]);
00383       vnl_vector_ref<double> y(matSize, &workd[iPntr[2]]);
00384       vnl_vector_ref<double> z(matSize, &workd[iPntr[3]]);  // z = Bx
00385 
00386       switch (ido)
00387       {
00388         case -1:
00389             // Performing y <- OP*x for the first time when mode != 2.
00390             if (mode != 2)
00391               B.mult(x, z);
00392             // no "break;" - initialization continues below
00393         case  1:
00394             // Performing y <- OP*w.
00395             if (mode != 2)
00396               opLU.solve(z, &y);
00397             else
00398               {
00399               A.mult(x, workVector);
00400               x.update(workVector);
00401               opLU.solve(x, &y);
00402               }
00403           break;
00404         case  2:
00405             B.mult(x, y);
00406           break;
00407         default:
00408             break;
00409       }
00410     }
00411   }
00412 
00413   long rvec   = 1;  // get the values and vectors
00414 
00415   // which Ritz vctors do we want?
00416   const int howMnyLength = 1;
00417   char howMny = 'A';  // all
00418 
00419   // selection vector for which Ritz vectors to calc.
00420   // we want them all, so allocate the space (dseupd uses it)
00421   vnl_vector<long> select(numberLanczosVecsL);
00422 
00423   // allocate eVals and eVecs
00424   nvalues = nconv;
00425   values = new double[nvalues];
00426   vectors = new vnl_vector<double>[nvalues];
00427 
00428   // hold the eigenvectors
00429   double *Z = new double[nvalues * matSize];
00430 
00431   v3p_netlib_dseupd_(&rvec, &howMny, select.data_block(), values, Z, &matSize, &sigma, &genEigProblem,
00432                      &matSize, which, &nEVL, &tolerance, resid, &numberLanczosVecsL, &V[1], &matSize, &iParam[1],
00433                      &iPntr[1], &workd[1], &workl[1], &lworkl, &info,
00434                      howMnyLength, genEigProblemLength, whichLength);
00435 
00436   // Copy the eigenvectors
00437   int evIx;
00438   for (evIx = 0; evIx < nvalues; evIx++)
00439   {
00440     vnl_vector_ref<double> tempEVec(matSize, &Z[evIx * matSize]);
00441     vectors[evIx] = tempEVec;
00442   }
00443 
00444   // Delete temporary space.
00445   delete[] Z;
00446   delete[] resid;
00447   delete[] V;
00448   delete[] workd;
00449   delete[] workl;
00450 
00451   return info;
00452 }
00453 
00454 
00455 //------------------------------------------------------------
00456 //: Callback from solver to calculate the product A p.
00457 int vnl_sparse_symmetric_eigensystem::CalculateProduct(int n, int m,
00458                                                        const double* p,
00459                                                        double* q)
00460 {
00461   // Call the special multiply method on the matrix.
00462   mat->mult(n,m,p,q);
00463 
00464   return 0;
00465 }
00466 
00467 //------------------------------------------------------------
00468 //: Callback to store vectors for dnlaso.
00469 int vnl_sparse_symmetric_eigensystem::SaveVectors(int n, int m,
00470                                                   const double* q,
00471                                                   int base)
00472 {
00473   // Store the contents of q.  Basically this is a fifo.  When a write
00474   // with base=0 is called, we start another fifo.
00475   if (base == 0) {
00476     for (unsigned i=0; i<temp_store.size(); ++i)
00477       delete temp_store[i];
00478     temp_store.clear();
00479   }
00480 
00481   double* temp = new double[n*m];
00482   vcl_memcpy(temp,q,n*m*sizeof(double));
00483 #ifdef DEBUG
00484     vcl_cout << "Save vectors " << base << ' ' << temp << '\n';
00485 #endif
00486 
00487   temp_store.push_back(temp);
00488 
00489   return 0;
00490 }
00491 
00492 //------------------------------------------------------------
00493 //: Callback to restore vectors for dnlaso.
00494 int vnl_sparse_symmetric_eigensystem::RestoreVectors(int n, int m,
00495                                                      double* q,
00496                                                      int base)
00497 {
00498   // Store the contents of q.  Basically this is a fifo.  When a read
00499   // with base=0 is called, we start another fifo.
00500   static int read_idx = 0;
00501   if (base == 0)
00502     read_idx = 0;
00503 
00504   double* temp = temp_store[read_idx];
00505   vcl_memcpy(q,temp,n*m*sizeof(double));
00506 #ifdef DEBUG
00507     vcl_cout << "Restore vectors " << base << ' ' << temp << '\n';
00508 #endif
00509 
00510   read_idx++;
00511   return 0;
00512 }
00513 
00514 //------------------------------------------------------------
00515 //: Return a calculated eigenvector.
00516 vnl_vector<double> vnl_sparse_symmetric_eigensystem::get_eigenvector(int i) const
00517 {
00518   assert(i>=0 && i<nvalues);
00519   return vectors[i];
00520 }
00521 
00522 double vnl_sparse_symmetric_eigensystem::get_eigenvalue(int i) const
00523 {
00524   assert(i>=0 && i<nvalues);
00525   return values[i];
00526 }