How to split a single Fortran file into separate subroutine and function files

Question

Another week, another silly question from someone trying not to smash their head against their desk. I have a Fortran90 file that does what I want it to do, no errors, and correct output. Now I need to separate this thing that works into its constituents, i.e. two subroutine files, two function files and a driver program. How do I do this without breaking it, because it is broken...The main issue is passing arrays...I think. Working single file code:

      program testnew
              implicit none
              integer,parameter :: p14r300 = SELECTED_REAL_KIND(14,300)
              integer,parameter :: k7 = SELECTED_INT_KIND(7)
              integer(kind=k7) :: n, ng
              real(kind=p14r300), dimension(:), allocatable :: xarr
              real(kind=p14r300), dimension(:), allocatable :: xabsc
              real(kind=p14r300), dimension(:), allocatable :: weight
              real(kind=p14r300) :: tol  

              do n=2,4
              allocate (xarr(n))
              allocate (xabsc(n))
              allocate (weight(n))
              call gauss_leg_int(ng, xabsc, weight)
              print *, ng, xabsc, weight
              deallocate (xarr)
              deallocate (xabsc)
              deallocate (weight)
              enddo

      return
      contains
             subroutine gauss_leg_int(ng, xabsc, weight)
      !==================================================================
      ! Subroutine that organizes the computations to find the abscissas
      ! and weights for Gauss-Legendre integration, where ng is the
      ! number of integration points(integer, input), and xabsc and
      ! weight are real arrays of length ng (output) that hold the
      ! abscissas and weights, respectively.
      !==================================================================
             integer(kind=k7) :: ng, i, iter
             real(kind=p14r300) :: x, w
       real(kind=p14r300), dimension(:), allocatable :: weight, xabsc

             do i=1,n
             call leg_root(n, tol, xarr)
             xabsc=xarr
             ng=n
             !do iter=1,n
             x=xabsc(i)
             print *,x
             w=2/((1-x**2)*leg_deriv(n, x)**2)
             !enddo
             weight(i)=w
             enddo
             end subroutine gauss_leg_int

             subroutine leg_root(n, tol, xarr)
      !==================================================================
      ! Subroutine that finds the set of roots of a Legendre polynomial,
      ! where n is the degree of the polynomial (input,integer), and tol
      ! is an absolute tolerance(input,real) for stopping the iteration
      ! when abs(P_l(x_i))<=tol.
      !==================================================================
             real(kind=p14r300) :: a, pi, x, y, pl, tol ! Declare real variables
             real(kind=p14r300), dimension(:), allocatable :: xarr ! Array
             integer(kind=k7) :: i, n, iter ! Declare integer variables
             a=1.0   ! Value to use on the next line
             pi=4*atan(a) ! Calculate Pi
             tol=1.d-14
             do i=1,n
             x=-cos(pi*(i-0.25)/(n+0.5))       ! Initial x value
             do iter=1,20    ! Set maximum number of iterations
             y=x-leg_poly(n, x)/leg_deriv(n, x)
             pl=leg_poly(n, y)-leg_poly(n, x)
             x=y     ! Once value of y is correct, make x the same
             if (abs(pl)<=tol) exit ! Once tolerance is reached, exit
             enddo
             !write (*,*) x
             xarr(i)=x
             !print *,xarr
             enddo
             !xarr(1,i*4)=x
             end subroutine leg_root

      function leg_poly(n, x) result(pn)
      !==================================================================
      ! Function for evaluating a given Legendre polynomial using the
      ! recurrence relation, where n is the degree of the
      ! polynomial(input, integer), and x is the location(input, real)
      ! in the interval -1<=x<=1 in which to evaluate the polynomial.
      ! The function result is the real value of P_n(x).
      !==================================================================
              real(kind=p14r300) :: pn, x, pln(0:n)
              integer(kind=k7) :: l, n

              pln(0)=1.0        ! First Legendre polynomial
              pln(1)=x          ! Second Legendre polynomial

              if (n<=1) then    ! Set the first two polynomials
                      pn=pln(n)
              else              ! Starts the recurrence to generate
                 do l=1,n-1     ! higher degree polynomials
                 pln(l+1)=((2.0*l+1.0)*x*pln(l)-l*pln(l-1))/(l+1)
                 enddo
                 pn=pln(n)
              endif
         end function leg_poly

       function leg_deriv(n, x) result(pdn)
       !=================================================================
       ! Function for evaluating the derivatives of a given Legendre
       ! polynomial using the recurrence relation, where n is the degree
       ! of the polynomial(input, integer), and x is the
       ! location(input, real) in the interval -1<=x<=1 in which to
       ! evaluate the derivative. The function result is the real value
       ! of Pd_n(x).
       !=================================================================
            real(kind=p14r300) :: pdn, x, pdln(0:n)
            integer(kind=k7) :: l, n

            pdln(0)=0        ! Derivative of first Legendre polynomial
            pdln(1)=1.0      ! Derivative of second Legendre polynomial

              if (n<=1) then   ! Set the first two Legendre polynomial
               pdn=pdln(n)     ! derivatives
               else            ! Starts the recurrence to generate
                  do l=1,n-1   ! higher degree polynomial derivatives
                   pdln(l+1)=((2.0*l+1.0)*x*pdln(l)-(l+1)*pdln(l-1))/l
                  enddo
                  pdn=pdln(n)
               endif
               end function leg_deriv
               end program

