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>

      7: program main
      8: #include <petsc/finclude/petscvec.h>
      9:   use petscvec
     10:   implicit none

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

     24:   PetscCallA(PetscInitialize(ierr))

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

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

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

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

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

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

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

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

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

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

     84:   !Return the value of the integral.

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

     91:   PetscCallA(PetscFinalize(ierr))

     93:   contains

     95:     function func(a)
     96: #include <petsc/finclude/petscvec.h>
     97:       use petscvec

     99:       implicit none
    100:       PetscScalar :: func
    101:       PetscScalar,INTENT(IN) :: a

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

    105:     end function func

    107: end program

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