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

Cleanup of first remapping version. To be tested.

Browse Source
remap-coarse
sfilippone 5 days ago
parent 23f84ce236
commit 8ac5b90454
46 changed files with 793 additions and 825 deletions
Show Stats Download Patch File Download Diff File

11
amgprec/amg_c_base_aggregator_mod.f90
Unescape Escape View File

@ -89,7 +89,7 @@ module amg_c_base_aggregator_mod
procedure, pass(ag) :: bld_tprol => amg_c_base_aggregator_build_tprol procedure, pass(ag) :: bld_tprol => amg_c_base_aggregator_build_tprol
procedure, pass(ag) :: mat_bld => amg_c_base_aggregator_mat_bld procedure, pass(ag) :: mat_bld => amg_c_base_aggregator_mat_bld
procedure, pass(ag) :: mat_asb => amg_c_base_aggregator_mat_asb procedure, pass(ag) :: mat_asb => amg_c_base_aggregator_mat_asb
procedure, pass(ag) :: bld_map => amg_c_base_aggregator_bld_map procedure, pass(ag) :: bld_linmap => amg_c_base_aggregator_bld_linmap
procedure, pass(ag) :: update_next => amg_c_base_aggregator_update_next procedure, pass(ag) :: update_next => amg_c_base_aggregator_update_next
procedure, pass(ag) :: clone => amg_c_base_aggregator_clone procedure, pass(ag) :: clone => amg_c_base_aggregator_clone
procedure, pass(ag) :: free => amg_c_base_aggregator_free procedure, pass(ag) :: free => amg_c_base_aggregator_free
@ -458,7 +458,7 @@ contains
end subroutine amg_c_base_aggregator_mat_asb end subroutine amg_c_base_aggregator_mat_asb
! !
!> Function bld_map !> Function bld_linmap
!! \memberof amg_c_base_aggregator_type !! \memberof amg_c_base_aggregator_type
!! \brief Build linear map between hierarchy levels !! \brief Build linear map between hierarchy levels
!! !!
@ -473,7 +473,7 @@ contains
!! \param map The output map !! \param map The output map
!! \param info Return code !! \param info Return code
!! !!
subroutine amg_c_base_aggregator_bld_map(ag,desc_a,desc_ac,ilaggr,nlaggr,& subroutine amg_c_base_aggregator_bld_linmap(ag,desc_a,desc_ac,ilaggr,nlaggr,&
& op_restr,op_prol,map,info) & op_restr,op_prol,map,info)
use psb_base_mod use psb_base_mod
implicit none implicit none
@ -484,7 +484,7 @@ contains
type(psb_clinmap_type), intent(out) :: map type(psb_clinmap_type), intent(out) :: map
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: err_act integer(psb_ipk_) :: err_act
character(len=20) :: name='c_base_aggregator_bld_map' character(len=20) :: name='c_base_aggregator_bld_linmap'
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -508,7 +508,6 @@ contains
9999 call psb_error_handler(err_act) 9999 call psb_error_handler(err_act)
return return
end subroutine amg_c_base_aggregator_bld_map end subroutine amg_c_base_aggregator_bld_linmap
end module amg_c_base_aggregator_mod end module amg_c_base_aggregator_mod

35
amgprec/amg_c_onelev_mod.f90
Unescape Escape View File

