From 2084fbd91d0e3b0e68b7ff8be4315cc0b0e6c9b6 Mon Sep 17 00:00:00 2001
From: Salvatore Filippone <salvatore.filippone@cranfield.ac.uk>
Date: Tue, 18 Mar 2014 15:24:22 +0000
Subject: [PATCH] psblas3:  base/modules/psb_c_csc_mat_mod.f90 
 base/modules/psb_c_csr_mat_mod.f90  base/modules/psb_c_mat_mod.f90 
 base/modules/psb_c_vect_mod.F90  base/modules/psb_const_mod.F90 
 base/modules/psb_d_csc_mat_mod.f90  base/modules/psb_d_csr_mat_mod.f90 
 base/modules/psb_d_mat_mod.f90  base/modules/psb_d_vect_mod.F90 
 base/modules/psb_i_base_vect_mod.f90  base/modules/psb_i_vect_mod.F90 
 base/modules/psb_s_csc_mat_mod.f90  base/modules/psb_s_csr_mat_mod.f90 
 base/modules/psb_s_mat_mod.f90  base/modules/psb_s_vect_mod.F90 
 base/modules/psb_z_csc_mat_mod.f90  base/modules/psb_z_csr_mat_mod.f90 
 base/modules/psb_z_mat_mod.f90  base/modules/psb_z_vect_mod.F90 
 base/serial/impl/psb_c_coo_impl.f90  base/serial/impl/psb_c_mat_impl.F90 
 base/serial/impl/psb_d_coo_impl.f90  base/serial/impl/psb_d_mat_impl.F90 
 base/serial/impl/psb_s_coo_impl.f90  base/serial/impl/psb_s_mat_impl.F90 
 base/serial/impl/psb_z_coo_impl.f90  base/serial/impl/psb_z_mat_impl.F90 
 base/serial/psb_cspspmm.f90  base/serial/psb_dspspmm.f90 
 base/serial/psb_sspspmm.f90  base/serial/psb_zspspmm.f90

New defaults for MOLD, inquiry and set functions.
New sparse-spare product version.
New is_by_rows/is_by_cols methods, as yet unused.
---
 base/modules/psb_c_csc_mat_mod.f90   | 10 ++++
 base/modules/psb_c_csr_mat_mod.f90   | 14 ++++-
 base/modules/psb_c_mat_mod.f90       | 79 ++++++++++++++++++++++++++-
 base/modules/psb_c_vect_mod.F90      | 69 ++++++++++++++++++++---
 base/modules/psb_const_mod.F90       |  5 +-
 base/modules/psb_d_csc_mat_mod.f90   | 10 ++++
 base/modules/psb_d_csr_mat_mod.f90   | 14 ++++-
 base/modules/psb_d_mat_mod.f90       | 79 ++++++++++++++++++++++++++-
 base/modules/psb_d_vect_mod.F90      | 69 ++++++++++++++++++++---
 base/modules/psb_i_base_vect_mod.f90 | 39 ++++++++++++-
 base/modules/psb_i_vect_mod.F90      | 82 +++++++++++++++++++++++-----
 base/modules/psb_s_csc_mat_mod.f90   | 10 ++++
 base/modules/psb_s_csr_mat_mod.f90   | 14 ++++-
 base/modules/psb_s_mat_mod.f90       | 79 ++++++++++++++++++++++++++-
 base/modules/psb_s_vect_mod.F90      | 69 ++++++++++++++++++++---
 base/modules/psb_z_csc_mat_mod.f90   | 10 ++++
 base/modules/psb_z_csr_mat_mod.f90   | 14 ++++-
 base/modules/psb_z_mat_mod.f90       | 79 ++++++++++++++++++++++++++-
 base/modules/psb_z_vect_mod.F90      | 69 ++++++++++++++++++++---
 base/serial/impl/psb_c_coo_impl.f90  |  2 +-
 base/serial/impl/psb_c_mat_impl.F90  | 16 +++++-
 base/serial/impl/psb_d_coo_impl.f90  |  2 +-
 base/serial/impl/psb_d_mat_impl.F90  | 16 +++++-
 base/serial/impl/psb_s_coo_impl.f90  |  2 +-
 base/serial/impl/psb_s_mat_impl.F90  | 16 +++++-
 base/serial/impl/psb_z_coo_impl.f90  |  2 +-
 base/serial/impl/psb_z_mat_impl.F90  | 16 +++++-
 base/serial/psb_cspspmm.f90          |  7 +++
 base/serial/psb_dspspmm.f90          |  7 +++
 base/serial/psb_sspspmm.f90          |  7 +++
 base/serial/psb_zspspmm.f90          |  7 +++
 31 files changed, 847 insertions(+), 67 deletions(-)

diff --git a/base/modules/psb_c_csc_mat_mod.f90 b/base/modules/psb_c_csc_mat_mod.f90
index 3770bcd3..17fbc83a 100644
--- a/base/modules/psb_c_csc_mat_mod.f90
+++ b/base/modules/psb_c_csc_mat_mod.f90
@@ -60,6 +60,7 @@ module psb_c_csc_mat_mod
     complex(psb_spk_), allocatable :: val(:)
 
   contains
+    procedure, pass(a) :: is_by_cols  => c_csc_is_by_cols
     procedure, pass(a) :: get_size    => c_csc_get_size
     procedure, pass(a) :: get_nzeros  => c_csc_get_nzeros
     procedure, nopass  :: get_fmt     => c_csc_get_fmt
@@ -517,6 +518,15 @@ contains
   ! == ===================================
 
   
+  function c_csc_is_by_cols(a) result(res)
+    implicit none 
+    class(psb_c_csc_sparse_mat), intent(in) :: a
+    logical  :: res
+    res = .true.
+     
+  end function c_csc_is_by_cols
+
+  
   function c_csc_sizeof(a) result(res)
     implicit none 
     class(psb_c_csc_sparse_mat), intent(in) :: a
diff --git a/base/modules/psb_c_csr_mat_mod.f90 b/base/modules/psb_c_csr_mat_mod.f90
index 4f093e96..9df858c7 100644
--- a/base/modules/psb_c_csr_mat_mod.f90
+++ b/base/modules/psb_c_csr_mat_mod.f90
@@ -61,6 +61,7 @@ module psb_c_csr_mat_mod
     complex(psb_spk_), allocatable :: val(:)
 
   contains
+    procedure, pass(a) :: is_by_rows  => c_csr_is_by_rows
     procedure, pass(a) :: get_size    => c_csr_get_size
     procedure, pass(a) :: get_nzeros  => c_csr_get_nzeros
     procedure, nopass  :: get_fmt     => c_csr_get_fmt
@@ -102,7 +103,8 @@ module psb_c_csr_mat_mod
   end type psb_c_csr_sparse_mat
 
   private :: c_csr_get_nzeros, c_csr_free,  c_csr_get_fmt, &
-       & c_csr_get_size, c_csr_sizeof, c_csr_get_nz_row
+       & c_csr_get_size, c_csr_sizeof, c_csr_get_nz_row, &
+       & c_csr_is_by_rows
 
   !> \memberof psb_c_csr_sparse_mat
   !| \see psb_base_mat_mod::psb_base_reallocate_nz
@@ -519,6 +521,16 @@ contains
   !
   ! == ===================================
 
+
+  
+  function c_csr_is_by_rows(a) result(res)
+    implicit none 
+    class(psb_c_csr_sparse_mat), intent(in) :: a
+    logical  :: res
+    res = .true.
+     
+  end function c_csr_is_by_rows
+
   
   function c_csr_sizeof(a) result(res)
     implicit none 
diff --git a/base/modules/psb_c_mat_mod.f90 b/base/modules/psb_c_mat_mod.f90
index 0b5452c7..78f5ffb3 100644
--- a/base/modules/psb_c_mat_mod.f90
+++ b/base/modules/psb_c_mat_mod.f90
@@ -92,6 +92,8 @@ module psb_c_mat_mod
     procedure, pass(a) :: is_upd      => psb_c_is_upd
     procedure, pass(a) :: is_asb      => psb_c_is_asb
     procedure, pass(a) :: is_sorted   => psb_c_is_sorted
+    procedure, pass(a) :: is_by_rows  => psb_c_is_by_rows
+    procedure, pass(a) :: is_by_cols  => psb_c_is_by_cols
     procedure, pass(a) :: is_upper    => psb_c_is_upper
     procedure, pass(a) :: is_lower    => psb_c_is_lower
     procedure, pass(a) :: is_triangle => psb_c_is_triangle
@@ -183,9 +185,21 @@ module psb_c_mat_mod
 
   private :: psb_c_get_nrows, psb_c_get_ncols, psb_c_get_nzeros, psb_c_get_size, &
        & psb_c_get_dupl, psb_c_is_null, psb_c_is_bld, &
-       & psb_c_is_upd, psb_c_is_asb, psb_c_is_sorted, psb_c_is_upper, &
+       & psb_c_is_upd, psb_c_is_asb, psb_c_is_sorted, &
+       & psb_c_is_by_rows, psb_c_is_by_cols, psb_c_is_upper, &
        & psb_c_is_lower, psb_c_is_triangle, psb_c_get_nz_row
 
