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

First D remapping version. Needs full testing an work on clone

Browse Source
remap-coarse
sfilippone 6 days ago
parent b7ac77da69
commit 42313bb843
24 changed files with 444 additions and 311 deletions
Show Stats Download Patch File Download Diff File

27
amgprec/amg_d_onelev_mod.f90
Unescape Escape View File

@ -183,6 +183,7 @@ module amg_d_onelev_mod
integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:)
contains
procedure, pass(rmp) :: clone => d_remap_data_clone
procedure, pass(rmp) :: move_alloc => d_remap_move_alloc
end type amg_d_remap_data_type
type amg_d_onelev_type
@ -706,6 +707,7 @@ contains
if (info == psb_success_) call psb_move_alloc(lv%tprol,b%tprol,info)
if (info == psb_success_) call psb_move_alloc(lv%desc_ac,b%desc_ac,info)
if (info == psb_success_) call psb_move_alloc(lv%linmap,b%linmap,info)
if (info == psb_success_) call lv%remap_data%move_alloc(b%remap_data,info)
b%base_a => lv%base_a
b%base_desc => lv%base_desc
@ -760,6 +762,7 @@ contains
info = psb_success_
nwv = lv%get_wrksz()
if (.not.allocated(lv%wrk)) allocate(lv%wrk,stat=info)
write(0,*) 'From allocate_wrk :',lv%remap_data%desc_ac_pre_remap%is_asb()
if (info == 0) then
if (lv%remap_data%desc_ac_pre_remap%is_asb()) then
!
@ -807,7 +810,8 @@ contains
info = psb_success_
call wk%free(info)
write(0,*) 'wrk_alloc D: "',trim(desc%get_fmt()),'"', present(desc2),desc%is_valid()
write(0,*) 'wrk_alloc D: "',trim(desc%get_fmt()),'"',&
& present(desc2),desc%is_valid()
allocate(wk%wv(nwv),stat=info)
if (present(desc2).and.(desc%is_valid())) then
@ -997,4 +1001,25 @@ contains
call psb_safe_ab_cpy(rmp%nrsrc,remap_out%nrsrc,info)
end subroutine d_remap_data_clone
subroutine d_remap_move_alloc(rmp, remap_out, info)
use psb_base_mod
implicit none
! Arguments
class(amg_d_remap_data_type), target, intent(inout) :: rmp
class(amg_d_remap_data_type), target, intent(inout) :: remap_out
integer(psb_ipk_), intent(out) :: info
!
integer(psb_ipk_) :: i
info = psb_success_
call psb_move_alloc(rmp%ac_pre_remap,remap_out%ac_pre_remap,info)
if (info == psb_success_) &
& call psb_move_alloc(rmp%desc_ac_pre_remap,remap_out%desc_ac_pre_remap,info)
remap_out%idest = rmp%idest
call move_alloc(rmp%isrc,remap_out%isrc)
call move_alloc(rmp%nrsrc,remap_out%nrsrc)
call move_alloc(rmp%naggr,remap_out%naggr)
end subroutine d_remap_move_alloc
end module amg_d_onelev_mod

47
amgprec/amg_d_prec_type.f90
Unescape Escape View File

@ -97,11 +97,11 @@ module amg_d_prec_type
! to keep track against what is put later in the multilevel array
!
integer(psb_ipk_) :: coarse_solver = -1
!
! The multilevel hierarchy
!
type(amg_d_onelev_type), allocatable :: precv(:)
integer(psb_ipk_) :: nlevs
contains
procedure, pass(prec) :: psb_d_apply2_vect => amg_d_apply2_vect
procedure, pass(prec) :: psb_d_apply1_vect => amg_d_apply1_vect
@ -119,6 +119,7 @@ module amg_d_prec_type
procedure, pass(prec) :: cmp_complexity => amg_d_cmp_compl
procedure, pass(prec) :: get_avg_cr => amg_d_get_avg_cr
procedure, pass(prec) :: cmp_avg_cr => amg_d_cmp_avg_cr
procedure, pass(prec) :: set_nlevs => amg_d_set_nlevs
procedure, pass(prec) :: get_nlevs => amg_d_get_nlevs
procedure, pass(prec) :: get_nzeros => amg_d_get_nzeros
procedure, pass(prec) :: sizeof => amg_dprec_sizeof
@ -141,6 +142,7 @@ module amg_d_prec_type
procedure, pass(prec) :: smoothers_free => amg_d_smoothers_free
procedure, pass(prec) :: descr => amg_dfile_prec_descr
procedure, pass(prec) :: memory_use => amg_dfile_prec_memory_use
end type amg_dprec_type
private :: amg_d_dump, amg_d_get_compl, amg_d_cmp_compl,&
@ -438,7 +440,17 @@ contains
if (allocated(prec%precv)) then
val = size(prec%precv)
end if
val = prec%nlevs
write(0,*) ' NLEVS: ',prec%nlevs, val,size(prec%precv)
end function amg_d_get_nlevs
subroutine amg_d_set_nlevs(prec,nl)
implicit none
class(amg_dprec_type), intent(inout) :: prec
integer(psb_ipk_) :: nl
prec%nlevs = nl
end subroutine amg_d_set_nlevs
!
! Function returning the size of the amg_prec_type data structure
! in bytes or in number of nonzeros of the operator(s) involved.
@ -509,7 +521,7 @@ contains
real(psb_dpk_) :: num, den, nmin
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: il
integer(psb_ipk_) :: il,nl
num = -done
den = done
@ -519,7 +531,10 @@ contains
num = prec%precv(il)%base_a%get_nzeros()
if (num >= dzero) then
den = num
do il=2,size(prec%precv)
nl = prec%get_nlevs()
write(0,*) 'Inside cmp_compl ',nl,size(prec%precv)
do il=2, nl
write(0,*) ' ',il,associated(prec%precv(il)%base_a)
num = num + max(0,prec%precv(il)%base_a%get_nzeros())
end do
end if
@ -562,13 +577,11 @@ contains
avgcr = dzero
ctxt = prec%ctxt
call psb_info(ctxt,iam,np)
if (allocated(prec%precv)) then
nl = size(prec%precv)
do il=2,nl
avgcr = avgcr + max(dzero,prec%precv(il)%szratio)
end do
avgcr = avgcr / (nl-1)
end if
nl = prec%get_nlevs()
do il=2,nl
avgcr = avgcr + max(dzero,prec%precv(il)%szratio)
end do
avgcr = avgcr / (nl-1)
call psb_sum(ctxt,avgcr)
prec%ag_data%avg_cr = avgcr/np
end subroutine amg_d_cmp_avg_cr
@ -667,7 +680,7 @@ contains
end if
if (allocated(prec%precv)) then
do i=1,size(prec%precv)
do i=1,prec%get_nlevs()
call prec%precv(i)%free_smoothers(info)
end do
end if
@ -857,7 +870,7 @@ contains
info = 0
ctxt = prec%ctxt
call psb_info(ctxt,iam,np)
iln = size(prec%precv)
iln = prec%get_nlevs()
if (present(istart)) then
il1 = max(1,istart)
else
@ -893,7 +906,7 @@ contains
info = psb_success_
if (allocated(prec%precv)) then
do i=1,size(prec%precv)
do i=1,prec%get_nlevs()
if (info == psb_success_ ) &
& call prec%precv(i)%cnv(info,amold=amold,vmold=vmold,imold=imold)
end do
@ -930,8 +943,9 @@ contains
pout%ctxt = prec%ctxt
pout%ag_data = prec%ag_data
pout%outer_sweeps = prec%outer_sweeps
pout%nlevs = prec%nlevs
if (allocated(prec%precv)) then
ln = size(prec%precv)
ln = prec%get_nlevs()
allocate(pout%precv(ln),stat=info)
if (info /= psb_success_) goto 9999
if (ln >= 1) then
@ -939,6 +953,7 @@ contains
end if
do lev=2, ln
if (info /= psb_success_) exit
write(0,*) 'Inner_clone must be checked and reimplemented! '
call prec%precv(lev)%clone(pout%precv(lev),info)
if (info == psb_success_) then
pout%precv(lev)%base_a => pout%precv(lev)%ac
@ -1018,7 +1033,7 @@ contains
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
nlev = size(prec%precv)
nlev = prec%get_nlevs()
level = 1
do level = 1, nlev
call prec%precv(level)%allocate_wrk(info,vmold=vmold)
@ -1058,7 +1073,7 @@ contains
end if
if (allocated(prec%precv)) then
nlev = size(prec%precv)
nlev = prec%get_nlevs()
do level = 1, nlev
call prec%precv(level)%free_wrk(info)
end do

