diff --git a/mlprec/impl/solver/mld_c_bwgs_solver_apply.f90 b/mlprec/impl/solver/mld_c_bwgs_solver_apply.f90
index fcd0f93f..e70bd789 100644
--- a/mlprec/impl/solver/mld_c_bwgs_solver_apply.f90
+++ b/mlprec/impl/solver/mld_c_bwgs_solver_apply.f90
@@ -52,7 +52,7 @@ subroutine mld_c_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,&
   character, intent(in), optional       :: init
   complex(psb_spk_),intent(inout), optional :: initu(:)
 
-  integer(psb_ipk_)  :: n_row,n_col, itx
+  integer(psb_ipk_)  :: n_row,n_col, itx, itxst
   complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   complex(psb_spk_), allocatable :: temp(:),wv(:),xit(:)
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
@@ -118,10 +118,13 @@ subroutine mld_c_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,&
   end if
 
   call psb_geasb(wv,desc_data,info) 
-  call psb_geasb(xit,desc_data,info) 
+  call psb_geasb(xit,desc_data,info)
+  itxst = 1
   select case (init_)
   case('Z') 
-    xit(:) = czero
+    call psb_geaxpby(cone,x,czero,wv,desc_data,info)
+    call psb_spsm(cone,sv%u,wv,czero,xit,desc_data,info)
+    itxst = 2
   case('Y')
     call psb_geaxpby(cone,y,czero,xit,desc_data,info)
   case('U')
@@ -144,7 +147,7 @@ subroutine mld_c_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,&
       ! Fixed number of iterations
       !
       !
-      do itx=1,sv%sweeps
+      do itx=itxst, sv%sweeps
         call psb_geaxpby(cone,x,czero,wv,desc_data,info)
         ! Update with L. The off-diagonal block is taken care
         ! from the Jacobi smoother, hence this is purely local. 
diff --git a/mlprec/impl/solver/mld_c_bwgs_solver_apply_vect.f90 b/mlprec/impl/solver/mld_c_bwgs_solver_apply_vect.f90
index 33d21fd1..41d645fc 100644
--- a/mlprec/impl/solver/mld_c_bwgs_solver_apply_vect.f90
+++ b/mlprec/impl/solver/mld_c_bwgs_solver_apply_vect.f90
@@ -53,7 +53,7 @@ subroutine mld_c_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
   character, intent(in), optional                :: init
   type(psb_c_vect_type),intent(inout), optional   :: initu
 
-  integer(psb_ipk_)   :: n_row,n_col, itx
+  integer(psb_ipk_)   :: n_row,n_col, itx, itxst
   complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   complex(psb_spk_), allocatable :: temp(:)
   integer(psb_ipk_)   :: ictxt,np,me,i, err_act
@@ -128,10 +128,12 @@ subroutine mld_c_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
   end if
 
   associate(tw => wv(1), xit => wv(2))
-
+    itxst = 1
     select case (init_)
     case('Z') 
-      call xit%zero()
+      call psb_geaxpby(cone,x,czero,tw,desc_data,info)
+      call psb_spsm(cone,sv%u,tw,czero,xit,desc_data,info)
+      itxst = 2
     case('Y')
       call psb_geaxpby(cone,y,czero,xit,desc_data,info)
     case('U')
@@ -154,7 +156,7 @@ subroutine mld_c_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
         ! Fixed number of iterations
         !
         !
