tminres
|
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 }