1: static char help[] = "Illustration of MatIS using a 1D Laplacian assembly\n\n";

  3: /*
  4:   MatIS means that the matrix is not assembled. The easiest way to think of this (for me) is that processes do not have
  5:   to hold full matrix rows. One process can hold part of row i, and another processes can hold another part. However, there
  6:   are still the same number of global rows. The local size here is not the size of the local IS block, which we call the
  7:   overlap size, since that is a property only of MatIS. It is the size of the local piece of the vector you multiply in
  8:   MatMult(). This allows PETSc to understand the parallel layout of the Vec, and how it matches the Mat. If you only know
  9:   the overlap size when assembling, it is best to use PETSC_DECIDE for the local size in the creation routine, so that PETSc
 10:   automatically partitions the unknowns.

 12:   Each P_1 element matrix for a cell will be

 14:     /  1 -1 \
 15:     \ -1  1 /

 17:   so that the assembled matrix has a tridiagonal [-1, 2, -1] pattern. We will use 1 cell per process for illustration,
 18:   and allow PETSc to decide the ownership.
 19: */

 21: #include <petscmat.h>

 23: int main(int argc, char **argv)
 24: {
 25:   MPI_Comm               comm;
 26:   Mat                    A;
 27:   Vec                    x, y;
 28:   ISLocalToGlobalMapping map;
 29:   PetscScalar            elemMat[4] = {1.0, -1.0, -1.0, 1.0};
 30:   PetscReal              error;
 31:   PetscInt               overlapSize = 2, globalIdx[2];
 32:   PetscMPIInt            rank, size;

 35:   PetscInitialize(&argc, &argv, NULL, help);
 36:   comm = PETSC_COMM_WORLD;
 37:   MPI_Comm_rank(comm, &rank);
 38:   MPI_Comm_size(comm, &size);
 39:   /* Create local-to-global map */
 40:   globalIdx[0] = rank;
 41:   globalIdx[1] = rank + 1;
 42:   ISLocalToGlobalMappingCreate(comm, 1, overlapSize, globalIdx, PETSC_COPY_VALUES, &map);
 43:   /* Create matrix */
 44:   MatCreateIS(comm, 1, PETSC_DECIDE, PETSC_DECIDE, size + 1, size + 1, map, map, &A);
 45:   PetscObjectSetName((PetscObject)A, "A");
 46:   ISLocalToGlobalMappingDestroy(&map);
 47:   MatISSetPreallocation(A, overlapSize, NULL, overlapSize, NULL);
 48:   MatSetValues(A, 2, globalIdx, 2, globalIdx, elemMat, ADD_VALUES);
 49:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 50:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 51:   /* Check that the constant vector is in the nullspace */
 52:   MatCreateVecs(A, &x, &y);
 53:   VecSet(x, 1.0);
 54:   PetscObjectSetName((PetscObject)x, "x");
 55:   VecViewFromOptions(x, NULL, "-x_view");
 56:   MatMult(A, x, y);
 57:   PetscObjectSetName((PetscObject)y, "y");
 58:   VecViewFromOptions(y, NULL, "-y_view");
 59:   VecNorm(y, NORM_2, &error);
 61:   /* Check that an interior unit vector gets mapped to something of 1-norm 4 */
 62:   if (size > 1) {
 63:     VecSet(x, 0.0);
 64:     VecSetValue(x, 1, 1.0, INSERT_VALUES);
 65:     VecAssemblyBegin(x);
 66:     VecAssemblyEnd(x);
 67:     MatMult(A, x, y);
 68:     VecNorm(y, NORM_1, &error);
 70:   }
 71:   /* Cleanup */
 72:   MatDestroy(&A);
 73:   VecDestroy(&x);
 74:   VecDestroy(&y);
 75:   PetscFinalize();
 76:   return 0;
 77: }

 79: /*TEST

 81:   test:
 82:     suffix: 0
 83:     requires:
 84:     args:

 86:   test:
 87:     suffix: 1
 88:     nsize: 3
 89:     args:

 91:   test:
 92:     suffix: 2
 93:     nsize: 7
 94:     args:

 96: TEST*/