!
!  This example shows how to avoid Fortran line lengths larger than 132 characters.
!  It avoids used of certain macros such as PetscCallA() and PetscCheckA() that
!  generate very long lines
!
!  We recommend starting from src/snes/tutorials/ex5f90.F90 instead of this example
!  because that does not have the restricted formatting that makes this version
!  more difficult to read
!
!  Description: This example solves a nonlinear system in parallel with SNES.
!  We solve the  Bratu (SFI - solid fuel ignition) problem in a 2D rectangular
!  domain, using distributed arrays (DMDAs) to partition the parallel grid.
!  The command line options include:
!    -par <param>, where <param> indicates the nonlinearity of the problem
!       problem SFI:  <parameter> = Bratu parameter (0 <= par <= 6.81)
!
!  --------------------------------------------------------------------------
!
!  Solid Fuel Ignition (SFI) problem.  This problem is modeled by
!  the partial differential equation
!
!          -Laplacian u - lambda*exp(u) = 0,  0 < x,y < 1,
!
!  with boundary conditions
!
!           u = 0  for  x = 0, x = 1, y = 0, y = 1.
!
!  A finite difference approximation with the usual 5-point stencil
!  is used to discretize the boundary value problem to obtain a nonlinear
!  system of equations.
!
!  --------------------------------------------------------------------------
#include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h>
module ex5f_mod
  use petscsnes
  use petscdmda
  implicit none
  PetscInt xs, xe, xm, gxs, gxe, gxm
  PetscInt ys, ye, ym, gys, gye, gym
  PetscInt mx, my
  PetscMPIInt rank, size
  PetscReal lambda
contains
! ---------------------------------------------------------------------
!
!  FormInitialGuess - Forms initial approximation.
!
!  Input Parameters:
!  X - vector
!
!  Output Parameter:
!  X - vector
!
!  Notes:
!  This routine serves as a wrapper for the lower-level routine
!  "ApplicationInitialGuess", where the actual computations are
!  done using the standard Fortran style of treating the local
!  vector data as a multidimensional array over the local mesh.
!  This routine merely handles ghost point scatters and accesses
!  the local vector data via VecGetArray() and VecRestoreArray().
!
  subroutine FormInitialGuess(X, ierr)

!  Input/output variables:
    Vec X
    PetscErrorCode ierr
!  Declarations for use with local arrays:
    PetscScalar, pointer :: lx_v(:)

    ierr = 0

!  Get a pointer 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.
!    - Note that the Fortran interface to VecGetArray() differs from the
!      C version.  See the users manual for details.

    call VecGetArray(X, lx_v, ierr)
    CHKERRQ(ierr)

!  Compute initial guess over the locally owned part of the grid

    call InitialGuessLocal(lx_v, ierr)
    CHKERRQ(ierr)

!  Restore vector

    call VecRestoreArray(X, lx_v, ierr)
    CHKERRQ(ierr)

  end

! ---------------------------------------------------------------------
!
!  InitialGuessLocal - Computes initial approximation, called by
!  the higher level routine FormInitialGuess().
!
!  Input Parameter:
!  x - local vector data
!
!  Output Parameters:
!  x - local vector data
!  ierr - error code
!
!  Notes:
!  This routine uses standard Fortran-style computations over a 2-dim array.
!
  subroutine InitialGuessLocal(x, ierr)

!  Input/output variables:
    PetscScalar x(xs:xe, ys:ye)
    PetscErrorCode ierr

!  Local variables:
    PetscInt i, j
    PetscReal temp1, temp, one, hx, hy

!  Set parameters

    ierr = 0
    one = 1.0
    hx = one/((real(mx) - 1))
    hy = one/((real(my) - 1))
    temp1 = lambda/(lambda + one)

    do j = ys, ye
      temp = (real(min(j - 1, my - j)))*hy
      do i = xs, xe
        if (i == 1 .or. j == 1 .or. i == mx .or. j == my) then
          x(i, j) = 0.0
        else
          x(i, j) = temp1*sqrt(min(real(min(i - 1, mx - i))*hx, (temp)))
        end if
      end do
    end do

  end

