Folded in new precinit/precset.

stopcriterion
Salvatore Filippone 18 years ago
parent a59f22b619
commit 206b19c1fc

@ -0,0 +1,199 @@
!!$
!!$
!!$ MD2P4
!!$ Multilevel Domain Decomposition Parallel Preconditioner Package for PSBLAS
!!$ for
!!$ Parallel Sparse BLAS v2.0
!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari University of Rome Tor Vergata
!!$ Daniela di Serafino Second University of Naples
!!$ Pasqua D'Ambra ICAR-CNR
!!$
!!$ 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 MD2P4 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 MD2P4 GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_dprecinit(p,ptype,info,nlev)
use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_dprecinit
implicit none
type(psb_dprec_type), intent(inout) :: p
character(len=*), intent(in) :: ptype
integer, intent(out) :: info
integer, optional, intent(in) :: nlev
integer :: nlev_, ilev_
info = 0
if (allocated(p%baseprecv)) then
call psb_precfree(p,info)
if (info /=0) then
! Do we want to do something?
endif
endif
select case(toupper(ptype(1:len_trim(ptype))))
case ('NONE','NOPREC')
nlev_ = 1
ilev_ = 1
allocate(p%baseprecv(nlev_),stat=info)
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(p_type_) = noprec_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_none_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
case ('DIAG')
nlev_ = 1
ilev_ = 1
allocate(p%baseprecv(nlev_),stat=info)
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(p_type_) = diag_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_none_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
case ('BJAC')
nlev_ = 1
ilev_ = 1
allocate(p%baseprecv(nlev_),stat=info)
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(p_type_) = bjac_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
case ('ASM','AS')
nlev_ = 1
ilev_ = 1
allocate(p%baseprecv(nlev_),stat=info)
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(p_type_) = asm_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_halo_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 1
p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
case ('MLD', 'ML')
if (present(nlev)) then
nlev_ = max(1,nlev)
else
nlev_ = 2
end if
if (nlev_ == 1) then
write(0,*) 'Warning: requested ML preconditioner with NLEV=1'
endif
ilev_ = 1
allocate(p%baseprecv(nlev_),stat=info)
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(p_type_) = asm_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_halo_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 1
p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
if (nlev_ == 1) return
do ilev_ = 2, nlev_ -1
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(p_type_) = bjac_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(ml_type_) = mult_ml_prec_
p%baseprecv(ilev_)%iprcparm(aggr_alg_) = loc_aggr_
p%baseprecv(ilev_)%iprcparm(smth_kind_) = smth_omg_
p%baseprecv(ilev_)%iprcparm(coarse_mat_) = mat_distr_
p%baseprecv(ilev_)%iprcparm(smth_pos_) = post_smooth_
p%baseprecv(ilev_)%iprcparm(glb_smth_) = 1
p%baseprecv(ilev_)%iprcparm(om_choice_) = lib_choice_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
p%baseprecv(ilev_)%dprcparm(smooth_omega_) = 4.d0/3.d0
end do
ilev_ = nlev_
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(p_type_) = bjac_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(ml_type_) = mult_ml_prec_
p%baseprecv(ilev_)%iprcparm(aggr_alg_) = loc_aggr_
p%baseprecv(ilev_)%iprcparm(smth_kind_) = smth_omg_
p%baseprecv(ilev_)%iprcparm(coarse_mat_) = mat_distr_
p%baseprecv(ilev_)%iprcparm(smth_pos_) = post_smooth_
p%baseprecv(ilev_)%iprcparm(glb_smth_) = 1
p%baseprecv(ilev_)%iprcparm(om_choice_) = lib_choice_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_umf_
p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 4
p%baseprecv(ilev_)%dprcparm(smooth_omega_) = 4.d0/3.d0
case default
write(0,*) 'Unknown preconditioner type request "',ptype,'"'
info = 2
end select
end subroutine psb_dprecinit

