1 static char help[] = "Test problems for Schur complement solvers.\n\n\n";
2
3 #include <petscsnes.h>
4
5 /*
6 Test 1:
7 I u = b
8
9 solution: u = b
10
11 Test 2:
12 / I 0 I \ / u_1 \ / b_1 \
13 | 0 I 0 | | u_2 | = | b_2 |
14 \ I 0 0 / \ u_3 / \ b_3 /
15
16 solution: u_1 = b_3, u_2 = b_2, u_3 = b_1 - b_3
17 */
18
ComputeFunctionLinear(SNES snes,Vec x,Vec f,PetscCtx ctx)19 PetscErrorCode ComputeFunctionLinear(SNES snes, Vec x, Vec f, PetscCtx ctx)
20 {
21 Mat A = (Mat)ctx;
22
23 PetscFunctionBeginUser;
24 PetscCall(MatMult(A, x, f));
25 PetscFunctionReturn(PETSC_SUCCESS);
26 }
27
ComputeJacobianLinear(SNES snes,Vec x,Mat A,Mat J,PetscCtx ctx)28 PetscErrorCode ComputeJacobianLinear(SNES snes, Vec x, Mat A, Mat J, PetscCtx ctx)
29 {
30 PetscFunctionBeginUser;
31 PetscFunctionReturn(PETSC_SUCCESS);
32 }
33
ConstructProblem1(Mat A,Vec b)34 PetscErrorCode ConstructProblem1(Mat A, Vec b)
35 {
36 PetscInt rStart, rEnd, row;
37
38 PetscFunctionBeginUser;
39 PetscCall(VecSet(b, -3.0));
40 PetscCall(MatGetOwnershipRange(A, &rStart, &rEnd));
41 for (row = rStart; row < rEnd; ++row) {
42 PetscScalar val = 1.0;
43
44 PetscCall(MatSetValues(A, 1, &row, 1, &row, &val, INSERT_VALUES));
45 }
46 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
47 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
48 PetscFunctionReturn(PETSC_SUCCESS);
49 }
50
CheckProblem1(Mat A,Vec b,Vec u)51 PetscErrorCode CheckProblem1(Mat A, Vec b, Vec u)
52 {
53 Vec errorVec;
54 PetscReal norm, error;
55
56 PetscFunctionBeginUser;
57 PetscCall(VecDuplicate(b, &errorVec));
58 PetscCall(VecWAXPY(errorVec, -1.0, b, u));
59 PetscCall(VecNorm(errorVec, NORM_2, &error));
60 PetscCall(VecNorm(b, NORM_2, &norm));
61 PetscCheck(error / norm <= 1000. * PETSC_MACHINE_EPSILON, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", (double)(error / norm));
62 PetscCall(VecDestroy(&errorVec));
63 PetscFunctionReturn(PETSC_SUCCESS);
64 }
65
ConstructProblem2(Mat A,Vec b)66 PetscErrorCode ConstructProblem2(Mat A, Vec b)
67 {
68 PetscInt N = 10, constraintSize = 4;
69 PetscInt row;
70
71 PetscFunctionBeginUser;
72 PetscCall(VecSet(b, -3.0));
73 for (row = 0; row < constraintSize; ++row) {
74 PetscScalar vals[2] = {1.0, 1.0};
75 PetscInt cols[2];
76
77 cols[0] = row;
78 cols[1] = row + N - constraintSize;
79 PetscCall(MatSetValues(A, 1, &row, 2, cols, vals, INSERT_VALUES));
80 }
81 for (row = constraintSize; row < N - constraintSize; ++row) {
82 PetscScalar val = 1.0;
83
84 PetscCall(MatSetValues(A, 1, &row, 1, &row, &val, INSERT_VALUES));
85 }
86 for (row = N - constraintSize; row < N; ++row) {
87 PetscInt col = row - (N - constraintSize);
88 PetscScalar val = 1.0;
89
90 PetscCall(MatSetValues(A, 1, &row, 1, &col, &val, INSERT_VALUES));
91 }
92 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
93 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
94 PetscFunctionReturn(PETSC_SUCCESS);
95 }
96
CheckProblem2(Mat A,Vec b,Vec u)97 PetscErrorCode CheckProblem2(Mat A, Vec b, Vec u)
98 {
99 PetscInt N = 10, constraintSize = 4, r;
100 PetscReal norm, error;
101 const PetscScalar *uArray, *bArray;
102
103 PetscFunctionBeginUser;
104 PetscCall(VecNorm(b, NORM_2, &norm));
105 PetscCall(VecGetArrayRead(u, &uArray));
106 PetscCall(VecGetArrayRead(b, &bArray));
107 error = 0.0;
108 for (r = 0; r < constraintSize; ++r) error += PetscRealPart(PetscSqr(uArray[r] - bArray[r + N - constraintSize]));
109
110 PetscCheck(error / norm <= 10000 * PETSC_MACHINE_EPSILON, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", (double)(error / norm));
111 error = 0.0;
112 for (r = constraintSize; r < N - constraintSize; ++r) error += PetscRealPart(PetscSqr(uArray[r] - bArray[r]));
113
114 PetscCheck(error / norm <= 10000 * PETSC_MACHINE_EPSILON, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", (double)(error / norm));
115 error = 0.0;
116 for (r = N - constraintSize; r < N; ++r) error += PetscRealPart(PetscSqr(uArray[r] - (bArray[r - (N - constraintSize)] - bArray[r])));
117
118 PetscCheck(error / norm <= 10000 * PETSC_MACHINE_EPSILON, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", (double)(error / norm));
119 PetscCall(VecRestoreArrayRead(u, &uArray));
120 PetscCall(VecRestoreArrayRead(b, &bArray));
121 PetscFunctionReturn(PETSC_SUCCESS);
122 }
123
main(int argc,char ** argv)124 int main(int argc, char **argv)
125 {
126 MPI_Comm comm;
127 SNES snes; /* nonlinear solver */
128 Vec u, r, b; /* solution, residual, and rhs vectors */
129 Mat A, J; /* Jacobian matrix */
130 PetscInt problem = 1, N = 10;
131
132 PetscFunctionBeginUser;
133 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
134 comm = PETSC_COMM_WORLD;
135 PetscCall(PetscOptionsGetInt(NULL, NULL, "-problem", &problem, NULL));
136 PetscCall(VecCreate(comm, &u));
137 PetscCall(VecSetSizes(u, PETSC_DETERMINE, N));
138 PetscCall(VecSetFromOptions(u));
139 PetscCall(VecDuplicate(u, &r));
140 PetscCall(VecDuplicate(u, &b));
141
142 PetscCall(MatCreate(comm, &A));
143 PetscCall(MatSetSizes(A, PETSC_DETERMINE, PETSC_DETERMINE, N, N));
144 PetscCall(MatSetFromOptions(A));
145 PetscCall(MatSeqAIJSetPreallocation(A, 5, NULL));
146 J = A;
147
148 switch (problem) {
149 case 1:
150 PetscCall(ConstructProblem1(A, b));
151 break;
152 case 2:
153 PetscCall(ConstructProblem2(A, b));
154 break;
155 default:
156 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid problem number %" PetscInt_FMT, problem);
157 }
158
159 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
160 PetscCall(SNESSetJacobian(snes, A, J, ComputeJacobianLinear, NULL));
161 PetscCall(SNESSetFunction(snes, r, ComputeFunctionLinear, A));
162 PetscCall(SNESSetFromOptions(snes));
163
164 PetscCall(SNESSolve(snes, b, u));
165 PetscCall(VecView(u, NULL));
166
167 switch (problem) {
168 case 1:
169 PetscCall(CheckProblem1(A, b, u));
170 break;
171 case 2:
172 PetscCall(CheckProblem2(A, b, u));
173 break;
174 default:
175 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid problem number %" PetscInt_FMT, problem);
176 }
177
178 if (A != J) PetscCall(MatDestroy(&A));
179 PetscCall(MatDestroy(&J));
180 PetscCall(VecDestroy(&u));
181 PetscCall(VecDestroy(&r));
182 PetscCall(VecDestroy(&b));
183 PetscCall(SNESDestroy(&snes));
184 PetscCall(PetscFinalize());
185 return 0;
186 }
187
188 /*TEST
189
190 test:
191 args: -snes_monitor
192
193 test:
194 suffix: 2
195 args: -problem 2 -pc_type jacobi -snes_monitor
196
197 TEST*/
198