152030a5eSPierre Jolivet static char help[] = "Illustration of MatIS using a 1D Laplacian assembly\n\n"; 2c4762a1bSJed Brown 3c4762a1bSJed Brown /* 4c4762a1bSJed Brown MatIS means that the matrix is not assembled. The easiest way to think of this (for me) is that processes do not have 5c4762a1bSJed Brown to hold full matrix rows. One process can hold part of row i, and another processes can hold another part. However, there 6c4762a1bSJed Brown 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 7c4762a1bSJed Brown overlap size, since that is a property only of MatIS. It is the size of the local piece of the vector you multiply in 8c4762a1bSJed Brown MatMult(). This allows PETSc to understand the parallel layout of the Vec, and how it matches the Mat. If you only know 9c4762a1bSJed Brown the overlap size when assembling, it is best to use PETSC_DECIDE for the local size in the creation routine, so that PETSc 10c4762a1bSJed Brown automatically partitions the unknowns. 11c4762a1bSJed Brown 12c4762a1bSJed Brown Each P_1 element matrix for a cell will be 13c4762a1bSJed Brown 14c4762a1bSJed Brown / 1 -1 \ 15c4762a1bSJed Brown \ -1 1 / 16c4762a1bSJed Brown 17c4762a1bSJed Brown so that the assembled matrix has a tridiagonal [-1, 2, -1] pattern. We will use 1 cell per process for illustration, 18c4762a1bSJed Brown and allow PETSc to decide the ownership. 19c4762a1bSJed Brown */ 20c4762a1bSJed Brown 21c4762a1bSJed Brown #include <petscmat.h> 22c4762a1bSJed Brown 23c4762a1bSJed Brown int main(int argc, char **argv) { 24c4762a1bSJed Brown MPI_Comm comm; 25c4762a1bSJed Brown Mat A; 26c4762a1bSJed Brown Vec x, y; 27c4762a1bSJed Brown ISLocalToGlobalMapping map; 28c4762a1bSJed Brown PetscScalar elemMat[4] = {1.0, -1.0, -1.0, 1.0}; 29c4762a1bSJed Brown PetscReal error; 30c4762a1bSJed Brown PetscInt overlapSize = 2, globalIdx[2]; 31c4762a1bSJed Brown PetscMPIInt rank, size; 32c4762a1bSJed Brown PetscErrorCode ierr; 33c4762a1bSJed Brown 34c4762a1bSJed Brown ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 35c4762a1bSJed Brown comm = PETSC_COMM_WORLD; 36*5f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_rank(comm, &rank)); 37*5f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_size(comm, &size)); 38c4762a1bSJed Brown /* Create local-to-global map */ 39c4762a1bSJed Brown globalIdx[0] = rank; 40c4762a1bSJed Brown globalIdx[1] = rank+1; 41*5f80ce2aSJacob Faibussowitsch CHKERRQ(ISLocalToGlobalMappingCreate(comm, 1, overlapSize, globalIdx, PETSC_COPY_VALUES, &map)); 42c4762a1bSJed Brown /* Create matrix */ 43*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatCreateIS(comm, 1, PETSC_DECIDE, PETSC_DECIDE, size+1, size+1, map, map, &A)); 44*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscObjectSetName((PetscObject) A, "A")); 45*5f80ce2aSJacob Faibussowitsch CHKERRQ(ISLocalToGlobalMappingDestroy(&map)); 46*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatISSetPreallocation(A, overlapSize, NULL, overlapSize, NULL)); 47*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetValues(A, 2, globalIdx, 2, globalIdx, elemMat, ADD_VALUES)); 48*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 49*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 50c4762a1bSJed Brown /* Check that the constant vector is in the nullspace */ 51*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatCreateVecs(A, &x, &y)); 52*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecSet(x, 1.0)); 53*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscObjectSetName((PetscObject) x, "x")); 54*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecViewFromOptions(x, NULL, "-x_view")); 55*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatMult(A, x, y)); 56*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscObjectSetName((PetscObject) y, "y")); 57*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecViewFromOptions(y, NULL, "-y_view")); 58*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecNorm(y, NORM_2, &error)); 592c71b3e2SJacob Faibussowitsch PetscCheckFalse(error > PETSC_SMALL,comm, PETSC_ERR_ARG_WRONG, "Invalid output, x should be in the nullspace of A"); 60c4762a1bSJed Brown /* Check that an interior unit vector gets mapped to something of 1-norm 4 */ 61c4762a1bSJed Brown if (size > 1) { 62*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecSet(x, 0.0)); 63*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecSetValue(x, 1, 1.0, INSERT_VALUES)); 64*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecAssemblyBegin(x)); 65*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecAssemblyEnd(x)); 66*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatMult(A, x, y)); 67*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecNorm(y, NORM_1, &error)); 682c71b3e2SJacob Faibussowitsch PetscCheckFalse(PetscAbsReal(error - 4) > PETSC_SMALL,comm, PETSC_ERR_ARG_WRONG, "Invalid output for matrix multiply"); 69c4762a1bSJed Brown } 70c4762a1bSJed Brown /* Cleanup */ 71*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatDestroy(&A)); 72*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&x)); 73*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&y)); 74c4762a1bSJed Brown ierr = PetscFinalize(); 75c4762a1bSJed Brown return ierr; 76c4762a1bSJed Brown } 77c4762a1bSJed Brown 78c4762a1bSJed Brown /*TEST 79c4762a1bSJed Brown 80c4762a1bSJed Brown test: 81c4762a1bSJed Brown suffix: 0 82c4762a1bSJed Brown requires: 83c4762a1bSJed Brown args: 84c4762a1bSJed Brown 85c4762a1bSJed Brown test: 86c4762a1bSJed Brown suffix: 1 87c4762a1bSJed Brown nsize: 3 88c4762a1bSJed Brown args: 89c4762a1bSJed Brown 90c4762a1bSJed Brown test: 91c4762a1bSJed Brown suffix: 2 92c4762a1bSJed Brown nsize: 7 93c4762a1bSJed Brown args: 94c4762a1bSJed Brown 95c4762a1bSJed Brown TEST*/ 96