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

Merge branch 'richardson' into unify_aggr_bld_richardson

Browse Source
unify_aggr_bld_richardson
Cirdans-Home 5 years ago
parent 43fd5268a1 be4e7b4837
commit 6bfd227a80
44 changed files with 2216 additions and 1992 deletions
Show Stats Download Patch File Download Diff File

292
mlprec/impl/smoother/mld_c_as_smoother_apply.f90
Unescape Escape View File

@ -105,116 +105,218 @@ subroutine mld_c_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
end if
endif
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
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in sub_aply Jacobi Sweeps = 1')
goto 9999
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)
!
! 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')
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
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='missing initu to smoother_apply')
& 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
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
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')
!
! 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/))
goto 9999
endif
do i=1, sweeps-1
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
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in sub_aply Jacobi Sweeps = 1')
goto 9999
endif
else if (sweeps >= 0) then
!
!
! Apply multiple sweeps of an AS solver
! to compute an approximate solution of a linear system.
!
! 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_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
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Y')
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
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
end do
if (info /= psb_success_) then
info=psb_err_internal_error_
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_
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_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
endif
endif
end if
if (.not.(4*isz <= size(work))) then
deallocate(aux,stat=info)
endif

335
mlprec/impl/smoother/mld_c_as_smoother_apply_vect.f90
Unescape Escape View File

@ -107,89 +107,140 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
end if
endif
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)
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 (sweeps >= 0) then
!
!
! Apply multiple sweeps of an AS solver
! to compute an approximate solution of a linear system.
!
!
if (size(wv) < 3) then
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
!
! 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()
goto 9999
end if
!
! 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)
! 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.
!
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,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
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='missing initu to smoother_apply')
& 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
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)
!
! 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/))
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,wv,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -197,47 +248,125 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
goto 9999
endif
do i=1, sweeps-1
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()
!
! 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.
! 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,wv(4:),info,init='Z')
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
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')
if (info /= psb_success_) exit
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)
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
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
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_)
!
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
!
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
end associate
if (info /= psb_success_) exit
else
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
endif
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
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)

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

@ -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_

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

@ -77,6 +77,7 @@ 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_

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

@ -61,6 +61,7 @@ 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

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

@ -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,146 +102,63 @@ subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
goto 9999
end if
endif
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
endif
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
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='wrong init to smoother_apply')
& a_err='missing initu 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)
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
!
! 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)
! Compute Y(j+1) = Y(j)+ M^(-1)*(X-A*Y(j)),
!
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%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')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
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_)
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
if (info /= psb_success_) exit
call sm%sv%apply(cone,tx,cone,ty,desc_data,trans_,aux,info,init='Y')
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_) exit
end do
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
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
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
deallocate(tx,ty,stat=info)
if (info /= psb_success_) then
info=psb_err_internal_error_
@ -250,6 +167,12 @@ 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_

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

@ -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,190 +108,90 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
end if
endif
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 (sweeps > 0) then
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)
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 (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 (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')
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')
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 default
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
& a_err='missing initu 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 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
end do
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
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
call sm%sv%apply(cone,tx,cone,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
if (info /= psb_success_) then
info=psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='subsolve with Jacobi sweeps > 1')
goto 9999
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
end associate
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_
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='invalid wv size in smoother_apply')
& a_err='subsolve with Jacobi sweeps > 1')
goto 9999
end if
associate(tx => wv(1), ty => wv(2))
!
! 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,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 associate
end select
else if (sweeps == 0) then
!
! 0 sweeps of smoother is the identity operator
!
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
else

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

@ -71,12 +71,11 @@ subroutine mld_c_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
nrow_a = a%get_nrows()
nztota = a%get_nzeros()
if( sm%checkres ) sm%pa => a
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)

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

@ -77,6 +77,7 @@ 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_

3
mlprec/impl/smoother/mld_c_l1_jac_smoother_bld.f90
Unescape Escape View File

