BlockKrylovSchur/BlockKrylovSchurEpetraExSVD.cpp

This is an example of how to use the Anasazi::BlockKrylovSchurSolMgr solver manager to compute an SVD, using Epetra data structures.

00001 //  This example shows how to use the block Krylov-Schur method to compute a few
00002 //  of the largest singular values (sigma) and corresponding right singular 
00003 //  vectors (v) for the matrix A by solving the symmetric problem:
00004 //
00005 //                             (A'*A)*v = sigma*v
00006 //
00007 //  where A is an m by n real matrix that is derived from the simplest finite
00008 //  difference discretization of the 2-dimensional kernel K(s,t)dt where
00009 //
00010 //                    K(s,t) = s(t-1)   if 0 <= s <= t <= 1
00011 //                             t(s-1)   if 0 <= t <= s <= 1
00012 //
00013 //  NOTE:  This example came from the ARPACK SVD driver dsvd.f
00014 //
00015 #include "AnasaziBlockKrylovSchurSolMgr.hpp"
00016 #include "AnasaziBasicEigenproblem.hpp"
00017 #include "AnasaziConfigDefs.hpp"
00018 #include "AnasaziEpetraAdapter.hpp"
00019 
00020 #include "Epetra_CrsMatrix.h"
00021 #include "Teuchos_LAPACK.hpp"
00022 
00023 #ifdef EPETRA_MPI
00024 #include "Epetra_MpiComm.h"
00025 #include <mpi.h>
00026 #else
00027 #include "Epetra_SerialComm.h"
00028 #endif
00029 #include "Epetra_Map.h"
00030 
00031 int main(int argc, char *argv[]) {
00032   int i, j, info;
00033   const double one = 1.0;
00034   const double zero = 0.0;
00035   Teuchos::LAPACK<int,double> lapack;
00036 
00037 #ifdef EPETRA_MPI
00038   // Initialize MPI
00039   MPI_Init(&argc,&argv);
00040   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00041 #else
00042   Epetra_SerialComm Comm;
00043 #endif
00044 
00045   int MyPID = Comm.MyPID();
00046 
00047   //  Dimension of the matrix
00048   int m = 500;
00049   int n = 100;
00050 
00051   // Construct a Map that puts approximately the same number of
00052   // equations on each processor.
00053 
00054   Epetra_Map RowMap(m, 0, Comm);
00055   Epetra_Map ColMap(n, 0, Comm);
00056 
00057   // Get update list and number of local equations from newly created Map.
00058 
00059   int NumMyRowElements = RowMap.NumMyElements();
00060   
00061   std::vector<int> MyGlobalRowElements(NumMyRowElements);
00062   RowMap.MyGlobalElements(&MyGlobalRowElements[0]);
00063 
00064   /* We are building an m by n matrix with entries
00065     
00066               A(i,j) = k*(si)*(tj - 1) if i <= j
00067                = k*(tj)*(si - 1) if i  > j
00068   
00069      where si = i/(m+1) and tj = j/(n+1) and k = 1/(n+1).
00070   */
00071 
00072   // Create an Epetra_Matrix
00073   Teuchos::RCP<Epetra_CrsMatrix> A = Teuchos::rcp( new Epetra_CrsMatrix(Copy, RowMap, n) );
00074 
00075   // Compute coefficients for discrete integral operator
00076   std::vector<double> Values(n);
00077   std::vector<int> Indices(n);
00078   double inv_mp1 = one/(m+1);
00079   double inv_np1 = one/(n+1);
00080   for (i=0; i<n; i++) { Indices[i] = i; }
00081   
00082   for (i=0; i<NumMyRowElements; i++) {
00083     //
00084     for (j=0; j<n; j++) {
00085       //
00086       if ( MyGlobalRowElements[i] <= j ) {
00087         Values[j] = inv_np1 * ( (MyGlobalRowElements[i]+one)*inv_mp1 ) * ( (j+one)*inv_np1 - one );  // k*(si)*(tj-1)
00088       }
00089       else {
00090         Values[j] = inv_np1 * ( (j+one)*inv_np1 ) * ( (MyGlobalRowElements[i]+one)*inv_mp1 - one );  // k*(tj)*(si-1)
00091       }
00092     }
00093     info = A->InsertGlobalValues(MyGlobalRowElements[i], n, &Values[0], &Indices[0]);
00094     assert( info==0 );
00095   }
00096 
00097   // Finish up
00098   info = A->FillComplete(ColMap, RowMap);
00099   assert( info==0 );
00100   info = A->OptimizeStorage();
00101   assert( info==0 );
00102   A->SetTracebackMode(1); // Shutdown Epetra Warning tracebacks
00103 
00104   //************************************
00105   // Start the block Arnoldi iteration
00106   //***********************************
00107   //
00108   //  Variables used for the Block Arnoldi Method
00109   //
00110   int nev = 4;
00111   int blockSize = 1;
00112   int numBlocks = 10;
00113   int maxRestarts = 20;
00114   int verbosity = Anasazi::Errors + Anasazi::Warnings + Anasazi::FinalSummary;
00115   double tol = lapack.LAMCH('E');
00116   std::string which = "LM";
00117   //
00118   // Create parameter list to pass into solver
00119   //
00120   Teuchos::ParameterList MyPL;
00121   MyPL.set( "Verbosity", verbosity );
00122   MyPL.set( "Which", which );
00123   MyPL.set( "Block Size", blockSize );
00124   MyPL.set( "Num Blocks", numBlocks );
00125   MyPL.set( "Maximum Restarts", maxRestarts );
00126   MyPL.set( "Convergence Tolerance", tol );
00127 
00128   typedef Anasazi::MultiVec<double> MV;
00129   typedef Anasazi::Operator<double> OP;
00130 
00131   // Create an Anasazi::EpetraMultiVec for an initial vector to start the solver. 
00132   // Note:  This needs to have the same number of columns as the blocksize.
00133   Teuchos::RCP<Anasazi::EpetraMultiVec> ivec = Teuchos::rcp( new Anasazi::EpetraMultiVec(ColMap, blockSize) );
00134   ivec->MvRandom();
00135 
00136   // Call the constructor for the (A^T*A) operator
00137   Teuchos::RCP<Anasazi::EpetraSymOp>  Amat = Teuchos::rcp( new Anasazi::EpetraSymOp(A) );  
00138   Teuchos::RCP<Anasazi::BasicEigenproblem<double, MV, OP> > MyProblem =
00139     Teuchos::rcp( new Anasazi::BasicEigenproblem<double, MV, OP>(Amat, ivec) );
00140 
00141   // Inform the eigenproblem that the matrix A is symmetric
00142   MyProblem->setHermitian(true);
00143 
00144   // Set the number of eigenvalues requested and the blocksize the solver should use
00145   MyProblem->setNEV( nev );
00146 
00147   // Inform the eigenproblem that you are finished passing it information
00148   bool boolret = MyProblem->setProblem();
00149   if (boolret != true) {
00150     if (MyPID == 0) {
00151       cout << "Anasazi::BasicEigenproblem::setProblem() returned with error." << endl;
00152     }
00153 #ifdef HAVE_MPI
00154     MPI_Finalize() ;
00155 #endif
00156     return -1;
00157   }
00158 
00159   // Initialize the Block Arnoldi solver
00160   Anasazi::BlockKrylovSchurSolMgr<double, MV, OP> MySolverMgr(MyProblem, MyPL);
00161   
00162   // Solve the problem to the specified tolerances or length
00163   Anasazi::ReturnType returnCode = MySolverMgr.solve();
00164   if (returnCode != Anasazi::Converged && MyPID==0) {
00165     cout << "Anasazi::EigensolverMgr::solve() returned unconverged." << endl;
00166   }
00167 
00168   // Get the eigenvalues and eigenvectors from the eigenproblem
00169   Anasazi::Eigensolution<double,MV> sol = MyProblem->getSolution();
00170   std::vector<Anasazi::Value<double> > evals = sol.Evals;
00171   int numev = sol.numVecs;
00172 
00173   if (numev > 0) {
00174     
00175     // Compute singular values/vectors and direct residuals.
00176     //
00177     // Compute singular values which are the square root of the eigenvalues
00178     if (MyPID==0) {
00179       cout<<"------------------------------------------------------"<<endl;
00180       cout<<"Computed Singular Values: "<<endl;
00181       cout<<"------------------------------------------------------"<<endl;
00182     }
00183     for (i=0; i<numev; i++) { evals[i].realpart = Teuchos::ScalarTraits<double>::squareroot( evals[i].realpart ); }
00184     //
00185     // Compute left singular vectors :  u = Av/sigma
00186     //
00187     std::vector<double> tempnrm(numev), directnrm(numev);
00188     std::vector<int> index(numev);
00189     for (i=0; i<numev; i++) { index[i] = i; }
00190     Anasazi::EpetraMultiVec Av(RowMap,numev), u(RowMap,numev);
00191     Anasazi::EpetraMultiVec* evecs = dynamic_cast<Anasazi::EpetraMultiVec* >(sol.Evecs->CloneView( index ));
00192     Teuchos::SerialDenseMatrix<int,double> S(numev,numev);
00193     A->Apply( *evecs, Av );
00194     Av.MvNorm( tempnrm );
00195     for (i=0; i<numev; i++) { S(i,i) = one/tempnrm[i]; };
00196     u.MvTimesMatAddMv( one, Av, S, zero );
00197     //
00198     // Compute direct residuals : || Av - sigma*u ||
00199     //
00200     for (i=0; i<numev; i++) { S(i,i) = evals[i].realpart; }
00201     Av.MvTimesMatAddMv( -one, u, S, one );
00202     Av.MvNorm( directnrm );
00203     if (MyPID==0) {
00204       cout.setf(std::ios_base::right, std::ios_base::adjustfield);
00205       cout<<std::setw(16)<<"Singular Value"
00206         <<std::setw(20)<<"Direct Residual"
00207         <<endl;
00208       cout<<"------------------------------------------------------"<<endl;
00209       for (i=0; i<numev; i++) {
00210         cout<<std::setw(16)<<evals[i].realpart
00211           <<std::setw(20)<<directnrm[i] 
00212           <<endl;
00213       }  
00214       cout<<"------------------------------------------------------"<<endl;
00215     }
00216   }
00217   
00218 #ifdef EPETRA_MPI
00219     MPI_Finalize() ;
00220 #endif
00221   
00222   return 0;
00223 }

Generated on Wed May 12 21:40:22 2010 for Anasazi by  doxygen 1.4.7