1 static const char help[] = "Attempts to solve for root of a function with multiple local minima.\n\ 2 With the proper initial guess, a backtracking line-search fails because Newton's method gets\n\ 3 stuck in a local minimum. However, a critical point line-search or Newton's method without a\n\ 4 line search succeeds.\n"; 5 6 /* Solve 1D problem f(x) = 8 * exp(-4 * (x - 2)^2) * (x - 2) + 2 * x = 0 7 8 This problem is based on the example given here: https://scicomp.stackexchange.com/a/2446/24756 9 Originally an optimization problem to find the minimum of the function 10 11 g(x) = x^2 - exp(-4 * (x - 2)^2) 12 13 it has been reformulated to solve dg(x)/dx = f(x) = 0. The reformulated problem has several local 14 minima that can cause problems for some global Newton root-finding methods. In this particular 15 example, an initial guess of x0 = 2.5 generates an initial search direction (-df/dx is positive) 16 away from the root and towards a local minimum in which a back-tracking line search gets trapped. 17 However, omitting a line-search or using a critical point line search, the solve is successful. 18 19 The test outputs the final result for x and f(x). 20 21 Example usage: 22 23 Get help: 24 ./ex99 -help 25 26 Monitor run (with default back-tracking line search; solve fails): 27 ./ex99 -snes_converged_reason -snes_monitor -snes_linesearch_monitor -ksp_converged_reason -ksp_monitor 28 29 Run without a line search; solve succeeds: 30 ./ex99 -snes_linesearch_type basic 31 32 Run with a critical point line search; solve succeeds: 33 ./ex99 -snes_linesearch_type cp 34 */ 35 36 #include <math.h> 37 #include <petscsnes.h> 38 39 extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*); 40 extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*); 41 42 int main(int argc,char **argv) 43 { 44 SNES snes; /* nonlinear solver context */ 45 KSP ksp; /* linear solver context */ 46 PC pc; /* preconditioner context */ 47 Vec x,r; /* solution, residual vectors */ 48 Mat J; /* Jacobian matrix */ 49 PetscErrorCode ierr; 50 PetscMPIInt size; 51 52 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 53 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 54 if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Example is only for sequential runs"); 55 56 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 57 Create nonlinear solver context 58 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 59 ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); 60 61 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 62 Create matrix and vector data structures; set corresponding routines 63 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 64 /* 65 Create vectors for solution and nonlinear function 66 */ 67 ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr); 68 ierr = VecSetSizes(x,PETSC_DECIDE,1);CHKERRQ(ierr); 69 ierr = VecSetFromOptions(x);CHKERRQ(ierr); 70 ierr = VecDuplicate(x,&r);CHKERRQ(ierr); 71 72 /* 73 Create Jacobian matrix data structure 74 */ 75 ierr = MatCreate(PETSC_COMM_WORLD,&J);CHKERRQ(ierr); 76 ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,1,1);CHKERRQ(ierr); 77 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 78 ierr = MatSetUp(J);CHKERRQ(ierr); 79 80 ierr = SNESSetFunction(snes,r,FormFunction,NULL);CHKERRQ(ierr); 81 ierr = SNESSetJacobian(snes,J,J,FormJacobian,NULL);CHKERRQ(ierr); 82 83 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 84 Customize nonlinear solver; set runtime options 85 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 86 /* 87 Set linear solver defaults for this problem. By extracting the 88 KSP and PC contexts from the SNES context, we can then 89 directly call any KSP and PC routines to set various options. 90 */ 91 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 92 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 93 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 94 ierr = KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);CHKERRQ(ierr); 95 96 /* 97 Set SNES/KSP/KSP/PC runtime options, e.g., 98 -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc> 99 These options will override those specified above as long as 100 SNESSetFromOptions() is called _after_ any other customization 101 routines. 102 */ 103 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 104 105 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 106 Evaluate initial guess; then solve nonlinear system 107 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 108 ierr = VecSet(x,2.5);CHKERRQ(ierr); 109 110 ierr = SNESSolve(snes,NULL,x);CHKERRQ(ierr); 111 112 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 113 Output x and f(x) 114 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 115 ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 116 ierr = VecView(r,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 117 118 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 119 Free work space. All PETSc objects should be destroyed when they 120 are no longer needed. 121 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 122 123 ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&r);CHKERRQ(ierr); 124 ierr = MatDestroy(&J);CHKERRQ(ierr); ierr = SNESDestroy(&snes);CHKERRQ(ierr); 125 ierr = PetscFinalize(); 126 return ierr; 127 } 128 129 PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx) 130 { 131 PetscErrorCode ierr; 132 const PetscScalar *xx; 133 PetscScalar *ff; 134 135 /* 136 Get pointers to vector data. 137 - For default PETSc vectors, VecGetArray() returns a pointer to 138 the data array. Otherwise, the routine is implementation dependent. 139 - You MUST call VecRestoreArray() when you no longer need access to 140 the array. 141 */ 142 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 143 ierr = VecGetArray(f,&ff);CHKERRQ(ierr); 144 145 /* Compute function */ 146 ff[0] = 8. * PetscExpScalar(-4. * (xx[0] - 2.) * (xx[0] - 2.)) * (xx[0] - 2.) + 2. * xx[0]; 147 148 /* Restore vectors */ 149 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 150 ierr = VecRestoreArray(f,&ff);CHKERRQ(ierr); 151 return 0; 152 } 153 154 PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy) 155 { 156 const PetscScalar *xx; 157 PetscScalar A[1]; 158 PetscErrorCode ierr; 159 PetscInt idx[1] = {0}; 160 161 /* 162 Get pointer to vector data 163 */ 164 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 165 166 /* 167 Compute Jacobian entries and insert into matrix. 168 - Since this is such a small problem, we set all entries for 169 the matrix at once. 170 */ 171 A[0] = 8. * ((xx[0] - 2.) * (PetscExpScalar(-4. * (xx[0] - 2.) * (xx[0] - 2.)) * -8. * (xx[0] - 2.)) 172 + PetscExpScalar(-4. * (xx[0] - 2.) * (xx[0] - 2.))) 173 + 2.; 174 175 ierr = MatSetValues(B,1,idx,1,idx,A,INSERT_VALUES);CHKERRQ(ierr); 176 177 /* 178 Restore vector 179 */ 180 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 181 182 /* 183 Assemble matrix 184 */ 185 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 186 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 187 if (jac != B) { 188 ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 189 ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 190 } 191 return 0; 192 } 193 194 /*TEST 195 196 test: 197 suffix: 1 198 args: -snes_linesearch_type cp 199 test: 200 suffix: 2 201 args: -snes_linesearch_type basic 202 test: 203 suffix: 3 204 205 TEST*/ 206