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