1 static char help[] = "Illustration of MatIS using a 1D Laplacian assembly\n\n"; 2 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. 11 12 Each P_1 element matrix for a cell will be 13 14 / 1 -1 \ 15 \ -1 1 / 16 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 */ 20 21 #include <petscmat.h> 22 23 int main(int argc, char **argv) { 24 MPI_Comm comm; 25 Mat A; 26 Vec x, y; 27 ISLocalToGlobalMapping map; 28 PetscScalar elemMat[4] = {1.0, -1.0, -1.0, 1.0}; 29 PetscReal error; 30 PetscInt overlapSize = 2, globalIdx[2]; 31 PetscMPIInt rank, size; 32 PetscErrorCode ierr; 33 34 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 35 comm = PETSC_COMM_WORLD; 36 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 37 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 38 /* Create local-to-global map */ 39 globalIdx[0] = rank; 40 globalIdx[1] = rank+1; 41 ierr = ISLocalToGlobalMappingCreate(comm, 1, overlapSize, globalIdx, PETSC_COPY_VALUES, &map);CHKERRQ(ierr); 42 /* Create matrix */ 43 ierr = MatCreateIS(comm, 1, PETSC_DECIDE, PETSC_DECIDE, size+1, size+1, map, map, &A);CHKERRQ(ierr); 44 ierr = PetscObjectSetName((PetscObject) A, "A");CHKERRQ(ierr); 45 ierr = ISLocalToGlobalMappingDestroy(&map);CHKERRQ(ierr); 46 ierr = MatISSetPreallocation(A, overlapSize, NULL, overlapSize, NULL);CHKERRQ(ierr); 47 ierr = MatSetValues(A, 2, globalIdx, 2, globalIdx, elemMat, ADD_VALUES);CHKERRQ(ierr); 48 ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 49 ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 50 /* Check that the constant vector is in the nullspace */ 51 ierr = MatCreateVecs(A, &x, &y);CHKERRQ(ierr); 52 ierr = VecSet(x, 1.0);CHKERRQ(ierr); 53 ierr = PetscObjectSetName((PetscObject) x, "x");CHKERRQ(ierr); 54 ierr = VecViewFromOptions(x, NULL, "-x_view");CHKERRQ(ierr); 55 ierr = MatMult(A, x, y);CHKERRQ(ierr); 56 ierr = PetscObjectSetName((PetscObject) y, "y");CHKERRQ(ierr); 57 ierr = VecViewFromOptions(y, NULL, "-y_view");CHKERRQ(ierr); 58 ierr = VecNorm(y, NORM_2, &error);CHKERRQ(ierr); 59 if (error > PETSC_SMALL) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid output, x should be in the nullspace of A"); 60 /* Check that an interior unit vector gets mapped to something of 1-norm 4 */ 61 if (size > 1) { 62 ierr = VecSet(x, 0.0);CHKERRQ(ierr); 63 ierr = VecSetValue(x, 1, 1.0, INSERT_VALUES);CHKERRQ(ierr); 64 ierr = VecAssemblyBegin(x);CHKERRQ(ierr); 65 ierr = VecAssemblyEnd(x);CHKERRQ(ierr); 66 ierr = MatMult(A, x, y);CHKERRQ(ierr); 67 ierr = VecNorm(y, NORM_1, &error);CHKERRQ(ierr); 68 if (PetscAbsReal(error - 4) > PETSC_SMALL) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid output for matrix multiply"); 69 } 70 /* Cleanup */ 71 ierr = MatDestroy(&A);CHKERRQ(ierr); 72 ierr = VecDestroy(&x);CHKERRQ(ierr); 73 ierr = VecDestroy(&y);CHKERRQ(ierr); 74 ierr = PetscFinalize(); 75 return ierr; 76 } 77 78 /*TEST 79 80 test: 81 suffix: 0 82 requires: 83 args: 84 85 test: 86 suffix: 1 87 nsize: 3 88 args: 89 90 test: 91 suffix: 2 92 nsize: 7 93 args: 94 95 TEST*/ 96