! ---------------------------------------------------------------------
!
!  FormFunctionLocal - Computes nonlinear function, called by
!  the higher level routine FormFunction().
!
!  Input Parameter:
!  x - local vector data
!
!  Output Parameters:
!  f - local vector data, f(x)
!  ierr - error code
!
!  Notes:
!  This routine uses standard Fortran-style computations over a 2-dim array.
!
!
  subroutine FormFunctionLocal(info, x, f, da, ierr)

    DM da

!  Input/output variables:
    DMDALocalInfo info
    PetscScalar x(gxs:gxe, gys:gye)
    PetscScalar f(xs:xe, ys:ye)
    PetscErrorCode ierr

!  Local variables:
    PetscScalar two, one, hx, hy
    PetscScalar hxdhy, hydhx, sc
    PetscScalar u, uxx, uyy
    PetscInt i, j

    xs = info%XS + 1
    xe = xs + info%XM - 1
    ys = info%YS + 1
    ye = ys + info%YM - 1
    mx = info%MX
    my = info%MY

    one = 1.0
    two = 2.0
    hx = one/(real(mx) - 1)
    hy = one/(real(my) - 1)
    sc = hx*hy*lambda
    hxdhy = hx/hy
    hydhx = hy/hx

!  Compute function over the locally owned part of the grid

    do j = ys, ye
      do i = xs, xe
        if (i == 1 .or. j == 1 .or. i == mx .or. j == my) then
          f(i, j) = x(i, j)
        else
          u = x(i, j)
          uxx = hydhx*(two*u - x(i - 1, j) - x(i + 1, j))
          uyy = hxdhy*(two*u - x(i, j - 1) - x(i, j + 1))
          f(i, j) = uxx + uyy - sc*exp(u)
        end if
      end do
    end do

    call PetscLogFlops(11.0d0*ym*xm, ierr)
    CHKERRQ(ierr)

  end

! ---------------------------------------------------------------------
!
!  FormJacobianLocal - Computes Jacobian matrix, called by
!  the higher level routine FormJacobian().
!
!  Input Parameters:
!  x        - local vector data
!
!  Output Parameters:
!  jac      - Jacobian matrix
!  jac_prec - optionally different matrix used to construct the preconditioner (not used here)
!  ierr     - error code
!
!  Notes:
!  This routine uses standard Fortran-style computations over a 2-dim array.
!
!  Notes:
!  Due to grid point reordering with DMDAs, we must always work
!  with the local grid points, and then transform them to the new
!  global numbering with the "ltog" mapping
!  We cannot work directly with the global numbers for the original
!  uniprocessor grid!
!
!  Two methods are available for imposing this transformation
!  when setting matrix entries:
!    (A) MatSetValuesLocal(), using the local ordering (including
!        ghost points!)
!          by calling MatSetValuesLocal()
!    (B) MatSetValues(), using the global ordering
!        - Use DMDAGetGlobalIndices() to extract the local-to-global map
!        - Then apply this map explicitly yourself
!        - Set matrix entries using the global ordering by calling
!          MatSetValues()
!  Option (A) seems cleaner/easier in many cases, and is the procedure
!  used in this example.
!
  subroutine FormJacobianLocal(info, x, A, jac, da, ierr)

    DM da

!  Input/output variables:
    PetscScalar x(gxs:gxe, gys:gye)
    Mat A, jac
    PetscErrorCode ierr
    DMDALocalInfo info

!  Local variables:
    PetscInt row, col(5), i, j, i1, i5
    PetscScalar two, one, hx, hy, v(5)
    PetscScalar hxdhy, hydhx, sc

!  Set parameters

    i1 = 1
    i5 = 5
    one = 1.0
    two = 2.0
    hx = one/(real(mx) - 1)
    hy = one/(real(my) - 1)
    sc = hx*hy
    hxdhy = hx/hy
    hydhx = hy/hx
! -Wmaybe-uninitialized
    v = 0.0
    col = 0

