Page 1 of 1

PSGEMR2D works for certain timesteps then fails

PostPosted: Tue Apr 15, 2008 6:45 pm
by papis
HI ppl,I would really like your help
I use scalapack to solve a linear system with LU factorization for each timestep of my application.
The matrices that i use don't change in size and are allocated in the first time step (deallocated at the last)

The values of the matrices have normal values (not NaN or INF)

The results of the solution are correct till the timestep that i have the problem

The strange thing is that a a certain timestep (it may change if i change the computer that i run the code) the code is stuck at the first psgemr2d (that distributes the main matrix).By watching gkrellm at that time i noticed that,
1) the cpus work around 90%
2) the seem to exchange data constantly(that may cause the cpu load)
3) the program stucks there(cpu's exchange data endlessly)

THe strange thing here is that the code works fine for a number of timesteps and with no apparent reason this happens

I have quoted the code i use to factorise the matrix A0 and then store the results in A0 but in the correct order (because the factorization uses pivoting).

Has anyone noticed similar behaviour with other routines?
Any idea of what is might happenning?

Code: Select all
subroutine lucmp(AO,ndim,Nmax)
Implicit none
integer, parameter :: wp = kind(1.)
Integer ,intent(in)           :: Nmax,ndim
real ,intent(in)::AO(ndim,ndim)
integer ::j,n,IA,JA,info,numroc,mloc_A,nloc_A,mloc_B,block_size,JB,IB,root
integer:: p,ierr,my_rank,nprow,npcol,icontxt,myrow,mycol,nproc,mypnum,i,nloc_b
real,allocatable :: A(:,:),B(:,:)
integer,dimension(9) :: descAO,descA,descb,descbO,descx
integer,dimension(Nmax+200) :: ipvt


block_size=64
call blacs_pinfo(mypnum,nproc)
nprow = nproc             
npcol = 1
CALL SL_INIT( icontxt, NPROW, NPCOL )
CALL BLACS_GRIDINFO( icontxt, NPROW, NPCOL, MYROW, MYCOL )
mloc_A = numroc(Nmax,block_size,myrow,0,nprow)
nloc_A = numroc(Nmax,block_size,mycol,0,npcol)
call descinit(descA,Nmax,Nmax,block_size,block_size,&
0,0,icontxt,mloc_A,info)

!write (*,*) 'nproc',Nmax,myrow,mycol,NPROW,NPCOL,mloc_a,mloc_b


allocate(A(mloc_A,nloc_A))
CALL SL_INIT( root, 1, 1 )
if  (mycol==0.and.myrow==0) then
call descinit(DESCAO,Nmax,Nmax,block_size,block_size,0,0,root,ndim,info)
else
   descAO(1:9) = 0
   descAO(2) = -1
endif
IA=1
JA=1!abs(Nmax*mycol-nloc_A*mycol+1)
IB=1
JB=1
call     PSGEMR2D(Nmax,Nmax, AO, IA, JA, DESCAO, A, IB, JB, DESCA,  &
                         icontxt )
if (myrow==0.and.mycol==0) then
call blacs_gridexit(root)
endif

call blacs_barrier(icontxt,'ALL')
call psgetrf(Nmax,Nmax,A,iA,jA,descA,ipvt,info)
if (info .ne. 0) print*,'pdgetrf info: ',info

call blacs_barrier(icontxt,'ALL')

call psmat( 'N', Nmax, Nmax, A, IA, JA, DESCA, IPVT,INFO)
CALL SL_INIT( root, 1, 1 )
if(myrow==0.and.mycol==0) then
call descinit(DESCAO,Nmax,Nmax,block_size,block_size,0,0,root,ndim,info)
else
   descAO(1:9) = 0
   descAO(2) = -1

endif

       CALL PSGEMR2D(Nmax,Nmax,A,1,1,descA,AO,&
                   1,1,DESCAO,icontxt)
if (myrow==0.and.mycol==0) then
endif
 deallocate(A)

call blacs_gridexit(icontxt)
!call blacs_exit(0)
return
end subroutine lucmp



PostPosted: Wed Apr 16, 2008 9:21 am
by papis
The main problem is tha i use mpi and scalapack.So slinit and gridinit may work wrong together with MPI_init.

Just updating the status :)