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

New Jacobi implementation.

Browse Source
stopcriterion
Salvatore Filippone 6 years ago
parent 19d33fae2c
commit b2047f95e0
24 changed files with 1076 additions and 544 deletions
Show Stats Download Patch File Download Diff File

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

@ -114,6 +114,66 @@ subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
endif endif
else if (sweeps >= 0) then else if (sweeps >= 0) then
if (associated(sm%pa)) then
!
! This means we are dealing with a pure Jacobi smoother/solver.
!
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 ! Apply multiple sweeps of a block-Jacobi solver
@ -180,6 +240,7 @@ subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
& a_err='subsolve with Jacobi sweeps > 1') & a_err='subsolve with Jacobi sweeps > 1')
goto 9999 goto 9999
end if end if
end if
deallocate(tx,ty,stat=info) deallocate(tx,ty,stat=info)
if (info /= psb_success_) then if (info /= psb_success_) then

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

@ -116,6 +116,68 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
endif endif
else if (sweeps >= 0) then else if (sweeps >= 0) then
if (associated(sm%pa)) then
!
! This means we are dealing with a pure Jacobi smoother/solver.
!
associate(tx => wv(1), ty => wv(2))
select case (init_)
case('Z')
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%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='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
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
else
! !
! !
! Apply multiple sweeps of a block-Jacobi solver ! Apply multiple sweeps of a block-Jacobi solver
@ -189,6 +251,7 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
end if end if
end associate end associate
end if
else else

16
mlprec/impl/smoother/mld_c_jac_smoother_bld.f90
Unescape Escape View File

@ -71,12 +71,14 @@ subroutine mld_c_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
nrow_a = a%get_nrows() nrow_a = a%get_nrows()
nztota = a%get_nzeros() nztota = a%get_nzeros()
select type (smsv => sm%sv) select type (smsv => sm%sv)
type is (mld_c_diag_solver_type) class is (mld_c_diag_solver_type)
call a%clip_diag(sm%nd,info) call sm%nd%free()
sm%pa => a
sm%nnz_nd_tot = nztota
class default class default
call a%csclip(sm%nd,info,& call a%csclip(sm%nd,info,&
& jmin=nrow_a+1,rscale=.false.,cscale=.false.) & jmin=nrow_a+1,rscale=.false.,cscale=.false.)
end select
if (info == psb_success_) then if (info == psb_success_) then
if (present(amold)) then if (present(amold)) then
call sm%nd%cscnv(info,& call sm%nd%cscnv(info,&
@ -86,11 +88,16 @@ subroutine mld_c_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
& type='csr',dupl=psb_dupl_add_) & type='csr',dupl=psb_dupl_add_)
endif endif
end if end if
sm%nnz_nd_tot = sm%nd%get_nzeros()
end select
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,& call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='clip & psb_spcnv csr 4') & a_err='clip & psb_spcnv csr 4')
goto 9999 goto 9999
end if end if
call psb_sum(ictxt,sm%nnz_nd_tot)
call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold) call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -98,9 +105,6 @@ subroutine mld_c_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
& a_err='solver build') & a_err='solver build')
goto 9999 goto 9999
end if end if
nzeros = sm%nd%get_nzeros()
call psb_sum(ictxt,nzeros)
sm%nnz_nd_tot = nzeros
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' end' & write(debug_unit,*) me,' ',trim(name),' end'

2
mlprec/impl/smoother/mld_c_jac_smoother_cnv.f90
Unescape Escape View File

