diff --git a/linsolve/Makefile b/linsolve/Makefile
index fa64b1b7e..40ad1e5bf 100644
--- a/linsolve/Makefile
+++ b/linsolve/Makefile
@@ -9,6 +9,8 @@ MODDIR=../modules
 MODOBJS= psb_base_linsolve_conv_mod.o \
 	psb_s_linsolve_conv_mod.o psb_c_linsolve_conv_mod.o \
 	psb_d_linsolve_conv_mod.o psb_z_linsolve_conv_mod.o \
+	psb_s_linsolve_mod.o psb_c_linsolve_mod.o \
+	psb_d_linsolve_mod.o psb_z_linsolve_mod.o \
 	psb_linsolve_mod.o
 
 OBJS=$(MODOBJS)
@@ -37,6 +39,9 @@ impld: $(OBJS)
 
 psb_s_linsolve_conv_mod.o psb_c_linsolve_conv_mod.o psb_d_linsolve_conv_mod.o psb_z_linsolve_conv_mod.o:  psb_base_linsolve_conv_mod.o
 psb_linsolve_conv_mod.o: psb_s_linsolve_conv_mod.o psb_c_linsolve_conv_mod.o psb_d_linsolve_conv_mod.o psb_z_linsolve_conv_mod.o 
+psb_linsolve_mod.o: psb_s_linsolve_mod.o psb_c_linsolve_mod.o psb_d_linsolve_mod.o psb_z_linsolve_mod.o 
+
+
 $(F90OBJS): $(MODOBJS)
 $(OBJS): $(MODDIR)/$(PRECMODNAME)$(.mod)  $(MODDIR)/$(BASEMODNAME)$(.mod)
 
diff --git a/linsolve/impl/psb_cbicg.f90 b/linsolve/impl/psb_cbicg.f90
index 246dcc10f..30d150aa6 100644
--- a/linsolve/impl/psb_cbicg.f90
+++ b/linsolve/impl/psb_cbicg.f90
@@ -95,7 +95,7 @@
 !
 
 subroutine psb_cbicg_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,istop)
+     & itmax,iter,err,itrace,istop,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_c_linsolve_conv_mod
@@ -111,6 +111,7 @@ subroutine psb_cbicg_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), optional, intent(in)        :: itmax, itrace, istop
   integer(psb_ipk_), optional, intent(out)       :: iter
   real(psb_spk_), optional, intent(out) :: err
+  type(psb_c_vect_type), intent(inout), optional   :: s1, s2
 ! !$   local data
   complex(psb_spk_), allocatable, target   :: aux(:)
   type(psb_c_vect_type), allocatable, target :: wwrk(:)
@@ -236,7 +237,8 @@ subroutine psb_cbicg_vect(a,prec,b,x,eps,desc_a,info,&
   itx   = 0
 
 
-  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
+  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
+       & desc_a,stopdat,info,s1=s1,s2=s2)
   if (info /= psb_success_) Then 
      call psb_errpush(psb_err_from_subroutine_non_,name)
      goto 9999
@@ -262,7 +264,7 @@ subroutine psb_cbicg_vect(a,prec,b,x,eps,desc_a,info,&
     rho = czero
     
     ! Perhaps we already satisfy the convergence criterion...
-    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
     if (info /= psb_success_) Then 
       call psb_errpush(psb_err_from_subroutine_non_,name)
       goto 9999
@@ -316,7 +318,7 @@ subroutine psb_cbicg_vect(a,prec,b,x,eps,desc_a,info,&
       call psb_geaxpby(-alpha,q,cone,r,desc_a,info)
       call psb_geaxpby(-alpha,qt,cone,rt,desc_a,info)
 
-      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
       if (info /= psb_success_) Then 
         call psb_errpush(psb_err_from_subroutine_non_,name)
         goto 9999
diff --git a/linsolve/impl/psb_ccg.F90 b/linsolve/impl/psb_ccg.F90
index 399246254..cf1ea1bea 100644
--- a/linsolve/impl/psb_ccg.F90
+++ b/linsolve/impl/psb_ccg.F90
@@ -96,7 +96,7 @@
 ! 
 !
 subroutine psb_ccg_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,istop,cond)
+     & itmax,iter,err,itrace,istop,cond,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_c_linsolve_conv_mod
@@ -112,6 +112,8 @@ subroutine psb_ccg_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   Real(psb_spk_), Optional, Intent(out) :: err,cond
+  type(psb_c_vect_type), intent(inout), optional   :: s1, s2
+      
 ! =   Local data
   complex(psb_spk_), allocatable, target   :: aux(:),td(:),tu(:),eig(:),ewrk(:)
   integer(psb_mpk_), allocatable :: ibl(:), ispl(:), iwrk(:)
@@ -245,7 +247,8 @@ subroutine psb_ccg_vect(a,prec,b,x,eps,desc_a,info,&
 
     rho = czero
     
-    call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
+    call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
+         &desc_a,stopdat,info,s1=s1,s2=s2)
     if (info /= psb_success_) Then 
       call psb_errpush(psb_err_from_subroutine_non_,name)
       goto 9999
@@ -289,7 +292,7 @@ subroutine psb_ccg_vect(a,prec,b,x,eps,desc_a,info,&
       call psb_geaxpby(alpha,p,cone,x,desc_a,info)
       call psb_geaxpby(-alpha,q,cone,r,desc_a,info)
 
-      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
       if (info /= psb_success_) Then 
         call psb_errpush(psb_err_from_subroutine_non_,name)
         goto 9999
diff --git a/linsolve/impl/psb_ccgs.f90 b/linsolve/impl/psb_ccgs.f90
index bb4ebd182..8783a0dad 100644
--- a/linsolve/impl/psb_ccgs.f90
+++ b/linsolve/impl/psb_ccgs.f90
@@ -93,7 +93,7 @@
 !                                         estimate of) residual. 
 !
 Subroutine psb_ccgs_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,istop)
+     & itmax,iter,err,itrace,istop,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_c_linsolve_conv_mod
@@ -109,6 +109,7 @@ Subroutine psb_ccgs_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace,istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   Real(psb_spk_), Optional, Intent(out) :: err
+  type(psb_c_vect_type), intent(inout), optional   :: s1, s2
 ! =   local data
   complex(psb_spk_), allocatable, target   :: aux(:)
   type(psb_c_vect_type), allocatable, target :: wwrk(:)
@@ -223,7 +224,8 @@ Subroutine psb_ccgs_vect(a,prec,b,x,eps,desc_a,info,&
 
   itx   = 0
 
-  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
+  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
+       & desc_a,stopdat,info,s1=s1,s2=s2)
   if (info /= psb_success_) Then 
      call psb_errpush(psb_err_from_subroutine_non_,name)
      goto 9999
@@ -246,7 +248,7 @@ Subroutine psb_ccgs_vect(a,prec,b,x,eps,desc_a,info,&
     
 
     ! Perhaps we already satisfy the convergence criterion...
-    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
     if (info /= psb_success_) Then 
       call psb_errpush(psb_err_from_subroutine_non_,name)
       goto 9999
@@ -320,7 +322,7 @@ Subroutine psb_ccgs_vect(a,prec,b,x,eps,desc_a,info,&
          goto 9999
       end if
 
-      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
       if (info /= psb_success_) Then 
         call psb_errpush(psb_err_from_subroutine_non_,name)
         goto 9999
diff --git a/linsolve/impl/psb_ccgstab.f90 b/linsolve/impl/psb_ccgstab.f90
index 4d5678a91..20c0a22a7 100644
--- a/linsolve/impl/psb_ccgstab.f90
+++ b/linsolve/impl/psb_ccgstab.f90
@@ -93,7 +93,7 @@
 !                                         where r is the (preconditioned, recursive
 !                                         estimate of) residual. 
 !
-Subroutine psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
+Subroutine psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_c_linsolve_conv_mod
@@ -109,6 +109,7 @@ Subroutine psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
   integer(psb_ipk_), Optional, Intent(in)      :: itmax, itrace, istop
   integer(psb_ipk_), Optional, Intent(out)     :: iter
   Real(psb_spk_), Optional, Intent(out) :: err
+  type(psb_c_vect_type), intent(inout), optional   :: s1, s2
 ! =   Local data
   complex(psb_spk_), allocatable, target   :: aux(:),wwrk(:,:)
   type(psb_c_vect_type) :: q, r, p, v, s, t, z, f
@@ -235,7 +236,8 @@ Subroutine psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
   End If
 
   itx   = 0
-  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
+  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
+       & desc_a,stopdat,info,s1=s1,s2=s2)
   if (psb_errstatus_fatal()) Then 
     call psb_errpush(psb_err_from_subroutine_non_,name)
     goto 9999
@@ -252,7 +254,7 @@ Subroutine psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
     call psb_geaxpby(cone,r,czero,q,desc_a,info)
 
     ! Perhaps we already satisfy the convergence criterion...
-    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
 
     if (psb_errstatus_fatal()) then 
       info=psb_err_from_subroutine_
@@ -372,7 +374,7 @@ Subroutine psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
       call psb_geaxpby(omega,z,cone,x,desc_a,info)
       call psb_geaxpby(cone,s,czero,r,desc_a,info)
       call psb_geaxpby(-omega,t,cone,r,desc_a,info)
-      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
 
       if (psb_errstatus_fatal()) Then 
         call psb_errpush(psb_err_from_subroutine_,name,a_err='X/R update ')
diff --git a/linsolve/impl/psb_ccgstabl.f90 b/linsolve/impl/psb_ccgstabl.f90
index aa81e00fd..b8898c3aa 100644
--- a/linsolve/impl/psb_ccgstabl.f90
+++ b/linsolve/impl/psb_ccgstabl.f90
@@ -104,7 +104,7 @@
 !
 !
 Subroutine psb_ccgstabl_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,irst,istop)
+     & itmax,iter,err,itrace,irst,istop,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_c_linsolve_conv_mod
@@ -120,6 +120,7 @@ Subroutine psb_ccgstabl_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, irst,istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   Real(psb_spk_), Optional, Intent(out) :: err
+  type(psb_c_vect_type), intent(inout), optional   :: s1, s2
 ! =   local data
   complex(psb_spk_), allocatable, target   :: aux(:), gamma(:),&
        & gamma1(:), gamma2(:), taum(:,:), sigma(:)
@@ -267,7 +268,8 @@ Subroutine psb_ccgstabl_vect(a,prec,b,x,eps,desc_a,info,&
   rt0 => wwrk(10)
   
 
-  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
+  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
+       & desc_a,stopdat,info,s1=s1,s2=s2)
   if (info /= psb_success_) Then 
      call psb_errpush(psb_err_from_subroutine_non_,name)
      goto 9999
@@ -305,7 +307,7 @@ Subroutine psb_ccgstabl_vect(a,prec,b,x,eps,desc_a,info,&
          & write(debug_unit,*) me,' ',trim(name),&
          & ' on entry to amax: b: ',b%get_nrows()
 
-    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
     if (info /= psb_success_) Then 
       call psb_errpush(psb_err_from_subroutine_non_,name)
       goto 9999
@@ -409,7 +411,7 @@ Subroutine psb_ccgstabl_vect(a,prec,b,x,eps,desc_a,info,&
         call psb_geaxpby(-gamma1(j),rh(j),cone,rh(0),desc_a,info)        
       enddo
       
-      if (psb_check_conv(methdname,itx,x,rh(0),desc_a,stopdat,info)) exit restart
+      if (psb_check_conv(methdname,itx,x,rh(0),desc_a,stopdat,info,s1=s1)) exit restart
       if (info /= psb_success_) Then 
         call psb_errpush(psb_err_from_subroutine_non_,name)
         goto 9999
diff --git a/linsolve/impl/psb_cfcg.F90 b/linsolve/impl/psb_cfcg.F90
index d5db7f947..252e19e8e 100644
--- a/linsolve/impl/psb_cfcg.F90
+++ b/linsolve/impl/psb_cfcg.F90
@@ -104,7 +104,7 @@
 ! 
 !
 subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,istop,cond)
+     & itmax,iter,err,itrace,istop,cond,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_c_linsolve_conv_mod
@@ -120,6 +120,7 @@ subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   real(psb_spk_), Optional, Intent(out) :: err,cond
+  type(psb_c_vect_type), intent(inout), optional   :: s1, s2
 ! =   Local data
   type(psb_c_vect_type)  :: v, w, d , q, r
   complex(psb_spk_) :: alpha, beta, delta, gamma, theta
@@ -227,7 +228,7 @@ subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,&
        & scratch=.true.,mold=x%v)
 
   call psb_init_conv(methdname,istop_,itrace_,itmax_,&
-       & a,x,b,eps,desc_a,stopdat,info)
+       & a,x,b,eps,desc_a,stopdat,info,s1=s1,s2=s2)
   itx = 0 
 
   restart: do 
@@ -246,7 +247,7 @@ subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,&
     end if
     
 
-    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) then
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) then
       if (debug.and.(me==0)) write(0,*) name,' Exit on  convergence from restart'
       exit restart
     end if
@@ -302,7 +303,7 @@ subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,&
 
       itx = itx + 1
 
-      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) then
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) then
         if (debug.and.(me==0)) write(0,*) name,' Exit on  convergence from iteration'
         exit restart
       end if
diff --git a/linsolve/impl/psb_cgcr.f90 b/linsolve/impl/psb_cgcr.f90
index 59129e161..2552003c6 100644
--- a/linsolve/impl/psb_cgcr.f90
+++ b/linsolve/impl/psb_cgcr.f90
@@ -106,7 +106,7 @@
 !
 
 subroutine psb_cgcr_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace, irst, istop)
+     & itmax,iter,err,itrace, irst, istop,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_c_linsolve_conv_mod
@@ -124,6 +124,7 @@ subroutine psb_cgcr_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, irst, istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   real(psb_spk_), Optional, Intent(out) :: err
+  type(psb_c_vect_type), intent(inout), optional   :: s1, s2
   ! =   local data
   complex(psb_spk_), allocatable   :: alpha(:), h(:,:)
   type(psb_c_vect_type), allocatable :: z(:), c(:), c_scale(:)
@@ -253,7 +254,8 @@ subroutine psb_cgcr_vect(a,prec,b,x,eps,desc_a,info,&
   itx = 0
   
   nrst = -1
-  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
+  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
+       & desc_a,stopdat,info,s1=s1,s2=s2)
   restart: do 
     if (itx>= itmax_) exit restart 
     h = czero
@@ -276,7 +278,7 @@ subroutine psb_cgcr_vect(a,prec,b,x,eps,desc_a,info,&
       goto 9999
     end if
     
-    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
         
     nrst = nrst + 1 
     
@@ -307,7 +309,7 @@ subroutine psb_cgcr_vect(a,prec,b,x,eps,desc_a,info,&
       call psb_geaxpby(cone, r, czero, r, desc_a, info)   
       call psb_geaxpby(-alpha(j), c_scale(j), cone, r, desc_a, info)   
       
-      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
       
       if (j >= irst) exit iteration
       
diff --git a/linsolve/impl/psb_ckrylov.f90 b/linsolve/impl/psb_ckrylov.f90
index 308e62aea..5f708ed56 100644
--- a/linsolve/impl/psb_ckrylov.f90
+++ b/linsolve/impl/psb_ckrylov.f90
@@ -80,7 +80,7 @@
 !                                           estimate of) residual 
 ! 
 Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,irst,istop,cond)
+     & itmax,iter,err,itrace,irst,istop,cond,s1,s2)
 
   use psb_base_mod
   use psb_prec_mod,only : psb_cprec_type
@@ -97,11 +97,12 @@ Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, irst,istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   Real(psb_spk_), Optional, Intent(out) :: err,cond
+  type(psb_c_vect_type), intent(inout), optional   :: s1, s2
 
 
   abstract interface
     subroutine psb_ckryl_vect(a,prec,b,x,eps,&
-         & desc_a,info,itmax,iter,err,itrace,istop)
+         & desc_a,info,itmax,iter,err,itrace,istop,s1,s2)
       import :: psb_ipk_, psb_spk_, psb_desc_type, &
            & psb_cspmat_type, psb_cprec_type, psb_c_vect_type
       type(psb_cspmat_type), intent(in)    :: a
