psblas-dev:

17 years ago · c1d96960d8
parent ffb96464ee
commit c1d96960d8
241 changed files with 5262 additions and 4977 deletions
--- a/28
+++ b/28
@ -69,7 +69,8 @@ prerequisites (see also SERIAL below):
 1. A working version of MPI
-2. The MPI version of the BLACS  from
+2. The MPI version of the BLACS; some systems (e.g. IBM SP) provide
   a specific version, otherwise you can find it at
   http://www.netlib.org/blacs/index.html 
 3. A version of the BLAS; if you don't have a specific version for
@ -78,14 +79,14 @@ prerequisites (see also SERIAL below):
 4. We have had good results with  the METIS library, from 
   http://www-users.cs.umn.edu/~karypis/metis/metis/main.html 
-   This is optional: it is  used in the util and test/fileread
+   This is optional; it is  used in the util and test/fileread
-   directories if you define the HAVE_METIS directive.
+   directories but only if you define the HAVE_METIS directive.
 We offer an experimental configure script: if everything works well,
 it will generate a Make.inc file suitable for building the library. 
 The script is capable of recognizing the needed libraries with their
 default names; if they are in unusual places consider adding the paths
-with --with-lib, or explicitly specify the names in --with-blas and
+with --with-lib, or explicitly specifying the names in --with-blas and
 friends.
 Please note that a common way for the configure script to fail is to
 specify inconsistent MPI vs. plain compilers, either directly or
@ -105,16 +106,18 @@ An annoying problem exists with some versions of MPICH: the configure
 script will set -DMPI_MOD, which is to say, the MPI call interfaces
 will be resolved by using the MPI Fortran module. However usage of the
 module may cause compilation to fail if coupled with the debugging
-option -g, because the compiler complains that it cannot find any
+option -ggdb, because the compiler complains that it cannot find any
 matching interface. 
-The solution: either take out the -g option, or, if you really need to
+The solution: either take out the -ggdb option, or, if you really need
-debug, force -DMPI_H in place of -DMPI_MOD.
+to debug, force -DMPI_H in place of -DMPI_MOD.
-As a backup alternative, you can  choose a Make.inc.XXX file in
+As a backup alternative, you can always choose a Make.inc.XXX file in
 directory Make that fits your compilers, modify the paths to libraries
 to  match your installation and copy it to Make.inc in the top
-directory. Then run 
+directory. 
 After you have Make.inc fixed,  run 
 make 
 to  compile the library; go to the test directory and its
 subdirectories to get test programs done.   
@ -131,6 +134,8 @@ SERIAL: We now provide an (experimental) option to run in serial
 TODO:
 Fix all reamining bugs. Bugs? We dont' have any ! ;-)
 Work on a single precision version if there is sufficient demand for
 it.
 The PSBLAS team. 
@ -145,8 +150,9 @@ Salvatore  Filippone
 Alfredo    Buttari
-In a sister package we are going to provide a library of  multilevel
+In a sister package called MLD2P4 we are going to provide a library of
-parallel preconditioners that were developed with the contribution of:  
+more sophisticated preconditioners, including parallel multilevel ones
 that were  developed with the contribution of:  
 Pasqua     D'Ambra
 Daniela    di Serafino
--- a/base/comm/psb_dgather.f90
+++ b/base/comm/psb_dgather.f90
@ -53,8 +53,8 @@ subroutine  psb_dgatherm(globx, locx, desc_a, info, iroot)
  use psb_penv_mod
  implicit none
-  real(kind(1.d0)), intent(in)    :: locx(:,:)
+  real(psb_dpk_), intent(in)    :: locx(:,:)
-  real(kind(1.d0)), intent(out)   :: globx(:,:)
+  real(psb_dpk_), intent(out)   :: globx(:,:)
  type(psb_desc_type), intent(in) :: desc_a
  integer, intent(out)            :: info
  integer, intent(in), optional   :: iroot
@ -225,8 +225,8 @@ subroutine  psb_dgatherv(globx, locx, desc_a, info, iroot)
  use psb_penv_mod
  implicit none
-  real(kind(1.d0)), intent(in)    :: locx(:)
+  real(psb_dpk_), intent(in)    :: locx(:)
-  real(kind(1.d0)), intent(out)   :: globx(:)
+  real(psb_dpk_), intent(out)   :: globx(:)
  type(psb_desc_type), intent(in) :: desc_a
  integer, intent(out)            :: info
  integer, intent(in), optional   :: iroot
--- a/base/comm/psb_dhalo.f90
+++ b/base/comm/psb_dhalo.f90
@ -64,11 +64,11 @@ subroutine  psb_dhalom(x,desc_a,info,alpha,jx,ik,work,tran,mode,data)
  use psb_penv_mod
  implicit none
-  real(kind(1.d0)), intent(inout), target   :: x(:,:)
+  real(psb_dpk_), intent(inout), target   :: x(:,:)
  type(psb_desc_type), intent(in)           :: desc_a
  integer, intent(out)                      :: info
-  real(kind(1.d0)), intent(in), optional    :: alpha
+  real(psb_dpk_), intent(in), optional    :: alpha
-  real(kind(1.d0)), optional, target        :: work(:)
+  real(psb_dpk_), optional, target        :: work(:)
  integer, intent(in), optional             :: mode,jx,ik,data
  character, intent(in), optional           :: tran
@ -76,7 +76,7 @@ subroutine  psb_dhalom(x,desc_a,info,alpha,jx,ik,work,tran,mode,data)
  integer                  :: ictxt, np, me,&
       & err_act, m, n, iix, jjx, ix, ijx, k, maxk, nrow, imode, i,&
       & err, liwork,data_
-  real(kind(1.d0)),pointer :: iwork(:), xp(:,:)
+  real(psb_dpk_),pointer :: iwork(:), xp(:,:)
  character                :: tran_
  character(len=20)        :: name, ch_err
  logical                  :: aliw
@ -291,11 +291,11 @@ subroutine  psb_dhalov(x,desc_a,info,alpha,work,tran,mode,data)
  use psb_penv_mod
  implicit none
-  real(kind(1.d0)), intent(inout)           :: x(:)
+  real(psb_dpk_), intent(inout)           :: x(:)
  type(psb_desc_type), intent(in)           :: desc_a
  integer, intent(out)                      :: info
-  real(kind(1.d0)), intent(in), optional    :: alpha
+  real(psb_dpk_), intent(in), optional    :: alpha
-  real(kind(1.d0)), target, optional        :: work(:)
+  real(psb_dpk_), target, optional        :: work(:)
  integer, intent(in), optional             :: mode,data
  character, intent(in), optional           :: tran
@ -303,7 +303,7 @@ subroutine  psb_dhalov(x,desc_a,info,alpha,work,tran,mode,data)
  integer                  :: ictxt, np, me,&
       & err_act, m, n, iix, jjx, ix, ijx, nrow, imode,&
       & err, liwork,data_
-  real(kind(1.d0)),pointer :: iwork(:)
+  real(psb_dpk_),pointer :: iwork(:)
  character                :: tran_
  character(len=20)        :: name, ch_err
  logical                  :: aliw
--- a/base/comm/psb_dovrl.f90
+++ b/base/comm/psb_dovrl.f90
@ -73,17 +73,17 @@ subroutine  psb_dovrlm(x,desc_a,info,jx,ik,work,update,mode)
  use psb_penv_mod
  implicit none
-  real(kind(1.d0)), intent(inout), target   :: x(:,:)
+  real(psb_dpk_), intent(inout), target   :: x(:,:)
  type(psb_desc_type), intent(in)           :: desc_a
  integer, intent(out)                      :: info
-  real(kind(1.d0)), optional, target        :: work(:)
+  real(psb_dpk_), optional, target        :: work(:)
  integer, intent(in), optional             :: update,jx,ik,mode
  ! locals
  integer                  :: ictxt, np, me, &
       & err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, maxk, update_,&
       & mode_, err, liwork
-  real(kind(1.d0)),pointer :: iwork(:), xp(:,:)
+  real(psb_dpk_),pointer :: iwork(:), xp(:,:)
  logical                  :: do_swap
  character(len=20)        :: name, ch_err
  logical                  :: aliw
@ -280,17 +280,17 @@ subroutine  psb_dovrlv(x,desc_a,info,work,update,mode)
  use psb_penv_mod
  implicit none
-  real(kind(1.d0)), intent(inout), target   :: x(:)
+  real(psb_dpk_), intent(inout), target   :: x(:)
  type(psb_desc_type), intent(in)           :: desc_a
  integer, intent(out)                      :: info
-  real(kind(1.d0)), optional, target        :: work(:)
+  real(psb_dpk_), optional, target        :: work(:)
  integer, intent(in), optional             :: update,mode
  ! locals
  integer                  :: ictxt, np, me, &
       & err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, update_,&
       & mode_, err, liwork
-  real(kind(1.d0)),pointer :: iwork(:)
+  real(psb_dpk_),pointer :: iwork(:)
  logical                  :: do_swap
  character(len=20)        :: name, ch_err
  logical                  :: aliw
--- a/base/comm/psb_dscatter.F90
+++ b/base/comm/psb_dscatter.F90
@ -57,8 +57,8 @@ subroutine  psb_dscatterm(globx, locx, desc_a, info, iroot)
  include 'mpif.h'
 #endif
-  real(kind(1.d0)), intent(out)    :: locx(:,:)
+  real(psb_dpk_), intent(out)    :: locx(:,:)
-  real(kind(1.d0)), intent(in)     :: globx(:,:)
+  real(psb_dpk_), intent(in)     :: globx(:,:)
  type(psb_desc_type), intent(in)  :: desc_a
  integer, intent(out)             :: info
  integer, intent(in), optional    :: iroot
@ -69,7 +69,7 @@ subroutine  psb_dscatterm(globx, locx, desc_a, info, iroot)
       & err_act, m, n, i, j, idx, nrow, iiroot, iglobx, jglobx,&
       & ilocx, jlocx, lda_locx, lda_globx, lock, globk, icomm, k, maxk, root, ilx,&
       & jlx, myrank, rootrank, c, pos
-  real(kind(1.d0)), allocatable  :: scatterv(:)
+  real(psb_dpk_), allocatable  :: scatterv(:)
  integer, allocatable           :: displ(:), l_t_g_all(:), all_dim(:)
  character(len=20)        :: name, ch_err
@ -286,8 +286,8 @@ subroutine  psb_dscatterv(globx, locx, desc_a, info, iroot)
  include 'mpif.h'
 #endif
-  real(kind(1.d0)), intent(out)    :: locx(:)
+  real(psb_dpk_), intent(out)    :: locx(:)
-  real(kind(1.d0)), intent(in)     :: globx(:)
+  real(psb_dpk_), intent(in)     :: globx(:)
  type(psb_desc_type), intent(in)  :: desc_a
  integer, intent(out)             :: info
  integer, intent(in), optional    :: iroot
@ -298,7 +298,7 @@ subroutine  psb_dscatterv(globx, locx, desc_a, info, iroot)
       & err_act, m, n, i, j, idx, nrow, iglobx, jglobx,&
       & ilocx, jlocx, lda_locx, lda_globx, root, k, icomm, myrank,&
       & rootrank, pos, ilx, jlx
-  real(kind(1.d0)), allocatable  :: scatterv(:)
+  real(psb_dpk_), allocatable  :: scatterv(:)
  integer, allocatable           :: displ(:), l_t_g_all(:), all_dim(:)
  character(len=20)        :: name, ch_err
--- a/base/comm/psb_ihalo.f90
+++ b/base/comm/psb_ihalo.f90
@ -68,7 +68,7 @@ subroutine  psb_ihalom(x,desc_a,info,alpha,jx,ik,work,tran,mode,data)
  integer, intent(inout), target           :: x(:,:)
  type(psb_desc_type), intent(in)          :: desc_a
  integer, intent(out)                     :: info
-  real(kind(1.d0)), intent(in), optional   :: alpha
+  real(psb_dpk_), intent(in), optional   :: alpha
  integer, intent(inout), optional, target :: work(:)
  integer, intent(in), optional            :: mode,jx,ik,data
  character, intent(in), optional          :: tran
@ -298,7 +298,7 @@ subroutine  psb_ihalov(x,desc_a,info,alpha,work,tran,mode,data)
  integer, intent(inout)                   :: x(:)
  type(psb_desc_type), intent(in)          :: desc_a
  integer, intent(out)                     :: info
-  real(kind(1.d0)), intent(in), optional   :: alpha
+  real(psb_dpk_), intent(in), optional   :: alpha
  integer, intent(inout), optional, target :: work(:)
  integer, intent(in), optional            :: mode,data
  character, intent(in), optional          :: tran
--- a/base/comm/psb_zgather.f90
+++ b/base/comm/psb_zgather.f90
@ -52,8 +52,8 @@ subroutine  psb_zgatherm(globx, locx, desc_a, info, iroot)
  use psb_penv_mod
  implicit none
-  complex(kind(1.d0)), intent(in)    :: locx(:,:)
+  complex(psb_dpk_), intent(in)    :: locx(:,:)
-  complex(kind(1.d0)), intent(out)   :: globx(:,:)
+  complex(psb_dpk_), intent(out)   :: globx(:,:)
  type(psb_desc_type), intent(in) :: desc_a
  integer, intent(out)            :: info
  integer, intent(in), optional   :: iroot
@ -228,8 +228,8 @@ subroutine  psb_zgatherv(globx, locx, desc_a, info, iroot)
  use psb_penv_mod
  implicit none
-  complex(kind(1.d0)), intent(in)    :: locx(:)
+  complex(psb_dpk_), intent(in)    :: locx(:)
-  complex(kind(1.d0)), intent(out)   :: globx(:)
+  complex(psb_dpk_), intent(out)   :: globx(:)
  type(psb_desc_type), intent(in) :: desc_a
  integer, intent(out)            :: info
  integer, intent(in), optional   :: iroot
--- a/base/comm/psb_zhalo.f90
+++ b/base/comm/psb_zhalo.f90
@ -64,11 +64,11 @@ subroutine  psb_zhalom(x,desc_a,info,alpha,jx,ik,work,tran,mode,data)
  use psb_penv_mod
  implicit none
-  complex(kind(1.d0)), intent(inout), target   :: x(:,:)
+  complex(psb_dpk_), intent(inout), target   :: x(:,:)
  type(psb_desc_type), intent(in)           :: desc_a
  integer, intent(out)                      :: info
-  complex(kind(1.d0)), intent(in), optional    :: alpha
+  complex(psb_dpk_), intent(in), optional    :: alpha
-  complex(kind(1.d0)), optional, target        :: work(:)
+  complex(psb_dpk_), optional, target        :: work(:)
  integer, intent(in), optional             :: mode,jx,ik,data
  character, intent(in), optional           :: tran
@ -76,7 +76,7 @@ subroutine  psb_zhalom(x,desc_a,info,alpha,jx,ik,work,tran,mode,data)
  integer                  :: ictxt, np, me, &
       & err_act, m, n, iix, jjx, ix, ijx, k, maxk, nrow, imode, i,&
       & err, liwork,data_
-  complex(kind(1.d0)),pointer :: iwork(:), xp(:,:)
+  complex(psb_dpk_),pointer :: iwork(:), xp(:,:)
  character                :: tran_
  character(len=20)        :: name, ch_err
  logical                  :: aliw
@ -292,18 +292,18 @@ subroutine  psb_zhalov(x,desc_a,info,alpha,work,tran,mode,data)
  use psb_penv_mod
  implicit none
-  complex(kind(1.d0)), intent(inout)           :: x(:)
+  complex(psb_dpk_), intent(inout)           :: x(:)
  type(psb_desc_type), intent(in)           :: desc_a
  integer, intent(out)                      :: info
-  complex(kind(1.d0)), intent(in), optional    :: alpha
+  complex(psb_dpk_), intent(in), optional    :: alpha
-  complex(kind(1.d0)), target, optional        :: work(:)
+  complex(psb_dpk_), target, optional        :: work(:)
  integer, intent(in), optional             :: mode,data
  character, intent(in), optional           :: tran
  ! locals
  integer                  :: ictxt, np, me, err_act, &
       & m, n, iix, jjx, ix, ijx, nrow, imode, err, liwork,data_
-  complex(kind(1.d0)),pointer :: iwork(:)
+  complex(psb_dpk_),pointer :: iwork(:)
  character                :: tran_
  character(len=20)        :: name, ch_err
  logical                  :: aliw
--- a/base/comm/psb_zovrl.f90
+++ b/base/comm/psb_zovrl.f90
@ -74,17 +74,17 @@ subroutine  psb_zovrlm(x,desc_a,info,jx,ik,work,update,mode)
  use psb_penv_mod
  implicit none
-  complex(kind(1.d0)), intent(inout), target   :: x(:,:)
+  complex(psb_dpk_), intent(inout), target   :: x(:,:)
  type(psb_desc_type), intent(in)           :: desc_a
  integer, intent(out)                      :: info
-  complex(kind(1.d0)), optional, target        :: work(:)
+  complex(psb_dpk_), optional, target        :: work(:)
  integer, intent(in), optional             :: update,jx,ik,mode
  ! locals
  integer                  :: ictxt, np, me, &
       & err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, maxk, update_,&
       & mode_, err, liwork
-  complex(kind(1.d0)),pointer :: iwork(:), xp(:,:)
+  complex(psb_dpk_),pointer :: iwork(:), xp(:,:)
  logical                  :: do_swap
  character(len=20)        :: name, ch_err
  logical                  :: aliw
@ -278,17 +278,17 @@ subroutine  psb_zovrlv(x,desc_a,info,work,update,mode)
  use psb_penv_mod
  implicit none
-  complex(kind(1.d0)), intent(inout), target :: x(:)
+  complex(psb_dpk_), intent(inout), target :: x(:)
  type(psb_desc_type), intent(in)            :: desc_a
  integer, intent(out)                       :: info
-  complex(kind(1.d0)), optional, target      :: work(:)
+  complex(psb_dpk_), optional, target      :: work(:)
  integer, intent(in), optional              :: update,mode
  ! locals
  integer                  :: ictxt, np, me, &
       & err_act, m, n, iix, jjx, ix, ijx, nrow, ncol, k, update_,&
       & mode_, err, liwork
-  complex(kind(1.d0)),pointer :: iwork(:)
+  complex(psb_dpk_),pointer :: iwork(:)
  logical                  :: do_swap
  character(len=20)        :: name, ch_err
  logical                  :: aliw
--- a/base/comm/psb_zscatter.F90
+++ b/base/comm/psb_zscatter.F90
@ -57,8 +57,8 @@ subroutine  psb_zscatterm(globx, locx, desc_a, info, iroot)
  include 'mpif.h'
 #endif
-  complex(kind(1.d0)), intent(out)    :: locx(:,:)
+  complex(psb_dpk_), intent(out)    :: locx(:,:)
-  complex(kind(1.d0)), intent(in)     :: globx(:,:)
+  complex(psb_dpk_), intent(in)     :: globx(:,:)
  type(psb_desc_type), intent(in)  :: desc_a
  integer, intent(out)             :: info
  integer, intent(in), optional    :: iroot
@ -69,7 +69,7 @@ subroutine  psb_zscatterm(globx, locx, desc_a, info, iroot)
       & err_act, m, n, i, j, idx, nrow, iiroot, iglobx, jglobx,&
       & ilocx, jlocx, lda_locx, lda_globx, lock, globk, icomm, k, maxk, root, ilx,&
       & jlx, myrank, rootrank, c, pos
-  complex(kind(1.d0)),allocatable  :: scatterv(:)
+  complex(psb_dpk_),allocatable  :: scatterv(:)
  integer, allocatable             :: displ(:), l_t_g_all(:), all_dim(:)
  character(len=20)        :: name, ch_err
@ -289,8 +289,8 @@ subroutine  psb_zscatterv(globx, locx, desc_a, info, iroot)
  include 'mpif.h'
 #endif
-  complex(kind(1.d0)), intent(out)    :: locx(:)
+  complex(psb_dpk_), intent(out)    :: locx(:)
-  complex(kind(1.d0)), intent(in)     :: globx(:)
+  complex(psb_dpk_), intent(in)     :: globx(:)
  type(psb_desc_type), intent(in)  :: desc_a
  integer, intent(out)             :: info
  integer, intent(in), optional    :: iroot
@ -301,7 +301,7 @@ subroutine  psb_zscatterv(globx, locx, desc_a, info, iroot)
       & err_act, m, n, i, j, idx, nrow, iglobx, jglobx,&
       & ilocx, jlocx, lda_locx, lda_globx, root, k, icomm, myrank,&
       & rootrank, pos, ilx, jlx
-  complex(kind(1.d0)), allocatable  :: scatterv(:)
+  complex(psb_dpk_), allocatable  :: scatterv(:)
  integer, allocatable              :: displ(:), l_t_g_all(:), all_dim(:)
  character(len=20)        :: name, ch_err
--- a/base/internals/psi_dswapdata.F90
+++ b/base/internals/psi_dswapdata.F90
@ -38,8 +38,8 @@
 !   it is capable of pruning empty exchanges, which are very likely in out 
 !   application environment. All the variants have the same structure 
 !   In all these subroutines X may be:    I    Integer
-!                                         D    real(kind(1.d0))
+!                                         D    real(psb_dpk_)
-!                                         Z    complex(kind(1.d0))
+!                                         Z    complex(psb_dpk_)
 !   Basically the operation is as follows: on each process, we identify 
 !   sections SND(Y) and RCV(Y); then we do a send on (PACK(SND(Y)));
 !   then we receive, and we do an update with Y = UNPACK(RCV(Y)) + BETA * Y 
@ -98,8 +98,8 @@ subroutine psi_dswapdatam(flag,n,beta,y,desc_a,work,info,data)
  integer, intent(in)      :: flag, n
  integer, intent(out)     :: info
-  real(kind(1.d0))         :: y(:,:), beta
+  real(psb_dpk_)         :: y(:,:), beta
-  real(kind(1.d0)), target :: work(:)
+  real(psb_dpk_), target :: work(:)
  type(psb_desc_type),target  :: desc_a
  integer, optional        :: data
@ -171,8 +171,8 @@ subroutine psi_dswapidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work
  integer, intent(in)      :: ictxt,icomm,flag,n
  integer, intent(out)     :: info
-  real(kind(1.d0))         :: y(:,:), beta
+  real(psb_dpk_)         :: y(:,:), beta
-  real(kind(1.d0)), target :: work(:)
+  real(psb_dpk_), target :: work(:)
  integer, intent(in)      :: idx(:),totxch,totsnd, totrcv
  ! locals
@ -187,7 +187,7 @@ subroutine psi_dswapidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work
       & albf,do_send,do_recv
  logical, parameter :: usersend=.false.
