diff --git a/linsolve/impl/psb_crgmres.f90 b/linsolve/impl/psb_crgmres.f90
index 6efc672a4..cf6e8dc67 100644
--- a/linsolve/impl/psb_crgmres.f90
+++ b/linsolve/impl/psb_crgmres.f90
@@ -134,7 +134,7 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
   type(psb_c_vect_type)              :: w, w1, xt
   real(psb_spk_) :: tmp 
   complex(psb_spk_) :: scal, gm, rti, rti1
-  integer(psb_ipk_) ::litmax, naux, it, k, itrace_,&
+  integer(psb_ipk_) ::itmax_, naux, it, k, itrace_,&
        & n_row, n_col, nl
   integer(psb_lpk_) :: mglob
   Logical, Parameter :: exchange=.True., noexchange=.False., use_srot=.true.
@@ -202,9 +202,9 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
   end select
 
   if (present(itmax)) then 
-    litmax = itmax
+    itmax_ = itmax
   else
-    litmax = 1000
+    itmax_ = 1000
   endif
 
   if (present(itrace)) then
@@ -369,7 +369,7 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
       goto 9999
     end if
     
-    if ((errnum <= eps*errden).or.(itx >= litmax)) exit restart  
+    if ((errnum <= eps*errden).or.(itx >= itmax_)) exit restart  
 
     if (itrace_ > 0) &
          & call log_conv(methdname,me,itx,itrace_,errnum,errden,deps)
@@ -391,26 +391,8 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
       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)
 
-      do k = 1, i
-        h(k,i) = psb_gedot(v(k),w,desc_a,info)
-        call psb_geaxpby(-h(k,i),v(k),cone,w,desc_a,info)
-      end do
-      h(i+1,i) = psb_genrm2(w,desc_a,info)
-      scal=cone/h(i+1,i)
-      call psb_geaxpby(scal,w,czero,v(i+1),desc_a,info)
-      do k=2,i
-        call crot(1,h(k-1,i),1,h(k,i),1,real(c(k-1),kind=psb_spk_),s(k-1))
-      enddo
-
-
-      rti  = h(i,i)
-      rti1 = h(i+1,i) 
-      call crotg(rti,rti1,tmp,s(i))
-      c(i) = cmplx(tmp,szero)
-      call crot(1,h(i,i),1,h(i+1,i),1,real(c(i),kind=psb_spk_),s(i))
-      h(i+1,i) = czero
-      call crot(1,rs(i),1,rs(i+1),1,real(c(i),kind=psb_spk_),s(i))
       
       select case(istop_)
       case(psb_istop_ani_)
@@ -418,16 +400,9 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
         ! build x and then compute the residual and its infinity norm
         !
         rst = rs
-        call w1%set(czero)
-        call ctrsm('l','u','n','n',i,1,cone,h,size(h,1),rst,size(rst,1))
-        if (debug_level >= psb_debug_ext_) &
-             & write(debug_unit,*) me,' ',trim(name),&
-             & ' Rebuild x-> RS:',rst(1:i)
-        do k=1, i
-          call psb_geaxpby(rst(k),v(k),cone,xt,desc_a,info)
-        end do
-        call prec%apply(xt,desc_a,info)
-        call psb_geaxpby(cone,x,cone,xt,desc_a,info)
+        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 psb_geaxpby(cone,b,czero,w1,desc_a,info)
         call psb_spmm(-cone,a,xt,cone,w1,desc_a,info,work=aux)
         rni = psb_geamax(w1,desc_a,info)
@@ -466,35 +441,9 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
         case(psb_istop_ani_)
           call psb_geaxpby(cone,xt,czero,x,desc_a,info)
 ! =          x = xt 
-        case(psb_istop_bn2_, psb_istop_rn2_abs_)
-          !
-          ! build x
-          !
-          call ctrsm('l','u','n','n',i,1,cone,h,size(h,1),rs,size(rs,1))
-          if (debug_level >= psb_debug_ext_) &
-               & write(debug_unit,*) me,' ',trim(name),&
-               & ' Rebuild x-> RS:',rs(1:i)
-          call w1%set(czero)
-          do k=1, i
-            call psb_geaxpby(rs(k),v(k),cone,w1,desc_a,info)
-          end do
-          call prec%apply(w1,w,desc_a,info)
-          call psb_geaxpby(cone,w,cone,x,desc_a,info)
-
-        case(psb_istop_rrn2_)
-          !
-          ! build x
-          !
-          call ctrsm('l','u','n','n',i,1,cone,h,size(h,1),rs,size(rs,1))
-          if (debug_level >= psb_debug_ext_) &
-               & write(debug_unit,*) me,' ',trim(name),&
-               & ' Rebuild x-> RS:',rs(1:i)
-          call w1%set(czero)
-          do k=1, i
-            call psb_geaxpby(rs(k),v(k),cone,w1,desc_a,info)
-          end do
-          call prec%apply(w1,w,desc_a,info)
-          call psb_geaxpby(cone,w,cone,x,desc_a,info)          
+        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)          !
+
         end select
 
         if (itrace_ > 0) &
@@ -512,38 +461,11 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
     case(psb_istop_ani_)
       call psb_geaxpby(cone,xt,czero,x,desc_a,info)!      x = xt 
 