@ -181,7 +181,8 @@ module amg_c_onelev_mod
integer(psb_ipk_) :: idest integer(psb_ipk_) :: idest
integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:) integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:)
contains contains
procedure, pass(rmp) :: clone => c_remap_data_clone procedure, pass(rmp) :: clone => c_remap_data_clone
procedure, pass(rmp) :: move_alloc => c_remap_move_alloc
end type amg_c_remap_data_type end type amg_c_remap_data_type
type amg_c_onelev_type type amg_c_onelev_type
@ -227,7 +228,7 @@ module amg_c_onelev_mod
procedure, pass(lv) :: get_wrksz => c_base_onelev_get_wrksize procedure, pass(lv) :: get_wrksz => c_base_onelev_get_wrksize
procedure, pass(lv) :: allocate_wrk => c_base_onelev_allocate_wrk procedure, pass(lv) :: allocate_wrk => c_base_onelev_allocate_wrk
procedure, pass(lv) :: free_wrk => c_base_onelev_free_wrk procedure, pass(lv) :: free_wrk => c_base_onelev_free_wrk
procedure, nopass :: stringval => amg_stringval procedure, nopass :: stringval => amg_stringval
procedure, pass(lv) :: move_alloc => c_base_onelev_move_alloc procedure, pass(lv) :: move_alloc => c_base_onelev_move_alloc
@ -639,7 +640,7 @@ contains
! Arguments ! Arguments
class(amg_c_onelev_type), target, intent(inout) :: lv class(amg_c_onelev_type), target, intent(inout) :: lv
class(amg_c_onelev_type), target, intent(inout) :: lvout class(amg_c_onelev_type), target, intent(inout) :: lvout
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
info = psb_success_ info = psb_success_
if (allocated(lv%sm)) then if (allocated(lv%sm)) then
@ -705,6 +706,7 @@ contains
if (info == psb_success_) call psb_move_alloc(lv%tprol,b%tprol,info) 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%desc_ac,b%desc_ac,info)
if (info == psb_success_) call psb_move_alloc(lv%linmap,b%linmap,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_a => lv%base_a
b%base_desc => lv%base_desc b%base_desc => lv%base_desc
@ -759,6 +761,7 @@ contains
info = psb_success_ info = psb_success_
nwv = lv%get_wrksz() nwv = lv%get_wrksz()
if (.not.allocated(lv%wrk)) allocate(lv%wrk,stat=info) 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 (info == 0) then
if (lv%remap_data%desc_ac_pre_remap%is_asb()) then if (lv%remap_data%desc_ac_pre_remap%is_asb()) then
! !
@ -806,11 +809,12 @@ contains
info = psb_success_ info = psb_success_
call wk%free(info) 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) allocate(wk%wv(nwv),stat=info)
if (present(desc2).and.(desc%is_valid())) then if (present(desc2).and.(desc%is_valid())) then
write(0,*) 'wrk_alloc D2:',desc2%get_fmt(),desc2%is_asb() !!$ write(0,*) 'wrk_alloc D2:',desc2%get_fmt(),desc2%is_asb()
if (desc2%get_local_cols()>desc%get_local_cols()) then if (desc2%get_local_cols()>desc%get_local_cols()) then
call inner_do_wrk_alloc(wk,nwv,desc2,vmold=vmold) call inner_do_wrk_alloc(wk,nwv,desc2,vmold=vmold)
else else
@ -996,4 +1000,25 @@ contains
call psb_safe_ab_cpy(rmp%nrsrc,remap_out%nrsrc,info) call psb_safe_ab_cpy(rmp%nrsrc,remap_out%nrsrc,info)
end subroutine c_remap_data_clone end subroutine c_remap_data_clone
subroutine c_remap_move_alloc(rmp, remap_out, info)
use psb_base_mod
implicit none
! Arguments
class(amg_c_remap_data_type), target, intent(inout) :: rmp
class(amg_c_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 c_remap_move_alloc
end module amg_c_onelev_mod end module amg_c_onelev_mod

67
amgprec/amg_c_prec_type.f90
Unescape Escape View File

@ -102,6 +102,7 @@ module amg_c_prec_type
! The multilevel hierarchy ! The multilevel hierarchy
! !
type(amg_c_onelev_type), allocatable :: precv(:) type(amg_c_onelev_type), allocatable :: precv(:)
integer(psb_ipk_) :: nlevs
contains contains
procedure, pass(prec) :: psb_c_apply2_vect => amg_c_apply2_vect procedure, pass(prec) :: psb_c_apply2_vect => amg_c_apply2_vect
procedure, pass(prec) :: psb_c_apply1_vect => amg_c_apply1_vect procedure, pass(prec) :: psb_c_apply1_vect => amg_c_apply1_vect
@ -119,6 +120,7 @@ module amg_c_prec_type
procedure, pass(prec) :: cmp_complexity => amg_c_cmp_compl procedure, pass(prec) :: cmp_complexity => amg_c_cmp_compl
procedure, pass(prec) :: get_avg_cr => amg_c_get_avg_cr procedure, pass(prec) :: get_avg_cr => amg_c_get_avg_cr
procedure, pass(prec) :: cmp_avg_cr => amg_c_cmp_avg_cr procedure, pass(prec) :: cmp_avg_cr => amg_c_cmp_avg_cr
procedure, pass(prec) :: set_nlevs => amg_c_set_nlevs
procedure, pass(prec) :: get_nlevs => amg_c_get_nlevs procedure, pass(prec) :: get_nlevs => amg_c_get_nlevs
procedure, pass(prec) :: get_nzeros => amg_c_get_nzeros procedure, pass(prec) :: get_nzeros => amg_c_get_nzeros
procedure, pass(prec) :: sizeof => amg_cprec_sizeof procedure, pass(prec) :: sizeof => amg_cprec_sizeof
@ -435,10 +437,19 @@ contains
class(amg_cprec_type), intent(in) :: prec class(amg_cprec_type), intent(in) :: prec
integer(psb_ipk_) :: val integer(psb_ipk_) :: val
val = 0 val = 0
if (allocated(prec%precv)) then !!$ if (allocated(prec%precv)) then
val = size(prec%precv) !!$ val = size(prec%precv)
end if !!$ end if
val = prec%nlevs
!!$ write(0,*) ' NLEVS: ',prec%nlevs, val,size(prec%precv)
end function amg_c_get_nlevs end function amg_c_get_nlevs
subroutine amg_c_set_nlevs(prec,nl)
implicit none
class(amg_cprec_type), intent(inout) :: prec
integer(psb_ipk_) :: nl
prec%nlevs = nl
end subroutine amg_c_set_nlevs
! !
! Function returning the size of the amg_prec_type data structure ! Function returning the size of the amg_prec_type data structure
! in bytes or in number of nonzeros of the operator(s) involved. ! in bytes or in number of nonzeros of the operator(s) involved.
@ -509,7 +520,7 @@ contains
real(psb_spk_) :: num, den, nmin real(psb_spk_) :: num, den, nmin
type(psb_ctxt_type) :: ctxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: il integer(psb_ipk_) :: il,nl
num = -sone num = -sone
den = sone den = sone
@ -519,7 +530,10 @@ contains
num = prec%precv(il)%base_a%get_nzeros() num = prec%precv(il)%base_a%get_nzeros()
if (num >= szero) then if (num >= szero) then
den = num 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()) num = num + max(0,prec%precv(il)%base_a%get_nzeros())
end do end do
end if end if
@ -550,7 +564,6 @@ contains
end function amg_c_get_avg_cr end function amg_c_get_avg_cr
subroutine amg_c_cmp_avg_cr(prec) subroutine amg_c_cmp_avg_cr(prec)
implicit none implicit none
class(amg_cprec_type), intent(inout) :: prec class(amg_cprec_type), intent(inout) :: prec
@ -558,17 +571,14 @@ contains
type(psb_ctxt_type) :: ctxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: il, nl, iam, np integer(psb_ipk_) :: il, nl, iam, np
avgcr = szero avgcr = szero
ctxt = prec%ctxt ctxt = prec%ctxt
call psb_info(ctxt,iam,np) call psb_info(ctxt,iam,np)
if (allocated(prec%precv)) then nl = prec%get_nlevs()
nl = size(prec%precv) do il=2,nl
do il=2,nl avgcr = avgcr + max(szero,prec%precv(il)%szratio)
avgcr = avgcr + max(szero,prec%precv(il)%szratio) end do
end do avgcr = avgcr / (nl-1)
avgcr = avgcr / (nl-1)
end if
call psb_sum(ctxt,avgcr) call psb_sum(ctxt,avgcr)
prec%ag_data%avg_cr = avgcr/np prec%ag_data%avg_cr = avgcr/np
end subroutine amg_c_cmp_avg_cr end subroutine amg_c_cmp_avg_cr
@ -586,9 +596,7 @@ contains
! error code. ! error code.
! !
subroutine amg_cprecfree(p,info) subroutine amg_cprecfree(p,info)
implicit none implicit none
! Arguments ! Arguments
type(amg_cprec_type), intent(inout) :: p type(amg_cprec_type), intent(inout) :: p
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -614,9 +622,7 @@ contains
end subroutine amg_cprecfree end subroutine amg_cprecfree
subroutine amg_c_prec_free(prec,info) subroutine amg_c_prec_free(prec,info)
implicit none implicit none
! Arguments ! Arguments
class(amg_cprec_type), intent(inout) :: prec class(amg_cprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -648,9 +654,7 @@ contains
end subroutine amg_c_prec_free end subroutine amg_c_prec_free
subroutine amg_c_smoothers_free(prec,info) subroutine amg_c_smoothers_free(prec,info)
implicit none implicit none
! Arguments ! Arguments
class(amg_cprec_type), intent(inout) :: prec class(amg_cprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -667,7 +671,7 @@ contains
end if end if
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
do i=1,size(prec%precv) do i=1,prec%get_nlevs()
call prec%precv(i)%free_smoothers(info) call prec%precv(i)%free_smoothers(info)
end do end do
end if end if
@ -681,9 +685,7 @@ contains
end subroutine amg_c_smoothers_free end subroutine amg_c_smoothers_free
subroutine amg_c_hierarchy_free(prec,info) subroutine amg_c_hierarchy_free(prec,info)
implicit none implicit none
! Arguments ! Arguments
class(amg_cprec_type), intent(inout) :: prec class(amg_cprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -710,7 +712,6 @@ contains
end subroutine amg_c_hierarchy_free end subroutine amg_c_hierarchy_free
! !
! Top level methods. ! Top level methods.
! !
@ -775,7 +776,6 @@ contains
end subroutine amg_c_apply1_vect end subroutine amg_c_apply1_vect
subroutine amg_c_apply2v(prec,x,y,desc_data,info,trans,work) subroutine amg_c_apply2v(prec,x,y,desc_data,info,trans,work)
implicit none implicit none
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
@ -840,7 +840,6 @@ contains
subroutine amg_c_dump(prec,info,istart,iend,iproc,prefix,head,& subroutine amg_c_dump(prec,info,istart,iend,iproc,prefix,head,&
& ac,rp,smoother,solver,tprol,& & ac,rp,smoother,solver,tprol,&
& global_num) & global_num)
implicit none implicit none
class(amg_cprec_type), intent(in) :: prec class(amg_cprec_type), intent(in) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -857,7 +856,7 @@ contains
info = 0 info = 0
ctxt = prec%ctxt ctxt = prec%ctxt
call psb_info(ctxt,iam,np) call psb_info(ctxt,iam,np)
iln = size(prec%precv) iln = prec%get_nlevs()
if (present(istart)) then if (present(istart)) then
il1 = max(1,istart) il1 = max(1,istart)
else else
@ -881,7 +880,6 @@ contains
end subroutine amg_c_dump end subroutine amg_c_dump
subroutine amg_c_cnv(prec,info,amold,vmold,imold) subroutine amg_c_cnv(prec,info,amold,vmold,imold)
implicit none implicit none
class(amg_cprec_type), intent(inout) :: prec class(amg_cprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -893,7 +891,7 @@ contains
info = psb_success_ info = psb_success_
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
do i=1,size(prec%precv) do i=1,prec%get_nlevs()
if (info == psb_success_ ) & if (info == psb_success_ ) &
& call prec%precv(i)%cnv(info,amold=amold,vmold=vmold,imold=imold) & call prec%precv(i)%cnv(info,amold=amold,vmold=vmold,imold=imold)
end do end do
@ -902,7 +900,6 @@ contains
end subroutine amg_c_cnv end subroutine amg_c_cnv
subroutine amg_c_clone(prec,precout,info) subroutine amg_c_clone(prec,precout,info)
implicit none implicit none
class(amg_cprec_type), intent(inout) :: prec class(amg_cprec_type), intent(inout) :: prec
class(psb_cprec_type), intent(inout) :: precout class(psb_cprec_type), intent(inout) :: precout
@ -914,7 +911,6 @@ contains
end subroutine amg_c_clone end subroutine amg_c_clone
subroutine amg_c_inner_clone(prec,precout,info) subroutine amg_c_inner_clone(prec,precout,info)
implicit none implicit none
class(amg_cprec_type), intent(inout) :: prec class(amg_cprec_type), intent(inout) :: prec
class(psb_cprec_type), target, intent(inout) :: precout class(psb_cprec_type), target, intent(inout) :: precout
@ -930,8 +926,9 @@ contains
pout%ctxt = prec%ctxt pout%ctxt = prec%ctxt
pout%ag_data = prec%ag_data pout%ag_data = prec%ag_data
pout%outer_sweeps = prec%outer_sweeps pout%outer_sweeps = prec%outer_sweeps
pout%nlevs = prec%nlevs
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
ln = size(prec%precv) ln = prec%get_nlevs()
allocate(pout%precv(ln),stat=info) allocate(pout%precv(ln),stat=info)
if (info /= psb_success_) goto 9999 if (info /= psb_success_) goto 9999
if (ln >= 1) then if (ln >= 1) then
@ -939,6 +936,7 @@ contains
end if end if
do lev=2, ln do lev=2, ln
if (info /= psb_success_) exit if (info /= psb_success_) exit
!!$ write(0,*) 'Inner_clone must be checked and reimplemented! '
call prec%precv(lev)%clone(pout%precv(lev),info) call prec%precv(lev)%clone(pout%precv(lev),info)
if (info == psb_success_) then if (info == psb_success_) then
pout%precv(lev)%base_a => pout%precv(lev)%ac pout%precv(lev)%base_a => pout%precv(lev)%ac
@ -1018,7 +1016,7 @@ contains
if (psb_errstatus_fatal()) then if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999 info = psb_err_internal_error_; goto 9999
end if end if
nlev = size(prec%precv) nlev = prec%get_nlevs()
level = 1 level = 1
do level = 1, nlev do level = 1, nlev
call prec%precv(level)%allocate_wrk(info,vmold=vmold) call prec%precv(level)%allocate_wrk(info,vmold=vmold)
@ -1041,7 +1039,6 @@ contains
subroutine amg_c_free_wrk(prec,info) subroutine amg_c_free_wrk(prec,info)
use psb_base_mod use psb_base_mod
implicit none implicit none
! Arguments ! Arguments
class(amg_cprec_type), intent(inout) :: prec class(amg_cprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -1058,7 +1055,7 @@ contains
end if end if
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
nlev = size(prec%precv) nlev = prec%get_nlevs()
do level = 1, nlev do level = 1, nlev
call prec%precv(level)%free_wrk(info) call prec%precv(level)%free_wrk(info)
end do end do

11
amgprec/amg_d_base_aggregator_mod.f90
Unescape Escape View File

@ -89,7 +89,7 @@ module amg_d_base_aggregator_mod
procedure, pass(ag) :: bld_tprol => amg_d_base_aggregator_build_tprol procedure, pass(ag) :: bld_tprol => amg_d_base_aggregator_build_tprol
procedure, pass(ag) :: mat_bld => amg_d_base_aggregator_mat_bld procedure, pass(ag) :: mat_bld => amg_d_base_aggregator_mat_bld
procedure, pass(ag) :: mat_asb => amg_d_base_aggregator_mat_asb procedure, pass(ag) :: mat_asb => amg_d_base_aggregator_mat_asb
procedure, pass(ag) :: bld_map => amg_d_base_aggregator_bld_map procedure, pass(ag) :: bld_linmap => amg_d_base_aggregator_bld_linmap
procedure, pass(ag) :: update_next => amg_d_base_aggregator_update_next procedure, pass(ag) :: update_next => amg_d_base_aggregator_update_next
procedure, pass(ag) :: clone => amg_d_base_aggregator_clone procedure, pass(ag) :: clone => amg_d_base_aggregator_clone
procedure, pass(ag) :: free => amg_d_base_aggregator_free procedure, pass(ag) :: free => amg_d_base_aggregator_free
@ -458,7 +458,7 @@ contains
end subroutine amg_d_base_aggregator_mat_asb end subroutine amg_d_base_aggregator_mat_asb
! !
!> Function bld_map !> Function bld_linmap
!! \memberof amg_d_base_aggregator_type !! \memberof amg_d_base_aggregator_type
!! \brief Build linear map between hierarchy levels !! \brief Build linear map between hierarchy levels
!! !!
@ -473,7 +473,7 @@ contains
!! \param map The output map !! \param map The output map
!! \param info Return code !! \param info Return code
!! !!
subroutine amg_d_base_aggregator_bld_map(ag,desc_a,desc_ac,ilaggr,nlaggr,& subroutine amg_d_base_aggregator_bld_linmap(ag,desc_a,desc_ac,ilaggr,nlaggr,&
& op_restr,op_prol,map,info) & op_restr,op_prol,map,info)
use psb_base_mod use psb_base_mod
implicit none implicit none
@ -484,7 +484,7 @@ contains
type(psb_dlinmap_type), intent(out) :: map type(psb_dlinmap_type), intent(out) :: map
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: err_act integer(psb_ipk_) :: err_act
character(len=20) :: name='d_base_aggregator_bld_map' character(len=20) :: name='d_base_aggregator_bld_linmap'
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -508,7 +508,6 @@ contains
9999 call psb_error_handler(err_act) 9999 call psb_error_handler(err_act)
return return
end subroutine amg_d_base_aggregator_bld_map end subroutine amg_d_base_aggregator_bld_linmap
end module amg_d_base_aggregator_mod end module amg_d_base_aggregator_mod

8
amgprec/amg_d_onelev_mod.f90
Unescape Escape View File

@ -182,8 +182,8 @@ module amg_d_onelev_mod
integer(psb_ipk_) :: idest integer(psb_ipk_) :: idest
integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:) integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:)
contains contains
procedure, pass(rmp) :: clone => d_remap_data_clone procedure, pass(rmp) :: clone => d_remap_data_clone
procedure, pass(rmp) :: move_alloc => d_remap_move_alloc procedure, pass(rmp) :: move_alloc => d_remap_move_alloc
end type amg_d_remap_data_type end type amg_d_remap_data_type
type amg_d_onelev_type type amg_d_onelev_type
@ -229,7 +229,7 @@ module amg_d_onelev_mod
procedure, pass(lv) :: get_wrksz => d_base_onelev_get_wrksize procedure, pass(lv) :: get_wrksz => d_base_onelev_get_wrksize
procedure, pass(lv) :: allocate_wrk => d_base_onelev_allocate_wrk procedure, pass(lv) :: allocate_wrk => d_base_onelev_allocate_wrk
procedure, pass(lv) :: free_wrk => d_base_onelev_free_wrk procedure, pass(lv) :: free_wrk => d_base_onelev_free_wrk
procedure, nopass :: stringval => amg_stringval procedure, nopass :: stringval => amg_stringval
procedure, pass(lv) :: move_alloc => d_base_onelev_move_alloc procedure, pass(lv) :: move_alloc => d_base_onelev_move_alloc
@ -641,7 +641,7 @@ contains
! Arguments ! Arguments
class(amg_d_onelev_type), target, intent(inout) :: lv class(amg_d_onelev_type), target, intent(inout) :: lv
class(amg_d_onelev_type), target, intent(inout) :: lvout class(amg_d_onelev_type), target, intent(inout) :: lvout
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
info = psb_success_ info = psb_success_
if (allocated(lv%sm)) then if (allocated(lv%sm)) then

8
amgprec/amg_d_parmatch_aggregator_mod.F90
Unescape Escape View File

@ -136,7 +136,7 @@ module amg_d_parmatch_aggregator_mod
procedure, pass(ag) :: mat_bld => amg_d_parmatch_aggregator_mat_bld procedure, pass(ag) :: mat_bld => amg_d_parmatch_aggregator_mat_bld
procedure, pass(ag) :: mat_asb => amg_d_parmatch_aggregator_mat_asb procedure, pass(ag) :: mat_asb => amg_d_parmatch_aggregator_mat_asb
procedure, pass(ag) :: inner_mat_asb => amg_d_parmatch_aggregator_inner_mat_asb procedure, pass(ag) :: inner_mat_asb => amg_d_parmatch_aggregator_inner_mat_asb
procedure, pass(ag) :: bld_map => amg_d_parmatch_aggregator_bld_map procedure, pass(ag) :: bld_linmap => amg_d_parmatch_aggregator_bld_linmap
procedure, pass(ag) :: csetc => amg_d_parmatch_aggr_csetc procedure, pass(ag) :: csetc => amg_d_parmatch_aggr_csetc
procedure, pass(ag) :: cseti => amg_d_parmatch_aggr_cseti procedure, pass(ag) :: cseti => amg_d_parmatch_aggr_cseti
procedure, pass(ag) :: default => amg_d_parmatch_aggr_set_default procedure, pass(ag) :: default => amg_d_parmatch_aggr_set_default
@ -643,7 +643,7 @@ contains
end select end select
end subroutine amg_d_parmatch_aggregator_clone end subroutine amg_d_parmatch_aggregator_clone
subroutine amg_d_parmatch_aggregator_bld_map(ag,desc_a,desc_ac,ilaggr,nlaggr,& subroutine amg_d_parmatch_aggregator_bld_linmap(ag,desc_a,desc_ac,ilaggr,nlaggr,&
& op_restr,op_prol,map,info) & op_restr,op_prol,map,info)
use psb_base_mod use psb_base_mod
implicit none implicit none
@ -654,7 +654,7 @@ contains
type(psb_dlinmap_type), intent(out) :: map type(psb_dlinmap_type), intent(out) :: map
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: err_act integer(psb_ipk_) :: err_act
character(len=20) :: name='d_parmatch_aggregator_bld_map' character(len=20) :: name='d_parmatch_aggregator_bld_linmap'
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -680,5 +680,5 @@ contains
9999 call psb_error_handler(err_act) 9999 call psb_error_handler(err_act)
return return
end subroutine amg_d_parmatch_aggregator_bld_map end subroutine amg_d_parmatch_aggregator_bld_linmap
end module amg_d_parmatch_aggregator_mod end module amg_d_parmatch_aggregator_mod

21
amgprec/amg_d_prec_type.f90
Unescape Escape View File

@ -97,6 +97,7 @@ module amg_d_prec_type
! to keep track against what is put later in the multilevel array ! to keep track against what is put later in the multilevel array
! !
integer(psb_ipk_) :: coarse_solver = -1 integer(psb_ipk_) :: coarse_solver = -1
! !
! The multilevel hierarchy ! The multilevel hierarchy
! !
@ -119,7 +120,7 @@ module amg_d_prec_type
procedure, pass(prec) :: cmp_complexity => amg_d_cmp_compl procedure, pass(prec) :: cmp_complexity => amg_d_cmp_compl
procedure, pass(prec) :: get_avg_cr => amg_d_get_avg_cr 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) :: cmp_avg_cr => amg_d_cmp_avg_cr
procedure, pass(prec) :: set_nlevs => amg_d_set_nlevs procedure, pass(prec) :: set_nlevs => amg_d_set_nlevs
procedure, pass(prec) :: get_nlevs => amg_d_get_nlevs procedure, pass(prec) :: get_nlevs => amg_d_get_nlevs
procedure, pass(prec) :: get_nzeros => amg_d_get_nzeros procedure, pass(prec) :: get_nzeros => amg_d_get_nzeros
procedure, pass(prec) :: sizeof => amg_dprec_sizeof procedure, pass(prec) :: sizeof => amg_dprec_sizeof
@ -142,7 +143,6 @@ module amg_d_prec_type
procedure, pass(prec) :: smoothers_free => amg_d_smoothers_free procedure, pass(prec) :: smoothers_free => amg_d_smoothers_free
procedure, pass(prec) :: descr => amg_dfile_prec_descr procedure, pass(prec) :: descr => amg_dfile_prec_descr
procedure, pass(prec) :: memory_use => amg_dfile_prec_memory_use procedure, pass(prec) :: memory_use => amg_dfile_prec_memory_use
end type amg_dprec_type end type amg_dprec_type
private :: amg_d_dump, amg_d_get_compl, amg_d_cmp_compl,& private :: amg_d_dump, amg_d_get_compl, amg_d_cmp_compl,&
@ -442,7 +442,6 @@ contains
!!$ end if !!$ end if
val = prec%nlevs val = prec%nlevs
!!$ write(0,*) ' NLEVS: ',prec%nlevs, val,size(prec%precv) !!$ write(0,*) ' NLEVS: ',prec%nlevs, val,size(prec%precv)
end function amg_d_get_nlevs end function amg_d_get_nlevs
subroutine amg_d_set_nlevs(prec,nl) subroutine amg_d_set_nlevs(prec,nl)
@ -565,7 +564,6 @@ contains
end function amg_d_get_avg_cr end function amg_d_get_avg_cr
subroutine amg_d_cmp_avg_cr(prec) subroutine amg_d_cmp_avg_cr(prec)
implicit none implicit none
class(amg_dprec_type), intent(inout) :: prec class(amg_dprec_type), intent(inout) :: prec
@ -573,7 +571,6 @@ contains
type(psb_ctxt_type) :: ctxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: il, nl, iam, np integer(psb_ipk_) :: il, nl, iam, np
avgcr = dzero avgcr = dzero
ctxt = prec%ctxt ctxt = prec%ctxt
call psb_info(ctxt,iam,np) call psb_info(ctxt,iam,np)
@ -599,9 +596,7 @@ contains
! error code. ! error code.
! !
subroutine amg_dprecfree(p,info) subroutine amg_dprecfree(p,info)
implicit none implicit none
! Arguments ! Arguments
type(amg_dprec_type), intent(inout) :: p type(amg_dprec_type), intent(inout) :: p
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -690,9 +685,7 @@ contains
end subroutine amg_d_smoothers_free end subroutine amg_d_smoothers_free
subroutine amg_d_hierarchy_free(prec,info) subroutine amg_d_hierarchy_free(prec,info)
implicit none implicit none
! Arguments ! Arguments
class(amg_dprec_type), intent(inout) :: prec class(amg_dprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -719,7 +712,6 @@ contains
end subroutine amg_d_hierarchy_free end subroutine amg_d_hierarchy_free
! !
! Top level methods. ! Top level methods.
! !
@ -784,7 +776,6 @@ contains
end subroutine amg_d_apply1_vect end subroutine amg_d_apply1_vect
subroutine amg_d_apply2v(prec,x,y,desc_data,info,trans,work) subroutine amg_d_apply2v(prec,x,y,desc_data,info,trans,work)
implicit none implicit none
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
@ -849,7 +840,6 @@ contains
subroutine amg_d_dump(prec,info,istart,iend,iproc,prefix,head,& subroutine amg_d_dump(prec,info,istart,iend,iproc,prefix,head,&
& ac,rp,smoother,solver,tprol,& & ac,rp,smoother,solver,tprol,&
& global_num) & global_num)
implicit none implicit none
class(amg_dprec_type), intent(in) :: prec class(amg_dprec_type), intent(in) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -890,7 +880,6 @@ contains
end subroutine amg_d_dump end subroutine amg_d_dump
subroutine amg_d_cnv(prec,info,amold,vmold,imold) subroutine amg_d_cnv(prec,info,amold,vmold,imold)
implicit none implicit none
class(amg_dprec_type), intent(inout) :: prec class(amg_dprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -911,7 +900,6 @@ contains
end subroutine amg_d_cnv end subroutine amg_d_cnv
subroutine amg_d_clone(prec,precout,info) subroutine amg_d_clone(prec,precout,info)
implicit none implicit none
class(amg_dprec_type), intent(inout) :: prec class(amg_dprec_type), intent(inout) :: prec
class(psb_dprec_type), intent(inout) :: precout class(psb_dprec_type), intent(inout) :: precout
@ -948,7 +936,7 @@ contains
end if end if
do lev=2, ln do lev=2, ln
if (info /= psb_success_) exit if (info /= psb_success_) exit
write(0,*) 'Inner_clone must be checked and reimplemented! ' !!$ write(0,*) 'Inner_clone must be checked and reimplemented! '
call prec%precv(lev)%clone(pout%precv(lev),info) call prec%precv(lev)%clone(pout%precv(lev),info)
if (info == psb_success_) then if (info == psb_success_) then
pout%precv(lev)%base_a => pout%precv(lev)%ac pout%precv(lev)%base_a => pout%precv(lev)%ac
@ -1007,6 +995,7 @@ contains
subroutine amg_d_allocate_wrk(prec,info,vmold,desc) subroutine amg_d_allocate_wrk(prec,info,vmold,desc)
use psb_base_mod use psb_base_mod
implicit none implicit none
! Arguments ! Arguments
class(amg_dprec_type), intent(inout) :: prec class(amg_dprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -1050,7 +1039,6 @@ contains
subroutine amg_d_free_wrk(prec,info) subroutine amg_d_free_wrk(prec,info)
use psb_base_mod use psb_base_mod
implicit none implicit none
! Arguments ! Arguments
class(amg_dprec_type), intent(inout) :: prec class(amg_dprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -1084,6 +1072,7 @@ contains
function amg_d_is_allocated_wrk(prec) result(res) function amg_d_is_allocated_wrk(prec) result(res)
use psb_base_mod use psb_base_mod
implicit none implicit none
! Arguments ! Arguments
class(amg_dprec_type), intent(in) :: prec class(amg_dprec_type), intent(in) :: prec
logical :: res logical :: res

11
amgprec/amg_s_base_aggregator_mod.f90
Unescape Escape View File

@ -89,7 +89,7 @@ module amg_s_base_aggregator_mod
procedure, pass(ag) :: bld_tprol => amg_s_base_aggregator_build_tprol procedure, pass(ag) :: bld_tprol => amg_s_base_aggregator_build_tprol
procedure, pass(ag) :: mat_bld => amg_s_base_aggregator_mat_bld procedure, pass(ag) :: mat_bld => amg_s_base_aggregator_mat_bld
procedure, pass(ag) :: mat_asb => amg_s_base_aggregator_mat_asb procedure, pass(ag) :: mat_asb => amg_s_base_aggregator_mat_asb
procedure, pass(ag) :: bld_map => amg_s_base_aggregator_bld_map procedure, pass(ag) :: bld_linmap => amg_s_base_aggregator_bld_linmap
procedure, pass(ag) :: update_next => amg_s_base_aggregator_update_next procedure, pass(ag) :: update_next => amg_s_base_aggregator_update_next
procedure, pass(ag) :: clone => amg_s_base_aggregator_clone procedure, pass(ag) :: clone => amg_s_base_aggregator_clone
procedure, pass(ag) :: free => amg_s_base_aggregator_free procedure, pass(ag) :: free => amg_s_base_aggregator_free
@ -458,7 +458,7 @@ contains
end subroutine amg_s_base_aggregator_mat_asb end subroutine amg_s_base_aggregator_mat_asb
! !
!> Function bld_map !> Function bld_linmap
!! \memberof amg_s_base_aggregator_type !! \memberof amg_s_base_aggregator_type
!! \brief Build linear map between hierarchy levels !! \brief Build linear map between hierarchy levels
!! !!
@ -473,7 +473,7 @@ contains
!! \param map The output map !! \param map The output map
!! \param info Return code !! \param info Return code
!! !!
subroutine amg_s_base_aggregator_bld_map(ag,desc_a,desc_ac,ilaggr,nlaggr,& subroutine amg_s_base_aggregator_bld_linmap(ag,desc_a,desc_ac,ilaggr,nlaggr,&
& op_restr,op_prol,map,info) & op_restr,op_prol,map,info)
use psb_base_mod use psb_base_mod
implicit none implicit none
@ -484,7 +484,7 @@ contains
type(psb_slinmap_type), intent(out) :: map type(psb_slinmap_type), intent(out) :: map
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: err_act integer(psb_ipk_) :: err_act
character(len=20) :: name='s_base_aggregator_bld_map' character(len=20) :: name='s_base_aggregator_bld_linmap'
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -508,7 +508,6 @@ contains
9999 call psb_error_handler(err_act) 9999 call psb_error_handler(err_act)
return return
end subroutine amg_s_base_aggregator_bld_map end subroutine amg_s_base_aggregator_bld_linmap
end module amg_s_base_aggregator_mod end module amg_s_base_aggregator_mod

35
amgprec/amg_s_onelev_mod.f90
Unescape Escape View File

@ -182,7 +182,8 @@ module amg_s_onelev_mod
integer(psb_ipk_) :: idest integer(psb_ipk_) :: idest
integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:) integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:)
contains contains
procedure, pass(rmp) :: clone => s_remap_data_clone procedure, pass(rmp) :: clone => s_remap_data_clone
procedure, pass(rmp) :: move_alloc => s_remap_move_alloc
end type amg_s_remap_data_type end type amg_s_remap_data_type
type amg_s_onelev_type type amg_s_onelev_type
@ -228,7 +229,7 @@ module amg_s_onelev_mod
procedure, pass(lv) :: get_wrksz => s_base_onelev_get_wrksize procedure, pass(lv) :: get_wrksz => s_base_onelev_get_wrksize
procedure, pass(lv) :: allocate_wrk => s_base_onelev_allocate_wrk procedure, pass(lv) :: allocate_wrk => s_base_onelev_allocate_wrk
procedure, pass(lv) :: free_wrk => s_base_onelev_free_wrk procedure, pass(lv) :: free_wrk => s_base_onelev_free_wrk
procedure, nopass :: stringval => amg_stringval procedure, nopass :: stringval => amg_stringval
procedure, pass(lv) :: move_alloc => s_base_onelev_move_alloc procedure, pass(lv) :: move_alloc => s_base_onelev_move_alloc
@ -640,7 +641,7 @@ contains
! Arguments ! Arguments
class(amg_s_onelev_type), target, intent(inout) :: lv class(amg_s_onelev_type), target, intent(inout) :: lv
class(amg_s_onelev_type), target, intent(inout) :: lvout class(amg_s_onelev_type), target, intent(inout) :: lvout
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
info = psb_success_ info = psb_success_
if (allocated(lv%sm)) then if (allocated(lv%sm)) then
@ -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%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%desc_ac,b%desc_ac,info)
if (info == psb_success_) call psb_move_alloc(lv%linmap,b%linmap,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_a => lv%base_a
b%base_desc => lv%base_desc b%base_desc => lv%base_desc
@ -760,6 +762,7 @@ contains
info = psb_success_ info = psb_success_
nwv = lv%get_wrksz() nwv = lv%get_wrksz()
if (.not.allocated(lv%wrk)) allocate(lv%wrk,stat=info) 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 (info == 0) then
if (lv%remap_data%desc_ac_pre_remap%is_asb()) then if (lv%remap_data%desc_ac_pre_remap%is_asb()) then
! !
@ -807,11 +810,12 @@ contains
info = psb_success_ info = psb_success_
call wk%free(info) 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) allocate(wk%wv(nwv),stat=info)
if (present(desc2).and.(desc%is_valid())) then if (present(desc2).and.(desc%is_valid())) then
write(0,*) 'wrk_alloc D2:',desc2%get_fmt(),desc2%is_asb() !!$ write(0,*) 'wrk_alloc D2:',desc2%get_fmt(),desc2%is_asb()
if (desc2%get_local_cols()>desc%get_local_cols()) then if (desc2%get_local_cols()>desc%get_local_cols()) then
call inner_do_wrk_alloc(wk,nwv,desc2,vmold=vmold) call inner_do_wrk_alloc(wk,nwv,desc2,vmold=vmold)
else else
@ -997,4 +1001,25 @@ contains
call psb_safe_ab_cpy(rmp%nrsrc,remap_out%nrsrc,info) call psb_safe_ab_cpy(rmp%nrsrc,remap_out%nrsrc,info)
end subroutine s_remap_data_clone end subroutine s_remap_data_clone
subroutine s_remap_move_alloc(rmp, remap_out, info)
use psb_base_mod
implicit none
! Arguments
class(amg_s_remap_data_type), target, intent(inout) :: rmp
class(amg_s_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 s_remap_move_alloc
end module amg_s_onelev_mod end module amg_s_onelev_mod

8
amgprec/amg_s_parmatch_aggregator_mod.F90
Unescape Escape View File

@ -136,7 +136,7 @@ module amg_s_parmatch_aggregator_mod
procedure, pass(ag) :: mat_bld => amg_s_parmatch_aggregator_mat_bld procedure, pass(ag) :: mat_bld => amg_s_parmatch_aggregator_mat_bld
procedure, pass(ag) :: mat_asb => amg_s_parmatch_aggregator_mat_asb procedure, pass(ag) :: mat_asb => amg_s_parmatch_aggregator_mat_asb
procedure, pass(ag) :: inner_mat_asb => amg_s_parmatch_aggregator_inner_mat_asb procedure, pass(ag) :: inner_mat_asb => amg_s_parmatch_aggregator_inner_mat_asb
procedure, pass(ag) :: bld_map => amg_s_parmatch_aggregator_bld_map procedure, pass(ag) :: bld_linmap => amg_s_parmatch_aggregator_bld_linmap
procedure, pass(ag) :: csetc => amg_s_parmatch_aggr_csetc procedure, pass(ag) :: csetc => amg_s_parmatch_aggr_csetc
procedure, pass(ag) :: cseti => amg_s_parmatch_aggr_cseti procedure, pass(ag) :: cseti => amg_s_parmatch_aggr_cseti
procedure, pass(ag) :: default => amg_s_parmatch_aggr_set_default procedure, pass(ag) :: default => amg_s_parmatch_aggr_set_default
@ -643,7 +643,7 @@ contains
end select end select
end subroutine amg_s_parmatch_aggregator_clone end subroutine amg_s_parmatch_aggregator_clone
subroutine amg_s_parmatch_aggregator_bld_map(ag,desc_a,desc_ac,ilaggr,nlaggr,& subroutine amg_s_parmatch_aggregator_bld_linmap(ag,desc_a,desc_ac,ilaggr,nlaggr,&
& op_restr,op_prol,map,info) & op_restr,op_prol,map,info)
use psb_base_mod use psb_base_mod
implicit none implicit none
@ -654,7 +654,7 @@ contains
type(psb_slinmap_type), intent(out) :: map type(psb_slinmap_type), intent(out) :: map
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: err_act integer(psb_ipk_) :: err_act
character(len=20) :: name='s_parmatch_aggregator_bld_map' character(len=20) :: name='s_parmatch_aggregator_bld_linmap'
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -680,5 +680,5 @@ contains
9999 call psb_error_handler(err_act) 9999 call psb_error_handler(err_act)
return return
end subroutine amg_s_parmatch_aggregator_bld_map end subroutine amg_s_parmatch_aggregator_bld_linmap
end module amg_s_parmatch_aggregator_mod end module amg_s_parmatch_aggregator_mod

67
amgprec/amg_s_prec_type.f90
Unescape Escape View File

@ -102,6 +102,7 @@ module amg_s_prec_type
! The multilevel hierarchy ! The multilevel hierarchy
! !
type(amg_s_onelev_type), allocatable :: precv(:) type(amg_s_onelev_type), allocatable :: precv(:)
integer(psb_ipk_) :: nlevs
contains contains
procedure, pass(prec) :: psb_s_apply2_vect => amg_s_apply2_vect procedure, pass(prec) :: psb_s_apply2_vect => amg_s_apply2_vect
procedure, pass(prec) :: psb_s_apply1_vect => amg_s_apply1_vect procedure, pass(prec) :: psb_s_apply1_vect => amg_s_apply1_vect
@ -119,6 +120,7 @@ module amg_s_prec_type
procedure, pass(prec) :: cmp_complexity => amg_s_cmp_compl procedure, pass(prec) :: cmp_complexity => amg_s_cmp_compl
procedure, pass(prec) :: get_avg_cr => amg_s_get_avg_cr procedure, pass(prec) :: get_avg_cr => amg_s_get_avg_cr
procedure, pass(prec) :: cmp_avg_cr => amg_s_cmp_avg_cr procedure, pass(prec) :: cmp_avg_cr => amg_s_cmp_avg_cr
procedure, pass(prec) :: set_nlevs => amg_s_set_nlevs
procedure, pass(prec) :: get_nlevs => amg_s_get_nlevs procedure, pass(prec) :: get_nlevs => amg_s_get_nlevs
procedure, pass(prec) :: get_nzeros => amg_s_get_nzeros procedure, pass(prec) :: get_nzeros => amg_s_get_nzeros
procedure, pass(prec) :: sizeof => amg_sprec_sizeof procedure, pass(prec) :: sizeof => amg_sprec_sizeof
@ -435,10 +437,19 @@ contains
class(amg_sprec_type), intent(in) :: prec class(amg_sprec_type), intent(in) :: prec
integer(psb_ipk_) :: val integer(psb_ipk_) :: val
val = 0 val = 0
if (allocated(prec%precv)) then !!$ if (allocated(prec%precv)) then
val = size(prec%precv) !!$ val = size(prec%precv)
end if !!$ end if
val = prec%nlevs
!!$ write(0,*) ' NLEVS: ',prec%nlevs, val,size(prec%precv)
end function amg_s_get_nlevs end function amg_s_get_nlevs
subroutine amg_s_set_nlevs(prec,nl)
implicit none
class(amg_sprec_type), intent(inout) :: prec
integer(psb_ipk_) :: nl
prec%nlevs = nl
end subroutine amg_s_set_nlevs
! !
! Function returning the size of the amg_prec_type data structure ! Function returning the size of the amg_prec_type data structure
! in bytes or in number of nonzeros of the operator(s) involved. ! in bytes or in number of nonzeros of the operator(s) involved.
@ -509,7 +520,7 @@ contains
real(psb_spk_) :: num, den, nmin real(psb_spk_) :: num, den, nmin
type(psb_ctxt_type) :: ctxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: il integer(psb_ipk_) :: il,nl
num = -sone num = -sone
den = sone den = sone
@ -519,7 +530,10 @@ contains
num = prec%precv(il)%base_a%get_nzeros() num = prec%precv(il)%base_a%get_nzeros()
if (num >= szero) then if (num >= szero) then
den = num 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()) num = num + max(0,prec%precv(il)%base_a%get_nzeros())
end do end do
end if end if
@ -550,7 +564,6 @@ contains
end function amg_s_get_avg_cr end function amg_s_get_avg_cr
subroutine amg_s_cmp_avg_cr(prec) subroutine amg_s_cmp_avg_cr(prec)
implicit none implicit none
class(amg_sprec_type), intent(inout) :: prec class(amg_sprec_type), intent(inout) :: prec
@ -558,17 +571,14 @@ contains
type(psb_ctxt_type) :: ctxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: il, nl, iam, np integer(psb_ipk_) :: il, nl, iam, np
avgcr = szero avgcr = szero
ctxt = prec%ctxt ctxt = prec%ctxt
call psb_info(ctxt,iam,np) call psb_info(ctxt,iam,np)
if (allocated(prec%precv)) then nl = prec%get_nlevs()
nl = size(prec%precv) do il=2,nl
do il=2,nl avgcr = avgcr + max(szero,prec%precv(il)%szratio)
avgcr = avgcr + max(szero,prec%precv(il)%szratio) end do
end do avgcr = avgcr / (nl-1)
avgcr = avgcr / (nl-1)
end if
call psb_sum(ctxt,avgcr) call psb_sum(ctxt,avgcr)
prec%ag_data%avg_cr = avgcr/np prec%ag_data%avg_cr = avgcr/np
end subroutine amg_s_cmp_avg_cr end subroutine amg_s_cmp_avg_cr
@ -586,9 +596,7 @@ contains
! error code. ! error code.
! !
subroutine amg_sprecfree(p,info) subroutine amg_sprecfree(p,info)
implicit none implicit none
! Arguments ! Arguments
type(amg_sprec_type), intent(inout) :: p type(amg_sprec_type), intent(inout) :: p
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -614,9 +622,7 @@ contains
end subroutine amg_sprecfree end subroutine amg_sprecfree
subroutine amg_s_prec_free(prec,info) subroutine amg_s_prec_free(prec,info)
implicit none implicit none
! Arguments ! Arguments
class(amg_sprec_type), intent(inout) :: prec class(amg_sprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -648,9 +654,7 @@ contains
end subroutine amg_s_prec_free end subroutine amg_s_prec_free
subroutine amg_s_smoothers_free(prec,info) subroutine amg_s_smoothers_free(prec,info)
implicit none implicit none
! Arguments ! Arguments
class(amg_sprec_type), intent(inout) :: prec class(amg_sprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -667,7 +671,7 @@ contains
end if end if
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
do i=1,size(prec%precv) do i=1,prec%get_nlevs()
call prec%precv(i)%free_smoothers(info) call prec%precv(i)%free_smoothers(info)
end do end do
end if end if
@ -681,9 +685,7 @@ contains
end subroutine amg_s_smoothers_free end subroutine amg_s_smoothers_free
subroutine amg_s_hierarchy_free(prec,info) subroutine amg_s_hierarchy_free(prec,info)
implicit none implicit none
! Arguments ! Arguments
class(amg_sprec_type), intent(inout) :: prec class(amg_sprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -710,7 +712,6 @@ contains
end subroutine amg_s_hierarchy_free end subroutine amg_s_hierarchy_free
! !
! Top level methods. ! Top level methods.
! !
@ -775,7 +776,6 @@ contains
end subroutine amg_s_apply1_vect end subroutine amg_s_apply1_vect
subroutine amg_s_apply2v(prec,x,y,desc_data,info,trans,work) subroutine amg_s_apply2v(prec,x,y,desc_data,info,trans,work)
implicit none implicit none
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
@ -840,7 +840,6 @@ contains
subroutine amg_s_dump(prec,info,istart,iend,iproc,prefix,head,& subroutine amg_s_dump(prec,info,istart,iend,iproc,prefix,head,&
& ac,rp,smoother,solver,tprol,& & ac,rp,smoother,solver,tprol,&
& global_num) & global_num)
implicit none implicit none
class(amg_sprec_type), intent(in) :: prec class(amg_sprec_type), intent(in) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -857,7 +856,7 @@ contains
info = 0 info = 0
ctxt = prec%ctxt ctxt = prec%ctxt
call psb_info(ctxt,iam,np) call psb_info(ctxt,iam,np)
iln = size(prec%precv) iln = prec%get_nlevs()
if (present(istart)) then if (present(istart)) then
il1 = max(1,istart) il1 = max(1,istart)
else else
@ -881,7 +880,6 @@ contains
end subroutine amg_s_dump end subroutine amg_s_dump
subroutine amg_s_cnv(prec,info,amold,vmold,imold) subroutine amg_s_cnv(prec,info,amold,vmold,imold)
implicit none implicit none
class(amg_sprec_type), intent(inout) :: prec class(amg_sprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -893,7 +891,7 @@ contains
info = psb_success_ info = psb_success_
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
do i=1,size(prec%precv) do i=1,prec%get_nlevs()
if (info == psb_success_ ) & if (info == psb_success_ ) &
& call prec%precv(i)%cnv(info,amold=amold,vmold=vmold,imold=imold) & call prec%precv(i)%cnv(info,amold=amold,vmold=vmold,imold=imold)
end do end do
@ -902,7 +900,6 @@ contains
end subroutine amg_s_cnv end subroutine amg_s_cnv
subroutine amg_s_clone(prec,precout,info) subroutine amg_s_clone(prec,precout,info)
implicit none implicit none
class(amg_sprec_type), intent(inout) :: prec class(amg_sprec_type), intent(inout) :: prec
class(psb_sprec_type), intent(inout) :: precout class(psb_sprec_type), intent(inout) :: precout
@ -914,7 +911,6 @@ contains
end subroutine amg_s_clone end subroutine amg_s_clone
subroutine amg_s_inner_clone(prec,precout,info) subroutine amg_s_inner_clone(prec,precout,info)
implicit none implicit none
class(amg_sprec_type), intent(inout) :: prec class(amg_sprec_type), intent(inout) :: prec
class(psb_sprec_type), target, intent(inout) :: precout class(psb_sprec_type), target, intent(inout) :: precout
@ -930,8 +926,9 @@ contains
pout%ctxt = prec%ctxt pout%ctxt = prec%ctxt
pout%ag_data = prec%ag_data pout%ag_data = prec%ag_data
pout%outer_sweeps = prec%outer_sweeps pout%outer_sweeps = prec%outer_sweeps
pout%nlevs = prec%nlevs
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
ln = size(prec%precv) ln = prec%get_nlevs()
allocate(pout%precv(ln),stat=info) allocate(pout%precv(ln),stat=info)
if (info /= psb_success_) goto 9999 if (info /= psb_success_) goto 9999
if (ln >= 1) then if (ln >= 1) then
@ -939,6 +936,7 @@ contains
end if end if
do lev=2, ln do lev=2, ln
if (info /= psb_success_) exit if (info /= psb_success_) exit
!!$ write(0,*) 'Inner_clone must be checked and reimplemented! '
call prec%precv(lev)%clone(pout%precv(lev),info) call prec%precv(lev)%clone(pout%precv(lev),info)
if (info == psb_success_) then if (info == psb_success_) then
pout%precv(lev)%base_a => pout%precv(lev)%ac pout%precv(lev)%base_a => pout%precv(lev)%ac
@ -1018,7 +1016,7 @@ contains
if (psb_errstatus_fatal()) then if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999 info = psb_err_internal_error_; goto 9999
end if end if
nlev = size(prec%precv) nlev = prec%get_nlevs()
level = 1 level = 1
do level = 1, nlev do level = 1, nlev
call prec%precv(level)%allocate_wrk(info,vmold=vmold) call prec%precv(level)%allocate_wrk(info,vmold=vmold)
@ -1041,7 +1039,6 @@ contains
subroutine amg_s_free_wrk(prec,info) subroutine amg_s_free_wrk(prec,info)
use psb_base_mod use psb_base_mod
implicit none implicit none
! Arguments ! Arguments
class(amg_sprec_type), intent(inout) :: prec class(amg_sprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -1058,7 +1055,7 @@ contains
end if end if
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
nlev = size(prec%precv) nlev = prec%get_nlevs()
do level = 1, nlev do level = 1, nlev
call prec%precv(level)%free_wrk(info) call prec%precv(level)%free_wrk(info)
end do end do

11
amgprec/amg_z_base_aggregator_mod.f90
Unescape Escape View File

@ -89,7 +89,7 @@ module amg_z_base_aggregator_mod
procedure, pass(ag) :: bld_tprol => amg_z_base_aggregator_build_tprol procedure, pass(ag) :: bld_tprol => amg_z_base_aggregator_build_tprol
procedure, pass(ag) :: mat_bld => amg_z_base_aggregator_mat_bld procedure, pass(ag) :: mat_bld => amg_z_base_aggregator_mat_bld
procedure, pass(ag) :: mat_asb => amg_z_base_aggregator_mat_asb procedure, pass(ag) :: mat_asb => amg_z_base_aggregator_mat_asb
procedure, pass(ag) :: bld_map => amg_z_base_aggregator_bld_map procedure, pass(ag) :: bld_linmap => amg_z_base_aggregator_bld_linmap
procedure, pass(ag) :: update_next => amg_z_base_aggregator_update_next procedure, pass(ag) :: update_next => amg_z_base_aggregator_update_next
procedure, pass(ag) :: clone => amg_z_base_aggregator_clone procedure, pass(ag) :: clone => amg_z_base_aggregator_clone
procedure, pass(ag) :: free => amg_z_base_aggregator_free procedure, pass(ag) :: free => amg_z_base_aggregator_free
@ -458,7 +458,7 @@ contains
end subroutine amg_z_base_aggregator_mat_asb end subroutine amg_z_base_aggregator_mat_asb
! !
!> Function bld_map !> Function bld_linmap
!! \memberof amg_z_base_aggregator_type !! \memberof amg_z_base_aggregator_type
!! \brief Build linear map between hierarchy levels !! \brief Build linear map between hierarchy levels
!! !!
@ -473,7 +473,7 @@ contains
!! \param map The output map !! \param map The output map
!! \param info Return code !! \param info Return code
!! !!
subroutine amg_z_base_aggregator_bld_map(ag,desc_a,desc_ac,ilaggr,nlaggr,& subroutine amg_z_base_aggregator_bld_linmap(ag,desc_a,desc_ac,ilaggr,nlaggr,&
& op_restr,op_prol,map,info) & op_restr,op_prol,map,info)
use psb_base_mod use psb_base_mod
implicit none implicit none
@ -484,7 +484,7 @@ contains
type(psb_zlinmap_type), intent(out) :: map type(psb_zlinmap_type), intent(out) :: map
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: err_act integer(psb_ipk_) :: err_act
character(len=20) :: name='z_base_aggregator_bld_map' character(len=20) :: name='z_base_aggregator_bld_linmap'
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -508,7 +508,6 @@ contains
9999 call psb_error_handler(err_act) 9999 call psb_error_handler(err_act)
return return
end subroutine amg_z_base_aggregator_bld_map end subroutine amg_z_base_aggregator_bld_linmap
end module amg_z_base_aggregator_mod end module amg_z_base_aggregator_mod

35
amgprec/amg_z_onelev_mod.f90
Unescape Escape View File

@ -181,7 +181,8 @@ module amg_z_onelev_mod
integer(psb_ipk_) :: idest integer(psb_ipk_) :: idest
integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:) integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:)
contains contains
procedure, pass(rmp) :: clone => z_remap_data_clone procedure, pass(rmp) :: clone => z_remap_data_clone
procedure, pass(rmp) :: move_alloc => z_remap_move_alloc
end type amg_z_remap_data_type end type amg_z_remap_data_type
type amg_z_onelev_type type amg_z_onelev_type
@ -227,7 +228,7 @@ module amg_z_onelev_mod
procedure, pass(lv) :: get_wrksz => z_base_onelev_get_wrksize procedure, pass(lv) :: get_wrksz => z_base_onelev_get_wrksize
procedure, pass(lv) :: allocate_wrk => z_base_onelev_allocate_wrk procedure, pass(lv) :: allocate_wrk => z_base_onelev_allocate_wrk
procedure, pass(lv) :: free_wrk => z_base_onelev_free_wrk procedure, pass(lv) :: free_wrk => z_base_onelev_free_wrk
procedure, nopass :: stringval => amg_stringval procedure, nopass :: stringval => amg_stringval
procedure, pass(lv) :: move_alloc => z_base_onelev_move_alloc procedure, pass(lv) :: move_alloc => z_base_onelev_move_alloc
@ -639,7 +640,7 @@ contains
! Arguments ! Arguments
class(amg_z_onelev_type), target, intent(inout) :: lv class(amg_z_onelev_type), target, intent(inout) :: lv
class(amg_z_onelev_type), target, intent(inout) :: lvout class(amg_z_onelev_type), target, intent(inout) :: lvout
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
info = psb_success_ info = psb_success_
if (allocated(lv%sm)) then if (allocated(lv%sm)) then
@ -705,6 +706,7 @@ contains
if (info == psb_success_) call psb_move_alloc(lv%tprol,b%tprol,info) 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%desc_ac,b%desc_ac,info)
if (info == psb_success_) call psb_move_alloc(lv%linmap,b%linmap,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_a => lv%base_a
b%base_desc => lv%base_desc b%base_desc => lv%base_desc
@ -759,6 +761,7 @@ contains
info = psb_success_ info = psb_success_
nwv = lv%get_wrksz() nwv = lv%get_wrksz()
if (.not.allocated(lv%wrk)) allocate(lv%wrk,stat=info) 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 (info == 0) then
if (lv%remap_data%desc_ac_pre_remap%is_asb()) then if (lv%remap_data%desc_ac_pre_remap%is_asb()) then
! !
@ -806,11 +809,12 @@ contains
info = psb_success_ info = psb_success_
call wk%free(info) 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) allocate(wk%wv(nwv),stat=info)
if (present(desc2).and.(desc%is_valid())) then if (present(desc2).and.(desc%is_valid())) then
write(0,*) 'wrk_alloc D2:',desc2%get_fmt(),desc2%is_asb() !!$ write(0,*) 'wrk_alloc D2:',desc2%get_fmt(),desc2%is_asb()
if (desc2%get_local_cols()>desc%get_local_cols()) then if (desc2%get_local_cols()>desc%get_local_cols()) then
call inner_do_wrk_alloc(wk,nwv,desc2,vmold=vmold) call inner_do_wrk_alloc(wk,nwv,desc2,vmold=vmold)
else else
@ -996,4 +1000,25 @@ contains
call psb_safe_ab_cpy(rmp%nrsrc,remap_out%nrsrc,info) call psb_safe_ab_cpy(rmp%nrsrc,remap_out%nrsrc,info)
end subroutine z_remap_data_clone end subroutine z_remap_data_clone
subroutine z_remap_move_alloc(rmp, remap_out, info)
use psb_base_mod
implicit none
! Arguments
class(amg_z_remap_data_type), target, intent(inout) :: rmp
class(amg_z_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 z_remap_move_alloc
end module amg_z_onelev_mod end module amg_z_onelev_mod

67
amgprec/amg_z_prec_type.f90
Unescape Escape View File

@ -102,6 +102,7 @@ module amg_z_prec_type
! The multilevel hierarchy ! The multilevel hierarchy
! !
type(amg_z_onelev_type), allocatable :: precv(:) type(amg_z_onelev_type), allocatable :: precv(:)
integer(psb_ipk_) :: nlevs
contains contains
procedure, pass(prec) :: psb_z_apply2_vect => amg_z_apply2_vect procedure, pass(prec) :: psb_z_apply2_vect => amg_z_apply2_vect
procedure, pass(prec) :: psb_z_apply1_vect => amg_z_apply1_vect procedure, pass(prec) :: psb_z_apply1_vect => amg_z_apply1_vect
@ -119,6 +120,7 @@ module amg_z_prec_type
procedure, pass(prec) :: cmp_complexity => amg_z_cmp_compl procedure, pass(prec) :: cmp_complexity => amg_z_cmp_compl
procedure, pass(prec) :: get_avg_cr => amg_z_get_avg_cr procedure, pass(prec) :: get_avg_cr => amg_z_get_avg_cr
procedure, pass(prec) :: cmp_avg_cr => amg_z_cmp_avg_cr procedure, pass(prec) :: cmp_avg_cr => amg_z_cmp_avg_cr
procedure, pass(prec) :: set_nlevs => amg_z_set_nlevs
procedure, pass(prec) :: get_nlevs => amg_z_get_nlevs procedure, pass(prec) :: get_nlevs => amg_z_get_nlevs
procedure, pass(prec) :: get_nzeros => amg_z_get_nzeros procedure, pass(prec) :: get_nzeros => amg_z_get_nzeros
procedure, pass(prec) :: sizeof => amg_zprec_sizeof procedure, pass(prec) :: sizeof => amg_zprec_sizeof
@ -435,10 +437,19 @@ contains
class(amg_zprec_type), intent(in) :: prec class(amg_zprec_type), intent(in) :: prec
integer(psb_ipk_) :: val integer(psb_ipk_) :: val
val = 0 val = 0
if (allocated(prec%precv)) then !!$ if (allocated(prec%precv)) then
val = size(prec%precv) !!$ val = size(prec%precv)
end if !!$ end if
val = prec%nlevs
!!$ write(0,*) ' NLEVS: ',prec%nlevs, val,size(prec%precv)
end function amg_z_get_nlevs end function amg_z_get_nlevs
subroutine amg_z_set_nlevs(prec,nl)
implicit none
class(amg_zprec_type), intent(inout) :: prec
integer(psb_ipk_) :: nl
prec%nlevs = nl
end subroutine amg_z_set_nlevs
! !
! Function returning the size of the amg_prec_type data structure ! Function returning the size of the amg_prec_type data structure
! in bytes or in number of nonzeros of the operator(s) involved. ! in bytes or in number of nonzeros of the operator(s) involved.
@ -509,7 +520,7 @@ contains
real(psb_dpk_) :: num, den, nmin real(psb_dpk_) :: num, den, nmin
type(psb_ctxt_type) :: ctxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: il integer(psb_ipk_) :: il,nl
num = -done num = -done
den = done den = done
@ -519,7 +530,10 @@ contains
num = prec%precv(il)%base_a%get_nzeros() num = prec%precv(il)%base_a%get_nzeros()
if (num >= dzero) then if (num >= dzero) then
den = num 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()) num = num + max(0,prec%precv(il)%base_a%get_nzeros())
end do end do
end if end if
@ -550,7 +564,6 @@ contains
end function amg_z_get_avg_cr end function amg_z_get_avg_cr
subroutine amg_z_cmp_avg_cr(prec) subroutine amg_z_cmp_avg_cr(prec)
implicit none implicit none
class(amg_zprec_type), intent(inout) :: prec class(amg_zprec_type), intent(inout) :: prec
@ -558,17 +571,14 @@ contains
type(psb_ctxt_type) :: ctxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: il, nl, iam, np integer(psb_ipk_) :: il, nl, iam, np
avgcr = dzero avgcr = dzero
ctxt = prec%ctxt ctxt = prec%ctxt
call psb_info(ctxt,iam,np) call psb_info(ctxt,iam,np)
if (allocated(prec%precv)) then nl = prec%get_nlevs()
nl = size(prec%precv) do il=2,nl
do il=2,nl avgcr = avgcr + max(dzero,prec%precv(il)%szratio)
avgcr = avgcr + max(dzero,prec%precv(il)%szratio) end do
end do avgcr = avgcr / (nl-1)
avgcr = avgcr / (nl-1)
end if
call psb_sum(ctxt,avgcr) call psb_sum(ctxt,avgcr)
prec%ag_data%avg_cr = avgcr/np prec%ag_data%avg_cr = avgcr/np
end subroutine amg_z_cmp_avg_cr end subroutine amg_z_cmp_avg_cr
@ -586,9 +596,7 @@ contains
! error code. ! error code.
! !
subroutine amg_zprecfree(p,info) subroutine amg_zprecfree(p,info)
implicit none implicit none
! Arguments ! Arguments
type(amg_zprec_type), intent(inout) :: p type(amg_zprec_type), intent(inout) :: p
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -614,9 +622,7 @@ contains
end subroutine amg_zprecfree end subroutine amg_zprecfree
subroutine amg_z_prec_free(prec,info) subroutine amg_z_prec_free(prec,info)
implicit none implicit none
! Arguments ! Arguments
class(amg_zprec_type), intent(inout) :: prec class(amg_zprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -648,9 +654,7 @@ contains
end subroutine amg_z_prec_free end subroutine amg_z_prec_free
subroutine amg_z_smoothers_free(prec,info) subroutine amg_z_smoothers_free(prec,info)
implicit none implicit none
! Arguments ! Arguments
class(amg_zprec_type), intent(inout) :: prec class(amg_zprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -667,7 +671,7 @@ contains
end if end if
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
do i=1,size(prec%precv) do i=1,prec%get_nlevs()
call prec%precv(i)%free_smoothers(info) call prec%precv(i)%free_smoothers(info)
end do end do
end if end if
@ -681,9 +685,7 @@ contains
end subroutine amg_z_smoothers_free end subroutine amg_z_smoothers_free
subroutine amg_z_hierarchy_free(prec,info) subroutine amg_z_hierarchy_free(prec,info)
implicit none implicit none
! Arguments ! Arguments
class(amg_zprec_type), intent(inout) :: prec class(amg_zprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -710,7 +712,6 @@ contains
end subroutine amg_z_hierarchy_free end subroutine amg_z_hierarchy_free
! !
! Top level methods. ! Top level methods.
! !
@ -775,7 +776,6 @@ contains
end subroutine amg_z_apply1_vect end subroutine amg_z_apply1_vect
subroutine amg_z_apply2v(prec,x,y,desc_data,info,trans,work) subroutine amg_z_apply2v(prec,x,y,desc_data,info,trans,work)
implicit none implicit none
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
@ -840,7 +840,6 @@ contains
subroutine amg_z_dump(prec,info,istart,iend,iproc,prefix,head,& subroutine amg_z_dump(prec,info,istart,iend,iproc,prefix,head,&
& ac,rp,smoother,solver,tprol,& & ac,rp,smoother,solver,tprol,&
& global_num) & global_num)
implicit none implicit none
class(amg_zprec_type), intent(in) :: prec class(amg_zprec_type), intent(in) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -857,7 +856,7 @@ contains
info = 0 info = 0
ctxt = prec%ctxt ctxt = prec%ctxt
call psb_info(ctxt,iam,np) call psb_info(ctxt,iam,np)
iln = size(prec%precv) iln = prec%get_nlevs()
if (present(istart)) then if (present(istart)) then
il1 = max(1,istart) il1 = max(1,istart)
else else
@ -881,7 +880,6 @@ contains
end subroutine amg_z_dump end subroutine amg_z_dump
subroutine amg_z_cnv(prec,info,amold,vmold,imold) subroutine amg_z_cnv(prec,info,amold,vmold,imold)
implicit none implicit none
class(amg_zprec_type), intent(inout) :: prec class(amg_zprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -893,7 +891,7 @@ contains
info = psb_success_ info = psb_success_
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
do i=1,size(prec%precv) do i=1,prec%get_nlevs()
if (info == psb_success_ ) & if (info == psb_success_ ) &
& call prec%precv(i)%cnv(info,amold=amold,vmold=vmold,imold=imold) & call prec%precv(i)%cnv(info,amold=amold,vmold=vmold,imold=imold)
end do end do
@ -902,7 +900,6 @@ contains
end subroutine amg_z_cnv end subroutine amg_z_cnv
subroutine amg_z_clone(prec,precout,info) subroutine amg_z_clone(prec,precout,info)
implicit none implicit none
class(amg_zprec_type), intent(inout) :: prec class(amg_zprec_type), intent(inout) :: prec
class(psb_zprec_type), intent(inout) :: precout class(psb_zprec_type), intent(inout) :: precout
@ -914,7 +911,6 @@ contains
end subroutine amg_z_clone end subroutine amg_z_clone
subroutine amg_z_inner_clone(prec,precout,info) subroutine amg_z_inner_clone(prec,precout,info)
implicit none implicit none
class(amg_zprec_type), intent(inout) :: prec class(amg_zprec_type), intent(inout) :: prec
class(psb_zprec_type), target, intent(inout) :: precout class(psb_zprec_type), target, intent(inout) :: precout
@ -930,8 +926,9 @@ contains
pout%ctxt = prec%ctxt pout%ctxt = prec%ctxt
pout%ag_data = prec%ag_data pout%ag_data = prec%ag_data
pout%outer_sweeps = prec%outer_sweeps pout%outer_sweeps = prec%outer_sweeps
pout%nlevs = prec%nlevs
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
ln = size(prec%precv) ln = prec%get_nlevs()
allocate(pout%precv(ln),stat=info) allocate(pout%precv(ln),stat=info)
if (info /= psb_success_) goto 9999 if (info /= psb_success_) goto 9999
if (ln >= 1) then if (ln >= 1) then
@ -939,6 +936,7 @@ contains
end if end if
do lev=2, ln do lev=2, ln
if (info /= psb_success_) exit if (info /= psb_success_) exit
!!$ write(0,*) 'Inner_clone must be checked and reimplemented! '
call prec%precv(lev)%clone(pout%precv(lev),info) call prec%precv(lev)%clone(pout%precv(lev),info)
if (info == psb_success_) then if (info == psb_success_) then
pout%precv(lev)%base_a => pout%precv(lev)%ac pout%precv(lev)%base_a => pout%precv(lev)%ac
@ -1018,7 +1016,7 @@ contains
if (psb_errstatus_fatal()) then if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999 info = psb_err_internal_error_; goto 9999
end if end if
nlev = size(prec%precv) nlev = prec%get_nlevs()
level = 1 level = 1
do level = 1, nlev do level = 1, nlev
call prec%precv(level)%allocate_wrk(info,vmold=vmold) call prec%precv(level)%allocate_wrk(info,vmold=vmold)
@ -1041,7 +1039,6 @@ contains
subroutine amg_z_free_wrk(prec,info) subroutine amg_z_free_wrk(prec,info)
use psb_base_mod use psb_base_mod
implicit none implicit none
! Arguments ! Arguments
class(amg_zprec_type), intent(inout) :: prec class(amg_zprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -1058,7 +1055,7 @@ contains
end if end if
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
nlev = size(prec%precv) nlev = prec%get_nlevs()
do level = 1, nlev do level = 1, nlev
call prec%precv(level)%free_wrk(info) call prec%precv(level)%free_wrk(info)
end do end do

322
amgprec/impl/amg_c_hierarchy_bld.F90
Unescape Escape View File

@ -64,11 +64,9 @@
! Error code. ! Error code.
! !
subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat) subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat)
use psb_base_mod use psb_base_mod
use amg_c_inner_mod use amg_c_inner_mod
use amg_c_prec_mod, amg_protect_name => amg_c_hierarchy_bld use amg_c_prec_mod, amg_protect_name => amg_c_hierarchy_bld
Implicit None Implicit None
! Arguments ! Arguments
@ -82,7 +80,7 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat)
type(psb_ctxt_type) :: ctxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: me,np integer(psb_ipk_) :: me,np
integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz,& integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz,&
& nplevs, mxplevs & nplevs, mxplevs, level
integer(psb_lpk_) :: iaggsize, casize, mncsize, mncszpp integer(psb_lpk_) :: iaggsize, casize, mncsize, mncszpp
real(psb_spk_) :: mnaggratio, sizeratio, athresh, aomega real(psb_spk_) :: mnaggratio, sizeratio, athresh, aomega
class(amg_c_base_smoother_type), allocatable :: coarse_sm, med_sm, & class(amg_c_base_smoother_type), allocatable :: coarse_sm, med_sm, &
@ -98,6 +96,7 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat)
character(len=40) :: ch_err character(len=40) :: ch_err
integer(psb_ipk_), save :: idx_bldtp=-1, idx_matasb=-1 integer(psb_ipk_), save :: idx_bldtp=-1, idx_matasb=-1
logical, parameter :: do_timings=.false. logical, parameter :: do_timings=.false.
logical :: stop_hierarchy_loop
type(psb_ctxt_type) :: lctxt type(psb_ctxt_type) :: lctxt
integer(psb_ipk_) :: lme,lnp integer(psb_ipk_) :: lme,lnp
@ -141,7 +140,7 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat)
mnaggratio = prec%ag_data%min_cr_ratio mnaggratio = prec%ag_data%min_cr_ratio
mncsize = prec%ag_data%min_coarse_size mncsize = prec%ag_data%min_coarse_size
mncszpp = prec%ag_data%min_coarse_size_per_process 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,iszv)
call psb_bcast(ctxt,mncsize) call psb_bcast(ctxt,mncsize)
call psb_bcast(ctxt,mncszpp) call psb_bcast(ctxt,mncszpp)
@ -167,7 +166,7 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat)
call psb_errpush(info,name,a_err='Inconsistent min_cr_ratio') call psb_errpush(info,name,a_err='Inconsistent min_cr_ratio')
goto 9999 goto 9999
end if end if
if (iszv /= size(prec%precv)) then if (iszv /= prec%get_nlevs()) then
info=psb_err_internal_error_ info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Inconsistent size of precv') call psb_errpush(info,name,a_err='Inconsistent size of precv')
goto 9999 goto 9999
@ -182,6 +181,7 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat)
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
endif endif
if (iszv == 1) then if (iszv == 1) then
! !
! This is OK, since it may be called by the user even if there ! This is OK, since it may be called by the user even if there
@ -229,7 +229,6 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat)
casize = mncsize casize = mncsize
end if end if
prec%ag_data%target_coarse_size = casize prec%ag_data%target_coarse_size = casize
nplevs = max(itwo,mxplevs) nplevs = max(itwo,mxplevs)
! !
@ -242,7 +241,7 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat)
goto 9999 goto 9999
end if end if
! !
! First set desired number of levels ! First set desired number of levels if different from default.
! !
if (iszv /= nplevs) then if (iszv /= nplevs) then
allocate(tprecv(nplevs),stat=info) allocate(tprecv(nplevs),stat=info)
@ -288,7 +287,8 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat)
call prec%precv(i)%free(info) call prec%precv(i)%free(info)
end do end do
call move_alloc(tprecv,prec%precv) call move_alloc(tprecv,prec%precv)
iszv = size(prec%precv) call prec%set_nlevs(nplevs)
iszv = prec%get_nlevs()
end if end if
! !
@ -304,16 +304,23 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat)
call psb_cd_renum_block(desc_a,prec%precv(1)%desc_ac,info) call psb_cd_renum_block(desc_a,prec%precv(1)%desc_ac,info)
prec%precv(1)%base_desc => prec%precv(1)%desc_ac prec%precv(1)%base_desc => prec%precv(1)%desc_ac
!
! Main build loop
!
newsz = 0 newsz = 0
stop_hierarchy_loop = .false.
array_build_loop: do i=2, iszv array_build_loop: do i=2, iszv
! !
! Check on the iprcparm contents: they should be the same ! Check on the iprcparm contents: they should be the same
! on all processes. ! on all processes.
! !
call psb_bcast(ctxt,prec%precv(i)%parms) call psb_bcast(ctxt,prec%precv(i)%parms)
!
! Get current context: might have performed remapping
!
lctxt = prec%precv(i-1)%base_desc%get_ctxt() lctxt = prec%precv(i-1)%base_desc%get_ctxt()
call psb_info(lctxt,lme,lnp) call psb_info(lctxt,lme,lnp)
write(0,*) 'Check at level',i,lme,lnp !!$ write(0,*) 'Check at level',i,lme,lnp
! !
! Sanity checks on the parameters ! Sanity checks on the parameters
! !
@ -329,8 +336,8 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat)
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Calling mlprcbld at level ',i & 'Calling mlprcbld at level ',i
! !
! Build the mapping between levels i-1 and i and the matrix ! Build the tentative mapping between levels i-1 and i
! at level i ! and the matrix at level i
! !
if (do_timings) call psb_tic(idx_bldtp) if (do_timings) call psb_tic(idx_bldtp)
if (info == psb_success_)& if (info == psb_success_)&
@ -352,62 +359,26 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat)
! Save op_prol just in case ! Save op_prol just in case
! !
call op_prol%clone(prec%precv(i)%tprol,info) call op_prol%clone(prec%precv(i)%tprol,info)
! !
! Check for early termination of aggregation loop. ! Check for early termination of aggregation loop.
! !
iaggsize = sum(nlaggr) if (i == 2) then
call amg_c_hierarchy_bld_cmp_newsz(i,iszv,&
sizeratio = iaggsize & desc_a%get_global_rows(),&
if (i==2) then & nlaggr,casize,mnaggratio,sizeratio,newsz)
sizeratio = desc_a%get_global_rows()/sizeratio
else else
sizeratio = sum(prec%precv(i-1)%linmap%naggr)/sizeratio call amg_c_hierarchy_bld_cmp_newsz(i,iszv,&
& sum(prec%precv(i-1)%linmap%naggr),&
& nlaggr,casize,mnaggratio,sizeratio,newsz)
end if end if
prec%precv(i)%szratio = sizeratio 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) call psb_bcast(ctxt,newsz)
!
! Handle reallocation, if needed, and then mat_asb to polish off the
! construction
!
if (newsz > 0) then if (newsz > 0) then
! !
! This is awkward, we are saving the aggregation parms, for the sake ! This is awkward, we are saving the aggregation parms, for the sake
@ -441,156 +412,112 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat)
& a_err=ch_err) & a_err=ch_err)
goto 9999 goto 9999
endif endif
if (amg_get_do_remap().and.(i>=2)) then !!$ write(0,*) ' Early exit of array_build_loop',i,iszv,info,&
block level = newsz
type(psb_ctxt_type) :: lctxt stop_hierarchy_loop = .true.
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
else 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(& if (info == psb_success_) call prec%precv(i)%mat_asb(&
& prec%precv(i-1)%base_a,prec%precv(i-1)%base_desc,& & prec%precv(i-1)%base_a,prec%precv(i-1)%base_desc,&
& ilaggr,nlaggr,op_prol,info) & ilaggr,nlaggr,op_prol,info)
if (do_timings) call psb_toc(idx_matasb) if (do_timings) call psb_toc(idx_matasb)
level = i
end if 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,& ! Do we want to remap onto a smaller subset of processes?
& a_err=ch_err) ! Will need a more sophisticated policy
goto 9999 !
endif block
if (i<iszv) call prec%precv(i)%update_aggr(prec%precv(i+1),info) type(psb_ctxt_type) :: lctxt
integer(psb_ipk_) :: lme,lnp
if (amg_get_do_remap().and.(i>=2)) then lctxt = prec%precv(level)%desc_ac%get_ctxt()
block call psb_info(lctxt,lme,lnp)
type(psb_ctxt_type) :: lctxt if (amg_c_policy_do_remap(lctxt,level,sum(nlaggr))) then
integer(psb_ipk_) :: lme,lnp !!$ write(0,*) ' Context on remapping ',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 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%desc_ac%clone(rmp%desc_ac_pre_remap,info)
call lv%ac%clone(rmp%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,& 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) & 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(),' ',& !!$ 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%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%naggr(:) = rmp%naggr(:)
lv%linmap%p_desc_V => rmp%desc_ac_pre_remap lv%linmap%p_desc_V => rmp%desc_ac_pre_remap
lv%base_a => lv%ac lv%base_a => lv%ac
lv%base_desc => lv%desc_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 associate
end if end if
end block !!$ write(0,*) 'Second Check on out remapping ',level,&
end if !!$ & prec%precv(level)%remap_data%desc_ac_pre_remap%is_asb(), newsz
write(0,*) ' End of array_build_loop',i,iszv,info end if
end do array_build_loop end block
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 if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,& write(ch_err,'(a,i7)') 'Mat asb fail @ level ',i
& a_err='prec reallocation') call psb_errpush(psb_err_internal_error_,name,&
& a_err=ch_err)
goto 9999 goto 9999
endif endif
do i=1,newsz if (stop_hierarchy_loop) then
call prec%precv(i)%move_alloc(tprecv(i),info) exit array_build_loop
end do else
do i=newsz+1, iszv if (i<iszv) call prec%precv(i)%update_aggr(prec%precv(i+1),info)
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 end if
do i=2, iszv end do array_build_loop
prec%precv(i)%base_a => prec%precv(i)%ac
prec%precv(i)%base_desc => prec%precv(i)%desc_ac !!$ write(0,*) ' Done array_build_loop',iszv,newsz,info,psb_errstatus_fatal()
! This is needed when the linmap object has been built
! reusing the base_desc descriptor through a pointer. if (newsz>0) then
! With PSBLAS 4 we will have a better solution !!$ do i=2,newsz
if (associated(prec%precv(i)%linmap%p_desc_U)) & !!$ write(0,*) me,'Newsz Out of array_build_loop ',i,':',&
& prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc !!$ & prec%precv(i)%remap_data%desc_ac_pre_remap%is_asb()
if (associated(prec%precv(i)%linmap%p_desc_V))& !!$ end do
& prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc !!$ write(0,*) 'Calling set_nlevs ',newsz
end do 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
end if end if
write(0,*) ' Done reallocating precv',iszv,newsz,info,psb_errstatus_fatal() iszv = prec%get_nlevs()
call psb_barrier(ctxt) call psb_barrier(ctxt)
!write(0,*) 'Should we remap? ' !!$ write(0,*) ' Done reallocating precv',iszv,newsz,info
if (.false.) then !!$
if (amg_get_do_remap().and.(np>=4)) then !!$ do i=2, iszv
!!$ write(0,*) 'Going for remapping ' !!$ write(0,*) me,'At end of hierarchy_bld level',i,':',&
if (.true.) then !!$ & prec%precv(i)%remap_data%desc_ac_pre_remap%is_asb()
associate(lv=>prec%precv(iszv), rmp => prec%precv(iszv)%remap_data) !!$ end do
call lv%desc_ac%clone(rmp%desc_ac_pre_remap,info) call psb_barrier(ctxt)
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 if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Internal hierarchy build' ) & a_err='Internal hierarchy build' )
goto 9999 goto 9999
endif 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_complexity()
call prec%cmp_avg_cr() call prec%cmp_avg_cr()
@ -730,4 +657,49 @@ contains
return return
end subroutine restore_smoothers end subroutine restore_smoothers
#endif #endif
function amg_c_policy_do_remap(ctxt,level,aggsize) result(res)
logical :: res
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: level
integer(psb_lpk_) :: aggsize
res = amg_get_do_remap().and.(level>=2)
!!$ res = .false.
end function amg_c_policy_do_remap
subroutine amg_c_hierarchy_bld_cmp_newsz(level,iszv,prevsize,&
& nlaggr,casize,mnratio,sizeratio,newsz)
implicit none
integer(psb_ipk_) :: level,iszv,newsz
integer(psb_lpk_) :: nlaggr(:)
integer(psb_lpk_) :: prevsize, casize
real(psb_spk_) :: mnratio, sizeratio
! ==============================
integer(psb_lpk_) :: iaggsize
newsz = 0
iaggsize = sum(nlaggr)
sizeratio = prevsize
sizeratio = sizeratio/iaggsize
!!$ write(0,*) 'From cmp_newsz: ',iaggsize,casize,&
!!$ & sizeratio,mnratio, level
if (iaggsize <= casize) newsz = level
if (level == iszv) newsz = level
if (level>2) then
if (sizeratio < mnratio) then
if (sizeratio > 1) then
newsz = level
else
!
! We are not gaining
!
newsz = level-1
end if
end if
end if
!!$ write(0,*) 'At end of cmp_newsz ',newsz
end subroutine amg_c_hierarchy_bld_cmp_newsz
end subroutine amg_c_hierarchy_bld end subroutine amg_c_hierarchy_bld

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

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

