Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
test/IC_LL/cxx_main.cpp
Go to the documentation of this file.
1
2/*@HEADER
3// ***********************************************************************
4//
5// Ifpack: Object-Oriented Algebraic Preconditioner Package
6// Copyright (2002) Sandia Corporation
7//
8// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9// license for use of this work by or on behalf of the U.S. Government.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39//
40// ***********************************************************************
41//@HEADER
42*/
43
44#include "Ifpack_ConfigDefs.h"
45
46#include "Epetra_ConfigDefs.h"
47#ifdef HAVE_MPI
48#include "mpi.h"
49#include "Epetra_MpiComm.h"
50#else
51#include "Epetra_SerialComm.h"
52#endif
53#include "Epetra_Comm.h"
54#include "Epetra_Map.h"
55#include "Epetra_Time.h"
56#include "Epetra_BlockMap.h"
57#include "Epetra_MultiVector.h"
58#include "Epetra_Vector.h"
59#include "Epetra_Export.h"
60#include "AztecOO.h"
61#include "Galeri_Maps.h"
62#include "Galeri_CrsMatrices.h"
63//#include "Ifpack_CrsRick.h"
64#include "Ifpack.h"
66#include "Teuchos_RefCountPtr.hpp"
67
68// function for fancy output
69
70std::string toString(const int& x) {
71 char s[100];
72 sprintf(s, "%d", x);
73 return std::string(s);
74}
75
76std::string toString(const double& x) {
77 char s[100];
78 sprintf(s, "%g", x);
79 return std::string(s);
80}
81
82// main driver
83
84int main(int argc, char *argv[]) {
85
86#ifdef HAVE_MPI
87 MPI_Init(&argc,&argv);
88 Epetra_MpiComm Comm (MPI_COMM_WORLD);
89#else
91#endif
92
93 // The problem is defined on a 2D grid, global size is nx * nx.
94 int nx = 30;
95 Teuchos::ParameterList GaleriList;
96 GaleriList.set("nx", nx);
97 GaleriList.set("ny", nx * Comm.NumProc());
98 GaleriList.set("mx", 1);
99 GaleriList.set("my", Comm.NumProc());
100 Teuchos::RefCountPtr<Epetra_Map> Map = Teuchos::rcp( Galeri::CreateMap64("Cartesian2D", Comm, GaleriList) );
101 Teuchos::RefCountPtr<Epetra_CrsMatrix> A = Teuchos::rcp( Galeri::CreateCrsMatrix("Laplace2D", &*Map, GaleriList) );
102 Teuchos::RefCountPtr<Epetra_MultiVector> LHS = Teuchos::rcp( new Epetra_MultiVector(*Map, 1) );
103 Teuchos::RefCountPtr<Epetra_MultiVector> RHS = Teuchos::rcp( new Epetra_MultiVector(*Map, 1) );
104 LHS->PutScalar(0.0); RHS->Random();
105
106 // ========================================= //
107 // Compare IC preconditioners to no precond. //
108 // ----------------------------------------- //
109
110 const double tol = 1e-5;
111 const int maxIter = 500;
112
113 // Baseline: No preconditioning
114 // Compute number of iterations, to compare to IC later.
115
116 // Here we create an AztecOO object
117 LHS->PutScalar(0.0);
118
119 AztecOO solver;
120 solver.SetUserMatrix(&*A);
121 solver.SetLHS(&*LHS);
122 solver.SetRHS(&*RHS);
123 solver.SetAztecOption(AZ_solver,AZ_cg);
124 //solver.SetPrecOperator(&*PrecDiag);
125 solver.SetAztecOption(AZ_output, 16);
126 solver.Iterate(maxIter, tol);
127
128 int Iters = solver.NumIters();
129 //cout << "No preconditioner iterations: " << Iters << endl;
130
131#if 0
132 // Not sure how to use Ifpack_CrsRick - leave out for now.
133 //
134 // I wanna test funky values to be sure that they have the same
135 // influence on the algorithms, both old and new
136 int LevelFill = 2;
137 double DropTol = 0.3333;
138 double Condest;
139
140 Teuchos::RefCountPtr<Ifpack_CrsRick> IC;
141 Ifpack_IlukGraph mygraph (A->Graph(), 0, 0);
142 IC = Teuchos::rcp( new Ifpack_CrsRick(*A, mygraph) );
143 IC->SetAbsoluteThreshold(0.00123);
144 IC->SetRelativeThreshold(0.9876);
145 // Init values from A
146 IC->InitValues(*A);
147 // compute the factors
148 IC->Factor();
149 // and now estimate the condition number
150 IC->Condest(false,Condest);
151
152 if( Comm.MyPID() == 0 ) {
153 cout << "Condition number estimate (level-of-fill = "
154 << LevelFill << ") = " << Condest << endl;
155 }
156
157 // Define label for printing out during the solve phase
158 std::string label = "Ifpack_CrsRick Preconditioner: LevelFill = " + toString(LevelFill) +
159 " Overlap = 0";
160 IC->SetLabel(label.c_str());
161
162 // Here we create an AztecOO object
163 LHS->PutScalar(0.0);
164
165 AztecOO solver;
166 solver.SetUserMatrix(&*A);
167 solver.SetLHS(&*LHS);
168 solver.SetRHS(&*RHS);
169 solver.SetAztecOption(AZ_solver,AZ_cg);
170 solver.SetPrecOperator(&*IC);
171 solver.SetAztecOption(AZ_output, 16);
172 solver.Iterate(maxIter, tol);
173
174 int RickIters = solver.NumIters();
175 //cout << "Ifpack_Rick iterations: " << RickIters << endl;
176
177 // Compare to no preconditioning
178 if (RickIters > Iters/2)
179 IFPACK_CHK_ERR(-1);
180
181#endif
182
184 // Same test with Ifpack_IC
185 // This is Crout threshold Cholesky, so different than IC(0)
186
187 Ifpack Factory;
188 Teuchos::RefCountPtr<Ifpack_Preconditioner> PrecIC = Teuchos::rcp( Factory.Create("IC", &*A) );
189
190 Teuchos::ParameterList List;
191 //List.get("fact: ict level-of-fill", 2.);
192 //List.get("fact: drop tolerance", 0.3333);
193 //List.get("fact: absolute threshold", 0.00123);
194 //List.get("fact: relative threshold", 0.9876);
195 //List.get("fact: relaxation value", 0.0);
196
197 IFPACK_CHK_ERR(PrecIC->SetParameters(List));
198 IFPACK_CHK_ERR(PrecIC->Compute());
199
200 // Here we create an AztecOO object
201 LHS->PutScalar(0.0);
202
203 //AztecOO solver;
204 solver.SetUserMatrix(&*A);
205 solver.SetLHS(&*LHS);
206 solver.SetRHS(&*RHS);
207 solver.SetAztecOption(AZ_solver,AZ_cg);
208 solver.SetPrecOperator(&*PrecIC);
209 solver.SetAztecOption(AZ_output, 16);
210 solver.Iterate(maxIter, tol);
211
212 int ICIters = solver.NumIters();
213 //cout << "Ifpack_IC iterations: " << ICIters << endl;
214
215 // Compare to no preconditioning
216 if (ICIters > Iters/2)
217 IFPACK_CHK_ERR(-1);
218
219#if 0
221 // Same test with Ifpack_ICT
222 // This is another threshold Cholesky
223
224 Teuchos::RefCountPtr<Ifpack_Preconditioner> PrecICT = Teuchos::rcp( Factory.Create("ICT", &*A) );
225
226 //Teuchos::ParameterList List;
227 //List.get("fact: level-of-fill", 2);
228 //List.get("fact: drop tolerance", 0.3333);
229 //List.get("fact: absolute threshold", 0.00123);
230 //List.get("fact: relative threshold", 0.9876);
231 //List.get("fact: relaxation value", 0.0);
232
233 IFPACK_CHK_ERR(PrecICT->SetParameters(List));
234 IFPACK_CHK_ERR(PrecICT->Compute());
235
236 // Here we create an AztecOO object
237 LHS->PutScalar(0.0);
238
239 solver.SetUserMatrix(&*A);
240 solver.SetLHS(&*LHS);
241 solver.SetRHS(&*RHS);
242 solver.SetAztecOption(AZ_solver,AZ_cg);
243 solver.SetPrecOperator(&*PrecICT);
244 solver.SetAztecOption(AZ_output, 16);
245 solver.Iterate(maxIter, tol);
246
247 int ICTIters = solver.NumIters();
248 //cout << "Ifpack_ICT iterations: " << ICTIters << endl;
249
250 // Compare to no preconditioning
251 if (ICTIters > Iters/2)
252 IFPACK_CHK_ERR(-1);
253#endif
254
255#ifdef HAVE_MPI
256 MPI_Finalize() ;
257#endif
258
259 return(EXIT_SUCCESS);
260}
#define IFPACK_CHK_ERR(ifpack_err)
#define RHS(a)
Definition: MatGenFD.c:60
int NumProc() const
int MyPID() const
Ifpack_CrsRick: A class for constructing and using an incomplete Cholesky (IC) factorization of a giv...
Ifpack_IlukGraph: A class for constructing level filled graphs for use with ILU(k) class precondition...
static Ifpack_Preconditioner * Create(EPrecType PrecType, Epetra_RowMatrix *Matrix, const int overlap=0, bool overrideSerialDefault=false)
Creates an instance of Ifpack_Preconditioner given the enum value of the preconditioner type (can not...
Definition: Ifpack.cpp:253
const double tol
int main(int argc, char *argv[])
std::string toString(const int &x)