Fix PR libfortran/21926 matmul does not deal with non-packed result

[Thomas Koenig]

On Tue, Jun 07, 2005 at 12:45:18PM +0200, Tobias Schlüter wrote:
> Thomas Koenig wrote:
> > This has been regression-tested on mainline.  OK?
> > 
> > This patch doesn't apply cleanly with 4.0, because of a cast
> > intrduced for cleanup, but the port is straightforward.
> > 
> > OK to apply the equivalent patch once the 4.0 branch reopens?
> 
> Both ok.

Committed to 4.1, thanks.

> I haven't understood how the 43 and 44 in the result came about.
> Can you explain that, please?

program main
  real, dimension(2,2) :: a
  real, dimension(4,2) :: c
  a = reshape((/ 1.0, 1.0, 0.0, 1.0/), shape(a))
  c = 42.

After this, c is, in Fortran order,

42. 42.
42. 42.
42. 42.
42. 42.

c(1:2,1:2) has a lower stride of 1 and an upper stride of 4.

In executing

  c(1:2,1:2) = matmul(a, transpose(a)

the matmul routines zeroed a contiguous amount of memory if the
lower stride was one, (which was in error), so we got


then

 0. 42.
 0. 42.
 0. 42.
 0. 42.

The result was then added to the correct elements, so we had

 1. 43.
 1. 44.
 0. 42.
 0. 42.

instead of the correct

 1.  1.
 1.  2.
42. 42.
42. 42.