diff --git a/tests/fileread/cf_sample.f90 b/tests/fileread/cf_sample.f90
index 920a190c..06dde492 100644
--- a/tests/fileread/cf_sample.f90
+++ b/tests/fileread/cf_sample.f90
@@ -264,7 +264,7 @@ program cf_sample
     call mld_precset(prec,mld_coarse_subsolve_, prec_choice%csbsolve,info)
     call mld_precset(prec,mld_coarse_mat_,      prec_choice%cmat,    info)
     call mld_precset(prec,mld_coarse_fillin_,   prec_choice%cfill,   info)
-    call mld_precset(prec,mld_coarse_iluthrs_,    prec_choice%cthres,  info)
+    call mld_precset(prec,mld_coarse_iluthrs_,  prec_choice%cthres,  info)
     call mld_precset(prec,mld_coarse_sweeps_,   prec_choice%cjswp,   info)
   end if
 
diff --git a/tests/fileread/df_sample.f90 b/tests/fileread/df_sample.f90
index acd4618e..c4aa58ed 100644
--- a/tests/fileread/df_sample.f90
+++ b/tests/fileread/df_sample.f90
@@ -264,7 +264,7 @@ program df_sample
     call mld_precset(prec,mld_coarse_subsolve_, prec_choice%csbsolve,info)
     call mld_precset(prec,mld_coarse_mat_,      prec_choice%cmat,    info)
     call mld_precset(prec,mld_coarse_fillin_,   prec_choice%cfill,   info)
-    call mld_precset(prec,mld_coarse_iluthrs_,    prec_choice%cthres,  info)
+    call mld_precset(prec,mld_coarse_iluthrs_,  prec_choice%cthres,  info)
     call mld_precset(prec,mld_coarse_sweeps_,   prec_choice%cjswp,   info)
   end if
 
diff --git a/tests/fileread/sf_sample.f90 b/tests/fileread/sf_sample.f90
index 4e09395e..a6bf71a6 100644
--- a/tests/fileread/sf_sample.f90
+++ b/tests/fileread/sf_sample.f90
@@ -264,7 +264,7 @@ program sf_sample
     call mld_precset(prec,mld_coarse_subsolve_, prec_choice%csbsolve,info)
     call mld_precset(prec,mld_coarse_mat_,      prec_choice%cmat,    info)
     call mld_precset(prec,mld_coarse_fillin_,   prec_choice%cfill,   info)
-    call mld_precset(prec,mld_coarse_iluthrs_,    prec_choice%cthres,  info)
+    call mld_precset(prec,mld_coarse_iluthrs_,  prec_choice%cthres,  info)
     call mld_precset(prec,mld_coarse_sweeps_,   prec_choice%cjswp,   info)
   end if
 
diff --git a/tests/fileread/zf_sample.f90 b/tests/fileread/zf_sample.f90
index 400e46ea..2bde9cc6 100644
--- a/tests/fileread/zf_sample.f90
+++ b/tests/fileread/zf_sample.f90
@@ -264,7 +264,7 @@ program zf_sample
     call mld_precset(prec,mld_coarse_subsolve_, prec_choice%csbsolve,info)
     call mld_precset(prec,mld_coarse_mat_,      prec_choice%cmat,    info)
     call mld_precset(prec,mld_coarse_fillin_,   prec_choice%cfill,   info)
-    call mld_precset(prec,mld_coarse_iluthrs_,    prec_choice%cthres,  info)
+    call mld_precset(prec,mld_coarse_iluthrs_,  prec_choice%cthres,  info)
     call mld_precset(prec,mld_coarse_sweeps_,   prec_choice%cjswp,   info)
   end if
 
diff --git a/tests/pdegen/ppde.f90 b/tests/pdegen/ppde.f90
index 24f043d9..7831170a 100644
--- a/tests/pdegen/ppde.f90
+++ b/tests/pdegen/ppde.f90
@@ -189,7 +189,7 @@ program ppde
     call mld_precset(prec,mld_coarse_subsolve_, prectype%csbsolve,info)
     call mld_precset(prec,mld_coarse_mat_,      prectype%cmat,    info)
     call mld_precset(prec,mld_coarse_fillin_,   prectype%cfill,   info)
-    call mld_precset(prec,mld_coarse_iluthrs_,    prectype%cthres,  info)
+    call mld_precset(prec,mld_coarse_iluthrs_,  prectype%cthres,  info)
     call mld_precset(prec,mld_coarse_sweeps_,   prectype%cjswp,   info)
   end if
   
