From daa48bbb8516bdbfda7680cb2ba60e5943e0b257 Mon Sep 17 00:00:00 2001
From: Salvatore Filippone <salvatore.filippone@cranfield.ac.uk>
Date: Wed, 23 Jul 2008 09:17:14 +0000
Subject: [PATCH] mld2p4:  examples/pdegen/mld_dexample_1lev.f90 
 examples/pdegen/mld_dexample_ml.f90  examples/pdegen/mld_sexample_1lev.f90 
 examples/pdegen/mld_sexample_ml.f90

Merged changes  to pdegen from Pasqua.
---
 examples/pdegen/mld_dexample_1lev.f90 | 53 +++++++++-------
 examples/pdegen/mld_dexample_ml.f90   | 57 +++++++++--------
 examples/pdegen/mld_sexample_1lev.f90 | 91 ++++++++++++++-------------
 examples/pdegen/mld_sexample_ml.f90   | 65 ++++++++++---------
 4 files changed, 143 insertions(+), 123 deletions(-)

diff --git a/examples/pdegen/mld_dexample_1lev.f90 b/examples/pdegen/mld_dexample_1lev.f90
index 7f43e870..dbc75816 100644
--- a/examples/pdegen/mld_dexample_1lev.f90
+++ b/examples/pdegen/mld_dexample_1lev.f90
@@ -292,7 +292,7 @@ contains
    ! local variables
     type(psb_dspmat_type)    :: a
     real(psb_dpk_)           :: zt(nbmax),glob_x,glob_y,glob_z
-    integer                  :: m,n,nnz,glob_row
+    integer                  :: m,n,nnz,glob_row,ipoints
     integer                  :: x,y,z,ia,indx_owner
     integer                  :: np, iam
     integer                  :: element
@@ -322,7 +322,8 @@ contains
     ! initialize array descriptor and sparse matrix storage; provide an
     ! estimate of the number of non zeroes 
 
-    m   = idim*idim*idim
+    ipoints=idim-2
+    m   = ipoints*ipoints*ipoints
     n   = m
     nnz = ((n*9)/(np))
     if(iam == psb_root_) write(0,'("Generating Matrix (size=",i0x,")...")')n
@@ -357,7 +358,6 @@ contains
     ! loop over rows belonging to current process in a block
     ! distribution.
 
-    !    icol(1)=1    
     do glob_row = 1, n
       call parts(glob_row,n,np,prv,nv)
       do inv = 1, nv
@@ -366,17 +366,17 @@ contains
           ! local matrix pointer 
           element=1
           ! compute gridpoint coordinates
-          if (mod(glob_row,(idim*idim)) == 0) then
-            x = glob_row/(idim*idim)
+          if (mod(glob_row,ipoints*ipoints) == 0) then
+            x = glob_row/(ipoints*ipoints)
           else
-            x = glob_row/(idim*idim)+1
+            x = glob_row/(ipoints*ipoints)+1
           endif
-          if (mod((glob_row-(x-1)*idim*idim),idim) == 0) then
-            y = (glob_row-(x-1)*idim*idim)/idim
+          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)*idim*idim)/idim+1
+            y = (glob_row-(x-1)*ipoints*ipoints)/ipoints+1
           endif
-          z = glob_row-(x-1)*idim*idim-(y-1)*idim
+          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
@@ -399,7 +399,7 @@ contains
                  & -a1(glob_x,glob_y,glob_z)
             val(element) = val(element)/(deltah*&
                  & deltah)
-            icol(element)=(x-2)*idim*idim+(y-1)*idim+(z)
+            icol(element)=(x-2)*ipoints*ipoints+(y-1)*ipoints+(z)
             element=element+1
           endif
           !  term depending on     (x,y-1,z)
@@ -408,13 +408,13 @@ contains
                  & -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))  
+            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)*idim*idim+(y-2)*idim+(z)
+            icol(element)=(x-1)*ipoints*ipoints+(y-2)*ipoints+(z)
             element=element+1
           endif
           !  term depending on     (x,y,z-1)