-  real(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf
+  real(psb_dpk_), pointer, dimension(:) :: sndbuf, rcvbuf
 #ifdef HAVE_VOLATILE
  volatile :: sndbuf, rcvbuf
 #endif
@ -529,8 +529,8 @@ end subroutine psi_dswapidxm
 !   it is capable of pruning empty exchanges, which are very likely in out 
 !   application environment. All the variants have the same structure 
 !   In all these subroutines X may be:    I    Integer
-!                                         D    real(kind(1.d0))
+!                                         D    real(psb_dpk_)
-!                                         Z    complex(kind(1.d0))
+!                                         Z    complex(psb_dpk_)
 !   Basically the operation is as follows: on each process, we identify 
 !   sections SND(Y) and RCV(Y); then we do a SEND(PACK(SND(Y)));
 !   then we receive, and we do an update with Y = UNPACK(RCV(Y)) + BETA * Y 
@ -587,8 +587,8 @@ subroutine psi_dswapdatav(flag,beta,y,desc_a,work,info,data)
  integer, intent(in)      :: flag
  integer, intent(out)     :: info
-  real(kind(1.d0))         :: y(:), beta
+  real(psb_dpk_)         :: y(:), beta
-  real(kind(1.d0)), target :: work(:)
+  real(psb_dpk_), target :: work(:)
  type(psb_desc_type),target    :: desc_a
  integer, optional        :: data
@ -661,8 +661,8 @@ subroutine psi_dswapidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,i
  integer, intent(in)      :: ictxt,icomm,flag
  integer, intent(out)     :: info
-  real(kind(1.d0))         :: y(:), beta
+  real(psb_dpk_)         :: y(:), beta
-  real(kind(1.d0)), target :: work(:)
+  real(psb_dpk_), target :: work(:)
  integer, intent(in)      :: idx(:),totxch,totsnd, totrcv
  ! locals
@ -678,7 +678,7 @@ subroutine psi_dswapidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,i
       & albf,do_send,do_recv
  logical, parameter :: usersend=.false.
-  real(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf
+  real(psb_dpk_), pointer, dimension(:) :: sndbuf, rcvbuf
 #ifdef HAVE_VOLATILE
  volatile :: sndbuf, rcvbuf
 #endif
--- a/base/internals/psi_dswaptran.F90
+++ b/base/internals/psi_dswaptran.F90
@ -42,8 +42,8 @@
 !   it is capable of pruning empty exchanges, which are very likely in out 
 !   application environment. All the variants have the same structure 
 !   In all these subroutines X may be:    I    Integer
-!                                         D    real(kind(1.d0))
+!                                         D    real(psb_dpk_)
-!                                         Z    complex(kind(1.d0))
+!                                         Z    complex(psb_dpk_)
 !   Basically the operation is as follows: on each process, we identify 
 !   sections SND(Y) and RCV(Y); then we do a SEND(PACK(SND(Y)));
 !   then we receive, and we do an update with Y = UNPACK(RCV(Y)) + BETA * Y 
@ -101,8 +101,8 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data)
  integer, intent(in)      :: flag, n
  integer, intent(out)     :: info
-  real(kind(1.d0))         :: y(:,:), beta
+  real(psb_dpk_)         :: y(:,:), beta
-  real(kind(1.d0)), target :: work(:)
+  real(psb_dpk_), target :: work(:)
  type(psb_desc_type),target      :: desc_a
  integer, optional        :: data
@ -175,8 +175,8 @@ subroutine psi_dtranidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work
  integer, intent(in)      :: ictxt,icomm,flag,n
  integer, intent(out)     :: info
-  real(kind(1.d0))         :: y(:,:), beta
+  real(psb_dpk_)         :: y(:,:), beta
-  real(kind(1.d0)), target :: work(:)
+  real(psb_dpk_), target :: work(:)
  integer, intent(in)      :: idx(:),totxch,totsnd, totrcv
  ! locals
@ -191,7 +191,7 @@ subroutine psi_dtranidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work
       & albf,do_send,do_recv
  logical, parameter :: usersend=.false.
-  real(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf
+  real(psb_dpk_), pointer, dimension(:) :: sndbuf, rcvbuf
 #ifdef HAVE_VOLATILE
  volatile :: sndbuf, rcvbuf
 #endif
@ -531,8 +531,8 @@ end subroutine psi_dtranidxm
 !   it is capable of pruning empty exchanges, which are very likely in out 
 !   application environment. All the variants have the same structure 
 !   In all these subroutines X may be:    I    Integer
-!                                         D    real(kind(1.d0))
+!                                         D    real(psb_dpk_)
-!                                         Z    complex(kind(1.d0))
+!                                         Z    complex(psb_dpk_)
 !   Basically the operation is as follows: on each process, we identify 
 !   sections SND(Y) and RCV(Y); then we do a SEND(PACK(SND(Y)));
 !   then we receive, and we do an update with Y = UNPACK(RCV(Y)) + BETA * Y 
@ -590,8 +590,8 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data)
  integer, intent(in)      :: flag
  integer, intent(out)     :: info
-  real(kind(1.d0))         :: y(:), beta
+  real(psb_dpk_)         :: y(:), beta
-  real(kind(1.d0)), target :: work(:)
+  real(psb_dpk_), target :: work(:)
  type(psb_desc_type),target  :: desc_a
  integer, optional    :: data
@ -665,8 +665,8 @@ subroutine psi_dtranidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,i
  integer, intent(in)      :: ictxt,icomm,flag
  integer, intent(out)     :: info
-  real(kind(1.d0))         :: y(:), beta
+  real(psb_dpk_)         :: y(:), beta
-  real(kind(1.d0)), target :: work(:)
+  real(psb_dpk_), target :: work(:)
  integer, intent(in)      :: idx(:),totxch,totsnd, totrcv
  ! locals
@ -681,7 +681,7 @@ subroutine psi_dtranidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,i
       & albf,do_send,do_recv
  logical, parameter :: usersend=.false.
-  real(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf
+  real(psb_dpk_), pointer, dimension(:) :: sndbuf, rcvbuf
  character(len=20)  :: name
  info = 0
--- a/base/internals/psi_iswapdata.F90
+++ b/base/internals/psi_iswapdata.F90
@ -38,8 +38,8 @@
 !   it is capable of pruning empty exchanges, which are very likely in out 
 !   application environment. All the variants have the same structure 
 !   In all these subroutines X may be:    I    Integer
-!                                         D    real(kind(1.d0))
+!                                         D    real(psb_dpk_)
-!                                         Z    complex(kind(1.d0))
+!                                         Z    complex(psb_dpk_)
 !   Basically the operation is as follows: on each process, we identify 
 !   sections SND(Y) and RCV(Y); then we do a send on (PACK(SND(Y)));
 !   then we receive, and we do an update with Y = UNPACK(RCV(Y)) + BETA * Y 
@ -528,8 +528,8 @@ end subroutine psi_iswapidxm
 !   it is capable of pruning empty exchanges, which are very likely in out 
 !   application environment. All the variants have the same structure 
 !   In all these subroutines X may be:    I    Integer
-!                                         D    real(kind(1.d0))
+!                                         D    real(psb_dpk_)
-!                                         Z    complex(kind(1.d0))
+!                                         Z    complex(psb_dpk_)
 !   Basically the operation is as follows: on each process, we identify 
 !   sections SND(Y) and RCV(Y); then we do a SEND(PACK(SND(Y)));
 !   then we receive, and we do an update with Y = UNPACK(RCV(Y)) + BETA * Y 
--- a/base/internals/psi_iswaptran.F90
+++ b/base/internals/psi_iswaptran.F90
@ -42,8 +42,8 @@
 !   it is capable of pruning empty exchanges, which are very likely in out 
 !   application environment. All the variants have the same structure 
 !   In all these subroutines X may be:    I    Integer
-!                                         D    real(kind(1.d0))
+!                                         D    real(psb_dpk_)
-!                                         Z    complex(kind(1.d0))
+!                                         Z    complex(psb_dpk_)
 !   Basically the operation is as follows: on each process, we identify 
 !   sections SND(Y) and RCV(Y); then we do a SEND(PACK(SND(Y)));
 !   then we receive, and we do an update with Y = UNPACK(RCV(Y)) + BETA * Y 
@ -531,8 +531,8 @@ end subroutine psi_itranidxm
 !   it is capable of pruning empty exchanges, which are very likely in out 
 !   application environment. All the variants have the same structure 
 !   In all these subroutines X may be:    I    Integer
-!                                         D    real(kind(1.d0))
+!                                         D    real(psb_dpk_)
-!                                         Z    complex(kind(1.d0))
+!                                         Z    complex(psb_dpk_)
 !   Basically the operation is as follows: on each process, we identify 
 !   sections SND(Y) and RCV(Y); then we do a SEND(PACK(SND(Y)));
 !   then we receive, and we do an update with Y = UNPACK(RCV(Y)) + BETA * Y 
--- a/base/internals/psi_zswapdata.F90
+++ b/base/internals/psi_zswapdata.F90
@ -38,8 +38,8 @@
 !   it is capable of pruning empty exchanges, which are very likely in out 
 !   application environment. All the variants have the same structure 
 !   In all these subroutines X may be:    I    Integer
-!                                         D    real(kind(1.d0))
+!                                         D    real(psb_dpk_)
-!                                         Z    complex(kind(1.d0))
+!                                         Z    complex(psb_dpk_)
 !   Basically the operation is as follows: on each process, we identify 
 !   sections SND(Y) and RCV(Y); then we do a send on (PACK(SND(Y)));
 !   then we receive, and we do an update with Y = UNPACK(RCV(Y)) + BETA * Y 
@ -97,8 +97,8 @@ subroutine psi_zswapdatam(flag,n,beta,y,desc_a,work,info,data)
  integer, intent(in)      :: flag, n
  integer, intent(out)     :: info
-  complex(kind(1.d0))         :: y(:,:), beta
+  complex(psb_dpk_)         :: y(:,:), beta
-  complex(kind(1.d0)), target :: work(:)
+  complex(psb_dpk_), target :: work(:)
  type(psb_desc_type),target      :: desc_a
  integer, optional        :: data
@ -170,8 +170,8 @@ subroutine psi_zswapidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work
  integer, intent(in)      :: ictxt,icomm,flag,n
  integer, intent(out)     :: info
-  complex(kind(1.d0))         :: y(:,:), beta
+  complex(psb_dpk_)         :: y(:,:), beta
-  complex(kind(1.d0)), target :: work(:)
+  complex(psb_dpk_), target :: work(:)
  integer, intent(in)      :: idx(:),totxch,totsnd, totrcv
  ! locals
@ -186,7 +186,7 @@ subroutine psi_zswapidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work
       & albf,do_send,do_recv
  logical, parameter :: usersend=.false.
-  complex(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf
+  complex(psb_dpk_), pointer, dimension(:) :: sndbuf, rcvbuf
 #ifdef HAVE_VOLATILE
  volatile :: sndbuf, rcvbuf
 #endif
@ -528,8 +528,8 @@ end subroutine psi_zswapidxm
 !   it is capable of pruning empty exchanges, which are very likely in out 
 !   application environment. All the variants have the same structure 
 !   In all these subroutines X may be:    I    Integer
-!                                         D    real(kind(1.d0))
+!                                         D    real(psb_dpk_)
-!                                         Z    complex(kind(1.d0))
+!                                         Z    complex(psb_dpk_)
 !   Basically the operation is as follows: on each process, we identify 
 !   sections SND(Y) and RCV(Y); then we do a SEND(PACK(SND(Y)));
 !   then we receive, and we do an update with Y = UNPACK(RCV(Y)) + BETA * Y 
@ -587,8 +587,8 @@ subroutine psi_zswapdatav(flag,beta,y,desc_a,work,info,data)
  integer, intent(in)      :: flag
  integer, intent(out)     :: info
-  complex(kind(1.d0))         :: y(:), beta
+  complex(psb_dpk_)         :: y(:), beta
-  complex(kind(1.d0)), target :: work(:)
+  complex(psb_dpk_), target :: work(:)
  type(psb_desc_type),target      :: desc_a
  integer, optional        :: data
@ -661,8 +661,8 @@ subroutine psi_zswapidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,i
  integer, intent(in)      :: ictxt,icomm,flag
  integer, intent(out)     :: info
-  complex(kind(1.d0))         :: y(:), beta
+  complex(psb_dpk_)         :: y(:), beta
-  complex(kind(1.d0)), target :: work(:)
+  complex(psb_dpk_), target :: work(:)
  integer, intent(in)      :: idx(:),totxch,totsnd, totrcv
  ! locals
@ -678,7 +678,7 @@ subroutine psi_zswapidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,i
       & albf,do_send,do_recv
  logical, parameter :: usersend=.false.
-  complex(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf
+  complex(psb_dpk_), pointer, dimension(:) :: sndbuf, rcvbuf
 #ifdef HAVE_VOLATILE
  volatile :: sndbuf, rcvbuf
 #endif
--- a/base/internals/psi_zswaptran.F90
+++ b/base/internals/psi_zswaptran.F90
@ -42,8 +42,8 @@
 !   it is capable of pruning empty exchanges, which are very likely in out 
 !   application environment. All the variants have the same structure 
 !   In all these subroutines X may be:    I    Integer
-!                                         D    real(kind(1.d0))
+!                                         D    real(psb_dpk_)
-!                                         Z    complex(kind(1.d0))
+!                                         Z    complex(psb_dpk_)
 !   Basically the operation is as follows: on each process, we identify 
 !   sections SND(Y) and RCV(Y); then we do a SEND(PACK(SND(Y)));
 !   then we receive, and we do an update with Y = UNPACK(RCV(Y)) + BETA * Y 
@ -101,8 +101,8 @@ subroutine psi_zswaptranm(flag,n,beta,y,desc_a,work,info,data)
  integer, intent(in)      :: flag, n
  integer, intent(out)     :: info
-  complex(kind(1.d0))         :: y(:,:), beta
+  complex(psb_dpk_)         :: y(:,:), beta
-  complex(kind(1.d0)), target :: work(:)
+  complex(psb_dpk_), target :: work(:)
  type(psb_desc_type),target       :: desc_a
  integer, optional         :: data
@ -175,8 +175,8 @@ subroutine psi_ztranidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work
  integer, intent(in)      :: ictxt,icomm,flag,n
  integer, intent(out)     :: info
-  complex(kind(1.d0))         :: y(:,:), beta
+  complex(psb_dpk_)         :: y(:,:), beta
-  complex(kind(1.d0)), target :: work(:)
+  complex(psb_dpk_), target :: work(:)
  integer, intent(in)      :: idx(:),totxch,totsnd, totrcv
  ! locals
@ -191,7 +191,7 @@ subroutine psi_ztranidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work
       & albf,do_send,do_recv
  logical, parameter :: usersend=.false.
-  complex(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf
+  complex(psb_dpk_), pointer, dimension(:) :: sndbuf, rcvbuf
 #ifdef HAVE_VOLATILE
  volatile :: sndbuf, rcvbuf
 #endif
@ -531,8 +531,8 @@ end subroutine psi_ztranidxm
 !   it is capable of pruning empty exchanges, which are very likely in out 
 !   application environment. All the variants have the same structure 
 !   In all these subroutines X may be:    I    Integer
-!                                         D    real(kind(1.d0))
+!                                         D    real(psb_dpk_)
-!                                         Z    complex(kind(1.d0))
+!                                         Z    complex(psb_dpk_)
 !   Basically the operation is as follows: on each process, we identify 
 !   sections SND(Y) and RCV(Y); then we do a SEND(PACK(SND(Y)));
 !   then we receive, and we do an update with Y = UNPACK(RCV(Y)) + BETA * Y 
@ -590,8 +590,8 @@ subroutine psi_zswaptranv(flag,beta,y,desc_a,work,info,data)
  integer, intent(in)      :: flag
  integer, intent(out)     :: info
-  complex(kind(1.d0))         :: y(:), beta
+  complex(psb_dpk_)         :: y(:), beta
-  complex(kind(1.d0)), target :: work(:)
+  complex(psb_dpk_), target :: work(:)
  type(psb_desc_type),target  :: desc_a
  integer, optional    :: data
@ -665,8 +665,8 @@ subroutine psi_ztranidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,i
  integer, intent(in)      :: ictxt,icomm,flag
  integer, intent(out)     :: info
-  complex(kind(1.d0))         :: y(:), beta
+  complex(psb_dpk_)         :: y(:), beta
-  complex(kind(1.d0)), target :: work(:)
+  complex(psb_dpk_), target :: work(:)
  integer, intent(in)      :: idx(:),totxch,totsnd, totrcv
  ! locals
@ -681,7 +681,7 @@ subroutine psi_ztranidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,i
       & albf,do_send,do_recv
  logical, parameter :: usersend=.false.
-  complex(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf
+  complex(psb_dpk_), pointer, dimension(:) :: sndbuf, rcvbuf
 #ifdef HAVE_VOLATILE
  volatile :: sndbuf, rcvbuf
 #endif
--- a/base/modules/Makefile
+++ b/base/modules/Makefile
@ -11,7 +11,7 @@ MODULES = psb_realloc_mod.o psb_string_mod.o psb_spmat_type.o \
 LIBMOD=psb_base_mod$(.mod)
 MPFOBJS=psb_penv_mod.o 
-OBJS = error.o psb_base_mod.o $(EXTRA_COBJS)
+OBJS = error.o psb_base_mod.o $(EXTRA_COBJS) cutil.o
 LIBDIR=..
 CINCLUDES=-I. 
 FINCLUDES=$(FMFLAG)$(LIBDIR) $(FMFLAG). $(FIFLAG).
@ -33,7 +33,7 @@ psb_desc_type.o: psb_const_mod.o psb_error_mod.o psb_penv_mod.o psb_realloc_mod.
 psb_inter_desc_type.o: psb_desc_type.o psb_spmat_type.o psb_error_mod.o psb_serial_mod.o psb_comm_mod.o
 psb_check_mod.o: psb_desc_type.o
 psb_serial_mod.o: psb_spmat_type.o psb_string_mod.o psb_sort_mod.o psi_serial_mod.o
-psb_sort_mod.o: psb_error_mod.o psb_realloc_mod.o
+psb_sort_mod.o: psb_error_mod.o psb_realloc_mod.o psb_const_mod.o
 psb_methd_mod.o: psb_serial_mod.o psb_desc_type.o psb_prec_type.o
 psb_tools_mod.o: psb_spmat_type.o psb_desc_type.o psi_mod.o psb_gps_mod.o psb_inter_desc_type.o
 psb_gps_mod.o:	psb_realloc_mod.o
--- a/base/modules/cutil.c
+++ b/base/modules/cutil.c
@ -0,0 +1,28 @@
 #include <sys/time.h>
 #include <stdio.h>
 #include <string.h>
 #ifdef LowerUnderscore
 #define psi_c_diffadd   psi_c_diffadd_
 #endif
 #ifdef LowerDoubleUnderscore
 #define psi_c_diffadd   psi_c_diffadd__
 #endif
 #ifdef LowerCase
 #define psi_c_diffadd   psi_c_diffadd
 #endif
 #ifdef UpperUnderscore
 #define psi_c_diffadd   PSI_C_DIFFADD_
 #endif
 #ifdef UpperDoubleUnderscore 
 #define psi_c_diffadd   PSI_C_DIFFADD__
 #endif
 #ifdef UpperCase
 #define psi_c_diffadd   PSI_C_DIFFADD
 #endif
 void psi_c_diffadd(void *p1, void *p2, int *ret)
 {
  *ret = p2-p1;
  return;
 }
--- a/base/modules/fakempi.c
+++ b/base/modules/fakempi.c
@ -2,7 +2,7 @@
 #include <stdio.h>
 #include <string.h>
-#ifdef LowerUndescore
+#ifdef LowerUnderscore
 #define mpi_wtime       mpi_wtime_
 #define mpi_send        mpi_send_
 #define mpi_irecv       mpi_irecv_
@ -12,7 +12,7 @@
 #define mpi_allgather   mpi_allgather_
 #define mpi_allgatherv  mpi_allgatherv_
 #endif
-#ifdef LowerDoubleUndescore
+#ifdef LowerDoubleUnderscore
 #define mpi_wtime       mpi_wtime__
 #define mpi_send        mpi_send__
 #define mpi_irecv       mpi_irecv__
@ -32,7 +32,7 @@
 #define mpi_allgather   mpi_allgather
 #define mpi_allgatherv  mpi_allgatherv
 #endif
-#ifdef UpperUndescore
+#ifdef UpperUnderscore
 #define mpi_wtime       MPI_WTIME_
 #define mpi_send        MPI_SEND_
 #define mpi_irecv       MPI_IRECV_
@ -42,7 +42,7 @@
 #define mpi_allgather   MPI_ALLGATHER_
 #define mpi_allgatherv  MPI_ALLGATHERV_
 #endif
-#ifdef UpperDoubleUndescore
+#ifdef UpperDoubleUnderscore 
 #define mpi_wtime       MPI_WTIME__
 #define mpi_send        MPI_SEND__
 #define mpi_irecv       MPI_IRECV__
--- a/base/modules/psb_comm_mod.f90
+++ b/base/modules/psb_comm_mod.f90
@ -34,18 +34,18 @@ module psb_comm_mod
  interface psb_ovrl
    subroutine  psb_dovrlm(x,desc_a,info,jx,ik,work,update,mode)
      use psb_descriptor_type
-      real(kind(1.d0)), intent(inout)           :: x(:,:)
+      real(psb_dpk_), intent(inout)           :: x(:,:)
      type(psb_desc_type), intent(in)           :: desc_a
      integer, intent(out)                      :: info
-      real(kind(1.d0)), intent(inout), optional :: work(:)
+      real(psb_dpk_), intent(inout), optional :: work(:)
      integer, intent(in), optional             :: update,jx,ik,mode
    end subroutine psb_dovrlm
    subroutine  psb_dovrlv(x,desc_a,info,work,update,mode)
      use psb_descriptor_type
-      real(kind(1.d0)), intent(inout)           :: x(:)
+      real(psb_dpk_), intent(inout)           :: x(:)
      type(psb_desc_type), intent(in)           :: desc_a
      integer, intent(out)                      :: info
-      real(kind(1.d0)), intent(inout), optional :: work(:)
+      real(psb_dpk_), intent(inout), optional :: work(:)
      integer, intent(in), optional             :: update,mode
    end subroutine psb_dovrlv
    subroutine  psb_iovrlm(x,desc_a,info,jx,ik,work,update,mode)
@ -53,7 +53,7 @@ module psb_comm_mod
      integer,          intent(inout)           :: x(:,:)
      type(psb_desc_type), intent(in)           :: desc_a
      integer, intent(out)                      :: info
-      real(kind(1.d0)), intent(inout), optional :: work(:)
+      real(psb_dpk_), intent(inout), optional :: work(:)
      integer, intent(in), optional             :: update,jx,ik,mode
    end subroutine psb_iovrlm
    subroutine  psb_iovrlv(x,desc_a,info,work,update,mode)
@ -61,23 +61,23 @@ module psb_comm_mod
      integer, intent(inout)                    :: x(:)
      type(psb_desc_type), intent(in)           :: desc_a
      integer, intent(out)                      :: info
-      real(kind(1.d0)), intent(inout), optional :: work(:)
+      real(psb_dpk_), intent(inout), optional :: work(:)
      integer, intent(in), optional             :: update,mode
    end subroutine psb_iovrlv
    subroutine  psb_zovrlm(x,desc_a,info,jx,ik,work,update,mode)
      use psb_descriptor_type
-      complex(kind(1.d0)), intent(inout)        :: x(:,:)
+      complex(psb_dpk_), intent(inout)        :: x(:,:)
      type(psb_desc_type), intent(in)           :: desc_a
      integer, intent(out)                      :: info
-      complex(kind(1.d0)), intent(inout), optional :: work(:)
+      complex(psb_dpk_), intent(inout), optional :: work(:)
      integer, intent(in), optional             :: update,jx,ik,mode
    end subroutine psb_zovrlm
    subroutine  psb_zovrlv(x,desc_a,info,work,update,mode)
      use psb_descriptor_type
-      complex(kind(1.d0)), intent(inout)        :: x(:)
+      complex(psb_dpk_), intent(inout)        :: x(:)
      type(psb_desc_type), intent(in)           :: desc_a
      integer, intent(out)                      :: info
-      complex(kind(1.d0)), intent(inout), optional :: work(:)
+      complex(psb_dpk_), intent(inout), optional :: work(:)
      integer, intent(in), optional             :: update,mode
    end subroutine psb_zovrlv
  end interface
@ -85,21 +85,21 @@ module psb_comm_mod
  interface psb_halo
    subroutine  psb_dhalom(x,desc_a,info,alpha,jx,ik,work,tran,mode,data)
      use psb_descriptor_type
-      real(kind(1.d0)), intent(inout)           :: x(:,:)
+      real(psb_dpk_), intent(inout)           :: x(:,:)
      type(psb_desc_type), intent(in)           :: desc_a
      integer, intent(out)                      :: info
-      real(kind(1.d0)), intent(in), optional    :: alpha
+      real(psb_dpk_), intent(in), optional    :: alpha
-      real(kind(1.d0)), target, optional        :: work(:)
+      real(psb_dpk_), target, optional        :: work(:)
      integer, intent(in), optional             :: mode,jx,ik,data
      character, intent(in), optional           :: tran
    end subroutine psb_dhalom
    subroutine  psb_dhalov(x,desc_a,info,alpha,work,tran,mode,data)
      use psb_descriptor_type
-      real(kind(1.d0)), intent(inout)           :: x(:)
+      real(psb_dpk_), intent(inout)           :: x(:)
      type(psb_desc_type), intent(in)           :: desc_a
      integer, intent(out)                      :: info
-      real(kind(1.d0)), intent(in), optional    :: alpha
+      real(psb_dpk_), intent(in), optional    :: alpha
-      real(kind(1.d0)), target, optional :: work(:)
+      real(psb_dpk_), target, optional :: work(:)
      integer, intent(in), optional             :: mode,data
      character, intent(in), optional           :: tran
    end subroutine psb_dhalov
@ -108,7 +108,7 @@ module psb_comm_mod
      integer, intent(inout) :: x(:,:)
      type(psb_desc_type), intent(in)        :: desc_a
      integer, intent(out)                   :: info
-      real(kind(1.d0)), intent(in), optional :: alpha
+      real(psb_dpk_), intent(in), optional :: alpha
      integer, intent(inout), optional       :: work(:)
      integer, intent(in), optional          :: mode,jx,ik,data
      character, intent(in), optional        :: tran
@ -118,28 +118,28 @@ module psb_comm_mod
      integer, intent(inout)                 :: x(:)
      type(psb_desc_type), intent(in)        :: desc_a
      integer, intent(out)                   :: info
-      real(kind(1.d0)), intent(in), optional :: alpha
+      real(psb_dpk_), intent(in), optional :: alpha
      integer, intent(inout), optional       :: work(:)
      integer, intent(in), optional          :: mode,data
      character, intent(in), optional        :: tran
    end subroutine psb_ihalov
    subroutine  psb_zhalom(x,desc_a,info,alpha,jx,ik,work,tran,mode,data)
      use psb_descriptor_type
-      complex(kind(1.d0)), intent(inout)        :: x(:,:)
+      complex(psb_dpk_), intent(inout)        :: x(:,:)
      type(psb_desc_type), intent(in)           :: desc_a
      integer, intent(out)                      :: info
-      complex(kind(1.d0)), intent(in), optional :: alpha
+      complex(psb_dpk_), intent(in), optional :: alpha
-      complex(kind(1.d0)), target, optional     :: work(:)
+      complex(psb_dpk_), target, optional     :: work(:)
      integer, intent(in), optional             :: mode,jx,ik,data
      character, intent(in), optional           :: tran
    end subroutine psb_zhalom
    subroutine  psb_zhalov(x,desc_a,info,alpha,work,tran,mode,data)
      use psb_descriptor_type
-      complex(kind(1.d0)), intent(inout)        :: x(:)
+      complex(psb_dpk_), intent(inout)        :: x(:)
      type(psb_desc_type), intent(in)           :: desc_a
      integer, intent(out)                      :: info
-      complex(kind(1.d0)), intent(in), optional :: alpha
+      complex(psb_dpk_), intent(in), optional :: alpha
-      complex(kind(1.d0)), target, optional     :: work(:)
+      complex(psb_dpk_), target, optional     :: work(:)
      integer, intent(in), optional             :: mode,data
      character, intent(in), optional           :: tran
    end subroutine psb_zhalov
@ -149,32 +149,32 @@ module psb_comm_mod
  interface psb_dscatter
    subroutine  psb_dscatterm(globx, locx, desc_a, info, root)
      use psb_descriptor_type
-      real(kind(1.d0)), intent(out)    :: locx(:,:)
+      real(psb_dpk_), intent(out)    :: locx(:,:)
-      real(kind(1.d0)), intent(in)     :: globx(:,:)
+      real(psb_dpk_), intent(in)     :: globx(:,:)
      type(psb_desc_type), intent(in)  :: desc_a
      integer, intent(out)             :: info
      integer, intent(in), optional    :: root
    end subroutine psb_dscatterm
    subroutine  psb_dscatterv(globx, locx, desc_a, info, root)
      use psb_descriptor_type
-      real(kind(1.d0)), intent(out)    :: locx(:)
+      real(psb_dpk_), intent(out)    :: locx(:)
-      real(kind(1.d0)), intent(in)     :: globx(:)
+      real(psb_dpk_), intent(in)     :: globx(:)
      type(psb_desc_type), intent(in)  :: desc_a
      integer, intent(out)             :: info
      integer, intent(in), optional    :: root
    end subroutine psb_dscatterv
    subroutine  psb_zscatterm(globx, locx, desc_a, info, root)
      use psb_descriptor_type
-      complex(kind(1.d0)), intent(out) :: locx(:,:)
+      complex(psb_dpk_), intent(out) :: locx(:,:)
-      complex(kind(1.d0)), intent(in)  :: globx(:,:)
+      complex(psb_dpk_), intent(in)  :: globx(:,:)
      type(psb_desc_type), intent(in)  :: desc_a
      integer, intent(out)             :: info
      integer, intent(in), optional    :: root
    end subroutine psb_zscatterm
    subroutine  psb_zscatterv(globx, locx, desc_a, info, root)
      use psb_descriptor_type
-      complex(kind(1.d0)), intent(out) :: locx(:)
+      complex(psb_dpk_), intent(out) :: locx(:)
-      complex(kind(1.d0)), intent(in)  :: globx(:)
+      complex(psb_dpk_), intent(in)  :: globx(:)
      type(psb_desc_type), intent(in)  :: desc_a
      integer, intent(out)             :: info
      integer, intent(in), optional    :: root
@ -216,32 +216,32 @@ module psb_comm_mod
    end subroutine psb_igatherv
    subroutine  psb_dgatherm(globx, locx, desc_a, info, root)
      use psb_descriptor_type
-      real(kind(1.d0)), intent(in)    :: locx(:,:)
+      real(psb_dpk_), intent(in)    :: locx(:,:)
-      real(kind(1.d0)), intent(out)   :: globx(:,:)
+      real(psb_dpk_), intent(out)   :: globx(:,:)
      type(psb_desc_type), intent(in) :: desc_a
      integer, intent(out)            :: info
      integer, intent(in), optional   :: root
    end subroutine psb_dgatherm
    subroutine  psb_dgatherv(globx, locx, desc_a, info, root)
      use psb_descriptor_type
-      real(kind(1.d0)), intent(in)    :: locx(:)
+      real(psb_dpk_), intent(in)    :: locx(:)
-      real(kind(1.d0)), intent(out)   :: globx(:)
+      real(psb_dpk_), intent(out)   :: globx(:)
      type(psb_desc_type), intent(in) :: desc_a
      integer, intent(out)            :: info
      integer, intent(in), optional   :: root
    end subroutine psb_dgatherv
    subroutine  psb_zgatherm(globx, locx, desc_a, info, root)
      use psb_descriptor_type
-      complex(kind(1.d0)), intent(in)  :: locx(:,:)
+      complex(psb_dpk_), intent(in)  :: locx(:,:)
-      complex(kind(1.d0)), intent(out) :: globx(:,:)
+      complex(psb_dpk_), intent(out) :: globx(:,:)
      type(psb_desc_type), intent(in)  :: desc_a
      integer, intent(out)             :: info
      integer, intent(in), optional    :: root
    end subroutine psb_zgatherm
    subroutine  psb_zgatherv(globx, locx, desc_a, info, root)
      use psb_descriptor_type
-      complex(kind(1.d0)), intent(in)  :: locx(:)
+      complex(psb_dpk_), intent(in)  :: locx(:)
-      complex(kind(1.d0)), intent(out) :: globx(:)
+      complex(psb_dpk_), intent(out) :: globx(:)
      type(psb_desc_type), intent(in)  :: desc_a
      integer, intent(out)             :: info
      integer, intent(in), optional    :: root
--- a/base/modules/psb_const_mod.f90
+++ b/base/modules/psb_const_mod.f90
@ -31,15 +31,21 @@
 !!$  
 module psb_const_mod
  !
  ! This must be the kind parameter corresponding to MPI_DOUBLE_PRECISION.
  !
  integer, parameter           :: psb_dpk_ = kind(1.d0)
  integer                      :: psb_sizeof_dp, psb_sizeof_int
  !
  !     Handy & miscellaneous constants
  !
-  integer, parameter             :: izero=0, ione=1
+  integer, parameter           :: izero=0, ione=1
-  integer, parameter             :: itwo=2, ithree=3,mone=-1, psb_root_=0
+  integer, parameter           :: itwo=2, ithree=3,mone=-1, psb_root_=0
-  real(kind(1.d0)), parameter    :: dzero=0.d0, done=1.d0
+  real(psb_dpk_), parameter    :: dzero=0.d0, done=1.d0
-  complex(kind(1.d0)), parameter :: zzero=(0.d0,0.0d0)
+  complex(psb_dpk_), parameter :: zzero=(0.d0,0.0d0)
-  complex(kind(1.d0)), parameter :: zone=(1.d0,0.0d0)
+  complex(psb_dpk_), parameter :: zone=(1.d0,0.0d0)
-  real(kind(1.d0)), parameter    :: epstol=1.d-32, psb_percent_=0.7
+  real(psb_dpk_), parameter    :: epstol=1.d-32, psb_percent_=0.7
-  character, parameter           :: psb_all_='A',  psb_topdef_=' '
+  character, parameter         :: psb_all_='A',  psb_topdef_=' '
 end module psb_const_mod
--- a/base/modules/psb_desc_type.f90
+++ b/base/modules/psb_desc_type.f90
@ -382,20 +382,20 @@ contains
    integer, external :: SizeofPairSearchTree
    val = 0
-    if (allocated(desc%matrix_data))  val = val + 4*size(desc%matrix_data)
+    if (allocated(desc%matrix_data))  val = val + psb_sizeof_int*size(desc%matrix_data)
-    if (allocated(desc%halo_index))   val = val + 4*size(desc%halo_index)
+    if (allocated(desc%halo_index))   val = val + psb_sizeof_int*size(desc%halo_index)
-    if (allocated(desc%ext_index))    val = val + 4*size(desc%ext_index)
+    if (allocated(desc%ext_index))    val = val + psb_sizeof_int*size(desc%ext_index)
-    if (allocated(desc%bnd_elem))     val = val + 4*size(desc%bnd_elem)
+    if (allocated(desc%bnd_elem))     val = val + psb_sizeof_int*size(desc%bnd_elem)
-    if (allocated(desc%ovrlap_index)) val = val + 4*size(desc%ovrlap_index)
+    if (allocated(desc%ovrlap_index)) val = val + psb_sizeof_int*size(desc%ovrlap_index)
-    if (allocated(desc%ovrlap_elem))  val = val + 4*size(desc%ovrlap_elem)
+    if (allocated(desc%ovrlap_elem))  val = val + psb_sizeof_int*size(desc%ovrlap_elem)
-    if (allocated(desc%ovr_mst_idx))  val = val + 4*size(desc%ovr_mst_idx)
+    if (allocated(desc%ovr_mst_idx))  val = val + psb_sizeof_int*size(desc%ovr_mst_idx)
-    if (allocated(desc%loc_to_glob))  val = val + 4*size(desc%loc_to_glob)    
+    if (allocated(desc%loc_to_glob))  val = val + psb_sizeof_int*size(desc%loc_to_glob)    
-    if (allocated(desc%glob_to_loc))  val = val + 4*size(desc%glob_to_loc)
+    if (allocated(desc%glob_to_loc))  val = val + psb_sizeof_int*size(desc%glob_to_loc)
-    if (allocated(desc%hashv))        val = val + 4*size(desc%hashv)
+    if (allocated(desc%hashv))        val = val + psb_sizeof_int*size(desc%hashv)
-    if (allocated(desc%glb_lc))       val = val + 4*size(desc%glb_lc)
+    if (allocated(desc%glb_lc))       val = val + psb_sizeof_int*size(desc%glb_lc)
-    if (allocated(desc%lprm))         val = val + 4*size(desc%lprm)
+    if (allocated(desc%lprm))         val = val + psb_sizeof_int*size(desc%lprm)
-    if (allocated(desc%idx_space))    val = val + 4*size(desc%idx_space)
+    if (allocated(desc%idx_space))    val = val + psb_sizeof_int*size(desc%idx_space)
-    if (allocated(desc%ptree))        val = val + 4*size(desc%ptree) +&
+    if (allocated(desc%ptree))        val = val + psb_sizeof_int*size(desc%ptree) +&
         &                                  SizeofPairSearchTree(desc%ptree)
    psb_cd_sizeof = val
--- a/base/modules/psb_error_mod.F90
+++ b/base/modules/psb_error_mod.F90
@ -62,7 +62,7 @@ module psb_error_mod
    character(len=20)        ::   routine=''       
    !  array of integer data to complete the error msg   
    integer,dimension(5)     ::   i_err_data=0     
-    !     real(kind(1.d0))(dim=10) ::   r_err_data=0.d0    !  array of real data to complete the error msg
+    !     real(psb_dpk_)(dim=10) ::   r_err_data=0.d0    !  array of real data to complete the error msg
    !     complex(dim=10)          ::   c_err_data=0.c0    !  array of complex data to complete the error msg
    !  array of character data to complete the error msg
    character(len=40)        ::   a_err_data=''      
--- a/base/modules/psb_inter_desc_type.f90
+++ b/base/modules/psb_inter_desc_type.f90
@ -479,14 +479,14 @@ contains
    use psi_mod
    implicit none 
    type(psb_inter_desc_type), intent(in) :: desc
-    real(kind(1.d0)), intent(in)     :: alpha,beta
+    real(psb_dpk_), intent(in)     :: alpha,beta
-    real(kind(1.d0)), intent(inout)  :: x(:)
+    real(psb_dpk_), intent(inout)  :: x(:)
-    real(kind(1.d0)), intent(out)    :: y(:)
+    real(psb_dpk_), intent(out)    :: y(:)
    integer, intent(out)             :: info 
-    real(kind(1.d0)), optional       :: work(:)
+    real(psb_dpk_), optional       :: work(:)
    !
-    real(kind(1.d0)), allocatable :: xt(:)
+    real(psb_dpk_), allocatable :: xt(:)
    integer                       :: itsz, i, j,totxch,totsnd,totrcv,&
         &  map_kind, map_data
    character(len=20), parameter  :: name='psb_forward_map'
@ -553,14 +553,14 @@ contains
    use psi_mod
    implicit none 
    type(psb_inter_desc_type), intent(in) :: desc
-    real(kind(1.d0)), intent(in)     :: alpha,beta
+    real(psb_dpk_), intent(in)     :: alpha,beta
-    real(kind(1.d0)), intent(inout)  :: x(:)
+    real(psb_dpk_), intent(inout)  :: x(:)
-    real(kind(1.d0)), intent(out)    :: y(:)
+    real(psb_dpk_), intent(out)    :: y(:)
    integer, intent(out)             :: info 
-    real(kind(1.d0)), optional       :: work(:)
+    real(psb_dpk_), optional       :: work(:)
    !
-    real(kind(1.d0)), allocatable :: xt(:)
+    real(psb_dpk_), allocatable :: xt(:)
    integer                       :: itsz, i, j,totxch,totsnd,totrcv,&
         & map_kind, map_data
    character(len=20), parameter  :: name='psb_backward_map'
@ -624,14 +624,14 @@ contains
    use psi_mod
    implicit none 
    type(psb_inter_desc_type), intent(in) :: desc
-    complex(kind(1.d0)), intent(in)       :: alpha,beta
+    complex(psb_dpk_), intent(in)       :: alpha,beta
-    complex(kind(1.d0)), intent(inout)    :: x(:)
+    complex(psb_dpk_), intent(inout)    :: x(:)
-    complex(kind(1.d0)), intent(out)      :: y(:)
+    complex(psb_dpk_), intent(out)      :: y(:)
    integer, intent(out)                  :: info 
-    complex(kind(1.d0)), optional         :: work(:)
+    complex(psb_dpk_), optional         :: work(:)
    !
-    complex(kind(1.d0)), allocatable :: xt(:)
+    complex(psb_dpk_), allocatable :: xt(:)
    integer                       :: itsz, i, j,totxch,totsnd,totrcv,&
         & map_kind, map_data
    character(len=20), parameter  :: name='psb_forward_map'
@ -695,14 +695,14 @@ contains
    use psi_mod
    implicit none 
    type(psb_inter_desc_type), intent(in) :: desc
-    complex(kind(1.d0)), intent(in)       :: alpha,beta
+    complex(psb_dpk_), intent(in)       :: alpha,beta
-    complex(kind(1.d0)), intent(inout)    :: x(:)
+    complex(psb_dpk_), intent(inout)    :: x(:)
-    complex(kind(1.d0)), intent(out)      :: y(:)
+    complex(psb_dpk_), intent(out)      :: y(:)
    integer, intent(out)                  :: info 
-    complex(kind(1.d0)), optional         :: work(:)
+    complex(psb_dpk_), optional         :: work(:)
    !
-    complex(kind(1.d0)), allocatable :: xt(:)
+    complex(psb_dpk_), allocatable :: xt(:)
    integer                       :: itsz, i, j,totxch,totsnd,totrcv,&
         & map_kind, map_data
    character(len=20), parameter  :: name='psb_backward_map'
--- a/base/modules/psb_penv_mod.F90
+++ b/base/modules/psb_penv_mod.F90
--- a/base/modules/psb_psblas_mod.f90
+++ b/base/modules/psb_psblas_mod.f90
@ -34,30 +34,30 @@ module psb_psblas_mod
  interface psb_gedot
    function psb_ddotv(x, y, desc_a,info) 
      use psb_descriptor_type
-      real(kind(1.d0))                   :: psb_ddotv
+      real(psb_dpk_)                   :: psb_ddotv
-      real(kind(1.d0)), intent(in)       :: x(:), y(:)
+      real(psb_dpk_), intent(in)       :: x(:), y(:)
      type(psb_desc_type), intent(in)    :: desc_a
      integer, intent(out)               :: info
    end function psb_ddotv
    function psb_ddot(x, y, desc_a, info, jx, jy) 
      use psb_descriptor_type
-      real(kind(1.d0))                   :: psb_ddot
+      real(psb_dpk_)                   :: psb_ddot
-      real(kind(1.d0)), intent(in)       :: x(:,:), y(:,:)
+      real(psb_dpk_), intent(in)       :: x(:,:), y(:,:)
      type(psb_desc_type), intent(in)    :: desc_a
      integer, optional, intent(in)      :: jx, jy
      integer, intent(out)               :: info
    end function psb_ddot
    function psb_zdotv(x, y, desc_a,info) 
      use psb_descriptor_type
-      complex(kind(1.d0))                :: psb_zdotv
+      complex(psb_dpk_)                :: psb_zdotv
-      complex(kind(1.d0)), intent(in)    :: x(:), y(:)
+      complex(psb_dpk_), intent(in)    :: x(:), y(:)
      type(psb_desc_type), intent(in)    :: desc_a
      integer, intent(out)               :: info
    end function psb_zdotv
    function psb_zdot(x, y, desc_a, info, jx, jy) 
      use psb_descriptor_type
-      complex(kind(1.d0))                :: psb_zdot
+      complex(psb_dpk_)                :: psb_zdot
-      complex(kind(1.d0)), intent(in)    :: x(:,:), y(:,:)
+      complex(psb_dpk_), intent(in)    :: x(:,:), y(:,:)
      type(psb_desc_type), intent(in)    :: desc_a
      integer, optional, intent(in)      :: jx, jy
      integer, intent(out)               :: info
@ -67,29 +67,29 @@ module psb_psblas_mod
  interface psb_gedots
    subroutine  psb_ddotvs(res,x, y, desc_a, info) 
      use psb_descriptor_type
-      real(kind(1.d0)), intent(out)      :: res
+      real(psb_dpk_), intent(out)      :: res
-      real(kind(1.d0)), intent(in)       :: x(:), y(:)
+      real(psb_dpk_), intent(in)       :: x(:), y(:)
      type(psb_desc_type), intent(in)    :: desc_a
      integer, intent(out)               :: info
    end subroutine psb_ddotvs
    subroutine  psb_dmdots(res,x, y, desc_a,info) 
      use psb_descriptor_type
-      real(kind(1.d0)), intent(out)      :: res(:)
+      real(psb_dpk_), intent(out)      :: res(:)
-      real(kind(1.d0)), intent(in)       :: x(:,:), y(:,:)
+      real(psb_dpk_), intent(in)       :: x(:,:), y(:,:)
      type(psb_desc_type), intent(in)    :: desc_a
      integer, intent(out)               :: info
    end subroutine psb_dmdots
    subroutine  psb_zdotvs(res,x, y, desc_a, info) 
      use psb_descriptor_type
-      complex(kind(1.d0)), intent(out)      :: res
+      complex(psb_dpk_), intent(out)      :: res
-      complex(kind(1.d0)), intent(in)       :: x(:), y(:)
+      complex(psb_dpk_), intent(in)       :: x(:), y(:)
      type(psb_desc_type), intent(in)    :: desc_a
      integer, intent(out)               :: info
    end subroutine psb_zdotvs
    subroutine  psb_zmdots(res,x, y, desc_a,info) 
      use psb_descriptor_type
-      complex(kind(1.d0)), intent(out)      :: res(:)
+      complex(psb_dpk_), intent(out)      :: res(:)
-      complex(kind(1.d0)), intent(in)       :: x(:,:), y(:,:)
+      complex(psb_dpk_), intent(in)       :: x(:,:), y(:,:)
      type(psb_desc_type), intent(in)    :: desc_a
      integer, intent(out)               :: info
    end subroutine psb_zmdots
@ -99,18 +99,18 @@ module psb_psblas_mod
     subroutine psb_daxpbyv(alpha, x, beta, y,&
 	  & desc_a, info)
       use psb_descriptor_type
-       real(kind(1.d0)), intent (in)       ::  x(:)
+       real(psb_dpk_), intent (in)       ::  x(:)
-       real(kind(1.d0)), intent (inout)    ::  y(:)
+       real(psb_dpk_), intent (inout)    ::  y(:)
-       real(kind(1.d0)), intent (in)       :: alpha, beta
+       real(psb_dpk_), intent (in)       :: alpha, beta
       type(psb_desc_type), intent (in)    :: desc_a
       integer, intent(out)                :: info
     end subroutine psb_daxpbyv
     subroutine psb_daxpby(alpha, x, beta, y,&
 	  & desc_a, info, n, jx, jy)
       use psb_descriptor_type
-       real(kind(1.d0)), intent (in)       ::  x(:,:)
+       real(psb_dpk_), intent (in)       ::  x(:,:)
-       real(kind(1.d0)), intent (inout)    ::  y(:,:)
+       real(psb_dpk_), intent (inout)    ::  y(:,:)
-       real(kind(1.d0)), intent (in)       ::  alpha, beta
+       real(psb_dpk_), intent (in)       ::  alpha, beta
       type(psb_desc_type), intent (in)    :: desc_a
       integer, optional :: n, jx, jy
       integer, intent(out)                :: info
@ -118,18 +118,18 @@ module psb_psblas_mod
     subroutine psb_zaxpbyv(alpha, x, beta, y,&
 	  & desc_a, info)
       use psb_descriptor_type
-       complex(kind(1.d0)), intent (in)       ::  x(:)
+       complex(psb_dpk_), intent (in)       ::  x(:)
-       complex(kind(1.d0)), intent (inout)    ::  y(:)
+       complex(psb_dpk_), intent (inout)    ::  y(:)
-       complex(kind(1.d0)), intent (in)       :: alpha, beta
+       complex(psb_dpk_), intent (in)       :: alpha, beta
       type(psb_desc_type), intent (in)    :: desc_a
       integer, intent(out)                :: info
     end subroutine psb_zaxpbyv
     subroutine psb_zaxpby(alpha, x, beta, y,&
 	  & desc_a, info, n, jx, jy)
       use psb_descriptor_type
-       complex(kind(1.d0)), intent (in)       ::  x(:,:)
+       complex(psb_dpk_), intent (in)       ::  x(:,:)
-       complex(kind(1.d0)), intent (inout)    ::  y(:,:)
+       complex(psb_dpk_), intent (inout)    ::  y(:,:)
-       complex(kind(1.d0)), intent (in)       ::  alpha, beta
+       complex(psb_dpk_), intent (in)       ::  alpha, beta
       type(psb_desc_type), intent (in)    :: desc_a
       integer, optional :: n, jx, jy
       integer, intent(out)                :: info
@ -139,31 +139,31 @@ module psb_psblas_mod
  interface psb_geamax
     function psb_damax(x, desc_a, info, jx)
       use psb_descriptor_type
-       real(kind(1.d0))   psb_damax
+       real(psb_dpk_)   psb_damax
-       real(kind(1.d0)), intent (in)       :: x(:,:)
+       real(psb_dpk_), intent (in)       :: x(:,:)
       type(psb_desc_type), intent (in)    :: desc_a
       integer, optional, intent (in)      :: jx
       integer, intent(out)                :: info
     end function psb_damax
     function psb_damaxv(x, desc_a,info)
       use psb_descriptor_type
-       real(kind(1.d0)) psb_damaxv
+       real(psb_dpk_) psb_damaxv
-       real(kind(1.d0)), intent (in)       :: x(:)
+       real(psb_dpk_), intent (in)       :: x(:)
       type(psb_desc_type), intent (in)    :: desc_a
       integer, intent(out)                :: info
     end function psb_damaxv
     function psb_zamax(x, desc_a, info, jx)
       use psb_descriptor_type
-       real(kind(1.d0))   psb_zamax
+       real(psb_dpk_)   psb_zamax
-       complex(kind(1.d0)), intent (in)       :: x(:,:)
+       complex(psb_dpk_), intent (in)       :: x(:,:)
       type(psb_desc_type), intent (in)    :: desc_a
       integer, optional, intent (in)      :: jx
       integer, intent(out)                :: info
     end function psb_zamax
     function psb_zamaxv(x, desc_a,info)
       use psb_descriptor_type
-       real(kind(1.d0)) psb_zamaxv
+       real(psb_dpk_) psb_zamaxv
-       complex(kind(1.d0)), intent (in)       :: x(:)
+       complex(psb_dpk_), intent (in)       :: x(:)
       type(psb_desc_type), intent (in)    :: desc_a
       integer, intent(out)                :: info
     end function psb_zamaxv
@ -172,30 +172,30 @@ module psb_psblas_mod
  interface psb_geamaxs
     subroutine  psb_damaxvs(res,x,desc_a,info)
       use psb_descriptor_type
-       real(kind(1.d0)), intent (out)      :: res
+       real(psb_dpk_), intent (out)      :: res
-       real(kind(1.d0)), intent (in)       :: x(:)
+       real(psb_dpk_), intent (in)       :: x(:)
       type(psb_desc_type), intent (in) :: desc_a
       integer, intent(out)               :: info
     end subroutine psb_damaxvs
     subroutine  psb_dmamax(res,x,desc_a,info,jx)
       use psb_descriptor_type
-       real(kind(1.d0)), intent (out)      :: res(:)
+       real(psb_dpk_), intent (out)      :: res(:)
-       real(kind(1.d0)), intent (in)       :: x(:,:)
+       real(psb_dpk_), intent (in)       :: x(:,:)
       type(psb_desc_type), intent (in)    :: desc_a
       integer, intent(out)                :: info
       integer, optional                   :: jx
     end subroutine psb_dmamax
     subroutine  psb_zamaxvs(res,x,desc_a,info)
       use psb_descriptor_type
-       real(kind(1.d0)), intent (out)      :: res
+       real(psb_dpk_), intent (out)      :: res
-       complex(kind(1.d0)), intent (in)       :: x(:)
+       complex(psb_dpk_), intent (in)       :: x(:)
       type(psb_desc_type), intent (in) :: desc_a
       integer, intent(out)               :: info
     end subroutine psb_zamaxvs
     subroutine  psb_zmamax(res,x,desc_a,info,jx)
       use psb_descriptor_type
-       real(kind(1.d0)), intent (out)      :: res(:)
+       real(psb_dpk_), intent (out)      :: res(:)
-       complex(kind(1.d0)), intent (in)       :: x(:,:)
+       complex(psb_dpk_), intent (in)       :: x(:,:)
       type(psb_desc_type), intent (in)    :: desc_a
       integer, intent(out)                :: info
       integer, optional                   :: jx
@ -205,31 +205,31 @@ module psb_psblas_mod
  interface psb_geasum
     function psb_dasum(x, desc_a, info, jx)
       use psb_descriptor_type
-       real(kind(1.d0))   psb_dasum
+       real(psb_dpk_)   psb_dasum
-       real(kind(1.d0)), intent (in)       :: x(:,:)
+       real(psb_dpk_), intent (in)       :: x(:,:)
       type(psb_desc_type), intent (in)    :: desc_a
       integer, optional, intent (in)      :: jx
       integer, intent(out)                :: info
     end function psb_dasum
     function psb_dasumv(x, desc_a, info)
       use psb_descriptor_type
-       real(kind(1.d0)) psb_dasumv
+       real(psb_dpk_) psb_dasumv
-       real(kind(1.d0)), intent (in)       :: x(:)
+       real(psb_dpk_), intent (in)       :: x(:)
       type(psb_desc_type), intent (in)    :: desc_a
       integer, intent(out)                :: info
     end function psb_dasumv
     function psb_zasum(x, desc_a, info, jx)
       use psb_descriptor_type
-       real(kind(1.d0))   psb_zasum
+       real(psb_dpk_)   psb_zasum
-       complex(kind(1.d0)), intent (in)       :: x(:,:)
+       complex(psb_dpk_), intent (in)       :: x(:,:)
       type(psb_desc_type), intent (in)    :: desc_a
       integer, optional, intent (in)      :: jx
       integer, intent(out)                :: info
     end function psb_zasum
     function psb_zasumv(x, desc_a, info)
       use psb_descriptor_type
-       real(kind(1.d0)) psb_zasumv
+       real(psb_dpk_) psb_zasumv
-       complex(kind(1.d0)), intent (in)       :: x(:)
+       complex(psb_dpk_), intent (in)       :: x(:)
       type(psb_desc_type), intent (in)    :: desc_a
       integer, intent(out)                :: info
     end function psb_zasumv
@ -238,29 +238,29 @@ module psb_psblas_mod
  interface psb_geasums
     subroutine  psb_dasumvs(res,x,desc_a,info)
       use psb_descriptor_type
-       real(kind(1.d0)), intent (out)      :: res
+       real(psb_dpk_), intent (out)      :: res
-       real(kind(1.d0)), intent (in)       :: x(:)
+       real(psb_dpk_), intent (in)       :: x(:)
       type(psb_desc_type), intent (in)    :: desc_a
       integer, intent(out)                :: info
     end subroutine psb_dasumvs
     subroutine  psb_dmasum(res,x,desc_a,info)
       use psb_descriptor_type
-       real(kind(1.d0)), intent (out)      :: res(:)
+       real(psb_dpk_), intent (out)      :: res(:)
-       real(kind(1.d0)), intent (in)       :: x(:,:)
+       real(psb_dpk_), intent (in)       :: x(:,:)
       type(psb_desc_type), intent (in)    :: desc_a
       integer, intent(out)                :: info
     end subroutine psb_dmasum
     subroutine  psb_zasumvs(res,x,desc_a,info)
       use psb_descriptor_type
-       real(kind(1.d0)), intent (out)      :: res
+       real(psb_dpk_), intent (out)      :: res
-       complex(kind(1.d0)), intent (in)       :: x(:)
+       complex(psb_dpk_), intent (in)       :: x(:)
       type(psb_desc_type), intent (in)    :: desc_a
       integer, intent(out)                :: info
     end subroutine psb_zasumvs
     subroutine  psb_zmasum(res,x,desc_a,info)
       use psb_descriptor_type
-       real(kind(1.d0)), intent (out)      :: res(:)
+       real(psb_dpk_), intent (out)      :: res(:)
-       complex(kind(1.d0)), intent (in)       :: x(:,:)
+       complex(psb_dpk_), intent (in)       :: x(:,:)
       type(psb_desc_type), intent (in)    :: desc_a
       integer, intent(out)                :: info
     end subroutine psb_zmasum
@ -270,31 +270,31 @@ module psb_psblas_mod
  interface psb_genrm2
     function psb_dnrm2(x, desc_a, info, jx)
       use psb_descriptor_type
-       real(kind(1.d0))   psb_dnrm2
+       real(psb_dpk_)   psb_dnrm2
-       real(kind(1.d0)), intent (in)       :: x(:,:)
+       real(psb_dpk_), intent (in)       :: x(:,:)
       type(psb_desc_type), intent (in)    :: desc_a
       integer, optional, intent (in)      :: jx
       integer, intent(out)                :: info
     end function psb_dnrm2
     function psb_dnrm2v(x, desc_a, info)
       use psb_descriptor_type
-       real(kind(1.d0)) psb_dnrm2v
+       real(psb_dpk_) psb_dnrm2v
-       real(kind(1.d0)), intent (in)       :: x(:)
+       real(psb_dpk_), intent (in)       :: x(:)
       type(psb_desc_type), intent (in)    :: desc_a
       integer, intent(out)                :: info
     end function psb_dnrm2v
     function psb_znrm2(x, desc_a, info, jx)
       use psb_descriptor_type
-       real(kind(1.d0))   psb_znrm2
+       real(psb_dpk_)   psb_znrm2
-       complex(kind(1.d0)), intent (in)       :: x(:,:)
+       complex(psb_dpk_), intent (in)       :: x(:,:)
       type(psb_desc_type), intent (in)    :: desc_a
       integer, optional, intent (in)      :: jx
       integer, intent(out)                :: info
     end function psb_znrm2
     function psb_znrm2v(x, desc_a, info)
       use psb_descriptor_type
-       real(kind(1.d0)) psb_znrm2v
+       real(psb_dpk_) psb_znrm2v
-       complex(kind(1.d0)), intent (in)       :: x(:)
+       complex(psb_dpk_), intent (in)       :: x(:)
       type(psb_desc_type), intent (in)    :: desc_a
       integer, intent(out)                :: info
     end function psb_znrm2v
@ -303,15 +303,15 @@ module psb_psblas_mod
  interface psb_genrm2s
     subroutine  psb_dnrm2vs(res,x,desc_a,info)
       use psb_descriptor_type
-       real(kind(1.d0)), intent (out)      :: res
+       real(psb_dpk_), intent (out)      :: res
-       real(kind(1.d0)), intent (in)       :: x(:)
+       real(psb_dpk_), intent (in)       :: x(:)
       type(psb_desc_type), intent (in)    :: desc_a
       integer, intent(out)                :: info
     end subroutine psb_dnrm2vs
     subroutine  psb_znrm2vs(res,x,desc_a,info)
       use psb_descriptor_type
-       real(kind(1.d0)), intent (out)      :: res
+       real(psb_dpk_), intent (out)      :: res
-       complex(kind(1.d0)), intent (in)       :: x(:)
+       complex(psb_dpk_), intent (in)       :: x(:)
       type(psb_desc_type), intent (in)    :: desc_a
       integer, intent(out)                :: info
     end subroutine psb_znrm2vs
@ -322,7 +322,7 @@ module psb_psblas_mod
     function psb_dnrmi(a, desc_a,info)
       use psb_serial_mod
       use psb_descriptor_type
-       real(kind(1.d0))                    :: psb_dnrmi
+       real(psb_dpk_)                    :: psb_dnrmi
       type (psb_dspmat_type), intent (in) :: a
       type (psb_desc_type), intent (in)   :: desc_a
       integer, intent(out)                :: info
@ -330,7 +330,7 @@ module psb_psblas_mod
     function psb_znrmi(a, desc_a,info)
       use psb_serial_mod
       use psb_descriptor_type
-       real(kind(1.d0))                    :: psb_znrmi
+       real(psb_dpk_)                    :: psb_znrmi
       type (psb_zspmat_type), intent (in) :: a
       type (psb_desc_type), intent (in)   :: desc_a
       integer, intent(out)                :: info
@ -343,12 +343,12 @@ module psb_psblas_mod
       use psb_serial_mod
       use psb_descriptor_type
       type (psb_dspmat_type), intent(in)   :: a
-       real(kind(1.d0)), intent(inout)      :: x(:,:)
+       real(psb_dpk_), intent(inout)      :: x(:,:)
-       real(kind(1.d0)), intent(inout)      :: y(:,:)
+       real(psb_dpk_), intent(inout)      :: y(:,:)
-       real(kind(1.d0)), intent(in)         :: alpha, beta
+       real(psb_dpk_), intent(in)         :: alpha, beta
       type(psb_desc_type), intent(in)      :: desc_a
       character, optional, intent(in)      :: trans
-       real(kind(1.d0)), optional, intent(inout),target :: work(:)
+       real(psb_dpk_), optional, intent(inout),target :: work(:)
       integer, optional, intent(in)        :: k, jx, jy
       logical, optional, intent(in)        :: doswap
       integer, intent(out)                 :: info
@ -358,12 +358,12 @@ module psb_psblas_mod
       use psb_serial_mod
       use psb_descriptor_type
       type (psb_dspmat_type), intent(in)   :: a
-       real(kind(1.d0)), intent(inout)      :: x(:)
+       real(psb_dpk_), intent(inout)      :: x(:)
-       real(kind(1.d0)), intent(inout)      :: y(:)
+       real(psb_dpk_), intent(inout)      :: y(:)
-       real(kind(1.d0)), intent(in)         :: alpha, beta
+       real(psb_dpk_), intent(in)         :: alpha, beta
       type(psb_desc_type), intent(in)      :: desc_a
       character, optional, intent(in)      :: trans
-       real(kind(1.d0)), optional, intent(inout),target :: work(:)
+       real(psb_dpk_), optional, intent(inout),target :: work(:)
       logical, optional, intent(in)        :: doswap
       integer, intent(out)                 :: info
     end subroutine psb_dspmv
@ -372,12 +372,12 @@ module psb_psblas_mod
       use psb_serial_mod
       use psb_descriptor_type
       type (psb_zspmat_type), intent(in)   :: a
-       complex(kind(1.d0)), intent(inout)      :: x(:,:)
+       complex(psb_dpk_), intent(inout)      :: x(:,:)
-       complex(kind(1.d0)), intent(inout)      :: y(:,:)
+       complex(psb_dpk_), intent(inout)      :: y(:,:)
-       complex(kind(1.d0)), intent(in)         :: alpha, beta
+       complex(psb_dpk_), intent(in)         :: alpha, beta
       type(psb_desc_type), intent(in)      :: desc_a
       character, optional, intent(in)      :: trans
-       complex(kind(1.d0)), optional, intent(inout),target :: work(:)
+       complex(psb_dpk_), optional, intent(inout),target :: work(:)
       integer, optional, intent(in)        :: k, jx, jy
       logical, optional, intent(in)        :: doswap
       integer, intent(out)                 :: info
@ -387,12 +387,12 @@ module psb_psblas_mod
       use psb_serial_mod
       use psb_descriptor_type
       type (psb_zspmat_type), intent(in)   :: a
-       complex(kind(1.d0)), intent(inout)      :: x(:)
+       complex(psb_dpk_), intent(inout)      :: x(:)
-       complex(kind(1.d0)), intent(inout)      :: y(:)
+       complex(psb_dpk_), intent(inout)      :: y(:)
-       complex(kind(1.d0)), intent(in)         :: alpha, beta
+       complex(psb_dpk_), intent(in)         :: alpha, beta
       type(psb_desc_type), intent(in)      :: desc_a
       character, optional, intent(in)      :: trans
-       complex(kind(1.d0)), optional, intent(inout),target :: work(:)
+       complex(psb_dpk_), optional, intent(inout),target :: work(:)
       logical, optional, intent(in)        :: doswap
       integer, intent(out)                 :: info
     end subroutine psb_zspmv
@ -405,14 +405,14 @@ module psb_psblas_mod
       use psb_serial_mod
       use psb_descriptor_type
       type (psb_dspmat_type), intent(in)     :: t
-       real(kind(1.d0)), intent(in)           :: x(:,:)
+       real(psb_dpk_), intent(in)           :: x(:,:)
-       real(kind(1.d0)), intent(inout)        :: y(:,:)
+       real(psb_dpk_), intent(inout)        :: y(:,:)
-       real(kind(1.d0)), intent(in)           :: alpha, beta
+       real(psb_dpk_), intent(in)           :: alpha, beta
       type(psb_desc_type), intent(in)     :: desc_a
       character, optional, intent(in)        :: trans, unit
       integer, optional, intent(in)          :: n, jx, jy
       integer, optional, intent(in)          :: choice
-       real(kind(1.d0)), optional, intent(in),target :: work(:), diag(:)
+       real(psb_dpk_), optional, intent(in),target :: work(:), diag(:)
       integer, intent(out)               :: info
     end subroutine psb_dspsm
     subroutine psb_dspsv(alpha, t, x, beta, y,&
@ -421,13 +421,13 @@ module psb_psblas_mod
       use psb_serial_mod
       use psb_descriptor_type
       type (psb_dspmat_type), intent(in)     :: t
-       real(kind(1.d0)), intent(in)           :: x(:)
+       real(psb_dpk_), intent(in)           :: x(:)
-       real(kind(1.d0)), intent(inout)        :: y(:)
+       real(psb_dpk_), intent(inout)        :: y(:)
-       real(kind(1.d0)), intent(in)           :: alpha, beta
+       real(psb_dpk_), intent(in)           :: alpha, beta
       type(psb_desc_type), intent(in)        :: desc_a
       character, optional, intent(in)        :: trans, unit
       integer, optional, intent(in)          :: choice
-       real(kind(1.d0)), optional, intent(in),target :: work(:), diag(:)
+       real(psb_dpk_), optional, intent(in),target :: work(:), diag(:)
       integer, intent(out)                   :: info
     end subroutine psb_dspsv
     subroutine psb_zspsm(alpha, t, x, beta, y,&
@ -436,14 +436,14 @@ module psb_psblas_mod
       use psb_serial_mod
       use psb_descriptor_type
       type (psb_zspmat_type), intent(in)     :: t
-       complex(kind(1.d0)), intent(in)           :: x(:,:)
+       complex(psb_dpk_), intent(in)           :: x(:,:)
-       complex(kind(1.d0)), intent(inout)        :: y(:,:)
+       complex(psb_dpk_), intent(inout)        :: y(:,:)
-       complex(kind(1.d0)), intent(in)           :: alpha, beta
+       complex(psb_dpk_), intent(in)           :: alpha, beta
       type(psb_desc_type), intent(in)     :: desc_a
       character, optional, intent(in)        :: trans, unit
       integer, optional, intent(in)          :: n, jx, jy
       integer, optional, intent(in)          :: choice
-       complex(kind(1.d0)), optional, intent(in),target :: work(:), diag(:)
+       complex(psb_dpk_), optional, intent(in),target :: work(:), diag(:)
       integer, intent(out)               :: info
     end subroutine psb_zspsm
     subroutine psb_zspsv(alpha, t, x, beta, y,&
@ -452,13 +452,13 @@ module psb_psblas_mod
       use psb_serial_mod
       use psb_descriptor_type
       type (psb_zspmat_type), intent(in)     :: t
-       complex(kind(1.d0)), intent(in)           :: x(:)
+       complex(psb_dpk_), intent(in)           :: x(:)
-       complex(kind(1.d0)), intent(inout)        :: y(:)
+       complex(psb_dpk_), intent(inout)        :: y(:)
-       complex(kind(1.d0)), intent(in)           :: alpha, beta
+       complex(psb_dpk_), intent(in)           :: alpha, beta
       type(psb_desc_type), intent(in)        :: desc_a
       character, optional, intent(in)        :: trans, unit
       integer, optional, intent(in)          :: choice
-       complex(kind(1.d0)), optional, intent(in),target :: work(:), diag(:)
+       complex(psb_dpk_), optional, intent(in),target :: work(:), diag(:)
       integer, intent(out)                   :: info
     end subroutine psb_zspsv
  end interface
--- a/base/modules/psb_realloc_mod.F90
+++ b/base/modules/psb_realloc_mod.F90
@ -192,8 +192,8 @@ Contains
    use psb_error_mod
    ! ...Subroutine Arguments  
-    real(kind(1.d0)), allocatable, intent(in)  :: vin(:)
+    real(psb_dpk_), allocatable, intent(in)  :: vin(:)
-    real(kind(1.d0)), allocatable, intent(out) :: vout(:)
+    real(psb_dpk_), allocatable, intent(out) :: vout(:)
    integer         :: info
    ! ...Local Variables
@ -242,8 +242,8 @@ Contains
    use psb_error_mod
    ! ...Subroutine Arguments  
-    real(kind(1.d0)), allocatable, intent(in)  :: vin(:,:)
+    real(psb_dpk_), allocatable, intent(in)  :: vin(:,:)
-    real(kind(1.d0)), allocatable, intent(out) :: vout(:,:)
+    real(psb_dpk_), allocatable, intent(out) :: vout(:,:)
    integer         :: info
    ! ...Local Variables
@ -294,8 +294,8 @@ Contains
    use psb_error_mod
    ! ...Subroutine Arguments  
-    complex(kind(1.d0)), allocatable, intent(in)  :: vin(:)
+    complex(psb_dpk_), allocatable, intent(in)  :: vin(:)
-    complex(kind(1.d0)), allocatable, intent(out) :: vout(:)
+    complex(psb_dpk_), allocatable, intent(out) :: vout(:)
    integer         :: info
    ! ...Local Variables
@ -344,8 +344,8 @@ Contains
    use psb_error_mod
    ! ...Subroutine Arguments  
-    complex(kind(1.d0)), allocatable, intent(in)  :: vin(:,:)
+    complex(psb_dpk_), allocatable, intent(in)  :: vin(:,:)
-    complex(kind(1.d0)), allocatable, intent(out) :: vout(:,:)
+    complex(psb_dpk_), allocatable, intent(out) :: vout(:,:)
    integer         :: info
    ! ...Local Variables
@ -495,8 +495,8 @@ Contains
    use psb_error_mod
    ! ...Subroutine Arguments  
-    real(kind(1.d0)), intent(in)               :: vin(:)
+    real(psb_dpk_), intent(in)               :: vin(:)
-    real(kind(1.d0)), allocatable, intent(out) :: vout(:)
+    real(psb_dpk_), allocatable, intent(out) :: vout(:)
    integer         :: info
    ! ...Local Variables
@ -543,8 +543,8 @@ Contains
    use psb_error_mod
    ! ...Subroutine Arguments  
-    real(kind(1.d0)), intent(in)               :: vin(:,:)
+    real(psb_dpk_), intent(in)               :: vin(:,:)
-    real(kind(1.d0)), allocatable, intent(out) :: vout(:,:)
+    real(psb_dpk_), allocatable, intent(out) :: vout(:,:)
    integer         :: info
    ! ...Local Variables
@ -594,8 +594,8 @@ Contains
    use psb_error_mod
    ! ...Subroutine Arguments  
-    complex(kind(1.d0)), intent(in)               :: vin(:)
+    complex(psb_dpk_), intent(in)               :: vin(:)
-    complex(kind(1.d0)), allocatable, intent(out) :: vout(:)
+    complex(psb_dpk_), allocatable, intent(out) :: vout(:)
    integer         :: info
    ! ...Local Variables
@ -642,8 +642,8 @@ Contains
    use psb_error_mod
    ! ...Subroutine Arguments  
-    complex(kind(1.d0)), intent(in)               :: vin(:,:)
+    complex(psb_dpk_), intent(in)               :: vin(:,:)
-    complex(kind(1.d0)), allocatable, intent(out) :: vout(:,:)
+    complex(psb_dpk_), allocatable, intent(out) :: vout(:,:)
    integer         :: info
    ! ...Local Variables
@ -719,7 +719,7 @@ Contains
  function psb_dsize1d(vin)
    integer :: psb_dsize1d
-    real(kind(1.d0)), allocatable, intent(in) :: vin(:)
+    real(psb_dpk_), allocatable, intent(in) :: vin(:)
    if (.not.allocated(vin)) then 
      psb_dsize1d = 0
@ -729,7 +729,7 @@ Contains
  end function psb_dsize1d
  function psb_dsize2d(vin,dim)
    integer :: psb_dsize2d
-    real(kind(1.d0)), allocatable, intent(in) :: vin(:,:)
+    real(psb_dpk_), allocatable, intent(in) :: vin(:,:)
    integer, optional :: dim
    integer :: dim_
@ -749,7 +749,7 @@ Contains
  function psb_zsize1d(vin)
    integer :: psb_zsize1d
-    complex(kind(1.d0)), allocatable, intent(in) :: vin(:)
+    complex(psb_dpk_), allocatable, intent(in) :: vin(:)
    if (.not.allocated(vin)) then 
      psb_zsize1d = 0
@ -760,7 +760,7 @@ Contains
  function psb_zsize2d(vin,dim)
    integer :: psb_zsize2d
-    complex(kind(1.d0)), allocatable, intent(in) :: vin(:,:)
+    complex(psb_dpk_), allocatable, intent(in) :: vin(:,:)
    integer, optional :: dim
    integer :: dim_
@ -841,10 +841,10 @@ Contains
    ! ...Subroutine Arguments  
    Integer,Intent(in)                 :: len
-    real(kind(1.d0)),allocatable, intent(inout) :: v(:)
+    real(psb_dpk_),allocatable, intent(inout) :: v(:)
    integer         :: info
    integer, optional, intent(in)          :: addsz,newsz
-    real(kind(1.d0)), optional, intent(in) :: pad
+    real(psb_dpk_), optional, intent(in) :: pad
    ! ...Local Variables
    character(len=20)  :: name
    logical, parameter :: debug=.false.
@ -900,10 +900,10 @@ Contains
    ! ...Subroutine Arguments  
    Integer,Intent(in)                             :: len
-    complex(kind(1.d0)),allocatable, intent(inout) :: v(:)
+    complex(psb_dpk_),allocatable, intent(inout) :: v(:)
    integer                                        :: info
    integer, optional, intent(in)                  :: addsz,newsz
-    complex(kind(1.d0)), optional, intent(in)      :: pad
+    complex(psb_dpk_), optional, intent(in)      :: pad
    ! ...Local Variables
    character(len=20)  :: name
    logical, parameter :: debug=.false.
@ -1038,13 +1038,13 @@ Contains
    ! ...Subroutine Arguments  
    Integer,Intent(in) :: len
-    Real(kind(1.d0)),allocatable, intent(inout) :: rrax(:)
+    Real(psb_dpk_),allocatable, intent(inout) :: rrax(:)
    integer :: info
-    real(kind(1.d0)), optional, intent(in) :: pad
+    real(psb_dpk_), optional, intent(in) :: pad
    integer, optional, intent(in) :: lb
    ! ...Local Variables
-    Real(kind(1.d0)),allocatable  :: tmp(:)
+    Real(psb_dpk_),allocatable  :: tmp(:)
    Integer :: dim,err_act,err, lb_, lbi,ub_
    character(len=20)  :: name
    logical, parameter :: debug=.false.
@ -1061,7 +1061,7 @@ Contains
    endif
    if ((len<0)) then 
      err=4025
-      call psb_errpush(err,name,i_err=(/len,0,0,0,0/),a_err='real(kind(1.d0))')
+      call psb_errpush(err,name,i_err=(/len,0,0,0,0/),a_err='real(psb_dpk_)')
      goto 9999
    end if
    ub_ = lb_ + len-1
@ -1073,7 +1073,7 @@ Contains
        Allocate(tmp(lb_:ub_),stat=info)
        if (info /= 0) then
          err=4025
-          call psb_errpush(err,name,i_err=(/len,0,0,0,0/),a_err='real(kind(1.d0))')
+          call psb_errpush(err,name,i_err=(/len,0,0,0,0/),a_err='real(psb_dpk_)')
          goto 9999
        end if
        tmp(lb_:lb_-1+min(len,dim))=rrax(lbi:lbi-1+min(len,dim))
@ -1084,7 +1084,7 @@ Contains
      Allocate(rrax(lb_:ub_),stat=info)
      if (info /= 0) then
        err=4025
-        call psb_errpush(err,name,i_err=(/len,0,0,0,0/),a_err='real(kind(1.d0))')
+        call psb_errpush(err,name,i_err=(/len,0,0,0,0/),a_err='real(psb_dpk_)')
        goto 9999
      end if
    endif
@ -1113,13 +1113,13 @@ Contains
    ! ...Subroutine Arguments  
    Integer,Intent(in) :: len
-    complex(kind(1.d0)),allocatable, intent(inout):: rrax(:)
+    complex(psb_dpk_),allocatable, intent(inout):: rrax(:)
    integer :: info
-    complex(kind(1.d0)), optional, intent(in) :: pad
+    complex(psb_dpk_), optional, intent(in) :: pad
    integer, optional, intent(in) :: lb
    ! ...Local Variables
-    complex(kind(1.d0)),allocatable  :: tmp(:)
+    complex(psb_dpk_),allocatable  :: tmp(:)
    Integer :: dim,err_act,err,lb_,ub_,lbi
    character(len=20)  :: name
    logical, parameter :: debug=.false.
@ -1135,7 +1135,7 @@ Contains
    endif
    if ((len<0)) then 
      err=4025
-      call psb_errpush(err,name,i_err=(/len,0,0,0,0/),a_err='complex(kind(1.d0))')
+      call psb_errpush(err,name,i_err=(/len,0,0,0,0/),a_err='complex(psb_dpk_)')
      goto 9999
    end if
    ub_ = lb_+len-1
@ -1147,7 +1147,7 @@ Contains
        Allocate(tmp(lb_:ub_),stat=info)
        if (info /= 0) then
          err=4025
-          call psb_errpush(err,name,i_err=(/len,0,0,0,0/),a_err='complex(kind(1.d0))')
+          call psb_errpush(err,name,i_err=(/len,0,0,0,0/),a_err='complex(psb_dpk_)')
          goto 9999
        end if
        tmp(lb_:lb_-1+min(len,dim))=rrax(lbi:lbi-1+min(len,dim))
@ -1158,7 +1158,7 @@ Contains
      Allocate(rrax(lb_:ub_),stat=info)
      if (info /= 0) then
        err=4025
-        call psb_errpush(err,name,i_err=(/len,0,0,0,0/),a_err='complex(kind(1.d0))')
+        call psb_errpush(err,name,i_err=(/len,0,0,0,0/),a_err='complex(psb_dpk_)')
        goto 9999
      end if
    endif
@ -1187,14 +1187,14 @@ Contains
    use psb_error_mod
    ! ...Subroutine Arguments  
    Integer,Intent(in) :: len1,len2
-    Real(kind(1.d0)),allocatable :: rrax(:,:)
+    Real(psb_dpk_),allocatable :: rrax(:,:)
    integer :: info
-    real(kind(1.d0)), optional, intent(in) :: pad
+    real(psb_dpk_), optional, intent(in) :: pad
    Integer,Intent(in), optional  :: lb1,lb2
    ! ...Local Variables
-    Real(kind(1.d0)),allocatable  :: tmp(:,:)
+    Real(psb_dpk_),allocatable  :: tmp(:,:)
    Integer :: dim,err_act,err, dim2,lb1_, lb2_, ub1_, ub2_,&
         & lbi1, lbi2
    character(len=20)  :: name
@ -1217,12 +1217,12 @@ Contains
    if (len1 < 0) then
      err=4025 
-      call psb_errpush(err,name,i_err=(/len1,0,0,0,0/),a_err='real(kind(1.d0))')
+      call psb_errpush(err,name,i_err=(/len1,0,0,0,0/),a_err='real(psb_dpk_)')
      goto 9999
    end if
    if (len2 < 0) then
      err=4025
-      call psb_errpush(err,name,i_err=(/len2,0,0,0,0/),a_err='real(kind(1.d0))')
+      call psb_errpush(err,name,i_err=(/len2,0,0,0,0/),a_err='real(psb_dpk_)')
      goto 9999
    end if
@ -1237,7 +1237,7 @@ Contains
        Allocate(tmp(lb1_:ub1_,lb2_:ub2_),stat=info)
        if (info /= 0) then
          err=4025
-          call psb_errpush(err,name,i_err=(/len1*len2,0,0,0,0/),a_err='real(kind(1.d0))')
+          call psb_errpush(err,name,i_err=(/len1*len2,0,0,0,0/),a_err='real(psb_dpk_)')
          goto 9999
        end if
        tmp(lb1_:lb1_-1+min(len1,dim),lb2_:lb2_-1+min(len2,dim2)) = &
@ -1250,7 +1250,7 @@ Contains
      Allocate(rrax(lb1_:ub1_,lb2_:ub2_),stat=info)
      if (info /= 0) then
        err=4025
-        call psb_errpush(err,name,i_err=(/len1*len2,0,0,0,0/),a_err='real(kind(1.d0))')
+        call psb_errpush(err,name,i_err=(/len1*len2,0,0,0,0/),a_err='real(psb_dpk_)')
        goto 9999
      end if
    endif
@ -1280,14 +1280,14 @@ Contains
    use psb_error_mod
    ! ...Subroutine Arguments  
    Integer,Intent(in) :: len1,len2
-    complex(kind(1.d0)),allocatable :: rrax(:,:)
+    complex(psb_dpk_),allocatable :: rrax(:,:)
    integer :: info
-    complex(kind(1.d0)), optional, intent(in) :: pad
+    complex(psb_dpk_), optional, intent(in) :: pad
    Integer,Intent(in), optional  :: lb1,lb2
    ! ...Local Variables
-    complex(kind(1.d0)),allocatable  :: tmp(:,:)
+    complex(psb_dpk_),allocatable  :: tmp(:,:)
    Integer :: dim,err_act,err,dim2,lb1_, lb2_, ub1_, ub2_,&
         & lbi1, lbi2
    character(len=20)  :: name
@ -1310,12 +1310,12 @@ Contains
    if (len1 < 0) then
      err=4025
-      call psb_errpush(err,name,i_err=(/len1,0,0,0,0/),a_err='complex(kind(1.d0))')
+      call psb_errpush(err,name,i_err=(/len1,0,0,0,0/),a_err='complex(psb_dpk_)')
      goto 9999
    end if
    if (len2 < 0) then
      err=4025
-      call psb_errpush(err,name,i_err=(/len2,0,0,0,0/),a_err='complex(kind(1.d0))')
+      call psb_errpush(err,name,i_err=(/len2,0,0,0,0/),a_err='complex(psb_dpk_)')
      goto 9999
    end if
@ -1330,7 +1330,7 @@ Contains
        Allocate(tmp(lb1_:ub1_,lb2_:ub2_),stat=info)
        if (info /= 0) then
          err=4025
-          call psb_errpush(err,name,i_err=(/len1*len2,0,0,0,0/),a_err='complex(kind(1.d0))')
+          call psb_errpush(err,name,i_err=(/len1*len2,0,0,0,0/),a_err='complex(psb_dpk_)')
          goto 9999
        end if
        tmp(lb1_:lb1_-1+min(len1,dim),lb2_:lb2_-1+min(len2,dim2)) = &
@ -1343,7 +1343,7 @@ Contains
      Allocate(rrax(lb1_:ub1_,lb2_:ub2_),stat=info)
      if (info /= 0) then
        err=4025
-        call psb_errpush(err,name,i_err=(/len1*len2,0,0,0,0/),a_err='complex(kind(1.d0))')
+        call psb_errpush(err,name,i_err=(/len1*len2,0,0,0,0/),a_err='complex(psb_dpk_)')
        goto 9999
      end if
    endif
@ -1513,7 +1513,7 @@ Contains
    ! ...Subroutine Arguments  
    Integer,Intent(in) :: len
    Integer,allocatable, intent(inout)  :: rrax(:),y(:)
-    Real(Kind(1.d0)),allocatable, intent(inout) :: z(:)
+    Real(psb_dpk_),allocatable, intent(inout) :: z(:)
    integer :: info
    character(len=20)  :: name
    integer :: err_act, err
@ -1562,7 +1562,7 @@ Contains
    ! ...Subroutine Arguments  
    Integer,Intent(in) :: len
    Integer,allocatable, intent(inout) :: rrax(:),y(:)
-    complex(Kind(1.d0)),allocatable, intent(inout) :: z(:)
+    complex(psb_dpk_),allocatable, intent(inout) :: z(:)
    integer :: info
    character(len=20)  :: name
    integer :: err_act, err
@ -1606,7 +1606,7 @@ Contains
  Subroutine psb_dtransfer1d(vin,vout,info)
    use psb_error_mod
-    real(kind(1.d0)), allocatable, intent(inout) :: vin(:),vout(:)
+    real(psb_dpk_), allocatable, intent(inout) :: vin(:),vout(:)
    integer, intent(out) :: info 
    !
    ! 
@ -1646,7 +1646,7 @@ Contains
  Subroutine psb_dtransfer2d(vin,vout,info)
    use psb_error_mod
-    real(kind(1.d0)), allocatable, intent(inout) :: vin(:,:),vout(:,:)
+    real(psb_dpk_), allocatable, intent(inout) :: vin(:,:),vout(:,:)
    integer, intent(out) :: info 
    !
    ! 
@ -1685,7 +1685,7 @@ Contains
  Subroutine psb_ztransfer1d(vin,vout,info)
    use psb_error_mod
-    complex(kind(1.d0)), allocatable, intent(inout) :: vin(:),vout(:)
+    complex(psb_dpk_), allocatable, intent(inout) :: vin(:),vout(:)
    integer, intent(out) :: info 
    !
    ! 
@ -1721,7 +1721,7 @@ Contains
  Subroutine psb_ztransfer2d(vin,vout,info)
    use psb_error_mod
-    complex(kind(1.d0)), allocatable, intent(inout) :: vin(:,:),vout(:,:)
+    complex(psb_dpk_), allocatable, intent(inout) :: vin(:,:),vout(:,:)
    integer, intent(out) :: info 
    !
    ! 
--- a/base/modules/psb_serial_mod.f90
+++ b/base/modules/psb_serial_mod.f90
@ -30,6 +30,7 @@
 !!$ 
 !!$  
 module psb_serial_mod
  use psb_const_mod
  use psb_spmat_type
  use psb_string_mod
  use psb_sort_mod
@ -62,14 +63,14 @@ module psb_serial_mod
    subroutine psb_dcsrws(rw,a,info,trans)
      use psb_spmat_type
      type(psb_dspmat_type) :: a
-      real(kind(1.d0)), allocatable   :: rw(:) 
+      real(psb_dpk_), allocatable   :: rw(:) 
      integer :: info
      character, optional :: trans
    end subroutine psb_dcsrws
    subroutine psb_zcsrws(rw,a,info,trans)
      use psb_spmat_type
      type(psb_zspmat_type) :: a
-      complex(kind(1.d0)), allocatable :: rw(:) 
+      complex(psb_dpk_), allocatable :: rw(:) 
      integer :: info
      character, optional :: trans
    end subroutine psb_zcsrws
@ -81,34 +82,34 @@ module psb_serial_mod
    subroutine psb_dcssm(alpha,t,b,beta,c,info,trans,unitd,d)
      use psb_spmat_type
      type(psb_dspmat_type) :: t
-      real(kind(1.d0)) :: alpha, beta, b(:,:), c(:,:)
+      real(psb_dpk_) :: alpha, beta, b(:,:), c(:,:)
      integer :: info
      character, optional :: trans, unitd
-      real(kind(1.d0)), optional, target :: d(:)
+      real(psb_dpk_), optional, target :: d(:)
    end subroutine psb_dcssm
    subroutine psb_dcssv(alpha,t,b,beta,c,info,trans,unitd,d)
      use psb_spmat_type
      type(psb_dspmat_type) :: t
-      real(kind(1.d0)) :: alpha, beta, b(:), c(:)
+      real(psb_dpk_) :: alpha, beta, b(:), c(:)
      integer :: info
      character, optional :: trans, unitd
-      real(kind(1.d0)), optional, target :: d(:)
+      real(psb_dpk_), optional, target :: d(:)
    end subroutine psb_dcssv
    subroutine psb_zcssm(alpha,t,b,beta,c,info,trans,unitd,d)
      use psb_spmat_type
      type(psb_zspmat_type) :: t
-      complex(kind(1.d0)) :: alpha, beta, b(:,:), c(:,:)
+      complex(psb_dpk_) :: alpha, beta, b(:,:), c(:,:)
      integer :: info
      character, optional :: trans, unitd
-      complex(kind(1.d0)), optional, target :: d(:)
+      complex(psb_dpk_), optional, target :: d(:)
    end subroutine psb_zcssm
    subroutine psb_zcssv(alpha,t,b,beta,c,info,trans,unitd,d)
      use psb_spmat_type
      type(psb_zspmat_type) :: t
-      complex(kind(1.d0)) :: alpha, beta, b(:), c(:)
+      complex(psb_dpk_) :: alpha, beta, b(:), c(:)
      integer :: info
      character, optional :: trans, unitd
-      complex(kind(1.d0)), optional, target :: d(:)
+      complex(psb_dpk_), optional, target :: d(:)
    end subroutine psb_zcssv
  end interface
@ -116,28 +117,28 @@ module psb_serial_mod
    subroutine psb_dcsmv(alpha,a,b,beta,c,info,trans)
      use psb_spmat_type
      type(psb_dspmat_type) :: a
-      real(kind(1.d0)) :: alpha, beta, b(:), c(:)
+      real(psb_dpk_) :: alpha, beta, b(:), c(:)
      integer :: info
      character, optional :: trans
    end subroutine psb_dcsmv
    subroutine psb_dcsmm(alpha,a,b,beta,c,info,trans)
      use psb_spmat_type
      type(psb_dspmat_type) :: a
-      real(kind(1.d0)) :: alpha, beta, b(:,:), c(:,:)
+      real(psb_dpk_) :: alpha, beta, b(:,:), c(:,:)
      integer :: info
      character, optional :: trans
    end subroutine psb_dcsmm
    subroutine psb_zcsmv(alpha,a,b,beta,c,info,trans)
      use psb_spmat_type
      type(psb_zspmat_type) :: a
-      complex(kind(1.d0)) :: alpha, beta, b(:), c(:)
+      complex(psb_dpk_) :: alpha, beta, b(:), c(:)
      integer :: info
      character, optional :: trans
    end subroutine psb_zcsmv
    subroutine psb_zcsmm(alpha,a,b,beta,c,info,trans)
      use psb_spmat_type
      type(psb_zspmat_type) :: a
-      complex(kind(1.d0)) :: alpha, beta, b(:,:), c(:,:)
+      complex(psb_dpk_) :: alpha, beta, b(:,:), c(:,:)
      integer :: info
      character, optional :: trans
    end subroutine psb_zcsmm
@ -291,7 +292,7 @@ module psb_serial_mod
      use psb_spmat_type
      integer, intent(in) :: nz, imin,imax,jmin,jmax
      integer, intent(in) :: ia(:),ja(:)
-      real(kind(1.d0)), intent(in) :: val(:)
+      real(psb_dpk_), intent(in) :: val(:)
      type(psb_dspmat_type), intent(inout) :: a
      integer, intent(out) :: info
      integer, intent(in), optional :: gtl(:)
@ -301,7 +302,7 @@ module psb_serial_mod
      use psb_spmat_type
      integer, intent(in) :: nz, imin,imax,jmin,jmax
      integer, intent(in) :: ia(:),ja(:)
-      complex(kind(1.d0)), intent(in) :: val(:)
+      complex(psb_dpk_), intent(in) :: val(:)
      type(psb_zspmat_type), intent(inout) :: a
      integer, intent(out) :: info
      integer, intent(in), optional :: gtl(:)
@ -378,13 +379,13 @@ module psb_serial_mod
  end interface
  interface psb_csnmi
-    real(kind(1.d0)) function psb_dcsnmi(a,info,trans)
+    real(psb_dpk_) function psb_dcsnmi(a,info,trans)
      use psb_spmat_type
      type(psb_dspmat_type), intent(in)  :: a
      integer, intent(out)       :: info
      character, optional        :: trans
    end function psb_dcsnmi
-    real(kind(1.d0)) function psb_zcsnmi(a,info,trans)
+    real(psb_dpk_) function psb_zcsnmi(a,info,trans)
      use psb_spmat_type
      type(psb_zspmat_type), intent(in)  :: a
      integer, intent(out)       :: info
@ -417,13 +418,13 @@ module psb_serial_mod
     subroutine psb_dspgtdiag(a,d,info)
       use psb_spmat_type
       type(psb_dspmat_type), intent(in)     :: a
-       real(kind(1.d0)), intent(inout) :: d(:) 
+       real(psb_dpk_), intent(inout) :: d(:) 
       integer, intent(out)  :: info
     end subroutine psb_dspgtdiag
     subroutine psb_zspgtdiag(a,d,info)
       use psb_spmat_type
       type(psb_zspmat_type), intent(in)     :: a
-       complex(kind(1.d0)), intent(inout) :: d(:) 
+       complex(psb_dpk_), intent(inout) :: d(:) 
       integer, intent(out)  :: info
     end subroutine psb_zspgtdiag
  end interface
@ -432,13 +433,13 @@ module psb_serial_mod
     subroutine psb_dspscal(a,d,info)
       use psb_spmat_type
       type(psb_dspmat_type), intent(inout) :: a
-       real(kind(1.d0)), intent(in) :: d(:) 
+       real(psb_dpk_), intent(in) :: d(:) 
       integer, intent(out)  :: info
     end subroutine psb_dspscal
     subroutine psb_zspscal(a,d,info)
       use psb_spmat_type
       type(psb_zspmat_type), intent(inout) :: a
-       complex(kind(1.d0)), intent(in) :: d(:) 
+       complex(psb_dpk_), intent(in) :: d(:) 
       integer, intent(out)  :: info
     end subroutine psb_zspscal
  end interface
@ -478,7 +479,7 @@ module psb_serial_mod
       integer, intent(in)                  :: irw
       integer, intent(out)                 :: nz
       integer, allocatable, intent(inout)  :: ia(:), ja(:)
-       real(kind(1.d0)), allocatable,  intent(inout)    :: val(:)
+       real(psb_dpk_), allocatable,  intent(inout)    :: val(:)
       integer,intent(out)                  :: info
       logical, intent(in), optional        :: append
       integer, intent(in), optional        :: iren(:)
@ -493,7 +494,7 @@ module psb_serial_mod
       integer, intent(in)                  :: irw
       integer, intent(out)                 :: nz
       integer, allocatable, intent(inout)  :: ia(:), ja(:)
-       complex(kind(1.d0)), allocatable,  intent(inout)    :: val(:)
+       complex(psb_dpk_), allocatable,  intent(inout)    :: val(:)
       integer,intent(out)                  :: info
       logical, intent(in), optional        :: append
       integer, intent(in), optional        :: iren(:)
@ -522,28 +523,32 @@ module psb_serial_mod
  interface psb_gelp
    ! 2-D version
    subroutine psb_dgelp(trans,iperm,x,info)
-      real(kind(1.d0)), intent(inout)      ::  x(:,:)
+      use psb_const_mod
      real(psb_dpk_), intent(inout)      ::  x(:,:)
      integer, intent(in)                  ::  iperm(:)
      integer, intent(out)                 ::  info
      character, intent(in)                :: trans
    end subroutine psb_dgelp
    ! 1-D version
    subroutine psb_dgelpv(trans,iperm,x,info)
-      real(kind(1.d0)), intent(inout)    ::  x(:)
+      use psb_const_mod
      real(psb_dpk_), intent(inout)    ::  x(:)
      integer, intent(in)                  ::  iperm(:)
      integer, intent(out)                 ::  info
      character, intent(in)              :: trans
    end subroutine psb_dgelpv
    ! 2-D version
    subroutine psb_zgelp(trans,iperm,x,info)
-      complex(kind(1.d0)), intent(inout)      ::  x(:,:)
+      use psb_const_mod
      complex(psb_dpk_), intent(inout)      ::  x(:,:)
      integer, intent(in)                  ::  iperm(:)
      integer, intent(out)                 ::  info
      character, intent(in)                :: trans
    end subroutine psb_zgelp
    ! 1-D version
    subroutine psb_zgelpv(trans,iperm,x,info)
-      complex(kind(1.d0)), intent(inout)    ::  x(:)
+      use psb_const_mod
      complex(psb_dpk_), intent(inout)    ::  x(:)
      integer, intent(in)                  ::  iperm(:)
      integer, intent(out)                 ::  info
      character, intent(in)              :: trans
--- a/base/modules/psb_sort_mod.f90
+++ b/base/modules/psb_sort_mod.f90
@ -42,7 +42,7 @@
 !  Addison-Wesley
 !
 module psb_sort_mod
-
+  use psb_const_mod
  ! 
  !  The up/down constant are defined in pairs having 
@ -66,12 +66,12 @@ module psb_sort_mod
  end type psb_int_idx_heap
  type psb_double_idx_heap
    integer              :: last, dir
-    real(kind(1.d0)), allocatable :: keys(:)
+    real(psb_dpk_), allocatable :: keys(:)
    integer, allocatable          :: idxs(:)
  end type psb_double_idx_heap
  type psb_dcomplex_idx_heap
    integer              :: last, dir
-    complex(kind(1.d0)), allocatable :: keys(:)
+    complex(psb_dpk_), allocatable :: keys(:)
    integer, allocatable             :: idxs(:)
  end type psb_dcomplex_idx_heap
@ -182,7 +182,7 @@ contains
  subroutine dmsort(x,ix,dir,flag)
    use psb_error_mod
    implicit none 
-    real(kind(1.d0)), intent(inout)  :: x(:) 
+    real(psb_dpk_), intent(inout)  :: x(:) 
    integer, optional, intent(in)    :: dir, flag
    integer, optional, intent(inout) :: ix(:)
@ -241,7 +241,7 @@ contains
  subroutine zamsort(x,ix,dir,flag)
    use psb_error_mod
    implicit none 
-    complex(kind(1.d0)), intent(inout)  :: x(:) 
+    complex(psb_dpk_), intent(inout)  :: x(:) 
    integer, optional, intent(in)    :: dir, flag
    integer, optional, intent(inout) :: ix(:)
@ -417,7 +417,7 @@ contains
  subroutine dqsort(x,ix,dir,flag)
    use psb_error_mod
    implicit none 
-    real(kind(1.d0)), intent(inout)  :: x(:) 
+    real(psb_dpk_), intent(inout)  :: x(:) 
    integer, optional, intent(in)    :: dir, flag
    integer, optional, intent(inout) :: ix(:)
@ -493,7 +493,7 @@ contains
  subroutine zqsort(x,ix,dir,flag)
    use psb_error_mod
    implicit none 
-    complex(kind(1.d0)), intent(inout)  :: x(:) 
+    complex(psb_dpk_), intent(inout)  :: x(:) 
    integer, optional, intent(in)    :: dir, flag
    integer, optional, intent(inout) :: ix(:)
@ -692,12 +692,12 @@ contains
  subroutine dhsort(x,ix,dir,flag)
    use psb_error_mod
    implicit none 
-    real(kind(1.d0)), intent(inout)  :: x(:) 
+    real(psb_dpk_), intent(inout)  :: x(:) 
    integer, optional, intent(in)    :: dir, flag
    integer, optional, intent(inout) :: ix(:)
    integer          :: dir_, flag_, n, i, l, err_act,info
-    real(kind(1.d0)) :: key
+    real(psb_dpk_) :: key
    integer          :: index
    character(len=20)  :: name
@ -800,12 +800,12 @@ contains
  subroutine zhsort(x,ix,dir,flag)
    use psb_error_mod
    implicit none 
-    complex(kind(1.d0)), intent(inout) :: x(:) 
+    complex(psb_dpk_), intent(inout) :: x(:) 
    integer, optional, intent(in)      :: dir, flag
    integer, optional, intent(inout)   :: ix(:)
    integer             :: dir_, flag_, n, i, l, err_act,info
-    complex(kind(1.d0)) :: key
+    complex(psb_dpk_) :: key
    integer             :: index
    character(len=20)  :: name
@ -1074,7 +1074,7 @@ contains
    use psb_realloc_mod
    implicit none 
-    real(kind(1.d0)), intent(in)      :: key
+    real(psb_dpk_), intent(in)      :: key
    integer, intent(in)               :: index
    type(psb_double_idx_heap), intent(inout) :: heap
    integer, intent(out)              :: info
@ -1106,7 +1106,7 @@ contains
    type(psb_double_idx_heap), intent(inout) :: heap
    integer, intent(out)              :: index,info
-    real(kind(1.d0)), intent(out)     :: key
+    real(psb_dpk_), intent(out)     :: key
    info = 0
@ -1295,7 +1295,7 @@ contains
    use psb_realloc_mod
    implicit none 
-    complex(kind(1.d0)), intent(in)            :: key
+    complex(psb_dpk_), intent(in)            :: key
    integer, intent(in)                        :: index
    type(psb_dcomplex_idx_heap), intent(inout) :: heap
    integer, intent(out)                       :: info
@ -1326,7 +1326,7 @@ contains
    type(psb_dcomplex_idx_heap), intent(inout) :: heap
    integer, intent(out)                       :: index,info
-    complex(kind(1.d0)), intent(out)           :: key
+    complex(psb_dpk_), intent(out)           :: key
    info = 0
@ -1582,13 +1582,13 @@ contains
    !   heap: the heap
    !   dir:  sorting direction
-    real(kind(1.d0)), intent(in)    :: key
+    real(psb_dpk_), intent(in)    :: key
    integer, intent(in)             :: dir
-    real(kind(1.d0)), intent(inout) :: heap(:)
+    real(psb_dpk_), intent(inout) :: heap(:)
    integer, intent(inout)          :: last
    integer, intent(out)            :: info
    integer                         :: i, i2
-    real(kind(1.d0))                :: temp
+    real(psb_dpk_)                :: temp
    info = 0
    if (last < 0) then 
@ -1680,14 +1680,14 @@ contains
  subroutine psi_double_heap_get_first(key,last,heap,dir,info)
    implicit none 
-    real(kind(1.d0)), intent(inout) :: key
+    real(psb_dpk_), intent(inout) :: key
    integer, intent(inout)          :: last
    integer, intent(in)             :: dir
-    real(kind(1.d0)), intent(inout) :: heap(:)
+    real(psb_dpk_), intent(inout) :: heap(:)
    integer, intent(out)            :: info
    integer                 :: i, j
-    real(kind(1.d0))        :: temp
+    real(psb_dpk_)        :: temp
    info = 0
@ -1812,13 +1812,13 @@ contains
    !   heap: the heap
    !   dir:  sorting direction
-    complex(kind(1.d0)), intent(in)    :: key
+    complex(psb_dpk_), intent(in)    :: key
    integer, intent(in)                :: dir
-    complex(kind(1.d0)), intent(inout) :: heap(:)
+    complex(psb_dpk_), intent(inout) :: heap(:)
    integer, intent(inout)             :: last
    integer, intent(out)               :: info
    integer                            :: i, i2
-    complex(kind(1.d0))                :: temp
+    complex(psb_dpk_)                :: temp
    info = 0
    if (last < 0) then 
@ -1910,14 +1910,14 @@ contains
  subroutine psi_dcomplex_heap_get_first(key,last,heap,dir,info)
    implicit none 
-    complex(kind(1.d0)), intent(inout) :: key
+    complex(psb_dpk_), intent(inout) :: key
    integer, intent(inout)             :: last
    integer, intent(in)                :: dir
-    complex(kind(1.d0)), intent(inout) :: heap(:)
+    complex(psb_dpk_), intent(inout) :: heap(:)
    integer, intent(out)               :: info
    integer                    :: i, j
-    complex(kind(1.d0))        :: temp
+    complex(psb_dpk_)        :: temp
    info = 0
@ -2301,13 +2301,13 @@ contains
    !   idxs: the indices
    !   dir:  sorting direction
-    real(kind(1.d0)), intent(in)     :: key
+    real(psb_dpk_), intent(in)     :: key
    integer, intent(in)              :: index,dir
-    real(kind(1.d0)), intent(inout)  :: heap(:)
+    real(psb_dpk_), intent(inout)  :: heap(:)
    integer, intent(inout)           :: idxs(:),last
    integer, intent(out)             :: info
    integer                          :: i, i2, itemp
-    real(kind(1.d0))                 :: temp 
+    real(psb_dpk_)                 :: temp 
    info = 0
    if (last < 0) then 
@ -2413,14 +2413,14 @@ contains
  subroutine psi_double_idx_heap_get_first(key,index,last,heap,idxs,dir,info)
    implicit none 
-    real(kind(1.d0)), intent(inout) :: heap(:)
+    real(psb_dpk_), intent(inout) :: heap(:)
    integer, intent(out)            :: index,info
    integer, intent(inout)          :: last,idxs(:)
    integer, intent(in)             :: dir
-    real(kind(1.d0)), intent(out)   :: key
+    real(psb_dpk_), intent(out)   :: key
    integer                         :: i, j,itemp
-    real(kind(1.d0))                :: temp
+    real(psb_dpk_)                :: temp
    info = 0
    if (last <= 0) then 
@ -2559,13 +2559,13 @@ contains
    !   idxs: the indices
    !   dir:  sorting direction
-    complex(kind(1.d0)), intent(in)    :: key
+    complex(psb_dpk_), intent(in)    :: key
    integer, intent(in)                :: index,dir
-    complex(kind(1.d0)), intent(inout) :: heap(:)
+    complex(psb_dpk_), intent(inout) :: heap(:)
    integer, intent(inout)             :: idxs(:),last
    integer, intent(out)               :: info
    integer                            :: i, i2, itemp
-    complex(kind(1.d0))                :: temp 
+    complex(psb_dpk_)                :: temp 
    info = 0
    if (last < 0) then 
@ -2671,14 +2671,14 @@ contains
  subroutine psi_dcomplex_idx_heap_get_first(key,index,last,heap,idxs,dir,info)
    implicit none 
-    complex(kind(1.d0)), intent(inout) :: heap(:)
+    complex(psb_dpk_), intent(inout) :: heap(:)
    integer, intent(out)               :: index,info
    integer, intent(inout)             :: last,idxs(:)
    integer, intent(in)                :: dir
-    complex(kind(1.d0)), intent(out)   :: key
+    complex(psb_dpk_), intent(out)   :: key
    integer                            :: i, j, itemp
-    complex(kind(1.d0))                :: temp
+    complex(psb_dpk_)                :: temp
    info = 0
    if (last <= 0) then 
--- a/base/modules/psb_spmat_type.f90
+++ b/base/modules/psb_spmat_type.f90
@ -61,7 +61,7 @@ module psb_spmat_type
  !|    character(len=5) :: fida         ! Storage format: CSR,COO etc.
  !|    character(len=11) :: descra      ! Matrix type: encodes general, triang. 
  !|    integer     :: infoa(psb_ifasize_)  ! Additional integer info 
-  !|    real(kind(1.d0)), allocatable  :: aspk(:)  ! Coefficients
+  !|    real(psb_dpk_), allocatable  :: aspk(:)  ! Coefficients
  !|    integer, allocatable :: ia1(:), ia2(:)     ! Row/column indices encoded
  !|    integer, allocatable :: pl(:), pr(:)       ! Row/column permutation
  !|  end type psb_dspmat_type
@ -70,7 +70,7 @@ module psb_spmat_type
  !|    character(len=5) :: fida
  !|    character(len=11) :: descra
  !|    integer     :: infoa(psb_ifasize_)
-  !|    complex(kind(1.d0)), allocatable  :: aspk(:)
+  !|    complex(psb_dpk_), allocatable  :: aspk(:)
  !|    integer, allocatable  :: ia1(:), ia2(:)
  !|    integer, allocatable  :: pl(:), pr(:)
  !|  end type psb_zspmat_type
@ -126,7 +126,7 @@ module psb_spmat_type
    character(len=5) :: fida
    character(len=11) :: descra
    integer     :: infoa(psb_ifasize_)
-    real(kind(1.d0)), allocatable  :: aspk(:)
+    real(psb_dpk_), allocatable  :: aspk(:)
    integer, allocatable :: ia1(:), ia2(:)
    integer, allocatable :: pl(:), pr(:)
  end type psb_dspmat_type
@ -135,7 +135,7 @@ module psb_spmat_type
    character(len=5) :: fida
    character(len=11) :: descra
    integer     :: infoa(psb_ifasize_)
-    complex(kind(1.d0)), allocatable  :: aspk(:)
+    complex(psb_dpk_), allocatable  :: aspk(:)
    integer, allocatable  :: ia1(:), ia2(:)
    integer, allocatable  :: pl(:), pr(:)
  end type psb_zspmat_type
@ -809,23 +809,23 @@ contains
    !locals
    integer :: val
-    val   = 4*size(a%infoa)
+    val   = psb_sizeof_int*size(a%infoa)
    if (allocated(a%aspk)) then 
-      val = val + 8 * size(a%aspk)
+      val = val + psb_sizeof_dp  * size(a%aspk)
    endif
    if (allocated(a%ia1)) then 
-      val = val + 4 * size(a%ia1)
+      val = val + psb_sizeof_int * size(a%ia1)
    endif
    if (allocated(a%ia2)) then 
-      val = val + 4 * size(a%ia2)
+      val = val + psb_sizeof_int * size(a%ia2)
    endif
    if (allocated(a%pl)) then 
-      val = val + 4 * size(a%pl)
+      val = val + psb_sizeof_int * size(a%pl)
    endif
    if (allocated(a%pr)) then 
-      val = val + 4 * size(a%pr)
+      val = val + psb_sizeof_int * size(a%pr)
    endif
@ -1296,23 +1296,23 @@ contains
    !locals
    integer :: val
-    val   = 4*size(a%infoa)
+    val   = psb_sizeof_int*size(a%infoa)
    if (allocated(a%aspk)) then 
-      val = val + 16 * size(a%aspk)
+      val = val + 2 * psb_sizeof_dp * size(a%aspk)
    endif
    if (allocated(a%ia1)) then 
-      val = val + 4 * size(a%ia1)
+      val = val + psb_sizeof_int * size(a%ia1)
    endif
    if (allocated(a%ia2)) then 
-      val = val + 4 * size(a%ia2)
+      val = val + psb_sizeof_int * size(a%ia2)
    endif
    if (allocated(a%pl)) then 
-      val = val + 4 * size(a%pl)
+      val = val + psb_sizeof_int * size(a%pl)
    endif
    if (allocated(a%pr)) then 
-      val = val + 4 * size(a%pr)
+      val = val + psb_sizeof_int * size(a%pr)
    endif
--- a/base/modules/psb_spsb_mod.f90
+++ b/base/modules/psb_spsb_mod.f90
@ -32,143 +32,155 @@
 module psb_spsb_mod
  interface
-      real(kind(1.d0)) function ddot(n,dx,incx,dy,incy)
+    real(psb_dpk_) function ddot(n,dx,incx,dy,incy)
-      real(kind(1.d0)) :: dx(*),dy(*),dtemp
+      use psb_const_mod
      real(psb_dpk_) :: dx(*),dy(*),dtemp
      integer          :: i,incx,incy,ix,iy,m,mp1,n
    end function ddot
- end interface
+  end interface
  interface
-      subroutine  daxpby(m, n, alpha, X, lldx, beta, Y, lldy, info)
+    subroutine  daxpby(m, n, alpha, X, lldx, beta, Y, lldy, info)
-      real(kind(1.d0)) :: X(lldx,*), Y(lldy,*)
+      use psb_const_mod
-      real(kind(1.d0)) :: alpha, beta
+      real(psb_dpk_) :: X(lldx,*), Y(lldy,*)
      real(psb_dpk_) :: alpha, beta
      integer          :: n, m, lldx, lldy, info
    end subroutine daxpby
- end interface
+  end interface
- interface
+  interface
-      subroutine dcsmm(trans,m,n,k,alpha,pl,fida,descra,a,ia1,ia2,&
+    subroutine dcsmm(trans,m,n,k,alpha,pl,fida,descra,a,ia1,ia2,&
-     &   infoa,pr,b,ldb,beta,c,ldc,work,lwork,ierror)
+         &   infoa,pr,b,ldb,beta,c,ldc,work,lwork,ierror)
      use psb_const_mod
      integer          ::  m,n,k,ldb,ldc,lwork, ierror
      character        ::  trans
-      real(kind(1.d0)) ::  alpha,beta
+      real(psb_dpk_) ::  alpha,beta
      integer          ::  ia1(*),ia2(*),infoa(*),pl(*),pr(*)
      character        ::  descra*11, fida*5
-      real(kind(1.d0)) ::  a(*),b(ldb,*),c(ldc,*),work(*)
+      real(psb_dpk_) ::  a(*),b(ldb,*),c(ldc,*),work(*)
    end subroutine dcsmm
- end interface
+  end interface
- interface
+  interface
-    real(kind(1.d0)) function dcsnmi(trans,m,n,fida,descra,a,ia1,ia2,&
+    real(psb_dpk_) function dcsnmi(trans,m,n,fida,descra,a,ia1,ia2,&
         &                 infoa,ierror)
      use psb_const_mod
      integer          ::  m,n, ierror
      character        ::  trans
      integer          ::  ia1(*),ia2(*),infoa(*)
      character        ::  descra*11, fida*5
-      real(kind(1.d0)) ::  a(*)
+      real(psb_dpk_) ::  a(*)
    end function dcsnmi
- end interface
+  end interface
- interface
+  interface
-      subroutine dcsrp(trans,m,n,fida,descra,ia1,ia2,infoa,&
+    subroutine dcsrp(trans,m,n,fida,descra,ia1,ia2,infoa,&
-     &  p,work,lwork,ierror)
+         &  p,work,lwork,ierror)
      use psb_const_mod
      implicit none                                                      
      integer          :: lwork, m, n, ierror
      character        :: trans
-      real(kind(1.d0)) :: work(lwork)
+      real(psb_dpk_) :: work(lwork)
      integer          :: ia1(*), ia2(*), infoa(*), p(*), int_val(5)
      character        :: descra*11, fida*5
    end subroutine dcsrp
- end interface
+  end interface
- interface
+  interface
-      subroutine dcssm(trans,m,n,alpha,unitd,d,&
+    subroutine dcssm(trans,m,n,alpha,unitd,d,&
-     &   pl,fidt,descrt,t,it1,it2,infot,pr,&
+         &   pl,fidt,descrt,t,it1,it2,infot,pr,&
-     &   b,ldb,beta,c,ldc,work,lwork,ierror)
+         &   b,ldb,beta,c,ldc,work,lwork,ierror)
-      real(kind(1.d0)) :: alpha, beta
+      use psb_const_mod
      real(psb_dpk_) :: alpha, beta
      integer          :: n, ldb, ldc, m, lwork, ierror
      character        :: unitd, trans
-      real(kind(1.d0)) :: t(*), b(ldb,*), c(ldc,*), d(*), work(*)
+      real(psb_dpk_) :: t(*), b(ldb,*), c(ldc,*), d(*), work(*)
      integer          :: it1(*), it2(*), infot(*), pl(*), pr(*)
      character        :: descrt*11, fidt*5
    end subroutine dcssm
- end interface
+  end interface
- interface
+  interface
    subroutine dcsupd(m, n, fida, descra, a, ia1,&
         &  ia2, infoa, ia, ja, fidh, descrh, h, ih1, ih2,&
         &  infoh, ih, jh, flag, glob_to_loc,&
         &  iwork, liwork, ierror)
      use psb_const_mod
      integer          :: ia, ja, ih, jh, m, n,&
           &  ierror, flag, liwork 
      integer          :: ia1(*),ia2(*),ih1(*),ih2(*),&
           &  infoa(*),infoh(*),iwork(*),&
           &  glob_to_loc(*)
      character        :: descra*11,descrh*11, fida*5, fidh*5
-      real(kind(1.d0)) :: a(*),h(*)    
+      real(psb_dpk_) :: a(*),h(*)    
    end subroutine dcsupd
- end interface
+  end interface
- interface
+  interface
    subroutine dgelp(trans,m,n,p,b,ldb,work,lwork,ierror)
      use psb_const_mod
      integer          :: ldb, m, n, lwork, ierror
      character        :: trans
-      real(kind(1.d0)) :: b(ldb,*), work(*)
+      real(psb_dpk_) :: b(ldb,*), work(*)
      integer          :: p(*)
    end subroutine dgelp
- end interface
+  end interface
- interface
+  interface
    subroutine dlpupd(m,n,perm,b,ldb,beta,c,ldc)
      use psb_const_mod
      integer          :: m, n, ldb, ldc
-      real(kind(1.d0)) :: beta
+      real(psb_dpk_) :: beta
      integer          :: perm(*)
-      real(kind(1.d0)) :: b(ldb,*), c(ldc,*)
+      real(psb_dpk_) :: b(ldb,*), c(ldc,*)
    end subroutine dlpupd
- end interface
+  end interface
- interface
+  interface
-      subroutine dswmm(trans,m,n,k,alpha,fida,descra,a,ia1,ia2,&
+    subroutine dswmm(trans,m,n,k,alpha,fida,descra,a,ia1,ia2,&
-           & infoa,b,ldb,beta,c,ldc,work,lwork,ierror)
+         & infoa,b,ldb,beta,c,ldc,work,lwork,ierror)
      use psb_const_mod
      integer          :: m,n,k,ldb,ldc,lwork,ierror
      character        :: trans
-      real(kind(1.d0)) :: alpha,beta
+      real(psb_dpk_) :: alpha,beta
      integer          :: ia1(*),ia2(*),infoa(*), int_val(5)
      character        :: descra*11, fida*5
-      real(kind(1.d0)) :: a(*),b(ldb,*),c(ldc,*),work(*)
+      real(psb_dpk_) :: a(*),b(ldb,*),c(ldc,*),work(*)
    end subroutine dswmm
- end interface
+  end interface
- interface
+  interface
    subroutine dswprt(m,n,fida,descra,a,ia1,ia2,infoa,title,&
         & iout,ierror)
      use psb_const_mod
      integer       :: m,n,iout,ierror
      integer       :: ia1(*),ia2(*),infoa(*)
      character     :: descra*11, fida*5, title*(*)
    end subroutine dswprt
- end interface
+  end interface
- interface
+  interface
-      subroutine dswsm(trans,m,n,alpha,unitd,d,fidt,descrt,t,it1,it2,& 
+    subroutine dswsm(trans,m,n,alpha,unitd,d,fidt,descrt,t,it1,it2,& 
-           & infot,b,ldb,beta,c,ldc,work,lwork,ierror)
+         & infot,b,ldb,beta,c,ldc,work,lwork,ierror)
      use psb_const_mod
      integer          :: m, n, ldb, ldc, lwork, ierror
      character        :: unitd, trans
-      real(kind(1.d0)) :: alpha, beta
+      real(psb_dpk_) :: alpha, beta
      integer          :: it1(*), it2(*), infot(*)
      character        :: descrt*11, fidt*5
-      real(kind(1.d0)) :: t(*), b(ldb,*), c(ldc,*), d(*), work(*)
+      real(psb_dpk_) :: t(*), b(ldb,*), c(ldc,*), d(*), work(*)
    end subroutine dswsm
- end interface
+  end interface
- interface
+  interface
-      subroutine symbmm (n, m, l, ia, ja, diaga, ib,&
+    subroutine symbmm (n, m, l, ia, ja, diaga, ib,&
-           & jb, diagb, ic, jc, diagc, index)
+         & jb, diagb, ic, jc, diagc, index)
      integer          :: ia(*), ja(*), diaga, ib(*),&
           & jb(*), diagb, diagc, index(*)
      integer, pointer :: ic(:),jc(:)
      integer          :: nze, info
    end subroutine symbmm
- end interface
+  end interface
 end module psb_spsb_mod
--- a/base/modules/psb_tools_mod.f90
+++ b/base/modules/psb_tools_mod.f90
@ -38,7 +38,7 @@ Module psb_tools_mod
    subroutine psb_dalloc(x, desc_a, info, n, lb)
      use psb_descriptor_type
      implicit none
-      real(kind(1.d0)), allocatable, intent(out) :: x(:,:)
+      real(psb_dpk_), allocatable, intent(out) :: x(:,:)
      type(psb_desc_type), intent(in) :: desc_a
      integer,intent(out)             :: info
      integer, optional, intent(in)   :: n, lb
@ -46,7 +46,7 @@ Module psb_tools_mod
    ! 1-D double precision version
    subroutine psb_dallocv(x, desc_a,info,n)
      use psb_descriptor_type
-      real(kind(1.d0)), allocatable, intent(out)       :: x(:)
+      real(psb_dpk_), allocatable, intent(out)       :: x(:)
      type(psb_desc_type), intent(in) :: desc_a
      integer,intent(out)             :: info
      integer, optional, intent(in)   :: n
@ -70,7 +70,7 @@ Module psb_tools_mod
    subroutine psb_zalloc(x, desc_a, info, n, lb)
      use psb_descriptor_type
      implicit none
-      complex(kind(1.d0)), allocatable, intent(out)    :: x(:,:)
+      complex(psb_dpk_), allocatable, intent(out)    :: x(:,:)
      type(psb_desc_type), intent(in) :: desc_a
      integer                         :: info
      integer, optional, intent(in)   :: n, lb
@ -78,7 +78,7 @@ Module psb_tools_mod
    ! 1-D double complex version
    subroutine psb_zallocv(x, desc_a,info,n)
      use psb_descriptor_type
-      complex(kind(1.d0)), allocatable, intent(out)    :: x(:)
+      complex(psb_dpk_), allocatable, intent(out)    :: x(:)
      type(psb_desc_type), intent(in) :: desc_a
      integer                         :: info
      integer, optional, intent(in)   :: n
@ -91,14 +91,14 @@ Module psb_tools_mod
    subroutine psb_dasb(x, desc_a, info)
      use psb_descriptor_type
      type(psb_desc_type), intent(in) ::  desc_a
-      real(kind(1.d0)), allocatable, intent(inout)       ::  x(:,:)
+      real(psb_dpk_), allocatable, intent(inout)       ::  x(:,:)
      integer, intent(out)            ::  info
    end subroutine psb_dasb
    ! 1-D double precision version
    subroutine psb_dasbv(x, desc_a, info)
      use psb_descriptor_type
      type(psb_desc_type), intent(in) ::  desc_a
-      real(kind(1.d0)), allocatable, intent(inout)   ::  x(:)
+      real(psb_dpk_), allocatable, intent(inout)   ::  x(:)
      integer, intent(out)        ::  info
    end subroutine psb_dasbv
    ! 2-D integer version
@ -119,14 +119,14 @@ Module psb_tools_mod
    subroutine psb_zasb(x, desc_a, info)
      use psb_descriptor_type
      type(psb_desc_type), intent(in) ::  desc_a
-      complex(kind(1.d0)), allocatable, intent(inout)       ::  x(:,:)
+      complex(psb_dpk_), allocatable, intent(inout)       ::  x(:,:)
      integer, intent(out)            ::  info
    end subroutine psb_zasb
    ! 1-D double precision version
    subroutine psb_zasbv(x, desc_a, info)
      use psb_descriptor_type
      type(psb_desc_type), intent(in) ::  desc_a
-      complex(kind(1.d0)), allocatable, intent(inout)   ::  x(:)
+      complex(psb_dpk_), allocatable, intent(inout)   ::  x(:)
      integer, intent(out)        ::  info
    end subroutine psb_zasbv
  end interface
@ -175,14 +175,14 @@ Module psb_tools_mod
    ! 2-D double precision version
    subroutine psb_dfree(x, desc_a, info)
      use psb_descriptor_type
-      real(kind(1.d0)),allocatable, intent(inout)        :: x(:,:)
+      real(psb_dpk_),allocatable, intent(inout)        :: x(:,:)
      type(psb_desc_type), intent(in) :: desc_a
      integer                         :: info
    end subroutine psb_dfree
    ! 1-D double precision version
    subroutine psb_dfreev(x, desc_a, info)
      use psb_descriptor_type
-      real(kind(1.d0)),allocatable, intent(inout)        :: x(:)
+      real(psb_dpk_),allocatable, intent(inout)        :: x(:)
      type(psb_desc_type), intent(in) :: desc_a
      integer                         :: info
    end subroutine psb_dfreev
@ -203,14 +203,14 @@ Module psb_tools_mod
    ! 2-D double precision version
    subroutine psb_zfree(x, desc_a, info)
      use psb_descriptor_type
-      complex(kind(1.d0)),allocatable, intent(inout)        :: x(:,:)
+      complex(psb_dpk_),allocatable, intent(inout)        :: x(:,:)
      type(psb_desc_type), intent(in) :: desc_a
      integer                         :: info
    end subroutine psb_zfree
    ! 1-D double precision version
    subroutine psb_zfreev(x, desc_a, info)
      use psb_descriptor_type
-      complex(kind(1.d0)),allocatable, intent(inout)        :: x(:)
+      complex(psb_dpk_),allocatable, intent(inout)        :: x(:)
      type(psb_desc_type), intent(in) :: desc_a
      integer                         :: info
    end subroutine psb_zfreev
@ -223,9 +223,9 @@ Module psb_tools_mod
      use psb_descriptor_type
      integer, intent(in)                ::  m
      type(psb_desc_type), intent(in)    ::  desc_a
-      real(kind(1.d0)),intent(inout)           ::  x(:,:)
+      real(psb_dpk_),intent(inout)           ::  x(:,:)
      integer, intent(in)                ::  irw(:)
-      real(kind(1.d0)), intent(in)       ::  val(:,:)
+      real(psb_dpk_), intent(in)       ::  val(:,:)
      integer, intent(out)               ::  info
      integer, optional, intent(in)      ::  dupl
    end subroutine psb_dinsi
@ -234,9 +234,9 @@ Module psb_tools_mod
      use psb_descriptor_type
      integer, intent(in)                ::  m
      type(psb_desc_type), intent(in)    ::  desc_a
-      real(kind(1.d0)),intent(inout)     ::  x(:)
+      real(psb_dpk_),intent(inout)     ::  x(:)
      integer, intent(in)                ::  irw(:)
-      real(kind(1.d0)), intent(in)       ::  val(:)
+      real(psb_dpk_), intent(in)       ::  val(:)
      integer, intent(out)               ::  info
      integer, optional, intent(in)      ::  dupl
    end subroutine psb_dinsvi
@ -267,9 +267,9 @@ Module psb_tools_mod
      use psb_descriptor_type
      integer, intent(in)              ::  m
      type(psb_desc_type), intent(in)  ::  desc_a
-      complex(kind(1.d0)),intent(inout)      ::  x(:,:)
+      complex(psb_dpk_),intent(inout)      ::  x(:,:)
      integer, intent(in)              ::  irw(:)
-      complex(kind(1.d0)), intent(in)  ::  val(:,:)
+      complex(psb_dpk_), intent(in)  ::  val(:,:)
      integer, intent(out)             ::  info
      integer, optional, intent(in)    ::  dupl
    end subroutine psb_zinsi
@ -278,9 +278,9 @@ Module psb_tools_mod
      use psb_descriptor_type
      integer, intent(in)              ::  m
      type(psb_desc_type), intent(in)  ::  desc_a
-      complex(kind(1.d0)),intent(inout)      ::  x(:)
+      complex(psb_dpk_),intent(inout)      ::  x(:)
      integer, intent(in)              ::  irw(:)
-      complex(kind(1.d0)), intent(in)  ::  val(:)
+      complex(psb_dpk_), intent(in)  ::  val(:)
      integer, intent(out)             ::  info
      integer, optional, intent(in)    ::  dupl
    end subroutine psb_zinsvi
@ -426,7 +426,7 @@ Module psb_tools_mod
      type(psb_desc_type), intent(inout)   :: desc_a
      type(psb_dspmat_type), intent(inout) :: a
      integer, intent(in)                  :: nz,ia(:),ja(:)
-      real(kind(1.d0)), intent(in)         :: val(:)
+      real(psb_dpk_), intent(in)         :: val(:)
      integer, intent(out)                 :: info
      logical, intent(in), optional        :: rebuild
    end subroutine psb_dspins
@ -446,7 +446,7 @@ Module psb_tools_mod
      type(psb_desc_type), intent(inout)   :: desc_a
      type(psb_zspmat_type), intent(inout) :: a
      integer, intent(in)                  :: nz,ia(:),ja(:)
-      complex(kind(1.d0)), intent(in)      :: val(:)
+      complex(psb_dpk_), intent(in)      :: val(:)
      integer, intent(out)                 :: info
      logical, intent(in), optional        :: rebuild
    end subroutine psb_zspins
--- a/base/modules/psi_mod.f90
+++ b/base/modules/psi_mod.f90
@ -100,8 +100,8 @@ module psi_mod
      use psb_descriptor_type
      integer, intent(in)  :: flag, n
      integer, intent(out) :: info
-      real(kind(1.d0))     :: y(:,:), beta
+      real(psb_dpk_)     :: y(:,:), beta
-      real(kind(1.d0)),target      :: work(:)
+      real(psb_dpk_),target      :: work(:)
      type(psb_desc_type), target  :: desc_a
      integer, optional    :: data
    end subroutine psi_dswapdatam
@ -109,23 +109,25 @@ module psi_mod
      use psb_descriptor_type
      integer, intent(in)  :: flag
      integer, intent(out) :: info
-      real(kind(1.d0))     :: y(:), beta 
+      real(psb_dpk_)     :: y(:), beta 
-      real(kind(1.d0)),target      :: work(:)
+      real(psb_dpk_),target      :: work(:)
      type(psb_desc_type), target  :: desc_a
      integer, optional    :: data
    end subroutine psi_dswapdatav
    subroutine psi_dswapidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work,info)
      use psb_const_mod
      integer, intent(in)  :: ictxt,icomm,flag, n
      integer, intent(out) :: info
-      real(kind(1.d0))     :: y(:,:), beta
+      real(psb_dpk_)     :: y(:,:), beta
-      real(kind(1.d0)),target :: work(:)
+      real(psb_dpk_),target :: work(:)
      integer, intent(in)  :: idx(:),totxch,totsnd,totrcv
    end subroutine psi_dswapidxm
    subroutine psi_dswapidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,info)
      use psb_const_mod
      integer, intent(in)  :: ictxt,icomm,flag
      integer, intent(out) :: info
-      real(kind(1.d0))     :: y(:), beta
+      real(psb_dpk_)     :: y(:), beta
-      real(kind(1.d0)),target :: work(:)
+      real(psb_dpk_),target :: work(:)
      integer, intent(in)  :: idx(:),totxch,totsnd,totrcv
    end subroutine psi_dswapidxv
    subroutine psi_iswapdatam(flag,n,beta,y,desc_a,work,info,data)
@ -147,6 +149,7 @@ module psi_mod
      integer, optional    :: data
    end subroutine psi_iswapdatav
    subroutine psi_iswapidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work,info)
      use psb_const_mod
      integer, intent(in)  :: ictxt,icomm,flag, n
      integer, intent(out) :: info
      integer              :: y(:,:), beta
@ -154,6 +157,7 @@ module psi_mod
      integer, intent(in)  :: idx(:),totxch,totsnd,totrcv
    end subroutine psi_iswapidxm
    subroutine psi_iswapidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,info)
      use psb_const_mod
      integer, intent(in)  :: ictxt,icomm,flag
      integer, intent(out) :: info
      integer              :: y(:), beta
@ -164,8 +168,8 @@ module psi_mod
      use psb_descriptor_type
      integer, intent(in)  :: flag, n
      integer, intent(out) :: info
-      complex(kind(1.d0))     :: y(:,:), beta
+      complex(psb_dpk_)     :: y(:,:), beta
-      complex(kind(1.d0)),target   :: work(:)
+      complex(psb_dpk_),target   :: work(:)
      type(psb_desc_type), target  :: desc_a
      integer, optional    :: data
    end subroutine psi_zswapdatam
@ -173,23 +177,25 @@ module psi_mod
      use psb_descriptor_type
      integer, intent(in)  :: flag
      integer, intent(out) :: info
-      complex(kind(1.d0))     :: y(:), beta
+      complex(psb_dpk_)     :: y(:), beta
-      complex(kind(1.d0)),target   :: work(:)
+      complex(psb_dpk_),target   :: work(:)
      type(psb_desc_type), target  :: desc_a
      integer, optional    :: data
    end subroutine psi_zswapdatav
    subroutine psi_zswapidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work,info)
      use psb_const_mod
      integer, intent(in)  :: ictxt,icomm,flag, n
      integer, intent(out) :: info
-      complex(kind(1.d0))     :: y(:,:), beta
+      complex(psb_dpk_)     :: y(:,:), beta
-      complex(kind(1.d0)),target :: work(:)
+      complex(psb_dpk_),target :: work(:)
      integer, intent(in)  :: idx(:),totxch,totsnd,totrcv
    end subroutine psi_zswapidxm
    subroutine psi_zswapidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,info)
      use psb_const_mod
      integer, intent(in)  :: ictxt,icomm,flag
      integer, intent(out) :: info
-      complex(kind(1.d0))     :: y(:), beta
+      complex(psb_dpk_)     :: y(:), beta
-      complex(kind(1.d0)),target :: work(:)
+      complex(psb_dpk_),target :: work(:)
      integer, intent(in)  :: idx(:),totxch,totsnd,totrcv
    end subroutine psi_zswapidxv
  end interface
@ -200,8 +206,8 @@ module psi_mod
      use psb_descriptor_type
      integer, intent(in)  :: flag, n
      integer, intent(out) :: info
-      real(kind(1.d0))     :: y(:,:), beta
+      real(psb_dpk_)     :: y(:,:), beta
-      real(kind(1.d0)),target     :: work(:)
+      real(psb_dpk_),target     :: work(:)
      type(psb_desc_type), target  :: desc_a
      integer, optional    :: data
    end subroutine psi_dswaptranm
@ -209,23 +215,25 @@ module psi_mod
      use psb_descriptor_type
      integer, intent(in)  :: flag
      integer, intent(out) :: info
-      real(kind(1.d0))     :: y(:), beta
+      real(psb_dpk_)     :: y(:), beta
-      real(kind(1.d0)),target     :: work(:)
+      real(psb_dpk_),target     :: work(:)
      type(psb_desc_type), target  :: desc_a
      integer, optional    :: data
    end subroutine psi_dswaptranv
    subroutine psi_dtranidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work,info)
      use psb_const_mod
      integer, intent(in)  :: ictxt,icomm,flag, n
      integer, intent(out) :: info
-      real(kind(1.d0))     :: y(:,:), beta
+      real(psb_dpk_)     :: y(:,:), beta
-      real(kind(1.d0)),target :: work(:)
+      real(psb_dpk_),target :: work(:)
      integer, intent(in)  :: idx(:),totxch,totsnd,totrcv
    end subroutine psi_dtranidxm
    subroutine psi_dtranidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,info)
      use psb_const_mod
      integer, intent(in)  :: ictxt,icomm,flag
      integer, intent(out) :: info
-      real(kind(1.d0))     :: y(:), beta
+      real(psb_dpk_)     :: y(:), beta
-      real(kind(1.d0)),target :: work(:)
+      real(psb_dpk_),target :: work(:)
      integer, intent(in)  :: idx(:),totxch,totsnd,totrcv
    end subroutine psi_dtranidxv
    subroutine psi_iswaptranm(flag,n,beta,y,desc_a,work,info,data)
@ -247,6 +255,7 @@ module psi_mod
      integer, optional    :: data
    end subroutine psi_iswaptranv
    subroutine psi_itranidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work,info)
      use psb_const_mod
      integer, intent(in)  :: ictxt,icomm,flag, n
      integer, intent(out) :: info
      integer              :: y(:,:), beta
@ -254,6 +263,7 @@ module psi_mod
      integer, intent(in)  :: idx(:),totxch,totsnd,totrcv
    end subroutine psi_itranidxm
    subroutine psi_itranidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,info)
      use psb_const_mod
      integer, intent(in)  :: ictxt,icomm,flag
      integer, intent(out) :: info
      integer              :: y(:), beta
