fdurastante/amg4psblas
1
0
Fork 0
mirror of https://github.com/sfilippone/amg4psblas.git synced 2026-10-06 22:55:12 +00:00
Code Issues Packages Projects Releases Wiki Activity

mld2p4:

Browse Source 
mlprec/mld_caggrmap_bld.f90
 mlprec/mld_cmlprec_bld.f90
 mlprec/mld_cprecinit.f90
 mlprec/mld_cprecset.f90
 mlprec/mld_daggrmap_bld.f90
 mlprec/mld_dmlprec_bld.f90
 mlprec/mld_dprecinit.f90
 mlprec/mld_dprecset.f90
 mlprec/mld_inner_mod.f90
 mlprec/mld_prec_type.f90
 mlprec/mld_saggrmap_bld.f90
 mlprec/mld_smlprec_bld.f90
 mlprec/mld_sprecinit.f90
 mlprec/mld_sprecset.f90
 mlprec/mld_zaggrmap_bld.f90
 mlprec/mld_zmlprec_bld.f90
 mlprec/mld_zprecinit.f90
 mlprec/mld_zprecset.f90
 test/fileread/cf_sample.f90
 test/fileread/df_bench.f90
 test/fileread/df_sample.f90
 test/fileread/runs/cfs.inp
 test/fileread/runs/dfs.inp
 test/fileread/runs/drt.sh
 test/fileread/runs/sfs.inp
 test/fileread/runs/zfs.inp
 test/fileread/sf_sample.f90
 test/fileread/zf_sample.f90