3
amgprec/impl/amg_cfile_prec_descr.f90
Unescape Escape View File

@ -129,7 +129,7 @@ subroutine amg_cfile_prec_descr(prec,info,iout,root, verbosity,prefix)
! ensured by amg_precbld). ! ensured by amg_precbld).
! !
if (me == root_) then if (me == root_) then
nlev = size(prec%precv) nlev = prec%get_nlevs()
do ilev = 1, nlev do ilev = 1, nlev
if (.not.allocated(prec%precv(ilev)%sm)) then if (.not.allocated(prec%precv(ilev)%sm)) then
info = 3111 info = 3111
@ -141,6 +141,7 @@ subroutine amg_cfile_prec_descr(prec,info,iout,root, verbosity,prefix)
write(iout_,*) write(iout_,*)
write(iout_,'(a,1x,a)') trim(prefix_),'Preconditioner description' write(iout_,'(a,1x,a)') trim(prefix_),'Preconditioner description'
write(iout_,*) 'At level :',1,' we have ',np,' processes'
if (nlev == 1) then if (nlev == 1) then
! !

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

@ -244,10 +244,10 @@ subroutine amg_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(p%precv) & ' Entry ', p%get_nlevs()
trans_ = psb_toupper(trans) trans_ = psb_toupper(trans)
nlev = size(p%precv) nlev = p%get_nlevs()
do_alloc_wrk = .not.allocated(p%precv(1)%wrk) do_alloc_wrk = .not.allocated(p%precv(1)%wrk)
@ -382,7 +382,7 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit() debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level() debug_level = psb_get_debug_level()
nlev = size(p%precv) nlev = p%get_nlevs()
if ((level < 1) .or. (level > nlev)) then if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_ml') & a_err='wrong call level to inner_ml')
@ -468,7 +468,7 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit() debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level() debug_level = psb_get_debug_level()
nlev = size(p%precv) nlev = p%get_nlevs()
if ((level < 1) .or. (level > nlev)) then if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_add') & a_err='wrong call level to inner_add')
@ -497,7 +497,8 @@ contains
if (allocated(p%precv(level)%sm2a)) then if (allocated(p%precv(level)%sm2a)) then
call psb_geaxpby(cone,vx2l,czero,vy2l,base_desc,info) call psb_geaxpby(cone,vx2l,czero,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 do k=1, sweeps
call p%precv(level)%sm%apply(cone,& call p%precv(level)%sm%apply(cone,&
& vy2l,czero,vty,& & vy2l,czero,vty,&
@ -593,7 +594,7 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit() debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level() debug_level = psb_get_debug_level()
nlev = size(p%precv) nlev = p%get_nlevs()
if ((level < 1) .or. (level > nlev)) then if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_mult') & a_err='wrong call level to inner_mult')
@ -604,7 +605,7 @@ contains
if(debug_level > 1) then if(debug_level > 1) then
write(debug_unit,*) me,' inner_mult at level ',level write(debug_unit,*) me,' inner_mult at level ',level
end if end if
write(debug_unit,*) me,' inner_mult at level (1):',level,np !!$ write(debug_unit,*) me,' inner_mult at level (1):',level,np
sweeps_post = p%precv(level)%parms%sweeps_post sweeps_post = p%precv(level)%parms%sweeps_post
sweeps_pre = p%precv(level)%parms%sweeps_pre sweeps_pre = p%precv(level)%parms%sweeps_pre
pre = ((sweeps_pre>0).and.(trans=='N')).or.((sweeps_post>0).and.(trans/='N')) pre = ((sweeps_pre>0).and.(trans=='N')).or.((sweeps_post>0).and.(trans/='N'))
@ -614,8 +615,8 @@ contains
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv) & wv => p%precv(level)%wrk%wv)
write(0,*) 'Inner mult at level (2):',level,' :',me,np,':',& !!$ write(0,*) 'Inner mult at level (2):',level,' :',me,np,':',&
& size(p%precv(level)%wrk%wv), allocated(p%precv(level)%wrk%wv) !!$ & size(p%precv(level)%wrk%wv), allocated(p%precv(level)%wrk%wv)
if (me >=0) then if (me >=0) then
if (level < nlev) then if (level < nlev) then
@ -831,7 +832,7 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit() debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level() debug_level = psb_get_debug_level()
nlev = size(p%precv) nlev = p%get_nlevs()
if ((level < 1) .or. (level > nlev)) then if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_add') & a_err='wrong call level to inner_add')

5
amgprec/impl/amg_cprecinit.F90
Unescape Escape View File

@ -214,9 +214,7 @@ subroutine amg_cprecinit(ctxt,prec,ptype,info)
allocate(amg_c_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) allocate(amg_c_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
call prec%precv(ilev_)%default() call prec%precv(ilev_)%default()
case ('ML') case ('ML')
nlev_ = prec%ag_data%max_levs nlev_ = prec%ag_data%max_levs
ilev_ = 1 ilev_ = 1
allocate(prec%precv(nlev_),stat=info) allocate(prec%precv(nlev_),stat=info)
@ -228,6 +226,8 @@ subroutine amg_cprecinit(ctxt,prec,ptype,info)
do ilev_ = 1, nlev_ do ilev_ = 1, nlev_
call prec%precv(ilev_)%default() call prec%precv(ilev_)%default()
end do end do
call prec%set_nlevs(nlev_)
call prec%set('ML_CYCLE','VCYCLE',info) call prec%set('ML_CYCLE','VCYCLE',info)
call prec%set('SMOOTHER_TYPE','FBGS',info) call prec%set('SMOOTHER_TYPE','FBGS',info)
#if defined(AMG_HAVE_MUMPS) #if defined(AMG_HAVE_MUMPS)
@ -241,7 +241,6 @@ subroutine amg_cprecinit(ctxt,prec,ptype,info)
write(psb_err_unit,*) name,& write(psb_err_unit,*) name,&
&': Warning: Unknown preconditioner type request "',ptype,'"' &': Warning: Unknown preconditioner type request "',ptype,'"'
info = psb_err_pivot_too_small_ info = psb_err_pivot_too_small_
end select end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)

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 mxplevs = p%ag_data%max_levs
mnaggratio = p%ag_data%min_cr_ratio mnaggratio = p%ag_data%min_cr_ratio
casize = p%ag_data%min_coarse_size casize = p%ag_data%min_coarse_size
iszv = p%get_nlevs() iszv = size(p%precv)
nprolv = size(prolv) nprolv = size(prolv)
nrestrv = size(restrv) nrestrv = size(restrv)
call psb_bcast(ctxt,iszv) 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') call psb_errpush(info,name,a_err='Inconsistent min_cr_ratio')
goto 9999 goto 9999
end if end if
if (iszv /= p%get_nlevs()) then if (iszv /= size(p%precv)) then
info=psb_err_internal_error_ info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Inconsistent size of precv') call psb_errpush(info,name,a_err='Inconsistent size of precv')
goto 9999 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) call p%precv(i)%free(info)
end do end do
call move_alloc(tprecv,p%precv) call move_alloc(tprecv,p%precv)
iszv = p%get_nlevs() iszv = size(p%precv)
end if end if
! !
! Finest level first; remember to fix base_a and base_desc ! 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 goto 9999
endif endif
iszv = p%get_nlevs() iszv = size(p%precv)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&

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