-    case(psb_istop_bn2_, psb_istop_rn2_abs_)
-      !
-      ! build x
-      !
-      call ctrsm('l','u','n','n',nl,1,cone,h,size(h,1),rs,size(rs,1))
-      if (debug_level >= psb_debug_ext_) &
-           & write(debug_unit,*) me,' ',trim(name),&
-           & ' Rebuild x-> RS:',rs(1:nl)
-      call w1%set(czero)
-      do k=1, nl
-        call psb_geaxpby(rs(k),v(k),cone,w1,desc_a,info)
-      end do
-      call prec%apply(w1,w,desc_a,info)
-      call psb_geaxpby(cone,w,cone,x,desc_a,info)
-      
-    case(psb_istop_rrn2_)
-      !
-      ! build x
-      !
-      call ctrsm('l','u','n','n',nl,1,cone,h,size(h,1),rs,size(rs,1))
-      if (debug_level >= psb_debug_ext_) &
-           & write(debug_unit,*) me,' ',trim(name),&
-           & ' Rebuild x-> RS:',rs(1:nl)
-      call w1%set(czero)
-      do k=1, nl
-        call psb_geaxpby(rs(k),v(k),cone,w1,desc_a,info)
-      end do
-      call prec%apply(w1,w,desc_a,info)
-      call psb_geaxpby(cone,w,cone,x,desc_a,info)
+    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)          !
     end select
     
-    if (itx >= litmax) then 
+    if (itx >= itmax_) then 
       if (itrace_ > 0) then 
         if (mod(itx,itrace_)/=0) &
              &  call log_conv(methdname,me,itx,ione,errnum,errden,deps)
@@ -573,6 +495,67 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
 
 9999 call psb_error_handler(err_act)
   return
+contains
+  !
+  ! Advance one step the Gram-Schmidt process, and apply
+  ! Givens rotations to keep H triangular
+  !
+  subroutine mgs(n,h,v,w,rs,c,s,desc_a,info)
+    complex(psb_spk_)   :: c(:), s(:), h(:,:), rs(:)
+    type(psb_c_vect_type) :: v(:), w
+    type(psb_desc_type) :: desc_a
+    complex(psb_spk_) :: scal, gm, rti, rti1
+    real(psb_spk_) :: tmp
+    integer(psb_ipk_) :: info
+    integer(psb_ipk_) :: k,n
+    !
+    
+    do k = 1, n
+      h(k,n) = psb_gedot(v(k),w,desc_a,info)
+      call psb_geaxpby(-h(k,n),v(k),cone,w,desc_a,info)
+    end do
+    h(n+1,n) = psb_genrm2(w,desc_a,info)
+    scal=cone/h(n+1,n)
+    call psb_geaxpby(scal,w,czero,v(n+1),desc_a,info)
+    do k=2,n
+      call crot(1,h(k-1:k-1,n),1,h(k:k,n),1,real(c(k-1),kind=psb_spk_),s(k-1))
+    enddo
+    
+    rti  = h(n,n)
+    rti1 = h(n+1,n) 
+    call crotg(rti,rti1,tmp,s(n))
+    c(n) = cmplx(tmp,szero)
+    call crot(1,h(n:n,n),1,h(n+1:n+1,n),1,real(c(n),kind=psb_spk_),s(n))
+    call crot(1,rs(n:n),1,rs(n+1:n+1),1,real(c(n),kind=psb_spk_),s(n))
+  end subroutine mgs
+
+  !
+  ! 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)
+    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
+    integer(psb_ipk_) :: info
+    integer(psb_ipk_) :: k,n
+
+    !
+    !
+    ! build x
+    !
+    call ctrsm('l','u','n','n',n,1,cone,h,size(h,1),rs,size(rs,1))
+    if (debug_level >= psb_debug_ext_) &
+         & write(debug_unit,*) me,' ',trim(name),&
+         & ' Rebuild x-> RS:',rs(1:n)
+    call w1%zero()
+    do k=1, n
+      call psb_geaxpby(rs(k),v(k),cone,w1,desc_a,info)
+    end do
+    call prec%apply(w1,w,desc_a,info)
+    call psb_geaxpby(cone,w,cone,x,desc_a,info)
+  end subroutine rebuildx
 
 end subroutine psb_crgmres_vect
 
diff --git a/linsolve/impl/psb_drgmres.f90 b/linsolve/impl/psb_drgmres.f90
index 2b16f1e9f..7cd4bb47b 100644
--- a/linsolve/impl/psb_drgmres.f90
+++ b/linsolve/impl/psb_drgmres.f90
@@ -134,7 +134,7 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
   type(psb_d_vect_type)              :: w, w1, xt
   real(psb_dpk_) :: tmp 
   real(psb_dpk_) :: scal, gm, rti, rti1
-  integer(psb_ipk_) ::litmax, naux, it, k, itrace_,&
+  integer(psb_ipk_) ::itmax_, naux, it, k, itrace_,&
        & n_row, n_col, nl
   integer(psb_lpk_) :: mglob
   Logical, Parameter :: exchange=.True., noexchange=.False., use_srot=.true.
@@ -202,9 +202,9 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
   end select
 
   if (present(itmax)) then 
-    litmax = itmax
+    itmax_ = itmax
   else
-    litmax = 1000
+    itmax_ = 1000
   endif
 
   if (present(itrace)) then
@@ -369,7 +369,7 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
       goto 9999
     end if
     
-    if ((errnum <= eps*errden).or.(itx >= litmax)) exit restart  
+    if ((errnum <= eps*errden).or.(itx >= itmax_)) exit restart  
 
     if (itrace_ > 0) &
          & call log_conv(methdname,me,itx,itrace_,errnum,errden,deps)
@@ -391,26 +391,8 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
       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)
 
-      do k = 1, i
-        h(k,i) = psb_gedot(v(k),w,desc_a,info)
-        call psb_geaxpby(-h(k,i),v(k),done,w,desc_a,info)
-      end do
-      h(i+1,i) = psb_genrm2(w,desc_a,info)
-      scal=done/h(i+1,i)
-      call psb_geaxpby(scal,w,dzero,v(i+1),desc_a,info)
-      do k=2,i
-        call drot(1,h(k-1,i),1,h(k,i),1,real(c(k-1),kind=psb_dpk_),s(k-1))
-      enddo
-
-
-      rti  = h(i,i)
-      rti1 = h(i+1,i) 
-      call drotg(rti,rti1,tmp,s(i))
-      c(i) = cmplx(tmp,dzero)
-      call drot(1,h(i,i),1,h(i+1,i),1,real(c(i),kind=psb_dpk_),s(i))
-      h(i+1,i) = dzero
-      call drot(1,rs(i),1,rs(i+1),1,real(c(i),kind=psb_dpk_),s(i))
       
       select case(istop_)
       case(psb_istop_ani_)
