2

I am writing a n-dimensional numerical solver in Fortran. I created a module for the same which is called from the main program. I used external for calling the unknown function when I wrote for first order ode. However on replicating the result for more than one dimension using external I get the following error

Error: EXTERNAL attribute conflicts with DIMENSION attribute in 'f_x' 

This means my compiler cannot process an external function with array output.

I tried using an interface block to predefine the function parameters with fixed dimensions and I end up with this error

Error: PROCEDURE attribute conflicts with INTENT attribute in 'f_x' 

Is there any way in Fortran to have an external function returning an array being initialized using a dummy function in a subroutine.

Here is the code that I used.

module int_adaptive
implicit none
contains

subroutine RK4nd(f_x, norder, nvar, x_0, t_0, t_f, x_out, t_out)
    INTERFACE
        FUNCTION f_x (x_0, t_0)
            integer, parameter  :: dp=kind(0.d0)
            REAL(dp), INTENT(IN)    :: x_0(2, 3), t_0
            REAL, intent(out)   :: f_x(2,3)
        END FUNCTION f_x
    END INTERFACE


    integer, parameter  :: dp=kind(0.d0)
    integer             :: norder, i, nvar
    external            :: f_x
    real(dp), intent(in):: x_0(norder, nvar), t_0, t_f
    real(dp)            :: x_out(norder, nvar), t_out, k1(norder, nvar), y1(norder, nvar), y(norder, nvar)
    real(dp)            :: k2(norder, nvar), k3(norder, nvar), k4(norder, nvar),&
                           y2(norder, nvar), y3(norder, nvar), y4(norder, nvar)!, f_x(norder, nvar)
    real(dp)            :: h_min, h_0, h_max, tol, err_est, h, t

    if (h_0<h) then
        h = h_0
    end if
    if ((t_f - t_0) < 0.0) then
        h = -h
    end if
    t = t_0
    y = x_0

    do while (t .NE. t_f)
        k1 = f_x(y, t)
        y1 = y + k1*h/2.0

        k2 = f_x(y1, t+h/2.0)
        y2 = y + k2*h/2.0

        k3 = f_x(y2, t+h/2.0)
        y3 = y + k3*h

        k4 = f_x(y3, t+h)
        y4 = y + (k1+ 2*k2 + 2*k3 +k4)*h/6.0

        t = t + h
        y = y4

    end do

    x_out = y
    t_out = t

end subroutine
end module

I could define a standard integrator with fixed function input inside a module and call it by name directly, but since I deal with ODE's of various orders it will be complicated each time I am changing the order I'll need to update this module as well.

francescalus
  • 30,576
  • 16
  • 61
  • 96
Astroynamicist
  • 193
  • 2
  • 11

2 Answers2

5

You have two distinct problems here.

The first is that your specification of the explicit interface for the function f_x is incorrect: the function result f_x must not have the intent(out) attribute. This covers the second error message.

The first error message is related to the later use of

external f_x

This external statement gives f_x the external attribute. Which is fine, because the dummy procedure is an external function. However, you've also given that procedure an explicit interface with the (now corrected) interface block. This interface block also states that the procedure has the external attribute.

Doing this, you've violated a constraint that an entity not be explicitly given the same attribute twice in one scoping block. To resolve this issue, you should remove one of the specifications. Because the function returns an array result, its interface must be explicit in the subroutine where it is referenced, so the interface block must be retained. That is, remove the external statement.

For clarity, you should also remove the commented out real(dp) f_x(norder, nvar) declaration. Such a declaration is erroneous for a function.

Community
  • 1
  • 1
francescalus
  • 30,576
  • 16
  • 61
  • 96
  • i got the point on removing the external block, but it was the intent part that did the trick. No compiler error. Need to test it. Thanks a lot – Astroynamicist Jan 17 '17 at 17:04
2

As francescalus already pointed out, a function should be either explicitly INTERFACEd, or declared as EXTERNAL - not both. Using EXTERNAL is old-fashion Fortran, since function interfacing was introduced in Fortran 90 to replace the need of EXTERNAL, which is a rather obsolete feature, still valid for backwards compatibility.

In addition, a function result cannot have an INTENT attribute. A function is meant to return a result by "its own name", as in y=f_x(...). Therefore, you should either declare f_x itself as REAL (without any INTENT attribute), or, even better, use the RESULT attribute, both in function declaration and its interface. Using the RESULT attribute is actually necessary in recursive functions, but it is a good practice to do so even if the function is not recursive. So, the interface for your function should be written as

INTERFACE
    FUNCTION f_x (x_0, t_0) RESULT(res)
        INTEGER, PARAMETER :: dp=kind(0.d0)
        REAL(kind=dp), DIMENSION(2,3), INTENT(IN) :: x_0
        REAL(kind=dp), INTENT(IN) :: t_0
        REAL, DIMENSION(2,3) :: res
    END FUNCTION f_x
END INTERFACE

In your actual f_x implementation, you should use res (or whatever you want to call the RESULT) as the variable holding the result this function should return.

Note that I also did a few modifications which are not strictly necessary, but highly recommended: declaring x_0 as you do, REAL(dp), INTENT(IN) :: x_0(2, 3) is also old-fashion Fortran (the DIMENSION attribute was introduced to make things clearer). Also some Fortran implementations, especially the F-standard, do not accept REAL(dp) and need REAL(kind=dp) instead.

Pap
  • 448
  • 1
  • 10
  • 15
  • I won't comment on the style point of preferring `result`, but there is something to say. Using `result` (not an attribute) isn't necessary in the interface block for a recursive function. The function result's name isn't a characteristic (of the function result or the function), and as there's no executable statement to reference the function rather than the function result it isn't needed to have a distinct name. The implementation needs to have a means of distinguishing between them but the names of the function result in the implementation and the interface needn't match, – francescalus Jan 17 '17 at 19:09
  • @francescalus : According to the F-standard, using `result` (which seems a function attribute to me, no matter how it is officially called) is imperative even for non-recursive functions. In this particular case, using `result` might be needed for another reason: the function's actual result is an array, and some compilers won't accept `real, dimension(2,3) :: f_x`. Of course, the names of the function result in the interface and the actual implementation don't need to match, so one could use `result(res)` in the interface and `result(foo)` in the function implementation. – Pap Jan 17 '17 at 19:27
  • @Pap thanks for the info, will look into result with more details, but I thought stating dimension explicitly was earlier practice as you need to declare scalar variables separately, which is why I switched to giving dimensions after stating the variable. Something new I learned. – Astroynamicist Jan 17 '17 at 20:08
  • @Astroynamicist : Indeed, using the DIMENSION attribute means you have to declare scalar variables separately (same for any other arrays of different dimensions). But the code is more clear when different entities are declared separately - I guess that was the reason the Fortran committee introduced the DIMENSION attribute. In any case, you don't have to use it. Your variables declarations are perfectly valid, the important part is how you define `f_x` and its interface. – Pap Jan 19 '17 at 09:26