Actual source code: ex37.c


  2: static char help[] = "Test MatGetMultiProcBlock() and MatCreateRedundantMatrix() \n\
  3: Reads a PETSc matrix and vector from a file and solves a linear system.\n\n";
  4: /*
  5:   Example:
  6:   mpiexec -n 4 ./ex37 -f <input_file> -nsubcomm 2 -psubcomm_view -subcomm_type <1 or 2>
  7: */

  9: #include <petscksp.h>
 10: #include <petscsys.h>

 12: int main(int argc, char **args)
 13: {
 14:   KSP          subksp;
 15:   Mat          A, subA;
 16:   Vec          x, b, u, subb, subx, subu;
 17:   PetscViewer  fd;
 18:   char         file[PETSC_MAX_PATH_LEN];
 19:   PetscBool    flg;
 20:   PetscInt     i, m, n, its;
 21:   PetscReal    norm;
 22:   PetscMPIInt  rank, size;
 23:   MPI_Comm     comm, subcomm;
 24:   PetscSubcomm psubcomm;
 25:   PetscInt     nsubcomm = 1, id;
 26:   PetscScalar *barray, *xarray, *uarray, *array, one = 1.0;
 27:   PetscInt     type = 1;

 30:   PetscInitialize(&argc, &args, (char *)0, help);
 31:   /* Load the matrix */
 32:   PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg);
 34:   PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd);

 36:   /* Load the matrix; then destroy the viewer.*/
 37:   MatCreate(PETSC_COMM_WORLD, &A);
 38:   MatLoad(A, fd);
 39:   PetscViewerDestroy(&fd);

 41:   PetscObjectGetComm((PetscObject)A, &comm);
 42:   MPI_Comm_size(comm, &size);
 43:   MPI_Comm_rank(comm, &rank);

 45:   /* Create rhs vector b */
 46:   MatGetLocalSize(A, &m, NULL);
 47:   VecCreate(PETSC_COMM_WORLD, &b);
 48:   VecSetSizes(b, m, PETSC_DECIDE);
 49:   VecSetFromOptions(b);
 50:   VecSet(b, one);

 52:   VecDuplicate(b, &x);
 53:   VecDuplicate(b, &u);
 54:   VecSet(x, 0.0);

 56:   /* Test MatGetMultiProcBlock() */
 57:   PetscOptionsGetInt(NULL, NULL, "-nsubcomm", &nsubcomm, NULL);
 58:   PetscOptionsGetInt(NULL, NULL, "-subcomm_type", &type, NULL);

 60:   PetscSubcommCreate(comm, &psubcomm);
 61:   PetscSubcommSetNumber(psubcomm, nsubcomm);
 62:   if (type == PETSC_SUBCOMM_GENERAL) { /* user provides color, subrank and duprank */
 63:     PetscMPIInt color, subrank, duprank, subsize;
 64:     duprank = size - 1 - rank;
 65:     subsize = size / nsubcomm;
 67:     color   = duprank / subsize;
 68:     subrank = duprank - color * subsize;
 69:     PetscSubcommSetTypeGeneral(psubcomm, color, subrank);
 70:   } else if (type == PETSC_SUBCOMM_CONTIGUOUS) {
 71:     PetscSubcommSetType(psubcomm, PETSC_SUBCOMM_CONTIGUOUS);
 72:   } else if (type == PETSC_SUBCOMM_INTERLACED) {
 73:     PetscSubcommSetType(psubcomm, PETSC_SUBCOMM_INTERLACED);
 74:   } else SETERRQ(psubcomm->parent, PETSC_ERR_SUP, "PetscSubcommType %" PetscInt_FMT " is not supported yet", type);
 75:   PetscSubcommSetFromOptions(psubcomm);
 76:   subcomm = PetscSubcommChild(psubcomm);

 78:   /* Test MatCreateRedundantMatrix() */
 79:   if (size > 1) {
 80:     PetscMPIInt subrank = -1, color = -1;
 81:     MPI_Comm    dcomm;

 83:     if (rank == 0) {
 84:       color   = 0;
 85:       subrank = 0;
 86:     } else if (rank == 1) {
 87:       color   = 0;
 88:       subrank = 1;
 89:     } else {
 90:       color   = 1;
 91:       subrank = 0;
 92:     }

 94:     PetscCommDuplicate(PETSC_COMM_WORLD, &dcomm, NULL);
 95:     MPI_Comm_split(dcomm, color, subrank, &subcomm);

 97:     MatCreate(subcomm, &subA);
 98:     MatSetSizes(subA, PETSC_DECIDE, PETSC_DECIDE, 10, 10);
 99:     MatSetFromOptions(subA);
100:     MatSetUp(subA);
101:     MatAssemblyBegin(subA, MAT_FINAL_ASSEMBLY);
102:     MatAssemblyEnd(subA, MAT_FINAL_ASSEMBLY);
103:     MatZeroEntries(subA);

105:     /* Test MatMult() */
106:     MatCreateVecs(subA, &subx, &subb);
107:     VecSet(subx, 1.0);
108:     MatMult(subA, subx, subb);

110:     VecDestroy(&subx);
111:     VecDestroy(&subb);
112:     MatDestroy(&subA);
113:     PetscCommDestroy(&dcomm);
114:   }