@@ -114,9 +115,10 @@ Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
       integer(psb_ipk_), optional, intent(in)        :: itmax, itrace,istop
       integer(psb_ipk_), optional, intent(out)       :: iter
       real(psb_spk_), optional, intent(out) :: err
+      type(psb_c_vect_type), intent(inout), optional   :: s1, s2
     end subroutine psb_ckryl_vect
     Subroutine psb_ckryl_rest_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err, itrace,irst,istop)
+         &itmax,iter,err, itrace,irst,istop,s1,s2)
       import :: psb_ipk_, psb_spk_, psb_desc_type, &
            & psb_cspmat_type, psb_cprec_type, psb_c_vect_type
       Type(psb_cspmat_type), Intent(in)    :: a
@@ -129,9 +131,10 @@ Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
       integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, irst,istop
       integer(psb_ipk_), Optional, Intent(out)       :: iter
       Real(psb_spk_), Optional, Intent(out) :: err
+      type(psb_c_vect_type), intent(inout), optional   :: s1, s2
     end subroutine psb_ckryl_rest_vect
     Subroutine psb_ckryl_cond_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err, itrace,istop,cond)
+         &itmax,iter,err, itrace,istop,cond,s1,s2)
       import :: psb_ipk_, psb_spk_, psb_desc_type, &
            & psb_cspmat_type, psb_cprec_type, psb_c_vect_type
       Type(psb_cspmat_type), Intent(in)    :: a
@@ -144,6 +147,7 @@ Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
       integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace,istop
       integer(psb_ipk_), Optional, Intent(out)       :: iter
       Real(psb_spk_), Optional, Intent(out) :: err, cond
+      type(psb_c_vect_type), intent(inout), optional   :: s1, s2
     end subroutine psb_ckryl_cond_vect
   end interface
 
@@ -180,37 +184,37 @@ Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
   select case(psb_toupper(method))
   case('CG') 
     call  psb_ccg_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond)
+         &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond,s1=s1,s2=s2)
   case('FCG') 
     call  psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond)
+         &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond,s1=s1,s2=s2)
   case('GCR') 
     call  psb_cgcr_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop)
+         &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
   case('CGS') 
     call  psb_ccgs_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop)
+         &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
   case('BICG') 
     call  psb_cbicg_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop)
+         &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
   case('BICGSTAB') 
     call  psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop)
+         &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
   case('RGMRES','GMRES')
     call  psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop)
+         &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop,s1=s1,s2=s2)
   case('MINRES','PMINRES')
     call  psb_cminres_vect(a,prec,b,x,eps,desc_a,info,&
          &itmax,iter,err,itrace=itrace_,istop=istop)
   case('BICGSTABL')
     call  psb_ccgstabl_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop)
+         &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop,s1=s1,s2=s2)
   case default
     if (me == 0) write(psb_err_unit,*) trim(name),&
          & ': Warning: Unknown method  ',method,&
          & ', defaulting to BiCGSTAB'
     call  psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop)
+         &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
   end select
 
   if ((info==psb_success_).and.do_alloc_wrk) call prec%free_wrk(info)
diff --git a/linsolve/impl/psb_crgmres.f90 b/linsolve/impl/psb_crgmres.f90
index ab32a8938..90f8588d1 100644
--- a/linsolve/impl/psb_crgmres.f90
+++ b/linsolve/impl/psb_crgmres.f90
@@ -102,13 +102,15 @@
 !                                            stopped when  |r| <= eps * (|a||x|+|b|)
 !                                         2: err =  |r|/|b|; here the iteration is
 !                                            stopped when  |r| <= eps * |b|
+!                                         3: Same as 2 but with X and B scaled
+!                                            by s1 and s2 
 !                                         where r is the (preconditioned, recursive
 !                                         estimate of) residual. 
 !    irst   -  integer(optional)          Input: restart parameter 
 !
 
 subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,irst,istop)
+     & itmax,iter,err,itrace,irst,istop,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_c_linsolve_conv_mod
@@ -124,6 +126,7 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, irst,istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   Real(psb_spk_), Optional, Intent(out) :: err
+  type(psb_c_vect_type), intent(inout), optional   :: s1, s2
 ! =   local data
   complex(psb_spk_), allocatable   :: aux(:)
   complex(psb_spk_), allocatable   :: c(:), s(:), h(:,:), rs(:), rst(:)
@@ -268,9 +271,20 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
   select case(istop_)
   case(psb_istop_ani_)
     ani = psb_spnrmi(a,desc_a,info)
-    bni = psb_geamax(b,desc_a,info)
+    if (present(s1)) then
+      call psb_gemlt(cone,s1,b,czero,v(1),desc_a,info)
+      bni = psb_geamax(v(1),desc_a,info)
+    else
+      bni = psb_geamax(b,desc_a,info)
+    end if
   case(psb_istop_bn2_)
-    bn2 = psb_genrm2(b,desc_a,info)    
+    if (present(s1)) then
+      call psb_gemlt(cone,s1,b,czero,v(1),desc_a,info)
+      bn2 = psb_genrm2(v(1),desc_a,info)
+    else
+      bn2 = psb_genrm2(b,desc_a,info)
+    end if
+      
   case(psb_istop_rn2_abs_)
     ! do nothing
   case(psb_istop_rrn2_)
@@ -282,6 +296,7 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
     end if
     
     call psb_spmm(-cone,a,x,cone,v(1),desc_a,info,work=aux)
+    if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info)
     if (info /= psb_success_) then 
       info=psb_err_from_subroutine_non_ 
       call psb_errpush(info,name)
@@ -323,7 +338,8 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
       call psb_errpush(info,name)
       goto 9999
     end if
-
+    if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info)
+    
     rs(1) = psb_genrm2(v(1),desc_a,info)
     rs(2:) = czero
     if (info /= psb_success_) then 
@@ -378,8 +394,14 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
     inner:  Do i=1,nl
       itx  = itx + 1
 
-      call prec%apply(v(i),w1,desc_a,info)
+      if (present(s2)) then
+        call psb_gediv(v(i),s2,w,desc_a,info)
+        call prec%apply(w,w1,desc_a,info)
+      else
+        call prec%apply(v(i),w1,desc_a,info)
+      end if
       call psb_spmm(cone,a,w1,czero,w,desc_a,info,work=aux)
+      if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
       !
       call mgs(i,h,v,w,rs,c,s,desc_a,info)
 
@@ -391,10 +413,11 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
         !
         rst = rs
         call psb_geaxpby(cone,x,czero,xt,desc_a,info)
-        call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info)
+        call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info,s2=s2)
 
         call psb_geaxpby(cone,b,czero,w1,desc_a,info)
         call psb_spmm(-cone,a,xt,cone,w1,desc_a,info,work=aux)
+        if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
         rni = psb_geamax(w1,desc_a,info)
         xni = psb_geamax(xt,desc_a,info)
         errnum = rni
@@ -432,7 +455,8 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
           call psb_geaxpby(cone,xt,czero,x,desc_a,info)
 ! =          x = xt 
         case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
-          call  rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info)          !
+          call  rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2)
+          !
 
         end select
 
@@ -452,7 +476,7 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
       call psb_geaxpby(cone,xt,czero,x,desc_a,info)!      x = xt 
 
     case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
-      call  rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info)          !
+      call  rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2)          !
     end select
     
     if (itx >= itmax_) then 
@@ -523,11 +547,12 @@ contains
   ! Rebuild solution X from the space V using the factor
   ! stored in R
   !
-  subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info)
+  subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2)
     complex(psb_spk_)   :: c(:), s(:), rs(:), h(:,:)
     type(psb_c_vect_type) :: v(:), w, w1, x
     type(psb_desc_type) :: desc_a
     class(psb_cprec_type) :: prec
+    type(psb_c_vect_type), intent(inout), optional   :: s2
     integer(psb_ipk_) :: info
     integer(psb_ipk_) :: k,n
 
@@ -539,12 +564,13 @@ contains
     if (debug_level >= psb_debug_ext_) &
          & write(debug_unit,*) me,' ',trim(name),&
          & ' Rebuild x-> RS:',rs(1:n)
-    call w1%zero()
+    call w%zero()
     do k=1, n
-      call psb_geaxpby(rs(k),v(k),cone,w1,desc_a,info)
+      call psb_geaxpby(rs(k),v(k),cone,w,desc_a,info)
     end do
-    call prec%apply(w1,w,desc_a,info)
-    call psb_geaxpby(cone,w,cone,x,desc_a,info)
+    if (present(s2)) call psb_gediv(s2,w,desc_a,info)
+    call prec%apply(w,w1,desc_a,info)
+    call psb_geaxpby(cone,w1,cone,x,desc_a,info)
   end subroutine rebuildx
 
 end subroutine psb_crgmres_vect
diff --git a/linsolve/impl/psb_dbicg.f90 b/linsolve/impl/psb_dbicg.f90
index f56ae6d3d..dc8018800 100644
--- a/linsolve/impl/psb_dbicg.f90
+++ b/linsolve/impl/psb_dbicg.f90
@@ -95,7 +95,7 @@
 !
 
 subroutine psb_dbicg_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,istop)
+     & itmax,iter,err,itrace,istop,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_d_linsolve_conv_mod
@@ -111,6 +111,7 @@ subroutine psb_dbicg_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), optional, intent(in)        :: itmax, itrace, istop
   integer(psb_ipk_), optional, intent(out)       :: iter
   real(psb_dpk_), optional, intent(out) :: err
+  type(psb_d_vect_type), intent(inout), optional   :: s1, s2
 ! !$   local data
   real(psb_dpk_), allocatable, target   :: aux(:)
   type(psb_d_vect_type), allocatable, target :: wwrk(:)
@@ -236,7 +237,8 @@ subroutine psb_dbicg_vect(a,prec,b,x,eps,desc_a,info,&
   itx   = 0
 
 
-  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
+  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
+       & desc_a,stopdat,info,s1=s1,s2=s2)
   if (info /= psb_success_) Then 
      call psb_errpush(psb_err_from_subroutine_non_,name)
      goto 9999
@@ -262,7 +264,7 @@ subroutine psb_dbicg_vect(a,prec,b,x,eps,desc_a,info,&
     rho = dzero
     
     ! Perhaps we already satisfy the convergence criterion...
-    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
     if (info /= psb_success_) Then 
       call psb_errpush(psb_err_from_subroutine_non_,name)
       goto 9999
@@ -316,7 +318,7 @@ subroutine psb_dbicg_vect(a,prec,b,x,eps,desc_a,info,&
       call psb_geaxpby(-alpha,q,done,r,desc_a,info)
       call psb_geaxpby(-alpha,qt,done,rt,desc_a,info)
 
-      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
       if (info /= psb_success_) Then 
         call psb_errpush(psb_err_from_subroutine_non_,name)
         goto 9999
diff --git a/linsolve/impl/psb_dcg.F90 b/linsolve/impl/psb_dcg.F90
index 767a4162b..a60742b58 100644
--- a/linsolve/impl/psb_dcg.F90
+++ b/linsolve/impl/psb_dcg.F90
@@ -96,7 +96,7 @@
 ! 
 !
 subroutine psb_dcg_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,istop,cond)
+     & itmax,iter,err,itrace,istop,cond,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_d_linsolve_conv_mod
@@ -112,6 +112,8 @@ subroutine psb_dcg_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   Real(psb_dpk_), Optional, Intent(out) :: err,cond
+  type(psb_d_vect_type), intent(inout), optional   :: s1, s2
+      
 ! =   Local data
   real(psb_dpk_), allocatable, target   :: aux(:),td(:),tu(:),eig(:),ewrk(:)
   integer(psb_mpk_), allocatable :: ibl(:), ispl(:), iwrk(:)
@@ -253,7 +255,8 @@ subroutine psb_dcg_vect(a,prec,b,x,eps,desc_a,info,&
 
     rho = dzero
     
-    call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
+    call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
+         &desc_a,stopdat,info,s1=s1,s2=s2)
     if (info /= psb_success_) Then 
       call psb_errpush(psb_err_from_subroutine_non_,name)
       goto 9999
@@ -306,7 +309,7 @@ subroutine psb_dcg_vect(a,prec,b,x,eps,desc_a,info,&
       call psb_geaxpby(alpha,p,done,x,desc_a,info)
       call psb_geaxpby(-alpha,q,done,r,desc_a,info)
 
-      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
       if (info /= psb_success_) Then 
         call psb_errpush(psb_err_from_subroutine_non_,name)
         goto 9999
diff --git a/linsolve/impl/psb_dcgs.f90 b/linsolve/impl/psb_dcgs.f90
index 85b1e73ba..60173a057 100644
--- a/linsolve/impl/psb_dcgs.f90
+++ b/linsolve/impl/psb_dcgs.f90
@@ -93,7 +93,7 @@
 !                                         estimate of) residual. 
 !
 Subroutine psb_dcgs_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,istop)
+     & itmax,iter,err,itrace,istop,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_d_linsolve_conv_mod
@@ -109,6 +109,7 @@ Subroutine psb_dcgs_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace,istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   Real(psb_dpk_), Optional, Intent(out) :: err
+  type(psb_d_vect_type), intent(inout), optional   :: s1, s2
 ! =   local data
   real(psb_dpk_), allocatable, target   :: aux(:)
   type(psb_d_vect_type), allocatable, target :: wwrk(:)
@@ -223,7 +224,8 @@ Subroutine psb_dcgs_vect(a,prec,b,x,eps,desc_a,info,&
 
   itx   = 0
 
-  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
+  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
+       & desc_a,stopdat,info,s1=s1,s2=s2)
   if (info /= psb_success_) Then 
      call psb_errpush(psb_err_from_subroutine_non_,name)
      goto 9999
@@ -246,7 +248,7 @@ Subroutine psb_dcgs_vect(a,prec,b,x,eps,desc_a,info,&
     
 
     ! Perhaps we already satisfy the convergence criterion...
-    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
     if (info /= psb_success_) Then 
       call psb_errpush(psb_err_from_subroutine_non_,name)
       goto 9999
@@ -320,7 +322,7 @@ Subroutine psb_dcgs_vect(a,prec,b,x,eps,desc_a,info,&
          goto 9999
       end if
 
-      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
       if (info /= psb_success_) Then 
         call psb_errpush(psb_err_from_subroutine_non_,name)
         goto 9999
diff --git a/linsolve/impl/psb_dcgstab.f90 b/linsolve/impl/psb_dcgstab.f90
index 65e20a760..95633379b 100644
--- a/linsolve/impl/psb_dcgstab.f90
+++ b/linsolve/impl/psb_dcgstab.f90
@@ -93,7 +93,7 @@
 !                                         where r is the (preconditioned, recursive
 !                                         estimate of) residual. 
 !
-Subroutine psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
+Subroutine psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_d_linsolve_conv_mod
@@ -109,6 +109,7 @@ Subroutine psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
   integer(psb_ipk_), Optional, Intent(in)      :: itmax, itrace, istop
   integer(psb_ipk_), Optional, Intent(out)     :: iter
   Real(psb_dpk_), Optional, Intent(out) :: err
+  type(psb_d_vect_type), intent(inout), optional   :: s1, s2
 ! =   Local data
   real(psb_dpk_), allocatable, target   :: aux(:),wwrk(:,:)
   type(psb_d_vect_type) :: q, r, p, v, s, t, z, f
@@ -235,7 +236,8 @@ Subroutine psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
   End If
 
   itx   = 0
-  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
+  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
+       & desc_a,stopdat,info,s1=s1,s2=s2)
   if (psb_errstatus_fatal()) Then 
     call psb_errpush(psb_err_from_subroutine_non_,name)
     goto 9999
@@ -252,7 +254,7 @@ Subroutine psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
     call psb_geaxpby(done,r,dzero,q,desc_a,info)
 
     ! Perhaps we already satisfy the convergence criterion...
-    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
 
     if (psb_errstatus_fatal()) then 
       info=psb_err_from_subroutine_
@@ -372,7 +374,7 @@ Subroutine psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
       call psb_geaxpby(omega,z,done,x,desc_a,info)
       call psb_geaxpby(done,s,dzero,r,desc_a,info)
       call psb_geaxpby(-omega,t,done,r,desc_a,info)
-      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
 
       if (psb_errstatus_fatal()) Then 
         call psb_errpush(psb_err_from_subroutine_,name,a_err='X/R update ')
diff --git a/linsolve/impl/psb_dcgstabl.f90 b/linsolve/impl/psb_dcgstabl.f90
index 893ae01df..152dffdc8 100644
--- a/linsolve/impl/psb_dcgstabl.f90
+++ b/linsolve/impl/psb_dcgstabl.f90
@@ -104,7 +104,7 @@
 !
 !
 Subroutine psb_dcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,irst,istop)