25
amgprec/impl/amg_cmlprec_aply.f90
Unescape Escape View File

@ -527,8 +527,7 @@ contains
! Apply the restriction
call p%precv(level+1)%map_rstr(cone,vx2l,&
& czero,p%precv(level+1)%wrk%vx2l,&
& info,work=work,&
& vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1))
& info,work=work,vtx=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
@ -547,8 +546,7 @@ contains
!
call p%precv(level+1)%map_prol(cone,&
& p%precv(level+1)%wrk%vy2l, cone,vy2l,&
& info,work=work,&
& vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1))
& info,work=work, vty=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
@ -667,8 +665,7 @@ contains
call p%precv(level+1)%map_rstr(cone,vty,&
& czero,p%precv(level+1)%wrk%vx2l,&
& info,work=work,&
& vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1))
& info,work=work,vtx=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
@ -678,8 +675,7 @@ contains
! Shortcut: just transfer x2l.
call p%precv(level+1)%map_rstr(cone,vx2l,&
& czero,p%precv(level+1)%wrk%vx2l,&
& info,work=work,&
& vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1))
& info,work=work,vtx=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
@ -694,8 +690,7 @@ contains
!
call p%precv(level+1)%map_prol(cone,&
& p%precv(level+1)%wrk%vy2l,cone,vy2l,&
& info,work=work,&
& vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1))
& info,work=work,vty=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
@ -714,7 +709,7 @@ contains
if (info == psb_success_) &
& call p%precv(level+1)%map_rstr(cone,vty,&
& czero,p%precv(level+1)%wrk%vx2l,info,work=work,&
& vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1))
& vtx=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during W-cycle restriction')
@ -725,8 +720,7 @@ contains
if (info == psb_success_) call p%precv(level+1)%map_prol(cone, &
& p%precv(level+1)%wrk%vy2l,cone,vy2l,&
& info,work=work,&
& vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1))
& info,work=work,vty=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -914,7 +908,7 @@ contains
call p%precv(level + 1)%map_rstr(cone,vty,&
& czero,p%precv(level + 1)%wrk%vx2l,&
&info,work=work,&
& vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1))
& vtx=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -949,8 +943,7 @@ contains
!
call p%precv(level+1)%map_prol(cone,&
& p%precv(level+1)%wrk%vy2l,cone,vy2l,&
& info,work=work,&
& vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1))
& info,work=work,vty=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&

8
amgprec/impl/amg_d_extprol_bld.F90
Unescape Escape View File

@ -164,7 +164,7 @@ subroutine amg_d_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold)
mxplevs = p%ag_data%max_levs
mnaggratio = p%ag_data%min_cr_ratio
casize = p%ag_data%min_coarse_size
iszv = size(p%precv)
iszv = p%get_nlevs()
nprolv = size(prolv)
nrestrv = size(restrv)
call psb_bcast(ctxt,iszv)
@ -188,7 +188,7 @@ subroutine amg_d_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold)
call psb_errpush(info,name,a_err='Inconsistent min_cr_ratio')
goto 9999
end if
if (iszv /= size(p%precv)) then
if (iszv /= p%get_nlevs()) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Inconsistent size of precv')
goto 9999
@ -283,7 +283,7 @@ subroutine amg_d_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold)
call p%precv(i)%free(info)
end do
call move_alloc(tprecv,p%precv)
iszv = size(p%precv)
iszv = p%get_nlevs()
end if
!
! Finest level first; remember to fix base_a and base_desc
@ -325,7 +325,7 @@ subroutine amg_d_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold)
goto 9999
endif
iszv = size(p%precv)
iszv = p%get_nlevs()
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&

328
amgprec/impl/amg_d_hierarchy_bld.F90
Unescape Escape View File

