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