mirror of
https://github.com/sfilippone/amg4psblas.git
synced 2026-10-07 07:04:59 +00:00
Compare commits
5
Commits
| Author | SHA1 | Date | |
|---|---|---|---|
 Salvatore Filippone
|
3f9abee095 | ||
 Salvatore Filippone
|
d9c5cfb8e9 | ||
|
Salvatore Filippone
|
4ba1232533 | ||
|
Salvatore Filippone
|
a3f8d802ff | ||
 Cirdans-Home
|
a704873923 |
@@ -29,7 +29,7 @@ extern "C" {
|
||||
psb_i_t mld_c_dprecfree(mld_c_dprec *ph);
|
||||
psb_i_t mld_c_dprecbld_opt(psb_c_dspmat *ah, psb_c_descriptor *cdh,
|
||||
mld_c_dprec *ph, const char *afmt);
|
||||
|
||||
psb_i_t mld_c_ddescr(mld_c_dprec *ph);
|
||||
|
||||
psb_i_t mld_c_dkrylov(const char *method, psb_c_dspmat *ah, mld_c_dprec *ph,
|
||||
psb_c_dvector *bh, psb_c_dvector *xh,
|
||||
|
||||
@@ -30,6 +30,7 @@ extern "C" {
|
||||
psb_i_t mld_c_zprecbld_opt(psb_c_zspmat *ah, psb_c_descriptor *cdh,
|
||||
mld_c_zprec *ph, const char *afmt);
|
||||
|
||||
psb_i_t mld_c_zdescr(mld_c_zprec *ph);
|
||||
|
||||
psb_i_t mld_c_zkrylov(const char *method, psb_c_zspmat *ah, mld_c_zprec *ph,
|
||||
psb_c_zvector *bh, psb_c_zvector *xh,
|
||||
|
||||
@@ -377,4 +377,32 @@ contains
|
||||
return
|
||||
end function mld_c_dprecfree
|
||||
|
||||
function mld_c_ddescr(ph) bind(c) result(res)
|
||||
use psb_base_mod
|
||||
use mld_prec_mod
|
||||
implicit none
|
||||
|
||||
integer(psb_c_ipk_) :: res
|
||||
type(psb_c_object_type) :: ph
|
||||
integer :: info
|
||||
type(mld_dprec_type), pointer :: precp
|
||||
|
||||
res = -1
|
||||
info = -1
|
||||
if (c_associated(ph%item)) then
|
||||
call c_f_pointer(ph%item,precp)
|
||||
else
|
||||
return
|
||||
end if
|
||||
|
||||
|
||||
call precp%descr()
|
||||
call flush(output_unit)
|
||||
|
||||
info = 0
|
||||
res = MLDC_ERR_FILTER(info)
|
||||
MLDC_ERR_HANDLE(res)
|
||||
return
|
||||
end function mld_c_ddescr
|
||||
|
||||
end module mld_dprec_cbind_mod
|
||||
|
||||
@@ -377,4 +377,32 @@ contains
|
||||
return
|
||||
end function mld_c_zprecfree
|
||||
|
||||
function mld_c_zdescr(ph) bind(c) result(res)
|
||||
use psb_base_mod
|
||||
use mld_prec_mod
|
||||
implicit none
|
||||
|
||||
integer(psb_c_ipk_) :: res
|
||||
type(psb_c_object_type) :: ph
|
||||
integer :: info
|
||||
type(mld_zprec_type), pointer :: precp
|
||||
|
||||
res = -1
|
||||
info = -1
|
||||
if (c_associated(ph%item)) then
|
||||
call c_f_pointer(ph%item,precp)
|
||||
else
|
||||
return
|
||||
end if
|
||||
|
||||
|
||||
call precp%descr()
|
||||
call flush(output_unit)
|
||||
|
||||
info = 0
|
||||
res = MLDC_ERR_FILTER(info)
|
||||
MLDC_ERR_HANDLE(res)
|
||||
return
|
||||
end function mld_c_zdescr
|
||||
|
||||
end module mld_zprec_cbind_mod
|
||||
|
||||
@@ -105,218 +105,116 @@ subroutine mld_c_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
end if
|
||||
endif
|
||||
|
||||
if (.false.) then
|
||||
if (sweeps > 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,sm%desc_data,info)
|
||||
call psb_geasb(ty,sm%desc_data,info)
|
||||
call psb_geasb(ww,sm%desc_data,info)
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
call sm%sv%apply(cone,tx,czero,ty,sm%desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%pa,y,cone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(cone,x,czero,ww,desc_data,info)
|
||||
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%pa,ty,cone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+D^(-1)*(X-A*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(cone,x,czero,ww,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%pa,ty,cone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
else
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,sm%desc_data,info)
|
||||
call psb_geasb(ty,sm%desc_data,info)
|
||||
call psb_geasb(ww,sm%desc_data,info)
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info)
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
endif
|
||||
end if
|
||||
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,sm%desc_data,info)
|
||||
call psb_geasb(ty,sm%desc_data,info)
|
||||
call psb_geasb(ww,sm%desc_data,info)
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
end if
|
||||
endif
|
||||
|
||||
|
||||
if (.not.(4*isz <= size(work))) then
|
||||
deallocate(aux,stat=info)
|
||||
endif
|
||||
|
||||
@@ -107,140 +107,89 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
end if
|
||||
endif
|
||||
|
||||
if (.true.) then
|
||||
|
||||
if (sweeps > 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 3) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
|
||||
|
||||
!
|
||||
! This is tricky. This smoother has a descriptor sm%desc_data
|
||||
! for an index space potentially different from
|
||||
! that of desc_data. Hence the size of the work vectors
|
||||
! could be wrong. We need to check and reallocate as needed.
|
||||
!
|
||||
do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
|
||||
|
||||
if (do_realloc_wv) then
|
||||
call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
|
||||
call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
end if
|
||||
|
||||
associate(tx => wv(1), ty => wv(2), ww => wv(3))
|
||||
|
||||
! Need to zero tx because of the apply_restr call.
|
||||
call tx%zero()
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
call sm%sv%apply(cone,tx,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(cone,x,czero,ww,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%pa,y,cone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(cone,x,czero,ww,desc_data,info)
|
||||
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%pa,ty,cone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+D^(-1)*(X-A*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(cone,x,czero,ww,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%pa,ty,cone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
|
||||
else
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 3) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! This is tricky. This smoother has a descriptor sm%desc_data
|
||||
! for an index space potentially different from
|
||||
! that of desc_data. Hence the size of the work vectors
|
||||
! could be wrong. We need to check and reallocate as needed.
|
||||
!
|
||||
do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
|
||||
|
||||
if (do_realloc_wv) then
|
||||
call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
|
||||
call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
end if
|
||||
|
||||
associate(tx => wv(1), ty => wv(2), ww => wv(3))
|
||||
|
||||
! Need to zero tx because of the apply_restr call.
|
||||
call tx%zero()
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info)
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
@@ -248,125 +197,47 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 3) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! This is tricky. This smoother has a descriptor sm%desc_data
|
||||
! for an index space potentially different from
|
||||
! that of desc_data. Hence the size of the work vectors
|
||||
! could be wrong. We need to check and reallocate as needed.
|
||||
!
|
||||
do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
|
||||
|
||||
if (do_realloc_wv) then
|
||||
call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
|
||||
call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
end if
|
||||
|
||||
associate(tx => wv(1), ty => wv(2), ww => wv(3))
|
||||
|
||||
! Need to zero tx because of the apply_restr call.
|
||||
call tx%zero()
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z')
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
end do
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
|
||||
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
endif
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
if (.not.(4*isz <= size(work))) then
|
||||
deallocate(aux,stat=info)
|
||||
|
||||
@@ -66,7 +66,7 @@ subroutine mld_c_as_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
|
||||
sm%pa => a
|
||||
|
||||
novr = sm%novr
|
||||
if (novr < 0) then
|
||||
info=psb_err_invalid_ovr_num_
|
||||
|
||||
@@ -77,7 +77,6 @@ subroutine mld_c_as_smoother_clone(sm,smout,info)
|
||||
allocate(smout%sv,mold=sm%sv,stat=info)
|
||||
if (info == psb_success_) call sm%sv%clone(smo%sv,info)
|
||||
end if
|
||||
smo%pa => sm%pa
|
||||
|
||||
class default
|
||||
info = psb_err_internal_error_
|
||||
|
||||
@@ -61,7 +61,6 @@ subroutine mld_c_as_smoother_free(sm,info)
|
||||
end if
|
||||
end if
|
||||
call sm%nd%free()
|
||||
sm%pa => null()
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
@@ -65,7 +65,7 @@ subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
|
||||
ictxt = desc_data%get_context()
|
||||
call psb_info(ictxt,me,np)
|
||||
|
||||
|
||||
|
||||
if (present(init)) then
|
||||
init_ = psb_toupper(init)
|
||||
else
|
||||
@@ -102,63 +102,146 @@ subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
|
||||
goto 9999
|
||||
end if
|
||||
endif
|
||||
|
||||
if (sweeps > 0) then
|
||||
call psb_geasb(tx,desc_data,info)
|
||||
call psb_geasb(ty,desc_data,info)
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
|
||||
! if .not.sv%is_iterative, there's no need to pass init
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,&
|
||||
& name,a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
endif
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
if (associated(sm%pa)) then
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ M^(-1)*(X-A*Y(j)),
|
||||
! This means we are dealing with a pure Jacobi smoother/solver.
|
||||
!
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call psb_geasb(tx,desc_data,info)
|
||||
call psb_geasb(ty,desc_data,info)
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
|
||||
! where is the diagonal and A the matrix.
|
||||
!
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(cone,tx,cone,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
end do
|
||||
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
else
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of a block-Jacobi solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,desc_data,info)
|
||||
call psb_geasb(ty,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one AXPBY and one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
call sm%sv%apply(cone,tx,cone,ty,desc_data,trans_,aux,info,init='Y')
|
||||
call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
end do
|
||||
case('Y')
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
|
||||
call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
|
||||
call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
end do
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
end if
|
||||
|
||||
|
||||
|
||||
deallocate(tx,ty,stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
@@ -167,12 +250,6 @@ subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
|
||||
@@ -49,7 +49,7 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
type(psb_c_vect_type),intent(inout) :: y
|
||||
complex(psb_spk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
integer(psb_ipk_), intent(in) :: sweeps
|
||||
integer(psb_ipk_), intent(in) :: sweeps
|
||||
complex(psb_spk_),target, intent(inout) :: work(:)
|
||||
type(psb_c_vect_type),intent(inout) :: wv(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
@@ -108,90 +108,190 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
end if
|
||||
endif
|
||||
|
||||
if (sweeps > 0) then
|
||||
if(sm%checkres) then
|
||||
call psb_geall(r,desc_data,info)
|
||||
call psb_geasb(r,desc_data,info)
|
||||
resdenum = psb_genrm2(x,desc_data,info)
|
||||
end if
|
||||
|
||||
if(sm%checkres) then
|
||||
call psb_geall(r,desc_data,info)
|
||||
call psb_geasb(r,desc_data,info)
|
||||
resdenum = psb_genrm2(x,desc_data,info)
|
||||
end if
|
||||
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
|
||||
! if .not.sv%is_iterative, there's no need to pass init
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
|
||||
|
||||
associate(tx => wv(1), ty => wv(2))
|
||||
select case (init_)
|
||||
case('Z')
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,&
|
||||
& name,a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
|
||||
else if (sweeps >= 0) then
|
||||
select type (smsv => sm%sv)
|
||||
class is (mld_c_diag_solver_type)
|
||||
!
|
||||
! This means we are dealing with a pure Jacobi smoother/solver.
|
||||
!
|
||||
associate(tx => wv(1), ty => wv(2))
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
case('Y')
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
|
||||
! where is the diagonal and A the matrix.
|
||||
!
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(cone,tx,cone,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
|
||||
call psb_geaxpby(cone,x,czero,r,r,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,r,desc_data,info)
|
||||
res = psb_genrm2(r,desc_data,info)
|
||||
if( sm%printres ) then
|
||||
call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
|
||||
end if
|
||||
if ( res < sm%tol*resdenum ) then
|
||||
if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) &
|
||||
& call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
|
||||
! where is the diagonal and A the matrix.
|
||||
!
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(cone,tx,cone,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
|
||||
call psb_geaxpby(cone,x,czero,r,r,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,r,desc_data,info)
|
||||
res = psb_genrm2(r,desc_data,info)
|
||||
if( sm%printres ) then
|
||||
call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
|
||||
end if
|
||||
if ( res < sm%tol*resdenum ) then
|
||||
if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) &
|
||||
& call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
class default
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of a block-Jacobi solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 2) then
|
||||
info = psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
associate(tx => wv(1), ty => wv(2))
|
||||
|
||||
end associate
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one AXPBY and one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
|
||||
call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
|
||||
call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
|
||||
call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
|
||||
call psb_geaxpby(cone,x,czero,r,r,desc_data,info)
|
||||
call psb_spmm(-cone,sm%pa,ty,cone,r,desc_data,info)
|
||||
res = psb_genrm2(r,desc_data,info)
|
||||
if( sm%printres ) then
|
||||
call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
|
||||
end if
|
||||
if (res < sm%tol*resdenum ) then
|
||||
if( (sm%printres).and.( mod(sm%printiter,sm%checkiter) /=0 ) ) &
|
||||
& call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
end associate
|
||||
end select
|
||||
|
||||
else
|
||||
|
||||
|
||||
@@ -71,11 +71,12 @@ subroutine mld_c_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
|
||||
nrow_a = a%get_nrows()
|
||||
nztota = a%get_nzeros()
|
||||
|
||||
sm%pa => a
|
||||
if( sm%checkres ) sm%pa => a
|
||||
|
||||
select type (smsv => sm%sv)
|
||||
class is (mld_c_diag_solver_type)
|
||||
call sm%nd%free()
|
||||
sm%pa => a
|
||||
sm%nd_nnz_tot = nztota
|
||||
call psb_sum(ictxt,sm%nd_nnz_tot)
|
||||
call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold)
|
||||
|
||||
@@ -77,7 +77,6 @@ subroutine mld_c_jac_smoother_clone(sm,smout,info)
|
||||
allocate(smout%sv,mold=sm%sv,stat=info)
|
||||
if (info == psb_success_) call sm%sv%clone(smo%sv,info)
|
||||
end if
|
||||
smo%pa => sm%pa
|
||||
|
||||
class default
|
||||
info = psb_err_internal_error_
|
||||
|
||||
@@ -72,11 +72,12 @@ subroutine mld_c_l1_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
|
||||
nrow_a = a%get_nrows()
|
||||
nztota = a%get_nzeros()
|
||||
|
||||
sm%pa => a
|
||||
if( sm%checkres ) sm%pa => a
|
||||
|
||||
select type (smsv => sm%sv)
|
||||
class is (mld_c_diag_solver_type)
|
||||
call sm%nd%free()
|
||||
sm%pa => a
|
||||
sm%nd_nnz_tot = nztota
|
||||
|
||||
call psb_sum(ictxt,sm%nd_nnz_tot)
|
||||
|
||||
@@ -105,218 +105,116 @@ subroutine mld_d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
end if
|
||||
endif
|
||||
|
||||
if (.false.) then
|
||||
if (sweeps > 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,sm%desc_data,info)
|
||||
call psb_geasb(ty,sm%desc_data,info)
|
||||
call psb_geasb(ww,sm%desc_data,info)
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
call sm%sv%apply(done,tx,dzero,ty,sm%desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%pa,y,done,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(done,x,dzero,ww,desc_data,info)
|
||||
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%pa,ty,done,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+D^(-1)*(X-A*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(done,x,dzero,ww,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%pa,ty,done,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
else
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,sm%desc_data,info)
|
||||
call psb_geasb(ty,sm%desc_data,info)
|
||||
call psb_geasb(ww,sm%desc_data,info)
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info)
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
endif
|
||||
end if
|
||||
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,sm%desc_data,info)
|
||||
call psb_geasb(ty,sm%desc_data,info)
|
||||
call psb_geasb(ww,sm%desc_data,info)
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
end if
|
||||
endif
|
||||
|
||||
|
||||
if (.not.(4*isz <= size(work))) then
|
||||
deallocate(aux,stat=info)
|
||||
endif
|
||||
|
||||
@@ -107,140 +107,89 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
end if
|
||||
endif
|
||||
|
||||
if (.true.) then
|
||||
|
||||
if (sweeps > 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 3) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
|
||||
|
||||
!
|
||||
! This is tricky. This smoother has a descriptor sm%desc_data
|
||||
! for an index space potentially different from
|
||||
! that of desc_data. Hence the size of the work vectors
|
||||
! could be wrong. We need to check and reallocate as needed.
|
||||
!
|
||||
do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
|
||||
|
||||
if (do_realloc_wv) then
|
||||
call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
|
||||
call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
end if
|
||||
|
||||
associate(tx => wv(1), ty => wv(2), ww => wv(3))
|
||||
|
||||
! Need to zero tx because of the apply_restr call.
|
||||
call tx%zero()
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
call sm%sv%apply(done,tx,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(done,x,dzero,ww,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%pa,y,done,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(done,x,dzero,ww,desc_data,info)
|
||||
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%pa,ty,done,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+D^(-1)*(X-A*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(done,x,dzero,ww,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%pa,ty,done,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
|
||||
else
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 3) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! This is tricky. This smoother has a descriptor sm%desc_data
|
||||
! for an index space potentially different from
|
||||
! that of desc_data. Hence the size of the work vectors
|
||||
! could be wrong. We need to check and reallocate as needed.
|
||||
!
|
||||
do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
|
||||
|
||||
if (do_realloc_wv) then
|
||||
call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
|
||||
call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
end if
|
||||
|
||||
associate(tx => wv(1), ty => wv(2), ww => wv(3))
|
||||
|
||||
! Need to zero tx because of the apply_restr call.
|
||||
call tx%zero()
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info)
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
@@ -248,125 +197,47 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 3) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! This is tricky. This smoother has a descriptor sm%desc_data
|
||||
! for an index space potentially different from
|
||||
! that of desc_data. Hence the size of the work vectors
|
||||
! could be wrong. We need to check and reallocate as needed.
|
||||
!
|
||||
do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
|
||||
|
||||
if (do_realloc_wv) then
|
||||
call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
|
||||
call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
end if
|
||||
|
||||
associate(tx => wv(1), ty => wv(2), ww => wv(3))
|
||||
|
||||
! Need to zero tx because of the apply_restr call.
|
||||
call tx%zero()
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z')
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
end do
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
|
||||
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
endif
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
if (.not.(4*isz <= size(work))) then
|
||||
deallocate(aux,stat=info)
|
||||
|
||||
@@ -66,7 +66,7 @@ subroutine mld_d_as_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
|
||||
sm%pa => a
|
||||
|
||||
novr = sm%novr
|
||||
if (novr < 0) then
|
||||
info=psb_err_invalid_ovr_num_
|
||||
|
||||
@@ -77,7 +77,6 @@ subroutine mld_d_as_smoother_clone(sm,smout,info)
|
||||
allocate(smout%sv,mold=sm%sv,stat=info)
|
||||
if (info == psb_success_) call sm%sv%clone(smo%sv,info)
|
||||
end if
|
||||
smo%pa => sm%pa
|
||||
|
||||
class default
|
||||
info = psb_err_internal_error_
|
||||
|
||||
@@ -61,7 +61,6 @@ subroutine mld_d_as_smoother_free(sm,info)
|
||||
end if
|
||||
end if
|
||||
call sm%nd%free()
|
||||
sm%pa => null()
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
@@ -65,7 +65,7 @@ subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
|
||||
ictxt = desc_data%get_context()
|
||||
call psb_info(ictxt,me,np)
|
||||
|
||||
|
||||
|
||||
if (present(init)) then
|
||||
init_ = psb_toupper(init)
|
||||
else
|
||||
@@ -102,63 +102,146 @@ subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
|
||||
goto 9999
|
||||
end if
|
||||
endif
|
||||
|
||||
if (sweeps > 0) then
|
||||
call psb_geasb(tx,desc_data,info)
|
||||
call psb_geasb(ty,desc_data,info)
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
|
||||
! if .not.sv%is_iterative, there's no need to pass init
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,&
|
||||
& name,a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
endif
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
if (associated(sm%pa)) then
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ M^(-1)*(X-A*Y(j)),
|
||||
! This means we are dealing with a pure Jacobi smoother/solver.
|
||||
!
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call psb_geasb(tx,desc_data,info)
|
||||
call psb_geasb(ty,desc_data,info)
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
|
||||
! where is the diagonal and A the matrix.
|
||||
!
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(done,tx,done,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
end do
|
||||
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
else
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of a block-Jacobi solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,desc_data,info)
|
||||
call psb_geasb(ty,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one AXPBY and one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
call sm%sv%apply(done,tx,done,ty,desc_data,trans_,aux,info,init='Y')
|
||||
call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
end do
|
||||
case('Y')
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
|
||||
call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
|
||||
call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
end do
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
end if
|
||||
|
||||
|
||||
|
||||
deallocate(tx,ty,stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
@@ -167,12 +250,6 @@ subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
|
||||
@@ -49,7 +49,7 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
type(psb_d_vect_type),intent(inout) :: y
|
||||
real(psb_dpk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
integer(psb_ipk_), intent(in) :: sweeps
|
||||
integer(psb_ipk_), intent(in) :: sweeps
|
||||
real(psb_dpk_),target, intent(inout) :: work(:)
|
||||
type(psb_d_vect_type),intent(inout) :: wv(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
@@ -108,90 +108,190 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
end if
|
||||
endif
|
||||
|
||||
if (sweeps > 0) then
|
||||
if(sm%checkres) then
|
||||
call psb_geall(r,desc_data,info)
|
||||
call psb_geasb(r,desc_data,info)
|
||||
resdenum = psb_genrm2(x,desc_data,info)
|
||||
end if
|
||||
|
||||
if(sm%checkres) then
|
||||
call psb_geall(r,desc_data,info)
|
||||
call psb_geasb(r,desc_data,info)
|
||||
resdenum = psb_genrm2(x,desc_data,info)
|
||||
end if
|
||||
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
|
||||
! if .not.sv%is_iterative, there's no need to pass init
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
|
||||
|
||||
associate(tx => wv(1), ty => wv(2))
|
||||
select case (init_)
|
||||
case('Z')
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,&
|
||||
& name,a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
|
||||
else if (sweeps >= 0) then
|
||||
select type (smsv => sm%sv)
|
||||
class is (mld_d_diag_solver_type)
|
||||
!
|
||||
! This means we are dealing with a pure Jacobi smoother/solver.
|
||||
!
|
||||
associate(tx => wv(1), ty => wv(2))
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
case('Y')
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
|
||||
! where is the diagonal and A the matrix.
|
||||
!
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(done,tx,done,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
|
||||
call psb_geaxpby(done,x,dzero,r,r,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,r,desc_data,info)
|
||||
res = psb_genrm2(r,desc_data,info)
|
||||
if( sm%printres ) then
|
||||
call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
|
||||
end if
|
||||
if ( res < sm%tol*resdenum ) then
|
||||
if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) &
|
||||
& call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
|
||||
! where is the diagonal and A the matrix.
|
||||
!
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(done,tx,done,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
|
||||
call psb_geaxpby(done,x,dzero,r,r,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,r,desc_data,info)
|
||||
res = psb_genrm2(r,desc_data,info)
|
||||
if( sm%printres ) then
|
||||
call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
|
||||
end if
|
||||
if ( res < sm%tol*resdenum ) then
|
||||
if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) &
|
||||
& call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
class default
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of a block-Jacobi solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 2) then
|
||||
info = psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
associate(tx => wv(1), ty => wv(2))
|
||||
|
||||
end associate
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one AXPBY and one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
|
||||
call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
|
||||
call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
|
||||
call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
|
||||
call psb_geaxpby(done,x,dzero,r,r,desc_data,info)
|
||||
call psb_spmm(-done,sm%pa,ty,done,r,desc_data,info)
|
||||
res = psb_genrm2(r,desc_data,info)
|
||||
if( sm%printres ) then
|
||||
call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
|
||||
end if
|
||||
if (res < sm%tol*resdenum ) then
|
||||
if( (sm%printres).and.( mod(sm%printiter,sm%checkiter) /=0 ) ) &
|
||||
& call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
end associate
|
||||
end select
|
||||
|
||||
else
|
||||
|
||||
|
||||
@@ -71,11 +71,12 @@ subroutine mld_d_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
|
||||
nrow_a = a%get_nrows()
|
||||
nztota = a%get_nzeros()
|
||||
|
||||
sm%pa => a
|
||||
if( sm%checkres ) sm%pa => a
|
||||
|
||||
select type (smsv => sm%sv)
|
||||
class is (mld_d_diag_solver_type)
|
||||
call sm%nd%free()
|
||||
sm%pa => a
|
||||
sm%nd_nnz_tot = nztota
|
||||
call psb_sum(ictxt,sm%nd_nnz_tot)
|
||||
call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold)
|
||||
|
||||
@@ -77,7 +77,6 @@ subroutine mld_d_jac_smoother_clone(sm,smout,info)
|
||||
allocate(smout%sv,mold=sm%sv,stat=info)
|
||||
if (info == psb_success_) call sm%sv%clone(smo%sv,info)
|
||||
end if
|
||||
smo%pa => sm%pa
|
||||
|
||||
class default
|
||||
info = psb_err_internal_error_
|
||||
|
||||
@@ -72,11 +72,12 @@ subroutine mld_d_l1_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
|
||||
nrow_a = a%get_nrows()
|
||||
nztota = a%get_nzeros()
|
||||
|
||||
sm%pa => a
|
||||
if( sm%checkres ) sm%pa => a
|
||||
|
||||
select type (smsv => sm%sv)
|
||||
class is (mld_d_diag_solver_type)
|
||||
call sm%nd%free()
|
||||
sm%pa => a
|
||||
sm%nd_nnz_tot = nztota
|
||||
|
||||
call psb_sum(ictxt,sm%nd_nnz_tot)
|
||||
|
||||
@@ -105,218 +105,116 @@ subroutine mld_s_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
end if
|
||||
endif
|
||||
|
||||
if (.false.) then
|
||||
if (sweeps > 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,sm%desc_data,info)
|
||||
call psb_geasb(ty,sm%desc_data,info)
|
||||
call psb_geasb(ww,sm%desc_data,info)
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
call sm%sv%apply(sone,tx,szero,ty,sm%desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%pa,y,sone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(sone,x,szero,ww,desc_data,info)
|
||||
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%pa,ty,sone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+D^(-1)*(X-A*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(sone,x,szero,ww,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%pa,ty,sone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
else
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,sm%desc_data,info)
|
||||
call psb_geasb(ty,sm%desc_data,info)
|
||||
call psb_geasb(ww,sm%desc_data,info)
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info)
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
endif
|
||||
end if
|
||||
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,sm%desc_data,info)
|
||||
call psb_geasb(ty,sm%desc_data,info)
|
||||
call psb_geasb(ww,sm%desc_data,info)
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
end if
|
||||
endif
|
||||
|
||||
|
||||
if (.not.(4*isz <= size(work))) then
|
||||
deallocate(aux,stat=info)
|
||||
endif
|
||||
|
||||
@@ -107,140 +107,89 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
end if
|
||||
endif
|
||||
|
||||
if (.true.) then
|
||||
|
||||
if (sweeps > 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 3) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
|
||||
|
||||
!
|
||||
! This is tricky. This smoother has a descriptor sm%desc_data
|
||||
! for an index space potentially different from
|
||||
! that of desc_data. Hence the size of the work vectors
|
||||
! could be wrong. We need to check and reallocate as needed.
|
||||
!
|
||||
do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
|
||||
|
||||
if (do_realloc_wv) then
|
||||
call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
|
||||
call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
end if
|
||||
|
||||
associate(tx => wv(1), ty => wv(2), ww => wv(3))
|
||||
|
||||
! Need to zero tx because of the apply_restr call.
|
||||
call tx%zero()
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
call sm%sv%apply(sone,tx,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(sone,x,szero,ww,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%pa,y,sone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(sone,x,szero,ww,desc_data,info)
|
||||
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%pa,ty,sone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+D^(-1)*(X-A*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(sone,x,szero,ww,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%pa,ty,sone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
|
||||
else
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 3) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! This is tricky. This smoother has a descriptor sm%desc_data
|
||||
! for an index space potentially different from
|
||||
! that of desc_data. Hence the size of the work vectors
|
||||
! could be wrong. We need to check and reallocate as needed.
|
||||
!
|
||||
do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
|
||||
|
||||
if (do_realloc_wv) then
|
||||
call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
|
||||
call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
end if
|
||||
|
||||
associate(tx => wv(1), ty => wv(2), ww => wv(3))
|
||||
|
||||
! Need to zero tx because of the apply_restr call.
|
||||
call tx%zero()
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info)
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
@@ -248,125 +197,47 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 3) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! This is tricky. This smoother has a descriptor sm%desc_data
|
||||
! for an index space potentially different from
|
||||
! that of desc_data. Hence the size of the work vectors
|
||||
! could be wrong. We need to check and reallocate as needed.
|
||||
!
|
||||
do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
|
||||
|
||||
if (do_realloc_wv) then
|
||||
call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
|
||||
call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
end if
|
||||
|
||||
associate(tx => wv(1), ty => wv(2), ww => wv(3))
|
||||
|
||||
! Need to zero tx because of the apply_restr call.
|
||||
call tx%zero()
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z')
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
end do
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
|
||||
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
endif
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
if (.not.(4*isz <= size(work))) then
|
||||
deallocate(aux,stat=info)
|
||||
|
||||
@@ -66,7 +66,7 @@ subroutine mld_s_as_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
|
||||
sm%pa => a
|
||||
|
||||
novr = sm%novr
|
||||
if (novr < 0) then
|
||||
info=psb_err_invalid_ovr_num_
|
||||
|
||||
@@ -77,7 +77,6 @@ subroutine mld_s_as_smoother_clone(sm,smout,info)
|
||||
allocate(smout%sv,mold=sm%sv,stat=info)
|
||||
if (info == psb_success_) call sm%sv%clone(smo%sv,info)
|
||||
end if
|
||||
smo%pa => sm%pa
|
||||
|
||||
class default
|
||||
info = psb_err_internal_error_
|
||||
|
||||
@@ -61,7 +61,6 @@ subroutine mld_s_as_smoother_free(sm,info)
|
||||
end if
|
||||
end if
|
||||
call sm%nd%free()
|
||||
sm%pa => null()
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
@@ -65,7 +65,7 @@ subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
|
||||
ictxt = desc_data%get_context()
|
||||
call psb_info(ictxt,me,np)
|
||||
|
||||
|
||||
|
||||
if (present(init)) then
|
||||
init_ = psb_toupper(init)
|
||||
else
|
||||
@@ -102,63 +102,146 @@ subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
|
||||
goto 9999
|
||||
end if
|
||||
endif
|
||||
|
||||
if (sweeps > 0) then
|
||||
call psb_geasb(tx,desc_data,info)
|
||||
call psb_geasb(ty,desc_data,info)
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
|
||||
! if .not.sv%is_iterative, there's no need to pass init
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,&
|
||||
& name,a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
endif
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
if (associated(sm%pa)) then
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ M^(-1)*(X-A*Y(j)),
|
||||
! This means we are dealing with a pure Jacobi smoother/solver.
|
||||
!
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call psb_geasb(tx,desc_data,info)
|
||||
call psb_geasb(ty,desc_data,info)
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
|
||||
! where is the diagonal and A the matrix.
|
||||
!
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(sone,tx,sone,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
end do
|
||||
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
else
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of a block-Jacobi solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,desc_data,info)
|
||||
call psb_geasb(ty,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one AXPBY and one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
call sm%sv%apply(sone,tx,sone,ty,desc_data,trans_,aux,info,init='Y')
|
||||
call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
end do
|
||||
case('Y')
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
|
||||
call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
|
||||
call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
end do
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
end if
|
||||
|
||||
|
||||
|
||||
deallocate(tx,ty,stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
@@ -167,12 +250,6 @@ subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
|
||||
@@ -49,7 +49,7 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
type(psb_s_vect_type),intent(inout) :: y
|
||||
real(psb_spk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
integer(psb_ipk_), intent(in) :: sweeps
|
||||
integer(psb_ipk_), intent(in) :: sweeps
|
||||
real(psb_spk_),target, intent(inout) :: work(:)
|
||||
type(psb_s_vect_type),intent(inout) :: wv(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
@@ -108,90 +108,190 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
end if
|
||||
endif
|
||||
|
||||
if (sweeps > 0) then
|
||||
if(sm%checkres) then
|
||||
call psb_geall(r,desc_data,info)
|
||||
call psb_geasb(r,desc_data,info)
|
||||
resdenum = psb_genrm2(x,desc_data,info)
|
||||
end if
|
||||
|
||||
if(sm%checkres) then
|
||||
call psb_geall(r,desc_data,info)
|
||||
call psb_geasb(r,desc_data,info)
|
||||
resdenum = psb_genrm2(x,desc_data,info)
|
||||
end if
|
||||
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
|
||||
! if .not.sv%is_iterative, there's no need to pass init
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
|
||||
|
||||
associate(tx => wv(1), ty => wv(2))
|
||||
select case (init_)
|
||||
case('Z')
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,&
|
||||
& name,a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
|
||||
else if (sweeps >= 0) then
|
||||
select type (smsv => sm%sv)
|
||||
class is (mld_s_diag_solver_type)
|
||||
!
|
||||
! This means we are dealing with a pure Jacobi smoother/solver.
|
||||
!
|
||||
associate(tx => wv(1), ty => wv(2))
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
case('Y')
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
|
||||
! where is the diagonal and A the matrix.
|
||||
!
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(sone,tx,sone,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
|
||||
call psb_geaxpby(sone,x,szero,r,r,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,r,desc_data,info)
|
||||
res = psb_genrm2(r,desc_data,info)
|
||||
if( sm%printres ) then
|
||||
call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
|
||||
end if
|
||||
if ( res < sm%tol*resdenum ) then
|
||||
if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) &
|
||||
& call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
|
||||
! where is the diagonal and A the matrix.
|
||||
!
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(sone,tx,sone,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
|
||||
call psb_geaxpby(sone,x,szero,r,r,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,r,desc_data,info)
|
||||
res = psb_genrm2(r,desc_data,info)
|
||||
if( sm%printres ) then
|
||||
call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
|
||||
end if
|
||||
if ( res < sm%tol*resdenum ) then
|
||||
if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) &
|
||||
& call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
class default
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of a block-Jacobi solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 2) then
|
||||
info = psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
associate(tx => wv(1), ty => wv(2))
|
||||
|
||||
end associate
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one AXPBY and one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
|
||||
call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
|
||||
call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
|
||||
call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
|
||||
call psb_geaxpby(sone,x,szero,r,r,desc_data,info)
|
||||
call psb_spmm(-sone,sm%pa,ty,sone,r,desc_data,info)
|
||||
res = psb_genrm2(r,desc_data,info)
|
||||
if( sm%printres ) then
|
||||
call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
|
||||
end if
|
||||
if (res < sm%tol*resdenum ) then
|
||||
if( (sm%printres).and.( mod(sm%printiter,sm%checkiter) /=0 ) ) &
|
||||
& call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
end associate
|
||||
end select
|
||||
|
||||
else
|
||||
|
||||
|
||||
@@ -71,11 +71,12 @@ subroutine mld_s_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
|
||||
nrow_a = a%get_nrows()
|
||||
nztota = a%get_nzeros()
|
||||
|
||||
sm%pa => a
|
||||
if( sm%checkres ) sm%pa => a
|
||||
|
||||
select type (smsv => sm%sv)
|
||||
class is (mld_s_diag_solver_type)
|
||||
call sm%nd%free()
|
||||
sm%pa => a
|
||||
sm%nd_nnz_tot = nztota
|
||||
call psb_sum(ictxt,sm%nd_nnz_tot)
|
||||
call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold)
|
||||
|
||||
@@ -77,7 +77,6 @@ subroutine mld_s_jac_smoother_clone(sm,smout,info)
|
||||
allocate(smout%sv,mold=sm%sv,stat=info)
|
||||
if (info == psb_success_) call sm%sv%clone(smo%sv,info)
|
||||
end if
|
||||
smo%pa => sm%pa
|
||||
|
||||
class default
|
||||
info = psb_err_internal_error_
|
||||
|
||||
@@ -72,11 +72,12 @@ subroutine mld_s_l1_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
|
||||
nrow_a = a%get_nrows()
|
||||
nztota = a%get_nzeros()
|
||||
|
||||
sm%pa => a
|
||||
if( sm%checkres ) sm%pa => a
|
||||
|
||||
select type (smsv => sm%sv)
|
||||
class is (mld_s_diag_solver_type)
|
||||
call sm%nd%free()
|
||||
sm%pa => a
|
||||
sm%nd_nnz_tot = nztota
|
||||
|
||||
call psb_sum(ictxt,sm%nd_nnz_tot)
|
||||
|
||||
@@ -105,218 +105,116 @@ subroutine mld_z_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
end if
|
||||
endif
|
||||
|
||||
if (.false.) then
|
||||
if (sweeps > 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,sm%desc_data,info)
|
||||
call psb_geasb(ty,sm%desc_data,info)
|
||||
call psb_geasb(ww,sm%desc_data,info)
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,sm%desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%pa,y,zone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(zone,x,zzero,ww,desc_data,info)
|
||||
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%pa,ty,zone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+D^(-1)*(X-A*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(zone,x,zzero,ww,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%pa,ty,zone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
else
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,sm%desc_data,info)
|
||||
call psb_geasb(ty,sm%desc_data,info)
|
||||
call psb_geasb(ww,sm%desc_data,info)
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info)
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
endif
|
||||
end if
|
||||
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,sm%desc_data,info)
|
||||
call psb_geasb(ty,sm%desc_data,info)
|
||||
call psb_geasb(ww,sm%desc_data,info)
|
||||
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
end if
|
||||
endif
|
||||
|
||||
|
||||
if (.not.(4*isz <= size(work))) then
|
||||
deallocate(aux,stat=info)
|
||||
endif
|
||||
|
||||
@@ -107,140 +107,89 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
end if
|
||||
endif
|
||||
|
||||
if (.true.) then
|
||||
|
||||
if (sweeps > 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 3) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
|
||||
|
||||
!
|
||||
! This is tricky. This smoother has a descriptor sm%desc_data
|
||||
! for an index space potentially different from
|
||||
! that of desc_data. Hence the size of the work vectors
|
||||
! could be wrong. We need to check and reallocate as needed.
|
||||
!
|
||||
do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
|
||||
|
||||
if (do_realloc_wv) then
|
||||
call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
|
||||
call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
end if
|
||||
|
||||
associate(tx => wv(1), ty => wv(2), ww => wv(3))
|
||||
|
||||
! Need to zero tx because of the apply_restr call.
|
||||
call tx%zero()
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(zone,x,zzero,ww,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%pa,y,zone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(zone,x,zzero,ww,desc_data,info)
|
||||
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%pa,ty,zone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+D^(-1)*(X-A*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(zone,x,zzero,ww,desc_data,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%pa,ty,zone,ww,desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_restr(ww,trans_,aux,info)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
|
||||
else
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 3) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! This is tricky. This smoother has a descriptor sm%desc_data
|
||||
! for an index space potentially different from
|
||||
! that of desc_data. Hence the size of the work vectors
|
||||
! could be wrong. We need to check and reallocate as needed.
|
||||
!
|
||||
do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
|
||||
|
||||
if (do_realloc_wv) then
|
||||
call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
|
||||
call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
end if
|
||||
|
||||
associate(tx => wv(1), ty => wv(2), ww => wv(3))
|
||||
|
||||
! Need to zero tx because of the apply_restr call.
|
||||
call tx%zero()
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info)
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
|
||||
!
|
||||
! Shortcut: in this case there is nothing else to be done.
|
||||
!
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
@@ -248,125 +197,47 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of an AS solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 3) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! This is tricky. This smoother has a descriptor sm%desc_data
|
||||
! for an index space potentially different from
|
||||
! that of desc_data. Hence the size of the work vectors
|
||||
! could be wrong. We need to check and reallocate as needed.
|
||||
!
|
||||
do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.&
|
||||
& (wv(3)%get_nrows() < sm%desc_data%get_local_cols())
|
||||
|
||||
if (do_realloc_wv) then
|
||||
call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v)
|
||||
call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v)
|
||||
end if
|
||||
|
||||
associate(tx => wv(1), ty => wv(2), ww => wv(3))
|
||||
|
||||
! Need to zero tx because of the apply_restr call.
|
||||
call tx%zero()
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
|
||||
if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
|
||||
select case (init_)
|
||||
case('Z')
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z')
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
|
||||
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
end do
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info)
|
||||
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
|
||||
& work=aux,trans=trans_)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
|
||||
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
|
||||
endif
|
||||
|
||||
end do
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
!
|
||||
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
|
||||
!
|
||||
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
call psb_errpush(info,name,&
|
||||
& i_err=(/itwo,sweeps,izero,izero,izero/))
|
||||
goto 9999
|
||||
|
||||
endif
|
||||
end if
|
||||
|
||||
if (.not.(4*isz <= size(work))) then
|
||||
deallocate(aux,stat=info)
|
||||
|
||||
@@ -66,7 +66,7 @@ subroutine mld_z_as_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
|
||||
sm%pa => a
|
||||
|
||||
novr = sm%novr
|
||||
if (novr < 0) then
|
||||
info=psb_err_invalid_ovr_num_
|
||||
|
||||
@@ -77,7 +77,6 @@ subroutine mld_z_as_smoother_clone(sm,smout,info)
|
||||
allocate(smout%sv,mold=sm%sv,stat=info)
|
||||
if (info == psb_success_) call sm%sv%clone(smo%sv,info)
|
||||
end if
|
||||
smo%pa => sm%pa
|
||||
|
||||
class default
|
||||
info = psb_err_internal_error_
|
||||
|
||||
@@ -61,7 +61,6 @@ subroutine mld_z_as_smoother_free(sm,info)
|
||||
end if
|
||||
end if
|
||||
call sm%nd%free()
|
||||
sm%pa => null()
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
@@ -65,7 +65,7 @@ subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
|
||||
ictxt = desc_data%get_context()
|
||||
call psb_info(ictxt,me,np)
|
||||
|
||||
|
||||
|
||||
if (present(init)) then
|
||||
init_ = psb_toupper(init)
|
||||
else
|
||||
@@ -102,63 +102,146 @@ subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
|
||||
goto 9999
|
||||
end if
|
||||
endif
|
||||
|
||||
if (sweeps > 0) then
|
||||
call psb_geasb(tx,desc_data,info)
|
||||
call psb_geasb(ty,desc_data,info)
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
|
||||
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
|
||||
! if .not.sv%is_iterative, there's no need to pass init
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,&
|
||||
& name,a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
endif
|
||||
|
||||
else if (sweeps >= 0) then
|
||||
if (associated(sm%pa)) then
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ M^(-1)*(X-A*Y(j)),
|
||||
! This means we are dealing with a pure Jacobi smoother/solver.
|
||||
!
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call psb_geasb(tx,desc_data,info)
|
||||
call psb_geasb(ty,desc_data,info)
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
|
||||
! where is the diagonal and A the matrix.
|
||||
!
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(zone,tx,zone,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
end do
|
||||
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
else
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of a block-Jacobi solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
call psb_geasb(tx,desc_data,info)
|
||||
call psb_geasb(ty,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one AXPBY and one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
call sm%sv%apply(zone,tx,zone,ty,desc_data,trans_,aux,info,init='Y')
|
||||
call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,info,init='Z')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
end do
|
||||
case('Y')
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
|
||||
call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
|
||||
call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
end do
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
end if
|
||||
|
||||
|
||||
|
||||
deallocate(tx,ty,stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
@@ -167,12 +250,6 @@ subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
|
||||
else
|
||||
|
||||
info = psb_err_iarg_neg_
|
||||
|
||||
@@ -49,7 +49,7 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
type(psb_z_vect_type),intent(inout) :: y
|
||||
complex(psb_dpk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
integer(psb_ipk_), intent(in) :: sweeps
|
||||
integer(psb_ipk_), intent(in) :: sweeps
|
||||
complex(psb_dpk_),target, intent(inout) :: work(:)
|
||||
type(psb_z_vect_type),intent(inout) :: wv(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
@@ -108,90 +108,190 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
||||
end if
|
||||
endif
|
||||
|
||||
if (sweeps > 0) then
|
||||
if(sm%checkres) then
|
||||
call psb_geall(r,desc_data,info)
|
||||
call psb_geasb(r,desc_data,info)
|
||||
resdenum = psb_genrm2(x,desc_data,info)
|
||||
end if
|
||||
|
||||
if(sm%checkres) then
|
||||
call psb_geall(r,desc_data,info)
|
||||
call psb_geasb(r,desc_data,info)
|
||||
resdenum = psb_genrm2(x,desc_data,info)
|
||||
end if
|
||||
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
|
||||
! if .not.sv%is_iterative, there's no need to pass init
|
||||
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
|
||||
|
||||
associate(tx => wv(1), ty => wv(2))
|
||||
select case (init_)
|
||||
case('Z')
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,&
|
||||
& name,a_err='Error in sub_aply Jacobi Sweeps = 1')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
|
||||
else if (sweeps >= 0) then
|
||||
select type (smsv => sm%sv)
|
||||
class is (mld_z_diag_solver_type)
|
||||
!
|
||||
! This means we are dealing with a pure Jacobi smoother/solver.
|
||||
!
|
||||
associate(tx => wv(1), ty => wv(2))
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
case('Y')
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
|
||||
! where is the diagonal and A the matrix.
|
||||
!
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(zone,tx,zone,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
|
||||
call psb_geaxpby(zone,x,zzero,r,r,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,r,desc_data,info)
|
||||
res = psb_genrm2(r,desc_data,info)
|
||||
if( sm%printres ) then
|
||||
call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
|
||||
end if
|
||||
if ( res < sm%tol*resdenum ) then
|
||||
if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) &
|
||||
& call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
|
||||
! where is the diagonal and A the matrix.
|
||||
!
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(zone,tx,zone,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
|
||||
call psb_geaxpby(zone,x,zzero,r,r,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,r,desc_data,info)
|
||||
res = psb_genrm2(r,desc_data,info)
|
||||
if( sm%printres ) then
|
||||
call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
|
||||
end if
|
||||
if ( res < sm%tol*resdenum ) then
|
||||
if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) &
|
||||
& call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
end associate
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
class default
|
||||
!
|
||||
!
|
||||
! Apply multiple sweeps of a block-Jacobi solver
|
||||
! to compute an approximate solution of a linear system.
|
||||
!
|
||||
!
|
||||
if (size(wv) < 2) then
|
||||
info = psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
& a_err='invalid wv size in smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
associate(tx => wv(1), ty => wv(2))
|
||||
|
||||
end associate
|
||||
!
|
||||
! Unroll the first iteration and fold it inside SELECT CASE
|
||||
! this will save one AXPBY and one SPMM when INIT=Z, and will be
|
||||
! significant when sweeps=1 (a common case)
|
||||
!
|
||||
select case (init_)
|
||||
case('Z')
|
||||
|
||||
else if (sweeps == 0) then
|
||||
!
|
||||
! 0 sweeps of smoother is the identity operator
|
||||
!
|
||||
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
|
||||
call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
|
||||
|
||||
case('Y')
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
|
||||
call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case('U')
|
||||
if (.not.present(initu)) then
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='missing initu to smoother_apply')
|
||||
goto 9999
|
||||
end if
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
|
||||
call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='wrong init to smoother_apply')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
do i=1, sweeps-1
|
||||
!
|
||||
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
|
||||
! block diagonal part and the remaining part of the local matrix
|
||||
! and Y(j) is the approximate solution at sweep j.
|
||||
!
|
||||
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
|
||||
call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
|
||||
|
||||
if (info /= psb_success_) exit
|
||||
|
||||
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
|
||||
call psb_geaxpby(zone,x,zzero,r,r,desc_data,info)
|
||||
call psb_spmm(-zone,sm%pa,ty,zone,r,desc_data,info)
|
||||
res = psb_genrm2(r,desc_data,info)
|
||||
if( sm%printres ) then
|
||||
call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
|
||||
end if
|
||||
if (res < sm%tol*resdenum ) then
|
||||
if( (sm%printres).and.( mod(sm%printiter,sm%checkiter) /=0 ) ) &
|
||||
& call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_internal_error_
|
||||
call psb_errpush(info,name,&
|
||||
& a_err='subsolve with Jacobi sweeps > 1')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
end associate
|
||||
end select
|
||||
|
||||
else
|
||||
|
||||
|
||||
@@ -71,11 +71,12 @@ subroutine mld_z_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
|
||||
nrow_a = a%get_nrows()
|
||||
nztota = a%get_nzeros()
|
||||
|
||||
sm%pa => a
|
||||
if( sm%checkres ) sm%pa => a
|
||||
|
||||
select type (smsv => sm%sv)
|
||||
class is (mld_z_diag_solver_type)
|
||||
call sm%nd%free()
|
||||
sm%pa => a
|
||||
sm%nd_nnz_tot = nztota
|
||||
call psb_sum(ictxt,sm%nd_nnz_tot)
|
||||
call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold)
|
||||
|
||||
@@ -77,7 +77,6 @@ subroutine mld_z_jac_smoother_clone(sm,smout,info)
|
||||
allocate(smout%sv,mold=sm%sv,stat=info)
|
||||
if (info == psb_success_) call sm%sv%clone(smo%sv,info)
|
||||
end if
|
||||
smo%pa => sm%pa
|
||||
|
||||
class default
|
||||
info = psb_err_internal_error_
|
||||
|
||||
@@ -72,11 +72,12 @@ subroutine mld_z_l1_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
|
||||
nrow_a = a%get_nrows()
|
||||
nztota = a%get_nzeros()
|
||||
|
||||
sm%pa => a
|
||||
if( sm%checkres ) sm%pa => a
|
||||
|
||||
select type (smsv => sm%sv)
|
||||
class is (mld_z_diag_solver_type)
|
||||
call sm%nd%free()
|
||||
sm%pa => a
|
||||
sm%nd_nnz_tot = nztota
|
||||
|
||||
call psb_sum(ictxt,sm%nd_nnz_tot)
|
||||
|
||||
@@ -65,7 +65,6 @@ module mld_c_as_smoother
|
||||
! parent type.
|
||||
! class(mld_c_base_solver_type), allocatable :: sv
|
||||
!
|
||||
type(psb_cspmat_type), pointer :: pa => null()
|
||||
type(psb_cspmat_type) :: nd
|
||||
type(psb_desc_type) :: desc_data
|
||||
integer(psb_ipk_) :: novr, restr, prol
|
||||
|
||||
@@ -65,7 +65,6 @@ module mld_d_as_smoother
|
||||
! parent type.
|
||||
! class(mld_d_base_solver_type), allocatable :: sv
|
||||
!
|
||||
type(psb_dspmat_type), pointer :: pa => null()
|
||||
type(psb_dspmat_type) :: nd
|
||||
type(psb_desc_type) :: desc_data
|
||||
integer(psb_ipk_) :: novr, restr, prol
|
||||
|
||||
@@ -65,7 +65,6 @@ module mld_s_as_smoother
|
||||
! parent type.
|
||||
! class(mld_s_base_solver_type), allocatable :: sv
|
||||
!
|
||||
type(psb_sspmat_type), pointer :: pa => null()
|
||||
type(psb_sspmat_type) :: nd
|
||||
type(psb_desc_type) :: desc_data
|
||||
integer(psb_ipk_) :: novr, restr, prol
|
||||
|
||||
@@ -65,7 +65,6 @@ module mld_z_as_smoother
|
||||
! parent type.
|
||||
! class(mld_z_base_solver_type), allocatable :: sv
|
||||
!
|
||||
type(psb_zspmat_type), pointer :: pa => null()
|
||||
type(psb_zspmat_type) :: nd
|
||||
type(psb_desc_type) :: desc_data
|
||||
integer(psb_ipk_) :: novr, restr, prol
|
||||
|
||||
Reference in New Issue
Block a user