diff --git a/examples/pdegen/mld_dexample_1lev.f90 b/examples/pdegen/mld_dexample_1lev.f90
index dbc75816..993ed972 100644
--- a/examples/pdegen/mld_dexample_1lev.f90
+++ b/examples/pdegen/mld_dexample_1lev.f90
@@ -81,32 +81,32 @@ program mld_dexample_1lev
   implicit none
 
 
-! sparse matrices
+  ! sparse matrices
   type(psb_dspmat_type) :: A
 
-! descriptor of sparse matrices
+  ! descriptor of sparse matrices
   type(psb_desc_type):: desc_A
 
-! preconditioner
+  ! preconditioner
   type(mld_dprec_type)  :: P
 
-! right-hand side, solution and residual vectors
+  ! right-hand side, solution and residual vectors
   real(psb_dpk_), allocatable , save  :: b(:), x(:), r(:)
 
-! solver parameters
+  ! solver parameters
   real(psb_dpk_)   :: tol, err
   integer :: itmax, iter, itrace, istop
 
-! parallel environment parameters
+  ! parallel environment parameters
   integer            :: ictxt, iam, np
 
-! other variables
+  ! other variables
   integer            :: i,info,j,amatsize,descsize,precsize
   integer            :: idim, nlev, ierr, ircode
   real(psb_dpk_) :: t1, t2, tprec, resmx, resmxp
   character(len=20)  :: name
 
-! initialize the parallel environment
+  ! initialize the parallel environment
   call psb_init(ictxt)
   call psb_info(ictxt,iam,np)
 
@@ -121,11 +121,11 @@ program mld_dexample_1lev
   info=0
   call psb_set_errverbosity(2)
 
-! get parameters
+  ! get parameters
 
   call get_parms(ictxt,idim,itmax,tol)
 
-! allocate and fill in the coefficient matrix, rhs and initial guess
+  ! allocate and fill in the coefficient matrix, rhs and initial guess
 
   call psb_barrier(ictxt)
   t1 = psb_wtime()
@@ -141,12 +141,12 @@ program mld_dexample_1lev
   if (iam == psb_root_) write(*,'("Overall matrix creation time : ",es10.4)')t2
   if (iam == psb_root_) write(*,'(" ")')
 
-! set RAS with overlap 2 and ILU(0) on the local blocks
+  ! set RAS with overlap 2 and ILU(0) on the local blocks
 
   call mld_precinit(P,'AS',info)
   call mld_precset(P,mld_sub_ovr_,2,info)
 
-! build the preconditioner
+  ! build the preconditioner
 
   call psb_barrier(ictxt)
   t1 = psb_wtime()
@@ -161,13 +161,13 @@ program mld_dexample_1lev
     goto 9999
   end if
 
-! set the initial guess
+  ! set the initial guess
 
   call psb_geall(x,desc_A,info)
   x(:) =0.0
   call psb_geasb(x,desc_A,info)
 
-! solve Ax=b with preconditioned BiCGSTAB
+  ! solve Ax=b with preconditioned BiCGSTAB
 
   call psb_barrier(ictxt)
   t1 = psb_wtime()
@@ -289,17 +289,15 @@ contains
         integer, intent(out) :: pv(*) 
       end subroutine parts
     end interface
-   ! local variables
+    ! local variables
     type(psb_dspmat_type)    :: a
     real(psb_dpk_)           :: zt(nbmax),glob_x,glob_y,glob_z
     integer                  :: m,n,nnz,glob_row,ipoints
     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(:)
     ! deltah dimension of each grid cell
     ! deltat discretization time
     real(psb_dpk_)         :: deltah
@@ -329,11 +327,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)
-    ! 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_spall(a,desc_a,info,nnz=nnz)
+    ! define  rhs from boundary conditions; also build initial guess 
+    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
       call psb_errpush(info,name)
       goto 9999
@@ -344,7 +342,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)
@@ -359,143 +357,140 @@ contains
     ! distribution.
 
     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,ipoints*ipoints) == 0) then