@@ -418,16 +400,9 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
         ! build x and then compute the residual and its infinity norm
         !
         rst = rs
-        call w1%set(dzero)
-        call dtrsm('l','u','n','n',i,1,done,h,size(h,1),rst,size(rst,1))
-        if (debug_level >= psb_debug_ext_) &
-             & write(debug_unit,*) me,' ',trim(name),&
-             & ' Rebuild x-> RS:',rst(1:i)
-        do k=1, i
-          call psb_geaxpby(rst(k),v(k),done,xt,desc_a,info)
-        end do
-        call prec%apply(xt,desc_a,info)
-        call psb_geaxpby(done,x,done,xt,desc_a,info)
+        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 psb_geaxpby(done,b,dzero,w1,desc_a,info)
         call psb_spmm(-done,a,xt,done,w1,desc_a,info,work=aux)
         rni = psb_geamax(w1,desc_a,info)
@@ -466,35 +441,9 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
         case(psb_istop_ani_)
           call psb_geaxpby(done,xt,dzero,x,desc_a,info)
 ! =          x = xt 
-        case(psb_istop_bn2_, psb_istop_rn2_abs_)
-          !
-          ! build x
-          !
-          call dtrsm('l','u','n','n',i,1,done,h,size(h,1),rs,size(rs,1))
-          if (debug_level >= psb_debug_ext_) &
-               & write(debug_unit,*) me,' ',trim(name),&
-               & ' Rebuild x-> RS:',rs(1:i)
-          call w1%set(dzero)
-          do k=1, i
-            call psb_geaxpby(rs(k),v(k),done,w1,desc_a,info)
-          end do
-          call prec%apply(w1,w,desc_a,info)
-          call psb_geaxpby(done,w,done,x,desc_a,info)
-
-        case(psb_istop_rrn2_)
-          !
-          ! build x
-          !
-          call dtrsm('l','u','n','n',i,1,done,h,size(h,1),rs,size(rs,1))
-          if (debug_level >= psb_debug_ext_) &
-               & write(debug_unit,*) me,' ',trim(name),&
-               & ' Rebuild x-> RS:',rs(1:i)
-          call w1%set(dzero)
-          do k=1, i
-            call psb_geaxpby(rs(k),v(k),done,w1,desc_a,info)
-          end do
-          call prec%apply(w1,w,desc_a,info)
-          call psb_geaxpby(done,w,done,x,desc_a,info)          
+        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)          !
+
         end select
 
         if (itrace_ > 0) &
@@ -512,38 +461,11 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
     case(psb_istop_ani_)
       call psb_geaxpby(done,xt,dzero,x,desc_a,info)!      x = xt 
 
-    case(psb_istop_bn2_, psb_istop_rn2_abs_)
-      !
-      ! build x
-      !
-      call dtrsm('l','u','n','n',nl,1,done,h,size(h,1),rs,size(rs,1))
-      if (debug_level >= psb_debug_ext_) &
-           & write(debug_unit,*) me,' ',trim(name),&
-           & ' Rebuild x-> RS:',rs(1:nl)
-      call w1%set(dzero)
-      do k=1, nl
-        call psb_geaxpby(rs(k),v(k),done,w1,desc_a,info)
-      end do
-      call prec%apply(w1,w,desc_a,info)
-      call psb_geaxpby(done,w,done,x,desc_a,info)
-      
-    case(psb_istop_rrn2_)
-      !
-      ! build x
-      !
-      call dtrsm('l','u','n','n',nl,1,done,h,size(h,1),rs,size(rs,1))
-      if (debug_level >= psb_debug_ext_) &
-           & write(debug_unit,*) me,' ',trim(name),&
-           & ' Rebuild x-> RS:',rs(1:nl)
-      call w1%set(dzero)
-      do k=1, nl
-        call psb_geaxpby(rs(k),v(k),done,w1,desc_a,info)
-      end do
-      call prec%apply(w1,w,desc_a,info)
-      call psb_geaxpby(done,w,done,x,desc_a,info)
+    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)          !
     end select
     
-    if (itx >= litmax) then 
+    if (itx >= itmax_) then 
       if (itrace_ > 0) then 
         if (mod(itx,itrace_)/=0) &
              &  call log_conv(methdname,me,itx,ione,errnum,errden,deps)
@@ -573,6 +495,67 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
 
 9999 call psb_error_handler(err_act)
   return
