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

mld2p4-2:

Browse Source 
mlprec/impl/smoother/mld_c_jac_smoother_apply.f90
 mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90
 mlprec/impl/smoother/mld_d_jac_smoother_apply.f90
 mlprec/impl/smoother/mld_d_jac_smoother_apply_vect.f90
 mlprec/impl/smoother/mld_s_jac_smoother_apply.f90
 mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90
 mlprec/impl/smoother/mld_z_jac_smoother_apply.f90
 mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90

Fixed application for various cases SWEEPS >=1.
stopcriterion
Salvatore Filippone 9 years ago
parent 4d3089901a
commit da6bde2130
8 changed files with 364 additions and 388 deletions
Show Stats Download Patch File Download Diff File

88
mlprec/impl/smoother/mld_c_jac_smoother_apply.f90
Unescape Escape View File

@ -55,7 +55,7 @@ subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
integer(psb_ipk_) :: n_row,n_col integer(psb_ipk_) :: n_row,n_col
complex(psb_spk_), allocatable :: tx(:),ty(:) complex(psb_spk_), allocatable :: tx(:),ty(:)
complex(psb_spk_), pointer :: ww(:), aux(:) complex(psb_spk_), pointer :: aux(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act integer(psb_ipk_) :: ictxt,np,me,i, err_act
character :: trans_, init_ character :: trans_, init_
character(len=20) :: name='c_jac_smoother_apply' character(len=20) :: name='c_jac_smoother_apply'
@ -89,71 +89,73 @@ subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
n_row = desc_data%get_local_rows() n_row = desc_data%get_local_rows()
n_col = desc_data%get_local_cols() n_col = desc_data%get_local_cols()
if (n_col <= size(work)) then if (4*n_col <= size(work)) then
ww => work(1:n_col) aux => work(:)
if ((4*n_col+n_col) <= size(work)) then
aux => work(n_col+1:)
else
allocate(aux(4*n_col),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,&
& i_err=(/4*n_col,izero,izero,izero,izero/),&
& a_err='complex(psb_spk_)')
goto 9999
end if
endif
else else
allocate(ww(n_col),aux(4*n_col),stat=info) allocate(aux(4*n_col),stat=info)
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_alloc_request_ info=psb_err_alloc_request_
call psb_errpush(info,name,& call psb_errpush(info,name,&
& i_err=(/5*n_col,izero,izero,izero,izero/),& & i_err=(/4*n_col,izero,izero,izero,izero/),&
& a_err='complex(psb_spk_)') & a_err='complex(psb_spk_)')
goto 9999 goto 9999
end if end if
endif endif
!!$ if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then
!!$ ! if .not.sv%is_iterative, there's no need to pass init ! 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) call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
!!$
!!$ if (info /= psb_success_) then if (info /= psb_success_) then
!!$ call psb_errpush(psb_err_internal_error_,& call psb_errpush(psb_err_internal_error_,&
!!$ & name,a_err='Error in sub_aply Jacobi Sweeps = 1') & name,a_err='Error in sub_aply Jacobi Sweeps = 1')
!!$ goto 9999 goto 9999
!!$ endif endif
!!$
!!$ else if (sweeps >= 1) then else if (sweeps >= 0) then
! !
! !
! Apply multiple sweeps of a block-Jacobi solver ! Apply multiple sweeps of a block-Jacobi solver
! to compute an approximate solution of a linear system. ! to compute an approximate solution of a linear system.
! !
! !
call psb_geasb(tx,desc_data,info) call psb_geasb(tx,desc_data,info)
call psb_geasb(ty,desc_data,info) call psb_geasb(ty,desc_data,info)
!
! 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_) select case (init_)
case('Z') case('Z')
ty(:) = czero
call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,info,init='Z')
case('Y') case('Y')
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
call psb_geaxpby(cone,y,czero,ty,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') case('U')
if (.not.present(initu)) then if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply') & a_err='missing initu to smoother_apply')
goto 9999 goto 9999
end if end if
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
call psb_geaxpby(cone,initu,czero,ty,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')
case default case default
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply') & a_err='wrong init to smoother_apply')
goto 9999 goto 9999
end select end select
do i=1, sweeps do i=1, sweeps-1
! !
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the ! 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 ! block diagonal part and the remaining part of the local matrix
@ -186,23 +188,17 @@ subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
goto 9999 goto 9999
end if end if
!!$ else else
!!$
!!$ info = psb_err_iarg_neg_
!!$ call psb_errpush(info,name,&
!!$ & i_err=(/itwo,sweeps,izero,izero,izero/))
!!$ goto 9999
!!$
!!$ endif
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
if (n_col <= size(work)) then endif
if ((4*n_col+n_col) <= size(work)) then
else if (.not.(4*n_col <= size(work))) then
deallocate(aux) deallocate(aux)
endif
else
deallocate(ww,aux)
endif endif
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)

92
mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90
Unescape Escape View File

@ -56,10 +56,10 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
integer(psb_ipk_) :: n_row,n_col integer(psb_ipk_) :: n_row,n_col
type(psb_c_vect_type) :: tx, ty type(psb_c_vect_type) :: tx, ty
complex(psb_spk_), pointer :: ww(:), aux(:) complex(psb_spk_), pointer :: aux(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act integer(psb_ipk_) :: ictxt,np,me,i, err_act
character :: trans_, init_ character :: trans_, init_
character(len=20) :: name='c_jac_smoother_apply' character(len=20) :: name='c_jac_smoother_apply_v'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -90,71 +90,74 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
n_row = desc_data%get_local_rows() n_row = desc_data%get_local_rows()
n_col = desc_data%get_local_cols() n_col = desc_data%get_local_cols()
if (n_col <= size(work)) then if (4*n_col <= size(work)) then
ww => work(1:n_col) aux => work(:)
if ((4*n_col+n_col) <= size(work)) then
aux => work(n_col+1:)
else
allocate(aux(4*n_col),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,&
& i_err=(/4*n_col,izero,izero,izero,izero/),&
& a_err='complex(psb_spk_)')
goto 9999
end if
endif
else else
allocate(ww(n_col),aux(4*n_col),stat=info) allocate(aux(4*n_col),stat=info)
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_alloc_request_ info=psb_err_alloc_request_
call psb_errpush(info,name,& call psb_errpush(info,name,&
& i_err=(/5*n_col,izero,izero,izero,izero/),& & i_err=(/4*n_col,izero,izero,izero,izero/),&
& a_err='complex(psb_spk_)') & a_err='complex(psb_spk_)')
goto 9999 goto 9999
end if end if
endif endif
!!$ if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then
!!$ ! if .not.sv%is_iterative, there's no need to pass init ! 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) call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
!!$
!!$ if (info /= psb_success_) then if (info /= psb_success_) then
!!$ call psb_errpush(psb_err_internal_error_,& call psb_errpush(psb_err_internal_error_,&
!!$ & name,a_err='Error in sub_aply Jacobi Sweeps = 1') & name,a_err='Error in sub_aply Jacobi Sweeps = 1')
!!$ goto 9999 goto 9999
!!$ endif endif
!!$
!!$ else if (sweeps >= 0) then else if (sweeps >= 0) then
! !
! !
! Apply multiple sweeps of a block-Jacobi solver ! Apply multiple sweeps of a block-Jacobi solver
! to compute an approximate solution of a linear system. ! to compute an approximate solution of a linear system.
! !
! !
call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.) call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.) call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.)
!
! 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_) select case (init_)
case('Z') case('Z')
call ty%zero()
call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,info,init='Z')
case('Y') case('Y')
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
call psb_geaxpby(cone,y,czero,ty,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') case('U')
if (.not.present(initu)) then if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply') & a_err='missing initu to smoother_apply')
goto 9999 goto 9999
end if end if
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
call psb_geaxpby(cone,initu,czero,ty,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')
case default case default
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply') & a_err='wrong init to smoother_apply')
goto 9999 goto 9999
end select end select
do i=1, sweeps do i=1, sweeps-1
! !
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the ! 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 ! block diagonal part and the remaining part of the local matrix
@ -188,22 +191,17 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
goto 9999 goto 9999
end if end if
!!$ else
!!$
!!$ info = psb_err_iarg_neg_
!!$ call psb_errpush(info,name,&
!!$ & i_err=(/itwo,sweeps,izero,izero,izero/))
!!$ goto 9999
!!$
!!$ endif
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
deallocate(aux)
endif
else else
deallocate(ww,aux)
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
endif
if (.not.(4*n_col <= size(work))) then
deallocate(aux)
endif endif
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)

88
mlprec/impl/smoother/mld_d_jac_smoother_apply.f90
Unescape Escape View File

@ -55,7 +55,7 @@ subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
integer(psb_ipk_) :: n_row,n_col integer(psb_ipk_) :: n_row,n_col
real(psb_dpk_), allocatable :: tx(:),ty(:) real(psb_dpk_), allocatable :: tx(:),ty(:)
real(psb_dpk_), pointer :: ww(:), aux(:) real(psb_dpk_), pointer :: aux(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act integer(psb_ipk_) :: ictxt,np,me,i, err_act
character :: trans_, init_ character :: trans_, init_
character(len=20) :: name='d_jac_smoother_apply' character(len=20) :: name='d_jac_smoother_apply'
@ -89,71 +89,73 @@ subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
n_row = desc_data%get_local_rows() n_row = desc_data%get_local_rows()
n_col = desc_data%get_local_cols() n_col = desc_data%get_local_cols()
if (n_col <= size(work)) then if (4*n_col <= size(work)) then
ww => work(1:n_col) aux => work(:)
if ((4*n_col+n_col) <= size(work)) then
aux => work(n_col+1:)
else
allocate(aux(4*n_col),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,&
& i_err=(/4*n_col,izero,izero,izero,izero/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
endif
else else
allocate(ww(n_col),aux(4*n_col),stat=info) allocate(aux(4*n_col),stat=info)
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_alloc_request_ info=psb_err_alloc_request_
call psb_errpush(info,name,& call psb_errpush(info,name,&
& i_err=(/5*n_col,izero,izero,izero,izero/),& & i_err=(/4*n_col,izero,izero,izero,izero/),&
& a_err='real(psb_dpk_)') & a_err='real(psb_dpk_)')
goto 9999 goto 9999
end if end if
endif endif
!!$ if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then
!!$ ! if .not.sv%is_iterative, there's no need to pass init ! 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) call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
!!$
!!$ if (info /= psb_success_) then if (info /= psb_success_) then
!!$ call psb_errpush(psb_err_internal_error_,& call psb_errpush(psb_err_internal_error_,&
!!$ & name,a_err='Error in sub_aply Jacobi Sweeps = 1') & name,a_err='Error in sub_aply Jacobi Sweeps = 1')
!!$ goto 9999 goto 9999
!!$ endif endif
!!$
!!$ else if (sweeps >= 1) then else if (sweeps >= 0) then
! !
! !
! Apply multiple sweeps of a block-Jacobi solver ! Apply multiple sweeps of a block-Jacobi solver
! to compute an approximate solution of a linear system. ! to compute an approximate solution of a linear system.
! !
! !
call psb_geasb(tx,desc_data,info) call psb_geasb(tx,desc_data,info)
call psb_geasb(ty,desc_data,info) call psb_geasb(ty,desc_data,info)
!
! 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_) select case (init_)
case('Z') case('Z')
ty(:) = dzero
call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,info,init='Z')
case('Y') case('Y')
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
call psb_geaxpby(done,y,dzero,ty,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') case('U')
if (.not.present(initu)) then if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply') & a_err='missing initu to smoother_apply')
goto 9999 goto 9999
end if end if
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
call psb_geaxpby(done,initu,dzero,ty,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')
case default case default
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply') & a_err='wrong init to smoother_apply')
goto 9999 goto 9999
end select end select
do i=1, sweeps do i=1, sweeps-1
! !
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the ! 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 ! block diagonal part and the remaining part of the local matrix
@ -186,23 +188,17 @@ subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
goto 9999 goto 9999
end if end if
!!$ else else
!!$
!!$ info = psb_err_iarg_neg_
!!$ call psb_errpush(info,name,&
!!$ & i_err=(/itwo,sweeps,izero,izero,izero/))
!!$ goto 9999
!!$
!!$ endif
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
if (n_col <= size(work)) then endif
if ((4*n_col+n_col) <= size(work)) then
else if (.not.(4*n_col <= size(work))) then
deallocate(aux) deallocate(aux)
endif
else
deallocate(ww,aux)
endif endif
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)

92
mlprec/impl/smoother/mld_d_jac_smoother_apply_vect.f90
Unescape Escape View File

@ -56,10 +56,10 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
integer(psb_ipk_) :: n_row,n_col integer(psb_ipk_) :: n_row,n_col
type(psb_d_vect_type) :: tx, ty type(psb_d_vect_type) :: tx, ty
real(psb_dpk_), pointer :: ww(:), aux(:) real(psb_dpk_), pointer :: aux(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act integer(psb_ipk_) :: ictxt,np,me,i, err_act
character :: trans_, init_ character :: trans_, init_
character(len=20) :: name='d_jac_smoother_apply' character(len=20) :: name='d_jac_smoother_apply_v'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -90,71 +90,74 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
n_row = desc_data%get_local_rows() n_row = desc_data%get_local_rows()
n_col = desc_data%get_local_cols() n_col = desc_data%get_local_cols()
if (n_col <= size(work)) then if (4*n_col <= size(work)) then
ww => work(1:n_col) aux => work(:)
if ((4*n_col+n_col) <= size(work)) then
aux => work(n_col+1:)
else
allocate(aux(4*n_col),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,&
& i_err=(/4*n_col,izero,izero,izero,izero/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
endif
else else
allocate(ww(n_col),aux(4*n_col),stat=info) allocate(aux(4*n_col),stat=info)
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_alloc_request_ info=psb_err_alloc_request_
call psb_errpush(info,name,& call psb_errpush(info,name,&
& i_err=(/5*n_col,izero,izero,izero,izero/),& & i_err=(/4*n_col,izero,izero,izero,izero/),&
& a_err='real(psb_dpk_)') & a_err='real(psb_dpk_)')
goto 9999 goto 9999
end if end if
endif endif
!!$ if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then
!!$ ! if .not.sv%is_iterative, there's no need to pass init ! 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) call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
!!$
!!$ if (info /= psb_success_) then if (info /= psb_success_) then
!!$ call psb_errpush(psb_err_internal_error_,& call psb_errpush(psb_err_internal_error_,&
!!$ & name,a_err='Error in sub_aply Jacobi Sweeps = 1') & name,a_err='Error in sub_aply Jacobi Sweeps = 1')
!!$ goto 9999 goto 9999
!!$ endif endif
!!$
!!$ else if (sweeps >= 0) then else if (sweeps >= 0) then
! !
! !
! Apply multiple sweeps of a block-Jacobi solver ! Apply multiple sweeps of a block-Jacobi solver
! to compute an approximate solution of a linear system. ! to compute an approximate solution of a linear system.
! !
! !
call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.) call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.) call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.)
!
! 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_) select case (init_)
case('Z') case('Z')
call ty%zero()
call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,info,init='Z')
case('Y') case('Y')
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
call psb_geaxpby(done,y,dzero,ty,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') case('U')
if (.not.present(initu)) then if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply') & a_err='missing initu to smoother_apply')
goto 9999 goto 9999
end if end if
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
call psb_geaxpby(done,initu,dzero,ty,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')
case default case default
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply') & a_err='wrong init to smoother_apply')
goto 9999 goto 9999
end select end select
do i=1, sweeps do i=1, sweeps-1
! !
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the ! 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 ! block diagonal part and the remaining part of the local matrix
@ -188,22 +191,17 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
goto 9999 goto 9999
end if end if
!!$ else
!!$
!!$ info = psb_err_iarg_neg_
!!$ call psb_errpush(info,name,&
!!$ & i_err=(/itwo,sweeps,izero,izero,izero/))
!!$ goto 9999
!!$
!!$ endif
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
deallocate(aux)
endif
else else
deallocate(ww,aux)
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
endif
if (.not.(4*n_col <= size(work))) then
deallocate(aux)
endif endif
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)

88
mlprec/impl/smoother/mld_s_jac_smoother_apply.f90
Unescape Escape View File

@ -55,7 +55,7 @@ subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
integer(psb_ipk_) :: n_row,n_col integer(psb_ipk_) :: n_row,n_col
real(psb_spk_), allocatable :: tx(:),ty(:) real(psb_spk_), allocatable :: tx(:),ty(:)
real(psb_spk_), pointer :: ww(:), aux(:) real(psb_spk_), pointer :: aux(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act integer(psb_ipk_) :: ictxt,np,me,i, err_act
character :: trans_, init_ character :: trans_, init_
character(len=20) :: name='s_jac_smoother_apply' character(len=20) :: name='s_jac_smoother_apply'
@ -89,71 +89,73 @@ subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
n_row = desc_data%get_local_rows() n_row = desc_data%get_local_rows()
n_col = desc_data%get_local_cols() n_col = desc_data%get_local_cols()
if (n_col <= size(work)) then if (4*n_col <= size(work)) then
ww => work(1:n_col) aux => work(:)
if ((4*n_col+n_col) <= size(work)) then
aux => work(n_col+1:)
else
allocate(aux(4*n_col),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,&
& i_err=(/4*n_col,izero,izero,izero,izero/),&
& a_err='real(psb_spk_)')
goto 9999
end if
endif
else else
allocate(ww(n_col),aux(4*n_col),stat=info) allocate(aux(4*n_col),stat=info)
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_alloc_request_ info=psb_err_alloc_request_
call psb_errpush(info,name,& call psb_errpush(info,name,&
& i_err=(/5*n_col,izero,izero,izero,izero/),& & i_err=(/4*n_col,izero,izero,izero,izero/),&
& a_err='real(psb_spk_)') & a_err='real(psb_spk_)')
goto 9999 goto 9999
end if end if
endif endif
!!$ if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then
!!$ ! if .not.sv%is_iterative, there's no need to pass init ! 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) call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
!!$
!!$ if (info /= psb_success_) then if (info /= psb_success_) then
!!$ call psb_errpush(psb_err_internal_error_,& call psb_errpush(psb_err_internal_error_,&
!!$ & name,a_err='Error in sub_aply Jacobi Sweeps = 1') & name,a_err='Error in sub_aply Jacobi Sweeps = 1')
!!$ goto 9999 goto 9999
!!$ endif endif
!!$
!!$ else if (sweeps >= 1) then else if (sweeps >= 0) then
! !
! !
! Apply multiple sweeps of a block-Jacobi solver ! Apply multiple sweeps of a block-Jacobi solver
! to compute an approximate solution of a linear system. ! to compute an approximate solution of a linear system.
! !
! !
call psb_geasb(tx,desc_data,info) call psb_geasb(tx,desc_data,info)
call psb_geasb(ty,desc_data,info) call psb_geasb(ty,desc_data,info)
!
! 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_) select case (init_)
case('Z') case('Z')
ty(:) = szero
call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,info,init='Z')
case('Y') case('Y')
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
call psb_geaxpby(sone,y,szero,ty,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') case('U')
if (.not.present(initu)) then if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply') & a_err='missing initu to smoother_apply')
goto 9999 goto 9999
end if end if
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
call psb_geaxpby(sone,initu,szero,ty,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')
case default case default
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply') & a_err='wrong init to smoother_apply')
goto 9999 goto 9999
end select end select
do i=1, sweeps do i=1, sweeps-1
! !
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the ! 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 ! block diagonal part and the remaining part of the local matrix
@ -186,23 +188,17 @@ subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
goto 9999 goto 9999
end if end if
!!$ else else
!!$
!!$ info = psb_err_iarg_neg_
!!$ call psb_errpush(info,name,&
!!$ & i_err=(/itwo,sweeps,izero,izero,izero/))
!!$ goto 9999
!!$
!!$ endif
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
if (n_col <= size(work)) then endif
if ((4*n_col+n_col) <= size(work)) then
else if (.not.(4*n_col <= size(work))) then
deallocate(aux) deallocate(aux)
endif
else
deallocate(ww,aux)
endif endif
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)

92
mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90
Unescape Escape View File

@ -56,10 +56,10 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
integer(psb_ipk_) :: n_row,n_col integer(psb_ipk_) :: n_row,n_col
type(psb_s_vect_type) :: tx, ty type(psb_s_vect_type) :: tx, ty
real(psb_spk_), pointer :: ww(:), aux(:) real(psb_spk_), pointer :: aux(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act integer(psb_ipk_) :: ictxt,np,me,i, err_act
character :: trans_, init_ character :: trans_, init_
character(len=20) :: name='s_jac_smoother_apply' character(len=20) :: name='s_jac_smoother_apply_v'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -90,71 +90,74 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
n_row = desc_data%get_local_rows() n_row = desc_data%get_local_rows()
n_col = desc_data%get_local_cols() n_col = desc_data%get_local_cols()
if (n_col <= size(work)) then if (4*n_col <= size(work)) then
ww => work(1:n_col) aux => work(:)
if ((4*n_col+n_col) <= size(work)) then
aux => work(n_col+1:)
else
allocate(aux(4*n_col),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,&
& i_err=(/4*n_col,izero,izero,izero,izero/),&
& a_err='real(psb_spk_)')
goto 9999
end if
endif
else else
allocate(ww(n_col),aux(4*n_col),stat=info) allocate(aux(4*n_col),stat=info)
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_alloc_request_ info=psb_err_alloc_request_
call psb_errpush(info,name,& call psb_errpush(info,name,&
& i_err=(/5*n_col,izero,izero,izero,izero/),& & i_err=(/4*n_col,izero,izero,izero,izero/),&
& a_err='real(psb_spk_)') & a_err='real(psb_spk_)')
goto 9999 goto 9999
end if end if
endif endif
!!$ if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then
!!$ ! if .not.sv%is_iterative, there's no need to pass init ! 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) call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
!!$
!!$ if (info /= psb_success_) then if (info /= psb_success_) then
!!$ call psb_errpush(psb_err_internal_error_,& call psb_errpush(psb_err_internal_error_,&
!!$ & name,a_err='Error in sub_aply Jacobi Sweeps = 1') & name,a_err='Error in sub_aply Jacobi Sweeps = 1')
!!$ goto 9999 goto 9999
!!$ endif endif
!!$
!!$ else if (sweeps >= 0) then else if (sweeps >= 0) then
! !
! !
! Apply multiple sweeps of a block-Jacobi solver ! Apply multiple sweeps of a block-Jacobi solver
! to compute an approximate solution of a linear system. ! to compute an approximate solution of a linear system.
! !
! !
call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.) call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.) call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.)
!
! 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_) select case (init_)
case('Z') case('Z')
call ty%zero()
call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,info,init='Z')
case('Y') case('Y')
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
call psb_geaxpby(sone,y,szero,ty,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') case('U')
if (.not.present(initu)) then if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply') & a_err='missing initu to smoother_apply')
goto 9999 goto 9999
end if end if
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
call psb_geaxpby(sone,initu,szero,ty,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')
case default case default
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply') & a_err='wrong init to smoother_apply')
goto 9999 goto 9999
end select end select
do i=1, sweeps do i=1, sweeps-1
! !
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the ! 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 ! block diagonal part and the remaining part of the local matrix
@ -188,22 +191,17 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
goto 9999 goto 9999
end if end if
!!$ else
!!$
!!$ info = psb_err_iarg_neg_
!!$ call psb_errpush(info,name,&
!!$ & i_err=(/itwo,sweeps,izero,izero,izero/))
!!$ goto 9999
!!$
!!$ endif
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
deallocate(aux)
endif
else else
deallocate(ww,aux)
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
endif
if (.not.(4*n_col <= size(work))) then
deallocate(aux)
endif endif
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)

88
mlprec/impl/smoother/mld_z_jac_smoother_apply.f90
Unescape Escape View File

@ -55,7 +55,7 @@ subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
integer(psb_ipk_) :: n_row,n_col integer(psb_ipk_) :: n_row,n_col
complex(psb_dpk_), allocatable :: tx(:),ty(:) complex(psb_dpk_), allocatable :: tx(:),ty(:)
complex(psb_dpk_), pointer :: ww(:), aux(:) complex(psb_dpk_), pointer :: aux(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act integer(psb_ipk_) :: ictxt,np,me,i, err_act
character :: trans_, init_ character :: trans_, init_
character(len=20) :: name='z_jac_smoother_apply' character(len=20) :: name='z_jac_smoother_apply'
@ -89,71 +89,73 @@ subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
n_row = desc_data%get_local_rows() n_row = desc_data%get_local_rows()
n_col = desc_data%get_local_cols() n_col = desc_data%get_local_cols()
if (n_col <= size(work)) then if (4*n_col <= size(work)) then
ww => work(1:n_col) aux => work(:)
if ((4*n_col+n_col) <= size(work)) then
aux => work(n_col+1:)
else
allocate(aux(4*n_col),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,&
& i_err=(/4*n_col,izero,izero,izero,izero/),&
& a_err='complex(psb_dpk_)')
goto 9999
end if
endif
else else
allocate(ww(n_col),aux(4*n_col),stat=info) allocate(aux(4*n_col),stat=info)
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_alloc_request_ info=psb_err_alloc_request_
call psb_errpush(info,name,& call psb_errpush(info,name,&
& i_err=(/5*n_col,izero,izero,izero,izero/),& & i_err=(/4*n_col,izero,izero,izero,izero/),&
& a_err='complex(psb_dpk_)') & a_err='complex(psb_dpk_)')
goto 9999 goto 9999
end if end if
endif endif
!!$ if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then
!!$ ! if .not.sv%is_iterative, there's no need to pass init ! 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) call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
!!$
!!$ if (info /= psb_success_) then if (info /= psb_success_) then
!!$ call psb_errpush(psb_err_internal_error_,& call psb_errpush(psb_err_internal_error_,&
!!$ & name,a_err='Error in sub_aply Jacobi Sweeps = 1') & name,a_err='Error in sub_aply Jacobi Sweeps = 1')
!!$ goto 9999 goto 9999
!!$ endif endif
!!$
!!$ else if (sweeps >= 1) then else if (sweeps >= 0) then
! !
! !
! Apply multiple sweeps of a block-Jacobi solver ! Apply multiple sweeps of a block-Jacobi solver
! to compute an approximate solution of a linear system. ! to compute an approximate solution of a linear system.
! !
! !
call psb_geasb(tx,desc_data,info) call psb_geasb(tx,desc_data,info)
call psb_geasb(ty,desc_data,info) call psb_geasb(ty,desc_data,info)
!
! 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_) select case (init_)
case('Z') case('Z')
ty(:) = zzero
call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,info,init='Z')
case('Y') case('Y')
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
call psb_geaxpby(zone,y,zzero,ty,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') case('U')
if (.not.present(initu)) then if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply') & a_err='missing initu to smoother_apply')
goto 9999 goto 9999
end if end if
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
call psb_geaxpby(zone,initu,zzero,ty,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')
case default case default
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply') & a_err='wrong init to smoother_apply')
goto 9999 goto 9999
end select end select
do i=1, sweeps do i=1, sweeps-1
! !
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the ! 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 ! block diagonal part and the remaining part of the local matrix
@ -186,23 +188,17 @@ subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
goto 9999 goto 9999
end if end if
!!$ else else
!!$
!!$ info = psb_err_iarg_neg_
!!$ call psb_errpush(info,name,&
!!$ & i_err=(/itwo,sweeps,izero,izero,izero/))
!!$ goto 9999
!!$
!!$ endif
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
if (n_col <= size(work)) then endif
if ((4*n_col+n_col) <= size(work)) then
else if (.not.(4*n_col <= size(work))) then
deallocate(aux) deallocate(aux)
endif
else
deallocate(ww,aux)
endif endif
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)

92
mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90
Unescape Escape View File

@ -56,10 +56,10 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
integer(psb_ipk_) :: n_row,n_col integer(psb_ipk_) :: n_row,n_col
type(psb_z_vect_type) :: tx, ty type(psb_z_vect_type) :: tx, ty
complex(psb_dpk_), pointer :: ww(:), aux(:) complex(psb_dpk_), pointer :: aux(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act integer(psb_ipk_) :: ictxt,np,me,i, err_act
character :: trans_, init_ character :: trans_, init_
character(len=20) :: name='z_jac_smoother_apply' character(len=20) :: name='z_jac_smoother_apply_v'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -90,71 +90,74 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
n_row = desc_data%get_local_rows() n_row = desc_data%get_local_rows()
n_col = desc_data%get_local_cols() n_col = desc_data%get_local_cols()
if (n_col <= size(work)) then if (4*n_col <= size(work)) then
ww => work(1:n_col) aux => work(:)
if ((4*n_col+n_col) <= size(work)) then
aux => work(n_col+1:)
else
allocate(aux(4*n_col),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,&
& i_err=(/4*n_col,izero,izero,izero,izero/),&
& a_err='complex(psb_dpk_)')
goto 9999
end if
endif
else else
allocate(ww(n_col),aux(4*n_col),stat=info) allocate(aux(4*n_col),stat=info)
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_alloc_request_ info=psb_err_alloc_request_
call psb_errpush(info,name,& call psb_errpush(info,name,&
& i_err=(/5*n_col,izero,izero,izero,izero/),& & i_err=(/4*n_col,izero,izero,izero,izero/),&
& a_err='complex(psb_dpk_)') & a_err='complex(psb_dpk_)')
goto 9999 goto 9999
end if end if
endif endif
!!$ if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then
!!$ ! if .not.sv%is_iterative, there's no need to pass init ! 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) call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
!!$
!!$ if (info /= psb_success_) then if (info /= psb_success_) then
!!$ call psb_errpush(psb_err_internal_error_,& call psb_errpush(psb_err_internal_error_,&
!!$ & name,a_err='Error in sub_aply Jacobi Sweeps = 1') & name,a_err='Error in sub_aply Jacobi Sweeps = 1')
!!$ goto 9999 goto 9999
!!$ endif endif
!!$
!!$ else if (sweeps >= 0) then else if (sweeps >= 0) then
! !
! !
! Apply multiple sweeps of a block-Jacobi solver ! Apply multiple sweeps of a block-Jacobi solver
! to compute an approximate solution of a linear system. ! to compute an approximate solution of a linear system.
! !
! !
call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.) call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.) call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.)
!
! 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_) select case (init_)
case('Z') case('Z')
call ty%zero()
call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,info,init='Z')
case('Y') case('Y')
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
call psb_geaxpby(zone,y,zzero,ty,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') case('U')
if (.not.present(initu)) then if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply') & a_err='missing initu to smoother_apply')
goto 9999 goto 9999
end if end if
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
call psb_geaxpby(zone,initu,zzero,ty,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')
case default case default
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply') & a_err='wrong init to smoother_apply')
goto 9999 goto 9999
end select end select
do i=1, sweeps do i=1, sweeps-1
! !
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the ! 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 ! block diagonal part and the remaining part of the local matrix
@ -188,22 +191,17 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
goto 9999 goto 9999
end if end if
!!$ else
!!$
!!$ info = psb_err_iarg_neg_
!!$ call psb_errpush(info,name,&
!!$ & i_err=(/itwo,sweeps,izero,izero,izero/))
!!$ goto 9999
!!$
!!$ endif
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
deallocate(aux)
endif
else else
deallocate(ww,aux)
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
endif
if (.not.(4*n_col <= size(work))) then
deallocate(aux)
endif endif
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)

Write Preview
Loading…
Cancel
Save