@@ -423,13 +423,13 @@ contains
                  & -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))  
+            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)*idim*idim+(y-1)*idim+(z-1)
+            icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z-1)
             element=element+1
           endif
           !  term depending on     (x,y,z)
@@ -441,40 +441,45 @@ contains
                & +a3(glob_x,glob_y,glob_z)
           val(element) = val(element)/(deltah*&
                & deltah)
-          icol(element)=(x-1)*idim*idim+(y-1)*idim+(z)
+          icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z)
           element=element+1                  
           !  term depending on     (x,y,z+1)
-          if (z==idim) then 
+          if (z==ipoints) 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))  
+            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)*idim*idim+(y-1)*idim+(z+1)
+            icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z+1)
             element=element+1
           endif
           !  term depending on     (x,y+1,z)
-          if (y==idim) then 
+          if (y==ipoints) 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))  
+           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)*idim*idim+(y)*idim+(z)
+            icol(element)=(x-1)*ipoints*ipoints+(y)*ipoints+(z)
             element=element+1
           endif
           !  term depending on     (x+1,y,z)
-          if (x<idim) then 
+          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)*idim*idim+(y-1)*idim+(z)
+            icol(element)=(x)*ipoints*ipoints+(y-1)*ipoints+(z)
             element=element+1
           endif
           irow(1:element-1)=glob_row
diff --git a/examples/pdegen/mld_dexample_ml.f90 b/examples/pdegen/mld_dexample_ml.f90
index d2607fff..4cd836e5 100644
--- a/examples/pdegen/mld_dexample_ml.f90
+++ b/examples/pdegen/mld_dexample_ml.f90
@@ -329,8 +329,8 @@ contains
     end interface
    ! local variables
     type(psb_dspmat_type)    :: a
-    real(psb_dpk_)         :: zt(nbmax),glob_x,glob_y,glob_z
-    integer                  :: m,n,nnz,glob_row
+    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
@@ -360,14 +360,15 @@ contains
     ! initialize array descriptor and sparse matrix storage; provide an
     ! estimate of the number of non zeroes 
 
-    m   = idim*idim*idim
+    ipoints=idim-2
+    m   = ipoints*ipoints*ipoints
     n   = m
     nnz = ((n*9)/(np))
     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 
+    ! 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
@@ -395,7 +396,6 @@ contains
     ! loop over rows belonging to current process in a block
     ! distribution.
 
-    !    icol(1)=1    
     do glob_row = 1, n
       call parts(glob_row,n,np,prv,nv)
       do inv = 1, nv
@@ -404,17 +404,17 @@ contains
           ! local matrix pointer 
           element=1
           ! compute gridpoint coordinates
-          if (mod(glob_row,(idim*idim)) == 0) then
-            x = glob_row/(idim*idim)
+          if (mod(glob_row,ipoints*ipoints) == 0) then
+            x = glob_row/(ipoints*ipoints)
           else
-            x = glob_row/(idim*idim)+1
+            x = glob_row/(ipoints*ipoints)+1
           endif
-          if (mod((glob_row-(x-1)*idim*idim),idim) == 0) then
-            y = (glob_row-(x-1)*idim*idim)/idim
+          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)*idim*idim)/idim+1
+            y = (glob_row-(x-1)*ipoints*ipoints)/ipoints+1
           endif
-          z = glob_row-(x-1)*idim*idim-(y-1)*idim
+          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
@@ -437,7 +437,7 @@ contains
                  & -a1(glob_x,glob_y,glob_z)
             val(element) = val(element)/(deltah*&
                  & deltah)
-            icol(element)=(x-2)*idim*idim+(y-1)*idim+(z)
+            icol(element)=(x-2)*ipoints*ipoints+(y-1)*ipoints+(z)
             element=element+1
           endif
           !  term depending on     (x,y-1,z)
@@ -446,13 +446,13 @@ contains
                  & -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))  
