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