mld2p4-smooth-2side:

mlprec/impl/level/mld_d_base_onelev_csetc.f90
 mlprec/impl/level/mld_d_base_onelev_cseti.f90
 mlprec/impl/level/mld_d_base_onelev_csetr.f90
 mlprec/impl/level/mld_d_base_onelev_setc.f90
 mlprec/impl/level/mld_d_base_onelev_seti.f90
 mlprec/impl/level/mld_d_base_onelev_setr.f90
 mlprec/impl/mld_dcprecset.F90
 mlprec/impl/mld_dprecset.F90
 mlprec/mld_d_prec_mod.f90
 tests/pdegen/ppde3d-gs.f90
 tests/pdegen/runs/ppde.inp

SET now works; next step will be some refactoring.
Note: the symmetrized ML for CG with FW/BW Gauss-Seidel does not seem
to work right now.
stopcriterion
Salvatore Filippone 9 years ago
parent d747bc9aae
commit e08492cdaf

@ -49,7 +49,8 @@ subroutine mld_d_base_onelev_csetc(lv,what,val,info,pos)
character(len=*), intent(in) :: val
integer(psb_ipk_), intent(out) :: info
character(len=*), optional, intent(in) :: pos
integer(psb_ipk_) :: err_act
! Local
integer(psb_ipk_) :: ipos_, err_act
character(len=20) :: name='d_base_onelev_csetc'
integer(psb_ipk_) :: ival
@ -61,9 +62,33 @@ subroutine mld_d_base_onelev_csetc(lv,what,val,info,pos)
if (ival >= 0) then
call lv%set(what,ival,info,pos=pos)
else
if (allocated(lv%sm)) then
call lv%sm%set(what,val,info)
if (present(pos)) then
select case(psb_toupper(trim(pos)))
case('PRE')
ipos_ = mld_pre_smooth_
case('POST')
ipos_ = mld_post_smooth_
case default
ipos_ = mld_pre_smooth_
end select
else
ipos_ = mld_pre_smooth_
end if
select case(ipos_)
case(mld_pre_smooth_)
if (allocated(lv%sm)) then
call lv%sm%set(what,val,info)
end if
case (mld_post_smooth_)
if (allocated(lv%sm2a)) then
call lv%sm2a%set(what,val,info)
end if
case default
! Impossible!!
info = psb_err_internal_error_
end select
end if

@ -49,12 +49,26 @@ subroutine mld_d_base_onelev_cseti(lv,what,val,info,pos)
integer(psb_ipk_), intent(in) :: val
integer(psb_ipk_), intent(out) :: info
character(len=*), optional, intent(in) :: pos
Integer(Psb_ipk_) :: err_act
! Local
integer(psb_ipk_) :: ipos_, err_act
character(len=20) :: name='d_base_onelev_cseti'
call psb_erractionsave(err_act)
info = psb_success_
if (present(pos)) then
select case(psb_toupper(trim(pos)))
case('PRE')
ipos_ = mld_pre_smooth_
case('POST')
ipos_ = mld_post_smooth_
case default
ipos_ = mld_pre_smooth_
end select
else
ipos_ = mld_pre_smooth_
end if
select case (psb_toupper(what))
case ('SMOOTHER_SWEEPS')
@ -99,9 +113,20 @@ subroutine mld_d_base_onelev_cseti(lv,what,val,info,pos)
lv%parms%coarse_solve = val
case default
if (allocated(lv%sm)) then
call lv%sm%set(what,val,info)
end if
select case(ipos_)
case(mld_pre_smooth_)
if (allocated(lv%sm)) then
call lv%sm%set(what,val,info)
end if
case (mld_post_smooth_)
if (allocated(lv%sm2a)) then
call lv%sm2a%set(what,val,info)
end if
case default
! Impossible!!
info = psb_err_internal_error_
end select
end select
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)

@ -49,7 +49,8 @@ subroutine mld_d_base_onelev_csetr(lv,what,val,info,pos)
real(psb_dpk_), intent(in) :: val
integer(psb_ipk_), intent(out) :: info
character(len=*), optional, intent(in) :: pos
integer(psb_ipk_) :: err_act
! Local
integer(psb_ipk_) :: ipos_, err_act
character(len=20) :: name='d_base_onelev_csetr'
call psb_erractionsave(err_act)
@ -69,9 +70,32 @@ subroutine mld_d_base_onelev_csetr(lv,what,val,info,pos)
lv%parms%aggr_scale = val
case default
if (allocated(lv%sm)) then
call lv%sm%set(what,val,info)
if (present(pos)) then
select case(psb_toupper(trim(pos)))
case('PRE')
ipos_ = mld_pre_smooth_
case('POST')
ipos_ = mld_post_smooth_
case default
ipos_ = mld_pre_smooth_
end select
else
ipos_ = mld_pre_smooth_
end if
select case(ipos_)
case(mld_pre_smooth_)
if (allocated(lv%sm)) then
call lv%sm%set(what,val,info)
end if
case (mld_post_smooth_)
if (allocated(lv%sm2a)) then
call lv%sm2a%set(what,val,info)
end if
case default
! Impossible!!
info = psb_err_internal_error_
end select
end select
if (info /= psb_success_) goto 9999

