Actual source code: ex221.c

  1: static char help[] = "Tests various routines for MATSHELL\n\n";

  3: #include <petscmat.h>

  5: typedef struct _n_User *User;
  6: struct _n_User {
  7:   Mat B;
  8: };

 10: static PetscErrorCode MatGetDiagonal_User(Mat A,Vec X)
 11: {
 12:   User           user;

 14:   MatShellGetContext(A,&user);
 15:   MatGetDiagonal(user->B,X);
 16:   return 0;
 17: }

 19: static PetscErrorCode MatMult_User(Mat A,Vec X,Vec Y)
 20: {
 21:   User           user;

 23:   MatShellGetContext(A,&user);
 24:   MatMult(user->B,X,Y);
 25:   return 0;
 26: }

 28: static PetscErrorCode MatMultTranspose_User(Mat A,Vec X,Vec Y)
 29: {
 30:   User           user;

 32:   MatShellGetContext(A,&user);
 33:   MatMultTranspose(user->B,X,Y);
 34:   return 0;
 35: }

 37: static PetscErrorCode MatCopy_User(Mat A,Mat X,MatStructure str)
 38: {
 39:   User           user,userX;

 41:   MatShellGetContext(A,&user);
 42:   MatShellGetContext(X,&userX);
 44:   PetscObjectReference((PetscObject)user->B);
 45:   return 0;
 46: }

 48: static PetscErrorCode MatDestroy_User(Mat A)
 49: {
 50:   User           user;

 52:   MatShellGetContext(A,&user);
 53:   PetscObjectDereference((PetscObject)user->B);
 54:   return 0;
 55: }

 57: int main(int argc,char **args)
 58: {
 59:   User           user;
 60:   Mat            A,S;
 61:   PetscScalar    *data,diag = 1.3;
 62:   PetscReal      tol = PETSC_SMALL;
 63:   PetscInt       i,j,m = PETSC_DECIDE,n = PETSC_DECIDE,M = 17,N = 15,s1,s2;
 64:   PetscInt       test, ntest = 2;
 65:   PetscMPIInt    rank,size;
 66:   PetscBool      nc = PETSC_FALSE, cong, flg;
 67:   PetscBool      ronl = PETSC_TRUE;
 68:   PetscBool      randomize = PETSC_FALSE, submat = PETSC_FALSE;
 69:   PetscBool      keep = PETSC_FALSE;
 70:   PetscBool      testzerorows = PETSC_TRUE, testdiagscale = PETSC_TRUE, testgetdiag = PETSC_TRUE, testsubmat = PETSC_TRUE;
 71:   PetscBool      testshift = PETSC_TRUE, testscale = PETSC_TRUE, testdup = PETSC_TRUE, testreset = PETSC_TRUE;
 72:   PetscBool      testaxpy = PETSC_TRUE, testaxpyd = PETSC_TRUE, testaxpyerr = PETSC_FALSE;

 74:   PetscInitialize(&argc,&args,(char*)0,help);
 75:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 76:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 77:   PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
 78:   PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
 79:   PetscOptionsGetInt(NULL,NULL,"-ml",&m,NULL);
 80:   PetscOptionsGetInt(NULL,NULL,"-nl",&n,NULL);
 81:   PetscOptionsGetBool(NULL,NULL,"-square_nc",&nc,NULL);
 82:   PetscOptionsGetBool(NULL,NULL,"-rows_only",&ronl,NULL);
 83:   PetscOptionsGetBool(NULL,NULL,"-randomize",&randomize,NULL);
 84:   PetscOptionsGetBool(NULL,NULL,"-submat",&submat,NULL);
 85:   PetscOptionsGetBool(NULL,NULL,"-test_zerorows",&testzerorows,NULL);
 86:   PetscOptionsGetBool(NULL,NULL,"-test_diagscale",&testdiagscale,NULL);
 87:   PetscOptionsGetBool(NULL,NULL,"-test_getdiag",&testgetdiag,NULL);
 88:   PetscOptionsGetBool(NULL,NULL,"-test_shift",&testshift,NULL);
 89:   PetscOptionsGetBool(NULL,NULL,"-test_scale",&testscale,NULL);
 90:   PetscOptionsGetBool(NULL,NULL,"-test_dup",&testdup,NULL);
 91:   PetscOptionsGetBool(NULL,NULL,"-test_reset",&testreset,NULL);
 92:   PetscOptionsGetBool(NULL,NULL,"-test_submat",&testsubmat,NULL);
 93:   PetscOptionsGetBool(NULL,NULL,"-test_axpy",&testaxpy,NULL);
 94:   PetscOptionsGetBool(NULL,NULL,"-test_axpy_different",&testaxpyd,NULL);
 95:   PetscOptionsGetBool(NULL,NULL,"-test_axpy_error",&testaxpyerr,NULL);
 96:   PetscOptionsGetInt(NULL,NULL,"-loop",&ntest,NULL);
 97:   PetscOptionsGetReal(NULL,NULL,"-tol",&tol,NULL);
 98:   PetscOptionsGetScalar(NULL,NULL,"-diag",&diag,NULL);
 99:   PetscOptionsGetBool(NULL,NULL,"-keep",&keep,NULL);
100:   /* This tests square matrices with different row/col layout */
101:   if (nc && size > 1) {
102:     M = PetscMax(PetscMax(N,M),1);
103:     N = M;
104:     m = n = 0;
105:     if (rank == 0) { m = M-1; n = 1; }
106:     else if (rank == 1) { m = 1; n = N-1; }
107:   }
108:   MatCreateDense(PETSC_COMM_WORLD,m,n,M,N,NULL,&A);
109:   MatGetLocalSize(A,&m,&n);
110:   MatGetSize(A,&M,&N);
111:   MatGetOwnershipRange(A,&s1,NULL);
112:   s2   = 1;
113:   while (s2 < M) s2 *= 10;
114:   MatDenseGetArray(A,&data);
115:   for (j = 0; j < N; j++) {
116:     for (i = 0; i < m; i++) {
117:       data[j*m + i] = s2*j + i + s1 + 1;
118:     }
119:   }
120:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
121:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

123:   if (submat) {
124:     Mat      A2;
125:     IS       r,c;
126:     PetscInt rst,ren,cst,cen;

128:     MatGetOwnershipRange(A,&rst,&ren);
129:     MatGetOwnershipRangeColumn(A,&cst,&cen);
130:     ISCreateStride(PetscObjectComm((PetscObject)A),(ren-rst)/2,rst,1,&r);
131:     ISCreateStride(PetscObjectComm((PetscObject)A),(cen-cst)/2,cst,1,&c);
132:     MatCreateSubMatrix(A,r,c,MAT_INITIAL_MATRIX,&A2);
133:     ISDestroy(&r);
134:     ISDestroy(&c);
135:     MatDestroy(&A);
136:     A = A2;
137:   }

139:   MatGetSize(A,&M,&N);
140:   MatGetLocalSize(A,&m,&n);
141:   MatHasCongruentLayouts(A,&cong);

143:   MatConvert(A,MATAIJ,MAT_INPLACE_MATRIX,&A);
144:   MatSetOption(A,MAT_KEEP_NONZERO_PATTERN,keep);
145:   PetscObjectSetName((PetscObject)A,"initial");
146:   MatViewFromOptions(A,NULL,"-view_mat");

148:   PetscNew(&user);
149:   MatCreateShell(PETSC_COMM_WORLD,m,n,M,N,user,&S);
150:   MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MatMult_User);
151:   MatShellSetOperation(S,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_User);
152:   if (cong) {
153:     MatShellSetOperation(S,MATOP_GET_DIAGONAL,(void (*)(void))MatGetDiagonal_User);
154:   }
155:   MatShellSetOperation(S,MATOP_COPY,(void (*)(void))MatCopy_User);
156:   MatShellSetOperation(S,MATOP_DESTROY,(void (*)(void))MatDestroy_User);
157:   MatDuplicate(A,MAT_COPY_VALUES,&user->B);

