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

mld2p4:

Browse Source 
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.
stopcriterion
Salvatore Filippone 18 years ago
parent e8be375189
commit daa48bbb85
4 changed files with 143 additions and 123 deletions
Show Stats Download Patch File Download Diff File

53
examples/pdegen/mld_dexample_1lev.f90
Unescape Escape View File

@ -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

53
examples/pdegen/mld_dexample_ml.f90
Unescape Escape View File

@ -330,7 +330,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
@ -360,7 +360,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
@ -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

61
examples/pdegen/mld_sexample_1lev.f90
Unescape Escape View File

@ -293,7 +293,7 @@ contains
! 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

61
examples/pdegen/mld_sexample_ml.f90
Unescape Escape View File

@ -331,7 +331,7 @@ contains
! 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
@ -356,12 +356,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
@ -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

Write Preview
Loading…
Cancel
Save