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