1: ! Computes the integral of 2*x/(1+x^2) from x=0..1
  2: ! This is equal to the ln(2).

  4: ! Contributed by Mike McCourt <mccomic@iit.edu> and Nathan Johnston <johnnat@iit.edu>
  5: ! Fortran translation by Arko Bhattacharjee <a.bhattacharjee@mpie.de>
  6: #include <petsc/finclude/petscvec.h>
  7: program main
  8:   use petscvec
  9:   implicit none

 11:   PetscErrorCode :: ierr
 12:   PetscMPIInt :: rank, size
 13:   PetscInt   ::  rstart, rend, i, k, N
 14:   PetscInt, parameter   ::   numPoints = 1000000
 15:   PetscScalar  ::  dummy
 16:   PetscScalar, parameter  :: h = 1.0/numPoints
 17:   PetscScalar, pointer, dimension(:)  :: xarray
 18:   PetscScalar :: myResult = 0
 19:   Vec x, xend
 20:   character(len=PETSC_MAX_PATH_LEN) :: output
 21:   PetscInt, parameter :: one = 1

 23:   PetscCallA(PetscInitialize(ierr))

 25:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
 26:   PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))

 28:   ! Create a parallel vector.
 29:   ! Here we set up our x vector which will be given values below.
 30:   ! The xend vector is a dummy vector to find the value of the
 31:   ! elements at the endpoints for use in the trapezoid rule.

 33:   PetscCallA(VecCreate(PETSC_COMM_WORLD, x, ierr))
 34:   PetscCallA(VecSetSizes(x, PETSC_DECIDE, numPoints, ierr))
 35:   PetscCallA(VecSetFromOptions(x, ierr))
 36:   PetscCallA(VecGetSize(x, N, ierr))
 37:   PetscCallA(VecSet(x, myResult, ierr))
 38:   PetscCallA(VecDuplicate(x, xend, ierr))
 39:   myResult = 0.5
 40:   if (rank == 0) then
 41:     i = 0
 42:     PetscCallA(VecSetValues(xend, one, [i], [myResult], INSERT_VALUES, ierr))
 43:   end if

 45:   if (rank == size - 1) then
 46:     i = N - 1
 47:     PetscCallA(VecSetValues(xend, one, [i], [myResult], INSERT_VALUES, ierr))
 48:   end if

 50:   ! Assemble vector, using the 2-step process:
 51:   ! VecAssemblyBegin(), VecAssemblyEnd()
 52:   ! Computations can be done while messages are in transition
 53:   ! by placing code between these two statements.

 55:   PetscCallA(VecAssemblyBegin(xend, ierr))
 56:   PetscCallA(VecAssemblyEnd(xend, ierr))

 58:   ! Set the x vector elements.
 59:   ! i*h will return 0 for i=0 and 1 for i=N-1.
 60:   ! The function evaluated (2x/(1+x^2)) is defined above.
 61:   ! Each evaluation is put into the local array of the vector without message passing.

 63:   PetscCallA(VecGetOwnershipRange(x, rstart, rend, ierr))
 64:   PetscCallA(VecGetArray(x, xarray, ierr))
 65:   k = 1
 66:   do i = rstart, rend - 1
 67:     xarray(k) = real(i)*h
 68:     xarray(k) = func(xarray(k))
 69:     k = k + 1
 70:   end do
 71:   PetscCallA(VecRestoreArray(x, xarray, ierr))

 73:   ! Evaluates the integral.  First the sum of all the points is taken.
 74:   ! That result is multiplied by the step size for the trapezoid rule.
 75:   ! Then half the value at each endpoint is subtracted,
 76:   ! this is part of the composite trapezoid rule.

 78:   PetscCallA(VecSum(x, myResult, ierr))
 79:   myResult = myResult*h
 80:   PetscCallA(VecDot(x, xend, dummy, ierr))
 81:   myResult = myResult - h*dummy

 83:   !Return the value of the integral.

 85:   write (output, '(a,f9.6,a)') 'ln(2) is', real(myResult), '\n'           ! PetscScalar might be complex
 86:   PetscCallA(PetscPrintf(PETSC_COMM_WORLD, trim(output), ierr))
 87:   PetscCallA(VecDestroy(x, ierr))
 88:   PetscCallA(VecDestroy(xend, ierr))

 90:   PetscCallA(PetscFinalize(ierr))

 92: contains

 94:   function func(a)
 95:     use petscvec

 97:     implicit none
 98:     PetscScalar :: func
 99:     PetscScalar, intent(IN) :: a

101:     func = 2.0*a/(1.0 + a*a)

103:   end function func

105: end program

107: !/*TEST
108: !
109: !     test:
110: !       nsize: 1
111: !       output_file: output/ex18_1.out
112: !
113: !     test:
114: !       nsize: 2
115: !       suffix: 2
116: !       output_file: output/ex18_1.out
117: !
118: !TEST*/