@ -59,6 +59,7 @@ subroutine mld_c_jac_smoother_cnv(sm,info,amold,vmold,imold)
if (info == psb_success_) then if (info == psb_success_) then
if (.not.associated(sm%pa)) then
if (present(amold)) then if (present(amold)) then
call sm%nd%cscnv(info,& call sm%nd%cscnv(info,&
& mold=amold,dupl=psb_dupl_add_) & mold=amold,dupl=psb_dupl_add_)
@ -67,6 +68,7 @@ subroutine mld_c_jac_smoother_cnv(sm,info,amold,vmold,imold)
& type='csr',dupl=psb_dupl_add_) & type='csr',dupl=psb_dupl_add_)
endif endif
end if end if
end if
if (allocated(sm%sv)) & if (allocated(sm%sv)) &
& call sm%sv%cnv(info,amold=amold,vmold=vmold,imold=imold) & call sm%sv%cnv(info,amold=amold,vmold=vmold,imold=imold)

1
mlprec/impl/smoother/mld_c_jac_smoother_descr.f90
Unescape Escape View File

@ -75,6 +75,7 @@ subroutine mld_c_jac_smoother_descr(sm,info,iout,coarse)
select type(smv=>sm%sv) select type(smv=>sm%sv)
class is (mld_c_diag_solver_type) class is (mld_c_diag_solver_type)
write(iout_,*) ' Point Jacobi ' write(iout_,*) ' Point Jacobi '
call smv%descr(info,iout_,coarse=coarse)
class is (mld_c_bwgs_solver_type) class is (mld_c_bwgs_solver_type)
write(iout_,*) ' Hybrid Backward Gauss-Seidel ' write(iout_,*) ' Hybrid Backward Gauss-Seidel '
class is (mld_c_gs_solver_type) class is (mld_c_gs_solver_type)

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

@ -114,6 +114,66 @@ subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
endif endif
else if (sweeps >= 0) then else if (sweeps >= 0) then
if (associated(sm%pa)) then
!
! This means we are dealing with a pure Jacobi smoother/solver.
!
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 ! Apply multiple sweeps of a block-Jacobi solver
@ -180,6 +240,7 @@ subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
& a_err='subsolve with Jacobi sweeps > 1') & a_err='subsolve with Jacobi sweeps > 1')
goto 9999 goto 9999
end if end if
end if
deallocate(tx,ty,stat=info) deallocate(tx,ty,stat=info)
if (info /= psb_success_) then if (info /= psb_success_) then

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

@ -116,6 +116,68 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
endif endif
else if (sweeps >= 0) then else if (sweeps >= 0) then
if (associated(sm%pa)) then
!
! This means we are dealing with a pure Jacobi smoother/solver.
!
associate(tx => wv(1), ty => wv(2))
select case (init_)
case('Z')
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%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='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
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
else
! !
! !
! Apply multiple sweeps of a block-Jacobi solver ! Apply multiple sweeps of a block-Jacobi solver
@ -189,6 +251,7 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
end if end if
end associate end associate
end if
else else

16
mlprec/impl/smoother/mld_d_jac_smoother_bld.f90
Unescape Escape View File

@ -71,12 +71,14 @@ subroutine mld_d_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
nrow_a = a%get_nrows() nrow_a = a%get_nrows()
nztota = a%get_nzeros() nztota = a%get_nzeros()
select type (smsv => sm%sv) select type (smsv => sm%sv)
type is (mld_d_diag_solver_type) class is (mld_d_diag_solver_type)
call a%clip_diag(sm%nd,info) call sm%nd%free()
sm%pa => a
sm%nnz_nd_tot = nztota
class default class default
call a%csclip(sm%nd,info,& call a%csclip(sm%nd,info,&
& jmin=nrow_a+1,rscale=.false.,cscale=.false.) & jmin=nrow_a+1,rscale=.false.,cscale=.false.)
end select
if (info == psb_success_) then if (info == psb_success_) then
if (present(amold)) then if (present(amold)) then
call sm%nd%cscnv(info,& call sm%nd%cscnv(info,&
@ -86,11 +88,16 @@ subroutine mld_d_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
& type='csr',dupl=psb_dupl_add_) & type='csr',dupl=psb_dupl_add_)
endif endif
end if end if
sm%nnz_nd_tot = sm%nd%get_nzeros()
end select
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,& call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='clip & psb_spcnv csr 4') & a_err='clip & psb_spcnv csr 4')
goto 9999 goto 9999
end if end if
call psb_sum(ictxt,sm%nnz_nd_tot)
call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold) call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -98,9 +105,6 @@ subroutine mld_d_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
& a_err='solver build') & a_err='solver build')
goto 9999 goto 9999
end if end if
nzeros = sm%nd%get_nzeros()
call psb_sum(ictxt,nzeros)
sm%nnz_nd_tot = nzeros
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' end' & write(debug_unit,*) me,' ',trim(name),' end'

