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 PetscFunctionBeginUser; 119 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B)); 120 121 PetscCall(PetscOptionsHasName(NULL, NULL, "-keep_nonzero_pattern", &keepnonzeropattern)); 122 if (keepnonzeropattern) PetscCall(MatSetOption(B, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE)); 123 124 PetscCall(MatZeroRowsIS(B, is, diag, 0, 0)); 125 PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD)); 126 PetscCall(MatDestroy(&B)); 127 PetscFunctionReturn(PETSC_SUCCESS); 128 } 129 130 PetscErrorCode TestMatZeroRows_with_no_allocation(Mat A, IS is, PetscScalar diag) 131 { 132 Mat B; 133 134 /* Now copy A into B, and test it with MatZeroRows() */ 135 PetscFunctionBeginUser; 136 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B)); 137 /* Set this flag after assembly. This way, it affects only MatZeroRows() */ 138 PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); 139 140 PetscCall(MatZeroRowsIS(B, is, diag, 0, 0)); 141 PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD)); 142 PetscCall(MatDestroy(&B)); 143 PetscFunctionReturn(PETSC_SUCCESS); 144 } 145 146 /*TEST 147 148 test: 149 nsize: 2 150 filter: grep -v " MPI process" 151 152 test: 153 suffix: 2 154 nsize: 3 155 args: -mat_type mpibaij -mat_block_size 3 156 filter: grep -v " MPI process" 157 158 test: 159 suffix: 3 160 nsize: 3 161 args: -mat_type mpiaij -keep_nonzero_pattern 162 filter: grep -v " MPI process" 163 164 test: 165 suffix: 4 166 nsize: 3 167 args: -keep_nonzero_pattern -mat_type mpibaij -mat_block_size 3 168 filter: grep -v " MPI process" 169 170 test: 171 requires: !defined(PETSC_HAVE_THREADSAFETY) 172 suffix: 5 173 nsize: 3 174 args: -mat_type mpibaij -mat_block_size 3 -mat_view 175 filter: grep -v " MPI process" 176 177 TEST*/ 178