diff --git a/mlprec/Makefile b/mlprec/Makefile
index b131bd4c..a29dbbdc 100644
--- a/mlprec/Makefile
+++ b/mlprec/Makefile
@@ -28,6 +28,7 @@ INNEROBJS= mld_dcoarse_bld.o  mld_dmlprec_bld.o	 mld_dslu_bld.o mld_dumf_bld.o \
 	mld_cilu0_fact.o   mld_ciluk_fact.o  mld_cilut_fact.o  mld_caggrmap_bld.o \
 	mld_cmlprec_aply.o mld_cslud_bld.o   mld_caggrmat_asb.o   \
 	mld_zcoarse_bld.o  mld_zmlprec_bld.o mld_zslu_bld.o    mld_zumf_bld.o \
+	mld_zilu0_fact.o   mld_ziluk_fact.o  mld_zilut_fact.o  mld_zaggrmap_bld.o \
 	$(MPFOBJS) 
 #
 
diff --git a/mlprec/mld_d_ilu_solver.f03 b/mlprec/mld_d_ilu_solver.f03
index dbb9449f..ea0bea18 100644
--- a/mlprec/mld_d_ilu_solver.f03
+++ b/mlprec/mld_d_ilu_solver.f03
@@ -48,7 +48,7 @@ module mld_d_ilu_solver
   use mld_d_prec_type
 
   type, extends(mld_d_base_solver_type) :: mld_d_ilu_solver_type
-    type(psb_dspmat_type)      :: l, u
+    type(psb_dspmat_type)       :: l, u
     real(psb_dpk_), allocatable :: d(:)
     integer                     :: fact_type, fill_in
     real(psb_dpk_)              :: thresh
diff --git a/mlprec/mld_zaggrmap_bld.f90 b/mlprec/mld_zaggrmap_bld.f90
index 668887a3..7974b0d6 100644
--- a/mlprec/mld_zaggrmap_bld.f90
+++ b/mlprec/mld_zaggrmap_bld.f90
@@ -121,19 +121,20 @@ subroutine mld_zaggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info)
     call mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info)    
 
   case (mld_sym_dec_aggr_)  
-    nr = psb_sp_get_nrows(a)
-    call psb_sp_clip(a,atmp,info,imax=nr,jmax=nr,&
+    nr = a%get_nrows()
+    call a%csclip(atmp,info,imax=nr,jmax=nr,&
          & rscale=.false.,cscale=.false.)
-    atmp%m=nr
-    atmp%k=nr
-    if (info == psb_success_) call psb_transp(atmp,atrans,fmt='COO')
+    call atmp%set_nrows(nr)
+    call atmp%set_ncols(nr)
+    if (info == psb_success_) call atrans%transp(atmp)
+    if (info == psb_success_) call atrans%cscnv(info,type='COO')
     if (info == psb_success_) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.) 
-    atmp%m=nr
-    atmp%k=nr
-    if (info == psb_success_) call psb_sp_free(atrans,info)
-    if (info == psb_success_) call psb_spcnv(atmp,info,afmt='csr')
+    call atmp%set_nrows(nr)
+    call atmp%set_ncols(nr)
+    if (info == psb_success_) call atrans%free()
+    if (info == psb_success_) call atmp%cscnv(info,type='CSR')
     if (info == psb_success_) call mld_dec_map_bld(theta,atmp,desc_a,nlaggr,ilaggr,info)    
-    if (info == psb_success_) call psb_sp_free(atmp,info)
+    if (info == psb_success_) call atmp%free()
 
   case default
 
@@ -198,7 +199,7 @@ contains
     nrow  = psb_cd_get_local_rows(desc_a)
     ncol  = psb_cd_get_local_cols(desc_a)
 
-    nr = a%m
+    nr = a%get_nrows()
     allocate(ilaggr(nr),neigh(nr),stat=info)
     if(info /= psb_success_) then
       info=psb_err_alloc_request_
@@ -214,7 +215,7 @@ contains
            & a_err='complex(psb_dpk_)')
       goto 9999
     end if
-    call psb_sp_getdiag(a,diag,info)
+    call a%get_diag(diag,info)
     if(info /= psb_success_) then
       info=psb_err_from_subroutine_
       call psb_errpush(info,name,a_err='psb_sp_getdiag')
@@ -247,7 +248,7 @@ contains
           naggr     = naggr + 1 
           ilaggr(i) = naggr
 
-          call psb_sp_getrow(i,a,nz,irow,icol,val,info)
+          call a%csget(i,i,nz,irow,icol,val,info)
           if (info /= psb_success_) then 
             info=psb_err_from_subroutine_
             call psb_errpush(info,name,a_err='psb_sp_getrow')
@@ -268,7 +269,7 @@ contains
           !
           ! 2. Untouched neighbours of these nodes are marked <0.
           !
-          call psb_neigh(a,i,neigh,n_ne,info,lev=2)
+          call a%get_neigh(i,neigh,n_ne,info,lev=2)
           if (info /= psb_success_) then 
             info=psb_err_from_subroutine_
             call psb_errpush(info,name,a_err='psb_neigh')
@@ -288,8 +289,7 @@ contains
     enddo
     if (debug_level >= psb_debug_outer_) then 
       write(debug_unit,*) me,' ',trim(name),&
-           & ' Check 1:',count(ilaggr == -(nr+1)),&
-           & (a%ia1(i),i=a%ia2(1),a%ia2(2)-1)
+           & ' Check 1:',count(ilaggr == -(nr+1))
     end if
 
     !
@@ -336,7 +336,7 @@ contains
         isz  = nr+1
         ia   = -1
         cpling = dzero
-        call psb_sp_getrow(i,a,nz,irow,icol,val,info)
+        call a%csget(i,i,nz,irow,icol,val,info)
         if (info /= psb_success_) then 
           info=psb_err_from_subroutine_
           call psb_errpush(info,name,a_err='psb_sp_getrow')
diff --git a/mlprec/mld_zilu0_fact.f90 b/mlprec/mld_zilu0_fact.f90
index 13a9ab7a..2c5ea6d5 100644
--- a/mlprec/mld_zilu0_fact.f90
+++ b/mlprec/mld_zilu0_fact.f90
@@ -99,10 +99,10 @@
 !               greater than 0. If the overlap is 0 or the matrix has been reordered
 !               (see mld_fact_bld), then blck is empty.
 !  
-subroutine mld_zilu0_fact(ialg,a,l,u,d,info,blck)
+subroutine mld_zilu0_fact(ialg,a,l,u,d,info,blck,upd)
 
   use psb_sparse_mod
-  use mld_inner_mod, mld_protect_name => mld_zilu0_fact
+  use mld_inner_mod!, mld_protect_name => mld_zilu0_fact
 
   implicit none
 
@@ -113,11 +113,14 @@ subroutine mld_zilu0_fact(ialg,a,l,u,d,info,blck)
   complex(psb_dpk_), intent(inout)     :: d(:)
   integer, intent(out)                :: info
   type(psb_zspmat_type),intent(in), optional, target :: blck
+  character, intent(in), optional      :: upd
 
   ! Local variables
-  integer   :: l1, l2,m,err_act
+  integer   :: l1, l2, m, err_act
   type(psb_zspmat_type), pointer  :: blck_
-  character(len=20)   :: name, ch_err
+  type(psb_z_csr_sparse_mat)       :: ll, uu
+  character                        :: upd_
+  character(len=20)    :: name, ch_err
 
   name='mld_zilu0_fact'
   info = psb_success_
@@ -130,28 +133,36 @@ subroutine mld_zilu0_fact(ialg,a,l,u,d,info,blck)
     blck_ => blck
   else
     allocate(blck_,stat=info) 
-    if (info /= psb_success_) then 
-      call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
-      goto 9999      
-    end if
-
-    call psb_nullify_sp(blck_)          ! Probably pointless.
-    call psb_sp_all(0,0,blck_,1,info)
-    if(info.ne.0) then
+    if (info == psb_success_) call blck_%csall(0,0,info,1) 
+    if (info /= psb_success_) then
       info=psb_err_from_subroutine_
-      ch_err='psb_sp_all'
+      ch_err='csall'
       call psb_errpush(info,name,a_err=ch_err)
       goto 9999
     end if
-
-    blck_%m=0
   endif
+  if (present(upd)) then 
+    upd_ = psb_toupper(upd) 
+  else
+    upd_ = 'F'
+  end if
+  
+  m = a%get_nrows() + blck_%get_nrows()
+  if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.&
+       & (m > size(d))    ) then 
+    write(0,*) 'Wrong allocation status for L,D,U? ',&
+         & l%get_nrows(),size(d),u%get_nrows()
+    info = -1
+    return
+  end if
 
+  call l%mv_to(ll)
+  call u%mv_to(uu)
   !
   ! Compute the ILU(0) or the MILU(0) factorization, depending on ialg
   !
-  call mld_zilu0_factint(ialg,m,a%m,a,blck_%m,blck_,&
-       & d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info)
+  call mld_zilu0_factint(ialg,a,blck_,&
+       & d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,upd_,info)
   if(info.ne.0) then
      info=psb_err_from_subroutine_
      ch_err='mld_zilu0_factint'
@@ -162,24 +173,22 @@ subroutine mld_zilu0_fact(ialg,a,l,u,d,info,blck)
   !
   ! Store information on the L and U sparse matrices
   !
-  l%infoa(1) = l1
-  l%fida     = 'CSR'
-  l%descra   = 'TLU'
-  u%infoa(1) = l2
-  u%fida     = 'CSR'
-  u%descra   = 'TUU'
-  l%m = m
-  l%k = m
-  u%m = m
-  u%k = m
-
+  call l%mv_from(ll)
+  call l%set_triangle()
+  call l%set_unit()
+  call l%set_lower()
+  call u%mv_from(uu)
+  call u%set_triangle()
+  call u%set_unit()
+  call u%set_upper()
+  
   !
   ! Nullify pointer / deallocate memory
   !
   if (present(blck)) then 
     blck_ => null() 
   else
-    call psb_sp_free(blck_,info)
+    call blck_%free()
     if(info.ne.0) then
       info=psb_err_from_subroutine_
       ch_err='psb_sp_free'
