1 static char help[] = "Newton's method to solve a two-variable system that comes from the Rosenbrock function.\n\n"; 2 3 /* 4 Include "petscsnes.h" so that we can use SNES solvers. Note that this 5 file automatically includes: 6 petscsys.h - base PETSc routines petscvec.h - vectors 7 petscmat.h - matrices 8 petscis.h - index sets petscksp.h - Krylov subspace methods 9 petscviewer.h - viewers petscpc.h - preconditioners 10 petscksp.h - linear solvers 11 */ 12 #include <petscsnes.h> 13 14 extern PetscErrorCode FormJacobian1(SNES, Vec, Mat, Mat, void *); 15 extern PetscErrorCode FormFunction1(SNES, Vec, Vec, void *); 16 17 int main(int argc, char **argv) 18 { 19 SNES snes; /* nonlinear solver context */ 20 Vec x, r; /* solution, residual vectors */ 21 Mat J; /* Jacobian matrix */ 22 PetscInt its; 23 PetscScalar *xx; 24 SNESConvergedReason reason; 25 26 PetscFunctionBeginUser; 27 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 28 29 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 30 Create nonlinear solver context 31 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 32 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 33 34 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 35 Create matrix and vector data structures; set corresponding routines 36 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 37 /* 38 Create vectors for solution and nonlinear function 39 */ 40 PetscCall(VecCreate(PETSC_COMM_WORLD, &x)); 41 PetscCall(VecSetSizes(x, PETSC_DECIDE, 2)); 42 PetscCall(VecSetFromOptions(x)); 43 PetscCall(VecDuplicate(x, &r)); 44 45 /* 46 Create Jacobian matrix data structure 47 */ 48 PetscCall(MatCreate(PETSC_COMM_WORLD, &J)); 49 PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, 2, 2)); 50 PetscCall(MatSetFromOptions(J)); 51 PetscCall(MatSetUp(J)); 52 53 /* 54 Set function evaluation routine and vector. 55 */ 56 PetscCall(SNESSetFunction(snes, r, FormFunction1, NULL)); 57 58 /* 59 Set Jacobian matrix data structure and Jacobian evaluation routine 60 */ 61 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian1, NULL)); 62 63 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 64 Customize nonlinear solver; set runtime options 65 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 66 PetscCall(SNESSetFromOptions(snes)); 67 68 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 69 Evaluate initial guess; then solve nonlinear system 70 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 71 PetscCall(VecGetArray(x, &xx)); 72 xx[0] = -1.2; 73 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)); 101 PetscCall(VecDestroy(&r)); 102 PetscCall(MatDestroy(&J)); 103 PetscCall(SNESDestroy(&snes)); 104 PetscCall(PetscFinalize()); 105 return 0; 106 } 107 /* ------------------------------------------------------------------- */ 108 /* 109 FormFunction1 - Evaluates nonlinear function, F(x). 110 111 Input Parameters: 112 . snes - the SNES context 113 . x - input vector 114 . ctx - optional user-defined context 115 116 Output Parameter: 117 . f - function vector 118 */ 119 PetscErrorCode FormFunction1(SNES snes, Vec x, Vec f, void *ctx) 120 { 121 PetscScalar *ff; 122 const PetscScalar *xx; 123 124 PetscFunctionBeginUser; 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 PetscCall(VecGetArrayRead(x, &xx)); 133 PetscCall(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 PetscCall(VecRestoreArrayRead(x, &xx)); 141 PetscCall(VecRestoreArray(f, &ff)); 142 PetscFunctionReturn(PETSC_SUCCESS); 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 PetscFunctionBeginUser; 165 /* 166 Get pointer to vector data 167 */ 168 PetscCall(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 PetscCall(MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES)); 180 181 /* 182 Restore vector 183 */ 184 PetscCall(VecRestoreArrayRead(x, &xx)); 185 186 /* 187 Assemble matrix 188 */ 189 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 190 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 191 if (jac != B) { 192 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 193 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 194 } 195 PetscFunctionReturn(PETSC_SUCCESS); 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