@ -82,7 +82,7 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info,cpymat)
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: me,np
integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz,&
& nplevs, mxplevs
& nplevs, mxplevs, level
integer(psb_lpk_) :: iaggsize, casize, mncsize, mncszpp
real(psb_dpk_) :: mnaggratio, sizeratio, athresh, aomega
class(amg_d_base_smoother_type), allocatable :: coarse_sm, med_sm, &
@ -98,6 +98,7 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info,cpymat)
character(len=40) :: ch_err
integer(psb_ipk_), save :: idx_bldtp=-1, idx_matasb=-1
logical, parameter :: do_timings=.false.
logical :: stop_hierarchy_loop
type(psb_ctxt_type) :: lctxt
integer(psb_ipk_) :: lme,lnp
@ -141,7 +142,7 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info,cpymat)
mnaggratio = prec%ag_data%min_cr_ratio
mncsize = prec%ag_data%min_coarse_size
mncszpp = prec%ag_data%min_coarse_size_per_process
iszv = size(prec%precv)
iszv = prec%get_nlevs()
call psb_bcast(ctxt,iszv)
call psb_bcast(ctxt,mncsize)
call psb_bcast(ctxt,mncszpp)
@ -167,7 +168,7 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info,cpymat)
call psb_errpush(info,name,a_err='Inconsistent min_cr_ratio')
goto 9999
end if
if (iszv /= size(prec%precv)) then
if (iszv /= prec%get_nlevs()) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Inconsistent size of precv')
goto 9999
@ -182,6 +183,7 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info,cpymat)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
if (iszv == 1) then
!
! This is OK, since it may be called by the user even if there
@ -229,7 +231,6 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info,cpymat)
casize = mncsize
end if
prec%ag_data%target_coarse_size = casize
nplevs = max(itwo,mxplevs)
!
@ -288,9 +289,9 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info,cpymat)
call prec%precv(i)%free(info)
end do
call move_alloc(tprecv,prec%precv)
iszv = size(prec%precv)
call prec%set_nlevs(nplevs)
iszv = prec%get_nlevs()
end if
!
! Finest level first; create a GEN_BLOCK
! copy of the descriptor.
@ -305,6 +306,7 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info,cpymat)
prec%precv(1)%base_desc => prec%precv(1)%desc_ac
newsz = 0
stop_hierarchy_loop = .false.
array_build_loop: do i=2, iszv
!
! Check on the iprcparm contents: they should be the same
@ -352,62 +354,23 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info,cpymat)
! Save op_prol just in case
!
call op_prol%clone(prec%precv(i)%tprol,info)
!
! Check for early termination of aggregation loop.
!
iaggsize = sum(nlaggr)
sizeratio = iaggsize
if (i==2) then
sizeratio = desc_a%get_global_rows()/sizeratio
!
if (i == 2) then
call amg_hierarchy_bld_newsz(i,iszv,&
& desc_a%get_global_rows(),nlaggr,mnaggratio,sizeratio,newsz)
else
sizeratio = sum(prec%precv(i-1)%linmap%naggr)/sizeratio
call amg_hierarchy_bld_newsz(i,iszv,&
& sum(prec%precv(i-1)%linmap%naggr),nlaggr,mnaggratio,sizeratio,newsz)
end if
prec%precv(i)%szratio = sizeratio
if (iaggsize <= casize) newsz = i
if (i == iszv) newsz = i
if (i>2) then
if (sizeratio < mnaggratio) then
if (sizeratio > 1) then
newsz = i
else
!
! We are not gaining
!
newsz = i-1
end if
end if
!
! TO BE REVIWED
!
write(0,*) me,lme,allocated(prec%precv(i-1)%linmap%naggr)
if (allocated(prec%precv(i-1)%linmap%naggr)) then
write(0,*) me,lme,size(prec%precv(i-1)%linmap%naggr), size(nlaggr)
end if
if (.false.) then
if (lme >=0) then
if (all(nlaggr == prec%precv(i-1)%linmap%naggr)) then
newsz=i-1
if (me == 0) then
write(debug_unit,*) trim(name),&
&': Warning: aggregates from level ',&
& newsz
write(debug_unit,*) trim(name),&
&': to level ',&
& iszv,' coincide.'
write(debug_unit,*) trim(name),&
&': Number of levels actually used :',newsz
write(debug_unit,*)
end if
end if
end if
end if
end if
call psb_bcast(ctxt,newsz)
! Handle reallocation, if needed, and then mat_asb to polish off the
! construction
if (newsz > 0) then
!
! This is awkward, we are saving the aggregation parms, for the sake
@ -441,156 +404,166 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info,cpymat)
& a_err=ch_err)
goto 9999
endif
if (amg_get_do_remap().and.(i>=2)) then
block
type(psb_ctxt_type) :: lctxt
integer(psb_ipk_) :: lme,lnp
lctxt = prec%precv(i)%desc_ac%get_ctxt()
call psb_info(lctxt,lme,lnp)
write(0,*) ' Context on remapping ',lme,lnp
if ((lme >=0).and.(lnp>=2)) then
associate(lv=>prec%precv(i), rmp => prec%precv(i)%remap_data)
call lv%desc_ac%clone(rmp%desc_ac_pre_remap,info)
call lv%ac%clone(rmp%ac_pre_remap,info)
write(0,*) ' Doing remapping ',lnp, lnp/2
call psb_remap(lnp/2,rmp%desc_ac_pre_remap,rmp%ac_pre_remap,&
& rmp%idest,rmp%isrc,rmp%nrsrc,rmp%naggr,lv%desc_ac,lv%ac,info)
!!$ write(0,*) me,' Out of remapping ',rmp%desc_ac_pre_remap%get_fmt(),' ',&
!!$ & lv%desc_ac%get_fmt(),sum(lv%linmap%naggr),sum(rmp%naggr)
write(0,*) 'Assignment ',size(lv%linmap%naggr),size(rmp%naggr)
lv%linmap%naggr(:) = rmp%naggr(:)
lv%linmap%p_desc_V => rmp%desc_ac_pre_remap
lv%base_a => lv%ac
lv%base_desc => lv%desc_ac
end associate
end if
end block
end if
exit array_build_loop
!!$ write(0,*) ' Early exit of array_build_loop',i,iszv,info,&
!!$ & prec%precv(i)%remap_data%desc_ac_pre_remap%is_asb()
!!$ exit array_build_loop
level = newsz
stop_hierarchy_loop = .true.
else
if (do_timings) call psb_tic(idx_matasb)
if (do_timings) call psb_tic(idx_matasb)
if (info == psb_success_) call prec%precv(i)%mat_asb(&
& prec%precv(i-1)%base_a,prec%precv(i-1)%base_desc,&
& ilaggr,nlaggr,op_prol,info)
if (do_timings) call psb_toc(idx_matasb)
level = i
end if
if (info /= psb_success_) then
write(ch_err,'(a,i7)') 'Mat asb fail @ level ',i
call psb_errpush(psb_err_internal_error_,name,&
& a_err=ch_err)
goto 9999
endif
if (i<iszv) call prec%precv(i)%update_aggr(prec%precv(i+1),info)
if (amg_get_do_remap().and.(i>=2)) then
!
! Do we want to remap onto a smaller subset of processes?
!
if (amg_policy_do_remap(level,sum(nlaggr))) then
block
type(psb_ctxt_type) :: lctxt
integer(psb_ipk_) :: lme,lnp
lctxt = prec%precv(i)%desc_ac%get_ctxt()
lctxt = prec%precv(level)%desc_ac%get_ctxt()
call psb_info(lctxt,lme,lnp)
write(0,*) ' Context on remapping ',lme,lnp
if ((lme >=0).and.(lnp>=2)) then
associate(lv=>prec%precv(i), rmp => prec%precv(i)%remap_data)
associate(lv=>prec%precv(level), rmp => prec%precv(level)%remap_data)
call lv%desc_ac%clone(rmp%desc_ac_pre_remap,info)
call lv%ac%clone(rmp%ac_pre_remap,info)
write(0,*) ' Doing remapping ',lnp, lnp/2
write(0,*) 'During first remapping desc_ac:',lv%desc_ac%is_asb(),&
& rmp%desc_ac_pre_remap%is_asb()
write(0,*) ' First Doing remapping ',lnp, lnp/2
call psb_remap(lnp/2,rmp%desc_ac_pre_remap,rmp%ac_pre_remap,&
& rmp%idest,rmp%isrc,rmp%nrsrc,rmp%naggr,lv%desc_ac,lv%ac,info)
!!$ write(0,*) me,' Out of remapping ',rmp%desc_ac_pre_remap%get_fmt(),' ',&
!!$ & lv%desc_ac%get_fmt(),sum(lv%linmap%naggr),sum(rmp%naggr)
write(0,*) 'Assignment ',size(lv%linmap%naggr),size(rmp%naggr)
write(0,*) 'First Assignment ',size(lv%linmap%naggr),size(rmp%naggr)
lv%linmap%naggr(:) = rmp%naggr(:)
lv%linmap%p_desc_V => rmp%desc_ac_pre_remap
lv%base_a => lv%ac
lv%base_desc => lv%desc_ac
block
integer(psb_ipk_) :: meu,npu,mev,npv
type(psb_ctxt_type) :: ct
ct = lv%linmap%p_desc_U%get_ctxt()
call psb_info(ct,meu,npu)
ct = lv%linmap%p_desc_V%get_ctxt()
call psb_info(ct,mev,npv)
write(0,*) 'First Check on out remapping ',i,&
& rmp%desc_ac_pre_remap%is_asb(),&
& ':',meu,npu,mev,npv
end block
end associate
end if
end block
write(0,*) 'Second Check on out remapping ',level,&
& prec%precv(level)%remap_data%desc_ac_pre_remap%is_asb()
end if
write(0,*) ' End of array_build_loop',i,iszv,info
end do array_build_loop
write(0,*) ' Done array_build_loop',iszv,newsz,info,psb_errstatus_fatal()
call psb_barrier(ctxt)
if (newsz > 0) then
!
! We exited early from the build loop, need to fix
! the size.
!
allocate(tprecv(newsz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='prec reallocation')
write(ch_err,'(a,i7)') 'Mat asb fail @ level ',i
call psb_errpush(psb_err_internal_error_,name,&
& a_err=ch_err)
goto 9999
endif
do i=1,newsz
call prec%precv(i)%move_alloc(tprecv(i),info)
if (stop_hierarchy_loop) then
exit array_build_loop
else
if (i<iszv) call prec%precv(i)%update_aggr(prec%precv(i+1),info)
end if
end do array_build_loop
write(0,*) ' Done array_build_loop',iszv,newsz,info,psb_errstatus_fatal()
if (newsz>0) then
do i=2,newsz
write(0,*) me,'Newsz Out of array_build_loop ',i,':',&
& prec%precv(i)%remap_data%desc_ac_pre_remap%is_asb()
end do
do i=newsz+1, iszv
call prec%precv(i)%free(info)
write(0,*) 'Calling set_nlevs ',newsz
call prec%set_nlevs(newsz)
else
do i=2, iszv
write(0,*) me,'Out of array_build_loop ',i,':',&
& prec%precv(i)%remap_data%desc_ac_pre_remap%is_asb()
end do
call move_alloc(tprecv,prec%precv)
! Ignore errors from transfer
info = psb_success_
!
! Restart
iszv = newsz
! Fix the pointers, but the level 1 should
! be treated differently
end if
iszv = prec%get_nlevs()
call psb_barrier(ctxt)
if (.false.) then
if (newsz > 0) then
!
! We exited early from the build loop, need to fix
! the size.
!
allocate(tprecv(newsz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='prec reallocation')
goto 9999
endif
do i=1,newsz
call prec%precv(i)%move_alloc(tprecv(i),info)
end do
do i=newsz+1, iszv
call prec%precv(i)%free(info)
end do
call move_alloc(tprecv,prec%precv)
! Ignore errors from transfer
info = psb_success_
!
! Restart
iszv = newsz
! Fix the pointers, but the level 1 should
! be treated differently
if (.not.associated(prec%precv(1)%base_a,a)) then
prec%precv(1)%base_a => prec%precv(1)%ac
end if
if (.not.associated(prec%precv(1)%base_desc,desc_a)) then
prec%precv(1)%base_desc => prec%precv(1)%desc_ac
end if
do i=2, iszv
prec%precv(i)%base_a => prec%precv(i)%ac
prec%precv(i)%base_desc => prec%precv(i)%desc_ac
! This is needed when the linmap object has been built
! reusing the base_desc descriptor through a pointer.
! With PSBLAS 4 we will have a better solution
if (associated(prec%precv(i)%linmap%p_desc_U)) &
& prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
if (associated(prec%precv(i)%linmap%p_desc_V))&
& prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
end do
end if
else
if (.not.associated(prec%precv(1)%base_a,a)) then
prec%precv(1)%base_a => prec%precv(1)%ac
end if
if (.not.associated(prec%precv(1)%base_desc,desc_a)) then
prec%precv(1)%base_desc => prec%precv(1)%desc_ac
end if
do i=2, iszv
prec%precv(i)%base_a => prec%precv(i)%ac
prec%precv(i)%base_desc => prec%precv(i)%desc_ac
! This is needed when the linmap object has been built
! reusing the base_desc descriptor through a pointer.
! With PSBLAS 4 we will have a better solution
if (associated(prec%precv(i)%linmap%p_desc_U)) &
& prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
if (associated(prec%precv(i)%linmap%p_desc_V))&
& prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
end do
end if
write(0,*) ' Done reallocating precv',iszv,newsz,info,psb_errstatus_fatal()
write(0,*) ' Done reallocating precv',iszv,newsz,info
do i=2, iszv
write(0,*) me,'At end of hierarchy_bld level',i,':',&
& prec%precv(i)%remap_data%desc_ac_pre_remap%is_asb()
end do
call psb_barrier(ctxt)
!write(0,*) 'Should we remap? '
if (.false.) then
if (amg_get_do_remap().and.(np>=4)) then
!!$ write(0,*) 'Going for remapping '
if (.true.) then
associate(lv=>prec%precv(iszv), rmp => prec%precv(iszv)%remap_data)
call lv%desc_ac%clone(rmp%desc_ac_pre_remap,info)
call lv%ac%clone(rmp%ac_pre_remap,info)
if (np >= 8) then
call psb_remap(np/4,rmp%desc_ac_pre_remap,rmp%ac_pre_remap,&
& rmp%idest,rmp%isrc,rmp%nrsrc,rmp%naggr,lv%desc_ac,lv%ac,info)
else
call psb_remap(np/2,rmp%desc_ac_pre_remap,rmp%ac_pre_remap,&
& rmp%idest,rmp%isrc,rmp%nrsrc,rmp%naggr,lv%desc_ac,lv%ac,info)
end if
!!$ write(0,*) me,' Out of remapping ',rmp%desc_ac_pre_remap%get_fmt(),' ',&
!!$ & lv%desc_ac%get_fmt(),sum(lv%linmap%naggr),sum(rmp%naggr)
lv%linmap%naggr(:) = rmp%naggr(:)
lv%linmap%p_desc_V => rmp%desc_ac_pre_remap
lv%base_a => lv%ac
lv%base_desc => lv%desc_ac
end associate
end if
end if
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Internal hierarchy build' )
goto 9999
endif
iszv = size(prec%precv)
iszv = prec%get_nlevs()
write(0,*) 'Going for cmp_complexity ',&
& allocated(prec%precv),iszv,size(prec%precv)
call prec%cmp_complexity()
call prec%cmp_avg_cr()
@ -730,4 +703,45 @@ contains
return
end subroutine restore_smoothers
#endif
function amg_policy_do_remap(level,aggsize) result(res)
logical :: res
integer(psb_ipk_) :: level
integer(psb_lpk_) :: aggsize
res = amg_get_do_remap().and.(level>=2)
!!$ res = .false.
end function amg_policy_do_remap
subroutine amg_hierarchy_bld_newsz(level,iszv,prevsize,&
& nlaggr,mnratio,sizeratio,newsz)
implicit none
integer(psb_ipk_) :: level,iszv,newsz
integer(psb_lpk_) :: nlaggr(:)
integer(psb_lpk_) :: prevsize
real(psb_dpk_) :: mnratio, sizeratio
! ==============================
integer(psb_lpk_) :: iaggsize, casize
newsz = 0
iaggsize = sum(nlaggr)
sizeratio = prevsize
sizeratio = sizeratio/iaggsize
if (iaggsize <= casize) newsz = level
if (level == iszv) newsz = level
if (level>2) then
if (sizeratio < mnaggratio) then
if (sizeratio > 1) then
newsz = level
else
!
! We are not gaining
!
newsz = level-1
end if
end if
end if
end subroutine amg_hierarchy_bld_newsz
end subroutine amg_d_hierarchy_bld

2
amgprec/impl/amg_d_hierarchy_rebld.f90
Unescape Escape View File

@ -122,7 +122,7 @@ subroutine amg_d_hierarchy_rebld(a,desc_a,prec,info)
end if
iszv = size(prec%precv)
iszv = prec%get_nlevs()
do i=2, iszv
call prec%precv(i-1)%base_a%cp_to(acsr)

4
amgprec/impl/amg_d_smoothers_bld.f90
Unescape Escape View File

@ -136,9 +136,9 @@ subroutine amg_d_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
!
! Check to ensure all procs have the same
!
iszv = size(prec%precv)
iszv = prec%get_nlevs()
call psb_bcast(ctxt,iszv)
if (iszv /= size(prec%precv)) then
if (iszv /= prec%get_nlevs()) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Inconsistent size of precv')
goto 9999

6
amgprec/impl/amg_dcprecset.F90
Unescape Escape View File

@ -129,7 +129,7 @@ subroutine amg_dcprecseti(p,what,val,info,ilev,ilmax,pos,idx)
return
endif
nlev_ = size(p%precv)
nlev_ = p%get_nlevs()
if (present(ilev)) then
ilev_ = ilev
@ -376,7 +376,7 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx)
return
endif
nlev_ = size(p%precv)
nlev_ = p%get_nlevs()
if (present(ilev)) then
ilev_ = ilev
@ -1057,7 +1057,7 @@ subroutine amg_dcprecsetr(p,what,val,info,ilev,ilmax,pos,idx)
info = 3111
return
endif
nlev_ = size(p%precv)
nlev_ = p%get_nlevs()
if (present(ilev)) then
ilev_ = ilev

4
amgprec/impl/amg_dfile_prec_descr.f90
Unescape Escape View File

@ -129,7 +129,9 @@ subroutine amg_dfile_prec_descr(prec,info,iout,root, verbosity,prefix)
! ensured by amg_precbld).
!
if (me == root_) then
nlev = size(prec%precv)
nlev = prec%get_nlevs()
write(0,*) 'From file_prec_descr: nlev ',nlev
flush(0)
do ilev = 1, nlev
if (.not.allocated(prec%precv(ilev)%sm)) then
info = 3111