@@ -277,24 +286,25 @@ contains
   !    info    -  integer, output.           
   !               Error code.
   !
-  subroutine mld_zilu0_factint(ialg,m,ma,a,mb,b,&
-       & d,laspk,lia1,lia2,uaspk,uia1,uia2,l1,l2,info)
+  subroutine mld_zilu0_factint(ialg,a,b,&
+       & d,lval,lja,lirp,uval,uja,uirp,l1,l2,upd,info)
 
     implicit none 
 
     ! Arguments
-    integer, intent(in)                              :: ialg
-    type(psb_zspmat_type),intent(in)                 :: a,b
-    integer,intent(inout)                            :: m,l1,l2,info
-    integer, intent(in)                              :: ma,mb
-    integer, dimension(:), intent(inout)             :: lia1,lia2,uia1,uia2
-    complex(psb_dpk_), dimension(:), intent(inout) :: laspk,uaspk,d
+    integer, intent(in)               :: ialg
+    type(psb_zspmat_type),intent(in)  :: a,b
+    integer,intent(inout)             :: l1,l2,info
+    integer, intent(inout)            :: lja(:),lirp(:),uja(:),uirp(:)
+    complex(psb_dpk_), intent(inout)  :: lval(:),uval(:),d(:)
+    character, intent(in)             :: upd
 
     ! Local variables
-    integer :: i,j,k,l,low1,low2,kk,jj,ll, ktrw,err_act
+    integer :: i,j,k,l,low1,low2,kk,jj,ll, ktrw,err_act, m
+    integer          :: ma,mb
     complex(psb_dpk_) :: dia,temp
     integer, parameter :: nrb=16
-    type(psb_zspmat_type) :: trw
+    type(psb_z_coo_sparse_mat) :: trw
     integer             :: int_err(5) 
     character(len=20)   :: name, ch_err
 
@@ -302,6 +312,8 @@ contains
     if(psb_get_errstatus().ne.0) return 
     info=psb_success_
     call psb_erractionsave(err_act)
+    ma = a%get_nrows()
+    mb = b%get_nrows()
 
     select case(ialg)
     case(mld_ilu_n_,mld_milu_n_)
@@ -312,154 +324,152 @@ contains
       goto 9999
     end select
 
-    call psb_nullify_sp(trw)
-    trw%m=0
-    trw%k=0
-
-    call psb_sp_all(trw,1,info)
-    if(info.ne.0) then
+    call trw%allocate(0,0,1)
+    if(info /= psb_success_) then
       info=psb_err_from_subroutine_
       ch_err='psb_sp_all'
       call psb_errpush(info,name,a_err=ch_err)
       goto 9999
     end if
-
-    lia2(1) = 1
-    uia2(1) = 1
-    l1      = 0
-    l2      = 0
     m = ma+mb
 
-    !
-    ! Cycle over the matrix rows
-    !
-    do i = 1, m
+    if (psb_toupper(upd) == 'F' ) then 
+      lirp(1) = 1
+      uirp(1) = 1
+      l1      = 0
+      l2      = 0
+
+      !
+      ! Cycle over the matrix rows
+      !
+      do i = 1, m
 
-      d(i) = zzero
+      d(i) = szero
 
-      if (i <= ma) then
-        !
-        ! Copy the i-th local row of the matrix, stored in a,
-        ! into laspk/d(i)/uaspk 
-        !
-        call ilu_copyin(i,ma,a,i,1,m,l1,lia1,laspk,&
-             & d(i),l2,uia1,uaspk,ktrw,trw)
-      else
-        !
-        ! Copy the i-th local row of the matrix, stored in b
-        ! (as (i-ma)-th row), into laspk/d(i)/uaspk 
-        !
-        call ilu_copyin(i-ma,mb,b,i,1,m,l1,lia1,laspk,&
-             & d(i),l2,uia1,uaspk,ktrw,trw)
-      endif
+        if (i <= ma) then
+          !
+          ! Copy the i-th local row of the matrix, stored in a,
+          ! into lval/d(i)/uval 
+          !
+          call ilu_copyin(i,ma,a,i,1,m,l1,lja,lval,&
+               & d(i),l2,uja,uval,ktrw,trw,upd)
+        else
+          !
+          ! Copy the i-th local row of the matrix, stored in b
+          ! (as (i-ma)-th row), into lval/d(i)/uval 
+          !
+          call ilu_copyin(i-ma,mb,b,i,1,m,l1,lja,lval,&
+               & d(i),l2,uja,uval,ktrw,trw,upd)
+        endif
 
-      lia2(i+1) = l1 + 1
-      uia2(i+1) = l2 + 1
+        lirp(i+1) = l1 + 1
+        uirp(i+1) = l2 + 1
 
-      dia = d(i)
-      do kk = lia2(i), lia2(i+1) - 1
-        !
-        ! Compute entry l(i,k) (lower factor L) of the incomplete
-        ! factorization
-        !
-        temp      = laspk(kk)
-        k         = lia1(kk)
-        laspk(kk) = temp*d(k)
-        !
-        ! Update the rest of row i (lower and upper factors L and U)
-        ! using l(i,k)
-        !
-        low1 = kk + 1
-        low2 = uia2(i)
-        !
-        updateloop: do  jj = uia2(k), uia2(k+1) - 1
+        dia = d(i)
+        do kk = lirp(i), lirp(i+1) - 1
           !
-          j = uia1(jj)
+          ! Compute entry l(i,k) (lower factor L) of the incomplete
+          ! factorization
           !
-          if (j < i) then
-            !
-            ! search l(i,*) (i-th row of L) for a matching index j
-            !
-            do  ll = low1, lia2(i+1) - 1
-              l = lia1(ll)
-              if (l > j) then
-                low1 = ll
-                exit
-              else if (l == j) then
-                laspk(ll) = laspk(ll) - temp*uaspk(jj)
-                low1 = ll + 1
-                cycle updateloop
-              end if
-            enddo
-
-          else if (j == i) then
+          temp      = lval(kk)
+          k         = lja(kk)
+          lval(kk) = temp*d(k)
+          !
+          ! Update the rest of row i (lower and upper factors L and U)
+          ! using l(i,k)
+          !
+          low1 = kk + 1
+          low2 = uirp(i)
+          !
+          updateloop: do  jj = uirp(k), uirp(k+1) - 1
             !
-            ! j=i: update the diagonal
+            j = uja(jj)
             !
-            dia = dia - temp*uaspk(jj)
-            cycle updateloop
+            if (j < i) then
+              !
+              ! search l(i,*) (i-th row of L) for a matching index j
+              !
+              do  ll = low1, lirp(i+1) - 1
+                l = lja(ll)
+                if (l > j) then
+                  low1 = ll
+                  exit
+                else if (l == j) then
+                  lval(ll) = lval(ll) - temp*uval(jj)
+                  low1 = ll + 1
+                  cycle updateloop
+                end if
+              enddo
+
+            else if (j == i) then
+              !
+              ! j=i: update the diagonal
+              !
+              dia = dia - temp*uval(jj)
+              cycle updateloop
+              !     
+            else if (j > i) then
+              !
+              ! search u(i,*) (i-th row of U) for a matching index j
+              !
+              do ll = low2, uirp(i+1) - 1
+                l = uja(ll)
+                if (l > j) then
+                  low2 = ll
+                  exit
+                else if (l == j) then
+                  uval(ll) = uval(ll) - temp*uval(jj)
+                  low2 = ll + 1
+                  cycle updateloop
+                end if
+              enddo
+            end if
             !     
-          else if (j > i) then
-            !
-            ! search u(i,*) (i-th row of U) for a matching index j
+            ! If we get here we missed the cycle updateloop, which means 
+            ! that this entry does not match; thus we accumulate on the
+            ! diagonal for MILU(0).
             !
-            do ll = low2, uia2(i+1) - 1
-              l = uia1(ll)
-              if (l > j) then
-                low2 = ll
-                exit
-              else if (l == j) then
-                uaspk(ll) = uaspk(ll) - temp*uaspk(jj)
-                low2 = ll + 1
-                cycle updateloop
-              end if
-            enddo
-          end if
+            if (ialg == mld_milu_n_) then 
+              dia = dia - temp*uval(jj)
+            end if
+          enddo updateloop
+        enddo
+        !     
+        ! Check the pivot size
+        !     
+        if (abs(dia) < s_epstol) then
+          !
+          ! Too small pivot: unstable factorization
           !     
-          ! If we get here we missed the cycle updateloop, which means 
-          ! that this entry does not match; thus we accumulate on the
-          ! diagonal for MILU(0).
+          info = psb_err_pivot_too_small_
+          int_err(1) = i
+          write(ch_err,'(g20.10)') abs(dia)
+          call psb_errpush(info,name,i_err=int_err,a_err=ch_err)
+          goto 9999
+        else
           !
-          if (ialg == mld_milu_n_) then 
-            dia = dia - temp*uaspk(jj)
-          end if
-        enddo updateloop
-      enddo
-      !     
-      ! Check the pivot size
-      !     
-      if (abs(dia) < d_epstol) then
-        !
-        ! Too small pivot: unstable factorization
-        !     
-        info = psb_err_pivot_too_small_
-        int_err(1) = i
-        write(ch_err,'(g20.10)') abs(dia)
-        call psb_errpush(info,name,i_err=int_err,a_err=ch_err)
-        goto 9999
-      else
+          ! Compute 1/pivot
+          !
+          dia = zone/dia
+        end if
+        d(i) = dia
         !
-        ! Compute 1/pivot
+        ! Scale row i of upper triangle
         !
-        dia = done/dia
-      end if
-      d(i) = dia
-      !
-      ! Scale row i of upper triangle
-      !
-      do  kk = uia2(i), uia2(i+1) - 1
-        uaspk(kk) = uaspk(kk)*dia
+        do  kk = uirp(i), uirp(i+1) - 1
+          uval(kk) = uval(kk)*dia
+        enddo
       enddo
