Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
test/ILU/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_CrsRiluk.h"
64#include "Ifpack.h"
65#include "Teuchos_RefCountPtr.hpp"
66
67// function for fancy output
68
69std::string toString(const int& x) {
70 char s[100];
71 sprintf(s, "%d", x);
72 return std::string(s);
73}
74
75std::string toString(const double& x) {
76 char s[100];
77 sprintf(s, "%g", x);
78 return std::string(s);
79}
80
81// main driver
82
83int main(int argc, char *argv[]) {
84
85#ifdef HAVE_MPI
86 MPI_Init(&argc,&argv);
87 Epetra_MpiComm Comm (MPI_COMM_WORLD);
88#else
90#endif
91
92 Teuchos::ParameterList GaleriList;
93 int nx = 30;
94
95 GaleriList.set("nx", nx);
96 GaleriList.set("ny", nx * Comm.NumProc());
97 GaleriList.set("mx", 1);
98 GaleriList.set("my", Comm.NumProc());
99 Teuchos::RefCountPtr<Epetra_Map> Map = Teuchos::rcp( Galeri::CreateMap("Cartesian2D", Comm, GaleriList) );
100 Teuchos::RefCountPtr<Epetra_CrsMatrix> A = Teuchos::rcp( Galeri::CreateCrsMatrix("Laplace2D", &*Map, GaleriList) );
101 Teuchos::RefCountPtr<Epetra_MultiVector> LHS = Teuchos::rcp( new Epetra_MultiVector(*Map, 1) );
102 Teuchos::RefCountPtr<Epetra_MultiVector> RHS = Teuchos::rcp( new Epetra_MultiVector(*Map, 1) );
103 LHS->PutScalar(0.0); RHS->Random();
104
105 // ============================ //
106 // Construct ILU preconditioner //
107 // ---------------------------- //
108
109 // I wanna test funky values to be sure that they have the same
110 // influence on the algorithms, both old and new
111 int LevelFill = 2;
112 double DropTol = 0.3333;
113 double Athresh = 0.0123;
114 double Rthresh = 0.9876;
115 double Relax = 0.1;
116 int Overlap = 2;
117
118 Teuchos::RefCountPtr<Ifpack_IlukGraph> Graph;
119 Teuchos::RefCountPtr<Ifpack_CrsRiluk> RILU;
120
121 Graph = Teuchos::rcp( new Ifpack_IlukGraph(A->Graph(), LevelFill, Overlap) );
122 int ierr;
123 ierr = Graph->ConstructFilledGraph();
124 IFPACK_CHK_ERR(ierr);
125
126 RILU = Teuchos::rcp( new Ifpack_CrsRiluk(*Graph) );
127 RILU->SetAbsoluteThreshold(Athresh);
128 RILU->SetRelativeThreshold(Rthresh);
129 RILU->SetRelaxValue(Relax);
130 int initerr = RILU->InitValues(*A);
131 if (initerr!=0) cout << Comm << "*ERR* InitValues = " << initerr;
132
133 RILU->Factor();
134
135 // Define label for printing out during the solve phase
136 std::string label = "Ifpack_CrsRiluk Preconditioner: LevelFill = " + toString(LevelFill) +
137 " Overlap = " + toString(Overlap) +
138 " Athresh = " + toString(Athresh) +
139 " Rthresh = " + toString(Rthresh);
140 // Here we create an AztecOO object
141 LHS->PutScalar(0.0);
142
143 int Niters = 1200;
144
145 AztecOO solver;
146 solver.SetUserMatrix(&*A);
147 solver.SetLHS(&*LHS);
148 solver.SetRHS(&*RHS);
149 solver.SetAztecOption(AZ_solver,AZ_gmres);
150 solver.SetPrecOperator(&*RILU);
151 solver.SetAztecOption(AZ_output, 16);
152 solver.Iterate(Niters, 5.0e-5);
153
154 int OldIters = solver.NumIters();
155
156 // now rebuild the same preconditioner using RILU, we expect the same
157 // number of iterations
158
159 Ifpack Factory;
160 Teuchos::RefCountPtr<Ifpack_Preconditioner> Prec = Teuchos::rcp( Factory.Create("ILU", &*A, Overlap) );
161
162 Teuchos::ParameterList List;
163 List.get("fact: level-of-fill", LevelFill);
164 List.get("fact: drop tolerance", DropTol);
165 List.get("fact: absolute threshold", Athresh);
166 List.get("fact: relative threshold", Rthresh);
167 List.get("fact: relax value", Relax);
168
169 IFPACK_CHK_ERR(Prec->SetParameters(List));
170 IFPACK_CHK_ERR(Prec->Compute());
171
172 // Here we create an AztecOO object
173 LHS->PutScalar(0.0);
174
175 solver.SetUserMatrix(&*A);
176 solver.SetLHS(&*LHS);
177 solver.SetRHS(&*RHS);
178 solver.SetAztecOption(AZ_solver,AZ_gmres);
179 solver.SetPrecOperator(&*Prec);
180 solver.SetAztecOption(AZ_output, 16);
181 solver.Iterate(Niters, 5.0e-5);
182
183 int NewIters = solver.NumIters();
184
185 if (OldIters != NewIters)
186 IFPACK_CHK_ERR(-1);
187
188
189#ifdef HAVE_IFPACK_SUPERLU
190 // Now test w/ SuperLU's ILU, if we've got it
191 Teuchos::RefCountPtr<Ifpack_Preconditioner> Prec2 = Teuchos::rcp( Factory.Create("SILU", &*A,0) );
192 Teuchos::ParameterList SList;
193 SList.set("fact: drop tolerance",1e-4);
194 SList.set("fact: zero pivot threshold",.1);
195 SList.set("fact: maximum fill factor",10.0);
196 // NOTE: There is a bug in SuperLU 4.0 which will crash the code if the maximum fill factor is set too low.
197 // This bug was reported to Sherry Li on 4/8/10.
198 SList.set("fact: silu drop rule",9);
199
200 IFPACK_CHK_ERR(Prec2->SetParameters(SList));
201 IFPACK_CHK_ERR(Prec2->Compute());
202
203 LHS->PutScalar(0.0);
204 solver.SetPrecOperator(&*Prec2);
205 solver.Iterate(Niters, 5.0e-5);
206 Prec2->Print(cout);
207
208#endif
209
210
211#ifdef HAVE_MPI
212 MPI_Finalize() ;
213#endif
214
215 return(EXIT_SUCCESS);
216}
#define IFPACK_CHK_ERR(ifpack_err)
adjacency_list< vecS, vecS, undirectedS, no_property, property< edge_weight_t, double > > Graph
#define RHS(a)
Definition: MatGenFD.c:60
int NumProc() const
Ifpack_CrsRiluk: A class for constructing and using an incomplete lower/upper (ILU) factorization of ...
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
int main(int argc, char *argv[])
std::string toString(const int &x)