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;

 25:   MatMult(A, x, f);
 26:   return(0);
 27: }

 29: PetscErrorCode ComputeJacobianLinear(SNES snes, Vec x, Mat A, Mat J, void *ctx)
 30: {
 32:   return(0);
 33: }

 35: PetscErrorCode ConstructProblem1(Mat A, Vec b)
 36: {
 37:   PetscInt       rStart, rEnd, row;

 41:   VecSet(b, -3.0);
 42:   MatGetOwnershipRange(A, &rStart, &rEnd);
 43:   for (row = rStart; row < rEnd; ++row) {
 44:     PetscScalar val = 1.0;

 46:     MatSetValues(A, 1, &row, 1, &row, &val, INSERT_VALUES);
 47:   }
 48:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 49:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 50:   return(0);
 51: }

 53: PetscErrorCode CheckProblem1(Mat A, Vec b, Vec u)
 54: {
 55:   Vec            errorVec;
 56:   PetscReal      norm, error;

 60:   VecDuplicate(b, &errorVec);
 61:   VecWAXPY(errorVec, -1.0, b, u);
 62:   VecNorm(errorVec, NORM_2, &error);
 63:   VecNorm(b, NORM_2, &norm);
 64:   if (error/norm > 1e-12) SETERRQ1(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", error/norm);
 65:   VecDestroy(&errorVec);
 66:   return(0);
 67: }

 69: PetscErrorCode ConstructProblem2(Mat A, Vec b)
 70: {
 71:   PetscInt       N = 10, constraintSize = 4;
 72:   PetscInt       row;

 76:   VecSet(b, -3.0);
 77:   for (row = 0; row < constraintSize; ++row) {
 78:     PetscScalar vals[2] = {1.0, 1.0};
 79:     PetscInt    cols[2];

 81:     cols[0] = row; cols[1] = row + N - constraintSize;
 82:     MatSetValues(A, 1, &row, 2, cols, vals, INSERT_VALUES);
 83:   }
 84:   for (row = constraintSize; row < N - constraintSize; ++row) {
 85:     PetscScalar val = 1.0;

 87:     MatSetValues(A, 1, &row, 1, &row, &val, INSERT_VALUES);
 88:   }
 89:   for (row = N - constraintSize; row < N; ++row) {
 90:     PetscInt    col = row - (N - constraintSize);
 91:     PetscScalar val = 1.0;

 93:     MatSetValues(A, 1, &row, 1, &col, &val, INSERT_VALUES);
 94:   }
 95:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 96:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 97:   return(0);
 98: }

100: PetscErrorCode CheckProblem2(Mat A, Vec b, Vec u)
101: {
102:   PetscInt          N = 10, constraintSize = 4, r;
103:   PetscReal         norm, error;
104:   const PetscScalar *uArray, *bArray;
105:   PetscErrorCode    ierr;

108:   VecNorm(b, NORM_2, &norm);
109:   VecGetArrayRead(u, &uArray);
110:   VecGetArrayRead(b, &bArray);
111:   error = 0.0;
112:   for (r = 0; r < constraintSize; ++r) error += PetscRealPart(PetscSqr(uArray[r] - bArray[r + N-constraintSize]));

114:   if (error/norm > 1e-12) SETERRQ1(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", error/norm);
115:   error = 0.0;
116:   for (r = constraintSize; r < N - constraintSize; ++r) error += PetscRealPart(PetscSqr(uArray[r] - bArray[r]));

118:   if (error/norm > 1e-12) SETERRQ1(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", error/norm);
119:   error = 0.0;
120:   for (r = N - constraintSize; r < N; ++r) error += PetscRealPart(PetscSqr(uArray[r] - (bArray[r - (N-constraintSize)] - bArray[r])));

122:   if (error/norm > 1e-12) SETERRQ1(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", error/norm);
123:   VecRestoreArrayRead(u, &uArray);
124:   VecRestoreArrayRead(b, &bArray);
125:   return(0);
126: }

128: int main(int argc, char **argv)
129: {
130:   MPI_Comm       comm;
131:   SNES           snes;                 /* nonlinear solver */
132:   Vec            u,r,b;                /* solution, residual, and rhs vectors */
133:   Mat            A,J;                  /* Jacobian matrix */
134:   PetscInt       problem = 1, N = 10;

137:   PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
138:   comm = PETSC_COMM_WORLD;
139:   PetscOptionsGetInt(NULL,NULL, "-problem", &problem, NULL);
140:   VecCreate(comm, &u);
141:   VecSetSizes(u, PETSC_DETERMINE, N);
142:   VecSetFromOptions(u);
143:   VecDuplicate(u, &r);
144:   VecDuplicate(u, &b);

146:   MatCreate(comm, &A);
147:   MatSetSizes(A, PETSC_DETERMINE, PETSC_DETERMINE, N, N);
148:   MatSetFromOptions(A);
149:   MatSeqAIJSetPreallocation(A, 5, NULL);
150:   J    = A;

152:   switch (problem) {
153:   case 1:
154:     ConstructProblem1(A, b);
155:     break;
156:   case 2:
157:     ConstructProblem2(A, b);
158:     break;
159:   default:
160:     SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid problem number %d", problem);
161:   }

163:   SNESCreate(PETSC_COMM_WORLD, &snes);
164:   SNESSetJacobian(snes, A, J, ComputeJacobianLinear, NULL);
165:   SNESSetFunction(snes, r, ComputeFunctionLinear, A);
166:   SNESSetFromOptions(snes);

168:   SNESSolve(snes, b, u);
169:   VecView(u, NULL);

171:   switch (problem) {
172:   case 1:
173:     CheckProblem1(A, b, u);
174:     break;
175:   case 2:
176:     CheckProblem2(A, b, u);
177:     break;
178:   default:
179:     SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid problem number %d", problem);
180:   }

182:   if (A != J) {
183:     MatDestroy(&A);
184:   }
185:   MatDestroy(&J);
186:   VecDestroy(&u);
187:   VecDestroy(&r);
188:   VecDestroy(&b);
189:   SNESDestroy(&snes);
190:   PetscFinalize();
191:   return ierr;
192: }