@ -264,8 +274,8 @@ module psi_mod
      use psb_descriptor_type
      integer, intent(in)  :: flag, n
      integer, intent(out) :: info
-      complex(kind(1.d0))     :: y(:,:), beta
+      complex(psb_dpk_)     :: y(:,:), beta
-      complex(kind(1.d0)),target   :: work(:)
+      complex(psb_dpk_),target   :: work(:)
      type(psb_desc_type), target  :: desc_a
      integer, optional    :: data
    end subroutine psi_zswaptranm
@ -273,23 +283,25 @@ module psi_mod
      use psb_descriptor_type
      integer, intent(in)  :: flag
      integer, intent(out) :: info
-      complex(kind(1.d0))     :: y(:), beta
+      complex(psb_dpk_)     :: y(:), beta
-      complex(kind(1.d0)),target   :: work(:)       
+      complex(psb_dpk_),target   :: work(:)       
      type(psb_desc_type), target  :: desc_a
      integer, optional    :: data
    end subroutine psi_zswaptranv
    subroutine psi_ztranidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work,info)
      use psb_const_mod
      integer, intent(in)  :: ictxt,icomm,flag, n
      integer, intent(out) :: info
-      complex(kind(1.d0))     :: y(:,:), beta
+      complex(psb_dpk_)     :: y(:,:), beta
-      complex(kind(1.d0)), target :: work(:)
+      complex(psb_dpk_), target :: work(:)
      integer, intent(in)  :: idx(:),totxch,totsnd,totrcv
    end subroutine psi_ztranidxm
    subroutine psi_ztranidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,info)
      use psb_const_mod
      integer, intent(in)  :: ictxt,icomm,flag
      integer, intent(out) :: info