2
amgprec/impl/amg_dfile_prec_memory_use.f90
Unescape Escape View File

@ -135,7 +135,7 @@ subroutine amg_dfile_prec_memory_use(prec,info,iout,root, verbosity,prefix,globa
write(iout_,*)
write(iout_,'(a,1x,a)') trim(prefix_),'Preconditioner memory usage'
end if
nlev = size(prec%precv)
nlev = prec%get_nlevs()
do ilev=1,nlev
call prec%precv(ilev)%memory_use(ilev,nlev,ilmin,info, &
& iout=iout_,verbosity=verbosity_,prefix=trim(prefix_),global=global)

65
amgprec/impl/amg_dmlprec_aply.f90
Unescape Escape View File

@ -244,10 +244,10 @@ subroutine amg_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(p%precv)
& ' Entry ', p%get_nlevs()
trans_ = psb_toupper(trans)
nlev = size(p%precv)
nlev = p%get_nlevs()
do_alloc_wrk = .not.allocated(p%precv(1)%wrk)
@ -382,7 +382,7 @@ contains
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nlev = size(p%precv)
nlev = p%get_nlevs()
if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_ml')
@ -468,7 +468,7 @@ contains
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nlev = size(p%precv)
nlev = p%get_nlevs()
if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_add')
@ -497,7 +497,8 @@ contains
if (allocated(p%precv(level)%sm2a)) then
call psb_geaxpby(done,vx2l,dzero,vy2l,base_desc,info)
sweeps = max(p%precv(level)%parms%sweeps_pre,p%precv(level)%parms%sweeps_post)
sweeps = max(p%precv(level)%parms%sweeps_pre,&
& p%precv(level)%parms%sweeps_post)
do k=1, sweeps
call p%precv(level)%sm%apply(done,&
& vy2l,dzero,vty,&
@ -527,8 +528,7 @@ contains
! Apply the restriction
call p%precv(level+1)%map_rstr(done,vx2l,&
& dzero,p%precv(level+1)%wrk%vx2l,&
& info,work=work,&
& vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1))
& info,work=work,vtx=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
@ -547,8 +547,7 @@ contains
!
call p%precv(level+1)%map_prol(done,&
& p%precv(level+1)%wrk%vy2l, done,vy2l,&
& info,work=work,&
& vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1))
& info,work=work, vty=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
@ -595,7 +594,7 @@ contains
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nlev = size(p%precv)
nlev = p%get_nlevs()
if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_mult')
@ -652,23 +651,23 @@ contains
!
if (pre) then
call psb_geaxpby(done,vx2l,&
& dzero,vty,&
& base_desc,info)
if (info == psb_success_) call psb_spmm(-done,base_a,&
& vy2l,done,vty,&
& base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
call psb_geaxpby(done,vx2l,&
& dzero,vty,&
& base_desc,info)
if (info == psb_success_) call psb_spmm(-done,base_a,&
& vy2l,done,vty,&
& base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
write(0,*) me,' Map_rstr from ', level,' to ',level+1
call p%precv(level+1)%map_rstr(done,vty,&
& dzero,p%precv(level+1)%wrk%vx2l,&
& info,work=work,&
& vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1))
& info,work=work,vtx=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
@ -678,8 +677,7 @@ contains
! Shortcut: just transfer x2l.
call p%precv(level+1)%map_rstr(done,vx2l,&
& dzero,p%precv(level+1)%wrk%vx2l,&
& info,work=work,&
& vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1))
& info,work=work,vtx=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
@ -694,8 +692,7 @@ contains
!
call p%precv(level+1)%map_prol(done,&
& p%precv(level+1)%wrk%vy2l,done,vy2l,&
& info,work=work,&
& vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1))
& info,work=work,vty=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
@ -714,7 +711,7 @@ contains
if (info == psb_success_) &
& call p%precv(level+1)%map_rstr(done,vty,&
& dzero,p%precv(level+1)%wrk%vx2l,info,work=work,&
& vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1))
& vtx=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during W-cycle restriction')
@ -725,8 +722,7 @@ contains
if (info == psb_success_) call p%precv(level+1)%map_prol(done, &
& p%precv(level+1)%wrk%vy2l,done,vy2l,&
& info,work=work,&
& vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1))
& info,work=work,vty=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -837,7 +833,7 @@ contains
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nlev = size(p%precv)
nlev = p%get_nlevs()
if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_add')
@ -914,7 +910,7 @@ contains
call p%precv(level + 1)%map_rstr(done,vty,&
& dzero,p%precv(level + 1)%wrk%vx2l,&
&info,work=work,&
& vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1))
& vtx=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -949,8 +945,7 @@ contains
!
call p%precv(level+1)%map_prol(done,&
& p%precv(level+1)%wrk%vy2l,done,vy2l,&
& info,work=work,&
& vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1))
& info,work=work,vty=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&

