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