Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
test/SupportGraph/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"
65#include "Teuchos_RefCountPtr.hpp"
66
67
68
69// function for fancy output
70
71std::string toString(const int& x) {
72 char s[100];
73 sprintf(s, "%d", x);
74 return std::string(s);
75}
76
77std::string toString(const double& x) {
78 char s[100];
79 sprintf(s, "%g", x);
80 return std::string(s);
81}
82
83// main driver
84
85int main(int argc, char *argv[]) {
86
87#ifdef HAVE_MPI
88 MPI_Init(&argc,&argv);
89 Epetra_MpiComm Comm (MPI_COMM_WORLD);
90#else
92#endif
93
94 int MyPID = Comm.MyPID();
95 bool verbose = false;
96 if (MyPID==0) verbose = true;
97
98
99 /*int npRows = -1;
100 int npCols = -1;
101 bool useTwoD = false;
102 int randomize = 1;
103 std::string matrix = "Laplacian";
104
105 Epetra_CrsMatrix *AK = NULL;
106 std::string filename = "email.mtx";
107 read_matrixmarket_file((char*) filename.c_str(), Comm, AK,
108 useTwoD, npRows, npCols,
109 randomize, false,
110 (matrix.find("Laplacian")!=std::string::npos));
111 Teuchos::RCP<Epetra_CrsMatrix> A(AK);
112 const Epetra_Map *AMap = &(AK->DomainMap());
113 Teuchos::RCP<const Epetra_Map> Map(AMap, false);*/
114
115 int nx = 30;
116 Teuchos::ParameterList GaleriList;
117 GaleriList.set("nx", nx);
118 GaleriList.set("ny", nx * Comm.NumProc());
119 GaleriList.set("mx", 1);
120 GaleriList.set("my", Comm.NumProc());
121 Teuchos::RefCountPtr<Epetra_Map> Map = Teuchos::rcp( Galeri::CreateMap("Cartesian2D", Comm, GaleriList) );
122 Teuchos::RefCountPtr<Epetra_CrsMatrix> A = Teuchos::rcp( Galeri::CreateCrsMatrix("Laplace2D", &*Map, GaleriList) );
123
124 Teuchos::RefCountPtr<Epetra_MultiVector> LHS = Teuchos::rcp( new Epetra_MultiVector(*Map, 1) );
125 Teuchos::RefCountPtr<Epetra_MultiVector> RHS = Teuchos::rcp( new Epetra_MultiVector(*Map, 1) );
126
127
128 LHS->PutScalar(0.0); RHS->Random();
129
130 // ==================================================== //
131 // Compare support graph preconditioners to no precond. //
132 // ---------------------------------------------------- //
133
134 const double tol = 1e-5;
135 const int maxIter = 500;
136
137 // Baseline: No preconditioning
138 // Compute number of iterations, to compare to IC later.
139
140 // Here we create an AztecOO object
141 LHS->PutScalar(0.0);
142
143 AztecOO solver;
144 solver.SetUserMatrix(&*A);
145 solver.SetLHS(&*LHS);
146 solver.SetRHS(&*RHS);
147 solver.SetAztecOption(AZ_solver,AZ_cg);
148 solver.SetAztecOption(AZ_output, 16);
149 solver.Iterate(maxIter, tol);
150
151 int Iters = solver.NumIters();
152
153
154 int SupportIters;
155 Ifpack Factory;
156 Teuchos::ParameterList List;
157
158#ifdef HAVE_IFPACK_AMESOS
160 // Same test with Ifpack_SupportGraph
161 // Factored with Amesos
162
163
164 Teuchos::RefCountPtr<Ifpack_Preconditioner> PrecSupportAmesos = Teuchos::rcp( Factory.Create("MSF Amesos", &*A) );
165 List.set("amesos: solver type","Klu");
166 List.set("MST: keep diagonal", 1.0);
167 List.set("MST: randomize", 1);
168 //List.set("fact: absolute threshold", 3.0);
169
170 IFPACK_CHK_ERR(PrecSupportAmesos->SetParameters(List));
171 IFPACK_CHK_ERR(PrecSupportAmesos->Initialize());
172 IFPACK_CHK_ERR(PrecSupportAmesos->Compute());
173
174
175 // Here we create an AztecOO object
176 LHS->PutScalar(0.0);
177
178 //AztecOO solver;
179 solver.SetUserMatrix(&*A);
180 solver.SetLHS(&*LHS);
181 solver.SetRHS(&*RHS);
182 solver.SetAztecOption(AZ_solver,AZ_cg);
183 solver.SetPrecOperator(&*PrecSupportAmesos);
184 solver.SetAztecOption(AZ_output, 16);
185 solver.Iterate(maxIter, tol);
186
187 SupportIters = solver.NumIters();
188
189
190
191
192
193 // Compare to no preconditioning
194 if (SupportIters > 2*Iters)
195 IFPACK_CHK_ERR(-1);
196
197#endif
198
200 // Same test with Ifpack_SupportGraph
201 // Factored with IC
202
203
204
205 Teuchos::RefCountPtr<Ifpack_Preconditioner> PrecSupportIC = Teuchos::rcp( Factory.Create("MSF IC", &*A) );
206
207
208
209 IFPACK_CHK_ERR(PrecSupportIC->SetParameters(List));
210 IFPACK_CHK_ERR(PrecSupportIC->Compute());
211
212
213 // Here we create an AztecOO object
214 LHS->PutScalar(0.0);
215
216 //AztecOO solver;
217 solver.SetUserMatrix(&*A);
218 solver.SetLHS(&*LHS);
219 solver.SetRHS(&*RHS);
220 solver.SetAztecOption(AZ_solver,AZ_cg);
221 solver.SetPrecOperator(&*PrecSupportIC);
222 solver.SetAztecOption(AZ_output, 16);
223 solver.Iterate(maxIter, tol);
224
225 SupportIters = solver.NumIters();
226
227 // Compare to no preconditioning
228 if (SupportIters > 2*Iters)
229 IFPACK_CHK_ERR(-1);
230
231
232
233
234
235#ifdef HAVE_MPI
236 MPI_Finalize() ;
237#endif
238
239 return(EXIT_SUCCESS);
240}
#define IFPACK_CHK_ERR(ifpack_err)
#define RHS(a)
Definition: MatGenFD.c:60
int NumProc() const
int MyPID() const
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)