fdurastante/amg4psblas
1
0
Fork 0
mirror of https://github.com/sfilippone/amg4psblas.git synced 2026-10-06 22:55:12 +00:00
Code Issues Packages Projects Releases Wiki Activity

Added use of WV in solve apply_vect

Browse Source
This commit is contained in:
Salvatore Filippone
2017-12-12 11:27:18 +00:00
parent 63233716c4
commit 54d1478e21
12 changed files with 668 additions and 588 deletions
Download Patch File Download Diff File Expand all files Collapse all files
mlprec/impl/solver/mld_c_bwgs_solver_apply_vect.f90
+61 -55
View File
@@ -54,7 +54,6 @@ subroutine mld_c_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
type(psb_c_vect_type),intent(inout), optional :: initu
integer(psb_ipk_) :: n_row,n_col, itx
type(psb_c_vect_type) :: tw, xit
complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
complex(psb_spk_), allocatable :: temp(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act
@@ -121,69 +120,76 @@ subroutine mld_c_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
goto 9999
end if
call psb_geasb(tw,desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(xit,desc_data,info,mold=x%v,scratch=.true.)
select case (init_)
case('Z')
call xit%zero()
case('Y')
call psb_geaxpby(cone,y,czero,xit,desc_data,info)
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,xit,desc_data,info)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
if (size(wv) < 2) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='invalid wv size')
goto 9999
end select
select case(trans_)
case('N')
if (sv%eps <=szero) then
!
! Fixed number of iterations
!
!
do itx=1,sv%sweeps
call psb_geaxpby(cone,x,czero,tw,desc_data,info)
! Update with L. The off-diagonal block is taken care
! from the Jacobi smoother, hence this is purely local.
call psb_spmm(-cone,sv%l,xit,cone,tw,desc_data,info,doswap=.false.)
call psb_spsm(cone,sv%u,tw,czero,xit,desc_data,info)
end do
call psb_geaxpby(alpha,xit,beta,y,desc_data,info)
end if
else
!
! Iterations to convergence, not implemented right now.
!
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='EPS>0 not implemented in GS subsolve')
associate(tw => wv(1), xit => wv(2))
select case (init_)
case('Z')
call xit%zero()
case('Y')
call psb_geaxpby(cone,y,czero,xit,desc_data,info)
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,xit,desc_data,info)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end if
end select
case default
select case(trans_)
case('N')
if (sv%eps <=szero) then
!
! Fixed number of iterations
!
!
do itx=1,sv%sweeps
call psb_geaxpby(cone,x,czero,tw,desc_data,info)
! Update with L. The off-diagonal block is taken care
! from the Jacobi smoother, hence this is purely local.
call psb_spmm(-cone,sv%l,xit,cone,tw,desc_data,info,doswap=.false.)
call psb_spsm(cone,sv%u,tw,czero,xit,desc_data,info)
end do
call psb_geaxpby(alpha,xit,beta,y,desc_data,info)
else
!
! Iterations to convergence, not implemented right now.
!
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='EPS>0 not implemented in GS subsolve')
goto 9999
end if
case default
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Invalid TRANS in GS subsolve')
goto 9999
end select
& a_err='Invalid TRANS in GS subsolve')
goto 9999
end select
if (info /= psb_success_) then
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in subsolve')
goto 9999
endif
call tw%free(info)
call xit%free(info)
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in subsolve')
goto 9999
endif
end associate
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
mlprec/impl/solver/mld_c_gs_solver_apply_vect.f90
+61 -55
View File
@@ -54,7 +54,6 @@ subroutine mld_c_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
type(psb_c_vect_type),intent(inout), optional :: initu
integer(psb_ipk_) :: n_row,n_col, itx
type(psb_c_vect_type) :: tw, xit
complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
complex(psb_spk_), allocatable :: temp(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act
@@ -121,69 +120,76 @@ subroutine mld_c_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
goto 9999
end if
call psb_geasb(tw,desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(xit,desc_data,info,mold=x%v,scratch=.true.)
select case (init_)
case('Z')
call xit%zero()
case('Y')
call psb_geaxpby(cone,y,czero,xit,desc_data,info)
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,xit,desc_data,info)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
if (size(wv) < 2) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='invalid wv size')
goto 9999
end select
select case(trans_)
case('N')
if (sv%eps <=szero) then
!
! Fixed number of iterations
!
!
do itx=1,sv%sweeps
call psb_geaxpby(cone,x,czero,tw,desc_data,info)
! Update with U. The off-diagonal block is taken care
! from the Jacobi smoother, hence this is purely local.
call psb_spmm(-cone,sv%u,xit,cone,tw,desc_data,info,doswap=.false.)
call psb_spsm(cone,sv%l,tw,czero,xit,desc_data,info)
end do
call psb_geaxpby(alpha,xit,beta,y,desc_data,info)
end if
else
!
! Iterations to convergence, not implemented right now.
!
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='EPS>0 not implemented in GS subsolve')
associate(tw => wv(1), xit => wv(2))
select case (init_)
case('Z')
call xit%zero()
case('Y')
call psb_geaxpby(cone,y,czero,xit,desc_data,info)
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,xit,desc_data,info)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end if
end select
case default
select case(trans_)
case('N')
if (sv%eps <=szero) then
!
! Fixed number of iterations
!
!
do itx=1,sv%sweeps
call psb_geaxpby(cone,x,czero,tw,desc_data,info)
! Update with U. The off-diagonal block is taken care
! from the Jacobi smoother, hence this is purely local.
call psb_spmm(-cone,sv%u,xit,cone,tw,desc_data,info,doswap=.false.)
call psb_spsm(cone,sv%l,tw,czero,xit,desc_data,info)
end do
call psb_geaxpby(alpha,xit,beta,y,desc_data,info)
else
!
! Iterations to convergence, not implemented right now.
!
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='EPS>0 not implemented in GS subsolve')
goto 9999
end if
case default
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Invalid TRANS in GS subsolve')
goto 9999
end select
& a_err='Invalid TRANS in GS subsolve')
goto 9999
end select
if (info /= psb_success_) then
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in subsolve')
goto 9999
endif
call tw%free(info)
call xit%free(info)
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in subsolve')
goto 9999
endif
end associate
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
mlprec/impl/solver/mld_c_ilu_solver_apply_vect.f90
+45 -37
View File
@@ -125,48 +125,56 @@ subroutine mld_c_ilu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
goto 9999
end if
call psb_geasb(tw,desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(tw1,desc_data,info,mold=x%v,scratch=.true.)
select case(trans_)
case('N')
call psb_spsm(cone,sv%l,x,czero,tw,desc_data,info,&
& trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux)
if (info == psb_success_) call psb_spsm(alpha,sv%u,tw,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_, work=aux)
case('T')
call psb_spsm(cone,sv%u,x,czero,tw,desc_data,info,&
& trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux)
if (info == psb_success_) call psb_spsm(alpha,sv%l,tw,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
case('C')
call psb_spsm(cone,sv%u,x,czero,tw,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
call tw1%mlt(cone,sv%dv,tw,czero,info,conjgx=trans_)
if (info == psb_success_) call psb_spsm(alpha,sv%l,tw1,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid TRANS in ILU subsolve')
if (size(wv) < 2) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='invalid wv size')
goto 9999
end select
end if
associate(tw => wv(1), tw1 => wv(2))
select case(trans_)
case('N')
call psb_spsm(cone,sv%l,x,czero,tw,desc_data,info,&
& trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux)
if (info == psb_success_) call psb_spsm(alpha,sv%u,tw,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_, work=aux)
case('T')
call psb_spsm(cone,sv%u,x,czero,tw,desc_data,info,&
& trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux)
if (info == psb_success_) call psb_spsm(alpha,sv%l,tw,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
case('C')
call psb_spsm(cone,sv%u,x,czero,tw,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
call tw1%mlt(cone,sv%dv,tw,czero,info,conjgx=trans_)
if (info == psb_success_) call psb_spsm(alpha,sv%l,tw1,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid TRANS in ILU subsolve')
goto 9999
end select
if (info /= psb_success_) then
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in subsolve')
goto 9999
endif
call tw%free(info)
call tw1%free(info)
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in subsolve')
goto 9999
endif
end associate
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
mlprec/impl/solver/mld_d_bwgs_solver_apply_vect.f90
+61 -55
View File
@@ -54,7 +54,6 @@ subroutine mld_d_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
type(psb_d_vect_type),intent(inout), optional :: initu
integer(psb_ipk_) :: n_row,n_col, itx
type(psb_d_vect_type) :: tw, xit
real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
real(psb_dpk_), allocatable :: temp(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act
@@ -121,69 +120,76 @@ subroutine mld_d_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
goto 9999
end if
call psb_geasb(tw,desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(xit,desc_data,info,mold=x%v,scratch=.true.)
select case (init_)
case('Z')
call xit%zero()
case('Y')
call psb_geaxpby(done,y,dzero,xit,desc_data,info)
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,xit,desc_data,info)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
if (size(wv) < 2) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='invalid wv size')
goto 9999
end select
select case(trans_)
case('N')
if (sv%eps <=dzero) then
!
! Fixed number of iterations
!
!
do itx=1,sv%sweeps
call psb_geaxpby(done,x,dzero,tw,desc_data,info)
! Update with L. The off-diagonal block is taken care
! from the Jacobi smoother, hence this is purely local.
call psb_spmm(-done,sv%l,xit,done,tw,desc_data,info,doswap=.false.)
call psb_spsm(done,sv%u,tw,dzero,xit,desc_data,info)
end do
call psb_geaxpby(alpha,xit,beta,y,desc_data,info)
end if
else
!
! Iterations to convergence, not implemented right now.
!
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='EPS>0 not implemented in GS subsolve')
associate(tw => wv(1), xit => wv(2))
select case (init_)
case('Z')
call xit%zero()
case('Y')
call psb_geaxpby(done,y,dzero,xit,desc_data,info)
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,xit,desc_data,info)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end if
end select
case default
select case(trans_)
case('N')
if (sv%eps <=dzero) then
!
! Fixed number of iterations
!
!
do itx=1,sv%sweeps
call psb_geaxpby(done,x,dzero,tw,desc_data,info)
! Update with L. The off-diagonal block is taken care
! from the Jacobi smoother, hence this is purely local.
call psb_spmm(-done,sv%l,xit,done,tw,desc_data,info,doswap=.false.)
call psb_spsm(done,sv%u,tw,dzero,xit,desc_data,info)
end do
call psb_geaxpby(alpha,xit,beta,y,desc_data,info)
else
!
! Iterations to convergence, not implemented right now.
!
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='EPS>0 not implemented in GS subsolve')
goto 9999
end if
case default
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Invalid TRANS in GS subsolve')
goto 9999
end select
& a_err='Invalid TRANS in GS subsolve')
goto 9999
end select
if (info /= psb_success_) then
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in subsolve')
goto 9999
endif
call tw%free(info)
call xit%free(info)
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in subsolve')
goto 9999
endif
end associate
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
mlprec/impl/solver/mld_d_gs_solver_apply_vect.f90
+61 -55
View File
@@ -54,7 +54,6 @@ subroutine mld_d_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
type(psb_d_vect_type),intent(inout), optional :: initu
integer(psb_ipk_) :: n_row,n_col, itx
type(psb_d_vect_type) :: tw, xit
real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
real(psb_dpk_), allocatable :: temp(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act
@@ -121,69 +120,76 @@ subroutine mld_d_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
goto 9999
end if
call psb_geasb(tw,desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(xit,desc_data,info,mold=x%v,scratch=.true.)
select case (init_)
case('Z')
call xit%zero()
case('Y')
call psb_geaxpby(done,y,dzero,xit,desc_data,info)
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,xit,desc_data,info)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
if (size(wv) < 2) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='invalid wv size')
goto 9999
end select
select case(trans_)
case('N')
if (sv%eps <=dzero) then
!
! Fixed number of iterations
!
!
do itx=1,sv%sweeps
call psb_geaxpby(done,x,dzero,tw,desc_data,info)
! Update with U. The off-diagonal block is taken care
! from the Jacobi smoother, hence this is purely local.
call psb_spmm(-done,sv%u,xit,done,tw,desc_data,info,doswap=.false.)
call psb_spsm(done,sv%l,tw,dzero,xit,desc_data,info)
end do
call psb_geaxpby(alpha,xit,beta,y,desc_data,info)
end if
else
!
! Iterations to convergence, not implemented right now.
!
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='EPS>0 not implemented in GS subsolve')
associate(tw => wv(1), xit => wv(2))
select case (init_)
case('Z')
call xit%zero()
case('Y')
call psb_geaxpby(done,y,dzero,xit,desc_data,info)
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,xit,desc_data,info)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end if
end select
case default
select case(trans_)
case('N')
if (sv%eps <=dzero) then
!
! Fixed number of iterations
!
!
do itx=1,sv%sweeps
call psb_geaxpby(done,x,dzero,tw,desc_data,info)
! Update with U. The off-diagonal block is taken care
! from the Jacobi smoother, hence this is purely local.
call psb_spmm(-done,sv%u,xit,done,tw,desc_data,info,doswap=.false.)
call psb_spsm(done,sv%l,tw,dzero,xit,desc_data,info)
end do
call psb_geaxpby(alpha,xit,beta,y,desc_data,info)
else
!
! Iterations to convergence, not implemented right now.
!
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='EPS>0 not implemented in GS subsolve')
goto 9999
end if
case default
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Invalid TRANS in GS subsolve')
goto 9999
end select
& a_err='Invalid TRANS in GS subsolve')
goto 9999
end select
if (info /= psb_success_) then
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in subsolve')
goto 9999
endif
call tw%free(info)
call xit%free(info)
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in subsolve')
goto 9999
endif
end associate
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
mlprec/impl/solver/mld_d_ilu_solver_apply_vect.f90
+45 -37
View File
@@ -125,48 +125,56 @@ subroutine mld_d_ilu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
goto 9999
end if
call psb_geasb(tw,desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(tw1,desc_data,info,mold=x%v,scratch=.true.)
select case(trans_)
case('N')
call psb_spsm(done,sv%l,x,dzero,tw,desc_data,info,&
& trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux)
if (info == psb_success_) call psb_spsm(alpha,sv%u,tw,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_, work=aux)
case('T')
call psb_spsm(done,sv%u,x,dzero,tw,desc_data,info,&
& trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux)
if (info == psb_success_) call psb_spsm(alpha,sv%l,tw,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
case('C')
call psb_spsm(done,sv%u,x,dzero,tw,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
call tw1%mlt(done,sv%dv,tw,dzero,info,conjgx=trans_)
if (info == psb_success_) call psb_spsm(alpha,sv%l,tw1,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid TRANS in ILU subsolve')
if (size(wv) < 2) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='invalid wv size')
goto 9999
end select
end if
associate(tw => wv(1), tw1 => wv(2))
select case(trans_)
case('N')
call psb_spsm(done,sv%l,x,dzero,tw,desc_data,info,&
& trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux)
if (info == psb_success_) call psb_spsm(alpha,sv%u,tw,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_, work=aux)
case('T')
call psb_spsm(done,sv%u,x,dzero,tw,desc_data,info,&
& trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux)
if (info == psb_success_) call psb_spsm(alpha,sv%l,tw,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
case('C')
call psb_spsm(done,sv%u,x,dzero,tw,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
call tw1%mlt(done,sv%dv,tw,dzero,info,conjgx=trans_)
if (info == psb_success_) call psb_spsm(alpha,sv%l,tw1,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid TRANS in ILU subsolve')
goto 9999
end select
if (info /= psb_success_) then
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in subsolve')
goto 9999
endif
call tw%free(info)
call tw1%free(info)
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in subsolve')
goto 9999
endif
end associate
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
mlprec/impl/solver/mld_s_bwgs_solver_apply_vect.f90
+61 -55
View File
@@ -54,7 +54,6 @@ subroutine mld_s_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
type(psb_s_vect_type),intent(inout), optional :: initu
integer(psb_ipk_) :: n_row,n_col, itx
type(psb_s_vect_type) :: tw, xit
real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
real(psb_spk_), allocatable :: temp(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act
@@ -121,69 +120,76 @@ subroutine mld_s_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
goto 9999
end if
call psb_geasb(tw,desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(xit,desc_data,info,mold=x%v,scratch=.true.)
select case (init_)
case('Z')
call xit%zero()
case('Y')
call psb_geaxpby(sone,y,szero,xit,desc_data,info)
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,xit,desc_data,info)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
if (size(wv) < 2) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='invalid wv size')
goto 9999
end select
select case(trans_)
case('N')
if (sv%eps <=szero) then
!
! Fixed number of iterations
!
!
do itx=1,sv%sweeps
call psb_geaxpby(sone,x,szero,tw,desc_data,info)
! Update with L. The off-diagonal block is taken care
! from the Jacobi smoother, hence this is purely local.
call psb_spmm(-sone,sv%l,xit,sone,tw,desc_data,info,doswap=.false.)
call psb_spsm(sone,sv%u,tw,szero,xit,desc_data,info)
end do
call psb_geaxpby(alpha,xit,beta,y,desc_data,info)
end if
else
!
! Iterations to convergence, not implemented right now.
!
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='EPS>0 not implemented in GS subsolve')
associate(tw => wv(1), xit => wv(2))
select case (init_)
case('Z')
call xit%zero()
case('Y')
call psb_geaxpby(sone,y,szero,xit,desc_data,info)
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,xit,desc_data,info)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end if
end select
case default
select case(trans_)
case('N')
if (sv%eps <=szero) then
!
! Fixed number of iterations
!
!
do itx=1,sv%sweeps
call psb_geaxpby(sone,x,szero,tw,desc_data,info)
! Update with L. The off-diagonal block is taken care
! from the Jacobi smoother, hence this is purely local.
call psb_spmm(-sone,sv%l,xit,sone,tw,desc_data,info,doswap=.false.)
call psb_spsm(sone,sv%u,tw,szero,xit,desc_data,info)
end do
call psb_geaxpby(alpha,xit,beta,y,desc_data,info)
else
!
! Iterations to convergence, not implemented right now.
!
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='EPS>0 not implemented in GS subsolve')
goto 9999
end if
case default
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Invalid TRANS in GS subsolve')
goto 9999
end select
& a_err='Invalid TRANS in GS subsolve')
goto 9999
end select
if (info /= psb_success_) then
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in subsolve')
goto 9999
endif
call tw%free(info)
call xit%free(info)
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in subsolve')
goto 9999
endif
end associate
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
mlprec/impl/solver/mld_s_gs_solver_apply_vect.f90
+61 -55
View File
@@ -54,7 +54,6 @@ subroutine mld_s_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
type(psb_s_vect_type),intent(inout), optional :: initu
integer(psb_ipk_) :: n_row,n_col, itx
type(psb_s_vect_type) :: tw, xit
real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
real(psb_spk_), allocatable :: temp(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act
@@ -121,69 +120,76 @@ subroutine mld_s_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
goto 9999
end if
call psb_geasb(tw,desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(xit,desc_data,info,mold=x%v,scratch=.true.)
select case (init_)
case('Z')
call xit%zero()
case('Y')
call psb_geaxpby(sone,y,szero,xit,desc_data,info)
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,xit,desc_data,info)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
if (size(wv) < 2) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='invalid wv size')
goto 9999
end select
select case(trans_)
case('N')
if (sv%eps <=szero) then
!
! Fixed number of iterations
!
!
do itx=1,sv%sweeps
call psb_geaxpby(sone,x,szero,tw,desc_data,info)
! Update with U. The off-diagonal block is taken care
! from the Jacobi smoother, hence this is purely local.
call psb_spmm(-sone,sv%u,xit,sone,tw,desc_data,info,doswap=.false.)
call psb_spsm(sone,sv%l,tw,szero,xit,desc_data,info)
end do
call psb_geaxpby(alpha,xit,beta,y,desc_data,info)
end if
else
!
! Iterations to convergence, not implemented right now.
!
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='EPS>0 not implemented in GS subsolve')
associate(tw => wv(1), xit => wv(2))
select case (init_)
case('Z')
call xit%zero()
case('Y')
call psb_geaxpby(sone,y,szero,xit,desc_data,info)
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,xit,desc_data,info)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end if
end select
case default
select case(trans_)
case('N')
if (sv%eps <=szero) then
!
! Fixed number of iterations
!
!
do itx=1,sv%sweeps
call psb_geaxpby(sone,x,szero,tw,desc_data,info)
! Update with U. The off-diagonal block is taken care
! from the Jacobi smoother, hence this is purely local.
call psb_spmm(-sone,sv%u,xit,sone,tw,desc_data,info,doswap=.false.)
call psb_spsm(sone,sv%l,tw,szero,xit,desc_data,info)
end do
call psb_geaxpby(alpha,xit,beta,y,desc_data,info)
else
!
! Iterations to convergence, not implemented right now.
!
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='EPS>0 not implemented in GS subsolve')
goto 9999
end if
case default
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Invalid TRANS in GS subsolve')
goto 9999
end select
& a_err='Invalid TRANS in GS subsolve')
goto 9999
end select
if (info /= psb_success_) then
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in subsolve')
goto 9999
endif
call tw%free(info)
call xit%free(info)
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in subsolve')
goto 9999
endif
end associate
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
mlprec/impl/solver/mld_s_ilu_solver_apply_vect.f90
+45 -37
View File
@@ -125,48 +125,56 @@ subroutine mld_s_ilu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
goto 9999
end if
call psb_geasb(tw,desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(tw1,desc_data,info,mold=x%v,scratch=.true.)
select case(trans_)
case('N')
call psb_spsm(sone,sv%l,x,szero,tw,desc_data,info,&
& trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux)
if (info == psb_success_) call psb_spsm(alpha,sv%u,tw,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_, work=aux)
case('T')
call psb_spsm(sone,sv%u,x,szero,tw,desc_data,info,&
& trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux)
if (info == psb_success_) call psb_spsm(alpha,sv%l,tw,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
case('C')
call psb_spsm(sone,sv%u,x,szero,tw,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
call tw1%mlt(sone,sv%dv,tw,szero,info,conjgx=trans_)
if (info == psb_success_) call psb_spsm(alpha,sv%l,tw1,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid TRANS in ILU subsolve')
if (size(wv) < 2) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='invalid wv size')
goto 9999
end select
end if
associate(tw => wv(1), tw1 => wv(2))
select case(trans_)
case('N')
call psb_spsm(sone,sv%l,x,szero,tw,desc_data,info,&
& trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux)
if (info == psb_success_) call psb_spsm(alpha,sv%u,tw,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_, work=aux)
case('T')
call psb_spsm(sone,sv%u,x,szero,tw,desc_data,info,&
& trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux)
if (info == psb_success_) call psb_spsm(alpha,sv%l,tw,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
case('C')
call psb_spsm(sone,sv%u,x,szero,tw,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
call tw1%mlt(sone,sv%dv,tw,szero,info,conjgx=trans_)
if (info == psb_success_) call psb_spsm(alpha,sv%l,tw1,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid TRANS in ILU subsolve')
goto 9999
end select
if (info /= psb_success_) then
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in subsolve')
goto 9999
endif
call tw%free(info)
call tw1%free(info)
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in subsolve')
goto 9999
endif
end associate
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
mlprec/impl/solver/mld_z_bwgs_solver_apply_vect.f90
+61 -55
View File
@@ -54,7 +54,6 @@ subroutine mld_z_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
type(psb_z_vect_type),intent(inout), optional :: initu
integer(psb_ipk_) :: n_row,n_col, itx
type(psb_z_vect_type) :: tw, xit
complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
complex(psb_dpk_), allocatable :: temp(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act
@@ -121,69 +120,76 @@ subroutine mld_z_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
goto 9999
end if
call psb_geasb(tw,desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(xit,desc_data,info,mold=x%v,scratch=.true.)
select case (init_)
case('Z')
call xit%zero()
case('Y')
call psb_geaxpby(zone,y,zzero,xit,desc_data,info)
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,xit,desc_data,info)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
if (size(wv) < 2) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='invalid wv size')
goto 9999
end select
select case(trans_)
case('N')
if (sv%eps <=dzero) then
!
! Fixed number of iterations
!
!
do itx=1,sv%sweeps
call psb_geaxpby(zone,x,zzero,tw,desc_data,info)
! Update with L. The off-diagonal block is taken care
! from the Jacobi smoother, hence this is purely local.
call psb_spmm(-zone,sv%l,xit,zone,tw,desc_data,info,doswap=.false.)
call psb_spsm(zone,sv%u,tw,zzero,xit,desc_data,info)
end do
call psb_geaxpby(alpha,xit,beta,y,desc_data,info)
end if
else
!
! Iterations to convergence, not implemented right now.
!
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='EPS>0 not implemented in GS subsolve')
associate(tw => wv(1), xit => wv(2))
select case (init_)
case('Z')
call xit%zero()
case('Y')
call psb_geaxpby(zone,y,zzero,xit,desc_data,info)
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,xit,desc_data,info)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end if
end select
case default
select case(trans_)
case('N')
if (sv%eps <=dzero) then
!
! Fixed number of iterations
!
!
do itx=1,sv%sweeps
call psb_geaxpby(zone,x,zzero,tw,desc_data,info)
! Update with L. The off-diagonal block is taken care
! from the Jacobi smoother, hence this is purely local.
call psb_spmm(-zone,sv%l,xit,zone,tw,desc_data,info,doswap=.false.)
call psb_spsm(zone,sv%u,tw,zzero,xit,desc_data,info)
end do
call psb_geaxpby(alpha,xit,beta,y,desc_data,info)
else
!
! Iterations to convergence, not implemented right now.
!
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='EPS>0 not implemented in GS subsolve')
goto 9999
end if
case default
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Invalid TRANS in GS subsolve')
goto 9999
end select
& a_err='Invalid TRANS in GS subsolve')
goto 9999
end select
if (info /= psb_success_) then
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in subsolve')
goto 9999
endif
call tw%free(info)
call xit%free(info)
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in subsolve')
goto 9999
endif
end associate
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
mlprec/impl/solver/mld_z_gs_solver_apply_vect.f90
+61 -55
View File
@@ -54,7 +54,6 @@ subroutine mld_z_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
type(psb_z_vect_type),intent(inout), optional :: initu
integer(psb_ipk_) :: n_row,n_col, itx
type(psb_z_vect_type) :: tw, xit
complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
complex(psb_dpk_), allocatable :: temp(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act
@@ -121,69 +120,76 @@ subroutine mld_z_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
goto 9999
end if
call psb_geasb(tw,desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(xit,desc_data,info,mold=x%v,scratch=.true.)
select case (init_)
case('Z')
call xit%zero()
case('Y')
call psb_geaxpby(zone,y,zzero,xit,desc_data,info)
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,xit,desc_data,info)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
if (size(wv) < 2) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='invalid wv size')
goto 9999
end select
select case(trans_)
case('N')
if (sv%eps <=dzero) then
!
! Fixed number of iterations
!
!
do itx=1,sv%sweeps
call psb_geaxpby(zone,x,zzero,tw,desc_data,info)
! Update with U. The off-diagonal block is taken care
! from the Jacobi smoother, hence this is purely local.
call psb_spmm(-zone,sv%u,xit,zone,tw,desc_data,info,doswap=.false.)
call psb_spsm(zone,sv%l,tw,zzero,xit,desc_data,info)
end do
call psb_geaxpby(alpha,xit,beta,y,desc_data,info)
end if
else
!
! Iterations to convergence, not implemented right now.
!
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='EPS>0 not implemented in GS subsolve')
associate(tw => wv(1), xit => wv(2))
select case (init_)
case('Z')
call xit%zero()
case('Y')
call psb_geaxpby(zone,y,zzero,xit,desc_data,info)
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,xit,desc_data,info)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end if
end select
case default
select case(trans_)
case('N')
if (sv%eps <=dzero) then
!
! Fixed number of iterations
!
!
do itx=1,sv%sweeps
call psb_geaxpby(zone,x,zzero,tw,desc_data,info)
! Update with U. The off-diagonal block is taken care
! from the Jacobi smoother, hence this is purely local.
call psb_spmm(-zone,sv%u,xit,zone,tw,desc_data,info,doswap=.false.)
call psb_spsm(zone,sv%l,tw,zzero,xit,desc_data,info)
end do
call psb_geaxpby(alpha,xit,beta,y,desc_data,info)
else
!
! Iterations to convergence, not implemented right now.
!
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='EPS>0 not implemented in GS subsolve')
goto 9999
end if
case default
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Invalid TRANS in GS subsolve')
goto 9999
end select
& a_err='Invalid TRANS in GS subsolve')
goto 9999
end select
if (info /= psb_success_) then
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in subsolve')
goto 9999
endif
call tw%free(info)
call xit%free(info)
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in subsolve')
goto 9999
endif
end associate
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
mlprec/impl/solver/mld_z_ilu_solver_apply_vect.f90
+45 -37
View File
@@ -125,48 +125,56 @@ subroutine mld_z_ilu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
goto 9999
end if
call psb_geasb(tw,desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(tw1,desc_data,info,mold=x%v,scratch=.true.)
select case(trans_)
case('N')
call psb_spsm(zone,sv%l,x,zzero,tw,desc_data,info,&
& trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux)
if (info == psb_success_) call psb_spsm(alpha,sv%u,tw,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_, work=aux)
case('T')
call psb_spsm(zone,sv%u,x,zzero,tw,desc_data,info,&
& trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux)
if (info == psb_success_) call psb_spsm(alpha,sv%l,tw,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
case('C')
call psb_spsm(zone,sv%u,x,zzero,tw,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
call tw1%mlt(zone,sv%dv,tw,zzero,info,conjgx=trans_)
if (info == psb_success_) call psb_spsm(alpha,sv%l,tw1,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid TRANS in ILU subsolve')
if (size(wv) < 2) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='invalid wv size')
goto 9999
end select
end if
associate(tw => wv(1), tw1 => wv(2))
select case(trans_)
case('N')
call psb_spsm(zone,sv%l,x,zzero,tw,desc_data,info,&
& trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux)
if (info == psb_success_) call psb_spsm(alpha,sv%u,tw,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_, work=aux)
case('T')
call psb_spsm(zone,sv%u,x,zzero,tw,desc_data,info,&
& trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux)
if (info == psb_success_) call psb_spsm(alpha,sv%l,tw,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
case('C')
call psb_spsm(zone,sv%u,x,zzero,tw,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
call tw1%mlt(zone,sv%dv,tw,zzero,info,conjgx=trans_)
if (info == psb_success_) call psb_spsm(alpha,sv%l,tw1,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid TRANS in ILU subsolve')
goto 9999
end select
if (info /= psb_success_) then
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in subsolve')
goto 9999
endif
call tw%free(info)
call tw1%free(info)
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in subsolve')
goto 9999
endif
end associate
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else