Actual source code: ex12.c


  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";

  5: #include <petscmat.h>

  7: extern PetscErrorCode TestMatZeroRows_Basic(Mat, IS, PetscScalar);
  8: extern PetscErrorCode TestMatZeroRows_with_no_allocation(Mat, IS, PetscScalar);

 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;

 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;

 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));

 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));

 61:   PetscCall(TestMatZeroRows_Basic(A, is, 0.0));
 62:   PetscCall(TestMatZeroRows_Basic(A, is, diag));

 64:   PetscCall(TestMatZeroRows_with_no_allocation(A, is, 0.0));
 65:   PetscCall(TestMatZeroRows_with_no_allocation(A, is, diag));

 67:   PetscCall(MatDestroy(&A));

 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));

103:   PetscCall(TestMatZeroRows_Basic(A, is, 0.0));
104:   PetscCall(TestMatZeroRows_Basic(A, is, diag));

106:   PetscCall(MatDestroy(&A));
107:   PetscCall(ISDestroy(&is));
108:   PetscCall(PetscFinalize());
109:   return 0;
110: }

112: PetscErrorCode TestMatZeroRows_Basic(Mat A, IS is, PetscScalar diag)
113: {
114:   Mat       B;
115:   PetscBool keepnonzeropattern;

117:   /* Now copy A into B, and test it with MatZeroRows() */
118:   PetscFunctionBeginUser;
119:   PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));

121:   PetscCall(PetscOptionsHasName(NULL, NULL, "-keep_nonzero_pattern", &keepnonzeropattern));
122:   if (keepnonzeropattern) PetscCall(MatSetOption(B, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE));

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: }

130: PetscErrorCode TestMatZeroRows_with_no_allocation(Mat A, IS is, PetscScalar diag)
131: {
132:   Mat B;

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));

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: }

146: /*TEST

148:    test:
149:       nsize: 2
150:       filter: grep -v " MPI process"

152:    test:
153:       suffix: 2
154:       nsize: 3
155:       args: -mat_type mpibaij -mat_block_size 3
156:       filter: grep -v " MPI process"

158:    test:
159:       suffix: 3
160:       nsize: 3
161:       args: -mat_type mpiaij -keep_nonzero_pattern
162:       filter: grep -v " MPI process"

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"

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"

177: TEST*/