Streamline GS first application step when INIT='Z'

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
@@ -119,9 +119,12 @@ subroutine mld_c_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,&
 
   call psb_geasb(wv,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. 

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. 

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
@@ -119,9 +119,12 @@ subroutine mld_c_gs_solver_apply(alpha,sv,x,beta,y,desc_data,&
 
   call psb_geasb(wv,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. 

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. 

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
@@ -119,9 +119,12 @@ subroutine mld_d_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,&
 
   call psb_geasb(wv,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. 

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. 

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
@@ -119,9 +119,12 @@ subroutine mld_d_gs_solver_apply(alpha,sv,x,beta,y,desc_data,&
 
   call psb_geasb(wv,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. 

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. 

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
@@ -119,9 +119,12 @@ subroutine mld_s_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,&
 
   call psb_geasb(wv,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. 

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. 

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
@@ -119,9 +119,12 @@ subroutine mld_s_gs_solver_apply(alpha,sv,x,beta,y,desc_data,&
 
   call psb_geasb(wv,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. 

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. 

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
@@ -119,9 +119,12 @@ subroutine mld_z_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,&
 
   call psb_geasb(wv,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. 

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. 

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
@@ -119,9 +119,12 @@ subroutine mld_z_gs_solver_apply(alpha,sv,x,beta,y,desc_data,&
 
   call psb_geasb(wv,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. 

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.