fdurastante
/
amg4psblas
mirror of https://github.com/sfilippone/amg4psblas.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

removed set moved to psblas

Browse Source
merge-amgext
Cirdans-Home 4 years ago
parent 32792507f5
commit b5cbfd5356
4 changed files with 584 additions and 584 deletions
Show Stats Download Patch File Download Diff File

292
amgprec/amg_c_prec_type.f90
Unescape Escape View File

@ -1,15 +1,15 @@
! !
! !
! AMG4PSBLAS version 1.0 ! AMG4PSBLAS version 1.0
! Algebraic Multigrid Package ! Algebraic Multigrid Package
! based on PSBLAS (Parallel Sparse BLAS version 3.5) ! based on PSBLAS (Parallel Sparse BLAS version 3.5)
! !
! (C) Copyright 2020 ! (C) Copyright 2020
! !
! Salvatore Filippone ! Salvatore Filippone
! Pasqua D'Ambra ! Pasqua D'Ambra
! Fabio Durastante ! Fabio Durastante
! !
! Redistribution and use in source and binary forms, with or without ! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions ! modification, are permitted provided that the following conditions
! are met: ! are met:
@ -21,7 +21,7 @@
! 3. The name of the AMG4PSBLAS group or the names of its contributors may ! 3. The name of the AMG4PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this ! not be used to endorse or promote products derived from this
! software without specific written permission. ! software without specific written permission.
! !
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
@ -33,21 +33,21 @@
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! 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 ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE. ! POSSIBILITY OF SUCH DAMAGE.
! !
! !
! File: amg_c_prec_type.f90 ! File: amg_c_prec_type.f90
! !
! Module: amg_c_prec_type ! Module: amg_c_prec_type
! !
! This module defines: ! This module defines:
! - the amg_c_prec_type data structure containing the preconditioner and related ! - the amg_c_prec_type data structure containing the preconditioner and related
! data structures; ! data structures;
! !
! It contains routines for ! It contains routines for
! - Building and applying; ! - Building and applying;
! - checking if the preconditioner is correctly defined; ! - checking if the preconditioner is correctly defined;
! - printing a description of the preconditioner; ! - printing a description of the preconditioner;
! - deallocating the preconditioner data structure. ! - deallocating the preconditioner data structure.
! !
module amg_c_prec_type module amg_c_prec_type
@ -70,25 +70,25 @@ module amg_c_prec_type
! It consists of an array of 'one-level' intermediate data structures ! It consists of an array of 'one-level' intermediate data structures
! of type amg_conelev_type, each containing the information needed to apply ! of type amg_conelev_type, each containing the information needed to apply
! the smoothing and the coarse-space correction at a generic level. RT is the ! the smoothing and the coarse-space correction at a generic level. RT is the
! real data type, i.e. S for both S and C, and D for both D and Z. ! real data type, i.e. S for both S and C, and D for both D and Z.
! !
! type amg_cprec_type ! type amg_cprec_type
! type(amg_conelev_type), allocatable :: precv(:) ! type(amg_conelev_type), allocatable :: precv(:)
! end type amg_cprec_type ! end type amg_cprec_type
! !
! Note that the levels are numbered in increasing order starting from ! Note that the levels are numbered in increasing order starting from
! the level 1 as the finest one, and the number of levels is given by ! the level 1 as the finest one, and the number of levels is given by
! size(precv(:)) which is the id of the coarsest level. ! size(precv(:)) which is the id of the coarsest level.
! In the multigrid literature many authors number the levels in the opposite ! In the multigrid literature many authors number the levels in the opposite
! order, with level 0 being the id of the coarsest level. ! order, with level 0 being the id of the coarsest level.
! !
! !
integer, parameter, private :: wv_size_=4 integer, parameter, private :: wv_size_=4
type, extends(psb_cprec_type) :: amg_cprec_type type, extends(psb_cprec_type) :: amg_cprec_type
type(amg_saggr_data) :: ag_data type(amg_saggr_data) :: ag_data
! !
! Number of outer sweeps. Sometimes 2 V-cycles may be better than 1 W-cycle. ! Number of outer sweeps. Sometimes 2 V-cycles may be better than 1 W-cycle.
! !
integer(psb_ipk_) :: outer_sweeps = 1 integer(psb_ipk_) :: outer_sweeps = 1
! !
@ -97,11 +97,11 @@ module amg_c_prec_type
! to keep track against what is put later in the multilevel array ! to keep track against what is put later in the multilevel array
! !
integer(psb_ipk_) :: coarse_solver = -1 integer(psb_ipk_) :: coarse_solver = -1
! !
! The multilevel hierarchy ! The multilevel hierarchy
! !
type(amg_c_onelev_type), allocatable :: precv(:) type(amg_c_onelev_type), allocatable :: precv(:)
contains contains
procedure, pass(prec) :: psb_c_apply2_vect => amg_c_apply2_vect procedure, pass(prec) :: psb_c_apply2_vect => amg_c_apply2_vect
procedure, pass(prec) :: psb_c_apply1_vect => amg_c_apply1_vect procedure, pass(prec) :: psb_c_apply1_vect => amg_c_apply1_vect
@ -127,7 +127,7 @@ module amg_c_prec_type
procedure, pass(prec) :: cseti => amg_ccprecseti procedure, pass(prec) :: cseti => amg_ccprecseti
procedure, pass(prec) :: csetc => amg_ccprecsetc procedure, pass(prec) :: csetc => amg_ccprecsetc
procedure, pass(prec) :: csetr => amg_ccprecsetr procedure, pass(prec) :: csetr => amg_ccprecsetr
generic, public :: set => cseti, csetc, csetr, setsm, setsv, setag generic, public :: set => setsm, setsv, setag
procedure, pass(prec) :: get_smoother => amg_c_get_smootherp procedure, pass(prec) :: get_smoother => amg_c_get_smootherp
procedure, pass(prec) :: get_solver => amg_c_get_solverp procedure, pass(prec) :: get_solver => amg_c_get_solverp
procedure, pass(prec) :: move_alloc => c_prec_move_alloc procedure, pass(prec) :: move_alloc => c_prec_move_alloc
@ -157,7 +157,7 @@ module amg_c_prec_type
interface amg_precdescr interface amg_precdescr
subroutine amg_cfile_prec_descr(prec,iout,root) subroutine amg_cfile_prec_descr(prec,iout,root)
import :: amg_cprec_type, psb_ipk_ import :: amg_cprec_type, psb_ipk_
implicit none implicit none
! Arguments ! Arguments
class(amg_cprec_type), intent(in) :: prec class(amg_cprec_type), intent(in) :: prec
integer(psb_ipk_), intent(in), optional :: iout integer(psb_ipk_), intent(in), optional :: iout
@ -211,7 +211,7 @@ module amg_c_prec_type
end subroutine amg_cprecaply1 end subroutine amg_cprecaply1
end interface end interface
interface interface
subroutine amg_cprecsetsm(prec,val,info,ilev,ilmax,pos) subroutine amg_cprecsetsm(prec,val,info,ilev,ilmax,pos)
import :: psb_cspmat_type, psb_desc_type, psb_spk_, & import :: psb_cspmat_type, psb_desc_type, psb_spk_, &
& amg_cprec_type, amg_c_base_smoother_type, psb_ipk_ & amg_cprec_type, amg_c_base_smoother_type, psb_ipk_
@ -243,7 +243,7 @@ module amg_c_prec_type
import :: psb_cspmat_type, psb_desc_type, psb_spk_, & import :: psb_cspmat_type, psb_desc_type, psb_spk_, &
& amg_cprec_type, psb_ipk_ & amg_cprec_type, psb_ipk_
class(amg_cprec_type), intent(inout) :: prec class(amg_cprec_type), intent(inout) :: prec
character(len=*), intent(in) :: what character(len=*), intent(in) :: what
integer(psb_ipk_), intent(in) :: val integer(psb_ipk_), intent(in) :: val
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx
@ -253,7 +253,7 @@ module amg_c_prec_type
import :: psb_cspmat_type, psb_desc_type, psb_spk_, & import :: psb_cspmat_type, psb_desc_type, psb_spk_, &
& amg_cprec_type, psb_ipk_ & amg_cprec_type, psb_ipk_
class(amg_cprec_type), intent(inout) :: prec class(amg_cprec_type), intent(inout) :: prec
character(len=*), intent(in) :: what character(len=*), intent(in) :: what
real(psb_spk_), intent(in) :: val real(psb_spk_), intent(in) :: val
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx
@ -263,7 +263,7 @@ module amg_c_prec_type
import :: psb_cspmat_type, psb_desc_type, psb_spk_, & import :: psb_cspmat_type, psb_desc_type, psb_spk_, &
& amg_cprec_type, psb_ipk_ & amg_cprec_type, psb_ipk_
class(amg_cprec_type), intent(inout) :: prec class(amg_cprec_type), intent(inout) :: prec
character(len=*), intent(in) :: what character(len=*), intent(in) :: what
character(len=*), intent(in) :: string character(len=*), intent(in) :: string
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx
@ -341,28 +341,28 @@ module amg_c_prec_type
! character, intent(in),optional :: upd ! character, intent(in),optional :: upd
end subroutine amg_c_smoothers_bld end subroutine amg_c_smoothers_bld
end interface amg_smoothers_bld end interface amg_smoothers_bld
contains contains
! !
! Function returning a pointer to the smoother ! Function returning a pointer to the smoother
! !
function amg_c_get_smootherp(prec,ilev) result(val) function amg_c_get_smootherp(prec,ilev) result(val)
implicit none implicit none
class(amg_cprec_type), target, intent(in) :: prec class(amg_cprec_type), target, intent(in) :: prec
integer(psb_ipk_), optional :: ilev integer(psb_ipk_), optional :: ilev
class(amg_c_base_smoother_type), pointer :: val class(amg_c_base_smoother_type), pointer :: val
integer(psb_ipk_) :: ilev_ integer(psb_ipk_) :: ilev_
val => null() val => null()
if (present(ilev)) then if (present(ilev)) then
ilev_ = ilev ilev_ = ilev
else else
! What is a good default? ! What is a good default?
ilev_ = 1 ilev_ = 1
end if end if
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
if ((1<=ilev_).and.(ilev_<=size(prec%precv))) then if ((1<=ilev_).and.(ilev_<=size(prec%precv))) then
if (allocated(prec%precv(ilev_)%sm)) then if (allocated(prec%precv(ilev_)%sm)) then
val => prec%precv(ilev_)%sm val => prec%precv(ilev_)%sm
end if end if
end if end if
@ -372,23 +372,23 @@ contains
! Function returning a pointer to the solver ! Function returning a pointer to the solver
! !
function amg_c_get_solverp(prec,ilev) result(val) function amg_c_get_solverp(prec,ilev) result(val)
implicit none implicit none
class(amg_cprec_type), target, intent(in) :: prec class(amg_cprec_type), target, intent(in) :: prec
integer(psb_ipk_), optional :: ilev integer(psb_ipk_), optional :: ilev
class(amg_c_base_solver_type), pointer :: val class(amg_c_base_solver_type), pointer :: val
integer(psb_ipk_) :: ilev_ integer(psb_ipk_) :: ilev_
val => null() val => null()
if (present(ilev)) then if (present(ilev)) then
ilev_ = ilev ilev_ = ilev
else else
! What is a good default? ! What is a good default?
ilev_ = 1 ilev_ = 1
end if end if
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
if ((1<=ilev_).and.(ilev_<=size(prec%precv))) then if ((1<=ilev_).and.(ilev_<=size(prec%precv))) then
if (allocated(prec%precv(ilev_)%sm)) then if (allocated(prec%precv(ilev_)%sm)) then
if (allocated(prec%precv(ilev_)%sm%sv)) then if (allocated(prec%precv(ilev_)%sm%sv)) then
val => prec%precv(ilev_)%sm%sv val => prec%precv(ilev_)%sm%sv
end if end if
end if end if
@ -399,25 +399,25 @@ contains
! Function returning the size of the precv(:) array ! Function returning the size of the precv(:) array
! !
function amg_c_get_nlevs(prec) result(val) function amg_c_get_nlevs(prec) result(val)
implicit none implicit none
class(amg_cprec_type), intent(in) :: prec class(amg_cprec_type), intent(in) :: prec
integer(psb_ipk_) :: val integer(psb_ipk_) :: val
val = 0 val = 0
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
val = size(prec%precv) val = size(prec%precv)
end if end if
end function amg_c_get_nlevs end function amg_c_get_nlevs
! !
! Function returning the size of the amg_prec_type data structure ! Function returning the size of the amg_prec_type data structure
! in bytes or in number of nonzeros of the operator(s) involved. ! in bytes or in number of nonzeros of the operator(s) involved.
! !
function amg_c_get_nzeros(prec) result(val) function amg_c_get_nzeros(prec) result(val)
implicit none implicit none
class(amg_cprec_type), intent(in) :: prec class(amg_cprec_type), intent(in) :: prec
integer(psb_epk_) :: val integer(psb_epk_) :: val
integer(psb_ipk_) :: i integer(psb_ipk_) :: i
val = 0 val = 0
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
do i=1, size(prec%precv) do i=1, size(prec%precv)
val = val + prec%precv(i)%get_nzeros() val = val + prec%precv(i)%get_nzeros()
end do end do
@ -425,13 +425,13 @@ contains
end function amg_c_get_nzeros end function amg_c_get_nzeros
function amg_cprec_sizeof(prec) result(val) function amg_cprec_sizeof(prec) result(val)
implicit none implicit none
class(amg_cprec_type), intent(in) :: prec class(amg_cprec_type), intent(in) :: prec
integer(psb_epk_) :: val integer(psb_epk_) :: val
integer(psb_ipk_) :: i integer(psb_ipk_) :: i
val = 0 val = 0
val = val + psb_sizeof_ip val = val + psb_sizeof_ip
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
do i=1, size(prec%precv) do i=1, size(prec%precv)
val = val + prec%precv(i)%sizeof() val = val + prec%precv(i)%sizeof()
end do end do
@ -444,40 +444,40 @@ contains
! various level to the nonzeroes at the fine level ! various level to the nonzeroes at the fine level
! (original matrix) ! (original matrix)
! !
function amg_c_get_compl(prec) result(val) function amg_c_get_compl(prec) result(val)
implicit none implicit none
class(amg_cprec_type), intent(in) :: prec class(amg_cprec_type), intent(in) :: prec
complex(psb_spk_) :: val complex(psb_spk_) :: val
val = prec%ag_data%op_complexity val = prec%ag_data%op_complexity
end function amg_c_get_compl end function amg_c_get_compl
subroutine amg_c_cmp_compl(prec)
implicit none subroutine amg_c_cmp_compl(prec)
implicit none
class(amg_cprec_type), intent(inout) :: prec class(amg_cprec_type), intent(inout) :: prec
real(psb_spk_) :: num, den, nmin real(psb_spk_) :: num, den, nmin
type(psb_ctxt_type) :: ctxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: il integer(psb_ipk_) :: il
num = -sone num = -sone
den = sone den = sone
ctxt = prec%ctxt ctxt = prec%ctxt
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
il = 1 il = 1
num = prec%precv(il)%base_a%get_nzeros() num = prec%precv(il)%base_a%get_nzeros()
if (num >= szero) then if (num >= szero) then
den = num den = num
do il=2,size(prec%precv) do il=2,size(prec%precv)
num = num + max(0,prec%precv(il)%base_a%get_nzeros()) num = num + max(0,prec%precv(il)%base_a%get_nzeros())
end do end do
end if end if
end if end if
nmin = num nmin = num
call psb_min(ctxt,nmin) call psb_min(ctxt,nmin)
if (nmin < szero) then if (nmin < szero) then
num = szero num = szero
den = sone den = sone
@ -487,25 +487,25 @@ contains
end if end if
prec%ag_data%op_complexity = num/den prec%ag_data%op_complexity = num/den
end subroutine amg_c_cmp_compl end subroutine amg_c_cmp_compl
! !
! Average coarsening ratio ! Average coarsening ratio
! !
function amg_c_get_avg_cr(prec) result(val) function amg_c_get_avg_cr(prec) result(val)
implicit none implicit none
class(amg_cprec_type), intent(in) :: prec class(amg_cprec_type), intent(in) :: prec
complex(psb_spk_) :: val complex(psb_spk_) :: val
val = prec%ag_data%avg_cr val = prec%ag_data%avg_cr
end function amg_c_get_avg_cr end function amg_c_get_avg_cr
subroutine amg_c_cmp_avg_cr(prec)
implicit none subroutine amg_c_cmp_avg_cr(prec)
implicit none
class(amg_cprec_type), intent(inout) :: prec class(amg_cprec_type), intent(inout) :: prec
real(psb_spk_) :: avgcr real(psb_spk_) :: avgcr
type(psb_ctxt_type) :: ctxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: il, nl, iam, np integer(psb_ipk_) :: il, nl, iam, np
@ -519,12 +519,12 @@ contains
do il=2,nl do il=2,nl
avgcr = avgcr + max(szero,prec%precv(il)%szratio) avgcr = avgcr + max(szero,prec%precv(il)%szratio)
end do end do
avgcr = avgcr / (nl-1) avgcr = avgcr / (nl-1)
end if end if
call psb_sum(ctxt,avgcr) call psb_sum(ctxt,avgcr)
prec%ag_data%avg_cr = avgcr/np prec%ag_data%avg_cr = avgcr/np
end subroutine amg_c_cmp_avg_cr end subroutine amg_c_cmp_avg_cr
! !
! Subroutines: amg_Tprec_free ! Subroutines: amg_Tprec_free
! Version: complex ! Version: complex
@ -538,74 +538,74 @@ contains
! error code. ! error code.
! !
subroutine amg_cprecfree(p,info) subroutine amg_cprecfree(p,info)
implicit none implicit none
! Arguments ! Arguments
type(amg_cprec_type), intent(inout) :: p type(amg_cprec_type), intent(inout) :: p
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
! Local variables ! Local variables
integer(psb_ipk_) :: me,err_act,i integer(psb_ipk_) :: me,err_act,i
character(len=20) :: name character(len=20) :: name
info=psb_success_ info=psb_success_
name = 'amg_cprecfree' name = 'amg_cprecfree'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; return info = psb_err_internal_error_; return
end if end if
me=-1 me=-1
call p%free(info) call p%free(info)
return return
end subroutine amg_cprecfree end subroutine amg_cprecfree
subroutine amg_c_prec_free(prec,info) subroutine amg_c_prec_free(prec,info)
implicit none implicit none
! Arguments ! Arguments
class(amg_cprec_type), intent(inout) :: prec class(amg_cprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
! Local variables ! Local variables
integer(psb_ipk_) :: me,err_act,i integer(psb_ipk_) :: me,err_act,i
character(len=20) :: name character(len=20) :: name
info=psb_success_ info=psb_success_
name = 'amg_cprecfree' name = 'amg_cprecfree'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999 info = psb_err_internal_error_; goto 9999
end if end if
me=-1 me=-1
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
do i=1,size(prec%precv) do i=1,size(prec%precv)
call prec%precv(i)%free(info) call prec%precv(i)%free(info)
end do end do
deallocate(prec%precv,stat=info) deallocate(prec%precv,stat=info)
end if end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
9999 call psb_error_handler(err_act) 9999 call psb_error_handler(err_act)
return return
end subroutine amg_c_prec_free end subroutine amg_c_prec_free
! !
! Top level methods. ! Top level methods.
! !
subroutine amg_c_apply2_vect(prec,x,y,desc_data,info,trans,work) subroutine amg_c_apply2_vect(prec,x,y,desc_data,info,trans,work)
implicit none implicit none
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
class(amg_cprec_type), intent(inout) :: prec class(amg_cprec_type), intent(inout) :: prec
type(psb_c_vect_type),intent(inout) :: x type(psb_c_vect_type),intent(inout) :: x
@ -618,13 +618,13 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
select type(prec) select type(prec)
type is (amg_cprec_type) type is (amg_cprec_type)
call amg_precapply(prec,x,y,desc_data,info,trans,work) call amg_precapply(prec,x,y,desc_data,info,trans,work)
class default class default
info = psb_err_missing_override_method_ info = psb_err_missing_override_method_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end select end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -636,7 +636,7 @@ contains
end subroutine amg_c_apply2_vect end subroutine amg_c_apply2_vect
subroutine amg_c_apply1_vect(prec,x,desc_data,info,trans,work) subroutine amg_c_apply1_vect(prec,x,desc_data,info,trans,work)
implicit none implicit none
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
class(amg_cprec_type), intent(inout) :: prec class(amg_cprec_type), intent(inout) :: prec
type(psb_c_vect_type),intent(inout) :: x type(psb_c_vect_type),intent(inout) :: x
@ -648,13 +648,13 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
select type(prec) select type(prec)
type is (amg_cprec_type) type is (amg_cprec_type)
call amg_precapply(prec,x,desc_data,info,trans,work) call amg_precapply(prec,x,desc_data,info,trans,work)
class default class default
info = psb_err_missing_override_method_ info = psb_err_missing_override_method_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end select end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -667,7 +667,7 @@ contains
subroutine amg_c_apply2v(prec,x,y,desc_data,info,trans,work) subroutine amg_c_apply2v(prec,x,y,desc_data,info,trans,work)
implicit none implicit none
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
class(amg_cprec_type), intent(inout) :: prec class(amg_cprec_type), intent(inout) :: prec
complex(psb_spk_),intent(inout) :: x(:) complex(psb_spk_),intent(inout) :: x(:)
@ -680,13 +680,13 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
select type(prec) select type(prec)
type is (amg_cprec_type) type is (amg_cprec_type)
call amg_precapply(prec,x,y,desc_data,info,trans,work) call amg_precapply(prec,x,y,desc_data,info,trans,work)
class default class default
info = psb_err_missing_override_method_ info = psb_err_missing_override_method_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end select end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -698,7 +698,7 @@ contains
end subroutine amg_c_apply2v end subroutine amg_c_apply2v
subroutine amg_c_apply1v(prec,x,desc_data,info,trans) subroutine amg_c_apply1v(prec,x,desc_data,info,trans)
implicit none implicit none
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
class(amg_cprec_type), intent(inout) :: prec class(amg_cprec_type), intent(inout) :: prec
complex(psb_spk_),intent(inout) :: x(:) complex(psb_spk_),intent(inout) :: x(:)
@ -709,13 +709,13 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
select type(prec) select type(prec)
type is (amg_cprec_type) type is (amg_cprec_type)
call amg_precapply(prec,x,desc_data,info,trans) call amg_precapply(prec,x,desc_data,info,trans)
class default class default
info = psb_err_missing_override_method_ info = psb_err_missing_override_method_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end select end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -730,8 +730,8 @@ contains
subroutine amg_c_dump(prec,info,istart,iend,iproc,prefix,head,& subroutine amg_c_dump(prec,info,istart,iend,iproc,prefix,head,&
& ac,rp,smoother,solver,tprol,& & ac,rp,smoother,solver,tprol,&
& global_num) & global_num)
implicit none implicit none
class(amg_cprec_type), intent(in) :: prec class(amg_cprec_type), intent(in) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: istart, iend, iproc integer(psb_ipk_), intent(in), optional :: istart, iend, iproc
@ -742,27 +742,27 @@ contains
integer(psb_ipk_) :: iam, np, iproc_ integer(psb_ipk_) :: iam, np, iproc_
character(len=80) :: prefix_ character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than character(len=120) :: fname ! len should be at least 20 more than
! len of prefix_ ! len of prefix_
info = 0 info = 0
icontxt = prec%ctxt icontxt = prec%ctxt
call psb_info(icontxt,iam,np) call psb_info(icontxt,iam,np)
iln = size(prec%precv) iln = size(prec%precv)
if (present(istart)) then if (present(istart)) then
il1 = max(1,istart) il1 = max(1,istart)
else else
il1 = min(2,iln) il1 = min(2,iln)
end if end if
if (present(iend)) then if (present(iend)) then
iln = min(iln, iend) iln = min(iln, iend)
end if end if
iproc_ = -1 iproc_ = -1
if (present(iproc)) then if (present(iproc)) then
iproc_ = iproc iproc_ = iproc
end if end if
if ((iproc_ == -1).or.(iproc_==iam)) then if ((iproc_ == -1).or.(iproc_==iam)) then
do lev=il1, iln do lev=il1, iln
call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,& call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,&
& ac=ac,smoother=smoother,solver=solver,rp=rp,tprol=tprol, & & ac=ac,smoother=smoother,solver=solver,rp=rp,tprol=tprol, &
@ -773,7 +773,7 @@ contains
subroutine amg_c_cnv(prec,info,amold,vmold,imold) subroutine amg_c_cnv(prec,info,amold,vmold,imold)
implicit none implicit none
class(amg_cprec_type), intent(inout) :: prec class(amg_cprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
class(psb_c_base_sparse_mat), intent(in), optional :: amold class(psb_c_base_sparse_mat), intent(in), optional :: amold
@ -781,7 +781,7 @@ contains
class(psb_i_base_vect_type), intent(in), optional :: imold class(psb_i_base_vect_type), intent(in), optional :: imold
integer(psb_ipk_) :: i integer(psb_ipk_) :: i
info = psb_success_ info = psb_success_
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
do i=1,size(prec%precv) do i=1,size(prec%precv)
@ -789,24 +789,24 @@ contains
& call prec%precv(i)%cnv(info,amold=amold,vmold=vmold,imold=imold) & call prec%precv(i)%cnv(info,amold=amold,vmold=vmold,imold=imold)
end do end do
end if end if
end subroutine amg_c_cnv end subroutine amg_c_cnv
subroutine amg_c_clone(prec,precout,info) subroutine amg_c_clone(prec,precout,info)
implicit none implicit none
class(amg_cprec_type), intent(inout) :: prec class(amg_cprec_type), intent(inout) :: prec
class(psb_cprec_type), intent(inout) :: precout class(psb_cprec_type), intent(inout) :: precout
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
call precout%free(info) call precout%free(info)
if (info == 0) call amg_c_inner_clone(prec,precout,info) if (info == 0) call amg_c_inner_clone(prec,precout,info)
end subroutine amg_c_clone end subroutine amg_c_clone
subroutine amg_c_inner_clone(prec,precout,info) subroutine amg_c_inner_clone(prec,precout,info)
implicit none implicit none
class(amg_cprec_type), intent(inout) :: prec class(amg_cprec_type), intent(inout) :: prec
class(psb_cprec_type), target, intent(inout) :: precout class(psb_cprec_type), target, intent(inout) :: precout
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -821,17 +821,17 @@ contains
pout%ctxt = prec%ctxt pout%ctxt = prec%ctxt
pout%ag_data = prec%ag_data pout%ag_data = prec%ag_data
pout%outer_sweeps = prec%outer_sweeps pout%outer_sweeps = prec%outer_sweeps
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
ln = size(prec%precv) ln = size(prec%precv)
allocate(pout%precv(ln),stat=info) allocate(pout%precv(ln),stat=info)
if (info /= psb_success_) goto 9999 if (info /= psb_success_) goto 9999
if (ln >= 1) then if (ln >= 1) then
call prec%precv(1)%clone(pout%precv(1),info) call prec%precv(1)%clone(pout%precv(1),info)
end if end if
do lev=2, ln do lev=2, ln
if (info /= psb_success_) exit if (info /= psb_success_) exit
call prec%precv(lev)%clone(pout%precv(lev),info) call prec%precv(lev)%clone(pout%precv(lev),info)
if (info == psb_success_) then if (info == psb_success_) then
pout%precv(lev)%base_a => pout%precv(lev)%ac pout%precv(lev)%base_a => pout%precv(lev)%ac
pout%precv(lev)%base_desc => pout%precv(lev)%desc_ac pout%precv(lev)%base_desc => pout%precv(lev)%desc_ac
pout%precv(lev)%map%p_desc_U => pout%precv(lev-1)%base_desc pout%precv(lev)%map%p_desc_U => pout%precv(lev-1)%base_desc
@ -842,7 +842,7 @@ contains
if (allocated(prec%precv(1)%wrk)) & if (allocated(prec%precv(1)%wrk)) &
& call pout%allocate_wrk(info,vmold=prec%precv(1)%wrk%vx2l%v) & call pout%allocate_wrk(info,vmold=prec%precv(1)%wrk%vx2l%v)
class default class default
write(0,*) 'Error: wrong out type' write(0,*) 'Error: wrong out type'
info = psb_err_invalid_input_ info = psb_err_invalid_input_
end select end select
@ -854,14 +854,14 @@ contains
implicit none implicit none
class(amg_cprec_type), intent(inout) :: prec class(amg_cprec_type), intent(inout) :: prec
class(amg_cprec_type), intent(inout), target :: b class(amg_cprec_type), intent(inout), target :: b
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i integer(psb_ipk_) :: i
if (same_type_as(prec,b)) then if (same_type_as(prec,b)) then
if (allocated(b%precv)) then if (allocated(b%precv)) then
! This might not be required if FINAL procedures are available. ! This might not be required if FINAL procedures are available.
call b%free(info) call b%free(info)
if (info /= psb_success_) then if (info /= psb_success_) then
!????? !?????
!!$ return !!$ return
endif endif
@ -869,7 +869,7 @@ contains
b%ctxt = prec%ctxt b%ctxt = prec%ctxt
b%ag_data = prec%ag_data b%ag_data = prec%ag_data
b%outer_sweeps = prec%outer_sweeps b%outer_sweeps = prec%outer_sweeps
call move_alloc(prec%precv,b%precv) call move_alloc(prec%precv,b%precv)
! Fix the pointers except on level 1. ! Fix the pointers except on level 1.
do i=2, size(b%precv) do i=2, size(b%precv)
@ -878,7 +878,7 @@ contains
b%precv(i)%map%p_desc_U => b%precv(i-1)%base_desc b%precv(i)%map%p_desc_U => b%precv(i-1)%base_desc
b%precv(i)%map%p_desc_V => b%precv(i)%base_desc b%precv(i)%map%p_desc_V => b%precv(i)%base_desc
end do end do
else else
write(0,*) 'Warning: PREC%move_alloc onto different type?' write(0,*) 'Warning: PREC%move_alloc onto different type?'
info = psb_err_internal_error_ info = psb_err_internal_error_
@ -888,7 +888,7 @@ contains
subroutine amg_c_allocate_wrk(prec,info,vmold,desc) subroutine amg_c_allocate_wrk(prec,info,vmold,desc)
use psb_base_mod use psb_base_mod
implicit none implicit none
! Arguments ! Arguments
class(amg_cprec_type), intent(inout) :: prec class(amg_cprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -896,37 +896,37 @@ contains
! !
! In MLD the DESC optional argument is ignored, since ! In MLD the DESC optional argument is ignored, since
! the necessary info is contained in the various entries of the ! the necessary info is contained in the various entries of the
! PRECV component. ! PRECV component.
type(psb_desc_type), intent(in), optional :: desc type(psb_desc_type), intent(in), optional :: desc
! Local variables ! Local variables
integer(psb_ipk_) :: me,err_act,i,j,level,nlev, nc2l integer(psb_ipk_) :: me,err_act,i,j,level,nlev, nc2l
character(len=20) :: name character(len=20) :: name
info=psb_success_ info=psb_success_
name = 'amg_c_allocate_wrk' name = 'amg_c_allocate_wrk'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999 info = psb_err_internal_error_; goto 9999
end if end if
nlev = size(prec%precv) nlev = size(prec%precv)
level = 1 level = 1
do level = 1, nlev do level = 1, nlev
call prec%precv(level)%allocate_wrk(info,vmold=vmold) call prec%precv(level)%allocate_wrk(info,vmold=vmold)
if (psb_errstatus_fatal()) then if (psb_errstatus_fatal()) then
nc2l = prec%precv(level)%base_desc%get_local_cols() nc2l = prec%precv(level)%base_desc%get_local_cols()
info=psb_err_alloc_request_ info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l/), a_err='complex(psb_spk_)') call psb_errpush(info,name,i_err=(/2*nc2l/), a_err='complex(psb_spk_)')
goto 9999 goto 9999
end if end if
end do end do
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
9999 call psb_error_handler(err_act) 9999 call psb_error_handler(err_act)
return return
end subroutine amg_c_allocate_wrk end subroutine amg_c_allocate_wrk
subroutine amg_c_free_wrk(prec,info) subroutine amg_c_free_wrk(prec,info)
@ -948,13 +948,13 @@ contains
info = psb_err_internal_error_; goto 9999 info = psb_err_internal_error_; goto 9999
end if end if
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
nlev = size(prec%precv) nlev = size(prec%precv)
do level = 1, nlev do level = 1, nlev
call prec%precv(level)%free_wrk(info) call prec%precv(level)%free_wrk(info)
end do end do
end if end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -966,7 +966,7 @@ contains
function amg_c_is_allocated_wrk(prec) result(res) function amg_c_is_allocated_wrk(prec) result(res)
use psb_base_mod use psb_base_mod
implicit none implicit none
! Arguments ! Arguments
class(amg_cprec_type), intent(in) :: prec class(amg_cprec_type), intent(in) :: prec
logical :: res logical :: res
@ -974,7 +974,7 @@ contains
res = .false. res = .false.
if (.not.allocated(prec%precv)) return if (.not.allocated(prec%precv)) return
res = allocated(prec%precv(1)%wrk) res = allocated(prec%precv(1)%wrk)
end function amg_c_is_allocated_wrk end function amg_c_is_allocated_wrk
end module amg_c_prec_type end module amg_c_prec_type

292
amgprec/amg_d_prec_type.f90
Unescape Escape View File

@ -1,15 +1,15 @@
! !
! !
! AMG4PSBLAS version 1.0 ! AMG4PSBLAS version 1.0
! Algebraic Multigrid Package ! Algebraic Multigrid Package
! based on PSBLAS (Parallel Sparse BLAS version 3.5) ! based on PSBLAS (Parallel Sparse BLAS version 3.5)
! !
! (C) Copyright 2020 ! (C) Copyright 2020
! !
! Salvatore Filippone ! Salvatore Filippone
! Pasqua D'Ambra ! Pasqua D'Ambra
! Fabio Durastante ! Fabio Durastante
! !
! Redistribution and use in source and binary forms, with or without ! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions ! modification, are permitted provided that the following conditions
! are met: ! are met:
@ -21,7 +21,7 @@
! 3. The name of the AMG4PSBLAS group or the names of its contributors may ! 3. The name of the AMG4PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this ! not be used to endorse or promote products derived from this
! software without specific written permission. ! software without specific written permission.
! !
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
@ -33,21 +33,21 @@
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! 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 ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE. ! POSSIBILITY OF SUCH DAMAGE.
! !
! !
! File: amg_d_prec_type.f90 ! File: amg_d_prec_type.f90
! !
! Module: amg_d_prec_type ! Module: amg_d_prec_type
! !
! This module defines: ! This module defines:
! - the amg_d_prec_type data structure containing the preconditioner and related ! - the amg_d_prec_type data structure containing the preconditioner and related
! data structures; ! data structures;
! !
! It contains routines for ! It contains routines for
! - Building and applying; ! - Building and applying;
! - checking if the preconditioner is correctly defined; ! - checking if the preconditioner is correctly defined;
! - printing a description of the preconditioner; ! - printing a description of the preconditioner;
! - deallocating the preconditioner data structure. ! - deallocating the preconditioner data structure.
! !
module amg_d_prec_type module amg_d_prec_type
@ -70,25 +70,25 @@ module amg_d_prec_type
! It consists of an array of 'one-level' intermediate data structures ! It consists of an array of 'one-level' intermediate data structures
! of type amg_donelev_type, each containing the information needed to apply ! of type amg_donelev_type, each containing the information needed to apply
! the smoothing and the coarse-space correction at a generic level. RT is the ! the smoothing and the coarse-space correction at a generic level. RT is the
! real data type, i.e. S for both S and C, and D for both D and Z. ! real data type, i.e. S for both S and C, and D for both D and Z.
! !
! type amg_dprec_type ! type amg_dprec_type
! type(amg_donelev_type), allocatable :: precv(:) ! type(amg_donelev_type), allocatable :: precv(:)
! end type amg_dprec_type ! end type amg_dprec_type
! !
! Note that the levels are numbered in increasing order starting from ! Note that the levels are numbered in increasing order starting from
! the level 1 as the finest one, and the number of levels is given by ! the level 1 as the finest one, and the number of levels is given by
! size(precv(:)) which is the id of the coarsest level. ! size(precv(:)) which is the id of the coarsest level.
! In the multigrid literature many authors number the levels in the opposite ! In the multigrid literature many authors number the levels in the opposite
! order, with level 0 being the id of the coarsest level. ! order, with level 0 being the id of the coarsest level.
! !
! !
integer, parameter, private :: wv_size_=4 integer, parameter, private :: wv_size_=4
type, extends(psb_dprec_type) :: amg_dprec_type type, extends(psb_dprec_type) :: amg_dprec_type
type(amg_daggr_data) :: ag_data type(amg_daggr_data) :: ag_data
! !
! Number of outer sweeps. Sometimes 2 V-cycles may be better than 1 W-cycle. ! Number of outer sweeps. Sometimes 2 V-cycles may be better than 1 W-cycle.
! !
integer(psb_ipk_) :: outer_sweeps = 1 integer(psb_ipk_) :: outer_sweeps = 1
! !
@ -97,11 +97,11 @@ module amg_d_prec_type
! to keep track against what is put later in the multilevel array ! to keep track against what is put later in the multilevel array
! !
integer(psb_ipk_) :: coarse_solver = -1 integer(psb_ipk_) :: coarse_solver = -1
! !
! The multilevel hierarchy ! The multilevel hierarchy
! !
type(amg_d_onelev_type), allocatable :: precv(:) type(amg_d_onelev_type), allocatable :: precv(:)
contains contains
procedure, pass(prec) :: psb_d_apply2_vect => amg_d_apply2_vect procedure, pass(prec) :: psb_d_apply2_vect => amg_d_apply2_vect
procedure, pass(prec) :: psb_d_apply1_vect => amg_d_apply1_vect procedure, pass(prec) :: psb_d_apply1_vect => amg_d_apply1_vect
@ -127,7 +127,7 @@ module amg_d_prec_type
procedure, pass(prec) :: cseti => amg_dcprecseti procedure, pass(prec) :: cseti => amg_dcprecseti
procedure, pass(prec) :: csetc => amg_dcprecsetc procedure, pass(prec) :: csetc => amg_dcprecsetc
procedure, pass(prec) :: csetr => amg_dcprecsetr procedure, pass(prec) :: csetr => amg_dcprecsetr
generic, public :: set => cseti, csetc, csetr, setsm, setsv, setag generic, public :: set => setsm, setsv, setag
procedure, pass(prec) :: get_smoother => amg_d_get_smootherp procedure, pass(prec) :: get_smoother => amg_d_get_smootherp
procedure, pass(prec) :: get_solver => amg_d_get_solverp procedure, pass(prec) :: get_solver => amg_d_get_solverp
procedure, pass(prec) :: move_alloc => d_prec_move_alloc procedure, pass(prec) :: move_alloc => d_prec_move_alloc
@ -157,7 +157,7 @@ module amg_d_prec_type
interface amg_precdescr interface amg_precdescr
subroutine amg_dfile_prec_descr(prec,iout,root) subroutine amg_dfile_prec_descr(prec,iout,root)
import :: amg_dprec_type, psb_ipk_ import :: amg_dprec_type, psb_ipk_
implicit none implicit none
! Arguments ! Arguments
class(amg_dprec_type), intent(in) :: prec class(amg_dprec_type), intent(in) :: prec
integer(psb_ipk_), intent(in), optional :: iout integer(psb_ipk_), intent(in), optional :: iout
@ -211,7 +211,7 @@ module amg_d_prec_type
end subroutine amg_dprecaply1 end subroutine amg_dprecaply1
end interface end interface
interface interface
subroutine amg_dprecsetsm(prec,val,info,ilev,ilmax,pos) subroutine amg_dprecsetsm(prec,val,info,ilev,ilmax,pos)
import :: psb_dspmat_type, psb_desc_type, psb_dpk_, & import :: psb_dspmat_type, psb_desc_type, psb_dpk_, &
& amg_dprec_type, amg_d_base_smoother_type, psb_ipk_ & amg_dprec_type, amg_d_base_smoother_type, psb_ipk_
@ -243,7 +243,7 @@ module amg_d_prec_type
import :: psb_dspmat_type, psb_desc_type, psb_dpk_, & import :: psb_dspmat_type, psb_desc_type, psb_dpk_, &
& amg_dprec_type, psb_ipk_ & amg_dprec_type, psb_ipk_
class(amg_dprec_type), intent(inout) :: prec class(amg_dprec_type), intent(inout) :: prec
character(len=*), intent(in) :: what character(len=*), intent(in) :: what
integer(psb_ipk_), intent(in) :: val integer(psb_ipk_), intent(in) :: val
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx
@ -253,7 +253,7 @@ module amg_d_prec_type
import :: psb_dspmat_type, psb_desc_type, psb_dpk_, & import :: psb_dspmat_type, psb_desc_type, psb_dpk_, &
& amg_dprec_type, psb_ipk_ & amg_dprec_type, psb_ipk_
class(amg_dprec_type), intent(inout) :: prec class(amg_dprec_type), intent(inout) :: prec
character(len=*), intent(in) :: what character(len=*), intent(in) :: what
real(psb_dpk_), intent(in) :: val real(psb_dpk_), intent(in) :: val
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx
@ -263,7 +263,7 @@ module amg_d_prec_type
import :: psb_dspmat_type, psb_desc_type, psb_dpk_, & import :: psb_dspmat_type, psb_desc_type, psb_dpk_, &
& amg_dprec_type, psb_ipk_ & amg_dprec_type, psb_ipk_
class(amg_dprec_type), intent(inout) :: prec class(amg_dprec_type), intent(inout) :: prec
character(len=*), intent(in) :: what character(len=*), intent(in) :: what
character(len=*), intent(in) :: string character(len=*), intent(in) :: string
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx
@ -341,28 +341,28 @@ module amg_d_prec_type
! character, intent(in),optional :: upd ! character, intent(in),optional :: upd
end subroutine amg_d_smoothers_bld end subroutine amg_d_smoothers_bld
end interface amg_smoothers_bld end interface amg_smoothers_bld
contains contains
! !
! Function returning a pointer to the smoother ! Function returning a pointer to the smoother
! !
function amg_d_get_smootherp(prec,ilev) result(val) function amg_d_get_smootherp(prec,ilev) result(val)
implicit none implicit none
class(amg_dprec_type), target, intent(in) :: prec class(amg_dprec_type), target, intent(in) :: prec
integer(psb_ipk_), optional :: ilev integer(psb_ipk_), optional :: ilev
class(amg_d_base_smoother_type), pointer :: val class(amg_d_base_smoother_type), pointer :: val
integer(psb_ipk_) :: ilev_ integer(psb_ipk_) :: ilev_
val => null() val => null()
if (present(ilev)) then if (present(ilev)) then
ilev_ = ilev ilev_ = ilev
else else
! What is a good default? ! What is a good default?
ilev_ = 1 ilev_ = 1
end if end if
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
if ((1<=ilev_).and.(ilev_<=size(prec%precv))) then if ((1<=ilev_).and.(ilev_<=size(prec%precv))) then
if (allocated(prec%precv(ilev_)%sm)) then if (allocated(prec%precv(ilev_)%sm)) then
val => prec%precv(ilev_)%sm val => prec%precv(ilev_)%sm
end if end if
end if end if
@ -372,23 +372,23 @@ contains
! Function returning a pointer to the solver ! Function returning a pointer to the solver
! !
function amg_d_get_solverp(prec,ilev) result(val) function amg_d_get_solverp(prec,ilev) result(val)
implicit none implicit none
class(amg_dprec_type), target, intent(in) :: prec class(amg_dprec_type), target, intent(in) :: prec
integer(psb_ipk_), optional :: ilev integer(psb_ipk_), optional :: ilev
class(amg_d_base_solver_type), pointer :: val class(amg_d_base_solver_type), pointer :: val
integer(psb_ipk_) :: ilev_ integer(psb_ipk_) :: ilev_
val => null() val => null()
if (present(ilev)) then if (present(ilev)) then
ilev_ = ilev ilev_ = ilev
else else
! What is a good default? ! What is a good default?
ilev_ = 1 ilev_ = 1
end if end if
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
if ((1<=ilev_).and.(ilev_<=size(prec%precv))) then if ((1<=ilev_).and.(ilev_<=size(prec%precv))) then
if (allocated(prec%precv(ilev_)%sm)) then if (allocated(prec%precv(ilev_)%sm)) then
if (allocated(prec%precv(ilev_)%sm%sv)) then if (allocated(prec%precv(ilev_)%sm%sv)) then
val => prec%precv(ilev_)%sm%sv val => prec%precv(ilev_)%sm%sv
end if end if
end if end if
@ -399,25 +399,25 @@ contains
! Function returning the size of the precv(:) array ! Function returning the size of the precv(:) array
! !
function amg_d_get_nlevs(prec) result(val) function amg_d_get_nlevs(prec) result(val)
implicit none implicit none
class(amg_dprec_type), intent(in) :: prec class(amg_dprec_type), intent(in) :: prec
integer(psb_ipk_) :: val integer(psb_ipk_) :: val
val = 0 val = 0
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
val = size(prec%precv) val = size(prec%precv)
end if end if
end function amg_d_get_nlevs end function amg_d_get_nlevs
! !
! Function returning the size of the amg_prec_type data structure ! Function returning the size of the amg_prec_type data structure
! in bytes or in number of nonzeros of the operator(s) involved. ! in bytes or in number of nonzeros of the operator(s) involved.
! !
function amg_d_get_nzeros(prec) result(val) function amg_d_get_nzeros(prec) result(val)
implicit none implicit none
class(amg_dprec_type), intent(in) :: prec class(amg_dprec_type), intent(in) :: prec
integer(psb_epk_) :: val integer(psb_epk_) :: val
integer(psb_ipk_) :: i integer(psb_ipk_) :: i
val = 0 val = 0
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
do i=1, size(prec%precv) do i=1, size(prec%precv)
val = val + prec%precv(i)%get_nzeros() val = val + prec%precv(i)%get_nzeros()
end do end do
@ -425,13 +425,13 @@ contains
end function amg_d_get_nzeros end function amg_d_get_nzeros
function amg_dprec_sizeof(prec) result(val) function amg_dprec_sizeof(prec) result(val)
implicit none implicit none
class(amg_dprec_type), intent(in) :: prec class(amg_dprec_type), intent(in) :: prec
integer(psb_epk_) :: val integer(psb_epk_) :: val
integer(psb_ipk_) :: i integer(psb_ipk_) :: i
val = 0 val = 0
val = val + psb_sizeof_ip val = val + psb_sizeof_ip
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
do i=1, size(prec%precv) do i=1, size(prec%precv)
val = val + prec%precv(i)%sizeof() val = val + prec%precv(i)%sizeof()
end do end do
@ -444,40 +444,40 @@ contains
! various level to the nonzeroes at the fine level ! various level to the nonzeroes at the fine level
! (original matrix) ! (original matrix)
! !
function amg_d_get_compl(prec) result(val) function amg_d_get_compl(prec) result(val)
implicit none implicit none
class(amg_dprec_type), intent(in) :: prec class(amg_dprec_type), intent(in) :: prec
real(psb_dpk_) :: val real(psb_dpk_) :: val
val = prec%ag_data%op_complexity val = prec%ag_data%op_complexity
end function amg_d_get_compl end function amg_d_get_compl
subroutine amg_d_cmp_compl(prec)
implicit none subroutine amg_d_cmp_compl(prec)
implicit none
class(amg_dprec_type), intent(inout) :: prec class(amg_dprec_type), intent(inout) :: prec
real(psb_dpk_) :: num, den, nmin real(psb_dpk_) :: num, den, nmin
type(psb_ctxt_type) :: ctxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: il integer(psb_ipk_) :: il
num = -done num = -done
den = done den = done
ctxt = prec%ctxt ctxt = prec%ctxt
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
il = 1 il = 1
num = prec%precv(il)%base_a%get_nzeros() num = prec%precv(il)%base_a%get_nzeros()
if (num >= dzero) then if (num >= dzero) then
den = num den = num
do il=2,size(prec%precv) do il=2,size(prec%precv)
num = num + max(0,prec%precv(il)%base_a%get_nzeros()) num = num + max(0,prec%precv(il)%base_a%get_nzeros())
end do end do
end if end if
end if end if
nmin = num nmin = num
call psb_min(ctxt,nmin) call psb_min(ctxt,nmin)
if (nmin < dzero) then if (nmin < dzero) then
num = dzero num = dzero
den = done den = done
@ -487,25 +487,25 @@ contains
end if end if
prec%ag_data%op_complexity = num/den prec%ag_data%op_complexity = num/den
end subroutine amg_d_cmp_compl end subroutine amg_d_cmp_compl
! !
! Average coarsening ratio ! Average coarsening ratio
! !
function amg_d_get_avg_cr(prec) result(val) function amg_d_get_avg_cr(prec) result(val)
implicit none implicit none
class(amg_dprec_type), intent(in) :: prec class(amg_dprec_type), intent(in) :: prec
real(psb_dpk_) :: val real(psb_dpk_) :: val
val = prec%ag_data%avg_cr val = prec%ag_data%avg_cr
end function amg_d_get_avg_cr end function amg_d_get_avg_cr
subroutine amg_d_cmp_avg_cr(prec)
implicit none subroutine amg_d_cmp_avg_cr(prec)
implicit none
class(amg_dprec_type), intent(inout) :: prec class(amg_dprec_type), intent(inout) :: prec
real(psb_dpk_) :: avgcr real(psb_dpk_) :: avgcr
type(psb_ctxt_type) :: ctxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: il, nl, iam, np integer(psb_ipk_) :: il, nl, iam, np
@ -519,12 +519,12 @@ contains
do il=2,nl do il=2,nl
avgcr = avgcr + max(dzero,prec%precv(il)%szratio) avgcr = avgcr + max(dzero,prec%precv(il)%szratio)
end do end do
avgcr = avgcr / (nl-1) avgcr = avgcr / (nl-1)
end if end if
call psb_sum(ctxt,avgcr) call psb_sum(ctxt,avgcr)
prec%ag_data%avg_cr = avgcr/np prec%ag_data%avg_cr = avgcr/np
end subroutine amg_d_cmp_avg_cr end subroutine amg_d_cmp_avg_cr
! !
! Subroutines: amg_Tprec_free ! Subroutines: amg_Tprec_free
! Version: real ! Version: real
@ -538,74 +538,74 @@ contains
! error code. ! error code.
! !
subroutine amg_dprecfree(p,info) subroutine amg_dprecfree(p,info)
implicit none implicit none
! Arguments ! Arguments
type(amg_dprec_type), intent(inout) :: p type(amg_dprec_type), intent(inout) :: p
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
! Local variables ! Local variables
integer(psb_ipk_) :: me,err_act,i integer(psb_ipk_) :: me,err_act,i
character(len=20) :: name character(len=20) :: name
info=psb_success_ info=psb_success_
name = 'amg_dprecfree' name = 'amg_dprecfree'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; return info = psb_err_internal_error_; return
end if end if
me=-1 me=-1
call p%free(info) call p%free(info)
return return
end subroutine amg_dprecfree end subroutine amg_dprecfree
subroutine amg_d_prec_free(prec,info) subroutine amg_d_prec_free(prec,info)
implicit none implicit none
! Arguments ! Arguments
class(amg_dprec_type), intent(inout) :: prec class(amg_dprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
! Local variables ! Local variables
integer(psb_ipk_) :: me,err_act,i integer(psb_ipk_) :: me,err_act,i
character(len=20) :: name character(len=20) :: name
info=psb_success_ info=psb_success_
name = 'amg_dprecfree' name = 'amg_dprecfree'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999 info = psb_err_internal_error_; goto 9999
end if end if
me=-1 me=-1
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
do i=1,size(prec%precv) do i=1,size(prec%precv)
call prec%precv(i)%free(info) call prec%precv(i)%free(info)
end do end do
deallocate(prec%precv,stat=info) deallocate(prec%precv,stat=info)
end if end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
9999 call psb_error_handler(err_act) 9999 call psb_error_handler(err_act)
return return
end subroutine amg_d_prec_free end subroutine amg_d_prec_free
! !
! Top level methods. ! Top level methods.
! !
subroutine amg_d_apply2_vect(prec,x,y,desc_data,info,trans,work) subroutine amg_d_apply2_vect(prec,x,y,desc_data,info,trans,work)
implicit none implicit none
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
class(amg_dprec_type), intent(inout) :: prec class(amg_dprec_type), intent(inout) :: prec
type(psb_d_vect_type),intent(inout) :: x type(psb_d_vect_type),intent(inout) :: x
@ -618,13 +618,13 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
select type(prec) select type(prec)
type is (amg_dprec_type) type is (amg_dprec_type)
call amg_precapply(prec,x,y,desc_data,info,trans,work) call amg_precapply(prec,x,y,desc_data,info,trans,work)
class default class default
info = psb_err_missing_override_method_ info = psb_err_missing_override_method_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end select end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -636,7 +636,7 @@ contains
end subroutine amg_d_apply2_vect end subroutine amg_d_apply2_vect
subroutine amg_d_apply1_vect(prec,x,desc_data,info,trans,work) subroutine amg_d_apply1_vect(prec,x,desc_data,info,trans,work)
implicit none implicit none
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
class(amg_dprec_type), intent(inout) :: prec class(amg_dprec_type), intent(inout) :: prec
type(psb_d_vect_type),intent(inout) :: x type(psb_d_vect_type),intent(inout) :: x
@ -648,13 +648,13 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
select type(prec) select type(prec)
type is (amg_dprec_type) type is (amg_dprec_type)
call amg_precapply(prec,x,desc_data,info,trans,work) call amg_precapply(prec,x,desc_data,info,trans,work)
class default class default
info = psb_err_missing_override_method_ info = psb_err_missing_override_method_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end select end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -667,7 +667,7 @@ contains
subroutine amg_d_apply2v(prec,x,y,desc_data,info,trans,work) subroutine amg_d_apply2v(prec,x,y,desc_data,info,trans,work)
implicit none implicit none
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
class(amg_dprec_type), intent(inout) :: prec class(amg_dprec_type), intent(inout) :: prec
real(psb_dpk_),intent(inout) :: x(:) real(psb_dpk_),intent(inout) :: x(:)
@ -680,13 +680,13 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
select type(prec) select type(prec)
type is (amg_dprec_type) type is (amg_dprec_type)
call amg_precapply(prec,x,y,desc_data,info,trans,work) call amg_precapply(prec,x,y,desc_data,info,trans,work)
class default class default
info = psb_err_missing_override_method_ info = psb_err_missing_override_method_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end select end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -698,7 +698,7 @@ contains
end subroutine amg_d_apply2v end subroutine amg_d_apply2v
subroutine amg_d_apply1v(prec,x,desc_data,info,trans) subroutine amg_d_apply1v(prec,x,desc_data,info,trans)
implicit none implicit none
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
class(amg_dprec_type), intent(inout) :: prec class(amg_dprec_type), intent(inout) :: prec
real(psb_dpk_),intent(inout) :: x(:) real(psb_dpk_),intent(inout) :: x(:)
@ -709,13 +709,13 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
select type(prec) select type(prec)
type is (amg_dprec_type) type is (amg_dprec_type)
call amg_precapply(prec,x,desc_data,info,trans) call amg_precapply(prec,x,desc_data,info,trans)
class default class default
info = psb_err_missing_override_method_ info = psb_err_missing_override_method_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end select end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -730,8 +730,8 @@ contains
subroutine amg_d_dump(prec,info,istart,iend,iproc,prefix,head,& subroutine amg_d_dump(prec,info,istart,iend,iproc,prefix,head,&
& ac,rp,smoother,solver,tprol,& & ac,rp,smoother,solver,tprol,&
& global_num) & global_num)
implicit none implicit none
class(amg_dprec_type), intent(in) :: prec class(amg_dprec_type), intent(in) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: istart, iend, iproc integer(psb_ipk_), intent(in), optional :: istart, iend, iproc
@ -742,27 +742,27 @@ contains
integer(psb_ipk_) :: iam, np, iproc_ integer(psb_ipk_) :: iam, np, iproc_
character(len=80) :: prefix_ character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than character(len=120) :: fname ! len should be at least 20 more than
! len of prefix_ ! len of prefix_
info = 0 info = 0
icontxt = prec%ctxt icontxt = prec%ctxt
call psb_info(icontxt,iam,np) call psb_info(icontxt,iam,np)
iln = size(prec%precv) iln = size(prec%precv)
if (present(istart)) then if (present(istart)) then
il1 = max(1,istart) il1 = max(1,istart)
else else
il1 = min(2,iln) il1 = min(2,iln)
end if end if
if (present(iend)) then if (present(iend)) then
iln = min(iln, iend) iln = min(iln, iend)
end if end if
iproc_ = -1 iproc_ = -1
if (present(iproc)) then if (present(iproc)) then
iproc_ = iproc iproc_ = iproc
end if end if
if ((iproc_ == -1).or.(iproc_==iam)) then if ((iproc_ == -1).or.(iproc_==iam)) then
do lev=il1, iln do lev=il1, iln
call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,& call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,&
& ac=ac,smoother=smoother,solver=solver,rp=rp,tprol=tprol, & & ac=ac,smoother=smoother,solver=solver,rp=rp,tprol=tprol, &
@ -773,7 +773,7 @@ contains
subroutine amg_d_cnv(prec,info,amold,vmold,imold) subroutine amg_d_cnv(prec,info,amold,vmold,imold)
implicit none implicit none
class(amg_dprec_type), intent(inout) :: prec class(amg_dprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
class(psb_d_base_sparse_mat), intent(in), optional :: amold class(psb_d_base_sparse_mat), intent(in), optional :: amold
@ -781,7 +781,7 @@ contains
class(psb_i_base_vect_type), intent(in), optional :: imold class(psb_i_base_vect_type), intent(in), optional :: imold
integer(psb_ipk_) :: i integer(psb_ipk_) :: i
info = psb_success_ info = psb_success_
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
do i=1,size(prec%precv) do i=1,size(prec%precv)
@ -789,24 +789,24 @@ contains
& call prec%precv(i)%cnv(info,amold=amold,vmold=vmold,imold=imold) & call prec%precv(i)%cnv(info,amold=amold,vmold=vmold,imold=imold)
end do end do
end if end if
end subroutine amg_d_cnv end subroutine amg_d_cnv
subroutine amg_d_clone(prec,precout,info) subroutine amg_d_clone(prec,precout,info)
implicit none implicit none
class(amg_dprec_type), intent(inout) :: prec class(amg_dprec_type), intent(inout) :: prec
class(psb_dprec_type), intent(inout) :: precout class(psb_dprec_type), intent(inout) :: precout
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
call precout%free(info) call precout%free(info)
if (info == 0) call amg_d_inner_clone(prec,precout,info) if (info == 0) call amg_d_inner_clone(prec,precout,info)
end subroutine amg_d_clone end subroutine amg_d_clone
subroutine amg_d_inner_clone(prec,precout,info) subroutine amg_d_inner_clone(prec,precout,info)
implicit none implicit none
class(amg_dprec_type), intent(inout) :: prec class(amg_dprec_type), intent(inout) :: prec
class(psb_dprec_type), target, intent(inout) :: precout class(psb_dprec_type), target, intent(inout) :: precout
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -821,17 +821,17 @@ contains
pout%ctxt = prec%ctxt pout%ctxt = prec%ctxt
pout%ag_data = prec%ag_data pout%ag_data = prec%ag_data
pout%outer_sweeps = prec%outer_sweeps pout%outer_sweeps = prec%outer_sweeps
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
ln = size(prec%precv) ln = size(prec%precv)
allocate(pout%precv(ln),stat=info) allocate(pout%precv(ln),stat=info)
if (info /= psb_success_) goto 9999 if (info /= psb_success_) goto 9999
if (ln >= 1) then if (ln >= 1) then
call prec%precv(1)%clone(pout%precv(1),info) call prec%precv(1)%clone(pout%precv(1),info)
end if end if
do lev=2, ln do lev=2, ln
if (info /= psb_success_) exit if (info /= psb_success_) exit
call prec%precv(lev)%clone(pout%precv(lev),info) call prec%precv(lev)%clone(pout%precv(lev),info)
if (info == psb_success_) then if (info == psb_success_) then
pout%precv(lev)%base_a => pout%precv(lev)%ac pout%precv(lev)%base_a => pout%precv(lev)%ac
pout%precv(lev)%base_desc => pout%precv(lev)%desc_ac pout%precv(lev)%base_desc => pout%precv(lev)%desc_ac
pout%precv(lev)%map%p_desc_U => pout%precv(lev-1)%base_desc pout%precv(lev)%map%p_desc_U => pout%precv(lev-1)%base_desc
@ -842,7 +842,7 @@ contains
if (allocated(prec%precv(1)%wrk)) & if (allocated(prec%precv(1)%wrk)) &
& call pout%allocate_wrk(info,vmold=prec%precv(1)%wrk%vx2l%v) & call pout%allocate_wrk(info,vmold=prec%precv(1)%wrk%vx2l%v)
class default class default
write(0,*) 'Error: wrong out type' write(0,*) 'Error: wrong out type'
info = psb_err_invalid_input_ info = psb_err_invalid_input_
end select end select
@ -854,14 +854,14 @@ contains
implicit none implicit none
class(amg_dprec_type), intent(inout) :: prec class(amg_dprec_type), intent(inout) :: prec
class(amg_dprec_type), intent(inout), target :: b class(amg_dprec_type), intent(inout), target :: b
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i integer(psb_ipk_) :: i
if (same_type_as(prec,b)) then if (same_type_as(prec,b)) then
if (allocated(b%precv)) then if (allocated(b%precv)) then
! This might not be required if FINAL procedures are available. ! This might not be required if FINAL procedures are available.
call b%free(info) call b%free(info)
if (info /= psb_success_) then if (info /= psb_success_) then
!????? !?????
!!$ return !!$ return
endif endif
@ -869,7 +869,7 @@ contains
b%ctxt = prec%ctxt b%ctxt = prec%ctxt
b%ag_data = prec%ag_data b%ag_data = prec%ag_data
b%outer_sweeps = prec%outer_sweeps b%outer_sweeps = prec%outer_sweeps
call move_alloc(prec%precv,b%precv) call move_alloc(prec%precv,b%precv)
! Fix the pointers except on level 1. ! Fix the pointers except on level 1.
do i=2, size(b%precv) do i=2, size(b%precv)
@ -878,7 +878,7 @@ contains
b%precv(i)%map%p_desc_U => b%precv(i-1)%base_desc b%precv(i)%map%p_desc_U => b%precv(i-1)%base_desc
b%precv(i)%map%p_desc_V => b%precv(i)%base_desc b%precv(i)%map%p_desc_V => b%precv(i)%base_desc
end do end do
else else
write(0,*) 'Warning: PREC%move_alloc onto different type?' write(0,*) 'Warning: PREC%move_alloc onto different type?'
info = psb_err_internal_error_ info = psb_err_internal_error_
@ -888,7 +888,7 @@ contains
subroutine amg_d_allocate_wrk(prec,info,vmold,desc) subroutine amg_d_allocate_wrk(prec,info,vmold,desc)
use psb_base_mod use psb_base_mod
implicit none implicit none
! Arguments ! Arguments
class(amg_dprec_type), intent(inout) :: prec class(amg_dprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -896,37 +896,37 @@ contains
! !
! In MLD the DESC optional argument is ignored, since ! In MLD the DESC optional argument is ignored, since
! the necessary info is contained in the various entries of the ! the necessary info is contained in the various entries of the
! PRECV component. ! PRECV component.
type(psb_desc_type), intent(in), optional :: desc type(psb_desc_type), intent(in), optional :: desc
! Local variables ! Local variables
integer(psb_ipk_) :: me,err_act,i,j,level,nlev, nc2l integer(psb_ipk_) :: me,err_act,i,j,level,nlev, nc2l
character(len=20) :: name character(len=20) :: name
info=psb_success_ info=psb_success_
name = 'amg_d_allocate_wrk' name = 'amg_d_allocate_wrk'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999 info = psb_err_internal_error_; goto 9999
end if end if
nlev = size(prec%precv) nlev = size(prec%precv)
level = 1 level = 1
do level = 1, nlev do level = 1, nlev
call prec%precv(level)%allocate_wrk(info,vmold=vmold) call prec%precv(level)%allocate_wrk(info,vmold=vmold)
if (psb_errstatus_fatal()) then if (psb_errstatus_fatal()) then
nc2l = prec%precv(level)%base_desc%get_local_cols() nc2l = prec%precv(level)%base_desc%get_local_cols()
info=psb_err_alloc_request_ info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l/), a_err='real(psb_dpk_)') call psb_errpush(info,name,i_err=(/2*nc2l/), a_err='real(psb_dpk_)')
goto 9999 goto 9999
end if end if
end do end do
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
9999 call psb_error_handler(err_act) 9999 call psb_error_handler(err_act)
return return
end subroutine amg_d_allocate_wrk end subroutine amg_d_allocate_wrk
subroutine amg_d_free_wrk(prec,info) subroutine amg_d_free_wrk(prec,info)
@ -948,13 +948,13 @@ contains
info = psb_err_internal_error_; goto 9999 info = psb_err_internal_error_; goto 9999
end if end if
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
nlev = size(prec%precv) nlev = size(prec%precv)
do level = 1, nlev do level = 1, nlev
call prec%precv(level)%free_wrk(info) call prec%precv(level)%free_wrk(info)
end do end do
end if end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -966,7 +966,7 @@ contains
function amg_d_is_allocated_wrk(prec) result(res) function amg_d_is_allocated_wrk(prec) result(res)
use psb_base_mod use psb_base_mod
implicit none implicit none
! Arguments ! Arguments
class(amg_dprec_type), intent(in) :: prec class(amg_dprec_type), intent(in) :: prec
logical :: res logical :: res
@ -974,7 +974,7 @@ contains
res = .false. res = .false.
if (.not.allocated(prec%precv)) return if (.not.allocated(prec%precv)) return
res = allocated(prec%precv(1)%wrk) res = allocated(prec%precv(1)%wrk)
end function amg_d_is_allocated_wrk end function amg_d_is_allocated_wrk
end module amg_d_prec_type end module amg_d_prec_type

292
amgprec/amg_s_prec_type.f90
Unescape Escape View File

@ -1,15 +1,15 @@
! !
! !
! AMG4PSBLAS version 1.0 ! AMG4PSBLAS version 1.0
! Algebraic Multigrid Package ! Algebraic Multigrid Package
! based on PSBLAS (Parallel Sparse BLAS version 3.5) ! based on PSBLAS (Parallel Sparse BLAS version 3.5)
! !
! (C) Copyright 2020 ! (C) Copyright 2020
! !
! Salvatore Filippone ! Salvatore Filippone
! Pasqua D'Ambra ! Pasqua D'Ambra
! Fabio Durastante ! Fabio Durastante
! !
! Redistribution and use in source and binary forms, with or without ! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions ! modification, are permitted provided that the following conditions
! are met: ! are met:
@ -21,7 +21,7 @@
! 3. The name of the AMG4PSBLAS group or the names of its contributors may ! 3. The name of the AMG4PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this ! not be used to endorse or promote products derived from this
! software without specific written permission. ! software without specific written permission.
! !
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
@ -33,21 +33,21 @@
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! 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 ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE. ! POSSIBILITY OF SUCH DAMAGE.
! !
! !
! File: amg_s_prec_type.f90 ! File: amg_s_prec_type.f90
! !
! Module: amg_s_prec_type ! Module: amg_s_prec_type
! !
! This module defines: ! This module defines:
! - the amg_s_prec_type data structure containing the preconditioner and related ! - the amg_s_prec_type data structure containing the preconditioner and related
! data structures; ! data structures;
! !
! It contains routines for ! It contains routines for
! - Building and applying; ! - Building and applying;
! - checking if the preconditioner is correctly defined; ! - checking if the preconditioner is correctly defined;
! - printing a description of the preconditioner; ! - printing a description of the preconditioner;
! - deallocating the preconditioner data structure. ! - deallocating the preconditioner data structure.
! !
module amg_s_prec_type module amg_s_prec_type
@ -70,25 +70,25 @@ module amg_s_prec_type
! It consists of an array of 'one-level' intermediate data structures ! It consists of an array of 'one-level' intermediate data structures
! of type amg_sonelev_type, each containing the information needed to apply ! of type amg_sonelev_type, each containing the information needed to apply
! the smoothing and the coarse-space correction at a generic level. RT is the ! the smoothing and the coarse-space correction at a generic level. RT is the
! real data type, i.e. S for both S and C, and D for both D and Z. ! real data type, i.e. S for both S and C, and D for both D and Z.
! !
! type amg_sprec_type ! type amg_sprec_type
! type(amg_sonelev_type), allocatable :: precv(:) ! type(amg_sonelev_type), allocatable :: precv(:)
! end type amg_sprec_type ! end type amg_sprec_type
! !
! Note that the levels are numbered in increasing order starting from ! Note that the levels are numbered in increasing order starting from
! the level 1 as the finest one, and the number of levels is given by ! the level 1 as the finest one, and the number of levels is given by
! size(precv(:)) which is the id of the coarsest level. ! size(precv(:)) which is the id of the coarsest level.
! In the multigrid literature many authors number the levels in the opposite ! In the multigrid literature many authors number the levels in the opposite
! order, with level 0 being the id of the coarsest level. ! order, with level 0 being the id of the coarsest level.
! !
! !
integer, parameter, private :: wv_size_=4 integer, parameter, private :: wv_size_=4
type, extends(psb_sprec_type) :: amg_sprec_type type, extends(psb_sprec_type) :: amg_sprec_type
type(amg_saggr_data) :: ag_data type(amg_saggr_data) :: ag_data
! !
! Number of outer sweeps. Sometimes 2 V-cycles may be better than 1 W-cycle. ! Number of outer sweeps. Sometimes 2 V-cycles may be better than 1 W-cycle.
! !
integer(psb_ipk_) :: outer_sweeps = 1 integer(psb_ipk_) :: outer_sweeps = 1
! !
@ -97,11 +97,11 @@ module amg_s_prec_type
! to keep track against what is put later in the multilevel array ! to keep track against what is put later in the multilevel array
! !
integer(psb_ipk_) :: coarse_solver = -1 integer(psb_ipk_) :: coarse_solver = -1
! !
! The multilevel hierarchy ! The multilevel hierarchy
! !
type(amg_s_onelev_type), allocatable :: precv(:) type(amg_s_onelev_type), allocatable :: precv(:)
contains contains
procedure, pass(prec) :: psb_s_apply2_vect => amg_s_apply2_vect procedure, pass(prec) :: psb_s_apply2_vect => amg_s_apply2_vect
procedure, pass(prec) :: psb_s_apply1_vect => amg_s_apply1_vect procedure, pass(prec) :: psb_s_apply1_vect => amg_s_apply1_vect
@ -127,7 +127,7 @@ module amg_s_prec_type
procedure, pass(prec) :: cseti => amg_scprecseti procedure, pass(prec) :: cseti => amg_scprecseti
procedure, pass(prec) :: csetc => amg_scprecsetc procedure, pass(prec) :: csetc => amg_scprecsetc
procedure, pass(prec) :: csetr => amg_scprecsetr procedure, pass(prec) :: csetr => amg_scprecsetr
generic, public :: set => cseti, csetc, csetr, setsm, setsv, setag generic, public :: set => setsm, setsv, setag
procedure, pass(prec) :: get_smoother => amg_s_get_smootherp procedure, pass(prec) :: get_smoother => amg_s_get_smootherp
procedure, pass(prec) :: get_solver => amg_s_get_solverp procedure, pass(prec) :: get_solver => amg_s_get_solverp
procedure, pass(prec) :: move_alloc => s_prec_move_alloc procedure, pass(prec) :: move_alloc => s_prec_move_alloc
@ -157,7 +157,7 @@ module amg_s_prec_type
interface amg_precdescr interface amg_precdescr
subroutine amg_sfile_prec_descr(prec,iout,root) subroutine amg_sfile_prec_descr(prec,iout,root)
import :: amg_sprec_type, psb_ipk_ import :: amg_sprec_type, psb_ipk_
implicit none implicit none
! Arguments ! Arguments
class(amg_sprec_type), intent(in) :: prec class(amg_sprec_type), intent(in) :: prec
integer(psb_ipk_), intent(in), optional :: iout integer(psb_ipk_), intent(in), optional :: iout
@ -211,7 +211,7 @@ module amg_s_prec_type
end subroutine amg_sprecaply1 end subroutine amg_sprecaply1
end interface end interface
interface interface
subroutine amg_sprecsetsm(prec,val,info,ilev,ilmax,pos) subroutine amg_sprecsetsm(prec,val,info,ilev,ilmax,pos)
import :: psb_sspmat_type, psb_desc_type, psb_spk_, & import :: psb_sspmat_type, psb_desc_type, psb_spk_, &
& amg_sprec_type, amg_s_base_smoother_type, psb_ipk_ & amg_sprec_type, amg_s_base_smoother_type, psb_ipk_
@ -243,7 +243,7 @@ module amg_s_prec_type
import :: psb_sspmat_type, psb_desc_type, psb_spk_, & import :: psb_sspmat_type, psb_desc_type, psb_spk_, &
& amg_sprec_type, psb_ipk_ & amg_sprec_type, psb_ipk_
class(amg_sprec_type), intent(inout) :: prec class(amg_sprec_type), intent(inout) :: prec
character(len=*), intent(in) :: what character(len=*), intent(in) :: what
integer(psb_ipk_), intent(in) :: val integer(psb_ipk_), intent(in) :: val
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx
@ -253,7 +253,7 @@ module amg_s_prec_type
import :: psb_sspmat_type, psb_desc_type, psb_spk_, & import :: psb_sspmat_type, psb_desc_type, psb_spk_, &
& amg_sprec_type, psb_ipk_ & amg_sprec_type, psb_ipk_
class(amg_sprec_type), intent(inout) :: prec class(amg_sprec_type), intent(inout) :: prec
character(len=*), intent(in) :: what character(len=*), intent(in) :: what
real(psb_spk_), intent(in) :: val real(psb_spk_), intent(in) :: val
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx
@ -263,7 +263,7 @@ module amg_s_prec_type
import :: psb_sspmat_type, psb_desc_type, psb_spk_, & import :: psb_sspmat_type, psb_desc_type, psb_spk_, &
& amg_sprec_type, psb_ipk_ & amg_sprec_type, psb_ipk_
class(amg_sprec_type), intent(inout) :: prec class(amg_sprec_type), intent(inout) :: prec
character(len=*), intent(in) :: what character(len=*), intent(in) :: what
character(len=*), intent(in) :: string character(len=*), intent(in) :: string
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx
@ -341,28 +341,28 @@ module amg_s_prec_type
! character, intent(in),optional :: upd ! character, intent(in),optional :: upd
end subroutine amg_s_smoothers_bld end subroutine amg_s_smoothers_bld
end interface amg_smoothers_bld end interface amg_smoothers_bld
contains contains
! !
! Function returning a pointer to the smoother ! Function returning a pointer to the smoother
! !
function amg_s_get_smootherp(prec,ilev) result(val) function amg_s_get_smootherp(prec,ilev) result(val)
implicit none implicit none
class(amg_sprec_type), target, intent(in) :: prec class(amg_sprec_type), target, intent(in) :: prec
integer(psb_ipk_), optional :: ilev integer(psb_ipk_), optional :: ilev
class(amg_s_base_smoother_type), pointer :: val class(amg_s_base_smoother_type), pointer :: val
integer(psb_ipk_) :: ilev_ integer(psb_ipk_) :: ilev_
val => null() val => null()
if (present(ilev)) then if (present(ilev)) then
ilev_ = ilev ilev_ = ilev
else else
! What is a good default? ! What is a good default?
ilev_ = 1 ilev_ = 1
end if end if
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
if ((1<=ilev_).and.(ilev_<=size(prec%precv))) then if ((1<=ilev_).and.(ilev_<=size(prec%precv))) then
if (allocated(prec%precv(ilev_)%sm)) then if (allocated(prec%precv(ilev_)%sm)) then
val => prec%precv(ilev_)%sm val => prec%precv(ilev_)%sm
end if end if
end if end if
@ -372,23 +372,23 @@ contains
! Function returning a pointer to the solver ! Function returning a pointer to the solver
! !
function amg_s_get_solverp(prec,ilev) result(val) function amg_s_get_solverp(prec,ilev) result(val)
implicit none implicit none
class(amg_sprec_type), target, intent(in) :: prec class(amg_sprec_type), target, intent(in) :: prec
integer(psb_ipk_), optional :: ilev integer(psb_ipk_), optional :: ilev
class(amg_s_base_solver_type), pointer :: val class(amg_s_base_solver_type), pointer :: val
integer(psb_ipk_) :: ilev_ integer(psb_ipk_) :: ilev_
val => null() val => null()
if (present(ilev)) then if (present(ilev)) then
ilev_ = ilev ilev_ = ilev
else else
! What is a good default? ! What is a good default?
ilev_ = 1 ilev_ = 1
end if end if
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
if ((1<=ilev_).and.(ilev_<=size(prec%precv))) then if ((1<=ilev_).and.(ilev_<=size(prec%precv))) then
if (allocated(prec%precv(ilev_)%sm)) then if (allocated(prec%precv(ilev_)%sm)) then
if (allocated(prec%precv(ilev_)%sm%sv)) then if (allocated(prec%precv(ilev_)%sm%sv)) then
val => prec%precv(ilev_)%sm%sv val => prec%precv(ilev_)%sm%sv
end if end if
end if end if
@ -399,25 +399,25 @@ contains
! Function returning the size of the precv(:) array ! Function returning the size of the precv(:) array
! !
function amg_s_get_nlevs(prec) result(val) function amg_s_get_nlevs(prec) result(val)
implicit none implicit none
class(amg_sprec_type), intent(in) :: prec class(amg_sprec_type), intent(in) :: prec
integer(psb_ipk_) :: val integer(psb_ipk_) :: val
val = 0 val = 0
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
val = size(prec%precv) val = size(prec%precv)
end if end if
end function amg_s_get_nlevs end function amg_s_get_nlevs
! !
! Function returning the size of the amg_prec_type data structure ! Function returning the size of the amg_prec_type data structure
! in bytes or in number of nonzeros of the operator(s) involved. ! in bytes or in number of nonzeros of the operator(s) involved.
! !
function amg_s_get_nzeros(prec) result(val) function amg_s_get_nzeros(prec) result(val)
implicit none implicit none
class(amg_sprec_type), intent(in) :: prec class(amg_sprec_type), intent(in) :: prec
integer(psb_epk_) :: val integer(psb_epk_) :: val
integer(psb_ipk_) :: i integer(psb_ipk_) :: i
val = 0 val = 0
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
do i=1, size(prec%precv) do i=1, size(prec%precv)
val = val + prec%precv(i)%get_nzeros() val = val + prec%precv(i)%get_nzeros()
end do end do
@ -425,13 +425,13 @@ contains
end function amg_s_get_nzeros end function amg_s_get_nzeros
function amg_sprec_sizeof(prec) result(val) function amg_sprec_sizeof(prec) result(val)
implicit none implicit none
class(amg_sprec_type), intent(in) :: prec class(amg_sprec_type), intent(in) :: prec
integer(psb_epk_) :: val integer(psb_epk_) :: val
integer(psb_ipk_) :: i integer(psb_ipk_) :: i
val = 0 val = 0
val = val + psb_sizeof_ip val = val + psb_sizeof_ip
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
do i=1, size(prec%precv) do i=1, size(prec%precv)
val = val + prec%precv(i)%sizeof() val = val + prec%precv(i)%sizeof()
end do end do
@ -444,40 +444,40 @@ contains
! various level to the nonzeroes at the fine level ! various level to the nonzeroes at the fine level
! (original matrix) ! (original matrix)
! !
function amg_s_get_compl(prec) result(val) function amg_s_get_compl(prec) result(val)
implicit none implicit none
class(amg_sprec_type), intent(in) :: prec class(amg_sprec_type), intent(in) :: prec
real(psb_spk_) :: val real(psb_spk_) :: val
val = prec%ag_data%op_complexity val = prec%ag_data%op_complexity
end function amg_s_get_compl end function amg_s_get_compl
subroutine amg_s_cmp_compl(prec)
implicit none subroutine amg_s_cmp_compl(prec)
implicit none
class(amg_sprec_type), intent(inout) :: prec class(amg_sprec_type), intent(inout) :: prec
real(psb_spk_) :: num, den, nmin real(psb_spk_) :: num, den, nmin
type(psb_ctxt_type) :: ctxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: il integer(psb_ipk_) :: il
num = -sone num = -sone
den = sone den = sone
ctxt = prec%ctxt ctxt = prec%ctxt
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
il = 1 il = 1
num = prec%precv(il)%base_a%get_nzeros() num = prec%precv(il)%base_a%get_nzeros()
if (num >= szero) then if (num >= szero) then
den = num den = num
do il=2,size(prec%precv) do il=2,size(prec%precv)
num = num + max(0,prec%precv(il)%base_a%get_nzeros()) num = num + max(0,prec%precv(il)%base_a%get_nzeros())
end do end do
end if end if
end if end if
nmin = num nmin = num
call psb_min(ctxt,nmin) call psb_min(ctxt,nmin)
if (nmin < szero) then if (nmin < szero) then
num = szero num = szero
den = sone den = sone
@ -487,25 +487,25 @@ contains
end if end if
prec%ag_data%op_complexity = num/den prec%ag_data%op_complexity = num/den
end subroutine amg_s_cmp_compl end subroutine amg_s_cmp_compl
! !
! Average coarsening ratio ! Average coarsening ratio
! !
function amg_s_get_avg_cr(prec) result(val) function amg_s_get_avg_cr(prec) result(val)
implicit none implicit none
class(amg_sprec_type), intent(in) :: prec class(amg_sprec_type), intent(in) :: prec
real(psb_spk_) :: val real(psb_spk_) :: val
val = prec%ag_data%avg_cr val = prec%ag_data%avg_cr
end function amg_s_get_avg_cr end function amg_s_get_avg_cr
subroutine amg_s_cmp_avg_cr(prec)
implicit none subroutine amg_s_cmp_avg_cr(prec)
implicit none
class(amg_sprec_type), intent(inout) :: prec class(amg_sprec_type), intent(inout) :: prec
real(psb_spk_) :: avgcr real(psb_spk_) :: avgcr
type(psb_ctxt_type) :: ctxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: il, nl, iam, np integer(psb_ipk_) :: il, nl, iam, np
@ -519,12 +519,12 @@ contains
do il=2,nl do il=2,nl
avgcr = avgcr + max(szero,prec%precv(il)%szratio) avgcr = avgcr + max(szero,prec%precv(il)%szratio)
end do end do
avgcr = avgcr / (nl-1) avgcr = avgcr / (nl-1)
end if end if
call psb_sum(ctxt,avgcr) call psb_sum(ctxt,avgcr)
prec%ag_data%avg_cr = avgcr/np prec%ag_data%avg_cr = avgcr/np
end subroutine amg_s_cmp_avg_cr end subroutine amg_s_cmp_avg_cr
! !
! Subroutines: amg_Tprec_free ! Subroutines: amg_Tprec_free
! Version: real ! Version: real
@ -538,74 +538,74 @@ contains
! error code. ! error code.
! !
subroutine amg_sprecfree(p,info) subroutine amg_sprecfree(p,info)
implicit none implicit none
! Arguments ! Arguments
type(amg_sprec_type), intent(inout) :: p type(amg_sprec_type), intent(inout) :: p
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
! Local variables ! Local variables
integer(psb_ipk_) :: me,err_act,i integer(psb_ipk_) :: me,err_act,i
character(len=20) :: name character(len=20) :: name
info=psb_success_ info=psb_success_
name = 'amg_sprecfree' name = 'amg_sprecfree'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; return info = psb_err_internal_error_; return
end if end if
me=-1 me=-1
call p%free(info) call p%free(info)
return return
end subroutine amg_sprecfree end subroutine amg_sprecfree
subroutine amg_s_prec_free(prec,info) subroutine amg_s_prec_free(prec,info)
implicit none implicit none
! Arguments ! Arguments
class(amg_sprec_type), intent(inout) :: prec class(amg_sprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
! Local variables ! Local variables
integer(psb_ipk_) :: me,err_act,i integer(psb_ipk_) :: me,err_act,i
character(len=20) :: name character(len=20) :: name
info=psb_success_ info=psb_success_
name = 'amg_sprecfree' name = 'amg_sprecfree'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999 info = psb_err_internal_error_; goto 9999
end if end if
me=-1 me=-1
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
do i=1,size(prec%precv) do i=1,size(prec%precv)
call prec%precv(i)%free(info) call prec%precv(i)%free(info)
end do end do
deallocate(prec%precv,stat=info) deallocate(prec%precv,stat=info)
end if end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
9999 call psb_error_handler(err_act) 9999 call psb_error_handler(err_act)
return return
end subroutine amg_s_prec_free end subroutine amg_s_prec_free
! !
! Top level methods. ! Top level methods.
! !
subroutine amg_s_apply2_vect(prec,x,y,desc_data,info,trans,work) subroutine amg_s_apply2_vect(prec,x,y,desc_data,info,trans,work)
implicit none implicit none
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
class(amg_sprec_type), intent(inout) :: prec class(amg_sprec_type), intent(inout) :: prec
type(psb_s_vect_type),intent(inout) :: x type(psb_s_vect_type),intent(inout) :: x
@ -618,13 +618,13 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
select type(prec) select type(prec)
type is (amg_sprec_type) type is (amg_sprec_type)
call amg_precapply(prec,x,y,desc_data,info,trans,work) call amg_precapply(prec,x,y,desc_data,info,trans,work)
class default class default
info = psb_err_missing_override_method_ info = psb_err_missing_override_method_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end select end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -636,7 +636,7 @@ contains
end subroutine amg_s_apply2_vect end subroutine amg_s_apply2_vect
subroutine amg_s_apply1_vect(prec,x,desc_data,info,trans,work) subroutine amg_s_apply1_vect(prec,x,desc_data,info,trans,work)
implicit none implicit none
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
class(amg_sprec_type), intent(inout) :: prec class(amg_sprec_type), intent(inout) :: prec
type(psb_s_vect_type),intent(inout) :: x type(psb_s_vect_type),intent(inout) :: x
@ -648,13 +648,13 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
select type(prec) select type(prec)
type is (amg_sprec_type) type is (amg_sprec_type)
call amg_precapply(prec,x,desc_data,info,trans,work) call amg_precapply(prec,x,desc_data,info,trans,work)
class default class default
info = psb_err_missing_override_method_ info = psb_err_missing_override_method_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end select end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -667,7 +667,7 @@ contains
subroutine amg_s_apply2v(prec,x,y,desc_data,info,trans,work) subroutine amg_s_apply2v(prec,x,y,desc_data,info,trans,work)
implicit none implicit none
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
class(amg_sprec_type), intent(inout) :: prec class(amg_sprec_type), intent(inout) :: prec
real(psb_spk_),intent(inout) :: x(:) real(psb_spk_),intent(inout) :: x(:)
@ -680,13 +680,13 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
select type(prec) select type(prec)
type is (amg_sprec_type) type is (amg_sprec_type)
call amg_precapply(prec,x,y,desc_data,info,trans,work) call amg_precapply(prec,x,y,desc_data,info,trans,work)
class default class default
info = psb_err_missing_override_method_ info = psb_err_missing_override_method_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end select end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -698,7 +698,7 @@ contains
end subroutine amg_s_apply2v end subroutine amg_s_apply2v
subroutine amg_s_apply1v(prec,x,desc_data,info,trans) subroutine amg_s_apply1v(prec,x,desc_data,info,trans)
implicit none implicit none
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
class(amg_sprec_type), intent(inout) :: prec class(amg_sprec_type), intent(inout) :: prec
real(psb_spk_),intent(inout) :: x(:) real(psb_spk_),intent(inout) :: x(:)
@ -709,13 +709,13 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
select type(prec) select type(prec)
type is (amg_sprec_type) type is (amg_sprec_type)
call amg_precapply(prec,x,desc_data,info,trans) call amg_precapply(prec,x,desc_data,info,trans)
class default class default
info = psb_err_missing_override_method_ info = psb_err_missing_override_method_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end select end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -730,8 +730,8 @@ contains
subroutine amg_s_dump(prec,info,istart,iend,iproc,prefix,head,& subroutine amg_s_dump(prec,info,istart,iend,iproc,prefix,head,&
& ac,rp,smoother,solver,tprol,& & ac,rp,smoother,solver,tprol,&
& global_num) & global_num)
implicit none implicit none
class(amg_sprec_type), intent(in) :: prec class(amg_sprec_type), intent(in) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: istart, iend, iproc integer(psb_ipk_), intent(in), optional :: istart, iend, iproc
@ -742,27 +742,27 @@ contains
integer(psb_ipk_) :: iam, np, iproc_ integer(psb_ipk_) :: iam, np, iproc_
character(len=80) :: prefix_ character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than character(len=120) :: fname ! len should be at least 20 more than
! len of prefix_ ! len of prefix_
info = 0 info = 0
icontxt = prec%ctxt icontxt = prec%ctxt
call psb_info(icontxt,iam,np) call psb_info(icontxt,iam,np)
iln = size(prec%precv) iln = size(prec%precv)
if (present(istart)) then if (present(istart)) then
il1 = max(1,istart) il1 = max(1,istart)
else else
il1 = min(2,iln) il1 = min(2,iln)
end if end if
if (present(iend)) then if (present(iend)) then
iln = min(iln, iend) iln = min(iln, iend)
end if end if
iproc_ = -1 iproc_ = -1
if (present(iproc)) then if (present(iproc)) then
iproc_ = iproc iproc_ = iproc
end if end if
if ((iproc_ == -1).or.(iproc_==iam)) then if ((iproc_ == -1).or.(iproc_==iam)) then
do lev=il1, iln do lev=il1, iln
call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,& call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,&
& ac=ac,smoother=smoother,solver=solver,rp=rp,tprol=tprol, & & ac=ac,smoother=smoother,solver=solver,rp=rp,tprol=tprol, &
@ -773,7 +773,7 @@ contains
subroutine amg_s_cnv(prec,info,amold,vmold,imold) subroutine amg_s_cnv(prec,info,amold,vmold,imold)
implicit none implicit none
class(amg_sprec_type), intent(inout) :: prec class(amg_sprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
class(psb_s_base_sparse_mat), intent(in), optional :: amold class(psb_s_base_sparse_mat), intent(in), optional :: amold
@ -781,7 +781,7 @@ contains
class(psb_i_base_vect_type), intent(in), optional :: imold class(psb_i_base_vect_type), intent(in), optional :: imold
integer(psb_ipk_) :: i integer(psb_ipk_) :: i
info = psb_success_ info = psb_success_
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
do i=1,size(prec%precv) do i=1,size(prec%precv)
@ -789,24 +789,24 @@ contains
& call prec%precv(i)%cnv(info,amold=amold,vmold=vmold,imold=imold) & call prec%precv(i)%cnv(info,amold=amold,vmold=vmold,imold=imold)
end do end do
end if end if
end subroutine amg_s_cnv end subroutine amg_s_cnv
subroutine amg_s_clone(prec,precout,info) subroutine amg_s_clone(prec,precout,info)
implicit none implicit none
class(amg_sprec_type), intent(inout) :: prec class(amg_sprec_type), intent(inout) :: prec
class(psb_sprec_type), intent(inout) :: precout class(psb_sprec_type), intent(inout) :: precout
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
call precout%free(info) call precout%free(info)
if (info == 0) call amg_s_inner_clone(prec,precout,info) if (info == 0) call amg_s_inner_clone(prec,precout,info)
end subroutine amg_s_clone end subroutine amg_s_clone
subroutine amg_s_inner_clone(prec,precout,info) subroutine amg_s_inner_clone(prec,precout,info)
implicit none implicit none
class(amg_sprec_type), intent(inout) :: prec class(amg_sprec_type), intent(inout) :: prec
class(psb_sprec_type), target, intent(inout) :: precout class(psb_sprec_type), target, intent(inout) :: precout
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -821,17 +821,17 @@ contains
pout%ctxt = prec%ctxt pout%ctxt = prec%ctxt
pout%ag_data = prec%ag_data pout%ag_data = prec%ag_data
pout%outer_sweeps = prec%outer_sweeps pout%outer_sweeps = prec%outer_sweeps
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
ln = size(prec%precv) ln = size(prec%precv)
allocate(pout%precv(ln),stat=info) allocate(pout%precv(ln),stat=info)
if (info /= psb_success_) goto 9999 if (info /= psb_success_) goto 9999
if (ln >= 1) then if (ln >= 1) then
call prec%precv(1)%clone(pout%precv(1),info) call prec%precv(1)%clone(pout%precv(1),info)
end if end if
do lev=2, ln do lev=2, ln
if (info /= psb_success_) exit if (info /= psb_success_) exit
call prec%precv(lev)%clone(pout%precv(lev),info) call prec%precv(lev)%clone(pout%precv(lev),info)
if (info == psb_success_) then if (info == psb_success_) then
pout%precv(lev)%base_a => pout%precv(lev)%ac pout%precv(lev)%base_a => pout%precv(lev)%ac
pout%precv(lev)%base_desc => pout%precv(lev)%desc_ac pout%precv(lev)%base_desc => pout%precv(lev)%desc_ac
pout%precv(lev)%map%p_desc_U => pout%precv(lev-1)%base_desc pout%precv(lev)%map%p_desc_U => pout%precv(lev-1)%base_desc
@ -842,7 +842,7 @@ contains
if (allocated(prec%precv(1)%wrk)) & if (allocated(prec%precv(1)%wrk)) &
& call pout%allocate_wrk(info,vmold=prec%precv(1)%wrk%vx2l%v) & call pout%allocate_wrk(info,vmold=prec%precv(1)%wrk%vx2l%v)
class default class default
write(0,*) 'Error: wrong out type' write(0,*) 'Error: wrong out type'
info = psb_err_invalid_input_ info = psb_err_invalid_input_
end select end select
@ -854,14 +854,14 @@ contains
implicit none implicit none
class(amg_sprec_type), intent(inout) :: prec class(amg_sprec_type), intent(inout) :: prec
class(amg_sprec_type), intent(inout), target :: b class(amg_sprec_type), intent(inout), target :: b
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i integer(psb_ipk_) :: i
if (same_type_as(prec,b)) then if (same_type_as(prec,b)) then
if (allocated(b%precv)) then if (allocated(b%precv)) then
! This might not be required if FINAL procedures are available. ! This might not be required if FINAL procedures are available.
call b%free(info) call b%free(info)
if (info /= psb_success_) then if (info /= psb_success_) then
!????? !?????
!!$ return !!$ return
endif endif
@ -869,7 +869,7 @@ contains
b%ctxt = prec%ctxt b%ctxt = prec%ctxt
b%ag_data = prec%ag_data b%ag_data = prec%ag_data
b%outer_sweeps = prec%outer_sweeps b%outer_sweeps = prec%outer_sweeps
call move_alloc(prec%precv,b%precv) call move_alloc(prec%precv,b%precv)
! Fix the pointers except on level 1. ! Fix the pointers except on level 1.
do i=2, size(b%precv) do i=2, size(b%precv)
@ -878,7 +878,7 @@ contains
b%precv(i)%map%p_desc_U => b%precv(i-1)%base_desc b%precv(i)%map%p_desc_U => b%precv(i-1)%base_desc
b%precv(i)%map%p_desc_V => b%precv(i)%base_desc b%precv(i)%map%p_desc_V => b%precv(i)%base_desc
end do end do
else else
write(0,*) 'Warning: PREC%move_alloc onto different type?' write(0,*) 'Warning: PREC%move_alloc onto different type?'
info = psb_err_internal_error_ info = psb_err_internal_error_
@ -888,7 +888,7 @@ contains
subroutine amg_s_allocate_wrk(prec,info,vmold,desc) subroutine amg_s_allocate_wrk(prec,info,vmold,desc)
use psb_base_mod use psb_base_mod
implicit none implicit none
! Arguments ! Arguments
class(amg_sprec_type), intent(inout) :: prec class(amg_sprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -896,37 +896,37 @@ contains
! !
! In MLD the DESC optional argument is ignored, since ! In MLD the DESC optional argument is ignored, since
! the necessary info is contained in the various entries of the ! the necessary info is contained in the various entries of the
! PRECV component. ! PRECV component.
type(psb_desc_type), intent(in), optional :: desc type(psb_desc_type), intent(in), optional :: desc
! Local variables ! Local variables
integer(psb_ipk_) :: me,err_act,i,j,level,nlev, nc2l integer(psb_ipk_) :: me,err_act,i,j,level,nlev, nc2l
character(len=20) :: name character(len=20) :: name
info=psb_success_ info=psb_success_
name = 'amg_s_allocate_wrk' name = 'amg_s_allocate_wrk'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999 info = psb_err_internal_error_; goto 9999
end if end if
nlev = size(prec%precv) nlev = size(prec%precv)
level = 1 level = 1
do level = 1, nlev do level = 1, nlev
call prec%precv(level)%allocate_wrk(info,vmold=vmold) call prec%precv(level)%allocate_wrk(info,vmold=vmold)
if (psb_errstatus_fatal()) then if (psb_errstatus_fatal()) then
nc2l = prec%precv(level)%base_desc%get_local_cols() nc2l = prec%precv(level)%base_desc%get_local_cols()
info=psb_err_alloc_request_ info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l/), a_err='real(psb_spk_)') call psb_errpush(info,name,i_err=(/2*nc2l/), a_err='real(psb_spk_)')
goto 9999 goto 9999
end if end if
end do end do
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
9999 call psb_error_handler(err_act) 9999 call psb_error_handler(err_act)
return return
end subroutine amg_s_allocate_wrk end subroutine amg_s_allocate_wrk
subroutine amg_s_free_wrk(prec,info) subroutine amg_s_free_wrk(prec,info)
@ -948,13 +948,13 @@ contains
info = psb_err_internal_error_; goto 9999 info = psb_err_internal_error_; goto 9999
end if end if
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
nlev = size(prec%precv) nlev = size(prec%precv)
do level = 1, nlev do level = 1, nlev
call prec%precv(level)%free_wrk(info) call prec%precv(level)%free_wrk(info)
end do end do
end if end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -966,7 +966,7 @@ contains
function amg_s_is_allocated_wrk(prec) result(res) function amg_s_is_allocated_wrk(prec) result(res)
use psb_base_mod use psb_base_mod
implicit none implicit none
! Arguments ! Arguments
class(amg_sprec_type), intent(in) :: prec class(amg_sprec_type), intent(in) :: prec
logical :: res logical :: res
@ -974,7 +974,7 @@ contains
res = .false. res = .false.
if (.not.allocated(prec%precv)) return if (.not.allocated(prec%precv)) return
res = allocated(prec%precv(1)%wrk) res = allocated(prec%precv(1)%wrk)
end function amg_s_is_allocated_wrk end function amg_s_is_allocated_wrk
end module amg_s_prec_type end module amg_s_prec_type

292
amgprec/amg_z_prec_type.f90
Unescape Escape View File

@ -1,15 +1,15 @@
! !
! !
! AMG4PSBLAS version 1.0 ! AMG4PSBLAS version 1.0
! Algebraic Multigrid Package ! Algebraic Multigrid Package
! based on PSBLAS (Parallel Sparse BLAS version 3.5) ! based on PSBLAS (Parallel Sparse BLAS version 3.5)
! !
! (C) Copyright 2020 ! (C) Copyright 2020
! !
! Salvatore Filippone ! Salvatore Filippone
! Pasqua D'Ambra ! Pasqua D'Ambra
! Fabio Durastante ! Fabio Durastante
! !
! Redistribution and use in source and binary forms, with or without ! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions ! modification, are permitted provided that the following conditions
! are met: ! are met:
@ -21,7 +21,7 @@
! 3. The name of the AMG4PSBLAS group or the names of its contributors may ! 3. The name of the AMG4PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this ! not be used to endorse or promote products derived from this
! software without specific written permission. ! software without specific written permission.
! !
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
@ -33,21 +33,21 @@
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! 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 ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE. ! POSSIBILITY OF SUCH DAMAGE.
! !
! !
! File: amg_z_prec_type.f90 ! File: amg_z_prec_type.f90
! !
! Module: amg_z_prec_type ! Module: amg_z_prec_type
! !
! This module defines: ! This module defines:
! - the amg_z_prec_type data structure containing the preconditioner and related ! - the amg_z_prec_type data structure containing the preconditioner and related
! data structures; ! data structures;
! !
! It contains routines for ! It contains routines for
! - Building and applying; ! - Building and applying;
! - checking if the preconditioner is correctly defined; ! - checking if the preconditioner is correctly defined;
! - printing a description of the preconditioner; ! - printing a description of the preconditioner;
! - deallocating the preconditioner data structure. ! - deallocating the preconditioner data structure.
! !
module amg_z_prec_type module amg_z_prec_type
@ -70,25 +70,25 @@ module amg_z_prec_type
! It consists of an array of 'one-level' intermediate data structures ! It consists of an array of 'one-level' intermediate data structures
! of type amg_zonelev_type, each containing the information needed to apply ! of type amg_zonelev_type, each containing the information needed to apply
! the smoothing and the coarse-space correction at a generic level. RT is the ! the smoothing and the coarse-space correction at a generic level. RT is the
! real data type, i.e. S for both S and C, and D for both D and Z. ! real data type, i.e. S for both S and C, and D for both D and Z.
! !
! type amg_zprec_type ! type amg_zprec_type
! type(amg_zonelev_type), allocatable :: precv(:) ! type(amg_zonelev_type), allocatable :: precv(:)
! end type amg_zprec_type ! end type amg_zprec_type
! !
! Note that the levels are numbered in increasing order starting from ! Note that the levels are numbered in increasing order starting from
! the level 1 as the finest one, and the number of levels is given by ! the level 1 as the finest one, and the number of levels is given by
! size(precv(:)) which is the id of the coarsest level. ! size(precv(:)) which is the id of the coarsest level.
! In the multigrid literature many authors number the levels in the opposite ! In the multigrid literature many authors number the levels in the opposite
! order, with level 0 being the id of the coarsest level. ! order, with level 0 being the id of the coarsest level.
! !
! !
integer, parameter, private :: wv_size_=4 integer, parameter, private :: wv_size_=4
type, extends(psb_zprec_type) :: amg_zprec_type type, extends(psb_zprec_type) :: amg_zprec_type
type(amg_daggr_data) :: ag_data type(amg_daggr_data) :: ag_data
! !
! Number of outer sweeps. Sometimes 2 V-cycles may be better than 1 W-cycle. ! Number of outer sweeps. Sometimes 2 V-cycles may be better than 1 W-cycle.
! !
integer(psb_ipk_) :: outer_sweeps = 1 integer(psb_ipk_) :: outer_sweeps = 1
! !
@ -97,11 +97,11 @@ module amg_z_prec_type
! to keep track against what is put later in the multilevel array ! to keep track against what is put later in the multilevel array
! !
integer(psb_ipk_) :: coarse_solver = -1 integer(psb_ipk_) :: coarse_solver = -1
! !
! The multilevel hierarchy ! The multilevel hierarchy
! !
type(amg_z_onelev_type), allocatable :: precv(:) type(amg_z_onelev_type), allocatable :: precv(:)
contains contains
procedure, pass(prec) :: psb_z_apply2_vect => amg_z_apply2_vect procedure, pass(prec) :: psb_z_apply2_vect => amg_z_apply2_vect
procedure, pass(prec) :: psb_z_apply1_vect => amg_z_apply1_vect procedure, pass(prec) :: psb_z_apply1_vect => amg_z_apply1_vect
@ -127,7 +127,7 @@ module amg_z_prec_type
procedure, pass(prec) :: cseti => amg_zcprecseti procedure, pass(prec) :: cseti => amg_zcprecseti
procedure, pass(prec) :: csetc => amg_zcprecsetc procedure, pass(prec) :: csetc => amg_zcprecsetc
procedure, pass(prec) :: csetr => amg_zcprecsetr procedure, pass(prec) :: csetr => amg_zcprecsetr
generic, public :: set => cseti, csetc, csetr, setsm, setsv, setag generic, public :: set => setsm, setsv, setag
procedure, pass(prec) :: get_smoother => amg_z_get_smootherp procedure, pass(prec) :: get_smoother => amg_z_get_smootherp
procedure, pass(prec) :: get_solver => amg_z_get_solverp procedure, pass(prec) :: get_solver => amg_z_get_solverp
procedure, pass(prec) :: move_alloc => z_prec_move_alloc procedure, pass(prec) :: move_alloc => z_prec_move_alloc
@ -157,7 +157,7 @@ module amg_z_prec_type
interface amg_precdescr interface amg_precdescr
subroutine amg_zfile_prec_descr(prec,iout,root) subroutine amg_zfile_prec_descr(prec,iout,root)
import :: amg_zprec_type, psb_ipk_ import :: amg_zprec_type, psb_ipk_
implicit none implicit none
! Arguments ! Arguments
class(amg_zprec_type), intent(in) :: prec class(amg_zprec_type), intent(in) :: prec
integer(psb_ipk_), intent(in), optional :: iout integer(psb_ipk_), intent(in), optional :: iout
@ -211,7 +211,7 @@ module amg_z_prec_type
end subroutine amg_zprecaply1 end subroutine amg_zprecaply1
end interface end interface
interface interface
subroutine amg_zprecsetsm(prec,val,info,ilev,ilmax,pos) subroutine amg_zprecsetsm(prec,val,info,ilev,ilmax,pos)
import :: psb_zspmat_type, psb_desc_type, psb_dpk_, & import :: psb_zspmat_type, psb_desc_type, psb_dpk_, &
& amg_zprec_type, amg_z_base_smoother_type, psb_ipk_ & amg_zprec_type, amg_z_base_smoother_type, psb_ipk_
@ -243,7 +243,7 @@ module amg_z_prec_type
import :: psb_zspmat_type, psb_desc_type, psb_dpk_, & import :: psb_zspmat_type, psb_desc_type, psb_dpk_, &
& amg_zprec_type, psb_ipk_ & amg_zprec_type, psb_ipk_
class(amg_zprec_type), intent(inout) :: prec class(amg_zprec_type), intent(inout) :: prec
character(len=*), intent(in) :: what character(len=*), intent(in) :: what
integer(psb_ipk_), intent(in) :: val integer(psb_ipk_), intent(in) :: val
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx
@ -253,7 +253,7 @@ module amg_z_prec_type
import :: psb_zspmat_type, psb_desc_type, psb_dpk_, & import :: psb_zspmat_type, psb_desc_type, psb_dpk_, &
& amg_zprec_type, psb_ipk_ & amg_zprec_type, psb_ipk_
class(amg_zprec_type), intent(inout) :: prec class(amg_zprec_type), intent(inout) :: prec
character(len=*), intent(in) :: what character(len=*), intent(in) :: what
real(psb_dpk_), intent(in) :: val real(psb_dpk_), intent(in) :: val
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx
@ -263,7 +263,7 @@ module amg_z_prec_type
import :: psb_zspmat_type, psb_desc_type, psb_dpk_, & import :: psb_zspmat_type, psb_desc_type, psb_dpk_, &
& amg_zprec_type, psb_ipk_ & amg_zprec_type, psb_ipk_
class(amg_zprec_type), intent(inout) :: prec class(amg_zprec_type), intent(inout) :: prec
character(len=*), intent(in) :: what character(len=*), intent(in) :: what
character(len=*), intent(in) :: string character(len=*), intent(in) :: string
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx
@ -341,28 +341,28 @@ module amg_z_prec_type
! character, intent(in),optional :: upd ! character, intent(in),optional :: upd
end subroutine amg_z_smoothers_bld end subroutine amg_z_smoothers_bld
end interface amg_smoothers_bld end interface amg_smoothers_bld
contains contains
! !
! Function returning a pointer to the smoother ! Function returning a pointer to the smoother
! !
function amg_z_get_smootherp(prec,ilev) result(val) function amg_z_get_smootherp(prec,ilev) result(val)
implicit none implicit none
class(amg_zprec_type), target, intent(in) :: prec class(amg_zprec_type), target, intent(in) :: prec
integer(psb_ipk_), optional :: ilev integer(psb_ipk_), optional :: ilev
class(amg_z_base_smoother_type), pointer :: val class(amg_z_base_smoother_type), pointer :: val
integer(psb_ipk_) :: ilev_ integer(psb_ipk_) :: ilev_
val => null() val => null()
if (present(ilev)) then if (present(ilev)) then
ilev_ = ilev ilev_ = ilev
else else
! What is a good default? ! What is a good default?
ilev_ = 1 ilev_ = 1
end if end if
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
if ((1<=ilev_).and.(ilev_<=size(prec%precv))) then if ((1<=ilev_).and.(ilev_<=size(prec%precv))) then
if (allocated(prec%precv(ilev_)%sm)) then if (allocated(prec%precv(ilev_)%sm)) then
val => prec%precv(ilev_)%sm val => prec%precv(ilev_)%sm
end if end if
end if end if
@ -372,23 +372,23 @@ contains
! Function returning a pointer to the solver ! Function returning a pointer to the solver
! !
function amg_z_get_solverp(prec,ilev) result(val) function amg_z_get_solverp(prec,ilev) result(val)
implicit none implicit none
class(amg_zprec_type), target, intent(in) :: prec class(amg_zprec_type), target, intent(in) :: prec
integer(psb_ipk_), optional :: ilev integer(psb_ipk_), optional :: ilev
class(amg_z_base_solver_type), pointer :: val class(amg_z_base_solver_type), pointer :: val
integer(psb_ipk_) :: ilev_ integer(psb_ipk_) :: ilev_
val => null() val => null()
if (present(ilev)) then if (present(ilev)) then
ilev_ = ilev ilev_ = ilev
else else
! What is a good default? ! What is a good default?
ilev_ = 1 ilev_ = 1
end if end if
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
if ((1<=ilev_).and.(ilev_<=size(prec%precv))) then if ((1<=ilev_).and.(ilev_<=size(prec%precv))) then
if (allocated(prec%precv(ilev_)%sm)) then if (allocated(prec%precv(ilev_)%sm)) then
if (allocated(prec%precv(ilev_)%sm%sv)) then if (allocated(prec%precv(ilev_)%sm%sv)) then
val => prec%precv(ilev_)%sm%sv val => prec%precv(ilev_)%sm%sv
end if end if
end if end if
@ -399,25 +399,25 @@ contains
! Function returning the size of the precv(:) array ! Function returning the size of the precv(:) array
! !
function amg_z_get_nlevs(prec) result(val) function amg_z_get_nlevs(prec) result(val)
implicit none implicit none
class(amg_zprec_type), intent(in) :: prec class(amg_zprec_type), intent(in) :: prec
integer(psb_ipk_) :: val integer(psb_ipk_) :: val
val = 0 val = 0
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
val = size(prec%precv) val = size(prec%precv)
end if end if
end function amg_z_get_nlevs end function amg_z_get_nlevs
! !
! Function returning the size of the amg_prec_type data structure ! Function returning the size of the amg_prec_type data structure
! in bytes or in number of nonzeros of the operator(s) involved. ! in bytes or in number of nonzeros of the operator(s) involved.
! !
function amg_z_get_nzeros(prec) result(val) function amg_z_get_nzeros(prec) result(val)
implicit none implicit none
class(amg_zprec_type), intent(in) :: prec class(amg_zprec_type), intent(in) :: prec
integer(psb_epk_) :: val integer(psb_epk_) :: val
integer(psb_ipk_) :: i integer(psb_ipk_) :: i
val = 0 val = 0
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
do i=1, size(prec%precv) do i=1, size(prec%precv)
val = val + prec%precv(i)%get_nzeros() val = val + prec%precv(i)%get_nzeros()
end do end do
@ -425,13 +425,13 @@ contains
end function amg_z_get_nzeros end function amg_z_get_nzeros
function amg_zprec_sizeof(prec) result(val) function amg_zprec_sizeof(prec) result(val)
implicit none implicit none
class(amg_zprec_type), intent(in) :: prec class(amg_zprec_type), intent(in) :: prec
integer(psb_epk_) :: val integer(psb_epk_) :: val
integer(psb_ipk_) :: i integer(psb_ipk_) :: i
val = 0 val = 0
val = val + psb_sizeof_ip val = val + psb_sizeof_ip
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
do i=1, size(prec%precv) do i=1, size(prec%precv)
val = val + prec%precv(i)%sizeof() val = val + prec%precv(i)%sizeof()
end do end do
@ -444,40 +444,40 @@ contains
! various level to the nonzeroes at the fine level ! various level to the nonzeroes at the fine level
! (original matrix) ! (original matrix)
! !
function amg_z_get_compl(prec) result(val) function amg_z_get_compl(prec) result(val)
implicit none implicit none
class(amg_zprec_type), intent(in) :: prec class(amg_zprec_type), intent(in) :: prec
complex(psb_dpk_) :: val complex(psb_dpk_) :: val
val = prec%ag_data%op_complexity val = prec%ag_data%op_complexity
end function amg_z_get_compl end function amg_z_get_compl
subroutine amg_z_cmp_compl(prec)
implicit none subroutine amg_z_cmp_compl(prec)
implicit none
class(amg_zprec_type), intent(inout) :: prec class(amg_zprec_type), intent(inout) :: prec
real(psb_dpk_) :: num, den, nmin real(psb_dpk_) :: num, den, nmin
type(psb_ctxt_type) :: ctxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: il integer(psb_ipk_) :: il
num = -done num = -done
den = done den = done
ctxt = prec%ctxt ctxt = prec%ctxt
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
il = 1 il = 1
num = prec%precv(il)%base_a%get_nzeros() num = prec%precv(il)%base_a%get_nzeros()
if (num >= dzero) then if (num >= dzero) then
den = num den = num
do il=2,size(prec%precv) do il=2,size(prec%precv)
num = num + max(0,prec%precv(il)%base_a%get_nzeros()) num = num + max(0,prec%precv(il)%base_a%get_nzeros())
end do end do
end if end if
end if end if
nmin = num nmin = num
call psb_min(ctxt,nmin) call psb_min(ctxt,nmin)
if (nmin < dzero) then if (nmin < dzero) then
num = dzero num = dzero
den = done den = done
@ -487,25 +487,25 @@ contains
end if end if
prec%ag_data%op_complexity = num/den prec%ag_data%op_complexity = num/den
end subroutine amg_z_cmp_compl end subroutine amg_z_cmp_compl
! !
! Average coarsening ratio ! Average coarsening ratio
! !
function amg_z_get_avg_cr(prec) result(val) function amg_z_get_avg_cr(prec) result(val)
implicit none implicit none
class(amg_zprec_type), intent(in) :: prec class(amg_zprec_type), intent(in) :: prec
complex(psb_dpk_) :: val complex(psb_dpk_) :: val
val = prec%ag_data%avg_cr val = prec%ag_data%avg_cr
end function amg_z_get_avg_cr end function amg_z_get_avg_cr
subroutine amg_z_cmp_avg_cr(prec)
implicit none subroutine amg_z_cmp_avg_cr(prec)
implicit none
class(amg_zprec_type), intent(inout) :: prec class(amg_zprec_type), intent(inout) :: prec
real(psb_dpk_) :: avgcr real(psb_dpk_) :: avgcr
type(psb_ctxt_type) :: ctxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: il, nl, iam, np integer(psb_ipk_) :: il, nl, iam, np
@ -519,12 +519,12 @@ contains
do il=2,nl do il=2,nl
avgcr = avgcr + max(dzero,prec%precv(il)%szratio) avgcr = avgcr + max(dzero,prec%precv(il)%szratio)
end do end do
avgcr = avgcr / (nl-1) avgcr = avgcr / (nl-1)
end if end if
call psb_sum(ctxt,avgcr) call psb_sum(ctxt,avgcr)
prec%ag_data%avg_cr = avgcr/np prec%ag_data%avg_cr = avgcr/np
end subroutine amg_z_cmp_avg_cr end subroutine amg_z_cmp_avg_cr
! !
! Subroutines: amg_Tprec_free ! Subroutines: amg_Tprec_free
! Version: complex ! Version: complex
@ -538,74 +538,74 @@ contains
! error code. ! error code.
! !
subroutine amg_zprecfree(p,info) subroutine amg_zprecfree(p,info)
implicit none implicit none
! Arguments ! Arguments
type(amg_zprec_type), intent(inout) :: p type(amg_zprec_type), intent(inout) :: p
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
! Local variables ! Local variables
integer(psb_ipk_) :: me,err_act,i integer(psb_ipk_) :: me,err_act,i
character(len=20) :: name character(len=20) :: name
info=psb_success_ info=psb_success_
name = 'amg_zprecfree' name = 'amg_zprecfree'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; return info = psb_err_internal_error_; return
end if end if
me=-1 me=-1
call p%free(info) call p%free(info)
return return
end subroutine amg_zprecfree end subroutine amg_zprecfree
subroutine amg_z_prec_free(prec,info) subroutine amg_z_prec_free(prec,info)
implicit none implicit none
! Arguments ! Arguments
class(amg_zprec_type), intent(inout) :: prec class(amg_zprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
! Local variables ! Local variables
integer(psb_ipk_) :: me,err_act,i integer(psb_ipk_) :: me,err_act,i
character(len=20) :: name character(len=20) :: name
info=psb_success_ info=psb_success_
name = 'amg_zprecfree' name = 'amg_zprecfree'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999 info = psb_err_internal_error_; goto 9999
end if end if
me=-1 me=-1
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
do i=1,size(prec%precv) do i=1,size(prec%precv)
call prec%precv(i)%free(info) call prec%precv(i)%free(info)
end do end do
deallocate(prec%precv,stat=info) deallocate(prec%precv,stat=info)
end if end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
9999 call psb_error_handler(err_act) 9999 call psb_error_handler(err_act)
return return
end subroutine amg_z_prec_free end subroutine amg_z_prec_free
! !
! Top level methods. ! Top level methods.
! !
subroutine amg_z_apply2_vect(prec,x,y,desc_data,info,trans,work) subroutine amg_z_apply2_vect(prec,x,y,desc_data,info,trans,work)
implicit none implicit none
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
class(amg_zprec_type), intent(inout) :: prec class(amg_zprec_type), intent(inout) :: prec
type(psb_z_vect_type),intent(inout) :: x type(psb_z_vect_type),intent(inout) :: x
@ -618,13 +618,13 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
select type(prec) select type(prec)
type is (amg_zprec_type) type is (amg_zprec_type)
call amg_precapply(prec,x,y,desc_data,info,trans,work) call amg_precapply(prec,x,y,desc_data,info,trans,work)
class default class default
info = psb_err_missing_override_method_ info = psb_err_missing_override_method_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end select end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -636,7 +636,7 @@ contains
end subroutine amg_z_apply2_vect end subroutine amg_z_apply2_vect
subroutine amg_z_apply1_vect(prec,x,desc_data,info,trans,work) subroutine amg_z_apply1_vect(prec,x,desc_data,info,trans,work)
implicit none implicit none
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
class(amg_zprec_type), intent(inout) :: prec class(amg_zprec_type), intent(inout) :: prec
type(psb_z_vect_type),intent(inout) :: x type(psb_z_vect_type),intent(inout) :: x
@ -648,13 +648,13 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
select type(prec) select type(prec)
type is (amg_zprec_type) type is (amg_zprec_type)
call amg_precapply(prec,x,desc_data,info,trans,work) call amg_precapply(prec,x,desc_data,info,trans,work)
class default class default
info = psb_err_missing_override_method_ info = psb_err_missing_override_method_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end select end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -667,7 +667,7 @@ contains
subroutine amg_z_apply2v(prec,x,y,desc_data,info,trans,work) subroutine amg_z_apply2v(prec,x,y,desc_data,info,trans,work)
implicit none implicit none
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
class(amg_zprec_type), intent(inout) :: prec class(amg_zprec_type), intent(inout) :: prec
complex(psb_dpk_),intent(inout) :: x(:) complex(psb_dpk_),intent(inout) :: x(:)
@ -680,13 +680,13 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
select type(prec) select type(prec)
type is (amg_zprec_type) type is (amg_zprec_type)
call amg_precapply(prec,x,y,desc_data,info,trans,work) call amg_precapply(prec,x,y,desc_data,info,trans,work)
class default class default
info = psb_err_missing_override_method_ info = psb_err_missing_override_method_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end select end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -698,7 +698,7 @@ contains
end subroutine amg_z_apply2v end subroutine amg_z_apply2v
subroutine amg_z_apply1v(prec,x,desc_data,info,trans) subroutine amg_z_apply1v(prec,x,desc_data,info,trans)
implicit none implicit none
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
class(amg_zprec_type), intent(inout) :: prec class(amg_zprec_type), intent(inout) :: prec
complex(psb_dpk_),intent(inout) :: x(:) complex(psb_dpk_),intent(inout) :: x(:)
@ -709,13 +709,13 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
select type(prec) select type(prec)
type is (amg_zprec_type) type is (amg_zprec_type)
call amg_precapply(prec,x,desc_data,info,trans) call amg_precapply(prec,x,desc_data,info,trans)
class default class default
info = psb_err_missing_override_method_ info = psb_err_missing_override_method_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end select end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -730,8 +730,8 @@ contains
subroutine amg_z_dump(prec,info,istart,iend,iproc,prefix,head,& subroutine amg_z_dump(prec,info,istart,iend,iproc,prefix,head,&
& ac,rp,smoother,solver,tprol,& & ac,rp,smoother,solver,tprol,&
& global_num) & global_num)
implicit none implicit none
class(amg_zprec_type), intent(in) :: prec class(amg_zprec_type), intent(in) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: istart, iend, iproc integer(psb_ipk_), intent(in), optional :: istart, iend, iproc
@ -742,27 +742,27 @@ contains
integer(psb_ipk_) :: iam, np, iproc_ integer(psb_ipk_) :: iam, np, iproc_
character(len=80) :: prefix_ character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than character(len=120) :: fname ! len should be at least 20 more than
! len of prefix_ ! len of prefix_
info = 0 info = 0
icontxt = prec%ctxt icontxt = prec%ctxt
call psb_info(icontxt,iam,np) call psb_info(icontxt,iam,np)
iln = size(prec%precv) iln = size(prec%precv)
if (present(istart)) then if (present(istart)) then
il1 = max(1,istart) il1 = max(1,istart)
else else
il1 = min(2,iln) il1 = min(2,iln)
end if end if
if (present(iend)) then if (present(iend)) then
iln = min(iln, iend) iln = min(iln, iend)
end if end if
iproc_ = -1 iproc_ = -1
if (present(iproc)) then if (present(iproc)) then
iproc_ = iproc iproc_ = iproc
end if end if
if ((iproc_ == -1).or.(iproc_==iam)) then if ((iproc_ == -1).or.(iproc_==iam)) then
do lev=il1, iln do lev=il1, iln
call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,& call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,&
& ac=ac,smoother=smoother,solver=solver,rp=rp,tprol=tprol, & & ac=ac,smoother=smoother,solver=solver,rp=rp,tprol=tprol, &
@ -773,7 +773,7 @@ contains
subroutine amg_z_cnv(prec,info,amold,vmold,imold) subroutine amg_z_cnv(prec,info,amold,vmold,imold)
implicit none implicit none
class(amg_zprec_type), intent(inout) :: prec class(amg_zprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
class(psb_z_base_sparse_mat), intent(in), optional :: amold class(psb_z_base_sparse_mat), intent(in), optional :: amold
@ -781,7 +781,7 @@ contains
class(psb_i_base_vect_type), intent(in), optional :: imold class(psb_i_base_vect_type), intent(in), optional :: imold
integer(psb_ipk_) :: i integer(psb_ipk_) :: i
info = psb_success_ info = psb_success_
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
do i=1,size(prec%precv) do i=1,size(prec%precv)
@ -789,24 +789,24 @@ contains
& call prec%precv(i)%cnv(info,amold=amold,vmold=vmold,imold=imold) & call prec%precv(i)%cnv(info,amold=amold,vmold=vmold,imold=imold)
end do end do
end if end if
end subroutine amg_z_cnv end subroutine amg_z_cnv
subroutine amg_z_clone(prec,precout,info) subroutine amg_z_clone(prec,precout,info)
implicit none implicit none
class(amg_zprec_type), intent(inout) :: prec class(amg_zprec_type), intent(inout) :: prec
class(psb_zprec_type), intent(inout) :: precout class(psb_zprec_type), intent(inout) :: precout
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
call precout%free(info) call precout%free(info)
if (info == 0) call amg_z_inner_clone(prec,precout,info) if (info == 0) call amg_z_inner_clone(prec,precout,info)
end subroutine amg_z_clone end subroutine amg_z_clone
subroutine amg_z_inner_clone(prec,precout,info) subroutine amg_z_inner_clone(prec,precout,info)
implicit none implicit none
class(amg_zprec_type), intent(inout) :: prec class(amg_zprec_type), intent(inout) :: prec
class(psb_zprec_type), target, intent(inout) :: precout class(psb_zprec_type), target, intent(inout) :: precout
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -821,17 +821,17 @@ contains
pout%ctxt = prec%ctxt pout%ctxt = prec%ctxt
pout%ag_data = prec%ag_data pout%ag_data = prec%ag_data
pout%outer_sweeps = prec%outer_sweeps pout%outer_sweeps = prec%outer_sweeps
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
ln = size(prec%precv) ln = size(prec%precv)
allocate(pout%precv(ln),stat=info) allocate(pout%precv(ln),stat=info)
if (info /= psb_success_) goto 9999 if (info /= psb_success_) goto 9999
if (ln >= 1) then if (ln >= 1) then
call prec%precv(1)%clone(pout%precv(1),info) call prec%precv(1)%clone(pout%precv(1),info)
end if end if
do lev=2, ln do lev=2, ln
if (info /= psb_success_) exit if (info /= psb_success_) exit
call prec%precv(lev)%clone(pout%precv(lev),info) call prec%precv(lev)%clone(pout%precv(lev),info)
if (info == psb_success_) then if (info == psb_success_) then
pout%precv(lev)%base_a => pout%precv(lev)%ac pout%precv(lev)%base_a => pout%precv(lev)%ac
pout%precv(lev)%base_desc => pout%precv(lev)%desc_ac pout%precv(lev)%base_desc => pout%precv(lev)%desc_ac
pout%precv(lev)%map%p_desc_U => pout%precv(lev-1)%base_desc pout%precv(lev)%map%p_desc_U => pout%precv(lev-1)%base_desc
@ -842,7 +842,7 @@ contains
if (allocated(prec%precv(1)%wrk)) & if (allocated(prec%precv(1)%wrk)) &
& call pout%allocate_wrk(info,vmold=prec%precv(1)%wrk%vx2l%v) & call pout%allocate_wrk(info,vmold=prec%precv(1)%wrk%vx2l%v)
class default class default
write(0,*) 'Error: wrong out type' write(0,*) 'Error: wrong out type'
info = psb_err_invalid_input_ info = psb_err_invalid_input_
end select end select
@ -854,14 +854,14 @@ contains
implicit none implicit none
class(amg_zprec_type), intent(inout) :: prec class(amg_zprec_type), intent(inout) :: prec
class(amg_zprec_type), intent(inout), target :: b class(amg_zprec_type), intent(inout), target :: b
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i integer(psb_ipk_) :: i
if (same_type_as(prec,b)) then if (same_type_as(prec,b)) then
if (allocated(b%precv)) then if (allocated(b%precv)) then
! This might not be required if FINAL procedures are available. ! This might not be required if FINAL procedures are available.
call b%free(info) call b%free(info)
if (info /= psb_success_) then if (info /= psb_success_) then
!????? !?????
!!$ return !!$ return
endif endif
@ -869,7 +869,7 @@ contains
b%ctxt = prec%ctxt b%ctxt = prec%ctxt
b%ag_data = prec%ag_data b%ag_data = prec%ag_data
b%outer_sweeps = prec%outer_sweeps b%outer_sweeps = prec%outer_sweeps
call move_alloc(prec%precv,b%precv) call move_alloc(prec%precv,b%precv)
! Fix the pointers except on level 1. ! Fix the pointers except on level 1.
do i=2, size(b%precv) do i=2, size(b%precv)
@ -878,7 +878,7 @@ contains
b%precv(i)%map%p_desc_U => b%precv(i-1)%base_desc b%precv(i)%map%p_desc_U => b%precv(i-1)%base_desc
b%precv(i)%map%p_desc_V => b%precv(i)%base_desc b%precv(i)%map%p_desc_V => b%precv(i)%base_desc
end do end do
else else
write(0,*) 'Warning: PREC%move_alloc onto different type?' write(0,*) 'Warning: PREC%move_alloc onto different type?'
info = psb_err_internal_error_ info = psb_err_internal_error_
@ -888,7 +888,7 @@ contains
subroutine amg_z_allocate_wrk(prec,info,vmold,desc) subroutine amg_z_allocate_wrk(prec,info,vmold,desc)
use psb_base_mod use psb_base_mod
implicit none implicit none
! Arguments ! Arguments
class(amg_zprec_type), intent(inout) :: prec class(amg_zprec_type), intent(inout) :: prec
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -896,37 +896,37 @@ contains
! !
! In MLD the DESC optional argument is ignored, since ! In MLD the DESC optional argument is ignored, since
! the necessary info is contained in the various entries of the ! the necessary info is contained in the various entries of the
! PRECV component. ! PRECV component.
type(psb_desc_type), intent(in), optional :: desc type(psb_desc_type), intent(in), optional :: desc
! Local variables ! Local variables
integer(psb_ipk_) :: me,err_act,i,j,level,nlev, nc2l integer(psb_ipk_) :: me,err_act,i,j,level,nlev, nc2l
character(len=20) :: name character(len=20) :: name
info=psb_success_ info=psb_success_
name = 'amg_z_allocate_wrk' name = 'amg_z_allocate_wrk'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999 info = psb_err_internal_error_; goto 9999
end if end if
nlev = size(prec%precv) nlev = size(prec%precv)
level = 1 level = 1
do level = 1, nlev do level = 1, nlev
call prec%precv(level)%allocate_wrk(info,vmold=vmold) call prec%precv(level)%allocate_wrk(info,vmold=vmold)
if (psb_errstatus_fatal()) then if (psb_errstatus_fatal()) then
nc2l = prec%precv(level)%base_desc%get_local_cols() nc2l = prec%precv(level)%base_desc%get_local_cols()
info=psb_err_alloc_request_ info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l/), a_err='complex(psb_dpk_)') call psb_errpush(info,name,i_err=(/2*nc2l/), a_err='complex(psb_dpk_)')
goto 9999 goto 9999
end if end if
end do end do
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
9999 call psb_error_handler(err_act) 9999 call psb_error_handler(err_act)
return return
end subroutine amg_z_allocate_wrk end subroutine amg_z_allocate_wrk
subroutine amg_z_free_wrk(prec,info) subroutine amg_z_free_wrk(prec,info)
@ -948,13 +948,13 @@ contains
info = psb_err_internal_error_; goto 9999 info = psb_err_internal_error_; goto 9999
end if end if
if (allocated(prec%precv)) then if (allocated(prec%precv)) then
nlev = size(prec%precv) nlev = size(prec%precv)
do level = 1, nlev do level = 1, nlev
call prec%precv(level)%free_wrk(info) call prec%precv(level)%free_wrk(info)
end do end do
end if end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -966,7 +966,7 @@ contains
function amg_z_is_allocated_wrk(prec) result(res) function amg_z_is_allocated_wrk(prec) result(res)
use psb_base_mod use psb_base_mod
implicit none implicit none
! Arguments ! Arguments
class(amg_zprec_type), intent(in) :: prec class(amg_zprec_type), intent(in) :: prec
logical :: res logical :: res
@ -974,7 +974,7 @@ contains
res = .false. res = .false.
if (.not.allocated(prec%precv)) return if (.not.allocated(prec%precv)) return
res = allocated(prec%precv(1)%wrk) res = allocated(prec%precv(1)%wrk)
end function amg_z_is_allocated_wrk end function amg_z_is_allocated_wrk
end module amg_z_prec_type end module amg_z_prec_type

Write Preview
Loading…
Cancel
Save