diff --git a/mlprec/Makefile b/mlprec/Makefile index 19f9f922..2f757ce9 100644 --- a/mlprec/Makefile +++ b/mlprec/Makefile @@ -6,7 +6,9 @@ HERE=. FINCLUDES=$(FMFLAG). $(FMFLAG)$(LIBDIR) $(FMFLAG)$(PSBLIBDIR) -MODOBJS=mld_prec_type.o mld_prec_mod.o mld_inner_mod.o mld_move_alloc_mod.o +MODOBJS=mld_base_prec_type.o \ + mld_s_prec_type.o mld_d_prec_type.o mld_c_prec_type.o mld_z_prec_type.o \ + mld_prec_type.o mld_prec_mod.o mld_inner_mod.o mld_move_alloc_mod.o MPFOBJS=mld_daggrmat_nosmth_asb.o mld_daggrmat_smth_asb.o mld_daggrmat_minnrg_asb.o MPCOBJS=mld_sslud_interface.o mld_dslud_interface.o mld_cslud_interface.o mld_zslud_interface.o INNEROBJS= mld_dcoarse_bld.o \ @@ -85,6 +87,8 @@ lib: mpobjs $(OBJS) /bin/cp -p $(LIBMOD) $(LOCAL_MODS) $(LIBDIR) $(F90OBJS) $(MPFOBJS): $(MODOBJS:.o=$(.mod)) +mld_s_prec_type.o mld_d_prec_type.o mld_c_prec_type.o mld_z_prec_type.o : mld_base_prec_type.o +mld_prec_type.o: mld_s_prec_type.o mld_d_prec_type.o mld_c_prec_type.o mld_z_prec_type.o mld_prec_mod.o mld_innner_mod.o: mld_prec_type.o mld_inner_mod.o: mld_move_alloc_mod.o mld_move_alloc_mod.o: mld_prec_type.o diff --git a/mlprec/mld_base_prec_type.f90 b/mlprec/mld_base_prec_type.f90 new file mode 100644 index 00000000..6f76b428 --- /dev/null +++ b/mlprec/mld_base_prec_type.f90 @@ -0,0 +1,851 @@ +!!$ +!!$ +!!$ MLD2P4 version 1.1 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 2.3.1) +!!$ +!!$ (C) Copyright 2008,2009 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +! File: mld_prec_type.f90 +! +! Module: mld_prec_type +! +! This module defines: +! - the mld_prec_type data structure containing the preconditioner and related +! data structures; +! - integer constants defining the preconditioner; +! - character constants describing the preconditioner (used by the routines +! printing out a preconditioner description); +! - the interfaces to the routines for the management of the preconditioner +! data structure (see below). +! +! It contains routines for +! - converting character constants defining the preconditioner into integer +! constants; +! - checking if the preconditioner is correctly defined; +! - printing a description of the preconditioner; +! - deallocating the preconditioner data structure. +! + +module mld_base_prec_type + + ! + ! This reduces the size of .mod file. Without the ONLY clause compilation + ! blows up on some systems. + ! + use psb_base_mod, only :& + & psb_d_sparse_mat, psb_z_sparse_mat,& + & psb_s_sparse_mat, psb_c_sparse_mat,& + & psb_desc_type,& + & psb_slinmap_type, psb_dlinmap_type,& + & psb_clinmap_type, psb_zlinmap_type, & + & psb_dpk_, psb_spk_, psb_long_int_k_, & + & psb_spfree, psb_cdfree, psb_halo_, psb_none_, psb_sum_, psb_avg_, & + & psb_nohalo_, psb_square_root_, psb_toupper, psb_root_,& + & psb_sizeof_int, psb_sizeof_long_int, psb_sizeof_sp, psb_sizeof_dp, psb_sizeof,& + & psb_cd_get_context, psb_info + use psb_prec_mod, only: psb_sprec_type, psb_dprec_type,& + & psb_cprec_type, psb_zprec_type, psb_d_base_prec_type + + + type mld_aux_onelev_map_type + integer :: naggr + integer, allocatable :: ilaggr(:) + end type mld_aux_onelev_map_type + type mld_aux_map_type + type(mld_aux_onelev_map_type), allocatable :: mapv(:) + end type mld_aux_map_type + + + ! + ! Entries in iprcparm + ! + ! These are in baseprec + ! + integer, parameter :: mld_smoother_type_ = 1 + integer, parameter :: mld_sub_solve_ = 2 + integer, parameter :: mld_sub_restr_ = 3 + integer, parameter :: mld_sub_prol_ = 4 + integer, parameter :: mld_sub_ren_ = 5 + integer, parameter :: mld_sub_ovr_ = 6 + integer, parameter :: mld_sub_fillin_ = 8 + integer, parameter :: mld_smoother_sweeps_ = 9 + !! 2 ints for 64 bit versions + integer, parameter :: mld_slu_ptr_ = 10 + integer, parameter :: mld_umf_symptr_ = 12 + integer, parameter :: mld_umf_numptr_ = 14 + integer, parameter :: mld_slud_ptr_ = 16 + integer, parameter :: mld_prec_status_ = 18 + ! + ! These are in onelev + ! + integer, parameter :: mld_ml_type_ = 20 + integer, parameter :: mld_smoother_pos_ = 21 + integer, parameter :: mld_aggr_kind_ = 22 + integer, parameter :: mld_aggr_alg_ = 23 + integer, parameter :: mld_aggr_omega_alg_ = 24 + integer, parameter :: mld_aggr_eig_ = 25 + integer, parameter :: mld_aggr_filter_ = 26 + integer, parameter :: mld_coarse_mat_ = 27 + integer, parameter :: mld_coarse_solve_ = 28 + integer, parameter :: mld_coarse_sweeps_ = 29 + integer, parameter :: mld_coarse_fillin_ = 30 + integer, parameter :: mld_coarse_subsolve_ = 31 + integer, parameter :: mld_ifpsz_ = 32 + + ! + ! Legal values for entry: mld_smoother_type_ + ! + integer, parameter :: mld_min_prec_=0, mld_noprec_=0, mld_diag_=1, mld_bjac_=2,& + & mld_as_=3, mld_pjac_=4, mld_max_prec_=4 + ! VERY IMPORTANT: we are relying on the following to be true: + ! mld_pjac_ == mld_diag_scale_ + ! mld_bjac_ == mld_milu_n_ (or mld_ilu_n_ would be fine) + ! mld_diag_scale_ < min(mld_slu_ mld_umf_, mld_sludist_) + ! This is a quick&dirty fix, but I have nothing better now... + ! + ! Legal values for entry: mld_sub_solve_ + ! + integer, parameter :: mld_f_none_=0,mld_ilu_n_=1,mld_milu_n_=2, mld_ilu_t_=3 + integer, parameter :: mld_diag_scale_=4, mld_slu_=5, mld_umf_=6, mld_sludist_=7 + integer, parameter :: mld_max_sub_solve_= 7 + ! + ! Legal values for entry: mld_sub_ren_ + ! + integer, parameter :: mld_renum_none_=0, mld_renum_glb_=1, mld_renum_gps_=2 + ! For the time being we are disabling GPS renumbering. + integer, parameter :: mld_max_renum_=1 + ! + ! Legal values for entry: mld_ml_type_ + ! + integer, parameter :: mld_no_ml_=0, mld_add_ml_=1, mld_mult_ml_=2 + integer, parameter :: mld_new_ml_prec_=3, mld_max_ml_type_=mld_mult_ml_ + ! + ! Legal values for entry: mld_smoother_pos_ + ! + integer, parameter :: mld_pre_smooth_=1, mld_post_smooth_=2,& + & mld_twoside_smooth_=3, mld_max_smooth_=mld_twoside_smooth_ + ! + ! Legal values for entry: mld_aggr_kind_ + ! + integer, parameter :: mld_no_smooth_=0, mld_smooth_prol_=1 + integer, parameter :: mld_min_energy_=2, mld_biz_prol_=3 + ! Disabling biz_prol for the time being. + integer, parameter :: mld_max_aggr_kind_=mld_min_energy_ + ! + ! Legal values for entry: mld_aggr_filter_ + ! + integer, parameter :: mld_no_filter_mat_=0, mld_filter_mat_=1 + integer, parameter :: mld_max_filter_mat_=mld_no_filter_mat_ + ! + ! Legal values for entry: mld_aggr_alg_ + ! + integer, parameter :: mld_dec_aggr_=0, mld_sym_dec_aggr_=1 + integer, parameter :: mld_glb_aggr_=2, mld_new_dec_aggr_=3 + integer, parameter :: mld_new_glb_aggr_=4 + integer, parameter :: mld_max_aggr_alg_=mld_dec_aggr_ + + ! + ! Legal values for entry: mld_aggr_omega_alg_ + ! + integer, parameter :: mld_eig_est_=0, mld_user_choice_=999 + ! + ! Legal values for entry: mld_aggr_eig_ + ! + integer, parameter :: mld_max_norm_=0 + ! + ! Legal values for entry: mld_coarse_mat_ + ! + integer, parameter :: mld_distr_mat_=0, mld_repl_mat_=1 + integer, parameter :: mld_max_coarse_mat_=mld_repl_mat_ + ! + ! Legal values for entry: mld_prec_status_ + ! + integer, parameter :: mld_prec_built_=98765 + + ! + ! Entries in rprcparm: ILU(k,t) threshold, smoothed aggregation omega + ! + integer, parameter :: mld_sub_iluthrs_ = 1 + integer, parameter :: mld_aggr_omega_val_ = 2 + integer, parameter :: mld_aggr_thresh_ = 3 + integer, parameter :: mld_coarse_iluthrs_ = 4 + integer, parameter :: mld_rfpsz_ = 8 + + ! + ! Fields for sparse matrices ensembles stored in av() + ! + integer, parameter :: mld_l_pr_=1, mld_u_pr_=2, mld_bp_ilu_avsz_=2 + integer, parameter :: mld_ap_nd_=3, mld_ac_=4, mld_sm_pr_t_=5, mld_sm_pr_=6 + integer, parameter :: mld_smth_avsz_=6, mld_max_avsz_=mld_smth_avsz_ + + ! + ! Character constants used by mld_file_prec_descr + ! + character(len=19), parameter, private :: & + & eigen_estimates(0:0)=(/'infinity norm '/) + character(len=19), parameter, private :: & + & smooth_names(1:3)=(/'pre-smoothing ','post-smoothing ',& + & 'pre/post-smoothing'/) + character(len=15), parameter, private :: & + & aggr_kinds(0:3)=(/'nonsmoothed ','smoothed ',& + & 'min energy ','bizr. smoothed'/) + character(len=15), parameter, private :: & + & matrix_names(0:1)=(/'distributed ','replicated '/) + character(len=18), parameter, private :: & + & aggr_names(0:4)=(/'local aggregation ','sym. local aggr. ',& + & 'global aggregation', 'new local aggr. ','new global aggr. '/) + character(len=6), parameter, private :: & + & restrict_names(0:4)=(/'none ','halo ',' ',' ',' '/) + character(len=12), parameter, private :: & + & prolong_names(0:3)=(/'none ','sum ','average ','square root'/) + character(len=15), parameter, private :: & + & ml_names(0:3)=(/'none ','additive ','multiplicative',& + & 'new ML '/) + character(len=15), parameter, private :: & + & fact_names(0:7)=(/'none ','ILU(n) ',& + & 'MILU(n) ','ILU(t,n) ',& + & 'SuperLU ','UMFPACK LU ',& + & 'SuperLU_Dist ','DiagSc-PntJac '/) + + interface mld_check_def + module procedure mld_icheck_def, mld_scheck_def, mld_dcheck_def + end interface + +contains + + ! + ! Subroutine: mld_stringval + ! + ! This routine converts the string contained into string into the corresponding + ! integer value. + ! + ! Arguments: + ! string - character(len=*), input. + ! The string to be converted. + ! val - integer, output. + ! The integer value corresponding to the string + ! info - integer, output. + ! Error code. + ! + subroutine mld_stringval(string,val,info) + implicit none + ! Arguments + character(len=*), intent(in) :: string + integer, intent(out) :: val, info + character(len=*), parameter :: name='mld_stringval' + + info = 0 + select case(psb_toupper(trim(string))) + case('NONE') + val = 0 + case('HALO') + val = psb_halo_ + case('SUM') + val = psb_sum_ + case('AVG') + val = psb_avg_ + case('FACT_NONE') + val = mld_f_none_ + case('ILU') + val = mld_ilu_n_ + case('MILU') + val = mld_milu_n_ + case('ILUT') + val = mld_ilu_t_ + case('UMF') + val = mld_umf_ + case('SLU') + val = mld_slu_ + case('SLUDIST') + val = mld_sludist_ + case('DSCALE') + val = mld_diag_scale_ + case('ADD') + val = mld_add_ml_ + case('MULT') + val = mld_mult_ml_ + case('DEC') + val = mld_dec_aggr_ + case('SYMDEC') + val = mld_sym_dec_aggr_ + case('GLB') + val = mld_glb_aggr_ + case('REPL') + val = mld_repl_mat_ + case('DIST') + val = mld_distr_mat_ + case('NONSMOOTHED') + val = mld_no_smooth_ + case('SMOOTHED') + val = mld_smooth_prol_ + case('MINENERGY') + val = mld_min_energy_ + case('PRE') + val = mld_pre_smooth_ + case('POST') + val = mld_post_smooth_ + case('TWOSIDE') + val = mld_twoside_smooth_ + case('NOPREC') + val = mld_noprec_ + case('DIAG') + val = mld_diag_ + case('BJAC') + val = mld_bjac_ + case('PJAC') + val = mld_pjac_ + case('AS') + val = mld_as_ + case('RENUM_NONE') + val = mld_renum_none_ + case('RENUM_GLOBAL') + val = mld_renum_glb_ + case('RENUM_GPS') + val = mld_renum_gps_ + case('A_NORMI') + val = mld_max_norm_ + case('USER_CHOICE') + val = mld_user_choice_ + case('EIG_EST') + val = mld_eig_est_ + case('FILTER') + val = mld_filter_mat_ + case('NO_FILTER') + val = mld_no_filter_mat_ + case default + val = -1 + info = -1 + end select + if (info /= 0) then + write(0,*) name,': Error: unknown request: "',trim(string),'"' + end if + end subroutine mld_stringval + + + ! + ! Routines printing out a description of the preconditioner + ! + + subroutine mld_base_prec_descr(iout,iprcparm, info,rprcparm,dprcparm) + implicit none + integer, intent(in) :: iprcparm(:),iout + integer, intent(out) :: info + real(psb_spk_), intent(in), optional :: rprcparm(:) + real(psb_dpk_), intent(in), optional :: dprcparm(:) + + info = 0 + if (count((/ present(rprcparm),present(dprcparm) /)) /= 1) then + info=581 +!!$ call psb_errpush(info,name,a_err=" rprcparm, dprcparm") + return + endif + + select case(iprcparm(mld_smoother_type_)) + case(mld_noprec_) + write(iout,*) ' No preconditioning' + case(mld_diag_) + write(iout,*) ' Diagonal scaling' + case(mld_pjac_) + write(iout,*) ' Point Jacobi ' + case(mld_bjac_) + write(iout,*) ' Block Jacobi with ',& + & fact_names(iprcparm(mld_sub_solve_)) + select case(iprcparm(mld_sub_solve_)) + case(mld_ilu_n_,mld_milu_n_) + write(iout,*) ' Fill level:',iprcparm(mld_sub_fillin_) + case(mld_ilu_t_) + write(iout,*) ' Fill level:',iprcparm(mld_sub_fillin_) + if (present(rprcparm)) then + write(iout,*) ' Fill threshold :',rprcparm(mld_sub_iluthrs_) + else + write(iout,*) ' Fill threshold :',dprcparm(mld_sub_iluthrs_) + end if + case(mld_slu_,mld_umf_,mld_sludist_,mld_diag_scale_) + case default + write(iout,*) ' Should never get here!' + end select + case(mld_as_) + write(iout,*) ' Additive Schwarz with ',& + & fact_names(iprcparm(mld_sub_solve_)) + select case(iprcparm(mld_sub_solve_)) + case(mld_ilu_n_,mld_milu_n_) + write(iout,*) ' Fill level:',iprcparm(mld_sub_fillin_) + case(mld_ilu_t_) + write(iout,*) ' Fill level:',iprcparm(mld_sub_fillin_) + if (present(rprcparm)) then + write(iout,*) ' Fill threshold :',rprcparm(mld_sub_iluthrs_) + else + write(iout,*) ' Fill threshold :',dprcparm(mld_sub_iluthrs_) + end if + case(mld_slu_,mld_umf_,mld_sludist_,mld_diag_scale_) + case default + write(iout,*) ' Should never get here!' + end select + write(iout,*) ' Overlap:',& + & iprcparm(mld_sub_ovr_) + write(iout,*) ' Restriction: ',& + & restrict_names(iprcparm(mld_sub_restr_)) + write(iout,*) ' Prolongation: ',& + & prolong_names(iprcparm(mld_sub_prol_)) + end select + return + end subroutine mld_base_prec_descr + + subroutine mld_ml_alg_descr(iout,ilev,iprcparm, info,rprcparm,dprcparm) + implicit none + integer, intent(in) :: iprcparm(:),iout,ilev + integer, intent(out) :: info + real(psb_spk_), intent(in), optional :: rprcparm(:) + real(psb_dpk_), intent(in), optional :: dprcparm(:) + + info = 0 + if (count((/ present(rprcparm),present(dprcparm) /)) /= 1) then + info=581 +!!$ call psb_errpush(info,name,a_err=" rprcparm, dprcparm") + return + endif + + if (iprcparm(mld_ml_type_)>mld_no_ml_) then + + write(iout,*) ' Multilevel type: ',& + & ml_names(iprcparm(mld_ml_type_)) + write(iout,*) ' Smoother position: ',& + & smooth_names(iprcparm(mld_smoother_pos_)) + write(iout,*) ' Aggregation: ', & + & aggr_names(iprcparm(mld_aggr_alg_)) + write(iout,*) ' Aggregation type: ', & + & aggr_kinds(iprcparm(mld_aggr_kind_)) + if (present(rprcparm)) then + write(iout,*) ' Aggregation threshold: ', & + & rprcparm(mld_aggr_thresh_) + else + write(iout,*) ' Aggregation threshold: ', & + & dprcparm(mld_aggr_thresh_) + end if + if (iprcparm(mld_aggr_kind_) /= mld_no_smooth_) then + if (iprcparm(mld_aggr_omega_alg_) == mld_eig_est_) then + write(iout,*) ' Damping omega computation: spectral radius estimate' + write(iout,*) ' Spectral radius estimate: ', & + & eigen_estimates(iprcparm(mld_aggr_eig_)) + else if (iprcparm(mld_aggr_omega_alg_) == mld_user_choice_) then + write(iout,*) ' Damping omega computation: user defined value.' + else + write(iout,*) ' Damping omega computation: unknown value in iprcparm!!' + end if + end if + end if + + return + end subroutine mld_ml_alg_descr + + subroutine mld_ml_level_descr(iout,ilev,iprcparm,nlaggr, info,rprcparm,dprcparm) + implicit none + integer, intent(in) :: iprcparm(:),iout,ilev + integer, intent(in), allocatable :: nlaggr(:) + integer, intent(out) :: info + real(psb_spk_), intent(in), optional :: rprcparm(:) + real(psb_dpk_), intent(in), optional :: dprcparm(:) + + info = 0 + if (count((/ present(rprcparm),present(dprcparm) /)) /= 1) then + info=581 +!!$ call psb_errpush(info,name,a_err=" rprcparm, dprcparm") + return + endif + + if (iprcparm(mld_ml_type_)>mld_no_ml_) then + write(iout,*) ' Level ',ilev + if (allocated(nlaggr)) then + write(iout,*) ' Size of coarse matrix: ', & + & sum(nlaggr(:)) + write(iout,*) ' Sizes of aggregates: ', & + & nlaggr(:) + end if + if (iprcparm(mld_aggr_kind_) /= mld_no_smooth_) then + if (present(rprcparm)) then + write(iout,*) ' Damping omega: ', & + & rprcparm(mld_aggr_omega_val_) + else + write(iout,*) ' Damping omega: ', & + & dprcparm(mld_aggr_omega_val_) + end if + end if + end if + + return + end subroutine mld_ml_level_descr + + subroutine mld_ml_coarse_descr(iout,ilev,iprcparm,iprcparm2,nlaggr,info,& + & rprcparm,dprcparm, rprcparm2,dprcparm2) + implicit none + integer, intent(in) :: iprcparm(:),iprcparm2(:),iout,ilev + integer, intent(in), allocatable :: nlaggr(:) + integer, intent(out) :: info + real(psb_spk_), intent(in), optional :: rprcparm(:), rprcparm2(:) + real(psb_dpk_), intent(in), optional :: dprcparm(:), dprcparm2(:) + + info = 0 + if (count((/ present(rprcparm),present(dprcparm) /)) /= 1) then + info=581 +!!$ call psb_errpush(info,name,a_err=" rprcparm, dprcparm") + return + endif + if (count((/ present(rprcparm2),present(dprcparm2) /)) /= 1) then + info=581 +!!$ call psb_errpush(info,name,a_err=" rprcparm, dprcparm") + return + endif + + if (iprcparm(mld_ml_type_)>mld_no_ml_) then + + write(iout,*) ' Level ',ilev,' (coarsest)' + write(iout,*) ' Coarsest matrix: ',& + & matrix_names(iprcparm(mld_coarse_mat_)) + if (allocated(nlaggr)) then + write(iout,*) ' Size of coarsest matrix: ', & + & sum( nlaggr(:)) + write(iout,*) ' Sizes of aggregates: ', & + & nlaggr(:) + end if + if (iprcparm(mld_aggr_kind_) /= mld_no_smooth_) then + + if (present(rprcparm)) then + write(iout,*) ' Damping omega: ', & + & rprcparm(mld_aggr_omega_val_) + else + write(iout,*) ' Damping omega: ', & + & dprcparm(mld_aggr_omega_val_) + end if + end if + if (iprcparm(mld_coarse_mat_) == mld_distr_mat_ .and. & + & iprcparm(mld_sub_solve_) /= mld_sludist_) then + write(iout,*) ' Coarsest matrix solver: block Jacobi with ', & + & fact_names(iprcparm2(mld_sub_solve_)) + write(iout,*) ' Number of Jacobi sweeps: ', & + & (iprcparm2(mld_smoother_sweeps_)) + else + write(iout,*) ' Coarsest matrix solver: ', & + & fact_names(iprcparm2(mld_sub_solve_)) + end if + select case(iprcparm2(mld_sub_solve_)) + case(mld_ilu_n_,mld_milu_n_) + write(iout,*) ' Fill level:',iprcparm2(mld_sub_fillin_) + case(mld_ilu_t_) + write(iout,*) ' Fill level:',iprcparm2(mld_sub_fillin_) + if (present(rprcparm2)) then + write(iout,*) ' Fill threshold :',rprcparm2(mld_sub_iluthrs_) + else if (present(dprcparm2)) then + write(iout,*) ' Fill threshold :',dprcparm2(mld_sub_iluthrs_) + end if + case(mld_slu_,mld_umf_,mld_sludist_,mld_diag_scale_) + case default + write(iout,*) ' Should never get here!' + end select + end if + + + return + end subroutine mld_ml_coarse_descr + + + ! + ! Functions/subroutines checking if the preconditioner is correctly defined + ! + + function is_legal_base_prec(ip) + implicit none + integer, intent(in) :: ip + logical :: is_legal_base_prec + + is_legal_base_prec = ((ip>=mld_noprec_).and.(ip<=mld_max_prec_)) + return + end function is_legal_base_prec + function is_legal_n_ovr(ip) + implicit none + integer, intent(in) :: ip + logical :: is_legal_n_ovr + + is_legal_n_ovr = (ip >= 0) + return + end function is_legal_n_ovr + function is_legal_renum(ip) + implicit none + integer, intent(in) :: ip + logical :: is_legal_renum + is_legal_renum = ((ip >= 0).and.(ip <= mld_max_renum_)) + return + end function is_legal_renum + function is_legal_jac_sweeps(ip) + implicit none + integer, intent(in) :: ip + logical :: is_legal_jac_sweeps + + is_legal_jac_sweeps = (ip >= 1) + return + end function is_legal_jac_sweeps + function is_legal_prolong(ip) + implicit none + integer, intent(in) :: ip + logical :: is_legal_prolong + is_legal_prolong = ((ip>=psb_none_).and.(ip<=psb_square_root_)) + return + end function is_legal_prolong + function is_legal_restrict(ip) + implicit none + integer, intent(in) :: ip + logical :: is_legal_restrict + is_legal_restrict = ((ip==psb_nohalo_).or.(ip==psb_halo_)) + return + end function is_legal_restrict + function is_legal_ml_type(ip) + implicit none + integer, intent(in) :: ip + logical :: is_legal_ml_type + + is_legal_ml_type = ((ip>=mld_no_ml_).and.(ip<=mld_max_ml_type_)) + return + end function is_legal_ml_type + function is_legal_ml_aggr_alg(ip) + implicit none + integer, intent(in) :: ip + logical :: is_legal_ml_aggr_alg + + is_legal_ml_aggr_alg = ((ip>=mld_dec_aggr_).and.(ip<=mld_max_aggr_alg_)) + return + end function is_legal_ml_aggr_alg + function is_legal_ml_aggr_omega_alg(ip) + implicit none + integer, intent(in) :: ip + logical :: is_legal_ml_aggr_omega_alg + + is_legal_ml_aggr_omega_alg = ((ip==mld_eig_est_).or.(ip==mld_user_choice_)) + return + end function is_legal_ml_aggr_omega_alg + function is_legal_ml_aggr_eig(ip) + implicit none + integer, intent(in) :: ip + logical :: is_legal_ml_aggr_eig + + is_legal_ml_aggr_eig = (ip==mld_max_norm_) + return + end function is_legal_ml_aggr_eig + function is_legal_ml_smooth_pos(ip) + implicit none + integer, intent(in) :: ip + logical :: is_legal_ml_smooth_pos + + is_legal_ml_smooth_pos = ((ip>=mld_pre_smooth_).and.(ip<=mld_max_smooth_)) + return + end function is_legal_ml_smooth_pos + function is_legal_ml_aggr_kind(ip) + implicit none + integer, intent(in) :: ip + logical :: is_legal_ml_aggr_kind + + is_legal_ml_aggr_kind = ((ip>=0).and.(ip<=mld_max_aggr_kind_)) + return + end function is_legal_ml_aggr_kind + function is_legal_ml_coarse_mat(ip) + implicit none + integer, intent(in) :: ip + logical :: is_legal_ml_coarse_mat + + is_legal_ml_coarse_mat = ((ip>=0).and.(ip<=mld_max_coarse_mat_)) + return + end function is_legal_ml_coarse_mat + function is_legal_aggr_filter(ip) + implicit none + integer, intent(in) :: ip + logical :: is_legal_aggr_filter + + is_legal_aggr_filter = ((ip>=0).and.(ip<=mld_max_filter_mat_)) + return + end function is_legal_aggr_filter + function is_distr_ml_coarse_mat(ip) + implicit none + integer, intent(in) :: ip + logical :: is_distr_ml_coarse_mat + + is_distr_ml_coarse_mat = (ip==mld_distr_mat_) + return + end function is_distr_ml_coarse_mat + function is_legal_ml_fact(ip) + implicit none + integer, intent(in) :: ip + logical :: is_legal_ml_fact + ! Here the minimum is really 1, mld_fact_none_ is not acceptable. + is_legal_ml_fact = ((ip>=1).and.(ip<=mld_max_sub_solve_)) + return + end function is_legal_ml_fact + function is_legal_ml_lev(ip) + implicit none + integer, intent(in) :: ip + logical :: is_legal_ml_lev + + is_legal_ml_lev = (ip >= 0) + return + end function is_legal_ml_lev + function is_legal_omega(ip) + implicit none + real(psb_dpk_), intent(in) :: ip + logical :: is_legal_omega + is_legal_omega = ((ip>=0.0d0).and.(ip<=2.0d0)) + return + end function is_legal_omega + function is_legal_fact_thrs(ip) + implicit none + real(psb_dpk_), intent(in) :: ip + logical :: is_legal_fact_thrs + + is_legal_fact_thrs = (ip>=0.0d0) + return + end function is_legal_fact_thrs + function is_legal_aggr_thrs(ip) + implicit none + real(psb_dpk_), intent(in) :: ip + logical :: is_legal_aggr_thrs + + is_legal_aggr_thrs = (ip>=0.0d0) + return + end function is_legal_aggr_thrs + + function is_legal_s_omega(ip) + implicit none + real(psb_spk_), intent(in) :: ip + logical :: is_legal_s_omega + is_legal_s_omega = ((ip>=0.0).and.(ip<=2.0)) + return + end function is_legal_s_omega + function is_legal_s_fact_thrs(ip) + implicit none + real(psb_spk_), intent(in) :: ip + logical :: is_legal_s_fact_thrs + + is_legal_s_fact_thrs = (ip>=0.0) + return + end function is_legal_s_fact_thrs + function is_legal_s_aggr_thrs(ip) + implicit none + real(psb_spk_), intent(in) :: ip + logical :: is_legal_s_aggr_thrs + + is_legal_s_aggr_thrs = (ip>=0.0) + return + end function is_legal_s_aggr_thrs + + + subroutine mld_icheck_def(ip,name,id,is_legal) + implicit none + integer, intent(inout) :: ip + integer, intent(in) :: id + character(len=*), intent(in) :: name + interface + function is_legal(i) + integer, intent(in) :: i + logical :: is_legal + end function is_legal + end interface + character(len=20), parameter :: rname='mld_check_def' + + if (.not.is_legal(ip)) then + write(0,*)trim(rname),': Error: Illegal value for ',& + & name,' :',ip, '. defaulting to ',id + ip = id + end if + end subroutine mld_icheck_def + + subroutine mld_scheck_def(ip,name,id,is_legal) + implicit none + real(psb_spk_), intent(inout) :: ip + real(psb_spk_), intent(in) :: id + character(len=*), intent(in) :: name + interface + function is_legal(i) + use psb_base_mod, only : psb_spk_ + real(psb_spk_), intent(in) :: i + logical :: is_legal + end function is_legal + end interface + character(len=20), parameter :: rname='mld_check_def' + + if (.not.is_legal(ip)) then + write(0,*)trim(rname),': Error: Illegal value for ',& + & name,' :',ip, '. defaulting to ',id + ip = id + end if + end subroutine mld_scheck_def + + subroutine mld_dcheck_def(ip,name,id,is_legal) + implicit none + real(psb_dpk_), intent(inout) :: ip + real(psb_dpk_), intent(in) :: id + character(len=*), intent(in) :: name + interface + function is_legal(i) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_), intent(in) :: i + logical :: is_legal + end function is_legal + end interface + character(len=20), parameter :: rname='mld_check_def' + + if (.not.is_legal(ip)) then + write(0,*)trim(rname),': Error: Illegal value for ',& + & name,' :',ip, '. defaulting to ',id + ip = id + end if + end subroutine mld_dcheck_def + + + function pr_to_str(iprec) + implicit none + + integer, intent(in) :: iprec + character(len=10) :: pr_to_str + + select case(iprec) + case(mld_noprec_) + pr_to_str='NOPREC' + case(mld_diag_) + pr_to_str='DIAG' + case(mld_bjac_) + pr_to_str='BJAC' + case(mld_as_) + pr_to_str='AS' + end select + + end function pr_to_str + + +end module mld_base_prec_type diff --git a/mlprec/mld_c_prec_type.f90 b/mlprec/mld_c_prec_type.f90 new file mode 100644 index 00000000..ba5e6a29 --- /dev/null +++ b/mlprec/mld_c_prec_type.f90 @@ -0,0 +1,688 @@ +!!$ +!!$ +!!$ MLD2P4 version 1.1 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 2.3.1) +!!$ +!!$ (C) Copyright 2008,2009 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +! File: mld_prec_type.f90 +! +! Module: mld_prec_type +! +! This module defines: +! - the mld_prec_type data structure containing the preconditioner and related +! data structures; +! - integer constants defining the preconditioner; +! - character constants describing the preconditioner (used by the routines +! printing out a preconditioner description); +! - the interfaces to the routines for the management of the preconditioner +! data structure (see below). +! +! It contains routines for +! - converting character constants defining the preconditioner into integer +! constants; +! - checking if the preconditioner is correctly defined; +! - printing a description of the preconditioner; +! - deallocating the preconditioner data structure. +! + +module mld_c_prec_type + + use mld_base_prec_type + + ! + ! Type: mld_Tprec_type. + ! + ! It is the data type containing all the information about the multilevel + ! preconditioner (here and in the following 'T' denotes 'd', 's', 'c' and + ! 'z', according to the real/complex, single/double precision version of + ! MLD2P4). It consists of an array of 'one-level' intermediate data structures + ! of type mld_Tonelev_type, each containing the information needed to apply + ! the smoothing and the coarse-space correction at a generic level. + ! + ! type mld_Tprec_type + ! type(mld_Tonelev_type), allocatable :: precv(:) + ! end type mld_Tprec_type + ! + ! Note that the levels are numbered in increasing order starting from + ! the finest one and the number of levels is given by size(precv(:)). + ! + ! + ! Type: mld_Tonelev_type. + ! + ! It is the data type containing the necessary items for the current + ! level (essentially, the base preconditioner, the current-level matrix + ! and the restriction and prolongation operators). + ! + ! type mld_Tonelev_type + ! type(mld_Tbaseprec_type) :: prec + ! integer, allocatable :: iprcparm(:) + ! real(psb_Tpk_), allocatable :: rprcparm(:) + ! type(psb_T_sparse_mat) :: ac + ! type(psb_desc_type) :: desc_ac + ! type(psb_T_sparse_mat), pointer :: base_a => null() + ! type(psb_desc_type), pointer :: base_desc => null() + ! type(psb_Tlinmap_type) :: map + ! end type mld_Tonelev_type + ! + ! Note that psb_Tpk denotes the kind of the real data type to be chosen + ! according to single/double precision version of MLD2P4. + ! + ! prec - type(mld_Tbaseprec_type). + ! The current level preconditioner (aka smoother). + ! iprcparm - integer, dimension(:), allocatable. + ! The integer parameters defining the multilevel strategy. + ! rprcparm - real(psb_Ypk_), dimension(:), allocatable. + ! The real parameters defining the multilevel strategy. + ! ac - The local part of the current-level matrix, built by + ! coarsening the previous-level matrix. + ! desc_ac - type(psb_desc_type). + ! The communication descriptor associated to the matrix + ! stored in ac. + ! base_a - type(psb_z_sparse_mat), pointer. + ! Pointer (really a pointer!) to the local part of the current + ! matrix (so we have a unified treatment of residuals). + ! We need this to avoid passing explicitly the current matrix + ! to the routine which applies the preconditioner. + ! base_desc - type(psb_desc_type), pointer. + ! Pointer to the communication descriptor associated to the + ! matrix pointed by base_a. + ! map - Stores the maps (restriction and prolongation) between the + ! vector spaces associated to the index spaces of the previous + ! and current levels. + ! + ! + ! Type: mld_Tbaseprec_type. + ! + ! It holds the smoother (base preconditioner) at a single level. + ! + ! type mld_Tbaseprec_type + ! type(psb_T_sparse_mat), allocatable :: av(:) + ! IntrType(psb_Tpk_), allocatable :: d(:) + ! type(psb_desc_type) :: desc_data + ! integer, allocatable :: iprcparm(:) + ! real(psb_Tpk_), allocatable :: rprcparm(:) + ! integer, allocatable :: perm(:), invperm(:) + ! end type mld_sbaseprec_type + ! + ! Note that IntrType denotes the real or complex data type, and psb_Tpk denotes + ! the kind of the real or complex type, according to the real/complex, single/double + ! precision version of MLD2P4. + ! + ! av - type(psb_T_sparse_mat), dimension(:), allocatable(:). + ! The sparse matrices needed to apply the preconditioner at + ! the current level ilev. + ! av(mld_l_pr_) - The L factor of the ILU factorization of the local + ! diagonal block of the current-level matrix A(ilev). + ! av(mld_u_pr_) - The U factor of the ILU factorization of the local + ! diagonal block of A(ilev), except its diagonal entries + ! (stored in d). + ! av(mld_ap_nd_) - The entries of the local part of A(ilev) outside + ! the diagonal block, for block-Jacobi sweeps. + ! d - real/complex(psb_Tpk_), dimension(:), allocatable. + ! The diagonal entries of the U factor in the ILU factorization + ! of A(ilev). + ! desc_data - type(psb_desc_type). + ! The communication descriptor associated to the base preconditioner, + ! i.e. to the sparse matrices needed to apply the base preconditioner + ! at the current level. + ! iprcparm - integer, dimension(:), allocatable. + ! The integer parameters defining the base preconditioner K(ilev) + ! (the iprcparm entries and values are specified below). + ! rprcparm - real(psb_Tpk_), dimension(:), allocatable. + ! The real parameters defining the base preconditioner K(ilev) + ! (the rprcparm entries and values are specified below). + ! perm - integer, dimension(:), allocatable. + ! The row and column permutations applied to the local part of + ! A(ilev) (defined only if iprcparm(mld_sub_ren_)>0). + ! invperm - integer, dimension(:), allocatable. + ! The inverse of the permutation stored in perm. + ! + ! Note that when the LU factorization of the (local part of the) matrix A(ilev) is + ! computed instead of the ILU one, by using UMFPACK, SuperLU or SuperLU_dist, the + ! corresponding L and U factors are stored in data structures provided by those + ! packages and pointed by prec%iprcparm(mld_umf_ptr), prec%iprcparm(mld_slu_ptr) + ! or prec%iprcparm(mld_slud_ptr). + ! + + type mld_cbaseprec_type + type(psb_c_sparse_mat), allocatable :: av(:) + complex(psb_spk_), allocatable :: d(:) + type(psb_desc_type) :: desc_data + integer, allocatable :: iprcparm(:) + real(psb_spk_), allocatable :: rprcparm(:) + integer, allocatable :: perm(:), invperm(:) + end type mld_cbaseprec_type + + type mld_conelev_type + type(mld_cbaseprec_type) :: prec + integer, allocatable :: iprcparm(:) + real(psb_spk_), allocatable :: rprcparm(:) + type(psb_c_sparse_mat) :: ac + type(psb_desc_type) :: desc_ac + type(psb_c_sparse_mat), pointer :: base_a => null() + type(psb_desc_type), pointer :: base_desc => null() + type(psb_clinmap_type) :: map + end type mld_conelev_type + + type, extends(psb_cprec_type) :: mld_cprec_type + type(mld_conelev_type), allocatable :: precv(:) + contains + procedure, pass(prec) :: c_apply2v => mld_c_apply2v + procedure, pass(prec) :: c_apply1v => mld_c_apply1v + end type mld_cprec_type + + ! + ! Interfaces to routines for checking the definition of the preconditioner, + ! for printing its description and for deallocating its data structure + ! + + interface mld_precfree + module procedure mld_cbase_precfree, mld_c_onelev_precfree, mld_cprec_free + end interface + + interface mld_nullify_baseprec + module procedure mld_nullify_cbaseprec + end interface + + interface mld_nullify_onelevprec + module procedure mld_nullify_c_onelevprec + end interface + + interface mld_precdescr + module procedure mld_cfile_prec_descr + end interface + + interface mld_sizeof + module procedure mld_cprec_sizeof, mld_cbaseprec_sizeof, mld_c_onelev_prec_sizeof + end interface + + interface mld_precaply + subroutine mld_cprecaply(prec,x,y,desc_data,info,trans,work) + use psb_base_mod, only : psb_c_sparse_mat, psb_desc_type, psb_spk_ + import mld_cprec_type + type(psb_desc_type),intent(in) :: desc_data + type(mld_cprec_type), intent(in) :: prec + complex(psb_spk_),intent(in) :: x(:) + complex(psb_spk_),intent(inout) :: y(:) + integer, intent(out) :: info + character(len=1), optional :: trans + complex(psb_spk_),intent(inout), optional, target :: work(:) + end subroutine mld_cprecaply + subroutine mld_cprecaply1(prec,x,desc_data,info,trans) + use psb_base_mod, only : psb_c_sparse_mat, psb_desc_type, psb_spk_ + import mld_cprec_type + type(psb_desc_type),intent(in) :: desc_data + type(mld_cprec_type), intent(in) :: prec + complex(psb_spk_),intent(inout) :: x(:) + integer, intent(out) :: info + character(len=1), optional :: trans + end subroutine mld_cprecaply1 + end interface + + +contains + + ! + ! Function returning the size of the mld_prec_type data structure + ! + + function mld_cprec_sizeof(prec) result(val) + implicit none + type(mld_cprec_type), intent(in) :: prec + integer(psb_long_int_k_) :: val + integer :: i + val = 0 + if (allocated(prec%precv)) then + do i=1, size(prec%precv) + val = val + mld_sizeof(prec%precv(i)) + end do + end if + end function mld_cprec_sizeof + + function mld_cbaseprec_sizeof(prec) result(val) + implicit none + type(mld_cbaseprec_type), intent(in) :: prec + integer(psb_long_int_k_) :: val + integer :: i + + val = 0 + if (allocated(prec%iprcparm)) then + val = val + psb_sizeof_int * size(prec%iprcparm) + if (prec%iprcparm(mld_prec_status_) == mld_prec_built_) then + select case(prec%iprcparm(mld_sub_solve_)) + case(mld_ilu_n_,mld_ilu_t_) + ! do nothing + case(mld_slu_) + case(mld_umf_) + case(mld_sludist_) + case default + end select + + end if + end if + if (allocated(prec%rprcparm)) val = val + psb_sizeof_sp * size(prec%rprcparm) + if (allocated(prec%d)) val = val + 2 * psb_sizeof_sp * size(prec%d) + if (allocated(prec%perm)) val = val + psb_sizeof_int * size(prec%perm) + if (allocated(prec%invperm)) val = val + psb_sizeof_int * size(prec%invperm) + val = val + psb_sizeof(prec%desc_data) + if (allocated(prec%av)) then + do i=1,size(prec%av) + val = val + psb_sizeof(prec%av(i)) + end do + end if + + end function mld_cbaseprec_sizeof + + function mld_c_onelev_prec_sizeof(prec) result(val) + implicit none + type(mld_conelev_type), intent(in) :: prec + integer(psb_long_int_k_) :: val + integer :: i + + val = mld_sizeof(prec%prec) + if (allocated(prec%iprcparm)) then + val = val + psb_sizeof_int * size(prec%iprcparm) + end if + if (allocated(prec%rprcparm)) val = val + psb_sizeof_sp * size(prec%rprcparm) + val = val + psb_sizeof(prec%desc_ac) + val = val + psb_sizeof(prec%ac) + val = val + psb_sizeof(prec%map) + + end function mld_c_onelev_prec_sizeof + + ! + ! Subroutine: mld_file_prec_descr + ! Version: real + ! + ! This routine prints a description of the preconditioner to the standard + ! output or to a file. It must be called after the preconditioner has been + ! built by mld_precbld. + ! + ! Arguments: + ! p - type(mld_Tprec_type), input. + ! The preconditioner data structure to be printed out. + ! info - integer, output. + ! error code. + ! iout - integer, input, optional. + ! The id of the file where the preconditioner description + ! will be printed. If iout is not present, then the standard + ! output is condidered. + ! + subroutine mld_cfile_prec_descr(p,info,iout) + implicit none + + ! Arguments + type(mld_cprec_type), intent(in) :: p + integer, intent(out) :: info + integer, intent(in), optional :: iout + + ! Local variables + integer :: ilev, nlev + integer :: ictxt, me, np + character(len=20), parameter :: name='mld_file_prec_descr' + integer :: iout_ + + info = 0 + if (present(iout)) then + iout_ = iout + else + iout_ = 6 + end if + if (iout_ < 0) iout_ = 6 + + if (allocated(p%precv)) then + ictxt = psb_cd_get_context(p%precv(1)%prec%desc_data) + + call psb_info(ictxt,me,np) + + ! + ! The preconditioner description is printed by processor psb_root_. + ! This agrees with the fact that all the parameters defining the + ! preconditioner have the same values on all the procs (this is + ! ensured by mld_precbld). + ! + if (me==psb_root_) then + write(iout_,*) + write(iout_,*) 'Preconditioner description' + nlev = size(p%precv) + if (nlev >= 1) then + ! + ! Print description of base preconditioner + ! + + write(iout_,*) ' ' + + if (nlev > 1) then + write(iout_,*) 'Multilevel Schwarz' + write(iout_,*) + write(iout_,*) 'Base preconditioner (smoother) details' + endif + + ilev = 1 + call mld_base_prec_descr(iout_,p%precv(ilev)%prec%iprcparm,info,& + & rprcparm=p%precv(ilev)%prec%rprcparm) + + end if + + if (nlev > 1) then + + ! + ! Print multilevel details + ! + write(iout_,*) + write(iout_,*) 'Multilevel details' + + do ilev = 2, nlev + if (.not.allocated(p%precv(ilev)%iprcparm)) then + info = 3111 + write(iout_,*) ' ',name,': error: inconsistent MLPREC part, should call MLD_PRECINIT' + return + endif + end do + + write(iout_,*) ' Number of levels: ',nlev + + ! + ! Currently, all the preconditioner parameters must have the same value at levels + ! 2,...,nlev-1, hence only the values at level 2 are printed + ! + + ilev=2 + call mld_ml_alg_descr(iout_,ilev,p%precv(ilev)%iprcparm, info,& + & rprcparm=p%precv(ilev)%rprcparm) + + ! + ! Coarse matrices are different at levels 2,...,nlev-1, hence related + ! info is printed separately + ! + write(iout_,*) + + do ilev = 2, nlev-1 + call mld_ml_level_descr(iout_,ilev,p%precv(ilev)%iprcparm,& + & p%precv(ilev)%map%naggr,info,& + & rprcparm=p%precv(ilev)%rprcparm) + end do + + ! + ! Print coarsest level details + ! + + ilev = nlev + write(iout_,*) + call mld_ml_coarse_descr(iout_,ilev,& + & p%precv(ilev)%iprcparm,p%precv(ilev)%prec%iprcparm,& + & p%precv(ilev)%map%naggr,info,& + & rprcparm=p%precv(ilev)%rprcparm,& + & rprcparm2=p%precv(ilev)%prec%rprcparm) + end if + + endif + write(iout_,*) + else + write(iout_,*) trim(name), & + & ': Error: no base preconditioner available, something is wrong!' + info = -2 + return + endif + + end subroutine mld_cfile_prec_descr + + ! + ! Subroutines: mld_Tbase_precfree, mld_T_onelev_precfree, mld_Tprec_free + ! Version: real/complex + ! + ! These routines deallocate the mld_Tbaseprec_type, mld_Tonelev_type and + ! mld_Tprec_type data structures. + ! + ! Arguments: + ! p - type(mld_Tbaseprec_type/mld_Tonelev_type/mld_Tprec_type), input. + ! The data structure to be deallocated. + ! info - integer, output. + ! error code. + ! + subroutine mld_cbase_precfree(p,info) + implicit none + type(mld_cbaseprec_type), intent(inout) :: p + integer, intent(out) :: info + integer :: i + + info = 0 + + if (allocated(p%d)) then + deallocate(p%d,stat=info) + end if + + if (allocated(p%av)) then + do i=1,size(p%av) + call psb_sp_free(p%av(i),info) + if (info /= 0) then + ! Actually, we don't care here about this. + ! Just let it go. + ! return + end if + enddo + deallocate(p%av,stat=info) + + end if + if (allocated(p%desc_data%matrix_data)) & + & call psb_cdfree(p%desc_data,info) + + if (allocated(p%rprcparm)) then + deallocate(p%rprcparm,stat=info) + end if + + if (allocated(p%perm)) then + deallocate(p%perm,stat=info) + endif + + if (allocated(p%invperm)) then + deallocate(p%invperm,stat=info) + endif + + if (allocated(p%iprcparm)) then + if (p%iprcparm(mld_prec_status_) == mld_prec_built_) then + if (p%iprcparm(mld_sub_solve_)==mld_slu_) then + call mld_cslu_free(p%iprcparm(mld_slu_ptr_),info) + end if + end if + deallocate(p%iprcparm,stat=info) + end if + call mld_nullify_baseprec(p) + end subroutine mld_cbase_precfree + + subroutine mld_c_onelev_precfree(p,info) + implicit none + + type(mld_conelev_type), intent(inout) :: p + integer, intent(out) :: info + integer :: i + + info = 0 + + ! Actually we might just deallocate the top level array, except + ! for the inner UMFPACK or SLU stuff + call mld_precfree(p%prec,info) + + call psb_sp_free(p%ac,info) + if (allocated(p%desc_ac%matrix_data)) & + & call psb_cdfree(p%desc_ac,info) + + if (allocated(p%rprcparm)) then + deallocate(p%rprcparm,stat=info) + end if + ! This is a pointer to something else, must not free it here. + nullify(p%base_a) + ! This is a pointer to something else, must not free it here. + nullify(p%base_desc) + + ! + ! free explicitly map??? + ! For now thanks to allocatable semantics + ! works anyway. + ! + + call mld_nullify_onelevprec(p) + end subroutine mld_c_onelev_precfree + + subroutine mld_nullify_c_onelevprec(p) + implicit none + + type(mld_conelev_type), intent(inout) :: p + + nullify(p%base_a) + nullify(p%base_desc) + + end subroutine mld_nullify_c_onelevprec + + subroutine mld_nullify_cbaseprec(p) + implicit none + + type(mld_cbaseprec_type), intent(inout) :: p + + + end subroutine mld_nullify_cbaseprec + + subroutine mld_cprec_free(p,info) + + use psb_base_mod + + implicit none + + ! Arguments + type(mld_cprec_type), intent(inout) :: p + integer, intent(out) :: info + + ! Local variables + integer :: me,err_act,i + character(len=20) :: name + + if(psb_get_errstatus().ne.0) return + info=0 + name = 'mld_dprecfree' + call psb_erractionsave(err_act) + + me=-1 + + if (allocated(p%precv)) then + do i=1,size(p%precv) + call mld_precfree(p%precv(i),info) + end do + deallocate(p%precv) + end if + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act.eq.psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine mld_cprec_free + + subroutine mld_c_apply2v(prec,x,y,desc_data,info,trans,work) + use psb_base_mod + + type(psb_desc_type),intent(in) :: desc_data + class(mld_cprec_type), intent(in) :: prec + complex(psb_spk_),intent(in) :: x(:) + complex(psb_spk_),intent(inout) :: y(:) + integer, intent(out) :: info + character(len=1), optional :: trans + complex(psb_spk_),intent(inout), optional, target :: work(:) + Integer :: err_act + character(len=20) :: name='s_prec_apply' + + call psb_erractionsave(err_act) + + select type(prec) +!!$ type is (mld_cprec_type) +!!$ call mld_precaply(prec,x,y,desc_data,info,trans,work) + class default + info = 700 + call psb_errpush(info,name) + goto 9999 + end select + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine mld_c_apply2v + subroutine mld_c_apply1v(prec,x,desc_data,info,trans) + use psb_base_mod + + type(psb_desc_type),intent(in) :: desc_data + class(mld_cprec_type), intent(in) :: prec + complex(psb_spk_),intent(inout) :: x(:) + integer, intent(out) :: info + character(len=1), optional :: trans + Integer :: err_act + character(len=20) :: name='c_prec_apply' + + call psb_erractionsave(err_act) + + select type(prec) +!!$ type is (mld_cprec_type) +!!$ call mld_precaply(prec,x,desc_data,info,trans) + class default + info = 700 + call psb_errpush(info,name) + goto 9999 + end select + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine mld_c_apply1v + +end module mld_c_prec_type diff --git a/mlprec/mld_d_prec_type.f90 b/mlprec/mld_d_prec_type.f90 new file mode 100644 index 00000000..882d5a32 --- /dev/null +++ b/mlprec/mld_d_prec_type.f90 @@ -0,0 +1,784 @@ +!!$ +!!$ +!!$ MLD2P4 version 1.1 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 2.3.1) +!!$ +!!$ (C) Copyright 2008,2009 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +! File: mld_prec_type.f90 +! +! Module: mld_prec_type +! +! This module defines: +! - the mld_prec_type data structure containing the preconditioner and related +! data structures; +! - integer constants defining the preconditioner; +! - character constants describing the preconditioner (used by the routines +! printing out a preconditioner description); +! - the interfaces to the routines for the management of the preconditioner +! data structure (see below). +! +! It contains routines for +! - converting character constants defining the preconditioner into integer +! constants; +! - checking if the preconditioner is correctly defined; +! - printing a description of the preconditioner; +! - deallocating the preconditioner data structure. +! + +module mld_d_prec_type + + use mld_base_prec_type + ! + ! Type: mld_Tprec_type. + ! + ! It is the data type containing all the information about the multilevel + ! preconditioner (here and in the following 'T' denotes 'd', 's', 'c' and + ! 'z', according to the real/complex, single/double precision version of + ! MLD2P4). It consists of an array of 'one-level' intermediate data structures + ! of type mld_Tonelev_type, each containing the information needed to apply + ! the smoothing and the coarse-space correction at a generic level. + ! + ! type mld_Tprec_type + ! type(mld_Tonelev_type), allocatable :: precv(:) + ! end type mld_Tprec_type + ! + ! Note that the levels are numbered in increasing order starting from + ! the finest one and the number of levels is given by size(precv(:)). + ! + ! + ! Type: mld_Tonelev_type. + ! + ! It is the data type containing the necessary items for the current + ! level (essentially, the base preconditioner, the current-level matrix + ! and the restriction and prolongation operators). + ! + ! type mld_Tonelev_type + ! type(mld_Tbaseprec_type) :: prec + ! integer, allocatable :: iprcparm(:) + ! real(psb_Tpk_), allocatable :: rprcparm(:) + ! type(psb_T_sparse_mat) :: ac + ! type(psb_desc_type) :: desc_ac + ! type(psb_T_sparse_mat), pointer :: base_a => null() + ! type(psb_desc_type), pointer :: base_desc => null() + ! type(psb_Tlinmap_type) :: map + ! end type mld_Tonelev_type + ! + ! Note that psb_Tpk denotes the kind of the real data type to be chosen + ! according to single/double precision version of MLD2P4. + ! + ! prec - type(mld_Tbaseprec_type). + ! The current level preconditioner (aka smoother). + ! iprcparm - integer, dimension(:), allocatable. + ! The integer parameters defining the multilevel strategy. + ! rprcparm - real(psb_Ypk_), dimension(:), allocatable. + ! The real parameters defining the multilevel strategy. + ! ac - The local part of the current-level matrix, built by + ! coarsening the previous-level matrix. + ! desc_ac - type(psb_desc_type). + ! The communication descriptor associated to the matrix + ! stored in ac. + ! base_a - type(psb_z_sparse_mat), pointer. + ! Pointer (really a pointer!) to the local part of the current + ! matrix (so we have a unified treatment of residuals). + ! We need this to avoid passing explicitly the current matrix + ! to the routine which applies the preconditioner. + ! base_desc - type(psb_desc_type), pointer. + ! Pointer to the communication descriptor associated to the + ! matrix pointed by base_a. + ! map - Stores the maps (restriction and prolongation) between the + ! vector spaces associated to the index spaces of the previous + ! and current levels. + ! + ! + ! Type: mld_Tbaseprec_type. + ! + ! It holds the smoother (base preconditioner) at a single level. + ! + ! type mld_Tbaseprec_type + ! type(psb_T_sparse_mat), allocatable :: av(:) + ! IntrType(psb_Tpk_), allocatable :: d(:) + ! type(psb_desc_type) :: desc_data + ! integer, allocatable :: iprcparm(:) + ! real(psb_Tpk_), allocatable :: rprcparm(:) + ! integer, allocatable :: perm(:), invperm(:) + ! end type mld_sbaseprec_type + ! + ! Note that IntrType denotes the real or complex data type, and psb_Tpk denotes + ! the kind of the real or complex type, according to the real/complex, single/double + ! precision version of MLD2P4. + ! + ! av - type(psb_T_sparse_mat), dimension(:), allocatable(:). + ! The sparse matrices needed to apply the preconditioner at + ! the current level ilev. + ! av(mld_l_pr_) - The L factor of the ILU factorization of the local + ! diagonal block of the current-level matrix A(ilev). + ! av(mld_u_pr_) - The U factor of the ILU factorization of the local + ! diagonal block of A(ilev), except its diagonal entries + ! (stored in d). + ! av(mld_ap_nd_) - The entries of the local part of A(ilev) outside + ! the diagonal block, for block-Jacobi sweeps. + ! d - real/complex(psb_Tpk_), dimension(:), allocatable. + ! The diagonal entries of the U factor in the ILU factorization + ! of A(ilev). + ! desc_data - type(psb_desc_type). + ! The communication descriptor associated to the base preconditioner, + ! i.e. to the sparse matrices needed to apply the base preconditioner + ! at the current level. + ! iprcparm - integer, dimension(:), allocatable. + ! The integer parameters defining the base preconditioner K(ilev) + ! (the iprcparm entries and values are specified below). + ! rprcparm - real(psb_Tpk_), dimension(:), allocatable. + ! The real parameters defining the base preconditioner K(ilev) + ! (the rprcparm entries and values are specified below). + ! perm - integer, dimension(:), allocatable. + ! The row and column permutations applied to the local part of + ! A(ilev) (defined only if iprcparm(mld_sub_ren_)>0). + ! invperm - integer, dimension(:), allocatable. + ! The inverse of the permutation stored in perm. + ! + ! Note that when the LU factorization of the (local part of the) matrix A(ilev) is + ! computed instead of the ILU one, by using UMFPACK, SuperLU or SuperLU_dist, the + ! corresponding L and U factors are stored in data structures provided by those + ! packages and pointed by prec%iprcparm(mld_umf_ptr), prec%iprcparm(mld_slu_ptr) + ! or prec%iprcparm(mld_slud_ptr). + ! + + type mld_d_base_solver_type +!!$ contains +!!$ procedure, pass(sv) :: build => d_base_solver_bld +!!$ procedure, pass(sv) :: apply => d_base_solver_apply + end type mld_d_base_solver_type + + type mld_d_base_smoother_type + class(mld_d_base_solver_type), allocatable :: sv + contains + procedure, pass(sm) :: build => d_base_smoother_bld + procedure, pass(sm) :: apply => d_base_smoother_apply + end type mld_d_base_smoother_type + + type, extends(psb_d_base_prec_type) :: mld_dbaseprec_type + class(mld_d_base_smoother_type), allocatable :: sm + type(psb_d_sparse_mat), allocatable :: av(:) + real(psb_dpk_), allocatable :: d(:) + type(psb_desc_type) :: desc_data + integer, allocatable :: iprcparm(:) + real(psb_dpk_), allocatable :: rprcparm(:) + integer, allocatable :: perm(:), invperm(:) + end type mld_dbaseprec_type + + type mld_donelev_type + type(mld_dbaseprec_type) :: prec + integer, allocatable :: iprcparm(:) + real(psb_dpk_), allocatable :: rprcparm(:) + type(psb_d_sparse_mat) :: ac + type(psb_desc_type) :: desc_ac + type(psb_d_sparse_mat), pointer :: base_a => null() + type(psb_desc_type), pointer :: base_desc => null() + type(psb_dlinmap_type) :: map + end type mld_donelev_type + + type, extends(psb_dprec_type) :: mld_dprec_type + type(mld_donelev_type), allocatable :: precv(:) + contains + procedure, pass(prec) :: d_apply2v => mld_d_apply2v + procedure, pass(prec) :: d_apply1v => mld_d_apply1v + end type mld_dprec_type + + ! + ! Interfaces to routines for checking the definition of the preconditioner, + ! for printing its description and for deallocating its data structure + ! + + interface mld_precfree + module procedure mld_dbase_precfree, mld_d_onelev_precfree, mld_dprec_free + end interface + + interface mld_nullify_baseprec + module procedure mld_nullify_dbaseprec + end interface + + interface mld_nullify_onelevprec + module procedure mld_nullify_d_onelevprec + end interface + + interface mld_precdescr + module procedure mld_dfile_prec_descr + end interface + + interface mld_sizeof + module procedure mld_dprec_sizeof, mld_dbaseprec_sizeof, mld_d_onelev_prec_sizeof + end interface + + interface mld_precaply + subroutine mld_dprecaply(prec,x,y,desc_data,info,trans,work) + use psb_base_mod, only : psb_d_sparse_mat, psb_desc_type, psb_dpk_ + import mld_dprec_type + type(psb_desc_type),intent(in) :: desc_data + type(mld_dprec_type), intent(in) :: prec + real(psb_dpk_),intent(in) :: x(:) + real(psb_dpk_),intent(inout) :: y(:) + integer, intent(out) :: info + character(len=1), optional :: trans + real(psb_dpk_),intent(inout), optional, target :: work(:) + end subroutine mld_dprecaply + subroutine mld_dprecaply1(prec,x,desc_data,info,trans) + use psb_base_mod, only : psb_d_sparse_mat, psb_desc_type, psb_dpk_ + import mld_dprec_type + type(psb_desc_type),intent(in) :: desc_data + type(mld_dprec_type), intent(in) :: prec + real(psb_dpk_),intent(inout) :: x(:) + integer, intent(out) :: info + character(len=1), optional :: trans + end subroutine mld_dprecaply1 + end interface + +contains + ! + ! Function returning the size of the mld_prec_type data structure + ! + + function mld_dprec_sizeof(prec) result(val) + implicit none + type(mld_dprec_type), intent(in) :: prec + integer(psb_long_int_k_) :: val + integer :: i + val = 0 + if (allocated(prec%precv)) then + do i=1, size(prec%precv) + val = val + mld_sizeof(prec%precv(i)) + end do + end if + end function mld_dprec_sizeof + + function mld_dbaseprec_sizeof(prec) result(val) + implicit none + type(mld_dbaseprec_type), intent(in) :: prec + integer(psb_long_int_k_) :: val + integer :: i + + val = 0 + if (allocated(prec%iprcparm)) then + val = val + psb_sizeof_int * size(prec%iprcparm) + if (prec%iprcparm(mld_prec_status_) == mld_prec_built_) then + select case(prec%iprcparm(mld_sub_solve_)) + case(mld_ilu_n_,mld_ilu_t_) + ! do nothing + case(mld_slu_) + case(mld_umf_) + case(mld_sludist_) + case default + end select + + end if + end if + if (allocated(prec%rprcparm)) val = val + psb_sizeof_dp * size(prec%rprcparm) + if (allocated(prec%d)) val = val + psb_sizeof_dp * size(prec%d) + if (allocated(prec%perm)) val = val + psb_sizeof_int * size(prec%perm) + if (allocated(prec%invperm)) val = val + psb_sizeof_int * size(prec%invperm) + val = val + psb_sizeof(prec%desc_data) + if (allocated(prec%av)) then + do i=1,size(prec%av) + val = val + psb_sizeof(prec%av(i)) + end do + end if + + + end function mld_dbaseprec_sizeof + + function mld_d_onelev_prec_sizeof(prec) result(val) + implicit none + type(mld_donelev_type), intent(in) :: prec + integer(psb_long_int_k_) :: val + integer :: i + + val = mld_sizeof(prec%prec) + if (allocated(prec%iprcparm)) & + & val = val + psb_sizeof_int * size(prec%iprcparm) +!!$ if (allocated(prec%ilaggr)) & +!!$ & val = val + psb_sizeof_int * size(prec%ilaggr) +!!$ if (allocated(prec%nlaggr)) & +!!$ & val = val + psb_sizeof_int * size(prec%nlaggr) + if (allocated(prec%rprcparm)) val = val + psb_sizeof_dp * size(prec%rprcparm) + val = val + psb_sizeof(prec%desc_ac) + val = val + psb_sizeof(prec%ac) + val = val + psb_sizeof(prec%map) + + end function mld_d_onelev_prec_sizeof + + ! + ! Subroutine: mld_file_prec_descr + ! Version: real + ! + ! This routine prints a description of the preconditioner to the standard + ! output or to a file. It must be called after the preconditioner has been + ! built by mld_precbld. + ! + ! Arguments: + ! p - type(mld_Tprec_type), input. + ! The preconditioner data structure to be printed out. + ! info - integer, output. + ! error code. + ! iout - integer, input, optional. + ! The id of the file where the preconditioner description + ! will be printed. If iout is not present, then the standard + ! output is condidered. + ! + subroutine mld_dfile_prec_descr(p,info,iout) + implicit none + ! Arguments + type(mld_dprec_type), intent(in) :: p + integer, intent(out) :: info + integer, intent(in), optional :: iout + + ! Local variables + integer :: ilev, nlev + integer :: ictxt, me, np + character(len=20), parameter :: name='mld_file_prec_descr' + integer :: iout_ + + info = 0 + if (present(iout)) then + iout_ = iout + else + iout_ = 6 + end if + if (iout_ < 0) iout_ = 6 + + if (allocated(p%precv)) then + ictxt = psb_cd_get_context(p%precv(1)%prec%desc_data) + + call psb_info(ictxt,me,np) + + ! + ! The preconditioner description is printed by processor psb_root_. + ! This agrees with the fact that all the parameters defining the + ! preconditioner have the same values on all the procs (this is + ! ensured by mld_precbld). + ! + if (me==psb_root_) then + + write(iout_,*) + write(iout_,'(a)') 'Preconditioner description' + nlev = size(p%precv) + if (nlev >= 1) then + ! + ! Print description of base preconditioner + ! + + write(iout_,*) ' ' + + if (nlev > 1) then + write(iout_,*) 'Multilevel Schwarz' + write(iout_,*) + write(iout_,*) 'Base preconditioner (smoother) details' + endif + + ilev = 1 + call mld_base_prec_descr(iout_,p%precv(ilev)%prec%iprcparm,info,& + & dprcparm=p%precv(ilev)%prec%rprcparm) + + end if + + if (nlev > 1) then + + ! + ! Print multilevel details + ! + write(iout_,*) + write(iout_,*) 'Multilevel details' + + do ilev = 2, nlev + if (.not.allocated(p%precv(ilev)%iprcparm)) then + info = 3111 + write(iout_,*) ' ',name,': error: inconsistent MLPREC part, should call MLD_PRECINIT' + return + endif + end do + + write(iout_,*) ' Number of levels: ',nlev + + ! + ! Currently, all the preconditioner parameters must have the same value at levels + ! 2,...,nlev-1, hence only the values at level 2 are printed + ! + + ilev=2 + call mld_ml_alg_descr(iout_,ilev,p%precv(ilev)%iprcparm, info,& + & dprcparm=p%precv(ilev)%rprcparm) + + ! + ! Coarse matrices are different at levels 2,...,nlev-1, hence related + ! info is printed separately + ! + write(iout_,*) + do ilev = 2, nlev-1 + call mld_ml_level_descr(iout_,ilev,p%precv(ilev)%iprcparm,& + & p%precv(ilev)%map%naggr,info,& + & dprcparm=p%precv(ilev)%rprcparm) + end do + + ! + ! Print coarsest level details + ! + + ilev = nlev + write(iout_,*) + call mld_ml_coarse_descr(iout_,ilev,& + & p%precv(ilev)%iprcparm,p%precv(ilev)%prec%iprcparm,& + & p%precv(ilev)%map%naggr,info,& + & dprcparm=p%precv(ilev)%rprcparm,& + & dprcparm2=p%precv(ilev)%prec%rprcparm) + end if + + endif + write(iout_,*) + else + write(iout_,*) trim(name), & + & ': Error: no base preconditioner available, something is wrong!' + info = -2 + return + endif + + + end subroutine mld_dfile_prec_descr + + ! + ! Subroutines: mld_Tbase_precfree, mld_T_onelev_precfree, mld_Tprec_free + ! Version: real/complex + ! + ! These routines deallocate the mld_Tbaseprec_type, mld_Tonelev_type and + ! mld_Tprec_type data structures. + ! + ! Arguments: + ! p - type(mld_Tbaseprec_type/mld_Tonelev_type/mld_Tprec_type), input. + ! The data structure to be deallocated. + ! info - integer, output. + ! error code. + ! + subroutine mld_dbase_precfree(p,info) + implicit none + + type(mld_dbaseprec_type), intent(inout) :: p + integer, intent(out) :: info + integer :: i + + info = 0 + + ! Actually we might just deallocate the top level array, except + ! for the inner UMFPACK or SLU stuff + + if (allocated(p%d)) then + deallocate(p%d,stat=info) + end if + + if (allocated(p%av)) then + do i=1,size(p%av) + call psb_sp_free(p%av(i),info) + if (info /= 0) then + ! Actually, we don't care here about this. + ! Just let it go. + ! return + end if + enddo + deallocate(p%av,stat=info) + end if + + if (allocated(p%desc_data%matrix_data)) & + & call psb_cdfree(p%desc_data,info) + + if (allocated(p%rprcparm)) then + deallocate(p%rprcparm,stat=info) + end if + + if (allocated(p%perm)) then + deallocate(p%perm,stat=info) + endif + + if (allocated(p%invperm)) then + deallocate(p%invperm,stat=info) + endif + + if (allocated(p%iprcparm)) then + if (p%iprcparm(mld_prec_status_) == mld_prec_built_) then + if (p%iprcparm(mld_sub_solve_)==mld_slu_) then + call mld_dslu_free(p%iprcparm(mld_slu_ptr_),info) + end if + if (p%iprcparm(mld_sub_solve_)==mld_sludist_) then + call mld_dsludist_free(p%iprcparm(mld_slud_ptr_),info) + end if + if (p%iprcparm(mld_sub_solve_)==mld_umf_) then + call mld_dumf_free(p%iprcparm(mld_umf_symptr_),& + & p%iprcparm(mld_umf_numptr_),info) + end if + end if + deallocate(p%iprcparm,stat=info) + end if + call mld_nullify_baseprec(p) + + end subroutine mld_dbase_precfree + + subroutine mld_d_onelev_precfree(p,info) + implicit none + + type(mld_donelev_type), intent(inout) :: p + integer, intent(out) :: info + integer :: i + + info = 0 + + ! Actually we might just deallocate the top level array, except + ! for the inner UMFPACK or SLU stuff + call mld_precfree(p%prec,info) + + call psb_sp_free(p%ac,info) + if (allocated(p%desc_ac%matrix_data)) & + & call psb_cdfree(p%desc_ac,info) + + if (allocated(p%rprcparm)) then + deallocate(p%rprcparm,stat=info) + end if + ! This is a pointer to something else, must not free it here. + nullify(p%base_a) + ! This is a pointer to something else, must not free it here. + nullify(p%base_desc) + + ! + ! free explicitly map??? + ! For now thanks to allocatable semantics + ! works anyway. + ! + + call mld_nullify_onelevprec(p) + end subroutine mld_d_onelev_precfree + + subroutine mld_nullify_dbaseprec(p) + implicit none + + type(mld_dbaseprec_type), intent(inout) :: p + + + end subroutine mld_nullify_dbaseprec + + subroutine mld_nullify_d_onelevprec(p) + implicit none + + type(mld_donelev_type), intent(inout) :: p + + nullify(p%base_a) + nullify(p%base_desc) + + end subroutine mld_nullify_d_onelevprec + + subroutine mld_dprec_free(p,info) + + use psb_base_mod + + implicit none + + ! Arguments + type(mld_dprec_type), intent(inout) :: p + integer, intent(out) :: info + + ! Local variables + integer :: me,err_act,i + character(len=20) :: name + + if(psb_get_errstatus().ne.0) return + info=0 + name = 'mld_dprecfree' + call psb_erractionsave(err_act) + + me=-1 + + if (allocated(p%precv)) then + do i=1,size(p%precv) + call mld_precfree(p%precv(i),info) + end do + deallocate(p%precv) + end if + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act.eq.psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine mld_dprec_free + + + + + subroutine d_base_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,work,info) + use psb_base_mod + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_base_smoother_type), intent(in) :: sm + real(psb_dpk_),intent(in) :: x(:) + real(psb_dpk_),intent(inout) :: y(:) + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='d_base_smoother_apply' + + call psb_erractionsave(err_act) + + info = 700 + call psb_errpush(info,name) + goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_base_smoother_apply + + subroutine d_base_smoother_bld(a,desc_a,sm,upd,info) + + use psb_base_mod + + Implicit None + + ! Arguments + type(psb_d_sparse_mat), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_d_base_smoother_type), intent(in) :: sm + + character, intent(in) :: upd + integer, intent(out) :: info + Integer :: err_act + character(len=20) :: name='d_base_smoother_bld' + + call psb_erractionsave(err_act) + + info = 700 + call psb_errpush(info,name) + goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine d_base_smoother_bld + + subroutine mld_d_apply2v(prec,x,y,desc_data,info,trans,work) + use psb_base_mod + type(psb_desc_type),intent(in) :: desc_data + class(mld_dprec_type), intent(in) :: prec + real(psb_dpk_),intent(in) :: x(:) + real(psb_dpk_),intent(inout) :: y(:) + integer, intent(out) :: info + character(len=1), optional :: trans + real(psb_dpk_),intent(inout), optional, target :: work(:) + Integer :: err_act + character(len=20) :: name='d_prec_apply' + + call psb_erractionsave(err_act) + + select type(prec) + type is (mld_dprec_type) + call mld_precaply(prec,x,y,desc_data,info,trans,work) + class default + info = 700 + call psb_errpush(info,name) + goto 9999 + end select + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine mld_d_apply2v + + subroutine mld_d_apply1v(prec,x,desc_data,info,trans) + use psb_base_mod + type(psb_desc_type),intent(in) :: desc_data + class(mld_dprec_type), intent(in) :: prec + real(psb_dpk_),intent(inout) :: x(:) + integer, intent(out) :: info + character(len=1), optional :: trans + Integer :: err_act + character(len=20) :: name='d_prec_apply' + + call psb_erractionsave(err_act) + + select type(prec) + type is (mld_dprec_type) + call mld_precaply(prec,x,desc_data,info,trans) + class default + info = 700 + call psb_errpush(info,name) + goto 9999 + end select + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine mld_d_apply1v + +end module mld_d_prec_type diff --git a/mlprec/mld_d_prec_type.mod b/mlprec/mld_d_prec_type.mod new file mode 100644 index 00000000..75ebcda4 --- /dev/null +++ b/mlprec/mld_d_prec_type.mod @@ -0,0 +1,152 @@ +VE 0 0 3 0 0 0 +MODULE MLD_D_PREC_TYPE,0 0 +FILE 0,mld_d_prec_type.f90 +USE MLD_BASE_PREC_TYPE 2 +REF PSB_D_BASE_PREC_TYPE(PSB_D_BASE_PREC_TYPE@PSB_D_PREC_TYPE),1 +REF PSB_DPREC_TYPE(PSB_DPREC_TYPE@PSB_D_PREC_TYPE),1 +TYPE MLD_D_BASE_SOLVER_TYPE,1,192,0 +ENDTYPE +TYPE MLD_D_BASE_SMOOTHER_TYPE,1,2272,12 + SV(0): 7,MLD_D_BASE_SOLVER_TYPE,b,0,0,400,1,1000000,4000010,0 +CONTAINS +TBP BUILD,A1,PASS(SM):D_BASE_SMOOTHER_BLD +TBP APPLY,A1,PASS(SM):D_BASE_SMOOTHER_APPLY +ENDTYPE +TYPE MLD_DBASEPREC_TYPE,1,2160,496 +(PSB_D_BASE_PREC_TYPE) + SM(0): 7,MLD_D_BASE_SMOOTHER_TYPE,b,0,0,400,1,1000000,4000010,0 + AV(12): 7,PSB_D_SPARSE_MAT,b,0,1,400,1 (1,8,0: 1,0,5,3),1000000,0,0 + D(32): 2,2,5,0,1,400,1 (1,8,0: 1,0,5,3),1000000,0,0 + DESC_DATA(52): 7,PSB_DESC_TYPE,b,0,0,0,1,0,0,0 + IPRCPARM(416): 1,3,3,0,1,400,1 (1,8,0: 1,0,5,3),1000000,0,0 + RPRCPARM(436): 2,2,5,0,1,400,1 (1,8,0: 1,0,5,3),1000000,0,0 + PERM(456): 1,3,3,0,1,400,1 (1,8,0: 1,0,5,3),1000000,0,0 + INVPERM(476): 1,3,3,0,1,400,1 (1,8,0: 1,0,5,3),1000000,0,0 +CONTAINS +ENDTYPE +TYPE MLD_DONELEV_TYPE,1,112,1740 + PREC(0): 7,MLD_DBASEPREC_TYPE,b,0,0,0,1,0,0,0 + IPRCPARM(496): 1,3,3,0,1,400,1 (1,8,0: 1,0,5,3),1000000,0,0 + RPRCPARM(516): 2,2,5,0,1,400,1 (1,8,0: 1,0,5,3),1000000,0,0 + AC(536): 7,PSB_D_SPARSE_MAT,b,0,0,0,1,0,0,0 + DESC_AC(548): 7,PSB_DESC_TYPE,b,0,0,0,1,0,0,0 + BASE_A(912): 7,PSB_D_SPARSE_MAT,b,0,0,20,1,0,0,0 += NULL + BASE_DESC(916): 7,PSB_DESC_TYPE,b,0,0,20,1,0,0,0 += NULL + MAP(920): 7,PSB_DLINMAP_TYPE,b,0,0,0,1,0,0,0 +ENDTYPE +TYPE MLD_DPREC_TYPE,1,2272,32 +(PSB_DPREC_TYPE) + PRECV(12): 7,MLD_DONELEV_TYPE,b,0,1,400,1 (1,8,0: 1,0,5,3),1000000,0,0 +CONTAINS +TBP D_APPLY2V,A1,PASS(PREC):MLD_D_APPLY2V +TBP D_APPLY1V,A1,PASS(PREC):MLD_D_APPLY1V +ENDTYPE +GMODPROC MLD_NULLIFY_ONELEVPREC: MLD_NULLIFY_D_ONELEVPREC +GMODPROC MLD_PRECFREE: MLD_DPREC_FREE +GMODPROC MLD_PRECFREE: MLD_D_ONELEV_PRECFREE +GMODPROC MLD_PRECFREE: MLD_DBASE_PRECFREE +GMODPROC MLD_SIZEOF: MLD_D_ONELEV_PREC_SIZEOF +GMODPROC MLD_SIZEOF: MLD_DBASEPREC_SIZEOF +GMODPROC MLD_SIZEOF: MLD_DPREC_SIZEOF +GMODPROC MLD_NULLIFY_BASEPREC: MLD_NULLIFY_DBASEPREC +GMODPROC MLD_PRECDESCR: MLD_DFILE_PREC_DESCR +PROC MLD_D_APPLY2V,7,8,0,17: 8,0,0,0,0,40000,1,0,80000,0 +USE PSB_BASE_MOD 2 +VAR PREC,3,,: 7,MLD_DPREC_TYPE,b,0,0,3,0,0,4000010,1 +VAR X,3,,: 2,2,5,0,1,3,0 (1,5,0: 1,0,5,3),0,1000,1 +VAR Y,3,,: 2,2,5,0,1,3,0 (1,5,0: 1,0,5,3),0,1000,3 +VAR DESC_DATA,3,,: 7,PSB_DESC_TYPE,b,0,0,3,0,0,1000,1 +VAR INFO,3,,: 1,3,3,0,0,83,0,0,1000,2 +VAR TRANS,3,,: 3,1,a,1,0,13,0,0,1000,0 +VAR WORK,3,,: 2,2,5,0,1,1013,0 (1,5,0: 1,0,5,3),0,1000,3 +ENDPROC +PROC D_BASE_SMOOTHER_BLD,5,8,0,17: 8,0,0,0,0,40000,1,0,80000,0 +USE PSB_BASE_MOD 2 +VAR A,3,,: 7,PSB_D_SPARSE_MAT,b,0,0,1003,0,0,0,1 +VAR DESC_A,3,,: 7,PSB_DESC_TYPE,b,0,0,3,0,0,0,1 +VAR SM,3,,: 7,MLD_D_BASE_SMOOTHER_TYPE,b,0,0,3,0,0,4000010,1 +VAR UPD,3,,: 3,1,a,1,0,3,0,0,0,1 +VAR INFO,3,,: 1,3,3,0,0,83,0,0,1000,2 +ENDPROC +PROC MLD_DPRECAPLY1,5,10,0,17: 8,0,0,0,0,40000,1,0,40000,0 +IMPORT MLD_DPREC_TYPE +USE PSB_BASE_MOD 2,ONLY:PSB_DPK_=>PSB_DPK_,PSB_DESC_TYPE=>PSB_DESC_TYPE,PSB_D_SPARSE_MAT=>PSB_D_SPARSE_MAT +VAR PREC,3,,: 7,MLD_DPREC_TYPE,b,0,0,3,0,0,0,1 +VAR X,3,,: 2,2,5,0,1,3,0 (1,5,0: 1,0,5,3),0,0,3 +VAR DESC_DATA,3,,: 7,PSB_DESC_TYPE,b,0,0,3,0,0,0,1 +VAR INFO,3,,: 1,3,3,0,0,3,0,0,0,2 +VAR TRANS,3,,: 3,1,a,1,0,13,0,0,0,0 +ENDPROC +PROC MLD_DFILE_PREC_DESCR,3,8,0,17: 8,0,0,0,0,40000,1,0,40000,0 +VAR P,3,,: 7,MLD_DPREC_TYPE,b,0,0,103,0,0,0,1 +VAR INFO,3,,: 1,3,3,0,0,83,0,0,1000,2 +VAR IOUT,3,,: 1,3,3,0,0,113,0,0,0,1 +ENDPROC +PROC MLD_D_APPLY1V,5,8,0,17: 8,0,0,0,0,40000,1,0,80000,0 +USE PSB_BASE_MOD 2 +VAR PREC,3,,: 7,MLD_DPREC_TYPE,b,0,0,3,0,0,4000010,1 +VAR X,3,,: 2,2,5,0,1,3,0 (1,5,0: 1,0,5,3),0,1000,3 +VAR DESC_DATA,3,,: 7,PSB_DESC_TYPE,b,0,0,3,0,0,1000,1 +VAR INFO,3,,: 1,3,3,0,0,83,0,0,1000,2 +VAR TRANS,3,,: 3,1,a,1,0,13,0,0,1000,0 +ENDPROC +PROC MLD_D_ONELEV_PREC_SIZEOF,1,8,0,17: 1,4,e,0,0,40181,1,0,40000,0 +RESULTVAR VAL,4,,MLD_D_ONELEV_PREC_SIZEOF: 1,4,e,0,0,181,0,0,0,0 +VAR PREC,3,,: 7,MLD_DONELEV_TYPE,b,0,0,103,0,0,0,1 +ENDPROC +PROC MLD_DPRECAPLY,7,10,0,17: 8,0,0,0,0,40000,1,0,40000,0 +IMPORT MLD_DPREC_TYPE +USE PSB_BASE_MOD 2,ONLY:PSB_DPK_=>PSB_DPK_,PSB_DESC_TYPE=>PSB_DESC_TYPE,PSB_D_SPARSE_MAT=>PSB_D_SPARSE_MAT +VAR PREC,3,,: 7,MLD_DPREC_TYPE,b,0,0,3,0,0,0,1 +VAR X,3,,: 2,2,5,0,1,3,0 (1,5,0: 1,0,5,3),0,0,1 +VAR Y,3,,: 2,2,5,0,1,3,0 (1,5,0: 1,0,5,3),0,0,3 +VAR DESC_DATA,3,,: 7,PSB_DESC_TYPE,b,0,0,3,0,0,0,1 +VAR INFO,3,,: 1,3,3,0,0,3,0,0,0,2 +VAR TRANS,3,,: 3,1,a,1,0,13,0,0,0,0 +VAR WORK,3,,: 2,2,5,0,1,1013,0 (1,5,0: 1,0,5,3),0,0,3 +ENDPROC +PROC MLD_DPREC_FREE,2,8,0,17: 8,0,0,0,0,40000,1,0,40000,0 +USE PSB_BASE_MOD 2 +VAR P,3,,: 7,MLD_DPREC_TYPE,b,0,0,183,0,0,0,3 +VAR INFO,3,,: 1,3,3,0,0,83,0,0,1000,2 +ENDPROC +PROC D_BASE_SMOOTHER_APPLY,9,8,0,17: 8,0,0,0,0,40000,1,0,80000,0 +USE PSB_BASE_MOD 2 +VAR ALPHA,3,,: 2,2,5,0,0,3,0,0,0,1 +VAR SM,3,,: 7,MLD_D_BASE_SMOOTHER_TYPE,b,0,0,3,0,0,4000010,1 +VAR X,3,,: 2,2,5,0,1,3,0 (1,5,0: 1,0,5,3),0,0,1 +VAR BETA,3,,: 2,2,5,0,0,3,0,0,0,1 +VAR Y,3,,: 2,2,5,0,1,3,0 (1,5,0: 1,0,5,3),0,0,3 +VAR DESC_DATA,3,,: 7,PSB_DESC_TYPE,b,0,0,3,0,0,0,1 +VAR TRANS,3,,: 3,1,a,1,0,3,0,0,0,1 +VAR WORK,3,,: 2,2,5,0,1,1003,0 (1,5,0: 1,0,5,3),0,0,3 +VAR INFO,3,,: 1,3,3,0,0,83,0,0,1000,2 +ENDPROC +PROC MLD_DBASEPREC_SIZEOF,1,8,0,17: 1,4,e,0,0,40181,1,0,40000,0 +RESULTVAR VAL,4,,MLD_DBASEPREC_SIZEOF: 1,4,e,0,0,181,0,0,0,0 +VAR PREC,3,,: 7,MLD_DBASEPREC_TYPE,b,0,0,103,0,0,0,1 +ENDPROC +PROC MLD_NULLIFY_D_ONELEVPREC,1,8,0,17: 8,0,0,0,0,40000,1,0,40000,0 +VAR P,3,,: 7,MLD_DONELEV_TYPE,b,0,0,83,0,0,0,3 +ENDPROC +PROC MLD_DBASE_PRECFREE,2,8,0,17: 8,0,0,0,0,40000,1,0,40000,0 +VAR P,3,,: 7,MLD_DBASEPREC_TYPE,b,0,0,183,0,0,1000,3 +VAR INFO,3,,: 1,3,3,0,0,183,0,0,1000,2 +ENDPROC +PROC MLD_DPREC_SIZEOF,1,8,0,17: 1,4,e,0,0,40181,1,0,40000,0 +RESULTVAR VAL,4,,MLD_DPREC_SIZEOF: 1,4,e,0,0,181,0,0,0,0 +VAR PREC,3,,: 7,MLD_DPREC_TYPE,b,0,0,103,0,0,0,1 +ENDPROC +PROC NULL,0,6,161,17: 0,0,0,118,0,20000,1,0,0,0 +ENDPROC +PROC MLD_NULLIFY_DBASEPREC,1,8,0,17: 8,0,0,0,0,40000,1,0,40000,0 +VAR P,3,,: 7,MLD_DBASEPREC_TYPE,b,0,0,3,0,0,0,3 +ENDPROC +PROC MLD_D_ONELEV_PRECFREE,2,8,0,17: 8,0,0,0,0,40000,1,0,40000,0 +VAR P,3,,: 7,MLD_DONELEV_TYPE,b,0,0,183,0,0,1000,3 +VAR INFO,3,,: 1,3,3,0,0,83,0,0,1000,2 +ENDPROC +GENERIC MLD_PRECAPLY: MLD_DPRECAPLY1,MLD_DPRECAPLY +END diff --git a/mlprec/mld_prec_type.f90 b/mlprec/mld_prec_type.f90 index 8e0f187d..d58442cc 100644 --- a/mlprec/mld_prec_type.f90 +++ b/mlprec/mld_prec_type.f90 @@ -59,2795 +59,10 @@ module mld_prec_type - ! - ! This reduces the size of .mod file. Without the ONLY clause compilation - ! blows up on some systems. - ! - use psb_base_mod, only :& - & psb_d_sparse_mat, psb_z_sparse_mat,& - & psb_s_sparse_mat, psb_c_sparse_mat,& - & psb_desc_type,& - & psb_slinmap_type, psb_dlinmap_type,& - & psb_clinmap_type, psb_zlinmap_type, & - & psb_dpk_, psb_spk_, psb_long_int_k_, & - & psb_spfree, psb_cdfree, psb_halo_, psb_none_, psb_sum_, psb_avg_, & - & psb_nohalo_, psb_square_root_, psb_toupper, psb_root_,& - & psb_sizeof_int, psb_sizeof_long_int, psb_sizeof_sp, psb_sizeof_dp, psb_sizeof,& - & psb_cd_get_context, psb_info - use psb_prec_mod, only: psb_sprec_type, psb_dprec_type,& - & psb_cprec_type, psb_zprec_type, psb_d_base_prec_type - - ! - ! Type: mld_Tprec_type. - ! - ! It is the data type containing all the information about the multilevel - ! preconditioner (here and in the following 'T' denotes 'd', 's', 'c' and - ! 'z', according to the real/complex, single/double precision version of - ! MLD2P4). It consists of an array of 'one-level' intermediate data structures - ! of type mld_Tonelev_type, each containing the information needed to apply - ! the smoothing and the coarse-space correction at a generic level. - ! - ! type mld_Tprec_type - ! type(mld_Tonelev_type), allocatable :: precv(:) - ! end type mld_Tprec_type - ! - ! Note that the levels are numbered in increasing order starting from - ! the finest one and the number of levels is given by size(precv(:)). - ! - ! - ! Type: mld_Tonelev_type. - ! - ! It is the data type containing the necessary items for the current - ! level (essentially, the base preconditioner, the current-level matrix - ! and the restriction and prolongation operators). - ! - ! type mld_Tonelev_type - ! type(mld_Tbaseprec_type) :: prec - ! integer, allocatable :: iprcparm(:) - ! real(psb_Tpk_), allocatable :: rprcparm(:) - ! type(psb_T_sparse_mat) :: ac - ! type(psb_desc_type) :: desc_ac - ! type(psb_T_sparse_mat), pointer :: base_a => null() - ! type(psb_desc_type), pointer :: base_desc => null() - ! type(psb_Tlinmap_type) :: map - ! end type mld_Tonelev_type - ! - ! Note that psb_Tpk denotes the kind of the real data type to be chosen - ! according to single/double precision version of MLD2P4. - ! - ! prec - type(mld_Tbaseprec_type). - ! The current level preconditioner (aka smoother). - ! iprcparm - integer, dimension(:), allocatable. - ! The integer parameters defining the multilevel strategy. - ! rprcparm - real(psb_Ypk_), dimension(:), allocatable. - ! The real parameters defining the multilevel strategy. - ! ac - The local part of the current-level matrix, built by - ! coarsening the previous-level matrix. - ! desc_ac - type(psb_desc_type). - ! The communication descriptor associated to the matrix - ! stored in ac. - ! base_a - type(psb_z_sparse_mat), pointer. - ! Pointer (really a pointer!) to the local part of the current - ! matrix (so we have a unified treatment of residuals). - ! We need this to avoid passing explicitly the current matrix - ! to the routine which applies the preconditioner. - ! base_desc - type(psb_desc_type), pointer. - ! Pointer to the communication descriptor associated to the - ! matrix pointed by base_a. - ! map - Stores the maps (restriction and prolongation) between the - ! vector spaces associated to the index spaces of the previous - ! and current levels. - ! - ! - ! Type: mld_Tbaseprec_type. - ! - ! It holds the smoother (base preconditioner) at a single level. - ! - ! type mld_Tbaseprec_type - ! type(psb_T_sparse_mat), allocatable :: av(:) - ! IntrType(psb_Tpk_), allocatable :: d(:) - ! type(psb_desc_type) :: desc_data - ! integer, allocatable :: iprcparm(:) - ! real(psb_Tpk_), allocatable :: rprcparm(:) - ! integer, allocatable :: perm(:), invperm(:) - ! end type mld_sbaseprec_type - ! - ! Note that IntrType denotes the real or complex data type, and psb_Tpk denotes - ! the kind of the real or complex type, according to the real/complex, single/double - ! precision version of MLD2P4. - ! - ! av - type(psb_T_sparse_mat), dimension(:), allocatable(:). - ! The sparse matrices needed to apply the preconditioner at - ! the current level ilev. - ! av(mld_l_pr_) - The L factor of the ILU factorization of the local - ! diagonal block of the current-level matrix A(ilev). - ! av(mld_u_pr_) - The U factor of the ILU factorization of the local - ! diagonal block of A(ilev), except its diagonal entries - ! (stored in d). - ! av(mld_ap_nd_) - The entries of the local part of A(ilev) outside - ! the diagonal block, for block-Jacobi sweeps. - ! d - real/complex(psb_Tpk_), dimension(:), allocatable. - ! The diagonal entries of the U factor in the ILU factorization - ! of A(ilev). - ! desc_data - type(psb_desc_type). - ! The communication descriptor associated to the base preconditioner, - ! i.e. to the sparse matrices needed to apply the base preconditioner - ! at the current level. - ! iprcparm - integer, dimension(:), allocatable. - ! The integer parameters defining the base preconditioner K(ilev) - ! (the iprcparm entries and values are specified below). - ! rprcparm - real(psb_Tpk_), dimension(:), allocatable. - ! The real parameters defining the base preconditioner K(ilev) - ! (the rprcparm entries and values are specified below). - ! perm - integer, dimension(:), allocatable. - ! The row and column permutations applied to the local part of - ! A(ilev) (defined only if iprcparm(mld_sub_ren_)>0). - ! invperm - integer, dimension(:), allocatable. - ! The inverse of the permutation stored in perm. - ! - ! Note that when the LU factorization of the (local part of the) matrix A(ilev) is - ! computed instead of the ILU one, by using UMFPACK, SuperLU or SuperLU_dist, the - ! corresponding L and U factors are stored in data structures provided by those - ! packages and pointed by prec%iprcparm(mld_umf_ptr), prec%iprcparm(mld_slu_ptr) - ! or prec%iprcparm(mld_slud_ptr). - ! - - type mld_sbaseprec_type - type(psb_s_sparse_mat), allocatable :: av(:) - real(psb_spk_), allocatable :: d(:) - type(psb_desc_type) :: desc_data - integer, allocatable :: iprcparm(:) - real(psb_spk_), allocatable :: rprcparm(:) - integer, allocatable :: perm(:), invperm(:) - end type mld_sbaseprec_type - - type mld_sonelev_type - type(mld_sbaseprec_type) :: prec - integer, allocatable :: iprcparm(:) - real(psb_spk_), allocatable :: rprcparm(:) - type(psb_s_sparse_mat) :: ac - type(psb_desc_type) :: desc_ac - type(psb_s_sparse_mat), pointer :: base_a => null() - type(psb_desc_type), pointer :: base_desc => null() - type(psb_slinmap_type) :: map - end type mld_sonelev_type - - type, extends(psb_sprec_type) :: mld_sprec_type - type(mld_sonelev_type), allocatable :: precv(:) - contains - procedure, pass(prec) :: s_apply2v => mld_s_apply2v - procedure, pass(prec) :: s_apply1v => mld_s_apply1v - end type mld_sprec_type - - type, extends(psb_d_base_prec_type) :: mld_dbaseprec_type - type(psb_d_sparse_mat), allocatable :: av(:) - real(psb_dpk_), allocatable :: d(:) - type(psb_desc_type) :: desc_data - integer, allocatable :: iprcparm(:) - real(psb_dpk_), allocatable :: rprcparm(:) - integer, allocatable :: perm(:), invperm(:) - end type mld_dbaseprec_type - - type mld_donelev_type - class(psb_d_base_prec_type), allocatable :: bprec - type(mld_dbaseprec_type) :: prec - integer, allocatable :: iprcparm(:) - real(psb_dpk_), allocatable :: rprcparm(:) - type(psb_d_sparse_mat) :: ac - type(psb_desc_type) :: desc_ac - type(psb_d_sparse_mat), pointer :: base_a => null() - type(psb_desc_type), pointer :: base_desc => null() - type(psb_dlinmap_type) :: map - end type mld_donelev_type - - type, extends(psb_dprec_type) :: mld_dprec_type - type(mld_donelev_type), allocatable :: precv(:) - contains - procedure, pass(prec) :: d_apply2v => mld_d_apply2v - procedure, pass(prec) :: d_apply1v => mld_d_apply1v - end type mld_dprec_type - - - type mld_cbaseprec_type - type(psb_c_sparse_mat), allocatable :: av(:) - complex(psb_spk_), allocatable :: d(:) - type(psb_desc_type) :: desc_data - integer, allocatable :: iprcparm(:) - real(psb_spk_), allocatable :: rprcparm(:) - integer, allocatable :: perm(:), invperm(:) - end type mld_cbaseprec_type - - type mld_conelev_type - type(mld_cbaseprec_type) :: prec - integer, allocatable :: iprcparm(:) - real(psb_spk_), allocatable :: rprcparm(:) - type(psb_c_sparse_mat) :: ac - type(psb_desc_type) :: desc_ac - type(psb_c_sparse_mat), pointer :: base_a => null() - type(psb_desc_type), pointer :: base_desc => null() - type(psb_clinmap_type) :: map - end type mld_conelev_type - - type, extends(psb_cprec_type) :: mld_cprec_type - type(mld_conelev_type), allocatable :: precv(:) - contains - procedure, pass(prec) :: c_apply2v => mld_c_apply2v - procedure, pass(prec) :: c_apply1v => mld_c_apply1v - end type mld_cprec_type - - type mld_zbaseprec_type - type(psb_z_sparse_mat), allocatable :: av(:) - complex(psb_dpk_), allocatable :: d(:) - type(psb_desc_type) :: desc_data - integer, allocatable :: iprcparm(:) - real(psb_dpk_), allocatable :: rprcparm(:) - integer, allocatable :: perm(:), invperm(:) - end type mld_zbaseprec_type - - type mld_zonelev_type - type(mld_zbaseprec_type) :: prec - integer, allocatable :: iprcparm(:) - real(psb_dpk_), allocatable :: rprcparm(:) - type(psb_z_sparse_mat) :: ac - type(psb_desc_type) :: desc_ac - type(psb_z_sparse_mat), pointer :: base_a => null() - type(psb_desc_type), pointer :: base_desc => null() - type(psb_zlinmap_type) :: map - end type mld_zonelev_type - - type, extends(psb_zprec_type) :: mld_zprec_type - type(mld_zonelev_type), allocatable :: precv(:) - contains - procedure, pass(prec) :: z_apply2v => mld_z_apply2v - procedure, pass(prec) :: z_apply1v => mld_z_apply1v - end type mld_zprec_type - - - ! - ! Entries in iprcparm - ! - ! These are in baseprec - ! - integer, parameter :: mld_smoother_type_ = 1 - integer, parameter :: mld_sub_solve_ = 2 - integer, parameter :: mld_sub_restr_ = 3 - integer, parameter :: mld_sub_prol_ = 4 - integer, parameter :: mld_sub_ren_ = 5 - integer, parameter :: mld_sub_ovr_ = 6 - integer, parameter :: mld_sub_fillin_ = 8 - integer, parameter :: mld_smoother_sweeps_ = 9 - !! 2 ints for 64 bit versions - integer, parameter :: mld_slu_ptr_ = 10 - integer, parameter :: mld_umf_symptr_ = 12 - integer, parameter :: mld_umf_numptr_ = 14 - integer, parameter :: mld_slud_ptr_ = 16 - integer, parameter :: mld_prec_status_ = 18 - ! - ! These are in onelev - ! - integer, parameter :: mld_ml_type_ = 20 - integer, parameter :: mld_smoother_pos_ = 21 - integer, parameter :: mld_aggr_kind_ = 22 - integer, parameter :: mld_aggr_alg_ = 23 - integer, parameter :: mld_aggr_omega_alg_ = 24 - integer, parameter :: mld_aggr_eig_ = 25 - integer, parameter :: mld_aggr_filter_ = 26 - integer, parameter :: mld_coarse_mat_ = 27 - integer, parameter :: mld_coarse_solve_ = 28 - integer, parameter :: mld_coarse_sweeps_ = 29 - integer, parameter :: mld_coarse_fillin_ = 30 - integer, parameter :: mld_coarse_subsolve_ = 31 - integer, parameter :: mld_ifpsz_ = 32 - - ! - ! Legal values for entry: mld_smoother_type_ - ! - integer, parameter :: mld_min_prec_=0, mld_noprec_=0, mld_diag_=1, mld_bjac_=2,& - & mld_as_=3, mld_max_prec_=3 - ! - ! Legal values for entry: mld_sub_solve_ - ! - integer, parameter :: mld_f_none_=0,mld_ilu_n_=1,mld_milu_n_=2, mld_ilu_t_=3 - integer, parameter :: mld_slu_=4, mld_umf_=5, mld_sludist_=6, mld_max_sub_solve_=mld_sludist_ - ! - ! Legal values for entry: mld_sub_ren_ - ! - integer, parameter :: mld_renum_none_=0, mld_renum_glb_=1, mld_renum_gps_=2 - ! For the time being we are disabling GPS renumbering. - integer, parameter :: mld_max_renum_=1 - ! - ! Legal values for entry: mld_ml_type_ - ! - integer, parameter :: mld_no_ml_=0, mld_add_ml_=1, mld_mult_ml_=2 - integer, parameter :: mld_new_ml_prec_=3, mld_max_ml_type_=mld_mult_ml_ - ! - ! Legal values for entry: mld_smoother_pos_ - ! - integer, parameter :: mld_pre_smooth_=1, mld_post_smooth_=2,& - & mld_twoside_smooth_=3, mld_max_smooth_=mld_twoside_smooth_ - ! - ! Legal values for entry: mld_aggr_kind_ - ! - integer, parameter :: mld_no_smooth_=0, mld_smooth_prol_=1 - integer, parameter :: mld_min_energy_=2, mld_biz_prol_=3 - ! Disabling biz_prol for the time being. - integer, parameter :: mld_max_aggr_kind_=mld_min_energy_ - ! - ! Legal values for entry: mld_aggr_filter_ - ! - integer, parameter :: mld_no_filter_mat_=0, mld_filter_mat_=1 - integer, parameter :: mld_max_filter_mat_=mld_no_filter_mat_ - ! - ! Legal values for entry: mld_aggr_alg_ - ! - integer, parameter :: mld_dec_aggr_=0, mld_sym_dec_aggr_=1 - integer, parameter :: mld_glb_aggr_=2, mld_new_dec_aggr_=3 - integer, parameter :: mld_new_glb_aggr_=4 - integer, parameter :: mld_max_aggr_alg_=mld_dec_aggr_ - - ! - ! Legal values for entry: mld_aggr_omega_alg_ - ! - integer, parameter :: mld_eig_est_=0, mld_user_choice_=999 - ! - ! Legal values for entry: mld_aggr_eig_ - ! - integer, parameter :: mld_max_norm_=0 - ! - ! Legal values for entry: mld_coarse_mat_ - ! - integer, parameter :: mld_distr_mat_=0, mld_repl_mat_=1 - integer, parameter :: mld_max_coarse_mat_=mld_repl_mat_ - ! - ! Legal values for entry: mld_prec_status_ - ! - integer, parameter :: mld_prec_built_=98765 - - ! - ! Entries in rprcparm: ILU(k,t) threshold, smoothed aggregation omega - ! - integer, parameter :: mld_sub_iluthrs_ = 1 - integer, parameter :: mld_aggr_omega_val_ = 2 - integer, parameter :: mld_aggr_thresh_ = 3 - integer, parameter :: mld_coarse_iluthrs_ = 4 - integer, parameter :: mld_rfpsz_ = 8 - - ! - ! Fields for sparse matrices ensembles stored in av() - ! - integer, parameter :: mld_l_pr_=1, mld_u_pr_=2, mld_bp_ilu_avsz_=2 - integer, parameter :: mld_ap_nd_=3, mld_ac_=4, mld_sm_pr_t_=5, mld_sm_pr_=6 - integer, parameter :: mld_smth_avsz_=6, mld_max_avsz_=mld_smth_avsz_ - - ! - ! Character constants used by mld_file_prec_descr - ! - character(len=19), parameter, private :: & - & eigen_estimates(0:0)=(/'infinity norm '/) - character(len=19), parameter, private :: & - & smooth_names(1:3)=(/'pre-smoothing ','post-smoothing ',& - & 'pre/post-smoothing'/) - character(len=15), parameter, private :: & - & aggr_kinds(0:3)=(/'nonsmoothed ','smoothed ',& - & 'min energy ','bizr. smoothed'/) - character(len=15), parameter, private :: & - & matrix_names(0:1)=(/'distributed ','replicated '/) - character(len=18), parameter, private :: & - & aggr_names(0:4)=(/'local aggregation ','sym. local aggr. ',& - & 'global aggregation', 'new local aggr. ','new global aggr. '/) - character(len=6), parameter, private :: & - & restrict_names(0:4)=(/'none ','halo ',' ',' ',' '/) - character(len=12), parameter, private :: & - & prolong_names(0:3)=(/'none ','sum ','average ','square root'/) - character(len=15), parameter, private :: & - & ml_names(0:3)=(/'none ','additive ','multiplicative',& - & 'new ML '/) - character(len=15), parameter, private :: & - & fact_names(0:6)=(/'none ','ILU(n) ',& - & 'MILU(n) ','ILU(t,n) ',& - & 'SuperLU ','UMFPACK LU ',& - & 'SuperLU_Dist '/) - - ! - ! Interfaces to routines for checking the definition of the preconditioner, - ! for printing its description and for deallocating its data structure - ! - - interface mld_precfree - module procedure mld_sbase_precfree, mld_cbase_precfree,& - & mld_dbase_precfree, mld_zbase_precfree, & - & mld_s_onelev_precfree, mld_d_onelev_precfree, & - & mld_c_onelev_precfree, mld_z_onelev_precfree, & - & mld_sprec_free, mld_dprec_free, & - & mld_cprec_free, mld_zprec_free - end interface - - interface mld_nullify_baseprec - module procedure mld_nullify_sbaseprec, mld_nullify_cbaseprec,& - & mld_nullify_dbaseprec, mld_nullify_zbaseprec - end interface - - interface mld_nullify_onelevprec - module procedure mld_nullify_s_onelevprec, mld_nullify_d_onelevprec,& - & mld_nullify_c_onelevprec, mld_nullify_z_onelevprec - end interface - - interface mld_check_def - module procedure mld_icheck_def, mld_scheck_def, mld_dcheck_def - end interface - - - interface mld_precdescr - module procedure mld_file_prec_descr, & - & mld_zfile_prec_descr,& - & mld_sfile_prec_descr,& - & mld_cfile_prec_descr - end interface - - interface mld_prec_short_descr - module procedure mld_prec_short_descr, mld_zprec_short_descr - end interface - - interface mld_sizeof - module procedure mld_sprec_sizeof, mld_cprec_sizeof, & - & mld_dprec_sizeof, mld_zprec_sizeof, & - & mld_sbaseprec_sizeof, mld_cbaseprec_sizeof,& - & mld_dbaseprec_sizeof, mld_zbaseprec_sizeof, & - & mld_s_onelev_prec_sizeof, mld_d_onelev_prec_sizeof,& - & mld_c_onelev_prec_sizeof, mld_z_onelev_prec_sizeof - end interface - - - - interface mld_precaply - subroutine mld_sprecaply(prec,x,y,desc_data,info,trans,work) - use psb_base_mod, only : psb_s_sparse_mat, psb_desc_type, psb_spk_ - import mld_sprec_type - type(psb_desc_type),intent(in) :: desc_data - type(mld_sprec_type), intent(in) :: prec - real(psb_spk_),intent(in) :: x(:) - real(psb_spk_),intent(inout) :: y(:) - integer, intent(out) :: info - character(len=1), optional :: trans - real(psb_spk_),intent(inout), optional, target :: work(:) - end subroutine mld_sprecaply - subroutine mld_sprecaply1(prec,x,desc_data,info,trans) - use psb_base_mod, only : psb_s_sparse_mat, psb_desc_type, psb_spk_ - import mld_sprec_type - type(psb_desc_type),intent(in) :: desc_data - type(mld_sprec_type), intent(in) :: prec - real(psb_spk_),intent(inout) :: x(:) - integer, intent(out) :: info - character(len=1), optional :: trans - end subroutine mld_sprecaply1 - subroutine mld_dprecaply(prec,x,y,desc_data,info,trans,work) - use psb_base_mod, only : psb_d_sparse_mat, psb_desc_type, psb_dpk_ - import mld_dprec_type - type(psb_desc_type),intent(in) :: desc_data - type(mld_dprec_type), intent(in) :: prec - real(psb_dpk_),intent(in) :: x(:) - real(psb_dpk_),intent(inout) :: y(:) - integer, intent(out) :: info - character(len=1), optional :: trans - real(psb_dpk_),intent(inout), optional, target :: work(:) - end subroutine mld_dprecaply - subroutine mld_dprecaply1(prec,x,desc_data,info,trans) - use psb_base_mod, only : psb_d_sparse_mat, psb_desc_type, psb_dpk_ - import mld_dprec_type - type(psb_desc_type),intent(in) :: desc_data - type(mld_dprec_type), intent(in) :: prec - real(psb_dpk_),intent(inout) :: x(:) - integer, intent(out) :: info - character(len=1), optional :: trans - end subroutine mld_dprecaply1 - subroutine mld_cprecaply(prec,x,y,desc_data,info,trans,work) - use psb_base_mod, only : psb_c_sparse_mat, psb_desc_type, psb_spk_ - import mld_cprec_type - type(psb_desc_type),intent(in) :: desc_data - type(mld_cprec_type), intent(in) :: prec - complex(psb_spk_),intent(in) :: x(:) - complex(psb_spk_),intent(inout) :: y(:) - integer, intent(out) :: info - character(len=1), optional :: trans - complex(psb_spk_),intent(inout), optional, target :: work(:) - end subroutine mld_cprecaply - subroutine mld_cprecaply1(prec,x,desc_data,info,trans) - use psb_base_mod, only : psb_c_sparse_mat, psb_desc_type, psb_spk_ - import mld_cprec_type - type(psb_desc_type),intent(in) :: desc_data - type(mld_cprec_type), intent(in) :: prec - complex(psb_spk_),intent(inout) :: x(:) - integer, intent(out) :: info - character(len=1), optional :: trans - end subroutine mld_cprecaply1 - subroutine mld_zprecaply(prec,x,y,desc_data,info,trans,work) - use psb_base_mod, only : psb_z_sparse_mat, psb_desc_type, psb_dpk_ - import mld_zprec_type - type(psb_desc_type),intent(in) :: desc_data - type(mld_zprec_type), intent(in) :: prec - complex(psb_dpk_),intent(in) :: x(:) - complex(psb_dpk_),intent(inout) :: y(:) - integer, intent(out) :: info - character(len=1), optional :: trans - complex(psb_dpk_),intent(inout), optional, target :: work(:) - end subroutine mld_zprecaply - subroutine mld_zprecaply1(prec,x,desc_data,info,trans) - use psb_base_mod, only : psb_z_sparse_mat, psb_desc_type, psb_dpk_ - import mld_zprec_type - type(psb_desc_type),intent(in) :: desc_data - type(mld_zprec_type), intent(in) :: prec - complex(psb_dpk_),intent(inout) :: x(:) - integer, intent(out) :: info - character(len=1), optional :: trans - end subroutine mld_zprecaply1 - end interface - - - - - - - -contains - - ! - ! Subroutine: mld_stringval - ! - ! This routine converts the string contained into string into the corresponding - ! integer value. - ! - ! Arguments: - ! string - character(len=*), input. - ! The string to be converted. - ! val - integer, output. - ! The integer value corresponding to the string - ! info - integer, output. - ! Error code. - ! - subroutine mld_stringval(string,val,info) - implicit none - ! Arguments - character(len=*), intent(in) :: string - integer, intent(out) :: val, info - character(len=*), parameter :: name='mld_stringval' - - info = 0 - select case(psb_toupper(trim(string))) - case('NONE') - val = 0 - case('HALO') - val = psb_halo_ - case('SUM') - val = psb_sum_ - case('AVG') - val = psb_avg_ - case('FACT_NONE') - val = mld_f_none_ - case('ILU') - val = mld_ilu_n_ - case('MILU') - val = mld_milu_n_ - case('ILUT') - val = mld_ilu_t_ - case('UMF') - val = mld_umf_ - case('SLU') - val = mld_slu_ - case('SLUDIST') - val = mld_sludist_ - case('ADD') - val = mld_add_ml_ - case('MULT') - val = mld_mult_ml_ - case('DEC') - val = mld_dec_aggr_ - case('SYMDEC') - val = mld_sym_dec_aggr_ - case('GLB') - val = mld_glb_aggr_ - case('REPL') - val = mld_repl_mat_ - case('DIST') - val = mld_distr_mat_ - case('NONSMOOTHED') - val = mld_no_smooth_ - case('SMOOTHED') - val = mld_smooth_prol_ - case('MINENERGY') - val = mld_min_energy_ - case('PRE') - val = mld_pre_smooth_ - case('POST') - val = mld_post_smooth_ - case('TWOSIDE') - val = mld_twoside_smooth_ - case('NOPREC') - val = mld_noprec_ - case('DIAG') - val = mld_diag_ - case('BJAC') - val = mld_bjac_ - case('AS') - val = mld_as_ - case('RENUM_NONE') - val = mld_renum_none_ - case('RENUM_GLOBAL') - val = mld_renum_glb_ - case('RENUM_GPS') - val = mld_renum_gps_ - case('A_NORMI') - val = mld_max_norm_ - case('USER_CHOICE') - val = mld_user_choice_ - case('EIG_EST') - val = mld_eig_est_ - case('FILTER') - val = mld_filter_mat_ - case('NO_FILTER') - val = mld_no_filter_mat_ - case default - val = -1 - info = -1 - end select - if (info /= 0) then - write(0,*) name,': Error: unknown request: "',trim(string),'"' - end if - end subroutine mld_stringval - - ! - ! Function returning the size of the mld_prec_type data structure - ! - - function mld_sprec_sizeof(prec) result(val) - implicit none - type(mld_sprec_type), intent(in) :: prec - integer(psb_long_int_k_) :: val - integer :: i - val = 0 - if (allocated(prec%precv)) then - do i=1, size(prec%precv) - val = val + mld_sizeof(prec%precv(i)) - end do - end if - end function mld_sprec_sizeof - - function mld_dprec_sizeof(prec) result(val) - implicit none - type(mld_dprec_type), intent(in) :: prec - integer(psb_long_int_k_) :: val - integer :: i - val = 0 - if (allocated(prec%precv)) then - do i=1, size(prec%precv) - val = val + mld_sizeof(prec%precv(i)) - end do - end if - end function mld_dprec_sizeof - - function mld_cprec_sizeof(prec) result(val) - implicit none - type(mld_cprec_type), intent(in) :: prec - integer(psb_long_int_k_) :: val - integer :: i - val = 0 - if (allocated(prec%precv)) then - do i=1, size(prec%precv) - val = val + mld_sizeof(prec%precv(i)) - end do - end if - end function mld_cprec_sizeof - - function mld_zprec_sizeof(prec) result(val) - implicit none - type(mld_zprec_type), intent(in) :: prec - integer(psb_long_int_k_) :: val - integer :: i - val = 0 - if (allocated(prec%precv)) then - do i=1, size(prec%precv) - val = val + mld_sizeof(prec%precv(i)) - end do - end if - end function mld_zprec_sizeof - - ! - ! Function returning the size of the mld_baseprec_type data structure - ! - - function mld_sbaseprec_sizeof(prec) result(val) - implicit none - type(mld_sbaseprec_type), intent(in) :: prec - integer(psb_long_int_k_) :: val - integer :: i - - val = 0 - if (allocated(prec%iprcparm)) then - val = val + psb_sizeof_int * size(prec%iprcparm) - if (prec%iprcparm(mld_prec_status_) == mld_prec_built_) then - select case(prec%iprcparm(mld_sub_solve_)) - case(mld_ilu_n_,mld_ilu_t_) - ! do nothing - case(mld_slu_) - case(mld_umf_) - case(mld_sludist_) - case default - end select - - end if - end if - if (allocated(prec%rprcparm)) val = val + psb_sizeof_sp * size(prec%rprcparm) - if (allocated(prec%d)) val = val + psb_sizeof_sp * size(prec%d) - if (allocated(prec%perm)) val = val + psb_sizeof_int * size(prec%perm) - if (allocated(prec%invperm)) val = val + psb_sizeof_int * size(prec%invperm) - val = val + psb_sizeof(prec%desc_data) - if (allocated(prec%av)) then - do i=1,size(prec%av) - val = val + psb_sizeof(prec%av(i)) - end do - end if - - end function mld_sbaseprec_sizeof - - function mld_dbaseprec_sizeof(prec) result(val) - implicit none - type(mld_dbaseprec_type), intent(in) :: prec - integer(psb_long_int_k_) :: val - integer :: i - - val = 0 - if (allocated(prec%iprcparm)) then - val = val + psb_sizeof_int * size(prec%iprcparm) - if (prec%iprcparm(mld_prec_status_) == mld_prec_built_) then - select case(prec%iprcparm(mld_sub_solve_)) - case(mld_ilu_n_,mld_ilu_t_) - ! do nothing - case(mld_slu_) - case(mld_umf_) - case(mld_sludist_) - case default - end select - - end if - end if - if (allocated(prec%rprcparm)) val = val + psb_sizeof_dp * size(prec%rprcparm) - if (allocated(prec%d)) val = val + psb_sizeof_dp * size(prec%d) - if (allocated(prec%perm)) val = val + psb_sizeof_int * size(prec%perm) - if (allocated(prec%invperm)) val = val + psb_sizeof_int * size(prec%invperm) - val = val + psb_sizeof(prec%desc_data) - if (allocated(prec%av)) then - do i=1,size(prec%av) - val = val + psb_sizeof(prec%av(i)) - end do - end if - - - end function mld_dbaseprec_sizeof - - function mld_cbaseprec_sizeof(prec) result(val) - implicit none - type(mld_cbaseprec_type), intent(in) :: prec - integer(psb_long_int_k_) :: val - integer :: i - - val = 0 - if (allocated(prec%iprcparm)) then - val = val + psb_sizeof_int * size(prec%iprcparm) - if (prec%iprcparm(mld_prec_status_) == mld_prec_built_) then - select case(prec%iprcparm(mld_sub_solve_)) - case(mld_ilu_n_,mld_ilu_t_) - ! do nothing - case(mld_slu_) - case(mld_umf_) - case(mld_sludist_) - case default - end select - - end if - end if - if (allocated(prec%rprcparm)) val = val + psb_sizeof_sp * size(prec%rprcparm) - if (allocated(prec%d)) val = val + 2 * psb_sizeof_sp * size(prec%d) - if (allocated(prec%perm)) val = val + psb_sizeof_int * size(prec%perm) - if (allocated(prec%invperm)) val = val + psb_sizeof_int * size(prec%invperm) - val = val + psb_sizeof(prec%desc_data) - if (allocated(prec%av)) then - do i=1,size(prec%av) - val = val + psb_sizeof(prec%av(i)) - end do - end if - - end function mld_cbaseprec_sizeof - - function mld_zbaseprec_sizeof(prec) result(val) - implicit none - type(mld_zbaseprec_type), intent(in) :: prec - integer(psb_long_int_k_) :: val - integer :: i - - val = 0 - if (allocated(prec%iprcparm)) then - val = val + psb_sizeof_int * size(prec%iprcparm) - if (prec%iprcparm(mld_prec_status_) == mld_prec_built_) then - select case(prec%iprcparm(mld_sub_solve_)) - case(mld_ilu_n_,mld_ilu_t_) - ! do nothing - case(mld_slu_) - case(mld_umf_) - case(mld_sludist_) - case default - end select - - end if - end if - if (allocated(prec%rprcparm)) val = val + psb_sizeof_dp * size(prec%rprcparm) - if (allocated(prec%d)) val = val + 2 * psb_sizeof_dp * size(prec%d) - if (allocated(prec%perm)) val = val + psb_sizeof_int * size(prec%perm) - if (allocated(prec%invperm)) val = val + psb_sizeof_int * size(prec%invperm) - val = val + psb_sizeof(prec%desc_data) - if (allocated(prec%av)) then - do i=1,size(prec%av) - val = val + psb_sizeof(prec%av(i)) - end do - end if - - end function mld_zbaseprec_sizeof - - function mld_s_onelev_prec_sizeof(prec) result(val) - implicit none - type(mld_sonelev_type), intent(in) :: prec - integer(psb_long_int_k_) :: val - integer :: i - - val = mld_sizeof(prec%prec) - if (allocated(prec%iprcparm)) then - val = val + psb_sizeof_int * size(prec%iprcparm) - end if - if (allocated(prec%rprcparm)) val = val + psb_sizeof_sp * size(prec%rprcparm) - val = val + psb_sizeof(prec%desc_ac) - val = val + psb_sizeof(prec%ac) - val = val + psb_sizeof(prec%map) - - end function mld_s_onelev_prec_sizeof - - function mld_d_onelev_prec_sizeof(prec) result(val) - implicit none - type(mld_donelev_type), intent(in) :: prec - integer(psb_long_int_k_) :: val - integer :: i - - val = mld_sizeof(prec%prec) - if (allocated(prec%iprcparm)) then - val = val + psb_sizeof_int * size(prec%iprcparm) - end if - if (allocated(prec%rprcparm)) val = val + psb_sizeof_dp * size(prec%rprcparm) - val = val + psb_sizeof(prec%desc_ac) - val = val + psb_sizeof(prec%ac) - val = val + psb_sizeof(prec%map) - - end function mld_d_onelev_prec_sizeof - - function mld_c_onelev_prec_sizeof(prec) result(val) - implicit none - type(mld_conelev_type), intent(in) :: prec - integer(psb_long_int_k_) :: val - integer :: i - - val = mld_sizeof(prec%prec) - if (allocated(prec%iprcparm)) then - val = val + psb_sizeof_int * size(prec%iprcparm) - end if - if (allocated(prec%rprcparm)) val = val + psb_sizeof_sp * size(prec%rprcparm) - val = val + psb_sizeof(prec%desc_ac) - val = val + psb_sizeof(prec%ac) - val = val + psb_sizeof(prec%map) - - end function mld_c_onelev_prec_sizeof - - function mld_z_onelev_prec_sizeof(prec) result(val) - implicit none - type(mld_zonelev_type), intent(in) :: prec - integer(psb_long_int_k_) :: val - integer :: i - - val = mld_sizeof(prec%prec) - if (allocated(prec%iprcparm)) then - val = val + psb_sizeof_int * size(prec%iprcparm) - end if - if (allocated(prec%rprcparm)) val = val + psb_sizeof_dp * size(prec%rprcparm) - val = val + psb_sizeof(prec%desc_ac) - val = val + psb_sizeof(prec%ac) - val = val + psb_sizeof(prec%map) - - end function mld_z_onelev_prec_sizeof - - ! - ! Routines printing out a description of the preconditioner - ! - - subroutine mld_base_prec_descr(iout,iprcparm, info,rprcparm,dprcparm) - implicit none - integer, intent(in) :: iprcparm(:),iout - integer, intent(out) :: info - real(psb_spk_), intent(in), optional :: rprcparm(:) - real(psb_dpk_), intent(in), optional :: dprcparm(:) - - info = 0 - if (count((/ present(rprcparm),present(dprcparm) /)) /= 1) then - info=581 -!!$ call psb_errpush(info,name,a_err=" rprcparm, dprcparm") - return - endif - - select case(iprcparm(mld_smoother_type_)) - case(mld_noprec_) - write(iout,*) ' No preconditioning' - case(mld_diag_) - write(iout,*) ' Diagonal scaling' - case(mld_bjac_) - write(iout,*) ' Block Jacobi with ',& - & fact_names(iprcparm(mld_sub_solve_)) - select case(iprcparm(mld_sub_solve_)) - case(mld_ilu_n_,mld_milu_n_) - write(iout,*) ' Fill level:',iprcparm(mld_sub_fillin_) - case(mld_ilu_t_) - write(iout,*) ' Fill level:',iprcparm(mld_sub_fillin_) - if (present(rprcparm)) then - write(iout,*) ' Fill threshold :',rprcparm(mld_sub_iluthrs_) - else - write(iout,*) ' Fill threshold :',dprcparm(mld_sub_iluthrs_) - end if - case(mld_slu_,mld_umf_,mld_sludist_) - case default - write(iout,*) ' Should never get here!' - end select - case(mld_as_) - write(iout,*) ' Additive Schwarz with ',& - & fact_names(iprcparm(mld_sub_solve_)) - select case(iprcparm(mld_sub_solve_)) - case(mld_ilu_n_,mld_milu_n_) - write(iout,*) ' Fill level:',iprcparm(mld_sub_fillin_) - case(mld_ilu_t_) - write(iout,*) ' Fill level:',iprcparm(mld_sub_fillin_) - if (present(rprcparm)) then - write(iout,*) ' Fill threshold :',rprcparm(mld_sub_iluthrs_) - else - write(iout,*) ' Fill threshold :',dprcparm(mld_sub_iluthrs_) - end if - case(mld_slu_,mld_umf_,mld_sludist_) - case default - write(iout,*) ' Should never get here!' - end select - write(iout,*) ' Overlap:',& - & iprcparm(mld_sub_ovr_) - write(iout,*) ' Restriction: ',& - & restrict_names(iprcparm(mld_sub_restr_)) - write(iout,*) ' Prolongation: ',& - & prolong_names(iprcparm(mld_sub_prol_)) - end select - return - end subroutine mld_base_prec_descr - - subroutine mld_ml_alg_descr(iout,ilev,iprcparm, info,rprcparm,dprcparm) - implicit none - integer, intent(in) :: iprcparm(:),iout,ilev - integer, intent(out) :: info - real(psb_spk_), intent(in), optional :: rprcparm(:) - real(psb_dpk_), intent(in), optional :: dprcparm(:) - - info = 0 - if (count((/ present(rprcparm),present(dprcparm) /)) /= 1) then - info=581 -!!$ call psb_errpush(info,name,a_err=" rprcparm, dprcparm") - return - endif - - if (iprcparm(mld_ml_type_)>mld_no_ml_) then - - write(iout,*) ' Multilevel type: ',& - & ml_names(iprcparm(mld_ml_type_)) - write(iout,*) ' Smoother position: ',& - & smooth_names(iprcparm(mld_smoother_pos_)) - write(iout,*) ' Aggregation: ', & - & aggr_names(iprcparm(mld_aggr_alg_)) - write(iout,*) ' Aggregation type: ', & - & aggr_kinds(iprcparm(mld_aggr_kind_)) - if (present(rprcparm)) then - write(iout,*) ' Aggregation threshold: ', & - & rprcparm(mld_aggr_thresh_) - else - write(iout,*) ' Aggregation threshold: ', & - & dprcparm(mld_aggr_thresh_) - end if - if (iprcparm(mld_aggr_kind_) /= mld_no_smooth_) then - if (iprcparm(mld_aggr_omega_alg_) == mld_eig_est_) then - write(iout,*) ' Damping omega computation: spectral radius estimate' - write(iout,*) ' Spectral radius estimate: ', & - & eigen_estimates(iprcparm(mld_aggr_eig_)) - else if (iprcparm(mld_aggr_omega_alg_) == mld_user_choice_) then - write(iout,*) ' Damping omega computation: user defined value.' - else - write(iout,*) ' Damping omega computation: unknown value in iprcparm!!' - end if - end if - end if - - return - end subroutine mld_ml_alg_descr - - subroutine mld_ml_level_descr(iout,ilev,iprcparm,nlaggr, info,rprcparm,dprcparm) - implicit none - integer, intent(in) :: iprcparm(:),iout,ilev - integer, intent(in), allocatable :: nlaggr(:) - integer, intent(out) :: info - real(psb_spk_), intent(in), optional :: rprcparm(:) - real(psb_dpk_), intent(in), optional :: dprcparm(:) - - info = 0 - if (count((/ present(rprcparm),present(dprcparm) /)) /= 1) then - info=581 -!!$ call psb_errpush(info,name,a_err=" rprcparm, dprcparm") - return - endif - - if (iprcparm(mld_ml_type_)>mld_no_ml_) then - write(iout,*) ' Level ',ilev - if (allocated(nlaggr)) then - write(iout,*) ' Size of coarse matrix: ', & - & sum(nlaggr(:)) - write(iout,*) ' Sizes of aggregates: ', & - & nlaggr(:) - end if - if (iprcparm(mld_aggr_kind_) /= mld_no_smooth_) then - if (present(rprcparm)) then - write(iout,*) ' Damping omega: ', & - & rprcparm(mld_aggr_omega_val_) - else - write(iout,*) ' Damping omega: ', & - & dprcparm(mld_aggr_omega_val_) - end if - end if - end if - - return - end subroutine mld_ml_level_descr - - subroutine mld_ml_coarse_descr(iout,ilev,iprcparm,iprcparm2,nlaggr,info,& - & rprcparm,dprcparm, rprcparm2,dprcparm2) - implicit none - integer, intent(in) :: iprcparm(:),iprcparm2(:),iout,ilev - integer, intent(in), allocatable :: nlaggr(:) - integer, intent(out) :: info - real(psb_spk_), intent(in), optional :: rprcparm(:), rprcparm2(:) - real(psb_dpk_), intent(in), optional :: dprcparm(:), dprcparm2(:) - - info = 0 - if (count((/ present(rprcparm),present(dprcparm) /)) /= 1) then - info=581 -!!$ call psb_errpush(info,name,a_err=" rprcparm, dprcparm") - return - endif - if (count((/ present(rprcparm2),present(dprcparm2) /)) /= 1) then - info=581 -!!$ call psb_errpush(info,name,a_err=" rprcparm, dprcparm") - return - endif - - if (iprcparm(mld_ml_type_)>mld_no_ml_) then - - write(iout,*) ' Level ',ilev,' (coarsest)' - write(iout,*) ' Coarsest matrix: ',& - & matrix_names(iprcparm(mld_coarse_mat_)) - if (allocated(nlaggr)) then - write(iout,*) ' Size of coarsest matrix: ', & - & sum( nlaggr(:)) - write(iout,*) ' Sizes of aggregates: ', & - & nlaggr(:) - end if - if (iprcparm(mld_aggr_kind_) /= mld_no_smooth_) then - - if (present(rprcparm)) then - write(iout,*) ' Damping omega: ', & - & rprcparm(mld_aggr_omega_val_) - else - write(iout,*) ' Damping omega: ', & - & dprcparm(mld_aggr_omega_val_) - end if - end if - if (iprcparm(mld_coarse_mat_) == mld_distr_mat_ .and. & - & iprcparm(mld_sub_solve_) /= mld_sludist_) then - write(iout,*) ' Coarsest matrix solver: block Jacobi with ', & - & fact_names(iprcparm2(mld_sub_solve_)) - write(iout,*) ' Number of Jacobi sweeps: ', & - & (iprcparm2(mld_smoother_sweeps_)) - else - write(iout,*) ' Coarsest matrix solver: ', & - & fact_names(iprcparm2(mld_sub_solve_)) - end if - select case(iprcparm2(mld_sub_solve_)) - case(mld_ilu_n_,mld_milu_n_) - write(iout,*) ' Fill level:',iprcparm2(mld_sub_fillin_) - case(mld_ilu_t_) - write(iout,*) ' Fill level:',iprcparm2(mld_sub_fillin_) - if (present(rprcparm2)) then - write(iout,*) ' Fill threshold :',rprcparm2(mld_sub_iluthrs_) - else if (present(dprcparm2)) then - write(iout,*) ' Fill threshold :',dprcparm2(mld_sub_iluthrs_) - end if - case(mld_slu_,mld_umf_,mld_sludist_) - case default - write(iout,*) ' Should never get here!' - end select - end if - - - return - end subroutine mld_ml_coarse_descr - - - ! - ! Subroutine: mld_file_prec_descr - ! Version: real - ! - ! This routine prints a description of the preconditioner to the standard - ! output or to a file. It must be called after the preconditioner has been - ! built by mld_precbld. - ! - ! Arguments: - ! p - type(mld_dprec_type), input. - ! The preconditioner data structure to be printed out. - ! info - integer, output. - ! error code. - ! iout - integer, input, optional. - ! The id of the file where the preconditioner description - ! will be printed. If iout is not present, then the standard - ! output is condidered. - ! - subroutine mld_file_prec_descr(p,info,iout) - implicit none - ! Arguments - type(mld_dprec_type), intent(in) :: p - integer, intent(out) :: info - integer, intent(in), optional :: iout - - ! Local variables - integer :: ilev, nlev - integer :: ictxt, me, np - character(len=20), parameter :: name='mld_file_prec_descr' - integer :: iout_ - - info = 0 - if (present(iout)) then - iout_ = iout - else - iout_ = 6 - end if - if (iout_ < 0) iout_ = 6 - - if (allocated(p%precv)) then - ictxt = psb_cd_get_context(p%precv(1)%prec%desc_data) - - call psb_info(ictxt,me,np) - - ! - ! The preconditioner description is printed by processor psb_root_. - ! This agrees with the fact that all the parameters defining the - ! preconditioner have the same values on all the procs (this is - ! ensured by mld_precbld). - ! - if (me==psb_root_) then - - write(iout_,*) - write(iout_,'(a)') 'Preconditioner description' - nlev = size(p%precv) - if (nlev >= 1) then - ! - ! Print description of base preconditioner - ! - - write(iout_,*) ' ' - - if (nlev > 1) then - write(iout_,*) 'Multilevel Schwarz' - write(iout_,*) - write(iout_,*) 'Base preconditioner (smoother) details' - endif - - ilev = 1 - call mld_base_prec_descr(iout_,p%precv(ilev)%prec%iprcparm,info,& - & dprcparm=p%precv(ilev)%prec%rprcparm) - - end if - - if (nlev > 1) then - - ! - ! Print multilevel details - ! - write(iout_,*) - write(iout_,*) 'Multilevel details' - - do ilev = 2, nlev - if (.not.allocated(p%precv(ilev)%iprcparm)) then - info = 3111 - write(iout_,*) ' ',name,': error: inconsistent MLPREC part, should call MLD_PRECINIT' - return - endif - end do - - write(iout_,*) ' Number of levels: ',nlev - - ! - ! Currently, all the preconditioner parameters must have the same value at levels - ! 2,...,nlev-1, hence only the values at level 2 are printed - ! - - ilev=2 - call mld_ml_alg_descr(iout_,ilev,p%precv(ilev)%iprcparm, info,& - & dprcparm=p%precv(ilev)%rprcparm) - - ! - ! Coarse matrices are different at levels 2,...,nlev-1, hence related - ! info is printed separately - ! - write(iout_,*) - do ilev = 2, nlev-1 - call mld_ml_level_descr(iout_,ilev,p%precv(ilev)%iprcparm,& - & p%precv(ilev)%map%naggr,info,& - & dprcparm=p%precv(ilev)%rprcparm) - end do - - ! - ! Print coarsest level details - ! - - ilev = nlev - write(iout_,*) - call mld_ml_coarse_descr(iout_,ilev,& - & p%precv(ilev)%iprcparm,p%precv(ilev)%prec%iprcparm,& - & p%precv(ilev)%map%naggr,info,& - & dprcparm=p%precv(ilev)%rprcparm,& - & dprcparm2=p%precv(ilev)%prec%rprcparm) - end if - - endif - write(iout_,*) - else - write(iout_,*) trim(name), & - & ': Error: no base preconditioner available, something is wrong!' - info = -2 - return - endif - - - end subroutine mld_file_prec_descr - - subroutine mld_sfile_prec_descr(p,info,iout) - implicit none - - ! Arguments - type(mld_sprec_type), intent(in) :: p - integer, intent(out) :: info - integer, intent(in), optional :: iout - - ! Local variables - integer :: ilev, nlev - integer :: ictxt, me, np - character(len=20), parameter :: name='mld_file_prec_descr' - integer :: iout_ - - info = 0 - if (present(iout)) then - iout_ = iout - else - iout_ = 6 - end if - if (iout_ < 0) iout_ = 6 - - if (allocated(p%precv)) then - ictxt = psb_cd_get_context(p%precv(1)%prec%desc_data) - - call psb_info(ictxt,me,np) - - ! - ! The preconditioner description is printed by processor psb_root_. - ! This agrees with the fact that all the parameters defining the - ! preconditioner have the same values on all the procs (this is - ! ensured by mld_precbld). - ! - if (me==psb_root_) then - - write(iout_,*) - write(iout_,*) 'Preconditioner description' - nlev = size(p%precv) - if (nlev >= 1) then - ! - ! Print description of base preconditioner - ! - - write(iout_,*) ' ' - - if (nlev > 1) then - write(iout_,*) 'Multilevel Schwarz' - write(iout_,*) - write(iout_,*) 'Base preconditioner (smoother) details' - endif - - ilev = 1 - call mld_base_prec_descr(iout_,p%precv(ilev)%prec%iprcparm,info,& - & rprcparm=p%precv(ilev)%prec%rprcparm) - - end if - - if (nlev > 1) then - - ! - ! Print multilevel details - ! - write(iout_,*) - write(iout_,*) 'Multilevel details' - - do ilev = 2, nlev - if (.not.allocated(p%precv(ilev)%iprcparm)) then - info = 3111 - write(iout_,*) ' ',name,': error: inconsistent MLPREC part, should call MLD_PRECINIT' - return - endif - end do - - write(iout_,*) ' Number of levels: ',nlev - - ! - ! Currently, all the preconditioner parameters must have the same value at levels - ! 2,...,nlev-1, hence only the values at level 2 are printed - ! - - ilev=2 - call mld_ml_alg_descr(iout_,ilev,p%precv(ilev)%iprcparm, info,& - & rprcparm=p%precv(ilev)%rprcparm) - - ! - ! Coarse matrices are different at levels 2,...,nlev-1, hence related - ! info is printed separately - ! - write(iout_,*) - do ilev = 2, nlev-1 - call mld_ml_level_descr(iout_,ilev,p%precv(ilev)%iprcparm,& - & p%precv(ilev)%map%naggr,info,& - & rprcparm=p%precv(ilev)%rprcparm) - end do - - ! - ! Print coarsest level details - ! - - ilev = nlev - write(iout_,*) - call mld_ml_coarse_descr(iout_,ilev,& - & p%precv(ilev)%iprcparm,p%precv(ilev)%prec%iprcparm,& - & p%precv(ilev)%map%naggr,info,& - & rprcparm=p%precv(ilev)%rprcparm, & - & rprcparm2=p%precv(ilev)%prec%rprcparm) - - end if - - endif - write(iout_,*) - else - - write(iout_,*) trim(name), & - & ': Error: no base preconditioner available, something is wrong!' - info = -2 - return - endif - - - end subroutine mld_sfile_prec_descr - - - function mld_prec_short_descr(p) - implicit none - type(mld_dprec_type), intent(in) :: p - character(len=20) :: mld_prec_short_descr - mld_prec_short_descr = ' ' - end function mld_prec_short_descr - - - ! - ! Subroutine: mld_zfile_prec_descr - ! Version: complex - ! - ! This routine prints to a file a description of the preconditioner. - ! - ! Arguments: - ! p - type(mld_zprec_type), input. - ! The preconditioner data structure to be printed out. - ! iout - integer, input. - ! The id of the file where the preconditioner description - ! will be printed. - ! - subroutine mld_zfile_prec_descr(p,info,iout) - implicit none - - ! Arguments - type(mld_zprec_type), intent(in) :: p - integer, intent(out) :: info - integer, intent(in), optional :: iout - - ! Local variables - integer :: ilev, nlev - integer :: ictxt, me, np - character(len=20), parameter :: name='mld_file_prec_descr' - integer :: iout_ - - info = 0 - if (present(iout)) then - iout_ = iout - else - iout_ = 6 - end if - if (iout_ < 0) iout_ = 6 - - if (allocated(p%precv)) then - ictxt = psb_cd_get_context(p%precv(1)%prec%desc_data) - - call psb_info(ictxt,me,np) - - ! - ! The preconditioner description is printed by processor psb_root_. - ! This agrees with the fact that all the parameters defining the - ! preconditioner have the same values on all the procs (this is - ! ensured by mld_precbld). - ! - if (me==psb_root_) then - - write(iout_,*) - write(iout_,*) 'Preconditioner description' - nlev = size(p%precv) - if (nlev >= 1) then - ! - ! Print description of base preconditioner - ! - - write(iout_,*) ' ' - - if (nlev > 1) then - write(iout_,*) 'Multilevel Schwarz' - write(iout_,*) - write(iout_,*) 'Base preconditioner (smoother) details' - endif - - ilev = 1 - call mld_base_prec_descr(iout_,p%precv(ilev)%prec%iprcparm,info,& - & dprcparm=p%precv(ilev)%prec%rprcparm) - - end if - - if (nlev > 1) then - - ! - ! Print multilevel details - ! - write(iout_,*) - write(iout_,*) 'Multilevel details' - - do ilev = 2, nlev - if (.not.allocated(p%precv(ilev)%iprcparm)) then - info = 3111 - write(iout_,*) ' ',name,': error: inconsistent MLPREC part, should call MLD_PRECINIT' - return - endif - end do - - write(iout_,*) ' Number of levels: ',nlev - - ! - ! Currently, all the preconditioner parameters must have the same value at levels - ! 2,...,nlev-1, hence only the values at level 2 are printed - ! - - ilev=2 - call mld_ml_alg_descr(iout_,ilev,p%precv(ilev)%iprcparm, info,& - & dprcparm=p%precv(ilev)%rprcparm) - - ! - ! Coarse matrices are different at levels 2,...,nlev-1, hence related - ! info is printed separately - ! - write(iout_,*) - do ilev = 2, nlev-1 - call mld_ml_level_descr(iout_,ilev,p%precv(ilev)%iprcparm,& - & p%precv(ilev)%map%naggr,info,& - & dprcparm=p%precv(ilev)%rprcparm) - end do - - ! - ! Print coarsest level details - ! - - ilev = nlev - write(iout_,*) - call mld_ml_coarse_descr(iout_,ilev,& - & p%precv(ilev)%iprcparm,p%precv(ilev)%prec%iprcparm,& - & p%precv(ilev)%map%naggr,info,& - & dprcparm=p%precv(ilev)%rprcparm,& - & dprcparm2=p%precv(ilev)%prec%rprcparm) - end if - - endif - write(iout_,*) - else - write(iout_,*) trim(name), & - & ': Error: no base preconditioner available, something is wrong!' - info = -2 - return - endif - - end subroutine mld_zfile_prec_descr - - subroutine mld_cfile_prec_descr(p,info,iout) - implicit none - - ! Arguments - type(mld_cprec_type), intent(in) :: p - integer, intent(out) :: info - integer, intent(in), optional :: iout - - ! Local variables - integer :: ilev, nlev - integer :: ictxt, me, np - character(len=20), parameter :: name='mld_file_prec_descr' - integer :: iout_ - - info = 0 - if (present(iout)) then - iout_ = iout - else - iout_ = 6 - end if - if (iout_ < 0) iout_ = 6 - - if (allocated(p%precv)) then - ictxt = psb_cd_get_context(p%precv(1)%prec%desc_data) - - call psb_info(ictxt,me,np) - - ! - ! The preconditioner description is printed by processor psb_root_. - ! This agrees with the fact that all the parameters defining the - ! preconditioner have the same values on all the procs (this is - ! ensured by mld_precbld). - ! - if (me==psb_root_) then - write(iout_,*) - write(iout_,*) 'Preconditioner description' - nlev = size(p%precv) - if (nlev >= 1) then - ! - ! Print description of base preconditioner - ! - - write(iout_,*) ' ' - - if (nlev > 1) then - write(iout_,*) 'Multilevel Schwarz' - write(iout_,*) - write(iout_,*) 'Base preconditioner (smoother) details' - endif - - ilev = 1 - call mld_base_prec_descr(iout_,p%precv(ilev)%prec%iprcparm,info,& - & rprcparm=p%precv(ilev)%prec%rprcparm) - - end if - - if (nlev > 1) then - - ! - ! Print multilevel details - ! - write(iout_,*) - write(iout_,*) 'Multilevel details' - - do ilev = 2, nlev - if (.not.allocated(p%precv(ilev)%iprcparm)) then - info = 3111 - write(iout_,*) ' ',name,': error: inconsistent MLPREC part, should call MLD_PRECINIT' - return - endif - end do - - write(iout_,*) ' Number of levels: ',nlev - - ! - ! Currently, all the preconditioner parameters must have the same value at levels - ! 2,...,nlev-1, hence only the values at level 2 are printed - ! - - ilev=2 - call mld_ml_alg_descr(iout_,ilev,p%precv(ilev)%iprcparm, info,& - & rprcparm=p%precv(ilev)%rprcparm) - - ! - ! Coarse matrices are different at levels 2,...,nlev-1, hence related - ! info is printed separately - ! - write(iout_,*) - - do ilev = 2, nlev-1 - call mld_ml_level_descr(iout_,ilev,p%precv(ilev)%iprcparm,& - & p%precv(ilev)%map%naggr,info,& - & rprcparm=p%precv(ilev)%rprcparm) - end do - - ! - ! Print coarsest level details - ! - - ilev = nlev - write(iout_,*) - call mld_ml_coarse_descr(iout_,ilev,& - & p%precv(ilev)%iprcparm,p%precv(ilev)%prec%iprcparm,& - & p%precv(ilev)%map%naggr,info,& - & rprcparm=p%precv(ilev)%rprcparm,& - & rprcparm2=p%precv(ilev)%prec%rprcparm) - end if - - endif - write(iout_,*) - else - write(iout_,*) trim(name), & - & ': Error: no base preconditioner available, something is wrong!' - info = -2 - return - endif - - end subroutine mld_cfile_prec_descr - - - function mld_zprec_short_descr(p) - implicit none - type(mld_zprec_type), intent(in) :: p - character(len=20) :: mld_zprec_short_descr - mld_zprec_short_descr = ' ' - - end function mld_zprec_short_descr - - - ! - ! Functions/subroutines checking if the preconditioner is correctly defined - ! - - function is_legal_base_prec(ip) - implicit none - integer, intent(in) :: ip - logical :: is_legal_base_prec - - is_legal_base_prec = ((ip>=mld_noprec_).and.(ip<=mld_max_prec_)) - return - end function is_legal_base_prec - function is_legal_n_ovr(ip) - implicit none - integer, intent(in) :: ip - logical :: is_legal_n_ovr - - is_legal_n_ovr = (ip >= 0) - return - end function is_legal_n_ovr - function is_legal_renum(ip) - implicit none - integer, intent(in) :: ip - logical :: is_legal_renum - is_legal_renum = ((ip >= 0).and.(ip <= mld_max_renum_)) - return - end function is_legal_renum - function is_legal_jac_sweeps(ip) - implicit none - integer, intent(in) :: ip - logical :: is_legal_jac_sweeps - - is_legal_jac_sweeps = (ip >= 1) - return - end function is_legal_jac_sweeps - function is_legal_prolong(ip) - implicit none - integer, intent(in) :: ip - logical :: is_legal_prolong - is_legal_prolong = ((ip>=psb_none_).and.(ip<=psb_square_root_)) - return - end function is_legal_prolong - function is_legal_restrict(ip) - implicit none - integer, intent(in) :: ip - logical :: is_legal_restrict - is_legal_restrict = ((ip==psb_nohalo_).or.(ip==psb_halo_)) - return - end function is_legal_restrict - function is_legal_ml_type(ip) - implicit none - integer, intent(in) :: ip - logical :: is_legal_ml_type - - is_legal_ml_type = ((ip>=mld_no_ml_).and.(ip<=mld_max_ml_type_)) - return - end function is_legal_ml_type - function is_legal_ml_aggr_alg(ip) - implicit none - integer, intent(in) :: ip - logical :: is_legal_ml_aggr_alg - - is_legal_ml_aggr_alg = ((ip>=mld_dec_aggr_).and.(ip<=mld_max_aggr_alg_)) - return - end function is_legal_ml_aggr_alg - function is_legal_ml_aggr_omega_alg(ip) - implicit none - integer, intent(in) :: ip - logical :: is_legal_ml_aggr_omega_alg - - is_legal_ml_aggr_omega_alg = ((ip==mld_eig_est_).or.(ip==mld_user_choice_)) - return - end function is_legal_ml_aggr_omega_alg - function is_legal_ml_aggr_eig(ip) - implicit none - integer, intent(in) :: ip - logical :: is_legal_ml_aggr_eig - - is_legal_ml_aggr_eig = (ip==mld_max_norm_) - return - end function is_legal_ml_aggr_eig - function is_legal_ml_smooth_pos(ip) - implicit none - integer, intent(in) :: ip - logical :: is_legal_ml_smooth_pos - - is_legal_ml_smooth_pos = ((ip>=mld_pre_smooth_).and.(ip<=mld_max_smooth_)) - return - end function is_legal_ml_smooth_pos - function is_legal_ml_aggr_kind(ip) - implicit none - integer, intent(in) :: ip - logical :: is_legal_ml_aggr_kind - - is_legal_ml_aggr_kind = ((ip>=0).and.(ip<=mld_max_aggr_kind_)) - return - end function is_legal_ml_aggr_kind - function is_legal_ml_coarse_mat(ip) - implicit none - integer, intent(in) :: ip - logical :: is_legal_ml_coarse_mat - - is_legal_ml_coarse_mat = ((ip>=0).and.(ip<=mld_max_coarse_mat_)) - return - end function is_legal_ml_coarse_mat - function is_legal_aggr_filter(ip) - implicit none - integer, intent(in) :: ip - logical :: is_legal_aggr_filter - - is_legal_aggr_filter = ((ip>=0).and.(ip<=mld_max_filter_mat_)) - return - end function is_legal_aggr_filter - function is_distr_ml_coarse_mat(ip) - implicit none - integer, intent(in) :: ip - logical :: is_distr_ml_coarse_mat - - is_distr_ml_coarse_mat = (ip==mld_distr_mat_) - return - end function is_distr_ml_coarse_mat - function is_legal_ml_fact(ip) - implicit none - integer, intent(in) :: ip - logical :: is_legal_ml_fact - ! Here the minimum is really 1, mld_fact_none_ is not acceptable. - is_legal_ml_fact = ((ip>=1).and.(ip<=mld_max_sub_solve_)) - return - end function is_legal_ml_fact - function is_legal_ml_lev(ip) - implicit none - integer, intent(in) :: ip - logical :: is_legal_ml_lev - - is_legal_ml_lev = (ip >= 0) - return - end function is_legal_ml_lev - function is_legal_omega(ip) - implicit none - real(psb_dpk_), intent(in) :: ip - logical :: is_legal_omega - is_legal_omega = ((ip>=0.0d0).and.(ip<=2.0d0)) - return - end function is_legal_omega - function is_legal_fact_thrs(ip) - implicit none - real(psb_dpk_), intent(in) :: ip - logical :: is_legal_fact_thrs - - is_legal_fact_thrs = (ip>=0.0d0) - return - end function is_legal_fact_thrs - function is_legal_aggr_thrs(ip) - implicit none - real(psb_dpk_), intent(in) :: ip - logical :: is_legal_aggr_thrs - - is_legal_aggr_thrs = (ip>=0.0d0) - return - end function is_legal_aggr_thrs - - function is_legal_s_omega(ip) - implicit none - real(psb_spk_), intent(in) :: ip - logical :: is_legal_s_omega - is_legal_s_omega = ((ip>=0.0).and.(ip<=2.0)) - return - end function is_legal_s_omega - function is_legal_s_fact_thrs(ip) - implicit none - real(psb_spk_), intent(in) :: ip - logical :: is_legal_s_fact_thrs - - is_legal_s_fact_thrs = (ip>=0.0) - return - end function is_legal_s_fact_thrs - function is_legal_s_aggr_thrs(ip) - implicit none - real(psb_spk_), intent(in) :: ip - logical :: is_legal_s_aggr_thrs - - is_legal_s_aggr_thrs = (ip>=0.0) - return - end function is_legal_s_aggr_thrs - - - subroutine mld_icheck_def(ip,name,id,is_legal) - implicit none - integer, intent(inout) :: ip - integer, intent(in) :: id - character(len=*), intent(in) :: name - interface - function is_legal(i) - integer, intent(in) :: i - logical :: is_legal - end function is_legal - end interface - character(len=20), parameter :: rname='mld_check_def' - - if (.not.is_legal(ip)) then - write(0,*)trim(rname),': Error: Illegal value for ',& - & name,' :',ip, '. defaulting to ',id - ip = id - end if - end subroutine mld_icheck_def - - subroutine mld_scheck_def(ip,name,id,is_legal) - implicit none - real(psb_spk_), intent(inout) :: ip - real(psb_spk_), intent(in) :: id - character(len=*), intent(in) :: name - interface - function is_legal(i) - use psb_base_mod, only : psb_spk_ - real(psb_spk_), intent(in) :: i - logical :: is_legal - end function is_legal - end interface - character(len=20), parameter :: rname='mld_check_def' - - if (.not.is_legal(ip)) then - write(0,*)trim(rname),': Error: Illegal value for ',& - & name,' :',ip, '. defaulting to ',id - ip = id - end if - end subroutine mld_scheck_def - - subroutine mld_dcheck_def(ip,name,id,is_legal) - implicit none - real(psb_dpk_), intent(inout) :: ip - real(psb_dpk_), intent(in) :: id - character(len=*), intent(in) :: name - interface - function is_legal(i) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_), intent(in) :: i - logical :: is_legal - end function is_legal - end interface - character(len=20), parameter :: rname='mld_check_def' - - if (.not.is_legal(ip)) then - write(0,*)trim(rname),': Error: Illegal value for ',& - & name,' :',ip, '. defaulting to ',id - ip = id - end if - end subroutine mld_dcheck_def - - ! - ! Subroutines: mld_Tbase_precfree, mld_T_onelev_precfree, mld_Tprec_free - ! Version: real/complex - ! - ! These routines deallocate the mld_Tbaseprec_type, mld_Tonelev_type and - ! mld_Tprec_type data structures. - ! - ! Arguments: - ! p - type(mld_Tbaseprec_type/mld_Tonelev_type/mld_Tprec_type), input. - ! The data structure to be deallocated. - ! info - integer, output. - ! error code. - ! - - subroutine mld_sbase_precfree(p,info) - implicit none - - type(mld_sbaseprec_type), intent(inout) :: p - integer, intent(out) :: info - integer :: i - - info = 0 - - ! Actually we might just deallocate the top level array, except - ! for the inner UMFPACK or SLU stuff - - if (allocated(p%d)) then - deallocate(p%d,stat=info) - end if - - if (allocated(p%av)) then - do i=1,size(p%av) - call p%av(i)%free() - enddo - deallocate(p%av,stat=info) - end if - - if (allocated(p%desc_data%matrix_data)) & - & call psb_cdfree(p%desc_data,info) - - if (allocated(p%rprcparm)) then - deallocate(p%rprcparm,stat=info) - end if - - - if (allocated(p%perm)) then - deallocate(p%perm,stat=info) - endif - - if (allocated(p%invperm)) then - deallocate(p%invperm,stat=info) - endif - - if (allocated(p%iprcparm)) then - if (p%iprcparm(mld_prec_status_) == mld_prec_built_) then - if (p%iprcparm(mld_sub_solve_)==mld_slu_) then - call mld_sslu_free(p%iprcparm(mld_slu_ptr_),info) - end if - end if - deallocate(p%iprcparm,stat=info) - end if - call mld_nullify_baseprec(p) - end subroutine mld_sbase_precfree - - - subroutine mld_s_onelev_precfree(p,info) - implicit none - - type(mld_sonelev_type), intent(inout) :: p - integer, intent(out) :: info - integer :: i - - info = 0 - - ! Actually we might just deallocate the top level array, except - ! for the inner UMFPACK or SLU stuff - call mld_precfree(p%prec,info) - - call p%ac%free() - if (allocated(p%desc_ac%matrix_data)) & - & call psb_cdfree(p%desc_ac,info) - - if (allocated(p%rprcparm)) then - deallocate(p%rprcparm,stat=info) - end if - ! This is a pointer to something else, must not free it here. - nullify(p%base_a) - ! This is a pointer to something else, must not free it here. - nullify(p%base_desc) - - ! - ! free explicitly map??? - ! For now thanks to allocatable semantics - ! works anyway. - ! - - call mld_nullify_onelevprec(p) - end subroutine mld_s_onelev_precfree - - - subroutine mld_nullify_s_onelevprec(p) - implicit none - - type(mld_sonelev_type), intent(inout) :: p - - nullify(p%base_a) - nullify(p%base_desc) - - end subroutine mld_nullify_s_onelevprec - - subroutine mld_nullify_sbaseprec(p) - implicit none - - type(mld_sbaseprec_type), intent(inout) :: p - - - end subroutine mld_nullify_sbaseprec - - - subroutine mld_dbase_precfree(p,info) - implicit none - - type(mld_dbaseprec_type), intent(inout) :: p - integer, intent(out) :: info - integer :: i - - info = 0 - - ! Actually we might just deallocate the top level array, except - ! for the inner UMFPACK or SLU stuff - - if (allocated(p%d)) then - deallocate(p%d,stat=info) - end if - - if (allocated(p%av)) then - do i=1,size(p%av) - call p%av(i)%free() - enddo - deallocate(p%av,stat=info) - end if - - if (allocated(p%desc_data%matrix_data)) & - & call psb_cdfree(p%desc_data,info) - - if (allocated(p%rprcparm)) then - deallocate(p%rprcparm,stat=info) - end if - - if (allocated(p%perm)) then - deallocate(p%perm,stat=info) - endif - - if (allocated(p%invperm)) then - deallocate(p%invperm,stat=info) - endif - - if (allocated(p%iprcparm)) then - if (p%iprcparm(mld_prec_status_) == mld_prec_built_) then - if (p%iprcparm(mld_sub_solve_)==mld_slu_) then - call mld_dslu_free(p%iprcparm(mld_slu_ptr_),info) - end if - if (p%iprcparm(mld_sub_solve_)==mld_sludist_) then - call mld_dsludist_free(p%iprcparm(mld_slud_ptr_),info) - end if - if (p%iprcparm(mld_sub_solve_)==mld_umf_) then - call mld_dumf_free(p%iprcparm(mld_umf_symptr_),& - & p%iprcparm(mld_umf_numptr_),info) - end if - end if - deallocate(p%iprcparm,stat=info) - end if - call mld_nullify_baseprec(p) - end subroutine mld_dbase_precfree - - subroutine mld_d_onelev_precfree(p,info) - implicit none - - type(mld_donelev_type), intent(inout) :: p - integer, intent(out) :: info - integer :: i - - info = 0 - - ! Actually we might just deallocate the top level array, except - ! for the inner UMFPACK or SLU stuff - call mld_precfree(p%prec,info) - - call p%ac%free() - if (allocated(p%desc_ac%matrix_data)) & - & call psb_cdfree(p%desc_ac,info) - - if (allocated(p%rprcparm)) then - deallocate(p%rprcparm,stat=info) - end if - ! This is a pointer to something else, must not free it here. - nullify(p%base_a) - ! This is a pointer to something else, must not free it here. - nullify(p%base_desc) - - ! - ! free explicitly map??? - ! For now thanks to allocatable semantics - ! works anyway. - ! - - call mld_nullify_onelevprec(p) - end subroutine mld_d_onelev_precfree - - subroutine mld_nullify_dbaseprec(p) - implicit none - - type(mld_dbaseprec_type), intent(inout) :: p - - - end subroutine mld_nullify_dbaseprec - - subroutine mld_nullify_d_onelevprec(p) - implicit none - - type(mld_donelev_type), intent(inout) :: p - - nullify(p%base_a) - nullify(p%base_desc) - - end subroutine mld_nullify_d_onelevprec - - subroutine mld_cbase_precfree(p,info) - implicit none - type(mld_cbaseprec_type), intent(inout) :: p - integer, intent(out) :: info - integer :: i - - info = 0 - - if (allocated(p%d)) then - deallocate(p%d,stat=info) - end if - - if (allocated(p%av)) then - do i=1,size(p%av) - call p%av(i)%free() - enddo - deallocate(p%av,stat=info) - - end if - if (allocated(p%desc_data%matrix_data)) & - & call psb_cdfree(p%desc_data,info) - - if (allocated(p%rprcparm)) then - deallocate(p%rprcparm,stat=info) - end if - - if (allocated(p%perm)) then - deallocate(p%perm,stat=info) - endif - - if (allocated(p%invperm)) then - deallocate(p%invperm,stat=info) - endif - - if (allocated(p%iprcparm)) then - if (p%iprcparm(mld_prec_status_) == mld_prec_built_) then - if (p%iprcparm(mld_sub_solve_)==mld_slu_) then - call mld_cslu_free(p%iprcparm(mld_slu_ptr_),info) - end if - end if - deallocate(p%iprcparm,stat=info) - end if - call mld_nullify_baseprec(p) - end subroutine mld_cbase_precfree - - subroutine mld_c_onelev_precfree(p,info) - implicit none - - type(mld_conelev_type), intent(inout) :: p - integer, intent(out) :: info - integer :: i - - info = 0 - - ! Actually we might just deallocate the top level array, except - ! for the inner UMFPACK or SLU stuff - call mld_precfree(p%prec,info) - - call p%ac%free() - if (allocated(p%desc_ac%matrix_data)) & - & call psb_cdfree(p%desc_ac,info) - - if (allocated(p%rprcparm)) then - deallocate(p%rprcparm,stat=info) - end if - ! This is a pointer to something else, must not free it here. - nullify(p%base_a) - ! This is a pointer to something else, must not free it here. - nullify(p%base_desc) - - ! - ! free explicitly map??? - ! For now thanks to allocatable semantics - ! works anyway. - ! - - call mld_nullify_onelevprec(p) - end subroutine mld_c_onelev_precfree - - subroutine mld_nullify_c_onelevprec(p) - implicit none - - type(mld_conelev_type), intent(inout) :: p - - nullify(p%base_a) - nullify(p%base_desc) - - end subroutine mld_nullify_c_onelevprec - - subroutine mld_nullify_cbaseprec(p) - implicit none - - type(mld_cbaseprec_type), intent(inout) :: p - - - end subroutine mld_nullify_cbaseprec - - subroutine mld_zbase_precfree(p,info) - implicit none - type(mld_zbaseprec_type), intent(inout) :: p - integer, intent(out) :: info - integer :: i - - info = 0 - - if (allocated(p%d)) then - deallocate(p%d,stat=info) - end if - - if (allocated(p%av)) then - do i=1,size(p%av) - call p%av(i)%free() - enddo - deallocate(p%av,stat=info) - - end if - if (allocated(p%desc_data%matrix_data)) & - & call psb_cdfree(p%desc_data,info) - - if (allocated(p%rprcparm)) then - deallocate(p%rprcparm,stat=info) - end if - - if (allocated(p%perm)) then - deallocate(p%perm,stat=info) - endif - - if (allocated(p%invperm)) then - deallocate(p%invperm,stat=info) - endif - - if (allocated(p%iprcparm)) then - if (p%iprcparm(mld_prec_status_) == mld_prec_built_) then - if (p%iprcparm(mld_sub_solve_)==mld_slu_) then - call mld_zslu_free(p%iprcparm(mld_slu_ptr_),info) - end if - if (p%iprcparm(mld_sub_solve_)==mld_umf_) then - call mld_zumf_free(p%iprcparm(mld_umf_symptr_),& - & p%iprcparm(mld_umf_numptr_),info) - end if - end if - deallocate(p%iprcparm,stat=info) - end if - call mld_nullify_baseprec(p) - end subroutine mld_zbase_precfree - - subroutine mld_z_onelev_precfree(p,info) - implicit none - - type(mld_zonelev_type), intent(inout) :: p - integer, intent(out) :: info - integer :: i - - info = 0 - - ! Actually we might just deallocate the top level array, except - ! for the inner UMFPACK or SLU stuff - call mld_precfree(p%prec,info) - - call p%ac%free() - if (allocated(p%desc_ac%matrix_data)) & - & call psb_cdfree(p%desc_ac,info) - - if (allocated(p%rprcparm)) then - deallocate(p%rprcparm,stat=info) - end if - ! This is a pointer to something else, must not free it here. - nullify(p%base_a) - ! This is a pointer to something else, must not free it here. - nullify(p%base_desc) - - ! - ! free explicitly map??? - ! For now thanks to allocatable semantics - ! works anyway. - ! - - call mld_nullify_onelevprec(p) - end subroutine mld_z_onelev_precfree - - subroutine mld_nullify_z_onelevprec(p) - implicit none - - type(mld_zonelev_type), intent(inout) :: p - - nullify(p%base_a) - nullify(p%base_desc) - - end subroutine mld_nullify_z_onelevprec - - - subroutine mld_nullify_zbaseprec(p) - implicit none - - type(mld_zbaseprec_type), intent(inout) :: p - - - end subroutine mld_nullify_zbaseprec - - - function pr_to_str(iprec) - implicit none - - integer, intent(in) :: iprec - character(len=10) :: pr_to_str - - select case(iprec) - case(mld_noprec_) - pr_to_str='NOPREC' - case(mld_diag_) - pr_to_str='DIAG' - case(mld_bjac_) - pr_to_str='BJAC' - case(mld_as_) - pr_to_str='AS' - end select - - end function pr_to_str - - subroutine mld_sprec_free(p,info) - - use psb_base_mod - - implicit none - - ! Arguments - type(mld_sprec_type), intent(inout) :: p - integer, intent(out) :: info - - ! Local variables - integer :: me,err_act,i - character(len=20) :: name - - if(psb_get_errstatus().ne.0) return - info=0 - name = 'mld_dprecfree' - call psb_erractionsave(err_act) - - me=-1 - - if (allocated(p%precv)) then - do i=1,size(p%precv) - call mld_precfree(p%precv(i),info) - end do - deallocate(p%precv) - end if - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act.eq.psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine mld_sprec_free - - subroutine mld_dprec_free(p,info) - - use psb_base_mod - - implicit none - - ! Arguments - type(mld_dprec_type), intent(inout) :: p - integer, intent(out) :: info - - ! Local variables - integer :: me,err_act,i - character(len=20) :: name - - if(psb_get_errstatus().ne.0) return - info=0 - name = 'mld_dprecfree' - call psb_erractionsave(err_act) - - me=-1 - - if (allocated(p%precv)) then - do i=1,size(p%precv) - call mld_precfree(p%precv(i),info) - end do - deallocate(p%precv) - end if - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act.eq.psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine mld_dprec_free - - subroutine mld_cprec_free(p,info) - - use psb_base_mod - - implicit none - - ! Arguments - type(mld_cprec_type), intent(inout) :: p - integer, intent(out) :: info - - ! Local variables - integer :: me,err_act,i - character(len=20) :: name - - if(psb_get_errstatus().ne.0) return - info=0 - name = 'mld_dprecfree' - call psb_erractionsave(err_act) - - me=-1 - - if (allocated(p%precv)) then - do i=1,size(p%precv) - call mld_precfree(p%precv(i),info) - end do - deallocate(p%precv) - end if - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act.eq.psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine mld_cprec_free - - - subroutine mld_zprec_free(p,info) - - use psb_base_mod - - implicit none - - ! Arguments - type(mld_zprec_type), intent(inout) :: p - integer, intent(out) :: info - - ! Local variables - integer :: me,err_act,i - character(len=20) :: name - - if(psb_get_errstatus().ne.0) return - info=0 - name = 'mld_dprecfree' - call psb_erractionsave(err_act) - - me=-1 - - if (allocated(p%precv)) then - do i=1,size(p%precv) - call mld_precfree(p%precv(i),info) - end do - deallocate(p%precv) - end if - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act.eq.psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine mld_zprec_free - - - - - - - subroutine mld_s_apply2v(prec,x,y,desc_data,info,trans,work) - use psb_base_mod - type(psb_desc_type),intent(in) :: desc_data - class(mld_sprec_type), intent(in) :: prec - real(psb_spk_),intent(in) :: x(:) - real(psb_spk_),intent(inout) :: y(:) - integer, intent(out) :: info - character(len=1), optional :: trans - real(psb_spk_),intent(inout), optional, target :: work(:) - Integer :: err_act - character(len=20) :: name='s_prec_apply' - - call psb_erractionsave(err_act) - - select type(prec) -!!$ type is (mld_sprec_type) -!!$ call mld_precaply(prec,x,y,desc_data,info,trans,work) - class default - info = 700 - call psb_errpush(info,name) - goto 9999 - end select - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine mld_s_apply2v - subroutine mld_s_apply1v(prec,x,desc_data,info,trans) - use psb_base_mod - type(psb_desc_type),intent(in) :: desc_data - class(mld_sprec_type), intent(in) :: prec - real(psb_spk_),intent(inout) :: x(:) - integer, intent(out) :: info - character(len=1), optional :: trans - Integer :: err_act - character(len=20) :: name='s_prec_apply' - - call psb_erractionsave(err_act) - - select type(prec) -!!$ type is (mld_sprec_type) -!!$ call mld_precaply(prec,x,desc_data,info,trans) - class default - info = 700 - call psb_errpush(info,name) - goto 9999 - end select - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - end subroutine mld_s_apply1v - - subroutine mld_d_apply2v(prec,x,y,desc_data,info,trans,work) - use psb_base_mod - type(psb_desc_type),intent(in) :: desc_data - class(mld_dprec_type), intent(in) :: prec - real(psb_dpk_),intent(in) :: x(:) - real(psb_dpk_),intent(inout) :: y(:) - integer, intent(out) :: info - character(len=1), optional :: trans - real(psb_dpk_),intent(inout), optional, target :: work(:) - Integer :: err_act - character(len=20) :: name='d_prec_apply' - - call psb_erractionsave(err_act) - - select type(prec) - type is (mld_dprec_type) - call mld_precaply(prec,x,y,desc_data,info,trans,work) - class default - info = 700 - call psb_errpush(info,name) - goto 9999 - end select - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine mld_d_apply2v - - subroutine mld_d_apply1v(prec,x,desc_data,info,trans) - use psb_base_mod - type(psb_desc_type),intent(in) :: desc_data - class(mld_dprec_type), intent(in) :: prec - real(psb_dpk_),intent(inout) :: x(:) - integer, intent(out) :: info - character(len=1), optional :: trans - Integer :: err_act - character(len=20) :: name='d_prec_apply' - - call psb_erractionsave(err_act) - - select type(prec) - type is (mld_dprec_type) - call mld_precaply(prec,x,desc_data,info,trans) - class default - info = 700 - call psb_errpush(info,name) - goto 9999 - end select - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine mld_d_apply1v - subroutine mld_c_apply2v(prec,x,y,desc_data,info,trans,work) - use psb_base_mod - - type(psb_desc_type),intent(in) :: desc_data - class(mld_cprec_type), intent(in) :: prec - complex(psb_spk_),intent(in) :: x(:) - complex(psb_spk_),intent(inout) :: y(:) - integer, intent(out) :: info - character(len=1), optional :: trans - complex(psb_spk_),intent(inout), optional, target :: work(:) - Integer :: err_act - character(len=20) :: name='s_prec_apply' - - call psb_erractionsave(err_act) - - select type(prec) -!!$ type is (mld_cprec_type) -!!$ call mld_precaply(prec,x,y,desc_data,info,trans,work) - class default - info = 700 - call psb_errpush(info,name) - goto 9999 - end select - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine mld_c_apply2v - subroutine mld_c_apply1v(prec,x,desc_data,info,trans) - use psb_base_mod - - type(psb_desc_type),intent(in) :: desc_data - class(mld_cprec_type), intent(in) :: prec - complex(psb_spk_),intent(inout) :: x(:) - integer, intent(out) :: info - character(len=1), optional :: trans - Integer :: err_act - character(len=20) :: name='c_prec_apply' - - call psb_erractionsave(err_act) - - select type(prec) -!!$ type is (mld_cprec_type) -!!$ call mld_precaply(prec,x,desc_data,info,trans) - class default - info = 700 - call psb_errpush(info,name) - goto 9999 - end select - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine mld_c_apply1v - - subroutine mld_z_apply2v(prec,x,y,desc_data,info,trans,work) - use psb_base_mod - type(psb_desc_type),intent(in) :: desc_data - class(mld_zprec_type), intent(in) :: prec - complex(psb_dpk_),intent(in) :: x(:) - complex(psb_dpk_),intent(inout) :: y(:) - integer, intent(out) :: info - character(len=1), optional :: trans - complex(psb_dpk_),intent(inout), optional, target :: work(:) - Integer :: err_act - character(len=20) :: name='z_prec_apply' - - call psb_erractionsave(err_act) - - select type(prec) -!!$ type is (mld_zprec_type) -!!$ call mld_precaply(prec,x,y,desc_data,info,trans,work) - class default - info = 700 - call psb_errpush(info,name) - goto 9999 - end select - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine mld_z_apply2v - subroutine mld_z_apply1v(prec,x,desc_data,info,trans) - use psb_base_mod - - type(psb_desc_type),intent(in) :: desc_data - class(mld_zprec_type), intent(in) :: prec - complex(psb_dpk_),intent(inout) :: x(:) - integer, intent(out) :: info - character(len=1), optional :: trans - Integer :: err_act - character(len=20) :: name='z_prec_apply' - - call psb_erractionsave(err_act) - - select type(prec) -!!$ type is (mld_zprec_type) -!!$ call mld_precaply(prec,x,desc_data,info,trans) - class default - info = 700 - call psb_errpush(info,name) - goto 9999 - end select - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine mld_z_apply1v + use mld_base_prec_type + use mld_s_prec_type + use mld_d_prec_type + use mld_c_prec_type + use mld_z_prec_type end module mld_prec_type diff --git a/mlprec/mld_s_prec_type.f90 b/mlprec/mld_s_prec_type.f90 new file mode 100644 index 00000000..4e14f5e4 --- /dev/null +++ b/mlprec/mld_s_prec_type.f90 @@ -0,0 +1,702 @@ +!!$ +!!$ +!!$ MLD2P4 version 1.1 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 2.3.1) +!!$ +!!$ (C) Copyright 2008,2009 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +! File: mld_prec_type.f90 +! +! Module: mld_prec_type +! +! This module defines: +! - the mld_prec_type data structure containing the preconditioner and related +! data structures; +! - integer constants defining the preconditioner; +! - character constants describing the preconditioner (used by the routines +! printing out a preconditioner description); +! - the interfaces to the routines for the management of the preconditioner +! data structure (see below). +! +! It contains routines for +! - converting character constants defining the preconditioner into integer +! constants; +! - checking if the preconditioner is correctly defined; +! - printing a description of the preconditioner; +! - deallocating the preconditioner data structure. +! + +module mld_s_prec_type + + use mld_base_prec_type + + ! + ! Type: mld_Tprec_type. + ! + ! It is the data type containing all the information about the multilevel + ! preconditioner (here and in the following 'T' denotes 'd', 's', 'c' and + ! 'z', according to the real/complex, single/double precision version of + ! MLD2P4). It consists of an array of 'one-level' intermediate data structures + ! of type mld_Tonelev_type, each containing the information needed to apply + ! the smoothing and the coarse-space correction at a generic level. + ! + ! type mld_Tprec_type + ! type(mld_Tonelev_type), allocatable :: precv(:) + ! end type mld_Tprec_type + ! + ! Note that the levels are numbered in increasing order starting from + ! the finest one and the number of levels is given by size(precv(:)). + ! + ! + ! Type: mld_Tonelev_type. + ! + ! It is the data type containing the necessary items for the current + ! level (essentially, the base preconditioner, the current-level matrix + ! and the restriction and prolongation operators). + ! + ! type mld_Tonelev_type + ! type(mld_Tbaseprec_type) :: prec + ! integer, allocatable :: iprcparm(:) + ! real(psb_Tpk_), allocatable :: rprcparm(:) + ! type(psb_T_sparse_mat) :: ac + ! type(psb_desc_type) :: desc_ac + ! type(psb_T_sparse_mat), pointer :: base_a => null() + ! type(psb_desc_type), pointer :: base_desc => null() + ! type(psb_Tlinmap_type) :: map + ! end type mld_Tonelev_type + ! + ! Note that psb_Tpk denotes the kind of the real data type to be chosen + ! according to single/double precision version of MLD2P4. + ! + ! prec - type(mld_Tbaseprec_type). + ! The current level preconditioner (aka smoother). + ! iprcparm - integer, dimension(:), allocatable. + ! The integer parameters defining the multilevel strategy. + ! rprcparm - real(psb_Ypk_), dimension(:), allocatable. + ! The real parameters defining the multilevel strategy. + ! ac - The local part of the current-level matrix, built by + ! coarsening the previous-level matrix. + ! desc_ac - type(psb_desc_type). + ! The communication descriptor associated to the matrix + ! stored in ac. + ! base_a - type(psb_z_sparse_mat), pointer. + ! Pointer (really a pointer!) to the local part of the current + ! matrix (so we have a unified treatment of residuals). + ! We need this to avoid passing explicitly the current matrix + ! to the routine which applies the preconditioner. + ! base_desc - type(psb_desc_type), pointer. + ! Pointer to the communication descriptor associated to the + ! matrix pointed by base_a. + ! map - Stores the maps (restriction and prolongation) between the + ! vector spaces associated to the index spaces of the previous + ! and current levels. + ! + ! + ! Type: mld_Tbaseprec_type. + ! + ! It holds the smoother (base preconditioner) at a single level. + ! + ! type mld_Tbaseprec_type + ! type(psb_T_sparse_mat), allocatable :: av(:) + ! IntrType(psb_Tpk_), allocatable :: d(:) + ! type(psb_desc_type) :: desc_data + ! integer, allocatable :: iprcparm(:) + ! real(psb_Tpk_), allocatable :: rprcparm(:) + ! integer, allocatable :: perm(:), invperm(:) + ! end type mld_sbaseprec_type + ! + ! Note that IntrType denotes the real or complex data type, and psb_Tpk denotes + ! the kind of the real or complex type, according to the real/complex, single/double + ! precision version of MLD2P4. + ! + ! av - type(psb_T_sparse_mat), dimension(:), allocatable(:). + ! The sparse matrices needed to apply the preconditioner at + ! the current level ilev. + ! av(mld_l_pr_) - The L factor of the ILU factorization of the local + ! diagonal block of the current-level matrix A(ilev). + ! av(mld_u_pr_) - The U factor of the ILU factorization of the local + ! diagonal block of A(ilev), except its diagonal entries + ! (stored in d). + ! av(mld_ap_nd_) - The entries of the local part of A(ilev) outside + ! the diagonal block, for block-Jacobi sweeps. + ! d - real/complex(psb_Tpk_), dimension(:), allocatable. + ! The diagonal entries of the U factor in the ILU factorization + ! of A(ilev). + ! desc_data - type(psb_desc_type). + ! The communication descriptor associated to the base preconditioner, + ! i.e. to the sparse matrices needed to apply the base preconditioner + ! at the current level. + ! iprcparm - integer, dimension(:), allocatable. + ! The integer parameters defining the base preconditioner K(ilev) + ! (the iprcparm entries and values are specified below). + ! rprcparm - real(psb_Tpk_), dimension(:), allocatable. + ! The real parameters defining the base preconditioner K(ilev) + ! (the rprcparm entries and values are specified below). + ! perm - integer, dimension(:), allocatable. + ! The row and column permutations applied to the local part of + ! A(ilev) (defined only if iprcparm(mld_sub_ren_)>0). + ! invperm - integer, dimension(:), allocatable. + ! The inverse of the permutation stored in perm. + ! + ! Note that when the LU factorization of the (local part of the) matrix A(ilev) is + ! computed instead of the ILU one, by using UMFPACK, SuperLU or SuperLU_dist, the + ! corresponding L and U factors are stored in data structures provided by those + ! packages and pointed by prec%iprcparm(mld_umf_ptr), prec%iprcparm(mld_slu_ptr) + ! or prec%iprcparm(mld_slud_ptr). + ! + + type mld_sbaseprec_type + type(psb_s_sparse_mat), allocatable :: av(:) + real(psb_spk_), allocatable :: d(:) + type(psb_desc_type) :: desc_data + integer, allocatable :: iprcparm(:) + real(psb_spk_), allocatable :: rprcparm(:) + integer, allocatable :: perm(:), invperm(:) + end type mld_sbaseprec_type + + type mld_sonelev_type + type(mld_sbaseprec_type) :: prec + integer, allocatable :: iprcparm(:) + real(psb_spk_), allocatable :: rprcparm(:) + type(psb_s_sparse_mat) :: ac + type(psb_desc_type) :: desc_ac + type(psb_s_sparse_mat), pointer :: base_a => null() + type(psb_desc_type), pointer :: base_desc => null() + type(psb_slinmap_type) :: map + end type mld_sonelev_type + + + type, extends(psb_sprec_type) :: mld_sprec_type + type(mld_sonelev_type), allocatable :: precv(:) + contains + procedure, pass(prec) :: s_apply2v => mld_s_apply2v + procedure, pass(prec) :: s_apply1v => mld_s_apply1v + end type mld_sprec_type + + + ! + ! Interfaces to routines for checking the definition of the preconditioner, + ! for printing its description and for deallocating its data structure + ! + + interface mld_precfree + module procedure mld_sbase_precfree, mld_s_onelev_precfree, mld_sprec_free + end interface + + interface mld_nullify_baseprec + module procedure mld_nullify_sbaseprec + end interface + + interface mld_nullify_onelevprec + module procedure mld_nullify_s_onelevprec + end interface + + + interface mld_precdescr + module procedure mld_sfile_prec_descr + end interface + + + interface mld_sizeof + module procedure mld_sprec_sizeof, mld_sbaseprec_sizeof, mld_s_onelev_prec_sizeof + end interface + + + interface mld_precaply + subroutine mld_sprecaply(prec,x,y,desc_data,info,trans,work) + use psb_base_mod, only : psb_s_sparse_mat, psb_desc_type, psb_spk_ + import mld_sprec_type + type(psb_desc_type),intent(in) :: desc_data + type(mld_sprec_type), intent(in) :: prec + real(psb_spk_),intent(in) :: x(:) + real(psb_spk_),intent(inout) :: y(:) + integer, intent(out) :: info + character(len=1), optional :: trans + real(psb_spk_),intent(inout), optional, target :: work(:) + end subroutine mld_sprecaply + subroutine mld_sprecaply1(prec,x,desc_data,info,trans) + use psb_base_mod, only : psb_s_sparse_mat, psb_desc_type, psb_spk_ + import mld_sprec_type + type(psb_desc_type),intent(in) :: desc_data + type(mld_sprec_type), intent(in) :: prec + real(psb_spk_),intent(inout) :: x(:) + integer, intent(out) :: info + character(len=1), optional :: trans + end subroutine mld_sprecaply1 + end interface + +contains + ! + ! Function returning the size of the mld_prec_type data structure + ! + + function mld_sprec_sizeof(prec) result(val) + implicit none + type(mld_sprec_type), intent(in) :: prec + integer(psb_long_int_k_) :: val + integer :: i + val = 0 + if (allocated(prec%precv)) then + do i=1, size(prec%precv) + val = val + mld_sizeof(prec%precv(i)) + end do + end if + end function mld_sprec_sizeof + + function mld_sbaseprec_sizeof(prec) result(val) + implicit none + type(mld_sbaseprec_type), intent(in) :: prec + integer(psb_long_int_k_) :: val + integer :: i + + val = 0 + if (allocated(prec%iprcparm)) then + val = val + psb_sizeof_int * size(prec%iprcparm) + if (prec%iprcparm(mld_prec_status_) == mld_prec_built_) then + select case(prec%iprcparm(mld_sub_solve_)) + case(mld_ilu_n_,mld_ilu_t_) + ! do nothing + case(mld_slu_) + case(mld_umf_) + case(mld_sludist_) + case default + end select + + end if + end if + if (allocated(prec%rprcparm)) val = val + psb_sizeof_sp * size(prec%rprcparm) + if (allocated(prec%d)) val = val + psb_sizeof_sp * size(prec%d) + if (allocated(prec%perm)) val = val + psb_sizeof_int * size(prec%perm) + if (allocated(prec%invperm)) val = val + psb_sizeof_int * size(prec%invperm) + val = val + psb_sizeof(prec%desc_data) + if (allocated(prec%av)) then + do i=1,size(prec%av) + val = val + psb_sizeof(prec%av(i)) + end do + end if + + end function mld_sbaseprec_sizeof + + function mld_s_onelev_prec_sizeof(prec) result(val) + implicit none + type(mld_sonelev_type), intent(in) :: prec + integer(psb_long_int_k_) :: val + integer :: i + + val = mld_sizeof(prec%prec) + if (allocated(prec%iprcparm)) then + val = val + psb_sizeof_int * size(prec%iprcparm) + end if + if (allocated(prec%rprcparm)) val = val + psb_sizeof_sp * size(prec%rprcparm) + val = val + psb_sizeof(prec%desc_ac) + val = val + psb_sizeof(prec%ac) + val = val + psb_sizeof(prec%map) + + end function mld_s_onelev_prec_sizeof + + ! + ! Subroutine: mld_file_prec_descr + ! Version: real + ! + ! This routine prints a description of the preconditioner to the standard + ! output or to a file. It must be called after the preconditioner has been + ! built by mld_precbld. + ! + ! Arguments: + ! p - type(mld_Tprec_type), input. + ! The preconditioner data structure to be printed out. + ! info - integer, output. + ! error code. + ! iout - integer, input, optional. + ! The id of the file where the preconditioner description + ! will be printed. If iout is not present, then the standard + ! output is condidered. + ! + subroutine mld_sfile_prec_descr(p,info,iout) + implicit none + + ! Arguments + type(mld_sprec_type), intent(in) :: p + integer, intent(out) :: info + integer, intent(in), optional :: iout + + ! Local variables + integer :: ilev, nlev + integer :: ictxt, me, np + character(len=20), parameter :: name='mld_file_prec_descr' + integer :: iout_ + + info = 0 + if (present(iout)) then + iout_ = iout + else + iout_ = 6 + end if + if (iout_ < 0) iout_ = 6 + + if (allocated(p%precv)) then + ictxt = psb_cd_get_context(p%precv(1)%prec%desc_data) + + call psb_info(ictxt,me,np) + + ! + ! The preconditioner description is printed by processor psb_root_. + ! This agrees with the fact that all the parameters defining the + ! preconditioner have the same values on all the procs (this is + ! ensured by mld_precbld). + ! + if (me==psb_root_) then + + write(iout_,*) + write(iout_,*) 'Preconditioner description' + nlev = size(p%precv) + if (nlev >= 1) then + ! + ! Print description of base preconditioner + ! + + write(iout_,*) ' ' + + if (nlev > 1) then + write(iout_,*) 'Multilevel Schwarz' + write(iout_,*) + write(iout_,*) 'Base preconditioner (smoother) details' + endif + + ilev = 1 + call mld_base_prec_descr(iout_,p%precv(ilev)%prec%iprcparm,info,& + & rprcparm=p%precv(ilev)%prec%rprcparm) + + end if + + if (nlev > 1) then + + ! + ! Print multilevel details + ! + write(iout_,*) + write(iout_,*) 'Multilevel details' + + do ilev = 2, nlev + if (.not.allocated(p%precv(ilev)%iprcparm)) then + info = 3111 + write(iout_,*) ' ',name,': error: inconsistent MLPREC part, should call MLD_PRECINIT' + return + endif + end do + + write(iout_,*) ' Number of levels: ',nlev + + ! + ! Currently, all the preconditioner parameters must have the same value at levels + ! 2,...,nlev-1, hence only the values at level 2 are printed + ! + + ilev=2 + call mld_ml_alg_descr(iout_,ilev,p%precv(ilev)%iprcparm, info,& + & rprcparm=p%precv(ilev)%rprcparm) + + ! + ! Coarse matrices are different at levels 2,...,nlev-1, hence related + ! info is printed separately + ! + write(iout_,*) + do ilev = 2, nlev-1 + call mld_ml_level_descr(iout_,ilev,p%precv(ilev)%iprcparm,& + & p%precv(ilev)%map%naggr,info,& + & rprcparm=p%precv(ilev)%rprcparm) + end do + + ! + ! Print coarsest level details + ! + + ilev = nlev + write(iout_,*) + call mld_ml_coarse_descr(iout_,ilev,& + & p%precv(ilev)%iprcparm,p%precv(ilev)%prec%iprcparm,& + & p%precv(ilev)%map%naggr,info,& + & rprcparm=p%precv(ilev)%rprcparm, & + & rprcparm2=p%precv(ilev)%prec%rprcparm) + + end if + + endif + write(iout_,*) + else + + write(iout_,*) trim(name), & + & ': Error: no base preconditioner available, something is wrong!' + info = -2 + return + endif + + + end subroutine mld_sfile_prec_descr + + + ! + ! Subroutines: mld_Tbase_precfree, mld_T_onelev_precfree, mld_Tprec_free + ! Version: real/complex + ! + ! These routines deallocate the mld_Tbaseprec_type, mld_Tonelev_type and + ! mld_Tprec_type data structures. + ! + ! Arguments: + ! p - type(mld_Tbaseprec_type/mld_Tonelev_type/mld_Tprec_type), input. + ! The data structure to be deallocated. + ! info - integer, output. + ! error code. + ! + subroutine mld_sbase_precfree(p,info) + implicit none + + type(mld_sbaseprec_type), intent(inout) :: p + integer, intent(out) :: info + integer :: i + + info = 0 + + ! Actually we might just deallocate the top level array, except + ! for the inner UMFPACK or SLU stuff + + if (allocated(p%d)) then + deallocate(p%d,stat=info) + end if + + if (allocated(p%av)) then + do i=1,size(p%av) + call psb_sp_free(p%av(i),info) + if (info /= 0) then + ! Actually, we don't care here about this. + ! Just let it go. + ! return + end if + enddo + deallocate(p%av,stat=info) + end if + + if (allocated(p%desc_data%matrix_data)) & + & call psb_cdfree(p%desc_data,info) + + if (allocated(p%rprcparm)) then + deallocate(p%rprcparm,stat=info) + end if + + + if (allocated(p%perm)) then + deallocate(p%perm,stat=info) + endif + + if (allocated(p%invperm)) then + deallocate(p%invperm,stat=info) + endif + + if (allocated(p%iprcparm)) then + if (p%iprcparm(mld_prec_status_) == mld_prec_built_) then + if (p%iprcparm(mld_sub_solve_)==mld_slu_) then + call mld_sslu_free(p%iprcparm(mld_slu_ptr_),info) + end if + end if + deallocate(p%iprcparm,stat=info) + end if + call mld_nullify_baseprec(p) + end subroutine mld_sbase_precfree + + + subroutine mld_s_onelev_precfree(p,info) + implicit none + + type(mld_sonelev_type), intent(inout) :: p + integer, intent(out) :: info + integer :: i + + info = 0 + + ! Actually we might just deallocate the top level array, except + ! for the inner UMFPACK or SLU stuff + call mld_precfree(p%prec,info) + + call psb_sp_free(p%ac,info) + if (allocated(p%desc_ac%matrix_data)) & + & call psb_cdfree(p%desc_ac,info) + + if (allocated(p%rprcparm)) then + deallocate(p%rprcparm,stat=info) + end if + ! This is a pointer to something else, must not free it here. + nullify(p%base_a) + ! This is a pointer to something else, must not free it here. + nullify(p%base_desc) + + ! + ! free explicitly map??? + ! For now thanks to allocatable semantics + ! works anyway. + ! + + call mld_nullify_onelevprec(p) + end subroutine mld_s_onelev_precfree + + + subroutine mld_nullify_s_onelevprec(p) + implicit none + + type(mld_sonelev_type), intent(inout) :: p + + nullify(p%base_a) + nullify(p%base_desc) + + end subroutine mld_nullify_s_onelevprec + + subroutine mld_nullify_sbaseprec(p) + implicit none + + type(mld_sbaseprec_type), intent(inout) :: p + + + end subroutine mld_nullify_sbaseprec + + + subroutine mld_sprec_free(p,info) + + use psb_base_mod + + implicit none + + ! Arguments + type(mld_sprec_type), intent(inout) :: p + integer, intent(out) :: info + + ! Local variables + integer :: me,err_act,i + character(len=20) :: name + + if(psb_get_errstatus().ne.0) return + info=0 + name = 'mld_dprecfree' + call psb_erractionsave(err_act) + + me=-1 + + if (allocated(p%precv)) then + do i=1,size(p%precv) + call mld_precfree(p%precv(i),info) + end do + deallocate(p%precv) + end if + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act.eq.psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine mld_sprec_free + + + + subroutine mld_s_apply2v(prec,x,y,desc_data,info,trans,work) + use psb_base_mod + type(psb_desc_type),intent(in) :: desc_data + class(mld_sprec_type), intent(in) :: prec + real(psb_spk_),intent(in) :: x(:) + real(psb_spk_),intent(inout) :: y(:) + integer, intent(out) :: info + character(len=1), optional :: trans + real(psb_spk_),intent(inout), optional, target :: work(:) + Integer :: err_act + character(len=20) :: name='s_prec_apply' + + call psb_erractionsave(err_act) + + select type(prec) + type is (mld_sprec_type) + call mld_precaply(prec,x,y,desc_data,info,trans,work) + class default + info = 700 + call psb_errpush(info,name) + goto 9999 + end select + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine mld_s_apply2v + subroutine mld_s_apply1v(prec,x,desc_data,info,trans) + use psb_base_mod + type(psb_desc_type),intent(in) :: desc_data + class(mld_sprec_type), intent(in) :: prec + real(psb_spk_),intent(inout) :: x(:) + integer, intent(out) :: info + character(len=1), optional :: trans + Integer :: err_act + character(len=20) :: name='s_prec_apply' + + call psb_erractionsave(err_act) + + select type(prec) + type is (mld_sprec_type) + call mld_precaply(prec,x,desc_data,info,trans) + class default + info = 700 + call psb_errpush(info,name) + goto 9999 + end select + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine mld_s_apply1v + +end module mld_s_prec_type diff --git a/mlprec/mld_z_prec_type.f90 b/mlprec/mld_z_prec_type.f90 new file mode 100644 index 00000000..81302bff --- /dev/null +++ b/mlprec/mld_z_prec_type.f90 @@ -0,0 +1,685 @@ +!!$ +!!$ +!!$ MLD2P4 version 1.1 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 2.3.1) +!!$ +!!$ (C) Copyright 2008,2009 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +! File: mld_prec_type.f90 +! +! Module: mld_prec_type +! +! This module defines: +! - the mld_prec_type data structure containing the preconditioner and related +! data structures; +! - integer constants defining the preconditioner; +! - character constants describing the preconditioner (used by the routines +! printing out a preconditioner description); +! - the interfaces to the routines for the management of the preconditioner +! data structure (see below). +! +! It contains routines for +! - converting character constants defining the preconditioner into integer +! constants; +! - checking if the preconditioner is correctly defined; +! - printing a description of the preconditioner; +! - deallocating the preconditioner data structure. +! + +module mld_z_prec_type + + use mld_base_prec_type + ! + ! Type: mld_Tprec_type. + ! + ! It is the data type containing all the information about the multilevel + ! preconditioner (here and in the following 'T' denotes 'd', 's', 'c' and + ! 'z', according to the real/complex, single/double precision version of + ! MLD2P4). It consists of an array of 'one-level' intermediate data structures + ! of type mld_Tonelev_type, each containing the information needed to apply + ! the smoothing and the coarse-space correction at a generic level. + ! + ! type mld_Tprec_type + ! type(mld_Tonelev_type), allocatable :: precv(:) + ! end type mld_Tprec_type + ! + ! Note that the levels are numbered in increasing order starting from + ! the finest one and the number of levels is given by size(precv(:)). + ! + ! + ! Type: mld_Tonelev_type. + ! + ! It is the data type containing the necessary items for the current + ! level (essentially, the base preconditioner, the current-level matrix + ! and the restriction and prolongation operators). + ! + ! type mld_Tonelev_type + ! type(mld_Tbaseprec_type) :: prec + ! integer, allocatable :: iprcparm(:) + ! real(psb_Tpk_), allocatable :: rprcparm(:) + ! type(psb_T_sparse_mat) :: ac + ! type(psb_desc_type) :: desc_ac + ! type(psb_T_sparse_mat), pointer :: base_a => null() + ! type(psb_desc_type), pointer :: base_desc => null() + ! type(psb_Tlinmap_type) :: map + ! end type mld_Tonelev_type + ! + ! Note that psb_Tpk denotes the kind of the real data type to be chosen + ! according to single/double precision version of MLD2P4. + ! + ! prec - type(mld_Tbaseprec_type). + ! The current level preconditioner (aka smoother). + ! iprcparm - integer, dimension(:), allocatable. + ! The integer parameters defining the multilevel strategy. + ! rprcparm - real(psb_Ypk_), dimension(:), allocatable. + ! The real parameters defining the multilevel strategy. + ! ac - The local part of the current-level matrix, built by + ! coarsening the previous-level matrix. + ! desc_ac - type(psb_desc_type). + ! The communication descriptor associated to the matrix + ! stored in ac. + ! base_a - type(psb_z_sparse_mat), pointer. + ! Pointer (really a pointer!) to the local part of the current + ! matrix (so we have a unified treatment of residuals). + ! We need this to avoid passing explicitly the current matrix + ! to the routine which applies the preconditioner. + ! base_desc - type(psb_desc_type), pointer. + ! Pointer to the communication descriptor associated to the + ! matrix pointed by base_a. + ! map - Stores the maps (restriction and prolongation) between the + ! vector spaces associated to the index spaces of the previous + ! and current levels. + ! + ! + ! Type: mld_Tbaseprec_type. + ! + ! It holds the smoother (base preconditioner) at a single level. + ! + ! type mld_Tbaseprec_type + ! type(psb_T_sparse_mat), allocatable :: av(:) + ! IntrType(psb_Tpk_), allocatable :: d(:) + ! type(psb_desc_type) :: desc_data + ! integer, allocatable :: iprcparm(:) + ! real(psb_Tpk_), allocatable :: rprcparm(:) + ! integer, allocatable :: perm(:), invperm(:) + ! end type mld_sbaseprec_type + ! + ! Note that IntrType denotes the real or complex data type, and psb_Tpk denotes + ! the kind of the real or complex type, according to the real/complex, single/double + ! precision version of MLD2P4. + ! + ! av - type(psb_T_sparse_mat), dimension(:), allocatable(:). + ! The sparse matrices needed to apply the preconditioner at + ! the current level ilev. + ! av(mld_l_pr_) - The L factor of the ILU factorization of the local + ! diagonal block of the current-level matrix A(ilev). + ! av(mld_u_pr_) - The U factor of the ILU factorization of the local + ! diagonal block of A(ilev), except its diagonal entries + ! (stored in d). + ! av(mld_ap_nd_) - The entries of the local part of A(ilev) outside + ! the diagonal block, for block-Jacobi sweeps. + ! d - real/complex(psb_Tpk_), dimension(:), allocatable. + ! The diagonal entries of the U factor in the ILU factorization + ! of A(ilev). + ! desc_data - type(psb_desc_type). + ! The communication descriptor associated to the base preconditioner, + ! i.e. to the sparse matrices needed to apply the base preconditioner + ! at the current level. + ! iprcparm - integer, dimension(:), allocatable. + ! The integer parameters defining the base preconditioner K(ilev) + ! (the iprcparm entries and values are specified below). + ! rprcparm - real(psb_Tpk_), dimension(:), allocatable. + ! The real parameters defining the base preconditioner K(ilev) + ! (the rprcparm entries and values are specified below). + ! perm - integer, dimension(:), allocatable. + ! The row and column permutations applied to the local part of + ! A(ilev) (defined only if iprcparm(mld_sub_ren_)>0). + ! invperm - integer, dimension(:), allocatable. + ! The inverse of the permutation stored in perm. + ! + ! Note that when the LU factorization of the (local part of the) matrix A(ilev) is + ! computed instead of the ILU one, by using UMFPACK, SuperLU or SuperLU_dist, the + ! corresponding L and U factors are stored in data structures provided by those + ! packages and pointed by prec%iprcparm(mld_umf_ptr), prec%iprcparm(mld_slu_ptr) + ! or prec%iprcparm(mld_slud_ptr). + ! + + type mld_zbaseprec_type + type(psb_z_sparse_mat), allocatable :: av(:) + complex(psb_dpk_), allocatable :: d(:) + type(psb_desc_type) :: desc_data + integer, allocatable :: iprcparm(:) + real(psb_dpk_), allocatable :: rprcparm(:) + integer, allocatable :: perm(:), invperm(:) + end type mld_zbaseprec_type + + type mld_zonelev_type + type(mld_zbaseprec_type) :: prec + integer, allocatable :: iprcparm(:) + real(psb_dpk_), allocatable :: rprcparm(:) + type(psb_z_sparse_mat) :: ac + type(psb_desc_type) :: desc_ac + type(psb_z_sparse_mat), pointer :: base_a => null() + type(psb_desc_type), pointer :: base_desc => null() + type(psb_zlinmap_type) :: map + end type mld_zonelev_type + + type, extends(psb_zprec_type) :: mld_zprec_type + type(mld_zonelev_type), allocatable :: precv(:) + contains + procedure, pass(prec) :: z_apply2v => mld_z_apply2v + procedure, pass(prec) :: z_apply1v => mld_z_apply1v + end type mld_zprec_type + + + ! + ! Interfaces to routines for checking the definition of the preconditioner, + ! for printing its description and for deallocating its data structure + ! + + interface mld_precfree + module procedure mld_zbase_precfree, mld_z_onelev_precfree, mld_zprec_free + end interface + + interface mld_nullify_baseprec + module procedure mld_nullify_zbaseprec + end interface + + interface mld_nullify_onelevprec + module procedure mld_nullify_z_onelevprec + end interface + + interface mld_precdescr + module procedure mld_zfile_prec_descr + end interface + + interface mld_sizeof + module procedure mld_zprec_sizeof, mld_zbaseprec_sizeof, mld_z_onelev_prec_sizeof + end interface + + interface mld_precaply + subroutine mld_zprecaply(prec,x,y,desc_data,info,trans,work) + use psb_base_mod, only : psb_z_sparse_mat, psb_desc_type, psb_dpk_ + import mld_zprec_type + type(psb_desc_type),intent(in) :: desc_data + type(mld_zprec_type), intent(in) :: prec + complex(psb_dpk_),intent(in) :: x(:) + complex(psb_dpk_),intent(inout) :: y(:) + integer, intent(out) :: info + character(len=1), optional :: trans + complex(psb_dpk_),intent(inout), optional, target :: work(:) + end subroutine mld_zprecaply + subroutine mld_zprecaply1(prec,x,desc_data,info,trans) + use psb_base_mod, only : psb_z_sparse_mat, psb_desc_type, psb_dpk_ + import mld_zprec_type + type(psb_desc_type),intent(in) :: desc_data + type(mld_zprec_type), intent(in) :: prec + complex(psb_dpk_),intent(inout) :: x(:) + integer, intent(out) :: info + character(len=1), optional :: trans + end subroutine mld_zprecaply1 + end interface + + +contains + ! + ! Function returning the size of the mld_prec_type data structure + ! + function mld_zprec_sizeof(prec) result(val) + implicit none + type(mld_zprec_type), intent(in) :: prec + integer(psb_long_int_k_) :: val + integer :: i + val = 0 + if (allocated(prec%precv)) then + do i=1, size(prec%precv) + val = val + mld_sizeof(prec%precv(i)) + end do + end if + end function mld_zprec_sizeof + + function mld_zbaseprec_sizeof(prec) result(val) + implicit none + type(mld_zbaseprec_type), intent(in) :: prec + integer(psb_long_int_k_) :: val + integer :: i + + val = 0 + if (allocated(prec%iprcparm)) then + val = val + psb_sizeof_int * size(prec%iprcparm) + if (prec%iprcparm(mld_prec_status_) == mld_prec_built_) then + select case(prec%iprcparm(mld_sub_solve_)) + case(mld_ilu_n_,mld_ilu_t_) + ! do nothing + case(mld_slu_) + case(mld_umf_) + case(mld_sludist_) + case default + end select + + end if + end if + if (allocated(prec%rprcparm)) val = val + psb_sizeof_dp * size(prec%rprcparm) + if (allocated(prec%d)) val = val + 2 * psb_sizeof_dp * size(prec%d) + if (allocated(prec%perm)) val = val + psb_sizeof_int * size(prec%perm) + if (allocated(prec%invperm)) val = val + psb_sizeof_int * size(prec%invperm) + val = val + psb_sizeof(prec%desc_data) + if (allocated(prec%av)) then + do i=1,size(prec%av) + val = val + psb_sizeof(prec%av(i)) + end do + end if + + end function mld_zbaseprec_sizeof + + function mld_z_onelev_prec_sizeof(prec) result(val) + implicit none + type(mld_zonelev_type), intent(in) :: prec + integer(psb_long_int_k_) :: val + integer :: i + + val = mld_sizeof(prec%prec) + if (allocated(prec%iprcparm)) then + val = val + psb_sizeof_int * size(prec%iprcparm) + end if + if (allocated(prec%rprcparm)) val = val + psb_sizeof_dp * size(prec%rprcparm) + val = val + psb_sizeof(prec%desc_ac) + val = val + psb_sizeof(prec%ac) + val = val + psb_sizeof(prec%map) + + end function mld_z_onelev_prec_sizeof + + ! + ! Subroutine: mld_zfile_prec_descr + ! Version: complex + ! + ! This routine prints to a file a description of the preconditioner. + ! + ! Arguments: + ! p - type(mld_Tprec_type), input. + ! The preconditioner data structure to be printed out. + ! iout - integer, input. + ! The id of the file where the preconditioner description + ! will be printed. + ! + subroutine mld_zfile_prec_descr(p,info,iout) + implicit none + + ! Arguments + type(mld_zprec_type), intent(in) :: p + integer, intent(out) :: info + integer, intent(in), optional :: iout + + ! Local variables + integer :: ilev, nlev + integer :: ictxt, me, np + character(len=20), parameter :: name='mld_file_prec_descr' + integer :: iout_ + + info = 0 + if (present(iout)) then + iout_ = iout + else + iout_ = 6 + end if + if (iout_ < 0) iout_ = 6 + + if (allocated(p%precv)) then + ictxt = psb_cd_get_context(p%precv(1)%prec%desc_data) + + call psb_info(ictxt,me,np) + + ! + ! The preconditioner description is printed by processor psb_root_. + ! This agrees with the fact that all the parameters defining the + ! preconditioner have the same values on all the procs (this is + ! ensured by mld_precbld). + ! + if (me==psb_root_) then + + write(iout_,*) + write(iout_,*) 'Preconditioner description' + nlev = size(p%precv) + if (nlev >= 1) then + ! + ! Print description of base preconditioner + ! + + write(iout_,*) ' ' + + if (nlev > 1) then + write(iout_,*) 'Multilevel Schwarz' + write(iout_,*) + write(iout_,*) 'Base preconditioner (smoother) details' + endif + + ilev = 1 + call mld_base_prec_descr(iout_,p%precv(ilev)%prec%iprcparm,info,& + & dprcparm=p%precv(ilev)%prec%rprcparm) + + end if + + if (nlev > 1) then + + ! + ! Print multilevel details + ! + write(iout_,*) + write(iout_,*) 'Multilevel details' + + do ilev = 2, nlev + if (.not.allocated(p%precv(ilev)%iprcparm)) then + info = 3111 + write(iout_,*) ' ',name,': error: inconsistent MLPREC part, should call MLD_PRECINIT' + return + endif + end do + + write(iout_,*) ' Number of levels: ',nlev + + ! + ! Currently, all the preconditioner parameters must have the same value at levels + ! 2,...,nlev-1, hence only the values at level 2 are printed + ! + + ilev=2 + call mld_ml_alg_descr(iout_,ilev,p%precv(ilev)%iprcparm, info,& + & dprcparm=p%precv(ilev)%rprcparm) + + ! + ! Coarse matrices are different at levels 2,...,nlev-1, hence related + ! info is printed separately + ! + write(iout_,*) + do ilev = 2, nlev-1 + call mld_ml_level_descr(iout_,ilev,p%precv(ilev)%iprcparm,& + & p%precv(ilev)%map%naggr,info,& + & dprcparm=p%precv(ilev)%rprcparm) + end do + + ! + ! Print coarsest level details + ! + + ilev = nlev + write(iout_,*) + call mld_ml_coarse_descr(iout_,ilev,& + & p%precv(ilev)%iprcparm,p%precv(ilev)%prec%iprcparm,& + & p%precv(ilev)%map%naggr,info,& + & dprcparm=p%precv(ilev)%rprcparm,& + & dprcparm2=p%precv(ilev)%prec%rprcparm) + end if + + endif + write(iout_,*) + else + write(iout_,*) trim(name), & + & ': Error: no base preconditioner available, something is wrong!' + info = -2 + return + endif + + end subroutine mld_zfile_prec_descr + + ! + ! Subroutines: mld_Tbase_precfree, mld_T_onelev_precfree, mld_Tprec_free + ! Version: real/complex + ! + ! These routines deallocate the mld_Tbaseprec_type, mld_Tonelev_type and + ! mld_Tprec_type data structures. + ! + ! Arguments: + ! p - type(mld_Tbaseprec_type/mld_Tonelev_type/mld_Tprec_type), input. + ! The data structure to be deallocated. + ! info - integer, output. + ! error code. + ! + subroutine mld_zbase_precfree(p,info) + implicit none + type(mld_zbaseprec_type), intent(inout) :: p + integer, intent(out) :: info + integer :: i + + info = 0 + + if (allocated(p%d)) then + deallocate(p%d,stat=info) + end if + + if (allocated(p%av)) then + do i=1,size(p%av) + call psb_sp_free(p%av(i),info) + if (info /= 0) then + ! Actually, we don't care here about this. + ! Just let it go. + ! return + end if + enddo + deallocate(p%av,stat=info) + + end if + if (allocated(p%desc_data%matrix_data)) & + & call psb_cdfree(p%desc_data,info) + + if (allocated(p%rprcparm)) then + deallocate(p%rprcparm,stat=info) + end if + + if (allocated(p%perm)) then + deallocate(p%perm,stat=info) + endif + + if (allocated(p%invperm)) then + deallocate(p%invperm,stat=info) + endif + + if (allocated(p%iprcparm)) then + if (p%iprcparm(mld_prec_status_) == mld_prec_built_) then + if (p%iprcparm(mld_sub_solve_)==mld_slu_) then + call mld_zslu_free(p%iprcparm(mld_slu_ptr_),info) + end if + if (p%iprcparm(mld_sub_solve_)==mld_umf_) then + call mld_zumf_free(p%iprcparm(mld_umf_symptr_),& + & p%iprcparm(mld_umf_numptr_),info) + end if + end if + deallocate(p%iprcparm,stat=info) + end if + call mld_nullify_baseprec(p) + end subroutine mld_zbase_precfree + + subroutine mld_z_onelev_precfree(p,info) + implicit none + + type(mld_zonelev_type), intent(inout) :: p + integer, intent(out) :: info + integer :: i + + info = 0 + + ! Actually we might just deallocate the top level array, except + ! for the inner UMFPACK or SLU stuff + call mld_precfree(p%prec,info) + + call psb_sp_free(p%ac,info) + if (allocated(p%desc_ac%matrix_data)) & + & call psb_cdfree(p%desc_ac,info) + + if (allocated(p%rprcparm)) then + deallocate(p%rprcparm,stat=info) + end if + ! This is a pointer to something else, must not free it here. + nullify(p%base_a) + ! This is a pointer to something else, must not free it here. + nullify(p%base_desc) + + ! + ! free explicitly map??? + ! For now thanks to allocatable semantics + ! works anyway. + ! + + call mld_nullify_onelevprec(p) + end subroutine mld_z_onelev_precfree + + subroutine mld_nullify_z_onelevprec(p) + implicit none + + type(mld_zonelev_type), intent(inout) :: p + + nullify(p%base_a) + nullify(p%base_desc) + + end subroutine mld_nullify_z_onelevprec + + + subroutine mld_nullify_zbaseprec(p) + implicit none + + type(mld_zbaseprec_type), intent(inout) :: p + + + end subroutine mld_nullify_zbaseprec + + subroutine mld_zprec_free(p,info) + + use psb_base_mod + + implicit none + + ! Arguments + type(mld_zprec_type), intent(inout) :: p + integer, intent(out) :: info + + ! Local variables + integer :: me,err_act,i + character(len=20) :: name + + if(psb_get_errstatus().ne.0) return + info=0 + name = 'mld_dprecfree' + call psb_erractionsave(err_act) + + me=-1 + + if (allocated(p%precv)) then + do i=1,size(p%precv) + call mld_precfree(p%precv(i),info) + end do + deallocate(p%precv) + end if + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act.eq.psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine mld_zprec_free + + subroutine mld_z_apply2v(prec,x,y,desc_data,info,trans,work) + use psb_base_mod + type(psb_desc_type),intent(in) :: desc_data + class(mld_zprec_type), intent(in) :: prec + complex(psb_dpk_),intent(in) :: x(:) + complex(psb_dpk_),intent(inout) :: y(:) + integer, intent(out) :: info + character(len=1), optional :: trans + complex(psb_dpk_),intent(inout), optional, target :: work(:) + Integer :: err_act + character(len=20) :: name='z_prec_apply' + + call psb_erractionsave(err_act) + + select type(prec) +!!$ type is (mld_zprec_type) +!!$ call mld_precaply(prec,x,y,desc_data,info,trans,work) + class default + info = 700 + call psb_errpush(info,name) + goto 9999 + end select + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine mld_z_apply2v + subroutine mld_z_apply1v(prec,x,desc_data,info,trans) + use psb_base_mod + + type(psb_desc_type),intent(in) :: desc_data + class(mld_zprec_type), intent(in) :: prec + complex(psb_dpk_),intent(inout) :: x(:) + integer, intent(out) :: info + character(len=1), optional :: trans + Integer :: err_act + character(len=20) :: name='z_prec_apply' + + call psb_erractionsave(err_act) + + select type(prec) +!!$ type is (mld_zprec_type) +!!$ call mld_precaply(prec,x,desc_data,info,trans) + class default + info = 700 + call psb_errpush(info,name) + goto 9999 + end select + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine mld_z_apply1v + +end module mld_z_prec_type