fdurastante
/
amg4psblas
mirror of https://github.com/sfilippone/amg4psblas.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

mld2p4-2:

Browse Source 
mlprec/Makefile
 mlprec/mld_c_dec_map_bld.F90
 mlprec/mld_c_inner_mod.f90
 mlprec/mld_caggrmap_bld.f90
 mlprec/mld_d_dec_map_bld.F90
 mlprec/mld_d_ilu_solver.f90
 mlprec/mld_d_inner_mod.f90
 mlprec/mld_daggrmap_bld.f90
 mlprec/mld_s_dec_map_bld.F90
 mlprec/mld_s_inner_mod.f90
 mlprec/mld_saggrmap_bld.f90
 mlprec/mld_z_dec_map_bld.F90
 mlprec/mld_z_inner_mod.f90
 mlprec/mld_zaggrmap_bld.f90

Separate the decmap_bld from its wrapper, make room for alternatives.
stopcriterion
Salvatore Filippone 14 years ago
parent 1f91926ff4
commit 7fa413283e
14 changed files with 1243 additions and 1210 deletions
Show Stats Download Patch File Download Diff File

8
mlprec/Makefile
Unescape Escape View File

@ -36,22 +36,22 @@ MPCOBJS=mld_sslud_interface.o mld_dslud_interface.o mld_cslud_interface.o mld_zs
SINNEROBJS= mld_scoarse_bld.o mld_smlprec_bld.o \
mld_silu0_fact.o mld_siluk_fact.o mld_silut_fact.o mld_saggrmap_bld.o \
mld_smlprec_aply.o mld_saggrmat_asb.o \
mld_s_dec_map_bld.o mld_smlprec_aply.o mld_saggrmat_asb.o \
$(SMPFOBJS)
DINNEROBJS= mld_dcoarse_bld.o mld_dmlprec_bld.o \
mld_dilu0_fact.o mld_diluk_fact.o mld_dilut_fact.o mld_daggrmap_bld.o \
mld_dmlprec_aply.o mld_daggrmat_asb.o \
mld_d_dec_map_bld.o mld_dmlprec_aply.o mld_daggrmat_asb.o \
$(DMPFOBJS)
CINNEROBJS= mld_ccoarse_bld.o mld_cmlprec_bld.o \
mld_cilu0_fact.o mld_ciluk_fact.o mld_cilut_fact.o mld_caggrmap_bld.o \
mld_cmlprec_aply.o mld_caggrmat_asb.o \
mld_c_dec_map_bld.o mld_cmlprec_aply.o mld_caggrmat_asb.o \
$(CMPFOBJS)
ZINNEROBJS= mld_zcoarse_bld.o mld_zmlprec_bld.o \
mld_zilu0_fact.o mld_ziluk_fact.o mld_zilut_fact.o mld_zaggrmap_bld.o \
mld_zmlprec_aply.o mld_zaggrmat_asb.o \
mld_z_dec_map_bld.o mld_zmlprec_aply.o mld_zaggrmat_asb.o \
$(ZMPFOBJS)
INNEROBJS= $(SINNEROBJS) $(DINNEROBJS) $(CINNEROBJS) $(ZINNEROBJS)

300
mlprec/mld_c_dec_map_bld.F90
Unescape Escape View File