+  class(psb_c_base_sparse_mat), allocatable, target, &
+       & save, private :: psb_c_base_mat_default
+
+  interface psb_set_mat_default
+    module procedure psb_c_set_mat_default
+  end interface
+
+  interface psb_get_mat_default
+    module procedure psb_c_get_mat_default
+  end interface
+
   interface psb_sizeof
     module procedure psb_c_sizeof
   end interface
@@ -812,6 +826,43 @@ module psb_c_mat_mod
 contains 
 
 
+  
+  subroutine  psb_c_set_mat_default(a) 
+    implicit none 
+    class(psb_c_base_sparse_mat), intent(in) :: a
+    
+    if (allocated(psb_c_base_mat_default)) then 
+      deallocate(psb_c_base_mat_default)
+    end if
+    allocate(psb_c_base_mat_default, mold=a)
+
+  end subroutine psb_c_set_mat_default
+  
+  function psb_c_get_mat_default(a) result(res)
+    implicit none 
+    class(psb_cspmat_type), intent(in) :: a
+    class(psb_c_base_sparse_mat), pointer :: res
+    
+    res => psb_c_get_base_mat_default()
+    
+  end function psb_c_get_mat_default
+
+  
+  function psb_c_get_base_mat_default() result(res)
+    implicit none 
+    class(psb_c_base_sparse_mat), pointer :: res
+    
+    if (.not.allocated(psb_c_base_mat_default)) then 
+      allocate(psb_c_csr_sparse_mat :: psb_c_base_mat_default)
+    end if
+
+    res => psb_c_base_mat_default
+    
+  end function psb_c_get_base_mat_default
+
+
+
+
   ! == ===================================
   !
   !
@@ -1007,6 +1058,32 @@ contains
 
   end function psb_c_is_sorted
 
+  function psb_c_is_by_rows(a) result(res)
+    implicit none 
+    class(psb_cspmat_type), intent(in) :: a
+    logical :: res
+
+    if (allocated(a%a)) then 
+      res = a%a%is_by_rows()
+    else
+      res = .false.
+    end if
+
+  end function psb_c_is_by_rows
+
+  function psb_c_is_by_cols(a) result(res)
+    implicit none 
+    class(psb_cspmat_type), intent(in) :: a
+    logical :: res
+
+    if (allocated(a%a)) then 
+      res = a%a%is_by_cols()
+    else
+      res = .false.
+    end if
+
+  end function psb_c_is_by_cols
+
 
 
   function psb_c_get_nzeros(a) result(res)
diff --git a/base/modules/psb_c_vect_mod.F90 b/base/modules/psb_c_vect_mod.F90
index 3848d9ff..83dca917 100644
--- a/base/modules/psb_c_vect_mod.F90
+++ b/base/modules/psb_c_vect_mod.F90
@@ -93,10 +93,56 @@ module psb_c_vect_mod
   interface psb_c_vect
     module procedure constructor, size_const
   end interface psb_c_vect
+ 
+  class(psb_c_base_vect_type), allocatable, target,&
+       & save, private :: psb_c_base_vect_default
+
+  interface psb_set_vect_default
+    module procedure psb_c_set_vect_default
+  end interface
+
+  interface psb_get_vect_default
+    module procedure psb_c_get_vect_default
+  end interface
+
 
 contains
 
   
+  subroutine  psb_c_set_vect_default(v) 
+    implicit none 
+    class(psb_c_base_vect_type), intent(in) :: v
+    
+    if (allocated(psb_c_base_vect_default)) then 
+      deallocate(psb_c_base_vect_default)
+    end if
+    allocate(psb_c_base_vect_default, mold=v)
+
+  end subroutine psb_c_set_vect_default
+  
+  function psb_c_get_vect_default(v) result(res)
+    implicit none 
+    class(psb_c_vect_type), intent(in) :: v
+    class(psb_c_base_vect_type), pointer :: res
+    
+    res => psb_c_get_base_vect_default()
+    
+  end function psb_c_get_vect_default
+
+  
+  function psb_c_get_base_vect_default() result(res)
+    implicit none 
+    class(psb_c_base_vect_type), pointer :: res
+    
+    if (.not.allocated(psb_c_base_vect_default)) then 
+      allocate(psb_c_base_vect_type :: psb_c_base_vect_default)
+    end if
+
+    res => psb_c_base_vect_default
+    
+  end function psb_c_get_base_vect_default
+
+  
   subroutine c_vect_clone(x,y,info)
     implicit none 
     class(psb_c_vect_type), intent(inout) :: x
@@ -115,6 +161,7 @@ contains
     class(psb_c_vect_type), intent(out) :: x
     class(psb_c_base_vect_type), intent(in), optional :: mold
     integer(psb_ipk_) :: info
+    class(psb_c_base_vect_type), pointer :: mld
 
     if (present(mold)) then 
 #ifdef HAVE_MOLD
@@ -123,7 +170,12 @@ contains
       call mold%mold(x%v,info)
 #endif
     else
-      allocate(psb_c_base_vect_type :: x%v,stat=info)
+#ifdef HAVE_MOLD
+      allocate(x%v,stat=info, mold=psb_c_get_base_vect_default())
+#else 
+      mld = psb_c_get_base_vect_default()
+      call mld%mold(x%v,info)
+#endif
     endif
 
     if (info == psb_success_) call x%v%bld(invect)
@@ -136,6 +188,7 @@ contains
     class(psb_c_vect_type), intent(out) :: x
     class(psb_c_base_vect_type), intent(in), optional :: mold
     integer(psb_ipk_) :: info
+    class(psb_c_base_vect_type), pointer :: mld
 
     if (present(mold)) then 
 #ifdef HAVE_MOLD
@@ -144,7 +197,12 @@ contains
       call mold%mold(x%v,info)
 #endif
     else
-      allocate(psb_c_base_vect_type :: x%v,stat=info)
+#ifdef HAVE_MOLD
+      allocate(x%v,stat=info, mold=psb_c_get_base_vect_default())
+#else 
+      mld = psb_c_get_base_vect_default()
+      call mld%mold(x%v,info)
+#endif
     endif
     if (info == psb_success_) call x%v%bld(n)
 
@@ -184,10 +242,7 @@ contains
     type(psb_c_vect_type) :: this
     integer(psb_ipk_) :: info
 
-    allocate(psb_c_base_vect_type :: this%v, stat=info)
-
-    if (info == 0) call this%v%bld(x)
-
+    call this%bld(x)
     call this%asb(size(x,kind=psb_ipk_),info)
 
   end function constructor
@@ -198,7 +253,7 @@ contains
     type(psb_c_vect_type) :: this
     integer(psb_ipk_) :: info
 
-    allocate(psb_c_base_vect_type :: this%v, stat=info)
+    call this%bld(n)
     call this%asb(n,info)
 
   end function size_const
diff --git a/base/modules/psb_const_mod.F90 b/base/modules/psb_const_mod.F90
index e305336a..77bd1e02 100644
--- a/base/modules/psb_const_mod.F90
+++ b/base/modules/psb_const_mod.F90
@@ -113,8 +113,9 @@ module psb_const_mod
   integer(psb_ipk_), parameter :: psb_iflag_=2, psb_ichk_=3
   integer(psb_ipk_), parameter :: psb_nnzt_=4, psb_zero_=5,psb_ipc_=6
 
-  integer(psb_ipk_), parameter :: psb_row_major_ = 0 
-  integer(psb_ipk_), parameter :: psb_col_major_ = 1
+  integer(psb_ipk_), parameter :: psb_unsorted_  = 0 
+  integer(psb_ipk_), parameter :: psb_row_major_ = 1 
+  integer(psb_ipk_), parameter :: psb_col_major_ = 2
   
   ! Duplicate coefficients handling
   ! These are usually set while calling spcnv as one of its
diff --git a/base/modules/psb_d_csc_mat_mod.f90 b/base/modules/psb_d_csc_mat_mod.f90
index 66aee0fe..e64e787a 100644
--- a/base/modules/psb_d_csc_mat_mod.f90
+++ b/base/modules/psb_d_csc_mat_mod.f90
@@ -60,6 +60,7 @@ module psb_d_csc_mat_mod
     real(psb_dpk_), allocatable :: val(:)
 
   contains
+    procedure, pass(a) :: is_by_cols  => d_csc_is_by_cols
     procedure, pass(a) :: get_size    => d_csc_get_size
     procedure, pass(a) :: get_nzeros  => d_csc_get_nzeros
     procedure, nopass  :: get_fmt     => d_csc_get_fmt
@@ -517,6 +518,15 @@ contains
   ! == ===================================
 
   
+  function d_csc_is_by_cols(a) result(res)
+    implicit none 
+    class(psb_d_csc_sparse_mat), intent(in) :: a
+    logical  :: res
+    res = .true.
+     
+  end function d_csc_is_by_cols
+
+  
   function d_csc_sizeof(a) result(res)
     implicit none 
     class(psb_d_csc_sparse_mat), intent(in) :: a
diff --git a/base/modules/psb_d_csr_mat_mod.f90 b/base/modules/psb_d_csr_mat_mod.f90
index f2d4fb53..25fa9c39 100644
--- a/base/modules/psb_d_csr_mat_mod.f90
+++ b/base/modules/psb_d_csr_mat_mod.f90
@@ -61,6 +61,7 @@ module psb_d_csr_mat_mod
     real(psb_dpk_), allocatable :: val(:)
 
   contains
