diff --git a/base/modules/psb_hash_mod.f90 b/base/modules/psb_hash_mod.f90
index 649537e3..603aefa5 100644
--- a/base/modules/psb_hash_mod.f90
+++ b/base/modules/psb_hash_mod.f90
@@ -30,26 +30,28 @@
 !!$ 
 !!$  
 !
-! package: psb_hash_mod
-! 
-!   This module implements a very simple minded hash table.
-!   The hash is based on the idea of open addressing with double hashing;
-!   the primary hash function h1(K) is simply the remainder modulo 2^N, while
-!   the secondary hash function is 1 if H1(k) == 0, otherwise IOR((2^N-H1(k)),1)
-!   (See Knuth: TAOCP, Vol. 3, sec. 6.4)
-!   These hash functions are not very smart; however they are very simple and fast.
-!   The intended usage of this hash table is to store indices of halo points, which
-!   are supposed to be few compared to the internal indices (which are stored elsewhere).
-!   Therefore, either the table has a very low occupancy, and this scheme will work,
-!   or we have lots more to worry about in parallel performance than the efficiency
-!   of this hashing scheme.
-!
-!
 !
 module psb_hash_mod
   use psb_const_mod
   
-  !
+  !> \class psb_hash_mod
+  !! \brief Simple hash module for storing integer keys. 
+  !! 
+  !!   This module implements a very simple minded hash table.
+  !!   The hash is based on the idea of open addressing with double hashing;
+  !!   the primary hash function h1(K) is simply the remainder modulo 2^N, while
+  !!   the secondary hash function is 1 if H1(k) == 0, otherwise IOR((2^N-H1(k)),1)
+  !!   (See Knuth: TAOCP, Vol. 3, sec. 6.4)
+  !!   
+  !!   These hash functions are not very smart; however they are very simple and fast.
+  !!   The intended usage of this hash table is to store indices of halo points, which
+  !!   are supposed to be few compared to the internal indices
+  !!   (which are stored elsewhere).
+  !!   Therefore, either the table has a very low occupancy, and this scheme will work,
+  !!   or we have lots more to worry about in parallel performance than the efficiency
+  !!   of this hashing scheme.
+  !!
+  !!
   ! For us a hash is a Nx2 table.
   ! Note: we are assuming that the keys are positive numbers.
   ! Allocatable scalars would be a nice solution...
diff --git a/test/pargen/ppde2d.f90 b/test/pargen/ppde2d.f90
index f9a845d7..922c4925 100644
--- a/test/pargen/ppde2d.f90
+++ b/test/pargen/ppde2d.f90
@@ -55,11 +55,62 @@
 ! then the corresponding vector is distributed according to a BLOCK
 ! data distribution.
 !
+module ppde2d_mod
+contains
+
+  !
+  ! functions parametrizing the differential equation 
+  !  
+  function b1(x,y)
+    use psb_base_mod, only : psb_dpk_
+    real(psb_dpk_) :: b1
+    real(psb_dpk_), intent(in) :: x,y
+    b1=1.d0/sqrt(2.d0)
+  end function b1
+  function b2(x,y)
+    use psb_base_mod, only : psb_dpk_
+    real(psb_dpk_) ::  b2
+    real(psb_dpk_), intent(in) :: x,y
+    b2=1.d0/sqrt(2.d0)
+  end function b2
+  function c(x,y)
+    use psb_base_mod, only : psb_dpk_
+    real(psb_dpk_) ::  c
+    real(psb_dpk_), intent(in) :: x,y
+    c=0.d0
+  end function c
+  function a1(x,y)
+    use psb_base_mod, only : psb_dpk_
+    real(psb_dpk_) ::  a1   
+    real(psb_dpk_), intent(in) :: x,y
+    a1=1.d0/80
+  end function a1
+  function a2(x,y)
+    use psb_base_mod, only : psb_dpk_
+    real(psb_dpk_) ::  a2
+    real(psb_dpk_), intent(in) :: x,y
+    a2=1.d0/80
+  end function a2
+  function g(x,y)
+    use psb_base_mod, only : psb_dpk_, done, dzero
+    real(psb_dpk_) ::  g
+    real(psb_dpk_), intent(in) :: x,y
+    g = dzero
+    if (x == done) then
+      g = done
+    else if (x == dzero) then 
+      g = exp(-y**2)
+    end if
+  end function g
+
+end module ppde2d_mod
+
 program ppde2d
   use psb_base_mod
   use psb_prec_mod
   use psb_krylov_mod
   use psb_util_mod
