|
|
|
@ -185,10 +185,6 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
integer :: debug_level, debug_unit
|
|
|
|
|
integer :: ismth, nlev, ilev, icm
|
|
|
|
|
character(len=20) :: name
|
|
|
|
|
type psb_mlprec_wrk_type
|
|
|
|
|
real(kind(1.d0)), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
|
|
|
|
|
end type psb_mlprec_wrk_type
|
|
|
|
|
type(psb_mlprec_wrk_type), allocatable :: mlprec_wrk(:)
|
|
|
|
|
|
|
|
|
|
name='mld_dmlprec_aply'
|
|
|
|
|
info = 0
|
|
|
|
@ -203,14 +199,6 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
& write(debug_unit,*) me,' ',trim(name),&
|
|
|
|
|
& ' Entry ', size(baseprecv)
|
|
|
|
|
|
|
|
|
|
nlev = size(baseprecv)
|
|
|
|
|
allocate(mlprec_wrk(nlev),stat=info)
|
|
|
|
|
if (info /= 0) then
|
|
|
|
|
call psb_errpush(4010,name,a_err='Allocate')
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
select case(baseprecv(2)%iprcparm(mld_ml_type_))
|
|
|
|
|
|
|
|
|
|
case(mld_no_ml_)
|
|
|
|
@ -223,6 +211,103 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
|
|
|
|
|
case(mld_add_ml_)
|
|
|
|
|
|
|
|
|
|
call add_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
|
|
|
|
|
case(mld_mult_ml_)
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
! Multiplicative multilevel (multiplicative among the levels, additive inside
|
|
|
|
|
! each level)
|
|
|
|
|
!
|
|
|
|
|
! Pre/post-smoothing versions
|
|
|
|
|
!
|
|
|
|
|
|
|
|
|
|
select case(baseprecv(2)%iprcparm(mld_smooth_pos_))
|
|
|
|
|
|
|
|
|
|
case(mld_post_smooth_)
|
|
|
|
|
|
|
|
|
|
call mlt_post_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
|
|
|
|
|
case(mld_pre_smooth_)
|
|
|
|
|
|
|
|
|
|
call mlt_pre_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
|
|
|
|
|
case(mld_twoside_smooth_)
|
|
|
|
|
|
|
|
|
|
call mlt_twoside_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
|
|
|
|
|
case default
|
|
|
|
|
info = 4013
|
|
|
|
|
call psb_errpush(info,name,a_err='invalid smooth_pos',&
|
|
|
|
|
& i_Err=(/baseprecv(2)%iprcparm(mld_smooth_pos_),0,0,0,0/))
|
|
|
|
|
goto 9999
|
|
|
|
|
|
|
|
|
|
end select
|
|
|
|
|
|
|
|
|
|
case default
|
|
|
|
|
info = 4013
|
|
|
|
|
call psb_errpush(info,name,a_err='invalid mltype',&
|
|
|
|
|
& i_Err=(/baseprecv(2)%iprcparm(mld_ml_type_),0,0,0,0/))
|
|
|
|
|
goto 9999
|
|
|
|
|
|
|
|
|
|
end select
|
|
|
|
|
|
|
|
|
|
call psb_erractionrestore(err_act)
|
|
|
|
|
return
|
|
|
|
|
|
|
|
|
|
9999 continue
|
|
|
|
|
call psb_erractionrestore(err_act)
|
|
|
|
|
if (err_act.eq.psb_act_abort_) then
|
|
|
|
|
call psb_error()
|
|
|
|
|
return
|
|
|
|
|
end if
|
|
|
|
|
return
|
|
|
|
|
|
|
|
|
|
contains
|
|
|
|
|
|
|
|
|
|
subroutine add_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
implicit none
|
|
|
|
|
! Arguments
|
|
|
|
|
type(psb_desc_type),intent(in) :: desc_data
|
|
|
|
|
type(mld_dbaseprc_type), intent(in) :: baseprecv(:)
|
|
|
|
|
real(kind(0.d0)),intent(in) :: alpha,beta
|
|
|
|
|
real(kind(0.d0)),intent(in) :: x(:)
|
|
|
|
|
real(kind(0.d0)),intent(inout) :: y(:)
|
|
|
|
|
character :: trans
|
|
|
|
|
real(kind(0.d0)),target :: work(:)
|
|
|
|
|
integer, intent(out) :: info
|
|
|
|
|
|
|
|
|
|
! Local variables
|
|
|
|
|
integer :: n_row,n_col
|
|
|
|
|
integer :: ictxt,np,me,i, nr2l,nc2l,err_act
|
|
|
|
|
integer :: debug_level, debug_unit
|
|
|
|
|
integer :: ismth, nlev, ilev, icm
|
|
|
|
|
character(len=20) :: name
|
|
|
|
|
type psb_mlprec_wrk_type
|
|
|
|
|
real(kind(1.d0)), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
|
|
|
|
|
end type psb_mlprec_wrk_type
|
|
|
|
|
type(psb_mlprec_wrk_type), allocatable :: mlprec_wrk(:)
|
|
|
|
|
|
|
|
|
|
name='add_ml_aply'
|
|
|
|
|
info = 0
|
|
|
|
|
call psb_erractionsave(err_act)
|
|
|
|
|
debug_unit = psb_get_debug_unit()
|
|
|
|
|
debug_level = psb_get_debug_level()
|
|
|
|
|
|
|
|
|
|
ictxt = psb_cd_get_context(desc_data)
|
|
|
|
|
call psb_info(ictxt, me, np)
|
|
|
|
|
|
|
|
|
|
if (debug_level >= psb_debug_inner_) &
|
|
|
|
|
& write(debug_unit,*) me,' ',trim(name),&
|
|
|
|
|
& ' Entry ', size(baseprecv)
|
|
|
|
|
|
|
|
|
|
nlev = size(baseprecv)
|
|
|
|
|
allocate(mlprec_wrk(nlev),stat=info)
|
|
|
|
|
if (info /= 0) then
|
|
|
|
|
call psb_errpush(4010,name,a_err='Allocate')
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
! Additive multilevel
|
|
|
|
|
!
|
|
|
|
@ -398,48 +483,105 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
call psb_errpush(4001,name,a_err='Error on final update')
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
deallocate(mlprec_wrk,stat=info)
|
|
|
|
|
if (info /= 0) then
|
|
|
|
|
call psb_errpush(4000,name)
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
case(mld_mult_ml_)
|
|
|
|
|
call psb_erractionrestore(err_act)
|
|
|
|
|
return
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
! Multiplicative multilevel (multiplicative among the levels, additive inside
|
|
|
|
|
! each level)
|
|
|
|
|
!
|
|
|
|
|
! Pre/post-smoothing versions
|
|
|
|
|
!
|
|
|
|
|
9999 continue
|
|
|
|
|
call psb_erractionrestore(err_act)
|
|
|
|
|
if (err_act.eq.psb_act_abort_) then
|
|
|
|
|
call psb_error()
|
|
|
|
|
return
|
|
|
|
|
end if
|
|
|
|
|
return
|
|
|
|
|
|
|
|
|
|
select case(baseprecv(2)%iprcparm(mld_smooth_pos_))
|
|
|
|
|
end subroutine add_ml_aply
|
|
|
|
|
|
|
|
|
|
case(mld_post_smooth_)
|
|
|
|
|
|
|
|
|
|
subroutine mlt_pre_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
implicit none
|
|
|
|
|
! Arguments
|
|
|
|
|
type(psb_desc_type),intent(in) :: desc_data
|
|
|
|
|
type(mld_dbaseprc_type), intent(in) :: baseprecv(:)
|
|
|
|
|
real(kind(0.d0)),intent(in) :: alpha,beta
|
|
|
|
|
real(kind(0.d0)),intent(in) :: x(:)
|
|
|
|
|
real(kind(0.d0)),intent(inout) :: y(:)
|
|
|
|
|
character :: trans
|
|
|
|
|
real(kind(0.d0)),target :: work(:)
|
|
|
|
|
integer, intent(out) :: info
|
|
|
|
|
|
|
|
|
|
! Local variables
|
|
|
|
|
integer :: n_row,n_col
|
|
|
|
|
integer :: ictxt,np,me,i, nr2l,nc2l,err_act
|
|
|
|
|
integer :: debug_level, debug_unit
|
|
|
|
|
integer :: ismth, nlev, ilev, icm
|
|
|
|
|
character(len=20) :: name
|
|
|
|
|
type psb_mlprec_wrk_type
|
|
|
|
|
real(kind(1.d0)), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
|
|
|
|
|
end type psb_mlprec_wrk_type
|
|
|
|
|
type(psb_mlprec_wrk_type), allocatable :: mlprec_wrk(:)
|
|
|
|
|
|
|
|
|
|
name='mlt_pre_ml_aply'
|
|
|
|
|
info = 0
|
|
|
|
|
call psb_erractionsave(err_act)
|
|
|
|
|
debug_unit = psb_get_debug_unit()
|
|
|
|
|
debug_level = psb_get_debug_level()
|
|
|
|
|
|
|
|
|
|
ictxt = psb_cd_get_context(desc_data)
|
|
|
|
|
call psb_info(ictxt, me, np)
|
|
|
|
|
|
|
|
|
|
if (debug_level >= psb_debug_inner_) &
|
|
|
|
|
& write(debug_unit,*) me,' ',trim(name),&
|
|
|
|
|
& ' Entry ', size(baseprecv)
|
|
|
|
|
|
|
|
|
|
nlev = size(baseprecv)
|
|
|
|
|
allocate(mlprec_wrk(nlev),stat=info)
|
|
|
|
|
if (info /= 0) then
|
|
|
|
|
call psb_errpush(4010,name,a_err='Allocate')
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
! Post-smoothing
|
|
|
|
|
! Pre-smoothing
|
|
|
|
|
!
|
|
|
|
|
! 1. X(1) = Xext
|
|
|
|
|
!
|
|
|
|
|
! 2. DO ilev=2, nlev
|
|
|
|
|
! 2. ! Apply the base preconditioner at the finest level.
|
|
|
|
|
! Y(1) = (K(1)^(-1))*X(1)
|
|
|
|
|
!
|
|
|
|
|
! ! Transfer X(ilev-1) to the next coarser level.
|
|
|
|
|
! X(ilev) = AV(ilev; sm_pr_t_)*X(ilev-1)
|
|
|
|
|
! 3. ! Compute the residual at the finest level.
|
|
|
|
|
! TX(1) = X(1) - A(1)*Y(1)
|
|
|
|
|
!
|
|
|
|
|
! ENDDO
|
|
|
|
|
! 4. DO ilev=2, nlev
|
|
|
|
|
!
|
|
|
|
|
! 3.! Apply the preconditioner at the coarsest level.
|
|
|
|
|
! Y(nlev) = (K(nlev)^(-1))*X(nlev)
|
|
|
|
|
! ! Transfer the residual to the current (coarser) level.
|
|
|
|
|
! X(ilev) = AV(ilev; sm_pr_t_)*TX(ilev-1)
|
|
|
|
|
!
|
|
|
|
|
! 4. DO ilev=nlev-1,1,-1
|
|
|
|
|
! ! Apply the base preconditioner at the current level.
|
|
|
|
|
! ! The sum over the subdomains is carried out in the
|
|
|
|
|
! ! application of K(ilev).
|
|
|
|
|
! Y(ilev) = (K(ilev)^(-1))*X(ilev)
|
|
|
|
|
!
|
|
|
|
|
! ! Transfer Y(ilev+1) to the next finer level.
|
|
|
|
|
! Y(ilev) = AV(ilev+1; sm_pr_)*Y(ilev+1)
|
|
|
|
|
! ! Compute the residual at the current level (except at
|
|
|
|
|
! ! the coarsest level).
|
|
|
|
|
! IF (ilev < nlev)
|
|
|
|
|
! TX(ilev) = (X(ilev)-A(ilev)*Y(ilev))
|
|
|
|
|
!
|
|
|
|
|
! ! Compute the residual at the current level and apply to it the
|
|
|
|
|
! ! base preconditioner. The sum over the subdomains is carried out
|
|
|
|
|
! ! in the application of K(ilev).
|
|
|
|
|
! Y(ilev) = Y(ilev) + (K(ilev)^(-1))*(X(ilev)-A(ilev)*Y(ilev))
|
|
|
|
|
! ENDDO
|
|
|
|
|
!
|
|
|
|
|
! 5. DO ilev=nlev-1,1,-1
|
|
|
|
|
!
|
|
|
|
|
! ! Transfer Y(ilev+1) to the next finer level
|
|
|
|
|
! Y(ilev) = Y(ilev) + AV(ilev+1; sm_pr_)*Y(ilev+1)
|
|
|
|
|
!
|
|
|
|
|
! ENDDO
|
|
|
|
|
!
|
|
|
|
|
! 5. Yext = beta*Yext + alpha*Y(1)
|
|
|
|
|
! 6. Yext = beta*Yext + alpha*Y(1)
|
|
|
|
|
!
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
@ -447,30 +589,55 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
!
|
|
|
|
|
! Copy the input vector X
|
|
|
|
|
!
|
|
|
|
|
if (debug_level >= psb_debug_inner_) &
|
|
|
|
|
& write(debug_unit,*) me,' ',trim(name),&
|
|
|
|
|
& ' desc_data status',allocated(desc_data%matrix_data)
|
|
|
|
|
|
|
|
|
|
n_col = psb_cd_get_local_cols(desc_data)
|
|
|
|
|
nc2l = psb_cd_get_local_cols(baseprecv(1)%desc_data)
|
|
|
|
|
|
|
|
|
|
allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), &
|
|
|
|
|
& mlprec_wrk(1)%tx(nc2l), stat=info)
|
|
|
|
|
mlprec_wrk(1)%x2l(:) = dzero
|
|
|
|
|
mlprec_wrk(1)%y2l(:) = dzero
|
|
|
|
|
mlprec_wrk(1)%tx(:) = dzero
|
|
|
|
|
if (info /= 0) then
|
|
|
|
|
info=4025
|
|
|
|
|
call psb_errpush(info,name,i_err=(/4*nc2l,0,0,0,0/),&
|
|
|
|
|
& a_err='real(kind(1.d0))')
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%tx,&
|
|
|
|
|
& baseprecv(1)%base_desc,info)
|
|
|
|
|
call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%x2l,&
|
|
|
|
|
& baseprecv(1)%base_desc,info)
|
|
|
|
|
mlprec_wrk(1)%y2l(:) = dzero
|
|
|
|
|
mlprec_wrk(1)%x2l(:) = x
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
! STEP 2
|
|
|
|
|
!
|
|
|
|
|
! Apply the base preconditioner at the finest level
|
|
|
|
|
!
|
|
|
|
|
call mld_baseprec_aply(done,baseprecv(1),mlprec_wrk(1)%x2l,&
|
|
|
|
|
& dzero,mlprec_wrk(1)%y2l,&
|
|
|
|
|
& baseprecv(1)%base_desc,&
|
|
|
|
|
& trans,work,info)
|
|
|
|
|
if (info /=0) then
|
|
|
|
|
call psb_errpush(4010,name,a_err=' baseprec_aply')
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
! STEP 3
|
|
|
|
|
!
|
|
|
|
|
! Compute the residual at the finest level
|
|
|
|
|
!
|
|
|
|
|
mlprec_wrk(1)%tx = mlprec_wrk(1)%x2l
|
|
|
|
|
|
|
|
|
|
call psb_spmm(-done,baseprecv(1)%base_a,mlprec_wrk(1)%y2l,&
|
|
|
|
|
& done,mlprec_wrk(1)%tx,baseprecv(1)%base_desc,info,work=work)
|
|
|
|
|
if (info /=0) then
|
|
|
|
|
call psb_errpush(4001,name,a_err=' fine level residual')
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
! STEP 4
|
|
|
|
|
!
|
|
|
|
|
! For each level but the finest one ...
|
|
|
|
|
!
|
|
|
|
|
do ilev=2, nlev
|
|
|
|
|
do ilev = 2, nlev
|
|
|
|
|
|
|
|
|
|
n_row = psb_cd_get_local_rows(baseprecv(ilev-1)%base_desc)
|
|
|
|
|
n_col = psb_cd_get_local_cols(baseprecv(ilev-1)%desc_data)
|
|
|
|
@ -479,15 +646,8 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
ismth = baseprecv(ilev)%iprcparm(mld_smooth_kind_)
|
|
|
|
|
icm = baseprecv(ilev)%iprcparm(mld_coarse_mat_)
|
|
|
|
|
|
|
|
|
|
if (debug_level >= psb_debug_inner_) &
|
|
|
|
|
& write(debug_unit,*) me,' ',trim(name), &
|
|
|
|
|
& ' starting up sweep ',&
|
|
|
|
|
& ilev,allocated(baseprecv(ilev)%iprcparm),n_row,n_col,&
|
|
|
|
|
& nc2l, nr2l,ismth
|
|
|
|
|
|
|
|
|
|
allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%y2l(nc2l),&
|
|
|
|
|
& mlprec_wrk(ilev)%x2l(nc2l), stat=info)
|
|
|
|
|
|
|
|
|
|
if (info /= 0) then
|
|
|
|
|
info=4025
|
|
|
|
|
call psb_errpush(info,name,i_err=(/4*nc2l,0,0,0,0/),&
|
|
|
|
@ -503,32 +663,29 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
!
|
|
|
|
|
! Apply the smoothed prolongator transpose
|
|
|
|
|
!
|
|
|
|
|
if (debug_level >= psb_debug_inner_) &
|
|
|
|
|
& write(debug_unit,*) me,' ',trim(name), ' up sweep ', ilev
|
|
|
|
|
|
|
|
|
|
call psb_halo(mlprec_wrk(ilev-1)%x2l,&
|
|
|
|
|
& baseprecv(ilev-1)%base_desc,info,work=work)
|
|
|
|
|
call psb_halo(mlprec_wrk(ilev-1)%tx,baseprecv(ilev-1)%base_desc,&
|
|
|
|
|
& info,work=work)
|
|
|
|
|
if (info == 0) call psb_csmm(done,baseprecv(ilev)%av(mld_sm_pr_t_),&
|
|
|
|
|
& mlprec_wrk(ilev-1)%x2l,dzero,mlprec_wrk(ilev)%x2l,info)
|
|
|
|
|
& mlprec_wrk(ilev-1)%tx,dzero,mlprec_wrk(ilev)%x2l,info)
|
|
|
|
|
else
|
|
|
|
|
!
|
|
|
|
|
! Apply the raw aggregation map transpose (take a shortcut)
|
|
|
|
|
!
|
|
|
|
|
mlprec_wrk(ilev)%x2l = dzero
|
|
|
|
|
do i=1,n_row
|
|
|
|
|
mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) = &
|
|
|
|
|
& mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) + &
|
|
|
|
|
& mlprec_wrk(ilev-1)%x2l(i)
|
|
|
|
|
& mlprec_wrk(ilev-1)%tx(i)
|
|
|
|
|
end do
|
|
|
|
|
|
|
|
|
|
end if
|
|
|
|
|
if (info /=0) then
|
|
|
|
|
call psb_errpush(4001,name,a_err='Error during restriction')
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
if (icm == mld_repl_mat_) Then
|
|
|
|
|
if (icm ==mld_repl_mat_) then
|
|
|
|
|
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nr2l))
|
|
|
|
|
else if (icm /= mld_distr_mat_) Then
|
|
|
|
|
else if (icm /= mld_distr_mat_) then
|
|
|
|
|
info = 4013
|
|
|
|
|
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_',&
|
|
|
|
|
& i_Err=(/icm,0,0,0,0/))
|
|
|
|
@ -536,46 +693,32 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
endif
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
! update x2l
|
|
|
|
|
! Apply the base preconditioner
|
|
|
|
|
!
|
|
|
|
|
call psb_geaxpby(done,mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%tx,&
|
|
|
|
|
& baseprecv(ilev)%base_desc,info)
|
|
|
|
|
if (info /= 0) then
|
|
|
|
|
call psb_errpush(4001,name,a_err='Error in update')
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
if (debug_level >= psb_debug_inner_) &
|
|
|
|
|
& write(debug_unit,*) me,' ',trim(name),&
|
|
|
|
|
& ' done up sweep ', ilev
|
|
|
|
|
|
|
|
|
|
enddo
|
|
|
|
|
call mld_baseprec_aply(done,baseprecv(ilev),mlprec_wrk(ilev)%x2l,&
|
|
|
|
|
& dzero,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%desc_data, 'N',work,info)
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
! STEP 3
|
|
|
|
|
!
|
|
|
|
|
! Apply the base preconditioner at the coarsest level
|
|
|
|
|
! Compute the residual (at all levels but the coarsest one)
|
|
|
|
|
!
|
|
|
|
|
call mld_baseprec_aply(done,baseprecv(nlev),mlprec_wrk(nlev)%x2l, &
|
|
|
|
|
& dzero, mlprec_wrk(nlev)%y2l,baseprecv(nlev)%desc_data,'N',work,info)
|
|
|
|
|
if (ilev < nlev) then
|
|
|
|
|
mlprec_wrk(ilev)%tx = mlprec_wrk(ilev)%x2l
|
|
|
|
|
if (info == 0) call psb_spmm(-done,baseprecv(ilev)%base_a,&
|
|
|
|
|
& mlprec_wrk(ilev)%y2l,done,mlprec_wrk(ilev)%tx,&
|
|
|
|
|
& baseprecv(ilev)%base_desc,info,work=work)
|
|
|
|
|
endif
|
|
|
|
|
if (info /=0) then
|
|
|
|
|
call psb_errpush(4010,name,a_err='baseprec_aply')
|
|
|
|
|
call psb_errpush(4001,name,a_err='Error on up sweep residual')
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
if (debug_level >= psb_debug_inner_) write(debug_unit,*) &
|
|
|
|
|
& me,' ',trim(name), ' done baseprec_aply ', nlev
|
|
|
|
|
enddo
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
! STEP 4
|
|
|
|
|
! STEP 5
|
|
|
|
|
!
|
|
|
|
|
! For each level but the coarsest one ...
|
|
|
|
|
!
|
|
|
|
|
do ilev=nlev-1, 1, -1
|
|
|
|
|
|
|
|
|
|
if (debug_level >= psb_debug_inner_) &
|
|
|
|
|
& write(debug_unit,*) me,' ',trim(name),&
|
|
|
|
|
& ' starting down sweep',ilev
|
|
|
|
|
do ilev = nlev-1, 1, -1
|
|
|
|
|
|
|
|
|
|
ismth = baseprecv(ilev+1)%iprcparm(mld_smooth_kind_)
|
|
|
|
|
n_row = psb_cd_get_local_rows(baseprecv(ilev)%base_desc)
|
|
|
|
@ -585,15 +728,14 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
! Apply the smoothed prolongator
|
|
|
|
|
!
|
|
|
|
|
if (ismth == mld_smooth_prol_) &
|
|
|
|
|
& call psb_halo(mlprec_wrk(ilev+1)%y2l,baseprecv(ilev+1)%desc_data,&
|
|
|
|
|
& info,work=work)
|
|
|
|
|
& call psb_halo(mlprec_wrk(ilev+1)%y2l,&
|
|
|
|
|
& baseprecv(ilev+1)%desc_data,info,work=work)
|
|
|
|
|
if (info == 0) call psb_csmm(done,baseprecv(ilev+1)%av(mld_sm_pr_),&
|
|
|
|
|
& mlprec_wrk(ilev+1)%y2l, dzero,mlprec_wrk(ilev)%y2l,info)
|
|
|
|
|
& mlprec_wrk(ilev+1)%y2l,done,mlprec_wrk(ilev)%y2l,info)
|
|
|
|
|
else
|
|
|
|
|
!
|
|
|
|
|
! Apply the raw aggregation map (take a shortcut)
|
|
|
|
|
!
|
|
|
|
|
mlprec_wrk(ilev)%y2l(:) = dzero
|
|
|
|
|
do i=1, n_row
|
|
|
|
|
mlprec_wrk(ilev)%y2l(i) = mlprec_wrk(ilev)%y2l(i) + &
|
|
|
|
|
& mlprec_wrk(ilev+1)%y2l(baseprecv(ilev+1)%mlia(i))
|
|
|
|
@ -603,78 +745,117 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
call psb_errpush(4001,name,a_err='Error during prolongation')
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
enddo
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
! Compute the residual
|
|
|
|
|
!
|
|
|
|
|
call psb_spmm(-done,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
|
|
|
|
|
& done,mlprec_wrk(ilev)%tx,baseprecv(ilev)%base_desc,info,work=work)
|
|
|
|
|
|
|
|
|
|
! STEP 6
|
|
|
|
|
!
|
|
|
|
|
! Apply the base preconditioner
|
|
|
|
|
! Compute the output vector Y
|
|
|
|
|
!
|
|
|
|
|
if (info == 0) call mld_baseprec_aply(done,baseprecv(ilev),mlprec_wrk(ilev)%tx,&
|
|
|
|
|
& done,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%base_desc, trans, work,info)
|
|
|
|
|
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,&
|
|
|
|
|
& baseprecv(1)%base_desc,info)
|
|
|
|
|
if (info /=0) then
|
|
|
|
|
call psb_errpush(4001,name,a_err=' spmm/baseprec_aply')
|
|
|
|
|
call psb_errpush(4001,name,a_err='Error on final update')
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
if (debug_level >= psb_debug_inner_) &
|
|
|
|
|
& write(debug_unit,*) me,' ',trim(name),&
|
|
|
|
|
& ' done down sweep',ilev
|
|
|
|
|
enddo
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
! STEP 5
|
|
|
|
|
!
|
|
|
|
|
! Compute the output vector Y
|
|
|
|
|
!
|
|
|
|
|
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,baseprecv(1)%base_desc,info)
|
|
|
|
|
|
|
|
|
|
if (info /=0) then
|
|
|
|
|
call psb_errpush(4001,name,a_err=' Final update')
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
deallocate(mlprec_wrk,stat=info)
|
|
|
|
|
if (info /= 0) then
|
|
|
|
|
call psb_errpush(4000,name)
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
case(mld_pre_smooth_)
|
|
|
|
|
call psb_erractionrestore(err_act)
|
|
|
|
|
return
|
|
|
|
|
|
|
|
|
|
9999 continue
|
|
|
|
|
call psb_erractionrestore(err_act)
|
|
|
|
|
if (err_act.eq.psb_act_abort_) then
|
|
|
|
|
call psb_error()
|
|
|
|
|
return
|
|
|
|
|
end if
|
|
|
|
|
return
|
|
|
|
|
|
|
|
|
|
end subroutine mlt_pre_ml_aply
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
subroutine mlt_post_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
implicit none
|
|
|
|
|
! Arguments
|
|
|
|
|
type(psb_desc_type),intent(in) :: desc_data
|
|
|
|
|
type(mld_dbaseprc_type), intent(in) :: baseprecv(:)
|
|
|
|
|
real(kind(0.d0)),intent(in) :: alpha,beta
|
|
|
|
|
real(kind(0.d0)),intent(in) :: x(:)
|
|
|
|
|
real(kind(0.d0)),intent(inout) :: y(:)
|
|
|
|
|
character :: trans
|
|
|
|
|
real(kind(0.d0)),target :: work(:)
|
|
|
|
|
integer, intent(out) :: info
|
|
|
|
|
|
|
|
|
|
! Local variables
|
|
|
|
|
integer :: n_row,n_col
|
|
|
|
|
integer :: ictxt,np,me,i, nr2l,nc2l,err_act
|
|
|
|
|
integer :: debug_level, debug_unit
|
|
|
|
|
integer :: ismth, nlev, ilev, icm
|
|
|
|
|
character(len=20) :: name
|
|
|
|
|
type psb_mlprec_wrk_type
|
|
|
|
|
real(kind(1.d0)), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
|
|
|
|
|
end type psb_mlprec_wrk_type
|
|
|
|
|
type(psb_mlprec_wrk_type), allocatable :: mlprec_wrk(:)
|
|
|
|
|
|
|
|
|
|
name='mlt_post_ml_aply'
|
|
|
|
|
info = 0
|
|
|
|
|
call psb_erractionsave(err_act)
|
|
|
|
|
debug_unit = psb_get_debug_unit()
|
|
|
|
|
debug_level = psb_get_debug_level()
|
|
|
|
|
|
|
|
|
|
ictxt = psb_cd_get_context(desc_data)
|
|
|
|
|
call psb_info(ictxt, me, np)
|
|
|
|
|
|
|
|
|
|
if (debug_level >= psb_debug_inner_) &
|
|
|
|
|
& write(debug_unit,*) me,' ',trim(name),&
|
|
|
|
|
& ' Entry ', size(baseprecv)
|
|
|
|
|
|
|
|
|
|
nlev = size(baseprecv)
|
|
|
|
|
allocate(mlprec_wrk(nlev),stat=info)
|
|
|
|
|
if (info /= 0) then
|
|
|
|
|
call psb_errpush(4010,name,a_err='Allocate')
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
! Pre-smoothing
|
|
|
|
|
! Post-smoothing
|
|
|
|
|
!
|
|
|
|
|
! 1. X(1) = Xext
|
|
|
|
|
!
|
|
|
|
|
! 2. ! Apply the base preconditioner at the finest level.
|
|
|
|
|
! Y(1) = (K(1)^(-1))*X(1)
|
|
|
|
|
!
|
|
|
|
|
! 3. ! Compute the residual at the finest level.
|
|
|
|
|
! TX(1) = X(1) - A(1)*Y(1)
|
|
|
|
|
!
|
|
|
|
|
! 4. DO ilev=2, nlev
|
|
|
|
|
!
|
|
|
|
|
! ! Transfer the residual to the current (coarser) level.
|
|
|
|
|
! X(ilev) = AV(ilev; sm_pr_t_)*TX(ilev-1)
|
|
|
|
|
!
|
|
|
|
|
! ! Apply the base preconditioner at the current level.
|
|
|
|
|
! ! The sum over the subdomains is carried out in the
|
|
|
|
|
! ! application of K(ilev).
|
|
|
|
|
! Y(ilev) = (K(ilev)^(-1))*X(ilev)
|
|
|
|
|
! 2. DO ilev=2, nlev
|
|
|
|
|
!
|
|
|
|
|
! ! Compute the residual at the current level (except at
|
|
|
|
|
! ! the coarsest level).
|
|
|
|
|
! IF (ilev < nlev)
|
|
|
|
|
! TX(ilev) = (X(ilev)-A(ilev)*Y(ilev))
|
|
|
|
|
! ! Transfer X(ilev-1) to the next coarser level.
|
|
|
|
|
! X(ilev) = AV(ilev; sm_pr_t_)*X(ilev-1)
|
|
|
|
|
!
|
|
|
|
|
! ENDDO
|
|
|
|
|
!
|
|
|
|
|
! 5. DO ilev=nlev-1,1,-1
|
|
|
|
|
! 3.! Apply the preconditioner at the coarsest level.
|
|
|
|
|
! Y(nlev) = (K(nlev)^(-1))*X(nlev)
|
|
|
|
|
!
|
|
|
|
|
! ! Transfer Y(ilev+1) to the next finer level
|
|
|
|
|
! Y(ilev) = Y(ilev) + AV(ilev+1; sm_pr_)*Y(ilev+1)
|
|
|
|
|
! 4. DO ilev=nlev-1,1,-1
|
|
|
|
|
!
|
|
|
|
|
! ! Transfer Y(ilev+1) to the next finer level.
|
|
|
|
|
! Y(ilev) = AV(ilev+1; sm_pr_)*Y(ilev+1)
|
|
|
|
|
!
|
|
|
|
|
! ! Compute the residual at the current level and apply to it the
|
|
|
|
|
! ! base preconditioner. The sum over the subdomains is carried out
|
|
|
|
|
! ! in the application of K(ilev).
|
|
|
|
|
! Y(ilev) = Y(ilev) + (K(ilev)^(-1))*(X(ilev)-A(ilev)*Y(ilev))
|
|
|
|
|
!
|
|
|
|
|
! ENDDO
|
|
|
|
|
!
|
|
|
|
|
! 6. Yext = beta*Yext + alpha*Y(1)
|
|
|
|
|
! 5. Yext = beta*Yext + alpha*Y(1)
|
|
|
|
|
!
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
@ -682,55 +863,30 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
!
|
|
|
|
|
! Copy the input vector X
|
|
|
|
|
!
|
|
|
|
|
if (debug_level >= psb_debug_inner_) &
|
|
|
|
|
& write(debug_unit,*) me,' ',trim(name),&
|
|
|
|
|
& ' desc_data status',allocated(desc_data%matrix_data)
|
|
|
|
|
|
|
|
|
|
n_col = psb_cd_get_local_cols(desc_data)
|
|
|
|
|
nc2l = psb_cd_get_local_cols(baseprecv(1)%desc_data)
|
|
|
|
|
|
|
|
|
|
allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), &
|
|
|
|
|
& mlprec_wrk(1)%tx(nc2l), stat=info)
|
|
|
|
|
if (info /= 0) then
|
|
|
|
|
info=4025
|
|
|
|
|
call psb_errpush(info,name,i_err=(/4*nc2l,0,0,0,0/),&
|
|
|
|
|
& a_err='real(kind(1.d0))')
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
mlprec_wrk(1)%x2l(:) = dzero
|
|
|
|
|
mlprec_wrk(1)%y2l(:) = dzero
|
|
|
|
|
mlprec_wrk(1)%x2l(:) = x
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
! STEP 2
|
|
|
|
|
!
|
|
|
|
|
! Apply the base preconditioner at the finest level
|
|
|
|
|
!
|
|
|
|
|
call mld_baseprec_aply(done,baseprecv(1),mlprec_wrk(1)%x2l,&
|
|
|
|
|
& dzero,mlprec_wrk(1)%y2l,&
|
|
|
|
|
& baseprecv(1)%base_desc,&
|
|
|
|
|
& trans,work,info)
|
|
|
|
|
if (info /=0) then
|
|
|
|
|
call psb_errpush(4010,name,a_err=' baseprec_aply')
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
! STEP 3
|
|
|
|
|
!
|
|
|
|
|
! Compute the residual at the finest level
|
|
|
|
|
!
|
|
|
|
|
mlprec_wrk(1)%tx = mlprec_wrk(1)%x2l
|
|
|
|
|
mlprec_wrk(1)%tx(:) = dzero
|
|
|
|
|
|
|
|
|
|
call psb_spmm(-done,baseprecv(1)%base_a,mlprec_wrk(1)%y2l,&
|
|
|
|
|
& done,mlprec_wrk(1)%tx,baseprecv(1)%base_desc,info,work=work)
|
|
|
|
|
if (info /=0) then
|
|
|
|
|
call psb_errpush(4001,name,a_err=' fine level residual')
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%tx,&
|
|
|
|
|
& baseprecv(1)%base_desc,info)
|
|
|
|
|
call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%x2l,&
|
|
|
|
|
& baseprecv(1)%base_desc,info)
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
! STEP 4
|
|
|
|
|
! STEP 2
|
|
|
|
|
!
|
|
|
|
|
! For each level but the finest one ...
|
|
|
|
|
!
|
|
|
|
|
do ilev = 2, nlev
|
|
|
|
|
do ilev=2, nlev
|
|
|
|
|
|
|
|
|
|
n_row = psb_cd_get_local_rows(baseprecv(ilev-1)%base_desc)
|
|
|
|
|
n_col = psb_cd_get_local_cols(baseprecv(ilev-1)%desc_data)
|
|
|
|
@ -739,8 +895,15 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
ismth = baseprecv(ilev)%iprcparm(mld_smooth_kind_)
|
|
|
|
|
icm = baseprecv(ilev)%iprcparm(mld_coarse_mat_)
|
|
|
|
|
|
|
|
|
|
if (debug_level >= psb_debug_inner_) &
|
|
|
|
|
& write(debug_unit,*) me,' ',trim(name), &
|
|
|
|
|
& ' starting up sweep ',&
|
|
|
|
|
& ilev,allocated(baseprecv(ilev)%iprcparm),n_row,n_col,&
|
|
|
|
|
& nc2l, nr2l,ismth
|
|
|
|
|
|
|
|
|
|
allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%y2l(nc2l),&
|
|
|
|
|
& mlprec_wrk(ilev)%x2l(nc2l), stat=info)
|
|
|
|
|
|
|
|
|
|
if (info /= 0) then
|
|
|
|
|
info=4025
|
|
|
|
|
call psb_errpush(info,name,i_err=(/4*nc2l,0,0,0,0/),&
|
|
|
|
@ -756,29 +919,32 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
!
|
|
|
|
|
! Apply the smoothed prolongator transpose
|
|
|
|
|
!
|
|
|
|
|
call psb_halo(mlprec_wrk(ilev-1)%tx,baseprecv(ilev-1)%base_desc,&
|
|
|
|
|
& info,work=work)
|
|
|
|
|
if (debug_level >= psb_debug_inner_) &
|
|
|
|
|
& write(debug_unit,*) me,' ',trim(name), ' up sweep ', ilev
|
|
|
|
|
|
|
|
|
|
call psb_halo(mlprec_wrk(ilev-1)%x2l,&
|
|
|
|
|
& baseprecv(ilev-1)%base_desc,info,work=work)
|
|
|
|
|
if (info == 0) call psb_csmm(done,baseprecv(ilev)%av(mld_sm_pr_t_),&
|
|
|
|
|
& mlprec_wrk(ilev-1)%tx,dzero,mlprec_wrk(ilev)%x2l,info)
|
|
|
|
|
& mlprec_wrk(ilev-1)%x2l,dzero,mlprec_wrk(ilev)%x2l,info)
|
|
|
|
|
else
|
|
|
|
|
!
|
|
|
|
|
! Apply the raw aggregation map transpose (take a shortcut)
|
|
|
|
|
!
|
|
|
|
|
mlprec_wrk(ilev)%x2l = dzero
|
|
|
|
|
do i=1,n_row
|
|
|
|
|
mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) = &
|
|
|
|
|
& mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) + &
|
|
|
|
|
& mlprec_wrk(ilev-1)%tx(i)
|
|
|
|
|
& mlprec_wrk(ilev-1)%x2l(i)
|
|
|
|
|
end do
|
|
|
|
|
|
|
|
|
|
end if
|
|
|
|
|
if (info /=0) then
|
|
|
|
|
call psb_errpush(4001,name,a_err='Error during restriction')
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
if (icm ==mld_repl_mat_) then
|
|
|
|
|
if (icm == mld_repl_mat_) Then
|
|
|
|
|
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nr2l))
|
|
|
|
|
else if (icm /= mld_distr_mat_) then
|
|
|
|
|
else if (icm /= mld_distr_mat_) Then
|
|
|
|
|
info = 4013
|
|
|
|
|
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_',&
|
|
|
|
|
& i_Err=(/icm,0,0,0,0/))
|
|
|
|
@ -786,32 +952,46 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
endif
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
! Apply the base preconditioner
|
|
|
|
|
! update x2l
|
|
|
|
|
!
|
|
|
|
|
call mld_baseprec_aply(done,baseprecv(ilev),mlprec_wrk(ilev)%x2l,&
|
|
|
|
|
& dzero,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%desc_data, 'N',work,info)
|
|
|
|
|
call psb_geaxpby(done,mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%tx,&
|
|
|
|
|
& baseprecv(ilev)%base_desc,info)
|
|
|
|
|
if (info /= 0) then
|
|
|
|
|
call psb_errpush(4001,name,a_err='Error in update')
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
if (debug_level >= psb_debug_inner_) &
|
|
|
|
|
& write(debug_unit,*) me,' ',trim(name),&
|
|
|
|
|
& ' done up sweep ', ilev
|
|
|
|
|
|
|
|
|
|
enddo
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
! Compute the residual (at all levels but the coarsest one)
|
|
|
|
|
! STEP 3
|
|
|
|
|
!
|
|
|
|
|
if (ilev < nlev) then
|
|
|
|
|
mlprec_wrk(ilev)%tx = mlprec_wrk(ilev)%x2l
|
|
|
|
|
if (info == 0) call psb_spmm(-done,baseprecv(ilev)%base_a,&
|
|
|
|
|
& mlprec_wrk(ilev)%y2l,done,mlprec_wrk(ilev)%tx,&
|
|
|
|
|
& baseprecv(ilev)%base_desc,info,work=work)
|
|
|
|
|
endif
|
|
|
|
|
! Apply the base preconditioner at the coarsest level
|
|
|
|
|
!
|
|
|
|
|
call mld_baseprec_aply(done,baseprecv(nlev),mlprec_wrk(nlev)%x2l, &
|
|
|
|
|
& dzero, mlprec_wrk(nlev)%y2l,baseprecv(nlev)%desc_data,'N',work,info)
|
|
|
|
|
if (info /=0) then
|
|
|
|
|
call psb_errpush(4001,name,a_err='Error on up sweep residual')
|
|
|
|
|
call psb_errpush(4010,name,a_err='baseprec_aply')
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
enddo
|
|
|
|
|
|
|
|
|
|
if (debug_level >= psb_debug_inner_) write(debug_unit,*) &
|
|
|
|
|
& me,' ',trim(name), ' done baseprec_aply ', nlev
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
! STEP 5
|
|
|
|
|
! STEP 4
|
|
|
|
|
!
|
|
|
|
|
! For each level but the coarsest one ...
|
|
|
|
|
!
|
|
|
|
|
do ilev = nlev-1, 1, -1
|
|
|
|
|
do ilev=nlev-1, 1, -1
|
|
|
|
|
|
|
|
|
|
if (debug_level >= psb_debug_inner_) &
|
|
|
|
|
& write(debug_unit,*) me,' ',trim(name),&
|
|
|
|
|
& ' starting down sweep',ilev
|
|
|
|
|
|
|
|
|
|
ismth = baseprecv(ilev+1)%iprcparm(mld_smooth_kind_)
|
|
|
|
|
n_row = psb_cd_get_local_rows(baseprecv(ilev)%base_desc)
|
|
|
|
@ -821,14 +1001,15 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
! Apply the smoothed prolongator
|
|
|
|
|
!
|
|
|
|
|
if (ismth == mld_smooth_prol_) &
|
|
|
|
|
& call psb_halo(mlprec_wrk(ilev+1)%y2l,&
|
|
|
|
|
& baseprecv(ilev+1)%desc_data,info,work=work)
|
|
|
|
|
& call psb_halo(mlprec_wrk(ilev+1)%y2l,baseprecv(ilev+1)%desc_data,&
|
|
|
|
|
& info,work=work)
|
|
|
|
|
if (info == 0) call psb_csmm(done,baseprecv(ilev+1)%av(mld_sm_pr_),&
|
|
|
|
|
& mlprec_wrk(ilev+1)%y2l,done,mlprec_wrk(ilev)%y2l,info)
|
|
|
|
|
& mlprec_wrk(ilev+1)%y2l, dzero,mlprec_wrk(ilev)%y2l,info)
|
|
|
|
|
else
|
|
|
|
|
!
|
|
|
|
|
! Apply the raw aggregation map (take a shortcut)
|
|
|
|
|
!
|
|
|
|
|
mlprec_wrk(ilev)%y2l(:) = dzero
|
|
|
|
|
do i=1, n_row
|
|
|
|
|
mlprec_wrk(ilev)%y2l(i) = mlprec_wrk(ilev)%y2l(i) + &
|
|
|
|
|
& mlprec_wrk(ilev+1)%y2l(baseprecv(ilev+1)%mlia(i))
|
|
|
|
@ -838,23 +1019,105 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
call psb_errpush(4001,name,a_err='Error during prolongation')
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
! Compute the residual
|
|
|
|
|
!
|
|
|
|
|
call psb_spmm(-done,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
|
|
|
|
|
& done,mlprec_wrk(ilev)%tx,baseprecv(ilev)%base_desc,info,work=work)
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
! Apply the base preconditioner
|
|
|
|
|
!
|
|
|
|
|
if (info == 0) call mld_baseprec_aply(done,baseprecv(ilev),mlprec_wrk(ilev)%tx,&
|
|
|
|
|
& done,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%base_desc, trans, work,info)
|
|
|
|
|
if (info /=0) then
|
|
|
|
|
call psb_errpush(4001,name,a_err=' spmm/baseprec_aply')
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
if (debug_level >= psb_debug_inner_) &
|
|
|
|
|
& write(debug_unit,*) me,' ',trim(name),&
|
|
|
|
|
& ' done down sweep',ilev
|
|
|
|
|
enddo
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
! STEP 6
|
|
|
|
|
! STEP 5
|
|
|
|
|
!
|
|
|
|
|
! Compute the output vector Y
|
|
|
|
|
!
|
|
|
|
|
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,&
|
|
|
|
|
& baseprecv(1)%base_desc,info)
|
|
|
|
|
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,baseprecv(1)%base_desc,info)
|
|
|
|
|
|
|
|
|
|
if (info /=0) then
|
|
|
|
|
call psb_errpush(4001,name,a_err='Error on final update')
|
|
|
|
|
call psb_errpush(4001,name,a_err=' Final update')
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
case(mld_twoside_smooth_)
|
|
|
|
|
deallocate(mlprec_wrk,stat=info)
|
|
|
|
|
if (info /= 0) then
|
|
|
|
|
call psb_errpush(4000,name)
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
call psb_erractionrestore(err_act)
|
|
|
|
|
return
|
|
|
|
|
|
|
|
|
|
9999 continue
|
|
|
|
|
call psb_erractionrestore(err_act)
|
|
|
|
|
if (err_act.eq.psb_act_abort_) then
|
|
|
|
|
call psb_error()
|
|
|
|
|
return
|
|
|
|
|
end if
|
|
|
|
|
return
|
|
|
|
|
|
|
|
|
|
end subroutine mlt_post_ml_aply
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
subroutine mlt_twoside_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
implicit none
|
|
|
|
|
! Arguments
|
|
|
|
|
type(psb_desc_type),intent(in) :: desc_data
|
|
|
|
|
type(mld_dbaseprc_type), intent(in) :: baseprecv(:)
|
|
|
|
|
real(kind(0.d0)),intent(in) :: alpha,beta
|
|
|
|
|
real(kind(0.d0)),intent(in) :: x(:)
|
|
|
|
|
real(kind(0.d0)),intent(inout) :: y(:)
|
|
|
|
|
character :: trans
|
|
|
|
|
real(kind(0.d0)),target :: work(:)
|
|
|
|
|
integer, intent(out) :: info
|
|
|
|
|
|
|
|
|
|
! Local variables
|
|
|
|
|
integer :: n_row,n_col
|
|
|
|
|
integer :: ictxt,np,me,i, nr2l,nc2l,err_act
|
|
|
|
|
integer :: debug_level, debug_unit
|
|
|
|
|
integer :: ismth, nlev, ilev, icm
|
|
|
|
|
character(len=20) :: name
|
|
|
|
|
type psb_mlprec_wrk_type
|
|
|
|
|
real(kind(1.d0)), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
|
|
|
|
|
end type psb_mlprec_wrk_type
|
|
|
|
|
type(psb_mlprec_wrk_type), allocatable :: mlprec_wrk(:)
|
|
|
|
|
|
|
|
|
|
name='mlt_twoside_ml_aply'
|
|
|
|
|
info = 0
|
|
|
|
|
call psb_erractionsave(err_act)
|
|
|
|
|
debug_unit = psb_get_debug_unit()
|
|
|
|
|
debug_level = psb_get_debug_level()
|
|
|
|
|
|
|
|
|
|
ictxt = psb_cd_get_context(desc_data)
|
|
|
|
|
call psb_info(ictxt, me, np)
|
|
|
|
|
|
|
|
|
|
if (debug_level >= psb_debug_inner_) &
|
|
|
|
|
& write(debug_unit,*) me,' ',trim(name),&
|
|
|
|
|
& ' Entry ', size(baseprecv)
|
|
|
|
|
|
|
|
|
|
nlev = size(baseprecv)
|
|
|
|
|
allocate(mlprec_wrk(nlev),stat=info)
|
|
|
|
|
if (info /= 0) then
|
|
|
|
|
call psb_errpush(4010,name,a_err='Allocate')
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
|
! Pre- and post-smoothing (symmetrized)
|
|
|
|
@ -1093,21 +1356,7 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
goto 9999
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
case default
|
|
|
|
|
info = 4013
|
|
|
|
|
call psb_errpush(info,name,a_err='invalid smooth_pos',&
|
|
|
|
|
& i_Err=(/baseprecv(2)%iprcparm(mld_smooth_pos_),0,0,0,0/))
|
|
|
|
|
goto 9999
|
|
|
|
|
|
|
|
|
|
end select
|
|
|
|
|
|
|
|
|
|
case default
|
|
|
|
|
info = 4013
|
|
|
|
|
call psb_errpush(info,name,a_err='invalid mltype',&
|
|
|
|
|
& i_Err=(/baseprecv(2)%iprcparm(mld_ml_type_),0,0,0,0/))
|
|
|
|
|
goto 9999
|
|
|
|
|
|
|
|
|
|
end select
|
|
|
|
|
|
|
|
|
|
deallocate(mlprec_wrk,stat=info)
|
|
|
|
|
if (info /= 0) then
|
|
|
|
@ -1126,5 +1375,10 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|
end if
|
|
|
|
|
return
|
|
|
|
|
|
|
|
|
|
end subroutine mlt_twoside_ml_aply
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
end subroutine mld_dmlprec_aply
|
|
|
|
|
|
|
|
|
|