+    procedure, pass(a) :: is_by_rows  => d_csr_is_by_rows
     procedure, pass(a) :: get_size    => d_csr_get_size
     procedure, pass(a) :: get_nzeros  => d_csr_get_nzeros
     procedure, nopass  :: get_fmt     => d_csr_get_fmt
@@ -102,7 +103,8 @@ module psb_d_csr_mat_mod
   end type psb_d_csr_sparse_mat
 
   private :: d_csr_get_nzeros, d_csr_free,  d_csr_get_fmt, &
-       & d_csr_get_size, d_csr_sizeof, d_csr_get_nz_row
+       & d_csr_get_size, d_csr_sizeof, d_csr_get_nz_row, &
+       & d_csr_is_by_rows
 
   !> \memberof psb_d_csr_sparse_mat
   !| \see psb_base_mat_mod::psb_base_reallocate_nz
@@ -519,6 +521,16 @@ contains
   !
   ! == ===================================
 
+
+  
+  function d_csr_is_by_rows(a) result(res)
+    implicit none 
+    class(psb_d_csr_sparse_mat), intent(in) :: a
+    logical  :: res
+    res = .true.
+     
+  end function d_csr_is_by_rows
+
   
   function d_csr_sizeof(a) result(res)
     implicit none 
diff --git a/base/modules/psb_d_mat_mod.f90 b/base/modules/psb_d_mat_mod.f90
index 1b9b8e3c..5f9eb689 100644
--- a/base/modules/psb_d_mat_mod.f90
+++ b/base/modules/psb_d_mat_mod.f90
@@ -92,6 +92,8 @@ module psb_d_mat_mod
     procedure, pass(a) :: is_upd      => psb_d_is_upd
     procedure, pass(a) :: is_asb      => psb_d_is_asb
     procedure, pass(a) :: is_sorted   => psb_d_is_sorted
+    procedure, pass(a) :: is_by_rows  => psb_d_is_by_rows
+    procedure, pass(a) :: is_by_cols  => psb_d_is_by_cols
     procedure, pass(a) :: is_upper    => psb_d_is_upper
     procedure, pass(a) :: is_lower    => psb_d_is_lower
     procedure, pass(a) :: is_triangle => psb_d_is_triangle
@@ -183,9 +185,21 @@ module psb_d_mat_mod
 
   private :: psb_d_get_nrows, psb_d_get_ncols, psb_d_get_nzeros, psb_d_get_size, &
        & psb_d_get_dupl, psb_d_is_null, psb_d_is_bld, &
-       & psb_d_is_upd, psb_d_is_asb, psb_d_is_sorted, psb_d_is_upper, &
+       & psb_d_is_upd, psb_d_is_asb, psb_d_is_sorted, &
+       & psb_d_is_by_rows, psb_d_is_by_cols, psb_d_is_upper, &
        & psb_d_is_lower, psb_d_is_triangle, psb_d_get_nz_row
 
+  class(psb_d_base_sparse_mat), allocatable, target, &
+       & save, private :: psb_d_base_mat_default
+
+  interface psb_set_mat_default
+    module procedure psb_d_set_mat_default
+  end interface
+
+  interface psb_get_mat_default
+    module procedure psb_d_get_mat_default
+  end interface
+
   interface psb_sizeof
     module procedure psb_d_sizeof
   end interface
@@ -812,6 +826,43 @@ module psb_d_mat_mod
 contains 
 
 
+  
+  subroutine  psb_d_set_mat_default(a) 
+    implicit none 
+    class(psb_d_base_sparse_mat), intent(in) :: a
+    
+    if (allocated(psb_d_base_mat_default)) then 
+      deallocate(psb_d_base_mat_default)
+    end if
+    allocate(psb_d_base_mat_default, mold=a)
+
+  end subroutine psb_d_set_mat_default
+  
+  function psb_d_get_mat_default(a) result(res)
+    implicit none 
+    class(psb_dspmat_type), intent(in) :: a
+    class(psb_d_base_sparse_mat), pointer :: res
+    
+    res => psb_d_get_base_mat_default()
+    
+  end function psb_d_get_mat_default
+
+  
+  function psb_d_get_base_mat_default() result(res)
+    implicit none 
+    class(psb_d_base_sparse_mat), pointer :: res
+    
+    if (.not.allocated(psb_d_base_mat_default)) then 
+      allocate(psb_d_csr_sparse_mat :: psb_d_base_mat_default)
+    end if
+
+    res => psb_d_base_mat_default
+    
+  end function psb_d_get_base_mat_default
+
+
+
+
   ! == ===================================
   !
   !
@@ -1007,6 +1058,32 @@ contains
 
   end function psb_d_is_sorted
 
+  function psb_d_is_by_rows(a) result(res)
+    implicit none 
+    class(psb_dspmat_type), intent(in) :: a
+    logical :: res
+
+    if (allocated(a%a)) then 
+      res = a%a%is_by_rows()
+    else
+      res = .false.
+    end if
+
+  end function psb_d_is_by_rows
+
+  function psb_d_is_by_cols(a) result(res)
+    implicit none 
+    class(psb_dspmat_type), intent(in) :: a
+    logical :: res
+
+    if (allocated(a%a)) then 
+      res = a%a%is_by_cols()
+    else
+      res = .false.
+    end if
+
+  end function psb_d_is_by_cols
+
 
 
   function psb_d_get_nzeros(a) result(res)
diff --git a/base/modules/psb_d_vect_mod.F90 b/base/modules/psb_d_vect_mod.F90
index 1b00c678..caf164d0 100644
--- a/base/modules/psb_d_vect_mod.F90
+++ b/base/modules/psb_d_vect_mod.F90
@@ -93,10 +93,56 @@ module psb_d_vect_mod
   interface psb_d_vect
     module procedure constructor, size_const
   end interface psb_d_vect
+ 
+  class(psb_d_base_vect_type), allocatable, target,&
+       & save, private :: psb_d_base_vect_default
+
+  interface psb_set_vect_default
+    module procedure psb_d_set_vect_default
+  end interface
+
+  interface psb_get_vect_default
+    module procedure psb_d_get_vect_default
+  end interface
+
 
 contains
 
   
+  subroutine  psb_d_set_vect_default(v) 
+    implicit none 
+    class(psb_d_base_vect_type), intent(in) :: v
+    
+    if (allocated(psb_d_base_vect_default)) then 
+      deallocate(psb_d_base_vect_default)
+    end if
+    allocate(psb_d_base_vect_default, mold=v)
+
+  end subroutine psb_d_set_vect_default
+  
+  function psb_d_get_vect_default(v) result(res)
+    implicit none 
+    class(psb_d_vect_type), intent(in) :: v
+    class(psb_d_base_vect_type), pointer :: res
+    
+    res => psb_d_get_base_vect_default()
+    
+  end function psb_d_get_vect_default
+
+  
+  function psb_d_get_base_vect_default() result(res)
+    implicit none 
+    class(psb_d_base_vect_type), pointer :: res
+    
+    if (.not.allocated(psb_d_base_vect_default)) then 
+      allocate(psb_d_base_vect_type :: psb_d_base_vect_default)
+    end if
+
+    res => psb_d_base_vect_default
+    
+  end function psb_d_get_base_vect_default
+
+  
   subroutine d_vect_clone(x,y,info)
     implicit none 
     class(psb_d_vect_type), intent(inout) :: x
@@ -115,6 +161,7 @@ contains
     class(psb_d_vect_type), intent(out) :: x
     class(psb_d_base_vect_type), intent(in), optional :: mold
     integer(psb_ipk_) :: info
+    class(psb_d_base_vect_type), pointer :: mld
 
     if (present(mold)) then 
 #ifdef HAVE_MOLD
@@ -123,7 +170,12 @@ contains
       call mold%mold(x%v,info)
 #endif
     else
-      allocate(psb_d_base_vect_type :: x%v,stat=info)
+#ifdef HAVE_MOLD
+      allocate(x%v,stat=info, mold=psb_d_get_base_vect_default())
+#else 
+      mld = psb_d_get_base_vect_default()
+      call mld%mold(x%v,info)
+#endif
     endif
 
     if (info == psb_success_) call x%v%bld(invect)
@@ -136,6 +188,7 @@ contains
     class(psb_d_vect_type), intent(out) :: x
     class(psb_d_base_vect_type), intent(in), optional :: mold
     integer(psb_ipk_) :: info
+    class(psb_d_base_vect_type), pointer :: mld
 
     if (present(mold)) then 
 #ifdef HAVE_MOLD
@@ -144,7 +197,12 @@ contains
       call mold%mold(x%v,info)
 #endif
     else
-      allocate(psb_d_base_vect_type :: x%v,stat=info)
+#ifdef HAVE_MOLD
+      allocate(x%v,stat=info, mold=psb_d_get_base_vect_default())
+#else 
+      mld = psb_d_get_base_vect_default()
+      call mld%mold(x%v,info)
+#endif
     endif
     if (info == psb_success_) call x%v%bld(n)
 
@@ -184,10 +242,7 @@ contains
     type(psb_d_vect_type) :: this
     integer(psb_ipk_) :: info
 
