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