@ -49,7 +49,8 @@ subroutine mld_d_base_onelev_setc(lv,what,val,info,pos)
character(len=*), intent(in) :: val
integer(psb_ipk_), intent(out) :: info
character(len=*), optional, intent(in) :: pos
integer(psb_ipk_) :: err_act
! Local
integer(psb_ipk_) :: ipos_, err_act
character(len=20) :: name='d_base_onelev_setc'
integer(psb_ipk_) :: ival
@ -61,9 +62,32 @@ subroutine mld_d_base_onelev_setc(lv,what,val,info,pos)
if (ival >= 0) then
call lv%set(what,ival,info,pos=pos)
else
if (allocated(lv%sm)) then
call lv%sm%set(what,val,info)
if (present(pos)) then
select case(psb_toupper(trim(pos)))
case('PRE')
ipos_ = mld_pre_smooth_
case('POST')
ipos_ = mld_post_smooth_
case default
ipos_ = mld_pre_smooth_
end select
else
ipos_ = mld_pre_smooth_
end if
select case(ipos_)
case(mld_pre_smooth_)
if (allocated(lv%sm)) then
call lv%sm%set(what,val,info)
end if
case (mld_post_smooth_)
if (allocated(lv%sm2a)) then
call lv%sm2a%set(what,val,info)
end if
case default
! Impossible!!
info = psb_err_internal_error_
end select
end if

@ -49,7 +49,8 @@ subroutine mld_d_base_onelev_seti(lv,what,val,info,pos)
integer(psb_ipk_), intent(in) :: val
integer(psb_ipk_), intent(out) :: info
character(len=*), optional, intent(in) :: pos
Integer(Psb_ipk_) :: err_act
! Local
integer(psb_ipk_) :: ipos_, err_act
character(len=20) :: name='d_base_onelev_seti'
call psb_erractionsave(err_act)
@ -99,9 +100,33 @@ subroutine mld_d_base_onelev_seti(lv,what,val,info,pos)
lv%parms%coarse_solve = val
case default
if (allocated(lv%sm)) then
call lv%sm%set(what,val,info)
if (present(pos)) then
select case(psb_toupper(trim(pos)))
case('PRE')
ipos_ = mld_pre_smooth_
case('POST')
ipos_ = mld_post_smooth_
case default
ipos_ = mld_pre_smooth_
end select
else
ipos_ = mld_pre_smooth_
end if
select case(ipos_)
case(mld_pre_smooth_)
if (allocated(lv%sm)) then
call lv%sm%set(what,val,info)
end if
case (mld_post_smooth_)
if (allocated(lv%sm2a)) then
call lv%sm2a%set(what,val,info)
end if
case default
! Impossible!!
info = psb_err_internal_error_
end select
end select
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)

@ -49,7 +49,8 @@ subroutine mld_d_base_onelev_setr(lv,what,val,info,pos)
real(psb_dpk_), intent(in) :: val
integer(psb_ipk_), intent(out) :: info
character(len=*), optional, intent(in) :: pos
integer(psb_ipk_) :: err_act
! Local
integer(psb_ipk_) :: ipos_, err_act
character(len=20) :: name='d_base_onelev_setr'
call psb_erractionsave(err_act)
@ -69,9 +70,32 @@ subroutine mld_d_base_onelev_setr(lv,what,val,info,pos)
lv%parms%aggr_scale = val
case default
if (allocated(lv%sm)) then
call lv%sm%set(what,val,info)
if (present(pos)) then
select case(psb_toupper(trim(pos)))
case('PRE')
ipos_ = mld_pre_smooth_
case('POST')
ipos_ = mld_post_smooth_
case default
ipos_ = mld_pre_smooth_
end select
else
ipos_ = mld_pre_smooth_
end if
select case(ipos_)
case(mld_pre_smooth_)
if (allocated(lv%sm)) then
call lv%sm%set(what,val,info)
end if
case (mld_post_smooth_)
if (allocated(lv%sm2a)) then
call lv%sm2a%set(what,val,info)
end if
case default
! Impossible!!
info = psb_err_internal_error_
end select
end select
if (info /= psb_success_) goto 9999