A module is what you want - is that what you have tried? Also note you are most likely losing precision because the kind of your real variables and the kind of your real constants is not consistent — Ian Bush, Feb 27 '20 at 07:10
Also it would be A LOT more useful if you posted the code that is not working and describe carefully how it is not working (doesn't compile? always wrong answers? sometimes wrong answers? ..) — Ian Bush, Feb 27 '20 at 07:32
@Ian_Bush I managed to get it done, I used a module for precision parameters. I am interested to know what you mean by losing precision. I will post my version here in a few... — CrazyIvan1978, Feb 27 '20 at 11:50
@CrazyIvan1978 I do not know what exactly Ian had in mind but in the subroutines you use single precision values a lot. For example `pi=4*atan(a)` or `x=-cos(pi*(i-0.25)/(n+0.5))`. — PetrH, Feb 27 '20 at 12:00
See https://stackoverflow.com/questions/6146005/is-there-a-better-double-precision-assignment-in-fortran-90/42622204#42622204 However, @PeterH there is no single precision in `pi=4*atan(a)`. In the latter case there are single precision constants there, but they are exactly representable. — Vladimir F Героям слава, Feb 27 '20 at 12:02
Thank you for your comments. @PetrH, so having these 0.25 and 0.5 and having a=1.0 messes up precision? — CrazyIvan1978, Feb 27 '20 at 13:41
@VladimirF Thanks, in the first case I meant that `4` is coerced to a single-precision real, or isn't it? I think that it is. This may not make a difference with the value `4` but would generally. In the second case, you get values like `3.75`, is this exactly representable? — PetrH, Feb 27 '20 at 14:42
@CrazyIvan1978 I would think so, although I am not sure about `a=1.0`, this is most likely OK as you declared `a` to be double precision. Generally, you want to use `0.25_p14r300`. — PetrH, Feb 27 '20 at 14:45
@PetrH No, 4 is converted to the kind of the rest of the real expression. It does not matter anyway as small integers are exactly representable in all reasonable floating point and fixed point (both binary and decadic) numeric formats. The same holds for numbers like 1/2 1/8 and their combinations like 3.75. 3.75 is exactly representable in reasonable binary and decadic formats. It is still better to use double precision literals for these values. Or, as I do, just divide by an integer 1,2,4,8... — Vladimir F Героям слава, Feb 27 '20 at 14:57
@PetrH Numbers that are not exactly representably are the like 0.1, 0.2, 0.3 and similar. These would only be exactly representable in decimal formats, not in binary. In the previous comment I should have used decimal, not decadic. — Vladimir F Героям слава, Feb 27 '20 at 15:00

score 0 · Answer 1 · answered Feb 27 '20 at 13:50

Here is how solved my conundrum, first I divided it all up, created a module for precision parameters, and edited each file to deal with the fact that arrays had to be passed: Precision MODULE:

  MODULE Precision
        !===========================================================
        ! Module to be used to declare precision parameters for any
        ! program using it.
        !===========================================================        

                IMPLICIT NONE
                INTEGER, PARAMETER :: p14r300=SELECTED_REAL_KIND(14,300)
                INTEGER, PARAMETER :: k7=SELECTED_INT_KIND(7)
        END MODULE Precision

Driver program:

program testnew
              USE Precision
              implicit none
              integer(kind=k7) :: n, ng
              real(kind=p14r300), allocatable :: xabsc(:)
              real(kind=p14r300), allocatable :: weight(:)

              do n=2,4     ! Set range based on provided table

        ! Allocate memory to dynamic arrays
              allocate (xabsc(n))
              allocate (weight(n))

        ! Call subroutine to obtain values of interest
              call gauss_leg_int(n, xabsc, weight)

        ! Print values to stdout
              print *, n
              print *, xabsc
              print *, weight

        ! Deallocate memory
              deallocate (xabsc)
              deallocate (weight)
              enddo       

       end program testnew

Functions:


      function leg_poly(n, x) result(pn)
      !==================================================================
      ! Function for evaluating a given Legendre polynomial using the
      ! recurrence relation, where n is the degree of the
      ! polynomial(input, integer), and x is the location(input, real)
      ! in the interval -1<=x<=1 in which to evaluate the polynomial.
      ! The function result is the real value of P_n(x).
      !==================================================================
              USE Precision
              IMPLICIT NONE
              real(kind=p14r300) :: pn, pln(0:n)
              real(kind=p14r300), intent(in) :: x
              integer(kind=k7), intent(in) :: n
              integer(kind=k7) :: l
              pln(0)=1.0        ! First Legendre polynomial
              pln(1)=x          ! Second Legendre polynomial

              if (n<=1) then    ! Set the first two polynomials
                      pn=pln(n)
              else              ! Starts the recurrence to generate
                 do l=1,n-1     ! higher degree polynomials
                 pln(l+1)=((2.0*l+1.0)*x*pln(l)-l*pln(l-1))/(l+1)
                 enddo
                 pn=pln(n)
              endif
         end function leg_poly


       function leg_deriv(n, x) result(pdn)
       !=================================================================
       ! Function for evaluating the derivatives of a given Legendre
       ! polynomial using the recurrence relation, where n is the degree
       ! of the polynomial(input, integer), and x is the
       ! location(input, real) in the interval -1<=x<=1 in which to
       ! evaluate the derivative. The function result is the real value
       ! of Pd_n(x).
       !=================================================================
            USE Precision
            IMPLICIT NONE
            REAL(KIND=p14r300) :: pdn, x, pdln(0:n)
            INTEGER(KIND=k7), INTENT(IN) :: n
            INTEGER(KIND=k7) :: l
            pdln(0)=0        ! Derivative of first Legendre polynomial
            pdln(1)=1.0      ! Derivative of second Legendre polynomial

              if (n<=1) then   ! Set the first two Legendre polynomial
               pdn=pdln(n)     ! derivatives
               else            ! Starts the recurrence to generate
                  do l=1,n-1   ! higher degree polynomial derivatives
                   pdln(l+1)=((2.0*l+1.0)*x*pdln(l)-(l+1)*pdln(l-1))/l
                  enddo
                  pdn=pdln(n)
               endif
               end function leg_deriv

Subroutines:


            subroutine leg_root(n, xarr)
      !==================================================================
      ! Subroutine that finds the set of roots of a Legendre polynomial,
      ! where n is the degree of the polynomial (input,integer), and tol
      ! is an absolute tolerance(input,real) for stopping the iteration
      ! when abs(P_l(x_i))<=tol.
      !==================================================================
             USE Precision
             IMPLICIT NONE
             real(kind=p14r300) :: a, pi, x, y, pl ! Declare real variables
             real(kind=p14r300) :: tol
             real(kind=p14r300), intent(out) :: xarr(n)
             integer(kind=k7) :: i, iter ! Declare integer variables
             integer(kind=k7), intent(in) :: n
             real(kind=p14r300),EXTERNAL :: leg_poly, leg_deriv


             a=1.0   ! Value to use on the next line
             pi=4*atan(a) ! Calculate Pi
             tol=1.d-14
             do i=1,n
             x=-cos(pi*(i-0.25)/(n+0.5))       ! Initial x value
             do iter=1,20    ! Set maximum number of iterations
             y=x-leg_poly(n, x)/leg_deriv(n, x)
             pl=leg_poly(n, y)-leg_poly(n, x)
             x=y     ! Make x the same as the computed y to repeat
             ! calculation 
             if (abs(pl)<=tol) exit ! Once tolerance is reached, exit
             enddo
             xarr(i)=x   ! Place x values into array
             enddo
             end subroutine leg_root            



            subroutine gauss_leg_int(ng, xabsc, weight)
      !==================================================================
      ! Subroutine that organizes the computations to find the abscissas
      ! and weights for Gauss-Legendre integration, where ng is the
      ! number of integration points(integer, input), and xabsc and
      ! weight are real arrays of length ng (output) that hold the
      ! abscissas and weights, respectively.
      !==================================================================
             USE Precision
             IMPLICIT NONE
             integer(kind=k7) :: i 
             integer(kind=k7), intent(in) :: ng
             real(kind=p14r300) :: x, w
             real(kind=p14r300), EXTERNAL :: leg_deriv
             real(kind=p14r300) :: xarr(ng)
             real(kind=p14r300), intent(out) :: weight(ng)
             real(kind=p14r300), intent(out) :: xabsc(ng)

             do i=1,ng
             call leg_root(ng, xabsc) ! Call subroutine to use xarr
             x=xabsc(i)      ! Loop over each x value per ng
             w=2/((1-x**2)*leg_deriv(ng, x)**2)  ! calculate weight
             weight(i)=w       ! Place weight values into array
             enddo
             end subroutine gauss_leg_int

It's not entirely clear, did you put the routines into a module ? — High Performance Mark, Feb 27 '20 at 13:54
Just to clarify Mark's comments - it is much better if you put the functions and subroutines in a module. In modern code external subprograms should be avoided. — Ian Bush, Feb 27 '20 at 14:46
I am learning FORTRAN90, I had to have it all broken up as part of my assignment, i.e. the module, main program, each subroutine and function had to be their own files. I would have rather had it all in one file, but alas. — CrazyIvan1978, Feb 27 '20 at 17:33
*each subroutine and function had to be their own files* If that is a requirement imposed on you by your educational institution ask for your money back and go elsewhere. Putting each routine in its own file, unless each is also in its own module, is worse than your original implementation with all the routines in a `contains` section in the same source as the `program`. Or perhaps you have misinterpreted the instructions and should have put all the routines in the same module as the parameter declarations. — High Performance Mark, Feb 27 '20 at 17:48
Oh come on Mark, no need to be so serious, it was to simplify grading, I think — CrazyIvan1978, Feb 27 '20 at 19:54

How to split a single Fortran file into separate subroutine and function files

1 Answers1