159:   /* Square and rows only scaling */
160:   ronl = cong ? ronl : PETSC_TRUE;

162:   for (test = 0; test < ntest; test++) {
163:     PetscReal err;

165:     MatMultAddEqual(A,S,10,&flg);
166:     if (!flg) {
167:       PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mult add\n",test);
168:     }
169:     MatMultTransposeAddEqual(A,S,10,&flg);
170:     if (!flg) {
171:       PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mult add (T)\n",test);
172:     }
173:     if (testzerorows) {
174:       Mat       ST,B,C,BT,BTT;
175:       IS        zr;
176:       Vec       x = NULL, b1 = NULL, b2 = NULL;
177:       PetscInt  *idxs = NULL, nr = 0;

179:       if (rank == (test%size)) {
180:         nr = 1;
181:         PetscMalloc1(nr,&idxs);
182:         if (test%2) {
183:           idxs[0] = (2*M - 1 - test/2)%M;
184:         } else {
185:           idxs[0] = (test/2)%M;
186:         }
187:         idxs[0] = PetscMax(idxs[0],0);
188:       }
189:       ISCreateGeneral(PETSC_COMM_WORLD,nr,idxs,PETSC_OWN_POINTER,&zr);
190:       PetscObjectSetName((PetscObject)zr,"ZR");
191:       ISViewFromOptions(zr,NULL,"-view_is");
192:       MatCreateVecs(A,&x,&b1);
193:       if (randomize) {
194:         VecSetRandom(x,NULL);
195:         VecSetRandom(b1,NULL);
196:       } else {
197:         VecSet(x,11.4);
198:         VecSet(b1,-14.2);
199:       }
200:       VecDuplicate(b1,&b2);
201:       VecCopy(b1,b2);
202:       PetscObjectSetName((PetscObject)b1,"A_B1");
203:       PetscObjectSetName((PetscObject)b2,"A_B2");
204:       if (size > 1 && !cong) { /* MATMPIAIJ ZeroRows and ZeroRowsColumns are buggy in this case */
205:         VecDestroy(&b1);
206:       }
207:       if (ronl) {
208:         MatZeroRowsIS(A,zr,diag,x,b1);
209:         MatZeroRowsIS(S,zr,diag,x,b2);
210:       } else {
211:         MatZeroRowsColumnsIS(A,zr,diag,x,b1);
212:         MatZeroRowsColumnsIS(S,zr,diag,x,b2);
213:         ISDestroy(&zr);
214:         /* Mix zerorows and zerorowscols */
215:         nr   = 0;
216:         idxs = NULL;
217:         if (rank == 0) {
218:           nr   = 1;
219:           PetscMalloc1(nr,&idxs);
220:           if (test%2) {
221:             idxs[0] = (3*M - 2 - test/2)%M;
222:           } else {
223:             idxs[0] = (test/2+1)%M;
224:           }
225:           idxs[0] = PetscMax(idxs[0],0);
226:         }
227:         ISCreateGeneral(PETSC_COMM_WORLD,nr,idxs,PETSC_OWN_POINTER,&zr);
228:         PetscObjectSetName((PetscObject)zr,"ZR2");
229:         ISViewFromOptions(zr,NULL,"-view_is");
230:         MatZeroRowsIS(A,zr,diag*2.0+PETSC_SMALL,NULL,NULL);
231:         MatZeroRowsIS(S,zr,diag*2.0+PETSC_SMALL,NULL,NULL);
232:       }
233:       ISDestroy(&zr);

235:       if (b1) {
236:         Vec b;

238:         VecViewFromOptions(b1,NULL,"-view_b");
239:         VecViewFromOptions(b2,NULL,"-view_b");
240:         VecDuplicate(b1,&b);
241:         VecCopy(b1,b);
242:         VecAXPY(b,-1.0,b2);
243:         VecNorm(b,NORM_INFINITY,&err);
244:         if (err >= tol) {
245:           PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error b %g\n",test,(double)err);
246:         }
247:         VecDestroy(&b);
248:       }
249:       VecDestroy(&b1);
250:       VecDestroy(&b2);
251:       VecDestroy(&x);
252:       MatConvert(S,MATDENSE,MAT_INITIAL_MATRIX,&B);

254:       MatCreateTranspose(S,&ST);
255:       MatComputeOperator(ST,MATDENSE,&BT);
256:       MatTranspose(BT,MAT_INITIAL_MATRIX,&BTT);
257:       PetscObjectSetName((PetscObject)B,"S");
258:       PetscObjectSetName((PetscObject)BTT,"STT");
259:       MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&C);
260:       PetscObjectSetName((PetscObject)C,"A");

262:       MatViewFromOptions(C,NULL,"-view_mat");
263:       MatViewFromOptions(B,NULL,"-view_mat");
264:       MatViewFromOptions(BTT,NULL,"-view_mat");

266:       MatAXPY(C,-1.0,B,SAME_NONZERO_PATTERN);
267:       MatNorm(C,NORM_FROBENIUS,&err);
268:       if (err >= tol) {
269:         PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat mult after %s %g\n",test,ronl ? "MatZeroRows" : "MatZeroRowsColumns",(double)err);
270:       }

272:       MatConvert(A,MATDENSE,MAT_REUSE_MATRIX,&C);
273:       MatAXPY(C,-1.0,BTT,SAME_NONZERO_PATTERN);
274:       MatNorm(C,NORM_FROBENIUS,&err);
275:       if (err >= tol) {
276:         PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat mult transpose after %s %g\n",test,ronl ? "MatZeroRows" : "MatZeroRowsColumns",(double)err);
277:       }

279:       MatDestroy(&ST);
280:       MatDestroy(&BTT);
281:       MatDestroy(&BT);
282:       MatDestroy(&B);
283:       MatDestroy(&C);
284:     }
285:     if (testdiagscale) { /* MatDiagonalScale() */
286:       Vec vr,vl;

288:       MatCreateVecs(A,&vr,&vl);
289:       if (randomize) {
290:         VecSetRandom(vr,NULL);
291:         VecSetRandom(vl,NULL);
292:       } else {
293:         VecSet(vr,test%2 ? 0.15 : 1.0/0.15);
294:         VecSet(vl,test%2 ? -1.2 : 1.0/-1.2);
295:       }
296:       MatDiagonalScale(A,vl,vr);
297:       MatDiagonalScale(S,vl,vr);
298:       VecDestroy(&vr);
299:       VecDestroy(&vl);
300:     }

