Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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"
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"
59#include "AztecOO.h"
62#include "Ifpack_IC.h"
63#include "Ifpack_ILU.h"
64#include "Ifpack_Amesos.h"
65
66bool verbose = false;
67
68bool 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
91 Teuchos::ParameterList List;
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// ======================================================================
220int 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
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}
Zero
#define IFPACK_CHK_ERR(ifpack_err)
#define IFPACK_ABS(x)
#define RHS(a)
Definition MatGenFD.c:60
int PutScalar(double ScalarConstant)
Ifpack_PointRelaxation: a class to define point relaxation preconditioners of for Epetra_RowMatrix's.
int main(int argc, char *argv[])
bool CompareWithAztecOO(Epetra_LinearProblem &Problem, const std::string what, int Overlap, int ival)