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 PetscMPIInt size; 50 51 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 52 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 53 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Example is only for sequential runs"); 54 55 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 56 Create nonlinear solver context 57 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 58 PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes)); 59 60 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 61 Create matrix and vector data structures; set corresponding routines 62 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 63 /* 64 Create vectors for solution and nonlinear function 65 */ 66 PetscCall(VecCreate(PETSC_COMM_WORLD,&x)); 67 PetscCall(VecSetSizes(x,PETSC_DECIDE,1)); 68 PetscCall(VecSetFromOptions(x)); 69 PetscCall(VecDuplicate(x,&r)); 70 71 /* 72 Create Jacobian matrix data structure 73 */ 74 PetscCall(MatCreate(PETSC_COMM_WORLD,&J)); 75 PetscCall(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,1,1)); 76 PetscCall(MatSetFromOptions(J)); 77 PetscCall(MatSetUp(J)); 78 79 PetscCall(SNESSetFunction(snes,r,FormFunction,NULL)); 80 PetscCall(SNESSetJacobian(snes,J,J,FormJacobian,NULL)); 81 82 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 83 Customize nonlinear solver; set runtime options 84 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 85 /* 86 Set linear solver defaults for this problem. By extracting the 87 KSP and PC contexts from the SNES context, we can then 88 directly call any KSP and PC routines to set various options. 89 */ 90 PetscCall(SNESGetKSP(snes,&ksp)); 91 PetscCall(KSPGetPC(ksp,&pc)); 92 PetscCall(PCSetType(pc,PCNONE)); 93 PetscCall(KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20)); 94 95 /* 96 Set SNES/KSP/KSP/PC runtime options, e.g., 97 -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc> 98 These options will override those specified above as long as 99 SNESSetFromOptions() is called _after_ any other customization 100 routines. 101 */ 102 PetscCall(SNESSetFromOptions(snes)); 103 104 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 105 Evaluate initial guess; then solve nonlinear system 106 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 107 PetscCall(VecSet(x,2.5)); 108 109 PetscCall(SNESSolve(snes,NULL,x)); 110 111 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 112 Output x and f(x) 113 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 114 PetscCall(VecView(x,PETSC_VIEWER_STDOUT_WORLD)); 115 PetscCall(VecView(r,PETSC_VIEWER_STDOUT_WORLD)); 116 117 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 118 Free work space. All PETSc objects should be destroyed when they 119 are no longer needed. 120 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 121 122 PetscCall(VecDestroy(&x)); PetscCall(VecDestroy(&r)); 123 PetscCall(MatDestroy(&J)); PetscCall(SNESDestroy(&snes)); 124 PetscCall(PetscFinalize()); 125 return 0; 126 } 127 128 PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx) 129 { 130 const PetscScalar *xx; 131 PetscScalar *ff; 132 133 /* 134 Get pointers to vector data. 135 - For default PETSc vectors, VecGetArray() returns a pointer to 136 the data array. Otherwise, the routine is implementation dependent. 137 - You MUST call VecRestoreArray() when you no longer need access to 138 the array. 139 */ 140 PetscCall(VecGetArrayRead(x,&xx)); 141 PetscCall(VecGetArray(f,&ff)); 142 143 /* Compute function */ 144 ff[0] = 8. * PetscExpScalar(-4. * (xx[0] - 2.) * (xx[0] - 2.)) * (xx[0] - 2.) + 2. * xx[0]; 145 146 /* Restore vectors */ 147 PetscCall(VecRestoreArrayRead(x,&xx)); 148 PetscCall(VecRestoreArray(f,&ff)); 149 return 0; 150 } 151 152 PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy) 153 { 154 const PetscScalar *xx; 155 PetscScalar A[1]; 156 PetscInt idx[1] = {0}; 157 158 /* 159 Get pointer to vector data 160 */ 161 PetscCall(VecGetArrayRead(x,&xx)); 162 163 /* 164 Compute Jacobian entries and insert into matrix. 165 - Since this is such a small problem, we set all entries for 166 the matrix at once. 167 */ 168 A[0] = 8. * ((xx[0] - 2.) * (PetscExpScalar(-4. * (xx[0] - 2.) * (xx[0] - 2.)) * -8. * (xx[0] - 2.)) 169 + PetscExpScalar(-4. * (xx[0] - 2.) * (xx[0] - 2.))) 170 + 2.; 171 172 PetscCall(MatSetValues(B,1,idx,1,idx,A,INSERT_VALUES)); 173 174 /* 175 Restore vector 176 */ 177 PetscCall(VecRestoreArrayRead(x,&xx)); 178 179 /* 180 Assemble matrix 181 */ 182 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 183 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 184 if (jac != B) { 185 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 186 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY)); 187 } 188 return 0; 189 } 190 191 /*TEST 192 193 test: 194 suffix: 1 195 args: -snes_linesearch_type cp 196 test: 197 suffix: 2 198 args: -snes_linesearch_type basic 199 test: 200 suffix: 3 201 test: 202 suffix: 4 203 args: -snes_type newtontrdc 204 test: 205 suffix: 5 206 args: -snes_type newtontrdc -snes_trdc_use_cauchy false 207 208 TEST*/ 209