@ -64,11 +64,9 @@
! Error code. ! Error code.
! !
subroutine amg_d_hierarchy_bld(a,desc_a,prec,info,cpymat) subroutine amg_d_hierarchy_bld(a,desc_a,prec,info,cpymat)
use psb_base_mod use psb_base_mod
use amg_d_inner_mod use amg_d_inner_mod
use amg_d_prec_mod, amg_protect_name => amg_d_hierarchy_bld use amg_d_prec_mod, amg_protect_name => amg_d_hierarchy_bld
Implicit None Implicit None
! Arguments ! Arguments
@ -306,7 +304,6 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info,cpymat)
call psb_cd_renum_block(desc_a,prec%precv(1)%desc_ac,info) call psb_cd_renum_block(desc_a,prec%precv(1)%desc_ac,info)
prec%precv(1)%base_desc => prec%precv(1)%desc_ac prec%precv(1)%base_desc => prec%precv(1)%desc_ac
! !
! Main build loop ! Main build loop
! !
@ -340,7 +337,7 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info,cpymat)
& 'Calling mlprcbld at level ',i & 'Calling mlprcbld at level ',i
! !
! Build the tentative mapping between levels i-1 and i ! Build the tentative mapping between levels i-1 and i
! and the matrixat level i ! and the matrix at level i
! !
if (do_timings) call psb_tic(idx_bldtp) if (do_timings) call psb_tic(idx_bldtp)
if (info == psb_success_)& if (info == psb_success_)&
@ -512,7 +509,6 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info,cpymat)
!!$ end do !!$ end do
call psb_barrier(ctxt) call psb_barrier(ctxt)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Internal hierarchy build' ) & a_err='Internal hierarchy build' )

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 end if
iszv = prec%get_nlevs() iszv = size(prec%precv)
do i=2, iszv do i=2, iszv
call prec%precv(i-1)%base_a%cp_to(acsr) call prec%precv(i-1)%base_a%cp_to(acsr)

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

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_,*)
write(iout_,'(a,1x,a)') trim(prefix_),'Preconditioner memory usage' write(iout_,'(a,1x,a)') trim(prefix_),'Preconditioner memory usage'
end if end if
nlev = prec%get_nlevs() nlev = size(prec%precv)
do ilev=1,nlev do ilev=1,nlev
call prec%precv(ilev)%memory_use(ilev,nlev,ilmin,info, & call prec%precv(ilev)%memory_use(ilev,nlev,ilmin,info, &
& iout=iout_,verbosity=verbosity_,prefix=trim(prefix_),global=global) & iout=iout_,verbosity=verbosity_,prefix=trim(prefix_),global=global)

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

