psblas3:

Make.inc.in
 base/Makefile
 base/modules/fakempi.c
 base/modules/psb_penv_mod.F90
 base/newserial/psbn_base_mat_mod.f03
 base/newserial/psbn_coo_mat.f03
 base/newserial/psbn_csr_mat.f03
 base/newserial/psbn_d_base_mat_mod.f03
 config/ax_f90_module_extension.m4
 config/ax_f90_module_flag.m4
 config/pac.m4
 configure.ac
 configure

Added  --with-serial-mpi.

Make.inc.in

@@ -15,6 +15,7 @@ F90=@F90@
 FC=@FC@
 CC=@CC@
 #F77=@F77@
+F03COPT=@F03COPT@
 F90COPT=@F90COPT@
 FCOPT=@FCOPT@
 CCOPT=@CCOPT@
@@ -37,6 +38,7 @@ BLAS=@BLAS_LIBS@
 BLACS=@BLACS_LIBS@
 METIS_LIB=@METIS_LIBS@
 LAPACK=@LAPACK_LIBS@
+EXTRA_COBJS=@FAKEMPI@
 
 PSBFDEFINES=@FDEFINES@
 PSBCDEFINES=@CDEFINES@

base/Makefile

@@ -6,11 +6,11 @@ LIBNAME=$(BASELIBNAME)
 LIBMOD=psb_base_mod$(.mod)
 lib:
 	(cd modules; make lib LIBNAME=$(BASELIBNAME) F90=$(MPF90) F90COPT="$(F90COPT) $(MPI_OPT)")
+	(cd serial; make lib LIBNAME=$(BASELIBNAME))
 	(cd newserial; make lib LIBNAME=$(BASELIBNAME) )
 	(cd comm; make lib LIBNAME=$(BASELIBNAME))
 	(cd internals; make lib LIBNAME=$(BASELIBNAME))
 	(cd tools; make lib LIBNAME=$(BASELIBNAME))
-	(cd serial; make lib LIBNAME=$(BASELIBNAME))
 	(cd psblas; make lib LIBNAME=$(BASELIBNAME))
 	/bin/cp -p $(HERE)/$(LIBNAME) $(LIBDIR)
 	/bin/cp -p $(LIBMOD) *$(.mod) $(LIBDIR)

base/modules/fakempi.c

@@ -64,8 +64,11 @@
 #endif
 
 #define mpi_integer        1
