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