Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
performance.cpp
Go to the documentation of this file.
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack: Object-Oriented Algebraic Preconditioner Package
5 // Copyright (2002) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 */
42 
43 #include <sys/resource.h>
44 #include "Ifpack_ConfigDefs.h"
45 
46 #ifdef HAVE_MPI
47 #include "Epetra_MpiComm.h"
48 #else
49 #include "Epetra_SerialComm.h"
50 #endif
51 #include "Epetra_CrsMatrix.h"
52 #include "Epetra_MultiVector.h"
53 #include "Epetra_Vector.h"
54 #include "Epetra_LinearProblem.h"
55 #include "Epetra_Map.h"
56 #include "Galeri_Maps.h"
57 #include "Galeri_CrsMatrices.h"
58 #include "Galeri_Utils.h"
59 #include "Teuchos_ParameterList.hpp"
60 #include "Teuchos_RefCountPtr.hpp"
61 #include "Ifpack_PointRelaxation.h"
62 #include "Ifpack_BlockRelaxation.h"
63 #include "Ifpack_SparseContainer.h"
64 
65 #include "Ifpack_Amesos.h"
66 #include "AztecOO.h"
67 
68 static bool verbose = false;
69 static bool SymmetricGallery = false;
70 static bool Solver = AZ_gmres;
71 const int NumVectors = 3;
72 
73 // ======================================================================
74 bool ComparePointAndBlock(std::string PrecType, const Teuchos::RefCountPtr<Epetra_RowMatrix>& A, int sweeps)
75 {
76  Epetra_MultiVector RHS(A->RowMatrixRowMap(), NumVectors);
77  Epetra_MultiVector LHS(A->RowMatrixRowMap(), NumVectors);
78  LHS.PutScalar(0.0); RHS.Random();
79 
80  Epetra_LinearProblem Problem(&*A, &LHS, &RHS);
81 
82  // Set up the list
84  List.set("relaxation: damping factor", 1.0);
85  List.set("relaxation: type", PrecType);
86  List.set("relaxation: sweeps",sweeps);
87  List.set("partitioner: type", "linear");
88  List.set("partitioner: local parts", A->NumMyRows());
89 
90  int ItersPoint, ItersBlock;
91 
92  // ================================================== //
93  // get the number of iterations with point relaxation //
94  // ================================================== //
95  {
96  RHS.PutScalar(1.0);
97  LHS.PutScalar(0.0);
98 
99  Ifpack_PointRelaxation Point(&*A);
100  Point.SetParameters(List);
101  Point.Compute();
102 
103  // set AztecOO solver object
104  AztecOO AztecOOSolver(Problem);
105  AztecOOSolver.SetAztecOption(AZ_solver,Solver);
106  AztecOOSolver.SetAztecOption(AZ_output,32);
107  AztecOOSolver.SetPrecOperator(&Point);
108 
109  AztecOOSolver.Iterate(1550,1e-2);
110 
111  double TrueResidual = AztecOOSolver.TrueResidual();
112  ItersPoint = AztecOOSolver.NumIters();
113  // some output
114  if (verbose && Problem.GetMatrix()->Comm().MyPID() == 0) {
115  cout << "Iterations = " << ItersPoint << endl;
116  cout << "Norm of the true residual = " << TrueResidual << endl;
117  }
118  }
119  if(verbose) printf(" point finished \n");
120  // ================================================== //
121  // get the number of iterations with block relaxation //
122  // ================================================== //
123  {
124 
125  RHS.PutScalar(1.0);
126  LHS.PutScalar(0.0);
127 
129  Block.SetParameters(List);
130  Block.Compute();
131 
132  // set AztecOO solver object
133  AztecOO AztecOOSolver(Problem);
134  AztecOOSolver.SetAztecOption(AZ_solver,Solver);
135  AztecOOSolver.SetAztecOption(AZ_output,32);
136  AztecOOSolver.SetPrecOperator(&Block);
137 
138  AztecOOSolver.Iterate(1550,1e-2);
139 
140  double TrueResidual = AztecOOSolver.TrueResidual();
141  ItersBlock = AztecOOSolver.NumIters();
142  // some output
143  if ( Problem.GetMatrix()->Comm().MyPID() == 0) {
144  cout << "Iterations " << ItersBlock << endl;
145  cout << "Norm of the true residual = " << TrueResidual << endl;
146  }
147  }
148 
149  if(verbose) printf(" point finished \n");
150 
151  int diff = ItersPoint - ItersBlock;
152  if (diff < 0) diff = -diff;
153 
154  if (diff > 10)
155  {
156  if(verbose) cout << "ComparePointandBlock TEST FAILED!" << endl;
157  return(false);
158  }
159  else {
160  if(verbose) cout << "ComparePointandBlock TEST PASSED" << endl;
161  return(true);
162  }
163 }
164 
165 
166 // ======================================================================
167 int main(int argc, char *argv[])
168 {
169 #ifdef HAVE_MPI
170  MPI_Init(&argc,&argv);
171  Epetra_MpiComm Comm( MPI_COMM_WORLD );
172 #else
173  Epetra_SerialComm Comm;
174 #endif
175 
176  verbose = (Comm.MyPID() == 0);
177 
178  int nx = 60;
179 
180  for (int i = 1 ; i < argc ; ++i) {
181  if (strcmp(argv[i],"-s") == 0) {
182  SymmetricGallery = true;
183  Solver = AZ_cg;
184  }
185  if(strcmp(argv[i],"-n") == 0 && i+1 < argc) {
186  i++;
187  nx = atoi(argv[i]);
188  }
189 
190  }
191 
192  // size of the global matrix.
193  Teuchos::ParameterList GaleriList;
194 
195  GaleriList.set("nx", nx);
196  GaleriList.set("ny", nx * Comm.NumProc());
197  GaleriList.set("mx", 1);
198  GaleriList.set("my", Comm.NumProc());
199  Teuchos::RefCountPtr<Epetra_Map> Map = Teuchos::rcp( Galeri::CreateMap("Cartesian2D", Comm, GaleriList) );
200  Teuchos::RefCountPtr<Epetra_CrsMatrix> A;
201  if (SymmetricGallery)
202  A = Teuchos::rcp( Galeri::CreateCrsMatrix("Laplace2D", &*Map, GaleriList) );
203  else
204  A = Teuchos::rcp( Galeri::CreateCrsMatrix("Recirc2D", &*Map, GaleriList) );
205 
206  // coordinates
207  Teuchos::RCP<Epetra_MultiVector> coord = Teuchos::rcp( Galeri::CreateCartesianCoordinates("2D",&*Map,GaleriList));
208 
209  // test the preconditioner
210  int TestPassed = true;
211  struct rusage usage;
212  //int ret;
213  //ret = getrusage(who, &usage);
214  struct timeval ru_utime;
215  // struct timeval ru_stime;
216  ru_utime = usage.ru_utime;
217 
218 
219  // ================================== //
220  // compare point and block relaxation //
221  // ================================== //
222 
223  TestPassed = TestPassed &&
224  ComparePointAndBlock("Jacobi",A,10);
225  if(verbose) printf(" Jacobi Finished \n");
226 
227  //ret = getrusage(who, &usage);
228  int sec = usage.ru_utime.tv_sec -ru_utime.tv_sec;
229  int usec = usage.ru_utime.tv_usec -ru_utime.tv_usec;
230  double tt = (double)sec + 1e-6*(double)usec;
231  ru_utime = usage.ru_utime;
232  if(verbose) printf(" Jacobi time %f \n",tt);
233 
234  TestPassed = TestPassed &&
235  ComparePointAndBlock("symmetric Gauss-Seidel",A,10);
236  if(verbose) printf(" sGS finished \n");
237 
238  //ret = getrusage(who, &usage);
239  sec = usage.ru_utime.tv_sec -ru_utime.tv_sec;
240  usec = usage.ru_utime.tv_usec -ru_utime.tv_usec;
241  tt = (double)sec + 1e-6*(double)usec;
242  ru_utime = usage.ru_utime;
243  if(verbose) printf(" sGS time %f \n",tt);
244 
245  if (!SymmetricGallery) {
246  TestPassed = TestPassed &&
247  ComparePointAndBlock("Gauss-Seidel",A,10);
248  //ret = getrusage(who, &usage);
249  sec = usage.ru_utime.tv_sec -ru_utime.tv_sec;
250  usec = usage.ru_utime.tv_usec -ru_utime.tv_usec;
251  tt = (double)sec + 1e-6*(double)usec;
252  ru_utime = usage.ru_utime;
253  if(verbose) printf(" GS time %f \n",tt);
254  if(verbose) printf(" GS Finished \n");
255  }
256 
257  if (!TestPassed) {
258  cout << "Test `Performance.exe' failed!" << endl;
259  exit(EXIT_FAILURE);
260  }
261 
262 #ifdef HAVE_MPI
263  MPI_Finalize();
264 #endif
265 
266  cout << endl;
267  cout << "Test `Performance.exe' passed!" << endl;
268  cout << endl;
269  return(EXIT_SUCCESS);
270 }
Ifpack_BlockRelaxation: a class to define block relaxation preconditioners of Epetra_RowMatrix&#39;s.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
const int NumVectors
Definition: performance.cpp:71
int MyPID() const
virtual int MyPID() const =0
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
static bool Solver
Definition: performance.cpp:70
virtual const Epetra_Comm & Comm() const =0
virtual int Compute()
Computes the preconditioner.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Ifpack_PointRelaxation: a class to define point relaxation preconditioners of for Epetra_RowMatrix&#39;s...
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
int main(int argc, char *argv[])
virtual int Compute()
Computes the preconditioners.
int NumProc() const
bool ComparePointAndBlock(std::string PrecType, const Teuchos::RefCountPtr< Epetra_RowMatrix > &A, int sweeps)
Definition: performance.cpp:74
Epetra_RowMatrix * GetMatrix() const
static bool SymmetricGallery
Definition: performance.cpp:69
#define RHS(a)
Definition: MatGenFD.c:60