diff --git a/krylov/psb_ckrylovsubspace_mod.F90 b/krylov/psb_ckrylovsubspace_mod.F90
index ee8ae829..dd50b1ae 100644
--- a/krylov/psb_ckrylovsubspace_mod.F90
+++ b/krylov/psb_ckrylovsubspace_mod.F90
@@ -44,6 +44,8 @@ module psb_ckrylovsubspace_mod
   use psb_util_mod
   use psb_krylov_mod
   implicit none
+  private
+  public psb_cpolkrylov
 
   type :: psb_cpolkrylov
     type(psb_c_vect_type), allocatable, dimension(:) :: v
@@ -180,7 +182,8 @@ contains
     type(psb_ctxt_type) :: ctxt
     integer(psb_ipk_) :: err_act, info
     character(len=20) :: name, ch_err
-    real(psb_spk_) :: nrm0,scal
+    real(psb_spk_) :: nrm0
+    complex(psb_spk_) :: scal
 
     info = psb_success_
     name = 'psb_carnoldi_kryl'
@@ -251,7 +254,8 @@ contains
         goto 9999
       end if
       scal = cone/kryl%h(i1,i)
-      call psb_gescal(kryl%v(i1),scal,kryl%v(i1),desc_a,info)
+      ! call psb_gescal(kryl%v(i1),scal,kryl%v(i1),desc_a,info)
+      call psb_geaxpby(scal,kryl%v(i1),czero,kryl%v(i1),desc_a,info)
       if (info /= psb_success_) then
         info=psb_err_from_subroutine_non_
         call psb_errpush(info,name)
@@ -295,7 +299,7 @@ contains
     call psb_info(ctxt,iam,np)
 
     if (iam == psb_root_) then
