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