Actual source code: ex18f.F90

      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: module ex18fmodule
      8:   use petscvec
      9:   implicit none

     11: contains

     13:   pure function func(a)
     14:     PetscScalar :: func
     15:     PetscScalar, intent(in) :: a

     17:     func = 2.0*a/(1.0 + a**2)

     19:   end function func

     21: end module ex18fmodule

     23: program main
     24:   use petscvec
     25:   use ex18fmodule

     27:   implicit none
     28:   PetscErrorCode :: ierr
     29:   PetscMPIInt :: rank, size
     30:   PetscInt   ::  rstart, rend, i, k, N
     31:   PetscInt, parameter ::numPoints = 1000000
     32:   PetscScalar  ::  dummy
     33:   PetscScalar, parameter  :: h = 1.0/numPoints
     34:   PetscScalar, pointer, dimension(:)  :: xarray
     35:   PetscScalar :: myResult = 0
     36:   Vec x, xend
     37:   character(len=PETSC_MAX_PATH_LEN) :: output

     39:   PetscCallA(PetscInitialize(ierr))

     41:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
     42:   PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))

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

     49:   PetscCallA(VecCreate(PETSC_COMM_WORLD, x, ierr))
     50:   PetscCallA(VecSetSizes(x, PETSC_DECIDE, numPoints, ierr))
     51:   PetscCallA(VecSetFromOptions(x, ierr))
     52:   PetscCallA(VecGetSize(x, N, ierr))
     53:   PetscCallA(VecSet(x, myResult, ierr))
     54:   PetscCallA(VecDuplicate(x, xend, ierr))
     55:   myResult = 0.5
     56:   if (rank == 0) then
     57:     i = 0
     58:     PetscCallA(VecSetValues(xend, 1_PETSC_INT_KIND, [i], [myResult], INSERT_VALUES, ierr))
     59:   end if

     61:   if (rank == size - 1) then
     62:     i = N - 1
     63:     PetscCallA(VecSetValues(xend, 1_PETSC_INT_KIND, [i], [myResult], INSERT_VALUES, ierr))
     64:   end if

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

     71:   PetscCallA(VecAssemblyBegin(xend, ierr))
     72:   PetscCallA(VecAssemblyEnd(xend, ierr))

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

     79:   PetscCallA(VecGetOwnershipRange(x, rstart, rend, ierr))
     80:   PetscCallA(VecGetArray(x, xarray, ierr))
     81:   k = 1
     82:   do i = rstart, rend - 1
     83:     xarray(k) = func(real(i)*h)
     84:     k = k + 1
     85:   end do
     86:   PetscCallA(VecRestoreArray(x, xarray, ierr))

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

     93:   PetscCallA(VecSum(x, myResult, ierr))
     94:   myResult = myResult*h
     95:   PetscCallA(VecDot(x, xend, dummy, ierr))
     96:   myResult = myResult - h*dummy

     98:   !Return the value of the integral.

    100:   write (output, '(a,f9.6,a)') 'ln(2) is', real(myResult), '\n'           ! PetscScalar might be complex
    101:   PetscCallA(PetscPrintf(PETSC_COMM_WORLD, trim(output), ierr))
    102:   PetscCallA(VecDestroy(x, ierr))
    103:   PetscCallA(VecDestroy(xend, ierr))

    105:   PetscCallA(PetscFinalize(ierr))
    106: end program

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