Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
test/Bugs_LL/Bug_5794_IndexBase_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 == 10.
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
22#include <stdio.h>
23
24#include "Epetra_ConfigDefs.h"
25
26#ifdef EPETRA_MPI
27#include <mpi.h>
28#include "Epetra_MpiComm.h"
29#define FINALIZE MPI_Finalize()
30#else
31#include "Epetra_SerialComm.h"
32#define FINALIZE
33#endif
34
35#include "Epetra_CrsMatrix.h"
36#include "Epetra_Map.h"
37#include "Epetra_Vector.h"
38
40// Build the following Laplacian matrix using Epetra or Epetra64.
41//
42// | 1 -1 |
43// |-1 2 -1 |
44// | -1 2 -1 |
45// | ... |
46// | -1 1|
47//
49
50#define MIN(a,b) ((a) < (b) ? (a) : (b))
51
52int main(int narg, char *arg[])
53{
54 using std::cout;
55
56#ifdef EPETRA_MPI
57 // Initialize MPI
58 MPI_Init(&narg,&arg);
59 Epetra_MpiComm comm(MPI_COMM_WORLD);
60#else
62#endif
63
64 int me = comm.MyPID();
65 int np = comm.NumProc();
66
67 ITYPE nGlobalRows = 10;
68 if (narg > 1)
69 nGlobalRows = (ITYPE) atol(arg[1]);
70
71 bool verbose = (nGlobalRows < 20);
72
73 // Linear map similar to Trilinos default,
74 // but want to allow adding OFFSET_EPETRA64 to the indices.
75 int nMyRows = (int) (nGlobalRows / np + (nGlobalRows % np > me));
76 ITYPE myFirstRow = (ITYPE)(me * (nGlobalRows / np) + MIN(nGlobalRows%np, me));
77 ITYPE *myGlobalRows = new ITYPE[nMyRows];
78 for (int i = 0; i < nMyRows; i++)
79 myGlobalRows[i] = (ITYPE)i + myFirstRow + OFFSET_EPETRA64;
80 Epetra_Map *rowMap = new Epetra_Map(-1, nMyRows, &myGlobalRows[0], 0, comm);
81 if (verbose) rowMap->Print(std::cout);
82
83 // Create an integer vector nnzPerRow that is used to build the Epetra Matrix.
84 // nnzPerRow[i] is the number of entries for the ith local equation
85 std::vector<int> nnzPerRow(nMyRows+1, 0);
86
87 // Also create lists of the nonzeros to be assigned to processors.
88 // To save programming time and complexity, these vectors are allocated
89 // bigger than they may actually be needed.
90 std::vector<ITYPE> iv(3*nMyRows+1);
91 std::vector<ITYPE> jv(3*nMyRows+1);
92 std::vector<double> vv(3*nMyRows+1);
93
94 // Generate the nonzeros for the Laplacian matrix.
95 ITYPE nMyNonzeros = 0;
96 for (ITYPE i = 0, myrowcnt = 0; i < nGlobalRows; i++) {
97 if (rowMap->MyGID(i+OFFSET_EPETRA64)) {
98 // This processor owns this row; add nonzeros.
99 if (i > 0) {
100 iv[nMyNonzeros] = i + OFFSET_EPETRA64;
101 jv[nMyNonzeros] = i-1 + OFFSET_EPETRA64;
102 vv[nMyNonzeros] = -1;
103 if (verbose)
104 std::cout << "(" << iv[nMyNonzeros] << "," << jv[nMyNonzeros] << ")="
105 << vv[nMyNonzeros] << " on processor " << me
106 << " in " << myrowcnt << std::endl;
107 nMyNonzeros++;
108 nnzPerRow[myrowcnt]++;
109 }
110
111 iv[nMyNonzeros] = i + OFFSET_EPETRA64;
112 jv[nMyNonzeros] = i + OFFSET_EPETRA64;
113 vv[nMyNonzeros] = ((i == 0 || i == nGlobalRows-1) ? 1. : 2.);
114 if (verbose)
115 std::cout << "(" << iv[nMyNonzeros] << "," << jv[nMyNonzeros] << ")="
116 << vv[nMyNonzeros] << " on processor " << me
117 << " in " << myrowcnt << std::endl;
118 nMyNonzeros++;
119 nnzPerRow[myrowcnt]++;
120
121 if (i < nGlobalRows - 1) {
122 iv[nMyNonzeros] = i + OFFSET_EPETRA64;
123 jv[nMyNonzeros] = i+1 + OFFSET_EPETRA64;
124 vv[nMyNonzeros] = -1;
125 if (verbose)
126 std::cout << "(" << iv[nMyNonzeros] << "," << jv[nMyNonzeros] << ")="
127 << vv[nMyNonzeros] << " on processor " << me
128 << " in " << myrowcnt << std::endl;
129 nMyNonzeros++;
130 nnzPerRow[myrowcnt]++;
131 }
132 myrowcnt++;
133 }
134 }
135
136 // Create an Epetra_Matrix
137 Epetra_CrsMatrix *A = new Epetra_CrsMatrix(Copy, *rowMap, &nnzPerRow[0], false);
138
139 // Insert the nonzeros.
140#ifndef NDEBUG
141 int info;
142#endif
143 ITYPE sum = 0;
144 for (int i=0; i < nMyRows; i++) {
145 if (nnzPerRow[i]) {
146 if (verbose) {
147 std::cout << "InsertGlobalValus row " << iv[sum]
148 << " count " << nnzPerRow[i]
149 << " cols " << jv[sum] << " " << jv[sum+1] << " ";
150 if (nnzPerRow[i] == 3) std::cout << jv[sum+2];
151 std::cout << std::endl;
152 }
153#ifndef NDEBUG
154 info =
155#endif
156 A->InsertGlobalValues(iv[sum],nnzPerRow[i],&vv[sum],&jv[sum]);
157 assert(info==0);
158 sum += nnzPerRow[i];
159 }
160 }
161
162 // Finish up
163#ifndef NDEBUG
164 info =
165#endif
166 A->FillComplete();
167 assert(info==0);
168 if (verbose) A->Print(std::cout);
169
170 // Sanity test: Product of matrix and vector of ones should have norm == 0
171 // and max/min/mean values of 0
172 Epetra_Vector sanity(A->RangeMap());
173 Epetra_Vector sanityres(A->DomainMap());
174 sanity.PutScalar(1.);
175 A->Multiply(false, sanity, sanityres);
176
177 double jjone, jjtwo, jjmax;
178 sanityres.Norm1(&jjone);
179 sanityres.Norm2(&jjtwo);
180 sanityres.NormInf(&jjmax);
181 if (me == 0)
182 std::cout << "SanityTest norms 1/2/inf: " << jjone << " "
183 << jjtwo << " " << jjmax << std::endl;
184
185 bool test_failed = (jjone != 0) || (jjtwo != 0) || (jjmax != 0);
186
187 sanityres.MinValue(&jjone);
188 sanityres.MeanValue(&jjtwo);
189 sanityres.MaxValue(&jjmax);
190 if (me == 0)
191 std::cout << "SanityTest values min/max/avg: " << jjone << " "
192 << jjmax << " " << jjtwo << std::endl;
193
194 test_failed = test_failed || (jjone != 0) || (jjtwo != 0) || (jjmax != 0);
195
196 if (me == 0) {
197 if(test_failed)
198 std::cout << "Bug_5794_IndexBase_LL tests FAILED" << std::endl;
199 }
200
201 delete A;
202 delete rowMap;
203 delete [] myGlobalRows;
204
205 FINALIZE;
206}
207
@ Copy
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[])