diff --git a/base/modules/Makefile b/base/modules/Makefile
index 399add1f..af7b66d3 100644
--- a/base/modules/Makefile
+++ b/base/modules/Makefile
@@ -60,6 +60,7 @@ SERIAL_MODS=serial/psb_s_serial_mod.o serial/psb_d_serial_mod.o \
 	auxil/psb_d_hsort_x_mod.o \
 	auxil/psb_c_hsort_x_mod.o \
 	auxil/psb_z_hsort_x_mod.o \
+        serial/psb_z_csrli_mat_mod.o \
 	serial/psb_base_mat_mod.o   serial/psb_mat_mod.o\
 	serial/psb_s_base_mat_mod.o serial/psb_s_csr_mat_mod.o serial/psb_s_csc_mat_mod.o serial/psb_s_mat_mod.o \
 	serial/psb_d_base_mat_mod.o serial/psb_d_csr_mat_mod.o serial/psb_d_csc_mat_mod.o serial/psb_d_mat_mod.o \
@@ -241,8 +242,8 @@ serial/psb_s_base_mat_mod.o serial/psb_d_base_mat_mod.o serial/psb_c_base_mat_mo
 serial/psb_s_base_mat_mod.o: serial/psb_s_base_vect_mod.o
 serial/psb_d_base_mat_mod.o: serial/psb_d_base_vect_mod.o 
 #serial/psb_ld_base_mat_mod.o: serial/psb_d_base_vect_mod.o 
-serial/psb_c_base_mat_mod.o: serial/psb_c_base_vect_mod.o
-serial/psb_z_base_mat_mod.o: serial/psb_z_base_vect_mod.o
+serial/psb_c_base_mat_mod.o: serial/psb_c_base_vect_mod.o serial/psb_s_base_mat_mod.o
+serial/psb_z_base_mat_mod.o: serial/psb_z_base_vect_mod.o serial/psb_d_base_mat_mod.o
 serial/psb_l_base_vect_mod.o: serial/psb_i_base_vect_mod.o 
 serial/psb_c_base_vect_mod.o serial/psb_s_base_vect_mod.o serial/psb_d_base_vect_mod.o serial/psb_z_base_vect_mod.o: serial/psb_i_base_vect_mod.o serial/psb_l_base_vect_mod.o
 serial/psb_i_base_vect_mod.o serial/psb_l_base_vect_mod.o serial/psb_c_base_vect_mod.o serial/psb_s_base_vect_mod.o serial/psb_d_base_vect_mod.o serial/psb_z_base_vect_mod.o: auxil/psi_serial_mod.o psb_realloc_mod.o
@@ -253,11 +254,11 @@ serial/psb_d_mat_mod.o: serial/psb_d_base_mat_mod.o serial/psb_d_csr_mat_mod.o s
 serial/psb_c_mat_mod.o: serial/psb_c_base_mat_mod.o serial/psb_c_csr_mat_mod.o serial/psb_c_csc_mat_mod.o serial/psb_c_vect_mod.o \
 	 serial/psb_i_vect_mod.o serial/psb_l_vect_mod.o
 serial/psb_z_mat_mod.o: serial/psb_z_base_mat_mod.o serial/psb_z_csr_mat_mod.o serial/psb_z_csc_mat_mod.o serial/psb_z_vect_mod.o \
-	 serial/psb_i_vect_mod.o serial/psb_l_vect_mod.o
+	 serial/psb_i_vect_mod.o serial/psb_l_vect_mod.o serial/psb_z_csrli_mat_mod.o 
 serial/psb_s_csc_mat_mod.o serial/psb_s_csr_mat_mod.o serial/psb_ls_csr_mat_mod.o: serial/psb_s_base_mat_mod.o 
 serial/psb_d_csc_mat_mod.o serial/psb_d_csr_mat_mod.o serial/psb_ld_csr_mat_mod.o: serial/psb_d_base_mat_mod.o 
 serial/psb_c_csc_mat_mod.o serial/psb_c_csr_mat_mod.o serial/psb_lc_csr_mat_mod.o: serial/psb_c_base_mat_mod.o 
-serial/psb_z_csc_mat_mod.o serial/psb_z_csr_mat_mod.o serial/psb_lz_csr_mat_mod.o: serial/psb_z_base_mat_mod.o
+serial/psb_z_csrli_mat_mod.o serial/psb_z_csc_mat_mod.o serial/psb_z_csr_mat_mod.o serial/psb_lz_csr_mat_mod.o: serial/psb_z_base_mat_mod.o
 
 serial/psb_mat_mod.o: serial/psb_vect_mod.o serial/psb_s_mat_mod.o serial/psb_d_mat_mod.o serial/psb_c_mat_mod.o serial/psb_z_mat_mod.o
 serial/psb_serial_mod.o: serial/psb_s_serial_mod.o serial/psb_d_serial_mod.o serial/psb_c_serial_mod.o serial/psb_z_serial_mod.o  auxil/psi_serial_mod.o
diff --git a/base/modules/serial/psb_c_base_mat_mod.F90 b/base/modules/serial/psb_c_base_mat_mod.F90
index 7a7bbb1f..5e85e68a 100644
--- a/base/modules/serial/psb_c_base_mat_mod.F90
+++ b/base/modules/serial/psb_c_base_mat_mod.F90
@@ -35,7 +35,7 @@ module psb_c_base_mat_mod
 
   use psb_base_mat_mod
   use psb_c_base_vect_mod
-
+  use psb_s_base_mat_mod
 
   !> \namespace  psb_base_mod  \class  psb_c_base_sparse_mat
   !! \extends psb_base_mat_mod::psb_base_sparse_mat
@@ -90,6 +90,8 @@ module psb_c_base_mat_mod
     procedure, pass(a) :: mv_from_lcoo   => psb_c_base_mv_from_lcoo
     procedure, pass(a) :: mv_to_lfmt     => psb_c_base_mv_to_lfmt
     procedure, pass(a) :: mv_from_lfmt   => psb_c_base_mv_from_lfmt
+    procedure, pass(a) :: cp_to_real     => psb_c_base_cp_to_real
+    procedure, pass(a) :: cp_from_real   => psb_c_base_cp_from_real
 
 
     !
@@ -182,6 +184,8 @@ module psb_c_base_mat_mod
     !
     procedure, pass(a) :: cp_to_lcoo     => psb_c_cp_coo_to_lcoo
     procedure, pass(a) :: cp_from_lcoo   => psb_c_cp_coo_from_lcoo
+    procedure, pass(a) :: cp_to_coo_real   => psb_c_cp_coo_to_coo_real
+    procedure, pass(a) :: cp_from_coo_real => psb_c_cp_coo_from_coo_real
 
     procedure, pass(a) :: csput_a      => psb_c_coo_csput_a
     procedure, pass(a) :: get_diag     => psb_c_coo_get_diag
@@ -1069,7 +1073,22 @@ module psb_c_base_mat_mod
       integer(psb_ipk_), intent(out)            :: info
     end subroutine psb_c_base_mv_from_lfmt
   end interface
-
+  interface
+    subroutine psb_c_base_cp_to_real(a,b,info)
+      import
+      class(psb_c_base_sparse_mat), intent(inout) :: a
+      class(psb_s_base_sparse_mat), intent(inout) :: b
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine psb_c_base_cp_to_real
+  end interface
+  interface
+    subroutine psb_c_base_cp_from_real(a,b,info)
+      import
+      class(psb_c_base_sparse_mat), intent(inout) :: a
+      class(psb_s_base_sparse_mat), intent(inout) :: b
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine psb_c_base_cp_from_real
+  end interface
 
   !
   !>
@@ -1954,6 +1973,23 @@ module psb_c_base_mat_mod
     end subroutine psb_c_cp_coo_from_lcoo
   end interface
 
+  interface
+    subroutine psb_c_cp_coo_to_coo_real(a,b,info)
+      import
+      class(psb_c_coo_sparse_mat), intent(inout) :: a
+      class(psb_s_coo_sparse_mat), intent(inout) :: b
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine psb_c_cp_coo_to_coo_real
+  end interface
+  interface
+    subroutine psb_c_cp_coo_from_coo_real(a,b,info)
+      import
+      class(psb_c_coo_sparse_mat), intent(inout) :: a
+      class(psb_s_coo_sparse_mat), intent(inout) :: b
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine psb_c_cp_coo_from_coo_real
+  end interface
+  
   !>
   !! \memberof  psb_c_coo_sparse_mat
   !! \see psb_c_base_mat_mod::psb_c_base_cp_from_coo
diff --git a/base/modules/serial/psb_d_base_mat_mod.F90 b/base/modules/serial/psb_d_base_mat_mod.F90
index eb49905d..3ef455d7 100644
--- a/base/modules/serial/psb_d_base_mat_mod.F90
+++ b/base/modules/serial/psb_d_base_mat_mod.F90
@@ -36,7 +36,6 @@ module psb_d_base_mat_mod
   use psb_base_mat_mod
   use psb_d_base_vect_mod
 
-
   !> \namespace  psb_base_mod  \class  psb_d_base_sparse_mat
   !! \extends psb_base_mat_mod::psb_base_sparse_mat
   !! The psb_d_base_sparse_mat type, extending psb_base_sparse_mat,
@@ -1070,7 +1069,6 @@ module psb_d_base_mat_mod
     end subroutine psb_d_base_mv_from_lfmt
   end interface
 
-
   !
   !>
   !! \memberof  psb_d_base_sparse_mat
@@ -1954,6 +1952,7 @@ module psb_d_base_mat_mod
     end subroutine psb_d_cp_coo_from_lcoo
   end interface
 
+  
   !>
   !! \memberof  psb_d_coo_sparse_mat
   !! \see psb_d_base_mat_mod::psb_d_base_cp_from_coo
diff --git a/base/modules/serial/psb_s_base_mat_mod.F90 b/base/modules/serial/psb_s_base_mat_mod.F90
index 79c8222b..30c84a23 100644
--- a/base/modules/serial/psb_s_base_mat_mod.F90
+++ b/base/modules/serial/psb_s_base_mat_mod.F90
@@ -36,7 +36,6 @@ module psb_s_base_mat_mod
   use psb_base_mat_mod
   use psb_s_base_vect_mod
 
-
   !> \namespace  psb_base_mod  \class  psb_s_base_sparse_mat
   !! \extends psb_base_mat_mod::psb_base_sparse_mat
   !! The psb_s_base_sparse_mat type, extending psb_base_sparse_mat,
