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));
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++){
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);
146 for(
int i = 0; i < Matrix->
NumMyRows(); i++){
149 Teuchos::Array<double> Values; Values.resize(entries);
150 Teuchos::Array<int> Indices; Indices.resize(entries);
153 for(
int j = 0; j < NumEntries; j++){
155 currVal[0] = Values[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]);
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){
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();
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);
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);
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){
286 printf(
"[%d]Failed matrix equivalency test.\n", MyPID);