-    enddo
-
-    call psb_sp_free(trw,info)
-    if(info.ne.0) then
-      info=psb_err_from_subroutine_
-      ch_err='psb_sp_free'
-      call psb_errpush(info,name,a_err=ch_err)
+    else
+      write(0,*) 'Update not implemented '
+      info = 31
+      call psb_errpush(info,name,i_err=(/13,0,0,0,0/),a_err=upd)
       goto 9999
+
     end if
 
+    call trw%free()
+
     call psb_erractionrestore(err_act)
     return
 
@@ -480,13 +490,13 @@ contains
   !  This routine copies a row of a sparse matrix A, stored in the psb_dspmat_type 
   !  data structure a, into the arrays laspk and uaspk and into the scalar variable
   !  dia, corresponding to the lower and upper triangles of A and to the diagonal
-  !  entry of the row, respectively. The entries in laspk and uaspk are stored
+  !  entry of the row, respectively. The entries in lval and uval are stored
   !  according to the CSR format; the corresponding column indices are stored in
-  !  the arrays lia1 and uia1.
+  !  the arrays lja and uja.
   !
   !  If the sparse matrix is in CSR format, a 'straight' copy is performed;
   !  otherwise psb_sp_getblk is used to extract a block of rows, which is then
-  !  copied into laspk, dia, uaspk row by row, through successive calls to
+  !  copied into lval, dia, uval row by row, through successive calls to
   !  ilu_copyin.
   !
   !  The routine is used by mld_zilu0_factint in the computation of the ILU(0)/MILU(0)
@@ -514,23 +524,23 @@ contains
   !               The output matrix will contain a clipped copy taken from
   !               a(1:m,jmin:jmax).
   !    l1      -  integer, input/output.
-  !               Pointer to the last occupied entry of laspk.
-  !    lia1    -  integer, dimension(:), input/output.
+  !               Pointer to the last occupied entry of lval.
+  !    lja    -  integer, dimension(:), input/output.
   !               The column indices of the nonzero entries of the lower triangle
-  !               copied in laspk row by row (see mld_zilu0_factint), according
+  !               copied in lval row by row (see mld_dilu0_factint), according
   !               to the CSR storage format.
-  !    laspk   -  complex(psb_dpk_), dimension(:), input/output.
+  !    lval   -  complex(psb_dpk_), dimension(:), input/output.
   !               The array where the entries of the row corresponding to the
   !               lower triangle are copied.
   !    dia     -  complex(psb_dpk_), output.
   !               The diagonal entry of the copied row.
   !    l2      -  integer, input/output.
-  !               Pointer to the last occupied entry of uaspk.
-  !    uia1    -  integer, dimension(:), input/output.
+  !               Pointer to the last occupied entry of uval.
+  !    uja    -  integer, dimension(:), input/output.
   !               The column indices of the nonzero entries of the upper triangle
-  !               copied in uaspk row by row (see mld_zilu0_factint), according
+  !               copied in uval row by row (see mld_dilu0_factint), according
   !               to the CSR storage format.
-  !    uaspk   -  complex(psb_dpk_), dimension(:), input/output.
+  !    uval   -  complex(psb_dpk_), dimension(:), input/output.
   !               The array where the entries of the row corresponding to the
   !               upper triangle are copied. 
   !    ktrw    -  integer, input/output.
@@ -544,8 +554,8 @@ contains
   !               until we empty the buffer. Thus we will make a call to psb_sp_getblk
   !               every nrb calls to copyin. If A is in CSR format it is unused.
   !
-  subroutine ilu_copyin(i,m,a,jd,jmin,jmax,l1,lia1,laspk,&
-       & dia,l2,uia1,uaspk,ktrw,trw)
+  subroutine ilu_copyin(i,m,a,jd,jmin,jmax,l1,lja,lval,&
+       & dia,l2,uja,uval,ktrw,trw,upd)
 
     use psb_sparse_mod
 
@@ -553,83 +563,95 @@ contains
 
     ! Arguments
     type(psb_zspmat_type), intent(in)    :: a
-    type(psb_zspmat_type), intent(inout) :: trw
+    type(psb_z_coo_sparse_mat), intent(inout) :: trw
     integer, intent(in)                  :: i,m,jd,jmin,jmax
     integer, intent(inout)               :: ktrw,l1,l2
-    integer, intent(inout)               :: lia1(:), uia1(:)
-    complex(psb_dpk_), intent(inout)   :: laspk(:), uaspk(:), dia
-
+    integer, intent(inout)               :: lja(:), uja(:)
+    complex(psb_dpk_), intent(inout)     :: lval(:), uval(:), dia
+    character, intent(in)                :: upd
     ! Local variables
-    integer               :: k,j,info,irb
-    integer, parameter    :: nrb=16
+    integer               :: k,j,info,irb, nz
+    integer, parameter    :: nrb=40
     character(len=20), parameter  :: name='ilu_copyin'
     character(len=20)             :: ch_err
 
     if (psb_get_errstatus() /= 0) return 
     info=psb_success_
     call psb_erractionsave(err_act)
+    if (psb_toupper(upd) == 'F') then 
 
-    if (psb_toupper(a%fida) == 'CSR') then
+      select type(aa => a%a) 
+      type is (psb_z_csr_sparse_mat) 
 
       !
       ! Take a fast shortcut if the matrix is stored in CSR format
       !
 
-      do j = a%ia2(i), a%ia2(i+1) - 1
-        k = a%ia1(j)
-        !           write(0,*)'KKKKK',k
-        if ((k < jd).and.(k >= jmin)) then
-          l1 = l1 + 1
-          laspk(l1) = a%aspk(j)
-          lia1(l1) = k
-        else if (k == jd) then
-          dia = a%aspk(j)
-        else if ((k > jd).and.(k <= jmax)) then
-          l2 = l2 + 1
-          uaspk(l2) = a%aspk(j)
-          uia1(l2) = k
-        end if
-      enddo
+        do j = aa%irp(i), aa%irp(i+1) - 1
+          k = aa%ja(j)
+          !           write(0,*)'KKKKK',k
+          if ((k < jd).and.(k >= jmin)) then
+            l1 = l1 + 1
+            lval(l1) = aa%val(j)
+            lja(l1) = k
+          else if (k == jd) then
+            dia = aa%val(j)
+          else if ((k > jd).and.(k <= jmax)) then
+            l2 = l2 + 1
+            uval(l2) = aa%val(j)
+            uja(l2) = k
+          end if
+        enddo
 
-    else
+      class default 
 
-      !
-      ! Otherwise use psb_sp_getblk, slower but able (in principle) of 
-      ! handling any format. In this case, a block of rows is extracted
-      ! instead of a single row, for performance reasons, and these
-      ! rows are copied one by one into laspk, dia, uaspk, through
-      ! successive calls to ilu_copyin.
-      !
+        !
+        ! Otherwise use psb_sp_getblk, slower but able (in principle) of 
+        ! handling any format. In this case, a block of rows is extracted
+        ! instead of a single row, for performance reasons, and these
+        ! rows are copied one by one into lval, dia, uval, through
+        ! successive calls to ilu_copyin.
+        !
 
-      if ((mod(i,nrb) == 1).or.(nrb == 1)) then 
-        irb = min(m-i+1,nrb)
-        call psb_sp_getblk(i,a,trw,info,lrw=i+irb-1)
-        if(info.ne.0) then
-          info=psb_err_from_subroutine_
-          ch_err='psb_sp_getblk'
-          call psb_errpush(info,name,a_err=ch_err)
-          goto 9999
-        end if
-        ktrw=1
-      end if
-
-      do 
-        if (ktrw > trw%infoa(psb_nnz_)) exit
-        if (trw%ia1(ktrw) > i) exit
-        k = trw%ia2(ktrw)
-        if ((k < jd).and.(k >= jmin)) then
-          l1 = l1 + 1
-          laspk(l1) = trw%aspk(ktrw)
-          lia1(l1) = k
-        else if (k == jd) then
-          dia = trw%aspk(ktrw)
-        else if ((k > jd).and.(k <= jmax)) then
-          l2 = l2 + 1
-          uaspk(l2) = trw%aspk(ktrw)
-          uia1(l2) = k
+        if ((mod(i,nrb) == 1).or.(nrb == 1)) then 
+          irb = min(m-i+1,nrb)
+          call aa%csget(i,i+irb-1,trw,info) 
+          if(info /= psb_success_) then
+            info=psb_err_from_subroutine_
+            ch_err='csget'
+            call psb_errpush(info,name,a_err=ch_err)
+            goto 9999
+          end if
+          ktrw=1
         end if
-        ktrw = ktrw + 1
-      enddo
+
+        nz = trw%get_nzeros()
+        do 
+          if (ktrw > nz) exit
+          if (trw%ia(ktrw) > i) exit
+          k = trw%ja(ktrw)
+          if ((k < jd).and.(k >= jmin)) then
+            l1 = l1 + 1
+            lval(l1) = trw%val(ktrw)
+            lja(l1) = k
+          else if (k == jd) then
+            dia = trw%val(ktrw)
+          else if ((k > jd).and.(k <= jmax)) then
+            l2 = l2 + 1
+            uval(l2) = trw%val(ktrw)
+            uja(l2) = k
+          end if
+          ktrw = ktrw + 1
+        enddo
+
+      end select
+
+    else
+
+      write(0,*) 'Update not implemented '
+      info = 31
+      call psb_errpush(info,name,i_err=(/13,0,0,0,0/),a_err=upd)
+      goto 9999
 
     end if
 
diff --git a/mlprec/mld_ziluk_fact.f90 b/mlprec/mld_ziluk_fact.f90
index 57207bf0..5f350da8 100644
--- a/mlprec/mld_ziluk_fact.f90
+++ b/mlprec/mld_ziluk_fact.f90
@@ -99,7 +99,7 @@
 subroutine mld_ziluk_fact(fill_in,ialg,a,l,u,d,info,blck)
 
   use psb_sparse_mod
-  use mld_inner_mod, mld_protect_name => mld_ziluk_fact
+  use mld_inner_mod!, mld_protect_name => mld_ziluk_fact
 
   implicit none
 
@@ -114,6 +114,7 @@ subroutine mld_ziluk_fact(fill_in,ialg,a,l,u,d,info,blck)
   integer   :: l1, l2, m, err_act
   
   type(psb_zspmat_type), pointer  :: blck_
