Actual source code: ex230.c
petsc-3.12.4 2020-02-04
Report Typos and Errors

      1: static char help[] = "Example of using MatPreallocator\n\n";

      3: /*T
      4:    Concepts: Mat^matrix preallocation
      5:    Processors: n
      6: T*/

      8:  #include <petscmat.h>

     10: PetscErrorCode ex1_nonsquare_bs1(void)
     11: {
     12:   Mat            A,preallocator;
     13:   PetscInt       M,N,m,n,bs;

     17:   /*
     18:      Create the Jacobian matrix
     19:   */
     20:   M = 10;
     21:   N = 8;
     22:   MatCreate(PETSC_COMM_WORLD,&A);
     23:   MatSetType(A,MATAIJ);
     24:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);
     25:   MatSetBlockSize(A,1);
     26:   MatSetFromOptions(A);

     28:   /*
     29:      Get the sizes of the jacobian.
     30:   */
     31:   MatGetLocalSize(A,&m,&n);
     32:   MatGetBlockSize(A,&bs);

     34:   /*
     35:      Create a preallocator matrix with sizes (local and global) matching the jacobian A.
     36:   */
     37:   MatCreate(PetscObjectComm((PetscObject)A),&preallocator);
     38:   MatSetType(preallocator,MATPREALLOCATOR);
     39:   MatSetSizes(preallocator,m,n,M,N);
     40:   MatSetBlockSize(preallocator,bs);
     41:   MatSetUp(preallocator);

     43:   /*
     44:      Insert non-zero pattern (e.g. perform a sweep over the grid).
     45:      You can use MatSetValues(), MatSetValuesBlocked() or MatSetValue().
     46:   */
     47:   {
     48:     PetscInt    ii,jj;
     49:     PetscScalar vv = 0.0;

     51:     ii = 3; jj = 3;
     52:     MatSetValues(preallocator,1,&ii,1,&jj,&vv,INSERT_VALUES);

     54:     ii = 7; jj = 4;
     55:     MatSetValues(preallocator,1,&ii,1,&jj,&vv,INSERT_VALUES);

     57:     ii = 9; jj = 7;
     58:     MatSetValue(preallocator,ii,jj,vv,INSERT_VALUES);
     59:   }
     60:   MatAssemblyBegin(preallocator,MAT_FINAL_ASSEMBLY);
     61:   MatAssemblyEnd(preallocator,MAT_FINAL_ASSEMBLY);

     63:   /*
     64:      Push the non-zero pattern defined within preallocator into A.
     65:      Internally this will call MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE).
     66:      The arg fill = PETSC_TRUE will insert zeros in the matrix A automatically.
     67:   */
     68:   MatPreallocatorPreallocate(preallocator,PETSC_TRUE,A);

     70:   /*
     71:      We no longer require the preallocator object so destroy it.
     72:   */
     73:   MatDestroy(&preallocator);

     75:   MatView(A,PETSC_VIEWER_STDOUT_WORLD);

     77:   /*
     78:      Insert non-zero values into A.
     79:   */
     80:   {
     81:     PetscInt    ii,jj;
     82:     PetscScalar vv;

     84:     ii = 3; jj = 3; vv = 0.3;
     85:     MatSetValue(A,ii,jj,vv,INSERT_VALUES);

     87:     ii = 7; jj = 4; vv = 3.3;
     88:     MatSetValues(A,1,&ii,1,&jj,&vv,INSERT_VALUES);

     90:     ii = 9; jj = 7; vv = 4.3;
     91:     MatSetValues(A,1,&ii,1,&jj,&vv,INSERT_VALUES);
     92:   }
     93:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
     94:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

     96:   MatView(A,PETSC_VIEWER_STDOUT_WORLD);

     98:   MatDestroy(&A);
     99:   return(0);
    100: }

    102: PetscErrorCode ex2_square_bsvariable(void)
    103: {
    104:   Mat            A,preallocator;
    105:   PetscInt       M,N,m,n,bs = 1;

    109:   PetscOptionsGetInt(NULL,NULL,"-block_size",&bs,NULL);

    111:   /*
    112:      Create the Jacobian matrix.
    113:   */
    114:   M = 10 * bs;
    115:   N = 10 * bs;
    116:   MatCreate(PETSC_COMM_WORLD,&A);
    117:   MatSetType(A,MATAIJ);
    118:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);
    119:   MatSetBlockSize(A,bs);
    120:   MatSetFromOptions(A);

    122:   /*
    123:      Get the sizes of the jacobian.
    124:   */
    125:   MatGetLocalSize(A,&m,&n);
    126:   MatGetBlockSize(A,&bs);

    128:   /*
    129:      Create a preallocator matrix with dimensions matching the jacobian A.
    130:   */
    131:   MatCreate(PetscObjectComm((PetscObject)A),&preallocator);
    132:   MatSetType(preallocator,MATPREALLOCATOR);
    133:   MatSetSizes(preallocator,m,n,M,N);
    134:   MatSetBlockSize(preallocator,bs);
    135:   MatSetUp(preallocator);

    137:   /*
    138:      Insert non-zero pattern (e.g. perform a sweep over the grid).
    139:      You can use MatSetValues(), MatSetValuesBlocked() or MatSetValue().
    140:   */
    141:   {
    142:     PetscInt  ii,jj;
    143:     PetscScalar *vv;

    145:     PetscCalloc1(bs*bs,&vv);

    147:     ii = 0; jj = 9;
    148:     MatSetValue(preallocator,ii,jj,vv[0],INSERT_VALUES);

    150:     ii = 0; jj = 0;
    151:     MatSetValuesBlocked(preallocator,1,&ii,1,&jj,vv,INSERT_VALUES);

    153:     ii = 2; jj = 4;
    154:     MatSetValuesBlocked(preallocator,1,&ii,1,&jj,vv,INSERT_VALUES);

    156:     ii = 4; jj = 2;
    157:     MatSetValuesBlocked(preallocator,1,&ii,1,&jj,vv,INSERT_VALUES);

    159:     ii = 4; jj = 4;
    160:     MatSetValuesBlocked(preallocator,1,&ii,1,&jj,vv,INSERT_VALUES);

    162:     ii = 9; jj = 9;
    163:     MatSetValuesBlocked(preallocator,1,&ii,1,&jj,vv,INSERT_VALUES);

    165:     PetscFree(vv);
    166:   }
    167:   MatAssemblyBegin(preallocator,MAT_FINAL_ASSEMBLY);
    168:   MatAssemblyEnd(preallocator,MAT_FINAL_ASSEMBLY);

    170:   /*
    171:      Push non-zero pattern defined from preallocator into A.
    172:      Internally this will call MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE).
    173:      The arg fill = PETSC_TRUE will insert zeros in the matrix A automatically.
    174:   */
    175:   MatPreallocatorPreallocate(preallocator,PETSC_TRUE,A);

    177:   /*
    178:      We no longer require the preallocator object so destroy it.
    179:   */
    180:   MatDestroy(&preallocator);

    182:   MatView(A,PETSC_VIEWER_STDOUT_WORLD);

    184:   {
    185:     PetscInt    ii,jj;
    186:     PetscScalar *vv;

    188:     PetscCalloc1(bs*bs,&vv);
    189:     for (ii=0; ii<bs*bs; ii++) {
    190:       vv[ii] = (PetscReal)(ii+1);
    191:     }

    193:     ii = 0; jj = 9;
    194:     MatSetValue(A,ii,jj,33.3,INSERT_VALUES);

    196:     ii = 0; jj = 0;
    197:     MatSetValuesBlocked(A,1,&ii,1,&jj,vv,INSERT_VALUES);

    199:     ii = 2; jj = 4;
    200:     MatSetValuesBlocked(A,1,&ii,1,&jj,vv,INSERT_VALUES);

    202:     ii = 4; jj = 2;
    203:     MatSetValuesBlocked(A,1,&ii,1,&jj,vv,INSERT_VALUES);

    205:     ii = 4; jj = 4;
    206:     MatSetValuesBlocked(A,1,&ii,1,&jj,vv,INSERT_VALUES);

    208:     ii = 9; jj = 9;
    209:     MatSetValuesBlocked(A,1,&ii,1,&jj,vv,INSERT_VALUES);

    211:     PetscFree(vv);
    212:   }
    213:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    214:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

    216:   MatView(A,PETSC_VIEWER_STDOUT_WORLD);

    218:   MatDestroy(&A);
    219:   return(0);
    220: }

    222: int main(int argc, char **args)
    223: {
    225:   PetscInt testid = 0;
    226:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
    227:   PetscOptionsGetInt(NULL,NULL,"-test_id",&testid,NULL);
    228:   switch (testid) {
    229:     case 0:
    230:       ex1_nonsquare_bs1();
    231:       break;
    232:     case 1:
    233:       ex2_square_bsvariable();
    234:       break;
    235:     default:
    236:       SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Invalid value for -test_id. Must be {0,1}");
    237:       break;
    238:   }
    239:   PetscFinalize();
    240:   return ierr;
    241: }

    243: /*TEST

    245:    test:
    246:      suffix: t0_a_aij
    247:      nsize: 1
    248:      args: -test_id 0 -mat_type aij

    250:    test:
    251:      suffix: t0_b_aij
    252:      nsize: 6
    253:      args: -test_id 0 -mat_type aij

    256:    test:
    257:      suffix: t1_a_aij
    258:      nsize: 1
    259:      args: -test_id 1 -mat_type aij

    261:    test:
    262:      suffix: t2_a_baij
    263:      nsize: 1
    264:      args: -test_id 1 -mat_type baij

    266:    test:
    267:      suffix: t3_a_sbaij
    268:      nsize: 1
    269:      args: -test_id 1 -mat_type sbaij

    273:    test:
    274:      suffix: t4_a_aij_bs3
    275:      nsize: 1
    276:      args: -test_id 1 -mat_type aij -block_size 3

    278:    test:
    279:      suffix: t5_a_baij_bs3
    280:      nsize: 1
    281:      args: -test_id 1 -mat_type baij -block_size 3

    283:    test:
    284:      suffix: t6_a_sbaij_bs3
    285:      nsize: 1
    286:      args: -test_id 1 -mat_type sbaij -block_size 3

    290:    test:
    291:      suffix: t4_b_aij_bs3
    292:      nsize: 6
    293:      args: -test_id 1 -mat_type aij -block_size 3

    295:    test:
    296:      suffix: t5_b_baij_bs3
    297:      nsize: 6
    298:      args: -test_id 1 -mat_type baij -block_size 3
    299:      filter: grep -v Mat_

    301:    test:
    302:      suffix: t6_b_sbaij_bs3
    303:      nsize: 6
    304:      args: -test_id 1 -mat_type sbaij -block_size 3
    305:      filter: grep -v Mat_

    307: TEST*/