@ -0,0 +1,300 @@
subroutine mld_c_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod
use mld_c_inner_mod, mld_protect_name => mld_c_dec_map_bld
implicit none
! Arguments
type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
real(psb_spk_), intent(in) :: theta
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
! Local variables
integer, allocatable :: ils(:), neigh(:), irow(:), icol(:)
complex(psb_spk_), allocatable :: val(:), diag(:)
integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz
real(psb_spk_) :: cpling, tcl
logical :: recovery
integer :: debug_level, debug_unit
integer :: ictxt,np,me,err_act
integer :: nrow, ncol, n_ne
character(len=20) :: name, ch_err
if(psb_get_errstatus() /= 0) return
info=psb_success_
name = 'mld_dec_map_bld'
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!
ictxt=psb_cd_get_context(desc_a)
call psb_info(ictxt,me,np)
nrow = psb_cd_get_local_rows(desc_a)
ncol = psb_cd_get_local_cols(desc_a)
nr = a%get_nrows()
allocate(ilaggr(nr),neigh(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nr,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
allocate(diag(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/nr,0,0,0,0/),&
& a_err='complex(psb_spk_)')
goto 9999
end if
call a%get_diag(diag,info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_sp_getdiag')
goto 9999
end if
do i=1, nr
ilaggr(i) = -(nr+1)
end do
! Note: -(nr+1) Untouched as yet
! -i 1<=i<=nr Adjacent to aggregate i
! i 1<=i<=nr Belonging to aggregate i
!
! Phase one: group nodes together.
! Very simple minded strategy.
!
naggr = 0
nlp = 0
do
icnt = 0
do i=1, nr
if (ilaggr(i) == -(nr+1)) then
!
! 1. Untouched nodes are marked >0 together
! with their neighbours
!
icnt = icnt + 1
naggr = naggr + 1
ilaggr(i) = naggr
call a%csget(i,i,nz,irow,icol,val,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='csget')
goto 9999
end if
do k=1, nz
j = icol(k)
if ((1<=j).and.(j<=nr).and.(i /= j)) then
if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
ilaggr(j) = naggr
else
ilaggr(j) = -naggr
endif
end if
enddo
!
! 2. Untouched neighbours of these nodes are marked <0.
!
call a%get_neigh(i,neigh,n_ne,info,lev=2)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_neigh')
goto 9999
end if
do n = 1, n_ne
m = neigh(n)
if ((1<=m).and.(m<=nr)) then
if (ilaggr(m) == -(nr+1)) ilaggr(m) = -naggr
endif
enddo
endif
enddo
nlp = nlp + 1
if (icnt == 0) exit
enddo
if (debug_level >= psb_debug_outer_) then
write(debug_unit,*) me,' ',trim(name),&
& ' Check 1:',count(ilaggr == -(nr+1))
end if
!
! Phase two: sweep over leftovers.
!
allocate(ils(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
do i=1, size(ils)
ils(i) = 0
end do
do i=1, nr
n = ilaggr(i)
if (n>0) then
if (n>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 1 ?')
goto 9999
else
ils(n) = ils(n) + 1
end if
end if
end do
if (debug_level >= psb_debug_outer_) then
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 1: number of aggregates ',naggr
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 1: nodes aggregated ',sum(ils)
end if
recovery=.false.
do i=1, nr
if (ilaggr(i) < 0) then
!
! Now some silly rule to break ties:
! Group with adjacent aggregate.
!
isz = nr+1
ia = -1
cpling = szero
call a%csget(i,i,nz,irow,icol,val,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_sp_getrow')
goto 9999
end if
do j=1, nz
k = icol(j)
if ((1<=k).and.(k<=nr).and.(k /= i)) then
tcl = abs(val(j)) / sqrt(abs(diag(i)*diag(k)))
if (abs(val(j)) > theta*sqrt(abs(diag(i)*diag(k)))) then
!!$ if (tcl > theta) then
n = ilaggr(k)
if (n>0) then
if (n>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?')
goto 9999
end if
if ((abs(val(j))>cpling) .or. &
& ((abs(val(j)) == cpling).and. (ils(n) < isz))) then
!!$ if ((tcl > cpling) .or. ((tcl == cpling).and. (ils(n) < isz))) then
ia = n
isz = ils(n)
cpling = abs(val(j))
!!$ cpling = tcl
endif
endif
endif
end if
enddo
if (ia == -1) then
! At this point, the easiest thing is to start a new aggregate
naggr = naggr + 1
ilaggr(i) = naggr
ils(ilaggr(i)) = 1
else
ilaggr(i) = ia
if (ia>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr ? ')
goto 9999
end if
ils(ia) = ils(ia) + 1
endif
end if
end do
if (debug_level >= psb_debug_outer_) then
if (recovery) then
write(debug_unit,*) me,' ',trim(name),&
& 'Had to recover from strange situation in loc_aggregate.'
write(debug_unit,*) me,' ',trim(name),&
& 'Perhaps an unsymmetric pattern?'
endif
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 2: number of aggregates ',naggr,sum(ils)
do i=1, naggr
write(debug_unit,*) me,' ',trim(name),&
& 'Size of aggregate ',i,' :',count(ilaggr == i), ils(i)
enddo
write(debug_unit,*) me,' ',trim(name),&
& maxval(ils(1:naggr))
write(debug_unit,*) me,' ',trim(name),&
& 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops'
end if
if (count(ilaggr<0) >0) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Fatal error: some leftovers')
goto 9999
endif
deallocate(ils,neigh,stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
if (naggr > ncol) then
write(0,*) name,'Error : naggr > ncol',naggr,ncol
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Fatal error: naggr>ncol')
goto 9999
end if
call psb_realloc(ncol,ilaggr,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
allocate(nlaggr(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/np,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
nlaggr(:) = 0
nlaggr(me+1) = naggr
call psb_sum(ictxt,nlaggr(1:np))
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine mld_c_dec_map_bld

12
mlprec/mld_c_inner_mod.f90
Unescape Escape View File

@ -102,6 +102,18 @@ module mld_c_inner_mod
end subroutine mld_caggrmap_bld
end interface mld_aggrmap_bld
interface mld_dec_map_bld
subroutine mld_c_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_
type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
real(psb_spk_), intent(in) :: theta
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
end subroutine mld_c_dec_map_bld
end interface mld_dec_map_bld
interface mld_aggrmat_asb
subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_

301
mlprec/mld_caggrmap_bld.f90
Unescape Escape View File

@ -160,305 +160,4 @@ subroutine mld_caggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info)
end if
return
contains
subroutine mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod
use mld_c_inner_mod !, mld_protect_name => mld_daggrmap_bld
implicit none
! Arguments
type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
real(psb_spk_), intent(in) :: theta
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
! Local variables
integer, allocatable :: ils(:), neigh(:), irow(:), icol(:)
complex(psb_spk_), allocatable :: val(:), diag(:)
integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz
real(psb_spk_) :: cpling, tcl
logical :: recovery
integer :: debug_level, debug_unit
integer :: ictxt,np,me,err_act
integer :: nrow, ncol, n_ne
character(len=20) :: name, ch_err
if(psb_get_errstatus() /= 0) return
info=psb_success_
name = 'mld_dec_map_bld'
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!
ictxt=psb_cd_get_context(desc_a)
call psb_info(ictxt,me,np)
nrow = psb_cd_get_local_rows(desc_a)
ncol = psb_cd_get_local_cols(desc_a)
nr = a%get_nrows()
allocate(ilaggr(nr),neigh(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nr,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
allocate(diag(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/nr,0,0,0,0/),&
& a_err='complex(psb_spk_)')
goto 9999
end if
call a%get_diag(diag,info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_sp_getdiag')
goto 9999
end if
do i=1, nr
ilaggr(i) = -(nr+1)
end do
! Note: -(nr+1) Untouched as yet
! -i 1<=i<=nr Adjacent to aggregate i
! i 1<=i<=nr Belonging to aggregate i
!
! Phase one: group nodes together.
! Very simple minded strategy.
!
naggr = 0
nlp = 0
do
icnt = 0
do i=1, nr
if (ilaggr(i) == -(nr+1)) then
!
! 1. Untouched nodes are marked >0 together
! with their neighbours
!
icnt = icnt + 1
naggr = naggr + 1
ilaggr(i) = naggr
call a%csget(i,i,nz,irow,icol,val,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='csget')
goto 9999
end if
do k=1, nz
j = icol(k)
if ((1<=j).and.(j<=nr).and.(i /= j)) then
if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
ilaggr(j) = naggr
else
ilaggr(j) = -naggr
endif
end if
enddo
!
! 2. Untouched neighbours of these nodes are marked <0.
!
call a%get_neigh(i,neigh,n_ne,info,lev=2)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_neigh')
goto 9999
end if
do n = 1, n_ne
m = neigh(n)
if ((1<=m).and.(m<=nr)) then
if (ilaggr(m) == -(nr+1)) ilaggr(m) = -naggr
endif
enddo
endif
enddo
nlp = nlp + 1
if (icnt == 0) exit
enddo
if (debug_level >= psb_debug_outer_) then
write(debug_unit,*) me,' ',trim(name),&
& ' Check 1:',count(ilaggr == -(nr+1))
end if
!
! Phase two: sweep over leftovers.
!
allocate(ils(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
do i=1, size(ils)
ils(i) = 0
end do
do i=1, nr
n = ilaggr(i)
if (n>0) then
if (n>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 1 ?')
goto 9999
else
ils(n) = ils(n) + 1
end if
end if
end do
if (debug_level >= psb_debug_outer_) then
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 1: number of aggregates ',naggr
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 1: nodes aggregated ',sum(ils)
end if
recovery=.false.
do i=1, nr
if (ilaggr(i) < 0) then
!
! Now some silly rule to break ties:
! Group with adjacent aggregate.
!
isz = nr+1
ia = -1
cpling = szero
call a%csget(i,i,nz,irow,icol,val,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_sp_getrow')
goto 9999
end if
do j=1, nz
k = icol(j)
if ((1<=k).and.(k<=nr).and.(k /= i)) then
tcl = abs(val(j)) / sqrt(abs(diag(i)*diag(k)))
if (abs(val(j)) > theta*sqrt(abs(diag(i)*diag(k)))) then
!!$ if (tcl > theta) then
n = ilaggr(k)
if (n>0) then
if (n>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?')
goto 9999
end if
if ((abs(val(j))>cpling) .or. &
& ((abs(val(j)) == cpling).and. (ils(n) < isz))) then
!!$ if ((tcl > cpling) .or. ((tcl == cpling).and. (ils(n) < isz))) then
ia = n
isz = ils(n)
cpling = abs(val(j))
!!$ cpling = tcl
endif
endif
endif
end if
enddo
if (ia == -1) then
! At this point, the easiest thing is to start a new aggregate
naggr = naggr + 1
ilaggr(i) = naggr
ils(ilaggr(i)) = 1
else
ilaggr(i) = ia
if (ia>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr ? ')
goto 9999
end if
ils(ia) = ils(ia) + 1
endif
end if
end do
if (debug_level >= psb_debug_outer_) then
if (recovery) then
write(debug_unit,*) me,' ',trim(name),&
& 'Had to recover from strange situation in loc_aggregate.'
write(debug_unit,*) me,' ',trim(name),&
& 'Perhaps an unsymmetric pattern?'
endif
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 2: number of aggregates ',naggr,sum(ils)
do i=1, naggr
write(debug_unit,*) me,' ',trim(name),&
& 'Size of aggregate ',i,' :',count(ilaggr == i), ils(i)
enddo
write(debug_unit,*) me,' ',trim(name),&
& maxval(ils(1:naggr))
write(debug_unit,*) me,' ',trim(name),&
& 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops'
end if
if (count(ilaggr<0) >0) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Fatal error: some leftovers')
goto 9999
endif
deallocate(ils,neigh,stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
if (naggr > ncol) then
write(0,*) name,'Error : naggr > ncol',naggr,ncol
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Fatal error: naggr>ncol')
goto 9999
end if
call psb_realloc(ncol,ilaggr,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
allocate(nlaggr(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/np,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
nlaggr(:) = 0
nlaggr(me+1) = naggr
call psb_sum(ictxt,nlaggr(1:np))
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine mld_dec_map_bld
end subroutine mld_caggrmap_bld

300
mlprec/mld_d_dec_map_bld.F90
Unescape Escape View File

@ -0,0 +1,300 @@
subroutine mld_d_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod
use mld_d_inner_mod, mld_protect_name => mld_d_dec_map_bld
implicit none
! Arguments
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
real(psb_dpk_), intent(in) :: theta
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
! Local variables
integer, allocatable :: ils(:), neigh(:), irow(:), icol(:)
real(psb_dpk_), allocatable :: val(:), diag(:)
integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz
real(psb_dpk_) :: cpling, tcl
logical :: recovery
integer :: debug_level, debug_unit
integer :: ictxt,np,me,err_act
integer :: nrow, ncol, n_ne
character(len=20) :: name, ch_err
if (psb_get_errstatus() /= 0) return
info=psb_success_
name = 'mld_dec_map_bld'
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!
ictxt=psb_cd_get_context(desc_a)
call psb_info(ictxt,me,np)
nrow = psb_cd_get_local_rows(desc_a)
ncol = psb_cd_get_local_cols(desc_a)
nr = a%get_nrows()
allocate(ilaggr(nr),neigh(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nr,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
allocate(diag(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/nr,0,0,0,0/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
call a%get_diag(diag,info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_sp_getdiag')
goto 9999
end if
do i=1, nr
ilaggr(i) = -(nr+1)
end do
! Note: -(nr+1) Untouched as yet
! -i 1<=i<=nr Adjacent to aggregate i
! i 1<=i<=nr Belonging to aggregate i
!
! Phase one: group nodes together.
! Very simple minded strategy.
!
naggr = 0
nlp = 0
do
icnt = 0
do i=1, nr
if (ilaggr(i) == -(nr+1)) then
!
! 1. Untouched nodes are marked >0 together
! with their neighbours
!
icnt = icnt + 1
naggr = naggr + 1
ilaggr(i) = naggr
call a%csget(i,i,nz,irow,icol,val,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='csget')
goto 9999
end if
do k=1, nz
j = icol(k)
if ((1<=j).and.(j<=nr).and.(i /= j)) then
if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
ilaggr(j) = naggr
else
ilaggr(j) = -naggr
endif
end if
enddo
!
! 2. Untouched neighbours of these nodes are marked <0.
!
call a%get_neigh(i,neigh,n_ne,info,lev=2)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_neigh')
goto 9999
end if
do n = 1, n_ne
m = neigh(n)
if ((1<=m).and.(m<=nr)) then
if (ilaggr(m) == -(nr+1)) ilaggr(m) = -naggr
endif
enddo
endif
enddo
nlp = nlp + 1
if (icnt == 0) exit
enddo
if (debug_level >= psb_debug_outer_) then
write(debug_unit,*) me,' ',trim(name),&
& ' Check 1:',count(ilaggr == -(nr+1))
end if
!
! Phase two: sweep over leftovers.
!
allocate(ils(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
do i=1, size(ils)
ils(i) = 0
end do
do i=1, nr
n = ilaggr(i)
if (n>0) then
if (n>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 1 ?')
goto 9999
else
ils(n) = ils(n) + 1
end if
end if
end do
if (debug_level >= psb_debug_outer_) then
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 1: number of aggregates ',naggr
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 1: nodes aggregated ',sum(ils)
end if
recovery=.false.
do i=1, nr
if (ilaggr(i) < 0) then
!
! Now some silly rule to break ties:
! Group with adjacent aggregate.
!
isz = nr+1
ia = -1
cpling = dzero
call a%csget(i,i,nz,irow,icol,val,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_sp_getrow')
goto 9999
end if
do j=1, nz
k = icol(j)
if ((1<=k).and.(k<=nr).and.(k /= i)) then
tcl = abs(val(j)) / sqrt(abs(diag(i)*diag(k)))
if (abs(val(j)) > theta*sqrt(abs(diag(i)*diag(k)))) then
!!$ if (tcl > theta) then
n = ilaggr(k)
if (n>0) then
if (n>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?')
goto 9999
end if
if ((abs(val(j))>cpling) .or. &
& ((abs(val(j)) == cpling).and. (ils(n) < isz))) then
!!$ if ((tcl > cpling) .or. ((tcl == cpling).and. (ils(n) < isz))) then
ia = n
isz = ils(n)
cpling = abs(val(j))
!!$ cpling = tcl
endif
endif
endif
end if
enddo
if (ia == -1) then
! At this point, the easiest thing is to start a new aggregate
naggr = naggr + 1
ilaggr(i) = naggr
ils(ilaggr(i)) = 1
else
ilaggr(i) = ia
if (ia>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr ? ')
goto 9999
end if
ils(ia) = ils(ia) + 1
endif
end if
end do
if (debug_level >= psb_debug_outer_) then
if (recovery) then
write(debug_unit,*) me,' ',trim(name),&
& 'Had to recover from strange situation in loc_aggregate.'
write(debug_unit,*) me,' ',trim(name),&
& 'Perhaps an unsymmetric pattern?'
endif
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 2: number of aggregates ',naggr,sum(ils)
do i=1, naggr
write(debug_unit,*) me,' ',trim(name),&
& 'Size of aggregate ',i,' :',count(ilaggr == i), ils(i)
enddo
write(debug_unit,*) me,' ',trim(name),&
& maxval(ils(1:naggr))
write(debug_unit,*) me,' ',trim(name),&
& 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops'
end if
if (count(ilaggr<0) >0) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Fatal error: some leftovers')
goto 9999
endif
deallocate(ils,neigh,stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
if (naggr > ncol) then
write(0,*) name,'Error : naggr > ncol',naggr,ncol
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Fatal error: naggr>ncol')
goto 9999
end if
call psb_realloc(ncol,ilaggr,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
allocate(nlaggr(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/np,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
nlaggr(:) = 0
nlaggr(me+1) = naggr
call psb_sum(ictxt,nlaggr(1:np))
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine mld_d_dec_map_bld

2
mlprec/mld_d_ilu_solver.f90
Unescape Escape View File

@ -407,8 +407,6 @@ contains
if (present(mold)) then
call sv%l%cscnv(info,mold=mold)
call sv%u%cscnv(info,mold=mold)
Write(0,*) 'Converted L into ',sv%l%get_fmt(),&
&' and U into ',sv%u%get_fmt()
end if
if (debug_level >= psb_debug_outer_) &

12
mlprec/mld_d_inner_mod.f90
Unescape Escape View File

@ -104,6 +104,18 @@ module mld_d_inner_mod
end subroutine mld_daggrmap_bld
end interface mld_aggrmap_bld
interface mld_dec_map_bld
subroutine mld_d_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
real(psb_dpk_), intent(in) :: theta
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
end subroutine mld_d_dec_map_bld
end interface mld_dec_map_bld
interface mld_aggrmat_asb
subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_

301
mlprec/mld_daggrmap_bld.f90
Unescape Escape View File

@ -160,305 +160,4 @@ subroutine mld_daggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info)
end if
return
contains
subroutine mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod
use mld_d_inner_mod !, mld_protect_name => mld_daggrmap_bld
implicit none
! Arguments
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
real(psb_dpk_), intent(in) :: theta
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
! Local variables
integer, allocatable :: ils(:), neigh(:), irow(:), icol(:)
real(psb_dpk_), allocatable :: val(:), diag(:)
integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz
real(psb_dpk_) :: cpling, tcl
logical :: recovery
integer :: debug_level, debug_unit
integer :: ictxt,np,me,err_act
integer :: nrow, ncol, n_ne
character(len=20) :: name, ch_err
if (psb_get_errstatus() /= 0) return
info=psb_success_
name = 'mld_dec_map_bld'
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!
ictxt=psb_cd_get_context(desc_a)
call psb_info(ictxt,me,np)
nrow = psb_cd_get_local_rows(desc_a)
ncol = psb_cd_get_local_cols(desc_a)
nr = a%get_nrows()
allocate(ilaggr(nr),neigh(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nr,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
allocate(diag(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/nr,0,0,0,0/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
call a%get_diag(diag,info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_sp_getdiag')
goto 9999
end if
do i=1, nr
ilaggr(i) = -(nr+1)
end do
! Note: -(nr+1) Untouched as yet
! -i 1<=i<=nr Adjacent to aggregate i
! i 1<=i<=nr Belonging to aggregate i
!
! Phase one: group nodes together.
! Very simple minded strategy.
!
naggr = 0
nlp = 0
do
icnt = 0
do i=1, nr
if (ilaggr(i) == -(nr+1)) then
!
! 1. Untouched nodes are marked >0 together
! with their neighbours
!
icnt = icnt + 1
naggr = naggr + 1
ilaggr(i) = naggr
call a%csget(i,i,nz,irow,icol,val,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='csget')
goto 9999
end if
do k=1, nz
j = icol(k)
if ((1<=j).and.(j<=nr).and.(i /= j)) then
if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
ilaggr(j) = naggr
else
ilaggr(j) = -naggr
endif
end if
enddo
!
! 2. Untouched neighbours of these nodes are marked <0.
!
call a%get_neigh(i,neigh,n_ne,info,lev=2)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_neigh')
goto 9999
end if
do n = 1, n_ne
m = neigh(n)
if ((1<=m).and.(m<=nr)) then
if (ilaggr(m) == -(nr+1)) ilaggr(m) = -naggr
endif
enddo
endif
enddo
nlp = nlp + 1
if (icnt == 0) exit
enddo
if (debug_level >= psb_debug_outer_) then
write(debug_unit,*) me,' ',trim(name),&
& ' Check 1:',count(ilaggr == -(nr+1))
end if
!
! Phase two: sweep over leftovers.
!
allocate(ils(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
do i=1, size(ils)
ils(i) = 0
end do
do i=1, nr
n = ilaggr(i)
if (n>0) then
if (n>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 1 ?')
goto 9999
else
ils(n) = ils(n) + 1
end if
end if
end do
if (debug_level >= psb_debug_outer_) then
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 1: number of aggregates ',naggr
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 1: nodes aggregated ',sum(ils)
end if
recovery=.false.
do i=1, nr
if (ilaggr(i) < 0) then
!
! Now some silly rule to break ties:
! Group with adjacent aggregate.
!
isz = nr+1
ia = -1
cpling = dzero
call a%csget(i,i,nz,irow,icol,val,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_sp_getrow')
goto 9999
end if
do j=1, nz
k = icol(j)
if ((1<=k).and.(k<=nr).and.(k /= i)) then
tcl = abs(val(j)) / sqrt(abs(diag(i)*diag(k)))
if (abs(val(j)) > theta*sqrt(abs(diag(i)*diag(k)))) then
!!$ if (tcl > theta) then
n = ilaggr(k)
if (n>0) then
if (n>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?')
goto 9999
end if
if ((abs(val(j))>cpling) .or. &
& ((abs(val(j)) == cpling).and. (ils(n) < isz))) then
!!$ if ((tcl > cpling) .or. ((tcl == cpling).and. (ils(n) < isz))) then
ia = n
isz = ils(n)
cpling = abs(val(j))
!!$ cpling = tcl
endif
endif
endif
end if
enddo
if (ia == -1) then
! At this point, the easiest thing is to start a new aggregate
naggr = naggr + 1
ilaggr(i) = naggr
ils(ilaggr(i)) = 1
else
ilaggr(i) = ia
if (ia>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr ? ')
goto 9999
end if
ils(ia) = ils(ia) + 1
endif
end if
end do
if (debug_level >= psb_debug_outer_) then
if (recovery) then
write(debug_unit,*) me,' ',trim(name),&
& 'Had to recover from strange situation in loc_aggregate.'
write(debug_unit,*) me,' ',trim(name),&
& 'Perhaps an unsymmetric pattern?'
endif
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 2: number of aggregates ',naggr,sum(ils)
do i=1, naggr
write(debug_unit,*) me,' ',trim(name),&
& 'Size of aggregate ',i,' :',count(ilaggr == i), ils(i)
enddo
write(debug_unit,*) me,' ',trim(name),&
& maxval(ils(1:naggr))
write(debug_unit,*) me,' ',trim(name),&
& 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops'
end if
if (count(ilaggr<0) >0) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Fatal error: some leftovers')
goto 9999
endif
deallocate(ils,neigh,stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
if (naggr > ncol) then
write(0,*) name,'Error : naggr > ncol',naggr,ncol
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Fatal error: naggr>ncol')
goto 9999
end if
call psb_realloc(ncol,ilaggr,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
allocate(nlaggr(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/np,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
nlaggr(:) = 0
nlaggr(me+1) = naggr
call psb_sum(ictxt,nlaggr(1:np))
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine mld_dec_map_bld
end subroutine mld_daggrmap_bld

296
mlprec/mld_s_dec_map_bld.F90
Unescape Escape View File

@ -0,0 +1,296 @@
subroutine mld_s_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod
use mld_s_inner_mod, mld_protect_name => mld_s_dec_map_bld
implicit none
! Arguments
type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
real(psb_spk_), intent(in) :: theta
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
! Local variables
integer, allocatable :: ils(:), neigh(:), irow(:), icol(:)
real(psb_spk_), allocatable :: val(:), diag(:)
integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz
real(psb_spk_) :: cpling, tcl
logical :: recovery
integer :: debug_level, debug_unit
integer :: ictxt,np,me,err_act
integer :: nrow, ncol, n_ne
character(len=20) :: name, ch_err
if(psb_get_errstatus() /= 0) return
info=psb_success_
name = 'mld_dec_map_bld'
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!
ictxt=psb_cd_get_context(desc_a)
call psb_info(ictxt,me,np)
nrow = psb_cd_get_local_rows(desc_a)
ncol = psb_cd_get_local_cols(desc_a)
nr = a%get_nrows()
allocate(ilaggr(nr),neigh(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nr,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
allocate(diag(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/nr,0,0,0,0/),&
& a_err='real(psb_spk_)')
goto 9999
end if
call a%get_diag(diag,info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_sp_getdiag')
goto 9999
end if
do i=1, nr
ilaggr(i) = -(nr+1)
end do
! Note: -(nr+1) Untouched as yet
! -i 1<=i<=nr Adjacent to aggregate i
! i 1<=i<=nr Belonging to aggregate i
!
! Phase one: group nodes together.
! Very simple minded strategy.
!
naggr = 0
nlp = 0
do
icnt = 0
do i=1, nr
if (ilaggr(i) == -(nr+1)) then
!
! 1. Untouched nodes are marked >0 together
! with their neighbours
!
icnt = icnt + 1
naggr = naggr + 1
ilaggr(i) = naggr
call a%csget(i,i,nz,irow,icol,val,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='csget')
goto 9999
end if
do k=1, nz
j = icol(k)
if ((1<=j).and.(j<=nr).and.(i /= j)) then
if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
ilaggr(j) = naggr
else
ilaggr(j) = -naggr
endif
end if
enddo
!
! 2. Untouched neighbours of these nodes are marked <0.
!
call a%get_neigh(i,neigh,n_ne,info,lev=2)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_neigh')
goto 9999
end if
do n = 1, n_ne
m = neigh(n)
if ((1<=m).and.(m<=nr)) then
if (ilaggr(m) == -(nr+1)) ilaggr(m) = -naggr
endif
enddo
endif
enddo
nlp = nlp + 1
if (icnt == 0) exit
enddo
if (debug_level >= psb_debug_outer_) then
write(debug_unit,*) me,' ',trim(name),&
& ' Check 1:',count(ilaggr == -(nr+1))
end if
!
! Phase two: sweep over leftovers.
!
allocate(ils(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
do i=1, size(ils)
ils(i) = 0
end do
do i=1, nr
n = ilaggr(i)
if (n>0) then
if (n>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 1 ?')
goto 9999
else
ils(n) = ils(n) + 1
end if
end if
end do
if (debug_level >= psb_debug_outer_) then
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 1: number of aggregates ',naggr
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 1: nodes aggregated ',sum(ils)
end if
recovery=.false.
do i=1, nr
if (ilaggr(i) < 0) then
!
! Now some silly rule to break ties:
! Group with adjacent aggregate.
!
isz = nr+1
ia = -1
cpling = szero
call a%csget(i,i,nz,irow,icol,val,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_sp_getrow')
goto 9999
end if
do j=1, nz
k = icol(j)
if ((1<=k).and.(k<=nr).and.(k /= i)) then
tcl = abs(val(j)) / sqrt(abs(diag(i)*diag(k)))
if (abs(val(j)) > theta*sqrt(abs(diag(i)*diag(k)))) then
n = ilaggr(k)
if (n>0) then
if (n>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?')
goto 9999
end if
if ((abs(val(j))>cpling) .or. &
& ((abs(val(j)) == cpling).and. (ils(n) < isz))) then
ia = n
isz = ils(n)
cpling = abs(val(j))
endif
endif
endif
end if
enddo
if (ia == -1) then
! At this point, the easiest thing is to start a new aggregate
naggr = naggr + 1
ilaggr(i) = naggr
ils(ilaggr(i)) = 1
else
ilaggr(i) = ia
if (ia>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr ? ')
goto 9999
end if
ils(ia) = ils(ia) + 1
endif
end if
end do
if (debug_level >= psb_debug_outer_) then
if (recovery) then
write(debug_unit,*) me,' ',trim(name),&
& 'Had to recover from strange situation in loc_aggregate.'
write(debug_unit,*) me,' ',trim(name),&
& 'Perhaps an unsymmetric pattern?'
endif
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 2: number of aggregates ',naggr,sum(ils)
do i=1, naggr
write(debug_unit,*) me,' ',trim(name),&
& 'Size of aggregate ',i,' :',count(ilaggr == i), ils(i)
enddo
write(debug_unit,*) me,' ',trim(name),&
& maxval(ils(1:naggr))
write(debug_unit,*) me,' ',trim(name),&
& 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops'
end if
if (count(ilaggr<0) >0) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Fatal error: some leftovers')
goto 9999
endif
deallocate(ils,neigh,stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
if (naggr > ncol) then
write(0,*) name,'Error : naggr > ncol',naggr,ncol
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Fatal error: naggr>ncol')
goto 9999
end if
call psb_realloc(ncol,ilaggr,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
allocate(nlaggr(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/np,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
nlaggr(:) = 0
nlaggr(me+1) = naggr
call psb_sum(ictxt,nlaggr(1:np))
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine mld_s_dec_map_bld

11
mlprec/mld_s_inner_mod.f90
Unescape Escape View File

@ -102,6 +102,17 @@ module mld_s_inner_mod
end subroutine mld_saggrmap_bld
end interface mld_aggrmap_bld
interface mld_dec_map_bld
subroutine mld_s_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_
type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
real(psb_spk_), intent(in) :: theta
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
end subroutine mld_s_dec_map_bld
end interface mld_dec_map_bld
interface mld_aggrmat_asb
subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_

301
mlprec/mld_saggrmap_bld.f90
Unescape Escape View File

@ -160,305 +160,4 @@ subroutine mld_saggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info)
end if
return
contains
subroutine mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod
use mld_s_inner_mod !, mld_protect_name => mld_daggrmap_bld
implicit none
! Arguments
type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
real(psb_spk_), intent(in) :: theta
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
! Local variables
integer, allocatable :: ils(:), neigh(:), irow(:), icol(:)
real(psb_spk_), allocatable :: val(:), diag(:)
integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz
real(psb_spk_) :: cpling, tcl
logical :: recovery
integer :: debug_level, debug_unit
integer :: ictxt,np,me,err_act
integer :: nrow, ncol, n_ne
character(len=20) :: name, ch_err
if(psb_get_errstatus() /= 0) return
info=psb_success_
name = 'mld_dec_map_bld'
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!
ictxt=psb_cd_get_context(desc_a)
call psb_info(ictxt,me,np)
nrow = psb_cd_get_local_rows(desc_a)
ncol = psb_cd_get_local_cols(desc_a)
nr = a%get_nrows()
allocate(ilaggr(nr),neigh(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nr,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
allocate(diag(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/nr,0,0,0,0/),&
& a_err='real(psb_spk_)')
goto 9999
end if
call a%get_diag(diag,info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_sp_getdiag')
goto 9999
end if
do i=1, nr
ilaggr(i) = -(nr+1)
end do
! Note: -(nr+1) Untouched as yet
! -i 1<=i<=nr Adjacent to aggregate i
! i 1<=i<=nr Belonging to aggregate i
!
! Phase one: group nodes together.
! Very simple minded strategy.
!
naggr = 0
nlp = 0
do
icnt = 0
do i=1, nr
if (ilaggr(i) == -(nr+1)) then
!
! 1. Untouched nodes are marked >0 together
! with their neighbours
!
icnt = icnt + 1
naggr = naggr + 1
ilaggr(i) = naggr
call a%csget(i,i,nz,irow,icol,val,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='csget')
goto 9999
end if
do k=1, nz
j = icol(k)
if ((1<=j).and.(j<=nr).and.(i /= j)) then
if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
ilaggr(j) = naggr
else
ilaggr(j) = -naggr
endif
end if
enddo
!
! 2. Untouched neighbours of these nodes are marked <0.
!
call a%get_neigh(i,neigh,n_ne,info,lev=2)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_neigh')
goto 9999
end if
do n = 1, n_ne
m = neigh(n)
if ((1<=m).and.(m<=nr)) then
if (ilaggr(m) == -(nr+1)) ilaggr(m) = -naggr
endif
enddo
endif
enddo
nlp = nlp + 1
if (icnt == 0) exit
enddo
if (debug_level >= psb_debug_outer_) then
write(debug_unit,*) me,' ',trim(name),&
& ' Check 1:',count(ilaggr == -(nr+1))
end if
!
! Phase two: sweep over leftovers.
!
allocate(ils(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
do i=1, size(ils)
ils(i) = 0
end do
do i=1, nr
n = ilaggr(i)
if (n>0) then
if (n>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 1 ?')
goto 9999
else
ils(n) = ils(n) + 1
end if
end if
end do
if (debug_level >= psb_debug_outer_) then
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 1: number of aggregates ',naggr
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 1: nodes aggregated ',sum(ils)
end if
recovery=.false.
do i=1, nr
if (ilaggr(i) < 0) then
!
! Now some silly rule to break ties:
! Group with adjacent aggregate.
!
isz = nr+1
ia = -1
cpling = szero
call a%csget(i,i,nz,irow,icol,val,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_sp_getrow')
goto 9999
end if
do j=1, nz
k = icol(j)
if ((1<=k).and.(k<=nr).and.(k /= i)) then
tcl = abs(val(j)) / sqrt(abs(diag(i)*diag(k)))
if (abs(val(j)) > theta*sqrt(abs(diag(i)*diag(k)))) then
!!$ if (tcl > theta) then
n = ilaggr(k)
if (n>0) then
if (n>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?')
goto 9999
end if
if ((abs(val(j))>cpling) .or. &
& ((abs(val(j)) == cpling).and. (ils(n) < isz))) then
!!$ if ((tcl > cpling) .or. ((tcl == cpling).and. (ils(n) < isz))) then
ia = n
isz = ils(n)
cpling = abs(val(j))
!!$ cpling = tcl
endif
endif
endif
end if
enddo
if (ia == -1) then
! At this point, the easiest thing is to start a new aggregate
naggr = naggr + 1
ilaggr(i) = naggr
ils(ilaggr(i)) = 1
else
ilaggr(i) = ia
if (ia>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr ? ')
goto 9999
end if
ils(ia) = ils(ia) + 1
endif
end if
end do
if (debug_level >= psb_debug_outer_) then
if (recovery) then
write(debug_unit,*) me,' ',trim(name),&
& 'Had to recover from strange situation in loc_aggregate.'
write(debug_unit,*) me,' ',trim(name),&
& 'Perhaps an unsymmetric pattern?'
endif
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 2: number of aggregates ',naggr,sum(ils)
do i=1, naggr
write(debug_unit,*) me,' ',trim(name),&
& 'Size of aggregate ',i,' :',count(ilaggr == i), ils(i)
enddo
write(debug_unit,*) me,' ',trim(name),&
& maxval(ils(1:naggr))
write(debug_unit,*) me,' ',trim(name),&
& 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops'
end if
if (count(ilaggr<0) >0) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Fatal error: some leftovers')
goto 9999
endif
deallocate(ils,neigh,stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
if (naggr > ncol) then
write(0,*) name,'Error : naggr > ncol',naggr,ncol
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Fatal error: naggr>ncol')
goto 9999
end if
call psb_realloc(ncol,ilaggr,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
allocate(nlaggr(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/np,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
nlaggr(:) = 0
nlaggr(me+1) = naggr
call psb_sum(ictxt,nlaggr(1:np))
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine mld_dec_map_bld
end subroutine mld_saggrmap_bld

296
mlprec/mld_z_dec_map_bld.F90
Unescape Escape View File

@ -0,0 +1,296 @@
subroutine mld_z_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod
use mld_z_inner_mod, mld_protect_name => mld_z_dec_map_bld
implicit none
! Arguments
type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
real(psb_dpk_), intent(in) :: theta
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
! Local variables
integer, allocatable :: ils(:), neigh(:), irow(:), icol(:)
complex(psb_dpk_), allocatable :: val(:), diag(:)
integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz
real(psb_dpk_) :: cpling, tcl
logical :: recovery
integer :: debug_level, debug_unit
integer :: ictxt,np,me,err_act
integer :: nrow, ncol, n_ne
character(len=20) :: name, ch_err
if(psb_get_errstatus() /= 0) return
info=psb_success_
name = 'mld_dec_map_bld'
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!
ictxt=psb_cd_get_context(desc_a)
call psb_info(ictxt,me,np)
nrow = psb_cd_get_local_rows(desc_a)
ncol = psb_cd_get_local_cols(desc_a)
nr = a%get_nrows()
allocate(ilaggr(nr),neigh(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nr,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
allocate(diag(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/nr,0,0,0,0/),&
& a_err='complex(psb_dpk_)')
goto 9999
end if
call a%get_diag(diag,info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_sp_getdiag')
goto 9999
end if
do i=1, nr
ilaggr(i) = -(nr+1)
end do
! Note: -(nr+1) Untouched as yet
! -i 1<=i<=nr Adjacent to aggregate i
! i 1<=i<=nr Belonging to aggregate i
!
! Phase one: group nodes together.
! Very simple minded strategy.
!
naggr = 0
nlp = 0
do
icnt = 0
do i=1, nr
if (ilaggr(i) == -(nr+1)) then
!
! 1. Untouched nodes are marked >0 together
! with their neighbours
!
icnt = icnt + 1
naggr = naggr + 1
ilaggr(i) = naggr
call a%csget(i,i,nz,irow,icol,val,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='csget')
goto 9999
end if
do k=1, nz
j = icol(k)
if ((1<=j).and.(j<=nr).and.(i /= j)) then
if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
ilaggr(j) = naggr
else
ilaggr(j) = -naggr
endif
end if
enddo
!
! 2. Untouched neighbours of these nodes are marked <0.
!
call a%get_neigh(i,neigh,n_ne,info,lev=2)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_neigh')
goto 9999
end if
do n = 1, n_ne
m = neigh(n)
if ((1<=m).and.(m<=nr)) then
if (ilaggr(m) == -(nr+1)) ilaggr(m) = -naggr
endif
enddo
endif
enddo
nlp = nlp + 1
if (icnt == 0) exit
enddo
if (debug_level >= psb_debug_outer_) then
write(debug_unit,*) me,' ',trim(name),&
& ' Check 1:',count(ilaggr == -(nr+1))
end if
!
! Phase two: sweep over leftovers.
!
allocate(ils(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
do i=1, size(ils)
ils(i) = 0
end do
do i=1, nr
n = ilaggr(i)
if (n>0) then
if (n>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 1 ?')
goto 9999
else
ils(n) = ils(n) + 1
end if
end if
end do
if (debug_level >= psb_debug_outer_) then
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 1: number of aggregates ',naggr
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 1: nodes aggregated ',sum(ils)
end if
recovery=.false.
do i=1, nr
if (ilaggr(i) < 0) then
!
! Now some silly rule to break ties:
! Group with adjacent aggregate.
!
isz = nr+1
ia = -1
cpling = dzero
call a%csget(i,i,nz,irow,icol,val,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_sp_getrow')
goto 9999
end if
do j=1, nz
k = icol(j)
if ((1<=k).and.(k<=nr).and.(k /= i)) then
tcl = abs(val(j)) / sqrt(abs(diag(i)*diag(k)))
if (abs(val(j)) > theta*sqrt(abs(diag(i)*diag(k)))) then
n = ilaggr(k)
if (n>0) then
if (n>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?')
goto 9999
end if
if ((abs(val(j))>cpling) .or. &
& ((abs(val(j)) == cpling).and. (ils(n) < isz))) then
ia = n
isz = ils(n)
cpling = abs(val(j))
endif
endif
endif
end if
enddo
if (ia == -1) then
! At this point, the easiest thing is to start a new aggregate
naggr = naggr + 1
ilaggr(i) = naggr
ils(ilaggr(i)) = 1
else
ilaggr(i) = ia
if (ia>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr ? ')
goto 9999
end if
ils(ia) = ils(ia) + 1
endif
end if
end do
if (debug_level >= psb_debug_outer_) then
if (recovery) then
write(debug_unit,*) me,' ',trim(name),&
& 'Had to recover from strange situation in loc_aggregate.'
write(debug_unit,*) me,' ',trim(name),&
& 'Perhaps an unsymmetric pattern?'
endif
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 2: number of aggregates ',naggr,sum(ils)
do i=1, naggr
write(debug_unit,*) me,' ',trim(name),&
& 'Size of aggregate ',i,' :',count(ilaggr == i), ils(i)
enddo
write(debug_unit,*) me,' ',trim(name),&
& maxval(ils(1:naggr))
write(debug_unit,*) me,' ',trim(name),&
& 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops'
end if
if (count(ilaggr<0) >0) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Fatal error: some leftovers')
goto 9999
endif
deallocate(ils,neigh,stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
if (naggr > ncol) then
write(0,*) name,'Error : naggr > ncol',naggr,ncol
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Fatal error: naggr>ncol')
goto 9999
end if
call psb_realloc(ncol,ilaggr,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
allocate(nlaggr(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/np,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
nlaggr(:) = 0
nlaggr(me+1) = naggr
call psb_sum(ictxt,nlaggr(1:np))
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine mld_z_dec_map_bld

12
mlprec/mld_z_inner_mod.f90
Unescape Escape View File

@ -102,6 +102,18 @@ module mld_z_inner_mod
end subroutine mld_zaggrmap_bld
end interface mld_aggrmap_bld
interface mld_dec_map_bld
subroutine mld_z_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_
type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
real(psb_dpk_), intent(in) :: theta
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
end subroutine mld_z_dec_map_bld
end interface mld_dec_map_bld
interface mld_aggrmat_asb
subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_

301
mlprec/mld_zaggrmap_bld.f90
Unescape Escape View File

@ -160,305 +160,4 @@ subroutine mld_zaggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info)
end if
return
contains
subroutine mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod
use mld_z_inner_mod !, mld_protect_name => mld_daggrmap_bld
implicit none
! Arguments
type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
real(psb_dpk_), intent(in) :: theta
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
! Local variables
integer, allocatable :: ils(:), neigh(:), irow(:), icol(:)
complex(psb_dpk_), allocatable :: val(:), diag(:)
integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz
real(psb_dpk_) :: cpling, tcl
logical :: recovery
integer :: debug_level, debug_unit
integer :: ictxt,np,me,err_act
integer :: nrow, ncol, n_ne
character(len=20) :: name, ch_err
if(psb_get_errstatus() /= 0) return
info=psb_success_
name = 'mld_dec_map_bld'
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!
ictxt=psb_cd_get_context(desc_a)
call psb_info(ictxt,me,np)
nrow = psb_cd_get_local_rows(desc_a)
ncol = psb_cd_get_local_cols(desc_a)
nr = a%get_nrows()
allocate(ilaggr(nr),neigh(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nr,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
allocate(diag(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/nr,0,0,0,0/),&
& a_err='complex(psb_dpk_)')
goto 9999
end if
call a%get_diag(diag,info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_sp_getdiag')
goto 9999
end if
do i=1, nr
ilaggr(i) = -(nr+1)
end do
! Note: -(nr+1) Untouched as yet
! -i 1<=i<=nr Adjacent to aggregate i
! i 1<=i<=nr Belonging to aggregate i
!
! Phase one: group nodes together.
! Very simple minded strategy.
!
naggr = 0
nlp = 0
do
icnt = 0
do i=1, nr
if (ilaggr(i) == -(nr+1)) then
!
! 1. Untouched nodes are marked >0 together
! with their neighbours
!
icnt = icnt + 1
naggr = naggr + 1
ilaggr(i) = naggr
call a%csget(i,i,nz,irow,icol,val,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='csget')
goto 9999
end if
do k=1, nz
j = icol(k)
if ((1<=j).and.(j<=nr).and.(i /= j)) then
if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
ilaggr(j) = naggr
else
ilaggr(j) = -naggr
endif
end if
enddo
!
! 2. Untouched neighbours of these nodes are marked <0.
!
call a%get_neigh(i,neigh,n_ne,info,lev=2)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_neigh')
goto 9999
end if
do n = 1, n_ne
m = neigh(n)
if ((1<=m).and.(m<=nr)) then
if (ilaggr(m) == -(nr+1)) ilaggr(m) = -naggr
endif
enddo
endif
enddo
nlp = nlp + 1
if (icnt == 0) exit
enddo
if (debug_level >= psb_debug_outer_) then
write(debug_unit,*) me,' ',trim(name),&
& ' Check 1:',count(ilaggr == -(nr+1))
end if
!
! Phase two: sweep over leftovers.
!
allocate(ils(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
do i=1, size(ils)
ils(i) = 0
end do
do i=1, nr
n = ilaggr(i)
if (n>0) then
if (n>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 1 ?')
goto 9999
else
ils(n) = ils(n) + 1
end if
end if
end do
if (debug_level >= psb_debug_outer_) then
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 1: number of aggregates ',naggr
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 1: nodes aggregated ',sum(ils)
end if
recovery=.false.
do i=1, nr
if (ilaggr(i) < 0) then
!
! Now some silly rule to break ties:
! Group with adjacent aggregate.
!
isz = nr+1
ia = -1
cpling = dzero
call a%csget(i,i,nz,irow,icol,val,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_sp_getrow')
goto 9999
end if
do j=1, nz
k = icol(j)
if ((1<=k).and.(k<=nr).and.(k /= i)) then
tcl = abs(val(j)) / sqrt(abs(diag(i)*diag(k)))
if (abs(val(j)) > theta*sqrt(abs(diag(i)*diag(k)))) then
!!$ if (tcl > theta) then
n = ilaggr(k)
if (n>0) then
if (n>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?')
goto 9999
end if
if ((abs(val(j))>cpling) .or. &
& ((abs(val(j)) == cpling).and. (ils(n) < isz))) then
!!$ if ((tcl > cpling) .or. ((tcl == cpling).and. (ils(n) < isz))) then
ia = n
isz = ils(n)
cpling = abs(val(j))
!!$ cpling = tcl
endif
endif
endif
end if
enddo
if (ia == -1) then
! At this point, the easiest thing is to start a new aggregate
naggr = naggr + 1
ilaggr(i) = naggr
ils(ilaggr(i)) = 1
else
ilaggr(i) = ia
if (ia>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr ? ')
goto 9999
end if
ils(ia) = ils(ia) + 1
endif
end if
end do
if (debug_level >= psb_debug_outer_) then
if (recovery) then
write(debug_unit,*) me,' ',trim(name),&
& 'Had to recover from strange situation in loc_aggregate.'
write(debug_unit,*) me,' ',trim(name),&
& 'Perhaps an unsymmetric pattern?'
endif
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 2: number of aggregates ',naggr,sum(ils)
do i=1, naggr
write(debug_unit,*) me,' ',trim(name),&
& 'Size of aggregate ',i,' :',count(ilaggr == i), ils(i)
enddo
write(debug_unit,*) me,' ',trim(name),&
& maxval(ils(1:naggr))
write(debug_unit,*) me,' ',trim(name),&
& 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops'
end if
if (count(ilaggr<0) >0) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Fatal error: some leftovers')
goto 9999
endif
deallocate(ils,neigh,stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
if (naggr > ncol) then
write(0,*) name,'Error : naggr > ncol',naggr,ncol
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Fatal error: naggr>ncol')
goto 9999
end if
call psb_realloc(ncol,ilaggr,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
allocate(nlaggr(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/np,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
nlaggr(:) = 0
nlaggr(me+1) = naggr
call psb_sum(ictxt,nlaggr(1:np))
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine mld_dec_map_bld
end subroutine mld_zaggrmap_bld

Write Preview
Loading…
Cancel
Save