!  Compute entries for the locally owned part of the Jacobian.
!   - Currently, all PETSc parallel matrix formats are partitioned by
!     contiguous chunks of rows across the processors.
!   - Each processor needs to insert only elements that it owns
!     locally (but any non-local elements will be sent to the
!     appropriate processor during matrix assembly).
!   - Here, we set all entries for a particular row at once.
!   - We can set matrix entries either using either
!     MatSetValuesLocal() or MatSetValues(), as discussed above.
!   - Note that MatSetValues() uses 0-based row and column numbers
!     in Fortran as well as in C.

    do j = ys, ye
      row = (j - gys)*gxm + xs - gxs - 1
      do i = xs, xe
        row = row + 1
!           boundary points
        if (i == 1 .or. j == 1 .or. i == mx .or. j == my) then
!       Some f90 compilers need 4th arg to be of same type in both calls
          col(1) = row
          v(1) = one
          call MatSetValuesLocal(jac, i1, [row], i1, [col], [v], INSERT_VALUES, ierr)
          CHKERRQ(ierr)
!           interior grid points
        else
          v(1) = -hxdhy
          v(2) = -hydhx
          v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(x(i, j))
          v(4) = -hydhx
          v(5) = -hxdhy
          col(1) = row - gxm
          col(2) = row - 1
          col(3) = row
          col(4) = row + 1
          col(5) = row + gxm
          call MatSetValuesLocal(jac, i1, [row], i5, [col], [v], INSERT_VALUES, ierr)
          CHKERRQ(ierr)
        end if
      end do
    end do
    call MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr)
    CHKERRQ(ierr)
    call MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr)
    CHKERRQ(ierr)
    if (A /= jac) then
      call MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)
      CHKERRQ(ierr)
      call MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)
      CHKERRQ(ierr)
    end if
  end

!
!     Simple convergence test based on the infinity norm of the residual being small
!
  subroutine MySNESConverged(snes, it, xnorm, snorm, fnorm, reason, dummy, ierr)

    SNES snes
    PetscInt it, dummy
    PetscReal xnorm, snorm, fnorm, nrm
    SNESConvergedReason reason
    Vec f
    PetscErrorCode ierr

    call SNESGetFunction(snes, f, PETSC_NULL_FUNCTION, dummy, ierr)
    CHKERRQ(ierr)
    call VecNorm(f, NORM_INFINITY, nrm, ierr)
    CHKERRQ(ierr)
    if (nrm <= 1.e-5) reason = SNES_CONVERGED_FNORM_ABS

  end

end module ex5f_mod

program main
  use ex5f_mod
  implicit none

! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!                   Variable declarations
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
!  Variables:
!     snes        - nonlinear solver
!     x, r        - solution, residual vectors
!     its         - iterations for convergence
!
!  See additional variable declarations in the file ex5f.h
!
  SNES snes
  Vec x, r
  PetscInt its, i1, i4
  PetscErrorCode ierr
  PetscReal lambda_max, lambda_min
  PetscBool flg
  DM da

! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  Initialize program
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  call PetscInitialize(ierr)
  CHKERRA(ierr)
  call MPI_Comm_size(PETSC_COMM_WORLD, size, ierr)
  CHKERRMPIA(ierr)
  call MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr)
  CHKERRMPIA(ierr)
!  Initialize problem parameters

  i1 = 1
  i4 = 4
  lambda_max = 6.81
  lambda_min = 0.0
  lambda = 6.0
  call PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-par', lambda, PETSC_NULL_BOOL, ierr)
  CHKERRA(ierr)

! this statement is split into multiple-lines to keep lines under 132 char limit - required by 'make check'
  if (lambda >= lambda_max .or. lambda <= lambda_min) then
    ierr = PETSC_ERR_ARG_OUTOFRANGE
    SETERRA(PETSC_COMM_WORLD, ierr, 'Lambda')
  end if

! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  Create nonlinear solver context
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  call SNESCreate(PETSC_COMM_WORLD, snes, ierr)
  CHKERRA(ierr)