302:     if (testscale) { /* MatScale() */
303:       MatScale(A,test%2 ? 1.4 : 1.0/1.4);
304:       MatScale(S,test%2 ? 1.4 : 1.0/1.4);
305:     }

307:     if (testshift && cong) { /* MatShift() : MATSHELL shift is broken when row/cols layout are not congruent and left/right scaling have been applied */
308:       MatShift(A,test%2 ? -77.5 : 77.5);
309:       MatShift(S,test%2 ? -77.5 : 77.5);
310:     }

312:     if (testgetdiag && cong) { /* MatGetDiagonal() */
313:       Vec dA,dS;

315:       MatCreateVecs(A,&dA,NULL);
316:       MatCreateVecs(S,&dS,NULL);
317:       MatGetDiagonal(A,dA);
318:       MatGetDiagonal(S,dS);
319:       VecAXPY(dA,-1.0,dS);
320:       VecNorm(dA,NORM_INFINITY,&err);
321:       if (err >= tol) {
322:         PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error diag %g\n",test,(double)err);
323:       }
324:       VecDestroy(&dA);
325:       VecDestroy(&dS);
326:     }

328:     if (testdup && !test) {
329:       Mat A2, S2;

331:       MatDuplicate(A,MAT_COPY_VALUES,&A2);
332:       MatDuplicate(S,MAT_COPY_VALUES,&S2);
333:       MatDestroy(&A);
334:       MatDestroy(&S);
335:       A = A2;
336:       S = S2;
337:     }