@@ -1070,7 +1069,6 @@ module psb_s_base_mat_mod
     end subroutine psb_s_base_mv_from_lfmt
   end interface
 
-
   !
   !>
   !! \memberof  psb_s_base_sparse_mat
@@ -1954,6 +1952,7 @@ module psb_s_base_mat_mod
     end subroutine psb_s_cp_coo_from_lcoo
   end interface
 
+  
   !>
   !! \memberof  psb_s_coo_sparse_mat
   !! \see psb_s_base_mat_mod::psb_s_base_cp_from_coo
diff --git a/base/modules/serial/psb_z_base_mat_mod.F90 b/base/modules/serial/psb_z_base_mat_mod.F90
index 5b6ca07b..066f5926 100644
--- a/base/modules/serial/psb_z_base_mat_mod.F90
+++ b/base/modules/serial/psb_z_base_mat_mod.F90
@@ -35,7 +35,7 @@ module psb_z_base_mat_mod
 
   use psb_base_mat_mod
   use psb_z_base_vect_mod
-
+  use psb_d_base_mat_mod
 
   !> \namespace  psb_base_mod  \class  psb_z_base_sparse_mat
   !! \extends psb_base_mat_mod::psb_base_sparse_mat
@@ -90,6 +90,8 @@ module psb_z_base_mat_mod
     procedure, pass(a) :: mv_from_lcoo   => psb_z_base_mv_from_lcoo
     procedure, pass(a) :: mv_to_lfmt     => psb_z_base_mv_to_lfmt
     procedure, pass(a) :: mv_from_lfmt   => psb_z_base_mv_from_lfmt
+    procedure, pass(a) :: cp_to_real     => psb_z_base_cp_to_real
+    procedure, pass(a) :: cp_from_real   => psb_z_base_cp_from_real
 
 
     !
@@ -182,6 +184,8 @@ module psb_z_base_mat_mod
     !
     procedure, pass(a) :: cp_to_lcoo     => psb_z_cp_coo_to_lcoo
     procedure, pass(a) :: cp_from_lcoo   => psb_z_cp_coo_from_lcoo
+    procedure, pass(a) :: cp_to_coo_real   => psb_z_cp_coo_to_coo_real
+    procedure, pass(a) :: cp_from_coo_real => psb_z_cp_coo_from_coo_real
 
     procedure, pass(a) :: csput_a      => psb_z_coo_csput_a
     procedure, pass(a) :: get_diag     => psb_z_coo_get_diag
@@ -1069,7 +1073,22 @@ module psb_z_base_mat_mod
       integer(psb_ipk_), intent(out)            :: info
     end subroutine psb_z_base_mv_from_lfmt
   end interface
-
+  interface
+    subroutine psb_z_base_cp_to_real(a,b,info)
+      import
+      class(psb_z_base_sparse_mat), intent(inout) :: a
+      class(psb_d_base_sparse_mat), intent(inout) :: b
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine psb_z_base_cp_to_real
+  end interface
+  interface
+    subroutine psb_z_base_cp_from_real(a,b,info)
+      import
+      class(psb_z_base_sparse_mat), intent(inout) :: a
+      class(psb_d_base_sparse_mat), intent(inout) :: b
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine psb_z_base_cp_from_real
+  end interface
 
   !
   !>
@@ -1954,6 +1973,23 @@ module psb_z_base_mat_mod
     end subroutine psb_z_cp_coo_from_lcoo
   end interface
 
+  interface
+    subroutine psb_z_cp_coo_to_coo_real(a,b,info)
+      import
+      class(psb_z_coo_sparse_mat), intent(inout) :: a
+      class(psb_d_coo_sparse_mat), intent(inout) :: b
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine psb_z_cp_coo_to_coo_real
+  end interface
+  interface
+    subroutine psb_z_cp_coo_from_coo_real(a,b,info)
+      import
+      class(psb_z_coo_sparse_mat), intent(inout) :: a
+      class(psb_d_coo_sparse_mat), intent(inout) :: b
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine psb_z_cp_coo_from_coo_real
+  end interface
+  
   !>
   !! \memberof  psb_z_coo_sparse_mat
   !! \see psb_z_base_mat_mod::psb_z_base_cp_from_coo
diff --git a/base/modules/serial/psb_z_csrli_mat_mod.f90 b/base/modules/serial/psb_z_csrli_mat_mod.f90
index 17b69f26..84fd6ac8 100644
--- a/base/modules/serial/psb_z_csrli_mat_mod.f90
+++ b/base/modules/serial/psb_z_csrli_mat_mod.f90
@@ -48,7 +48,7 @@ module psb_z_csrli_mat_mod
 
   type, extends(psb_z_csr_sparse_mat) :: psb_z_csrli_sparse_mat
 
-    complex(psb_dpk_) :: lambda
+    complex(psb_dpk_) :: lambda=zzero
     
   contains
     procedure, nopass  :: get_fmt     => z_csrli_get_fmt
@@ -57,15 +57,13 @@ module psb_z_csrli_mat_mod
     procedure, pass(a) :: inner_cssm  => psb_z_csrli_cssm
     procedure, pass(a) :: inner_cssv  => psb_z_csrli_cssv
     procedure, pass(a) :: scals       => psb_z_csrli_scals
-    !procedure, pass(a) :: scalv       => psb_z_csrli_scal
+    procedure, pass(a) :: scalv       => psb_z_csrli_scal
     procedure, pass(a) :: maxval      => psb_z_csrli_maxval
     procedure, pass(a) :: spnmi       => psb_z_csrli_csnmi
     procedure, pass(a) :: rowsum      => psb_z_csrli_rowsum
     procedure, pass(a) :: arwsum      => psb_z_csrli_arwsum
     procedure, pass(a) :: colsum      => psb_z_csrli_colsum
     procedure, pass(a) :: aclsum      => psb_z_csrli_aclsum
-!!$    procedure, pass(a) :: reallocate_nz => psb_z_csrli_reallocate_nz
-!!$    procedure, pass(a) :: allocate_mnnz => psb_z_csrli_allocate_mnnz
     procedure, pass(a) :: tril          => psb_z_csrli_tril
     procedure, pass(a) :: triu          => psb_z_csrli_triu
     procedure, pass(a) :: cp_to_coo   => psb_z_cp_csrli_to_coo
@@ -83,6 +81,9 @@ module psb_z_csrli_mat_mod
     procedure, pass(a) :: free        => z_csrli_free
     procedure, pass(a) :: mold        => psb_z_csrli_mold
 
+    procedure, pass(a) :: set_lambda  => z_csrli_set_lambda
+    procedure, pass(a) :: get_lambda  => z_csrli_get_lambda
+    
   end type psb_z_csrli_sparse_mat
 
   private :: z_csrli_get_nzeros, z_csrli_free,  z_csrli_get_fmt, &
@@ -109,14 +110,14 @@ module psb_z_csrli_mat_mod
     end subroutine psb_z_csrli_reinit
   end interface
 
-  !> \memberof psb_z_csrli_sparse_mat
-  !| \see psb_base_mat_mod::psb_base_trim
-  interface
-    subroutine  psb_z_csrli_trim(a)
-      import
-      class(psb_z_csrli_sparse_mat), intent(inout) :: a
-    end subroutine psb_z_csrli_trim
-  end interface
+!!$  !> \memberof psb_z_csrli_sparse_mat
+!!$  !| \see psb_base_mat_mod::psb_base_trim
+!!$  interface
+!!$    subroutine  psb_z_csrli_trim(a)
+!!$      import
+!!$      class(psb_z_csrli_sparse_mat), intent(inout) :: a
+!!$    end subroutine psb_z_csrli_trim
+!!$  end interface
 
 
   !> \memberof psb_z_csrli_sparse_mat
@@ -246,14 +247,14 @@ module psb_z_csrli_mat_mod
   !! \memberof  psb_z_csrli_sparse_mat
   !! \see psb_z_base_mat_mod::psb_z_base_clean_zeros
   !
-  interface
-    subroutine  psb_z_csrli_clean_zeros(a, info)
-      import
-      class(psb_z_csrli_sparse_mat), intent(inout) :: a
-      integer(psb_ipk_), intent(out)              :: info
-    end subroutine psb_z_csrli_clean_zeros
-  end interface
-
+!!$  interface
+!!$    subroutine  psb_z_csrli_clean_zeros(a, info)
+!!$      import
+!!$      class(psb_z_csrli_sparse_mat), intent(inout) :: a
+!!$      integer(psb_ipk_), intent(out)              :: info
+!!$    end subroutine psb_z_csrli_clean_zeros
+!!$  end interface
+!!$
   !> \memberof psb_z_csrli_sparse_mat
   !! \see psb_z_base_mat_mod::psb_z_base_cp_to_coo
   interface
@@ -363,36 +364,36 @@ module psb_z_csrli_mat_mod
   end interface
 
 
