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 ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); 38 39 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 40 Create matrix and vector data structures; set corresponding routines 41 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 42 /* 43 Create vectors for solution and nonlinear function 44 */ 45 ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr); 46 ierr = VecSetSizes(x,PETSC_DECIDE,2);CHKERRQ(ierr); 47 ierr = VecSetFromOptions(x);CHKERRQ(ierr); 48 ierr = VecDuplicate(x,&r);CHKERRQ(ierr); 49 50 /* 51 Create Jacobian matrix data structure 52 */ 53 ierr = MatCreate(PETSC_COMM_WORLD,&J);CHKERRQ(ierr); 54 ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);CHKERRQ(ierr); 55 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 56 ierr = MatSetUp(J);CHKERRQ(ierr); 57 58 /* 59 Set function evaluation routine and vector. 60 */ 61 ierr = SNESSetFunction(snes,r,FormFunction1,NULL);CHKERRQ(ierr); 62 63 /* 64 Set Jacobian matrix data structure and Jacobian evaluation routine 65 */ 66 ierr = SNESSetJacobian(snes,J,J,FormJacobian1,NULL);CHKERRQ(ierr); 67 68 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 69 Customize nonlinear solver; set runtime options 70 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 71 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 72 73 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 74 Evaluate initial guess; then solve nonlinear system 75 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 76 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 77 xx[0] = -1.2; xx[1] = 1.0; 78 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 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 ierr = SNESSolve(snes,NULL,x);CHKERRQ(ierr); 88 ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 89 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 90 ierr = SNESGetConvergedReason(snes,&reason);CHKERRQ(ierr); 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 ierr = PetscPrintf(PETSC_COMM_WORLD,"%s number of SNES iterations = %D\n",reason>0 ? "CONVERGED" : SNESConvergedReasons[reason],its);CHKERRQ(ierr); 98 99 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 100 Free work space. All PETSc objects should be destroyed when they 101 are no longer needed. 102 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 103 104 ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&r);CHKERRQ(ierr); 105 ierr = MatDestroy(&J);CHKERRQ(ierr); ierr = SNESDestroy(&snes);CHKERRQ(ierr); 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 PetscErrorCode ierr; 124 PetscScalar *ff; 125 const PetscScalar *xx; 126 127 /* 128 Get pointers to vector data. 129 - For default PETSc vectors, VecGetArray() returns a pointer to 130 the data array. Otherwise, the routine is implementation dependent. 131 - You MUST call VecRestoreArray() when you no longer need access to 132 the array. 133 */ 134 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 135 ierr = VecGetArray(f,&ff);CHKERRQ(ierr); 136 137 /* Compute function */ 138 ff[0] = -2.0 + 2.0*xx[0] + 400.0*xx[0]*xx[0]*xx[0] - 400.0*xx[0]*xx[1]; 139 ff[1] = -200.0*xx[0]*xx[0] + 200.0*xx[1]; 140 141 /* Restore vectors */ 142 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 143 ierr = VecRestoreArray(f,&ff);CHKERRQ(ierr); 144 return 0; 145 } 146 /* ------------------------------------------------------------------- */ 147 /* 148 FormJacobian1 - Evaluates Jacobian matrix. 149 150 Input Parameters: 151 . snes - the SNES context 152 . x - input vector 153 . dummy - optional user-defined context (not used here) 154 155 Output Parameters: 156 . jac - Jacobian matrix 157 . B - optionally different preconditioning matrix 158 . flag - flag indicating matrix structure 159 */ 160 PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat jac,Mat B,void *dummy) 161 { 162 const PetscScalar *xx; 163 PetscScalar A[4]; 164 PetscErrorCode ierr; 165 PetscInt idx[2] = {0,1}; 166 167 /* 168 Get pointer to vector data 169 */ 170 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 171 172 /* 173 Compute Jacobian entries and insert into matrix. 174 - Since this is such a small problem, we set all entries for 175 the matrix at once. 176 */ 177 A[0] = 2.0 + 1200.0*xx[0]*xx[0] - 400.0*xx[1]; 178 A[1] = -400.0*xx[0]; 179 A[2] = -400.0*xx[0]; 180 A[3] = 200; 181 ierr = MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES);CHKERRQ(ierr); 182 183 /* 184 Restore vector 185 */ 186 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 187 188 /* 189 Assemble matrix 190 */ 191 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 192 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 193 if (jac != B) { 194 ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 195 ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 196 } 197 return 0; 198 } 199 200 /*TEST 201 202 test: 203 args: -snes_monitor_short -snes_max_it 1000 204 requires: !single 205 206 TEST*/ 207