10
amgprec/impl/amg_dmlprec_aply_a.f90
Unescape Escape View File

@ -249,11 +249,11 @@ subroutine amg_dmlprec_aply_a(alpha,p,x,beta,y,desc_data,trans,work,info)
if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(p%precv)
& ' Entry ', p%get_nlevs()
trans_ = psb_toupper(trans)
nlev = size(p%precv)
nlev = p%get_nlevs()
allocate(mlwrk(nlev),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
@ -362,7 +362,7 @@ contains
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nlev = size(p%precv)
nlev = p%get_nlevs()
if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_ml')
@ -445,7 +445,7 @@ contains
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nlev = size(p%precv)
nlev = p%get_nlevs()
if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_add')
@ -550,7 +550,7 @@ contains
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nlev = size(p%precv)
nlev = p%get_nlevs()
if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_mult')

2
amgprec/impl/amg_dprecinit.F90
Unescape Escape View File

@ -255,8 +255,8 @@ subroutine amg_dprecinit(ctxt,prec,ptype,info)
write(psb_err_unit,*) name,&
&': Warning: Unknown preconditioner type request "',ptype,'"'
info = psb_err_pivot_too_small_
end select
call prec%set_nlevs(nlev_)
call psb_erractionrestore(err_act)
return

25
amgprec/impl/amg_smlprec_aply.f90
Unescape Escape View File

@ -527,8 +527,7 @@ contains
! Apply the restriction
call p%precv(level+1)%map_rstr(sone,vx2l,&
& szero,p%precv(level+1)%wrk%vx2l,&
& info,work=work,&
& vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1))
& info,work=work,vtx=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
@ -547,8 +546,7 @@ contains
!
call p%precv(level+1)%map_prol(sone,&
& p%precv(level+1)%wrk%vy2l, sone,vy2l,&
& info,work=work,&
& vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1))
& info,work=work, vty=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
@ -667,8 +665,7 @@ contains
call p%precv(level+1)%map_rstr(sone,vty,&
& szero,p%precv(level+1)%wrk%vx2l,&
& info,work=work,&
& vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1))
& info,work=work,vtx=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
@ -678,8 +675,7 @@ contains
! Shortcut: just transfer x2l.
call p%precv(level+1)%map_rstr(sone,vx2l,&
& szero,p%precv(level+1)%wrk%vx2l,&
& info,work=work,&
& vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1))
& info,work=work,vtx=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
@ -694,8 +690,7 @@ contains
!
call p%precv(level+1)%map_prol(sone,&
& p%precv(level+1)%wrk%vy2l,sone,vy2l,&
& info,work=work,&
& vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1))
& info,work=work,vty=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
@ -714,7 +709,7 @@ contains
if (info == psb_success_) &
& call p%precv(level+1)%map_rstr(sone,vty,&
& szero,p%precv(level+1)%wrk%vx2l,info,work=work,&
& vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1))
& vtx=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during W-cycle restriction')
@ -725,8 +720,7 @@ contains
if (info == psb_success_) call p%precv(level+1)%map_prol(sone, &
& p%precv(level+1)%wrk%vy2l,sone,vy2l,&
& info,work=work,&
& vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1))
& info,work=work,vty=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -914,7 +908,7 @@ contains
call p%precv(level + 1)%map_rstr(sone,vty,&
& szero,p%precv(level + 1)%wrk%vx2l,&
&info,work=work,&
& vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1))
& vtx=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -949,8 +943,7 @@ contains
!
call p%precv(level+1)%map_prol(sone,&
& p%precv(level+1)%wrk%vy2l,sone,vy2l,&
& info,work=work,&
& vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1))
& info,work=work,vty=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&