116:   /* Create subA */
117:   MatGetMultiProcBlock(A, subcomm, MAT_INITIAL_MATRIX, &subA);
118:   MatGetMultiProcBlock(A, subcomm, MAT_REUSE_MATRIX, &subA);

120:   /* Create sub vectors without arrays. Place b's and x's local arrays into subb and subx */
121:   MatGetLocalSize(subA, &m, &n);
122:   VecCreateMPIWithArray(subcomm, 1, m, PETSC_DECIDE, NULL, &subb);
123:   VecCreateMPIWithArray(subcomm, 1, n, PETSC_DECIDE, NULL, &subx);
124:   VecCreateMPIWithArray(subcomm, 1, n, PETSC_DECIDE, NULL, &subu);

126:   VecGetArray(b, &barray);
127:   VecGetArray(x, &xarray);
128:   VecGetArray(u, &uarray);
129:   VecPlaceArray(subb, barray);
130:   VecPlaceArray(subx, xarray);
131:   VecPlaceArray(subu, uarray);

133:   /* Create linear solvers associated with subA */
134:   KSPCreate(subcomm, &subksp);
135:   KSPSetOperators(subksp, subA, subA);
136:   KSPSetFromOptions(subksp);

138:   /* Solve sub systems */
139:   KSPSolve(subksp, subb, subx);
140:   KSPGetIterationNumber(subksp, &its);

142:   /* check residual */
143:   MatMult(subA, subx, subu);
144:   VecAXPY(subu, -1.0, subb);
145:   VecNorm(u, NORM_2, &norm);
146:   if (norm > 1.e-4 && rank == 0) {
147:     PetscPrintf(PETSC_COMM_WORLD, "[%d]  Number of iterations = %3" PetscInt_FMT "\n", rank, its);
148:     PetscPrintf(PETSC_COMM_WORLD, "Error: Residual norm of each block |subb - subA*subx |= %g\n", (double)norm);
149:   }
150:   VecResetArray(subb);
151:   VecResetArray(subx);
152:   VecResetArray(subu);

154:   PetscOptionsGetInt(NULL, NULL, "-subvec_view", &id, &flg);
155:   if (flg && rank == id) {
156:     PetscPrintf(PETSC_COMM_SELF, "[%d] subb:\n", rank);
157:     VecGetArray(subb, &array);
158:     for (i = 0; i < m; i++) PetscPrintf(PETSC_COMM_SELF, "%g\n", (double)PetscRealPart(array[i]));
159:     VecRestoreArray(subb, &array);
160:     PetscPrintf(PETSC_COMM_SELF, "[%d] subx:\n", rank);
161:     VecGetArray(subx, &array);
162:     for (i = 0; i < m; i++) PetscPrintf(PETSC_COMM_SELF, "%g\n", (double)PetscRealPart(array[i]));
163:     VecRestoreArray(subx, &array);
164:   }

166:   VecRestoreArray(x, &xarray);
167:   VecRestoreArray(b, &barray);
168:   VecRestoreArray(u, &uarray);
169:   MatDestroy(&subA);
170:   VecDestroy(&subb);
171:   VecDestroy(&subx);
172:   VecDestroy(&subu);
173:   KSPDestroy(&subksp);
174:   PetscSubcommDestroy(&psubcomm);
175:   if (size > 1) MPI_Comm_free(&subcomm);
176:   MatDestroy(&A);
177:   VecDestroy(&b);
178:   VecDestroy(&u);
179:   VecDestroy(&x);

181:   PetscFinalize();
182:   return 0;
183: }

185: /*TEST

187:     test:
188:       args: -f ${DATAFILESPATH}/matrices/small -nsubcomm 1
189:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
190:       output_file: output/ex37.out

192:     test:
193:       suffix: 2
194:       args: -f ${DATAFILESPATH}/matrices/small -nsubcomm 2
195:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
196:       nsize: 4
197:       output_file: output/ex37.out

199:     test:
200:       suffix: mumps
201:       args: -f ${DATAFILESPATH}/matrices/small -nsubcomm 2 -pc_factor_mat_solver_type mumps -pc_type lu
202:       requires: datafilespath  mumps !complex double !defined(PETSC_USE_64BIT_INDICES)
203:       nsize: 4
204:       output_file: output/ex37.out

206:     test:
207:       suffix: 3
208:       nsize: 4
209:       args: -f ${DATAFILESPATH}/matrices/small -nsubcomm 2 -subcomm_type 0
210:       requires: datafilespath  !complex double !defined(PETSC_USE_64BIT_INDICES)
211:       output_file: output/ex37.out

213:     test:
214:       suffix: 4
215:       nsize: 4
216:       args: -f ${DATAFILESPATH}/matrices/small -nsubcomm 2 -subcomm_type 1
217:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
218:       output_file: output/ex37.out

220:     test:
221:       suffix: 5
222:       nsize: 4
223:       args: -f ${DATAFILESPATH}/matrices/small -nsubcomm 2 -subcomm_type 2
224:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
225:       output_file: output/ex37.out

227: TEST*/