+            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)*idim*idim+(y-2)*idim+(z)
+            icol(element)=(x-1)*ipoints*ipoints+(y-2)*ipoints+(z)
             element=element+1
           endif
           !  term depending on     (x,y,z-1)
@@ -461,13 +461,13 @@ contains
                  & -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))  
+            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)*idim*idim+(y-1)*idim+(z-1)
+            icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z-1)
             element=element+1
           endif
           !  term depending on     (x,y,z)
@@ -479,40 +479,45 @@ contains
                & +a3(glob_x,glob_y,glob_z)
           val(element) = val(element)/(deltah*&
                & deltah)
-          icol(element)=(x-1)*idim*idim+(y-1)*idim+(z)
+          icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z)
           element=element+1                  
           !  term depending on     (x,y,z+1)
-          if (z==idim) then 
+          if (z==ipoints) 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))  
+            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)*idim*idim+(y-1)*idim+(z+1)
+            icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z+1)
             element=element+1
           endif
           !  term depending on     (x,y+1,z)
-          if (y==idim) then 
+          if (y==ipoints) 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))  
+           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)*idim*idim+(y)*idim+(z)
+            icol(element)=(x-1)*ipoints*ipoints+(y)*ipoints+(z)
             element=element+1
           endif
           !  term depending on     (x+1,y,z)
-          if (x<idim) then 
+          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)*idim*idim+(y-1)*idim+(z)
+            icol(element)=(x)*ipoints*ipoints+(y-1)*ipoints+(z)
             element=element+1
           endif
           irow(1:element-1)=glob_row
diff --git a/examples/pdegen/mld_sexample_1lev.f90 b/examples/pdegen/mld_sexample_1lev.f90
index 5f885a57..f3a1241a 100644
--- a/examples/pdegen/mld_sexample_1lev.f90
+++ b/examples/pdegen/mld_sexample_1lev.f90
@@ -81,33 +81,33 @@ program mld_sexample_1lev
   implicit none
 
 
-! sparse matrices
+  ! sparse matrices
   type(psb_sspmat_type) :: A
 
-! descriptor of sparse matrices
+  ! descriptor of sparse matrices
   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 parameters
+  ! solver parameters
   real(psb_spk_)   :: 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
   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)
 
@@ -122,11 +122,11 @@ program mld_sexample_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()
@@ -142,12 +142,12 @@ program mld_sexample_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()
@@ -162,13 +162,13 @@ program mld_sexample_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()
@@ -290,10 +290,10 @@ 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
+    integer                  :: m,n,nnz,glob_row,ipoints
     integer                  :: x,y,z,ia,indx_owner
     integer                  :: np, iam
     integer                  :: element
@@ -318,12 +318,13 @@ contains
 
     call psb_info(ictxt, iam, np)
 
-    deltah = 1.e0/(idim-1)
+    deltah = 1.d0/(idim-1)
 
     ! initialize array descriptor and sparse matrix storage; provide an
     ! estimate of the number of non zeroes 
 
-    m   = idim*idim*idim
+    ipoints=idim-2
+    m   = ipoints*ipoints*ipoints
     n   = m
     nnz = ((n*9)/(np))
     if(iam == psb_root_) write(0,'("Generating Matrix (size=",i0x,")...")')n
@@ -351,14 +352,13 @@ contains
       goto 9999
     endif
 
-    tins = 0.e0
+    tins = 0.d0
     call psb_barrier(ictxt)
     t1 = psb_wtime()
 
     ! loop over rows belonging to current process in a block
     ! distribution.
 
-    !    icol(1)=1    
     do glob_row = 1, n
       call parts(glob_row,n,np,prv,nv)
       do inv = 1, nv
@@ -367,24 +367,24 @@ contains
           ! local matrix pointer 
           element=1
           ! compute gridpoint coordinates
-          if (mod(glob_row,(idim*idim)) == 0) then
-            x = glob_row/(idim*idim)
+          if (mod(glob_row,ipoints*ipoints) == 0) then
+            x = glob_row/(ipoints*ipoints)
           else
-            x = glob_row/(idim*idim)+1
+            x = glob_row/(ipoints*ipoints)+1
           endif
