Actual source code: ex221.c

petsc-3.12.2 2019-11-22
Report Typos and Errors
  1: static char help[] = "Tests MatZeroRows and MatZeroRowsColumns 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;

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

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

 27:   MatShellGetContext(A,&user);
 28:   MatMult(user->B,X,Y);
 29:   return(0);
 30: }

 32: static PetscErrorCode MatMultTranspose_User(Mat A,Vec X,Vec Y)
 33: {
 34:   User           user;

 38:   MatShellGetContext(A,&user);
 39:   MatMultTranspose(user->B,X,Y);
 40:   return(0);
 41: }

 43: int main(int argc,char **args)
 44: {
 45:   User           user;
 46:   Mat            A,S;
 47:   PetscScalar    *data,diag = 1.3;
 48:   PetscReal      tol = PETSC_SMALL;
 49:   PetscInt       i,j,m = PETSC_DECIDE,n = PETSC_DECIDE,M = 17,N = 15,s1,s2;
 50:   PetscInt       test, ntest = 2;
 51:   PetscMPIInt    rank,size;
 52:   PetscBool      nc = PETSC_FALSE, cong;
 53:   PetscBool      ronl = PETSC_TRUE;
 54:   PetscBool      randomize = PETSC_FALSE;
 55:   PetscBool      keep = PETSC_FALSE;
 56:   PetscBool      testzerorows = PETSC_TRUE, testdiagscale = PETSC_TRUE, testgetdiag = PETSC_TRUE;
 57:   PetscBool      testshift = PETSC_TRUE, testscale = PETSC_TRUE, testdup = PETSC_TRUE, testreset = PETSC_TRUE;

 60:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 61:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 62:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 63:   PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
 64:   PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
 65:   PetscOptionsGetInt(NULL,NULL,"-ml",&m,NULL);
 66:   PetscOptionsGetInt(NULL,NULL,"-nl",&n,NULL);
 67:   PetscOptionsGetBool(NULL,NULL,"-square_nc",&nc,NULL);
 68:   PetscOptionsGetBool(NULL,NULL,"-rows_only",&ronl,NULL);
 69:   PetscOptionsGetBool(NULL,NULL,"-randomize",&randomize,NULL);
 70:   PetscOptionsGetBool(NULL,NULL,"-test_zerorows",&testzerorows,NULL);
 71:   PetscOptionsGetBool(NULL,NULL,"-test_diagscale",&testdiagscale,NULL);
 72:   PetscOptionsGetBool(NULL,NULL,"-test_getdiag",&testgetdiag,NULL);
 73:   PetscOptionsGetBool(NULL,NULL,"-test_shift",&testshift,NULL);
 74:   PetscOptionsGetBool(NULL,NULL,"-test_scale",&testscale,NULL);
 75:   PetscOptionsGetBool(NULL,NULL,"-test_dup",&testdup,NULL);
 76:   PetscOptionsGetBool(NULL,NULL,"-test_reset",&testreset,NULL);
 77:   PetscOptionsGetInt(NULL,NULL,"-loop",&ntest,NULL);
 78:   PetscOptionsGetReal(NULL,NULL,"-tol",&tol,NULL);
 79:   PetscOptionsGetScalar(NULL,NULL,"-diag",&diag,NULL);
 80:   PetscOptionsGetBool(NULL,NULL,"-keep",&keep,NULL);
 81:   /* This tests square matrices with different row/col layout */
 82:   if (nc && size > 1) {
 83:     M = PetscMax(PetscMax(N,M),1);
 84:     N = M;
 85:     m = n = 0;
 86:     if (rank == 0) { m = M-1; n = 1; }
 87:     else if (rank == 1) { m = 1; n = N-1; }
 88:   }
 89:   MatCreateDense(PETSC_COMM_WORLD,m,n,M,N,NULL,&A);
 90:   MatGetLocalSize(A,&m,&n);
 91:   MatGetSize(A,&M,&N);
 92:   MatHasCongruentLayouts(A,&cong);
 93:   MatGetOwnershipRange(A,&s1,NULL);
 94:   s2   = 1;
 95:   while (s2 < M) s2 *= 10;
 96:   MatDenseGetArray(A,&data);
 97:   for (j = 0; j < N; j++) {
 98:     for (i = 0; i < m; i++) {
 99:       data[j*m + i] = s2*j + i + s1 + 1;
100:     }
101:   }
102:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
103:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

105:   MatConvert(A,MATAIJ,MAT_INPLACE_MATRIX,&A);
106:   MatSetOption(A,MAT_KEEP_NONZERO_PATTERN,keep);
107:   PetscObjectSetName((PetscObject)A,"initial");
108:   MatViewFromOptions(A,NULL,"-view_mat");

110:   PetscNew(&user);
111:   MatCreateShell(PETSC_COMM_WORLD,m,n,M,N,user,&S);
112:   MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MatMult_User);
113:   MatShellSetOperation(S,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_User);
114:   if (cong) {
115:     MatShellSetOperation(S,MATOP_GET_DIAGONAL,(void (*)(void))MatGetDiagonal_User);
116:   }
117:   MatDuplicate(A,MAT_COPY_VALUES,&user->B);

