diff --git a/base/comm/psb_cgather.f90 b/base/comm/psb_cgather.f90
index 538fe4e17..6f2f6e902 100644
--- a/base/comm/psb_cgather.f90
+++ b/base/comm/psb_cgather.f90
@@ -58,8 +58,8 @@ subroutine  psb_cgather_vect(globx, locx, desc_a, info, iroot)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
-  integer(psb_ipk_) :: ierr(5), err_act, n, ilocx, iglobx, jlocx,&
-       & jglobx, lda_locx, lda_globx, m, k, jlx, ilx, i, idx
+  integer(psb_ipk_) :: ierr(5), err_act, jlx, ilx, lda_locx, lda_globx, i
+  integer(psb_lpk_) :: m, n, k, ilocx,  jlocx, idx, iglobx, jglobx
   complex(psb_spk_), allocatable :: llocx(:)
   character(len=20)        :: name, ch_err
 
@@ -171,8 +171,8 @@ subroutine  psb_cgather_multivect(globx, locx, desc_a, info, iroot)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
-  integer(psb_ipk_) :: ierr(5), err_act, n, ilocx, iglobx, jlocx,&
-       & jglobx, lda_locx, lda_globx, m, k, jlx, ilx, i, idx
+  integer(psb_ipk_) :: ierr(5), err_act, jlx, ilx, lda_locx, lda_globx, i 
+  integer(psb_lpk_) :: m, n, k, ilocx, jlocx, idx, iglobx, jglobx
   complex(psb_spk_), allocatable :: llocx(:,:)
   character(len=20)        :: name, ch_err
 
diff --git a/base/comm/psb_cgather_a.f90 b/base/comm/psb_cgather_a.f90
index de1949719..fcfd78eda 100644
--- a/base/comm/psb_cgather_a.f90
+++ b/base/comm/psb_cgather_a.f90
@@ -58,9 +58,9 @@ subroutine  psb_cgatherm(globx, locx, desc_a, info, iroot)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
-  integer(psb_ipk_) :: ierr(5), err_act, n, ilocx, iglobx, jlocx,&
-       & jglobx, lda_locx, lda_globx, m, lock, globk, maxk, k, jlx, &
-       & ilx, i, j, idx
+  integer(psb_ipk_) :: ierr(5), err_act, lda_locx, lda_globx, lock, globk,&
+       & maxk, k, jlx, ilx, i, j
+  integer(psb_lpk_) :: m, n, ilocx,  jlocx, idx, iglobx, jglobx
 
   character(len=20)        :: name, ch_err
 
@@ -230,8 +230,9 @@ subroutine  psb_cgatherv(globx, locx, desc_a, info, iroot)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
-  integer(psb_ipk_) :: ierr(5), err_act, n, ilocx, iglobx, jlocx,&
-       & jglobx, lda_locx, lda_globx, m, k, jlx, ilx, i, idx
+  integer(psb_ipk_) :: ierr(5), err_act, lda_locx, lda_globx, lock, globk,&
+       & maxk, k, jlx, ilx, i, j
+  integer(psb_lpk_) :: m, n, ilocx,  jlocx, idx, iglobx, jglobx
 
   character(len=20)        :: name, ch_err
 
diff --git a/base/comm/psb_chalo.f90 b/base/comm/psb_chalo.f90
index f245d51e2..07827aa7a 100644
--- a/base/comm/psb_chalo.f90
+++ b/base/comm/psb_chalo.f90
@@ -66,9 +66,9 @@ subroutine  psb_chalo_vect(x,desc_a,info,work,tran,mode,data)
   character, intent(in), optional         :: tran
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me,&
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, imode,&
-       & err, liwork,data_
+  integer(psb_ipk_) :: ictxt, np, me, err_act, iix, jjx, &
+       & nrow, imode, err, liwork,data_
+  integer(psb_lpk_) :: m, n, ix, ijx
   complex(psb_spk_),pointer :: iwork(:)
   character                 :: tran_
   character(len=20)         :: name, ch_err
@@ -119,7 +119,7 @@ subroutine  psb_chalo_vect(x,desc_a,info,work,tran,mode,data)
   endif
 
   ! check vector correctness
