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;
20: PetscInitialize(&argc, &args, NULL, help);
21: MPI_Comm_size(PETSC_COMM_WORLD, &size);
22: PetscMalloc1(found, &outer);
23: PetscOptionsGetIntArray(NULL, NULL, "-outer_fieldsplit_sizes", outer, &found, &flg);
24: PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
25: if (flg) {
27: j = 0;
28: for (i = 0; i < found; ++i) j += outer[i];
30: }
31: 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: 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: 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: MatCreateBAIJ(PETSC_COMM_WORLD, 2, PETSC_DECIDE, PETSC_DECIDE, size, size, 1, NULL, 0, NULL, array + (Q + 1) * i);
41: }
42: MatAssemblyBegin(array[(Q + 1) * i], MAT_FINAL_ASSEMBLY);
43: MatAssemblyEnd(array[(Q + 1) * i], MAT_FINAL_ASSEMBLY);
44: MatShift(array[(Q + 1) * i], 100 + i + 1);
45: if (i == 3) {
46: MatDuplicate(array[(Q + 1) * i], MAT_COPY_VALUES, &a);
47: MatDestroy(array + (Q + 1) * i);
48: MatCreateHermitianTranspose(a, array + (Q + 1) * i);
49: MatDestroy(&a);
50: }
51: size *= 2;
52: }
53: MatGetSize(array[0], &M, NULL);
54: for (i = 2; i < Q; ++i) {
55: MatGetSize(array[(Q + 1) * i], NULL, &N);
56: if (i != Q - 1) {
57: 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: MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, M, N, NULL, array + i);
60: }
61: MatAssemblyBegin(array[i], MAT_FINAL_ASSEMBLY);
62: MatAssemblyEnd(array[i], MAT_FINAL_ASSEMBLY);
63: MatSetRandom(array[i], rctx);
64: if (i == 3) {
65: MatDuplicate(array[i], MAT_COPY_VALUES, &a);
66: MatDestroy(array + i);
67: MatCreateHermitianTranspose(a, array + i);
68: MatDestroy(&a);
69: }
70: }
71: MatGetSize(array[0], NULL, &N);
72: for (i = 2; i < Q; i += 2) {
73: MatGetSize(array[(Q + 1) * i], &M, NULL);
74: if (i != Q - 1) {
75: MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, M, N, 2, NULL, 2, NULL, array + Q * i);
76: } else {
77: MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, N, M, NULL, array + Q * i);
78: }
79: MatAssemblyBegin(array[Q * i], MAT_FINAL_ASSEMBLY);
80: MatAssemblyEnd(array[Q * i], MAT_FINAL_ASSEMBLY);
81: MatSetRandom(array[Q * i], rctx);
82: if (i == Q - 1) {
83: MatDuplicate(array[Q * i], MAT_COPY_VALUES, &a);
84: MatDestroy(array + Q * i);
85: MatCreateHermitianTranspose(a, array + Q * i);
86: MatDestroy(&a);
87: }
88: }
89: MatGetSize(array[(Q + 1) * 3], &M, NULL);
90: for (i = 1; i < 3; ++i) {
91: MatGetSize(array[(Q + 1) * i], NULL, &N);
92: MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, M, N, 2, NULL, 2, NULL, array + Q * 3 + i);
93: MatAssemblyBegin(array[Q * 3 + i], MAT_FINAL_ASSEMBLY);
94: MatAssemblyEnd(array[Q * 3 + i], MAT_FINAL_ASSEMBLY);
95: MatSetRandom(array[Q * 3 + i], rctx);
96: }
97: MatGetSize(array[(Q + 1) * 1], NULL, &N);
98: MatGetSize(array[(Q + 1) * (Q - 1)], &M, NULL);
99: MatCreateBAIJ(PETSC_COMM_WORLD, 2, PETSC_DECIDE, PETSC_DECIDE, M, N, 0, NULL, 0, NULL, &a);
100: MatAssemblyBegin(a, MAT_FINAL_ASSEMBLY);
101: MatAssemblyEnd(a, MAT_FINAL_ASSEMBLY);
102: MatCreateHermitianTranspose(a, array + Q + Q - 1);
103: MatDestroy(&a);
104: MatDestroy(array + Q * Q - 1);
105: MatCreateNest(PETSC_COMM_WORLD, Q, NULL, Q, NULL, array, &A);
106: for (i = 0; i < Q; ++i) MatDestroy(array + (Q + 1) * i);
107: for (i = 2; i < Q; ++i) {
108: MatDestroy(array + i);
109: MatDestroy(array + Q * i);
110: }
111: for (i = 1; i < 3; ++i) MatDestroy(array + Q * 3 + i);
112: MatDestroy(array + Q + Q - 1);
113: KSPSetOperators(ksp, A, A);
114: MatNestGetISs(A, rows, NULL);
115: KSPGetPC(ksp, &pc);
116: PCSetType(pc, PCFIELDSPLIT);
117: M = 0;
118: for (j = 0; j < found; ++j) {
119: IS expand1, expand2;
120: ISDuplicate(rows[M], &expand1);
121: for (i = 1; i < outer[j]; ++i) {
122: ISExpand(expand1, rows[M + i], &expand2);
123: ISDestroy(&expand1);
124: expand1 = expand2;
125: }
126: M += outer[j];
127: PCFieldSplitSetIS(pc, NULL, expand1);
128: ISDestroy(&expand1);
129: }
130: KSPSetFromOptions(ksp);
131: flg = PETSC_FALSE;
132: PetscOptionsGetBool(NULL, NULL, "-test_matmult", &flg, NULL);
133: if (flg) {
134: Mat D, E, F, H;
135: MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &D);
136: MatMultEqual(A, D, 10, &flg);
138: MatZeroEntries(D);
139: MatConvert(A, MATDENSE, MAT_REUSE_MATRIX, &D);
140: MatMultEqual(A, D, 10, &flg);
142: MatConvert(D, MATAIJ, MAT_INITIAL_MATRIX, &E);
143: MatMultEqual(E, D, 10, &flg);
145: MatZeroEntries(E);
146: MatConvert(D, MATAIJ, MAT_REUSE_MATRIX, &E);
147: MatMultEqual(E, D, 10, &flg);
149: MatConvert(E, MATDENSE, MAT_INPLACE_MATRIX, &E);
150: MatMultEqual(A, E, 10, &flg);
152: MatConvert(D, MATAIJ, MAT_INPLACE_MATRIX, &D);
153: MatMultEqual(A, D, 10, &flg);
155: MatDestroy(&E);
156: MatCreateHermitianTranspose(D, &E);
157: MatConvert(E, MATAIJ, MAT_INITIAL_MATRIX, &F);
158: MatMultEqual(E, F, 10, &flg);
160: MatZeroEntries(F);
161: MatConvert(E, MATAIJ, MAT_REUSE_MATRIX, &F);
162: MatMultEqual(E, F, 10, &flg);
164: MatDestroy(&F);
165: MatConvert(E, MATAIJ, MAT_INPLACE_MATRIX, &E);
166: MatCreateHermitianTranspose(D, &H);
167: MatMultEqual(E, H, 10, &flg);
169: MatDestroy(&H);
170: MatDestroy(&E);
171: MatDestroy(&D);
172: }
173: KSPSetUp(ksp);
174: MatCreateVecs(A, &b, &x);
175: VecSetRandom(b, rctx);
176: VecGetSubVector(b, rows[Q - 1], &sub);
177: VecSet(sub, 0.0);
178: VecRestoreSubVector(b, rows[Q - 1], &sub);
179: KSPSolve(ksp, b, x);
180: VecDestroy(&b);
181: VecDestroy(&x);
182: PetscRandomDestroy(&rctx);
183: MatDestroy(&A);
184: KSPDestroy(&ksp);
185: PetscFree(outer);
186: PetscFinalize();
187: return 0;
188: }
190: /*TEST
192: test:
193: nsize: {{1 3}shared output}
194: suffix: wo_explicit_schur
195: requires: !complex
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*/