+     & itmax,iter,err,itrace,irst,istop,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_d_linsolve_conv_mod
@@ -120,6 +120,7 @@ Subroutine psb_dcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, irst,istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   Real(psb_dpk_), Optional, Intent(out) :: err
+  type(psb_d_vect_type), intent(inout), optional   :: s1, s2
 ! =   local data
   real(psb_dpk_), allocatable, target   :: aux(:), gamma(:),&
        & gamma1(:), gamma2(:), taum(:,:), sigma(:)
@@ -267,7 +268,8 @@ Subroutine psb_dcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
   rt0 => wwrk(10)
   
 
-  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
+  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
+       & desc_a,stopdat,info,s1=s1,s2=s2)
   if (info /= psb_success_) Then 
      call psb_errpush(psb_err_from_subroutine_non_,name)
      goto 9999
@@ -305,7 +307,7 @@ Subroutine psb_dcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
          & write(debug_unit,*) me,' ',trim(name),&
          & ' on entry to amax: b: ',b%get_nrows()
 
-    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
     if (info /= psb_success_) Then 
       call psb_errpush(psb_err_from_subroutine_non_,name)
       goto 9999
@@ -409,7 +411,7 @@ Subroutine psb_dcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
         call psb_geaxpby(-gamma1(j),rh(j),done,rh(0),desc_a,info)        
       enddo
       
-      if (psb_check_conv(methdname,itx,x,rh(0),desc_a,stopdat,info)) exit restart
+      if (psb_check_conv(methdname,itx,x,rh(0),desc_a,stopdat,info,s1=s1)) exit restart
       if (info /= psb_success_) Then 
         call psb_errpush(psb_err_from_subroutine_non_,name)
         goto 9999
diff --git a/linsolve/impl/psb_dfcg.F90 b/linsolve/impl/psb_dfcg.F90
index 0885ac3f3..9e0d62213 100644
--- a/linsolve/impl/psb_dfcg.F90
+++ b/linsolve/impl/psb_dfcg.F90
@@ -104,7 +104,7 @@
 ! 
 !
 subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,istop,cond)
+     & itmax,iter,err,itrace,istop,cond,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_d_linsolve_conv_mod
@@ -120,6 +120,7 @@ subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   real(psb_dpk_), Optional, Intent(out) :: err,cond
+  type(psb_d_vect_type), intent(inout), optional   :: s1, s2
 ! =   Local data
   type(psb_d_vect_type)  :: v, w, d , q, r
   real(psb_dpk_) :: alpha, beta, delta, gamma, theta
@@ -227,7 +228,7 @@ subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,&
        & scratch=.true.,mold=x%v)
 
   call psb_init_conv(methdname,istop_,itrace_,itmax_,&
-       & a,x,b,eps,desc_a,stopdat,info)
+       & a,x,b,eps,desc_a,stopdat,info,s1=s1,s2=s2)
   itx = 0 
 
   restart: do 
@@ -246,7 +247,7 @@ subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,&
     end if
     
 
-    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) then
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) then
       if (debug.and.(me==0)) write(0,*) name,' Exit on  convergence from restart'
       exit restart
     end if
@@ -302,7 +303,7 @@ subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,&
 
       itx = itx + 1
 
-      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) then
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) then
         if (debug.and.(me==0)) write(0,*) name,' Exit on  convergence from iteration'
         exit restart
       end if
diff --git a/linsolve/impl/psb_dgcr.f90 b/linsolve/impl/psb_dgcr.f90
index f26503ef7..43400a9aa 100644
--- a/linsolve/impl/psb_dgcr.f90
+++ b/linsolve/impl/psb_dgcr.f90
@@ -106,7 +106,7 @@
 !
 
 subroutine psb_dgcr_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace, irst, istop)
+     & itmax,iter,err,itrace, irst, istop,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_d_linsolve_conv_mod
@@ -124,6 +124,7 @@ subroutine psb_dgcr_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, irst, istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   real(psb_dpk_), Optional, Intent(out) :: err
+  type(psb_d_vect_type), intent(inout), optional   :: s1, s2
   ! =   local data
   real(psb_dpk_), allocatable   :: alpha(:), h(:,:)
   type(psb_d_vect_type), allocatable :: z(:), c(:), c_scale(:)
@@ -253,7 +254,8 @@ subroutine psb_dgcr_vect(a,prec,b,x,eps,desc_a,info,&
   itx = 0
   
   nrst = -1
-  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
+  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
+       & desc_a,stopdat,info,s1=s1,s2=s2)
   restart: do 
     if (itx>= itmax_) exit restart 
     h = dzero
@@ -276,7 +278,7 @@ subroutine psb_dgcr_vect(a,prec,b,x,eps,desc_a,info,&
       goto 9999
     end if
     
-    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
         
     nrst = nrst + 1 
     
@@ -307,7 +309,7 @@ subroutine psb_dgcr_vect(a,prec,b,x,eps,desc_a,info,&
       call psb_geaxpby(done, r, dzero, r, desc_a, info)   
       call psb_geaxpby(-alpha(j), c_scale(j), done, r, desc_a, info)   
       
-      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
       
       if (j >= irst) exit iteration
       
diff --git a/linsolve/impl/psb_dkrylov.f90 b/linsolve/impl/psb_dkrylov.f90
index 37daca4f9..e9aa21563 100644
--- a/linsolve/impl/psb_dkrylov.f90
+++ b/linsolve/impl/psb_dkrylov.f90
@@ -80,7 +80,7 @@
 !                                           estimate of) residual 
 ! 
 Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,irst,istop,cond)
+     & itmax,iter,err,itrace,irst,istop,cond,s1,s2)
 
   use psb_base_mod
   use psb_prec_mod,only : psb_dprec_type
@@ -97,11 +97,12 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, irst,istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   Real(psb_dpk_), Optional, Intent(out) :: err,cond
+  type(psb_d_vect_type), intent(inout), optional   :: s1, s2
 
 
   abstract interface
     subroutine psb_dkryl_vect(a,prec,b,x,eps,&
-         & desc_a,info,itmax,iter,err,itrace,istop)
+         & desc_a,info,itmax,iter,err,itrace,istop,s1,s2)
       import :: psb_ipk_, psb_dpk_, psb_desc_type, &
            & psb_dspmat_type, psb_dprec_type, psb_d_vect_type
       type(psb_dspmat_type), intent(in)    :: a
@@ -114,9 +115,10 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
       integer(psb_ipk_), optional, intent(in)        :: itmax, itrace,istop
       integer(psb_ipk_), optional, intent(out)       :: iter
       real(psb_dpk_), optional, intent(out) :: err
+      type(psb_d_vect_type), intent(inout), optional   :: s1, s2
     end subroutine psb_dkryl_vect
     Subroutine psb_dkryl_rest_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err, itrace,irst,istop)
+         &itmax,iter,err, itrace,irst,istop,s1,s2)
       import :: psb_ipk_, psb_dpk_, psb_desc_type, &
            & psb_dspmat_type, psb_dprec_type, psb_d_vect_type
       Type(psb_dspmat_type), Intent(in)    :: a
@@ -129,9 +131,10 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
       integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, irst,istop
       integer(psb_ipk_), Optional, Intent(out)       :: iter
       Real(psb_dpk_), Optional, Intent(out) :: err
+      type(psb_d_vect_type), intent(inout), optional   :: s1, s2
     end subroutine psb_dkryl_rest_vect
     Subroutine psb_dkryl_cond_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err, itrace,istop,cond)
+         &itmax,iter,err, itrace,istop,cond,s1,s2)
       import :: psb_ipk_, psb_dpk_, psb_desc_type, &
            & psb_dspmat_type, psb_dprec_type, psb_d_vect_type
       Type(psb_dspmat_type), Intent(in)    :: a
@@ -144,6 +147,7 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
       integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace,istop
       integer(psb_ipk_), Optional, Intent(out)       :: iter
       Real(psb_dpk_), Optional, Intent(out) :: err, cond
+      type(psb_d_vect_type), intent(inout), optional   :: s1, s2
     end subroutine psb_dkryl_cond_vect
   end interface
 
@@ -180,37 +184,37 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
   select case(psb_toupper(method))
   case('CG') 
     call  psb_dcg_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond)
+         &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond,s1=s1,s2=s2)
   case('FCG') 
     call  psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond)
+         &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond,s1=s1,s2=s2)
   case('GCR') 
     call  psb_dgcr_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop)
+         &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
   case('CGS') 
     call  psb_dcgs_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop)
+         &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
   case('BICG') 
     call  psb_dbicg_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop)
+         &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
   case('BICGSTAB') 
     call  psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop)
+         &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
   case('RGMRES','GMRES')
     call  psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop)
+         &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop,s1=s1,s2=s2)
   case('MINRES','PMINRES')
     call  psb_dminres_vect(a,prec,b,x,eps,desc_a,info,&
          &itmax,iter,err,itrace=itrace_,istop=istop)
   case('BICGSTABL')
     call  psb_dcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop)
+         &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop,s1=s1,s2=s2)
   case default
     if (me == 0) write(psb_err_unit,*) trim(name),&
          & ': Warning: Unknown method  ',method,&
          & ', defaulting to BiCGSTAB'
     call  psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop)
+         &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
   end select
 
   if ((info==psb_success_).and.do_alloc_wrk) call prec%free_wrk(info)
diff --git a/linsolve/impl/psb_drgmres.f90 b/linsolve/impl/psb_drgmres.f90
index 8da587173..32ab3b290 100644
--- a/linsolve/impl/psb_drgmres.f90
+++ b/linsolve/impl/psb_drgmres.f90
@@ -102,13 +102,15 @@
 !                                            stopped when  |r| <= eps * (|a||x|+|b|)
 !                                         2: err =  |r|/|b|; here the iteration is
 !                                            stopped when  |r| <= eps * |b|
+!                                         3: Same as 2 but with X and B scaled
+!                                            by s1 and s2 
 !                                         where r is the (preconditioned, recursive
 !                                         estimate of) residual. 
 !    irst   -  integer(optional)          Input: restart parameter 
 !
 
 subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,irst,istop)
+     & itmax,iter,err,itrace,irst,istop,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_d_linsolve_conv_mod
@@ -124,6 +126,7 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, irst,istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   Real(psb_dpk_), Optional, Intent(out) :: err
+  type(psb_d_vect_type), intent(inout), optional   :: s1, s2
 ! =   local data
   real(psb_dpk_), allocatable   :: aux(:)
   real(psb_dpk_), allocatable   :: c(:), s(:), h(:,:), rs(:), rst(:)
@@ -268,9 +271,20 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
   select case(istop_)
   case(psb_istop_ani_)
     ani = psb_spnrmi(a,desc_a,info)
-    bni = psb_geamax(b,desc_a,info)
+    if (present(s1)) then
+      call psb_gemlt(done,s1,b,dzero,v(1),desc_a,info)
+      bni = psb_geamax(v(1),desc_a,info)
+    else
+      bni = psb_geamax(b,desc_a,info)
+    end if
   case(psb_istop_bn2_)
-    bn2 = psb_genrm2(b,desc_a,info)    
+    if (present(s1)) then
+      call psb_gemlt(done,s1,b,dzero,v(1),desc_a,info)
+      bn2 = psb_genrm2(v(1),desc_a,info)
+    else
+      bn2 = psb_genrm2(b,desc_a,info)
+    end if
+      
   case(psb_istop_rn2_abs_)
     ! do nothing
   case(psb_istop_rrn2_)
@@ -282,6 +296,7 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
     end if
     
     call psb_spmm(-done,a,x,done,v(1),desc_a,info,work=aux)
+    if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info)
     if (info /= psb_success_) then 
       info=psb_err_from_subroutine_non_ 
       call psb_errpush(info,name)
@@ -323,7 +338,8 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
       call psb_errpush(info,name)
       goto 9999
     end if
-
+    if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info)
+    
     rs(1) = psb_genrm2(v(1),desc_a,info)
     rs(2:) = dzero
     if (info /= psb_success_) then 
@@ -378,8 +394,14 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
     inner:  Do i=1,nl
       itx  = itx + 1
 
-      call prec%apply(v(i),w1,desc_a,info)
+      if (present(s2)) then
+        call psb_gediv(v(i),s2,w,desc_a,info)
+        call prec%apply(w,w1,desc_a,info)
+      else
+        call prec%apply(v(i),w1,desc_a,info)
+      end if
       call psb_spmm(done,a,w1,dzero,w,desc_a,info,work=aux)
+      if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
       !
       call mgs(i,h,v,w,rs,c,s,desc_a,info)
 
@@ -391,10 +413,11 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
         !
         rst = rs
         call psb_geaxpby(done,x,dzero,xt,desc_a,info)
-        call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info)
+        call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info,s2=s2)
 
         call psb_geaxpby(done,b,dzero,w1,desc_a,info)
         call psb_spmm(-done,a,xt,done,w1,desc_a,info,work=aux)
+        if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
         rni = psb_geamax(w1,desc_a,info)
         xni = psb_geamax(xt,desc_a,info)
         errnum = rni
@@ -432,7 +455,8 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
           call psb_geaxpby(done,xt,dzero,x,desc_a,info)
 ! =          x = xt 
         case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
-          call  rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info)          !
+          call  rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2)
+          !
 
         end select
 
@@ -452,7 +476,7 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
       call psb_geaxpby(done,xt,dzero,x,desc_a,info)!      x = xt 
 
     case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
-      call  rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info)          !
+      call  rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2)          !
     end select
     
     if (itx >= itmax_) then 
@@ -523,11 +547,12 @@ contains
   ! Rebuild solution X from the space V using the factor
   ! stored in R
   !
-  subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info)
+  subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2)
     real(psb_dpk_)   :: c(:), s(:), rs(:), h(:,:)
     type(psb_d_vect_type) :: v(:), w, w1, x
     type(psb_desc_type) :: desc_a
     class(psb_dprec_type) :: prec
+    type(psb_d_vect_type), intent(inout), optional   :: s2
     integer(psb_ipk_) :: info
     integer(psb_ipk_) :: k,n
 
@@ -539,12 +564,13 @@ contains
     if (debug_level >= psb_debug_ext_) &
          & write(debug_unit,*) me,' ',trim(name),&
          & ' Rebuild x-> RS:',rs(1:n)
-    call w1%zero()
+    call w%zero()
     do k=1, n
-      call psb_geaxpby(rs(k),v(k),done,w1,desc_a,info)
+      call psb_geaxpby(rs(k),v(k),done,w,desc_a,info)
     end do
-    call prec%apply(w1,w,desc_a,info)
-    call psb_geaxpby(done,w,done,x,desc_a,info)
+    if (present(s2)) call psb_gediv(s2,w,desc_a,info)
+    call prec%apply(w,w1,desc_a,info)
+    call psb_geaxpby(done,w1,done,x,desc_a,info)
   end subroutine rebuildx
 
 end subroutine psb_drgmres_vect
diff --git a/linsolve/impl/psb_sbicg.f90 b/linsolve/impl/psb_sbicg.f90
index 168ee891d..88d81d977 100644
--- a/linsolve/impl/psb_sbicg.f90
+++ b/linsolve/impl/psb_sbicg.f90
@@ -95,7 +95,7 @@
 !
 
 subroutine psb_sbicg_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,istop)
+     & itmax,iter,err,itrace,istop,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_s_linsolve_conv_mod
@@ -111,6 +111,7 @@ subroutine psb_sbicg_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), optional, intent(in)        :: itmax, itrace, istop
   integer(psb_ipk_), optional, intent(out)       :: iter
   real(psb_spk_), optional, intent(out) :: err
