Re: TRANSPOSE and MATMUL interaction

Just for ino, the Matmul and Manual results, and all three sets of results
are all identical using Windows version of G95, Linux version of Intel Fortran 9.1
and PathScale 3.0 and AIX 5.2 xlf90 version 8.1.0.

Skip

On Thu, 30 Aug 2007 19:38:40 -0400, Paul van Delst <Paul.vanDelst@xxxxxxxx> wrote:

-|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

.