+  use ppde2d_mod
   implicit none
 
   ! input parameters
@@ -321,52 +372,6 @@ contains
     write(iout,*)'               >= 1 do tracing every itrace'
     write(iout,*)'               iterations ' 
   end subroutine pr_usage
-
-  !
-  ! functions parametrizing the differential equation 
-  !  
-  function b1(x,y)
-    use psb_base_mod, only : psb_dpk_
-    real(psb_dpk_) :: b1
-    real(psb_dpk_), intent(in) :: x,y
-    b1=1.d0/sqrt(2.d0)
-  end function b1
-  function b2(x,y)
-    use psb_base_mod, only : psb_dpk_
-    real(psb_dpk_) ::  b2
-    real(psb_dpk_), intent(in) :: x,y
-    b2=1.d0/sqrt(2.d0)
-  end function b2
-  function c(x,y)
-    use psb_base_mod, only : psb_dpk_
-    real(psb_dpk_) ::  c
-    real(psb_dpk_), intent(in) :: x,y
-    c=0.d0
-  end function c
-  function a1(x,y)
-    use psb_base_mod, only : psb_dpk_
-    real(psb_dpk_) ::  a1   
-    real(psb_dpk_), intent(in) :: x,y
-    a1=1.d0/80
-  end function a1
-  function a2(x,y)
-    use psb_base_mod, only : psb_dpk_
-    real(psb_dpk_) ::  a2
-    real(psb_dpk_), intent(in) :: x,y
-    a2=1.d0/80
-  end function a2
-  function g(x,y)
-    use psb_base_mod, only : psb_dpk_, done
-    real(psb_dpk_) ::  g
-    real(psb_dpk_), intent(in) :: x,y
-    g = dzero
-    if (x == done) then
-      g = done
-    else if (x == dzero) then 
-      g = exp(-y**2)
-    end if
-  end function g
-
 end program ppde2d
 
 
diff --git a/test/pargen/ppde3d.f90 b/test/pargen/ppde3d.f90
index 2b300db8..a8ad8f1d 100644
--- a/test/pargen/ppde3d.f90
+++ b/test/pargen/ppde3d.f90
@@ -55,11 +55,73 @@
 ! then the corresponding vector is distributed according to a BLOCK
 ! data distribution.
 !
+module ppde3d_mod
+contains
+  !
+  ! functions parametrizing the differential equation 
+  !  
+  function b1(x,y,z)
+    use psb_base_mod, only : psb_dpk_
+    real(psb_dpk_) :: b1
+    real(psb_dpk_), intent(in) :: x,y,z
+    b1=1.d0/sqrt(3.d0)
+  end function b1
+  function b2(x,y,z)
+    use psb_base_mod, only : psb_dpk_
+    real(psb_dpk_) ::  b2
+    real(psb_dpk_), intent(in) :: x,y,z
+    b2=1.d0/sqrt(3.d0)
+  end function b2
+  function b3(x,y,z)
+    use psb_base_mod, only : psb_dpk_
+    real(psb_dpk_) ::  b3
+    real(psb_dpk_), intent(in) :: x,y,z      
+    b3=1.d0/sqrt(3.d0)
+  end function b3
+  function c(x,y,z)
+    use psb_base_mod, only : psb_dpk_
+    real(psb_dpk_) ::  c
+    real(psb_dpk_), intent(in) :: x,y,z      
+    c=0.d0
+  end function c
+  function a1(x,y,z)
+    use psb_base_mod, only : psb_dpk_
+    real(psb_dpk_) ::  a1   
+    real(psb_dpk_), intent(in) :: x,y,z
+    a1=1.d0/80
+  end function a1
+  function a2(x,y,z)
+    use psb_base_mod, only : psb_dpk_
+    real(psb_dpk_) ::  a2
+    real(psb_dpk_), intent(in) :: x,y,z
+    a2=1.d0/80
+  end function a2
+  function a3(x,y,z)
+    use psb_base_mod, only : psb_dpk_
+    real(psb_dpk_) ::  a3
+    real(psb_dpk_), intent(in) :: x,y,z
+    a3=1.d0/80
+  end function a3
+  function g(x,y,z)
+    use psb_base_mod, only : psb_dpk_, done, dzero
+    real(psb_dpk_) ::  g
+    real(psb_dpk_), intent(in) :: x,y,z
+    g = dzero
+    if (x == done) then
+      g = done
+    else if (x == dzero) then 
+      g = exp(y**2-z**2)
+    end if
+  end function g
+
+end module ppde3d_mod
+
 program ppde3d
   use psb_base_mod
   use psb_prec_mod
   use psb_krylov_mod
   use psb_util_mod
