1 2 static char help[] = "Tests SNESLinesearch handling of Inf/Nan.\n\n"; 3 4 /* 5 Include "petscsnes.h" so that we can use SNES solvers. Note that this 6 file automatically includes: 7 petscsys.h - base PETSc routines petscvec.h - vectors 8 petscmat.h - matrices 9 petscis.h - index sets petscksp.h - Krylov subspace methods 10 petscviewer.h - viewers petscpc.h - preconditioners 11 petscksp.h - linear solvers 12 */ 13 /*F 14 This examples solves 15 \begin{equation} 16 F\genfrac{(}{)}{0pt}{}{x_0}{x_1} = \genfrac{(}{)}{0pt}{}{\sin(3 x_0) + x_0}{x_1} 17 \end{equation} 18 F*/ 19 #include <petscsnes.h> 20 21 /* 22 User-defined routines 23 */ 24 extern PetscErrorCode FormJacobian2(SNES, Vec, Mat, Mat, void *); 25 extern PetscErrorCode FormFunction2(SNES, Vec, Vec, void *); 26 extern PetscErrorCode FormObjective(SNES, Vec, PetscReal *, void *); 27 28 /* 29 This is a very hacking way to trigger the objective function generating an infinity at a particular count to the call FormObjective(). 30 Different line searches evaluate the full step at different counts. For l2 it is the third call (infatcount == 2) while for bt it is the second call. 31 */ 32 PetscInt infatcount = 0; 33 34 int main(int argc, char **argv) 35 { 36 SNES snes; /* nonlinear solver context */ 37 KSP ksp; /* linear solver context */ 38 PC pc; /* preconditioner context */ 39 Vec x, r; /* solution, residual vectors */ 40 Mat J; /* Jacobian matrix */ 41 PetscInt its; 42 PetscMPIInt size; 43 PetscScalar *xx; 44 PetscBool flg; 45 char type[256]; 46 47 PetscFunctionBeginUser; 48 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 49 PetscCall(PetscOptionsGetString(NULL, NULL, "-snes_linesearch_type", type, sizeof(type), &flg)); 50 if (flg) { 51 PetscCall(PetscStrcmp(type, SNESLINESEARCHBT, &flg)); 52 if (flg) infatcount = 1; 53 PetscCall(PetscStrcmp(type, SNESLINESEARCHL2, &flg)); 54 if (flg) infatcount = 2; 55 } 56 57 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 58 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Example is only for sequential runs"); 59 60 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 61 Create nonlinear solver context 62 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 63 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 64 65 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 66 Create matrix and vector data structures; set corresponding routines 67 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 68 /* 69 Create vectors for solution and nonlinear function 70 */ 71 PetscCall(VecCreate(PETSC_COMM_WORLD, &x)); 72 PetscCall(VecSetSizes(x, PETSC_DECIDE, 2)); 73 PetscCall(VecSetFromOptions(x)); 74 PetscCall(VecDuplicate(x, &r)); 75 76 /* 77 Create Jacobian matrix data structure 78 */ 79 PetscCall(MatCreate(PETSC_COMM_WORLD, &J)); 80 PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, 2, 2)); 81 PetscCall(MatSetFromOptions(J)); 82 PetscCall(MatSetUp(J)); 83 84 PetscCall(SNESSetFunction(snes, r, FormFunction2, NULL)); 85 PetscCall(SNESSetObjective(snes, FormObjective, NULL)); 86 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian2, NULL)); 87 88 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 89 Customize nonlinear solver; set runtime options 90 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 91 /* 92 Set linear solver defaults for this problem. By extracting the 93 KSP and PC contexts from the SNES context, we can then 94 directly call any KSP and PC routines to set various options. 95 */ 96 PetscCall(SNESGetKSP(snes, &ksp)); 97 PetscCall(KSPGetPC(ksp, &pc)); 98 PetscCall(PCSetType(pc, PCNONE)); 99 PetscCall(KSPSetTolerances(ksp, 1.e-4, PETSC_DEFAULT, PETSC_DEFAULT, 20)); 100 101 /* 102 Set SNES/KSP/KSP/PC runtime options, e.g., 103 -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc> 104 These options will override those specified above as long as 105 SNESSetFromOptions() is called _after_ any other customization 106 routines. 107 */ 108 PetscCall(SNESSetFromOptions(snes)); 109 110 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 111 Evaluate initial guess; then solve nonlinear system 112 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 113 PetscCall(VecGetArray(x, &xx)); 114 xx[0] = 2.0; 115 xx[1] = 3.0; 116 PetscCall(VecRestoreArray(x, &xx)); 117 118 /* 119 Note: The user should initialize the vector, x, with the initial guess 120 for the nonlinear solver prior to calling SNESSolve(). In particular, 121 to employ an initial guess of zero, the user should explicitly set 122 this vector to zero by calling VecSet(). 123 */ 124 125 PetscCall(SNESSolve(snes, NULL, x)); 126 PetscCall(SNESGetIterationNumber(snes, &its)); 127 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its)); 128 129 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 130 Free work space. All PETSc objects should be destroyed when they 131 are no longer needed. 132 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 133 134 PetscCall(VecDestroy(&x)); 135 PetscCall(VecDestroy(&r)); 136 PetscCall(MatDestroy(&J)); 137 PetscCall(SNESDestroy(&snes)); 138 PetscCall(PetscFinalize()); 139 return 0; 140 } 141 142 PetscErrorCode FormObjective(SNES snes, Vec x, PetscReal *f, void *dummy) 143 { 144 Vec F; 145 static PetscInt cnt = 0; 146 const PetscReal one = 1.0, zero = 0.0; 147 PetscReal inf; 148 149 PetscFunctionBeginUser; 150 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 151 inf = one / zero; 152 PetscCall(PetscFPTrapPop()); 153 if (cnt++ == infatcount) *f = inf; 154 else { 155 PetscCall(VecDuplicate(x, &F)); 156 PetscCall(FormFunction2(snes, x, F, dummy)); 157 PetscCall(VecNorm(F, NORM_2, f)); 158 PetscCall(VecDestroy(&F)); 159 *f = (*f) * (*f); 160 } 161 PetscFunctionReturn(PETSC_SUCCESS); 162 } 163 164 PetscErrorCode FormFunction2(SNES snes, Vec x, Vec f, void *dummy) 165 { 166 const PetscScalar *xx; 167 PetscScalar *ff; 168 169 PetscFunctionBeginUser; 170 /* 171 Get pointers to vector data. 172 - For default PETSc vectors, VecGetArray() returns a pointer to 173 the data array. Otherwise, the routine is implementation dependent. 174 - You MUST call VecRestoreArray() when you no longer need access to 175 the array. 176 */ 177 PetscCall(VecGetArrayRead(x, &xx)); 178 PetscCall(VecGetArray(f, &ff)); 179 180 /* 181 Compute function 182 */ 183 ff[0] = PetscSinScalar(3.0 * xx[0]) + xx[0]; 184 ff[1] = xx[1]; 185 186 /* 187 Restore vectors 188 */ 189 PetscCall(VecRestoreArrayRead(x, &xx)); 190 PetscCall(VecRestoreArray(f, &ff)); 191 PetscFunctionReturn(PETSC_SUCCESS); 192 } 193 194 PetscErrorCode FormJacobian2(SNES snes, Vec x, Mat jac, Mat B, void *dummy) 195 { 196 const PetscScalar *xx; 197 PetscScalar A[4]; 198 PetscInt idx[2] = {0, 1}; 199 200 PetscFunctionBeginUser; 201 /* 202 Get pointer to vector data 203 */ 204 PetscCall(VecGetArrayRead(x, &xx)); 205 206 /* 207 Compute Jacobian entries and insert into matrix. 208 - Since this is such a small problem, we set all entries for 209 the matrix at once. 210 */ 211 A[0] = 3.0 * PetscCosScalar(3.0 * xx[0]) + 1.0; 212 A[1] = 0.0; 213 A[2] = 0.0; 214 A[3] = 1.0; 215 PetscCall(MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES)); 216 217 /* 218 Restore vector 219 */ 220 PetscCall(VecRestoreArrayRead(x, &xx)); 221 222 /* 223 Assemble matrix 224 */ 225 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 226 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 227 if (jac != B) { 228 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 229 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 230 } 231 PetscFunctionReturn(PETSC_SUCCESS); 232 } 233 234 /*TEST 235 236 test: 237 args: -snes_converged_reason -snes_linesearch_monitor -snes_linesearch_type l2 238 filter: grep Inf 239 240 TEST*/ 241