@ -72,12 +72,11 @@ subroutine mld_c_l1_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
nrow_a = a%get_nrows()
nztota = a%get_nzeros()
if( sm%checkres ) sm%pa => a
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)

292
mlprec/impl/smoother/mld_d_as_smoother_apply.f90
Unescape Escape View File

@ -105,116 +105,218 @@ subroutine mld_d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
end if
endif
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
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in sub_aply Jacobi Sweeps = 1')
goto 9999
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)
!
! 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')
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
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='missing initu to smoother_apply')
& 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
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
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')
!
! 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/))
goto 9999
endif
do i=1, sweeps-1
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
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in sub_aply Jacobi Sweeps = 1')
goto 9999
endif
else if (sweeps >= 0) then
!
!
! Apply multiple sweeps of an AS solver
! to compute an approximate solution of a linear system.
!
! 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_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
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Y')
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
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
end do
if (info /= psb_success_) then
info=psb_err_internal_error_
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_
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_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
endif
endif
end if
if (.not.(4*isz <= size(work))) then
deallocate(aux,stat=info)
endif

335
mlprec/impl/smoother/mld_d_as_smoother_apply_vect.f90
Unescape Escape View File

@ -107,89 +107,140 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
end if
endif
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)
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 (sweeps >= 0) then
!
!
! Apply multiple sweeps of an AS solver
! to compute an approximate solution of a linear system.
!
!
if (size(wv) < 3) then
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
!
! 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()
goto 9999
end if
!
! 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)
! 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.
!
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,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
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='missing initu to smoother_apply')
& 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
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)
!
! 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/))
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,wv,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -197,47 +248,125 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
goto 9999
endif
do i=1, sweeps-1
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()
!
! 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.
! 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,wv(4:),info,init='Z')
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
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')
if (info /= psb_success_) exit
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)
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
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
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_)
!
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
!
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
end associate
if (info /= psb_success_) exit
else
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
endif
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
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)

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

@ -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_

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

@ -77,6 +77,7 @@ 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_

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

@ -61,6 +61,7 @@ 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

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

@ -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,146 +102,63 @@ subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
goto 9999
end if
endif
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
endif
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
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='wrong init to smoother_apply')
& a_err='missing initu 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)
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
!
! 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)
! Compute Y(j+1) = Y(j)+ M^(-1)*(X-A*Y(j)),
!
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%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')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
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_)
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
if (info /= psb_success_) exit
call sm%sv%apply(done,tx,done,ty,desc_data,trans_,aux,info,init='Y')
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_) exit
end do
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
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
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
deallocate(tx,ty,stat=info)
if (info /= psb_success_) then
info=psb_err_internal_error_
@ -250,6 +167,12 @@ 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_

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

@ -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,190 +108,90 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
end if
endif
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 (sweeps > 0) then
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)
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 (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 (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')
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')
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 default
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
& a_err='missing initu 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 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
end do
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
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
call sm%sv%apply(done,tx,done,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
if (info /= psb_success_) then
info=psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='subsolve with Jacobi sweeps > 1')
goto 9999
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
end associate
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_
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='invalid wv size in smoother_apply')
& a_err='subsolve with Jacobi sweeps > 1')
goto 9999
end if
associate(tx => wv(1), ty => wv(2))
!
! 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,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 associate
end select
else if (sweeps == 0) then
!
! 0 sweeps of smoother is the identity operator
!
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
else

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

@ -71,12 +71,11 @@ subroutine mld_d_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
nrow_a = a%get_nrows()
nztota = a%get_nzeros()
if( sm%checkres ) sm%pa => a
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)

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

@ -77,6 +77,7 @@ 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_

3
mlprec/impl/smoother/mld_d_l1_jac_smoother_bld.f90
Unescape Escape View File

@ -72,12 +72,11 @@ subroutine mld_d_l1_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
nrow_a = a%get_nrows()
nztota = a%get_nzeros()
if( sm%checkres ) sm%pa => a
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)

292
mlprec/impl/smoother/mld_s_as_smoother_apply.f90
Unescape Escape View File

