From b0d7272f11f717a10dfa75048f56455d24d01ae2 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Fri, 15 Jul 2016 11:22:00 +0000 Subject: [PATCH] mld2p4-2: mlprec/impl/Makefile mlprec/impl/level/mld_c_base_onelev_dump.f90 mlprec/impl/level/mld_d_base_onelev_dump.f90 mlprec/impl/level/mld_s_base_onelev_dump.f90 mlprec/impl/level/mld_z_base_onelev_dump.f90 mlprec/impl/mld_c_extprol_bld.f90 mlprec/impl/mld_s_extprol_bld.f90 mlprec/impl/mld_z_extprol_bld.f90 mlprec/mld_c_onelev_mod.f90 mlprec/mld_c_prec_mod.f90 mlprec/mld_c_prec_type.f90 mlprec/mld_d_onelev_mod.f90 mlprec/mld_d_prec_mod.f90 mlprec/mld_d_prec_type.f90 mlprec/mld_s_onelev_mod.f90 mlprec/mld_s_prec_mod.f90 mlprec/mld_s_prec_type.f90 mlprec/mld_z_onelev_mod.f90 mlprec/mld_z_prec_mod.f90 mlprec/mld_z_prec_type.f90 New option for dump with global numbering. New method to build hierarchy with externally supplied restrictors and prolongators. --- mlprec/impl/Makefile | 6 +- mlprec/impl/level/mld_c_base_onelev_dump.f90 | 56 +- mlprec/impl/level/mld_d_base_onelev_dump.f90 | 56 +- mlprec/impl/level/mld_s_base_onelev_dump.f90 | 56 +- mlprec/impl/level/mld_z_base_onelev_dump.f90 | 56 +- mlprec/impl/mld_c_extprol_bld.f90 | 524 +++++++++++++++++++ mlprec/impl/mld_s_extprol_bld.f90 | 524 +++++++++++++++++++ mlprec/impl/mld_z_extprol_bld.f90 | 524 +++++++++++++++++++ mlprec/mld_c_onelev_mod.f90 | 5 +- mlprec/mld_c_prec_mod.f90 | 38 +- mlprec/mld_c_prec_type.f90 | 6 +- mlprec/mld_d_onelev_mod.f90 | 5 +- mlprec/mld_d_prec_mod.f90 | 7 +- mlprec/mld_d_prec_type.f90 | 6 +- mlprec/mld_s_onelev_mod.f90 | 5 +- mlprec/mld_s_prec_mod.f90 | 38 +- mlprec/mld_s_prec_type.f90 | 6 +- mlprec/mld_z_onelev_mod.f90 | 5 +- mlprec/mld_z_prec_mod.f90 | 38 +- mlprec/mld_z_prec_type.f90 | 6 +- 20 files changed, 1854 insertions(+), 113 deletions(-) create mode 100644 mlprec/impl/mld_c_extprol_bld.f90 create mode 100644 mlprec/impl/mld_s_extprol_bld.f90 create mode 100644 mlprec/impl/mld_z_extprol_bld.f90 diff --git a/mlprec/impl/Makefile b/mlprec/impl/Makefile index 37258583..c9c0e713 100644 --- a/mlprec/impl/Makefile +++ b/mlprec/impl/Makefile @@ -32,19 +32,19 @@ SINNEROBJS= mld_scoarse_bld.o mld_smlprec_bld.o mld_s_bld_mlhier_aggsize.o mld mld_s_ml_prec_bld.o mld_s_hierarchy_bld.o \ mld_silu0_fact.o mld_siluk_fact.o mld_silut_fact.o mld_saggrmap_bld.o \ mld_s_dec_map_bld.o mld_smlprec_aply.o mld_saggrmat_asb.o \ - $(SMPFOBJS) + $(SMPFOBJS) mld_s_extprol_bld.o ZINNEROBJS= mld_zcoarse_bld.o mld_zmlprec_bld.o mld_z_bld_mlhier_aggsize.o mld_z_bld_mlhier_array.o \ mld_z_ml_prec_bld.o mld_z_hierarchy_bld.o \ mld_zilu0_fact.o mld_ziluk_fact.o mld_zilut_fact.o mld_zaggrmap_bld.o \ mld_z_dec_map_bld.o mld_zmlprec_aply.o mld_zaggrmat_asb.o \ - $(ZMPFOBJS) + $(ZMPFOBJS) mld_z_extprol_bld.o CINNEROBJS= mld_ccoarse_bld.o mld_cmlprec_bld.o mld_c_bld_mlhier_aggsize.o mld_c_bld_mlhier_array.o \ mld_c_ml_prec_bld.o mld_c_hierarchy_bld.o \ mld_cilu0_fact.o mld_ciluk_fact.o mld_cilut_fact.o mld_caggrmap_bld.o \ mld_c_dec_map_bld.o mld_cmlprec_aply.o mld_caggrmat_asb.o \ - $(CMPFOBJS) + $(CMPFOBJS) mld_c_extprol_bld.o INNEROBJS= $(SINNEROBJS) $(DINNEROBJS) $(CINNEROBJS) $(ZINNEROBJS) diff --git a/mlprec/impl/level/mld_c_base_onelev_dump.f90 b/mlprec/impl/level/mld_c_base_onelev_dump.f90 index 545071f0..50c78d4f 100644 --- a/mlprec/impl/level/mld_c_base_onelev_dump.f90 +++ b/mlprec/impl/level/mld_c_base_onelev_dump.f90 @@ -36,8 +36,9 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -subroutine mld_c_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,solver) - +subroutine mld_c_base_onelev_dump(lv,level,info,prefix,head,ac,rp,& + & smoother,solver,global_num) + use psb_base_mod use mld_c_onelev_mod, mld_protect_name => mld_c_base_onelev_dump implicit none @@ -45,13 +46,14 @@ subroutine mld_c_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,solve integer(psb_ipk_), intent(in) :: level integer(psb_ipk_), intent(out) :: info character(len=*), intent(in), optional :: prefix, head - logical, optional, intent(in) :: ac, rp, smoother, solver + logical, optional, intent(in) :: ac, rp, smoother, solver,global_num + ! Local variables integer(psb_ipk_) :: i, j, il1, iln, lname, lev integer(psb_ipk_) :: icontxt,iam, np character(len=80) :: prefix_ character(len=120) :: fname ! len should be at least 20 more than - logical :: ac_, rp_ - ! len of prefix_ + logical :: ac_, rp_,global_num_ + integer(psb_ipk_), allocatable :: ivr(:), ivc(:) info = 0 @@ -78,23 +80,47 @@ subroutine mld_c_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,solve else rp_ = .false. end if + if (present(global_num)) then + global_num_ = global_num + else + global_num_ = .false. + end if lname = len_trim(prefix_) fname = trim(prefix_) write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam lname = lname + 5 - if (level >= 2) then - if (ac_) then - write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_ac.mtx' - call lv%ac%print(fname,head=head) + if (global_num_) then + if (level >= 2) then + if (ac_) then + ivr = lv%desc_ac%get_global_indices(owned=.false.) + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_ac.mtx' + call lv%ac%print(fname,head=head,iv=ivr) + end if + if (rp_) then + ivr = lv%map%p_desc_X%get_global_indices(owned=.false.) + ivc = lv%map%p_desc_Y%get_global_indices(owned=.false.) + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_r.mtx' + call lv%map%map_X2Y%print(fname,head=head,ivr=ivc,ivc=ivr) + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_p.mtx' + call lv%map%map_Y2X%print(fname,head=head,ivr=ivr,ivc=ivc) + end if end if - if (rp_) then - write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_r.mtx' - call lv%map%map_X2Y%print(fname,head=head) - write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_p.mtx' - call lv%map%map_Y2X%print(fname,head=head) + else + if (level >= 2) then + if (ac_) then + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_ac.mtx' + call lv%ac%print(fname,head=head) + end if + if (rp_) then + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_r.mtx' + call lv%map%map_X2Y%print(fname,head=head) + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_p.mtx' + call lv%map%map_Y2X%print(fname,head=head) + end if end if end if + if (allocated(lv%sm)) & & call lv%sm%dump(icontxt,level,info,smoother=smoother, & & solver=solver,prefix=prefix) @@ -103,5 +129,5 @@ subroutine mld_c_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,solve call lv%sm2a%dump(icontxt,level,info,smoother=smoother, & & solver=solver,prefix=prefix_) end if - + end subroutine mld_c_base_onelev_dump diff --git a/mlprec/impl/level/mld_d_base_onelev_dump.f90 b/mlprec/impl/level/mld_d_base_onelev_dump.f90 index 7254e535..6dae5045 100644 --- a/mlprec/impl/level/mld_d_base_onelev_dump.f90 +++ b/mlprec/impl/level/mld_d_base_onelev_dump.f90 @@ -36,8 +36,9 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -subroutine mld_d_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,solver) - +subroutine mld_d_base_onelev_dump(lv,level,info,prefix,head,ac,rp,& + & smoother,solver,global_num) + use psb_base_mod use mld_d_onelev_mod, mld_protect_name => mld_d_base_onelev_dump implicit none @@ -45,13 +46,14 @@ subroutine mld_d_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,solve integer(psb_ipk_), intent(in) :: level integer(psb_ipk_), intent(out) :: info character(len=*), intent(in), optional :: prefix, head - logical, optional, intent(in) :: ac, rp, smoother, solver + logical, optional, intent(in) :: ac, rp, smoother, solver,global_num + ! Local variables integer(psb_ipk_) :: i, j, il1, iln, lname, lev integer(psb_ipk_) :: icontxt,iam, np character(len=80) :: prefix_ character(len=120) :: fname ! len should be at least 20 more than - logical :: ac_, rp_ - ! len of prefix_ + logical :: ac_, rp_,global_num_ + integer(psb_ipk_), allocatable :: ivr(:), ivc(:) info = 0 @@ -78,23 +80,47 @@ subroutine mld_d_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,solve else rp_ = .false. end if + if (present(global_num)) then + global_num_ = global_num + else + global_num_ = .false. + end if lname = len_trim(prefix_) fname = trim(prefix_) write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam lname = lname + 5 - if (level >= 2) then - if (ac_) then - write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_ac.mtx' - call lv%ac%print(fname,head=head) + if (global_num_) then + if (level >= 2) then + if (ac_) then + ivr = lv%desc_ac%get_global_indices(owned=.false.) + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_ac.mtx' + call lv%ac%print(fname,head=head,iv=ivr) + end if + if (rp_) then + ivr = lv%map%p_desc_X%get_global_indices(owned=.false.) + ivc = lv%map%p_desc_Y%get_global_indices(owned=.false.) + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_r.mtx' + call lv%map%map_X2Y%print(fname,head=head,ivr=ivc,ivc=ivr) + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_p.mtx' + call lv%map%map_Y2X%print(fname,head=head,ivr=ivr,ivc=ivc) + end if end if - if (rp_) then - write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_r.mtx' - call lv%map%map_X2Y%print(fname,head=head) - write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_p.mtx' - call lv%map%map_Y2X%print(fname,head=head) + else + if (level >= 2) then + if (ac_) then + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_ac.mtx' + call lv%ac%print(fname,head=head) + end if + if (rp_) then + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_r.mtx' + call lv%map%map_X2Y%print(fname,head=head) + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_p.mtx' + call lv%map%map_Y2X%print(fname,head=head) + end if end if end if + if (allocated(lv%sm)) & & call lv%sm%dump(icontxt,level,info,smoother=smoother, & & solver=solver,prefix=prefix) @@ -103,5 +129,5 @@ subroutine mld_d_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,solve call lv%sm2a%dump(icontxt,level,info,smoother=smoother, & & solver=solver,prefix=prefix_) end if - + end subroutine mld_d_base_onelev_dump diff --git a/mlprec/impl/level/mld_s_base_onelev_dump.f90 b/mlprec/impl/level/mld_s_base_onelev_dump.f90 index 9b020017..21a1cdac 100644 --- a/mlprec/impl/level/mld_s_base_onelev_dump.f90 +++ b/mlprec/impl/level/mld_s_base_onelev_dump.f90 @@ -36,8 +36,9 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -subroutine mld_s_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,solver) - +subroutine mld_s_base_onelev_dump(lv,level,info,prefix,head,ac,rp,& + & smoother,solver,global_num) + use psb_base_mod use mld_s_onelev_mod, mld_protect_name => mld_s_base_onelev_dump implicit none @@ -45,13 +46,14 @@ subroutine mld_s_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,solve integer(psb_ipk_), intent(in) :: level integer(psb_ipk_), intent(out) :: info character(len=*), intent(in), optional :: prefix, head - logical, optional, intent(in) :: ac, rp, smoother, solver + logical, optional, intent(in) :: ac, rp, smoother, solver,global_num + ! Local variables integer(psb_ipk_) :: i, j, il1, iln, lname, lev integer(psb_ipk_) :: icontxt,iam, np character(len=80) :: prefix_ character(len=120) :: fname ! len should be at least 20 more than - logical :: ac_, rp_ - ! len of prefix_ + logical :: ac_, rp_,global_num_ + integer(psb_ipk_), allocatable :: ivr(:), ivc(:) info = 0 @@ -78,23 +80,47 @@ subroutine mld_s_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,solve else rp_ = .false. end if + if (present(global_num)) then + global_num_ = global_num + else + global_num_ = .false. + end if lname = len_trim(prefix_) fname = trim(prefix_) write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam lname = lname + 5 - if (level >= 2) then - if (ac_) then - write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_ac.mtx' - call lv%ac%print(fname,head=head) + if (global_num_) then + if (level >= 2) then + if (ac_) then + ivr = lv%desc_ac%get_global_indices(owned=.false.) + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_ac.mtx' + call lv%ac%print(fname,head=head,iv=ivr) + end if + if (rp_) then + ivr = lv%map%p_desc_X%get_global_indices(owned=.false.) + ivc = lv%map%p_desc_Y%get_global_indices(owned=.false.) + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_r.mtx' + call lv%map%map_X2Y%print(fname,head=head,ivr=ivc,ivc=ivr) + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_p.mtx' + call lv%map%map_Y2X%print(fname,head=head,ivr=ivr,ivc=ivc) + end if end if - if (rp_) then - write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_r.mtx' - call lv%map%map_X2Y%print(fname,head=head) - write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_p.mtx' - call lv%map%map_Y2X%print(fname,head=head) + else + if (level >= 2) then + if (ac_) then + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_ac.mtx' + call lv%ac%print(fname,head=head) + end if + if (rp_) then + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_r.mtx' + call lv%map%map_X2Y%print(fname,head=head) + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_p.mtx' + call lv%map%map_Y2X%print(fname,head=head) + end if end if end if + if (allocated(lv%sm)) & & call lv%sm%dump(icontxt,level,info,smoother=smoother, & & solver=solver,prefix=prefix) @@ -103,5 +129,5 @@ subroutine mld_s_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,solve call lv%sm2a%dump(icontxt,level,info,smoother=smoother, & & solver=solver,prefix=prefix_) end if - + end subroutine mld_s_base_onelev_dump diff --git a/mlprec/impl/level/mld_z_base_onelev_dump.f90 b/mlprec/impl/level/mld_z_base_onelev_dump.f90 index fab0c9ff..68227d61 100644 --- a/mlprec/impl/level/mld_z_base_onelev_dump.f90 +++ b/mlprec/impl/level/mld_z_base_onelev_dump.f90 @@ -36,8 +36,9 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -subroutine mld_z_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,solver) - +subroutine mld_z_base_onelev_dump(lv,level,info,prefix,head,ac,rp,& + & smoother,solver,global_num) + use psb_base_mod use mld_z_onelev_mod, mld_protect_name => mld_z_base_onelev_dump implicit none @@ -45,13 +46,14 @@ subroutine mld_z_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,solve integer(psb_ipk_), intent(in) :: level integer(psb_ipk_), intent(out) :: info character(len=*), intent(in), optional :: prefix, head - logical, optional, intent(in) :: ac, rp, smoother, solver + logical, optional, intent(in) :: ac, rp, smoother, solver,global_num + ! Local variables integer(psb_ipk_) :: i, j, il1, iln, lname, lev integer(psb_ipk_) :: icontxt,iam, np character(len=80) :: prefix_ character(len=120) :: fname ! len should be at least 20 more than - logical :: ac_, rp_ - ! len of prefix_ + logical :: ac_, rp_,global_num_ + integer(psb_ipk_), allocatable :: ivr(:), ivc(:) info = 0 @@ -78,23 +80,47 @@ subroutine mld_z_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,solve else rp_ = .false. end if + if (present(global_num)) then + global_num_ = global_num + else + global_num_ = .false. + end if lname = len_trim(prefix_) fname = trim(prefix_) write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam lname = lname + 5 - if (level >= 2) then - if (ac_) then - write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_ac.mtx' - call lv%ac%print(fname,head=head) + if (global_num_) then + if (level >= 2) then + if (ac_) then + ivr = lv%desc_ac%get_global_indices(owned=.false.) + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_ac.mtx' + call lv%ac%print(fname,head=head,iv=ivr) + end if + if (rp_) then + ivr = lv%map%p_desc_X%get_global_indices(owned=.false.) + ivc = lv%map%p_desc_Y%get_global_indices(owned=.false.) + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_r.mtx' + call lv%map%map_X2Y%print(fname,head=head,ivr=ivc,ivc=ivr) + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_p.mtx' + call lv%map%map_Y2X%print(fname,head=head,ivr=ivr,ivc=ivc) + end if end if - if (rp_) then - write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_r.mtx' - call lv%map%map_X2Y%print(fname,head=head) - write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_p.mtx' - call lv%map%map_Y2X%print(fname,head=head) + else + if (level >= 2) then + if (ac_) then + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_ac.mtx' + call lv%ac%print(fname,head=head) + end if + if (rp_) then + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_r.mtx' + call lv%map%map_X2Y%print(fname,head=head) + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_p.mtx' + call lv%map%map_Y2X%print(fname,head=head) + end if end if end if + if (allocated(lv%sm)) & & call lv%sm%dump(icontxt,level,info,smoother=smoother, & & solver=solver,prefix=prefix) @@ -103,5 +129,5 @@ subroutine mld_z_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,solve call lv%sm2a%dump(icontxt,level,info,smoother=smoother, & & solver=solver,prefix=prefix_) end if - + end subroutine mld_z_base_onelev_dump diff --git a/mlprec/impl/mld_c_extprol_bld.f90 b/mlprec/impl/mld_c_extprol_bld.f90 new file mode 100644 index 00000000..da912024 --- /dev/null +++ b/mlprec/impl/mld_c_extprol_bld.f90 @@ -0,0 +1,524 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +! File: mld_c_hierarchy_bld.f90 +! +! Subroutine: mld_c_hierarchy_bld +! Version: real +! +! This routine builds the preconditioner according to the requirements made by +! the user trough the subroutines mld_precinit and mld_precset. +! +! A multilevel preconditioner is regarded as an array of 'one-level' data structures, +! each containing the part of the preconditioner associated to a certain level, +! (for more details see the description of mld_Tonelev_type in mld_prec_type.f90). +! The levels are numbered in increasing order starting from the finest one, i.e. +! level 1 is the finest level. No transfer operators are associated to level 1. +! +! +! Arguments: +! a - type(psb_cspmat_type). +! The sparse matrix structure containing the local part of the +! matrix to be preconditioned. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of a. +! p - type(mld_cprec_type), input/output. +! The preconditioner data structure containing the local part +! of the preconditioner to be built. +! info - integer, output. +! Error code. +! +! amold - class(psb_c_base_sparse_mat), input, optional +! Mold for the inner format of matrices contained in the +! preconditioner +! +! +! vmold - class(psb_c_base_vect_type), input, optional +! Mold for the inner format of vectors contained in the +! preconditioner +! +! +! +subroutine mld_c_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold) + + use psb_base_mod + use mld_c_inner_mod + use mld_c_prec_mod, mld_protect_name => mld_c_extprol_bld + + Implicit None + + ! Arguments + type(psb_cspmat_type),intent(in), target :: a + type(psb_cspmat_type),intent(inout), target :: prolv(:) + type(psb_cspmat_type),intent(inout), target :: restrv(:) + type(psb_desc_type), intent(inout), target :: desc_a + type(mld_cprec_type),intent(inout),target :: p + integer(psb_ipk_), intent(out) :: info + class(psb_c_base_sparse_mat), intent(in), optional :: amold + class(psb_c_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + ! !$ character, intent(in), optional :: upd + + ! Local Variables + integer(psb_ipk_) :: ictxt, me,np + integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize, nplevs, mxplevs + integer(psb_ipk_) :: nprolv, nrestrv + real(psb_spk_) :: mnaggratio + integer(psb_ipk_) :: ipv(mld_ifpsz_), val + class(mld_c_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm + type(mld_sml_parms) :: baseparms, medparms, coarseparms + type(mld_c_onelev_type), allocatable :: tprecv(:) + integer(psb_ipk_) :: int_err(5) + character :: upd_ + integer(psb_ipk_) :: debug_level, debug_unit + character(len=20) :: name, ch_err + logical, parameter :: debug=.false. + + if (psb_get_errstatus().ne.0) return + info=psb_success_ + err=0 + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + name = 'mld_c_extprol_bld' + info = psb_success_ + int_err(1) = 0 + ictxt = desc_a%get_context() + call psb_info(ictxt, me, np) + p%ictxt = ictxt + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Entering ' + ! + ! For the time being we are commenting out the UPDATE argument + ! we plan to resurrect it later. + ! !$ if (present(upd)) then + ! !$ if (debug_level >= psb_debug_outer_) & + ! !$ & write(debug_unit,*) me,' ',trim(name),'UPD ', upd + ! !$ + ! !$ if ((psb_toupper(upd).eq.'F').or.(psb_toupper(upd).eq.'T')) then + ! !$ upd_=psb_toupper(upd) + ! !$ else + ! !$ upd_='F' + ! !$ endif + ! !$ else + ! !$ upd_='F' + ! !$ endif + upd_ = 'F' + + if (.not.allocated(p%precv)) then + !! Error: should have called mld_cprecinit + info=3111 + call psb_errpush(info,name) + goto 9999 + end if + + ! + ! Check to ensure all procs have the same + ! + newsz = -1 + casize = p%coarse_aggr_size + nplevs = p%n_prec_levs + mxplevs = p%max_prec_levs + mnaggratio = p%min_aggr_ratio + casize = p%coarse_aggr_size + iszv = size(p%precv) + nprolv = size(prolv) + nrestrv = size(restrv) + call psb_bcast(ictxt,iszv) + call psb_bcast(ictxt,casize) + call psb_bcast(ictxt,nplevs) + call psb_bcast(ictxt,mxplevs) + call psb_bcast(ictxt,mnaggratio) + call psb_bcast(ictxt,nprolv) + call psb_bcast(ictxt,nrestrv) + if (casize /= p%coarse_aggr_size) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent coarse_aggr_size') + goto 9999 + end if + if (nplevs /= p%n_prec_levs) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent n_prec_levs') + goto 9999 + end if + if (mxplevs /= p%max_prec_levs) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent max_prec_levs') + goto 9999 + end if + if (mnaggratio /= p%min_aggr_ratio) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent min_aggr_ratio') + goto 9999 + end if + if (iszv /= size(p%precv)) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent size of precv') + goto 9999 + end if + if (nprolv /= size(prolv)) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent size of prolv') + goto 9999 + end if + if (nrestrv /= size(restrv)) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent size of restrv') + goto 9999 + end if + if (nrestrv /= nprolv) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent size prolv vs restrv') + goto 9999 + end if + + if (iszv <= 1) then + ! We should only ever get here for multilevel. + info=psb_err_from_subroutine_ + ch_err='size bpv' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + + if (nrestrv < 1) then + ! We should only ever get here for multilevel. + info=psb_err_from_subroutine_ + ch_err='size restrv' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + + ! + nplevs = nrestrv + 1 + p%n_prec_levs = nplevs + + ! + ! Fixed number of levels. + ! + nplevs = max(itwo,min(nplevs,mxplevs)) + + coarseparms = p%precv(iszv)%parms + baseparms = p%precv(1)%parms + medparms = p%precv(2)%parms + + allocate(coarse_sm, source=p%precv(iszv)%sm,stat=info) + if (info == psb_success_) & + & allocate(med_sm, source=p%precv(2)%sm,stat=info) + if (info == psb_success_) & + & allocate(base_sm, source=p%precv(1)%sm,stat=info) + if (info /= psb_success_) then + write(0,*) 'Error in saving smoothers',info + call psb_errpush(psb_err_internal_error_,name,a_err='Base level precbuild.') + goto 9999 + end if + + if (iszv /= nplevs) then + allocate(tprecv(nplevs),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='prec reallocation') + goto 9999 + endif + tprecv(1)%parms = baseparms + allocate(tprecv(1)%sm,source=base_sm,stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='prec reallocation') + goto 9999 + endif + do i=2,nplevs-1 + tprecv(i)%parms = medparms + allocate(tprecv(i)%sm,source=med_sm,stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='prec reallocation') + goto 9999 + endif + end do + tprecv(nplevs)%parms = coarseparms + allocate(tprecv(nplevs)%sm,source=coarse_sm,stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='prec reallocation') + goto 9999 + endif + do i=1,iszv + call p%precv(i)%free(info) + end do + call move_alloc(tprecv,p%precv) + iszv = size(p%precv) + end if + ! + ! Finest level first; remember to fix base_a and base_desc + ! + p%precv(1)%base_a => a + p%precv(1)%base_desc => desc_a + newsz = 0 + array_build_loop: do i=2, iszv + + ! + ! Sanity checks on the parameters + ! + if (i p%precv(i)%ac + p%precv(i)%base_desc => p%precv(i)%desc_ac + + + if (i>2) then + if (all(p%precv(i)%map%naggr == p%precv(i-1)%map%naggr)) then + newsz=i-1 + end if + call psb_bcast(ictxt,newsz) + if (newsz > 0) exit array_build_loop + end if + end do array_build_loop + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Internal hierarchy build' ) + goto 9999 + endif + + iszv = size(p%precv) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Exiting with',iszv,' levels' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +contains + + subroutine mld_c_extaggr_bld(a,desc_a,p,op_restr,op_prol,info) + use psb_base_mod + use mld_c_inner_mod + + implicit none + + ! Arguments + type(psb_cspmat_type), intent(in), target :: a + type(psb_cspmat_type), intent(inout) :: op_restr,op_prol + type(psb_desc_type), intent(in), target :: desc_a + type(mld_c_onelev_type), intent(inout),target :: p + integer(psb_ipk_), intent(out) :: info + + ! Local variables + character(len=20) :: name + integer(psb_mpik_) :: ictxt, np, me, ncol + integer(psb_ipk_) :: err_act,ntaggr,nzl + integer(psb_ipk_), allocatable :: ilaggr(:), nlaggr(:) + type(psb_cspmat_type) :: ac, am3, am4 + type(psb_c_coo_sparse_mat) :: acoo, bcoo + type(psb_c_csr_sparse_mat) :: acsr1 + logical, parameter :: debug=.false. + + name='mld_c_extaggr_bld' + if (psb_get_errstatus().ne.0) return + call psb_erractionsave(err_act) + info = psb_success_ + ictxt = desc_a%get_context() + call psb_info(ictxt,me,np) + allocate(nlaggr(np),ilaggr(1)) + nlaggr = 0 + ilaggr = 0 + p%parms%aggr_alg = mld_ext_aggr_ + call mld_check_def(p%parms%ml_type,'Multilevel type',& + & mld_mult_ml_,is_legal_ml_type) + call mld_check_def(p%parms%coarse_mat,'Coarse matrix',& + & mld_distr_mat_,is_legal_ml_coarse_mat) + + nlaggr(me+1) = op_restr%get_nrows() + if (op_restr%get_nrows() /= op_prol%get_ncols()) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent restr/prol sizes') + goto 9999 + end if + call psb_sum(ictxt,nlaggr) + ntaggr = sum(nlaggr) + ncol = desc_a%get_local_cols() + if (debug) write(0,*)me,' Sizes:',op_restr%get_nrows(),op_restr%get_ncols(),& + & op_prol%get_nrows(),op_prol%get_ncols(), a%get_nrows(),a%get_ncols() + ! + ! Compute local part of AC + ! + call psb_spspmm(a,op_prol,am3,info) + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 2') + goto 9999 + end if + call psb_sphalo(am3,desc_a,am4,info,& + & colcnv=.false.,rowscale=.true.) + if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4) + if (info == psb_success_) call am4%free() + if(info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3') + goto 9999 + end if + call psb_spspmm(op_restr,am3,ac,info) + if (info == psb_success_) call am3%free() + if (info == psb_success_) call ac%cscnv(info,type='csr',dupl=psb_dupl_add_) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,a_err='Build ac = op_restr x am3') + goto 9999 + end if + + select case(p%parms%coarse_mat) + + case(mld_distr_mat_) + + call ac%mv_to(bcoo) + if (p%parms%clean_zeros) call bcoo%clean_zeros(info) + nzl = bcoo%get_nzeros() + + if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1)) + if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info) + if (info == psb_success_) call psb_cdasb(p%desc_ac,info) + if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I') + if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I') + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Creating p%desc_ac and converting ac') + goto 9999 + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Assembld aux descr. distr.' + call p%ac%mv_from(bcoo) + + call p%ac%set_nrows(p%desc_ac%get_local_rows()) + call p%ac%set_ncols(p%desc_ac%get_local_cols()) + call p%ac%set_asb() + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free') + goto 9999 + end if + + if (np>1) then + call op_prol%mv_to(acsr1) + nzl = acsr1%get_nzeros() + call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I') + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc') + goto 9999 + end if + call op_prol%mv_from(acsr1) + endif + call op_prol%set_ncols(p%desc_ac%get_local_cols()) + + if (np>1) then + call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_) + call op_restr%mv_to(acoo) + nzl = acoo%get_nzeros() + if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I') + call acoo%set_dupl(psb_dupl_add_) + if (info == psb_success_) call op_restr%mv_from(acoo) + if (info == psb_success_) call op_restr%cscnv(info,type='csr') + if(info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Converting op_restr to local') + goto 9999 + end if + end if + call op_restr%set_nrows(p%desc_ac%get_local_cols()) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Done ac ' + + case(mld_repl_mat_) + ! + ! + call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.) + if (info == psb_success_) call psb_cdasb(p%desc_ac,info) + if ((info == psb_success_).and.p%parms%clean_zeros) call ac%clean_zeros(info) + if (info == psb_success_) & + & call psb_gather(p%ac,ac,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.) + + if (info /= psb_success_) goto 9999 + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='invalid mld_coarse_mat_') + goto 9999 + end select + + call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv') + goto 9999 + end if + + ! + ! Copy the prolongation/restriction matrices into the descriptor map. + ! op_restr => PR^T i.e. restriction operator + ! op_prol => PR i.e. prolongation operator + ! + + p%map = psb_linmap(psb_map_aggr_,desc_a,& + & p%desc_ac,op_restr,op_prol,ilaggr,nlaggr) + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + end subroutine mld_c_extaggr_bld + +end subroutine mld_c_extprol_bld diff --git a/mlprec/impl/mld_s_extprol_bld.f90 b/mlprec/impl/mld_s_extprol_bld.f90 new file mode 100644 index 00000000..ed429863 --- /dev/null +++ b/mlprec/impl/mld_s_extprol_bld.f90 @@ -0,0 +1,524 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +! File: mld_s_hierarchy_bld.f90 +! +! Subroutine: mld_s_hierarchy_bld +! Version: real +! +! This routine builds the preconditioner according to the requirements made by +! the user trough the subroutines mld_precinit and mld_precset. +! +! A multilevel preconditioner is regarded as an array of 'one-level' data structures, +! each containing the part of the preconditioner associated to a certain level, +! (for more details see the description of mld_Tonelev_type in mld_prec_type.f90). +! The levels are numbered in increasing order starting from the finest one, i.e. +! level 1 is the finest level. No transfer operators are associated to level 1. +! +! +! Arguments: +! a - type(psb_sspmat_type). +! The sparse matrix structure containing the local part of the +! matrix to be preconditioned. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of a. +! p - type(mld_sprec_type), input/output. +! The preconditioner data structure containing the local part +! of the preconditioner to be built. +! info - integer, output. +! Error code. +! +! amold - class(psb_s_base_sparse_mat), input, optional +! Mold for the inner format of matrices contained in the +! preconditioner +! +! +! vmold - class(psb_s_base_vect_type), input, optional +! Mold for the inner format of vectors contained in the +! preconditioner +! +! +! +subroutine mld_s_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold) + + use psb_base_mod + use mld_s_inner_mod + use mld_s_prec_mod, mld_protect_name => mld_s_extprol_bld + + Implicit None + + ! Arguments + type(psb_sspmat_type),intent(in), target :: a + type(psb_sspmat_type),intent(inout), target :: prolv(:) + type(psb_sspmat_type),intent(inout), target :: restrv(:) + type(psb_desc_type), intent(inout), target :: desc_a + type(mld_sprec_type),intent(inout),target :: p + integer(psb_ipk_), intent(out) :: info + class(psb_s_base_sparse_mat), intent(in), optional :: amold + class(psb_s_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + ! !$ character, intent(in), optional :: upd + + ! Local Variables + integer(psb_ipk_) :: ictxt, me,np + integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize, nplevs, mxplevs + integer(psb_ipk_) :: nprolv, nrestrv + real(psb_spk_) :: mnaggratio + integer(psb_ipk_) :: ipv(mld_ifpsz_), val + class(mld_s_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm + type(mld_sml_parms) :: baseparms, medparms, coarseparms + type(mld_s_onelev_type), allocatable :: tprecv(:) + integer(psb_ipk_) :: int_err(5) + character :: upd_ + integer(psb_ipk_) :: debug_level, debug_unit + character(len=20) :: name, ch_err + logical, parameter :: debug=.false. + + if (psb_get_errstatus().ne.0) return + info=psb_success_ + err=0 + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + name = 'mld_s_extprol_bld' + info = psb_success_ + int_err(1) = 0 + ictxt = desc_a%get_context() + call psb_info(ictxt, me, np) + p%ictxt = ictxt + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Entering ' + ! + ! For the time being we are commenting out the UPDATE argument + ! we plan to resurrect it later. + ! !$ if (present(upd)) then + ! !$ if (debug_level >= psb_debug_outer_) & + ! !$ & write(debug_unit,*) me,' ',trim(name),'UPD ', upd + ! !$ + ! !$ if ((psb_toupper(upd).eq.'F').or.(psb_toupper(upd).eq.'T')) then + ! !$ upd_=psb_toupper(upd) + ! !$ else + ! !$ upd_='F' + ! !$ endif + ! !$ else + ! !$ upd_='F' + ! !$ endif + upd_ = 'F' + + if (.not.allocated(p%precv)) then + !! Error: should have called mld_sprecinit + info=3111 + call psb_errpush(info,name) + goto 9999 + end if + + ! + ! Check to ensure all procs have the same + ! + newsz = -1 + casize = p%coarse_aggr_size + nplevs = p%n_prec_levs + mxplevs = p%max_prec_levs + mnaggratio = p%min_aggr_ratio + casize = p%coarse_aggr_size + iszv = size(p%precv) + nprolv = size(prolv) + nrestrv = size(restrv) + call psb_bcast(ictxt,iszv) + call psb_bcast(ictxt,casize) + call psb_bcast(ictxt,nplevs) + call psb_bcast(ictxt,mxplevs) + call psb_bcast(ictxt,mnaggratio) + call psb_bcast(ictxt,nprolv) + call psb_bcast(ictxt,nrestrv) + if (casize /= p%coarse_aggr_size) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent coarse_aggr_size') + goto 9999 + end if + if (nplevs /= p%n_prec_levs) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent n_prec_levs') + goto 9999 + end if + if (mxplevs /= p%max_prec_levs) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent max_prec_levs') + goto 9999 + end if + if (mnaggratio /= p%min_aggr_ratio) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent min_aggr_ratio') + goto 9999 + end if + if (iszv /= size(p%precv)) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent size of precv') + goto 9999 + end if + if (nprolv /= size(prolv)) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent size of prolv') + goto 9999 + end if + if (nrestrv /= size(restrv)) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent size of restrv') + goto 9999 + end if + if (nrestrv /= nprolv) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent size prolv vs restrv') + goto 9999 + end if + + if (iszv <= 1) then + ! We should only ever get here for multilevel. + info=psb_err_from_subroutine_ + ch_err='size bpv' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + + if (nrestrv < 1) then + ! We should only ever get here for multilevel. + info=psb_err_from_subroutine_ + ch_err='size restrv' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + + ! + nplevs = nrestrv + 1 + p%n_prec_levs = nplevs + + ! + ! Fixed number of levels. + ! + nplevs = max(itwo,min(nplevs,mxplevs)) + + coarseparms = p%precv(iszv)%parms + baseparms = p%precv(1)%parms + medparms = p%precv(2)%parms + + allocate(coarse_sm, source=p%precv(iszv)%sm,stat=info) + if (info == psb_success_) & + & allocate(med_sm, source=p%precv(2)%sm,stat=info) + if (info == psb_success_) & + & allocate(base_sm, source=p%precv(1)%sm,stat=info) + if (info /= psb_success_) then + write(0,*) 'Error in saving smoothers',info + call psb_errpush(psb_err_internal_error_,name,a_err='Base level precbuild.') + goto 9999 + end if + + if (iszv /= nplevs) then + allocate(tprecv(nplevs),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='prec reallocation') + goto 9999 + endif + tprecv(1)%parms = baseparms + allocate(tprecv(1)%sm,source=base_sm,stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='prec reallocation') + goto 9999 + endif + do i=2,nplevs-1 + tprecv(i)%parms = medparms + allocate(tprecv(i)%sm,source=med_sm,stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='prec reallocation') + goto 9999 + endif + end do + tprecv(nplevs)%parms = coarseparms + allocate(tprecv(nplevs)%sm,source=coarse_sm,stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='prec reallocation') + goto 9999 + endif + do i=1,iszv + call p%precv(i)%free(info) + end do + call move_alloc(tprecv,p%precv) + iszv = size(p%precv) + end if + ! + ! Finest level first; remember to fix base_a and base_desc + ! + p%precv(1)%base_a => a + p%precv(1)%base_desc => desc_a + newsz = 0 + array_build_loop: do i=2, iszv + + ! + ! Sanity checks on the parameters + ! + if (i p%precv(i)%ac + p%precv(i)%base_desc => p%precv(i)%desc_ac + + + if (i>2) then + if (all(p%precv(i)%map%naggr == p%precv(i-1)%map%naggr)) then + newsz=i-1 + end if + call psb_bcast(ictxt,newsz) + if (newsz > 0) exit array_build_loop + end if + end do array_build_loop + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Internal hierarchy build' ) + goto 9999 + endif + + iszv = size(p%precv) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Exiting with',iszv,' levels' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +contains + + subroutine mld_s_extaggr_bld(a,desc_a,p,op_restr,op_prol,info) + use psb_base_mod + use mld_s_inner_mod + + implicit none + + ! Arguments + type(psb_sspmat_type), intent(in), target :: a + type(psb_sspmat_type), intent(inout) :: op_restr,op_prol + type(psb_desc_type), intent(in), target :: desc_a + type(mld_s_onelev_type), intent(inout),target :: p + integer(psb_ipk_), intent(out) :: info + + ! Local variables + character(len=20) :: name + integer(psb_mpik_) :: ictxt, np, me, ncol + integer(psb_ipk_) :: err_act,ntaggr,nzl + integer(psb_ipk_), allocatable :: ilaggr(:), nlaggr(:) + type(psb_sspmat_type) :: ac, am3, am4 + type(psb_s_coo_sparse_mat) :: acoo, bcoo + type(psb_s_csr_sparse_mat) :: acsr1 + logical, parameter :: debug=.false. + + name='mld_s_extaggr_bld' + if (psb_get_errstatus().ne.0) return + call psb_erractionsave(err_act) + info = psb_success_ + ictxt = desc_a%get_context() + call psb_info(ictxt,me,np) + allocate(nlaggr(np),ilaggr(1)) + nlaggr = 0 + ilaggr = 0 + p%parms%aggr_alg = mld_ext_aggr_ + call mld_check_def(p%parms%ml_type,'Multilevel type',& + & mld_mult_ml_,is_legal_ml_type) + call mld_check_def(p%parms%coarse_mat,'Coarse matrix',& + & mld_distr_mat_,is_legal_ml_coarse_mat) + + nlaggr(me+1) = op_restr%get_nrows() + if (op_restr%get_nrows() /= op_prol%get_ncols()) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent restr/prol sizes') + goto 9999 + end if + call psb_sum(ictxt,nlaggr) + ntaggr = sum(nlaggr) + ncol = desc_a%get_local_cols() + if (debug) write(0,*)me,' Sizes:',op_restr%get_nrows(),op_restr%get_ncols(),& + & op_prol%get_nrows(),op_prol%get_ncols(), a%get_nrows(),a%get_ncols() + ! + ! Compute local part of AC + ! + call psb_spspmm(a,op_prol,am3,info) + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 2') + goto 9999 + end if + call psb_sphalo(am3,desc_a,am4,info,& + & colcnv=.false.,rowscale=.true.) + if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4) + if (info == psb_success_) call am4%free() + if(info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3') + goto 9999 + end if + call psb_spspmm(op_restr,am3,ac,info) + if (info == psb_success_) call am3%free() + if (info == psb_success_) call ac%cscnv(info,type='csr',dupl=psb_dupl_add_) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,a_err='Build ac = op_restr x am3') + goto 9999 + end if + + select case(p%parms%coarse_mat) + + case(mld_distr_mat_) + + call ac%mv_to(bcoo) + if (p%parms%clean_zeros) call bcoo%clean_zeros(info) + nzl = bcoo%get_nzeros() + + if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1)) + if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info) + if (info == psb_success_) call psb_cdasb(p%desc_ac,info) + if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I') + if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I') + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Creating p%desc_ac and converting ac') + goto 9999 + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Assembld aux descr. distr.' + call p%ac%mv_from(bcoo) + + call p%ac%set_nrows(p%desc_ac%get_local_rows()) + call p%ac%set_ncols(p%desc_ac%get_local_cols()) + call p%ac%set_asb() + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free') + goto 9999 + end if + + if (np>1) then + call op_prol%mv_to(acsr1) + nzl = acsr1%get_nzeros() + call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I') + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc') + goto 9999 + end if + call op_prol%mv_from(acsr1) + endif + call op_prol%set_ncols(p%desc_ac%get_local_cols()) + + if (np>1) then + call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_) + call op_restr%mv_to(acoo) + nzl = acoo%get_nzeros() + if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I') + call acoo%set_dupl(psb_dupl_add_) + if (info == psb_success_) call op_restr%mv_from(acoo) + if (info == psb_success_) call op_restr%cscnv(info,type='csr') + if(info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Converting op_restr to local') + goto 9999 + end if + end if + call op_restr%set_nrows(p%desc_ac%get_local_cols()) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Done ac ' + + case(mld_repl_mat_) + ! + ! + call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.) + if (info == psb_success_) call psb_cdasb(p%desc_ac,info) + if ((info == psb_success_).and.p%parms%clean_zeros) call ac%clean_zeros(info) + if (info == psb_success_) & + & call psb_gather(p%ac,ac,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.) + + if (info /= psb_success_) goto 9999 + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='invalid mld_coarse_mat_') + goto 9999 + end select + + call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv') + goto 9999 + end if + + ! + ! Copy the prolongation/restriction matrices into the descriptor map. + ! op_restr => PR^T i.e. restriction operator + ! op_prol => PR i.e. prolongation operator + ! + + p%map = psb_linmap(psb_map_aggr_,desc_a,& + & p%desc_ac,op_restr,op_prol,ilaggr,nlaggr) + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + end subroutine mld_s_extaggr_bld + +end subroutine mld_s_extprol_bld diff --git a/mlprec/impl/mld_z_extprol_bld.f90 b/mlprec/impl/mld_z_extprol_bld.f90 new file mode 100644 index 00000000..f6a63d85 --- /dev/null +++ b/mlprec/impl/mld_z_extprol_bld.f90 @@ -0,0 +1,524 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +! File: mld_z_hierarchy_bld.f90 +! +! Subroutine: mld_z_hierarchy_bld +! Version: real +! +! This routine builds the preconditioner according to the requirements made by +! the user trough the subroutines mld_precinit and mld_precset. +! +! A multilevel preconditioner is regarded as an array of 'one-level' data structures, +! each containing the part of the preconditioner associated to a certain level, +! (for more details see the description of mld_Tonelev_type in mld_prec_type.f90). +! The levels are numbered in increasing order starting from the finest one, i.e. +! level 1 is the finest level. No transfer operators are associated to level 1. +! +! +! Arguments: +! a - type(psb_zspmat_type). +! The sparse matrix structure containing the local part of the +! matrix to be preconditioned. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of a. +! p - type(mld_zprec_type), input/output. +! The preconditioner data structure containing the local part +! of the preconditioner to be built. +! info - integer, output. +! Error code. +! +! amold - class(psb_z_base_sparse_mat), input, optional +! Mold for the inner format of matrices contained in the +! preconditioner +! +! +! vmold - class(psb_z_base_vect_type), input, optional +! Mold for the inner format of vectors contained in the +! preconditioner +! +! +! +subroutine mld_z_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold) + + use psb_base_mod + use mld_z_inner_mod + use mld_z_prec_mod, mld_protect_name => mld_z_extprol_bld + + Implicit None + + ! Arguments + type(psb_zspmat_type),intent(in), target :: a + type(psb_zspmat_type),intent(inout), target :: prolv(:) + type(psb_zspmat_type),intent(inout), target :: restrv(:) + type(psb_desc_type), intent(inout), target :: desc_a + type(mld_zprec_type),intent(inout),target :: p + integer(psb_ipk_), intent(out) :: info + class(psb_z_base_sparse_mat), intent(in), optional :: amold + class(psb_z_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + ! !$ character, intent(in), optional :: upd + + ! Local Variables + integer(psb_ipk_) :: ictxt, me,np + integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize, nplevs, mxplevs + integer(psb_ipk_) :: nprolv, nrestrv + real(psb_dpk_) :: mnaggratio + integer(psb_ipk_) :: ipv(mld_ifpsz_), val + class(mld_z_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm + type(mld_dml_parms) :: baseparms, medparms, coarseparms + type(mld_z_onelev_type), allocatable :: tprecv(:) + integer(psb_ipk_) :: int_err(5) + character :: upd_ + integer(psb_ipk_) :: debug_level, debug_unit + character(len=20) :: name, ch_err + logical, parameter :: debug=.false. + + if (psb_get_errstatus().ne.0) return + info=psb_success_ + err=0 + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + name = 'mld_z_extprol_bld' + info = psb_success_ + int_err(1) = 0 + ictxt = desc_a%get_context() + call psb_info(ictxt, me, np) + p%ictxt = ictxt + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Entering ' + ! + ! For the time being we are commenting out the UPDATE argument + ! we plan to resurrect it later. + ! !$ if (present(upd)) then + ! !$ if (debug_level >= psb_debug_outer_) & + ! !$ & write(debug_unit,*) me,' ',trim(name),'UPD ', upd + ! !$ + ! !$ if ((psb_toupper(upd).eq.'F').or.(psb_toupper(upd).eq.'T')) then + ! !$ upd_=psb_toupper(upd) + ! !$ else + ! !$ upd_='F' + ! !$ endif + ! !$ else + ! !$ upd_='F' + ! !$ endif + upd_ = 'F' + + if (.not.allocated(p%precv)) then + !! Error: should have called mld_zprecinit + info=3111 + call psb_errpush(info,name) + goto 9999 + end if + + ! + ! Check to ensure all procs have the same + ! + newsz = -1 + casize = p%coarse_aggr_size + nplevs = p%n_prec_levs + mxplevs = p%max_prec_levs + mnaggratio = p%min_aggr_ratio + casize = p%coarse_aggr_size + iszv = size(p%precv) + nprolv = size(prolv) + nrestrv = size(restrv) + call psb_bcast(ictxt,iszv) + call psb_bcast(ictxt,casize) + call psb_bcast(ictxt,nplevs) + call psb_bcast(ictxt,mxplevs) + call psb_bcast(ictxt,mnaggratio) + call psb_bcast(ictxt,nprolv) + call psb_bcast(ictxt,nrestrv) + if (casize /= p%coarse_aggr_size) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent coarse_aggr_size') + goto 9999 + end if + if (nplevs /= p%n_prec_levs) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent n_prec_levs') + goto 9999 + end if + if (mxplevs /= p%max_prec_levs) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent max_prec_levs') + goto 9999 + end if + if (mnaggratio /= p%min_aggr_ratio) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent min_aggr_ratio') + goto 9999 + end if + if (iszv /= size(p%precv)) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent size of precv') + goto 9999 + end if + if (nprolv /= size(prolv)) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent size of prolv') + goto 9999 + end if + if (nrestrv /= size(restrv)) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent size of restrv') + goto 9999 + end if + if (nrestrv /= nprolv) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent size prolv vs restrv') + goto 9999 + end if + + if (iszv <= 1) then + ! We should only ever get here for multilevel. + info=psb_err_from_subroutine_ + ch_err='size bpv' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + + if (nrestrv < 1) then + ! We should only ever get here for multilevel. + info=psb_err_from_subroutine_ + ch_err='size restrv' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + + ! + nplevs = nrestrv + 1 + p%n_prec_levs = nplevs + + ! + ! Fixed number of levels. + ! + nplevs = max(itwo,min(nplevs,mxplevs)) + + coarseparms = p%precv(iszv)%parms + baseparms = p%precv(1)%parms + medparms = p%precv(2)%parms + + allocate(coarse_sm, source=p%precv(iszv)%sm,stat=info) + if (info == psb_success_) & + & allocate(med_sm, source=p%precv(2)%sm,stat=info) + if (info == psb_success_) & + & allocate(base_sm, source=p%precv(1)%sm,stat=info) + if (info /= psb_success_) then + write(0,*) 'Error in saving smoothers',info + call psb_errpush(psb_err_internal_error_,name,a_err='Base level precbuild.') + goto 9999 + end if + + if (iszv /= nplevs) then + allocate(tprecv(nplevs),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='prec reallocation') + goto 9999 + endif + tprecv(1)%parms = baseparms + allocate(tprecv(1)%sm,source=base_sm,stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='prec reallocation') + goto 9999 + endif + do i=2,nplevs-1 + tprecv(i)%parms = medparms + allocate(tprecv(i)%sm,source=med_sm,stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='prec reallocation') + goto 9999 + endif + end do + tprecv(nplevs)%parms = coarseparms + allocate(tprecv(nplevs)%sm,source=coarse_sm,stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='prec reallocation') + goto 9999 + endif + do i=1,iszv + call p%precv(i)%free(info) + end do + call move_alloc(tprecv,p%precv) + iszv = size(p%precv) + end if + ! + ! Finest level first; remember to fix base_a and base_desc + ! + p%precv(1)%base_a => a + p%precv(1)%base_desc => desc_a + newsz = 0 + array_build_loop: do i=2, iszv + + ! + ! Sanity checks on the parameters + ! + if (i p%precv(i)%ac + p%precv(i)%base_desc => p%precv(i)%desc_ac + + + if (i>2) then + if (all(p%precv(i)%map%naggr == p%precv(i-1)%map%naggr)) then + newsz=i-1 + end if + call psb_bcast(ictxt,newsz) + if (newsz > 0) exit array_build_loop + end if + end do array_build_loop + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Internal hierarchy build' ) + goto 9999 + endif + + iszv = size(p%precv) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Exiting with',iszv,' levels' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +contains + + subroutine mld_z_extaggr_bld(a,desc_a,p,op_restr,op_prol,info) + use psb_base_mod + use mld_z_inner_mod + + implicit none + + ! Arguments + type(psb_zspmat_type), intent(in), target :: a + type(psb_zspmat_type), intent(inout) :: op_restr,op_prol + type(psb_desc_type), intent(in), target :: desc_a + type(mld_z_onelev_type), intent(inout),target :: p + integer(psb_ipk_), intent(out) :: info + + ! Local variables + character(len=20) :: name + integer(psb_mpik_) :: ictxt, np, me, ncol + integer(psb_ipk_) :: err_act,ntaggr,nzl + integer(psb_ipk_), allocatable :: ilaggr(:), nlaggr(:) + type(psb_zspmat_type) :: ac, am3, am4 + type(psb_z_coo_sparse_mat) :: acoo, bcoo + type(psb_z_csr_sparse_mat) :: acsr1 + logical, parameter :: debug=.false. + + name='mld_z_extaggr_bld' + if (psb_get_errstatus().ne.0) return + call psb_erractionsave(err_act) + info = psb_success_ + ictxt = desc_a%get_context() + call psb_info(ictxt,me,np) + allocate(nlaggr(np),ilaggr(1)) + nlaggr = 0 + ilaggr = 0 + p%parms%aggr_alg = mld_ext_aggr_ + call mld_check_def(p%parms%ml_type,'Multilevel type',& + & mld_mult_ml_,is_legal_ml_type) + call mld_check_def(p%parms%coarse_mat,'Coarse matrix',& + & mld_distr_mat_,is_legal_ml_coarse_mat) + + nlaggr(me+1) = op_restr%get_nrows() + if (op_restr%get_nrows() /= op_prol%get_ncols()) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Inconsistent restr/prol sizes') + goto 9999 + end if + call psb_sum(ictxt,nlaggr) + ntaggr = sum(nlaggr) + ncol = desc_a%get_local_cols() + if (debug) write(0,*)me,' Sizes:',op_restr%get_nrows(),op_restr%get_ncols(),& + & op_prol%get_nrows(),op_prol%get_ncols(), a%get_nrows(),a%get_ncols() + ! + ! Compute local part of AC + ! + call psb_spspmm(a,op_prol,am3,info) + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 2') + goto 9999 + end if + call psb_sphalo(am3,desc_a,am4,info,& + & colcnv=.false.,rowscale=.true.) + if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4) + if (info == psb_success_) call am4%free() + if(info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3') + goto 9999 + end if + call psb_spspmm(op_restr,am3,ac,info) + if (info == psb_success_) call am3%free() + if (info == psb_success_) call ac%cscnv(info,type='csr',dupl=psb_dupl_add_) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,a_err='Build ac = op_restr x am3') + goto 9999 + end if + + select case(p%parms%coarse_mat) + + case(mld_distr_mat_) + + call ac%mv_to(bcoo) + if (p%parms%clean_zeros) call bcoo%clean_zeros(info) + nzl = bcoo%get_nzeros() + + if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1)) + if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info) + if (info == psb_success_) call psb_cdasb(p%desc_ac,info) + if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I') + if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I') + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Creating p%desc_ac and converting ac') + goto 9999 + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Assembld aux descr. distr.' + call p%ac%mv_from(bcoo) + + call p%ac%set_nrows(p%desc_ac%get_local_rows()) + call p%ac%set_ncols(p%desc_ac%get_local_cols()) + call p%ac%set_asb() + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free') + goto 9999 + end if + + if (np>1) then + call op_prol%mv_to(acsr1) + nzl = acsr1%get_nzeros() + call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I') + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc') + goto 9999 + end if + call op_prol%mv_from(acsr1) + endif + call op_prol%set_ncols(p%desc_ac%get_local_cols()) + + if (np>1) then + call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_) + call op_restr%mv_to(acoo) + nzl = acoo%get_nzeros() + if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I') + call acoo%set_dupl(psb_dupl_add_) + if (info == psb_success_) call op_restr%mv_from(acoo) + if (info == psb_success_) call op_restr%cscnv(info,type='csr') + if(info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Converting op_restr to local') + goto 9999 + end if + end if + call op_restr%set_nrows(p%desc_ac%get_local_cols()) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Done ac ' + + case(mld_repl_mat_) + ! + ! + call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.) + if (info == psb_success_) call psb_cdasb(p%desc_ac,info) + if ((info == psb_success_).and.p%parms%clean_zeros) call ac%clean_zeros(info) + if (info == psb_success_) & + & call psb_gather(p%ac,ac,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.) + + if (info /= psb_success_) goto 9999 + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='invalid mld_coarse_mat_') + goto 9999 + end select + + call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv') + goto 9999 + end if + + ! + ! Copy the prolongation/restriction matrices into the descriptor map. + ! op_restr => PR^T i.e. restriction operator + ! op_prol => PR i.e. prolongation operator + ! + + p%map = psb_linmap(psb_map_aggr_,desc_a,& + & p%desc_ac,op_restr,op_prol,ilaggr,nlaggr) + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + end subroutine mld_z_extaggr_bld + +end subroutine mld_z_extprol_bld diff --git a/mlprec/mld_c_onelev_mod.f90 b/mlprec/mld_c_onelev_mod.f90 index 2c39c3d6..26959820 100644 --- a/mlprec/mld_c_onelev_mod.f90 +++ b/mlprec/mld_c_onelev_mod.f90 @@ -354,7 +354,8 @@ module mld_c_onelev_mod end interface interface - subroutine mld_c_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,solver) + subroutine mld_c_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,& + & solver,global_num) import :: psb_cspmat_type, psb_c_vect_type, psb_c_base_vect_type, & & psb_clinmap_type, psb_spk_, mld_c_onelev_type, & & psb_ipk_, psb_long_int_k_, psb_desc_type @@ -363,7 +364,7 @@ module mld_c_onelev_mod integer(psb_ipk_), intent(in) :: level integer(psb_ipk_), intent(out) :: info character(len=*), intent(in), optional :: prefix, head - logical, optional, intent(in) :: ac, rp, smoother, solver + logical, optional, intent(in) :: ac, rp, smoother, solver, global_num end subroutine mld_c_base_onelev_dump end interface diff --git a/mlprec/mld_c_prec_mod.f90 b/mlprec/mld_c_prec_mod.f90 index 259e8167..5d0fea10 100644 --- a/mlprec/mld_c_prec_mod.f90 +++ b/mlprec/mld_c_prec_mod.f90 @@ -52,7 +52,7 @@ module mld_c_prec_mod use mld_c_diag_solver use mld_c_ilu_solver use mld_c_gs_solver - + interface mld_precinit subroutine mld_cprecinit(p,ptype,info,nlev) import :: psb_cspmat_type, psb_desc_type, psb_spk_, & @@ -62,13 +62,13 @@ module mld_c_prec_mod integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: nlev end subroutine mld_cprecinit - end interface + end interface mld_precinit interface mld_precset module procedure mld_c_iprecsetsm, mld_c_iprecsetsv, & & mld_c_iprecseti, mld_c_iprecsetc, mld_c_iprecsetr, & & mld_c_cprecseti, mld_c_cprecsetc, mld_c_cprecsetr - end interface + end interface mld_precset interface mld_precbld subroutine mld_cprecbld(a,desc_a,prec,info,amold,vmold,imold) @@ -83,9 +83,9 @@ module mld_c_prec_mod class(psb_c_base_sparse_mat), intent(in), optional :: amold class(psb_c_base_vect_type), intent(in), optional :: vmold class(psb_i_base_vect_type), intent(in), optional :: imold - ! character, intent(in),optional :: upd + ! character, intent(in),optional :: upd end subroutine mld_cprecbld - end interface + end interface mld_precbld interface mld_hierarchy_bld subroutine mld_c_hierarchy_bld(a,desc_a,prec,info,amold,vmold,imold) @@ -100,10 +100,30 @@ module mld_c_prec_mod class(psb_c_base_sparse_mat), intent(in), optional :: amold class(psb_c_base_vect_type), intent(in), optional :: vmold class(psb_i_base_vect_type), intent(in), optional :: imold - ! character, intent(in),optional :: upd + ! character, intent(in),optional :: upd end subroutine mld_c_hierarchy_bld - end interface + end interface mld_hierarchy_bld + + interface mld_extprol_bld + subroutine mld_c_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold) + import :: psb_cspmat_type, psb_desc_type, psb_spk_, & + & psb_c_base_sparse_mat, psb_c_base_vect_type, & + & psb_i_base_vect_type, mld_cprec_type, psb_ipk_ + ! Arguments + type(psb_cspmat_type),intent(in), target :: a + type(psb_cspmat_type),intent(inout), target :: prolv(:) + type(psb_cspmat_type),intent(inout), target :: restrv(:) + type(psb_desc_type), intent(inout), target :: desc_a + type(mld_cprec_type),intent(inout),target :: p + integer(psb_ipk_), intent(out) :: info + class(psb_c_base_sparse_mat), intent(in), optional :: amold + class(psb_c_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + ! !$ character, intent(in), optional :: upd + end subroutine mld_c_extprol_bld + end interface mld_extprol_bld + interface mld_ml_prec_bld subroutine mld_c_ml_prec_bld(a,desc_a,prec,info,amold,vmold,imold) import :: psb_cspmat_type, psb_desc_type, psb_spk_, & @@ -117,9 +137,9 @@ module mld_c_prec_mod class(psb_c_base_sparse_mat), intent(in), optional :: amold class(psb_c_base_vect_type), intent(in), optional :: vmold class(psb_i_base_vect_type), intent(in), optional :: imold - ! character, intent(in),optional :: upd + ! character, intent(in),optional :: upd end subroutine mld_c_ml_prec_bld - end interface + end interface mld_ml_prec_bld contains diff --git a/mlprec/mld_c_prec_type.f90 b/mlprec/mld_c_prec_type.f90 index ca296236..7bb565c8 100644 --- a/mlprec/mld_c_prec_type.f90 +++ b/mlprec/mld_c_prec_type.f90 @@ -744,14 +744,14 @@ contains end subroutine mld_c_apply1v - subroutine mld_c_dump(prec,info,istart,iend,prefix,head,ac,rp,smoother,solver) + subroutine mld_c_dump(prec,info,istart,iend,prefix,head,ac,rp,smoother,solver,global_num) implicit none class(mld_cprec_type), intent(in) :: prec integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: istart, iend character(len=*), intent(in), optional :: prefix, head - logical, optional, intent(in) :: smoother, solver,ac, rp + logical, optional, intent(in) :: smoother, solver,ac, rp, global_num integer(psb_ipk_) :: i, j, il1, iln, lname, lev integer(psb_ipk_) :: icontxt,iam, np character(len=80) :: prefix_ @@ -772,7 +772,7 @@ contains do lev=il1, iln call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,& - & ac=ac,smoother=smoother,solver=solver,rp=rp) + & ac=ac,smoother=smoother,solver=solver,rp=rp,global_num=global_num) end do end subroutine mld_c_dump diff --git a/mlprec/mld_d_onelev_mod.f90 b/mlprec/mld_d_onelev_mod.f90 index 8dc2e240..bbc6f8a4 100644 --- a/mlprec/mld_d_onelev_mod.f90 +++ b/mlprec/mld_d_onelev_mod.f90 @@ -354,7 +354,8 @@ module mld_d_onelev_mod end interface interface - subroutine mld_d_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,solver) + subroutine mld_d_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,& + & solver,global_num) import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, & & psb_dlinmap_type, psb_dpk_, mld_d_onelev_type, & & psb_ipk_, psb_long_int_k_, psb_desc_type @@ -363,7 +364,7 @@ module mld_d_onelev_mod integer(psb_ipk_), intent(in) :: level integer(psb_ipk_), intent(out) :: info character(len=*), intent(in), optional :: prefix, head - logical, optional, intent(in) :: ac, rp, smoother, solver + logical, optional, intent(in) :: ac, rp, smoother, solver, global_num end subroutine mld_d_base_onelev_dump end interface diff --git a/mlprec/mld_d_prec_mod.f90 b/mlprec/mld_d_prec_mod.f90 index 48919fc3..25a6ba3c 100644 --- a/mlprec/mld_d_prec_mod.f90 +++ b/mlprec/mld_d_prec_mod.f90 @@ -109,20 +109,21 @@ module mld_d_prec_mod import :: psb_dspmat_type, psb_desc_type, psb_dpk_, & & psb_d_base_sparse_mat, psb_d_base_vect_type, & & psb_i_base_vect_type, mld_dprec_type, psb_ipk_ - implicit none + ! Arguments type(psb_dspmat_type),intent(in), target :: a type(psb_dspmat_type),intent(inout), target :: prolv(:) type(psb_dspmat_type),intent(inout), target :: restrv(:) - type(psb_desc_type), intent(inout), target :: desc_a + type(psb_desc_type), intent(inout), target :: desc_a type(mld_dprec_type),intent(inout),target :: p integer(psb_ipk_), intent(out) :: info class(psb_d_base_sparse_mat), intent(in), optional :: amold class(psb_d_base_vect_type), intent(in), optional :: vmold class(psb_i_base_vect_type), intent(in), optional :: imold + ! !$ character, intent(in), optional :: upd end subroutine mld_d_extprol_bld end interface mld_extprol_bld - + interface mld_ml_prec_bld subroutine mld_d_ml_prec_bld(a,desc_a,prec,info,amold,vmold,imold) import :: psb_dspmat_type, psb_desc_type, psb_dpk_, & diff --git a/mlprec/mld_d_prec_type.f90 b/mlprec/mld_d_prec_type.f90 index a5f3f583..f4f594fb 100644 --- a/mlprec/mld_d_prec_type.f90 +++ b/mlprec/mld_d_prec_type.f90 @@ -744,14 +744,14 @@ contains end subroutine mld_d_apply1v - subroutine mld_d_dump(prec,info,istart,iend,prefix,head,ac,rp,smoother,solver) + subroutine mld_d_dump(prec,info,istart,iend,prefix,head,ac,rp,smoother,solver,global_num) implicit none class(mld_dprec_type), intent(in) :: prec integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: istart, iend character(len=*), intent(in), optional :: prefix, head - logical, optional, intent(in) :: smoother, solver,ac, rp + logical, optional, intent(in) :: smoother, solver,ac, rp, global_num integer(psb_ipk_) :: i, j, il1, iln, lname, lev integer(psb_ipk_) :: icontxt,iam, np character(len=80) :: prefix_ @@ -772,7 +772,7 @@ contains do lev=il1, iln call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,& - & ac=ac,smoother=smoother,solver=solver,rp=rp) + & ac=ac,smoother=smoother,solver=solver,rp=rp,global_num=global_num) end do end subroutine mld_d_dump diff --git a/mlprec/mld_s_onelev_mod.f90 b/mlprec/mld_s_onelev_mod.f90 index 3855f2e2..4cdc378b 100644 --- a/mlprec/mld_s_onelev_mod.f90 +++ b/mlprec/mld_s_onelev_mod.f90 @@ -354,7 +354,8 @@ module mld_s_onelev_mod end interface interface - subroutine mld_s_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,solver) + subroutine mld_s_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,& + & solver,global_num) import :: psb_sspmat_type, psb_s_vect_type, psb_s_base_vect_type, & & psb_slinmap_type, psb_spk_, mld_s_onelev_type, & & psb_ipk_, psb_long_int_k_, psb_desc_type @@ -363,7 +364,7 @@ module mld_s_onelev_mod integer(psb_ipk_), intent(in) :: level integer(psb_ipk_), intent(out) :: info character(len=*), intent(in), optional :: prefix, head - logical, optional, intent(in) :: ac, rp, smoother, solver + logical, optional, intent(in) :: ac, rp, smoother, solver, global_num end subroutine mld_s_base_onelev_dump end interface diff --git a/mlprec/mld_s_prec_mod.f90 b/mlprec/mld_s_prec_mod.f90 index 2a52def1..da3d546d 100644 --- a/mlprec/mld_s_prec_mod.f90 +++ b/mlprec/mld_s_prec_mod.f90 @@ -52,7 +52,7 @@ module mld_s_prec_mod use mld_s_diag_solver use mld_s_ilu_solver use mld_s_gs_solver - + interface mld_precinit subroutine mld_sprecinit(p,ptype,info,nlev) import :: psb_sspmat_type, psb_desc_type, psb_spk_, & @@ -62,13 +62,13 @@ module mld_s_prec_mod integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: nlev end subroutine mld_sprecinit - end interface + end interface mld_precinit interface mld_precset module procedure mld_s_iprecsetsm, mld_s_iprecsetsv, & & mld_s_iprecseti, mld_s_iprecsetc, mld_s_iprecsetr, & & mld_s_cprecseti, mld_s_cprecsetc, mld_s_cprecsetr - end interface + end interface mld_precset interface mld_precbld subroutine mld_sprecbld(a,desc_a,prec,info,amold,vmold,imold) @@ -83,9 +83,9 @@ module mld_s_prec_mod class(psb_s_base_sparse_mat), intent(in), optional :: amold class(psb_s_base_vect_type), intent(in), optional :: vmold class(psb_i_base_vect_type), intent(in), optional :: imold - ! character, intent(in),optional :: upd + ! character, intent(in),optional :: upd end subroutine mld_sprecbld - end interface + end interface mld_precbld interface mld_hierarchy_bld subroutine mld_s_hierarchy_bld(a,desc_a,prec,info,amold,vmold,imold) @@ -100,10 +100,30 @@ module mld_s_prec_mod class(psb_s_base_sparse_mat), intent(in), optional :: amold class(psb_s_base_vect_type), intent(in), optional :: vmold class(psb_i_base_vect_type), intent(in), optional :: imold - ! character, intent(in),optional :: upd + ! character, intent(in),optional :: upd end subroutine mld_s_hierarchy_bld - end interface + end interface mld_hierarchy_bld + + interface mld_extprol_bld + subroutine mld_s_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold) + import :: psb_sspmat_type, psb_desc_type, psb_spk_, & + & psb_s_base_sparse_mat, psb_s_base_vect_type, & + & psb_i_base_vect_type, mld_sprec_type, psb_ipk_ + ! Arguments + type(psb_sspmat_type),intent(in), target :: a + type(psb_sspmat_type),intent(inout), target :: prolv(:) + type(psb_sspmat_type),intent(inout), target :: restrv(:) + type(psb_desc_type), intent(inout), target :: desc_a + type(mld_sprec_type),intent(inout),target :: p + integer(psb_ipk_), intent(out) :: info + class(psb_s_base_sparse_mat), intent(in), optional :: amold + class(psb_s_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + ! !$ character, intent(in), optional :: upd + end subroutine mld_s_extprol_bld + end interface mld_extprol_bld + interface mld_ml_prec_bld subroutine mld_s_ml_prec_bld(a,desc_a,prec,info,amold,vmold,imold) import :: psb_sspmat_type, psb_desc_type, psb_spk_, & @@ -117,9 +137,9 @@ module mld_s_prec_mod class(psb_s_base_sparse_mat), intent(in), optional :: amold class(psb_s_base_vect_type), intent(in), optional :: vmold class(psb_i_base_vect_type), intent(in), optional :: imold - ! character, intent(in),optional :: upd + ! character, intent(in),optional :: upd end subroutine mld_s_ml_prec_bld - end interface + end interface mld_ml_prec_bld contains diff --git a/mlprec/mld_s_prec_type.f90 b/mlprec/mld_s_prec_type.f90 index 384d587d..623e1698 100644 --- a/mlprec/mld_s_prec_type.f90 +++ b/mlprec/mld_s_prec_type.f90 @@ -744,14 +744,14 @@ contains end subroutine mld_s_apply1v - subroutine mld_s_dump(prec,info,istart,iend,prefix,head,ac,rp,smoother,solver) + subroutine mld_s_dump(prec,info,istart,iend,prefix,head,ac,rp,smoother,solver,global_num) implicit none class(mld_sprec_type), intent(in) :: prec integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: istart, iend character(len=*), intent(in), optional :: prefix, head - logical, optional, intent(in) :: smoother, solver,ac, rp + logical, optional, intent(in) :: smoother, solver,ac, rp, global_num integer(psb_ipk_) :: i, j, il1, iln, lname, lev integer(psb_ipk_) :: icontxt,iam, np character(len=80) :: prefix_ @@ -772,7 +772,7 @@ contains do lev=il1, iln call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,& - & ac=ac,smoother=smoother,solver=solver,rp=rp) + & ac=ac,smoother=smoother,solver=solver,rp=rp,global_num=global_num) end do end subroutine mld_s_dump diff --git a/mlprec/mld_z_onelev_mod.f90 b/mlprec/mld_z_onelev_mod.f90 index e04fd13a..238afa08 100644 --- a/mlprec/mld_z_onelev_mod.f90 +++ b/mlprec/mld_z_onelev_mod.f90 @@ -354,7 +354,8 @@ module mld_z_onelev_mod end interface interface - subroutine mld_z_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,solver) + subroutine mld_z_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,& + & solver,global_num) import :: psb_zspmat_type, psb_z_vect_type, psb_z_base_vect_type, & & psb_zlinmap_type, psb_dpk_, mld_z_onelev_type, & & psb_ipk_, psb_long_int_k_, psb_desc_type @@ -363,7 +364,7 @@ module mld_z_onelev_mod integer(psb_ipk_), intent(in) :: level integer(psb_ipk_), intent(out) :: info character(len=*), intent(in), optional :: prefix, head - logical, optional, intent(in) :: ac, rp, smoother, solver + logical, optional, intent(in) :: ac, rp, smoother, solver, global_num end subroutine mld_z_base_onelev_dump end interface diff --git a/mlprec/mld_z_prec_mod.f90 b/mlprec/mld_z_prec_mod.f90 index 40c25ef3..cd1aa321 100644 --- a/mlprec/mld_z_prec_mod.f90 +++ b/mlprec/mld_z_prec_mod.f90 @@ -52,7 +52,7 @@ module mld_z_prec_mod use mld_z_diag_solver use mld_z_ilu_solver use mld_z_gs_solver - + interface mld_precinit subroutine mld_zprecinit(p,ptype,info,nlev) import :: psb_zspmat_type, psb_desc_type, psb_dpk_, & @@ -62,13 +62,13 @@ module mld_z_prec_mod integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: nlev end subroutine mld_zprecinit - end interface + end interface mld_precinit interface mld_precset module procedure mld_z_iprecsetsm, mld_z_iprecsetsv, & & mld_z_iprecseti, mld_z_iprecsetc, mld_z_iprecsetr, & & mld_z_cprecseti, mld_z_cprecsetc, mld_z_cprecsetr - end interface + end interface mld_precset interface mld_precbld subroutine mld_zprecbld(a,desc_a,prec,info,amold,vmold,imold) @@ -83,9 +83,9 @@ module mld_z_prec_mod class(psb_z_base_sparse_mat), intent(in), optional :: amold class(psb_z_base_vect_type), intent(in), optional :: vmold class(psb_i_base_vect_type), intent(in), optional :: imold - ! character, intent(in),optional :: upd + ! character, intent(in),optional :: upd end subroutine mld_zprecbld - end interface + end interface mld_precbld interface mld_hierarchy_bld subroutine mld_z_hierarchy_bld(a,desc_a,prec,info,amold,vmold,imold) @@ -100,10 +100,30 @@ module mld_z_prec_mod class(psb_z_base_sparse_mat), intent(in), optional :: amold class(psb_z_base_vect_type), intent(in), optional :: vmold class(psb_i_base_vect_type), intent(in), optional :: imold - ! character, intent(in),optional :: upd + ! character, intent(in),optional :: upd end subroutine mld_z_hierarchy_bld - end interface + end interface mld_hierarchy_bld + + interface mld_extprol_bld + subroutine mld_z_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold) + import :: psb_zspmat_type, psb_desc_type, psb_dpk_, & + & psb_z_base_sparse_mat, psb_z_base_vect_type, & + & psb_i_base_vect_type, mld_zprec_type, psb_ipk_ + ! Arguments + type(psb_zspmat_type),intent(in), target :: a + type(psb_zspmat_type),intent(inout), target :: prolv(:) + type(psb_zspmat_type),intent(inout), target :: restrv(:) + type(psb_desc_type), intent(inout), target :: desc_a + type(mld_zprec_type),intent(inout),target :: p + integer(psb_ipk_), intent(out) :: info + class(psb_z_base_sparse_mat), intent(in), optional :: amold + class(psb_z_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + ! !$ character, intent(in), optional :: upd + end subroutine mld_z_extprol_bld + end interface mld_extprol_bld + interface mld_ml_prec_bld subroutine mld_z_ml_prec_bld(a,desc_a,prec,info,amold,vmold,imold) import :: psb_zspmat_type, psb_desc_type, psb_dpk_, & @@ -117,9 +137,9 @@ module mld_z_prec_mod class(psb_z_base_sparse_mat), intent(in), optional :: amold class(psb_z_base_vect_type), intent(in), optional :: vmold class(psb_i_base_vect_type), intent(in), optional :: imold - ! character, intent(in),optional :: upd + ! character, intent(in),optional :: upd end subroutine mld_z_ml_prec_bld - end interface + end interface mld_ml_prec_bld contains diff --git a/mlprec/mld_z_prec_type.f90 b/mlprec/mld_z_prec_type.f90 index b563cce4..c75c6a1a 100644 --- a/mlprec/mld_z_prec_type.f90 +++ b/mlprec/mld_z_prec_type.f90 @@ -744,14 +744,14 @@ contains end subroutine mld_z_apply1v - subroutine mld_z_dump(prec,info,istart,iend,prefix,head,ac,rp,smoother,solver) + subroutine mld_z_dump(prec,info,istart,iend,prefix,head,ac,rp,smoother,solver,global_num) implicit none class(mld_zprec_type), intent(in) :: prec integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: istart, iend character(len=*), intent(in), optional :: prefix, head - logical, optional, intent(in) :: smoother, solver,ac, rp + logical, optional, intent(in) :: smoother, solver,ac, rp, global_num integer(psb_ipk_) :: i, j, il1, iln, lname, lev integer(psb_ipk_) :: icontxt,iam, np character(len=80) :: prefix_ @@ -772,7 +772,7 @@ contains do lev=il1, iln call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,& - & ac=ac,smoother=smoother,solver=solver,rp=rp) + & ac=ac,smoother=smoother,solver=solver,rp=rp,global_num=global_num) end do end subroutine mld_z_dump