+  type(psb_s_vect_type), intent(inout), optional   :: s1, s2
 ! !$   local data
   real(psb_spk_), allocatable, target   :: aux(:)
   type(psb_s_vect_type), allocatable, target :: wwrk(:)
@@ -236,7 +237,8 @@ subroutine psb_sbicg_vect(a,prec,b,x,eps,desc_a,info,&
   itx   = 0
 
 
-  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
+  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
+       & desc_a,stopdat,info,s1=s1,s2=s2)
   if (info /= psb_success_) Then 
      call psb_errpush(psb_err_from_subroutine_non_,name)
      goto 9999
@@ -262,7 +264,7 @@ subroutine psb_sbicg_vect(a,prec,b,x,eps,desc_a,info,&
     rho = szero
     
     ! Perhaps we already satisfy the convergence criterion...
-    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
     if (info /= psb_success_) Then 
       call psb_errpush(psb_err_from_subroutine_non_,name)
       goto 9999
@@ -316,7 +318,7 @@ subroutine psb_sbicg_vect(a,prec,b,x,eps,desc_a,info,&
       call psb_geaxpby(-alpha,q,sone,r,desc_a,info)
       call psb_geaxpby(-alpha,qt,sone,rt,desc_a,info)
 
-      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
       if (info /= psb_success_) Then 
         call psb_errpush(psb_err_from_subroutine_non_,name)
         goto 9999
diff --git a/linsolve/impl/psb_scg.F90 b/linsolve/impl/psb_scg.F90
index 6a1c750b4..b93603660 100644
--- a/linsolve/impl/psb_scg.F90
+++ b/linsolve/impl/psb_scg.F90
@@ -96,7 +96,7 @@
 ! 
 !
 subroutine psb_scg_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,istop,cond)
+     & itmax,iter,err,itrace,istop,cond,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_s_linsolve_conv_mod
@@ -112,6 +112,8 @@ subroutine psb_scg_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   Real(psb_spk_), Optional, Intent(out) :: err,cond
+  type(psb_s_vect_type), intent(inout), optional   :: s1, s2
+      
 ! =   Local data
   real(psb_spk_), allocatable, target   :: aux(:),td(:),tu(:),eig(:),ewrk(:)
   integer(psb_mpk_), allocatable :: ibl(:), ispl(:), iwrk(:)
@@ -253,7 +255,8 @@ subroutine psb_scg_vect(a,prec,b,x,eps,desc_a,info,&
 
     rho = szero
     
-    call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
+    call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
+         &desc_a,stopdat,info,s1=s1,s2=s2)
     if (info /= psb_success_) Then 
       call psb_errpush(psb_err_from_subroutine_non_,name)
       goto 9999
@@ -306,7 +309,7 @@ subroutine psb_scg_vect(a,prec,b,x,eps,desc_a,info,&
       call psb_geaxpby(alpha,p,sone,x,desc_a,info)
       call psb_geaxpby(-alpha,q,sone,r,desc_a,info)
 
-      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
       if (info /= psb_success_) Then 
         call psb_errpush(psb_err_from_subroutine_non_,name)
         goto 9999
diff --git a/linsolve/impl/psb_scgs.f90 b/linsolve/impl/psb_scgs.f90
index 8b2b74832..1e13d39a4 100644
--- a/linsolve/impl/psb_scgs.f90
+++ b/linsolve/impl/psb_scgs.f90
@@ -93,7 +93,7 @@
 !                                         estimate of) residual. 
 !
 Subroutine psb_scgs_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,istop)
+     & itmax,iter,err,itrace,istop,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_s_linsolve_conv_mod
@@ -109,6 +109,7 @@ Subroutine psb_scgs_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace,istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   Real(psb_spk_), Optional, Intent(out) :: err
+  type(psb_s_vect_type), intent(inout), optional   :: s1, s2
 ! =   local data
   real(psb_spk_), allocatable, target   :: aux(:)
   type(psb_s_vect_type), allocatable, target :: wwrk(:)
@@ -223,7 +224,8 @@ Subroutine psb_scgs_vect(a,prec,b,x,eps,desc_a,info,&
 
   itx   = 0
 
-  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
+  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
+       & desc_a,stopdat,info,s1=s1,s2=s2)
   if (info /= psb_success_) Then 
      call psb_errpush(psb_err_from_subroutine_non_,name)
      goto 9999
@@ -246,7 +248,7 @@ Subroutine psb_scgs_vect(a,prec,b,x,eps,desc_a,info,&
     
 
     ! Perhaps we already satisfy the convergence criterion...
-    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
     if (info /= psb_success_) Then 
       call psb_errpush(psb_err_from_subroutine_non_,name)
       goto 9999
@@ -320,7 +322,7 @@ Subroutine psb_scgs_vect(a,prec,b,x,eps,desc_a,info,&
          goto 9999
       end if
 
-      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
       if (info /= psb_success_) Then 
         call psb_errpush(psb_err_from_subroutine_non_,name)
         goto 9999
diff --git a/linsolve/impl/psb_scgstab.f90 b/linsolve/impl/psb_scgstab.f90
index fe6e0e2cf..f1e72b73b 100644
--- a/linsolve/impl/psb_scgstab.f90
+++ b/linsolve/impl/psb_scgstab.f90
@@ -93,7 +93,7 @@
 !                                         where r is the (preconditioned, recursive
 !                                         estimate of) residual. 
 !
-Subroutine psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
+Subroutine psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_s_linsolve_conv_mod
@@ -109,6 +109,7 @@ Subroutine psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
   integer(psb_ipk_), Optional, Intent(in)      :: itmax, itrace, istop
   integer(psb_ipk_), Optional, Intent(out)     :: iter
   Real(psb_spk_), Optional, Intent(out) :: err
+  type(psb_s_vect_type), intent(inout), optional   :: s1, s2
 ! =   Local data
   real(psb_spk_), allocatable, target   :: aux(:),wwrk(:,:)
   type(psb_s_vect_type) :: q, r, p, v, s, t, z, f
@@ -235,7 +236,8 @@ Subroutine psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
   End If
 
   itx   = 0
-  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
+  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
+       & desc_a,stopdat,info,s1=s1,s2=s2)
   if (psb_errstatus_fatal()) Then 
     call psb_errpush(psb_err_from_subroutine_non_,name)
     goto 9999
@@ -252,7 +254,7 @@ Subroutine psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
     call psb_geaxpby(sone,r,szero,q,desc_a,info)
 
     ! Perhaps we already satisfy the convergence criterion...
-    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
 
     if (psb_errstatus_fatal()) then 
       info=psb_err_from_subroutine_
@@ -372,7 +374,7 @@ Subroutine psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
       call psb_geaxpby(omega,z,sone,x,desc_a,info)
       call psb_geaxpby(sone,s,szero,r,desc_a,info)
       call psb_geaxpby(-omega,t,sone,r,desc_a,info)
-      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
 
       if (psb_errstatus_fatal()) Then 
         call psb_errpush(psb_err_from_subroutine_,name,a_err='X/R update ')
diff --git a/linsolve/impl/psb_scgstabl.f90 b/linsolve/impl/psb_scgstabl.f90
index ab504cb83..a9d5e783a 100644
--- a/linsolve/impl/psb_scgstabl.f90
+++ b/linsolve/impl/psb_scgstabl.f90
@@ -104,7 +104,7 @@
 !
 !
 Subroutine psb_scgstabl_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,irst,istop)
+     & itmax,iter,err,itrace,irst,istop,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_s_linsolve_conv_mod
@@ -120,6 +120,7 @@ Subroutine psb_scgstabl_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, irst,istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   Real(psb_spk_), Optional, Intent(out) :: err
+  type(psb_s_vect_type), intent(inout), optional   :: s1, s2
 ! =   local data
   real(psb_spk_), allocatable, target   :: aux(:), gamma(:),&
        & gamma1(:), gamma2(:), taum(:,:), sigma(:)
@@ -267,7 +268,8 @@ Subroutine psb_scgstabl_vect(a,prec,b,x,eps,desc_a,info,&
   rt0 => wwrk(10)
   
 
-  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
+  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
+       & desc_a,stopdat,info,s1=s1,s2=s2)
   if (info /= psb_success_) Then 
      call psb_errpush(psb_err_from_subroutine_non_,name)
      goto 9999
@@ -305,7 +307,7 @@ Subroutine psb_scgstabl_vect(a,prec,b,x,eps,desc_a,info,&
          & write(debug_unit,*) me,' ',trim(name),&
          & ' on entry to amax: b: ',b%get_nrows()
 
-    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
     if (info /= psb_success_) Then 
       call psb_errpush(psb_err_from_subroutine_non_,name)
       goto 9999
@@ -409,7 +411,7 @@ Subroutine psb_scgstabl_vect(a,prec,b,x,eps,desc_a,info,&
         call psb_geaxpby(-gamma1(j),rh(j),sone,rh(0),desc_a,info)        
       enddo
       
-      if (psb_check_conv(methdname,itx,x,rh(0),desc_a,stopdat,info)) exit restart
+      if (psb_check_conv(methdname,itx,x,rh(0),desc_a,stopdat,info,s1=s1)) exit restart
       if (info /= psb_success_) Then 
         call psb_errpush(psb_err_from_subroutine_non_,name)
         goto 9999
diff --git a/linsolve/impl/psb_sfcg.F90 b/linsolve/impl/psb_sfcg.F90
index a18470e02..e55cd6a8b 100644
--- a/linsolve/impl/psb_sfcg.F90
+++ b/linsolve/impl/psb_sfcg.F90
@@ -104,7 +104,7 @@
 ! 
 !
 subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,istop,cond)
+     & itmax,iter,err,itrace,istop,cond,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_s_linsolve_conv_mod
@@ -120,6 +120,7 @@ subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   real(psb_spk_), Optional, Intent(out) :: err,cond
+  type(psb_s_vect_type), intent(inout), optional   :: s1, s2
 ! =   Local data
   type(psb_s_vect_type)  :: v, w, d , q, r
   real(psb_spk_) :: alpha, beta, delta, gamma, theta
@@ -227,7 +228,7 @@ subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,&
        & scratch=.true.,mold=x%v)
 
   call psb_init_conv(methdname,istop_,itrace_,itmax_,&
-       & a,x,b,eps,desc_a,stopdat,info)
+       & a,x,b,eps,desc_a,stopdat,info,s1=s1,s2=s2)
   itx = 0 
 
   restart: do 
@@ -246,7 +247,7 @@ subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,&
     end if
     
 
-    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) then
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) then
       if (debug.and.(me==0)) write(0,*) name,' Exit on  convergence from restart'
       exit restart
     end if
@@ -302,7 +303,7 @@ subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,&
 
       itx = itx + 1
 
-      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) then
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) then
         if (debug.and.(me==0)) write(0,*) name,' Exit on  convergence from iteration'
         exit restart
       end if
diff --git a/linsolve/impl/psb_sgcr.f90 b/linsolve/impl/psb_sgcr.f90
index aeccfcff9..8966c280d 100644
--- a/linsolve/impl/psb_sgcr.f90
+++ b/linsolve/impl/psb_sgcr.f90
@@ -106,7 +106,7 @@
 !
 
 subroutine psb_sgcr_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace, irst, istop)
+     & itmax,iter,err,itrace, irst, istop,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_s_linsolve_conv_mod
@@ -124,6 +124,7 @@ subroutine psb_sgcr_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, irst, istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   real(psb_spk_), Optional, Intent(out) :: err
+  type(psb_s_vect_type), intent(inout), optional   :: s1, s2
   ! =   local data
   real(psb_spk_), allocatable   :: alpha(:), h(:,:)
   type(psb_s_vect_type), allocatable :: z(:), c(:), c_scale(:)
@@ -253,7 +254,8 @@ subroutine psb_sgcr_vect(a,prec,b,x,eps,desc_a,info,&
   itx = 0
   
   nrst = -1
-  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
+  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
+       & desc_a,stopdat,info,s1=s1,s2=s2)
   restart: do 
     if (itx>= itmax_) exit restart 
     h = szero
@@ -276,7 +278,7 @@ subroutine psb_sgcr_vect(a,prec,b,x,eps,desc_a,info,&
       goto 9999
     end if
     
-    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
         
     nrst = nrst + 1 
     
@@ -307,7 +309,7 @@ subroutine psb_sgcr_vect(a,prec,b,x,eps,desc_a,info,&
       call psb_geaxpby(sone, r, szero, r, desc_a, info)   
       call psb_geaxpby(-alpha(j), c_scale(j), sone, r, desc_a, info)   
       
-      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
       
       if (j >= irst) exit iteration
       
diff --git a/linsolve/impl/psb_skrylov.f90 b/linsolve/impl/psb_skrylov.f90
index cce35cf9a..c7bf798f9 100644
--- a/linsolve/impl/psb_skrylov.f90
+++ b/linsolve/impl/psb_skrylov.f90
@@ -80,7 +80,7 @@
 !                                           estimate of) residual 
 ! 
 Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,irst,istop,cond)
+     & itmax,iter,err,itrace,irst,istop,cond,s1,s2)
 
   use psb_base_mod
   use psb_prec_mod,only : psb_sprec_type
@@ -97,11 +97,12 @@ Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, irst,istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   Real(psb_spk_), Optional, Intent(out) :: err,cond
+  type(psb_s_vect_type), intent(inout), optional   :: s1, s2
 
 
   abstract interface
     subroutine psb_skryl_vect(a,prec,b,x,eps,&
-         & desc_a,info,itmax,iter,err,itrace,istop)
+         & desc_a,info,itmax,iter,err,itrace,istop,s1,s2)
       import :: psb_ipk_, psb_spk_, psb_desc_type, &
            & psb_sspmat_type, psb_sprec_type, psb_s_vect_type
       type(psb_sspmat_type), intent(in)    :: a
@@ -114,9 +115,10 @@ Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
       integer(psb_ipk_), optional, intent(in)        :: itmax, itrace,istop
       integer(psb_ipk_), optional, intent(out)       :: iter
       real(psb_spk_), optional, intent(out) :: err
+      type(psb_s_vect_type), intent(inout), optional   :: s1, s2
     end subroutine psb_skryl_vect
     Subroutine psb_skryl_rest_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err, itrace,irst,istop)
+         &itmax,iter,err, itrace,irst,istop,s1,s2)
       import :: psb_ipk_, psb_spk_, psb_desc_type, &
            & psb_sspmat_type, psb_sprec_type, psb_s_vect_type
       Type(psb_sspmat_type), Intent(in)    :: a
@@ -129,9 +131,10 @@ Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
       integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, irst,istop
       integer(psb_ipk_), Optional, Intent(out)       :: iter
       Real(psb_spk_), Optional, Intent(out) :: err
+      type(psb_s_vect_type), intent(inout), optional   :: s1, s2
     end subroutine psb_skryl_rest_vect
     Subroutine psb_skryl_cond_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err, itrace,istop,cond)
+         &itmax,iter,err, itrace,istop,cond,s1,s2)
       import :: psb_ipk_, psb_spk_, psb_desc_type, &
            & psb_sspmat_type, psb_sprec_type, psb_s_vect_type
       Type(psb_sspmat_type), Intent(in)    :: a
@@ -144,6 +147,7 @@ Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
       integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace,istop
       integer(psb_ipk_), Optional, Intent(out)       :: iter
       Real(psb_spk_), Optional, Intent(out) :: err, cond
+      type(psb_s_vect_type), intent(inout), optional   :: s1, s2
     end subroutine psb_skryl_cond_vect
   end interface
 
@@ -180,37 +184,37 @@ Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
   select case(psb_toupper(method))
   case('CG') 
     call  psb_scg_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond)
+         &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond,s1=s1,s2=s2)
   case('FCG') 
     call  psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond)
+         &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond,s1=s1,s2=s2)
   case('GCR') 
     call  psb_sgcr_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop)
+         &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
   case('CGS') 
     call  psb_scgs_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop)
+         &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
   case('BICG') 
     call  psb_sbicg_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop)
+         &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
   case('BICGSTAB') 
     call  psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop)
+         &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
   case('RGMRES','GMRES')
     call  psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop)
+         &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop,s1=s1,s2=s2)
   case('MINRES','PMINRES')
     call  psb_sminres_vect(a,prec,b,x,eps,desc_a,info,&
          &itmax,iter,err,itrace=itrace_,istop=istop)
   case('BICGSTABL')
     call  psb_scgstabl_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop)