@ -105,116 +105,218 @@ subroutine mld_s_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
end if
endif
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
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in sub_aply Jacobi Sweeps = 1')
goto 9999
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)
!
! 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')
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
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='missing initu to smoother_apply')
& 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
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
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')
!
! 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/))
goto 9999
endif
do i=1, sweeps-1
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
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in sub_aply Jacobi Sweeps = 1')
goto 9999
endif
else if (sweeps >= 0) then
!
!
! Apply multiple sweeps of an AS solver
! to compute an approximate solution of a linear system.
!
! 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_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
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Y')
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
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
end do
if (info /= psb_success_) then
info=psb_err_internal_error_
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_
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_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
endif
endif
end if
if (.not.(4*isz <= size(work))) then
deallocate(aux,stat=info)
endif

335
mlprec/impl/smoother/mld_s_as_smoother_apply_vect.f90
Unescape Escape View File

@ -107,89 +107,140 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
end if
endif
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)
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 (sweeps >= 0) then
!
!
! Apply multiple sweeps of an AS solver
! to compute an approximate solution of a linear system.
!
!
if (size(wv) < 3) then
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
!
! 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()
goto 9999
end if
!
! 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)
! 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.
!
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,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
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='missing initu to smoother_apply')
& 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
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)
!
! 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/))
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,wv,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -197,47 +248,125 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
goto 9999
endif
do i=1, sweeps-1
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()
!
! 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.
! 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,wv(4:),info,init='Z')
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
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')
if (info /= psb_success_) exit
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)
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
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
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_)
!
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
!
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
end associate
if (info /= psb_success_) exit
else
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
endif
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
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)

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

@ -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_

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

@ -77,6 +77,7 @@ 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_

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

@ -61,6 +61,7 @@ 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

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

@ -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,146 +102,63 @@ subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
goto 9999
end if
endif
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
endif
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
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='wrong init to smoother_apply')
& a_err='missing initu 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)
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
!
! 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)
! Compute Y(j+1) = Y(j)+ M^(-1)*(X-A*Y(j)),
!
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%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')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
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_)
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
if (info /= psb_success_) exit
call sm%sv%apply(sone,tx,sone,ty,desc_data,trans_,aux,info,init='Y')
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_) exit
end do
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
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
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
deallocate(tx,ty,stat=info)
if (info /= psb_success_) then
info=psb_err_internal_error_
@ -250,6 +167,12 @@ 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_

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

@ -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,190 +108,90 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
end if
endif
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 (sweeps > 0) then
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)
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 (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 (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')
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')
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 default
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
& a_err='missing initu 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 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
end do
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
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
call sm%sv%apply(sone,tx,sone,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
if (info /= psb_success_) then
info=psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='subsolve with Jacobi sweeps > 1')
goto 9999
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
end associate
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_
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='invalid wv size in smoother_apply')
& a_err='subsolve with Jacobi sweeps > 1')
goto 9999
end if
associate(tx => wv(1), ty => wv(2))
!
! 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,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 associate
end select
else if (sweeps == 0) then
!
! 0 sweeps of smoother is the identity operator
!
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
else

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

@ -71,12 +71,11 @@ subroutine mld_s_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
nrow_a = a%get_nrows()
nztota = a%get_nzeros()
if( sm%checkres ) sm%pa => a
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)

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

@ -77,6 +77,7 @@ 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_

3
mlprec/impl/smoother/mld_s_l1_jac_smoother_bld.f90
Unescape Escape View File

@ -72,12 +72,11 @@ subroutine mld_s_l1_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
nrow_a = a%get_nrows()
nztota = a%get_nzeros()
if( sm%checkres ) sm%pa => a
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)

292
mlprec/impl/smoother/mld_z_as_smoother_apply.f90
Unescape Escape View File