+contains
+  !
+  ! Advance one step the Gram-Schmidt process, and apply
+  ! Givens rotations to keep H triangular
+  !
+  subroutine mgs(n,h,v,w,rs,c,s,desc_a,info)
+    real(psb_dpk_)   :: c(:), s(:), h(:,:), rs(:)
+    type(psb_d_vect_type) :: v(:), w
+    type(psb_desc_type) :: desc_a
+    real(psb_dpk_) :: scal, gm, rti, rti1
+    real(psb_dpk_) :: tmp
+    integer(psb_ipk_) :: info
+    integer(psb_ipk_) :: k,n
+    !
+    
+    do k = 1, n
+      h(k,n) = psb_gedot(v(k),w,desc_a,info)
+      call psb_geaxpby(-h(k,n),v(k),done,w,desc_a,info)
+    end do
+    h(n+1,n) = psb_genrm2(w,desc_a,info)
+    scal=done/h(n+1,n)
+    call psb_geaxpby(scal,w,dzero,v(n+1),desc_a,info)
+    do k=2,n
+      call drot(1,h(k-1:k-1,n),1,h(k:k,n),1,real(c(k-1),kind=psb_dpk_),s(k-1))
+    enddo
+    
+    rti  = h(n,n)
+    rti1 = h(n+1,n) 
+    call drotg(rti,rti1,tmp,s(n))
+    c(n) = cmplx(tmp,dzero)
+    call drot(1,h(n:n,n),1,h(n+1:n+1,n),1,real(c(n),kind=psb_dpk_),s(n))
+    call drot(1,rs(n:n),1,rs(n+1:n+1),1,real(c(n),kind=psb_dpk_),s(n))
+  end subroutine mgs
+
+  !
+  ! 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)
+    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
+    integer(psb_ipk_) :: info
+    integer(psb_ipk_) :: k,n
+
+    !
+    !
+    ! build x
+    !
+    call dtrsm('l','u','n','n',n,1,done,h,size(h,1),rs,size(rs,1))
+    if (debug_level >= psb_debug_ext_) &
+         & write(debug_unit,*) me,' ',trim(name),&
+         & ' Rebuild x-> RS:',rs(1:n)
+    call w1%zero()
+    do k=1, n
+      call psb_geaxpby(rs(k),v(k),done,w1,desc_a,info)
+    end do
+    call prec%apply(w1,w,desc_a,info)
+    call psb_geaxpby(done,w,done,x,desc_a,info)
+  end subroutine rebuildx
 
 end subroutine psb_drgmres_vect
 
diff --git a/linsolve/impl/psb_srgmres.f90 b/linsolve/impl/psb_srgmres.f90
index 998d5b714..fc2dfca11 100644
--- a/linsolve/impl/psb_srgmres.f90
+++ b/linsolve/impl/psb_srgmres.f90
@@ -134,7 +134,7 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
   type(psb_s_vect_type)              :: w, w1, xt
   real(psb_spk_) :: tmp 
   real(psb_spk_) :: scal, gm, rti, rti1
-  integer(psb_ipk_) ::litmax, naux, it, k, itrace_,&
+  integer(psb_ipk_) ::itmax_, naux, it, k, itrace_,&
        & n_row, n_col, nl
   integer(psb_lpk_) :: mglob
   Logical, Parameter :: exchange=.True., noexchange=.False., use_srot=.true.
@@ -202,9 +202,9 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
   end select
 
   if (present(itmax)) then 
-    litmax = itmax
+    itmax_ = itmax
   else
-    litmax = 1000
+    itmax_ = 1000
   endif
 
   if (present(itrace)) then
@@ -369,7 +369,7 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
       goto 9999
     end if
     
-    if ((errnum <= eps*errden).or.(itx >= litmax)) exit restart  
+    if ((errnum <= eps*errden).or.(itx >= itmax_)) exit restart  
 
     if (itrace_ > 0) &
          & call log_conv(methdname,me,itx,itrace_,errnum,errden,deps)
@@ -391,26 +391,8 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
       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)
 
-      do k = 1, i
-        h(k,i) = psb_gedot(v(k),w,desc_a,info)
-        call psb_geaxpby(-h(k,i),v(k),sone,w,desc_a,info)
-      end do
-      h(i+1,i) = psb_genrm2(w,desc_a,info)
-      scal=sone/h(i+1,i)
-      call psb_geaxpby(scal,w,szero,v(i+1),desc_a,info)
-      do k=2,i
-        call srot(1,h(k-1,i),1,h(k,i),1,real(c(k-1),kind=psb_spk_),s(k-1))
-      enddo
-
-
-      rti  = h(i,i)
-      rti1 = h(i+1,i) 
-      call srotg(rti,rti1,tmp,s(i))
-      c(i) = cmplx(tmp,szero)
-      call srot(1,h(i,i),1,h(i+1,i),1,real(c(i),kind=psb_spk_),s(i))
-      h(i+1,i) = szero
-      call srot(1,rs(i),1,rs(i+1),1,real(c(i),kind=psb_spk_),s(i))
       
       select case(istop_)
       case(psb_istop_ani_)
@@ -418,16 +400,9 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
         ! build x and then compute the residual and its infinity norm
         !
         rst = rs
-        call w1%set(szero)
-        call strsm('l','u','n','n',i,1,sone,h,size(h,1),rst,size(rst,1))
-        if (debug_level >= psb_debug_ext_) &
-             & write(debug_unit,*) me,' ',trim(name),&
-             & ' Rebuild x-> RS:',rst(1:i)
-        do k=1, i
-          call psb_geaxpby(rst(k),v(k),sone,xt,desc_a,info)
-        end do
-        call prec%apply(xt,desc_a,info)
-        call psb_geaxpby(sone,x,sone,xt,desc_a,info)
+        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 psb_geaxpby(sone,b,szero,w1,desc_a,info)
         call psb_spmm(-sone,a,xt,sone,w1,desc_a,info,work=aux)
         rni = psb_geamax(w1,desc_a,info)
@@ -466,35 +441,9 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
         case(psb_istop_ani_)
           call psb_geaxpby(sone,xt,szero,x,desc_a,info)
 ! =          x = xt 