-    allocate(psb_d_base_vect_type :: this%v, stat=info)
-
-    if (info == 0) call this%v%bld(x)
-
+    call this%bld(x)
     call this%asb(size(x,kind=psb_ipk_),info)
 
   end function constructor
@@ -198,7 +253,7 @@ contains
     type(psb_d_vect_type) :: this
     integer(psb_ipk_) :: info
 
-    allocate(psb_d_base_vect_type :: this%v, stat=info)
+    call this%bld(n)
     call this%asb(n,info)
 
   end function size_const
diff --git a/base/modules/psb_i_base_vect_mod.f90 b/base/modules/psb_i_base_vect_mod.f90
index 9e81d16b..c188e0b7 100644
--- a/base/modules/psb_i_base_vect_mod.f90
+++ b/base/modules/psb_i_base_vect_mod.f90
@@ -141,9 +141,11 @@ module psb_i_base_vect_mod
     !
     procedure, pass(x) :: gthab    => i_base_gthab
     procedure, pass(x) :: gthzv    => i_base_gthzv
-    generic, public    :: gth      => gthab, gthzv
+    procedure, pass(x) :: gthzv_x  => i_base_gthzv_x
+    generic, public    :: gth      => gthab, gthzv, gthzv_x
     procedure, pass(y) :: sctb     => i_base_sctb
-    generic, public    :: sct      => sctb
+    procedure, pass(y) :: sctb_x   => i_base_sctb_x
+    generic, public    :: sct      => sctb, sctb_x
   end type psb_i_base_vect_type
 
   public  :: psb_i_base_vect
@@ -574,7 +576,7 @@ contains
   !    
   function  i_base_get_vect(x) result(res)
     class(psb_i_base_vect_type), intent(inout) :: x
-    integer(psb_ipk_), allocatable             :: res(:)
+    integer(psb_ipk_), allocatable                 :: res(:)
     integer(psb_ipk_) :: info
     
     if (.not.allocated(x%v)) return 
@@ -1042,6 +1044,26 @@ contains
     call psi_gth(n,idx,alpha,x%v,beta,y)
 
   end subroutine i_base_gthab
+  !
+  ! shortcut alpha=1 beta=0
+  ! 
+  !> Function  base_gthzv
+  !! \memberof  psb_i_base_vect_type
+  !! \brief gather into an array special alpha=1 beta=0
+  !!    Y = X(IDX(:))
+  !! \param n  how many entries to consider
+  !! \param idx(:) indices
+  subroutine i_base_gthzv_x(i,n,idx,x,y)
+    use psi_serial_mod
+    integer(psb_ipk_) :: i,n
+    class(psb_i_base_vect_type) :: idx
+    integer(psb_ipk_) ::  y(:)
+    class(psb_i_base_vect_type) :: x
+    
+    call x%gth(n,idx%v(i:),y)
+
+  end subroutine i_base_gthzv_x
+
   !
   ! shortcut alpha=1 beta=0
   ! 
@@ -1087,4 +1109,15 @@ contains
 
   end subroutine i_base_sctb
 
+  subroutine i_base_sctb_x(i,n,idx,x,beta,y)
+    use psi_serial_mod
+    integer(psb_ipk_) :: i, n
+    class(psb_i_base_vect_type) :: idx
+    integer(psb_ipk_) :: beta, x(:)
+    class(psb_i_base_vect_type) :: y
+    
+    call y%sct(n,idx%v(i:),x,beta)
+
+  end subroutine i_base_sctb_x
+
 end module psb_i_base_vect_mod
diff --git a/base/modules/psb_i_vect_mod.F90 b/base/modules/psb_i_vect_mod.F90
index 95000951..4b660383 100644
--- a/base/modules/psb_i_vect_mod.F90
+++ b/base/modules/psb_i_vect_mod.F90
@@ -93,10 +93,56 @@ module psb_i_vect_mod
   interface psb_i_vect
     module procedure constructor, size_const
   end interface psb_i_vect
+ 
+  class(psb_i_base_vect_type), allocatable, target,&
+       & save, private :: psb_i_base_vect_default
+
+  interface psb_set_vect_default
+    module procedure psb_i_set_vect_default
+  end interface
+
+  interface psb_get_vect_default
+    module procedure psb_i_get_vect_default
+  end interface
+
 
 contains
 
   
+  subroutine  psb_i_set_vect_default(v) 
+    implicit none 
+    class(psb_i_base_vect_type), intent(in) :: v
+    
+    if (allocated(psb_i_base_vect_default)) then 
+      deallocate(psb_i_base_vect_default)
+    end if
+    allocate(psb_i_base_vect_default, mold=v)
+
+  end subroutine psb_i_set_vect_default
+  
+  function psb_i_get_vect_default(v) result(res)
+    implicit none 
+    class(psb_i_vect_type), intent(in) :: v
+    class(psb_i_base_vect_type), pointer :: res
+    
+    res => psb_i_get_base_vect_default()
+    
+  end function psb_i_get_vect_default
+
+  
+  function psb_i_get_base_vect_default() result(res)
+    implicit none 
+    class(psb_i_base_vect_type), pointer :: res
+    
+    if (.not.allocated(psb_i_base_vect_default)) then 
+      allocate(psb_i_base_vect_type :: psb_i_base_vect_default)
+    end if
+
+    res => psb_i_base_vect_default
+    
+  end function psb_i_get_base_vect_default
+
+  
   subroutine i_vect_clone(x,y,info)
     implicit none 
     class(psb_i_vect_type), intent(inout) :: x
@@ -115,6 +161,7 @@ contains
     class(psb_i_vect_type), intent(out) :: x
     class(psb_i_base_vect_type), intent(in), optional :: mold
     integer(psb_ipk_) :: info
+    class(psb_i_base_vect_type), pointer :: mld
 
     if (present(mold)) then 
 #ifdef HAVE_MOLD
@@ -123,7 +170,12 @@ contains
       call mold%mold(x%v,info)
 #endif
     else
-      allocate(psb_i_base_vect_type :: x%v,stat=info)
+#ifdef HAVE_MOLD
+      allocate(x%v,stat=info, mold=psb_i_get_base_vect_default())
+#else 
+      mld = psb_i_get_base_vect_default()
+      call mld%mold(x%v,info)
+#endif
     endif
 
     if (info == psb_success_) call x%v%bld(invect)
@@ -136,6 +188,7 @@ contains
     class(psb_i_vect_type), intent(out) :: x
     class(psb_i_base_vect_type), intent(in), optional :: mold
     integer(psb_ipk_) :: info
+    class(psb_i_base_vect_type), pointer :: mld
 
     if (present(mold)) then 
 #ifdef HAVE_MOLD
@@ -144,7 +197,12 @@ contains
       call mold%mold(x%v,info)
 #endif
     else
-      allocate(psb_i_base_vect_type :: x%v,stat=info)
+#ifdef HAVE_MOLD
+      allocate(x%v,stat=info, mold=psb_i_get_base_vect_default())
+#else 
+      mld = psb_i_get_base_vect_default()
+      call mld%mold(x%v,info)
+#endif
     endif
     if (info == psb_success_) call x%v%bld(n)
 
@@ -152,7 +210,7 @@ contains
 
   function  i_vect_get_vect(x) result(res)
     class(psb_i_vect_type), intent(inout)  :: x
-    integer(psb_ipk_), allocatable         :: res(:)
+    integer(psb_ipk_), allocatable                 :: res(:)
     integer(psb_ipk_) :: info
 
     if (allocated(x%v)) then
@@ -184,10 +242,7 @@ contains
     type(psb_i_vect_type) :: this
     integer(psb_ipk_) :: info
 
-    allocate(psb_i_base_vect_type :: this%v, stat=info)
-
-    if (info == 0) call this%v%bld(x)
-
+    call this%bld(x)
     call this%asb(size(x,kind=psb_ipk_),info)
 
   end function constructor
@@ -198,7 +253,7 @@ contains
     type(psb_i_vect_type) :: this
     integer(psb_ipk_) :: info
 
-    allocate(psb_i_base_vect_type :: this%v, stat=info)
+    call this%bld(n)
     call this%asb(n,info)
 
   end function size_const
@@ -558,10 +613,10 @@ contains
     use psi_serial_mod
     implicit none 
     class(psb_i_vect_type), intent(inout)  :: x
-    integer(psb_ipk_), intent(in)          :: n, dupl
-    integer(psb_ipk_), intent(in)          :: irl(:)
-    integer(psb_ipk_), intent(in)          :: val(:)
-    integer(psb_ipk_), intent(out)         :: info
+    integer(psb_ipk_), intent(in)               :: n, dupl
+    integer(psb_ipk_), intent(in)               :: irl(:)
+    integer(psb_ipk_), intent(in)        :: val(:)
+    integer(psb_ipk_), intent(out)              :: info
 
     integer(psb_ipk_) :: i
 
@@ -580,8 +635,7 @@ contains
     class(psb_i_vect_type), intent(inout) :: x
     class(psb_i_base_vect_type), intent(in), optional :: mold
     class(psb_i_base_vect_type), allocatable :: tmp
-    integer(psb_ipk_)                        :: info
-
+    integer(psb_ipk_) :: info
 
     if (present(mold)) then 
 #ifdef HAVE_MOLD