-        do itx=1,sv%sweeps
+        do itx=itxst,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. 
diff --git a/mlprec/impl/solver/mld_c_gs_solver_apply.f90 b/mlprec/impl/solver/mld_c_gs_solver_apply.f90
index c5b15321..2f053b08 100644
--- a/mlprec/impl/solver/mld_c_gs_solver_apply.f90
+++ b/mlprec/impl/solver/mld_c_gs_solver_apply.f90
@@ -52,7 +52,7 @@ subroutine mld_c_gs_solver_apply(alpha,sv,x,beta,y,desc_data,&
   character, intent(in), optional       :: init
   complex(psb_spk_),intent(inout), optional :: initu(:)
 
-  integer(psb_ipk_)  :: n_row,n_col, itx
+  integer(psb_ipk_)  :: n_row,n_col, itx, itxst
   complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   complex(psb_spk_), allocatable :: temp(:),wv(:),xit(:)
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
@@ -118,10 +118,13 @@ subroutine mld_c_gs_solver_apply(alpha,sv,x,beta,y,desc_data,&
   end if
 
   call psb_geasb(wv,desc_data,info) 
-  call psb_geasb(xit,desc_data,info) 
+  call psb_geasb(xit,desc_data,info)
+  itxst = 1
   select case (init_)
   case('Z') 
-    xit(:) = czero
+    call psb_geaxpby(cone,x,czero,wv,desc_data,info)
+    call psb_spsm(cone,sv%l,wv,czero,xit,desc_data,info)
+    itxst = 2 
   case('Y')
     call psb_geaxpby(cone,y,czero,xit,desc_data,info)
   case('U')
@@ -144,7 +147,7 @@ subroutine mld_c_gs_solver_apply(alpha,sv,x,beta,y,desc_data,&
       ! Fixed number of iterations
       !
       !
-      do itx=1,sv%sweeps
+      do itx=itxst,sv%sweeps
         call psb_geaxpby(cone,x,czero,wv,desc_data,info)
         ! Update with U. The off-diagonal block is taken care
         ! from the Jacobi smoother, hence this is purely local. 
diff --git a/mlprec/impl/solver/mld_c_gs_solver_apply_vect.f90 b/mlprec/impl/solver/mld_c_gs_solver_apply_vect.f90
index d0c9c326..e322c164 100644
--- a/mlprec/impl/solver/mld_c_gs_solver_apply_vect.f90
+++ b/mlprec/impl/solver/mld_c_gs_solver_apply_vect.f90
@@ -53,7 +53,7 @@ subroutine mld_c_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
   character, intent(in), optional                :: init
   type(psb_c_vect_type),intent(inout), optional   :: initu
 
-  integer(psb_ipk_)   :: n_row,n_col, itx
+  integer(psb_ipk_)   :: n_row,n_col, itx, itxst
   complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   complex(psb_spk_), allocatable :: temp(:)
   integer(psb_ipk_)   :: ictxt,np,me,i, err_act
@@ -128,10 +128,12 @@ subroutine mld_c_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
   end if
 
   associate(tw => wv(1), xit => wv(2))
-
+    itxst = 1 
     select case (init_)
     case('Z') 
-      call xit%zero()
+      call psb_geaxpby(cone,x,czero,tw,desc_data,info)
+      call psb_spsm(cone,sv%l,tw,czero,xit,desc_data,info)
+      itxst = 2
     case('Y')
       call psb_geaxpby(cone,y,czero,xit,desc_data,info)
     case('U')
@@ -154,7 +156,7 @@ subroutine mld_c_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
         ! Fixed number of iterations
         !
         !
-        do itx=1,sv%sweeps
+        do itx=itxst,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. 
diff --git a/mlprec/impl/solver/mld_d_bwgs_solver_apply.f90 b/mlprec/impl/solver/mld_d_bwgs_solver_apply.f90
index 4ce02fea..6eaece03 100644
--- a/mlprec/impl/solver/mld_d_bwgs_solver_apply.f90
+++ b/mlprec/impl/solver/mld_d_bwgs_solver_apply.f90
@@ -52,7 +52,7 @@ subroutine mld_d_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,&
   character, intent(in), optional       :: init
   real(psb_dpk_),intent(inout), optional :: initu(:)
 
-  integer(psb_ipk_)  :: n_row,n_col, itx
+  integer(psb_ipk_)  :: n_row,n_col, itx, itxst
   real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   real(psb_dpk_), allocatable :: temp(:),wv(:),xit(:)
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
@@ -118,10 +118,13 @@ subroutine mld_d_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,&
   end if
 
   call psb_geasb(wv,desc_data,info) 
-  call psb_geasb(xit,desc_data,info) 
+  call psb_geasb(xit,desc_data,info)
+  itxst = 1
   select case (init_)
   case('Z') 
-    xit(:) = dzero
+    call psb_geaxpby(done,x,dzero,wv,desc_data,info)
+    call psb_spsm(done,sv%u,wv,dzero,xit,desc_data,info)
+    itxst = 2
   case('Y')
     call psb_geaxpby(done,y,dzero,xit,desc_data,info)
   case('U')
@@ -144,7 +147,7 @@ subroutine mld_d_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,&
       ! Fixed number of iterations
       !
       !
-      do itx=1,sv%sweeps
+      do itx=itxst, sv%sweeps
         call psb_geaxpby(done,x,dzero,wv,desc_data,info)
         ! Update with L. The off-diagonal block is taken care
         ! from the Jacobi smoother, hence this is purely local. 
diff --git a/mlprec/impl/solver/mld_d_bwgs_solver_apply_vect.f90 b/mlprec/impl/solver/mld_d_bwgs_solver_apply_vect.f90
index fe4fdd7c..6ee62fef 100644
--- a/mlprec/impl/solver/mld_d_bwgs_solver_apply_vect.f90
+++ b/mlprec/impl/solver/mld_d_bwgs_solver_apply_vect.f90
@@ -53,7 +53,7 @@ subroutine mld_d_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
   character, intent(in), optional                :: init
   type(psb_d_vect_type),intent(inout), optional   :: initu
 
-  integer(psb_ipk_)   :: n_row,n_col, itx
+  integer(psb_ipk_)   :: n_row,n_col, itx, itxst
   real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   real(psb_dpk_), allocatable :: temp(:)
   integer(psb_ipk_)   :: ictxt,np,me,i, err_act
@@ -128,10 +128,12 @@ subroutine mld_d_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
   end if
 
   associate(tw => wv(1), xit => wv(2))
-
+    itxst = 1
     select case (init_)
     case('Z') 
-      call xit%zero()
+      call psb_geaxpby(done,x,dzero,tw,desc_data,info)
+      call psb_spsm(done,sv%u,tw,dzero,xit,desc_data,info)
+      itxst = 2
     case('Y')
       call psb_geaxpby(done,y,dzero,xit,desc_data,info)
     case('U')
@@ -154,7 +156,7 @@ subroutine mld_d_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
         ! Fixed number of iterations
         !
         !
-        do itx=1,sv%sweeps
+        do itx=itxst,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. 
diff --git a/mlprec/impl/solver/mld_d_gs_solver_apply.f90 b/mlprec/impl/solver/mld_d_gs_solver_apply.f90
index 7e5e2da8..0459a523 100644
--- a/mlprec/impl/solver/mld_d_gs_solver_apply.f90
+++ b/mlprec/impl/solver/mld_d_gs_solver_apply.f90
@@ -52,7 +52,7 @@ subroutine mld_d_gs_solver_apply(alpha,sv,x,beta,y,desc_data,&
   character, intent(in), optional       :: init
   real(psb_dpk_),intent(inout), optional :: initu(:)
 
-  integer(psb_ipk_)  :: n_row,n_col, itx
+  integer(psb_ipk_)  :: n_row,n_col, itx, itxst
   real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   real(psb_dpk_), allocatable :: temp(:),wv(:),xit(:)
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
@@ -118,10 +118,13 @@ subroutine mld_d_gs_solver_apply(alpha,sv,x,beta,y,desc_data,&
   end if
 
   call psb_geasb(wv,desc_data,info) 
-  call psb_geasb(xit,desc_data,info) 
+  call psb_geasb(xit,desc_data,info)
+  itxst = 1
   select case (init_)
   case('Z') 
-    xit(:) = dzero
+    call psb_geaxpby(done,x,dzero,wv,desc_data,info)
+    call psb_spsm(done,sv%l,wv,dzero,xit,desc_data,info)
+    itxst = 2 
   case('Y')
     call psb_geaxpby(done,y,dzero,xit,desc_data,info)
   case('U')
@@ -144,7 +147,7 @@ subroutine mld_d_gs_solver_apply(alpha,sv,x,beta,y,desc_data,&
       ! Fixed number of iterations
       !
       !
-      do itx=1,sv%sweeps
+      do itx=itxst,sv%sweeps
         call psb_geaxpby(done,x,dzero,wv,desc_data,info)
         ! Update with U. The off-diagonal block is taken care
         ! from the Jacobi smoother, hence this is purely local. 
diff --git a/mlprec/impl/solver/mld_d_gs_solver_apply_vect.f90 b/mlprec/impl/solver/mld_d_gs_solver_apply_vect.f90
index 5ed87f5b..a5e30f44 100644
--- a/mlprec/impl/solver/mld_d_gs_solver_apply_vect.f90
+++ b/mlprec/impl/solver/mld_d_gs_solver_apply_vect.f90
@@ -53,7 +53,7 @@ subroutine mld_d_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
   character, intent(in), optional                :: init
   type(psb_d_vect_type),intent(inout), optional   :: initu
 
-  integer(psb_ipk_)   :: n_row,n_col, itx
+  integer(psb_ipk_)   :: n_row,n_col, itx, itxst
   real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   real(psb_dpk_), allocatable :: temp(:)
   integer(psb_ipk_)   :: ictxt,np,me,i, err_act
@@ -128,10 +128,12 @@ subroutine mld_d_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
   end if
 
   associate(tw => wv(1), xit => wv(2))
-
+    itxst = 1 
     select case (init_)
     case('Z') 
-      call xit%zero()
+      call psb_geaxpby(done,x,dzero,tw,desc_data,info)
+      call psb_spsm(done,sv%l,tw,dzero,xit,desc_data,info)
+      itxst = 2
     case('Y')
       call psb_geaxpby(done,y,dzero,xit,desc_data,info)
     case('U')
@@ -154,7 +156,7 @@ subroutine mld_d_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
         ! Fixed number of iterations
         !
         !
-        do itx=1,sv%sweeps
+        do itx=itxst,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. 
diff --git a/mlprec/impl/solver/mld_s_bwgs_solver_apply.f90 b/mlprec/impl/solver/mld_s_bwgs_solver_apply.f90
index f6955fbd..e232e393 100644
--- a/mlprec/impl/solver/mld_s_bwgs_solver_apply.f90
+++ b/mlprec/impl/solver/mld_s_bwgs_solver_apply.f90
@@ -52,7 +52,7 @@ subroutine mld_s_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,&
   character, intent(in), optional       :: init
   real(psb_spk_),intent(inout), optional :: initu(:)
 
-  integer(psb_ipk_)  :: n_row,n_col, itx
+  integer(psb_ipk_)  :: n_row,n_col, itx, itxst
   real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   real(psb_spk_), allocatable :: temp(:),wv(:),xit(:)
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
@@ -118,10 +118,13 @@ subroutine mld_s_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,&
   end if
 
   call psb_geasb(wv,desc_data,info) 
-  call psb_geasb(xit,desc_data,info) 
+  call psb_geasb(xit,desc_data,info)
+  itxst = 1
   select case (init_)
   case('Z') 
-    xit(:) = szero
+    call psb_geaxpby(sone,x,szero,wv,desc_data,info)
+    call psb_spsm(sone,sv%u,wv,szero,xit,desc_data,info)
+    itxst = 2
   case('Y')
     call psb_geaxpby(sone,y,szero,xit,desc_data,info)
   case('U')
@@ -144,7 +147,7 @@ subroutine mld_s_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,&
       ! Fixed number of iterations
       !
       !
-      do itx=1,sv%sweeps
+      do itx=itxst, sv%sweeps
         call psb_geaxpby(sone,x,szero,wv,desc_data,info)
         ! Update with L. The off-diagonal block is taken care
         ! from the Jacobi smoother, hence this is purely local. 
diff --git a/mlprec/impl/solver/mld_s_bwgs_solver_apply_vect.f90 b/mlprec/impl/solver/mld_s_bwgs_solver_apply_vect.f90
index a4565edb..7f043d74 100644
--- a/mlprec/impl/solver/mld_s_bwgs_solver_apply_vect.f90
+++ b/mlprec/impl/solver/mld_s_bwgs_solver_apply_vect.f90
@@ -53,7 +53,7 @@ subroutine mld_s_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
   character, intent(in), optional                :: init
   type(psb_s_vect_type),intent(inout), optional   :: initu
 
-  integer(psb_ipk_)   :: n_row,n_col, itx
+  integer(psb_ipk_)   :: n_row,n_col, itx, itxst
   real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   real(psb_spk_), allocatable :: temp(:)
   integer(psb_ipk_)   :: ictxt,np,me,i, err_act
@@ -128,10 +128,12 @@ subroutine mld_s_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
   end if
 
   associate(tw => wv(1), xit => wv(2))
-
+    itxst = 1
     select case (init_)
     case('Z') 
-      call xit%zero()
+      call psb_geaxpby(sone,x,szero,tw,desc_data,info)
+      call psb_spsm(sone,sv%u,tw,szero,xit,desc_data,info)
+      itxst = 2
     case('Y')
       call psb_geaxpby(sone,y,szero,xit,desc_data,info)
     case('U')
@@ -154,7 +156,7 @@ subroutine mld_s_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
         ! Fixed number of iterations
         !
         !
-        do itx=1,sv%sweeps
+        do itx=itxst,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. 
diff --git a/mlprec/impl/solver/mld_s_gs_solver_apply.f90 b/mlprec/impl/solver/mld_s_gs_solver_apply.f90
index adfd751d..0148d03c 100644
--- a/mlprec/impl/solver/mld_s_gs_solver_apply.f90
+++ b/mlprec/impl/solver/mld_s_gs_solver_apply.f90
@@ -52,7 +52,7 @@ subroutine mld_s_gs_solver_apply(alpha,sv,x,beta,y,desc_data,&
   character, intent(in), optional       :: init
   real(psb_spk_),intent(inout), optional :: initu(:)
 
-  integer(psb_ipk_)  :: n_row,n_col, itx
+  integer(psb_ipk_)  :: n_row,n_col, itx, itxst
   real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   real(psb_spk_), allocatable :: temp(:),wv(:),xit(:)
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
@@ -118,10 +118,13 @@ subroutine mld_s_gs_solver_apply(alpha,sv,x,beta,y,desc_data,&
   end if
 
   call psb_geasb(wv,desc_data,info) 
-  call psb_geasb(xit,desc_data,info) 
+  call psb_geasb(xit,desc_data,info)
+  itxst = 1
   select case (init_)
   case('Z') 
-    xit(:) = szero
+    call psb_geaxpby(sone,x,szero,wv,desc_data,info)
+    call psb_spsm(sone,sv%l,wv,szero,xit,desc_data,info)
+    itxst = 2 
   case('Y')
     call psb_geaxpby(sone,y,szero,xit,desc_data,info)
   case('U')
@@ -144,7 +147,7 @@ subroutine mld_s_gs_solver_apply(alpha,sv,x,beta,y,desc_data,&
       ! Fixed number of iterations
       !
       !
-      do itx=1,sv%sweeps
+      do itx=itxst,sv%sweeps
         call psb_geaxpby(sone,x,szero,wv,desc_data,info)
         ! Update with U. The off-diagonal block is taken care
         ! from the Jacobi smoother, hence this is purely local. 
diff --git a/mlprec/impl/solver/mld_s_gs_solver_apply_vect.f90 b/mlprec/impl/solver/mld_s_gs_solver_apply_vect.f90
index 74bfe835..a52d4edf 100644
--- a/mlprec/impl/solver/mld_s_gs_solver_apply_vect.f90
+++ b/mlprec/impl/solver/mld_s_gs_solver_apply_vect.f90
@@ -53,7 +53,7 @@ subroutine mld_s_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
   character, intent(in), optional                :: init
   type(psb_s_vect_type),intent(inout), optional   :: initu
 
-  integer(psb_ipk_)   :: n_row,n_col, itx
+  integer(psb_ipk_)   :: n_row,n_col, itx, itxst
   real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   real(psb_spk_), allocatable :: temp(:)
   integer(psb_ipk_)   :: ictxt,np,me,i, err_act
@@ -128,10 +128,12 @@ subroutine mld_s_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
   end if
 
   associate(tw => wv(1), xit => wv(2))
-
+    itxst = 1 
     select case (init_)
     case('Z') 
-      call xit%zero()
+      call psb_geaxpby(sone,x,szero,tw,desc_data,info)
+      call psb_spsm(sone,sv%l,tw,szero,xit,desc_data,info)
+      itxst = 2
     case('Y')
       call psb_geaxpby(sone,y,szero,xit,desc_data,info)
     case('U')
@@ -154,7 +156,7 @@ subroutine mld_s_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
         ! Fixed number of iterations
         !
         !
-        do itx=1,sv%sweeps
+        do itx=itxst,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. 
diff --git a/mlprec/impl/solver/mld_z_bwgs_solver_apply.f90 b/mlprec/impl/solver/mld_z_bwgs_solver_apply.f90
index 02dce73c..89f870c8 100644
--- a/mlprec/impl/solver/mld_z_bwgs_solver_apply.f90
+++ b/mlprec/impl/solver/mld_z_bwgs_solver_apply.f90
@@ -52,7 +52,7 @@ subroutine mld_z_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,&
   character, intent(in), optional       :: init
   complex(psb_dpk_),intent(inout), optional :: initu(:)
 
-  integer(psb_ipk_)  :: n_row,n_col, itx
+  integer(psb_ipk_)  :: n_row,n_col, itx, itxst
   complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   complex(psb_dpk_), allocatable :: temp(:),wv(:),xit(:)
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
@@ -118,10 +118,13 @@ subroutine mld_z_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,&
   end if
 
   call psb_geasb(wv,desc_data,info) 
-  call psb_geasb(xit,desc_data,info) 
+  call psb_geasb(xit,desc_data,info)
+  itxst = 1
   select case (init_)
   case('Z') 
-    xit(:) = zzero
+    call psb_geaxpby(zone,x,zzero,wv,desc_data,info)
+    call psb_spsm(zone,sv%u,wv,zzero,xit,desc_data,info)
+    itxst = 2
   case('Y')
     call psb_geaxpby(zone,y,zzero,xit,desc_data,info)
   case('U')
@@ -144,7 +147,7 @@ subroutine mld_z_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,&
       ! Fixed number of iterations
       !
       !
-      do itx=1,sv%sweeps
+      do itx=itxst, sv%sweeps
         call psb_geaxpby(zone,x,zzero,wv,desc_data,info)
         ! Update with L. The off-diagonal block is taken care
         ! from the Jacobi smoother, hence this is purely local. 
diff --git a/mlprec/impl/solver/mld_z_bwgs_solver_apply_vect.f90 b/mlprec/impl/solver/mld_z_bwgs_solver_apply_vect.f90
index 9454046f..5e167ba4 100644
--- a/mlprec/impl/solver/mld_z_bwgs_solver_apply_vect.f90
+++ b/mlprec/impl/solver/mld_z_bwgs_solver_apply_vect.f90
@@ -53,7 +53,7 @@ subroutine mld_z_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
   character, intent(in), optional                :: init
   type(psb_z_vect_type),intent(inout), optional   :: initu
 
-  integer(psb_ipk_)   :: n_row,n_col, itx
+  integer(psb_ipk_)   :: n_row,n_col, itx, itxst
   complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   complex(psb_dpk_), allocatable :: temp(:)
   integer(psb_ipk_)   :: ictxt,np,me,i, err_act
@@ -128,10 +128,12 @@ subroutine mld_z_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
   end if
 
   associate(tw => wv(1), xit => wv(2))
-
+    itxst = 1
     select case (init_)
     case('Z') 
-      call xit%zero()
+      call psb_geaxpby(zone,x,zzero,tw,desc_data,info)
+      call psb_spsm(zone,sv%u,tw,zzero,xit,desc_data,info)
+      itxst = 2
     case('Y')
       call psb_geaxpby(zone,y,zzero,xit,desc_data,info)
     case('U')
@@ -154,7 +156,7 @@ subroutine mld_z_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
         ! Fixed number of iterations
         !
         !
-        do itx=1,sv%sweeps
+        do itx=itxst,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. 
diff --git a/mlprec/impl/solver/mld_z_gs_solver_apply.f90 b/mlprec/impl/solver/mld_z_gs_solver_apply.f90
index 4acfff15..679a8867 100644
--- a/mlprec/impl/solver/mld_z_gs_solver_apply.f90
+++ b/mlprec/impl/solver/mld_z_gs_solver_apply.f90
@@ -52,7 +52,7 @@ subroutine mld_z_gs_solver_apply(alpha,sv,x,beta,y,desc_data,&
   character, intent(in), optional       :: init
   complex(psb_dpk_),intent(inout), optional :: initu(:)
 
-  integer(psb_ipk_)  :: n_row,n_col, itx
+  integer(psb_ipk_)  :: n_row,n_col, itx, itxst
   complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   complex(psb_dpk_), allocatable :: temp(:),wv(:),xit(:)
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
@@ -118,10 +118,13 @@ subroutine mld_z_gs_solver_apply(alpha,sv,x,beta,y,desc_data,&
   end if
 
   call psb_geasb(wv,desc_data,info) 
-  call psb_geasb(xit,desc_data,info) 
+  call psb_geasb(xit,desc_data,info)
+  itxst = 1
   select case (init_)
   case('Z') 
-    xit(:) = zzero
+    call psb_geaxpby(zone,x,zzero,wv,desc_data,info)
+    call psb_spsm(zone,sv%l,wv,zzero,xit,desc_data,info)
+    itxst = 2 
   case('Y')
     call psb_geaxpby(zone,y,zzero,xit,desc_data,info)
   case('U')
@@ -144,7 +147,7 @@ subroutine mld_z_gs_solver_apply(alpha,sv,x,beta,y,desc_data,&
       ! Fixed number of iterations
       !
       !
-      do itx=1,sv%sweeps
+      do itx=itxst,sv%sweeps
         call psb_geaxpby(zone,x,zzero,wv,desc_data,info)
         ! Update with U. The off-diagonal block is taken care
         ! from the Jacobi smoother, hence this is purely local. 
diff --git a/mlprec/impl/solver/mld_z_gs_solver_apply_vect.f90 b/mlprec/impl/solver/mld_z_gs_solver_apply_vect.f90
index 5174ddc9..905c2c66 100644
--- a/mlprec/impl/solver/mld_z_gs_solver_apply_vect.f90
+++ b/mlprec/impl/solver/mld_z_gs_solver_apply_vect.f90
@@ -53,7 +53,7 @@ subroutine mld_z_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
   character, intent(in), optional                :: init
   type(psb_z_vect_type),intent(inout), optional   :: initu
 
-  integer(psb_ipk_)   :: n_row,n_col, itx
+  integer(psb_ipk_)   :: n_row,n_col, itx, itxst
   complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   complex(psb_dpk_), allocatable :: temp(:)
   integer(psb_ipk_)   :: ictxt,np,me,i, err_act
@@ -128,10 +128,12 @@ subroutine mld_z_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
   end if
 
   associate(tw => wv(1), xit => wv(2))
-
+    itxst = 1 
     select case (init_)
     case('Z') 
-      call xit%zero()
+      call psb_geaxpby(zone,x,zzero,tw,desc_data,info)
+      call psb_spsm(zone,sv%l,tw,zzero,xit,desc_data,info)
+      itxst = 2
     case('Y')
       call psb_geaxpby(zone,y,zzero,xit,desc_data,info)
     case('U')
@@ -154,7 +156,7 @@ subroutine mld_z_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
         ! Fixed number of iterations
         !
         !
-        do itx=1,sv%sweeps
+        do itx=itxst,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.