-#define mpi_double         3
+#define mpi_integer8       2
-#define mpi_double_complex 5
+#define mpi_real           3
+#define mpi_double         4
+#define mpi_complex        5
+#define mpi_double_complex 6
 
 double mpi_wtime() 
 {
@@ -104,9 +107,18 @@ void mpi_alltoall(void* sdb, int* sdc, int* sdt,
   if (*sdt == mpi_integer) {
     memcpy(rvb,sdb, (*sdc)*sizeof(int));
   }
+  if (*sdt == mpi_integer8) {
+    memcpy(rvb,sdb, (*sdc)*2*sizeof(int));
+  }
+  if (*sdt == mpi_real) {
+    memcpy(rvb,sdb, (*sdc)*sizeof(float));
+  } 
   if (*sdt == mpi_double) {
     memcpy(rvb,sdb, (*sdc)*sizeof(double));
   } 
+  if (*sdt == mpi_complex) {
+    memcpy(rvb,sdb, (*sdc)*2*sizeof(float));
+  }
   if (*sdt == mpi_double_complex) {
     memcpy(rvb,sdb, (*sdc)*2*sizeof(double));
   }
@@ -123,10 +135,22 @@ void mpi_alltoallv(void* sdb, int* sdc, int* sdspl, int* sdt,
     memcpy((rvb+rdspl[0]*sizeof(int)),
 	   (sdb+sdspl[0]*sizeof(int)),(*sdc)*sizeof(int));
   }
+  if (*sdt == mpi_integer8) {
+    memcpy((rvb+rdspl[0]*2*sizeof(int)),
+	   (sdb+sdspl[0]*2*sizeof(int)),(*sdc)*2*sizeof(int));
+  }
+  if (*sdt == mpi_real) {
+    memcpy((rvb+rdspl[0]*sizeof(float)),
+	   (sdb+sdspl[0]*sizeof(float)),(*sdc)*sizeof(float));
+  } 
   if (*sdt == mpi_double) {
     memcpy((rvb+rdspl[0]*sizeof(double)),
 	   (sdb+sdspl[0]*sizeof(double)),(*sdc)*sizeof(double));
   } 
+  if (*sdt == mpi_complex) {
+    memcpy((rvb+rdspl[0]*2*sizeof(float)),
+	   (sdb+sdspl[0]*2*sizeof(float)),(*sdc)*2*sizeof(float));
+  }
   if (*sdt == mpi_double_complex) {
     memcpy((rvb+rdspl[0]*2*sizeof(double)),
 	   (sdb+sdspl[0]*2*sizeof(double)),(*sdc)*2*sizeof(double));
@@ -143,9 +167,18 @@ void mpi_allgather(void* sdb, int* sdc, int* sdt,
   if (*sdt == mpi_integer) {
     memcpy(rvb,sdb, (*sdc)*sizeof(int));
   }
+  if (*sdt == mpi_integer8) {
+    memcpy(rvb,sdb, (*sdc)*2*sizeof(int));
+  }
+  if (*sdt == mpi_real) {
+    memcpy(rvb,sdb, (*sdc)*sizeof(float));
+  } 
   if (*sdt == mpi_double) {
     memcpy(rvb,sdb, (*sdc)*sizeof(double));
   } 
+  if (*sdt == mpi_complex) {
+    memcpy(rvb,sdb, (*sdc)*2*sizeof(float));
+  }
   if (*sdt == mpi_double_complex) {
     memcpy(rvb,sdb, (*sdc)*2*sizeof(double));
   }
@@ -158,23 +191,26 @@ void mpi_allgatherv(void* sdb, int* sdc, int* sdt,
 {
   int i,j,k; 
   
-  if (*sdt == mpi_integer) {
-    memcpy(rvb,sdb, (*sdc)*sizeof(int));
-  }
-  if (*sdt == mpi_double) {
-    memcpy(rvb,sdb, (*sdc)*sizeof(double));
-  } 
-  if (*sdt == mpi_double_complex) {
-    memcpy(rvb,sdb, (*sdc)*2*sizeof(double));
-  }
   if (*sdt == mpi_integer) {
     memcpy((rvb+rdspl[0]*sizeof(int)),
 	   (sdb),(*sdc)*sizeof(int));
   }
+  if (*sdt == mpi_integer8) {
+    memcpy((rvb+rdspl[0]*2*sizeof(int)),
+	   (sdb),(*sdc)*2*sizeof(int));
+  }
+  if (*sdt == mpi_real) {
+    memcpy((rvb+rdspl[0]*sizeof(float)),
+	   (sdb),(*sdc)*sizeof(float));
+  } 
   if (*sdt == mpi_double) {
     memcpy((rvb+rdspl[0]*sizeof(double)),
 	   (sdb),(*sdc)*sizeof(double));
   } 
+  if (*sdt == mpi_complex) {
+    memcpy((rvb+rdspl[0]*2*sizeof(float)),
+	   (sdb),(*sdc)*2*sizeof(float));
+  }
   if (*sdt == mpi_double_complex) {
     memcpy((rvb+rdspl[0]*2*sizeof(double)),
 	   (sdb),(*sdc)*2*sizeof(double));

base/modules/psb_penv_mod.F90

@@ -32,11 +32,16 @@
 #if defined(SERIAL_MPI)
 ! Provide a fake mpi module just to keep the compiler(s) happy.
 module mpi
+  use psb_const_mod
   integer, parameter :: mpi_success=0
   integer, parameter :: mpi_request_null=0
   integer, parameter :: mpi_status_size=1
-  integer, parameter :: mpi_integer=1, mpi_double_precision=3
+  integer, parameter :: mpi_integer          = 1
-  integer, parameter :: mpi_double_complex=5 
+  integer, parameter :: mpi_integer8         = 2
+  integer, parameter :: mpi_real             = 3
+  integer, parameter :: mpi_double_precision = 4
+  integer, parameter :: mpi_complex          = 5   
+  integer, parameter :: mpi_double_complex   = 6 
   real(psb_dpk_), external :: mpi_wtime
 end module mpi
 #endif    
@@ -144,9 +149,7 @@ module psb_penv_mod
 
 #if defined(SERIAL_MPI)
   integer, private, save :: nctxt=0
-#endif
+#else 
-
-
 #if defined(HAVE_KSENDID)
   interface 
     integer function krecvid(contxt,proc_to_comm,myrow)
@@ -159,6 +162,8 @@ module psb_penv_mod
     end function ksendid
   end interface
 #endif
+#endif
+  
   private psi_get_sizes
 contains 
 
@@ -2351,6 +2356,7 @@ contains
     else
       root_ = -1
     endif
+#if !defined(SERIAL_MPI)
     call psb_info(ictxt,iam,np)
     call psb_get_mpicomm(ictxt,icomm)
     mpi_int8_type = mpi_integer8
@@ -2367,6 +2373,7 @@ contains
         call mpi_reduce(dat,dat_,isz,mpi_int8_type,mpi_sum,root_,icomm,info)
       end if
     endif
+#endif
 
   end subroutine psb_i8sumv
 
@@ -2391,6 +2398,7 @@ contains
     else
       root_ = -1
     endif
+#if !defined(SERIAL_MPI)
     call psb_info(ictxt,iam,np)
     call psb_get_mpicomm(ictxt,icomm)
     mpi_int8_type = mpi_integer8
@@ -2405,7 +2413,7 @@ contains
         call mpi_reduce(dat,dat_,1,mpi_int8_type,mpi_sum,root_,icomm,info)
       end if
     endif
-
+#endif
   end subroutine psb_i8sums
 
   subroutine psb_i8amx_mpi_user(inv, outv,len,type) 
@@ -3502,7 +3510,7 @@ contains
   subroutine psb_set_coher(ictxt,isvch)
     integer :: ictxt, isvch
     ! Ensure global repeatability for convergence checks.
-#if !defined(HAVE_ESSL_BLACS)
+#if (!defined(HAVE_ESSL_BLACS)) &&(!defined(SERIAL_MPI))
     Call blacs_get(ictxt,15,isvch)
     Call blacs_set(ictxt,15,1)
 #else
@@ -3513,7 +3521,7 @@ contains
   subroutine psb_restore_coher(ictxt,isvch)
     integer :: ictxt, isvch
     ! Ensure global coherence for convergence checks.
-#if !defined(HAVE_ESSL_BLACS)
+#if (!defined(HAVE_ESSL_BLACS)) &&(!defined(SERIAL_MPI))
     Call blacs_set(ictxt,15,isvch)
 #else
     ! Do nothing: ESSL does coherence by default,

base/newserial/psbn_base_mat_mod.f03

@@ -1,5 +1,7 @@
 module psbn_base_mat_mod
   
+  use psb_const_mod 
+
   integer, parameter :: psbn_spmat_null_=0, psbn_spmat_bld_=1
   integer, parameter :: psbn_spmat_asb_=2, psbn_spmat_upd_=4
 
@@ -24,117 +26,117 @@ module psbn_base_mat_mod
 
   type  :: psbn_base_sparse_mat
     integer              :: m, n
-    integer, private     :: state 
+    integer, private     :: state, duplicate 
     logical, private     :: triangle, unitd, upper, sorted
   contains 
-    procedure, pass(a) :: base_get_nrows
+    procedure, pass(a) :: get_nrows
-    procedure, pass(a) :: base_get_ncols
+    procedure, pass(a) :: get_ncols
-    procedure, pass(a) :: base_get_nzeros
+    procedure, pass(a) :: get_nzeros
-    procedure, pass(a) :: base_get_size
+    procedure, pass(a) :: get_size
-    procedure, pass(a) :: base_get_state
+    procedure, pass(a) :: get_state
-    procedure, pass(a) :: base_is_bld
+    procedure, pass(a) :: get_dupl
-    procedure, pass(a) :: base_is_upd
+    procedure, pass(a) :: is_null
-    procedure, pass(a) :: base_is_asb
+    procedure, pass(a) :: is_bld
-    procedure, pass(a) :: base_is_sorted
+    procedure, pass(a) :: is_upd
-    procedure, pass(a) :: base_is_upper
+    procedure, pass(a) :: is_asb
-    procedure, pass(a) :: base_is_lower
+    procedure, pass(a) :: is_sorted
-    procedure, pass(a) :: base_is_triangle
+    procedure, pass(a) :: is_upper
-    procedure, pass(a) :: base_is_unit
+    procedure, pass(a) :: is_lower
-    procedure, pass(a) :: base_get_neigh
+    procedure, pass(a) :: is_triangle
-    procedure, pass(a) :: base_allocate_mn
+    procedure, pass(a) :: is_unit
-    procedure, pass(a) :: base_allocate_mnnz
+    procedure, pass(a) :: get_neigh
-    procedure, pass(a) :: base_reallocate_nz
+    procedure, pass(a) :: allocate_mn
-    procedure, pass(a) :: base_free
+    procedure, pass(a) :: allocate_mnnz
-    generic, public    :: allocate => base_allocate_mn, base_allocate_mnnz
+    procedure, pass(a) :: reallocate_nz
-    generic, public    :: reallocate => base_reallocate_nz
+    procedure, pass(a) :: free
-    generic, public    :: get_nrows => base_get_nrows
+    generic,   public  :: allocate => allocate_mn, allocate_mnnz
-    generic, public    :: get_ncols => base_get_ncols
+    generic,   public  :: reallocate => reallocate_nz
-    generic, public    :: get_nzeros => base_get_nzeros
-    generic, public    :: get_size  => base_get_size
-    generic, public    :: get_state => base_get_state
-    generic, public    :: is_triangle => base_is_triangle
-    generic, public    :: is_unit => base_is_unit
-    generic, public    :: is_upper => base_is_upper
-    generic, public    :: is_lower => base_is_lower
-    generic, public    :: is_bld => base_is_bld
-    generic, public    :: is_upd => base_is_upd
-    generic, public    :: is_asb => base_is_asb
-    generic, public    :: is_sorted => base_is_sorted
-    generic, public    :: get_neigh => base_get_neigh
-    generic, public    :: free => base_free
     
   end type psbn_base_sparse_mat
   
 contains
  
-  function base_get_state(a) result(res)
+  function get_dupl(a) result(res)
+    class(psbn_base_sparse_mat), intent(in) :: a
+    integer :: res
+    res = a%duplicate
+  end function get_dupl
+ 
+ 
+  function get_state(a) result(res)
     class(psbn_base_sparse_mat), intent(in) :: a
     integer :: res
     res = a%state
-  end function base_get_state
+  end function get_state
  
-  function base_get_nrows(a) result(res)
+  function get_nrows(a) result(res)
     class(psbn_base_sparse_mat), intent(in) :: a
     integer :: res
     res = a%m
-  end function base_get_nrows
+  end function get_nrows
 
-  function base_get_ncols(a) result(res)
+  function get_ncols(a) result(res)
     class(psbn_base_sparse_mat), intent(in) :: a
     integer :: res
     res = a%n
-  end function base_get_ncols
+  end function get_ncols
 
-  function base_is_triangle(a) result(res)
+  function is_triangle(a) result(res)
     class(psbn_base_sparse_mat), intent(in) :: a
     logical :: res
     res = a%triangle
-  end function base_is_triangle
+  end function is_triangle
 
-  function base_is_unit(a) result(res)
+  function is_unit(a) result(res)
     class(psbn_base_sparse_mat), intent(in) :: a
     logical :: res
     res = a%unitd
-  end function base_is_unit
+  end function is_unit
 
-  function base_is_upper(a) result(res)
+  function is_upper(a) result(res)
     class(psbn_base_sparse_mat), intent(in) :: a
     logical :: res
     res = a%upper
-  end function base_is_upper
+  end function is_upper
 
-  function base_is_lower(a) result(res)
+  function is_lower(a) result(res)
     class(psbn_base_sparse_mat), intent(in) :: a
     logical :: res
     res = .not.a%upper
-  end function base_is_lower
+  end function is_lower
+
+  function is_null(a) result(res)
+    class(psbn_base_sparse_mat), intent(in) :: a
+    logical :: res
+    res = (a%state == psbn_spmat_null_)
+  end function is_null
 
-  function base_is_bld(a) result(res)
+  function is_bld(a) result(res)
     class(psbn_base_sparse_mat), intent(in) :: a
     logical :: res
     res = (a%state == psbn_spmat_bld_)
-  end function base_is_bld
+  end function is_bld
 
-  function base_is_upd(a) result(res)
+  function is_upd(a) result(res)
     class(psbn_base_sparse_mat), intent(in) :: a
     logical :: res
     res = (a%state == psbn_spmat_upd_)
-  end function base_is_upd
+  end function is_upd
 
-  function base_is_asb(a) result(res)
+  function is_asb(a) result(res)
     class(psbn_base_sparse_mat), intent(in) :: a
     logical :: res
     res = (a%state == psbn_spmat_asb_)
-  end function base_is_asb
+  end function is_asb
 
-  function base_is_sorted(a) result(res)
+  function is_sorted(a) result(res)
     class(psbn_base_sparse_mat), intent(in) :: a
     logical :: res
     res = a%sorted
-  end function base_is_sorted
+  end function is_sorted
 
 
-  function base_get_nzeros(a) result(res)
+  function get_nzeros(a) result(res)
     use psb_error_mod
     class(psbn_base_sparse_mat), intent(in) :: a
     integer :: res
@@ -155,15 +157,15 @@ contains
     end if
     return
 
-  end function base_get_nzeros
+  end function get_nzeros
 
-  function base_get_size(a) result(res)
+  function get_size(a) result(res)
     use psb_error_mod
     class(psbn_base_sparse_mat), intent(in) :: a
     integer :: res
     
     Integer :: err_act
-    character(len=20)  :: name='base_get_size'
+    character(len=20)  :: name='get_size'
     logical, parameter :: debug=.false.
 
     call psb_erractionsave(err_act)
@@ -178,10 +180,10 @@ contains
     end if
     return
 
-  end function base_get_size
+  end function get_size
 
 
-  subroutine base_get_neigh(a,idx,neigh,n,info,lev)
+  subroutine get_neigh(a,idx,neigh,n,info,lev)
     use psb_error_mod
     class(psbn_base_sparse_mat), intent(in) :: a   
     integer, intent(in)                :: idx 
@@ -191,7 +193,7 @@ contains
     integer, optional, intent(in)      :: lev 
     
     Integer :: err_act
-    character(len=20)  :: name='base_get_neigh'
+    character(len=20)  :: name='get_neigh'
     logical, parameter :: debug=.false.
 
     call psb_erractionsave(err_act)
@@ -206,15 +208,15 @@ contains
     end if
     return
 
-  end subroutine base_get_neigh
+  end subroutine get_neigh
 
-  subroutine  base_allocate_mn(m,n,a) 
+  subroutine  allocate_mn(m,n,a) 
     use psb_error_mod
     integer, intent(in) :: m,n
     class(psbn_base_sparse_mat), intent(inout) :: a
 
     Integer :: err_act
-    character(len=20)  :: name='base_allocate_mn'
+    character(len=20)  :: name='allocate_mn'
     logical, parameter :: debug=.false.
 
     call psb_erractionsave(err_act)
@@ -228,14 +230,14 @@ contains
     end if
     return
 
-  end subroutine base_allocate_mn
+  end subroutine allocate_mn
 
-  subroutine  base_allocate_mnnz(m,n,nz,a) 
+  subroutine  allocate_mnnz(m,n,nz,a) 
     use psb_error_mod
     integer, intent(in) :: m,n,nz
     class(psbn_base_sparse_mat), intent(inout) :: a
     Integer :: err_act
-    character(len=20)  :: name='base_allocate_mnz'
+    character(len=20)  :: name='allocate_mnz'
     logical, parameter :: debug=.false.
 
     call psb_erractionsave(err_act)
@@ -249,14 +251,14 @@ contains
     end if
     return
 
-  end subroutine base_allocate_mnnz
+  end subroutine allocate_mnnz
 
-  subroutine  base_reallocate_nz(nz,a) 
+  subroutine  reallocate_nz(nz,a) 
     use psb_error_mod
     integer, intent(in) :: nz
     class(psbn_base_sparse_mat), intent(inout) :: a
     Integer :: err_act
-    character(len=20)  :: name='base_reallocate_nz'
+    character(len=20)  :: name='reallocate_nz'
     logical, parameter :: debug=.false.
 
     call psb_erractionsave(err_act)
@@ -270,13 +272,13 @@ contains
     end if
     return
 
-  end subroutine base_reallocate_nz
+  end subroutine reallocate_nz
 
-  subroutine  base_free(a) 
+  subroutine  free(a) 
     use psb_error_mod
     class(psbn_base_sparse_mat), intent(inout) :: a
     Integer :: err_act
-    character(len=20)  :: name='base_free'
+    character(len=20)  :: name='free'
     logical, parameter :: debug=.false.
 
     call psb_erractionsave(err_act)
@@ -290,7 +292,7 @@ contains
     end if
     return
 
-  end subroutine base_free
+  end subroutine free
 
 end module psbn_base_mat_mod
 

base/newserial/psbn_coo_mat.f03

@@ -3,20 +3,23 @@ module psbn_d_coo_sparse_mat_mod
   use psbn_d_base_mat_mod
 
   type, extends(psbn_d_base_sparse_mat) :: psbn_d_coo_sparse_mat
-    integer              :: nnz,  state
+
+    integer              :: nnz
     logical              :: sorted
     
     integer, allocatable :: ia(:), ja(:)
-    real(kind(1.d0)), allocatable :: val(:)
+    real(psb_dpk_), allocatable :: val(:)
+
   contains
-    procedure, pass(a)  :: d_coo_get_nzeros
+
+    procedure, pass(a)  :: get_nzeros => d_coo_get_nzeros
+    procedure, pass(a)  :: set_nzeros => d_coo_set_nzeros
     procedure, pass(a)  :: d_base_csmm => d_coo_csmm
     procedure, pass(a)  :: d_base_csmv => d_coo_csmv
-    generic, public     :: base_get_nzeros => d_coo_get_nzeros
     procedure, pass(a)  :: d_base_cssm => d_coo_cssm
     procedure, pass(a)  :: d_base_cssv => d_coo_cssv
-    procedure, pass(a)  :: d_base_csins => d_coo_csins
+    procedure, pass(a)  :: csins => d_coo_csins
-    procedure, pass(a)  :: base_reallocate_nz => d_coo_reallocate_nz
+    procedure, pass(a)  :: reallocate_nz => d_coo_reallocate_nz
 
   end type psbn_d_coo_sparse_mat
 
@@ -24,7 +27,6 @@ contains
 
 
   subroutine  d_coo_reallocate_nz(nz,a) 
-    use psb_const_mod
     use psb_error_mod
     use psb_realloc_mod
     integer, intent(in) :: nz
@@ -60,23 +62,33 @@ contains
   function d_coo_get_nzeros(a) result(res)
     class(psbn_d_coo_sparse_mat), intent(in) :: a
     integer :: res
-    res = a%nnz
+    res  = a%nnz
   end function d_coo_get_nzeros
   
 
-  subroutine d_coo_csins(nz,val,ia,ja,a,info) 
+  subroutine  d_coo_set_nzeros(nz,a)
-    use psb_const_mod
+    integer, intent(in) :: nz
+    class(psbn_d_coo_sparse_mat), intent(inout) :: a
+
+    a%nnz = nz
+
+  end subroutine d_coo_set_nzeros
+  
+
+  subroutine d_coo_csins(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl) 
     use psb_error_mod
     use psb_realloc_mod
-    class(psbn_d_coo_sparse_mat), intent(in) :: a
+    class(psbn_d_coo_sparse_mat), intent(inout) :: a
-    real(psb_dpk_), intent(in)    :: val(:)
+    real(psb_dpk_), intent(in)      :: val(:)
-    integer, intent(in)           :: nz, ia(:), ja(:)
+    integer, intent(in)             :: nz, ia(:), ja(:), imin,imax,jmin,jmax
     integer, intent(out)            :: info
+    integer, intent(in), optional   :: gtl(:)
 
-    Integer :: err_act
+
+    Integer            :: err_act
     character(len=20)  :: name='d_coo_csins'
     logical, parameter :: debug=.false.
-    integer :: nza, i,j,k, nzl, isza, int_err(5)
+    integer            :: nza, i,j,k, nzl, isza, int_err(5)
 
     call psb_erractionsave(err_act)
     info = 0
@@ -109,20 +121,43 @@ contains
 
     if (nz == 0) return
 
+
+    nza  = a%get_nzeros()
+    isza = a%get_size()
     if (a%is_bld()) then 
       ! Build phase. Must handle reallocations in a sensible way.
-      nza  = a%get_nzeros()
-      isza = a%get_size()
       if (isza < (nza+nz)) then 
-        
+        call a%reallocate(max(nza+nz,int(1.5*isza)))
-        
+        isza = a%get_size()
       endif
 
+      call psb_inner_ins(nz,ia,ja,val,nza,a%ia,a%ja,a%val,isza,&
+           & imin,imax,jmin,jmax,info,gtl)
+      call a%set_nzeros(nz+nza)
+
     else  if (a%is_upd()) then 
+      if (a%is_sorted()) then 
+
+
+!!$#ifdef FIXED_NAG_SEGV
+!!$        call  d_coo_srch_upd(nz,ia,ja,val,a,&
+!!$             & imin,imax,jmin,jmax,info,gtl)
+!!$#else 
+        call  d_coo_srch_upd(nz,ia,ja,val,&
+             & a%ia,a%ja,a%val,&
+             & a%get_dupl(),a%get_nzeros(),a%get_nrows(),&
+             & info,gtl)
+!!$#endif
+
+      else
+        info = 1121
+      end if
 
     else 
       ! State is wrong.
       info = 1121
+    end if
+    if (info /= 0) then
       call psb_errpush(info,name)
       goto 9999
     end if
@@ -141,86 +176,83 @@ contains
 
 
   contains
-
+!!$
-    subroutine psb_inner_upd(nz,ia,ja,val,nza,aspk,maxsz,&
+!!$    subroutine psb_inner_upd(nz,ia,ja,val,nza,aspk,maxsz,&
-         & imin,imax,jmin,jmax,nzl,info,gtl,ng)
+!!$         & imin,imax,jmin,jmax,info,gtl)
-      implicit none 
+!!$      implicit none 
-
+!!$
-      integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,maxsz
+!!$      integer, intent(in) :: nz, imin,imax,jmin,jmax,maxsz
-      integer, intent(in) :: ia(:),ja(:)
+!!$      integer, intent(in) :: ia(:),ja(:)
-      integer, intent(inout) :: nza
+!!$      integer, intent(inout) :: nza
-      real(psb_dpk_), intent(in) :: val(:)
+!!$      real(psb_dpk_), intent(in) :: val(:)
-      real(psb_dpk_), intent(inout) :: aspk(:)
+!!$      real(psb_dpk_), intent(inout) :: aspk(:)
-      integer, intent(out) :: info
+!!$      integer, intent(out) :: info
-      integer, intent(in), optional  :: ng,gtl(:)
+!!$      integer, intent(in), optional  :: gtl(:)
-      integer  :: i,ir,ic
+!!$      integer  :: i,ir,ic, ng,nzl
-      character(len=20)    :: name, ch_err
+!!$      character(len=20)    :: name, ch_err
-
+!!$
-
+!!$
-      name='psb_inner_upd'
+!!$      name='psb_inner_upd'
-
+!!$      nzl = 0
-      if (present(gtl)) then 
+!!$      if (present(gtl)) then 
-        if (.not.present(ng)) then 
+!!$        ng = size(gtl) 
-          info = -1
+!!$        if ((nza > nzl)) then 
-          return
+!!$          do i=1, nz 
-        endif
+!!$            nza = nza + 1 
-        if ((nza > nzl)) then 
+!!$            if (nza>maxsz) then 
-          do i=1, nz 
+!!$              call psb_errpush(50,name,i_err=(/7,maxsz,5,0,nza /))
-            nza = nza + 1 
+!!$              info = -71
-            if (nza>maxsz) then 
+!!$              return
-              call psb_errpush(50,name,i_err=(/7,maxsz,5,0,nza /))
+!!$            endif
-              info = -71
+!!$            aspk(nza) = val(i)
-              return
+!!$          end do
-            endif
+!!$        else
-            aspk(nza) = val(i)
+!!$          do i=1, nz 
-          end do
+!!$            ir = ia(i)
-        else
+!!$            ic = ja(i) 
-          do i=1, nz 
+!!$            if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then 
-            ir = ia(i)
+!!$              ir = gtl(ir)
-            ic = ja(i) 
+!!$              ic = gtl(ic) 
-            if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then 
+!!$              if ((ir >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then 
-              ir = gtl(ir)
+!!$                nza = nza + 1 
-              ic = gtl(ic) 
+!!$                if (nza>maxsz) then 
-              if ((ir >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then 
+!!$                  info = -72
-                nza = nza + 1 
+!!$                  return
-                if (nza>maxsz) then 
+!!$                endif
-                  info = -72
+!!$                aspk(nza) = val(i)
-                  return
+!!$              end if
-                endif
+!!$            end if
-                aspk(nza) = val(i)
+!!$          end do
-              end if
+!!$        end if
-            end if
+!!$      else
-          end do
+!!$        if ((nza >= nzl)) then 
-        end if
+!!$          do i=1, nz 
-      else
+!!$            nza = nza + 1 
-        if ((nza >= nzl)) then 
+!!$            if (nza>maxsz) then 
-          do i=1, nz 
+!!$              info = -73
-            nza = nza + 1 
+!!$              return
-            if (nza>maxsz) then 
+!!$            endif
-              info = -73
+!!$            aspk(nza) = val(i)
-              return
+!!$          end do
-            endif
+!!$        else
-            aspk(nza) = val(i)
+!!$          do i=1, nz 
-          end do
+!!$            ir = ia(i)
-        else
+!!$            ic = ja(i) 
-          do i=1, nz 
+!!$            if ((ir >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then 
-            ir = ia(i)
+!!$              nza = nza + 1 
-            ic = ja(i) 
+!!$              if (nza>maxsz) then 
-            if ((ir >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then 
+!!$                info = -74
-              nza = nza + 1 
+!!$                return
-              if (nza>maxsz) then 
+!!$              endif
-                info = -74
+!!$              aspk(nza) = val(i)
-                return
+!!$            end if
-              endif
+!!$          end do
-              aspk(nza) = val(i)
+!!$        end if
-            end if
+!!$      end if
-          end do
+!!$    end subroutine psb_inner_upd
-        end if
-      end if
-    end subroutine psb_inner_upd
 
     subroutine psb_inner_ins(nz,ia,ja,val,nza,ia1,ia2,aspk,maxsz,&
-         & imin,imax,jmin,jmax,info,gtl,ng)
+         & imin,imax,jmin,jmax,info,gtl)
       implicit none 
 
       integer, intent(in) :: nz, imin,imax,jmin,jmax,maxsz
@@ -229,15 +261,13 @@ contains
       real(psb_dpk_), intent(in) :: val(:)
       real(psb_dpk_), intent(inout) :: aspk(:)
       integer, intent(out) :: info
-      integer, intent(in), optional  :: ng,gtl(:)
+      integer, intent(in), optional  :: gtl(:)
-      integer :: i,ir,ic
+      integer :: i,ir,ic,ng
 
       info = 0
       if (present(gtl)) then 
-        if (.not.present(ng)) then 
+        ng = size(gtl) 
-          info = -1
+
-          return
-        endif
         do i=1, nz 
           ir = ia(i)
           ic = ja(i) 
@@ -277,12 +307,448 @@ contains
     end subroutine psb_inner_ins
 
 
+!!$#ifdef FIXED_NAG_SEGV
+!!$    subroutine d_coo_srch_upd(nz,ia,ja,val,a,&
+!!$         & imin,imax,jmin,jmax,info,gtl)
+!!$      
+!!$      use psb_const_mod
+!!$      use psb_realloc_mod
+!!$      use psb_string_mod
+!!$      use psb_serial_mod
+!!$      implicit none 
+!!$      
+!!$      class(psbn_d_coo_sparse_mat), intent(inout) :: a
+!!$      integer, intent(in) :: nz, imin,imax,jmin,jmax
+!!$      integer, intent(in) :: ia(:),ja(:)
+!!$      real(psb_dpk_), intent(in) :: val(:)
+!!$      integer, intent(out) :: info
+!!$      integer, intent(in), optional  :: gtl(:)
+!!$      integer  :: i,ir,ic, ilr, ilc, ip, &
+!!$           & i1,i2,nc,nnz,dupl,ng
+!!$      integer              :: debug_level, debug_unit
+!!$      character(len=20)    :: name='d_coo_srch_upd'
+!!$      
+!!$      info = 0
+!!$      debug_unit  = psb_get_debug_unit()
+!!$      debug_level = psb_get_debug_level()
+!!$
+!!$      dupl = a%get_dupl()
+!!$      
+!!$      if (.not.a%is_sorted()) then 
+!!$        info = -4
+!!$        return
+!!$      end if
+!!$      
+!!$      ilr = -1 
+!!$      ilc = -1 
+!!$      nnz = a%get_nzeros()
+!!$      
+!!$      
+!!$      if (present(gtl)) then
+!!$        ng = size(gtl)
+!!$        
+!!$        select case(dupl)
+!!$        case(psbn_dupl_ovwrt_,psbn_dupl_err_)
+!!$          ! Overwrite.
+!!$          ! Cannot test for error, should have been caught earlier.
+!!$          do i=1, nz
+!!$            ir = ia(i)
+!!$            ic = ja(i) 
+!!$            if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then 
+!!$              ir = gtl(ir)
+!!$              if ((ir > 0).and.(ir <= a%m)) then 
+!!$                ic = gtl(ic) 
+!!$                if (ir /= ilr) then 
+!!$                  i1 = psb_ibsrch(ir,nnz,a%ia)
+!!$                  i2 = i1
+!!$                  do 
+!!$                    if (i2+1 > nnz) exit
+!!$                    if (a%ia(i2+1) /= a%ia(i2)) exit
+!!$                    i2 = i2 + 1
+!!$                  end do
+!!$                  do 
+!!$                    if (i1-1 < 1) exit
+!!$                    if (a%ia(i1-1) /= a%ia(i1)) exit
+!!$                    i1 = i1 - 1
+!!$                  end do
+!!$                  ilr = ir
+!!$                else
+!!$                  i1 = 1
+!!$                  i2 = 1
+!!$                end if
+!!$                nc = i2-i1+1
+!!$                ip = psb_issrch(ic,nc,a%ja(i1:i2))
+!!$                if (ip>0) then 
+!!$                  a%val(i1+ip-1) = val(i)
+!!$                else
+!!$                  info = i 
+!!$                  return
+!!$                end if
+!!$              else
+!!$                if (debug_level >= psb_debug_serial_) &
+!!$                     & write(debug_unit,*) trim(name),&
+!!$                     & ': Discarding row that does not belong to us.'
+!!$              endif
+!!$            end if
+!!$          end do
+!!$        case(psbn_dupl_add_)
+!!$          ! Add
+!!$          do i=1, nz
+!!$            ir = ia(i)
+!!$            ic = ja(i) 
+!!$            if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then 
+!!$              ir = gtl(ir)
+!!$              ic = gtl(ic) 
+!!$              if ((ir > 0).and.(ir <= a%m)) then 
+!!$                
+!!$                if (ir /= ilr) then 
+!!$                  i1 = psb_ibsrch(ir,nnz,a%ia)
+!!$                  i2 = i1
+!!$                  do 
+!!$                    if (i2+1 > nnz) exit
+!!$                    if (a%ia(i2+1) /= a%ia(i2)) exit
+!!$                    i2 = i2 + 1
+!!$                  end do
+!!$                  do 
+!!$                    if (i1-1 < 1) exit
+!!$                    if (a%ia(i1-1) /= a%ia(i1)) exit
+!!$                    i1 = i1 - 1
+!!$                  end do
+!!$                  ilr = ir
+!!$                else
+!!$                  i1 = 1
+!!$                  i2 = 1
+!!$                end if
+!!$                nc = i2-i1+1
+!!$                ip = psb_issrch(ic,nc,a%ja(i1:i2))
+!!$                if (ip>0) then 
+!!$                  a%val(i1+ip-1) = a%val(i1+ip-1) + val(i)
+!!$                else
+!!$                  info = i 
+!!$                  return
+!!$                end if
+!!$              else
+!!$                if (debug_level >= psb_debug_serial_) &
+!!$                     & write(debug_unit,*) trim(name),&
+!!$                     & ': Discarding row that does not belong to us.'              
+!!$              end if
+!!$            end if
+!!$          end do
+!!$          
+!!$        case default
+!!$          info = -3
+!!$          if (debug_level >= psb_debug_serial_) &
+!!$               & write(debug_unit,*) trim(name),&
+!!$               & ': Duplicate handling: ',dupl
+!!$        end select
+!!$        
+!!$      else
+!!$        
+!!$        select case(dupl)
+!!$        case(psbn_dupl_ovwrt_,psbn_dupl_err_)
+!!$          ! Overwrite.
+!!$          ! Cannot test for error, should have been caught earlier.
+!!$          do i=1, nz
+!!$            ir = ia(i)
+!!$            ic = ja(i) 
+!!$            if ((ir > 0).and.(ir <= a%m)) then 
+!!$
+!!$              if (ir /= ilr) then 
+!!$                i1 = psb_ibsrch(ir,nnz,a%ia)
+!!$                i2 = i1
+!!$                do 
+!!$                  if (i2+1 > nnz) exit
+!!$                  if (a%ia(i2+1) /= a%ia(i2)) exit
+!!$                  i2 = i2 + 1
+!!$                end do
+!!$                do 
+!!$                  if (i1-1 < 1) exit
+!!$                  if (a%ia(i1-1) /= a%ia(i1)) exit
+!!$                  i1 = i1 - 1
+!!$                end do
+!!$                ilr = ir
+!!$              else
+!!$                i1 = 1
+!!$                i2 = 1
+!!$              end if
+!!$              nc = i2-i1+1
+!!$              ip = psb_issrch(ic,nc,a%ja(i1:i2))
+!!$              if (ip>0) then 
+!!$                a%val(i1+ip-1) = val(i)
+!!$              else
+!!$                info = i 
+!!$                return
+!!$              end if
+!!$            end if
+!!$          end do
+!!$
+!!$        case(psbn_dupl_add_)
+!!$          ! Add
+!!$          do i=1, nz
+!!$            ir = ia(i)
+!!$            ic = ja(i) 
+!!$            if ((ir > 0).and.(ir <= a%m)) then 
+!!$
+!!$              if (ir /= ilr) then 
+!!$                i1 = psb_ibsrch(ir,nnz,a%ia)
+!!$                i2 = i1
+!!$                do 
+!!$                  if (i2+1 > nnz) exit
+!!$                  if (a%ia(i2+1) /= a%ia(i2)) exit
+!!$                  i2 = i2 + 1
+!!$                end do
+!!$                do 
+!!$                  if (i1-1 < 1) exit
+!!$                  if (a%ia(i1-1) /= a%ia(i1)) exit
+!!$                  i1 = i1 - 1
+!!$                end do
+!!$                ilr = ir
+!!$              else
+!!$                i1 = 1
+!!$                i2 = 1
+!!$              end if
+!!$              nc = i2-i1+1
+!!$              ip = psb_issrch(ic,nc,a%ja(i1:i2))
+!!$              if (ip>0) then 
+!!$                a%val(i1+ip-1) = a%val(i1+ip-1) + val(i)
+!!$              else
+!!$                info = i 
+!!$                return
+!!$              end if
+!!$            end if
+!!$          end do
+!!$
+!!$        case default
+!!$          info = -3
+!!$          if (debug_level >= psb_debug_serial_) &
+!!$               & write(debug_unit,*) trim(name),&
+!!$               & ': Duplicate handling: ',dupl
+!!$        end select
+!!$
+!!$      end if
+!!$
+!!$    end subroutine d_coo_srch_upd
+!!$
+!!$#else
+    subroutine d_coo_srch_upd(nz,ia,ja,val,&
+         & aia,aja,aval,dupl,nza,nra,&
+         & info,gtl)
+
+      use psb_error_mod
+      use psb_sort_mod
+      implicit none 
+
+      integer, intent(inout) :: aia(:),aja(:)
+      real(psb_dpk_), intent(inout) :: aval(:)
+      integer, intent(in) :: nz, dupl,nza, nra
+      integer, intent(in) :: ia(:),ja(:)
+      real(psb_dpk_), intent(in) :: val(:)
+      integer, intent(out) :: info
+      integer, intent(in), optional  :: gtl(:)
+      integer  :: i,ir,ic, ilr, ilc, ip, &
+           & i1,i2,nc,ng
+      integer              :: debug_level, debug_unit
+      character(len=20)    :: name='d_coo_srch_upd'
+
+      info = 0
+      debug_unit  = psb_get_debug_unit()
+      debug_level = psb_get_debug_level()
+
+
+      ilr = -1 
+      ilc = -1 
+
+
+      if (present(gtl)) then
+        ng = size(gtl)
+        
+        select case(dupl)
+
+        case(psbn_dupl_ovwrt_,psbn_dupl_err_)
+          ! Overwrite.
+          ! Cannot test for error, should have been caught earlier.
+          do i=1, nz
+            ir = ia(i)
+            ic = ja(i) 
+            if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then 
+              ir = gtl(ir)
+              if ((ir > 0).and.(ir <= nra)) then 
+                ic = gtl(ic) 
+                if (ir /= ilr) then 
+                  i1 = psb_ibsrch(ir,nza,aia)
+                  i2 = i1
+                  do 
+                    if (i2+1 > nza) exit
+                    if (aia(i2+1) /= aia(i2)) exit
+                    i2 = i2 + 1
+                  end do
+                  do 
+                    if (i1-1 < 1) exit
+                    if (aia(i1-1) /= aia(i1)) exit
+                    i1 = i1 - 1
+                  end do
+                  ilr = ir
+                else
+                  i1 = 1
+                  i2 = 1
+                end if
+                nc = i2-i1+1
+                ip = psb_issrch(ic,nc,aja(i1:i2))
+                if (ip>0) then 
+                  aval(i1+ip-1) = val(i)
+                else
+                  info = i 
+                  return
+                end if
+              else
+                if (debug_level >= psb_debug_serial_) &
+                     & write(debug_unit,*) trim(name),&
+                     & ': Discarding row that does not belong to us.'
+              endif
+            end if
+          end do
+        case(psbn_dupl_add_)
+          ! Add
+          do i=1, nz
+            ir = ia(i)
+            ic = ja(i) 
+            if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then 
+              ir = gtl(ir)
+              ic = gtl(ic) 
+              if ((ir > 0).and.(ir <= nra)) then 
+                
+                if (ir /= ilr) then 
+                  i1 = psb_ibsrch(ir,nza,aia)
+                  i2 = i1
+                  do 
+                    if (i2+1 > nza) exit
+                    if (aia(i2+1) /= aia(i2)) exit
+                    i2 = i2 + 1
+                  end do
+                  do 
+                    if (i1-1 < 1) exit
+                    if (aia(i1-1) /= aia(i1)) exit
+                    i1 = i1 - 1
+                  end do
+                  ilr = ir
+                else
+                  i1 = 1
+                  i2 = 1
+                end if
+                nc = i2-i1+1
+                ip = psb_issrch(ic,nc,aja(i1:i2))
+                if (ip>0) then 
+                  aval(i1+ip-1) = aval(i1+ip-1) + val(i)
+                else
+                  info = i 
+                  return
+                end if
+              else
+                if (debug_level >= psb_debug_serial_) &
+                     & write(debug_unit,*) trim(name),&
+                     & ': Discarding row that does not belong to us.'              
+              end if
+            end if
+          end do
           
+        case default
+          info = -3
+          if (debug_level >= psb_debug_serial_) &
+               & write(debug_unit,*) trim(name),&
+               & ': Duplicate handling: ',dupl
+        end select
+        
+      else
+        
+        select case(dupl)
+        case(psbn_dupl_ovwrt_,psbn_dupl_err_)
+          ! Overwrite.
+          ! Cannot test for error, should have been caught earlier.
+          do i=1, nz
+            ir = ia(i)
+            ic = ja(i) 
+            if ((ir > 0).and.(ir <= nra)) then 
+
+              if (ir /= ilr) then 
+                i1 = psb_ibsrch(ir,nza,aia)
+                i2 = i1
+                do 
+                  if (i2+1 > nza) exit
+                  if (aia(i2+1) /= aia(i2)) exit
+                  i2 = i2 + 1
+                end do
+                do 
+                  if (i1-1 < 1) exit
+                  if (aia(i1-1) /= aia(i1)) exit
+                  i1 = i1 - 1
+                end do
+                ilr = ir
+              else
+                i1 = 1
+                i2 = 1
+              end if
+              nc = i2-i1+1
+              ip = psb_issrch(ic,nc,aja(i1:i2))
+              if (ip>0) then 
+                aval(i1+ip-1) = val(i)
+              else
+                info = i 
+                return
+              end if
+            end if
+          end do
+
+        case(psbn_dupl_add_)
+          ! Add
+          do i=1, nz
+            ir = ia(i)
+            ic = ja(i) 
+            if ((ir > 0).and.(ir <= nra)) then 
+
+              if (ir /= ilr) then 
+                i1 = psb_ibsrch(ir,nza,aia)
+                i2 = i1
+                do 
+                  if (i2+1 > nza) exit
+                  if (aia(i2+1) /= aia(i2)) exit
+                  i2 = i2 + 1
+                end do
+                do 
+                  if (i1-1 < 1) exit
+                  if (aia(i1-1) /= aia(i1)) exit
+                  i1 = i1 - 1
+                end do
+                ilr = ir
+              else
+                i1 = 1
+                i2 = 1
+              end if
+              nc = i2-i1+1
+              ip = psb_issrch(ic,nc,aja(i1:i2))
+              if (ip>0) then 
+                aval(i1+ip-1) = aval(i1+ip-1) + val(i)
+              else
+                info = i 
+                return
+              end if
+            end if
+          end do
+
+        case default
+          info = -3
+          if (debug_level >= psb_debug_serial_) &
+               & write(debug_unit,*) trim(name),&
+               & ': Duplicate handling: ',dupl
+        end select
+
+      end if
+
+    end subroutine d_coo_srch_upd
+!!$#endif
   end subroutine d_coo_csins
 
 
   subroutine d_coo_csmv(alpha,a,x,beta,y,info,trans) 
-    use psb_const_mod
     use psb_error_mod
     class(psbn_d_coo_sparse_mat), intent(in) :: a
     real(psb_dpk_), intent(in)          :: alpha, beta, x(:)
@@ -433,7 +899,6 @@ contains
   end subroutine d_coo_csmv
 
   subroutine d_coo_csmm(alpha,a,x,beta,y,info,trans) 
-    use psb_const_mod
     use psb_error_mod
     class(psbn_d_coo_sparse_mat), intent(in) :: a
     real(psb_dpk_), intent(in)          :: alpha, beta, x(:,:)
@@ -593,7 +1058,6 @@ contains
   end subroutine d_coo_csmm
 
   subroutine d_coo_cssv(alpha,a,x,beta,y,info,trans) 
-    use psb_const_mod
     use psb_error_mod
     class(psbn_d_coo_sparse_mat), intent(in) :: a
     real(psb_dpk_), intent(in)          :: alpha, beta, x(:)
@@ -689,7 +1153,7 @@ contains
   contains 
 
     subroutine inner_coosv(tra,a,x,y,info) 
-      use psb_const_mod
+
       logical, intent(in)                 :: tra  
       class(psbn_d_coo_sparse_mat), intent(in) :: a
       real(psb_dpk_), intent(in)          :: x(:)
@@ -850,7 +1314,6 @@ contains
 
 
   subroutine d_coo_cssm(alpha,a,x,beta,y,info,trans) 
-    use psb_const_mod
     use psb_error_mod
     class(psbn_d_coo_sparse_mat), intent(in) :: a
     real(psb_dpk_), intent(in)          :: alpha, beta, x(:,:)
@@ -949,7 +1412,6 @@ contains
   contains 
 
     subroutine inner_coosm(tra,a,x,y,info) 
-      use psb_const_mod
       logical, intent(in)                 :: tra  
       class(psbn_d_coo_sparse_mat), intent(in) :: a
       real(psb_dpk_), intent(in)          :: x(:,:)
@@ -1113,7 +1575,5 @@ contains
 
   end subroutine d_coo_cssm
   
-
-  
 end module psbn_d_coo_sparse_mat_mod
 

base/newserial/psbn_csr_mat.f03

@@ -3,22 +3,58 @@ module psbn_d_csr_sparse_mat_mod
   use psbn_d_base_mat_mod
 
   type, extends(psbn_d_base_sparse_mat) :: psbn_d_csr_sparse_mat
+
     logical              :: sorted
     integer, allocatable :: irp(:), ja(:)
-    real(kind(1.d0)), allocatable :: val(:)
+    real(psb_dpk_), allocatable :: val(:)
+
   contains
-    procedure, pass(a)  :: d_csr_get_nzeros
+    procedure, pass(a)  :: get_nzeros => d_csr_get_nzeros
     procedure, pass(a)  :: d_base_csmm => d_csr_csmm
     procedure, pass(a)  :: d_base_csmv => d_csr_csmv
-    generic, public     :: base_get_nzeros => d_csr_get_nzeros
     procedure, pass(a)  :: d_base_cssm => d_csr_cssm
     procedure, pass(a)  :: d_base_cssv => d_csr_cssv
+    procedure, pass(a)  :: reallocate_nz => d_csr_reallocate_nz
+    procedure, pass(a)  :: csins => d_csr_csins
 
 
   end type psbn_d_csr_sparse_mat
 
 contains 
 
+  subroutine  d_csr_reallocate_nz(nz,a) 
+    use psb_error_mod
+    use psb_realloc_mod
+    integer, intent(in) :: nz
+    class(psbn_d_csr_sparse_mat), intent(inout) :: a
+    Integer :: err_act, info
+    character(len=20)  :: name='d_csr_reallocate_nz'
+    logical, parameter :: debug=.false.
+
+    call psb_erractionsave(err_act)
+
+    call psb_realloc(nx,a%ja,info)
+    if (info == 0) call psb_realloc(nx,a%val,info)
+    if (info == 0) call psb_realloc(max(nx,a%m+1,a%n+1),a%irp,info)
+    if (info /= 0) then 
+      call psb_errpush(4000,name)
+      goto 9999
+    end if
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+
+  end subroutine d_csr_reallocate_nz
+
   function d_csr_get_nzeros(a) result(res)
     class(psbn_d_csr_sparse_mat), intent(in) :: a
     integer :: res
@@ -26,8 +62,516 @@ contains
   end function d_csr_get_nzeros
   
 
-  subroutine d_csr_csmv(alpha,a,x,beta,y,info,trans) 
+  subroutine d_csr_csins(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl) 
     use psb_const_mod
+    use psb_error_mod
+    class(psbn_d_csr_sparse_mat), intent(inout) :: a
+    real(psb_dpk_), intent(in)      :: val(:)
+    integer, intent(in)             :: nz, ia(:), ja(:), imin,imax,jmin,jmax
+    integer, intent(out)            :: info
+    integer, intent(in), optional   :: gtl(:)
+
+
+    Integer            :: err_act
+    character(len=20)  :: name='d_csr_csins'
+    logical, parameter :: debug=.false.
+    integer            :: nza, i,j,k, nzl, isza, int_err(5)
+
+    call psb_erractionsave(err_act)
+    info = 0
+
+    if (nz <= 0) then 
+      info = 10
+      int_err(1)=1
+      call psb_errpush(info,name,i_err=int_err)
+      goto 9999
+    end if
+    if (size(ia) < nz) then 
+      info = 35
+      int_err(1)=2
+      call psb_errpush(info,name,i_err=int_err)
+      goto 9999
+    end if
+
+    if (size(ja) < nz) then 
+      info = 35
+      int_err(1)=3
+      call psb_errpush(info,name,i_err=int_err)
+      goto 9999
+    end if
+    if (size(val) < nz) then 
+      info = 35
+      int_err(1)=4
+      call psb_errpush(info,name,i_err=int_err)
+      goto 9999
+    end if
+
+    if (nz == 0) return
+
+
+    nza  = a%get_nzeros()
+
+    if (a%is_bld()) then 
+      ! Build phase should only ever be in COO
+      info = 1121
+
+    else  if (a%is_upd()) then 
+      if (a%is_sorted()) then 
+
+
+!!$#ifdef FIXED_NAG_SEGV
+!!$        call  d_csr_srch_upd(nz,ia,ja,val,a,&
+!!$             & imin,imax,jmin,jmax,info,gtl)
+!!$#else 
+        call  d_csr_srch_upd(nz,ia,ja,val,&
+             & a%irp,a%ja,a%val,&
+             & a%get_dupl(),a%get_nzeros(),a%get_nrows(),&
+             & info,gtl)
+!!$#endif
+
+      else
+        info = 1121
+      end if
+
+    else 
+      ! State is wrong.
+      info = 1121
+    end if
+    if (info /= 0) then
+      call psb_errpush(info,name)
+      goto 9999
+    end if
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+
+
+  contains
+
+
+!!$#ifdef FIXED_NAG_SEGV
+!!$    subroutine d_csr_srch_upd(nz,ia,ja,val,a,&
+!!$         & imin,imax,jmin,jmax,info,gtl)
+!!$      
+!!$      use psb_const_mod
+!!$      use psb_realloc_mod
+!!$      use psb_string_mod
+!!$      use psb_serial_mod
+!!$      implicit none 
+!!$      
+!!$      class(psbn_d_csr_sparse_mat), intent(inout) :: a
+!!$      integer, intent(in) :: nz, imin,imax,jmin,jmax
+!!$      integer, intent(in) :: ia(:),ja(:)
+!!$      real(psb_dpk_), intent(in) :: val(:)
+!!$      integer, intent(out) :: info
+!!$      integer, intent(in), optional  :: gtl(:)
+!!$      integer  :: i,ir,ic, ilr, ilc, ip, &
+!!$           & i1,i2,nc,nnz,dupl,ng
+!!$      integer              :: debug_level, debug_unit
+!!$      character(len=20)    :: name='d_csr_srch_upd'
+!!$      
+!!$      info = 0
+!!$      debug_unit  = psb_get_debug_unit()
+!!$      debug_level = psb_get_debug_level()
+!!$
+!!$      dupl = a%get_dupl()
+!!$      
+!!$      if (.not.a%is_sorted()) then 
+!!$        info = -4
+!!$        return
+!!$      end if
+!!$      
+!!$      ilr = -1 
+!!$      ilc = -1 
+!!$      nnz = a%get_nzeros()
+!!$      
+!!$      
+!!$      if (present(gtl)) then
+!!$        ng = size(gtl)
+!!$        
+!!$        select case(dupl)
+!!$        case(psbn_dupl_ovwrt_,psbn_dupl_err_)
+!!$          ! Overwrite.
+!!$          ! Cannot test for error, should have been caught earlier.
+!!$          do i=1, nz
+!!$            ir = ia(i)
+!!$            ic = ja(i) 
+!!$            if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then 
+!!$              ir = gtl(ir)
+!!$              if ((ir > 0).and.(ir <= a%m)) then 
+!!$                ic = gtl(ic) 
+!!$                if (ir /= ilr) then 
+!!$                  i1 = psb_ibsrch(ir,nnz,a%ia)
+!!$                  i2 = i1
+!!$                  do 
+!!$                    if (i2+1 > nnz) exit
+!!$                    if (a%ia(i2+1) /= a%ia(i2)) exit
+!!$                    i2 = i2 + 1
+!!$                  end do
+!!$                  do 
+!!$                    if (i1-1 < 1) exit
+!!$                    if (a%ia(i1-1) /= a%ia(i1)) exit
+!!$                    i1 = i1 - 1
+!!$                  end do
+!!$                  ilr = ir
+!!$                else
+!!$                  i1 = 1
+!!$                  i2 = 1
+!!$                end if
+!!$                nc = i2-i1+1
+!!$                ip = psb_issrch(ic,nc,a%ja(i1:i2))
+!!$                if (ip>0) then 
+!!$                  a%val(i1+ip-1) = val(i)
+!!$                else
+!!$                  info = i 
+!!$                  return
+!!$                end if
+!!$              else
+!!$                if (debug_level >= psb_debug_serial_) &
+!!$                     & write(debug_unit,*) trim(name),&
+!!$                     & ': Discarding row that does not belong to us.'
+!!$              endif
+!!$            end if
+!!$          end do
+!!$        case(psbn_dupl_add_)
+!!$          ! Add
+!!$          do i=1, nz
+!!$            ir = ia(i)
+!!$            ic = ja(i) 
+!!$            if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then 
+!!$              ir = gtl(ir)
+!!$              ic = gtl(ic) 
+!!$              if ((ir > 0).and.(ir <= a%m)) then 
+!!$                
+!!$                if (ir /= ilr) then 
+!!$                  i1 = psb_ibsrch(ir,nnz,a%ia)
+!!$                  i2 = i1
+!!$                  do 
+!!$                    if (i2+1 > nnz) exit
+!!$                    if (a%ia(i2+1) /= a%ia(i2)) exit
+!!$                    i2 = i2 + 1
+!!$                  end do
+!!$                  do 
+!!$                    if (i1-1 < 1) exit
+!!$                    if (a%ia(i1-1) /= a%ia(i1)) exit
+!!$                    i1 = i1 - 1
+!!$                  end do
+!!$                  ilr = ir
+!!$                else
+!!$                  i1 = 1
+!!$                  i2 = 1
+!!$                end if
+!!$                nc = i2-i1+1
+!!$                ip = psb_issrch(ic,nc,a%ja(i1:i2))
+!!$                if (ip>0) then 
+!!$                  a%val(i1+ip-1) = a%val(i1+ip-1) + val(i)
+!!$                else
+!!$                  info = i 
+!!$                  return
+!!$                end if
+!!$              else
+!!$                if (debug_level >= psb_debug_serial_) &
+!!$                     & write(debug_unit,*) trim(name),&
+!!$                     & ': Discarding row that does not belong to us.'              
+!!$              end if
+!!$            end if
+!!$          end do
+!!$          
+!!$        case default
+!!$          info = -3
+!!$          if (debug_level >= psb_debug_serial_) &
+!!$               & write(debug_unit,*) trim(name),&
+!!$               & ': Duplicate handling: ',dupl
+!!$        end select
+!!$        
+!!$      else
+!!$        
+!!$        select case(dupl)
+!!$        case(psbn_dupl_ovwrt_,psbn_dupl_err_)
+!!$          ! Overwrite.
+!!$          ! Cannot test for error, should have been caught earlier.
+!!$          do i=1, nz
+!!$            ir = ia(i)
+!!$            ic = ja(i) 
+!!$            if ((ir > 0).and.(ir <= a%m)) then 
+!!$
+!!$              if (ir /= ilr) then 
+!!$                i1 = psb_ibsrch(ir,nnz,a%ia)
+!!$                i2 = i1
+!!$                do 
+!!$                  if (i2+1 > nnz) exit
+!!$                  if (a%ia(i2+1) /= a%ia(i2)) exit
+!!$                  i2 = i2 + 1
+!!$                end do
+!!$                do 
+!!$                  if (i1-1 < 1) exit
+!!$                  if (a%ia(i1-1) /= a%ia(i1)) exit
+!!$                  i1 = i1 - 1
+!!$                end do
+!!$                ilr = ir
+!!$              else
+!!$                i1 = 1
+!!$                i2 = 1
+!!$              end if
+!!$              nc = i2-i1+1
+!!$              ip = psb_issrch(ic,nc,a%ja(i1:i2))
+!!$              if (ip>0) then 
+!!$                a%val(i1+ip-1) = val(i)
+!!$              else
+!!$                info = i 
+!!$                return
+!!$              end if
+!!$            end if
+!!$          end do
+!!$
+!!$        case(psbn_dupl_add_)
+!!$          ! Add
+!!$          do i=1, nz
+!!$            ir = ia(i)
+!!$            ic = ja(i) 
+!!$            if ((ir > 0).and.(ir <= a%m)) then 
+!!$
+!!$              if (ir /= ilr) then 
+!!$                i1 = psb_ibsrch(ir,nnz,a%ia)
+!!$                i2 = i1
+!!$                do 
+!!$                  if (i2+1 > nnz) exit
+!!$                  if (a%ia(i2+1) /= a%ia(i2)) exit
+!!$                  i2 = i2 + 1
+!!$                end do
+!!$                do 
+!!$                  if (i1-1 < 1) exit
+!!$                  if (a%ia(i1-1) /= a%ia(i1)) exit
+!!$                  i1 = i1 - 1
+!!$                end do
+!!$                ilr = ir
+!!$              else
+!!$                i1 = 1
+!!$                i2 = 1
+!!$              end if
+!!$              nc = i2-i1+1
+!!$              ip = psb_issrch(ic,nc,a%ja(i1:i2))
+!!$              if (ip>0) then 
+!!$                a%val(i1+ip-1) = a%val(i1+ip-1) + val(i)
+!!$              else
+!!$                info = i 
+!!$                return
+!!$              end if
+!!$            end if
+!!$          end do
+!!$
+!!$        case default
+!!$          info = -3
+!!$          if (debug_level >= psb_debug_serial_) &
+!!$               & write(debug_unit,*) trim(name),&
+!!$               & ': Duplicate handling: ',dupl
+!!$        end select
+!!$
+!!$      end if
+!!$
+!!$    end subroutine d_csr_srch_upd
+!!$
+!!$#else
+    subroutine d_csr_srch_upd(nz,ia,ja,val,&
+         & airp,aja,aval,dupl,nza,nra,&
+         & info,gtl)
+
+      use psb_error_mod
+      use psb_sort_mod
+      implicit none 
+
+      integer, intent(inout) :: airp(:),aja(:)
+      real(psb_dpk_), intent(inout) :: aval(:)
+      integer, intent(in) :: nz, dupl,nza, nra
+      integer, intent(in) :: ia(:),ja(:)
+      real(psb_dpk_), intent(in) :: val(:)
+      integer, intent(out) :: info
+      integer, intent(in), optional  :: gtl(:)
+
+      integer  :: debug_level, debug_unit
+      character(len=20)    :: name='d_csr_srch_upd'
+      integer  :: i,ir,ic, ilr, ilc, ip, &
+           & i1,i2,nc,lb,ub,m, ng
+
+      info = 0
+      debug_unit  = psb_get_debug_unit()
+      debug_level = psb_get_debug_level()
+
+
+      if (present(gtl)) then 
+        ng = size(gtl)
+
+        select case(dupl)
+        case(psbn_dupl_ovwrt_,psbn_dupl_err_)
+          ! Overwrite.
+          ! Cannot test for error, should have been caught earlier.
+
+          ilr = -1 
+          ilc = -1 
+          do i=1, nz
+            ir = ia(i)
+            ic = ja(i) 
+            if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then 
+              ir = gtl(ir)
+              ic = gtl(ic)
+              if ((ir > 0).and.(ir <= a%m)) then 
+                i1 = airp(ir)
+                i2 = airp(ir+1)
+                nc=i2-i1
+
+                ip = psb_ibsrch(ic,nc,aja(i1:i2-1))    
+                if (ip>0) then 
+                  aval(i1+ip-1) = val(i)
+                else
+                  if (debug_level >= psb_debug_serial_) &
+                       & write(debug_unit,*) trim(name),&
+                       & ': Was searching ',ic,' in: ',i1,i2,&
+                       & ' : ',aja(i1:i2-1)
+                  info = i
+                  return
+                end if
+
+              else
+
+                if (debug_level >= psb_debug_serial_) &
+                     & write(debug_unit,*) trim(name),&
+                     & ': Discarding row that does not belong to us.'
+              end if
+            end if
+          end do
+
+        case(psbn_dupl_add_)
+          ! Add
+          ilr = -1 
+          ilc = -1 
+          do i=1, nz
+            ir = ia(i)
+            ic = ja(i) 
+            if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then 
+              ir = gtl(ir)
+              ic = gtl(ic)
+              if ((ir > 0).and.(ir <= a%m)) then 
+                i1 = airp(ir)
+                i2 = airp(ir+1)
+                nc = i2-i1
+                ip = psb_ibsrch(ic,nc,aja(i1:i2-1))
+                if (ip>0) then 
+                  aval(i1+ip-1) = aval(i1+ip-1) + val(i)
+                else
+                  if (debug_level >= psb_debug_serial_) &
+                       & write(debug_unit,*) trim(name),&
+                       & ': Was searching ',ic,' in: ',i1,i2,&
+                       & ' : ',aja(i1:i2-1)
+                  info = i
+                  return
+                end if
+              else
+                if (debug_level >= psb_debug_serial_) &
+                     & write(debug_unit,*) trim(name),&
+                     & ': Discarding row that does not belong to us.'
+              end if
+
+            end if
+          end do
+
+        case default
+          info = -3
+          if (debug_level >= psb_debug_serial_) &
+               & write(debug_unit,*) trim(name),&
+               & ': Duplicate handling: ',dupl
+        end select
+
+      else
+
+        select case(dupl)
+        case(psbn_dupl_ovwrt_,psbn_dupl_err_)
+          ! Overwrite.
+          ! Cannot test for error, should have been caught earlier.
+
+          ilr = -1 
+          ilc = -1 
+          do i=1, nz
+            ir = ia(i)
+            ic = ja(i) 
+
+            if ((ir > 0).and.(ir <= a%m)) then 
+
+              i1 = airp(ir)
+              i2 = airp(ir+1)
+              nc=i2-i1
+
+              ip = psb_ibsrch(ic,nc,aja(i1:i2-1))    
+              if (ip>0) then 
+                aval(i1+ip-1) = val(i)
+              else
+                if (debug_level >= psb_debug_serial_) &
+                     & write(debug_unit,*) trim(name),&
+                     & ': Was searching ',ic,' in: ',i1,i2,&
+                     & ' : ',aja(i1:i2-1)
+                info = i
+                return
+              end if
+
+            else
+              if (debug_level >= psb_debug_serial_) &
+                   & write(debug_unit,*) trim(name),&
+                   & ': Discarding row that does not belong to us.'
+            end if
+
+          end do
+
+        case(psbn_dupl_add_)
+          ! Add
+          ilr = -1 
+          ilc = -1 
+          do i=1, nz
+            ir = ia(i)
+            ic = ja(i) 
+            if ((ir > 0).and.(ir <= a%m)) then 
+              i1 = airp(ir)
+              i2 = airp(ir+1)
+              nc = i2-i1
+              ip = psb_ibsrch(ic,nc,aja(i1:i2-1))
+              if (ip>0) then 
+                aval(i1+ip-1) = aval(i1+ip-1) + val(i)
+              else
+                info = i
+                return
+              end if
+            else
+              if (debug_level >= psb_debug_serial_) &
+                   & write(debug_unit,*) trim(name),&
+                   & ': Discarding row that does not belong to us.'
+            end if
+          end do
+
+        case default
+          info = -3
+          if (debug_level >= psb_debug_serial_) &
+               & write(debug_unit,*) trim(name),&
+               & ': Duplicate handling: ',dupl
+        end select
+        
+      end if
+    end subroutine d_csr_srch_upd
+!!$#endif
+  end subroutine d_csr_csins
+
+
+
+
+
+
+  subroutine d_csr_csmv(alpha,a,x,beta,y,info,trans) 
     use psb_error_mod
     class(psbn_d_csr_sparse_mat), intent(in) :: a
     real(psb_dpk_), intent(in)          :: alpha, beta, x(:)
@@ -290,7 +834,6 @@ contains
   end subroutine d_csr_csmv
 
   subroutine d_csr_csmm(alpha,a,x,beta,y,info,trans) 
-    use psb_const_mod
     use psb_error_mod
     class(psbn_d_csr_sparse_mat), intent(in) :: a
     real(psb_dpk_), intent(in)          :: alpha, beta, x(:,:)
@@ -558,7 +1101,6 @@ contains
 
 
   subroutine d_csr_cssv(alpha,a,x,beta,y,info,trans) 
-    use psb_const_mod
     use psb_error_mod
     class(psbn_d_csr_sparse_mat), intent(in) :: a
     real(psb_dpk_), intent(in)          :: alpha, beta, x(:)
@@ -644,7 +1186,6 @@ contains
   contains 
 
     subroutine inner_csrsv(tra,a,x,y) 
-      use psb_const_mod
       logical, intent(in)                 :: tra  
       class(psbn_d_csr_sparse_mat), intent(in) :: a
       real(psb_dpk_), intent(in)          :: x(:)
@@ -750,7 +1291,6 @@ contains
 
 
   subroutine d_csr_cssm(alpha,a,x,beta,y,info,trans) 
-    use psb_const_mod
     use psb_error_mod
     class(psbn_d_csr_sparse_mat), intent(in) :: a
     real(psb_dpk_), intent(in)          :: alpha, beta, x(:,:)
@@ -849,7 +1389,6 @@ contains
   contains 
 
     subroutine inner_csrsm(tra,a,x,y,info) 
-      use psb_const_mod
       logical, intent(in)                 :: tra  
       class(psbn_d_csr_sparse_mat), intent(in) :: a
       real(psb_dpk_), intent(in)          :: x(:,:)

base/newserial/psbn_d_base_mat_mod.f03

@@ -2,6 +2,7 @@
 module psbn_d_base_mat_mod
 
   use psbn_base_mat_mod
+
   type, extends(psbn_base_sparse_mat) :: psbn_d_base_sparse_mat
   contains
     procedure, pass(a) :: d_base_csmv
@@ -10,23 +11,23 @@ module psbn_d_base_mat_mod
     procedure, pass(a) :: d_base_cssv
     procedure, pass(a) :: d_base_cssm
     generic, public    :: psbn_cssm => d_base_cssm, d_base_cssv
-    procedure, pass(a) :: d_base_csins
+    procedure, pass(a) :: csins
-    generic, public    :: csins => d_base_csins
     
   end type psbn_d_base_sparse_mat
 
 contains 
 
-  subroutine d_base_csins(nz,val,ia,ja,a,info) 
+  subroutine csins(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl) 
-    use psb_const_mod
     use psb_error_mod
-    class(psbn_d_base_sparse_mat), intent(in) :: a
+    use psb_realloc_mod
-    real(psb_dpk_), intent(in)    :: val(:)
+    class(psbn_d_base_sparse_mat), intent(inout) :: a
-    integer, intent(in)           :: nz, ia(:), ja(:)
+    real(psb_dpk_), intent(in)      :: val(:)
+    integer, intent(in)             :: nz, ia(:), ja(:), imin,imax,jmin,jmax
     integer, intent(out)            :: info
+    integer, intent(in), optional   :: gtl(:)
 
     Integer :: err_act
-    character(len=20)  :: name='d_base_csins'
+    character(len=20)  :: name='csins'
     logical, parameter :: debug=.false.
 
     call psb_erractionsave(err_act)
@@ -41,10 +42,9 @@ contains
     end if
     return
 
-  end subroutine d_base_csins
+  end subroutine csins
 
   subroutine d_base_csmm(alpha,a,x,beta,y,info,trans) 
-    use psb_const_mod
     use psb_error_mod
     class(psbn_d_base_sparse_mat), intent(in) :: a
     real(psb_dpk_), intent(in)    :: alpha, beta, x(:,:)
@@ -60,7 +60,8 @@ contains
     ! This is the base version. If we get here
     ! it means the derived class is incomplete,
     ! so we throw an error.
-    call psb_errpush(700,name)
+    info = 700
+    call psb_errpush(info,name)
          
     if (err_act /= psb_act_ret_) then
       call psb_error()
@@ -85,7 +86,8 @@ contains
     ! This is the base version. If we get here
     ! it means the derived class is incomplete,
     ! so we throw an error.
-    call psb_errpush(700,name)
+    info = 700
+    call psb_errpush(info,name)
           
     if (err_act /= psb_act_ret_) then
       call psb_error()
@@ -108,11 +110,11 @@ contains
     logical, parameter :: debug=.false.
 
     call psb_erractionsave(err_act)
-    info = 700
     ! This is the base version. If we get here
     ! it means the derived class is incomplete,
     ! so we throw an error.
-    call psb_errpush(700,name)
+    info = 700
+    call psb_errpush(info,name)
           
     if (err_act /= psb_act_ret_) then
       call psb_error()
@@ -134,11 +136,11 @@ contains
     logical, parameter :: debug=.false.
 
     call psb_erractionsave(err_act)
-    info = 700
     ! This is the base version. If we get here
     ! it means the derived class is incomplete,
     ! so we throw an error.
-    call psb_errpush(700,name)
+    info = 700
+    call psb_errpush(info,name)
           
     if (err_act /= psb_act_ret_) then
       call psb_error()

config/ax_f90_module_extension.m4

@@ -13,7 +13,7 @@ dnl make it work even if ac_fc_ext="f"
 
 AC_DEFUN([AX_F90_MODULE_EXTENSION],[
 AC_CACHE_CHECK([fortran 90 modules extension],
-ax_f90_modext,
+ax_cv_f90_modext,
 [AC_LANG_PUSH(Fortran)
 i=0
 while test \( -f tmpdir_$i \) -o \( -d tmpdir_$i \) ; do
@@ -29,16 +29,16 @@ AC_COMPILE_IFELSE([
        end subroutine conftest_routine
       end module conftest_module
   ],
-  [ax_f90_modext=`ls | sed -n 's,conftest_module\.,,p'`
+  [ax_cv_f90_modext=`ls | sed -n 's,conftest_module\.,,p'`
-   if test x$ax_f90_modext = x ; then
+   if test x$ax_cv_f90_modext = x ; then
 dnl Some F90 compilers put module filename in uppercase letters
-     ax_f90_modext=`ls | sed -n 's,CONFTEST_MODULE\.,,p'`
+     ax_cv_f90_modext=`ls | sed -n 's,CONFTEST_MODULE\.,,p'`
-     if test x$ax_f90_modext = x ; then
+     if test x$ax_cv_f90_modext = x ; then
-       ax_f90_modext=unknown
+       ax_cv_f90_modext=unknown
      fi
    fi
   ],
-  [ax_f90_modext=unknown])
+  [ax_cv_f90_modext=unknown])
 cd ..
 rm -fr tmpdir_$i
 AC_LANG_POP(Fortran)

config/ax_f90_module_flag.m4

@@ -18,7 +18,7 @@ dnl make it work even if ac_fc_ext="f"
 
 AC_DEFUN([AX_F90_MODULE_FLAG],[
 AC_CACHE_CHECK([fortran 90 modules inclusion flag],
-ax_f90_modflag,
+ax_cv_f90_modflag,
 [AC_LANG_PUSH(Fortran)
 i=0
 while test \( -f tmpdir_$i \) -o \( -d tmpdir_$i \) ; do
@@ -35,9 +35,9 @@ AC_COMPILE_IFELSE([
       end module conftest_module
   ],[],[])
 cd ..
-ax_f90_modflag="not found"
+ax_cv_f90_modflag="not found"
 for ax_flag in "-I " "-M" "-p"; do
-  if test "$ax_f90_modflag" = "not found" ; then
+  if test "$ax_cv_f90_modflag" = "not found" ; then
     ax_save_FCFLAGS="$FCFLAGS"
     FCFLAGS="$ax_save_FCFLAGS ${ax_flag}tmpdir_$i"
     AC_COMPILE_IFELSE([
@@ -45,7 +45,7 @@ for ax_flag in "-I " "-M" "-p"; do
        use conftest_module
        call conftest_routine
        end program conftest_program
-      ],[ax_f90_modflag="$ax_flag"],[])
+      ],[ax_cv_f90_modflag="$ax_flag"],[])
     FCFLAGS="$ax_save_FCFLAGS"
   fi
 done

config/pac.m4

@@ -285,6 +285,38 @@ AC_MSG_RESULT(no)
 ]
 )
 
+
+dnl @synopsis PAC_ARG_SERIAL_MPI
+dnl
+dnl Test for --with-serial-mpi={yes|no}
+dnl 
+dnl 
+dnl
+dnl Example use:
+dnl
+dnl
+dnl @author Salvatore Filippone <salvatore.filippone@uniroma2.it>
+dnl
+AC_DEFUN([PAC_ARG_SERIAL_MPI],
+[
+AC_MSG_CHECKING([whether we want serial (fake) mpi])
+AC_ARG_WITH(serial-mpi,
+AC_HELP_STRING([--with-serial-mpi], 
+[Specify whether to enable a fake mpi library to run in serial mode.
+  --with-serial-mpi={yes|no}]),
+[
+pac_cv_serial_mpi="${withval}";
+],
+[pac_cv_serial_mpi="no";]
+)
+if test x"$pac_cv_serial_mpi" == x"yes" ; then
+   AC_MSG_RESULT([yes.])
+else
+   AC_MSG_RESULT([no.])
+fi
+]
+)
+
 dnl @synopsis PAC_FORTRAN_HAVE_PSBLAS( [ACTION-IF-FOUND [, ACTION-IF-NOT-FOUND]])
 dnl
 dnl Will try to compile and link a program using the PSBLAS library

configure.ac

@@ -109,14 +109,22 @@ fi
 if test "X$CC" == "X" ; then
 	AC_MSG_ERROR([Problem : No C compiler specified nor found!])
 fi
+
 ###############################################################################
 #			Suitable MPI compilers detection
 ###############################################################################
 # Note: Someday we will contemplate a fake MPI - configured version of PSBLAS
 ###############################################################################
+# First check whether the user required our serial (fake) mpi.
+PAC_ARG_SERIAL_MPI
 
 #Note : we miss the name of the Intel C compiler
-
+if test x"$pac_cv_serial_mpi" == x"yes" ; then
+   FAKEMPI="fakempi.o";
+   MPIFC="$FC";
+   MPIF77="$F77";
+   MPICC="$CC";
+else 
 AC_LANG([C])
 if test "X$MPICC" = "X" ; then
     # This is our MPICC compiler preference: it will override ACX_MPI's first try.
@@ -139,7 +147,7 @@ if test "X$MPIFC" = "X" ; then
 fi
 
 ACX_MPI([], [AC_MSG_ERROR([[Cannot find any suitable MPI implementation for Fortran]])])
-
+fi
 # We leave a default language for the next checks.
 dnl AC_LANG([Fortran 77])
 AC_LANG([C])
@@ -391,7 +399,7 @@ if test "X$FCOPT" == "X" ; then
 fi
 if test "X$psblas_cv_fc" == X"nag" ; then
    # Add needed options 
-   FCOPT="$FCOPT -mismatch -dcfuns"
+   FCOPT="$FCOPT -mismatch -dcfuns -f2003"
 fi
 FFLAGS="${FCOPT}"
 
@@ -424,9 +432,13 @@ else
 fi
 if test "X$psblas_cv_fc" == X"nag" ; then
    # Add needed options 
-     F90COPT="$F90COPT -mismatch -dcfuns"
+     F90COPT="$F90COPT -mismatch -dcfuns -f2003"
      EXTRA_OPT="-mismatch_all"
+     F03COPT="${F90COPT}"
+else
+     F03COPT="${F90COPT}"
 fi
+
 FCFLAGS="${F90COPT}"
 
 # COPT,FCOPT, F90COPT are aliases for FFLAGS,CFLAGS,FCFLAGS .
@@ -441,7 +453,6 @@ if test "X$psblas_cv_fc" == X"xlf" ; then
   #  - it is not said that mpxlf95 gets chosen by the configure script.
   F90="xlf95   -qsuffix=f=f90:cpp=F90"
  #  F90="xlf95"
-  MPF90="mpxlf95 -qsuffix=f=f90:cpp=F90"
 #    FC="xlf"
 
 # Note : this gives problems in base/serial/aux/isaperm.f
@@ -450,11 +461,17 @@ if test "X$psblas_cv_fc" == X"xlf" ; then
 # Note : this is the cure
      FC="xlf     -qsuffix=f=f:cpp=F"
      # Note : maybe we will want xlf -qsuffix=cpp=F
-
     F77="xlf"
-  MPF77="mpxlf95 -qfixed  -qsuffix=f=f:cpp=F"
      CC="xlc"
-   MPCC="mpcc"
+  if test x"$pac_cv_serial_mpi" == x"yes" ; then
+     MPF90="xlf95 -qsuffix=f=f90:cpp=F90"
+     MPF77="xlf95 -qfixed  -qsuffix=f=f:cpp=F"
+     MPCC="xlc"
+  else 
+     MPF90="mpxlf95 -qsuffix=f=f90:cpp=F90"
+     MPF77="mpxlf95 -qfixed  -qsuffix=f=f:cpp=F"
+     MPCC="mpcc"
+  fi
  #MPFCC="mpxlc"
   # Note : -qfixed should be not specified in the environment FFLAGS or things will break.
   #        This fact should be documented somewhere.
@@ -487,8 +504,8 @@ fi
 ##############################################################################
 AX_F90_MODULE_EXTENSION
 AX_F90_MODULE_FLAG
-MODEXT=".$ax_f90_modext"
+MODEXT=".$ax_cv_f90_modext"
-FMFLAG="${ax_f90_modflag%%[ ]*}"
+FMFLAG="${ax_cv_f90_modflag%%[ ]*}"
 FIFLAG=-I
 
 ###############################################################################
@@ -504,11 +521,14 @@ FIFLAG=-I
 
 ###############################################################################
 # Custom test : do we have a module or include for MPI Fortran interface?
-
+if test x"$pac_cv_serial_mpi" == x"yes" ; then
-PAC_FORTRAN_CHECK_HAVE_MPI_MOD(
+   FDEFINES="$psblas_cv_define_prepend-DSERIAL_MPI $psblas_cv_define_prepend-DMPI_MOD $FDEFINES";
-	[FDEFINES="$psblas_cv_define_prepend-DMPI_MOD $FDEFINES"],
+else 
+ PAC_FORTRAN_CHECK_HAVE_MPI_MOD(
+ 	[FDEFINES="$psblas_cv_define_prepend-DMPI_MOD $FDEFINES"],
 	[FDEFINES="$psblas_cv_define_prepend-DMPI_H $FDEFINES"]
-)
+ )
+fi
 
 # Custom test : do we have move_alloc ?
 PAC_FORTRAN_HAVE_MOVE_ALLOC(
@@ -577,6 +597,7 @@ ACX_LAPACK(
 ###############################################################################
 #		BLACS library presence checks
 ###############################################################################
+if test x"$pac_cv_serial_mpi" == x"no" ; then
 save_FC="$FC";
 save_CC="$CC";
 FC="$MPIFC";
@@ -584,6 +605,8 @@ CC="$MPICC";
 PAC_CHECK_BLACS
 FC="$save_FC";
 CC="$save_CC";
+fi
+
 PAC_MAKE_IS_GNUMAKE
 
 ###############################################################################
@@ -644,12 +667,14 @@ AC_SUBST(CDEFINES)
 AC_SUBST(BASELIBNAME)
 AC_SUBST(F90)
 AC_SUBST(F90COPT)
+AC_SUBST(F03COPT)
 AC_SUBST(MPF90)
 AC_SUBST(MPF77)
 AC_SUBST(MPCC)
 AC_SUBST(FCOPT)
 AC_SUBST(CCOPT)
 AC_SUBST(EXTRA_OPT)
+AC_SUBST(FAKEMPI)
 AC_SUBST(FIFLAG)
 AC_SUBST(FMFLAG)
 AC_SUBST(MODEXT)
@@ -710,9 +735,9 @@ $(.mod).o:
 %$(.mod): %.f03
 	$(F90) $(F90COPT) $(FINCLUDES)  -c $<
 %.o: %.F03
-	$(F90) $(F90COPT) $(FINCLUDES) $(FDEFINES)  -c $<
+	$(F90) $(F03COPT) $(FINCLUDES) $(FDEFINES)  -c $<
 %$(.mod): %.F03
-	$(F90) $(F90COPT) $(FINCLUDES) $(FDEFINES) -c $<'
+	$(F90) $(F03COPT) $(FINCLUDES) $(FDEFINES) -c $<'
 
 else
 
@@ -748,9 +773,9 @@ $(.mod).o:
 .f03.o:
 	$(F90) $(F90COPT) $(FINCLUDES) -c $<
 .F03.o:
-	$(F90) $(F90COPT) $(FINCLUDES) $(FDEFINES) -c $<
+	$(F90) $(F03COPT) $(FINCLUDES) $(FDEFINES) -c $<
 .F03$(.mod):
-	$(F90) $(F90COPT) $(FINCLUDES) $(FDEFINES) -c $<'
+	$(F90) $(F03COPT) $(FINCLUDES) $(FDEFINES) -c $<'
 
 fi
 AC_SUBST(PSBLASRULES)
@@ -785,6 +810,7 @@ dnl	FFLAGS                : ${FFLAGS}
 dnl	FCFLAGS               : ${FCFLAGS}
 	MODEXT                : ${MODEXT}
 	FMFLAG                : ${FMFLAG}
+   	F03COPT               : ${F03COPT}
    	F90COPT               : ${F90COPT}
    	FCOPT                 : ${FCOPT}
    	CCOPT                 : ${CCOPT}