Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
example/petra_power_method/cxx_main.cpp
Go to the documentation of this file.
1//@HEADER
2// ************************************************************************
3//
4// Epetra: Linear Algebra Services Package
5// Copyright 2011 Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
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 <cstdio>
44#include <cstdlib>
45#include <cassert>
46#include <string>
47#include <cmath>
48#include <vector>
49#include "Epetra_Map.h"
50#include "Epetra_Time.h"
51#include "Epetra_MultiVector.h"
52#include "Epetra_Util.h"
53#include "Epetra_Vector.h"
54#include "Epetra_CrsMatrix.h"
55#ifdef EPETRA_MPI
56#include "mpi.h"
57#include "Epetra_MpiComm.h"
58#endif
59#ifndef __cplusplus
60#define __cplusplus
61#endif
62#include "Epetra_Comm.h"
63#include "Epetra_SerialComm.h"
64#include "Epetra_Version.h"
65
66// prototype
67int power_method(Epetra_CrsMatrix& A, double & lambda, int niters, double tolerance,
68 bool verbose);
69
70
71int main(int argc, char *argv[])
72{
73 int ierr = 0, i;
74
75#ifdef EPETRA_MPI
76
77 // Initialize MPI
78
79 MPI_Init(&argc,&argv);
80
81 Epetra_MpiComm Comm( MPI_COMM_WORLD );
82
83#else
84
86
87#endif
88
89 int MyPID = Comm.MyPID();
90 int NumProc = Comm.NumProc();
91 bool verbose = (MyPID==0);
92
93 if (verbose)
94 std::cout << Epetra_Version() << std::endl << std::endl;
95
96 std::cout << Comm << std::endl;
97
98 // Get the number of local equations from the command line
99 if (argc!=2)
100 {
101 if (verbose)
102 std::cout << "Usage: " << argv[0] << " number_of_equations" << std::endl;
103 std::exit(1);
104 }
105 int NumGlobalElements = std::atoi(argv[1]);
106
107 if (NumGlobalElements < NumProc)
108 {
109 if (verbose)
110 std::cout << "numGlobalBlocks = " << NumGlobalElements
111 << " cannot be < number of processors = " << NumProc << std::endl;
112 std::exit(1);
113 }
114
115 // Construct a Map that puts approximately the same number of
116 // equations on each processor.
117
118 Epetra_Map Map(NumGlobalElements, 0, Comm);
119
120 // Get update list and number of local equations from newly created Map.
121
122 int NumMyElements = Map.NumMyElements();
123
124 std::vector<int> MyGlobalElements(NumMyElements);
125 Map.MyGlobalElements(&MyGlobalElements[0]);
126
127 // Create an integer vector NumNz that is used to build the Petra Matrix.
128 // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation
129 // on this processor
130
131 std::vector<int> NumNz(NumMyElements);
132
133 // We are building a tridiagonal matrix where each row has (-1 2 -1)
134 // So we need 2 off-diagonal terms (except for the first and last equation)
135
136 for (i=0; i<NumMyElements; i++)
137 if (MyGlobalElements[i]==0 || MyGlobalElements[i] == NumGlobalElements-1)
138 NumNz[i] = 2;
139 else
140 NumNz[i] = 3;
141
142 // Create a Epetra_Matrix
143
144 Epetra_CrsMatrix A(Copy, Map, &NumNz[0]);
145
146 // Add rows one-at-a-time
147 // Need some vectors to help
148 // Off diagonal Values will always be -1
149
150
151 std::vector<double> Values(2);
152 Values[0] = -1.0; Values[1] = -1.0;
153 std::vector<int> Indices(2);
154 double two = 2.0;
155 int NumEntries;
156
157 for (i=0; i<NumMyElements; i++)
158 {
159 if (MyGlobalElements[i]==0)
160 {
161 Indices[0] = 1;
162 NumEntries = 1;
163 }
164 else if (MyGlobalElements[i] == NumGlobalElements-1)
165 {
166 Indices[0] = NumGlobalElements-2;
167 NumEntries = 1;
168 }
169 else
170 {
171 Indices[0] = MyGlobalElements[i]-1;
172 Indices[1] = MyGlobalElements[i]+1;
173 NumEntries = 2;
174 }
175 ierr = A.InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
176 assert(ierr==0);
177 // Put in the diagonal entry
178 ierr = A.InsertGlobalValues(MyGlobalElements[i], 1, &two, &MyGlobalElements[i]);
179 assert(ierr==0);
180 }
181
182 // Finish up
183 ierr = A.FillComplete();
184 assert(ierr==0);
185
186 // Create vectors for Power method
187
188
189 // variable needed for iteration
190 double lambda = 0.0;
191 int niters = NumGlobalElements*10;
192 double tolerance = 1.0e-2;
193
194 // Iterate
195 Epetra_Flops counter;
196 A.SetFlopCounter(counter);
197 Epetra_Time timer(Comm);
198 ierr += power_method(A, lambda, niters, tolerance, verbose);
199 double elapsed_time = timer.ElapsedTime();
200 double total_flops =counter.Flops();
201 double MFLOPs = total_flops/elapsed_time/1000000.0;
202
203 if (verbose)
204 std::cout << "\n\nTotal MFLOPs for first solve = " << MFLOPs << std::endl<< std::endl;
205
206 // Increase diagonal dominance
207 if (verbose)
208 std::cout << "\nIncreasing magnitude of first diagonal term, solving again\n\n"
209 << std::endl;
210
211 if (A.MyGlobalRow(0)) {
212 int numvals = A.NumGlobalEntries(0);
213 std::vector<double> Rowvals(numvals);
214 std::vector<int> Rowinds(numvals);
215 A.ExtractGlobalRowCopy(0, numvals, numvals, Epetra_Util_data_ptr(Rowvals), Epetra_Util_data_ptr(Rowinds)); // Get A[0,0]
216 for (i=0; i<numvals; i++) if (Rowinds[i] == 0) Rowvals[i] *= 10.0;
217
218 A.ReplaceGlobalValues(0, numvals, Epetra_Util_data_ptr(Rowvals), Epetra_Util_data_ptr(Rowinds));
219 }
220
221 // Iterate (again)
222 lambda = 0.0;
223 timer.ResetStartTime();
224 counter.ResetFlops();
225 ierr += power_method(A, lambda, niters, tolerance, verbose);
226 elapsed_time = timer.ElapsedTime();
227 total_flops = counter.Flops();
228 MFLOPs = total_flops/elapsed_time/1000000.0;
229
230 if (verbose)
231 std::cout << "\n\nTotal MFLOPs for second solve = " << MFLOPs << std::endl<< std::endl;
232
233
234 // Release all objects
235#ifdef EPETRA_MPI
236 MPI_Finalize() ;
237#endif
238
239/* end main
240*/
241return ierr ;
242}
243
244int power_method(Epetra_CrsMatrix& A, double &lambda, int niters,
245 double tolerance, bool verbose) {
246
247 Epetra_Vector q(A.RowMap());
248 Epetra_Vector z(A.RowMap());
249 Epetra_Vector resid(A.RowMap());
250
251 Epetra_Flops * counter = A.GetFlopCounter();
252 if (counter!=0) {
253 q.SetFlopCounter(A);
254 z.SetFlopCounter(A);
255 resid.SetFlopCounter(A);
256 }
257
258 // Fill z with random Numbers
259 z.Random();
260
261 // variable needed for iteration
262 double normz, residual;
263
264 int ierr = 1;
265
266 for (int iter = 0; iter < niters; iter++)
267 {
268 z.Norm2(&normz); // Compute 2-norm of z
269 q.Scale(1.0/normz, z);
270 A.Multiply(false, q, z); // Compute z = A*q
271 q.Dot(z, &lambda); // Approximate maximum eigenvalue
272 if (iter%100==0 || iter+1==niters)
273 {
274 resid.Update(1.0, z, -lambda, q, 0.0); // Compute A*q - lambda*q
275 resid.Norm2(&residual);
276 if (verbose) std::cout << "Iter = " << iter << " Lambda = " << lambda
277 << " Residual of A*q - lambda*q = "
278 << residual << std::endl;
279 }
280 if (residual < tolerance) {
281 ierr = 0;
282 break;
283 }
284 }
285 return(ierr);
286}
T * Epetra_Util_data_ptr(std::vector< T > &vec)
Function that returns either a pointer to the first entry in the vector or, if the vector is empty,...
std::string Epetra_Version()
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
int NumMyElements() const
Number of elements on the calling processor.
Epetra_Flops * GetFlopCounter() const
Get the pointer to the Epetra_Flops() object associated with this object, returns 0 if none.
void SetFlopCounter(const Epetra_Flops &FlopCounter_in)
Set the internal Epetra_Flops() pointer.
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
int NumGlobalEntries(long long Row) const
Returns the current number of nonzero entries in specified global row on this processor.
const Epetra_Map & RowMap() const
Returns the Epetra_Map object associated with the rows of this matrix.
bool MyGlobalRow(int GID) const
Returns true of GID is owned by the calling processor, otherwise it returns false.
int FillComplete(bool OptimizeDataStorage=true)
Signal that data entry is complete. Perform transformations to local index space.
virtual int ReplaceGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Replace specified existing values with this list of entries for a given global row of the matrix.
int ExtractGlobalRowCopy(int GlobalRow, int Length, int &NumEntries, double *Values, int *Indices) const
Returns a copy of the specified global row in user-provided arrays.
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
Returns the result of a Epetra_CrsMatrix multiplied by a Epetra_Vector x in y.
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Insert a list of elements in a given global row of the matrix.
Epetra_Flops: The Epetra Floating Point Operations Class.
void ResetFlops()
Resets the number of floating point operations to zero for this multi-vector.
double Flops() const
Returns the number of floating point operations with this object and resets the count.
Epetra_Map: A class for partitioning vectors and matrices.
Definition Epetra_Map.h:119
Epetra_MpiComm: The Epetra MPI Communication Class.
int NumProc() const
Returns total number of processes.
int MyPID() const
Return my process ID.
int Scale(double ScalarValue)
Scale the current values of a multi-vector, this = ScalarValue*this.
int Dot(const Epetra_MultiVector &A, double *Result) const
Computes dot product of each corresponding pair of vectors.
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
Update multi-vector values with scaled values of A, this = ScalarThis*this + ScalarA*A.
int Random()
Set multi-vector values to random numbers.
int Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
Epetra_SerialComm: The Epetra Serial Communication Class.
Epetra_Time: The Epetra Timing Class.
Definition Epetra_Time.h:75
void ResetStartTime(void)
Epetra_Time function to reset the start time for a timer object.
double ElapsedTime(void) const
Epetra_Time elapsed time function.
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
int main(int argc, char *argv[])
int power_method(Epetra_CrsMatrix &A, double &lambda, int niters, double tolerance, bool verbose)