From 69cd86f079aa176231bb16df2d279c1da910af96 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Mon, 16 Oct 2017 16:54:39 +0100 Subject: [PATCH] New example programs --- examples/fileread/mld_dexample_1lev.f90 | 15 +- examples/pdegen/mld_dexample_1lev.f90 | 12 +- examples/pdegen/mld_dexample_ml.f90 | 5 +- examples/pdegen/mld_dpde_mod.f90 | 348 +++++++++++++++++++++++ examples/pdegen/mld_sexample_1lev.f90 | 12 +- examples/pdegen/mld_sexample_ml.f90 | 5 +- examples/pdegen/mld_spde_mod.f90 | 348 +++++++++++++++++++++++ tests/fileread/runs/dfs.inp | 26 +- tests/pdegen/mld_d_pde2d.f90 | 310 ++++++++++++++++++++- tests/pdegen/mld_d_pde3d.f90 | 351 +++++++++++++++++++++++- tests/pdegen/mld_s_pde2d.f90 | 310 ++++++++++++++++++++- tests/pdegen/mld_s_pde3d.f90 | 351 +++++++++++++++++++++++- 12 files changed, 2052 insertions(+), 41 deletions(-) diff --git a/examples/fileread/mld_dexample_1lev.f90 b/examples/fileread/mld_dexample_1lev.f90 index 8d290e8b..f5cc89ea 100644 --- a/examples/fileread/mld_dexample_1lev.f90 +++ b/examples/fileread/mld_dexample_1lev.f90 @@ -193,12 +193,15 @@ program mld_dexample_1lev ! set RAS - call P%init('AS',info) - - ! set number of overlaps - - call P%set('SUB_OVR',2,info) - +!!$ call P%init('AS',info) +!!$ +!!$ ! set number of overlaps +!!$ +!!$ call P%set('SUB_OVR',2,info) +!!$ + call P%init('FBGS',info) + !call prec%set('smoother_sweeps', p_choice%outer_sweeps, info) + ! build the preconditioner t1 = psb_wtime() diff --git a/examples/pdegen/mld_dexample_1lev.f90 b/examples/pdegen/mld_dexample_1lev.f90 index 12a4b2eb..fda9e6e8 100644 --- a/examples/pdegen/mld_dexample_1lev.f90 +++ b/examples/pdegen/mld_dexample_1lev.f90 @@ -93,7 +93,7 @@ program mld_dexample_1lev real(psb_dpk_) :: resmx, resmxp real(psb_dpk_) :: t1, t2, tprec character(len=5) :: afmt='CSR' - character(len=20) :: name + character(len=20) :: name, kmethod ! initialize the parallel environment call psb_init(ictxt) @@ -125,7 +125,7 @@ program mld_dexample_1lev call psb_barrier(ictxt) t1 = psb_wtime() - call psb_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,& + call mld_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 @@ -166,18 +166,19 @@ program mld_dexample_1lev call x%zero() call psb_geasb(x,desc_A,info) - ! solve Ax=b with preconditioned BiCGSTAB + ! solve Ax=b with preconditioned Krylov method: BiCGSTAB + kmethod = 'BiCGSTAB' call psb_barrier(ictxt) t1 = psb_wtime() - call psb_krylov('BICGSTAB',A,P,b,x,tol,desc_A,info,itmax,iter,err,itrace=1,istop=2) + call psb_krylov(kmethod,A,P,b,x,tol,desc_A,info,itmax,iter,err,itrace=1,istop=2) t2 = psb_wtime() - t1 call psb_amx(ictxt,t2) call psb_geall(r,desc_A,info) - call r%set(dzero) + call r%zero() 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) @@ -197,6 +198,7 @@ program mld_dexample_1lev write(*,'(" ")') write(*,'("Matrix from PDE example")') write(*,'("Computed solution on ",i8," processors")')np + write(*,'("Krylov method : ",a)') kmethod write(*,'("Iterations to convergence : ",i6)')iter write(*,'("Error estimate on exit : ",es12.5)')err write(*,'("Time to build prec. : ",es12.5)')tprec diff --git a/examples/pdegen/mld_dexample_ml.f90 b/examples/pdegen/mld_dexample_ml.f90 index 3640ad15..7c2614c7 100644 --- a/examples/pdegen/mld_dexample_ml.f90 +++ b/examples/pdegen/mld_dexample_ml.f90 @@ -146,7 +146,7 @@ program mld_dexample_ml call psb_barrier(ictxt) t1 = psb_wtime() - call psb_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,& + call mld_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 @@ -219,7 +219,7 @@ program mld_dexample_ml call x%zero() call psb_geasb(x,desc_A,info) - ! solve Ax=b with preconditioned CG + ! solve Ax=b with preconditioned Krylov method call psb_barrier(ictxt) t1 = psb_wtime() @@ -250,6 +250,7 @@ program mld_dexample_ml write(*,'(" ")') write(*,'("Matrix from PDE example")') write(*,'("Computed solution on ",i8," processors")')np + write(*,'("Krylov method : ",a)') kmethod write(*,'("Iterations to convergence : ",i6)')iter write(*,'("Error estimate on exit : ",es12.5)')err write(*,'("Time to build prec. : ",es12.5)')tprec diff --git a/examples/pdegen/mld_dpde_mod.f90 b/examples/pdegen/mld_dpde_mod.f90 index 9c1a0179..b9ae4911 100644 --- a/examples/pdegen/mld_dpde_mod.f90 +++ b/examples/pdegen/mld_dpde_mod.f90 @@ -35,7 +35,355 @@ ! POSSIBILITY OF SUCH DAMAGE. ! module mld_d_pde_mod + use psb_base_mod, only : psb_dpk_, psb_ipk_, psb_desc_type,& + & psb_dspmat_type, psb_d_vect_type, dzero,& + & psb_d_base_sparse_mat, psb_d_base_vect_type, psb_i_base_vect_type + + interface + function d_func_3d(x,y,z) result(val) + import :: psb_dpk_ + real(psb_dpk_), intent(in) :: x,y,z + real(psb_dpk_) :: val + end function d_func_3d + end interface + + interface mld_gen_pde3d + module procedure mld_d_gen_pde3d + end interface mld_gen_pde3d + + contains + + function d_null_func_3d(x,y,z) result(val) + + real(psb_dpk_), intent(in) :: x,y,z + real(psb_dpk_) :: val + + val = dzero + + end function d_null_func_3d + + ! + ! subroutine to allocate and fill in the coefficient matrix and + ! the rhs. + ! + subroutine mld_d_gen_pde3d(ictxt,idim,a,bv,xv,desc_a,afmt,& + & a1,a2,a3,b1,b2,b3,c,g,info,f,amold,vmold,imold,nrl,iv) + use psb_base_mod + ! + ! Discretizes the partial differential equation + ! + ! 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 + ! u = g + ! + ! on the unit cube 0<=x,y,z<=1. + ! + ! + ! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation. + ! + implicit none + procedure(d_func_3d) :: b1,b2,b3,c,a1,a2,a3,g + integer(psb_ipk_) :: idim + type(psb_dspmat_type) :: a + type(psb_d_vect_type) :: xv,bv + type(psb_desc_type) :: desc_a + integer(psb_ipk_) :: ictxt, info + character(len=*) :: afmt + procedure(d_func_3d), optional :: f + class(psb_d_base_sparse_mat), optional :: amold + class(psb_d_base_vect_type), optional :: vmold + class(psb_i_base_vect_type), optional :: imold + integer(psb_ipk_), optional :: nrl,iv(:) + + ! Local variables. + + integer(psb_ipk_), parameter :: nb=20 + type(psb_d_csc_sparse_mat) :: acsc + type(psb_d_coo_sparse_mat) :: acoo + type(psb_d_csr_sparse_mat) :: acsr + real(psb_dpk_) :: zt(nb),x,y,z + integer(psb_ipk_) :: m,n,nnz,glob_row,nlr,i,ii,ib,k + integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner + integer(psb_ipk_) :: np, iam, nr, nt + integer(psb_ipk_) :: icoeff + integer(psb_ipk_), allocatable :: irow(:),icol(:),myidx(:) + real(psb_dpk_), allocatable :: val(:) + ! deltah dimension of each grid cell + ! deltat discretization time + real(psb_dpk_) :: deltah, sqdeltah, deltah2 + real(psb_dpk_), parameter :: rhs=dzero,one=done,zero=dzero + real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb + integer(psb_ipk_) :: err_act + procedure(d_func_3d), pointer :: f_ + character(len=20) :: name, ch_err,tmpfmt + + info = psb_success_ + name = 'create_matrix' + call psb_erractionsave(err_act) + + call psb_info(ictxt, iam, np) + + + if (present(f)) then + f_ => f + else + f_ => d_null_func_3d + end if + + deltah = 1.d0/(idim+2) + sqdeltah = deltah*deltah + deltah2 = 2.d0* deltah + + ! initialize array descriptor and sparse matrix storage. provide an + ! estimate of the number of non zeroes + + m = idim*idim*idim + n = m + nnz = ((n*9)/(np)) + if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n + + if (.not.present(iv)) then + if (present(nrl)) then + nr = nrl + else + ! + ! Using a simple BLOCK distribution. + ! + nt = (m+np-1)/np + nr = max(0,min(nt,m-(iam*nt))) + end if + + nt = nr + call psb_sum(ictxt,nt) + if (nt /= m) then + write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m + info = -1 + call psb_barrier(ictxt) + call psb_abort(ictxt) + return + end if + else + if (size(iv) /= m) then + write(psb_err_unit,*) iam, 'Initialization error IV',size(iv),m + info = -1 + call psb_barrier(ictxt) + call psb_abort(ictxt) + return + end if + + end if + call psb_barrier(ictxt) + t0 = psb_wtime() + if (present(iv)) then + call psb_cdall(ictxt,desc_a,info,vg=iv) + else + call psb_cdall(ictxt,desc_a,info,nl=nr) + end if + + 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(xv,desc_a,info) + if (info == psb_success_) call psb_geall(bv,desc_a,info) + + 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),stat=info) + if (info /= psb_success_ ) then + info=psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + endif + + myidx = desc_a%get_global_indices() + nlr = size(myidx) + + ! 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) + icoeff = 1 + do k=1,ib + i=ii+k-1 + ! local matrix pointer + glob_row=myidx(i) + ! compute gridpoint coordinates + if (mod(glob_row,(idim*idim)) == 0) then + ix = glob_row/(idim*idim) + else + ix = glob_row/(idim*idim)+1 + endif + if (mod((glob_row-(ix-1)*idim*idim),idim) == 0) then + iy = (glob_row-(ix-1)*idim*idim)/idim + else + iy = (glob_row-(ix-1)*idim*idim)/idim+1 + endif + iz = glob_row-(ix-1)*idim*idim-(iy-1)*idim + ! x, y, x coordinates + x = (ix-1)*deltah + y = (iy-1)*deltah + z = (iz-1)*deltah + zt(k) = f_(x,y,z) + ! internal point: build discretization + ! + ! term depending on (x-1,y,z) + ! + val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2 + if (ix == 1) then + zt(k) = g(dzero,y,z)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-2)*idim*idim+(iy-1)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y-1,z) + val(icoeff) = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2 + if (iy == 1) then + zt(k) = g(x,dzero,z)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim*idim+(iy-2)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y,z-1) + val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2 + if (iz == 1) then + zt(k) = g(x,y,dzero)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz-1) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + ! term depending on (x,y,z) + val(icoeff)=2.d0*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah & + & + c(x,y,z) + icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + ! term depending on (x,y,z+1) + val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2 + if (iz == idim) then + zt(k) = g(x,y,done)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz+1) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y+1,z) + val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2 + if (iy == idim) then + zt(k) = g(x,done,z)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim*idim+(iy)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x+1,y,z) + val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2 + if (ix==idim) then + zt(k) = g(done,y,z)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix)*idim*idim+(iy-1)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + end do + call psb_spins(icoeff-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),bv,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_ + ch_err='insert rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + deallocate(val,irow,icol) + + call psb_barrier(ictxt) + t1 = psb_wtime() + call psb_cdasb(desc_a,info,mold=imold) + tcdasb = psb_wtime()-t1 + call psb_barrier(ictxt) + t1 = psb_wtime() + if (info == psb_success_) then + if (present(amold)) then + call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,mold=amold) + else + call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) + end if + end if + call psb_barrier(ictxt) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='asb rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + if (info == psb_success_) call psb_geasb(xv,desc_a,info,mold=vmold) + if (info == psb_success_) call psb_geasb(bv,desc_a,info,mold=vmold) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='asb rout.' + call psb_errpush(info,name,a_err=ch_err) + 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 + tmpfmt = a%get_fmt() + write(psb_out_unit,'("The matrix has been generated and assembled in ",a3," format.")')& + & tmpfmt + write(psb_out_unit,'("-allocation time : ",es12.5)') talc + write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen + write(psb_out_unit,'("-desc asbly time : ",es12.5)') tcdasb + write(psb_out_unit,'("- mat asbly time : ",es12.5)') tasb + write(psb_out_unit,'("-total time : ",es12.5)') ttot + + end if + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + end subroutine mld_d_gen_pde3d + + ! ! functions parametrizing the differential equation ! diff --git a/examples/pdegen/mld_sexample_1lev.f90 b/examples/pdegen/mld_sexample_1lev.f90 index fa8eb2f1..3f7e5547 100644 --- a/examples/pdegen/mld_sexample_1lev.f90 +++ b/examples/pdegen/mld_sexample_1lev.f90 @@ -93,7 +93,7 @@ program mld_sexample_1lev real(psb_spk_) :: resmx, resmxp real(psb_dpk_) :: t1, t2, tprec character(len=5) :: afmt='CSR' - character(len=20) :: name + character(len=20) :: name, kmethod ! initialize the parallel environment call psb_init(ictxt) @@ -125,7 +125,7 @@ program mld_sexample_1lev call psb_barrier(ictxt) t1 = psb_wtime() - call psb_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,& + call mld_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 @@ -166,18 +166,19 @@ program mld_sexample_1lev call x%zero() call psb_geasb(x,desc_A,info) - ! solve Ax=b with preconditioned BiCGSTAB + ! solve Ax=b with preconditioned Krylov method: BiCGSTAB + kmethod = 'BiCGSTAB' call psb_barrier(ictxt) t1 = psb_wtime() - call psb_krylov('BICGSTAB',A,P,b,x,tol,desc_A,info,itmax,iter,err,itrace=1,istop=2) + call psb_krylov(kmethod,A,P,b,x,tol,desc_A,info,itmax,iter,err,itrace=1,istop=2) t2 = psb_wtime() - t1 call psb_amx(ictxt,t2) call psb_geall(r,desc_A,info) - call r%set(szero) + call r%zero() 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) @@ -197,6 +198,7 @@ program mld_sexample_1lev write(*,'(" ")') write(*,'("Matrix from PDE example")') write(*,'("Computed solution on ",i8," processors")')np + write(*,'("Krylov method : ",a)') kmethod write(*,'("Iterations to convergence : ",i6)')iter write(*,'("Error estimate on exit : ",es12.5)')err write(*,'("Time to build prec. : ",es12.5)')tprec diff --git a/examples/pdegen/mld_sexample_ml.f90 b/examples/pdegen/mld_sexample_ml.f90 index ae05dd2f..d1cb4ff4 100644 --- a/examples/pdegen/mld_sexample_ml.f90 +++ b/examples/pdegen/mld_sexample_ml.f90 @@ -146,7 +146,7 @@ program mld_sexample_ml call psb_barrier(ictxt) t1 = psb_wtime() - call psb_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,& + call mld_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 @@ -219,7 +219,7 @@ program mld_sexample_ml call x%zero() call psb_geasb(x,desc_A,info) - ! solve Ax=b with preconditioned CG + ! solve Ax=b with preconditioned Krylov method call psb_barrier(ictxt) t1 = psb_wtime() @@ -250,6 +250,7 @@ program mld_sexample_ml write(*,'(" ")') write(*,'("Matrix from PDE example")') write(*,'("Computed solution on ",i8," processors")')np + write(*,'("Krylov method : ",a)') kmethod write(*,'("Iterations to convergence : ",i6)')iter write(*,'("Error estimate on exit : ",es12.5)')err write(*,'("Time to build prec. : ",es12.5)')tprec diff --git a/examples/pdegen/mld_spde_mod.f90 b/examples/pdegen/mld_spde_mod.f90 index a4899828..94e46e21 100644 --- a/examples/pdegen/mld_spde_mod.f90 +++ b/examples/pdegen/mld_spde_mod.f90 @@ -35,7 +35,355 @@ ! POSSIBILITY OF SUCH DAMAGE. ! module mld_s_pde_mod + use psb_base_mod, only : psb_spk_, psb_ipk_, psb_desc_type,& + & psb_sspmat_type, psb_s_vect_type, szero,& + & psb_s_base_sparse_mat, psb_s_base_vect_type, psb_i_base_vect_type + + interface + function s_func_3d(x,y,z) result(val) + import :: psb_spk_ + real(psb_spk_), intent(in) :: x,y,z + real(psb_spk_) :: val + end function s_func_3d + end interface + + interface mld_gen_pde3d + module procedure mld_s_gen_pde3d + end interface mld_gen_pde3d + + contains + + function s_null_func_3d(x,y,z) result(val) + + real(psb_spk_), intent(in) :: x,y,z + real(psb_spk_) :: val + + val = szero + + end function s_null_func_3d + + ! + ! subroutine to allocate and fill in the coefficient matrix and + ! the rhs. + ! + subroutine mld_s_gen_pde3d(ictxt,idim,a,bv,xv,desc_a,afmt,& + & a1,a2,a3,b1,b2,b3,c,g,info,f,amold,vmold,imold,nrl,iv) + use psb_base_mod + ! + ! Discretizes the partial differential equation + ! + ! 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 + ! u = g + ! + ! on the unit cube 0<=x,y,z<=1. + ! + ! + ! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation. + ! + implicit none + procedure(s_func_3d) :: b1,b2,b3,c,a1,a2,a3,g + integer(psb_ipk_) :: idim + type(psb_sspmat_type) :: a + type(psb_s_vect_type) :: xv,bv + type(psb_desc_type) :: desc_a + integer(psb_ipk_) :: ictxt, info + character(len=*) :: afmt + procedure(s_func_3d), optional :: f + class(psb_s_base_sparse_mat), optional :: amold + class(psb_s_base_vect_type), optional :: vmold + class(psb_i_base_vect_type), optional :: imold + integer(psb_ipk_), optional :: nrl,iv(:) + + ! Local variables. + + integer(psb_ipk_), parameter :: nb=20 + type(psb_s_csc_sparse_mat) :: acsc + type(psb_s_coo_sparse_mat) :: acoo + type(psb_s_csr_sparse_mat) :: acsr + real(psb_spk_) :: zt(nb),x,y,z + integer(psb_ipk_) :: m,n,nnz,glob_row,nlr,i,ii,ib,k + integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner + integer(psb_ipk_) :: np, iam, nr, nt + integer(psb_ipk_) :: icoeff + integer(psb_ipk_), allocatable :: irow(:),icol(:),myidx(:) + real(psb_spk_), allocatable :: val(:) + ! deltah dimension of each grid cell + ! deltat discretization time + real(psb_spk_) :: deltah, sqdeltah, deltah2 + real(psb_spk_), parameter :: rhs=szero,one=sone,zero=szero + real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb + integer(psb_ipk_) :: err_act + procedure(s_func_3d), pointer :: f_ + character(len=20) :: name, ch_err,tmpfmt + + info = psb_success_ + name = 'create_matrix' + call psb_erractionsave(err_act) + + call psb_info(ictxt, iam, np) + + + if (present(f)) then + f_ => f + else + f_ => s_null_func_3d + end if + + deltah = 1.d0/(idim+2) + sqdeltah = deltah*deltah + deltah2 = 2.d0* deltah + + ! initialize array descriptor and sparse matrix storage. provide an + ! estimate of the number of non zeroes + + m = idim*idim*idim + n = m + nnz = ((n*9)/(np)) + if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n + + if (.not.present(iv)) then + if (present(nrl)) then + nr = nrl + else + ! + ! Using a simple BLOCK distribution. + ! + nt = (m+np-1)/np + nr = max(0,min(nt,m-(iam*nt))) + end if + + nt = nr + call psb_sum(ictxt,nt) + if (nt /= m) then + write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m + info = -1 + call psb_barrier(ictxt) + call psb_abort(ictxt) + return + end if + else + if (size(iv) /= m) then + write(psb_err_unit,*) iam, 'Initialization error IV',size(iv),m + info = -1 + call psb_barrier(ictxt) + call psb_abort(ictxt) + return + end if + + end if + call psb_barrier(ictxt) + t0 = psb_wtime() + if (present(iv)) then + call psb_cdall(ictxt,desc_a,info,vg=iv) + else + call psb_cdall(ictxt,desc_a,info,nl=nr) + end if + + 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(xv,desc_a,info) + if (info == psb_success_) call psb_geall(bv,desc_a,info) + + 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),stat=info) + if (info /= psb_success_ ) then + info=psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + endif + + myidx = desc_a%get_global_indices() + nlr = size(myidx) + + ! 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) + icoeff = 1 + do k=1,ib + i=ii+k-1 + ! local matrix pointer + glob_row=myidx(i) + ! compute gridpoint coordinates + if (mod(glob_row,(idim*idim)) == 0) then + ix = glob_row/(idim*idim) + else + ix = glob_row/(idim*idim)+1 + endif + if (mod((glob_row-(ix-1)*idim*idim),idim) == 0) then + iy = (glob_row-(ix-1)*idim*idim)/idim + else + iy = (glob_row-(ix-1)*idim*idim)/idim+1 + endif + iz = glob_row-(ix-1)*idim*idim-(iy-1)*idim + ! x, y, x coordinates + x = (ix-1)*deltah + y = (iy-1)*deltah + z = (iz-1)*deltah + zt(k) = f_(x,y,z) + ! internal point: build discretization + ! + ! term depending on (x-1,y,z) + ! + val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2 + if (ix == 1) then + zt(k) = g(szero,y,z)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-2)*idim*idim+(iy-1)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y-1,z) + val(icoeff) = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2 + if (iy == 1) then + zt(k) = g(x,szero,z)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim*idim+(iy-2)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y,z-1) + val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2 + if (iz == 1) then + zt(k) = g(x,y,szero)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz-1) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + ! term depending on (x,y,z) + val(icoeff)=2.d0*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah & + & + c(x,y,z) + icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + ! term depending on (x,y,z+1) + val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2 + if (iz == idim) then + zt(k) = g(x,y,sone)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz+1) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y+1,z) + val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2 + if (iy == idim) then + zt(k) = g(x,sone,z)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim*idim+(iy)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x+1,y,z) + val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2 + if (ix==idim) then + zt(k) = g(sone,y,z)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix)*idim*idim+(iy-1)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + end do + call psb_spins(icoeff-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),bv,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_ + ch_err='insert rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + deallocate(val,irow,icol) + + call psb_barrier(ictxt) + t1 = psb_wtime() + call psb_cdasb(desc_a,info,mold=imold) + tcdasb = psb_wtime()-t1 + call psb_barrier(ictxt) + t1 = psb_wtime() + if (info == psb_success_) then + if (present(amold)) then + call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,mold=amold) + else + call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) + end if + end if + call psb_barrier(ictxt) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='asb rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + if (info == psb_success_) call psb_geasb(xv,desc_a,info,mold=vmold) + if (info == psb_success_) call psb_geasb(bv,desc_a,info,mold=vmold) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='asb rout.' + call psb_errpush(info,name,a_err=ch_err) + 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 + tmpfmt = a%get_fmt() + write(psb_out_unit,'("The matrix has been generated and assembled in ",a3," format.")')& + & tmpfmt + write(psb_out_unit,'("-allocation time : ",es12.5)') talc + write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen + write(psb_out_unit,'("-desc asbly time : ",es12.5)') tcdasb + write(psb_out_unit,'("- mat asbly time : ",es12.5)') tasb + write(psb_out_unit,'("-total time : ",es12.5)') ttot + + end if + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + end subroutine mld_s_gen_pde3d + + ! ! functions parametrizing the differential equation ! diff --git a/tests/fileread/runs/dfs.inp b/tests/fileread/runs/dfs.inp index 23853f9d..56938d6a 100644 --- a/tests/fileread/runs/dfs.inp +++ b/tests/fileread/runs/dfs.inp @@ -1,11 +1,11 @@ -mld_mat.mtx ! This matrix (and others) from: http://math.nist.gov/MatrixMarket/ or -mld_rhs.mtx ! rhs | http://www.cise.ufl.edu/research/sparse/matrices/index.html +poisson_200x200.mtx ! This matrix (and others) from: http://math.nist.gov/MatrixMarket/ or +NONE ! rhs | http://www.cise.ufl.edu/research/sparse/matrices/index.html NONE ! Initial guess -NONE !mld_sol.mtx ! Reference solution +NONE !mld_sol.mtx ! Reference solution MM ! File format: MatrixMarket or Harwell-Boeing CSR ! Storage format: CSR COO JAD GRAPH ! PART (partition method): BLOCK GRAPH -BICGSTAB ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES +CG ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES 2 ! ISTOPC 00500 ! ITMAX 1 ! ITRACE @@ -14,35 +14,35 @@ BICGSTAB ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FC ML-VCYCLE-FBGS-R-UMF ! Longer descriptive name for preconditioner (up to 20 chars) ML ! Preconditioner type: NONE JACOBI BJAC AS ML 1 ! Number of sweeps for 1-level, outer sweeps for ML; lines below ignored for non-ML -WCYCLE ! Type of multilevel CYCLE: VCYCLE WCYCLE KCYCLE MULT ADD +VCYCLE ! Type of multilevel CYCLE: VCYCLE WCYCLE KCYCLE MULT ADD -3 ! Max Number of levels in a multilevel preconditioner; if <0, lib default --3 ! Target coarse matrix size; if <0, lib default +-500 ! Target coarse matrix size; if <0, lib default SMOOTHED ! Type of aggregation: SMOOTHED UNSMOOTHED -SYMDEC ! Parallel aggregation: DEC, SYMDEC +DEC ! Parallel aggregation: DEC, SYMDEC NATURAL ! Ordering of aggregation NATURAL DEGREE NOFILTER ! Filtering of matrix ? FILTER NOFILTER -1.5 ! Coarsening ratio, if < 0 use library default -2 ! Number of thresholds in vector, next line ignored if <= 0 +-2 ! Number of thresholds in vector, next line ignored if <= 0 0.05 0.025 ! Thresholds -0.0100d0 ! Smoothed aggregation threshold: >= 0.0 +0.0d0 ! Smoothed aggregation threshold: >= 0.0 BJAC ! Smoother type JACOBI GS BGS BJAC AS; ignored for non-ML 1 ! Number of sweeps for (pre) smoother 0 ! Number of overlap layers for AS preconditioner HALO ! AS restriction operator: NONE HALO NONE ! AS prolongation operator: NONE SUM AVG -ILU ! Subdomain solver: JACOBI BS BGS ILU ILUT MILU MUMPS SLU UMF +GS ! Subdomain solver: JACOBI BS BGS ILU ILUT MILU MUMPS SLU UMF 0 ! Fill level P for ILU(P) and ILU(T,P) 1.d-4 ! Threshold T for ILU(T,P) -NONE ! Second (post) smoother, ignored if NONE +BJAC ! Second (post) smoother, ignored if NONE 1 ! Number of sweeps for (post) smoother 0 ! Number of overlap layers for AS preconditioner HALO ! AS restriction operator: NONE HALO NONE ! AS prolongation operator: NONE SUM AVG -JACOBI ! Subdomain solver: JACOBI BS BGS ILU ILUT MILU MUMPS SLU UMF +BWGS ! Subdomain solver: JACOBI BS BGS ILU ILUT MILU MUMPS SLU UMF 0 ! Fill level P for ILU(P) and ILU(T,P) 1.d-4 ! Threshold T for ILU(T,P) DEFLT ! Coarsest-level solver: MUMPS UMF SLU SLUDIST JACOBI GS BJAC -ILU ! Coarsest-level subsolver: ILU UMF SLU MUMPS SLUDIST +ILU ! Coarsest-level subsolver: ILU UMF SLU MUMPS SLUDIST DIST ! Coarsest-level matrix distribution: DIST REPL 1 ! Coarsest-level fillin P for ILU(P) and ILU(T,P) 1.d-4 ! Coarsest-level threshold T for ILU(T,P) diff --git a/tests/pdegen/mld_d_pde2d.f90 b/tests/pdegen/mld_d_pde2d.f90 index 2e2074b6..3260204e 100644 --- a/tests/pdegen/mld_d_pde2d.f90 +++ b/tests/pdegen/mld_d_pde2d.f90 @@ -62,8 +62,314 @@ ! data distribution. ! module mld_d_pde2d_mod + use psb_base_mod, only : psb_dpk_, psb_ipk_, psb_desc_type,& + & psb_dspmat_type, psb_d_vect_type, dzero,& + & psb_d_base_sparse_mat, psb_d_base_vect_type, psb_i_base_vect_type + + interface + function d_func_2d(x,y) result(val) + import :: psb_dpk_ + real(psb_dpk_), intent(in) :: x,y + real(psb_dpk_) :: val + end function d_func_2d + end interface + + interface mld_gen_pde2d + module procedure mld_d_gen_pde2d + end interface mld_gen_pde2d contains + function d_null_func_2d(x,y) result(val) + + real(psb_dpk_), intent(in) :: x,y + real(psb_dpk_) :: val + + val = dzero + + end function d_null_func_2d + + + ! + ! subroutine to allocate and fill in the coefficient matrix and + ! the rhs. + ! + subroutine mld_d_gen_pde2d(ictxt,idim,a,bv,xv,desc_a,afmt,& + & a1,a2,b1,b2,c,g,info,f,amold,vmold,imold,nrl) + use psb_base_mod + ! + ! Discretizes the partial differential equation + ! + ! a1 dd(u) a2 dd(u) b1 d(u) b2 d(u) + ! - ------ - ------ + ----- + ------ + c u = f + ! dxdx dydy dx dy + ! + ! with Dirichlet boundary conditions + ! u = g + ! + ! on the unit square 0<=x,y<=1. + ! + ! + ! Note that if b1=b2=c=0., the PDE is the Laplace equation. + ! + implicit none + procedure(d_func_2d) :: b1,b2,c,a1,a2,g + integer(psb_ipk_) :: idim + type(psb_dspmat_type) :: a + type(psb_d_vect_type) :: xv,bv + type(psb_desc_type) :: desc_a + integer(psb_ipk_) :: ictxt, info + character(len=*) :: afmt + procedure(d_func_2d), optional :: f + class(psb_d_base_sparse_mat), optional :: amold + class(psb_d_base_vect_type), optional :: vmold + class(psb_i_base_vect_type), optional :: imold + integer(psb_ipk_), optional :: nrl + + ! Local variables. + + integer(psb_ipk_), parameter :: nb=20 + type(psb_d_csc_sparse_mat) :: acsc + type(psb_d_coo_sparse_mat) :: acoo + type(psb_d_csr_sparse_mat) :: acsr + real(psb_dpk_) :: zt(nb),x,y,z + integer(psb_ipk_) :: m,n,nnz,glob_row,nlr,i,ii,ib,k + integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner + integer(psb_ipk_) :: np, iam, nr, nt + integer(psb_ipk_) :: icoeff + integer(psb_ipk_), allocatable :: irow(:),icol(:),myidx(:) + real(psb_dpk_), allocatable :: val(:) + ! deltah dimension of each grid cell + ! deltat discretization time + real(psb_dpk_) :: deltah, sqdeltah, deltah2 + real(psb_dpk_), parameter :: rhs=0.e0,one=1.e0,zero=0.e0 + real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb + integer(psb_ipk_) :: err_act + procedure(d_func_2d), pointer :: f_ + character(len=20) :: name, ch_err,tmpfmt + + info = psb_success_ + name = 'create_matrix' + call psb_erractionsave(err_act) + + call psb_info(ictxt, iam, np) + + + if (present(f)) then + f_ => f + else + f_ => d_null_func_2d + end if + + deltah = 1.e0/(idim+2) + sqdeltah = deltah*deltah + deltah2 = 2.e0* deltah + + ! initialize array descriptor and sparse matrix storage. provide an + ! estimate of the number of non zeroes + + m = idim*idim + n = m + nnz = ((n*7)/(np)) + if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n + + if (present(nrl)) then + nr = nrl + else + ! + ! Using a simple BLOCK distribution. + ! + nt = (m+np-1)/np + nr = max(0,min(nt,m-(iam*nt))) + end if + + nt = nr + call psb_sum(ictxt,nt) + if (nt /= m) then + write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m + info = -1 + call psb_barrier(ictxt) + call psb_abort(ictxt) + return + end if + 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(xv,desc_a,info) + if (info == psb_success_) call psb_geall(bv,desc_a,info) + 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),stat=info) + if (info /= psb_success_ ) then + info=psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + endif + + + myidx = desc_a%get_global_indices() + nlr = size(myidx) + + ! 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) + icoeff = 1 + do k=1,ib + i=ii+k-1 + ! local matrix pointer + glob_row=myidx(i) + ! compute gridpoint coordinates + if (mod(glob_row,(idim)) == 0) then + ix = glob_row/(idim) + else + ix = glob_row/(idim)+1 + endif + iy = (glob_row-(ix-1)*idim) + ! x, y + x = (ix-1)*deltah + y = (iy-1)*deltah + + zt(k) = f_(x,y) + ! internal point: build discretization + ! + ! term depending on (x-1,y) + ! + val(icoeff) = -a1(x,y)/sqdeltah-b1(x,y)/deltah2 + if (ix == 1) then + zt(k) = g(dzero,y)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-2)*idim+iy + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y-1) + val(icoeff) = -a2(x,y)/sqdeltah-b2(x,y)/deltah2 + if (iy == 1) then + zt(k) = g(x,dzero)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim+(iy-1) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + ! term depending on (x,y) + val(icoeff)=2.e0*(a1(x,y) + a2(x,y))/sqdeltah + c(x,y) + icol(icoeff) = (ix-1)*idim+iy + irow(icoeff) = glob_row + icoeff = icoeff+1 + ! term depending on (x,y+1) + val(icoeff)=-a2(x,y)/sqdeltah+b2(x,y)/deltah2 + if (iy == idim) then + zt(k) = g(x,done)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim+(iy+1) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x+1,y) + val(icoeff)=-a1(x,y)/sqdeltah+b1(x,y)/deltah2 + if (ix==idim) then + zt(k) = g(done,y)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix)*idim+(iy) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + end do + call psb_spins(icoeff-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),bv,desc_a,info) + if(info /= psb_success_) exit + zt(:)=0.e0 + 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_ + ch_err='insert rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + deallocate(val,irow,icol) + + call psb_barrier(ictxt) + t1 = psb_wtime() + call psb_cdasb(desc_a,info,mold=imold) + tcdasb = psb_wtime()-t1 + call psb_barrier(ictxt) + t1 = psb_wtime() + if (info == psb_success_) then + if (present(amold)) then + call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,mold=amold) + else + call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) + end if + end if + call psb_barrier(ictxt) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='asb rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + if (info == psb_success_) call psb_geasb(xv,desc_a,info,mold=vmold) + if (info == psb_success_) call psb_geasb(bv,desc_a,info,mold=vmold) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='asb rout.' + call psb_errpush(info,name,a_err=ch_err) + 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 + tmpfmt = a%get_fmt() + write(psb_out_unit,'("The matrix has been generated and assembled in ",a3," format.")')& + & tmpfmt + write(psb_out_unit,'("-allocation time : ",es12.5)') talc + write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen + write(psb_out_unit,'("-desc asbly time : ",es12.5)') tcdasb + write(psb_out_unit,'("- mat asbly time : ",es12.5)') tasb + write(psb_out_unit,'("-total time : ",es12.5)') ttot + + end if + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + end subroutine mld_d_gen_pde2d + + ! ! functions parametrizing the differential equation ! @@ -249,13 +555,13 @@ program mld_d_pde2d ! call psb_barrier(ictxt) t1 = psb_wtime() - call psb_gen_pde2d(ictxt,idim,a,b,x,desc_a,afmt,& + call mld_gen_pde2d(ictxt,idim,a,b,x,desc_a,afmt,& & a1,a2,b1,b2,c,g,info) call psb_barrier(ictxt) t2 = psb_wtime() - t1 if(info /= psb_success_) then info=psb_err_from_subroutine_ - ch_err='psb_gen_pde2d' + ch_err='mld_gen_pde2d' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if diff --git a/tests/pdegen/mld_d_pde3d.f90 b/tests/pdegen/mld_d_pde3d.f90 index 1c63fcd2..773f3d5f 100644 --- a/tests/pdegen/mld_d_pde3d.f90 +++ b/tests/pdegen/mld_d_pde3d.f90 @@ -63,7 +63,354 @@ ! data distribution. ! module mld_d_pde3d_mod + use psb_base_mod, only : psb_dpk_, psb_ipk_, psb_desc_type,& + & psb_dspmat_type, psb_d_vect_type, dzero,& + & psb_d_base_sparse_mat, psb_d_base_vect_type, psb_i_base_vect_type + + interface + function d_func_3d(x,y,z) result(val) + import :: psb_dpk_ + real(psb_dpk_), intent(in) :: x,y,z + real(psb_dpk_) :: val + end function d_func_3d + end interface + + interface mld_gen_pde3d + module procedure mld_d_gen_pde3d + end interface mld_gen_pde3d + + contains + + function d_null_func_3d(x,y,z) result(val) + + real(psb_dpk_), intent(in) :: x,y,z + real(psb_dpk_) :: val + + val = dzero + + end function d_null_func_3d + + ! + ! subroutine to allocate and fill in the coefficient matrix and + ! the rhs. + ! + subroutine mld_d_gen_pde3d(ictxt,idim,a,bv,xv,desc_a,afmt,& + & a1,a2,a3,b1,b2,b3,c,g,info,f,amold,vmold,imold,nrl,iv) + use psb_base_mod + ! + ! Discretizes the partial differential equation + ! + ! 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 + ! u = g + ! + ! on the unit cube 0<=x,y,z<=1. + ! + ! + ! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation. + ! + implicit none + procedure(d_func_3d) :: b1,b2,b3,c,a1,a2,a3,g + integer(psb_ipk_) :: idim + type(psb_dspmat_type) :: a + type(psb_d_vect_type) :: xv,bv + type(psb_desc_type) :: desc_a + integer(psb_ipk_) :: ictxt, info + character(len=*) :: afmt + procedure(d_func_3d), optional :: f + class(psb_d_base_sparse_mat), optional :: amold + class(psb_d_base_vect_type), optional :: vmold + class(psb_i_base_vect_type), optional :: imold + integer(psb_ipk_), optional :: nrl,iv(:) + + ! Local variables. + + integer(psb_ipk_), parameter :: nb=20 + type(psb_d_csc_sparse_mat) :: acsc + type(psb_d_coo_sparse_mat) :: acoo + type(psb_d_csr_sparse_mat) :: acsr + real(psb_dpk_) :: zt(nb),x,y,z + integer(psb_ipk_) :: m,n,nnz,glob_row,nlr,i,ii,ib,k + integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner + integer(psb_ipk_) :: np, iam, nr, nt + integer(psb_ipk_) :: icoeff + integer(psb_ipk_), allocatable :: irow(:),icol(:),myidx(:) + real(psb_dpk_), allocatable :: val(:) + ! deltah dimension of each grid cell + ! deltat discretization time + real(psb_dpk_) :: deltah, sqdeltah, deltah2 + real(psb_dpk_), parameter :: rhs=dzero,one=done,zero=dzero + real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb + integer(psb_ipk_) :: err_act + procedure(d_func_3d), pointer :: f_ + character(len=20) :: name, ch_err,tmpfmt + + info = psb_success_ + name = 'create_matrix' + call psb_erractionsave(err_act) + + call psb_info(ictxt, iam, np) + + + if (present(f)) then + f_ => f + else + f_ => d_null_func_3d + end if + + deltah = 1.d0/(idim+2) + sqdeltah = deltah*deltah + deltah2 = 2.d0* deltah + + ! initialize array descriptor and sparse matrix storage. provide an + ! estimate of the number of non zeroes + + m = idim*idim*idim + n = m + nnz = ((n*9)/(np)) + if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n + + if (.not.present(iv)) then + if (present(nrl)) then + nr = nrl + else + ! + ! Using a simple BLOCK distribution. + ! + nt = (m+np-1)/np + nr = max(0,min(nt,m-(iam*nt))) + end if + + nt = nr + call psb_sum(ictxt,nt) + if (nt /= m) then + write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m + info = -1 + call psb_barrier(ictxt) + call psb_abort(ictxt) + return + end if + else + if (size(iv) /= m) then + write(psb_err_unit,*) iam, 'Initialization error IV',size(iv),m + info = -1 + call psb_barrier(ictxt) + call psb_abort(ictxt) + return + end if + + end if + call psb_barrier(ictxt) + t0 = psb_wtime() + if (present(iv)) then + call psb_cdall(ictxt,desc_a,info,vg=iv) + else + call psb_cdall(ictxt,desc_a,info,nl=nr) + end if + + 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(xv,desc_a,info) + if (info == psb_success_) call psb_geall(bv,desc_a,info) + + 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),stat=info) + if (info /= psb_success_ ) then + info=psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + endif + + myidx = desc_a%get_global_indices() + nlr = size(myidx) + + ! 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) + icoeff = 1 + do k=1,ib + i=ii+k-1 + ! local matrix pointer + glob_row=myidx(i) + ! compute gridpoint coordinates + if (mod(glob_row,(idim*idim)) == 0) then + ix = glob_row/(idim*idim) + else + ix = glob_row/(idim*idim)+1 + endif + if (mod((glob_row-(ix-1)*idim*idim),idim) == 0) then + iy = (glob_row-(ix-1)*idim*idim)/idim + else + iy = (glob_row-(ix-1)*idim*idim)/idim+1 + endif + iz = glob_row-(ix-1)*idim*idim-(iy-1)*idim + ! x, y, x coordinates + x = (ix-1)*deltah + y = (iy-1)*deltah + z = (iz-1)*deltah + zt(k) = f_(x,y,z) + ! internal point: build discretization + ! + ! term depending on (x-1,y,z) + ! + val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2 + if (ix == 1) then + zt(k) = g(dzero,y,z)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-2)*idim*idim+(iy-1)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y-1,z) + val(icoeff) = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2 + if (iy == 1) then + zt(k) = g(x,dzero,z)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim*idim+(iy-2)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y,z-1) + val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2 + if (iz == 1) then + zt(k) = g(x,y,dzero)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz-1) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + ! term depending on (x,y,z) + val(icoeff)=2.d0*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah & + & + c(x,y,z) + icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + ! term depending on (x,y,z+1) + val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2 + if (iz == idim) then + zt(k) = g(x,y,done)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz+1) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y+1,z) + val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2 + if (iy == idim) then + zt(k) = g(x,done,z)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim*idim+(iy)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x+1,y,z) + val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2 + if (ix==idim) then + zt(k) = g(done,y,z)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix)*idim*idim+(iy-1)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + end do + call psb_spins(icoeff-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),bv,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_ + ch_err='insert rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + deallocate(val,irow,icol) + + call psb_barrier(ictxt) + t1 = psb_wtime() + call psb_cdasb(desc_a,info,mold=imold) + tcdasb = psb_wtime()-t1 + call psb_barrier(ictxt) + t1 = psb_wtime() + if (info == psb_success_) then + if (present(amold)) then + call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,mold=amold) + else + call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) + end if + end if + call psb_barrier(ictxt) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='asb rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + if (info == psb_success_) call psb_geasb(xv,desc_a,info,mold=vmold) + if (info == psb_success_) call psb_geasb(bv,desc_a,info,mold=vmold) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='asb rout.' + call psb_errpush(info,name,a_err=ch_err) + 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 + tmpfmt = a%get_fmt() + write(psb_out_unit,'("The matrix has been generated and assembled in ",a3," format.")')& + & tmpfmt + write(psb_out_unit,'("-allocation time : ",es12.5)') talc + write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen + write(psb_out_unit,'("-desc asbly time : ",es12.5)') tcdasb + write(psb_out_unit,'("- mat asbly time : ",es12.5)') tasb + write(psb_out_unit,'("-total time : ",es12.5)') ttot + + end if + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + end subroutine mld_d_gen_pde3d + ! ! functions parametrizing the differential equation ! @@ -262,13 +609,13 @@ program mld_d_pde3d call psb_barrier(ictxt) t1 = psb_wtime() - call psb_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,& + call mld_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 info=psb_err_from_subroutine_ - ch_err='psb_gen_pdeed' + ch_err='mld_gen_pde3d' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if diff --git a/tests/pdegen/mld_s_pde2d.f90 b/tests/pdegen/mld_s_pde2d.f90 index 960cc129..f1412ce5 100644 --- a/tests/pdegen/mld_s_pde2d.f90 +++ b/tests/pdegen/mld_s_pde2d.f90 @@ -62,8 +62,314 @@ ! data distribution. ! module mld_s_pde2d_mod + use psb_base_mod, only : psb_spk_, psb_ipk_, psb_desc_type,& + & psb_sspmat_type, psb_s_vect_type, szero,& + & psb_s_base_sparse_mat, psb_s_base_vect_type, psb_i_base_vect_type + + interface + function s_func_2d(x,y) result(val) + import :: psb_spk_ + real(psb_spk_), intent(in) :: x,y + real(psb_spk_) :: val + end function s_func_2d + end interface + + interface mld_gen_pde2d + module procedure mld_s_gen_pde2d + end interface mld_gen_pde2d contains + function s_null_func_2d(x,y) result(val) + + real(psb_spk_), intent(in) :: x,y + real(psb_spk_) :: val + + val = szero + + end function s_null_func_2d + + + ! + ! subroutine to allocate and fill in the coefficient matrix and + ! the rhs. + ! + subroutine mld_s_gen_pde2d(ictxt,idim,a,bv,xv,desc_a,afmt,& + & a1,a2,b1,b2,c,g,info,f,amold,vmold,imold,nrl) + use psb_base_mod + ! + ! Discretizes the partial differential equation + ! + ! a1 dd(u) a2 dd(u) b1 d(u) b2 d(u) + ! - ------ - ------ + ----- + ------ + c u = f + ! dxdx dydy dx dy + ! + ! with Dirichlet boundary conditions + ! u = g + ! + ! on the unit square 0<=x,y<=1. + ! + ! + ! Note that if b1=b2=c=0., the PDE is the Laplace equation. + ! + implicit none + procedure(s_func_2d) :: b1,b2,c,a1,a2,g + integer(psb_ipk_) :: idim + type(psb_sspmat_type) :: a + type(psb_s_vect_type) :: xv,bv + type(psb_desc_type) :: desc_a + integer(psb_ipk_) :: ictxt, info + character(len=*) :: afmt + procedure(s_func_2d), optional :: f + class(psb_s_base_sparse_mat), optional :: amold + class(psb_s_base_vect_type), optional :: vmold + class(psb_i_base_vect_type), optional :: imold + integer(psb_ipk_), optional :: nrl + + ! Local variables. + + integer(psb_ipk_), parameter :: nb=20 + type(psb_s_csc_sparse_mat) :: acsc + type(psb_s_coo_sparse_mat) :: acoo + type(psb_s_csr_sparse_mat) :: acsr + real(psb_spk_) :: zt(nb),x,y,z + integer(psb_ipk_) :: m,n,nnz,glob_row,nlr,i,ii,ib,k + integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner + integer(psb_ipk_) :: np, iam, nr, nt + integer(psb_ipk_) :: icoeff + integer(psb_ipk_), allocatable :: irow(:),icol(:),myidx(:) + real(psb_spk_), allocatable :: val(:) + ! deltah dimension of each grid cell + ! deltat discretization time + real(psb_spk_) :: deltah, sqdeltah, deltah2 + real(psb_spk_), parameter :: rhs=0.e0,one=1.e0,zero=0.e0 + real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb + integer(psb_ipk_) :: err_act + procedure(s_func_2d), pointer :: f_ + character(len=20) :: name, ch_err,tmpfmt + + info = psb_success_ + name = 'create_matrix' + call psb_erractionsave(err_act) + + call psb_info(ictxt, iam, np) + + + if (present(f)) then + f_ => f + else + f_ => s_null_func_2d + end if + + deltah = 1.e0/(idim+2) + sqdeltah = deltah*deltah + deltah2 = 2.e0* deltah + + ! initialize array descriptor and sparse matrix storage. provide an + ! estimate of the number of non zeroes + + m = idim*idim + n = m + nnz = ((n*7)/(np)) + if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n + + if (present(nrl)) then + nr = nrl + else + ! + ! Using a simple BLOCK distribution. + ! + nt = (m+np-1)/np + nr = max(0,min(nt,m-(iam*nt))) + end if + + nt = nr + call psb_sum(ictxt,nt) + if (nt /= m) then + write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m + info = -1 + call psb_barrier(ictxt) + call psb_abort(ictxt) + return + end if + 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(xv,desc_a,info) + if (info == psb_success_) call psb_geall(bv,desc_a,info) + 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),stat=info) + if (info /= psb_success_ ) then + info=psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + endif + + + myidx = desc_a%get_global_indices() + nlr = size(myidx) + + ! 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) + icoeff = 1 + do k=1,ib + i=ii+k-1 + ! local matrix pointer + glob_row=myidx(i) + ! compute gridpoint coordinates + if (mod(glob_row,(idim)) == 0) then + ix = glob_row/(idim) + else + ix = glob_row/(idim)+1 + endif + iy = (glob_row-(ix-1)*idim) + ! x, y + x = (ix-1)*deltah + y = (iy-1)*deltah + + zt(k) = f_(x,y) + ! internal point: build discretization + ! + ! term depending on (x-1,y) + ! + val(icoeff) = -a1(x,y)/sqdeltah-b1(x,y)/deltah2 + if (ix == 1) then + zt(k) = g(szero,y)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-2)*idim+iy + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y-1) + val(icoeff) = -a2(x,y)/sqdeltah-b2(x,y)/deltah2 + if (iy == 1) then + zt(k) = g(x,szero)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim+(iy-1) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + ! term depending on (x,y) + val(icoeff)=2.e0*(a1(x,y) + a2(x,y))/sqdeltah + c(x,y) + icol(icoeff) = (ix-1)*idim+iy + irow(icoeff) = glob_row + icoeff = icoeff+1 + ! term depending on (x,y+1) + val(icoeff)=-a2(x,y)/sqdeltah+b2(x,y)/deltah2 + if (iy == idim) then + zt(k) = g(x,sone)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim+(iy+1) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x+1,y) + val(icoeff)=-a1(x,y)/sqdeltah+b1(x,y)/deltah2 + if (ix==idim) then + zt(k) = g(sone,y)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix)*idim+(iy) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + end do + call psb_spins(icoeff-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),bv,desc_a,info) + if(info /= psb_success_) exit + zt(:)=0.e0 + 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_ + ch_err='insert rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + deallocate(val,irow,icol) + + call psb_barrier(ictxt) + t1 = psb_wtime() + call psb_cdasb(desc_a,info,mold=imold) + tcdasb = psb_wtime()-t1 + call psb_barrier(ictxt) + t1 = psb_wtime() + if (info == psb_success_) then + if (present(amold)) then + call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,mold=amold) + else + call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) + end if + end if + call psb_barrier(ictxt) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='asb rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + if (info == psb_success_) call psb_geasb(xv,desc_a,info,mold=vmold) + if (info == psb_success_) call psb_geasb(bv,desc_a,info,mold=vmold) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='asb rout.' + call psb_errpush(info,name,a_err=ch_err) + 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 + tmpfmt = a%get_fmt() + write(psb_out_unit,'("The matrix has been generated and assembled in ",a3," format.")')& + & tmpfmt + write(psb_out_unit,'("-allocation time : ",es12.5)') talc + write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen + write(psb_out_unit,'("-desc asbly time : ",es12.5)') tcdasb + write(psb_out_unit,'("- mat asbly time : ",es12.5)') tasb + write(psb_out_unit,'("-total time : ",es12.5)') ttot + + end if + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + end subroutine mld_s_gen_pde2d + + ! ! functions parametrizing the differential equation ! @@ -249,13 +555,13 @@ program mld_s_pde2d ! call psb_barrier(ictxt) t1 = psb_wtime() - call psb_gen_pde2d(ictxt,idim,a,b,x,desc_a,afmt,& + call mld_gen_pde2d(ictxt,idim,a,b,x,desc_a,afmt,& & a1,a2,b1,b2,c,g,info) call psb_barrier(ictxt) t2 = psb_wtime() - t1 if(info /= psb_success_) then info=psb_err_from_subroutine_ - ch_err='psb_gen_pde2d' + ch_err='mld_gen_pde2d' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if diff --git a/tests/pdegen/mld_s_pde3d.f90 b/tests/pdegen/mld_s_pde3d.f90 index 4102f46f..06f633f6 100644 --- a/tests/pdegen/mld_s_pde3d.f90 +++ b/tests/pdegen/mld_s_pde3d.f90 @@ -63,7 +63,354 @@ ! data distribution. ! module mld_s_pde3d_mod + use psb_base_mod, only : psb_spk_, psb_ipk_, psb_desc_type,& + & psb_sspmat_type, psb_s_vect_type, szero,& + & psb_s_base_sparse_mat, psb_s_base_vect_type, psb_i_base_vect_type + + interface + function s_func_3d(x,y,z) result(val) + import :: psb_spk_ + real(psb_spk_), intent(in) :: x,y,z + real(psb_spk_) :: val + end function s_func_3d + end interface + + interface mld_gen_pde3d + module procedure mld_s_gen_pde3d + end interface mld_gen_pde3d + + contains + + function s_null_func_3d(x,y,z) result(val) + + real(psb_spk_), intent(in) :: x,y,z + real(psb_spk_) :: val + + val = szero + + end function s_null_func_3d + + ! + ! subroutine to allocate and fill in the coefficient matrix and + ! the rhs. + ! + subroutine mld_s_gen_pde3d(ictxt,idim,a,bv,xv,desc_a,afmt,& + & a1,a2,a3,b1,b2,b3,c,g,info,f,amold,vmold,imold,nrl,iv) + use psb_base_mod + ! + ! Discretizes the partial differential equation + ! + ! 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 + ! u = g + ! + ! on the unit cube 0<=x,y,z<=1. + ! + ! + ! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation. + ! + implicit none + procedure(s_func_3d) :: b1,b2,b3,c,a1,a2,a3,g + integer(psb_ipk_) :: idim + type(psb_sspmat_type) :: a + type(psb_s_vect_type) :: xv,bv + type(psb_desc_type) :: desc_a + integer(psb_ipk_) :: ictxt, info + character(len=*) :: afmt + procedure(s_func_3d), optional :: f + class(psb_s_base_sparse_mat), optional :: amold + class(psb_s_base_vect_type), optional :: vmold + class(psb_i_base_vect_type), optional :: imold + integer(psb_ipk_), optional :: nrl,iv(:) + + ! Local variables. + + integer(psb_ipk_), parameter :: nb=20 + type(psb_s_csc_sparse_mat) :: acsc + type(psb_s_coo_sparse_mat) :: acoo + type(psb_s_csr_sparse_mat) :: acsr + real(psb_spk_) :: zt(nb),x,y,z + integer(psb_ipk_) :: m,n,nnz,glob_row,nlr,i,ii,ib,k + integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner + integer(psb_ipk_) :: np, iam, nr, nt + integer(psb_ipk_) :: icoeff + integer(psb_ipk_), allocatable :: irow(:),icol(:),myidx(:) + real(psb_spk_), allocatable :: val(:) + ! deltah dimension of each grid cell + ! deltat discretization time + real(psb_spk_) :: deltah, sqdeltah, deltah2 + real(psb_spk_), parameter :: rhs=szero,one=sone,zero=szero + real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb + integer(psb_ipk_) :: err_act + procedure(s_func_3d), pointer :: f_ + character(len=20) :: name, ch_err,tmpfmt + + info = psb_success_ + name = 'create_matrix' + call psb_erractionsave(err_act) + + call psb_info(ictxt, iam, np) + + + if (present(f)) then + f_ => f + else + f_ => s_null_func_3d + end if + + deltah = 1.d0/(idim+2) + sqdeltah = deltah*deltah + deltah2 = 2.d0* deltah + + ! initialize array descriptor and sparse matrix storage. provide an + ! estimate of the number of non zeroes + + m = idim*idim*idim + n = m + nnz = ((n*9)/(np)) + if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n + + if (.not.present(iv)) then + if (present(nrl)) then + nr = nrl + else + ! + ! Using a simple BLOCK distribution. + ! + nt = (m+np-1)/np + nr = max(0,min(nt,m-(iam*nt))) + end if + + nt = nr + call psb_sum(ictxt,nt) + if (nt /= m) then + write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m + info = -1 + call psb_barrier(ictxt) + call psb_abort(ictxt) + return + end if + else + if (size(iv) /= m) then + write(psb_err_unit,*) iam, 'Initialization error IV',size(iv),m + info = -1 + call psb_barrier(ictxt) + call psb_abort(ictxt) + return + end if + + end if + call psb_barrier(ictxt) + t0 = psb_wtime() + if (present(iv)) then + call psb_cdall(ictxt,desc_a,info,vg=iv) + else + call psb_cdall(ictxt,desc_a,info,nl=nr) + end if + + 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(xv,desc_a,info) + if (info == psb_success_) call psb_geall(bv,desc_a,info) + + 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),stat=info) + if (info /= psb_success_ ) then + info=psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + endif + + myidx = desc_a%get_global_indices() + nlr = size(myidx) + + ! 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) + icoeff = 1 + do k=1,ib + i=ii+k-1 + ! local matrix pointer + glob_row=myidx(i) + ! compute gridpoint coordinates + if (mod(glob_row,(idim*idim)) == 0) then + ix = glob_row/(idim*idim) + else + ix = glob_row/(idim*idim)+1 + endif + if (mod((glob_row-(ix-1)*idim*idim),idim) == 0) then + iy = (glob_row-(ix-1)*idim*idim)/idim + else + iy = (glob_row-(ix-1)*idim*idim)/idim+1 + endif + iz = glob_row-(ix-1)*idim*idim-(iy-1)*idim + ! x, y, x coordinates + x = (ix-1)*deltah + y = (iy-1)*deltah + z = (iz-1)*deltah + zt(k) = f_(x,y,z) + ! internal point: build discretization + ! + ! term depending on (x-1,y,z) + ! + val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2 + if (ix == 1) then + zt(k) = g(szero,y,z)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-2)*idim*idim+(iy-1)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y-1,z) + val(icoeff) = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2 + if (iy == 1) then + zt(k) = g(x,szero,z)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim*idim+(iy-2)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y,z-1) + val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2 + if (iz == 1) then + zt(k) = g(x,y,szero)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz-1) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + ! term depending on (x,y,z) + val(icoeff)=2.d0*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah & + & + c(x,y,z) + icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + ! term depending on (x,y,z+1) + val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2 + if (iz == idim) then + zt(k) = g(x,y,sone)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz+1) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y+1,z) + val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2 + if (iy == idim) then + zt(k) = g(x,sone,z)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim*idim+(iy)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x+1,y,z) + val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2 + if (ix==idim) then + zt(k) = g(sone,y,z)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix)*idim*idim+(iy-1)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + end do + call psb_spins(icoeff-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),bv,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_ + ch_err='insert rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + deallocate(val,irow,icol) + + call psb_barrier(ictxt) + t1 = psb_wtime() + call psb_cdasb(desc_a,info,mold=imold) + tcdasb = psb_wtime()-t1 + call psb_barrier(ictxt) + t1 = psb_wtime() + if (info == psb_success_) then + if (present(amold)) then + call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,mold=amold) + else + call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) + end if + end if + call psb_barrier(ictxt) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='asb rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + if (info == psb_success_) call psb_geasb(xv,desc_a,info,mold=vmold) + if (info == psb_success_) call psb_geasb(bv,desc_a,info,mold=vmold) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='asb rout.' + call psb_errpush(info,name,a_err=ch_err) + 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 + tmpfmt = a%get_fmt() + write(psb_out_unit,'("The matrix has been generated and assembled in ",a3," format.")')& + & tmpfmt + write(psb_out_unit,'("-allocation time : ",es12.5)') talc + write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen + write(psb_out_unit,'("-desc asbly time : ",es12.5)') tcdasb + write(psb_out_unit,'("- mat asbly time : ",es12.5)') tasb + write(psb_out_unit,'("-total time : ",es12.5)') ttot + + end if + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + end subroutine mld_s_gen_pde3d + ! ! functions parametrizing the differential equation ! @@ -262,13 +609,13 @@ program mld_s_pde3d call psb_barrier(ictxt) t1 = psb_wtime() - call psb_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,& + call mld_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 info=psb_err_from_subroutine_ - ch_err='psb_gen_pdeed' + ch_err='mld_gen_pde3d' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if