Actual source code: ex9.c

      2: static char help[] = "Tests MatCreateComposite()\n\n";

      4: /*T
      5:    Concepts: Mat^composite matrices
      6:    Processors: n
      7: T*/

      9: /*
     10:   Include "petscmat.h" so that we can use matrices.
     11:   automatically includes:
     12:      petscsys.h       - base PETSc routines   petscvec.h    - vectors
     13:      petscmat.h    - matrices
     14:      petscis.h     - index sets            petscviewer.h - viewers
     15: */
     16: #include <petscmat.h>

     18: int main(int argc,char **args)
     19: {
     20:   Mat              *A,B;           /* matrix */
     21:   PetscErrorCode   ierr;
     22:   Vec              x,y,v,v2,z,z2;
     23:   PetscReal        rnorm;
     24:   PetscInt         n = 20;         /* size of the matrix */
     25:   PetscInt         nmat = 3;       /* number of matrices */
     26:   PetscInt         i;
     27:   PetscRandom      rctx;
     28:   MatCompositeType type;
     29:   PetscScalar      scalings[5]={2,3,4,5,6};

     31:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
     32:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
     33:   PetscOptionsGetInt(NULL,NULL,"-nmat",&nmat,NULL);

     35:   /*
     36:      Create random matrices
     37:   */
     38:   PetscMalloc1(nmat+3,&A);
     39:   PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
     40:   MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n/2,3,NULL,3,NULL,&A[0]);
     41:   for (i = 1; i < nmat+1; i++) {
     42:     MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,3,NULL,3,NULL,&A[i]);
     43:   }
     44:   MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n/2,n,3,NULL,3,NULL,&A[nmat+1]);
     45:   for (i = 0; i < nmat+2; i++) {
     46:     MatSetRandom(A[i],rctx);
     47:   }

     49:   MatCreateVecs(A[1],&x,&y);
     50:   VecDuplicate(y,&z);
     51:   VecDuplicate(z,&z2);
     52:   MatCreateVecs(A[0],&v,NULL);
     53:   VecDuplicate(v,&v2);

     55:   /* Test MatMult of an ADDITIVE MatComposite B made up of A[1],A[2],A[3] with separate scalings */

     57:   /* Do MatMult with A[1],A[2],A[3] by hand and store the result in z */
     58:   VecSet(x,1.0);
     59:   MatMult(A[1],x,z);
     60:   VecScale(z,scalings[1]);
     61:   for (i = 2; i < nmat+1; i++) {
     62:     MatMult(A[i],x,z2);
     63:     VecAXPY(z,scalings[i],z2);
     64:   }

     66:   /* Do MatMult using MatComposite and store the result in y */
     67:   VecSet(y,0.0);
     68:   MatCreateComposite(PETSC_COMM_WORLD,nmat,A+1,&B);
     69:   MatSetFromOptions(B);
     70:   MatCompositeSetScalings(B,&scalings[1]);
     71:   MatMultAdd(B,x,y,y);

     73:   /* Diff y and z */
     74:   VecAXPY(y,-1.0,z);
     75:   VecNorm(y,NORM_2,&rnorm);
     76:   if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
     77:     PetscPrintf(PETSC_COMM_WORLD,"Error with composite add %g\n",(double)rnorm);
     78:   }

     80:   /* Test MatCompositeMerge on ADDITIVE MatComposite */
     81:   MatCompositeSetMatStructure(B,DIFFERENT_NONZERO_PATTERN); /* default */
     82:   MatCompositeMerge(B);
     83:   MatMult(B,x,y);
     84:   MatDestroy(&B);
     85:   VecAXPY(y,-1.0,z);
     86:   VecNorm(y,NORM_2,&rnorm);
     87:   if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
     88:     PetscPrintf(PETSC_COMM_WORLD,"Error with composite add after merge %g\n",(double)rnorm);
     89:   }

     91:   /*
     92:      Test n x n/2 multiplicative composite B made up of A[0],A[1],A[2] with separate scalings
     93:   */

     95:   /* Do MatMult with A[0],A[1],A[2] by hand and store the result in z */
     96:   VecSet(v,1.0);
     97:   MatMult(A[0],v,z);
     98:   VecScale(z,scalings[0]);
     99:   for (i = 1; i < nmat; i++) {
    100:     MatMult(A[i],z,y);
    101:     VecScale(y,scalings[i]);
    102:     VecCopy(y,z);
    103:   }

    105:   /* Do MatMult using MatComposite and store the result in y */
    106:   MatCreateComposite(PETSC_COMM_WORLD,nmat,A,&B);
    107:   MatCompositeSetType(B,MAT_COMPOSITE_MULTIPLICATIVE);
    108:   MatCompositeSetMergeType(B,MAT_COMPOSITE_MERGE_LEFT);
    109:   MatSetFromOptions(B);
    110:   MatCompositeSetScalings(B,&scalings[0]);
    111:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    112:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); /* do MatCompositeMerge() if -mat_composite_merge 1 */
    113:   MatMult(B,v,y);
    114:   MatDestroy(&B);

    116:   /* Diff y and z */
    117:   VecAXPY(y,-1.0,z);
    118:   VecNorm(y,NORM_2,&rnorm);
    119:   if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
    120:     PetscPrintf(PETSC_COMM_WORLD,"Error with composite multiplicative %g\n",(double)rnorm);
    121:   }

    123:   /*
    124:      Test n/2 x n multiplicative composite B made up of A[2], A[3], A[4] without separate scalings
    125:   */
    126:   VecSet(x,1.0);
    127:   MatMult(A[2],x,z);
    128:   for (i = 3; i < nmat+1; i++) {
    129:     MatMult(A[i],z,y);
    130:     VecCopy(y,z);
    131:   }
    132:   MatMult(A[nmat+1],z,v);

    134:   MatCreateComposite(PETSC_COMM_WORLD,nmat,A+2,&B);
    135:   MatCompositeSetType(B,MAT_COMPOSITE_MULTIPLICATIVE);
    136:   MatSetFromOptions(B);
    137:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    138:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); /* do MatCompositeMerge() if -mat_composite_merge 1 */
    139:   MatMult(B,x,v2);
    140:   MatDestroy(&B);

    142:   VecAXPY(v2,-1.0,v);
    143:   VecNorm(v2,NORM_2,&rnorm);
    144:   if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
    145:     PetscPrintf(PETSC_COMM_WORLD,"Error with composite multiplicative %g\n",(double)rnorm);
    146:   }

    148:   /*
    149:      Test get functions
    150:   */
    151:   MatCreateComposite(PETSC_COMM_WORLD,nmat,A,&B);
    152:   MatCompositeGetNumberMat(B,&n);
    153:   if (nmat != n) {
    154:     PetscPrintf(PETSC_COMM_WORLD,"Error with GetNumberMat %D != %D\n",nmat,n);
    155:   }
    156:   MatCompositeGetMat(B,0,&A[nmat+2]);
    157:   if (A[0] != A[nmat+2]) {
    158:     PetscPrintf(PETSC_COMM_WORLD,"Error with GetMat\n");
    159:   }
    160:   MatCompositeGetType(B,&type);
    161:   if (type != MAT_COMPOSITE_ADDITIVE) {
    162:     PetscPrintf(PETSC_COMM_WORLD,"Error with GetType\n");
    163:   }
    164:   MatDestroy(&B);

    166:   /*
    167:      Free work space.  All PETSc objects should be destroyed when they
    168:      are no longer needed.
    169:   */
    170:   VecDestroy(&x);
    171:   VecDestroy(&y);
    172:   VecDestroy(&v);
    173:   VecDestroy(&v2);
    174:   VecDestroy(&z);
    175:   VecDestroy(&z2);
    176:   PetscRandomDestroy(&rctx);
    177:   for (i = 0; i < nmat+2; i++) {
    178:     MatDestroy(&A[i]);
    179:   }
    180:   PetscFree(A);

    182:   PetscFinalize();
    183:   return ierr;
    184: }

    186: /*TEST

    188:    test:
    189:       nsize: 2
    190:       requires: double
    191:       args: -mat_composite_merge {{0 1}shared output} -mat_composite_merge_mvctx {{0 1}shared output}

    193: TEST*/