+  use ppde3d_mod
   implicit none
 
   ! input parameters
@@ -325,63 +387,6 @@ contains
     write(iout,*)'               iterations ' 
   end subroutine pr_usage
 
-  !
-  ! functions parametrizing the differential equation 
-  !  
-  function b1(x,y,z)
-    use psb_base_mod, only : psb_dpk_
-    real(psb_dpk_) :: b1
-    real(psb_dpk_), intent(in) :: x,y,z
-    b1=1.d0/sqrt(3.d0)
-  end function b1
-  function b2(x,y,z)
-    use psb_base_mod, only : psb_dpk_
-    real(psb_dpk_) ::  b2
-    real(psb_dpk_), intent(in) :: x,y,z
-    b2=1.d0/sqrt(3.d0)
-  end function b2
-  function b3(x,y,z)
-    use psb_base_mod, only : psb_dpk_
-    real(psb_dpk_) ::  b3
-    real(psb_dpk_), intent(in) :: x,y,z      
-    b3=1.d0/sqrt(3.d0)
-  end function b3
-  function c(x,y,z)
-    use psb_base_mod, only : psb_dpk_
-    real(psb_dpk_) ::  c
-    real(psb_dpk_), intent(in) :: x,y,z      
-    c=0.d0
-  end function c
-  function a1(x,y,z)
-    use psb_base_mod, only : psb_dpk_
-    real(psb_dpk_) ::  a1   
-    real(psb_dpk_), intent(in) :: x,y,z
-    a1=1.d0/80
-  end function a1
-  function a2(x,y,z)
-    use psb_base_mod, only : psb_dpk_
-    real(psb_dpk_) ::  a2
-    real(psb_dpk_), intent(in) :: x,y,z
-    a2=1.d0/80
-  end function a2
-  function a3(x,y,z)
-    use psb_base_mod, only : psb_dpk_
-    real(psb_dpk_) ::  a3
-    real(psb_dpk_), intent(in) :: x,y,z
-    a3=1.d0/80
-  end function a3
-  function g(x,y,z)
-    use psb_base_mod, only : psb_dpk_, done
-    real(psb_dpk_) ::  g
-    real(psb_dpk_), intent(in) :: x,y,z
-    g = dzero
-    if (x == done) then
-      g = done
-    else if (x == dzero) then 
-      g = exp(y**2-z**2)
-    end if
-  end function g
-
 end program ppde3d
 
 
diff --git a/test/pargen/spde2d.f90 b/test/pargen/spde2d.f90
index 4d7ce462..ffa5d60d 100644
--- a/test/pargen/spde2d.f90
+++ b/test/pargen/spde2d.f90
@@ -55,11 +55,61 @@
 ! then the corresponding vector is distributed according to a BLOCK
 ! data distribution.
 !
