1 #include <petsc.h>
2
3 static char help[] = "Solves a linear system with a MatNest and nested fields.\n\n";
4
5 #define Q 5 /* everything is hardwired for a 5x5 MatNest for now */
6
main(int argc,char ** args)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;
18
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 }
189
190 /*TEST
191
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
198
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
205
206 TEST*/
207