43 #include "Teuchos_UnitTestHarness.hpp" 44 #include "Teuchos_Assert.hpp" 45 #include "Teuchos_Array.hpp" 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){
110 Teuchos::Array<int> cols; cols.resize(
N);
111 Teuchos::Array<double> values; values.resize(
N);
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){
124 Teuchos::Array<int> cols; cols.resize(
N);
125 Teuchos::Array<double> values; values.resize(
N);
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));
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++){
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);
174 for(
int i = 0; i < Matrix.
NumMyRows(); i++){
177 Teuchos::Array<double> Values; Values.resize(entries);
178 Teuchos::Array<int> Indices; Indices.resize(entries);
181 for(
int j = 0; j < NumEntries; j++){
183 currVal[0] = Values[j];
199 printf(
"Multivectors do not have same number of vectors.\n");
203 for(
int j = 0; j < num_vectors; j++){
205 printf(
"Vectors are not same size on local processor.\n");
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]);
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++){
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);
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);
254 Teuchos::Array<int> vals; vals.resize(HypreMatrix.
Comm().
NumProc());
255 int my_vals[1]; my_vals[0] = (int)retVal;
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());
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.
virtual int GatherAll(double *MyVals, double *AllVals, int Count) const=0
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
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)
virtual const Epetra_Map & RowMatrixRowMap() const
virtual int NumMyRowEntries(int MyRow, int &NumEntries) const=0
virtual const Epetra_Comm & Comm() const=0
Epetra_CrsMatrix * newCrsMatrix(EpetraExt_HypreIJMatrix &Matrix)
virtual int NumMyRows() const=0
virtual int NumProc() const=0
virtual const Epetra_Map & RowMatrixColMap() const
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
EpetraExt_HypreIJMatrix * newHypreMatrix(const int N, const int type)
virtual int NumMyRows() const
virtual int NumGlobalRows() const