EpetraExt Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
hypre_Helpers.cpp
Go to the documentation of this file.
1//@HEADER
2// ***********************************************************************
3//
4// EpetraExt: Epetra Extended - 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#include "hypre_Helpers.hpp"
43#include "Teuchos_UnitTestHarness.hpp"
44#include "Teuchos_Assert.hpp"
45#include "Teuchos_Array.hpp"
46#include <iostream>
47#include <fstream>
48
49#include <math.h>
50#include <cstdlib>
51#include <time.h>
52#include <vector>
53#include <map>
54
55EpetraExt_HypreIJMatrix* newHypreMatrix(const int N, const int type)
56{
57 HYPRE_IJMatrix Matrix;
58 int ierr = 0;
59 int i;
60 int numprocs, rank;
61 MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
62 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
63
64 int ilower = (N/numprocs)*rank;
65 int iupper = (N/numprocs)*(rank+1);
66 if(rank == numprocs-1){ /*Last processor */
67 iupper = N;
68 }
69
70 // Create
71 ierr += HYPRE_IJMatrixCreate(MPI_COMM_WORLD, ilower, iupper-1, ilower, iupper-1, &Matrix);
72 ierr += HYPRE_IJMatrixSetObjectType(Matrix,HYPRE_PARCSR);
73 // Initialize
74 ierr += HYPRE_IJMatrixInitialize(Matrix);
75
76 if(Matrix == NULL){
77 printf("Error! Matrix is NULL!\n");
78 std::cin >> ierr;
79 }
80
81 if(type == 0){
82 // Set values
83 int rows[1];
84 int cols[1];
85 double values[1];
86 int ncols = 1;
87 for(i = ilower; i < iupper; i++) {
88 rows[0] = i;
89 cols[0] = i;
90 values[0] = 1.0;
91 ierr += HYPRE_IJMatrixSetValues(Matrix, 1, &ncols, rows, cols, values);
92 }
93 } else if(type == 1){
94 srand(time(NULL));
95 // Set values
96 int rows[1];
97 int cols[1];
98 double values[1];
99 int ncols = 1;
100 for(i = ilower; i < iupper; i++) {
101 rows[0] = i;
102 cols[0] = i;
103 values[0] = ( (double)rand()/(double)RAND_MAX ) * 100;
104 ierr += HYPRE_IJMatrixSetValues(Matrix, 1, &ncols, rows, cols, values);
105 }
106
107 } else if(type == 2){
108 // Set values
109 int rows[1];
110 Teuchos::Array<int> cols; cols.resize(N);
111 Teuchos::Array<double> values; values.resize(N);
112 int ncols = N;
113 for(i = ilower; i < iupper; i++) {
114 for(int j = 0; j < N; j++){
115 cols[j] = j;
116 values[j] = j;
117 }
118 rows[0] = i;
119 ierr += HYPRE_IJMatrixSetValues(Matrix, 1, &ncols, rows, &cols[0], &values[0]);
120 }
121 } else if(type == 3){
122 srand(time(NULL));
123 int rows[1];
124 Teuchos::Array<int> cols; cols.resize(N);
125 Teuchos::Array<double> values; values.resize(N);
126 int ncols = N;
127 for(i = ilower; i < iupper; i++) {
128 for(int j = 0; j < N; j++){
129 cols[j] = j;
130 double currVal = ( (double)rand()/(double)RAND_MAX ) * 100;
131 values[j] = currVal;
132 }
133 rows[0] = i;
134 ierr += HYPRE_IJMatrixSetValues(Matrix, 1, &ncols, rows, &cols[0], &values[0]);
135 }
136 } else {
137 srand(time(NULL));
138 int rows[1];
139 for(i = ilower; i < iupper; i++) {
140 int ncols = (int)(1+( (double)rand()/(double)RAND_MAX ) * 0.5*(N-1));
141 TEUCHOS_TEST_FOR_EXCEPTION(ncols <= 0, std::logic_error, "ncols is negative");
142 Teuchos::Array<int> cols; cols.resize(ncols);
143 Teuchos::Array<double> values; values.resize(ncols);
144 for(int j = 0; j < ncols; j++){
145 int index = 0;
146 if(i-(ncols/2) >= 0 && i+(ncols/2) < N){
147 index = j + i - (ncols/2);
148 } else if (i-(ncols/2) < 0){
149 index = j + i;
150 } else{
151 index = j + i - (ncols-1);
152 }
153
154 cols[j] = index;
155 double currVal = ( (double)rand()/(double)RAND_MAX ) * 100;
156 values[j] = currVal;
157 }
158 rows[0] = i;
159 ierr += HYPRE_IJMatrixSetValues(Matrix, 1, &ncols, rows, &cols[0], &values[0]);
160 }
161 }
162 // Assemble
163 ierr += HYPRE_IJMatrixAssemble(Matrix);
165 //Don't HYPRE_IJMatrixDestroy(Matrix);
166 return RetMat;
167}
168
170{
171 int N = Matrix.NumGlobalRows();
172 Epetra_CrsMatrix* TestMat = new Epetra_CrsMatrix(Copy, Matrix.RowMatrixRowMap(), Matrix.RowMatrixColMap(), N, false);
173
174 for(int i = 0; i < Matrix.NumMyRows(); i++){
175 int entries;
176 Matrix.NumMyRowEntries(i,entries);
177 Teuchos::Array<double> Values; Values.resize(entries);
178 Teuchos::Array<int> Indices; Indices.resize(entries);
179 int NumEntries;
180 Matrix.ExtractMyRowCopy(i,entries,NumEntries,&Values[0], &Indices[0]);
181 for(int j = 0; j < NumEntries; j++){
182 double currVal[1];
183 currVal[0] = Values[j];
184 int currInd[1];
185 currInd[0] = Matrix.RowMatrixColMap().GID(Indices[j]);
186 TestMat->InsertGlobalValues(Matrix.RowMatrixRowMap().GID(i), 1, currVal, currInd);
187 }
188 }
189 TestMat->FillComplete();
190 return TestMat;
191}
192
194
195 bool retVal = true;
196
197 int num_vectors = Y1.NumVectors();
198 if(Y2.NumVectors() != num_vectors){
199 printf("Multivectors do not have same number of vectors.\n");
200 return false;
201 }
202
203 for(int j = 0; j < num_vectors; j++){
204 if(Y1.MyLength() != Y2.MyLength()){
205 printf("Vectors are not same size on local processor.\n");
206 return false;
207 }
208 Teuchos::Array<double> Y1_vals; Y1_vals.resize(Y1.MyLength());
209 Teuchos::Array<double> Y2_vals; Y2_vals.resize(Y2.MyLength());
210 (*Y1(j)).ExtractCopy(&Y1_vals[0]);
211 (*Y2(j)).ExtractCopy(&Y2_vals[0]);
212
213 for(int i = 0; i < Y1.MyLength(); i++){
214 if(fabs(Y1_vals[i] - Y2_vals[i]) > tol){
215 printf("Vector number[%d] ", j);
216 printf("Val1[%d] = %f != Val2[%d] = %f\n", i, Y1_vals[i], i, Y2_vals[i]);
217 retVal = false;
218 }
219 }
220 }
221 if(retVal == false){
222 printf("Failed equivalent vectors.\n");
223 }
224 return retVal;
225}
226
227
228
229bool EquivalentMatrices(Epetra_RowMatrix &HypreMatrix, Epetra_RowMatrix &CrsMatrix, const double tol){
230 bool retVal = true;
231 for(int j = 0; j < HypreMatrix.NumMyRows(); j++){
232
233 int NumEntries;
234 int entries1;
235 int entries2;
236 HypreMatrix.NumMyRowEntries(j,NumEntries);
237
238 Teuchos::Array<double> Y1_vals; Y1_vals.resize(NumEntries);
239 Teuchos::Array<double> Y2_vals; Y2_vals.resize(NumEntries);
240 Teuchos::Array<int> indices1; indices1.resize(NumEntries);
241 Teuchos::Array<int> indices2; indices2.resize(NumEntries);
242
243 HypreMatrix.ExtractMyRowCopy(j,NumEntries, entries1, &Y1_vals[0], &indices1[0]);
244 CrsMatrix.ExtractMyRowCopy(j,NumEntries, entries2, &Y2_vals[0], &indices2[0]);
245
246 std::map<int,double> Y1map;
247 std::map<int,double> Y2map;
248 for (int i=0; i < NumEntries ; ++i) {
249 Y1map[indices1[i]] = Y1_vals[i];
250 Y2map[indices2[i]] = Y2_vals[i];
251 }
252 retVal = retVal && (Y1map == Y2map);
253 }
254 Teuchos::Array<int> vals; vals.resize(HypreMatrix.Comm().NumProc());
255 int my_vals[1]; my_vals[0] = (int)retVal;
256 HypreMatrix.Comm().GatherAll(my_vals, &vals[0], 1);
257 for(int i = 0; i < HypreMatrix.Comm().NumProc(); i++){
258 if(vals[i] == false){
259 retVal = false;
260 }
261 }
262 if(retVal == false){
263 printf("[%d]Failed matrix equivalency test.\n", HypreMatrix.Comm().MyPID());
264 if(HypreMatrix.Comm().MyPID() == 0){
265 //std::ofstream outfile("Matrices.txt");
266
267 //HypreMatrix.Print(outfile);
268 //CrsMatrix.Print(outfile);
269 //outfile.close();
270 }
271 }
272 return retVal;
273}
Copy
int NumMyRowEntries(int MyRow, int &NumEntries) const
Return the current number of values stored for the specified local row.
int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
Returns a copy of the specified local row in user-provided arrays.
virtual const Epetra_Map & RowMatrixColMap() const
virtual int NumMyRows() const
virtual int NumGlobalRows() const
virtual const Epetra_Map & RowMatrixRowMap() const
int GID(int LID) const
virtual int NumProc() const=0
virtual int GatherAll(double *MyVals, double *AllVals, int Count) const=0
virtual int MyPID() const=0
int FillComplete(bool OptimizeDataStorage=true)
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
int NumVectors() const
int MyLength() const
virtual const Epetra_Comm & Comm() const=0
virtual int NumMyRows() const=0
virtual int NumMyRowEntries(int MyRow, int &NumEntries) const=0
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
bool EquivalentVectors(Epetra_MultiVector &Y1, Epetra_MultiVector &Y2, const double tol)
Epetra_CrsMatrix * newCrsMatrix(EpetraExt_HypreIJMatrix &Matrix)
EpetraExt_HypreIJMatrix * newHypreMatrix(const int N, const int type)
bool EquivalentMatrices(Epetra_RowMatrix &HypreMatrix, Epetra_RowMatrix &CrsMatrix, const double tol)