@ -105,116 +105,218 @@ subroutine mld_z_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
end if
endif
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
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in sub_aply Jacobi Sweeps = 1')
goto 9999
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)
!
! 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')
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
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='missing initu to smoother_apply')
& 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
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
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')
!
! 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/))
goto 9999
endif
do i=1, sweeps-1
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
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in sub_aply Jacobi Sweeps = 1')
goto 9999
endif
else if (sweeps >= 0) then
!
!
! Apply multiple sweeps of an AS solver
! to compute an approximate solution of a linear system.
!
! 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_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
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Y')
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
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
end do
if (info /= psb_success_) then
info=psb_err_internal_error_
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_
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_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
endif
endif
end if
if (.not.(4*isz <= size(work))) then
deallocate(aux,stat=info)
endif

335
mlprec/impl/smoother/mld_z_as_smoother_apply_vect.f90
Unescape Escape View File

@ -107,89 +107,140 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
end if
endif
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)
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 (sweeps >= 0) then
!
!
! Apply multiple sweeps of an AS solver
! to compute an approximate solution of a linear system.
!
!
if (size(wv) < 3) then
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
!
! 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()
goto 9999
end if
!
! 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)
! 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.
!
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,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
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='missing initu to smoother_apply')
& 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
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)
!
! 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/))
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,wv,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -197,47 +248,125 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
goto 9999
endif
do i=1, sweeps-1
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()
!
! 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.
! 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,wv(4:),info,init='Z')
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
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')
if (info /= psb_success_) exit
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)
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
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
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_)
!
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
!
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
end associate
if (info /= psb_success_) exit
else
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y')
endif
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
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)

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

@ -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_

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

@ -77,6 +77,7 @@ 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_

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

@ -61,6 +61,7 @@ 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

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

@ -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,146 +102,63 @@ subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
goto 9999
end if
endif
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
endif
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
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='wrong init to smoother_apply')
& a_err='missing initu 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)
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
!
! 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)
! Compute Y(j+1) = Y(j)+ M^(-1)*(X-A*Y(j)),
!
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%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')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
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_)
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
if (info /= psb_success_) exit
call sm%sv%apply(zone,tx,zone,ty,desc_data,trans_,aux,info,init='Y')
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_) exit
end do
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
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
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
deallocate(tx,ty,stat=info)
if (info /= psb_success_) then
info=psb_err_internal_error_
@ -250,6 +167,12 @@ 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_

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

@ -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,190 +108,90 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
end if
endif
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 (sweeps > 0) then
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)
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 (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 (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')
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')
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 default
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
& a_err='missing initu 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 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
end do
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
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
call sm%sv%apply(zone,tx,zone,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
if (info /= psb_success_) then
info=psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='subsolve with Jacobi sweeps > 1')
goto 9999
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
end associate
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_
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='invalid wv size in smoother_apply')
& a_err='subsolve with Jacobi sweeps > 1')
goto 9999
end if
associate(tx => wv(1), ty => wv(2))
!
! 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,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 associate
end select
else if (sweeps == 0) then
!
! 0 sweeps of smoother is the identity operator
!
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
else

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

@ -71,12 +71,11 @@ subroutine mld_z_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
nrow_a = a%get_nrows()
nztota = a%get_nzeros()
if( sm%checkres ) sm%pa => a
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)

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

@ -77,6 +77,7 @@ 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_

3
mlprec/impl/smoother/mld_z_l1_jac_smoother_bld.f90
Unescape Escape View File

@ -72,12 +72,11 @@ subroutine mld_z_l1_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
nrow_a = a%get_nrows()
nztota = a%get_nzeros()
if( sm%checkres ) sm%pa => a
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)

1
mlprec/mld_c_as_smoother.f90
Unescape Escape View File

@ -65,6 +65,7 @@ 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

1
mlprec/mld_d_as_smoother.f90
Unescape Escape View File

@ -65,6 +65,7 @@ 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

1
mlprec/mld_s_as_smoother.f90
Unescape Escape View File

@ -65,6 +65,7 @@ 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

1
mlprec/mld_z_as_smoother.f90
Unescape Escape View File

@ -65,6 +65,7 @@ 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

Write Preview
Loading…
Cancel
Save