1 2 static char help[] = "Tests the use of MatZeroRows() for parallel matrices.\n\ 3 This example also tests the use of MatDuplicate() for both MPIAIJ and MPIBAIJ matrices"; 4 5 #include <petscmat.h> 6 7 extern PetscErrorCode TestMatZeroRows_Basic(Mat, IS, PetscScalar); 8 extern PetscErrorCode TestMatZeroRows_with_no_allocation(Mat, IS, PetscScalar); 9 10 int main(int argc, char **args) 11 { 12 Mat A; 13 PetscInt i, j, m = 3, n, Ii, J, Imax; 14 PetscMPIInt rank, size; 15 PetscScalar v, diag = -4.0; 16 IS is; 17 18 PetscFunctionBeginUser; 19 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 20 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 21 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 22 n = 2 * size; 23 24 /* create A Square matrix for the five point stencil,YET AGAIN*/ 25 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 26 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n)); 27 PetscCall(MatSetFromOptions(A)); 28 PetscCall(MatSetUp(A)); 29 for (i = 0; i < m; i++) { 30 for (j = 2 * rank; j < 2 * rank + 2; j++) { 31 v = -1.0; 32 Ii = j + n * i; 33 if (i > 0) { 34 J = Ii - n; 35 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 36 } 37 if (i < m - 1) { 38 J = Ii + n; 39 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 40 } 41 if (j > 0) { 42 J = Ii - 1; 43 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 44 } 45 if (j < n - 1) { 46 J = Ii + 1; 47 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 48 } 49 v = 4.0; 50 PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES)); 51 } 52 } 53 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 54 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 55 56 /* Create AN IS required by MatZeroRows() */ 57 Imax = n * rank; 58 if (Imax >= n * m - m - 1) Imax = m * n - m - 1; 59 PetscCall(ISCreateStride(PETSC_COMM_SELF, m, Imax, 1, &is)); 60 61 PetscCall(TestMatZeroRows_Basic(A, is, 0.0)); 62 PetscCall(TestMatZeroRows_Basic(A, is, diag)); 63 64 PetscCall(TestMatZeroRows_with_no_allocation(A, is, 0.0)); 65 PetscCall(TestMatZeroRows_with_no_allocation(A, is, diag)); 66 67 PetscCall(MatDestroy(&A)); 68 69 /* Now Create a rectangular matrix with five point stencil (app) 70 n+size is used so that this dimension is always divisible by size. 71 This way, we can always use bs = size for any number of procs */ 72 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 73 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n, m * (n + size))); 74 PetscCall(MatSetFromOptions(A)); 75 PetscCall(MatSetUp(A)); 76 for (i = 0; i < m; i++) { 77 for (j = 2 * rank; j < 2 * rank + 2; j++) { 78 v = -1.0; 79 Ii = j + n * i; 80 if (i > 0) { 81 J = Ii - n; 82 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 83 } 84 if (i < m - 1) { 85 J = Ii + n; 86 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 87 } 88 if (j > 0) { 89 J = Ii - 1; 90 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 91 } 92 if (j < n + size - 1) { 93 J = Ii + 1; 94 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 95 } 96 v = 4.0; 97 PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES)); 98 } 99 } 100 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 101 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 102 103 PetscCall(TestMatZeroRows_Basic(A, is, 0.0)); 104 PetscCall(TestMatZeroRows_Basic(A, is, diag)); 105 106 PetscCall(MatDestroy(&A)); 107 PetscCall(ISDestroy(&is)); 108 PetscCall(PetscFinalize()); 109 return 0; 110 } 111 112 PetscErrorCode TestMatZeroRows_Basic(Mat A, IS is, PetscScalar diag) 113 { 114 Mat B; 115 PetscBool keepnonzeropattern; 116 117 /* Now copy A into B, and test it with MatZeroRows() */ 118 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B)); 119 120 PetscCall(PetscOptionsHasName(NULL, NULL, "-keep_nonzero_pattern", &keepnonzeropattern)); 121 if (keepnonzeropattern) PetscCall(MatSetOption(B, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE)); 122 123 PetscCall(MatZeroRowsIS(B, is, diag, 0, 0)); 124 PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD)); 125 PetscCall(MatDestroy(&B)); 126 return 0; 127 } 128 129 PetscErrorCode TestMatZeroRows_with_no_allocation(Mat A, IS is, PetscScalar diag) 130 { 131 Mat B; 132 133 /* Now copy A into B, and test it with MatZeroRows() */ 134 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B)); 135 /* Set this flag after assembly. This way, it affects only MatZeroRows() */ 136 PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); 137 138 PetscCall(MatZeroRowsIS(B, is, diag, 0, 0)); 139 PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD)); 140 PetscCall(MatDestroy(&B)); 141 return 0; 142 } 143 144 /*TEST 145 146 test: 147 nsize: 2 148 filter: grep -v " MPI process" 149 150 test: 151 suffix: 2 152 nsize: 3 153 args: -mat_type mpibaij -mat_block_size 3 154 filter: grep -v " MPI process" 155 156 test: 157 suffix: 3 158 nsize: 3 159 args: -mat_type mpiaij -keep_nonzero_pattern 160 filter: grep -v " MPI process" 161 162 test: 163 suffix: 4 164 nsize: 3 165 args: -keep_nonzero_pattern -mat_type mpibaij -mat_block_size 3 166 filter: grep -v " MPI process" 167 168 TEST*/ 169