+  type(psb_z_csr_sparse_mat)       :: ll, uu
   character(len=20)   :: name, ch_err
 
   name='mld_ziluk_fact'
@@ -127,26 +128,32 @@ subroutine mld_ziluk_fact(fill_in,ialg,a,l,u,d,info,blck)
     blck_ => blck
   else
     allocate(blck_,stat=info) 
-    if (info /= psb_success_) then 
-      call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
-      goto 9999      
-    end if
-
-    call psb_sp_all(0,0,blck_,1,info)
+    if (info == psb_success_) call blck_%csall(0,0,info,1) 
     if (info /= psb_success_) then
-       info=psb_err_from_subroutine_
-       ch_err='psb_sp_all'
-       call psb_errpush(info,name,a_err=ch_err)
-       goto 9999
+      info=psb_err_from_subroutine_
+      ch_err='csall'
+      call psb_errpush(info,name,a_err=ch_err)
+      goto 9999
     end if
-
   endif
 
+  m = a%get_nrows() + blck_%get_nrows()
+  if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.&
+       & (m > size(d))    ) then 
+    write(0,*) 'Wrong allocation status for L,D,U? ',&
+         & l%get_nrows(),size(d),u%get_nrows()
+    info = -1
+    return
+  end if
+
+  call l%mv_to(ll)
+  call u%mv_to(uu)
+
   !
   ! Compute the ILU(k) or the MILU(k) factorization, depending on ialg
   !
-  call mld_ziluk_factint(fill_in,ialg,m,a,blck_,&
-       & d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info)
+  call mld_ziluk_factint(fill_in,ialg,a,blck_,&
+       & d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,info)
   if (info /= psb_success_) then
      info=psb_err_from_subroutine_
      ch_err='mld_ziluk_factint'
@@ -157,33 +164,32 @@ subroutine mld_ziluk_fact(fill_in,ialg,a,l,u,d,info,blck)
   !
   ! Store information on the L and U sparse matrices
   !
-  l%infoa(1) = l1
-  l%fida     = 'CSR'
-  l%descra   = 'TLU'
-  u%infoa(1) = l2
-  u%fida     = 'CSR'
-  u%descra   = 'TUU'
-  l%m = m
-  l%k = m
-  u%m = m
-  u%k = m
-
+  call l%mv_from(ll)
+  call l%set_triangle()
+  call l%set_unit()
+  call l%set_lower()
+  call u%mv_from(uu)
+  call u%set_triangle()
+  call u%set_unit()
+  call u%set_upper()
+  
   !
-  ! Nullify the pointer / deallocate the memory
+  ! Nullify pointer / deallocate memory
   !
   if (present(blck)) then 
     blck_ => null() 
   else
-    call psb_sp_free(blck_,info)
-    if (info /= psb_success_) then
-       info=psb_err_from_subroutine_
-       ch_err='psb_sp_free'
-       call psb_errpush(info,name,a_err=ch_err)
-       goto 9999
+    call blck_%free()
+    deallocate(blck_,stat=info) 
+    if(info.ne.0) then
+      info=psb_err_from_subroutine_
+      ch_err='psb_sp_free'
+      call psb_errpush(info,name,a_err=ch_err)
+      goto 9999
     end if
-    deallocate(blck_) 
   endif
 
+
   call psb_erractionrestore(err_act)
   return
 
@@ -248,43 +254,43 @@ contains
   !    lia2    -  integer, dimension(:), input/output.
   !               The indices identifying the first nonzero entry of each row
   !               of the L factor in laspk, according to the CSR storage format. 
-  !    uaspk   -  complex(psb_dpk_), dimension(:), input/output.
+  !    uval   -  complex(psb_dpk_), dimension(:), input/output.
   !               The U factor in the incomplete factorization.
   !               The entries of U are stored according to the CSR format.
-  !    uia1    -  integer, dimension(:), input/output.
+  !    uja    -  integer, dimension(:), input/output.
   !               The column indices of the nonzero entries of the U factor,
   !               according to the CSR storage format.
-  !    uia2    -  integer, dimension(:), input/output.
+  !    uirp    -  integer, dimension(:), input/output.
   !               The indices identifying the first nonzero entry of each row
-  !               of the U factor in uaspk, according to the CSR storage format. 
+  !               of the U factor in uval, according to the CSR storage format. 
   !    l1      -  integer, output
   !               The number of nonzero entries in laspk.
   !    l2      -  integer, output
-  !               The number of nonzero entries in uaspk.
+  !               The number of nonzero entries in uval.
   !    info    -  integer, output.           
   !               Error code.
   !
-  subroutine mld_ziluk_factint(fill_in,ialg,m,a,b,&
-       & d,laspk,lia1,lia2,uaspk,uia1,uia2,l1,l2,info)
+  subroutine mld_ziluk_factint(fill_in,ialg,a,b,&
+       & d,lval,lja,lirp,uval,uja,uirp,l1,l2,info)
 
     use psb_sparse_mod
 
     implicit none
 
   ! Arguments 
-    integer, intent(in)                             :: fill_in, ialg
-    type(psb_zspmat_type), intent(in)               :: a,b
-    integer, intent(inout)                          :: m,l1,l2,info
-    integer, allocatable, intent(inout)             :: lia1(:),lia2(:),uia1(:),uia2(:)
-    complex(psb_dpk_), allocatable, intent(inout) :: laspk(:),uaspk(:)
+    integer, intent(in)                        :: fill_in, ialg
+    type(psb_zspmat_type),intent(in)           :: a,b
+    integer,intent(inout)                      :: l1,l2,info
+    integer, allocatable, intent(inout)        :: lja(:),lirp(:),uja(:),uirp(:)
+    complex(psb_dpk_), allocatable, intent(inout) :: lval(:),uval(:)
     complex(psb_dpk_), intent(inout)              :: d(:)
 
   ! Local variables
-    integer :: ma,mb,i, ktrw,err_act,nidx
-    integer, allocatable            :: uplevs(:), rowlevs(:),idxs(:)
-    complex(psb_dpk_), allocatable :: row(:)
+    integer :: ma,mb,i, ktrw,err_act,nidx, m
+    integer, allocatable          :: uplevs(:), rowlevs(:),idxs(:)
+    complex(psb_dpk_), allocatable   :: row(:)
     type(psb_int_heap) :: heap
-    type(psb_zspmat_type) :: trw
+    type(psb_z_coo_sparse_mat) :: trw
     character(len=20), parameter  :: name='mld_ziluk_factint'
     character(len=20)             :: ch_err
 
@@ -292,6 +298,7 @@ contains
     info=psb_success_
     call psb_erractionsave(err_act)
 
+
     select case(ialg)
     case(mld_ilu_n_,mld_milu_n_)
       ! Ok 
@@ -306,16 +313,17 @@ contains
       goto 9999
     end if
 
-    ma = a%m
-    mb = b%m
+    ma = a%get_nrows()
+    mb = b%get_nrows()
     m  = ma+mb
 
     !
     ! Allocate a temporary buffer for the iluk_copyin function 
     !
-    call psb_sp_all(0,0,trw,1,info)
-    if (info == psb_success_) call psb_ensure_size(m+1,lia2,info)
-    if (info == psb_success_) call psb_ensure_size(m+1,uia2,info)
+
+    call trw%allocate(0,0,1)
+    if (info == psb_success_) call psb_ensure_size(m+1,lirp,info)
+    if (info == psb_success_) call psb_ensure_size(m+1,uirp,info)
 
     if (info /= psb_success_) then
       info=psb_err_from_subroutine_
@@ -325,14 +333,14 @@ contains
     
     l1=0
     l2=0
-    lia2(1) = 1
-    uia2(1) = 1
+    lirp(1) = 1
+    uirp(1) = 1
 
     !
     ! Allocate memory to hold the entries of a row and the corresponding
     ! fill levels
     !
-    allocate(uplevs(size(uaspk)),rowlevs(m),row(m),stat=info)
+    allocate(uplevs(size(uval)),rowlevs(m),row(m),stat=info)
     if (info /= psb_success_) then
       info=psb_err_from_subroutine_
       call psb_errpush(info,name,a_err='Allocate')
@@ -375,12 +383,12 @@ contains
       ! do not have a lowlevs variable.
       !
       if (info == psb_success_) call iluk_fact(fill_in,i,row,rowlevs,heap,&
-           & d,uia1,uia2,uaspk,uplevs,nidx,idxs,info)
+           & d,uja,uirp,uval,uplevs,nidx,idxs,info)
       !
-      ! Copy the row into laspk/d(i)/uaspk
+      ! Copy the row into lval/d(i)/uval
       ! 
       if (info == psb_success_) call iluk_copyout(fill_in,ialg,i,m,row,rowlevs,nidx,idxs,&
-           & l1,l2,lia1,lia2,laspk,d,uia1,uia2,uaspk,uplevs,info)
+           & l1,l2,lja,lirp,lval,d,uja,uirp,uval,uplevs,info)
       if (info /= psb_success_) then
         info=psb_err_internal_error_
         call psb_errpush(info,name,a_err='Copy/factor loop')
@@ -397,7 +405,7 @@ contains
       call psb_errpush(info,name,a_err='Deallocate')
       goto 9999
     end if
-    if (info == psb_success_) call psb_sp_free(trw,info)
+    if (info == psb_success_) call trw%free()
     if (info /= psb_success_) then
       info=psb_err_from_subroutine_
       ch_err='psb_sp_free'
@@ -473,7 +481,7 @@ contains
   !    ktrw    -  integer, input/output.
   !               The index identifying the last entry taken from the
   !               staging buffer trw. See below.
-  !    trw     -  type(psb_dspmat_type), input/output.
+  !    trw     -  type(psb_zspmat_type), input/output.
   !               A staging buffer. If the matrix A is not in CSR format, we use
   !               the psb_sp_getblk routine and store its output in trw; when we 
   !               need to call psb_sp_getblk we do it for a block of rows, and then
@@ -489,7 +497,7 @@ contains
   
   ! Arguments 
     type(psb_zspmat_type), intent(in)    :: a
