Actual source code: ex81a.c
1: #include <petsc.h>
3: static char help[] = "Solves a linear system with a MatNest and PCFIELDSPLIT with fields defined from the command line.\n\n";
4: /* similar to ex81.c except the PCFIELDSPLIT fields are defined from the command line with fields instead of hardwired IS from the MATNEST */
6: #define Q 5 /* everything is hardwired for a 5x5 MatNest for now */
8: int main(int argc, char **args)
9: {
10: KSP ksp;
11: PC pc;
12: Mat array[Q * Q], A, a;
13: Vec b, x, sub;
14: IS rows[Q];
15: PetscInt i, M, N;
16: PetscMPIInt size;
17: PetscRandom rctx;
19: PetscFunctionBeginUser;
20: PetscCall(PetscInitialize(&argc, &args, NULL, help));
21: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
22: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
23: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
24: size = PetscMax(3, size);
25: for (i = 0; i < Q * Q; ++i) array[i] = NULL;
26: for (i = 0; i < Q; ++i) {
27: if (i == 0) {
28: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, size, size, 1, NULL, 0, NULL, array + (Q + 1) * i));
29: } else if (i == 1 || i == 3) {
30: PetscCall(MatCreateSBAIJ(PETSC_COMM_WORLD, 2, PETSC_DECIDE, PETSC_DECIDE, size, size, 1, NULL, 0, NULL, array + (Q + 1) * i));
31: } else if (i == 2 || i == 4) {
32: PetscCall(MatCreateBAIJ(PETSC_COMM_WORLD, 2, PETSC_DECIDE, PETSC_DECIDE, size, size, 1, NULL, 0, NULL, array + (Q + 1) * i));
33: }
34: PetscCall(MatAssemblyBegin(array[(Q + 1) * i], MAT_FINAL_ASSEMBLY));
35: PetscCall(MatAssemblyEnd(array[(Q + 1) * i], MAT_FINAL_ASSEMBLY));
36: PetscCall(MatShift(array[(Q + 1) * i], 100 + i + 1));
37: if (i == 3) {
38: PetscCall(MatDuplicate(array[(Q + 1) * i], MAT_COPY_VALUES, &a));
39: PetscCall(MatDestroy(array + (Q + 1) * i));
40: PetscCall(MatCreateHermitianTranspose(a, array + (Q + 1) * i));
41: PetscCall(MatDestroy(&a));
42: }
43: size *= 2;
44: }
45: PetscCall(MatGetSize(array[0], &M, NULL));
46: for (i = 2; i < Q; ++i) {
47: PetscCall(MatGetSize(array[(Q + 1) * i], NULL, &N));
48: if (i != Q - 1) {
49: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, i == 3 ? N : M, i == 3 ? M : N, 0, NULL, 0, NULL, array + i));
50: } else {
51: PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, M, N, NULL, array + i));
52: }
53: PetscCall(MatAssemblyBegin(array[i], MAT_FINAL_ASSEMBLY));
54: PetscCall(MatAssemblyEnd(array[i], MAT_FINAL_ASSEMBLY));
55: PetscCall(MatSetRandom(array[i], rctx));
56: if (i == 3) {
57: PetscCall(MatDuplicate(array[i], MAT_COPY_VALUES, &a));
58: PetscCall(MatDestroy(array + i));
59: PetscCall(MatCreateHermitianTranspose(a, array + i));
60: PetscCall(MatDestroy(&a));
61: }
62: }
63: PetscCall(MatGetSize(array[0], NULL, &N));
64: for (i = 2; i < Q; i += 2) {
65: PetscCall(MatGetSize(array[(Q + 1) * i], &M, NULL));
66: if (i != Q - 1) {
67: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, M, N, 2, NULL, 2, NULL, array + Q * i));
68: } else {
69: PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, N, M, NULL, array + Q * i));
70: }
71: PetscCall(MatAssemblyBegin(array[Q * i], MAT_FINAL_ASSEMBLY));
72: PetscCall(MatAssemblyEnd(array[Q * i], MAT_FINAL_ASSEMBLY));
73: PetscCall(MatSetRandom(array[Q * i], rctx));
74: if (i == Q - 1) {
75: PetscCall(MatDuplicate(array[Q * i], MAT_COPY_VALUES, &a));
76: PetscCall(MatDestroy(array + Q * i));
77: PetscCall(MatCreateHermitianTranspose(a, array + Q * i));
78: PetscCall(MatDestroy(&a));
79: }
80: }
81: PetscCall(MatGetSize(array[(Q + 1) * 3], &M, NULL));
82: for (i = 1; i < 3; ++i) {
83: PetscCall(MatGetSize(array[(Q + 1) * i], NULL, &N));
84: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, M, N, 2, NULL, 2, NULL, array + Q * 3 + i));
85: PetscCall(MatAssemblyBegin(array[Q * 3 + i], MAT_FINAL_ASSEMBLY));
86: PetscCall(MatAssemblyEnd(array[Q * 3 + i], MAT_FINAL_ASSEMBLY));
87: PetscCall(MatSetRandom(array[Q * 3 + i], rctx));
88: }
89: PetscCall(MatGetSize(array[(Q + 1) * 1], NULL, &N));
90: PetscCall(MatGetSize(array[(Q + 1) * (Q - 1)], &M, NULL));
91: PetscCall(MatCreateBAIJ(PETSC_COMM_WORLD, 2, PETSC_DECIDE, PETSC_DECIDE, M, N, 0, NULL, 0, NULL, &a));
92: PetscCall(MatAssemblyBegin(a, MAT_FINAL_ASSEMBLY));
93: PetscCall(MatAssemblyEnd(a, MAT_FINAL_ASSEMBLY));
94: PetscCall(MatCreateHermitianTranspose(a, array + Q + Q - 1));
95: PetscCall(MatDestroy(&a));
96: PetscCall(MatDestroy(array + Q * Q - 1));
97: PetscCall(MatCreateNest(PETSC_COMM_WORLD, Q, NULL, Q, NULL, array, &A));
98: for (i = 0; i < Q; ++i) PetscCall(MatDestroy(array + (Q + 1) * i));
99: for (i = 2; i < Q; ++i) {
100: PetscCall(MatDestroy(array + i));
101: PetscCall(MatDestroy(array + Q * i));
102: }
103: for (i = 1; i < 3; ++i) PetscCall(MatDestroy(array + Q * 3 + i));
104: PetscCall(MatDestroy(array + Q + Q - 1));
105: PetscCall(KSPSetOperators(ksp, A, A));
106: PetscCall(MatNestGetISs(A, rows, NULL));
107: PetscCall(KSPGetPC(ksp, &pc));
108: PetscCall(PCSetType(pc, PCFIELDSPLIT));
109: PetscCall(KSPSetFromOptions(ksp));
110: PetscCall(MatCreateVecs(A, &b, &x));
111: PetscCall(VecSetRandom(b, rctx));
112: PetscCall(VecGetSubVector(b, rows[Q - 1], &sub));
113: PetscCall(VecSet(sub, 0.0));
114: PetscCall(VecRestoreSubVector(b, rows[Q - 1], &sub));
115: PetscCall(KSPSolve(ksp, b, x));
116: PetscCall(VecDestroy(&b));
117: PetscCall(VecDestroy(&x));
118: PetscCall(PetscRandomDestroy(&rctx));
119: PetscCall(MatDestroy(&A));
120: PetscCall(KSPDestroy(&ksp));
121: PetscCall(PetscFinalize());
122: return 0;
123: }
125: /*TEST
127: test:
128: nsize: 3
129: suffix: 1
130: requires: !complex !single
131: filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g" -e "s/iterations [4-5]/iterations 4/g" -e "s/hermitiantranspose/transpose/g"
132: args: -pc_fieldsplit_0_fields 0,1 -pc_fieldsplit_1_fields 2,3 -pc_fieldsplit_2_fields 4 -pc_type fieldsplit -ksp_converged_reason -fieldsplit_pc_type jacobi -ksp_view
134: test:
135: suffix: 2
136: nsize: 3
137: requires: !complex !single
138: filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g" -e "s/iterations [4-5]/iterations 4/g" -e "s/hermitiantranspose/transpose/g"
139: args: -pc_type fieldsplit -ksp_converged_reason -fieldsplit_pc_type jacobi -ksp_view
141: TEST*/