-          if (mod((glob_row-(x-1)*idim*idim),idim) == 0) then
-            y = (glob_row-(x-1)*idim*idim)/idim
+          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)*idim*idim)/idim+1
+            y = (glob_row-(x-1)*ipoints*ipoints)/ipoints+1
           endif
-          z = glob_row-(x-1)*idim*idim-(y-1)*idim
+          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.e0
+          zt(1) = 0.d0
           ! internal point: build discretization
           !   
           !  term depending on   (x-1,y,z)
@@ -400,7 +400,7 @@ contains
                  & -a1(glob_x,glob_y,glob_z)
             val(element) = val(element)/(deltah*&
                  & deltah)
-            icol(element)=(x-2)*idim*idim+(y-1)*idim+(z)
+            icol(element)=(x-2)*ipoints*ipoints+(y-1)*ipoints+(z)
             element=element+1
           endif
           !  term depending on     (x,y-1,z)
@@ -409,13 +409,13 @@ contains
                  & -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))  
+            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)*idim*idim+(y-2)*idim+(z)
+            icol(element)=(x-1)*ipoints*ipoints+(y-2)*ipoints+(z)
             element=element+1
           endif
           !  term depending on     (x,y,z-1)
@@ -424,13 +424,13 @@ contains
                  & -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))  
+            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)*idim*idim+(y-1)*idim+(z-1)
+            icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z-1)
             element=element+1
           endif
           !  term depending on     (x,y,z)
@@ -442,40 +442,45 @@ contains
                & +a3(glob_x,glob_y,glob_z)
           val(element) = val(element)/(deltah*&
                & deltah)
-          icol(element)=(x-1)*idim*idim+(y-1)*idim+(z)
+          icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z)
           element=element+1                  
           !  term depending on     (x,y,z+1)
-          if (z==idim) then 
+          if (z==ipoints) 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))  
+            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)*idim*idim+(y-1)*idim+(z+1)
+            icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z+1)
             element=element+1
           endif
           !  term depending on     (x,y+1,z)
-          if (y==idim) then 
+          if (y==ipoints) 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))  
+            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)*idim*idim+(y)*idim+(z)
+            icol(element)=(x-1)*ipoints*ipoints+(y)*ipoints+(z)
             element=element+1
           endif
           !  term depending on     (x+1,y,z)
-          if (x<idim) then 
+          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)*idim*idim+(y-1)*idim+(z)
+            icol(element)=(x)*ipoints*ipoints+(y-1)*ipoints+(z)
             element=element+1
           endif
           irow(1:element-1)=glob_row
@@ -487,7 +492,7 @@ contains
           tins = tins + (psb_wtime()-t3)
           call psb_geins(1,(/ia/),zt(1:1),b,desc_a,info)
           if(info /= 0) exit
-          zt(1)=0.e0
+          zt(1)=0.d0
           call psb_geins(1,(/ia/),zt(1:1),xv,desc_a,info)
           if(info /= 0) exit
         end if
diff --git a/examples/pdegen/mld_sexample_ml.f90 b/examples/pdegen/mld_sexample_ml.f90
index e97f5d46..6f91dab0 100644
--- a/examples/pdegen/mld_sexample_ml.f90
+++ b/examples/pdegen/mld_sexample_ml.f90
@@ -330,8 +330,8 @@ contains
     end interface
    ! local variables
     type(psb_sspmat_type)    :: a
-    real(psb_spk_)         :: zt(nbmax),glob_x,glob_y,glob_z
-    integer                  :: m,n,nnz,glob_row
+    real(psb_spk_)           :: 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
@@ -356,19 +356,20 @@ contains
 
     call psb_info(ictxt, iam, np)
 
-    deltah = 1.e0/(idim-1)
+    deltah = 1.d0/(idim-1)
 
     ! initialize array descriptor and sparse matrix storage; provide an
     ! estimate of the number of non zeroes 
 
