1 2 static char help[] = "Newton's method to solve a two-variable system that comes from the Rosenbrock function.\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 #include <petscsnes.h> 14 15 extern PetscErrorCode FormJacobian1(SNES,Vec,Mat,Mat,void*); 16 extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*); 17 18 int main(int argc,char **argv) 19 { 20 SNES snes; /* nonlinear solver context */ 21 Vec x,r; /* solution, residual vectors */ 22 Mat J; /* Jacobian matrix */ 23 PetscInt its; 24 PetscScalar *xx; 25 SNESConvergedReason reason; 26 27 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 28 29 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 30 Create nonlinear solver context 31 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 32 PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes)); 33 34 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 35 Create matrix and vector data structures; set corresponding routines 36 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 37 /* 38 Create vectors for solution and nonlinear function 39 */ 40 PetscCall(VecCreate(PETSC_COMM_WORLD,&x)); 41 PetscCall(VecSetSizes(x,PETSC_DECIDE,2)); 42 PetscCall(VecSetFromOptions(x)); 43 PetscCall(VecDuplicate(x,&r)); 44 45 /* 46 Create Jacobian matrix data structure 47 */ 48 PetscCall(MatCreate(PETSC_COMM_WORLD,&J)); 49 PetscCall(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2)); 50 PetscCall(MatSetFromOptions(J)); 51 PetscCall(MatSetUp(J)); 52 53 /* 54 Set function evaluation routine and vector. 55 */ 56 PetscCall(SNESSetFunction(snes,r,FormFunction1,NULL)); 57 58 /* 59 Set Jacobian matrix data structure and Jacobian evaluation routine 60 */ 61 PetscCall(SNESSetJacobian(snes,J,J,FormJacobian1,NULL)); 62 63 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 64 Customize nonlinear solver; set runtime options 65 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 66 PetscCall(SNESSetFromOptions(snes)); 67 68 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 69 Evaluate initial guess; then solve nonlinear system 70 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 71 PetscCall(VecGetArray(x,&xx)); 72 xx[0] = -1.2; xx[1] = 1.0; 73 PetscCall(VecRestoreArray(x,&xx)); 74 75 /* 76 Note: The user should initialize the vector, x, with the initial guess 77 for the nonlinear solver prior to calling SNESSolve(). In particular, 78 to employ an initial guess of zero, the user should explicitly set 79 this vector to zero by calling VecSet(). 80 */ 81 82 PetscCall(SNESSolve(snes,NULL,x)); 83 PetscCall(VecView(x,PETSC_VIEWER_STDOUT_WORLD)); 84 PetscCall(SNESGetIterationNumber(snes,&its)); 85 PetscCall(SNESGetConvergedReason(snes,&reason)); 86 /* 87 Some systems computes a residual that is identically zero, thus converging 88 due to FNORM_ABS, others converge due to FNORM_RELATIVE. Here, we only 89 report the reason if the iteration did not converge so that the tests are 90 reproducible. 91 */ 92 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"%s number of SNES iterations = %" PetscInt_FMT "\n",reason>0 ? "CONVERGED" : SNESConvergedReasons[reason],its)); 93 94 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 95 Free work space. All PETSc objects should be destroyed when they 96 are no longer needed. 97 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 98 99 PetscCall(VecDestroy(&x)); PetscCall(VecDestroy(&r)); 100 PetscCall(MatDestroy(&J)); PetscCall(SNESDestroy(&snes)); 101 PetscCall(PetscFinalize()); 102 return 0; 103 } 104 /* ------------------------------------------------------------------- */ 105 /* 106 FormFunction1 - Evaluates nonlinear function, F(x). 107 108 Input Parameters: 109 . snes - the SNES context 110 . x - input vector 111 . ctx - optional user-defined context 112 113 Output Parameter: 114 . f - function vector 115 */ 116 PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx) 117 { 118 PetscScalar *ff; 119 const PetscScalar *xx; 120 121 /* 122 Get pointers to vector data. 123 - For default PETSc vectors, VecGetArray() returns a pointer to 124 the data array. Otherwise, the routine is implementation dependent. 125 - You MUST call VecRestoreArray() when you no longer need access to 126 the array. 127 */ 128 PetscCall(VecGetArrayRead(x,&xx)); 129 PetscCall(VecGetArray(f,&ff)); 130 131 /* Compute function */ 132 ff[0] = -2.0 + 2.0*xx[0] + 400.0*xx[0]*xx[0]*xx[0] - 400.0*xx[0]*xx[1]; 133 ff[1] = -200.0*xx[0]*xx[0] + 200.0*xx[1]; 134 135 /* Restore vectors */ 136 PetscCall(VecRestoreArrayRead(x,&xx)); 137 PetscCall(VecRestoreArray(f,&ff)); 138 return 0; 139 } 140 /* ------------------------------------------------------------------- */ 141 /* 142 FormJacobian1 - Evaluates Jacobian matrix. 143 144 Input Parameters: 145 . snes - the SNES context 146 . x - input vector 147 . dummy - optional user-defined context (not used here) 148 149 Output Parameters: 150 . jac - Jacobian matrix 151 . B - optionally different preconditioning matrix 152 . flag - flag indicating matrix structure 153 */ 154 PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat jac,Mat B,void *dummy) 155 { 156 const PetscScalar *xx; 157 PetscScalar A[4]; 158 PetscInt idx[2] = {0,1}; 159 160 /* 161 Get pointer to vector data 162 */ 163 PetscCall(VecGetArrayRead(x,&xx)); 164 165 /* 166 Compute Jacobian entries and insert into matrix. 167 - Since this is such a small problem, we set all entries for 168 the matrix at once. 169 */ 170 A[0] = 2.0 + 1200.0*xx[0]*xx[0] - 400.0*xx[1]; 171 A[1] = -400.0*xx[0]; 172 A[2] = -400.0*xx[0]; 173 A[3] = 200; 174 PetscCall(MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES)); 175 176 /* 177 Restore vector 178 */ 179 PetscCall(VecRestoreArrayRead(x,&xx)); 180 181 /* 182 Assemble matrix 183 */ 184 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 185 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 186 if (jac != B) { 187 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 188 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY)); 189 } 190 return 0; 191 } 192 193 /*TEST 194 195 test: 196 suffix: 1 197 args: -snes_monitor_short -snes_max_it 1000 198 requires: !single 199 200 test: 201 suffix: 2 202 args: -snes_monitor_short -snes_max_it 1000 -snes_type newtontrdc -snes_trdc_use_cauchy false 203 requires: !single 204 205 TEST*/ 206