+module spde2d_mod
+contains
+  !
+  ! functions parametrizing the differential equation 
+  !  
+  function b1(x,y)
+    use psb_base_mod, only : psb_spk_
+    real(psb_spk_) :: b1
+    real(psb_spk_), intent(in) :: x,y
+    b1=1.e0/sqrt(2.e0)
+  end function b1
+  function b2(x,y)
+    use psb_base_mod, only : psb_spk_
+    real(psb_spk_) ::  b2
+    real(psb_spk_), intent(in) :: x,y
+    b2=1.e0/sqrt(2.e0)
+  end function b2
+  function c(x,y)
+    use psb_base_mod, only : psb_spk_
+    real(psb_spk_) ::  c
+    real(psb_spk_), intent(in) :: x,y
+    c=0.e0
+  end function c
+  function a1(x,y)
+    use psb_base_mod, only : psb_spk_
+    real(psb_spk_) ::  a1   
+    real(psb_spk_), intent(in) :: x,y
+    a1=1.e0/80
+  end function a1
+  function a2(x,y)
+    use psb_base_mod, only : psb_spk_
+    real(psb_spk_) ::  a2
+    real(psb_spk_), intent(in) :: x,y
+    a2=1.e0/80
+  end function a2
+  function g(x,y)
+    use psb_base_mod, only : psb_spk_, sone, szero
+    real(psb_spk_) ::  g
+    real(psb_spk_), intent(in) :: x,y
+    g = szero
+    if (x == sone) then
+      g = sone
+    else if (x == szero) then 
+      g = exp(-y**2)
+    end if
+  end function g
+
+end module spde2d_mod
+
 program spde2d
   use psb_base_mod
   use psb_prec_mod
   use psb_krylov_mod
   use psb_util_mod
+  use spde2d_mod
   implicit none
 
   ! input parameters
@@ -322,51 +372,6 @@ contains
     write(iout,*)'               iterations ' 
   end subroutine pr_usage
 
-  !
-  ! functions parametrizing the differential equation 
-  !  
-  function b1(x,y)
-    use psb_base_mod, only : psb_spk_
-    real(psb_spk_) :: b1
-    real(psb_spk_), intent(in) :: x,y
-    b1=1.e0/sqrt(2.e0)
-  end function b1
-  function b2(x,y)
-    use psb_base_mod, only : psb_spk_
-    real(psb_spk_) ::  b2
-    real(psb_spk_), intent(in) :: x,y
-    b2=1.e0/sqrt(2.e0)
-  end function b2
-  function c(x,y)
-    use psb_base_mod, only : psb_spk_
-    real(psb_spk_) ::  c
-    real(psb_spk_), intent(in) :: x,y
-    c=0.e0
-  end function c
-  function a1(x,y)
-    use psb_base_mod, only : psb_spk_
-    real(psb_spk_) ::  a1   
-    real(psb_spk_), intent(in) :: x,y
-    a1=1.e0/80
-  end function a1
-  function a2(x,y)
-    use psb_base_mod, only : psb_spk_
-    real(psb_spk_) ::  a2
-    real(psb_spk_), intent(in) :: x,y
-    a2=1.e0/80
-  end function a2
-  function g(x,y)
-    use psb_base_mod, only : psb_spk_, sone
-    real(psb_spk_) ::  g
-    real(psb_spk_), intent(in) :: x,y
-    g = szero
-    if (x == sone) then
-      g = sone
-    else if (x == szero) then 
-      g = exp(-y**2)
-    end if
-  end function g
-
 end program spde2d
 
 
diff --git a/test/pargen/spde3d.f90 b/test/pargen/spde3d.f90
index 9eadb27f..554ca997 100644
--- a/test/pargen/spde3d.f90
+++ b/test/pargen/spde3d.f90
@@ -56,11 +56,73 @@
 ! data distribution.
 !
 !
