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;

 19:   PetscInitialize(&argc, &args, (char *)0, help);
 20:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 21:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 22:   n = 2 * size;

 24:   /* create A Square matrix for the five point stencil,YET AGAIN*/
 25:   MatCreate(PETSC_COMM_WORLD, &A);
 26:   MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n);
 27:   MatSetFromOptions(A);
 28:   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:         MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 36:       }
 37:       if (i < m - 1) {
 38:         J = Ii + n;
 39:         MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 40:       }
 41:       if (j > 0) {
 42:         J = Ii - 1;
 43:         MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 44:       }
 45:       if (j < n - 1) {
 46:         J = Ii + 1;
 47:         MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 48:       }
 49:       v = 4.0;
 50:       MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES);
 51:     }
 52:   }
 53:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 54:   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:   ISCreateStride(PETSC_COMM_SELF, m, Imax, 1, &is);

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

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

 67:   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:   MatCreate(PETSC_COMM_WORLD, &A);
 73:   MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n, m * (n + size));
 74:   MatSetFromOptions(A);
 75:   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:         MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 83:       }
 84:       if (i < m - 1) {
 85:         J = Ii + n;
 86:         MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 87:       }
 88:       if (j > 0) {
 89:         J = Ii - 1;
 90:         MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 91:       }
 92:       if (j < n + size - 1) {
 93:         J = Ii + 1;
 94:         MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 95:       }
 96:       v = 4.0;
 97:       MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES);
 98:     }
 99:   }
100:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
101:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

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

106:   MatDestroy(&A);
107:   ISDestroy(&is);
108:   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:   MatDuplicate(A, MAT_COPY_VALUES, &B);

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

123:   MatZeroRowsIS(B, is, diag, 0, 0);
124:   MatView(B, PETSC_VIEWER_STDOUT_WORLD);
125:   MatDestroy(&B);
126:   return 0;
127: }

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

133:   /* Now copy A into B, and test it with MatZeroRows() */
134:   MatDuplicate(A, MAT_COPY_VALUES, &B);
135:   /* Set this flag after assembly. This way, it affects only MatZeroRows() */
136:   MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);

138:   MatZeroRowsIS(B, is, diag, 0, 0);
139:   MatView(B, PETSC_VIEWER_STDOUT_WORLD);
140:   MatDestroy(&B);
141:   return 0;
142: }

144: /*TEST

146:    test:
147:       nsize: 2
148:       filter: grep -v " MPI process"

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

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

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"

168: TEST*/