tminres
TrilinosExample/ex1.cpp
Go to the documentation of this file.
00001 // tminres is free software; you can redistribute it and/or modify it under the
00002 // terms of the GNU Lesser General Public License (as published by the Free
00003 // Software Foundation) version 2.1 dated February 1999.
00004 //
00005 // Authors:
00006 // - Umberto Villa, Emory University - uvilla@emory.edu
00007 // - Michael Saunders, Stanford University
00008 // - Santiago Akle, Stanford University
00009 
00016 #include "EpetraVectorAdapter.hpp"
00017 #include "EpetraOperatorAdapter.hpp"
00018 #include "tminres.hpp"
00019 
00020 #include <mpi.h>
00021 
00022 #include <Epetra_ConfigDefs.h>
00023 #include <Epetra_MpiComm.h>
00024 #include <Epetra_Map.h>
00025 #include <Epetra_Vector.h>
00026 #include <Epetra_CrsMatrix.h>
00027 #include <EpetraExt_CrsMatrixIn.h>
00028 #include <ml_MultiLevelPreconditioner.h>
00029 #include <Teuchos_ParameterList.hpp>
00030 
00073 int main(int argc,char * argv[])
00074 {
00075 
00076         //(1) MPI initialization and allocation of the Epetra Communicator
00077         MPI_Init(&argc, &argv);
00078         std::cout<< "MPI Initialization\n";
00079 
00080         Epetra_MpiComm * comm(new Epetra_MpiComm(MPI_COMM_WORLD));
00081 
00082         // get process information
00083         int numProc = comm->NumProc();
00084         int myPID   = comm->MyPID();
00085 
00086         bool verbose( 0 == myPID );
00087 
00088         // Open a new scope: all Epetra Object must be deleted before the Epetra_Comm is deallocated
00089         // and the MPI session is closed
00090         {
00091                    // Define the map for the parallel vectors and operators.
00092                    int dim(2312);                                       //global dimension of the problem.
00093                    Epetra_Map map(dim, 1, *comm);       //define a linear map
00094 
00095                    //(1) Read the Stokes FE Matrix and Preconditioner from file.
00096                    Epetra_CrsMatrix * A = NULL;
00097                    Epetra_CrsMatrix * P = NULL;
00098 
00099                    if( verbose )
00100                            std::cout << "Reading matrix market file (Operator)" << std::endl;
00101 
00102                    int ierr(0);
00103                    ierr = EpetraExt::MatrixMarketFileToCrsMatrix("StokesMatrix.mtx",map,map,map,A);
00104 
00105                    if(ierr==0 && verbose)
00106                            std::cout<< "Read success \n";
00107                    if(ierr != 0 && verbose)
00108                            std::cout<<" Reading error " << ierr << "\n";
00109 
00110                    if( verbose )
00111                            std::cout << "Reading matrix market file (Preconditioner)" << std::endl;
00112 
00113                    ierr = EpetraExt::MatrixMarketFileToCrsMatrix("StokesPreconditioner.mtx",map,map,map,P);
00114 
00115                    if(ierr==0 && verbose)
00116                            std::cout<< "Read success \n";
00117                    if(ierr != 0 && verbose)
00118                            std::cout<<" Reading error " << ierr << "\n";
00119 
00120                    //(2) Create the MultiLevel preconditioner for the matrix P
00121                    Teuchos::ParameterList MLList;               //List of options for ml
00122                    ML_Epetra::SetDefaults("SA",MLList); //Use the Smoothed Aggregation defaults
00123 
00124                    ML_Epetra::MultiLevelPreconditioner * MLPrec(new ML_Epetra::MultiLevelPreconditioner(*P, MLList) ); //Create the actual preconditioner
00125 
00126                    //(3) Wrap the Epetra_Operator A and MLPrec with the adapters to be used in MINRES
00127                    EpetraOperatorAdapter op(*A);
00128                    EpetraPreconditionerAdapter * prec = new EpetraPreconditionerAdapter(*MLPrec);
00129 
00130                    //(4) Define the exact solution, the rhsm and the computed solution vector (and their wrapper)
00131                    Epetra_Vector ev_sol(map), ev_rhs(map), ev_x(map);
00132                    EpetraVectorAdapter sol(ev_sol), rhs(ev_rhs), x(ev_x);
00133 
00134                    sol.Randomize( 1 );  //Pick a random exact solution ||sol||_2 = 1
00135 
00136                    op.Apply(sol, rhs);  // rhs = op * sol
00137                    double rhsNorm2(0);
00138                    rhsNorm2 = InnerProduct(rhs, rhs);
00139                    if(verbose)
00140                            std::cout<<"|| rhs ||_2 = " << sqrt(rhsNorm2) << "\n";
00141 
00142                    //(5) Set up minres parameters
00143                    x = 0.0;                             //Let the initial guess to be identically 0
00144                    double shift(0);             //No shift in minres
00145                    int max_iter(500);   //Set the maximum number of iterations
00146                    double tol(1.e-10);  //Set the tolerance of the stopping criterion
00147 
00148                    //(6) Solve the problem with minres
00149                    MINRES(op, x, rhs, prec, shift, max_iter, tol, verbose);
00150 
00151                    //(7) Compare against the exact solution
00152                    subtract(x, sol, x);
00153                    double err2 = InnerProduct(x,x);
00154 
00155                    if(verbose)
00156                            std::cout<< "|| x_ex - x_n || = " << sqrt(err2) << "\n";
00157 
00158                    delete prec;
00159                    delete MLPrec;
00160                    delete P;
00161                    delete A;
00162 
00163 
00164         }
00165 
00166         //(8) Dellacate the comm and finalize mpi
00167         delete comm;
00168         MPI_Finalize();
00169 
00170         return 0;
00171 }
 All Classes Files Functions Friends