Re: TRANSPOSE and MATMUL interaction

Peter wrote:
On Aug 30, 2:50 pm, Paul van Delst <Paul.vanDe...@xxxxxxxx> wrote:
What platform and compiler are you using? How are the results different? You wouldn't be
doing this on an AIX machine with xlf by any chance, would you?

Windows XP Pro/Visual Studio 2005/Intel Visual Fortran 9.1. The result
I believe correct is 1.2e-3, and the incorrect result is 6e-304, which
is the same as an uninitialized variable, and wasn't changed by the
MATMUL assignment.


Ah. That's about as from AIX xlf95 as you can get I think :o)

If it's easy for you to do, what does the output from the test program pasted below look like for your system?

It's a test program for a system I use (AIX xlf95 v10.1) that replicates a compiler bug related to matmul. On our system I've found that the following "combinations" work,

b = matmul(A,x)
b(1:ni) = matmul(A,x)
b(1:ni) = matmul(A,x(1:nj))
b = matmul(A,x(1:nj))
b = matmul(A(1:ni,1:nj),x(1:nj))

but if you do either of these:

b(1:ni) = matmul(A(1:ni,1:nj),x)
b(1:ni) = matmul(A(1:ni,1:nj),x(1:nj))

our results in "b" will be wrong.

IBM found the bug and, apparently, fixed it; but the patch still hasn't trickled down to our system yet. :o(


<-----code begins here------>
program tmatmul

implicit none
integer, parameter :: nX=4, nY=2
real, dimension(nX,nX) :: A
real, dimension(nX) :: x, b
integer :: i

write(*,'(/5x,"The following produces incorrect results...")')
A = 0.0
x = 0.0
b = 0.0
A(1:nX,1:nY) = RESHAPE( (/(real(i)/314.15927,i=1,nX*nY)/), (/nX,nY/) )
x(1:nX) = (/(real(i),i=1,nX)/)
call multiplybusted1(A(1:nX,1:nY),x(1:nY),b(1:nX))

write(*,'(/5x,"The following produces correct results...")')
A = 0.0
x = 0.0
b = 0.0
A(1:nX,1:nY) = RESHAPE( (/(real(i)/314.15927,i=1,nX*nY)/), (/nX,nY/) )
x(1:nX) = (/(real(i),i=1,nX)/)
call multiply(A(1:nX,1:nY),x(1:nY),b(1:nX))

write(*,'(/5x,"Calling the busted routine again...")')
A = 0.0
x = 0.0
b = 0.0
A(1:nX,1:nY) = RESHAPE( (/(real(i)/314.15927,i=1,nX*nY)/), (/nX,nY/) )
x(1:nX) = (/(real(i),i=1,nX)/)
call multiplybusted1(A(1:nX,1:nY),x(1:nY),b(1:nX))

contains

subroutine multiply(A,x,b)
real, dimension(:,:), intent(in) :: A
real, dimension(:), intent(in) :: x
real, dimension(:), intent(out) :: b
integer :: i, j, ni, nj
real :: bruteforce

ni=size(A,dim=1)
nj=size(A,dim=2)

b = matmul(A,x)

do i=1,ni
bruteforce=0.
do j=1,nj
bruteforce=bruteforce+(A(i,j)*x(j))
end do
print *, '===',i,' Matmul result:',b(i),' Manual mult:',bruteforce
end do
end subroutine multiply

subroutine multiplybusted1(A,x,b)
real, dimension(:,:), intent(in) :: A
real, dimension(:), intent(in) :: x
real, dimension(:), intent(out) :: b
integer :: i, j, ni, nj
real :: bruteforce

ni=size(A,dim=1)
nj=size(A,dim=2)

b(1:ni) = matmul(A(1:ni,1:nj),x(1:nj))

do i=1,ni
bruteforce=0.
do j=1,nj
bruteforce=bruteforce+(A(i,j)*x(j))
end do
print *, '===',i,' Matmul result:',b(i),' Manual mult:',bruteforce
end do
end subroutine multiplybusted1

end program tmatmul
.