-  !> \memberof psb_z_csrli_sparse_mat
-  !! \see psb_z_base_mat_mod::psb_z_base_csput_a
-  interface
-    subroutine psb_z_csrli_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info)
-      import
-      class(psb_z_csrli_sparse_mat), intent(inout) :: a
-      complex(psb_dpk_), intent(in)      :: val(:)
-      integer(psb_ipk_), intent(in)             :: nz,ia(:), ja(:),&
-           &  imin,imax,jmin,jmax
-      integer(psb_ipk_), intent(out)            :: info
-    end subroutine psb_z_csrli_csput_a
-  end interface
-
-  !> \memberof psb_z_csrli_sparse_mat
-  !! \see psb_base_mat_mod::psb_base_csgetptn
-  interface
-    subroutine psb_z_csrli_csgetptn(imin,imax,a,nz,ia,ja,info,&
-         & jmin,jmax,iren,append,nzin,rscale,cscale)
-      import
-      class(psb_z_csrli_sparse_mat), intent(in) :: a
-      integer(psb_ipk_), intent(in)                  :: imin,imax
-      integer(psb_ipk_), intent(out)                 :: nz
-      integer(psb_ipk_), allocatable, intent(inout)  :: ia(:), ja(:)
-      integer(psb_ipk_),intent(out)                  :: info
-      logical, intent(in), optional        :: append
-      integer(psb_ipk_), intent(in), optional        :: iren(:)
-      integer(psb_ipk_), intent(in), optional        :: jmin,jmax, nzin
-      logical, intent(in), optional        :: rscale,cscale
-    end subroutine psb_z_csrli_csgetptn
-  end interface
+!!$  !> \memberof psb_z_csrli_sparse_mat
+!!$  !! \see psb_z_base_mat_mod::psb_z_base_csput_a
+!!$  interface
+!!$    subroutine psb_z_csrli_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info)
+!!$      import
+!!$      class(psb_z_csrli_sparse_mat), intent(inout) :: a
+!!$      complex(psb_dpk_), intent(in)      :: val(:)
+!!$      integer(psb_ipk_), intent(in)             :: nz,ia(:), ja(:),&
+!!$           &  imin,imax,jmin,jmax
+!!$      integer(psb_ipk_), intent(out)            :: info
+!!$    end subroutine psb_z_csrli_csput_a
+!!$  end interface
+!!$
+!!$  !> \memberof psb_z_csrli_sparse_mat
+!!$  !! \see psb_base_mat_mod::psb_base_csgetptn
+!!$  interface
+!!$    subroutine psb_z_csrli_csgetptn(imin,imax,a,nz,ia,ja,info,&
+!!$         & jmin,jmax,iren,append,nzin,rscale,cscale)
+!!$      import
+!!$      class(psb_z_csrli_sparse_mat), intent(in) :: a
+!!$      integer(psb_ipk_), intent(in)                  :: imin,imax
+!!$      integer(psb_ipk_), intent(out)                 :: nz
+!!$      integer(psb_ipk_), allocatable, intent(inout)  :: ia(:), ja(:)
+!!$      integer(psb_ipk_),intent(out)                  :: info
+!!$      logical, intent(in), optional        :: append
+!!$      integer(psb_ipk_), intent(in), optional        :: iren(:)
+!!$      integer(psb_ipk_), intent(in), optional        :: jmin,jmax, nzin
+!!$      logical, intent(in), optional        :: rscale,cscale
+!!$    end subroutine psb_z_csrli_csgetptn
+!!$  end interface
 
   !> \memberof psb_z_csrli_sparse_mat
   !! \see psb_z_base_mat_mod::psb_z_base_csgetrow
@@ -580,6 +581,21 @@ contains
     res = 'CSRLI'
   end function z_csrli_get_fmt
 
+
+  function z_csrli_get_lambda(a) result(res)
+    implicit none
+    class(psb_z_csrli_sparse_mat), intent(in) :: a
+    complex(psb_dpk_) :: res
+    res = a%lambda
+  end function z_csrli_get_lambda
+
+  subroutine z_csrli_set_lambda(a,val) 
+    implicit none
+    class(psb_z_csrli_sparse_mat), intent(inout) :: a
+    complex(psb_dpk_), intent(in) :: val
+    a%lambda = val
+  end subroutine z_csrli_set_lambda
+
   ! == ===================================
   !
   !
diff --git a/base/modules/serial/psb_z_mat_mod.F90 b/base/modules/serial/psb_z_mat_mod.F90
index c534cad5..fa0358be 100644
--- a/base/modules/serial/psb_z_mat_mod.F90
+++ b/base/modules/serial/psb_z_mat_mod.F90
@@ -81,7 +81,8 @@ module psb_z_mat_mod
   use psb_z_base_mat_mod
   use psb_z_csr_mat_mod,  only : psb_z_csr_sparse_mat, psb_lz_csr_sparse_mat
   use psb_z_csc_mat_mod,  only : psb_z_csc_sparse_mat, psb_lz_csc_sparse_mat
-
+  use psb_z_csrli_mat_mod
+  
   type :: psb_zspmat_type
 
     class(psb_z_base_sparse_mat), allocatable  :: a   
diff --git a/base/serial/impl/Makefile b/base/serial/impl/Makefile
index 00088741..87c26ebe 100644
--- a/base/serial/impl/Makefile
+++ b/base/serial/impl/Makefile
@@ -16,7 +16,7 @@ DOBJS=psb_d_csr_impl.o psb_d_coo_impl.o psb_d_csc_impl.o psb_d_mat_impl.o
 COBJS=psb_c_csr_impl.o psb_c_coo_impl.o psb_c_csc_impl.o psb_c_mat_impl.o
 #\
  psb_c_lcoo_impl.o  psb_c_lcsr_impl.o 
-ZOBJS=psb_z_csr_impl.o psb_z_coo_impl.o psb_z_csc_impl.o psb_z_mat_impl.o
+ZOBJS=psb_z_csr_impl.o psb_z_coo_impl.o psb_z_csc_impl.o psb_z_mat_impl.o psb_z_csrli_impl.o
 #\
  psb_z_lcoo_impl.o  psb_z_lcsr_impl.o 
 
diff --git a/base/serial/impl/psb_c_base_mat_impl.F90 b/base/serial/impl/psb_c_base_mat_impl.F90
index 17f2cdc8..35acb595 100644
--- a/base/serial/impl/psb_c_base_mat_impl.F90
+++ b/base/serial/impl/psb_c_base_mat_impl.F90
@@ -2507,6 +2507,81 @@ subroutine psb_c_base_mv_from_lfmt(a,b,info)
 
 end subroutine psb_c_base_mv_from_lfmt
 
+subroutine psb_c_base_cp_to_real(a,b,info)
+  use psb_c_base_mat_mod, psb_protect_name => psb_c_base_cp_to_real
+  use psb_s_base_mat_mod
+  use psb_error_mod
+  use psb_realloc_mod
+  implicit none
+  
+  class(psb_c_base_sparse_mat), intent(inout) :: a
+  class(psb_s_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_), intent(out)            :: info
+
+  !
+  integer(psb_ipk_)  :: err_act
+  character(len=20)  :: name='to_real'
+  logical, parameter :: debug=.false.
+  
+  type(psb_s_coo_sparse_mat) ::  rtemp
+  type(psb_c_coo_sparse_mat)  ::  ctemp
+  
+  !
+  ! Default implementation
+  !
+  info  = psb_success_
+  call psb_erractionsave(err_act)
+  call a%sync()
+  call a%cp_to_coo(ctemp,info)
+  call ctemp%cp_to_coo_real(rtemp,info)
+  call ctemp%free()
+  call rtemp%cp_to_fmt(b,info)
+  call rtemp%free()
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end subroutine psb_c_base_cp_to_real
+
+
+subroutine psb_c_base_cp_from_real(a,b,info)
+  use psb_c_base_mat_mod, psb_protect_name => psb_c_base_cp_from_real
+  use psb_s_base_mat_mod
+  use psb_error_mod
+  use psb_realloc_mod
+  implicit none
+  class(psb_c_base_sparse_mat), intent(inout) :: a
+  class(psb_s_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_), intent(out)            :: info
+  !
+  integer(psb_ipk_)  :: err_act
+  character(len=20)  :: name='to_real'
+  logical, parameter :: debug=.false.
+  
+  type(psb_s_coo_sparse_mat) ::  rtemp
+  type(psb_c_coo_sparse_mat)  ::  ctemp
+  
+  !
+  ! Default implementation
+  !
+  info  = psb_success_
+  call psb_erractionsave(err_act)
+
+  call rtemp%cp_from_fmt(b,info)
+  call ctemp%cp_from_coo_real(rtemp,info)
+  call rtemp%free()
+  call a%mv_from_coo(ctemp,info)
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end subroutine psb_c_base_cp_from_real
+
 !
 !
 ! lc implementation
diff --git a/base/serial/impl/psb_c_coo_impl.F90 b/base/serial/impl/psb_c_coo_impl.F90
index 306aefd9..65f9f2e7 100644
--- a/base/serial/impl/psb_c_coo_impl.F90
+++ b/base/serial/impl/psb_c_coo_impl.F90
@@ -5114,6 +5114,92 @@ subroutine psb_c_cp_coo_from_lcoo(a,b,info)
 
 end subroutine psb_c_cp_coo_from_lcoo
 
+subroutine psb_c_cp_coo_to_coo_real(a,b,info)
+  use psb_error_mod
+  use psb_c_base_mat_mod, psb_protect_name => psb_c_cp_coo_to_coo_real
+  implicit none
+
+  class(psb_c_coo_sparse_mat), intent(inout) :: a
+  class(psb_s_coo_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_), intent(out)            :: info
+
+
+  integer(psb_ipk_)  :: err_act
+  character(len=20)  :: name='to_coo_real'
+  logical, parameter :: debug=.false.
+  integer(psb_ipk_) :: m,n,nz
+
+  call psb_erractionsave(err_act)
+  info = psb_success_
+
+  call a%sync()
+  m  = a%get_nrows()
+  n  = a%get_ncols()
+  nz = a%get_nzeros()
+  call b%set_nrows(m)
+  call b%set_ncols(n)
+  call b%reallocate(nz)
+  b%ia(1:nz)  = a%ia(1:nz)
+  b%ja(1:nz)  = a%ja(1:nz)
+  b%val(1:nz)  = real(a%val(1:nz))
+  call b%set_nzeros(nz)
+  call b%set_sorted(a%is_sorted())
+  call b%set_host()
+  
+  call psb_erractionrestore(err_act)
+  return
+
+9999 continue
+  call psb_errpush(info,name)
+
+  call psb_error_handler(err_act)
+
+  return
+
+end subroutine psb_c_cp_coo_to_coo_real
+
+subroutine psb_c_cp_coo_from_coo_real(a,b,info)
+  use psb_error_mod
+  use psb_c_base_mat_mod, psb_protect_name => psb_c_cp_coo_from_coo_real
+  implicit none
+  class(psb_c_coo_sparse_mat), intent(inout) :: a
+  class(psb_s_coo_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_), intent(out)            :: info
+
+
+  integer(psb_ipk_)  :: err_act
+  character(len=20)  :: name='from_coo_real'
+  logical, parameter :: debug=.false.
+  integer(psb_ipk_) :: m,n,nz
+
+  call psb_erractionsave(err_act)
+  info = psb_success_
+
+  call b%sync()
+  m  = b%get_nrows()
+  n  = b%get_ncols()
+  nz = b%get_nzeros()
+  call a%set_nrows(m)
+  call a%set_ncols(n)
+  call a%reallocate(nz)
+  a%ia(1:nz)  = b%ia(1:nz)
+  a%ja(1:nz)  = b%ja(1:nz)
+  a%val(1:nz)  = b%val(1:nz)
+  call a%set_nzeros(nz)
+  call a%set_sorted(b%is_sorted())
+  call a%set_host()
+  
+  call psb_erractionrestore(err_act)
+  return
+
+9999 continue
+  call psb_errpush(info,name)
+
+  call psb_error_handler(err_act)
+
+  return
+
+end subroutine psb_c_cp_coo_from_coo_real
 
 !
 !
