xref: /petsc/src/mat/tests/ex115.c (revision 732aec7a18f2199fb53bb9a2f3aef439a834ce31)
1 static char help[] = "Tests MatHYPRE\n";
2 
3 #include <petscmathypre.h>
4 
main(int argc,char ** args)5 int main(int argc, char **args)
6 {
7   Mat                 A, B, C, D;
8   Mat                 pAB, CD;
9   hypre_ParCSRMatrix *parcsr;
10   PetscReal           err;
11   PetscInt            i, j, N = 6, M = 6;
12   PetscBool           flg, testptap = PETSC_TRUE, testmatmatmult = PETSC_TRUE;
13   PetscReal           norm;
14   char                file[256];
15   MatType             mtype = MATAIJ;
16 
17   PetscFunctionBeginUser;
18   PetscCall(PetscInitialize(&argc, &args, NULL, help));
19   PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg));
20 #if defined(PETSC_USE_COMPLEX)
21   testptap       = PETSC_FALSE;
22   testmatmatmult = PETSC_FALSE;
23   PetscCall(PetscOptionsInsertString(NULL, "-options_left 0"));
24 #endif
25   PetscCall(PetscOptionsGetBool(NULL, NULL, "-ptap", &testptap, NULL));
26   PetscCall(PetscOptionsGetBool(NULL, NULL, "-matmatmult", &testmatmatmult, NULL));
27   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
28 #if PetscDefined(HAVE_HYPRE_DEVICE)
29   #if PetscDefined(HAVE_HIP)
30   mtype = MATAIJHIPSPARSE;
31   #elif PetscDefined(HAVE_CUDA)
32   mtype = MATAIJCUSPARSE;
33   #endif
34 #endif
35 
36   if (!flg) { /* Create a matrix and test MatSetValues */
37     PetscMPIInt size;
38 
39     PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
40     PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
41     PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL));
42     PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N));
43     PetscCall(MatSetType(A, mtype));
44     PetscCall(MatSeqAIJSetPreallocation(A, 9, NULL));
45     PetscCall(MatMPIAIJSetPreallocation(A, 9, NULL, 9, NULL));
46     PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
47     PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, M, N));
48 #if PetscDefined(HAVE_HYPRE_DEVICE)
49     PetscCall(MatSetType(B, mtype));
50 #else
51     PetscCall(MatSetType(B, MATHYPRE));
52 #endif
53     PetscCall(MatSeqAIJSetPreallocation(B, 9, NULL));
54     PetscCall(MatMPIAIJSetPreallocation(B, 9, NULL, 9, NULL));
55     if (M == N) {
56       PetscCall(MatHYPRESetPreallocation(B, 9, NULL, 9, NULL));
57     } else {
58       PetscCall(MatHYPRESetPreallocation(B, 6, NULL, 6, NULL));
59     }
60     if (M == N) {
61       for (i = 0; i < M; i++) {
62         PetscInt    cols[] = {0, 1, 2, 3, 4, 5};
63         PetscScalar vals[] = {0, 1. / size, 2. / size, 3. / size, 4. / size, 5. / size};
64         for (j = i - 2; j < i + 1; j++) {
65           if (j >= N) {
66             PetscCall(MatSetValue(A, i, N - 1, (1. * j * N + i) / (3. * N * size), ADD_VALUES));
67             PetscCall(MatSetValue(B, i, N - 1, (1. * j * N + i) / (3. * N * size), ADD_VALUES));
68           } else if (i > j) {
69             PetscCall(MatSetValue(A, i, PetscMin(j, N - 1), (1. * j * N + i) / (2. * N * size), ADD_VALUES));
70             PetscCall(MatSetValue(B, i, PetscMin(j, N - 1), (1. * j * N + i) / (2. * N * size), ADD_VALUES));
71           } else {
72             PetscCall(MatSetValue(A, i, PetscMin(j, N - 1), -1. - (1. * j * N + i) / (4. * N * size), ADD_VALUES));
73             PetscCall(MatSetValue(B, i, PetscMin(j, N - 1), -1. - (1. * j * N + i) / (4. * N * size), ADD_VALUES));
74           }
75         }
76         PetscCall(MatSetValues(A, 1, &i, PetscMin(6, N), cols, vals, ADD_VALUES));
77         PetscCall(MatSetValues(B, 1, &i, PetscMin(6, N), cols, vals, ADD_VALUES));
78       }
79     } else {
80       PetscInt  rows[2];
81       PetscBool test_offproc = PETSC_FALSE;
82 
83       PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_offproc", &test_offproc, NULL));
84       if (test_offproc) {
85         const PetscInt *ranges;
86         PetscMPIInt     rank;
87 
88         PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
89         PetscCall(MatGetOwnershipRanges(A, &ranges));
90         rows[0] = ranges[(rank + 1) % size];
91         rows[1] = ranges[(rank + 1) % size + 1];
92       } else {
93         PetscCall(MatGetOwnershipRange(A, &rows[0], &rows[1]));
94       }
95       for (i = rows[0]; i < rows[1]; i++) {
96         PetscInt    cols[] = {0, 1, 2, 3, 4, 5};
97         PetscScalar vals[] = {-1, 1, -2, 2, -3, 3};
98 
99         PetscCall(MatSetValues(A, 1, &i, PetscMin(6, N), cols, vals, INSERT_VALUES));
100         PetscCall(MatSetValues(B, 1, &i, PetscMin(6, N), cols, vals, INSERT_VALUES));
101       }
102     }
103     /* MAT_FLUSH_ASSEMBLY currently not supported */
104     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
105     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
106     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
107     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
108 #if PetscDefined(HAVE_HYPRE_DEVICE)
109     PetscCall(MatConvert(B, MATHYPRE, MAT_INPLACE_MATRIX, &B));
110 #endif
111 
112 #if defined(PETSC_USE_COMPLEX)
113     /* make the matrix imaginary */
114     PetscCall(MatScale(A, PETSC_i));
115     PetscCall(MatScale(B, PETSC_i));
116 #endif
117 
118 #if !PetscDefined(HAVE_HYPRE_DEVICE)
119     /* MatAXPY further exercises MatSetValues_HYPRE */
120     PetscCall(MatAXPY(B, -1., A, DIFFERENT_NONZERO_PATTERN));
121     PetscCall(MatConvert(B, MATMPIAIJ, MAT_INITIAL_MATRIX, &C));
122     PetscCall(MatNorm(C, NORM_INFINITY, &err));
123     PetscCheck(err <= PETSC_SMALL, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error MatSetValues %g", err);
124     PetscCall(MatDestroy(&C));
125 #endif
126     PetscCall(MatDestroy(&B));
127   } else {
128     PetscViewer viewer;
129 
130     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &viewer));
131     PetscCall(MatSetFromOptions(A));
132     PetscCall(MatLoad(A, viewer));
133     PetscCall(MatSetType(A, mtype));
134     PetscCall(PetscViewerDestroy(&viewer));
135     PetscCall(MatGetSize(A, &M, &N));
136   }
137 
138   /* check conversion routines */
139   PetscCall(MatViewFromOptions(A, NULL, "-view_A"));
140   PetscCall(MatConvert(A, MATHYPRE, MAT_INITIAL_MATRIX, &B));
141   PetscCall(MatViewFromOptions(B, NULL, "-view_convert"));
142   PetscCall(MatMultEqual(B, A, 4, &flg));
143   PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error Mat HYPRE init");
144   PetscCall(MatConvert(A, MATHYPRE, MAT_REUSE_MATRIX, &B));
145   PetscCall(MatViewFromOptions(B, NULL, "-view_convert"));
146   PetscCall(MatMultEqual(B, A, 4, &flg));
147   PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error Mat HYPRE reuse");
148   PetscCall(MatConvert(B, MATIS, MAT_INITIAL_MATRIX, &D));
149   PetscCall(MatConvert(B, MATIS, MAT_REUSE_MATRIX, &D));
150   PetscCall(MatMultEqual(D, A, 4, &flg));
151   PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error Mat IS");
152   PetscCall(MatConvert(B, MATAIJ, MAT_INITIAL_MATRIX, &C));
153   PetscCall(MatConvert(B, MATAIJ, MAT_REUSE_MATRIX, &C));
154   PetscCall(MatMultEqual(C, A, 4, &flg));
155   PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error Mat AIJ");
156   PetscCall(MatAXPY(C, -1., A, SAME_NONZERO_PATTERN));
157   PetscCall(MatNorm(C, NORM_INFINITY, &err));
158   PetscCheck(err <= PETSC_SMALL, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error Mat AIJ %g", err);
159   PetscCall(MatDestroy(&C));
160   PetscCall(MatConvert(D, MATAIJ, MAT_INITIAL_MATRIX, &C));
161   PetscCall(MatAXPY(C, -1., A, SAME_NONZERO_PATTERN));
162   PetscCall(MatNorm(C, NORM_INFINITY, &err));
163   PetscCheck(err <= PETSC_SMALL, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error Mat IS %g", err);
164   PetscCall(MatDestroy(&C));
165   PetscCall(MatDestroy(&D));
166 
167   /* check MatCreateFromParCSR */
168   PetscCall(MatHYPREGetParCSR(B, &parcsr));
169   PetscCall(MatCreateFromParCSR(parcsr, MATAIJ, PETSC_COPY_VALUES, &D));
170   PetscCall(MatDestroy(&D));
171   PetscCall(MatCreateFromParCSR(parcsr, MATHYPRE, PETSC_USE_POINTER, &C));
172 
173   /* check MatMult operations */
174   PetscCall(MatMultEqual(A, B, 4, &flg));
175   PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error MatMult B");
176   PetscCall(MatMultEqual(A, C, 4, &flg));
177   PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error MatMult C");
178   PetscCall(MatMultAddEqual(A, B, 4, &flg));
179   PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error MatMultAdd B");
180   PetscCall(MatMultAddEqual(A, C, 4, &flg));
181   PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error MatMultAdd C");
182   PetscCall(MatMultTransposeEqual(A, B, 4, &flg));
183   PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error MatMultTranspose B");
184   PetscCall(MatMultTransposeEqual(A, C, 4, &flg));
185   PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error MatMultTranspose C");
186   PetscCall(MatMultTransposeAddEqual(A, B, 4, &flg));
187   PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error MatMultTransposeAdd B");
188   PetscCall(MatMultTransposeAddEqual(A, C, 4, &flg));
189   PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error MatMultTransposeAdd C");
190 
191   /* check PtAP */
192   if (testptap && M == N) {
193     Mat pP, hP;
194 
195     /* PETSc MatPtAP -> output is a MatAIJ
196        It uses HYPRE functions when -matptap_via hypre is specified at command line */
197     PetscCall(MatPtAP(A, A, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &pP));
198     PetscCall(MatPtAP(A, A, MAT_REUSE_MATRIX, PETSC_DETERMINE, &pP));
199     PetscCall(MatNorm(pP, NORM_INFINITY, &norm));
200     PetscCall(MatPtAPMultEqual(A, A, pP, 10, &flg));
201     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatPtAP_MatAIJ");
202 
203     /* MatPtAP_HYPRE_HYPRE -> output is a MatHYPRE */
204     PetscCall(MatPtAP(C, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &hP));
205     PetscCall(MatPtAP(C, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &hP));
206     PetscCall(MatPtAPMultEqual(C, B, hP, 10, &flg));
207     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatPtAP_HYPRE_HYPRE");
208 
209     /* Test MatAXPY_Basic() */
210 #if !PetscDefined(HAVE_HYPRE_DEVICE)
211     PetscCall(MatAXPY(hP, -1., pP, DIFFERENT_NONZERO_PATTERN));
212     PetscCall(MatHasOperation(hP, MATOP_NORM, &flg));
213     if (!flg) { /* TODO add MatNorm_HYPRE */
214       PetscCall(MatConvert(hP, MATAIJ, MAT_INPLACE_MATRIX, &hP));
215     }
216     PetscCall(MatNorm(hP, NORM_INFINITY, &err));
217     PetscCheck(err / norm <= PETSC_SMALL, PetscObjectComm((PetscObject)hP), PETSC_ERR_PLIB, "Error MatPtAP %g %g", err, norm);
218 #endif
219     PetscCall(MatDestroy(&pP));
220     PetscCall(MatDestroy(&hP));
221 
222     /* MatPtAP_AIJ_HYPRE -> output can be decided at runtime with -matptap_hypre_outtype */
223 #if !PetscDefined(HAVE_HYPRE_DEVICE)
224     PetscCall(MatPtAP(A, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &hP));
225     PetscCall(MatPtAP(A, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &hP));
226     PetscCall(MatPtAPMultEqual(A, B, hP, 10, &flg));
227     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatPtAP_AIJ_HYPRE");
228     PetscCall(MatDestroy(&hP));
229 #endif
230   }
231   PetscCall(MatDestroy(&C));
232   PetscCall(MatDestroy(&B));
233 
234   /* check MatMatMult */
235   if (testmatmatmult) {
236     PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &B));
237     PetscCall(MatConvert(A, MATHYPRE, MAT_INITIAL_MATRIX, &C));
238     PetscCall(MatConvert(B, MATHYPRE, MAT_INITIAL_MATRIX, &D));
239 
240     /* PETSc MatMatMult -> output is a MatAIJ
241        It uses HYPRE functions when -matmatmult_via hypre is specified at command line */
242     PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &pAB));
243     PetscCall(MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &pAB));
244     PetscCall(MatNorm(pAB, NORM_INFINITY, &norm));
245     PetscCall(MatMatMultEqual(A, B, pAB, 10, &flg));
246     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatMatMult_AIJ_AIJ");
247 
248     /* MatMatMult_HYPRE_HYPRE -> output is a MatHYPRE */
249     PetscCall(MatMatMult(C, D, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &CD));
250     PetscCall(MatMatMult(C, D, MAT_REUSE_MATRIX, PETSC_DETERMINE, &CD));
251     PetscCall(MatMatMultEqual(C, D, CD, 10, &flg));
252     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatMatMult_HYPRE_HYPRE");
253 
254     /* Test MatAXPY_Basic() */
255 #if !PetscDefined(HAVE_HYPRE_DEVICE)
256     PetscCall(MatAXPY(CD, -1., pAB, DIFFERENT_NONZERO_PATTERN));
257 
258     PetscCall(MatHasOperation(CD, MATOP_NORM, &flg));
259     if (!flg) { /* TODO add MatNorm_HYPRE */
260       PetscCall(MatConvert(CD, MATAIJ, MAT_INPLACE_MATRIX, &CD));
261     }
262     PetscCall(MatNorm(CD, NORM_INFINITY, &err));
263     PetscCheck((err / norm) <= PETSC_SMALL, PetscObjectComm((PetscObject)CD), PETSC_ERR_PLIB, "Error MatMatMult %g %g", err, norm);
264 #endif
265 
266     PetscCall(MatDestroy(&C));
267     PetscCall(MatDestroy(&D));
268     PetscCall(MatDestroy(&pAB));
269     PetscCall(MatDestroy(&CD));
270 
271     /* When configured with HYPRE, MatMatMatMult is available for the triplet transpose(aij)-aij-aij */
272 #if !PetscDefined(HAVE_HYPRE_DEVICE)
273     Mat CAB;
274     PetscCall(MatCreateTranspose(A, &C));
275     PetscCall(MatMatMatMult(C, A, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &CAB));
276     PetscCall(MatDestroy(&C));
277     PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &C));
278     PetscCall(MatMatMult(C, A, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &D));
279     PetscCall(MatDestroy(&C));
280     PetscCall(MatMatMult(D, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C));
281     PetscCall(MatNorm(C, NORM_INFINITY, &norm));
282     PetscCall(MatAXPY(C, -1., CAB, DIFFERENT_NONZERO_PATTERN));
283     PetscCall(MatNorm(C, NORM_INFINITY, &err));
284     PetscCheck((err / norm) <= PETSC_SMALL, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error MatMatMatMult %g %g", err, norm);
285     PetscCall(MatDestroy(&C));
286     PetscCall(MatDestroy(&D));
287     PetscCall(MatDestroy(&CAB));
288 #endif
289     PetscCall(MatDestroy(&B));
290   }
291 
292   /* Check MatView */
293   PetscCall(MatConvert(A, MATHYPRE, MAT_INITIAL_MATRIX, &B));
294   PetscCall(MatViewFromOptions(B, NULL, "-view_B"));
295 
296   /* Check MatDuplicate/MatCopy */
297   for (j = 0; j < 3; j++) {
298     MatDuplicateOption dop;
299 
300     dop = MAT_COPY_VALUES;
301     if (j == 1) dop = MAT_DO_NOT_COPY_VALUES;
302     if (j == 2) dop = MAT_SHARE_NONZERO_PATTERN;
303 
304     for (i = 0; i < 3; i++) {
305       MatStructure str;
306 
307       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Dup/Copy tests: %" PetscInt_FMT " %" PetscInt_FMT "\n", j, i));
308 
309       str = DIFFERENT_NONZERO_PATTERN;
310       if (i == 1) str = SAME_NONZERO_PATTERN;
311       if (i == 2) str = SUBSET_NONZERO_PATTERN;
312 
313       PetscCall(MatDuplicate(A, dop, &C));
314       PetscCall(MatDuplicate(B, dop, &D));
315       if (dop != MAT_COPY_VALUES) {
316         PetscCall(MatCopy(A, C, str));
317         PetscCall(MatCopy(B, D, str));
318       }
319       /* AXPY with AIJ and HYPRE */
320 #if !PetscDefined(HAVE_HYPRE_DEVICE)
321       PetscCall(MatAXPY(C, -1.0, D, str));
322       PetscCall(MatNorm(C, NORM_INFINITY, &err));
323 #else
324       err = 0.0;
325 #endif
326       if (err > PETSC_SMALL) {
327         PetscCall(MatViewFromOptions(A, NULL, "-view_duplicate_diff"));
328         PetscCall(MatViewFromOptions(B, NULL, "-view_duplicate_diff"));
329         PetscCall(MatViewFromOptions(C, NULL, "-view_duplicate_diff"));
330         PetscCall(MatViewFromOptions(D, NULL, "-view_duplicate_diff"));
331         SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error test 1 MatDuplicate/MatCopy %g (%" PetscInt_FMT ",%" PetscInt_FMT ")", err, j, i);
332       }
333       /* AXPY with HYPRE and HYPRE */
334       PetscCall(MatAXPY(D, -1.0, B, str));
335       if (err > PETSC_SMALL) {
336         PetscCall(MatViewFromOptions(A, NULL, "-view_duplicate_diff"));
337         PetscCall(MatViewFromOptions(B, NULL, "-view_duplicate_diff"));
338         PetscCall(MatViewFromOptions(D, NULL, "-view_duplicate_diff"));
339         SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error test 2 MatDuplicate/MatCopy %g (%" PetscInt_FMT ",%" PetscInt_FMT ")", err, j, i);
340       }
341 #if !PetscDefined(HAVE_HYPRE_DEVICE)
342       /* Copy from HYPRE to AIJ */
343       PetscCall(MatCopy(B, C, str));
344       /* Copy from AIJ to HYPRE */
345       PetscCall(MatCopy(A, D, str));
346       /* AXPY with HYPRE and AIJ */
347       PetscCall(MatAXPY(D, -1.0, C, str));
348       PetscCall(MatHasOperation(D, MATOP_NORM, &flg));
349       if (!flg) { /* TODO add MatNorm_HYPRE */
350         PetscCall(MatConvert(D, MATAIJ, MAT_INPLACE_MATRIX, &D));
351       }
352       PetscCall(MatNorm(D, NORM_INFINITY, &err));
353       if (err > PETSC_SMALL) {
354         PetscCall(MatViewFromOptions(A, NULL, "-view_duplicate_diff"));
355         PetscCall(MatViewFromOptions(C, NULL, "-view_duplicate_diff"));
356         PetscCall(MatViewFromOptions(D, NULL, "-view_duplicate_diff"));
357         SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error test 3 MatDuplicate/MatCopy %g (%" PetscInt_FMT ",%" PetscInt_FMT ")", err, j, i);
358       }
359 #endif
360       PetscCall(MatDestroy(&C));
361       PetscCall(MatDestroy(&D));
362     }
363   }
364   PetscCall(MatDestroy(&B));
365 
366   PetscCall(MatHasCongruentLayouts(A, &flg));
367   if (flg) {
368     Vec y, y2;
369 
370     PetscCall(MatConvert(A, MATHYPRE, MAT_INITIAL_MATRIX, &B));
371     PetscCall(MatCreateVecs(A, NULL, &y));
372     PetscCall(MatCreateVecs(B, NULL, &y2));
373     PetscCall(MatGetDiagonal(A, y));
374     PetscCall(MatGetDiagonal(B, y2));
375     PetscCall(VecAXPY(y2, -1.0, y));
376     PetscCall(VecNorm(y2, NORM_INFINITY, &err));
377     if (err > PETSC_SMALL) {
378       PetscCall(VecViewFromOptions(y, NULL, "-view_diagonal_diff"));
379       PetscCall(VecViewFromOptions(y2, NULL, "-view_diagonal_diff"));
380       SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error MatGetDiagonal %g", err);
381     }
382     PetscCall(MatDestroy(&B));
383     PetscCall(VecDestroy(&y));
384     PetscCall(VecDestroy(&y2));
385   }
386 
387   PetscCall(MatDestroy(&A));
388 
389   PetscCall(PetscFinalize());
390   return 0;
391 }
392 
393 /*TEST
394 
395    build:
396       requires: hypre
397 
398    test:
399       suffix: 1
400       args: -N 11 -M 11
401       output_file: output/ex115_1.out
402 
403    test:
404       suffix: 2
405       nsize: 3
406       args: -N 13 -M 13 -matmatmult_via hypre -options_left 0
407       output_file: output/ex115_1.out
408 
409    test:
410       suffix: 3
411       nsize: 4
412       args: -M 13 -N 7 -matmatmult_via hypre -options_left 0
413       output_file: output/ex115_1.out
414 
415    test:
416       suffix: 4
417       nsize: 2
418       args: -M 12 -N 19
419       output_file: output/ex115_1.out
420 
421    test:
422       suffix: 5
423       nsize: 3
424       args: -M 13 -N 13 -options_left 0 -matptap_via hypre -matptap_hypre_outtype hypre
425       output_file: output/ex115_1.out
426 
427    test:
428       requires: !defined(PETSC_HAVE_HYPRE_DEVICE)
429       suffix: 6
430       nsize: 3
431       args: -M 12 -N 19 -test_offproc
432       output_file: output/ex115_1.out
433 
434    test:
435       requires: !defined(PETSC_HAVE_HYPRE_DEVICE)
436       suffix: 7
437       nsize: 3
438       args: -M 19 -N 12 -test_offproc -view_B ::ascii_info_detail
439       output_file: output/ex115_7.out
440 
441    test:
442       requires: !defined(PETSC_HAVE_HYPRE_DEVICE)
443       suffix: 8
444       nsize: 3
445       args: -M 1 -N 12 -test_offproc
446       output_file: output/ex115_1.out
447 
448    test:
449       requires: !defined(PETSC_HAVE_HYPRE_DEVICE)
450       suffix: 9
451       nsize: 3
452       args: -M 1 -N 2 -test_offproc
453       output_file: output/ex115_1.out
454 
455 TEST*/
456