psblas3:

test/pargen/ppde.f90

Changed matrix generation routine in sample program. Must be done in a
better way!

test/pargen/ppde.f90

@@ -39,7 +39,7 @@
 ! The PDE is a general second order equation in 3d
 !
 !   b1 dd(u)  b2 dd(u)   b3 dd(u)  a1 d(u)    a2 d(u)    a3 d(u)  
-! -   ------ -  ------ -  ------ -  -----  -  ------  -  ------ + a4 u  = 0
+! -   ------ -  ------ -  ------ +  -----  +  ------  +  ------ + a4 u  = 0
 !      dxdx     dydy       dzdz       dx       dy         dz   
 !
 ! with Dirichlet boundary conditions, on the unit cube  0<=x,y,z<=1.
@@ -148,7 +148,7 @@ program ppde
 !!$  call psb_cdcpy(desc_a,desc_b,info)
 !!$  call psb_set_debug_level(9999)
 
-  call psb_cdbldext(a,desc_a,2,desc_b,info,extype=psb_ovt_asov_)
+!!$  call psb_cdbldext(a,desc_a,2,desc_b,info,extype=psb_ovt_asov_)
   if (info /= 0) then 
     write(0,*) 'Error from bldext'
     call psb_abort(ictxt)
@@ -365,7 +365,7 @@ contains
     !   discretize the partial diferential equation
     ! 
     !   b1 dd(u)  b2 dd(u)    b3 dd(u)    a1 d(u)   a2 d(u)  a3 d(u)  
-    ! -   ------ -  ------ -  ------ -  -----  -  ------  -  ------ + a4 u 
+    ! -   ------ -  ------ -  ------ +  -----  +  ------  +  ------ + a4 u 
     !      dxdx     dydy       dzdz        dx       dy         dz   
     !
     ! with Dirichlet boundary conditions, on the unit cube  0<=x,y,z<=1.
@@ -399,7 +399,7 @@ contains
     real(psb_dpk_), allocatable :: val(:)
     ! deltah dimension of each grid cell
     ! deltat discretization time
-    real(psb_dpk_)            :: deltah, deltah2
+    real(psb_dpk_)            :: deltah, sqdeltah, deltah2
     real(psb_dpk_), parameter :: rhs=0.d0,one=1.d0,zero=0.d0
     real(psb_dpk_)   :: t0, t1, t2, t3, tasb, talc, ttot, tgen 
     real(psb_dpk_)   :: a1, a2, a3, a4, b1, b2, b3 
@@ -415,7 +415,8 @@ contains
     call psb_info(ictxt, iam, np)
 
     deltah   = 1.d0/(idim-1)
-    deltah2 = deltah*deltah
+    sqdeltah = deltah*deltah
+    deltah2  = 2.d0* deltah
 
     ! initialize array descriptor and sparse matrix storage. provide an
     ! estimate of the number of non zeroes 
@@ -507,71 +508,72 @@ contains
         !  term depending on   (x-1,y,z)
         !
         if (ix == 1) then 
-          val(element) = -b1(x,y,z)/deltah2-a1(x,y,z)/deltah
+          val(element) = -b1(x,y,z)/sqdeltah-a1(x,y,z)/deltah2
           zt(k) = exp(-x**2-y**2-z**2)*(-val(element))
         else
-          val(element)  = -b1(x,y,z)/deltah2-a1(x,y,z)/deltah
+          val(element)  = -b1(x,y,z)/sqdeltah-a1(x,y,z)/deltah2
           icol(element) = (ix-2)*idim*idim+(iy-1)*idim+(iz)
           irow(element) = glob_row
           element       = element+1
         endif
         !  term depending on     (x,y-1,z)
         if (iy == 1) then 
-          val(element)  = -b2(x,y,z)/deltah2-a2(x,y,z)/deltah
+          val(element)  = -b2(x,y,z)/sqdeltah-a2(x,y,z)/deltah2
           zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element))  
         else
-          val(element)  = -b2(x,y,z)/deltah2-a2(x,y,z)/deltah
+          val(element)  = -b2(x,y,z)/sqdeltah-a2(x,y,z)/deltah2
           icol(element) = (ix-1)*idim*idim+(iy-2)*idim+(iz)
           irow(element) = glob_row
           element       = element+1
         endif
         !  term depending on     (x,y,z-1)
         if (iz == 1) then 
-          val(element)=-b3(x,y,z)/deltah2-a3(x,y,z)/deltah
+          val(element)=-b3(x,y,z)/sqdeltah-a3(x,y,z)/deltah2
           zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element))  
         else
