diff --git a/mlprec/impl/Makefile b/mlprec/impl/Makefile index 030c3be5..4105d252 100644 --- a/mlprec/impl/Makefile +++ b/mlprec/impl/Makefile @@ -22,22 +22,22 @@ MPFOBJS=$(SMPFOBJS) $(DMPFOBJS) $(CMPFOBJS) $(ZMPFOBJS) MPCOBJS=mld_sslud_interface.o mld_dslud_interface.o mld_cslud_interface.o mld_zslud_interface.o -DINNEROBJS= mld_dcoarse_bld.o mld_dmlprec_bld.o mld_d_bld_mlhier_aggsize.o mld_d_bld_mlhier_array.o \ +DINNEROBJS= mld_dcoarse_bld.o mld_dmlprec_bld.o mld_d_bld_mlhier_aggsize.o mld_d_bld_mlhier_array.o \ mld_dilu0_fact.o mld_diluk_fact.o mld_dilut_fact.o mld_daggrmap_bld.o \ mld_d_dec_map_bld.o mld_dmlprec_aply.o mld_daggrmat_asb.o \ $(DMPFOBJS) -SINNEROBJS= mld_scoarse_bld.o mld_smlprec_bld.o \ +SINNEROBJS= mld_scoarse_bld.o mld_smlprec_bld.o mld_s_bld_mlhier_aggsize.o mld_s_bld_mlhier_array.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) -ZINNEROBJS= mld_zcoarse_bld.o mld_zmlprec_bld.o \ +ZINNEROBJS= mld_zcoarse_bld.o mld_zmlprec_bld.o mld_z_bld_mlhier_aggsize.o mld_z_bld_mlhier_array.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) -CINNEROBJS= mld_ccoarse_bld.o mld_cmlprec_bld.o \ +CINNEROBJS= mld_ccoarse_bld.o mld_cmlprec_bld.o mld_c_bld_mlhier_aggsize.o mld_c_bld_mlhier_array.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) diff --git a/mlprec/impl/level/Makefile b/mlprec/impl/level/Makefile index 87b23f0a..d11b1b81 100644 --- a/mlprec/impl/level/Makefile +++ b/mlprec/impl/level/Makefile @@ -8,7 +8,8 @@ HERE=../.. FINCLUDES=$(FMFLAG)../.. $(FMFLAG)$(INCDIR) $(FMFLAG)$(PSBINCDIR) $(FMFLAG)$(PSBLIBDIR) -OBJS=mld_c_base_onelev_check.o \ +OBJS=mld_c_base_onelev_build.o \ +mld_c_base_onelev_check.o \ mld_c_base_onelev_cnv.o \ mld_c_base_onelev_csetc.o \ mld_c_base_onelev_cseti.o \ @@ -35,6 +36,7 @@ mld_d_base_onelev_seti.o \ mld_d_base_onelev_setr.o \ mld_d_base_onelev_setsm.o \ mld_d_base_onelev_setsv.o \ +mld_s_base_onelev_build.o \ mld_s_base_onelev_check.o \ mld_s_base_onelev_cnv.o \ mld_s_base_onelev_csetc.o \ @@ -48,6 +50,7 @@ mld_s_base_onelev_seti.o \ mld_s_base_onelev_setr.o \ mld_s_base_onelev_setsm.o \ mld_s_base_onelev_setsv.o \ +mld_z_base_onelev_build.o \ mld_z_base_onelev_check.o \ mld_z_base_onelev_cnv.o \ mld_z_base_onelev_csetc.o \ diff --git a/mlprec/impl/level/mld_c_base_onelev_build.f90 b/mlprec/impl/level/mld_c_base_onelev_build.f90 new file mode 100644 index 00000000..e1f64075 --- /dev/null +++ b/mlprec/impl/level/mld_c_base_onelev_build.f90 @@ -0,0 +1,125 @@ +!!$ +!!$ +!!$ 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. +!!$ +!!$ +subroutine mld_c_base_onelev_build(lv,info,amold,vmold,imold) + use psb_base_mod + use mld_c_onelev_mod, mld_protect_name => mld_c_base_onelev_build + implicit none + class(mld_c_onelev_type), target, intent(inout) :: lv + 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 + ! Local + integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize + integer(psb_ipk_) :: ictxt, me, np + integer(psb_ipk_) :: ipv(mld_ifpsz_), val + integer(psb_ipk_) :: int_err(5) + integer(psb_ipk_) :: debug_level, debug_unit + character(len=20) :: name, ch_err + + name = 'mld_onelev_build' + 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() + if (.not.associated(lv%base_desc)) then + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Unassociated base DESC') + goto 9999 + end if + info = psb_success_ + int_err(1) = 0 + ictxt = lv%base_desc%get_ctxt() + call psb_info(ictxt,me,np) + + if (.not.allocated(lv%sm)) then + !! Error: should have called mld_dprecinit + info=3111 + call psb_errpush(info,name) + goto 9999 + end if + if (.not.allocated(lv%sm%sv)) then + !! Error: should have called mld_dprecinit + info=3111 + call psb_errpush(info,name) + goto 9999 + end if + lv%ac_nz_loc = lv%ac%get_nzeros() + lv%ac_nz_tot = lv%ac_nz_loc + call psb_sum(ictxt,lv%ac_nz_tot) + ! end do + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Calling mlprcbld at level ',i + call mld_check_def(lv%parms%sweeps,& + & 'Jacobi sweeps',ione,is_int_positive) + call mld_check_def(lv%parms%sweeps_pre,& + & 'Jacobi sweeps',ione,is_int_positive) + call mld_check_def(lv%parms%sweeps_post,& + & 'Jacobi sweeps',ione,is_int_positive) + + call lv%sm%build(lv%base_a,lv%base_desc,& + & 'F',info,amold=amold,vmold=vmold,imold=imold) + if (info == 0) then + if (allocated(lv%sm2a)) then + call lv%sm2a%build(lv%base_a,lv%base_desc,'F',info,& + & amold=amold,vmold=vmold,imold=imold) + lv%sm2 => lv%sm2a + else + lv%sm2 => lv%sm + end if + end if + if (info /=0) then + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Smoother bld error') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_c_base_onelev_build diff --git a/mlprec/impl/level/mld_s_base_onelev_build.f90 b/mlprec/impl/level/mld_s_base_onelev_build.f90 new file mode 100644 index 00000000..223bc6e9 --- /dev/null +++ b/mlprec/impl/level/mld_s_base_onelev_build.f90 @@ -0,0 +1,125 @@ +!!$ +!!$ +!!$ 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. +!!$ +!!$ +subroutine mld_s_base_onelev_build(lv,info,amold,vmold,imold) + use psb_base_mod + use mld_s_onelev_mod, mld_protect_name => mld_s_base_onelev_build + implicit none + class(mld_s_onelev_type), target, intent(inout) :: lv + 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 + ! Local + integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize + integer(psb_ipk_) :: ictxt, me, np + integer(psb_ipk_) :: ipv(mld_ifpsz_), val + integer(psb_ipk_) :: int_err(5) + integer(psb_ipk_) :: debug_level, debug_unit + character(len=20) :: name, ch_err + + name = 'mld_onelev_build' + 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() + if (.not.associated(lv%base_desc)) then + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Unassociated base DESC') + goto 9999 + end if + info = psb_success_ + int_err(1) = 0 + ictxt = lv%base_desc%get_ctxt() + call psb_info(ictxt,me,np) + + if (.not.allocated(lv%sm)) then + !! Error: should have called mld_dprecinit + info=3111 + call psb_errpush(info,name) + goto 9999 + end if + if (.not.allocated(lv%sm%sv)) then + !! Error: should have called mld_dprecinit + info=3111 + call psb_errpush(info,name) + goto 9999 + end if + lv%ac_nz_loc = lv%ac%get_nzeros() + lv%ac_nz_tot = lv%ac_nz_loc + call psb_sum(ictxt,lv%ac_nz_tot) + ! end do + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Calling mlprcbld at level ',i + call mld_check_def(lv%parms%sweeps,& + & 'Jacobi sweeps',ione,is_int_positive) + call mld_check_def(lv%parms%sweeps_pre,& + & 'Jacobi sweeps',ione,is_int_positive) + call mld_check_def(lv%parms%sweeps_post,& + & 'Jacobi sweeps',ione,is_int_positive) + + call lv%sm%build(lv%base_a,lv%base_desc,& + & 'F',info,amold=amold,vmold=vmold,imold=imold) + if (info == 0) then + if (allocated(lv%sm2a)) then + call lv%sm2a%build(lv%base_a,lv%base_desc,'F',info,& + & amold=amold,vmold=vmold,imold=imold) + lv%sm2 => lv%sm2a + else + lv%sm2 => lv%sm + end if + end if + if (info /=0) then + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Smoother bld error') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_s_base_onelev_build diff --git a/mlprec/impl/level/mld_z_base_onelev_build.f90 b/mlprec/impl/level/mld_z_base_onelev_build.f90 new file mode 100644 index 00000000..6937c690 --- /dev/null +++ b/mlprec/impl/level/mld_z_base_onelev_build.f90 @@ -0,0 +1,125 @@ +!!$ +!!$ +!!$ 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. +!!$ +!!$ +subroutine mld_z_base_onelev_build(lv,info,amold,vmold,imold) + use psb_base_mod + use mld_z_onelev_mod, mld_protect_name => mld_z_base_onelev_build + implicit none + class(mld_z_onelev_type), target, intent(inout) :: lv + 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 + ! Local + integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize + integer(psb_ipk_) :: ictxt, me, np + integer(psb_ipk_) :: ipv(mld_ifpsz_), val + integer(psb_ipk_) :: int_err(5) + integer(psb_ipk_) :: debug_level, debug_unit + character(len=20) :: name, ch_err + + name = 'mld_onelev_build' + 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() + if (.not.associated(lv%base_desc)) then + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Unassociated base DESC') + goto 9999 + end if + info = psb_success_ + int_err(1) = 0 + ictxt = lv%base_desc%get_ctxt() + call psb_info(ictxt,me,np) + + if (.not.allocated(lv%sm)) then + !! Error: should have called mld_dprecinit + info=3111 + call psb_errpush(info,name) + goto 9999 + end if + if (.not.allocated(lv%sm%sv)) then + !! Error: should have called mld_dprecinit + info=3111 + call psb_errpush(info,name) + goto 9999 + end if + lv%ac_nz_loc = lv%ac%get_nzeros() + lv%ac_nz_tot = lv%ac_nz_loc + call psb_sum(ictxt,lv%ac_nz_tot) + ! end do + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Calling mlprcbld at level ',i + call mld_check_def(lv%parms%sweeps,& + & 'Jacobi sweeps',ione,is_int_positive) + call mld_check_def(lv%parms%sweeps_pre,& + & 'Jacobi sweeps',ione,is_int_positive) + call mld_check_def(lv%parms%sweeps_post,& + & 'Jacobi sweeps',ione,is_int_positive) + + call lv%sm%build(lv%base_a,lv%base_desc,& + & 'F',info,amold=amold,vmold=vmold,imold=imold) + if (info == 0) then + if (allocated(lv%sm2a)) then + call lv%sm2a%build(lv%base_a,lv%base_desc,'F',info,& + & amold=amold,vmold=vmold,imold=imold) + lv%sm2 => lv%sm2a + else + lv%sm2 => lv%sm + end if + end if + if (info /=0) then + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Smoother bld error') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_z_base_onelev_build diff --git a/mlprec/impl/mld_c_bld_mlhier_aggsize.f90 b/mlprec/impl/mld_c_bld_mlhier_aggsize.f90 new file mode 100644 index 00000000..fde2d1a1 --- /dev/null +++ b/mlprec/impl/mld_c_bld_mlhier_aggsize.f90 @@ -0,0 +1,235 @@ +!!$ +!!$ +!!$ 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. +!!$ +!!$ +! +! Build an aggregation hierarchy with a target aggregation size +! +! +subroutine mld_c_bld_mlhier_aggsize(casize,a,desc_a,iszv,precv,info) + use psb_base_mod + use mld_c_inner_mod, mld_protect_name => mld_c_bld_mlhier_aggsize + use mld_c_prec_mod + implicit none + integer(psb_ipk_), intent(in) :: casize + integer(psb_ipk_), intent(inout) :: iszv + type(psb_cspmat_type),intent(in), target :: a + type(psb_desc_type), intent(inout), target :: desc_a + type(mld_c_onelev_type), allocatable, target, intent(inout) :: precv(:) + integer(psb_ipk_), intent(out) :: info + ! Local + integer(psb_ipk_) :: ictxt, me,np + integer(psb_ipk_) :: err,i,k, err_act, newsz + integer(psb_ipk_) :: ipv(mld_ifpsz_), val + integer(psb_ipk_) :: int_err(5) + character :: upd_ + class(mld_c_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm + type(mld_sml_parms) :: baseparms, medparms, coarseparms + type(mld_c_onelev_node), pointer :: head, tail, newnode, current + integer(psb_ipk_) :: debug_level, debug_unit + character(len=20) :: name, ch_err + name = 'mld_bld_mlhier_aggsize' + 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() + ictxt = desc_a%get_ctxt() + call psb_info(ictxt,me,np) + ! + ! New strategy to build according to coarse size. + ! + coarseparms = precv(iszv)%parms + baseparms = precv(1)%parms + medparms = precv(2)%parms + + allocate(coarse_sm, source=precv(iszv)%sm,stat=info) + if (info == psb_success_) & + & allocate(med_sm, source=precv(2)%sm,stat=info) + if (info == psb_success_) & + & allocate(base_sm, source=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 + ! + ! Replicated matrix should only ever happen at coarse level. + ! + call mld_check_def(baseparms%coarse_mat,'Coarse matrix',& + & mld_distr_mat_,is_distr_ml_coarse_mat) + ! + ! Now build a doubly linked list + ! + allocate(newnode,stat=info) + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='List start'); goto 9999 + end if + head => newnode + tail => newnode + newnode%item%base_a => a + newnode%item%base_desc => desc_a + newnode%item%parms = baseparms + newsz = 1 + current => head + list_build_loop: do + allocate(newnode,stat=info) + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='List start'); goto 9999 + end if + current%next => newnode + newnode%prev => current + newsz = newsz + 1 + newnode%item%parms = medparms + newnode%item%parms%aggr_thresh = current%item%parms%aggr_thresh/2 + call mld_coarse_bld(current%item%base_a, current%item%base_desc, & + & newnode%item,info) + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='build next level'); goto 9999 + end if + + if (newsz>2) then + if (all(current%item%map%naggr == newnode%item%map%naggr)) then + ! + ! We are not gaining anything + ! + newsz = newsz-1 + current%next => null() + call newnode%item%free(info) + if (info == psb_success_) deallocate(newnode,stat=info) + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='Deallocate at list end'); goto 9999 + end if + exit list_build_loop + end if + end if + + current => current%next + tail => current + if (sum(current%item%map%naggr) <= casize) then + ! + ! Target reached; but we may need to rebuild. + ! + exit list_build_loop + end if + end do list_build_loop + ! + ! At this point, we are at the list tail, + ! and it needs to be rebuilt in case the parms were + ! different. + ! + ! But the threshold has to be fixed before rebuliding + coarseparms%aggr_thresh = current%item%parms%aggr_thresh + current%item%parms = coarseparms + call mld_coarse_bld(current%prev%item%base_a,& + & current%prev%item%base_desc, & + & current%item,info) + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='build next level'); goto 9999 + end if + + ! + ! Ok, now allocate the output vector and fix items. + ! + do i=1,iszv + if (info == psb_success_) call precv(i)%free(info) + end do + if (info == psb_success_) deallocate(precv,stat=info) + if (info == psb_success_) allocate(precv(newsz),stat=info) + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='Reallocate precv'); goto 9999 + end if + newnode => head + do i=1, newsz + current => newnode + ! First do a move_alloc. + ! This handles the AC, DESC_AC and MAP fields + if (info == psb_success_) & + & call current%item%move_alloc(precv(i),info) + ! Now set the smoother/solver parts. + if (info == psb_success_) then + if (i ==1) then + ! This is a workaround for a bug in gfortran 4.7.2 + allocate(precv(i)%sm,source=base_sm,stat=info) + else if (i < newsz) then + allocate(precv(i)%sm,source=med_sm,stat=info) + else + allocate(precv(i)%sm,source=coarse_sm,stat=info) + end if + end if + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='list cpy'); goto 9999 + end if + if (i == 1) then + precv(i)%base_a => a + precv(i)%base_desc => desc_a + else + precv(i)%base_a => precv(i)%ac + precv(i)%base_desc => precv(i)%desc_ac + precv(i)%map%p_desc_X => precv(i-1)%base_desc + precv(i)%map%p_desc_Y => precv(i)%base_desc + end if + + newnode => current%next + deallocate(current) + end do + call base_sm%free(info) + if (info == psb_success_) call med_sm%free(info) + if (info == psb_success_) call coarse_sm%free(info) + if (info == psb_success_) deallocate(coarse_sm,med_sm,base_sm,stat=info) + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='final cleanup'); goto 9999 + end if + iszv = newsz + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_c_bld_mlhier_aggsize diff --git a/mlprec/impl/mld_c_bld_mlhier_array.f90 b/mlprec/impl/mld_c_bld_mlhier_array.f90 new file mode 100644 index 00000000..226973de --- /dev/null +++ b/mlprec/impl/mld_c_bld_mlhier_array.f90 @@ -0,0 +1,188 @@ +!!$ +!!$ +!!$ 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. +!!$ +!!$ + +subroutine mld_c_bld_mlhier_array(a,desc_a,iszv,precv,info) + use psb_base_mod + use mld_c_inner_mod, mld_protect_name => mld_c_bld_mlhier_array + use mld_c_prec_mod + implicit none + integer(psb_ipk_), intent(inout) :: iszv + type(psb_cspmat_type),intent(in), target :: a + type(psb_desc_type), intent(inout), target :: desc_a + type(mld_c_onelev_type),intent(inout), allocatable, target :: precv(:) + integer(psb_ipk_), intent(out) :: info + ! Local + integer(psb_ipk_) :: ictxt, me,np + integer(psb_ipk_) :: err,i,k, err_act, newsz + integer(psb_ipk_) :: ipv(mld_ifpsz_), val + type(mld_c_onelev_type), allocatable :: tprecv(:) + integer(psb_ipk_) :: int_err(5) + integer(psb_ipk_) :: debug_level, debug_unit + character(len=20) :: name, ch_err + name = 'mld_bld_array_hierarchy' + 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() + ictxt = desc_a%get_ctxt() + call psb_info(ictxt,me,np) + + ! + ! + ! Build the matrix and the transfer operators corresponding + ! to the remaining levels + ! + ! + ! Check on the iprcparm contents: they should be the same + ! on all processes. + ! + call psb_bcast(ictxt,precv(1)%parms) + ! + ! Finest level first; remember to fix base_a and base_desc + ! + precv(1)%base_a => a + precv(1)%base_desc => desc_a + iszv = size(precv) + + array_build_loop: do i=2, iszv + ! + ! Check on the iprcparm contents: they should be the same + ! on all processes. + ! + call psb_bcast(ictxt,precv(i)%parms) + + ! + ! Sanity checks on the parameters + ! + if (i= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Calling mlprcbld at level ',i + ! + ! Build the mapping between levels i-1 and i and the matrix + ! at level i + ! + if (info == psb_success_) call mld_coarse_bld(precv(i-1)%base_a,& + & precv(i-1)%base_desc, precv(i),info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Init upper level preconditioner') + goto 9999 + endif + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Return from ',i,' call to mlprcbld ',info + + if (i>2) then + if (all(precv(i)%map%naggr == 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 (newsz > 0) then + 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 + allocate(tprecv(newsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='prec reallocation') + goto 9999 + endif + do i=1,newsz-1 + call precv(i)%move_alloc(tprecv(i),info) + end do + call precv(iszv)%move_alloc(tprecv(newsz),info) + do i=newsz+1, iszv + call precv(i)%free(info) + end do + call move_alloc(tprecv,precv) + ! Ignore errors from transfer + info = psb_success_ + ! + ! Restart + iszv = newsz + ! Fix the pointers, but the level 1 should + ! be already OK + do i=2, iszv - 1 + precv(i)%base_a => precv(i)%ac + precv(i)%base_desc => precv(i)%desc_ac + precv(i)%map%p_desc_X => precv(i-1)%base_desc + precv(i)%map%p_desc_Y => precv(i)%base_desc + end do + + i = iszv + if (info == psb_success_) call mld_coarse_bld(precv(i-1)%base_a,& + & precv(i-1)%base_desc, precv(i),info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='coarse rebuild') + goto 9999 + endif + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_c_bld_mlhier_array diff --git a/mlprec/impl/mld_cmlprec_bld.f90 b/mlprec/impl/mld_cmlprec_bld.f90 index 7b6a3d7e..4521b09d 100644 --- a/mlprec/impl/mld_cmlprec_bld.f90 +++ b/mlprec/impl/mld_cmlprec_bld.f90 @@ -93,15 +93,11 @@ subroutine mld_cmlprec_bld(a,desc_a,p,info,amold,vmold,imold) !!$ character, intent(in), optional :: upd ! Local Variables - type(mld_cprec_type) :: t_prec integer(psb_ipk_) :: ictxt, me,np integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize integer(psb_ipk_) :: ipv(mld_ifpsz_), val integer(psb_ipk_) :: int_err(5) character :: upd_ - class(mld_c_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm - type(mld_sml_parms) :: baseparms, medparms, coarseparms - type(mld_c_onelev_node), pointer :: head, tail, newnode, current integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name, ch_err @@ -124,18 +120,18 @@ subroutine mld_cmlprec_bld(a,desc_a,p,info,amold,vmold,imold) ! ! 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 + ! !$ 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 @@ -148,9 +144,9 @@ subroutine mld_cmlprec_bld(a,desc_a,p,info,amold,vmold,imold) ! ! Check to ensure all procs have the same ! - newsz = -1 + newsz = -1 casize = p%coarse_aggr_size - iszv = size(p%precv) + iszv = size(p%precv) call psb_bcast(ictxt,iszv) call psb_bcast(ictxt,casize) if (casize > 0) then @@ -177,292 +173,20 @@ subroutine mld_cmlprec_bld(a,desc_a,p,info,amold,vmold,imold) - if (casize>0) then - ! - ! New strategy to build according to coarse size. - ! - 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 - ! - ! Replicated matrix should only ever happen at coarse level. - ! - call mld_check_def(baseparms%coarse_mat,'Coarse matrix',& - & mld_distr_mat_,is_distr_ml_coarse_mat) - ! - ! Now build a doubly linked list - ! - allocate(newnode,stat=info) - if (info /= psb_success_) then - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='List start'); goto 9999 - end if - head => newnode - tail => newnode - newnode%item%base_a => a - newnode%item%base_desc => desc_a - newnode%item%parms = baseparms - newsz = 1 - current => head - list_build_loop: do - allocate(newnode,stat=info) - if (info /= psb_success_) then - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='List start'); goto 9999 - end if - current%next => newnode - newnode%prev => current - newsz = newsz + 1 - newnode%item%parms = medparms - newnode%item%parms%aggr_thresh = current%item%parms%aggr_thresh/2 - call mld_coarse_bld(current%item%base_a, current%item%base_desc, & - & newnode%item,info) - if (info /= psb_success_) then - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='build next level'); goto 9999 - end if - - if (newsz>2) then - if (all(current%item%map%naggr == newnode%item%map%naggr)) then - ! - ! We are not gaining anything - ! - newsz = newsz-1 - current%next => null() - call newnode%item%free(info) - if (info == psb_success_) deallocate(newnode,stat=info) - if (info /= psb_success_) then - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='Deallocate at list end'); goto 9999 - end if - exit list_build_loop - end if - end if - - current => current%next - tail => current - if (sum(newnode%item%map%naggr) <= casize) then - ! - ! Target reached; but we may need to rebuild. - ! - exit list_build_loop - end if - end do list_build_loop - ! - ! At this point, we are at the list tail, - ! and it needs to be rebuilt in case the parms were - ! different. - ! - ! But the threshold has to be fixed before rebuliding - coarseparms%aggr_thresh = current%item%parms%aggr_thresh - current%item%parms = coarseparms - call mld_coarse_bld(current%prev%item%base_a,& - & current%prev%item%base_desc, & - & current%item,info) - if (info /= psb_success_) then - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='build next level'); goto 9999 - end if - + if (casize>0) then ! - ! Ok, now allocate the output vector and fix items. + ! This should become the default strategy, we specify a target aggregation size. ! - do i=1,iszv - if (info == psb_success_) call p%precv(i)%free(info) - end do - if (info == psb_success_) deallocate(p%precv,stat=info) - if (info == psb_success_) allocate(p%precv(newsz),stat=info) - if (info /= psb_success_) then - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='Reallocate precv'); goto 9999 - end if - newnode => head - do i=1, newsz - current => newnode - ! First do a move_alloc. - ! This handles the AC, DESC_AC and MAP fields - if (info == psb_success_) & - & call current%item%move_alloc(p%precv(i),info) - ! Now set the smoother/solver parts. - if (info == psb_success_) then - if (i ==1) then - ! This is a workaround for a bug in gfortran 4.7.2 - call doallc(i,p%precv,base_sm,info) - ! !$ allocate(p%precv(i)%sm,source=base_sm,stat=info) - else if (i < newsz) then - call doallc(i,p%precv,med_sm,info) - ! !$ allocate(p%precv(i)%sm,source=med_sm,stat=info) - else - call doallc(i,p%precv,coarse_sm,info) - ! !$ allocate(p%precv(i)%sm,source=coarse_sm,stat=info) - end if - end if - if (info /= psb_success_) then - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='list cpy'); goto 9999 - end if - if (i == 1) then - p%precv(i)%base_a => a - p%precv(i)%base_desc => desc_a - else - p%precv(i)%base_a => p%precv(i)%ac - p%precv(i)%base_desc => p%precv(i)%desc_ac - p%precv(i)%map%p_desc_X => p%precv(i-1)%base_desc - p%precv(i)%map%p_desc_Y => p%precv(i)%base_desc - end if - - newnode => current%next - deallocate(current) - end do - call base_sm%free(info) - if (info == psb_success_) call med_sm%free(info) - if (info == psb_success_) call coarse_sm%free(info) - if (info == psb_success_) deallocate(coarse_sm,med_sm,base_sm,stat=info) - if (info /= psb_success_) then - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='final cleanup'); goto 9999 - end if - iszv = newsz - + call mld_bld_mlhier_aggsize(casize,a,desc_a,iszv,p%precv,info) else ! - ! Default, oldstyle - ! - - ! - ! Build the matrix and the transfer operators corresponding - ! to the remaining levels - ! + ! Oldstyle with fixed number of levels. ! - ! Check on the iprcparm contents: they should be the same - ! on all processes. - ! - call psb_bcast(ictxt,p%precv(1)%parms) - ! - ! Finest level first; remember to fix base_a and base_desc - ! - p%precv(1)%base_a => a - p%precv(1)%base_desc => desc_a - - - do i=2, iszv - ! - ! Check on the iprcparm contents: they should be the same - ! on all processes. - ! - call psb_bcast(ictxt,p%precv(i)%parms) - - ! - ! Sanity checks on the parameters - ! - if (i= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Calling mlprcbld at level ',i - ! - ! Build the mapping between levels i-1 and i and the matrix - ! at level i - ! - if (info == psb_success_) call mld_coarse_bld(p%precv(i-1)%base_a,& - & p%precv(i-1)%base_desc, p%precv(i),info) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Init upper level preconditioner') - goto 9999 - endif - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Return from ',i,' call to mlprcbld ',info - - 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 - end if - end do - - if (newsz > 0) then - 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 - allocate(t_prec%precv(newsz),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='prec reallocation') - goto 9999 - endif - do i=1,newsz-1 - call p%precv(i)%move_alloc(t_prec%precv(i),info) - end do - call p%precv(iszv)%move_alloc(t_prec%precv(newsz),info) - do i=newsz+1, iszv - call p%precv(i)%free(info) - end do - call t_prec%move_alloc(p,info) - ! Ignore errors from transfer - info = psb_success_ - ! - ! Restart - iszv = newsz - ! Fix the pointers, but the level 1 should - ! be already OK - do i=2, iszv - 1 - p%precv(i)%base_a => p%precv(i)%ac - p%precv(i)%base_desc => p%precv(i)%desc_ac - p%precv(i)%map%p_desc_X => p%precv(i-1)%base_desc - p%precv(i)%map%p_desc_Y => p%precv(i)%base_desc - end do - - i = iszv - if (info == psb_success_) call mld_coarse_bld(p%precv(i-1)%base_a,& - & p%precv(i-1)%base_desc, p%precv(i),info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='coarse rebuild') - goto 9999 - endif - end if + call mld_bld_mlhier_array(a,desc_a,p%precv,info) end if - ! Fix nzeros - do i=1,iszv - p%precv(i)%ac_nz_loc = p%precv(i)%ac%get_nzeros() - p%precv(i)%ac_nz_tot = p%precv(i)%ac_nz_loc - call psb_sum(ictxt,p%precv(i)%ac_nz_tot) - end do - - ! - ! The coarse space hierarchy has been built. + ! ! Now do the preconditioner build. ! @@ -471,51 +195,20 @@ subroutine mld_cmlprec_bld(a,desc_a,p,info,amold,vmold,imold) ! ! build the base preconditioner at level i ! - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Calling mlprcbld at level ',i - call mld_check_def(p%precv(i)%parms%sweeps,& - & 'Jacobi sweeps',ione,is_int_positive) - call mld_check_def(p%precv(i)%parms%sweeps_pre,& - & 'Jacobi sweeps',ione,is_int_positive) - call mld_check_def(p%precv(i)%parms%sweeps_post,& - & 'Jacobi sweeps',ione,is_int_positive) - if (.not.allocated(p%precv(i)%sm)) then - !! Error: should have called mld_dprecinit - info=3111 - call psb_errpush(info,name) - goto 9999 - end if - if (.not.allocated(p%precv(i)%sm%sv)) then - !! Error: should have called mld_dprecinit - info=3111 - call psb_errpush(info,name) - goto 9999 - end if - - call p%precv(i)%sm%build(p%precv(i)%base_a,p%precv(i)%base_desc,& - & 'F',info,amold=amold,vmold=vmold,imold=imold) - if (info == 0) then - if (allocated(p%precv(i)%sm2a)) then - call p%precv(i)%sm2a%build(a,desc_a,upd_,info,& - & amold=amold,vmold=vmold,imold=imold) - p%precv(i)%sm2 => p%precv(i)%sm2a - else - p%precv(i)%sm2 => p%precv(i)%sm - end if - end if - + call p%precv(i)%bld(info,amold=amold,vmold=vmold,imold=imold) + if (info /= psb_success_) then + write(ch_err,'(a,i7)') 'Error @ level',i call psb_errpush(psb_err_internal_error_,name,& - & a_err='One level preconditioner build.') + & a_err=ch_err) goto 9999 endif - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Return from ',i,' call to mlprcbld ',info end do + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Exiting with',iszv,' levels' call psb_erractionrestore(err_act) return @@ -523,18 +216,5 @@ subroutine mld_cmlprec_bld(a,desc_a,p,info,amold,vmold,imold) 9999 call psb_error_handler(err_act) return - -contains - - subroutine doallc(i,v,src,info) - type(mld_c_onelev_type), intent(inout) :: v(:) - integer(psb_ipk_), intent(in) :: i - class(mld_c_base_smoother_type), intent(in) :: src - integer(psb_ipk_), intent(out) :: info - - if ((lbound(v,1)<=i).and.(i<=ubound(v,1))) then - allocate(v(i)%sm,source=src,stat=info) - end if - end subroutine doallc end subroutine mld_cmlprec_bld diff --git a/mlprec/impl/mld_dmlprec_bld.f90 b/mlprec/impl/mld_dmlprec_bld.f90 index f8b93d4c..94894ce6 100644 --- a/mlprec/impl/mld_dmlprec_bld.f90 +++ b/mlprec/impl/mld_dmlprec_bld.f90 @@ -93,7 +93,6 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info,amold,vmold,imold) !!$ character, intent(in), optional :: upd ! Local Variables - type(mld_dprec_type) :: t_prec integer(psb_ipk_) :: ictxt, me,np integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize integer(psb_ipk_) :: ipv(mld_ifpsz_), val @@ -211,8 +210,6 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info,amold,vmold,imold) & write(debug_unit,*) me,' ',trim(name),& & 'Exiting with',iszv,' levels' - - call psb_erractionrestore(err_act) return diff --git a/mlprec/impl/mld_s_bld_mlhier_aggsize.f90 b/mlprec/impl/mld_s_bld_mlhier_aggsize.f90 new file mode 100644 index 00000000..f655155a --- /dev/null +++ b/mlprec/impl/mld_s_bld_mlhier_aggsize.f90 @@ -0,0 +1,235 @@ +!!$ +!!$ +!!$ 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. +!!$ +!!$ +! +! Build an aggregation hierarchy with a target aggregation size +! +! +subroutine mld_s_bld_mlhier_aggsize(casize,a,desc_a,iszv,precv,info) + use psb_base_mod + use mld_s_inner_mod, mld_protect_name => mld_s_bld_mlhier_aggsize + use mld_s_prec_mod + implicit none + integer(psb_ipk_), intent(in) :: casize + integer(psb_ipk_), intent(inout) :: iszv + type(psb_sspmat_type),intent(in), target :: a + type(psb_desc_type), intent(inout), target :: desc_a + type(mld_s_onelev_type), allocatable, target, intent(inout) :: precv(:) + integer(psb_ipk_), intent(out) :: info + ! Local + integer(psb_ipk_) :: ictxt, me,np + integer(psb_ipk_) :: err,i,k, err_act, newsz + integer(psb_ipk_) :: ipv(mld_ifpsz_), val + integer(psb_ipk_) :: int_err(5) + character :: upd_ + class(mld_s_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm + type(mld_sml_parms) :: baseparms, medparms, coarseparms + type(mld_s_onelev_node), pointer :: head, tail, newnode, current + integer(psb_ipk_) :: debug_level, debug_unit + character(len=20) :: name, ch_err + name = 'mld_bld_mlhier_aggsize' + 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() + ictxt = desc_a%get_ctxt() + call psb_info(ictxt,me,np) + ! + ! New strategy to build according to coarse size. + ! + coarseparms = precv(iszv)%parms + baseparms = precv(1)%parms + medparms = precv(2)%parms + + allocate(coarse_sm, source=precv(iszv)%sm,stat=info) + if (info == psb_success_) & + & allocate(med_sm, source=precv(2)%sm,stat=info) + if (info == psb_success_) & + & allocate(base_sm, source=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 + ! + ! Replicated matrix should only ever happen at coarse level. + ! + call mld_check_def(baseparms%coarse_mat,'Coarse matrix',& + & mld_distr_mat_,is_distr_ml_coarse_mat) + ! + ! Now build a doubly linked list + ! + allocate(newnode,stat=info) + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='List start'); goto 9999 + end if + head => newnode + tail => newnode + newnode%item%base_a => a + newnode%item%base_desc => desc_a + newnode%item%parms = baseparms + newsz = 1 + current => head + list_build_loop: do + allocate(newnode,stat=info) + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='List start'); goto 9999 + end if + current%next => newnode + newnode%prev => current + newsz = newsz + 1 + newnode%item%parms = medparms + newnode%item%parms%aggr_thresh = current%item%parms%aggr_thresh/2 + call mld_coarse_bld(current%item%base_a, current%item%base_desc, & + & newnode%item,info) + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='build next level'); goto 9999 + end if + + if (newsz>2) then + if (all(current%item%map%naggr == newnode%item%map%naggr)) then + ! + ! We are not gaining anything + ! + newsz = newsz-1 + current%next => null() + call newnode%item%free(info) + if (info == psb_success_) deallocate(newnode,stat=info) + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='Deallocate at list end'); goto 9999 + end if + exit list_build_loop + end if + end if + + current => current%next + tail => current + if (sum(current%item%map%naggr) <= casize) then + ! + ! Target reached; but we may need to rebuild. + ! + exit list_build_loop + end if + end do list_build_loop + ! + ! At this point, we are at the list tail, + ! and it needs to be rebuilt in case the parms were + ! different. + ! + ! But the threshold has to be fixed before rebuliding + coarseparms%aggr_thresh = current%item%parms%aggr_thresh + current%item%parms = coarseparms + call mld_coarse_bld(current%prev%item%base_a,& + & current%prev%item%base_desc, & + & current%item,info) + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='build next level'); goto 9999 + end if + + ! + ! Ok, now allocate the output vector and fix items. + ! + do i=1,iszv + if (info == psb_success_) call precv(i)%free(info) + end do + if (info == psb_success_) deallocate(precv,stat=info) + if (info == psb_success_) allocate(precv(newsz),stat=info) + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='Reallocate precv'); goto 9999 + end if + newnode => head + do i=1, newsz + current => newnode + ! First do a move_alloc. + ! This handles the AC, DESC_AC and MAP fields + if (info == psb_success_) & + & call current%item%move_alloc(precv(i),info) + ! Now set the smoother/solver parts. + if (info == psb_success_) then + if (i ==1) then + ! This is a workaround for a bug in gfortran 4.7.2 + allocate(precv(i)%sm,source=base_sm,stat=info) + else if (i < newsz) then + allocate(precv(i)%sm,source=med_sm,stat=info) + else + allocate(precv(i)%sm,source=coarse_sm,stat=info) + end if + end if + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='list cpy'); goto 9999 + end if + if (i == 1) then + precv(i)%base_a => a + precv(i)%base_desc => desc_a + else + precv(i)%base_a => precv(i)%ac + precv(i)%base_desc => precv(i)%desc_ac + precv(i)%map%p_desc_X => precv(i-1)%base_desc + precv(i)%map%p_desc_Y => precv(i)%base_desc + end if + + newnode => current%next + deallocate(current) + end do + call base_sm%free(info) + if (info == psb_success_) call med_sm%free(info) + if (info == psb_success_) call coarse_sm%free(info) + if (info == psb_success_) deallocate(coarse_sm,med_sm,base_sm,stat=info) + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='final cleanup'); goto 9999 + end if + iszv = newsz + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_s_bld_mlhier_aggsize diff --git a/mlprec/impl/mld_s_bld_mlhier_array.f90 b/mlprec/impl/mld_s_bld_mlhier_array.f90 new file mode 100644 index 00000000..fd4864ce --- /dev/null +++ b/mlprec/impl/mld_s_bld_mlhier_array.f90 @@ -0,0 +1,188 @@ +!!$ +!!$ +!!$ 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. +!!$ +!!$ + +subroutine mld_s_bld_mlhier_array(a,desc_a,iszv,precv,info) + use psb_base_mod + use mld_s_inner_mod, mld_protect_name => mld_s_bld_mlhier_array + use mld_s_prec_mod + implicit none + integer(psb_ipk_), intent(inout) :: iszv + type(psb_sspmat_type),intent(in), target :: a + type(psb_desc_type), intent(inout), target :: desc_a + type(mld_s_onelev_type),intent(inout), allocatable, target :: precv(:) + integer(psb_ipk_), intent(out) :: info + ! Local + integer(psb_ipk_) :: ictxt, me,np + integer(psb_ipk_) :: err,i,k, err_act, newsz + integer(psb_ipk_) :: ipv(mld_ifpsz_), val + type(mld_s_onelev_type), allocatable :: tprecv(:) + integer(psb_ipk_) :: int_err(5) + integer(psb_ipk_) :: debug_level, debug_unit + character(len=20) :: name, ch_err + name = 'mld_bld_array_hierarchy' + 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() + ictxt = desc_a%get_ctxt() + call psb_info(ictxt,me,np) + + ! + ! + ! Build the matrix and the transfer operators corresponding + ! to the remaining levels + ! + ! + ! Check on the iprcparm contents: they should be the same + ! on all processes. + ! + call psb_bcast(ictxt,precv(1)%parms) + ! + ! Finest level first; remember to fix base_a and base_desc + ! + precv(1)%base_a => a + precv(1)%base_desc => desc_a + iszv = size(precv) + + array_build_loop: do i=2, iszv + ! + ! Check on the iprcparm contents: they should be the same + ! on all processes. + ! + call psb_bcast(ictxt,precv(i)%parms) + + ! + ! Sanity checks on the parameters + ! + if (i= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Calling mlprcbld at level ',i + ! + ! Build the mapping between levels i-1 and i and the matrix + ! at level i + ! + if (info == psb_success_) call mld_coarse_bld(precv(i-1)%base_a,& + & precv(i-1)%base_desc, precv(i),info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Init upper level preconditioner') + goto 9999 + endif + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Return from ',i,' call to mlprcbld ',info + + if (i>2) then + if (all(precv(i)%map%naggr == 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 (newsz > 0) then + 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 + allocate(tprecv(newsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='prec reallocation') + goto 9999 + endif + do i=1,newsz-1 + call precv(i)%move_alloc(tprecv(i),info) + end do + call precv(iszv)%move_alloc(tprecv(newsz),info) + do i=newsz+1, iszv + call precv(i)%free(info) + end do + call move_alloc(tprecv,precv) + ! Ignore errors from transfer + info = psb_success_ + ! + ! Restart + iszv = newsz + ! Fix the pointers, but the level 1 should + ! be already OK + do i=2, iszv - 1 + precv(i)%base_a => precv(i)%ac + precv(i)%base_desc => precv(i)%desc_ac + precv(i)%map%p_desc_X => precv(i-1)%base_desc + precv(i)%map%p_desc_Y => precv(i)%base_desc + end do + + i = iszv + if (info == psb_success_) call mld_coarse_bld(precv(i-1)%base_a,& + & precv(i-1)%base_desc, precv(i),info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='coarse rebuild') + goto 9999 + endif + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_s_bld_mlhier_array diff --git a/mlprec/impl/mld_smlprec_bld.f90 b/mlprec/impl/mld_smlprec_bld.f90 index 8880e9d9..78a159f1 100644 --- a/mlprec/impl/mld_smlprec_bld.f90 +++ b/mlprec/impl/mld_smlprec_bld.f90 @@ -93,15 +93,11 @@ subroutine mld_smlprec_bld(a,desc_a,p,info,amold,vmold,imold) !!$ character, intent(in), optional :: upd ! Local Variables - type(mld_sprec_type) :: t_prec integer(psb_ipk_) :: ictxt, me,np integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize integer(psb_ipk_) :: ipv(mld_ifpsz_), val integer(psb_ipk_) :: int_err(5) character :: upd_ - class(mld_s_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm - type(mld_sml_parms) :: baseparms, medparms, coarseparms - type(mld_s_onelev_node), pointer :: head, tail, newnode, current integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name, ch_err @@ -124,18 +120,18 @@ subroutine mld_smlprec_bld(a,desc_a,p,info,amold,vmold,imold) ! ! 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 + ! !$ 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 @@ -148,9 +144,9 @@ subroutine mld_smlprec_bld(a,desc_a,p,info,amold,vmold,imold) ! ! Check to ensure all procs have the same ! - newsz = -1 + newsz = -1 casize = p%coarse_aggr_size - iszv = size(p%precv) + iszv = size(p%precv) call psb_bcast(ictxt,iszv) call psb_bcast(ictxt,casize) if (casize > 0) then @@ -177,292 +173,20 @@ subroutine mld_smlprec_bld(a,desc_a,p,info,amold,vmold,imold) - if (casize>0) then - ! - ! New strategy to build according to coarse size. - ! - 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 - ! - ! Replicated matrix should only ever happen at coarse level. - ! - call mld_check_def(baseparms%coarse_mat,'Coarse matrix',& - & mld_distr_mat_,is_distr_ml_coarse_mat) - ! - ! Now build a doubly linked list - ! - allocate(newnode,stat=info) - if (info /= psb_success_) then - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='List start'); goto 9999 - end if - head => newnode - tail => newnode - newnode%item%base_a => a - newnode%item%base_desc => desc_a - newnode%item%parms = baseparms - newsz = 1 - current => head - list_build_loop: do - allocate(newnode,stat=info) - if (info /= psb_success_) then - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='List start'); goto 9999 - end if - current%next => newnode - newnode%prev => current - newsz = newsz + 1 - newnode%item%parms = medparms - newnode%item%parms%aggr_thresh = current%item%parms%aggr_thresh/2 - call mld_coarse_bld(current%item%base_a, current%item%base_desc, & - & newnode%item,info) - if (info /= psb_success_) then - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='build next level'); goto 9999 - end if - - if (newsz>2) then - if (all(current%item%map%naggr == newnode%item%map%naggr)) then - ! - ! We are not gaining anything - ! - newsz = newsz-1 - current%next => null() - call newnode%item%free(info) - if (info == psb_success_) deallocate(newnode,stat=info) - if (info /= psb_success_) then - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='Deallocate at list end'); goto 9999 - end if - exit list_build_loop - end if - end if - - current => current%next - tail => current - if (sum(newnode%item%map%naggr) <= casize) then - ! - ! Target reached; but we may need to rebuild. - ! - exit list_build_loop - end if - end do list_build_loop - ! - ! At this point, we are at the list tail, - ! and it needs to be rebuilt in case the parms were - ! different. - ! - ! But the threshold has to be fixed before rebuliding - coarseparms%aggr_thresh = current%item%parms%aggr_thresh - current%item%parms = coarseparms - call mld_coarse_bld(current%prev%item%base_a,& - & current%prev%item%base_desc, & - & current%item,info) - if (info /= psb_success_) then - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='build next level'); goto 9999 - end if - + if (casize>0) then ! - ! Ok, now allocate the output vector and fix items. + ! This should become the default strategy, we specify a target aggregation size. ! - do i=1,iszv - if (info == psb_success_) call p%precv(i)%free(info) - end do - if (info == psb_success_) deallocate(p%precv,stat=info) - if (info == psb_success_) allocate(p%precv(newsz),stat=info) - if (info /= psb_success_) then - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='Reallocate precv'); goto 9999 - end if - newnode => head - do i=1, newsz - current => newnode - ! First do a move_alloc. - ! This handles the AC, DESC_AC and MAP fields - if (info == psb_success_) & - & call current%item%move_alloc(p%precv(i),info) - ! Now set the smoother/solver parts. - if (info == psb_success_) then - if (i ==1) then - ! This is a workaround for a bug in gfortran 4.7.2 - call doallc(i,p%precv,base_sm,info) - ! !$ allocate(p%precv(i)%sm,source=base_sm,stat=info) - else if (i < newsz) then - call doallc(i,p%precv,med_sm,info) - ! !$ allocate(p%precv(i)%sm,source=med_sm,stat=info) - else - call doallc(i,p%precv,coarse_sm,info) - ! !$ allocate(p%precv(i)%sm,source=coarse_sm,stat=info) - end if - end if - if (info /= psb_success_) then - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='list cpy'); goto 9999 - end if - if (i == 1) then - p%precv(i)%base_a => a - p%precv(i)%base_desc => desc_a - else - p%precv(i)%base_a => p%precv(i)%ac - p%precv(i)%base_desc => p%precv(i)%desc_ac - p%precv(i)%map%p_desc_X => p%precv(i-1)%base_desc - p%precv(i)%map%p_desc_Y => p%precv(i)%base_desc - end if - - newnode => current%next - deallocate(current) - end do - call base_sm%free(info) - if (info == psb_success_) call med_sm%free(info) - if (info == psb_success_) call coarse_sm%free(info) - if (info == psb_success_) deallocate(coarse_sm,med_sm,base_sm,stat=info) - if (info /= psb_success_) then - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='final cleanup'); goto 9999 - end if - iszv = newsz - + call mld_bld_mlhier_aggsize(casize,a,desc_a,iszv,p%precv,info) else ! - ! Default, oldstyle - ! - - ! - ! Build the matrix and the transfer operators corresponding - ! to the remaining levels - ! + ! Oldstyle with fixed number of levels. ! - ! Check on the iprcparm contents: they should be the same - ! on all processes. - ! - call psb_bcast(ictxt,p%precv(1)%parms) - ! - ! Finest level first; remember to fix base_a and base_desc - ! - p%precv(1)%base_a => a - p%precv(1)%base_desc => desc_a - - - do i=2, iszv - ! - ! Check on the iprcparm contents: they should be the same - ! on all processes. - ! - call psb_bcast(ictxt,p%precv(i)%parms) - - ! - ! Sanity checks on the parameters - ! - if (i= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Calling mlprcbld at level ',i - ! - ! Build the mapping between levels i-1 and i and the matrix - ! at level i - ! - if (info == psb_success_) call mld_coarse_bld(p%precv(i-1)%base_a,& - & p%precv(i-1)%base_desc, p%precv(i),info) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Init upper level preconditioner') - goto 9999 - endif - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Return from ',i,' call to mlprcbld ',info - - 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 - end if - end do - - if (newsz > 0) then - 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 - allocate(t_prec%precv(newsz),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='prec reallocation') - goto 9999 - endif - do i=1,newsz-1 - call p%precv(i)%move_alloc(t_prec%precv(i),info) - end do - call p%precv(iszv)%move_alloc(t_prec%precv(newsz),info) - do i=newsz+1, iszv - call p%precv(i)%free(info) - end do - call t_prec%move_alloc(p,info) - ! Ignore errors from transfer - info = psb_success_ - ! - ! Restart - iszv = newsz - ! Fix the pointers, but the level 1 should - ! be already OK - do i=2, iszv - 1 - p%precv(i)%base_a => p%precv(i)%ac - p%precv(i)%base_desc => p%precv(i)%desc_ac - p%precv(i)%map%p_desc_X => p%precv(i-1)%base_desc - p%precv(i)%map%p_desc_Y => p%precv(i)%base_desc - end do - - i = iszv - if (info == psb_success_) call mld_coarse_bld(p%precv(i-1)%base_a,& - & p%precv(i-1)%base_desc, p%precv(i),info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='coarse rebuild') - goto 9999 - endif - end if + call mld_bld_mlhier_array(a,desc_a,p%precv,info) end if - ! Fix nzeros - do i=1,iszv - p%precv(i)%ac_nz_loc = p%precv(i)%ac%get_nzeros() - p%precv(i)%ac_nz_tot = p%precv(i)%ac_nz_loc - call psb_sum(ictxt,p%precv(i)%ac_nz_tot) - end do - - ! - ! The coarse space hierarchy has been built. + ! ! Now do the preconditioner build. ! @@ -471,51 +195,20 @@ subroutine mld_smlprec_bld(a,desc_a,p,info,amold,vmold,imold) ! ! build the base preconditioner at level i ! - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Calling mlprcbld at level ',i - call mld_check_def(p%precv(i)%parms%sweeps,& - & 'Jacobi sweeps',ione,is_int_positive) - call mld_check_def(p%precv(i)%parms%sweeps_pre,& - & 'Jacobi sweeps',ione,is_int_positive) - call mld_check_def(p%precv(i)%parms%sweeps_post,& - & 'Jacobi sweeps',ione,is_int_positive) - if (.not.allocated(p%precv(i)%sm)) then - !! Error: should have called mld_dprecinit - info=3111 - call psb_errpush(info,name) - goto 9999 - end if - if (.not.allocated(p%precv(i)%sm%sv)) then - !! Error: should have called mld_dprecinit - info=3111 - call psb_errpush(info,name) - goto 9999 - end if - - call p%precv(i)%sm%build(p%precv(i)%base_a,p%precv(i)%base_desc,& - & 'F',info,amold=amold,vmold=vmold,imold=imold) - if (info == 0) then - if (allocated(p%precv(i)%sm2a)) then - call p%precv(i)%sm2a%build(a,desc_a,upd_,info,& - & amold=amold,vmold=vmold,imold=imold) - p%precv(i)%sm2 => p%precv(i)%sm2a - else - p%precv(i)%sm2 => p%precv(i)%sm - end if - end if - + call p%precv(i)%bld(info,amold=amold,vmold=vmold,imold=imold) + if (info /= psb_success_) then + write(ch_err,'(a,i7)') 'Error @ level',i call psb_errpush(psb_err_internal_error_,name,& - & a_err='One level preconditioner build.') + & a_err=ch_err) goto 9999 endif - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Return from ',i,' call to mlprcbld ',info end do + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Exiting with',iszv,' levels' call psb_erractionrestore(err_act) return @@ -523,18 +216,5 @@ subroutine mld_smlprec_bld(a,desc_a,p,info,amold,vmold,imold) 9999 call psb_error_handler(err_act) return - -contains - - subroutine doallc(i,v,src,info) - type(mld_s_onelev_type), intent(inout) :: v(:) - integer(psb_ipk_), intent(in) :: i - class(mld_s_base_smoother_type), intent(in) :: src - integer(psb_ipk_), intent(out) :: info - - if ((lbound(v,1)<=i).and.(i<=ubound(v,1))) then - allocate(v(i)%sm,source=src,stat=info) - end if - end subroutine doallc end subroutine mld_smlprec_bld diff --git a/mlprec/impl/mld_z_bld_mlhier_aggsize.f90 b/mlprec/impl/mld_z_bld_mlhier_aggsize.f90 new file mode 100644 index 00000000..cce4d6e0 --- /dev/null +++ b/mlprec/impl/mld_z_bld_mlhier_aggsize.f90 @@ -0,0 +1,235 @@ +!!$ +!!$ +!!$ 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. +!!$ +!!$ +! +! Build an aggregation hierarchy with a target aggregation size +! +! +subroutine mld_z_bld_mlhier_aggsize(casize,a,desc_a,iszv,precv,info) + use psb_base_mod + use mld_z_inner_mod, mld_protect_name => mld_z_bld_mlhier_aggsize + use mld_z_prec_mod + implicit none + integer(psb_ipk_), intent(in) :: casize + integer(psb_ipk_), intent(inout) :: iszv + type(psb_zspmat_type),intent(in), target :: a + type(psb_desc_type), intent(inout), target :: desc_a + type(mld_z_onelev_type), allocatable, target, intent(inout) :: precv(:) + integer(psb_ipk_), intent(out) :: info + ! Local + integer(psb_ipk_) :: ictxt, me,np + integer(psb_ipk_) :: err,i,k, err_act, newsz + integer(psb_ipk_) :: ipv(mld_ifpsz_), val + integer(psb_ipk_) :: int_err(5) + character :: upd_ + class(mld_z_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm + type(mld_dml_parms) :: baseparms, medparms, coarseparms + type(mld_z_onelev_node), pointer :: head, tail, newnode, current + integer(psb_ipk_) :: debug_level, debug_unit + character(len=20) :: name, ch_err + name = 'mld_bld_mlhier_aggsize' + 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() + ictxt = desc_a%get_ctxt() + call psb_info(ictxt,me,np) + ! + ! New strategy to build according to coarse size. + ! + coarseparms = precv(iszv)%parms + baseparms = precv(1)%parms + medparms = precv(2)%parms + + allocate(coarse_sm, source=precv(iszv)%sm,stat=info) + if (info == psb_success_) & + & allocate(med_sm, source=precv(2)%sm,stat=info) + if (info == psb_success_) & + & allocate(base_sm, source=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 + ! + ! Replicated matrix should only ever happen at coarse level. + ! + call mld_check_def(baseparms%coarse_mat,'Coarse matrix',& + & mld_distr_mat_,is_distr_ml_coarse_mat) + ! + ! Now build a doubly linked list + ! + allocate(newnode,stat=info) + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='List start'); goto 9999 + end if + head => newnode + tail => newnode + newnode%item%base_a => a + newnode%item%base_desc => desc_a + newnode%item%parms = baseparms + newsz = 1 + current => head + list_build_loop: do + allocate(newnode,stat=info) + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='List start'); goto 9999 + end if + current%next => newnode + newnode%prev => current + newsz = newsz + 1 + newnode%item%parms = medparms + newnode%item%parms%aggr_thresh = current%item%parms%aggr_thresh/2 + call mld_coarse_bld(current%item%base_a, current%item%base_desc, & + & newnode%item,info) + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='build next level'); goto 9999 + end if + + if (newsz>2) then + if (all(current%item%map%naggr == newnode%item%map%naggr)) then + ! + ! We are not gaining anything + ! + newsz = newsz-1 + current%next => null() + call newnode%item%free(info) + if (info == psb_success_) deallocate(newnode,stat=info) + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='Deallocate at list end'); goto 9999 + end if + exit list_build_loop + end if + end if + + current => current%next + tail => current + if (sum(current%item%map%naggr) <= casize) then + ! + ! Target reached; but we may need to rebuild. + ! + exit list_build_loop + end if + end do list_build_loop + ! + ! At this point, we are at the list tail, + ! and it needs to be rebuilt in case the parms were + ! different. + ! + ! But the threshold has to be fixed before rebuliding + coarseparms%aggr_thresh = current%item%parms%aggr_thresh + current%item%parms = coarseparms + call mld_coarse_bld(current%prev%item%base_a,& + & current%prev%item%base_desc, & + & current%item,info) + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='build next level'); goto 9999 + end if + + ! + ! Ok, now allocate the output vector and fix items. + ! + do i=1,iszv + if (info == psb_success_) call precv(i)%free(info) + end do + if (info == psb_success_) deallocate(precv,stat=info) + if (info == psb_success_) allocate(precv(newsz),stat=info) + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='Reallocate precv'); goto 9999 + end if + newnode => head + do i=1, newsz + current => newnode + ! First do a move_alloc. + ! This handles the AC, DESC_AC and MAP fields + if (info == psb_success_) & + & call current%item%move_alloc(precv(i),info) + ! Now set the smoother/solver parts. + if (info == psb_success_) then + if (i ==1) then + ! This is a workaround for a bug in gfortran 4.7.2 + allocate(precv(i)%sm,source=base_sm,stat=info) + else if (i < newsz) then + allocate(precv(i)%sm,source=med_sm,stat=info) + else + allocate(precv(i)%sm,source=coarse_sm,stat=info) + end if + end if + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='list cpy'); goto 9999 + end if + if (i == 1) then + precv(i)%base_a => a + precv(i)%base_desc => desc_a + else + precv(i)%base_a => precv(i)%ac + precv(i)%base_desc => precv(i)%desc_ac + precv(i)%map%p_desc_X => precv(i-1)%base_desc + precv(i)%map%p_desc_Y => precv(i)%base_desc + end if + + newnode => current%next + deallocate(current) + end do + call base_sm%free(info) + if (info == psb_success_) call med_sm%free(info) + if (info == psb_success_) call coarse_sm%free(info) + if (info == psb_success_) deallocate(coarse_sm,med_sm,base_sm,stat=info) + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='final cleanup'); goto 9999 + end if + iszv = newsz + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_z_bld_mlhier_aggsize diff --git a/mlprec/impl/mld_z_bld_mlhier_array.f90 b/mlprec/impl/mld_z_bld_mlhier_array.f90 new file mode 100644 index 00000000..7c84d857 --- /dev/null +++ b/mlprec/impl/mld_z_bld_mlhier_array.f90 @@ -0,0 +1,188 @@ +!!$ +!!$ +!!$ 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. +!!$ +!!$ + +subroutine mld_z_bld_mlhier_array(a,desc_a,iszv,precv,info) + use psb_base_mod + use mld_z_inner_mod, mld_protect_name => mld_z_bld_mlhier_array + use mld_z_prec_mod + implicit none + integer(psb_ipk_), intent(inout) :: iszv + type(psb_zspmat_type),intent(in), target :: a + type(psb_desc_type), intent(inout), target :: desc_a + type(mld_z_onelev_type),intent(inout), allocatable, target :: precv(:) + integer(psb_ipk_), intent(out) :: info + ! Local + integer(psb_ipk_) :: ictxt, me,np + integer(psb_ipk_) :: err,i,k, err_act, newsz + integer(psb_ipk_) :: ipv(mld_ifpsz_), val + type(mld_z_onelev_type), allocatable :: tprecv(:) + integer(psb_ipk_) :: int_err(5) + integer(psb_ipk_) :: debug_level, debug_unit + character(len=20) :: name, ch_err + name = 'mld_bld_array_hierarchy' + 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() + ictxt = desc_a%get_ctxt() + call psb_info(ictxt,me,np) + + ! + ! + ! Build the matrix and the transfer operators corresponding + ! to the remaining levels + ! + ! + ! Check on the iprcparm contents: they should be the same + ! on all processes. + ! + call psb_bcast(ictxt,precv(1)%parms) + ! + ! Finest level first; remember to fix base_a and base_desc + ! + precv(1)%base_a => a + precv(1)%base_desc => desc_a + iszv = size(precv) + + array_build_loop: do i=2, iszv + ! + ! Check on the iprcparm contents: they should be the same + ! on all processes. + ! + call psb_bcast(ictxt,precv(i)%parms) + + ! + ! Sanity checks on the parameters + ! + if (i= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Calling mlprcbld at level ',i + ! + ! Build the mapping between levels i-1 and i and the matrix + ! at level i + ! + if (info == psb_success_) call mld_coarse_bld(precv(i-1)%base_a,& + & precv(i-1)%base_desc, precv(i),info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Init upper level preconditioner') + goto 9999 + endif + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Return from ',i,' call to mlprcbld ',info + + if (i>2) then + if (all(precv(i)%map%naggr == 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 (newsz > 0) then + 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 + allocate(tprecv(newsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='prec reallocation') + goto 9999 + endif + do i=1,newsz-1 + call precv(i)%move_alloc(tprecv(i),info) + end do + call precv(iszv)%move_alloc(tprecv(newsz),info) + do i=newsz+1, iszv + call precv(i)%free(info) + end do + call move_alloc(tprecv,precv) + ! Ignore errors from transfer + info = psb_success_ + ! + ! Restart + iszv = newsz + ! Fix the pointers, but the level 1 should + ! be already OK + do i=2, iszv - 1 + precv(i)%base_a => precv(i)%ac + precv(i)%base_desc => precv(i)%desc_ac + precv(i)%map%p_desc_X => precv(i-1)%base_desc + precv(i)%map%p_desc_Y => precv(i)%base_desc + end do + + i = iszv + if (info == psb_success_) call mld_coarse_bld(precv(i-1)%base_a,& + & precv(i-1)%base_desc, precv(i),info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='coarse rebuild') + goto 9999 + endif + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_z_bld_mlhier_array diff --git a/mlprec/impl/mld_zmlprec_bld.f90 b/mlprec/impl/mld_zmlprec_bld.f90 index f4212f60..f31bbf3c 100644 --- a/mlprec/impl/mld_zmlprec_bld.f90 +++ b/mlprec/impl/mld_zmlprec_bld.f90 @@ -93,15 +93,11 @@ subroutine mld_zmlprec_bld(a,desc_a,p,info,amold,vmold,imold) !!$ character, intent(in), optional :: upd ! Local Variables - type(mld_zprec_type) :: t_prec integer(psb_ipk_) :: ictxt, me,np integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize integer(psb_ipk_) :: ipv(mld_ifpsz_), val integer(psb_ipk_) :: int_err(5) character :: upd_ - class(mld_z_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm - type(mld_dml_parms) :: baseparms, medparms, coarseparms - type(mld_z_onelev_node), pointer :: head, tail, newnode, current integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name, ch_err @@ -124,18 +120,18 @@ subroutine mld_zmlprec_bld(a,desc_a,p,info,amold,vmold,imold) ! ! 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 + ! !$ 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 @@ -148,9 +144,9 @@ subroutine mld_zmlprec_bld(a,desc_a,p,info,amold,vmold,imold) ! ! Check to ensure all procs have the same ! - newsz = -1 + newsz = -1 casize = p%coarse_aggr_size - iszv = size(p%precv) + iszv = size(p%precv) call psb_bcast(ictxt,iszv) call psb_bcast(ictxt,casize) if (casize > 0) then @@ -177,292 +173,20 @@ subroutine mld_zmlprec_bld(a,desc_a,p,info,amold,vmold,imold) - if (casize>0) then - ! - ! New strategy to build according to coarse size. - ! - 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 - ! - ! Replicated matrix should only ever happen at coarse level. - ! - call mld_check_def(baseparms%coarse_mat,'Coarse matrix',& - & mld_distr_mat_,is_distr_ml_coarse_mat) - ! - ! Now build a doubly linked list - ! - allocate(newnode,stat=info) - if (info /= psb_success_) then - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='List start'); goto 9999 - end if - head => newnode - tail => newnode - newnode%item%base_a => a - newnode%item%base_desc => desc_a - newnode%item%parms = baseparms - newsz = 1 - current => head - list_build_loop: do - allocate(newnode,stat=info) - if (info /= psb_success_) then - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='List start'); goto 9999 - end if - current%next => newnode - newnode%prev => current - newsz = newsz + 1 - newnode%item%parms = medparms - newnode%item%parms%aggr_thresh = current%item%parms%aggr_thresh/2 - call mld_coarse_bld(current%item%base_a, current%item%base_desc, & - & newnode%item,info) - if (info /= psb_success_) then - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='build next level'); goto 9999 - end if - - if (newsz>2) then - if (all(current%item%map%naggr == newnode%item%map%naggr)) then - ! - ! We are not gaining anything - ! - newsz = newsz-1 - current%next => null() - call newnode%item%free(info) - if (info == psb_success_) deallocate(newnode,stat=info) - if (info /= psb_success_) then - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='Deallocate at list end'); goto 9999 - end if - exit list_build_loop - end if - end if - - current => current%next - tail => current - if (sum(newnode%item%map%naggr) <= casize) then - ! - ! Target reached; but we may need to rebuild. - ! - exit list_build_loop - end if - end do list_build_loop - ! - ! At this point, we are at the list tail, - ! and it needs to be rebuilt in case the parms were - ! different. - ! - ! But the threshold has to be fixed before rebuliding - coarseparms%aggr_thresh = current%item%parms%aggr_thresh - current%item%parms = coarseparms - call mld_coarse_bld(current%prev%item%base_a,& - & current%prev%item%base_desc, & - & current%item,info) - if (info /= psb_success_) then - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='build next level'); goto 9999 - end if - + if (casize>0) then ! - ! Ok, now allocate the output vector and fix items. + ! This should become the default strategy, we specify a target aggregation size. ! - do i=1,iszv - if (info == psb_success_) call p%precv(i)%free(info) - end do - if (info == psb_success_) deallocate(p%precv,stat=info) - if (info == psb_success_) allocate(p%precv(newsz),stat=info) - if (info /= psb_success_) then - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='Reallocate precv'); goto 9999 - end if - newnode => head - do i=1, newsz - current => newnode - ! First do a move_alloc. - ! This handles the AC, DESC_AC and MAP fields - if (info == psb_success_) & - & call current%item%move_alloc(p%precv(i),info) - ! Now set the smoother/solver parts. - if (info == psb_success_) then - if (i ==1) then - ! This is a workaround for a bug in gfortran 4.7.2 - call doallc(i,p%precv,base_sm,info) - ! !$ allocate(p%precv(i)%sm,source=base_sm,stat=info) - else if (i < newsz) then - call doallc(i,p%precv,med_sm,info) - ! !$ allocate(p%precv(i)%sm,source=med_sm,stat=info) - else - call doallc(i,p%precv,coarse_sm,info) - ! !$ allocate(p%precv(i)%sm,source=coarse_sm,stat=info) - end if - end if - if (info /= psb_success_) then - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='list cpy'); goto 9999 - end if - if (i == 1) then - p%precv(i)%base_a => a - p%precv(i)%base_desc => desc_a - else - p%precv(i)%base_a => p%precv(i)%ac - p%precv(i)%base_desc => p%precv(i)%desc_ac - p%precv(i)%map%p_desc_X => p%precv(i-1)%base_desc - p%precv(i)%map%p_desc_Y => p%precv(i)%base_desc - end if - - newnode => current%next - deallocate(current) - end do - call base_sm%free(info) - if (info == psb_success_) call med_sm%free(info) - if (info == psb_success_) call coarse_sm%free(info) - if (info == psb_success_) deallocate(coarse_sm,med_sm,base_sm,stat=info) - if (info /= psb_success_) then - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='final cleanup'); goto 9999 - end if - iszv = newsz - + call mld_bld_mlhier_aggsize(casize,a,desc_a,iszv,p%precv,info) else ! - ! Default, oldstyle - ! - - ! - ! Build the matrix and the transfer operators corresponding - ! to the remaining levels - ! + ! Oldstyle with fixed number of levels. ! - ! Check on the iprcparm contents: they should be the same - ! on all processes. - ! - call psb_bcast(ictxt,p%precv(1)%parms) - ! - ! Finest level first; remember to fix base_a and base_desc - ! - p%precv(1)%base_a => a - p%precv(1)%base_desc => desc_a - - - do i=2, iszv - ! - ! Check on the iprcparm contents: they should be the same - ! on all processes. - ! - call psb_bcast(ictxt,p%precv(i)%parms) - - ! - ! Sanity checks on the parameters - ! - if (i= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Calling mlprcbld at level ',i - ! - ! Build the mapping between levels i-1 and i and the matrix - ! at level i - ! - if (info == psb_success_) call mld_coarse_bld(p%precv(i-1)%base_a,& - & p%precv(i-1)%base_desc, p%precv(i),info) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Init upper level preconditioner') - goto 9999 - endif - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Return from ',i,' call to mlprcbld ',info - - 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 - end if - end do - - if (newsz > 0) then - 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 - allocate(t_prec%precv(newsz),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='prec reallocation') - goto 9999 - endif - do i=1,newsz-1 - call p%precv(i)%move_alloc(t_prec%precv(i),info) - end do - call p%precv(iszv)%move_alloc(t_prec%precv(newsz),info) - do i=newsz+1, iszv - call p%precv(i)%free(info) - end do - call t_prec%move_alloc(p,info) - ! Ignore errors from transfer - info = psb_success_ - ! - ! Restart - iszv = newsz - ! Fix the pointers, but the level 1 should - ! be already OK - do i=2, iszv - 1 - p%precv(i)%base_a => p%precv(i)%ac - p%precv(i)%base_desc => p%precv(i)%desc_ac - p%precv(i)%map%p_desc_X => p%precv(i-1)%base_desc - p%precv(i)%map%p_desc_Y => p%precv(i)%base_desc - end do - - i = iszv - if (info == psb_success_) call mld_coarse_bld(p%precv(i-1)%base_a,& - & p%precv(i-1)%base_desc, p%precv(i),info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='coarse rebuild') - goto 9999 - endif - end if + call mld_bld_mlhier_array(a,desc_a,p%precv,info) end if - ! Fix nzeros - do i=1,iszv - p%precv(i)%ac_nz_loc = p%precv(i)%ac%get_nzeros() - p%precv(i)%ac_nz_tot = p%precv(i)%ac_nz_loc - call psb_sum(ictxt,p%precv(i)%ac_nz_tot) - end do - - ! - ! The coarse space hierarchy has been built. + ! ! Now do the preconditioner build. ! @@ -471,51 +195,20 @@ subroutine mld_zmlprec_bld(a,desc_a,p,info,amold,vmold,imold) ! ! build the base preconditioner at level i ! - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Calling mlprcbld at level ',i - call mld_check_def(p%precv(i)%parms%sweeps,& - & 'Jacobi sweeps',ione,is_int_positive) - call mld_check_def(p%precv(i)%parms%sweeps_pre,& - & 'Jacobi sweeps',ione,is_int_positive) - call mld_check_def(p%precv(i)%parms%sweeps_post,& - & 'Jacobi sweeps',ione,is_int_positive) - if (.not.allocated(p%precv(i)%sm)) then - !! Error: should have called mld_dprecinit - info=3111 - call psb_errpush(info,name) - goto 9999 - end if - if (.not.allocated(p%precv(i)%sm%sv)) then - !! Error: should have called mld_dprecinit - info=3111 - call psb_errpush(info,name) - goto 9999 - end if - - call p%precv(i)%sm%build(p%precv(i)%base_a,p%precv(i)%base_desc,& - & 'F',info,amold=amold,vmold=vmold,imold=imold) - if (info == 0) then - if (allocated(p%precv(i)%sm2a)) then - call p%precv(i)%sm2a%build(a,desc_a,upd_,info,& - & amold=amold,vmold=vmold,imold=imold) - p%precv(i)%sm2 => p%precv(i)%sm2a - else - p%precv(i)%sm2 => p%precv(i)%sm - end if - end if - + call p%precv(i)%bld(info,amold=amold,vmold=vmold,imold=imold) + if (info /= psb_success_) then + write(ch_err,'(a,i7)') 'Error @ level',i call psb_errpush(psb_err_internal_error_,name,& - & a_err='One level preconditioner build.') + & a_err=ch_err) goto 9999 endif - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Return from ',i,' call to mlprcbld ',info end do + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Exiting with',iszv,' levels' call psb_erractionrestore(err_act) return @@ -523,18 +216,5 @@ subroutine mld_zmlprec_bld(a,desc_a,p,info,amold,vmold,imold) 9999 call psb_error_handler(err_act) return - -contains - - subroutine doallc(i,v,src,info) - type(mld_z_onelev_type), intent(inout) :: v(:) - integer(psb_ipk_), intent(in) :: i - class(mld_z_base_smoother_type), intent(in) :: src - integer(psb_ipk_), intent(out) :: info - - if ((lbound(v,1)<=i).and.(i<=ubound(v,1))) then - allocate(v(i)%sm,source=src,stat=info) - end if - end subroutine doallc end subroutine mld_zmlprec_bld diff --git a/mlprec/mld_c_inner_mod.f90 b/mlprec/mld_c_inner_mod.f90 index a1aa0253..269bb7ab 100644 --- a/mlprec/mld_c_inner_mod.f90 +++ b/mlprec/mld_c_inner_mod.f90 @@ -108,6 +108,33 @@ module mld_c_inner_mod end subroutine mld_ccoarse_bld end interface mld_coarse_bld + + interface mld_bld_mlhier_aggsize + subroutine mld_c_bld_mlhier_aggsize(casize,a,desc_a,iszv,precv,info) + use psb_base_mod, only : psb_ipk_, psb_cspmat_type, psb_desc_type + use mld_c_prec_type, only : mld_c_onelev_type + implicit none + integer(psb_ipk_), intent(in) :: casize + integer(psb_ipk_), intent(inout) :: iszv + type(psb_cspmat_type),intent(in), target :: a + type(psb_desc_type), intent(inout), target :: desc_a + type(mld_c_onelev_type), allocatable, target, intent(inout) :: precv(:) + integer(psb_ipk_), intent(out) :: info + end subroutine mld_c_bld_mlhier_aggsize + end interface mld_bld_mlhier_aggsize + + interface mld_bld_mlhier_array + subroutine mld_c_bld_mlhier_array(a,desc_a,precv,info) + use psb_base_mod, only : psb_ipk_, psb_cspmat_type, psb_desc_type + use mld_c_prec_type, only : mld_c_onelev_type + implicit none + type(psb_cspmat_type),intent(in), target :: a + type(psb_desc_type), intent(inout), target :: desc_a + type(mld_c_onelev_type), allocatable, target, intent(inout) :: precv(:) + integer(psb_ipk_), intent(out) :: info + end subroutine mld_c_bld_mlhier_array + end interface mld_bld_mlhier_array + interface mld_aggrmap_bld subroutine mld_caggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_, psb_ipk_ diff --git a/mlprec/mld_c_onelev_mod.f90 b/mlprec/mld_c_onelev_mod.f90 index 566a4df2..2c39c3d6 100644 --- a/mlprec/mld_c_onelev_mod.f90 +++ b/mlprec/mld_c_onelev_mod.f90 @@ -131,6 +131,7 @@ module mld_c_onelev_mod type(psb_desc_type), pointer :: base_desc => null() type(psb_clinmap_type) :: map contains + procedure, pass(lv) :: bld => mld_c_base_onelev_build procedure, pass(lv) :: clone => c_base_onelev_clone procedure, pass(lv) :: cnv => mld_c_base_onelev_cnv procedure, pass(lv) :: descr => mld_c_base_onelev_descr @@ -166,6 +167,20 @@ module mld_c_onelev_mod + interface + subroutine mld_c_base_onelev_build(lv,info,amold,vmold,imold) + import :: psb_c_base_sparse_mat, psb_c_base_vect_type, & + & psb_i_base_vect_type, psb_spk_, mld_c_onelev_type, & + & psb_ipk_, psb_long_int_k_, psb_desc_type + implicit none + class(mld_c_onelev_type), target, intent(inout) :: lv + 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 + end subroutine mld_c_base_onelev_build + end interface + interface subroutine mld_c_base_onelev_descr(lv,il,nl,ilmin,info,iout) import :: psb_cspmat_type, psb_c_vect_type, psb_c_base_vect_type, & diff --git a/mlprec/mld_d_inner_mod.f90 b/mlprec/mld_d_inner_mod.f90 index 3f4d4e49..668543a2 100644 --- a/mlprec/mld_d_inner_mod.f90 +++ b/mlprec/mld_d_inner_mod.f90 @@ -64,6 +64,7 @@ module mld_d_inner_mod end subroutine mld_dmlprec_bld end interface mld_mlprec_bld + interface mld_mlprec_aply subroutine mld_dmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info) use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_ @@ -94,6 +95,7 @@ module mld_d_inner_mod end subroutine mld_dmlprec_aply_vect end interface mld_mlprec_aply + interface mld_coarse_bld subroutine mld_dcoarse_bld(a,desc_a,p,info) use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_ @@ -106,6 +108,7 @@ module mld_d_inner_mod end subroutine mld_dcoarse_bld end interface mld_coarse_bld + interface mld_bld_mlhier_aggsize subroutine mld_d_bld_mlhier_aggsize(casize,a,desc_a,iszv,precv,info) use psb_base_mod, only : psb_ipk_, psb_dspmat_type, psb_desc_type @@ -146,6 +149,7 @@ module mld_d_inner_mod end subroutine mld_daggrmap_bld end interface mld_aggrmap_bld + interface mld_dec_map_bld subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_ @@ -189,6 +193,7 @@ module mld_d_inner_mod end subroutine mld_daggrmat_var_asb end interface + procedure(mld_daggrmat_var_asb) :: mld_daggrmat_nosmth_asb, & & mld_daggrmat_smth_asb, mld_daggrmat_minnrg_asb, & & mld_daggrmat_biz_asb diff --git a/mlprec/mld_d_onelev_mod.f90 b/mlprec/mld_d_onelev_mod.f90 index 875a61ed..8dc2e240 100644 --- a/mlprec/mld_d_onelev_mod.f90 +++ b/mlprec/mld_d_onelev_mod.f90 @@ -56,8 +56,8 @@ module mld_d_onelev_mod use mld_base_prec_type use mld_d_base_smoother_mod - use psb_base_mod, only : psb_dspmat_type, psb_d_base_sparse_mat,& - & psb_d_vect_type, psb_d_base_vect_type, psb_dlinmap_type, psb_dpk_, & + use psb_base_mod, only : psb_dspmat_type, psb_d_vect_type, & + & psb_d_base_vect_type, psb_dlinmap_type, psb_dpk_, & & psb_ipk_, psb_long_int_k_, psb_desc_type, psb_i_base_vect_type, & & psb_erractionsave, psb_error_handler ! @@ -165,6 +165,8 @@ module mld_d_onelev_mod & d_base_onelev_nullify, d_base_onelev_get_nzeros, & & d_base_onelev_clone, d_base_onelev_move_alloc + + interface subroutine mld_d_base_onelev_build(lv,info,amold,vmold,imold) import :: psb_d_base_sparse_mat, psb_d_base_vect_type, & @@ -178,7 +180,7 @@ module mld_d_onelev_mod class(psb_i_base_vect_type), intent(in), optional :: imold end subroutine mld_d_base_onelev_build end interface - + interface subroutine mld_d_base_onelev_descr(lv,il,nl,ilmin,info,iout) import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, & @@ -495,7 +497,7 @@ contains subroutine d_base_onelev_move_alloc(lv, b,info) - use psb_base_mod, only : psb_move_alloc + use psb_base_mod implicit none class(mld_d_onelev_type), target, intent(inout) :: lv, b integer(psb_ipk_), intent(out) :: info diff --git a/mlprec/mld_s_inner_mod.f90 b/mlprec/mld_s_inner_mod.f90 index 74cbcdfb..017496e4 100644 --- a/mlprec/mld_s_inner_mod.f90 +++ b/mlprec/mld_s_inner_mod.f90 @@ -108,6 +108,33 @@ module mld_s_inner_mod end subroutine mld_scoarse_bld end interface mld_coarse_bld + + interface mld_bld_mlhier_aggsize + subroutine mld_s_bld_mlhier_aggsize(casize,a,desc_a,iszv,precv,info) + use psb_base_mod, only : psb_ipk_, psb_sspmat_type, psb_desc_type + use mld_s_prec_type, only : mld_s_onelev_type + implicit none + integer(psb_ipk_), intent(in) :: casize + integer(psb_ipk_), intent(inout) :: iszv + type(psb_sspmat_type),intent(in), target :: a + type(psb_desc_type), intent(inout), target :: desc_a + type(mld_s_onelev_type), allocatable, target, intent(inout) :: precv(:) + integer(psb_ipk_), intent(out) :: info + end subroutine mld_s_bld_mlhier_aggsize + end interface mld_bld_mlhier_aggsize + + interface mld_bld_mlhier_array + subroutine mld_s_bld_mlhier_array(a,desc_a,precv,info) + use psb_base_mod, only : psb_ipk_, psb_sspmat_type, psb_desc_type + use mld_s_prec_type, only : mld_s_onelev_type + implicit none + type(psb_sspmat_type),intent(in), target :: a + type(psb_desc_type), intent(inout), target :: desc_a + type(mld_s_onelev_type), allocatable, target, intent(inout) :: precv(:) + integer(psb_ipk_), intent(out) :: info + end subroutine mld_s_bld_mlhier_array + end interface mld_bld_mlhier_array + interface mld_aggrmap_bld subroutine mld_saggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_, psb_ipk_ diff --git a/mlprec/mld_s_onelev_mod.f90 b/mlprec/mld_s_onelev_mod.f90 index bc47ee5c..3855f2e2 100644 --- a/mlprec/mld_s_onelev_mod.f90 +++ b/mlprec/mld_s_onelev_mod.f90 @@ -131,6 +131,7 @@ module mld_s_onelev_mod type(psb_desc_type), pointer :: base_desc => null() type(psb_slinmap_type) :: map contains + procedure, pass(lv) :: bld => mld_s_base_onelev_build procedure, pass(lv) :: clone => s_base_onelev_clone procedure, pass(lv) :: cnv => mld_s_base_onelev_cnv procedure, pass(lv) :: descr => mld_s_base_onelev_descr @@ -166,6 +167,20 @@ module mld_s_onelev_mod + interface + subroutine mld_s_base_onelev_build(lv,info,amold,vmold,imold) + import :: psb_s_base_sparse_mat, psb_s_base_vect_type, & + & psb_i_base_vect_type, psb_spk_, mld_s_onelev_type, & + & psb_ipk_, psb_long_int_k_, psb_desc_type + implicit none + class(mld_s_onelev_type), target, intent(inout) :: lv + 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 + end subroutine mld_s_base_onelev_build + end interface + interface subroutine mld_s_base_onelev_descr(lv,il,nl,ilmin,info,iout) import :: psb_sspmat_type, psb_s_vect_type, psb_s_base_vect_type, & diff --git a/mlprec/mld_z_inner_mod.f90 b/mlprec/mld_z_inner_mod.f90 index e9c3d360..a90ea2a6 100644 --- a/mlprec/mld_z_inner_mod.f90 +++ b/mlprec/mld_z_inner_mod.f90 @@ -108,6 +108,33 @@ module mld_z_inner_mod end subroutine mld_zcoarse_bld end interface mld_coarse_bld + + interface mld_bld_mlhier_aggsize + subroutine mld_z_bld_mlhier_aggsize(casize,a,desc_a,iszv,precv,info) + use psb_base_mod, only : psb_ipk_, psb_zspmat_type, psb_desc_type + use mld_z_prec_type, only : mld_z_onelev_type + implicit none + integer(psb_ipk_), intent(in) :: casize + integer(psb_ipk_), intent(inout) :: iszv + type(psb_zspmat_type),intent(in), target :: a + type(psb_desc_type), intent(inout), target :: desc_a + type(mld_z_onelev_type), allocatable, target, intent(inout) :: precv(:) + integer(psb_ipk_), intent(out) :: info + end subroutine mld_z_bld_mlhier_aggsize + end interface mld_bld_mlhier_aggsize + + interface mld_bld_mlhier_array + subroutine mld_z_bld_mlhier_array(a,desc_a,precv,info) + use psb_base_mod, only : psb_ipk_, psb_zspmat_type, psb_desc_type + use mld_z_prec_type, only : mld_z_onelev_type + implicit none + type(psb_zspmat_type),intent(in), target :: a + type(psb_desc_type), intent(inout), target :: desc_a + type(mld_z_onelev_type), allocatable, target, intent(inout) :: precv(:) + integer(psb_ipk_), intent(out) :: info + end subroutine mld_z_bld_mlhier_array + end interface mld_bld_mlhier_array + interface mld_aggrmap_bld subroutine mld_zaggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_, psb_ipk_ diff --git a/mlprec/mld_z_onelev_mod.f90 b/mlprec/mld_z_onelev_mod.f90 index 105cf8a5..e04fd13a 100644 --- a/mlprec/mld_z_onelev_mod.f90 +++ b/mlprec/mld_z_onelev_mod.f90 @@ -131,6 +131,7 @@ module mld_z_onelev_mod type(psb_desc_type), pointer :: base_desc => null() type(psb_zlinmap_type) :: map contains + procedure, pass(lv) :: bld => mld_z_base_onelev_build procedure, pass(lv) :: clone => z_base_onelev_clone procedure, pass(lv) :: cnv => mld_z_base_onelev_cnv procedure, pass(lv) :: descr => mld_z_base_onelev_descr @@ -166,6 +167,20 @@ module mld_z_onelev_mod + interface + subroutine mld_z_base_onelev_build(lv,info,amold,vmold,imold) + import :: psb_z_base_sparse_mat, psb_z_base_vect_type, & + & psb_i_base_vect_type, psb_dpk_, mld_z_onelev_type, & + & psb_ipk_, psb_long_int_k_, psb_desc_type + implicit none + class(mld_z_onelev_type), target, intent(inout) :: lv + 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 + end subroutine mld_z_base_onelev_build + end interface + interface subroutine mld_z_base_onelev_descr(lv,il,nl,ilmin,info,iout) import :: psb_zspmat_type, psb_z_vect_type, psb_z_base_vect_type, &