2
mlprec/impl/smoother/mld_d_jac_smoother_cnv.f90
Unescape Escape View File

@ -59,6 +59,7 @@ subroutine mld_d_jac_smoother_cnv(sm,info,amold,vmold,imold)
if (info == psb_success_) then if (info == psb_success_) then
if (.not.associated(sm%pa)) then
if (present(amold)) then if (present(amold)) then
call sm%nd%cscnv(info,& call sm%nd%cscnv(info,&
& mold=amold,dupl=psb_dupl_add_) & mold=amold,dupl=psb_dupl_add_)
@ -67,6 +68,7 @@ subroutine mld_d_jac_smoother_cnv(sm,info,amold,vmold,imold)
& type='csr',dupl=psb_dupl_add_) & type='csr',dupl=psb_dupl_add_)
endif endif
end if end if
end if
if (allocated(sm%sv)) & if (allocated(sm%sv)) &
& call sm%sv%cnv(info,amold=amold,vmold=vmold,imold=imold) & call sm%sv%cnv(info,amold=amold,vmold=vmold,imold=imold)

1
mlprec/impl/smoother/mld_d_jac_smoother_descr.f90
Unescape Escape View File

@ -75,6 +75,7 @@ subroutine mld_d_jac_smoother_descr(sm,info,iout,coarse)
select type(smv=>sm%sv) select type(smv=>sm%sv)
class is (mld_d_diag_solver_type) class is (mld_d_diag_solver_type)
write(iout_,*) ' Point Jacobi ' write(iout_,*) ' Point Jacobi '
call smv%descr(info,iout_,coarse=coarse)
class is (mld_d_bwgs_solver_type) class is (mld_d_bwgs_solver_type)
write(iout_,*) ' Hybrid Backward Gauss-Seidel ' write(iout_,*) ' Hybrid Backward Gauss-Seidel '
class is (mld_d_gs_solver_type) class is (mld_d_gs_solver_type)

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

@ -114,6 +114,66 @@ subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
endif endif
else if (sweeps >= 0) then else if (sweeps >= 0) then
if (associated(sm%pa)) then
!
! This means we are dealing with a pure Jacobi smoother/solver.
!
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 ! Apply multiple sweeps of a block-Jacobi solver
@ -180,6 +240,7 @@ subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
& a_err='subsolve with Jacobi sweeps > 1') & a_err='subsolve with Jacobi sweeps > 1')
goto 9999 goto 9999
end if end if
end if
deallocate(tx,ty,stat=info) deallocate(tx,ty,stat=info)
if (info /= psb_success_) then if (info /= psb_success_) then

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

@ -116,6 +116,68 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
endif endif
else if (sweeps >= 0) then else if (sweeps >= 0) then
if (associated(sm%pa)) then
!
! This means we are dealing with a pure Jacobi smoother/solver.
!
associate(tx => wv(1), ty => wv(2))
select case (init_)
case('Z')
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%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='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
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
else
! !
! !
! Apply multiple sweeps of a block-Jacobi solver ! Apply multiple sweeps of a block-Jacobi solver
@ -189,6 +251,7 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
end if end if
end associate end associate
end if
else else

16
mlprec/impl/smoother/mld_s_jac_smoother_bld.f90
Unescape Escape View File

