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