25
amgprec/impl/amg_zmlprec_aply.f90
Unescape Escape View File

@ -527,8 +527,7 @@ contains
! Apply the restriction
call p%precv(level+1)%map_rstr(zone,vx2l,&
& zzero,p%precv(level+1)%wrk%vx2l,&
& info,work=work,&
& vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1))
& info,work=work,vtx=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
@ -547,8 +546,7 @@ contains
!
call p%precv(level+1)%map_prol(zone,&
& p%precv(level+1)%wrk%vy2l, zone,vy2l,&
& info,work=work,&
& vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1))
& info,work=work, vty=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
@ -667,8 +665,7 @@ contains
call p%precv(level+1)%map_rstr(zone,vty,&
& zzero,p%precv(level+1)%wrk%vx2l,&
& info,work=work,&
& vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1))
& info,work=work,vtx=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
@ -678,8 +675,7 @@ contains
! Shortcut: just transfer x2l.
call p%precv(level+1)%map_rstr(zone,vx2l,&
& zzero,p%precv(level+1)%wrk%vx2l,&
& info,work=work,&
& vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1))
& info,work=work,vtx=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
@ -694,8 +690,7 @@ contains
!
call p%precv(level+1)%map_prol(zone,&
& p%precv(level+1)%wrk%vy2l,zone,vy2l,&
& info,work=work,&
& vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1))
& info,work=work,vty=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
@ -714,7 +709,7 @@ contains
if (info == psb_success_) &
& call p%precv(level+1)%map_rstr(zone,vty,&
& zzero,p%precv(level+1)%wrk%vx2l,info,work=work,&
& vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1))
& vtx=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during W-cycle restriction')
@ -725,8 +720,7 @@ contains
if (info == psb_success_) call p%precv(level+1)%map_prol(zone, &
& p%precv(level+1)%wrk%vy2l,zone,vy2l,&
& info,work=work,&
& vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1))
& info,work=work,vty=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -914,7 +908,7 @@ contains
call p%precv(level + 1)%map_rstr(zone,vty,&
& zzero,p%precv(level + 1)%wrk%vx2l,&
&info,work=work,&
& vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1))
& vtx=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -949,8 +943,7 @@ contains
!
call p%precv(level+1)%map_prol(zone,&
& p%precv(level+1)%wrk%vy2l,zone,vy2l,&
& info,work=work,&
& vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1))
& info,work=work,vty=wv(1))
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&

12
amgprec/impl/level/amg_c_base_onelev_map_prol.F90
Unescape Escape View File

@ -46,6 +46,14 @@ subroutine amg_c_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt
integer(psb_ipk_), intent(out) :: info
complex(psb_spk_), optional :: work(:)
type(psb_c_vect_type), optional, target, intent(inout) :: vtx,vty
type(psb_c_vect_type), pointer :: vtx_
!!$ write(0,*) 'New map_rstr',lv%remap_data%ac_pre_remap%is_asb()
if (present(vtx)) then
vtx_ => vtx
else
vtx_ => lv%wrk%wv(1)
end if
!!$ write(0,*) 'New map_prol',lv%remap_data%ac_pre_remap%is_asb()
if (lv%remap_data%ac_pre_remap%is_asb()) then
@ -98,7 +106,7 @@ subroutine amg_c_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt
call psb_rcv(ctxt,tv%v%v(1:nrl),idest)
call tv%set_host()
call lv%linmap%map_V2U(alpha,tv,beta,vect_u,info,&
& work=work,vtx=vtx,vty=vty)
& work=work,vtx=vtx_,vty=vty)
end associate
!!$ write(0,*) me, ' Prolongator with remap done '
!!$ flush(0)
@ -107,7 +115,7 @@ subroutine amg_c_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt
else
! Default transfer
call lv%linmap%map_V2U(alpha,vect_v,beta,vect_u,info,&
& work=work,vtx=vtx,vty=vty)
& work=work,vtx=vtx_,vty=vty)
end if
end subroutine amg_c_base_onelev_map_prol_v

22
amgprec/impl/level/amg_c_base_onelev_map_rstr.F90
Unescape Escape View File

@ -47,8 +47,15 @@ subroutine amg_c_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
integer(psb_ipk_), intent(out) :: info
complex(psb_spk_), optional :: work(:)
type(psb_c_vect_type), optional, target, intent(inout) :: vtx,vty
type(psb_c_vect_type), pointer :: vty_
!!$ write(0,*) 'New map_rstr',lv%remap_data%ac_pre_remap%is_asb()
if (present(vty)) then
vty_ => vty
else
vty_ => lv%wrk%wv(1)
end if
if (lv%remap_data%ac_pre_remap%is_asb()) then
!
! Remap has happened, deal with it
@ -74,11 +81,11 @@ subroutine amg_c_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
nrl = lv%remap_data%desc_ac_pre_remap%get_local_rows()
call psb_geall(tv,lv%remap_data%desc_ac_pre_remap,info)
call psb_geasb(tv,lv%remap_data%desc_ac_pre_remap,info,mold=vect_u%v)
write(0,*) me,' map_rstr calling U2V: ',me,np,rme,rnp,tv%get_nrows()
write(0,*) me,' remap map_rstr calling U2V: ',me,np,rme,rnp,tv%get_nrows()
flush(0)
call psb_barrier(ctxt)
call lv%linmap%map_U2V(alpha,vect_u,beta,tv,info,&
& work=work,vtx=vtx,vty=vty)
& work=work,vtx=vtx,vty=vty_)
call tv%sync()
!rsnd = tv%get_vect()
!call psb_snd(ctxt,rsnd(1:nrl),idest)
@ -103,8 +110,15 @@ subroutine amg_c_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
else
! Default transfer
call lv%linmap%map_U2V(alpha,vect_u,beta,vect_v,info,&
& work=work,vtx=vtx,vty=vty)
block
type(psb_ctxt_type) :: ctxt, rctxt
integer(psb_mpk_) :: me, np
ctxt = lv%linmap%p_desc_U%get_ctxt()
call psb_info(ctxt,me,np)
write(0,*) me,' map_rstr calling U2V: ',me,np
call lv%linmap%map_U2V(alpha,vect_u,beta,vect_v,info,&
& work=work,vtx=vtx,vty=vty_)
end block
end if
end subroutine amg_c_base_onelev_map_rstr_v

