Actual source code: ex9.c

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

      3: /*
      4:   Include "petscmat.h" so that we can use matrices.
      5:   automatically includes:
      6:      petscsys.h       - base PETSc routines   petscvec.h    - vectors
      7:      petscmat.h    - matrices
      8:      petscis.h     - index sets            petscviewer.h - viewers
      9: */
     10: #include <petscmat.h>

     12: int main(int argc, char **args)
     13: {
     14:   Mat             *A, B; /* matrix */
     15:   Vec              x, y, v, v2, z, z2;
     16:   PetscReal        rnorm;
     17:   PetscInt         n    = 20; /* size of the matrix */
     18:   PetscInt         nmat = 3;  /* number of matrices */
     19:   PetscInt         i;
     20:   PetscRandom      rctx;
     21:   MatCompositeType type;
     22:   PetscScalar      scalings[5] = {2, 3, 4, 5, 6};

     24:   PetscFunctionBeginUser;
     25:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
     26:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
     27:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nmat", &nmat, NULL));

     29:   /*
     30:      Create random matrices
     31:   */
     32:   PetscCall(PetscMalloc1(nmat + 3, &A));
     33:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
     34:   PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, n / 2, 3, NULL, 3, NULL, &A[0]));
     35:   for (i = 1; i < nmat + 1; i++) PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, n, 3, NULL, 3, NULL, &A[i]));
     36:   PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n / 2, n, 3, NULL, 3, NULL, &A[nmat + 1]));
     37:   for (i = 0; i < nmat + 2; i++) PetscCall(MatSetRandom(A[i], rctx));

     39:   PetscCall(MatCreateVecs(A[1], &x, &y));
     40:   PetscCall(VecDuplicate(y, &z));
     41:   PetscCall(VecDuplicate(z, &z2));
     42:   PetscCall(MatCreateVecs(A[0], &v, NULL));
     43:   PetscCall(VecDuplicate(v, &v2));

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

     47:   /* Do MatMult with A[1],A[2],A[3] by hand and store the result in z */
     48:   PetscCall(VecSet(x, 1.0));
     49:   PetscCall(MatMult(A[1], x, z));
     50:   PetscCall(VecScale(z, scalings[1]));
     51:   for (i = 2; i < nmat + 1; i++) {
     52:     PetscCall(MatMult(A[i], x, z2));
     53:     PetscCall(VecAXPY(z, scalings[i], z2));
     54:   }

     56:   /* Do MatMult using MatComposite and store the result in y */
     57:   PetscCall(VecSet(y, 0.0));
     58:   PetscCall(MatCreateComposite(PETSC_COMM_WORLD, nmat, A + 1, &B));
     59:   PetscCall(MatSetFromOptions(B));
     60:   PetscCall(MatCompositeSetScalings(B, &scalings[1]));
     61:   PetscCall(MatMultAdd(B, x, y, y));

     63:   /* Diff y and z */
     64:   PetscCall(VecAXPY(y, -1.0, z));
     65:   PetscCall(VecNorm(y, NORM_2, &rnorm));
     66:   if (rnorm > 10000.0 * PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with composite add %g\n", (double)rnorm));

     68:   /* Test MatCompositeMerge on ADDITIVE MatComposite */
     69:   PetscCall(MatCompositeSetMatStructure(B, DIFFERENT_NONZERO_PATTERN)); /* default */
     70:   PetscCall(MatCompositeMerge(B));
     71:   PetscCall(MatMult(B, x, y));
     72:   PetscCall(MatDestroy(&B));
     73:   PetscCall(VecAXPY(y, -1.0, z));
     74:   PetscCall(VecNorm(y, NORM_2, &rnorm));
     75:   if (rnorm > 10000.0 * PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with composite add after merge %g\n", (double)rnorm));

     77:   /*
     78:      Test n x n/2 multiplicative composite B made up of A[0],A[1],A[2] with separate scalings
     79:   */

     81:   /* Do MatMult with A[0],A[1],A[2] by hand and store the result in z */
     82:   PetscCall(VecSet(v, 1.0));
     83:   PetscCall(MatMult(A[0], v, z));
     84:   PetscCall(VecScale(z, scalings[0]));
     85:   for (i = 1; i < nmat; i++) {
     86:     PetscCall(MatMult(A[i], z, y));
     87:     PetscCall(VecScale(y, scalings[i]));
     88:     PetscCall(VecCopy(y, z));
     89:   }

     91:   /* Do MatMult using MatComposite and store the result in y */
     92:   PetscCall(MatCreateComposite(PETSC_COMM_WORLD, nmat, A, &B));
     93:   PetscCall(MatCompositeSetType(B, MAT_COMPOSITE_MULTIPLICATIVE));
     94:   PetscCall(MatCompositeSetMergeType(B, MAT_COMPOSITE_MERGE_LEFT));
     95:   PetscCall(MatSetFromOptions(B));
     96:   PetscCall(MatCompositeSetScalings(B, &scalings[0]));
     97:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
     98:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); /* do MatCompositeMerge() if -mat_composite_merge 1 */
     99:   PetscCall(MatMult(B, v, y));
    100:   PetscCall(MatDestroy(&B));

    102:   /* Diff y and z */
    103:   PetscCall(VecAXPY(y, -1.0, z));
    104:   PetscCall(VecNorm(y, NORM_2, &rnorm));
    105:   if (rnorm > 10000.0 * PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with composite multiplicative %g\n", (double)rnorm));

    107:   /*
    108:      Test n/2 x n multiplicative composite B made up of A[2], A[3], A[4] without separate scalings
    109:   */
    110:   PetscCall(VecSet(x, 1.0));
    111:   PetscCall(MatMult(A[2], x, z));
    112:   for (i = 3; i < nmat + 1; i++) {
    113:     PetscCall(MatMult(A[i], z, y));
    114:     PetscCall(VecCopy(y, z));
    115:   }
    116:   PetscCall(MatMult(A[nmat + 1], z, v));

    118:   PetscCall(MatCreateComposite(PETSC_COMM_WORLD, nmat, A + 2, &B));
    119:   PetscCall(MatCompositeSetType(B, MAT_COMPOSITE_MULTIPLICATIVE));
    120:   PetscCall(MatSetFromOptions(B));
    121:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    122:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); /* do MatCompositeMerge() if -mat_composite_merge 1 */
    123:   PetscCall(MatMult(B, x, v2));
    124:   PetscCall(MatDestroy(&B));

    126:   PetscCall(VecAXPY(v2, -1.0, v));
    127:   PetscCall(VecNorm(v2, NORM_2, &rnorm));
    128:   if (rnorm > 10000.0 * PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with composite multiplicative %g\n", (double)rnorm));

    130:   /*
    131:      Test get functions
    132:   */
    133:   PetscCall(MatCreateComposite(PETSC_COMM_WORLD, nmat, A, &B));
    134:   PetscCall(MatCompositeGetNumberMat(B, &n));
    135:   if (nmat != n) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with GetNumberMat %" PetscInt_FMT " != %" PetscInt_FMT "\n", nmat, n));
    136:   PetscCall(MatCompositeGetMat(B, 0, &A[nmat + 2]));
    137:   if (A[0] != A[nmat + 2]) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with GetMat\n"));
    138:   PetscCall(MatCompositeGetType(B, &type));
    139:   if (type != MAT_COMPOSITE_ADDITIVE) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with GetType\n"));
    140:   PetscCall(MatDestroy(&B));

    142:   /*
    143:      Free work space.  All PETSc objects should be destroyed when they
    144:      are no longer needed.
    145:   */
    146:   PetscCall(VecDestroy(&x));
    147:   PetscCall(VecDestroy(&y));
    148:   PetscCall(VecDestroy(&v));
    149:   PetscCall(VecDestroy(&v2));
    150:   PetscCall(VecDestroy(&z));
    151:   PetscCall(VecDestroy(&z2));
    152:   PetscCall(PetscRandomDestroy(&rctx));
    153:   for (i = 0; i < nmat + 2; i++) PetscCall(MatDestroy(&A[i]));
    154:   PetscCall(PetscFree(A));

    156:   PetscCall(PetscFinalize());
    157:   return 0;
    158: }

    160: /*TEST

    162:    test:
    163:       nsize: 2
    164:       requires: double
    165:       args: -mat_composite_merge {{0 1}shared output} -mat_composite_merge_mvctx {{0 1}shared output}

    167: TEST*/