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*/