diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90
index a7a4202f..b248a460 100644
--- a/amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90
+++ b/amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90
@@ -108,11 +108,6 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
       goto 9999
     end if
   endif
-!!$  if (me == 0) write(0,*) name,' Unimplemented apply_vect '
-!!$  info =psb_err_internal_error_
-!!$  call psb_errpush(info,&
-!!$       & name,a_err='Error in sub_aply Polynomial not implemented')
-!!$  goto 9999
 
   if (size(wv) < 4) then
     info = psb_err_internal_error_
@@ -150,8 +145,8 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
           else
             call psb_spmm(-done,sm%pa,tz,done,r,desc_data,info,work=aux,trans=trans_)
           end if
-!!$        res  = psb_genrm2(r,desc_data,info)
-!!$        write(0,*) 'Polynomial smoother ',i,res
+!!$          res  = psb_genrm2(r,desc_data,info)
+!!$          write(0,*) 'Polynomial smoother LOTTES ',i,res
           ! x_k =  x_{k-1} + z_k
         end do
       end block
@@ -186,8 +181,8 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
           else
             call psb_spmm(-done,sm%pa,tz,done,r,desc_data,info,work=aux,trans=trans_)
           end if
-!!$        res  = psb_genrm2(r,desc_data,info)
-!!$        write(0,*) 'Polynomial smoother ',i,res
+!!$          res  = psb_genrm2(r,desc_data,info)
+!!$          write(0,*) 'Polynomial smoother LOTTES_BETA ',i,res
           ! x_k =  x_{k-1} + z_k
         end do
       end block
@@ -221,7 +216,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
           rho = done/(2*sigma - rho_old)
           call psb_geaxpby((2*rho/delta),r,(rho*rho_old),tz,desc_data,info)
 !!$          res  = psb_genrm2(r,desc_data,info)
-!!$          write(0,*) 'Polynomial smoother ',i,res
+!!$          write(0,*) 'Polynomial smoother NEW ',i,res
           ! x_k =  x_{k-1} + z_k
         end do
       end block
@@ -244,205 +239,10 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
     end if
   end associate
 