119:   /* Square and rows only scaling */
120:   ronl = cong ? ronl : PETSC_TRUE;

122:   for (test = 0; test < ntest; test++) {
123:     PetscReal err;

125:     if (testzerorows) {
126:       Mat       ST,B,C,BT,BTT;
127:       IS        zr;
128:       Vec       x = NULL, b1 = NULL, b2 = NULL;
129:       PetscInt  *idxs = NULL, nr = 0;

131:       if (rank == (test%size)) {
132:         nr = 1;
133:         PetscMalloc1(nr,&idxs);
134:         if (test%2) {
135:           idxs[0] = (2*M - 1 - test/2)%M;
136:         } else {
137:           idxs[0] = (test/2)%M;
138:         }
139:         idxs[0] = PetscMax(idxs[0],0);
140:       }
141:       ISCreateGeneral(PETSC_COMM_WORLD,nr,idxs,PETSC_OWN_POINTER,&zr);
142:       PetscObjectSetName((PetscObject)zr,"ZR");
143:       ISViewFromOptions(zr,NULL,"-view_is");
144:       MatCreateVecs(A,&x,&b1);
145:       if (randomize) {
146:         VecSetRandom(x,NULL);
147:         VecSetRandom(b1,NULL);
148:       } else {
149:         VecSet(x,11.4);
150:         VecSet(b1,-14.2);
151:       }
152:       VecDuplicate(b1,&b2);
153:       VecCopy(b1,b2);
154:       PetscObjectSetName((PetscObject)b1,"A_B1");
155:       PetscObjectSetName((PetscObject)b2,"A_B2");
156:       if (size > 1 && !cong) { /* MATMPIAIJ ZeroRows and ZeroRowsColumns are buggy in this case */
157:         VecDestroy(&b1);
158:       }
159:       if (ronl) {
160:         MatZeroRowsIS(A,zr,diag,x,b1);
161:         MatZeroRowsIS(S,zr,diag,x,b2);
162:       } else {
163:         MatZeroRowsColumnsIS(A,zr,diag,x,b1);
164:         MatZeroRowsColumnsIS(S,zr,diag,x,b2);
165:         ISDestroy(&zr);
166:         /* Mix zerorows and zerorowscols */
167:         nr   = 0;
168:         idxs = NULL;
169:         if (!rank) {
170:           nr   = 1;
171:           PetscMalloc1(nr,&idxs);
172:           if (test%2) {
173:             idxs[0] = (3*M - 2 - test/2)%M;
174:           } else {
175:             idxs[0] = (test/2+1)%M;
176:           }
177:           idxs[0] = PetscMax(idxs[0],0);
178:         }
179:         ISCreateGeneral(PETSC_COMM_WORLD,nr,idxs,PETSC_OWN_POINTER,&zr);
180:         PetscObjectSetName((PetscObject)zr,"ZR2");
181:         ISViewFromOptions(zr,NULL,"-view_is");
182:         MatZeroRowsIS(A,zr,diag*2.0+PETSC_SMALL,NULL,NULL);
183:         MatZeroRowsIS(S,zr,diag*2.0+PETSC_SMALL,NULL,NULL);
184:       }
185:       ISDestroy(&zr);

187:       if (b1) {
188:         Vec b;

190:         VecViewFromOptions(b1,NULL,"-view_b");
191:         VecViewFromOptions(b2,NULL,"-view_b");
192:         VecDuplicate(b1,&b);
193:         VecCopy(b1,b);
194:         VecAXPY(b,-1.0,b2);
195:         VecNorm(b,NORM_INFINITY,&err);
196:         if (err >= tol) {
197:           PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error b %g\n",test,(double)err);
198:         }
199:         VecDestroy(&b);
200:       }
201:       VecDestroy(&b1);
202:       VecDestroy(&b2);
203:       VecDestroy(&x);
204:       MatConvert(S,MATDENSE,MAT_INITIAL_MATRIX,&B);

206:       MatCreateTranspose(S,&ST);
207:       MatComputeOperator(ST,MATDENSE,&BT);
208:       MatTranspose(BT,MAT_INITIAL_MATRIX,&BTT);
209:       PetscObjectSetName((PetscObject)B,"S");
210:       PetscObjectSetName((PetscObject)BTT,"STT");
211:       MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&C);
212:       PetscObjectSetName((PetscObject)C,"A");

214:       MatViewFromOptions(C,NULL,"-view_mat");
215:       MatViewFromOptions(B,NULL,"-view_mat");
216:       MatViewFromOptions(BTT,NULL,"-view_mat");

218:       MatAXPY(C,-1.0,B,SAME_NONZERO_PATTERN);
219:       MatNorm(C,NORM_FROBENIUS,&err);
220:       if (err >= tol) {
221:         PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error mat mult %g\n",test,(double)err);
222:       }

224:       MatConvert(A,MATDENSE,MAT_REUSE_MATRIX,&C);
225:       MatAXPY(C,-1.0,BTT,SAME_NONZERO_PATTERN);
226:       MatNorm(C,NORM_FROBENIUS,&err);
227:       if (err >= tol) {
228:         PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error mat mult transpose %g\n",test,(double)err);
229:       }

231:       MatDestroy(&ST);
232:       MatDestroy(&BTT);
233:       MatDestroy(&BT);
234:       MatDestroy(&B);
235:       MatDestroy(&C);
236:     }
237:     if (testdiagscale) { /* MatDiagonalScale() */
238:       Vec vr,vl;

240:       MatCreateVecs(A,&vr,&vl);
241:       if (randomize) {
242:         VecSetRandom(vr,NULL);
243:         VecSetRandom(vl,NULL);
244:       } else {
245:         VecSet(vr,test%2 ? 0.15 : 1.0/0.15);
246:         VecSet(vl,test%2 ? -1.2 : 1.0/-1.2);
247:       }
248:       MatDiagonalScale(A,vl,vr);
249:       MatDiagonalScale(S,vl,vr);
250:       VecDestroy(&vr);
251:       VecDestroy(&vl);
252:     }

