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, void *ctx)
20: {
21: Mat A = (Mat) ctx;
24: MatMult(A, x, f);
25: return 0;
26: }
28: PetscErrorCode ComputeJacobianLinear(SNES snes, Vec x, Mat A, Mat J, void *ctx)
29: {
31: return 0;
32: }
34: PetscErrorCode ConstructProblem1(Mat A, Vec b)
35: {
36: PetscInt rStart, rEnd, row;
39: VecSet(b, -3.0);
40: MatGetOwnershipRange(A, &rStart, &rEnd);
41: for (row = rStart; row < rEnd; ++row) {
42: PetscScalar val = 1.0;
44: MatSetValues(A, 1, &row, 1, &row, &val, INSERT_VALUES);
45: }
46: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
47: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
48: return 0;
49: }
51: PetscErrorCode CheckProblem1(Mat A, Vec b, Vec u)
52: {
53: Vec errorVec;
54: PetscReal norm, error;
57: VecDuplicate(b, &errorVec);
58: VecWAXPY(errorVec, -1.0, b, u);
59: VecNorm(errorVec, NORM_2, &error);
60: VecNorm(b, NORM_2, &norm);
62: VecDestroy(&errorVec);
63: return 0;
64: }
66: PetscErrorCode ConstructProblem2(Mat A, Vec b)
67: {
68: PetscInt N = 10, constraintSize = 4;
69: PetscInt row;
72: VecSet(b, -3.0);
73: for (row = 0; row < constraintSize; ++row) {
74: PetscScalar vals[2] = {1.0, 1.0};
75: PetscInt cols[2];
77: cols[0] = row; cols[1] = row + N - constraintSize;
78: MatSetValues(A, 1, &row, 2, cols, vals, INSERT_VALUES);
79: }
80: for (row = constraintSize; row < N - constraintSize; ++row) {
81: PetscScalar val = 1.0;
83: MatSetValues(A, 1, &row, 1, &row, &val, INSERT_VALUES);
84: }
85: for (row = N - constraintSize; row < N; ++row) {
86: PetscInt col = row - (N - constraintSize);
87: PetscScalar val = 1.0;
89: MatSetValues(A, 1, &row, 1, &col, &val, INSERT_VALUES);
90: }
91: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
92: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
93: return 0;
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;
103: VecNorm(b, NORM_2, &norm);
104: VecGetArrayRead(u, &uArray);
105: VecGetArrayRead(b, &bArray);
106: error = 0.0;
107: for (r = 0; r < constraintSize; ++r) error += PetscRealPart(PetscSqr(uArray[r] - bArray[r + N-constraintSize]));
110: error = 0.0;
111: for (r = constraintSize; r < N - constraintSize; ++r) error += PetscRealPart(PetscSqr(uArray[r] - bArray[r]));
114: error = 0.0;
115: for (r = N - constraintSize; r < N; ++r) error += PetscRealPart(PetscSqr(uArray[r] - (bArray[r - (N-constraintSize)] - bArray[r])));
118: VecRestoreArrayRead(u, &uArray);
119: VecRestoreArrayRead(b, &bArray);
120: return 0;
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: PetscInitialize(&argc, &argv, NULL,help);
132: comm = PETSC_COMM_WORLD;
133: PetscOptionsGetInt(NULL,NULL, "-problem", &problem, NULL);
134: VecCreate(comm, &u);
135: VecSetSizes(u, PETSC_DETERMINE, N);
136: VecSetFromOptions(u);
137: VecDuplicate(u, &r);
138: VecDuplicate(u, &b);
140: MatCreate(comm, &A);
141: MatSetSizes(A, PETSC_DETERMINE, PETSC_DETERMINE, N, N);
142: MatSetFromOptions(A);
143: MatSeqAIJSetPreallocation(A, 5, NULL);
144: J = A;
146: switch (problem) {
147: case 1:
148: ConstructProblem1(A, b);
149: break;
150: case 2:
151: ConstructProblem2(A, b);
152: break;
153: default:
154: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid problem number %d", problem);
155: }
157: SNESCreate(PETSC_COMM_WORLD, &snes);
158: SNESSetJacobian(snes, A, J, ComputeJacobianLinear, NULL);
159: SNESSetFunction(snes, r, ComputeFunctionLinear, A);
160: SNESSetFromOptions(snes);
162: SNESSolve(snes, b, u);
163: VecView(u, NULL);
165: switch (problem) {
166: case 1:
167: CheckProblem1(A, b, u);
168: break;
169: case 2:
170: CheckProblem2(A, b, u);
171: break;
172: default:
173: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid problem number %d", problem);
174: }
176: if (A != J) {
177: MatDestroy(&A);
178: }
179: MatDestroy(&J);
180: VecDestroy(&u);
181: VecDestroy(&r);
182: VecDestroy(&b);
183: SNESDestroy(&snes);
184: PetscFinalize();
185: return 0;
186: }
188: /*TEST
190: test:
191: args: -snes_monitor
193: test:
194: suffix: 2
195: args: -problem 2 -pc_type jacobi -snes_monitor
197: TEST*/