How are these arrays used in this Fortran algorithm?
I am writing some routines in Fortran90 to do some numerical computation. However, as part of this, I need to include some code from the netlib template library that is written in Fortran77. I am having a hard time getting them to work, especially understanding the use of arrays.
For example, I need to use the GMRES subroutine. 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 )
Note how arrays are defined WORK (*), WORK2 (*). Therefore, in my opinion, these are 1D arrays of arbitrary length. But then in the description of the argument it seems that they are 2D arrays, matrices, with WORK size (LDW, 5). So are they 1D or 2D?
In addition, in the GMRES algorithm, these WORK arrays are used as follows:
CALL MATVEC(SCLR1, WORK(NDX1), SCLR2, WORK(NDX2))
So, if the WORK arrays are 2D, doesn't the above rank mismatch mean? What does it mean to access a two-dimensional array with only one index? Should I define WORK arrays as 2D or 1D?
Edit
The GMRES routine requires the matvec routine. In the GMRES code, it is referred to as
CALL MATVEC(SCLR1, X, SCLR2, WORK(NDX2))
and also how
CALL MATVEC(SCLR1, WORK(NDX1), SCLR2, WORK(NDX2))
My MATVEC routine that I am trying to provide looks like this:
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 is make_Jac
returning my matrix for the problem I am working on, along with its size n. The blas routine then dgemv
processes the matrix-vector product.
source to share
The same Fortran array can be manipulated as one dimensional, two dimensional, etc. It is stored in contiguous memory anyway.
Let's say you have
double precision x(3, 2)
call somefunc (x)
It can be accessed inside somefunc like y (6).
Array elements are stored in "top column 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 compute linear access using simple arithmetic. Alas, this "relaxed type" is also a common source of errors.
source to share
In addition to wallyk's answer, the dummy argument work(*)
is an assumed sized array of rank 1. With an array like this
Rank and extents can be different for effective and bogus arguments; only the size of the effective argument is accepted by the dummy argument.
This means that it is perfectly acceptable for a structure like
double precision work(LDW,5)
call GMRES(..., work, ...)
! ...
end
subroutine GMRES(..., work, ...)
double precision work(*)
! ...
end subroutine
Indeed, with a dummy argument as a rank-1 array, it is not allowed to refer to it as a rank-2 array. The ranged size of the array will look something like this:
double precision work(LDW,*)
where then, of course, work(ndx1)
it would be bad.
Coming to the comment from roygvib, the line later appears in the Netlib source code
call GMRESREVCOM(..., work, ...)
where this subroutine has a dummy argument
double precision work(LDW,*)
There is probably a chance that the user of the code will first provide an actual argument of rank 2.
What all this means is that it doesn't matter what rank the actual argument is passed work
in GMRES
, as long as it has at least no LDW*5
elements. I would be careful to label the dummy argument as "arbitrary length", although work(LDW*5+1)
it would be wrong as a reference (according to my first example). The size of the dummy array is exactly the size of the passed array.
The later call is matvec
not a concern for another reason. This subroutine has four arguments, the first and third of which are scalar. The second and fourth are again arrays of estimated rank-1 size. We've established that arrays of assumed sizes don't care about the rank of the effective / actual argument, but you're probably wondering why you could pass a scalar argument to work(ndx1)
this array of rank-1.
The answer to this question is called sequence association. This means that when your actual argument is a pointer to an array element and the dummy argument is an array argument argument, then that dummy argument is an argument associated with the sequence of elements from the actual argument, starting at the specified element of the element.
That is, you have an array of rank-1 [work(ndx1), work(ndx1+1), ...]
, like your array x
in matvec
.
This is all fine as long as you are not trying to reference outside of your actual argument.
source to share