!==============================================================================!
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!==============================================================================!
!==============================================================================!
! |
! WSMP_SETUP Oct. 2008 |
! WSMP_CLEANUP Oct. 2008 |
! |
!==============================================================================!
!==============================================================================!
subroutine wsmp_setup (ndof,vo,wsmp,osolve,params,threadinfo,istep,iter)
use definitions
use threads
implicit none
integer ndof
type (void) vo
type (topology),dimension(:),allocatable::tpl
type (wsmp_data) wsmp
type (octreesolve) osolve
type (parameters) params
type (thread) threadinfo
integer istep,iter
integer icol(100000),kk,nn,n0,jproc,idg_locp,nheightp,iatemp
integer ie,k1,k2,j,i2,i1,idof1,idof2,k,err,i3,k3,i,iloc
integer iproc,nproc,ierr
character*72 :: shift
INCLUDE 'mpif.h'
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
shift=' '
!==============================================================================!
!==============================================================================!
! wsmp%n is the size of the system/matrix
! wsmp%nz is the total number of non-zero terms in the global matrix
! wsmp%n_iproc is the number of columns on given thread iproc
! wsmp%n_iproc_st is the first column on given thread iproc, between 1 and n
! wsmp%n_iproc_end is the last column on given thread iproc, between 1 and n
! wsmp%nz_loc is the total number of non-zero terms on given thread iproc
! wsmp%ia is an integer array of size n_iproc+1 (see wsmp manual)
! wsmp%ja is an interer array of size nz_loc
! wsmp%avals is a double precision array of size nz_loc containing the values
! wsmp%b is a double precision array of size n_iproc containing the rhs
!==============================================================================!
!=====[find matrix size]=======================================================!
!==============================================================================!
wsmp%n = vo%nnode * ndof
!==============================================================================!
!=====[allocate memory for topology]===========================================!
!==============================================================================!
allocate (tpl(wsmp%n),stat=err)
tpl%nheight=0
if (ndof==3) then
tpl%nheightmax=17*ndof
elseif(ndof==1) then
tpl%nheightmax=27*ndof
else
stop 'pb in wsmp_Setup'
end if
do i=1,wsmp%n
allocate (tpl(i)%icol(tpl(i)%nheightmax),stat=err)
tpl(i)%icol(:)=0
enddo
call heap (threadinfo,'tpl(:)%icol(:)','wsmp_setup',wsmp%n*tpl%nheightmax,'int',+1)
!==============================================================================!
!=====[find_connectivity_dimension]============================================!
!==============================================================================!
! the list of non-zero terms is computed for every column
! i1 is a column index
! i2 is a row index
! for symmetric systems, we only build and store the lower triangular part
! of the matrix, i.e. i2.ge.i1
do ie=1,osolve%nleaves
if (vo%leaf(ie).eq.0) then
do k1=1,params%mpe
do idof1=1,ndof
i1=(vo%ftr(osolve%icon(k1,ie))-1)*ndof+idof1
do k2=1,params%mpe
do idof2=1,ndof
i2=(vo%ftr(osolve%icon(k2,ie))-1)*ndof+idof2
if (i1.gt.i2.and.ndof.gt.1) goto 2
do j=1,tpl(i1)%nheight
if (i2.eq.tpl(i1)%icol(j)) goto 2
enddo
if (tpl(i1)%nheight.eq.tpl(i1)%nheightmax) then
do k=1,tpl(i1)%nheight
icol(k)=tpl(i1)%icol(k)
enddo
if (tpl(i1)%nheight.gt.0) deallocate (tpl(i1)%icol)
tpl(i1)%nheightmax=tpl(i1)%nheightmax+ndof
allocate (tpl(i1)%icol(tpl(i1)%nheightmax),stat=err)
tpl(i1)%icol(:)=0
if (err.ne.0) call stop_run ('Error alloc tpl(i1)%icol$')
do k=1,tpl(i1)%nheight
tpl(i1)%icol(k)=icol(k)
enddo
endif
tpl(i1)%nheight=tpl(i1)%nheight+1
tpl(i1)%icol(tpl(i1)%nheight)=i2
2 continue
enddo
enddo
enddo
enddo
endif
enddo
!print *,iproc,'nface=',osolve%nface
do ie=1,osolve%nface
if (vo%face(ie).eq.0) then
do k1=1,9
do idof1=1,ndof
i1=(vo%ftr(osolve%iface(k1,ie))-1)*ndof+idof1
do k2=1,9
do idof2=1,ndof
i2=(vo%ftr(osolve%iface(k2,ie))-1)*ndof+idof2
if (i1.gt.i2.and.ndof.gt.1) goto 3 !!!!!!NEW
do j=1,tpl(i1)%nheight
if (i2.eq.tpl(i1)%icol(j)) goto 3
enddo
if (tpl(i1)%nheight.eq.tpl(i1)%nheightmax) then
do k=1,tpl(i1)%nheight
icol(k)=tpl(i1)%icol(k)
enddo
if (tpl(i1)%nheight.gt.0) deallocate (tpl(i1)%icol)
tpl(i1)%nheightmax=tpl(i1)%nheightmax+ndof
allocate (tpl(i1)%icol(tpl(i1)%nheightmax),stat=err)
if (err.ne.0) call stop_run ('Error alloc tpl(i1)%icol$')
do k=1,tpl(i1)%nheight
tpl(i1)%icol(k)=icol(k)
enddo
endif
tpl(i1)%nheight=tpl(i1)%nheight+1
tpl(i1)%icol(tpl(i1)%nheight)=i2
3 continue
enddo
enddo
enddo
enddo
endif
enddo
!print *,iproc,'min/max tpl%nheight',minval(tpl(:)%nheight),maxval(tpl(:)%nheight)
!print *,iproc,'avrg tpl%nheight',sum(tpl(:)%nheight)/dble(wsmp%n)
!==============================================================================!
!=====[computing total number of non-zero terms]===============================!
!==============================================================================!
! the total number of non-zero terms is the sum of the height of each column
wsmp%nz=sum(tpl(1:vo%nnode*ndof)%nheight)
if (iproc.eq.0) write(*,'(a,i9)') shift//'n =',wsmp%n
if (iproc.eq.0) write(*,'(a,i9)') shift//'nz=',wsmp%nz
!==============================================================================!
!=====[sorting icol's]=========================================================!
!==============================================================================!
! wsmp requires the indices of the row contained in each column to be sorted
do i=1,wsmp%n
call iqsort_s (tpl(i)%icol,tpl(i)%nheight)
end do
!==============================================================================!
!=====[allocating memory]======================================================!
!==============================================================================!
allocate(wsmp%irn(wsmp%nz)) ; call heap (threadinfo,'wsmp%irn','wsmp_setup',size(wsmp%irn),'int',+1)
allocate(wsmp%jcn(wsmp%nz)) ; call heap (threadinfo,'wsmp%jcn','wsmp_setup',size(wsmp%jcn),'int',+1)
allocate(wsmp%idg(wsmp%n+1)) ; call heap (threadinfo,'wsmp%idg','wsmp_setup',size(wsmp%idg),'int',+1)
allocate(wsmp%iproc_col(wsmp%n)) ; call heap (threadinfo,'wsmp%iproc_col','wsmp_setup',size(wsmp%iproc_col),'bool',+1)
wsmp%irn=0
wsmp%jcn=0
wsmp%idg=0
!==============================================================================!
!=====[find_connectivity]======================================================!
!==============================================================================!
! the coordinates (irn,jcn) of the non zero terms in the global matrix are
! computed, as well as the indices of the diagonal terms in the list of dof's.
!if (ndof==1) then
!print *,iproc,'irn size',size(wsmp%irn)
!print *,iproc,'jcn size',size(wsmp%jcn)
!print *,iproc,'nz ',wsmp%nz
!end if
kk=0
do i=1,wsmp%n
do j=1,tpl(i)%nheight
kk=kk+1
if (kk.gt.wsmp%nz) call stop_run ('error in find_connectivity$')
wsmp%irn(kk)=tpl(i)%icol(j)
wsmp%jcn(kk)=i
enddo
enddo
!wsmp%idg(1)=1
!do i=2,wsmp%n+1
! wsmp%idg(i)=wsmp%idg(i-1)+tpl(i-1)%nheight
!enddo
!==============================================================================!
!=====[find_processors]========================================================!
!==============================================================================!
nn =(wsmp%n-1)/nproc+1
wsmp%n_iproc_st = nn*iproc+1
wsmp%n_iproc_end = min(wsmp%n_iproc_st+nn-1,wsmp%n)
wsmp%n_iproc = wsmp%n_iproc_end - wsmp%n_iproc_st +1
wsmp%iproc_col=.false.
wsmp%iproc_col(wsmp%n_iproc_st:wsmp%n_iproc_end)=.true.
wsmp%nz_loc=sum(tpl(:)%nheight,wsmp%iproc_col)
!==============================================================================!
!=====[allocating memory]======================================================!
!==============================================================================!
allocate(wsmp%ja(wsmp%nz_loc)) ; call heap (threadinfo,'wsmp%ja','wsmp_setup',size(wsmp%ja),'int',+1)
allocate(wsmp%ia(wsmp%n_iproc+1)) ; call heap (threadinfo,'wsmp%ia','wsmp_setup',size(wsmp%ia),'int',+1)
!allocate(wsmp%irn_loc(wsmp%nz_loc))
!allocate(wsmp%jcn_loc(wsmp%nz_loc))
!allocate(wsmp%idg_loc(wsmp%n+1))
!wsmp%irn_loc=0
!wsmp%jcn_loc=0
!wsmp%idg_loc=0
!==============================================================================!
!=====[find_connectivity_local]================================================!
!==============================================================================!
!idg_locp=0
!do i=1,wsmp%n
! if (wsmp%iproc_col(i)) then
! if (idg_locp.eq.0) then
! wsmp%idg_loc(i)=1
! else
! wsmp%idg_loc(i)=idg_locp+nheightp
! endif
! idg_locp=wsmp%idg_loc(i)
! nheightp=tpl(i)%nheight
! endif
!enddo
!wsmp%idg_loc(wsmp%n+1)=idg_locp+nheightp
iloc=1
do i=1,wsmp%n
if (wsmp%iproc_col(i)) then
do k=1,tpl(i)%nheight
! wsmp%irn_loc(iloc+k-1)=wsmp%irn(wsmp%idg(i)+k-1)
! wsmp%irn_loc(iloc+k-1)=tpl(i)%icol(k)
wsmp%ja(iloc+k-1)=tpl(i)%icol(k)
! wsmp%jcn_loc(iloc+k-1)=wsmp%jcn(wsmp%idg(i)+k-1)
! wsmp%jcn_loc(iloc+k-1)=i
enddo
iloc=iloc+tpl(i)%nheight
endif
enddo
!==============================================================================!
!=====[outputs visual reprenstation of matrix]=================================!
!==============================================================================!
call visualise_matrix (params%visualise_matrix,wsmp%nz,wsmp%irn,wsmp%jcn,wsmp%n,istep,iter,ndof)
!==============================================================================!
!=====[allocate memory for matrix values]======================================!
!==============================================================================!
allocate(wsmp%avals(wsmp%nz_loc)) ; call heap (threadinfo,'wsmp%avals','wsmp_setup',size(wsmp%avals),'dp',+1)
wsmp%avals=0.d0
wsmp%nrhs = 1
wsmp%ldb = wsmp%n_iproc
wsmp%naux = 0
wsmp%mrp=0
allocate(wsmp%b (wsmp%ldb,wsmp%nrhs)) ; call heap (threadinfo,'wsmp%b','wsmp_setup',size(wsmp%b),'dp',+1)
wsmp%b=0.d0
if (ndof==1) then
allocate(wsmp%rmisc(wsmp%ldb,wsmp%nrhs)) ; call heap (threadinfo,'wsmp%rmisc','wsmp_setup',size(wsmp%rmisc),'dp',+1)
wsmp%rmisc=0.d0
end if
if (ndof==3) then
allocate(wsmp%perm (wsmp%n)) ; call heap (threadinfo,'wsmp%perm','wsmp_setup',size(wsmp%perm),'int',+1)
allocate(wsmp%invp (wsmp%n)) ; call heap (threadinfo,'wsmp%invp','wsmp_setup',size(wsmp%invp),'int',+1)
wsmp%perm=0
wsmp%invp=0
end if
!==============================================================================!
!=====[build ia,ja]============================================================!
!==============================================================================!
iatemp=0
k=0
do i=wsmp%n_iproc_st,wsmp%n_iproc_end
k=k+1
if (iatemp==0) then
wsmp%ia(k)=1
else
wsmp%ia(k)=iatemp+nheightp
end if
iatemp=wsmp%ia(k)
nheightp=tpl(i)%nheight
end do
k=k+1
wsmp%ia(k)=iatemp+nheightp
!if (ndof==1) then
!write(*,*) iproc,'n , nz',wsmp%n,wsmp%nz
!write(*,*) iproc,'col st/end ',wsmp%n_iproc_st,wsmp%n_iproc_end
!write(*,*) iproc,'n_iproc',wsmp%n_iproc
!write(*,*) iproc,'nz_loc',wsmp%nz_loc
!write(*,*) iproc,'min/max ia',minval(wsmp%ia),maxval(wsmp%ia)
!write(*,*) iproc,'min/max ja',minval(wsmp%ja),maxval(wsmp%ja)
!end if
!==============================================================================!
!=====[deallocate memory]======================================================!
!==============================================================================!
call heap (threadinfo,'wsmp%irn','wsmp_setup',size(wsmp%irn),'int',-1) ; deallocate (wsmp%irn)
call heap (threadinfo,'wsmp%jcn','wsmp_setup',size(wsmp%jcn),'int',-1) ; deallocate (wsmp%jcn)
call heap (threadinfo,'wsmp%idg','wsmp_setup',size(wsmp%idg),'int',-1) ; deallocate (wsmp%idg)
call heap (threadinfo,'tpl(:)%icol(:)','wsmp_setup',wsmp%n*tpl%nheightmax,'int',-1)
do i=1,wsmp%n
deallocate (tpl(i)%icol)
enddo
deallocate (tpl)
end subroutine wsmp_setup
!==============================================================================!
!==============================================================================!
! |
! WSMP_CLEANUP Oct. 2008 |
! |
!==============================================================================!
!==============================================================================!
subroutine wsmp_cleanup (wsmp,threadinfo)
use definitions
use threads
implicit none
type(wsmp_data) wsmp
type (thread) threadinfo
call heap (threadinfo,'wsmp%ia','wsmp_cleanup',size(wsmp%ia),'int',-1) ; deallocate(wsmp%ia)
call heap (threadinfo,'wsmp%ja','wsmp_cleanup',size(wsmp%ja),'int',-1) ; deallocate(wsmp%ja)
call heap (threadinfo,'wsmp%iproc_col','wsmp_cleanup',size(wsmp%iproc_col),'bool',-1) ; deallocate(wsmp%iproc_col)
call heap (threadinfo,'wsmp%avals','wsmp_cleanup',size(wsmp%avals),'dp',-1) ; deallocate(wsmp%avals)
call heap (threadinfo,'wsmp%b','wsmp_cleanup',size(wsmp%b),'dp',-1) ; deallocate(wsmp%b)
if (allocated(wsmp%perm)) then
call heap (threadinfo,'wsmp%perm','wsmp_cleanup',size(wsmp%perm),'int',-1)
deallocate(wsmp%perm)
end if
if (allocated(wsmp%invp)) then
call heap (threadinfo,'wsmp%invp','wsmp_cleanup',size(wsmp%invp),'int',-1)
deallocate(wsmp%invp)
end if
if (allocated(wsmp%rmisc)) then
call heap (threadinfo,'wsmp%rmisc','wsmp_cleanup',size(wsmp%rmisc),'dp',-1)
deallocate(wsmp%rmisc)
end if
!call pwsmp_clear()
end subroutine wsmp_cleanup
!==============================================================================!
!==============================================================================!