Vectors and Matrices

I'm drawing a blank about how to solve this problem. I have a
subroutine with a matrix dummy argument, stored in fortran as a 2D
array with assumed shape. I have a vector stored in fortran as a 1D
array. I want to pass that vector to the subroutine, and tell the
subroutine that there is only one column.

With f77 and storage association, this was a straightforward thing
to do. But with assumed shape arrays, the only way I can think of
doing this is to make a local copy, pass the copy as the actual
argument, and then copy the result back into the vector.

Is there a better way of doing this? One that avoids that
unnecessary allocate-copy-deallocate overhead?

The subroutine is not pure, so I cannot declare it as elemental so
that it works on both 1D and 2D arrays. Also, it doesn't work with
scalars anyway, so this doesn't really make sense.

RESHAPE would work if the dummy array were INTENT(IN), but that is
not the case. The dummy argument needs to be INTENT(INOUT).

Is there a standard-conforming way to get a 2D pointer to point to
the elements of the 1D array? If so, I could then use the 2D
pointer as the actual argument. I remember something like this in
one of the recent language standards (2003 or 2008), but the last I
checked that feature is not yet commonly supported.

Another option would be to declare the vector as a 2D array with one
column. However, this vector is used in multiple routines, so that
just pushes the problem back one level, and I don't want to make
that kind of global change to the code.

$.02 -Ron Shepard

Ron Shepard <ron-shepard@nospam.comcast.net> wrote:

> I'm drawing a blank about how to solve this problem. I have a
> subroutine with a matrix dummy argument, stored in fortran as a 2D
> array with assumed shape. I have a vector stored in fortran as a 1D
> array. I want to pass that vector to the subroutine, and tell the
> subroutine that there is only one column.

> With f77 and storage association, this was a straightforward thing
> to do. But with assumed shape arrays, the only way I can think of
> doing this is to make a local copy, pass the copy as the actual
> argument, and then copy the result back into the vector.

> Is there a better way of doing this? One that avoids that
> unnecessary allocate-copy-deallocate overhead?

I would probably go for the allocate/copy/deallocate unless it was
called an especially large number of times with different sizes,
or was especially large.

My first thought, though, is to pass a C_PTR from C_LOC to
C_F_PTR to generate an appropriate pointer to a 2D array.
That does assume that the vector is in contiguous storage, though
that would be usual.

I believe that requires only F2003 features, and that many compilers
now implement those.

-- glen

Ron Shepard <ron-shepard@NOSPAM.comcast.net> wrote:

> I'm drawing a blank about how to solve this problem. I have a
> subroutine with a matrix dummy argument, stored in fortran as a 2D
> array with assumed shape. I have a vector stored in fortran as a 1D
> array. I want to pass that vector to the subroutine, and tell the
> subroutine that there is only one column.
>
> With f77 and storage association, this was a straightforward thing
> to do. But with assumed shape arrays, the only way I can think of
> doing this is to make a local copy, pass the copy as the actual
> argument, and then copy the result back into the vector.
>
> Is there a better way of doing this? One that avoids that
> unnecessary allocate-copy-deallocate overhead?
....
> Is there a standard-conforming way to get a 2D pointer to point to
> the elements of the 1D array? If so, I could then use the 2D
> pointer as the actual argument. I remember something like this in
> one of the recent language standards (2003 or 2008), but the last I
> checked that feature is not yet commonly supported.

Yep. But you are right that it is an f2003 feature.

P(1,n) => T

I think there might be possibilities of using a subroutine to hack such
a thing. Make the subroutine have an explicit-shape dummy argument with
dimensions (1,n). You can still use sequence association in f90/f95 -
just not with an assumed shape dummy. So make a different subroutine
with an explicit-shape dummy just for the purpose of getting a pointer
to it. Something like