@@ -356,7 +356,6 @@ contains
     return
 
   end subroutine get_parms
-
   !
   !  print an error message 
   !  
@@ -421,7 +420,6 @@ contains
     integer                  :: x,y,z,ia,indx_owner
     integer                  :: np, iam
     integer                  :: element
-    integer                  :: nv, inv
     integer, allocatable     :: irow(:),icol(:)
     real(psb_dpk_), allocatable :: val(:)
     integer, allocatable     :: prv(:)
@@ -454,11 +452,11 @@ contains
     if(iam == psb_root_) write(0,'("Generating Matrix (size=",i0x,")...")')n
 
     call psb_cdall(ictxt,desc_a,info,mg=n,parts=parts)
-    call psb_spall(a,desc_a,info,nnz=nnz)
+    if (info == 0) call psb_spall(a,desc_a,info,nnz=nnz)
     ! define  rhs from boundary conditions; also build initial guess 
-    call psb_geall(b,desc_a,info)
-    call psb_geall(xv,desc_a,info)
-    if(info /= 0) then
+    if (info == 0) call psb_geall(b,desc_a,info)
+    if (info == 0) call psb_geall(xv,desc_a,info)
+    if (info /= 0) then
       info=4010
       ch_err='allocation rout.'
       call psb_errpush(info,name,a_err=ch_err)
@@ -470,7 +468,7 @@ contains
     ! a bunch of rows per call. 
     ! 
     allocate(val(20*nbmax),irow(20*nbmax),&
-         &icol(20*nbmax),prv(np),stat=info)
+         &icol(20*nbmax),stat=info)
     if (info /= 0 ) then 
       info=4000
       call psb_errpush(info,name)
@@ -486,138 +484,135 @@ contains
 
     !    icol(1)=1    
     do glob_row = 1, n
