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