Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
test/Bugs_LL/Bug_5797_Sort_MakeColMap_LL/cxx_main.cpp
Go to the documentation of this file.
1// Program for testing Epetra64 implementation.
2// Builds a Laplacian matrx using either Epetra or Epetra64.
3// Multiplies it by the vector of all ones and checks norms.
4//
5// To run:
6// mpirun -np # test64.exe [n]
7// where n is an optional argument specifying the number of matrix rows.
8// Default n == 25.
9//
10// Two macro definitions below control behavior:
11// ITYPE: int --> use Epetra
12// long long --> use Epetra64
13// OFFSET_EPETRA64: Add this value to each row/column index. Resulting
14// matrix rows/columns are indexed from
15// OFFSET_EPETRA64 to OFFSET_EPETRA64+n-1.
16
17#include <limits.h>
18#define ITYPE long long
19//#define OFFSET_EPETRA64 ((long long)(2)*(long long)INT_MAX)
20//#define OFFSET_EPETRA64 (INT_MAX-5)
21#define OFFSET_EPETRA64 0
22
23#include <stdio.h>
24
25#include "Epetra_ConfigDefs.h"
26
27#ifdef EPETRA_MPI
28#include <mpi.h>
29#include "Epetra_MpiComm.h"
30#define FINALIZE MPI_Finalize()
31#else
32#include "Epetra_SerialComm.h"
33#define FINALIZE
34#endif
35
36#include "Epetra_CrsMatrix.h"
37#include "Epetra_Map.h"
38#include "Epetra_Vector.h"
39
41// Build the following Laplacian matrix using Epetra or Epetra64
42//
43// | 2 -1 ... -1... |
44// |-1 3 -1 ... -1... |
45// | ... |
46// |-1 -1 4 -1... -1... |
47// | ... |
48// | -1 -1 3 -1|
49// | -1 -1 2|
50//
52
53#define MIN(a,b) ((a) < (b) ? (a) : (b))
54
55int main(int narg, char *arg[])
56{
57 using std::cout;
58
59#ifdef EPETRA_MPI
60 // Initialize MPI
61 MPI_Init(&narg,&arg);
62 Epetra_MpiComm comm(MPI_COMM_WORLD);
63#else
65#endif
66
67 int me = comm.MyPID();
68 int np = comm.NumProc();
69
70 ITYPE nGlobalRows = 25;
71 if (narg > 1)
72 nGlobalRows = (ITYPE) atol(arg[1]);
73
74 bool verbose = (nGlobalRows < 26);
75
76 // Linear map similar to Trilinos default,
77 // but want to allow adding OFFSET_EPETRA64 to the indices.
78 int nMyRows = (int) (nGlobalRows / np + (nGlobalRows % np > me));
79 ITYPE myFirstRow = (ITYPE)(me * (nGlobalRows / np) + MIN(nGlobalRows%np, me));
80 ITYPE *myGlobalRows = new ITYPE[nMyRows];
81 for (int i = 0; i < nMyRows; i++)
82 myGlobalRows[i] = (ITYPE)i + myFirstRow + OFFSET_EPETRA64;
83 Epetra_Map *rowMap = new Epetra_Map(-1, nMyRows, &myGlobalRows[0], 0, comm);
84 if (verbose) rowMap->Print(std::cout);
85
86 // Create an integer vector nnzPerRow that is used to build the Epetra Matrix.
87 // nnzPerRow[i] is the number of entries for the ith local equation
88 std::vector<int> nnzPerRow(nMyRows+1, 0);
89
90 // Also create lists of the nonzeros to be assigned to processors.
91 // To save programming time and complexity, these vectors are allocated
92 // bigger than they may actually be needed.
93 std::vector<ITYPE> iv(5*nMyRows+1);
94 std::vector<ITYPE> jv(5*nMyRows+1);
95 std::vector<double> vv(5*nMyRows+1);
96
97 // Generate the nonzeros for the Laplacian matrix.
98 ITYPE nMyNonzeros = 0;
99 for (ITYPE i = 0, myrowcnt = 0; i < nGlobalRows; i++) {
100 if (rowMap->MyGID(i+OFFSET_EPETRA64)) {
101 int idegree = 0;
102 // This processor owns this row; add nonzeros.
103 if (i > 0) {
104 iv[nMyNonzeros] = i + OFFSET_EPETRA64;
105 jv[nMyNonzeros] = i-1 + OFFSET_EPETRA64;
106 vv[nMyNonzeros] = -1;
107 if (verbose)
108 std::cout << "(" << iv[nMyNonzeros] << "," << jv[nMyNonzeros] << ")="
109 << vv[nMyNonzeros] << " on processor " << me
110 << " in " << myrowcnt << std::endl;
111 nMyNonzeros++;
112 nnzPerRow[myrowcnt]++;
113 idegree++;
114 }
115
116 if (i < nGlobalRows - 1) {
117 iv[nMyNonzeros] = i + OFFSET_EPETRA64;
118 jv[nMyNonzeros] = i+1 + OFFSET_EPETRA64;
119 vv[nMyNonzeros] = -1;
120 if (verbose)
121 std::cout << "(" << iv[nMyNonzeros] << "," << jv[nMyNonzeros] << ")="
122 << vv[nMyNonzeros] << " on processor " << me
123 << " in " << myrowcnt << std::endl;
124 nMyNonzeros++;
125 nnzPerRow[myrowcnt]++;
126 idegree++;
127 }
128
129 if (i+5 < nGlobalRows) {
130 iv[nMyNonzeros] = i + OFFSET_EPETRA64;
131 jv[nMyNonzeros] = i+5 + OFFSET_EPETRA64;
132 vv[nMyNonzeros] = -1;
133 if (verbose)
134 std::cout << "(" << iv[nMyNonzeros] << "," << jv[nMyNonzeros] << ")="
135 << vv[nMyNonzeros] << " on processor " << me
136 << " in " << myrowcnt << std::endl;
137 nMyNonzeros++;
138 nnzPerRow[myrowcnt]++;
139 idegree++;
140 }
141
142 if (i-5 >= 0) {
143 iv[nMyNonzeros] = i + OFFSET_EPETRA64;
144 jv[nMyNonzeros] = i-5 + OFFSET_EPETRA64;
145 vv[nMyNonzeros] = -1;
146 if (verbose)
147 std::cout << "(" << iv[nMyNonzeros] << "," << jv[nMyNonzeros] << ")="
148 << vv[nMyNonzeros] << " on processor " << me
149 << " in " << myrowcnt << std::endl;
150 nMyNonzeros++;
151 nnzPerRow[myrowcnt]++;
152 idegree++;
153 }
154
155 iv[nMyNonzeros] = i + OFFSET_EPETRA64;
156 jv[nMyNonzeros] = i + OFFSET_EPETRA64;
157 vv[nMyNonzeros] = idegree;
158 if (verbose)
159 std::cout << "(" << iv[nMyNonzeros] << "," << jv[nMyNonzeros] << ")="
160 << vv[nMyNonzeros] << " on processor " << me
161 << " in " << myrowcnt << std::endl;
162 nMyNonzeros++;
163 nnzPerRow[myrowcnt]++;
164
165 myrowcnt++;
166 }
167 }
168
169 // Create an Epetra_Matrix
170 // TODO: Should be able to use StaticProfile == true, but there appears to
171 // TODO: be a bug in Epetra64 (or in my usage of it). I filed bug 5791
172 // TODO: and will use StaticProfile == false for now. When the bug is fixed,
173 // TODO: we should switch back to StaticProfile == true, as it will be more
174 // TODO: efficient.
175 // #warning Once bug 5791 is fixed, change StaticProfile back to true.
176 // -- Done. Chetan.
177 Epetra_CrsMatrix *A = new Epetra_CrsMatrix(Copy, *rowMap, &nnzPerRow[0], true);
178
179 // Insert the nonzeros.
180#ifndef NDEBUG
181 int info;
182#endif
183 ITYPE sum = 0;
184 for (int i=0; i < nMyRows; i++) {
185 if (nnzPerRow[i]) {
186 if (verbose) {
187 std::cout << "InsertGlobalValus row " << iv[sum]
188 << " count " << nnzPerRow[i]
189 << " cols " << jv[sum] << " " << jv[sum+1] << " ";
190 if (nnzPerRow[i] == 3) std::cout << jv[sum+2];
191 std::cout << std::endl;
192 }
193#ifndef NDEBUG
194 info =
195#endif
196 A->InsertGlobalValues(iv[sum],nnzPerRow[i],&vv[sum],&jv[sum]);
197 assert(info==0);
198 sum += nnzPerRow[i];
199 }
200 }
201
202 // Finish up
203#ifndef NDEBUG
204 info =
205#endif
206 A->FillComplete();
207 assert(info==0);
208 if (verbose) A->Print(std::cout);
209
210 // Sanity test: Product of matrix and vector of ones should have norm == 0
211 // and max/min/mean values of 0
212 Epetra_Vector sanity(A->RangeMap());
213 Epetra_Vector sanityres(A->DomainMap());
214 sanity.PutScalar(1.);
215 A->Multiply(false, sanity, sanityres);
216
217 double jjone, jjtwo, jjmax;
218 sanityres.Norm1(&jjone);
219 sanityres.Norm2(&jjtwo);
220 sanityres.NormInf(&jjmax);
221 if (me == 0)
222 std::cout << "SanityTest norms 1/2/inf: " << jjone << " "
223 << jjtwo << " " << jjmax << std::endl;
224
225 bool test_failed = (jjone != 0) || (jjtwo != 0) || (jjmax != 0);
226
227 sanityres.MinValue(&jjone);
228 sanityres.MeanValue(&jjtwo);
229 sanityres.MaxValue(&jjmax);
230 if (me == 0)
231 std::cout << "SanityTest values min/max/avg: " << jjone << " "
232 << jjmax << " " << jjtwo << std::endl;
233
234 test_failed = test_failed || (jjone != 0) || (jjtwo != 0) || (jjmax != 0);
235
236 if (me == 0) {
237 if(test_failed)
238 std::cout << "Bug_5797_Sort_MakeColMap_LL tests FAILED" << std::endl;
239 }
240
241 delete A;
242 delete rowMap;
243 delete [] myGlobalRows;
244
245 FINALIZE;
246}
247
virtual void Print(std::ostream &os) const
Print object to an output stream.
bool MyGID(int GID_in) const
Returns true if the GID passed in belongs to the calling processor in this map, otherwise returns fal...
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
int FillComplete(bool OptimizeDataStorage=true)
Signal that data entry is complete. Perform transformations to local index space.
virtual void Print(std::ostream &os) const
Print method.
const Epetra_Map & RangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
const Epetra_Map & DomainMap() const
Returns the Epetra_Map object associated with the domain of this matrix operator.
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_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 NormInf(double *Result) const
Compute Inf-norm of each vector in multi-vector.
int MinValue(double *Result) const
Compute minimum value of each vector in multi-vector.
int Norm1(double *Result) const
Compute 1-norm of each vector in multi-vector.
int MeanValue(double *Result) const
Compute mean (average) value of each vector in multi-vector.
int Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
int MaxValue(double *Result) const
Compute maximum value of each vector in multi-vector.
int PutScalar(double ScalarConstant)
Initialize all values in a multi-vector with constant value.
Epetra_SerialComm: The Epetra Serial Communication Class.
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
int main(int narg, char *arg[])