diff --git a/base/serial/impl/psb_z_base_mat_impl.F90 b/base/serial/impl/psb_z_base_mat_impl.F90
index 404027c5..ea2b66a9 100644
--- a/base/serial/impl/psb_z_base_mat_impl.F90
+++ b/base/serial/impl/psb_z_base_mat_impl.F90
@@ -2507,6 +2507,81 @@ subroutine psb_z_base_mv_from_lfmt(a,b,info)
 
 end subroutine psb_z_base_mv_from_lfmt
 
+subroutine psb_z_base_cp_to_real(a,b,info)
+  use psb_z_base_mat_mod, psb_protect_name => psb_z_base_cp_to_real
+  use psb_d_base_mat_mod
+  use psb_error_mod
+  use psb_realloc_mod
+  implicit none
+  
+  class(psb_z_base_sparse_mat), intent(inout) :: a
+  class(psb_d_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_), intent(out)            :: info
+
+  !
+  integer(psb_ipk_)  :: err_act
+  character(len=20)  :: name='to_real'
+  logical, parameter :: debug=.false.
+  
+  type(psb_d_coo_sparse_mat) ::  rtemp
+  type(psb_z_coo_sparse_mat)  ::  ctemp
+  
+  !
+  ! Default implementation
+  !
+  info  = psb_success_
+  call psb_erractionsave(err_act)
+  call a%sync()
+  call a%cp_to_coo(ctemp,info)
+  call ctemp%cp_to_coo_real(rtemp,info)
+  call ctemp%free()
+  call rtemp%cp_to_fmt(b,info)
+  call rtemp%free()
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end subroutine psb_z_base_cp_to_real
+
+
+subroutine psb_z_base_cp_from_real(a,b,info)
+  use psb_z_base_mat_mod, psb_protect_name => psb_z_base_cp_from_real
+  use psb_d_base_mat_mod
+  use psb_error_mod
+  use psb_realloc_mod
+  implicit none
+  class(psb_z_base_sparse_mat), intent(inout) :: a
+  class(psb_d_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_), intent(out)            :: info
+  !
+  integer(psb_ipk_)  :: err_act
+  character(len=20)  :: name='to_real'
+  logical, parameter :: debug=.false.
+  
+  type(psb_d_coo_sparse_mat) ::  rtemp
+  type(psb_z_coo_sparse_mat)  ::  ctemp
+  
+  !
+  ! Default implementation
+  !
+  info  = psb_success_
+  call psb_erractionsave(err_act)
+
+  call rtemp%cp_from_fmt(b,info)
+  call ctemp%cp_from_coo_real(rtemp,info)
+  call rtemp%free()
+  call a%mv_from_coo(ctemp,info)
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end subroutine psb_z_base_cp_from_real
+
 !
 !
 ! lz implementation
diff --git a/base/serial/impl/psb_z_coo_impl.F90 b/base/serial/impl/psb_z_coo_impl.F90
index b9dd7e59..1af133c8 100644
--- a/base/serial/impl/psb_z_coo_impl.F90
+++ b/base/serial/impl/psb_z_coo_impl.F90
@@ -5114,6 +5114,92 @@ subroutine psb_z_cp_coo_from_lcoo(a,b,info)
 
 end subroutine psb_z_cp_coo_from_lcoo
 
+subroutine psb_z_cp_coo_to_coo_real(a,b,info)
+  use psb_error_mod
+  use psb_z_base_mat_mod, psb_protect_name => psb_z_cp_coo_to_coo_real
+  implicit none
+
+  class(psb_z_coo_sparse_mat), intent(inout) :: a
+  class(psb_d_coo_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_), intent(out)            :: info
+
+
+  integer(psb_ipk_)  :: err_act
+  character(len=20)  :: name='to_coo_real'
+  logical, parameter :: debug=.false.
+  integer(psb_ipk_) :: m,n,nz
+
+  call psb_erractionsave(err_act)
+  info = psb_success_
+
+  call a%sync()
+  m  = a%get_nrows()
+  n  = a%get_ncols()
+  nz = a%get_nzeros()
+  call b%set_nrows(m)
+  call b%set_ncols(n)
+  call b%reallocate(nz)
+  b%ia(1:nz)  = a%ia(1:nz)
+  b%ja(1:nz)  = a%ja(1:nz)
+  b%val(1:nz)  = real(a%val(1:nz))
+  call b%set_nzeros(nz)
+  call b%set_sorted(a%is_sorted())
+  call b%set_host()
+  
+  call psb_erractionrestore(err_act)
+  return
+
+9999 continue
+  call psb_errpush(info,name)
+
+  call psb_error_handler(err_act)
+
+  return
+
+end subroutine psb_z_cp_coo_to_coo_real
+
+subroutine psb_z_cp_coo_from_coo_real(a,b,info)
+  use psb_error_mod
+  use psb_z_base_mat_mod, psb_protect_name => psb_z_cp_coo_from_coo_real
+  implicit none
+  class(psb_z_coo_sparse_mat), intent(inout) :: a
+  class(psb_d_coo_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_), intent(out)            :: info
+
+
+  integer(psb_ipk_)  :: err_act
+  character(len=20)  :: name='from_coo_real'
+  logical, parameter :: debug=.false.
+  integer(psb_ipk_) :: m,n,nz
+
+  call psb_erractionsave(err_act)
+  info = psb_success_
+
+  call b%sync()
+  m  = b%get_nrows()
+  n  = b%get_ncols()
+  nz = b%get_nzeros()
+  call a%set_nrows(m)
+  call a%set_ncols(n)
+  call a%reallocate(nz)
+  a%ia(1:nz)  = b%ia(1:nz)
+  a%ja(1:nz)  = b%ja(1:nz)
+  a%val(1:nz)  = b%val(1:nz)
+  call a%set_nzeros(nz)
+  call a%set_sorted(b%is_sorted())
+  call a%set_host()
+  
+  call psb_erractionrestore(err_act)
+  return
+
+9999 continue
+  call psb_errpush(info,name)
+
+  call psb_error_handler(err_act)
+
+  return
+
+end subroutine psb_z_cp_coo_from_coo_real
 
 !
 !
diff --git a/test/pargen/Makefile b/test/pargen/Makefile
index 20a95c0b..07943a43 100644
--- a/test/pargen/Makefile
+++ b/test/pargen/Makefile
@@ -16,11 +16,15 @@ FINCLUDES=$(FMFLAG)$(MODDIR) $(FMFLAG).
 
 EXEDIR=./runs
 
-all: runsd psb_d_pde3d psb_s_pde3d psb_d_pde2d psb_s_pde2d
+all: runsd psb_tzcsrli psb_d_pde3d psb_s_pde3d psb_d_pde2d psb_s_pde2d
 
 runsd:
 	(if test ! -d runs ; then mkdir runs; fi)
 
+psb_tzcsrli: psb_tzcsrli.o 
+	$(FLINK)  psb_tzcsrli.o -o psb_tzcsrli $(PSBLAS_LIB) $(LDLIBS)
+	/bin/mv psb_tzcsrli $(EXEDIR)
+
 psb_d_pde3d: psb_d_pde3d.o 
 	$(FLINK)  psb_d_pde3d.o -o psb_d_pde3d $(PSBLAS_LIB) $(LDLIBS)
 	/bin/mv psb_d_pde3d $(EXEDIR)
@@ -41,7 +45,7 @@ psb_s_pde2d: psb_s_pde2d.o
 
 
 clean: 
-	/bin/rm -f psb_d_pde3d.o psb_s_pde3d.o  psb_d_pde2d.o psb_s_pde2d.o  *$(.mod) \
+	/bin/rm -f psb_tzcsrli.o psb_d_pde3d.o psb_s_pde3d.o  psb_d_pde2d.o psb_s_pde2d.o  *$(.mod) \
 	$(EXEDIR)/psb_d_pde3d $(EXEDIR)/psb_s_pde3d $(EXEDIR)/psb_d_pde2d $(EXEDIR)/psb_s_pde2d 
 verycleanlib: 
 	(cd ../..; make veryclean)
