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 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 48 PetscCall(PetscOptionsGetString(NULL,NULL,"-snes_linesearch_type",type,sizeof(type),&flg)); 49 if (flg) { 50 PetscCall(PetscStrcmp(type,SNESLINESEARCHBT,&flg)); 51 if (flg) infatcount = 1; 52 PetscCall(PetscStrcmp(type,SNESLINESEARCHL2,&flg)); 53 if (flg) infatcount = 2; 54 } 55 56 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 57 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Example is only for sequential runs"); 58 59 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 60 Create nonlinear solver context 61 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 62 PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes)); 63 64 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 65 Create matrix and vector data structures; set corresponding routines 66 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 67 /* 68 Create vectors for solution and nonlinear function 69 */ 70 PetscCall(VecCreate(PETSC_COMM_WORLD,&x)); 71 PetscCall(VecSetSizes(x,PETSC_DECIDE,2)); 72 PetscCall(VecSetFromOptions(x)); 73 PetscCall(VecDuplicate(x,&r)); 74 75 /* 76 Create Jacobian matrix data structure 77 */ 78 PetscCall(MatCreate(PETSC_COMM_WORLD,&J)); 79 PetscCall(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2)); 80 PetscCall(MatSetFromOptions(J)); 81 PetscCall(MatSetUp(J)); 82 83 PetscCall(SNESSetFunction(snes,r,FormFunction2,NULL)); 84 PetscCall(SNESSetObjective(snes,FormObjective,NULL)); 85 PetscCall(SNESSetJacobian(snes,J,J,FormJacobian2,NULL)); 86 87 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 88 Customize nonlinear solver; set runtime options 89 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 90 /* 91 Set linear solver defaults for this problem. By extracting the 92 KSP and PC contexts from the SNES context, we can then 93 directly call any KSP and PC routines to set various options. 94 */ 95 PetscCall(SNESGetKSP(snes,&ksp)); 96 PetscCall(KSPGetPC(ksp,&pc)); 97 PetscCall(PCSetType(pc,PCNONE)); 98 PetscCall(KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20)); 99 100 /* 101 Set SNES/KSP/KSP/PC runtime options, e.g., 102 -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc> 103 These options will override those specified above as long as 104 SNESSetFromOptions() is called _after_ any other customization 105 routines. 106 */ 107 PetscCall(SNESSetFromOptions(snes)); 108 109 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 110 Evaluate initial guess; then solve nonlinear system 111 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 112 PetscCall(VecGetArray(x,&xx)); 113 xx[0] = 2.0; xx[1] = 3.0; 114 PetscCall(VecRestoreArray(x,&xx)); 115 116 /* 117 Note: The user should initialize the vector, x, with the initial guess 118 for the nonlinear solver prior to calling SNESSolve(). In particular, 119 to employ an initial guess of zero, the user should explicitly set 120 this vector to zero by calling VecSet(). 121 */ 122 123 PetscCall(SNESSolve(snes,NULL,x)); 124 PetscCall(SNESGetIterationNumber(snes,&its)); 125 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %" PetscInt_FMT "\n",its)); 126 127 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 128 Free work space. All PETSc objects should be destroyed when they 129 are no longer needed. 130 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 131 132 PetscCall(VecDestroy(&x)); PetscCall(VecDestroy(&r)); 133 PetscCall(MatDestroy(&J)); PetscCall(SNESDestroy(&snes)); 134 PetscCall(PetscFinalize()); 135 return 0; 136 } 137 138 PetscErrorCode FormObjective(SNES snes,Vec x,PetscReal *f,void *dummy) 139 { 140 Vec F; 141 static PetscInt cnt = 0; 142 143 if (cnt++ == infatcount) *f = INFINITY; 144 else { 145 PetscCall(VecDuplicate(x,&F)); 146 PetscCall(FormFunction2(snes,x,F,dummy)); 147 PetscCall(VecNorm(F,NORM_2,f)); 148 PetscCall(VecDestroy(&F)); 149 *f = (*f)*(*f); 150 } 151 return 0; 152 } 153 154 /* ------------------------------------------------------------------- */ 155 PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy) 156 { 157 const PetscScalar *xx; 158 PetscScalar *ff; 159 160 /* 161 Get pointers to vector data. 162 - For default PETSc vectors, VecGetArray() returns a pointer to 163 the data array. Otherwise, the routine is implementation dependent. 164 - You MUST call VecRestoreArray() when you no longer need access to 165 the array. 166 */ 167 PetscCall(VecGetArrayRead(x,&xx)); 168 PetscCall(VecGetArray(f,&ff)); 169 170 /* 171 Compute function 172 */ 173 ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0]; 174 ff[1] = xx[1]; 175 176 /* 177 Restore vectors 178 */ 179 PetscCall(VecRestoreArrayRead(x,&xx)); 180 PetscCall(VecRestoreArray(f,&ff)); 181 return 0; 182 } 183 /* ------------------------------------------------------------------- */ 184 PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat jac,Mat B,void *dummy) 185 { 186 const PetscScalar *xx; 187 PetscScalar A[4]; 188 PetscInt idx[2] = {0,1}; 189 190 /* 191 Get pointer to vector data 192 */ 193 PetscCall(VecGetArrayRead(x,&xx)); 194 195 /* 196 Compute Jacobian entries and insert into matrix. 197 - Since this is such a small problem, we set all entries for 198 the matrix at once. 199 */ 200 A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0; 201 A[2] = 0.0; A[3] = 1.0; 202 PetscCall(MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES)); 203 204 /* 205 Restore vector 206 */ 207 PetscCall(VecRestoreArrayRead(x,&xx)); 208 209 /* 210 Assemble matrix 211 */ 212 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 213 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 214 if (jac != B) { 215 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 216 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY)); 217 } 218 return 0; 219 } 220 221 /*TEST 222 223 build: 224 requires: infinity 225 226 test: 227 args: -snes_converged_reason -snes_linesearch_monitor -snes_linesearch_type l2 228 filter: grep Inf 229 230 TEST*/ 231