-            x = glob_row/(ipoints*ipoints)
-          else
-            x = glob_row/(ipoints*ipoints)+1
-          endif
-          if (mod((glob_row-(x-1)*ipoints*ipoints),ipoints) == 0) then
-            y = (glob_row-(x-1)*ipoints*ipoints)/ipoints
-          else
-            y = (glob_row-(x-1)*ipoints*ipoints)/ipoints+1
-          endif
-          z = glob_row-(x-1)*ipoints*ipoints-(y-1)*ipoints
-          ! 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)*ipoints*ipoints+(y-1)*ipoints+(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_x**2-glob_z**2)*(-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)*ipoints*ipoints+(y-2)*ipoints+(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_x**2-glob_y**2)*(-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)*ipoints*ipoints+(y-1)*ipoints+(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,ipoints*ipoints) == 0) then
+          x = glob_row/(ipoints*ipoints)
+        else
+          x = glob_row/(ipoints*ipoints)+1
+        endif
+        if (mod((glob_row-(x-1)*ipoints*ipoints),ipoints) == 0) then
+          y = (glob_row-(x-1)*ipoints*ipoints)/ipoints
+        else
+          y = (glob_row-(x-1)*ipoints*ipoints)/ipoints+1
+        endif
+        z = glob_row-(x-1)*ipoints*ipoints-(y-1)*ipoints
+        ! 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)*ipoints*ipoints+(y-1)*ipoints+(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_x**2-glob_z**2)*(-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)*ipoints*ipoints+(y-2)*ipoints+(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_x**2-glob_y**2)*(-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)*ipoints*ipoints+(y-1)*ipoints+(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)*ipoints*ipoints+(y-1)*ipoints+(z)
+        element=element+1                  
+        !  term depending on     (x,y,z+1)
+        if (z==ipoints) then 
+          val(element)=-b1(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          zt(1) = exp(-glob_x**2-glob_y**2)*exp(-glob_z)*(-val(element))  
+        else
+          val(element)=-b1(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z+1)
+          element=element+1
+        endif
+        !  term depending on     (x,y+1,z)
+        if (y==ipoints) then 
+          val(element)=-b2(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          zt(1) = exp(-glob_x**2-glob_z**2)*exp(-glob_y)*(-val(element))  
+        else
+          val(element)=-b2(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          icol(element)=(x-1)*ipoints*ipoints+(y)*ipoints+(z)
+          element=element+1
+        endif
+        !  term depending on     (x+1,y,z)
+        if (x==ipoints) then 
+          val(element)=-b3(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)
           val(element) = val(element)/(deltah*&
                & deltah)
-          icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z)
-          element=element+1                  
-          !  term depending on     (x,y,z+1)
-          if (z==ipoints) then 
-            val(element)=-b1(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            zt(1) = exp(-glob_x**2-glob_y**2)*exp(-glob_z)*(-val(element))  
-          else
-            val(element)=-b1(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z+1)
-            element=element+1
-          endif
-          !  term depending on     (x,y+1,z)
-          if (y==ipoints) then 
-            val(element)=-b2(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-           zt(1) = exp(-glob_x**2-glob_z**2)*exp(-glob_y)*(-val(element))  
-          else
-            val(element)=-b2(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            icol(element)=(x-1)*ipoints*ipoints+(y)*ipoints+(z)
-            element=element+1
-          endif
-          !  term depending on     (x+1,y,z)
-          if (x==ipoints) then 
-            val(element)=-b3(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)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            icol(element)=(x)*ipoints*ipoints+(y-1)*ipoints+(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)*ipoints*ipoints+(y-1)*ipoints+(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)    
@@ -525,7 +520,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/examples/pdegen/mld_dexample_ml.f90 b/examples/pdegen/mld_dexample_ml.f90
index 4cd836e5..09d8c59b 100644
--- a/examples/pdegen/mld_dexample_ml.f90
+++ b/examples/pdegen/mld_dexample_ml.f90
@@ -80,36 +80,36 @@ program mld_dexample_ml
 
   implicit none
 
-! input parameters
+  ! input parameters
 
-! sparse matrices
+  ! sparse matrices
   type(psb_dspmat_type) :: A
 
-! sparse matrices descriptor
+  ! sparse matrices descriptor
   type(psb_desc_type):: desc_A
 
-! preconditioner
+  ! preconditioner
   type(mld_dprec_type)  :: P
 
-! right-hand side, solution and residual vectors
+  ! right-hand side, solution and residual vectors
   real(psb_dpk_), allocatable , save  :: b(:), x(:), r(:)
 
-! solver and preconditioner parameters
+  ! solver and preconditioner parameters
   real(psb_dpk_)   :: tol, err
   integer          :: itmax, iter, istop
   integer          :: nlev
 
-! parallel environment parameters
+  ! parallel environment parameters
   integer            :: ictxt, iam, np
 
-! other variables
+  ! other variables
   integer            :: choice       
   integer            :: i,info,j,amatsize,descsize,precsize
   integer            :: idim, ierr, ircode
   real(psb_dpk_)     :: t1, t2, tprec, resmx, resmxp
   character(len=20)  :: name
 
-! initialize the parallel environment
+  ! initialize the parallel environment
 
   call psb_init(ictxt)
   call psb_info(ictxt,iam,np)
@@ -125,11 +125,11 @@ program mld_dexample_ml
   info=0
   call psb_set_errverbosity(2)
 
-! get parameters
+  ! get parameters
 
   call get_parms(ictxt,choice,idim,itmax,tol)
 
-!  allocate and fill in the coefficient matrix, rhs and initial guess
+  !  allocate and fill in the coefficient matrix, rhs and initial guess
 
   call psb_barrier(ictxt)
   t1 = psb_wtime()
@@ -149,40 +149,40 @@ program mld_dexample_ml
 
   case(1)
 
-! initialize the default multi-level preconditioner, i.e. hybrid
-! Schwarz, using RAS (with overlap 1 and ILU(0) on the blocks)
-! as post-smoother and 4 block-Jacobi sweeps (with UMFPACK LU
-! on the blocks) as distributed coarse-level solver
+    ! initialize the default multi-level preconditioner, i.e. hybrid
+    ! Schwarz, using RAS (with overlap 1 and ILU(0) on the blocks)
+    ! as post-smoother and 4 block-Jacobi sweeps (with UMFPACK LU
+    ! on the blocks) as distributed coarse-level solver
 
-  call mld_precinit(P,'ML',info)
+    call mld_precinit(P,'ML',info)
 
   case(2)
 
-! set a three-level hybrid Schwarz preconditioner, which uses
-! block Jacobi (with ILU(0) on the blocks) as post-smoother,
-! a coarsest matrix replicated on the processors, and the
-! LU factorization from UMFPACK as coarse-level solver
+    ! set a three-level hybrid Schwarz preconditioner, which uses
+    ! block Jacobi (with ILU(0) on the blocks) as post-smoother,
+    ! a coarsest matrix replicated on the processors, and the
+    ! LU factorization from UMFPACK as coarse-level solver
 
-  call mld_precinit(P,'ML',info,nlev=3)
-  call mld_precset(P,mld_smoother_type_,'BJAC',info)
-  call mld_precset(P,mld_coarse_mat_,'REPL',info)
-  call mld_precset(P,mld_coarse_solve_,'UMF',info)
+    call mld_precinit(P,'ML',info,nlev=3)
+    call mld_precset(P,mld_smoother_type_,'BJAC',info)
+    call mld_precset(P,mld_coarse_mat_,'REPL',info)
+    call mld_precset(P,mld_coarse_solve_,'UMF',info)
 
   case(3)
 
-! set a three-level additive Schwarz preconditioner, which uses
-! RAS (with overlap 1 and ILU(0) on the blocks) as pre- and
-! post-smoother, and 5 block-Jacobi sweeps (with UMFPACK LU
-! on the blocks) as distributed coarsest-level solver
+    ! set a three-level additive Schwarz preconditioner, which uses
+    ! RAS (with overlap 1 and ILU(0) on the blocks) as pre- and
+    ! post-smoother, and 5 block-Jacobi sweeps (with UMFPACK LU
+    ! on the blocks) as distributed coarsest-level solver
 
-  call mld_precinit(P,'ML',info,nlev=3)
-  call mld_precset(P,mld_ml_type_,'ADD',info)
-  call mld_precset(P,mld_smoother_pos_,'TWOSIDE',info)
-  call mld_precset(P,mld_coarse_sweeps_,5,info)
+    call mld_precinit(P,'ML',info,nlev=3)
+    call mld_precset(P,mld_ml_type_,'ADD',info)
+    call mld_precset(P,mld_smoother_pos_,'TWOSIDE',info)
+    call mld_precset(P,mld_coarse_sweeps_,5,info)
 
   end select
 
-! build the preconditioner
+  ! build the preconditioner
 
   call psb_barrier(ictxt)
   t1 = psb_wtime()
@@ -197,13 +197,13 @@ program mld_dexample_ml
     goto 9999
   end if
 
-! set the solver parameters and the initial guess
+  ! set the solver parameters and the initial guess
 
   call psb_geall(x,desc_A,info)
   x(:) =0.0
   call psb_geasb(x,desc_A,info)
 
-! solve Ax=b with preconditioned BiCGSTAB
+  ! solve Ax=b with preconditioned BiCGSTAB
 
   call psb_barrier(ictxt)
   t1 = psb_wtime()
@@ -327,17 +327,15 @@ contains
         integer, intent(out) :: pv(*) 
       end subroutine parts
     end interface
-   ! local variables
+    ! local variables
     type(psb_dspmat_type)    :: a
     real(psb_dpk_)           :: zt(nbmax),glob_x,glob_y,glob_z
     integer                  :: m,n,nnz,glob_row,ipoints
     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(:)
     ! deltah dimension of each grid cell
     ! deltat discretization time
     real(psb_dpk_)         :: deltah
@@ -367,11 +365,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)
-    ! 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_spall(a,desc_a,info,nnz=nnz)
+    ! define  rhs from boundary conditions; also build initial guess 
+    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
       call psb_errpush(info,name)
       goto 9999
@@ -382,7 +380,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)
@@ -397,143 +395,140 @@ contains
     ! distribution.
 
     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,ipoints*ipoints) == 0) then
-            x = glob_row/(ipoints*ipoints)
-          else
-            x = glob_row/(ipoints*ipoints)+1
-          endif
-          if (mod((glob_row-(x-1)*ipoints*ipoints),ipoints) == 0) then
-            y = (glob_row-(x-1)*ipoints*ipoints)/ipoints
-          else
-            y = (glob_row-(x-1)*ipoints*ipoints)/ipoints+1
-          endif
-          z = glob_row-(x-1)*ipoints*ipoints-(y-1)*ipoints
-          ! 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)*ipoints*ipoints+(y-1)*ipoints+(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_x**2-glob_z**2)*(-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)*ipoints*ipoints+(y-2)*ipoints+(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_x**2-glob_y**2)*(-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)*ipoints*ipoints+(y-1)*ipoints+(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,ipoints*ipoints) == 0) then
+          x = glob_row/(ipoints*ipoints)
+        else
+          x = glob_row/(ipoints*ipoints)+1
+        endif
+        if (mod((glob_row-(x-1)*ipoints*ipoints),ipoints) == 0) then
+          y = (glob_row-(x-1)*ipoints*ipoints)/ipoints
+        else
+          y = (glob_row-(x-1)*ipoints*ipoints)/ipoints+1
+        endif
+        z = glob_row-(x-1)*ipoints*ipoints-(y-1)*ipoints
+        ! 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)*ipoints*ipoints+(y-1)*ipoints+(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_x**2-glob_z**2)*(-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)*ipoints*ipoints+(y-2)*ipoints+(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_x**2-glob_y**2)*(-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)*ipoints*ipoints+(y-1)*ipoints+(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)*ipoints*ipoints+(y-1)*ipoints+(z)
+        element=element+1                  
+        !  term depending on     (x,y,z+1)
+        if (z==ipoints) then 
+          val(element)=-b1(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          zt(1) = exp(-glob_x**2-glob_y**2)*exp(-glob_z)*(-val(element))  
+        else
+          val(element)=-b1(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z+1)
+          element=element+1
+        endif
+        !  term depending on     (x,y+1,z)
+        if (y==ipoints) then 
+          val(element)=-b2(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          zt(1) = exp(-glob_x**2-glob_z**2)*exp(-glob_y)*(-val(element))  
+        else
+          val(element)=-b2(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          icol(element)=(x-1)*ipoints*ipoints+(y)*ipoints+(z)
+          element=element+1
+        endif
+        !  term depending on     (x+1,y,z)
+        if (x==ipoints) then 
+          val(element)=-b3(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)
           val(element) = val(element)/(deltah*&
                & deltah)
-          icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z)
-          element=element+1                  
-          !  term depending on     (x,y,z+1)
-          if (z==ipoints) then 
-            val(element)=-b1(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            zt(1) = exp(-glob_x**2-glob_y**2)*exp(-glob_z)*(-val(element))  
-          else
-            val(element)=-b1(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z+1)
-            element=element+1
-          endif
-          !  term depending on     (x,y+1,z)
-          if (y==ipoints) then 
-            val(element)=-b2(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-           zt(1) = exp(-glob_x**2-glob_z**2)*exp(-glob_y)*(-val(element))  
-          else
-            val(element)=-b2(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            icol(element)=(x-1)*ipoints*ipoints+(y)*ipoints+(z)
-            element=element+1
-          endif
-          !  term depending on     (x+1,y,z)
-          if (x==ipoints) then 
-            val(element)=-b3(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)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            icol(element)=(x)*ipoints*ipoints+(y-1)*ipoints+(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)*ipoints*ipoints+(y-1)*ipoints+(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)    
@@ -563,7 +558,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/examples/pdegen/mld_sexample_1lev.f90 b/examples/pdegen/mld_sexample_1lev.f90
index f3a1241a..8c449bae 100644
--- a/examples/pdegen/mld_sexample_1lev.f90
+++ b/examples/pdegen/mld_sexample_1lev.f90
@@ -300,7 +300,6 @@ contains
     integer                  :: nv, inv
     integer, allocatable     :: irow(:),icol(:)
     real(psb_spk_), allocatable :: val(:)
-    integer, allocatable     :: prv(:)
     ! deltah dimension of each grid cell
     ! deltat discretization time
     real(psb_spk_)         :: deltah
@@ -330,11 +329,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)
-    ! 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_spall(a,desc_a,info,nnz=nnz)
+    ! define  rhs from boundary conditions; also build initial guess 
+    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
       call psb_errpush(info,name)
       goto 9999
@@ -345,7 +344,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)
@@ -360,143 +359,140 @@ contains
     ! distribution.
 
     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,ipoints*ipoints) == 0) then
-            x = glob_row/(ipoints*ipoints)
-          else
-            x = glob_row/(ipoints*ipoints)+1
-          endif
-          if (mod((glob_row-(x-1)*ipoints*ipoints),ipoints) == 0) then
-            y = (glob_row-(x-1)*ipoints*ipoints)/ipoints
-          else
-            y = (glob_row-(x-1)*ipoints*ipoints)/ipoints+1
-          endif
-          z = glob_row-(x-1)*ipoints*ipoints-(y-1)*ipoints
-          ! 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)*ipoints*ipoints+(y-1)*ipoints+(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_x**2-glob_z**2)*(-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)*ipoints*ipoints+(y-2)*ipoints+(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_x**2-glob_y**2)*(-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)*ipoints*ipoints+(y-1)*ipoints+(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,ipoints*ipoints) == 0) then
+          x = glob_row/(ipoints*ipoints)
+        else
+          x = glob_row/(ipoints*ipoints)+1
+        endif
+        if (mod((glob_row-(x-1)*ipoints*ipoints),ipoints) == 0) then
+          y = (glob_row-(x-1)*ipoints*ipoints)/ipoints
+        else
+          y = (glob_row-(x-1)*ipoints*ipoints)/ipoints+1
+        endif
+        z = glob_row-(x-1)*ipoints*ipoints-(y-1)*ipoints
+        ! 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)*ipoints*ipoints+(y-1)*ipoints+(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_x**2-glob_z**2)*(-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)*ipoints*ipoints+(y-2)*ipoints+(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_x**2-glob_y**2)*(-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)*ipoints*ipoints+(y-1)*ipoints+(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)*ipoints*ipoints+(y-1)*ipoints+(z)
+        element=element+1                  
+        !  term depending on     (x,y,z+1)
+        if (z==ipoints) then 
+          val(element)=-b1(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          zt(1) = exp(-glob_x**2-glob_y**2)*exp(-glob_z)*(-val(element))  
+        else
+          val(element)=-b1(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z+1)
+          element=element+1
+        endif
+        !  term depending on     (x,y+1,z)
+        if (y==ipoints) then 
+          val(element)=-b2(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          zt(1) = exp(-glob_x**2-glob_z**2)*exp(-glob_y)*(-val(element))  
+        else
+          val(element)=-b2(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          icol(element)=(x-1)*ipoints*ipoints+(y)*ipoints+(z)
+          element=element+1
+        endif
+        !  term depending on     (x+1,y,z)
+        if (x==ipoints) then 
+          val(element)=-b3(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)
           val(element) = val(element)/(deltah*&
                & deltah)
-          icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z)
-          element=element+1                  
-          !  term depending on     (x,y,z+1)
-          if (z==ipoints) then 
-            val(element)=-b1(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            zt(1) = exp(-glob_x**2-glob_y**2)*exp(-glob_z)*(-val(element))  
-          else
-            val(element)=-b1(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z+1)
-            element=element+1
-          endif
-          !  term depending on     (x,y+1,z)
-          if (y==ipoints) then 
-            val(element)=-b2(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            zt(1) = exp(-glob_x**2-glob_z**2)*exp(-glob_y)*(-val(element))  
-          else
-            val(element)=-b2(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            icol(element)=(x-1)*ipoints*ipoints+(y)*ipoints+(z)
-            element=element+1
-          endif
-          !  term depending on     (x+1,y,z)
-          if (x==ipoints) then 
-            val(element)=-b3(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)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            icol(element)=(x)*ipoints*ipoints+(y-1)*ipoints+(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)*ipoints*ipoints+(y-1)*ipoints+(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)    
@@ -526,7 +522,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/examples/pdegen/mld_sexample_ml.f90 b/examples/pdegen/mld_sexample_ml.f90
index 6f91dab0..a96c3368 100644
--- a/examples/pdegen/mld_sexample_ml.f90
+++ b/examples/pdegen/mld_sexample_ml.f90
@@ -80,29 +80,29 @@ program mld_sexample_ml
 
   implicit none
 
-! input parameters
+  ! input parameters
 
-! sparse matrices
+  ! sparse matrices
   type(psb_sspmat_type) :: A
 
-! sparse matrices descriptor
+  ! sparse matrices descriptor
   type(psb_desc_type):: desc_A
 
-! preconditioner
+  ! preconditioner
   type(mld_sprec_type)  :: P
 
-! right-hand side, solution and residual vectors
+  ! right-hand side, solution and residual vectors
   real(psb_spk_), allocatable , save  :: b(:), x(:), r(:)
 
-! solver and preconditioner parameters
+  ! solver and preconditioner parameters
   real(psb_spk_)   :: tol, err
   integer          :: itmax, iter, istop
   integer          :: nlev
 
-! parallel environment parameters
+  ! parallel environment parameters
   integer            :: ictxt, iam, np
 
-! other variables
+  ! other variables
   integer            :: choice       
   integer            :: i,info,j,amatsize,descsize,precsize
   integer            :: idim, ierr, ircode
@@ -110,7 +110,7 @@ program mld_sexample_ml
   real(psb_spk_)     :: resmx, resmxp
   character(len=20)  :: name
 
-! initialize the parallel environment
+  ! initialize the parallel environment
 
   call psb_init(ictxt)
   call psb_info(ictxt,iam,np)
@@ -126,11 +126,11 @@ program mld_sexample_ml
   info=0
   call psb_set_errverbosity(2)
 
-! get parameters
+  ! get parameters
 
   call get_parms(ictxt,choice,idim,itmax,tol)
 
-!  allocate and fill in the coefficient matrix, rhs and initial guess
+  !  allocate and fill in the coefficient matrix, rhs and initial guess
 
   call psb_barrier(ictxt)
   t1 = psb_wtime()
@@ -150,40 +150,40 @@ program mld_sexample_ml
 
   case(1)
 
-! initialize the default multi-level preconditioner, i.e. hybrid
-! Schwarz, using RAS (with overlap 1 and ILU(0) on the blocks)
-! as post-smoother and 4 block-Jacobi sweeps (with UMFPACK LU
-! on the blocks) as distributed coarse-level solver
+    ! initialize the default multi-level preconditioner, i.e. hybrid
+    ! Schwarz, using RAS (with overlap 1 and ILU(0) on the blocks)
+    ! as post-smoother and 4 block-Jacobi sweeps (with UMFPACK LU
+    ! on the blocks) as distributed coarse-level solver
 
-  call mld_precinit(P,'ML',info)
+    call mld_precinit(P,'ML',info)
 
   case(2)
 
-! set a three-level hybrid Schwarz preconditioner, which uses
-! block Jacobi (with ILU(0) on the blocks) as post-smoother,
-! a coarsest matrix replicated on the processors, and the
-! LU factorization from UMFPACK as coarse-level solver
+    ! set a three-level hybrid Schwarz preconditioner, which uses
+    ! block Jacobi (with ILU(0) on the blocks) as post-smoother,
+    ! a coarsest matrix replicated on the processors, and the
+    ! LU factorization from UMFPACK as coarse-level solver
 
-  call mld_precinit(P,'ML',info,nlev=3)
-  call mld_precset(P,mld_smoother_type_,'BJAC',info)
-  call mld_precset(P,mld_coarse_mat_,'REPL',info)
-  call mld_precset(P,mld_coarse_solve_,'UMF',info)
+    call mld_precinit(P,'ML',info,nlev=3)
+    call mld_precset(P,mld_smoother_type_,'BJAC',info)
+    call mld_precset(P,mld_coarse_mat_,'REPL',info)
+    call mld_precset(P,mld_coarse_solve_,'UMF',info)
 
   case(3)
 
-! set a three-level additive Schwarz preconditioner, which uses
-! RAS (with overlap 1 and ILU(0) on the blocks) as pre- and
-! post-smoother, and 5 block-Jacobi sweeps (with UMFPACK LU
-! on the blocks) as distributed coarsest-level solver
+    ! set a three-level additive Schwarz preconditioner, which uses
+    ! RAS (with overlap 1 and ILU(0) on the blocks) as pre- and
+    ! post-smoother, and 5 block-Jacobi sweeps (with UMFPACK LU
+    ! on the blocks) as distributed coarsest-level solver
 
-  call mld_precinit(P,'ML',info,nlev=3)
-  call mld_precset(P,mld_ml_type_,'ADD',info)
-  call mld_precset(P,mld_smoother_pos_,'TWOSIDE',info)
-  call mld_precset(P,mld_coarse_sweeps_,5,info)
+    call mld_precinit(P,'ML',info,nlev=3)
+    call mld_precset(P,mld_ml_type_,'ADD',info)
+    call mld_precset(P,mld_smoother_pos_,'TWOSIDE',info)
+    call mld_precset(P,mld_coarse_sweeps_,5,info)
 
   end select
 
-! build the preconditioner
+  ! build the preconditioner
 
   call psb_barrier(ictxt)
   t1 = psb_wtime()
@@ -198,13 +198,13 @@ program mld_sexample_ml
     goto 9999
   end if
 
-! set the solver parameters and the initial guess
+  ! set the solver parameters and the initial guess
 
   call psb_geall(x,desc_A,info)
   x(:) =0.0
   call psb_geasb(x,desc_A,info)
 
-! solve Ax=b with preconditioned BiCGSTAB
+  ! solve Ax=b with preconditioned BiCGSTAB
 
   call psb_barrier(ictxt)
   t1 = psb_wtime()
@@ -328,7 +328,7 @@ contains
         integer, intent(out) :: pv(*) 
       end subroutine parts
     end interface
-   ! local variables
+    ! local variables
     type(psb_sspmat_type)    :: a
     real(psb_spk_)           :: zt(nbmax),glob_x,glob_y,glob_z
     integer                  :: m,n,nnz,glob_row,ipoints
@@ -338,7 +338,6 @@ contains
     integer                  :: nv, inv
     integer, allocatable     :: irow(:),icol(:)
     real(psb_spk_), allocatable :: val(:)
-    integer, allocatable     :: prv(:)
     ! deltah dimension of each grid cell
     ! deltat discretization time
     real(psb_spk_)         :: deltah
@@ -368,11 +367,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)
-    ! 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_spall(a,desc_a,info,nnz=nnz)
+    ! define  rhs from boundary conditions; also build initial guess 
+    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
       call psb_errpush(info,name)
       goto 9999
@@ -383,7 +382,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)
@@ -398,143 +397,140 @@ contains
     ! distribution.
 
     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,ipoints*ipoints) == 0) then
-            x = glob_row/(ipoints*ipoints)
-          else
-            x = glob_row/(ipoints*ipoints)+1
-          endif
-          if (mod((glob_row-(x-1)*ipoints*ipoints),ipoints) == 0) then
-            y = (glob_row-(x-1)*ipoints*ipoints)/ipoints
-          else
-            y = (glob_row-(x-1)*ipoints*ipoints)/ipoints+1
-          endif
-          z = glob_row-(x-1)*ipoints*ipoints-(y-1)*ipoints
-          ! 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)*ipoints*ipoints+(y-1)*ipoints+(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_x**2-glob_z**2)*(-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)*ipoints*ipoints+(y-2)*ipoints+(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_x**2-glob_y**2)*(-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)*ipoints*ipoints+(y-1)*ipoints+(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,ipoints*ipoints) == 0) then
+          x = glob_row/(ipoints*ipoints)
+        else
+          x = glob_row/(ipoints*ipoints)+1
+        endif
+        if (mod((glob_row-(x-1)*ipoints*ipoints),ipoints) == 0) then
+          y = (glob_row-(x-1)*ipoints*ipoints)/ipoints
+        else
+          y = (glob_row-(x-1)*ipoints*ipoints)/ipoints+1
+        endif
+        z = glob_row-(x-1)*ipoints*ipoints-(y-1)*ipoints
+        ! 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)*ipoints*ipoints+(y-1)*ipoints+(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_x**2-glob_z**2)*(-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)*ipoints*ipoints+(y-2)*ipoints+(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_x**2-glob_y**2)*(-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)*ipoints*ipoints+(y-1)*ipoints+(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)*ipoints*ipoints+(y-1)*ipoints+(z)
+        element=element+1                  
+        !  term depending on     (x,y,z+1)
+        if (z==ipoints) then 
+          val(element)=-b1(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          zt(1) = exp(-glob_x**2-glob_y**2)*exp(-glob_z)*(-val(element))  
+        else
+          val(element)=-b1(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z+1)
+          element=element+1
+        endif
+        !  term depending on     (x,y+1,z)
+        if (y==ipoints) then 
+          val(element)=-b2(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          zt(1) = exp(-glob_x**2-glob_z**2)*exp(-glob_y)*(-val(element))  
+        else
+          val(element)=-b2(glob_x,glob_y,glob_z)
+          val(element) = val(element)/(deltah*&
+               & deltah)
+          icol(element)=(x-1)*ipoints*ipoints+(y)*ipoints+(z)
+          element=element+1
+        endif
+        !  term depending on     (x+1,y,z)
+        if (x==ipoints) then 
+          val(element)=-b3(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)
           val(element) = val(element)/(deltah*&
                & deltah)
-          icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z)
-          element=element+1                  
-          !  term depending on     (x,y,z+1)
-          if (z==ipoints) then 
-            val(element)=-b1(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            zt(1) = exp(-glob_x**2-glob_y**2)*exp(-glob_z)*(-val(element))  
-          else
-            val(element)=-b1(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z+1)
-            element=element+1
-          endif
-          !  term depending on     (x,y+1,z)
-          if (y==ipoints) then 
-            val(element)=-b2(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-           zt(1) = exp(-glob_x**2-glob_z**2)*exp(-glob_y)*(-val(element))  
-          else
-            val(element)=-b2(glob_x,glob_y,glob_z)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            icol(element)=(x-1)*ipoints*ipoints+(y)*ipoints+(z)
-            element=element+1
-          endif
-          !  term depending on     (x+1,y,z)
-          if (x==ipoints) then 
-            val(element)=-b3(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)
-            val(element) = val(element)/(deltah*&
-                 & deltah)
-            icol(element)=(x)*ipoints*ipoints+(y-1)*ipoints+(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)*ipoints*ipoints+(y-1)*ipoints+(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)    
@@ -564,7 +560,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/ppde.f90 b/tests/pdegen/ppde.f90
index 7831170a..8b553158 100644
--- a/tests/pdegen/ppde.f90
+++ b/tests/pdegen/ppde.f90
@@ -1,4 +1,5 @@
-!!$  
+!!$ 
+!!$ 
 !!$                           MLD2P4  version 1.0
 !!$  MultiLevel Domain Decomposition Parallel Preconditioners Package
 !!$             based on PSBLAS (Parallel Sparse BLAS version 2.2)
@@ -18,14 +19,14 @@
 !!$    2. Redistributions in binary form must reproduce the above copyright
 !!$       notice, this list of conditions, and the following disclaimer in the
 !!$       documentation and/or other materials provided with the distribution.
-!!$    3. The name of the PSBLAS group or the names of its contributors may
+!!$    3. The name of the MLD2P4 group or the names of its contributors may
 !!$       not be used to endorse or promote products derived from this
 !!$       software without specific written permission.
 !!$ 
 !!$  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 !!$  ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
 !!$  TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
-!!$  PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
+!!$  PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
 !!$  BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 !!$  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 !!$  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
@@ -265,7 +266,7 @@ program ppde
 
 contains
   !
-  ! get iteration parameters from the command line
+  ! get iteration parameters from standard input
   !
   subroutine  get_parms(ictxt,kmethd,prectype,afmt,idim,istopc,itmax,itrace,irst)
     integer           :: ictxt
@@ -389,14 +390,15 @@ contains
     ! -   ------ -  ------ -  ------ -  -----  -  ------  -  ------ + a4 u 
     !      dxdx     dydy       dzdz        dx       dy         dz   
     !
-    !  = 0 
-    ! 
-    ! boundary condition: dirichlet
-    !    0< x,y,z<1
-    !  
-    !  u(x,y,z)(2b1+2b2+2b3+a1+a2+a3)+u(x-1,y,z)(-b1-a1)+u(x,y-1,z)(-b2-a2)+
-    !  + u(x,y,z-1)(-b3-a3)-u(x+1,y,z)b1-u(x,y+1,z)b2-u(x,y,z+1)b3
-
+    ! with Dirichlet boundary conditions, on the unit cube  0<=x,y,z<=1.
+    !
+    ! Boundary conditions are set in a very simple way, by adding 
+    ! equations of the form
+    !
+    !   u(x,y) = exp(-x^2-y^2-z^2)
+    !
+    ! Note that if a1=a2=a3=a4=0., the PDE is the well-known Laplace equation.
+    !
     use psb_base_mod
     implicit none
     integer                        :: idim
@@ -422,7 +424,6 @@ contains
     integer                  :: element
     integer, allocatable     :: irow(:),icol(:)
     real(psb_dpk_), allocatable :: val(:)
-    integer, allocatable     :: prv(:)
     ! deltah dimension of each grid cell
     ! deltat discretization time
     real(psb_dpk_)         :: deltah
diff --git a/tests/pdegen/spde.f90 b/tests/pdegen/spde.f90
index 82474467..d3e1cda5 100644
--- a/tests/pdegen/spde.f90
+++ b/tests/pdegen/spde.f90
@@ -1,4 +1,5 @@
-!!$  
+!!$ 
+!!$ 
 !!$                           MLD2P4  version 1.0
 !!$  MultiLevel Domain Decomposition Parallel Preconditioners Package
 !!$             based on PSBLAS (Parallel Sparse BLAS version 2.2)
@@ -18,14 +19,14 @@
 !!$    2. Redistributions in binary form must reproduce the above copyright
 !!$       notice, this list of conditions, and the following disclaimer in the
 !!$       documentation and/or other materials provided with the distribution.
-!!$    3. The name of the PSBLAS group or the names of its contributors may
+!!$    3. The name of the MLD2P4 group or the names of its contributors may
 !!$       not be used to endorse or promote products derived from this
 !!$       software without specific written permission.
 !!$ 
 !!$  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 !!$  ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
 !!$  TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
-!!$  PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
+!!$  PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
 !!$  BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 !!$  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 !!$  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
@@ -265,7 +266,7 @@ program spde
 
 contains
   !
-  ! get iteration parameters from the command line
+  ! get iteration parameters from standard input
   !
   subroutine  get_parms(ictxt,kmethd,prectype,afmt,idim,istopc,itmax,itrace,irst)
     integer           :: ictxt
@@ -389,14 +390,15 @@ contains
     ! -   ------ -  ------ -  ------ -  -----  -  ------  -  ------ + a4 u 
     !      dxdx     dydy       dzdz        dx       dy         dz   
     !
-    !  = 0 
-    ! 
-    ! boundary condition: dirichlet
-    !    0< x,y,z<1
-    !  
-    !  u(x,y,z)(2b1+2b2+2b3+a1+a2+a3)+u(x-1,y,z)(-b1-a1)+u(x,y-1,z)(-b2-a2)+
-    !  + u(x,y,z-1)(-b3-a3)-u(x+1,y,z)b1-u(x,y+1,z)b2-u(x,y,z+1)b3
-
+    ! with Dirichlet boundary conditions, on the unit cube  0<=x,y,z<=1.
+    !
+    ! Boundary conditions are set in a very simple way, by adding 
+    ! equations of the form
+    !
+    !   u(x,y) = exp(-x^2-y^2-z^2)
+    !
+    ! Note that if a1=a2=a3=a4=0., the PDE is the well-known Laplace equation.
+    !
     use psb_base_mod
     implicit none
     integer                        :: idim
@@ -422,7 +424,6 @@ contains
     integer                  :: element
     integer, allocatable     :: irow(:),icol(:)
     real(psb_spk_), allocatable :: val(:)
-    integer, allocatable     :: prv(:)
     ! deltah dimension of each grid cell
     ! deltat discretization time
     real(psb_spk_)         :: deltah