@ -71,12 +71,14 @@ subroutine mld_s_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
nrow_a = a%get_nrows() nrow_a = a%get_nrows()
nztota = a%get_nzeros() nztota = a%get_nzeros()
select type (smsv => sm%sv) select type (smsv => sm%sv)
type is (mld_s_diag_solver_type) class is (mld_s_diag_solver_type)
call a%clip_diag(sm%nd,info) call sm%nd%free()
sm%pa => a
sm%nnz_nd_tot = nztota
class default class default
call a%csclip(sm%nd,info,& call a%csclip(sm%nd,info,&
& jmin=nrow_a+1,rscale=.false.,cscale=.false.) & jmin=nrow_a+1,rscale=.false.,cscale=.false.)
end select
if (info == psb_success_) then if (info == psb_success_) then
if (present(amold)) then if (present(amold)) then
call sm%nd%cscnv(info,& call sm%nd%cscnv(info,&
@ -86,11 +88,16 @@ subroutine mld_s_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
& type='csr',dupl=psb_dupl_add_) & type='csr',dupl=psb_dupl_add_)
endif endif
end if end if
sm%nnz_nd_tot = sm%nd%get_nzeros()
end select
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,& call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='clip & psb_spcnv csr 4') & a_err='clip & psb_spcnv csr 4')
goto 9999 goto 9999
end if end if
call psb_sum(ictxt,sm%nnz_nd_tot)
call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold) call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -98,9 +105,6 @@ subroutine mld_s_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
& a_err='solver build') & a_err='solver build')
goto 9999 goto 9999
end if end if
nzeros = sm%nd%get_nzeros()
call psb_sum(ictxt,nzeros)
sm%nnz_nd_tot = nzeros
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' end' & write(debug_unit,*) me,' ',trim(name),' end'

2
mlprec/impl/smoother/mld_s_jac_smoother_cnv.f90
Unescape Escape View File

@ -59,6 +59,7 @@ subroutine mld_s_jac_smoother_cnv(sm,info,amold,vmold,imold)
if (info == psb_success_) then if (info == psb_success_) then
if (.not.associated(sm%pa)) then
if (present(amold)) then if (present(amold)) then
call sm%nd%cscnv(info,& call sm%nd%cscnv(info,&
& mold=amold,dupl=psb_dupl_add_) & mold=amold,dupl=psb_dupl_add_)
@ -67,6 +68,7 @@ subroutine mld_s_jac_smoother_cnv(sm,info,amold,vmold,imold)
& type='csr',dupl=psb_dupl_add_) & type='csr',dupl=psb_dupl_add_)
endif endif
end if end if
end if
if (allocated(sm%sv)) & if (allocated(sm%sv)) &
& call sm%sv%cnv(info,amold=amold,vmold=vmold,imold=imold) & call sm%sv%cnv(info,amold=amold,vmold=vmold,imold=imold)

1
mlprec/impl/smoother/mld_s_jac_smoother_descr.f90
Unescape Escape View File

@ -75,6 +75,7 @@ subroutine mld_s_jac_smoother_descr(sm,info,iout,coarse)
select type(smv=>sm%sv) select type(smv=>sm%sv)
class is (mld_s_diag_solver_type) class is (mld_s_diag_solver_type)
write(iout_,*) ' Point Jacobi ' write(iout_,*) ' Point Jacobi '
call smv%descr(info,iout_,coarse=coarse)
class is (mld_s_bwgs_solver_type) class is (mld_s_bwgs_solver_type)
write(iout_,*) ' Hybrid Backward Gauss-Seidel ' write(iout_,*) ' Hybrid Backward Gauss-Seidel '
class is (mld_s_gs_solver_type) class is (mld_s_gs_solver_type)

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