-      complex(kind(1.d0))     :: y(:), beta
+      complex(psb_dpk_)     :: y(:), beta
-      complex(kind(1.d0)), target :: work(:)
+      complex(psb_dpk_), target :: work(:)
      integer, intent(in)  :: idx(:),totxch,totsnd,totrcv
    end subroutine psi_ztranidxv
  end interface
@ -678,7 +690,7 @@ contains
    use psb_penv_mod
    implicit none
-    real(kind(1.d0)), intent(inout), target :: x(:)
+    real(psb_dpk_), intent(inout), target :: x(:)
    type(psb_desc_type), intent(in)         :: desc_a
    integer, intent(in)                     :: update
    integer, intent(out)                    :: info
@ -750,7 +762,7 @@ contains
    use psb_penv_mod
    implicit none
-    real(kind(1.d0)), intent(inout), target :: x(:,:)
+    real(psb_dpk_), intent(inout), target :: x(:,:)
    type(psb_desc_type), intent(in)         :: desc_a
    integer, intent(in)                     :: update
    integer, intent(out)                    :: info
@ -821,7 +833,7 @@ contains
    use psb_penv_mod
    implicit none
-    complex(kind(1.d0)), intent(inout), target :: x(:)
+    complex(psb_dpk_), intent(inout), target :: x(:)
    type(psb_desc_type), intent(in)         :: desc_a
    integer, intent(in)                     :: update
    integer, intent(out)                    :: info
