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:   const PetscScalar *uArray, *bArray;
117:   PetscErrorCode    ierr;

120:   VecNorm(b, NORM_2, &norm);
121:   VecGetArrayRead(u, &uArray);
122:   VecGetArrayRead(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:   VecRestoreArrayRead(u, &uArray);
136:   VecRestoreArrayRead(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,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: }