Actual source code: ex92.c
petsc-3.4.2 2013-07-02

      2: static char help[] = "Tests MatIncreaseOverlap(), MatGetSubMatrices() for parallel MatSBAIJ format.\n";
      3: /* Example of usage:
      4:       mpiexec -n 2 ./ex92 -nd 2 -ov 3 -mat_block_size 2 -view_id 0 -test_overlap -test_submat
      5: */
      6: #include <petscmat.h>

     10: int main(int argc,char **args)
     11: {
     12:   Mat            A,Atrans,sA,*submatA,*submatsA;
     14:   PetscMPIInt    size,rank;
     15:   PetscInt       bs=1,mbs=10,ov=1,i,j,k,*rows,*cols,nd=2,*idx,rstart,rend,sz,M,N,Mbs;
     16:   PetscScalar    *vals,rval,one=1.0;
     17:   IS             *is1,*is2;
     18:   PetscRandom    rand;
     19:   PetscBool      flg,TestOverlap,TestSubMat,TestAllcols;
     20:   PetscLogStage  stages[2];
     21:   PetscInt       vid = -1;

     23:   PetscInitialize(&argc,&args,(char*)0,help);
     24:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
     25:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

     27:   PetscOptionsGetInt(NULL,"-mat_block_size",&bs,NULL);
     28:   PetscOptionsGetInt(NULL,"-mat_mbs",&mbs,NULL);
     29:   PetscOptionsGetInt(NULL,"-ov",&ov,NULL);
     30:   PetscOptionsGetInt(NULL,"-nd",&nd,NULL);
     31:   PetscOptionsGetInt(NULL,"-view_id",&vid,NULL);
     32:   PetscOptionsHasName(NULL, "-test_overlap", &TestOverlap);
     33:   PetscOptionsHasName(NULL, "-test_submat", &TestSubMat);
     34:   PetscOptionsHasName(NULL, "-test_allcols", &TestAllcols);

     36:   MatCreate(PETSC_COMM_WORLD,&A);
     37:   MatSetSizes(A,mbs*bs,mbs*bs,PETSC_DECIDE,PETSC_DECIDE);
     38:   MatSetType(A,MATBAIJ);
     39:   MatSeqBAIJSetPreallocation(A,bs,PETSC_DEFAULT,NULL);
     40:   MatMPIBAIJSetPreallocation(A,bs,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);

     42:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
     43:   PetscRandomSetFromOptions(rand);

     45:   MatGetOwnershipRange(A,&rstart,&rend);
     46:   MatGetSize(A,&M,&N);
     47:   Mbs  = M/bs;

     49:   PetscMalloc(bs*sizeof(PetscInt),&rows);
     50:   PetscMalloc(bs*sizeof(PetscInt),&cols);
     51:   PetscMalloc(bs*bs*sizeof(PetscScalar),&vals);
     52:   PetscMalloc(M*sizeof(PetscScalar),&idx);

     54:   /* Now set blocks of values */
     55:   for (j=0; j<bs*bs; j++) vals[j] = 0.0;
     56:   for (i=0; i<Mbs; i++) {
     57:     cols[0] = i*bs; rows[0] = i*bs;
     58:     for (j=1; j<bs; j++) {
     59:       rows[j] = rows[j-1]+1;
     60:       cols[j] = cols[j-1]+1;
     61:     }
     62:     MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
     63:   }
     64:   /* second, add random blocks */
     65:   for (i=0; i<20*bs; i++) {
     66:     PetscRandomGetValue(rand,&rval);
     67:     cols[0] = bs*(PetscInt)(PetscRealPart(rval)*Mbs);
     68:     PetscRandomGetValue(rand,&rval);
     69:     rows[0] = rstart + bs*(PetscInt)(PetscRealPart(rval)*mbs);
     70:     for (j=1; j<bs; j++) {
     71:       rows[j] = rows[j-1]+1;
     72:       cols[j] = cols[j-1]+1;
     73:     }

     75:     for (j=0; j<bs*bs; j++) {
     76:       PetscRandomGetValue(rand,&rval);
     77:       vals[j] = rval;
     78:     }
     79:     MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
     80:   }

     82:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
     83:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

     85:   /* make A a symmetric matrix: A <- A^T + A */
     86:   MatTranspose(A,MAT_INITIAL_MATRIX, &Atrans);
     87:   MatAXPY(A,one,Atrans,DIFFERENT_NONZERO_PATTERN);
     88:   MatDestroy(&Atrans);
     89:   MatTranspose(A,MAT_INITIAL_MATRIX, &Atrans);
     90:   MatEqual(A, Atrans, &flg);
     91:   if (flg) {
     92:     MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
     93:   } else SETERRQ(PETSC_COMM_SELF,1,"A+A^T is non-symmetric");
     94:   MatDestroy(&Atrans);

     96:   /* create a SeqSBAIJ matrix sA (= A) */
     97:   MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA);
     98:   if (vid >= 0 && vid < size) {
     99:     if (!rank) printf("A: \n");
    100:     MatView(A,PETSC_VIEWER_STDOUT_WORLD);
    101:     if (!rank) printf("sA: \n");
    102:     MatView(sA,PETSC_VIEWER_STDOUT_WORLD);
    103:   }

    105:   /* Test sA==A through MatMult() */
    106:   MatMultEqual(A,sA,10,&flg);
    107:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Error in MatConvert(): A != sA");

    109:   /* Test MatIncreaseOverlap() */
    110:   PetscMalloc(nd*sizeof(IS **),&is1);
    111:   PetscMalloc(nd*sizeof(IS **),&is2);

    113:   for (i=0; i<nd; i++) {
    114:     if (!TestAllcols) {
    115:       PetscRandomGetValue(rand,&rval);
    116:       sz   = (PetscInt)((0.5+0.2*PetscRealPart(rval))*mbs); /* 0.5*mbs < sz < 0.7*mbs */

    118:       for (j=0; j<sz; j++) {
    119:         PetscRandomGetValue(rand,&rval);
    120:         idx[j*bs] = bs*(PetscInt)(PetscRealPart(rval)*Mbs);
    121:         for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
    122:       }
    123:       ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is1+i);
    124:       ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is2+i);
    125:       if (rank == vid) {
    126:         PetscPrintf(PETSC_COMM_SELF," [%d] IS sz[%d]: %d\n",rank,i,sz);
    127:         ISView(is2[i],PETSC_VIEWER_STDOUT_SELF);
    128:       }
    129:     } else { /* Test all rows and colums */
    130:       sz   = M;
    131:       ISCreateStride(PETSC_COMM_SELF,sz,0,1,is1+i);
    132:       ISCreateStride(PETSC_COMM_SELF,sz,0,1,is2+i);

    134:       if (rank == vid) {
    135:         PetscBool colflag;
    136:         ISIdentity(is2[i],&colflag);
    137:         printf("[%d] is2[%d], colflag %d\n",rank,(int)i,(int)colflag);
    138:         ISView(is2[i],PETSC_VIEWER_STDOUT_SELF);
    139:       }
    140:     }
    141:   }

    143:   PetscLogStageRegister("MatOv_SBAIJ",&stages[0]);
    144:   PetscLogStageRegister("MatOv_BAIJ",&stages[1]);

    146:   /* Test MatIncreaseOverlap */
    147:   if (TestOverlap) {
    148:     PetscLogStagePush(stages[0]);
    149:     MatIncreaseOverlap(sA,nd,is2,ov);
    150:     PetscLogStagePop();

    152:     PetscLogStagePush(stages[1]);
    153:     MatIncreaseOverlap(A,nd,is1,ov);
    154:     PetscLogStagePop();

    156:     if (rank == vid) {
    157:       printf("\n[%d] IS from BAIJ:\n",rank);
    158:       ISView(is1[0],PETSC_VIEWER_STDOUT_SELF);
    159:       printf("\n[%d] IS from SBAIJ:\n",rank);
    160:       ISView(is2[0],PETSC_VIEWER_STDOUT_SELF);
    161:     }

    163:     for (i=0; i<nd; ++i) {
    164:       ISEqual(is1[i],is2[i],&flg);
    165:       if (!flg) {
    166:         if (!rank) {
    167:           ISSort(is1[i]);
    168:           /* ISView(is1[i],PETSC_VIEWER_STDOUT_SELF); */
    169:           ISSort(is2[i]);
    170:           /* ISView(is2[i],PETSC_VIEWER_STDOUT_SELF); */
    171:         }
    172:         SETERRQ1(PETSC_COMM_SELF,1,"i=%D, is1 != is2",i);
    173:       }
    174:     }
    175:   }

    177:   /* Test MatGetSubmatrices */
    178:   if (TestSubMat) {
    179:     for (i = 0; i < nd; ++i) {
    180:       ISSort(is1[i]);
    181:       ISSort(is2[i]);
    182:     }
    183:     MatGetSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
    184:     MatGetSubMatrices(sA,nd,is2,is2,MAT_INITIAL_MATRIX,&submatsA);

    186:     MatMultEqual(A,sA,10,&flg);
    187:     if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"A != sA");

    189:     /* Now test MatGetSubmatrices with MAT_REUSE_MATRIX option */
    190:     MatGetSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
    191:     MatGetSubMatrices(sA,nd,is2,is2,MAT_REUSE_MATRIX,&submatsA);
    192:     MatMultEqual(A,sA,10,&flg);
    193:     if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"MatGetSubmatrices(): A != sA");

    195:     for (i=0; i<nd; ++i) {
    196:       MatDestroy(&submatA[i]);
    197:       MatDestroy(&submatsA[i]);
    198:     }
    199:     PetscFree(submatA);
    200:     PetscFree(submatsA);
    201:   }

    203:   /* Free allocated memory */
    204:   for (i=0; i<nd; ++i) {
    205:     ISDestroy(&is1[i]);
    206:     ISDestroy(&is2[i]);
    207:   }
    208:   PetscFree(is1);
    209:   PetscFree(is2);
    210:   PetscFree(idx);
    211:   PetscFree(rows);
    212:   PetscFree(cols);
    213:   PetscFree(vals);
    214:   MatDestroy(&A);
    215:   MatDestroy(&sA);
    216:   PetscRandomDestroy(&rand);
    217:   PetscFinalize();
    218:   return 0;
    219: }