@ -893,7 +905,7 @@ contains
    use psb_penv_mod
    implicit none
-    complex(kind(1.d0)), intent(inout), target :: x(:,:)
+    complex(psb_dpk_), intent(inout), target :: x(:,:)
    type(psb_desc_type), intent(in)         :: desc_a
    integer, intent(in)                     :: update
    integer, intent(out)                    :: info
@ -1112,8 +1124,8 @@ contains
    use psb_penv_mod
    implicit none
-    real(kind(1.d0)), intent(inout)  :: x(:)
+    real(psb_dpk_), intent(inout)  :: x(:)
-    real(kind(1.d0)), allocatable    :: xs(:)
+    real(psb_dpk_), allocatable    :: xs(:)
    type(psb_desc_type), intent(in)  :: desc_a
    integer, intent(out)             :: info
@ -1166,8 +1178,8 @@ contains
    use psb_penv_mod
    implicit none
-    real(kind(1.d0)), intent(inout)  :: x(:)
+    real(psb_dpk_), intent(inout)  :: x(:)
-    real(kind(1.d0))                 :: xs(:)
+    real(psb_dpk_)                 :: xs(:)
    type(psb_desc_type), intent(in)  :: desc_a
    integer, intent(out)             :: info
@ -1216,8 +1228,8 @@ contains
    use psb_penv_mod
    implicit none
-    real(kind(1.d0)), intent(inout)  :: x(:,:)
+    real(psb_dpk_), intent(inout)  :: x(:,:)
-    real(kind(1.d0)), allocatable    :: xs(:,:)
+    real(psb_dpk_), allocatable    :: xs(:,:)
    type(psb_desc_type), intent(in)  :: desc_a
    integer, intent(out)             :: info
@ -1271,8 +1283,8 @@ contains
    use psb_penv_mod
    implicit none
-    real(kind(1.d0)), intent(inout)  :: x(:,:)
+    real(psb_dpk_), intent(inout)  :: x(:,:)
-    real(kind(1.d0))                 :: xs(:,:)
+    real(psb_dpk_)                 :: xs(:,:)
    type(psb_desc_type), intent(in)  :: desc_a
    integer, intent(out)             :: info
@ -1329,8 +1341,8 @@ contains
    use psb_penv_mod
    implicit none
-    complex(kind(1.d0)), intent(inout)  :: x(:)
+    complex(psb_dpk_), intent(inout)  :: x(:)
-    complex(kind(1.d0)), allocatable    :: xs(:)
+    complex(psb_dpk_), allocatable    :: xs(:)
    type(psb_desc_type), intent(in)  :: desc_a
    integer, intent(out)             :: info
@ -1383,8 +1395,8 @@ contains
    use psb_penv_mod
    implicit none
-    complex(kind(1.d0)), intent(inout)  :: x(:)
+    complex(psb_dpk_), intent(inout)  :: x(:)
-    complex(kind(1.d0))                 :: xs(:)
+    complex(psb_dpk_)                 :: xs(:)
    type(psb_desc_type), intent(in)  :: desc_a
    integer, intent(out)             :: info
@ -1433,8 +1445,8 @@ contains
    use psb_penv_mod
    implicit none
-    complex(kind(1.d0)), intent(inout)  :: x(:,:)
+    complex(psb_dpk_), intent(inout)  :: x(:,:)
-    complex(kind(1.d0)), allocatable    :: xs(:,:)
+    complex(psb_dpk_), allocatable    :: xs(:,:)
    type(psb_desc_type), intent(in)  :: desc_a
    integer, intent(out)             :: info
@ -1488,8 +1500,8 @@ contains
    use psb_penv_mod
    implicit none
-    complex(kind(1.d0)), intent(inout)  :: x(:,:)
+    complex(psb_dpk_), intent(inout)  :: x(:,:)
-    complex(kind(1.d0))                 :: xs(:,:)
+    complex(psb_dpk_)                 :: xs(:,:)
    type(psb_desc_type), intent(in)  :: desc_a
    integer, intent(out)             :: info
@ -1759,7 +1771,7 @@ contains
 !!$    implicit none
 !!$
 !!$    integer :: n, k, idx(:)
-!!$    real(kind(1.d0)) :: x(:,:), y(:)
+!!$    real(psb_dpk_) :: x(:,:), y(:)
 !!$
 !!$    ! Locals
 !!$    integer :: i, j, pt
@ -1780,7 +1792,7 @@ contains
 !!$    implicit none
 !!$
 !!$    integer :: n, idx(:)
-!!$    real(kind(1.d0)) :: x(:), y(:)
+!!$    real(psb_dpk_) :: x(:), y(:)
 !!$
 !!$    ! Locals
 !!$    integer :: i
@ -1798,7 +1810,7 @@ contains
 !!$    implicit none
 !!$
 !!$    integer :: n, k, idx(:)
-!!$    real(kind(1.d0)) :: beta, x(:), y(:,:)
+!!$    real(psb_dpk_) :: beta, x(:), y(:,:)
 !!$
 !!$    ! Locals
 !!$    integer :: i, j, pt
@ -1836,7 +1848,7 @@ contains
 !!$    implicit none
 !!$
 !!$    integer :: n, idx(:)
-!!$    real(kind(1.d0)) :: beta, x(:), y(:)
+!!$    real(psb_dpk_) :: beta, x(:), y(:)
 !!$
 !!$    ! Locals
 !!$    integer :: i
@ -1969,7 +1981,7 @@ contains
 !!$    implicit none
 !!$
 !!$    integer :: n, k, idx(:)
-!!$    complex(kind(1.d0)) :: x(:,:), y(:)
+!!$    complex(psb_dpk_) :: x(:,:), y(:)
 !!$
 !!$    ! Locals
 !!$    integer :: i, j, pt
@ -1991,7 +2003,7 @@ contains
 !!$    implicit none
 !!$
 !!$    integer :: n, idx(:)
-!!$    complex(kind(1.d0)) :: x(:), y(:)
+!!$    complex(psb_dpk_) :: x(:), y(:)
 !!$
 !!$    ! Locals
 !!$    integer :: i
@ -2008,7 +2020,7 @@ contains
 !!$    implicit none
 !!$
 !!$    integer :: n, k, idx(:)
-!!$    complex(kind(1.d0)) :: beta, x(:), y(:,:)
+!!$    complex(psb_dpk_) :: beta, x(:), y(:,:)
 !!$
 !!$    ! Locals
 !!$    integer :: i, j, pt
@ -2047,7 +2059,7 @@ contains
 !!$    implicit none
 !!$
 !!$    integer :: n, idx(:)
