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

mld2p4:

Browse Source 
tests/fileread/cf_sample.f90
 tests/fileread/df_sample.f90
 tests/fileread/sf_sample.f90
 tests/fileread/zf_sample.f90
 tests/pdegen/ppde.f90
 tests/pdegen/spde.f90

Sample programs: formatting, comments, ownership function in matrix generation.
stopcriterion
Salvatore Filippone 18 years ago
parent 0a436e66d8
commit 8782eb20df
6 changed files with 272 additions and 282 deletions
Show Stats Download Patch File Download Diff File

2
tests/fileread/cf_sample.f90
Unescape Escape View File

@ -264,7 +264,7 @@ program cf_sample
call mld_precset(prec,mld_coarse_subsolve_, prec_choice%csbsolve,info) call mld_precset(prec,mld_coarse_subsolve_, prec_choice%csbsolve,info)
call mld_precset(prec,mld_coarse_mat_, prec_choice%cmat, info) call mld_precset(prec,mld_coarse_mat_, prec_choice%cmat, info)
call mld_precset(prec,mld_coarse_fillin_, prec_choice%cfill, info) call mld_precset(prec,mld_coarse_fillin_, prec_choice%cfill, info)
call mld_precset(prec,mld_coarse_iluthrs_, prec_choice%cthres, info) call mld_precset(prec,mld_coarse_iluthrs_, prec_choice%cthres, info)
call mld_precset(prec,mld_coarse_sweeps_, prec_choice%cjswp, info) call mld_precset(prec,mld_coarse_sweeps_, prec_choice%cjswp, info)
end if end if

2
tests/fileread/df_sample.f90
Unescape Escape View File

@ -264,7 +264,7 @@ program df_sample
call mld_precset(prec,mld_coarse_subsolve_, prec_choice%csbsolve,info) call mld_precset(prec,mld_coarse_subsolve_, prec_choice%csbsolve,info)
call mld_precset(prec,mld_coarse_mat_, prec_choice%cmat, info) call mld_precset(prec,mld_coarse_mat_, prec_choice%cmat, info)
call mld_precset(prec,mld_coarse_fillin_, prec_choice%cfill, info) call mld_precset(prec,mld_coarse_fillin_, prec_choice%cfill, info)
call mld_precset(prec,mld_coarse_iluthrs_, prec_choice%cthres, info) call mld_precset(prec,mld_coarse_iluthrs_, prec_choice%cthres, info)
call mld_precset(prec,mld_coarse_sweeps_, prec_choice%cjswp, info) call mld_precset(prec,mld_coarse_sweeps_, prec_choice%cjswp, info)
end if end if

2
tests/fileread/sf_sample.f90
Unescape Escape View File

@ -264,7 +264,7 @@ program sf_sample
call mld_precset(prec,mld_coarse_subsolve_, prec_choice%csbsolve,info) call mld_precset(prec,mld_coarse_subsolve_, prec_choice%csbsolve,info)
call mld_precset(prec,mld_coarse_mat_, prec_choice%cmat, info) call mld_precset(prec,mld_coarse_mat_, prec_choice%cmat, info)
call mld_precset(prec,mld_coarse_fillin_, prec_choice%cfill, info) call mld_precset(prec,mld_coarse_fillin_, prec_choice%cfill, info)
call mld_precset(prec,mld_coarse_iluthrs_, prec_choice%cthres, info) call mld_precset(prec,mld_coarse_iluthrs_, prec_choice%cthres, info)
call mld_precset(prec,mld_coarse_sweeps_, prec_choice%cjswp, info) call mld_precset(prec,mld_coarse_sweeps_, prec_choice%cjswp, info)
end if end if

2
tests/fileread/zf_sample.f90
Unescape Escape View File

