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