fdurastante
/
amg4psblas
mirror of https://github.com/sfilippone/amg4psblas.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Aligned matrix generation with that from base PSBLAS library.

Browse Source
stopcriterion
Salvatore Filippone 7 years ago
parent 11f25be218
commit 77624b2518
4 changed files with 290 additions and 266 deletions
Show Stats Download Patch File Download Diff File

123
tests/pdegen/mld_d_pde2d.f90
Unescape Escape View File

@ -90,13 +90,64 @@ contains
end function d_null_func_2d
!
! functions parametrizing the differential equation
!
function b1(x,y)
use psb_base_mod, only : psb_dpk_, done, dzero
implicit none
real(psb_dpk_) :: b1
real(psb_dpk_), intent(in) :: x,y
b1=done/sqrt((2*done))
end function b1
function b2(x,y)
use psb_base_mod, only : psb_dpk_, done, dzero
implicit none
real(psb_dpk_) :: b2
real(psb_dpk_), intent(in) :: x,y
b2=done/sqrt((2*done))
end function b2
function c(x,y)
use psb_base_mod, only : psb_dpk_, done, dzero
implicit none
real(psb_dpk_) :: c
real(psb_dpk_), intent(in) :: x,y
c=0.d0
end function c
function a1(x,y)
use psb_base_mod, only : psb_dpk_, done, dzero
implicit none
real(psb_dpk_) :: a1
real(psb_dpk_), intent(in) :: x,y
a1=done/80
end function a1
function a2(x,y)
use psb_base_mod, only : psb_dpk_, done, dzero
implicit none
real(psb_dpk_) :: a2
real(psb_dpk_), intent(in) :: x,y
a2=done/80
end function a2
function g(x,y)
use psb_base_mod, only : psb_dpk_, done, dzero
implicit none
real(psb_dpk_) :: g
real(psb_dpk_), intent(in) :: x,y
g = dzero
if (x == done) then
g = done
else if (x == dzero) then
g = exp(-y**2)
end if
end function g
!
! subroutine to allocate and fill in the coefficient matrix and
! the rhs.
!
subroutine mld_d_gen_pde2d(ictxt,idim,a,bv,xv,desc_a,afmt,&
& a1,a2,b1,b2,c,g,info,f,amold,vmold,imold,partition,nrl,iv)
subroutine mld_d_gen_pde2d(ictxt,idim,a,bv,xv,desc_a,afmt,info,&
& f,amold,vmold,imold,partition,nrl,iv)
use psb_base_mod
use psb_util_mod
!
@ -115,7 +166,6 @@ contains
! Note that if b1=b2=c=0., the PDE is the Laplace equation.
!
implicit none
procedure(d_func_2d) :: b1,b2,c,a1,a2,g
integer(psb_ipk_) :: idim
type(psb_dspmat_type) :: a
type(psb_d_vect_type) :: xv,bv
@ -148,7 +198,7 @@ contains
! deltah dimension of each grid cell
! deltat discretization time
real(psb_dpk_) :: deltah, sqdeltah, deltah2
real(psb_dpk_), parameter :: rhs=0.e0,one=1.e0,zero=0.e0
real(psb_dpk_), parameter :: rhs=dzero,one=done,zero=dzero
real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb
integer(psb_ipk_) :: err_act
procedure(d_func_2d), pointer :: f_
@ -167,9 +217,9 @@ contains
f_ => d_null_func_2d
end if
deltah = 1.e0/(idim+2)
deltah = done/(idim+2)
sqdeltah = deltah*deltah
deltah2 = 2.e0* deltah
deltah2 = (2*done)* deltah
if (present(partition)) then
if ((1<= partition).and.(partition <= 3)) then
@ -355,7 +405,7 @@ contains
if (ix == 1) then
zt(k) = g(dzero,y)*(-val(icoeff)) + zt(k)
else
icol(icoeff) = (ix-2)*idim+iy
call ijk2idx(icol(icoeff),ix-1,iy,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
@ -364,14 +414,14 @@ contains
if (iy == 1) then
zt(k) = g(x,dzero)*(-val(icoeff)) + zt(k)
else
icol(icoeff) = (ix-1)*idim+(iy-1)
call ijk2idx(icol(icoeff),ix,iy-1,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
! term depending on (x,y)
val(icoeff)=2.e0*(a1(x,y) + a2(x,y))/sqdeltah + c(x,y)
icol(icoeff) = (ix-1)*idim+iy
val(icoeff)=(2*done)*(a1(x,y) + a2(x,y))/sqdeltah + c(x,y)
call ijk2idx(icol(icoeff),ix,iy,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
! term depending on (x,y+1)
@ -379,7 +429,7 @@ contains
if (iy == idim) then
zt(k) = g(x,done)*(-val(icoeff)) + zt(k)
else
icol(icoeff) = (ix-1)*idim+(iy+1)
call ijk2idx(icol(icoeff),ix,iy+1,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
@ -388,7 +438,7 @@ contains
if (ix==idim) then
zt(k) = g(done,y)*(-val(icoeff)) + zt(k)
else
icol(icoeff) = (ix)*idim+(iy)
call ijk2idx(icol(icoeff),ix+1,iy,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
@ -398,7 +448,7 @@ contains
if(info /= psb_success_) exit
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info)
if(info /= psb_success_) exit
zt(:)=0.e0
zt(:)=dzero
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
if(info /= psb_success_) exit
end do
@ -469,50 +519,6 @@ contains
end subroutine mld_d_gen_pde2d
!
! functions parametrizing the differential equation
!
function b1(x,y)
use psb_base_mod, only : psb_dpk_,done,dzero
real(psb_dpk_) :: b1
real(psb_dpk_), intent(in) :: x,y
b1=dzero
end function b1
function b2(x,y)
use psb_base_mod, only : psb_dpk_,done,dzero
real(psb_dpk_) :: b2
real(psb_dpk_), intent(in) :: x,y
b2=dzero
end function b2
function c(x,y)
use psb_base_mod, only : psb_dpk_,done,dzero
real(psb_dpk_) :: c
real(psb_dpk_), intent(in) :: x,y
c=dzero
end function c
function a1(x,y)
use psb_base_mod, only : psb_dpk_,done,dzero
real(psb_dpk_) :: a1
real(psb_dpk_), intent(in) :: x,y
a1=done
end function a1
function a2(x,y)
use psb_base_mod, only : psb_dpk_,done,dzero
real(psb_dpk_) :: a2
real(psb_dpk_), intent(in) :: x,y
a2=done
end function a2
function g(x,y)
use psb_base_mod, only : psb_dpk_, done, dzero
real(psb_dpk_) :: g
real(psb_dpk_), intent(in) :: x,y
g = dzero
if (x == done) then
g = done
else if (x == dzero) then
g = exp(-y**2)
end if
end function g
end module mld_d_pde2d_mod
program mld_d_pde2d
@ -654,8 +660,7 @@ program mld_d_pde2d
!
call psb_barrier(ictxt)
t1 = psb_wtime()
call mld_gen_pde2d(ictxt,idim,a,b,x,desc_a,afmt,&
& a1,a2,b1,b2,c,g,info)
call mld_gen_pde2d(ictxt,idim,a,b,x,desc_a,afmt,info)
call psb_barrier(ictxt)
t2 = psb_wtime() - t1
if(info /= psb_success_) then

153
tests/pdegen/mld_d_pde3d.f90
Unescape Escape View File

@ -92,13 +92,78 @@ contains
val = dzero
end function d_null_func_3d
!
! functions parametrizing the differential equation
!
function b1(x,y,z)
use psb_base_mod, only : psb_dpk_, done, dzero
implicit none
real(psb_dpk_) :: b1
real(psb_dpk_), intent(in) :: x,y,z
b1=done/sqrt((3*done))
end function b1
function b2(x,y,z)
use psb_base_mod, only : psb_dpk_, done, dzero
implicit none
real(psb_dpk_) :: b2
real(psb_dpk_), intent(in) :: x,y,z
b2=done/sqrt((3*done))
end function b2
function b3(x,y,z)
use psb_base_mod, only : psb_dpk_, done, dzero
implicit none
real(psb_dpk_) :: b3
real(psb_dpk_), intent(in) :: x,y,z
b3=done/sqrt((3*done))
end function b3
function c(x,y,z)
use psb_base_mod, only : psb_dpk_, done, dzero
implicit none
real(psb_dpk_) :: c
real(psb_dpk_), intent(in) :: x,y,z
c=dzero
end function c
function a1(x,y,z)
use psb_base_mod, only : psb_dpk_, done, dzero
implicit none
real(psb_dpk_) :: a1
real(psb_dpk_), intent(in) :: x,y,z
a1=done/80
end function a1
function a2(x,y,z)
use psb_base_mod, only : psb_dpk_, done, dzero
implicit none
real(psb_dpk_) :: a2
real(psb_dpk_), intent(in) :: x,y,z
a2=done/80
end function a2
function a3(x,y,z)
use psb_base_mod, only : psb_dpk_, done, dzero
implicit none
real(psb_dpk_) :: a3
real(psb_dpk_), intent(in) :: x,y,z
a3=done/80
end function a3
function g(x,y,z)
use psb_base_mod, only : psb_dpk_, done, dzero
implicit none
real(psb_dpk_) :: g
real(psb_dpk_), intent(in) :: x,y,z
g = dzero
if (x == done) then
g = done
else if (x == dzero) then
g = exp(y**2-z**2)
end if
end function g
!
! subroutine to allocate and fill in the coefficient matrix and
! the rhs.
!
subroutine mld_d_gen_pde3d(ictxt,idim,a,bv,xv,desc_a,afmt,&
& a1,a2,a3,b1,b2,b3,c,g,info,f,amold,vmold,imold,partition,nrl,iv)
subroutine mld_d_gen_pde3d(ictxt,idim,a,bv,xv,desc_a,afmt,info,&
& f,amold,vmold,imold,partition,nrl,iv)
use psb_base_mod
use psb_util_mod
!
@ -117,7 +182,6 @@ contains
! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation.
!
implicit none
procedure(d_func_3d) :: b1,b2,b3,c,a1,a2,a3,g
integer(psb_ipk_) :: idim
type(psb_dspmat_type) :: a
type(psb_d_vect_type) :: xv,bv
@ -169,9 +233,9 @@ contains
f_ => d_null_func_3d
end if
deltah = 1.d0/(idim+2)
deltah = done/(idim+2)
sqdeltah = deltah*deltah
deltah2 = 2.d0* deltah
deltah2 = (2*done)* deltah
if (present(partition)) then
if ((1<= partition).and.(partition <= 3)) then
@ -189,7 +253,7 @@ contains
m = idim*idim*idim
n = m
nnz = ((n*9)/(np))
nnz = ((n*7)/(np))
if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n
t0 = psb_wtime()
select case(partition_)
@ -362,7 +426,7 @@ contains
if (ix == 1) then
zt(k) = g(dzero,y,z)*(-val(icoeff)) + zt(k)
else
icol(icoeff) = (ix-2)*idim*idim+(iy-1)*idim+(iz)
call ijk2idx(icol(icoeff),ix-1,iy,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
@ -371,7 +435,7 @@ contains
if (iy == 1) then
zt(k) = g(x,dzero,z)*(-val(icoeff)) + zt(k)
else
icol(icoeff) = (ix-1)*idim*idim+(iy-2)*idim+(iz)
call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
@ -380,15 +444,15 @@ contains
if (iz == 1) then
zt(k) = g(x,y,dzero)*(-val(icoeff)) + zt(k)
else
icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz-1)
call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
! term depending on (x,y,z)
val(icoeff)=2.d0*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah &
val(icoeff)=(2*done)*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah &
& + c(x,y,z)
icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz)
call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
! term depending on (x,y,z+1)
@ -396,7 +460,7 @@ contains
if (iz == idim) then
zt(k) = g(x,y,done)*(-val(icoeff)) + zt(k)
else
icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz+1)
call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
@ -405,7 +469,7 @@ contains
if (iy == idim) then
zt(k) = g(x,done,z)*(-val(icoeff)) + zt(k)
else
icol(icoeff) = (ix-1)*idim*idim+(iy)*idim+(iz)
call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
@ -414,7 +478,7 @@ contains
if (ix==idim) then
zt(k) = g(done,y,z)*(-val(icoeff)) + zt(k)
else
icol(icoeff) = (ix)*idim*idim+(iy-1)*idim+(iz)
call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
@ -424,7 +488,7 @@ contains
if(info /= psb_success_) exit
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info)
if(info /= psb_success_) exit
zt(:)=0.d0
zt(:)=dzero
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
if(info /= psb_success_) exit
end do
@ -494,62 +558,6 @@ contains
return
end subroutine mld_d_gen_pde3d
!
! functions parametrizing the differential equation
!
function b1(x,y,z)
use psb_base_mod, only : psb_dpk_,done,dzero
real(psb_dpk_) :: b1
real(psb_dpk_), intent(in) :: x,y,z
b1=dzero
end function b1
function b2(x,y,z)
use psb_base_mod, only : psb_dpk_,done,dzero
real(psb_dpk_) :: b2
real(psb_dpk_), intent(in) :: x,y,z
b2=dzero
end function b2
function b3(x,y,z)
use psb_base_mod, only : psb_dpk_,done,dzero
real(psb_dpk_) :: b3
real(psb_dpk_), intent(in) :: x,y,z
b3=dzero
end function b3
function c(x,y,z)
use psb_base_mod, only : psb_dpk_,done,dzero
real(psb_dpk_) :: c
real(psb_dpk_), intent(in) :: x,y,z
c=dzero
end function c
function a1(x,y,z)
use psb_base_mod, only : psb_dpk_,done,dzero
real(psb_dpk_) :: a1
real(psb_dpk_), intent(in) :: x,y,z
a1=done
end function a1
function a2(x,y,z)
use psb_base_mod, only : psb_dpk_,done,dzero
real(psb_dpk_) :: a2
real(psb_dpk_), intent(in) :: x,y,z
a2=done
end function a2
function a3(x,y,z)
use psb_base_mod, only : psb_dpk_,done,dzero
real(psb_dpk_) :: a3
real(psb_dpk_), intent(in) :: x,y,z
a3=done
end function a3
function g(x,y,z)
use psb_base_mod, only : psb_dpk_,done,dzero
real(psb_dpk_) :: g
real(psb_dpk_), intent(in) :: x,y,z
g = dzero
if (x == done) then
g = done
else if (x == dzero) then
g = exp(y**2-z**2)
end if
end function g
end module mld_d_pde3d_mod
program mld_d_pde3d
@ -692,8 +700,7 @@ program mld_d_pde3d
call psb_barrier(ictxt)
t1 = psb_wtime()
call mld_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,&
& a1,a2,a3,b1,b2,b3,c,g,info)
call mld_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,info)
call psb_barrier(ictxt)
t2 = psb_wtime() - t1
if(info /= psb_success_) then