diff --git a/base/modules/psb_s_csc_mat_mod.f90 b/base/modules/psb_s_csc_mat_mod.f90
index 29f5de43..2508bb49 100644
--- a/base/modules/psb_s_csc_mat_mod.f90
+++ b/base/modules/psb_s_csc_mat_mod.f90
@@ -60,6 +60,7 @@ module psb_s_csc_mat_mod
     real(psb_spk_), allocatable :: val(:)
 
   contains
+    procedure, pass(a) :: is_by_cols  => s_csc_is_by_cols
     procedure, pass(a) :: get_size    => s_csc_get_size
     procedure, pass(a) :: get_nzeros  => s_csc_get_nzeros
     procedure, nopass  :: get_fmt     => s_csc_get_fmt
@@ -517,6 +518,15 @@ contains
   ! == ===================================
 
   
+  function s_csc_is_by_cols(a) result(res)
+    implicit none 
+    class(psb_s_csc_sparse_mat), intent(in) :: a
+    logical  :: res
+    res = .true.
+     
+  end function s_csc_is_by_cols
+
+  
   function s_csc_sizeof(a) result(res)
     implicit none 
     class(psb_s_csc_sparse_mat), intent(in) :: a
diff --git a/base/modules/psb_s_csr_mat_mod.f90 b/base/modules/psb_s_csr_mat_mod.f90
index 47bf8266..a1d57e3c 100644
--- a/base/modules/psb_s_csr_mat_mod.f90
+++ b/base/modules/psb_s_csr_mat_mod.f90
@@ -61,6 +61,7 @@ module psb_s_csr_mat_mod
     real(psb_spk_), allocatable :: val(:)
 
   contains
+    procedure, pass(a) :: is_by_rows  => s_csr_is_by_rows
     procedure, pass(a) :: get_size    => s_csr_get_size
     procedure, pass(a) :: get_nzeros  => s_csr_get_nzeros
     procedure, nopass  :: get_fmt     => s_csr_get_fmt
@@ -102,7 +103,8 @@ module psb_s_csr_mat_mod
   end type psb_s_csr_sparse_mat
 
   private :: s_csr_get_nzeros, s_csr_free,  s_csr_get_fmt, &
-       & s_csr_get_size, s_csr_sizeof, s_csr_get_nz_row
+       & s_csr_get_size, s_csr_sizeof, s_csr_get_nz_row, &
+       & s_csr_is_by_rows
 
   !> \memberof psb_s_csr_sparse_mat
   !| \see psb_base_mat_mod::psb_base_reallocate_nz
@@ -519,6 +521,16 @@ contains
   !
   ! == ===================================
 
+
+  
+  function s_csr_is_by_rows(a) result(res)
+    implicit none 
+    class(psb_s_csr_sparse_mat), intent(in) :: a
+    logical  :: res
+    res = .true.
+     
+  end function s_csr_is_by_rows
+
   
   function s_csr_sizeof(a) result(res)
     implicit none 
diff --git a/base/modules/psb_s_mat_mod.f90 b/base/modules/psb_s_mat_mod.f90
index 79fecc37..d84db53a 100644
--- a/base/modules/psb_s_mat_mod.f90
+++ b/base/modules/psb_s_mat_mod.f90
@@ -92,6 +92,8 @@ module psb_s_mat_mod
     procedure, pass(a) :: is_upd      => psb_s_is_upd
     procedure, pass(a) :: is_asb      => psb_s_is_asb
     procedure, pass(a) :: is_sorted   => psb_s_is_sorted
+    procedure, pass(a) :: is_by_rows  => psb_s_is_by_rows
+    procedure, pass(a) :: is_by_cols  => psb_s_is_by_cols
     procedure, pass(a) :: is_upper    => psb_s_is_upper
     procedure, pass(a) :: is_lower    => psb_s_is_lower
     procedure, pass(a) :: is_triangle => psb_s_is_triangle
@@ -183,9 +185,21 @@ module psb_s_mat_mod
 
   private :: psb_s_get_nrows, psb_s_get_ncols, psb_s_get_nzeros, psb_s_get_size, &
        & psb_s_get_dupl, psb_s_is_null, psb_s_is_bld, &
-       & psb_s_is_upd, psb_s_is_asb, psb_s_is_sorted, psb_s_is_upper, &
+       & psb_s_is_upd, psb_s_is_asb, psb_s_is_sorted, &
+       & psb_s_is_by_rows, psb_s_is_by_cols, psb_s_is_upper, &
        & psb_s_is_lower, psb_s_is_triangle, psb_s_get_nz_row
 
+  class(psb_s_base_sparse_mat), allocatable, target, &
+       & save, private :: psb_s_base_mat_default
+
+  interface psb_set_mat_default
+    module procedure psb_s_set_mat_default
+  end interface
+
+  interface psb_get_mat_default
+    module procedure psb_s_get_mat_default
+  end interface
+
   interface psb_sizeof
     module procedure psb_s_sizeof
   end interface
@@ -812,6 +826,43 @@ module psb_s_mat_mod
 contains 
 
 
+  
+  subroutine  psb_s_set_mat_default(a) 
+    implicit none 
+    class(psb_s_base_sparse_mat), intent(in) :: a
+    
+    if (allocated(psb_s_base_mat_default)) then 
+      deallocate(psb_s_base_mat_default)
+    end if
+    allocate(psb_s_base_mat_default, mold=a)
+
+  end subroutine psb_s_set_mat_default
+  
+  function psb_s_get_mat_default(a) result(res)
+    implicit none 
+    class(psb_sspmat_type), intent(in) :: a
+    class(psb_s_base_sparse_mat), pointer :: res
+    
+    res => psb_s_get_base_mat_default()
+    
+  end function psb_s_get_mat_default
+
+  
+  function psb_s_get_base_mat_default() result(res)
+    implicit none 
+    class(psb_s_base_sparse_mat), pointer :: res
+    
+    if (.not.allocated(psb_s_base_mat_default)) then 
+      allocate(psb_s_csr_sparse_mat :: psb_s_base_mat_default)
+    end if
+
+    res => psb_s_base_mat_default
+    
+  end function psb_s_get_base_mat_default
+
+
+
+
   ! == ===================================
   !
   !
@@ -1007,6 +1058,32 @@ contains
 
   end function psb_s_is_sorted
 
+  function psb_s_is_by_rows(a) result(res)
+    implicit none 
+    class(psb_sspmat_type), intent(in) :: a
+    logical :: res
+
+    if (allocated(a%a)) then 
+      res = a%a%is_by_rows()
+    else
+      res = .false.
+    end if
+
+  end function psb_s_is_by_rows
+
+  function psb_s_is_by_cols(a) result(res)
+    implicit none 
+    class(psb_sspmat_type), intent(in) :: a
+    logical :: res
+
+    if (allocated(a%a)) then 
+      res = a%a%is_by_cols()
+    else
+      res = .false.
+    end if
+
+  end function psb_s_is_by_cols
+
 
 
   function psb_s_get_nzeros(a) result(res)
diff --git a/base/modules/psb_s_vect_mod.F90 b/base/modules/psb_s_vect_mod.F90
index fe418b1a..f5e82bf2 100644
--- a/base/modules/psb_s_vect_mod.F90
+++ b/base/modules/psb_s_vect_mod.F90
@@ -93,10 +93,56 @@ module psb_s_vect_mod
   interface psb_s_vect
     module procedure constructor, size_const
   end interface psb_s_vect
+ 
+  class(psb_s_base_vect_type), allocatable, target,&
+       & save, private :: psb_s_base_vect_default
+
+  interface psb_set_vect_default
+    module procedure psb_s_set_vect_default
+  end interface
+
+  interface psb_get_vect_default
+    module procedure psb_s_get_vect_default
+  end interface
+
 
 contains
 
   
+  subroutine  psb_s_set_vect_default(v) 
+    implicit none 
+    class(psb_s_base_vect_type), intent(in) :: v
+    
+    if (allocated(psb_s_base_vect_default)) then 
+      deallocate(psb_s_base_vect_default)
+    end if
+    allocate(psb_s_base_vect_default, mold=v)
+
+  end subroutine psb_s_set_vect_default
+  
+  function psb_s_get_vect_default(v) result(res)
+    implicit none 
+    class(psb_s_vect_type), intent(in) :: v
+    class(psb_s_base_vect_type), pointer :: res
+    
+    res => psb_s_get_base_vect_default()
+    
+  end function psb_s_get_vect_default
+
+  
+  function psb_s_get_base_vect_default() result(res)
+    implicit none 
+    class(psb_s_base_vect_type), pointer :: res
+    
+    if (.not.allocated(psb_s_base_vect_default)) then 
+      allocate(psb_s_base_vect_type :: psb_s_base_vect_default)
+    end if
+
+    res => psb_s_base_vect_default
+    
+  end function psb_s_get_base_vect_default
+
+  
   subroutine s_vect_clone(x,y,info)
     implicit none 
     class(psb_s_vect_type), intent(inout) :: x
@@ -115,6 +161,7 @@ contains
     class(psb_s_vect_type), intent(out) :: x
     class(psb_s_base_vect_type), intent(in), optional :: mold
     integer(psb_ipk_) :: info
