EpetraExt Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
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 
55 EpetraExt_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);
164  EpetraExt_HypreIJMatrix* RetMat = new EpetraExt_HypreIJMatrix(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 
229 bool 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 }
virtual const Epetra_Map & RowMatrixRowMap() const
virtual int GatherAll(double *MyVals, double *AllVals, int Count) const =0
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
bool EquivalentVectors(Epetra_MultiVector &Y1, Epetra_MultiVector &Y2, const double tol)
bool EquivalentMatrices(Epetra_RowMatrix &HypreMatrix, Epetra_RowMatrix &CrsMatrix, const double tol)
virtual int MyPID() const =0
int FillComplete(bool OptimizeDataStorage=true)
int NumMyRowEntries(int MyRow, int &NumEntries) const
Return the current number of values stored for the specified local row.
virtual int NumMyRowEntries(int MyRow, int &NumEntries) const =0
const int N
virtual const Epetra_Comm & Comm() const =0
Epetra_CrsMatrix * newCrsMatrix(EpetraExt_HypreIJMatrix &Matrix)
int GID(int LID) const
const double tol
virtual int NumGlobalRows() const
virtual int NumMyRows() const =0
virtual const Epetra_Map & RowMatrixColMap() const
void resize(size_type new_size, const value_type &x=value_type())
virtual int NumMyRows() const
virtual int NumProc() const =0
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const =0
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.
EpetraExt_HypreIJMatrix * newHypreMatrix(const int N, const int type)