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