-    type(psb_zspmat_type), intent(inout) :: trw
+    type(psb_z_coo_sparse_mat), intent(inout) :: trw
     integer, intent(in)                  :: i,m,jmin,jmax
     integer, intent(inout)               :: ktrw,info
     integer, intent(inout)               :: rowlevs(:)
@@ -497,8 +505,8 @@ contains
     type(psb_int_heap), intent(inout)    :: heap
 
   ! Local variables
-    integer               :: k,j,irb,err_act
-    integer, parameter    :: nrb=16
+    integer               :: k,j,irb,err_act,nz
+    integer, parameter    :: nrb=40
     character(len=20), parameter  :: name='iluk_copyin'
     character(len=20)             :: ch_err
 
@@ -507,22 +515,22 @@ contains
     call psb_erractionsave(err_act)
     call psb_init_heap(heap,info) 
 
-    if (psb_toupper(a%fida) == 'CSR') then
-
+    select type (aa=> a%a) 
+    type is (psb_z_csr_sparse_mat) 
       !
       ! Take a fast shortcut if the matrix is stored in CSR format
       !
       
-      do j = a%ia2(i), a%ia2(i+1) - 1
-        k          = a%ia1(j)
+      do j = aa%irp(i), aa%irp(i+1) - 1
+        k          = aa%ja(j)
         if ((jmin<=k).and.(k<=jmax)) then 
-          row(k)     = a%aspk(j)
+          row(k)     = aa%val(j)
           rowlevs(k) = 0
           call psb_insert_heap(k,heap,info)
         end if
       end do
 
-    else
+    class default
 
       !
       ! Otherwise use psb_sp_getblk, slower but able (in principle) of 
@@ -534,7 +542,7 @@ contains
 
       if ((mod(i,nrb) == 1).or.(nrb == 1)) then 
         irb = min(m-i+1,nrb)
-        call psb_sp_getblk(i,a,trw,info,lrw=i+irb-1)
+        call aa%csget(i,i+irb-1,trw,info)
         if (info /= psb_success_) then
           info=psb_err_from_subroutine_
           ch_err='psb_sp_getblk'
@@ -543,19 +551,19 @@ contains
         end if
         ktrw=1
       end if
-      
+      nz = trw%get_nzeros()
       do 
-        if (ktrw > trw%infoa(psb_nnz_)) exit
-        if (trw%ia1(ktrw) > i) exit
-        k          = trw%ia2(ktrw)
+        if (ktrw > nz) exit
+        if (trw%ia(ktrw) > i) exit
+        k          = trw%ja(ktrw)
         if ((jmin<=k).and.(k<=jmax)) then 
-          row(k)     = trw%aspk(ktrw)
+          row(k)     = trw%val(ktrw)
           rowlevs(k) = 0
           call psb_insert_heap(k,heap,info)
         end if
         ktrw       = ktrw + 1
       enddo
-    end if
+    end select
     call psb_erractionrestore(err_act)
     return
 
@@ -611,17 +619,17 @@ contains
   !    d       -  complex(psb_dpk_), input.
   !               The inverse of the diagonal entries of the part of the U factor
   !               above the current row (see iluk_copyout).
-  !    uia1    -  integer, dimension(:), input.
+  !    uja    -  integer, dimension(:), input.
   !               The column indices of the nonzero entries of the part of the U