-  call psb_chkvect(m,ione,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
@@ -205,9 +205,9 @@ subroutine  psb_chalo_multivect(x,desc_a,info,work,tran,mode,data)
   character, intent(in), optional         :: tran
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me,&
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, imode,&
-       & err, liwork,data_
+  integer(psb_ipk_) :: ictxt, np, me, err_act, iix, jjx, &
+       & nrow, imode, err, liwork,data_
+  integer(psb_lpk_) :: m, n, ix, ijx
   complex(psb_spk_),pointer :: iwork(:)
   character                 :: tran_
   character(len=20)         :: name, ch_err
@@ -258,7 +258,7 @@ subroutine  psb_chalo_multivect(x,desc_a,info,work,tran,mode,data)
   endif
 
   ! check vector correctness
-  call psb_chkvect(m,ione,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
diff --git a/base/comm/psb_chalo_a.f90 b/base/comm/psb_chalo_a.f90
index e28902f91..3ccfcd378 100644
--- a/base/comm/psb_chalo_a.f90
+++ b/base/comm/psb_chalo_a.f90
@@ -66,8 +66,9 @@ subroutine  psb_chalom(x,desc_a,info,jx,ik,work,tran,mode,data)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me
-  integer(psb_ipk_) :: err_act, m, n, iix, jjx, ix, ijx, k, maxk, nrow, imode, i,&
+  integer(psb_ipk_) :: err_act, iix, jjx, k, maxk, nrow, imode, i,&
        & err, liwork,data_, ldx
+  integer(psb_lpk_) :: m, n, ix, ijx
   complex(psb_spk_),pointer :: iwork(:), xp(:,:)
   character                :: tran_
   character(len=20)        :: name, ch_err
@@ -129,7 +130,7 @@ subroutine  psb_chalom(x,desc_a,info,jx,ik,work,tran,mode,data)
   endif
   ldx = size(x,1)
   ! check vector correctness
-  call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
@@ -273,8 +274,8 @@ subroutine  psb_chalov(x,desc_a,info,work,tran,mode,data)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me
-  integer(psb_ipk_) :: err_act, ldx, &
-       & m, n, iix, jjx, ix, ijx, nrow, imode, err, liwork,data_
+  integer(psb_ipk_) :: err_act, ldx, iix, jjx, nrow, imode, err, liwork,data_
+  integer(psb_lpk_) :: m, n, ix, ijx
   complex(psb_spk_),pointer :: iwork(:)
   character                :: tran_
   character(len=20)        :: name, ch_err
@@ -319,7 +320,7 @@ subroutine  psb_chalov(x,desc_a,info,work,tran,mode,data)
   endif
   ldx = size(x,1)
   ! check vector correctness
-  call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
diff --git a/base/comm/psb_covrl.f90 b/base/comm/psb_covrl.f90
index 58a5c09d8..7a64ab092 100644
--- a/base/comm/psb_covrl.f90
+++ b/base/comm/psb_covrl.f90
@@ -75,9 +75,9 @@ subroutine  psb_covrl_vect(x,desc_a,info,work,update,mode)
   integer(psb_ipk_), intent(in), optional          :: update,mode
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me, &
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, update_,&
-       & mode_, err, liwork,ldx
+  integer(psb_ipk_) :: ictxt, np, me, err_act, k, iix, jjx, &
+       & nrow, imode, err, liwork,data_, update_, mode_, ncol
+  integer(psb_lpk_) :: m, n, ix, ijx
   complex(psb_spk_),pointer :: iwork(:)
   logical                  :: do_swap
   character(len=20)        :: name, ch_err
@@ -127,7 +127,7 @@ subroutine  psb_covrl_vect(x,desc_a,info,work,update,mode)
   do_swap = (mode_ /= 0)
 
   ! check vector correctness
-  call psb_chkvect(m,ione,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
@@ -200,9 +200,9 @@ subroutine  psb_covrl_multivect(x,desc_a,info,work,update,mode)
   integer(psb_ipk_), intent(in), optional          :: update,mode
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me, &
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, update_,&
-       & mode_, err, liwork,ldx
+  integer(psb_ipk_) :: ictxt, np, me, err_act, k, iix, jjx, &
+       & nrow, imode, err, liwork,data_, update_, mode_, ncol
+  integer(psb_lpk_) :: m, n, ix, ijx
   complex(psb_spk_),pointer :: iwork(:)
   logical                  :: do_swap
   character(len=20)        :: name, ch_err
@@ -252,7 +252,7 @@ subroutine  psb_covrl_multivect(x,desc_a,info,work,update,mode)
   do_swap = (mode_ /= 0)
 
   ! check vector correctness
-  call psb_chkvect(m,ione,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
diff --git a/base/comm/psb_covrl_a.f90 b/base/comm/psb_covrl_a.f90
index 6b6366853..b64a3137c 100644
--- a/base/comm/psb_covrl_a.f90
+++ b/base/comm/psb_covrl_a.f90
@@ -77,8 +77,9 @@ subroutine  psb_covrlm(x,desc_a,info,jx,ik,work,update,mode)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me
-  integer(psb_ipk_) :: err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, maxk, update_,&
+  integer(psb_ipk_) :: err_act, iix, jjx, nrow, ncol, k, maxk, update_,&
        & mode_, err, liwork, ldx
+  integer(psb_lpk_) :: m, n, ix, ijx
   complex(psb_spk_),pointer :: iwork(:), xp(:,:)
   logical                  :: do_swap
   character(len=20)        :: name, ch_err
@@ -137,7 +138,7 @@ subroutine  psb_covrlm(x,desc_a,info,jx,ik,work,update,mode)
   do_swap = (mode_ /= 0)
   ldx = size(x,1)
   ! check vector correctness
-  call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
@@ -271,9 +272,9 @@ subroutine  psb_covrlv(x,desc_a,info,work,update,mode)
   integer(psb_ipk_), intent(in), optional             :: update,mode
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me, &
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, update_,&
-       & mode_, err, liwork, ldx
+  integer(psb_ipk_) :: ictxt, np, me, err_act, iix, jjx, nrow, ncol, &
+       & k, update_, mode_, err, liwork, ldx
+  integer(psb_lpk_) :: m, n, ix, ijx
   complex(psb_spk_),pointer :: iwork(:)
   logical                  :: do_swap
   character(len=20)        :: name, ch_err
@@ -318,7 +319,7 @@ subroutine  psb_covrlv(x,desc_a,info,work,update,mode)
   do_swap = (mode_ /= 0)
   ldx = size(x,1) 
   ! check vector correctness
-  call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
diff --git a/base/comm/psb_cspgather.F90 b/base/comm/psb_cspgather.F90
index 7f6b5a4ee..e9d3780f3 100644
--- a/base/comm/psb_cspgather.F90
+++ b/base/comm/psb_cspgather.F90
@@ -51,6 +51,7 @@ subroutine  psb_csp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
   logical, intent(in), optional   :: keepnum,keeploc
 
   type(psb_c_coo_sparse_mat)      :: loc_coo, glob_coo
+  integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:)
   integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg
   integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl
   logical :: keepnum_, keeploc_
@@ -101,14 +102,19 @@ subroutine  psb_csp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
       call loca%mv_to(loc_coo)
     end if
     nzl = loc_coo%get_nzeros()
-    call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I')
-    call psb_loc_to_glob(loc_coo%ja(1:nzl),desc_a,info,iact='I')
+    call psb_realloc(nzl,locia,info)
+    call psb_realloc(nzl,locja,info)    
+    call psb_loc_to_glob(loc_coo%ia(1:nzl),locia(1:nzl),desc_a,info,iact='I')
+    call psb_loc_to_glob(loc_coo%ja(1:nzl),locja(1:nzl),desc_a,info,iact='I')
     nzbr(:) = 0
     nzbr(me+1) = nzl
     call psb_sum(ictxt,nzbr(1:np))
     nzg = sum(nzbr)
+    if (info == psb_success_) call psb_realloc(nzg,glbia,info)
+    if (info == psb_success_) call psb_realloc(nzg,glbja,info)    
     if (info == psb_success_) call glob_coo%allocate(nrg,ncg,nzg)
     if (info /= psb_success_) goto 9999
+
     do ip=1,np
       idisp(ip) = sum(nzbr(1:ip-1))
     enddo
@@ -117,24 +123,33 @@ subroutine  psb_csp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
          & glob_coo%val,nzbr,idisp,&
          & psb_mpi_c_spk_,icomm,minfo)
     if (minfo == psb_success_) call &
-         & mpi_allgatherv(loc_coo%ia,ndx,psb_mpi_ipk_,&
-         & glob_coo%ia,nzbr,idisp,&
-         & psb_mpi_ipk_,icomm,minfo)
+         & mpi_allgatherv(locia,ndx,psb_mpi_lpk_,&
+         & glbia,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
     if (minfo == psb_success_) call &
-         & mpi_allgatherv(loc_coo%ja,ndx,psb_mpi_ipk_,&
-         & glob_coo%ja,nzbr,idisp,&
-         & psb_mpi_ipk_,icomm,minfo)
+         & mpi_allgatherv(locja,ndx,psb_mpi_lpk_,&
+         & glbja,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
     
     if (minfo /= psb_success_) then 
       info  = minfo
       call psb_errpush(psb_err_internal_error_,name,a_err=' from mpi_allgatherv')
       goto 9999
-    end if
-    
+    end if    
     call loc_coo%free()
+    deallocate(locia,locja, stat=info)
+    !
+    ! Is the code below safe? For very large cases
+    ! the indices in glob_coo will overflow. But then,
+    ! for very large cases it does not make sense to
+    ! gather the matrix on a single procecss anyway...
+    !
+    glob_coo%ia(1:nzg) = glbia(1:nzg)
+    glob_coo%ja(1:nzg) = glbja(1:nzg)
     call glob_coo%set_nzeros(nzg)
     if (present(dupl)) call glob_coo%set_dupl(dupl)
     call globa%mv_from(glob_coo)
+    deallocate(glbia,glbja, stat=info)
 
   else
     write(psb_err_unit,*) 'SP_ALLGATHER: Not implemented yet with keepnum ',keepnum_
diff --git a/base/comm/psb_dgather.f90 b/base/comm/psb_dgather.f90
index 1af953c67..f7ec78da5 100644
--- a/base/comm/psb_dgather.f90
+++ b/base/comm/psb_dgather.f90
@@ -58,8 +58,8 @@ subroutine  psb_dgather_vect(globx, locx, desc_a, info, iroot)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
-  integer(psb_ipk_) :: ierr(5), err_act, n, ilocx, iglobx, jlocx,&
-       & jglobx, lda_locx, lda_globx, m, k, jlx, ilx, i, idx
+  integer(psb_ipk_) :: ierr(5), err_act, jlx, ilx, lda_locx, lda_globx, i
+  integer(psb_lpk_) :: m, n, k, ilocx,  jlocx, idx, iglobx, jglobx
   real(psb_dpk_), allocatable :: llocx(:)
   character(len=20)        :: name, ch_err
 
@@ -171,8 +171,8 @@ subroutine  psb_dgather_multivect(globx, locx, desc_a, info, iroot)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
-  integer(psb_ipk_) :: ierr(5), err_act, n, ilocx, iglobx, jlocx,&
-       & jglobx, lda_locx, lda_globx, m, k, jlx, ilx, i, idx
+  integer(psb_ipk_) :: ierr(5), err_act, jlx, ilx, lda_locx, lda_globx, i 
+  integer(psb_lpk_) :: m, n, k, ilocx, jlocx, idx, iglobx, jglobx
   real(psb_dpk_), allocatable :: llocx(:,:)
   character(len=20)        :: name, ch_err
 
diff --git a/base/comm/psb_dgather_a.f90 b/base/comm/psb_dgather_a.f90
index ed4af9f98..5dd202952 100644
--- a/base/comm/psb_dgather_a.f90
+++ b/base/comm/psb_dgather_a.f90
@@ -58,9 +58,9 @@ subroutine  psb_dgatherm(globx, locx, desc_a, info, iroot)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
-  integer(psb_ipk_) :: ierr(5), err_act, n, ilocx, iglobx, jlocx,&
-       & jglobx, lda_locx, lda_globx, m, lock, globk, maxk, k, jlx, &
-       & ilx, i, j, idx
+  integer(psb_ipk_) :: ierr(5), err_act, lda_locx, lda_globx, lock, globk,&
+       & maxk, k, jlx, ilx, i, j
+  integer(psb_lpk_) :: m, n, ilocx,  jlocx, idx, iglobx, jglobx
 
   character(len=20)        :: name, ch_err
 
@@ -230,8 +230,9 @@ subroutine  psb_dgatherv(globx, locx, desc_a, info, iroot)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
-  integer(psb_ipk_) :: ierr(5), err_act, n, ilocx, iglobx, jlocx,&
-       & jglobx, lda_locx, lda_globx, m, k, jlx, ilx, i, idx
+  integer(psb_ipk_) :: ierr(5), err_act, lda_locx, lda_globx, lock, globk,&
+       & maxk, k, jlx, ilx, i, j
+  integer(psb_lpk_) :: m, n, ilocx,  jlocx, idx, iglobx, jglobx
 
   character(len=20)        :: name, ch_err
 
diff --git a/base/comm/psb_dhalo.f90 b/base/comm/psb_dhalo.f90
index 058305822..4f48007b0 100644
--- a/base/comm/psb_dhalo.f90
+++ b/base/comm/psb_dhalo.f90
@@ -66,9 +66,9 @@ subroutine  psb_dhalo_vect(x,desc_a,info,work,tran,mode,data)
   character, intent(in), optional         :: tran
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me,&
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, imode,&
-       & err, liwork,data_
+  integer(psb_ipk_) :: ictxt, np, me, err_act, iix, jjx, &
+       & nrow, imode, err, liwork,data_
+  integer(psb_lpk_) :: m, n, ix, ijx
   real(psb_dpk_),pointer :: iwork(:)
   character                 :: tran_
   character(len=20)         :: name, ch_err
@@ -119,7 +119,7 @@ subroutine  psb_dhalo_vect(x,desc_a,info,work,tran,mode,data)
   endif
 
   ! check vector correctness
-  call psb_chkvect(m,ione,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
@@ -205,9 +205,9 @@ subroutine  psb_dhalo_multivect(x,desc_a,info,work,tran,mode,data)
   character, intent(in), optional         :: tran
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me,&
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, imode,&
-       & err, liwork,data_
+  integer(psb_ipk_) :: ictxt, np, me, err_act, iix, jjx, &
+       & nrow, imode, err, liwork,data_
+  integer(psb_lpk_) :: m, n, ix, ijx
   real(psb_dpk_),pointer :: iwork(:)
   character                 :: tran_
   character(len=20)         :: name, ch_err
@@ -258,7 +258,7 @@ subroutine  psb_dhalo_multivect(x,desc_a,info,work,tran,mode,data)
   endif
 
   ! check vector correctness
-  call psb_chkvect(m,ione,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
diff --git a/base/comm/psb_dhalo_a.f90 b/base/comm/psb_dhalo_a.f90
index 497667f1a..89bde9aea 100644
--- a/base/comm/psb_dhalo_a.f90
+++ b/base/comm/psb_dhalo_a.f90
@@ -66,8 +66,9 @@ subroutine  psb_dhalom(x,desc_a,info,jx,ik,work,tran,mode,data)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me
-  integer(psb_ipk_) :: err_act, m, n, iix, jjx, ix, ijx, k, maxk, nrow, imode, i,&
+  integer(psb_ipk_) :: err_act, iix, jjx, k, maxk, nrow, imode, i,&
        & err, liwork,data_, ldx
+  integer(psb_lpk_) :: m, n, ix, ijx
   real(psb_dpk_),pointer :: iwork(:), xp(:,:)
   character                :: tran_
   character(len=20)        :: name, ch_err
@@ -129,7 +130,7 @@ subroutine  psb_dhalom(x,desc_a,info,jx,ik,work,tran,mode,data)
   endif
   ldx = size(x,1)
   ! check vector correctness
-  call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
@@ -273,8 +274,8 @@ subroutine  psb_dhalov(x,desc_a,info,work,tran,mode,data)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me
-  integer(psb_ipk_) :: err_act, ldx, &
-       & m, n, iix, jjx, ix, ijx, nrow, imode, err, liwork,data_
+  integer(psb_ipk_) :: err_act, ldx, iix, jjx, nrow, imode, err, liwork,data_
+  integer(psb_lpk_) :: m, n, ix, ijx
   real(psb_dpk_),pointer :: iwork(:)
   character                :: tran_
   character(len=20)        :: name, ch_err
@@ -319,7 +320,7 @@ subroutine  psb_dhalov(x,desc_a,info,work,tran,mode,data)
   endif
   ldx = size(x,1)
   ! check vector correctness
-  call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
diff --git a/base/comm/psb_dovrl.f90 b/base/comm/psb_dovrl.f90
index 13ac175d8..b255511a8 100644
--- a/base/comm/psb_dovrl.f90
+++ b/base/comm/psb_dovrl.f90
@@ -75,9 +75,9 @@ subroutine  psb_dovrl_vect(x,desc_a,info,work,update,mode)
   integer(psb_ipk_), intent(in), optional          :: update,mode
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me, &
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, update_,&
-       & mode_, err, liwork,ldx
+  integer(psb_ipk_) :: ictxt, np, me, err_act, k, iix, jjx, &
+       & nrow, imode, err, liwork,data_, update_, mode_, ncol
+  integer(psb_lpk_) :: m, n, ix, ijx
   real(psb_dpk_),pointer :: iwork(:)
   logical                  :: do_swap
   character(len=20)        :: name, ch_err
@@ -127,7 +127,7 @@ subroutine  psb_dovrl_vect(x,desc_a,info,work,update,mode)
   do_swap = (mode_ /= 0)
 
   ! check vector correctness
-  call psb_chkvect(m,ione,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
@@ -200,9 +200,9 @@ subroutine  psb_dovrl_multivect(x,desc_a,info,work,update,mode)
   integer(psb_ipk_), intent(in), optional          :: update,mode
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me, &
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, update_,&
-       & mode_, err, liwork,ldx
+  integer(psb_ipk_) :: ictxt, np, me, err_act, k, iix, jjx, &
+       & nrow, imode, err, liwork,data_, update_, mode_, ncol
+  integer(psb_lpk_) :: m, n, ix, ijx
   real(psb_dpk_),pointer :: iwork(:)
   logical                  :: do_swap
   character(len=20)        :: name, ch_err
@@ -252,7 +252,7 @@ subroutine  psb_dovrl_multivect(x,desc_a,info,work,update,mode)
   do_swap = (mode_ /= 0)
 
   ! check vector correctness
-  call psb_chkvect(m,ione,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
diff --git a/base/comm/psb_dovrl_a.f90 b/base/comm/psb_dovrl_a.f90
index e38dfaa07..12732460d 100644
--- a/base/comm/psb_dovrl_a.f90
+++ b/base/comm/psb_dovrl_a.f90
@@ -77,8 +77,9 @@ subroutine  psb_dovrlm(x,desc_a,info,jx,ik,work,update,mode)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me
-  integer(psb_ipk_) :: err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, maxk, update_,&
+  integer(psb_ipk_) :: err_act, iix, jjx, nrow, ncol, k, maxk, update_,&
        & mode_, err, liwork, ldx
+  integer(psb_lpk_) :: m, n, ix, ijx
   real(psb_dpk_),pointer :: iwork(:), xp(:,:)
   logical                  :: do_swap
   character(len=20)        :: name, ch_err
@@ -137,7 +138,7 @@ subroutine  psb_dovrlm(x,desc_a,info,jx,ik,work,update,mode)
   do_swap = (mode_ /= 0)
   ldx = size(x,1)
   ! check vector correctness
-  call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
@@ -271,9 +272,9 @@ subroutine  psb_dovrlv(x,desc_a,info,work,update,mode)
   integer(psb_ipk_), intent(in), optional             :: update,mode
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me, &
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, update_,&
-       & mode_, err, liwork, ldx
+  integer(psb_ipk_) :: ictxt, np, me, err_act, iix, jjx, nrow, ncol, &
+       & k, update_, mode_, err, liwork, ldx
+  integer(psb_lpk_) :: m, n, ix, ijx
   real(psb_dpk_),pointer :: iwork(:)
   logical                  :: do_swap
   character(len=20)        :: name, ch_err
@@ -318,7 +319,7 @@ subroutine  psb_dovrlv(x,desc_a,info,work,update,mode)
   do_swap = (mode_ /= 0)
   ldx = size(x,1) 
   ! check vector correctness
-  call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
diff --git a/base/comm/psb_dspgather.F90 b/base/comm/psb_dspgather.F90
index afebc345a..61416cb1e 100644
--- a/base/comm/psb_dspgather.F90
+++ b/base/comm/psb_dspgather.F90
@@ -51,6 +51,7 @@ subroutine  psb_dsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
   logical, intent(in), optional   :: keepnum,keeploc
 
   type(psb_d_coo_sparse_mat)      :: loc_coo, glob_coo
+  integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:)
   integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg
   integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl
   logical :: keepnum_, keeploc_
@@ -101,14 +102,19 @@ subroutine  psb_dsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
       call loca%mv_to(loc_coo)
     end if
     nzl = loc_coo%get_nzeros()
-    call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I')
-    call psb_loc_to_glob(loc_coo%ja(1:nzl),desc_a,info,iact='I')
+    call psb_realloc(nzl,locia,info)
+    call psb_realloc(nzl,locja,info)    
+    call psb_loc_to_glob(loc_coo%ia(1:nzl),locia(1:nzl),desc_a,info,iact='I')
+    call psb_loc_to_glob(loc_coo%ja(1:nzl),locja(1:nzl),desc_a,info,iact='I')
     nzbr(:) = 0
     nzbr(me+1) = nzl
     call psb_sum(ictxt,nzbr(1:np))
     nzg = sum(nzbr)
+    if (info == psb_success_) call psb_realloc(nzg,glbia,info)
+    if (info == psb_success_) call psb_realloc(nzg,glbja,info)    
     if (info == psb_success_) call glob_coo%allocate(nrg,ncg,nzg)
     if (info /= psb_success_) goto 9999
+
     do ip=1,np
       idisp(ip) = sum(nzbr(1:ip-1))
     enddo
@@ -117,24 +123,33 @@ subroutine  psb_dsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
          & glob_coo%val,nzbr,idisp,&
          & psb_mpi_r_dpk_,icomm,minfo)
     if (minfo == psb_success_) call &
-         & mpi_allgatherv(loc_coo%ia,ndx,psb_mpi_ipk_,&
-         & glob_coo%ia,nzbr,idisp,&
-         & psb_mpi_ipk_,icomm,minfo)
+         & mpi_allgatherv(locia,ndx,psb_mpi_lpk_,&
+         & glbia,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
     if (minfo == psb_success_) call &
-         & mpi_allgatherv(loc_coo%ja,ndx,psb_mpi_ipk_,&
-         & glob_coo%ja,nzbr,idisp,&
-         & psb_mpi_ipk_,icomm,minfo)
+         & mpi_allgatherv(locja,ndx,psb_mpi_lpk_,&
+         & glbja,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
     
     if (minfo /= psb_success_) then 
       info  = minfo
       call psb_errpush(psb_err_internal_error_,name,a_err=' from mpi_allgatherv')
       goto 9999
-    end if
-    
+    end if    
     call loc_coo%free()
+    deallocate(locia,locja, stat=info)
+    !
+    ! Is the code below safe? For very large cases
+    ! the indices in glob_coo will overflow. But then,
+    ! for very large cases it does not make sense to
+    ! gather the matrix on a single procecss anyway...
+    !
+    glob_coo%ia(1:nzg) = glbia(1:nzg)
+    glob_coo%ja(1:nzg) = glbja(1:nzg)
     call glob_coo%set_nzeros(nzg)
     if (present(dupl)) call glob_coo%set_dupl(dupl)
     call globa%mv_from(glob_coo)
+    deallocate(glbia,glbja, stat=info)
 
   else
     write(psb_err_unit,*) 'SP_ALLGATHER: Not implemented yet with keepnum ',keepnum_
diff --git a/base/comm/psb_egather_a.f90 b/base/comm/psb_egather_a.f90
index c10e55658..77d76dd27 100644
--- a/base/comm/psb_egather_a.f90
+++ b/base/comm/psb_egather_a.f90
@@ -58,9 +58,9 @@ subroutine  psb_egatherm(globx, locx, desc_a, info, iroot)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
-  integer(psb_ipk_) :: ierr(5), err_act, n, ilocx, iglobx, jlocx,&
-       & jglobx, lda_locx, lda_globx, m, lock, globk, maxk, k, jlx, &
-       & ilx, i, j, idx
+  integer(psb_ipk_) :: ierr(5), err_act, lda_locx, lda_globx, lock, globk,&
+       & maxk, k, jlx, ilx, i, j
+  integer(psb_lpk_) :: m, n, ilocx,  jlocx, idx, iglobx, jglobx
 
   character(len=20)        :: name, ch_err
 
@@ -230,8 +230,9 @@ subroutine  psb_egatherv(globx, locx, desc_a, info, iroot)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
-  integer(psb_ipk_) :: ierr(5), err_act, n, ilocx, iglobx, jlocx,&
-       & jglobx, lda_locx, lda_globx, m, k, jlx, ilx, i, idx
+  integer(psb_ipk_) :: ierr(5), err_act, lda_locx, lda_globx, lock, globk,&
+       & maxk, k, jlx, ilx, i, j
+  integer(psb_lpk_) :: m, n, ilocx,  jlocx, idx, iglobx, jglobx
 
   character(len=20)        :: name, ch_err
 
diff --git a/base/comm/psb_ehalo_a.f90 b/base/comm/psb_ehalo_a.f90
index 171495456..b9921c0db 100644
--- a/base/comm/psb_ehalo_a.f90
+++ b/base/comm/psb_ehalo_a.f90
@@ -66,8 +66,9 @@ subroutine  psb_ehalom(x,desc_a,info,jx,ik,work,tran,mode,data)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me
-  integer(psb_ipk_) :: err_act, m, n, iix, jjx, ix, ijx, k, maxk, nrow, imode, i,&
+  integer(psb_ipk_) :: err_act, iix, jjx, k, maxk, nrow, imode, i,&
        & err, liwork,data_, ldx
+  integer(psb_lpk_) :: m, n, ix, ijx
   integer(psb_epk_),pointer :: iwork(:), xp(:,:)
   character                :: tran_
   character(len=20)        :: name, ch_err
@@ -129,7 +130,7 @@ subroutine  psb_ehalom(x,desc_a,info,jx,ik,work,tran,mode,data)
   endif
   ldx = size(x,1)
   ! check vector correctness
-  call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
@@ -273,8 +274,8 @@ subroutine  psb_ehalov(x,desc_a,info,work,tran,mode,data)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me
-  integer(psb_ipk_) :: err_act, ldx, &
-       & m, n, iix, jjx, ix, ijx, nrow, imode, err, liwork,data_
+  integer(psb_ipk_) :: err_act, ldx, iix, jjx, nrow, imode, err, liwork,data_
+  integer(psb_lpk_) :: m, n, ix, ijx
   integer(psb_epk_),pointer :: iwork(:)
   character                :: tran_
   character(len=20)        :: name, ch_err
@@ -319,7 +320,7 @@ subroutine  psb_ehalov(x,desc_a,info,work,tran,mode,data)
   endif
   ldx = size(x,1)
   ! check vector correctness
-  call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
diff --git a/base/comm/psb_eovrl_a.f90 b/base/comm/psb_eovrl_a.f90
index 9da7708f3..3c5eccc4d 100644
--- a/base/comm/psb_eovrl_a.f90
+++ b/base/comm/psb_eovrl_a.f90
@@ -77,8 +77,9 @@ subroutine  psb_eovrlm(x,desc_a,info,jx,ik,work,update,mode)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me
-  integer(psb_ipk_) :: err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, maxk, update_,&
+  integer(psb_ipk_) :: err_act, iix, jjx, nrow, ncol, k, maxk, update_,&
        & mode_, err, liwork, ldx
+  integer(psb_lpk_) :: m, n, ix, ijx
   integer(psb_epk_),pointer :: iwork(:), xp(:,:)
   logical                  :: do_swap
   character(len=20)        :: name, ch_err
@@ -137,7 +138,7 @@ subroutine  psb_eovrlm(x,desc_a,info,jx,ik,work,update,mode)
   do_swap = (mode_ /= 0)
   ldx = size(x,1)
   ! check vector correctness
-  call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
@@ -271,9 +272,9 @@ subroutine  psb_eovrlv(x,desc_a,info,work,update,mode)
   integer(psb_ipk_), intent(in), optional             :: update,mode
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me, &
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, update_,&
-       & mode_, err, liwork, ldx
+  integer(psb_ipk_) :: ictxt, np, me, err_act, iix, jjx, nrow, ncol, &
+       & k, update_, mode_, err, liwork, ldx
+  integer(psb_lpk_) :: m, n, ix, ijx
   integer(psb_epk_),pointer :: iwork(:)
   logical                  :: do_swap
   character(len=20)        :: name, ch_err
@@ -318,7 +319,7 @@ subroutine  psb_eovrlv(x,desc_a,info,work,update,mode)
   do_swap = (mode_ /= 0)
   ldx = size(x,1) 
   ! check vector correctness
-  call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
diff --git a/base/comm/psb_igather.f90 b/base/comm/psb_igather.f90
index dcef40dac..8b13af798 100644
--- a/base/comm/psb_igather.f90
+++ b/base/comm/psb_igather.f90
@@ -58,8 +58,8 @@ subroutine  psb_igather_vect(globx, locx, desc_a, info, iroot)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
-  integer(psb_ipk_) :: ierr(5), err_act, n, ilocx, iglobx, jlocx,&
-       & jglobx, lda_locx, lda_globx, m, k, jlx, ilx, i, idx
+  integer(psb_ipk_) :: ierr(5), err_act, jlx, ilx, lda_locx, lda_globx, i
+  integer(psb_lpk_) :: m, n, k, ilocx,  jlocx, idx, iglobx, jglobx
   integer(psb_ipk_), allocatable :: llocx(:)
   character(len=20)        :: name, ch_err
 
@@ -171,8 +171,8 @@ subroutine  psb_igather_multivect(globx, locx, desc_a, info, iroot)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
-  integer(psb_ipk_) :: ierr(5), err_act, n, ilocx, iglobx, jlocx,&
-       & jglobx, lda_locx, lda_globx, m, k, jlx, ilx, i, idx
+  integer(psb_ipk_) :: ierr(5), err_act, jlx, ilx, lda_locx, lda_globx, i 
+  integer(psb_lpk_) :: m, n, k, ilocx, jlocx, idx, iglobx, jglobx
   integer(psb_ipk_), allocatable :: llocx(:,:)
   character(len=20)        :: name, ch_err
 
diff --git a/base/comm/psb_ihalo.f90 b/base/comm/psb_ihalo.f90
index cbcd15b04..ea401ae39 100644
--- a/base/comm/psb_ihalo.f90
+++ b/base/comm/psb_ihalo.f90
@@ -66,9 +66,9 @@ subroutine  psb_ihalo_vect(x,desc_a,info,work,tran,mode,data)
   character, intent(in), optional         :: tran
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me,&
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, imode,&
-       & err, liwork,data_
+  integer(psb_ipk_) :: ictxt, np, me, err_act, iix, jjx, &
+       & nrow, imode, err, liwork,data_
+  integer(psb_lpk_) :: m, n, ix, ijx
   integer(psb_ipk_),pointer :: iwork(:)
   character                 :: tran_
   character(len=20)         :: name, ch_err
@@ -119,7 +119,7 @@ subroutine  psb_ihalo_vect(x,desc_a,info,work,tran,mode,data)
   endif
 
   ! check vector correctness
-  call psb_chkvect(m,ione,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
@@ -205,9 +205,9 @@ subroutine  psb_ihalo_multivect(x,desc_a,info,work,tran,mode,data)
   character, intent(in), optional         :: tran
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me,&
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, imode,&
-       & err, liwork,data_
+  integer(psb_ipk_) :: ictxt, np, me, err_act, iix, jjx, &
+       & nrow, imode, err, liwork,data_
+  integer(psb_lpk_) :: m, n, ix, ijx
   integer(psb_ipk_),pointer :: iwork(:)
   character                 :: tran_
   character(len=20)         :: name, ch_err
@@ -258,7 +258,7 @@ subroutine  psb_ihalo_multivect(x,desc_a,info,work,tran,mode,data)
   endif
 
   ! check vector correctness
-  call psb_chkvect(m,ione,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
diff --git a/base/comm/psb_iovrl.f90 b/base/comm/psb_iovrl.f90
index 6cc917044..dcc97b874 100644
--- a/base/comm/psb_iovrl.f90
+++ b/base/comm/psb_iovrl.f90
@@ -75,9 +75,9 @@ subroutine  psb_iovrl_vect(x,desc_a,info,work,update,mode)
   integer(psb_ipk_), intent(in), optional          :: update,mode
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me, &
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, update_,&
-       & mode_, err, liwork,ldx
+  integer(psb_ipk_) :: ictxt, np, me, err_act, k, iix, jjx, &
+       & nrow, imode, err, liwork,data_, update_, mode_, ncol
+  integer(psb_lpk_) :: m, n, ix, ijx
   integer(psb_ipk_),pointer :: iwork(:)
   logical                  :: do_swap
   character(len=20)        :: name, ch_err
@@ -127,7 +127,7 @@ subroutine  psb_iovrl_vect(x,desc_a,info,work,update,mode)
   do_swap = (mode_ /= 0)
 
   ! check vector correctness
-  call psb_chkvect(m,ione,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
@@ -200,9 +200,9 @@ subroutine  psb_iovrl_multivect(x,desc_a,info,work,update,mode)
   integer(psb_ipk_), intent(in), optional          :: update,mode
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me, &
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, update_,&
-       & mode_, err, liwork,ldx
+  integer(psb_ipk_) :: ictxt, np, me, err_act, k, iix, jjx, &
+       & nrow, imode, err, liwork,data_, update_, mode_, ncol
+  integer(psb_lpk_) :: m, n, ix, ijx
   integer(psb_ipk_),pointer :: iwork(:)
   logical                  :: do_swap
   character(len=20)        :: name, ch_err
@@ -252,7 +252,7 @@ subroutine  psb_iovrl_multivect(x,desc_a,info,work,update,mode)
   do_swap = (mode_ /= 0)
 
   ! check vector correctness
-  call psb_chkvect(m,ione,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
diff --git a/base/comm/psb_ispgather.F90 b/base/comm/psb_ispgather.F90
index 2ef3a5f50..08d4eea54 100644
--- a/base/comm/psb_ispgather.F90
+++ b/base/comm/psb_ispgather.F90
@@ -51,6 +51,7 @@ subroutine  psb_isp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
   logical, intent(in), optional   :: keepnum,keeploc
 
   type(psb_i_coo_sparse_mat)      :: loc_coo, glob_coo
+  integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:)
   integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg
   integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl
   logical :: keepnum_, keeploc_
@@ -101,14 +102,19 @@ subroutine  psb_isp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
       call loca%mv_to(loc_coo)
     end if
     nzl = loc_coo%get_nzeros()
-    call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I')
-    call psb_loc_to_glob(loc_coo%ja(1:nzl),desc_a,info,iact='I')
+    call psb_realloc(nzl,locia,info)
+    call psb_realloc(nzl,locja,info)    
+    call psb_loc_to_glob(loc_coo%ia(1:nzl),locia(1:nzl),desc_a,info,iact='I')
+    call psb_loc_to_glob(loc_coo%ja(1:nzl),locja(1:nzl),desc_a,info,iact='I')
     nzbr(:) = 0
     nzbr(me+1) = nzl
     call psb_sum(ictxt,nzbr(1:np))
     nzg = sum(nzbr)
+    if (info == psb_success_) call psb_realloc(nzg,glbia,info)
+    if (info == psb_success_) call psb_realloc(nzg,glbja,info)    
     if (info == psb_success_) call glob_coo%allocate(nrg,ncg,nzg)
     if (info /= psb_success_) goto 9999
+
     do ip=1,np
       idisp(ip) = sum(nzbr(1:ip-1))
     enddo
@@ -117,24 +123,33 @@ subroutine  psb_isp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
          & glob_coo%val,nzbr,idisp,&
          & psb_mpi_ipk_,icomm,minfo)
     if (minfo == psb_success_) call &
-         & mpi_allgatherv(loc_coo%ia,ndx,psb_mpi_ipk_,&
-         & glob_coo%ia,nzbr,idisp,&
-         & psb_mpi_ipk_,icomm,minfo)
+         & mpi_allgatherv(locia,ndx,psb_mpi_lpk_,&
+         & glbia,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
     if (minfo == psb_success_) call &
-         & mpi_allgatherv(loc_coo%ja,ndx,psb_mpi_ipk_,&
-         & glob_coo%ja,nzbr,idisp,&
-         & psb_mpi_ipk_,icomm,minfo)
+         & mpi_allgatherv(locja,ndx,psb_mpi_lpk_,&
+         & glbja,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
     
     if (minfo /= psb_success_) then 
       info  = minfo
       call psb_errpush(psb_err_internal_error_,name,a_err=' from mpi_allgatherv')
       goto 9999
-    end if
-    
+    end if    
     call loc_coo%free()
+    deallocate(locia,locja, stat=info)
+    !
+    ! Is the code below safe? For very large cases
+    ! the indices in glob_coo will overflow. But then,
+    ! for very large cases it does not make sense to
+    ! gather the matrix on a single procecss anyway...
+    !
+    glob_coo%ia(1:nzg) = glbia(1:nzg)
+    glob_coo%ja(1:nzg) = glbja(1:nzg)
     call glob_coo%set_nzeros(nzg)
     if (present(dupl)) call glob_coo%set_dupl(dupl)
     call globa%mv_from(glob_coo)
+    deallocate(glbia,glbja, stat=info)
 
   else
     write(psb_err_unit,*) 'SP_ALLGATHER: Not implemented yet with keepnum ',keepnum_
diff --git a/base/comm/psb_lgather.f90 b/base/comm/psb_lgather.f90
index 161a0a334..9b6f0fcad 100644
--- a/base/comm/psb_lgather.f90
+++ b/base/comm/psb_lgather.f90
@@ -58,8 +58,8 @@ subroutine  psb_lgather_vect(globx, locx, desc_a, info, iroot)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
-  integer(psb_ipk_) :: ierr(5), err_act, n, ilocx, iglobx, jlocx,&
-       & jglobx, lda_locx, lda_globx, m, k, jlx, ilx, i, idx
+  integer(psb_ipk_) :: ierr(5), err_act, jlx, ilx, lda_locx, lda_globx, i
+  integer(psb_lpk_) :: m, n, k, ilocx,  jlocx, idx, iglobx, jglobx
   integer(psb_lpk_), allocatable :: llocx(:)
   character(len=20)        :: name, ch_err
 
@@ -171,8 +171,8 @@ subroutine  psb_lgather_multivect(globx, locx, desc_a, info, iroot)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
-  integer(psb_ipk_) :: ierr(5), err_act, n, ilocx, iglobx, jlocx,&
-       & jglobx, lda_locx, lda_globx, m, k, jlx, ilx, i, idx
+  integer(psb_ipk_) :: ierr(5), err_act, jlx, ilx, lda_locx, lda_globx, i 
+  integer(psb_lpk_) :: m, n, k, ilocx, jlocx, idx, iglobx, jglobx
   integer(psb_lpk_), allocatable :: llocx(:,:)
   character(len=20)        :: name, ch_err
 
diff --git a/base/comm/psb_lhalo.f90 b/base/comm/psb_lhalo.f90
index e9a46f434..f08bdf704 100644
--- a/base/comm/psb_lhalo.f90
+++ b/base/comm/psb_lhalo.f90
@@ -66,9 +66,9 @@ subroutine  psb_lhalo_vect(x,desc_a,info,work,tran,mode,data)
   character, intent(in), optional         :: tran
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me,&
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, imode,&
-       & err, liwork,data_
+  integer(psb_ipk_) :: ictxt, np, me, err_act, iix, jjx, &
+       & nrow, imode, err, liwork,data_
+  integer(psb_lpk_) :: m, n, ix, ijx
   integer(psb_lpk_),pointer :: iwork(:)
   character                 :: tran_
   character(len=20)         :: name, ch_err
@@ -119,7 +119,7 @@ subroutine  psb_lhalo_vect(x,desc_a,info,work,tran,mode,data)
   endif
 
   ! check vector correctness
-  call psb_chkvect(m,ione,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
@@ -205,9 +205,9 @@ subroutine  psb_lhalo_multivect(x,desc_a,info,work,tran,mode,data)
   character, intent(in), optional         :: tran
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me,&
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, imode,&
-       & err, liwork,data_
+  integer(psb_ipk_) :: ictxt, np, me, err_act, iix, jjx, &
+       & nrow, imode, err, liwork,data_
+  integer(psb_lpk_) :: m, n, ix, ijx
   integer(psb_lpk_),pointer :: iwork(:)
   character                 :: tran_
   character(len=20)         :: name, ch_err
@@ -258,7 +258,7 @@ subroutine  psb_lhalo_multivect(x,desc_a,info,work,tran,mode,data)
   endif
 
   ! check vector correctness
-  call psb_chkvect(m,ione,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
diff --git a/base/comm/psb_lovrl.f90 b/base/comm/psb_lovrl.f90
index edb4049a6..94adc24d4 100644
--- a/base/comm/psb_lovrl.f90
+++ b/base/comm/psb_lovrl.f90
@@ -75,9 +75,9 @@ subroutine  psb_lovrl_vect(x,desc_a,info,work,update,mode)
   integer(psb_ipk_), intent(in), optional          :: update,mode
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me, &
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, update_,&
-       & mode_, err, liwork,ldx
+  integer(psb_ipk_) :: ictxt, np, me, err_act, k, iix, jjx, &
+       & nrow, imode, err, liwork,data_, update_, mode_, ncol
+  integer(psb_lpk_) :: m, n, ix, ijx
   integer(psb_lpk_),pointer :: iwork(:)
   logical                  :: do_swap
   character(len=20)        :: name, ch_err
@@ -127,7 +127,7 @@ subroutine  psb_lovrl_vect(x,desc_a,info,work,update,mode)
   do_swap = (mode_ /= 0)
 
   ! check vector correctness
-  call psb_chkvect(m,ione,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
@@ -200,9 +200,9 @@ subroutine  psb_lovrl_multivect(x,desc_a,info,work,update,mode)
   integer(psb_ipk_), intent(in), optional          :: update,mode
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me, &
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, update_,&
-       & mode_, err, liwork,ldx
+  integer(psb_ipk_) :: ictxt, np, me, err_act, k, iix, jjx, &
+       & nrow, imode, err, liwork,data_, update_, mode_, ncol
+  integer(psb_lpk_) :: m, n, ix, ijx
   integer(psb_lpk_),pointer :: iwork(:)
   logical                  :: do_swap
   character(len=20)        :: name, ch_err
@@ -252,7 +252,7 @@ subroutine  psb_lovrl_multivect(x,desc_a,info,work,update,mode)
   do_swap = (mode_ /= 0)
 
   ! check vector correctness
-  call psb_chkvect(m,ione,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
diff --git a/base/comm/psb_lspgather.F90 b/base/comm/psb_lspgather.F90
index 344d9c982..e3c1eb051 100644
--- a/base/comm/psb_lspgather.F90
+++ b/base/comm/psb_lspgather.F90
@@ -51,6 +51,7 @@ subroutine  psb_lsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
   logical, intent(in), optional   :: keepnum,keeploc
 
   type(psb_l_coo_sparse_mat)      :: loc_coo, glob_coo
+  integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:)
   integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg
   integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl
   logical :: keepnum_, keeploc_
@@ -101,14 +102,19 @@ subroutine  psb_lsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
       call loca%mv_to(loc_coo)
     end if
     nzl = loc_coo%get_nzeros()
-    call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I')
-    call psb_loc_to_glob(loc_coo%ja(1:nzl),desc_a,info,iact='I')
+    call psb_realloc(nzl,locia,info)
+    call psb_realloc(nzl,locja,info)    
+    call psb_loc_to_glob(loc_coo%ia(1:nzl),locia(1:nzl),desc_a,info,iact='I')
+    call psb_loc_to_glob(loc_coo%ja(1:nzl),locja(1:nzl),desc_a,info,iact='I')
     nzbr(:) = 0
     nzbr(me+1) = nzl
     call psb_sum(ictxt,nzbr(1:np))
     nzg = sum(nzbr)
+    if (info == psb_success_) call psb_realloc(nzg,glbia,info)
+    if (info == psb_success_) call psb_realloc(nzg,glbja,info)    
     if (info == psb_success_) call glob_coo%allocate(nrg,ncg,nzg)
     if (info /= psb_success_) goto 9999
+
     do ip=1,np
       idisp(ip) = sum(nzbr(1:ip-1))
     enddo
@@ -117,24 +123,33 @@ subroutine  psb_lsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
          & glob_coo%val,nzbr,idisp,&
          & psb_mpi_lpk_,icomm,minfo)
     if (minfo == psb_success_) call &
-         & mpi_allgatherv(loc_coo%ia,ndx,psb_mpi_ipk_,&
-         & glob_coo%ia,nzbr,idisp,&
-         & psb_mpi_ipk_,icomm,minfo)
+         & mpi_allgatherv(locia,ndx,psb_mpi_lpk_,&
+         & glbia,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
     if (minfo == psb_success_) call &
-         & mpi_allgatherv(loc_coo%ja,ndx,psb_mpi_ipk_,&
-         & glob_coo%ja,nzbr,idisp,&
-         & psb_mpi_ipk_,icomm,minfo)
+         & mpi_allgatherv(locja,ndx,psb_mpi_lpk_,&
+         & glbja,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
     
     if (minfo /= psb_success_) then 
       info  = minfo
       call psb_errpush(psb_err_internal_error_,name,a_err=' from mpi_allgatherv')
       goto 9999
-    end if
-    
+    end if    
     call loc_coo%free()
+    deallocate(locia,locja, stat=info)
+    !
+    ! Is the code below safe? For very large cases
+    ! the indices in glob_coo will overflow. But then,
+    ! for very large cases it does not make sense to
+    ! gather the matrix on a single procecss anyway...
+    !
+    glob_coo%ia(1:nzg) = glbia(1:nzg)
+    glob_coo%ja(1:nzg) = glbja(1:nzg)
     call glob_coo%set_nzeros(nzg)
     if (present(dupl)) call glob_coo%set_dupl(dupl)
     call globa%mv_from(glob_coo)
+    deallocate(glbia,glbja, stat=info)
 
   else
     write(psb_err_unit,*) 'SP_ALLGATHER: Not implemented yet with keepnum ',keepnum_
diff --git a/base/comm/psb_mgather_a.f90 b/base/comm/psb_mgather_a.f90
index 2c760794c..6b9467ff6 100644
--- a/base/comm/psb_mgather_a.f90
+++ b/base/comm/psb_mgather_a.f90
@@ -58,9 +58,9 @@ subroutine  psb_mgatherm(globx, locx, desc_a, info, iroot)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
-  integer(psb_ipk_) :: ierr(5), err_act, n, ilocx, iglobx, jlocx,&
-       & jglobx, lda_locx, lda_globx, m, lock, globk, maxk, k, jlx, &
-       & ilx, i, j, idx
+  integer(psb_ipk_) :: ierr(5), err_act, lda_locx, lda_globx, lock, globk,&
+       & maxk, k, jlx, ilx, i, j
+  integer(psb_lpk_) :: m, n, ilocx,  jlocx, idx, iglobx, jglobx
 
   character(len=20)        :: name, ch_err
 
@@ -230,8 +230,9 @@ subroutine  psb_mgatherv(globx, locx, desc_a, info, iroot)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
-  integer(psb_ipk_) :: ierr(5), err_act, n, ilocx, iglobx, jlocx,&
-       & jglobx, lda_locx, lda_globx, m, k, jlx, ilx, i, idx
+  integer(psb_ipk_) :: ierr(5), err_act, lda_locx, lda_globx, lock, globk,&
+       & maxk, k, jlx, ilx, i, j
+  integer(psb_lpk_) :: m, n, ilocx,  jlocx, idx, iglobx, jglobx
 
   character(len=20)        :: name, ch_err
 
diff --git a/base/comm/psb_mhalo_a.f90 b/base/comm/psb_mhalo_a.f90
index 67d1c3cc1..6762633af 100644
--- a/base/comm/psb_mhalo_a.f90
+++ b/base/comm/psb_mhalo_a.f90
@@ -66,8 +66,9 @@ subroutine  psb_mhalom(x,desc_a,info,jx,ik,work,tran,mode,data)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me
-  integer(psb_ipk_) :: err_act, m, n, iix, jjx, ix, ijx, k, maxk, nrow, imode, i,&
+  integer(psb_ipk_) :: err_act, iix, jjx, k, maxk, nrow, imode, i,&
        & err, liwork,data_, ldx
+  integer(psb_lpk_) :: m, n, ix, ijx
   integer(psb_mpk_),pointer :: iwork(:), xp(:,:)
   character                :: tran_
   character(len=20)        :: name, ch_err
@@ -129,7 +130,7 @@ subroutine  psb_mhalom(x,desc_a,info,jx,ik,work,tran,mode,data)
   endif
   ldx = size(x,1)
   ! check vector correctness
-  call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
@@ -273,8 +274,8 @@ subroutine  psb_mhalov(x,desc_a,info,work,tran,mode,data)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me
-  integer(psb_ipk_) :: err_act, ldx, &
-       & m, n, iix, jjx, ix, ijx, nrow, imode, err, liwork,data_
+  integer(psb_ipk_) :: err_act, ldx, iix, jjx, nrow, imode, err, liwork,data_
+  integer(psb_lpk_) :: m, n, ix, ijx
   integer(psb_mpk_),pointer :: iwork(:)
   character                :: tran_
   character(len=20)        :: name, ch_err
@@ -319,7 +320,7 @@ subroutine  psb_mhalov(x,desc_a,info,work,tran,mode,data)
   endif
   ldx = size(x,1)
   ! check vector correctness
-  call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
diff --git a/base/comm/psb_movrl_a.f90 b/base/comm/psb_movrl_a.f90
index 3df061c4c..95f8f0343 100644
--- a/base/comm/psb_movrl_a.f90
+++ b/base/comm/psb_movrl_a.f90
@@ -77,8 +77,9 @@ subroutine  psb_movrlm(x,desc_a,info,jx,ik,work,update,mode)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me
-  integer(psb_ipk_) :: err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, maxk, update_,&
+  integer(psb_ipk_) :: err_act, iix, jjx, nrow, ncol, k, maxk, update_,&
        & mode_, err, liwork, ldx
+  integer(psb_lpk_) :: m, n, ix, ijx
   integer(psb_mpk_),pointer :: iwork(:), xp(:,:)
   logical                  :: do_swap
   character(len=20)        :: name, ch_err
@@ -137,7 +138,7 @@ subroutine  psb_movrlm(x,desc_a,info,jx,ik,work,update,mode)
   do_swap = (mode_ /= 0)
   ldx = size(x,1)
   ! check vector correctness
-  call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
@@ -271,9 +272,9 @@ subroutine  psb_movrlv(x,desc_a,info,work,update,mode)
   integer(psb_ipk_), intent(in), optional             :: update,mode
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me, &
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, update_,&
-       & mode_, err, liwork, ldx
+  integer(psb_ipk_) :: ictxt, np, me, err_act, iix, jjx, nrow, ncol, &
+       & k, update_, mode_, err, liwork, ldx
+  integer(psb_lpk_) :: m, n, ix, ijx
   integer(psb_mpk_),pointer :: iwork(:)
   logical                  :: do_swap
   character(len=20)        :: name, ch_err
@@ -318,7 +319,7 @@ subroutine  psb_movrlv(x,desc_a,info,work,update,mode)
   do_swap = (mode_ /= 0)
   ldx = size(x,1) 
   ! check vector correctness
-  call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
diff --git a/base/comm/psb_sgather.f90 b/base/comm/psb_sgather.f90
index bb1c0390e..13192b99e 100644
--- a/base/comm/psb_sgather.f90
+++ b/base/comm/psb_sgather.f90
@@ -58,8 +58,8 @@ subroutine  psb_sgather_vect(globx, locx, desc_a, info, iroot)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
-  integer(psb_ipk_) :: ierr(5), err_act, n, ilocx, iglobx, jlocx,&
-       & jglobx, lda_locx, lda_globx, m, k, jlx, ilx, i, idx
+  integer(psb_ipk_) :: ierr(5), err_act, jlx, ilx, lda_locx, lda_globx, i
+  integer(psb_lpk_) :: m, n, k, ilocx,  jlocx, idx, iglobx, jglobx
   real(psb_spk_), allocatable :: llocx(:)
   character(len=20)        :: name, ch_err
 
@@ -171,8 +171,8 @@ subroutine  psb_sgather_multivect(globx, locx, desc_a, info, iroot)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
-  integer(psb_ipk_) :: ierr(5), err_act, n, ilocx, iglobx, jlocx,&
-       & jglobx, lda_locx, lda_globx, m, k, jlx, ilx, i, idx
+  integer(psb_ipk_) :: ierr(5), err_act, jlx, ilx, lda_locx, lda_globx, i 
+  integer(psb_lpk_) :: m, n, k, ilocx, jlocx, idx, iglobx, jglobx
   real(psb_spk_), allocatable :: llocx(:,:)
   character(len=20)        :: name, ch_err
 
diff --git a/base/comm/psb_sgather_a.f90 b/base/comm/psb_sgather_a.f90
index 7d6ba531c..62146841b 100644
--- a/base/comm/psb_sgather_a.f90
+++ b/base/comm/psb_sgather_a.f90
@@ -58,9 +58,9 @@ subroutine  psb_sgatherm(globx, locx, desc_a, info, iroot)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
-  integer(psb_ipk_) :: ierr(5), err_act, n, ilocx, iglobx, jlocx,&
-       & jglobx, lda_locx, lda_globx, m, lock, globk, maxk, k, jlx, &
-       & ilx, i, j, idx
+  integer(psb_ipk_) :: ierr(5), err_act, lda_locx, lda_globx, lock, globk,&
+       & maxk, k, jlx, ilx, i, j
+  integer(psb_lpk_) :: m, n, ilocx,  jlocx, idx, iglobx, jglobx
 
   character(len=20)        :: name, ch_err
 
@@ -230,8 +230,9 @@ subroutine  psb_sgatherv(globx, locx, desc_a, info, iroot)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
-  integer(psb_ipk_) :: ierr(5), err_act, n, ilocx, iglobx, jlocx,&
-       & jglobx, lda_locx, lda_globx, m, k, jlx, ilx, i, idx
+  integer(psb_ipk_) :: ierr(5), err_act, lda_locx, lda_globx, lock, globk,&
+       & maxk, k, jlx, ilx, i, j
+  integer(psb_lpk_) :: m, n, ilocx,  jlocx, idx, iglobx, jglobx
 
   character(len=20)        :: name, ch_err
 
diff --git a/base/comm/psb_shalo.f90 b/base/comm/psb_shalo.f90
index 657d27372..444c43fe9 100644
--- a/base/comm/psb_shalo.f90
+++ b/base/comm/psb_shalo.f90
@@ -66,9 +66,9 @@ subroutine  psb_shalo_vect(x,desc_a,info,work,tran,mode,data)
   character, intent(in), optional         :: tran
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me,&
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, imode,&
-       & err, liwork,data_
+  integer(psb_ipk_) :: ictxt, np, me, err_act, iix, jjx, &
+       & nrow, imode, err, liwork,data_
+  integer(psb_lpk_) :: m, n, ix, ijx
   real(psb_spk_),pointer :: iwork(:)
   character                 :: tran_
   character(len=20)         :: name, ch_err
@@ -119,7 +119,7 @@ subroutine  psb_shalo_vect(x,desc_a,info,work,tran,mode,data)
   endif
 
   ! check vector correctness
-  call psb_chkvect(m,ione,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
@@ -205,9 +205,9 @@ subroutine  psb_shalo_multivect(x,desc_a,info,work,tran,mode,data)
   character, intent(in), optional         :: tran
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me,&
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, imode,&
-       & err, liwork,data_
+  integer(psb_ipk_) :: ictxt, np, me, err_act, iix, jjx, &
+       & nrow, imode, err, liwork,data_
+  integer(psb_lpk_) :: m, n, ix, ijx
   real(psb_spk_),pointer :: iwork(:)
   character                 :: tran_
   character(len=20)         :: name, ch_err
@@ -258,7 +258,7 @@ subroutine  psb_shalo_multivect(x,desc_a,info,work,tran,mode,data)
   endif
 
   ! check vector correctness
-  call psb_chkvect(m,ione,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
diff --git a/base/comm/psb_shalo_a.f90 b/base/comm/psb_shalo_a.f90
index 99f89d266..5e9b83545 100644
--- a/base/comm/psb_shalo_a.f90
+++ b/base/comm/psb_shalo_a.f90
@@ -66,8 +66,9 @@ subroutine  psb_shalom(x,desc_a,info,jx,ik,work,tran,mode,data)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me
-  integer(psb_ipk_) :: err_act, m, n, iix, jjx, ix, ijx, k, maxk, nrow, imode, i,&
+  integer(psb_ipk_) :: err_act, iix, jjx, k, maxk, nrow, imode, i,&
        & err, liwork,data_, ldx
+  integer(psb_lpk_) :: m, n, ix, ijx
   real(psb_spk_),pointer :: iwork(:), xp(:,:)
   character                :: tran_
   character(len=20)        :: name, ch_err
@@ -129,7 +130,7 @@ subroutine  psb_shalom(x,desc_a,info,jx,ik,work,tran,mode,data)
   endif
   ldx = size(x,1)
   ! check vector correctness
-  call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
@@ -273,8 +274,8 @@ subroutine  psb_shalov(x,desc_a,info,work,tran,mode,data)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me
-  integer(psb_ipk_) :: err_act, ldx, &
-       & m, n, iix, jjx, ix, ijx, nrow, imode, err, liwork,data_
+  integer(psb_ipk_) :: err_act, ldx, iix, jjx, nrow, imode, err, liwork,data_
+  integer(psb_lpk_) :: m, n, ix, ijx
   real(psb_spk_),pointer :: iwork(:)
   character                :: tran_
   character(len=20)        :: name, ch_err
@@ -319,7 +320,7 @@ subroutine  psb_shalov(x,desc_a,info,work,tran,mode,data)
   endif
   ldx = size(x,1)
   ! check vector correctness
-  call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
diff --git a/base/comm/psb_sovrl.f90 b/base/comm/psb_sovrl.f90
index 489fa0eb7..83f9dcc5a 100644
--- a/base/comm/psb_sovrl.f90
+++ b/base/comm/psb_sovrl.f90
@@ -75,9 +75,9 @@ subroutine  psb_sovrl_vect(x,desc_a,info,work,update,mode)
   integer(psb_ipk_), intent(in), optional          :: update,mode
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me, &
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, update_,&
-       & mode_, err, liwork,ldx
+  integer(psb_ipk_) :: ictxt, np, me, err_act, k, iix, jjx, &
+       & nrow, imode, err, liwork,data_, update_, mode_, ncol
+  integer(psb_lpk_) :: m, n, ix, ijx
   real(psb_spk_),pointer :: iwork(:)
   logical                  :: do_swap
   character(len=20)        :: name, ch_err
@@ -127,7 +127,7 @@ subroutine  psb_sovrl_vect(x,desc_a,info,work,update,mode)
   do_swap = (mode_ /= 0)
 
   ! check vector correctness
-  call psb_chkvect(m,ione,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
@@ -200,9 +200,9 @@ subroutine  psb_sovrl_multivect(x,desc_a,info,work,update,mode)
   integer(psb_ipk_), intent(in), optional          :: update,mode
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me, &
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, update_,&
-       & mode_, err, liwork,ldx
+  integer(psb_ipk_) :: ictxt, np, me, err_act, k, iix, jjx, &
+       & nrow, imode, err, liwork,data_, update_, mode_, ncol
+  integer(psb_lpk_) :: m, n, ix, ijx
   real(psb_spk_),pointer :: iwork(:)
   logical                  :: do_swap
   character(len=20)        :: name, ch_err
@@ -252,7 +252,7 @@ subroutine  psb_sovrl_multivect(x,desc_a,info,work,update,mode)
   do_swap = (mode_ /= 0)
 
   ! check vector correctness
-  call psb_chkvect(m,ione,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
diff --git a/base/comm/psb_sovrl_a.f90 b/base/comm/psb_sovrl_a.f90
index 9dd994985..e12a4c6a2 100644
--- a/base/comm/psb_sovrl_a.f90
+++ b/base/comm/psb_sovrl_a.f90
@@ -77,8 +77,9 @@ subroutine  psb_sovrlm(x,desc_a,info,jx,ik,work,update,mode)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me
-  integer(psb_ipk_) :: err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, maxk, update_,&
+  integer(psb_ipk_) :: err_act, iix, jjx, nrow, ncol, k, maxk, update_,&
        & mode_, err, liwork, ldx
+  integer(psb_lpk_) :: m, n, ix, ijx
   real(psb_spk_),pointer :: iwork(:), xp(:,:)
   logical                  :: do_swap
   character(len=20)        :: name, ch_err
@@ -137,7 +138,7 @@ subroutine  psb_sovrlm(x,desc_a,info,jx,ik,work,update,mode)
   do_swap = (mode_ /= 0)
   ldx = size(x,1)
   ! check vector correctness
-  call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
@@ -271,9 +272,9 @@ subroutine  psb_sovrlv(x,desc_a,info,work,update,mode)
   integer(psb_ipk_), intent(in), optional             :: update,mode
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me, &
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, update_,&
-       & mode_, err, liwork, ldx
+  integer(psb_ipk_) :: ictxt, np, me, err_act, iix, jjx, nrow, ncol, &
+       & k, update_, mode_, err, liwork, ldx
+  integer(psb_lpk_) :: m, n, ix, ijx
   real(psb_spk_),pointer :: iwork(:)
   logical                  :: do_swap
   character(len=20)        :: name, ch_err
@@ -318,7 +319,7 @@ subroutine  psb_sovrlv(x,desc_a,info,work,update,mode)
   do_swap = (mode_ /= 0)
   ldx = size(x,1) 
   ! check vector correctness
-  call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
diff --git a/base/comm/psb_sspgather.F90 b/base/comm/psb_sspgather.F90
index 62b7ba445..055ebb049 100644
--- a/base/comm/psb_sspgather.F90
+++ b/base/comm/psb_sspgather.F90
@@ -51,6 +51,7 @@ subroutine  psb_ssp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
   logical, intent(in), optional   :: keepnum,keeploc
 
   type(psb_s_coo_sparse_mat)      :: loc_coo, glob_coo
+  integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:)
   integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg
   integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl
   logical :: keepnum_, keeploc_
@@ -101,14 +102,19 @@ subroutine  psb_ssp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
       call loca%mv_to(loc_coo)
     end if
     nzl = loc_coo%get_nzeros()
-    call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I')
-    call psb_loc_to_glob(loc_coo%ja(1:nzl),desc_a,info,iact='I')
+    call psb_realloc(nzl,locia,info)
+    call psb_realloc(nzl,locja,info)    
+    call psb_loc_to_glob(loc_coo%ia(1:nzl),locia(1:nzl),desc_a,info,iact='I')
+    call psb_loc_to_glob(loc_coo%ja(1:nzl),locja(1:nzl),desc_a,info,iact='I')
     nzbr(:) = 0
     nzbr(me+1) = nzl
     call psb_sum(ictxt,nzbr(1:np))
     nzg = sum(nzbr)
+    if (info == psb_success_) call psb_realloc(nzg,glbia,info)
+    if (info == psb_success_) call psb_realloc(nzg,glbja,info)    
     if (info == psb_success_) call glob_coo%allocate(nrg,ncg,nzg)
     if (info /= psb_success_) goto 9999
+
     do ip=1,np
       idisp(ip) = sum(nzbr(1:ip-1))
     enddo
@@ -117,24 +123,33 @@ subroutine  psb_ssp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
          & glob_coo%val,nzbr,idisp,&
          & psb_mpi_r_spk_,icomm,minfo)
     if (minfo == psb_success_) call &
-         & mpi_allgatherv(loc_coo%ia,ndx,psb_mpi_ipk_,&
-         & glob_coo%ia,nzbr,idisp,&
-         & psb_mpi_ipk_,icomm,minfo)
+         & mpi_allgatherv(locia,ndx,psb_mpi_lpk_,&
+         & glbia,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
     if (minfo == psb_success_) call &
-         & mpi_allgatherv(loc_coo%ja,ndx,psb_mpi_ipk_,&
-         & glob_coo%ja,nzbr,idisp,&
-         & psb_mpi_ipk_,icomm,minfo)
+         & mpi_allgatherv(locja,ndx,psb_mpi_lpk_,&
+         & glbja,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
     
     if (minfo /= psb_success_) then 
       info  = minfo
       call psb_errpush(psb_err_internal_error_,name,a_err=' from mpi_allgatherv')
       goto 9999
-    end if
-    
+    end if    
     call loc_coo%free()
+    deallocate(locia,locja, stat=info)
+    !
+    ! Is the code below safe? For very large cases
+    ! the indices in glob_coo will overflow. But then,
+    ! for very large cases it does not make sense to
+    ! gather the matrix on a single procecss anyway...
+    !
+    glob_coo%ia(1:nzg) = glbia(1:nzg)
+    glob_coo%ja(1:nzg) = glbja(1:nzg)
     call glob_coo%set_nzeros(nzg)
     if (present(dupl)) call glob_coo%set_dupl(dupl)
     call globa%mv_from(glob_coo)
+    deallocate(glbia,glbja, stat=info)
 
   else
     write(psb_err_unit,*) 'SP_ALLGATHER: Not implemented yet with keepnum ',keepnum_
diff --git a/base/comm/psb_zgather.f90 b/base/comm/psb_zgather.f90
index b4e919a45..59b6df405 100644
--- a/base/comm/psb_zgather.f90
+++ b/base/comm/psb_zgather.f90
@@ -58,8 +58,8 @@ subroutine  psb_zgather_vect(globx, locx, desc_a, info, iroot)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
-  integer(psb_ipk_) :: ierr(5), err_act, n, ilocx, iglobx, jlocx,&
-       & jglobx, lda_locx, lda_globx, m, k, jlx, ilx, i, idx
+  integer(psb_ipk_) :: ierr(5), err_act, jlx, ilx, lda_locx, lda_globx, i
+  integer(psb_lpk_) :: m, n, k, ilocx,  jlocx, idx, iglobx, jglobx
   complex(psb_dpk_), allocatable :: llocx(:)
   character(len=20)        :: name, ch_err
 
@@ -171,8 +171,8 @@ subroutine  psb_zgather_multivect(globx, locx, desc_a, info, iroot)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
-  integer(psb_ipk_) :: ierr(5), err_act, n, ilocx, iglobx, jlocx,&
-       & jglobx, lda_locx, lda_globx, m, k, jlx, ilx, i, idx
+  integer(psb_ipk_) :: ierr(5), err_act, jlx, ilx, lda_locx, lda_globx, i 
+  integer(psb_lpk_) :: m, n, k, ilocx, jlocx, idx, iglobx, jglobx
   complex(psb_dpk_), allocatable :: llocx(:,:)
   character(len=20)        :: name, ch_err
 
diff --git a/base/comm/psb_zgather_a.f90 b/base/comm/psb_zgather_a.f90
index 71b4f840f..014564872 100644
--- a/base/comm/psb_zgather_a.f90
+++ b/base/comm/psb_zgather_a.f90
@@ -58,9 +58,9 @@ subroutine  psb_zgatherm(globx, locx, desc_a, info, iroot)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
-  integer(psb_ipk_) :: ierr(5), err_act, n, ilocx, iglobx, jlocx,&
-       & jglobx, lda_locx, lda_globx, m, lock, globk, maxk, k, jlx, &
-       & ilx, i, j, idx
+  integer(psb_ipk_) :: ierr(5), err_act, lda_locx, lda_globx, lock, globk,&
+       & maxk, k, jlx, ilx, i, j
+  integer(psb_lpk_) :: m, n, ilocx,  jlocx, idx, iglobx, jglobx
 
   character(len=20)        :: name, ch_err
 
@@ -230,8 +230,9 @@ subroutine  psb_zgatherv(globx, locx, desc_a, info, iroot)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
-  integer(psb_ipk_) :: ierr(5), err_act, n, ilocx, iglobx, jlocx,&
-       & jglobx, lda_locx, lda_globx, m, k, jlx, ilx, i, idx
+  integer(psb_ipk_) :: ierr(5), err_act, lda_locx, lda_globx, lock, globk,&
+       & maxk, k, jlx, ilx, i, j
+  integer(psb_lpk_) :: m, n, ilocx,  jlocx, idx, iglobx, jglobx
 
   character(len=20)        :: name, ch_err
 
diff --git a/base/comm/psb_zhalo.f90 b/base/comm/psb_zhalo.f90
index a9f4f06c0..3d9db06ca 100644
--- a/base/comm/psb_zhalo.f90
+++ b/base/comm/psb_zhalo.f90
@@ -66,9 +66,9 @@ subroutine  psb_zhalo_vect(x,desc_a,info,work,tran,mode,data)
   character, intent(in), optional         :: tran
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me,&
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, imode,&
-       & err, liwork,data_
+  integer(psb_ipk_) :: ictxt, np, me, err_act, iix, jjx, &
+       & nrow, imode, err, liwork,data_
+  integer(psb_lpk_) :: m, n, ix, ijx
   complex(psb_dpk_),pointer :: iwork(:)
   character                 :: tran_
   character(len=20)         :: name, ch_err
@@ -119,7 +119,7 @@ subroutine  psb_zhalo_vect(x,desc_a,info,work,tran,mode,data)
   endif
 
   ! check vector correctness
-  call psb_chkvect(m,ione,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
@@ -205,9 +205,9 @@ subroutine  psb_zhalo_multivect(x,desc_a,info,work,tran,mode,data)
   character, intent(in), optional         :: tran
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me,&
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, imode,&
-       & err, liwork,data_
+  integer(psb_ipk_) :: ictxt, np, me, err_act, iix, jjx, &
+       & nrow, imode, err, liwork,data_
+  integer(psb_lpk_) :: m, n, ix, ijx
   complex(psb_dpk_),pointer :: iwork(:)
   character                 :: tran_
   character(len=20)         :: name, ch_err
@@ -258,7 +258,7 @@ subroutine  psb_zhalo_multivect(x,desc_a,info,work,tran,mode,data)
   endif
 
   ! check vector correctness
-  call psb_chkvect(m,ione,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
diff --git a/base/comm/psb_zhalo_a.f90 b/base/comm/psb_zhalo_a.f90
index 7dd7fbfe3..6fe5d70af 100644
--- a/base/comm/psb_zhalo_a.f90
+++ b/base/comm/psb_zhalo_a.f90
@@ -66,8 +66,9 @@ subroutine  psb_zhalom(x,desc_a,info,jx,ik,work,tran,mode,data)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me
-  integer(psb_ipk_) :: err_act, m, n, iix, jjx, ix, ijx, k, maxk, nrow, imode, i,&
+  integer(psb_ipk_) :: err_act, iix, jjx, k, maxk, nrow, imode, i,&
        & err, liwork,data_, ldx
+  integer(psb_lpk_) :: m, n, ix, ijx
   complex(psb_dpk_),pointer :: iwork(:), xp(:,:)
   character                :: tran_
   character(len=20)        :: name, ch_err
@@ -129,7 +130,7 @@ subroutine  psb_zhalom(x,desc_a,info,jx,ik,work,tran,mode,data)
   endif
   ldx = size(x,1)
   ! check vector correctness
-  call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
@@ -273,8 +274,8 @@ subroutine  psb_zhalov(x,desc_a,info,work,tran,mode,data)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me
-  integer(psb_ipk_) :: err_act, ldx, &
-       & m, n, iix, jjx, ix, ijx, nrow, imode, err, liwork,data_
+  integer(psb_ipk_) :: err_act, ldx, iix, jjx, nrow, imode, err, liwork,data_
+  integer(psb_lpk_) :: m, n, ix, ijx
   complex(psb_dpk_),pointer :: iwork(:)
   character                :: tran_
   character(len=20)        :: name, ch_err
@@ -319,7 +320,7 @@ subroutine  psb_zhalov(x,desc_a,info,work,tran,mode,data)
   endif
   ldx = size(x,1)
   ! check vector correctness
-  call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
diff --git a/base/comm/psb_zovrl.f90 b/base/comm/psb_zovrl.f90
index 8a1782b49..033bd0179 100644
--- a/base/comm/psb_zovrl.f90
+++ b/base/comm/psb_zovrl.f90
@@ -75,9 +75,9 @@ subroutine  psb_zovrl_vect(x,desc_a,info,work,update,mode)
   integer(psb_ipk_), intent(in), optional          :: update,mode
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me, &
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, update_,&
-       & mode_, err, liwork,ldx
+  integer(psb_ipk_) :: ictxt, np, me, err_act, k, iix, jjx, &
+       & nrow, imode, err, liwork,data_, update_, mode_, ncol
+  integer(psb_lpk_) :: m, n, ix, ijx
   complex(psb_dpk_),pointer :: iwork(:)
   logical                  :: do_swap
   character(len=20)        :: name, ch_err
@@ -127,7 +127,7 @@ subroutine  psb_zovrl_vect(x,desc_a,info,work,update,mode)
   do_swap = (mode_ /= 0)
 
   ! check vector correctness
-  call psb_chkvect(m,ione,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
@@ -200,9 +200,9 @@ subroutine  psb_zovrl_multivect(x,desc_a,info,work,update,mode)
   integer(psb_ipk_), intent(in), optional          :: update,mode
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me, &
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, update_,&
-       & mode_, err, liwork,ldx
+  integer(psb_ipk_) :: ictxt, np, me, err_act, k, iix, jjx, &
+       & nrow, imode, err, liwork,data_, update_, mode_, ncol
+  integer(psb_lpk_) :: m, n, ix, ijx
   complex(psb_dpk_),pointer :: iwork(:)
   logical                  :: do_swap
   character(len=20)        :: name, ch_err
@@ -252,7 +252,7 @@ subroutine  psb_zovrl_multivect(x,desc_a,info,work,update,mode)
   do_swap = (mode_ /= 0)
 
   ! check vector correctness
-  call psb_chkvect(m,ione,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
diff --git a/base/comm/psb_zovrl_a.f90 b/base/comm/psb_zovrl_a.f90
index 2e03c72f5..00c084a62 100644
--- a/base/comm/psb_zovrl_a.f90
+++ b/base/comm/psb_zovrl_a.f90
@@ -77,8 +77,9 @@ subroutine  psb_zovrlm(x,desc_a,info,jx,ik,work,update,mode)
 
   ! locals
   integer(psb_mpk_) :: ictxt, np, me
-  integer(psb_ipk_) :: err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, maxk, update_,&
+  integer(psb_ipk_) :: err_act, iix, jjx, nrow, ncol, k, maxk, update_,&
        & mode_, err, liwork, ldx
+  integer(psb_lpk_) :: m, n, ix, ijx
   complex(psb_dpk_),pointer :: iwork(:), xp(:,:)
   logical                  :: do_swap
   character(len=20)        :: name, ch_err
@@ -137,7 +138,7 @@ subroutine  psb_zovrlm(x,desc_a,info,jx,ik,work,update,mode)
   do_swap = (mode_ /= 0)
   ldx = size(x,1)
   ! check vector correctness
-  call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
@@ -271,9 +272,9 @@ subroutine  psb_zovrlv(x,desc_a,info,work,update,mode)
   integer(psb_ipk_), intent(in), optional             :: update,mode
 
   ! locals
-  integer(psb_ipk_) :: ictxt, np, me, &
-       & err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, update_,&
-       & mode_, err, liwork, ldx
+  integer(psb_ipk_) :: ictxt, np, me, err_act, iix, jjx, nrow, ncol, &
+       & k, update_, mode_, err, liwork, ldx
+  integer(psb_lpk_) :: m, n, ix, ijx
   complex(psb_dpk_),pointer :: iwork(:)
   logical                  :: do_swap
   character(len=20)        :: name, ch_err
@@ -318,7 +319,7 @@ subroutine  psb_zovrlv(x,desc_a,info,work,update,mode)
   do_swap = (mode_ /= 0)
   ldx = size(x,1) 
   ! check vector correctness
-  call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
+  call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
     ch_err='psb_chkvect'
diff --git a/base/comm/psb_zspgather.F90 b/base/comm/psb_zspgather.F90
index 777bf71fb..5e62017b9 100644
--- a/base/comm/psb_zspgather.F90
+++ b/base/comm/psb_zspgather.F90
@@ -51,6 +51,7 @@ subroutine  psb_zsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
   logical, intent(in), optional   :: keepnum,keeploc
 
   type(psb_z_coo_sparse_mat)      :: loc_coo, glob_coo
+  integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:)
   integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg
   integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl
   logical :: keepnum_, keeploc_
@@ -101,14 +102,19 @@ subroutine  psb_zsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
       call loca%mv_to(loc_coo)
     end if
     nzl = loc_coo%get_nzeros()
-    call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I')
-    call psb_loc_to_glob(loc_coo%ja(1:nzl),desc_a,info,iact='I')
+    call psb_realloc(nzl,locia,info)
+    call psb_realloc(nzl,locja,info)    
+    call psb_loc_to_glob(loc_coo%ia(1:nzl),locia(1:nzl),desc_a,info,iact='I')
+    call psb_loc_to_glob(loc_coo%ja(1:nzl),locja(1:nzl),desc_a,info,iact='I')
     nzbr(:) = 0
     nzbr(me+1) = nzl
     call psb_sum(ictxt,nzbr(1:np))
     nzg = sum(nzbr)
+    if (info == psb_success_) call psb_realloc(nzg,glbia,info)
+    if (info == psb_success_) call psb_realloc(nzg,glbja,info)    
     if (info == psb_success_) call glob_coo%allocate(nrg,ncg,nzg)
     if (info /= psb_success_) goto 9999
+
     do ip=1,np
       idisp(ip) = sum(nzbr(1:ip-1))
     enddo
@@ -117,24 +123,33 @@ subroutine  psb_zsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
          & glob_coo%val,nzbr,idisp,&
          & psb_mpi_c_dpk_,icomm,minfo)
     if (minfo == psb_success_) call &
-         & mpi_allgatherv(loc_coo%ia,ndx,psb_mpi_ipk_,&
-         & glob_coo%ia,nzbr,idisp,&
-         & psb_mpi_ipk_,icomm,minfo)
+         & mpi_allgatherv(locia,ndx,psb_mpi_lpk_,&
+         & glbia,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
     if (minfo == psb_success_) call &
-         & mpi_allgatherv(loc_coo%ja,ndx,psb_mpi_ipk_,&
-         & glob_coo%ja,nzbr,idisp,&
-         & psb_mpi_ipk_,icomm,minfo)
+         & mpi_allgatherv(locja,ndx,psb_mpi_lpk_,&
+         & glbja,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
     
     if (minfo /= psb_success_) then 
       info  = minfo
       call psb_errpush(psb_err_internal_error_,name,a_err=' from mpi_allgatherv')
       goto 9999
-    end if
-    
+    end if    
     call loc_coo%free()
+    deallocate(locia,locja, stat=info)
+    !
+    ! Is the code below safe? For very large cases
+    ! the indices in glob_coo will overflow. But then,
+    ! for very large cases it does not make sense to
+    ! gather the matrix on a single procecss anyway...
+    !
+    glob_coo%ia(1:nzg) = glbia(1:nzg)
+    glob_coo%ja(1:nzg) = glbja(1:nzg)
     call glob_coo%set_nzeros(nzg)
     if (present(dupl)) call glob_coo%set_dupl(dupl)
     call globa%mv_from(glob_coo)
+    deallocate(glbia,glbja, stat=info)
 
   else
     write(psb_err_unit,*) 'SP_ALLGATHER: Not implemented yet with keepnum ',keepnum_