!  Set convergence test routine if desired

  call PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-my_snes_convergence', flg, ierr)
  CHKERRA(ierr)
  if (flg) then
    call SNESSetConvergenceTest(snes, MySNESConverged, 0, PETSC_NULL_FUNCTION, ierr)
    CHKERRA(ierr)
  end if

! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  Create vector data structures; set function evaluation routine
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

!  Create distributed array (DMDA) to manage parallel grid and vectors

!     This really needs only the star-type stencil, but we use the box stencil

  call DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, i4, i4, PETSC_DECIDE, PETSC_DECIDE, &
                    i1, i1, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, da, ierr)
  CHKERRA(ierr)
  call DMSetFromOptions(da, ierr)
  CHKERRA(ierr)
  call DMSetUp(da, ierr)
  CHKERRA(ierr)

!  Extract global and local vectors from DMDA; then duplicate for remaining
!  vectors that are the same types

  call DMCreateGlobalVector(da, x, ierr)
  CHKERRA(ierr)
  call VecDuplicate(x, r, ierr)
  CHKERRA(ierr)

!  Get local grid boundaries (for 2-dimensional DMDA)

  call DMDAGetInfo(da, PETSC_NULL_INTEGER, mx, my, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, &
                   PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, &
                   PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMDASTENCILTYPE, ierr)
  CHKERRA(ierr)
  call DMDAGetCorners(da, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr)
  CHKERRA(ierr)
  call DMDAGetGhostCorners(da, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr)
  CHKERRA(ierr)

!  Here we shift the starting indices up by one so that we can easily
!  use the Fortran convention of 1-based indices (rather 0-based indices).

  xs = xs + 1
  ys = ys + 1
  gxs = gxs + 1
  gys = gys + 1

  ye = ys + ym - 1
  xe = xs + xm - 1
  gye = gys + gym - 1
  gxe = gxs + gxm - 1

!  Set function evaluation routine and vector

  call DMDASNESSetFunctionLocal(da, INSERT_VALUES, FormFunctionLocal, da, ierr)
  CHKERRA(ierr)
  call DMDASNESSetJacobianLocal(da, FormJacobianLocal, da, ierr)
  CHKERRA(ierr)
  call SNESSetDM(snes, da, ierr)
  CHKERRA(ierr)

! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  Customize nonlinear solver; set runtime options
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

!  Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)

  call SNESSetFromOptions(snes, ierr)
  CHKERRA(ierr)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  Evaluate initial guess; then solve nonlinear system.
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

!  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().

  call FormInitialGuess(x, ierr)
  CHKERRA(ierr)
  call SNESSolve(snes, PETSC_NULL_VEC, x, ierr)
  CHKERRA(ierr)
  call SNESGetIterationNumber(snes, its, ierr)
  CHKERRA(ierr)
  if (rank == 0) then
    write (6, 100) its
  end if
100 format('Number of SNES iterations = ', i5)

! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  Free work space.  All PETSc objects should be destroyed when they
!  are no longer needed.
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  call VecDestroy(x, ierr)
  CHKERRA(ierr)
  call VecDestroy(r, ierr)
  CHKERRA(ierr)
  call SNESDestroy(snes, ierr)
  CHKERRA(ierr)
  call DMDestroy(da, ierr)
  CHKERRA(ierr)
  call PetscFinalize(ierr)
  CHKERRA(ierr)
end
!/*TEST
!
!   build:
!      requires: !complex !single
!
!   test:
!      nsize: 4
!      args: -snes_mf -pc_type none -da_processors_x 4 -da_processors_y 1 -snes_monitor_short \
!            -ksp_gmres_cgs_refinement_type refine_always
!
!   test:
!      suffix: 2
!      nsize: 4
!      args: -da_processors_x 2 -da_processors_y 2 -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
!
!   test:
!      suffix: 3
!      nsize: 3
!      args: -snes_fd -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
!
!   test:
!      suffix: 6
!      nsize: 1
!      args: -snes_monitor_short -my_snes_convergence
!
!TEST*/