-  !               factor above the current row, stored in uaspk row by row (see
+  !               factor above the current row, stored in uval row by row (see
   !               iluk_copyout, called by mld_ziluk_factint), according to the CSR
   !               storage format.
-  !    uia2    -  integer, dimension(:), input.
+  !    uirp    -  integer, dimension(:), input.
   !               The indices identifying the first nonzero entry of each row of
-  !               the U factor above the current row, stored in uaspk row by row
+  !               the U factor above the current row, stored in uval row by row
   !               (see iluk_copyout, called by mld_ziluk_factint), according to
   !               the CSR storage format.
-  !    uaspk   -  complex(psb_dpk_), dimension(:), input.
+  !    uval   -  complex(psb_dpk_), dimension(:), input.
   !               The entries of the U factor above the current row (except the
   !               diagonal ones), stored according to the CSR format.
   !    uplevs  -  integer, dimension(:), input.
@@ -638,7 +646,7 @@ contains
   !               Note: this argument is intent(inout) and not only intent(out)
   !               to retain its allocation, done by this routine.
   !
-  subroutine iluk_fact(fill_in,i,row,rowlevs,heap,d,uia1,uia2,uaspk,uplevs,nidx,idxs,info)
+  subroutine iluk_fact(fill_in,i,row,rowlevs,heap,d,uja,uirp,uval,uplevs,nidx,idxs,info)
 
     use psb_sparse_mod
 
@@ -650,8 +658,8 @@ contains
     integer, intent(inout)               :: nidx,info
     integer, intent(inout)               :: rowlevs(:)
     integer, allocatable, intent(inout)  :: idxs(:)
-    integer, intent(inout)               :: uia1(:),uia2(:),uplevs(:)
-    complex(psb_dpk_), intent(inout)   :: row(:), uaspk(:),d(:)
+    integer, intent(inout)               :: uja(:),uirp(:),uplevs(:)
+    complex(psb_dpk_), intent(inout)        :: row(:), uval(:),d(:)
 
     ! Local variables
     integer               :: k,j,lrwk,jj,lastk, iret
@@ -694,8 +702,8 @@ contains
         row(k) = row(k) * d(k)    ! d(k) == 1/a(k,k)          
         lrwk   = rowlevs(k)
           
-        do jj=uia2(k),uia2(k+1)-1
-          j = uia1(jj)
+        do jj=uirp(k),uirp(k+1)-1
+          j = uja(jj)
           if (j<=k) then 
             info = -i
             return
@@ -715,7 +723,7 @@ contains
           !
           ! Update row(j) and the corresponding fill level
           !
-          row(j)     = row(j) - rwk * uaspk(jj)
+          row(j)     = row(j) - rwk * uval(jj)
           rowlevs(j) = min(rowlevs(j),lrwk+uplevs(jj)+1)
         end do
 
@@ -730,19 +738,19 @@ contains
   ! Note: internal subroutine of mld_ziluk_fact
   !
   !  This routine copies a matrix row, computed by iluk_fact by applying an
-  !  elimination step of the ILU(k) factorization, into the arrays laspk, uaspk,
+  !  elimination step of the ILU(k) factorization, into the arrays lval, uval,
   !  d, corresponding to the L factor, the U factor and the diagonal of U,
   !  respectively.
   !
   !  Note that
-  !  - the part of the row stored into uaspk is scaled by the corresponding diagonal
+  !  - the part of the row stored into uval is scaled by the corresponding diagonal
   !    entry, according to the LDU form of the incomplete factorization;
   !  - the inverse of the diagonal entries of U is actually stored into d; this is
   !    then managed in the solve stage associated to the ILU(k)/MILU(k) factorization;
   !  - if the MILU(k) factorization has been required (ialg == mld_milu_n_), the
   !    row entries discarded because their fill levels are too high are added to
   !    the diagonal entry of the row;
-  !  - the row entries are stored in laspk and uaspk according to the CSR format;
+  !  - the row entries are stored in lval and uval according to the CSR format;
   !  - the arrays row and rowlevs are re-initialized for future use in mld_iluk_fact
   !    (see also iluk_copyin and iluk_fact).
   !
@@ -780,32 +788,32 @@ contains
   !               examined during the elimination step carried out by the routine
   !               iluk_fact.
   !    l1      -  integer, input/output.
-  !               Pointer to the last occupied entry of laspk.
+  !               Pointer to the last occupied entry of lval.
   !    l2      -  integer, input/output.
-  !               Pointer to the last occupied entry of uaspk.
-  !    lia1    -  integer, dimension(:), input/output.
+  !               Pointer to the last occupied entry of uval.
+  !    lja    -  integer, dimension(:), input/output.
   !               The column indices of the nonzero entries of the L factor,
-  !               copied in laspk row by row (see mld_ziluk_factint), according
+  !               copied in lval row by row (see mld_ziluk_factint), according
   !               to the CSR storage format.
-  !    lia2    -  integer, dimension(:), input/output.
+  !    lirp    -  integer, dimension(:), input/output.
   !               The indices identifying the first nonzero entry of each row
-  !               of the L factor, copied in laspk row by row (see 
+  !               of the L factor, copied in lval row by row (see 
   !               mld_ziluk_factint), according to the CSR storage format.
-  !    laspk   -  complex(psb_dpk_), dimension(:), input/output.
+  !    lval   -  complex(psb_dpk_), dimension(:), input/output.
   !               The array where the entries of the row corresponding to the
   !               L factor are copied.
   !    d       -  complex(psb_dpk_), dimension(:), input/output.
   !               The array where the inverse of the diagonal entry of the
   !               row is copied (only d(i) is used by the routine). 
-  !    uia1    -  integer, dimension(:), input/output.
+  !    uja    -  integer, dimension(:), input/output.
   !               The column indices of the nonzero entries of the U factor
-  !               copied in uaspk row by row (see mld_ziluk_factint), according
+  !               copied in uval row by row (see mld_ziluk_factint), according
   !               to the CSR storage format.
-  !    uia2    -  integer, dimension(:), input/output.
+  !    uirp    -  integer, dimension(:), input/output.
   !               The indices identifying the first nonzero entry of each row
-  !               of the U factor copied in uaspk row by row (see
+  !               of the U factor copied in uval row by row (see
   !               mld_zilu_fctint), according to the CSR storage format.
-  !    uaspk   -  complex(psb_dpk_), dimension(:), input/output.
+  !    uval   -  complex(psb_dpk_), dimension(:), input/output.
   !               The array where the entries of the row corresponding to the
   !               U factor are copied.
   !    uplevs  -  integer, dimension(:), input.
@@ -813,19 +821,19 @@ contains
   !               U factor above the current row.
   !
   subroutine iluk_copyout(fill_in,ialg,i,m,row,rowlevs,nidx,idxs,&
-       &  l1,l2,lia1,lia2,laspk,d,uia1,uia2,uaspk,uplevs,info)
+       &  l1,l2,lja,lirp,lval,d,uja,uirp,uval,uplevs,info)
 
     use psb_sparse_mod
 
     implicit none 
 
     ! Arguments
-    integer, intent(in)                             :: fill_in, ialg, i, m, nidx
-    integer, intent(inout)                          :: l1, l2, info
-    integer, intent(inout)                          :: rowlevs(:), idxs(:)
-    integer, allocatable, intent(inout)             :: uia1(:), uia2(:), lia1(:), lia2(:),uplevs(:)
-    complex(psb_dpk_), allocatable, intent(inout) :: uaspk(:), laspk(:)
-    complex(psb_dpk_), intent(inout)              :: row(:), d(:)
+    integer, intent(in)                  :: fill_in, ialg, i, m, nidx
+    integer, intent(inout)               :: l1, l2, info
+    integer, intent(inout)               :: rowlevs(:), idxs(:)
+    integer, allocatable, intent(inout)  :: uja(:), uirp(:), lja(:), lirp(:),uplevs(:)
+    complex(psb_dpk_), allocatable, intent(inout) :: uval(:), lval(:)
+    complex(psb_dpk_), intent(inout)     :: row(:), d(:)
 
     ! Local variables
     integer               :: j,isz,err_act,int_err(5),idxp
@@ -848,21 +856,21 @@ contains
         !  
         if (rowlevs(j) <= fill_in) then 
           l1     = l1 + 1 
-          if (size(laspk) < l1) then
+          if (size(lval) < l1) then
             ! 
             ! Figure out a good reallocation size!
             ! 
             isz  = (max((l1/i)*m,int(1.2*l1),l1+100))
-            call psb_realloc(isz,laspk,info) 
-            if (info == psb_success_) call psb_realloc(isz,lia1,info) 
+            call psb_realloc(isz,lval,info) 
+            if (info == psb_success_) call psb_realloc(isz,lja,info) 
             if (info /= psb_success_) then 
               info=psb_err_from_subroutine_
               call psb_errpush(info,name,a_err='Allocate')
               goto 9999
             end if
           end if
-          lia1(l1)   = j
-          laspk(l1)  = row(j)
+          lja(l1)   = j
+          lval(l1)  = row(j)
         else if (ialg == mld_milu_n_) then
           !
           ! MILU(k): add discarded entries to the diagonal one
@@ -890,13 +898,13 @@ contains
         ! 
         if (rowlevs(j) <= fill_in) then 
           l2     = l2 + 1 
-          if (size(uaspk) < l2) then 
+          if (size(uval) < l2) then 
             ! 
             ! Figure out a good reallocation size!
             !
             isz  = max((l2/i)*m,int(1.2*l2),l2+100)
-            call psb_realloc(isz,uaspk,info) 
-            if (info == psb_success_) call psb_realloc(isz,uia1,info) 
+            call psb_realloc(isz,uval,info) 
+            if (info == psb_success_) call psb_realloc(isz,uja,info) 
             if (info == psb_success_) call psb_realloc(isz,uplevs,info,pad=(m+1))
             if (info /= psb_success_) then 
               info=psb_err_from_subroutine_
@@ -904,8 +912,8 @@ contains
               goto 9999
             end if
           end if
-          uia1(l2)   = j
-          uaspk(l2)  = row(j)
+          uja(l2)   = j
+          uval(l2)  = row(j)
           uplevs(l2) = rowlevs(j) 
         else if (ialg == mld_milu_n_) then
           !
@@ -923,10 +931,10 @@ contains
 
     !
     ! Store the pointers to the first non occupied entry of in
-    ! laspk and uaspk
+    ! lval  and uval
     !
-    lia2(i+1) = l1 + 1
-    uia2(i+1) = l2 + 1
+    lirp(i+1) = l1 + 1
+    uirp(i+1) = l2 + 1
 
     !     
     ! Check the pivot size
@@ -950,8 +958,8 @@ contains
     !
     ! Scale the upper part
     !
-    do j=uia2(i), uia2(i+1)-1
-      uaspk(j) = d(i)*uaspk(j)
+    do j=uirp(i), uirp(i+1)-1
+      uval(j) = d(i)*uval(j)
     end do
 
     call psb_erractionrestore(err_act)
diff --git a/mlprec/mld_zilut_fact.f90 b/mlprec/mld_zilut_fact.f90
index bb50d159..a01d8879 100644
--- a/mlprec/mld_zilut_fact.f90
+++ b/mlprec/mld_zilut_fact.f90
@@ -95,7 +95,7 @@
 subroutine mld_zilut_fact(fill_in,thres,a,l,u,d,info,blck)
   
   use psb_sparse_mod
-  use mld_inner_mod, mld_protect_name => mld_zilut_fact
+  use mld_inner_mod!, mld_protect_name => mld_zilut_fact
 
   implicit none
 
@@ -112,6 +112,7 @@ subroutine mld_zilut_fact(fill_in,thres,a,l,u,d,info,blck)
   integer   ::  l1, l2, m, err_act
   
   type(psb_zspmat_type), pointer  :: blck_
+  type(psb_z_csr_sparse_mat)       :: ll, uu
   character(len=20)   :: name, ch_err
 
   name='mld_zilut_fact'
@@ -130,26 +131,32 @@ subroutine mld_zilut_fact(fill_in,thres,a,l,u,d,info,blck)
     blck_ => blck
   else
     allocate(blck_,stat=info) 
-    if (info /= psb_success_) then 
-      call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
-      goto 9999      
-    end if
-
-    call psb_sp_all(0,0,blck_,1,info)
+    if (info == psb_success_) call blck_%csall(0,0,info,1) 
     if (info /= psb_success_) then
-       info=psb_err_from_subroutine_
-       ch_err='psb_sp_all'
-       call psb_errpush(info,name,a_err=ch_err)
-       goto 9999
+      info=psb_err_from_subroutine_
+      ch_err='csall'
+      call psb_errpush(info,name,a_err=ch_err)
+      goto 9999
     end if
-
   endif
 
+  m = a%get_nrows() + blck_%get_nrows()
+  if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.&
+       & (m > size(d))    ) then 
+    write(0,*) 'Wrong allocation status for L,D,U? ',&
+         & l%get_nrows(),size(d),u%get_nrows()
+    info = -1
+    return
+  end if
+
+  call l%mv_to(ll)
+  call u%mv_to(uu)
+
   !
   ! Compute the ILU(k,t) factorization
   !
-  call mld_zilut_factint(fill_in,thres,m,a,blck_,&
-       & d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info)
+  call mld_zilut_factint(fill_in,thres,a,blck_,&
+       & d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,info)
   if (info /= psb_success_) then
      info=psb_err_from_subroutine_
      ch_err='mld_zilut_factint'
@@ -160,31 +167,29 @@ subroutine mld_zilut_fact(fill_in,thres,a,l,u,d,info,blck)
   !
   ! Store information on the L and U sparse matrices
   !
-  l%infoa(1) = l1
-  l%fida     = 'CSR'
-  l%descra   = 'TLU'
-  u%infoa(1) = l2
-  u%fida     = 'CSR'
-  u%descra   = 'TUU'
-  l%m = m
-  l%k = m
-  u%m = m
-  u%k = m
-
+  call l%mv_from(ll)
+  call l%set_triangle()
+  call l%set_unit()
+  call l%set_lower()
+  call u%mv_from(uu)
+  call u%set_triangle()
+  call u%set_unit()
+  call u%set_upper()
+  
   !
-  ! Nullify the pointer / deallocate the memory
+  ! Nullify pointer / deallocate memory
   !
   if (present(blck)) then 
     blck_ => null() 
   else
-    call psb_sp_free(blck_,info)
-    if (info /= psb_success_) then
-       info=psb_err_from_subroutine_
-       ch_err='psb_sp_free'
-       call psb_errpush(info,name,a_err=ch_err)
-       goto 9999
+    call blck_%free()
+    deallocate(blck_,stat=info) 
+    if(info.ne.0) then
+      info=psb_err_from_subroutine_
+      ch_err='psb_sp_free'
+      call psb_errpush(info,name,a_err=ch_err)
+      goto 9999
     end if
-    deallocate(blck_) 
   endif
 
   call psb_erractionrestore(err_act)
@@ -241,53 +246,53 @@ contains
   !    d       -  complex(psb_dpk_), dimension(:), output.
   !               The inverse of the diagonal entries of the U factor in the incomplete
   !               factorization.
-  !    laspk   -  complex(psb_dpk_), dimension(:), input/output.
+  !    lval   -  complex(psb_dpk_), dimension(:), input/output.
   !               The L factor in the incomplete factorization.
   !    lia1    -  integer, dimension(:), input/output.
   !               The column indices of the nonzero entries of the L factor,
   !               according to the CSR storage format.
-  !    lia2    -  integer, dimension(:), input/output.
+  !    lirp    -  integer, dimension(:), input/output.
   !               The indices identifying the first nonzero entry of each row
-  !               of the L factor in laspk, according to the CSR storage format. 
-  !    uaspk   -  complex(psb_dpk_), dimension(:), input/output.
+  !               of the L factor in lval, according to the CSR storage format. 
+  !    uval   -  complex(psb_dpk_), dimension(:), input/output.
   !               The U factor in the incomplete factorization.
   !               The entries of U are stored according to the CSR format.
-  !    uia1    -  integer, dimension(:), input/output.
+  !    uja    -  integer, dimension(:), input/output.
   !               The column indices of the nonzero entries of the U factor,
   !               according to the CSR storage format.
-  !    uia2    -  integer, dimension(:), input/output.
+  !    uirp    -  integer, dimension(:), input/output.
   !               The indices identifying the first nonzero entry of each row
-  !               of the U factor in uaspk, according to the CSR storage format. 
+  !               of the U factor in uval, according to the CSR storage format. 
   !    l1      -  integer, output
-  !               The number of nonzero entries in laspk.
+  !               The number of nonzero entries in lval.
   !    l2      -  integer, output
-  !               The number of nonzero entries in uaspk.
+  !               The number of nonzero entries in uval.
   !    info    -  integer, output.           
   !               Error code.
   !
-  subroutine mld_zilut_factint(fill_in,thres,m,a,b,&
-       & d,laspk,lia1,lia2,uaspk,uia1,uia2,l1,l2,info)
+  subroutine mld_zilut_factint(fill_in,thres,a,b,&
+       & d,lval,lja,lirp,uval,uja,uirp,l1,l2,info)
 
     use psb_sparse_mod
 
     implicit none 
 
   ! Arguments
-    integer, intent(in)                           :: fill_in
-    real(psb_dpk_), intent(in)                    :: thres
-    type(psb_zspmat_type), intent(in)             :: a,b
-    integer, intent(inout)                        :: m,l1,l2,info
-    integer, allocatable, intent(inout)           :: lia1(:),lia2(:),uia1(:),uia2(:)
-    complex(psb_dpk_), allocatable, intent(inout) :: laspk(:),uaspk(:)
+    integer, intent(in)                        :: fill_in
+    real(psb_dpk_), intent(in)                 :: thres
+    type(psb_zspmat_type),intent(in)           :: a,b
+    integer,intent(inout)                      :: l1,l2,info
+    integer, allocatable, intent(inout)        :: lja(:),lirp(:),uja(:),uirp(:)
+    complex(psb_dpk_), allocatable, intent(inout) :: lval(:),uval(:)
     complex(psb_dpk_), intent(inout)              :: d(:)
 
     ! Local Variables
-    integer :: i, ktrw,err_act,nidx,nlw,nup,jmaxup, ma, mb
+    integer :: i, ktrw,err_act,nidx,nlw,nup,jmaxup, ma, mb, m
     real(psb_dpk_) ::  nrmi
     integer, allocatable          :: idxs(:)
     complex(psb_dpk_), allocatable :: row(:)
     type(psb_int_heap) :: heap
-    type(psb_zspmat_type) :: trw
+    type(psb_z_coo_sparse_mat) :: trw
     character(len=20), parameter  :: name='mld_zilut_factint'
     character(len=20)             :: ch_err
 
@@ -296,16 +301,16 @@ contains
     call psb_erractionsave(err_act)
 
 
-    ma = a%m
-    mb = b%m
+    ma = a%get_nrows()
+    mb = b%get_nrows()
     m  = ma+mb
 
     !
     ! Allocate a temporary buffer for the ilut_copyin function 
     !
-    call psb_sp_all(0,0,trw,1,info)
-    if (info == psb_success_) call psb_ensure_size(m+1,lia2,info)
-    if (info == psb_success_) call psb_ensure_size(m+1,uia2,info)
+    call trw%allocate(0,0,1)
+    if (info == psb_success_) call psb_ensure_size(m+1,lirp,info)
+    if (info == psb_success_) call psb_ensure_size(m+1,uirp,info)
 
     if (info /= psb_success_) then
       info=psb_err_from_subroutine_
@@ -315,8 +320,8 @@ contains
     
     l1=0
     l2=0
-    lia2(1) = 1
-    uia2(1) = 1
+    lirp(1) = 1
+    uirp(1) = 1
 
     !
     ! Allocate memory to hold the entries of a row
@@ -354,12 +359,12 @@ contains
       ! Do an elimination step on current row
       !
       if (info == psb_success_) call ilut_fact(thres,i,nrmi,row,heap,&
-           & d,uia1,uia2,uaspk,nidx,idxs,info)
+           & d,uja,uirp,uval,nidx,idxs,info)
       !
