Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
test/IC/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 
70 std::string toString(const int& x) {
71  char s[100];
72  sprintf(s, "%d", x);
73  return std::string(s);
74 }
75 
76 std::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 
84 int main(int argc, char *argv[]) {
85 
86 #ifdef HAVE_MPI
87  MPI_Init(&argc,&argv);
88  Epetra_MpiComm Comm (MPI_COMM_WORLD);
89 #else
90  Epetra_SerialComm Comm;
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::CreateMap("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 
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
220  // 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 }
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
std::string toString(const int &x)
Ifpack_IlukGraph: A class for constructing level filled graphs for use with ILU(k) class precondition...
int MyPID() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
const double tol
int main(int argc, char *argv[])
Ifpack_CrsRick: A class for constructing and using an incomplete Cholesky (IC) factorization of a giv...
int NumProc() const
Ifpack: a function class to define Ifpack preconditioners.
Definition: Ifpack.h:138
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
#define IFPACK_CHK_ERR(ifpack_err)
#define RHS(a)
Definition: MatGenFD.c:60