BlockKrylovSchur/BlockKrylovSchurEpetraEx.cpp

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

00001 //  This example computes the specified eigenvalues of the discretized 2D Convection-Diffusion
00002 //  equation using the block Krylov-Schur method.  This discretized operator is constructed as an
00003 //  Epetra matrix, then passed into the Anasazi::EpetraOp to be used in the construction of the
00004 //  Krylov decomposition.  The specifics of the block Krylov-Schur method can be set by the user.
00005 #include "AnasaziBlockKrylovSchurSolMgr.hpp"
00006 #include "AnasaziBasicEigenproblem.hpp"
00007 #include "AnasaziConfigDefs.hpp"
00008 #include "AnasaziEpetraAdapter.hpp"
00009 #include "AnasaziBasicSort.hpp"
00010 #include "Epetra_CrsMatrix.h"
00011 
00012 #include "Teuchos_LAPACK.hpp"
00013 #include "Teuchos_CommandLineProcessor.hpp"
00014 
00015 #ifdef EPETRA_MPI
00016 #include "Epetra_MpiComm.h"
00017 #else
00018 #include "Epetra_SerialComm.h"
00019 #endif
00020 #include "Epetra_Map.h"
00021 
00022 int main(int argc, char *argv[]) {
00023 
00024   using std::cout;
00025 
00026 #ifdef EPETRA_MPI
00027   // Initialize MPI
00028   MPI_Init(&argc,&argv);
00029   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00030 #else
00031   Epetra_SerialComm Comm;
00032 #endif
00033 
00034   bool boolret;
00035   int MyPID = Comm.MyPID();
00036 
00037   bool verbose = true;
00038   bool debug = false;
00039   std::string which("SM");
00040 
00041   Teuchos::CommandLineProcessor cmdp(false,true);
00042   cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
00043   cmdp.setOption("debug","nodebug",&debug,"Print debugging information.");
00044   cmdp.setOption("sort",&which,"Targetted eigenvalues (SM,LM,SR,LR,SI,or LI).");
00045   if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
00046 #ifdef HAVE_MPI
00047     MPI_Finalize();
00048 #endif
00049     return -1;
00050   }
00051 
00052   typedef double ScalarType;
00053   typedef Teuchos::ScalarTraits<ScalarType>          SCT;
00054   typedef SCT::magnitudeType               MagnitudeType;
00055   typedef Epetra_MultiVector                          MV;
00056   typedef Epetra_Operator                             OP;
00057   typedef Anasazi::MultiVecTraits<ScalarType,MV>     MVT;
00058   typedef Anasazi::OperatorTraits<ScalarType,MV,OP>  OPT;
00059 
00060   //  Dimension of the matrix
00061   int nx = 10;        // Discretization points in any one direction.
00062   int NumGlobalElements = nx*nx;  // Size of matrix nx*nx
00063 
00064   // Construct a Map that puts approximately the same number of
00065   // equations on each processor.
00066 
00067   Epetra_Map Map(NumGlobalElements, 0, Comm);
00068 
00069   // Get update list and number of local equations from newly created Map.
00070 
00071   int NumMyElements = Map.NumMyElements();
00072 
00073   std::vector<int> MyGlobalElements(NumMyElements);
00074   Map.MyGlobalElements(&MyGlobalElements[0]);
00075 
00076   // Create an integer vector NumNz that is used to build the Petra Matrix.
00077   // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation
00078   // on this processor
00079   std::vector<int> NumNz(NumMyElements);
00080 
00081   /* We are building a matrix of block structure:
00082   
00083       | T -I          |
00084       |-I  T -I       |
00085       |   -I  T       |
00086       |        ...  -I|
00087       |           -I T|
00088 
00089    where each block is dimension nx by nx and the matrix is on the order of
00090    nx*nx.  The block T is a tridiagonal matrix. 
00091   */
00092 
00093   for (int i=0; i<NumMyElements; i++) {
00094     if (MyGlobalElements[i] == 0 || MyGlobalElements[i] == NumGlobalElements-1 || 
00095         MyGlobalElements[i] == nx-1 || MyGlobalElements[i] == nx*(nx-1) ) {
00096       NumNz[i] = 3;
00097     }
00098     else if (MyGlobalElements[i] < nx || MyGlobalElements[i] > nx*(nx-1) || 
00099              MyGlobalElements[i]%nx == 0 || (MyGlobalElements[i]+1)%nx == 0) {
00100       NumNz[i] = 4;
00101     }
00102     else {
00103       NumNz[i] = 5;
00104     }
00105   }
00106 
00107   // Create an Epetra_Matrix
00108 
00109   Teuchos::RCP<Epetra_CrsMatrix> A = Teuchos::rcp( new Epetra_CrsMatrix(Copy, Map, &NumNz[0]) );
00110 
00111   // Diffusion coefficient, can be set by user.
00112   // When rho*h/2 <= 1, the discrete convection-diffusion operator has real eigenvalues.
00113   // When rho*h/2 > 1, the operator has complex eigenvalues.
00114   double rho = 2*nx+1;
00115   
00116   // Compute coefficients for discrete convection-diffution operator
00117   const double one = 1.0;
00118   std::vector<double> Values(4);
00119   std::vector<int> Indices(4);
00120   double h = one /(nx+1);
00121   double h2 = h*h;
00122   double c = 5.0e-01*rho/ h;
00123   Values[0] = -one/h2 - c; Values[1] = -one/h2 + c; Values[2] = -one/h2; Values[3]= -one/h2;
00124   double diag = 4.0 / h2;
00125   int NumEntries, info;
00126 
00127   for (int i=0; i<NumMyElements; i++)
00128   {
00129     if (MyGlobalElements[i]==0)
00130     {
00131       Indices[0] = 1;
00132       Indices[1] = nx;
00133       NumEntries = 2;
00134       info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[1], &Indices[0]);
00135       assert( info==0 );
00136     }
00137     else if (MyGlobalElements[i] == nx*(nx-1))
00138     {
00139       Indices[0] = nx*(nx-1)+1;
00140       Indices[1] = nx*(nx-2);
00141       NumEntries = 2;
00142       info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[1], &Indices[0]);
00143       assert( info==0 );
00144     }
00145     else if (MyGlobalElements[i] == nx-1)
00146     {
00147       Indices[0] = nx-2;
00148       NumEntries = 1;
00149       info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
00150       assert( info==0 );
00151       Indices[0] = 2*nx-1;
00152       info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[2], &Indices[0]);
00153       assert( info==0 );
00154     }
00155     else if (MyGlobalElements[i] == NumGlobalElements-1)
00156     {
00157       Indices[0] = NumGlobalElements-2;
00158       NumEntries = 1;
00159       info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
00160       assert( info==0 );
00161       Indices[0] = nx*(nx-1)-1;
00162       info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[2], &Indices[0]);
00163       assert( info==0 );
00164     }
00165     else if (MyGlobalElements[i] < nx)
00166     {
00167       Indices[0] = MyGlobalElements[i]-1;
00168       Indices[1] = MyGlobalElements[i]+1;
00169       Indices[2] = MyGlobalElements[i]+nx;
00170       NumEntries = 3;
00171       info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
00172       assert( info==0 );
00173     }
00174     else if (MyGlobalElements[i] > nx*(nx-1))
00175     {
00176       Indices[0] = MyGlobalElements[i]-1;
00177       Indices[1] = MyGlobalElements[i]+1;
00178       Indices[2] = MyGlobalElements[i]-nx;
00179       NumEntries = 3;
00180       info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
00181       assert( info==0 );
00182     }
00183     else if (MyGlobalElements[i]%nx == 0)
00184     {
00185       Indices[0] = MyGlobalElements[i]+1;
00186       Indices[1] = MyGlobalElements[i]-nx;
00187       Indices[2] = MyGlobalElements[i]+nx;
00188       NumEntries = 3;
00189       info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[1], &Indices[0]);
00190       assert( info==0 );
00191     }
00192     else if ((MyGlobalElements[i]+1)%nx == 0)
00193     {
00194       Indices[0] = MyGlobalElements[i]-nx;
00195       Indices[1] = MyGlobalElements[i]+nx;
00196       NumEntries = 2;
00197       info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[2], &Indices[0]);
00198       assert( info==0 );
00199       Indices[0] = MyGlobalElements[i]-1;
00200       NumEntries = 1;
00201       info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
00202       assert( info==0 );
00203     }
00204     else
00205     {
00206       Indices[0] = MyGlobalElements[i]-1;
00207       Indices[1] = MyGlobalElements[i]+1;
00208       Indices[2] = MyGlobalElements[i]-nx;
00209       Indices[3] = MyGlobalElements[i]+nx;
00210       NumEntries = 4;
00211       info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
00212       assert( info==0 );
00213     }
00214     // Put in the diagonal entry
00215     info = A->InsertGlobalValues(MyGlobalElements[i], 1, &diag, &MyGlobalElements[i]);
00216     assert( info==0 );
00217   }
00218 
00219   // Finish up
00220   info = A->FillComplete();
00221   assert( info==0 );
00222   A->SetTracebackMode(1); // Shutdown Epetra Warning tracebacks
00223 
00224   //************************************
00225   // Start the block Arnoldi iteration
00226   //***********************************
00227   //
00228   //  Variables used for the Block Krylov Schur Method
00229   //    
00230   int nev = 4;
00231   int blockSize = 1;
00232   int numBlocks = 20;
00233   int maxRestarts = 500;
00234   //int stepSize = 5;
00235   double tol = 1e-8;
00236 
00237   // Create a sort manager to pass into the block Krylov-Schur solver manager
00238   // -->  Make sure the reference-counted pointer is of type Anasazi::SortManager<>
00239   // -->  The block Krylov-Schur solver manager uses Anasazi::BasicSort<> by default,
00240   //      so you can also pass in the parameter "Which", instead of a sort manager.
00241   Teuchos::RCP<Anasazi::SortManager<MagnitudeType> > MySort =     
00242     Teuchos::rcp( new Anasazi::BasicSort<MagnitudeType>( which ) );
00243 
00244   // Set verbosity level
00245   int verbosity = Anasazi::Errors + Anasazi::Warnings;
00246   if (verbose) {
00247     verbosity += Anasazi::FinalSummary + Anasazi::TimingDetails;
00248   }
00249   if (debug) {
00250     verbosity += Anasazi::Debug;
00251   }
00252   //
00253   // Create parameter list to pass into solver manager
00254   //
00255   Teuchos::ParameterList MyPL;
00256   MyPL.set( "Verbosity", verbosity );
00257   MyPL.set( "Sort Manager", MySort );
00258   //MyPL.set( "Which", which );  
00259   MyPL.set( "Block Size", blockSize );
00260   MyPL.set( "Num Blocks", numBlocks );
00261   MyPL.set( "Maximum Restarts", maxRestarts );
00262   //MyPL.set( "Step Size", stepSize );
00263   MyPL.set( "Convergence Tolerance", tol );
00264 
00265   // Create an Epetra_MultiVector for an initial vector to start the solver.
00266   // Note:  This needs to have the same number of columns as the blocksize.
00267   Teuchos::RCP<Epetra_MultiVector> ivec = Teuchos::rcp( new Epetra_MultiVector(Map, blockSize) );
00268   ivec->Random();
00269 
00270   // Create the eigenproblem.
00271   Teuchos::RCP<Anasazi::BasicEigenproblem<double, MV, OP> > MyProblem =
00272     Teuchos::rcp( new Anasazi::BasicEigenproblem<double, MV, OP>(A, ivec) );
00273   
00274   // Inform the eigenproblem that the operator A is symmetric
00275   MyProblem->setHermitian(rho==0.0); 
00276   
00277   // Set the number of eigenvalues requested
00278   MyProblem->setNEV( nev );
00279   
00280   // Inform the eigenproblem that you are finishing passing it information
00281   boolret = MyProblem->setProblem();
00282   if (boolret != true) {
00283     if (verbose && MyPID == 0) {
00284       cout << "Anasazi::BasicEigenproblem::setProblem() returned with error." << endl;
00285     }
00286 #ifdef HAVE_MPI
00287     MPI_Finalize() ;
00288 #endif
00289     return -1;
00290   }
00291   
00292   // Initialize the Block Arnoldi solver
00293   Anasazi::BlockKrylovSchurSolMgr<double, MV, OP> MySolverMgr(MyProblem, MyPL);
00294   
00295   // Solve the problem to the specified tolerances or length
00296   Anasazi::ReturnType returnCode = MySolverMgr.solve();
00297   if (returnCode != Anasazi::Converged && MyPID==0 && verbose) {
00298     cout << "Anasazi::EigensolverMgr::solve() returned unconverged." << endl;
00299   }
00300 
00301   // Get the Ritz values from the eigensolver
00302   std::vector<Anasazi::Value<double> > ritzValues = MySolverMgr.getRitzValues();
00303   
00304   // Output computed eigenvalues and their direct residuals
00305   if (verbose && MyPID==0) {
00306     int numritz = (int)ritzValues.size();
00307     cout.setf(std::ios_base::right, std::ios_base::adjustfield);
00308     cout<<endl<< "Computed Ritz Values"<< endl;
00309     if (MyProblem->isHermitian()) {
00310       cout<< std::setw(16) << "Real Part"
00311         << endl;
00312       cout<<"-----------------------------------------------------------"<<endl;
00313       for (int i=0; i<numritz; i++) {
00314         cout<< std::setw(16) << ritzValues[i].realpart 
00315           << endl;
00316       }  
00317       cout<<"-----------------------------------------------------------"<<endl;
00318     } 
00319     else {
00320       cout<< std::setw(16) << "Real Part"
00321         << std::setw(16) << "Imag Part"
00322         << endl;
00323       cout<<"-----------------------------------------------------------"<<endl;
00324       for (int i=0; i<numritz; i++) {
00325         cout<< std::setw(16) << ritzValues[i].realpart 
00326           << std::setw(16) << ritzValues[i].imagpart 
00327           << endl;
00328       }  
00329       cout<<"-----------------------------------------------------------"<<endl;
00330       }  
00331     }
00332 
00333   // Get the eigenvalues and eigenvectors from the eigenproblem
00334   Anasazi::Eigensolution<ScalarType,MV> sol = MyProblem->getSolution();
00335   std::vector<Anasazi::Value<ScalarType> > evals = sol.Evals;
00336   Teuchos::RCP<MV> evecs = sol.Evecs;
00337   std::vector<int> index = sol.index;
00338   int numev = sol.numVecs;
00339   
00340   if (numev > 0) {
00341     // Compute residuals.
00342     Teuchos::LAPACK<int,double> lapack;
00343     std::vector<double> normA(numev);
00344     
00345     if (MyProblem->isHermitian()) {
00346       // Get storage
00347       Epetra_MultiVector Aevecs(Map,numev);
00348       Teuchos::SerialDenseMatrix<int,double> B(numev,numev);
00349       B.putScalar(0.0); 
00350       for (int i=0; i<numev; i++) {B(i,i) = evals[i].realpart;}
00351       
00352       // Compute A*evecs
00353       OPT::Apply( *A, *evecs, Aevecs );
00354       
00355       // Compute A*evecs - lambda*evecs and its norm
00356       MVT::MvTimesMatAddMv( -1.0, *evecs, B, 1.0, Aevecs );
00357       MVT::MvNorm( Aevecs, normA );
00358       
00359       // Scale the norms by the eigenvalue
00360       for (int i=0; i<numev; i++) {
00361         normA[i] /= Teuchos::ScalarTraits<double>::magnitude( evals[i].realpart );
00362       }
00363     } else {
00364       // The problem is non-Hermitian.
00365       int i=0;
00366       std::vector<int> curind(1);
00367       std::vector<double> resnorm(1), tempnrm(1);
00368       Teuchos::RCP<MV> evecr, eveci, tempAevec;
00369       Epetra_MultiVector Aevec(Map,numev);
00370       
00371       // Compute A*evecs
00372       OPT::Apply( *A, *evecs, Aevec );
00373       
00374       Teuchos::SerialDenseMatrix<int,double> Breal(1,1), Bimag(1,1);
00375       while (i<numev) {
00376         if (index[i]==0) {
00377           // Get a view of the current eigenvector (evecr)
00378           curind[0] = i;
00379           evecr = MVT::CloneView( *evecs, curind );
00380 
00381           // Get a copy of A*evecr
00382           tempAevec = MVT::CloneCopy( Aevec, curind );
00383 
00384           // Compute A*evecr - lambda*evecr
00385           Breal(0,0) = evals[i].realpart;
00386           MVT::MvTimesMatAddMv( -1.0, *evecr, Breal, 1.0, *tempAevec );
00387 
00388           // Compute the norm of the residual and increment counter
00389           MVT::MvNorm( *tempAevec, resnorm );
00390           normA[i] = resnorm[0]/Teuchos::ScalarTraits<MagnitudeType>::magnitude( evals[i].realpart );
00391           i++;
00392         } else {
00393           // Get a view of the real part of the eigenvector (evecr)
00394           curind[0] = i;
00395           evecr = MVT::CloneView( *evecs, curind );
00396 
00397           // Get a copy of A*evecr
00398           tempAevec = MVT::CloneCopy( Aevec, curind );
00399 
00400           // Get a view of the imaginary part of the eigenvector (eveci)
00401           curind[0] = i+1;
00402           eveci = MVT::CloneView( *evecs, curind );
00403 
00404           // Set the eigenvalue into Breal and Bimag
00405           Breal(0,0) = evals[i].realpart;
00406           Bimag(0,0) = evals[i].imagpart;
00407 
00408           // Compute A*evecr - evecr*lambdar + eveci*lambdai
00409           MVT::MvTimesMatAddMv( -1.0, *evecr, Breal, 1.0, *tempAevec );
00410           MVT::MvTimesMatAddMv( 1.0, *eveci, Bimag, 1.0, *tempAevec );
00411           MVT::MvNorm( *tempAevec, tempnrm );
00412 
00413           // Get a copy of A*eveci
00414           tempAevec = MVT::CloneCopy( Aevec, curind );
00415 
00416           // Compute A*eveci - eveci*lambdar - evecr*lambdai
00417           MVT::MvTimesMatAddMv( -1.0, *evecr, Bimag, 1.0, *tempAevec );
00418           MVT::MvTimesMatAddMv( -1.0, *eveci, Breal, 1.0, *tempAevec );
00419           MVT::MvNorm( *tempAevec, resnorm );
00420 
00421           // Compute the norms and scale by magnitude of eigenvalue
00422           normA[i] = lapack.LAPY2( tempnrm[i], resnorm[i] ) /
00423             lapack.LAPY2( evals[i].realpart, evals[i].imagpart );
00424           normA[i+1] = normA[i];
00425 
00426           i=i+2;
00427         }
00428       }
00429     }
00430 
00431     // Output computed eigenvalues and their direct residuals
00432     if (verbose && MyPID==0) {
00433       cout.setf(std::ios_base::right, std::ios_base::adjustfield);
00434       cout<<endl<< "Actual Residuals"<<endl;
00435       if (MyProblem->isHermitian()) {
00436         cout<< std::setw(16) << "Real Part"
00437           << std::setw(20) << "Direct Residual"<< endl;
00438         cout<<"-----------------------------------------------------------"<<endl;
00439         for (int i=0; i<numev; i++) {
00440           cout<< std::setw(16) << evals[i].realpart 
00441             << std::setw(20) << normA[i] << endl;
00442         }  
00443         cout<<"-----------------------------------------------------------"<<endl;
00444       } 
00445       else {
00446         cout<< std::setw(16) << "Real Part"
00447           << std::setw(16) << "Imag Part"
00448           << std::setw(20) << "Direct Residual"<< endl;
00449         cout<<"-----------------------------------------------------------"<<endl;
00450         for (int i=0; i<numev; i++) {
00451           cout<< std::setw(16) << evals[i].realpart 
00452             << std::setw(16) << evals[i].imagpart 
00453             << std::setw(20) << normA[i] << endl;
00454         }  
00455         cout<<"-----------------------------------------------------------"<<endl;
00456       }  
00457     }
00458   }
00459 
00460 #ifdef EPETRA_MPI
00461     MPI_Finalize();
00462 #endif
00463 
00464   return 0;
00465 }

Generated on Tue Jul 13 09:22:46 2010 for Anasazi by  doxygen 1.4.7