From 80c02a507eafbcc4674a43e7cc367b2b8ad932b0 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Wed, 15 Feb 2012 18:37:13 +0000 Subject: [PATCH] psblas: prec/impl/psb_d_nullprec_impl.f90 prec/psb_d_nullprec.f90 prec/psb_s_nullprec.f90 prec/psb_z_nullprec.f90 test/fileread/runs/dfs.inp test/pargen/ppde.f90 test/pargen/runs/ppde.inp util/Makefile util/psb_d_genmat_impl.f90 util/psb_d_genmat_mod.f90 util/psb_util_mod.f90 Fixed silly bug in nullprec. Started split genmat from ppde test, should move to have shared mat generation for 3D and 2D. --- prec/impl/psb_d_nullprec_impl.f90 | 3 - prec/psb_d_nullprec.f90 | 4 +- prec/psb_s_nullprec.f90 | 4 +- prec/psb_z_nullprec.f90 | 4 +- test/fileread/runs/dfs.inp | 2 +- test/pargen/ppde.f90 | 179 +++++++++---------- test/pargen/runs/ppde.inp | 4 +- util/Makefile | 7 +- util/psb_d_genmat_impl.f90 | 284 ++++++++++++++++++++++++++++++ util/psb_d_genmat_mod.f90 | 41 +++++ util/psb_util_mod.f90 | 1 + 11 files changed, 426 insertions(+), 107 deletions(-) create mode 100644 util/psb_d_genmat_impl.f90 create mode 100644 util/psb_d_genmat_mod.f90 diff --git a/prec/impl/psb_d_nullprec_impl.f90 b/prec/impl/psb_d_nullprec_impl.f90 index cc47dc51..6b6412bc 100644 --- a/prec/impl/psb_d_nullprec_impl.f90 +++ b/prec/impl/psb_d_nullprec_impl.f90 @@ -15,9 +15,6 @@ subroutine psb_d_null_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) call psb_erractionsave(err_act) - ! - ! This is the base version and we should throw an error. - ! Or should it be the NULL preonditioner??? ! info = psb_success_ diff --git a/prec/psb_d_nullprec.f90 b/prec/psb_d_nullprec.f90 index d956f195..f3473553 100644 --- a/prec/psb_d_nullprec.f90 +++ b/prec/psb_d_nullprec.f90 @@ -4,8 +4,8 @@ module psb_d_nullprec type, extends(psb_d_base_prec_type) :: psb_d_null_prec_type contains - procedure, pass(prec) :: c_apply_v => psb_d_null_apply_vect - procedure, pass(prec) :: c_apply => psb_d_null_apply + procedure, pass(prec) :: d_apply_v => psb_d_null_apply_vect + procedure, pass(prec) :: d_apply => psb_d_null_apply procedure, pass(prec) :: precbld => psb_d_null_precbld procedure, pass(prec) :: precinit => psb_d_null_precinit procedure, pass(prec) :: precseti => psb_d_null_precseti diff --git a/prec/psb_s_nullprec.f90 b/prec/psb_s_nullprec.f90 index e460b46d..0d998bdc 100644 --- a/prec/psb_s_nullprec.f90 +++ b/prec/psb_s_nullprec.f90 @@ -4,8 +4,8 @@ module psb_s_nullprec type, extends(psb_s_base_prec_type) :: psb_s_null_prec_type contains - procedure, pass(prec) :: c_apply_v => psb_s_null_apply_vect - procedure, pass(prec) :: c_apply => psb_s_null_apply + procedure, pass(prec) :: s_apply_v => psb_s_null_apply_vect + procedure, pass(prec) :: s_apply => psb_s_null_apply procedure, pass(prec) :: precbld => psb_s_null_precbld procedure, pass(prec) :: precinit => psb_s_null_precinit procedure, pass(prec) :: precseti => psb_s_null_precseti diff --git a/prec/psb_z_nullprec.f90 b/prec/psb_z_nullprec.f90 index f7bb8208..f2b10c0f 100644 --- a/prec/psb_z_nullprec.f90 +++ b/prec/psb_z_nullprec.f90 @@ -4,8 +4,8 @@ module psb_z_nullprec type, extends(psb_z_base_prec_type) :: psb_z_null_prec_type contains - procedure, pass(prec) :: c_apply_v => psb_z_null_apply_vect - procedure, pass(prec) :: c_apply => psb_z_null_apply + procedure, pass(prec) :: z_apply_v => psb_z_null_apply_vect + procedure, pass(prec) :: z_apply => psb_z_null_apply procedure, pass(prec) :: precbld => psb_z_null_precbld procedure, pass(prec) :: precinit => psb_z_null_precinit procedure, pass(prec) :: precseti => psb_z_null_precseti diff --git a/test/fileread/runs/dfs.inp b/test/fileread/runs/dfs.inp index bf3d1a1e..0acdcc26 100644 --- a/test/fileread/runs/dfs.inp +++ b/test/fileread/runs/dfs.inp @@ -7,7 +7,7 @@ BJAC Preconditioner NONE DIAG BJAC CSR Storage format CSR COO JAD 3 IPART: Partition method 0: BLK 2: graph (with Metis) 2 ISTOPC -00800 ITMAX +02100 ITMAX -1 ITRACE 30 IRST (restart for RGMRES and BiCGSTABL) 1.d-6 EPS diff --git a/test/pargen/ppde.f90 b/test/pargen/ppde.f90 index 505074f1..87ebbe0c 100644 --- a/test/pargen/ppde.f90 +++ b/test/pargen/ppde.f90 @@ -81,11 +81,10 @@ program ppde type(psb_dspmat_type) :: a type(psb_dprec_type) :: prec ! descriptor - type(psb_desc_type) :: desc_a, desc_b - ! dense matrices - type(psb_d_vect_type) :: xxv,bv, vtst - real(psb_dpk_), allocatable :: tst(:) - ! blacs parameters + type(psb_desc_type) :: desc_a + ! dense vectors + type(psb_d_vect_type) :: xxv,bv + ! parallel environment integer(psb_ipk_) :: ictxt, iam, np ! solver parameters @@ -95,12 +94,12 @@ program ppde ! other variables integer(psb_ipk_) :: info, i - character(len=20) :: name,ch_err - character(len=40) :: fname + character(len=20) :: name,ch_err + character(len=40) :: fname info=psb_success_ - + call psb_init(ictxt) call psb_info(ictxt,iam,np) @@ -129,7 +128,9 @@ program ppde ! call psb_barrier(ictxt) t1 = psb_wtime() - call create_matrix(idim,a,bv,xxv,desc_a,ictxt,afmt,info) + call gen_prob3d(ictxt,idim,a,bv,xxv,desc_a,afmt,& + & a1,a2,a3,b1,b2,b3,c,g,info) +!!$ call create_matrix(idim,a,bv,xxv,desc_a,ictxt,afmt,info) call psb_barrier(ictxt) t2 = psb_wtime() - t1 if(info /= psb_success_) then @@ -140,19 +141,6 @@ program ppde end if if (iam == psb_root_) write(psb_out_unit,'("Overall matrix creation time : ",es12.5)')t2 if (iam == psb_root_) write(psb_out_unit,'(" ")') -!!$ write(fname,'(a,i0,a)') 'pde-',idim,'.hb' -!!$ call hb_write(a,info,filename=fname,rhs=b,key='PDEGEN',mtitle='MLD2P4 pdegen Test matrix ') -!!$ write(fname,'(a,i2.2,a,i2.2,a)') 'amat-',iam,'-',np,'.mtx' -!!$ call a%print(fname) -!!$ call psb_cdprt(20+iam,desc_a,short=.false.) -!!$ call psb_cdcpy(desc_a,desc_b,info) -!!$ call psb_set_debug_level(9999) - - call psb_cdbldext(a,desc_a,2,desc_b,info,extype=psb_ovt_asov_) - if (info /= 0) then - write(0,*) 'Error from bldext' - call psb_abort(ictxt) - end if ! ! prepare the preconditioner. ! @@ -213,23 +201,9 @@ program ppde write(psb_out_unit,'("Total memory occupation for PREC: ",i12)')precsize write(psb_out_unit,'("Total memory occupation for DESC_A: ",i12)')descsize write(psb_out_unit,'("Storage type for DESC_A: ",a)') desc_a%indxmap%get_fmt() - write(psb_out_unit,'("Storage type for DESC_B: ",a)') desc_b%indxmap%get_fmt() - end if - - ! - if (.false.) then - call psb_geall(tst,desc_b, info) - call psb_geall(vtst,desc_b, info) - vtst%v%v = iam+1 - call psb_geasb(vtst,desc_b,info) - tst = vtst%get_vect() - call psb_geasb(tst,desc_b,info) - call psb_ovrl(vtst,desc_b,info,update=psb_avg_) - call psb_ovrl(tst,desc_b,info,update=psb_avg_) - write(0,*) iam,' After ovrl:',vtst%v%v - write(0,*) iam,' After ovrl:',tst end if + ! ! cleanup storage and exit ! @@ -370,18 +344,13 @@ contains ! ! 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 - use psb_mat_mod implicit none integer(psb_ipk_) :: idim - integer(psb_ipk_), parameter :: nb=20 + integer(psb_ipk_), parameter :: nb=20 type(psb_d_vect_type) :: xxv,bv type(psb_desc_type) :: desc_a integer(psb_ipk_) :: ictxt, info @@ -401,9 +370,9 @@ contains ! deltat discretization time real(psb_dpk_) :: deltah, sqdeltah, 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_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb real(psb_dpk_) :: a1, a2, a3, a4, b1, b2, b3 - external :: a1, a2, a3, a4, b1, b2, b3 + external :: a1, a2, a3, a4, b1, b2, b3 integer(psb_ipk_) :: err_act character(len=20) :: name, ch_err,tmpfmt @@ -414,7 +383,7 @@ contains call psb_info(ictxt, iam, np) - deltah = 1.d0/(idim-1) + deltah = 1.d0/(idim+2) sqdeltah = deltah*deltah deltah2 = 2.d0* deltah @@ -500,7 +469,15 @@ contains x = ix*deltah y = iy*deltah z = iz*deltah - + if (glob_row == 1) then + write(0,*) 'Starting from ',ix,iy,iz,x,y,z,deltah + end if + if (glob_row == nt) then + write(0,*) 'Ending at ',ix,iy,iz,x,y,z,deltah + end if + if (i == nlr) then + write(0,*) 'Ending at ',ix,iy,iz,x,y,z,deltah + end if ! check on boundary points zt(k) = 0.d0 ! internal point: build discretization @@ -596,6 +573,9 @@ contains call psb_barrier(ictxt) t1 = psb_wtime() call psb_cdasb(desc_a,info) + tcdasb = psb_wtime()-t1 + call psb_barrier(ictxt) + t1 = psb_wtime() if (info == psb_success_) & & call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) call psb_barrier(ictxt) @@ -627,7 +607,8 @@ contains & tmpfmt write(psb_out_unit,'("-allocation time : ",es12.5)') talc write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen - write(psb_out_unit,'("-assembly time : ",es12.5)') tasb + 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 @@ -642,51 +623,63 @@ contains end if return end subroutine create_matrix + ! + ! 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.414d0 + 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.414d0 + 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.414d0 + 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 + end function g + end program ppde -! -! 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.414d0 -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=1.414d0 -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.414d0 -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=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/80 -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/80 -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/80 -end function b3 diff --git a/test/pargen/runs/ppde.inp b/test/pargen/runs/ppde.inp index 410aa3ef..60315ead 100644 --- a/test/pargen/runs/ppde.inp +++ b/test/pargen/runs/ppde.inp @@ -2,9 +2,9 @@ BICGSTAB Iterative method BICGSTAB CGS BICG BICGSTABL RGMRES BJAC Preconditioner NONE DIAG BJAC CSR Storage format for matrix A: CSR COO JAD -040 Domain size (acutal system is this**3) +080 Domain size (acutal system is this**3) 2 Stopping criterion -0200 MAXIT +1000 MAXIT -2 ITRACE 02 IRST restart for RGMRES and BiCGSTABL diff --git a/util/Makefile b/util/Makefile index 629b4cba..864ece2c 100644 --- a/util/Makefile +++ b/util/Makefile @@ -7,7 +7,7 @@ HERE=. BASEOBJS= psb_blockpart_mod.o psb_metispart_mod.o \ psb_hbio_mod.o psb_mmio_mod.o psb_mat_dist_mod.o \ - psb_renum_mod.o psb_gps_mod.o + psb_renum_mod.o psb_gps_mod.o psb_d_genmat_mod.o IMPLOBJS= psb_s_hbio_impl.o psb_d_hbio_impl.o \ psb_c_hbio_impl.o psb_z_hbio_impl.o \ psb_s_mmio_impl.o psb_d_mmio_impl.o \ @@ -15,7 +15,8 @@ IMPLOBJS= psb_s_hbio_impl.o psb_d_hbio_impl.o \ psb_s_mat_dist_impl.o psb_d_mat_dist_impl.o \ psb_c_mat_dist_impl.o psb_z_mat_dist_impl.o \ psb_s_renum_impl.o psb_d_renum_impl.o \ - psb_c_renum_impl.o psb_z_renum_impl.o + psb_c_renum_impl.o psb_z_renum_impl.o \ + psb_d_genmat_impl.o MODOBJS=psb_util_mod.o $(BASEOBJS) COBJS=psb_amd_order.o OBJS=$(MODOBJS) $(IMPLOBJS) $(COBJS) @@ -35,6 +36,8 @@ $(HERE)/$(LIBNAME): $(OBJS) $(OBJS): $(LIBDIR)/$(BASEMODNAME)$(.mod) psb_util_mod.o: $(BASEOBJS) +$(IMPLOBJS): $(BASEOBJS) + veryclean: clean /bin/rm -f $(HERE)/$(LIBNAME) diff --git a/util/psb_d_genmat_impl.f90 b/util/psb_d_genmat_impl.f90 new file mode 100644 index 00000000..21322bf5 --- /dev/null +++ b/util/psb_d_genmat_impl.f90 @@ -0,0 +1,284 @@ +! +! subroutine to allocate and fill in the coefficient matrix and +! the rhs. +! +subroutine gen_prob3d(ictxt,idim,a,bv,xv,desc_a,afmt,a1,a2,a3,b1,b2,b3,c,g,info) + use psb_base_mod + use psb_d_genmat_mod, psb_protect_name => gen_prob3d + ! + ! 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 = 0 + ! 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 :: afmt*5 + + ! 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.d0,one=1.d0,zero=0.d0 + real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb + integer(psb_ipk_) :: err_act + + character(len=20) :: name, ch_err,tmpfmt + + info = psb_success_ + name = 'create_matrix' + call psb_erractionsave(err_act) + + call psb_info(ictxt, iam, np) + + 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 + + ! + ! 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(xv,desc_a,info) + if (info == psb_success_) call psb_geall(bv,desc_a,info) + nlr = desc_a%get_local_rows() + 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) + 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*deltah + y = iy*deltah + z = iz*deltah + zt(k) = 0.d0 + ! 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)) + 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)) + 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)) + 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*a1(x,y,z) + 2*a2(x,y,z) + 2*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)) + 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)) + 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)) + 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) + tcdasb = psb_wtime()-t1 + call psb_barrier(ictxt) + t1 = psb_wtime() + if (info == psb_success_) & + & call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) + 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) + if (info == psb_success_) call psb_geasb(bv,desc_a,info) + 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 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + return +end subroutine gen_prob3d diff --git a/util/psb_d_genmat_mod.f90 b/util/psb_d_genmat_mod.f90 new file mode 100644 index 00000000..e66bde9f --- /dev/null +++ b/util/psb_d_genmat_mod.f90 @@ -0,0 +1,41 @@ +module psb_d_genmat_mod + use psb_base_mod + 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 + subroutine gen_prob3d(ictxt,idim,a,bv,xv,desc_a,afmt,a1,a2,a3,b1,b2,b3,c,g,info) + ! + ! 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 = 0 + ! dxdx dydy dzdz dx dy dz + ! + ! with Dirichlet boundary conditions + ! u = g + ! + ! on the unit cube 0<=x,y,z<=1. + ! + ! + ! Note that if a1=a2=a3=c=0., the PDE is the Laplace equation. + ! + import :: psb_ipk_, psb_desc_type, psb_dspmat_type, psb_d_vect_type, d_func_3d + implicit none + procedure(d_func_3d) :: a1,a2,a3,c,b1,b2,b3,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 :: afmt*5 + end subroutine gen_prob3d + end interface + + +end module psb_d_genmat_mod diff --git a/util/psb_util_mod.f90 b/util/psb_util_mod.f90 index 141b9f0b..4b78019c 100644 --- a/util/psb_util_mod.f90 +++ b/util/psb_util_mod.f90 @@ -38,5 +38,6 @@ module psb_util_mod use psb_mmio_mod use psb_mat_dist_mod use psb_renum_mod + use psb_d_genmat_mod end module psb_util_mod