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