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