Ifpack 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// Didasko Tutorial Package
5// Copyright (2005) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
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 about Didasko? Contact Marzio Sala (marzio.sala _AT_ gmail.com)
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 "Teuchos_RCP.hpp"
47#include <iostream>
48#include <fstream>
49#include "Galeri_Maps.h"
50#include "Galeri_CrsMatrices.h"
51#include "Galeri_Utils.h"
52#include "Ifpack_Hypre.h"
53
54#include <math.h>
55#include <stdlib.h>
56#include <time.h>
57#include <vector>
58#include <sstream>
59#include <map>
60
61using Teuchos::RCP;
62using Teuchos::rcp;
64{
65 HYPRE_IJMatrix Matrix;
66 int ierr = 0;
67 int i;
68 int numprocs, rank;
69 MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
70 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
71
72 int ilower = (N/numprocs)*rank;
73 int iupper = (N/numprocs)*(rank+1);
74 if(rank == numprocs-1){ /*Last processor */
75 iupper = N;
76 }
77
78 // Create
79 ierr += HYPRE_IJMatrixCreate(MPI_COMM_WORLD, ilower, iupper-1, ilower, iupper-1, &Matrix);
80 ierr += HYPRE_IJMatrixSetObjectType(Matrix,HYPRE_PARCSR);
81 // Initialize
82 ierr += HYPRE_IJMatrixInitialize(Matrix);
83
84 if(Matrix == NULL){
85 printf("Error! Matrix is NULL!\n");
86 std::cin >> ierr;
87 }
88
89 srand(time(NULL));
90 int rows[1];
91 for(i = ilower; i < iupper; i++) {
92 int ncols = (int)(1+( (double)rand()/(double)RAND_MAX ) * 0.5*(N-1));
93 TEUCHOS_TEST_FOR_EXCEPTION(ncols <= 0, std::logic_error, "ncols is negative");
94 Teuchos::Array<int> cols; cols.resize(ncols);
95 Teuchos::Array<double> values; values.resize(ncols);
96 for(int j = 0; j < ncols; j++){
97 int index = 0;
98 if(i-(ncols/2) >= 0 && i+(ncols/2) < N){
99 index = j + i - (ncols/2);
100 } else if (i-(ncols/2) < 0){
101 index = j + i;
102 } else{
103 index = j + i - (ncols-1);
104 }
105
106 cols[j] = index;
107 double currVal = ( (double)rand()/(double)RAND_MAX ) * 100;
108 values[j] = currVal;
109 }
110 rows[0] = i;
111 ierr += HYPRE_IJMatrixAddToValues(Matrix, 1, &ncols, rows, &cols[0], &values[0]);
112 }
113 // Assemble
114 ierr += HYPRE_IJMatrixAssemble(Matrix);
116 //Don't HYPRE_IJMatrixDestroy(Matrix);
117 return RetMat;
118}
119
121
122 Epetra_MpiComm Comm(MPI_COMM_WORLD);
123 // pointer to the matrix to be created
124 Epetra_CrsMatrix* Matrix;
125 //Epetra_CrsMatrix* Matrix;
126 // container for parameters
127 Teuchos::ParameterList GaleriList;
128 int nx = N * Comm.NumProc();
129 int ny = N * Comm.NumProc();
130 GaleriList.set("nx", nx);
131 GaleriList.set("ny", ny);
132
133 // amk Nov 24, 2015: This map described a non-contiguous distribution before.
134 // hypre didn't like that at all, so I changed it
135 Epetra_Map Map(nx*ny,0,Comm);
136
137 Matrix = Galeri::CreateCrsMatrix("Laplace2D", &Map, GaleriList);
138 return Matrix;
139}
140
142{
143 int N = Matrix->NumGlobalRows();
144 Epetra_CrsMatrix* TestMat = new Epetra_CrsMatrix(Copy, Matrix->RowMatrixRowMap(), Matrix->RowMatrixColMap(), N, false);
145
146 for(int i = 0; i < Matrix->NumMyRows(); i++){
147 int entries;
148 Matrix->NumMyRowEntries(i,entries);
149 Teuchos::Array<double> Values; Values.resize(entries);
150 Teuchos::Array<int> Indices; Indices.resize(entries);
151 int NumEntries;
152 Matrix->ExtractMyRowCopy(i,entries,NumEntries,&Values[0], &Indices[0]);
153 for(int j = 0; j < NumEntries; j++){
154 double currVal[1];
155 currVal[0] = Values[j];
156 int currInd[1];
157 currInd[0] = Matrix->RowMatrixColMap().GID(Indices[j]);
158 TestMat->InsertGlobalValues(Matrix->RowMatrixRowMap().GID(i), 1, currVal, currInd);
159 }
160 }
161 TestMat->FillComplete();
162 return TestMat;
163}
164
166
167 bool retVal = true;
168
169 int num_vectors = Y1.NumVectors();
170 if(Y2.NumVectors() != num_vectors){
171 printf("Multivectors do not have same number of vectors.\n");
172 return false;
173 }
174
175 for(int j = 0; j < num_vectors; j++){
176 if(Y1.MyLength() != Y2.MyLength()){
177 printf("Vectors are not same size on local processor.\n");
178 return false;
179 }
180 for(int i = 0; i < Y1.GlobalLength(); i++){
181 int Y1_local = Y1.Map().LID(i);
182 int Y2_local = Y2.Map().LID(i);
183 if(Y1_local < 0 || Y2_local < 0){
184 continue;
185 }
186 if(fabs((*Y1(j))[Y1_local] - (*Y2(j))[Y2_local]) > tol){
187 printf("Vector number[%d] ", j);
188 printf("Val1[%d] = %f != Val2[%d] = %f\n", i, (*Y1(j))[Y1_local], i, (*Y2(j))[Y2_local]);
189 retVal = false;
190 }
191 }
192 }
193 Teuchos::Array<int> vals; vals.resize(Y1.Comm().NumProc());
194 int my_vals[1]; my_vals[0] = (int)retVal;
195 Y1.Comm().GatherAll(my_vals, &vals[0], 1);
196 for(int i = 0; i < Y1.Comm().NumProc(); i++){
197 if(vals[i] == false){
198 retVal = false;
199 }
200 }
201 if(retVal == false){
202 printf("[%d]Failed vector equivalency test.\n", Y1.Comm().MyPID());
203 }
204 return retVal;
205}
206
207
208
209bool EquivalentMatrices(Epetra_RowMatrix &HypreMatrix, Epetra_RowMatrix &CrsMatrix, const double tol){
210 bool retVal = true;
211 int MyPID = HypreMatrix.Comm().MyPID();
212 if(HypreMatrix.NumMyRows() != CrsMatrix.NumMyRows()){
213 printf("Different number of local rows.");
214 return false;
215 }
216 for(int j = 0; j < HypreMatrix.NumGlobalRows(); j++){
217 int hyp_j = HypreMatrix.RowMatrixRowMap().LID(j);
218 int crs_j = CrsMatrix.RowMatrixRowMap().LID(j);
219 if(hyp_j < 0 || crs_j < 0){
220 continue;
221 }
222
223 int NumEntries = HypreMatrix.NumMyCols();
224 int entries1;
225 int entries2;
226
227 Teuchos::Array<double> Y1_vals; Y1_vals.resize(NumEntries);
228 Teuchos::Array<double> Y2_vals; Y2_vals.resize(NumEntries);
229 Teuchos::Array<int> indices1; indices1.resize(NumEntries);
230 Teuchos::Array<int> indices2; indices2.resize(NumEntries);
231
232 HypreMatrix.ExtractMyRowCopy(hyp_j,NumEntries, entries1, &Y1_vals[0], &indices1[0]);
233 CrsMatrix.ExtractMyRowCopy(crs_j,NumEntries, entries2, &Y2_vals[0], &indices2[0]);
234 if(entries1 != entries2){
235 printf("Row[%d], Differing number of entries %d != %d\n", j, entries1, entries2);
236 }
237 for(int i = 0; i < entries1; i++){
238 indices1[i] = HypreMatrix.RowMatrixColMap().GID(indices1[i]);
239 indices2[i] = CrsMatrix.RowMatrixColMap().GID(indices2[i]);
240 }
241 for(int i = 1; i < entries1; ++i){
242 int value = indices1[i];
243 double my_val = Y1_vals[i];
244 int jj = i-1;
245 while(jj >= 0 && indices1[jj] > value){
246 indices1[jj+1] = indices1[jj];
247 Y1_vals[jj+1] = Y1_vals[jj];
248 jj = jj-1;
249 }
250 indices1[jj+1] = value;
251 Y1_vals[jj+1] = my_val;
252 }
253 for(int i = 1; i < entries2; ++i){
254 int value = indices2[i];
255 double my_val = Y2_vals[i];
256 int jj = i-1;
257 while(jj >= 0 && indices2[jj] > value){
258 indices2[jj+1] = indices2[jj];
259 Y2_vals[jj+1] = Y2_vals[jj];
260 jj = jj-1;
261 }
262 indices2[jj+1] = value;
263 Y2_vals[jj+1] = my_val;
264 }
265
266 for(int i = 0; i < entries1; i++){
267 if(indices1[i] != indices2[i]){
268 printf("indices[%d], %d != %d\n", i, indices1[i], indices2[i]);
269 retVal = false;
270 }
271 if(fabs(Y1_vals[i] - Y2_vals[i]) > tol){
272 printf("Failed at %d\n", i);
273 retVal = false;
274 }
275 }
276 }
277 Teuchos::Array<int> vals; vals.resize(HypreMatrix.Comm().NumProc());
278 int my_vals[1]; my_vals[0] = (int)retVal;
279 HypreMatrix.Comm().GatherAll(my_vals, &vals[0], 1);
280 for(int i = 0; i < HypreMatrix.Comm().NumProc(); i++){
281 if(vals[i] == false){
282 retVal = false;
283 }
284 }
285 if(retVal == false){
286 printf("[%d]Failed matrix equivalency test.\n", MyPID);
287 }
288 return retVal;
289}
Copy
const double tol
const int N
int NumMyRowEntries(int MyRow, int &NumEntries) const
int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
virtual const Epetra_Map & RowMatrixColMap() const
virtual int NumMyRows() const
virtual int NumGlobalRows() const
virtual const Epetra_Map & RowMatrixRowMap() const
int NumProc() const
Epetra_CrsMatrix * GetCrsMatrix(EpetraExt_HypreIJMatrix *Matrix)
bool EquivalentVectors(Epetra_MultiVector &Y1, Epetra_MultiVector &Y2, const double tol)
EpetraExt_HypreIJMatrix * newHypreMatrix(const int N)
bool EquivalentMatrices(Epetra_RowMatrix &HypreMatrix, Epetra_RowMatrix &CrsMatrix, const double tol)
Epetra_CrsMatrix * newCrsMatrix(int N)