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_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
66 int power_method(Epetra_CrsMatrix& A, double & lambda, int niters, double tolerance,
67  bool verbose);
68 
69 
70 int 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 
84  Epetra_SerialComm Comm;
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  int 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, 0, 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<int> 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<int> 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 = 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<int> 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 */
240 return ierr ;
241 }
242 
243 int 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 }
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.
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.