254:     if (testscale) { /* MatScale() */
255:       MatScale(A,test%2 ? 1.4 : 1.0/1.4);
256:       MatScale(S,test%2 ? 1.4 : 1.0/1.4);
257:     }

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

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

267:       MatCreateVecs(A,&dA,NULL);
268:       MatCreateVecs(S,&dS,NULL);
269:       MatGetDiagonal(A,dA);
270:       MatGetDiagonal(S,dS);
271:       VecAXPY(dA,-1.0,dS);
272:       VecNorm(dA,NORM_INFINITY,&err);
273:       if (err >= tol) {
274:         PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error diag %g\n",test,(double)err);
275:       }
276:       VecDestroy(&dA);
277:       VecDestroy(&dS);
278:     }

280:     if (testdup && !test) {
281:       Mat A2, S2;

283:       MatDuplicate(A,MAT_COPY_VALUES,&A2);
284:       MatDuplicate(S,MAT_COPY_VALUES,&S2);
285:       MatDestroy(&A);
286:       MatDestroy(&S);
287:       A = A2;
288:       S = S2;
289:     }

291:     if (testreset && (ntest == 1 || test == ntest-2)) {
292:       /* reset MATSHELL */
293:       MatAssemblyBegin(S,MAT_FINAL_ASSEMBLY);
294:       MatAssemblyEnd(S,MAT_FINAL_ASSEMBLY);
295:       /* reset A */
296:       MatCopy(user->B,A,DIFFERENT_NONZERO_PATTERN);
297:     }
298:   }

300:   MatDestroy(&A);
301:   MatDestroy(&user->B);
302:   MatDestroy(&S);
303:   PetscFree(user);
304:   PetscFinalize();
305:   return ierr;
306: }


309: /*TEST

311:    testset:
312:      suffix: rect
313:      requires: !single
314:      output_file: output/ex221_1.out
315:      nsize: {{1 3}}
316:      args: -keep {{0 1}} -M {{12 19}} -N {{19 12}}

318:    testset:
319:      suffix: square
320:      requires: !single
321:      output_file: output/ex221_1.out
322:      nsize: {{1 3}}
323:      args: -M 21 -N 21 -loop 4 -rows_only {{0 1}} -keep {{0 1}}
324: TEST*/