@ -264,7 +264,7 @@ program zf_sample
call mld_precset(prec,mld_coarse_subsolve_, prec_choice%csbsolve,info) call mld_precset(prec,mld_coarse_subsolve_, prec_choice%csbsolve,info)
call mld_precset(prec,mld_coarse_mat_, prec_choice%cmat, info) call mld_precset(prec,mld_coarse_mat_, prec_choice%cmat, info)
call mld_precset(prec,mld_coarse_fillin_, prec_choice%cfill, info) call mld_precset(prec,mld_coarse_fillin_, prec_choice%cfill, info)
call mld_precset(prec,mld_coarse_iluthrs_, prec_choice%cthres, info) call mld_precset(prec,mld_coarse_iluthrs_, prec_choice%cthres, info)
call mld_precset(prec,mld_coarse_sweeps_, prec_choice%cjswp, info) call mld_precset(prec,mld_coarse_sweeps_, prec_choice%cjswp, info)
end if end if

273
tests/pdegen/ppde.f90
Unescape Escape View File

@ -189,7 +189,7 @@ program ppde
call mld_precset(prec,mld_coarse_subsolve_, prectype%csbsolve,info) call mld_precset(prec,mld_coarse_subsolve_, prectype%csbsolve,info)
call mld_precset(prec,mld_coarse_mat_, prectype%cmat, info) call mld_precset(prec,mld_coarse_mat_, prectype%cmat, info)
call mld_precset(prec,mld_coarse_fillin_, prectype%cfill, info) call mld_precset(prec,mld_coarse_fillin_, prectype%cfill, info)
call mld_precset(prec,mld_coarse_iluthrs_, prectype%cthres, info) call mld_precset(prec,mld_coarse_iluthrs_, prectype%cthres, info)
call mld_precset(prec,mld_coarse_sweeps_, prectype%cjswp, info) call mld_precset(prec,mld_coarse_sweeps_, prectype%cjswp, info)
end if end if
@ -356,7 +356,6 @@ contains
return return
end subroutine get_parms end subroutine get_parms
! !
! print an error message ! print an error message
! !
@ -421,7 +420,6 @@ contains
integer :: x,y,z,ia,indx_owner integer :: x,y,z,ia,indx_owner
integer :: np, iam integer :: np, iam
integer :: element integer :: element
integer :: nv, inv
integer, allocatable :: irow(:),icol(:) integer, allocatable :: irow(:),icol(:)
real(psb_dpk_), allocatable :: val(:) real(psb_dpk_), allocatable :: val(:)
integer, allocatable :: prv(:) integer, allocatable :: prv(:)
@ -454,11 +452,11 @@ contains
if(iam == psb_root_) write(0,'("Generating Matrix (size=",i0x,")...")')n if(iam == psb_root_) write(0,'("Generating Matrix (size=",i0x,")...")')n
call psb_cdall(ictxt,desc_a,info,mg=n,parts=parts) call psb_cdall(ictxt,desc_a,info,mg=n,parts=parts)
call psb_spall(a,desc_a,info,nnz=nnz) if (info == 0) 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) if (info == 0) call psb_geall(b,desc_a,info)
call psb_geall(xv,desc_a,info) if (info == 0) call psb_geall(xv,desc_a,info)
if(info /= 0) then if (info /= 0) then
info=4010 info=4010
ch_err='allocation rout.' ch_err='allocation rout.'
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
@ -470,7 +468,7 @@ contains
! a bunch of rows per call. ! a bunch of rows per call.
! !
allocate(val(20*nbmax),irow(20*nbmax),& allocate(val(20*nbmax),irow(20*nbmax),&
&icol(20*nbmax),prv(np),stat=info) &icol(20*nbmax),stat=info)
if (info /= 0 ) then if (info /= 0 ) then
info=4000 info=4000
call psb_errpush(info,name) call psb_errpush(info,name)
@ -486,138 +484,135 @@ contains
! icol(1)=1 ! icol(1)=1
do glob_row = 1, n do glob_row = 1, n
call parts(glob_row,n,np,prv,nv) ! Figure out which rows are local to the current process:
do inv = 1, nv if (psb_is_owned(glob_row,desc_a)) then
indx_owner = prv(inv) ! local matrix pointer
if (indx_owner == iam) then element=1
! local matrix pointer ! compute gridpoint coordinates
element=1 if (mod(glob_row,(idim*idim)) == 0) then
! compute gridpoint coordinates x = glob_row/(idim*idim)
if (mod(glob_row,(idim*idim)) == 0) then else
x = glob_row/(idim*idim) x = glob_row/(idim*idim)+1
else endif
x = glob_row/(idim*idim)+1 if (mod((glob_row-(x-1)*idim*idim),idim) == 0) then
endif y = (glob_row-(x-1)*idim*idim)/idim
if (mod((glob_row-(x-1)*idim*idim),idim) == 0) then else
y = (glob_row-(x-1)*idim*idim)/idim y = (glob_row-(x-1)*idim*idim)/idim+1
else endif
y = (glob_row-(x-1)*idim*idim)/idim+1 z = glob_row-(x-1)*idim*idim-(y-1)*idim
endif ! glob_x, glob_y, glob_x coordinates
z = glob_row-(x-1)*idim*idim-(y-1)*idim glob_x=x*deltah
! glob_x, glob_y, glob_x coordinates glob_y=y*deltah
glob_x=x*deltah glob_z=z*deltah
glob_y=y*deltah
glob_z=z*deltah ! check on boundary points
zt(1) = 0.d0
! check on boundary points ! internal point: build discretization
zt(1) = 0.d0 !
! internal point: build discretization ! term depending on (x-1,y,z)
! !
! term depending on (x-1,y,z) if (x==1) then
! val(element)=-b1(glob_x,glob_y,glob_z)&
if (x==1) then & -a1(glob_x,glob_y,glob_z)
val(element)=-b1(glob_x,glob_y,glob_z)& val(element) = val(element)/(deltah*&
& -a1(glob_x,glob_y,glob_z) & deltah)
val(element) = val(element)/(deltah*& zt(1) = exp(-glob_y**2-glob_z**2)*(-val(element))
& deltah) else
zt(1) = exp(-glob_y**2-glob_z**2)*(-val(element)) val(element)=-b1(glob_x,glob_y,glob_z)&
else & -a1(glob_x,glob_y,glob_z)
val(element)=-b1(glob_x,glob_y,glob_z)& val(element) = val(element)/(deltah*&
& -a1(glob_x,glob_y,glob_z) & deltah)
val(element) = val(element)/(deltah*& icol(element)=(x-2)*idim*idim+(y-1)*idim+(z)
& deltah) element=element+1
icol(element)=(x-2)*idim*idim+(y-1)*idim+(z) endif
element=element+1 ! term depending on (x,y-1,z)
endif if (y==1) then
! term depending on (x,y-1,z) val(element)=-b2(glob_x,glob_y,glob_z)&
if (y==1) then & -a2(glob_x,glob_y,glob_z)
val(element)=-b2(glob_x,glob_y,glob_z)& val(element) = val(element)/(deltah*&
& -a2(glob_x,glob_y,glob_z) & deltah)
val(element) = val(element)/(deltah*& zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))
& deltah) else
zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element)) val(element)=-b2(glob_x,glob_y,glob_z)&
else & -a2(glob_x,glob_y,glob_z)
val(element)=-b2(glob_x,glob_y,glob_z)& val(element) = val(element)/(deltah*&
& -a2(glob_x,glob_y,glob_z) & deltah)
val(element) = val(element)/(deltah*& icol(element)=(x-1)*idim*idim+(y-2)*idim+(z)
& deltah) element=element+1
icol(element)=(x-1)*idim*idim+(y-2)*idim+(z) endif
element=element+1 ! term depending on (x,y,z-1)
endif if (z==1) then
! term depending on (x,y,z-1) val(element)=-b3(glob_x,glob_y,glob_z)&
if (z==1) then & -a3(glob_x,glob_y,glob_z)
val(element)=-b3(glob_x,glob_y,glob_z)& val(element) = val(element)/(deltah*&
& -a3(glob_x,glob_y,glob_z) & deltah)
val(element) = val(element)/(deltah*& zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))
& deltah) else
zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element)) val(element)=-b3(glob_x,glob_y,glob_z)&
else & -a3(glob_x,glob_y,glob_z)
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)
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*& val(element) = val(element)/(deltah*&
& deltah) & deltah)
icol(element)=(x-1)*idim*idim+(y-1)*idim+(z) icol(element)=(x-1)*idim*idim+(y-1)*idim+(z-1)
element=element+1 element=element+1
! term depending on (x,y,z+1) endif
if (z==idim) then ! term depending on (x,y,z)
val(element)=-b1(glob_x,glob_y,glob_z) val(element)=2*b1(glob_x,glob_y,glob_z)&
val(element) = val(element)/(deltah*& & +2*b2(glob_x,glob_y,glob_z)&
& deltah) & +2*b3(glob_x,glob_y,glob_z)&
zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element)) & +a1(glob_x,glob_y,glob_z)&
else & +a2(glob_x,glob_y,glob_z)&
val(element)=-b1(glob_x,glob_y,glob_z) & +a3(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*& val(element) = val(element)/(deltah*&
& deltah) & deltah)
icol(element)=(x-1)*idim*idim+(y-1)*idim+(z+1) icol(element)=(x-1)*idim*idim+(y-1)*idim+(z)
element=element+1 element=element+1
endif ! term depending on (x,y,z+1)
! term depending on (x,y+1,z) if (z==idim) then
if (y==idim) then val(element)=-b1(glob_x,glob_y,glob_z)
val(element)=-b2(glob_x,glob_y,glob_z) val(element) = val(element)/(deltah*&
val(element) = val(element)/(deltah*& & deltah)
& deltah) zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))
zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element)) else
else val(element)=-b1(glob_x,glob_y,glob_z)
val(element)=-b2(glob_x,glob_y,glob_z) val(element) = val(element)/(deltah*&
val(element) = val(element)/(deltah*& & deltah)
& deltah) icol(element)=(x-1)*idim*idim+(y-1)*idim+(z+1)
icol(element)=(x-1)*idim*idim+(y)*idim+(z) element=element+1
element=element+1 endif
endif ! term depending on (x,y+1,z)
! term depending on (x+1,y,z) if (y==idim) then
if (x<idim) then val(element)=-b2(glob_x,glob_y,glob_z)
val(element)=-b3(glob_x,glob_y,glob_z) val(element) = val(element)/(deltah*&
val(element) = val(element)/(deltah*& & deltah)
& deltah) zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))
icol(element)=(x)*idim*idim+(y-1)*idim+(z) else
element=element+1 val(element)=-b2(glob_x,glob_y,glob_z)
endif val(element) = val(element)/(deltah*&
irow(1:element-1)=glob_row & deltah)
ia=glob_row icol(element)=(x-1)*idim*idim+(y)*idim+(z)
element=element+1
t3 = psb_wtime() endif
call psb_spins(element-1,irow,icol,val,a,desc_a,info) ! term depending on (x+1,y,z)
if(info /= 0) exit if (x<idim) then
tins = tins + (psb_wtime()-t3) val(element)=-b3(glob_x,glob_y,glob_z)
call psb_geins(1,(/ia/),zt(1:1),b,desc_a,info) val(element) = val(element)/(deltah*&
if(info /= 0) exit & deltah)
zt(1)=0.d0 icol(element)=(x)*idim*idim+(y-1)*idim+(z)
call psb_geins(1,(/ia/),zt(1:1),xv,desc_a,info) element=element+1
if(info /= 0) exit endif
end if irow(1:element-1)=glob_row
end do 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 end do
call psb_barrier(ictxt) call psb_barrier(ictxt)
@ -649,7 +644,7 @@ contains
call psb_amx(ictxt,tasb) call psb_amx(ictxt,tasb)
if(iam == psb_root_) then 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) & a%fida(1:3)
write(*,'("-pspins time : ",es10.4)')tins write(*,'("-pspins time : ",es10.4)')tins
write(*,'("-insert time : ",es10.4)')t2 write(*,'("-insert time : ",es10.4)')t2