-!!$    complex(kind(1.d0)) :: beta, x(:), y(:)
+!!$    complex(psb_dpk_) :: beta, x(:), y(:)
 !!$
 !!$    ! Locals
 !!$    integer :: i
--- a/base/modules/psi_serial_mod.f90
+++ b/base/modules/psi_serial_mod.f90
@ -105,7 +105,7 @@ contains
    implicit none
    integer :: n, idx(:)
-    real(kind(1.d0)) :: x(:), y(:), alpha, beta
+    real(psb_dpk_) :: x(:), y(:), alpha, beta
    ! Locals
    integer :: i
@ -155,7 +155,7 @@ contains
    implicit none
    integer :: n, idx(:)
-    complex(kind(1.d0)) :: x(:), y(:),alpha,beta
+    complex(psb_dpk_) :: x(:), y(:),alpha,beta
    ! Locals
    integer :: i
@ -208,7 +208,7 @@ contains
    implicit none
    integer :: n, k, idx(:)
-    real(kind(1.d0)) :: x(:,:), y(:)
+    real(psb_dpk_) :: x(:,:), y(:)
    ! Locals
    integer :: i, j, pt
@ -252,7 +252,7 @@ contains
    implicit none
    integer :: n, k, idx(:)
-    complex(kind(1.d0)) :: x(:,:), y(:)
+    complex(psb_dpk_) :: x(:,:), y(:)
    ! Locals
    integer :: i, j, pt
@ -273,7 +273,7 @@ contains
    implicit none
    integer :: n, idx(:)
-    real(kind(1.d0)) :: x(:), y(:)
+    real(psb_dpk_) :: x(:), y(:)
    ! Locals
    integer :: i
@ -307,7 +307,7 @@ contains
    implicit none
    integer :: n, idx(:)
-    complex(kind(1.d0)) :: x(:), y(:)
+    complex(psb_dpk_) :: x(:), y(:)
    ! Locals
    integer :: i
@ -325,7 +325,7 @@ contains
    implicit none
    integer :: n, k, idx(:)
-    real(kind(1.d0)) :: beta, x(:), y(:,:)
+    real(psb_dpk_) :: beta, x(:), y(:,:)
    ! Locals
    integer :: i, j, pt
@ -363,7 +363,7 @@ contains
    implicit none
    integer :: n, idx(:)
-    real(kind(1.d0)) :: beta, x(:), y(:)
+    real(psb_dpk_) :: beta, x(:), y(:)
    ! Locals
    integer :: i
@ -456,7 +456,7 @@ contains
    implicit none
    integer :: n, k, idx(:)
-    complex(kind(1.d0)) :: beta, x(:), y(:,:)
+    complex(psb_dpk_) :: beta, x(:), y(:,:)
    ! Locals
    integer :: i, j, pt
@ -495,7 +495,7 @@ contains
    implicit none
    integer :: n, idx(:)
-    complex(kind(1.d0)) :: beta, x(:), y(:)
+    complex(psb_dpk_) :: beta, x(:), y(:)
    ! Locals
    integer :: i
--- a/base/psblas/psb_damax.f90
+++ b/base/psblas/psb_damax.f90
@ -52,16 +52,16 @@ function psb_damax (x,desc_a, info, jx)
  use psb_error_mod
  implicit none
-  real(kind(1.d0)), intent(in)      :: x(:,:)
+  real(psb_dpk_), intent(in)      :: x(:,:)
  type(psb_desc_type), intent(in)   :: desc_a
  integer, intent(out)              :: info
  integer, optional, intent(in)     :: jx
-  real(kind(1.d0))                  :: psb_damax
+  real(psb_dpk_)                  :: psb_damax
  ! locals
  integer                  :: ictxt, np, me,&
       & err_act, iix, jjx, ix, ijx, m, imax, idamax
-  real(kind(1.d0))         :: amax
+  real(psb_dpk_)         :: amax
  character(len=20)        :: name, ch_err
  name='psb_damax'
@ -180,15 +180,15 @@ function psb_damaxv (x,desc_a, info)
  use psb_error_mod
  implicit none
-  real(kind(1.d0)), intent(in)      :: x(:)
+  real(psb_dpk_), intent(in)      :: x(:)
  type(psb_desc_type), intent(in)   :: desc_a
  integer, intent(out)              :: info
-  real(kind(1.d0))                  :: psb_damaxv
+  real(psb_dpk_)                  :: psb_damaxv
  ! locals
  integer                  :: ictxt, np, me,&
       & err_act, iix, jjx, jx, ix, m, imax, idamax
-  real(kind(1.d0))         :: amax
+  real(psb_dpk_)         :: amax
  character(len=20)        :: name, ch_err
  name='psb_damaxv'
@ -304,15 +304,15 @@ subroutine psb_damaxvs (res,x,desc_a, info)
  use psb_error_mod
  implicit none
-  real(kind(1.d0)), intent(in)      :: x(:)
+  real(psb_dpk_), intent(in)      :: x(:)
  type(psb_desc_type), intent(in)   :: desc_a
  integer, intent(out)              :: info
-  real(kind(1.D0)), intent(out)     :: res
+  real(psb_dpk_), intent(out)     :: res
  ! locals
  integer                  :: ictxt, np, me,&
       & err_act, iix, jjx, ix, ijx, m, imax, idamax
-  real(kind(1.d0))         :: amax
+  real(psb_dpk_)         :: amax
  character(len=20)        :: name, ch_err
  name='psb_damaxvs'
@ -426,16 +426,16 @@ subroutine psb_dmamaxs (res,x,desc_a, info,jx)
  use psb_error_mod
  implicit none
-  real(kind(1.d0)), intent(in)      :: x(:,:)
+  real(psb_dpk_), intent(in)      :: x(:,:)
  type(psb_desc_type), intent(in)   :: desc_a
  integer, intent(out)              :: info
  integer, optional, intent(in)     :: jx
-  real(kind(1.d0)), intent(out) :: res(:)
+  real(psb_dpk_), intent(out) :: res(:)
  ! locals
  integer                  :: ictxt, np, me,&
       & err_act, iix, jjx, ix, ijx, m, imax, i, k, idamax
-  real(kind(1.d0))         :: amax
+  real(psb_dpk_)         :: amax
  character(len=20)        :: name, ch_err
  name='psb_dmamaxs'
--- a/base/psblas/psb_dasum.f90
+++ b/base/psblas/psb_dasum.f90
@ -53,16 +53,16 @@ function psb_dasum (x,desc_a, info, jx)
  use psb_penv_mod
  implicit none
-  real(kind(1.d0)), intent(in)      :: x(:,:)
+  real(psb_dpk_), intent(in)      :: x(:,:)
  type(psb_desc_type), intent(in)   :: desc_a
  integer, intent(out)              :: info
  integer, optional, intent(in)     :: jx
-  real(kind(1.d0))                  :: psb_dasum
+  real(psb_dpk_)                  :: psb_dasum
  ! locals
  integer           :: ictxt, np, me, err_act, &
       &      iix, jjx, ix, ijx, m, i, idx, ndm
-  real(kind(1.d0))  :: asum, dasum
+  real(psb_dpk_)  :: asum, dasum
  character(len=20) :: name, ch_err
  name='psb_dasum'
@ -197,14 +197,14 @@ function psb_dasumv (x,desc_a, info)
  use psb_penv_mod
  implicit none
-  real(kind(1.d0)), intent(in)      :: x(:)
+  real(psb_dpk_), intent(in)      :: x(:)
  type(psb_desc_type), intent(in)   :: desc_a
  integer, intent(out)              :: info
-  real(kind(1.d0))                  :: psb_dasumv
+  real(psb_dpk_)                  :: psb_dasumv
  ! locals
  integer            :: ictxt, np, me, err_act, iix, jjx, jx, ix, m, i, idx, ndm
-  real(kind(1.d0))   :: asum, dasum
+  real(psb_dpk_)   :: asum, dasum
  character(len=20)  :: name, ch_err
  name='psb_dasumv'
@ -334,14 +334,14 @@ subroutine psb_dasumvs(res,x,desc_a, info)
  use psb_penv_mod
  implicit none
-  real(kind(1.d0)), intent(in)      :: x(:)
+  real(psb_dpk_), intent(in)      :: x(:)
-  real(kind(1.d0)), intent(out)     :: res
+  real(psb_dpk_), intent(out)     :: res
  type(psb_desc_type), intent(in)   :: desc_a
  integer, intent(out)              :: info
  ! locals
  integer            :: ictxt, np, me, err_act, iix, jjx, ix, jx, m, i, idx, ndm
-  real(kind(1.d0))   :: asum, dasum
+  real(psb_dpk_)   :: asum, dasum
  character(len=20)  :: name, ch_err
  name='psb_dasumvs'
--- a/base/psblas/psb_daxpby.f90
+++ b/base/psblas/psb_daxpby.f90
@ -61,9 +61,9 @@ subroutine  psb_daxpby(alpha, x, beta,y,desc_a,info, n, jx, jy)
  integer, intent(in), optional   :: n, jx, jy
  integer, intent(out)            :: info
  type(psb_desc_type), intent(in) :: desc_a
-  real(kind(1.D0)), intent(in)    :: alpha, beta
+  real(psb_dpk_), intent(in)    :: alpha, beta
-  real(kind(1.D0)), intent(in)    :: x(:,:)
+  real(psb_dpk_), intent(in)    :: x(:,:)
-  real(kind(1.D0)), intent(inout) :: y(:,:)
+  real(psb_dpk_), intent(inout) :: y(:,:)
  ! locals
  integer                  :: ictxt, np, me,&
@ -215,9 +215,9 @@ subroutine  psb_daxpbyv(alpha, x, beta,y,desc_a,info)
  integer, intent(out)            :: info
  type(psb_desc_type), intent(in) :: desc_a
-  real(kind(1.D0)), intent(in)    :: alpha, beta
+  real(psb_dpk_), intent(in)    :: alpha, beta
-  real(kind(1.D0)), intent(in)    :: x(:)
+  real(psb_dpk_), intent(in)    :: x(:)
-  real(kind(1.D0)), intent(inout) :: y(:)
+  real(psb_dpk_), intent(inout) :: y(:)
  ! locals
  integer                  :: ictxt, np, me,&
--- a/base/psblas/psb_ddot.f90
+++ b/base/psblas/psb_ddot.f90
@ -55,17 +55,17 @@ function psb_ddot(x, y,desc_a, info, jx, jy)
  use psb_penv_mod
  implicit none
-  real(kind(1.d0)), intent(in)     :: x(:,:), y(:,:)
+  real(psb_dpk_), intent(in)     :: x(:,:), y(:,:)
  type(psb_desc_type), intent(in)  :: desc_a
  integer, intent(in), optional    :: jx, jy
  integer, intent(out)             :: info
-  real(kind(1.D0))                 :: psb_ddot
+  real(psb_dpk_)                 :: psb_ddot
  ! locals
  integer                  :: ictxt, np, me, idx, ndm,&
       & err_act, iix, jjx, ix, ijx, iy, ijy, iiy, jjy, i, m
-  real(kind(1.D0))         :: dot_local
+  real(psb_dpk_)         :: dot_local
-  real(kind(1.d0))         :: ddot
+  real(psb_dpk_)         :: ddot
  character(len=20)        :: name, ch_err
  name='psb_ddot'
@ -208,16 +208,16 @@ function psb_ddotv(x, y,desc_a, info)
  use psb_penv_mod
  implicit none
-  real(kind(1.d0)), intent(in)     :: x(:), y(:)
+  real(psb_dpk_), intent(in)     :: x(:), y(:)
  type(psb_desc_type), intent(in)  :: desc_a
  integer, intent(out)             :: info
-  real(kind(1.D0))                 :: psb_ddotv
+  real(psb_dpk_)                 :: psb_ddotv
  ! locals
  integer                  :: ictxt, np, me, idx, ndm,&
       & err_act, iix, jjx, ix, jx, iy, jy, iiy, jjy, i, m
-  real(kind(1.D0))         :: dot_local
+  real(psb_dpk_)         :: dot_local
-  real(kind(1.d0))         :: ddot
+  real(psb_dpk_)         :: ddot
  character(len=20)        :: name, ch_err
  name='psb_ddot'
@ -345,16 +345,16 @@ subroutine psb_ddotvs(res, x, y,desc_a, info)
  use psb_penv_mod
  implicit none
-  real(kind(1.d0)), intent(in)     :: x(:), y(:)
+  real(psb_dpk_), intent(in)     :: x(:), y(:)
-  real(kind(1.d0)), intent(out)    :: res
+  real(psb_dpk_), intent(out)    :: res
  type(psb_desc_type), intent(in)  :: desc_a
  integer, intent(out)             :: info
  ! locals
  integer                  :: ictxt, np, me, idx, ndm,&
       & err_act, iix, jjx, ix, iy,  iiy, jjy, i, m
-  real(kind(1.D0))         :: dot_local
+  real(psb_dpk_)         :: dot_local
-  real(kind(1.d0))         :: ddot
+  real(psb_dpk_)         :: ddot
  character(len=20)        :: name, ch_err
  name='psb_ddot'
@ -480,16 +480,16 @@ subroutine psb_dmdots(res, x, y, desc_a, info)
  use psb_penv_mod
  implicit none
-  real(kind(1.d0)), intent(in)     :: x(:,:), y(:,:)
+  real(psb_dpk_), intent(in)     :: x(:,:), y(:,:)
-  real(kind(1.d0)), intent(out)    :: res(:)
+  real(psb_dpk_), intent(out)    :: res(:)
  type(psb_desc_type), intent(in)  :: desc_a
  integer, intent(out)             :: info
  ! locals
  integer                  :: ictxt, np, me, idx, ndm,&
       & err_act, iix, jjx, ix, iy, iiy, jjy, i, m, j, k
-  real(kind(1.d0)),allocatable  :: dot_local(:)
+  real(psb_dpk_),allocatable  :: dot_local(:)
-  real(kind(1.d0))         :: ddot
+  real(psb_dpk_)         :: ddot
  character(len=20)        :: name, ch_err
  name='psb_dmdots'
--- a/base/psblas/psb_dnrm2.f90
+++ b/base/psblas/psb_dnrm2.f90
@ -51,16 +51,16 @@ function psb_dnrm2(x, desc_a, info, jx)
  use psb_penv_mod
  implicit none
-  real(kind(1.d0)), intent(in)      ::  x(:,:)
+  real(psb_dpk_), intent(in)      ::  x(:,:)
  type(psb_desc_type), intent(in)   :: desc_a
  integer, intent(in), optional     :: jx
  integer, intent(out)              :: info
-  real(kind(1.D0))                  :: psb_dnrm2
+  real(psb_dpk_)                  :: psb_dnrm2
  ! locals
  integer                  :: ictxt, np, me,&
       & err_act, iix, jjx, ndim, ix, ijx, i, m, id, idx, ndm 
-  real(kind(1.d0))         :: nrm2, dnrm2, dd
+  real(psb_dpk_)         :: nrm2, dnrm2, dd
  external dcombnrm2
  character(len=20)        :: name, ch_err
@ -186,15 +186,15 @@ function psb_dnrm2v(x, desc_a, info)
  use psb_penv_mod
  implicit none
-  real(kind(1.d0)), intent(in)      :: x(:)
+  real(psb_dpk_), intent(in)      :: x(:)
  type(psb_desc_type), intent(in)   :: desc_a
  integer, intent(out)              :: info
-  real(kind(1.D0))                  :: psb_dnrm2v
+  real(psb_dpk_)                  :: psb_dnrm2v
  ! locals
  integer                  :: ictxt, np, me,&
       & err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm
-  real(kind(1.d0))         :: nrm2, dnrm2, dd
+  real(psb_dpk_)         :: nrm2, dnrm2, dd
  external dcombnrm2
  character(len=20)        :: name, ch_err
@ -317,15 +317,15 @@ subroutine psb_dnrm2vs(res, x, desc_a, info)
  use psb_penv_mod
  implicit none
-  real(kind(1.d0)), intent(in)      :: x(:)
+  real(psb_dpk_), intent(in)      :: x(:)
-  real(kind(1.d0)), intent(out)     :: res
+  real(psb_dpk_), intent(out)     :: res
  type(psb_desc_type), intent(in)   :: desc_a
  integer, intent(out)              :: info
  ! locals
  integer                  :: ictxt, np, me,&
       & err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm
-  real(kind(1.d0))         :: nrm2, dnrm2, dd
+  real(psb_dpk_)         :: nrm2, dnrm2, dd
  external dcombnrm2
  character(len=20)        :: name, ch_err
--- a/base/psblas/psb_dnrmi.f90
+++ b/base/psblas/psb_dnrmi.f90
@ -52,12 +52,12 @@ function psb_dnrmi(a,desc_a,info)
  type(psb_dspmat_type), intent(in)   :: a
  integer, intent(out)                :: info
  type(psb_desc_type), intent(in)     :: desc_a
-  real(kind(1.d0))                    :: psb_dnrmi
+  real(psb_dpk_)                    :: psb_dnrmi
  ! locals
  integer                  :: ictxt, np, me,&
       & err_act, n, iia, jja, ia, ja, mdim, ndim, m
-  real(kind(1.d0))         :: nrmi, dcsnmi
+  real(psb_dpk_)         :: nrmi, dcsnmi
  character(len=20)        :: name, ch_err
  name='psb_dnrmi'
--- a/base/psblas/psb_dspmm.f90
+++ b/base/psblas/psb_dspmm.f90
@ -76,13 +76,13 @@ subroutine  psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
  use psb_penv_mod
  implicit none
-  real(kind(1.D0)), intent(in)             :: alpha, beta
+  real(psb_dpk_), intent(in)             :: alpha, beta
-  real(kind(1.d0)), intent(inout), target  :: x(:,:)
+  real(psb_dpk_), intent(inout), target  :: x(:,:)
-  real(kind(1.d0)), intent(inout), target  :: y(:,:)
+  real(psb_dpk_), intent(inout), target  :: y(:,:)
  type(psb_dspmat_type), intent(in)        :: a
  type(psb_desc_type), intent(in)          :: desc_a
  integer, intent(out)                     :: info
-  real(kind(1.d0)), optional, target       :: work(:)
+  real(psb_dpk_), optional, target       :: work(:)
  character, intent(in), optional          :: trans
  integer, intent(in), optional            :: k, jx, jy
  logical, intent(in), optional            :: doswap
@ -93,8 +93,8 @@ subroutine  psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
       & m, nrow, ncol, lldx, lldy, liwork, iiy, jjy,&
       & i, ib, ib1, ip, idx
  integer, parameter       :: nb=4
-  real(kind(1.d0)), pointer     :: xp(:,:), yp(:,:), iwork(:)
+  real(psb_dpk_), pointer     :: xp(:,:), yp(:,:), iwork(:)
-  real(kind(1.d0)), allocatable :: xvsave(:,:)
+  real(psb_dpk_), allocatable :: xvsave(:,:)
  character                :: trans_
  character(len=20)        :: name, ch_err
  logical                  :: aliw, doswap_
@ -421,13 +421,13 @@ subroutine  psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
  use psb_penv_mod
  implicit none
-  real(kind(1.D0)), intent(in)             :: alpha, beta
+  real(psb_dpk_), intent(in)             :: alpha, beta
-  real(kind(1.d0)), intent(inout), target  :: x(:)
+  real(psb_dpk_), intent(inout), target  :: x(:)
-  real(kind(1.d0)), intent(inout), target  :: y(:)
+  real(psb_dpk_), intent(inout), target  :: y(:)
  type(psb_dspmat_type), intent(in)        :: a
  type(psb_desc_type), intent(in)          :: desc_a
  integer, intent(out)                     :: info
-  real(kind(1.d0)), optional, target       :: work(:)
+  real(psb_dpk_), optional, target       :: work(:)
  character, intent(in), optional          :: trans
  logical, intent(in), optional            :: doswap
@ -437,8 +437,8 @@ subroutine  psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
       & m, nrow, ncol, lldx, lldy, liwork, jx, jy, iiy, jjy,&
       & ib, ip, idx
  integer, parameter       :: nb=4
-  real(kind(1.d0)), pointer :: iwork(:), xp(:), yp(:)
+  real(psb_dpk_), pointer :: iwork(:), xp(:), yp(:)
-  real(kind(1.d0)), allocatable :: xvsave(:)
+  real(psb_dpk_), allocatable :: xvsave(:)
  character                :: trans_
  character(len=20)        :: name, ch_err
  logical                  :: aliw, doswap_