@ -114,6 +114,66 @@ subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
endif endif
else if (sweeps >= 0) then else if (sweeps >= 0) then
if (associated(sm%pa)) then
!
! This means we are dealing with a pure Jacobi smoother/solver.
!
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 ! Apply multiple sweeps of a block-Jacobi solver
@ -180,6 +240,7 @@ subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
& a_err='subsolve with Jacobi sweeps > 1') & a_err='subsolve with Jacobi sweeps > 1')
goto 9999 goto 9999
end if end if
end if
deallocate(tx,ty,stat=info) deallocate(tx,ty,stat=info)
if (info /= psb_success_) then if (info /= psb_success_) then

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

@ -116,6 +116,68 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
endif endif
else if (sweeps >= 0) then else if (sweeps >= 0) then
if (associated(sm%pa)) then
!
! This means we are dealing with a pure Jacobi smoother/solver.
!
associate(tx => wv(1), ty => wv(2))
select case (init_)
case('Z')
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%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='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
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
else
! !
! !
! Apply multiple sweeps of a block-Jacobi solver ! Apply multiple sweeps of a block-Jacobi solver
@ -189,6 +251,7 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
end if end if
end associate end associate
end if
else else

16
mlprec/impl/smoother/mld_z_jac_smoother_bld.f90
Unescape Escape View File

@ -71,12 +71,14 @@ subroutine mld_z_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
nrow_a = a%get_nrows() nrow_a = a%get_nrows()
nztota = a%get_nzeros() nztota = a%get_nzeros()
select type (smsv => sm%sv) select type (smsv => sm%sv)
type is (mld_z_diag_solver_type) class is (mld_z_diag_solver_type)
call a%clip_diag(sm%nd,info) call sm%nd%free()
sm%pa => a
sm%nnz_nd_tot = nztota
class default class default
call a%csclip(sm%nd,info,& call a%csclip(sm%nd,info,&
& jmin=nrow_a+1,rscale=.false.,cscale=.false.) & jmin=nrow_a+1,rscale=.false.,cscale=.false.)
end select
if (info == psb_success_) then if (info == psb_success_) then
if (present(amold)) then if (present(amold)) then
call sm%nd%cscnv(info,& call sm%nd%cscnv(info,&
@ -86,11 +88,16 @@ subroutine mld_z_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
& type='csr',dupl=psb_dupl_add_) & type='csr',dupl=psb_dupl_add_)
endif endif
end if end if
sm%nnz_nd_tot = sm%nd%get_nzeros()
end select
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,& call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='clip & psb_spcnv csr 4') & a_err='clip & psb_spcnv csr 4')
goto 9999 goto 9999
end if end if
call psb_sum(ictxt,sm%nnz_nd_tot)
call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold) call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -98,9 +105,6 @@ subroutine mld_z_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
& a_err='solver build') & a_err='solver build')
goto 9999 goto 9999
end if end if
nzeros = sm%nd%get_nzeros()
call psb_sum(ictxt,nzeros)
sm%nnz_nd_tot = nzeros
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' end' & write(debug_unit,*) me,' ',trim(name),' end'

2
mlprec/impl/smoother/mld_z_jac_smoother_cnv.f90
Unescape Escape View File

@ -59,6 +59,7 @@ subroutine mld_z_jac_smoother_cnv(sm,info,amold,vmold,imold)
if (info == psb_success_) then if (info == psb_success_) then
if (.not.associated(sm%pa)) then
if (present(amold)) then if (present(amold)) then
call sm%nd%cscnv(info,& call sm%nd%cscnv(info,&
& mold=amold,dupl=psb_dupl_add_) & mold=amold,dupl=psb_dupl_add_)
@ -67,6 +68,7 @@ subroutine mld_z_jac_smoother_cnv(sm,info,amold,vmold,imold)
& type='csr',dupl=psb_dupl_add_) & type='csr',dupl=psb_dupl_add_)
endif endif
end if end if
end if
if (allocated(sm%sv)) & if (allocated(sm%sv)) &
& call sm%sv%cnv(info,amold=amold,vmold=vmold,imold=imold) & call sm%sv%cnv(info,amold=amold,vmold=vmold,imold=imold)