@ -651,20 +651,19 @@ contains
! !
if (pre) then if (pre) then
call psb_geaxpby(done,vx2l,& call psb_geaxpby(done,vx2l,&
& dzero,vty,& & dzero,vty,&
& base_desc,info) & base_desc,info)
if (info == psb_success_) call psb_spmm(-done,base_a,& if (info == psb_success_) call psb_spmm(-done,base_a,&
& vy2l,done,vty,& & vy2l,done,vty,&
& base_desc,info,work=work,trans=trans) & base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue') & a_err='Error during residue')
goto 9999 goto 9999
end if end if
!!$ write(0,*) me,' Map_rstr from ', level,' to ',level+1
call p%precv(level+1)%map_rstr(done,vty,& call p%precv(level+1)%map_rstr(done,vty,&
& dzero,p%precv(level+1)%wrk%vx2l,& & dzero,p%precv(level+1)%wrk%vx2l,&
& info,work=work,vtx=wv(1)) & info,work=work,vtx=wv(1))

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

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

@ -226,9 +226,7 @@ subroutine amg_dprecinit(ctxt,prec,ptype,info)
allocate(amg_d_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) allocate(amg_d_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
call prec%precv(ilev_)%default() call prec%precv(ilev_)%default()
case ('ML') case ('ML')
nlev_ = prec%ag_data%max_levs nlev_ = prec%ag_data%max_levs
ilev_ = 1 ilev_ = 1
allocate(prec%precv(nlev_),stat=info) allocate(prec%precv(nlev_),stat=info)

322
amgprec/impl/amg_s_hierarchy_bld.F90
Unescape Escape View File

@ -64,11 +64,9 @@
! Error code. ! Error code.
! !
subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat) subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat)
use psb_base_mod use psb_base_mod
use amg_s_inner_mod use amg_s_inner_mod
use amg_s_prec_mod, amg_protect_name => amg_s_hierarchy_bld use amg_s_prec_mod, amg_protect_name => amg_s_hierarchy_bld
Implicit None Implicit None
! Arguments ! Arguments
@ -82,7 +80,7 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat)
type(psb_ctxt_type) :: ctxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: me,np integer(psb_ipk_) :: me,np
integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz,& integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz,&
& nplevs, mxplevs & nplevs, mxplevs, level
integer(psb_lpk_) :: iaggsize, casize, mncsize, mncszpp integer(psb_lpk_) :: iaggsize, casize, mncsize, mncszpp
real(psb_spk_) :: mnaggratio, sizeratio, athresh, aomega real(psb_spk_) :: mnaggratio, sizeratio, athresh, aomega
class(amg_s_base_smoother_type), allocatable :: coarse_sm, med_sm, & class(amg_s_base_smoother_type), allocatable :: coarse_sm, med_sm, &
@ -98,6 +96,7 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat)
character(len=40) :: ch_err character(len=40) :: ch_err
integer(psb_ipk_), save :: idx_bldtp=-1, idx_matasb=-1 integer(psb_ipk_), save :: idx_bldtp=-1, idx_matasb=-1
logical, parameter :: do_timings=.false. logical, parameter :: do_timings=.false.
logical :: stop_hierarchy_loop
type(psb_ctxt_type) :: lctxt type(psb_ctxt_type) :: lctxt
integer(psb_ipk_) :: lme,lnp integer(psb_ipk_) :: lme,lnp
@ -141,7 +140,7 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat)
mnaggratio = prec%ag_data%min_cr_ratio mnaggratio = prec%ag_data%min_cr_ratio
mncsize = prec%ag_data%min_coarse_size mncsize = prec%ag_data%min_coarse_size
mncszpp = prec%ag_data%min_coarse_size_per_process 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,iszv)
call psb_bcast(ctxt,mncsize) call psb_bcast(ctxt,mncsize)
call psb_bcast(ctxt,mncszpp) call psb_bcast(ctxt,mncszpp)
@ -167,7 +166,7 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat)
call psb_errpush(info,name,a_err='Inconsistent min_cr_ratio') call psb_errpush(info,name,a_err='Inconsistent min_cr_ratio')
goto 9999 goto 9999
end if end if
if (iszv /= size(prec%precv)) then if (iszv /= prec%get_nlevs()) then
info=psb_err_internal_error_ info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Inconsistent size of precv') call psb_errpush(info,name,a_err='Inconsistent size of precv')
goto 9999 goto 9999
@ -182,6 +181,7 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat)
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
endif endif
if (iszv == 1) then if (iszv == 1) then
! !
! This is OK, since it may be called by the user even if there ! This is OK, since it may be called by the user even if there
@ -229,7 +229,6 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat)
casize = mncsize casize = mncsize
end if end if
prec%ag_data%target_coarse_size = casize prec%ag_data%target_coarse_size = casize
nplevs = max(itwo,mxplevs) nplevs = max(itwo,mxplevs)
! !
@ -242,7 +241,7 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat)
goto 9999 goto 9999
end if end if
! !
! First set desired number of levels ! First set desired number of levels if different from default.
! !
if (iszv /= nplevs) then if (iszv /= nplevs) then
allocate(tprecv(nplevs),stat=info) allocate(tprecv(nplevs),stat=info)
@ -288,7 +287,8 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat)
call prec%precv(i)%free(info) call prec%precv(i)%free(info)
end do end do
call move_alloc(tprecv,prec%precv) call move_alloc(tprecv,prec%precv)
iszv = size(prec%precv) call prec%set_nlevs(nplevs)
iszv = prec%get_nlevs()
end if end if
! !
@ -304,16 +304,23 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat)
call psb_cd_renum_block(desc_a,prec%precv(1)%desc_ac,info) call psb_cd_renum_block(desc_a,prec%precv(1)%desc_ac,info)
prec%precv(1)%base_desc => prec%precv(1)%desc_ac prec%precv(1)%base_desc => prec%precv(1)%desc_ac
!
! Main build loop
!
newsz = 0 newsz = 0
stop_hierarchy_loop = .false.
array_build_loop: do i=2, iszv array_build_loop: do i=2, iszv
! !
! Check on the iprcparm contents: they should be the same ! Check on the iprcparm contents: they should be the same
! on all processes. ! on all processes.
! !
call psb_bcast(ctxt,prec%precv(i)%parms) call psb_bcast(ctxt,prec%precv(i)%parms)
!
! Get current context: might have performed remapping
!
lctxt = prec%precv(i-1)%base_desc%get_ctxt() lctxt = prec%precv(i-1)%base_desc%get_ctxt()
call psb_info(lctxt,lme,lnp) call psb_info(lctxt,lme,lnp)
write(0,*) 'Check at level',i,lme,lnp !!$ write(0,*) 'Check at level',i,lme,lnp
! !
! Sanity checks on the parameters ! Sanity checks on the parameters
! !
@ -329,8 +336,8 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat)
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Calling mlprcbld at level ',i & 'Calling mlprcbld at level ',i
! !
! Build the mapping between levels i-1 and i and the matrix ! Build the tentative mapping between levels i-1 and i
! at level i ! and the matrix at level i
! !
if (do_timings) call psb_tic(idx_bldtp) if (do_timings) call psb_tic(idx_bldtp)
if (info == psb_success_)& if (info == psb_success_)&
@ -352,62 +359,26 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat)
! Save op_prol just in case ! Save op_prol just in case
! !
call op_prol%clone(prec%precv(i)%tprol,info) call op_prol%clone(prec%precv(i)%tprol,info)
! !
! Check for early termination of aggregation loop. ! Check for early termination of aggregation loop.
! !
iaggsize = sum(nlaggr) if (i == 2) then
call amg_s_hierarchy_bld_cmp_newsz(i,iszv,&
sizeratio = iaggsize & desc_a%get_global_rows(),&
if (i==2) then & nlaggr,casize,mnaggratio,sizeratio,newsz)
sizeratio = desc_a%get_global_rows()/sizeratio
else else
sizeratio = sum(prec%precv(i-1)%linmap%naggr)/sizeratio call amg_s_hierarchy_bld_cmp_newsz(i,iszv,&
& sum(prec%precv(i-1)%linmap%naggr),&
& nlaggr,casize,mnaggratio,sizeratio,newsz)
end if end if
prec%precv(i)%szratio = sizeratio 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) call psb_bcast(ctxt,newsz)
!
! Handle reallocation, if needed, and then mat_asb to polish off the
! construction
!
if (newsz > 0) then if (newsz > 0) then
! !
! This is awkward, we are saving the aggregation parms, for the sake ! This is awkward, we are saving the aggregation parms, for the sake
@ -441,156 +412,112 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat)
& a_err=ch_err) & a_err=ch_err)
goto 9999 goto 9999
endif endif
if (amg_get_do_remap().and.(i>=2)) then !!$ write(0,*) ' Early exit of array_build_loop',i,iszv,info,&
block level = newsz
type(psb_ctxt_type) :: lctxt stop_hierarchy_loop = .true.
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
else 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(& if (info == psb_success_) call prec%precv(i)%mat_asb(&
& prec%precv(i-1)%base_a,prec%precv(i-1)%base_desc,& & prec%precv(i-1)%base_a,prec%precv(i-1)%base_desc,&
& ilaggr,nlaggr,op_prol,info) & ilaggr,nlaggr,op_prol,info)
if (do_timings) call psb_toc(idx_matasb) if (do_timings) call psb_toc(idx_matasb)
level = i
end if 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,& ! Do we want to remap onto a smaller subset of processes?
& a_err=ch_err) ! Will need a more sophisticated policy
goto 9999 !
endif block
if (i<iszv) call prec%precv(i)%update_aggr(prec%precv(i+1),info) type(psb_ctxt_type) :: lctxt
integer(psb_ipk_) :: lme,lnp
if (amg_get_do_remap().and.(i>=2)) then lctxt = prec%precv(level)%desc_ac%get_ctxt()
block call psb_info(lctxt,lme,lnp)
type(psb_ctxt_type) :: lctxt if (amg_s_policy_do_remap(lctxt,level,sum(nlaggr))) then
integer(psb_ipk_) :: lme,lnp !!$ write(0,*) ' Context on remapping ',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 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%desc_ac%clone(rmp%desc_ac_pre_remap,info)
call lv%ac%clone(rmp%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,& 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) & 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(),' ',& !!$ 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%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%naggr(:) = rmp%naggr(:)
lv%linmap%p_desc_V => rmp%desc_ac_pre_remap lv%linmap%p_desc_V => rmp%desc_ac_pre_remap
lv%base_a => lv%ac lv%base_a => lv%ac
lv%base_desc => lv%desc_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 associate
end if end if
end block !!$ write(0,*) 'Second Check on out remapping ',level,&
end if !!$ & prec%precv(level)%remap_data%desc_ac_pre_remap%is_asb(), newsz
write(0,*) ' End of array_build_loop',i,iszv,info end if
end do array_build_loop end block
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 if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,& write(ch_err,'(a,i7)') 'Mat asb fail @ level ',i
& a_err='prec reallocation') call psb_errpush(psb_err_internal_error_,name,&
& a_err=ch_err)
goto 9999 goto 9999
endif endif
do i=1,newsz if (stop_hierarchy_loop) then
call prec%precv(i)%move_alloc(tprecv(i),info) exit array_build_loop
end do else
do i=newsz+1, iszv if (i<iszv) call prec%precv(i)%update_aggr(prec%precv(i+1),info)
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 end if
do i=2, iszv end do array_build_loop
prec%precv(i)%base_a => prec%precv(i)%ac
prec%precv(i)%base_desc => prec%precv(i)%desc_ac !!$ write(0,*) ' Done array_build_loop',iszv,newsz,info,psb_errstatus_fatal()
! This is needed when the linmap object has been built
! reusing the base_desc descriptor through a pointer. if (newsz>0) then
! With PSBLAS 4 we will have a better solution !!$ do i=2,newsz
if (associated(prec%precv(i)%linmap%p_desc_U)) & !!$ write(0,*) me,'Newsz Out of array_build_loop ',i,':',&
& prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc !!$ & prec%precv(i)%remap_data%desc_ac_pre_remap%is_asb()
if (associated(prec%precv(i)%linmap%p_desc_V))& !!$ end do
& prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc !!$ write(0,*) 'Calling set_nlevs ',newsz
end do 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
end if end if
write(0,*) ' Done reallocating precv',iszv,newsz,info,psb_errstatus_fatal() iszv = prec%get_nlevs()
call psb_barrier(ctxt) call psb_barrier(ctxt)
!write(0,*) 'Should we remap? ' !!$ write(0,*) ' Done reallocating precv',iszv,newsz,info
if (.false.) then !!$
if (amg_get_do_remap().and.(np>=4)) then !!$ do i=2, iszv
!!$ write(0,*) 'Going for remapping ' !!$ write(0,*) me,'At end of hierarchy_bld level',i,':',&
if (.true.) then !!$ & prec%precv(i)%remap_data%desc_ac_pre_remap%is_asb()
associate(lv=>prec%precv(iszv), rmp => prec%precv(iszv)%remap_data) !!$ end do
call lv%desc_ac%clone(rmp%desc_ac_pre_remap,info) call psb_barrier(ctxt)
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 if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Internal hierarchy build' ) & a_err='Internal hierarchy build' )
goto 9999 goto 9999
endif 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_complexity()
call prec%cmp_avg_cr() call prec%cmp_avg_cr()
@ -730,4 +657,49 @@ contains
return return
end subroutine restore_smoothers end subroutine restore_smoothers
#endif #endif
function amg_s_policy_do_remap(ctxt,level,aggsize) result(res)
logical :: res
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: level
integer(psb_lpk_) :: aggsize
res = amg_get_do_remap().and.(level>=2)
!!$ res = .false.
end function amg_s_policy_do_remap
subroutine amg_s_hierarchy_bld_cmp_newsz(level,iszv,prevsize,&
& nlaggr,casize,mnratio,sizeratio,newsz)
implicit none
integer(psb_ipk_) :: level,iszv,newsz
integer(psb_lpk_) :: nlaggr(:)
integer(psb_lpk_) :: prevsize, casize
real(psb_spk_) :: mnratio, sizeratio
! ==============================
integer(psb_lpk_) :: iaggsize
newsz = 0
iaggsize = sum(nlaggr)
sizeratio = prevsize
sizeratio = sizeratio/iaggsize
!!$ write(0,*) 'From cmp_newsz: ',iaggsize,casize,&
!!$ & sizeratio,mnratio, level
if (iaggsize <= casize) newsz = level
if (level == iszv) newsz = level
if (level>2) then
if (sizeratio < mnratio) then
if (sizeratio > 1) then
newsz = level
else
!
! We are not gaining
!
newsz = level-1
end if
end if
end if
!!$ write(0,*) 'At end of cmp_newsz ',newsz
end subroutine amg_s_hierarchy_bld_cmp_newsz
end subroutine amg_s_hierarchy_bld end subroutine amg_s_hierarchy_bld

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

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

