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 33*b122ec5aSJacob Faibussowitsch CHKERRQ(PetscInitialize(&argc, &argv, NULL, help)); 34c4762a1bSJed Brown comm = PETSC_COMM_WORLD; 355f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_rank(comm, &rank)); 365f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_size(comm, &size)); 37c4762a1bSJed Brown /* Create local-to-global map */ 38c4762a1bSJed Brown globalIdx[0] = rank; 39c4762a1bSJed Brown globalIdx[1] = rank+1; 405f80ce2aSJacob Faibussowitsch CHKERRQ(ISLocalToGlobalMappingCreate(comm, 1, overlapSize, globalIdx, PETSC_COPY_VALUES, &map)); 41c4762a1bSJed Brown /* Create matrix */ 425f80ce2aSJacob Faibussowitsch CHKERRQ(MatCreateIS(comm, 1, PETSC_DECIDE, PETSC_DECIDE, size+1, size+1, map, map, &A)); 435f80ce2aSJacob Faibussowitsch CHKERRQ(PetscObjectSetName((PetscObject) A, "A")); 445f80ce2aSJacob Faibussowitsch CHKERRQ(ISLocalToGlobalMappingDestroy(&map)); 455f80ce2aSJacob Faibussowitsch CHKERRQ(MatISSetPreallocation(A, overlapSize, NULL, overlapSize, NULL)); 465f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetValues(A, 2, globalIdx, 2, globalIdx, elemMat, ADD_VALUES)); 475f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 485f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 49c4762a1bSJed Brown /* Check that the constant vector is in the nullspace */ 505f80ce2aSJacob Faibussowitsch CHKERRQ(MatCreateVecs(A, &x, &y)); 515f80ce2aSJacob Faibussowitsch CHKERRQ(VecSet(x, 1.0)); 525f80ce2aSJacob Faibussowitsch CHKERRQ(PetscObjectSetName((PetscObject) x, "x")); 535f80ce2aSJacob Faibussowitsch CHKERRQ(VecViewFromOptions(x, NULL, "-x_view")); 545f80ce2aSJacob Faibussowitsch CHKERRQ(MatMult(A, x, y)); 555f80ce2aSJacob Faibussowitsch CHKERRQ(PetscObjectSetName((PetscObject) y, "y")); 565f80ce2aSJacob Faibussowitsch CHKERRQ(VecViewFromOptions(y, NULL, "-y_view")); 575f80ce2aSJacob Faibussowitsch CHKERRQ(VecNorm(y, NORM_2, &error)); 582c71b3e2SJacob Faibussowitsch PetscCheckFalse(error > PETSC_SMALL,comm, PETSC_ERR_ARG_WRONG, "Invalid output, x should be in the nullspace of A"); 59c4762a1bSJed Brown /* Check that an interior unit vector gets mapped to something of 1-norm 4 */ 60c4762a1bSJed Brown if (size > 1) { 615f80ce2aSJacob Faibussowitsch CHKERRQ(VecSet(x, 0.0)); 625f80ce2aSJacob Faibussowitsch CHKERRQ(VecSetValue(x, 1, 1.0, INSERT_VALUES)); 635f80ce2aSJacob Faibussowitsch CHKERRQ(VecAssemblyBegin(x)); 645f80ce2aSJacob Faibussowitsch CHKERRQ(VecAssemblyEnd(x)); 655f80ce2aSJacob Faibussowitsch CHKERRQ(MatMult(A, x, y)); 665f80ce2aSJacob Faibussowitsch CHKERRQ(VecNorm(y, NORM_1, &error)); 672c71b3e2SJacob Faibussowitsch PetscCheckFalse(PetscAbsReal(error - 4) > PETSC_SMALL,comm, PETSC_ERR_ARG_WRONG, "Invalid output for matrix multiply"); 68c4762a1bSJed Brown } 69c4762a1bSJed Brown /* Cleanup */ 705f80ce2aSJacob Faibussowitsch CHKERRQ(MatDestroy(&A)); 715f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&x)); 725f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&y)); 73*b122ec5aSJacob Faibussowitsch CHKERRQ(PetscFinalize()); 74*b122ec5aSJacob Faibussowitsch return 0; 75c4762a1bSJed Brown } 76c4762a1bSJed Brown 77c4762a1bSJed Brown /*TEST 78c4762a1bSJed Brown 79c4762a1bSJed Brown test: 80c4762a1bSJed Brown suffix: 0 81c4762a1bSJed Brown requires: 82c4762a1bSJed Brown args: 83c4762a1bSJed Brown 84c4762a1bSJed Brown test: 85c4762a1bSJed Brown suffix: 1 86c4762a1bSJed Brown nsize: 3 87c4762a1bSJed Brown args: 88c4762a1bSJed Brown 89c4762a1bSJed Brown test: 90c4762a1bSJed Brown suffix: 2 91c4762a1bSJed Brown nsize: 7 92c4762a1bSJed Brown args: 93c4762a1bSJed Brown 94c4762a1bSJed Brown TEST*/ 95