Actual source code: ex167.c


  2: static char help[] = "Extract submatrices using unsorted indices. For SEQSBAIJ either sort both rows and columns, or sort none.\n\n";
  3: /*
  4:    Take a 4x4 grid and form a 5-point stencil graph Laplacian over it.
  5:    Partition the grid into two subdomains by splitting into two in the j-direction (slowest varying).
  6:    Impose an overlap of 1 and order the subdomains with the j-direction varying fastest.
  7:    Extract the subdomain submatrices, one per rank.
  8: */
  9: /* Results:
 10:     Sequential:
 11:     - seqaij:   will error out, if rows or columns are unsorted
 12:     - seqbaij:  will extract submatrices correctly even for unsorted row or column indices
 13:     - seqsbaij: will extract submatrices correctly even for unsorted row and column indices (both must be sorted or not);
 14:                 CANNOT automatically report inversions, because MatGetRow is not available.
 15:     MPI:
 16:     - mpiaij:   will error out, if columns are unsorted
 17:     - mpibaij:  will error out, if columns are unsorted.
 18:     - mpisbaij: will error out, if columns are unsorted; even with unsorted rows will produce correct submatrices;
 19:                 CANNOT automatically report inversions, because MatGetRow is not available.
 20: */

 22: #include <petscmat.h>
 23: #include <petscis.h>

 25: int main(int argc, char **args)
 26: {
 27:   Mat             A, *S;
 28:   IS              rowis[2], colis[2];
 29:   PetscInt        n, N, i, j, k, l, nsub, Jlow[2] = {0, 1}, *jlow, Jhigh[2] = {3, 4}, *jhigh, row, col, *subindices, ncols;
 30:   const PetscInt *cols;
 31:   PetscScalar     v;
 32:   PetscMPIInt     rank, size, p, inversions, total_inversions;
 33:   PetscBool       sort_rows, sort_cols, show_inversions;

 36:   PetscInitialize(&argc, &args, (char *)0, help);
 37:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 38:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

 41:   MatCreate(PETSC_COMM_WORLD, &A);
 42:   if (size > 1) {
 43:     n = 8;
 44:     N = 16;
 45:   } else {
 46:     n = 16;
 47:     N = 16;
 48:   }
 49:   MatSetSizes(A, n, n, N, N);
 50:   MatSetFromOptions(A);
 51:   MatSetUp(A);

 53:   /* Don't care if the entries are set multiple times by different procs. */
 54:   for (i = 0; i < 4; ++i) {
 55:     for (j = 0; j < 4; ++j) {
 56:       row = j * 4 + i;
 57:       v   = -1.0;
 58:       if (i > 0) {
 59:         col = row - 1;
 60:         MatSetValues(A, 1, &row, 1, &col, &v, INSERT_VALUES);
 61:       }
 62:       if (i < 3) {
 63:         col = row + 1;
 64:         MatSetValues(A, 1, &row, 1, &col, &v, INSERT_VALUES);
 65:       }
 66:       if (j > 0) {
 67:         col = row - 4;
 68:         MatSetValues(A, 1, &row, 1, &col, &v, INSERT_VALUES);
 69:       }
 70:       if (j < 3) {
 71:         col = row + 4;
 72:         MatSetValues(A, 1, &row, 1, &col, &v, INSERT_VALUES);
 73:       }
 74:       v = 4.0;
 75:       MatSetValues(A, 1, &row, 1, &row, &v, INSERT_VALUES);
 76:     }
 77:   }
 78:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 79:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 80:   PetscPrintf(PETSC_COMM_WORLD, "Original matrix\n");
 81:   MatView(A, PETSC_VIEWER_STDOUT_WORLD);

 83:   if (size > 1) {
 84:     nsub = 1; /* one subdomain per rank */
 85:   } else {
 86:     nsub = 2; /* both subdomains on rank 0 */
 87:   }
 88:   if (rank) {
 89:     jlow  = Jlow + 1;
 90:     jhigh = Jhigh + 1;
 91:   } else {
 92:     jlow  = Jlow;
 93:     jhigh = Jhigh;
 94:   }
 95:   sort_rows = PETSC_FALSE;
 96:   PetscOptionsGetBool(NULL, NULL, "-sort_rows", &sort_rows, NULL);
 97:   sort_cols = PETSC_FALSE;
 98:   PetscOptionsGetBool(NULL, NULL, "-sort_cols", &sort_cols, NULL);
 99:   for (l = 0; l < nsub; ++l) {
100:     PetscMalloc1(12, &subindices);
101:     k = 0;
102:     for (i = 0; i < 4; ++i) {
103:       for (j = jlow[l]; j < jhigh[l]; ++j) {
104:         subindices[k] = j * 4 + i;
105:         k++;
106:       }
107:     }
108:     ISCreateGeneral(PETSC_COMM_SELF, 12, subindices, PETSC_OWN_POINTER, rowis + l);
109:     if ((sort_rows && !sort_cols) || (!sort_rows && sort_cols)) {
110:       ISDuplicate(rowis[l], colis + l);
111:     } else {
112:       PetscObjectReference((PetscObject)rowis[l]);
113:       colis[l] = rowis[l];
114:     }
115:     if (sort_rows) ISSort(rowis[l]);
116:     if (sort_cols) ISSort(colis[l]);
117:   }

119:   MatCreateSubMatrices(A, nsub, rowis, colis, MAT_INITIAL_MATRIX, &S);

121:   show_inversions = PETSC_FALSE;

123:   PetscOptionsGetBool(NULL, NULL, "-show_inversions", &show_inversions, NULL);

125:   inversions = 0;
126:   for (p = 0; p < size; ++p) {
127:     if (p == rank) {
128:       PetscPrintf(PETSC_COMM_SELF, "[%" PetscInt_FMT ":%" PetscInt_FMT "]: Number of subdomains: %" PetscInt_FMT ":\n", rank, size, nsub);
129:       for (l = 0; l < nsub; ++l) {
130:         PetscInt i0, i1;
131:         PetscPrintf(PETSC_COMM_SELF, "[%" PetscInt_FMT ":%" PetscInt_FMT "]: Subdomain row IS %" PetscInt_FMT ":\n", rank, size, l);
132:         ISView(rowis[l], PETSC_VIEWER_STDOUT_SELF);
133:         PetscPrintf(PETSC_COMM_SELF, "[%" PetscInt_FMT ":%" PetscInt_FMT "]: Subdomain col IS %" PetscInt_FMT ":\n", rank, size, l);
134:         ISView(colis[l], PETSC_VIEWER_STDOUT_SELF);
135:         PetscPrintf(PETSC_COMM_SELF, "[%" PetscInt_FMT ":%" PetscInt_FMT "]: Submatrix %" PetscInt_FMT ":\n", rank, size, l);
136:         MatView(S[l], PETSC_VIEWER_STDOUT_SELF);
137:         if (show_inversions) {
138:           MatGetOwnershipRange(S[l], &i0, &i1);
139:           for (i = i0; i < i1; ++i) {
140:             MatGetRow(S[l], i, &ncols, &cols, NULL);
141:             for (j = 1; j < ncols; ++j) {
142:               if (cols[j] < cols[j - 1]) {
143:                 PetscPrintf(PETSC_COMM_SELF, "***Inversion in row %" PetscInt_FMT ": col[%" PetscInt_FMT "] = %" PetscInt_FMT " < %" PetscInt_FMT " = col[%" PetscInt_FMT "]\n", i, j, cols[j], cols[j - 1], j - 1);
144:                 inversions++;
145:               }
146:             }
147:             MatRestoreRow(S[l], i, &ncols, &cols, NULL);
148:           }
149:         }
150:       }
151:     }
152:     MPI_Barrier(PETSC_COMM_WORLD);
153:   }
154:   if (show_inversions) {
155:     MPI_Reduce(&inversions, &total_inversions, 1, MPIU_INT, MPI_SUM, 0, PETSC_COMM_WORLD);
156:     PetscPrintf(PETSC_COMM_WORLD, "*Total inversions: %" PetscInt_FMT "\n", total_inversions);
157:   }
158:   MatDestroy(&A);

160:   for (l = 0; l < nsub; ++l) {
161:     ISDestroy(&(rowis[l]));
162:     ISDestroy(&(colis[l]));
163:   }
164:   MatDestroySubMatrices(nsub, &S);
165:   PetscFinalize();
166:   return 0;
167: }