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 { 21 Mat A = (Mat) ctx; 22 23 PetscFunctionBeginUser; 24 PetscCall(MatMult(A, x, f)); 25 PetscFunctionReturn(0); 26 } 27 28 PetscErrorCode ComputeJacobianLinear(SNES snes, Vec x, Mat A, Mat J, void *ctx) 29 { 30 PetscFunctionBeginUser; 31 PetscFunctionReturn(0); 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(0); 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", error/norm); 62 PetscCall(VecDestroy(&errorVec)); 63 PetscFunctionReturn(0); 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; cols[1] = row + N - constraintSize; 78 PetscCall(MatSetValues(A, 1, &row, 2, cols, vals, INSERT_VALUES)); 79 } 80 for (row = constraintSize; row < N - constraintSize; ++row) { 81 PetscScalar val = 1.0; 82 83 PetscCall(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; 88 89 PetscCall(MatSetValues(A, 1, &row, 1, &col, &val, INSERT_VALUES)); 90 } 91 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 92 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 93 PetscFunctionReturn(0); 94 } 95 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; 101 102 PetscFunctionBeginUser; 103 PetscCall(VecNorm(b, NORM_2, &norm)); 104 PetscCall(VecGetArrayRead(u, &uArray)); 105 PetscCall(VecGetArrayRead(b, &bArray)); 106 error = 0.0; 107 for (r = 0; r < constraintSize; ++r) error += PetscRealPart(PetscSqr(uArray[r] - bArray[r + N-constraintSize])); 108 109 PetscCheck(error/norm <= 10000*PETSC_MACHINE_EPSILON,PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", error/norm); 110 error = 0.0; 111 for (r = constraintSize; r < N - constraintSize; ++r) error += PetscRealPart(PetscSqr(uArray[r] - bArray[r])); 112 113 PetscCheck(error/norm <= 10000*PETSC_MACHINE_EPSILON,PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", error/norm); 114 error = 0.0; 115 for (r = N - constraintSize; r < N; ++r) error += PetscRealPart(PetscSqr(uArray[r] - (bArray[r - (N-constraintSize)] - bArray[r]))); 116 117 PetscCheck(error/norm <= 10000*PETSC_MACHINE_EPSILON,PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", error/norm); 118 PetscCall(VecRestoreArrayRead(u, &uArray)); 119 PetscCall(VecRestoreArrayRead(b, &bArray)); 120 PetscFunctionReturn(0); 121 } 122 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; 130 131 PetscCall(PetscInitialize(&argc, &argv, NULL,help)); 132 comm = PETSC_COMM_WORLD; 133 PetscCall(PetscOptionsGetInt(NULL,NULL, "-problem", &problem, NULL)); 134 PetscCall(VecCreate(comm, &u)); 135 PetscCall(VecSetSizes(u, PETSC_DETERMINE, N)); 136 PetscCall(VecSetFromOptions(u)); 137 PetscCall(VecDuplicate(u, &r)); 138 PetscCall(VecDuplicate(u, &b)); 139 140 PetscCall(MatCreate(comm, &A)); 141 PetscCall(MatSetSizes(A, PETSC_DETERMINE, PETSC_DETERMINE, N, N)); 142 PetscCall(MatSetFromOptions(A)); 143 PetscCall(MatSeqAIJSetPreallocation(A, 5, NULL)); 144 J = A; 145 146 switch (problem) { 147 case 1: 148 PetscCall(ConstructProblem1(A, b)); 149 break; 150 case 2: 151 PetscCall(ConstructProblem2(A, b)); 152 break; 153 default: 154 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid problem number %d", problem); 155 } 156 157 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 158 PetscCall(SNESSetJacobian(snes, A, J, ComputeJacobianLinear, NULL)); 159 PetscCall(SNESSetFunction(snes, r, ComputeFunctionLinear, A)); 160 PetscCall(SNESSetFromOptions(snes)); 161 162 PetscCall(SNESSolve(snes, b, u)); 163 PetscCall(VecView(u, NULL)); 164 165 switch (problem) { 166 case 1: 167 PetscCall(CheckProblem1(A, b, u)); 168 break; 169 case 2: 170 PetscCall(CheckProblem2(A, b, u)); 171 break; 172 default: 173 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid problem number %d", problem); 174 } 175 176 if (A != J) { 177 PetscCall(MatDestroy(&A)); 178 } 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