3
amgprec/impl/amg_sfile_prec_descr.f90
Unescape Escape View File

@ -129,7 +129,7 @@ subroutine amg_sfile_prec_descr(prec,info,iout,root, verbosity,prefix)
! ensured by amg_precbld). ! ensured by amg_precbld).
! !
if (me == root_) then if (me == root_) then
nlev = size(prec%precv) nlev = prec%get_nlevs()
do ilev = 1, nlev do ilev = 1, nlev
if (.not.allocated(prec%precv(ilev)%sm)) then if (.not.allocated(prec%precv(ilev)%sm)) then
info = 3111 info = 3111
@ -141,6 +141,7 @@ subroutine amg_sfile_prec_descr(prec,info,iout,root, verbosity,prefix)
write(iout_,*) write(iout_,*)
write(iout_,'(a,1x,a)') trim(prefix_),'Preconditioner description' write(iout_,'(a,1x,a)') trim(prefix_),'Preconditioner description'
write(iout_,*) 'At level :',1,' we have ',np,' processes'
if (nlev == 1) then if (nlev == 1) then
! !

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

@ -244,10 +244,10 @@ subroutine amg_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(p%precv) & ' Entry ', p%get_nlevs()
trans_ = psb_toupper(trans) trans_ = psb_toupper(trans)
nlev = size(p%precv) nlev = p%get_nlevs()
do_alloc_wrk = .not.allocated(p%precv(1)%wrk) do_alloc_wrk = .not.allocated(p%precv(1)%wrk)
@ -382,7 +382,7 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit() debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level() debug_level = psb_get_debug_level()
nlev = size(p%precv) nlev = p%get_nlevs()
if ((level < 1) .or. (level > nlev)) then if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_ml') & a_err='wrong call level to inner_ml')
@ -468,7 +468,7 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit() debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level() debug_level = psb_get_debug_level()
nlev = size(p%precv) nlev = p%get_nlevs()
if ((level < 1) .or. (level > nlev)) then if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_add') & a_err='wrong call level to inner_add')
@ -497,7 +497,8 @@ contains
if (allocated(p%precv(level)%sm2a)) then if (allocated(p%precv(level)%sm2a)) then
call psb_geaxpby(sone,vx2l,szero,vy2l,base_desc,info) call psb_geaxpby(sone,vx2l,szero,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 do k=1, sweeps
call p%precv(level)%sm%apply(sone,& call p%precv(level)%sm%apply(sone,&
& vy2l,szero,vty,& & vy2l,szero,vty,&
@ -593,7 +594,7 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit() debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level() debug_level = psb_get_debug_level()
nlev = size(p%precv) nlev = p%get_nlevs()
if ((level < 1) .or. (level > nlev)) then if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_mult') & a_err='wrong call level to inner_mult')
@ -604,7 +605,7 @@ contains
if(debug_level > 1) then if(debug_level > 1) then
write(debug_unit,*) me,' inner_mult at level ',level write(debug_unit,*) me,' inner_mult at level ',level
end if end if
write(debug_unit,*) me,' inner_mult at level (1):',level,np !!$ write(debug_unit,*) me,' inner_mult at level (1):',level,np
sweeps_post = p%precv(level)%parms%sweeps_post sweeps_post = p%precv(level)%parms%sweeps_post
sweeps_pre = p%precv(level)%parms%sweeps_pre sweeps_pre = p%precv(level)%parms%sweeps_pre
pre = ((sweeps_pre>0).and.(trans=='N')).or.((sweeps_post>0).and.(trans/='N')) pre = ((sweeps_pre>0).and.(trans=='N')).or.((sweeps_post>0).and.(trans/='N'))
@ -614,8 +615,8 @@ contains
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv) & wv => p%precv(level)%wrk%wv)
write(0,*) 'Inner mult at level (2):',level,' :',me,np,':',& !!$ write(0,*) 'Inner mult at level (2):',level,' :',me,np,':',&
& size(p%precv(level)%wrk%wv), allocated(p%precv(level)%wrk%wv) !!$ & size(p%precv(level)%wrk%wv), allocated(p%precv(level)%wrk%wv)
if (me >=0) then if (me >=0) then
if (level < nlev) then if (level < nlev) then
@ -831,7 +832,7 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit() debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level() debug_level = psb_get_debug_level()
nlev = size(p%precv) nlev = p%get_nlevs()
if ((level < 1) .or. (level > nlev)) then if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_add') & a_err='wrong call level to inner_add')

5
amgprec/impl/amg_sprecinit.F90
Unescape Escape View File

@ -223,9 +223,7 @@ subroutine amg_sprecinit(ctxt,prec,ptype,info)
allocate(amg_s_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) allocate(amg_s_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
call prec%precv(ilev_)%default() call prec%precv(ilev_)%default()
case ('ML') case ('ML')
nlev_ = prec%ag_data%max_levs nlev_ = prec%ag_data%max_levs
ilev_ = 1 ilev_ = 1
allocate(prec%precv(nlev_),stat=info) allocate(prec%precv(nlev_),stat=info)
@ -237,6 +235,8 @@ subroutine amg_sprecinit(ctxt,prec,ptype,info)
do ilev_ = 1, nlev_ do ilev_ = 1, nlev_
call prec%precv(ilev_)%default() call prec%precv(ilev_)%default()
end do end do
call prec%set_nlevs(nlev_)
call prec%set('ML_CYCLE','VCYCLE',info) call prec%set('ML_CYCLE','VCYCLE',info)
call prec%set('SMOOTHER_TYPE','FBGS',info) call prec%set('SMOOTHER_TYPE','FBGS',info)
#if defined(AMG_HAVE_MUMPS) #if defined(AMG_HAVE_MUMPS)
@ -250,7 +250,6 @@ subroutine amg_sprecinit(ctxt,prec,ptype,info)
write(psb_err_unit,*) name,& write(psb_err_unit,*) name,&
&': Warning: Unknown preconditioner type request "',ptype,'"' &': Warning: Unknown preconditioner type request "',ptype,'"'
info = psb_err_pivot_too_small_ info = psb_err_pivot_too_small_
end select end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)

322
amgprec/impl/amg_z_hierarchy_bld.F90
Unescape Escape View File

@ -64,11 +64,9 @@
! Error code. ! Error code.
! !
subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat) subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat)
use psb_base_mod use psb_base_mod
use amg_z_inner_mod use amg_z_inner_mod
use amg_z_prec_mod, amg_protect_name => amg_z_hierarchy_bld use amg_z_prec_mod, amg_protect_name => amg_z_hierarchy_bld
Implicit None Implicit None
! Arguments ! Arguments
@ -82,7 +80,7 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat)
type(psb_ctxt_type) :: ctxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: me,np integer(psb_ipk_) :: me,np
integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz,& integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz,&
& nplevs, mxplevs & nplevs, mxplevs, level
integer(psb_lpk_) :: iaggsize, casize, mncsize, mncszpp integer(psb_lpk_) :: iaggsize, casize, mncsize, mncszpp
real(psb_dpk_) :: mnaggratio, sizeratio, athresh, aomega real(psb_dpk_) :: mnaggratio, sizeratio, athresh, aomega
class(amg_z_base_smoother_type), allocatable :: coarse_sm, med_sm, & class(amg_z_base_smoother_type), allocatable :: coarse_sm, med_sm, &
@ -98,6 +96,7 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat)
character(len=40) :: ch_err character(len=40) :: ch_err
integer(psb_ipk_), save :: idx_bldtp=-1, idx_matasb=-1 integer(psb_ipk_), save :: idx_bldtp=-1, idx_matasb=-1
logical, parameter :: do_timings=.false. logical, parameter :: do_timings=.false.
logical :: stop_hierarchy_loop
type(psb_ctxt_type) :: lctxt type(psb_ctxt_type) :: lctxt
integer(psb_ipk_) :: lme,lnp integer(psb_ipk_) :: lme,lnp
@ -141,7 +140,7 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat)
mnaggratio = prec%ag_data%min_cr_ratio mnaggratio = prec%ag_data%min_cr_ratio
mncsize = prec%ag_data%min_coarse_size mncsize = prec%ag_data%min_coarse_size
mncszpp = prec%ag_data%min_coarse_size_per_process 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,iszv)
call psb_bcast(ctxt,mncsize) call psb_bcast(ctxt,mncsize)
call psb_bcast(ctxt,mncszpp) call psb_bcast(ctxt,mncszpp)
@ -167,7 +166,7 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat)
call psb_errpush(info,name,a_err='Inconsistent min_cr_ratio') call psb_errpush(info,name,a_err='Inconsistent min_cr_ratio')
goto 9999 goto 9999
end if end if
if (iszv /= size(prec%precv)) then if (iszv /= prec%get_nlevs()) then
info=psb_err_internal_error_ info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Inconsistent size of precv') call psb_errpush(info,name,a_err='Inconsistent size of precv')
goto 9999 goto 9999
@ -182,6 +181,7 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat)
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
endif endif
if (iszv == 1) then if (iszv == 1) then
! !
! This is OK, since it may be called by the user even if there ! This is OK, since it may be called by the user even if there
@ -229,7 +229,6 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat)
casize = mncsize casize = mncsize
end if end if
prec%ag_data%target_coarse_size = casize prec%ag_data%target_coarse_size = casize
nplevs = max(itwo,mxplevs) nplevs = max(itwo,mxplevs)
! !
@ -242,7 +241,7 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat)
goto 9999 goto 9999
end if end if
! !
! First set desired number of levels ! First set desired number of levels if different from default.
! !
if (iszv /= nplevs) then if (iszv /= nplevs) then
allocate(tprecv(nplevs),stat=info) allocate(tprecv(nplevs),stat=info)
@ -288,7 +287,8 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat)
call prec%precv(i)%free(info) call prec%precv(i)%free(info)
end do end do
call move_alloc(tprecv,prec%precv) call move_alloc(tprecv,prec%precv)
iszv = size(prec%precv) call prec%set_nlevs(nplevs)
iszv = prec%get_nlevs()
end if end if
! !
@ -304,16 +304,23 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat)
call psb_cd_renum_block(desc_a,prec%precv(1)%desc_ac,info) call psb_cd_renum_block(desc_a,prec%precv(1)%desc_ac,info)
prec%precv(1)%base_desc => prec%precv(1)%desc_ac prec%precv(1)%base_desc => prec%precv(1)%desc_ac
!
! Main build loop
!
newsz = 0 newsz = 0
stop_hierarchy_loop = .false.
array_build_loop: do i=2, iszv array_build_loop: do i=2, iszv
! !
! Check on the iprcparm contents: they should be the same ! Check on the iprcparm contents: they should be the same
! on all processes. ! on all processes.
! !
call psb_bcast(ctxt,prec%precv(i)%parms) call psb_bcast(ctxt,prec%precv(i)%parms)
!
! Get current context: might have performed remapping
!
lctxt = prec%precv(i-1)%base_desc%get_ctxt() lctxt = prec%precv(i-1)%base_desc%get_ctxt()
call psb_info(lctxt,lme,lnp) call psb_info(lctxt,lme,lnp)
write(0,*) 'Check at level',i,lme,lnp !!$ write(0,*) 'Check at level',i,lme,lnp
! !
! Sanity checks on the parameters ! Sanity checks on the parameters
! !
@ -329,8 +336,8 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat)
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Calling mlprcbld at level ',i & 'Calling mlprcbld at level ',i
! !
! Build the mapping between levels i-1 and i and the matrix ! Build the tentative mapping between levels i-1 and i
! at level i ! and the matrix at level i
! !
if (do_timings) call psb_tic(idx_bldtp) if (do_timings) call psb_tic(idx_bldtp)
if (info == psb_success_)& if (info == psb_success_)&
@ -352,62 +359,26 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat)
! Save op_prol just in case ! Save op_prol just in case
! !
call op_prol%clone(prec%precv(i)%tprol,info) call op_prol%clone(prec%precv(i)%tprol,info)
! !
! Check for early termination of aggregation loop. ! Check for early termination of aggregation loop.
! !
iaggsize = sum(nlaggr) if (i == 2) then
call amg_z_hierarchy_bld_cmp_newsz(i,iszv,&
sizeratio = iaggsize & desc_a%get_global_rows(),&
if (i==2) then & nlaggr,casize,mnaggratio,sizeratio,newsz)
sizeratio = desc_a%get_global_rows()/sizeratio
else else
sizeratio = sum(prec%precv(i-1)%linmap%naggr)/sizeratio call amg_z_hierarchy_bld_cmp_newsz(i,iszv,&
& sum(prec%precv(i-1)%linmap%naggr),&
& nlaggr,casize,mnaggratio,sizeratio,newsz)
end if end if
prec%precv(i)%szratio = sizeratio 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) call psb_bcast(ctxt,newsz)
!
! Handle reallocation, if needed, and then mat_asb to polish off the
! construction
!
if (newsz > 0) then if (newsz > 0) then
! !
! This is awkward, we are saving the aggregation parms, for the sake ! This is awkward, we are saving the aggregation parms, for the sake
@ -441,156 +412,112 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat)
& a_err=ch_err) & a_err=ch_err)
goto 9999 goto 9999
endif endif
if (amg_get_do_remap().and.(i>=2)) then !!$ write(0,*) ' Early exit of array_build_loop',i,iszv,info,&
block level = newsz
type(psb_ctxt_type) :: lctxt stop_hierarchy_loop = .true.
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
else 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(& if (info == psb_success_) call prec%precv(i)%mat_asb(&
& prec%precv(i-1)%base_a,prec%precv(i-1)%base_desc,& & prec%precv(i-1)%base_a,prec%precv(i-1)%base_desc,&
& ilaggr,nlaggr,op_prol,info) & ilaggr,nlaggr,op_prol,info)
if (do_timings) call psb_toc(idx_matasb) if (do_timings) call psb_toc(idx_matasb)
level = i
end if 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,& ! Do we want to remap onto a smaller subset of processes?
& a_err=ch_err) ! Will need a more sophisticated policy
goto 9999 !
endif block
if (i<iszv) call prec%precv(i)%update_aggr(prec%precv(i+1),info) type(psb_ctxt_type) :: lctxt
integer(psb_ipk_) :: lme,lnp
if (amg_get_do_remap().and.(i>=2)) then lctxt = prec%precv(level)%desc_ac%get_ctxt()
block call psb_info(lctxt,lme,lnp)
type(psb_ctxt_type) :: lctxt if (amg_z_policy_do_remap(lctxt,level,sum(nlaggr))) then
integer(psb_ipk_) :: lme,lnp !!$ write(0,*) ' Context on remapping ',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 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%desc_ac%clone(rmp%desc_ac_pre_remap,info)
call lv%ac%clone(rmp%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,& 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) & 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(),' ',& !!$ 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%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%naggr(:) = rmp%naggr(:)
lv%linmap%p_desc_V => rmp%desc_ac_pre_remap lv%linmap%p_desc_V => rmp%desc_ac_pre_remap
lv%base_a => lv%ac lv%base_a => lv%ac
lv%base_desc => lv%desc_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 associate
end if end if
end block !!$ write(0,*) 'Second Check on out remapping ',level,&
end if !!$ & prec%precv(level)%remap_data%desc_ac_pre_remap%is_asb(), newsz
write(0,*) ' End of array_build_loop',i,iszv,info end if
end do array_build_loop end block
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 if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,& write(ch_err,'(a,i7)') 'Mat asb fail @ level ',i
& a_err='prec reallocation') call psb_errpush(psb_err_internal_error_,name,&
& a_err=ch_err)
goto 9999 goto 9999
endif endif
do i=1,newsz if (stop_hierarchy_loop) then
call prec%precv(i)%move_alloc(tprecv(i),info) exit array_build_loop
end do else
do i=newsz+1, iszv if (i<iszv) call prec%precv(i)%update_aggr(prec%precv(i+1),info)
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 end if
do i=2, iszv end do array_build_loop
prec%precv(i)%base_a => prec%precv(i)%ac
prec%precv(i)%base_desc => prec%precv(i)%desc_ac !!$ write(0,*) ' Done array_build_loop',iszv,newsz,info,psb_errstatus_fatal()
! This is needed when the linmap object has been built
! reusing the base_desc descriptor through a pointer. if (newsz>0) then
! With PSBLAS 4 we will have a better solution !!$ do i=2,newsz
if (associated(prec%precv(i)%linmap%p_desc_U)) & !!$ write(0,*) me,'Newsz Out of array_build_loop ',i,':',&
& prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc !!$ & prec%precv(i)%remap_data%desc_ac_pre_remap%is_asb()
if (associated(prec%precv(i)%linmap%p_desc_V))& !!$ end do
& prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc !!$ write(0,*) 'Calling set_nlevs ',newsz
end do 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
end if end if
write(0,*) ' Done reallocating precv',iszv,newsz,info,psb_errstatus_fatal() iszv = prec%get_nlevs()
call psb_barrier(ctxt) call psb_barrier(ctxt)
!write(0,*) 'Should we remap? ' !!$ write(0,*) ' Done reallocating precv',iszv,newsz,info
if (.false.) then !!$
if (amg_get_do_remap().and.(np>=4)) then !!$ do i=2, iszv
!!$ write(0,*) 'Going for remapping ' !!$ write(0,*) me,'At end of hierarchy_bld level',i,':',&
if (.true.) then !!$ & prec%precv(i)%remap_data%desc_ac_pre_remap%is_asb()
associate(lv=>prec%precv(iszv), rmp => prec%precv(iszv)%remap_data) !!$ end do
call lv%desc_ac%clone(rmp%desc_ac_pre_remap,info) call psb_barrier(ctxt)
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 if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Internal hierarchy build' ) & a_err='Internal hierarchy build' )
goto 9999 goto 9999
endif 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_complexity()
call prec%cmp_avg_cr() call prec%cmp_avg_cr()
@ -730,4 +657,49 @@ contains
return return
end subroutine restore_smoothers end subroutine restore_smoothers
#endif #endif
function amg_z_policy_do_remap(ctxt,level,aggsize) result(res)
logical :: res
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: level
integer(psb_lpk_) :: aggsize
res = amg_get_do_remap().and.(level>=2)
!!$ res = .false.
end function amg_z_policy_do_remap
subroutine amg_z_hierarchy_bld_cmp_newsz(level,iszv,prevsize,&
& nlaggr,casize,mnratio,sizeratio,newsz)
implicit none
integer(psb_ipk_) :: level,iszv,newsz
integer(psb_lpk_) :: nlaggr(:)
integer(psb_lpk_) :: prevsize, casize
real(psb_dpk_) :: mnratio, sizeratio
! ==============================
integer(psb_lpk_) :: iaggsize
newsz = 0
iaggsize = sum(nlaggr)
sizeratio = prevsize
sizeratio = sizeratio/iaggsize
!!$ write(0,*) 'From cmp_newsz: ',iaggsize,casize,&
!!$ & sizeratio,mnratio, level
if (iaggsize <= casize) newsz = level
if (level == iszv) newsz = level
if (level>2) then
if (sizeratio < mnratio) then
if (sizeratio > 1) then
newsz = level
else
!
! We are not gaining
!
newsz = level-1
end if
end if
end if
!!$ write(0,*) 'At end of cmp_newsz ',newsz
end subroutine amg_z_hierarchy_bld_cmp_newsz
end subroutine amg_z_hierarchy_bld end subroutine amg_z_hierarchy_bld

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

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