-          val(element)=-b3(x,y,z)/deltah2-a3(x,y,z)/deltah
+          val(element)=-b3(x,y,z)/sqdeltah-a3(x,y,z)/deltah2
           icol(element) = (ix-1)*idim*idim+(iy-1)*idim+(iz-1)
           irow(element) = glob_row
           element       = element+1
         endif
         !  term depending on     (x,y,z)
-        val(element)=(2*b1(x,y,z) + 2*b2(x,y,z) + 2*b3(x,y,z))/deltah2&
+        val(element)=(2*b1(x,y,z) + 2*b2(x,y,z) + 2*b3(x,y,z))/sqdeltah&
-             & + (a1(x,y,z) + a2(x,y,z) + a3(x,y,z)+ a4(x,y,z))/deltah
+             & +a4(x,y,z)
         icol(element) = (ix-1)*idim*idim+(iy-1)*idim+(iz)
         irow(element) = glob_row
         element       = element+1                  
         !  term depending on     (x,y,z+1)
         if (iz == idim) then 
-          val(element)=-b1(x,y,z)/deltah2
+          val(element)=-b3(x,y,z)/sqdeltah+a3(x,y,z)/deltah2
           zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element))  
         else
-          val(element)=-b1(x,y,z)/deltah2
+          val(element)=-b3(x,y,z)/sqdeltah+a3(x,y,z)/deltah2
           icol(element) = (ix-1)*idim*idim+(iy-1)*idim+(iz+1)
           irow(element) = glob_row
           element       = element+1
         endif
         !  term depending on     (x,y+1,z)
         if (iy == idim) then 
-          val(element)=-b2(x,y,z)/deltah2
+          val(element)=-b2(x,y,z)/sqdeltah+a2(x,y,z)/deltah2
           zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element))  
         else
-          val(element)=-b2(x,y,z)/deltah2
+          val(element)=-b2(x,y,z)/sqdeltah+a2(x,y,z)/deltah2
           icol(element) = (ix-1)*idim*idim+(iy)*idim+(iz)
           irow(element) = glob_row
           element       = element+1
         endif
         !  term depending on     (x+1,y,z)
         if (ix==idim) then 
-          val(element)=-b3(x,y,z)/deltah2
+          val(element)=-b1(x,y,z)/sqdeltah+a1(x,y,z)/deltah2
           zt(k) = exp(-y**2-z**2)*exp(-x)*(-val(element))  
         else
-          val(element)=-b3(x,y,z)/deltah2
+          val(element)=-b1(x,y,z)/sqdeltah+a1(x,y,z)/deltah2
           icol(element) = (ix)*idim*idim+(iy-1)*idim+(iz)
           irow(element) = glob_row
           element       = element+1
         endif
 
+
       end do
       call psb_spins(element-1,irow,icol,val,a,desc_a,info)
       if(info /= psb_success_) exit
@@ -649,43 +651,43 @@ function a1(x,y,z)
   use psb_base_mod, only : psb_dpk_
   real(psb_dpk_) :: a1
   real(psb_dpk_) :: x,y,z
-  a1=1.d0
+  a1=1.414d0
 end function a1
 function a2(x,y,z)
   use psb_base_mod, only : psb_dpk_
   real(psb_dpk_) ::  a2
   real(psb_dpk_) :: x,y,z
-  a2=2.d1*y
+  a2=1.414d0
 end function a2
 function a3(x,y,z)
   use psb_base_mod, only : psb_dpk_
   real(psb_dpk_) ::  a3
   real(psb_dpk_) :: x,y,z      
-  a3=1.d0
+  a3=1.414d0
 end function a3
 function a4(x,y,z)
   use psb_base_mod, only : psb_dpk_
   real(psb_dpk_) ::  a4
   real(psb_dpk_) :: x,y,z      
-  a4=1.d0
+  a4=0.d0
 end function a4
 function b1(x,y,z)
   use psb_base_mod, only : psb_dpk_
   real(psb_dpk_) ::  b1   
   real(psb_dpk_) :: x,y,z
-  b1=1.d0
+  b1=1.d0/80
 end function b1
 function b2(x,y,z)
   use psb_base_mod, only : psb_dpk_
   real(psb_dpk_) ::  b2
   real(psb_dpk_) :: x,y,z
-  b2=1.d0
+  b2=1.d0/80
 end function b2
 function b3(x,y,z)
   use psb_base_mod, only : psb_dpk_
   real(psb_dpk_) ::  b3
   real(psb_dpk_) :: x,y,z
-  b3=1.d0
+  b3=1.d0/80
 end function b3