+         &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop,s1=s1,s2=s2)
   case default
     if (me == 0) write(psb_err_unit,*) trim(name),&
          & ': Warning: Unknown method  ',method,&
          & ', defaulting to BiCGSTAB'
     call  psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop)
+         &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
   end select
 
   if ((info==psb_success_).and.do_alloc_wrk) call prec%free_wrk(info)
diff --git a/linsolve/impl/psb_srgmres.f90 b/linsolve/impl/psb_srgmres.f90
index 208665dbe..b986c6cc5 100644
--- a/linsolve/impl/psb_srgmres.f90
+++ b/linsolve/impl/psb_srgmres.f90
@@ -102,13 +102,15 @@
 !                                            stopped when  |r| <= eps * (|a||x|+|b|)
 !                                         2: err =  |r|/|b|; here the iteration is
 !                                            stopped when  |r| <= eps * |b|
+!                                         3: Same as 2 but with X and B scaled
+!                                            by s1 and s2 
 !                                         where r is the (preconditioned, recursive
 !                                         estimate of) residual. 
 !    irst   -  integer(optional)          Input: restart parameter 
 !
 
 subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,irst,istop)
+     & itmax,iter,err,itrace,irst,istop,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_s_linsolve_conv_mod
@@ -124,6 +126,7 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, irst,istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   Real(psb_spk_), Optional, Intent(out) :: err
+  type(psb_s_vect_type), intent(inout), optional   :: s1, s2
 ! =   local data
   real(psb_spk_), allocatable   :: aux(:)
   real(psb_spk_), allocatable   :: c(:), s(:), h(:,:), rs(:), rst(:)
@@ -268,9 +271,20 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
   select case(istop_)
   case(psb_istop_ani_)
     ani = psb_spnrmi(a,desc_a,info)
-    bni = psb_geamax(b,desc_a,info)
+    if (present(s1)) then
+      call psb_gemlt(sone,s1,b,szero,v(1),desc_a,info)
+      bni = psb_geamax(v(1),desc_a,info)
+    else
+      bni = psb_geamax(b,desc_a,info)
+    end if
   case(psb_istop_bn2_)
-    bn2 = psb_genrm2(b,desc_a,info)    
+    if (present(s1)) then
+      call psb_gemlt(sone,s1,b,szero,v(1),desc_a,info)
+      bn2 = psb_genrm2(v(1),desc_a,info)
+    else
+      bn2 = psb_genrm2(b,desc_a,info)
+    end if
+      
   case(psb_istop_rn2_abs_)
     ! do nothing
   case(psb_istop_rrn2_)
@@ -282,6 +296,7 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
     end if
     
     call psb_spmm(-sone,a,x,sone,v(1),desc_a,info,work=aux)
+    if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info)
     if (info /= psb_success_) then 
       info=psb_err_from_subroutine_non_ 
       call psb_errpush(info,name)
@@ -323,7 +338,8 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
       call psb_errpush(info,name)
       goto 9999
     end if
-
+    if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info)
+    
     rs(1) = psb_genrm2(v(1),desc_a,info)
     rs(2:) = szero
     if (info /= psb_success_) then 
@@ -378,8 +394,14 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
     inner:  Do i=1,nl
       itx  = itx + 1
 
-      call prec%apply(v(i),w1,desc_a,info)
+      if (present(s2)) then
+        call psb_gediv(v(i),s2,w,desc_a,info)
+        call prec%apply(w,w1,desc_a,info)
+      else
+        call prec%apply(v(i),w1,desc_a,info)
+      end if
       call psb_spmm(sone,a,w1,szero,w,desc_a,info,work=aux)
+      if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
       !
       call mgs(i,h,v,w,rs,c,s,desc_a,info)
 
@@ -391,10 +413,11 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
         !
         rst = rs
         call psb_geaxpby(sone,x,szero,xt,desc_a,info)
-        call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info)
+        call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info,s2=s2)
 
         call psb_geaxpby(sone,b,szero,w1,desc_a,info)
         call psb_spmm(-sone,a,xt,sone,w1,desc_a,info,work=aux)
+        if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
         rni = psb_geamax(w1,desc_a,info)
         xni = psb_geamax(xt,desc_a,info)
         errnum = rni
@@ -432,7 +455,8 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
           call psb_geaxpby(sone,xt,szero,x,desc_a,info)
 ! =          x = xt 
         case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
-          call  rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info)          !
+          call  rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2)
+          !
 
         end select
 
@@ -452,7 +476,7 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
       call psb_geaxpby(sone,xt,szero,x,desc_a,info)!      x = xt 
 
     case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
-      call  rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info)          !
+      call  rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2)          !
     end select
     
     if (itx >= itmax_) then 
@@ -523,11 +547,12 @@ contains
   ! Rebuild solution X from the space V using the factor
   ! stored in R
   !
-  subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info)
+  subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2)
     real(psb_spk_)   :: c(:), s(:), rs(:), h(:,:)
     type(psb_s_vect_type) :: v(:), w, w1, x
     type(psb_desc_type) :: desc_a
     class(psb_sprec_type) :: prec
+    type(psb_s_vect_type), intent(inout), optional   :: s2
     integer(psb_ipk_) :: info
     integer(psb_ipk_) :: k,n
 
@@ -539,12 +564,13 @@ contains
     if (debug_level >= psb_debug_ext_) &
          & write(debug_unit,*) me,' ',trim(name),&
          & ' Rebuild x-> RS:',rs(1:n)
-    call w1%zero()
+    call w%zero()
     do k=1, n
-      call psb_geaxpby(rs(k),v(k),sone,w1,desc_a,info)
+      call psb_geaxpby(rs(k),v(k),sone,w,desc_a,info)
     end do
-    call prec%apply(w1,w,desc_a,info)
-    call psb_geaxpby(sone,w,sone,x,desc_a,info)
+    if (present(s2)) call psb_gediv(s2,w,desc_a,info)
+    call prec%apply(w,w1,desc_a,info)
+    call psb_geaxpby(sone,w1,sone,x,desc_a,info)
   end subroutine rebuildx
 
 end subroutine psb_srgmres_vect
diff --git a/linsolve/impl/psb_zbicg.f90 b/linsolve/impl/psb_zbicg.f90
index 0396fe06c..55e482423 100644
--- a/linsolve/impl/psb_zbicg.f90
+++ b/linsolve/impl/psb_zbicg.f90
@@ -95,7 +95,7 @@
 !
 
 subroutine psb_zbicg_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,istop)
+     & itmax,iter,err,itrace,istop,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_z_linsolve_conv_mod
@@ -111,6 +111,7 @@ subroutine psb_zbicg_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), optional, intent(in)        :: itmax, itrace, istop
   integer(psb_ipk_), optional, intent(out)       :: iter
   real(psb_dpk_), optional, intent(out) :: err
+  type(psb_z_vect_type), intent(inout), optional   :: s1, s2
 ! !$   local data
   complex(psb_dpk_), allocatable, target   :: aux(:)
   type(psb_z_vect_type), allocatable, target :: wwrk(:)
@@ -236,7 +237,8 @@ subroutine psb_zbicg_vect(a,prec,b,x,eps,desc_a,info,&
   itx   = 0
 
 
-  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
+  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
+       & desc_a,stopdat,info,s1=s1,s2=s2)
   if (info /= psb_success_) Then 
      call psb_errpush(psb_err_from_subroutine_non_,name)
      goto 9999
@@ -262,7 +264,7 @@ subroutine psb_zbicg_vect(a,prec,b,x,eps,desc_a,info,&
     rho = zzero
     
     ! Perhaps we already satisfy the convergence criterion...
-    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
     if (info /= psb_success_) Then 
       call psb_errpush(psb_err_from_subroutine_non_,name)
       goto 9999
@@ -316,7 +318,7 @@ subroutine psb_zbicg_vect(a,prec,b,x,eps,desc_a,info,&
       call psb_geaxpby(-alpha,q,zone,r,desc_a,info)
       call psb_geaxpby(-alpha,qt,zone,rt,desc_a,info)
 
-      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
       if (info /= psb_success_) Then 
         call psb_errpush(psb_err_from_subroutine_non_,name)
         goto 9999
diff --git a/linsolve/impl/psb_zcg.F90 b/linsolve/impl/psb_zcg.F90
index aaf654ea3..cb34f85d1 100644
--- a/linsolve/impl/psb_zcg.F90
+++ b/linsolve/impl/psb_zcg.F90
@@ -96,7 +96,7 @@
 ! 
 !
 subroutine psb_zcg_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,istop,cond)
+     & itmax,iter,err,itrace,istop,cond,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_z_linsolve_conv_mod
@@ -112,6 +112,8 @@ subroutine psb_zcg_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   Real(psb_dpk_), Optional, Intent(out) :: err,cond
+  type(psb_z_vect_type), intent(inout), optional   :: s1, s2
+      
 ! =   Local data
   complex(psb_dpk_), allocatable, target   :: aux(:),td(:),tu(:),eig(:),ewrk(:)
   integer(psb_mpk_), allocatable :: ibl(:), ispl(:), iwrk(:)
@@ -245,7 +247,8 @@ subroutine psb_zcg_vect(a,prec,b,x,eps,desc_a,info,&
 
     rho = zzero
     
-    call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
+    call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
+         &desc_a,stopdat,info,s1=s1,s2=s2)
     if (info /= psb_success_) Then 
       call psb_errpush(psb_err_from_subroutine_non_,name)
       goto 9999
@@ -289,7 +292,7 @@ subroutine psb_zcg_vect(a,prec,b,x,eps,desc_a,info,&
       call psb_geaxpby(alpha,p,zone,x,desc_a,info)
       call psb_geaxpby(-alpha,q,zone,r,desc_a,info)
 
-      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
       if (info /= psb_success_) Then 
         call psb_errpush(psb_err_from_subroutine_non_,name)
         goto 9999
diff --git a/linsolve/impl/psb_zcgs.f90 b/linsolve/impl/psb_zcgs.f90
index 969bf1bf4..208976112 100644
--- a/linsolve/impl/psb_zcgs.f90
+++ b/linsolve/impl/psb_zcgs.f90
@@ -93,7 +93,7 @@
 !                                         estimate of) residual. 
 !
 Subroutine psb_zcgs_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,istop)
+     & itmax,iter,err,itrace,istop,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_z_linsolve_conv_mod
@@ -109,6 +109,7 @@ Subroutine psb_zcgs_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace,istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   Real(psb_dpk_), Optional, Intent(out) :: err
+  type(psb_z_vect_type), intent(inout), optional   :: s1, s2
 ! =   local data
   complex(psb_dpk_), allocatable, target   :: aux(:)
   type(psb_z_vect_type), allocatable, target :: wwrk(:)
@@ -223,7 +224,8 @@ Subroutine psb_zcgs_vect(a,prec,b,x,eps,desc_a,info,&
 
   itx   = 0
 
-  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
+  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
+       & desc_a,stopdat,info,s1=s1,s2=s2)
   if (info /= psb_success_) Then 
      call psb_errpush(psb_err_from_subroutine_non_,name)
      goto 9999
@@ -246,7 +248,7 @@ Subroutine psb_zcgs_vect(a,prec,b,x,eps,desc_a,info,&
     
 
     ! Perhaps we already satisfy the convergence criterion...
-    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
     if (info /= psb_success_) Then 
       call psb_errpush(psb_err_from_subroutine_non_,name)
       goto 9999
@@ -320,7 +322,7 @@ Subroutine psb_zcgs_vect(a,prec,b,x,eps,desc_a,info,&
          goto 9999
       end if
 
-      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
       if (info /= psb_success_) Then 
         call psb_errpush(psb_err_from_subroutine_non_,name)
         goto 9999
diff --git a/linsolve/impl/psb_zcgstab.f90 b/linsolve/impl/psb_zcgstab.f90
index 09383bfdd..5777fa209 100644
--- a/linsolve/impl/psb_zcgstab.f90
+++ b/linsolve/impl/psb_zcgstab.f90
@@ -93,7 +93,7 @@
 !                                         where r is the (preconditioned, recursive
 !                                         estimate of) residual. 
 !
-Subroutine psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
+Subroutine psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_z_linsolve_conv_mod
@@ -109,6 +109,7 @@ Subroutine psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
   integer(psb_ipk_), Optional, Intent(in)      :: itmax, itrace, istop
   integer(psb_ipk_), Optional, Intent(out)     :: iter
   Real(psb_dpk_), Optional, Intent(out) :: err
+  type(psb_z_vect_type), intent(inout), optional   :: s1, s2
 ! =   Local data
   complex(psb_dpk_), allocatable, target   :: aux(:),wwrk(:,:)
   type(psb_z_vect_type) :: q, r, p, v, s, t, z, f
@@ -235,7 +236,8 @@ Subroutine psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
   End If
 
   itx   = 0
-  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
+  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
+       & desc_a,stopdat,info,s1=s1,s2=s2)
   if (psb_errstatus_fatal()) Then 
     call psb_errpush(psb_err_from_subroutine_non_,name)
     goto 9999
@@ -252,7 +254,7 @@ Subroutine psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
     call psb_geaxpby(zone,r,zzero,q,desc_a,info)
 
     ! Perhaps we already satisfy the convergence criterion...
-    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
 
     if (psb_errstatus_fatal()) then 
       info=psb_err_from_subroutine_
@@ -372,7 +374,7 @@ Subroutine psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
       call psb_geaxpby(omega,z,zone,x,desc_a,info)
       call psb_geaxpby(zone,s,zzero,r,desc_a,info)
       call psb_geaxpby(-omega,t,zone,r,desc_a,info)
-      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
 
       if (psb_errstatus_fatal()) Then 
         call psb_errpush(psb_err_from_subroutine_,name,a_err='X/R update ')
diff --git a/linsolve/impl/psb_zcgstabl.f90 b/linsolve/impl/psb_zcgstabl.f90
index 8bbf4ca6b..c99cf7d20 100644
--- a/linsolve/impl/psb_zcgstabl.f90
+++ b/linsolve/impl/psb_zcgstabl.f90
@@ -104,7 +104,7 @@
 !
 !
 Subroutine psb_zcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,irst,istop)
+     & itmax,iter,err,itrace,irst,istop,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_z_linsolve_conv_mod
@@ -120,6 +120,7 @@ Subroutine psb_zcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, irst,istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   Real(psb_dpk_), Optional, Intent(out) :: err
+  type(psb_z_vect_type), intent(inout), optional   :: s1, s2
 ! =   local data
   complex(psb_dpk_), allocatable, target   :: aux(:), gamma(:),&
        & gamma1(:), gamma2(:), taum(:,:), sigma(:)
@@ -267,7 +268,8 @@ Subroutine psb_zcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
   rt0 => wwrk(10)
   
 
-  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
+  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
+       & desc_a,stopdat,info,s1=s1,s2=s2)
   if (info /= psb_success_) Then 
      call psb_errpush(psb_err_from_subroutine_non_,name)
      goto 9999
@@ -305,7 +307,7 @@ Subroutine psb_zcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
          & write(debug_unit,*) me,' ',trim(name),&
          & ' on entry to amax: b: ',b%get_nrows()
 
-    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
     if (info /= psb_success_) Then 
       call psb_errpush(psb_err_from_subroutine_non_,name)
       goto 9999
@@ -409,7 +411,7 @@ Subroutine psb_zcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
         call psb_geaxpby(-gamma1(j),rh(j),zone,rh(0),desc_a,info)        
       enddo
       
-      if (psb_check_conv(methdname,itx,x,rh(0),desc_a,stopdat,info)) exit restart
+      if (psb_check_conv(methdname,itx,x,rh(0),desc_a,stopdat,info,s1=s1)) exit restart
       if (info /= psb_success_) Then 
         call psb_errpush(psb_err_from_subroutine_non_,name)
         goto 9999
diff --git a/linsolve/impl/psb_zfcg.F90 b/linsolve/impl/psb_zfcg.F90
index 3dfda5500..7f2806129 100644
--- a/linsolve/impl/psb_zfcg.F90
+++ b/linsolve/impl/psb_zfcg.F90
@@ -104,7 +104,7 @@
 ! 
 !
 subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,istop,cond)
