Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
example/petra_power_method_LL/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_Vector.h"
53#include "Epetra_CrsMatrix.h"
54#ifdef EPETRA_MPI
55#include "mpi.h"
56#include "Epetra_MpiComm.h"
57#endif
58#ifndef __cplusplus
59#define __cplusplus
60#endif
61#include "Epetra_Comm.h"
62#include "Epetra_SerialComm.h"
63#include "Epetra_Version.h"
64
65// prototype
66int power_method(Epetra_CrsMatrix& A, double & lambda, int niters, double tolerance,
67 bool verbose);
68
69
70int main(int argc, char *argv[])
71{
72 int ierr = 0, i;
73
74#ifdef EPETRA_MPI
75
76 // Initialize MPI
77
78 MPI_Init(&argc,&argv);
79
80 Epetra_MpiComm Comm( MPI_COMM_WORLD );
81
82#else
83
85
86#endif
87
88 int MyPID = Comm.MyPID();
89 int NumProc = Comm.NumProc();
90 bool verbose = (MyPID==0);
91
92 if (verbose)
93 std::cout << Epetra_Version() << std::endl << std::endl;
94
95 std::cout << Comm << std::endl;
96
97 // Get the number of local equations from the command line
98 if (argc!=2)
99 {
100 if (verbose)
101 std::cout << "Usage: " << argv[0] << " number_of_equations" << std::endl;
102 std::exit(1);
103 }
104 long long NumGlobalElements = std::atoi(argv[1]);
105
106 if (NumGlobalElements < NumProc)
107 {
108 if (verbose)
109 std::cout << "numGlobalBlocks = " << NumGlobalElements
110 << " cannot be < number of processors = " << NumProc << std::endl;
111 std::exit(1);
112 }
113
114 // Construct a Map that puts approximately the same number of
115 // equations on each processor.
116
117 Epetra_Map Map(NumGlobalElements, 0LL, Comm);
118
119 // Get update list and number of local equations from newly created Map.
120
121 int NumMyElements = Map.NumMyElements();
122
123 std::vector<long long> MyGlobalElements(NumMyElements);
124 Map.MyGlobalElements(&MyGlobalElements[0]);
125
126 // Create an integer vector NumNz that is used to build the Petra Matrix.
127 // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation
128 // on this processor
129
130 std::vector<int> NumNz(NumMyElements);
131
132 // We are building a tridiagonal matrix where each row has (-1 2 -1)
133 // So we need 2 off-diagonal terms (except for the first and last equation)
134
135 for (i=0; i<NumMyElements; i++)
136 if (MyGlobalElements[i]==0 || MyGlobalElements[i] == NumGlobalElements-1)
137 NumNz[i] = 2;
138 else
139 NumNz[i] = 3;
140
141 // Create a Epetra_Matrix
142
143 Epetra_CrsMatrix A(Copy, Map, &NumNz[0]);
144
145 // Add rows one-at-a-time
146 // Need some vectors to help
147 // Off diagonal Values will always be -1
148
149
150 std::vector<double> Values(2);
151 Values[0] = -1.0; Values[1] = -1.0;
152 std::vector<long long> Indices(2);
153 double two = 2.0;
154 int NumEntries;
155
156 for (i=0; i<NumMyElements; i++)
157 {
158 if (MyGlobalElements[i]==0)
159 {
160 Indices[0] = 1;
161 NumEntries = 1;
162 }
163 else if (MyGlobalElements[i] == NumGlobalElements-1)
164 {
165 Indices[0] = NumGlobalElements-2;
166 NumEntries = 1;
167 }
168 else
169 {
170 Indices[0] = MyGlobalElements[i]-1;
171 Indices[1] = MyGlobalElements[i]+1;
172 NumEntries = 2;
173 }
174 ierr = A.InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
175 assert(ierr==0);
176 // Put in the diagonal entry
177 ierr = A.InsertGlobalValues(MyGlobalElements[i], 1, &two, &MyGlobalElements[i]);
178 assert(ierr==0);
179 }
180
181 // Finish up
182 ierr = A.FillComplete();
183 assert(ierr==0);
184
185 // Create vectors for Power method
186
187
188 // variable needed for iteration
189 double lambda = 0.0;
190 int niters = (int) NumGlobalElements*10;
191 double tolerance = 1.0e-2;
192
193 // Iterate
194 Epetra_Flops counter;
195 A.SetFlopCounter(counter);
196 Epetra_Time timer(Comm);
197 ierr += power_method(A, lambda, niters, tolerance, verbose);
198 double elapsed_time = timer.ElapsedTime();
199 double total_flops =counter.Flops();
200 double MFLOPs = total_flops/elapsed_time/1000000.0;
201
202 if (verbose)
203 std::cout << "\n\nTotal MFLOPs for first solve = " << MFLOPs << std::endl<< std::endl;
204
205 // Increase diagonal dominance
206 if (verbose)
207 std::cout << "\nIncreasing magnitude of first diagonal term, solving again\n\n"
208 << std::endl;
209
210 if (A.MyGlobalRow(0)) {
211 int numvals = A.NumGlobalEntries(0);
212 std::vector<double> Rowvals(numvals);
213 std::vector<long long> Rowinds(numvals);
214 A.ExtractGlobalRowCopy(0, numvals, numvals, &Rowvals[0], &Rowinds[0]); // Get A[0,0]
215 for (i=0; i<numvals; i++) if (Rowinds[i] == 0) Rowvals[i] *= 10.0;
216
217 A.ReplaceGlobalValues(0, numvals, &Rowvals[0], &Rowinds[0]);
218 }
219
220 // Iterate (again)
221 lambda = 0.0;
222 timer.ResetStartTime();
223 counter.ResetFlops();
224 ierr += power_method(A, lambda, niters, tolerance, verbose);
225 elapsed_time = timer.ElapsedTime();
226 total_flops = counter.Flops();
227 MFLOPs = total_flops/elapsed_time/1000000.0;
228
229 if (verbose)
230 std::cout << "\n\nTotal MFLOPs for second solve = " << MFLOPs << std::endl<< std::endl;
231
232
233 // Release all objects
234#ifdef EPETRA_MPI
235 MPI_Finalize() ;
236#endif
237
238/* end main
239*/
240return ierr ;
241}
242
243int power_method(Epetra_CrsMatrix& A, double &lambda, int niters,
244 double tolerance, bool verbose) {
245
246 Epetra_Vector q(A.RowMap());
247 Epetra_Vector z(A.RowMap());
248 Epetra_Vector resid(A.RowMap());
249
250 Epetra_Flops * counter = A.GetFlopCounter();
251 if (counter!=0) {
252 q.SetFlopCounter(A);
253 z.SetFlopCounter(A);
254 resid.SetFlopCounter(A);
255 }
256
257 // Fill z with random Numbers
258 z.Random();
259
260 // variable needed for iteration
261 double normz, residual;
262
263 int ierr = 1;
264
265 for (int iter = 0; iter < niters; iter++)
266 {
267 z.Norm2(&normz); // Compute 2-norm of z
268 q.Scale(1.0/normz, z);
269 A.Multiply(false, q, z); // Compute z = A*q
270 q.Dot(z, &lambda); // Approximate maximum eigenvalue
271 if (iter%100==0 || iter+1==niters)
272 {
273 resid.Update(1.0, z, -lambda, q, 0.0); // Compute A*q - lambda*q
274 resid.Norm2(&residual);
275 if (verbose) std::cout << "Iter = " << iter << " Lambda = " << lambda
276 << " Residual of A*q - lambda*q = "
277 << residual << std::endl;
278 }
279 if (residual < tolerance) {
280 ierr = 0;
281 break;
282 }
283 }
284 return(ierr);
285}
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)