Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
71 std::string toString(const int& x) {
72  char s[100];
73  sprintf(s, "%d", x);
74  return std::string(s);
75 }
76 
77 std::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 
85 int main(int argc, char *argv[]) {
86 
87 #ifdef HAVE_MPI
88  MPI_Init(&argc,&argv);
89  Epetra_MpiComm Comm (MPI_COMM_WORLD);
90 #else
91  Epetra_SerialComm Comm;
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;
157 
158 #ifdef HAVE_IFPACK_AMESOS
159  // 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 }
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)
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[])
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