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