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, PetscCtx ctx) 20 { 21 Mat A = (Mat)ctx; 22 23 PetscFunctionBeginUser; 24 PetscCall(MatMult(A, x, f)); 25 PetscFunctionReturn(PETSC_SUCCESS); 26 } 27 28 PetscErrorCode ComputeJacobianLinear(SNES snes, Vec x, Mat A, Mat J, PetscCtx ctx) 29 { 30 PetscFunctionBeginUser; 31 PetscFunctionReturn(PETSC_SUCCESS); 32 } 33 34 PetscErrorCode ConstructProblem1(Mat A, Vec b) 35 { 36 PetscInt rStart, rEnd, row; 37 38 PetscFunctionBeginUser; 39 PetscCall(VecSet(b, -3.0)); 40 PetscCall(MatGetOwnershipRange(A, &rStart, &rEnd)); 41 for (row = rStart; row < rEnd; ++row) { 42 PetscScalar val = 1.0; 43 44 PetscCall(MatSetValues(A, 1, &row, 1, &row, &val, INSERT_VALUES)); 45 } 46 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 47 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 48 PetscFunctionReturn(PETSC_SUCCESS); 49 } 50 51 PetscErrorCode CheckProblem1(Mat A, Vec b, Vec u) 52 { 53 Vec errorVec; 54 PetscReal norm, error; 55 56 PetscFunctionBeginUser; 57 PetscCall(VecDuplicate(b, &errorVec)); 58 PetscCall(VecWAXPY(errorVec, -1.0, b, u)); 59 PetscCall(VecNorm(errorVec, NORM_2, &error)); 60 PetscCall(VecNorm(b, NORM_2, &norm)); 61 PetscCheck(error / norm <= 1000. * PETSC_MACHINE_EPSILON, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", (double)(error / norm)); 62 PetscCall(VecDestroy(&errorVec)); 63 PetscFunctionReturn(PETSC_SUCCESS); 64 } 65 66 PetscErrorCode ConstructProblem2(Mat A, Vec b) 67 { 68 PetscInt N = 10, constraintSize = 4; 69 PetscInt row; 70 71 PetscFunctionBeginUser; 72 PetscCall(VecSet(b, -3.0)); 73 for (row = 0; row < constraintSize; ++row) { 74 PetscScalar vals[2] = {1.0, 1.0}; 75 PetscInt cols[2]; 76 77 cols[0] = row; 78 cols[1] = row + N - constraintSize; 79 PetscCall(MatSetValues(A, 1, &row, 2, cols, vals, INSERT_VALUES)); 80 } 81 for (row = constraintSize; row < N - constraintSize; ++row) { 82 PetscScalar val = 1.0; 83 84 PetscCall(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; 89 90 PetscCall(MatSetValues(A, 1, &row, 1, &col, &val, INSERT_VALUES)); 91 } 92 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 93 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 94 PetscFunctionReturn(PETSC_SUCCESS); 95 } 96 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; 102 103 PetscFunctionBeginUser; 104 PetscCall(VecNorm(b, NORM_2, &norm)); 105 PetscCall(VecGetArrayRead(u, &uArray)); 106 PetscCall(VecGetArrayRead(b, &bArray)); 107 error = 0.0; 108 for (r = 0; r < constraintSize; ++r) error += PetscRealPart(PetscSqr(uArray[r] - bArray[r + N - constraintSize])); 109 110 PetscCheck(error / norm <= 10000 * PETSC_MACHINE_EPSILON, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", (double)(error / norm)); 111 error = 0.0; 112 for (r = constraintSize; r < N - constraintSize; ++r) error += PetscRealPart(PetscSqr(uArray[r] - bArray[r])); 113 114 PetscCheck(error / norm <= 10000 * PETSC_MACHINE_EPSILON, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", (double)(error / norm)); 115 error = 0.0; 116 for (r = N - constraintSize; r < N; ++r) error += PetscRealPart(PetscSqr(uArray[r] - (bArray[r - (N - constraintSize)] - bArray[r]))); 117 118 PetscCheck(error / norm <= 10000 * PETSC_MACHINE_EPSILON, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", (double)(error / norm)); 119 PetscCall(VecRestoreArrayRead(u, &uArray)); 120 PetscCall(VecRestoreArrayRead(b, &bArray)); 121 PetscFunctionReturn(PETSC_SUCCESS); 122 } 123 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; 131 132 PetscFunctionBeginUser; 133 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 134 comm = PETSC_COMM_WORLD; 135 PetscCall(PetscOptionsGetInt(NULL, NULL, "-problem", &problem, NULL)); 136 PetscCall(VecCreate(comm, &u)); 137 PetscCall(VecSetSizes(u, PETSC_DETERMINE, N)); 138 PetscCall(VecSetFromOptions(u)); 139 PetscCall(VecDuplicate(u, &r)); 140 PetscCall(VecDuplicate(u, &b)); 141 142 PetscCall(MatCreate(comm, &A)); 143 PetscCall(MatSetSizes(A, PETSC_DETERMINE, PETSC_DETERMINE, N, N)); 144 PetscCall(MatSetFromOptions(A)); 145 PetscCall(MatSeqAIJSetPreallocation(A, 5, NULL)); 146 J = A; 147 148 switch (problem) { 149 case 1: 150 PetscCall(ConstructProblem1(A, b)); 151 break; 152 case 2: 153 PetscCall(ConstructProblem2(A, b)); 154 break; 155 default: 156 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid problem number %" PetscInt_FMT, problem); 157 } 158 159 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 160 PetscCall(SNESSetJacobian(snes, A, J, ComputeJacobianLinear, NULL)); 161 PetscCall(SNESSetFunction(snes, r, ComputeFunctionLinear, A)); 162 PetscCall(SNESSetFromOptions(snes)); 163 164 PetscCall(SNESSolve(snes, b, u)); 165 PetscCall(VecView(u, NULL)); 166 167 switch (problem) { 168 case 1: 169 PetscCall(CheckProblem1(A, b, u)); 170 break; 171 case 2: 172 PetscCall(CheckProblem2(A, b, u)); 173 break; 174 default: 175 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid problem number %" PetscInt_FMT, problem); 176 } 177 178 if (A != J) PetscCall(MatDestroy(&A)); 179 PetscCall(MatDestroy(&J)); 180 PetscCall(VecDestroy(&u)); 181 PetscCall(VecDestroy(&r)); 182 PetscCall(VecDestroy(&b)); 183 PetscCall(SNESDestroy(&snes)); 184 PetscCall(PetscFinalize()); 185 return 0; 186 } 187 188 /*TEST 189 190 test: 191 args: -snes_monitor 192 193 test: 194 suffix: 2 195 args: -problem 2 -pc_type jacobi -snes_monitor 196 197 TEST*/ 198