3
amgprec/impl/amg_zfile_prec_descr.f90
Unescape Escape View File

@ -129,7 +129,7 @@ subroutine amg_zfile_prec_descr(prec,info,iout,root, verbosity,prefix)
! ensured by amg_precbld). ! ensured by amg_precbld).
! !
if (me == root_) then if (me == root_) then
nlev = size(prec%precv) nlev = prec%get_nlevs()
do ilev = 1, nlev do ilev = 1, nlev
if (.not.allocated(prec%precv(ilev)%sm)) then if (.not.allocated(prec%precv(ilev)%sm)) then
info = 3111 info = 3111
@ -141,6 +141,7 @@ subroutine amg_zfile_prec_descr(prec,info,iout,root, verbosity,prefix)
write(iout_,*) write(iout_,*)
write(iout_,'(a,1x,a)') trim(prefix_),'Preconditioner description' write(iout_,'(a,1x,a)') trim(prefix_),'Preconditioner description'
write(iout_,*) 'At level :',1,' we have ',np,' processes'
if (nlev == 1) then if (nlev == 1) then
! !

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

@ -244,10 +244,10 @@ subroutine amg_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(p%precv) & ' Entry ', p%get_nlevs()
trans_ = psb_toupper(trans) trans_ = psb_toupper(trans)
nlev = size(p%precv) nlev = p%get_nlevs()
do_alloc_wrk = .not.allocated(p%precv(1)%wrk) do_alloc_wrk = .not.allocated(p%precv(1)%wrk)
@ -382,7 +382,7 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit() debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level() debug_level = psb_get_debug_level()
nlev = size(p%precv) nlev = p%get_nlevs()
if ((level < 1) .or. (level > nlev)) then if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_ml') & a_err='wrong call level to inner_ml')
@ -468,7 +468,7 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit() debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level() debug_level = psb_get_debug_level()
nlev = size(p%precv) nlev = p%get_nlevs()
if ((level < 1) .or. (level > nlev)) then if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_add') & a_err='wrong call level to inner_add')
@ -497,7 +497,8 @@ contains
if (allocated(p%precv(level)%sm2a)) then if (allocated(p%precv(level)%sm2a)) then
call psb_geaxpby(zone,vx2l,zzero,vy2l,base_desc,info) call psb_geaxpby(zone,vx2l,zzero,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 do k=1, sweeps
call p%precv(level)%sm%apply(zone,& call p%precv(level)%sm%apply(zone,&
& vy2l,zzero,vty,& & vy2l,zzero,vty,&
@ -593,7 +594,7 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit() debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level() debug_level = psb_get_debug_level()
nlev = size(p%precv) nlev = p%get_nlevs()
if ((level < 1) .or. (level > nlev)) then if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_mult') & a_err='wrong call level to inner_mult')
@ -604,7 +605,7 @@ contains
if(debug_level > 1) then if(debug_level > 1) then
write(debug_unit,*) me,' inner_mult at level ',level write(debug_unit,*) me,' inner_mult at level ',level
end if end if
write(debug_unit,*) me,' inner_mult at level (1):',level,np !!$ write(debug_unit,*) me,' inner_mult at level (1):',level,np
sweeps_post = p%precv(level)%parms%sweeps_post sweeps_post = p%precv(level)%parms%sweeps_post
sweeps_pre = p%precv(level)%parms%sweeps_pre sweeps_pre = p%precv(level)%parms%sweeps_pre
pre = ((sweeps_pre>0).and.(trans=='N')).or.((sweeps_post>0).and.(trans/='N')) pre = ((sweeps_pre>0).and.(trans=='N')).or.((sweeps_post>0).and.(trans/='N'))
@ -614,8 +615,8 @@ contains
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv) & wv => p%precv(level)%wrk%wv)
write(0,*) 'Inner mult at level (2):',level,' :',me,np,':',& !!$ write(0,*) 'Inner mult at level (2):',level,' :',me,np,':',&
& size(p%precv(level)%wrk%wv), allocated(p%precv(level)%wrk%wv) !!$ & size(p%precv(level)%wrk%wv), allocated(p%precv(level)%wrk%wv)
if (me >=0) then if (me >=0) then
if (level < nlev) then if (level < nlev) then
@ -831,7 +832,7 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit() debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level() debug_level = psb_get_debug_level()
nlev = size(p%precv) nlev = p%get_nlevs()
if ((level < 1) .or. (level > nlev)) then if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_add') & a_err='wrong call level to inner_add')

5
amgprec/impl/amg_zprecinit.F90
Unescape Escape View File

@ -217,9 +217,7 @@ subroutine amg_zprecinit(ctxt,prec,ptype,info)
allocate(amg_z_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) allocate(amg_z_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
call prec%precv(ilev_)%default() call prec%precv(ilev_)%default()
case ('ML') case ('ML')
nlev_ = prec%ag_data%max_levs nlev_ = prec%ag_data%max_levs
ilev_ = 1 ilev_ = 1
allocate(prec%precv(nlev_),stat=info) allocate(prec%precv(nlev_),stat=info)
@ -231,6 +229,8 @@ subroutine amg_zprecinit(ctxt,prec,ptype,info)
do ilev_ = 1, nlev_ do ilev_ = 1, nlev_
call prec%precv(ilev_)%default() call prec%precv(ilev_)%default()
end do end do
call prec%set_nlevs(nlev_)
call prec%set('ML_CYCLE','VCYCLE',info) call prec%set('ML_CYCLE','VCYCLE',info)
call prec%set('SMOOTHER_TYPE','FBGS',info) call prec%set('SMOOTHER_TYPE','FBGS',info)
#if defined(AMG_HAVE_UMF) #if defined(AMG_HAVE_UMF)
@ -246,7 +246,6 @@ subroutine amg_zprecinit(ctxt,prec,ptype,info)
write(psb_err_unit,*) name,& write(psb_err_unit,*) name,&
&': Warning: Unknown preconditioner type request "',ptype,'"' &': Warning: Unknown preconditioner type request "',ptype,'"'
info = psb_err_pivot_too_small_ info = psb_err_pivot_too_small_
end select end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)

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

@ -46,6 +46,7 @@ subroutine amg_c_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
complex(psb_spk_), optional :: work(:) complex(psb_spk_), optional :: work(:)
type(psb_c_vect_type), optional, target, intent(inout) :: vtx,vty type(psb_c_vect_type), optional, target, intent(inout) :: vtx,vty
type(psb_c_vect_type), pointer :: vtx_ type(psb_c_vect_type), pointer :: vtx_
!!$ 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()

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

@ -48,6 +48,7 @@ subroutine amg_c_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
complex(psb_spk_), optional :: work(:) complex(psb_spk_), optional :: work(:)
type(psb_c_vect_type), optional, target, intent(inout) :: vtx,vty type(psb_c_vect_type), optional, target, intent(inout) :: vtx,vty
type(psb_c_vect_type), pointer :: vty_ type(psb_c_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 if (present(vty)) then
@ -64,7 +65,7 @@ subroutine amg_c_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
block block
type(psb_ctxt_type) :: ctxt, rctxt type(psb_ctxt_type) :: ctxt, rctxt
integer(psb_mpk_) :: i,j,ip, idest, nsrc, nrl, kp integer(psb_mpk_) :: i,j,ip, idest, nsrc, nrl, kp
integer(psb_mpk_) :: me, np, rme, rnp integer(psb_mpk_) :: rme, rnp
complex(psb_spk_), allocatable :: rsnd(:), rrcv(:) complex(psb_spk_), allocatable :: rsnd(:), rrcv(:)
type(psb_c_vect_type) :: tv type(psb_c_vect_type) :: tv
@ -72,7 +73,7 @@ subroutine amg_c_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
call psb_info(ctxt,me,np) call psb_info(ctxt,me,np)
rctxt = lv%desc_ac%get_ctxt() rctxt = lv%desc_ac%get_ctxt()
call psb_info(rctxt,rme,rnp) 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 idest = lv%remap_data%idest
associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc) associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc)
!!$ write(0,*) 'Should apply maps, then send data from ',me,' to ',idest !!$ write(0,*) 'Should apply maps, then send data from ',me,' to ',idest
@ -81,15 +82,16 @@ 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() 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_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) call psb_geasb(tv,lv%remap_data%desc_ac_pre_remap,info,mold=vect_u%v)
write(0,*) me,' remap 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) !!$ & psb_errstatus_fatal()
!!$ flush(0)
call psb_barrier(ctxt) call psb_barrier(ctxt)
call lv%linmap%map_U2V(alpha,vect_u,beta,tv,info,& 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() call tv%sync()
!rsnd = tv%get_vect() !rsnd = tv%get_vect()
!call psb_snd(ctxt,rsnd(1:nrl),idest) !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) call psb_snd(ctxt,tv%v%v(1:nrl),idest)
if (rme >=0) then if (rme >=0) then
allocate(rrcv(sum(nrsrc))) allocate(rrcv(sum(nrsrc)))
@ -98,29 +100,28 @@ subroutine amg_c_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
do i = 1,size(isrc) do i = 1,size(isrc)
ip = isrc(i) ip = isrc(i)
nrl = nrsrc(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) call psb_rcv(ctxt,rrcv(kp+1:kp+nrl),ip)
kp = kp + nrl kp = kp + nrl
end do end do
call vect_v%set_vect(rrcv) call vect_v%set_vect(rrcv)
end if end if
end associate end associate
!!$ write(0,*) me, ' Restrictor with remap done ' !!$ write(0,*) me, ' Restrictor with remap done ',psb_errstatus_fatal()
end block end block
else else
! Default transfer ! Default transfer
block block
type(psb_ctxt_type) :: ctxt, rctxt type(psb_ctxt_type) :: ctxt, rctxt
integer(psb_mpk_) :: me, np
ctxt = lv%linmap%p_desc_U%get_ctxt() ctxt = lv%linmap%p_desc_U%get_ctxt()
call psb_info(ctxt,me,np) call psb_info(ctxt,me,np)
write(0,*) me,' map_rstr calling U2V: ',me,np !!$ write(0,*) me,' map_rstr calling U2V: ',me,np
call lv%linmap%map_U2V(alpha,vect_u,beta,vect_v,info,& call lv%linmap%map_U2V(alpha,vect_u,beta,vect_v,info,&
& work=work,vtx=vtx,vty=vty_) & work=work,vtx=vtx,vty=vty_)
end block end block
end if end if
!!$ write(0,*) me, 'End of restriction ',info,psb_errstatus_fatal()
end subroutine amg_c_base_onelev_map_rstr_v end subroutine amg_c_base_onelev_map_rstr_v
subroutine amg_c_base_onelev_map_rstr_a(lv,alpha,u,beta,v,info,work) subroutine amg_c_base_onelev_map_rstr_a(lv,alpha,u,beta,v,info,work)