-
-
-
-!!$  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 (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 (amg_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')
-!!$
-!!$        case default
-!!$          call psb_errpush(psb_err_internal_error_,name,&
-!!$               & a_err='wrong  init to smoother_apply')
-!!$          goto 9999
-!!$        end select
-!!$
-!!$        do i=1, sweeps-1
-!!$          !
-!!$          ! Compute Y(j+1) =  Y(j)+ D^(-1)*(X-A*Y(j)),
-!!$          !  where is the diagonal and  A the matrix.
-!!$          !
-!!$          call psb_geaxpby(done,x,dzero,tx,desc_data,info)
-!!$          call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
-!!$
-!!$          if (info /= psb_success_) exit
-!!$
-!!$          call sm%sv%apply(done,tx,done,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
-!!$
-!!$          if (info /= psb_success_) exit
-!!$
-!!$          if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
-!!$            call psb_geaxpby(done,x,dzero,r,r,desc_data,info)
-!!$            call psb_spmm(-done,sm%pa,ty,done,r,desc_data,info)
-!!$            res  = psb_genrm2(r,desc_data,info)
-!!$            if( sm%printres ) then
-!!$              call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
-!!$            end if
-!!$            if ( res < sm%tol*resdenum ) then
-!!$              if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) &
-!!$               & call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
-!!$              exit
-!!$            end if
-!!$          end if
-!!$
-!!$        end do
-!!$
-!!$
-!!$        if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
-!!$
-!!$        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
-!!$
-!!$    class default
-!!$      !
-!!$      !
-!!$      ! Apply multiple sweeps of a block-Jacobi solver
-!!$      ! to compute an approximate solution of a linear system.
-!!$      !
-!!$      !
-!!$      if (size(wv) < 2) then
-!!$        info = psb_err_internal_error_
-!!$        call psb_errpush(info,name,&
-!!$             & a_err='invalid wv size in smoother_apply')
-!!$        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 select
-!!$
-!!$  else
-!!$
-!!$    info = psb_err_iarg_neg_
-!!$    call psb_errpush(info,name,&
-!!$         & i_err=(/itwo,sweeps,izero,izero,izero/))
-!!$    goto 9999
-!!$
-!!$  endif
-!!$
   if (.not.(4*n_col <= size(work))) then
     deallocate(aux)
   endif
 
-!!$  if(sm%checkres) then
-!!$    call psb_gefree(r,desc_data,info)
-!!$  end if
-
   call psb_erractionrestore(err_act)
   return
 
diff --git a/amgprec/impl/smoother/amg_s_poly_smoother_apply_vect.f90 b/amgprec/impl/smoother/amg_s_poly_smoother_apply_vect.f90
index 76be3e99..835e7eae 100644
--- a/amgprec/impl/smoother/amg_s_poly_smoother_apply_vect.f90
+++ b/amgprec/impl/smoother/amg_s_poly_smoother_apply_vect.f90
@@ -108,11 +108,6 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
       goto 9999
     end if
   endif
-!!$  if (me == 0) write(0,*) name,' Unimplemented apply_vect '
-!!$  info =psb_err_internal_error_
-!!$  call psb_errpush(info,&
-!!$       & name,a_err='Error in sub_aply Polynomial not implemented')
-!!$  goto 9999
 
   if (size(wv) < 4) then
     info = psb_err_internal_error_
@@ -150,8 +145,8 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
           else
             call psb_spmm(-sone,sm%pa,tz,sone,r,desc_data,info,work=aux,trans=trans_)
           end if
-!!$        res  = psb_genrm2(r,desc_data,info)
-!!$        write(0,*) 'Polynomial smoother ',i,res
+!!$          res  = psb_genrm2(r,desc_data,info)
+!!$          write(0,*) 'Polynomial smoother LOTTES ',i,res
           ! x_k =  x_{k-1} + z_k
         end do
       end block
@@ -186,8 +181,8 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
           else
             call psb_spmm(-sone,sm%pa,tz,sone,r,desc_data,info,work=aux,trans=trans_)
           end if
-!!$        res  = psb_genrm2(r,desc_data,info)
-!!$        write(0,*) 'Polynomial smoother ',i,res
+!!$          res  = psb_genrm2(r,desc_data,info)
+!!$          write(0,*) 'Polynomial smoother LOTTES_BETA ',i,res
           ! x_k =  x_{k-1} + z_k
         end do
       end block
@@ -221,7 +216,7 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
           rho = sone/(2*sigma - rho_old)
           call psb_geaxpby((2*rho/delta),r,(rho*rho_old),tz,desc_data,info)
 !!$          res  = psb_genrm2(r,desc_data,info)
-!!$          write(0,*) 'Polynomial smoother ',i,res
+!!$          write(0,*) 'Polynomial smoother NEW ',i,res
           ! x_k =  x_{k-1} + z_k
         end do
       end block
@@ -244,205 +239,10 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
     end if
   end associate
 
-
-
-
-!!$  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 (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 (amg_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')
-!!$
-!!$        case default
-!!$          call psb_errpush(psb_err_internal_error_,name,&
-!!$               & a_err='wrong  init to smoother_apply')
-!!$          goto 9999
-!!$        end select
-!!$
-!!$        do i=1, sweeps-1
-!!$          !
-!!$          ! Compute Y(j+1) =  Y(j)+ D^(-1)*(X-A*Y(j)),
-!!$          !  where is the diagonal and  A the matrix.
-!!$          !
-!!$          call psb_geaxpby(sone,x,szero,tx,desc_data,info)
-!!$          call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
-!!$
-!!$          if (info /= psb_success_) exit
-!!$
-!!$          call sm%sv%apply(sone,tx,sone,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
-!!$
-!!$          if (info /= psb_success_) exit
-!!$
-!!$          if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
-!!$            call psb_geaxpby(sone,x,szero,r,r,desc_data,info)
-!!$            call psb_spmm(-sone,sm%pa,ty,sone,r,desc_data,info)
-!!$            res  = psb_genrm2(r,desc_data,info)
-!!$            if( sm%printres ) then
-!!$              call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
-!!$            end if
-!!$            if ( res < sm%tol*resdenum ) then
-!!$              if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) &
-!!$               & call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
-!!$              exit
-!!$            end if
-!!$          end if
-!!$
-!!$        end do
-!!$
-!!$
-!!$        if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
-!!$
-!!$        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
-!!$
-!!$    class default
-!!$      !
-!!$      !
-!!$      ! Apply multiple sweeps of a block-Jacobi solver
-!!$      ! to compute an approximate solution of a linear system.
-!!$      !
-!!$      !
-!!$      if (size(wv) < 2) then
-!!$        info = psb_err_internal_error_
-!!$        call psb_errpush(info,name,&
-!!$             & a_err='invalid wv size in smoother_apply')
-!!$        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 select
-!!$
-!!$  else
-!!$
-!!$    info = psb_err_iarg_neg_
-!!$    call psb_errpush(info,name,&
-!!$         & i_err=(/itwo,sweeps,izero,izero,izero/))
-!!$    goto 9999
-!!$
-!!$  endif
-!!$
   if (.not.(4*n_col <= size(work))) then
     deallocate(aux)
   endif
 
-!!$  if(sm%checkres) then
-!!$    call psb_gefree(r,desc_data,info)
-!!$  end if
-
   call psb_erractionrestore(err_act)
   return