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; 74 xx[1] = 1.0; 75 PetscCall(VecRestoreArray(x, &xx)); 76 77 /* 78 Note: The user should initialize the vector, x, with the initial guess 79 for the nonlinear solver prior to calling SNESSolve(). In particular, 80 to employ an initial guess of zero, the user should explicitly set 81 this vector to zero by calling VecSet(). 82 */ 83 84 PetscCall(SNESSolve(snes, NULL, x)); 85 PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD)); 86 PetscCall(SNESGetIterationNumber(snes, &its)); 87 PetscCall(SNESGetConvergedReason(snes, &reason)); 88 /* 89 Some systems computes a residual that is identically zero, thus converging 90 due to FNORM_ABS, others converge due to FNORM_RELATIVE. Here, we only 91 report the reason if the iteration did not converge so that the tests are 92 reproducible. 93 */ 94 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s number of SNES iterations = %" PetscInt_FMT "\n", reason > 0 ? "CONVERGED" : SNESConvergedReasons[reason], its)); 95 96 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 97 Free work space. All PETSc objects should be destroyed when they 98 are no longer needed. 99 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 100 101 PetscCall(VecDestroy(&x)); 102 PetscCall(VecDestroy(&r)); 103 PetscCall(MatDestroy(&J)); 104 PetscCall(SNESDestroy(&snes)); 105 PetscCall(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 PetscFunctionBeginUser; 126 /* 127 Get pointers to vector data. 128 - For default PETSc vectors, VecGetArray() returns a pointer to 129 the data array. Otherwise, the routine is implementation dependent. 130 - You MUST call VecRestoreArray() when you no longer need access to 131 the array. 132 */ 133 PetscCall(VecGetArrayRead(x, &xx)); 134 PetscCall(VecGetArray(f, &ff)); 135 136 /* Compute function */ 137 ff[0] = -2.0 + 2.0 * xx[0] + 400.0 * xx[0] * xx[0] * xx[0] - 400.0 * xx[0] * xx[1]; 138 ff[1] = -200.0 * xx[0] * xx[0] + 200.0 * xx[1]; 139 140 /* Restore vectors */ 141 PetscCall(VecRestoreArrayRead(x, &xx)); 142 PetscCall(VecRestoreArray(f, &ff)); 143 PetscFunctionReturn(PETSC_SUCCESS); 144 } 145 /* ------------------------------------------------------------------- */ 146 /* 147 FormJacobian1 - Evaluates Jacobian matrix. 148 149 Input Parameters: 150 . snes - the SNES context 151 . x - input vector 152 . dummy - optional user-defined context (not used here) 153 154 Output Parameters: 155 . jac - Jacobian matrix 156 . B - optionally different preconditioning matrix 157 . flag - flag indicating matrix structure 158 */ 159 PetscErrorCode FormJacobian1(SNES snes, Vec x, Mat jac, Mat B, void *dummy) 160 { 161 const PetscScalar *xx; 162 PetscScalar A[4]; 163 PetscInt idx[2] = {0, 1}; 164 165 PetscFunctionBeginUser; 166 /* 167 Get pointer to vector data 168 */ 169 PetscCall(VecGetArrayRead(x, &xx)); 170 171 /* 172 Compute Jacobian entries and insert into matrix. 173 - Since this is such a small problem, we set all entries for 174 the matrix at once. 175 */ 176 A[0] = 2.0 + 1200.0 * xx[0] * xx[0] - 400.0 * xx[1]; 177 A[1] = -400.0 * xx[0]; 178 A[2] = -400.0 * xx[0]; 179 A[3] = 200; 180 PetscCall(MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES)); 181 182 /* 183 Restore vector 184 */ 185 PetscCall(VecRestoreArrayRead(x, &xx)); 186 187 /* 188 Assemble matrix 189 */ 190 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 191 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 192 if (jac != B) { 193 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 194 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 195 } 196 PetscFunctionReturn(PETSC_SUCCESS); 197 } 198 199 /*TEST 200 201 test: 202 suffix: 1 203 args: -snes_monitor_short -snes_max_it 1000 204 requires: !single 205 206 test: 207 suffix: 2 208 args: -snes_monitor_short -snes_max_it 1000 -snes_type newtontrdc -snes_trdc_use_cauchy false 209 requires: !single 210 211 TEST*/ 212