43 #include "Teuchos_UnitTestHarness.hpp" 44 #include "Teuchos_Assert.hpp" 45 #include "Teuchos_Array.hpp" 46 #include "Teuchos_RCP.hpp" 49 #include "Galeri_Maps.h" 50 #include "Galeri_CrsMatrices.h" 51 #include "Galeri_Utils.h" 65 HYPRE_IJMatrix Matrix;
69 MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
70 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
72 int ilower = (
N/numprocs)*rank;
73 int iupper = (
N/numprocs)*(rank+1);
74 if(rank == numprocs-1){
79 ierr += HYPRE_IJMatrixCreate(MPI_COMM_WORLD, ilower, iupper-1, ilower, iupper-1, &Matrix);
80 ierr += HYPRE_IJMatrixSetObjectType(Matrix,HYPRE_PARCSR);
82 ierr += HYPRE_IJMatrixInitialize(Matrix);
85 printf(
"Error! Matrix is NULL!\n");
91 for(i = ilower; i < iupper; i++) {
92 int ncols = (int)(1+( (
double)rand()/(
double)RAND_MAX ) * 0.5*(
N-1));
96 for(
int j = 0; j < ncols; j++){
98 if(i-(ncols/2) >= 0 && i+(ncols/2) <
N){
99 index = j + i - (ncols/2);
100 }
else if (i-(ncols/2) < 0){
103 index = j + i - (ncols-1);
107 double currVal = ( (double)rand()/(double)RAND_MAX ) * 100;
111 ierr += HYPRE_IJMatrixAddToValues(Matrix, 1, &ncols, rows, &cols[0], &values[0]);
114 ierr += HYPRE_IJMatrixAssemble(Matrix);
122 Epetra_MpiComm Comm(MPI_COMM_WORLD);
124 Epetra_CrsMatrix* Matrix;
128 int nx =
N * Comm.NumProc();
129 int ny =
N * Comm.NumProc();
130 GaleriList.
set(
"nx", nx);
131 GaleriList.
set(
"ny", ny);
135 Epetra_Map Map(nx*ny,0,Comm);
137 Matrix = Galeri::CreateCrsMatrix(
"Laplace2D", &Map, GaleriList);
143 int N = Matrix->NumGlobalRows();
144 Epetra_CrsMatrix* TestMat =
new Epetra_CrsMatrix(
Copy, Matrix->RowMatrixRowMap(), Matrix->RowMatrixColMap(),
N,
false);
146 for(
int i = 0; i < Matrix->NumMyRows(); i++){
153 for(
int j = 0; j < NumEntries; j++){
155 currVal[0] = Values[j];
157 currInd[0] = Matrix->RowMatrixColMap().GID(Indices[j]);
158 TestMat->InsertGlobalValues(Matrix->RowMatrixRowMap().GID(i), 1, currVal, currInd);
161 TestMat->FillComplete();
169 int num_vectors = Y1.NumVectors();
170 if(Y2.NumVectors() != num_vectors){
171 printf(
"Multivectors do not have same number of vectors.\n");
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");
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){
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]);
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){
202 printf(
"[%d]Failed vector equivalency test.\n", Y1.Comm().MyPID());
211 int MyPID = HypreMatrix.Comm().MyPID();
212 if(HypreMatrix.NumMyRows() != CrsMatrix.NumMyRows()){
213 printf(
"Different number of local rows.");
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){
223 int NumEntries = HypreMatrix.NumMyCols();
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);
237 for(
int i = 0; i < entries1; i++){
238 indices1[i] = HypreMatrix.RowMatrixColMap().GID(indices1[i]);
239 indices2[i] = CrsMatrix.RowMatrixColMap().GID(indices2[i]);
241 for(
int i = 1; i < entries1; ++i){
242 int value = indices1[i];
243 double my_val = Y1_vals[i];
245 while(jj >= 0 && indices1[jj] > value){
246 indices1[jj+1] = indices1[jj];
247 Y1_vals[jj+1] = Y1_vals[jj];
250 indices1[jj+1] = value;
251 Y1_vals[jj+1] = my_val;
253 for(
int i = 1; i < entries2; ++i){
254 int value = indices2[i];
255 double my_val = Y2_vals[i];
257 while(jj >= 0 && indices2[jj] > value){
258 indices2[jj+1] = indices2[jj];
259 Y2_vals[jj+1] = Y2_vals[jj];
262 indices2[jj+1] = value;
263 Y2_vals[jj+1] = my_val;
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]);
271 if(fabs(Y1_vals[i] - Y2_vals[i]) >
tol){
272 printf(
"Failed at %d\n", i);
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){
286 printf(
"[%d]Failed matrix equivalency test.\n", MyPID);
int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
int NumMyRowEntries(int MyRow, int &NumEntries) const
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Epetra_CrsMatrix * GetCrsMatrix(EpetraExt_HypreIJMatrix *Matrix)
bool EquivalentVectors(Epetra_MultiVector &Y1, Epetra_MultiVector &Y2, const double tol)
bool EquivalentMatrices(Epetra_RowMatrix &HypreMatrix, Epetra_RowMatrix &CrsMatrix, const double tol)
EpetraExt_HypreIJMatrix * newHypreMatrix(const int N)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void resize(size_type new_size, const value_type &x=value_type())
Epetra_CrsMatrix * newCrsMatrix(int N)