Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
test/CompareWithAztecOO/cxx_main.cpp
Go to the documentation of this file.
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack: Object-Oriented Algebraic Preconditioner Package
5 // Copyright (2002) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 */
42 
43 #include "Ifpack_ConfigDefs.h"
44 
45 #ifdef HAVE_MPI
46 #include "Epetra_MpiComm.h"
47 #else
48 #include "Epetra_SerialComm.h"
49 #endif
50 #include "Epetra_CrsMatrix.h"
51 #include "Epetra_Vector.h"
52 #include "Epetra_LinearProblem.h"
53 #include "Epetra_Time.h"
54 #include "Galeri_Maps.h"
55 #include "Galeri_CrsMatrices.h"
56 #include "Teuchos_ParameterList.hpp"
57 #include "Teuchos_RefCountPtr.hpp"
58 #include "Ifpack_AdditiveSchwarz.h"
59 #include "AztecOO.h"
61 #include "Ifpack_PointRelaxation.h"
62 #include "Ifpack_IC.h"
63 #include "Ifpack_ILU.h"
64 #include "Ifpack_Amesos.h"
65 
66 bool verbose = false;
67 
68 bool CompareWithAztecOO(Epetra_LinearProblem& Problem, const std::string what,
69  int Overlap, int ival)
70 {
71  using std::cout;
72  using std::endl;
73 
74  AztecOO AztecOOSolver(Problem);
75  AztecOOSolver.SetAztecOption(AZ_solver,AZ_gmres);
76  AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
77  AztecOOSolver.SetAztecOption(AZ_overlap,Overlap);
78  AztecOOSolver.SetAztecOption(AZ_graph_fill,ival);
79  AztecOOSolver.SetAztecOption(AZ_poly_ord, ival);
80  AztecOOSolver.SetAztecParam(AZ_drop, 0.0);
81  AztecOOSolver.SetAztecParam(AZ_athresh, 0.0);
82  AztecOOSolver.SetAztecParam(AZ_rthresh, 0.0);
83 
84  Epetra_MultiVector& RHS = *(Problem.GetRHS());
85  Epetra_MultiVector& LHS = *(Problem.GetLHS());
86  Teuchos::RefCountPtr<Epetra_RowMatrix> A = Teuchos::rcp(Problem.GetMatrix(), false);
87 
88  LHS.Random();
89  A->Multiply(false,LHS,RHS);
90 
92  List.set("fact: level-of-fill", ival);
93  List.set("relaxation: sweeps", ival);
94  List.set("relaxation: damping factor", 1.0);
95  List.set("relaxation: zero starting solution", true);
96 
97  //default combine mode is as for AztecOO
98  List.set("schwarz: combine mode", Zero);
99 
100  Epetra_Time Time(A->Comm());
101 
102  Teuchos::RefCountPtr<Ifpack_Preconditioner> Prec;
103 
104  if (what == "Jacobi") {
105  Prec = Teuchos::rcp( new Ifpack_PointRelaxation(&*A) );
106  List.set("relaxation: type", "Jacobi");
107  AztecOOSolver.SetAztecOption(AZ_precond,AZ_Jacobi);
108  AztecOOSolver.SetAztecOption(AZ_reorder,0);
109  }
110  else if (what == "IC no reord") {
111  Prec = Teuchos::rcp( new Ifpack_AdditiveSchwarz<Ifpack_IC>(&*A,Overlap) );
112  AztecOOSolver.SetAztecOption(AZ_precond,AZ_dom_decomp);
113  AztecOOSolver.SetAztecOption(AZ_subdomain_solve,AZ_icc);
114  AztecOOSolver.SetAztecOption(AZ_reorder,0);
115  }
116  else if (what == "IC reord") {
117  Prec = Teuchos::rcp( new Ifpack_AdditiveSchwarz<Ifpack_IC>(&*A,Overlap) );
118  List.set("schwarz: use reordering", true);
119  AztecOOSolver.SetAztecOption(AZ_precond,AZ_dom_decomp);
120  AztecOOSolver.SetAztecOption(AZ_subdomain_solve,AZ_icc);
121  AztecOOSolver.SetAztecOption(AZ_reorder,1);
122  }
123  else if (what == "ILU no reord") {
124  Prec = Teuchos::rcp( new Ifpack_AdditiveSchwarz<Ifpack_ILU>(&*A,Overlap) );
125  AztecOOSolver.SetAztecOption(AZ_precond,AZ_dom_decomp);
126  AztecOOSolver.SetAztecOption(AZ_subdomain_solve,AZ_ilu);
127  AztecOOSolver.SetAztecOption(AZ_reorder,0);
128  }
129  else if (what == "ILU reord") {
130  Prec = Teuchos::rcp( new Ifpack_AdditiveSchwarz<Ifpack_ILU>(&*A,Overlap) );
131  List.set("schwarz: use reordering", true);
132  AztecOOSolver.SetAztecOption(AZ_precond,AZ_dom_decomp);
133  AztecOOSolver.SetAztecOption(AZ_subdomain_solve,AZ_ilu);
134  AztecOOSolver.SetAztecOption(AZ_reorder,1);
135  }
136 #ifdef HAVE_IFPACK_AMESOS
137  else if (what == "LU") {
138  Prec = Teuchos::rcp( new Ifpack_AdditiveSchwarz<Ifpack_Amesos>(&*A,Overlap) );
139  List.set("amesos: solver type", "Klu");
140  AztecOOSolver.SetAztecOption(AZ_precond,AZ_dom_decomp);
141  AztecOOSolver.SetAztecOption(AZ_subdomain_solve,AZ_lu);
142  }
143 #endif
144  else {
145  cerr << "Option not recognized" << endl;
146  exit(EXIT_FAILURE);
147  }
148 
149  // ==================================== //
150  // Solve with AztecOO's preconditioners //
151  // ==================================== //
152 
153  LHS.PutScalar(0.0);
154 
155  Time.ResetStartTime();
156  AztecOOSolver.Iterate(150,1e-5);
157 
158  if (verbose) {
159  cout << endl;
160  cout << "==================================================" << endl;
161  cout << "Testing `" << what << "', Overlap = "
162  << Overlap << ", ival = " << ival << endl;
163  cout << endl;
164  cout << "[AztecOO] Total time = " << Time.ElapsedTime() << " (s)" << endl;
165  cout << "[AztecOO] Residual = " << AztecOOSolver.TrueResidual() << " (s)" << endl;
166  cout << "[AztecOO] Iterations = " << AztecOOSolver.NumIters() << endl;
167  cout << endl;
168  }
169 
170  int AztecOOPrecIters = AztecOOSolver.NumIters();
171 
172  // =========================================== //
173  // Create the IFPACK preconditioner and solver //
174  // =========================================== //
175 
176  Epetra_Time Time2(A->Comm());
177  assert(Prec != Teuchos::null);
178  IFPACK_CHK_ERR(Prec->SetParameters(List));
179 
180  Time.ResetStartTime();
181  IFPACK_CHK_ERR(Prec->Initialize());
182  if (verbose)
183  cout << "[IFPACK] Time for Initialize() = "
184  << Time.ElapsedTime() << " (s)" << endl;
185 
186  Time.ResetStartTime();
187  IFPACK_CHK_ERR(Prec->Compute());
188  if (verbose)
189  cout << "[IFPACK] Time for Compute() = "
190  << Time.ElapsedTime() << " (s)" << endl;
191 
192 
193  AztecOOSolver.SetPrecOperator(&*Prec);
194 
195  LHS.PutScalar(0.0);
196 
197  Time.ResetStartTime();
198  AztecOOSolver.Iterate(150,1e-5);
199 
200  if (verbose) {
201  cout << "[IFPACK] Total time = " << Time2.ElapsedTime() << " (s)" << endl;
202  cout << "[IFPACK] Residual = " << AztecOOSolver.TrueResidual() << " (s)" << endl;
203  cout << "[IFPACK] Iterations = " << AztecOOSolver.NumIters() << endl;
204  cout << endl;
205  }
206 
207  int IFPACKPrecIters = AztecOOSolver.NumIters();
208 
209  if (IFPACK_ABS(AztecOOPrecIters - IFPACKPrecIters) > 3) {
210  cerr << "TEST FAILED (" << AztecOOPrecIters << " != "
211  << IFPACKPrecIters << ")" << endl;
212  return(false);
213  }
214  else
215  return(true);
216 
217 }
218 
219 // ======================================================================
220 int main(int argc, char *argv[])
221 {
222 
223 #ifdef HAVE_MPI
224  MPI_Init(&argc,&argv);
225  Epetra_MpiComm Comm(MPI_COMM_WORLD);
226 #else
227  Epetra_SerialComm Comm;
228 #endif
229 
230  int nx = 30;
231  Teuchos::ParameterList GaleriList;
232  GaleriList.set("n", nx * nx);
233  GaleriList.set("nx", nx);
234  GaleriList.set("ny", nx);
235 
236  Teuchos::RefCountPtr<Epetra_Map> Map = Teuchos::rcp( Galeri::CreateMap("Linear", Comm, GaleriList) );
237  Teuchos::RefCountPtr<Epetra_RowMatrix> A = Teuchos::rcp( Galeri::CreateCrsMatrix("Laplace2D", &*Map, GaleriList) );
238  Epetra_Vector LHS(*Map);
239  Epetra_Vector RHS(*Map);
240  Epetra_LinearProblem Problem(&*A, &LHS, &RHS);
241 
242  int TestPassed = true;
243 
244  // Jacobi as in AztecOO (no overlap)
245  for (int ival = 1 ; ival < 10 ; ival += 3) {
246  TestPassed = TestPassed &&
247  CompareWithAztecOO(Problem,"Jacobi",0,ival);
248  }
249 
250 #if 0
251  // AztecOO with IC and overlap complains, also with
252  // large fill-ins (in parallel)
253  TestPassed = TestPassed &&
254  CompareWithAztecOO(Problem,"IC no reord",0,0);
255  TestPassed = TestPassed &&
256  CompareWithAztecOO(Problem,"IC reord",0,0);
257 
258  vector<std::string> Tests;
259  // now test solvers that accept overlap
260  Tests.push_back("ILU no reord");
261  Tests.push_back("ILU reord");
262  // following requires --enable-aztecoo-azlu
263 #ifdef HAVE_IFPACK_AMESOS
264  //Tests.push_back("LU");
265 #endif
266 
267  for (unsigned int i = 0 ; i < Tests.size() ; ++i) {
268  for (int overlap = 0 ; overlap < 1 ; overlap += 2) {
269  for (int ival = 0 ; ival < 10 ; ival += 4)
270  TestPassed = TestPassed &&
271  CompareWithAztecOO(Problem,Tests[i],overlap,ival);
272  }
273  }
274 #endif
275 
276  if (!TestPassed) {
277  cerr << "Test `CompareWithAztecOO.exe' FAILED!" << endl;
278  exit(EXIT_FAILURE);
279  }
280 
281 #ifdef HAVE_MPI
282  MPI_Finalize() ;
283 #endif
284  cout << "Test `CompareWithAztecOO.exe' passed!" << endl;
285 
286  exit(EXIT_SUCCESS);
287 }
Epetra_MultiVector * GetLHS() const
Epetra_MultiVector * GetRHS() const
bool CompareWithAztecOO(Epetra_LinearProblem &Problem, const std::string what, int Overlap, int ival)
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Ifpack_AdditiveSchwarz: a class to define Additive Schwarz preconditioners of Epetra_RowMatrix&#39;s.
Ifpack_PointRelaxation: a class to define point relaxation preconditioners of for Epetra_RowMatrix&#39;s...
int main(int argc, char *argv[])
Epetra_RowMatrix * GetMatrix() const
#define IFPACK_ABS(x)
#define IFPACK_CHK_ERR(ifpack_err)
#define RHS(a)
Definition: MatGenFD.c:60