Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
main.cc
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 <stdio.h>
44#include <stdlib.h>
45#include <ctype.h>
46#include <assert.h>
47#include <string.h>
48#include <math.h>
49#include "Petra_Comm.h"
50#include "Petra_Map.h"
51#include "Petra_RDP_MultiVector.h"
52#include "Petra_RDP_Vector.h"
53#include "Petra_RDP_DCRS_Matrix.h"
54#ifdef PETRA_MPI
55#include "mpi.h"
56#endif
57#ifndef __cplusplus
58#define __cplusplus
59#endif
60
61int main(int argc, char *argv[])
62{
63 int i;
64
65#ifdef PETRA_MPI
66 MPI_Init(&argc,&argv);
67#endif
68
69 // get number of processors and the name of this processor
70
71#ifdef PETRA_MPI
72 Petra_Comm& comm = *new Petra_Comm(MPI_COMM_WORLD);
73#else
74 Petra_Comm& comm = *new Petra_Comm();
75#endif
76
77 int NumProc = comm.getNumProc();
78 int MyPID = comm.getMyPID();
79 cout << "Processor " << MyPID << " of " << NumProc << " is alive." << endl;
80
81 // Get the number of local equations from the command line
82 if (argc!=2)
83 {
84 if (MyPID==0) cout << "Usage: " << argv[0] << " number_of_equations" << endl;
85 exit(1);
86 }
87 int numGlobalEquations = atoi(argv[1]);
88
89 if (numGlobalEquations < NumProc)
90 {
91 if (MyPID==0)
92 cout << "numGlobalBlocks = " << numGlobalEquations
93 << " cannot be < number of processors = " << NumProc << endl;
94 exit(1);
95 }
96
97 // Construct a map that puts approximately the same number of equations on each processor
98
99 Petra_Map& map = *new Petra_Map(numGlobalEquations, comm);
100
101 // Get update list and number of local equations from newly created map
102 int * UpdateList = map.getUpdateList();
103 int numLocalEquations = map.numLocalEquations();
104
105 // Create an integer vector numNz that is used to build the Petra Matrix.
106 // numNz[i] is the number of OFF-DIAGONAL term for the ith global equation on this processor
107
108 int * numNz = new int[numLocalEquations];
109
110 // We are building a tridiagonal matrix where each row has (-1 2 -1)
111 // So we need 2 off-diagonal terms (except for the first and last equation)
112
113 for (i=0; i<numLocalEquations; i++)
114 if (UpdateList[i]==0 || UpdateList[i] == numGlobalEquations-1)
115 numNz[i] = 1;
116 else
117 numNz[i] = 2;
118
119 // Create a Petra_Matrix
120
121 Petra_RDP_DCRS_Matrix& A = *new Petra_RDP_DCRS_Matrix(map);
122
123 // Allocate space using numNz
124
125 assert(A.allocate(numNz)==0);
126
127 // Add rows one-at-a-time
128 // Need some vectors to help
129 // Off diagonal values will always be -1
130
131
132 double *values = new double[2];
133 values[0] = -1.0; values[1] = -1.0;
134 int *indices = new int[2];
135 double two = 2.0;
136 int numEntries;
137
138 for (i=0; i<numLocalEquations; i++)
139 {
140 if (UpdateList[i]==0)
141 {
142 indices[0] = 1;
143 numEntries = 1;
144 }
145 else if (UpdateList[i] == numGlobalEquations-1)
146 {
147 indices[0] = numGlobalEquations-2;
148 numEntries = 1;
149 }
150 else
151 {
152 indices[0] = UpdateList[i]-1;
153 indices[1] = UpdateList[i]+1;
154 numEntries = 2;
155 }
156 assert(A.putRow(UpdateList[i], numEntries, values, indices)==0);
157 assert(A.putRow(UpdateList[i], 1, &two, UpdateList+i)==0); // Put in the diagonal entry
158 }
159
160 // Finish up
161 assert(A.fillComplete()==0);
162
163 // Create vectors for Power method
164
165 Petra_RDP_Vector& q = *new Petra_RDP_Vector(map);
166 Petra_RDP_Vector& z = *new Petra_RDP_Vector(map);
167 Petra_RDP_Vector& resid = *new Petra_RDP_Vector(map);
168
169 // Fill z with random numbers
170 z.random();
171
172 // variable needed for iteration
173 double normz, lambda, residual;
174
175 // Iterate
176 int niters = 500*numGlobalEquations;
177 double tolerance = 1.0e-10;
178 for (int iter = 0; iter < niters; iter++)
179 {
180 z.norm2(&normz); // Compute 2-norm of z
181 q.scaleCopy(z, 1.0/normz);
182 A.matvec(q, z); // Compute z = A*q
183 q.dotProd(z, &lambda); // Approximate maximum eigenvaluE
184 if (iter%100==0 || iter+1==niters)
185 {
186 resid.linComb(z, -lambda, q); // Compute A*q - lambda*q
187 resid.norm2(&residual);
188 if (MyPID==0) cout << "Iter = " << iter << " Lambda = " << lambda
189 << " Residual of A*q - lambda*q = " << residual << endl;
190 }
191 if (residual < tolerance) break;
192 }
193
194 // Release all objects
195
196 delete [] numNz;
197 delete [] values;
198 delete [] indices;
199
200 delete &resid;
201 delete &z;
202 delete &q;
203 delete &A;
204 delete &map;
205 delete &comm;
206
207#ifdef PETRA_MPI
208 MPI_Finalize() ;
209#endif
210
211/* end main
212*/
213return 0 ;
214}
int main(int argc, char *argv[])
Definition main.cc:61