-        case(psb_istop_bn2_, psb_istop_rn2_abs_)
-          !
-          ! build x
-          !
-          call strsm('l','u','n','n',i,1,sone,h,size(h,1),rs,size(rs,1))
-          if (debug_level >= psb_debug_ext_) &
-               & write(debug_unit,*) me,' ',trim(name),&
-               & ' Rebuild x-> RS:',rs(1:i)
-          call w1%set(szero)
-          do k=1, i
-            call psb_geaxpby(rs(k),v(k),sone,w1,desc_a,info)
-          end do
-          call prec%apply(w1,w,desc_a,info)
-          call psb_geaxpby(sone,w,sone,x,desc_a,info)
-
-        case(psb_istop_rrn2_)
-          !
-          ! build x
-          !
-          call strsm('l','u','n','n',i,1,sone,h,size(h,1),rs,size(rs,1))
-          if (debug_level >= psb_debug_ext_) &
-               & write(debug_unit,*) me,' ',trim(name),&
-               & ' Rebuild x-> RS:',rs(1:i)
-          call w1%set(szero)
-          do k=1, i
-            call psb_geaxpby(rs(k),v(k),sone,w1,desc_a,info)
-          end do
-          call prec%apply(w1,w,desc_a,info)
-          call psb_geaxpby(sone,w,sone,x,desc_a,info)          
+        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)          !
+
         end select
 
         if (itrace_ > 0) &
@@ -512,38 +461,11 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
     case(psb_istop_ani_)
       call psb_geaxpby(sone,xt,szero,x,desc_a,info)!      x = xt 
 
-    case(psb_istop_bn2_, psb_istop_rn2_abs_)
-      !
-      ! build x
-      !
-      call strsm('l','u','n','n',nl,1,sone,h,size(h,1),rs,size(rs,1))
-      if (debug_level >= psb_debug_ext_) &
-           & write(debug_unit,*) me,' ',trim(name),&
-           & ' Rebuild x-> RS:',rs(1:nl)
-      call w1%set(szero)
-      do k=1, nl
-        call psb_geaxpby(rs(k),v(k),sone,w1,desc_a,info)
-      end do
-      call prec%apply(w1,w,desc_a,info)
-      call psb_geaxpby(sone,w,sone,x,desc_a,info)
-      
-    case(psb_istop_rrn2_)
-      !
-      ! build x
-      !
-      call strsm('l','u','n','n',nl,1,sone,h,size(h,1),rs,size(rs,1))
-      if (debug_level >= psb_debug_ext_) &
-           & write(debug_unit,*) me,' ',trim(name),&
-           & ' Rebuild x-> RS:',rs(1:nl)
-      call w1%set(szero)
-      do k=1, nl
-        call psb_geaxpby(rs(k),v(k),sone,w1,desc_a,info)
-      end do
-      call prec%apply(w1,w,desc_a,info)
-      call psb_geaxpby(sone,w,sone,x,desc_a,info)
+    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)          !
     end select
     
-    if (itx >= litmax) then 
+    if (itx >= itmax_) then 
       if (itrace_ > 0) then 
         if (mod(itx,itrace_)/=0) &
              &  call log_conv(methdname,me,itx,ione,errnum,errden,deps)
@@ -573,6 +495,67 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
 
 9999 call psb_error_handler(err_act)
   return
+contains
+  !
+  ! Advance one step the Gram-Schmidt process, and apply
+  ! Givens rotations to keep H triangular
+  !
+  subroutine mgs(n,h,v,w,rs,c,s,desc_a,info)
+    real(psb_spk_)   :: c(:), s(:), h(:,:), rs(:)
+    type(psb_s_vect_type) :: v(:), w
+    type(psb_desc_type) :: desc_a
+    real(psb_spk_) :: scal, gm, rti, rti1
+    real(psb_spk_) :: tmp
+    integer(psb_ipk_) :: info
+    integer(psb_ipk_) :: k,n
+    !
+    
+    do k = 1, n
+      h(k,n) = psb_gedot(v(k),w,desc_a,info)
+      call psb_geaxpby(-h(k,n),v(k),sone,w,desc_a,info)
+    end do
+    h(n+1,n) = psb_genrm2(w,desc_a,info)
+    scal=sone/h(n+1,n)
+    call psb_geaxpby(scal,w,szero,v(n+1),desc_a,info)
+    do k=2,n
+      call srot(1,h(k-1:k-1,n),1,h(k:k,n),1,real(c(k-1),kind=psb_spk_),s(k-1))
+    enddo
+    
+    rti  = h(n,n)
+    rti1 = h(n+1,n) 
+    call srotg(rti,rti1,tmp,s(n))
+    c(n) = cmplx(tmp,szero)
+    call srot(1,h(n:n,n),1,h(n+1:n+1,n),1,real(c(n),kind=psb_spk_),s(n))
+    call srot(1,rs(n:n),1,rs(n+1:n+1),1,real(c(n),kind=psb_spk_),s(n))
+  end subroutine mgs
+
+  !
+  ! 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)
+    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
+    integer(psb_ipk_) :: info
+    integer(psb_ipk_) :: k,n
+
+    !
+    !
+    ! build x
+    !
+    call strsm('l','u','n','n',n,1,sone,h,size(h,1),rs,size(rs,1))
+    if (debug_level >= psb_debug_ext_) &
+         & write(debug_unit,*) me,' ',trim(name),&
+         & ' Rebuild x-> RS:',rs(1:n)
+    call w1%zero()
+    do k=1, n
+      call psb_geaxpby(rs(k),v(k),sone,w1,desc_a,info)
+    end do
+    call prec%apply(w1,w,desc_a,info)
+    call psb_geaxpby(sone,w,sone,x,desc_a,info)
+  end subroutine rebuildx
 
 end subroutine psb_srgmres_vect
 
diff --git a/linsolve/impl/psb_zrgmres.f90 b/linsolve/impl/psb_zrgmres.f90
index 0b1f25624..63ca823d1 100644
--- a/linsolve/impl/psb_zrgmres.f90
+++ b/linsolve/impl/psb_zrgmres.f90
@@ -134,7 +134,7 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
   type(psb_z_vect_type)              :: w, w1, xt
   real(psb_dpk_) :: tmp 
   complex(psb_dpk_) :: scal, gm, rti, rti1
