Actual source code: ex221.c
petsc-3.12.4 2020-02-04
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*/