Actual source code: ex81.c
1: #include <petsc.h>
3: static char help[] = "Solves a linear system with a MatNest and nested fields.\n\n";
5: #define Q 5 /* everything is hardwired for a 5x5 MatNest for now */
7: int main(int argc, char **args)
8: {
9: KSP ksp;
10: PC pc;
11: Mat array[Q * Q], A, a;
12: Vec b, x, sub;
13: IS rows[Q];
14: PetscInt i, j, *outer, M, N, found = Q;
15: PetscMPIInt size;
16: PetscBool flg;
17: PetscRandom rctx;
19: PetscFunctionBeginUser;
20: PetscCall(PetscInitialize(&argc, &args, NULL, help));
21: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
22: PetscCall(PetscMalloc1(found, &outer));
23: PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-outer_fieldsplit_sizes", outer, &found, &flg));
24: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
25: if (flg) {
26: PetscCheck(found != 1, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must supply more than one field");
27: j = 0;
28: for (i = 0; i < found; ++i) j += outer[i];
29: PetscCheck(j == Q, PETSC_COMM_WORLD, PETSC_ERR_USER, "Sum of outer fieldsplit sizes (%" PetscInt_FMT ") greater than number of blocks in MatNest (%d)", j, Q);
30: }
31: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
32: size = PetscMax(3, size);
33: for (i = 0; i < Q * Q; ++i) array[i] = NULL;
34: for (i = 0; i < Q; ++i) {
35: if (i == 0) {
36: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, size, size, 1, NULL, 0, NULL, array + (Q + 1) * i));
37: } else if (i == 1 || i == 3) {
38: PetscCall(MatCreateSBAIJ(PETSC_COMM_WORLD, 2, PETSC_DECIDE, PETSC_DECIDE, size, size, 1, NULL, 0, NULL, array + (Q + 1) * i));
39: } else if (i == 2 || i == 4) {
40: PetscCall(MatCreateBAIJ(PETSC_COMM_WORLD, 2, PETSC_DECIDE, PETSC_DECIDE, size, size, 1, NULL, 0, NULL, array + (Q + 1) * i));
41: }
42: PetscCall(MatAssemblyBegin(array[(Q + 1) * i], MAT_FINAL_ASSEMBLY));
43: PetscCall(MatAssemblyEnd(array[(Q + 1) * i], MAT_FINAL_ASSEMBLY));
44: PetscCall(MatShift(array[(Q + 1) * i], 100 + i + 1));
45: if (i == 3) {
46: PetscCall(MatDuplicate(array[(Q + 1) * i], MAT_COPY_VALUES, &a));
47: PetscCall(MatDestroy(array + (Q + 1) * i));
48: PetscCall(MatCreateHermitianTranspose(a, array + (Q + 1) * i));
49: PetscCall(MatDestroy(&a));
50: }
51: size *= 2;
52: }
53: PetscCall(MatGetSize(array[0], &M, NULL));
54: for (i = 2; i < Q; ++i) {
55: PetscCall(MatGetSize(array[(Q + 1) * i], NULL, &N));
56: if (i != Q - 1) {
57: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, i == 3 ? N : M, i == 3 ? M : N, 0, NULL, 0, NULL, array + i));
58: } else {
59: PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, M, N, NULL, array + i));
60: }
61: PetscCall(MatAssemblyBegin(array[i], MAT_FINAL_ASSEMBLY));
62: PetscCall(MatAssemblyEnd(array[i], MAT_FINAL_ASSEMBLY));
63: PetscCall(MatSetRandom(array[i], rctx));
64: if (i == 3) {
65: PetscCall(MatDuplicate(array[i], MAT_COPY_VALUES, &a));
66: PetscCall(MatDestroy(array + i));
67: PetscCall(MatCreateHermitianTranspose(a, array + i));
68: PetscCall(MatDestroy(&a));
69: }
70: }
71: PetscCall(MatGetSize(array[0], NULL, &N));
72: for (i = 2; i < Q; i += 2) {
73: PetscCall(MatGetSize(array[(Q + 1) * i], &M, NULL));
74: if (i != Q - 1) {
75: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, M, N, 2, NULL, 2, NULL, array + Q * i));
76: } else {
77: PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, N, M, NULL, array + Q * i));
78: }
79: PetscCall(MatAssemblyBegin(array[Q * i], MAT_FINAL_ASSEMBLY));
80: PetscCall(MatAssemblyEnd(array[Q * i], MAT_FINAL_ASSEMBLY));
81: PetscCall(MatSetRandom(array[Q * i], rctx));
82: if (i == Q - 1) {
83: PetscCall(MatDuplicate(array[Q * i], MAT_COPY_VALUES, &a));
84: PetscCall(MatDestroy(array + Q * i));
85: PetscCall(MatCreateHermitianTranspose(a, array + Q * i));
86: PetscCall(MatDestroy(&a));
87: }
88: }
89: PetscCall(MatGetSize(array[(Q + 1) * 3], &M, NULL));
90: for (i = 1; i < 3; ++i) {
91: PetscCall(MatGetSize(array[(Q + 1) * i], NULL, &N));
92: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, M, N, 2, NULL, 2, NULL, array + Q * 3 + i));
93: PetscCall(MatAssemblyBegin(array[Q * 3 + i], MAT_FINAL_ASSEMBLY));
94: PetscCall(MatAssemblyEnd(array[Q * 3 + i], MAT_FINAL_ASSEMBLY));
95: PetscCall(MatSetRandom(array[Q * 3 + i], rctx));
96: }
97: PetscCall(MatGetSize(array[(Q + 1) * 1], NULL, &N));
98: PetscCall(MatGetSize(array[(Q + 1) * (Q - 1)], &M, NULL));
99: PetscCall(MatCreateBAIJ(PETSC_COMM_WORLD, 2, PETSC_DECIDE, PETSC_DECIDE, M, N, 0, NULL, 0, NULL, &a));
100: PetscCall(MatAssemblyBegin(a, MAT_FINAL_ASSEMBLY));
101: PetscCall(MatAssemblyEnd(a, MAT_FINAL_ASSEMBLY));
102: PetscCall(MatCreateHermitianTranspose(a, array + Q + Q - 1));
103: PetscCall(MatDestroy(&a));
104: PetscCall(MatDestroy(array + Q * Q - 1));
105: PetscCall(MatCreateNest(PETSC_COMM_WORLD, Q, NULL, Q, NULL, array, &A));
106: for (i = 0; i < Q; ++i) PetscCall(MatDestroy(array + (Q + 1) * i));
107: for (i = 2; i < Q; ++i) {
108: PetscCall(MatDestroy(array + i));
109: PetscCall(MatDestroy(array + Q * i));
110: }
111: for (i = 1; i < 3; ++i) PetscCall(MatDestroy(array + Q * 3 + i));
112: PetscCall(MatDestroy(array + Q + Q - 1));
113: PetscCall(KSPSetOperators(ksp, A, A));
114: PetscCall(MatNestGetISs(A, rows, NULL));
115: PetscCall(KSPGetPC(ksp, &pc));
116: PetscCall(PCSetType(pc, PCFIELDSPLIT));
117: M = 0;
118: for (j = 0; j < found; ++j) {
119: IS expand1, expand2;
120: PetscCall(ISDuplicate(rows[M], &expand1));
121: for (i = 1; i < outer[j]; ++i) {
122: PetscCall(ISExpand(expand1, rows[M + i], &expand2));
123: PetscCall(ISDestroy(&expand1));
124: expand1 = expand2;
125: }
126: M += outer[j];
127: PetscCall(PCFieldSplitSetIS(pc, NULL, expand1));
128: PetscCall(ISDestroy(&expand1));
129: }
130: PetscCall(KSPSetFromOptions(ksp));
131: flg = PETSC_FALSE;
132: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matmult", &flg, NULL));
133: if (flg) {
134: Mat D, E, F, H;
135: PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &D));
136: PetscCall(MatMultEqual(A, D, 10, &flg));
137: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatDense != MatNest");
138: PetscCall(MatZeroEntries(D));
139: PetscCall(MatConvert(A, MATDENSE, MAT_REUSE_MATRIX, &D));
140: PetscCall(MatMultEqual(A, D, 10, &flg));
141: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatDense != MatNest");
142: PetscCall(MatConvert(D, MATAIJ, MAT_INITIAL_MATRIX, &E));
143: PetscCall(MatMultEqual(E, D, 10, &flg));
144: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatDense != MatAIJ");
145: PetscCall(MatZeroEntries(E));
146: PetscCall(MatConvert(D, MATAIJ, MAT_REUSE_MATRIX, &E));
147: PetscCall(MatMultEqual(E, D, 10, &flg));
148: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatDense != MatAIJ");
149: PetscCall(MatConvert(E, MATDENSE, MAT_INPLACE_MATRIX, &E));
150: PetscCall(MatMultEqual(A, E, 10, &flg));
151: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatAIJ != MatNest");
152: PetscCall(MatConvert(D, MATAIJ, MAT_INPLACE_MATRIX, &D));
153: PetscCall(MatMultEqual(A, D, 10, &flg));
154: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatDense != MatNest");
155: PetscCall(MatDestroy(&E));
156: PetscCall(MatCreateHermitianTranspose(D, &E));
157: PetscCall(MatConvert(E, MATAIJ, MAT_INITIAL_MATRIX, &F));
158: PetscCall(MatMultEqual(E, F, 10, &flg));
159: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatHermitianTranspose != MatAIJ");
160: PetscCall(MatZeroEntries(F));
161: PetscCall(MatConvert(E, MATAIJ, MAT_REUSE_MATRIX, &F));
162: PetscCall(MatMultEqual(E, F, 10, &flg));
163: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatHermitianTranspose != MatAIJ");
164: PetscCall(MatDestroy(&F));
165: PetscCall(MatConvert(E, MATAIJ, MAT_INPLACE_MATRIX, &E));
166: PetscCall(MatCreateHermitianTranspose(D, &H));
167: PetscCall(MatMultEqual(E, H, 10, &flg));
168: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatHermitianTranspose != MatAIJ");
169: PetscCall(MatDestroy(&H));
170: PetscCall(MatDestroy(&E));
171: PetscCall(MatDestroy(&D));
172: }
173: PetscCall(KSPSetUp(ksp));
174: PetscCall(MatCreateVecs(A, &b, &x));
175: PetscCall(VecSetRandom(b, rctx));
176: PetscCall(VecGetSubVector(b, rows[Q - 1], &sub));
177: PetscCall(VecSet(sub, 0.0));
178: PetscCall(VecRestoreSubVector(b, rows[Q - 1], &sub));
179: PetscCall(KSPSolve(ksp, b, x));
180: PetscCall(VecDestroy(&b));
181: PetscCall(VecDestroy(&x));
182: PetscCall(PetscRandomDestroy(&rctx));
183: PetscCall(MatDestroy(&A));
184: PetscCall(KSPDestroy(&ksp));
185: PetscCall(PetscFree(outer));
186: PetscCall(PetscFinalize());
187: return 0;
188: }
190: /*TEST
192: test:
193: nsize: {{1 3}shared output}
194: suffix: wo_explicit_schur
195: requires: !complex !single
196: filter: sed -e "s/seq/mpi/g" -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g" -e "s/ 1 MPI process/ 3 MPI processes/g" -e "s/iterations [4-5]/iterations 4/g" -e "s/hermitiantranspose/transpose/g"
197: args: -outer_fieldsplit_sizes {{1,2,2 2,1,2 2,2,1}separate output} -ksp_view_mat -pc_type fieldsplit -ksp_converged_reason -fieldsplit_pc_type jacobi -test_matmult
199: test:
200: nsize: {{1 3}shared output}
201: suffix: w_explicit_schur
202: requires: !complex
203: filter: sed -e "s/seq/mpi/g" -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g" -e "s/ 1 MPI process/ 3 MPI processes/g" -e "s/iterations [1-2]/iterations 1/g" -e "s/hermitiantranspose/transpose/g"
204: args: -outer_fieldsplit_sizes {{1,4 2,3 3,2 4,1}separate output} -ksp_view_mat -pc_type fieldsplit -ksp_converged_reason -fieldsplit_pc_type jacobi -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition full
206: TEST*/