--- a/base/psblas/psb_dspsm.f90
+++ b/base/psblas/psb_dspsm.f90
@ -88,14 +88,14 @@ subroutine  psb_dspsm(alpha,a,x,beta,y,desc_a,info,&
  use psb_penv_mod
  implicit none
-  real(kind(1.D0)), intent(in)              :: alpha, beta
+  real(psb_dpk_), intent(in)              :: alpha, beta
-  real(kind(1.d0)), intent(in), target      :: x(:,:)
+  real(psb_dpk_), intent(in), target      :: x(:,:)
-  real(kind(1.d0)), intent(inout), target   :: y(:,:)
+  real(psb_dpk_), intent(inout), target   :: y(:,:)
  type (psb_dspmat_type), intent(in)        :: a
  type(psb_desc_type), intent(in)           :: desc_a
  integer, intent(out)                      :: info
-  real(kind(1.d0)), intent(in), optional, target      :: diag(:)
+  real(psb_dpk_), intent(in), optional, target      :: diag(:)
-  real(kind(1.d0)), optional, target       :: work(:)
+  real(psb_dpk_), optional, target       :: work(:)
  character, intent(in), optional           :: trans, unitd
  integer, intent(in), optional             :: choice
  integer, intent(in), optional             :: k, jx, jy
@ -108,7 +108,7 @@ subroutine  psb_dspsm(alpha,a,x,beta,y,desc_a,info,&
  character                :: lunitd
  integer, parameter       :: nb=4
-  real(kind(1.d0)),pointer :: iwork(:), xp(:,:), yp(:,:), id(:)
+  real(psb_dpk_),pointer :: iwork(:), xp(:,:), yp(:,:), id(:)
  character                :: itrans
  character(len=20)        :: name, ch_err
  logical                  :: aliw
@ -376,14 +376,14 @@ subroutine  psb_dspsv(alpha,a,x,beta,y,desc_a,info,&
  use psb_penv_mod
  implicit none 
-  real(kind(1.D0)), intent(in)              :: alpha, beta
+  real(psb_dpk_), intent(in)              :: alpha, beta
-  real(kind(1.d0)), intent(in), target      :: x(:)
+  real(psb_dpk_), intent(in), target      :: x(:)
-  real(kind(1.d0)), intent(inout), target   :: y(:)
+  real(psb_dpk_), intent(inout), target   :: y(:)
  type(psb_dspmat_type), intent(in)         :: a
  type(psb_desc_type), intent(in)           :: desc_a
  integer, intent(out)                      :: info
-  real(kind(1.d0)), intent(in), optional, target    :: diag(:)
+  real(psb_dpk_), intent(in), optional, target    :: diag(:)
-  real(kind(1.d0)), optional, target        :: work(:)
+  real(psb_dpk_), optional, target        :: work(:)
  character, intent(in), optional           :: trans, unitd
  integer, intent(in), optional             :: choice
@ -395,7 +395,7 @@ subroutine  psb_dspsv(alpha,a,x,beta,y,desc_a,info,&
  character                :: lunitd
  integer, parameter       :: nb=4
-  real(kind(1.d0)),pointer :: iwork(:), xp(:), yp(:), id(:)
+  real(psb_dpk_),pointer :: iwork(:), xp(:), yp(:), id(:)
  character                :: itrans
  character(len=20)        :: name, ch_err
  logical                  :: aliw
--- a/base/psblas/psb_zamax.f90
+++ b/base/psblas/psb_zamax.f90
@ -52,16 +52,16 @@ function psb_zamax (x,desc_a, info, jx)
  use psb_error_mod
  implicit none
-  complex(kind(1.d0)), intent(in)      :: x(:,:)
+  complex(psb_dpk_), intent(in)      :: x(:,:)
  type(psb_desc_type), intent(in)   :: desc_a
  integer, intent(out)              :: info
  integer, optional, intent(in)     :: jx
-  real(kind(1.d0))                  :: psb_zamax
+  real(psb_dpk_)                  :: psb_zamax
  ! locals
  integer                  :: ictxt, np, me,&
       & err_act, iix, jjx, ix, ijx, m, imax, izamax
-  real(kind(1.d0))         :: amax
+  real(psb_dpk_)         :: amax
  character(len=20)        :: name, ch_err
  double complex    ::         zdum
  double precision  ::  cabs1
@ -183,16 +183,16 @@ function psb_zamaxv (x,desc_a, info)
  use psb_error_mod
  implicit none
-  complex(kind(1.d0)), intent(in)   :: x(:)
+  complex(psb_dpk_), intent(in)   :: x(:)
  type(psb_desc_type), intent(in)   :: desc_a
  integer, intent(out)              :: info
-  real(kind(1.d0))                  :: psb_zamaxv
+  real(psb_dpk_)                  :: psb_zamaxv
  ! locals
  integer                  :: ictxt, np, me,&
       & err_act, iix, jjx, jx, ix, m, imax, izamax
-  real(kind(1.d0))         :: amax
+  real(psb_dpk_)         :: amax
-  complex(kind(1.d0))      :: cmax
+  complex(psb_dpk_)      :: cmax
  character(len=20)        :: name, ch_err
  double complex    ::         zdum
  double precision  ::  cabs1
@ -312,17 +312,17 @@ subroutine psb_zamaxvs(res,x,desc_a, info)
  use psb_error_mod
  implicit none
-  complex(kind(1.d0)), intent(in)   :: x(:)
+  complex(psb_dpk_), intent(in)   :: x(:)
  type(psb_desc_type), intent(in)   :: desc_a
  integer, intent(out)              :: info
-  real(kind(1.D0)), intent(out)     :: res
+  real(psb_dpk_), intent(out)     :: res
  ! locals
  integer                  :: ictxt, np, me,&
       & err_act, iix, jjx, ix, ijx, m, imax, izamax
-  real(kind(1.d0))         :: amax
+  real(psb_dpk_)         :: amax
  character(len=20)        :: name, ch_err
-  complex(kind(1.d0))      :: cmax
+  complex(psb_dpk_)      :: cmax
  double complex    ::         zdum
  double precision  ::  cabs1
  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
@ -438,18 +438,18 @@ subroutine psb_zmamaxs(res,x,desc_a, info,jx)
  use psb_error_mod
  implicit none
-  complex(kind(1.d0)), intent(in)   :: x(:,:)
+  complex(psb_dpk_), intent(in)   :: x(:,:)
  type(psb_desc_type), intent(in)   :: desc_a
  integer, intent(out)              :: info
  integer, optional, intent(in)     :: jx
-  real(kind(1.d0)), intent(out)     :: res(:)
+  real(psb_dpk_), intent(out)     :: res(:)
  ! locals
  integer                  :: ictxt, np, me,&
       & err_act, iix, jjx, ix, ijx, m, imax, i, k, izamax
-  real(kind(1.d0))         :: amax
+  real(psb_dpk_)         :: amax
  character(len=20)        :: name, ch_err
-  complex(kind(1.d0))      :: cmax
+  complex(psb_dpk_)      :: cmax
  double complex    ::         zdum
  double precision  ::  cabs1
  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
--- a/base/psblas/psb_zasum.f90
+++ b/base/psblas/psb_zasum.f90
@ -53,18 +53,18 @@ function psb_zasum (x,desc_a, info, jx)
  use psb_penv_mod
  implicit none
-  complex(kind(1.d0)), intent(in)   :: x(:,:)
+  complex(psb_dpk_), intent(in)   :: x(:,:)
  type(psb_desc_type), intent(in)   :: desc_a
  integer, intent(out)              :: info
  integer, optional, intent(in)     :: jx
-  real(kind(1.d0))                  :: psb_zasum
+  real(psb_dpk_)                  :: psb_zasum
  ! locals
  integer                  :: ictxt, np, me, &
       & err_act, iix, jjx, ix, ijx, m, i, idx, ndm
-  real(kind(1.d0))         :: asum, dzasum
+  real(psb_dpk_)         :: asum, dzasum
  character(len=20)        :: name, ch_err
-  complex(kind(1.d0))      :: cmax
+  complex(psb_dpk_)      :: cmax
  double complex    ::         zdum
  double precision  ::  cabs1
  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
@ -201,17 +201,17 @@ function psb_zasumv(x,desc_a, info)
  use psb_penv_mod
  implicit none
-  complex(kind(1.d0)), intent(in)   :: x(:)
+  complex(psb_dpk_), intent(in)   :: x(:)
  type(psb_desc_type), intent(in)   :: desc_a
  integer, intent(out)              :: info
-  real(kind(1.d0))                  :: psb_zasumv
+  real(psb_dpk_)                  :: psb_zasumv
  ! locals
  integer                  :: ictxt, np, me,&
       & err_act, iix, jjx, jx, ix, m, i, idx, ndm
-  real(kind(1.d0))         :: asum, dzasum
+  real(psb_dpk_)         :: asum, dzasum
  character(len=20)        :: name, ch_err
-  complex(kind(1.d0))      :: cmax
+  complex(psb_dpk_)      :: cmax
  double complex    ::         zdum
  double precision  ::  cabs1
  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
@ -344,17 +344,17 @@ subroutine psb_zasumvs(res,x,desc_a, info)
  use psb_penv_mod
  implicit none
-  complex(kind(1.d0)), intent(in)   :: x(:)
+  complex(psb_dpk_), intent(in)   :: x(:)
-  real(kind(1.d0)), intent(out)     :: res
+  real(psb_dpk_), intent(out)     :: res
  type(psb_desc_type), intent(in)   :: desc_a
  integer, intent(out)              :: info
  ! locals
  integer                  :: ictxt, np, me,&
       & err_act, iix, jjx, ix, jx, m, i, idx, ndm
-  real(kind(1.d0))         :: asum, dzasum
+  real(psb_dpk_)         :: asum, dzasum
  character(len=20)        :: name, ch_err
-  complex(kind(1.d0))      :: cmax
+  complex(psb_dpk_)      :: cmax
  double complex    ::         zdum
  double precision  ::  cabs1
  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
--- a/base/psblas/psb_zaxpby.f90
+++ b/base/psblas/psb_zaxpby.f90
@ -61,9 +61,9 @@ subroutine  psb_zaxpby(alpha, x, beta,y,desc_a,info, n, jx, jy)
  integer, intent(in), optional   :: n, jx, jy
  integer, intent(out)            :: info
  type(psb_desc_type), intent(in) :: desc_a
-  complex(kind(1.D0)), intent(in)    :: alpha, beta
+  complex(psb_dpk_), intent(in)    :: alpha, beta
-  complex(kind(1.D0)), intent(in)    :: x(:,:)
+  complex(psb_dpk_), intent(in)    :: x(:,:)
-  complex(kind(1.D0)), intent(inout) :: y(:,:)
+  complex(psb_dpk_), intent(inout) :: y(:,:)
  ! locals
  integer                  :: ictxt, np, me,&
@ -214,9 +214,9 @@ subroutine  psb_zaxpbyv(alpha, x, beta,y,desc_a,info)
  integer, intent(out)            :: info
  type(psb_desc_type), intent(in) :: desc_a
-  complex(kind(1.D0)), intent(in)    :: alpha, beta
+  complex(psb_dpk_), intent(in)    :: alpha, beta
-  complex(kind(1.D0)), intent(in)    :: x(:)
+  complex(psb_dpk_), intent(in)    :: x(:)
-  complex(kind(1.D0)), intent(inout) :: y(:)
+  complex(psb_dpk_), intent(inout) :: y(:)
  ! locals
  integer                  :: ictxt, np, me,&
--- a/base/psblas/psb_zdot.f90
+++ b/base/psblas/psb_zdot.f90
@ -55,17 +55,17 @@ function psb_zdot(x, y,desc_a, info, jx, jy)
  use psb_penv_mod
  implicit none
-  complex(kind(1.d0)), intent(in)     :: x(:,:), y(:,:)
+  complex(psb_dpk_), intent(in)     :: x(:,:), y(:,:)
  type(psb_desc_type), intent(in)  :: desc_a
  integer, intent(in), optional    :: jx, jy
  integer, intent(out)             :: info
-  complex(kind(1.D0))              :: psb_zdot
+  complex(psb_dpk_)              :: psb_zdot
  ! locals
  integer                  :: ictxt, np, me, idx, ndm,&
       & err_act, iix, jjx, ix, ijx, iy, ijy, iiy, jjy, i, m
-  complex(kind(1.D0))         :: dot_local
+  complex(psb_dpk_)         :: dot_local
-  complex(kind(1.d0))         :: zdotc
+  complex(psb_dpk_)         :: zdotc
  character(len=20)        :: name, ch_err
  name='psb_zdot'
@ -208,16 +208,16 @@ function psb_zdotv(x, y,desc_a, info)
  use psb_penv_mod
  implicit none
-  complex(kind(1.d0)), intent(in)  :: x(:), y(:)
+  complex(psb_dpk_), intent(in)  :: x(:), y(:)
  type(psb_desc_type), intent(in)  :: desc_a
  integer, intent(out)             :: info
-  complex(kind(1.D0))              :: psb_zdotv
+  complex(psb_dpk_)              :: psb_zdotv
  ! locals
  integer                  :: ictxt, np, me, idx, ndm,&
       & err_act, iix, jjx, ix, jx, iy, jy, iiy, jjy, i, m
-  complex(kind(1.D0))         :: dot_local
+  complex(psb_dpk_)         :: dot_local
-  complex(kind(1.d0))         :: zdotc
+  complex(psb_dpk_)         :: zdotc
  character(len=20)        :: name, ch_err
  name='psb_zdot'
@ -345,16 +345,16 @@ subroutine psb_zdotvs(res, x, y,desc_a, info)
  use psb_penv_mod
  implicit none
-  complex(kind(1.d0)), intent(in)     :: x(:), y(:)
+  complex(psb_dpk_), intent(in)     :: x(:), y(:)
-  complex(kind(1.d0)), intent(out)    :: res
+  complex(psb_dpk_), intent(out)    :: res
  type(psb_desc_type), intent(in)  :: desc_a
  integer, intent(out)             :: info
  ! locals
  integer                  :: ictxt, np, me, idx, ndm,&
       & err_act, iix, jjx, ix, iy, iiy, jjy, i, m
-  complex(kind(1.D0))         :: dot_local
+  complex(psb_dpk_)         :: dot_local
-  complex(kind(1.d0))         :: zdotc
+  complex(psb_dpk_)         :: zdotc
  character(len=20)        :: name, ch_err
  name='psb_zdot'
@ -480,16 +480,16 @@ subroutine psb_zmdots(res, x, y, desc_a, info)
  use psb_penv_mod
  implicit none
-  complex(kind(1.d0)), intent(in)     :: x(:,:), y(:,:)
+  complex(psb_dpk_), intent(in)     :: x(:,:), y(:,:)
-  complex(kind(1.d0)), intent(out)    :: res(:)
+  complex(psb_dpk_), intent(out)    :: res(:)
  type(psb_desc_type), intent(in)  :: desc_a
  integer, intent(out)             :: info
  ! locals
  integer                  :: ictxt, np, me, idx, ndm,&
       & err_act, iix, jjx, ix, iy, iiy, jjy, i, m, j, k
-  complex(kind(1.d0)),allocatable  :: dot_local(:)
+  complex(psb_dpk_),allocatable  :: dot_local(:)
-  complex(kind(1.d0))         :: zdotc
+  complex(psb_dpk_)         :: zdotc
  character(len=20)        :: name, ch_err
  name='psb_zmdots'
--- a/base/psblas/psb_znrm2.f90
+++ b/base/psblas/psb_znrm2.f90
@ -51,16 +51,16 @@ function psb_znrm2(x, desc_a, info, jx)
  use psb_penv_mod
  implicit none
-  complex(kind(1.d0)), intent(in)      ::  x(:,:)
+  complex(psb_dpk_), intent(in)      ::  x(:,:)
  type(psb_desc_type), intent(in)   :: desc_a
  integer, intent(in), optional     :: jx
  integer, intent(out)              :: info
-  real(kind(1.D0))                  :: psb_znrm2
+  real(psb_dpk_)                  :: psb_znrm2
  ! locals
  integer                  :: ictxt, np, me,&
       & err_act, iix, jjx, ndim, ix, ijx, i, m, id, idx, ndm
-  real(kind(1.d0))         :: nrm2, dznrm2, dd
+  real(psb_dpk_)         :: nrm2, dznrm2, dd
  external dcombnrm2
  character(len=20)        :: name, ch_err
@ -187,15 +187,15 @@ function psb_znrm2v(x, desc_a, info)
  use psb_penv_mod
  implicit none
-  complex(kind(1.d0)), intent(in)   :: x(:)
+  complex(psb_dpk_), intent(in)   :: x(:)
  type(psb_desc_type), intent(in)   :: desc_a
  integer, intent(out)              :: info
-  real(kind(1.D0))                  :: psb_znrm2v
+  real(psb_dpk_)                  :: psb_znrm2v
  ! locals
  integer                  :: ictxt, np, me,&
       & err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm
-  real(kind(1.d0))         :: nrm2, dznrm2, dd
+  real(psb_dpk_)         :: nrm2, dznrm2, dd
  external dcombnrm2
  character(len=20)        :: name, ch_err
@ -319,15 +319,15 @@ subroutine psb_znrm2vs(res, x, desc_a, info)
  use psb_penv_mod
  implicit none
-  complex(kind(1.d0)), intent(in)   :: x(:)
+  complex(psb_dpk_), intent(in)   :: x(:)
-  real(kind(1.d0)), intent(out)     :: res
+  real(psb_dpk_), intent(out)     :: res
  type(psb_desc_type), intent(in)   :: desc_a
  integer, intent(out)              :: info
  ! locals
  integer                  :: ictxt, np, me,&
       & err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm
-  real(kind(1.d0))         :: nrm2, dznrm2, dd
+  real(psb_dpk_)         :: nrm2, dznrm2, dd
  external dcombnrm2
  character(len=20)        :: name, ch_err
--- a/base/psblas/psb_znrmi.f90
+++ b/base/psblas/psb_znrmi.f90
@ -52,12 +52,12 @@ function psb_znrmi(a,desc_a,info)
  type(psb_zspmat_type), intent(in)   :: a
  integer, intent(out)                :: info
  type(psb_desc_type), intent(in)     :: desc_a
-  real(kind(1.d0))                    :: psb_znrmi
+  real(psb_dpk_)                    :: psb_znrmi
  ! locals
  integer                  :: ictxt, np, me,&
       & err_act, n, iia, jja, ia, ja, mdim, ndim, m
-  real(kind(1.d0))         :: nrmi, zcsnmi
+  real(psb_dpk_)         :: nrmi, zcsnmi
  character(len=20)        :: name, ch_err
  name='psb_znrmi'
--- a/base/psblas/psb_zspmm.f90
+++ b/base/psblas/psb_zspmm.f90
@ -76,13 +76,13 @@ subroutine  psb_zspmm(alpha,a,x,beta,y,desc_a,info,&
  use psb_penv_mod
  implicit none
-  complex(kind(1.D0)), intent(in)             :: alpha, beta
+  complex(psb_dpk_), intent(in)             :: alpha, beta
-  complex(kind(1.d0)), intent(inout), target  :: x(:,:)
+  complex(psb_dpk_), intent(inout), target  :: x(:,:)
-  complex(kind(1.d0)), intent(inout), target  :: y(:,:)
+  complex(psb_dpk_), intent(inout), target  :: y(:,:)
  type(psb_zspmat_type), intent(in)        :: a
  type(psb_desc_type), intent(in)          :: desc_a
  integer, intent(out)                     :: info
-  complex(kind(1.d0)), optional, target      :: work(:)
+  complex(psb_dpk_), optional, target      :: work(:)
  character, intent(in), optional          :: trans
  integer, intent(in), optional            :: k, jx, jy
  logical, intent(in), optional            :: doswap
@ -93,8 +93,8 @@ subroutine  psb_zspmm(alpha,a,x,beta,y,desc_a,info,&
       & m, nrow, ncol, lldx, lldy, liwork, iiy, jjy,&
       & i, ib, ib1, ip, idx
  integer, parameter               :: nb=4
-  complex(kind(1.d0)), pointer     :: xp(:,:), yp(:,:), iwork(:)
+  complex(psb_dpk_), pointer     :: xp(:,:), yp(:,:), iwork(:)
-  complex(kind(1.d0)), allocatable :: xvsave(:,:)
+  complex(psb_dpk_), allocatable :: xvsave(:,:)
  character                        :: trans_
  character(len=20)                :: name, ch_err
  logical                          :: aliw, doswap_
@ -421,13 +421,13 @@ subroutine  psb_zspmv(alpha,a,x,beta,y,desc_a,info,&
  use psb_penv_mod
  implicit none
-  complex(kind(1.D0)), intent(in)             :: alpha, beta
+  complex(psb_dpk_), intent(in)             :: alpha, beta
-  complex(kind(1.d0)), intent(inout), target  :: x(:)
+  complex(psb_dpk_), intent(inout), target  :: x(:)
-  complex(kind(1.d0)), intent(inout), target  :: y(:)
+  complex(psb_dpk_), intent(inout), target  :: y(:)
  type(psb_zspmat_type), intent(in)        :: a
  type(psb_desc_type), intent(in)          :: desc_a
  integer, intent(out)                     :: info
-  complex(kind(1.d0)), optional, target    :: work(:)
+  complex(psb_dpk_), optional, target    :: work(:)
  character, intent(in), optional          :: trans
  logical, intent(in), optional            :: doswap
@ -437,8 +437,8 @@ subroutine  psb_zspmv(alpha,a,x,beta,y,desc_a,info,&
       & m, nrow, ncol, lldx, lldy, liwork, jx, jy, iiy, jjy,&
       & ib, ip, idx
  integer, parameter           :: nb=4
-  complex(kind(1.d0)), pointer :: iwork(:), xp(:), yp(:)
+  complex(psb_dpk_), pointer :: iwork(:), xp(:), yp(:)
-  complex(kind(1.d0)), allocatable :: xvsave(:)  
+  complex(psb_dpk_), allocatable :: xvsave(:)  
  character                    :: trans_
  character(len=20)            :: name, ch_err
  logical                      :: aliw, doswap_
--- a/base/psblas/psb_zspsm.f90
+++ b/base/psblas/psb_zspsm.f90
@ -87,14 +87,14 @@ subroutine  psb_zspsm(alpha,a,x,beta,y,desc_a,info,&
  use psb_penv_mod
  implicit none
-  complex(kind(1.D0)), intent(in)              :: alpha, beta
+  complex(psb_dpk_), intent(in)              :: alpha, beta
-  complex(kind(1.d0)), intent(in), target      :: x(:,:)
+  complex(psb_dpk_), intent(in), target      :: x(:,:)
-  complex(kind(1.d0)), intent(inout), target   :: y(:,:)
+  complex(psb_dpk_), intent(inout), target   :: y(:,:)
  type (psb_zspmat_type), intent(in)        :: a
  type(psb_desc_type), intent(in)           :: desc_a
  integer, intent(out)                      :: info
-  complex(kind(1.d0)), intent(in), optional, target      :: diag(:)
+  complex(psb_dpk_), intent(in), optional, target      :: diag(:)
-  complex(kind(1.d0)), optional, target       :: work(:)
+  complex(psb_dpk_), optional, target       :: work(:)
  character, intent(in), optional           :: trans, unitd
  integer, intent(in), optional             :: choice
  integer, intent(in), optional             :: k, jx, jy
@ -107,7 +107,7 @@ subroutine  psb_zspsm(alpha,a,x,beta,y,desc_a,info,&
  character                :: lunitd
  integer, parameter       :: nb=4
-  complex(kind(1.d0)),pointer :: iwork(:), xp(:,:), yp(:,:), id(:)
+  complex(psb_dpk_),pointer :: iwork(:), xp(:,:), yp(:,:), id(:)
  character                :: itrans
  character(len=20)        :: name, ch_err
  logical                  :: aliw
@ -375,14 +375,14 @@ subroutine  psb_zspsv(alpha,a,x,beta,y,desc_a,info,&
  use psb_penv_mod
  implicit none 
-  complex(kind(1.D0)), intent(in)              :: alpha, beta
+  complex(psb_dpk_), intent(in)              :: alpha, beta
-  complex(kind(1.d0)), intent(in), target      :: x(:)
+  complex(psb_dpk_), intent(in), target      :: x(:)
-  complex(kind(1.d0)), intent(inout), target   :: y(:)
+  complex(psb_dpk_), intent(inout), target   :: y(:)
  type(psb_zspmat_type), intent(in)         :: a
  type(psb_desc_type), intent(in)           :: desc_a
  integer, intent(out)                      :: info
-  complex(kind(1.d0)), intent(in), optional, target    :: diag(:)
+  complex(psb_dpk_), intent(in), optional, target    :: diag(:)
-  complex(kind(1.d0)), optional, target      :: work(:)
+  complex(psb_dpk_), optional, target      :: work(:)
  character, intent(in), optional           :: trans, unitd
  integer, intent(in), optional             :: choice
@ -394,7 +394,7 @@ subroutine  psb_zspsv(alpha,a,x,beta,y,desc_a,info,&
  character                :: lunitd
  integer, parameter       :: nb=4
-  complex(kind(1.d0)),pointer :: iwork(:), xp(:), yp(:), id(:)
+  complex(psb_dpk_),pointer :: iwork(:), xp(:), yp(:), id(:)
  character                :: itrans
  character(len=20)        :: name, ch_err
  logical                  :: aliw
--- a/base/serial/aux/dasr.f90
+++ b/base/serial/aux/dasr.f90
@ -39,10 +39,10 @@ subroutine dasr(n,x,dir)
  !
  !     .. Scalar Arguments ..
  integer, intent(in) :: n, dir 
-  real(kind(1.d0)) ::  x(n)
+  real(psb_dpk_) ::  x(n)
  !     ..
  !     .. Local Scalars ..
-  real(kind(1.d0)) :: piv, xt, xk
+  real(psb_dpk_) :: piv, xt, xk
  integer i, j, ilx, iux, istp, lpiv
  integer n1, n2
@ -314,9 +314,9 @@ contains
  subroutine disr_up(n,x)
    implicit none
    integer :: n
-    real(kind(1.d0)) :: x(n)
+    real(psb_dpk_) :: x(n)
    integer :: i,j
-    real(kind(1.d0)) :: xx,xax
+    real(psb_dpk_) :: xx,xax
    do j=n-1,1,-1
      if (abs(x(j+1)) < abs(x(j))) then
@ -337,9 +337,9 @@ contains
  subroutine disr_dw(n,x)
    implicit none
    integer :: n
-    real(kind(1.d0)) :: x(n)
+    real(psb_dpk_) :: x(n)
    integer :: i,j
-    real(kind(1.d0)) :: xx,xax
+    real(psb_dpk_) :: xx,xax
    do j=n-1,1,-1
      if (abs(x(j+1)) > abs(x(j))) then
--- a/base/serial/aux/dasrx.f90
+++ b/base/serial/aux/dasrx.f90
@ -38,11 +38,11 @@ subroutine dasrx(n,x,indx,dir,flag)
  !
  !     .. Scalar Arguments ..
  integer, intent(in) :: n, dir, flag
-  real(kind(1.d0)) ::  x(n)
+  real(psb_dpk_) ::  x(n)
  integer ::  indx(n)
  !     ..
  !     .. Local Scalars ..
-  real(kind(1.d0)) :: piv, xt, xk
+  real(psb_dpk_) :: piv, xt, xk
  integer i, j, ilx, iux, istp, lpiv
  integer ixt, n1, n2
@ -354,10 +354,10 @@ contains
  subroutine idasrx_up(n,x,indx)
    implicit none
    integer :: n
-    real(kind(1.d0)) :: x(n)
+    real(psb_dpk_) :: x(n)
    integer :: indx(n)
    integer :: i,j,ix
-    real(kind(1.d0)) :: xx,xax
+    real(psb_dpk_) :: xx,xax
    do j=n-1,1,-1
      if (abs(x(j+1)) < abs(x(j))) then
@ -381,10 +381,10 @@ contains
  subroutine idasrx_dw(n,x,indx)
    implicit none
    integer :: n
-    real(kind(1.d0)) :: x(n)
+    real(psb_dpk_) :: x(n)
    integer :: indx(n)
    integer :: i,j,ix
-    real(kind(1.d0)) :: xx,xax
+    real(psb_dpk_) :: xx,xax
    do j=n-1,1,-1
      if (abs(x(j+1)) > abs(x(j))) then
--- a/base/serial/aux/dmsort_dw.f90
+++ b/base/serial/aux/dmsort_dw.f90
@ -68,9 +68,10 @@
 !
 !
 subroutine dmsort_dw(n,k,l,iret)
  use psb_const_mod
  implicit none
  integer n, iret
-  real(kind(1.d0)) :: k(n)
+  real(psb_dpk_) :: k(n)
  integer l(0:n+1)
  !
  integer p,q,s,t
--- a/base/serial/aux/dmsort_up.f90
+++ b/base/serial/aux/dmsort_up.f90
@ -68,9 +68,10 @@
 !
 !
 subroutine dmsort_up(n,k,l,iret)
  use psb_const_mod
  implicit none
  integer n, iret
-  real(kind(1.d0))  k(n)
+  real(psb_dpk_)  k(n)
  integer l(0:n+1)
  !
  integer p,q,s,t
--- a/base/serial/aux/dmsr.f90
+++ b/base/serial/aux/dmsr.f90
@ -37,13 +37,13 @@ subroutine dmsr(n,x,idir)
  implicit none
  integer :: n, idir
-  real(kind(1.d0)) :: x(n)
+  real(psb_dpk_) :: x(n)
  integer, allocatable :: iaux(:)
  integer :: lswap, iret, info, lp, k
-  real(kind(1.d0)) :: swap
+  real(psb_dpk_) :: swap
  if (n<0) then 
 !!$    write(0,*) 'Error: IMSR: N<0'
--- a/base/serial/aux/dmsrx.f90
+++ b/base/serial/aux/dmsrx.f90
@ -36,13 +36,13 @@ subroutine dmsrx(n,x,indx,idir,flag)
  use psb_serial_mod
  implicit none
  integer :: n,idir,flag
-  real(kind(1.d0)) :: x(n)
+  real(psb_dpk_) :: x(n)
  integer :: indx(n)
  integer, allocatable :: iaux(:)
  integer :: iret, info, lp, k,lswap, ixswap
-  real(kind(1.d0)) :: swap
+  real(psb_dpk_) :: swap
  if (n<0) then 
    write(0,*) 'Error: DMSRX: N<0'
--- a/base/serial/aux/dsr.f90
+++ b/base/serial/aux/dsr.f90
@ -39,10 +39,10 @@ subroutine dsr(n,x,dir)
  !
  !     .. Scalar Arguments ..
  integer, intent(in) :: n, dir 
-  real(kind(1.d0)) ::  x(n)
+  real(psb_dpk_) ::  x(n)
  !     ..
  !     .. Local Scalars ..
-  real(kind(1.d0)) :: piv, xt, xk
+  real(psb_dpk_) :: piv, xt, xk
  integer i, j, ilx, iux, istp, lpiv
  integer n1, n2
@ -314,9 +314,9 @@ contains
  subroutine disr_up(n,x)
    implicit none
    integer :: n
-    real(kind(1.d0)) :: x(n)
+    real(psb_dpk_) :: x(n)
    integer :: i,j
-    real(kind(1.d0)) :: xx
+    real(psb_dpk_) :: xx
    do j=n-1,1,-1
      if (x(j+1) < x(j)) then
@ -336,9 +336,9 @@ contains
  subroutine disr_dw(n,x)
    implicit none
    integer :: n
-    real(kind(1.d0)) :: x(n)
+    real(psb_dpk_) :: x(n)
    integer :: i,j
-    real(kind(1.d0)) :: xx
+    real(psb_dpk_) :: xx
    do j=n-1,1,-1
      if (x(j+1) > x(j)) then
--- a/base/serial/aux/dsrx.f90
+++ b/base/serial/aux/dsrx.f90
@ -38,11 +38,11 @@ subroutine dsrx(n,x,indx,dir,flag)
  !
  !     .. Scalar Arguments ..
  integer, intent(in) :: n, dir, flag
-  real(kind(1.d0)) ::  x(n)
+  real(psb_dpk_) ::  x(n)
  integer ::  indx(n)
  !     ..
  !     .. Local Scalars ..
-  real(kind(1.d0)) :: piv, xk, xt
+  real(psb_dpk_) :: piv, xk, xt
  integer i, j, ilx, iux, istp, lpiv
  integer ixt, n1, n2
@ -356,10 +356,10 @@ contains
  subroutine idsrx_up(n,x,indx)
    implicit none
    integer :: n
-    real(kind(1.d0)) :: x(n)
+    real(psb_dpk_) :: x(n)
    integer :: indx(n)
    integer :: i,j,ix
-    real(kind(1.d0)) :: xx
+    real(psb_dpk_) :: xx
    do j=n-1,1,-1
      if (x(j+1) < x(j)) then
@ -382,10 +382,10 @@ contains
  subroutine idsrx_dw(n,x,indx)
    implicit none
    integer :: n
-    real(kind(1.d0)) :: x(n)
+    real(psb_dpk_) :: x(n)
    integer :: indx(n)
    integer :: i,j,ix
-    real(kind(1.d0)) :: xx
+    real(psb_dpk_) :: xx
    do j=n-1,1,-1
      if (x(j+1) > x(j)) then
--- a/base/serial/aux/zacmp_mod.f90
+++ b/base/serial/aux/zacmp_mod.f90
@ -30,6 +30,7 @@
 !!$ 
 !!$  
 module zacmp_mod
  use psb_const_mod
  interface operator(<)
    module procedure zalt
  end interface
@ -46,26 +47,30 @@ module zacmp_mod
 contains
  function zalt(a,b)
-    complex(kind(1.d0)), intent(in) :: a,b
+    use psb_const_mod
    complex(psb_dpk_), intent(in) :: a,b
    logical :: zalt
    zalt = (abs(a) < abs(b))
  end function zalt
  function zale(a,b)
-    complex(kind(1.d0)), intent(in) :: a,b
+    use psb_const_mod
    complex(psb_dpk_), intent(in) :: a,b
    logical :: zale
    zale = (abs(a) <= abs(b))
  end function zale
  function zagt(a,b)
-    complex(kind(1.d0)), intent(in) :: a,b
+    use psb_const_mod
    complex(psb_dpk_), intent(in) :: a,b
    logical :: zagt
    zagt = (abs(a) > abs(b))
  end function zagt
  function zage(a,b)
-    complex(kind(1.d0)), intent(in) :: a,b
+    use psb_const_mod
    complex(psb_dpk_), intent(in) :: a,b
    logical :: zage
    zage = (abs(a) >= abs(b))
--- a/base/serial/aux/zalcmp_mod.f90
+++ b/base/serial/aux/zalcmp_mod.f90
@ -30,6 +30,7 @@
 !!$ 
 !!$  
 module zalcmp_mod
  use psb_const_mod
  interface operator(<)
    module procedure zallt
  end interface
@ -46,14 +47,16 @@ module zalcmp_mod
 contains
  function zallt(a,b)
-    complex(kind(1.d0)), intent(in) :: a,b
+    use psb_const_mod
    complex(psb_dpk_), intent(in) :: a,b
    logical :: zallt
    zallt = (abs(real(a))<abs(real(b))).or. &
         & ((abs(real(a))==abs(real(b))).and.(abs(aimag(a))<abs(aimag(b))))
  end function zallt
  function zalle(a,b)
-    complex(kind(1.d0)), intent(in) :: a,b
+    use psb_const_mod
    complex(psb_dpk_), intent(in) :: a,b
    logical :: zalle
    zalle = (abs(real(a))<abs(real(b))).or. &
@ -61,14 +64,16 @@ contains
  end function zalle
  function zalgt(a,b)
-    complex(kind(1.d0)), intent(in) :: a,b
+    use psb_const_mod
    complex(psb_dpk_), intent(in) :: a,b
    logical :: zalgt
    zalgt = (abs(real(a))>abs(real(b))).or. &
         & ((abs(real(a))==abs(real(b))).and.(abs(aimag(a))>abs(aimag(b))))
  end function zalgt
  function zalge(a,b)
-    complex(kind(1.d0)), intent(in) :: a,b
+    use psb_const_mod
    complex(psb_dpk_), intent(in) :: a,b
    logical :: zalge
    zalge = (abs(real(a))>abs(real(b))).or. &
--- a/base/serial/aux/zalsr.f90
+++ b/base/serial/aux/zalsr.f90
@ -40,10 +40,10 @@ subroutine zalsr(n,x,dir)
  !
  !     .. Scalar Arguments ..
  integer, intent(in) :: n, dir 
-  complex(kind(1.d0)) ::  x(n)
+  complex(psb_dpk_) ::  x(n)
  !     ..
  !     .. Local Scalars ..
-  complex(kind(1.d0)) :: xk, piv, xt
+  complex(psb_dpk_) :: xk, piv, xt
  integer i, j, ilx, iux, istp, lpiv
  integer n1, n2
@ -316,9 +316,9 @@ contains
    use zalcmp_mod
    implicit none
    integer :: n
-    complex(kind(1.d0)) :: x(n)
+    complex(psb_dpk_) :: x(n)
    integer :: i,j
-    complex(kind(1.d0)) :: xx
+    complex(psb_dpk_) :: xx
    do j=n-1,1,-1
      if (x(j+1) < x(j)) then
@ -339,9 +339,9 @@ contains
    use zalcmp_mod
    implicit none
    integer :: n
-    complex(kind(1.d0)) :: x(n)
+    complex(psb_dpk_) :: x(n)
    integer :: i,j
-    complex(kind(1.d0)) :: xx
+    complex(psb_dpk_) :: xx
    do j=n-1,1,-1
      if (x(j+1) > x(j)) then
--- a/base/serial/aux/zalsrx.f90
+++ b/base/serial/aux/zalsrx.f90
@ -39,11 +39,11 @@ subroutine zalsrx(n,x,indx,dir,flag)
  !
  !     .. Scalar Arguments ..
  integer, intent(in) :: n, dir, flag
-  complex(kind(1.d0)) ::  x(n)
+  complex(psb_dpk_) ::  x(n)
  integer ::  indx(n)
  !     ..
  !     .. Local Scalars ..
-  complex(kind(1.d0)) :: piv, xk, xt
+  complex(psb_dpk_) :: piv, xk, xt
  integer i, j, ilx, iux, istp, lpiv
  integer ixt, n1, n2
@ -359,10 +359,10 @@ contains
    use zalcmp_mod
    implicit none
    integer :: n
-    complex(kind(1.d0)) :: x(n)
+    complex(psb_dpk_) :: x(n)
    integer :: indx(n)
    integer :: i,j,ix
-    complex(kind(1.d0)) :: xx
+    complex(psb_dpk_) :: xx
    do j=n-1,1,-1
      if (x(j+1) < x(j)) then
@ -386,10 +386,10 @@ contains
    use zalcmp_mod
    implicit none
    integer :: n
-    complex(kind(1.d0)) :: x(n)
+    complex(psb_dpk_) :: x(n)
    integer :: indx(n)
    integer :: i,j,ix
-    complex(kind(1.d0)) :: xx
+    complex(psb_dpk_) :: xx
    do j=n-1,1,-1
      if (x(j+1) > x(j)) then
--- a/base/serial/aux/zamsort_dw.f90
+++ b/base/serial/aux/zamsort_dw.f90
@ -71,7 +71,7 @@ subroutine zamsort_dw(n,k,l,iret)
  use zacmp_mod
  implicit none
  integer n, iret
-  complex(kind(1.d0)) :: k(n)
+  complex(psb_dpk_) :: k(n)
  integer l(0:n+1)
  !
  integer p,q,s,t
--- a/base/serial/aux/zamsort_up.f90
+++ b/base/serial/aux/zamsort_up.f90
@ -71,7 +71,7 @@ subroutine zamsort_up(n,k,l,iret)
  use zacmp_mod
  implicit none
  integer n, iret
-  complex(kind(1.d0))  k(n)
+  complex(psb_dpk_)  k(n)
  integer l(0:n+1)
  !
  integer p,q,s,t
--- a/base/serial/aux/zamsr.f90
+++ b/base/serial/aux/zamsr.f90
@ -37,13 +37,13 @@ subroutine zamsr(n,x,idir)
  implicit none
  integer :: n, idir
-  complex(kind(1.d0)) :: x(n)
+  complex(psb_dpk_) :: x(n)
  integer, allocatable :: iaux(:)
  integer :: lswap, iret, info, lp, k
-  complex(kind(1.d0)) :: swap
+  complex(psb_dpk_) :: swap
  if (n<0) then 
 !!$    write(0,*) 'Error: IMSR: N<0'
--- a/base/serial/aux/zamsrx.f90
+++ b/base/serial/aux/zamsrx.f90
@ -36,13 +36,13 @@ subroutine zamsrx(n,x,indx,idir,flag)
  use psb_serial_mod
  implicit none
  integer :: n,idir,flag
-  complex(kind(1.d0)) :: x(n)
+  complex(psb_dpk_) :: x(n)
  integer :: indx(n)
  integer, allocatable :: iaux(:)
  integer :: iret, info, lp, k,lswap, ixswap
-  complex(kind(1.d0)) :: swap
+  complex(psb_dpk_) :: swap
  if (n<0) then 
    write(0,*) 'Error: DMSRX: N<0'
--- a/base/serial/aux/zasr.f90
+++ b/base/serial/aux/zasr.f90
@ -40,10 +40,10 @@ subroutine zasr(n,x,dir)
  !
  !     .. Scalar Arguments ..
  integer, intent(in) :: n, dir 
-  complex(kind(1.d0)) ::  x(n)
+  complex(psb_dpk_) ::  x(n)
  !     ..
  !     .. Local Scalars ..
-  complex(kind(1.d0)) :: xk, piv, xt
+  complex(psb_dpk_) :: xk, piv, xt
  integer i, j, ilx, iux, istp, lpiv
  integer n1, n2
@ -316,9 +316,9 @@ contains
    use zacmp_mod
    implicit none
    integer :: n
-    complex(kind(1.d0)) :: x(n)
+    complex(psb_dpk_) :: x(n)
    integer :: i,j
-    complex(kind(1.d0)) :: xx
+    complex(psb_dpk_) :: xx
    do j=n-1,1,-1
      if (x(j+1) < x(j)) then
@ -339,9 +339,9 @@ contains
    use zacmp_mod
    implicit none
    integer :: n
-    complex(kind(1.d0)) :: x(n)
+    complex(psb_dpk_) :: x(n)
    integer :: i,j
-    complex(kind(1.d0)) :: xx
+    complex(psb_dpk_) :: xx
    do j=n-1,1,-1
      if (x(j+1) > x(j)) then
--- a/base/serial/aux/zasrx.f90
+++ b/base/serial/aux/zasrx.f90
@ -39,11 +39,11 @@ subroutine zasrx(n,x,indx,dir,flag)
  !
  !     .. Scalar Arguments ..
  integer, intent(in) :: n, dir, flag
-  complex(kind(1.d0)) ::  x(n)
+  complex(psb_dpk_) ::  x(n)
  integer ::  indx(n)
  !     ..
  !     .. Local Scalars ..
-  complex(kind(1.d0)) :: piv, xk, xt
+  complex(psb_dpk_) :: piv, xk, xt
  integer i, j, ilx, iux, istp, lpiv
  integer ixt, n1, n2
@ -359,10 +359,10 @@ contains
    use zacmp_mod
    implicit none
    integer :: n
-    complex(kind(1.d0)) :: x(n)
+    complex(psb_dpk_) :: x(n)
    integer :: indx(n)
    integer :: i,j,ix
-    complex(kind(1.d0)) :: xx
+    complex(psb_dpk_) :: xx
    do j=n-1,1,-1
      if (x(j+1) < x(j)) then
@ -386,10 +386,10 @@ contains
    use zacmp_mod
    implicit none
    integer :: n
-    complex(kind(1.d0)) :: x(n)
+    complex(psb_dpk_) :: x(n)
    integer :: indx(n)
    integer :: i,j,ix
-    complex(kind(1.d0)) :: xx
+    complex(psb_dpk_) :: xx
    do j=n-1,1,-1
      if (x(j+1) > x(j)) then
--- a/base/serial/aux/zlcmp_mod.f90
+++ b/base/serial/aux/zlcmp_mod.f90
@ -30,6 +30,7 @@
 !!$ 
 !!$  
 module zlcmp_mod
  use psb_const_mod
  interface operator(<)
    module procedure zllt
  end interface
@ -46,26 +47,30 @@ module zlcmp_mod
 contains
  function zllt(a,b)
-    complex(kind(1.d0)), intent(in) :: a,b
+    use psb_const_mod
    complex(psb_dpk_), intent(in) :: a,b
    logical :: zllt
    zllt = (real(a)<real(b)).or.((real(a)==real(b)).and.(aimag(a)<aimag(b)))
  end function zllt
  function zlle(a,b)
-    complex(kind(1.d0)), intent(in) :: a,b
+    use psb_const_mod
    complex(psb_dpk_), intent(in) :: a,b
    logical :: zlle
    zlle = (real(a)<real(b)).or.((real(a)==real(b)).and.(aimag(a)<=aimag(b)))
  end function zlle
  function zlgt(a,b)
-    complex(kind(1.d0)), intent(in) :: a,b
+    use psb_const_mod
    complex(psb_dpk_), intent(in) :: a,b
    logical :: zlgt
    zlgt = (real(a)>real(b)).or.((real(a)==real(b)).and.(aimag(a)>aimag(b)))
  end function zlgt
  function zlge(a,b)
-    complex(kind(1.d0)), intent(in) :: a,b
+    use psb_const_mod
    complex(psb_dpk_), intent(in) :: a,b
    logical :: zlge
    zlge = (real(a)>real(b)).or.((real(a)==real(b)).and.(aimag(a)>=aimag(b)))
--- a/base/serial/aux/zlsr.f90
+++ b/base/serial/aux/zlsr.f90
@ -40,10 +40,10 @@ subroutine zlsr(n,x,dir)
  !
  !     .. Scalar Arguments ..
  integer, intent(in) :: n, dir 
-  complex(kind(1.d0)) ::  x(n)
+  complex(psb_dpk_) ::  x(n)
  !     ..
  !     .. Local Scalars ..
-  complex(kind(1.d0)) :: xk, piv, xt
+  complex(psb_dpk_) :: xk, piv, xt
  integer i, j, ilx, iux, istp, lpiv
  integer n1, n2
@ -316,9 +316,9 @@ contains
    use zlcmp_mod
    implicit none
    integer :: n
-    complex(kind(1.d0)) :: x(n)
+    complex(psb_dpk_) :: x(n)
    integer :: i,j
-    complex(kind(1.d0)) :: xx
+    complex(psb_dpk_) :: xx
    do j=n-1,1,-1
      if (x(j+1) < x(j)) then
@ -339,9 +339,9 @@ contains
    use zlcmp_mod
    implicit none
    integer :: n
-    complex(kind(1.d0)) :: x(n)
+    complex(psb_dpk_) :: x(n)
    integer :: i,j
-    complex(kind(1.d0)) :: xx
+    complex(psb_dpk_) :: xx
    do j=n-1,1,-1
      if (x(j+1) > x(j)) then
--- a/base/serial/aux/zlsrx.f90
+++ b/base/serial/aux/zlsrx.f90
@ -39,11 +39,11 @@ subroutine zlsrx(n,x,indx,dir,flag)
  !
  !     .. Scalar Arguments ..
  integer, intent(in) :: n, dir, flag
-  complex(kind(1.d0)) ::  x(n)
+  complex(psb_dpk_) ::  x(n)
  integer ::  indx(n)
  !     ..
  !     .. Local Scalars ..
-  complex(kind(1.d0)) :: piv, xk, xt
+  complex(psb_dpk_) :: piv, xk, xt
  integer i, j, ilx, iux, istp, lpiv
  integer ixt, n1, n2
@ -359,10 +359,10 @@ contains
    use zlcmp_mod
    implicit none
    integer :: n
-    complex(kind(1.d0)) :: x(n)
+    complex(psb_dpk_) :: x(n)
    integer :: indx(n)
    integer :: i,j,ix
-    complex(kind(1.d0)) :: xx
+    complex(psb_dpk_) :: xx
    do j=n-1,1,-1
      if (x(j+1) < x(j)) then
@ -386,10 +386,10 @@ contains
    use zlcmp_mod
    implicit none
    integer :: n
-    complex(kind(1.d0)) :: x(n)
+    complex(psb_dpk_) :: x(n)
    integer :: indx(n)
    integer :: i,j,ix
-    complex(kind(1.d0)) :: xx
+    complex(psb_dpk_) :: xx
    do j=n-1,1,-1
      if (x(j+1) > x(j)) then
--- a/base/serial/coo/dcoomm.f
+++ b/base/serial/coo/dcoomm.f
@ -37,14 +37,15 @@ c
 *
      SUBROUTINE DCOOMM(TRANSA,M,K,N,ALPHA,DESCRA,AR,
     *   IA,JA,INFOA,B,LDB,BETA,C,LDC,WORK,LWORK,IERROR)
      use psb_const_mod
 C
 C
 C     .. Scalar Arguments ..
-      DOUBLE PRECISION  ALPHA, BETA
+      real(psb_dpk_)  ALPHA, BETA
      INTEGER           K, LDB, LDC, M, N, LWORK, IERROR
      CHARACTER         TRANSA
 C     .. Array Arguments ..
-      DOUBLE PRECISION  AR(*), B(LDB,*), C(LDC,*),  WORK(*)
+      real(psb_dpk_)  AR(*), B(LDB,*), C(LDC,*),  WORK(*)
      INTEGER           IA(*), JA(*),INFOA(*)
      CHARACTER         DESCRA*11
 C     .. Local Scalars ..
--- a/base/serial/coo/dcoomv.f
+++ b/base/serial/coo/dcoomv.f
@ -167,18 +167,19 @@ C
 ***********************************************************************
      SUBROUTINE DCOOMV (TRANS,DIAG,M,N,ALPHA,AS,IA,JA,INFOA,X,
     +  BETA,Y,WORK,IERROR)
      use psb_const_mod
 C     .. Parameters ..
-      DOUBLE PRECISION  ONE, ZERO
+      real(psb_dpk_)  ONE, ZERO
      PARAMETER         (ONE=1.0D0,ZERO=0.0D0)
 C     .. Scalar Arguments ..
-      DOUBLE PRECISION  ALPHA, BETA
+      real(psb_dpk_)  ALPHA, BETA
      INTEGER           M, N, IERROR
      CHARACTER         DIAG, TRANS
 C     .. Array Arguments ..
-      DOUBLE PRECISION  AS(*), WORK(*), X(*), Y(*)
+      real(psb_dpk_)  AS(*), WORK(*), X(*), Y(*)
      INTEGER           IA(*), JA(*),infoa(*)
 C     .. Local Scalars ..
-      DOUBLE PRECISION  ACC,  TX
+      real(psb_dpk_)  ACC,  TX
      INTEGER           I, J, K, NNZ, IR, JC
      LOGICAL           SYM, TRA, UNI
 C     .. Executable Statements ..
--- a/base/serial/coo/dcoonrmi.f
+++ b/base/serial/coo/dcoonrmi.f
@ -30,8 +30,9 @@ C POSSIBILITY OF SUCH DAMAGE.
 C
 C 
 C     ... Compute Infinity norm for sparse matrix in CSR Format ...
-      DOUBLE PRECISION FUNCTION DCOONRMI(TRANS,M,N,DESCRA,A,IA1,IA2,
+      real(psb_dpk_) FUNCTION DCOONRMI(TRANS,M,N,DESCRA,A,IA1,IA2,
     +   INFOA,IERROR)
      use psb_const_mod
      IMPLICIT NONE
 C     .. Scalar Arguments ..
      INTEGER           M,N, IERROR
@ -39,10 +40,10 @@ C     .. Scalar Arguments ..
 C     .. Array Arguments ..
      INTEGER           IA1(*),IA2(*),INFOA(*)
      CHARACTER         DESCRA*11
-      DOUBLE PRECISION  A(*)
+      real(psb_dpk_)  A(*)
 C     .. Local scalars ..
      INTEGER I, J, K, NNZ
-      DOUBLE PRECISION NRMI, SUM
+      real(psb_dpk_) NRMI, SUM
      IERROR=0
      NRMI = 0.0
--- a/base/serial/coo/dcoorws.f
+++ b/base/serial/coo/dcoorws.f
@ -31,6 +31,7 @@ C
 C 
      SUBROUTINE DCOORWS(TRANS,M,N,DESCRA,A,IA1,IA2,
     &   INFOA,ROWSUM,IERROR)
      use psb_const_mod
      IMPLICIT NONE
 C     .. Scalar Arguments ..
      INTEGER           M,N, IERROR
@ -38,10 +39,10 @@ C     .. Scalar Arguments ..
 C     .. Array Arguments ..
      INTEGER           IA1(*),IA2(*),INFOA(*)
      CHARACTER         DESCRA*11
-      DOUBLE PRECISION  A(*), ROWSUM(*)
+      real(psb_dpk_)  A(*), ROWSUM(*)
 C     .. Local scalars ..
      INTEGER I, J, NNZ, K
-      DOUBLE PRECISION  SUM
+      real(psb_dpk_)  SUM
      logical psb_lsame
      external psb_lsame
--- a/base/serial/coo/dcoosm.f
+++ b/base/serial/coo/dcoosm.f
@ -32,11 +32,12 @@ C
      SUBROUTINE DCOOSM(TRANST,M,N,UNITD,D,ALPHA,DESCRA,A,IA,JA,INFOA,
     *                  B,LDB,BETA,C,LDC,WORK,LWORK,IERROR)
      use psb_error_mod
      use psb_const_mod
      IMPLICIT NONE
-      DOUBLE PRECISION  ALPHA, BETA
+      real(psb_dpk_)  ALPHA, BETA
      INTEGER           LDB, LDC, LWORK, M, N, IERROR
      CHARACTER         UNITD, TRANST
-      DOUBLE PRECISION  A(*), B(LDB,*), C(LDC,*), D(*), WORK(*)
+      real(psb_dpk_)  A(*), B(LDB,*), C(LDC,*), D(*), WORK(*)
      INTEGER           IA(*), JA(*), INFOA(*), INT_VAL(5)
      CHARACTER         DESCRA*11
      INTEGER           I, K, ERR_ACT
--- a/base/serial/coo/dcoosv.f
+++ b/base/serial/coo/dcoosv.f
@ -35,14 +35,15 @@ C  "right" place, i.e. the last in its row for Lower and the first
 C  for Upper. 
 C
      SUBROUTINE DCOOSV (UPLO,TRANS,DIAG,N,AS,IA,JA,INFOA,B,X,IERROR)
      use psb_const_mod
      use psb_error_mod
-      DOUBLE PRECISION   ZERO
+      real(psb_dpk_)   ZERO
      PARAMETER         (ZERO=0.0D0)
      INTEGER           N,IERROR
      CHARACTER         DIAG, TRANS, UPLO
-      DOUBLE PRECISION  AS(*), B(*), X(*)
+      real(psb_dpk_)  AS(*), B(*), X(*)
      INTEGER           IA(*), JA(*),INFOA(*)
-      DOUBLE PRECISION  ACC
+      real(psb_dpk_)  ACC
      INTEGER           I, J, K, NNZ, II
      LOGICAL           LOW, TRA, UNI
      integer           :: debug_level, debug_unit
--- a/base/serial/coo/zcoomm.f
+++ b/base/serial/coo/zcoomm.f
@ -36,14 +36,15 @@ c
 *
      SUBROUTINE ZCOOMM(TRANSA,M,K,N,ALPHA,DESCRA,AR,
     *   IA,JA,INFOA,B,LDB,BETA,C,LDC,WORK,LWORK)
      use psb_const_mod
 C
 C
 C     .. Scalar Arguments ..
-      COMPLEX*16        ALPHA, BETA
+      complex(psb_dpk_)        ALPHA, BETA
      INTEGER           K, LDB, LDC, M, N, LWORK
      CHARACTER         TRANSA
 C     .. Array Arguments ..
-      COMPLEX*16        AR(*), B(LDB,*), C(LDC,*),  WORK(*)
+      complex(psb_dpk_)        AR(*), B(LDB,*), C(LDC,*),  WORK(*)
      INTEGER           IA(*), JA(*),INFOA(*)
      CHARACTER         DESCRA*11
 C     .. Local Scalars ..
--- a/base/serial/coo/zcoomv.f
+++ b/base/serial/coo/zcoomv.f
@ -166,18 +166,19 @@ C
 ***********************************************************************
      SUBROUTINE ZCOOMV (TRANS,DIAG,M,N,ALPHA,AS,IA,JA,INFOA,X,
     +   BETA,Y,WORK)
      use psb_const_mod
 C     .. Parameters ..
-      COMPLEX*16         ONE, ZERO
+      complex(psb_dpk_)         ONE, ZERO
      PARAMETER         (ONE=(1.0D0,0.0D0),ZERO=(0.0D0,0.0D0))
 C     .. Scalar Arguments ..
-      COMPLEX*16       ALPHA, BETA
+      complex(psb_dpk_)       ALPHA, BETA
      INTEGER           M, N
      CHARACTER         DIAG, TRANS
 C     .. Array Arguments ..
-      COMPLEX*16        AS(*), WORK(*), X(*), Y(*)
+      complex(psb_dpk_)        AS(*), WORK(*), X(*), Y(*)
      INTEGER           IA(*), JA(*),infoa(*)
 C     .. Local Scalars ..
-      COMPLEX*16        ACC,  TX
+      complex(psb_dpk_)        ACC,  TX
      INTEGER           I, J, K, NCOLA, NROWA,DUM, NNZ, IR, JC
      LOGICAL           SYM, TRA, UNI
 C     .. Executable Statements ..
--- a/base/serial/coo/zcoonrmi.f
+++ b/base/serial/coo/zcoonrmi.f
@ -29,8 +29,9 @@ C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 C POSSIBILITY OF SUCH DAMAGE.
 C
 C     ... Compute Infinity norm for sparse matrix in CSR Format ...
-      DOUBLE PRECISION FUNCTION ZCOONRMI(TRANS,M,N,DESCRA,A,IA1,IA2,
+      real(psb_dpk_) FUNCTION ZCOONRMI(TRANS,M,N,DESCRA,A,IA1,IA2,
     +   INFOA,IERROR)
      use psb_const_mod
      IMPLICIT NONE
 C     .. Scalar Arguments ..
      INTEGER           M,N, IERROR
@ -38,10 +39,10 @@ C     .. Scalar Arguments ..
 C     .. Array Arguments ..
      INTEGER           IA1(*),IA2(*),INFOA(*)
      CHARACTER         DESCRA*11
-      COMPLEX*16        A(*)
+      complex(psb_dpk_)        A(*)
 C     .. Local scalars ..
      INTEGER I, J, K, NNZ
-      DOUBLE PRECISION NRMI, SUM
+      real(psb_dpk_) NRMI, SUM
      NRMI = 0.0
      NNZ  = INFOA(1)
--- a/base/serial/coo/zcoorws.f
+++ b/base/serial/coo/zcoorws.f
@ -30,6 +30,7 @@ C POSSIBILITY OF SUCH DAMAGE.
 C
      SUBROUTINE ZCOORWS(TRANS,M,N,DESCRA,A,IA1,IA2,
     &  INFOA,ROWSUM,IERROR)
      use psb_const_mod
      IMPLICIT NONE
 C     .. Scalar Arguments ..
      INTEGER           M,N, IERROR
@ -37,10 +38,10 @@ C     .. Scalar Arguments ..
 C     .. Array Arguments ..
      INTEGER           IA1(*),IA2(*),INFOA(*)
      CHARACTER         DESCRA*11
-      COMPLEX*16        A(*), ROWSUM(*)
+      complex(psb_dpk_)        A(*), ROWSUM(*)
 C     .. Local scalars ..
      INTEGER I, J, NNZ, K
-      DOUBLE PRECISION  SUM
+      real(psb_dpk_)  SUM
      logical psb_lsame
      external psb_lsame
--- a/base/serial/coo/zcoosm.f
+++ b/base/serial/coo/zcoosm.f
@ -31,12 +31,13 @@ C
 C 
      SUBROUTINE ZCOOSM(TRANST,M,N,UNITD,D,ALPHA,DESCRA,A,IA,JA,INFOA,
     *                  B,LDB,BETA,C,LDC,WORK,LWORK,IERROR)
      use psb_const_mod
      use psb_error_mod
      IMPLICIT NONE
-      COMPLEX*16  ALPHA, BETA
+      complex(psb_dpk_)  ALPHA, BETA
      INTEGER           LDB, LDC, LWORK, M, N, IERROR
      CHARACTER         UNITD, TRANST
-      COMPLEX*16        A(*), B(LDB,*), C(LDC,*), D(*), WORK(*)
+      complex(psb_dpk_)        A(*), B(LDB,*), C(LDC,*), D(*), WORK(*)
      INTEGER           IA(*), JA(*), INFOA(*), INT_VAL(5)
      CHARACTER         DESCRA*11
      INTEGER           I, K, ERR_ACT
--- a/base/serial/coo/zcoosv.f
+++ b/base/serial/coo/zcoosv.f
@ -34,14 +34,15 @@ C  "right" place, i.e. the last in its row for Lower and the first
 C  for Upper. 
 C
      SUBROUTINE ZCOOSV (UPLO,TRANS,DIAG,N,AS,IA,JA,INFOA,B,X,IERROR)
      use psb_const_mod
      use psb_error_mod
-      COMPLEX*16         ZERO
+      complex(psb_dpk_)         ZERO
      PARAMETER         (ZERO=(0.0D0,0.0D0))
      INTEGER           N,IERROR
      CHARACTER         DIAG, TRANS, UPLO
-      COMPLEX*16        AS(*), B(*), X(*)
+      complex(psb_dpk_)        AS(*), B(*), X(*)
      INTEGER           IA(*), JA(*),INFOA(*)
-      COMPLEX*16        ACC
+      complex(psb_dpk_)        ACC
      INTEGER           I, J, K, NNZ, II, JJ
      LOGICAL           LOW, TRA, UNI
      integer           :: debug_level, debug_unit
--- a/base/serial/csr/dcrnrmi.f
+++ b/base/serial/csr/dcrnrmi.f
@ -30,8 +30,9 @@ C POSSIBILITY OF SUCH DAMAGE.
 C
 C 
 C     ... Compute Infinity norm for sparse matrix in CSR Format ...
-      DOUBLE PRECISION FUNCTION DCRNRMI(TRANS,M,N,DESCRA,A,IA1,IA2,
+      real(psb_dpk_) FUNCTION DCRNRMI(TRANS,M,N,DESCRA,A,IA1,IA2,
     +  INFOA,IERROR)
      use psb_const_mod
      IMPLICIT NONE
 C     .. Scalar Arguments ..
      INTEGER           M,N, IERROR
@ -39,10 +40,10 @@ C     .. Scalar Arguments ..
 C     .. Array Arguments ..
      INTEGER           IA1(*),IA2(*),INFOA(*)
      CHARACTER         DESCRA*11
-      DOUBLE PRECISION  A(*)
+      real(psb_dpk_)  A(*)
 C     .. Local scalars ..
      INTEGER I, J
-      DOUBLE PRECISION NRMI, SUM
+      real(psb_dpk_) NRMI, SUM
      IERROR=0
      NRMI = 0.0
--- a/base/serial/csr/dcsrck.f
+++ b/base/serial/csr/dcsrck.f
@ -89,12 +89,13 @@ C             IERROR > 0 error in integrity check
      SUBROUTINE DCSRCK(TRANS,M,N,DESCRA,A,IA1,IA2,                       
     +   WORK,LWORK,IERROR)
      use psb_const_mod
      IMPLICIT NONE                                                     
 C     .. Scalar Arguments ..
      INTEGER          LWORK,M, N, IERROR
      CHARACTER        TRANS
 C     .. Array Arguments ..
-      DOUBLE PRECISION A(*), WORK(*)
+      real(psb_dpk_) A(*), WORK(*)
      INTEGER          IA1(*), IA2(*)
      CHARACTER        DESCRA*11
 C     .. Local Scalars ..
--- a/base/serial/csr/dcsrmm.f
+++ b/base/serial/csr/dcsrmm.f
@ -37,14 +37,13 @@ c
 *     
      SUBROUTINE DCSRMM(TRANSA,M,K,N,ALPHA,DESCRA,AR,
     *  JA,IA,B,LDB,BETA,C,LDC,WORK,LWORK,IERROR)
-C     
+      use psb_const_mod
 C     
 C     .. Scalar Arguments ..
-      DOUBLE PRECISION  ALPHA, BETA
+      real(psb_dpk_)  ALPHA, BETA
      INTEGER           K, LDB, LDC, M, N, LWORK,IERROR
      CHARACTER         TRANSA
 C     .. Array Arguments ..
-      DOUBLE PRECISION  AR(*), B(LDB,*), C(LDC,*),  WORK(*)
+      real(psb_dpk_)  AR(*), B(LDB,*), C(LDC,*),  WORK(*)
      INTEGER           IA(*), JA(*)
      CHARACTER         DESCRA*11
 C     .. Local Scalars ..
--- a/base/serial/csr/dcsrmv.f
+++ b/base/serial/csr/dcsrmv.f
@ -167,18 +167,19 @@ C
 ***********************************************************************
      SUBROUTINE DCSRMV(TRANS,DIAG,M,N,ALPHA,AS,JA,IA,X,BETA,Y,
     +  WORK,LWORK,IERROR)
      use psb_const_mod
 C     .. Parameters ..
-      DOUBLE PRECISION  ONE, ZERO
+      real(psb_dpk_)  ONE, ZERO
      PARAMETER         (ONE=1.0D0,ZERO=0.0D0)
 C     .. Scalar Arguments ..
-      DOUBLE PRECISION  ALPHA, BETA
+      real(psb_dpk_)  ALPHA, BETA
      INTEGER           M, N,LWORK,IERROR
      CHARACTER         DIAG, TRANS
 C     .. Array Arguments ..
-      DOUBLE PRECISION  AS(*), WORK(*), X(*), Y(*)
+      real(psb_dpk_)  AS(*), WORK(*), X(*), Y(*)
      INTEGER           IA(*), JA(*)
 C     .. Local Scalars ..
-      DOUBLE PRECISION  ACC
+      real(psb_dpk_)  ACC
      INTEGER           I, J, K, NCOLA, NROWA
      LOGICAL           SYM, TRA, UNI
 C     .. Executable Statements ..
--- a/base/serial/csr/dcsrmv2.f
+++ b/base/serial/csr/dcsrmv2.f
@ -167,20 +167,21 @@ C
 ***********************************************************************
      SUBROUTINE DCSRMV2(TRANS,DIAG,M,N,ALPHA,AS,JA,IA,X,LDX,
     +  BETA,Y,LDY, WORK,LWORK,IERROR)
      use psb_const_mod
      integer  nb
      parameter (nb=2)
 C     .. Parameters ..
-      DOUBLE PRECISION  ONE, ZERO
+      real(psb_dpk_)  ONE, ZERO
      PARAMETER         (ONE=1.0D0,ZERO=0.0D0)
 C     .. Scalar Arguments ..
-      DOUBLE PRECISION  ALPHA, BETA
+      real(psb_dpk_)  ALPHA, BETA
      INTEGER           M, N,LWORK,IERROR,ldx,ldy
      CHARACTER         DIAG, TRANS
 C     .. Array Arguments ..
-      DOUBLE PRECISION  AS(*), WORK(*), X(LDX,NB), Y(LDY,NB)
+      real(psb_dpk_)  AS(*), WORK(*), X(LDX,NB), Y(LDY,NB)
      INTEGER           IA(*), JA(*)
 C     .. Local Scalars ..
-      DOUBLE PRECISION  ACC(nb)
+      real(psb_dpk_)  ACC(nb)
      INTEGER           I, J, K, NCOLA, NROWA
      LOGICAL           SYM, TRA, UNI
 C     .. Executable Statements ..
--- a/base/serial/csr/dcsrmv3.f
+++ b/base/serial/csr/dcsrmv3.f
@ -167,20 +167,21 @@ C
 ***********************************************************************
      SUBROUTINE DCSRMV3(TRANS,DIAG,M,N,ALPHA,AS,JA,IA,X,LDX,
     +  BETA,Y,LDY, WORK,LWORK,IERROR)
      use psb_const_mod
      integer  nb
      parameter (nb=3)
 C     .. Parameters ..
-      DOUBLE PRECISION  ONE, ZERO
+      real(psb_dpk_)  ONE, ZERO
      PARAMETER         (ONE=1.0D0,ZERO=0.0D0)
 C     .. Scalar Arguments ..
-      DOUBLE PRECISION  ALPHA, BETA
+      real(psb_dpk_)  ALPHA, BETA
      INTEGER           M, N,LWORK,IERROR,ldx,ldy
      CHARACTER         DIAG, TRANS
 C     .. Array Arguments ..
-      DOUBLE PRECISION  AS(*), WORK(*), X(LDX,NB), Y(LDY,NB)
+      real(psb_dpk_)  AS(*), WORK(*), X(LDX,NB), Y(LDY,NB)
      INTEGER           IA(*), JA(*)
 C     .. Local Scalars ..
-      DOUBLE PRECISION  ACC(nb)
+      real(psb_dpk_)  ACC(nb)
      INTEGER           I, J, K, NCOLA, NROWA
      LOGICAL           SYM, TRA, UNI
 C     .. Executable Statements ..
--- a/base/serial/csr/dcsrmv4.f
+++ b/base/serial/csr/dcsrmv4.f
@ -167,20 +167,21 @@ C
 ***********************************************************************
      SUBROUTINE DCSRMV4(TRANS,DIAG,M,N,ALPHA,AS,JA,IA,X,LDX,
     +  BETA,Y,LDY, WORK,LWORK,IERROR)
      use psb_const_mod
      integer  nb
      parameter (nb=4)
 C     .. Parameters ..
-      DOUBLE PRECISION  ONE, ZERO
+      real(psb_dpk_)  ONE, ZERO
      PARAMETER         (ONE=1.0D0,ZERO=0.0D0)
 C     .. Scalar Arguments ..
-      DOUBLE PRECISION  ALPHA, BETA
+      real(psb_dpk_)  ALPHA, BETA
      INTEGER           M, N,LWORK,IERROR,ldx,ldy
      CHARACTER         DIAG, TRANS
 C     .. Array Arguments ..
-      DOUBLE PRECISION  AS(*), WORK(*), X(LDX,NB), Y(LDY,NB)
+      real(psb_dpk_)  AS(*), WORK(*), X(LDX,NB), Y(LDY,NB)
      INTEGER           IA(*), JA(*)
 C     .. Local Scalars ..
-      DOUBLE PRECISION  ACC(nb)
+      real(psb_dpk_)  ACC(nb)
      INTEGER           I, J, K, NCOLA, NROWA
      LOGICAL           SYM, TRA, UNI
 C     .. Executable Statements ..
--- a/base/serial/csr/dcsrrws.f
+++ b/base/serial/csr/dcsrrws.f
@ -31,6 +31,7 @@ C
 C 
      SUBROUTINE DCSRRWS(TRANS,M,N,DESCRA,A,IA1,IA2,
     &  INFOA,ROWSUM,IERROR)
      use psb_const_mod
      IMPLICIT NONE
 C     .. Scalar Arguments ..
      INTEGER           M,N, IERROR
@ -38,7 +39,7 @@ C     .. Scalar Arguments ..
 C     .. Array Arguments ..
      INTEGER           IA1(*),IA2(*),INFOA(*)
      CHARACTER         DESCRA*11
-      DOUBLE PRECISION  A(*), ROWSUM(*)
+      real(psb_dpk_)  A(*), ROWSUM(*)
 C     .. Local scalars ..
      INTEGER I, J
--- a/base/serial/csr/dcsrsm.f
+++ b/base/serial/csr/dcsrsm.f
@ -31,10 +31,11 @@ C
 C 
      SUBROUTINE DCSRSM(TRANST,M,N,UNITD,D,ALPHA,DESCRA,A,JA,IA,
     *                  B,LDB,BETA,C,LDC,WORK,LWORK,IERROR)
-      DOUBLE PRECISION  ALPHA, BETA
+      use psb_const_mod
      real(psb_dpk_)  ALPHA, BETA
      INTEGER           LDB, LDC, LWORK, M, N, IERROR
      CHARACTER         UNITD, TRANST
-      DOUBLE PRECISION  A(*), B(LDB,*), C(LDC,*), D(*), WORK(*)
+      real(psb_dpk_)  A(*), B(LDB,*), C(LDC,*), D(*), WORK(*)
      INTEGER           IA(*), JA(*)
      CHARACTER         DESCRA*11
      INTEGER           I, K
--- a/base/serial/csr/dcsrsv.f
+++ b/base/serial/csr/dcsrsv.f
@ -30,13 +30,14 @@ C POSSIBILITY OF SUCH DAMAGE.
 C
 C 
      SUBROUTINE DCSRSV(UPLO,TRANS,DIAG,N,AS,JA,IA,B,X,IERROR)
-      DOUBLE PRECISION   ZERO
+      use psb_const_mod
      real(psb_dpk_)   ZERO
      PARAMETER         (ZERO=0.0D0)
      INTEGER           N, IERROR
      CHARACTER         DIAG, TRANS, UPLO
-      DOUBLE PRECISION  AS(*), B(*), X(*)
+      real(psb_dpk_)  AS(*), B(*), X(*)
      INTEGER           IA(*), JA(*)
-      DOUBLE PRECISION  ACC
+      real(psb_dpk_)  ACC
      INTEGER           I, J, K
      LOGICAL           LOW, TRA, UNI
      UNI = (DIAG.EQ.'U')
--- a/base/serial/csr/zcrnrmi.f
+++ b/base/serial/csr/zcrnrmi.f
@ -28,8 +28,9 @@ C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 C POSSIBILITY OF SUCH DAMAGE.
 C
-      DOUBLE PRECISION FUNCTION ZCRNRMI(TRANS,M,N,DESCRA,A,IA1,IA2,      
+      real(psb_dpk_) FUNCTION ZCRNRMI(TRANS,M,N,DESCRA,A,IA1,IA2,      
     &   INFOA,IERROR)
      use psb_const_mod
      IMPLICIT NONE
 C     .. Scalar Arguments ..
      INTEGER           M,N, IERROR
@ -37,10 +38,10 @@ C     .. Scalar Arguments ..
 C     .. Array Arguments ..
      INTEGER    IA1(*),IA2(*),INFOA(*)
      CHARACTER  DESCRA*11
-      COMPLEX*16 A(*)
+      complex(psb_dpk_) A(*)
 C     .. Local scalars ..
      INTEGER I, J
-      DOUBLE PRECISION NRMI, SUM
+      real(psb_dpk_) NRMI, SUM
      NRMI = 0.D0
      DO I = 1, M
--- a/base/serial/csr/zcsrck.f
+++ b/base/serial/csr/zcsrck.f
@ -89,12 +89,13 @@ C             IERROR > 0 error in integrity check
      SUBROUTINE ZCSRCK(TRANS,M,N,DESCRA,A,IA1,IA2,                       
     +   WORK,LWORK,IERROR)
      use psb_const_mod
      IMPLICIT NONE                                                     
 C     .. Scalar Arguments ..
      INTEGER          LWORK,M, N, IERROR
      CHARACTER        TRANS
 C     .. Array Arguments ..
-      complex(kind(1.d0)) A(*), WORK(*)
+      complex(psb_dpk_) A(*), WORK(*)
      INTEGER          IA1(*), IA2(*)
      CHARACTER        DESCRA*11
 C     .. Local Scalars ..
--- a/base/serial/csr/zcsrmm.f
+++ b/base/serial/csr/zcsrmm.f
@ -136,14 +136,13 @@ C             Unchanged on exit.
 C
      SUBROUTINE ZCSRMM(TRANSA,M,K,N,ALPHA,DESCRA,AR,
     *                  JA,IA,B,LDB,BETA,C,LDC,WORK,LWORK)
-C
+      use psb_const_mod
 C
 C     .. Scalar Arguments ..
-      COMPLEX*16 ALPHA, BETA
+      complex(psb_dpk_) ALPHA, BETA
      INTEGER    K, LDB, LDC, M, N, LWORK
      CHARACTER  TRANSA
 C     .. Array Arguments ..
-      COMPLEX*16 AR(*), B(LDB,*), C(LDC,*),  WORK(*)
+      complex(psb_dpk_) AR(*), B(LDB,*), C(LDC,*),  WORK(*)
      INTEGER    IA(*), JA(*)
      CHARACTER  DESCRA*11
 C     .. Local Scalars ..
--- a/base/serial/csr/zcsrrws.f
+++ b/base/serial/csr/zcsrrws.f
@ -30,6 +30,7 @@ C POSSIBILITY OF SUCH DAMAGE.
 C
      SUBROUTINE ZCSRRWS(TRANS,M,N,DESCRA,A,IA1,IA2,
     &  INFOA,ROWSUM,IERROR)
      use psb_const_mod
      IMPLICIT NONE
 C     .. Scalar Arguments ..
      INTEGER           M,N, IERROR
@ -37,10 +38,10 @@ C     .. Scalar Arguments ..
 C     .. Array Arguments ..
      INTEGER           IA1(*),IA2(*),INFOA(*)
      CHARACTER         DESCRA*11
-      COMPLEX*16        A(*), ROWSUM(*)
+      complex(psb_dpk_)        A(*), ROWSUM(*)
 C     .. Local scalars ..
      INTEGER I, J
-      DOUBLE PRECISION  SUM
+      real(psb_dpk_)  SUM
      IF (TRANS.EQ.'N') THEN
        DO I = 1, M
--- a/base/serial/csr/zcsrsm.f
+++ b/base/serial/csr/zcsrsm.f
@ -31,10 +31,11 @@ C
 C 
      SUBROUTINE ZCSRSM(TRANST,M,N,UNITD,D,ALPHA,DESCRA,A,JA,IA,
     *                  B,LDB,BETA,C,LDC,WORK,LWORK)
-      COMPLEX*16 ALPHA, BETA
+      use psb_const_mod
      complex(psb_dpk_) ALPHA, BETA
      INTEGER    LDB, LDC, LWORK, M, N
      CHARACTER  UNITD, TRANST
-      COMPLEX*16 A(*), B(LDB,*), C(LDC,*), D(*), WORK(*)
+      complex(psb_dpk_) A(*), B(LDB,*), C(LDC,*), D(*), WORK(*)
      INTEGER    IA(*), JA(*)
      CHARACTER  DESCRA*11
      INTEGER    I, K
--- a/base/serial/csr/zsrmv.f
+++ b/base/serial/csr/zsrmv.f
@ -178,18 +178,19 @@ C
 ***********************************************************************
 C
      SUBROUTINE ZSRMV (TRANS,DIAG,M,N,ALPHA,AS,JA,IA,X,BETA,Y,WORK)
      use psb_const_mod
 C     .. Parameters ..
-      COMPLEX*16 ONE, ZERO
+      complex(psb_dpk_) ONE, ZERO
      PARAMETER (ONE=(1.0D0, 0.0D0), ZERO=(0.0D0, 0.0D0))
 C     .. Scalar Arguments ..
-      COMPLEX*16 ALPHA, BETA
+      complex(psb_dpk_) ALPHA, BETA
      INTEGER    M, N
      CHARACTER  DIAG, TRANS
 C     .. Array Arguments ..
-      COMPLEX*16 AS(*), WORK(*), X(*), Y(*)
+      complex(psb_dpk_) AS(*), WORK(*), X(*), Y(*)
      INTEGER    IA(*), JA(*)
 C     .. Local Scalars ..
-      COMPLEX*16 ACC
+      complex(psb_dpk_) ACC
      INTEGER    I, J, K, NCOLA, NROWA,DUM
      LOGICAL    SYM, TRA, COTRA, UNI
 C     .. Executable Statements ..
--- a/base/serial/csr/zsrsv.f
+++ b/base/serial/csr/zsrsv.f
@ -29,13 +29,14 @@ C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 C POSSIBILITY OF SUCH DAMAGE.
 C
      SUBROUTINE ZSRSV (UPLO,TRANS,DIAG,N,AS,JA,IA,B,X)
-      COMPLEX*16 ZERO
+      use psb_const_mod
      complex(psb_dpk_) ZERO
      PARAMETER  (ZERO = (0.0D0, 0.0D0))
      INTEGER    N
      CHARACTER  DIAG, TRANS, UPLO
-      COMPLEX*16 AS(*), B(*), X(*)
+      complex(psb_dpk_) AS(*), B(*), X(*)
      INTEGER    IA(*), JA(*)
-      COMPLEX*16 ACC
+      complex(psb_dpk_) ACC
      INTEGER    I, J, K
      LOGICAL    LOW, TRA, COTRA, UNI
      UNI = (DIAG.EQ.'U')
--- a/Show More
+++ b/Show More