Actual source code: ex184.c

  1: static char help[] = "Example of inverting a block diagonal matrix.\n"
  2:                      "\n";

  4: #include <petscmat.h>

  6: int main(int argc, char **args)
  7: {
  8:   Mat          A, A_inv;
  9:   PetscMPIInt  rank, size;
 10:   PetscInt     M, m, bs, rstart, rend, j, x, y;
 11:   PetscInt    *dnnz;
 12:   PetscScalar *v;
 13:   Vec          X, Y;
 14:   PetscReal    norm;

 16:   PetscFunctionBeginUser;
 17:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 18:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 19:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

 21:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex184", "Mat");
 22:   M = 8;
 23:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_size", &M, NULL));
 24:   bs = 3;
 25:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_block_size", &bs, NULL));
 26:   PetscOptionsEnd();

 28:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 29:   PetscCall(MatSetFromOptions(A));
 30:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M * bs, M * bs));
 31:   PetscCall(MatSetBlockSize(A, bs));
 32:   PetscCall(MatSetUp(A)); /* called so that MatGetLocalSize() will work */
 33:   PetscCall(MatGetLocalSize(A, &m, NULL));
 34:   PetscCall(PetscMalloc1(m / bs, &dnnz));
 35:   for (j = 0; j < m / bs; j++) dnnz[j] = 1;
 36:   PetscCall(MatXAIJSetPreallocation(A, bs, dnnz, NULL, NULL, NULL));
 37:   PetscCall(PetscFree(dnnz));

 39:   PetscCall(PetscMalloc1(bs * bs, &v));
 40:   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
 41:   for (j = rstart / bs; j < rend / bs; j++) {
 42:     for (x = 0; x < bs; x++) {
 43:       for (y = 0; y < bs; y++) {
 44:         if (x == y) {
 45:           v[y + bs * x] = 2 * bs;
 46:         } else {
 47:           v[y + bs * x] = -1 * (x < y) - 2 * (x > y);
 48:         }
 49:       }
 50:     }
 51:     PetscCall(MatSetValuesBlocked(A, 1, &j, 1, &j, v, INSERT_VALUES));
 52:   }
 53:   PetscCall(PetscFree(v));
 54:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 55:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 57:   /* check that A  = inv(inv(A)) */
 58:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A_inv));
 59:   PetscCall(MatSetFromOptions(A_inv));
 60:   PetscCall(MatInvertBlockDiagonalMat(A, A_inv));

 62:   /* Test A_inv * A on a random vector */
 63:   PetscCall(MatCreateVecs(A, &X, &Y));
 64:   PetscCall(VecSetRandom(X, NULL));
 65:   PetscCall(MatMult(A, X, Y));
 66:   PetscCall(VecScale(X, -1));
 67:   PetscCall(MatMultAdd(A_inv, Y, X, X));
 68:   PetscCall(VecNorm(X, NORM_MAX, &norm));
 69:   if (norm > PETSC_SMALL) {
 70:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error exceeds tolerance.\nInverse of block diagonal A\n"));
 71:     PetscCall(MatView(A_inv, PETSC_VIEWER_STDOUT_WORLD));
 72:   }

 74:   PetscCall(MatDestroy(&A));
 75:   PetscCall(MatDestroy(&A_inv));
 76:   PetscCall(VecDestroy(&X));
 77:   PetscCall(VecDestroy(&Y));

 79:   PetscCall(PetscFinalize());
 80:   return 0;
 81: }

 83: /*TEST
 84:   test:
 85:     suffix: seqaij
 86:     args: -mat_type seqaij -mat_size 12 -mat_block_size 3
 87:     nsize: 1
 88:   test:
 89:     suffix: seqbaij
 90:     args: -mat_type seqbaij -mat_size 12 -mat_block_size 3
 91:     nsize: 1
 92:   test:
 93:     suffix: mpiaij
 94:     args: -mat_type mpiaij -mat_size 12 -mat_block_size 3
 95:     nsize: 2
 96:   test:
 97:     suffix: mpibaij
 98:     args: -mat_type mpibaij -mat_size 12 -mat_block_size 3
 99:     nsize: 2
100: TEST*/