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