-      ! Copy the row into laspk/d(i)/uaspk
+      ! Copy the row into lval/d(i)/uval
       ! 
       if (info == psb_success_) call ilut_copyout(fill_in,thres,i,m,nlw,nup,jmaxup,nrmi,row,nidx,idxs,&
-           & l1,l2,lia1,lia2,laspk,d,uia1,uia2,uaspk,info)
+           & l1,l2,lja,lirp,lval,d,uja,uirp,uval,info)
 
       if (info /= psb_success_) then
         info=psb_err_internal_error_
@@ -378,7 +383,7 @@ contains
       call psb_errpush(info,name,a_err='Deallocate')
       goto 9999
     end if
-    if (info == psb_success_) call psb_sp_free(trw,info)
+    if (info == psb_success_) call trw%free()
     if (info /= psb_success_) then
       info=psb_err_from_subroutine_
       ch_err='psb_sp_free'
@@ -482,17 +487,17 @@ contains
   subroutine ilut_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,nrmi,row,heap,ktrw,trw,info)
     use psb_sparse_mod
     implicit none 
-    type(psb_zspmat_type), intent(in)    :: a
-    type(psb_zspmat_type), intent(inout) :: trw
-    integer, intent(in)                  :: i, m,jmin,jmax,jd
-    integer, intent(inout)               :: ktrw,nlw,nup,jmaxup,info
-    real(psb_dpk_), intent(inout)        :: nrmi
-    complex(psb_dpk_), intent(inout)     :: row(:)
-    type(psb_int_heap), intent(inout)    :: heap
+    type(psb_zspmat_type), intent(in)         :: a
+    type(psb_z_coo_sparse_mat), intent(inout) :: trw
+    integer, intent(in)                       :: i, m,jmin,jmax,jd
+    integer, intent(inout)                    :: ktrw,nlw,nup,jmaxup,info
+    real(psb_dpk_), intent(inout)             :: nrmi
+    complex(psb_dpk_), intent(inout)          :: row(:)
+    type(psb_int_heap), intent(inout)         :: heap
     
     integer               :: k,j,irb,kin,nz
-    integer, parameter    :: nrb=16
-    real(psb_dpk_)      :: dmaxup
+    integer, parameter    :: nrb=40
+    real(psb_dpk_)        :: dmaxup
     real(psb_dpk_), external    :: dznrm2
     character(len=20), parameter  :: name='mld_zilut_factint'
 
@@ -517,25 +522,21 @@ contains
     nlw    = 0
     nup    = 0
     jmaxup = 0
-    dmaxup = dzero
-    nrmi   = dzero
-
-    if (psb_toupper(a%fida) == 'CSR') then
-
+    dmaxup = szero
+    nrmi   = szero
+    
+    select type (aa=> a%a) 
+    type is (psb_z_csr_sparse_mat) 
       !
       ! Take a fast shortcut if the matrix is stored in CSR format
-      ! 
-
-      do j = a%ia2(i), a%ia2(i+1) - 1
-        k          = a%ia1(j)
+      !
+      
+      do j = aa%irp(i), aa%irp(i+1) - 1
+        k          = aa%ja(j)
         if ((jmin<=k).and.(k<=jmax)) then 
-          row(k)     = a%aspk(j)
+          row(k)     = aa%val(j)
           call psb_insert_heap(k,heap,info)
-          if (info /= psb_success_) then
-            info=psb_err_from_subroutine_
-            call psb_errpush(info,name,a_err='psb_insert_heap')
-            goto 9999
-          end if
+          if (info /= psb_success_) exit
         end if
         if (k<jd) nlw = nlw + 1 
         if (k>jd) then 
@@ -546,9 +547,17 @@ contains
           end if
         end if
       end do
-      nz   = a%ia2(i+1) - a%ia2(i)
-      nrmi = dznrm2(nz,a%aspk(a%ia2(i)),ione)
-    else
+      if (info /= psb_success_) then 
+        info=psb_err_from_subroutine_
+        call psb_errpush(info,name,a_err='psb_insert_heap')
+        goto 9999
+      end if
+      
+      nz   = aa%irp(i+1) - aa%irp(i)
+      nrmi = dznrm2(nz,aa%val(aa%irp(i)),ione)
+
+
+    class default
 
       !
       ! Otherwise use psb_sp_getblk, slower but able (in principle) of 
@@ -560,7 +569,7 @@ contains
 
       if ((mod(i,nrb) == 1).or.(nrb == 1)) then 
         irb = min(m-i+1,nrb)
-        call psb_sp_getblk(i,a,trw,info,lrw=i+irb-1)
+        call aa%csget(i,i+irb-1,trw,info)
         if (info /= psb_success_) then
           info=psb_err_from_subroutine_
           call psb_errpush(info,name,a_err='psb_sp_getblk')
@@ -570,18 +579,16 @@ contains
       end if
       
       kin = ktrw
+      nz = trw%get_nzeros()
       do 
-        if (ktrw > trw%infoa(psb_nnz_)) exit
-        if (trw%ia1(ktrw) > i) exit
-        k          = trw%ia2(ktrw)
+        if (ktrw > nz) exit
+        if (trw%ia(ktrw) > i) exit
+        k          = trw%ja(ktrw)
         if ((jmin<=k).and.(k<=jmax)) then 
-          row(k)     = trw%aspk(ktrw)
+          row(k)     = trw%val(ktrw)
           call psb_insert_heap(k,heap,info)
-          if (info /= psb_success_) then
-            info=psb_err_from_subroutine_
-            call psb_errpush(info,name,a_err='psb_insert_heap')
-            goto 9999
-          end if
+          if (info /= psb_success_) exit
+          
         end if
         if (k<jd) nlw = nlw + 1 
         if (k>jd) then 
@@ -594,8 +601,9 @@ contains
         ktrw       = ktrw + 1
       enddo
       nz = ktrw - kin
-      nrmi = dznrm2(nz,trw%aspk(kin),ione)
-    end if
+      nrmi = dznrm2(nz,trw%val(kin),ione)
+    end select
+
     call psb_erractionrestore(err_act)
     return
 
@@ -645,17 +653,17 @@ contains
   !    d       -  complex(psb_dpk_), input.
   !               The inverse of the diagonal entries of the part of the U factor
   !               above the current row (see ilut_copyout).
-  !    uia1    -  integer, dimension(:), input.
+  !    uja    -  integer, dimension(:), input.
   !               The column indices of the nonzero entries of the part of the U
