Epetra Package Browser (Single Doxygen Collection)  Development
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
67 int power_method(Epetra_CrsMatrix& A, double & lambda, int niters, double tolerance,
68  bool verbose);
69 
70 
71 int 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 
85  Epetra_SerialComm Comm;
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 */
241 return ierr ;
242 }
243 
244 int 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 }
const Epetra_Map & RowMap() const
Returns the Epetra_Map object associated with the rows of this matrix.
void SetFlopCounter(const Epetra_Flops &FlopCounter_in)
Set the internal Epetra_Flops() pointer.
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:119
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.
int main(int argc, char *argv[])
Epetra_Flops * GetFlopCounter() const
Get the pointer to the Epetra_Flops() object associated with this object, returns 0 if none...
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.
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...
Definition: Epetra_Util.h:413
int power_method(Epetra_CrsMatrix &A, double &lambda, int niters, double tolerance, bool verbose)
bool MyGlobalRow(int GID) const
Returns true of GID is owned by the calling processor, otherwise it returns false.
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
Epetra_MpiComm: The Epetra MPI Communication Class.
std::string Epetra_Version()
int FillComplete(bool OptimizeDataStorage=true)
Signal that data entry is complete. Perform transformations to local index space. ...
Epetra_Time: The Epetra Timing Class.
Definition: Epetra_Time.h:75
double Flops() const
Returns the number of floating point operations with this object and resets the count.
Definition: Epetra_Flops.h:74
int NumProc() const
Returns total number of processes (always returns 1 for SerialComm).
int NumMyElements() const
Number of elements on the calling processor.
int NumGlobalEntries(long long Row) const
Returns the current number of nonzero entries in specified global row on this processor.
Epetra_SerialComm: The Epetra Serial Communication Class.
void ResetFlops()
Resets the number of floating point operations to zero for this multi-vector.
Definition: Epetra_Flops.h:77
Epetra_Flops: The Epetra Floating Point Operations Class.
Definition: Epetra_Flops.h:58
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 MyPID() const
Return my process ID.
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
void ResetStartTime(void)
Epetra_Time function to reset the start time for a timer object.
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...
double ElapsedTime(void) const
Epetra_Time elapsed time function.