Re: Matrix multiplication using subarrays?

Peter <sherwood@xxxxxxxxxxxxxx> wrote:

I'm having difficulty with matrix multiplication, which I think is due
to a misunderstanding regarding subarrays. I have a matrix A and
vector b declared as

REAL, DIMENSION(mMax,nMax) :: A
REAL, DIMENSION(nMax) :: b

mMax and nMax are the maximum dimensions, but normally the actual
dimensions m and n will be less. The actual data is put into elements
of A with first subscript 1,2,...,m and second subscript 1,2,...,n.

Then, I want to do a matrix multiplication of m x n A and n b, so I
write

MATMUL(A(1:m,1:n), b(1:n))

This does not give the result I expect, apparently because MATMUL
doesn't know where to find the non-contiguous elements of A that are
in use.

How should this kind of array handling be done?

I think you are misdiagnosing the problem. Since you haven't given
enough data for us to do the diagnosis, there isn't much help I can
offer.

Rule number 1 of asking for debugging help is to show raw data rather
than your conclusions. (Well, ok, that's rule number 2. Rule number 1 is
to show at least something. We do get people asking things like "My
program doesn't work. Why?" and providing no more data than that.) In
your case, you have given us your conclusion, namely "apparently because
MATMUL doesn't know where to find the non-contiguous elements". That
conclusion is unlikely to be correct. You haven't given us the raw data
from which we can make our own conclusion. Note that telling us what
values A and B have and what result MATMUL gives would not count as raw
data either. That would be processed data telling us what values you
think A and B have and what values you think MATMUL gave. Odds are
reasonably high that what you thought about them isn't so.

What you have shown should work fine. If it doesn't give the results you
expect, then you'll need to show us a sufficiently complete working
program so that we can see why not. The obvious possibilities include
such things as

1. Your expectation might be wrong. Ok, odds are that isn't it, butit is
at least a case to check for.

2. A and B don't have the values you thought they had. There are *LOTS*
of ways that could happen. You don't show any data at all on that, so I
can't guess.

3. On the other end, perhaps MATMUL is giving the right answer, but you
aren't correctly using the result of matmul, so that whatever you end up
looking at isn't actually the result of matmul. Hmm. Come to think of
it, this is a very real possibility. Real enough that I think I'll make
it my number 1 guess. That's an awfully weak guess; I won't put much
confidence in it. But I'll go out on a limb and make it my guess. It
would be easy to mess up and your question doesn't even acknowledge it
as an issue, which makes it seem likely that you didn't think of it.

The matmul that you showed should give a result of rank 1 and size m.
What are you doing with it? You can't just have matmul as a standalone
statement; it has to appear in some context. As one example, if you do

x = MATMUL(A(1:m,1:n), b(1:n))

then x had better be of rank 1 and have size m. If it has, for example,
size mMax instead of n, that won't work (unless mMax happens to equal
m). You'd need to instead have x(1:m) on the left-hand side.

4. Infinite numbers of other things.

As an aside, since you are presumably using f90/f95 anyway (given the
array slice notation and the MATMUL intrinsic), I often find it more
convenient to use an allocatable array for cases like this. Then you
don't have to have the concept of the "used" part of the array being
different from the dimensions. Balancing that, you do have to fiddle
with allocating the array. And it gets messier if you have to start
filling the array before you know how big it will end up. So there are
tradeoffs. Perhaps using the slice might be best in your situation (I
don't have enough data to tell). But I do suggest at least considering
the question.

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

Relevant Pages

  • Re: GNAT compiler switches and optimization
    ... people asked for from a dual processor machine. ... This shows that, for those with dual-processor machines, easy-to-create Ada is faster than easy-to-create FORTRAN. ... But Ada with explicit loops and FORTRAN with matmul are hardly equivalent. ... It looks like a matmul implemented in C with manual array subscripting logic.. ...
    (comp.lang.ada)
  • Re: Intrinsic functions in f90
    ... You can use MATMUL for matrix vector ... The MASK must be conformable with the input array argument, ...
    (comp.lang.fortran)