@ -34,156 +34,126 @@
!!$ POSSIBILITY OF SUCH DAMAGE. !!$ POSSIBILITY OF SUCH DAMAGE.
!!$ !!$
!!$ !!$
subroutine psb_dprecset(p,ptype,info,iv,rs,rv,ilev,nlev) subroutine psb_dprecseti(p,what,val,info,ilev)
use psb_base_mod use psb_base_mod
use psb_prec_mod, mld_protect_name => psb_dprecset use psb_prec_mod, psb_protect_name => psb_dprecseti
implicit none implicit none
type(psb_dprec_type), intent(inout) :: p type(psb_dprec_type), intent(inout) :: p
character(len=*), intent(in) :: ptype integer, intent(in) :: what
integer, intent(in) :: val
integer, intent(out) :: info integer, intent(out) :: info
integer, optional, intent(in) :: iv(:) integer, optional, intent(in) :: ilev
integer, optional, intent(in) :: nlev,ilev integer :: ilev_, nlev_
real(kind(1.d0)), optional, intent(in) :: rs
real(kind(1.d0)), optional, intent(in) :: rv(:)
character(len=len(ptype)) :: typeup
integer :: isz, err, nlev_, ilev_, i
info = 0 info = 0
if (present(ilev)) then if (present(ilev)) then
ilev_ = max(1, ilev) ilev_ = ilev
else else
ilev_ = 1 ilev_ = 1
end if end if
if (present(nlev)) then
if (allocated(p%baseprecv)) then
write(0,*) 'Warning: NLEV is ignored when P is already allocated'
end if
nlev_ = max(1, nlev)
else
nlev_ = 1
end if
if (.not.allocated(p%baseprecv)) then if (.not.allocated(p%baseprecv)) then
allocate(p%baseprecv(nlev_),stat=err) write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
else info = -1
nlev_ = size(p%baseprecv) return
endif endif
nlev_ = size(p%baseprecv)
if ((ilev_<1).or.(ilev_ > nlev_)) then if ((ilev_<1).or.(ilev_ > nlev_)) then
write(0,*) 'PRECSET ERRROR: ilev out of bounds' write(0,*) 'PRECSET ERRROR: ilev out of bounds'
info = -1 info = -1
return return
endif endif
if (.not.allocated(p%baseprecv(ilev_)%iprcparm)) then
write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
info = -1
return
endif
call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return if (ilev_ == 1) then
p%baseprecv(ilev_)%iprcparm(:) = 0 ! Rules for fine level are slightly different.
select case(what)
select case(toupper(ptype(1:len_trim(ptype)))) case(p_type_,f_type_,restr_,prol_,iren_,n_ovr_,ilu_fill_in_,jac_sweeps_)
case ('NONE','NOPREC') p%baseprecv(ilev_)%iprcparm(what) = val
p%baseprecv(ilev_)%iprcparm(:) = 0 case default
p%baseprecv(ilev_)%iprcparm(p_type_) = noprec_ write(0,*) 'Error: trying to call PRECSET with an invalid WHAT'
p%baseprecv(ilev_)%iprcparm(f_type_) = f_none_ info = -2
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_ end select
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_ else if (ilev_ > 1) then
p%baseprecv(ilev_)%iprcparm(iren_) = 0 select case(what)
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0 case(p_type_,f_type_,restr_,prol_,iren_,n_ovr_,ilu_fill_in_,jac_sweeps_,&
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1 & ml_type_,aggr_alg_,smth_kind_,coarse_mat_,smth_pos_,glb_smth_,&
& om_choice_)
case ('DIAG') p%baseprecv(ilev_)%iprcparm(what) = val
p%baseprecv(ilev_)%iprcparm(:) = 0 case default
p%baseprecv(ilev_)%iprcparm(p_type_) = diag_ write(0,*) 'Error: trying to call PRECSET with an invalid WHAT'
p%baseprecv(ilev_)%iprcparm(f_type_) = f_none_ info = -2
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_ end select
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
case ('BJAC')
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(p_type_) = bjac_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
case ('ASM','AS')
p%baseprecv(ilev_)%iprcparm(:) = 0
! Defaults first
p%baseprecv(ilev_)%iprcparm(p_type_) = asm_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_halo_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 1
p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
if (present(iv)) then
isz = size(iv)
if (isz >= 1) p%baseprecv(ilev_)%iprcparm(n_ovr_) = iv(1)
if (isz >= 2) p%baseprecv(ilev_)%iprcparm(restr_) = iv(2)
if (isz >= 3) p%baseprecv(ilev_)%iprcparm(prol_) = iv(3)
if (isz >= 4) p%baseprecv(ilev_)%iprcparm(f_type_) = iv(4)
if (isz >= 5) p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = iv(5)
! Do not consider renum for the time being.
!!$ if (isz >= 5) p%baseprecv(ilev_)%iprcparm(iren_) = iv(5)
endif endif
end subroutine psb_dprecseti
subroutine psb_dprecsetd(p,what,val,info,ilev)
use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_dprecsetd
implicit none
type(psb_dprec_type), intent(inout) :: p
integer, intent(in) :: what
real(kind(1.d0)), intent(in) :: val
integer, intent(out) :: info
integer, optional, intent(in) :: ilev
case ('ML', '2L', '2LEV') integer :: ilev_,nlev_
info = 0
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(p_type_) = bjac_ if (present(ilev)) then
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_ ilev_ = ilev
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_ else
p%baseprecv(ilev_)%iprcparm(iren_) = 0 ilev_ = 1
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(ml_type_) = mult_ml_prec_
p%baseprecv(ilev_)%iprcparm(aggr_alg_) = loc_aggr_
p%baseprecv(ilev_)%iprcparm(smth_kind_) = smth_omg_
p%baseprecv(ilev_)%iprcparm(coarse_mat_) = mat_distr_
p%baseprecv(ilev_)%iprcparm(smth_pos_) = post_smooth_
p%baseprecv(ilev_)%iprcparm(glb_smth_) = 1
p%baseprecv(ilev_)%iprcparm(om_choice_) = lib_choice_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(ilev_)%dprcparm(smooth_omega_) = 4.d0/3.d0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
if (present(iv)) then
isz = size(iv)
if (isz >= 1) p%baseprecv(ilev_)%iprcparm(ml_type_) = iv(1)
if (isz >= 2) p%baseprecv(ilev_)%iprcparm(aggr_alg_) = iv(2)
if (isz >= 3) p%baseprecv(ilev_)%iprcparm(coarse_mat_) = iv(3)
if (isz >= 4) p%baseprecv(ilev_)%iprcparm(smth_pos_) = iv(4)
if (isz >= 5) p%baseprecv(ilev_)%iprcparm(f_type_) = iv(5)
if (isz >= 6) p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = iv(6)
if (isz >= 7) p%baseprecv(ilev_)%iprcparm(smth_kind_) = iv(7)
end if end if
if (present(rs)) then if (.not.allocated(p%baseprecv)) then
p%baseprecv(ilev_)%iprcparm(om_choice_) = user_choice_ write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
p%baseprecv(ilev_)%dprcparm(smooth_omega_) = rs info = -1
return
endif endif
nlev_ = size(p%baseprecv)
if ((ilev_<1).or.(ilev_ > nlev_)) then
write(0,*) 'PRECSET ERRROR: ilev out of bounds'
info = -1
return
endif
if (.not.allocated(p%baseprecv(ilev_)%dprcparm)) then
write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
info = -1
return
endif
if (ilev_ == 1) then
! Rules for fine level are slightly different.
select case(what)
! Right now we don't have any at base level. Will change when
! we implement F_ILU_E_
case default case default
write(0,*) 'Unknown preconditioner type request "',ptype,'"' write(0,*) 'Error: trying to call PRECSET with an invalid WHAT'
err = 2 info = -2
end select end select
else if (ilev_ > 1) then
select case(what)
case(smooth_omega_)
p%baseprecv(ilev_)%dprcparm(what) = val
case default
write(0,*) 'Error: trying to call PRECSET with an invalid WHAT'
info = -2
end select
endif
info = err end subroutine psb_dprecsetd
end subroutine psb_dprecset