@ -108,7 +108,7 @@ subroutine mld_dcprecseti(p,what,val,info,ilev,pos)
integer(psb_ipk_), intent(in) :: val
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: ilev
character(len=*), optional, intent(in) :: pos
character(len=*), optional, intent(in) :: pos
! Local variables
integer(psb_ipk_) :: ilev_, nlev_
@ -367,91 +367,123 @@ subroutine mld_dcprecseti(p,what,val,info,ilev,pos)
contains
subroutine onelev_set_smoother(level,val,info,pos)
type(mld_d_onelev_type), intent(inout) :: level
class(mld_d_onelev_type), intent(inout) :: level
integer(psb_ipk_), intent(in) :: val
integer(psb_ipk_), intent(out) :: info
character(len=*), optional, intent(in) :: pos
character(len=*), optional, intent(in) :: pos
! Local
integer(psb_ipk_) :: ipos_
info = psb_success_
if (present(pos)) then
select case(psb_toupper(trim(pos)))
case('PRE')
ipos_ = mld_pre_smooth_
case('POST')
ipos_ = mld_post_smooth_
case default
ipos_ = mld_pre_smooth_
end select
else
ipos_ = mld_pre_smooth_
end if
select case(ipos_)
case(mld_pre_smooth_)
call inner_set_smoother(level%sm,val,info)
case (mld_post_smooth_)
call inner_set_smoother(level%sm2a,val,info)
case default
! Impossible!!
info = psb_err_internal_error_
end select
end subroutine onelev_set_smoother
subroutine inner_set_smoother(sm,val,info)
class(mld_d_base_smoother_type), allocatable, intent(inout) :: sm
integer(psb_ipk_), intent(in) :: val
integer(psb_ipk_), intent(out) :: info
!
! This here requires a bit more attention.
!
select case (val)
case (mld_noprec_)
if (allocated(level%sm)) then
select type (sm => level%sm)
if (allocated(sm)) then
select type (sms => sm)
type is (mld_d_base_smoother_type)
! do nothing
class default
call level%sm%free(info)
if (info == 0) deallocate(level%sm)
call sm%free(info)
if (info == 0) deallocate(sm)
if (info == 0) allocate(mld_d_base_smoother_type ::&
& level%sm, stat=info)
& sm, stat=info)
if (info == 0) allocate(mld_d_id_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
end select
else
allocate(mld_d_base_smoother_type ::&
& level%sm, stat=info)
& sm, stat=info)
if (info ==0) allocate(mld_d_id_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
endif
case (mld_jac_)
if (allocated(level%sm)) then
select type (sm => level%sm)
if (allocated(sm)) then
select type (sms => sm)
class is (mld_d_jac_smoother_type)
! do nothing
! do nothing
class default
call level%sm%free(info)
if (info == 0) deallocate(level%sm)
call sm%free(info)
if (info == 0) deallocate(sm)
if (info == 0) allocate(mld_d_jac_smoother_type :: &
& level%sm, stat=info)
& sm, stat=info)
if (info == 0) allocate(mld_d_diag_solver_type :: &
& level%sm%sv, stat=info)
& sm%sv, stat=info)
end select
else
allocate(mld_d_jac_smoother_type :: level%sm, stat=info)
allocate(mld_d_jac_smoother_type :: sm, stat=info)
if (info == 0) allocate(mld_d_diag_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
endif
case (mld_bjac_)
if (allocated(level%sm)) then
select type (sm => level%sm)
if (allocated(sm)) then
select type (sms => sm)
class is (mld_d_jac_smoother_type)
! do nothing
! do nothing
class default
call level%sm%free(info)
if (info == 0) deallocate(level%sm)
call sm%free(info)
if (info == 0) deallocate(sm)
if (info == 0) allocate(mld_d_jac_smoother_type ::&
& level%sm, stat=info)
& sm, stat=info)
if (info == 0) allocate(mld_d_ilu_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
end select
else
allocate(mld_d_jac_smoother_type :: level%sm, stat=info)
allocate(mld_d_jac_smoother_type :: sm, stat=info)
if (info == 0) allocate(mld_d_ilu_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
endif
case (mld_as_)
if (allocated(level%sm)) then
select type (sm => level%sm)
if (allocated(sm)) then
select type (sms => sm)
class is (mld_d_as_smoother_type)
! do nothing
! do nothing
class default
call level%sm%free(info)
if (info == 0) deallocate(level%sm)
call sm%free(info)
if (info == 0) deallocate(sm)
if (info == 0) allocate(mld_d_as_smoother_type ::&
& level%sm, stat=info)
& sm, stat=info)
if (info == 0) allocate(mld_d_ilu_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
end select
else
allocate(mld_d_as_smoother_type :: level%sm, stat=info)
allocate(mld_d_as_smoother_type :: sm, stat=info)
if (info == 0) allocate(mld_d_ilu_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
endif
case default
@ -459,40 +491,72 @@ contains
! Do nothing and hope for the best :)
!
end select
if (allocated(level%sm)) &
& call level%sm%default()
if (allocated(sm)) &
& call sm%default()
end subroutine inner_set_smoother
end subroutine onelev_set_smoother
subroutine onelev_set_solver(level,val,info,pos)
type(mld_d_onelev_type), intent(inout) :: level
class(mld_d_onelev_type), intent(inout) :: level
integer(psb_ipk_), intent(in) :: val
integer(psb_ipk_), intent(out) :: info
character(len=*), optional, intent(in) :: pos
character(len=*), optional, intent(in) :: pos
! Local
integer(psb_ipk_) :: ipos_
info = psb_success_
if (present(pos)) then
select case(psb_toupper(trim(pos)))
case('PRE')
ipos_ = mld_pre_smooth_
case('POST')
ipos_ = mld_post_smooth_
case default
ipos_ = mld_pre_smooth_
end select
else
ipos_ = mld_pre_smooth_
end if
select case(ipos_)
case(mld_pre_smooth_)
call inner_set_solver(level%sm,val,info)
case (mld_post_smooth_)
call inner_set_solver(level%sm2a,val,info)
case default
! Impossible!!
info = psb_err_internal_error_
end select
end subroutine onelev_set_solver
subroutine inner_set_solver(sm,val,info)
class(mld_d_base_smoother_type), allocatable, intent(inout) :: sm
integer(psb_ipk_), intent(in) :: val
integer(psb_ipk_), intent(out) :: info
!
! This here requires a bit more attention.
! Yes, the first argument is a smoother, to catch the case where
! user is trying to set a solver on a not-yet-allocated smoother.
!
select case (val)
case (mld_f_none_)
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) then
select type (sv => level%sm%sv)
if (allocated(sm)) then
if (allocated(sm%sv)) then
select type (sv => sm%sv)
class is (mld_d_id_solver_type)
! do nothing
class default
call level%sm%sv%free(info)
if (info == 0) deallocate(level%sm%sv)
call sm%sv%free(info)
if (info == 0) deallocate(sm%sv)
if (info == 0) allocate(mld_d_id_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
end select
else
allocate(mld_d_id_solver_type :: level%sm%sv, stat=info)
allocate(mld_d_id_solver_type :: sm%sv, stat=info)
endif
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) &
& call level%sm%sv%default()
if (allocated(sm)) then
if (allocated(sm%sv)) &
& call sm%sv%default()
end if
else
write(0,*) 'Calling set_solver without a smoother?'
@ -500,23 +564,23 @@ contains
end if
case (mld_diag_scale_)
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) then
select type (sv => level%sm%sv)
if (allocated(sm)) then
if (allocated(sm%sv)) then
select type (sv => sm%sv)
class is (mld_d_diag_solver_type)
! do nothing
class default
call level%sm%sv%free(info)
if (info == 0) deallocate(level%sm%sv)
call sm%sv%free(info)
if (info == 0) deallocate(sm%sv)
if (info == 0) allocate(mld_d_diag_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
end select
else
allocate(mld_d_diag_solver_type :: level%sm%sv, stat=info)
allocate(mld_d_diag_solver_type :: sm%sv, stat=info)
endif
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) &
& call level%sm%sv%default()
if (allocated(sm)) then
if (allocated(sm%sv)) &
& call sm%sv%default()
end if
else
write(0,*) 'Calling set_solver without a smoother?'
@ -524,22 +588,22 @@ contains
end if
case (mld_gs_)
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) then
select type (sv => level%sm%sv)
class is (mld_d_gs_solver_type)
! do nothing
class default
call level%sm%sv%free(info)
if (info == 0) deallocate(level%sm%sv)
if (allocated(sm)) then
if (allocated(sm%sv)) then
select type (sv => sm%sv)
class is (mld_d_gs_solver_type)
! do nothing
class default
call sm%sv%free(info)
if (info == 0) deallocate(sm%sv)
if (info == 0) allocate(mld_d_gs_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
end select
else
allocate(mld_d_gs_solver_type :: level%sm%sv, stat=info)
allocate(mld_d_gs_solver_type :: sm%sv, stat=info)
endif
if (allocated(level%sm%sv)) then
call level%sm%sv%default()
if (allocated(sm%sv)) then
call sm%sv%default()
else
endif
@ -549,25 +613,25 @@ contains
end if
case (mld_ilu_n_,mld_milu_n_,mld_ilu_t_)
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) then
select type (sv => level%sm%sv)
if (allocated(sm)) then
if (allocated(sm%sv)) then
select type (sv => sm%sv)
class is (mld_d_ilu_solver_type)
! do nothing
class default
call level%sm%sv%free(info)
if (info == 0) deallocate(level%sm%sv)
call sm%sv%free(info)
if (info == 0) deallocate(sm%sv)
if (info == 0) allocate(mld_d_ilu_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
end select
else
allocate(mld_d_ilu_solver_type :: level%sm%sv, stat=info)
allocate(mld_d_ilu_solver_type :: sm%sv, stat=info)
endif
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) &
& call level%sm%sv%default()
if (allocated(sm)) then
if (allocated(sm%sv)) &
& call sm%sv%default()
end if
call level%sm%sv%set('SUB_SOLVE',val,info)
call sm%sv%set('SUB_SOLVE',val,info)
else
write(0,*) 'Calling set_solver without a smoother?'
info = -5
@ -575,23 +639,23 @@ contains
#ifdef HAVE_SLU_
case (mld_slu_)
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) then
select type (sv => level%sm%sv)
if (allocated(sm)) then
if (allocated(sm%sv)) then
select type (sv => sm%sv)
class is (mld_d_slu_solver_type)
! do nothing
class default
call level%sm%sv%free(info)
if (info == 0) deallocate(level%sm%sv)
call sm%sv%free(info)
if (info == 0) deallocate(sm%sv)
if (info == 0) allocate(mld_d_slu_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
end select
else
allocate(mld_d_slu_solver_type :: level%sm%sv, stat=info)
allocate(mld_d_slu_solver_type :: sm%sv, stat=info)
endif
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) &
& call level%sm%sv%default()
if (allocated(sm)) then
if (allocated(sm%sv)) &
& call sm%sv%default()
end if
else
write(0,*) 'Calling set_solver without a smoother?'
@ -600,44 +664,44 @@ contains
#endif
#ifdef HAVE_MUMPS_
case (mld_mumps_)
if (allocated(level%sm%sv)) then
select type (sv => level%sm%sv)
if (allocated(sm%sv)) then
select type (sv => sm%sv)
class is (mld_d_mumps_solver_type)
! do nothing
! do nothing
class default
call level%sm%sv%free(info)
if (info == 0) deallocate(level%sm%sv)
call sm%sv%free(info)
if (info == 0) deallocate(sm%sv)
if (info == 0) allocate(mld_d_mumps_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
end select
else
allocate(mld_d_mumps_solver_type :: level%sm%sv, stat=info)
allocate(mld_d_mumps_solver_type :: sm%sv, stat=info)
endif
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) &
& call level%sm%sv%default()
if (allocated(sm)) then
if (allocated(sm%sv)) &
& call sm%sv%default()
end if
#endif
#ifdef HAVE_UMF_
case (mld_umf_)
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) then
select type (sv => level%sm%sv)
if (allocated(sm)) then
if (allocated(sm%sv)) then
select type (sv => sm%sv)
class is (mld_d_umf_solver_type)
! do nothing
class default
call level%sm%sv%free(info)
if (info == 0) deallocate(level%sm%sv)
call sm%sv%free(info)
if (info == 0) deallocate(sm%sv)
if (info == 0) allocate(mld_d_umf_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
end select
else
allocate(mld_d_umf_solver_type :: level%sm%sv, stat=info)
allocate(mld_d_umf_solver_type :: sm%sv, stat=info)
endif
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) &
& call level%sm%sv%default()
if (allocated(sm)) then
if (allocated(sm%sv)) &
& call sm%sv%default()
end if
else
write(0,*) 'Calling set_solver without a smoother?'
@ -646,23 +710,23 @@ contains
#endif
#ifdef HAVE_SLUDIST_
case (mld_sludist_)
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) then
select type (sv => level%sm%sv)
if (allocated(sm)) then
if (allocated(sm%sv)) then
select type (sv => sm%sv)
class is (mld_d_sludist_solver_type)
! do nothing
class default
call level%sm%sv%free(info)
if (info == 0) deallocate(level%sm%sv)
call sm%sv%free(info)
if (info == 0) deallocate(sm%sv)
if (info == 0) allocate(mld_d_sludist_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
end select
else
allocate(mld_d_sludist_solver_type :: level%sm%sv, stat=info)
allocate(mld_d_sludist_solver_type :: sm%sv, stat=info)
endif
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) &
& call level%sm%sv%default()
if (allocated(sm)) then
if (allocated(sm%sv)) &
& call sm%sv%default()
end if
else
write(0,*) 'Calling set_solver without a smoother?'
@ -674,9 +738,7 @@ contains
! Do nothing and hope for the best :)
!
end select
end subroutine onelev_set_solver
end subroutine inner_set_solver
end subroutine mld_dcprecseti

