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