@ -61,31 +61,64 @@ module psb_prec_mod
end subroutine psb_zprecbld end subroutine psb_zprecbld
end interface end interface
interface psb_precset interface psb_precinit
subroutine psb_dprecset(prec,ptype,info,iv,rs,rv,ilev,nlev) subroutine psb_dprecinit(p,ptype,info,nlev)
use psb_base_mod use psb_base_mod
use psb_prec_type use psb_prec_type
implicit none type(psb_dprec_type), intent(inout) :: p
type(psb_dprec_type), intent(inout) :: prec
character(len=*), intent(in) :: ptype character(len=*), intent(in) :: ptype
integer, intent(out) :: info integer, intent(out) :: info
integer, optional, intent(in) :: iv(:) integer, optional, intent(in) :: nlev
integer, optional, intent(in) :: nlev,ilev end subroutine psb_dprecinit
real(kind(1.d0)), optional, intent(in) :: rs subroutine psb_zprecinit(p,ptype,info,nlev)
real(kind(1.d0)), optional, intent(in) :: rv(:)
end subroutine psb_dprecset
subroutine psb_zprecset(prec,ptype,info,iv,rs,rv,ilev,nlev)
use psb_base_mod use psb_base_mod
use psb_prec_type use psb_prec_type
implicit none type(psb_zprec_type), intent(inout) :: p
type(psb_zprec_type), intent(inout) :: prec
character(len=*), intent(in) :: ptype character(len=*), intent(in) :: ptype
integer, intent(out) :: info integer, intent(out) :: info
integer, optional, intent(in) :: iv(:) integer, optional, intent(in) :: nlev
real(kind(1.d0)), optional, intent(in) :: rs end subroutine psb_zprecinit
real(kind(1.d0)), optional, intent(in) :: rv(:) end interface
integer, optional, intent(in) :: nlev,ilev
end subroutine psb_zprecset interface psb_precset
subroutine psb_dprecseti(p,what,val,info,ilev)
use psb_base_mod
use psb_prec_type
type(psb_dprec_type), intent(inout) :: p
integer, intent(in) :: what
integer, intent(in) :: val
integer, intent(out) :: info
integer, optional, intent(in) :: ilev
end subroutine psb_dprecseti
subroutine psb_dprecsetd(p,what,val,info,ilev)
use psb_base_mod
use psb_prec_type
type(psb_dprec_type), intent(inout) :: p
integer, intent(in) :: what
real(kind(1.d0)), intent(in) :: val
integer, intent(out) :: info
integer, optional, intent(in) :: ilev
end subroutine psb_dprecsetd
subroutine psb_zprecseti(p,what,val,info,ilev)
use psb_base_mod
use psb_prec_type
type(psb_zprec_type), intent(inout) :: p
integer, intent(in) :: what
integer, intent(in) :: val
integer, intent(out) :: info
integer, optional, intent(in) :: ilev
end subroutine psb_zprecseti
subroutine psb_zprecsetd(p,what,val,info,ilev)
use psb_base_mod
use psb_prec_type
type(psb_zprec_type), intent(inout) :: p
integer, intent(in) :: what
real(kind(1.d0)), intent(in) :: val
integer, intent(out) :: info
integer, optional, intent(in) :: ilev
end subroutine psb_zprecsetd
end interface end interface

