Actual source code: ex9.c


  2: static char help[] = "Tests MPI parallel matrix creation. Test MatCreateRedundantMatrix() \n\n";

  4: #include <petscmat.h>

  6: int main(int argc, char **args)
  7: {
  8:   Mat         C, Credundant;
  9:   MatInfo     info;
 10:   PetscMPIInt rank, size, subsize;
 11:   PetscInt    i, j, m = 3, n = 2, low, high, iglobal;
 12:   PetscInt    Ii, J, ldim, nsubcomms;
 13:   PetscBool   flg_info, flg_mat;
 14:   PetscScalar v, one = 1.0;
 15:   Vec         x, y;
 16:   MPI_Comm    subcomm;

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

 25:   MatCreate(PETSC_COMM_WORLD, &C);
 26:   MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n);
 27:   MatSetFromOptions(C);
 28:   MatSetUp(C);

 30:   /* Create the matrix for the five point stencil, YET AGAIN */
 31:   for (i = 0; i < m; i++) {
 32:     for (j = 2 * rank; j < 2 * rank + 2; j++) {
 33:       v  = -1.0;
 34:       Ii = j + n * i;
 35:       if (i > 0) {
 36:         J = Ii - n;
 37:         MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 38:       }
 39:       if (i < m - 1) {
 40:         J = Ii + n;
 41:         MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 42:       }
 43:       if (j > 0) {
 44:         J = Ii - 1;
 45:         MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 46:       }
 47:       if (j < n - 1) {
 48:         J = Ii + 1;
 49:         MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 50:       }
 51:       v = 4.0;
 52:       MatSetValues(C, 1, &Ii, 1, &Ii, &v, INSERT_VALUES);
 53:     }
 54:   }

 56:   /* Add extra elements (to illustrate variants of MatGetInfo) */
 57:   Ii = n;
 58:   J  = n - 2;
 59:   v  = 100.0;
 60:   MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 61:   Ii = n - 2;
 62:   J  = n;
 63:   v  = 100.0;
 64:   MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);

 66:   MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
 67:   MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);

 69:   /* Form vectors */
 70:   MatCreateVecs(C, &x, &y);
 71:   VecGetLocalSize(x, &ldim);
 72:   VecGetOwnershipRange(x, &low, &high);
 73:   for (i = 0; i < ldim; i++) {
 74:     iglobal = i + low;
 75:     v       = one * ((PetscReal)i) + 100.0 * rank;
 76:     VecSetValues(x, 1, &iglobal, &v, INSERT_VALUES);
 77:   }
 78:   VecAssemblyBegin(x);
 79:   VecAssemblyEnd(x);

 81:   MatMult(C, x, y);

 83:   PetscOptionsHasName(NULL, NULL, "-view_info", &flg_info);
 84:   if (flg_info) {
 85:     PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO);
 86:     MatView(C, PETSC_VIEWER_STDOUT_WORLD);

 88:     MatGetInfo(C, MAT_GLOBAL_SUM, &info);
 89:     PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "matrix information (global sums):\nnonzeros = %" PetscInt_FMT ", allocated nonzeros = %" PetscInt_FMT "\n", (PetscInt)info.nz_used, (PetscInt)info.nz_allocated);
 90:     MatGetInfo(C, MAT_GLOBAL_MAX, &info);
 91:     PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "matrix information (global max):\nnonzeros = %" PetscInt_FMT ", allocated nonzeros = %" PetscInt_FMT "\n", (PetscInt)info.nz_used, (PetscInt)info.nz_allocated);
 92:   }

 94:   PetscOptionsHasName(NULL, NULL, "-view_mat", &flg_mat);
 95:   if (flg_mat) MatView(C, PETSC_VIEWER_STDOUT_WORLD);

 97:   /* Test MatCreateRedundantMatrix() */
 98:   nsubcomms = size;
 99:   PetscOptionsGetInt(NULL, NULL, "-nsubcomms", &nsubcomms, NULL);
100:   MatCreateRedundantMatrix(C, nsubcomms, MPI_COMM_NULL, MAT_INITIAL_MATRIX, &Credundant);
101:   MatCreateRedundantMatrix(C, nsubcomms, MPI_COMM_NULL, MAT_REUSE_MATRIX, &Credundant);

103:   PetscObjectGetComm((PetscObject)Credundant, &subcomm);
104:   MPI_Comm_size(subcomm, &subsize);

106:   if (subsize == 2 && flg_mat) {
107:     PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(subcomm), "\n[%d] Credundant:\n", rank);
108:     MatView(Credundant, PETSC_VIEWER_STDOUT_(subcomm));
109:   }
110:   MatDestroy(&Credundant);

112:   /* Test MatCreateRedundantMatrix() with user-provided subcomm */
113:   {
114:     PetscSubcomm psubcomm;

116:     PetscSubcommCreate(PETSC_COMM_WORLD, &psubcomm);
117:     PetscSubcommSetNumber(psubcomm, nsubcomms);
118:     PetscSubcommSetType(psubcomm, PETSC_SUBCOMM_CONTIGUOUS);
119:     /* enable runtime switch of psubcomm type, e.g., '-psubcomm_type interlaced */
120:     PetscSubcommSetFromOptions(psubcomm);

122:     MatCreateRedundantMatrix(C, nsubcomms, PetscSubcommChild(psubcomm), MAT_INITIAL_MATRIX, &Credundant);
123:     MatCreateRedundantMatrix(C, nsubcomms, PetscSubcommChild(psubcomm), MAT_REUSE_MATRIX, &Credundant);

125:     PetscSubcommDestroy(&psubcomm);
126:     MatDestroy(&Credundant);
127:   }

129:   VecDestroy(&x);
130:   VecDestroy(&y);
131:   MatDestroy(&C);
132:   PetscFinalize();
133:   return 0;
134: }

136: /*TEST

138:    test:
139:       nsize: 3
140:       args: -view_info

142:    test:
143:       suffix: 2
144:       nsize: 3
145:       args: -nsubcomms 2 -view_mat -psubcomm_type interlaced

147:    test:
148:       suffix: 3
149:       nsize: 3
150:       args: -nsubcomms 2 -view_mat -psubcomm_type contiguous

152:    test:
153:       suffix: 3_baij
154:       nsize: 3
155:       args: -mat_type baij -nsubcomms 2 -view_mat

157:    test:
158:       suffix: 3_sbaij
159:       nsize: 3
160:       args: -mat_type sbaij -nsubcomms 2 -view_mat

162:    test:
163:       suffix: 3_dense
164:       nsize: 3
165:       args: -mat_type dense -nsubcomms 2 -view_mat

167:    test:
168:       suffix: 4_baij
169:       nsize: 3
170:       args: -mat_type baij -nsubcomms 2 -view_mat -psubcomm_type interlaced

172:    test:
173:       suffix: 4_sbaij
174:       nsize: 3
175:       args: -mat_type sbaij -nsubcomms 2 -view_mat -psubcomm_type interlaced

177:    test:
178:       suffix: 4_dense
179:       nsize: 3
180:       args: -mat_type dense -nsubcomms 2 -view_mat -psubcomm_type interlaced

182: TEST*/