339:     if (testsubmat) {
340:       Mat      sA,sS,dA,dS,At,St;
341:       IS       r,c;
342:       PetscInt rst,ren,cst,cen;

344:       MatGetOwnershipRange(A,&rst,&ren);
345:       MatGetOwnershipRangeColumn(A,&cst,&cen);
346:       ISCreateStride(PetscObjectComm((PetscObject)A),(ren-rst)/2,rst,1,&r);
347:       ISCreateStride(PetscObjectComm((PetscObject)A),(cen-cst)/2,cst,1,&c);
348:       MatCreateSubMatrix(A,r,c,MAT_INITIAL_MATRIX,&sA);
349:       MatCreateSubMatrix(S,r,c,MAT_INITIAL_MATRIX,&sS);
350:       MatMultAddEqual(sA,sS,10,&flg);
351:       if (!flg) {
352:         PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix mult add\n",test);
353:       }
354:       MatMultTransposeAddEqual(sA,sS,10,&flg);
355:       if (!flg) {
356:         PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix mult add (T)\n",test);
357:       }
358:       MatConvert(sA,MATDENSE,MAT_INITIAL_MATRIX,&dA);
359:       MatConvert(sS,MATDENSE,MAT_INITIAL_MATRIX,&dS);
360:       MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN);
361:       MatNorm(dA,NORM_FROBENIUS,&err);
362:       if (err >= tol) {
363:         PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat submatrix %g\n",test,(double)err);
364:       }
365:       MatDestroy(&sA);
366:       MatDestroy(&sS);
367:       MatDestroy(&dA);
368:       MatDestroy(&dS);
369:       MatCreateTranspose(A,&At);
370:       MatCreateTranspose(S,&St);
371:       MatCreateSubMatrix(At,c,r,MAT_INITIAL_MATRIX,&sA);
372:       MatCreateSubMatrix(St,c,r,MAT_INITIAL_MATRIX,&sS);
373:       MatMultAddEqual(sA,sS,10,&flg);
374:       if (!flg) {
375:         PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix (T) mult add\n",test);
376:       }
377:       MatMultTransposeAddEqual(sA,sS,10,&flg);
378:       if (!flg) {
379:         PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix (T) mult add (T)\n",test);
380:       }
381:       MatConvert(sA,MATDENSE,MAT_INITIAL_MATRIX,&dA);
382:       MatConvert(sS,MATDENSE,MAT_INITIAL_MATRIX,&dS);
383:       MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN);
384:       MatNorm(dA,NORM_FROBENIUS,&err);
385:       if (err >= tol) {
386:         PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat submatrix (T) %g\n",test,(double)err);
387:       }
388:       MatDestroy(&sA);
389:       MatDestroy(&sS);
390:       MatDestroy(&dA);
391:       MatDestroy(&dS);
392:       MatDestroy(&At);
393:       MatDestroy(&St);
394:       ISDestroy(&r);
395:       ISDestroy(&c);
396:     }