subroutine make_pointer(p,t,n)
real, target :: t(1,n)
real, pointer :: p(:,
p => t
end subroutine make_pointer

But whether the standard guarantees this to work is too subtle a
question for me at this hour (or perhaps even at other hours). The
treatment of target dummy arguments in the standard is quite messy and
one of the parts of the standard that I've never found easy. Anyway,
this *MIGHT* be ok. It might also depend on details of the actual
argument for T; I rather suspect it does. Is T itself an assumed-shape
one-dimensional dummy at the next higher level? If so, I bet it won't
work, as you might end up with a pointer to a temporary copy. But I'm
basing that just on thinking about what an implementation might do
instead of on reading the messy words about it in the standard. There's
supposed to be a correlation between those two viewpoints, but it can be
a tricky one.

> Another option would be to declare the vector as a 2D array with one
> column. However, this vector is used in multiple routines, so that
> just pushes the problem back one level, and I don't want to make
> that kind of global change to the code.

I suppose I should also mention the possibility of equivalence. But
that's pretty limitting and so likely not to be a good choice.

--
Richard Maine | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle | -- Mark Twain

glen herrmannsfeldt <gah@ugcs.caltech.edu> wrote:

> Ron Shepard <ron-shepard@nospam.comcast.net> wrote:
>
> > I'm drawing a blank about how to solve this problem. I have a
> > subroutine with a matrix dummy argument, stored in fortran as a 2D
> > array with assumed shape. I have a vector stored in fortran as a 1D
> > array. I want to pass that vector to the subroutine, and tell the
> > subroutine that there is only one column.

> My first thought, though, is to pass a C_PTR from C_LOC to
> C_F_PTR to generate an appropriate pointer to a 2D array.
> That does assume that the vector is in contiguous storage, though
> that would be usual.
>
> I believe that requires only F2003 features, and that many compilers
> now implement those.

Oh yes. I forgot about that one. Of course, the pointer remapping that
Ron said didn't yet work with all his compilers is also f2003. However,
compilers were pretty fast to implement much of the C interop stuff in
f2003, so availability of that solution might still be better.

--
Richard Maine | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle | -- Mark Twain

On Apr 14, 4:59*pm, Ron Shepard <ron-shep...@NOSPAM.comcast.net>
wrote:
> I'm drawing a blank about how to solve this problem. *I have a
> subroutine with a matrix dummy argument, stored in fortran as a 2D
> array with assumed shape. *I have a vector stored in fortran as a 1D
> array. I want to pass that vector to the subroutine, and tell the
> subroutine that there is only one column.
>
> With f77 and storage association, this was a straightforward thing
> to do. *But with assumed shape arrays, the only way I can think of
> doing this is to make a local copy, pass the copy as the actual
> argument, and then copy the result back into the vector.
>
> Is there a better way of doing this? *One that avoids that
> unnecessary allocate-copy-deallocate overhead?
>
> The subroutine is not pure, so I cannot declare it as elemental so
> that it works on both 1D and 2D arrays. Also, it doesn't work with
> scalars anyway, so this doesn't really make sense.
>
> RESHAPE would work if the dummy array were INTENT(IN), but that is
> not the case. *The dummy argument needs to be INTENT(INOUT).
>
> Is there a standard-conforming way to get a 2D pointer to point to
> the elements of the 1D array? *If so, I could then use the 2D
> pointer as the actual argument. *I remember something like this in
> one of the recent language standards (2003 or 2008), but the last I
> checked that feature is not yet commonly supported.
>
> Another option would be to declare the vector as a 2D array with one
> column. *However, this vector is used in multiple routines, so that
> just pushes the problem back one level, and I don't want to make
> that kind of global change to the code.

You could use generic subroutines: one that accepts a vector,
the other accepts a matrix.

Alternatively, you could include two arguments in the calls,
one optional, of course.

In article <1kiir0w.1cp5poqhf0p6oN%nospam@see.signature>,
nospam@see.signature (Richard Maine) wrote:

> glen herrmannsfeldt <gah@ugcs.caltech.edu> wrote:
>
> > Ron Shepard <ron-shepard@nospam.comcast.net> wrote:
> >
> > > I'm drawing a blank about how to solve this problem. I have a
> > > subroutine with a matrix dummy argument, stored in fortran as a 2D
> > > array with assumed shape. I have a vector stored in fortran as a 1D
> > > array. I want to pass that vector to the subroutine, and tell the
> > > subroutine that there is only one column.
>
> > My first thought, though, is to pass a C_PTR from C_LOC to
> > C_F_PTR to generate an appropriate pointer to a 2D array.
> > That does assume that the vector is in contiguous storage, though
> > that would be usual.
> >
> > I believe that requires only F2003 features, and that many compilers
> > now implement those.
>
> Oh yes. I forgot about that one. Of course, the pointer remapping that
> Ron said didn't yet work with all his compilers is also f2003. However,
> compilers were pretty fast to implement much of the C interop stuff in
> f2003, so availability of that solution might still be better.

Thanks to everyone for the comments. I think your explicit shape
trick will work:

subroutine make_pointer(p,t,n)
real, target :: t(1,n)
real, pointer :: p(:,
p => t
end subroutine make_pointer

I need t(n,1) column vectors rather than row vectors in the matrix,
but that should not matter. I have done this in the past, but I had
forgotten that the target needs to be explicit shape to get around
the argument rank checking.

I agree about your comments about whether this is strickly standard
conforming. I remember testing this with LOC() on a bunch of
compilers about 10 years ago, and it did seem to work correctly
(i.e. it avoided the copy and referenced the correct elements), but
I also could not verify that it was required to work for all
possible combinations of targets (strided input, derived types as
well as integers and reals, etc.).

F90 and later provides a simple and straightforward way to access a
column of a matrix as a scalar, matrix(:,k), but it does not provide
a similarly straightforward way to do the opposite. When combined
with the TKR matching of actual and dummy arguments, this causes the
kind of problem I had here.

This all results from using a set of library routines that work with
matrices, and matching them to another set of library routines that
work with vectors. Going one way is simple, going the other way
requires some trickery. I think I can localize the above to just a
few interface routines, so I think this is quite acceptable
esthetically.

Anyway, thanks again to everyone for the comments.

$.02 -Ron Shepard

Richard Maine <nospam@see.signature> wrote:
> Ron Shepard <ron-shepard@NOSPAM.comcast.net> wrote:

>> I'm drawing a blank about how to solve this problem. I have a
>> subroutine with a matrix dummy argument, stored in fortran as a 2D
>> array with assumed shape. I have a vector stored in fortran as a 1D
>> array. I want to pass that vector to the subroutine, and tell the
>> subroutine that there is only one column.

(snip)

> I think there might be possibilities of using a subroutine to hack such
> a thing. Make the subroutine have an explicit-shape dummy argument with
> dimensions (1,n). You can still use sequence association in f90/f95 -
> just not with an assumed shape dummy. So make a different subroutine
> with an explicit-shape dummy just for the purpose of getting a pointer
> to it. Something like

> subroutine make_pointer(p,t,n)
> real, target :: t(1,n)
> real, pointer :: p(:,
> p => t
> end subroutine make_pointer

> But whether the standard guarantees this to work is too subtle a
> question for me at this hour (or perhaps even at other hours). The
> treatment of target dummy arguments in the standard is quite messy and
> one of the parts of the standard that I've never found easy.

Easier, though, is to just call the original subroutine with the
new dummy argument.

subroutine onedee(t,n)
real t(1,n)
call twodee(t)
return
end

(Either t(1,n), or t(n,1), as appropriate.)

For extra complication, have a subroutine that accepts assumed
shape and then calls that one:

subroutine onedeeas(t)
real t(
call onedee(t,size(t))
return
end

-- glen

Ron Shepard <ron-shepard@NOSPAM.comcast.net> wrote:

> I need t(n,1) column vectors rather than row vectors in the matrix,
> but that should not matter. I have done this in the past, but I had
> forgotten that the target needs to be explicit shape to get around
> the argument rank checking.
....
> F90 and later provides a simple and straightforward way to access a
> column of a matrix as a scalar, matrix(:,k), but it does not provide
> a similarly straightforward way to do the opposite. When combined
> with the TKR matching of actual and dummy arguments, this causes the
> kind of problem I had here.

Note that TKR matching is *NOT* required in general. In particular,
that's not your problem in this case (at least as described). TKR
matching is for generics.

Assumed shape needs rank matching, but that's for a different reason,
which is pretty much unrelated to TKR. It needs rank matching because it
is assuming a shape, so there has to be a shape of the appropriate rank
to assume.

--
Richard Maine | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle | -- Mark Twain