+    class(psb_s_base_vect_type), pointer :: mld
 
     if (present(mold)) then 
 #ifdef HAVE_MOLD
@@ -123,7 +170,12 @@ contains
       call mold%mold(x%v,info)
 #endif
     else
-      allocate(psb_s_base_vect_type :: x%v,stat=info)
+#ifdef HAVE_MOLD
+      allocate(x%v,stat=info, mold=psb_s_get_base_vect_default())
+#else 
+      mld = psb_s_get_base_vect_default()
+      call mld%mold(x%v,info)
+#endif
     endif
 
     if (info == psb_success_) call x%v%bld(invect)
@@ -136,6 +188,7 @@ contains
     class(psb_s_vect_type), intent(out) :: x
     class(psb_s_base_vect_type), intent(in), optional :: mold
     integer(psb_ipk_) :: info
+    class(psb_s_base_vect_type), pointer :: mld
 
     if (present(mold)) then 
 #ifdef HAVE_MOLD
@@ -144,7 +197,12 @@ contains
       call mold%mold(x%v,info)
 #endif
     else
-      allocate(psb_s_base_vect_type :: x%v,stat=info)
+#ifdef HAVE_MOLD
+      allocate(x%v,stat=info, mold=psb_s_get_base_vect_default())
+#else 
+      mld = psb_s_get_base_vect_default()
+      call mld%mold(x%v,info)
+#endif
     endif
     if (info == psb_success_) call x%v%bld(n)
 
@@ -184,10 +242,7 @@ contains
     type(psb_s_vect_type) :: this
     integer(psb_ipk_) :: info
 
-    allocate(psb_s_base_vect_type :: this%v, stat=info)
-
-    if (info == 0) call this%v%bld(x)
-
+    call this%bld(x)
     call this%asb(size(x,kind=psb_ipk_),info)
 
   end function constructor
@@ -198,7 +253,7 @@ contains
     type(psb_s_vect_type) :: this
     integer(psb_ipk_) :: info
 
-    allocate(psb_s_base_vect_type :: this%v, stat=info)
+    call this%bld(n)
     call this%asb(n,info)
 
   end function size_const
diff --git a/base/modules/psb_z_csc_mat_mod.f90 b/base/modules/psb_z_csc_mat_mod.f90
index 3de58cbc..dc5dbf1c 100644
--- a/base/modules/psb_z_csc_mat_mod.f90
+++ b/base/modules/psb_z_csc_mat_mod.f90
@@ -60,6 +60,7 @@ module psb_z_csc_mat_mod
     complex(psb_dpk_), allocatable :: val(:)
 
   contains
+    procedure, pass(a) :: is_by_cols  => z_csc_is_by_cols
     procedure, pass(a) :: get_size    => z_csc_get_size
     procedure, pass(a) :: get_nzeros  => z_csc_get_nzeros
     procedure, nopass  :: get_fmt     => z_csc_get_fmt
@@ -517,6 +518,15 @@ contains
   ! == ===================================
 
   
+  function z_csc_is_by_cols(a) result(res)
+    implicit none 
+    class(psb_z_csc_sparse_mat), intent(in) :: a
+    logical  :: res
+    res = .true.
+     
+  end function z_csc_is_by_cols
+
+  
   function z_csc_sizeof(a) result(res)
     implicit none 
     class(psb_z_csc_sparse_mat), intent(in) :: a
diff --git a/base/modules/psb_z_csr_mat_mod.f90 b/base/modules/psb_z_csr_mat_mod.f90
index c5c9ff55..65b36351 100644
--- a/base/modules/psb_z_csr_mat_mod.f90
+++ b/base/modules/psb_z_csr_mat_mod.f90
@@ -61,6 +61,7 @@ module psb_z_csr_mat_mod
     complex(psb_dpk_), allocatable :: val(:)
 
   contains
+    procedure, pass(a) :: is_by_rows  => z_csr_is_by_rows
     procedure, pass(a) :: get_size    => z_csr_get_size
     procedure, pass(a) :: get_nzeros  => z_csr_get_nzeros
     procedure, nopass  :: get_fmt     => z_csr_get_fmt
@@ -102,7 +103,8 @@ module psb_z_csr_mat_mod
   end type psb_z_csr_sparse_mat
 
   private :: z_csr_get_nzeros, z_csr_free,  z_csr_get_fmt, &
-       & z_csr_get_size, z_csr_sizeof, z_csr_get_nz_row
+       & z_csr_get_size, z_csr_sizeof, z_csr_get_nz_row, &
+       & z_csr_is_by_rows
 
   !> \memberof psb_z_csr_sparse_mat
   !| \see psb_base_mat_mod::psb_base_reallocate_nz
@@ -519,6 +521,16 @@ contains
   !
   ! == ===================================
 
+
+  
+  function z_csr_is_by_rows(a) result(res)
+    implicit none 
+    class(psb_z_csr_sparse_mat), intent(in) :: a
+    logical  :: res
+    res = .true.
+     
+  end function z_csr_is_by_rows
+
   
   function z_csr_sizeof(a) result(res)
     implicit none 
diff --git a/base/modules/psb_z_mat_mod.f90 b/base/modules/psb_z_mat_mod.f90
index 28f352f8..887b2f36 100644
--- a/base/modules/psb_z_mat_mod.f90
+++ b/base/modules/psb_z_mat_mod.f90
@@ -92,6 +92,8 @@ module psb_z_mat_mod
     procedure, pass(a) :: is_upd      => psb_z_is_upd
     procedure, pass(a) :: is_asb      => psb_z_is_asb
     procedure, pass(a) :: is_sorted   => psb_z_is_sorted
+    procedure, pass(a) :: is_by_rows  => psb_z_is_by_rows
+    procedure, pass(a) :: is_by_cols  => psb_z_is_by_cols
     procedure, pass(a) :: is_upper    => psb_z_is_upper
     procedure, pass(a) :: is_lower    => psb_z_is_lower
     procedure, pass(a) :: is_triangle => psb_z_is_triangle
@@ -183,9 +185,21 @@ module psb_z_mat_mod
 
   private :: psb_z_get_nrows, psb_z_get_ncols, psb_z_get_nzeros, psb_z_get_size, &
        & psb_z_get_dupl, psb_z_is_null, psb_z_is_bld, &
-       & psb_z_is_upd, psb_z_is_asb, psb_z_is_sorted, psb_z_is_upper, &
+       & psb_z_is_upd, psb_z_is_asb, psb_z_is_sorted, &
+       & psb_z_is_by_rows, psb_z_is_by_cols, psb_z_is_upper, &
        & psb_z_is_lower, psb_z_is_triangle, psb_z_get_nz_row
 
+  class(psb_z_base_sparse_mat), allocatable, target, &
+       & save, private :: psb_z_base_mat_default
+
+  interface psb_set_mat_default
+    module procedure psb_z_set_mat_default
+  end interface
+
+  interface psb_get_mat_default
+    module procedure psb_z_get_mat_default
+  end interface
+
   interface psb_sizeof
     module procedure psb_z_sizeof
   end interface
@@ -812,6 +826,43 @@ module psb_z_mat_mod
 contains 
 
 
+  
+  subroutine  psb_z_set_mat_default(a) 
+    implicit none 
+    class(psb_z_base_sparse_mat), intent(in) :: a
+    
+    if (allocated(psb_z_base_mat_default)) then 
+      deallocate(psb_z_base_mat_default)
+    end if
+    allocate(psb_z_base_mat_default, mold=a)
+
+  end subroutine psb_z_set_mat_default
+  
+  function psb_z_get_mat_default(a) result(res)
+    implicit none 
+    class(psb_zspmat_type), intent(in) :: a
+    class(psb_z_base_sparse_mat), pointer :: res
+    
+    res => psb_z_get_base_mat_default()
+    
+  end function psb_z_get_mat_default
+
+  
+  function psb_z_get_base_mat_default() result(res)
+    implicit none 
+    class(psb_z_base_sparse_mat), pointer :: res
+    
+    if (.not.allocated(psb_z_base_mat_default)) then 
+      allocate(psb_z_csr_sparse_mat :: psb_z_base_mat_default)
+    end if
+
+    res => psb_z_base_mat_default
+    
+  end function psb_z_get_base_mat_default
+
+
+
+
   ! == ===================================
   !
   !
@@ -1007,6 +1058,32 @@ contains
 
   end function psb_z_is_sorted
 
+  function psb_z_is_by_rows(a) result(res)
+    implicit none 
+    class(psb_zspmat_type), intent(in) :: a
+    logical :: res
+
+    if (allocated(a%a)) then 
+      res = a%a%is_by_rows()
+    else
+      res = .false.
+    end if
+
+  end function psb_z_is_by_rows
+
+  function psb_z_is_by_cols(a) result(res)
+    implicit none 
+    class(psb_zspmat_type), intent(in) :: a
+    logical :: res
+
+    if (allocated(a%a)) then 
+      res = a%a%is_by_cols()
+    else
+      res = .false.
+    end if
+
+  end function psb_z_is_by_cols
+
 
 
   function psb_z_get_nzeros(a) result(res)