273
tests/pdegen/spde.f90
Unescape Escape View File

@ -189,7 +189,7 @@ program spde
call mld_precset(prec,mld_coarse_subsolve_, prectype%csbsolve,info) call mld_precset(prec,mld_coarse_subsolve_, prectype%csbsolve,info)
call mld_precset(prec,mld_coarse_mat_, prectype%cmat, info) call mld_precset(prec,mld_coarse_mat_, prectype%cmat, info)
call mld_precset(prec,mld_coarse_fillin_, prectype%cfill, info) call mld_precset(prec,mld_coarse_fillin_, prectype%cfill, info)
call mld_precset(prec,mld_coarse_iluthrs_, prectype%cthres, info) call mld_precset(prec,mld_coarse_iluthrs_, prectype%cthres, info)
call mld_precset(prec,mld_coarse_sweeps_, prectype%cjswp, info) call mld_precset(prec,mld_coarse_sweeps_, prectype%cjswp, info)
end if end if
@ -356,7 +356,6 @@ contains
return return
end subroutine get_parms end subroutine get_parms
! !
! print an error message ! print an error message
! !
@ -421,7 +420,6 @@ contains
integer :: x,y,z,ia,indx_owner integer :: x,y,z,ia,indx_owner
integer :: np, iam integer :: np, iam
integer :: element integer :: element
integer :: nv, inv
integer, allocatable :: irow(:),icol(:) integer, allocatable :: irow(:),icol(:)
real(psb_spk_), allocatable :: val(:) real(psb_spk_), allocatable :: val(:)
integer, allocatable :: prv(:) integer, allocatable :: prv(:)
@ -454,11 +452,11 @@ contains
if(iam == psb_root_) write(0,'("Generating Matrix (size=",i0x,")...")')n if(iam == psb_root_) write(0,'("Generating Matrix (size=",i0x,")...")')n
call psb_cdall(ictxt,desc_a,info,mg=n,parts=parts) call psb_cdall(ictxt,desc_a,info,mg=n,parts=parts)
call psb_spall(a,desc_a,info,nnz=nnz) if (info == 0) 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) if (info == 0) call psb_geall(b,desc_a,info)
call psb_geall(xv,desc_a,info) if (info == 0) call psb_geall(xv,desc_a,info)
if(info /= 0) then if (info /= 0) then
info=4010 info=4010
ch_err='allocation rout.' ch_err='allocation rout.'
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
@ -470,7 +468,7 @@ contains
! a bunch of rows per call. ! a bunch of rows per call.
! !
allocate(val(20*nbmax),irow(20*nbmax),& allocate(val(20*nbmax),irow(20*nbmax),&
&icol(20*nbmax),prv(np),stat=info) &icol(20*nbmax),stat=info)
if (info /= 0 ) then if (info /= 0 ) then
info=4000 info=4000
call psb_errpush(info,name) call psb_errpush(info,name)
@ -486,138 +484,135 @@ contains
! icol(1)=1 ! icol(1)=1
do glob_row = 1, n do glob_row = 1, n
call parts(glob_row,n,np,prv,nv) ! Figure out which rows are local to the current process:
do inv = 1, nv if (psb_is_owned(glob_row,desc_a)) then
indx_owner = prv(inv) ! local matrix pointer
if (indx_owner == iam) then element=1
! local matrix pointer ! compute gridpoint coordinates
element=1 if (mod(glob_row,(idim*idim)) == 0) then
! compute gridpoint coordinates x = glob_row/(idim*idim)
if (mod(glob_row,(idim*idim)) == 0) then else
x = glob_row/(idim*idim) x = glob_row/(idim*idim)+1
else endif
x = glob_row/(idim*idim)+1 if (mod((glob_row-(x-1)*idim*idim),idim) == 0) then
endif y = (glob_row-(x-1)*idim*idim)/idim
if (mod((glob_row-(x-1)*idim*idim),idim) == 0) then else
y = (glob_row-(x-1)*idim*idim)/idim y = (glob_row-(x-1)*idim*idim)/idim+1
else endif
y = (glob_row-(x-1)*idim*idim)/idim+1 z = glob_row-(x-1)*idim*idim-(y-1)*idim
endif ! glob_x, glob_y, glob_x coordinates
z = glob_row-(x-1)*idim*idim-(y-1)*idim glob_x=x*deltah
! glob_x, glob_y, glob_x coordinates glob_y=y*deltah
glob_x=x*deltah glob_z=z*deltah
glob_y=y*deltah
glob_z=z*deltah ! check on boundary points
zt(1) = 0.d0
! check on boundary points ! internal point: build discretization
zt(1) = 0.d0 !
! internal point: build discretization ! term depending on (x-1,y,z)
! !
! term depending on (x-1,y,z) if (x==1) then
! val(element)=-b1(glob_x,glob_y,glob_z)&
if (x==1) then & -a1(glob_x,glob_y,glob_z)
val(element)=-b1(glob_x,glob_y,glob_z)& val(element) = val(element)/(deltah*&
& -a1(glob_x,glob_y,glob_z) & deltah)
val(element) = val(element)/(deltah*& zt(1) = exp(-glob_y**2-glob_z**2)*(-val(element))
& deltah) else
zt(1) = exp(-glob_y**2-glob_z**2)*(-val(element)) val(element)=-b1(glob_x,glob_y,glob_z)&
else & -a1(glob_x,glob_y,glob_z)
val(element)=-b1(glob_x,glob_y,glob_z)& val(element) = val(element)/(deltah*&
& -a1(glob_x,glob_y,glob_z) & deltah)
val(element) = val(element)/(deltah*& icol(element)=(x-2)*idim*idim+(y-1)*idim+(z)
& deltah) element=element+1
icol(element)=(x-2)*idim*idim+(y-1)*idim+(z) endif
element=element+1 ! term depending on (x,y-1,z)
endif if (y==1) then
! term depending on (x,y-1,z) val(element)=-b2(glob_x,glob_y,glob_z)&
if (y==1) then & -a2(glob_x,glob_y,glob_z)
val(element)=-b2(glob_x,glob_y,glob_z)& val(element) = val(element)/(deltah*&
& -a2(glob_x,glob_y,glob_z) & deltah)
val(element) = val(element)/(deltah*& zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))
& deltah) else
zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element)) val(element)=-b2(glob_x,glob_y,glob_z)&
else & -a2(glob_x,glob_y,glob_z)
val(element)=-b2(glob_x,glob_y,glob_z)& val(element) = val(element)/(deltah*&
& -a2(glob_x,glob_y,glob_z) & deltah)
val(element) = val(element)/(deltah*& icol(element)=(x-1)*idim*idim+(y-2)*idim+(z)
& deltah) element=element+1
icol(element)=(x-1)*idim*idim+(y-2)*idim+(z) endif
element=element+1 ! term depending on (x,y,z-1)
endif if (z==1) then
! term depending on (x,y,z-1) val(element)=-b3(glob_x,glob_y,glob_z)&
if (z==1) then & -a3(glob_x,glob_y,glob_z)
val(element)=-b3(glob_x,glob_y,glob_z)& val(element) = val(element)/(deltah*&
& -a3(glob_x,glob_y,glob_z) & deltah)
val(element) = val(element)/(deltah*& zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))
& deltah) else
zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element)) val(element)=-b3(glob_x,glob_y,glob_z)&
else & -a3(glob_x,glob_y,glob_z)
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)
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*& val(element) = val(element)/(deltah*&
& deltah) & deltah)
icol(element)=(x-1)*idim*idim+(y-1)*idim+(z) icol(element)=(x-1)*idim*idim+(y-1)*idim+(z-1)
element=element+1 element=element+1
! term depending on (x,y,z+1) endif
if (z==idim) then ! term depending on (x,y,z)
val(element)=-b1(glob_x,glob_y,glob_z) val(element)=2*b1(glob_x,glob_y,glob_z)&
val(element) = val(element)/(deltah*& & +2*b2(glob_x,glob_y,glob_z)&
& deltah) & +2*b3(glob_x,glob_y,glob_z)&
zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element)) & +a1(glob_x,glob_y,glob_z)&
else & +a2(glob_x,glob_y,glob_z)&
val(element)=-b1(glob_x,glob_y,glob_z) & +a3(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*& val(element) = val(element)/(deltah*&
& deltah) & deltah)
icol(element)=(x-1)*idim*idim+(y-1)*idim+(z+1) icol(element)=(x-1)*idim*idim+(y-1)*idim+(z)
element=element+1 element=element+1
endif ! term depending on (x,y,z+1)
! term depending on (x,y+1,z) if (z==idim) then
if (y==idim) then val(element)=-b1(glob_x,glob_y,glob_z)
val(element)=-b2(glob_x,glob_y,glob_z) val(element) = val(element)/(deltah*&
val(element) = val(element)/(deltah*& & deltah)
& deltah) zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))
zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element)) else
else val(element)=-b1(glob_x,glob_y,glob_z)
val(element)=-b2(glob_x,glob_y,glob_z) val(element) = val(element)/(deltah*&
val(element) = val(element)/(deltah*& & deltah)
& deltah) icol(element)=(x-1)*idim*idim+(y-1)*idim+(z+1)
icol(element)=(x-1)*idim*idim+(y)*idim+(z) element=element+1
element=element+1 endif
endif ! term depending on (x,y+1,z)
! term depending on (x+1,y,z) if (y==idim) then
if (x<idim) then val(element)=-b2(glob_x,glob_y,glob_z)
val(element)=-b3(glob_x,glob_y,glob_z) val(element) = val(element)/(deltah*&
val(element) = val(element)/(deltah*& & deltah)
& deltah) zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))
icol(element)=(x)*idim*idim+(y-1)*idim+(z) else
element=element+1 val(element)=-b2(glob_x,glob_y,glob_z)
endif val(element) = val(element)/(deltah*&
irow(1:element-1)=glob_row & deltah)
ia=glob_row icol(element)=(x-1)*idim*idim+(y)*idim+(z)
element=element+1
t3 = psb_wtime() endif
call psb_spins(element-1,irow,icol,val,a,desc_a,info) ! term depending on (x+1,y,z)
if(info /= 0) exit if (x<idim) then
tins = tins + (psb_wtime()-t3) val(element)=-b3(glob_x,glob_y,glob_z)
call psb_geins(1,(/ia/),zt(1:1),b,desc_a,info) val(element) = val(element)/(deltah*&
if(info /= 0) exit & deltah)
zt(1)=0.d0 icol(element)=(x)*idim*idim+(y-1)*idim+(z)
call psb_geins(1,(/ia/),zt(1:1),xv,desc_a,info) element=element+1
if(info /= 0) exit endif
end if irow(1:element-1)=glob_row
end do 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 end do
call psb_barrier(ictxt) call psb_barrier(ictxt)
@ -649,7 +644,7 @@ contains
call psb_amx(ictxt,tasb) call psb_amx(ictxt,tasb)
if(iam == psb_root_) then 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) & a%fida(1:3)
write(*,'("-pspins time : ",es10.4)')tins write(*,'("-pspins time : ",es10.4)')tins
write(*,'("-insert time : ",es10.4)')t2 write(*,'("-insert time : ",es10.4)')t2

Write Preview
Loading…
Cancel
Save