13
amgprec/impl/level/amg_d_base_onelev_map_prol.F90
Unescape Escape View File

@ -46,6 +46,15 @@ subroutine amg_d_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt
integer(psb_ipk_), intent(out) :: info
real(psb_dpk_), optional :: work(:)
type(psb_d_vect_type), optional, target, intent(inout) :: vtx,vty
type(psb_d_vect_type), pointer :: vtx_
!!$ write(0,*) 'New map_rstr',lv%remap_data%ac_pre_remap%is_asb()
if (present(vtx)) then
vtx_ => vtx
else
vtx_ => lv%wrk%wv(1)
end if
!!$ write(0,*) 'New map_prol',lv%remap_data%ac_pre_remap%is_asb()
if (lv%remap_data%ac_pre_remap%is_asb()) then
@ -98,7 +107,7 @@ subroutine amg_d_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt
call psb_rcv(ctxt,tv%v%v(1:nrl),idest)
call tv%set_host()
call lv%linmap%map_V2U(alpha,tv,beta,vect_u,info,&
& work=work,vtx=vtx,vty=vty)
& work=work,vtx=vtx_,vty=vty)
end associate
!!$ write(0,*) me, ' Prolongator with remap done '
!!$ flush(0)
@ -107,7 +116,7 @@ subroutine amg_d_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt
else
! Default transfer
call lv%linmap%map_V2U(alpha,vect_v,beta,vect_u,info,&
& work=work,vtx=vtx,vty=vty)
& work=work,vtx=vtx_,vty=vty)
end if
end subroutine amg_d_base_onelev_map_prol_v

46
amgprec/impl/level/amg_d_base_onelev_map_rstr.F90
Unescape Escape View File

@ -35,7 +35,6 @@
! POSSIBILITY OF SUCH DAMAGE.
!
!
subroutine amg_d_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
& work,vtx,vty)
use psb_base_mod
@ -47,17 +46,25 @@ subroutine amg_d_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
integer(psb_ipk_), intent(out) :: info
real(psb_dpk_), optional :: work(:)
type(psb_d_vect_type), optional, target, intent(inout) :: vtx,vty
type(psb_d_vect_type), pointer :: vty_
integer(psb_mpk_) :: me, np
!!$ write(0,*) 'New map_rstr',lv%remap_data%ac_pre_remap%is_asb()
write(0,*) 'New map_rstr ',lv%remap_data%ac_pre_remap%is_asb()
if (present(vty)) then
vty_ => vty
else
vty_ => lv%wrk%wv(1)
end if
if (lv%remap_data%ac_pre_remap%is_asb()) then
!
! Remap has happened, deal with it
!
!!$ write(0,*) 'Remap handling not implemented yet '
write(0,*) 'Level map_rstr with remapping '
block
type(psb_ctxt_type) :: ctxt, rctxt
integer(psb_mpk_) :: i,j,ip, idest, nsrc, nrl, kp
integer(psb_mpk_) :: me, np, rme, rnp
integer(psb_mpk_) :: rme, rnp
real(psb_dpk_), allocatable :: rsnd(:), rrcv(:)
type(psb_d_vect_type) :: tv
@ -65,24 +72,25 @@ subroutine amg_d_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
call psb_info(ctxt,me,np)
rctxt = lv%desc_ac%get_ctxt()
call psb_info(rctxt,rme,rnp)
!!$ write(0,*) 'New context ',rme,rnp
write(0,*) 'New context map rstr',rme,rnp,me,np
idest = lv%remap_data%idest
associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc)
!!$ write(0,*) 'Should apply maps, then send data from ',me,' to ',idest
!!$ if (rme >= 0) write(0,*) rme, ' Receiving data from ',isrc(:)
write(0,*) 'Should apply maps, then send data from ',me,' to ',idest
if (rme >= 0) write(0,*) rme, ' Receiving data from ',isrc(:)
nsrc = size(isrc)
nrl = lv%remap_data%desc_ac_pre_remap%get_local_rows()
call psb_geall(tv,lv%remap_data%desc_ac_pre_remap,info)
call psb_geasb(tv,lv%remap_data%desc_ac_pre_remap,info,mold=vect_u%v)
write(0,*) me,' map_rstr calling U2V: ',me,np,rme,rnp,tv%get_nrows()
write(0,*) me,' remap map_rstr calling U2V: ',me,np,rme,rnp,tv%get_nrows(),&
& psb_errstatus_fatal()
flush(0)
call psb_barrier(ctxt)
call lv%linmap%map_U2V(alpha,vect_u,beta,tv,info,&
& work=work,vtx=vtx,vty=vty)
call tv%sync()
& work=work,vtx=vtx,vty=vty_)
!call tv%sync()
!rsnd = tv%get_vect()
!call psb_snd(ctxt,rsnd(1:nrl),idest)
write(0,*) me,' map_rstr sending ',me,idest
write(0,*) me,' map_rstr sending ',me,idest,psb_errstatus_fatal()
call psb_snd(ctxt,tv%v%v(1:nrl),idest)
if (rme >=0) then
allocate(rrcv(sum(nrsrc)))
@ -91,22 +99,28 @@ subroutine amg_d_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
do i = 1,size(isrc)
ip = isrc(i)
nrl = nrsrc(i)
write(0,*) me,' map_rstr receiving',rme,i,ip
write(0,*) me,' map_rstr receiving',rme,ip,psb_errstatus_fatal()
call psb_rcv(ctxt,rrcv(kp+1:kp+nrl),ip)
kp = kp + nrl
end do
call vect_v%set_vect(rrcv)
end if
end associate
!!$ write(0,*) me, ' Restrictor with remap done '
write(0,*) me, ' Restrictor with remap done ',psb_errstatus_fatal()
end block
else
! Default transfer
call lv%linmap%map_U2V(alpha,vect_u,beta,vect_v,info,&
& work=work,vtx=vtx,vty=vty)
block
type(psb_ctxt_type) :: ctxt, rctxt
ctxt = lv%linmap%p_desc_U%get_ctxt()
call psb_info(ctxt,me,np)
write(0,*) me,' map_rstr calling U2V: ',me,np
call lv%linmap%map_U2V(alpha,vect_u,beta,vect_v,info,&
& work=work,vtx=vtx,vty=vty_)
end block
end if
write(0,*) me, 'End of restriction ',info,psb_errstatus_fatal()
end subroutine amg_d_base_onelev_map_rstr_v
subroutine amg_d_base_onelev_map_rstr_a(lv,alpha,u,beta,v,info,work)

13
amgprec/impl/level/amg_s_base_onelev_map_prol.F90
Unescape Escape View File

@ -46,6 +46,15 @@ subroutine amg_s_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt
integer(psb_ipk_), intent(out) :: info
real(psb_spk_), optional :: work(:)
type(psb_s_vect_type), optional, target, intent(inout) :: vtx,vty
type(psb_s_vect_type), pointer :: vtx_
!!$ write(0,*) 'New map_rstr',lv%remap_data%ac_pre_remap%is_asb()
if (present(vtx)) then
vtx_ => vtx
else
vtx_ => lv%wrk%wv(1)
end if
!!$ write(0,*) 'New map_prol',lv%remap_data%ac_pre_remap%is_asb()
if (lv%remap_data%ac_pre_remap%is_asb()) then
@ -98,7 +107,7 @@ subroutine amg_s_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt
call psb_rcv(ctxt,tv%v%v(1:nrl),idest)
call tv%set_host()
call lv%linmap%map_V2U(alpha,tv,beta,vect_u,info,&
& work=work,vtx=vtx,vty=vty)
& work=work,vtx=vtx_,vty=vty)
end associate
!!$ write(0,*) me, ' Prolongator with remap done '
!!$ flush(0)
@ -107,7 +116,7 @@ subroutine amg_s_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt
else
! Default transfer
call lv%linmap%map_V2U(alpha,vect_v,beta,vect_u,info,&
& work=work,vtx=vtx,vty=vty)
& work=work,vtx=vtx_,vty=vty)
end if
end subroutine amg_s_base_onelev_map_prol_v