diff --git a/base/modules/psb_z_vect_mod.F90 b/base/modules/psb_z_vect_mod.F90
index 8e985c18..b901ba43 100644
--- a/base/modules/psb_z_vect_mod.F90
+++ b/base/modules/psb_z_vect_mod.F90
@@ -93,10 +93,56 @@ module psb_z_vect_mod
   interface psb_z_vect
     module procedure constructor, size_const
   end interface psb_z_vect
+ 
+  class(psb_z_base_vect_type), allocatable, target,&
+       & save, private :: psb_z_base_vect_default
+
+  interface psb_set_vect_default
+    module procedure psb_z_set_vect_default
+  end interface
+
+  interface psb_get_vect_default
+    module procedure psb_z_get_vect_default
+  end interface
+
 
 contains
 
   
+  subroutine  psb_z_set_vect_default(v) 
+    implicit none 
+    class(psb_z_base_vect_type), intent(in) :: v
+    
+    if (allocated(psb_z_base_vect_default)) then 
+      deallocate(psb_z_base_vect_default)
+    end if
+    allocate(psb_z_base_vect_default, mold=v)
+
+  end subroutine psb_z_set_vect_default
+  
+  function psb_z_get_vect_default(v) result(res)
+    implicit none 
+    class(psb_z_vect_type), intent(in) :: v
+    class(psb_z_base_vect_type), pointer :: res
+    
+    res => psb_z_get_base_vect_default()
+    
+  end function psb_z_get_vect_default
+
+  
+  function psb_z_get_base_vect_default() result(res)
+    implicit none 
+    class(psb_z_base_vect_type), pointer :: res
+    
+    if (.not.allocated(psb_z_base_vect_default)) then 
+      allocate(psb_z_base_vect_type :: psb_z_base_vect_default)
+    end if
+
+    res => psb_z_base_vect_default
+    
+  end function psb_z_get_base_vect_default
+
+  
   subroutine z_vect_clone(x,y,info)
     implicit none 
     class(psb_z_vect_type), intent(inout) :: x
@@ -115,6 +161,7 @@ contains
     class(psb_z_vect_type), intent(out) :: x
     class(psb_z_base_vect_type), intent(in), optional :: mold
     integer(psb_ipk_) :: info
+    class(psb_z_base_vect_type), pointer :: mld
 
     if (present(mold)) then 
 #ifdef HAVE_MOLD
@@ -123,7 +170,12 @@ contains
       call mold%mold(x%v,info)
 #endif
     else
-      allocate(psb_z_base_vect_type :: x%v,stat=info)
+#ifdef HAVE_MOLD
+      allocate(x%v,stat=info, mold=psb_z_get_base_vect_default())
+#else 
+      mld = psb_z_get_base_vect_default()
+      call mld%mold(x%v,info)
+#endif
     endif
 
     if (info == psb_success_) call x%v%bld(invect)
@@ -136,6 +188,7 @@ contains
     class(psb_z_vect_type), intent(out) :: x
     class(psb_z_base_vect_type), intent(in), optional :: mold
     integer(psb_ipk_) :: info
+    class(psb_z_base_vect_type), pointer :: mld
 
     if (present(mold)) then 
 #ifdef HAVE_MOLD
@@ -144,7 +197,12 @@ contains
       call mold%mold(x%v,info)
 #endif
     else
-      allocate(psb_z_base_vect_type :: x%v,stat=info)
+#ifdef HAVE_MOLD
+      allocate(x%v,stat=info, mold=psb_z_get_base_vect_default())
+#else 
+      mld = psb_z_get_base_vect_default()
+      call mld%mold(x%v,info)
+#endif
     endif
     if (info == psb_success_) call x%v%bld(n)
 
@@ -184,10 +242,7 @@ contains
     type(psb_z_vect_type) :: this
     integer(psb_ipk_) :: info
 
-    allocate(psb_z_base_vect_type :: this%v, stat=info)
-
-    if (info == 0) call this%v%bld(x)
-
+    call this%bld(x)
     call this%asb(size(x,kind=psb_ipk_),info)
 
   end function constructor
@@ -198,7 +253,7 @@ contains
     type(psb_z_vect_type) :: this
     integer(psb_ipk_) :: info
 
-    allocate(psb_z_base_vect_type :: this%v, stat=info)
+    call this%bld(n)
     call this%asb(n,info)
 
   end function size_const
diff --git a/base/serial/impl/psb_c_coo_impl.f90 b/base/serial/impl/psb_c_coo_impl.f90
index 4b328ef3..f20998fd 100644
--- a/base/serial/impl/psb_c_coo_impl.f90
+++ b/base/serial/impl/psb_c_coo_impl.f90
@@ -3375,7 +3375,7 @@ subroutine psb_c_fix_coo(a,info,idir)
   if (present(idir)) then 
     idir_ = idir
   else
-    idir_ = 0
+    idir_ = psb_row_major_
   endif
   
   nra = a%get_nrows()
diff --git a/base/serial/impl/psb_c_mat_impl.F90 b/base/serial/impl/psb_c_mat_impl.F90
index 55fbb02d..02535f2d 100644
--- a/base/serial/impl/psb_c_mat_impl.F90
+++ b/base/serial/impl/psb_c_mat_impl.F90
@@ -1164,6 +1164,7 @@ subroutine psb_c_cscnv(a,b,info,type,mold,upd,dupl)
 
 
   class(psb_c_base_sparse_mat), allocatable  :: altmp
+  class(psb_c_base_sparse_mat), pointer      :: mld
   integer(psb_ipk_) :: err_act
   character(len=20)  :: name='cscnv'
   logical, parameter :: debug=.false.
@@ -1205,7 +1206,12 @@ subroutine psb_c_cscnv(a,b,info,type,mold,upd,dupl)
       goto 9999
     end select
   else
-    allocate(psb_c_csr_sparse_mat :: altmp, stat=info) 
+#if defined(HAVE_MOLD)
+    allocate(altmp, mold=psb_get_mat_default(a),stat=info) 
+#else
+    mld = psb_get_mat_default(a)
+    call mld%mold(altmp,info)
+#endif
   end if
 
   if (info /= psb_success_) then 
@@ -1265,6 +1271,7 @@ subroutine psb_c_cscnv_ip(a,info,type,mold,dupl)
 
 
   class(psb_c_base_sparse_mat), allocatable  :: altmp
+  class(psb_c_base_sparse_mat), pointer      :: mld
   integer(psb_ipk_) :: err_act
   character(len=20)  :: name='cscnv_ip'
   logical, parameter :: debug=.false.
@@ -1313,7 +1320,12 @@ subroutine psb_c_cscnv_ip(a,info,type,mold,dupl)
       goto 9999
     end select
   else
-    allocate(psb_c_csr_sparse_mat :: altmp, stat=info) 
+#if defined(HAVE_MOLD)
+    allocate(altmp, mold=psb_get_mat_default(a),stat=info) 
+#else
+    mld = psb_get_mat_default(a)
+    call mld%mold(altmp,info)
+#endif
   end if
 
   if (info /= psb_success_) then 
diff --git a/base/serial/impl/psb_d_coo_impl.f90 b/base/serial/impl/psb_d_coo_impl.f90
index 428cc9f8..91059b93 100644
--- a/base/serial/impl/psb_d_coo_impl.f90
+++ b/base/serial/impl/psb_d_coo_impl.f90
@@ -3375,7 +3375,7 @@ subroutine psb_d_fix_coo(a,info,idir)
   if (present(idir)) then 
     idir_ = idir
   else
-    idir_ = 0
+    idir_ = psb_row_major_
   endif
   
   nra = a%get_nrows()
diff --git a/base/serial/impl/psb_d_mat_impl.F90 b/base/serial/impl/psb_d_mat_impl.F90
index b675735a..629e1f04 100644
--- a/base/serial/impl/psb_d_mat_impl.F90
+++ b/base/serial/impl/psb_d_mat_impl.F90
@@ -1164,6 +1164,7 @@ subroutine psb_d_cscnv(a,b,info,type,mold,upd,dupl)
 
 
   class(psb_d_base_sparse_mat), allocatable  :: altmp
+  class(psb_d_base_sparse_mat), pointer      :: mld
   integer(psb_ipk_) :: err_act
   character(len=20)  :: name='cscnv'
   logical, parameter :: debug=.false.
@@ -1205,7 +1206,12 @@ subroutine psb_d_cscnv(a,b,info,type,mold,upd,dupl)
       goto 9999
     end select
   else
-    allocate(psb_d_csr_sparse_mat :: altmp, stat=info) 
+#if defined(HAVE_MOLD)
+    allocate(altmp, mold=psb_get_mat_default(a),stat=info) 
+#else
+    mld = psb_get_mat_default(a)
+    call mld%mold(altmp,info)
+#endif
   end if
 
   if (info /= psb_success_) then 
@@ -1265,6 +1271,7 @@ subroutine psb_d_cscnv_ip(a,info,type,mold,dupl)
 
 
   class(psb_d_base_sparse_mat), allocatable  :: altmp
+  class(psb_d_base_sparse_mat), pointer      :: mld
   integer(psb_ipk_) :: err_act
   character(len=20)  :: name='cscnv_ip'
   logical, parameter :: debug=.false.
@@ -1313,7 +1320,12 @@ subroutine psb_d_cscnv_ip(a,info,type,mold,dupl)
       goto 9999
     end select
   else