Added an aggregation threshold to account for anisotropies.
This commit is contained in:
Salvatore Filippone
2008-06-06 08:54:48 +00:00
parent b7d9b380d5
commit b9901f1f8e
28 changed files with 848 additions and 433 deletions
Download Patch File Download Diff File Expand all files Collapse all files
mlprec/mld_caggrmap_bld.f90
+175 -97
View File
@@ -74,7 +74,7 @@
! info - integer, output.
! Error code.
!
subroutine mld_caggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
subroutine mld_caggrmap_bld(aggr_type,theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod
use mld_inner_mod, mld_protect_name => mld_caggrmap_bld
@@ -83,7 +83,8 @@ subroutine mld_caggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
! Arguments
integer, intent(in) :: aggr_type
type(psb_cspmat_type), intent(in), target :: a
real(psb_spk_), intent(in) :: theta
type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
@@ -91,9 +92,7 @@ subroutine mld_caggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
! Local variables
integer, allocatable :: ils(:), neigh(:)
integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m
type(psb_cspmat_type), target :: atmp, atrans
type(psb_cspmat_type), pointer :: apnt
type(psb_cspmat_type) :: atmp, atrans
logical :: recovery
integer :: debug_level, debug_unit
integer :: ictxt,np,me,err_act
@@ -117,8 +116,93 @@ subroutine mld_caggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
ncol = psb_cd_get_local_cols(desc_a)
select case (aggr_type)
case (mld_dec_aggr_,mld_sym_dec_aggr_)
case (mld_dec_aggr_)
call mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info)
case (mld_sym_dec_aggr_)
nr = psb_sp_get_nrows(a)
call psb_sp_clip(a,atmp,info,imax=nr,jmax=nr,&
& rscale=.false.,cscale=.false.)
atmp%m=nr
atmp%k=nr
if (info == 0) call psb_transp(atmp,atrans,fmt='COO')
if (info == 0) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.)
atmp%m=nr
atmp%k=nr
if (info == 0) call psb_sp_free(atrans,info)
if (info == 0) call psb_ipcoo2csr(atmp,info)
if (info == 0) call mld_dec_map_bld(theta,atmp,desc_a,nlaggr,ilaggr,info)
if (info == 0) call psb_sp_free(atmp,info)
case default
info = -1
call psb_errpush(30,name,i_err=(/1,aggr_type,0,0,0/))
goto 9999
end select
if (info/=0) then
info=4001
call psb_errpush(info,name,a_err='dec_map_bld')
goto 9999
end if
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
contains
subroutine mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod
use mld_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=0
name = 'mld_dec_map_bld'
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!
! Note. At the time being we are ignoring aggr_type so
! that we only have decoupled aggregation. This might
! change in the future.
!
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%m
allocate(ilaggr(nr),neigh(nr),stat=info)
if(info /= 0) then
@@ -128,36 +212,28 @@ subroutine mld_caggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
goto 9999
end if
allocate(diag(nr),stat=info)
if(info /= 0) then
info=4025
call psb_errpush(info,name,i_err=(/nr,0,0,0,0/),&
& a_err='complex(psb_spk_)')
goto 9999
end if
call psb_sp_getdiag(a,diag,info)
if(info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_sp_getdiag')
goto 9999
end if
do i=1, nr
ilaggr(i) = -(nr+1)
end do
if (aggr_type == mld_dec_aggr_) then
apnt => a
else
call psb_sp_clip(a,atmp,info,imax=nr,jmax=nr,&
& rscale=.false.,cscale=.false.)
atmp%m=nr
atmp%k=nr
if (info == 0) call psb_transp(atmp,atrans,fmt='COO')
if (info == 0) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.)
atmp%m=nr
atmp%k=nr
if (info == 0) call psb_sp_free(atrans,info)
if (info == 0) call psb_ipcoo2csr(atmp,info)
apnt => atmp
if (info/=0) then
info=4001
call psb_errpush(info,name,a_err='init apnt')
goto 9999
end if
end if
! 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.
@@ -172,26 +248,32 @@ subroutine mld_caggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
! 1. Untouched nodes are marked >0 together
! with their neighbours
!
icnt = icnt + 1
naggr = naggr + 1
icnt = icnt + 1
naggr = naggr + 1
ilaggr(i) = naggr
call psb_neigh(apnt,i,neigh,n_ne,info,lev=1)
call psb_sp_getrow(i,a,nz,irow,icol,val,info)
if (info/=0) then
info=4010
call psb_errpush(info,name,a_err='psb_neigh')
call psb_errpush(info,name,a_err='psb_sp_getrow')
goto 9999
end if
do k=1, n_ne
j = neigh(k)
if ((1<=j).and.(j<=nr)) then
ilaggr(j) = naggr
endif
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 psb_neigh(apnt,i,neigh,n_ne,info,lev=2)
call psb_neigh(a,i,neigh,n_ne,info,lev=2)
if (info/=0) then
info=4010
call psb_errpush(info,name,a_err='psb_neigh')
@@ -218,7 +300,7 @@ subroutine mld_caggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
!
! Phase two: sweep over leftovers.
!
allocate(ils(naggr+10),stat=info)
allocate(ils(nr),stat=info)
if(info /= 0) then
info=4025
call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),&
@@ -254,57 +336,55 @@ subroutine mld_caggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
if (ilaggr(i) < 0) then
!
! Now some silly rule to break ties:
! Group with smallest adjacent aggregate.
! Group with adjacent aggregate.
!
isz = nr+1
ia = -1
call psb_neigh(apnt,i,neigh,n_ne,info,lev=1)
isz = nr+1
ia = -1
cpling = 0.0
call psb_sp_getrow(i,a,nz,irow,icol,val,info)
if (info/=0) then
info=4010
call psb_errpush(info,name,a_err='psb_neigh')
call psb_errpush(info,name,a_err='psb_sp_getrow')
goto 9999
end if
do j=1, n_ne
k = neigh(j)
if ((1<=k).and.(k<=nr)) then
n = ilaggr(k)
if (n>0) then
if (n>naggr) then
info=4001
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?')
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=4001
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?')
goto 9999
end if
if (ils(n) < isz) then
ia = n
isz = ils(n)
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
endif
end if
enddo
if (ia == -1) then
if (ilaggr(i) > -(nr+1)) then
ilaggr(i) = abs(ilaggr(i))
if (ilaggr(I)>naggr) then
info=4001
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 3 ?')
goto 9999
end if
ils(ilaggr(i)) = ils(ilaggr(i)) + 1
!
! This might happen if the pattern is non symmetric.
! Need a better handling.
!
recovery = .true.
else
info=4001
call psb_errpush(info,name,a_err='Unrecoverable error !!')
goto 9999
endif
! 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=4001
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 4? ')
@@ -312,8 +392,9 @@ subroutine mld_caggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
end if
ils(ia) = ils(ia) + 1
endif
end if
enddo
end do
if (debug_level >= psb_debug_outer_) then
if (recovery) then
write(debug_unit,*) me,' ',trim(name),&
@@ -345,6 +426,12 @@ subroutine mld_caggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
call psb_errpush(info,name)
goto 9999
end if
if (naggr > ncol) then
write(0,*) name,'Error : naggr > ncol',naggr,ncol
info=4001
call psb_errpush(info,name,a_err='Fatal error: naggr>ncol')
goto 9999
end if
call psb_realloc(ncol,ilaggr,info)
if (info/=0) then
@@ -366,26 +453,17 @@ subroutine mld_caggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
nlaggr(me+1) = naggr
call psb_sum(ictxt,nlaggr(1:np))
if (aggr_type == mld_sym_dec_aggr_) then
call psb_sp_free(atmp,info)
end if
case default
info = -1
call psb_errpush(30,name,i_err=(/1,aggr_type,0,0,0/))
goto 9999
end select
call psb_erractionrestore(err_act)
return
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()
return
end if
return
end if
return
end subroutine mld_dec_map_bld
end subroutine mld_caggrmap_bld
mlprec/mld_cmlprec_bld.f90
+3 -1
View File
@@ -109,6 +109,7 @@ subroutine mld_cmlprec_bld(a,desc_a,p,info)
call mld_check_def(p%rprcparm(mld_fact_thrs_),'Eps',szero,is_legal_s_fact_thrs)
end select
call mld_check_def(p%rprcparm(mld_aggr_damp_),'Omega',szero,is_legal_s_omega)
call mld_check_def(p%rprcparm(mld_aggr_thresh_),'Aggr_Thresh',szero,is_legal_s_aggr_thrs)
call mld_check_def(p%iprcparm(mld_smooth_sweeps_),'Jacobi sweeps',&
& 1,is_legal_jac_sweeps)
@@ -118,7 +119,8 @@ subroutine mld_cmlprec_bld(a,desc_a,p,info)
! aggregation algorithm. This also defines a tentative prolongator from
! the coarse to the fine level.
!
call mld_aggrmap_bld(p%iprcparm(mld_aggr_alg_),a,desc_a,p%nlaggr,p%mlia,info)
call mld_aggrmap_bld(p%iprcparm(mld_aggr_alg_),p%rprcparm(mld_aggr_thresh_),&
& a,desc_a,p%nlaggr,p%mlia,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='mld_aggrmap_bld')
goto 9999
mlprec/mld_cprecinit.f90
+6 -3
View File
@@ -189,7 +189,8 @@ subroutine mld_cprecinit(p,ptype,info,nlev)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%rprcparm(:) = szero
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_as_
p%baseprecv(ilev_)%iprcparm(mld_sub_solve_) = mld_ilu_n_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_halo_
@@ -204,7 +205,8 @@ subroutine mld_cprecinit(p,ptype,info,nlev)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%rprcparm(:) = szero
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_bjac_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_
@@ -225,7 +227,8 @@ subroutine mld_cprecinit(p,ptype,info,nlev)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%rprcparm(:) = szero
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_bjac_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_
mlprec/mld_cprecset.f90
+12 -3
View File
@@ -595,7 +595,7 @@ subroutine mld_cprecsetr(p,what,val,info,ilev)
else if (ilev_ > 1) then
select case(what)
case(mld_aggr_damp_,mld_fact_thrs_)
case(mld_aggr_damp_,mld_aggr_thresh_,mld_fact_thrs_)
p%baseprecv(ilev_)%rprcparm(what) = val
case default
write(0,*) name,': Error: invalid WHAT'
@@ -610,7 +610,7 @@ subroutine mld_cprecsetr(p,what,val,info,ilev)
select case(what)
case(mld_fact_thrs_)
do ilev_=1,nlev_-1
do ilev_=1,nlev_
if (.not.allocated(p%baseprecv(ilev_)%rprcparm)) then
write(0,*) name,': Error: Uninitialized preconditioner component, should call MLD_PRECINIT'
info = -1
@@ -619,7 +619,16 @@ subroutine mld_cprecsetr(p,what,val,info,ilev)
p%baseprecv(ilev_)%rprcparm(what) = val
end do
case(mld_aggr_damp_)
do ilev_=2,nlev_-1
do ilev_=2,nlev_
if (.not.allocated(p%baseprecv(ilev_)%rprcparm)) then
write(0,*) name,': Error: Uninitialized preconditioner component, should call MLD_PRECINIT'
info = -1
return
endif
p%baseprecv(ilev_)%rprcparm(what) = val
end do
case(mld_aggr_thresh_)
do ilev_=2,nlev_
if (.not.allocated(p%baseprecv(ilev_)%rprcparm)) then
write(0,*) name,': Error: Uninitialized preconditioner component, should call MLD_PRECINIT'
info = -1
mlprec/mld_daggrmap_bld.f90
+176 -97
View File
@@ -74,16 +74,17 @@
! info - integer, output.
! Error code.
!
subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
subroutine mld_daggrmap_bld(aggr_type,theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod
use mld_inner_mod, mld_protect_name => mld_daggrmap_bld
implicit none
! Arguments
integer, intent(in) :: aggr_type
type(psb_dspmat_type), intent(in), target :: a
real(psb_dpk_), intent(in) :: theta
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
@@ -91,8 +92,7 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
! Local variables
integer, allocatable :: ils(:), neigh(:)
integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m
type(psb_dspmat_type), target :: atmp, atrans
type(psb_dspmat_type), pointer :: apnt
type(psb_dspmat_type) :: atmp, atrans
logical :: recovery
integer :: debug_level, debug_unit
integer :: ictxt,np,me,err_act
@@ -116,8 +116,93 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
ncol = psb_cd_get_local_cols(desc_a)
select case (aggr_type)
case (mld_dec_aggr_,mld_sym_dec_aggr_)
case (mld_dec_aggr_)
call mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info)
case (mld_sym_dec_aggr_)
nr = psb_sp_get_nrows(a)
call psb_sp_clip(a,atmp,info,imax=nr,jmax=nr,&
& rscale=.false.,cscale=.false.)
atmp%m=nr
atmp%k=nr
if (info == 0) call psb_transp(atmp,atrans,fmt='COO')
if (info == 0) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.)
atmp%m=nr
atmp%k=nr
if (info == 0) call psb_sp_free(atrans,info)
if (info == 0) call psb_ipcoo2csr(atmp,info)
if (info == 0) call mld_dec_map_bld(theta,atmp,desc_a,nlaggr,ilaggr,info)
if (info == 0) call psb_sp_free(atmp,info)
case default
info = -1
call psb_errpush(30,name,i_err=(/1,aggr_type,0,0,0/))
goto 9999
end select
if (info/=0) then
info=4001
call psb_errpush(info,name,a_err='dec_map_bld')
goto 9999
end if
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
contains
subroutine mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod
use mld_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=0
name = 'mld_dec_map_bld'
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!
! Note. At the time being we are ignoring aggr_type so
! that we only have decoupled aggregation. This might
! change in the future.
!
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%m
allocate(ilaggr(nr),neigh(nr),stat=info)
if(info /= 0) then
@@ -127,36 +212,28 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
goto 9999
end if
allocate(diag(nr),stat=info)
if(info /= 0) then
info=4025
call psb_errpush(info,name,i_err=(/nr,0,0,0,0/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
call psb_sp_getdiag(a,diag,info)
if(info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_sp_getdiag')
goto 9999
end if
do i=1, nr
ilaggr(i) = -(nr+1)
end do
if (aggr_type == mld_dec_aggr_) then
apnt => a
else
call psb_sp_clip(a,atmp,info,imax=nr,jmax=nr,&
& rscale=.false.,cscale=.false.)
atmp%m=nr
atmp%k=nr
if (info == 0) call psb_transp(atmp,atrans,fmt='COO')
if (info == 0) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.)
atmp%m=nr
atmp%k=nr
if (info == 0) call psb_sp_free(atrans,info)
if (info == 0) call psb_ipcoo2csr(atmp,info)
apnt => atmp
if (info/=0) then
info=4001
call psb_errpush(info,name,a_err='init apnt')
goto 9999
end if
end if
! 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.
@@ -171,26 +248,32 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
! 1. Untouched nodes are marked >0 together
! with their neighbours
!
icnt = icnt + 1
naggr = naggr + 1
icnt = icnt + 1
naggr = naggr + 1
ilaggr(i) = naggr
call psb_neigh(apnt,i,neigh,n_ne,info,lev=1)
call psb_sp_getrow(i,a,nz,irow,icol,val,info)
if (info/=0) then
info=4010
call psb_errpush(info,name,a_err='psb_neigh')
call psb_errpush(info,name,a_err='psb_sp_getrow')
goto 9999
end if
do k=1, n_ne
j = neigh(k)
if ((1<=j).and.(j<=nr)) then
ilaggr(j) = naggr
endif
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 psb_neigh(apnt,i,neigh,n_ne,info,lev=2)
call psb_neigh(a,i,neigh,n_ne,info,lev=2)
if (info/=0) then
info=4010
call psb_errpush(info,name,a_err='psb_neigh')
@@ -217,7 +300,7 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
!
! Phase two: sweep over leftovers.
!
allocate(ils(naggr+10),stat=info)
allocate(ils(nr),stat=info)
if(info /= 0) then
info=4025
call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),&
@@ -253,57 +336,55 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
if (ilaggr(i) < 0) then
!
! Now some silly rule to break ties:
! Group with smallest adjacent aggregate.
! Group with adjacent aggregate.
!
isz = nr+1
ia = -1
call psb_neigh(apnt,i,neigh,n_ne,info,lev=1)
isz = nr+1
ia = -1
cpling = 0.0d0
call psb_sp_getrow(i,a,nz,irow,icol,val,info)
if (info/=0) then
info=4010
call psb_errpush(info,name,a_err='psb_neigh')
call psb_errpush(info,name,a_err='psb_sp_getrow')
goto 9999
end if
do j=1, n_ne
k = neigh(j)
if ((1<=k).and.(k<=nr)) then
n = ilaggr(k)
if (n>0) then
if (n>naggr) then
info=4001
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?')
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=4001
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?')
goto 9999
end if
if (ils(n) < isz) then
ia = n
isz = ils(n)
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
endif
end if
enddo
if (ia == -1) then
if (ilaggr(i) > -(nr+1)) then
ilaggr(i) = abs(ilaggr(i))
if (ilaggr(I)>naggr) then
info=4001
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 3 ?')
goto 9999
end if
ils(ilaggr(i)) = ils(ilaggr(i)) + 1
!
! This might happen if the pattern is non symmetric.
! Need a better handling.
!
recovery = .true.
else
info=4001
call psb_errpush(info,name,a_err='Unrecoverable error !!')
goto 9999
endif
! 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=4001
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 4? ')
@@ -311,8 +392,9 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
end if
ils(ia) = ils(ia) + 1
endif
end if
enddo
end do
if (debug_level >= psb_debug_outer_) then
if (recovery) then
write(debug_unit,*) me,' ',trim(name),&
@@ -344,6 +426,12 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
call psb_errpush(info,name)
goto 9999
end if
if (naggr > ncol) then
write(0,*) name,'Error : naggr > ncol',naggr,ncol
info=4001
call psb_errpush(info,name,a_err='Fatal error: naggr>ncol')
goto 9999
end if
call psb_realloc(ncol,ilaggr,info)
if (info/=0) then
@@ -365,26 +453,17 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
nlaggr(me+1) = naggr
call psb_sum(ictxt,nlaggr(1:np))
if (aggr_type == mld_sym_dec_aggr_) then
call psb_sp_free(atmp,info)
end if
case default
info = -1
call psb_errpush(30,name,i_err=(/1,aggr_type,0,0,0/))
goto 9999
end select
call psb_erractionrestore(err_act)
return
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()
return
end if
return
end if
return
end subroutine mld_dec_map_bld
end subroutine mld_daggrmap_bld
mlprec/mld_dmlprec_bld.f90
+3 -1
View File
@@ -109,6 +109,7 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info)
call mld_check_def(p%rprcparm(mld_fact_thrs_),'Eps',dzero,is_legal_fact_thrs)
end select
call mld_check_def(p%rprcparm(mld_aggr_damp_),'Omega',dzero,is_legal_omega)
call mld_check_def(p%rprcparm(mld_aggr_thresh_),'Aggr_Thresh',dzero,is_legal_aggr_thrs)
call mld_check_def(p%iprcparm(mld_smooth_sweeps_),'Jacobi sweeps',&
& 1,is_legal_jac_sweeps)
@@ -118,7 +119,8 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info)
! aggregation algorithm. This also defines a tentative prolongator from
! the coarse to the fine level.
!
call mld_aggrmap_bld(p%iprcparm(mld_aggr_alg_),a,desc_a,p%nlaggr,p%mlia,info)
call mld_aggrmap_bld(p%iprcparm(mld_aggr_alg_),p%rprcparm(mld_aggr_thresh_),&
& a,desc_a,p%nlaggr,p%mlia,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='mld_aggrmap_bld')
goto 9999
mlprec/mld_dprecinit.f90
+6 -3
View File
@@ -189,7 +189,8 @@ subroutine mld_dprecinit(p,ptype,info,nlev)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%rprcparm(:) = dzero
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_as_
p%baseprecv(ilev_)%iprcparm(mld_sub_solve_) = mld_ilu_n_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_halo_
@@ -204,7 +205,8 @@ subroutine mld_dprecinit(p,ptype,info,nlev)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%rprcparm(:) = dzero
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_bjac_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_
@@ -225,7 +227,8 @@ subroutine mld_dprecinit(p,ptype,info,nlev)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%rprcparm(:) = dzero
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_bjac_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_
mlprec/mld_dprecset.f90
+12 -3
View File
@@ -595,7 +595,7 @@ subroutine mld_dprecsetr(p,what,val,info,ilev)
else if (ilev_ > 1) then
select case(what)
case(mld_aggr_damp_,mld_fact_thrs_)
case(mld_aggr_damp_,mld_aggr_thresh_,mld_fact_thrs_)
p%baseprecv(ilev_)%rprcparm(what) = val
case default
write(0,*) name,': Error: invalid WHAT'
@@ -610,7 +610,7 @@ subroutine mld_dprecsetr(p,what,val,info,ilev)
select case(what)
case(mld_fact_thrs_)
do ilev_=1,nlev_-1
do ilev_=1,nlev_
if (.not.allocated(p%baseprecv(ilev_)%rprcparm)) then
write(0,*) name,': Error: Uninitialized preconditioner component, should call MLD_PRECINIT'
info = -1
@@ -619,7 +619,16 @@ subroutine mld_dprecsetr(p,what,val,info,ilev)
p%baseprecv(ilev_)%rprcparm(what) = val
end do
case(mld_aggr_damp_)
do ilev_=2,nlev_-1
do ilev_=2,nlev_
if (.not.allocated(p%baseprecv(ilev_)%rprcparm)) then
write(0,*) name,': Error: Uninitialized preconditioner component, should call MLD_PRECINIT'
info = -1
return
endif
p%baseprecv(ilev_)%rprcparm(what) = val
end do
case(mld_aggr_thresh_)
do ilev_=2,nlev_
if (.not.allocated(p%baseprecv(ilev_)%rprcparm)) then
write(0,*) name,': Error: Uninitialized preconditioner component, should call MLD_PRECINIT'
info = -1
mlprec/mld_inner_mod.f90
+12 -8
View File
@@ -384,38 +384,42 @@ module mld_inner_mod
end interface
interface mld_aggrmap_bld
subroutine mld_saggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
subroutine mld_saggrmap_bld(aggr_type,theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_
use mld_prec_type, only : mld_sbaseprc_type
integer, intent(in) :: aggr_type
type(psb_sspmat_type), intent(in), target :: a
real(psb_spk_), intent(in) :: theta
type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
end subroutine mld_saggrmap_bld
subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
subroutine mld_daggrmap_bld(aggr_type,theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_
use mld_prec_type, only : mld_dbaseprc_type
integer, intent(in) :: aggr_type
type(psb_dspmat_type), intent(in), target :: a
real(psb_dpk_), intent(in) :: theta
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
end subroutine mld_daggrmap_bld
subroutine mld_caggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
subroutine mld_caggrmap_bld(aggr_type,theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_
use mld_prec_type, only : mld_cbaseprc_type
integer, intent(in) :: aggr_type
type(psb_cspmat_type), intent(in), target :: a
real(psb_spk_), intent(in) :: theta
type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
end subroutine mld_caggrmap_bld
subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
subroutine mld_zaggrmap_bld(aggr_type,theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_
use mld_prec_type, only : mld_zbaseprc_type
integer, intent(in) :: aggr_type
type(psb_zspmat_type), intent(in), target :: a
real(psb_dpk_), intent(in) :: theta
type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
mlprec/mld_prec_type.f90
+27 -3
View File
@@ -66,7 +66,9 @@ module mld_prec_type
& psb_dspmat_type, psb_zspmat_type,&
& psb_sspmat_type, psb_cspmat_type,&
& psb_desc_type, psb_inter_desc_type, psb_sizeof, psb_dpk_, psb_spk_,&
& psb_sp_free, psb_cdfree
& psb_sp_free, psb_cdfree, psb_halo_, psb_none_, &
& psb_nohalo_, psb_square_root_, &
& psb_sizeof_int, psb_sizeof_sp, psb_sizeof_dp, psb_sizeof
!
! Type: mld_dprec_type, mld_zprec_type
@@ -718,6 +720,8 @@ contains
& aggr_names(p%baseprecv(ilev)%iprcparm(mld_aggr_alg_))
write(iout,*) 'Aggregation smoothing: ', &
& aggr_kinds(p%baseprecv(ilev)%iprcparm(mld_aggr_kind_))
write(iout,*) 'Aggregation threshold: ', &
& p%baseprecv(ilev)%rprcparm(mld_aggr_thresh_)
if (p%baseprecv(ilev)%iprcparm(mld_aggr_kind_) /= mld_no_smooth_) then
write(iout,*) 'Damping omega: ', &
& p%baseprecv(ilev)%rprcparm(mld_aggr_damp_)
@@ -822,6 +826,8 @@ contains
& aggr_names(p%baseprecv(ilev)%iprcparm(mld_aggr_alg_))
write(iout,*) 'Aggregation smoothing: ', &
& aggr_kinds(p%baseprecv(ilev)%iprcparm(mld_aggr_kind_))
write(iout,*) 'Aggregation threshold: ', &
& p%baseprecv(ilev)%rprcparm(mld_aggr_thresh_)
if (p%baseprecv(ilev)%iprcparm(mld_aggr_kind_) /= mld_no_smooth_) then
write(iout,*) 'Damping omega: ', &
& p%baseprecv(ilev)%rprcparm(mld_aggr_damp_)
@@ -944,8 +950,10 @@ contains
if (p%baseprecv(ilev)%iprcparm(mld_ml_type_)>mld_no_ml_) then
write(iout,*) 'Multilevel aggregation: ', &
& aggr_names(p%baseprecv(ilev)%iprcparm(mld_aggr_alg_))
write(iout,*) 'Smoother: ', &
write(iout,*) 'Aggregation smoothing: ', &
& aggr_kinds(p%baseprecv(ilev)%iprcparm(mld_aggr_kind_))
write(iout,*) 'Aggregation threshold: ', &
& p%baseprecv(ilev)%rprcparm(mld_aggr_thresh_)
if (p%baseprecv(ilev)%iprcparm(mld_aggr_kind_) /= mld_no_smooth_) then
write(iout,*) 'Smoothing omega: ', &
& p%baseprecv(ilev)%rprcparm(mld_aggr_damp_)
@@ -1047,8 +1055,10 @@ contains
if (p%baseprecv(ilev)%iprcparm(mld_ml_type_)>mld_no_ml_) then
write(iout,*) 'Multilevel aggregation: ', &
& aggr_names(p%baseprecv(ilev)%iprcparm(mld_aggr_alg_))
write(iout,*) 'Smoother: ', &
write(iout,*) 'Aggregation smoothing: ', &
& aggr_kinds(p%baseprecv(ilev)%iprcparm(mld_aggr_kind_))
write(iout,*) 'Aggregation threshold: ', &
& p%baseprecv(ilev)%rprcparm(mld_aggr_thresh_)
if (p%baseprecv(ilev)%iprcparm(mld_aggr_kind_) /= mld_no_smooth_) then
write(iout,*) 'Smoothing omega: ', &
& p%baseprecv(ilev)%rprcparm(mld_aggr_damp_)
@@ -1205,6 +1215,13 @@ contains
is_legal_fact_thrs = (ip>=0.0d0)
return
end function is_legal_fact_thrs
function is_legal_aggr_thrs(ip)
real(psb_dpk_), intent(in) :: ip
logical :: is_legal_aggr_thrs
is_legal_aggr_thrs = (ip>=0.0d0)
return
end function is_legal_aggr_thrs
function is_legal_s_omega(ip)
real(psb_spk_), intent(in) :: ip
@@ -1219,6 +1236,13 @@ contains
is_legal_s_fact_thrs = (ip>=0.0)
return
end function is_legal_s_fact_thrs
function is_legal_s_aggr_thrs(ip)
real(psb_spk_), intent(in) :: ip
logical :: is_legal_s_aggr_thrs
is_legal_s_aggr_thrs = (ip>=0.0)
return
end function is_legal_s_aggr_thrs
subroutine mld_icheck_def(ip,name,id,is_legal)
mlprec/mld_saggrmap_bld.f90
+175 -96
View File
@@ -74,7 +74,7 @@
! info - integer, output.
! Error code.
!
subroutine mld_saggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
subroutine mld_saggrmap_bld(aggr_type,theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod
use mld_inner_mod, mld_protect_name => mld_saggrmap_bld
@@ -83,7 +83,8 @@ subroutine mld_saggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
! Arguments
integer, intent(in) :: aggr_type
type(psb_sspmat_type), intent(in), target :: a
real(psb_spk_), intent(in) :: theta
type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
@@ -91,8 +92,7 @@ subroutine mld_saggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
! Local variables
integer, allocatable :: ils(:), neigh(:)
integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m
type(psb_sspmat_type), target :: atmp, atrans
type(psb_sspmat_type), pointer :: apnt
type(psb_sspmat_type) :: atmp, atrans
logical :: recovery
integer :: debug_level, debug_unit
integer :: ictxt,np,me,err_act
@@ -116,8 +116,93 @@ subroutine mld_saggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
ncol = psb_cd_get_local_cols(desc_a)
select case (aggr_type)
case (mld_dec_aggr_,mld_sym_dec_aggr_)
case (mld_dec_aggr_)
call mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info)
case (mld_sym_dec_aggr_)
nr = psb_sp_get_nrows(a)
call psb_sp_clip(a,atmp,info,imax=nr,jmax=nr,&
& rscale=.false.,cscale=.false.)
atmp%m=nr
atmp%k=nr
if (info == 0) call psb_transp(atmp,atrans,fmt='COO')
if (info == 0) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.)
atmp%m=nr
atmp%k=nr
if (info == 0) call psb_sp_free(atrans,info)
if (info == 0) call psb_ipcoo2csr(atmp,info)
if (info == 0) call mld_dec_map_bld(theta,atmp,desc_a,nlaggr,ilaggr,info)
if (info == 0) call psb_sp_free(atmp,info)
case default
info = -1
call psb_errpush(30,name,i_err=(/1,aggr_type,0,0,0/))
goto 9999
end select
if (info/=0) then
info=4001
call psb_errpush(info,name,a_err='dec_map_bld')
goto 9999
end if
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
contains
subroutine mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod
use mld_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=0
name = 'mld_dec_map_bld'
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!
! Note. At the time being we are ignoring aggr_type so
! that we only have decoupled aggregation. This might
! change in the future.
!
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%m
allocate(ilaggr(nr),neigh(nr),stat=info)
if(info /= 0) then
@@ -127,36 +212,28 @@ subroutine mld_saggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
goto 9999
end if
allocate(diag(nr),stat=info)
if(info /= 0) then
info=4025
call psb_errpush(info,name,i_err=(/nr,0,0,0,0/),&
& a_err='real(psb_spk_)')
goto 9999
end if
call psb_sp_getdiag(a,diag,info)
if(info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_sp_getdiag')
goto 9999
end if
do i=1, nr
ilaggr(i) = -(nr+1)
end do
if (aggr_type == mld_dec_aggr_) then
apnt => a
else
call psb_sp_clip(a,atmp,info,imax=nr,jmax=nr,&
& rscale=.false.,cscale=.false.)
atmp%m=nr
atmp%k=nr
if (info == 0) call psb_transp(atmp,atrans,fmt='COO')
if (info == 0) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.)
atmp%m=nr
atmp%k=nr
if (info == 0) call psb_sp_free(atrans,info)
if (info == 0) call psb_ipcoo2csr(atmp,info)
apnt => atmp
if (info/=0) then
info=4001
call psb_errpush(info,name,a_err='init apnt')
goto 9999
end if
end if
! 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.
@@ -171,26 +248,32 @@ subroutine mld_saggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
! 1. Untouched nodes are marked >0 together
! with their neighbours
!
icnt = icnt + 1
naggr = naggr + 1
icnt = icnt + 1
naggr = naggr + 1
ilaggr(i) = naggr
call psb_neigh(apnt,i,neigh,n_ne,info,lev=1)
call psb_sp_getrow(i,a,nz,irow,icol,val,info)
if (info/=0) then
info=4010
call psb_errpush(info,name,a_err='psb_neigh')
call psb_errpush(info,name,a_err='psb_sp_getrow')
goto 9999
end if
do k=1, n_ne
j = neigh(k)
if ((1<=j).and.(j<=nr)) then
ilaggr(j) = naggr
endif
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 psb_neigh(apnt,i,neigh,n_ne,info,lev=2)
call psb_neigh(a,i,neigh,n_ne,info,lev=2)
if (info/=0) then
info=4010
call psb_errpush(info,name,a_err='psb_neigh')
@@ -217,7 +300,7 @@ subroutine mld_saggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
!
! Phase two: sweep over leftovers.
!
allocate(ils(naggr+10),stat=info)
allocate(ils(nr),stat=info)
if(info /= 0) then
info=4025
call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),&
@@ -253,57 +336,55 @@ subroutine mld_saggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
if (ilaggr(i) < 0) then
!
! Now some silly rule to break ties:
! Group with smallest adjacent aggregate.
! Group with adjacent aggregate.
!
isz = nr+1
ia = -1
call psb_neigh(apnt,i,neigh,n_ne,info,lev=1)
isz = nr+1
ia = -1
cpling = 0.0
call psb_sp_getrow(i,a,nz,irow,icol,val,info)
if (info/=0) then
info=4010
call psb_errpush(info,name,a_err='psb_neigh')
call psb_errpush(info,name,a_err='psb_sp_getrow')
goto 9999
end if
do j=1, n_ne
k = neigh(j)
if ((1<=k).and.(k<=nr)) then
n = ilaggr(k)
if (n>0) then
if (n>naggr) then
info=4001
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?')
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=4001
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?')
goto 9999
end if
if (ils(n) < isz) then
ia = n
isz = ils(n)
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
endif
end if
enddo
if (ia == -1) then
if (ilaggr(i) > -(nr+1)) then
ilaggr(i) = abs(ilaggr(i))
if (ilaggr(I)>naggr) then
info=4001
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 3 ?')
goto 9999
end if
ils(ilaggr(i)) = ils(ilaggr(i)) + 1
!
! This might happen if the pattern is non symmetric.
! Need a better handling.
!
recovery = .true.
else
info=4001
call psb_errpush(info,name,a_err='Unrecoverable error !!')
goto 9999
endif
! 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=4001
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 4? ')
@@ -311,8 +392,9 @@ subroutine mld_saggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
end if
ils(ia) = ils(ia) + 1
endif
end if
enddo
end do
if (debug_level >= psb_debug_outer_) then
if (recovery) then
write(debug_unit,*) me,' ',trim(name),&
@@ -344,6 +426,12 @@ subroutine mld_saggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
call psb_errpush(info,name)
goto 9999
end if
if (naggr > ncol) then
write(0,*) name,'Error : naggr > ncol',naggr,ncol
info=4001
call psb_errpush(info,name,a_err='Fatal error: naggr>ncol')
goto 9999
end if
call psb_realloc(ncol,ilaggr,info)
if (info/=0) then
@@ -365,26 +453,17 @@ subroutine mld_saggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
nlaggr(me+1) = naggr
call psb_sum(ictxt,nlaggr(1:np))
if (aggr_type == mld_sym_dec_aggr_) then
call psb_sp_free(atmp,info)
end if
case default
info = -1
call psb_errpush(30,name,i_err=(/1,aggr_type,0,0,0/))
goto 9999
end select
call psb_erractionrestore(err_act)
return
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()
return
end if
return
end if
return
end subroutine mld_dec_map_bld
end subroutine mld_saggrmap_bld
mlprec/mld_smlprec_bld.f90
+3 -1
View File
@@ -109,6 +109,7 @@ subroutine mld_smlprec_bld(a,desc_a,p,info)
call mld_check_def(p%rprcparm(mld_fact_thrs_),'Eps',szero,is_legal_s_fact_thrs)
end select
call mld_check_def(p%rprcparm(mld_aggr_damp_),'Omega',szero,is_legal_s_omega)
call mld_check_def(p%rprcparm(mld_aggr_thresh_),'Aggr_Thresh',szero,is_legal_s_aggr_thrs)
call mld_check_def(p%iprcparm(mld_smooth_sweeps_),'Jacobi sweeps',&
& 1,is_legal_jac_sweeps)
@@ -118,7 +119,8 @@ subroutine mld_smlprec_bld(a,desc_a,p,info)
! aggregation algorithm. This also defines a tentative prolongator from
! the coarse to the fine level.
!
call mld_aggrmap_bld(p%iprcparm(mld_aggr_alg_),a,desc_a,p%nlaggr,p%mlia,info)
call mld_aggrmap_bld(p%iprcparm(mld_aggr_alg_),p%rprcparm(mld_aggr_thresh_),&
& a,desc_a,p%nlaggr,p%mlia,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='mld_aggrmap_bld')
goto 9999
mlprec/mld_sprecinit.f90
+6 -3
View File
@@ -189,7 +189,8 @@ subroutine mld_sprecinit(p,ptype,info,nlev)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%rprcparm(:) = szero
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_as_
p%baseprecv(ilev_)%iprcparm(mld_sub_solve_) = mld_ilu_n_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_halo_
@@ -204,7 +205,8 @@ subroutine mld_sprecinit(p,ptype,info,nlev)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%rprcparm(:) = szero
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_bjac_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_
@@ -225,7 +227,8 @@ subroutine mld_sprecinit(p,ptype,info,nlev)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%rprcparm(:) = szero
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_bjac_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_
mlprec/mld_sprecset.f90
+12 -3
View File
@@ -595,7 +595,7 @@ subroutine mld_sprecsetr(p,what,val,info,ilev)
else if (ilev_ > 1) then
select case(what)
case(mld_aggr_damp_,mld_fact_thrs_)
case(mld_aggr_damp_,mld_aggr_thresh_,mld_fact_thrs_)
p%baseprecv(ilev_)%rprcparm(what) = val
case default
write(0,*) name,': Error: invalid WHAT'
@@ -610,7 +610,7 @@ subroutine mld_sprecsetr(p,what,val,info,ilev)
select case(what)
case(mld_fact_thrs_)
do ilev_=1,nlev_-1
do ilev_=1,nlev_
if (.not.allocated(p%baseprecv(ilev_)%rprcparm)) then
write(0,*) name,': Error: Uninitialized preconditioner component, should call MLD_PRECINIT'
info = -1
@@ -619,7 +619,16 @@ subroutine mld_sprecsetr(p,what,val,info,ilev)
p%baseprecv(ilev_)%rprcparm(what) = val
end do
case(mld_aggr_damp_)
do ilev_=2,nlev_-1
do ilev_=2,nlev_
if (.not.allocated(p%baseprecv(ilev_)%rprcparm)) then
write(0,*) name,': Error: Uninitialized preconditioner component, should call MLD_PRECINIT'
info = -1
return
endif
p%baseprecv(ilev_)%rprcparm(what) = val
end do
case(mld_aggr_thresh_)
do ilev_=2,nlev_
if (.not.allocated(p%baseprecv(ilev_)%rprcparm)) then
write(0,*) name,': Error: Uninitialized preconditioner component, should call MLD_PRECINIT'
info = -1
mlprec/mld_zaggrmap_bld.f90
+175 -97
View File
@@ -74,7 +74,7 @@
! info - integer, output.
! Error code.
!
subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
subroutine mld_zaggrmap_bld(aggr_type,theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod
use mld_inner_mod, mld_protect_name => mld_zaggrmap_bld
@@ -83,7 +83,8 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
! Arguments
integer, intent(in) :: aggr_type
type(psb_zspmat_type), intent(in), target :: a
real(psb_dpk_), intent(in) :: theta
type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
@@ -91,9 +92,7 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
! Local variables
integer, allocatable :: ils(:), neigh(:)
integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m
type(psb_zspmat_type), target :: atmp, atrans
type(psb_zspmat_type), pointer :: apnt
type(psb_zspmat_type) :: atmp, atrans
logical :: recovery
integer :: debug_level, debug_unit
integer :: ictxt,np,me,err_act
@@ -117,8 +116,93 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
ncol = psb_cd_get_local_cols(desc_a)
select case (aggr_type)
case (mld_dec_aggr_,mld_sym_dec_aggr_)
case (mld_dec_aggr_)
call mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info)
case (mld_sym_dec_aggr_)
nr = psb_sp_get_nrows(a)
call psb_sp_clip(a,atmp,info,imax=nr,jmax=nr,&
& rscale=.false.,cscale=.false.)
atmp%m=nr
atmp%k=nr
if (info == 0) call psb_transp(atmp,atrans,fmt='COO')
if (info == 0) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.)
atmp%m=nr
atmp%k=nr
if (info == 0) call psb_sp_free(atrans,info)
if (info == 0) call psb_ipcoo2csr(atmp,info)
if (info == 0) call mld_dec_map_bld(theta,atmp,desc_a,nlaggr,ilaggr,info)
if (info == 0) call psb_sp_free(atmp,info)
case default
info = -1
call psb_errpush(30,name,i_err=(/1,aggr_type,0,0,0/))
goto 9999
end select
if (info/=0) then
info=4001
call psb_errpush(info,name,a_err='dec_map_bld')
goto 9999
end if
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
contains
subroutine mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod
use mld_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=0
name = 'mld_dec_map_bld'
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!
! Note. At the time being we are ignoring aggr_type so
! that we only have decoupled aggregation. This might
! change in the future.
!
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%m
allocate(ilaggr(nr),neigh(nr),stat=info)
if(info /= 0) then
@@ -128,36 +212,28 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
goto 9999
end if
allocate(diag(nr),stat=info)
if(info /= 0) then
info=4025
call psb_errpush(info,name,i_err=(/nr,0,0,0,0/),&
& a_err='complex(psb_dpk_)')
goto 9999
end if
call psb_sp_getdiag(a,diag,info)
if(info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_sp_getdiag')
goto 9999
end if
do i=1, nr
ilaggr(i) = -(nr+1)
end do
if (aggr_type == mld_dec_aggr_) then
apnt => a
else
call psb_sp_clip(a,atmp,info,imax=nr,jmax=nr,&
& rscale=.false.,cscale=.false.)
atmp%m=nr
atmp%k=nr
if (info == 0) call psb_transp(atmp,atrans,fmt='COO')
if (info == 0) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.)
atmp%m=nr
atmp%k=nr
if (info == 0) call psb_sp_free(atrans,info)
if (info == 0) call psb_ipcoo2csr(atmp,info)
apnt => atmp
if (info/=0) then
info=4001
call psb_errpush(info,name,a_err='init apnt')
goto 9999
end if
end if
! 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.
@@ -172,26 +248,32 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
! 1. Untouched nodes are marked >0 together
! with their neighbours
!
icnt = icnt + 1
naggr = naggr + 1
icnt = icnt + 1
naggr = naggr + 1
ilaggr(i) = naggr
call psb_neigh(apnt,i,neigh,n_ne,info,lev=1)
call psb_sp_getrow(i,a,nz,irow,icol,val,info)
if (info/=0) then
info=4010
call psb_errpush(info,name,a_err='psb_neigh')
call psb_errpush(info,name,a_err='psb_sp_getrow')
goto 9999
end if
do k=1, n_ne
j = neigh(k)
if ((1<=j).and.(j<=nr)) then
ilaggr(j) = naggr
endif
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 psb_neigh(apnt,i,neigh,n_ne,info,lev=2)
call psb_neigh(a,i,neigh,n_ne,info,lev=2)
if (info/=0) then
info=4010
call psb_errpush(info,name,a_err='psb_neigh')
@@ -218,7 +300,7 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
!
! Phase two: sweep over leftovers.
!
allocate(ils(naggr+10),stat=info)
allocate(ils(nr),stat=info)
if(info /= 0) then
info=4025
call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),&
@@ -254,57 +336,55 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
if (ilaggr(i) < 0) then
!
! Now some silly rule to break ties:
! Group with smallest adjacent aggregate.
! Group with adjacent aggregate.
!
isz = nr+1
ia = -1
call psb_neigh(apnt,i,neigh,n_ne,info,lev=1)
isz = nr+1
ia = -1
cpling = 0.0d0
call psb_sp_getrow(i,a,nz,irow,icol,val,info)
if (info/=0) then
info=4010
call psb_errpush(info,name,a_err='psb_neigh')
call psb_errpush(info,name,a_err='psb_sp_getrow')
goto 9999
end if
do j=1, n_ne
k = neigh(j)
if ((1<=k).and.(k<=nr)) then
n = ilaggr(k)
if (n>0) then
if (n>naggr) then
info=4001
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?')
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=4001
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?')
goto 9999
end if
if (ils(n) < isz) then
ia = n
isz = ils(n)
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
endif
end if
enddo
if (ia == -1) then
if (ilaggr(i) > -(nr+1)) then
ilaggr(i) = abs(ilaggr(i))
if (ilaggr(I)>naggr) then
info=4001
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 3 ?')
goto 9999
end if
ils(ilaggr(i)) = ils(ilaggr(i)) + 1
!
! This might happen if the pattern is non symmetric.
! Need a better handling.
!
recovery = .true.
else
info=4001
call psb_errpush(info,name,a_err='Unrecoverable error !!')
goto 9999
endif
! 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=4001
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 4? ')
@@ -312,8 +392,9 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
end if
ils(ia) = ils(ia) + 1
endif
end if
enddo
end do
if (debug_level >= psb_debug_outer_) then
if (recovery) then
write(debug_unit,*) me,' ',trim(name),&
@@ -345,6 +426,12 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
call psb_errpush(info,name)
goto 9999
end if
if (naggr > ncol) then
write(0,*) name,'Error : naggr > ncol',naggr,ncol
info=4001
call psb_errpush(info,name,a_err='Fatal error: naggr>ncol')
goto 9999
end if
call psb_realloc(ncol,ilaggr,info)
if (info/=0) then
@@ -366,26 +453,17 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
nlaggr(me+1) = naggr
call psb_sum(ictxt,nlaggr(1:np))
if (aggr_type == mld_sym_dec_aggr_) then
call psb_sp_free(atmp,info)
end if
case default
info = -1
call psb_errpush(30,name,i_err=(/1,aggr_type,0,0,0/))
goto 9999
end select
call psb_erractionrestore(err_act)
return
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()
return
end if
return
end if
return
end subroutine mld_dec_map_bld
end subroutine mld_zaggrmap_bld
mlprec/mld_zmlprec_bld.f90
+3 -1
View File
@@ -109,6 +109,7 @@ subroutine mld_zmlprec_bld(a,desc_a,p,info)
call mld_check_def(p%rprcparm(mld_fact_thrs_),'Eps',dzero,is_legal_fact_thrs)
end select
call mld_check_def(p%rprcparm(mld_aggr_damp_),'Omega',dzero,is_legal_omega)
call mld_check_def(p%rprcparm(mld_aggr_thresh_),'Aggr_Thresh',dzero,is_legal_aggr_thrs)
call mld_check_def(p%iprcparm(mld_smooth_sweeps_),'Jacobi sweeps',&
& 1,is_legal_jac_sweeps)
@@ -118,7 +119,8 @@ subroutine mld_zmlprec_bld(a,desc_a,p,info)
! aggregation algorithm. This also defines a tentative prolongator from
! the coarse to the fine level.
!
call mld_aggrmap_bld(p%iprcparm(mld_aggr_alg_),a,desc_a,p%nlaggr,p%mlia,info)
call mld_aggrmap_bld(p%iprcparm(mld_aggr_alg_),p%rprcparm(mld_aggr_thresh_),&
& a,desc_a,p%nlaggr,p%mlia,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='mld_aggrmap_bld')
goto 9999
mlprec/mld_zprecinit.f90
+6 -3
View File
@@ -189,7 +189,8 @@ subroutine mld_zprecinit(p,ptype,info,nlev)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%rprcparm(:) = dzero
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_as_
p%baseprecv(ilev_)%iprcparm(mld_sub_solve_) = mld_ilu_n_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_halo_
@@ -204,7 +205,8 @@ subroutine mld_zprecinit(p,ptype,info,nlev)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%rprcparm(:) = dzero
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_bjac_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_
@@ -225,7 +227,8 @@ subroutine mld_zprecinit(p,ptype,info,nlev)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%rprcparm(:) = dzero
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_bjac_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_
mlprec/mld_zprecset.f90
+12 -3
View File
@@ -595,7 +595,7 @@ subroutine mld_zprecsetr(p,what,val,info,ilev)
else if (ilev_ > 1) then
select case(what)
case(mld_aggr_damp_,mld_fact_thrs_)
case(mld_aggr_damp_,mld_aggr_thresh_,mld_fact_thrs_)
p%baseprecv(ilev_)%rprcparm(what) = val
case default
write(0,*) name,': Error: invalid WHAT'
@@ -610,7 +610,7 @@ subroutine mld_zprecsetr(p,what,val,info,ilev)
select case(what)
case(mld_fact_thrs_)
do ilev_=1,nlev_-1
do ilev_=1,nlev_
if (.not.allocated(p%baseprecv(ilev_)%rprcparm)) then
write(0,*) name,': Error: Uninitialized preconditioner component, should call MLD_PRECINIT'
info = -1
@@ -619,7 +619,16 @@ subroutine mld_zprecsetr(p,what,val,info,ilev)
p%baseprecv(ilev_)%rprcparm(what) = val
end do
case(mld_aggr_damp_)
do ilev_=2,nlev_-1
do ilev_=2,nlev_
if (.not.allocated(p%baseprecv(ilev_)%rprcparm)) then
write(0,*) name,': Error: Uninitialized preconditioner component, should call MLD_PRECINIT'
info = -1
return
endif
p%baseprecv(ilev_)%rprcparm(what) = val
end do
case(mld_aggr_thresh_)
do ilev_=2,nlev_
if (.not.allocated(p%baseprecv(ilev_)%rprcparm)) then
write(0,*) name,': Error: Uninitialized preconditioner component, should call MLD_PRECINIT'
info = -1
test/fileread/cf_sample.f90
+4 -1
View File
@@ -67,6 +67,7 @@ program cf_sample
real(psb_spk_) :: cthres ! Threshold for fact. 1 ILU(T)
integer :: cjswp ! Jacobi sweeps
real(psb_spk_) :: omega ! smoother omega
real(psb_spk_) :: athres ! smoother aggregation threshold
end type precdata
type(precdata) :: prec_choice
@@ -242,7 +243,7 @@ program cf_sample
call mld_precset(prec,mld_coarse_solve_,prec_choice%csolve,info)
call mld_precset(prec,mld_sub_fill_in_,prec_choice%cfill,info,ilev=nlv)
call mld_precset(prec,mld_fact_thrs_,prec_choice%cthres,info,ilev=nlv)
call mld_precset(prec,mld_smooth_sweeps_,prec_choice%cjswp,info,ilev=nlv)
call mld_precset(prec,mld_aggr_thresh_,prec_choice%athres,info)
call mld_precset(prec,mld_smooth_sweeps_,prec_choice%cjswp,info,ilev=nlv)
if (prec_choice%omega>=0.0) then
call mld_precset(prec,mld_aggr_damp_,prec_choice%omega,info,ilev=nlv)
@@ -392,6 +393,7 @@ contains
call read_data(prec%cthres,5) ! Threshold for fact. 1 ILU(T)
call read_data(prec%cjswp,5) ! Jacobi sweeps
call read_data(prec%omega,5) ! smoother omega
call read_data(prec%athres,5) ! smoother aggr thresh
end if
end if
@@ -426,6 +428,7 @@ contains
call psb_bcast(icontxt,prec%cthres) ! Threshold for fact. 1 ILU(T)
call psb_bcast(icontxt,prec%cjswp) ! Jacobi sweeps
call psb_bcast(icontxt,prec%omega) ! smoother omega
call psb_bcast(icontxt,prec%athres) ! smoother aggr thresh
end if
end subroutine get_parms
test/fileread/df_bench.f90
+1
View File
@@ -29,6 +29,7 @@ program df_bench
real(psb_dpk_) :: thr2 ! Threshold for fact. 1 ILU(T)
integer :: jswp ! Jacobi sweeps
real(psb_dpk_) :: omega ! smoother omega
real(psb_dpk_) :: athres ! smoother aggregation
character(len=40) :: descr ! verbose description of the prec
end type precdata
type(precdata), allocatable :: precs(:)
test/fileread/df_sample.f90
+4 -1
View File
@@ -67,6 +67,7 @@ program df_sample
real(psb_dpk_) :: cthres ! Threshold for fact. 1 ILU(T)
integer :: cjswp ! Jacobi sweeps
real(psb_dpk_) :: omega ! smoother omega
real(psb_dpk_) :: athres ! smoother aggregation threshold
end type precdata
type(precdata) :: prec_choice
@@ -242,7 +243,7 @@ program df_sample
call mld_precset(prec,mld_coarse_solve_,prec_choice%csolve,info)
call mld_precset(prec,mld_sub_fill_in_,prec_choice%cfill,info,ilev=nlv)
call mld_precset(prec,mld_fact_thrs_,prec_choice%cthres,info,ilev=nlv)
call mld_precset(prec,mld_smooth_sweeps_,prec_choice%cjswp,info,ilev=nlv)
call mld_precset(prec,mld_aggr_thresh_,prec_choice%athres,info)
call mld_precset(prec,mld_smooth_sweeps_,prec_choice%cjswp,info,ilev=nlv)
if (prec_choice%omega>=0.0) then
call mld_precset(prec,mld_aggr_damp_,prec_choice%omega,info,ilev=nlv)
@@ -392,6 +393,7 @@ contains
call read_data(prec%cthres,5) ! Threshold for fact. 1 ILU(T)
call read_data(prec%cjswp,5) ! Jacobi sweeps
call read_data(prec%omega,5) ! smoother omega
call read_data(prec%athres,5) ! smoother aggr thresh
end if
end if
@@ -426,6 +428,7 @@ contains
call psb_bcast(icontxt,prec%cthres) ! Threshold for fact. 1 ILU(T)
call psb_bcast(icontxt,prec%cjswp) ! Jacobi sweeps
call psb_bcast(icontxt,prec%omega) ! smoother omega
call psb_bcast(icontxt,prec%athres) ! smoother aggr thresh
end if
end subroutine get_parms
test/fileread/runs/cfs.inp
+1
View File
@@ -27,3 +27,4 @@ ILU ! Coarse level: solver ILU ILUT UMF SLU SLUDIST
1.d-4 ! Coarse level: Threshold T for ILU(T,P)
4 ! Coarse level: Number of Jacobi sweeps
-1.0d0 ! Smoother Omega: if < 0 means library choice.
0.01d0 ! Smoother Aggregation Threshold: >= 0.0
test/fileread/runs/dfs.inp
+3 -2
View File
@@ -1,5 +1,5 @@
thm1000x600.mtx ! les_t4.mtx ! young1r.mtx This (and others) from: http://math.nist.gov/MatrixMarket/ or
NONE !les_t4.rhs ! rhs.mtx http://www.cise.ufl.edu/research/sparse/matrices/index.html
matphi_140x33x45.mtx !A_1M_gps.mtx !thm1000x600.mtx ! les_t4.mtx ! young1r.mtx This (and others) from: http://math.nist.gov/MatrixMarket/ or
tnoto_phi.mtx !NONE !les_t4.rhs ! rhs.mtx http://www.cise.ufl.edu/research/sparse/matrices/index.html
BICGSTAB ! Iterative method: BiCGSTAB BiCG CGS RGMRES BiCGSTABL CG
CSR ! Storage format CSR COO JAD
0 ! IPART: Partition method 0: BLK 2: graph (with Metis)
@@ -27,3 +27,4 @@ ILU ! Coarse level: solver ILU ILUT UMF SLU SLUDIST
1.d-4 ! Coarse level: Threshold T for ILU(T,P)
4 ! Coarse level: Number of Jacobi sweeps
-1.0d0 ! Smoother Omega: if < 0 means library choice.
0.01d0 ! Smoother Aggregation Threshold: >= 0.0
test/fileread/runs/drt.sh
+1 -1
View File
@@ -16,7 +16,7 @@ do
# 3rd batch: 14bis,64bis 4 sweeps
echo "mpirun -np $np -machinefile locm df_bench >>log.part$part.ren$renum.${np}p"
#/usr/local/mpich-gcc42/bin/mpirun -np $np -machinefile locm df_bench <<EOF
/usr/local/mpich-gcc42/bin/mpirun -np $np -machinefile locm df_bench >>run.gfc.kiva.log.${np}p.$date 2>err.gfc.kiva.log.${np}p.$date <<EOF
mpirun -np $np -machinefile locm df_bench >>run.gfc.kiva.log.${np}p.$date 2>err.gfc.kiva.log.${np}p.$date <<EOF
out.${np}p.$date Out file 1: summary
stat.${np}p.$date Out file 2: detailed for statistics
BICGSTAB iterative method to use
test/fileread/runs/sfs.inp
+1
View File
@@ -27,3 +27,4 @@ ILU ! Coarse level: solver ILU ILUT UMF SLU SLUDIST
1.d-4 ! Coarse level: Threshold T for ILU(T,P)
4 ! Coarse level: Number of Jacobi sweeps
-1.0d0 ! Smoother Omega: if < 0 means library choice.
0.01d0 ! Smoother Aggregation Threshold: >= 0.0
test/fileread/runs/zfs.inp
+1
View File
@@ -27,3 +27,4 @@ ILU ! Coarse level: solver ILU ILUT UMF SLU SLUDIST
1.d-4 ! Coarse level: Threshold T for ILU(T,P)
4 ! Coarse level: Number of Jacobi sweeps
-1.0d0 ! Smoother Omega: if < 0 means library choice.
0.01d0 ! Smoother Aggregation Threshold: >= 0.0
test/fileread/sf_sample.f90
+4 -1
View File
@@ -67,6 +67,7 @@ program sf_sample
real(psb_spk_) :: cthres ! Threshold for fact. 1 ILU(T)
integer :: cjswp ! Jacobi sweeps
real(psb_spk_) :: omega ! smoother omega
real(psb_spk_) :: athres ! smoother aggregation threshold
end type precdata
type(precdata) :: prec_choice
@@ -242,7 +243,7 @@ program sf_sample
call mld_precset(prec,mld_coarse_solve_,prec_choice%csolve,info)
call mld_precset(prec,mld_sub_fill_in_,prec_choice%cfill,info,ilev=nlv)
call mld_precset(prec,mld_fact_thrs_,prec_choice%cthres,info,ilev=nlv)
call mld_precset(prec,mld_smooth_sweeps_,prec_choice%cjswp,info,ilev=nlv)
call mld_precset(prec,mld_aggr_thresh_,prec_choice%athres,info)
call mld_precset(prec,mld_smooth_sweeps_,prec_choice%cjswp,info,ilev=nlv)
if (prec_choice%omega>=0.0) then
call mld_precset(prec,mld_aggr_damp_,prec_choice%omega,info,ilev=nlv)
@@ -392,6 +393,7 @@ contains
call read_data(prec%cthres,5) ! Threshold for fact. 1 ILU(T)
call read_data(prec%cjswp,5) ! Jacobi sweeps
call read_data(prec%omega,5) ! smoother omega
call read_data(prec%athres,5) ! smoother aggr thresh
end if
end if
@@ -426,6 +428,7 @@ contains
call psb_bcast(icontxt,prec%cthres) ! Threshold for fact. 1 ILU(T)
call psb_bcast(icontxt,prec%cjswp) ! Jacobi sweeps
call psb_bcast(icontxt,prec%omega) ! smoother omega
call psb_bcast(icontxt,prec%athres) ! smoother aggr thresh
end if
end subroutine get_parms
test/fileread/zf_sample.f90
+4 -1
View File
@@ -67,6 +67,7 @@ program zf_sample
real(psb_dpk_) :: cthres ! Threshold for fact. 1 ILU(T)
integer :: cjswp ! Jacobi sweeps
real(psb_dpk_) :: omega ! smoother omega
real(psb_dpk_) :: athres ! smoother aggregation threshold
end type precdata
type(precdata) :: prec_choice
@@ -242,7 +243,7 @@ program zf_sample
call mld_precset(prec,mld_coarse_solve_,prec_choice%csolve,info)
call mld_precset(prec,mld_sub_fill_in_,prec_choice%cfill,info,ilev=nlv)
call mld_precset(prec,mld_fact_thrs_,prec_choice%cthres,info,ilev=nlv)
call mld_precset(prec,mld_smooth_sweeps_,prec_choice%cjswp,info,ilev=nlv)
call mld_precset(prec,mld_aggr_thresh_,prec_choice%athres,info)
call mld_precset(prec,mld_smooth_sweeps_,prec_choice%cjswp,info,ilev=nlv)
if (prec_choice%omega>=0.0) then
call mld_precset(prec,mld_aggr_damp_,prec_choice%omega,info,ilev=nlv)
@@ -392,6 +393,7 @@ contains
call read_data(prec%cthres,5) ! Threshold for fact. 1 ILU(T)
call read_data(prec%cjswp,5) ! Jacobi sweeps
call read_data(prec%omega,5) ! smoother omega
call read_data(prec%athres,5) ! smoother aggr thresh
end if
end if
@@ -426,6 +428,7 @@ contains
call psb_bcast(icontxt,prec%cthres) ! Threshold for fact. 1 ILU(T)
call psb_bcast(icontxt,prec%cjswp) ! Jacobi sweeps
call psb_bcast(icontxt,prec%omega) ! smoother omega
call psb_bcast(icontxt,prec%athres) ! smoother aggr thresh
end if
end subroutine get_parms