How to output 2 or more arrays in a fortran's function?

Question

I am writing a program which computes the LU decomposition of a matrix, with partial pivoting, and I would like the function to output several (2 or 3) matrices without running the program several times to output each one individually, which is a waste of time since it gets me everything I want in one run. Is there a way of doing this? For example, here is my function using Doolittle's algorithm, for square matrix which don't need pivoting. I want my output to be matrix l and u at once, but I know no means of doing that.

function lu_d(aa) result(l)

real, dimension (:,:) :: aa !input matrix
real, dimension (size(aa,1), size(aa,2)) :: a !keeping input variable intact
real, dimension (size(a,1), size(a,2)) :: l , u !lower and upper matrices
integer :: i,j,k !index
real :: s !auxiliar variable

a=aa

do j=1 , size(a,2)
  u(1,j)=a(1,j)
end do

l(1,1)=1

do j=2, size(a,2)
  l(1,j)=0
end do

do i=2, size(a,1)

  l(i,1)=a(i,1)/u(1,1)
  u(i,1)=0

  do j=2, i-1

    s=0
    u(i,j)=0

    do k=1, j-1
      s=s+l(i,k)*u(k,j)
    end do

    l(i,j)=(a(i,j)-s)/u(j,j)

  end do

  l(i,i)=1

  do j=i, size(a,2)

    s=0
    l(i,j)=0

    do k=1, i-1
      s=s+l(i,k)*u(k,j)
    end do

    u(i,j)=a(i,j)-s

  end do

end do

end function

Always use tag [tag:fortran]. Add a specific version where necessary to distinguish. Only a very small number of people follow [tag:fortran95]. — Vladimir F Героям слава, Sep 12 '16 at 09:26
Of course there is a way to do this sort of things, but you should show some effort of yours. How does your code look like? Which problems have you find? — Vladimir F Героям слава, Sep 12 '16 at 09:27
I look up for some information online, but I haven't found anything about it. I will edit the post, and add my code! Just a moment. — António Carneiro, Sep 12 '16 at 09:39
Oh, by output you mean return as a function result? I can't see any array `m` in your code so I am still not sure what you want exactly, but I think it will be easier with a subroutine instead of the function. — Vladimir F Героям слава, Sep 12 '16 at 10:37
sorr, matrix u, fixed it. So, is it impossible to have 2 results for one function? Something like a "vector of matrices". — António Carneiro, Sep 12 '16 at 10:41

score 1 · Accepted Answer · answered Sep 12 '16 at 10:49

You could switch from using a function to using a subroutine. This way you can output values for multiple arrays in the arguments list. Additionally using the INTENT definition when declaring variables in the subroutine, e.g.:

REAL,INTENT(IN)::a declares a and does not allow its values to be altered inside the subroutine/function

REAL,INTENT(OUT)::b declares b and disregards any values it has coming into the subroutine/function

REAL,INTENT(INOUT)::c this is the case by default, if you don't write anything.

I will assume you need the output to be l and u (rather than m), in which case the structure would look something like the one below. Note that l and m should either be declared in the main program and their size defined with respect to aa (as in the first case shown below) OR declared with an allocatable size in the main program, passed to the subroutine without being allocated and allocated within the subroutine (second example). The latter may require you to put the subroutine in a module so that the interfaces are handled properly.

First example:

SUBROUTINE lu_d(aa,l,m)
implicit none
real,intent(in):: a(:,:)
real,intent(out):: l(:,:), m(:,:)
integer:: i,j,k
real:: s

<operations>

RETURN
END SUBROUTINE lud_d

Second example:

SUBROUTINE lu_d(aa,l,m)
implicit none
real,intent(in):: a(:,:)
real,allocatable,intent(out):: l(:,:), m(:,:)
integer:: i,j,k,size_a1,size_a2
real:: s

size_a1=size(aa,1)
size_a2=size(aa,2)
allocate( l(size_a1,size_a2), m(size_a1,size_a2))

<operations>

RETURN
END SUBROUTINE lud_d

Surprisingly I do not do this... But can a function return a structure? If so then it could return a structure with the two arrays, but the subroutine seems easier... Which is probably why I do these things that way. — Holmz, Sep 12 '16 at 11:29
That is a fair point, technically you could also output a matrix that holds `l` and `m` with the function approach, but then you would have to separate them afterwards. Declaring the result of the function as e.g. `A(size(aa,1),size(aa,2),2)` and then having `A(:,:,1)=l; A(:,:,2)=m` should do it, but you can see the inconvenience. — ptev, Sep 12 '16 at 11:37
@Holmz You can do that but without anonymous structures and without syntactic sugar for easy unpacking that is very inconvenient in Fortran. — Vladimir F Героям слава, Sep 12 '16 at 14:14
@ptev note that the RETURN is superfluous and may lead some to thing that it is necessary to have RETURN and STOP everywhere. I would definitely not use the RETURN here. — Vladimir F Героям слава, Sep 12 '16 at 14:15
Noted - need to no some research about the proper use myself as I have inherited a number of codes which have it. — ptev, Sep 12 '16 at 15:32

How to output 2 or more arrays in a fortran's function?

1 Answers1