@ -360,91 +360,123 @@ subroutine mld_dprecseti(p,what,val,info,ilev,pos)
contains
subroutine onelev_set_smoother(level,val,info,pos)
type(mld_d_onelev_type), intent(inout) :: level
class(mld_d_onelev_type), intent(inout) :: level
integer(psb_ipk_), intent(in) :: val
integer(psb_ipk_), intent(out) :: info
character(len=*), optional, intent(in) :: pos
character(len=*), optional, intent(in) :: pos
! Local
integer(psb_ipk_) :: ipos_
info = psb_success_
if (present(pos)) then
select case(psb_toupper(trim(pos)))
case('PRE')
ipos_ = mld_pre_smooth_
case('POST')
ipos_ = mld_post_smooth_
case default
ipos_ = mld_pre_smooth_
end select
else
ipos_ = mld_pre_smooth_
end if
select case(ipos_)
case(mld_pre_smooth_)
call inner_set_smoother(level%sm,val,info)
case (mld_post_smooth_)
call inner_set_smoother(level%sm2a,val,info)
case default
! Impossible!!
info = psb_err_internal_error_
end select
end subroutine onelev_set_smoother
subroutine inner_set_smoother(sm,val,info)
class(mld_d_base_smoother_type), allocatable, intent(inout) :: sm
integer(psb_ipk_), intent(in) :: val
integer(psb_ipk_), intent(out) :: info
!
! This here requires a bit more attention.
!
select case (val)
case (mld_noprec_)
if (allocated(level%sm)) then
select type (sm => level%sm)
if (allocated(sm)) then
select type (sms => sm)
type is (mld_d_base_smoother_type)
! do nothing
class default
call level%sm%free(info)
if (info == 0) deallocate(level%sm)
call sm%free(info)
if (info == 0) deallocate(sm)
if (info == 0) allocate(mld_d_base_smoother_type ::&
& level%sm, stat=info)
& sm, stat=info)
if (info == 0) allocate(mld_d_id_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
end select
else
allocate(mld_d_base_smoother_type ::&
& level%sm, stat=info)
& sm, stat=info)
if (info ==0) allocate(mld_d_id_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
endif
case (mld_jac_)
if (allocated(level%sm)) then
select type (sm => level%sm)
if (allocated(sm)) then
select type (sms => sm)
class is (mld_d_jac_smoother_type)
! do nothing
! do nothing
class default
call level%sm%free(info)
if (info == 0) deallocate(level%sm)
call sm%free(info)
if (info == 0) deallocate(sm)
if (info == 0) allocate(mld_d_jac_smoother_type :: &
& level%sm, stat=info)
& sm, stat=info)
if (info == 0) allocate(mld_d_diag_solver_type :: &
& level%sm%sv, stat=info)
& sm%sv, stat=info)
end select
else
allocate(mld_d_jac_smoother_type :: level%sm, stat=info)
allocate(mld_d_jac_smoother_type :: sm, stat=info)
if (info == 0) allocate(mld_d_diag_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
endif
case (mld_bjac_)
if (allocated(level%sm)) then
select type (sm => level%sm)
if (allocated(sm)) then
select type (sms => sm)
class is (mld_d_jac_smoother_type)
! do nothing
! do nothing
class default
call level%sm%free(info)
if (info == 0) deallocate(level%sm)
call sm%free(info)
if (info == 0) deallocate(sm)
if (info == 0) allocate(mld_d_jac_smoother_type ::&
& level%sm, stat=info)
& sm, stat=info)
if (info == 0) allocate(mld_d_ilu_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
end select
else
allocate(mld_d_jac_smoother_type :: level%sm, stat=info)
allocate(mld_d_jac_smoother_type :: sm, stat=info)
if (info == 0) allocate(mld_d_ilu_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
endif
case (mld_as_)
if (allocated(level%sm)) then
select type (sm => level%sm)
if (allocated(sm)) then
select type (sms => sm)
class is (mld_d_as_smoother_type)
! do nothing
! do nothing
class default
call level%sm%free(info)
if (info == 0) deallocate(level%sm)
call sm%free(info)
if (info == 0) deallocate(sm)
if (info == 0) allocate(mld_d_as_smoother_type ::&
& level%sm, stat=info)
& sm, stat=info)
if (info == 0) allocate(mld_d_ilu_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
end select
else
allocate(mld_d_as_smoother_type :: level%sm, stat=info)
allocate(mld_d_as_smoother_type :: sm, stat=info)
if (info == 0) allocate(mld_d_ilu_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
endif
case default
@ -452,65 +484,96 @@ contains
! Do nothing and hope for the best :)
!
end select
if (allocated(level%sm)) &
& call level%sm%default()
if (allocated(sm)) &
& call sm%default()
end subroutine inner_set_smoother
end subroutine onelev_set_smoother
subroutine onelev_set_solver(level,val,info,pos)
type(mld_d_onelev_type), intent(inout) :: level
class(mld_d_onelev_type), intent(inout) :: level
integer(psb_ipk_), intent(in) :: val
integer(psb_ipk_), intent(out) :: info
character(len=*), optional, intent(in) :: pos
! Local
integer(psb_ipk_) :: ipos_
info = psb_success_
if (present(pos)) then
select case(psb_toupper(trim(pos)))
case('PRE')
ipos_ = mld_pre_smooth_
case('POST')
ipos_ = mld_post_smooth_
case default
ipos_ = mld_pre_smooth_
end select
else
ipos_ = mld_pre_smooth_
end if
select case(ipos_)
case(mld_pre_smooth_)
call inner_set_solver(level%sm,val,info)
case (mld_post_smooth_)
call inner_set_solver(level%sm2a,val,info)
case default
! Impossible!!
info = psb_err_internal_error_
end select
end subroutine onelev_set_solver
subroutine inner_set_solver(sm,val,info)
class(mld_d_base_smoother_type), allocatable, intent(inout) :: sm
integer(psb_ipk_), intent(in) :: val
integer(psb_ipk_), intent(out) :: info
!
! This here requires a bit more attention.
! Yes, the first argument is a smoother, to catch the case where
! user is trying to set a solver on a not-yet-allocated smoother.
!
select case (val)
case (mld_f_none_)
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) then
select type (sv => level%sm%sv)
class is (mld_d_id_solver_type)
! do nothing
class default
call level%sm%sv%free(info)
if (info == 0) deallocate(level%sm%sv)
if (allocated(sm)) then
if (allocated(sm%sv)) then
select type (sv => sm%sv)
class is (mld_d_id_solver_type)
! do nothing
class default
call sm%sv%free(info)
if (info == 0) deallocate(sm%sv)
if (info == 0) allocate(mld_d_id_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
end select
else
allocate(mld_d_id_solver_type :: level%sm%sv, stat=info)
allocate(mld_d_id_solver_type :: sm%sv, stat=info)
endif
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) &
& call level%sm%sv%default()
if (allocated(sm)) then
if (allocated(sm%sv)) &
& call sm%sv%default()
end if
else
write(0,*) 'Calling set_solver without a smoother?'
info = -5
end if
case (mld_diag_scale_)
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) then
select type (sv => level%sm%sv)
class is (mld_d_diag_solver_type)
! do nothing
class default
call level%sm%sv%free(info)
if (info == 0) deallocate(level%sm%sv)
if (allocated(sm)) then
if (allocated(sm%sv)) then
select type (sv => sm%sv)
class is (mld_d_diag_solver_type)
! do nothing
class default
call sm%sv%free(info)
if (info == 0) deallocate(sm%sv)
if (info == 0) allocate(mld_d_diag_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
end select
else
allocate(mld_d_diag_solver_type :: level%sm%sv, stat=info)
allocate(mld_d_diag_solver_type :: sm%sv, stat=info)
endif
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) &
& call level%sm%sv%default()
if (allocated(sm)) then
if (allocated(sm%sv)) &
& call sm%sv%default()
end if
else
write(0,*) 'Calling set_solver without a smoother?'
@ -518,23 +581,24 @@ contains
end if
case (mld_gs_)
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) then
select type (sv => level%sm%sv)
class is (mld_d_gs_solver_type)
! do nothing
class default
call level%sm%sv%free(info)
if (info == 0) deallocate(level%sm%sv)
if (allocated(sm)) then
if (allocated(sm%sv)) then
select type (sv => sm%sv)
class is (mld_d_gs_solver_type)
! do nothing
class default
call sm%sv%free(info)
if (info == 0) deallocate(sm%sv)
if (info == 0) allocate(mld_d_gs_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
end select
else
allocate(mld_d_gs_solver_type :: level%sm%sv, stat=info)
allocate(mld_d_gs_solver_type :: sm%sv, stat=info)
endif
if (allocated(sm%sv)) then
call sm%sv%default()
else
endif
if (allocated(level%sm%sv)) then
call level%sm%sv%default()
end if
else
write(0,*) 'Calling set_solver without a smoother?'
@ -542,133 +606,132 @@ contains
end if
case (mld_ilu_n_,mld_milu_n_,mld_ilu_t_)
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) then
select type (sv => level%sm%sv)
class is (mld_d_ilu_solver_type)
! do nothing
class default
call level%sm%sv%free(info)
if (info == 0) deallocate(level%sm%sv)
if (allocated(sm)) then
if (allocated(sm%sv)) then
select type (sv => sm%sv)
class is (mld_d_ilu_solver_type)
! do nothing
class default
call sm%sv%free(info)
if (info == 0) deallocate(sm%sv)
if (info == 0) allocate(mld_d_ilu_solver_type ::&
& level%sm%sv, stat=info)
& sm%sv, stat=info)
end select
else
allocate(mld_d_ilu_solver_type :: level%sm%sv, stat=info)
allocate(mld_d_ilu_solver_type :: sm%sv, stat=info)
endif
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) &
& call level%sm%sv%default()
if (allocated(sm)) then
if (allocated(sm%sv)) &
& call sm%sv%default()
end if
call level%sm%sv%set(mld_sub_solve_,val,info)
call sm%sv%set('SUB_SOLVE',val,info)
else
write(0,*) 'Calling set_solver without a smoother?'
info = -5
end if
#ifdef HAVE_UMF_
case (mld_umf_)
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) then
select type (sv => level%sm%sv)
class is (mld_d_umf_solver_type)
#ifdef HAVE_SLU_
case (mld_slu_)
if (allocated(sm)) then
if (allocated(sm%sv)) then
select type (sv => sm%sv)
class is (mld_d_slu_solver_type)
! do nothing
class default
call level%sm%sv%free(info)
if (info == 0) deallocate(level%sm%sv)
if (info == 0) allocate(mld_d_umf_solver_type ::&
& level%sm%sv, stat=info)
call sm%sv%free(info)
if (info == 0) deallocate(sm%sv)
if (info == 0) allocate(mld_d_slu_solver_type ::&
& sm%sv, stat=info)
end select
else
allocate(mld_d_umf_solver_type :: level%sm%sv, stat=info)
allocate(mld_d_slu_solver_type :: sm%sv, stat=info)
endif
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) &
& call level%sm%sv%default()
if (allocated(sm)) then
if (allocated(sm%sv)) &
& call sm%sv%default()
end if
else
write(0,*) 'Calling set_solver without a smoother?'
info = -5
end if
#endif
#ifdef HAVE_SLUDIST_
case (mld_sludist_)
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) then
select type (sv => level%sm%sv)
class is (mld_d_sludist_solver_type)
#ifdef HAVE_MUMPS_
case (mld_mumps_)
if (allocated(sm%sv)) then
select type (sv => sm%sv)
class is (mld_d_mumps_solver_type)
! do nothing
class default
call sm%sv%free(info)
if (info == 0) deallocate(sm%sv)
if (info == 0) allocate(mld_d_mumps_solver_type ::&
& sm%sv, stat=info)
end select
else
allocate(mld_d_mumps_solver_type :: sm%sv, stat=info)
endif
if (allocated(sm)) then
if (allocated(sm%sv)) &
& call sm%sv%default()
end if
#endif
#ifdef HAVE_UMF_
case (mld_umf_)
if (allocated(sm)) then
if (allocated(sm%sv)) then
select type (sv => sm%sv)
class is (mld_d_umf_solver_type)
! do nothing
class default
call level%sm%sv%free(info)
if (info == 0) deallocate(level%sm%sv)
if (info == 0) allocate(mld_d_sludist_solver_type ::&
& level%sm%sv, stat=info)
call sm%sv%free(info)
if (info == 0) deallocate(sm%sv)
if (info == 0) allocate(mld_d_umf_solver_type ::&
& sm%sv, stat=info)
end select
else
allocate(mld_d_sludist_solver_type :: level%sm%sv, stat=info)
allocate(mld_d_umf_solver_type :: sm%sv, stat=info)
endif
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) &
& call level%sm%sv%default()
if (allocated(sm)) then
if (allocated(sm%sv)) &
& call sm%sv%default()
end if
else
write(0,*) 'Calling set_solver without a smoother?'
info = -5
end if
#endif
#ifdef HAVE_SLU_
case (mld_slu_)
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) then
select type (sv => level%sm%sv)
class is (mld_d_slu_solver_type)
#ifdef HAVE_SLUDIST_
case (mld_sludist_)
if (allocated(sm)) then
if (allocated(sm%sv)) then
select type (sv => sm%sv)
class is (mld_d_sludist_solver_type)
! do nothing
class default
call level%sm%sv%free(info)
if (info == 0) deallocate(level%sm%sv)
if (info == 0) allocate(mld_d_slu_solver_type ::&
& level%sm%sv, stat=info)
call sm%sv%free(info)
if (info == 0) deallocate(sm%sv)
if (info == 0) allocate(mld_d_sludist_solver_type ::&
& sm%sv, stat=info)
end select
else
allocate(mld_d_slu_solver_type :: level%sm%sv, stat=info)
allocate(mld_d_sludist_solver_type :: sm%sv, stat=info)
endif
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) &
& call level%sm%sv%default()
if (allocated(sm)) then
if (allocated(sm%sv)) &
& call sm%sv%default()
end if
else
write(0,*) 'Calling set_solver without a smoother?'
info = -5
end if
#endif
#ifdef HAVE_MUMPS_
case (mld_mumps_)
if (allocated(level%sm%sv)) then
select type (sv => level%sm%sv)
class is (mld_d_mumps_solver_type)
! do nothing
class default
call level%sm%sv%free(info)
if (info == 0) deallocate(level%sm%sv)
if (info == 0) allocate(mld_d_mumps_solver_type ::&
& level%sm%sv, stat=info)
end select
else
allocate(mld_d_mumps_solver_type :: level%sm%sv, stat=info)
endif
if (allocated(level%sm)) then
if (allocated(level%sm%sv)) then
call level%sm%sv%default()
end if
end if
#endif
case default
!
! Do nothing and hope for the best :)
!
end select
end subroutine onelev_set_solver
end subroutine inner_set_solver
end subroutine mld_dprecseti
@ -734,32 +797,35 @@ subroutine mld_dprecsetsm(p,val,info,ilev,pos)
case(mld_pre_smooth_)
do ilev_ = ilmin, ilmax
if (allocated(p%precv(ilev_)%sm)) then
if (allocated(p%precv(ilev_)%sm%sv)) then
deallocate(p%precv(ilev_)%sm%sv)
endif
deallocate(p%precv(ilev_)%sm)
end if
if (.not.same_type_as(p%precv(ilev_)%sm,val)) then
call p%precv(ilev_)%sm%free(info)
deallocate(p%precv(ilev_)%sm, stat=info)
end if
endif
if (.not.allocated(p%precv(ilev_)%sm)) then
#ifdef HAVE_MOLD
allocate(p%precv(ilev_)%sm,mold=val)
allocate(p%precv(ilev_)%sm,mold=val)
#else
allocate(p%precv(ilev_)%sm,source=val)
allocate(p%precv(ilev_)%sm,source=val)
#endif
end if
call p%precv(ilev_)%sm%default()
p%precv(ilev_)%sm2 => p%precv(ilev_)%sm
end do
case(mld_post_smooth_)
do ilev_ = ilmin, ilmax
if (allocated(p%precv(ilev_)%sm2a)) then
if (allocated(p%precv(ilev_)%sm2a%sv)) then
deallocate(p%precv(ilev_)%sm2a%sv)
endif
deallocate(p%precv(ilev_)%sm2a)
end if
if (.not.same_type_as(p%precv(ilev_)%sm2a,val)) then
call p%precv(ilev_)%sm2a%free(info)
deallocate(p%precv(ilev_)%sm2a, stat=info)
endif
if (.not.allocated(p%precv(ilev_)%sm2a)) then
#ifdef HAVE_MOLD
allocate(p%precv(ilev_)%sm2a,mold=val)
allocate(p%precv(ilev_)%sm2a,mold=val)
#else
allocate(p%precv(ilev_)%sm2a,source=val)
allocate(p%precv(ilev_)%sm2a,source=val)
#endif
end if
call p%precv(ilev_)%sm2a%default()
p%precv(ilev_)%sm2 => p%precv(ilev_)%sm2a
end do
@ -834,6 +900,7 @@ subroutine mld_dprecsetsv(p,val,info,ilev,pos)
if (allocated(p%precv(ilev_)%sm)) then
if (allocated(p%precv(ilev_)%sm%sv)) then
if (.not.same_type_as(p%precv(ilev_)%sm%sv,val)) then
call p%precv(ilev_)%sm%sv%free(info)
deallocate(p%precv(ilev_)%sm%sv,stat=info)
if (info /= 0) then
info = 3111
@ -870,6 +937,7 @@ subroutine mld_dprecsetsv(p,val,info,ilev,pos)
if (allocated(p%precv(ilev_)%sm2a)) then
if (allocated(p%precv(ilev_)%sm2a%sv)) then
if (.not.same_type_as(p%precv(ilev_)%sm2a%sv,val)) then
call p%precv(ilev_)%sm2a%sv%free(info)
deallocate(p%precv(ilev_)%sm2a%sv,stat=info)
if (info /= 0) then
info = 3111

