Error: Statement Function is recursive

fortran, fortran90, module, recursion, subroutine

Solution

As is so often the case, `implicit none` is your friend.

Let's take the errors one at a time:

B(i)=B(i)-B(k)*C       !error: Statement Function is recursive

the compiler doesn't recognize `B` in this context; there is no variable declared `B` here(*), so the best it can do is assume that it's a real-valued statement function defining B as a function of I. Statement functions are space saving but confusing ways of defining a function inline, but amongst other things they can't be recursive; here you'd be defining `B(i)` in terms of `B(i)` which obviously fails.

(*) But! you cry. `B` is an array field in my type `LAE`! Yes, but we're not in the context of LAE here; and indeed, in the context of this function there are no variables of type `LAE` to even use the B values of. This is because the procedure was defined `nopass`; you need to have a variable which is the object being operated on which is of class `LAE` so that we can access those fields. That looks something like this:

    type LAE
        !...
    contains
        procedure::RowReduction
    end type LAE

contains
    subroutine RowReduction(self)
        class(LAE), intent(InOut) :: self
        double precision::C
        integer :: i, j, k
        do k= 1, self%N
            do i= k+1, self%N
                if( self%A(k,k) /= 0 ) then
        !....

Note that we have to define self to be of `class(LAE)` rather than type; class is a superset of type and is needed when dealing with extensible objects which include those with (re-)allocatable components. Note too we've added implicit none which would tell you right away that `B` wasn't defined, and thus specified the integer indices i,j, and k.

Once `N`, `A`, and `B` are being correctly referenced as fields of `self`, then most of the rest of the program is correct. Note that you have to `reshape` your `LAE1%A` array:

LAE1%A=reshape((/2,0,0,0,2,0,0,0,2/), (/N, N/))

But otherwise things seem fine.

module LinearSolution
implicit none

    type LAE
        integer::N
        double precision,dimension(:,:),allocatable::A
        double precision,dimension(  :),allocatable::B
    contains
        procedure::RowReduction
    end type LAE

contains
    subroutine RowReduction(self)
        class(LAE), intent(InOut) :: self
        double precision::C
        integer :: i, j, k
        do k= 1, self%N
            do i= k+1, self%N
                if( self%A(k,k) /= 0 ) then
                    C = self%A(i,k) / self%A(k,k)
                    self%B(i) = self%B(i)- self%B(k)*C
                    do j=k+1, self%N
                         self%A(i,j) = self%A(i,j) - self%A(k,j)*C
                    end do
                end if
            end do
        end do

        do k = self%N,1,-1
            do i=k-1,1,-1
                if( self%A(k,k)/=0) then
                    C= self%A(i,k)/ self%A(k,k)
                    self%B(i)= self%B(i)- self%B(k)*C
                end if
            end do
        end do

        do k=1, self%N
            if( self%A(k,k)/=0 ) then
                self%B(k) = self%B(k) / self%A(k,k)
            end if
        end do
    end subroutine RowReduction
end module LinearSolution

program TestLAE
    use LinearSolution
    implicit none

    integer, parameter :: N = 3
    type(LAE)::LAE1
    LAE1%N=N
    allocate(LAE1%B(1:N))
    allocate(LAE1%A(1:N,1:N))
    LAE1%B=(/1,1,1/)
    LAE1%A=reshape((/2,0,0,0,2,0,0,0,2/), (/N, N/))
    call LAE1%RowReduction
    print*, LAE1%B(1),LAE1%B(2),LAE1%B(3)
end program

Running gives:

$ gfortran -o lae lae.f90 -Wall -std=f2003
$ ./lae
  0.50000000000000000       0.50000000000000000       0.50000000000000000     

Problem

Original source