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; xx[1] = 3.0; 115 PetscCall(VecRestoreArray(x,&xx)); 116 117 /* 118 Note: The user should initialize the vector, x, with the initial guess 119 for the nonlinear solver prior to calling SNESSolve(). In particular, 120 to employ an initial guess of zero, the user should explicitly set 121 this vector to zero by calling VecSet(). 122 */ 123 124 PetscCall(SNESSolve(snes,NULL,x)); 125 PetscCall(SNESGetIterationNumber(snes,&its)); 126 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %" PetscInt_FMT "\n",its)); 127 128 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 129 Free work space. All PETSc objects should be destroyed when they 130 are no longer needed. 131 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 132 133 PetscCall(VecDestroy(&x)); PetscCall(VecDestroy(&r)); 134 PetscCall(MatDestroy(&J)); PetscCall(SNESDestroy(&snes)); 135 PetscCall(PetscFinalize()); 136 return 0; 137 } 138 139 PetscErrorCode FormObjective(SNES snes,Vec x,PetscReal *f,void *dummy) 140 { 141 Vec F; 142 static PetscInt cnt = 0; 143 144 if (cnt++ == infatcount) *f = INFINITY; 145 else { 146 PetscCall(VecDuplicate(x,&F)); 147 PetscCall(FormFunction2(snes,x,F,dummy)); 148 PetscCall(VecNorm(F,NORM_2,f)); 149 PetscCall(VecDestroy(&F)); 150 *f = (*f)*(*f); 151 } 152 return 0; 153 } 154 155 /* ------------------------------------------------------------------- */ 156 PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy) 157 { 158 const PetscScalar *xx; 159 PetscScalar *ff; 160 161 /* 162 Get pointers to vector data. 163 - For default PETSc vectors, VecGetArray() returns a pointer to 164 the data array. Otherwise, the routine is implementation dependent. 165 - You MUST call VecRestoreArray() when you no longer need access to 166 the array. 167 */ 168 PetscCall(VecGetArrayRead(x,&xx)); 169 PetscCall(VecGetArray(f,&ff)); 170 171 /* 172 Compute function 173 */ 174 ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0]; 175 ff[1] = xx[1]; 176 177 /* 178 Restore vectors 179 */ 180 PetscCall(VecRestoreArrayRead(x,&xx)); 181 PetscCall(VecRestoreArray(f,&ff)); 182 return 0; 183 } 184 /* ------------------------------------------------------------------- */ 185 PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat jac,Mat B,void *dummy) 186 { 187 const PetscScalar *xx; 188 PetscScalar A[4]; 189 PetscInt idx[2] = {0,1}; 190 191 /* 192 Get pointer to vector data 193 */ 194 PetscCall(VecGetArrayRead(x,&xx)); 195 196 /* 197 Compute Jacobian entries and insert into matrix. 198 - Since this is such a small problem, we set all entries for 199 the matrix at once. 200 */ 201 A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0; 202 A[2] = 0.0; A[3] = 1.0; 203 PetscCall(MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES)); 204 205 /* 206 Restore vector 207 */ 208 PetscCall(VecRestoreArrayRead(x,&xx)); 209 210 /* 211 Assemble matrix 212 */ 213 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 214 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 215 if (jac != B) { 216 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 217 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY)); 218 } 219 return 0; 220 } 221 222 /*TEST 223 224 build: 225 requires: infinity 226 227 test: 228 args: -snes_converged_reason -snes_linesearch_monitor -snes_linesearch_type l2 229 filter: grep Inf 230 231 TEST*/ 232