-      call parts(glob_row,n,np,prv,nv)
-      do inv = 1, nv
-        indx_owner = prv(inv)
-        if (indx_owner == iam) then
-          ! local matrix pointer 
-          element=1
-          ! compute gridpoint coordinates
-          if (mod(glob_row,(idim*idim)) == 0) then
-            x = glob_row/(idim*idim)
-          else
-            x = glob_row/(idim*idim)+1
-          endif
-          if (mod((glob_row-(x-1)*idim*idim),idim) == 0) then
-            y = (glob_row-(x-1)*idim*idim)/idim
-          else
-            y = (glob_row-(x-1)*idim*idim)/idim+1
-          endif
-          z = glob_row-(x-1)*idim*idim-(y-1)*idim
-          ! glob_x, glob_y, glob_x coordinates
-          glob_x=x*deltah
-          glob_y=y*deltah
-          glob_z=z*deltah
-
-          ! check on boundary points 
-          zt(1) = 0.d0
-          ! internal point: build discretization
-          !   
-          !  term depending on   (x-1,y,z)
-          !
-          if (x==1) then 
-            val(element)=-b1(glob_x,glob_y,glob_z)&
-                 & -a1(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            zt(1) = exp(-glob_y**2-glob_z**2)*(-val(element))
-          else
-            val(element)=-b1(glob_x,glob_y,glob_z)&
-                 & -a1(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            icol(element)=(x-2)*idim*idim+(y-1)*idim+(z)
-            element=element+1
-          endif
-          !  term depending on     (x,y-1,z)
-          if (y==1) then 
-            val(element)=-b2(glob_x,glob_y,glob_z)&
-                 & -a2(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))  
-          else
-            val(element)=-b2(glob_x,glob_y,glob_z)&
-                 & -a2(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            icol(element)=(x-1)*idim*idim+(y-2)*idim+(z)
-            element=element+1
-          endif
-          !  term depending on     (x,y,z-1)
-          if (z==1) then 
-            val(element)=-b3(glob_x,glob_y,glob_z)&
-                 & -a3(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))  
-          else
-            val(element)=-b3(glob_x,glob_y,glob_z)&
-                 & -a3(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            icol(element)=(x-1)*idim*idim+(y-1)*idim+(z-1)
-            element=element+1
-          endif
-          !  term depending on     (x,y,z)
-          val(element)=2*b1(glob_x,glob_y,glob_z)&
-               & +2*b2(glob_x,glob_y,glob_z)&
-               & +2*b3(glob_x,glob_y,glob_z)&
-               & +a1(glob_x,glob_y,glob_z)&
-               & +a2(glob_x,glob_y,glob_z)&
-               & +a3(glob_x,glob_y,glob_z)
+      ! Figure out which rows are local  to the current process:
+      if (psb_is_owned(glob_row,desc_a)) then
+        ! local matrix pointer 
+        element=1
+        ! compute gridpoint coordinates
+        if (mod(glob_row,(idim*idim)) == 0) then
+          x = glob_row/(idim*idim)
+        else
+          x = glob_row/(idim*idim)+1
+        endif
+        if (mod((glob_row-(x-1)*idim*idim),idim) == 0) then
+          y = (glob_row-(x-1)*idim*idim)/idim
+        else
+          y = (glob_row-(x-1)*idim*idim)/idim+1
+        endif
+        z = glob_row-(x-1)*idim*idim-(y-1)*idim
+        ! glob_x, glob_y, glob_x coordinates
+        glob_x=x*deltah
+        glob_y=y*deltah
+        glob_z=z*deltah
+
+        ! check on boundary points 
+        zt(1) = 0.d0
+        ! internal point: build discretization
+        !   
+        !  term depending on   (x-1,y,z)
+        !
+        if (x==1) then 
+          val(element)=-b1(glob_x,glob_y,glob_z)&
+               & -a1(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          zt(1) = exp(-glob_y**2-glob_z**2)*(-val(element))
+        else
+          val(element)=-b1(glob_x,glob_y,glob_z)&
+               & -a1(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          icol(element)=(x-2)*idim*idim+(y-1)*idim+(z)
+          element=element+1
+        endif
+        !  term depending on     (x,y-1,z)
+        if (y==1) then 
+          val(element)=-b2(glob_x,glob_y,glob_z)&
+               & -a2(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))  
+        else
+          val(element)=-b2(glob_x,glob_y,glob_z)&
+               & -a2(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          icol(element)=(x-1)*idim*idim+(y-2)*idim+(z)
+          element=element+1
+        endif
+        !  term depending on     (x,y,z-1)
+        if (z==1) then 
+          val(element)=-b3(glob_x,glob_y,glob_z)&
+               & -a3(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))  
+        else
+          val(element)=-b3(glob_x,glob_y,glob_z)&
+               & -a3(glob_x,glob_y,glob_z)
           val(element) = val(element)/(deltah*&
                & deltah)
-          icol(element)=(x-1)*idim*idim+(y-1)*idim+(z)
-          element=element+1                  
-          !  term depending on     (x,y,z+1)
-          if (z==idim) then 
-            val(element)=-b1(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))  
-          else
-            val(element)=-b1(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            icol(element)=(x-1)*idim*idim+(y-1)*idim+(z+1)
-            element=element+1
-          endif
-          !  term depending on     (x,y+1,z)
-          if (y==idim) then 
-            val(element)=-b2(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))  
-          else
-            val(element)=-b2(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            icol(element)=(x-1)*idim*idim+(y)*idim+(z)
-            element=element+1
-          endif
-          !  term depending on     (x+1,y,z)
-          if (x<idim) then 
-            val(element)=-b3(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            icol(element)=(x)*idim*idim+(y-1)*idim+(z)
-            element=element+1
-          endif
-          irow(1:element-1)=glob_row
-          ia=glob_row
-
-          t3 = psb_wtime()
-          call psb_spins(element-1,irow,icol,val,a,desc_a,info)
-          if(info /= 0) exit
-          tins = tins + (psb_wtime()-t3)
-          call psb_geins(1,(/ia/),zt(1:1),b,desc_a,info)
-          if(info /= 0) exit
-          zt(1)=0.d0
-          call psb_geins(1,(/ia/),zt(1:1),xv,desc_a,info)
-          if(info /= 0) exit
-        end if
-      end do
+          icol(element)=(x-1)*idim*idim+(y-1)*idim+(z-1)
+          element=element+1
+        endif
+        !  term depending on     (x,y,z)
+        val(element)=2*b1(glob_x,glob_y,glob_z)&
+             & +2*b2(glob_x,glob_y,glob_z)&
+             & +2*b3(glob_x,glob_y,glob_z)&
+             & +a1(glob_x,glob_y,glob_z)&
+             & +a2(glob_x,glob_y,glob_z)&
+             & +a3(glob_x,glob_y,glob_z)
+        val(element) = val(element)/(deltah*&
+             & deltah)
+        icol(element)=(x-1)*idim*idim+(y-1)*idim+(z)
+        element=element+1                  
+        !  term depending on     (x,y,z+1)
+        if (z==idim) then 
+          val(element)=-b1(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))  
+        else
+          val(element)=-b1(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          icol(element)=(x-1)*idim*idim+(y-1)*idim+(z+1)
+          element=element+1
+        endif
+        !  term depending on     (x,y+1,z)
+        if (y==idim) then 
+          val(element)=-b2(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))  
+        else
+          val(element)=-b2(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          icol(element)=(x-1)*idim*idim+(y)*idim+(z)
+          element=element+1
+        endif
+        !  term depending on     (x+1,y,z)
+        if (x<idim) then 
+          val(element)=-b3(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          icol(element)=(x)*idim*idim+(y-1)*idim+(z)
+          element=element+1
+        endif
+        irow(1:element-1)=glob_row
+        ia=glob_row
+
+        t3 = psb_wtime()
+        call psb_spins(element-1,irow,icol,val,a,desc_a,info)
+        if(info /= 0) exit
+        tins = tins + (psb_wtime()-t3)
+        call psb_geins(1,(/ia/),zt(1:1),b,desc_a,info)
+        if(info /= 0) exit
+        zt(1)=0.d0
+        call psb_geins(1,(/ia/),zt(1:1),xv,desc_a,info)
+        if(info /= 0) exit
+      end if
     end do
 
     call psb_barrier(ictxt)    
@@ -649,7 +644,7 @@ contains
     call psb_amx(ictxt,tasb)
 
     if(iam == psb_root_) then
-      write(*,'("The matrix has been generated and assembeld in ",a3," format.")')&
+      write(*,'("The matrix has been generated and assembled in ",a3," format.")')&
            &   a%fida(1:3)
       write(*,'("-pspins time   : ",es10.4)')tins
       write(*,'("-insert time   : ",es10.4)')t2
diff --git a/tests/pdegen/spde.f90 b/tests/pdegen/spde.f90
index b6d203e0..82474467 100644
--- a/tests/pdegen/spde.f90
+++ b/tests/pdegen/spde.f90
@@ -189,7 +189,7 @@ program spde
     call mld_precset(prec,mld_coarse_subsolve_, prectype%csbsolve,info)
     call mld_precset(prec,mld_coarse_mat_,      prectype%cmat,    info)
     call mld_precset(prec,mld_coarse_fillin_,   prectype%cfill,   info)
-    call mld_precset(prec,mld_coarse_iluthrs_,    prectype%cthres,  info)
+    call mld_precset(prec,mld_coarse_iluthrs_,  prectype%cthres,  info)
     call mld_precset(prec,mld_coarse_sweeps_,   prectype%cjswp,   info)
   end if
   
@@ -356,7 +356,6 @@ contains
     return
 
   end subroutine get_parms
-
   !
   !  print an error message 
   !  
@@ -421,7 +420,6 @@ contains
     integer                  :: x,y,z,ia,indx_owner
     integer                  :: np, iam
     integer                  :: element
-    integer                  :: nv, inv
     integer, allocatable     :: irow(:),icol(:)
     real(psb_spk_), allocatable :: val(:)
     integer, allocatable     :: prv(:)
@@ -454,11 +452,11 @@ contains
     if(iam == psb_root_) write(0,'("Generating Matrix (size=",i0x,")...")')n
 
     call psb_cdall(ictxt,desc_a,info,mg=n,parts=parts)
-    call psb_spall(a,desc_a,info,nnz=nnz)
+    if (info == 0) call psb_spall(a,desc_a,info,nnz=nnz)
     ! define  rhs from boundary conditions; also build initial guess 
-    call psb_geall(b,desc_a,info)
-    call psb_geall(xv,desc_a,info)
-    if(info /= 0) then
+    if (info == 0) call psb_geall(b,desc_a,info)
+    if (info == 0) call psb_geall(xv,desc_a,info)
+    if (info /= 0) then
       info=4010
       ch_err='allocation rout.'
       call psb_errpush(info,name,a_err=ch_err)
@@ -470,7 +468,7 @@ contains
     ! a bunch of rows per call. 
     ! 
     allocate(val(20*nbmax),irow(20*nbmax),&
-         &icol(20*nbmax),prv(np),stat=info)
+         &icol(20*nbmax),stat=info)
     if (info /= 0 ) then 
       info=4000
       call psb_errpush(info,name)
@@ -486,138 +484,135 @@ contains
 
     !    icol(1)=1    
     do glob_row = 1, n
-      call parts(glob_row,n,np,prv,nv)
-      do inv = 1, nv
-        indx_owner = prv(inv)
-        if (indx_owner == iam) then
-          ! local matrix pointer 
-          element=1
-          ! compute gridpoint coordinates
-          if (mod(glob_row,(idim*idim)) == 0) then
-            x = glob_row/(idim*idim)
-          else
-            x = glob_row/(idim*idim)+1
-          endif
-          if (mod((glob_row-(x-1)*idim*idim),idim) == 0) then
-            y = (glob_row-(x-1)*idim*idim)/idim
-          else
-            y = (glob_row-(x-1)*idim*idim)/idim+1
-          endif
-          z = glob_row-(x-1)*idim*idim-(y-1)*idim
-          ! glob_x, glob_y, glob_x coordinates
-          glob_x=x*deltah
-          glob_y=y*deltah
-          glob_z=z*deltah
-
-          ! check on boundary points 
-          zt(1) = 0.d0
-          ! internal point: build discretization
-          !   
-          !  term depending on   (x-1,y,z)
-          !
-          if (x==1) then 
-            val(element)=-b1(glob_x,glob_y,glob_z)&
-                 & -a1(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            zt(1) = exp(-glob_y**2-glob_z**2)*(-val(element))
-          else
-            val(element)=-b1(glob_x,glob_y,glob_z)&
-                 & -a1(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            icol(element)=(x-2)*idim*idim+(y-1)*idim+(z)
-            element=element+1
-          endif
-          !  term depending on     (x,y-1,z)
-          if (y==1) then 
-            val(element)=-b2(glob_x,glob_y,glob_z)&
-                 & -a2(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))  
-          else
-            val(element)=-b2(glob_x,glob_y,glob_z)&
-                 & -a2(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            icol(element)=(x-1)*idim*idim+(y-2)*idim+(z)
-            element=element+1
-          endif
-          !  term depending on     (x,y,z-1)
-          if (z==1) then 
-            val(element)=-b3(glob_x,glob_y,glob_z)&
-                 & -a3(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))  
-          else
-            val(element)=-b3(glob_x,glob_y,glob_z)&
-                 & -a3(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            icol(element)=(x-1)*idim*idim+(y-1)*idim+(z-1)
-            element=element+1
-          endif
-          !  term depending on     (x,y,z)
-          val(element)=2*b1(glob_x,glob_y,glob_z)&
-               & +2*b2(glob_x,glob_y,glob_z)&
-               & +2*b3(glob_x,glob_y,glob_z)&
-               & +a1(glob_x,glob_y,glob_z)&
-               & +a2(glob_x,glob_y,glob_z)&
-               & +a3(glob_x,glob_y,glob_z)
+      ! Figure out which rows are local  to the current process:
+      if (psb_is_owned(glob_row,desc_a)) then
+        ! local matrix pointer 
+        element=1
+        ! compute gridpoint coordinates
+        if (mod(glob_row,(idim*idim)) == 0) then
+          x = glob_row/(idim*idim)
+        else
+          x = glob_row/(idim*idim)+1
+        endif
+        if (mod((glob_row-(x-1)*idim*idim),idim) == 0) then
+          y = (glob_row-(x-1)*idim*idim)/idim
+        else
+          y = (glob_row-(x-1)*idim*idim)/idim+1
+        endif
+        z = glob_row-(x-1)*idim*idim-(y-1)*idim
+        ! glob_x, glob_y, glob_x coordinates
+        glob_x=x*deltah
+        glob_y=y*deltah
+        glob_z=z*deltah
+
+        ! check on boundary points 
+        zt(1) = 0.d0
+        ! internal point: build discretization
+        !   
+        !  term depending on   (x-1,y,z)
+        !
+        if (x==1) then 
+          val(element)=-b1(glob_x,glob_y,glob_z)&
+               & -a1(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          zt(1) = exp(-glob_y**2-glob_z**2)*(-val(element))
+        else
+          val(element)=-b1(glob_x,glob_y,glob_z)&
+               & -a1(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          icol(element)=(x-2)*idim*idim+(y-1)*idim+(z)
+          element=element+1
+        endif
+        !  term depending on     (x,y-1,z)
+        if (y==1) then 
+          val(element)=-b2(glob_x,glob_y,glob_z)&
+               & -a2(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))  
+        else
+          val(element)=-b2(glob_x,glob_y,glob_z)&
+               & -a2(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          icol(element)=(x-1)*idim*idim+(y-2)*idim+(z)
+          element=element+1
+        endif
+        !  term depending on     (x,y,z-1)
+        if (z==1) then 
+          val(element)=-b3(glob_x,glob_y,glob_z)&
+               & -a3(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))  
+        else
+          val(element)=-b3(glob_x,glob_y,glob_z)&
+               & -a3(glob_x,glob_y,glob_z)
           val(element) = val(element)/(deltah*&
                & deltah)
-          icol(element)=(x-1)*idim*idim+(y-1)*idim+(z)
-          element=element+1                  
-          !  term depending on     (x,y,z+1)
-          if (z==idim) then 
-            val(element)=-b1(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))  
-          else
-            val(element)=-b1(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            icol(element)=(x-1)*idim*idim+(y-1)*idim+(z+1)
-            element=element+1
-          endif
-          !  term depending on     (x,y+1,z)
-          if (y==idim) then 
-            val(element)=-b2(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))  
-          else
-            val(element)=-b2(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            icol(element)=(x-1)*idim*idim+(y)*idim+(z)
-            element=element+1
-          endif
-          !  term depending on     (x+1,y,z)
-          if (x<idim) then 
-            val(element)=-b3(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            icol(element)=(x)*idim*idim+(y-1)*idim+(z)
-            element=element+1
-          endif
-          irow(1:element-1)=glob_row
-          ia=glob_row
-
-          t3 = psb_wtime()
-          call psb_spins(element-1,irow,icol,val,a,desc_a,info)
-          if(info /= 0) exit
-          tins = tins + (psb_wtime()-t3)
-          call psb_geins(1,(/ia/),zt(1:1),b,desc_a,info)
-          if(info /= 0) exit
-          zt(1)=0.d0
-          call psb_geins(1,(/ia/),zt(1:1),xv,desc_a,info)
-          if(info /= 0) exit
-        end if
-      end do
+          icol(element)=(x-1)*idim*idim+(y-1)*idim+(z-1)
+          element=element+1
+        endif
+        !  term depending on     (x,y,z)
+        val(element)=2*b1(glob_x,glob_y,glob_z)&
+             & +2*b2(glob_x,glob_y,glob_z)&
+             & +2*b3(glob_x,glob_y,glob_z)&
+             & +a1(glob_x,glob_y,glob_z)&
+             & +a2(glob_x,glob_y,glob_z)&
+             & +a3(glob_x,glob_y,glob_z)
+        val(element) = val(element)/(deltah*&
+             & deltah)
+        icol(element)=(x-1)*idim*idim+(y-1)*idim+(z)
+        element=element+1                  
+        !  term depending on     (x,y,z+1)
+        if (z==idim) then 
+          val(element)=-b1(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))  
+        else
+          val(element)=-b1(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          icol(element)=(x-1)*idim*idim+(y-1)*idim+(z+1)
+          element=element+1
+        endif
+        !  term depending on     (x,y+1,z)
+        if (y==idim) then 
+          val(element)=-b2(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))  
+        else
+          val(element)=-b2(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          icol(element)=(x-1)*idim*idim+(y)*idim+(z)
+          element=element+1
+        endif
+        !  term depending on     (x+1,y,z)
+        if (x<idim) then 
+          val(element)=-b3(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          icol(element)=(x)*idim*idim+(y-1)*idim+(z)
+          element=element+1
+        endif
+        irow(1:element-1)=glob_row
+        ia=glob_row
+
+        t3 = psb_wtime()
+        call psb_spins(element-1,irow,icol,val,a,desc_a,info)
+        if(info /= 0) exit
+        tins = tins + (psb_wtime()-t3)
+        call psb_geins(1,(/ia/),zt(1:1),b,desc_a,info)
+        if(info /= 0) exit
+        zt(1)=0.d0
+        call psb_geins(1,(/ia/),zt(1:1),xv,desc_a,info)
+        if(info /= 0) exit
+      end if
     end do
 
     call psb_barrier(ictxt)    
@@ -649,7 +644,7 @@ contains
     call psb_amx(ictxt,tasb)
 
     if(iam == psb_root_) then
-      write(*,'("The matrix has been generated and assembeld in ",a3," format.")')&
+      write(*,'("The matrix has been generated and assembled in ",a3," format.")')&
            &   a%fida(1:3)
       write(*,'("-pspins time   : ",es10.4)')tins
       write(*,'("-insert time   : ",es10.4)')t2