How are these arrays being used in this Fortran algorithm?

Question

I'm writing some subroutines in Fortran90 to perform some numerical computations. However, as part of this, I need to include some codes from the netlib templates library that are written in Fortran77. I'm having a hard time getting them to work - specifically understanding the usage of arrays.

For instance, I need to use a subroutine called GMRES. Here are the arguments:

  SUBROUTINE GMRES( N, B, X, RESTRT, WORK, LDW, WORK2, LDW2, 
 $                  ITER, RESID, MATVEC, PSOLVE, INFO )

*     .. Scalar Arguments ..
  INTEGER            N, RESTRT, LDW, LDW2, ITER, INFO
  DOUBLE PRECISION   RESID
*     ..
*     .. Array Arguments ..
  DOUBLE PRECISION   B( * ), X( * ), WORK( * ), WORK2( * )
*     ..
*     .. Function Arguments ..
*
  EXTERNAL           MATVEC, PSOLVE
*
*
*  Arguments
*  =========
*
*  N       (input) INTEGER. 
*          On entry, the dimension of the matrix.
*          Unchanged on exit.
* 
*  B       (input) DOUBLE PRECISION array, dimension N.
*          On entry, right hand side vector B.
*          Unchanged on exit.
*
*  X       (input/output) DOUBLE PRECISION array, dimension N.
*          On input, the initial guess; on exit, the iterated       solution.
*
*  RESTRT  (input) INTEGER
*          Restart parameter, <= N. This parameter controls the amount
*          of memory required for matrix WORK2.
*
*  WORK    (workspace) DOUBLE PRECISION array, dimension (LDW,5).
*          Note that if the initial guess is the zero vector, then 
*          storing the initial residual is not necessary.
*
*  LDW     (input) INTEGER
*          The leading dimension of the array WORK. LDW >= max(1,N).
*
*  WORK2   (workspace) DOUBLE PRECISION array, dimension     (LDW2,2*RESTRT+2).
*          This workspace is used for constructing and storing the
*          upper Hessenberg matrix. The two extra columns are used to
*          store the Givens rotation matrices.
*
*  LDW2    (input) INTEGER
*          The leading dimension of the array WORK2.
*          LDW2 >= max(1,RESTRT).
*
*  ITER    (input/output) INTEGER
*          On input, the maximum iterations to be performed.
*          On output, actual number of iterations performed.
*
*  RESID   (input/output) DOUBLE PRECISION
*          On input, the allowable error tolerance.
*          On ouput, the norm of the residual vector if solution
*          approximated to tolerance, otherwise reset to input
*          tolerance.
*
*  INFO    (output) INTEGER
*          =  0:  successful exit
*          =  1:  maximum number of iterations performed;
*                 convergence not achieved.
*
*  MATVEC  (external subroutine)
*          The user must provide a subroutine to perform the
*          matrix-vector product A*x = y.
*          Vector x must remain unchanged. The solution is
*          over-written on vector y.
*
*          The call is:
*
*             CALL MATVEC( X, Y )

Notice how the arrays are defined WORK( * ), WORK2( * ). So to my mind these are 1D arrays of arbitrary length. But then in the argument description, it seems to be suggesting that they are 2D arrays, matrices, with dimension WORK(LDW, 5). So are they 1D or 2D?

Also, in the GMRES algorithm, these WORK arrays are used like this:

CALL MATVEC(SCLR1, WORK(NDX1), SCLR2, WORK(NDX2))

So if the WORK arrays are 2D, wouldn't the above give a rank mismatch? What does it mean to access a 2D array with just one index like that? Should I define the WORK arrays as 2D or 1D?

Edit

The GMRES routine requires a matvec routine to be supplied. In the GMRES code it is being called like

CALL MATVEC(SCLR1, X, SCLR2, WORK(NDX2))

and also like

CALL MATVEC(SCLR1, WORK(NDX1), SCLR2, WORK(NDX2))

My subroutine MATVEC that I'm trying to supply looks like:

subroutine matvec(alpha, x, beta, y)

    real(dp), intent(in) :: alpha, beta
    real(dp), dimension(*), intent(in) :: x
    real(dp), dimension(*), intent(inout) :: y
    real(dp), dimension(*,*) :: A
    integer :: n

    call make_Jac(n,A)
    call dgemv('notranspose', n, n, alpha, A, n, x, 1, beta, y, 1)

end subroutine matvec

Where make_Jac returns my matrix for the problem I'm working on, along with its dimension n. The blas routine dgemv then handles the matrix-vector product.

Answer 1

WORK( * ) declares an assumed size array, which can be two-dimensional. See here.

The compiler will not complain if you feed a one-dimensional array to the subroutine, but weird things (up to and including a segmentation fault) might happen.

Better use arrays matching the specifications.

Answer 2

The same Fortran array can be managed as one dimensional, two dimensional, etc. It is stored in contiguous memory in any case.

Let's say you have

double precision x(3, 2)

call somefunc (x)

This can be accessed, inside somefunc, as y (6).

The array elements are stored in "column major order" which means

x(1, 1) is y(1)
x(2, 1) is y(2)
x(3, 1) is y(3)
x(1, 2) is y(4)
x(2, 2) is y(5)
x(3, 2) is y(6)

As long as the function knows the limits of each dimension, it can calculate linear access by simple arithmetic. Alas, this "relaxed type" is also a frequent source of bugs.

Answer 3

Further to the answer by wallyk the dummy argument work(*) is a rank-1 assumed size array. With such an array

The rank and extents may differ for the effective and dummy arguments; only the size of the effective argument is assumed by the dummy argument.

This means that it is quite acceptable for such a structure as

double precision work(LDW,5)
call GMRES(..., work, ...)
! ...

end

subroutine GMRES(..., work, ...)
  double precision work(*)
  ! ...
end subroutine

Indeed, with the dummy argument as a rank-1 array it isn't allowed to reference it as a rank-2 array. A rank-2 assumed size array would look something like

double precision work(LDW,*)

where then, of course, work(ndx1) would be bad.

Coming to the comment by roygvib, later on in the Netlib source code there is the line

call GMRESREVCOM(..., work, ...)

where that subroutine has the dummy argument

double precision work(LDW,*)

There is probably, then, an expectation that the user of the code will provide initially a rank-2 actual argument.

What all of this means is that it doesn't matter what rank the actual argument passed to work in GMRES as long as it has at least LDW*5 elements. I'd be careful calling the dummy argument as being of "arbitrary length", though, as referencing work(LDW*5+1) (according to my first example) would be wrong. The size of the dummy array is exactly the size of the passed array.

The later call to matvec is not troublesome for yet another reason. This subroutine has four arguments, the first and third of which are scalar. The second and fourth are again assumed-size arrays of rank-1. We've established that assumed-size arrays don't care about the rank of the effective/actual argument, but you're likely wondering why you can pass the scalar argument work(ndx1) to this rank-1 array.

The answer to that is something called sequence association. This means that when your actual argument is an array element designator and the dummy argument is an array dummy argument then that dummy argument is argument associated with a sequence of elements from the actual argument, starting with the element element designated.

That is, you have a rank-1 array like [work(ndx1), work(ndx1+1), ...] as your array x in matvec.

This is all fine, as long as you don't attempt to reference beyond the extent of your actual argument.