+     & itmax,iter,err,itrace,istop,cond,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_z_linsolve_conv_mod
@@ -120,6 +120,7 @@ subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   real(psb_dpk_), Optional, Intent(out) :: err,cond
+  type(psb_z_vect_type), intent(inout), optional   :: s1, s2
 ! =   Local data
   type(psb_z_vect_type)  :: v, w, d , q, r
   complex(psb_dpk_) :: alpha, beta, delta, gamma, theta
@@ -227,7 +228,7 @@ subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,&
        & scratch=.true.,mold=x%v)
 
   call psb_init_conv(methdname,istop_,itrace_,itmax_,&
-       & a,x,b,eps,desc_a,stopdat,info)
+       & a,x,b,eps,desc_a,stopdat,info,s1=s1,s2=s2)
   itx = 0 
 
   restart: do 
@@ -246,7 +247,7 @@ subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,&
     end if
     
 
-    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) then
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) then
       if (debug.and.(me==0)) write(0,*) name,' Exit on  convergence from restart'
       exit restart
     end if
@@ -302,7 +303,7 @@ subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,&
 
       itx = itx + 1
 
-      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) then
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) then
         if (debug.and.(me==0)) write(0,*) name,' Exit on  convergence from iteration'
         exit restart
       end if
diff --git a/linsolve/impl/psb_zgcr.f90 b/linsolve/impl/psb_zgcr.f90
index 812143be8..e2c64d91a 100644
--- a/linsolve/impl/psb_zgcr.f90
+++ b/linsolve/impl/psb_zgcr.f90
@@ -106,7 +106,7 @@
 !
 
 subroutine psb_zgcr_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace, irst, istop)
+     & itmax,iter,err,itrace, irst, istop,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_z_linsolve_conv_mod
@@ -124,6 +124,7 @@ subroutine psb_zgcr_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, irst, istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   real(psb_dpk_), Optional, Intent(out) :: err
+  type(psb_z_vect_type), intent(inout), optional   :: s1, s2
   ! =   local data
   complex(psb_dpk_), allocatable   :: alpha(:), h(:,:)
   type(psb_z_vect_type), allocatable :: z(:), c(:), c_scale(:)
@@ -253,7 +254,8 @@ subroutine psb_zgcr_vect(a,prec,b,x,eps,desc_a,info,&
   itx = 0
   
   nrst = -1
-  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
+  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,&
+       & desc_a,stopdat,info,s1=s1,s2=s2)
   restart: do 
     if (itx>= itmax_) exit restart 
     h = zzero
@@ -276,7 +278,7 @@ subroutine psb_zgcr_vect(a,prec,b,x,eps,desc_a,info,&
       goto 9999
     end if
     
-    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
         
     nrst = nrst + 1 
     
@@ -307,7 +309,7 @@ subroutine psb_zgcr_vect(a,prec,b,x,eps,desc_a,info,&
       call psb_geaxpby(zone, r, zzero, r, desc_a, info)   
       call psb_geaxpby(-alpha(j), c_scale(j), zone, r, desc_a, info)   
       
-      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info,s1=s1)) exit restart
       
       if (j >= irst) exit iteration
       
diff --git a/linsolve/impl/psb_zkrylov.f90 b/linsolve/impl/psb_zkrylov.f90
index 3baff8c34..04a5f839b 100644
--- a/linsolve/impl/psb_zkrylov.f90
+++ b/linsolve/impl/psb_zkrylov.f90
@@ -80,7 +80,7 @@
 !                                           estimate of) residual 
 ! 
 Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,irst,istop,cond)
+     & itmax,iter,err,itrace,irst,istop,cond,s1,s2)
 
   use psb_base_mod
   use psb_prec_mod,only : psb_zprec_type
@@ -97,11 +97,12 @@ Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, irst,istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   Real(psb_dpk_), Optional, Intent(out) :: err,cond
+  type(psb_z_vect_type), intent(inout), optional   :: s1, s2
 
 
   abstract interface
     subroutine psb_zkryl_vect(a,prec,b,x,eps,&
-         & desc_a,info,itmax,iter,err,itrace,istop)
+         & desc_a,info,itmax,iter,err,itrace,istop,s1,s2)
       import :: psb_ipk_, psb_dpk_, psb_desc_type, &
            & psb_zspmat_type, psb_zprec_type, psb_z_vect_type
       type(psb_zspmat_type), intent(in)    :: a
@@ -114,9 +115,10 @@ Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
       integer(psb_ipk_), optional, intent(in)        :: itmax, itrace,istop
       integer(psb_ipk_), optional, intent(out)       :: iter
       real(psb_dpk_), optional, intent(out) :: err
+      type(psb_z_vect_type), intent(inout), optional   :: s1, s2
     end subroutine psb_zkryl_vect
     Subroutine psb_zkryl_rest_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err, itrace,irst,istop)
+         &itmax,iter,err, itrace,irst,istop,s1,s2)
       import :: psb_ipk_, psb_dpk_, psb_desc_type, &
            & psb_zspmat_type, psb_zprec_type, psb_z_vect_type
       Type(psb_zspmat_type), Intent(in)    :: a
@@ -129,9 +131,10 @@ Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
       integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, irst,istop
       integer(psb_ipk_), Optional, Intent(out)       :: iter
       Real(psb_dpk_), Optional, Intent(out) :: err
+      type(psb_z_vect_type), intent(inout), optional   :: s1, s2
     end subroutine psb_zkryl_rest_vect
     Subroutine psb_zkryl_cond_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err, itrace,istop,cond)
+         &itmax,iter,err, itrace,istop,cond,s1,s2)
       import :: psb_ipk_, psb_dpk_, psb_desc_type, &
            & psb_zspmat_type, psb_zprec_type, psb_z_vect_type
       Type(psb_zspmat_type), Intent(in)    :: a
@@ -144,6 +147,7 @@ Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
       integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace,istop
       integer(psb_ipk_), Optional, Intent(out)       :: iter
       Real(psb_dpk_), Optional, Intent(out) :: err, cond
+      type(psb_z_vect_type), intent(inout), optional   :: s1, s2
     end subroutine psb_zkryl_cond_vect
   end interface
 
@@ -180,37 +184,37 @@ Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
   select case(psb_toupper(method))
   case('CG') 
     call  psb_zcg_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond)
+         &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond,s1=s1,s2=s2)
   case('FCG') 
     call  psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond)
+         &itmax,iter,err,itrace=itrace_,istop=istop,cond=cond,s1=s1,s2=s2)
   case('GCR') 
     call  psb_zgcr_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop)
+         &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
   case('CGS') 
     call  psb_zcgs_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop)
+         &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
   case('BICG') 
     call  psb_zbicg_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop)
+         &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
   case('BICGSTAB') 
     call  psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop)
+         &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
   case('RGMRES','GMRES')
     call  psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop)
+         &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop,s1=s1,s2=s2)
   case('MINRES','PMINRES')
     call  psb_zminres_vect(a,prec,b,x,eps,desc_a,info,&
          &itmax,iter,err,itrace=itrace_,istop=istop)
   case('BICGSTABL')
     call  psb_zcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop)
+         &itmax,iter,err,itrace=itrace_,irst=irst,istop=istop,s1=s1,s2=s2)
   case default
     if (me == 0) write(psb_err_unit,*) trim(name),&
          & ': Warning: Unknown method  ',method,&
          & ', defaulting to BiCGSTAB'
     call  psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,&
-         &itmax,iter,err,itrace=itrace_,istop=istop)
+         &itmax,iter,err,itrace=itrace_,istop=istop,s1=s1,s2=s2)
   end select
 
   if ((info==psb_success_).and.do_alloc_wrk) call prec%free_wrk(info)
diff --git a/linsolve/impl/psb_zrgmres.f90 b/linsolve/impl/psb_zrgmres.f90
index 357bb7f10..35b05f7b8 100644
--- a/linsolve/impl/psb_zrgmres.f90
+++ b/linsolve/impl/psb_zrgmres.f90
@@ -102,13 +102,15 @@
 !                                            stopped when  |r| <= eps * (|a||x|+|b|)
 !                                         2: err =  |r|/|b|; here the iteration is
 !                                            stopped when  |r| <= eps * |b|
+!                                         3: Same as 2 but with X and B scaled
+!                                            by s1 and s2 
 !                                         where r is the (preconditioned, recursive
 !                                         estimate of) residual. 
 !    irst   -  integer(optional)          Input: restart parameter 
 !
 
 subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
-     & itmax,iter,err,itrace,irst,istop)
+     & itmax,iter,err,itrace,irst,istop,s1,s2)
   use psb_base_mod
   use psb_prec_mod
   use psb_z_linsolve_conv_mod
@@ -124,6 +126,7 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(in)        :: itmax, itrace, irst,istop
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   Real(psb_dpk_), Optional, Intent(out) :: err
+  type(psb_z_vect_type), intent(inout), optional   :: s1, s2
 ! =   local data
   complex(psb_dpk_), allocatable   :: aux(:)
   complex(psb_dpk_), allocatable   :: c(:), s(:), h(:,:), rs(:), rst(:)
@@ -268,9 +271,20 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
   select case(istop_)
   case(psb_istop_ani_)
     ani = psb_spnrmi(a,desc_a,info)
-    bni = psb_geamax(b,desc_a,info)
+    if (present(s1)) then
+      call psb_gemlt(zone,s1,b,zzero,v(1),desc_a,info)
+      bni = psb_geamax(v(1),desc_a,info)
+    else
+      bni = psb_geamax(b,desc_a,info)
+    end if
   case(psb_istop_bn2_)
-    bn2 = psb_genrm2(b,desc_a,info)    
+    if (present(s1)) then
+      call psb_gemlt(zone,s1,b,zzero,v(1),desc_a,info)
+      bn2 = psb_genrm2(v(1),desc_a,info)
+    else
+      bn2 = psb_genrm2(b,desc_a,info)
+    end if
+      
   case(psb_istop_rn2_abs_)
     ! do nothing
   case(psb_istop_rrn2_)
@@ -282,6 +296,7 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
     end if
     
     call psb_spmm(-zone,a,x,zone,v(1),desc_a,info,work=aux)
+    if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info)
     if (info /= psb_success_) then 
       info=psb_err_from_subroutine_non_ 
       call psb_errpush(info,name)
@@ -323,7 +338,8 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
       call psb_errpush(info,name)
       goto 9999
     end if
-
+    if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info)
+    
     rs(1) = psb_genrm2(v(1),desc_a,info)
     rs(2:) = zzero
     if (info /= psb_success_) then 
@@ -378,8 +394,14 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
     inner:  Do i=1,nl
       itx  = itx + 1
 
-      call prec%apply(v(i),w1,desc_a,info)
+      if (present(s2)) then
+        call psb_gediv(v(i),s2,w,desc_a,info)
+        call prec%apply(w,w1,desc_a,info)
+      else
+        call prec%apply(v(i),w1,desc_a,info)
+      end if
       call psb_spmm(zone,a,w1,zzero,w,desc_a,info,work=aux)
+      if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
       !
       call mgs(i,h,v,w,rs,c,s,desc_a,info)
 
@@ -391,10 +413,11 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
         !
         rst = rs
         call psb_geaxpby(zone,x,zzero,xt,desc_a,info)
-        call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info)
+        call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info,s2=s2)
 
         call psb_geaxpby(zone,b,zzero,w1,desc_a,info)
         call psb_spmm(-zone,a,xt,zone,w1,desc_a,info,work=aux)
+        if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
         rni = psb_geamax(w1,desc_a,info)
         xni = psb_geamax(xt,desc_a,info)
         errnum = rni
@@ -432,7 +455,8 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
           call psb_geaxpby(zone,xt,zzero,x,desc_a,info)
 ! =          x = xt 
         case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
-          call  rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info)          !
+          call  rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2)
+          !
 
         end select
 
@@ -452,7 +476,7 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
       call psb_geaxpby(zone,xt,zzero,x,desc_a,info)!      x = xt 
 
     case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
-      call  rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info)          !
+      call  rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2)          !
     end select
     
     if (itx >= itmax_) then 
@@ -523,11 +547,12 @@ contains
   ! Rebuild solution X from the space V using the factor
   ! stored in R
   !
-  subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info)
+  subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2)
     complex(psb_dpk_)   :: c(:), s(:), rs(:), h(:,:)
     type(psb_z_vect_type) :: v(:), w, w1, x
     type(psb_desc_type) :: desc_a
     class(psb_zprec_type) :: prec
+    type(psb_z_vect_type), intent(inout), optional   :: s2
     integer(psb_ipk_) :: info
     integer(psb_ipk_) :: k,n
 
@@ -539,12 +564,13 @@ contains
     if (debug_level >= psb_debug_ext_) &
          & write(debug_unit,*) me,' ',trim(name),&
          & ' Rebuild x-> RS:',rs(1:n)
-    call w1%zero()
+    call w%zero()
     do k=1, n
-      call psb_geaxpby(rs(k),v(k),zone,w1,desc_a,info)
+      call psb_geaxpby(rs(k),v(k),zone,w,desc_a,info)
     end do
-    call prec%apply(w1,w,desc_a,info)
-    call psb_geaxpby(zone,w,zone,x,desc_a,info)
+    if (present(s2)) call psb_gediv(s2,w,desc_a,info)
+    call prec%apply(w,w1,desc_a,info)
+    call psb_geaxpby(zone,w1,zone,x,desc_a,info)
   end subroutine rebuildx
 
 end subroutine psb_zrgmres_vect
diff --git a/linsolve/psb_c_linsolve_conv_mod.f90 b/linsolve/psb_c_linsolve_conv_mod.f90
index 474acfaa6..840410813 100644
--- a/linsolve/psb_c_linsolve_conv_mod.f90
+++ b/linsolve/psb_c_linsolve_conv_mod.f90
@@ -83,15 +83,18 @@ contains
     stopdat%controls(psb_ik_itmax_) = itmax
 
     select case(stopdat%controls(psb_ik_stopc_))
-    case (1) 
+    case (psb_istop_ani_) 
       stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info)
       if (info == psb_success_)&
            & stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info)
 
-    case (2) 
+    case (psb_istop_bn2_) 
       stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info)
-
-    case (3)
+      
+    case (psb_istop_rn2_abs_) 
+      ! Do nothing
+      
+    case (psb_istop_rrn2_)
       call psb_geall(r,desc_a,info)
       call psb_geaxpby(cone,b,czero,r,desc_a,info)
       call psb_spmm(-cone,a,x,cone,r,desc_a,info)
@@ -108,8 +111,8 @@ contains
     end if
 
     stopdat%values(psb_ik_eps_)    = eps
-    stopdat%values(psb_ik_errnum_) = dzero
-    stopdat%values(psb_ik_errden_) = done
+    stopdat%values(psb_ik_errnum_) = szero
+    stopdat%values(psb_ik_errden_) = sone
 
     if ((stopdat%controls(psb_ik_trace_) > 0).and. (me == 0))&
          &  call log_header(methdname) 
@@ -123,7 +126,6 @@ contains
 
   end subroutine psb_c_init_conv
 
-
   function psb_c_check_conv(methdname,it,x,r,desc_a,stopdat,info) result(res)
     use psb_base_mod
     implicit none 
@@ -149,19 +151,26 @@ contains
     res = .false. 
 
     select case(stopdat%controls(psb_ik_stopc_)) 
-    case(1)
+    case(psb_istop_ani_)
       stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info)
-      if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
+      if (info == psb_success_) &
+           & stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
       stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_)
       stopdat%values(psb_ik_errden_) =&
            & (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)&
            &  +stopdat%values(psb_ik_bni_))
-    case(2)
+      
+    case(psb_istop_bn2_)
       stopdat%values(psb_ik_rn2_)  = psb_genrm2(r,desc_a,info)
       stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
       stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_)
 
-    case(3)
+    case (psb_istop_rn2_abs_) 
+      stopdat%values(psb_ik_rn2_abs_) = psb_genrm2(r,desc_a,info)
+      stopdat%values(psb_ik_errnum_)  = stopdat%values(psb_ik_rn2_abs_)
+      stopdat%values(psb_ik_errden_)  = sone
+
+    case(psb_istop_rrn2_)
       stopdat%values(psb_ik_rn2_)  = psb_genrm2(r,desc_a,info)
       stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
       stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_r0n2_)