123
tests/pdegen/mld_s_pde2d.f90
Unescape Escape View File

@ -90,13 +90,64 @@ contains
end function s_null_func_2d
!
! functions parametrizing the differential equation
!
function b1(x,y)
use psb_base_mod, only : psb_spk_, sone, szero
implicit none
real(psb_spk_) :: b1
real(psb_spk_), intent(in) :: x,y
b1=sone/sqrt((2*sone))
end function b1
function b2(x,y)
use psb_base_mod, only : psb_spk_, sone, szero
implicit none
real(psb_spk_) :: b2
real(psb_spk_), intent(in) :: x,y
b2=sone/sqrt((2*sone))
end function b2
function c(x,y)
use psb_base_mod, only : psb_spk_, sone, szero
implicit none
real(psb_spk_) :: c
real(psb_spk_), intent(in) :: x,y
c=0.d0
end function c
function a1(x,y)
use psb_base_mod, only : psb_spk_, sone, szero
implicit none
real(psb_spk_) :: a1
real(psb_spk_), intent(in) :: x,y
a1=sone/80
end function a1
function a2(x,y)
use psb_base_mod, only : psb_spk_, sone, szero
implicit none
real(psb_spk_) :: a2
real(psb_spk_), intent(in) :: x,y
a2=sone/80
end function a2
function g(x,y)
use psb_base_mod, only : psb_spk_, sone, szero
implicit none
real(psb_spk_) :: g
real(psb_spk_), intent(in) :: x,y
g = szero
if (x == sone) then
g = sone
else if (x == szero) then
g = exp(-y**2)
end if
end function g
!
! subroutine to allocate and fill in the coefficient matrix and
! the rhs.
!
subroutine mld_s_gen_pde2d(ictxt,idim,a,bv,xv,desc_a,afmt,&
& a1,a2,b1,b2,c,g,info,f,amold,vmold,imold,partition,nrl,iv)
subroutine mld_s_gen_pde2d(ictxt,idim,a,bv,xv,desc_a,afmt,info,&
& f,amold,vmold,imold,partition,nrl,iv)
use psb_base_mod
use psb_util_mod
!
@ -115,7 +166,6 @@ contains
! Note that if b1=b2=c=0., the PDE is the Laplace equation.
!
implicit none
procedure(s_func_2d) :: b1,b2,c,a1,a2,g
integer(psb_ipk_) :: idim
type(psb_sspmat_type) :: a
type(psb_s_vect_type) :: xv,bv
@ -148,7 +198,7 @@ contains
! deltah dimension of each grid cell
! deltat discretization time
real(psb_spk_) :: deltah, sqdeltah, deltah2
real(psb_spk_), parameter :: rhs=0.e0,one=1.e0,zero=0.e0
real(psb_spk_), parameter :: rhs=szero,one=sone,zero=szero
real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb
integer(psb_ipk_) :: err_act
procedure(s_func_2d), pointer :: f_
@ -167,9 +217,9 @@ contains
f_ => s_null_func_2d
end if
deltah = 1.e0/(idim+2)
deltah = sone/(idim+2)
sqdeltah = deltah*deltah
deltah2 = 2.e0* deltah
deltah2 = (2*sone)* deltah
if (present(partition)) then
if ((1<= partition).and.(partition <= 3)) then
@ -355,7 +405,7 @@ contains
if (ix == 1) then
zt(k) = g(szero,y)*(-val(icoeff)) + zt(k)
else
icol(icoeff) = (ix-2)*idim+iy
call ijk2idx(icol(icoeff),ix-1,iy,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
@ -364,14 +414,14 @@ contains
if (iy == 1) then
zt(k) = g(x,szero)*(-val(icoeff)) + zt(k)
else
icol(icoeff) = (ix-1)*idim+(iy-1)
call ijk2idx(icol(icoeff),ix,iy-1,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
! term depending on (x,y)
val(icoeff)=2.e0*(a1(x,y) + a2(x,y))/sqdeltah + c(x,y)
icol(icoeff) = (ix-1)*idim+iy
val(icoeff)=(2*sone)*(a1(x,y) + a2(x,y))/sqdeltah + c(x,y)
call ijk2idx(icol(icoeff),ix,iy,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
! term depending on (x,y+1)
@ -379,7 +429,7 @@ contains
if (iy == idim) then
zt(k) = g(x,sone)*(-val(icoeff)) + zt(k)
else
icol(icoeff) = (ix-1)*idim+(iy+1)
call ijk2idx(icol(icoeff),ix,iy+1,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
@ -388,7 +438,7 @@ contains
if (ix==idim) then
zt(k) = g(sone,y)*(-val(icoeff)) + zt(k)
else
icol(icoeff) = (ix)*idim+(iy)
call ijk2idx(icol(icoeff),ix+1,iy,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
@ -398,7 +448,7 @@ contains
if(info /= psb_success_) exit
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info)
if(info /= psb_success_) exit
zt(:)=0.e0
zt(:)=szero
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
if(info /= psb_success_) exit
end do
@ -469,50 +519,6 @@ contains
end subroutine mld_s_gen_pde2d
!
! functions parametrizing the differential equation
!
function b1(x,y)
use psb_base_mod, only : psb_spk_,sone,szero
real(psb_spk_) :: b1
real(psb_spk_), intent(in) :: x,y
b1=szero
end function b1
function b2(x,y)
use psb_base_mod, only : psb_spk_,sone,szero
real(psb_spk_) :: b2
real(psb_spk_), intent(in) :: x,y
b2=szero
end function b2
function c(x,y)
use psb_base_mod, only : psb_spk_,sone,szero
real(psb_spk_) :: c
real(psb_spk_), intent(in) :: x,y
c=szero
end function c
function a1(x,y)
use psb_base_mod, only : psb_spk_,sone,szero
real(psb_spk_) :: a1
real(psb_spk_), intent(in) :: x,y
a1=sone
end function a1
function a2(x,y)
use psb_base_mod, only : psb_spk_,sone,szero
real(psb_spk_) :: a2
real(psb_spk_), intent(in) :: x,y
a2=sone
end function a2
function g(x,y)
use psb_base_mod, only : psb_spk_, sone, szero
real(psb_spk_) :: g
real(psb_spk_), intent(in) :: x,y
g = szero
if (x == sone) then
g = sone
else if (x == szero) then
g = exp(-y**2)
end if
end function g
end module mld_s_pde2d_mod
program mld_s_pde2d
@ -654,8 +660,7 @@ program mld_s_pde2d
!
call psb_barrier(ictxt)
t1 = psb_wtime()
call mld_gen_pde2d(ictxt,idim,a,b,x,desc_a,afmt,&
& a1,a2,b1,b2,c,g,info)
call mld_gen_pde2d(ictxt,idim,a,b,x,desc_a,afmt,info)
call psb_barrier(ictxt)
t2 = psb_wtime() - t1
if(info /= psb_success_) then

153
tests/pdegen/mld_s_pde3d.f90
Unescape Escape View File

@ -92,13 +92,78 @@ contains
val = szero
end function s_null_func_3d
!
! functions parametrizing the differential equation
!
function b1(x,y,z)
use psb_base_mod, only : psb_spk_, sone, szero
implicit none
real(psb_spk_) :: b1
real(psb_spk_), intent(in) :: x,y,z
b1=sone/sqrt((3*sone))
end function b1
function b2(x,y,z)
use psb_base_mod, only : psb_spk_, sone, szero
implicit none
real(psb_spk_) :: b2
real(psb_spk_), intent(in) :: x,y,z
b2=sone/sqrt((3*sone))
end function b2
function b3(x,y,z)
use psb_base_mod, only : psb_spk_, sone, szero
implicit none
real(psb_spk_) :: b3
real(psb_spk_), intent(in) :: x,y,z
b3=sone/sqrt((3*sone))
end function b3
function c(x,y,z)
use psb_base_mod, only : psb_spk_, sone, szero
implicit none
real(psb_spk_) :: c
real(psb_spk_), intent(in) :: x,y,z
c=szero
end function c
function a1(x,y,z)
use psb_base_mod, only : psb_spk_, sone, szero
implicit none
real(psb_spk_) :: a1
real(psb_spk_), intent(in) :: x,y,z
a1=sone/80
end function a1
function a2(x,y,z)
use psb_base_mod, only : psb_spk_, sone, szero
implicit none
real(psb_spk_) :: a2
real(psb_spk_), intent(in) :: x,y,z
a2=sone/80
end function a2
function a3(x,y,z)
use psb_base_mod, only : psb_spk_, sone, szero
implicit none
real(psb_spk_) :: a3
real(psb_spk_), intent(in) :: x,y,z
a3=sone/80
end function a3
function g(x,y,z)
use psb_base_mod, only : psb_spk_, sone, szero
implicit none
real(psb_spk_) :: g
real(psb_spk_), intent(in) :: x,y,z
g = szero
if (x == sone) then
g = sone
else if (x == szero) then
g = exp(y**2-z**2)
end if
end function g
!
! subroutine to allocate and fill in the coefficient matrix and
! the rhs.
!
subroutine mld_s_gen_pde3d(ictxt,idim,a,bv,xv,desc_a,afmt,&
& a1,a2,a3,b1,b2,b3,c,g,info,f,amold,vmold,imold,partition,nrl,iv)
subroutine mld_s_gen_pde3d(ictxt,idim,a,bv,xv,desc_a,afmt,info,&
& f,amold,vmold,imold,partition,nrl,iv)
use psb_base_mod
use psb_util_mod
!
@ -117,7 +182,6 @@ contains
! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation.
!
implicit none
procedure(s_func_3d) :: b1,b2,b3,c,a1,a2,a3,g
integer(psb_ipk_) :: idim
type(psb_sspmat_type) :: a
type(psb_s_vect_type) :: xv,bv
@ -169,9 +233,9 @@ contains
f_ => s_null_func_3d
end if
deltah = 1.d0/(idim+2)
deltah = sone/(idim+2)
sqdeltah = deltah*deltah
deltah2 = 2.d0* deltah
deltah2 = (2*sone)* deltah
if (present(partition)) then
if ((1<= partition).and.(partition <= 3)) then
@ -189,7 +253,7 @@ contains
m = idim*idim*idim
n = m
nnz = ((n*9)/(np))
nnz = ((n*7)/(np))
if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n
t0 = psb_wtime()
select case(partition_)
@ -362,7 +426,7 @@ contains
if (ix == 1) then
zt(k) = g(szero,y,z)*(-val(icoeff)) + zt(k)
else
icol(icoeff) = (ix-2)*idim*idim+(iy-1)*idim+(iz)
call ijk2idx(icol(icoeff),ix-1,iy,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
@ -371,7 +435,7 @@ contains
if (iy == 1) then
zt(k) = g(x,szero,z)*(-val(icoeff)) + zt(k)
else
icol(icoeff) = (ix-1)*idim*idim+(iy-2)*idim+(iz)
call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
@ -380,15 +444,15 @@ contains
if (iz == 1) then
zt(k) = g(x,y,szero)*(-val(icoeff)) + zt(k)
else
icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz-1)
call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
! term depending on (x,y,z)
val(icoeff)=2.d0*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah &
val(icoeff)=(2*sone)*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah &
& + c(x,y,z)
icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz)
call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
! term depending on (x,y,z+1)
@ -396,7 +460,7 @@ contains
if (iz == idim) then
zt(k) = g(x,y,sone)*(-val(icoeff)) + zt(k)
else
icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz+1)
call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
@ -405,7 +469,7 @@ contains
if (iy == idim) then
zt(k) = g(x,sone,z)*(-val(icoeff)) + zt(k)
else
icol(icoeff) = (ix-1)*idim*idim+(iy)*idim+(iz)
call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
@ -414,7 +478,7 @@ contains
if (ix==idim) then
zt(k) = g(sone,y,z)*(-val(icoeff)) + zt(k)
else
icol(icoeff) = (ix)*idim*idim+(iy-1)*idim+(iz)
call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
@ -424,7 +488,7 @@ contains
if(info /= psb_success_) exit
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info)
if(info /= psb_success_) exit
zt(:)=0.d0
zt(:)=szero
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
if(info /= psb_success_) exit
end do
@ -494,62 +558,6 @@ contains
return
end subroutine mld_s_gen_pde3d
!
! functions parametrizing the differential equation
!
function b1(x,y,z)
use psb_base_mod, only : psb_spk_,sone,szero
real(psb_spk_) :: b1
real(psb_spk_), intent(in) :: x,y,z
b1=szero
end function b1
function b2(x,y,z)
use psb_base_mod, only : psb_spk_,sone,szero
real(psb_spk_) :: b2
real(psb_spk_), intent(in) :: x,y,z
b2=szero
end function b2
function b3(x,y,z)
use psb_base_mod, only : psb_spk_,sone,szero
real(psb_spk_) :: b3
real(psb_spk_), intent(in) :: x,y,z
b3=szero
end function b3
function c(x,y,z)
use psb_base_mod, only : psb_spk_,sone,szero
real(psb_spk_) :: c
real(psb_spk_), intent(in) :: x,y,z
c=szero
end function c
function a1(x,y,z)
use psb_base_mod, only : psb_spk_,sone,szero
real(psb_spk_) :: a1
real(psb_spk_), intent(in) :: x,y,z
a1=sone
end function a1
function a2(x,y,z)
use psb_base_mod, only : psb_spk_,sone,szero
real(psb_spk_) :: a2
real(psb_spk_), intent(in) :: x,y,z
a2=sone
end function a2
function a3(x,y,z)
use psb_base_mod, only : psb_spk_,sone,szero
real(psb_spk_) :: a3
real(psb_spk_), intent(in) :: x,y,z
a3=sone
end function a3
function g(x,y,z)
use psb_base_mod, only : psb_spk_,sone,szero
real(psb_spk_) :: g
real(psb_spk_), intent(in) :: x,y,z
g = szero
if (x == sone) then
g = sone
else if (x == szero) then
g = exp(y**2-z**2)
end if
end function g
end module mld_s_pde3d_mod
program mld_s_pde3d
@ -692,8 +700,7 @@ program mld_s_pde3d
call psb_barrier(ictxt)
t1 = psb_wtime()
call mld_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,&
& a1,a2,a3,b1,b2,b3,c,g,info)
call mld_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,info)
call psb_barrier(ictxt)
t2 = psb_wtime() - t1
if(info /= psb_success_) then

Write Preview
Loading…
Cancel
Save