diff --git a/examples/pdegen/mld_dexample_1lev.f90 b/examples/pdegen/mld_dexample_1lev.f90 index 0a552c4d..5c1fb0f1 100644 --- a/examples/pdegen/mld_dexample_1lev.f90 +++ b/examples/pdegen/mld_dexample_1lev.f90 @@ -46,30 +46,20 @@ ! - choice = 2, hybrid three-level Schwarz preconditioner (Sec. 6.1, Fig. 3) ! - choice = 3, additive three-level Schwarz preconditioner (Sec. 6.1, Fig. 4) ! +! ! The PDE is a general second order equation in 3d ! -! b1 dd(u) b2 dd(u) b3 dd(u) a1 d(u) a2 d(u) a3 d(u) -! - ------ - ------ - ------ - ----- - ------ - ------ + a4 u = 0 +! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u) +! - ------ - ------ - ------ + ----- + ------ + ------ + c u = f ! dxdx dydy dzdz dx dy dz ! -! with Dirichlet boundary conditions, on the unit cube 0<=x,y,z<=1. -! -! Example taken from: -! C.T.Kelley -! Iterative Methods for Linear and Nonlinear Equations -! SIAM 1995 -! -! In this sample program the index space of the discretized -! computational domain is first numbered sequentially in a standard way, -! then the corresponding vector is distributed according to a BLOCK -! data distribution. +! with Dirichlet boundary conditions +! u = g ! -! Boundary conditions are set in a very simple way, by adding -! equations of the form +! on the unit cube 0<=x,y,z<=1. ! -! u(x,y) = exp(-x^2-y^2-z^2) ! -! Note that if a1=a2=a3=a4=0., the PDE is the well-known Laplace equation. +! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation. ! program mld_dexample_1lev use psb_base_mod @@ -91,7 +81,7 @@ program mld_dexample_1lev type(mld_dprec_type) :: P ! right-hand side, solution and residual vectors - real(psb_dpk_), allocatable , save :: b(:), x(:), r(:) + type(psb_d_vect_type) :: x, b, r ! solver parameters real(psb_dpk_) :: tol, err @@ -105,6 +95,7 @@ program mld_dexample_1lev integer(psb_long_int_k_) :: amatsize, precsize, descsize integer :: idim, nlev, ierr, ircode real(psb_dpk_) :: t1, t2, tprec, resmx, resmxp + character(len=5) :: afmt='CSR' character(len=20) :: name ! initialize the parallel environment @@ -137,7 +128,8 @@ program mld_dexample_1lev call psb_barrier(ictxt) t1 = psb_wtime() - call create_matrix(idim,a,b,x,desc_a,ictxt,info) + call psb_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,& + & a1,a2,a3,b1,b2,b3,c,g,info) call psb_barrier(ictxt) t2 = psb_wtime() - t1 if(info /= psb_success_) then @@ -172,7 +164,7 @@ program mld_dexample_1lev ! set the initial guess call psb_geall(x,desc_A,info) - x(:) =0.0 + call x%set(dzero) call psb_geasb(x,desc_A,info) ! solve Ax=b with preconditioned BiCGSTAB @@ -186,12 +178,12 @@ program mld_dexample_1lev call psb_amx(ictxt,t2) call psb_geall(r,desc_A,info) - r(:) =0.0 + call r%set(dzero) call psb_geasb(r,desc_A,info) call psb_geaxpby(done,b,dzero,r,desc_A,info) call psb_spmm(-done,A,x,done,r,desc_A,info) - call psb_genrm2s(resmx,r,desc_A,info) - call psb_geamaxs(resmxp,r,desc_A,info) + resmx = psb_genrm2(r,desc_A,info) + resmxp = psb_geamax(r,desc_A,info) amatsize = a%sizeof() descsize = desc_a%sizeof() @@ -259,332 +251,60 @@ contains call psb_bcast(ictxt,tol) end subroutine get_parms - - ! - ! subroutine to allocate and fill in the coefficient matrix and - ! the rhs ! - subroutine create_matrix(idim,a,b,xv,desc_a,ictxt,info) - ! - ! Discretize the partial diferential equation - ! - ! b1 dd(u) b2 dd(u) b3 dd(u) a1 d(u) a2 d(u) a3 d(u) - ! - ------ - ------ - ------ - ----- - ------ - ------ + a4 u = 0 - ! dxdx dydy dzdz dx dy dz - ! - ! with Dirichlet boundary conditions, on the unit cube 0<=x,y,z<=1. - ! - ! Boundary conditions are set in a very simple way, by adding - ! equations of the form - ! - ! u(x,y) = exp(-x^2-y^2-z^2) - ! - ! Note that if a1=a2=a3=a4=0., the PDE is the well-known Laplace equation. - ! - use psb_base_mod - implicit none - integer :: idim - integer, parameter :: nb=20 - real(psb_dpk_), allocatable :: b(:),xv(:) - type(psb_desc_type) :: desc_a - integer :: ictxt, info - character :: afmt*5 - type(psb_dspmat_type) :: a - real(psb_dpk_) :: zt(nb),x,y,z - integer :: m,n,nnz,glob_row,nlr,i,ii,ib,k - integer :: ix,iy,iz,ia,indx_owner, ipoints - integer :: np, iam, nr, nt - integer :: element - integer, allocatable :: irow(:),icol(:),myidx(:) - real(psb_dpk_), allocatable :: val(:) - ! deltah dimension of each grid cell - ! deltat discretization time - real(psb_dpk_) :: deltah, deltah2 - real(psb_dpk_),parameter :: rhs=0.d0,one=1.d0,zero=0.d0 - real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen - real(psb_dpk_) :: a1, a2, a3, a4, b1, b2, b3 - external :: a1, a2, a3, a4, b1, b2, b3 - integer :: err_act - - character(len=20) :: name, ch_err - - info = psb_success_ - name = 'create_matrix' - call psb_erractionsave(err_act) - - call psb_info(ictxt, iam, np) - - deltah = 1.d0/(idim-1) - deltah2 = deltah*deltah - - ! initialize array descriptor and sparse matrix storage. provide an - ! estimate of the number of non zeroes - - ipoints=idim-2 - m = ipoints*ipoints*ipoints - n = m - nnz = ((n*9)/(np)) - if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n - - ! - ! Using a simple BLOCK distribution. - ! - nt = (m+np-1)/np - nr = max(0,min(nt,m-(iam*nt))) - - nt = nr - call psb_sum(ictxt,nt) - if (nt /= m) write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m - call psb_barrier(ictxt) - t0 = psb_wtime() - call psb_cdall(ictxt,desc_a,info,nl=nr) - if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz) - ! define rhs from boundary conditions; also build initial guess - if (info == psb_success_) call psb_geall(b,desc_a,info) - if (info == psb_success_) call psb_geall(xv,desc_a,info) - nlr = psb_cd_get_local_rows(desc_a) - call psb_barrier(ictxt) - talc = psb_wtime()-t0 - - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='allocation rout.' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - ! we build an auxiliary matrix consisting of one row at a - ! time; just a small matrix. might be extended to generate - ! a bunch of rows per call. - ! - allocate(val(20*nb),irow(20*nb),& - &icol(20*nb),myidx(nlr),stat=info) - if (info /= psb_success_ ) then - info=psb_err_alloc_dealloc_ - call psb_errpush(info,name) - goto 9999 - endif - - do i=1,nlr - myidx(i) = i - end do - - - call psb_loc_to_glob(myidx,desc_a,info) - - ! loop over rows belonging to current process in a block - ! distribution. - - call psb_barrier(ictxt) - t1 = psb_wtime() - do ii=1, nlr,nb - ib = min(nb,nlr-ii+1) - element = 1 - do k=1,ib - i=ii+k-1 - ! local matrix pointer - glob_row=myidx(i) - ! compute gridpoint coordinates - if (mod(glob_row,ipoints*ipoints) == 0) then - ix = glob_row/(ipoints*ipoints) - else - ix = glob_row/(ipoints*ipoints)+1 - endif - if (mod((glob_row-(ix-1)*ipoints*ipoints),ipoints) == 0) then - iy = (glob_row-(ix-1)*ipoints*ipoints)/ipoints - else - iy = (glob_row-(ix-1)*ipoints*ipoints)/ipoints+1 - endif - iz = glob_row-(ix-1)*ipoints*ipoints-(iy-1)*ipoints - ! x, y, x coordinates - x=ix*deltah - y=iy*deltah - z=iz*deltah - - ! check on boundary points - zt(k) = 0.d0 - ! internal point: build discretization - ! - ! term depending on (x-1,y,z) - ! - if (ix == 1) then - val(element) = -b1(x,y,z)/deltah2-a1(x,y,z)/deltah - zt(k) = exp(-x**2-y**2-z**2)*(-val(element)) - else - val(element) = -b1(x,y,z)/deltah2-a1(x,y,z)/deltah - icol(element) = (ix-2)*ipoints*ipoints+(iy-1)*ipoints+(iz) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x,y-1,z) - if (iy == 1) then - val(element) = -b2(x,y,z)/deltah2-a2(x,y,z)/deltah - zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element) = -b2(x,y,z)/deltah2-a2(x,y,z)/deltah - icol(element) = (ix-1)*ipoints*ipoints+(iy-2)*ipoints+(iz) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x,y,z-1) - if (iz == 1) then - val(element)=-b3(x,y,z)/deltah2-a3(x,y,z)/deltah - zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element)=-b3(x,y,z)/deltah2-a3(x,y,z)/deltah - icol(element) = (ix-1)*ipoints*ipoints+(iy-1)*ipoints+(iz-1) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x,y,z) - val(element)=(2*b1(x,y,z) + 2*b2(x,y,z) + 2*b3(x,y,z))/deltah2& - & + (a1(x,y,z) + a2(x,y,z) + a3(x,y,z)+ a4(x,y,z))/deltah - icol(element) = (ix-1)*ipoints*ipoints+(iy-1)*ipoints+(iz) - irow(element) = glob_row - element = element+1 - ! term depending on (x,y,z+1) - if (iz == ipoints) then - val(element)=-b1(x,y,z)/deltah2 - zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element)=-b1(x,y,z)/deltah2 - icol(element) = (ix-1)*ipoints*ipoints+(iy-1)*ipoints+(iz+1) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x,y+1,z) - if (iy == ipoints) then - val(element)=-b2(x,y,z)/deltah2 - zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element)=-b2(x,y,z)/deltah2 - icol(element) = (ix-1)*ipoints*ipoints+(iy)*ipoints+(iz) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x+1,y,z) - if (ix==ipoints) then - val(element)=-b3(x,y,z)/deltah2 - zt(k) = exp(-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element)=-b3(x,y,z)/deltah2 - icol(element) = (ix)*ipoints*ipoints+(iy-1)*ipoints+(iz) - irow(element) = glob_row - element = element+1 - endif - - end do - call psb_spins(element-1,irow,icol,val,a,desc_a,info) - if(info /= psb_success_) exit - call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),b,desc_a,info) - if(info /= psb_success_) exit - zt(:)=0.d0 - call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info) - if(info /= psb_success_) exit - end do - - tgen = psb_wtime()-t1 - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name) - goto 9999 + ! functions parametrizing the differential equation + ! + function b1(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: b1 + real(psb_dpk_), intent(in) :: x,y,z + b1=1.d0/sqrt(3.d0) + end function b1 + function b2(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: b2 + real(psb_dpk_), intent(in) :: x,y,z + b2=1.d0/sqrt(3.d0) + end function b2 + function b3(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: b3 + real(psb_dpk_), intent(in) :: x,y,z + b3=1.d0/sqrt(3.d0) + end function b3 + function c(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: c + real(psb_dpk_), intent(in) :: x,y,z + c=0.d0 + end function c + function a1(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: a1 + real(psb_dpk_), intent(in) :: x,y,z + a1=1.d0/80 + end function a1 + function a2(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: a2 + real(psb_dpk_), intent(in) :: x,y,z + a2=1.d0/80 + end function a2 + function a3(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: a3 + real(psb_dpk_), intent(in) :: x,y,z + a3=1.d0/80 + end function a3 + function g(x,y,z) + use psb_base_mod, only : psb_dpk_, done + 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 - - deallocate(val,irow,icol) - - call psb_barrier(ictxt) - t1 = psb_wtime() - call psb_cdasb(desc_a,info) - if (info == psb_success_) & - & call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_) - call psb_barrier(ictxt) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name) - goto 9999 - end if - call psb_geasb(b,desc_a,info) - call psb_geasb(xv,desc_a,info) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name) - goto 9999 - end if - tasb = psb_wtime()-t1 - call psb_barrier(ictxt) - ttot = psb_wtime() - t0 - - call psb_amx(ictxt,talc) - call psb_amx(ictxt,tgen) - call psb_amx(ictxt,tasb) - call psb_amx(ictxt,ttot) - if(iam == psb_root_) then - write(*,'("The matrix has been generated and assembled in ",a3," format.")')& - & a%get_fmt() - write(*,'("-allocation time : ",es12.5)') talc - write(*,'("-coeff. gen. time : ",es12.5)') tgen - write(*,'("-assembly time : ",es12.5)') tasb - write(*,'("-total time : ",es12.5)') ttot - - end if - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error(ictxt) - return - end if - return - end subroutine create_matrix + end function g end program mld_dexample_1lev -! -! functions parametrizing the differential equation -! -function a1(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: a1 - real(psb_dpk_) :: x,y,z - !a1=1.d0 - a1=0.d0 -end function a1 -function a2(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: a2 - real(psb_dpk_) :: x,y,z - !a2=2.d1*y - a2=0.d0 -end function a2 -function a3(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: a3 - real(psb_dpk_) :: x,y,z - !a3=1.d0 - a3=0.d0 -end function a3 -function a4(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: a4 - real(psb_dpk_) :: x,y,z - !a4=1.d0 - a4=0.d0 -end function a4 -function b1(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: b1 - real(psb_dpk_) :: x,y,z - b1=1.d0 -end function b1 -function b2(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: b2 - real(psb_dpk_) :: x,y,z - b2=1.d0 -end function b2 -function b3(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: b3 - real(psb_dpk_) :: x,y,z - b3=1.d0 -end function b3 diff --git a/examples/pdegen/mld_dexample_ml.f90 b/examples/pdegen/mld_dexample_ml.f90 index 79f5d682..38c196ff 100644 --- a/examples/pdegen/mld_dexample_ml.f90 +++ b/examples/pdegen/mld_dexample_ml.f90 @@ -48,29 +48,23 @@ ! ! The PDE is a general second order equation in 3d ! -! b1 dd(u) b2 dd(u) b3 dd(u) a1 d(u) a2 d(u) a3 d(u) -! - ------ - ------ - ------ - ----- - ------ - ------ + a4 u = 0 +! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u) +! - ------ - ------ - ------ + ----- + ------ + ------ + c u = f ! dxdx dydy dzdz dx dy dz ! -! with Dirichlet boundary conditions, on the unit cube 0<=x,y,z<=1. +! with Dirichlet boundary conditions +! u = g ! -! Example taken from: -! C.T.Kelley -! Iterative Methods for Linear and Nonlinear Equations -! SIAM 1995 +! on the unit cube 0<=x,y,z<=1. +! +! +! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation. ! ! In this sample program the index space of the discretized ! computational domain is first numbered sequentially in a standard way, ! then the corresponding vector is distributed according to a BLOCK ! data distribution. ! -! Boundary conditions are set in a very simple way, by adding -! equations of the form -! -! u(x,y) = exp(-x^2-y^2-z^2) -! -! Note that if a1=a2=a3=a4=0., the PDE is the well-known Laplace equation. -! program mld_dexample_ml use psb_base_mod use mld_prec_mod @@ -92,7 +86,7 @@ program mld_dexample_ml type(mld_dprec_type) :: P ! right-hand side, solution and residual vectors - real(psb_dpk_), allocatable , save :: b(:), x(:), r(:) + type(psb_d_vect_type) :: x, b, r ! solver and preconditioner parameters real(psb_dpk_) :: tol, err @@ -108,6 +102,7 @@ program mld_dexample_ml integer(psb_long_int_k_) :: amatsize, precsize, descsize integer :: idim, ierr, ircode real(psb_dpk_) :: t1, t2, tprec, resmx, resmxp + character(len=5) :: afmt='CSR' character(len=20) :: name ! initialize the parallel environment @@ -141,7 +136,8 @@ program mld_dexample_ml call psb_barrier(ictxt) t1 = psb_wtime() - call create_matrix(idim,a,b,x,desc_a,ictxt,info) + call psb_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,& + & a1,a2,a3,b1,b2,b3,c,g,info) call psb_barrier(ictxt) t2 = psb_wtime() - t1 if(info /= psb_success_) then @@ -208,7 +204,7 @@ program mld_dexample_ml ! set the solver parameters and the initial guess call psb_geall(x,desc_A,info) - x(:) =0.0 + call x%set(dzero) call psb_geasb(x,desc_A,info) ! solve Ax=b with preconditioned BiCGSTAB @@ -222,12 +218,12 @@ program mld_dexample_ml call psb_amx(ictxt,t2) call psb_geall(r,desc_A,info) - r(:) =0.0 + call r%set(dzero) call psb_geasb(r,desc_A,info) call psb_geaxpby(done,b,dzero,r,desc_A,info) call psb_spmm(-done,A,x,done,r,desc_A,info) - call psb_genrm2s(resmx,r,desc_A,info) - call psb_geamaxs(resmxp,r,desc_A,info) + resmx = psb_genrm2(r,desc_A,info) + resmxp = psb_geamax(r,desc_A,info) amatsize = a%sizeof() descsize = desc_a%sizeof() @@ -298,331 +294,61 @@ contains end subroutine get_parms - ! - ! subroutine to allocate and fill in the coefficient matrix and - ! the rhs - ! - subroutine create_matrix(idim,a,b,xv,desc_a,ictxt,info) - ! - ! Discretize the partial diferential equation - ! - ! b1 dd(u) b2 dd(u) b3 dd(u) a1 d(u) a2 d(u) a3 d(u) - ! - ------ - ------ - ------ - ----- - ------ - ------ + a4 u = 0 - ! dxdx dydy dzdz dx dy dz - ! - ! with Dirichlet boundary conditions, on the unit cube 0<=x,y,z<=1. - ! - ! Boundary conditions are set in a very simple way, by adding - ! equations of the form - ! - ! u(x,y) = exp(-x^2-y^2-z^2) - ! - ! Note that if a1=a2=a3=a4=0., the PDE is the well-known Laplace equation. - ! - use psb_base_mod - implicit none - integer :: idim - integer, parameter :: nb=20 - real(psb_dpk_), allocatable :: b(:),xv(:) - type(psb_desc_type) :: desc_a - integer :: ictxt, info - character :: afmt*5 - type(psb_dspmat_type) :: a - real(psb_dpk_) :: zt(nb),x,y,z - integer :: m,n,nnz,glob_row,nlr,i,ii,ib,k - integer :: ix,iy,iz,ia,indx_owner, ipoints - integer :: np, iam, nr, nt - integer :: element - integer, allocatable :: irow(:),icol(:),myidx(:) - real(psb_dpk_), allocatable :: val(:) - ! deltah dimension of each grid cell - ! deltat discretization time - real(psb_dpk_) :: deltah, deltah2 - real(psb_dpk_),parameter :: rhs=0.d0,one=1.d0,zero=0.d0 - real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen - real(psb_dpk_) :: a1, a2, a3, a4, b1, b2, b3 - external :: a1, a2, a3, a4, b1, b2, b3 - integer :: err_act - - character(len=20) :: name, ch_err - - info = psb_success_ - name = 'create_matrix' - call psb_erractionsave(err_act) - - call psb_info(ictxt, iam, np) - - deltah = 1.d0/(idim-1) - deltah2 = deltah*deltah - - ! initialize array descriptor and sparse matrix storage. provide an - ! estimate of the number of non zeroes - - ipoints=idim-2 - m = ipoints*ipoints*ipoints - n = m - nnz = ((n*9)/(np)) - if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n - - ! - ! Using a simple BLOCK distribution. - ! - nt = (m+np-1)/np - nr = max(0,min(nt,m-(iam*nt))) - - nt = nr - call psb_sum(ictxt,nt) - if (nt /= m) write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m - call psb_barrier(ictxt) - t0 = psb_wtime() - call psb_cdall(ictxt,desc_a,info,nl=nr) - if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz) - ! define rhs from boundary conditions; also build initial guess - if (info == psb_success_) call psb_geall(b,desc_a,info) - if (info == psb_success_) call psb_geall(xv,desc_a,info) - nlr = psb_cd_get_local_rows(desc_a) - call psb_barrier(ictxt) - talc = psb_wtime()-t0 - - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='allocation rout.' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - ! we build an auxiliary matrix consisting of one row at a - ! time; just a small matrix. might be extended to generate - ! a bunch of rows per call. - ! - allocate(val(20*nb),irow(20*nb),& - &icol(20*nb),myidx(nlr),stat=info) - if (info /= psb_success_ ) then - info=psb_err_alloc_dealloc_ - call psb_errpush(info,name) - goto 9999 - endif - - do i=1,nlr - myidx(i) = i - end do - - - call psb_loc_to_glob(myidx,desc_a,info) - - ! loop over rows belonging to current process in a block - ! distribution. - - call psb_barrier(ictxt) - t1 = psb_wtime() - do ii=1, nlr,nb - ib = min(nb,nlr-ii+1) - element = 1 - do k=1,ib - i=ii+k-1 - ! local matrix pointer - glob_row=myidx(i) - ! compute gridpoint coordinates - if (mod(glob_row,ipoints*ipoints) == 0) then - ix = glob_row/(ipoints*ipoints) - else - ix = glob_row/(ipoints*ipoints)+1 - endif - if (mod((glob_row-(ix-1)*ipoints*ipoints),ipoints) == 0) then - iy = (glob_row-(ix-1)*ipoints*ipoints)/ipoints - else - iy = (glob_row-(ix-1)*ipoints*ipoints)/ipoints+1 - endif - iz = glob_row-(ix-1)*ipoints*ipoints-(iy-1)*ipoints - ! x, y, x coordinates - x=ix*deltah - y=iy*deltah - z=iz*deltah - - ! check on boundary points - zt(k) = 0.d0 - ! internal point: build discretization - ! - ! term depending on (x-1,y,z) - ! - if (ix == 1) then - val(element) = -b1(x,y,z)/deltah2-a1(x,y,z)/deltah - zt(k) = exp(-x**2-y**2-z**2)*(-val(element)) - else - val(element) = -b1(x,y,z)/deltah2-a1(x,y,z)/deltah - icol(element) = (ix-2)*ipoints*ipoints+(iy-1)*ipoints+(iz) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x,y-1,z) - if (iy == 1) then - val(element) = -b2(x,y,z)/deltah2-a2(x,y,z)/deltah - zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element) = -b2(x,y,z)/deltah2-a2(x,y,z)/deltah - icol(element) = (ix-1)*ipoints*ipoints+(iy-2)*ipoints+(iz) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x,y,z-1) - if (iz == 1) then - val(element)=-b3(x,y,z)/deltah2-a3(x,y,z)/deltah - zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element)=-b3(x,y,z)/deltah2-a3(x,y,z)/deltah - icol(element) = (ix-1)*ipoints*ipoints+(iy-1)*ipoints+(iz-1) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x,y,z) - val(element)=(2*b1(x,y,z) + 2*b2(x,y,z) + 2*b3(x,y,z))/deltah2& - & + (a1(x,y,z) + a2(x,y,z) + a3(x,y,z)+ a4(x,y,z))/deltah - icol(element) = (ix-1)*ipoints*ipoints+(iy-1)*ipoints+(iz) - irow(element) = glob_row - element = element+1 - ! term depending on (x,y,z+1) - if (iz == ipoints) then - val(element)=-b1(x,y,z)/deltah2 - zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element)=-b1(x,y,z)/deltah2 - icol(element) = (ix-1)*ipoints*ipoints+(iy-1)*ipoints+(iz+1) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x,y+1,z) - if (iy == ipoints) then - val(element)=-b2(x,y,z)/deltah2 - zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element)=-b2(x,y,z)/deltah2 - icol(element) = (ix-1)*ipoints*ipoints+(iy)*ipoints+(iz) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x+1,y,z) - if (ix==ipoints) then - val(element)=-b3(x,y,z)/deltah2 - zt(k) = exp(-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element)=-b3(x,y,z)/deltah2 - icol(element) = (ix)*ipoints*ipoints+(iy-1)*ipoints+(iz) - irow(element) = glob_row - element = element+1 - endif - - end do - call psb_spins(element-1,irow,icol,val,a,desc_a,info) - if(info /= psb_success_) exit - call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),b,desc_a,info) - if(info /= psb_success_) exit - zt(:)=0.d0 - call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info) - if(info /= psb_success_) exit - end do - - tgen = psb_wtime()-t1 - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name) - goto 9999 - end if - - deallocate(val,irow,icol) - - call psb_barrier(ictxt) - t1 = psb_wtime() - call psb_cdasb(desc_a,info) - if (info == psb_success_) & - & call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_) - call psb_barrier(ictxt) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name) - goto 9999 - end if - call psb_geasb(b,desc_a,info) - call psb_geasb(xv,desc_a,info) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name) - goto 9999 - end if - tasb = psb_wtime()-t1 - call psb_barrier(ictxt) - ttot = psb_wtime() - t0 - - call psb_amx(ictxt,talc) - call psb_amx(ictxt,tgen) - call psb_amx(ictxt,tasb) - call psb_amx(ictxt,ttot) - if(iam == psb_root_) then - write(*,'("The matrix has been generated and assembled in ",a3," format.")')& - & a%get_fmt() - write(*,'("-allocation time : ",es12.5)') talc - write(*,'("-coeff. gen. time : ",es12.5)') tgen - write(*,'("-assembly time : ",es12.5)') tasb - write(*,'("-total time : ",es12.5)') ttot - - end if - call psb_erractionrestore(err_act) - return -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error(ictxt) - return + ! + ! functions parametrizing the differential equation + ! + function b1(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: b1 + real(psb_dpk_), intent(in) :: x,y,z + b1=1.d0/sqrt(3.d0) + end function b1 + function b2(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: b2 + real(psb_dpk_), intent(in) :: x,y,z + b2=1.d0/sqrt(3.d0) + end function b2 + function b3(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: b3 + real(psb_dpk_), intent(in) :: x,y,z + b3=1.d0/sqrt(3.d0) + end function b3 + function c(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: c + real(psb_dpk_), intent(in) :: x,y,z + c=0.d0 + end function c + function a1(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: a1 + real(psb_dpk_), intent(in) :: x,y,z + a1=1.d0/80 + end function a1 + function a2(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: a2 + real(psb_dpk_), intent(in) :: x,y,z + a2=1.d0/80 + end function a2 + function a3(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: a3 + real(psb_dpk_), intent(in) :: x,y,z + a3=1.d0/80 + end function a3 + function g(x,y,z) + use psb_base_mod, only : psb_dpk_, done + 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 - return - end subroutine create_matrix + end function g end program mld_dexample_ml -! -! functions parametrizing the differential equation -! -function a1(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: a1 - real(psb_dpk_) :: x,y,z -! a1=1.d0 - a1=0.d0 -end function a1 -function a2(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: a2 - real(psb_dpk_) :: x,y,z -! a2=2.d1*y - a2=0.d0 -end function a2 -function a3(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: a3 - real(psb_dpk_) :: x,y,z -! a3=1.d0 - a3=0.d0 -end function a3 -function a4(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: a4 - real(psb_dpk_) :: x,y,z -! a4=1.d0 - a4=0.d0 -end function a4 -function b1(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: b1 - real(psb_dpk_) :: x,y,z - b1=1.d0 -end function b1 -function b2(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: b2 - real(psb_dpk_) :: x,y,z - b2=1.d0 -end function b2 -function b3(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: b3 - real(psb_dpk_) :: x,y,z - b3=1.d0 -end function b3 diff --git a/examples/pdegen/mld_sexample_1lev.f90 b/examples/pdegen/mld_sexample_1lev.f90 index 97b36820..725ed4db 100644 --- a/examples/pdegen/mld_sexample_1lev.f90 +++ b/examples/pdegen/mld_sexample_1lev.f90 @@ -48,29 +48,23 @@ ! ! The PDE is a general second order equation in 3d ! -! b1 dd(u) b2 dd(u) b3 dd(u) a1 d(u) a2 d(u) a3 d(u) -! - ------ - ------ - ------ - ----- - ------ - ------ + a4 u = 0 +! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u) +! - ------ - ------ - ------ + ----- + ------ + ------ + c u = f ! dxdx dydy dzdz dx dy dz ! -! with Dirichlet boundary conditions, on the unit cube 0<=x,y,z<=1. +! with Dirichlet boundary conditions +! u = g ! -! Example taken from: -! C.T.Kelley -! Iterative Methods for Linear and Nonlinear Equations -! SIAM 1995 +! on the unit cube 0<=x,y,z<=1. +! +! +! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation. ! ! In this sample program the index space of the discretized ! computational domain is first numbered sequentially in a standard way, ! then the corresponding vector is distributed according to a BLOCK ! data distribution. ! -! Boundary conditions are set in a very simple way, by adding -! equations of the form -! -! u(x,y) = exp(-x^2-y^2-z^2) -! -! Note that if a1=a2=a3=a4=0., the PDE is the well-known Laplace equation. -! program mld_sexample_1lev use psb_base_mod use mld_prec_mod @@ -91,7 +85,7 @@ program mld_sexample_1lev type(mld_sprec_type) :: P ! right-hand side, solution and residual vectors - real(psb_spk_), allocatable , save :: b(:), x(:), r(:) + type(psb_s_vect_type) :: x, b, r ! solver parameters real(psb_spk_) :: tol, err @@ -106,6 +100,7 @@ program mld_sexample_1lev integer :: idim, nlev, ierr, ircode real(psb_dpk_) :: t1, t2, tprec real(psb_spk_) :: resmx, resmxp + character(len=5) :: afmt='CSR' character(len=20) :: name ! initialize the parallel environment @@ -138,7 +133,8 @@ program mld_sexample_1lev call psb_barrier(ictxt) t1 = psb_wtime() - call create_matrix(idim,a,b,x,desc_a,ictxt,info) + call psb_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,& + & a1,a2,a3,b1,b2,b3,c,g,info) call psb_barrier(ictxt) t2 = psb_wtime() - t1 if(info /= psb_success_) then @@ -173,7 +169,7 @@ program mld_sexample_1lev ! set the initial guess call psb_geall(x,desc_A,info) - x(:) =0.0 + call x%set(szero) call psb_geasb(x,desc_A,info) ! solve Ax=b with preconditioned BiCGSTAB @@ -187,12 +183,12 @@ program mld_sexample_1lev call psb_amx(ictxt,t2) call psb_geall(r,desc_A,info) - r(:) =0.0 + call r%set(szero) call psb_geasb(r,desc_A,info) call psb_geaxpby(sone,b,szero,r,desc_A,info) call psb_spmm(-sone,A,x,sone,r,desc_A,info) - call psb_genrm2s(resmx,r,desc_A,info) - call psb_geamaxs(resmxp,r,desc_A,info) + resmx = psb_genrm2(r,desc_A,info) + resmxp = psb_geamax(r,desc_A,info) amatsize = a%sizeof() descsize = desc_a%sizeof() @@ -260,332 +256,60 @@ contains call psb_bcast(ictxt,tol) end subroutine get_parms - - ! - ! subroutine to allocate and fill in the coefficient matrix and - ! the rhs ! - subroutine create_matrix(idim,a,b,xv,desc_a,ictxt,info) - ! - ! Discretize the partial diferential equation - ! - ! b1 dd(u) b2 dd(u) b3 dd(u) a1 d(u) a2 d(u) a3 d(u) - ! - ------ - ------ - ------ - ----- - ------ - ------ + a4 u = 0 - ! dxdx dydy dzdz dx dy dz - ! - ! with Dirichlet boundary conditions, on the unit cube 0<=x,y,z<=1. - ! - ! Boundary conditions are set in a very simple way, by adding - ! equations of the form - ! - ! u(x,y) = exp(-x^2-y^2-z^2) - ! - ! Note that if a1=a2=a3=a4=0., the PDE is the well-known Laplace equation. - ! - use psb_base_mod - implicit none - integer :: idim - integer, parameter :: nb=20 - real(psb_spk_), allocatable :: b(:),xv(:) - type(psb_desc_type) :: desc_a - integer :: ictxt, info - character :: afmt*5 - type(psb_sspmat_type) :: a - real(psb_spk_) :: zt(nb),x,y,z - integer :: m,n,nnz,glob_row,nlr,i,ii,ib,k - integer :: ix,iy,iz,ia,indx_owner, ipoints - integer :: np, iam, nr, nt - integer :: element - integer, allocatable :: irow(:),icol(:),myidx(:) - real(psb_spk_), allocatable :: val(:) - ! deltah dimension of each grid cell - ! deltat discretization time - real(psb_spk_) :: deltah, deltah2 - real(psb_spk_),parameter :: rhs=0.0,one=1.0,zero=0.0 - real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen - real(psb_spk_) :: a1, a2, a3, a4, b1, b2, b3 - external :: a1, a2, a3, a4, b1, b2, b3 - integer :: err_act - - character(len=20) :: name, ch_err - - info = psb_success_ - name = 'create_matrix' - call psb_erractionsave(err_act) - - call psb_info(ictxt, iam, np) - - deltah = 1.d0/(idim-1) - deltah2 = deltah*deltah - - ! initialize array descriptor and sparse matrix storage. provide an - ! estimate of the number of non zeroes - - ipoints=idim-2 - m = ipoints*ipoints*ipoints - n = m - nnz = ((n*9)/(np)) - if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n - - ! - ! Using a simple BLOCK distribution. - ! - nt = (m+np-1)/np - nr = max(0,min(nt,m-(iam*nt))) - - nt = nr - call psb_sum(ictxt,nt) - if (nt /= m) write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m - call psb_barrier(ictxt) - t0 = psb_wtime() - call psb_cdall(ictxt,desc_a,info,nl=nr) - if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz) - ! define rhs from boundary conditions; also build initial guess - if (info == psb_success_) call psb_geall(b,desc_a,info) - if (info == psb_success_) call psb_geall(xv,desc_a,info) - nlr = psb_cd_get_local_rows(desc_a) - call psb_barrier(ictxt) - talc = psb_wtime()-t0 - - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='allocation rout.' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - ! we build an auxiliary matrix consisting of one row at a - ! time; just a small matrix. might be extended to generate - ! a bunch of rows per call. - ! - allocate(val(20*nb),irow(20*nb),& - &icol(20*nb),myidx(nlr),stat=info) - if (info /= psb_success_ ) then - info=psb_err_alloc_dealloc_ - call psb_errpush(info,name) - goto 9999 - endif - - do i=1,nlr - myidx(i) = i - end do - - - call psb_loc_to_glob(myidx,desc_a,info) - - ! loop over rows belonging to current process in a block - ! distribution. - - call psb_barrier(ictxt) - t1 = psb_wtime() - do ii=1, nlr,nb - ib = min(nb,nlr-ii+1) - element = 1 - do k=1,ib - i=ii+k-1 - ! local matrix pointer - glob_row=myidx(i) - ! compute gridpoint coordinates - if (mod(glob_row,ipoints*ipoints) == 0) then - ix = glob_row/(ipoints*ipoints) - else - ix = glob_row/(ipoints*ipoints)+1 - endif - if (mod((glob_row-(ix-1)*ipoints*ipoints),ipoints) == 0) then - iy = (glob_row-(ix-1)*ipoints*ipoints)/ipoints - else - iy = (glob_row-(ix-1)*ipoints*ipoints)/ipoints+1 - endif - iz = glob_row-(ix-1)*ipoints*ipoints-(iy-1)*ipoints - ! x, y, x coordinates - x=ix*deltah - y=iy*deltah - z=iz*deltah - - ! check on boundary points - zt(k) = 0.d0 - ! internal point: build discretization - ! - ! term depending on (x-1,y,z) - ! - if (ix == 1) then - val(element) = -b1(x,y,z)/deltah2-a1(x,y,z)/deltah - zt(k) = exp(-x**2-y**2-z**2)*(-val(element)) - else - val(element) = -b1(x,y,z)/deltah2-a1(x,y,z)/deltah - icol(element) = (ix-2)*ipoints*ipoints+(iy-1)*ipoints+(iz) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x,y-1,z) - if (iy == 1) then - val(element) = -b2(x,y,z)/deltah2-a2(x,y,z)/deltah - zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element) = -b2(x,y,z)/deltah2-a2(x,y,z)/deltah - icol(element) = (ix-1)*ipoints*ipoints+(iy-2)*ipoints+(iz) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x,y,z-1) - if (iz == 1) then - val(element)=-b3(x,y,z)/deltah2-a3(x,y,z)/deltah - zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element)=-b3(x,y,z)/deltah2-a3(x,y,z)/deltah - icol(element) = (ix-1)*ipoints*ipoints+(iy-1)*ipoints+(iz-1) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x,y,z) - val(element)=(2*b1(x,y,z) + 2*b2(x,y,z) + 2*b3(x,y,z))/deltah2& - & + (a1(x,y,z) + a2(x,y,z) + a3(x,y,z)+ a4(x,y,z))/deltah - icol(element) = (ix-1)*ipoints*ipoints+(iy-1)*ipoints+(iz) - irow(element) = glob_row - element = element+1 - ! term depending on (x,y,z+1) - if (iz == ipoints) then - val(element)=-b1(x,y,z)/deltah2 - zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element)=-b1(x,y,z)/deltah2 - icol(element) = (ix-1)*ipoints*ipoints+(iy-1)*ipoints+(iz+1) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x,y+1,z) - if (iy == ipoints) then - val(element)=-b2(x,y,z)/deltah2 - zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element)=-b2(x,y,z)/deltah2 - icol(element) = (ix-1)*ipoints*ipoints+(iy)*ipoints+(iz) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x+1,y,z) - if (ix==ipoints) then - val(element)=-b3(x,y,z)/deltah2 - zt(k) = exp(-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element)=-b3(x,y,z)/deltah2 - icol(element) = (ix)*ipoints*ipoints+(iy-1)*ipoints+(iz) - irow(element) = glob_row - element = element+1 - endif - - end do - call psb_spins(element-1,irow,icol,val,a,desc_a,info) - if(info /= psb_success_) exit - call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),b,desc_a,info) - if(info /= psb_success_) exit - zt(:)=0.d0 - call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info) - if(info /= psb_success_) exit - end do - - tgen = psb_wtime()-t1 - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name) - goto 9999 - end if - - deallocate(val,irow,icol) - - call psb_barrier(ictxt) - t1 = psb_wtime() - call psb_cdasb(desc_a,info) - if (info == psb_success_) & - & call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_) - call psb_barrier(ictxt) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name) - goto 9999 - end if - call psb_geasb(b,desc_a,info) - call psb_geasb(xv,desc_a,info) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name) - goto 9999 - end if - tasb = psb_wtime()-t1 - call psb_barrier(ictxt) - ttot = psb_wtime() - t0 - - call psb_amx(ictxt,talc) - call psb_amx(ictxt,tgen) - call psb_amx(ictxt,tasb) - call psb_amx(ictxt,ttot) - if(iam == psb_root_) then - write(*,'("The matrix has been generated and assembled in ",a3," format.")')& - & a%get_fmt() - write(*,'("-allocation time : ",es12.5)') talc - write(*,'("-coeff. gen. time : ",es12.5)') tgen - write(*,'("-assembly time : ",es12.5)') tasb - write(*,'("-total time : ",es12.5)') ttot - + ! functions parametrizing the differential equation + ! + function b1(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: b1 + real(psb_spk_), intent(in) :: x,y,z + b1=1.e0/sqrt(3.e0) + end function b1 + function b2(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: b2 + real(psb_spk_), intent(in) :: x,y,z + b2=1.e0/sqrt(3.e0) + end function b2 + function b3(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: b3 + real(psb_spk_), intent(in) :: x,y,z + b3=1.e0/sqrt(3.e0) + end function b3 + function c(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: c + real(psb_spk_), intent(in) :: x,y,z + c=0.e0 + end function c + function a1(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: a1 + real(psb_spk_), intent(in) :: x,y,z + a1=1.e0/80 + end function a1 + function a2(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: a2 + real(psb_spk_), intent(in) :: x,y,z + a2=1.e0/80 + end function a2 + function a3(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: a3 + real(psb_spk_), intent(in) :: x,y,z + a3=1.e0/80 + end function a3 + function g(x,y,z) + use psb_base_mod, only : psb_spk_, sone + 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 - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error(ictxt) - return - end if - return - end subroutine create_matrix + end function g end program mld_sexample_1lev -! -! functions parametrizing the differential equation -! -function a1(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: a1 - real(psb_spk_) :: x,y,z - !a1=1.e0 - a1=0.e0 -end function a1 -function a2(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: a2 - real(psb_spk_) :: x,y,z - !a2=2.e1*y - a2=0.e0 -end function a2 -function a3(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: a3 - real(psb_spk_) :: x,y,z - !a3=1.e0 - a3=0.e0 -end function a3 -function a4(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: a4 - real(psb_spk_) :: x,y,z - !a4=1.e0 - a4=0.e0 -end function a4 -function b1(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: b1 - real(psb_spk_) :: x,y,z - b1=1.e0 -end function b1 -function b2(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: b2 - real(psb_spk_) :: x,y,z - b2=1.e0 -end function b2 -function b3(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: b3 - real(psb_spk_) :: x,y,z - b3=1.e0 -end function b3 diff --git a/examples/pdegen/mld_sexample_ml.f90 b/examples/pdegen/mld_sexample_ml.f90 index bbe07863..4a2e3d1c 100644 --- a/examples/pdegen/mld_sexample_ml.f90 +++ b/examples/pdegen/mld_sexample_ml.f90 @@ -48,29 +48,23 @@ ! ! The PDE is a general second order equation in 3d ! -! b1 dd(u) b2 dd(u) b3 dd(u) a1 d(u) a2 d(u) a3 d(u) -! - ------ - ------ - ------ - ----- - ------ - ------ + a4 u = 0 +! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u) +! - ------ - ------ - ------ + ----- + ------ + ------ + c u = f ! dxdx dydy dzdz dx dy dz ! -! with Dirichlet boundary conditions, on the unit cube 0<=x,y,z<=1. +! with Dirichlet boundary conditions +! u = g ! -! Example taken from: -! C.T.Kelley -! Iterative Methods for Linear and Nonlinear Equations -! SIAM 1995 +! on the unit cube 0<=x,y,z<=1. +! +! +! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation. ! ! In this sample program the index space of the discretized ! computational domain is first numbered sequentially in a standard way, ! then the corresponding vector is distributed according to a BLOCK ! data distribution. ! -! Boundary conditions are set in a very simple way, by adding -! equations of the form -! -! u(x,y) = exp(-x^2-y^2-z^2) -! -! Note that if a1=a2=a3=a4=0., the PDE is the well-known Laplace equation. -! program mld_sexample_ml use psb_base_mod use mld_prec_mod @@ -92,7 +86,7 @@ program mld_sexample_ml type(mld_sprec_type) :: P ! right-hand side, solution and residual vectors - real(psb_spk_), allocatable , save :: b(:), x(:), r(:) + type(psb_s_vect_type) :: x, b, r ! solver and preconditioner parameters real(psb_spk_) :: tol, err @@ -109,6 +103,7 @@ program mld_sexample_ml integer :: idim, ierr, ircode real(psb_dpk_) :: t1, t2, tprec real(psb_spk_) :: resmx, resmxp + character(len=5) :: afmt='CSR' character(len=20) :: name ! initialize the parallel environment @@ -142,7 +137,8 @@ program mld_sexample_ml call psb_barrier(ictxt) t1 = psb_wtime() - call create_matrix(idim,a,b,x,desc_a,ictxt,info) + call psb_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,& + & a1,a2,a3,b1,b2,b3,c,g,info) call psb_barrier(ictxt) t2 = psb_wtime() - t1 if(info /= psb_success_) then @@ -209,7 +205,7 @@ program mld_sexample_ml ! set the solver parameters and the initial guess call psb_geall(x,desc_A,info) - x(:) =0.0 + call x%set(szero) call psb_geasb(x,desc_A,info) ! solve Ax=b with preconditioned BiCGSTAB @@ -223,12 +219,12 @@ program mld_sexample_ml call psb_amx(ictxt,t2) call psb_geall(r,desc_A,info) - r(:) =0.0 + call r%set(szero) call psb_geasb(r,desc_A,info) call psb_geaxpby(sone,b,szero,r,desc_A,info) call psb_spmm(-sone,A,x,sone,r,desc_A,info) - call psb_genrm2s(resmx,r,desc_A,info) - call psb_geamaxs(resmxp,r,desc_A,info) + resmx = psb_genrm2(r,desc_A,info) + resmxp = psb_geamax(r,desc_A,info) amatsize = a%sizeof() descsize = desc_a%sizeof() @@ -298,332 +294,61 @@ contains call psb_bcast(ictxt,tol) end subroutine get_parms - - ! - ! subroutine to allocate and fill in the coefficient matrix and - ! the rhs ! - subroutine create_matrix(idim,a,b,xv,desc_a,ictxt,info) - ! - ! Discretize the partial diferential equation - ! - ! b1 dd(u) b2 dd(u) b3 dd(u) a1 d(u) a2 d(u) a3 d(u) - ! - ------ - ------ - ------ - ----- - ------ - ------ + a4 u = 0 - ! dxdx dydy dzdz dx dy dz - ! - ! with Dirichlet boundary conditions, on the unit cube 0<=x,y,z<=1. - ! - ! Boundary conditions are set in a very simple way, by adding - ! equations of the form - ! - ! u(x,y) = exp(-x^2-y^2-z^2) - ! - ! Note that if a1=a2=a3=a4=0., the PDE is the well-known Laplace equation. - ! - use psb_base_mod - implicit none - integer :: idim - integer, parameter :: nb=20 - real(psb_spk_), allocatable :: b(:),xv(:) - type(psb_desc_type) :: desc_a - integer :: ictxt, info - character :: afmt*5 - type(psb_sspmat_type) :: a - real(psb_spk_) :: zt(nb),x,y,z - integer :: m,n,nnz,glob_row,nlr,i,ii,ib,k - integer :: ix,iy,iz,ia,indx_owner, ipoints - integer :: np, iam, nr, nt - integer :: element - integer, allocatable :: irow(:),icol(:),myidx(:) - real(psb_spk_), allocatable :: val(:) - ! deltah dimension of each grid cell - ! deltat discretization time - real(psb_spk_) :: deltah, deltah2 - real(psb_spk_),parameter :: rhs=0.0,one=1.0,zero=0.0 - real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen - real(psb_spk_) :: a1, a2, a3, a4, b1, b2, b3 - external :: a1, a2, a3, a4, b1, b2, b3 - integer :: err_act - - character(len=20) :: name, ch_err - - info = psb_success_ - name = 'create_matrix' - call psb_erractionsave(err_act) - - call psb_info(ictxt, iam, np) - - deltah = 1.d0/(idim-1) - deltah2 = deltah*deltah - - ! initialize array descriptor and sparse matrix storage. provide an - ! estimate of the number of non zeroes - - ipoints=idim-2 - m = ipoints*ipoints*ipoints - n = m - nnz = ((n*9)/(np)) - if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n - - ! - ! Using a simple BLOCK distribution. - ! - nt = (m+np-1)/np - nr = max(0,min(nt,m-(iam*nt))) - - nt = nr - call psb_sum(ictxt,nt) - if (nt /= m) write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m - call psb_barrier(ictxt) - t0 = psb_wtime() - call psb_cdall(ictxt,desc_a,info,nl=nr) - if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz) - ! define rhs from boundary conditions; also build initial guess - if (info == psb_success_) call psb_geall(b,desc_a,info) - if (info == psb_success_) call psb_geall(xv,desc_a,info) - nlr = psb_cd_get_local_rows(desc_a) - call psb_barrier(ictxt) - talc = psb_wtime()-t0 - - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='allocation rout.' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - ! we build an auxiliary matrix consisting of one row at a - ! time; just a small matrix. might be extended to generate - ! a bunch of rows per call. - ! - allocate(val(20*nb),irow(20*nb),& - &icol(20*nb),myidx(nlr),stat=info) - if (info /= psb_success_ ) then - info=psb_err_alloc_dealloc_ - call psb_errpush(info,name) - goto 9999 - endif - - do i=1,nlr - myidx(i) = i - end do - - - call psb_loc_to_glob(myidx,desc_a,info) - - ! loop over rows belonging to current process in a block - ! distribution. - - call psb_barrier(ictxt) - t1 = psb_wtime() - do ii=1, nlr,nb - ib = min(nb,nlr-ii+1) - element = 1 - do k=1,ib - i=ii+k-1 - ! local matrix pointer - glob_row=myidx(i) - ! compute gridpoint coordinates - if (mod(glob_row,ipoints*ipoints) == 0) then - ix = glob_row/(ipoints*ipoints) - else - ix = glob_row/(ipoints*ipoints)+1 - endif - if (mod((glob_row-(ix-1)*ipoints*ipoints),ipoints) == 0) then - iy = (glob_row-(ix-1)*ipoints*ipoints)/ipoints - else - iy = (glob_row-(ix-1)*ipoints*ipoints)/ipoints+1 - endif - iz = glob_row-(ix-1)*ipoints*ipoints-(iy-1)*ipoints - ! x, y, x coordinates - x=ix*deltah - y=iy*deltah - z=iz*deltah - - ! check on boundary points - zt(k) = 0.d0 - ! internal point: build discretization - ! - ! term depending on (x-1,y,z) - ! - if (ix == 1) then - val(element) = -b1(x,y,z)/deltah2-a1(x,y,z)/deltah - zt(k) = exp(-x**2-y**2-z**2)*(-val(element)) - else - val(element) = -b1(x,y,z)/deltah2-a1(x,y,z)/deltah - icol(element) = (ix-2)*ipoints*ipoints+(iy-1)*ipoints+(iz) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x,y-1,z) - if (iy == 1) then - val(element) = -b2(x,y,z)/deltah2-a2(x,y,z)/deltah - zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element) = -b2(x,y,z)/deltah2-a2(x,y,z)/deltah - icol(element) = (ix-1)*ipoints*ipoints+(iy-2)*ipoints+(iz) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x,y,z-1) - if (iz == 1) then - val(element)=-b3(x,y,z)/deltah2-a3(x,y,z)/deltah - zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element)=-b3(x,y,z)/deltah2-a3(x,y,z)/deltah - icol(element) = (ix-1)*ipoints*ipoints+(iy-1)*ipoints+(iz-1) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x,y,z) - val(element)=(2*b1(x,y,z) + 2*b2(x,y,z) + 2*b3(x,y,z))/deltah2& - & + (a1(x,y,z) + a2(x,y,z) + a3(x,y,z)+ a4(x,y,z))/deltah - icol(element) = (ix-1)*ipoints*ipoints+(iy-1)*ipoints+(iz) - irow(element) = glob_row - element = element+1 - ! term depending on (x,y,z+1) - if (iz == ipoints) then - val(element)=-b1(x,y,z)/deltah2 - zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element)=-b1(x,y,z)/deltah2 - icol(element) = (ix-1)*ipoints*ipoints+(iy-1)*ipoints+(iz+1) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x,y+1,z) - if (iy == ipoints) then - val(element)=-b2(x,y,z)/deltah2 - zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element)=-b2(x,y,z)/deltah2 - icol(element) = (ix-1)*ipoints*ipoints+(iy)*ipoints+(iz) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x+1,y,z) - if (ix==ipoints) then - val(element)=-b3(x,y,z)/deltah2 - zt(k) = exp(-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element)=-b3(x,y,z)/deltah2 - icol(element) = (ix)*ipoints*ipoints+(iy-1)*ipoints+(iz) - irow(element) = glob_row - element = element+1 - endif - - end do - call psb_spins(element-1,irow,icol,val,a,desc_a,info) - if(info /= psb_success_) exit - call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),b,desc_a,info) - if(info /= psb_success_) exit - zt(:)=0.d0 - call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info) - if(info /= psb_success_) exit - end do - - tgen = psb_wtime()-t1 - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name) - goto 9999 - end if - - deallocate(val,irow,icol) - - call psb_barrier(ictxt) - t1 = psb_wtime() - call psb_cdasb(desc_a,info) - if (info == psb_success_) & - & call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_) - call psb_barrier(ictxt) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name) - goto 9999 - end if - call psb_geasb(b,desc_a,info) - call psb_geasb(xv,desc_a,info) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name) - goto 9999 - end if - tasb = psb_wtime()-t1 - call psb_barrier(ictxt) - ttot = psb_wtime() - t0 - - call psb_amx(ictxt,talc) - call psb_amx(ictxt,tgen) - call psb_amx(ictxt,tasb) - call psb_amx(ictxt,ttot) - if(iam == psb_root_) then - write(*,'("The matrix has been generated and assembled in ",a3," format.")')& - & a%get_fmt() - write(*,'("-allocation time : ",es12.5)') talc - write(*,'("-coeff. gen. time : ",es12.5)') tgen - write(*,'("-assembly time : ",es12.5)') tasb - write(*,'("-total time : ",es12.5)') ttot - + ! functions parametrizing the differential equation + ! + function b1(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: b1 + real(psb_spk_), intent(in) :: x,y,z + b1=1.e0/sqrt(3.e0) + end function b1 + function b2(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: b2 + real(psb_spk_), intent(in) :: x,y,z + b2=1.e0/sqrt(3.e0) + end function b2 + function b3(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: b3 + real(psb_spk_), intent(in) :: x,y,z + b3=1.e0/sqrt(3.e0) + end function b3 + function c(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: c + real(psb_spk_), intent(in) :: x,y,z + c=0.e0 + end function c + function a1(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: a1 + real(psb_spk_), intent(in) :: x,y,z + a1=1.e0/80 + end function a1 + function a2(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: a2 + real(psb_spk_), intent(in) :: x,y,z + a2=1.e0/80 + end function a2 + function a3(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: a3 + real(psb_spk_), intent(in) :: x,y,z + a3=1.e0/80 + end function a3 + function g(x,y,z) + use psb_base_mod, only : psb_spk_, sone + 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 - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error(ictxt) - return - end if - return - end subroutine create_matrix + end function g + end program mld_sexample_ml -! -! functions parametrizing the differential equation -! -function a1(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: a1 - real(psb_spk_) :: x,y,z -! a1=1.e0 - a1=0.e0 -end function a1 -function a2(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: a2 - real(psb_spk_) :: x,y,z -! a2=2.e1*y - a2=0.e0 -end function a2 -function a3(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: a3 - real(psb_spk_) :: x,y,z -! a3=1.e0 - a3=0.e0 -end function a3 -function a4(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: a4 - real(psb_spk_) :: x,y,z -! a4=1.e0 - a4=0.e0 -end function a4 -function b1(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: b1 - real(psb_spk_) :: x,y,z - b1=1.e0 -end function b1 -function b2(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: b2 - real(psb_spk_) :: x,y,z - b2=1.e0 -end function b2 -function b3(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: b3 - real(psb_spk_) :: x,y,z - b3=1.e0 -end function b3 diff --git a/mlprec/mld_base_prec_type.F90 b/mlprec/mld_base_prec_type.F90 index c49ba04f..a920c19f 100644 --- a/mlprec/mld_base_prec_type.F90 +++ b/mlprec/mld_base_prec_type.F90 @@ -195,11 +195,14 @@ module mld_base_prec_type ! ! Legal values for entry: mld_ilu_scale_ ! - integer, parameter :: mld_ilu_scale_none_ = 0 - integer, parameter :: mld_ilu_scale_maxval_ = 1 - integer, parameter :: mld_ilu_scale_diag_ = 2 + integer, parameter :: mld_ilu_scale_none_ = 0 + integer, parameter :: mld_ilu_scale_maxval_ = 1 + integer, parameter :: mld_ilu_scale_diag_ = 2 + integer, parameter :: mld_ilu_scale_arwsum_ = 3 + integer, parameter :: mld_ilu_scale_aclsum_ = 4 + integer, parameter :: mld_ilu_scale_dabsum_ = 5 ! For the time being enable only maxval scale - integer, parameter :: mld_max_ilu_scale_ = 1 + integer, parameter :: mld_max_ilu_scale_ = 1 ! ! Legal values for entry: mld_ml_type_ !