22
amgprec/impl/level/amg_s_base_onelev_map_rstr.F90
Unescape Escape View File

@ -47,8 +47,15 @@ subroutine amg_s_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
integer(psb_ipk_), intent(out) :: info
real(psb_spk_), optional :: work(:)
type(psb_s_vect_type), optional, target, intent(inout) :: vtx,vty
type(psb_s_vect_type), pointer :: vty_
!!$ write(0,*) 'New map_rstr',lv%remap_data%ac_pre_remap%is_asb()
if (present(vty)) then
vty_ => vty
else
vty_ => lv%wrk%wv(1)
end if
if (lv%remap_data%ac_pre_remap%is_asb()) then
!
! Remap has happened, deal with it
@ -74,11 +81,11 @@ subroutine amg_s_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
nrl = lv%remap_data%desc_ac_pre_remap%get_local_rows()
call psb_geall(tv,lv%remap_data%desc_ac_pre_remap,info)
call psb_geasb(tv,lv%remap_data%desc_ac_pre_remap,info,mold=vect_u%v)
write(0,*) me,' map_rstr calling U2V: ',me,np,rme,rnp,tv%get_nrows()
write(0,*) me,' remap map_rstr calling U2V: ',me,np,rme,rnp,tv%get_nrows()
flush(0)
call psb_barrier(ctxt)
call lv%linmap%map_U2V(alpha,vect_u,beta,tv,info,&
& work=work,vtx=vtx,vty=vty)
& work=work,vtx=vtx,vty=vty_)
call tv%sync()
!rsnd = tv%get_vect()
!call psb_snd(ctxt,rsnd(1:nrl),idest)
@ -103,8 +110,15 @@ subroutine amg_s_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
else
! Default transfer
call lv%linmap%map_U2V(alpha,vect_u,beta,vect_v,info,&
& work=work,vtx=vtx,vty=vty)
block
type(psb_ctxt_type) :: ctxt, rctxt
integer(psb_mpk_) :: me, np
ctxt = lv%linmap%p_desc_U%get_ctxt()
call psb_info(ctxt,me,np)
write(0,*) me,' map_rstr calling U2V: ',me,np
call lv%linmap%map_U2V(alpha,vect_u,beta,vect_v,info,&
& work=work,vtx=vtx,vty=vty_)
end block
end if
end subroutine amg_s_base_onelev_map_rstr_v

13
amgprec/impl/level/amg_z_base_onelev_map_prol.F90
Unescape Escape View File

@ -46,6 +46,15 @@ subroutine amg_z_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt
integer(psb_ipk_), intent(out) :: info
complex(psb_dpk_), optional :: work(:)
type(psb_z_vect_type), optional, target, intent(inout) :: vtx,vty
type(psb_z_vect_type), pointer :: vtx_
!!$ write(0,*) 'New map_rstr',lv%remap_data%ac_pre_remap%is_asb()
if (present(vtx)) then
vtx_ => vtx
else
vtx_ => lv%wrk%wv(1)
end if
!!$ write(0,*) 'New map_prol',lv%remap_data%ac_pre_remap%is_asb()
if (lv%remap_data%ac_pre_remap%is_asb()) then
@ -98,7 +107,7 @@ subroutine amg_z_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt
call psb_rcv(ctxt,tv%v%v(1:nrl),idest)
call tv%set_host()
call lv%linmap%map_V2U(alpha,tv,beta,vect_u,info,&
& work=work,vtx=vtx,vty=vty)
& work=work,vtx=vtx_,vty=vty)
end associate
!!$ write(0,*) me, ' Prolongator with remap done '
!!$ flush(0)
@ -107,7 +116,7 @@ subroutine amg_z_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt
else
! Default transfer
call lv%linmap%map_V2U(alpha,vect_v,beta,vect_u,info,&
& work=work,vtx=vtx,vty=vty)
& work=work,vtx=vtx_,vty=vty)
end if
end subroutine amg_z_base_onelev_map_prol_v

22
amgprec/impl/level/amg_z_base_onelev_map_rstr.F90
Unescape Escape View File

@ -47,8 +47,15 @@ subroutine amg_z_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
integer(psb_ipk_), intent(out) :: info
complex(psb_dpk_), optional :: work(:)
type(psb_z_vect_type), optional, target, intent(inout) :: vtx,vty
type(psb_z_vect_type), pointer :: vty_
!!$ write(0,*) 'New map_rstr',lv%remap_data%ac_pre_remap%is_asb()
if (present(vty)) then
vty_ => vty
else
vty_ => lv%wrk%wv(1)
end if
if (lv%remap_data%ac_pre_remap%is_asb()) then
!
! Remap has happened, deal with it
@ -74,11 +81,11 @@ subroutine amg_z_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
nrl = lv%remap_data%desc_ac_pre_remap%get_local_rows()
call psb_geall(tv,lv%remap_data%desc_ac_pre_remap,info)
call psb_geasb(tv,lv%remap_data%desc_ac_pre_remap,info,mold=vect_u%v)
write(0,*) me,' map_rstr calling U2V: ',me,np,rme,rnp,tv%get_nrows()
write(0,*) me,' remap map_rstr calling U2V: ',me,np,rme,rnp,tv%get_nrows()
flush(0)
call psb_barrier(ctxt)
call lv%linmap%map_U2V(alpha,vect_u,beta,tv,info,&
& work=work,vtx=vtx,vty=vty)
& work=work,vtx=vtx,vty=vty_)
call tv%sync()
!rsnd = tv%get_vect()
!call psb_snd(ctxt,rsnd(1:nrl),idest)
@ -103,8 +110,15 @@ subroutine amg_z_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
else
! Default transfer
call lv%linmap%map_U2V(alpha,vect_u,beta,vect_v,info,&
& work=work,vtx=vtx,vty=vty)
block
type(psb_ctxt_type) :: ctxt, rctxt
integer(psb_mpk_) :: me, np
ctxt = lv%linmap%p_desc_U%get_ctxt()
call psb_info(ctxt,me,np)
write(0,*) me,' map_rstr calling U2V: ',me,np
call lv%linmap%map_U2V(alpha,vect_u,beta,vect_v,info,&
& work=work,vtx=vtx,vty=vty_)
end block
end if
end subroutine amg_z_base_onelev_map_rstr_v

12
samples/advanced/remap/amg_d_pde3d.F90
Unescape Escape View File

@ -472,6 +472,12 @@ program amg_d_pde3d
call psb_errpush(psb_err_from_subroutine_,name,a_err='amg_hierarchy_bld')
goto 9999
end if
do i=2, size(prec%precv)
write(0,*) iam,'Between hier and smoothers_bld level',i,':',&
& prec%precv(i)%remap_data%desc_ac_pre_remap%is_asb()
end do
call psb_barrier(ctxt)
t1 = psb_wtime()
call prec%smoothers_build(a,desc_a,info)
@ -491,6 +497,12 @@ program amg_d_pde3d
write(psb_out_unit,'(" ")')
end if
do i=2, size(prec%precv)
write(0,*) iam,'After smoothers_bld level',i,':',&
& prec%precv(i)%remap_data%desc_ac_pre_remap%is_asb()
end do
call prec%descr(info,iout=psb_out_unit)
if (p_choice%dump) then
call prec%dump(info,istart=p_choice%dlmin,iend=p_choice%dlmax,&

Write Preview
Loading…
Cancel
Save