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