57 HYPRE_IJMatrix Matrix;
61 MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
62 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
64 int ilower = (
N/numprocs)*rank;
65 int iupper = (
N/numprocs)*(rank+1);
66 if(rank == numprocs-1){
71 ierr += HYPRE_IJMatrixCreate(MPI_COMM_WORLD, ilower, iupper-1, ilower, iupper-1, &Matrix);
72 ierr += HYPRE_IJMatrixSetObjectType(Matrix,HYPRE_PARCSR);
74 ierr += HYPRE_IJMatrixInitialize(Matrix);
77 printf(
"Error! Matrix is NULL!\n");
87 for(i = ilower; i < iupper; i++) {
91 ierr += HYPRE_IJMatrixSetValues(Matrix, 1, &ncols, rows, cols, values);
100 for(i = ilower; i < iupper; i++) {
103 values[0] = ( (double)rand()/(double)RAND_MAX ) * 100;
104 ierr += HYPRE_IJMatrixSetValues(Matrix, 1, &ncols, rows, cols, values);
107 }
else if(type == 2){
113 for(i = ilower; i < iupper; i++) {
114 for(
int j = 0; j <
N; j++){
119 ierr += HYPRE_IJMatrixSetValues(Matrix, 1, &ncols, rows, &cols[0], &values[0]);
121 }
else if(type == 3){
127 for(i = ilower; i < iupper; i++) {
128 for(
int j = 0; j <
N; j++){
130 double currVal = ( (double)rand()/(double)RAND_MAX ) * 100;
134 ierr += HYPRE_IJMatrixSetValues(Matrix, 1, &ncols, rows, &cols[0], &values[0]);
139 for(i = ilower; i < iupper; i++) {
140 int ncols = (int)(1+( (
double)rand()/(
double)RAND_MAX ) * 0.5*(
N-1));
144 for(
int j = 0; j < ncols; j++){
146 if(i-(ncols/2) >= 0 && i+(ncols/2) <
N){
147 index = j + i - (ncols/2);
148 }
else if (i-(ncols/2) < 0){
151 index = j + i - (ncols-1);
155 double currVal = ( (double)rand()/(double)RAND_MAX ) * 100;
159 ierr += HYPRE_IJMatrixSetValues(Matrix, 1, &ncols, rows, &cols[0], &values[0]);
163 ierr += HYPRE_IJMatrixAssemble(Matrix);
171 int N = Matrix.NumGlobalRows();
172 Epetra_CrsMatrix* TestMat =
new Epetra_CrsMatrix(
Copy, Matrix.RowMatrixRowMap(), Matrix.RowMatrixColMap(),
N,
false);
174 for(
int i = 0; i < Matrix.NumMyRows(); i++){
181 for(
int j = 0; j < NumEntries; j++){
183 currVal[0] = Values[j];
185 currInd[0] = Matrix.RowMatrixColMap().GID(Indices[j]);
186 TestMat->InsertGlobalValues(Matrix.RowMatrixRowMap().GID(i), 1, currVal, currInd);
189 TestMat->FillComplete();
197 int num_vectors = Y1.NumVectors();
198 if(Y2.NumVectors() != num_vectors){
199 printf(
"Multivectors do not have same number of vectors.\n");
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");
210 (*Y1(j)).ExtractCopy(&Y1_vals[0]);
211 (*Y2(j)).ExtractCopy(&Y2_vals[0]);
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]);
222 printf(
"Failed equivalent vectors.\n");
231 for(
int j = 0; j < HypreMatrix.NumMyRows(); j++){
236 HypreMatrix.NumMyRowEntries(j,NumEntries);
243 HypreMatrix.ExtractMyRowCopy(j,NumEntries, entries1, &Y1_vals[0], &indices1[0]);
244 CrsMatrix.ExtractMyRowCopy(j,NumEntries, entries2, &Y2_vals[0], &indices2[0]);
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];
252 retVal = retVal && (Y1map == Y2map);
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){
263 printf(
"[%d]Failed matrix equivalency test.\n", HypreMatrix.Comm().MyPID());
264 if(HypreMatrix.Comm().MyPID() == 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.
int NumMyRowEntries(int MyRow, int &NumEntries) const
Return the current number of values stored for the specified local row.
#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)
Epetra_CrsMatrix * newCrsMatrix(EpetraExt_HypreIJMatrix &Matrix)
void resize(size_type new_size, const value_type &x=value_type())
EpetraExt_HypreIJMatrix * newHypreMatrix(const int N, const int type)