@ -115,7 +115,7 @@ contains
integer(psb_ipk_), intent(out) :: info
character(len=*), optional, intent(in) :: pos
call p%set(what,val,info)
call p%set(what,val,info,pos=pos)
end subroutine mld_d_iprecseti
subroutine mld_d_iprecsetr(p,what,val,info,pos)
@ -125,7 +125,7 @@ contains
integer(psb_ipk_), intent(out) :: info
character(len=*), optional, intent(in) :: pos
call p%set(what,val,info)
call p%set(what,val,info,pos=pos)
end subroutine mld_d_iprecsetr
subroutine mld_d_iprecsetc(p,what,val,info,pos)
@ -135,7 +135,7 @@ contains
integer(psb_ipk_), intent(out) :: info
character(len=*), optional, intent(in) :: pos
call p%set(what,val,info)
call p%set(what,val,info,pos=pos)
end subroutine mld_d_iprecsetc
subroutine mld_d_cprecseti(p,what,val,info,pos)
@ -145,7 +145,7 @@ contains
integer(psb_ipk_), intent(out) :: info
character(len=*), optional, intent(in) :: pos
call p%set(what,val,info)
call p%set(what,val,info,pos=pos)
end subroutine mld_d_cprecseti
subroutine mld_d_cprecsetr(p,what,val,info,pos)
@ -155,7 +155,7 @@ contains
integer(psb_ipk_), intent(out) :: info
character(len=*), optional, intent(in) :: pos
call p%set(what,val,info)
call p%set(what,val,info,pos=pos)
end subroutine mld_d_cprecsetr
subroutine mld_d_cprecsetc(p,what,val,info,pos)
@ -165,7 +165,7 @@ contains
integer(psb_ipk_), intent(out) :: info
character(len=*), optional, intent(in) :: pos
call p%set(what,val,info)
call p%set(what,val,info,pos=pos)
end subroutine mld_d_cprecsetc
end module mld_d_prec_mod