@@ -201,8 +210,8 @@ contains
 
   end function psb_c_check_conv
 
-
-  subroutine psb_c_init_conv_vect(methdname,stopc,trace,itmax,a,x,b,eps,desc_a,stopdat,info)
+  subroutine psb_c_init_conv_vect(methdname,stopc,trace,itmax,&
+       & a,x,b,eps,desc_a,stopdat,info,s1,s2)
     use psb_base_mod
     implicit none 
     character(len=*), intent(in)      :: methdname
@@ -213,6 +222,7 @@ contains
     type(psb_desc_type), intent(in)   :: desc_a
     type(psb_itconv_type)             :: stopdat
     integer(psb_ipk_), intent(out)              :: info
+    type(psb_c_vect_type), optional :: s1, s2
 
     type(psb_ctxt_type) :: ctxt
     integer(psb_ipk_) :: me, np, err_act
@@ -236,15 +246,18 @@ contains
     stopdat%controls(psb_ik_itmax_) = itmax
 
     select case(stopdat%controls(psb_ik_stopc_))
-    case (1) 
+    case (psb_istop_ani_) 
       stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info)
       if (info == psb_success_)&
            & stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info)
 
-    case (2) 
+    case (psb_istop_bn2_) 
       stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info)
+      
+    case (psb_istop_rn2_abs_) 
+      ! Do nothing
 
-    case (3)
+    case (psb_istop_rrn2_)
       call psb_geasb(r,desc_a,info,scratch=.true.)
       call psb_geaxpby(cone,b,czero,r,desc_a,info)
       call psb_spmm(-cone,a,x,cone,r,desc_a,info)
@@ -261,8 +274,8 @@ contains
     end if
 
     stopdat%values(psb_ik_eps_)    = eps
-    stopdat%values(psb_ik_errnum_) = dzero
-    stopdat%values(psb_ik_errden_) = done
+    stopdat%values(psb_ik_errnum_) = szero
+    stopdat%values(psb_ik_errden_) = sone
 
     if ((stopdat%controls(psb_ik_trace_) > 0).and. (me == 0))&
          &  call log_header(methdname) 
@@ -276,7 +289,8 @@ contains
 
   end subroutine psb_c_init_conv_vect
 
-  function psb_c_check_conv_vect(methdname,it,x,r,desc_a,stopdat,info) result(res)
+  function psb_c_check_conv_vect(methdname,it,x,r,&
+       & desc_a,stopdat,info,s1,s2) result(res)
     use psb_base_mod
     implicit none 
     character(len=*), intent(in)     :: methdname
@@ -286,6 +300,7 @@ contains
     type(psb_itconv_type)            :: stopdat
     logical                          :: res
     integer(psb_ipk_), intent(out)             :: info
+    type(psb_c_vect_type), optional :: s1, s2
 
     type(psb_ctxt_type) :: ctxt
     integer(psb_ipk_) :: me, np, err_act
@@ -303,19 +318,26 @@ contains
 
 
     select case(stopdat%controls(psb_ik_stopc_)) 
-    case(1)
+    case(psb_istop_ani_)
       stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info)
-      if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
+      if (info == psb_success_) &
+           & stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
       stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_)
       stopdat%values(psb_ik_errden_) = &
            & (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)&
            & +stopdat%values(psb_ik_bni_))
-    case(2)
+      
+    case(psb_istop_bn2_)
       stopdat%values(psb_ik_rn2_)    = psb_genrm2(r,desc_a,info)
       stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
       stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_)
 
-    case(3)
+    case (psb_istop_rn2_abs_) 
+      stopdat%values(psb_ik_rn2_abs_) = psb_genrm2(r,desc_a,info)
+      stopdat%values(psb_ik_errnum_)  = stopdat%values(psb_ik_rn2_abs_)
+      stopdat%values(psb_ik_errden_)  = sone
+
+    case(psb_istop_rrn2_)
       stopdat%values(psb_ik_rn2_)  = psb_genrm2(r,desc_a,info)
       stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
       stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_r0n2_)
diff --git a/linsolve/psb_d_linsolve_conv_mod.f90 b/linsolve/psb_d_linsolve_conv_mod.f90
index f92edce39..64e5c3010 100644
--- a/linsolve/psb_d_linsolve_conv_mod.f90
+++ b/linsolve/psb_d_linsolve_conv_mod.f90
@@ -83,15 +83,18 @@ contains
     stopdat%controls(psb_ik_itmax_) = itmax
 
     select case(stopdat%controls(psb_ik_stopc_))
-    case (1) 
+    case (psb_istop_ani_) 
       stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info)
       if (info == psb_success_)&
            & stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info)
 
-    case (2) 
+    case (psb_istop_bn2_) 
       stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info)
-
-    case (3)
+      
+    case (psb_istop_rn2_abs_) 
+      ! Do nothing
+      
+    case (psb_istop_rrn2_)
       call psb_geall(r,desc_a,info)
       call psb_geaxpby(done,b,dzero,r,desc_a,info)
       call psb_spmm(-done,a,x,done,r,desc_a,info)
@@ -123,7 +126,6 @@ contains
 
   end subroutine psb_d_init_conv
 
-
   function psb_d_check_conv(methdname,it,x,r,desc_a,stopdat,info) result(res)
     use psb_base_mod
     implicit none 
@@ -149,19 +151,26 @@ contains
     res = .false. 
 
     select case(stopdat%controls(psb_ik_stopc_)) 
-    case(1)
+    case(psb_istop_ani_)
       stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info)
-      if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
+      if (info == psb_success_) &
+           & stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
       stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_)
       stopdat%values(psb_ik_errden_) =&
            & (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)&
            &  +stopdat%values(psb_ik_bni_))
-    case(2)
+      
+    case(psb_istop_bn2_)
       stopdat%values(psb_ik_rn2_)  = psb_genrm2(r,desc_a,info)
       stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
       stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_)
 
-    case(3)
+    case (psb_istop_rn2_abs_) 
+      stopdat%values(psb_ik_rn2_abs_) = psb_genrm2(r,desc_a,info)
+      stopdat%values(psb_ik_errnum_)  = stopdat%values(psb_ik_rn2_abs_)
+      stopdat%values(psb_ik_errden_)  = done
+
+    case(psb_istop_rrn2_)
       stopdat%values(psb_ik_rn2_)  = psb_genrm2(r,desc_a,info)
       stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
       stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_r0n2_)
@@ -201,8 +210,8 @@ contains
 
   end function psb_d_check_conv
 
-
-  subroutine psb_d_init_conv_vect(methdname,stopc,trace,itmax,a,x,b,eps,desc_a,stopdat,info)
+  subroutine psb_d_init_conv_vect(methdname,stopc,trace,itmax,&
+       & a,x,b,eps,desc_a,stopdat,info,s1,s2)
     use psb_base_mod
     implicit none 
     character(len=*), intent(in)      :: methdname
@@ -213,6 +222,7 @@ contains
     type(psb_desc_type), intent(in)   :: desc_a
     type(psb_itconv_type)             :: stopdat
     integer(psb_ipk_), intent(out)              :: info
+    type(psb_d_vect_type), optional :: s1, s2
 
     type(psb_ctxt_type) :: ctxt
     integer(psb_ipk_) :: me, np, err_act
@@ -236,15 +246,18 @@ contains
     stopdat%controls(psb_ik_itmax_) = itmax
 
     select case(stopdat%controls(psb_ik_stopc_))
-    case (1) 
+    case (psb_istop_ani_) 
       stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info)
       if (info == psb_success_)&
            & stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info)
 
-    case (2) 
+    case (psb_istop_bn2_) 
       stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info)
+      
+    case (psb_istop_rn2_abs_) 
+      ! Do nothing
 
-    case (3)
+    case (psb_istop_rrn2_)
       call psb_geasb(r,desc_a,info,scratch=.true.)
       call psb_geaxpby(done,b,dzero,r,desc_a,info)
       call psb_spmm(-done,a,x,done,r,desc_a,info)
@@ -276,7 +289,8 @@ contains
 
   end subroutine psb_d_init_conv_vect
 
-  function psb_d_check_conv_vect(methdname,it,x,r,desc_a,stopdat,info) result(res)
+  function psb_d_check_conv_vect(methdname,it,x,r,&
+       & desc_a,stopdat,info,s1,s2) result(res)
     use psb_base_mod
     implicit none 
     character(len=*), intent(in)     :: methdname
@@ -286,6 +300,7 @@ contains
     type(psb_itconv_type)            :: stopdat
     logical                          :: res
     integer(psb_ipk_), intent(out)             :: info
+    type(psb_d_vect_type), optional :: s1, s2
 
     type(psb_ctxt_type) :: ctxt
     integer(psb_ipk_) :: me, np, err_act
@@ -303,19 +318,26 @@ contains
 
 
     select case(stopdat%controls(psb_ik_stopc_)) 
-    case(1)
+    case(psb_istop_ani_)
       stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info)
-      if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
+      if (info == psb_success_) &
+           & stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
       stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_)
       stopdat%values(psb_ik_errden_) = &
            & (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)&
            & +stopdat%values(psb_ik_bni_))
-    case(2)
+      
+    case(psb_istop_bn2_)
       stopdat%values(psb_ik_rn2_)    = psb_genrm2(r,desc_a,info)
       stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
       stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_)
 
-    case(3)
+    case (psb_istop_rn2_abs_) 
+      stopdat%values(psb_ik_rn2_abs_) = psb_genrm2(r,desc_a,info)
+      stopdat%values(psb_ik_errnum_)  = stopdat%values(psb_ik_rn2_abs_)
+      stopdat%values(psb_ik_errden_)  = done
+
+    case(psb_istop_rrn2_)
       stopdat%values(psb_ik_rn2_)  = psb_genrm2(r,desc_a,info)
       stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
       stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_r0n2_)
diff --git a/linsolve/psb_linsolve_mod.f90 b/linsolve/psb_linsolve_mod.f90
index 1cd324622..7d72a9342 100644
--- a/linsolve/psb_linsolve_mod.f90
+++ b/linsolve/psb_linsolve_mod.f90
@@ -36,179 +36,9 @@
 Module psb_linsolve_mod
 
   use psb_const_mod
-  public 
-
-  interface psb_krylov
-    
-    Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
-         & itmax,iter,err,itrace,irst,istop,cond)
-      
-      use psb_base_mod, only  : psb_ipk_, psb_desc_type, psb_sspmat_type, &
-           & psb_spk_, psb_s_vect_type
-      use psb_prec_mod, only : psb_sprec_type
-      
-      character(len=*)                      :: method
-      Type(psb_sspmat_type), Intent(in)     :: a
-      Type(psb_desc_type), Intent(in)       :: desc_a
-      class(psb_sprec_type), intent(inout)  :: prec 
-      type(psb_s_vect_type), Intent(inout)  :: b
-      type(psb_s_vect_type), Intent(inout)  :: x
-      Real(psb_spk_), Intent(in)            :: eps
-      integer(psb_ipk_), intent(out)                  :: info
-      integer(psb_ipk_), Optional, Intent(in)         :: itmax, itrace, irst,istop
-      integer(psb_ipk_), Optional, Intent(out)        :: iter
-      Real(psb_spk_), Optional, Intent(out) :: err,cond
-
-    end Subroutine psb_skrylov_vect
-
-    Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
-         & itmax,iter,err,itrace,irst,istop,cond)
-      
-      use psb_base_mod, only  : psb_ipk_, psb_desc_type, psb_cspmat_type, &
-           & psb_spk_, psb_c_vect_type
-      use psb_prec_mod, only : psb_cprec_type
-      
-      character(len=*)                      :: method
-      Type(psb_cspmat_type), Intent(in)     :: a
-      Type(psb_desc_type), Intent(in)       :: desc_a
-      class(psb_cprec_type), intent(inout)  :: prec 
-      type(psb_c_vect_type), Intent(inout)  :: b
-      type(psb_c_vect_type), Intent(inout)  :: x
-      Real(psb_spk_), Intent(in)            :: eps
-      integer(psb_ipk_), intent(out)                  :: info
-      integer(psb_ipk_), Optional, Intent(in)         :: itmax, itrace, irst,istop
-      integer(psb_ipk_), Optional, Intent(out)        :: iter
-      Real(psb_spk_), Optional, Intent(out) :: err,cond
-
-    end Subroutine psb_ckrylov_vect
-
-    Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
-         & itmax,iter,err,itrace,irst,istop,cond)
-      
-      use psb_base_mod, only  : psb_ipk_, psb_desc_type, psb_dspmat_type, &
-           & psb_dpk_, psb_d_vect_type
-      use psb_prec_mod, only : psb_dprec_type
-      
-      character(len=*)                      :: method
-      Type(psb_dspmat_type), Intent(in)     :: a
-      Type(psb_desc_type), Intent(in)       :: desc_a
-      class(psb_dprec_type), intent(inout)  :: prec 
-      type(psb_d_vect_type), Intent(inout)  :: b
-      type(psb_d_vect_type), Intent(inout)  :: x
-      Real(psb_dpk_), Intent(in)            :: eps
-      integer(psb_ipk_), intent(out)                  :: info
-      integer(psb_ipk_), Optional, Intent(in)         :: itmax, itrace, irst,istop
-      integer(psb_ipk_), Optional, Intent(out)        :: iter
-      Real(psb_dpk_), Optional, Intent(out) :: err,cond
-
-    end Subroutine psb_dkrylov_vect
-
-    Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
-         & itmax,iter,err,itrace,irst,istop,cond)
-      
-      use psb_base_mod, only  : psb_ipk_, psb_desc_type, psb_zspmat_type, &
-           & psb_dpk_, psb_z_vect_type
-      use psb_prec_mod, only : psb_zprec_type
-      
-      character(len=*)                      :: method
-      Type(psb_zspmat_type), Intent(in)     :: a
-      Type(psb_desc_type), Intent(in)       :: desc_a
-      class(psb_zprec_type), intent(inout)  :: prec 
-      type(psb_z_vect_type), Intent(inout)  :: b
-      type(psb_z_vect_type), Intent(inout)  :: x
-      Real(psb_dpk_), Intent(in)            :: eps
-      integer(psb_ipk_), intent(out)                  :: info
-      integer(psb_ipk_), Optional, Intent(in)         :: itmax, itrace, irst,istop
-      integer(psb_ipk_), Optional, Intent(out)        :: iter
-      Real(psb_dpk_), Optional, Intent(out) :: err,cond
-
-    end Subroutine psb_zkrylov_vect
-
-  end interface
-
-
-  interface psb_richardson
-    
-    Subroutine psb_srichardson_vect(a,prec,b,x,eps,desc_a,info,&
-         & itmax,iter,err,itrace,istop)
-      
-      use psb_base_mod, only  : psb_ipk_, psb_desc_type, psb_sspmat_type, &
-           & psb_spk_, psb_s_vect_type
-      use psb_prec_mod, only : psb_sprec_type
-      
-      Type(psb_sspmat_type), Intent(in)     :: a
-      Type(psb_desc_type), Intent(in)       :: desc_a
-      class(psb_sprec_type), intent(inout)  :: prec 
-      type(psb_s_vect_type), Intent(inout)  :: b
-      type(psb_s_vect_type), Intent(inout)  :: x
-      Real(psb_spk_), Intent(in)            :: eps
-      integer(psb_ipk_), intent(out)                  :: info
-      integer(psb_ipk_), Optional, Intent(in)         :: itmax, itrace, istop
-      integer(psb_ipk_), Optional, Intent(out)        :: iter
-      Real(psb_spk_), Optional, Intent(out) :: err
-
-    end Subroutine psb_srichardson_vect
-
-    Subroutine psb_crichardson_vect(a,prec,b,x,eps,desc_a,info,&
-         & itmax,iter,err,itrace,istop)
-      
-      use psb_base_mod, only  : psb_ipk_, psb_desc_type, psb_cspmat_type, &
-           & psb_spk_, psb_c_vect_type
-      use psb_prec_mod, only : psb_cprec_type
-      
-      Type(psb_cspmat_type), Intent(in)     :: a
-      Type(psb_desc_type), Intent(in)       :: desc_a
-      class(psb_cprec_type), intent(inout)  :: prec 
-      type(psb_c_vect_type), Intent(inout)  :: b
-      type(psb_c_vect_type), Intent(inout)  :: x
-      Real(psb_spk_), Intent(in)            :: eps
-      integer(psb_ipk_), intent(out)                  :: info
-      integer(psb_ipk_), Optional, Intent(in)         :: itmax, itrace, istop
-      integer(psb_ipk_), Optional, Intent(out)        :: iter
-      Real(psb_spk_), Optional, Intent(out) :: err
-
-    end Subroutine psb_crichardson_vect
-
-    Subroutine psb_drichardson_vect(a,prec,b,x,eps,desc_a,info,&
-         & itmax,iter,err,itrace,istop)
-      
-      use psb_base_mod, only  : psb_ipk_, psb_desc_type, psb_dspmat_type, &
-           & psb_dpk_, psb_d_vect_type
-      use psb_prec_mod, only : psb_dprec_type
-      
-      Type(psb_dspmat_type), Intent(in)     :: a
-      Type(psb_desc_type), Intent(in)       :: desc_a
-      class(psb_dprec_type), intent(inout)  :: prec 
-      type(psb_d_vect_type), Intent(inout)  :: b
-      type(psb_d_vect_type), Intent(inout)  :: x
-      Real(psb_dpk_), Intent(in)            :: eps
-      integer(psb_ipk_), intent(out)                  :: info
-      integer(psb_ipk_), Optional, Intent(in)         :: itmax, itrace, istop
-      integer(psb_ipk_), Optional, Intent(out)        :: iter
-      Real(psb_dpk_), Optional, Intent(out) :: err
-
-    end Subroutine psb_drichardson_vect
-
-    Subroutine psb_zrichardson_vect(a,prec,b,x,eps,desc_a,info,&
-         & itmax,iter,err,itrace,istop)
-      
-      use psb_base_mod, only  : psb_ipk_, psb_desc_type, psb_zspmat_type, &
-           & psb_dpk_, psb_z_vect_type
-      use psb_prec_mod, only : psb_zprec_type
-      
-      Type(psb_zspmat_type), Intent(in)     :: a
-      Type(psb_desc_type), Intent(in)       :: desc_a
-      class(psb_zprec_type), intent(inout)  :: prec 
-      type(psb_z_vect_type), Intent(inout)  :: b
-      type(psb_z_vect_type), Intent(inout)  :: x
-      Real(psb_dpk_), Intent(in)            :: eps
-      integer(psb_ipk_), intent(out)                  :: info
-      integer(psb_ipk_), Optional, Intent(in)         :: itmax, itrace, istop
-      integer(psb_ipk_), Optional, Intent(out)        :: iter
-      Real(psb_dpk_), Optional, Intent(out) :: err
-
-    end Subroutine psb_zrichardson_vect
-
-  end interface
+  use psb_s_linsolve_mod
+  use psb_d_linsolve_mod
+  use psb_c_linsolve_mod
+  use psb_z_linsolve_mod
   
 end module psb_linsolve_mod