-    m   = idim*idim*idim
+    ipoints=idim-2
+    m   = ipoints*ipoints*ipoints
     n   = m
     nnz = ((n*9)/(np))
     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 
+    ! 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
@@ -389,14 +390,13 @@ contains
       goto 9999
     endif
 
-    tins = 0.e0
+    tins = 0.d0
     call psb_barrier(ictxt)
     t1 = psb_wtime()
 
     ! loop over rows belonging to current process in a block
     ! distribution.
 
-    !    icol(1)=1    
     do glob_row = 1, n
       call parts(glob_row,n,np,prv,nv)
       do inv = 1, nv
@@ -405,24 +405,24 @@ contains
           ! local matrix pointer 
           element=1
           ! compute gridpoint coordinates
-          if (mod(glob_row,(idim*idim)) == 0) then
-            x = glob_row/(idim*idim)
+          if (mod(glob_row,ipoints*ipoints) == 0) then
+            x = glob_row/(ipoints*ipoints)
           else
-            x = glob_row/(idim*idim)+1
+            x = glob_row/(ipoints*ipoints)+1
           endif
-          if (mod((glob_row-(x-1)*idim*idim),idim) == 0) then
-            y = (glob_row-(x-1)*idim*idim)/idim
+          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)*idim*idim)/idim+1
+            y = (glob_row-(x-1)*ipoints*ipoints)/ipoints+1
           endif
-          z = glob_row-(x-1)*idim*idim-(y-1)*idim
+          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.e0
+          zt(1) = 0.d0
           ! internal point: build discretization
           !   
           !  term depending on   (x-1,y,z)
@@ -438,7 +438,7 @@ contains
                  & -a1(glob_x,glob_y,glob_z)
             val(element) = val(element)/(deltah*&
                  & deltah)
-            icol(element)=(x-2)*idim*idim+(y-1)*idim+(z)
+            icol(element)=(x-2)*ipoints*ipoints+(y-1)*ipoints+(z)
             element=element+1
           endif
           !  term depending on     (x,y-1,z)
@@ -447,13 +447,13 @@ contains
                  & -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))  
+            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)*idim*idim+(y-2)*idim+(z)
+            icol(element)=(x-1)*ipoints*ipoints+(y-2)*ipoints+(z)
             element=element+1
           endif
           !  term depending on     (x,y,z-1)
@@ -462,13 +462,13 @@ contains
                  & -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))  
+            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)*idim*idim+(y-1)*idim+(z-1)
+            icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z-1)
             element=element+1
           endif
           !  term depending on     (x,y,z)
@@ -480,40 +480,45 @@ contains
                & +a3(glob_x,glob_y,glob_z)
           val(element) = val(element)/(deltah*&
                & deltah)
-          icol(element)=(x-1)*idim*idim+(y-1)*idim+(z)
+          icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z)
           element=element+1                  
           !  term depending on     (x,y,z+1)
-          if (z==idim) then 
+          if (z==ipoints) 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))  
+            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)*idim*idim+(y-1)*idim+(z+1)
+            icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z+1)
             element=element+1
           endif
           !  term depending on     (x,y+1,z)
-          if (y==idim) then 
+          if (y==ipoints) 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))  
+           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)*idim*idim+(y)*idim+(z)
+            icol(element)=(x-1)*ipoints*ipoints+(y)*ipoints+(z)
             element=element+1
           endif
           !  term depending on     (x+1,y,z)
-          if (x<idim) then 
+          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)*idim*idim+(y-1)*idim+(z)
+            icol(element)=(x)*ipoints*ipoints+(y-1)*ipoints+(z)
             element=element+1
           endif
           irow(1:element-1)=glob_row
@@ -525,7 +530,7 @@ contains
           tins = tins + (psb_wtime()-t3)
           call psb_geins(1,(/ia/),zt(1:1),b,desc_a,info)
           if(info /= 0) exit
-          zt(1)=0.e0
+          zt(1)=0.d0
           call psb_geins(1,(/ia/),zt(1:1),xv,desc_a,info)
           if(info /= 0) exit
         end if