Problem with scalapack routine pdgemm

[rosenbe2]

I have been developing a code using pdgemm. It seems to work if the c-matrix is initialized to zero, but not when it is initialized to anything else. I have tested the program out on three different machines and I obtain the same results - pdgemm yields incorrect answeres when the c-matrix is not identically zero. dgemm works fine. I have included the fortran source code below. Can anybody tell me what I am doing wrong.

program pmm
integer m,n,mb,nb,ntimer,na
parameter(m=1500,n=m,mb=64,nb=mb,ntimer=3,na=330)
parameter(mc=465,nc=mc)
real, dimension(m,n) :: a,b,c,e,f,s
real la(na,na),lb(na,na),lc(na,na),d(mc,nc)
real*8 times(10),diff(10)

integer iam,nprocs,ictxt,myrow,mycol
integer desca(9),descb(9),descc(9)
integer rsrc,csrc,llda,lldb,lldc,info
integer lm,ln
integer i,j,iloc,jloc

integer ipr,npr,ierr_mpi
include 'mpif.h'

call MPI_INIT(ierr_mpi)
call MPI_COMM_RANK(MPI_COMM_WORLD,ipr,ierr_mpi)
call MPI_COMM_SIZE(MPI_COMM_WORLD,npr,ierr_mpi)

print*,'%pmm, ',ipr,' out of ',npr


fac=1.d-6
do i=1,mc
do j=1,nc
a(i,j)=(10.d0*i+j )*fac
b(i,j)=(a(i,j)+5.d0)*fac
s(i,j)=1.d0 *fac
end do
end do

times=0.d0

call slboot()
c
c...... by hand "c"
c
c=s
call sltimer(1)
if (ipr.eq.0) then
do i=1,mc
do k=1,nc
d(k,i)=a(i,k)
end do
end do
do j=1,nc
do i=1,mc
t=0.d0
do k=1,nc
t=t+d(k,i)*b(k,j)
end do
c(i,j)=c(i,j)+t
end do
end do
end if
call sltimer(1)
c
c...... by lapack "e"
c
e=s
call sltimer(2)
if (ipr.eq.0) call dgemm('n','n',mc,nc,nc,1.d0,a,m,b,m,1.d0,e,m)
call sltimer(2)
c
c...... by scalapack "f"
c
f=s
call MPI_BARRIER(MPI_COMM_WORLD,ierr_mpi)
call sltimer(3)
call blacs_pinfo(iam,nprocs)
nprow=int(sqrt(float(nprocs)))
do while(mod(nprocs,nprow).ne.0)
nprow=nprow-1
end do
npcol=nprocs/nprow
if (iam.eq.0) print'(a,i3,a,i3)','%pmm, processor grid: ',nprow,' x ',npcol

call blacs_get(-1,0,ictxt)
call blacs_gridinit(ictxt,'r',nprow,npcol)
call blacs_gridinfo(ictxt, nprow,npcol,myrow,mycol)

lm=numroc(mc,mb,myrow,0,nprow)
ln=numroc(nc,nb,mycol,0,npcol)

call descinit(desca,mc,nc,mb,nb,0,0,ictxt,na,info)

do iloc=1,lm
do jloc=1,ln
i=indxl2g(iloc,mb,myrow,0,nprow)
j=indxl2g(jloc,nb,mycol,0,npcol)
la(iloc,jloc)=a(i,j)
lb(iloc,jloc)=b(i,j)
lc(iloc,jloc)=f(i,j)
end do
end do

call pdgemm('n','n',mc,nc,nc,1.d0,la,1,1,desca,lb,1,1,desca,1.d0,lc,1,1,desca)

do iloc=1,lm
do jloc=1,ln
i=indxl2g(iloc,mb,myrow,0,nprow)
j=indxl2g(jloc,nb,mycol,0,npcol)
f(i,j)=lc(iloc,jloc)
end do
end do

call dgsum2d(ictxt,'All',' ',mc,nc,f,m,-1,-1)
call sltimer(3)

if (ipr.eq.0) then
diff=0.0
do i=1,mc
do j=1,nc
diff(1)=diff(1)+abs(e(i,j)-c(i,j))/(e(i,j)+c(i,j))*2.0
diff(2)=diff(2)+abs(f(i,j)-c(i,j))/(f(i,j)+c(i,j))*2.0
end do
end do
print*,'%pmm, diff dgemm pdgemm '
print'(a,10e10.3)','%pmm, ',diff(1:ntimer-1)/(mc*nc)
end if

call blacs_barrier(ictxt,'All')
call slcombine(ictxt,'All','>','W',ntimer,1,times)
if (iam.eq.0) then
print*,'%pmm, hand dgemm pdgemm'
print'(a,10f8.3)','%pmm, abs ',times(1:ntimer)
print'(a,10f8.3)','%pmm, rel ',times(1:ntimer)/times(ntimer)
end if

call blacs_gridexit(ictxt)
call blacs_exit(1)

call MPI_BARRIER(MPI_COMM_WORLD,ierr_mpi)
call MPI_Finalize(ierr_mpi)

stop
end