@ -132,7 +132,7 @@ psb_dumf_factor_(int *n, int *nnz,
if ( *info == UMFPACK_OK ) { if ( *info == UMFPACK_OK ) {
*info = 0; *info = 0;
} else { } else {
printf("umfpack_di_symbolic() error returns INFO= %d\n", *info); printf("umfpack_di_symbolic() returned INFO= %d\n", *info);
*info = -11; *info = -11;
*numptr = (fptr) NULL; *numptr = (fptr) NULL;
return; return;
@ -148,7 +148,7 @@ psb_dumf_factor_(int *n, int *nnz,
*info = 0; *info = 0;
*numptr = (fptr) Numeric; *numptr = (fptr) Numeric;
} else { } else {
printf("umfpack_di_numeric() error returns INFO= %d\n", *info); printf("umfpack_di_numeric() returned INFO= %d\n", *info);
*info = -12; *info = -12;
*numptr = (fptr) NULL; *numptr = (fptr) NULL;
} }

@ -0,0 +1,199 @@
!!$
!!$
!!$ MD2P4
!!$ Multilevel Domain Decomposition Parallel Preconditioner Package for PSBLAS
!!$ for
!!$ Parallel Sparse BLAS v2.0
!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari University of Rome Tor Vergata
!!$ Daniela di Serafino Second University of Naples
!!$ Pasqua D'Ambra ICAR-CNR
!!$
!!$ 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 MD2P4 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 MD2P4 GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_zprecinit(p,ptype,info,nlev)
use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_zprecinit
implicit none
type(psb_zprec_type), intent(inout) :: p
character(len=*), intent(in) :: ptype
integer, intent(out) :: info
integer, optional, intent(in) :: nlev
integer :: nlev_, ilev_
info = 0
if (allocated(p%baseprecv)) then
call psb_precfree(p,info)
if (info /=0) then
! Do we want to do something?
endif
endif
select case(toupper(ptype(1:len_trim(ptype))))
case ('NONE','NOPREC')
nlev_ = 1
ilev_ = 1
allocate(p%baseprecv(nlev_),stat=info)
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(p_type_) = noprec_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_none_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
case ('DIAG')
nlev_ = 1
ilev_ = 1
allocate(p%baseprecv(nlev_),stat=info)
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(p_type_) = diag_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_none_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
case ('BJAC')
nlev_ = 1
ilev_ = 1
allocate(p%baseprecv(nlev_),stat=info)
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(p_type_) = bjac_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
case ('ASM','AS')
nlev_ = 1
ilev_ = 1
allocate(p%baseprecv(nlev_),stat=info)
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(p_type_) = asm_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_halo_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 1
p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
case ('MLD', 'ML')
if (present(nlev)) then
nlev_ = max(1,nlev)
else
nlev_ = 2
end if
if (nlev_ == 1) then
write(0,*) 'Warning: requested ML preconditioner with NLEV=1'
endif
ilev_ = 1
allocate(p%baseprecv(nlev_),stat=info)
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(p_type_) = asm_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_halo_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 1
p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
if (nlev_ == 1) return
do ilev_ = 2, nlev_ -1
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(p_type_) = bjac_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(ml_type_) = mult_ml_prec_
p%baseprecv(ilev_)%iprcparm(aggr_alg_) = loc_aggr_
p%baseprecv(ilev_)%iprcparm(smth_kind_) = smth_omg_
p%baseprecv(ilev_)%iprcparm(coarse_mat_) = mat_distr_
p%baseprecv(ilev_)%iprcparm(smth_pos_) = post_smooth_
p%baseprecv(ilev_)%iprcparm(glb_smth_) = 1
p%baseprecv(ilev_)%iprcparm(om_choice_) = lib_choice_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
p%baseprecv(ilev_)%dprcparm(smooth_omega_) = 4.d0/3.d0
end do
ilev_ = nlev_
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(p_type_) = bjac_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(ml_type_) = mult_ml_prec_
p%baseprecv(ilev_)%iprcparm(aggr_alg_) = loc_aggr_
p%baseprecv(ilev_)%iprcparm(smth_kind_) = smth_omg_
p%baseprecv(ilev_)%iprcparm(coarse_mat_) = mat_distr_
p%baseprecv(ilev_)%iprcparm(smth_pos_) = post_smooth_
p%baseprecv(ilev_)%iprcparm(glb_smth_) = 1
p%baseprecv(ilev_)%iprcparm(om_choice_) = lib_choice_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_umf_
p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 4
p%baseprecv(ilev_)%dprcparm(smooth_omega_) = 4.d0/3.d0
case default
write(0,*) 'Unknown preconditioner type request "',ptype,'"'
info = 2
end select
end subroutine psb_zprecinit

@ -34,156 +34,126 @@
!!$ POSSIBILITY OF SUCH DAMAGE. !!$ POSSIBILITY OF SUCH DAMAGE.
!!$ !!$
!!$ !!$
subroutine psb_zprecset(p,ptype,info,iv,rs,rv,ilev,nlev) subroutine psb_zprecseti(p,what,val,info,ilev)
use psb_base_mod use psb_base_mod
use psb_prec_mod, mld_protect_name => psb_zprecset use psb_prec_mod, psb_protect_name => psb_zprecseti
implicit none implicit none
type(psb_zprec_type), intent(inout) :: p type(psb_zprec_type), intent(inout) :: p
character(len=*), intent(in) :: ptype integer, intent(in) :: what
integer, intent(in) :: val
integer, intent(out) :: info integer, intent(out) :: info
integer, optional, intent(in) :: iv(:) integer, optional, intent(in) :: ilev
integer, optional, intent(in) :: nlev,ilev integer :: ilev_, nlev_
real(kind(1.d0)), optional, intent(in) :: rs
real(kind(1.d0)), optional, intent(in) :: rv(:)
character(len=len(ptype)) :: typeup
integer :: isz, err, nlev_, ilev_, i
info = 0 info = 0
if (present(ilev)) then if (present(ilev)) then
ilev_ = max(1, ilev) ilev_ = ilev
else else
ilev_ = 1 ilev_ = 1
end if end if
if (present(nlev)) then
if (allocated(p%baseprecv)) then
write(0,*) 'Warning: NLEV is ignored when P is already allocated'
end if
nlev_ = max(1, nlev)
else
nlev_ = 1
end if
if (.not.allocated(p%baseprecv)) then if (.not.allocated(p%baseprecv)) then
allocate(p%baseprecv(nlev_),stat=err) write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
else info = -1
nlev_ = size(p%baseprecv) return
endif endif
nlev_ = size(p%baseprecv)
if ((ilev_<1).or.(ilev_ > nlev_)) then if ((ilev_<1).or.(ilev_ > nlev_)) then
write(0,*) 'PRECSET ERRROR: ilev out of bounds' write(0,*) 'PRECSET ERRROR: ilev out of bounds'
info = -1 info = -1
return return
endif endif
if (.not.allocated(p%baseprecv(ilev_)%iprcparm)) then
write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
info = -1
return
endif
call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info) if (ilev_ == 1) then
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info) ! Rules for fine level are slightly different.
if (info /= 0) return select case(what)
p%baseprecv(ilev_)%iprcparm(:) = 0 case(p_type_,f_type_,restr_,prol_,iren_,n_ovr_,ilu_fill_in_,jac_sweeps_)
p%baseprecv(ilev_)%iprcparm(what) = val
select case(toupper(ptype(1:len_trim(ptype)))) case default
case ('NONE','NOPREC') write(0,*) 'Error: trying to call PRECSET with an invalid WHAT'
p%baseprecv(ilev_)%iprcparm(:) = 0 info = -2
p%baseprecv(ilev_)%iprcparm(p_type_) = noprec_ end select
p%baseprecv(ilev_)%iprcparm(f_type_) = f_none_ else if (ilev_ > 1) then
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_ select case(what)
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_ case(p_type_,f_type_,restr_,prol_,iren_,n_ovr_,ilu_fill_in_,jac_sweeps_,&
p%baseprecv(ilev_)%iprcparm(iren_) = 0 & ml_type_,aggr_alg_,smth_kind_,coarse_mat_,smth_pos_,glb_smth_,&
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0 & om_choice_)
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1 p%baseprecv(ilev_)%iprcparm(what) = val
case default
case ('DIAG') write(0,*) 'Error: trying to call PRECSET with an invalid WHAT'
p%baseprecv(ilev_)%iprcparm(:) = 0 info = -2
p%baseprecv(ilev_)%iprcparm(p_type_) = diag_ end select
p%baseprecv(ilev_)%iprcparm(f_type_) = f_none_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
case ('BJAC')
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(p_type_) = bjac_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
case ('ASM','AS')
p%baseprecv(ilev_)%iprcparm(:) = 0
! Defaults first
p%baseprecv(ilev_)%iprcparm(p_type_) = asm_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_halo_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 1
p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
if (present(iv)) then
isz = size(iv)
if (isz >= 1) p%baseprecv(ilev_)%iprcparm(n_ovr_) = iv(1)
if (isz >= 2) p%baseprecv(ilev_)%iprcparm(restr_) = iv(2)
if (isz >= 3) p%baseprecv(ilev_)%iprcparm(prol_) = iv(3)
if (isz >= 4) p%baseprecv(ilev_)%iprcparm(f_type_) = iv(4)
if (isz >= 5) p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = iv(5)
! Do not consider renum for the time being.
!!$ if (isz >= 5) p%baseprecv(ilev_)%iprcparm(iren_) = iv(5)
endif endif
end subroutine psb_zprecseti
subroutine psb_zprecsetd(p,what,val,info,ilev)
use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_zprecsetd
implicit none
type(psb_zprec_type), intent(inout) :: p
integer, intent(in) :: what
real(kind(1.d0)), intent(in) :: val
integer, intent(out) :: info
integer, optional, intent(in) :: ilev
integer :: ilev_,nlev_
info = 0
case ('ML', '2L', '2LEV') if (present(ilev)) then
ilev_ = ilev
else
p%baseprecv(ilev_)%iprcparm(:) = 0 ilev_ = 1
p%baseprecv(ilev_)%iprcparm(p_type_) = bjac_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(ml_type_) = mult_ml_prec_
p%baseprecv(ilev_)%iprcparm(aggr_alg_) = loc_aggr_
p%baseprecv(ilev_)%iprcparm(smth_kind_) = smth_omg_
p%baseprecv(ilev_)%iprcparm(coarse_mat_) = mat_distr_
p%baseprecv(ilev_)%iprcparm(smth_pos_) = post_smooth_
p%baseprecv(ilev_)%iprcparm(glb_smth_) = 1
p%baseprecv(ilev_)%iprcparm(om_choice_) = lib_choice_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(ilev_)%dprcparm(smooth_omega_) = 4.d0/3.d0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
if (present(iv)) then
isz = size(iv)
if (isz >= 1) p%baseprecv(ilev_)%iprcparm(ml_type_) = iv(1)
if (isz >= 2) p%baseprecv(ilev_)%iprcparm(aggr_alg_) = iv(2)
if (isz >= 3) p%baseprecv(ilev_)%iprcparm(coarse_mat_) = iv(3)
if (isz >= 4) p%baseprecv(ilev_)%iprcparm(smth_pos_) = iv(4)
if (isz >= 5) p%baseprecv(ilev_)%iprcparm(f_type_) = iv(5)
if (isz >= 6) p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = iv(6)
if (isz >= 7) p%baseprecv(ilev_)%iprcparm(smth_kind_) = iv(7)
end if end if
if (present(rs)) then if (.not.allocated(p%baseprecv)) then
p%baseprecv(ilev_)%iprcparm(om_choice_) = user_choice_ write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
p%baseprecv(ilev_)%dprcparm(smooth_omega_) = rs info = -1
return
endif endif
nlev_ = size(p%baseprecv)
if ((ilev_<1).or.(ilev_ > nlev_)) then
write(0,*) 'PRECSET ERRROR: ilev out of bounds'
info = -1
return
endif
if (.not.allocated(p%baseprecv(ilev_)%dprcparm)) then
write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
info = -1
return
endif
if (ilev_ == 1) then
! Rules for fine level are slightly different.
select case(what)
! Right now we don't have any at base level. Will change when
! we implement F_ILU_E_
case default case default
write(0,*) 'Unknown preconditioner type request "',ptype,'"' write(0,*) 'Error: trying to call PRECSET with an invalid WHAT'
err = 2 info = -2
end select end select
else if (ilev_ > 1) then
select case(what)
case(smooth_omega_)
p%baseprecv(ilev_)%dprcparm(what) = val
case default
write(0,*) 'Error: trying to call PRECSET with an invalid WHAT'
info = -2
end select
endif
info = err end subroutine psb_zprecsetd
end subroutine psb_zprecset

Loading…
Cancel
Save