-  !               factor above the current row, stored in uaspk row by row (see
+  !               factor above the current row, stored in uval row by row (see
   !               ilut_copyout, called by mld_zilut_factint), according to the CSR
   !               storage format.
-  !    uia2    -  integer, dimension(:), input.
+  !    uirp    -  integer, dimension(:), input.
   !               The indices identifying the first nonzero entry of each row of
-  !               the U factor above the current row, stored in uaspk row by row
+  !               the U factor above the current row, stored in uval row by row
   !               (see ilut_copyout, called by mld_zilut_factint), according to
   !               the CSR storage format.
-  !    uaspk   -  complex(psb_dpk_), dimension(:), input.
+  !    uval   -  complex(psb_dpk_), dimension(:), input.
   !               The entries of the U factor above the current row (except the
   !               diagonal ones), stored according to the CSR format.
   !    nidx    -  integer, output.
@@ -669,7 +677,7 @@ contains
   !               Note: this argument is intent(inout) and not only intent(out)
   !               to retain its allocation, done by this routine.
   !
-  subroutine ilut_fact(thres,i,nrmi,row,heap,d,uia1,uia2,uaspk,nidx,idxs,info)
+  subroutine ilut_fact(thres,i,nrmi,row,heap,d,uja,uirp,uval,nidx,idxs,info)
 
     use psb_sparse_mod
 
@@ -681,8 +689,8 @@ contains
     integer, intent(inout)              :: nidx,info
     real(psb_dpk_), intent(in)          :: thres,nrmi
     integer, allocatable, intent(inout) :: idxs(:)
-    integer, intent(inout)              :: uia1(:),uia2(:)
-    complex(psb_dpk_), intent(inout)    :: row(:), uaspk(:),d(:)
+    integer, intent(inout)              :: uja(:),uirp(:)
+    complex(psb_dpk_), intent(inout)    :: row(:), uval(:),d(:)
 
     ! Local Variables
     integer               :: k,j,jj,lastk, iret
@@ -726,8 +734,8 @@ contains
           ! Note: since U is scaled while copying it out (see ilut_copyout),
           ! we can use rwk in the update below.
           !           
-          do jj=uia2(k),uia2(k+1)-1
-            j = uia1(jj)
+          do jj=uirp(k),uirp(k+1)-1
+            j = uja(jj)
             if (j<=k) then 
               info = -i 
               return
@@ -736,7 +744,7 @@ contains
             ! Update row(j) and, if it is not to be discarded, insert
             ! its index into the heap for further processing.
             !
-            row(j)     = row(j) - rwk * uaspk(jj)
+            row(j)     = row(j) - rwk * uval(jj)
             if (abs(row(j)) < thres*nrmi) then
               ! 
               ! Drop the entry.
@@ -771,8 +779,8 @@ contains
   ! Note: internal subroutine of mld_zilut_fact
   !
   !  This routine copies a matrix row, computed by ilut_fact by applying an
-  !  elimination step of the ILU(k,t) factorization, into the arrays laspk,
-  !  uaspk, d, corresponding to the L factor, the U factor and the diagonal
+  !  elimination step of the ILU(k,t) factorization, into the arrays lval,
+  !  uval, d, corresponding to the L factor, the U factor and the diagonal
   !  of U, respectively.
   !
   !  Note that
@@ -781,11 +789,11 @@ contains
   !    the 'lower part' of the row, and the nup+k ones in the 'upper part';
   !  - the entry in the upper part of the row which has maximum absolute value
   !    in the original matrix is included in the above nup+k entries anyway;
-  !  - the part of the row stored into uaspk is scaled by the corresponding
+  !  - the part of the row stored into uval is scaled by the corresponding
   !    diagonal entry, according to the LDU form of the incomplete factorization;
   !  - the inverse of the diagonal entries of U is actually stored into d; this
   !    is then managed in the solve stage associated to the ILU(k,t) factorization;
-  !  - the row entries are stored in laspk and uaspk according to the CSR format;
+  !  - the row entries are stored in lval and uval according to the CSR format;
   !  - the array row is re-initialized for future use in mld_ilut_fact(see also
   !    ilut_copyin and ilut_fact).
   !
@@ -825,49 +833,49 @@ contains
   !               examined during the elimination step carried out by the routine
   !               ilut_fact.
   !    l1      -  integer, input/output.
-  !               Pointer to the last occupied entry of laspk.
+  !               Pointer to the last occupied entry of lval.
   !    l2      -  integer, input/output.
-  !               Pointer to the last occupied entry of uaspk.
-  !    lia1    -  integer, dimension(:), input/output.
+  !               Pointer to the last occupied entry of uval.
+  !    lja    -  integer, dimension(:), input/output.
   !               The column indices of the nonzero entries of the L factor,
-  !               copied in laspk row by row (see mld_zilut_factint), according
+  !               copied in lval row by row (see mld_zilut_factint), according
   !               to the CSR storage format.
-  !    lia2    -  integer, dimension(:), input/output.
+  !    lirp    -  integer, dimension(:), input/output.
   !               The indices identifying the first nonzero entry of each row
-  !               of the L factor, copied in laspk row by row (see 
+  !               of the L factor, copied in lval row by row (see 
   !               mld_zilut_factint), according to the CSR storage format.
-  !    laspk   -  complex(psb_dpk_), dimension(:), input/output.
+  !    lval   -  complex(psb_dpk_), dimension(:), input/output.
   !               The array where the entries of the row corresponding to the
   !               L factor are copied.
   !    d       -  complex(psb_dpk_), dimension(:), input/output.
   !               The array where the inverse of the diagonal entry of the
   !               row is copied (only d(i) is used by the routine). 
-  !    uia1    -  integer, dimension(:), input/output.
+  !    uja    -  integer, dimension(:), input/output.
   !               The column indices of the nonzero entries of the U factor
-  !               copied in uaspk row by row (see mld_zilut_factint), according
+  !               copied in uval row by row (see mld_zilut_factint), according
   !               to the CSR storage format.
-  !    uia2    -  integer, dimension(:), input/output.
+  !    uirp    -  integer, dimension(:), input/output.
   !               The indices identifying the first nonzero entry of each row
-  !               of the U factor copied in uaspk row by row (see
+  !               of the U factor copied in uval row by row (see
   !               mld_zilu_fctint), according to the CSR storage format.
-  !    uaspk   -  complex(psb_dpk_), dimension(:), input/output.
+  !    uval   -  complex(psb_dpk_), dimension(:), input/output.
   !               The array where the entries of the row corresponding to the
   !               U factor are copied.
   !
   subroutine ilut_copyout(fill_in,thres,i,m,nlw,nup,jmaxup,nrmi,row, &
-       & nidx,idxs,l1,l2,lia1,lia2,laspk,d,uia1,uia2,uaspk,info)
+       & nidx,idxs,l1,l2,lja,lirp,lval,d,uja,uirp,uval,info)
 
     use psb_sparse_mod
 
     implicit none 
 
     ! Arguments
-    integer, intent(in)                          :: fill_in,i,m,nidx,nlw,nup,jmaxup
-    integer, intent(in)                          :: idxs(:)
-    integer, intent(inout)                       :: l1,l2, info
-    integer, allocatable, intent(inout)          :: uia1(:),uia2(:), lia1(:),lia2(:)
-    real(psb_dpk_), intent(in)                   :: thres,nrmi
-    complex(psb_dpk_),allocatable, intent(inout) :: uaspk(:), laspk(:)
+    integer, intent(in)                       :: fill_in,i,m,nidx,nlw,nup,jmaxup
+    integer, intent(in)                       :: idxs(:)
+    integer, intent(inout)                    :: l1,l2, info
+    integer, allocatable, intent(inout)       :: uja(:),uirp(:), lja(:),lirp(:)
+    real(psb_dpk_), intent(in)                :: thres,nrmi
+    complex(psb_dpk_),allocatable, intent(inout) :: uval(:), lval(:)
     complex(psb_dpk_), intent(inout)             :: row(:), d(:)
 
     ! Local variables
@@ -966,21 +974,21 @@ contains
     !
     do k=1,nz
       l1     = l1 + 1 
-      if (size(laspk) < l1) then
+      if (size(lval) < l1) then
         ! 
         ! Figure out a good reallocation size!
         ! 
         isz  = (max((l1/i)*m,int(1.2*l1),l1+100))
-        call psb_realloc(isz,laspk,info) 
-        if (info == psb_success_) call psb_realloc(isz,lia1,info) 
+        call psb_realloc(isz,lval,info) 
+        if (info == psb_success_) call psb_realloc(isz,lja,info) 
         if (info /= psb_success_) then 
           info=psb_err_from_subroutine_
           call psb_errpush(info,name,a_err='Allocate')
           goto 9999
         end if
       end if
-      lia1(l1)   = xwid(k)
-      laspk(l1)  = xw(indx(k))
+      lja(l1)   = xwid(k)
+      lval(l1)  = xw(indx(k))
     end do
     
     !
@@ -1022,7 +1030,7 @@ contains
           !
           ! Compute 1/pivot
           !
-          d(i) = done/d(i)
+          d(i) = zone/d(i)
         end if
       end if
     end if
@@ -1112,21 +1120,21 @@ contains
     !
     do k=1,nz
       l2     = l2 + 1 
-      if (size(uaspk) < l2) then
+      if (size(uval) < l2) then
         ! 
         ! Figure out a good reallocation size!
         ! 
         isz  = max((l2/i)*m,int(1.2*l2),l2+100)
-        call psb_realloc(isz,uaspk,info) 
-        if (info == psb_success_) call psb_realloc(isz,uia1,info) 
+        call psb_realloc(isz,uval,info) 
+        if (info == psb_success_) call psb_realloc(isz,uja,info) 
         if (info /= psb_success_) then 
           info=psb_err_from_subroutine_
           call psb_errpush(info,name,a_err='Allocate')
           goto 9999
         end if
       end if
-      uia1(l2)   = xwid(k)
-      uaspk(l2)  = d(i)*xw(indx(k))
+      uja(l2)   = xwid(k)
+      uval(l2)  = d(i)*xw(indx(k))
     end do
 
     !
@@ -1138,10 +1146,10 @@ contains
 
     !
     ! Store the pointers to the first non occupied entry of in
-    ! laspk and uaspk
+    ! lval and uval
     !
-    lia2(i+1) = l1 + 1
-    uia2(i+1) = l2 + 1
+    lirp(i+1) = l1 + 1
+    uirp(i+1) = l2 + 1
 
     call psb_erractionrestore(err_act)
     return