-      call mm_array_write(kryl%h,"Projected Matrix",info,filename=filename//"_h")
+      call mm_array_write(kryl%h,"Projected Matrix",info,filename=trim(filename)//"_h")
       if (info /= psb_success_) then
         info=psb_err_from_subroutine_non_
         call psb_errpush(info,name)
@@ -315,7 +319,7 @@ contains
       vglobal(:,i) = w
     end do
     if (iam == psb_root_) then
-      call mm_array_write(vglobal,"Krylov basis",info,filename=filename//"_v")
+      call mm_array_write(vglobal,"Krylov basis",info,filename=trim(filename)//"_v")
       if (info /= psb_success_) then
         info=psb_err_from_subroutine_non_
         call psb_errpush(info,name)
diff --git a/krylov/psb_dkrylovsubspace_mod.F90 b/krylov/psb_dkrylovsubspace_mod.F90
index 9cc28b16..50b4c418 100644
--- a/krylov/psb_dkrylovsubspace_mod.F90
+++ b/krylov/psb_dkrylovsubspace_mod.F90
@@ -44,6 +44,8 @@ module psb_dkrylovsubspace_mod
   use psb_util_mod
   use psb_krylov_mod
   implicit none
+  private
+  public psb_dpolkrylov
 
   type :: psb_dpolkrylov
     type(psb_d_vect_type), allocatable, dimension(:) :: v
@@ -180,7 +182,8 @@ contains
     type(psb_ctxt_type) :: ctxt
     integer(psb_ipk_) :: err_act, info
     character(len=20) :: name, ch_err
-    real(psb_dpk_) :: nrm0,scal
+    real(psb_dpk_) :: nrm0
+    real(psb_dpk_) :: scal
 
     info = psb_success_
     name = 'psb_darnoldi_kryl'
@@ -251,7 +254,8 @@ contains
         goto 9999
       end if
       scal = done/kryl%h(i1,i)
-      call psb_gescal(kryl%v(i1),scal,kryl%v(i1),desc_a,info)
+      ! call psb_gescal(kryl%v(i1),scal,kryl%v(i1),desc_a,info)
+      call psb_geaxpby(scal,kryl%v(i1),dzero,kryl%v(i1),desc_a,info)
       if (info /= psb_success_) then
         info=psb_err_from_subroutine_non_
         call psb_errpush(info,name)
@@ -295,7 +299,7 @@ contains
     call psb_info(ctxt,iam,np)
 
     if (iam == psb_root_) then
-      call mm_array_write(kryl%h,"Projected Matrix",info,filename=filename//"_h")
+      call mm_array_write(kryl%h,"Projected Matrix",info,filename=trim(filename)//"_h")
       if (info /= psb_success_) then
         info=psb_err_from_subroutine_non_
         call psb_errpush(info,name)
@@ -315,7 +319,7 @@ contains
       vglobal(:,i) = w
     end do
     if (iam == psb_root_) then
-      call mm_array_write(vglobal,"Krylov basis",info,filename=filename//"_v")
+      call mm_array_write(vglobal,"Krylov basis",info,filename=trim(filename)//"_v")
       if (info /= psb_success_) then
         info=psb_err_from_subroutine_non_
         call psb_errpush(info,name)
diff --git a/krylov/psb_skrylovsubspace_mod.F90 b/krylov/psb_skrylovsubspace_mod.F90
index 13233c96..25314284 100644
--- a/krylov/psb_skrylovsubspace_mod.F90
+++ b/krylov/psb_skrylovsubspace_mod.F90
@@ -44,6 +44,8 @@ module psb_skrylovsubspace_mod
   use psb_util_mod
   use psb_krylov_mod
   implicit none
+  private
+  public psb_spolkrylov
 
   type :: psb_spolkrylov
     type(psb_s_vect_type), allocatable, dimension(:) :: v
@@ -180,7 +182,8 @@ contains
     type(psb_ctxt_type) :: ctxt
     integer(psb_ipk_) :: err_act, info
     character(len=20) :: name, ch_err
-    real(psb_spk_) :: nrm0,scal
+    real(psb_spk_) :: nrm0
+    real(psb_spk_) :: scal
 
     info = psb_success_
     name = 'psb_sarnoldi_kryl'
@@ -251,7 +254,8 @@ contains
         goto 9999
       end if
       scal = sone/kryl%h(i1,i)
-      call psb_gescal(kryl%v(i1),scal,kryl%v(i1),desc_a,info)
+      ! call psb_gescal(kryl%v(i1),scal,kryl%v(i1),desc_a,info)
+      call psb_geaxpby(scal,kryl%v(i1),szero,kryl%v(i1),desc_a,info)
       if (info /= psb_success_) then
         info=psb_err_from_subroutine_non_
         call psb_errpush(info,name)
@@ -295,7 +299,7 @@ contains
     call psb_info(ctxt,iam,np)
 
     if (iam == psb_root_) then
-      call mm_array_write(kryl%h,"Projected Matrix",info,filename=filename//"_h")
+      call mm_array_write(kryl%h,"Projected Matrix",info,filename=trim(filename)//"_h")
       if (info /= psb_success_) then
         info=psb_err_from_subroutine_non_
         call psb_errpush(info,name)
@@ -315,7 +319,7 @@ contains
       vglobal(:,i) = w
     end do
     if (iam == psb_root_) then
-      call mm_array_write(vglobal,"Krylov basis",info,filename=filename//"_v")
+      call mm_array_write(vglobal,"Krylov basis",info,filename=trim(filename)//"_v")
       if (info /= psb_success_) then
         info=psb_err_from_subroutine_non_
         call psb_errpush(info,name)
diff --git a/krylov/psb_zkrylovsubspace_mod.F90 b/krylov/psb_zkrylovsubspace_mod.F90
index 62930f76..bba262a7 100644
--- a/krylov/psb_zkrylovsubspace_mod.F90
+++ b/krylov/psb_zkrylovsubspace_mod.F90
@@ -44,6 +44,8 @@ module psb_zkrylovsubspace_mod
   use psb_util_mod
   use psb_krylov_mod
   implicit none
+  private
+  public psb_zpolkrylov
 
   type :: psb_zpolkrylov
     type(psb_z_vect_type), allocatable, dimension(:) :: v
@@ -180,7 +182,8 @@ contains
     type(psb_ctxt_type) :: ctxt
     integer(psb_ipk_) :: err_act, info
     character(len=20) :: name, ch_err
-    real(psb_dpk_) :: nrm0,scal
+    real(psb_dpk_) :: nrm0
+    complex(psb_dpk_) :: scal
 
     info = psb_success_
     name = 'psb_zarnoldi_kryl'
@@ -251,7 +254,8 @@ contains
         goto 9999
       end if
       scal = zone/kryl%h(i1,i)
-      call psb_gescal(kryl%v(i1),scal,kryl%v(i1),desc_a,info)
+      ! call psb_gescal(kryl%v(i1),scal,kryl%v(i1),desc_a,info)
+      call psb_geaxpby(scal,kryl%v(i1),zzero,kryl%v(i1),desc_a,info)
       if (info /= psb_success_) then
         info=psb_err_from_subroutine_non_
         call psb_errpush(info,name)
@@ -295,7 +299,7 @@ contains
     call psb_info(ctxt,iam,np)
 
     if (iam == psb_root_) then
-      call mm_array_write(kryl%h,"Projected Matrix",info,filename=filename//"_h")
+      call mm_array_write(kryl%h,"Projected Matrix",info,filename=trim(filename)//"_h")
       if (info /= psb_success_) then
         info=psb_err_from_subroutine_non_
         call psb_errpush(info,name)
@@ -315,7 +319,7 @@ contains
       vglobal(:,i) = w
     end do
     if (iam == psb_root_) then
-      call mm_array_write(vglobal,"Krylov basis",info,filename=filename//"_v")
+      call mm_array_write(vglobal,"Krylov basis",info,filename=trim(filename)//"_v")
       if (info /= psb_success_) then
         info=psb_err_from_subroutine_non_
         call psb_errpush(info,name)
diff --git a/test/krylovspace/Makefile b/test/krylovspace/Makefile
new file mode 100644
index 00000000..132a2e15
--- /dev/null
+++ b/test/krylovspace/Makefile
@@ -0,0 +1,52 @@
+INSTALLDIR=../..
+INCDIR=$(INSTALLDIR)/include/
+MODDIR=$(INSTALLDIR)/modules/
+include $(INCDIR)/Make.inc.psblas
+#
+# Libraries used
+#
+LIBDIR=$(INSTALLDIR)/lib/
+PSBLAS_LIB= -L$(LIBDIR) -lpsb_util -lpsb_krylov -lpsb_prec -lpsb_base
+LDLIBS=$(PSBLDLIBS)
+
+FINCLUDES=$(FMFLAG)$(MODDIR) $(FMFLAG).
+
+SFOBJS=getp.o  psb_sf_sample.o   
+DFOBJS=getp.o  psb_df_sample.o   
+CFOBJS=getp.o  psb_cf_sample.o   
+ZFOBJS=getp.o  psb_zf_sample.o   
+
+EXEDIR=./runs
+
+all: runsd psb_sf_sample psb_df_sample psb_cf_sample psb_zf_sample 
+
+runsd:
+	(if test ! -d runs ; then mkdir runs; fi)
+
+psb_sf_sample.o psb_df_sample.o psb_cf_sample.o psb_zf_sample.o: getp.o
+
+psb_sf_sample: $(SFOBJS)
+	$(FLINK) $(LOPT) $(SFOBJS) -o psb_sf_sample $(PSBLAS_LIB) $(LDLIBS) 
+	/bin/mv psb_sf_sample $(EXEDIR) 
+psb_df_sample: $(DFOBJS)
+	$(FLINK) $(LOPT) $(DFOBJS) -o psb_df_sample $(PSBLAS_LIB) $(LDLIBS) 
+	/bin/mv psb_df_sample $(EXEDIR) 
+psb_cf_sample: $(CFOBJS)
+	$(FLINK) $(LOPT) $(CFOBJS) -o psb_cf_sample $(PSBLAS_LIB) $(LDLIBS)
+	/bin/mv psb_cf_sample $(EXEDIR) 
+psb_zf_sample: $(ZFOBJS)
+	$(FLINK) $(LOPT) $(ZFOBJS) -o psb_zf_sample $(PSBLAS_LIB) $(LDLIBS)
+	/bin/mv psb_zf_sample $(EXEDIR) 
+
+.f90.o:
+	$(MPFC) $(FCOPT) $(FINCLUDES) $(FDEFINES) -c $<
+
+clean: 
+	/bin/rm -f $(DFOBJS) $(ZFOBJS) $(SFOBJS) $(CFOBJS)\
+        *$(.mod) $(EXEDIR)/psb_*f_sample 
+
+lib:
+	(cd ../../; make library)
+verycleanlib:
+	(cd ../../; make veryclean)
+
diff --git a/test/krylovspace/getp.f90 b/test/krylovspace/getp.f90
new file mode 100644
index 00000000..82bef762
--- /dev/null
+++ b/test/krylovspace/getp.f90
@@ -0,0 +1,281 @@
+!   
+!                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.
+!   
+!    
+Module getp
+  interface get_parms
+    module procedure get_dparms, get_sparms
+  end interface
+
+contains
+  !
+  ! Get iteration parameters from the command line
+  !
+  subroutine  get_dparms(ctxt,mtrx_file,rhs_file,filefmt,kmethd,ptype,part,&
+       & afmt,istopc,itmax,itrace,irst,eps)
+    use psb_base_mod
+    type(psb_ctxt_type) :: ctxt
+    character(len=2)  :: filefmt
+    character(len=40) :: kmethd, mtrx_file, rhs_file, ptype
+    character(len=20) :: part
+    integer(psb_ipk_) :: iret, istopc,itmax,itrace,irst
+    character(len=40) :: charbuf
+    real(psb_dpk_) :: eps
+    character    :: afmt*5
+    integer(psb_ipk_) :: np, iam
+    integer(psb_ipk_) :: inparms(40), ip, inp_unit
+    character(len=1024)   :: filename 
+
+    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 Input Parameters
+      read(inp_unit,*) ip
+      if (ip >= 5) then
+        read(inp_unit,*) mtrx_file
+        read(inp_unit,*) rhs_file
+        read(inp_unit,*) filefmt
+        read(inp_unit,*) kmethd
+        read(inp_unit,*) ptype
+        read(inp_unit,*) afmt
+        read(inp_unit,*) part
+
+
+        call psb_bcast(ctxt,mtrx_file)
+        call psb_bcast(ctxt,rhs_file)
+        call psb_bcast(ctxt,filefmt)
+        call psb_bcast(ctxt,kmethd)
+        call psb_bcast(ctxt,ptype)
+        call psb_bcast(ctxt,afmt)
+        call psb_bcast(ctxt,part)
+
+        if (ip >= 7) then
+          read(inp_unit,*) istopc
+        else
+          istopc=1        
+        endif
+        if (ip >= 8) then
+          read(inp_unit,*) itmax
+        else
+          itmax=500
+        endif
+        if (ip >= 9) then
+          read(inp_unit,*) itrace
+        else
+          itrace=-1
+        endif
+        if (ip >= 10) then
+          read(inp_unit,*) irst
+        else
+          irst  = 1
+        endif
+        if (ip >= 11) then
+          read(inp_unit,*) eps
+        else
+          eps=1.d-6
+        endif
+        inparms(1) = istopc
+        inparms(2) = itmax
+        inparms(3) = itrace
+        inparms(4) = irst
+        call psb_bcast(ctxt,inparms(1:4))
+        call psb_bcast(ctxt,eps)
+
+        write(psb_out_unit,'("Solving matrix       : ",a)')  mtrx_file      
+        write(psb_out_unit,'("Number of processors : ",i3)') np
+        write(psb_out_unit,'("Data distribution    : ",a)') part
+        write(psb_out_unit,'("Iterative method     : ",a)')  kmethd
+        write(psb_out_unit,'("Preconditioner       : ",a)')  ptype
+        write(psb_out_unit,'("Restart parameter    : ",i2)') irst
+        write(psb_out_unit,'("Storage format       : ",a)')  afmt(1:3)
+        write(psb_out_unit,'(" ")')
+      else
+        write(psb_err_unit,*) 'Wrong format for input file'
+        call psb_abort(ctxt)
+        stop 1
+      end if
+      if (inp_unit /= psb_inp_unit) then
+        close(inp_unit)
+      end if
+    else
+      ! Receive Parameters
+      call psb_bcast(ctxt,mtrx_file)
+      call psb_bcast(ctxt,rhs_file)
+      call psb_bcast(ctxt,filefmt)
+      call psb_bcast(ctxt,kmethd)
+      call psb_bcast(ctxt,ptype)
+      call psb_bcast(ctxt,afmt)
+      call psb_bcast(ctxt,part)
+
+      call psb_bcast(ctxt,inparms(1:4))
+      istopc =  inparms(1) 
+      itmax  =  inparms(2) 
+      itrace =  inparms(3) 
+      irst   =  inparms(4) 
+      call psb_bcast(ctxt,eps)
+
+    end if
+
+  end subroutine get_dparms
+  
+  subroutine  get_sparms(ctxt,mtrx_file,rhs_file,filefmt,kmethd,ptype,part,&
+       & afmt,istopc,itmax,itrace,irst,eps)
+    use psb_base_mod
+    type(psb_ctxt_type) :: ctxt
+    character(len=2)  :: filefmt
+    character(len=40) :: kmethd, mtrx_file, rhs_file, ptype
+    character(len=20) :: part
+    integer(psb_ipk_) :: iret, istopc,itmax,itrace,irst
+    character(len=40) :: charbuf
+    real(psb_spk_) :: eps
+    character    :: afmt*5
+    integer(psb_ipk_) :: np, iam
+    integer(psb_ipk_) :: inparms(40), ip, inp_unit
+    character(len=1024)   :: filename
+
+    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=inp_unit
+      end if
+      ! Read Input Parameters
+      read(inp_unit,*) ip
+      if (ip >= 5) then
+        read(inp_unit,*) mtrx_file
+        read(inp_unit,*) rhs_file
+        read(inp_unit,*) filefmt
+        read(inp_unit,*) kmethd
+        read(inp_unit,*) ptype
+        read(inp_unit,*) afmt
+        read(inp_unit,*) ipart
+
+
+        call psb_bcast(ctxt,mtrx_file)
+        call psb_bcast(ctxt,rhs_file)
+        call psb_bcast(ctxt,filefmt)
+        call psb_bcast(ctxt,kmethd)
+        call psb_bcast(ctxt,ptype)
+        call psb_bcast(ctxt,afmt)
+        call psb_bcast(ctxt,part)
+
+        if (ip >= 7) then
+          read(inp_unit,*) istopc
+        else
+          istopc=1        
+        endif
+        if (ip >= 8) then
+          read(inp_unit,*) itmax
+        else
+          itmax=500
+        endif
+        if (ip >= 9) then
+          read(inp_unit,*) itrace
+        else
+          itrace=-1
+        endif
+        if (ip >= 10) then
+          read(inp_unit,*) irst
+        else
+          irst  = 1
+        endif
+        if (ip >= 11) then
+          read(inp_unit,*) eps
+        else
+          eps=1.d-6
+        endif
+        inparms(1) = istopc
+        inparms(2) = itmax
+        inparms(3) = itrace
+        inparms(4) = irst
+        call psb_bcast(ctxt,inparms(1:4))
+        call psb_bcast(ctxt,eps)
+
+        write(psb_out_unit,'("Solving matrix       : ",a)')  mtrx_file      
+        write(psb_out_unit,'("Number of processors : ",i3)') np
+        write(psb_out_unit,'("Data distribution    : ",a)') part
+        write(psb_out_unit,'("Iterative method     : ",a)')  kmethd
+        write(psb_out_unit,'("Preconditioner       : ",a)')  ptype
+        write(psb_out_unit,'("Restart parameter    : ",i2)') irst
+        write(psb_out_unit,'("Storage format       : ",a)')  afmt(1:3)
+        write(psb_out_unit,'(" ")')
+      else
+        write(psb_err_unit,*) 'Wrong format for input file'
+        call psb_abort(ctxt)
+        stop 1
+      end if
+      if (inp_unit /= psb_inp_unit) then
+        close(inp_unit)
+      end if
+    else
+      ! Receive Parameters
+      call psb_bcast(ctxt,mtrx_file)
+      call psb_bcast(ctxt,rhs_file)
+      call psb_bcast(ctxt,filefmt)
+      call psb_bcast(ctxt,kmethd)
+      call psb_bcast(ctxt,ptype)
+      call psb_bcast(ctxt,afmt)
+      call psb_bcast(ctxt,part)
+
+      call psb_bcast(ctxt,inparms(1:4))
+      istopc =  inparms(1) 
+      itmax  =  inparms(2) 
+      itrace =  inparms(3) 
+      irst   =  inparms(4) 
+      call psb_bcast(ctxt,eps)
+
+    end if
+
+  end subroutine get_sparms
+
+end module getp
diff --git a/test/krylovspace/psb_cf_sample.f90 b/test/krylovspace/psb_cf_sample.f90
new file mode 100644
index 00000000..bf37e0d8
--- /dev/null
+++ b/test/krylovspace/psb_cf_sample.f90
@@ -0,0 +1,250 @@
+!
+!                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.
+!
+!
+program psb_ckrylov_sample
+  use psb_base_mod
+  use psb_prec_mod
+  use psb_krylov_mod
+  use psb_util_mod
+  use psb_ckrylovsubspace_mod
+  use getp
+  implicit none
+
+  ! input parameters
+  character(len=40) :: kmethd, ptype, mtrx_file, rhs_file
+
+  ! sparse matrices
+  type(psb_cspmat_type)  :: a
+  type(psb_lcspmat_type) :: aux_a
+  ! Krylov space
+  type(psb_cpolkrylov) :: kryl
+
+  ! dense matrices
+  complex(psb_spk_), allocatable, target ::  aux_b(:,:), d(:)
+  complex(psb_spk_), allocatable , save  :: x_col_glob(:), r_col_glob(:)
+  complex(psb_spk_), pointer  :: b_col_glob(:)
+  type(psb_c_vect_type)    :: b_col, x_col, r_col
+
+  ! communications data structure
+  type(psb_desc_type):: desc_a
+
+  type(psb_ctxt_type) :: ctxt
+  integer(psb_ipk_) :: iam, np
+  integer(psb_lpk_) :: lnp
+  ! solver paramters
+  integer(psb_ipk_) :: iter, itmax, ierr, itrace, ircode,&
+       & methd, istopc, irst
+  integer(psb_epk_) :: amatsize, descsize
+  real(psb_spk_)   :: err, eps, cond
+
+  character(len=5)   :: afmt
+  character(len=20)  :: name, part
+  character(len=2)   :: filefmt
+  integer(psb_ipk_), parameter :: iunit=12
+  integer(psb_ipk_) :: iparm(20)
+
+  ! other variables
+  integer(psb_ipk_) :: i,info,j,m_problem, err_act
+  integer(psb_ipk_) :: internal, m,ii,nnzero
+  real(psb_dpk_) :: t1, t2, tprec
+  real(psb_spk_) :: r_amax, b_amax, scale,resmx,resmxp
+  integer(psb_ipk_) :: nrhs, nrow, n_row, dim, ne, nv
+  integer(psb_ipk_), allocatable :: ivg(:)
+  integer(psb_ipk_), allocatable :: ipv(:)
+  character(len=40)  :: fname, fnout
+
+
+  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
+
+
+  name='psb_ckrylov_sample'
+  if(psb_errstatus_fatal()) goto 9999
+  info=psb_success_
+  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,mtrx_file,rhs_file,filefmt,kmethd,ptype,&
+       & part,afmt,istopc,itmax,itrace,irst,eps)
+
+  call psb_barrier(ctxt)
+  t1 = psb_wtime()
+  ! read the input matrix to be processed and (possibly) the rhs
+  nrhs = 1
+
+  if (iam == psb_root_) then
+    select case(psb_toupper(filefmt))
+    case('MM')
+      ! For Matrix Market we have an input file for the matrix
+      ! and an (optional) second file for the RHS.
+      call mm_mat_read(aux_a,info,iunit=iunit,filename=mtrx_file)
+      if (info == psb_success_) then
+        if (rhs_file /= 'NONE') then
+          call mm_array_read(aux_b,info,iunit=iunit,filename=rhs_file)
+        end if
+      end if
+
+    case ('HB')
+      ! For Harwell-Boeing we have a single file which may or may not
+      ! contain an RHS.
+      call hb_read(aux_a,info,iunit=iunit,b=aux_b,filename=mtrx_file)
+
+    case default
+      info = -1
+      write(psb_err_unit,*) 'Wrong choice for fileformat ', filefmt
+    end select
+    if (info /= psb_success_) then
+      write(psb_err_unit,*) 'Error while reading input matrix '
+      call psb_abort(ctxt)
+    end if
+
+    m_problem = aux_a%get_nrows()
+    call psb_bcast(ctxt,m_problem)
+
+    ! At this point aux_b may still be unallocated
+    if (size(aux_b,dim=1) == m_problem) then
+      ! if any rhs were present, broadcast the first one
+      write(psb_err_unit,'("Ok, got an rhs ")')
+      b_col_glob =>aux_b(:,1)
+    else
+      write(psb_out_unit,'("Generating an rhs...")')
+      write(psb_out_unit,'(" ")')
+      call psb_realloc(m_problem,1,aux_b,ircode)
+      if (ircode /= 0) then
+         call psb_errpush(psb_err_alloc_dealloc_,name)
+         goto 9999
+      endif
+
+      b_col_glob => aux_b(:,1)
+      do i=1, m_problem
+         b_col_glob(i) = cone
+      enddo
+    endif
+
+  else
+    call psb_bcast(ctxt,m_problem)
+
+  end if
+
+  ! switch over different partition types
+  select case(psb_toupper(part))
+  case('BLOCK')
+    if (iam == psb_root_) write(psb_out_unit,'("Partition type: block")')
+    call psb_matdist(aux_a, a,  ctxt,desc_a,info,fmt=afmt,parts=part_block)
+
+  case('GRAPH')
+    if (iam == psb_root_) then
+      write(psb_out_unit,'("Partition type: graph vector")')
+      write(psb_out_unit,'(" ")')
+      !      write(psb_err_unit,'("Build type: graph")')
+      call aux_a%cscnv(info,type='csr')
+      lnp = np
+      call build_mtpart(aux_a,lnp)
+
+    endif
+    call psb_barrier(ctxt)
+    call distr_mtpart(psb_root_,ctxt)
+    call getv_mtpart(ivg)
+    call psb_matdist(aux_a, a, ctxt,desc_a,info,fmt=afmt,vg=ivg)
+
+  case default
+    if (iam == psb_root_) write(psb_out_unit,'("Partition type: block")')
+    call psb_matdist(aux_a, a,  ctxt,desc_a,info,fmt=afmt,parts=part_block)
+  end select
+
+  call psb_scatter(b_col_glob,b_col,desc_a,info,root=psb_root_)
+  call psb_geall(x_col,desc_a,info)
+  call x_col%zero()
+  call psb_geasb(x_col,desc_a,info)
+  call psb_geall(r_col,desc_a,info)
+  call r_col%zero()
+  call psb_geasb(r_col,desc_a,info)
+  t2 = psb_wtime() - t1
+
+
+  call psb_amx(ctxt, t2)
+
+  if (iam == psb_root_) then
+     write(psb_out_unit,'(" ")')
+     write(psb_out_unit,'("Time to read and partition matrix : ",es12.5)')t2
+     write(psb_out_unit,'(" ")')
+  end if
+
+  ! Allocate the Krylov space
+  call kryl%allocate(desc_a,itmax,info)
+  ! Build
+  call kryl%arnoldi(a,desc_a,b_col,istopc)
+  ! Write
+  call kryl%write(desc_a,mtrx_file)
+
+  amatsize = a%sizeof()
+  descsize = desc_a%sizeof()
+  call psb_sum(ctxt,amatsize)
+  call psb_sum(ctxt,descsize)
+  if (iam == psb_root_) then
+    write(psb_out_unit,'("Matrix: ",a)')mtrx_file
+    write(psb_out_unit,'("Computed space on ",i8," processors")')np
+    write(psb_out_unit,'("Time to build space      : ",es12.5)')t2
+    write(psb_out_unit,'("Time per iteration       : ",es12.5)')t2/(iter)
+    write(psb_out_unit,'("Total memory occupation for A:      ",i12)')amatsize
+    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
+
+  call psb_gefree(b_col, desc_a,info)
+  call psb_gefree(x_col, desc_a,info)
+  call psb_spfree(a, desc_a,info)
+  call kryl%free(desc_a,info)
+  call psb_cdfree(desc_a,info)
+  call psb_exit(ctxt)
+  stop
+
+9999 call psb_error(ctxt)
+
+  stop
+end program psb_ckrylov_sample
diff --git a/test/krylovspace/psb_df_sample.f90 b/test/krylovspace/psb_df_sample.f90
new file mode 100644
index 00000000..1c1dad91
--- /dev/null
+++ b/test/krylovspace/psb_df_sample.f90
@@ -0,0 +1,250 @@
+!
+!                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.
+!
+!
+program psb_dkrylov_sample
+  use psb_base_mod
+  use psb_prec_mod
+  use psb_krylov_mod
+  use psb_util_mod
+  use psb_dkrylovsubspace_mod
+  use getp
+  implicit none
+
+  ! input parameters
+  character(len=40) :: kmethd, ptype, mtrx_file, rhs_file
+
+  ! sparse matrices
+  type(psb_dspmat_type)  :: a
+  type(psb_ldspmat_type) :: aux_a
+  ! Krylov space
+  type(psb_dpolkrylov) :: kryl
+
+  ! dense matrices
+  real(psb_dpk_), allocatable, target ::  aux_b(:,:), d(:)
+  real(psb_dpk_), allocatable , save  :: x_col_glob(:), r_col_glob(:)
+  real(psb_dpk_), pointer  :: b_col_glob(:)
+  type(psb_d_vect_type)    :: b_col, x_col, r_col
+
+  ! communications data structure
+  type(psb_desc_type):: desc_a
+
+  type(psb_ctxt_type) :: ctxt
+  integer(psb_ipk_) :: iam, np
+  integer(psb_lpk_) :: lnp
+  ! solver paramters
+  integer(psb_ipk_) :: iter, itmax, ierr, itrace, ircode,&
+       & methd, istopc, irst
+  integer(psb_epk_) :: amatsize, descsize
+  real(psb_dpk_)   :: err, eps, cond
+
+  character(len=5)   :: afmt
+  character(len=20)  :: name, part
+  character(len=2)   :: filefmt
+  integer(psb_ipk_), parameter :: iunit=12
+  integer(psb_ipk_) :: iparm(20)
+
+  ! other variables
+  integer(psb_ipk_) :: i,info,j,m_problem, err_act
+  integer(psb_ipk_) :: internal, m,ii,nnzero
+  real(psb_dpk_) :: t1, t2, tprec
+  real(psb_dpk_) :: r_amax, b_amax, scale,resmx,resmxp
+  integer(psb_ipk_) :: nrhs, nrow, n_row, dim, ne, nv
+  integer(psb_ipk_), allocatable :: ivg(:)
+  integer(psb_ipk_), allocatable :: ipv(:)
+  character(len=40)  :: fname, fnout
+
+
+  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
+
+
+  name='psb_dkrylov_sample'
+  if(psb_errstatus_fatal()) goto 9999
+  info=psb_success_
+  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,mtrx_file,rhs_file,filefmt,kmethd,ptype,&
+       & part,afmt,istopc,itmax,itrace,irst,eps)
+
+  call psb_barrier(ctxt)
+  t1 = psb_wtime()
+  ! read the input matrix to be processed and (possibly) the rhs
+  nrhs = 1
+
+  if (iam == psb_root_) then
+    select case(psb_toupper(filefmt))
+    case('MM')
+      ! For Matrix Market we have an input file for the matrix
+      ! and an (optional) second file for the RHS.
+      call mm_mat_read(aux_a,info,iunit=iunit,filename=mtrx_file)
+      if (info == psb_success_) then
+        if (rhs_file /= 'NONE') then
+          call mm_array_read(aux_b,info,iunit=iunit,filename=rhs_file)
+        end if
+      end if
+
+    case ('HB')
+      ! For Harwell-Boeing we have a single file which may or may not
+      ! contain an RHS.
+      call hb_read(aux_a,info,iunit=iunit,b=aux_b,filename=mtrx_file)
+
+    case default
+      info = -1
+      write(psb_err_unit,*) 'Wrong choice for fileformat ', filefmt
+    end select
+    if (info /= psb_success_) then
+      write(psb_err_unit,*) 'Error while reading input matrix '
+      call psb_abort(ctxt)
+    end if
+
+    m_problem = aux_a%get_nrows()
+    call psb_bcast(ctxt,m_problem)
+
+    ! At this point aux_b may still be unallocated
+    if (size(aux_b,dim=1) == m_problem) then
+      ! if any rhs were present, broadcast the first one
+      write(psb_err_unit,'("Ok, got an rhs ")')
+      b_col_glob =>aux_b(:,1)
+    else
+      write(psb_out_unit,'("Generating an rhs...")')
+      write(psb_out_unit,'(" ")')
+      call psb_realloc(m_problem,1,aux_b,ircode)
+      if (ircode /= 0) then
+         call psb_errpush(psb_err_alloc_dealloc_,name)
+         goto 9999
+      endif
+
+      b_col_glob => aux_b(:,1)
+      do i=1, m_problem
+         b_col_glob(i) = done
+      enddo
+    endif
+
+  else
+    call psb_bcast(ctxt,m_problem)
+
+  end if
+
+  ! switch over different partition types
+  select case(psb_toupper(part))
+  case('BLOCK')
+    if (iam == psb_root_) write(psb_out_unit,'("Partition type: block")')
+    call psb_matdist(aux_a, a,  ctxt,desc_a,info,fmt=afmt,parts=part_block)
+
+  case('GRAPH')
+    if (iam == psb_root_) then
+      write(psb_out_unit,'("Partition type: graph vector")')
+      write(psb_out_unit,'(" ")')
+      !      write(psb_err_unit,'("Build type: graph")')
+      call aux_a%cscnv(info,type='csr')
+      lnp = np
+      call build_mtpart(aux_a,lnp)
+
+    endif
+    call psb_barrier(ctxt)
+    call distr_mtpart(psb_root_,ctxt)
+    call getv_mtpart(ivg)
+    call psb_matdist(aux_a, a, ctxt,desc_a,info,fmt=afmt,vg=ivg)
+
+  case default
+    if (iam == psb_root_) write(psb_out_unit,'("Partition type: block")')
+    call psb_matdist(aux_a, a,  ctxt,desc_a,info,fmt=afmt,parts=part_block)
+  end select
+
+  call psb_scatter(b_col_glob,b_col,desc_a,info,root=psb_root_)
+  call psb_geall(x_col,desc_a,info)
+  call x_col%zero()
+  call psb_geasb(x_col,desc_a,info)
+  call psb_geall(r_col,desc_a,info)
+  call r_col%zero()
+  call psb_geasb(r_col,desc_a,info)
+  t2 = psb_wtime() - t1
+
+
+  call psb_amx(ctxt, t2)
+
+  if (iam == psb_root_) then
+     write(psb_out_unit,'(" ")')
+     write(psb_out_unit,'("Time to read and partition matrix : ",es12.5)')t2
+     write(psb_out_unit,'(" ")')
+  end if
+
+  ! Allocate the Krylov space
+  call kryl%allocate(desc_a,itmax,info)
+  ! Build
+  call kryl%arnoldi(a,desc_a,b_col,istopc)
+  ! Write
+  call kryl%write(desc_a,mtrx_file)
+
+  amatsize = a%sizeof()
+  descsize = desc_a%sizeof()
+  call psb_sum(ctxt,amatsize)
+  call psb_sum(ctxt,descsize)
+  if (iam == psb_root_) then
+    write(psb_out_unit,'("Matrix: ",a)')mtrx_file
+    write(psb_out_unit,'("Computed space on ",i8," processors")')np
+    write(psb_out_unit,'("Time to build space      : ",es12.5)')t2
+    write(psb_out_unit,'("Time per iteration       : ",es12.5)')t2/(iter)
+    write(psb_out_unit,'("Total memory occupation for A:      ",i12)')amatsize
+    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
+
+  call psb_gefree(b_col, desc_a,info)
+  call psb_gefree(x_col, desc_a,info)
+  call psb_spfree(a, desc_a,info)
+  call kryl%free(desc_a,info)
+  call psb_cdfree(desc_a,info)
+  call psb_exit(ctxt)
+  stop
+
+9999 call psb_error(ctxt)
+
+  stop
+end program psb_dkrylov_sample
diff --git a/test/krylovspace/psb_sf_sample.f90 b/test/krylovspace/psb_sf_sample.f90
new file mode 100644
index 00000000..63ea31da
--- /dev/null
+++ b/test/krylovspace/psb_sf_sample.f90
@@ -0,0 +1,250 @@
+!
+!                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.
+!
+!
+program psb_skrylov_sample
+  use psb_base_mod
+  use psb_prec_mod
+  use psb_krylov_mod
+  use psb_util_mod
+  use psb_skrylovsubspace_mod
+  use getp
+  implicit none
+
+  ! input parameters
+  character(len=40) :: kmethd, ptype, mtrx_file, rhs_file
+
+  ! sparse matrices
+  type(psb_sspmat_type)  :: a
+  type(psb_lsspmat_type) :: aux_a
+  ! Krylov space
+  type(psb_spolkrylov) :: kryl
+
+  ! dense matrices
+  real(psb_spk_), allocatable, target ::  aux_b(:,:), d(:)
+  real(psb_spk_), allocatable , save  :: x_col_glob(:), r_col_glob(:)
+  real(psb_spk_), pointer  :: b_col_glob(:)
+  type(psb_s_vect_type)    :: b_col, x_col, r_col
+
+  ! communications data structure
+  type(psb_desc_type):: desc_a
+
+  type(psb_ctxt_type) :: ctxt
+  integer(psb_ipk_) :: iam, np
+  integer(psb_lpk_) :: lnp
+  ! solver paramters
+  integer(psb_ipk_) :: iter, itmax, ierr, itrace, ircode,&
+       & methd, istopc, irst
+  integer(psb_epk_) :: amatsize, descsize
+  real(psb_spk_)   :: err, eps, cond
+
+  character(len=5)   :: afmt
+  character(len=20)  :: name, part
+  character(len=2)   :: filefmt
+  integer(psb_ipk_), parameter :: iunit=12
+  integer(psb_ipk_) :: iparm(20)
+
+  ! other variables
+  integer(psb_ipk_) :: i,info,j,m_problem, err_act
+  integer(psb_ipk_) :: internal, m,ii,nnzero
+  real(psb_dpk_) :: t1, t2, tprec
+  real(psb_spk_) :: r_amax, b_amax, scale,resmx,resmxp
+  integer(psb_ipk_) :: nrhs, nrow, n_row, dim, ne, nv
+  integer(psb_ipk_), allocatable :: ivg(:)
+  integer(psb_ipk_), allocatable :: ipv(:)
+  character(len=40)  :: fname, fnout
+
+
+  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
+
+
+  name='psb_skrylov_sample'
+  if(psb_errstatus_fatal()) goto 9999
+  info=psb_success_
+  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,mtrx_file,rhs_file,filefmt,kmethd,ptype,&
+       & part,afmt,istopc,itmax,itrace,irst,eps)
+
+  call psb_barrier(ctxt)
+  t1 = psb_wtime()
+  ! read the input matrix to be processed and (possibly) the rhs
+  nrhs = 1
+
+  if (iam == psb_root_) then
+    select case(psb_toupper(filefmt))
+    case('MM')
+      ! For Matrix Market we have an input file for the matrix
+      ! and an (optional) second file for the RHS.
+      call mm_mat_read(aux_a,info,iunit=iunit,filename=mtrx_file)
+      if (info == psb_success_) then
+        if (rhs_file /= 'NONE') then
+          call mm_array_read(aux_b,info,iunit=iunit,filename=rhs_file)
+        end if
+      end if
+
+    case ('HB')
+      ! For Harwell-Boeing we have a single file which may or may not
+      ! contain an RHS.
+      call hb_read(aux_a,info,iunit=iunit,b=aux_b,filename=mtrx_file)
+
+    case default
+      info = -1
+      write(psb_err_unit,*) 'Wrong choice for fileformat ', filefmt
+    end select
+    if (info /= psb_success_) then
+      write(psb_err_unit,*) 'Error while reading input matrix '
+      call psb_abort(ctxt)
+    end if
+
+    m_problem = aux_a%get_nrows()
+    call psb_bcast(ctxt,m_problem)
+
+    ! At this point aux_b may still be unallocated
+    if (size(aux_b,dim=1) == m_problem) then
+      ! if any rhs were present, broadcast the first one
+      write(psb_err_unit,'("Ok, got an rhs ")')
+      b_col_glob =>aux_b(:,1)
+    else
+      write(psb_out_unit,'("Generating an rhs...")')
+      write(psb_out_unit,'(" ")')
+      call psb_realloc(m_problem,1,aux_b,ircode)
+      if (ircode /= 0) then
+         call psb_errpush(psb_err_alloc_dealloc_,name)
+         goto 9999
+      endif
+
+      b_col_glob => aux_b(:,1)
+      do i=1, m_problem
+         b_col_glob(i) = sone
+      enddo
+    endif
+
+  else
+    call psb_bcast(ctxt,m_problem)
+
+  end if
+
+  ! switch over different partition types
+  select case(psb_toupper(part))
+  case('BLOCK')
+    if (iam == psb_root_) write(psb_out_unit,'("Partition type: block")')
+    call psb_matdist(aux_a, a,  ctxt,desc_a,info,fmt=afmt,parts=part_block)
+
+  case('GRAPH')
+    if (iam == psb_root_) then
+      write(psb_out_unit,'("Partition type: graph vector")')
+      write(psb_out_unit,'(" ")')
+      !      write(psb_err_unit,'("Build type: graph")')
+      call aux_a%cscnv(info,type='csr')
+      lnp = np
+      call build_mtpart(aux_a,lnp)
+
+    endif
+    call psb_barrier(ctxt)
+    call distr_mtpart(psb_root_,ctxt)
+    call getv_mtpart(ivg)
+    call psb_matdist(aux_a, a, ctxt,desc_a,info,fmt=afmt,vg=ivg)
+
+  case default
+    if (iam == psb_root_) write(psb_out_unit,'("Partition type: block")')
+    call psb_matdist(aux_a, a,  ctxt,desc_a,info,fmt=afmt,parts=part_block)
+  end select
+
+  call psb_scatter(b_col_glob,b_col,desc_a,info,root=psb_root_)
+  call psb_geall(x_col,desc_a,info)
+  call x_col%zero()
+  call psb_geasb(x_col,desc_a,info)
+  call psb_geall(r_col,desc_a,info)
+  call r_col%zero()
+  call psb_geasb(r_col,desc_a,info)
+  t2 = psb_wtime() - t1
+
+
+  call psb_amx(ctxt, t2)
+
+  if (iam == psb_root_) then
+     write(psb_out_unit,'(" ")')
+     write(psb_out_unit,'("Time to read and partition matrix : ",es12.5)')t2
+     write(psb_out_unit,'(" ")')
+  end if
+
+  ! Allocate the Krylov space
+  call kryl%allocate(desc_a,itmax,info)
+  ! Build
+  call kryl%arnoldi(a,desc_a,b_col,istopc)
+  ! Write
+  call kryl%write(desc_a,mtrx_file)
+
+  amatsize = a%sizeof()
+  descsize = desc_a%sizeof()
+  call psb_sum(ctxt,amatsize)
+  call psb_sum(ctxt,descsize)
+  if (iam == psb_root_) then
+    write(psb_out_unit,'("Matrix: ",a)')mtrx_file
+    write(psb_out_unit,'("Computed space on ",i8," processors")')np
+    write(psb_out_unit,'("Time to build space      : ",es12.5)')t2
+    write(psb_out_unit,'("Time per iteration       : ",es12.5)')t2/(iter)
+    write(psb_out_unit,'("Total memory occupation for A:      ",i12)')amatsize
+    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
+
+  call psb_gefree(b_col, desc_a,info)
+  call psb_gefree(x_col, desc_a,info)
+  call psb_spfree(a, desc_a,info)
+  call kryl%free(desc_a,info)
+  call psb_cdfree(desc_a,info)
+  call psb_exit(ctxt)
+  stop
+
+9999 call psb_error(ctxt)
+
+  stop
+end program psb_skrylov_sample
diff --git a/test/krylovspace/psb_zf_sample.f90 b/test/krylovspace/psb_zf_sample.f90
new file mode 100644
index 00000000..09d61c26
--- /dev/null
+++ b/test/krylovspace/psb_zf_sample.f90
@@ -0,0 +1,250 @@
+!
+!                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.
+!
+!
+program psb_zkrylov_sample
+  use psb_base_mod
+  use psb_prec_mod
+  use psb_krylov_mod
+  use psb_util_mod
+  use psb_zkrylovsubspace_mod
+  use getp
+  implicit none
+
+  ! input parameters
+  character(len=40) :: kmethd, ptype, mtrx_file, rhs_file
+
+  ! sparse matrices
+  type(psb_zspmat_type)  :: a
+  type(psb_lzspmat_type) :: aux_a
+  ! Krylov space
+  type(psb_zpolkrylov) :: kryl
+
+  ! dense matrices
+  complex(psb_dpk_), allocatable, target ::  aux_b(:,:), d(:)
+  complex(psb_dpk_), allocatable , save  :: x_col_glob(:), r_col_glob(:)
+  complex(psb_dpk_), pointer  :: b_col_glob(:)
+  type(psb_z_vect_type)    :: b_col, x_col, r_col
+
+  ! communications data structure
+  type(psb_desc_type):: desc_a
+
+  type(psb_ctxt_type) :: ctxt
+  integer(psb_ipk_) :: iam, np
+  integer(psb_lpk_) :: lnp
+  ! solver paramters
+  integer(psb_ipk_) :: iter, itmax, ierr, itrace, ircode,&
+       & methd, istopc, irst
+  integer(psb_epk_) :: amatsize, descsize
+  real(psb_dpk_)   :: err, eps, cond
+
+  character(len=5)   :: afmt
+  character(len=20)  :: name, part
+  character(len=2)   :: filefmt
+  integer(psb_ipk_), parameter :: iunit=12
+  integer(psb_ipk_) :: iparm(20)
+
+  ! other variables
+  integer(psb_ipk_) :: i,info,j,m_problem, err_act
+  integer(psb_ipk_) :: internal, m,ii,nnzero
+  real(psb_dpk_) :: t1, t2, tprec
+  real(psb_dpk_) :: r_amax, b_amax, scale,resmx,resmxp
+  integer(psb_ipk_) :: nrhs, nrow, n_row, dim, ne, nv
+  integer(psb_ipk_), allocatable :: ivg(:)
+  integer(psb_ipk_), allocatable :: ipv(:)
+  character(len=40)  :: fname, fnout
+
+
+  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
+
+
+  name='psb_zkrylov_sample'
+  if(psb_errstatus_fatal()) goto 9999
+  info=psb_success_
+  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,mtrx_file,rhs_file,filefmt,kmethd,ptype,&
+       & part,afmt,istopc,itmax,itrace,irst,eps)
+
+  call psb_barrier(ctxt)
+  t1 = psb_wtime()
+  ! read the input matrix to be processed and (possibly) the rhs
+  nrhs = 1
+
+  if (iam == psb_root_) then
+    select case(psb_toupper(filefmt))
+    case('MM')
+      ! For Matrix Market we have an input file for the matrix
+      ! and an (optional) second file for the RHS.
+      call mm_mat_read(aux_a,info,iunit=iunit,filename=mtrx_file)
+      if (info == psb_success_) then
+        if (rhs_file /= 'NONE') then
+          call mm_array_read(aux_b,info,iunit=iunit,filename=rhs_file)
+        end if
+      end if
+
+    case ('HB')
+      ! For Harwell-Boeing we have a single file which may or may not
+      ! contain an RHS.
+      call hb_read(aux_a,info,iunit=iunit,b=aux_b,filename=mtrx_file)
+
+    case default
+      info = -1
+      write(psb_err_unit,*) 'Wrong choice for fileformat ', filefmt
+    end select
+    if (info /= psb_success_) then
+      write(psb_err_unit,*) 'Error while reading input matrix '
+      call psb_abort(ctxt)
+    end if
+
+    m_problem = aux_a%get_nrows()
+    call psb_bcast(ctxt,m_problem)
+
+    ! At this point aux_b may still be unallocated
+    if (size(aux_b,dim=1) == m_problem) then
+      ! if any rhs were present, broadcast the first one
+      write(psb_err_unit,'("Ok, got an rhs ")')
+      b_col_glob =>aux_b(:,1)
+    else
+      write(psb_out_unit,'("Generating an rhs...")')
+      write(psb_out_unit,'(" ")')
+      call psb_realloc(m_problem,1,aux_b,ircode)
+      if (ircode /= 0) then
+         call psb_errpush(psb_err_alloc_dealloc_,name)
+         goto 9999
+      endif
+
+      b_col_glob => aux_b(:,1)
+      do i=1, m_problem
+         b_col_glob(i) = zone
+      enddo
+    endif
+
+  else
+    call psb_bcast(ctxt,m_problem)
+
+  end if
+
+  ! switch over different partition types
+  select case(psb_toupper(part))
+  case('BLOCK')
+    if (iam == psb_root_) write(psb_out_unit,'("Partition type: block")')
+    call psb_matdist(aux_a, a,  ctxt,desc_a,info,fmt=afmt,parts=part_block)
+
+  case('GRAPH')
+    if (iam == psb_root_) then
+      write(psb_out_unit,'("Partition type: graph vector")')
+      write(psb_out_unit,'(" ")')
+      !      write(psb_err_unit,'("Build type: graph")')
+      call aux_a%cscnv(info,type='csr')
+      lnp = np
+      call build_mtpart(aux_a,lnp)
+
+    endif
+    call psb_barrier(ctxt)
+    call distr_mtpart(psb_root_,ctxt)
+    call getv_mtpart(ivg)
+    call psb_matdist(aux_a, a, ctxt,desc_a,info,fmt=afmt,vg=ivg)
+
+  case default
+    if (iam == psb_root_) write(psb_out_unit,'("Partition type: block")')
+    call psb_matdist(aux_a, a,  ctxt,desc_a,info,fmt=afmt,parts=part_block)
+  end select
+
+  call psb_scatter(b_col_glob,b_col,desc_a,info,root=psb_root_)
+  call psb_geall(x_col,desc_a,info)
+  call x_col%zero()
+  call psb_geasb(x_col,desc_a,info)
+  call psb_geall(r_col,desc_a,info)
+  call r_col%zero()
+  call psb_geasb(r_col,desc_a,info)
+  t2 = psb_wtime() - t1
+
+
+  call psb_amx(ctxt, t2)
+
+  if (iam == psb_root_) then
+     write(psb_out_unit,'(" ")')
+     write(psb_out_unit,'("Time to read and partition matrix : ",es12.5)')t2
+     write(psb_out_unit,'(" ")')
+  end if
+
+  ! Allocate the Krylov space
+  call kryl%allocate(desc_a,itmax,info)
+  ! Build
+  call kryl%arnoldi(a,desc_a,b_col,istopc)
+  ! Write
+  call kryl%write(desc_a,mtrx_file)
+
+  amatsize = a%sizeof()
+  descsize = desc_a%sizeof()
+  call psb_sum(ctxt,amatsize)
+  call psb_sum(ctxt,descsize)
+  if (iam == psb_root_) then
+    write(psb_out_unit,'("Matrix: ",a)')mtrx_file
+    write(psb_out_unit,'("Computed space on ",i8," processors")')np
+    write(psb_out_unit,'("Time to build space      : ",es12.5)')t2
+    write(psb_out_unit,'("Time per iteration       : ",es12.5)')t2/(iter)
+    write(psb_out_unit,'("Total memory occupation for A:      ",i12)')amatsize
+    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
+
+  call psb_gefree(b_col, desc_a,info)
+  call psb_gefree(x_col, desc_a,info)
+  call psb_spfree(a, desc_a,info)
+  call kryl%free(desc_a,info)
+  call psb_cdfree(desc_a,info)
+  call psb_exit(ctxt)
+  stop
+
+9999 call psb_error(ctxt)
+
+  stop
+end program psb_zkrylov_sample