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