Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 }
void SetFlopCounter(const Epetra_Flops &FlopCounter_in)
Set the internal Epetra_Flops() pointer.
bool MyGlobalRow(int GID) const
Returns true of GID is owned by the calling processor, otherwise it returns false.
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:127
int NumGlobalEntries(long long Row) const
Returns the current number of nonzero entries in specified global row on this processor.
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 MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
double ElapsedTime(void) const
Epetra_Time elapsed time function.
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:430
int power_method(Epetra_CrsMatrix &A, double &lambda, int niters, double tolerance, bool verbose)
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
int MyPID() const
Return my process ID.
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:83
const Epetra_Map & RowMap() const
Returns the Epetra_Map object associated with the rows of this matrix.
int NumMyElements() const
Number of elements on the calling processor.
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 NumProc() const
Returns total number of processes (always returns 1 for SerialComm).
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:85
Epetra_Flops: The Epetra Floating Point Operations Class.
Definition: Epetra_Flops.h:66
double Flops() const
Returns the number of floating point operations with this object and resets the count.
Definition: Epetra_Flops.h:82
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
int main(int argc, char *argv[])
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...
Epetra_Flops * GetFlopCounter() const
Get the pointer to the Epetra_Flops() object associated with this object, returns 0 if none...