-    allocate(psb_d_csr_sparse_mat :: altmp, stat=info) 
+#if defined(HAVE_MOLD)
+    allocate(altmp, mold=psb_get_mat_default(a),stat=info) 
+#else
+    mld = psb_get_mat_default(a)
+    call mld%mold(altmp,info)
+#endif
   end if
 
   if (info /= psb_success_) then 
diff --git a/base/serial/impl/psb_s_coo_impl.f90 b/base/serial/impl/psb_s_coo_impl.f90
index 4f16f90e..7e166348 100644
--- a/base/serial/impl/psb_s_coo_impl.f90
+++ b/base/serial/impl/psb_s_coo_impl.f90
@@ -3375,7 +3375,7 @@ subroutine psb_s_fix_coo(a,info,idir)
   if (present(idir)) then 
     idir_ = idir
   else
-    idir_ = 0
+    idir_ = psb_row_major_
   endif
   
   nra = a%get_nrows()
diff --git a/base/serial/impl/psb_s_mat_impl.F90 b/base/serial/impl/psb_s_mat_impl.F90
index 67dd27e8..eaeec360 100644
--- a/base/serial/impl/psb_s_mat_impl.F90
+++ b/base/serial/impl/psb_s_mat_impl.F90
@@ -1164,6 +1164,7 @@ subroutine psb_s_cscnv(a,b,info,type,mold,upd,dupl)
 
 
   class(psb_s_base_sparse_mat), allocatable  :: altmp
+  class(psb_s_base_sparse_mat), pointer      :: mld
   integer(psb_ipk_) :: err_act
   character(len=20)  :: name='cscnv'
   logical, parameter :: debug=.false.
@@ -1205,7 +1206,12 @@ subroutine psb_s_cscnv(a,b,info,type,mold,upd,dupl)
       goto 9999
     end select
   else
-    allocate(psb_s_csr_sparse_mat :: altmp, stat=info) 
+#if defined(HAVE_MOLD)
+    allocate(altmp, mold=psb_get_mat_default(a),stat=info) 
+#else
+    mld = psb_get_mat_default(a)
+    call mld%mold(altmp,info)
+#endif
   end if
 
   if (info /= psb_success_) then 
@@ -1265,6 +1271,7 @@ subroutine psb_s_cscnv_ip(a,info,type,mold,dupl)
 
 
   class(psb_s_base_sparse_mat), allocatable  :: altmp
+  class(psb_s_base_sparse_mat), pointer      :: mld
   integer(psb_ipk_) :: err_act
   character(len=20)  :: name='cscnv_ip'
   logical, parameter :: debug=.false.
@@ -1313,7 +1320,12 @@ subroutine psb_s_cscnv_ip(a,info,type,mold,dupl)
       goto 9999
     end select
   else
-    allocate(psb_s_csr_sparse_mat :: altmp, stat=info) 
+#if defined(HAVE_MOLD)
+    allocate(altmp, mold=psb_get_mat_default(a),stat=info) 
+#else
+    mld = psb_get_mat_default(a)
+    call mld%mold(altmp,info)
+#endif
   end if
 
   if (info /= psb_success_) then 
diff --git a/base/serial/impl/psb_z_coo_impl.f90 b/base/serial/impl/psb_z_coo_impl.f90
index 2af0afce..d9927f8b 100644
--- a/base/serial/impl/psb_z_coo_impl.f90
+++ b/base/serial/impl/psb_z_coo_impl.f90
@@ -3375,7 +3375,7 @@ subroutine psb_z_fix_coo(a,info,idir)
   if (present(idir)) then 
     idir_ = idir
   else
-    idir_ = 0
+    idir_ = psb_row_major_
   endif
   
   nra = a%get_nrows()
diff --git a/base/serial/impl/psb_z_mat_impl.F90 b/base/serial/impl/psb_z_mat_impl.F90
index 83a757ca..65e2c121 100644
--- a/base/serial/impl/psb_z_mat_impl.F90
+++ b/base/serial/impl/psb_z_mat_impl.F90
@@ -1164,6 +1164,7 @@ subroutine psb_z_cscnv(a,b,info,type,mold,upd,dupl)
 
 
   class(psb_z_base_sparse_mat), allocatable  :: altmp
+  class(psb_z_base_sparse_mat), pointer      :: mld
   integer(psb_ipk_) :: err_act
   character(len=20)  :: name='cscnv'
   logical, parameter :: debug=.false.
@@ -1205,7 +1206,12 @@ subroutine psb_z_cscnv(a,b,info,type,mold,upd,dupl)
       goto 9999
     end select
   else
-    allocate(psb_z_csr_sparse_mat :: altmp, stat=info) 
+#if defined(HAVE_MOLD)
+    allocate(altmp, mold=psb_get_mat_default(a),stat=info) 
+#else
+    mld = psb_get_mat_default(a)
+    call mld%mold(altmp,info)
+#endif
   end if
 
   if (info /= psb_success_) then 
@@ -1265,6 +1271,7 @@ subroutine psb_z_cscnv_ip(a,info,type,mold,dupl)
 
 
   class(psb_z_base_sparse_mat), allocatable  :: altmp
+  class(psb_z_base_sparse_mat), pointer      :: mld
   integer(psb_ipk_) :: err_act
   character(len=20)  :: name='cscnv_ip'
   logical, parameter :: debug=.false.
@@ -1313,7 +1320,12 @@ subroutine psb_z_cscnv_ip(a,info,type,mold,dupl)
       goto 9999
     end select
   else
-    allocate(psb_z_csr_sparse_mat :: altmp, stat=info) 
+#if defined(HAVE_MOLD)
+    allocate(altmp, mold=psb_get_mat_default(a),stat=info) 
+#else
+    mld = psb_get_mat_default(a)
+    call mld%mold(altmp,info)
+#endif
   end if
 
   if (info /= psb_success_) then 
diff --git a/base/serial/psb_cspspmm.f90 b/base/serial/psb_cspspmm.f90
index 628cf129..8d39ea72 100644
--- a/base/serial/psb_cspspmm.f90
+++ b/base/serial/psb_cspspmm.f90
@@ -56,6 +56,10 @@ subroutine psb_cspspmm(a,b,c,info)
     goto 9999
   endif
 
+
+  ! 
+  ! Shortcuts for special cases
+  !
   done_spmm = .false. 
   select type(aa=>a%a)
   class is (psb_c_csr_sparse_mat) 
@@ -90,6 +94,9 @@ subroutine psb_cspspmm(a,b,c,info)
 
   end select
   
+  !
+  ! General code
+  !
   if (.not.done_spmm) then 
     call psb_symbmm(a,b,c,info)
     if (info == psb_success_) call psb_numbmm(a,b,c)
diff --git a/base/serial/psb_dspspmm.f90 b/base/serial/psb_dspspmm.f90
index d418ba6d..251a2d2d 100644
--- a/base/serial/psb_dspspmm.f90
+++ b/base/serial/psb_dspspmm.f90
@@ -56,6 +56,10 @@ subroutine psb_dspspmm(a,b,c,info)
     goto 9999
   endif
 
+
+  ! 
+  ! Shortcuts for special cases
+  !
   done_spmm = .false. 
   select type(aa=>a%a)
   class is (psb_d_csr_sparse_mat) 
@@ -90,6 +94,9 @@ subroutine psb_dspspmm(a,b,c,info)
 
   end select
   
+  !
+  ! General code
+  !
   if (.not.done_spmm) then 
     call psb_symbmm(a,b,c,info)
     if (info == psb_success_) call psb_numbmm(a,b,c)
diff --git a/base/serial/psb_sspspmm.f90 b/base/serial/psb_sspspmm.f90
index 55e79a7d..9c31bc17 100644
--- a/base/serial/psb_sspspmm.f90
+++ b/base/serial/psb_sspspmm.f90
@@ -56,6 +56,10 @@ subroutine psb_sspspmm(a,b,c,info)
     goto 9999
   endif
 
+
+  ! 
+  ! Shortcuts for special cases
+  !
   done_spmm = .false. 
   select type(aa=>a%a)
   class is (psb_s_csr_sparse_mat) 
@@ -90,6 +94,9 @@ subroutine psb_sspspmm(a,b,c,info)
 
   end select
   
+  !
+  ! General code
+  !
   if (.not.done_spmm) then 
     call psb_symbmm(a,b,c,info)
     if (info == psb_success_) call psb_numbmm(a,b,c)
diff --git a/base/serial/psb_zspspmm.f90 b/base/serial/psb_zspspmm.f90
index f3afe7fe..7b7fd1a1 100644
--- a/base/serial/psb_zspspmm.f90
+++ b/base/serial/psb_zspspmm.f90
@@ -56,6 +56,10 @@ subroutine psb_zspspmm(a,b,c,info)
     goto 9999
   endif
 
+
+  ! 
+  ! Shortcuts for special cases
+  !
   done_spmm = .false. 
   select type(aa=>a%a)
   class is (psb_z_csr_sparse_mat) 
@@ -90,6 +94,9 @@ subroutine psb_zspspmm(a,b,c,info)
 
   end select
   
+  !
+  ! General code
+  !
   if (.not.done_spmm) then 
     call psb_symbmm(a,b,c,info)
     if (info == psb_success_) call psb_numbmm(a,b,c)