Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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"
64
65#include "Ifpack_Amesos.h"
66#include "AztecOO.h"
67
68static bool verbose = false;
69static bool SymmetricGallery = false;
70static bool Solver = AZ_gmres;
71const int NumVectors = 3;
72
73// ======================================================================
74bool 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
83 Teuchos::ParameterList 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// ======================================================================
167int main(int argc, char *argv[])
168{
169#ifdef HAVE_MPI
170 MPI_Init(&argc,&argv);
171 Epetra_MpiComm Comm( MPI_COMM_WORLD );
172#else
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;
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}
#define RHS(a)
Definition: MatGenFD.c:60
virtual int MyPID() const=0
Epetra_RowMatrix * GetMatrix() const
int NumProc() const
int MyPID() const
int PutScalar(double ScalarConstant)
virtual const Epetra_Comm & Comm() const=0
Ifpack_BlockRelaxation: a class to define block relaxation preconditioners of Epetra_RowMatrix's.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
virtual int Compute()
Computes the preconditioner.
Ifpack_PointRelaxation: a class to define point relaxation preconditioners of for Epetra_RowMatrix's.
virtual int Compute()
Computes the preconditioners.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
int main(int argc, char *argv[])
static bool Solver
Definition: performance.cpp:70
const int NumVectors
Definition: performance.cpp:71
static bool verbose
Definition: performance.cpp:68
bool ComparePointAndBlock(std::string PrecType, const Teuchos::RefCountPtr< Epetra_RowMatrix > &A, int sweeps)
Definition: performance.cpp:74
static bool SymmetricGallery
Definition: performance.cpp:69