1
mlprec/impl/smoother/mld_z_jac_smoother_descr.f90
Unescape Escape View File

@ -75,6 +75,7 @@ subroutine mld_z_jac_smoother_descr(sm,info,iout,coarse)
select type(smv=>sm%sv) select type(smv=>sm%sv)
class is (mld_z_diag_solver_type) class is (mld_z_diag_solver_type)
write(iout_,*) ' Point Jacobi ' write(iout_,*) ' Point Jacobi '
call smv%descr(info,iout_,coarse=coarse)
class is (mld_z_bwgs_solver_type) class is (mld_z_bwgs_solver_type)
write(iout_,*) ' Hybrid Backward Gauss-Seidel ' write(iout_,*) ' Hybrid Backward Gauss-Seidel '
class is (mld_z_gs_solver_type) class is (mld_z_gs_solver_type)

2
mlprec/mld_c_jac_smoother.f90
Unescape Escape View File

@ -58,6 +58,7 @@ module mld_c_jac_smoother
! parent type. ! parent type.
! class(mld_c_base_solver_type), allocatable :: sv ! class(mld_c_base_solver_type), allocatable :: sv
! !
type(psb_cspmat_type), pointer :: pa => null()
type(psb_cspmat_type) :: nd type(psb_cspmat_type) :: nd
integer(psb_ipk_) :: nnz_nd_tot integer(psb_ipk_) :: nnz_nd_tot
contains contains
@ -217,6 +218,7 @@ contains
end if end if
end if end if
call sm%nd%free() call sm%nd%free()
sm%pa => null()
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

2
mlprec/mld_d_jac_smoother.f90
Unescape Escape View File

@ -58,6 +58,7 @@ module mld_d_jac_smoother
! parent type. ! parent type.
! class(mld_d_base_solver_type), allocatable :: sv ! class(mld_d_base_solver_type), allocatable :: sv
! !
type(psb_dspmat_type), pointer :: pa => null()
type(psb_dspmat_type) :: nd type(psb_dspmat_type) :: nd
integer(psb_ipk_) :: nnz_nd_tot integer(psb_ipk_) :: nnz_nd_tot
contains contains
@ -217,6 +218,7 @@ contains
end if end if
end if end if
call sm%nd%free() call sm%nd%free()
sm%pa => null()
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

2
mlprec/mld_s_jac_smoother.f90
Unescape Escape View File

@ -58,6 +58,7 @@ module mld_s_jac_smoother
! parent type. ! parent type.
! class(mld_s_base_solver_type), allocatable :: sv ! class(mld_s_base_solver_type), allocatable :: sv
! !
type(psb_sspmat_type), pointer :: pa => null()
type(psb_sspmat_type) :: nd type(psb_sspmat_type) :: nd
integer(psb_ipk_) :: nnz_nd_tot integer(psb_ipk_) :: nnz_nd_tot
contains contains
@ -217,6 +218,7 @@ contains
end if end if
end if end if
call sm%nd%free() call sm%nd%free()
sm%pa => null()
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

2
mlprec/mld_z_jac_smoother.f90
Unescape Escape View File

@ -58,6 +58,7 @@ module mld_z_jac_smoother
! parent type. ! parent type.
! class(mld_z_base_solver_type), allocatable :: sv ! class(mld_z_base_solver_type), allocatable :: sv
! !
type(psb_zspmat_type), pointer :: pa => null()
type(psb_zspmat_type) :: nd type(psb_zspmat_type) :: nd
integer(psb_ipk_) :: nnz_nd_tot integer(psb_ipk_) :: nnz_nd_tot
contains contains
@ -217,6 +218,7 @@ contains
end if end if
end if end if
call sm%nd%free() call sm%nd%free()
sm%pa => null()
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

Write Preview
Loading…
Cancel
Save