Actual source code: ex68.c
petsc-3.5.4 2015-05-23
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: */
21: PetscErrorCode ComputeFunctionLinear(SNES snes, Vec x, Vec f, void *ctx)
22: {
23: Mat A = (Mat) ctx;
27: MatMult(A, x, f);
28: return(0);
29: }
33: PetscErrorCode ComputeJacobianLinear(SNES snes, Vec x, Mat A, Mat J, void *ctx)
34: {
36: return(0);
37: }
41: PetscErrorCode ConstructProblem1(Mat A, Vec b)
42: {
43: PetscInt rStart, rEnd, row;
47: VecSet(b, -3.0);
48: MatGetOwnershipRange(A, &rStart, &rEnd);
49: for (row = rStart; row < rEnd; ++row) {
50: PetscScalar val = 1.0;
52: MatSetValues(A, 1, &row, 1, &row, &val, INSERT_VALUES);
53: }
54: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
55: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
56: return(0);
57: }
61: PetscErrorCode CheckProblem1(Mat A, Vec b, Vec u)
62: {
63: Vec errorVec;
64: PetscReal norm, error;
68: VecDuplicate(b, &errorVec);
69: VecWAXPY(errorVec, -1.0, b, u);
70: VecNorm(errorVec, NORM_2, &error);
71: VecNorm(b, NORM_2, &norm);
72: if (error/norm > 1e-12) SETERRQ1(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", error/norm);
73: VecDestroy(&errorVec);
74: return(0);
75: }
79: PetscErrorCode ConstructProblem2(Mat A, Vec b)
80: {
81: PetscInt N = 10, constraintSize = 4;
82: PetscInt row;
86: VecSet(b, -3.0);
87: for (row = 0; row < constraintSize; ++row) {
88: PetscScalar vals[2] = {1.0, 1.0};
89: PetscInt cols[2];
91: cols[0] = row; cols[1] = row + N - constraintSize;
92: MatSetValues(A, 1, &row, 2, cols, vals, INSERT_VALUES);
93: }
94: for (row = constraintSize; row < N - constraintSize; ++row) {
95: PetscScalar val = 1.0;
97: MatSetValues(A, 1, &row, 1, &row, &val, INSERT_VALUES);
98: }
99: for (row = N - constraintSize; row < N; ++row) {
100: PetscInt col = row - (N - constraintSize);
101: PetscScalar val = 1.0;
103: MatSetValues(A, 1, &row, 1, &col, &val, INSERT_VALUES);
104: }
105: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
106: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
107: return(0);
108: }
112: PetscErrorCode CheckProblem2(Mat A, Vec b, Vec u)
113: {
114: PetscInt N = 10, constraintSize = 4, r;
115: PetscReal norm, error;
116: PetscScalar *uArray, *bArray;
120: VecNorm(b, NORM_2, &norm);
121: VecGetArray(u, &uArray);
122: VecGetArray(b, &bArray);
123: error = 0.0;
124: for (r = 0; r < constraintSize; ++r) error += PetscRealPart(PetscSqr(uArray[r] - bArray[r + N-constraintSize]));
126: if (error/norm > 1e-12) SETERRQ1(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", error/norm);
127: error = 0.0;
128: for (r = constraintSize; r < N - constraintSize; ++r) error += PetscRealPart(PetscSqr(uArray[r] - bArray[r]));
130: if (error/norm > 1e-12) SETERRQ1(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", error/norm);
131: error = 0.0;
132: for (r = N - constraintSize; r < N; ++r) error += PetscRealPart(PetscSqr(uArray[r] - (bArray[r - (N-constraintSize)] - bArray[r])));
134: if (error/norm > 1e-12) SETERRQ1(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", error/norm);
135: VecRestoreArray(u, &uArray);
136: VecRestoreArray(b, &bArray);
137: return(0);
138: }
142: int main(int argc, char **argv)
143: {
144: MPI_Comm comm;
145: SNES snes; /* nonlinear solver */
146: Vec u,r,b; /* solution, residual, and rhs vectors */
147: Mat A,J; /* Jacobian matrix */
148: PetscInt problem = 1, N = 10;
151: PetscInitialize(&argc, &argv, NULL, help);
152: comm = PETSC_COMM_WORLD;
153: PetscOptionsGetInt(NULL, "-problem", &problem, NULL);
154: VecCreate(comm, &u);
155: VecSetSizes(u, PETSC_DETERMINE, N);
156: VecSetFromOptions(u);
157: VecDuplicate(u, &r);
158: VecDuplicate(u, &b);
160: MatCreate(comm, &A);
161: MatSetSizes(A, PETSC_DETERMINE, PETSC_DETERMINE, N, N);
162: MatSetFromOptions(A);
163: MatSeqAIJSetPreallocation(A, 5, NULL);
164: J = A;
166: switch (problem) {
167: case 1:
168: ConstructProblem1(A, b);
169: break;
170: case 2:
171: ConstructProblem2(A, b);
172: break;
173: default:
174: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid problem number %d", problem);
175: }
177: SNESCreate(PETSC_COMM_WORLD, &snes);
178: SNESSetJacobian(snes, A, J, ComputeJacobianLinear, NULL);
179: SNESSetFunction(snes, r, ComputeFunctionLinear, A);
180: SNESSetFromOptions(snes);
182: SNESSolve(snes, b, u);
183: VecView(u, NULL);
185: switch (problem) {
186: case 1:
187: CheckProblem1(A, b, u);
188: break;
189: case 2:
190: CheckProblem2(A, b, u);
191: break;
192: default:
193: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid problem number %d", problem);
194: }
196: if (A != J) {
197: MatDestroy(&A);
198: }
199: MatDestroy(&J);
200: VecDestroy(&u);
201: VecDestroy(&r);
202: VecDestroy(&b);
203: SNESDestroy(&snes);
204: PetscFinalize();
205: return 0;
206: }