diff --git a/linsolve/psb_s_linsolve_conv_mod.f90 b/linsolve/psb_s_linsolve_conv_mod.f90
index 8cac66a41..3d7028052 100644
--- a/linsolve/psb_s_linsolve_conv_mod.f90
+++ b/linsolve/psb_s_linsolve_conv_mod.f90
@@ -83,15 +83,18 @@ contains
     stopdat%controls(psb_ik_itmax_) = itmax
 
     select case(stopdat%controls(psb_ik_stopc_))
-    case (1) 
+    case (psb_istop_ani_) 
       stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info)
       if (info == psb_success_)&
            & stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info)
 
-    case (2) 
+    case (psb_istop_bn2_) 
       stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info)
-
-    case (3)
+      
+    case (psb_istop_rn2_abs_) 
+      ! Do nothing
+      
+    case (psb_istop_rrn2_)
       call psb_geall(r,desc_a,info)
       call psb_geaxpby(sone,b,szero,r,desc_a,info)
       call psb_spmm(-sone,a,x,sone,r,desc_a,info)
@@ -108,8 +111,8 @@ contains
     end if
 
     stopdat%values(psb_ik_eps_)    = eps
-    stopdat%values(psb_ik_errnum_) = dzero
-    stopdat%values(psb_ik_errden_) = done
+    stopdat%values(psb_ik_errnum_) = szero
+    stopdat%values(psb_ik_errden_) = sone
 
     if ((stopdat%controls(psb_ik_trace_) > 0).and. (me == 0))&
          &  call log_header(methdname) 
@@ -123,7 +126,6 @@ contains
 
   end subroutine psb_s_init_conv
 
-
   function psb_s_check_conv(methdname,it,x,r,desc_a,stopdat,info) result(res)
     use psb_base_mod
     implicit none 
@@ -149,19 +151,26 @@ contains
     res = .false. 
 
     select case(stopdat%controls(psb_ik_stopc_)) 
-    case(1)
+    case(psb_istop_ani_)
       stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info)
-      if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
+      if (info == psb_success_) &
+           & stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
       stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_)
       stopdat%values(psb_ik_errden_) =&
            & (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)&
            &  +stopdat%values(psb_ik_bni_))
-    case(2)
+      
+    case(psb_istop_bn2_)
       stopdat%values(psb_ik_rn2_)  = psb_genrm2(r,desc_a,info)
       stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
       stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_)
 
-    case(3)
+    case (psb_istop_rn2_abs_) 
+      stopdat%values(psb_ik_rn2_abs_) = psb_genrm2(r,desc_a,info)
+      stopdat%values(psb_ik_errnum_)  = stopdat%values(psb_ik_rn2_abs_)
+      stopdat%values(psb_ik_errden_)  = sone
+
+    case(psb_istop_rrn2_)
       stopdat%values(psb_ik_rn2_)  = psb_genrm2(r,desc_a,info)
       stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
       stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_r0n2_)
@@ -201,8 +210,8 @@ contains
 
   end function psb_s_check_conv
 
-
-  subroutine psb_s_init_conv_vect(methdname,stopc,trace,itmax,a,x,b,eps,desc_a,stopdat,info)
+  subroutine psb_s_init_conv_vect(methdname,stopc,trace,itmax,&
+       & a,x,b,eps,desc_a,stopdat,info,s1,s2)
     use psb_base_mod
     implicit none 
     character(len=*), intent(in)      :: methdname
@@ -213,6 +222,7 @@ contains
     type(psb_desc_type), intent(in)   :: desc_a
     type(psb_itconv_type)             :: stopdat
     integer(psb_ipk_), intent(out)              :: info
+    type(psb_s_vect_type), optional :: s1, s2
 
     type(psb_ctxt_type) :: ctxt
     integer(psb_ipk_) :: me, np, err_act
@@ -236,15 +246,18 @@ contains
     stopdat%controls(psb_ik_itmax_) = itmax
 
     select case(stopdat%controls(psb_ik_stopc_))
-    case (1) 
+    case (psb_istop_ani_) 
       stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info)
       if (info == psb_success_)&
            & stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info)
 
-    case (2) 
+    case (psb_istop_bn2_) 
       stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info)
+      
+    case (psb_istop_rn2_abs_) 
+      ! Do nothing
 
-    case (3)
+    case (psb_istop_rrn2_)
       call psb_geasb(r,desc_a,info,scratch=.true.)
       call psb_geaxpby(sone,b,szero,r,desc_a,info)
       call psb_spmm(-sone,a,x,sone,r,desc_a,info)
@@ -261,8 +274,8 @@ contains
     end if
 
     stopdat%values(psb_ik_eps_)    = eps
-    stopdat%values(psb_ik_errnum_) = dzero
-    stopdat%values(psb_ik_errden_) = done
+    stopdat%values(psb_ik_errnum_) = szero
+    stopdat%values(psb_ik_errden_) = sone
 
     if ((stopdat%controls(psb_ik_trace_) > 0).and. (me == 0))&
          &  call log_header(methdname) 
@@ -276,7 +289,8 @@ contains
 
   end subroutine psb_s_init_conv_vect
 
-  function psb_s_check_conv_vect(methdname,it,x,r,desc_a,stopdat,info) result(res)
+  function psb_s_check_conv_vect(methdname,it,x,r,&
+       & desc_a,stopdat,info,s1,s2) result(res)
     use psb_base_mod
     implicit none 
     character(len=*), intent(in)     :: methdname
@@ -286,6 +300,7 @@ contains
     type(psb_itconv_type)            :: stopdat
     logical                          :: res
     integer(psb_ipk_), intent(out)             :: info
+    type(psb_s_vect_type), optional :: s1, s2
 
     type(psb_ctxt_type) :: ctxt
     integer(psb_ipk_) :: me, np, err_act
@@ -303,19 +318,26 @@ contains
 
 
     select case(stopdat%controls(psb_ik_stopc_)) 
-    case(1)
+    case(psb_istop_ani_)
       stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info)
-      if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
+      if (info == psb_success_) &
+           & stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
       stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_)
       stopdat%values(psb_ik_errden_) = &
            & (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)&
            & +stopdat%values(psb_ik_bni_))
-    case(2)
+      
+    case(psb_istop_bn2_)
       stopdat%values(psb_ik_rn2_)    = psb_genrm2(r,desc_a,info)
       stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
       stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_)
 
-    case(3)
+    case (psb_istop_rn2_abs_) 
+      stopdat%values(psb_ik_rn2_abs_) = psb_genrm2(r,desc_a,info)
+      stopdat%values(psb_ik_errnum_)  = stopdat%values(psb_ik_rn2_abs_)
+      stopdat%values(psb_ik_errden_)  = sone
+
+    case(psb_istop_rrn2_)
       stopdat%values(psb_ik_rn2_)  = psb_genrm2(r,desc_a,info)
       stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
       stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_r0n2_)
diff --git a/linsolve/psb_z_linsolve_conv_mod.f90 b/linsolve/psb_z_linsolve_conv_mod.f90
index d6082262c..efe60f425 100644
--- a/linsolve/psb_z_linsolve_conv_mod.f90
+++ b/linsolve/psb_z_linsolve_conv_mod.f90
@@ -83,15 +83,18 @@ contains
     stopdat%controls(psb_ik_itmax_) = itmax
 
     select case(stopdat%controls(psb_ik_stopc_))
-    case (1) 
+    case (psb_istop_ani_) 
       stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info)
       if (info == psb_success_)&
            & stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info)
 
-    case (2) 
+    case (psb_istop_bn2_) 
       stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info)
-
-    case (3)
+      
+    case (psb_istop_rn2_abs_) 
+      ! Do nothing
+      
+    case (psb_istop_rrn2_)
       call psb_geall(r,desc_a,info)
       call psb_geaxpby(zone,b,zzero,r,desc_a,info)
       call psb_spmm(-zone,a,x,zone,r,desc_a,info)
@@ -123,7 +126,6 @@ contains
 
   end subroutine psb_z_init_conv
 
-
   function psb_z_check_conv(methdname,it,x,r,desc_a,stopdat,info) result(res)
     use psb_base_mod
     implicit none 
@@ -149,19 +151,26 @@ contains
     res = .false. 
 
     select case(stopdat%controls(psb_ik_stopc_)) 
-    case(1)
+    case(psb_istop_ani_)
       stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info)
-      if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
+      if (info == psb_success_) &
+           & stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
       stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_)
       stopdat%values(psb_ik_errden_) =&
            & (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)&
            &  +stopdat%values(psb_ik_bni_))
-    case(2)
+      
+    case(psb_istop_bn2_)
       stopdat%values(psb_ik_rn2_)  = psb_genrm2(r,desc_a,info)
       stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
       stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_)
 
-    case(3)
+    case (psb_istop_rn2_abs_) 
+      stopdat%values(psb_ik_rn2_abs_) = psb_genrm2(r,desc_a,info)
+      stopdat%values(psb_ik_errnum_)  = stopdat%values(psb_ik_rn2_abs_)
+      stopdat%values(psb_ik_errden_)  = done
+
+    case(psb_istop_rrn2_)
       stopdat%values(psb_ik_rn2_)  = psb_genrm2(r,desc_a,info)
       stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
       stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_r0n2_)
@@ -201,8 +210,8 @@ contains
 
   end function psb_z_check_conv
 
-
-  subroutine psb_z_init_conv_vect(methdname,stopc,trace,itmax,a,x,b,eps,desc_a,stopdat,info)
+  subroutine psb_z_init_conv_vect(methdname,stopc,trace,itmax,&
+       & a,x,b,eps,desc_a,stopdat,info,s1,s2)
     use psb_base_mod
     implicit none 
     character(len=*), intent(in)      :: methdname
@@ -213,6 +222,7 @@ contains
     type(psb_desc_type), intent(in)   :: desc_a
     type(psb_itconv_type)             :: stopdat
     integer(psb_ipk_), intent(out)              :: info
+    type(psb_z_vect_type), optional :: s1, s2
 
     type(psb_ctxt_type) :: ctxt
     integer(psb_ipk_) :: me, np, err_act
@@ -236,15 +246,18 @@ contains
     stopdat%controls(psb_ik_itmax_) = itmax
 
     select case(stopdat%controls(psb_ik_stopc_))
-    case (1) 
+    case (psb_istop_ani_) 
       stopdat%values(psb_ik_ani_) = psb_spnrmi(a,desc_a,info)
       if (info == psb_success_)&
            & stopdat%values(psb_ik_bni_) = psb_geamax(b,desc_a,info)
 
-    case (2) 
+    case (psb_istop_bn2_) 
       stopdat%values(psb_ik_bn2_) = psb_genrm2(b,desc_a,info)
+      
+    case (psb_istop_rn2_abs_) 
+      ! Do nothing
 
-    case (3)
+    case (psb_istop_rrn2_)
       call psb_geasb(r,desc_a,info,scratch=.true.)
       call psb_geaxpby(zone,b,zzero,r,desc_a,info)
       call psb_spmm(-zone,a,x,zone,r,desc_a,info)
@@ -276,7 +289,8 @@ contains
 
   end subroutine psb_z_init_conv_vect
 
-  function psb_z_check_conv_vect(methdname,it,x,r,desc_a,stopdat,info) result(res)
+  function psb_z_check_conv_vect(methdname,it,x,r,&
+       & desc_a,stopdat,info,s1,s2) result(res)
     use psb_base_mod
     implicit none 
     character(len=*), intent(in)     :: methdname
@@ -286,6 +300,7 @@ contains
     type(psb_itconv_type)            :: stopdat
     logical                          :: res
     integer(psb_ipk_), intent(out)             :: info
+    type(psb_z_vect_type), optional :: s1, s2
 
     type(psb_ctxt_type) :: ctxt
     integer(psb_ipk_) :: me, np, err_act
@@ -303,19 +318,26 @@ contains
 
 
     select case(stopdat%controls(psb_ik_stopc_)) 
-    case(1)
+    case(psb_istop_ani_)
       stopdat%values(psb_ik_rni_) = psb_geamax(r,desc_a,info)
-      if (info == psb_success_) stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
+      if (info == psb_success_) &
+           & stopdat%values(psb_ik_xni_) = psb_geamax(x,desc_a,info)
       stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rni_)
       stopdat%values(psb_ik_errden_) = &
            & (stopdat%values(psb_ik_ani_)*stopdat%values(psb_ik_xni_)&
            & +stopdat%values(psb_ik_bni_))
-    case(2)
+      
+    case(psb_istop_bn2_)
       stopdat%values(psb_ik_rn2_)    = psb_genrm2(r,desc_a,info)
       stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
       stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_bn2_)
 
-    case(3)
+    case (psb_istop_rn2_abs_) 
+      stopdat%values(psb_ik_rn2_abs_) = psb_genrm2(r,desc_a,info)
+      stopdat%values(psb_ik_errnum_)  = stopdat%values(psb_ik_rn2_abs_)
+      stopdat%values(psb_ik_errden_)  = done
+
+    case(psb_istop_rrn2_)
       stopdat%values(psb_ik_rn2_)  = psb_genrm2(r,desc_a,info)
       stopdat%values(psb_ik_errnum_) = stopdat%values(psb_ik_rn2_)
       stopdat%values(psb_ik_errden_) = stopdat%values(psb_ik_r0n2_)