-  integer(psb_ipk_) ::litmax, naux, it, k, itrace_,&
+  integer(psb_ipk_) ::itmax_, naux, it, k, itrace_,&
        & n_row, n_col, nl
   integer(psb_lpk_) :: mglob
   Logical, Parameter :: exchange=.True., noexchange=.False., use_srot=.true.
@@ -202,9 +202,9 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
   end select
 
   if (present(itmax)) then 
-    litmax = itmax
+    itmax_ = itmax
   else
-    litmax = 1000
+    itmax_ = 1000
   endif
 
   if (present(itrace)) then
@@ -369,7 +369,7 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
       goto 9999
     end if
     
-    if ((errnum <= eps*errden).or.(itx >= litmax)) exit restart  
+    if ((errnum <= eps*errden).or.(itx >= itmax_)) exit restart  
 
     if (itrace_ > 0) &
          & call log_conv(methdname,me,itx,itrace_,errnum,errden,deps)
@@ -391,26 +391,8 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
       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)
 
-      do k = 1, i
-        h(k,i) = psb_gedot(v(k),w,desc_a,info)
-        call psb_geaxpby(-h(k,i),v(k),zone,w,desc_a,info)
-      end do
-      h(i+1,i) = psb_genrm2(w,desc_a,info)
-      scal=zone/h(i+1,i)
-      call psb_geaxpby(scal,w,zzero,v(i+1),desc_a,info)
-      do k=2,i
-        call zrot(1,h(k-1,i),1,h(k,i),1,real(c(k-1),kind=psb_dpk_),s(k-1))
-      enddo
-
-
-      rti  = h(i,i)
-      rti1 = h(i+1,i) 
-      call zrotg(rti,rti1,tmp,s(i))
-      c(i) = cmplx(tmp,dzero)
-      call zrot(1,h(i,i),1,h(i+1,i),1,real(c(i),kind=psb_dpk_),s(i))
-      h(i+1,i) = zzero
-      call zrot(1,rs(i),1,rs(i+1),1,real(c(i),kind=psb_dpk_),s(i))
       
       select case(istop_)
       case(psb_istop_ani_)
@@ -418,16 +400,9 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
         ! build x and then compute the residual and its infinity norm
         !
         rst = rs
-        call w1%set(zzero)
-        call ztrsm('l','u','n','n',i,1,zone,h,size(h,1),rst,size(rst,1))
-        if (debug_level >= psb_debug_ext_) &
-             & write(debug_unit,*) me,' ',trim(name),&
-             & ' Rebuild x-> RS:',rst(1:i)
-        do k=1, i
-          call psb_geaxpby(rst(k),v(k),zone,xt,desc_a,info)
-        end do
-        call prec%apply(xt,desc_a,info)
-        call psb_geaxpby(zone,x,zone,xt,desc_a,info)
+        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 psb_geaxpby(zone,b,zzero,w1,desc_a,info)
         call psb_spmm(-zone,a,xt,zone,w1,desc_a,info,work=aux)
         rni = psb_geamax(w1,desc_a,info)
@@ -466,35 +441,9 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
         case(psb_istop_ani_)
           call psb_geaxpby(zone,xt,zzero,x,desc_a,info)
 ! =          x = xt 
-        case(psb_istop_bn2_, psb_istop_rn2_abs_)
-          !
-          ! build x
-          !
-          call ztrsm('l','u','n','n',i,1,zone,h,size(h,1),rs,size(rs,1))
-          if (debug_level >= psb_debug_ext_) &
-               & write(debug_unit,*) me,' ',trim(name),&
-               & ' Rebuild x-> RS:',rs(1:i)
-          call w1%set(zzero)
-          do k=1, i
-            call psb_geaxpby(rs(k),v(k),zone,w1,desc_a,info)
-          end do
-          call prec%apply(w1,w,desc_a,info)
-          call psb_geaxpby(zone,w,zone,x,desc_a,info)
-
-        case(psb_istop_rrn2_)
-          !
-          ! build x
-          !
-          call ztrsm('l','u','n','n',i,1,zone,h,size(h,1),rs,size(rs,1))
-          if (debug_level >= psb_debug_ext_) &
-               & write(debug_unit,*) me,' ',trim(name),&
-               & ' Rebuild x-> RS:',rs(1:i)
-          call w1%set(zzero)
-          do k=1, i
-            call psb_geaxpby(rs(k),v(k),zone,w1,desc_a,info)
-          end do
-          call prec%apply(w1,w,desc_a,info)
-          call psb_geaxpby(zone,w,zone,x,desc_a,info)          
+        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)          !
+
         end select
 
         if (itrace_ > 0) &
@@ -512,38 +461,11 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
     case(psb_istop_ani_)
       call psb_geaxpby(zone,xt,zzero,x,desc_a,info)!      x = xt 
 
-    case(psb_istop_bn2_, psb_istop_rn2_abs_)
-      !
-      ! build x
-      !
-      call ztrsm('l','u','n','n',nl,1,zone,h,size(h,1),rs,size(rs,1))
-      if (debug_level >= psb_debug_ext_) &
-           & write(debug_unit,*) me,' ',trim(name),&
-           & ' Rebuild x-> RS:',rs(1:nl)
-      call w1%set(zzero)
-      do k=1, nl
-        call psb_geaxpby(rs(k),v(k),zone,w1,desc_a,info)
-      end do
-      call prec%apply(w1,w,desc_a,info)
-      call psb_geaxpby(zone,w,zone,x,desc_a,info)
-      
-    case(psb_istop_rrn2_)
-      !
-      ! build x
-      !
-      call ztrsm('l','u','n','n',nl,1,zone,h,size(h,1),rs,size(rs,1))
-      if (debug_level >= psb_debug_ext_) &
-           & write(debug_unit,*) me,' ',trim(name),&
-           & ' Rebuild x-> RS:',rs(1:nl)
-      call w1%set(zzero)
-      do k=1, nl
-        call psb_geaxpby(rs(k),v(k),zone,w1,desc_a,info)
-      end do
-      call prec%apply(w1,w,desc_a,info)
-      call psb_geaxpby(zone,w,zone,x,desc_a,info)
+    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)          !
     end select
     