diff --git a/test/pargen/psb_tzcsrli.F90 b/test/pargen/psb_tzcsrli.F90
new file mode 100644
index 00000000..d2697623
--- /dev/null
+++ b/test/pargen/psb_tzcsrli.F90
@@ -0,0 +1,1139 @@
+!
+!                Parallel Sparse BLAS  version 3.5
+!      (C) Copyright 2006-2018
+!        Salvatore Filippone
+!        Alfredo Buttari
+!
+!    Redistribution and use in source and binary forms, with or without
+!    modification, are permitted provided that the following conditions
+!    are met:
+!      1. Redistributions of source code must retain the above copyright
+!         notice, this list of conditions and the following disclaimer.
+!      2. Redistributions in binary form must reproduce the above copyright
+!         notice, this list of conditions, and the following disclaimer in the
+!         documentation and/or other materials provided with the distribution.
+!      3. The name of the PSBLAS group or the names of its contributors may
+!         not be used to endorse or promote products derived from this
+!         software without specific written permission.
+!
+!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
+!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+!    POSSIBILITY OF SUCH DAMAGE.
+!
+!
+! File: psb_d_pde3d.f90
+!
+! Program: psb_d_pde3d
+! This sample program solves a linear system obtained by discretizing a
+! PDE with Dirichlet BCs.
+!
+!
+! The PDE is a general second order equation in 3d
+!
+!   a1 dd(u)  a2 dd(u)    a3 dd(u)    b1 d(u)   b2 d(u)  b3 d(u)
+! -   ------ -  ------ -  ------ +  -----  +  ------  +  ------ + c u = f
+!      dxdx     dydy       dzdz        dx       dy         dz
+!
+! with Dirichlet boundary conditions
+!   u = g
+!
+!  on the unit cube  0<=x,y,z<=1.
+!
+!
+! Note that if b1=b2=b3=c=0., the PDE is the  Laplace equation.
+!
+! There are three choices available for data distribution:
+! 1. A simple BLOCK distribution
+! 2. A ditribution based on arbitrary assignment of indices to processes,
+!    typically from a graph partitioner
+! 3. A 3D distribution in which the unit cube is partitioned
+!    into subcubes, each one assigned to a process.
+!
+!
+module psb_d_pde3d_mod
+
+
+  use psb_base_mod, only : psb_dpk_, psb_ipk_, psb_lpk_, psb_desc_type,&
+       &  psb_dspmat_type, psb_d_vect_type, dzero,&
+       &  psb_d_base_sparse_mat, psb_d_base_vect_type, &
+       &  psb_i_base_vect_type, psb_l_base_vect_type
+
+  interface
+    function d_func_3d(x,y,z) result(val)
+      import :: psb_dpk_
+      real(psb_dpk_), intent(in) :: x,y,z
+      real(psb_dpk_) :: val
+    end function d_func_3d
+  end interface
+
+  interface psb_gen_pde3d
+    module procedure  psb_d_gen_pde3d
+  end interface psb_gen_pde3d
+
+contains
+
+  function d_null_func_3d(x,y,z) result(val)
+
+    real(psb_dpk_), intent(in) :: x,y,z
+    real(psb_dpk_) :: val
+
+    val = dzero
+
+  end function d_null_func_3d
+  !
+  ! functions parametrizing the differential equation
+  !
+
+  !
+  ! Note: b1, b2 and b3 are the coefficients of the first
+  ! derivative of the unknown function. The default
+  ! we apply here is to have them zero, so that the resulting
+  ! matrix is symmetric/hermitian and suitable for
+  ! testing with CG and FCG.
+  ! When testing methods for non-hermitian matrices you can
+  ! change the B1/B2/B3 functions to e.g. done/sqrt((3*done))
+  !
+  function b1(x,y,z)
+    use psb_base_mod, only : psb_dpk_, done, dzero
+    implicit none
+    real(psb_dpk_) :: b1
+    real(psb_dpk_), intent(in) :: x,y,z
+    b1=dzero
+  end function b1
+  function b2(x,y,z)
+    use psb_base_mod, only : psb_dpk_, done, dzero
+    implicit none
+    real(psb_dpk_) ::  b2
+    real(psb_dpk_), intent(in) :: x,y,z
+    b2=dzero
+  end function b2
+  function b3(x,y,z)
+    use psb_base_mod, only : psb_dpk_, done, dzero
+    implicit none
+    real(psb_dpk_) ::  b3
+    real(psb_dpk_), intent(in) :: x,y,z
+    b3=dzero
+  end function b3
+  function c(x,y,z)
+    use psb_base_mod, only : psb_dpk_, done, dzero
+    implicit none
+    real(psb_dpk_) ::  c
+    real(psb_dpk_), intent(in) :: x,y,z
+    c=dzero
+  end function c
+  function a1(x,y,z)
+    use psb_base_mod, only : psb_dpk_, done, dzero
+    implicit none
+    real(psb_dpk_) ::  a1
+    real(psb_dpk_), intent(in) :: x,y,z
+    a1=done/80
+  end function a1
+  function a2(x,y,z)
+    use psb_base_mod, only : psb_dpk_, done, dzero
+    implicit none
+    real(psb_dpk_) ::  a2
+    real(psb_dpk_), intent(in) :: x,y,z
+    a2=done/80
+  end function a2
+  function a3(x,y,z)
+    use psb_base_mod, only : psb_dpk_, done, dzero
+    implicit none
+    real(psb_dpk_) ::  a3
+    real(psb_dpk_), intent(in) :: x,y,z
+    a3=done/80
+  end function a3
+  function g(x,y,z)
+    use psb_base_mod, only : psb_dpk_, done, dzero
+    implicit none
+    real(psb_dpk_) ::  g
+    real(psb_dpk_), intent(in) :: x,y,z
+    g = dzero
+    if (x == done) then
+      g = done
+    else if (x == dzero) then
+      g = exp(y**2-z**2)
+    end if
+  end function g
+
+
+  !
+  !  subroutine to allocate and fill in the coefficient matrix and
+  !  the rhs.
+  !
+  subroutine psb_d_gen_pde3d(ctxt,idim,a,bv,xv,desc_a,afmt,info,&
+       & f,amold,vmold,imold,partition,nrl,iv)
+    use psb_base_mod
+    use psb_util_mod
+    !
+    !   Discretizes the partial differential equation
+    !
+    !   a1 dd(u)  a2 dd(u)    a3 dd(u)    b1 d(u)   b2 d(u)  b3 d(u)
+    ! -   ------ -  ------ -  ------ +  -----  +  ------  +  ------ + c u = f
+    !      dxdx     dydy       dzdz        dx       dy         dz
+    !
+    ! with Dirichlet boundary conditions
+    !   u = g
+    !
+    !  on the unit cube  0<=x,y,z<=1.
+    !
+    !
+    ! Note that if b1=b2=b3=c=0., the PDE is the  Laplace equation.
+    !
+    implicit none
+    integer(psb_ipk_)     :: idim
+    type(psb_dspmat_type) :: a
+    type(psb_d_vect_type) :: xv,bv
+    type(psb_desc_type)   :: desc_a
+    type(psb_ctxt_type)   :: ctxt
+    integer(psb_ipk_)     :: info
+    character(len=*)      :: afmt
+    procedure(d_func_3d), optional :: f
+    class(psb_d_base_sparse_mat), optional :: amold
+    class(psb_d_base_vect_type), optional :: vmold
+    class(psb_i_base_vect_type), optional :: imold
+    integer(psb_ipk_), optional :: partition, nrl,iv(:)
+
+    ! Local variables.
+
+    integer(psb_ipk_), parameter :: nb=20
+    type(psb_d_csc_sparse_mat)  :: acsc
+    type(psb_d_coo_sparse_mat)  :: acoo
+    type(psb_d_csr_sparse_mat)  :: acsr
+    real(psb_dpk_)           :: zt(nb),x,y,z
+    integer(psb_ipk_) :: nnz,nr,nlr,i,j,ii,ib,k, partition_, mysz
+    integer(psb_lpk_) :: m,n,glob_row,nt
+    integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
+    ! For 3D partition
+    ! Note: integer control variables going directly into an MPI call
+    ! must be 4 bytes, i.e. psb_mpk_
+    integer(psb_mpk_) :: npdims(3), npp, minfo
+    integer(psb_ipk_) :: npx,npy,npz, iamx,iamy,iamz,mynx,myny,mynz
+    integer(psb_ipk_), allocatable :: bndx(:),bndy(:),bndz(:)
+    ! Process grid
+    integer(psb_ipk_) :: np, iam
+    integer(psb_ipk_) :: icoeff
+    integer(psb_lpk_), allocatable     :: irow(:),icol(:),myidx(:)
+    real(psb_dpk_), allocatable :: val(:)
+    ! deltah dimension of each grid cell
+    ! deltat discretization time
+    real(psb_dpk_)            :: deltah, sqdeltah, deltah2
+    real(psb_dpk_), parameter :: rhs=dzero,one=done,zero=dzero
+    real(psb_dpk_)    :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb
+    integer(psb_ipk_) :: err_act
+    procedure(d_func_3d), pointer :: f_
+    character(len=20)  :: name, ch_err,tmpfmt
+
+    info = psb_success_
+    name = 'create_matrix'
+    call psb_erractionsave(err_act)
+
+    call psb_info(ctxt, iam, np)
+
+
+    if (present(f)) then
+      f_ => f
+    else
+      f_ => d_null_func_3d
+    end if
+
+    deltah   = done/(idim+1)
+    sqdeltah = deltah*deltah
+    deltah2  = (2*done)* deltah
+
+    if (present(partition)) then
+      if ((1<= partition).and.(partition <= 3)) then
+        partition_ = partition
+      else
+        write(*,*) 'Invalid partition choice ',partition,' defaulting to 3'
+        partition_ = 3
+      end if
+    else
+      partition_ = 3
+    end if
+
+    ! initialize array descriptor and sparse matrix storage. provide an
+    ! estimate of the number of non zeroes
+
+    m   = (1_psb_lpk_*idim)*idim*idim
+    n   = m
+    nnz = ((n*7)/(np))
+    if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n
+    t0 = psb_wtime()
+    select case(partition_)
+    case(1)
+      ! A BLOCK partition
+      if (present(nrl)) then
+        nr = nrl
+      else
+        !
+        ! Using a simple BLOCK distribution.
+        !
+        nt = (m+np-1)/np
+        nr = max(0,min(nt,m-(iam*nt)))
+      end if
+
+      nt = nr
+      call psb_sum(ctxt,nt)
+      if (nt /= m) then
+        write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m
+        info = -1
+        call psb_barrier(ctxt)
+        call psb_abort(ctxt)
+        return
+      end if
+
+      !
+      ! First example  of use of CDALL: specify for each process a number of
+      ! contiguous rows
+      !
+      call psb_cdall(ctxt,desc_a,info,nl=nr)
+      myidx = desc_a%get_global_indices()
+      nlr = size(myidx)
+
+    case(2)
+      ! A  partition  defined by the user through IV
+
+      if (present(iv)) then
+        if (size(iv) /= m) then
+          write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m
+          info = -1
+          call psb_barrier(ctxt)
+          call psb_abort(ctxt)
+          return
+        end if
+      else
+        write(psb_err_unit,*) iam, 'Initialization error: IV not present'
+        info = -1
+        call psb_barrier(ctxt)
+        call psb_abort(ctxt)
+        return
+      end if
+
+      !
+      ! Second example  of use of CDALL: specify for each row the
+      ! process that owns it
+      !
+      call psb_cdall(ctxt,desc_a,info,vg=iv)
+      myidx = desc_a%get_global_indices()
+      nlr = size(myidx)
+
+    case(3)
+      ! A 3-dimensional partition
+
+      ! A nifty MPI function will split the process list
+      npdims = 0
+#if defined(SERIAL_MPI)
+      npdims = 1
+#else
+      call mpi_dims_create(np,3,npdims,info)
+#endif
+      npx = npdims(1)
+      npy = npdims(2)
+      npz = npdims(3)
+
+      allocate(bndx(0:npx),bndy(0:npy),bndz(0:npz))
+      ! We can reuse idx2ijk for process indices as well.
+      call idx2ijk(iamx,iamy,iamz,iam,npx,npy,npz,base=0)
+      ! Now let's split the 3D cube in hexahedra
+      call dist1Didx(bndx,idim,npx)
+      mynx = bndx(iamx+1)-bndx(iamx)
+      call dist1Didx(bndy,idim,npy)
+      myny = bndy(iamy+1)-bndy(iamy)
+      call dist1Didx(bndz,idim,npz)
+      mynz = bndz(iamz+1)-bndz(iamz)
+
+      ! How many indices do I own?
+      nlr = mynx*myny*mynz
+      allocate(myidx(nlr))
+      ! Now, let's generate the list of indices I own
+      nr = 0
+      do i=bndx(iamx),bndx(iamx+1)-1
+        do j=bndy(iamy),bndy(iamy+1)-1
+          do k=bndz(iamz),bndz(iamz+1)-1
+            nr = nr + 1
+            call ijk2idx(myidx(nr),i,j,k,idim,idim,idim)
+          end do
+        end do
+      end do
+      if (nr /= nlr) then
+        write(psb_err_unit,*) iam,iamx,iamy,iamz, 'Initialization error: NR vs NLR ',&
+             & nr,nlr,mynx,myny,mynz
+        info = -1
+        call psb_barrier(ctxt)
+        call psb_abort(ctxt)
+      end if
+
+      !
+      ! Third example  of use of CDALL: specify for each process
+      ! the set of global indices it owns.
+      !
+      call psb_cdall(ctxt,desc_a,info,vl=myidx)
+
+
+      !
+      ! Specify process topology
+      !
+      block
+        !
+        ! Use adjcncy methods 
+        ! 
+        integer(psb_mpk_), allocatable :: neighbours(:)
+        integer(psb_mpk_) :: cnt
+        logical, parameter :: debug_adj=.true.
+        if (debug_adj.and.(np > 1)) then 
+          cnt = 0
+          allocate(neighbours(np))
+          if (iamx < npx-1) then
+            cnt = cnt + 1 
+            call ijk2idx(neighbours(cnt),iamx+1,iamy,iamz,npx,npy,npz,base=0)
+          end if
+          if (iamy < npy-1) then
+            cnt = cnt + 1 
+            call ijk2idx(neighbours(cnt),iamx,iamy+1,iamz,npx,npy,npz,base=0)
+          end if
+          if (iamz < npz-1) then
+            cnt = cnt + 1 
+            call ijk2idx(neighbours(cnt),iamx,iamy,iamz+1,npx,npy,npz,base=0)
+          end if
+          if (iamx >0) then
+            cnt = cnt + 1 
+            call ijk2idx(neighbours(cnt),iamx-1,iamy,iamz,npx,npy,npz,base=0)
+          end if
+          if (iamy >0) then
+            cnt = cnt + 1 
+            call ijk2idx(neighbours(cnt),iamx,iamy-1,iamz,npx,npy,npz,base=0)
+          end if
+          if (iamz >0) then
+            cnt = cnt + 1 
+            call ijk2idx(neighbours(cnt),iamx,iamy,iamz-1,npx,npy,npz,base=0)
+          end if
+          call psb_realloc(cnt, neighbours,info)
+          call desc_a%set_p_adjcncy(neighbours)
+          !write(0,*) iam,' Check on neighbours: ',desc_a%get_p_adjcncy()
+        end if
+      end block
+
+    case default
+      write(psb_err_unit,*) iam, 'Initialization error: should not get here'
+      info = -1
+      call psb_barrier(ctxt)
+      call psb_abort(ctxt)
+      return
+    end select
+
+
+    if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz, &
+         & bldmode=psb_matbld_remote_,dupl=psb_dupl_add_)
+    ! define  rhs from boundary conditions; also build initial guess
+    if (info == psb_success_) call psb_geall(xv,desc_a,info)
+    if (info == psb_success_) call psb_geall(bv,desc_a,info,&
+         & bldmode=psb_matbld_remote_,dupl=psb_dupl_add_)
+
+    call psb_barrier(ctxt)
+    talc = psb_wtime()-t0
+
+    if (info /= psb_success_) then
+      info=psb_err_from_subroutine_
+      ch_err='allocation rout.'
+      call psb_errpush(info,name,a_err=ch_err)
+      goto 9999
+    end if
+
+    ! we build an auxiliary matrix consisting of one row at a
+    ! time; just a small matrix. might be extended to generate
+    ! a bunch of rows per call.
+    !
+    allocate(val(20*nb),irow(20*nb),&
+         &icol(20*nb),stat=info)
+    if (info /= psb_success_ ) then
+      info=psb_err_alloc_dealloc_
+      call psb_errpush(info,name)
+      goto 9999
+    endif
+
+
+    ! loop over rows belonging to current process in a block
+    ! distribution.
+
+    call psb_barrier(ctxt)
+    t1 = psb_wtime()
+    do ii=1, nlr, nb
+      ib = min(nb,nlr-ii+1)
+      !ib = min(nb,mysz-ii+1)
+      icoeff = 1
+      do k=1,ib
+        i=ii+k-1
+        ! local matrix pointer
+        glob_row=myidx(i)
+        ! compute gridpoint coordinates
+        call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
+        ! x, y, z coordinates
+        x = (ix-1)*deltah
+        y = (iy-1)*deltah
+        z = (iz-1)*deltah
+        zt(k) = f_(x,y,z)
+        ! internal point: build discretization
+        !
+        !  term depending on   (x-1,y,z)
+        !
+        val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2
+        if (ix == 1) then
+          zt(k) = g(dzero,y,z)*(-val(icoeff)) + zt(k)
+        else
+          call ijk2idx(icol(icoeff),ix-1,iy,iz,idim,idim,idim)
+          irow(icoeff) = glob_row
+          icoeff       = icoeff+1
+        endif
+        !  term depending on     (x,y-1,z)
+        val(icoeff)  = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2
+        if (iy == 1) then
+          zt(k) = g(x,dzero,z)*(-val(icoeff))   + zt(k)
+        else
+          call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim)
+          irow(icoeff) = glob_row
+          icoeff       = icoeff+1
+        endif
+        !  term depending on     (x,y,z-1)
+        val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2
+        if (iz == 1) then
+          zt(k) = g(x,y,dzero)*(-val(icoeff))   + zt(k)
+        else
+          call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim)
+          irow(icoeff) = glob_row
+          icoeff       = icoeff+1
+        endif
+
+        !  term depending on     (x,y,z)
+        val(icoeff)=(2*done)*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah &
+             & + c(x,y,z)
+        call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim)
+        irow(icoeff) = glob_row
+        icoeff       = icoeff+1
+        !  term depending on     (x,y,z+1)
+        val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2
+        if (iz == idim) then
+          zt(k) = g(x,y,done)*(-val(icoeff))   + zt(k)
+        else
+          call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim)
+          irow(icoeff) = glob_row
+          icoeff       = icoeff+1
+        endif
+        !  term depending on     (x,y+1,z)
+        val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2
+        if (iy == idim) then
+          zt(k) = g(x,done,z)*(-val(icoeff))   + zt(k)
+        else
+          call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim)
+          irow(icoeff) = glob_row
+          icoeff       = icoeff+1
+        endif
+        !  term depending on     (x+1,y,z)
+        val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2
+        if (ix==idim) then
+          zt(k) = g(done,y,z)*(-val(icoeff))   + zt(k)
+        else
+          call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim)
+          irow(icoeff) = glob_row
+          icoeff       = icoeff+1
+        endif
+
+      end do
+      call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info)
+      if(info /= psb_success_) exit
+      call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info)
+      if(info /= psb_success_) exit
+      zt(:)=dzero
+      call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
+      if(info /= psb_success_) exit
+    end do
+
+    tgen = psb_wtime()-t1
+    if(info /= psb_success_) then
+      info=psb_err_from_subroutine_
+      ch_err='insert rout.'
+      call psb_errpush(info,name,a_err=ch_err)
+      goto 9999
+    end if
+
+
+    call psb_barrier(ctxt)
+    t1 = psb_wtime()
+    call psb_cdasb(desc_a,info,mold=imold)
+    tcdasb = psb_wtime()-t1
+
+    if (.false.) then 
+      !
+      !  Add extra rows to test remote build. 
+      !
+      block
+        integer(psb_ipk_) :: ks, i
+        ks = desc_a%get_local_cols()-desc_a%get_local_rows()
+        if (ks > 0) ks = max(1,ks / 10)
+        mysz = nlr+ks
+        call psb_realloc(mysz,myidx,info)
+        do i=nlr+1, mysz
+          myidx(i) = i
+        end do
+        call desc_a%l2gv1(myidx(nlr+1:mysz),info)
+        !write(0,*) iam,' Check on extra nodes ',nlr,mysz,':',myidx(nlr+1:mysz)
+        do ii= nlr+1, mysz, nb
+          ib = min(nb,mysz-ii+1)
+          icoeff = 1
+          do k=1,ib
+            i=ii+k-1
+            ! local matrix pointer
+            glob_row=myidx(i)
+            ! compute gridpoint coordinates
+            call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
+            ! x, y, z coordinates
+            x = (ix-1)*deltah
+            y = (iy-1)*deltah
+            z = (iz-1)*deltah
+            zt(k) = f_(x,y,z)
+            ! internal point: build discretization
+            !
+            !  term depending on   (x-1,y,z)
+            !
+            val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2
+            if (ix == 1) then
+              zt(k) = g(dzero,y,z)*(-val(icoeff)) + zt(k)
+            else
+              call ijk2idx(icol(icoeff),ix-1,iy,iz,idim,idim,idim)
+              irow(icoeff) = glob_row
+              icoeff       = icoeff+1
+            endif
+            !  term depending on     (x,y-1,z)
+            val(icoeff)  = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2
+            if (iy == 1) then
+              zt(k) = g(x,dzero,z)*(-val(icoeff))   + zt(k)
+            else
+              call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim)
+              irow(icoeff) = glob_row
+              icoeff       = icoeff+1
+            endif
+            !  term depending on     (x,y,z-1)
+            val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2
+            if (iz == 1) then
+              zt(k) = g(x,y,dzero)*(-val(icoeff))   + zt(k)
+            else
+              call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim)
+              irow(icoeff) = glob_row
+              icoeff       = icoeff+1
+            endif
+
+            !  term depending on     (x,y,z)
+            val(icoeff)=(2*done)*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah &
+                 & + c(x,y,z)
+            call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim)
+            irow(icoeff) = glob_row
+            icoeff       = icoeff+1
+            !  term depending on     (x,y,z+1)
+            val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2
+            if (iz == idim) then
+              zt(k) = g(x,y,done)*(-val(icoeff))   + zt(k)
+            else
+              call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim)
+              irow(icoeff) = glob_row
+              icoeff       = icoeff+1
+            endif
+            !  term depending on     (x,y+1,z)
+            val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2
+            if (iy == idim) then
+              zt(k) = g(x,done,z)*(-val(icoeff))   + zt(k)
+            else
+              call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim)
+              irow(icoeff) = glob_row
+              icoeff       = icoeff+1
+            endif
+            !  term depending on     (x+1,y,z)
+            val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2
+            if (ix==idim) then
+              zt(k) = g(done,y,z)*(-val(icoeff))   + zt(k)
+            else
+              call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim)
+              irow(icoeff) = glob_row
+              icoeff       = icoeff+1
+            endif
+
+          end do
+          call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info)
+          if(info /= psb_success_) exit
+          call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info)
+          if(info /= psb_success_) exit
+          zt(:)=dzero
+          call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
+          if(info /= psb_success_) exit
+        end do
+
+      end block
+    end if
+
+    call psb_barrier(ctxt)
+    t1 = psb_wtime()
+    if (info == psb_success_) then
+      if (present(amold)) then
+        call psb_spasb(a,desc_a,info,mold=amold)
+      else
+        call psb_spasb(a,desc_a,info,afmt=afmt)
+      end if
+    end if
+    call psb_barrier(ctxt)
+    if(info /= psb_success_) then
+      info=psb_err_from_subroutine_
+      ch_err='asb rout.'
+      call psb_errpush(info,name,a_err=ch_err)
+      goto 9999
+    end if
+    if (info == psb_success_) call psb_geasb(xv,desc_a,info,mold=vmold)
+    if (info == psb_success_) call psb_geasb(bv,desc_a,info,mold=vmold)
+    if(info /= psb_success_) then
+      info=psb_err_from_subroutine_
+      ch_err='asb rout.'
+      call psb_errpush(info,name,a_err=ch_err)
+      goto 9999
+    end if
+    tasb = psb_wtime()-t1
+    call psb_barrier(ctxt)
+    ttot = psb_wtime() - t0
+
+    call psb_amx(ctxt,talc)
+    call psb_amx(ctxt,tgen)
+    call psb_amx(ctxt,tasb)
+    call psb_amx(ctxt,ttot)
+    if(iam == psb_root_) then
+      tmpfmt = a%get_fmt()
+      write(psb_out_unit,'("The matrix has been generated and assembled in ",a3," format.")')&
+           &   tmpfmt
+      write(psb_out_unit,'("-allocation  time : ",es12.5)') talc
+      write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen
+      write(psb_out_unit,'("-desc asbly  time : ",es12.5)') tcdasb
+      write(psb_out_unit,'("- mat asbly  time : ",es12.5)') tasb
+      write(psb_out_unit,'("-total       time : ",es12.5)') ttot
+
+    end if
+    deallocate(val,irow,icol)
+    call psb_erractionrestore(err_act)
+    return
+
+9999 call psb_error_handler(ctxt,err_act)
+
+    return
+  end subroutine psb_d_gen_pde3d
+  function outside(i,j,k,bndx,bndy,bndz,iamx,iamy,iamz) result(res)
+    logical :: res
+    integer(psb_ipk_), intent(in) :: i,j,k,iamx,iamy,iamz
+    integer(psb_ipk_), intent(in) :: bndx(0:),bndy(0:),bndz(0:)
+
+    res = (i<bndx(iamx)).or.(i>=bndx(iamx+1)) &
+         & .or.(j<bndy(iamy)).or.(j>=bndy(iamy+1)) &
+         & .or.(k<bndz(iamz)).or.(k>=bndz(iamz+1))
+  end function outside
+end module psb_d_pde3d_mod
+
+program psb_tzcsrli
+  use psb_base_mod
+  use psb_prec_mod
+  use psb_krylov_mod
+  use psb_util_mod
+  use psb_d_pde3d_mod
+  implicit none
+
+  ! input parameters
+  character(len=20) :: kmethd, ptype
+  character(len=5)  :: afmt
+  integer(psb_ipk_) :: idim
+  integer(psb_epk_) :: system_size
+
+  ! miscellaneous
+  real(psb_dpk_), parameter :: one = done
+  real(psb_dpk_) :: t1, t2, tprec
+
+  ! sparse matrix and preconditioner
+  type(psb_dspmat_type) :: a
+  type(psb_zspmat_type) :: za
+  type(psb_z_csrli_sparse_mat) :: azcsrli
+  type(psb_dprec_type)  :: prec
+  ! descriptor
+  type(psb_desc_type)   :: desc_a
+  ! dense vectors
+  type(psb_d_vect_type) :: xxv,bv
+  ! parallel environment
+  type(psb_ctxt_type) :: ctxt
+  integer(psb_ipk_)   :: iam, np
+ 
+  
+  ! solver parameters
+  integer(psb_ipk_) :: iter, itmax,itrace, istopc, irst, ipart
+  integer(psb_epk_) :: amatsize, precsize, descsize, d2size
+  real(psb_dpk_)   :: err, eps
+
+  ! Parameters for solvers in Block-Jacobi preconditioner
+  type ainvparms
+    character(len=12) :: alg, orth_alg, ilu_alg, ilut_scale
+    integer(psb_ipk_) :: fill, inv_fill
+    real(psb_dpk_)    :: thresh, inv_thresh
+  end type ainvparms
+  type(ainvparms)     :: parms
+
+  ! other variables
+  integer(psb_ipk_) :: info, i
+  character(len=20) :: name,ch_err
+  character(len=40) :: fname
+
+  info=psb_success_
+
+
+  call psb_init(ctxt)
+  call psb_info(ctxt,iam,np)
+
+  if (iam < 0) then
+    ! This should not happen, but just in case
+    call psb_exit(ctxt)
+    stop
+  endif
+  if(psb_errstatus_fatal()) goto 9999
+  name='pde3d90'
+  call psb_set_errverbosity(itwo)
+  !
+  ! Hello world
+  !
+  if (iam == psb_root_) then
+    write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_
+    write(*,*) 'This is the ',trim(name),' sample program'
+  end if
+  !
+  !  get parameters
+  !
+  call get_parms(ctxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst,ipart,parms)
+  !
+  !  allocate and fill in the coefficient matrix, rhs and initial guess
+  !
+  call psb_barrier(ctxt)
+  t1 = psb_wtime()
+  call psb_gen_pde3d(ctxt,idim,a,bv,xxv,desc_a,afmt,info,partition=ipart)
+  call psb_barrier(ctxt)
+  t2 = psb_wtime() - t1
+  if(info /= psb_success_) then
+    info=psb_err_from_subroutine_
+    ch_err='psb_gen_pde3d'
+    call psb_errpush(info,name,a_err=ch_err)
+    goto 9999
+  end if
+  if (iam == psb_root_) write(psb_out_unit,'("Overall matrix creation time : ",es12.5)')t2
+  if (iam == psb_root_) write(psb_out_unit,'(" ")')
+  !
+  !  prepare the preconditioner.
+  !
+  if(iam == psb_root_) write(psb_out_unit,'("Setting preconditioner to : ",a)')ptype
+  call prec%init(ctxt,ptype,info)
+  !
+  ! Set the options for the BJAC preconditioner
+  !
+  if (psb_toupper(ptype) == "BJAC") then
+      call prec%set('sub_solve',       parms%alg,   info)
+      select case (psb_toupper(parms%alg))
+      case ("ILU")
+        call prec%set('sub_fillin',      parms%fill,       info)
+        call prec%set('ilu_alg',         parms%ilu_alg,    info)
+      case ("ILUT")
+        call prec%set('sub_fillin',      parms%fill,       info)
+        call prec%set('sub_iluthrs',     parms%thresh,     info)
+        call prec%set('ilut_scale',      parms%ilut_scale, info)
+      case ("AINV")
+        call prec%set('inv_thresh',      parms%inv_thresh, info)
+        call prec%set('inv_fillin',      parms%inv_fill,   info)
+        call prec%set('ilut_scale',      parms%ilut_scale, info)
+        call prec%set('ainv_alg',        parms%orth_alg,   info)
+      case ("INVK")
+        call prec%set('sub_fillin',      parms%fill,       info)
+        call prec%set('inv_fillin',      parms%inv_fill,   info)
+        call prec%set('ilut_scale',      parms%ilut_scale, info)
+      case ("INVT")
+        call prec%set('sub_fillin',      parms%fill,       info)
+        call prec%set('inv_fillin',      parms%inv_fill,   info)
+        call prec%set('sub_iluthrs',     parms%thresh,     info)
+        call prec%set('inv_thresh',      parms%inv_thresh, info)
+        call prec%set('ilut_scale',      parms%ilut_scale, info)
+      case default
+        ! Do nothing, use default setting in the init routine
+      end select
+  else
+    ! nothing to set for NONE or DIAG preconditioner
+  end if
+
+  call psb_barrier(ctxt)
+  t1 = psb_wtime()
+  call prec%build(a,desc_a,info)
+  if(info /= psb_success_) then
+    info=psb_err_from_subroutine_
+    ch_err='psb_precbld'
+    call psb_errpush(info,name,a_err=ch_err)
+    goto 9999
+  end if
+
+  tprec = psb_wtime()-t1
+
+  call psb_amx(ctxt,tprec)
+
+  if (iam == psb_root_) write(psb_out_unit,'("Preconditioner time : ",es12.5)')tprec
+  if (iam == psb_root_) write(psb_out_unit,'(" ")')
+  call prec%descr(info)
+  !
+  ! iterative method parameters
+  !
+  if(iam == psb_root_) write(psb_out_unit,'("Calling iterative method ",a)')kmethd
+  call psb_barrier(ctxt)
+  t1 = psb_wtime()
+  eps   = 1.d-6
+  call psb_krylov(kmethd,a,prec,bv,xxv,eps,desc_a,info,&
+       & itmax=itmax,iter=iter,err=err,itrace=itrace,istop=istopc,irst=irst)
+
+  if(info /= psb_success_) then
+    info=psb_err_from_subroutine_
+    ch_err='solver routine'
+    call psb_errpush(info,name,a_err=ch_err)
+    goto 9999
+  end if
+
+  call psb_barrier(ctxt)
+  t2 = psb_wtime() - t1
+  call psb_amx(ctxt,t2)
+  amatsize = a%sizeof()
+  descsize = desc_a%sizeof()
+  precsize = prec%sizeof()
+  system_size = desc_a%get_global_rows()
+  call psb_sum(ctxt,amatsize)
+  call psb_sum(ctxt,descsize)
+  call psb_sum(ctxt,precsize)
+
+  if (iam == psb_root_) then
+    write(psb_out_unit,'(" ")')
+    write(psb_out_unit,'("Number of processes           : ",i12)')np
+    write(psb_out_unit,'("Linear system size            : ",i12)') system_size
+    write(psb_out_unit,'("Time to solve system          : ",es12.5)')t2
+    write(psb_out_unit,'("Time per iteration            : ",es12.5)')t2/iter
+    write(psb_out_unit,'("Number of iterations          : ",i12)')iter
+    write(psb_out_unit,'("Convergence indicator on exit : ",es12.5)')err
+    write(psb_out_unit,'("Info  on exit                 : ",i12)')info
+    write(psb_out_unit,'("Total memory occupation for      A: ",i12)')amatsize
+    write(psb_out_unit,'("Total memory occupation for   PREC: ",i12)')precsize
+    write(psb_out_unit,'("Total memory occupation for DESC_A: ",i12)')descsize
+    write(psb_out_unit,'("Storage format for               A: ",a)') a%get_fmt()
+    write(psb_out_unit,'("Storage format for          DESC_A: ",a)') desc_a%get_fmt()
+  end if
+
+
+  !
+  !  cleanup storage and exit
+  !
+  call psb_gefree(bv,desc_a,info)
+  call psb_gefree(xxv,desc_a,info)
+  call psb_spfree(a,desc_a,info)
+  call prec%free(info)
+  call psb_cdfree(desc_a,info)
+  if(info /= psb_success_) then
+    info=psb_err_from_subroutine_
+    ch_err='free routine'
+    call psb_errpush(info,name,a_err=ch_err)
+    goto 9999
+  end if
+
+  call psb_exit(ctxt)
+  stop
+
+9999 call psb_error(ctxt)
+
+  stop
+
+contains
+  !
+  ! get iteration parameters from standard input
+  !
+  subroutine  get_parms(ctxt,kmethd,ptype,afmt,idim,istopc,&
+       & itmax,itrace,irst,ipart,parms)
+    type(psb_ctxt_type) :: ctxt
+    character(len=*) :: kmethd, ptype, afmt
+    integer(psb_ipk_) :: idim, istopc,itmax,itrace,irst,ipart
+    integer(psb_ipk_) :: np, iam
+    integer(psb_ipk_) :: ip, inp_unit
+    character(len=1024)   :: filename
+    type(ainvparms)   :: parms
+
+    call psb_info(ctxt, iam, np)
+
+    if (iam == 0) then
+      if (command_argument_count()>0) then
+        call get_command_argument(1,filename)
+        inp_unit = 30
+        open(inp_unit,file=filename,action='read',iostat=info)
+        if (info /= 0) then
+          write(psb_err_unit,*) 'Could not open file ',filename,' for input'
+          call psb_abort(ctxt)
+          stop
+        else
+          write(psb_err_unit,*) 'Opened file ',trim(filename),' for input'
+        end if
+      else
+        inp_unit=psb_inp_unit
+      end if
+      read(inp_unit,*) ip
+      if (ip >= 3) then
+        read(inp_unit,*) kmethd
+        read(inp_unit,*) ptype
+        read(inp_unit,*) afmt
+
+        read(inp_unit,*) idim
+        if (ip >= 4) then
+          read(inp_unit,*) ipart
+        else
+          ipart = 3
+        endif
+        if (ip >= 5) then
+          read(inp_unit,*) istopc
+        else
+          istopc=1
+        endif
+        if (ip >= 6) then
+          read(inp_unit,*) itmax
+        else
+          itmax=500
+        endif
+        if (ip >= 7) then
+          read(inp_unit,*) itrace
+        else
+          itrace=-1
+        endif
+        if (ip >= 8) then
+          read(inp_unit,*) irst
+        else
+          irst=1
+        endif
+        if (ip >= 9) then
+          read(inp_unit,*) parms%alg
+          read(inp_unit,*) parms%ilu_alg
+          read(inp_unit,*) parms%ilut_scale
+          read(inp_unit,*) parms%fill
+          read(inp_unit,*) parms%inv_fill
+          read(inp_unit,*) parms%thresh
+          read(inp_unit,*) parms%inv_thresh
+          read(inp_unit,*) parms%orth_alg
+        else
+          parms%alg =  'ILU' ! Block Solver ILU,ILUT,INVK,AINVT,AORTH
+          parms%ilu_alg = 'NONE' ! If ILU : MILU or NONE othewise ignored
+          parms%ilut_scale = 'NONE' ! If ILUT: NONE, MAXVAL, DIAG, ARWSUM, ACLSUM, ARCSUM
+          parms%fill = 0     ! Level of fill for forward factorization
+          parms%inv_fill = 1 ! Level of fill for inverse factorization (only INVK)
+          parms%thresh = 1E-1_psb_dpk_ ! Threshold for forward factorization
+          parms%inv_thresh = 1E-1_psb_dpk_ ! Threshold for inverse factorization
+          parms%orth_alg = 'LLK'  ! What orthogonalization algorithm?
+        endif
+
+        write(psb_out_unit,'("Solving matrix       : ell1")')
+        write(psb_out_unit,&
+             & '("Grid dimensions      : ",i4," x ",i4," x ",i4)') &
+             & idim,idim,idim
+        write(psb_out_unit,'("Number of processors : ",i0)')np
+        select case(ipart)
+        case(1)
+          write(psb_out_unit,'("Data distribution    : BLOCK")')
+        case(3)
+          write(psb_out_unit,'("Data distribution    : 3D")')
+        case default
+          ipart = 3
+          write(psb_out_unit,'("Unknown data distrbution, defaulting to 3D")')
+        end select
+        write(psb_out_unit,'("Preconditioner       : ",a)') ptype
+        if( psb_toupper(ptype) == "BJAC" ) then
+          write(psb_out_unit,'("Block subsolver      : ",a)') parms%alg
+          select case (psb_toupper(parms%alg))
+          case ('ILU')
+            write(psb_out_unit,'("Fill in              : ",i0)') parms%fill
+            write(psb_out_unit,'("MILU                 : ",a)') parms%ilu_alg
+          case ('ILUT')
+            write(psb_out_unit,'("Fill in              : ",i0)') parms%fill
+            write(psb_out_unit,'("Threshold            : ",es12.5)') parms%thresh
+            write(psb_out_unit,'("Scaling              : ",a)') parms%ilut_scale
+          case ('INVK')
+            write(psb_out_unit,'("Fill in              : ",i0)') parms%fill
+            write(psb_out_unit,'("Invese Fill in       : ",i0)') parms%inv_fill
+            write(psb_out_unit,'("Scaling              : ",a)') parms%ilut_scale
+          case ('INVT')
+            write(psb_out_unit,'("Fill in              : ",i0)') parms%fill
+            write(psb_out_unit,'("Threshold            : ",es12.5)') parms%thresh
+            write(psb_out_unit,'("Invese Fill in       : ",i0)') parms%inv_fill
+            write(psb_out_unit,'("Inverse Threshold    : ",es12.5)') parms%inv_thresh
+            write(psb_out_unit,'("Scaling              : ",a)') parms%ilut_scale
+          case ('AINV','AORTH')
+            write(psb_out_unit,'("Inverse Threshold    : ",es12.5)') parms%inv_thresh
+            write(psb_out_unit,'("Invese Fill in       : ",i0)') parms%inv_fill
+            write(psb_out_unit,'("Orthogonalization    : ",a)') parms%orth_alg
+            write(psb_out_unit,'("Scaling              : ",a)') parms%ilut_scale
+          case default
+              write(psb_out_unit,'("Unknown diagonal solver")')
+          end select
+        end if
+        write(psb_out_unit,'("Iterative method     : ",a)') kmethd
+        write(psb_out_unit,'(" ")')
+      else
+        ! wrong number of parameter, print an error message and exit
+        call pr_usage(izero)
+        call psb_abort(ctxt)
+        stop 1
+      endif
+      if (inp_unit /= psb_inp_unit) then
+        close(inp_unit)
+      end if
+
+    end if
+    ! broadcast parameters to all processors
+    call psb_bcast(ctxt,kmethd)
+    call psb_bcast(ctxt,afmt)
+    call psb_bcast(ctxt,ptype)
+    call psb_bcast(ctxt,idim)
+    call psb_bcast(ctxt,ipart)
+    call psb_bcast(ctxt,istopc)
+    call psb_bcast(ctxt,itmax)
+    call psb_bcast(ctxt,itrace)
+    call psb_bcast(ctxt,irst)
+    call psb_bcast(ctxt,parms%alg)
+    call psb_bcast(ctxt,parms%fill)
+    call psb_bcast(ctxt,parms%inv_fill)
+    call psb_bcast(ctxt,parms%thresh)
+    call psb_bcast(ctxt,parms%inv_thresh)
+    call psb_bcast(ctxt,parms%orth_alg)
+    call psb_bcast(ctxt,parms%ilut_scale)
+
+    return
+
+  end subroutine get_parms
+  !
+  !  print an error message
+  !
+  subroutine pr_usage(iout)
+    integer(psb_ipk_) :: iout
+    write(iout,*)'incorrect parameter(s) found'
+    write(iout,*)' usage:  pde3d90 methd prec dim &
+         &[istop itmax itrace]'
+    write(iout,*)' where:'
+    write(iout,*)'     methd:    cgstab cgs rgmres bicgstabl'
+    write(iout,*)'     prec :    bjac diag none'
+    write(iout,*)'     dim       number of points along each axis'
+    write(iout,*)'               the size of the resulting linear '
+    write(iout,*)'               system is dim**3'
+    write(iout,*)'     ipart     data partition  1  3      '
+    write(iout,*)'     istop     stopping criterion  1, 2  '
+    write(iout,*)'     itmax     maximum number of iterations [500] '
+    write(iout,*)'     itrace    <=0  (no tracing, default) or '
+    write(iout,*)'               >= 1 do tracing every itrace'
+    write(iout,*)'               iterations '
+  end subroutine pr_usage
+
+end program psb_tzcsrli