xref: /petsc/src/mat/tests/ex221.c (revision c4762a1b19cd2af06abeed90e8f9d34fb975dd94) 
1 static char help[] = "Tests various routines for MATSHELL\n\n";
2 
3 #include <petscmat.h>
4 
5 typedef struct _n_User *User;
6 struct _n_User {
7   Mat B;
8 };
9 
10 static PetscErrorCode MatGetDiagonal_User(Mat A,Vec X)
11 {
12   User           user;
13   PetscErrorCode ierr;
14 
15   PetscFunctionBegin;
16   ierr = MatShellGetContext(A,&user);CHKERRQ(ierr);
17   ierr = MatGetDiagonal(user->B,X);CHKERRQ(ierr);
18   PetscFunctionReturn(0);
19 }
20 
21 static PetscErrorCode MatMult_User(Mat A,Vec X,Vec Y)
22 {
23   User           user;
24   PetscErrorCode ierr;
25 
26   PetscFunctionBegin;
27   ierr = MatShellGetContext(A,&user);CHKERRQ(ierr);
28   ierr = MatMult(user->B,X,Y);CHKERRQ(ierr);
29   PetscFunctionReturn(0);
30 }
31 
32 static PetscErrorCode MatMultTranspose_User(Mat A,Vec X,Vec Y)
33 {
34   User           user;
35   PetscErrorCode ierr;
36 
37   PetscFunctionBegin;
38   ierr = MatShellGetContext(A,&user);CHKERRQ(ierr);
39   ierr = MatMultTranspose(user->B,X,Y);CHKERRQ(ierr);
40   PetscFunctionReturn(0);
41 }
42 
43 static PetscErrorCode MatCopy_User(Mat A,Mat X,MatStructure str)
44 {
45   User           user,userX;
46   PetscErrorCode ierr;
47 
48   PetscFunctionBegin;
49   ierr = MatShellGetContext(A,&user);CHKERRQ(ierr);
50   ierr = MatShellGetContext(X,&userX);CHKERRQ(ierr);
51   if (user != userX) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"This should not happen");
52   ierr = PetscObjectReference((PetscObject)user->B);CHKERRQ(ierr);
53   PetscFunctionReturn(0);
54 }
55 
56 static PetscErrorCode MatDestroy_User(Mat A)
57 {
58   User           user;
59   PetscErrorCode ierr;
60 
61   PetscFunctionBegin;
62   ierr = MatShellGetContext(A,&user);CHKERRQ(ierr);
63   ierr = PetscObjectDereference((PetscObject)user->B);CHKERRQ(ierr);
64   PetscFunctionReturn(0);
65 }
66 
67 int main(int argc,char **args)
68 {
69   User           user;
70   Mat            A,S;
71   PetscScalar    *data,diag = 1.3;
72   PetscReal      tol = PETSC_SMALL;
73   PetscInt       i,j,m = PETSC_DECIDE,n = PETSC_DECIDE,M = 17,N = 15,s1,s2;
74   PetscInt       test, ntest = 2;
75   PetscMPIInt    rank,size;
76   PetscBool      nc = PETSC_FALSE, cong, flg;
77   PetscBool      ronl = PETSC_TRUE;
78   PetscBool      randomize = PETSC_FALSE, submat = PETSC_FALSE;
79   PetscBool      keep = PETSC_FALSE;
80   PetscBool      testzerorows = PETSC_TRUE, testdiagscale = PETSC_TRUE, testgetdiag = PETSC_TRUE, testsubmat = PETSC_TRUE;
81   PetscBool      testshift = PETSC_TRUE, testscale = PETSC_TRUE, testdup = PETSC_TRUE, testreset = PETSC_TRUE;
82   PetscBool      testaxpy = PETSC_TRUE, testaxpyd = PETSC_TRUE, testaxpyerr = PETSC_FALSE;
83   PetscErrorCode ierr;
84 
85   ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
86   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
87   ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
88   ierr = PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);CHKERRQ(ierr);
89   ierr = PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);CHKERRQ(ierr);
90   ierr = PetscOptionsGetInt(NULL,NULL,"-ml",&m,NULL);CHKERRQ(ierr);
91   ierr = PetscOptionsGetInt(NULL,NULL,"-nl",&n,NULL);CHKERRQ(ierr);
92   ierr = PetscOptionsGetBool(NULL,NULL,"-square_nc",&nc,NULL);CHKERRQ(ierr);
93   ierr = PetscOptionsGetBool(NULL,NULL,"-rows_only",&ronl,NULL);CHKERRQ(ierr);
94   ierr = PetscOptionsGetBool(NULL,NULL,"-randomize",&randomize,NULL);CHKERRQ(ierr);
95   ierr = PetscOptionsGetBool(NULL,NULL,"-submat",&submat,NULL);CHKERRQ(ierr);
96   ierr = PetscOptionsGetBool(NULL,NULL,"-test_zerorows",&testzerorows,NULL);CHKERRQ(ierr);
97   ierr = PetscOptionsGetBool(NULL,NULL,"-test_diagscale",&testdiagscale,NULL);CHKERRQ(ierr);
98   ierr = PetscOptionsGetBool(NULL,NULL,"-test_getdiag",&testgetdiag,NULL);CHKERRQ(ierr);
99   ierr = PetscOptionsGetBool(NULL,NULL,"-test_shift",&testshift,NULL);CHKERRQ(ierr);
100   ierr = PetscOptionsGetBool(NULL,NULL,"-test_scale",&testscale,NULL);CHKERRQ(ierr);
101   ierr = PetscOptionsGetBool(NULL,NULL,"-test_dup",&testdup,NULL);CHKERRQ(ierr);
102   ierr = PetscOptionsGetBool(NULL,NULL,"-test_reset",&testreset,NULL);CHKERRQ(ierr);
103   ierr = PetscOptionsGetBool(NULL,NULL,"-test_submat",&testsubmat,NULL);CHKERRQ(ierr);
104   ierr = PetscOptionsGetBool(NULL,NULL,"-test_axpy",&testaxpy,NULL);CHKERRQ(ierr);
105   ierr = PetscOptionsGetBool(NULL,NULL,"-test_axpy_different",&testaxpyd,NULL);CHKERRQ(ierr);
106   ierr = PetscOptionsGetBool(NULL,NULL,"-test_axpy_error",&testaxpyerr,NULL);CHKERRQ(ierr);
107   ierr = PetscOptionsGetInt(NULL,NULL,"-loop",&ntest,NULL);CHKERRQ(ierr);
108   ierr = PetscOptionsGetReal(NULL,NULL,"-tol",&tol,NULL);CHKERRQ(ierr);
109   ierr = PetscOptionsGetScalar(NULL,NULL,"-diag",&diag,NULL);CHKERRQ(ierr);
110   ierr = PetscOptionsGetBool(NULL,NULL,"-keep",&keep,NULL);CHKERRQ(ierr);
111   /* This tests square matrices with different row/col layout */
112   if (nc && size > 1) {
113     M = PetscMax(PetscMax(N,M),1);
114     N = M;
115     m = n = 0;
116     if (rank == 0) { m = M-1; n = 1; }
117     else if (rank == 1) { m = 1; n = N-1; }
118   }
119   ierr = MatCreateDense(PETSC_COMM_WORLD,m,n,M,N,NULL,&A);CHKERRQ(ierr);
120   ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr);
121   ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
122   ierr = MatGetOwnershipRange(A,&s1,NULL);CHKERRQ(ierr);
123   s2   = 1;
124   while (s2 < M) s2 *= 10;
125   ierr = MatDenseGetArray(A,&data);CHKERRQ(ierr);
126   for (j = 0; j < N; j++) {
127     for (i = 0; i < m; i++) {
128       data[j*m + i] = s2*j + i + s1 + 1;
129     }
130   }
131   ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
132   ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
133 
134   if (submat) {
135     Mat      A2;
136     IS       r,c;
137     PetscInt rst,ren,cst,cen;
138 
139     ierr = MatGetOwnershipRange(A,&rst,&ren);CHKERRQ(ierr);
140     ierr = MatGetOwnershipRangeColumn(A,&cst,&cen);CHKERRQ(ierr);
141     ierr = ISCreateStride(PetscObjectComm((PetscObject)A),(ren-rst)/2,rst,1,&r);CHKERRQ(ierr);
142     ierr = ISCreateStride(PetscObjectComm((PetscObject)A),(cen-cst)/2,cst,1,&c);CHKERRQ(ierr);
143     ierr = MatCreateSubMatrix(A,r,c,MAT_INITIAL_MATRIX,&A2);CHKERRQ(ierr);
144     ierr = ISDestroy(&r);CHKERRQ(ierr);
145     ierr = ISDestroy(&c);CHKERRQ(ierr);
146     ierr = MatDestroy(&A);CHKERRQ(ierr);
147     A = A2;
148   }
149 
150   ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
151   ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr);
152   ierr = MatHasCongruentLayouts(A,&cong);CHKERRQ(ierr);
153 
154   ierr = MatConvert(A,MATAIJ,MAT_INPLACE_MATRIX,&A);CHKERRQ(ierr);
155   ierr = MatSetOption(A,MAT_KEEP_NONZERO_PATTERN,keep);CHKERRQ(ierr);
156   ierr = PetscObjectSetName((PetscObject)A,"initial");CHKERRQ(ierr);
157   ierr = MatViewFromOptions(A,NULL,"-view_mat");CHKERRQ(ierr);
158 
159   ierr = PetscNew(&user);CHKERRQ(ierr);
160   ierr = MatCreateShell(PETSC_COMM_WORLD,m,n,M,N,user,&S);CHKERRQ(ierr);
161   ierr = MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MatMult_User);CHKERRQ(ierr);
162   ierr = MatShellSetOperation(S,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_User);CHKERRQ(ierr);
163   if (cong) {
164     ierr = MatShellSetOperation(S,MATOP_GET_DIAGONAL,(void (*)(void))MatGetDiagonal_User);CHKERRQ(ierr);
165   }
166   ierr = MatShellSetOperation(S,MATOP_COPY,(void (*)(void))MatCopy_User);CHKERRQ(ierr);
167   ierr = MatShellSetOperation(S,MATOP_DESTROY,(void (*)(void))MatDestroy_User);CHKERRQ(ierr);
168   ierr = MatDuplicate(A,MAT_COPY_VALUES,&user->B);CHKERRQ(ierr);
169 
170   /* Square and rows only scaling */
171   ronl = cong ? ronl : PETSC_TRUE;
172 
173   for (test = 0; test < ntest; test++) {
174     PetscReal err;
175 
176     ierr = MatMultAddEqual(A,S,10,&flg);CHKERRQ(ierr);
177     if (!flg) {
178       ierr = PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error mult add\n",test);CHKERRQ(ierr);
179     }
180     ierr = MatMultTransposeAddEqual(A,S,10,&flg);CHKERRQ(ierr);
181     if (!flg) {
182       ierr = PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error mult add (T)\n",test);CHKERRQ(ierr);
183     }
184     if (testzerorows) {
185       Mat       ST,B,C,BT,BTT;
186       IS        zr;
187       Vec       x = NULL, b1 = NULL, b2 = NULL;
188       PetscInt  *idxs = NULL, nr = 0;
189 
190       if (rank == (test%size)) {
191         nr = 1;
192         ierr = PetscMalloc1(nr,&idxs);CHKERRQ(ierr);
193         if (test%2) {
194           idxs[0] = (2*M - 1 - test/2)%M;
195         } else {
196           idxs[0] = (test/2)%M;
197         }
198         idxs[0] = PetscMax(idxs[0],0);
199       }
200       ierr = ISCreateGeneral(PETSC_COMM_WORLD,nr,idxs,PETSC_OWN_POINTER,&zr);CHKERRQ(ierr);
201       ierr = PetscObjectSetName((PetscObject)zr,"ZR");CHKERRQ(ierr);
202       ierr = ISViewFromOptions(zr,NULL,"-view_is");CHKERRQ(ierr);
203       ierr = MatCreateVecs(A,&x,&b1);CHKERRQ(ierr);
204       if (randomize) {
205         ierr = VecSetRandom(x,NULL);CHKERRQ(ierr);
206         ierr = VecSetRandom(b1,NULL);CHKERRQ(ierr);
207       } else {
208         ierr = VecSet(x,11.4);CHKERRQ(ierr);
209         ierr = VecSet(b1,-14.2);CHKERRQ(ierr);
210       }
211       ierr = VecDuplicate(b1,&b2);CHKERRQ(ierr);
212       ierr = VecCopy(b1,b2);CHKERRQ(ierr);
213       ierr = PetscObjectSetName((PetscObject)b1,"A_B1");CHKERRQ(ierr);
214       ierr = PetscObjectSetName((PetscObject)b2,"A_B2");CHKERRQ(ierr);
215       if (size > 1 && !cong) { /* MATMPIAIJ ZeroRows and ZeroRowsColumns are buggy in this case */
216         ierr = VecDestroy(&b1);CHKERRQ(ierr);
217       }
218       if (ronl) {
219         ierr = MatZeroRowsIS(A,zr,diag,x,b1);CHKERRQ(ierr);
220         ierr = MatZeroRowsIS(S,zr,diag,x,b2);CHKERRQ(ierr);
221       } else {
222         ierr = MatZeroRowsColumnsIS(A,zr,diag,x,b1);CHKERRQ(ierr);
223         ierr = MatZeroRowsColumnsIS(S,zr,diag,x,b2);CHKERRQ(ierr);
224         ierr = ISDestroy(&zr);CHKERRQ(ierr);
225         /* Mix zerorows and zerorowscols */
226         nr   = 0;
227         idxs = NULL;
228         if (!rank) {
229           nr   = 1;
230           ierr = PetscMalloc1(nr,&idxs);CHKERRQ(ierr);
231           if (test%2) {
232             idxs[0] = (3*M - 2 - test/2)%M;
233           } else {
234             idxs[0] = (test/2+1)%M;
235           }
236           idxs[0] = PetscMax(idxs[0],0);
237         }
238         ierr = ISCreateGeneral(PETSC_COMM_WORLD,nr,idxs,PETSC_OWN_POINTER,&zr);CHKERRQ(ierr);
239         ierr = PetscObjectSetName((PetscObject)zr,"ZR2");CHKERRQ(ierr);
240         ierr = ISViewFromOptions(zr,NULL,"-view_is");CHKERRQ(ierr);
241         ierr = MatZeroRowsIS(A,zr,diag*2.0+PETSC_SMALL,NULL,NULL);CHKERRQ(ierr);
242         ierr = MatZeroRowsIS(S,zr,diag*2.0+PETSC_SMALL,NULL,NULL);CHKERRQ(ierr);
243       }
244       ierr = ISDestroy(&zr);CHKERRQ(ierr);
245 
246       if (b1) {
247         Vec b;
248 
249         ierr = VecViewFromOptions(b1,NULL,"-view_b");CHKERRQ(ierr);
250         ierr = VecViewFromOptions(b2,NULL,"-view_b");CHKERRQ(ierr);
251         ierr = VecDuplicate(b1,&b);CHKERRQ(ierr);
252         ierr = VecCopy(b1,b);CHKERRQ(ierr);
253         ierr = VecAXPY(b,-1.0,b2);CHKERRQ(ierr);
254         ierr = VecNorm(b,NORM_INFINITY,&err);CHKERRQ(ierr);
255         if (err >= tol) {
256           ierr = PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error b %g\n",test,(double)err);CHKERRQ(ierr);
257         }
258         ierr = VecDestroy(&b);CHKERRQ(ierr);
259       }
260       ierr = VecDestroy(&b1);CHKERRQ(ierr);
261       ierr = VecDestroy(&b2);CHKERRQ(ierr);
262       ierr = VecDestroy(&x);CHKERRQ(ierr);
263       ierr = MatConvert(S,MATDENSE,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr);
264 
265       ierr = MatCreateTranspose(S,&ST);CHKERRQ(ierr);
266       ierr = MatComputeOperator(ST,MATDENSE,&BT);CHKERRQ(ierr);
267       ierr = MatTranspose(BT,MAT_INITIAL_MATRIX,&BTT);CHKERRQ(ierr);
268       ierr = PetscObjectSetName((PetscObject)B,"S");CHKERRQ(ierr);
269       ierr = PetscObjectSetName((PetscObject)BTT,"STT");CHKERRQ(ierr);
270       ierr = MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr);
271       ierr = PetscObjectSetName((PetscObject)C,"A");CHKERRQ(ierr);
272 
273       ierr = MatViewFromOptions(C,NULL,"-view_mat");CHKERRQ(ierr);
274       ierr = MatViewFromOptions(B,NULL,"-view_mat");CHKERRQ(ierr);
275       ierr = MatViewFromOptions(BTT,NULL,"-view_mat");CHKERRQ(ierr);
276 
277       ierr = MatAXPY(C,-1.0,B,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
278       ierr = MatNorm(C,NORM_FROBENIUS,&err);CHKERRQ(ierr);
279       if (err >= tol) {
280         ierr = PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error mat mult after %s %g\n",test,ronl ? "MatZeroRows" : "MatZeroRowsColumns",(double)err);CHKERRQ(ierr);
281       }
282 
283       ierr = MatConvert(A,MATDENSE,MAT_REUSE_MATRIX,&C);CHKERRQ(ierr);
284       ierr = MatAXPY(C,-1.0,BTT,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
285       ierr = MatNorm(C,NORM_FROBENIUS,&err);CHKERRQ(ierr);
286       if (err >= tol) {
287         ierr = PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error mat mult transpose after %s %g\n",test,ronl ? "MatZeroRows" : "MatZeroRowsColumns",(double)err);CHKERRQ(ierr);
288       }
289 
290       ierr = MatDestroy(&ST);CHKERRQ(ierr);
291       ierr = MatDestroy(&BTT);CHKERRQ(ierr);
292       ierr = MatDestroy(&BT);CHKERRQ(ierr);
293       ierr = MatDestroy(&B);CHKERRQ(ierr);
294       ierr = MatDestroy(&C);CHKERRQ(ierr);
295     }
296     if (testdiagscale) { /* MatDiagonalScale() */
297       Vec vr,vl;
298 
299       ierr = MatCreateVecs(A,&vr,&vl);CHKERRQ(ierr);
300       if (randomize) {
301         ierr = VecSetRandom(vr,NULL);CHKERRQ(ierr);
302         ierr = VecSetRandom(vl,NULL);CHKERRQ(ierr);
303       } else {
304         ierr = VecSet(vr,test%2 ? 0.15 : 1.0/0.15);CHKERRQ(ierr);
305         ierr = VecSet(vl,test%2 ? -1.2 : 1.0/-1.2);CHKERRQ(ierr);
306       }
307       ierr = MatDiagonalScale(A,vl,vr);CHKERRQ(ierr);
308       ierr = MatDiagonalScale(S,vl,vr);CHKERRQ(ierr);
309       ierr = VecDestroy(&vr);CHKERRQ(ierr);
310       ierr = VecDestroy(&vl);CHKERRQ(ierr);
311     }
312 
313     if (testscale) { /* MatScale() */
314       ierr = MatScale(A,test%2 ? 1.4 : 1.0/1.4);CHKERRQ(ierr);
315       ierr = MatScale(S,test%2 ? 1.4 : 1.0/1.4);CHKERRQ(ierr);
316     }
317 
318     if (testshift && cong) { /* MatShift() : MATSHELL shift is broken when row/cols layout are not congruent and left/right scaling have been applied */
319       ierr = MatShift(A,test%2 ? -77.5 : 77.5);CHKERRQ(ierr);
320       ierr = MatShift(S,test%2 ? -77.5 : 77.5);CHKERRQ(ierr);
321     }
322 
323     if (testgetdiag && cong) { /* MatGetDiagonal() */
324       Vec dA,dS;
325 
326       ierr = MatCreateVecs(A,&dA,NULL);CHKERRQ(ierr);
327       ierr = MatCreateVecs(S,&dS,NULL);CHKERRQ(ierr);
328       ierr = MatGetDiagonal(A,dA);CHKERRQ(ierr);
329       ierr = MatGetDiagonal(S,dS);CHKERRQ(ierr);
330       ierr = VecAXPY(dA,-1.0,dS);CHKERRQ(ierr);
331       ierr = VecNorm(dA,NORM_INFINITY,&err);CHKERRQ(ierr);
332       if (err >= tol) {
333         ierr = PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error diag %g\n",test,(double)err);CHKERRQ(ierr);
334       }
335       ierr = VecDestroy(&dA);CHKERRQ(ierr);
336       ierr = VecDestroy(&dS);CHKERRQ(ierr);
337     }
338 
339     if (testdup && !test) {
340       Mat A2, S2;
341 
342       ierr = MatDuplicate(A,MAT_COPY_VALUES,&A2);CHKERRQ(ierr);
343       ierr = MatDuplicate(S,MAT_COPY_VALUES,&S2);CHKERRQ(ierr);
344       ierr = MatDestroy(&A);CHKERRQ(ierr);
345       ierr = MatDestroy(&S);CHKERRQ(ierr);
346       A = A2;
347       S = S2;
348     }
349 
350     if (testsubmat) {
351       Mat      sA,sS,dA,dS,At,St;
352       IS       r,c;
353       PetscInt rst,ren,cst,cen;
354 
355       ierr = MatGetOwnershipRange(A,&rst,&ren);CHKERRQ(ierr);
356       ierr = MatGetOwnershipRangeColumn(A,&cst,&cen);CHKERRQ(ierr);
357       ierr = ISCreateStride(PetscObjectComm((PetscObject)A),(ren-rst)/2,rst,1,&r);CHKERRQ(ierr);
358       ierr = ISCreateStride(PetscObjectComm((PetscObject)A),(cen-cst)/2,cst,1,&c);CHKERRQ(ierr);
359       ierr = MatCreateSubMatrix(A,r,c,MAT_INITIAL_MATRIX,&sA);CHKERRQ(ierr);
360       ierr = MatCreateSubMatrix(S,r,c,MAT_INITIAL_MATRIX,&sS);CHKERRQ(ierr);
361       ierr = MatMultAddEqual(sA,sS,10,&flg);CHKERRQ(ierr);
362       if (!flg) {
363         ierr = PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error submatrix mult add\n",test);CHKERRQ(ierr);
364       }
365       ierr = MatMultTransposeAddEqual(sA,sS,10,&flg);CHKERRQ(ierr);
366       if (!flg) {
367         ierr = PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error submatrix mult add (T)\n",test);CHKERRQ(ierr);
368       }
369       ierr = MatConvert(sA,MATDENSE,MAT_INITIAL_MATRIX,&dA);CHKERRQ(ierr);
370       ierr = MatConvert(sS,MATDENSE,MAT_INITIAL_MATRIX,&dS);CHKERRQ(ierr);
371       ierr = MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
372       ierr = MatNorm(dA,NORM_FROBENIUS,&err);CHKERRQ(ierr);
373       if (err >= tol) {
374         ierr = PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error mat submatrix %g\n",test,(double)err);CHKERRQ(ierr);
375       }
376       ierr = MatDestroy(&sA);CHKERRQ(ierr);
377       ierr = MatDestroy(&sS);CHKERRQ(ierr);
378       ierr = MatDestroy(&dA);CHKERRQ(ierr);
379       ierr = MatDestroy(&dS);CHKERRQ(ierr);
380       ierr = MatCreateTranspose(A,&At);CHKERRQ(ierr);
381       ierr = MatCreateTranspose(S,&St);CHKERRQ(ierr);
382       ierr = MatCreateSubMatrix(At,c,r,MAT_INITIAL_MATRIX,&sA);CHKERRQ(ierr);
383       ierr = MatCreateSubMatrix(St,c,r,MAT_INITIAL_MATRIX,&sS);CHKERRQ(ierr);
384       ierr = MatMultAddEqual(sA,sS,10,&flg);CHKERRQ(ierr);
385       if (!flg) {
386         ierr = PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error submatrix (T) mult add\n",test);CHKERRQ(ierr);
387       }
388       ierr = MatMultTransposeAddEqual(sA,sS,10,&flg);CHKERRQ(ierr);
389       if (!flg) {
390         ierr = PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error submatrix (T) mult add (T)\n",test);CHKERRQ(ierr);
391       }
392       ierr = MatConvert(sA,MATDENSE,MAT_INITIAL_MATRIX,&dA);CHKERRQ(ierr);
393       ierr = MatConvert(sS,MATDENSE,MAT_INITIAL_MATRIX,&dS);CHKERRQ(ierr);
394       ierr = MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
395       ierr = MatNorm(dA,NORM_FROBENIUS,&err);CHKERRQ(ierr);
396       if (err >= tol) {
397         ierr = PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error mat submatrix (T) %g\n",test,(double)err);CHKERRQ(ierr);
398       }
399       ierr = MatDestroy(&sA);CHKERRQ(ierr);
400       ierr = MatDestroy(&sS);CHKERRQ(ierr);
401       ierr = MatDestroy(&dA);CHKERRQ(ierr);
402       ierr = MatDestroy(&dS);CHKERRQ(ierr);
403       ierr = MatDestroy(&At);CHKERRQ(ierr);
404       ierr = MatDestroy(&St);CHKERRQ(ierr);
405       ierr = ISDestroy(&r);CHKERRQ(ierr);
406       ierr = ISDestroy(&c);CHKERRQ(ierr);
407     }
408 
409     if (testaxpy) {
410       Mat          tA,tS,dA,dS;
411       MatStructure str[3] = { SAME_NONZERO_PATTERN, SUBSET_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN };
412 
413       ierr = MatDuplicate(A,MAT_COPY_VALUES,&tA);CHKERRQ(ierr);
414       if (testaxpyd && !(test%2)) {
415         ierr = PetscObjectReference((PetscObject)tA);CHKERRQ(ierr);
416         tS   = tA;
417       } else {
418         ierr = PetscObjectReference((PetscObject)S);CHKERRQ(ierr);
419         tS   = S;
420       }
421       ierr = MatAXPY(A,0.5,tA,str[test%3]);CHKERRQ(ierr);
422       ierr = MatAXPY(S,0.5,tS,str[test%3]);CHKERRQ(ierr);
423       /* this will trigger an error the next MatMult or MatMultTranspose call for S */
424       if (testaxpyerr) { ierr = MatScale(tA,0);CHKERRQ(ierr); }
425       ierr = MatDestroy(&tA);CHKERRQ(ierr);
426       ierr = MatDestroy(&tS);CHKERRQ(ierr);
427       ierr = MatMultAddEqual(A,S,10,&flg);CHKERRQ(ierr);
428       if (!flg) {
429         ierr = PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error axpy mult add\n",test);CHKERRQ(ierr);
430       }
431       ierr = MatMultTransposeAddEqual(A,S,10,&flg);CHKERRQ(ierr);
432       if (!flg) {
433         ierr = PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error axpy mult add (T)\n",test);CHKERRQ(ierr);
434       }
435       ierr = MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&dA);CHKERRQ(ierr);
436       ierr = MatConvert(S,MATDENSE,MAT_INITIAL_MATRIX,&dS);CHKERRQ(ierr);
437       ierr = MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
438       ierr = MatNorm(dA,NORM_FROBENIUS,&err);CHKERRQ(ierr);
439       if (err >= tol) {
440         ierr = PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error mat submatrix %g\n",test,(double)err);CHKERRQ(ierr);
441       }
442       ierr = MatDestroy(&dA);CHKERRQ(ierr);
443       ierr = MatDestroy(&dS);CHKERRQ(ierr);
444     }
445 
446     if (testreset && (ntest == 1 || test == ntest-2)) {
447       /* reset MATSHELL */
448       ierr = MatAssemblyBegin(S,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
449       ierr = MatAssemblyEnd(S,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
450       /* reset A */
451       ierr = MatCopy(user->B,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
452     }
453   }
454 
455   ierr = MatDestroy(&A);CHKERRQ(ierr);
456   ierr = MatDestroy(&S);CHKERRQ(ierr);
457   ierr = PetscFree(user);CHKERRQ(ierr);
458   ierr = PetscFinalize();
459   return ierr;
460 }
461 
462 
463 /*TEST
464 
465    testset:
466      suffix: rect
467      requires: !single
468      output_file: output/ex221_1.out
469      nsize: {{1 3}}
470      args: -loop 3 -keep {{0 1}} -M {{12 19}} -N {{19 12}} -submat {{0 1}} -test_axpy_different {{0 1}}
471 
472    testset:
473      suffix: square
474      requires: !single
475      output_file: output/ex221_1.out
476      nsize: {{1 3}}
477      args: -M 21 -N 21 -loop 4 -rows_only {{0 1}} -keep {{0 1}} -submat {{0 1}} -test_axpy_different {{0 1}}
478 TEST*/
479