-    if (itx >= litmax) then 
+    if (itx >= itmax_) then 
       if (itrace_ > 0) then 
         if (mod(itx,itrace_)/=0) &
              &  call log_conv(methdname,me,itx,ione,errnum,errden,deps)
@@ -573,6 +495,67 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
 
 9999 call psb_error_handler(err_act)
   return
+contains
+  !
+  ! Advance one step the Gram-Schmidt process, and apply
+  ! Givens rotations to keep H triangular
+  !
+  subroutine mgs(n,h,v,w,rs,c,s,desc_a,info)
+    complex(psb_dpk_)   :: c(:), s(:), h(:,:), rs(:)
+    type(psb_z_vect_type) :: v(:), w
+    type(psb_desc_type) :: desc_a
+    complex(psb_dpk_) :: scal, gm, rti, rti1
+    real(psb_dpk_) :: tmp
+    integer(psb_ipk_) :: info
+    integer(psb_ipk_) :: k,n
+    !
+    
+    do k = 1, n
+      h(k,n) = psb_gedot(v(k),w,desc_a,info)
+      call psb_geaxpby(-h(k,n),v(k),zone,w,desc_a,info)
+    end do
+    h(n+1,n) = psb_genrm2(w,desc_a,info)
+    scal=zone/h(n+1,n)
+    call psb_geaxpby(scal,w,zzero,v(n+1),desc_a,info)
+    do k=2,n
+      call zrot(1,h(k-1:k-1,n),1,h(k:k,n),1,real(c(k-1),kind=psb_dpk_),s(k-1))
+    enddo
+    
+    rti  = h(n,n)
+    rti1 = h(n+1,n) 
+    call zrotg(rti,rti1,tmp,s(n))
+    c(n) = cmplx(tmp,dzero)
+    call zrot(1,h(n:n,n),1,h(n+1:n+1,n),1,real(c(n),kind=psb_dpk_),s(n))
+    call zrot(1,rs(n:n),1,rs(n+1:n+1),1,real(c(n),kind=psb_dpk_),s(n))
+  end subroutine mgs
+
+  !
+  ! 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)
+    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
+    integer(psb_ipk_) :: info
+    integer(psb_ipk_) :: k,n
+
+    !
+    !
+    ! build x
+    !
+    call ztrsm('l','u','n','n',n,1,zone,h,size(h,1),rs,size(rs,1))
+    if (debug_level >= psb_debug_ext_) &
+         & write(debug_unit,*) me,' ',trim(name),&
+         & ' Rebuild x-> RS:',rs(1:n)
+    call w1%zero()
+    do k=1, n
+      call psb_geaxpby(rs(k),v(k),zone,w1,desc_a,info)
+    end do
+    call prec%apply(w1,w,desc_a,info)
+    call psb_geaxpby(zone,w,zone,x,desc_a,info)
+  end subroutine rebuildx
 
 end subroutine psb_zrgmres_vect
 
diff --git a/util/psb_d_mmio_impl.f90 b/util/psb_d_mmio_impl.f90
index 374eeed92..c5423f02d 100644
--- a/util/psb_d_mmio_impl.f90
+++ b/util/psb_d_mmio_impl.f90
@@ -39,8 +39,10 @@ subroutine mm_dvet_read(b, info, iunit, filename)
   integer(psb_ipk_) :: nrow, ncol, i,root, np,  me,  ircode, j, infile
   character            :: mmheader*15, fmt*15, object*10, type*10, sym*15,&
        & line*1024
+  logical :: opened
 
   info = psb_success_
+  opened = .false.
   if (present(filename)) then
     if (filename == '-') then 
       infile=5
@@ -51,6 +53,7 @@ subroutine mm_dvet_read(b, info, iunit, filename)
         infile=99
       endif
       open(infile,file=filename, status='OLD', err=901, action='READ')
+      opened = .true.
     endif
   else 
     if (present(iunit)) then 
@@ -83,7 +86,7 @@ subroutine mm_dvet_read(b, info, iunit, filename)
     end do
 
   end if      ! read right hand sides
-  if (infile /= 5) close(infile)
+  if (opened) close(infile)
 
   return 
   ! open failed
@@ -111,8 +114,10 @@ subroutine mm_dvet2_read(b, info, iunit, filename)
   integer(psb_ipk_) :: nrow, ncol, i,root, np,  me,  ircode, j, infile
   character            :: mmheader*15, fmt*15, object*10, type*10, sym*15,&
        & line*1024
+  logical :: opened
 
   info = psb_success_
+  opened = .false.
   if (present(filename)) then
     if (filename == '-') then 
       infile=5
@@ -123,6 +128,7 @@ subroutine mm_dvet2_read(b, info, iunit, filename)
         infile=99
       endif
       open(infile,file=filename, status='OLD', err=901, action='READ')
+      opened = .true.
     endif
   else 
     if (present(iunit)) then 
@@ -153,7 +159,7 @@ subroutine mm_dvet2_read(b, info, iunit, filename)
     read(infile,fmt=*,end=902) ((b(i,j), i=1,nrow),j=1,ncol)
 
   end if      ! read right hand sides
-  if (infile /= 5) close(infile)
+  if (opened) close(infile)
 
   return 
   ! open failed
@@ -182,8 +188,10 @@ subroutine mm_dvet2_write(b, header, info, iunit, filename)
   integer(psb_ipk_) :: nrow, ncol, i,root, np,  me,  ircode, j, outfile
 
   character(len=80)                 :: frmtv 
+  logical :: opened
 
   info = psb_success_