@ -270,6 +270,7 @@ program ppde3d
call mld_precset(prec,'coarse_aggr_size', prectype%csize, info)
call prec%set(dbsmth,info,pos='post')
call prec%set(dbwgs,info,pos='post')
call mld_precset(prec,'solver_sweeps', prectype%svsweeps, info, pos='post')
else
nlv = 1
call mld_precinit(prec,prectype%prec, info, nlev=nlv)
@ -281,6 +282,9 @@ program ppde3d
call mld_precset(prec,'sub_fillin', prectype%fill1, info)
call mld_precset(prec,'solver_sweeps', prectype%svsweeps, info)
call mld_precset(prec,'sub_iluthrs', prectype%thr1, info)
call prec%set(dbsmth,info,pos='post')
call prec%set(dbwgs,info,pos='post')
call mld_precset(prec,'solver_sweeps', prectype%svsweeps, info, pos='post')
end if
call psb_barrier(ictxt)

@ -1,21 +1,21 @@
RGMRES ! Iterative method: BiCGSTAB BiCG CGS RGMRES BiCGSTABL CG
CG ! Iterative method: BiCGSTAB BiCG CGS RGMRES BiCGSTABL CG
CSR ! Storage format CSR COO JAD
0030 ! IDIM; domain size is idim**3
2 ! ISTOPC
0100 ! ITMAX
0500 ! ITMAX
10 ! ITRACE
30 ! IRST (restart for RGMRES and BiCGSTABL)
1.d-6 ! EPS
3L-MUL-RAS-BJAC4-ILU ! Descriptive name for preconditioner (up to 40 chars)
ML ! Preconditioner NONE JACOBI BJAC AS ML
1 ! Number of overlap layers for AS preconditioner at finest level
0 ! Number of overlap layers for AS preconditioner at finest level
HALO ! Restriction operator NONE HALO
NONE ! Prolongation operator NONE SUM AVG
GS ! Subdomain solver DSCALE ILU MILU ILUT UMF SLU
2 ! sweeps for GS
4 ! Solver sweeps for GS
0 ! Level-set N for ILU(N), and P for ILUT
1.d-4 ! Threshold T for ILU(T,P)
1 ! Smoother/Jacobi sweeps
2 ! Smoother/Jacobi sweeps
BJAC ! Smoother type JACOBI BJAC AS; ignored for non-ML
2 ! Number of levels in a multilevel preconditioner
SMOOTHED ! Kind of aggregation: SMOOTHED, NONSMOOTHED, MINENERGY
@ -25,7 +25,7 @@ MULT ! Type of multilevel correction: ADD MULT
TWOSIDE ! Side of correction PRE POST TWOSIDE (ignored for ADD)
DIST ! Coarse level: matrix distribution DIST REPL
BJAC ! Coarse level: solver JACOBI BJAC UMF SLU SLUDIST MUMPS
GS ! Coarse level: subsolver DSCALE ILU UMF SLU SLUDIST MUMPS
ILU ! Coarse level: subsolver DSCALE ILU UMF SLU SLUDIST MUMPS
1 ! Coarse level: Level-set N for ILU(N)
1.d-4 ! Coarse level: Threshold T for ILU(T,P)
4 ! Coarse level: Number of Jacobi sweeps

Loading…
Cancel
Save