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

      2: static char help[] = "Tests MatMult(), MatMultAdd(), MatMultTranspose().\n\
      3: Also MatMultTransposeAdd(), MatScale(), MatGetDiagonal(), and MatDiagonalScale().\n\n";

      5:  #include <petscmat.h>

      7: int main(int argc,char **args)
      8: {
      9:   Mat            C;
     10:   Vec            s,u,w,x,y,z;
     12:   PetscInt       i,j,m = 8,n,rstart,rend,vstart,vend;
     13:   PetscScalar    one = 1.0,negone = -1.0,v,alpha=0.1;
     14:   PetscReal      norm, tol = PETSC_SQRT_MACHINE_EPSILON;
     15:   PetscBool      flg;

     17:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
     18:   PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_COMMON);
     19:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
     20:   n    = m;
     21:   PetscOptionsHasName(NULL,NULL,"-rectA",&flg);
     22:   if (flg) n += 2;
     23:   PetscOptionsHasName(NULL,NULL,"-rectB",&flg);
     24:   if (flg) n -= 2;

     26:   /* ---------- Assemble matrix and vectors ----------- */

     28:   MatCreate(PETSC_COMM_WORLD,&C);
     29:   MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m,n);
     30:   MatSetFromOptions(C);
     31:   MatSetUp(C);
     32:   MatGetOwnershipRange(C,&rstart,&rend);
     33:   VecCreate(PETSC_COMM_WORLD,&x);
     34:   VecSetSizes(x,PETSC_DECIDE,m);
     35:   VecSetFromOptions(x);
     36:   VecDuplicate(x,&z);
     37:   VecDuplicate(x,&w);
     38:   VecCreate(PETSC_COMM_WORLD,&y);
     39:   VecSetSizes(y,PETSC_DECIDE,n);
     40:   VecSetFromOptions(y);
     41:   VecDuplicate(y,&u);
     42:   VecDuplicate(y,&s);
     43:   VecGetOwnershipRange(y,&vstart,&vend);

     45:   /* Assembly */
     46:   for (i=rstart; i<rend; i++) {
     47:     v    = 100*(i+1);
     48:     VecSetValues(z,1,&i,&v,INSERT_VALUES);
     49:     for (j=0; j<n; j++) {
     50:       v    = 10*(i+1)+j+1;
     51:       MatSetValues(C,1,&i,1,&j,&v,INSERT_VALUES);
     52:     }
     53:   }

     55:   /* Flush off proc Vec values and do more assembly */
     56:   VecAssemblyBegin(z);
     57:   for (i=vstart; i<vend; i++) {
     58:     v    = one*((PetscReal)i);
     59:     VecSetValues(y,1,&i,&v,INSERT_VALUES);
     60:     v    = 100.0*i;
     61:     VecSetValues(u,1,&i,&v,INSERT_VALUES);
     62:   }

     64:   /* Flush off proc Mat values and do more assembly */
     65:   MatAssemblyBegin(C,MAT_FLUSH_ASSEMBLY);
     66:   for (i=rstart; i<rend; i++) {
     67:     for (j=0; j<n; j++) {
     68:       v    = 10*(i+1)+j+1;
     69:       MatSetValues(C,1,&i,1,&j,&v,INSERT_VALUES);
     70:     }
     71:   }
     72:   /* Try overlap Coomunication with the next stage XXXSetValues */
     73:   VecAssemblyEnd(z);

     75:   MatAssemblyEnd(C,MAT_FLUSH_ASSEMBLY);
     76:   CHKMEMQ;
     77:   /* The Assembly for the second Stage */
     78:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
     79:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
     80:   VecAssemblyBegin(y);
     81:   VecAssemblyEnd(y);
     82:   MatScale(C,alpha);
     83:   VecAssemblyBegin(u);
     84:   VecAssemblyEnd(u);
     85:   PetscPrintf(PETSC_COMM_WORLD,"testing MatMult()\n");
     86:   CHKMEMQ;
     87:   MatMult(C,y,x);
     88:   CHKMEMQ;
     89:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);
     90:   PetscPrintf(PETSC_COMM_WORLD,"testing MatMultAdd()\n");
     91:   MatMultAdd(C,y,z,w);
     92:   VecAXPY(x,one,z);
     93:   VecAXPY(x,negone,w);
     94:   VecNorm(x,NORM_2,&norm);
     95:   if (norm > tol) {
     96:     PetscPrintf(PETSC_COMM_WORLD,"Norm of error difference = %g\n",(double)norm);
     97:   }

     99:   /* ------- Test MatMultTranspose(), MatMultTransposeAdd() ------- */

    101:   for (i=rstart; i<rend; i++) {
    102:     v    = one*((PetscReal)i);
    103:     VecSetValues(x,1,&i,&v,INSERT_VALUES);
    104:   }
    105:   VecAssemblyBegin(x);
    106:   VecAssemblyEnd(x);
    107:   PetscPrintf(PETSC_COMM_WORLD,"testing MatMultTranspose()\n");
    108:   MatMultTranspose(C,x,y);
    109:   VecView(y,PETSC_VIEWER_STDOUT_WORLD);

    111:   PetscPrintf(PETSC_COMM_WORLD,"testing MatMultTransposeAdd()\n");
    112:   MatMultTransposeAdd(C,x,u,s);
    113:   VecAXPY(y,one,u);
    114:   VecAXPY(y,negone,s);
    115:   VecNorm(y,NORM_2,&norm);
    116:   if (norm > tol) {
    117:     PetscPrintf(PETSC_COMM_WORLD,"Norm of error difference = %g\n",(double)norm);
    118:   }

    120:   /* -------------------- Test MatGetDiagonal() ------------------ */

    122:   PetscPrintf(PETSC_COMM_WORLD,"testing MatGetDiagonal(), MatDiagonalScale()\n");
    123:   MatView(C,PETSC_VIEWER_STDOUT_WORLD);
    124:   VecSet(x,one);
    125:   MatGetDiagonal(C,x);
    126:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);
    127:   for (i=vstart; i<vend; i++) {
    128:     v    = one*((PetscReal)(i+1));
    129:     VecSetValues(y,1,&i,&v,INSERT_VALUES);
    130:   }

    132:   /* -------------------- Test () MatDiagonalScale ------------------ */
    133:   PetscOptionsHasName(NULL,NULL,"-test_diagonalscale",&flg);
    134:   if (flg) {
    135:     MatDiagonalScale(C,x,y);
    136:     MatView(C,PETSC_VIEWER_STDOUT_WORLD);
    137:   }
    138:   /* Free data structures */
    139:   VecDestroy(&u); VecDestroy(&s);
    140:   VecDestroy(&w); VecDestroy(&x);
    141:   VecDestroy(&y); VecDestroy(&z);
    142:   MatDestroy(&C);

    144:   PetscFinalize();
    145:   return ierr;
    146: }

    149: /*TEST

    152:    test:
    153:       suffix: 11_A
    154:       args: -mat_type seqaij -rectA
    155:       filter: grep -v type

    157:    test:
    158:       args: -mat_type seqdense -rectA
    159:       suffix: 12_A

    161:    test:
    162:       args: -mat_type seqaij -rectB
    163:       suffix: 11_B
    164:       filter: grep -v type

    166:    test:
    167:       args: -mat_type seqdense -rectB
    168:       suffix: 12_B

    170:    test:
    171:       suffix: 21
    172:       args: -mat_type mpiaij
    173:       filter: grep -v type

    175:    test:
    176:       suffix: 22
    177:       args: -mat_type mpidense

    179:    test:
    180:       suffix: 23
    181:       nsize: 3
    182:       args: -mat_type mpiaij
    183:       filter: grep -v type

    185:    test:
    186:       suffix: 24
    187:       nsize: 3
    188:       args: -mat_type mpidense

    190:    test:
    191:       suffix: 2_aijcusparse_1
    192:       args: -mat_type mpiaijcusparse -vec_type cuda
    193:       filter: grep -v type
    194:       output_file: output/ex5_21.out
    195:       requires: cuda

    197:    test:
    198:       suffix: 2_aijcusparse_2
    199:       nsize: 3
    200:       args: -mat_type mpiaijcusparse -vec_type cuda
    201:       filter: grep -v type
    202:       output_file: output/ex5_23.out
    203:       requires: cuda

    205:    test:
    206:       suffix: 2_aijcusparse_3
    207:       nsize: 3
    208:       args: -mat_type mpiaijcusparse -vec_type cuda -vecscatter_type sf
    209:       filter: grep -v type
    210:       output_file: output/ex5_23.out
    211:       requires: cuda

    213:    test:
    214:       suffix: 31
    215:       args: -mat_type mpiaij -test_diagonalscale
    216:       filter: grep -v type

    218:    test:
    219:       suffix: 32
    220:       args: -mat_type mpibaij -test_diagonalscale
    221:       filter: grep -v Mat_

    223:    test:
    224:       suffix: 33
    225:       nsize: 3
    226:       args: -mat_type mpiaij -test_diagonalscale
    227:       filter: grep -v type

    229:    test:
    230:       suffix: 34
    231:       nsize: 3
    232:       args: -mat_type mpibaij -test_diagonalscale
    233:       filter: grep -v Mat_

    235:    test:
    236:       suffix: 3_aijcusparse_1
    237:       args: -mat_type mpiaijcusparse -vec_type cuda -test_diagonalscale
    238:       filter: grep -v type
    239:       output_file: output/ex5_31.out
    240:       requires: cuda

    242:    test:
    243:       suffix: 3_aijcusparse_2
    244:       nsize: 3
    245:       args: -mat_type mpiaijcusparse -vec_type cuda -test_diagonalscale
    246:       filter: grep -v type
    247:       output_file: output/ex5_33.out
    248:       requires: cuda

    250:    test:
    251:       suffix: aijcusparse_1
    252:       args: -mat_type seqaijcusparse -vec_type cuda -rectA
    253:       filter: grep -v type
    254:       output_file: output/ex5_11_A.out
    255:       requires: cuda

    257:    test:
    258:       suffix: aijcusparse_2
    259:       args: -mat_type seqaijcusparse -vec_type cuda -rectB
    260:       filter: grep -v type
    261:       output_file: output/ex5_11_B.out
    262:       requires: cuda

    264:    test:
    265:       suffix: sell_1
    266:       args: -mat_type sell
    267:       output_file: output/ex5_41.out

    269:    test:
    270:       suffix: sell_2
    271:       nsize: 3
    272:       args: -mat_type sell
    273:       output_file: output/ex5_43.out

    275:    test:
    276:       suffix: sell_3
    277:       args: -mat_type sell -test_diagonalscale
    278:       output_file: output/ex5_51.out

    280:    test:
    281:       suffix: sell_4
    282:       nsize: 3
    283:       args: -mat_type sell -test_diagonalscale
    284:       output_file: output/ex5_53.out

    286: TEST*/