+  opened = .false.
   if (present(filename)) then
     if (filename == '-') then 
       outfile=6
@@ -194,6 +202,7 @@ subroutine mm_dvet2_write(b, header, info, iunit, filename)
         outfile=99
       endif
       open(outfile,file=filename, err=901, action='WRITE')
+      opened = .true.
     endif
   else 
     if (present(iunit)) then 
@@ -212,7 +221,7 @@ subroutine mm_dvet2_write(b, header, info, iunit, filename)
 
   write(outfile,fmt='(es26.18,1x)') ((b(i,j), i=1,nrow),j=1,ncol)
 
-  if (outfile /= 6) close(outfile)
+  if (opened) close(outfile)
 
   return 
   ! open failed
@@ -234,8 +243,10 @@ subroutine mm_dvet1_write(b, header, info, iunit, filename)
   integer(psb_ipk_) :: nrow, ncol, i,root, np,  me,  ircode, j, outfile
 
   character(len=80)                 :: frmtv 
+  logical :: opened
 
   info = psb_success_
+  opened = .false.
   if (present(filename)) then
     if (filename == '-') then 
       outfile=6
@@ -246,6 +257,7 @@ subroutine mm_dvet1_write(b, header, info, iunit, filename)
         outfile=99
       endif
       open(outfile,file=filename, err=901, action='WRITE')
+      opened = .true.
     endif
   else 
     if (present(iunit)) then 
@@ -268,7 +280,7 @@ subroutine mm_dvet1_write(b, header, info, iunit, filename)
     write(outfile,frmtv) b(i)
   end do
 
-  if (outfile /= 6) close(outfile)
+  if (opened) close(outfile)
 
   return 
   ! open failed
@@ -324,9 +336,10 @@ subroutine dmm_mat_read(a, info, iunit, filename)
   integer(psb_ipk_) :: nrow, ncol, nnzero
   integer(psb_ipk_) :: ircode, i,nzr,infile
   type(psb_d_coo_sparse_mat), allocatable :: acoo
+  logical :: opened
 
   info = psb_success_
-
+  opened = .false.
   if (present(filename)) then
     if (filename == '-') then 
       infile=5
@@ -337,6 +350,7 @@ subroutine dmm_mat_read(a, info, iunit, filename)
         infile=99
       endif
       open(infile,file=filename, status='OLD', err=901, action='READ')
+      opened = .true.
     endif
   else 
     if (present(iunit)) then 
@@ -421,7 +435,7 @@ subroutine dmm_mat_read(a, info, iunit, filename)
     call a%mv_from(acoo)
   end if
 
-  if (infile /= 5) close(infile)
+  if (opened) close(infile)
   return 
 
   ! open failed
@@ -449,10 +463,11 @@ subroutine dmm_mat_write(a,mtitle,info,iunit,filename)
   integer(psb_ipk_), optional, intent(in)          :: iunit
   character(len=*), optional, intent(in) :: filename
   integer(psb_ipk_) :: iout
+  logical :: opened
 
 
   info = psb_success_
-
+  opened = .false.
   if (present(filename)) then 
     if (filename == '-') then 
       iout=6
@@ -463,6 +478,7 @@ subroutine dmm_mat_write(a,mtitle,info,iunit,filename)
         iout=99
       endif
       open(iout,file=filename, err=901, action='WRITE')
+      opened = .true.
     endif
   else 
     if (present(iunit)) then 
@@ -474,7 +490,7 @@ subroutine dmm_mat_write(a,mtitle,info,iunit,filename)
 
   call a%print(iout,head=mtitle)
 
-  if (iout /= 6) close(iout)
+  if (opened) close(iout)
 
 
   return
@@ -498,9 +514,10 @@ subroutine ldmm_mat_read(a, info, iunit, filename)
   integer(psb_lpk_) :: nrow, ncol, nnzero, i, nzr
   integer(psb_ipk_) :: ircode, infile
   type(psb_ld_coo_sparse_mat), allocatable :: acoo
+  logical :: opened
 
   info = psb_success_
-
+  opened = .false.
   if (present(filename)) then
     if (filename == '-') then 
       infile=5
@@ -511,6 +528,7 @@ subroutine ldmm_mat_read(a, info, iunit, filename)
         infile=99
       endif
       open(infile,file=filename, status='OLD', err=901, action='READ')
+      opened = .true.
     endif
   else 
     if (present(iunit)) then 
@@ -595,7 +613,7 @@ subroutine ldmm_mat_read(a, info, iunit, filename)
     call a%mv_from(acoo)
   end if
 
-  if (infile /= 5) close(infile)
+  if (opened) close(infile)
   return 
 
   ! open failed
@@ -623,10 +641,10 @@ subroutine ldmm_mat_write(a,mtitle,info,iunit,filename)
   integer(psb_ipk_), optional, intent(in)          :: iunit
   character(len=*), optional, intent(in) :: filename
   integer(psb_ipk_) :: iout
-
+  logical :: opened
 
   info = psb_success_
-
+  opened = .false.
   if (present(filename)) then 
     if (filename == '-') then 
       iout=6
@@ -637,6 +655,7 @@ subroutine ldmm_mat_write(a,mtitle,info,iunit,filename)
         iout=99
       endif
       open(iout,file=filename, err=901, action='WRITE')
+      opened = .true.
     endif
   else 
     if (present(iunit)) then 
@@ -648,7 +667,7 @@ subroutine ldmm_mat_write(a,mtitle,info,iunit,filename)
 
   call a%print(iout,head=mtitle)
 
-  if (iout /= 6) close(iout)
+  if (opened) close(iout)
 
 
   return
@@ -658,5 +677,3 @@ subroutine ldmm_mat_write(a,mtitle,info,iunit,filename)
   write(psb_err_unit,*) 'Error while opening ',filename
   return
 end subroutine ldmm_mat_write
-
-