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;
 78:     cols[1] = row + N - constraintSize;
 79:     MatSetValues(A, 1, &row, 2, cols, vals, INSERT_VALUES);
 80:   }
 81:   for (row = constraintSize; row < N - constraintSize; ++row) {
 82:     PetscScalar val = 1.0;

 84:     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;

 90:     MatSetValues(A, 1, &row, 1, &col, &val, INSERT_VALUES);
 91:   }
 92:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 93:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 94:   return 0;
 95: }

 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;

104:   VecNorm(b, NORM_2, &norm);
105:   VecGetArrayRead(u, &uArray);
106:   VecGetArrayRead(b, &bArray);
107:   error = 0.0;
108:   for (r = 0; r < constraintSize; ++r) error += PetscRealPart(PetscSqr(uArray[r] - bArray[r + N - constraintSize]));

111:   error = 0.0;
112:   for (r = constraintSize; r < N - constraintSize; ++r) error += PetscRealPart(PetscSqr(uArray[r] - bArray[r]));

115:   error = 0.0;
116:   for (r = N - constraintSize; r < N; ++r) error += PetscRealPart(PetscSqr(uArray[r] - (bArray[r - (N - constraintSize)] - bArray[r])));

119:   VecRestoreArrayRead(u, &uArray);
120:   VecRestoreArrayRead(b, &bArray);
121:   return 0;
122: }

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;

133:   PetscInitialize(&argc, &argv, NULL, help);
134:   comm = PETSC_COMM_WORLD;
135:   PetscOptionsGetInt(NULL, NULL, "-problem", &problem, NULL);
136:   VecCreate(comm, &u);
137:   VecSetSizes(u, PETSC_DETERMINE, N);
138:   VecSetFromOptions(u);
139:   VecDuplicate(u, &r);
140:   VecDuplicate(u, &b);

142:   MatCreate(comm, &A);
143:   MatSetSizes(A, PETSC_DETERMINE, PETSC_DETERMINE, N, N);
144:   MatSetFromOptions(A);
145:   MatSeqAIJSetPreallocation(A, 5, NULL);
146:   J = A;

148:   switch (problem) {
149:   case 1:
150:     ConstructProblem1(A, b);
151:     break;
152:   case 2:
153:     ConstructProblem2(A, b);
154:     break;
155:   default:
156:     SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid problem number %" PetscInt_FMT, problem);
157:   }

159:   SNESCreate(PETSC_COMM_WORLD, &snes);
160:   SNESSetJacobian(snes, A, J, ComputeJacobianLinear, NULL);
161:   SNESSetFunction(snes, r, ComputeFunctionLinear, A);
162:   SNESSetFromOptions(snes);

164:   SNESSolve(snes, b, u);
165:   VecView(u, NULL);

167:   switch (problem) {
168:   case 1:
169:     CheckProblem1(A, b, u);
170:     break;
171:   case 2:
172:     CheckProblem2(A, b, u);
173:     break;
174:   default:
175:     SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid problem number %" PetscInt_FMT, problem);
176:   }

178:   if (A != J) MatDestroy(&A);
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*/