1 static char help[] = "Test problems for Schur complement solvers.\n\n\n"; 2 3 #include <petscsnes.h> 4 5 /* 6 Test 1: 7 I u = b 8 9 solution: u = b 10 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 / 15 16 solution: u_1 = b_3, u_2 = b_2, u_3 = b_1 - b_3 17 */ 18 19 PetscErrorCode ComputeFunctionLinear(SNES snes, Vec x, Vec f, void *ctx) { 20 Mat A = (Mat)ctx; 21 22 PetscFunctionBeginUser; 23 PetscCall(MatMult(A, x, f)); 24 PetscFunctionReturn(0); 25 } 26 27 PetscErrorCode ComputeJacobianLinear(SNES snes, Vec x, Mat A, Mat J, void *ctx) { 28 PetscFunctionBeginUser; 29 PetscFunctionReturn(0); 30 } 31 32 PetscErrorCode ConstructProblem1(Mat A, Vec b) { 33 PetscInt rStart, rEnd, row; 34 35 PetscFunctionBeginUser; 36 PetscCall(VecSet(b, -3.0)); 37 PetscCall(MatGetOwnershipRange(A, &rStart, &rEnd)); 38 for (row = rStart; row < rEnd; ++row) { 39 PetscScalar val = 1.0; 40 41 PetscCall(MatSetValues(A, 1, &row, 1, &row, &val, INSERT_VALUES)); 42 } 43 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 44 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 45 PetscFunctionReturn(0); 46 } 47 48 PetscErrorCode CheckProblem1(Mat A, Vec b, Vec u) { 49 Vec errorVec; 50 PetscReal norm, error; 51 52 PetscFunctionBeginUser; 53 PetscCall(VecDuplicate(b, &errorVec)); 54 PetscCall(VecWAXPY(errorVec, -1.0, b, u)); 55 PetscCall(VecNorm(errorVec, NORM_2, &error)); 56 PetscCall(VecNorm(b, NORM_2, &norm)); 57 PetscCheck(error / norm <= 1000. * PETSC_MACHINE_EPSILON, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", (double)(error / norm)); 58 PetscCall(VecDestroy(&errorVec)); 59 PetscFunctionReturn(0); 60 } 61 62 PetscErrorCode ConstructProblem2(Mat A, Vec b) { 63 PetscInt N = 10, constraintSize = 4; 64 PetscInt row; 65 66 PetscFunctionBeginUser; 67 PetscCall(VecSet(b, -3.0)); 68 for (row = 0; row < constraintSize; ++row) { 69 PetscScalar vals[2] = {1.0, 1.0}; 70 PetscInt cols[2]; 71 72 cols[0] = row; 73 cols[1] = row + N - constraintSize; 74 PetscCall(MatSetValues(A, 1, &row, 2, cols, vals, INSERT_VALUES)); 75 } 76 for (row = constraintSize; row < N - constraintSize; ++row) { 77 PetscScalar val = 1.0; 78 79 PetscCall(MatSetValues(A, 1, &row, 1, &row, &val, INSERT_VALUES)); 80 } 81 for (row = N - constraintSize; row < N; ++row) { 82 PetscInt col = row - (N - constraintSize); 83 PetscScalar val = 1.0; 84 85 PetscCall(MatSetValues(A, 1, &row, 1, &col, &val, INSERT_VALUES)); 86 } 87 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 88 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 89 PetscFunctionReturn(0); 90 } 91 92 PetscErrorCode CheckProblem2(Mat A, Vec b, Vec u) { 93 PetscInt N = 10, constraintSize = 4, r; 94 PetscReal norm, error; 95 const PetscScalar *uArray, *bArray; 96 97 PetscFunctionBeginUser; 98 PetscCall(VecNorm(b, NORM_2, &norm)); 99 PetscCall(VecGetArrayRead(u, &uArray)); 100 PetscCall(VecGetArrayRead(b, &bArray)); 101 error = 0.0; 102 for (r = 0; r < constraintSize; ++r) error += PetscRealPart(PetscSqr(uArray[r] - bArray[r + N - constraintSize])); 103 104 PetscCheck(error / norm <= 10000 * PETSC_MACHINE_EPSILON, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", (double)(error / norm)); 105 error = 0.0; 106 for (r = constraintSize; r < N - constraintSize; ++r) error += PetscRealPart(PetscSqr(uArray[r] - bArray[r])); 107 108 PetscCheck(error / norm <= 10000 * PETSC_MACHINE_EPSILON, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", (double)(error / norm)); 109 error = 0.0; 110 for (r = N - constraintSize; r < N; ++r) error += PetscRealPart(PetscSqr(uArray[r] - (bArray[r - (N - constraintSize)] - bArray[r]))); 111 112 PetscCheck(error / norm <= 10000 * PETSC_MACHINE_EPSILON, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", (double)(error / norm)); 113 PetscCall(VecRestoreArrayRead(u, &uArray)); 114 PetscCall(VecRestoreArrayRead(b, &bArray)); 115 PetscFunctionReturn(0); 116 } 117 118 int main(int argc, char **argv) { 119 MPI_Comm comm; 120 SNES snes; /* nonlinear solver */ 121 Vec u, r, b; /* solution, residual, and rhs vectors */ 122 Mat A, J; /* Jacobian matrix */ 123 PetscInt problem = 1, N = 10; 124 125 PetscFunctionBeginUser; 126 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 127 comm = PETSC_COMM_WORLD; 128 PetscCall(PetscOptionsGetInt(NULL, NULL, "-problem", &problem, NULL)); 129 PetscCall(VecCreate(comm, &u)); 130 PetscCall(VecSetSizes(u, PETSC_DETERMINE, N)); 131 PetscCall(VecSetFromOptions(u)); 132 PetscCall(VecDuplicate(u, &r)); 133 PetscCall(VecDuplicate(u, &b)); 134 135 PetscCall(MatCreate(comm, &A)); 136 PetscCall(MatSetSizes(A, PETSC_DETERMINE, PETSC_DETERMINE, N, N)); 137 PetscCall(MatSetFromOptions(A)); 138 PetscCall(MatSeqAIJSetPreallocation(A, 5, NULL)); 139 J = A; 140 141 switch (problem) { 142 case 1: PetscCall(ConstructProblem1(A, b)); break; 143 case 2: PetscCall(ConstructProblem2(A, b)); break; 144 default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid problem number %" PetscInt_FMT, problem); 145 } 146 147 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 148 PetscCall(SNESSetJacobian(snes, A, J, ComputeJacobianLinear, NULL)); 149 PetscCall(SNESSetFunction(snes, r, ComputeFunctionLinear, A)); 150 PetscCall(SNESSetFromOptions(snes)); 151 152 PetscCall(SNESSolve(snes, b, u)); 153 PetscCall(VecView(u, NULL)); 154 155 switch (problem) { 156 case 1: PetscCall(CheckProblem1(A, b, u)); break; 157 case 2: PetscCall(CheckProblem2(A, b, u)); break; 158 default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid problem number %" PetscInt_FMT, problem); 159 } 160 161 if (A != J) { PetscCall(MatDestroy(&A)); } 162 PetscCall(MatDestroy(&J)); 163 PetscCall(VecDestroy(&u)); 164 PetscCall(VecDestroy(&r)); 165 PetscCall(VecDestroy(&b)); 166 PetscCall(SNESDestroy(&snes)); 167 PetscCall(PetscFinalize()); 168 return 0; 169 } 170 171 /*TEST 172 173 test: 174 args: -snes_monitor 175 176 test: 177 suffix: 2 178 args: -problem 2 -pc_type jacobi -snes_monitor 179 180 TEST*/ 181