Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
test/CrsRiluk/cxx_main.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 "Ifpack_ConfigDefs.h"
44#include <stdio.h>
45#include <stdlib.h>
46#include <ctype.h>
47#include <assert.h>
48#include <string.h>
49#include <math.h>
50#include "Epetra_Comm.h"
51#include "Epetra_Map.h"
52#include "Epetra_Time.h"
53#include "Epetra_CrsMatrix.h"
54#include "Epetra_Vector.h"
55#include "Epetra_Object.h"
56#include "Ifpack_Version.h"
57#ifdef EPETRA_MPI
58#include "mpi.h"
59#include "Epetra_MpiComm.h"
60#else
61#include "Epetra_SerialComm.h"
62#endif
63
64
65// prototype
69 Epetra_Vector& resid,
70 double * lambda, int niters, double tolerance,
71 bool verbose);
72
73
74int main(int argc, char *argv[])
75{
76 using std::cout;
77 using std::endl;
78
79 int ierr = 0, i;
80
81#ifdef EPETRA_MPI
82
83 // Initialize MPI
84
85 MPI_Init(&argc,&argv);
86 int rank; // My process ID
87
88 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
89 Epetra_MpiComm Comm( MPI_COMM_WORLD );
90
91#else
92
93 int rank = 0;
95
96#endif
97
98 bool verbose = false;
99
100 // Check if we should print results to standard out
101 if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
102
103 // char tmp;
104 // if (rank==0) cout << "Press any key to continue..."<< endl;
105 // if (rank==0) cin >> tmp;
106 // Comm.Barrier();
107
108 int MyPID = Comm.MyPID();
109 int NumProc = Comm.NumProc();
110
111 if (verbose && MyPID==0)
112 cout << Ifpack_Version() << endl << endl;
113
114 if (verbose) cout << "Processor "<<MyPID<<" of "<< NumProc
115 << " is alive."<<endl;
116
117 bool verbose1 = verbose;
118
119 // Redefine verbose to only print on PE 0
120 if (verbose && rank!=0) verbose = false;
121
122 int NumMyEquations = 10000;
123 int NumGlobalEquations = NumMyEquations*NumProc+EPETRA_MIN(NumProc,3);
124 if (MyPID < 3) NumMyEquations++;
125
126 // Construct a Map that puts approximately the same Number of equations on each processor
127
128 Epetra_Map Map(NumGlobalEquations, NumMyEquations, 0, Comm);
129
130 // Get update list and number of local equations from newly created Map
131 int * MyGlobalElements = new int[Map.NumMyElements()];
132 Map.MyGlobalElements(MyGlobalElements);
133
134 // Create an integer vector NumNz that is used to build the Petra Matrix.
135 // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation on this processor
136
137 int * NumNz = new int[NumMyEquations];
138
139 // We are building a tridiagonal matrix where each row has (-1 2 -1)
140 // So we need 2 off-diagonal terms (except for the first and last equation)
141
142 for (i=0; i<NumMyEquations; i++)
143 if (MyGlobalElements[i]==0 || MyGlobalElements[i] == NumGlobalEquations-1)
144 NumNz[i] = 1;
145 else
146 NumNz[i] = 2;
147
148 // Create a Epetra_Matrix
149
150 Epetra_CrsMatrix A(Copy, Map, NumNz);
151
152 // Add rows one-at-a-time
153 // Need some vectors to help
154 // Off diagonal Values will always be -1
155
156
157 double *Values = new double[2];
158 Values[0] = -1.0; Values[1] = -1.0;
159 int *Indices = new int[2];
160 double two = 2.0;
161 int NumEntries;
162
163 for (i=0; i<NumMyEquations; i++)
164 {
165 if (MyGlobalElements[i]==0)
166 {
167 Indices[0] = 1;
168 NumEntries = 1;
169 }
170 else if (MyGlobalElements[i] == NumGlobalEquations-1)
171 {
172 Indices[0] = NumGlobalEquations-2;
173 NumEntries = 1;
174 }
175 else
176 {
177 Indices[0] = MyGlobalElements[i]-1;
178 Indices[1] = MyGlobalElements[i]+1;
179 NumEntries = 2;
180 }
181 int ierr2;
182 ierr2 = A.InsertGlobalValues(MyGlobalElements[i], NumEntries, Values, Indices);
183 IFPACK_CHK_ERR(ierr2);
184 ierr2 = A.InsertGlobalValues(MyGlobalElements[i], 1, &two, MyGlobalElements+i); // Put in the diagonal entry
185 IFPACK_CHK_ERR(ierr2);
186 }
187
188 // Finish up
189 A.FillComplete();
190
191 // Create vectors for Power method
192
193 Epetra_Vector q(Map);
194 Epetra_Vector z(Map);
195 Epetra_Vector resid(Map);
196
197 // variable needed for iteration
198 double lambda = 0.0;
199 int niters = 10000;
200 double tolerance = 1.0e-3;
201
202 // Iterate
203 Epetra_Time timer(Comm);
204 ierr += power_method(A, q, z, resid, &lambda, niters, tolerance, verbose);
205 double elapsed_time = timer.ElapsedTime();
206 double total_flops = A.Flops() + q.Flops() + z.Flops() + resid.Flops();
207 double MFLOPs = total_flops/elapsed_time/1000000.0;
208
209 if (verbose) cout << "\n\nTotal MFLOPs for first solve = " << MFLOPs << endl<< endl;
210
211 // Increase diagonal dominance
212 if (verbose) cout << "\n\nIncreasing the magnitude of first diagonal term and solving again\n\n"
213 << endl;
214
215
216 if (A.MyGlobalRow(0)) {
217 int numvals = A.NumGlobalEntries(0);
218 double * Rowvals = new double [numvals];
219 int * Rowinds = new int [numvals];
220 A.ExtractGlobalRowCopy(0, numvals, numvals, Rowvals, Rowinds); // Get A[0,0]
221
222 for (i=0; i<numvals; i++) if (Rowinds[i] == 0) Rowvals[i] *= 10.0;
223
224 A.ReplaceGlobalValues(0, numvals, Rowvals, Rowinds);
225 }
226 // Iterate (again)
227 lambda = 0.0;
228 timer.ResetStartTime();
229 A.ResetFlops(); q.ResetFlops(); z.ResetFlops(); resid.ResetFlops();
230 ierr += power_method(A, q, z, resid, &lambda, niters, tolerance, verbose);
231 elapsed_time = timer.ElapsedTime();
232 total_flops = A.Flops() + q.Flops() + z.Flops() + resid.Flops();
233 MFLOPs = total_flops/elapsed_time/1000000.0;
234
235 if (verbose) cout << "\n\nTotal MFLOPs for second solve = " << MFLOPs << endl<< endl;
236
237
238 // Release all objects
239 delete [] NumNz;
240 delete [] Values;
241 delete [] Indices;
242 delete [] MyGlobalElements;
243
244
245
246 if (verbose1) {
247 // Test ostream << operator (if verbose1)
248 // Construct a Map that puts 2 equations on each PE
249
250 int NumMyElements1 = 2;
251 int NumMyEquations1 = NumMyElements1;
252 int NumGlobalEquations1 = NumMyEquations1*NumProc;
253
254 Epetra_Map Map1(-1, NumMyElements1, 0, Comm);
255
256 // Get update list and number of local equations from newly created Map
257 int * MyGlobalElements1 = new int[Map1.NumMyElements()];
258 Map1.MyGlobalElements(MyGlobalElements1);
259
260 // Create an integer vector NumNz that is used to build the Petra Matrix.
261 // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation on this processor
262
263 int * NumNz1 = new int[NumMyEquations1];
264
265 // We are building a tridiagonal matrix where each row has (-1 2 -1)
266 // So we need 2 off-diagonal terms (except for the first and last equation)
267
268 for (i=0; i<NumMyEquations1; i++)
269 if (MyGlobalElements1[i]==0 || MyGlobalElements1[i] == NumGlobalEquations1-1)
270 NumNz1[i] = 1;
271 else
272 NumNz1[i] = 2;
273
274 // Create a Epetra_Matrix
275
276 Epetra_CrsMatrix A1(Copy, Map1, NumNz1);
277
278 // Add rows one-at-a-time
279 // Need some vectors to help
280 // Off diagonal Values will always be -1
281
282
283 double *Values1 = new double[2];
284 Values1[0] = -1.0; Values1[1] = -1.0;
285 int *Indices1 = new int[2];
286 double two1 = 2.0;
287 int NumEntries1;
288
289 for (i=0; i<NumMyEquations1; i++)
290 {
291 if (MyGlobalElements1[i]==0)
292 {
293 Indices1[0] = 1;
294 NumEntries1 = 1;
295 }
296 else if (MyGlobalElements1[i] == NumGlobalEquations1-1)
297 {
298 Indices1[0] = NumGlobalEquations1-2;
299 NumEntries1 = 1;
300 }
301 else
302 {
303 Indices1[0] = MyGlobalElements1[i]-1;
304 Indices1[1] = MyGlobalElements1[i]+1;
305 NumEntries1 = 2;
306 }
307 int ierr2;
308 ierr2 = A1.InsertGlobalValues(MyGlobalElements1[i], NumEntries1, Values1, Indices1);
309 IFPACK_CHK_ERR(ierr2);
310 ierr2 = A1.InsertGlobalValues(MyGlobalElements1[i], 1, &two1, MyGlobalElements1+i); // Put in the diagonal entry
311 IFPACK_CHK_ERR(ierr2);
312 }
313
314 // Finish up
315 A1.FillComplete();
316
317 if (verbose) cout << "\n\nPrint out tridiagonal matrix, each part on each processor.\n\n" << endl;
318 cout << A1 << endl;
319
320 // Release all objects
321 delete [] NumNz1;
322 delete [] Values1;
323 delete [] Indices1;
324 delete [] MyGlobalElements1;
325
326 }
327
328#ifdef EPETRA_MPI
329 MPI_Finalize() ;
330#endif
331
332/* end main
333*/
334return ierr ;
335}
336
338 Epetra_Vector& q,
339 Epetra_Vector& z,
340 Epetra_Vector& resid,
341 double * lambda, int niters, double tolerance,
342 bool verbose) {
343 using std::cout;
344 using std::endl;
345
346 // Fill z with random Numbers
347 z.Random();
348
349 // variable needed for iteration
350 double normz, residual;
351
352 int ierr = 1;
353
354 for (int iter = 0; iter < niters; iter++)
355 {
356 z.Norm2(&normz); // Compute 2-norm of z
357 q.Scale(1.0/normz, z);
358 A.Multiply(false, q, z); // Compute z = A*q
359 q.Dot(z, lambda); // Approximate maximum eigenvaluE
360 if (iter%100==0 || iter+1==niters)
361 {
362 resid.Update(1.0, z, -(*lambda), q, 0.0); // Compute A*q - lambda*q
363 resid.Norm2(&residual);
364 if (verbose) cout << "Iter = " << iter << " Lambda = " << *lambda
365 << " Residual of A*q - lambda*q = " << residual << endl;
366 }
367 if (residual < tolerance) {
368 ierr = 0;
369 break;
370 }
371 }
372 return(ierr);
373}
#define EPETRA_MIN(x, y)
Copy
#define IFPACK_CHK_ERR(ifpack_err)
std::string Ifpack_Version()
int MyGlobalElements(int *MyGlobalElementList) const
int NumMyElements() const
void ResetFlops() const
double Flops() const
int NumGlobalEntries(long long Row) const
bool MyGlobalRow(int GID) const
int FillComplete(bool OptimizeDataStorage=true)
virtual int ReplaceGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
int ExtractGlobalRowCopy(int_type Row, int Length, int &NumEntries, double *values, int_type *Indices) const
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
int NumProc() const
int MyPID() const
int Scale(double ScalarValue)
int Dot(const Epetra_MultiVector &A, double *Result) const
void ResetStartTime(void)
double ElapsedTime(void) const
int main(int argc, char *argv[])
int power_method(Epetra_CrsMatrix &A, Epetra_Vector &q, Epetra_Vector &z, Epetra_Vector &resid, double *lambda, int niters, double tolerance, bool verbose)