xref: /petsc/src/mat/tests/ex70.c (revision d5b43468fb8780a8feea140ccd6fa3e6a50411cc)
1 #include <petscmat.h>
2 
3 static char help[] = "Tests MatMat operations with MAT_REUSE_MATRIX and already allocated dense result.\n\n";
4 
5 static PetscScalar MAGIC_NUMBER = 12345;
6 
7 static PetscErrorCode CheckLocal(Mat A, Mat B, PetscScalar *a, PetscScalar *b)
8 {
9   PetscBool wA = PETSC_FALSE, wB = PETSC_FALSE;
10   PetscBool wAv = PETSC_FALSE, wBv = PETSC_FALSE;
11   PetscInt  lda, i, j, m, n;
12 
13   PetscFunctionBegin;
14   if (a) {
15     const PetscScalar *Aa;
16     PetscCall(MatDenseGetArrayRead(A, &Aa));
17     wA = (PetscBool)(a != Aa);
18     PetscCall(MatDenseGetLDA(A, &lda));
19     PetscCall(MatGetLocalSize(A, &m, &n));
20     for (j = 0; j < n; j++) {
21       for (i = m; i < lda; i++) {
22         if (Aa[j * lda + i] != MAGIC_NUMBER) wAv = PETSC_TRUE;
23       }
24     }
25     PetscCall(MatDenseRestoreArrayRead(A, &Aa));
26   }
27   if (b) {
28     const PetscScalar *Bb;
29     PetscCall(MatDenseGetArrayRead(B, &Bb));
30     wB = (PetscBool)(b != Bb);
31     PetscCall(MatDenseGetLDA(B, &lda));
32     PetscCall(MatGetLocalSize(B, &m, &n));
33     for (j = 0; j < n; j++) {
34       for (i = m; i < lda; i++) {
35         if (Bb[j * lda + i] != MAGIC_NUMBER) wBv = PETSC_TRUE;
36       }
37     }
38     PetscCall(MatDenseRestoreArrayRead(B, &Bb));
39   }
40   PetscCheck(!wA && !wB, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Wrong array in first Mat? %d, Wrong array in second Mat? %d", wA, wB);
41   PetscCheck(!wAv && !wBv, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Wrong data in first Mat? %d, Wrong data in second Mat? %d", wAv, wBv);
42   PetscFunctionReturn(0);
43 }
44 
45 typedef struct {
46   Mat A;
47   Mat P;
48   Mat R;
49 } proj_data;
50 
51 PetscErrorCode proj_destroy(void *ctx)
52 {
53   proj_data *userdata = (proj_data *)ctx;
54 
55   PetscFunctionBegin;
56   PetscCheck(userdata, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing userdata");
57   PetscCall(MatDestroy(&userdata->A));
58   PetscCall(MatDestroy(&userdata->P));
59   PetscCall(MatDestroy(&userdata->R));
60   PetscCall(PetscFree(userdata));
61   PetscFunctionReturn(0);
62 }
63 
64 PetscErrorCode proj_mult(Mat S, Vec X, Vec Y)
65 {
66   Mat        A, R, P;
67   Vec        Ax, Ay;
68   Vec        Px, Py;
69   proj_data *userdata;
70 
71   PetscFunctionBegin;
72   PetscCall(MatShellGetContext(S, &userdata));
73   PetscCheck(userdata, PetscObjectComm((PetscObject)S), PETSC_ERR_PLIB, "Missing userdata");
74   A = userdata->A;
75   R = userdata->R;
76   P = userdata->P;
77   PetscCheck(A, PetscObjectComm((PetscObject)S), PETSC_ERR_PLIB, "Missing matrix");
78   PetscCheck(R || P, PetscObjectComm((PetscObject)S), PETSC_ERR_PLIB, "Missing projectors");
79   PetscCheck(!R || !P, PetscObjectComm((PetscObject)S), PETSC_ERR_PLIB, "Both projectors");
80   PetscCall(MatCreateVecs(A, &Ax, &Ay));
81   if (R) {
82     PetscCall(MatCreateVecs(R, &Py, &Px));
83   } else {
84     PetscCall(MatCreateVecs(P, &Px, &Py));
85   }
86   PetscCall(VecCopy(X, Px));
87   if (P) {
88     PetscCall(MatMult(P, Px, Py));
89   } else {
90     PetscCall(MatMultTranspose(R, Px, Py));
91   }
92   PetscCall(VecCopy(Py, Ax));
93   PetscCall(MatMult(A, Ax, Ay));
94   PetscCall(VecCopy(Ay, Py));
95   if (P) {
96     PetscCall(MatMultTranspose(P, Py, Px));
97   } else {
98     PetscCall(MatMult(R, Py, Px));
99   }
100   PetscCall(VecCopy(Px, Y));
101   PetscCall(VecDestroy(&Px));
102   PetscCall(VecDestroy(&Py));
103   PetscCall(VecDestroy(&Ax));
104   PetscCall(VecDestroy(&Ay));
105   PetscFunctionReturn(0);
106 }
107 
108 PetscErrorCode MyPtShellPMultSymbolic(Mat S, Mat P, Mat PtAP, void **ctx)
109 {
110   proj_data *userdata;
111 
112   PetscFunctionBegin;
113   PetscCall(PetscNew(&userdata));
114   PetscCall(MatShellSetContext(PtAP, userdata));
115   *ctx = (void *)userdata;
116   PetscFunctionReturn(0);
117 }
118 
119 PetscErrorCode MyPtShellPMultNumeric(Mat S, Mat P, Mat PtAP, void *ctx)
120 {
121   Mat        A;
122   proj_data *userdata = (proj_data *)ctx;
123 
124   PetscFunctionBegin;
125   PetscCall(MatShellGetContext(S, &A));
126   PetscCall(PetscObjectReference((PetscObject)A));
127   PetscCall(PetscObjectReference((PetscObject)P));
128   PetscCall(MatDestroy(&userdata->A));
129   PetscCall(MatDestroy(&userdata->P));
130   PetscCall(MatDestroy(&userdata->R));
131   userdata->A = A;
132   userdata->P = P;
133   PetscCall(MatShellSetOperation(PtAP, MATOP_MULT, (void (*)(void))proj_mult));
134   PetscCall(MatSetUp(PtAP));
135   PetscCall(MatAssemblyBegin(PtAP, MAT_FINAL_ASSEMBLY));
136   PetscCall(MatAssemblyEnd(PtAP, MAT_FINAL_ASSEMBLY));
137   PetscFunctionReturn(0);
138 }
139 
140 PetscErrorCode MyRShellRtMultSymbolic(Mat S, Mat R, Mat RARt, void **ctx)
141 {
142   proj_data *userdata;
143 
144   PetscFunctionBegin;
145   PetscCall(PetscNew(&userdata));
146   PetscCall(MatShellSetContext(RARt, userdata));
147   *ctx = (void *)userdata;
148   PetscFunctionReturn(0);
149 }
150 
151 PetscErrorCode MyRShellRtMultNumeric(Mat S, Mat R, Mat RARt, void *ctx)
152 {
153   Mat        A;
154   proj_data *userdata = (proj_data *)ctx;
155 
156   PetscFunctionBegin;
157   PetscCall(MatShellGetContext(S, &A));
158   PetscCall(PetscObjectReference((PetscObject)A));
159   PetscCall(PetscObjectReference((PetscObject)R));
160   PetscCall(MatDestroy(&userdata->A));
161   PetscCall(MatDestroy(&userdata->P));
162   PetscCall(MatDestroy(&userdata->R));
163   userdata->A = A;
164   userdata->R = R;
165   PetscCall(MatShellSetOperation(RARt, MATOP_MULT, (void (*)(void))proj_mult));
166   PetscCall(MatSetUp(RARt));
167   PetscCall(MatAssemblyBegin(RARt, MAT_FINAL_ASSEMBLY));
168   PetscCall(MatAssemblyEnd(RARt, MAT_FINAL_ASSEMBLY));
169   PetscFunctionReturn(0);
170 }
171 
172 PetscErrorCode MyMatShellMatMultNumeric(Mat S, Mat B, Mat C, void *ctx)
173 {
174   Mat A;
175 
176   PetscFunctionBegin;
177   PetscCall(MatShellGetContext(S, &A));
178   PetscCall(MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C));
179   PetscFunctionReturn(0);
180 }
181 
182 PetscErrorCode MyMatTransposeShellMatMultNumeric(Mat S, Mat B, Mat C, void *ctx)
183 {
184   Mat A;
185 
186   PetscFunctionBegin;
187   PetscCall(MatShellGetContext(S, &A));
188   PetscCall(MatTransposeMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C));
189   PetscFunctionReturn(0);
190 }
191 
192 PetscErrorCode MyMatShellMatTransposeMultNumeric(Mat S, Mat B, Mat C, void *ctx)
193 {
194   Mat A;
195 
196   PetscFunctionBegin;
197   PetscCall(MatShellGetContext(S, &A));
198   PetscCall(MatMatTransposeMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C));
199   PetscFunctionReturn(0);
200 }
201 
202 int main(int argc, char **args)
203 {
204   Mat          X, B, A, Bt, T, T2, PtAP = NULL, RARt = NULL, R = NULL;
205   Vec          r, l, rs, ls;
206   PetscInt     m, n, k, M = 10, N = 10, K = 5, ldx = 3, ldb = 5, ldr = 4;
207   const char  *deft = MATAIJ;
208   char         mattype[256];
209   PetscBool    flg, symm = PETSC_FALSE, testtt = PETSC_TRUE, testnest = PETSC_TRUE, testtranspose = PETSC_TRUE, testcircular = PETSC_FALSE, local = PETSC_TRUE;
210   PetscBool    testhtranspose = PETSC_FALSE; /* Hermitian transpose is not handled correctly and generates an error */
211   PetscBool    xgpu = PETSC_FALSE, bgpu = PETSC_FALSE, testshellops = PETSC_FALSE, testproj = PETSC_TRUE, testrart = PETSC_TRUE, testmatmatt = PETSC_TRUE, testmattmat = PETSC_TRUE;
212   PetscScalar *dataX = NULL, *dataB = NULL, *dataR = NULL, *dataBt = NULL;
213   PetscScalar *aX, *aB, *aBt;
214   PetscReal    err;
215 
216   PetscFunctionBeginUser;
217   PetscCall(PetscInitialize(&argc, &args, NULL, help));
218   PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL));
219   PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
220   PetscCall(PetscOptionsGetInt(NULL, NULL, "-K", &K, NULL));
221   PetscCall(PetscOptionsGetBool(NULL, NULL, "-symm", &symm, NULL));
222   PetscCall(PetscOptionsGetBool(NULL, NULL, "-local", &local, NULL));
223   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ldx", &ldx, NULL));
224   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ldb", &ldb, NULL));
225   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ldr", &ldr, NULL));
226   PetscCall(PetscOptionsGetBool(NULL, NULL, "-testtranspose", &testtranspose, NULL));
227   PetscCall(PetscOptionsGetBool(NULL, NULL, "-testnest", &testnest, NULL));
228   PetscCall(PetscOptionsGetBool(NULL, NULL, "-testtt", &testtt, NULL));
229   PetscCall(PetscOptionsGetBool(NULL, NULL, "-testcircular", &testcircular, NULL));
230   PetscCall(PetscOptionsGetBool(NULL, NULL, "-testshellops", &testshellops, NULL));
231   PetscCall(PetscOptionsGetBool(NULL, NULL, "-testproj", &testproj, NULL));
232   PetscCall(PetscOptionsGetBool(NULL, NULL, "-testrart", &testrart, NULL));
233   PetscCall(PetscOptionsGetBool(NULL, NULL, "-testmatmatt", &testmatmatt, NULL));
234   PetscCall(PetscOptionsGetBool(NULL, NULL, "-testmattmat", &testmattmat, NULL));
235   PetscCall(PetscOptionsGetBool(NULL, NULL, "-xgpu", &xgpu, NULL));
236   PetscCall(PetscOptionsGetBool(NULL, NULL, "-bgpu", &bgpu, NULL));
237   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-magic_number", &MAGIC_NUMBER, NULL));
238   if (M != N) testproj = PETSC_FALSE;
239 
240   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
241   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N));
242   PetscCall(MatSetType(A, MATAIJ));
243   PetscCall(MatSetUp(A));
244   PetscCall(MatSetRandom(A, NULL));
245   if (M == N && symm) {
246     Mat AT;
247 
248     PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &AT));
249     PetscCall(MatAXPY(A, 1.0, AT, DIFFERENT_NONZERO_PATTERN));
250     PetscCall(MatDestroy(&AT));
251     PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
252   }
253   PetscCall(MatViewFromOptions(A, NULL, "-A_init_view"));
254   PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "");
255   PetscCall(PetscOptionsFList("-A_mat_type", "Matrix type", "MatSetType", MatList, deft, mattype, 256, &flg));
256   PetscOptionsEnd();
257   if (flg) {
258     Mat A2;
259 
260     PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
261     PetscCall(MatConvert(A, mattype, MAT_INPLACE_MATRIX, &A));
262     PetscCall(MatMultEqual(A, A2, 10, &flg));
263     if (!flg) {
264       Mat AE, A2E;
265 
266       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with convert\n"));
267       PetscCall(MatComputeOperator(A, MATDENSE, &AE));
268       PetscCall(MatComputeOperator(A2, MATDENSE, &A2E));
269       PetscCall(MatView(AE, NULL));
270       PetscCall(MatView(A2E, NULL));
271       PetscCall(MatAXPY(A2E, -1.0, A, SAME_NONZERO_PATTERN));
272       PetscCall(MatView(A2E, NULL));
273       PetscCall(MatDestroy(&A2E));
274       PetscCall(MatDestroy(&AE));
275     }
276     PetscCall(MatDestroy(&A2));
277   }
278   PetscCall(MatViewFromOptions(A, NULL, "-A_view"));
279 
280   PetscCall(MatGetLocalSize(A, &m, &n));
281   if (local) {
282     PetscInt i;
283 
284     PetscCall(PetscMalloc1((m + ldx) * K, &dataX));
285     PetscCall(PetscMalloc1((n + ldb) * K, &dataB));
286     for (i = 0; i < (m + ldx) * K; i++) dataX[i] = MAGIC_NUMBER;
287     for (i = 0; i < (n + ldb) * K; i++) dataB[i] = MAGIC_NUMBER;
288   }
289   PetscCall(MatCreateDense(PETSC_COMM_WORLD, n, PETSC_DECIDE, N, K, dataB, &B));
290   PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, M, K, dataX, &X));
291   if (local) {
292     PetscCall(MatDenseSetLDA(X, m + ldx));
293     PetscCall(MatDenseSetLDA(B, n + ldb));
294   }
295   PetscCall(MatGetLocalSize(B, NULL, &k));
296   if (local) {
297     PetscInt i;
298 
299     PetscCall(PetscMalloc1((k + ldr) * N, &dataBt));
300     for (i = 0; i < (k + ldr) * N; i++) dataBt[i] = MAGIC_NUMBER;
301   }
302   PetscCall(MatCreateDense(PETSC_COMM_WORLD, k, n, K, N, dataBt, &Bt));
303   if (local) PetscCall(MatDenseSetLDA(Bt, k + ldr));
304 
305   /* store pointer to dense data for testing */
306   PetscCall(MatDenseGetArrayRead(B, (const PetscScalar **)&dataB));
307   PetscCall(MatDenseGetArrayRead(X, (const PetscScalar **)&dataX));
308   PetscCall(MatDenseGetArrayRead(Bt, (const PetscScalar **)&dataBt));
309   aX  = dataX;
310   aB  = dataB;
311   aBt = dataBt;
312   PetscCall(MatDenseRestoreArrayRead(Bt, (const PetscScalar **)&dataBt));
313   PetscCall(MatDenseRestoreArrayRead(B, (const PetscScalar **)&dataB));
314   PetscCall(MatDenseRestoreArrayRead(X, (const PetscScalar **)&dataX));
315   if (local) {
316     dataX  = aX;
317     dataB  = aB;
318     dataBt = aBt;
319   }
320 
321   PetscCall(MatSetRandom(X, NULL));
322   PetscCall(MatSetRandom(B, NULL));
323   PetscCall(MatSetRandom(Bt, NULL));
324   PetscCall(CheckLocal(X, NULL, aX, NULL));
325   PetscCall(CheckLocal(Bt, B, aBt, aB));
326 
327   /* convert to CUDA if needed */
328   if (bgpu) {
329     PetscCall(MatConvert(B, MATDENSECUDA, MAT_INPLACE_MATRIX, &B));
330     PetscCall(MatConvert(Bt, MATDENSECUDA, MAT_INPLACE_MATRIX, &Bt));
331   }
332   if (xgpu) PetscCall(MatConvert(X, MATDENSECUDA, MAT_INPLACE_MATRIX, &X));
333   PetscCall(CheckLocal(B, X, aB, aX));
334 
335   /* Test MatDenseGetSubMatrix */
336   {
337     Mat B2, T3, T4;
338 
339     PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
340     PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &T4));
341     PetscCall(MatSetRandom(T4, NULL));
342     PetscCall(MatAXPY(B2, 1.0, T4, SAME_NONZERO_PATTERN));
343     PetscCall(MatDenseGetSubMatrix(B, PETSC_DECIDE, PETSC_DECIDE, PetscMin(1, K - 1), PetscMin(2, K), &T));
344     PetscCall(MatDenseGetSubMatrix(T4, PETSC_DECIDE, PETSC_DECIDE, PetscMin(1, K - 1), PetscMin(2, K), &T2));
345     PetscCall(MatDenseGetSubMatrix(B2, PETSC_DECIDE, PETSC_DECIDE, PetscMin(1, K - 1), PetscMin(2, K), &T3));
346     PetscCall(MatAXPY(T, 1.0, T2, SAME_NONZERO_PATTERN));
347     PetscCall(MatAXPY(T3, -1.0, T, SAME_NONZERO_PATTERN));
348     PetscCall(MatNorm(T3, NORM_FROBENIUS, &err));
349     if (err > PETSC_SMALL) {
350       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with MatDenseGetSubMatrix\n"));
351       PetscCall(MatView(T3, NULL));
352     }
353     PetscCall(MatDenseRestoreSubMatrix(B, &T));
354     PetscCall(MatDenseRestoreSubMatrix(T4, &T2));
355     PetscCall(MatDenseRestoreSubMatrix(B2, &T3));
356     PetscCall(CheckLocal(B, NULL, aB, NULL));
357     PetscCall(MatDestroy(&B2));
358     PetscCall(MatDestroy(&T4));
359     if (N >= 2) {
360       PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
361       PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &T4));
362       PetscCall(MatSetRandom(T4, NULL));
363       PetscCall(MatAXPY(B2, 1.0, T4, SAME_NONZERO_PATTERN));
364       PetscCall(MatDenseGetSubMatrix(B, N - 2, PETSC_DECIDE, PetscMin(1, K - 1), PetscMin(2, K), &T));
365       PetscCall(MatDenseGetSubMatrix(T4, N - 2, PETSC_DECIDE, PetscMin(1, K - 1), PetscMin(2, K), &T2));
366       PetscCall(MatDenseGetSubMatrix(B2, N - 2, PETSC_DECIDE, PetscMin(1, K - 1), PetscMin(2, K), &T3));
367       PetscCall(MatAXPY(T, 1.0, T2, SAME_NONZERO_PATTERN));
368       PetscCall(MatAXPY(T3, -1.0, T, SAME_NONZERO_PATTERN));
369       PetscCall(MatNorm(T3, NORM_FROBENIUS, &err));
370       if (err > PETSC_SMALL) {
371         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with MatDenseGetSubMatrix\n"));
372         PetscCall(MatView(T3, NULL));
373       }
374       PetscCall(MatDenseRestoreSubMatrix(B, &T));
375       PetscCall(MatDenseRestoreSubMatrix(T4, &T2));
376       PetscCall(MatDenseRestoreSubMatrix(B2, &T3));
377       PetscCall(CheckLocal(B, NULL, aB, NULL));
378       PetscCall(MatDestroy(&B2));
379       PetscCall(MatDestroy(&T4));
380       PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
381       PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &T4));
382       PetscCall(MatSetRandom(T4, NULL));
383       PetscCall(MatAXPY(B2, 1.0, T4, SAME_NONZERO_PATTERN));
384       PetscCall(MatDenseGetSubMatrix(B, PETSC_DECIDE, 2, PetscMin(1, K - 1), PetscMin(2, K), &T));
385       PetscCall(MatDenseGetSubMatrix(T4, PETSC_DECIDE, 2, PetscMin(1, K - 1), PetscMin(2, K), &T2));
386       PetscCall(MatDenseGetSubMatrix(B2, PETSC_DECIDE, 2, PetscMin(1, K - 1), PetscMin(2, K), &T3));
387       PetscCall(MatAXPY(T, 1.0, T2, SAME_NONZERO_PATTERN));
388       PetscCall(MatAXPY(T3, -1.0, T, SAME_NONZERO_PATTERN));
389       PetscCall(MatNorm(T3, NORM_FROBENIUS, &err));
390       if (err > PETSC_SMALL) {
391         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with MatDenseGetSubMatrix\n"));
392         PetscCall(MatView(T3, NULL));
393       }
394       PetscCall(MatDenseRestoreSubMatrix(B, &T));
395       PetscCall(MatDenseRestoreSubMatrix(T4, &T2));
396       PetscCall(MatDenseRestoreSubMatrix(B2, &T3));
397       PetscCall(CheckLocal(B, NULL, aB, NULL));
398       PetscCall(MatDestroy(&B2));
399       PetscCall(MatDestroy(&T4));
400     }
401   }
402 
403   /* Test reusing a previously allocated dense buffer */
404   PetscCall(MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &X));
405   PetscCall(CheckLocal(B, X, aB, aX));
406   PetscCall(MatMatMultEqual(A, B, X, 10, &flg));
407   if (!flg) {
408     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with reusage\n"));
409     PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T));
410     PetscCall(MatAXPY(T, -1.0, X, SAME_NONZERO_PATTERN));
411     PetscCall(MatView(T, NULL));
412     PetscCall(MatDestroy(&T));
413   }
414 
415   /* Test MatTransposeMat and MatMatTranspose */
416   if (testmattmat) {
417     PetscCall(MatTransposeMatMult(A, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &B));
418     PetscCall(CheckLocal(B, X, aB, aX));
419     PetscCall(MatTransposeMatMultEqual(A, X, B, 10, &flg));
420     if (!flg) {
421       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with reusage (MatTransposeMat)\n"));
422       PetscCall(MatTransposeMatMult(A, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &B));
423       PetscCall(MatAXPY(T, -1.0, B, SAME_NONZERO_PATTERN));
424       PetscCall(MatView(T, NULL));
425       PetscCall(MatDestroy(&T));
426     }
427   }
428   if (testmatmatt) {
429     PetscCall(MatMatTransposeMult(A, Bt, MAT_REUSE_MATRIX, PETSC_DEFAULT, &X));
430     PetscCall(CheckLocal(Bt, X, aBt, aX));
431     PetscCall(MatMatTransposeMultEqual(A, Bt, X, 10, &flg));
432     if (!flg) {
433       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with reusage (MatMatTranspose)\n"));
434       PetscCall(MatMatTransposeMult(A, Bt, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T));
435       PetscCall(MatAXPY(T, -1.0, X, SAME_NONZERO_PATTERN));
436       PetscCall(MatView(T, NULL));
437       PetscCall(MatDestroy(&T));
438     }
439   }
440 
441   /* Test projection operations (PtAP and RARt) */
442   if (testproj) {
443     PetscCall(MatPtAP(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &PtAP));
444     PetscCall(MatPtAPMultEqual(A, B, PtAP, 10, &flg));
445     if (!flg) {
446       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with PtAP\n"));
447       PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T));
448       PetscCall(MatTransposeMatMult(B, T, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T2));
449       PetscCall(MatAXPY(T2, -1.0, PtAP, SAME_NONZERO_PATTERN));
450       PetscCall(MatView(T2, NULL));
451       PetscCall(MatDestroy(&T2));
452       PetscCall(MatDestroy(&T));
453     }
454     PetscCall(PetscMalloc1((k + ldr) * M, &dataR));
455     PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, m, K, M, dataR, &R));
456     PetscCall(MatDenseSetLDA(R, k + ldr));
457     PetscCall(MatSetRandom(R, NULL));
458     if (testrart) { /* fails for AIJCUSPARSE because RA operation is not defined */
459       PetscCall(MatRARt(A, R, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &RARt));
460       PetscCall(MatRARtMultEqual(A, R, RARt, 10, &flg));
461       if (!flg) {
462         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with RARt\n"));
463         PetscCall(MatMatTransposeMult(A, R, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T));
464         PetscCall(MatMatMult(R, T, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T2));
465         PetscCall(MatAXPY(T2, -1.0, RARt, SAME_NONZERO_PATTERN));
466         PetscCall(MatView(T2, NULL));
467         PetscCall(MatDestroy(&T2));
468         PetscCall(MatDestroy(&T));
469       }
470     }
471   }
472 
473   /* Test MatDenseGetColumnVec and friends */
474   PetscCall(MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &X));
475   PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T));
476   PetscCall(MatDuplicate(T, MAT_DO_NOT_COPY_VALUES, &T2));
477   for (k = 0; k < K; k++) {
478     Vec Xv, Tv, T2v;
479 
480     PetscCall(MatDenseGetColumnVecRead(X, k, &Xv));
481     PetscCall(MatDenseGetColumnVec(T, k, &Tv));
482     PetscCall(MatDenseGetColumnVecWrite(T2, k, &T2v));
483     PetscCall(VecCopy(Xv, T2v));
484     PetscCall(VecAXPY(Tv, -1., Xv));
485     PetscCall(MatDenseRestoreColumnVecRead(X, k, &Xv));
486     PetscCall(MatDenseRestoreColumnVec(T, k, &Tv));
487     PetscCall(MatDenseRestoreColumnVecWrite(T2, k, &T2v));
488   }
489   PetscCall(MatNorm(T, NORM_FROBENIUS, &err));
490   if (err > PETSC_SMALL) {
491     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with MatDenseGetColumnVec\n"));
492     PetscCall(MatView(T, NULL));
493   }
494   PetscCall(MatAXPY(T2, -1., X, SAME_NONZERO_PATTERN));
495   PetscCall(MatNorm(T2, NORM_FROBENIUS, &err));
496   if (err > PETSC_SMALL) {
497     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with MatDenseGetColumnVecWrite\n"));
498     PetscCall(MatView(T2, NULL));
499   }
500   PetscCall(MatDestroy(&T));
501   PetscCall(MatDestroy(&T2));
502 
503   /* Test with MatShell */
504   PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &T));
505   PetscCall(MatConvert(T, MATSHELL, MAT_INITIAL_MATRIX, &T2));
506   PetscCall(MatDestroy(&T));
507 
508   /* scale matrix */
509   PetscCall(MatScale(A, 2.0));
510   PetscCall(MatScale(T2, 2.0));
511   PetscCall(MatCreateVecs(A, &r, &l));
512   PetscCall(VecSetRandom(r, NULL));
513   PetscCall(VecSetRandom(l, NULL));
514   PetscCall(MatCreateVecs(T2, &rs, &ls));
515   PetscCall(VecCopy(r, rs));
516   PetscCall(VecCopy(l, ls));
517   if (testproj) {
518     PetscCall(MatDiagonalScale(A, r, r));
519     PetscCall(MatDiagonalScale(T2, rs, rs));
520   } else {
521     PetscCall(MatDiagonalScale(A, l, r));
522     PetscCall(MatDiagonalScale(T2, ls, rs));
523   }
524   PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &T));
525   PetscCall(MatAXPY(A, 4.5, T, SAME_NONZERO_PATTERN));
526   PetscCall(MatAXPY(T2, 4.5, T, DIFFERENT_NONZERO_PATTERN));
527   PetscCall(MatMultEqual(T2, A, 10, &flg));
528   if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with MATSHELL (MatMult)\n"));
529   PetscCall(MatMultTransposeEqual(T2, A, 10, &flg));
530   if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with MATSHELL (MatMultTranspose)\n"));
531   PetscCall(MatDestroy(&T));
532   PetscCall(VecDestroy(&ls));
533   PetscCall(VecDestroy(&rs));
534   PetscCall(VecDestroy(&l));
535   PetscCall(VecDestroy(&r));
536 
537   /* recompute projections, test reusage */
538   if (PtAP) PetscCall(MatPtAP(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &PtAP));
539   if (RARt) PetscCall(MatRARt(A, R, MAT_REUSE_MATRIX, PETSC_DEFAULT, &RARt));
540   if (testshellops) { /* test callbacks for user defined MatProducts */
541     PetscCall(MatShellSetMatProductOperation(T2, MATPRODUCT_AB, NULL, MyMatShellMatMultNumeric, NULL, MATDENSE, MATDENSE));
542     PetscCall(MatShellSetMatProductOperation(T2, MATPRODUCT_AB, NULL, MyMatShellMatMultNumeric, NULL, MATDENSECUDA, MATDENSECUDA));
543     PetscCall(MatShellSetMatProductOperation(T2, MATPRODUCT_AtB, NULL, MyMatTransposeShellMatMultNumeric, NULL, MATDENSE, MATDENSE));
544     PetscCall(MatShellSetMatProductOperation(T2, MATPRODUCT_AtB, NULL, MyMatTransposeShellMatMultNumeric, NULL, MATDENSECUDA, MATDENSECUDA));
545     PetscCall(MatShellSetMatProductOperation(T2, MATPRODUCT_ABt, NULL, MyMatShellMatTransposeMultNumeric, NULL, MATDENSE, MATDENSE));
546     PetscCall(MatShellSetMatProductOperation(T2, MATPRODUCT_ABt, NULL, MyMatShellMatTransposeMultNumeric, NULL, MATDENSECUDA, MATDENSECUDA));
547     if (testproj) {
548       PetscCall(MatShellSetMatProductOperation(T2, MATPRODUCT_PtAP, MyPtShellPMultSymbolic, MyPtShellPMultNumeric, proj_destroy, MATDENSE, MATSHELL));
549       PetscCall(MatShellSetMatProductOperation(T2, MATPRODUCT_PtAP, MyPtShellPMultSymbolic, MyPtShellPMultNumeric, proj_destroy, MATDENSECUDA, MATSHELL));
550       PetscCall(MatShellSetMatProductOperation(T2, MATPRODUCT_RARt, MyRShellRtMultSymbolic, MyRShellRtMultNumeric, proj_destroy, MATDENSE, MATSHELL));
551       PetscCall(MatShellSetMatProductOperation(T2, MATPRODUCT_RARt, MyRShellRtMultSymbolic, MyRShellRtMultNumeric, proj_destroy, MATDENSECUDA, MATSHELL));
552     }
553   }
554   PetscCall(CheckLocal(B, X, aB, aX));
555   /* we either use the shell operations or the loop over columns code, applying the operator */
556   PetscCall(MatMatMult(T2, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &X));
557   PetscCall(CheckLocal(B, X, aB, aX));
558   PetscCall(MatMatMultEqual(T2, B, X, 10, &flg));
559   if (!flg) {
560     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with reusage (MATSHELL)\n"));
561     PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T));
562     PetscCall(MatAXPY(T, -1.0, X, SAME_NONZERO_PATTERN));
563     PetscCall(MatView(T, NULL));
564     PetscCall(MatDestroy(&T));
565   }
566   if (testproj) {
567     PetscCall(MatPtAPMultEqual(T2, B, PtAP, 10, &flg));
568     if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with PtAP (MATSHELL)\n"));
569     if (testshellops) { /* projections fail if the product operations are not specified */
570       PetscCall(MatPtAP(T2, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T));
571       PetscCall(MatPtAP(T2, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &T));
572       PetscCall(MatPtAPMultEqual(T2, B, T, 10, &flg));
573       if (!flg) {
574         Mat TE;
575 
576         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with PtAP (MATSHELL user defined)\n"));
577         PetscCall(MatComputeOperator(T, MATDENSE, &TE));
578         PetscCall(MatView(TE, NULL));
579         PetscCall(MatView(PtAP, NULL));
580         PetscCall(MatAXPY(TE, -1.0, PtAP, SAME_NONZERO_PATTERN));
581         PetscCall(MatView(TE, NULL));
582         PetscCall(MatDestroy(&TE));
583       }
584       PetscCall(MatDestroy(&T));
585     }
586     if (RARt) {
587       PetscCall(MatRARtMultEqual(T2, R, RARt, 10, &flg));
588       if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with RARt (MATSHELL)\n"));
589     }
590     if (testshellops) {
591       PetscCall(MatRARt(T2, R, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T));
592       PetscCall(MatRARt(T2, R, MAT_REUSE_MATRIX, PETSC_DEFAULT, &T));
593       PetscCall(MatRARtMultEqual(T2, R, T, 10, &flg));
594       if (!flg) {
595         Mat TE;
596 
597         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with RARt (MATSHELL user defined)\n"));
598         PetscCall(MatComputeOperator(T, MATDENSE, &TE));
599         PetscCall(MatView(TE, NULL));
600         if (RARt) {
601           PetscCall(MatView(RARt, NULL));
602           PetscCall(MatAXPY(TE, -1.0, RARt, SAME_NONZERO_PATTERN));
603           PetscCall(MatView(TE, NULL));
604         }
605         PetscCall(MatDestroy(&TE));
606       }
607       PetscCall(MatDestroy(&T));
608     }
609   }
610 
611   if (testmattmat) { /* we either use the shell operations or the loop over columns code applying the transposed operator */
612     PetscCall(MatTransposeMatMult(T2, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &B));
613     PetscCall(CheckLocal(B, X, aB, aX));
614     PetscCall(MatTransposeMatMultEqual(T2, X, B, 10, &flg));
615     if (!flg) {
616       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with reusage (MatTranspose, MATSHELL)\n"));
617       PetscCall(MatTransposeMatMult(A, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T));
618       PetscCall(MatAXPY(T, -1.0, B, SAME_NONZERO_PATTERN));
619       PetscCall(MatView(T, NULL));
620       PetscCall(MatDestroy(&T));
621     }
622   }
623   if (testmatmatt && testshellops) { /* only when shell operations are set */
624     PetscCall(MatMatTransposeMult(T2, Bt, MAT_REUSE_MATRIX, PETSC_DEFAULT, &X));
625     PetscCall(CheckLocal(Bt, X, aBt, aX));
626     PetscCall(MatMatTransposeMultEqual(T2, Bt, X, 10, &flg));
627     if (!flg) {
628       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with reusage (MatMatTranspose, MATSHELL)\n"));
629       PetscCall(MatMatTransposeMult(A, Bt, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T));
630       PetscCall(MatAXPY(T, -1.0, X, SAME_NONZERO_PATTERN));
631       PetscCall(MatView(T, NULL));
632       PetscCall(MatDestroy(&T));
633     }
634   }
635   PetscCall(MatDestroy(&T2));
636 
637   if (testnest) { /* test with MatNest */
638     Mat NA;
639 
640     PetscCall(MatCreateNest(PETSC_COMM_WORLD, 1, NULL, 1, NULL, &A, &NA));
641     PetscCall(MatViewFromOptions(NA, NULL, "-NA_view"));
642     PetscCall(MatMatMult(NA, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &X));
643     PetscCall(CheckLocal(B, X, aB, aX));
644     PetscCall(MatMatMultEqual(NA, B, X, 10, &flg));
645     if (!flg) {
646       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with Nest\n"));
647       PetscCall(MatMatMult(NA, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T));
648       PetscCall(MatAXPY(T, -1.0, X, SAME_NONZERO_PATTERN));
649       PetscCall(MatView(T, NULL));
650       PetscCall(MatDestroy(&T));
651     }
652     PetscCall(MatDestroy(&NA));
653   }
654 
655   if (testtranspose) { /* test with Transpose */
656     Mat TA;
657 
658     PetscCall(MatCreateTranspose(A, &TA));
659     PetscCall(MatMatMult(TA, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &B));
660     PetscCall(CheckLocal(B, X, aB, aX));
661     PetscCall(MatMatMultEqual(TA, X, B, 10, &flg));
662     if (!flg) {
663       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with Transpose\n"));
664       PetscCall(MatMatMult(TA, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T));
665       PetscCall(MatAXPY(T, -1.0, B, SAME_NONZERO_PATTERN));
666       PetscCall(MatView(T, NULL));
667       PetscCall(MatDestroy(&T));
668     }
669     PetscCall(MatDestroy(&TA));
670   }
671 
672   if (testhtranspose) { /* test with Hermitian Transpose */
673     Mat TA;
674 
675     PetscCall(MatCreateHermitianTranspose(A, &TA));
676     PetscCall(MatMatMult(TA, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &B));
677     PetscCall(CheckLocal(B, X, aB, aX));
678     PetscCall(MatMatMultEqual(TA, X, B, 10, &flg));
679     if (!flg) {
680       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with Transpose\n"));
681       PetscCall(MatMatMult(TA, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T));
682       PetscCall(MatAXPY(T, -1.0, B, SAME_NONZERO_PATTERN));
683       PetscCall(MatView(T, NULL));
684       PetscCall(MatDestroy(&T));
685     }
686     PetscCall(MatDestroy(&TA));
687   }
688 
689   if (testtt) { /* test with Transpose(Transpose) */
690     Mat TA, TTA;
691 
692     PetscCall(MatCreateTranspose(A, &TA));
693     PetscCall(MatCreateTranspose(TA, &TTA));
694     PetscCall(MatMatMult(TTA, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &X));
695     PetscCall(CheckLocal(B, X, aB, aX));
696     PetscCall(MatMatMultEqual(TTA, B, X, 10, &flg));
697     if (!flg) {
698       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with Transpose(Transpose)\n"));
699       PetscCall(MatMatMult(TTA, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T));
700       PetscCall(MatAXPY(T, -1.0, X, SAME_NONZERO_PATTERN));
701       PetscCall(MatView(T, NULL));
702       PetscCall(MatDestroy(&T));
703     }
704     PetscCall(MatDestroy(&TA));
705     PetscCall(MatDestroy(&TTA));
706   }
707 
708   if (testcircular) { /* test circular */
709     Mat AB;
710 
711     PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AB));
712     PetscCall(MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &X));
713     PetscCall(CheckLocal(B, X, aB, aX));
714     if (M == N && N == K) {
715       PetscCall(MatMatMult(A, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &B));
716     } else {
717       PetscCall(MatTransposeMatMult(A, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &B));
718     }
719     PetscCall(CheckLocal(B, X, aB, aX));
720     PetscCall(MatDestroy(&AB));
721   }
722 
723   /* Test by Pierre Jolivet */
724   {
725     Mat C, D, D2, AtA;
726     PetscCall(MatCreateNormal(A, &AtA));
727     PetscCall(MatDuplicate(X, MAT_DO_NOT_COPY_VALUES, &C));
728     PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &D));
729     PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &D2));
730     PetscCall(MatSetRandom(B, NULL));
731     PetscCall(MatSetRandom(C, NULL));
732     PetscCall(MatSetRandom(D, NULL));
733     PetscCall(MatSetRandom(D2, NULL));
734     PetscCall(MatProductCreateWithMat(A, B, NULL, C));
735     PetscCall(MatProductSetType(C, MATPRODUCT_AB));
736     PetscCall(MatProductSetFromOptions(C));
737     PetscCall(MatProductSymbolic(C));
738     PetscCall(MatProductCreateWithMat(A, C, NULL, D));
739     PetscCall(MatProductSetType(D, MATPRODUCT_AtB));
740     PetscCall(MatProductSetFromOptions(D));
741     PetscCall(MatProductSymbolic(D));
742     PetscCall(MatProductNumeric(C));
743     PetscCall(MatMatMultEqual(A, B, C, 10, &flg));
744     if (!flg) {
745       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with Normal (AB != C)\n"));
746       PetscCall(MatView(A, NULL));
747       PetscCall(MatView(B, NULL));
748       PetscCall(MatView(C, NULL));
749     }
750     PetscCall(MatProductNumeric(D));
751     PetscCall(MatMatMultEqual(AtA, B, D, 10, &flg));
752     if (!flg) {
753       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with Normal (2)\n"));
754       PetscCall(MatMatMult(AtA, C, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T));
755       PetscCall(MatView(D, NULL));
756       PetscCall(MatView(T, NULL));
757       PetscCall(MatAXPY(T, -1.0, D, SAME_NONZERO_PATTERN));
758       PetscCall(MatView(T, NULL));
759       PetscCall(MatDestroy(&T));
760     }
761     PetscCall(MatDestroy(&C));
762     PetscCall(MatDestroy(&D));
763     PetscCall(MatDestroy(&D2));
764     PetscCall(MatDestroy(&AtA));
765   }
766 
767   PetscCall(MatDestroy(&X));
768   PetscCall(MatDestroy(&Bt));
769   PetscCall(MatDestroy(&B));
770   PetscCall(MatDestroy(&A));
771   PetscCall(MatDestroy(&R));
772   PetscCall(MatDestroy(&PtAP));
773   PetscCall(MatDestroy(&RARt));
774   PetscCall(PetscFree(dataX));
775   PetscCall(PetscFree(dataB));
776   PetscCall(PetscFree(dataR));
777   PetscCall(PetscFree(dataBt));
778   PetscCall(PetscFinalize());
779   return 0;
780 }
781 
782 /*TEST
783 
784   test:
785     suffix: 1
786     args: -local {{0 1}} -testshellops
787 
788   test:
789     output_file: output/ex70_1.out
790     requires: cuda
791     suffix: 1_cuda
792     args: -local {{0 1}} -xgpu {{0 1}} -bgpu {{0 1}} -A_mat_type {{seqaijcusparse seqaij}} -testshellops {{0 1}}
793 
794   test:
795     output_file: output/ex70_1.out
796     nsize: 2
797     suffix: 1_par
798     args: -local {{0 1}} -testmatmatt 0
799 
800   test:
801     output_file: output/ex70_1.out
802     requires: cuda
803     nsize: 2
804     suffix: 1_par_cuda
805     args: -local {{0 1}} -xgpu {{0 1}} -bgpu {{0 1}} -A_mat_type {{mpiaijcusparse mpiaij}} -testnest 0 -testmatmatt 0 -matmatmult_Bbn 3
806 
807   test:
808     output_file: output/ex70_1.out
809     suffix: 2
810     nsize: 1
811     args: -M {{7 11}} -N {{12 9}} -K {{1 3}} -local {{0 1}}
812 
813   testset:
814     requires: cuda
815     output_file: output/ex70_1.out
816     nsize: 1
817     args: -M 7 -N 9 -K 2 -local {{0 1}} -testnest 0 -A_mat_type {{seqdensecuda seqdense}} -xgpu {{0 1}} -bgpu {{0 1}}
818     test:
819       requires: !complex
820       suffix: 2_cuda_real
821     test:
822       # complex+single gives a little bigger error in the MatDenseGetColumnVec test
823       requires: complex !single
824       suffix: 2_cuda_complex
825 
826   test:
827     output_file: output/ex70_1.out
828     suffix: 2_par
829     nsize: 2
830     args: -M {{7 11}} -N {{12 9}} -K {{1 3}} -local {{0 1}} -testcircular -testmatmatt 0
831 
832   test:
833     requires: cuda
834     output_file: output/ex70_1.out
835     suffix: 2_par_cuda
836     nsize: 2
837     args: -M 11 -N 9 -K 1 -local {{0 1}} -testcircular 0 -A_mat_type mpiaijcusparse -xgpu -bgpu -testnest 0 -testmatmatt 0
838 
839   test:
840     output_file: output/ex70_1.out
841     suffix: 3
842     nsize: {{1 3}}
843     args: -M 13 -N 13 -K {{1 3}} -local {{0 1}} -A_mat_type sbaij -symm -testproj 0 -testmatmatt 0
844 
845   test:
846     output_file: output/ex70_1.out
847     suffix: 4
848     nsize: 1
849     args: -M 3 -N 3 -K 3 -local {{0 1}} -testcircular
850 
851   test:
852     output_file: output/ex70_1.out
853     suffix: 5
854     nsize: {{2 4}}
855     args: -M 3 -N 3 -K 3 -local {{0 1}} -testcircular -testmatmatt 0
856 
857   test:
858     output_file: output/ex70_1.out
859     suffix: 6
860     nsize: 1
861     args: -M {{1 3}} -N {{2 5}} -K {{1 2}} -local {{0 1}} -testcircular
862 
863   test:
864     output_file: output/ex70_1.out
865     suffix: 7
866     nsize: 1
867     args: -M 13 -N 13 -K {{1 3}} -local {{0 1}} -A_mat_type dense -testnest -testcircular
868 
869 TEST*/
870