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: end if
46: if (rank == size - 1) then
47: i = N - 1
48: PetscCallA(VecSetValues(xend, one, [i], [myResult], INSERT_VALUES, ierr))
49: end if
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(VecGetArray(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(VecRestoreArray(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*/