xref: /petsc/src/mat/tests/ex70.c (revision ebead697dbf761eb322f829370bbe90b3bd93fa3)
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   PetscFunctionBeginUser;
217   PetscCall(PetscInitialize(&argc,&args,NULL,help));
218   PetscCall(PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL));
219   PetscCall(PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL));
220   PetscCall(PetscOptionsGetInt(NULL,NULL,"-K",&K,NULL));
221   PetscCall(PetscOptionsGetBool(NULL,NULL,"-symm",&symm,NULL));
222   PetscCall(PetscOptionsGetBool(NULL,NULL,"-local",&local,NULL));
223   PetscCall(PetscOptionsGetInt(NULL,NULL,"-ldx",&ldx,NULL));
224   PetscCall(PetscOptionsGetInt(NULL,NULL,"-ldb",&ldb,NULL));
225   PetscCall(PetscOptionsGetInt(NULL,NULL,"-ldr",&ldr,NULL));
226   PetscCall(PetscOptionsGetBool(NULL,NULL,"-testtranspose",&testtranspose,NULL));
227   PetscCall(PetscOptionsGetBool(NULL,NULL,"-testnest",&testnest,NULL));
228   PetscCall(PetscOptionsGetBool(NULL,NULL,"-testtt",&testtt,NULL));
229   PetscCall(PetscOptionsGetBool(NULL,NULL,"-testcircular",&testcircular,NULL));
230   PetscCall(PetscOptionsGetBool(NULL,NULL,"-testshellops",&testshellops,NULL));
231   PetscCall(PetscOptionsGetBool(NULL,NULL,"-testproj",&testproj,NULL));
232   PetscCall(PetscOptionsGetBool(NULL,NULL,"-testrart",&testrart,NULL));
233   PetscCall(PetscOptionsGetBool(NULL,NULL,"-testmatmatt",&testmatmatt,NULL));
234   PetscCall(PetscOptionsGetBool(NULL,NULL,"-testmattmat",&testmattmat,NULL));
235   PetscCall(PetscOptionsGetBool(NULL,NULL,"-xgpu",&xgpu,NULL));
236   PetscCall(PetscOptionsGetBool(NULL,NULL,"-bgpu",&bgpu,NULL));
237   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-magic_number",&MAGIC_NUMBER,NULL));
238   if (M != N) testproj = PETSC_FALSE;
239 
240   PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
241   PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N));
242   PetscCall(MatSetType(A,MATAIJ));
243   PetscCall(MatSetUp(A));
244   PetscCall(MatSetRandom(A,NULL));
245   if (M==N && symm) {
246     Mat AT;
247 
248     PetscCall(MatTranspose(A,MAT_INITIAL_MATRIX,&AT));
249     PetscCall(MatAXPY(A,1.0,AT,DIFFERENT_NONZERO_PATTERN));
250     PetscCall(MatDestroy(&AT));
251     PetscCall(MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE));
252   }
253   PetscCall(MatViewFromOptions(A,NULL,"-A_init_view"));
254   PetscOptionsBegin(PETSC_COMM_WORLD,"","","");
255   PetscCall(PetscOptionsFList("-A_mat_type","Matrix type","MatSetType",MatList,deft,mattype,256,&flg));
256   PetscOptionsEnd();
257   if (flg) {
258     Mat A2;
259 
260     PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&A2));
261     PetscCall(MatConvert(A,mattype,MAT_INPLACE_MATRIX,&A));
262     PetscCall(MatMultEqual(A,A2,10,&flg));
263     if (!flg) {
264       Mat AE,A2E;
265 
266       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error with convert\n"));
267       PetscCall(MatComputeOperator(A,MATDENSE,&AE));
268       PetscCall(MatComputeOperator(A2,MATDENSE,&A2E));
269       PetscCall(MatView(AE,NULL));
270       PetscCall(MatView(A2E,NULL));
271       PetscCall(MatAXPY(A2E,-1.0,A,SAME_NONZERO_PATTERN));
272       PetscCall(MatView(A2E,NULL));
273       PetscCall(MatDestroy(&A2E));
274       PetscCall(MatDestroy(&AE));
275     }
276     PetscCall(MatDestroy(&A2));
277   }
278   PetscCall(MatViewFromOptions(A,NULL,"-A_view"));
279 
280   PetscCall(MatGetLocalSize(A,&m,&n));
281   if (local) {
282     PetscInt i;
283 
284     PetscCall(PetscMalloc1((m+ldx)*K,&dataX));
285     PetscCall(PetscMalloc1((n+ldb)*K,&dataB));
286     for (i=0;i<(m+ldx)*K;i++) dataX[i] = MAGIC_NUMBER;
287     for (i=0;i<(n+ldb)*K;i++) dataB[i] = MAGIC_NUMBER;
288   }
289   PetscCall(MatCreateDense(PETSC_COMM_WORLD,n,PETSC_DECIDE,N,K,dataB,&B));
290   PetscCall(MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,M,K,dataX,&X));
291   if (local) {
292     PetscCall(MatDenseSetLDA(X,m+ldx));
293     PetscCall(MatDenseSetLDA(B,n+ldb));
294   }
295   PetscCall(MatGetLocalSize(B,NULL,&k));
296   if (local) {
297     PetscInt i;
298 
299     PetscCall(PetscMalloc1((k+ldr)*N,&dataBt));
300     for (i=0;i<(k+ldr)*N;i++) dataBt[i] = MAGIC_NUMBER;
301   }
302   PetscCall(MatCreateDense(PETSC_COMM_WORLD,k,n,K,N,dataBt,&Bt));
303   if (local) PetscCall(MatDenseSetLDA(Bt,k+ldr));
304 
305   /* store pointer to dense data for testing */
306   PetscCall(MatDenseGetArrayRead(B,(const PetscScalar**)&dataB));
307   PetscCall(MatDenseGetArrayRead(X,(const PetscScalar**)&dataX));
308   PetscCall(MatDenseGetArrayRead(Bt,(const PetscScalar**)&dataBt));
309   aX   = dataX;
310   aB   = dataB;
311   aBt  = dataBt;
312   PetscCall(MatDenseRestoreArrayRead(Bt,(const PetscScalar**)&dataBt));
313   PetscCall(MatDenseRestoreArrayRead(B,(const PetscScalar**)&dataB));
314   PetscCall(MatDenseRestoreArrayRead(X,(const PetscScalar**)&dataX));
315   if (local) {
316     dataX  = aX;
317     dataB  = aB;
318     dataBt = aBt;
319   }
320 
321   PetscCall(MatSetRandom(X,NULL));
322   PetscCall(MatSetRandom(B,NULL));
323   PetscCall(MatSetRandom(Bt,NULL));
324   PetscCall(CheckLocal(X,NULL,aX,NULL));
325   PetscCall(CheckLocal(Bt,B,aBt,aB));
326 
327   /* convert to CUDA if needed */
328   if (bgpu) {
329     PetscCall(MatConvert(B,MATDENSECUDA,MAT_INPLACE_MATRIX,&B));
330     PetscCall(MatConvert(Bt,MATDENSECUDA,MAT_INPLACE_MATRIX,&Bt));
331   }
332   if (xgpu) {
333     PetscCall(MatConvert(X,MATDENSECUDA,MAT_INPLACE_MATRIX,&X));
334   }
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_DEFAULT,&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_DEFAULT,&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_DEFAULT,&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_DEFAULT,&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_DEFAULT,&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_DEFAULT,&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_DEFAULT,&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_DEFAULT,&T));
450       PetscCall(MatTransposeMatMult(B,T,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&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_DEFAULT,&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_DEFAULT,&T));
466         PetscCall(MatMatMult(R,T,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&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_DEFAULT,&X));
477   PetscCall(MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&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) {
531     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error with MATSHELL (MatMult)\n"));
532   }
533   PetscCall(MatMultTransposeEqual(T2,A,10,&flg));
534   if (!flg) {
535     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error with MATSHELL (MatMultTranspose)\n"));
536   }
537   PetscCall(MatDestroy(&T));
538   PetscCall(VecDestroy(&ls));
539   PetscCall(VecDestroy(&rs));
540   PetscCall(VecDestroy(&l));
541   PetscCall(VecDestroy(&r));
542 
543   /* recompute projections, test reusage */
544   if (PtAP) PetscCall(MatPtAP(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&PtAP));
545   if (RARt) PetscCall(MatRARt(A,R,MAT_REUSE_MATRIX,PETSC_DEFAULT,&RARt));
546   if (testshellops) { /* test callbacks for user defined MatProducts */
547     PetscCall(MatShellSetMatProductOperation(T2,MATPRODUCT_AB,NULL,MyMatShellMatMultNumeric,NULL,MATDENSE,MATDENSE));
548     PetscCall(MatShellSetMatProductOperation(T2,MATPRODUCT_AB,NULL,MyMatShellMatMultNumeric,NULL,MATDENSECUDA,MATDENSECUDA));
549     PetscCall(MatShellSetMatProductOperation(T2,MATPRODUCT_AtB,NULL,MyMatTransposeShellMatMultNumeric,NULL,MATDENSE,MATDENSE));
550     PetscCall(MatShellSetMatProductOperation(T2,MATPRODUCT_AtB,NULL,MyMatTransposeShellMatMultNumeric,NULL,MATDENSECUDA,MATDENSECUDA));
551     PetscCall(MatShellSetMatProductOperation(T2,MATPRODUCT_ABt,NULL,MyMatShellMatTransposeMultNumeric,NULL,MATDENSE,MATDENSE));
552     PetscCall(MatShellSetMatProductOperation(T2,MATPRODUCT_ABt,NULL,MyMatShellMatTransposeMultNumeric,NULL,MATDENSECUDA,MATDENSECUDA));
553     if (testproj) {
554       PetscCall(MatShellSetMatProductOperation(T2,MATPRODUCT_PtAP,MyPtShellPMultSymbolic,MyPtShellPMultNumeric,proj_destroy,MATDENSE,MATSHELL));
555       PetscCall(MatShellSetMatProductOperation(T2,MATPRODUCT_PtAP,MyPtShellPMultSymbolic,MyPtShellPMultNumeric,proj_destroy,MATDENSECUDA,MATSHELL));
556       PetscCall(MatShellSetMatProductOperation(T2,MATPRODUCT_RARt,MyRShellRtMultSymbolic,MyRShellRtMultNumeric,proj_destroy,MATDENSE,MATSHELL));
557       PetscCall(MatShellSetMatProductOperation(T2,MATPRODUCT_RARt,MyRShellRtMultSymbolic,MyRShellRtMultNumeric,proj_destroy,MATDENSECUDA,MATSHELL));
558     }
559   }
560   PetscCall(CheckLocal(B,X,aB,aX));
561   /* we either use the shell operations or the loop over columns code, applying the operator */
562   PetscCall(MatMatMult(T2,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X));
563   PetscCall(CheckLocal(B,X,aB,aX));
564   PetscCall(MatMatMultEqual(T2,B,X,10,&flg));
565   if (!flg) {
566     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error with reusage (MATSHELL)\n"));
567     PetscCall(MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T));
568     PetscCall(MatAXPY(T,-1.0,X,SAME_NONZERO_PATTERN));
569     PetscCall(MatView(T,NULL));
570     PetscCall(MatDestroy(&T));
571   }
572   if (testproj) {
573     PetscCall(MatPtAPMultEqual(T2,B,PtAP,10,&flg));
574     if (!flg) {
575       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error with PtAP (MATSHELL)\n"));
576     }
577     if (testshellops) { /* projections fail if the product operations are not specified */
578       PetscCall(MatPtAP(T2,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T));
579       PetscCall(MatPtAP(T2,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&T));
580       PetscCall(MatPtAPMultEqual(T2,B,T,10,&flg));
581       if (!flg) {
582         Mat TE;
583 
584         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error with PtAP (MATSHELL user defined)\n"));
585         PetscCall(MatComputeOperator(T,MATDENSE,&TE));
586         PetscCall(MatView(TE,NULL));
587         PetscCall(MatView(PtAP,NULL));
588         PetscCall(MatAXPY(TE,-1.0,PtAP,SAME_NONZERO_PATTERN));
589         PetscCall(MatView(TE,NULL));
590         PetscCall(MatDestroy(&TE));
591       }
592       PetscCall(MatDestroy(&T));
593     }
594     if (RARt) {
595       PetscCall(MatRARtMultEqual(T2,R,RARt,10,&flg));
596       if (!flg) {
597         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error with RARt (MATSHELL)\n"));
598       }
599     }
600     if (testshellops) {
601       PetscCall(MatRARt(T2,R,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T));
602       PetscCall(MatRARt(T2,R,MAT_REUSE_MATRIX,PETSC_DEFAULT,&T));
603       PetscCall(MatRARtMultEqual(T2,R,T,10,&flg));
604       if (!flg) {
605         Mat TE;
606 
607         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error with RARt (MATSHELL user defined)\n"));
608         PetscCall(MatComputeOperator(T,MATDENSE,&TE));
609         PetscCall(MatView(TE,NULL));
610         if (RARt) {
611           PetscCall(MatView(RARt,NULL));
612           PetscCall(MatAXPY(TE,-1.0,RARt,SAME_NONZERO_PATTERN));
613           PetscCall(MatView(TE,NULL));
614         }
615         PetscCall(MatDestroy(&TE));
616       }
617       PetscCall(MatDestroy(&T));
618     }
619   }
620 
621   if (testmattmat) { /* we either use the shell operations or the loop over columns code applying the transposed operator */
622     PetscCall(MatTransposeMatMult(T2,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B));
623     PetscCall(CheckLocal(B,X,aB,aX));
624     PetscCall(MatTransposeMatMultEqual(T2,X,B,10,&flg));
625     if (!flg) {
626       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error with reusage (MatTranspose, MATSHELL)\n"));
627       PetscCall(MatTransposeMatMult(A,X,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T));
628       PetscCall(MatAXPY(T,-1.0,B,SAME_NONZERO_PATTERN));
629       PetscCall(MatView(T,NULL));
630       PetscCall(MatDestroy(&T));
631     }
632   }
633   if (testmatmatt && testshellops) { /* only when shell operations are set */
634     PetscCall(MatMatTransposeMult(T2,Bt,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X));
635     PetscCall(CheckLocal(Bt,X,aBt,aX));
636     PetscCall(MatMatTransposeMultEqual(T2,Bt,X,10,&flg));
637     if (!flg) {
638       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error with reusage (MatMatTranspose, MATSHELL)\n"));
639       PetscCall(MatMatTransposeMult(A,Bt,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T));
640       PetscCall(MatAXPY(T,-1.0,X,SAME_NONZERO_PATTERN));
641       PetscCall(MatView(T,NULL));
642       PetscCall(MatDestroy(&T));
643     }
644   }
645   PetscCall(MatDestroy(&T2));
646 
647   if (testnest) { /* test with MatNest */
648     Mat NA;
649 
650     PetscCall(MatCreateNest(PETSC_COMM_WORLD,1,NULL,1,NULL,&A,&NA));
651     PetscCall(MatViewFromOptions(NA,NULL,"-NA_view"));
652     PetscCall(MatMatMult(NA,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X));
653     PetscCall(CheckLocal(B,X,aB,aX));
654     PetscCall(MatMatMultEqual(NA,B,X,10,&flg));
655     if (!flg) {
656       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error with Nest\n"));
657       PetscCall(MatMatMult(NA,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T));
658       PetscCall(MatAXPY(T,-1.0,X,SAME_NONZERO_PATTERN));
659       PetscCall(MatView(T,NULL));
660       PetscCall(MatDestroy(&T));
661     }
662     PetscCall(MatDestroy(&NA));
663   }
664 
665   if (testtranspose) { /* test with Transpose */
666     Mat TA;
667 
668     PetscCall(MatCreateTranspose(A,&TA));
669     PetscCall(MatMatMult(TA,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B));
670     PetscCall(CheckLocal(B,X,aB,aX));
671     PetscCall(MatMatMultEqual(TA,X,B,10,&flg));
672     if (!flg) {
673       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error with Transpose\n"));
674       PetscCall(MatMatMult(TA,X,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T));
675       PetscCall(MatAXPY(T,-1.0,B,SAME_NONZERO_PATTERN));
676       PetscCall(MatView(T,NULL));
677       PetscCall(MatDestroy(&T));
678     }
679     PetscCall(MatDestroy(&TA));
680   }
681 
682   if (testhtranspose) { /* test with Hermitian Transpose */
683     Mat TA;
684 
685     PetscCall(MatCreateHermitianTranspose(A,&TA));
686     PetscCall(MatMatMult(TA,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B));
687     PetscCall(CheckLocal(B,X,aB,aX));
688     PetscCall(MatMatMultEqual(TA,X,B,10,&flg));
689     if (!flg) {
690       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error with Transpose\n"));
691       PetscCall(MatMatMult(TA,X,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T));
692       PetscCall(MatAXPY(T,-1.0,B,SAME_NONZERO_PATTERN));
693       PetscCall(MatView(T,NULL));
694       PetscCall(MatDestroy(&T));
695     }
696     PetscCall(MatDestroy(&TA));
697   }
698 
699   if (testtt) { /* test with Transpose(Transpose) */
700     Mat TA, TTA;
701 
702     PetscCall(MatCreateTranspose(A,&TA));
703     PetscCall(MatCreateTranspose(TA,&TTA));
704     PetscCall(MatMatMult(TTA,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X));
705     PetscCall(CheckLocal(B,X,aB,aX));
706     PetscCall(MatMatMultEqual(TTA,B,X,10,&flg));
707     if (!flg) {
708       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error with Transpose(Transpose)\n"));
709       PetscCall(MatMatMult(TTA,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T));
710       PetscCall(MatAXPY(T,-1.0,X,SAME_NONZERO_PATTERN));
711       PetscCall(MatView(T,NULL));
712       PetscCall(MatDestroy(&T));
713     }
714     PetscCall(MatDestroy(&TA));
715     PetscCall(MatDestroy(&TTA));
716   }
717 
718   if (testcircular) { /* test circular */
719     Mat AB;
720 
721     PetscCall(MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AB));
722     PetscCall(MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X));
723     PetscCall(CheckLocal(B,X,aB,aX));
724     if (M == N && N == K) {
725       PetscCall(MatMatMult(A,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B));
726     } else {
727       PetscCall(MatTransposeMatMult(A,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B));
728     }
729     PetscCall(CheckLocal(B,X,aB,aX));
730     PetscCall(MatDestroy(&AB));
731   }
732 
733   /* Test by Pierre Jolivet */
734   {
735     Mat C,D,D2,AtA;
736     PetscCall(MatCreateNormal(A,&AtA));
737     PetscCall(MatDuplicate(X,MAT_DO_NOT_COPY_VALUES,&C));
738     PetscCall(MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&D));
739     PetscCall(MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&D2));
740     PetscCall(MatSetRandom(B,NULL));
741     PetscCall(MatSetRandom(C,NULL));
742     PetscCall(MatSetRandom(D,NULL));
743     PetscCall(MatSetRandom(D2,NULL));
744     PetscCall(MatProductCreateWithMat(A,B,NULL,C));
745     PetscCall(MatProductSetType(C,MATPRODUCT_AB));
746     PetscCall(MatProductSetFromOptions(C));
747     PetscCall(MatProductSymbolic(C));
748     PetscCall(MatProductCreateWithMat(A,C,NULL,D));
749     PetscCall(MatProductSetType(D, MATPRODUCT_AtB));
750     PetscCall(MatProductSetFromOptions(D));
751     PetscCall(MatProductSymbolic(D));
752     PetscCall(MatProductNumeric(C));
753     PetscCall(MatProductNumeric(D));
754     PetscCall(MatMatMultEqual(AtA,B,D,10,&flg));
755     if (!flg) {
756       PetscCall(MatMatMult(AtA,C,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T));
757       PetscCall(MatAXPY(T,-1.0,D,SAME_NONZERO_PATTERN));
758       PetscCall(MatView(T,NULL));
759       PetscCall(MatDestroy(&T));
760     }
761     PetscCall(MatDestroy(&C));
762     PetscCall(MatDestroy(&D));
763     PetscCall(MatDestroy(&D2));
764     PetscCall(MatDestroy(&AtA));
765   }
766 
767   PetscCall(MatDestroy(&X));
768   PetscCall(MatDestroy(&Bt));
769   PetscCall(MatDestroy(&B));
770   PetscCall(MatDestroy(&A));
771   PetscCall(MatDestroy(&R));
772   PetscCall(MatDestroy(&PtAP));
773   PetscCall(MatDestroy(&RARt));
774   PetscCall(PetscFree(dataX));
775   PetscCall(PetscFree(dataB));
776   PetscCall(PetscFree(dataR));
777   PetscCall(PetscFree(dataBt));
778   PetscCall(PetscFinalize());
779   return 0;
780 }
781 
782 /*TEST
783 
784   test:
785     suffix: 1
786     args: -local {{0 1}} -testshellops
787 
788   test:
789     output_file: output/ex70_1.out
790     requires: cuda
791     suffix: 1_cuda
792     args: -local {{0 1}} -xgpu {{0 1}} -bgpu {{0 1}} -A_mat_type {{seqaijcusparse seqaij}} -testshellops {{0 1}}
793 
794   test:
795     output_file: output/ex70_1.out
796     nsize: 2
797     suffix: 1_par
798     args: -local {{0 1}} -testmatmatt 0
799 
800   test:
801     output_file: output/ex70_1.out
802     requires: cuda
803     nsize: 2
804     suffix: 1_par_cuda
805     args: -local {{0 1}} -xgpu {{0 1}} -bgpu {{0 1}} -A_mat_type {{mpiaijcusparse mpiaij}} -testnest 0 -testmatmatt 0 -matmatmult_Bbn 3
806 
807   test:
808     output_file: output/ex70_1.out
809     suffix: 2
810     nsize: 1
811     args: -M {{7 11}} -N {{12 9}} -K {{1 3}} -local {{0 1}}
812 
813   testset:
814     requires: cuda
815     output_file: output/ex70_1.out
816     nsize: 1
817     args: -M 7 -N 9 -K 2 -local {{0 1}} -testnest 0 -A_mat_type {{seqdensecuda seqdense}} -xgpu {{0 1}} -bgpu {{0 1}}
818     test:
819       requires: !complex
820       suffix: 2_cuda_real
821     test:
822       # complex+single gives a little bigger error in the MatDenseGetColumnVec test
823       requires: complex !single
824       suffix: 2_cuda_complex
825 
826   test:
827     output_file: output/ex70_1.out
828     suffix: 2_par
829     nsize: 2
830     args: -M {{7 11}} -N {{12 9}} -K {{1 3}} -local {{0 1}} -testcircular -testmatmatt 0
831 
832   test:
833     requires: cuda
834     output_file: output/ex70_1.out
835     suffix: 2_par_cuda
836     nsize: 2
837     args: -M 11 -N 9 -K 1 -local {{0 1}} -testcircular 0 -A_mat_type mpiaijcusparse -xgpu -bgpu -testnest 0 -testmatmatt 0
838 
839   test:
840     output_file: output/ex70_1.out
841     suffix: 3
842     nsize: {{1 3}}
843     args: -M 13 -N 13 -K {{1 3}} -local {{0 1}} -A_mat_type sbaij -symm -testproj 0 -testmatmatt 0
844 
845   test:
846     output_file: output/ex70_1.out
847     suffix: 4
848     nsize: 1
849     args: -M 3 -N 3 -K 3 -local {{0 1}} -testcircular
850 
851   test:
852     output_file: output/ex70_1.out
853     suffix: 5
854     nsize: {{2 4}}
855     args: -M 3 -N 3 -K 3 -local {{0 1}} -testcircular -testmatmatt 0
856 
857   test:
858     output_file: output/ex70_1.out
859     suffix: 6
860     nsize: 1
861     args: -M {{1 3}} -N {{2 5}} -K {{1 2}} -local {{0 1}} -testcircular
862 
863   test:
864     output_file: output/ex70_1.out
865     suffix: 7
866     nsize: 1
867     args: -M 13 -N 13 -K {{1 3}} -local {{0 1}} -A_mat_type dense -testnest -testcircular
868 
869 TEST*/
870