Actual source code: ex9.c
1: static char help[] = "Tests MPI parallel matrix creation. Test MatCreateRedundantMatrix() \n\n";
3: #include <petscmat.h>
5: int main(int argc, char **args)
6: {
7: Mat C, Credundant;
8: MatInfo info;
9: PetscMPIInt rank, size, subsize;
10: PetscInt i, j, m = 3, n = 2, low, high, iglobal;
11: PetscInt Ii, J, ldim, nsubcomms;
12: PetscBool flg_info, flg_mat;
13: PetscScalar v, one = 1.0;
14: Vec x, y;
15: MPI_Comm subcomm;
17: PetscFunctionBeginUser;
18: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
19: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
20: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
21: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
22: n = 2 * size;
24: PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
25: PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n));
26: PetscCall(MatSetFromOptions(C));
27: PetscCall(MatSetUp(C));
29: /* Create the matrix for the five point stencil, YET AGAIN */
30: for (i = 0; i < m; i++) {
31: for (j = 2 * rank; j < 2 * rank + 2; j++) {
32: v = -1.0;
33: Ii = j + n * i;
34: if (i > 0) {
35: J = Ii - n;
36: PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
37: }
38: if (i < m - 1) {
39: J = Ii + n;
40: PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
41: }
42: if (j > 0) {
43: J = Ii - 1;
44: PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
45: }
46: if (j < n - 1) {
47: J = Ii + 1;
48: PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
49: }
50: v = 4.0;
51: PetscCall(MatSetValues(C, 1, &Ii, 1, &Ii, &v, INSERT_VALUES));
52: }
53: }
55: /* Add extra elements (to illustrate variants of MatGetInfo) */
56: Ii = n;
57: J = n - 2;
58: v = 100.0;
59: PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
60: Ii = n - 2;
61: J = n;
62: v = 100.0;
63: PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
65: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
66: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
68: /* Form vectors */
69: PetscCall(MatCreateVecs(C, &x, &y));
70: PetscCall(VecGetLocalSize(x, &ldim));
71: PetscCall(VecGetOwnershipRange(x, &low, &high));
72: for (i = 0; i < ldim; i++) {
73: iglobal = i + low;
74: v = one * ((PetscReal)i) + 100.0 * rank;
75: PetscCall(VecSetValues(x, 1, &iglobal, &v, INSERT_VALUES));
76: }
77: PetscCall(VecAssemblyBegin(x));
78: PetscCall(VecAssemblyEnd(x));
80: PetscCall(MatMult(C, x, y));
82: PetscCall(PetscOptionsHasName(NULL, NULL, "-view_info", &flg_info));
83: if (flg_info) {
84: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO));
85: PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
87: PetscCall(MatGetInfo(C, MAT_GLOBAL_SUM, &info));
88: PetscCall(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));
89: PetscCall(MatGetInfo(C, MAT_GLOBAL_MAX, &info));
90: PetscCall(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));
91: }
93: PetscCall(PetscOptionsHasName(NULL, NULL, "-view_mat", &flg_mat));
94: if (flg_mat) PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
96: /* Test MatCreateRedundantMatrix() */
97: nsubcomms = size;
98: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nsubcomms", &nsubcomms, NULL));
99: PetscCall(MatCreateRedundantMatrix(C, nsubcomms, MPI_COMM_NULL, MAT_INITIAL_MATRIX, &Credundant));
100: PetscCall(MatCreateRedundantMatrix(C, nsubcomms, MPI_COMM_NULL, MAT_REUSE_MATRIX, &Credundant));
102: PetscCall(PetscObjectGetComm((PetscObject)Credundant, &subcomm));
103: PetscCallMPI(MPI_Comm_size(subcomm, &subsize));
105: if (subsize == 2 && flg_mat) {
106: PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(subcomm), "\n[%d] Credundant:\n", rank));
107: PetscCall(MatView(Credundant, PETSC_VIEWER_STDOUT_(subcomm)));
108: }
109: PetscCall(MatDestroy(&Credundant));
111: /* Test MatCreateRedundantMatrix() with user-provided subcomm */
112: {
113: PetscSubcomm psubcomm;
115: PetscCall(PetscSubcommCreate(PETSC_COMM_WORLD, &psubcomm));
116: PetscCall(PetscSubcommSetNumber(psubcomm, nsubcomms));
117: PetscCall(PetscSubcommSetType(psubcomm, PETSC_SUBCOMM_CONTIGUOUS));
118: /* enable runtime switch of psubcomm type, e.g., '-psubcomm_type interlaced */
119: PetscCall(PetscSubcommSetFromOptions(psubcomm));
121: PetscCall(MatCreateRedundantMatrix(C, nsubcomms, PetscSubcommChild(psubcomm), MAT_INITIAL_MATRIX, &Credundant));
122: PetscCall(MatCreateRedundantMatrix(C, nsubcomms, PetscSubcommChild(psubcomm), MAT_REUSE_MATRIX, &Credundant));
124: PetscCall(PetscSubcommDestroy(&psubcomm));
125: PetscCall(MatDestroy(&Credundant));
126: }
128: PetscCall(VecDestroy(&x));
129: PetscCall(VecDestroy(&y));
130: PetscCall(MatDestroy(&C));
131: PetscCall(PetscFinalize());
132: return 0;
133: }
135: /*TEST
137: test:
138: nsize: 3
139: args: -view_info
141: test:
142: suffix: 2
143: nsize: 3
144: args: -nsubcomms 2 -view_mat -psubcomm_type interlaced
146: test:
147: suffix: 3
148: nsize: 3
149: args: -nsubcomms 2 -view_mat -psubcomm_type contiguous
151: test:
152: suffix: 3_baij
153: nsize: 3
154: args: -mat_type baij -nsubcomms 2 -view_mat
156: test:
157: suffix: 3_sbaij
158: nsize: 3
159: args: -mat_type sbaij -nsubcomms 2 -view_mat
161: test:
162: suffix: 3_dense
163: nsize: 3
164: args: -mat_type dense -nsubcomms 2 -view_mat
166: test:
167: suffix: 4_baij
168: nsize: 3
169: args: -mat_type baij -nsubcomms 2 -view_mat -psubcomm_type interlaced
171: test:
172: suffix: 4_sbaij
173: nsize: 3
174: args: -mat_type sbaij -nsubcomms 2 -view_mat -psubcomm_type interlaced
176: test:
177: suffix: 4_dense
178: nsize: 3
179: args: -mat_type dense -nsubcomms 2 -view_mat -psubcomm_type interlaced
181: TEST*/