static char help[] = "Newton's method to solve a two-variable system that comes from the Rosenbrock function.\n\n"; /*T Concepts: SNES^basic example T*/ /* Include "petscsnes.h" so that we can use SNES solvers. Note that this file automatically includes: petscsys.h - base PETSc routines petscvec.h - vectors petscmat.h - matrices petscis.h - index sets petscksp.h - Krylov subspace methods petscviewer.h - viewers petscpc.h - preconditioners petscksp.h - linear solvers */ #include extern PetscErrorCode FormJacobian1(SNES,Vec,Mat,Mat,void*); extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*); int main(int argc,char **argv) { SNES snes; /* nonlinear solver context */ Vec x,r; /* solution, residual vectors */ Mat J; /* Jacobian matrix */ PetscErrorCode ierr; PetscInt its; PetscScalar *xx; SNESConvergedReason reason; ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Create nonlinear solver context - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Create matrix and vector data structures; set corresponding routines - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* Create vectors for solution and nonlinear function */ ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr); ierr = VecSetSizes(x,PETSC_DECIDE,2);CHKERRQ(ierr); ierr = VecSetFromOptions(x);CHKERRQ(ierr); ierr = VecDuplicate(x,&r);CHKERRQ(ierr); /* Create Jacobian matrix data structure */ ierr = MatCreate(PETSC_COMM_WORLD,&J);CHKERRQ(ierr); ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);CHKERRQ(ierr); ierr = MatSetFromOptions(J);CHKERRQ(ierr); ierr = MatSetUp(J);CHKERRQ(ierr); /* Set function evaluation routine and vector. */ ierr = SNESSetFunction(snes,r,FormFunction1,NULL);CHKERRQ(ierr); /* Set Jacobian matrix data structure and Jacobian evaluation routine */ ierr = SNESSetJacobian(snes,J,J,FormJacobian1,NULL);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Customize nonlinear solver; set runtime options - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Evaluate initial guess; then solve nonlinear system - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = VecGetArray(x,&xx);CHKERRQ(ierr); xx[0] = -1.2; xx[1] = 1.0; ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); /* Note: The user should initialize the vector, x, with the initial guess for the nonlinear solver prior to calling SNESSolve(). In particular, to employ an initial guess of zero, the user should explicitly set this vector to zero by calling VecSet(). */ ierr = SNESSolve(snes,NULL,x);CHKERRQ(ierr); ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); ierr = SNESGetConvergedReason(snes,&reason);CHKERRQ(ierr); /* Some systems computes a residual that is identically zero, thus converging due to FNORM_ABS, others converge due to FNORM_RELATIVE. Here, we only report the reason if the iteration did not converge so that the tests are reproducible. */ ierr = PetscPrintf(PETSC_COMM_WORLD,"%s number of SNES iterations = %D\n",reason>0 ? "CONVERGED" : SNESConvergedReasons[reason],its);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Free work space. All PETSc objects should be destroyed when they are no longer needed. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&r);CHKERRQ(ierr); ierr = MatDestroy(&J);CHKERRQ(ierr); ierr = SNESDestroy(&snes);CHKERRQ(ierr); ierr = PetscFinalize(); return ierr; } /* ------------------------------------------------------------------- */ /* FormFunction1 - Evaluates nonlinear function, F(x). Input Parameters: . snes - the SNES context . x - input vector . ctx - optional user-defined context Output Parameter: . f - function vector */ PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx) { PetscErrorCode ierr; PetscScalar *ff; const PetscScalar *xx; /* Get pointers to vector data. - For default PETSc vectors, VecGetArray() returns a pointer to the data array. Otherwise, the routine is implementation dependent. - You MUST call VecRestoreArray() when you no longer need access to the array. */ ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); ierr = VecGetArray(f,&ff);CHKERRQ(ierr); /* Compute function */ ff[0] = -2.0 + 2.0*xx[0] + 400.0*xx[0]*xx[0]*xx[0] - 400.0*xx[0]*xx[1]; ff[1] = -200.0*xx[0]*xx[0] + 200.0*xx[1]; /* Restore vectors */ ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); ierr = VecRestoreArray(f,&ff);CHKERRQ(ierr); return 0; } /* ------------------------------------------------------------------- */ /* FormJacobian1 - Evaluates Jacobian matrix. Input Parameters: . snes - the SNES context . x - input vector . dummy - optional user-defined context (not used here) Output Parameters: . jac - Jacobian matrix . B - optionally different preconditioning matrix . flag - flag indicating matrix structure */ PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat jac,Mat B,void *dummy) { const PetscScalar *xx; PetscScalar A[4]; PetscErrorCode ierr; PetscInt idx[2] = {0,1}; /* Get pointer to vector data */ ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); /* Compute Jacobian entries and insert into matrix. - Since this is such a small problem, we set all entries for the matrix at once. */ A[0] = 2.0 + 1200.0*xx[0]*xx[0] - 400.0*xx[1]; A[1] = -400.0*xx[0]; A[2] = -400.0*xx[0]; A[3] = 200; ierr = MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES);CHKERRQ(ierr); /* Restore vector */ ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); /* Assemble matrix */ ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); if (jac != B) { ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); } return 0; } /*TEST test: args: -snes_monitor_short -snes_max_it 1000 requires: !single TEST*/