398:     if (testaxpy) {
399:       Mat          tA,tS,dA,dS;
400:       MatStructure str[3] = { SAME_NONZERO_PATTERN, SUBSET_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN };

402:       MatDuplicate(A,MAT_COPY_VALUES,&tA);
403:       if (testaxpyd && !(test%2)) {
404:         PetscObjectReference((PetscObject)tA);
405:         tS   = tA;
406:       } else {
407:         PetscObjectReference((PetscObject)S);
408:         tS   = S;
409:       }
410:       MatAXPY(A,0.5,tA,str[test%3]);
411:       MatAXPY(S,0.5,tS,str[test%3]);
412:       /* this will trigger an error the next MatMult or MatMultTranspose call for S */
413:       if (testaxpyerr) MatScale(tA,0);
414:       MatDestroy(&tA);
415:       MatDestroy(&tS);
416:       MatMultAddEqual(A,S,10,&flg);
417:       if (!flg) {
418:         PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error axpy mult add\n",test);
419:       }
420:       MatMultTransposeAddEqual(A,S,10,&flg);
421:       if (!flg) {
422:         PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error axpy mult add (T)\n",test);
423:       }
424:       MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&dA);
425:       MatConvert(S,MATDENSE,MAT_INITIAL_MATRIX,&dS);
426:       MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN);
427:       MatNorm(dA,NORM_FROBENIUS,&err);
428:       if (err >= tol) {
429:         PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat submatrix %g\n",test,(double)err);
430:       }
431:       MatDestroy(&dA);
432:       MatDestroy(&dS);
433:     }

435:     if (testreset && (ntest == 1 || test == ntest-2)) {
436:       /* reset MATSHELL */
437:       MatAssemblyBegin(S,MAT_FINAL_ASSEMBLY);
438:       MatAssemblyEnd(S,MAT_FINAL_ASSEMBLY);
439:       /* reset A */
440:       MatCopy(user->B,A,DIFFERENT_NONZERO_PATTERN);
441:     }
442:   }

444:   MatDestroy(&A);
445:   MatDestroy(&S);
446:   PetscFree(user);
447:   PetscFinalize();
448:   return 0;
449: }

451: /*TEST

453:    testset:
454:      suffix: rect
455:      requires: !single
456:      output_file: output/ex221_1.out
457:      nsize: {{1 3}}
458:      args: -loop 3 -keep {{0 1}} -M {{12 19}} -N {{19 12}} -submat {{0 1}} -test_axpy_different {{0 1}}

460:    testset:
461:      suffix: square
462:      requires: !single
463:      output_file: output/ex221_1.out
464:      nsize: {{1 3}}
465:      args: -M 21 -N 21 -loop 4 -rows_only {{0 1}} -keep {{0 1}} -submat {{0 1}} -test_axpy_different {{0 1}}
466: TEST*/