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, inf = one / zero; 147 148 PetscFunctionBeginUser; 149 if (cnt++ == infatcount) *f = inf; 150 else { 151 PetscCall(VecDuplicate(x, &F)); 152 PetscCall(FormFunction2(snes, x, F, dummy)); 153 PetscCall(VecNorm(F, NORM_2, f)); 154 PetscCall(VecDestroy(&F)); 155 *f = (*f) * (*f); 156 } 157 PetscFunctionReturn(PETSC_SUCCESS); 158 } 159 160 PetscErrorCode FormFunction2(SNES snes, Vec x, Vec f, void *dummy) 161 { 162 const PetscScalar *xx; 163 PetscScalar *ff; 164 165 PetscFunctionBeginUser; 166 /* 167 Get pointers to vector data. 168 - For default PETSc vectors, VecGetArray() returns a pointer to 169 the data array. Otherwise, the routine is implementation dependent. 170 - You MUST call VecRestoreArray() when you no longer need access to 171 the array. 172 */ 173 PetscCall(VecGetArrayRead(x, &xx)); 174 PetscCall(VecGetArray(f, &ff)); 175 176 /* 177 Compute function 178 */ 179 ff[0] = PetscSinScalar(3.0 * xx[0]) + xx[0]; 180 ff[1] = xx[1]; 181 182 /* 183 Restore vectors 184 */ 185 PetscCall(VecRestoreArrayRead(x, &xx)); 186 PetscCall(VecRestoreArray(f, &ff)); 187 PetscFunctionReturn(PETSC_SUCCESS); 188 } 189 190 PetscErrorCode FormJacobian2(SNES snes, Vec x, Mat jac, Mat B, void *dummy) 191 { 192 const PetscScalar *xx; 193 PetscScalar A[4]; 194 PetscInt idx[2] = {0, 1}; 195 196 PetscFunctionBeginUser; 197 /* 198 Get pointer to vector data 199 */ 200 PetscCall(VecGetArrayRead(x, &xx)); 201 202 /* 203 Compute Jacobian entries and insert into matrix. 204 - Since this is such a small problem, we set all entries for 205 the matrix at once. 206 */ 207 A[0] = 3.0 * PetscCosScalar(3.0 * xx[0]) + 1.0; 208 A[1] = 0.0; 209 A[2] = 0.0; 210 A[3] = 1.0; 211 PetscCall(MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES)); 212 213 /* 214 Restore vector 215 */ 216 PetscCall(VecRestoreArrayRead(x, &xx)); 217 218 /* 219 Assemble matrix 220 */ 221 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 222 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 223 if (jac != B) { 224 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 225 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 226 } 227 PetscFunctionReturn(PETSC_SUCCESS); 228 } 229 230 /*TEST 231 232 test: 233 args: -snes_converged_reason -snes_linesearch_monitor -snes_linesearch_type l2 234 filter: grep Inf 235 236 TEST*/ 237