Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
test/CompareWithAztecOO_LL/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 
72  AztecOO AztecOOSolver(Problem);
73  AztecOOSolver.SetAztecOption(AZ_solver,AZ_gmres);
74  AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
75  AztecOOSolver.SetAztecOption(AZ_overlap,Overlap);
76  AztecOOSolver.SetAztecOption(AZ_graph_fill,ival);
77  AztecOOSolver.SetAztecOption(AZ_poly_ord, ival);
78  AztecOOSolver.SetAztecParam(AZ_drop, 0.0);
79  AztecOOSolver.SetAztecParam(AZ_athresh, 0.0);
80  AztecOOSolver.SetAztecParam(AZ_rthresh, 0.0);
81 
82  Epetra_MultiVector& RHS = *(Problem.GetRHS());
83  Epetra_MultiVector& LHS = *(Problem.GetLHS());
84  Teuchos::RefCountPtr<Epetra_RowMatrix> A = Teuchos::rcp(Problem.GetMatrix(), false);
85 
86  LHS.Random();
87  A->Multiply(false,LHS,RHS);
88 
90  List.set("fact: level-of-fill", ival);
91  List.set("relaxation: sweeps", ival);
92  List.set("relaxation: damping factor", 1.0);
93  List.set("relaxation: zero starting solution", true);
94 
95  //default combine mode is as for AztecOO
96  List.set("schwarz: combine mode", Zero);
97 
98  Epetra_Time Time(A->Comm());
99 
100  Teuchos::RefCountPtr<Ifpack_Preconditioner> Prec;
101 
102  if (what == "Jacobi") {
103  Prec = Teuchos::rcp( new Ifpack_PointRelaxation(&*A) );
104  List.set("relaxation: type", "Jacobi");
105  AztecOOSolver.SetAztecOption(AZ_precond,AZ_Jacobi);
106  AztecOOSolver.SetAztecOption(AZ_reorder,0);
107  }
108  else if (what == "IC no reord") {
109  Prec = Teuchos::rcp( new Ifpack_AdditiveSchwarz<Ifpack_IC>(&*A,Overlap) );
110  AztecOOSolver.SetAztecOption(AZ_precond,AZ_dom_decomp);
111  AztecOOSolver.SetAztecOption(AZ_subdomain_solve,AZ_icc);
112  AztecOOSolver.SetAztecOption(AZ_reorder,0);
113  }
114  else if (what == "IC reord") {
115  Prec = Teuchos::rcp( new Ifpack_AdditiveSchwarz<Ifpack_IC>(&*A,Overlap) );
116  List.set("schwarz: use reordering", true);
117  AztecOOSolver.SetAztecOption(AZ_precond,AZ_dom_decomp);
118  AztecOOSolver.SetAztecOption(AZ_subdomain_solve,AZ_icc);
119  AztecOOSolver.SetAztecOption(AZ_reorder,1);
120  }
121  else if (what == "ILU no reord") {
122  Prec = Teuchos::rcp( new Ifpack_AdditiveSchwarz<Ifpack_ILU>(&*A,Overlap) );
123  AztecOOSolver.SetAztecOption(AZ_precond,AZ_dom_decomp);
124  AztecOOSolver.SetAztecOption(AZ_subdomain_solve,AZ_ilu);
125  AztecOOSolver.SetAztecOption(AZ_reorder,0);
126  }
127  else if (what == "ILU reord") {
128  Prec = Teuchos::rcp( new Ifpack_AdditiveSchwarz<Ifpack_ILU>(&*A,Overlap) );
129  List.set("schwarz: use reordering", true);
130  AztecOOSolver.SetAztecOption(AZ_precond,AZ_dom_decomp);
131  AztecOOSolver.SetAztecOption(AZ_subdomain_solve,AZ_ilu);
132  AztecOOSolver.SetAztecOption(AZ_reorder,1);
133  }
134 #ifdef HAVE_IFPACK_AMESOS
135  else if (what == "LU") {
136  Prec = Teuchos::rcp( new Ifpack_AdditiveSchwarz<Ifpack_Amesos>(&*A,Overlap) );
137  List.set("amesos: solver type", "Klu");
138  AztecOOSolver.SetAztecOption(AZ_precond,AZ_dom_decomp);
139  AztecOOSolver.SetAztecOption(AZ_subdomain_solve,AZ_lu);
140  }
141 #endif
142  else {
143  cerr << "Option not recognized" << endl;
144  exit(EXIT_FAILURE);
145  }
146 
147  // ==================================== //
148  // Solve with AztecOO's preconditioners //
149  // ==================================== //
150 
151  LHS.PutScalar(0.0);
152 
153  Time.ResetStartTime();
154  AztecOOSolver.Iterate(150,1e-5);
155 
156  if (verbose) {
157  cout << endl;
158  cout << "==================================================" << endl;
159  cout << "Testing `" << what << "', Overlap = "
160  << Overlap << ", ival = " << ival << endl;
161  cout << endl;
162  cout << "[AztecOO] Total time = " << Time.ElapsedTime() << " (s)" << endl;
163  cout << "[AztecOO] Residual = " << AztecOOSolver.TrueResidual() << " (s)" << endl;
164  cout << "[AztecOO] Iterations = " << AztecOOSolver.NumIters() << endl;
165  cout << endl;
166  }
167 
168  int AztecOOPrecIters = AztecOOSolver.NumIters();
169 
170  // =========================================== //
171  // Create the IFPACK preconditioner and solver //
172  // =========================================== //
173 
174  Epetra_Time Time2(A->Comm());
175  assert(Prec != Teuchos::null);
176  IFPACK_CHK_ERR(Prec->SetParameters(List));
177 
178  Time.ResetStartTime();
179  IFPACK_CHK_ERR(Prec->Initialize());
180  if (verbose)
181  cout << "[IFPACK] Time for Initialize() = "
182  << Time.ElapsedTime() << " (s)" << endl;
183 
184  Time.ResetStartTime();
185  IFPACK_CHK_ERR(Prec->Compute());
186  if (verbose)
187  cout << "[IFPACK] Time for Compute() = "
188  << Time.ElapsedTime() << " (s)" << endl;
189 
190 
191  AztecOOSolver.SetPrecOperator(&*Prec);
192 
193  LHS.PutScalar(0.0);
194 
195  Time.ResetStartTime();
196  AztecOOSolver.Iterate(150,1e-5);
197 
198  if (verbose) {
199  cout << "[IFPACK] Total time = " << Time2.ElapsedTime() << " (s)" << endl;
200  cout << "[IFPACK] Residual = " << AztecOOSolver.TrueResidual() << " (s)" << endl;
201  cout << "[IFPACK] Iterations = " << AztecOOSolver.NumIters() << endl;
202  cout << endl;
203  }
204 
205  int IFPACKPrecIters = AztecOOSolver.NumIters();
206 
207  if (IFPACK_ABS(AztecOOPrecIters - IFPACKPrecIters) > 3) {
208  cerr << "TEST FAILED (" << AztecOOPrecIters << " != "
209  << IFPACKPrecIters << ")" << endl;
210  return(false);
211  }
212  else
213  return(true);
214 
215 }
216 
217 // ======================================================================
218 int main(int argc, char *argv[])
219 {
220 
221 #ifdef HAVE_MPI
222  MPI_Init(&argc,&argv);
223  Epetra_MpiComm Comm(MPI_COMM_WORLD);
224 #else
225  Epetra_SerialComm Comm;
226 #endif
227 
228  int nx = 30;
229  Teuchos::ParameterList GaleriList;
230  GaleriList.set("n", nx * nx);
231  GaleriList.set("nx", nx);
232  GaleriList.set("ny", nx);
233 
234  Teuchos::RefCountPtr<Epetra_Map> Map = Teuchos::rcp( Galeri::CreateMap64("Linear", Comm, GaleriList) );
235  Teuchos::RefCountPtr<Epetra_RowMatrix> A = Teuchos::rcp( Galeri::CreateCrsMatrix("Laplace2D", &*Map, GaleriList) );
236  Epetra_Vector LHS(*Map);
237  Epetra_Vector RHS(*Map);
238  Epetra_LinearProblem Problem(&*A, &LHS, &RHS);
239 
240  int TestPassed = true;
241 
242  // Jacobi as in AztecOO (no overlap)
243  for (int ival = 1 ; ival < 10 ; ival += 3) {
244  TestPassed = TestPassed &&
245  CompareWithAztecOO(Problem,"Jacobi",0,ival);
246  }
247 
248 #if 0
249  // AztecOO with IC and overlap complains, also with
250  // large fill-ins (in parallel)
251  TestPassed = TestPassed &&
252  CompareWithAztecOO(Problem,"IC no reord",0,0);
253  TestPassed = TestPassed &&
254  CompareWithAztecOO(Problem,"IC reord",0,0);
255 
256  vector<std::string> Tests;
257  // now test solvers that accept overlap
258  Tests.push_back("ILU no reord");
259  Tests.push_back("ILU reord");
260  // following requires --enable-aztecoo-azlu
261 #ifdef HAVE_IFPACK_AMESOS
262  //Tests.push_back("LU");
263 #endif
264 
265  for (unsigned int i = 0 ; i < Tests.size() ; ++i) {
266  for (int overlap = 0 ; overlap < 1 ; overlap += 2) {
267  for (int ival = 0 ; ival < 10 ; ival += 4)
268  TestPassed = TestPassed &&
269  CompareWithAztecOO(Problem,Tests[i],overlap,ival);
270  }
271  }
272 #endif
273 
274  if (!TestPassed) {
275  cerr << "Test `CompareWithAztecOO.exe' FAILED!" << endl;
276  exit(EXIT_FAILURE);
277  }
278 
279 #ifdef HAVE_MPI
280  MPI_Finalize() ;
281 #endif
282  cout << "Test `CompareWithAztecOO.exe' passed!" << endl;
283 
284  exit(EXIT_SUCCESS);
285 }
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