fdurastante/amg4psblas
1
0
Fork 0
mirror of https://github.com/sfilippone/amg4psblas.git synced 2026-10-06 22:55:12 +00:00
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
 tests/pdegen/ppde.f90
 tests/pdegen/spde.f90

Clean matrix generation code for pdegen examples.
This commit is contained in:
Salvatore Filippone
2008-08-26 19:56:04 +00:00
parent 8782eb20df
commit f25cb7fb4d
6 changed files with 663 additions and 679 deletions
Download Patch File Download Diff File Expand all files Collapse all files
examples/pdegen/mld_dexample_1lev.f90
+152 -157
View File
@@ -81,32 +81,32 @@ program mld_dexample_1lev
implicit none
! sparse matrices
! sparse matrices
type(psb_dspmat_type) :: A
! descriptor of sparse matrices
! descriptor of sparse matrices
type(psb_desc_type):: desc_A
! preconditioner
! preconditioner
type(mld_dprec_type) :: P
! right-hand side, solution and residual vectors
! right-hand side, solution and residual vectors
real(psb_dpk_), allocatable , save :: b(:), x(:), r(:)
! solver parameters
! solver parameters
real(psb_dpk_) :: 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, resmx, resmxp
character(len=20) :: name
! initialize the parallel environment
! initialize the parallel environment
call psb_init(ictxt)
call psb_info(ictxt,iam,np)
@@ -121,11 +121,11 @@ program mld_dexample_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()
@@ -141,12 +141,12 @@ program mld_dexample_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()
@@ -161,13 +161,13 @@ program mld_dexample_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()
@@ -289,17 +289,15 @@ contains
integer, intent(out) :: pv(*)
end subroutine parts
end interface
! local variables
! local variables
type(psb_dspmat_type) :: a
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
integer :: nv, inv
integer, allocatable :: irow(:),icol(:)
real(psb_dpk_), allocatable :: val(:)
integer, allocatable :: prv(:)
! deltah dimension of each grid cell
! deltat discretization time
real(psb_dpk_) :: deltah
@@ -329,11 +327,11 @@ contains
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
call psb_geall(b,desc_a,info)
call psb_geall(xv,desc_a,info)
if(info /= 0) then
if (info == 0) call psb_spall(a,desc_a,info,nnz=nnz)
! define rhs from boundary conditions; also build initial guess
if (info == 0) call psb_geall(b,desc_a,info)
if (info == 0) call psb_geall(xv,desc_a,info)
if (info /= 0) then
info=4010
call psb_errpush(info,name)
goto 9999
@@ -344,7 +342,7 @@ contains
! a bunch of rows per call.
!
allocate(val(20*nbmax),irow(20*nbmax),&
&icol(20*nbmax),prv(np),stat=info)
&icol(20*nbmax),stat=info)
if (info /= 0 ) then
info=4000
call psb_errpush(info,name)
@@ -359,143 +357,140 @@ contains
! distribution.
do glob_row = 1, n
call parts(glob_row,n,np,prv,nv)
do inv = 1, nv
indx_owner = prv(inv)
if (indx_owner == iam) then
! local matrix pointer
element=1
! compute gridpoint coordinates
if (mod(glob_row,ipoints*ipoints) == 0) then
x = glob_row/(ipoints*ipoints)
else
x = glob_row/(ipoints*ipoints)+1
endif
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)*ipoints*ipoints)/ipoints+1
endif
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
! Figure out which rows are local to the current process:
if (psb_is_owned(glob_row,desc_a)) then
! local matrix pointer
element=1
! compute gridpoint coordinates
if (mod(glob_row,ipoints*ipoints) == 0) then
x = glob_row/(ipoints*ipoints)
else
x = glob_row/(ipoints*ipoints)+1
endif
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)*ipoints*ipoints)/ipoints+1
endif
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.d0
! internal point: build discretization
!
! term depending on (x-1,y,z)
!
if (x==1) then
val(element)=-b1(glob_x,glob_y,glob_z)&
& -a1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
zt(1) = exp(-glob_y**2-glob_z**2)*(-val(element))
else
val(element)=-b1(glob_x,glob_y,glob_z)&
& -a1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element)=(x-2)*ipoints*ipoints+(y-1)*ipoints+(z)
element=element+1
endif
! term depending on (x,y-1,z)
if (y==1) then
val(element)=-b2(glob_x,glob_y,glob_z)&
& -a2(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y-2)*ipoints+(z)
element=element+1
endif
! term depending on (x,y,z-1)
if (z==1) then
val(element)=-b3(glob_x,glob_y,glob_z)&
& -a3(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y-1)*ipoints+(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)
! check on boundary points
zt(1) = 0.d0
! internal point: build discretization
!
! term depending on (x-1,y,z)
!
if (x==1) then
val(element)=-b1(glob_x,glob_y,glob_z)&
& -a1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z)
element=element+1
! term depending on (x,y,z+1)
if (z==ipoints) then
val(element)=-b1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y-1)*ipoints+(z+1)
element=element+1
endif
! term depending on (x,y+1,z)
if (y==ipoints) then
val(element)=-b2(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y)*ipoints+(z)
element=element+1
endif
! term depending on (x+1,y,z)
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)*ipoints*ipoints+(y-1)*ipoints+(z)
element=element+1
endif
irow(1:element-1)=glob_row
ia=glob_row
zt(1) = exp(-glob_y**2-glob_z**2)*(-val(element))
else
val(element)=-b1(glob_x,glob_y,glob_z)&
& -a1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element)=(x-2)*ipoints*ipoints+(y-1)*ipoints+(z)
element=element+1
endif
! term depending on (x,y-1,z)
if (y==1) then
val(element)=-b2(glob_x,glob_y,glob_z)&
& -a2(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y-2)*ipoints+(z)
element=element+1
endif
! term depending on (x,y,z-1)
if (z==1) then
val(element)=-b3(glob_x,glob_y,glob_z)&
& -a3(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y-1)*ipoints+(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*&
& deltah)
icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z)
element=element+1
! term depending on (x,y,z+1)
if (z==ipoints) then
val(element)=-b1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y-1)*ipoints+(z+1)
element=element+1
endif
! term depending on (x,y+1,z)
if (y==ipoints) then
val(element)=-b2(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y)*ipoints+(z)
element=element+1
endif
! term depending on (x+1,y,z)
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)*ipoints*ipoints+(y-1)*ipoints+(z)
element=element+1
endif
irow(1:element-1)=glob_row
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
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
call psb_barrier(ictxt)
@@ -525,7 +520,7 @@ contains
call psb_amx(ictxt,tasb)
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)
write(*,'("-pspins time : ",es10.4)')tins
write(*,'("-insert time : ",es10.4)')t2
examples/pdegen/mld_dexample_ml.f90
+173 -178
View File
@@ -80,36 +80,36 @@ program mld_dexample_ml
implicit none
! input parameters
! input parameters
! sparse matrices
! sparse matrices
type(psb_dspmat_type) :: A
! sparse matrices descriptor
! sparse matrices descriptor
type(psb_desc_type):: desc_A
! preconditioner
! preconditioner
type(mld_dprec_type) :: P
! right-hand side, solution and residual vectors
! right-hand side, solution and residual vectors
real(psb_dpk_), allocatable , save :: b(:), x(:), r(:)
! solver and preconditioner parameters
! solver and preconditioner parameters
real(psb_dpk_) :: tol, err
integer :: itmax, iter, istop
integer :: nlev
! parallel environment parameters
! parallel environment parameters
integer :: ictxt, iam, np
! other variables
! other variables
integer :: choice
integer :: i,info,j,amatsize,descsize,precsize
integer :: idim, ierr, ircode
real(psb_dpk_) :: t1, t2, tprec, resmx, resmxp
character(len=20) :: name
! initialize the parallel environment
! initialize the parallel environment
call psb_init(ictxt)
call psb_info(ictxt,iam,np)
@@ -125,11 +125,11 @@ program mld_dexample_ml
info=0
call psb_set_errverbosity(2)
! get parameters
! get parameters
call get_parms(ictxt,choice,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()
@@ -149,40 +149,40 @@ program mld_dexample_ml
case(1)
! initialize the default multi-level preconditioner, i.e. hybrid
! Schwarz, using RAS (with overlap 1 and ILU(0) on the blocks)
! as post-smoother and 4 block-Jacobi sweeps (with UMFPACK LU
! on the blocks) as distributed coarse-level solver
! initialize the default multi-level preconditioner, i.e. hybrid
! Schwarz, using RAS (with overlap 1 and ILU(0) on the blocks)
! as post-smoother and 4 block-Jacobi sweeps (with UMFPACK LU
! on the blocks) as distributed coarse-level solver
call mld_precinit(P,'ML',info)
call mld_precinit(P,'ML',info)
case(2)
! set a three-level hybrid Schwarz preconditioner, which uses
! block Jacobi (with ILU(0) on the blocks) as post-smoother,
! a coarsest matrix replicated on the processors, and the
! LU factorization from UMFPACK as coarse-level solver
! set a three-level hybrid Schwarz preconditioner, which uses
! block Jacobi (with ILU(0) on the blocks) as post-smoother,
! a coarsest matrix replicated on the processors, and the
! LU factorization from UMFPACK as coarse-level solver
call mld_precinit(P,'ML',info,nlev=3)
call mld_precset(P,mld_smoother_type_,'BJAC',info)
call mld_precset(P,mld_coarse_mat_,'REPL',info)
call mld_precset(P,mld_coarse_solve_,'UMF',info)
call mld_precinit(P,'ML',info,nlev=3)
call mld_precset(P,mld_smoother_type_,'BJAC',info)
call mld_precset(P,mld_coarse_mat_,'REPL',info)
call mld_precset(P,mld_coarse_solve_,'UMF',info)
case(3)
! set a three-level additive Schwarz preconditioner, which uses
! RAS (with overlap 1 and ILU(0) on the blocks) as pre- and
! post-smoother, and 5 block-Jacobi sweeps (with UMFPACK LU
! on the blocks) as distributed coarsest-level solver
! set a three-level additive Schwarz preconditioner, which uses
! RAS (with overlap 1 and ILU(0) on the blocks) as pre- and
! post-smoother, and 5 block-Jacobi sweeps (with UMFPACK LU
! on the blocks) as distributed coarsest-level solver
call mld_precinit(P,'ML',info,nlev=3)
call mld_precset(P,mld_ml_type_,'ADD',info)
call mld_precset(P,mld_smoother_pos_,'TWOSIDE',info)
call mld_precset(P,mld_coarse_sweeps_,5,info)
call mld_precinit(P,'ML',info,nlev=3)
call mld_precset(P,mld_ml_type_,'ADD',info)
call mld_precset(P,mld_smoother_pos_,'TWOSIDE',info)
call mld_precset(P,mld_coarse_sweeps_,5,info)
end select
! build the preconditioner
! build the preconditioner
call psb_barrier(ictxt)
t1 = psb_wtime()
@@ -197,13 +197,13 @@ program mld_dexample_ml
goto 9999
end if
! set the solver parameters and the initial guess
! set the solver parameters and 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()
@@ -327,17 +327,15 @@ contains
integer, intent(out) :: pv(*)
end subroutine parts
end interface
! local variables
! local variables
type(psb_dspmat_type) :: a
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
integer :: nv, inv
integer, allocatable :: irow(:),icol(:)
real(psb_dpk_), allocatable :: val(:)
integer, allocatable :: prv(:)
! deltah dimension of each grid cell
! deltat discretization time
real(psb_dpk_) :: deltah
@@ -367,11 +365,11 @@ contains
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
call psb_geall(b,desc_a,info)
call psb_geall(xv,desc_a,info)
if(info /= 0) then
if (info == 0) call psb_spall(a,desc_a,info,nnz=nnz)
! define rhs from boundary conditions; also build initial guess
if (info == 0) call psb_geall(b,desc_a,info)
if (info == 0) call psb_geall(xv,desc_a,info)
if (info /= 0) then
info=4010
call psb_errpush(info,name)
goto 9999
@@ -382,7 +380,7 @@ contains
! a bunch of rows per call.
!
allocate(val(20*nbmax),irow(20*nbmax),&
&icol(20*nbmax),prv(np),stat=info)
&icol(20*nbmax),stat=info)
if (info /= 0 ) then
info=4000
call psb_errpush(info,name)
@@ -397,143 +395,140 @@ contains
! distribution.
do glob_row = 1, n
call parts(glob_row,n,np,prv,nv)
do inv = 1, nv
indx_owner = prv(inv)
if (indx_owner == iam) then
! local matrix pointer
element=1
! compute gridpoint coordinates
if (mod(glob_row,ipoints*ipoints) == 0) then
x = glob_row/(ipoints*ipoints)
else
x = glob_row/(ipoints*ipoints)+1
endif
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)*ipoints*ipoints)/ipoints+1
endif
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
! Figure out which rows are local to the current process:
if (psb_is_owned(glob_row,desc_a)) then
! local matrix pointer
element=1
! compute gridpoint coordinates
if (mod(glob_row,ipoints*ipoints) == 0) then
x = glob_row/(ipoints*ipoints)
else
x = glob_row/(ipoints*ipoints)+1
endif
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)*ipoints*ipoints)/ipoints+1
endif
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.d0
! internal point: build discretization
!
! term depending on (x-1,y,z)
!
if (x==1) then
val(element)=-b1(glob_x,glob_y,glob_z)&
& -a1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
zt(1) = exp(-glob_y**2-glob_z**2)*(-val(element))
else
val(element)=-b1(glob_x,glob_y,glob_z)&
& -a1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element)=(x-2)*ipoints*ipoints+(y-1)*ipoints+(z)
element=element+1
endif
! term depending on (x,y-1,z)
if (y==1) then
val(element)=-b2(glob_x,glob_y,glob_z)&
& -a2(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y-2)*ipoints+(z)
element=element+1
endif
! term depending on (x,y,z-1)
if (z==1) then
val(element)=-b3(glob_x,glob_y,glob_z)&
& -a3(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y-1)*ipoints+(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)
! check on boundary points
zt(1) = 0.d0
! internal point: build discretization
!
! term depending on (x-1,y,z)
!
if (x==1) then
val(element)=-b1(glob_x,glob_y,glob_z)&
& -a1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z)
element=element+1
! term depending on (x,y,z+1)
if (z==ipoints) then
val(element)=-b1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y-1)*ipoints+(z+1)
element=element+1
endif
! term depending on (x,y+1,z)
if (y==ipoints) then
val(element)=-b2(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y)*ipoints+(z)
element=element+1
endif
! term depending on (x+1,y,z)
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)*ipoints*ipoints+(y-1)*ipoints+(z)
element=element+1
endif
irow(1:element-1)=glob_row
ia=glob_row
zt(1) = exp(-glob_y**2-glob_z**2)*(-val(element))
else
val(element)=-b1(glob_x,glob_y,glob_z)&
& -a1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element)=(x-2)*ipoints*ipoints+(y-1)*ipoints+(z)
element=element+1
endif
! term depending on (x,y-1,z)
if (y==1) then
val(element)=-b2(glob_x,glob_y,glob_z)&
& -a2(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y-2)*ipoints+(z)
element=element+1
endif
! term depending on (x,y,z-1)
if (z==1) then
val(element)=-b3(glob_x,glob_y,glob_z)&
& -a3(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y-1)*ipoints+(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*&
& deltah)
icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z)
element=element+1
! term depending on (x,y,z+1)
if (z==ipoints) then
val(element)=-b1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y-1)*ipoints+(z+1)
element=element+1
endif
! term depending on (x,y+1,z)
if (y==ipoints) then
val(element)=-b2(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y)*ipoints+(z)
element=element+1
endif
! term depending on (x+1,y,z)
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)*ipoints*ipoints+(y-1)*ipoints+(z)
element=element+1
endif
irow(1:element-1)=glob_row
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
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
call psb_barrier(ictxt)
@@ -563,7 +558,7 @@ contains
call psb_amx(ictxt,tasb)
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)
write(*,'("-pspins time : ",es10.4)')tins
write(*,'("-insert time : ",es10.4)')t2
examples/pdegen/mld_sexample_1lev.f90
+137 -141
View File
@@ -300,7 +300,6 @@ contains
integer :: nv, inv
integer, allocatable :: irow(:),icol(:)
real(psb_spk_), allocatable :: val(:)
integer, allocatable :: prv(:)
! deltah dimension of each grid cell
! deltat discretization time
real(psb_spk_) :: deltah
@@ -330,11 +329,11 @@ contains
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
call psb_geall(b,desc_a,info)
call psb_geall(xv,desc_a,info)
if(info /= 0) then
if (info == 0) call psb_spall(a,desc_a,info,nnz=nnz)
! define rhs from boundary conditions; also build initial guess
if (info == 0) call psb_geall(b,desc_a,info)
if (info == 0) call psb_geall(xv,desc_a,info)
if (info /= 0) then
info=4010
call psb_errpush(info,name)
goto 9999
@@ -345,7 +344,7 @@ contains
! a bunch of rows per call.
!
allocate(val(20*nbmax),irow(20*nbmax),&
&icol(20*nbmax),prv(np),stat=info)
&icol(20*nbmax),stat=info)
if (info /= 0 ) then
info=4000
call psb_errpush(info,name)
@@ -360,143 +359,140 @@ contains
! distribution.
do glob_row = 1, n
call parts(glob_row,n,np,prv,nv)
do inv = 1, nv
indx_owner = prv(inv)
if (indx_owner == iam) then
! local matrix pointer
element=1
! compute gridpoint coordinates
if (mod(glob_row,ipoints*ipoints) == 0) then
x = glob_row/(ipoints*ipoints)
else
x = glob_row/(ipoints*ipoints)+1
endif
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)*ipoints*ipoints)/ipoints+1
endif
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
! Figure out which rows are local to the current process:
if (psb_is_owned(glob_row,desc_a)) then
! local matrix pointer
element=1
! compute gridpoint coordinates
if (mod(glob_row,ipoints*ipoints) == 0) then
x = glob_row/(ipoints*ipoints)
else
x = glob_row/(ipoints*ipoints)+1
endif
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)*ipoints*ipoints)/ipoints+1
endif
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.d0
! internal point: build discretization
!
! term depending on (x-1,y,z)
!
if (x==1) then
val(element)=-b1(glob_x,glob_y,glob_z)&
& -a1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
zt(1) = exp(-glob_y**2-glob_z**2)*(-val(element))
else
val(element)=-b1(glob_x,glob_y,glob_z)&
& -a1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element)=(x-2)*ipoints*ipoints+(y-1)*ipoints+(z)
element=element+1
endif
! term depending on (x,y-1,z)
if (y==1) then
val(element)=-b2(glob_x,glob_y,glob_z)&
& -a2(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y-2)*ipoints+(z)
element=element+1
endif
! term depending on (x,y,z-1)
if (z==1) then
val(element)=-b3(glob_x,glob_y,glob_z)&
& -a3(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y-1)*ipoints+(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)
! check on boundary points
zt(1) = 0.d0
! internal point: build discretization
!
! term depending on (x-1,y,z)
!
if (x==1) then
val(element)=-b1(glob_x,glob_y,glob_z)&
& -a1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z)
element=element+1
! term depending on (x,y,z+1)
if (z==ipoints) then
val(element)=-b1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y-1)*ipoints+(z+1)
element=element+1
endif
! term depending on (x,y+1,z)
if (y==ipoints) then
val(element)=-b2(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y)*ipoints+(z)
element=element+1
endif
! term depending on (x+1,y,z)
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)*ipoints*ipoints+(y-1)*ipoints+(z)
element=element+1
endif
irow(1:element-1)=glob_row
ia=glob_row
zt(1) = exp(-glob_y**2-glob_z**2)*(-val(element))
else
val(element)=-b1(glob_x,glob_y,glob_z)&
& -a1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element)=(x-2)*ipoints*ipoints+(y-1)*ipoints+(z)
element=element+1
endif
! term depending on (x,y-1,z)
if (y==1) then
val(element)=-b2(glob_x,glob_y,glob_z)&
& -a2(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y-2)*ipoints+(z)
element=element+1
endif
! term depending on (x,y,z-1)
if (z==1) then
val(element)=-b3(glob_x,glob_y,glob_z)&
& -a3(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y-1)*ipoints+(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*&
& deltah)
icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z)
element=element+1
! term depending on (x,y,z+1)
if (z==ipoints) then
val(element)=-b1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y-1)*ipoints+(z+1)
element=element+1
endif
! term depending on (x,y+1,z)
if (y==ipoints) then
val(element)=-b2(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y)*ipoints+(z)
element=element+1
endif
! term depending on (x+1,y,z)
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)*ipoints*ipoints+(y-1)*ipoints+(z)
element=element+1
endif
irow(1:element-1)=glob_row
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
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
call psb_barrier(ictxt)
@@ -526,7 +522,7 @@ contains
call psb_amx(ictxt,tasb)
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)
write(*,'("-pspins time : ",es10.4)')tins
write(*,'("-insert time : ",es10.4)')t2
examples/pdegen/mld_sexample_ml.f90
+173 -177
View File
@@ -80,29 +80,29 @@ program mld_sexample_ml
implicit none
! input parameters
! input parameters
! sparse matrices
! sparse matrices
type(psb_sspmat_type) :: A
! sparse matrices descriptor
! sparse matrices descriptor
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 and preconditioner parameters
! solver and preconditioner parameters
real(psb_spk_) :: tol, err
integer :: itmax, iter, istop
integer :: nlev
! parallel environment parameters
! parallel environment parameters
integer :: ictxt, iam, np
! other variables
! other variables
integer :: choice
integer :: i,info,j,amatsize,descsize,precsize
integer :: idim, ierr, ircode
@@ -110,7 +110,7 @@ program mld_sexample_ml
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)
@@ -126,11 +126,11 @@ program mld_sexample_ml
info=0
call psb_set_errverbosity(2)
! get parameters
! get parameters
call get_parms(ictxt,choice,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()
@@ -150,40 +150,40 @@ program mld_sexample_ml
case(1)
! initialize the default multi-level preconditioner, i.e. hybrid
! Schwarz, using RAS (with overlap 1 and ILU(0) on the blocks)
! as post-smoother and 4 block-Jacobi sweeps (with UMFPACK LU
! on the blocks) as distributed coarse-level solver
! initialize the default multi-level preconditioner, i.e. hybrid
! Schwarz, using RAS (with overlap 1 and ILU(0) on the blocks)
! as post-smoother and 4 block-Jacobi sweeps (with UMFPACK LU
! on the blocks) as distributed coarse-level solver
call mld_precinit(P,'ML',info)
call mld_precinit(P,'ML',info)
case(2)
! set a three-level hybrid Schwarz preconditioner, which uses
! block Jacobi (with ILU(0) on the blocks) as post-smoother,
! a coarsest matrix replicated on the processors, and the
! LU factorization from UMFPACK as coarse-level solver
! set a three-level hybrid Schwarz preconditioner, which uses
! block Jacobi (with ILU(0) on the blocks) as post-smoother,
! a coarsest matrix replicated on the processors, and the
! LU factorization from UMFPACK as coarse-level solver
call mld_precinit(P,'ML',info,nlev=3)
call mld_precset(P,mld_smoother_type_,'BJAC',info)
call mld_precset(P,mld_coarse_mat_,'REPL',info)
call mld_precset(P,mld_coarse_solve_,'UMF',info)
call mld_precinit(P,'ML',info,nlev=3)
call mld_precset(P,mld_smoother_type_,'BJAC',info)
call mld_precset(P,mld_coarse_mat_,'REPL',info)
call mld_precset(P,mld_coarse_solve_,'UMF',info)
case(3)
! set a three-level additive Schwarz preconditioner, which uses
! RAS (with overlap 1 and ILU(0) on the blocks) as pre- and
! post-smoother, and 5 block-Jacobi sweeps (with UMFPACK LU
! on the blocks) as distributed coarsest-level solver
! set a three-level additive Schwarz preconditioner, which uses
! RAS (with overlap 1 and ILU(0) on the blocks) as pre- and
! post-smoother, and 5 block-Jacobi sweeps (with UMFPACK LU
! on the blocks) as distributed coarsest-level solver
call mld_precinit(P,'ML',info,nlev=3)
call mld_precset(P,mld_ml_type_,'ADD',info)
call mld_precset(P,mld_smoother_pos_,'TWOSIDE',info)
call mld_precset(P,mld_coarse_sweeps_,5,info)
call mld_precinit(P,'ML',info,nlev=3)
call mld_precset(P,mld_ml_type_,'ADD',info)
call mld_precset(P,mld_smoother_pos_,'TWOSIDE',info)
call mld_precset(P,mld_coarse_sweeps_,5,info)
end select
! build the preconditioner
! build the preconditioner
call psb_barrier(ictxt)
t1 = psb_wtime()
@@ -198,13 +198,13 @@ program mld_sexample_ml
goto 9999
end if
! set the solver parameters and the initial guess
! set the solver parameters and 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()
@@ -328,7 +328,7 @@ 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,ipoints
@@ -338,7 +338,6 @@ contains
integer :: nv, inv
integer, allocatable :: irow(:),icol(:)
real(psb_spk_), allocatable :: val(:)
integer, allocatable :: prv(:)
! deltah dimension of each grid cell
! deltat discretization time
real(psb_spk_) :: deltah
@@ -368,11 +367,11 @@ contains
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
call psb_geall(b,desc_a,info)
call psb_geall(xv,desc_a,info)
if(info /= 0) then
if (info == 0) call psb_spall(a,desc_a,info,nnz=nnz)
! define rhs from boundary conditions; also build initial guess
if (info == 0) call psb_geall(b,desc_a,info)
if (info == 0) call psb_geall(xv,desc_a,info)
if (info /= 0) then
info=4010
call psb_errpush(info,name)
goto 9999
@@ -383,7 +382,7 @@ contains
! a bunch of rows per call.
!
allocate(val(20*nbmax),irow(20*nbmax),&
&icol(20*nbmax),prv(np),stat=info)
&icol(20*nbmax),stat=info)
if (info /= 0 ) then
info=4000
call psb_errpush(info,name)
@@ -398,143 +397,140 @@ contains
! distribution.
do glob_row = 1, n
call parts(glob_row,n,np,prv,nv)
do inv = 1, nv
indx_owner = prv(inv)
if (indx_owner == iam) then
! local matrix pointer
element=1
! compute gridpoint coordinates
if (mod(glob_row,ipoints*ipoints) == 0) then
x = glob_row/(ipoints*ipoints)
else
x = glob_row/(ipoints*ipoints)+1
endif
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)*ipoints*ipoints)/ipoints+1
endif
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
! Figure out which rows are local to the current process:
if (psb_is_owned(glob_row,desc_a)) then
! local matrix pointer
element=1
! compute gridpoint coordinates
if (mod(glob_row,ipoints*ipoints) == 0) then
x = glob_row/(ipoints*ipoints)
else
x = glob_row/(ipoints*ipoints)+1
endif
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)*ipoints*ipoints)/ipoints+1
endif
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.d0
! internal point: build discretization
!
! term depending on (x-1,y,z)
!
if (x==1) then
val(element)=-b1(glob_x,glob_y,glob_z)&
& -a1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
zt(1) = exp(-glob_y**2-glob_z**2)*(-val(element))
else
val(element)=-b1(glob_x,glob_y,glob_z)&
& -a1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element)=(x-2)*ipoints*ipoints+(y-1)*ipoints+(z)
element=element+1
endif
! term depending on (x,y-1,z)
if (y==1) then
val(element)=-b2(glob_x,glob_y,glob_z)&
& -a2(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y-2)*ipoints+(z)
element=element+1
endif
! term depending on (x,y,z-1)
if (z==1) then
val(element)=-b3(glob_x,glob_y,glob_z)&
& -a3(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y-1)*ipoints+(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)
! check on boundary points
zt(1) = 0.d0
! internal point: build discretization
!
! term depending on (x-1,y,z)
!
if (x==1) then
val(element)=-b1(glob_x,glob_y,glob_z)&
& -a1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z)
element=element+1
! term depending on (x,y,z+1)
if (z==ipoints) then
val(element)=-b1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y-1)*ipoints+(z+1)
element=element+1
endif
! term depending on (x,y+1,z)
if (y==ipoints) then
val(element)=-b2(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y)*ipoints+(z)
element=element+1
endif
! term depending on (x+1,y,z)
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)*ipoints*ipoints+(y-1)*ipoints+(z)
element=element+1
endif
irow(1:element-1)=glob_row
ia=glob_row
zt(1) = exp(-glob_y**2-glob_z**2)*(-val(element))
else
val(element)=-b1(glob_x,glob_y,glob_z)&
& -a1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element)=(x-2)*ipoints*ipoints+(y-1)*ipoints+(z)
element=element+1
endif
! term depending on (x,y-1,z)
if (y==1) then
val(element)=-b2(glob_x,glob_y,glob_z)&
& -a2(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y-2)*ipoints+(z)
element=element+1
endif
! term depending on (x,y,z-1)
if (z==1) then
val(element)=-b3(glob_x,glob_y,glob_z)&
& -a3(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y-1)*ipoints+(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*&
& deltah)
icol(element)=(x-1)*ipoints*ipoints+(y-1)*ipoints+(z)
element=element+1
! term depending on (x,y,z+1)
if (z==ipoints) then
val(element)=-b1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y-1)*ipoints+(z+1)
element=element+1
endif
! term depending on (x,y+1,z)
if (y==ipoints) then
val(element)=-b2(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
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)*ipoints*ipoints+(y)*ipoints+(z)
element=element+1
endif
! term depending on (x+1,y,z)
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)*ipoints*ipoints+(y-1)*ipoints+(z)
element=element+1
endif
irow(1:element-1)=glob_row
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
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
call psb_barrier(ictxt)
@@ -564,7 +560,7 @@ contains
call psb_amx(ictxt,tasb)
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)
write(*,'("-pspins time : ",es10.4)')tins
write(*,'("-insert time : ",es10.4)')t2