+module spde3d_mod
+contains
+  
+  !
+  ! functions parametrizing the differential equation 
+  !  
+  function b1(x,y,z)
+    use psb_base_mod, only : psb_spk_
+    real(psb_spk_) :: b1
+    real(psb_spk_), intent(in) :: x,y,z
+    b1=1.e0/sqrt(3.e0)
+  end function b1
+  function b2(x,y,z)
+    use psb_base_mod, only : psb_spk_
+    real(psb_spk_) ::  b2
+    real(psb_spk_), intent(in) :: x,y,z
+    b2=1.e0/sqrt(3.e0)
+  end function b2
+  function b3(x,y,z)
+    use psb_base_mod, only : psb_spk_
+    real(psb_spk_) ::  b3
+    real(psb_spk_), intent(in) :: x,y,z      
+    b3=1.e0/sqrt(3.e0)
+  end function b3
+  function c(x,y,z)
+    use psb_base_mod, only : psb_spk_
+    real(psb_spk_) ::  c
+    real(psb_spk_), intent(in) :: x,y,z      
+    c=0.e0
+  end function c
+  function a1(x,y,z)
+    use psb_base_mod, only : psb_spk_
+    real(psb_spk_) ::  a1   
+    real(psb_spk_), intent(in) :: x,y,z
+    a1=1.e0/80
+  end function a1
+  function a2(x,y,z)
+    use psb_base_mod, only : psb_spk_
+    real(psb_spk_) ::  a2
+    real(psb_spk_), intent(in) :: x,y,z
+    a2=1.e0/80
+  end function a2
+  function a3(x,y,z)
+    use psb_base_mod, only : psb_spk_
+    real(psb_spk_) ::  a3
+    real(psb_spk_), intent(in) :: x,y,z
+    a3=1.e0/80
+  end function a3
+  function g(x,y,z)
+    use psb_base_mod, only : psb_spk_, sone, szero
+    real(psb_spk_) ::  g
+    real(psb_spk_), intent(in) :: x,y,z
+    g = szero
+    if (x == sone) then
+      g = sone
+    else if (x == szero) then 
+      g = exp(y**2-z**2)
+    end if
+  end function g
+end module spde3d_mod
+
 program spde3d
   use psb_base_mod
   use psb_prec_mod
   use psb_krylov_mod
   use psb_util_mod
+  use spde3d_mod
   implicit none
 
   ! input parameters
@@ -326,63 +388,6 @@ contains
     write(iout,*)'               iterations ' 
   end subroutine pr_usage
 
-  !
-  ! functions parametrizing the differential equation 
-  !  
-  function b1(x,y,z)
-    use psb_base_mod, only : psb_spk_
-    real(psb_spk_) :: b1
-    real(psb_spk_), intent(in) :: x,y,z
-    b1=1.e0/sqrt(3.e0)
-  end function b1
-  function b2(x,y,z)
-    use psb_base_mod, only : psb_spk_
-    real(psb_spk_) ::  b2
-    real(psb_spk_), intent(in) :: x,y,z
-    b2=1.e0/sqrt(3.e0)
-  end function b2
-  function b3(x,y,z)
-    use psb_base_mod, only : psb_spk_
-    real(psb_spk_) ::  b3
-    real(psb_spk_), intent(in) :: x,y,z      
-    b3=1.e0/sqrt(3.e0)
-  end function b3
-  function c(x,y,z)
-    use psb_base_mod, only : psb_spk_
-    real(psb_spk_) ::  c
-    real(psb_spk_), intent(in) :: x,y,z      
-    c=0.e0
-  end function c
-  function a1(x,y,z)
-    use psb_base_mod, only : psb_spk_
-    real(psb_spk_) ::  a1   
-    real(psb_spk_), intent(in) :: x,y,z
-    a1=1.e0/80
-  end function a1
-  function a2(x,y,z)
-    use psb_base_mod, only : psb_spk_
-    real(psb_spk_) ::  a2
-    real(psb_spk_), intent(in) :: x,y,z
-    a2=1.e0/80
-  end function a2
-  function a3(x,y,z)
-    use psb_base_mod, only : psb_spk_
-    real(psb_spk_) ::  a3
-    real(psb_spk_), intent(in) :: x,y,z
-    a3=1.e0/80
-  end function a3
-  function g(x,y,z)
-    use psb_base_mod, only : psb_spk_, sone
-    real(psb_spk_) ::  g
-    real(psb_spk_), intent(in) :: x,y,z
-    g = szero
-    if (x == sone) then
-      g = sone
-    else if (x == szero) then 
-      g = exp(y**2-z**2)
-    end if
-  end function g
-
 end program spde3d