2
amgprec/impl/level/amg_c_base_onelev_mat_asb.f90
Unescape Escape View File

@ -168,7 +168,7 @@ subroutine amg_c_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,t_prol,info)
if (do_timings) call psb_tic(idx_mapbld) if (do_timings) call psb_tic(idx_mapbld)
if (info == psb_success_) call lv%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) if (info == psb_success_) call lv%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info == psb_success_) call lv%aggr%bld_map(desc_a, lv%desc_ac,& if (info == psb_success_) call lv%aggr%bld_linmap(desc_a, lv%desc_ac,&
& ilaggr,nlaggr,op_restr,op_prol,lv%linmap,info) & ilaggr,nlaggr,op_restr,op_prol,lv%linmap,info)
if (do_timings) call psb_toc(idx_mapbld) if (do_timings) call psb_toc(idx_mapbld)
if(info /= psb_success_) then if(info /= psb_success_) then

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

@ -35,6 +35,7 @@
! POSSIBILITY OF SUCH DAMAGE. ! POSSIBILITY OF SUCH DAMAGE.
! !
! !
subroutine amg_d_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& subroutine amg_d_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
& work,vtx,vty) & work,vtx,vty)
use psb_base_mod use psb_base_mod
@ -47,9 +48,9 @@ subroutine amg_d_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
real(psb_dpk_), optional :: work(:) real(psb_dpk_), optional :: work(:)
type(psb_d_vect_type), optional, target, intent(inout) :: vtx,vty type(psb_d_vect_type), optional, target, intent(inout) :: vtx,vty
type(psb_d_vect_type), pointer :: vty_ type(psb_d_vect_type), pointer :: vty_
integer(psb_mpk_) :: me, np 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 if (present(vty)) then
vty_ => vty vty_ => vty
else else
@ -60,7 +61,7 @@ subroutine amg_d_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
! !
! Remap has happened, deal with it ! Remap has happened, deal with it
! !
!!$ write(0,*) 'Level map_rstr with remapping ' !!$ write(0,*) 'Remap handling not implemented yet '
block block
type(psb_ctxt_type) :: ctxt, rctxt type(psb_ctxt_type) :: ctxt, rctxt
integer(psb_mpk_) :: i,j,ip, idest, nsrc, nrl, kp integer(psb_mpk_) :: i,j,ip, idest, nsrc, nrl, kp
@ -87,7 +88,7 @@ subroutine amg_d_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
call psb_barrier(ctxt) call psb_barrier(ctxt)
call lv%linmap%map_U2V(alpha,vect_u,beta,tv,info,& 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() call tv%sync()
!rsnd = tv%get_vect() !rsnd = tv%get_vect()
!call psb_snd(ctxt,rsnd(1:nrl),idest) !call psb_snd(ctxt,rsnd(1:nrl),idest)
!!$ write(0,*) me,' map_rstr sending ',me,idest,psb_errstatus_fatal() !!$ write(0,*) me,' map_rstr sending ',me,idest,psb_errstatus_fatal()

2
amgprec/impl/level/amg_d_base_onelev_mat_asb.f90
Unescape Escape View File

@ -168,7 +168,7 @@ subroutine amg_d_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,t_prol,info)
if (do_timings) call psb_tic(idx_mapbld) if (do_timings) call psb_tic(idx_mapbld)
if (info == psb_success_) call lv%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) if (info == psb_success_) call lv%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info == psb_success_) call lv%aggr%bld_map(desc_a, lv%desc_ac,& if (info == psb_success_) call lv%aggr%bld_linmap(desc_a, lv%desc_ac,&
& ilaggr,nlaggr,op_restr,op_prol,lv%linmap,info) & ilaggr,nlaggr,op_restr,op_prol,lv%linmap,info)
if (do_timings) call psb_toc(idx_mapbld) if (do_timings) call psb_toc(idx_mapbld)
if(info /= psb_success_) then if(info /= psb_success_) then

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

@ -48,6 +48,7 @@ subroutine amg_s_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
real(psb_spk_), optional :: work(:) real(psb_spk_), optional :: work(:)
type(psb_s_vect_type), optional, target, intent(inout) :: vtx,vty type(psb_s_vect_type), optional, target, intent(inout) :: vtx,vty
type(psb_s_vect_type), pointer :: vty_ type(psb_s_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 if (present(vty)) then
@ -64,7 +65,7 @@ subroutine amg_s_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
block block
type(psb_ctxt_type) :: ctxt, rctxt type(psb_ctxt_type) :: ctxt, rctxt
integer(psb_mpk_) :: i,j,ip, idest, nsrc, nrl, kp integer(psb_mpk_) :: i,j,ip, idest, nsrc, nrl, kp
integer(psb_mpk_) :: me, np, rme, rnp integer(psb_mpk_) :: rme, rnp
real(psb_spk_), allocatable :: rsnd(:), rrcv(:) real(psb_spk_), allocatable :: rsnd(:), rrcv(:)
type(psb_s_vect_type) :: tv type(psb_s_vect_type) :: tv
@ -72,7 +73,7 @@ subroutine amg_s_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
call psb_info(ctxt,me,np) call psb_info(ctxt,me,np)
rctxt = lv%desc_ac%get_ctxt() rctxt = lv%desc_ac%get_ctxt()
call psb_info(rctxt,rme,rnp) 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 idest = lv%remap_data%idest
associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc) associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc)
!!$ write(0,*) 'Should apply maps, then send data from ',me,' to ',idest !!$ write(0,*) 'Should apply maps, then send data from ',me,' to ',idest
@ -81,15 +82,16 @@ 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() 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_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) call psb_geasb(tv,lv%remap_data%desc_ac_pre_remap,info,mold=vect_u%v)
write(0,*) me,' remap 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) !!$ & psb_errstatus_fatal()
!!$ flush(0)
call psb_barrier(ctxt) call psb_barrier(ctxt)
call lv%linmap%map_U2V(alpha,vect_u,beta,tv,info,& 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() call tv%sync()
!rsnd = tv%get_vect() !rsnd = tv%get_vect()
!call psb_snd(ctxt,rsnd(1:nrl),idest) !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) call psb_snd(ctxt,tv%v%v(1:nrl),idest)
if (rme >=0) then if (rme >=0) then
allocate(rrcv(sum(nrsrc))) allocate(rrcv(sum(nrsrc)))
@ -98,29 +100,28 @@ subroutine amg_s_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
do i = 1,size(isrc) do i = 1,size(isrc)
ip = isrc(i) ip = isrc(i)
nrl = nrsrc(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) call psb_rcv(ctxt,rrcv(kp+1:kp+nrl),ip)
kp = kp + nrl kp = kp + nrl
end do end do
call vect_v%set_vect(rrcv) call vect_v%set_vect(rrcv)
end if end if
end associate end associate
!!$ write(0,*) me, ' Restrictor with remap done ' !!$ write(0,*) me, ' Restrictor with remap done ',psb_errstatus_fatal()
end block end block
else else
! Default transfer ! Default transfer
block block
type(psb_ctxt_type) :: ctxt, rctxt type(psb_ctxt_type) :: ctxt, rctxt
integer(psb_mpk_) :: me, np
ctxt = lv%linmap%p_desc_U%get_ctxt() ctxt = lv%linmap%p_desc_U%get_ctxt()
call psb_info(ctxt,me,np) call psb_info(ctxt,me,np)
write(0,*) me,' map_rstr calling U2V: ',me,np !!$ write(0,*) me,' map_rstr calling U2V: ',me,np
call lv%linmap%map_U2V(alpha,vect_u,beta,vect_v,info,& call lv%linmap%map_U2V(alpha,vect_u,beta,vect_v,info,&
& work=work,vtx=vtx,vty=vty_) & work=work,vtx=vtx,vty=vty_)
end block end block
end if end if
!!$ write(0,*) me, 'End of restriction ',info,psb_errstatus_fatal()
end subroutine amg_s_base_onelev_map_rstr_v end subroutine amg_s_base_onelev_map_rstr_v
subroutine amg_s_base_onelev_map_rstr_a(lv,alpha,u,beta,v,info,work) subroutine amg_s_base_onelev_map_rstr_a(lv,alpha,u,beta,v,info,work)

2
amgprec/impl/level/amg_s_base_onelev_mat_asb.f90
Unescape Escape View File

@ -168,7 +168,7 @@ subroutine amg_s_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,t_prol,info)
if (do_timings) call psb_tic(idx_mapbld) if (do_timings) call psb_tic(idx_mapbld)
if (info == psb_success_) call lv%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) if (info == psb_success_) call lv%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info == psb_success_) call lv%aggr%bld_map(desc_a, lv%desc_ac,& if (info == psb_success_) call lv%aggr%bld_linmap(desc_a, lv%desc_ac,&
& ilaggr,nlaggr,op_restr,op_prol,lv%linmap,info) & ilaggr,nlaggr,op_restr,op_prol,lv%linmap,info)
if (do_timings) call psb_toc(idx_mapbld) if (do_timings) call psb_toc(idx_mapbld)
if(info /= psb_success_) then if(info /= psb_success_) then

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

@ -48,6 +48,7 @@ subroutine amg_z_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
complex(psb_dpk_), optional :: work(:) complex(psb_dpk_), optional :: work(:)
type(psb_z_vect_type), optional, target, intent(inout) :: vtx,vty type(psb_z_vect_type), optional, target, intent(inout) :: vtx,vty
type(psb_z_vect_type), pointer :: vty_ type(psb_z_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 if (present(vty)) then
@ -64,7 +65,7 @@ subroutine amg_z_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
block block
type(psb_ctxt_type) :: ctxt, rctxt type(psb_ctxt_type) :: ctxt, rctxt
integer(psb_mpk_) :: i,j,ip, idest, nsrc, nrl, kp integer(psb_mpk_) :: i,j,ip, idest, nsrc, nrl, kp
integer(psb_mpk_) :: me, np, rme, rnp integer(psb_mpk_) :: rme, rnp
complex(psb_dpk_), allocatable :: rsnd(:), rrcv(:) complex(psb_dpk_), allocatable :: rsnd(:), rrcv(:)
type(psb_z_vect_type) :: tv type(psb_z_vect_type) :: tv
@ -72,7 +73,7 @@ subroutine amg_z_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
call psb_info(ctxt,me,np) call psb_info(ctxt,me,np)
rctxt = lv%desc_ac%get_ctxt() rctxt = lv%desc_ac%get_ctxt()
call psb_info(rctxt,rme,rnp) 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 idest = lv%remap_data%idest
associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc) associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc)
!!$ write(0,*) 'Should apply maps, then send data from ',me,' to ',idest !!$ write(0,*) 'Should apply maps, then send data from ',me,' to ',idest
@ -81,15 +82,16 @@ 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() 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_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) call psb_geasb(tv,lv%remap_data%desc_ac_pre_remap,info,mold=vect_u%v)
write(0,*) me,' remap 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) !!$ & psb_errstatus_fatal()
!!$ flush(0)
call psb_barrier(ctxt) call psb_barrier(ctxt)
call lv%linmap%map_U2V(alpha,vect_u,beta,tv,info,& 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() call tv%sync()
!rsnd = tv%get_vect() !rsnd = tv%get_vect()
!call psb_snd(ctxt,rsnd(1:nrl),idest) !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) call psb_snd(ctxt,tv%v%v(1:nrl),idest)
if (rme >=0) then if (rme >=0) then
allocate(rrcv(sum(nrsrc))) allocate(rrcv(sum(nrsrc)))
@ -98,29 +100,28 @@ subroutine amg_z_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
do i = 1,size(isrc) do i = 1,size(isrc)
ip = isrc(i) ip = isrc(i)
nrl = nrsrc(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) call psb_rcv(ctxt,rrcv(kp+1:kp+nrl),ip)
kp = kp + nrl kp = kp + nrl
end do end do
call vect_v%set_vect(rrcv) call vect_v%set_vect(rrcv)
end if end if
end associate end associate
!!$ write(0,*) me, ' Restrictor with remap done ' !!$ write(0,*) me, ' Restrictor with remap done ',psb_errstatus_fatal()
end block end block
else else
! Default transfer ! Default transfer
block block
type(psb_ctxt_type) :: ctxt, rctxt type(psb_ctxt_type) :: ctxt, rctxt
integer(psb_mpk_) :: me, np
ctxt = lv%linmap%p_desc_U%get_ctxt() ctxt = lv%linmap%p_desc_U%get_ctxt()
call psb_info(ctxt,me,np) call psb_info(ctxt,me,np)
write(0,*) me,' map_rstr calling U2V: ',me,np !!$ write(0,*) me,' map_rstr calling U2V: ',me,np
call lv%linmap%map_U2V(alpha,vect_u,beta,vect_v,info,& call lv%linmap%map_U2V(alpha,vect_u,beta,vect_v,info,&
& work=work,vtx=vtx,vty=vty_) & work=work,vtx=vtx,vty=vty_)
end block end block
end if end if
!!$ write(0,*) me, 'End of restriction ',info,psb_errstatus_fatal()
end subroutine amg_z_base_onelev_map_rstr_v end subroutine amg_z_base_onelev_map_rstr_v
subroutine amg_z_base_onelev_map_rstr_a(lv,alpha,u,beta,v,info,work) subroutine amg_z_base_onelev_map_rstr_a(lv,alpha,u,beta,v,info,work)

2
amgprec/impl/level/amg_z_base_onelev_mat_asb.f90
Unescape Escape View File

@ -168,7 +168,7 @@ subroutine amg_z_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,t_prol,info)
if (do_timings) call psb_tic(idx_mapbld) if (do_timings) call psb_tic(idx_mapbld)
if (info == psb_success_) call lv%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) if (info == psb_success_) call lv%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info == psb_success_) call lv%aggr%bld_map(desc_a, lv%desc_ac,& if (info == psb_success_) call lv%aggr%bld_linmap(desc_a, lv%desc_ac,&
& ilaggr,nlaggr,op_restr,op_prol,lv%linmap,info) & ilaggr,nlaggr,op_restr,op_prol,lv%linmap,info)
if (do_timings) call psb_toc(idx_mapbld) if (do_timings) call psb_toc(idx_mapbld)
if(info /= psb_success_) then if(info /= psb_success_) then

Write Preview
Loading…
Cancel
Save