diff --git a/tests/pdegen/Makefile b/tests/pdegen/Makefile index 6e90a96a..0986f454 100644 --- a/tests/pdegen/Makefile +++ b/tests/pdegen/Makefile @@ -1,4 +1,4 @@ -AMGDIR=../.. +AMGDIR=../../.. AMGINCDIR=$(AMGDIR)/include include $(AMGINCDIR)/Make.inc.amg4psblas AMGMODDIR=$(AMGDIR)/modules @@ -9,42 +9,44 @@ FINCLUDES=$(FMFLAG). $(FMFLAG)$(AMGMODDIR) $(FMFLAG)$(AMGINCDIR) $(PSBLAS_INCLUD LINKOPT= EXEDIR=./runs -all: amg_s_pde3d amg_d_pde3d amg_s_pde2d amg_d_pde2d +all: amg_s_pde3d amg_d_pde3d amg_s_pde2d amg_d_pde2d -amg_d_pde3d: amg_d_pde3d.o data_input.o - $(FLINK) $(LINKOPT) amg_d_pde3d.o data_input.o -o amg_d_pde3d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS) +amg_d_pde3d: amg_d_pde3d.o amg_d_genpde_mod.o amg_d_pde3d_base_mod.o amg_d_pde3d_exp_mod.o amg_d_pde3d_gauss_mod.o data_input.o + $(FLINK) $(LINKOPT) amg_d_pde3d.o amg_d_genpde_mod.o amg_d_pde3d_base_mod.o amg_d_pde3d_exp_mod.o amg_d_pde3d_gauss_mod.o data_input.o -o amg_d_pde3d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS) /bin/mv amg_d_pde3d $(EXEDIR) -amg_s_pde3d: amg_s_pde3d.o data_input.o - $(FLINK) $(LINKOPT) amg_s_pde3d.o data_input.o -o amg_s_pde3d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS) +amg_s_pde3d: amg_s_pde3d.o amg_s_genpde_mod.o amg_s_pde3d_base_mod.o amg_s_pde3d_exp_mod.o amg_s_pde3d_gauss_mod.o data_input.o + $(FLINK) $(LINKOPT) amg_s_pde3d.o amg_s_genpde_mod.o amg_s_pde3d_base_mod.o amg_s_pde3d_exp_mod.o amg_s_pde3d_gauss_mod.o data_input.o -o amg_s_pde3d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS) /bin/mv amg_s_pde3d $(EXEDIR) -amg_d_pde2d: amg_d_pde2d.o data_input.o - $(FLINK) $(LINKOPT) amg_d_pde2d.o data_input.o -o amg_d_pde2d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS) +amg_d_pde2d: amg_d_pde2d.o amg_d_genpde_mod.o amg_d_pde2d_base_mod.o amg_d_pde2d_exp_mod.o amg_d_pde2d_box_mod.o data_input.o + $(FLINK) $(LINKOPT) amg_d_pde2d.o amg_d_genpde_mod.o amg_d_pde2d_base_mod.o amg_d_pde2d_exp_mod.o amg_d_pde2d_box_mod.o data_input.o -o amg_d_pde2d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS) /bin/mv amg_d_pde2d $(EXEDIR) -amg_s_pde2d: amg_s_pde2d.o data_input.o - $(FLINK) $(LINKOPT) amg_s_pde2d.o data_input.o -o amg_s_pde2d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS) +amg_s_pde2d: amg_s_pde2d.o amg_s_genpde_mod.o amg_s_pde2d_base_mod.o amg_s_pde2d_exp_mod.o amg_s_pde2d_box_mod.o data_input.o + $(FLINK) $(LINKOPT) amg_s_pde2d.o amg_s_genpde_mod.o amg_s_pde2d_base_mod.o amg_s_pde2d_exp_mod.o amg_s_pde2d_box_mod.o data_input.o -o amg_s_pde2d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS) /bin/mv amg_s_pde2d $(EXEDIR) amg_d_pde3d_rebld: amg_d_pde3d_rebld.o data_input.o $(FLINK) $(LINKOPT) amg_d_pde3d_rebld.o data_input.o -o amg_d_pde3d_rebld $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS) /bin/mv amg_d_pde3d_rebld $(EXEDIR) -amg_d_pde3d.o amg_s_pde3d.o amg_d_pde2d.o amg_s_pde2d.o: data_input.o +amg_d_pde3d.o amg_s_pde3d.o amg_d_pde2d.o amg_s_pde2d.o: data_input.o + +amg_d_pde3d.o: amg_d_genpde_mod.o amg_d_pde3d_base_mod.o amg_d_pde3d_exp_mod.o amg_d_pde3d_gauss_mod.o +amg_s_pde3d.o: amg_s_genpde_mod.o amg_s_pde3d_base_mod.o amg_s_pde3d_exp_mod.o amg_s_pde3d_gauss_mod.o +amg_d_pde2d.o: amg_d_genpde_mod.o amg_d_pde2d_base_mod.o amg_d_pde2d_exp_mod.o amg_d_pde2d_box_mod.o +amg_s_pde2d.o: amg_s_genpde_mod.o amg_s_pde2d_base_mod.o amg_s_pde2d_exp_mod.o amg_s_pde2d_box_mod.o check: all cd runs && ./amg_d_pde2d f + else + f_ => d_null_func_3d + end if + + if (present(partition)) then + if ((1<= partition).and.(partition <= 3)) then + partition_ = partition + else + write(*,*) 'Invalid partition choice ',partition,' defaulting to 3' + partition_ = 3 + end if + else + partition_ = 3 + end if + deltah = done/(idim+2) + sqdeltah = deltah*deltah + deltah2 = 2.0_psb_dpk_* deltah + + if (present(partition)) then + if ((1<= partition).and.(partition <= 3)) then + partition_ = partition + else + write(*,*) 'Invalid partition choice ',partition,' defaulting to 3' + partition_ = 3 + end if + else + partition_ = 3 + end if + + ! initialize array descriptor and sparse matrix storage. provide an + ! estimate of the number of non zeroes + + m = (1_psb_lpk_*idim)*idim*idim + n = m + nnz = 7*((n+np-1)/np) + if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n + t0 = psb_wtime() + select case(partition_) + case(1) + ! A BLOCK partition + 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(ctxt,nt) + if (nt /= m) then + write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + + ! + ! First example of use of CDALL: specify for each process a number of + ! contiguous rows + ! + call psb_cdall(ctxt,desc_a,info,nl=nr) + myidx = desc_a%get_global_indices() + nlr = size(myidx) + + case(2) + ! A partition defined by the user through IV + + if (present(iv)) then + if (size(iv) /= m) then + write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + else + write(psb_err_unit,*) iam, 'Initialization error: IV not present' + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + + ! + ! Second example of use of CDALL: specify for each row the + ! process that owns it + ! + call psb_cdall(ctxt,desc_a,info,vg=iv) + myidx = desc_a%get_global_indices() + nlr = size(myidx) + + case(3) + ! A 3-dimensional partition + + ! A nifty MPI function will split the process list + npdims = 0 + call mpi_dims_create(np,3,npdims,info) + npx = npdims(1) + npy = npdims(2) + npz = npdims(3) + + allocate(bndx(0:npx),bndy(0:npy),bndz(0:npz)) + ! We can reuse idx2ijk for process indices as well. + call idx2ijk(iamx,iamy,iamz,iam,npx,npy,npz,base=0) + ! Now let's split the 3D cube in hexahedra + call dist1Didx(bndx,idim,npx) + mynx = bndx(iamx+1)-bndx(iamx) + call dist1Didx(bndy,idim,npy) + myny = bndy(iamy+1)-bndy(iamy) + call dist1Didx(bndz,idim,npz) + mynz = bndz(iamz+1)-bndz(iamz) + + ! How many indices do I own? + nlr = mynx*myny*mynz + allocate(myidx(nlr)) + ! Now, let's generate the list of indices I own + nr = 0 + do i=bndx(iamx),bndx(iamx+1)-1 + do j=bndy(iamy),bndy(iamy+1)-1 + do k=bndz(iamz),bndz(iamz+1)-1 + nr = nr + 1 + call ijk2idx(myidx(nr),i,j,k,idim,idim,idim) + end do + end do + end do + if (nr /= nlr) then + write(psb_err_unit,*) iam,iamx,iamy,iamz, 'Initialization error: NR vs NLR ',& + & nr,nlr,mynx,myny,mynz + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + end if + + ! + ! Third example of use of CDALL: specify for each process + ! the set of global indices it owns. + ! + call psb_cdall(ctxt,desc_a,info,vl=myidx) + + case default + write(psb_err_unit,*) iam, 'Initialization error: should not get here' + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end select + + + 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(ctxt) + 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 + + + ! loop over rows belonging to current process in a block + ! distribution. + + call psb_barrier(ctxt) + 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 + call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim) + ! x, y, z 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 + call ijk2idx(icol(icoeff),ix-1,iy,iz,idim,idim,idim) + 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 + call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim) + 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 + call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + ! term depending on (x,y,z) + val(icoeff)=(2*done)*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah & + & + c(x,y,z) + call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim) + 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 + call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim) + 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 + call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim) + 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 + call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim) + 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(:)=dzero + 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(ctxt) + t1 = psb_wtime() + call psb_cdasb(desc_a,info) + tcdasb = psb_wtime()-t1 + call psb_barrier(ctxt) + 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(ctxt) + 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(ctxt) + ttot = psb_wtime() - t0 + + call psb_amx(ctxt,talc) + call psb_amx(ctxt,tgen) + call psb_amx(ctxt,tasb) + call psb_amx(ctxt,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(ctxt) + return + end if + return + end subroutine amg_d_gen_pde3d + + + + ! + ! subroutine to allocate and fill in the coefficient matrix and + ! the rhs. + ! + subroutine amg_d_gen_pde2d(ctxt,idim,a,bv,xv,desc_a,afmt,& + & a1,a2,b1,b2,c,g,info,f,amold,vmold,partition, nrl,iv) + use psb_base_mod + use psb_util_mod + ! + ! Discretizes the partial differential equation + ! + ! d d(u) d d(u) b1 d(u) b2 d(u) + ! - -- a1 ---- - -- a1 ---- + ----- + ------ + c u = f + ! dx dx dy dy 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_) :: info + type(psb_ctxt_type) :: ctxt + character :: afmt*5 + procedure(d_func_2d), optional :: f + class(psb_d_base_sparse_mat), optional :: amold + class(psb_d_base_vect_type), optional :: vmold + integer(psb_ipk_), optional :: partition, 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,xph,xmh,yph,ymh,zph,zmh + integer(psb_ipk_) :: nnz,nr,nlr,i,j,ii,ib,k, partition_ + integer(psb_lpk_) :: m,n,glob_row,nt + integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner + ! For 2D partition + ! Note: integer control variables going directly into an MPI call + ! must be 4 bytes, i.e. psb_mpk_ + integer(psb_mpk_) :: npdims(2), npp, minfo + integer(psb_ipk_) :: npx,npy,iamx,iamy,mynx,myny + integer(psb_ipk_), allocatable :: bndx(:),bndy(:) + ! Process grid + integer(psb_ipk_) :: np, iam + integer(psb_ipk_) :: icoeff + integer(psb_lpk_), allocatable :: irow(:),icol(:),myidx(:) + real(psb_dpk_), allocatable :: val(:) + ! deltah dimension of each grid cell + ! deltat discretization time + real(psb_dpk_) :: deltah, sqdeltah, deltah2, dd + real(psb_dpk_), parameter :: rhs=0.d0,one=done,zero=0.d0 + 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(ctxt, iam, np) + + + if (present(f)) then + f_ => f + else + f_ => d_null_func_2d + end if + + deltah = done/(idim+2) + sqdeltah = deltah*deltah + deltah2 = 2.0_psb_dpk_* deltah + + + if (present(partition)) then + if ((1<= partition).and.(partition <= 3)) then + partition_ = partition + else + write(*,*) 'Invalid partition choice ',partition,' defaulting to 3' + partition_ = 3 + end if + else + partition_ = 3 + end if + + ! initialize array descriptor and sparse matrix storage. provide an + ! estimate of the number of non zeroes + + m = (1_psb_lpk_)*idim*idim + n = m + nnz = 7*((n+np-1)/np) + if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n + t0 = psb_wtime() + select case(partition_) + case(1) + ! A BLOCK partition + 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(ctxt,nt) + if (nt /= m) then + write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + + ! + ! First example of use of CDALL: specify for each process a number of + ! contiguous rows + ! + call psb_cdall(ctxt,desc_a,info,nl=nr) + myidx = desc_a%get_global_indices() + nlr = size(myidx) + + case(2) + ! A partition defined by the user through IV + + if (present(iv)) then + if (size(iv) /= m) then + write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + else + write(psb_err_unit,*) iam, 'Initialization error: IV not present' + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + + ! + ! Second example of use of CDALL: specify for each row the + ! process that owns it + ! + call psb_cdall(ctxt,desc_a,info,vg=iv) + myidx = desc_a%get_global_indices() + nlr = size(myidx) + + case(3) + ! A 2-dimensional partition + + ! A nifty MPI function will split the process list + npdims = 0 + call mpi_dims_create(np,2,npdims,info) + npx = npdims(1) + npy = npdims(2) + + allocate(bndx(0:npx),bndy(0:npy)) + ! We can reuse idx2ijk for process indices as well. + call idx2ijk(iamx,iamy,iam,npx,npy,base=0) + ! Now let's split the 2D square in rectangles + call dist1Didx(bndx,idim,npx) + mynx = bndx(iamx+1)-bndx(iamx) + call dist1Didx(bndy,idim,npy) + myny = bndy(iamy+1)-bndy(iamy) + + ! How many indices do I own? + nlr = mynx*myny + allocate(myidx(nlr)) + ! Now, let's generate the list of indices I own + nr = 0 + do i=bndx(iamx),bndx(iamx+1)-1 + do j=bndy(iamy),bndy(iamy+1)-1 + nr = nr + 1 + call ijk2idx(myidx(nr),i,j,idim,idim) + end do + end do + if (nr /= nlr) then + write(psb_err_unit,*) iam,iamx,iamy, 'Initialization error: NR vs NLR ',& + & nr,nlr,mynx,myny + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + end if + + ! + ! Third example of use of CDALL: specify for each process + ! the set of global indices it owns. + ! + call psb_cdall(ctxt,desc_a,info,vl=myidx) + + case default + write(psb_err_unit,*) iam, 'Initialization error: should not get here' + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end select + + + 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(ctxt) + 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 + + + ! loop over rows belonging to current process in a block + ! distribution. + + call psb_barrier(ctxt) + 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 + call idx2ijk(ix,iy,glob_row,idim,idim) + ! x, y coordinates + 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 + call ijk2idx(icol(icoeff),ix-1,iy,idim,idim) + 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 + call ijk2idx(icol(icoeff),ix,iy-1,idim,idim) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + ! term depending on (x,y) + val(icoeff)=(2*done)*(a1(x,y) + a2(x,y))/sqdeltah + c(x,y) + call ijk2idx(icol(icoeff),ix,iy,idim,idim) + 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 + call ijk2idx(icol(icoeff),ix,iy+1,idim,idim) + 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 + call ijk2idx(icol(icoeff),ix+1,iy,idim,idim) + 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(:)=dzero + 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(ctxt) + t1 = psb_wtime() + call psb_cdasb(desc_a,info) + tcdasb = psb_wtime()-t1 + call psb_barrier(ctxt) + 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(ctxt) + 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(ctxt) + ttot = psb_wtime() - t0 + + call psb_amx(ctxt,talc) + call psb_amx(ctxt,tgen) + call psb_amx(ctxt,tasb) + call psb_amx(ctxt,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(ctxt) + return + end if + return + end subroutine amg_d_gen_pde2d +end module amg_d_genpde_mod diff --git a/tests/pdegen/amg_d_pde2d.f90 b/tests/pdegen/amg_d_pde2d.f90 index 524c55ab..061df38e 100644 --- a/tests/pdegen/amg_d_pde2d.f90 +++ b/tests/pdegen/amg_d_pde2d.f90 @@ -1,15 +1,15 @@ -! -! +! +! ! AMG4PSBLAS version 1.0 ! Algebraic Multigrid Package ! based on PSBLAS (Parallel Sparse BLAS version 3.5) -! -! (C) Copyright 2020 -! -! Salvatore Filippone -! Pasqua D'Ambra -! Fabio Durastante -! +! +! (C) Copyright 2020 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -21,7 +21,7 @@ ! 3. The name of the AMG4PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -33,23 +33,23 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! +! +! ! File: amg_d_pde2d.f90 ! ! Program: amg_d_pde2d ! This sample program solves a linear system obtained by discretizing a -! PDE with Dirichlet BCs. -! +! PDE with Dirichlet BCs. +! ! ! The PDE is a general second order equation in 2d ! -! a1 dd(u) a2 dd(u) b1 d(u) b2 d(u) +! a1 dd(u) a2 dd(u) b1 d(u) b2 d(u) ! - ------ - ------ ----- + ------ + c u = f -! dxdx dydy dx dy +! dxdx dydy dx dy ! ! with Dirichlet boundary conditions -! u = g +! u = g ! ! on the unit square 0<=x,y<=1. ! @@ -63,495 +63,25 @@ ! 3. A 2D distribution in which the unit square is partitioned ! into rectangles, each one assigned to a process. ! -module amg_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 amg_gen_pde2d - module procedure amg_d_gen_pde2d - end interface amg_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 - - ! - ! functions parametrizing the differential equation - ! - - ! - ! Note: b1 and b2 are the coefficients of the first - ! derivative of the unknown function. The default - ! we apply here is to have them zero, so that the resulting - ! matrix is symmetric/hermitian and suitable for - ! testing with CG and FCG. - ! When testing methods for non-hermitian matrices you can - ! change the B1/B2 functions to e.g. done/sqrt((2*done)) - ! - function b1(x,y) - use psb_base_mod, only : psb_dpk_, done, dzero - implicit none - real(psb_dpk_) :: b1 - real(psb_dpk_), intent(in) :: x,y - b1=dzero - end function b1 - function b2(x,y) - use psb_base_mod, only : psb_dpk_, done, dzero - implicit none - real(psb_dpk_) :: b2 - real(psb_dpk_), intent(in) :: x,y - b2=dzero - end function b2 - function c(x,y) - use psb_base_mod, only : psb_dpk_, done, dzero - implicit none - real(psb_dpk_) :: c - real(psb_dpk_), intent(in) :: x,y - c=0.d0 - end function c - function a1(x,y) - use psb_base_mod, only : psb_dpk_, done, dzero - implicit none - real(psb_dpk_) :: a1 - real(psb_dpk_), intent(in) :: x,y - a1=done/80 - end function a1 - function a2(x,y) - use psb_base_mod, only : psb_dpk_, done, dzero - implicit none - real(psb_dpk_) :: a2 - real(psb_dpk_), intent(in) :: x,y - a2=done/80 - end function a2 - function g(x,y) - use psb_base_mod, only : psb_dpk_, done, dzero - implicit none - real(psb_dpk_) :: g - real(psb_dpk_), intent(in) :: x,y - g = dzero - if (x == done) then - g = done - else if (x == dzero) then - g = exp(-y**2) - end if - end function g - - - ! - ! subroutine to allocate and fill in the coefficient matrix and - ! the rhs. - ! - subroutine amg_d_gen_pde2d(ctxt,idim,a,bv,xv,desc_a,afmt,info,& - & f,amold,vmold,imold,partition,nrl,iv) - use psb_base_mod - use psb_util_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 - integer(psb_ipk_) :: idim - type(psb_dspmat_type) :: a - type(psb_d_vect_type) :: xv,bv - type(psb_desc_type) :: desc_a - type(psb_ctxt_type) :: ctxt - integer(psb_ipk_) :: 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 :: partition, 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_) :: nnz,nr,nlr,i,j,ii,ib,k, partition_ - integer(psb_lpk_) :: m,n,glob_row,nt - integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner - ! For 2D partition - ! Note: integer control variables going directly into an MPI call - ! must be 4 bytes, i.e. psb_mpk_ - integer(psb_mpk_) :: npdims(2), npp, minfo - integer(psb_ipk_) :: npx,npy,iamx,iamy,mynx,myny - integer(psb_ipk_), allocatable :: bndx(:),bndy(:) - ! Process grid - integer(psb_ipk_) :: np, iam - integer(psb_ipk_) :: icoeff - integer(psb_lpk_), 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_2d), pointer :: f_ - character(len=20) :: name, ch_err,tmpfmt - - info = psb_success_ - name = 'create_matrix' - call psb_erractionsave(err_act) - - call psb_info(ctxt, iam, np) - - - if (present(f)) then - f_ => f - else - f_ => d_null_func_2d - end if - - deltah = done/(idim+1) - sqdeltah = deltah*deltah - deltah2 = (2*done)* deltah - - if (present(partition)) then - if ((1<= partition).and.(partition <= 3)) then - partition_ = partition - else - write(*,*) 'Invalid partition choice ',partition,' defaulting to 3' - partition_ = 3 - end if - else - partition_ = 3 - end if - - ! initialize array descriptor and sparse matrix storage. provide an - ! estimate of the number of non zeroes - - m = (1_psb_lpk_)*idim*idim - n = m - nnz = 7*((n+np-1)/np) - if (iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n - t0 = psb_wtime() - select case(partition_) - case(1) - ! A BLOCK partition - 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(ctxt,nt) - if (nt /= m) then - write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - return - end if - - ! - ! First example of use of CDALL: specify for each process a number of - ! contiguous rows - ! - call psb_cdall(ctxt,desc_a,info,nl=nr) - myidx = desc_a%get_global_indices() - nlr = size(myidx) - - case(2) - ! A partition defined by the user through IV - - if (present(iv)) then - if (size(iv) /= m) then - write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - return - end if - else - write(psb_err_unit,*) iam, 'Initialization error: IV not present' - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - return - end if - - ! - ! Second example of use of CDALL: specify for each row the - ! process that owns it - ! - call psb_cdall(ctxt,desc_a,info,vg=iv) - myidx = desc_a%get_global_indices() - nlr = size(myidx) - - case(3) - ! A 2-dimensional partition - - ! A nifty MPI function will split the process list - npdims = 0 - call mpi_dims_create(np,2,npdims,info) - npx = npdims(1) - npy = npdims(2) - - allocate(bndx(0:npx),bndy(0:npy)) - ! We can reuse idx2ijk for process indices as well. - call idx2ijk(iamx,iamy,iam,npx,npy,base=0) - ! Now let's split the 2D square in rectangles - call dist1Didx(bndx,idim,npx) - mynx = bndx(iamx+1)-bndx(iamx) - call dist1Didx(bndy,idim,npy) - myny = bndy(iamy+1)-bndy(iamy) - - ! How many indices do I own? - nlr = mynx*myny - allocate(myidx(nlr)) - ! Now, let's generate the list of indices I own - nr = 0 - do i=bndx(iamx),bndx(iamx+1)-1 - do j=bndy(iamy),bndy(iamy+1)-1 - nr = nr + 1 - call ijk2idx(myidx(nr),i,j,idim,idim) - end do - end do - if (nr /= nlr) then - write(psb_err_unit,*) iam,iamx,iamy, 'Initialization error: NR vs NLR ',& - & nr,nlr,mynx,myny - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - end if - - ! - ! Third example of use of CDALL: specify for each process - ! the set of global indices it owns. - ! - call psb_cdall(ctxt,desc_a,info,vl=myidx) - - case default - write(psb_err_unit,*) iam, 'Initialization error: should not get here' - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - return - end select - - - 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(ctxt) - 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 - - - ! loop over rows belonging to current process in a block - ! distribution. - - call psb_barrier(ctxt) - 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 - call idx2ijk(ix,iy,glob_row,idim,idim) - ! x, y coordinates - 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 - call ijk2idx(icol(icoeff),ix-1,iy,idim,idim) - 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 - call ijk2idx(icol(icoeff),ix,iy-1,idim,idim) - irow(icoeff) = glob_row - icoeff = icoeff+1 - endif - - ! term depending on (x,y) - val(icoeff)=(2*done)*(a1(x,y) + a2(x,y))/sqdeltah + c(x,y) - call ijk2idx(icol(icoeff),ix,iy,idim,idim) - 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 - call ijk2idx(icol(icoeff),ix,iy+1,idim,idim) - 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 - call ijk2idx(icol(icoeff),ix+1,iy,idim,idim) - 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(:)=dzero - 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(ctxt) - t1 = psb_wtime() - call psb_cdasb(desc_a,info,mold=imold) - tcdasb = psb_wtime()-t1 - call psb_barrier(ctxt) - 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(ctxt) - 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(ctxt) - ttot = psb_wtime() - t0 - - call psb_amx(ctxt,talc) - call psb_amx(ctxt,tgen) - call psb_amx(ctxt,tasb) - call psb_amx(ctxt,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(ctxt,err_act) - - return - end subroutine amg_d_gen_pde2d - -end module amg_d_pde2d_mod - - program amg_d_pde2d use psb_base_mod use amg_prec_mod use psb_krylov_mod use psb_util_mod use data_input - use amg_d_pde2d_mod + use amg_d_pde2d_base_mod + use amg_d_pde2d_exp_mod + use amg_d_pde2d_box_mod + use amg_d_genpde_mod implicit none ! input parameters character(len=20) :: kmethd, ptype - character(len=5) :: afmt + character(len=5) :: afmt, pdecoeff integer(psb_ipk_) :: idim integer(psb_epk_) :: system_size - ! miscellaneous + ! miscellaneous real(psb_dpk_) :: t1, t2, tprec, thier, tslv ! sparse matrix and preconditioner @@ -651,7 +181,7 @@ program amg_d_pde2d call psb_init(ctxt) call psb_info(ctxt,iam,np) - if (iam < 0) then + if (iam < 0) then ! This should not happen, but just in case call psb_exit(ctxt) stop @@ -662,22 +192,37 @@ program amg_d_pde2d ! ! Hello world ! - if (iam == psb_root_) then - write(*,*) 'Welcome to MLD2P4 version: ',amg_version_string_ + if (iam == psb_root_) then + write(*,*) 'Welcome to AMG4PSBLAS version: ',amg_version_string_ write(*,*) 'This is the ',trim(name),' sample program' end if ! ! get parameters ! - call get_parms(ctxt,afmt,idim,s_choice,p_choice) + call get_parms(ctxt,afmt,idim,s_choice,p_choice,pdecoeff) ! - ! allocate and fill in the coefficient matrix, rhs and initial guess + ! allocate and fill in the coefficient matrix, rhs and initial guess ! call psb_barrier(ctxt) t1 = psb_wtime() - call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,info) + select case(psb_toupper(trim(pdecoeff))) + case("CONST") + call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,& + & a1,a2,b1,b2,c,g,info) + case("EXP") + call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,& + & a1_exp,a2_exp,b1_exp,b2_exp,c_exp,g_exp,info) + case("BOX") + call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,& + & a1_box,a2_box,b1_box,b2_box,c_box,g_box,info) + case default + info=psb_err_from_subroutine_ + ch_err='amg_gen_pdecoeff' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end select call psb_barrier(ctxt) t2 = psb_wtime() - t1 if(info /= psb_success_) then @@ -687,6 +232,8 @@ program amg_d_pde2d goto 9999 end if + if (iam == psb_root_) & + & write(psb_out_unit,'("PDE Coefficients : ",a)')pdecoeff if (iam == psb_root_) & & write(psb_out_unit,'("Overall matrix creation time : ",es12.5)')t2 if (iam == psb_root_) & @@ -702,7 +249,7 @@ program amg_d_pde2d case ('JACOBI','L1-JACOBI','GS','FWGS','FBGS') ! 1-level sweeps from "outer_sweeps" call prec%set('smoother_sweeps', p_choice%jsweeps, info) - + case ('BJAC') call prec%set('smoother_sweeps', p_choice%jsweeps, info) call prec%set('sub_solve', p_choice%solve, info) @@ -717,8 +264,8 @@ program amg_d_pde2d call prec%set('sub_solve', p_choice%solve, info) call prec%set('sub_fillin', p_choice%fill, info) call prec%set('sub_iluthrs', p_choice%thr, info) - - case ('ML') + + case ('ML') ! multilevel preconditioner call prec%set('ml_cycle', p_choice%mlcycle, info) @@ -783,7 +330,7 @@ program amg_d_pde2d call prec%set('coarse_sweeps', p_choice%cjswp, info) end select - + ! build the preconditioner call psb_barrier(ctxt) t1 = psb_wtime() @@ -813,7 +360,7 @@ program amg_d_pde2d end if ! - ! iterative method parameters + ! iterative method parameters ! call psb_barrier(ctxt) t1 = psb_wtime() @@ -853,9 +400,10 @@ program amg_d_pde2d call psb_sum(ctxt,descsize) call psb_sum(ctxt,precsize) call prec%descr(iout=psb_out_unit) - if (iam == psb_root_) then + if (iam == psb_root_) then write(psb_out_unit,'("Computed solution on ",i8," processors")') np write(psb_out_unit,'("Linear system size : ",i12)') system_size + write(psb_out_unit,'("PDE Coefficients : ",a)') trim(pdecoeff) write(psb_out_unit,'("Krylov method : ",a)') trim(s_choice%kmethd) write(psb_out_unit,'("Preconditioner : ",a)') trim(p_choice%descr) write(psb_out_unit,'("Iterations to convergence : ",i12)') iter @@ -877,7 +425,7 @@ program amg_d_pde2d end if - ! + ! ! cleanup storage and exit ! call psb_gefree(b,desc_a,info) @@ -904,7 +452,7 @@ contains ! ! get iteration parameters from standard input ! - subroutine get_parms(ctxt,afmt,idim,solve,prec) + subroutine get_parms(ctxt,afmt,idim,solve,prec,pdecoeff) implicit none @@ -913,6 +461,7 @@ contains character(len=*) :: afmt type(solverdata) :: solve type(precdata) :: prec + character(len=*) :: pdecoeff integer(psb_ipk_) :: iam, nm, np, inp_unit character(len=1024) :: filename @@ -937,6 +486,7 @@ contains ! call read_data(afmt,inp_unit) ! matrix storage format call read_data(idim,inp_unit) ! Discretization grid size + call read_data(pdecoeff,inp_unit) ! PDE Coefficients ! Krylov solver data call read_data(solve%kmethd,inp_unit) ! Krylov solver call read_data(solve%istopc,inp_unit) ! stopping criterion @@ -998,6 +548,7 @@ contains call psb_bcast(ctxt,afmt) call psb_bcast(ctxt,idim) + call psb_bcast(ctxt,pdecoeff) call psb_bcast(ctxt,solve%kmethd) call psb_bcast(ctxt,solve%istopc) @@ -1010,7 +561,7 @@ contains call psb_bcast(ctxt,prec%ptype) ! broadcast first (pre-)smoother / 1-lev prec data - call psb_bcast(ctxt,prec%smther) + call psb_bcast(ctxt,prec%smther) call psb_bcast(ctxt,prec%jsweeps) call psb_bcast(ctxt,prec%novr) call psb_bcast(ctxt,prec%restr) @@ -1018,7 +569,7 @@ contains call psb_bcast(ctxt,prec%solve) call psb_bcast(ctxt,prec%fill) call psb_bcast(ctxt,prec%thr) - ! broadcast second (post-)smoother + ! broadcast second (post-)smoother call psb_bcast(ctxt,prec%smther2) call psb_bcast(ctxt,prec%jsweeps2) call psb_bcast(ctxt,prec%novr2) @@ -1027,12 +578,12 @@ contains call psb_bcast(ctxt,prec%solve2) call psb_bcast(ctxt,prec%fill2) call psb_bcast(ctxt,prec%thr2) - + ! broadcast AMG parameters call psb_bcast(ctxt,prec%mlcycle) call psb_bcast(ctxt,prec%outer_sweeps) call psb_bcast(ctxt,prec%maxlevs) - + call psb_bcast(ctxt,prec%aggr_prol) call psb_bcast(ctxt,prec%par_aggr_alg) call psb_bcast(ctxt,prec%aggr_ord) @@ -1044,7 +595,7 @@ contains call psb_bcast(ctxt,prec%athresv) end if call psb_bcast(ctxt,prec%athres) - + call psb_bcast(ctxt,prec%csize) call psb_bcast(ctxt,prec%cmat) call psb_bcast(ctxt,prec%csolve) diff --git a/tests/pdegen/amg_d_pde2d_base_mod.f90 b/tests/pdegen/amg_d_pde2d_base_mod.f90 new file mode 100644 index 00000000..c716d789 --- /dev/null +++ b/tests/pdegen/amg_d_pde2d_base_mod.f90 @@ -0,0 +1,53 @@ +module amg_d_pde2d_base_mod + use psb_base_mod, only : psb_dpk_, dzero, done + real(psb_dpk_), save, private :: epsilon=done/80 +contains + subroutine pde_set_parm(dat) + real(psb_dpk_), intent(in) :: dat + epsilon = dat + end subroutine pde_set_parm + ! + ! functions parametrizing the differential equation + ! + function b1(x,y) + use psb_base_mod, only : psb_dpk_, dzero, done + real(psb_dpk_) :: b1 + real(psb_dpk_), intent(in) :: x,y + b1 = dzero/1.414_psb_dpk_ + end function b1 + function b2(x,y) + use psb_base_mod, only : psb_dpk_, dzero, done + real(psb_dpk_) :: b2 + real(psb_dpk_), intent(in) :: x,y + b2 = dzero/1.414_psb_dpk_ + end function b2 + function c(x,y) + use psb_base_mod, only : psb_dpk_, dzero, done + real(psb_dpk_) :: c + real(psb_dpk_), intent(in) :: x,y + c = dzero + end function c + function a1(x,y) + use psb_base_mod, only : psb_dpk_, dzero, done + real(psb_dpk_) :: a1 + real(psb_dpk_), intent(in) :: x,y + a1=done*epsilon + end function a1 + function a2(x,y) + use psb_base_mod, only : psb_dpk_, dzero, done + real(psb_dpk_) :: a2 + real(psb_dpk_), intent(in) :: x,y + a2=done*epsilon + end function a2 + function g(x,y) + use psb_base_mod, only : psb_dpk_, dzero, done + real(psb_dpk_) :: g + real(psb_dpk_), intent(in) :: x,y + g = dzero + if (x == done) then + g = done + else if (x == dzero) then + g = done + end if + end function g +end module amg_d_pde2d_base_mod diff --git a/tests/pdegen/amg_d_pde2d_box_mod.f90 b/tests/pdegen/amg_d_pde2d_box_mod.f90 new file mode 100644 index 00000000..b5518be2 --- /dev/null +++ b/tests/pdegen/amg_d_pde2d_box_mod.f90 @@ -0,0 +1,53 @@ +module amg_d_pde2d_box_mod + use psb_base_mod, only : psb_dpk_, dzero, done + real(psb_dpk_), save, private :: epsilon=done/80 +contains + subroutine pde_set_parm(dat) + real(psb_dpk_), intent(in) :: dat + epsilon = dat + end subroutine pde_set_parm + ! + ! functions parametrizing the differential equation + ! + function b1_box(x,y) + use psb_base_mod, only : psb_dpk_, dzero, done + real(psb_dpk_) :: b1_box + real(psb_dpk_), intent(in) :: x,y + b1_box = done/1.414_psb_dpk_ + end function b1_box + function b2_box(x,y) + use psb_base_mod, only : psb_dpk_, dzero, done + real(psb_dpk_) :: b2_box + real(psb_dpk_), intent(in) :: x,y + b2_box = done/1.414_psb_dpk_ + end function b2_box + function c_box(x,y) + use psb_base_mod, only : psb_dpk_, dzero, done + real(psb_dpk_) :: c_box + real(psb_dpk_), intent(in) :: x,y + c_box = dzero + end function c_box + function a1_box(x,y) + use psb_base_mod, only : psb_dpk_, dzero, done + real(psb_dpk_) :: a1_box + real(psb_dpk_), intent(in) :: x,y + a1_box=done*epsilon + end function a1_box + function a2_box(x,y) + use psb_base_mod, only : psb_dpk_, dzero, done + real(psb_dpk_) :: a2_box + real(psb_dpk_), intent(in) :: x,y + a2_box=done*epsilon + end function a2_box + function g_box(x,y) + use psb_base_mod, only : psb_dpk_, dzero, done + real(psb_dpk_) :: g_box + real(psb_dpk_), intent(in) :: x,y + g_box = dzero + if (x == done) then + g_box = done + else if (x == dzero) then + g_box = done + end if + end function g_box +end module amg_d_pde2d_box_mod diff --git a/tests/pdegen/amg_d_pde2d_exp_mod.f90 b/tests/pdegen/amg_d_pde2d_exp_mod.f90 new file mode 100644 index 00000000..1929d471 --- /dev/null +++ b/tests/pdegen/amg_d_pde2d_exp_mod.f90 @@ -0,0 +1,53 @@ +module amg_d_pde2d_exp_mod + use psb_base_mod, only : psb_dpk_, done, dzero + real(psb_dpk_), save, private :: epsilon=done/80 +contains + subroutine pde_set_parm(dat) + real(psb_dpk_), intent(in) :: dat + epsilon = dat + end subroutine pde_set_parm + ! + ! functions parametrizing the differential equation + ! + function b1_exp(x,y) + use psb_base_mod, only : psb_dpk_, done, dzero + real(psb_dpk_) :: b1_exp + real(psb_dpk_), intent(in) :: x,y + b1_exp = dzero + end function b1_exp + function b2_exp(x,y) + use psb_base_mod, only : psb_dpk_, done, dzero + real(psb_dpk_) :: b2_exp + real(psb_dpk_), intent(in) :: x,y + b2_exp = dzero + end function b2_exp + function c_exp(x,y) + use psb_base_mod, only : psb_dpk_, done, dzero + real(psb_dpk_) :: c_exp + real(psb_dpk_), intent(in) :: x,y + c_exp = dzero + end function c_exp + function a1_exp(x,y) + use psb_base_mod, only : psb_dpk_, done, dzero + real(psb_dpk_) :: a1_exp + real(psb_dpk_), intent(in) :: x,y + a1=done*epsilon*exp(-(x+y)) + end function a1_exp + function a2_exp(x,y) + use psb_base_mod, only : psb_dpk_, done, dzero + real(psb_dpk_) :: a2_exp + real(psb_dpk_), intent(in) :: x,y + a2=done*epsilon*exp(-(x+y)) + end function a2_exp + function g_exp(x,y) + use psb_base_mod, only : psb_dpk_, done, dzero + real(psb_dpk_) :: g_exp + real(psb_dpk_), intent(in) :: x,y + g_exp = dzero + if (x == done) then + g_exp = done + else if (x == dzero) then + g_exp = done + end if + end function g_exp +end module amg_d_pde2d_exp_mod diff --git a/tests/pdegen/amg_d_pde3d.f90 b/tests/pdegen/amg_d_pde3d.f90 index 43e3591a..8c04d686 100644 --- a/tests/pdegen/amg_d_pde3d.f90 +++ b/tests/pdegen/amg_d_pde3d.f90 @@ -1,15 +1,15 @@ -! -! +! +! ! AMG4PSBLAS version 1.0 ! Algebraic Multigrid Package ! based on PSBLAS (Parallel Sparse BLAS version 3.5) -! -! (C) Copyright 2020 -! -! Salvatore Filippone -! Pasqua D'Ambra -! Fabio Durastante -! +! +! (C) Copyright 2020 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -21,7 +21,7 @@ ! 3. The name of the AMG4PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -33,24 +33,24 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! +! +! ! ! File: amg_d_pde3d.f90 ! ! Program: amg_d_pde3d ! This sample program solves a linear system obtained by discretizing a -! PDE with Dirichlet BCs. -! +! PDE with Dirichlet BCs. +! ! ! The PDE is a general second order equation in 3d ! -! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u) +! 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 +! dxdx dydy dzdz dx dy dz ! ! with Dirichlet boundary conditions -! u = g +! u = g ! ! on the unit cube 0<=x,y,z<=1. ! @@ -64,534 +64,25 @@ ! 3. A 3D distribution in which the unit cube is partitioned ! into subcubes, each one assigned to a process. ! -module amg_d_pde3d_mod - use psb_base_mod, only : psb_dpk_, psb_ipk_, psb_lpk_, 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, psb_l_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 amg_gen_pde3d - module procedure amg_d_gen_pde3d - end interface amg_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 - ! - ! functions parametrizing the differential equation - ! - ! - ! Note: b1, b2 and b3 are the coefficients of the first - ! derivative of the unknown function. The default - ! we apply here is to have them zero, so that the resulting - ! matrix is symmetric/hermitian and suitable for - ! testing with CG and FCG. - ! When testing methods for non-hermitian matrices you can - ! change the B1/B2/B3 functions to e.g. done/sqrt((3*done)) - ! - function b1(x,y,z) - use psb_base_mod, only : psb_dpk_, done, dzero - implicit none - real(psb_dpk_) :: b1 - real(psb_dpk_), intent(in) :: x,y,z - b1=dzero - end function b1 - function b2(x,y,z) - use psb_base_mod, only : psb_dpk_, done, dzero - implicit none - real(psb_dpk_) :: b2 - real(psb_dpk_), intent(in) :: x,y,z - b2=dzero - end function b2 - function b3(x,y,z) - use psb_base_mod, only : psb_dpk_, done, dzero - implicit none - real(psb_dpk_) :: b3 - real(psb_dpk_), intent(in) :: x,y,z - - b3=dzero - end function b3 - function c(x,y,z) - use psb_base_mod, only : psb_dpk_, done, dzero - implicit none - real(psb_dpk_) :: c - real(psb_dpk_), intent(in) :: x,y,z - c=dzero - end function c - function a1(x,y,z) - use psb_base_mod, only : psb_dpk_, done, dzero - implicit none - real(psb_dpk_) :: a1 - real(psb_dpk_), intent(in) :: x,y,z - a1=done/80 - end function a1 - function a2(x,y,z) - use psb_base_mod, only : psb_dpk_, done, dzero - implicit none - real(psb_dpk_) :: a2 - real(psb_dpk_), intent(in) :: x,y,z - a2=done/80 - end function a2 - function a3(x,y,z) - use psb_base_mod, only : psb_dpk_, done, dzero - implicit none - real(psb_dpk_) :: a3 - real(psb_dpk_), intent(in) :: x,y,z - a3=done/80 - end function a3 - function g(x,y,z) - use psb_base_mod, only : psb_dpk_, done, dzero - implicit none - 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 - - - ! - ! subroutine to allocate and fill in the coefficient matrix and - ! the rhs. - ! - subroutine amg_d_gen_pde3d(ctxt,idim,a,bv,xv,desc_a,afmt,info,& - & f,amold,vmold,imold,partition,nrl,iv) - use psb_base_mod - use psb_util_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 - integer(psb_ipk_) :: idim - type(psb_dspmat_type) :: a - type(psb_d_vect_type) :: xv,bv - type(psb_desc_type) :: desc_a - type(psb_ctxt_type) :: ctxt - integer(psb_ipk_) :: 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 :: partition, 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_) :: nnz,nr,nlr,i,j,ii,ib,k, partition_ - integer(psb_lpk_) :: m,n,glob_row,nt - integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner - ! For 3D partition - ! Note: integer control variables going directly into an MPI call - ! must be 4 bytes, i.e. psb_mpk_ - integer(psb_mpk_) :: npdims(3), npp, minfo - integer(psb_ipk_) :: npx,npy,npz, iamx,iamy,iamz,mynx,myny,mynz - integer(psb_ipk_), allocatable :: bndx(:),bndy(:),bndz(:) - ! Process grid - integer(psb_ipk_) :: np, iam - integer(psb_ipk_) :: icoeff - integer(psb_lpk_), 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(ctxt, iam, np) - - - if (present(f)) then - f_ => f - else - f_ => d_null_func_3d - end if - - deltah = done/(idim+1) - sqdeltah = deltah*deltah - deltah2 = (2*done)* deltah - - if (present(partition)) then - if ((1<= partition).and.(partition <= 3)) then - partition_ = partition - else - write(*,*) 'Invalid partition choice ',partition,' defaulting to 3' - partition_ = 3 - end if - else - partition_ = 3 - end if - - ! initialize array descriptor and sparse matrix storage. provide an - ! estimate of the number of non zeroes - - m = (1_psb_lpk_*idim)*idim*idim - n = m - nnz = 7*((n+np-1)/np) - if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n - t0 = psb_wtime() - select case(partition_) - case(1) - ! A BLOCK partition - 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(ctxt,nt) - if (nt /= m) then - write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - return - end if - - ! - ! First example of use of CDALL: specify for each process a number of - ! contiguous rows - ! - call psb_cdall(ctxt,desc_a,info,nl=nr) - myidx = desc_a%get_global_indices() - nlr = size(myidx) - - case(2) - ! A partition defined by the user through IV - - if (present(iv)) then - if (size(iv) /= m) then - write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - return - end if - else - write(psb_err_unit,*) iam, 'Initialization error: IV not present' - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - return - end if - - ! - ! Second example of use of CDALL: specify for each row the - ! process that owns it - ! - call psb_cdall(ctxt,desc_a,info,vg=iv) - myidx = desc_a%get_global_indices() - nlr = size(myidx) - - case(3) - ! A 3-dimensional partition - - ! A nifty MPI function will split the process list - npdims = 0 - call mpi_dims_create(np,3,npdims,info) - npx = npdims(1) - npy = npdims(2) - npz = npdims(3) - - allocate(bndx(0:npx),bndy(0:npy),bndz(0:npz)) - ! We can reuse idx2ijk for process indices as well. - call idx2ijk(iamx,iamy,iamz,iam,npx,npy,npz,base=0) - ! Now let's split the 3D cube in hexahedra - call dist1Didx(bndx,idim,npx) - mynx = bndx(iamx+1)-bndx(iamx) - call dist1Didx(bndy,idim,npy) - myny = bndy(iamy+1)-bndy(iamy) - call dist1Didx(bndz,idim,npz) - mynz = bndz(iamz+1)-bndz(iamz) - - ! How many indices do I own? - nlr = mynx*myny*mynz - allocate(myidx(nlr)) - ! Now, let's generate the list of indices I own - nr = 0 - do i=bndx(iamx),bndx(iamx+1)-1 - do j=bndy(iamy),bndy(iamy+1)-1 - do k=bndz(iamz),bndz(iamz+1)-1 - nr = nr + 1 - call ijk2idx(myidx(nr),i,j,k,idim,idim,idim) - end do - end do - end do - if (nr /= nlr) then - write(psb_err_unit,*) iam,iamx,iamy,iamz, 'Initialization error: NR vs NLR ',& - & nr,nlr,mynx,myny,mynz - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - end if - - ! - ! Third example of use of CDALL: specify for each process - ! the set of global indices it owns. - ! - call psb_cdall(ctxt,desc_a,info,vl=myidx) - - case default - write(psb_err_unit,*) iam, 'Initialization error: should not get here' - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - return - end select - - - 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(ctxt) - 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 - - - ! loop over rows belonging to current process in a block - ! distribution. - - call psb_barrier(ctxt) - 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 - call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim) - ! x, y, z 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 - call ijk2idx(icol(icoeff),ix-1,iy,iz,idim,idim,idim) - 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 - call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim) - 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 - call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim) - irow(icoeff) = glob_row - icoeff = icoeff+1 - endif - - ! term depending on (x,y,z) - val(icoeff)=(2*done)*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah & - & + c(x,y,z) - call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim) - 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 - call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim) - 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 - call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim) - 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 - call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim) - 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(:)=dzero - 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(ctxt) - t1 = psb_wtime() - call psb_cdasb(desc_a,info,mold=imold) - tcdasb = psb_wtime()-t1 - call psb_barrier(ctxt) - 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(ctxt) - 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(ctxt) - ttot = psb_wtime() - t0 - - call psb_amx(ctxt,talc) - call psb_amx(ctxt,tgen) - call psb_amx(ctxt,tasb) - call psb_amx(ctxt,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(ctxt,err_act) - - return - end subroutine amg_d_gen_pde3d - -end module amg_d_pde3d_mod - program amg_d_pde3d use psb_base_mod use amg_prec_mod use psb_krylov_mod use psb_util_mod use data_input - use amg_d_pde3d_mod + use amg_d_pde3d_base_mod + use amg_d_pde3d_exp_mod + use amg_d_pde3d_gauss_mod + use amg_d_genpde_mod implicit none ! input parameters character(len=20) :: kmethd, ptype - character(len=5) :: afmt + character(len=5) :: afmt, pdecoeff integer(psb_ipk_) :: idim integer(psb_epk_) :: system_size - ! miscellaneous + ! miscellaneous real(psb_dpk_) :: t1, t2, tprec, thier, tslv ! sparse matrix and preconditioner @@ -691,7 +182,7 @@ program amg_d_pde3d call psb_init(ctxt) call psb_info(ctxt,iam,np) - if (iam < 0) then + if (iam < 0) then ! This should not happen, but just in case call psb_exit(ctxt) stop @@ -702,23 +193,40 @@ program amg_d_pde3d ! ! Hello world ! - if (iam == psb_root_) then - write(*,*) 'Welcome to MLD2P4 version: ',amg_version_string_ + if (iam == psb_root_) then + write(*,*) 'Welcome to AMG4PSBLAS version: ',amg_version_string_ write(*,*) 'This is the ',trim(name),' sample program' end if ! ! get parameters ! - call get_parms(ctxt,afmt,idim,s_choice,p_choice) + call get_parms(ctxt,afmt,idim,s_choice,p_choice,pdecoeff) ! - ! allocate and fill in the coefficient matrix, rhs and initial guess + ! allocate and fill in the coefficient matrix, rhs and initial guess ! call psb_barrier(ctxt) t1 = psb_wtime() - call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,info) + select case(psb_toupper(trim(pdecoeff))) + case("CONST") + call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,& + & a1,a2,a3,b1,b2,b3,c,g,info) + case("EXP") + call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,& + & a1_exp,a2_exp,a3_exp,b1_exp,b2_exp,b3_exp,c_exp,g_exp,info) + case("GAUSS") + call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,& + & a1_gauss,a2_gauss,a3_gauss,b1_gauss,b2_gauss,b3_gauss,c_gauss,g_gauss,info) + case default + info=psb_err_from_subroutine_ + ch_err='amg_gen_pdecoeff' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end select + + call psb_barrier(ctxt) t2 = psb_wtime() - t1 if(info /= psb_success_) then @@ -728,6 +236,8 @@ program amg_d_pde3d goto 9999 end if + if (iam == psb_root_) & + & write(psb_out_unit,'("PDE Coefficients : ",a)')pdecoeff if (iam == psb_root_) & & write(psb_out_unit,'("Overall matrix creation time : ",es12.5)')t2 if (iam == psb_root_) & @@ -743,7 +253,7 @@ program amg_d_pde3d case ('JACOBI','L1-JACOBI','GS','FWGS','FBGS') ! 1-level sweeps from "outer_sweeps" call prec%set('smoother_sweeps', p_choice%jsweeps, info) - + case ('BJAC') call prec%set('smoother_sweeps', p_choice%jsweeps, info) call prec%set('sub_solve', p_choice%solve, info) @@ -758,8 +268,8 @@ program amg_d_pde3d call prec%set('sub_solve', p_choice%solve, info) call prec%set('sub_fillin', p_choice%fill, info) call prec%set('sub_iluthrs', p_choice%thr, info) - - case ('ML') + + case ('ML') ! multilevel preconditioner call prec%set('ml_cycle', p_choice%mlcycle, info) @@ -824,7 +334,7 @@ program amg_d_pde3d call prec%set('coarse_sweeps', p_choice%cjswp, info) end select - + ! build the preconditioner call psb_barrier(ctxt) t1 = psb_wtime() @@ -854,7 +364,7 @@ program amg_d_pde3d end if ! - ! iterative method parameters + ! iterative method parameters ! call psb_barrier(ctxt) t1 = psb_wtime() @@ -894,9 +404,10 @@ program amg_d_pde3d call psb_sum(ctxt,descsize) call psb_sum(ctxt,precsize) call prec%descr(iout=psb_out_unit) - if (iam == psb_root_) then + if (iam == psb_root_) then write(psb_out_unit,'("Computed solution on ",i8," processors")') np write(psb_out_unit,'("Linear system size : ",i12)') system_size + write(psb_out_unit,'("PDE Coefficients : ",a)') trim(pdecoeff) write(psb_out_unit,'("Krylov method : ",a)') trim(s_choice%kmethd) write(psb_out_unit,'("Preconditioner : ",a)') trim(p_choice%descr) write(psb_out_unit,'("Iterations to convergence : ",i12)') iter @@ -918,7 +429,7 @@ program amg_d_pde3d end if - ! + ! ! cleanup storage and exit ! call psb_gefree(b,desc_a,info) @@ -945,7 +456,7 @@ contains ! ! get iteration parameters from standard input ! - subroutine get_parms(ctxt,afmt,idim,solve,prec) + subroutine get_parms(ctxt,afmt,idim,solve,prec,pdecoeff) implicit none @@ -954,6 +465,7 @@ contains character(len=*) :: afmt type(solverdata) :: solve type(precdata) :: prec + character(len=*) :: pdecoeff integer(psb_ipk_) :: iam, nm, np, inp_unit character(len=1024) :: filename @@ -978,6 +490,7 @@ contains ! call read_data(afmt,inp_unit) ! matrix storage format call read_data(idim,inp_unit) ! Discretization grid size + call read_data(pdecoeff,inp_unit) ! PDE Coefficients ! Krylov solver data call read_data(solve%kmethd,inp_unit) ! Krylov solver call read_data(solve%istopc,inp_unit) ! stopping criterion @@ -1039,6 +552,7 @@ contains call psb_bcast(ctxt,afmt) call psb_bcast(ctxt,idim) + call psb_bcast(ctxt,pdecoeff) call psb_bcast(ctxt,solve%kmethd) call psb_bcast(ctxt,solve%istopc) @@ -1051,7 +565,7 @@ contains call psb_bcast(ctxt,prec%ptype) ! broadcast first (pre-)smoother / 1-lev prec data - call psb_bcast(ctxt,prec%smther) + call psb_bcast(ctxt,prec%smther) call psb_bcast(ctxt,prec%jsweeps) call psb_bcast(ctxt,prec%novr) call psb_bcast(ctxt,prec%restr) @@ -1059,7 +573,7 @@ contains call psb_bcast(ctxt,prec%solve) call psb_bcast(ctxt,prec%fill) call psb_bcast(ctxt,prec%thr) - ! broadcast second (post-)smoother + ! broadcast second (post-)smoother call psb_bcast(ctxt,prec%smther2) call psb_bcast(ctxt,prec%jsweeps2) call psb_bcast(ctxt,prec%novr2) @@ -1068,12 +582,12 @@ contains call psb_bcast(ctxt,prec%solve2) call psb_bcast(ctxt,prec%fill2) call psb_bcast(ctxt,prec%thr2) - + ! broadcast AMG parameters call psb_bcast(ctxt,prec%mlcycle) call psb_bcast(ctxt,prec%outer_sweeps) call psb_bcast(ctxt,prec%maxlevs) - + call psb_bcast(ctxt,prec%aggr_prol) call psb_bcast(ctxt,prec%par_aggr_alg) call psb_bcast(ctxt,prec%aggr_ord) @@ -1085,7 +599,7 @@ contains call psb_bcast(ctxt,prec%athresv) end if call psb_bcast(ctxt,prec%athres) - + call psb_bcast(ctxt,prec%csize) call psb_bcast(ctxt,prec%cmat) call psb_bcast(ctxt,prec%csolve) diff --git a/tests/pdegen/amg_d_pde3d_base_mod.f90 b/tests/pdegen/amg_d_pde3d_base_mod.f90 new file mode 100644 index 00000000..76aad5cd --- /dev/null +++ b/tests/pdegen/amg_d_pde3d_base_mod.f90 @@ -0,0 +1,65 @@ +module amg_d_pde3d_base_mod + use psb_base_mod, only : psb_dpk_, done + real(psb_dpk_), save, private :: epsilon=done/80 +contains + subroutine pde_set_parm(dat) + real(psb_dpk_), intent(in) :: dat + epsilon = dat + end subroutine pde_set_parm + ! + ! functions parametrizing the differential equation + ! + function b1(x,y,z) + use psb_base_mod, only : psb_dpk_, done + real(psb_dpk_) :: b1 + real(psb_dpk_), intent(in) :: x,y,z + b1=done/sqrt(3.0_psb_dpk_) + end function b1 + function b2(x,y,z) + use psb_base_mod, only : psb_dpk_, done + real(psb_dpk_) :: b2 + real(psb_dpk_), intent(in) :: x,y,z + b2=done/sqrt(3.0_psb_dpk_) + end function b2 + function b3(x,y,z) + use psb_base_mod, only : psb_dpk_, done + real(psb_dpk_) :: b3 + real(psb_dpk_), intent(in) :: x,y,z + b3=done/sqrt(3.0_psb_dpk_) + end function b3 + function c(x,y,z) + use psb_base_mod, only : psb_dpk_, done + real(psb_dpk_) :: c + real(psb_dpk_), intent(in) :: x,y,z + c=dzero + 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=epsilon + 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=epsilon + 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=epsilon + end function a3 + function g(x,y,z) + use psb_base_mod, only : psb_dpk_, done, dzero + 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 = done + end if + end function g +end module amg_d_pde3d_base_mod diff --git a/tests/pdegen/amg_d_pde3d_exp_mod.f90 b/tests/pdegen/amg_d_pde3d_exp_mod.f90 new file mode 100644 index 00000000..fdbb2970 --- /dev/null +++ b/tests/pdegen/amg_d_pde3d_exp_mod.f90 @@ -0,0 +1,65 @@ +module amg_d_pde3d_exp_mod + use psb_base_mod, only : psb_dpk_, done + real(psb_dpk_), save, private :: epsilon=done/160 +contains + subroutine pde_set_parm(dat) + real(psb_dpk_), intent(in) :: dat + epsilon = dat + end subroutine pde_set_parm + ! + ! functions parametrizing the differential equation + ! + function b1_exp(x,y,z) + use psb_base_mod, only : psb_dpk_, dzero + real(psb_dpk_) :: b1_exp + real(psb_dpk_), intent(in) :: x,y,z + b1_exp=dzero/sqrt(3.0_psb_dpk_) + end function b1_exp + function b2_exp(x,y,z) + use psb_base_mod, only : psb_dpk_, dzero + real(psb_dpk_) :: b2_exp + real(psb_dpk_), intent(in) :: x,y,z + b2_exp=dzero/sqrt(3.0_psb_dpk_) + end function b2_exp + function b3_exp(x,y,z) + use psb_base_mod, only : psb_dpk_, dzero + real(psb_dpk_) :: b3_exp + real(psb_dpk_), intent(in) :: x,y,z + b3_exp=dzero/sqrt(3.0_psb_dpk_) + end function b3_exp + function c_exp(x,y,z) + use psb_base_mod, only : psb_dpk_, dzero + real(psb_dpk_) :: c_exp + real(psb_dpk_), intent(in) :: x,y,z + c_exp=dzero + end function c_exp + function a1_exp(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: a1_exp + real(psb_dpk_), intent(in) :: x,y,z + a1_exp=epsilon*exp(-(x+y+z)) + end function a1_exp + function a2_exp(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: a2_exp + real(psb_dpk_), intent(in) :: x,y,z + a2_exp=epsilon*exp(-(x+y+z)) + end function a2_exp + function a3_exp(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: a3_exp + real(psb_dpk_), intent(in) :: x,y,z + a3_exp=epsilon*exp(-(x+y+z)) + end function a3_exp + function g_exp(x,y,z) + use psb_base_mod, only : psb_dpk_, done, dzero + real(psb_dpk_) :: g_exp + real(psb_dpk_), intent(in) :: x,y,z + g_exp = dzero + if (x == done) then + g_exp = done + else if (x == dzero) then + g_exp = done + end if + end function g_exp +end module amg_d_pde3d_exp_mod diff --git a/tests/pdegen/amg_d_pde3d_gauss_mod.f90 b/tests/pdegen/amg_d_pde3d_gauss_mod.f90 new file mode 100644 index 00000000..3787c76a --- /dev/null +++ b/tests/pdegen/amg_d_pde3d_gauss_mod.f90 @@ -0,0 +1,65 @@ +module amg_d_pde3d_gauss_mod + use psb_base_mod, only : psb_dpk_, done + real(psb_dpk_), save, private :: epsilon=done/80 +contains + subroutine pde_set_parm(dat) + real(psb_dpk_), intent(in) :: dat + epsilon = dat + end subroutine pde_set_parm + ! + ! functions parametrizing the differential equation + ! + function b1_gauss(x,y,z) + use psb_base_mod, only : psb_dpk_, done + real(psb_dpk_) :: b1_gauss + real(psb_dpk_), intent(in) :: x,y,z + b1_gauss=done/sqrt(3.0_psb_dpk_)-2*x*exp(-(x**2+y**2+z**2)) + end function b1_gauss + function b2_gauss(x,y,z) + use psb_base_mod, only : psb_dpk_, done + real(psb_dpk_) :: b2_gauss + real(psb_dpk_), intent(in) :: x,y,z + b2_gauss=done/sqrt(3.0_psb_dpk_)-2*y*exp(-(x**2+y**2+z**2)) + end function b2_gauss + function b3_gauss(x,y,z) + use psb_base_mod, only : psb_dpk_, done + real(psb_dpk_) :: b3_gauss + real(psb_dpk_), intent(in) :: x,y,z + b3_gauss=done/sqrt(3.0_psb_dpk_)-2*z*exp(-(x**2+y**2+z**2)) + end function b3_gauss + function c_gauss(x,y,z) + use psb_base_mod, only : psb_dpk_, dzero + real(psb_dpk_) :: c_gauss + real(psb_dpk_), intent(in) :: x,y,z + c=dzero + end function c_gauss + function a1_gauss(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: a1_gauss + real(psb_dpk_), intent(in) :: x,y,z + a1_gauss=epsilon*exp(-(x**2+y**2+z**2)) + end function a1_gauss + function a2_gauss(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: a2_gauss + real(psb_dpk_), intent(in) :: x,y,z + a2_gauss=epsilon*exp(-(x**2+y**2+z**2)) + end function a2_gauss + function a3_gauss(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: a3_gauss + real(psb_dpk_), intent(in) :: x,y,z + a3_gauss=epsilon*exp(-(x**2+y**2+z**2)) + end function a3_gauss + function g_gauss(x,y,z) + use psb_base_mod, only : psb_dpk_, done, dzero + real(psb_dpk_) :: g_gauss + real(psb_dpk_), intent(in) :: x,y,z + g_gauss = dzero + if (x == done) then + g_gauss = done + else if (x == dzero) then + g_gauss = done + end if + end function g_gauss +end module amg_d_pde3d_gauss_mod diff --git a/tests/pdegen/amg_s_genpde_mod.f90 b/tests/pdegen/amg_s_genpde_mod.f90 new file mode 100644 index 00000000..a448590f --- /dev/null +++ b/tests/pdegen/amg_s_genpde_mod.f90 @@ -0,0 +1,857 @@ +module amg_s_genpde_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 amg_gen_pde3d + module procedure amg_s_gen_pde3d + end interface amg_gen_pde3d + + 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 amg_gen_pde2d + module procedure amg_s_gen_pde2d + end interface amg_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 + + 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 amg_s_gen_pde3d(ctxt,idim,a,bv,xv,desc_a,afmt,& + & a1,a2,a3,b1,b2,b3,c,g,info,f,amold,vmold,partition, nrl,iv) + use psb_base_mod + use psb_util_mod + ! + ! Discretizes the partial differential equation + ! + ! d a1 d(u) d a1 d(u) d a1 d(u) b1 d(u) b2 d(u) b3 d(u) + ! - ------ - ------ - ------ + ----- + ------ + ------ + c u = f + ! dx dx dy dy dz dz 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_) :: info + type(psb_ctxt_type) :: ctxt + character :: afmt*5 + procedure(s_func_3d), optional :: f + class(psb_s_base_sparse_mat), optional :: amold + class(psb_s_base_vect_type), optional :: vmold + integer(psb_ipk_), optional :: partition, 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,xph,xmh,yph,ymh,zph,zmh + integer(psb_ipk_) :: nnz,nr,nlr,i,j,ii,ib,k, partition_ + integer(psb_lpk_) :: m,n,glob_row,nt + integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner + ! For 3D partition + ! Note: integer control variables going directly into an MPI call + ! must be 4 bytes, i.e. psb_mpk_ + integer(psb_mpk_) :: npdims(3), npp, minfo + integer(psb_ipk_) :: npx,npy,npz, iamx,iamy,iamz,mynx,myny,mynz + integer(psb_ipk_), allocatable :: bndx(:),bndy(:),bndz(:) + ! Process grid + integer(psb_ipk_) :: np, iam + integer(psb_ipk_) :: icoeff + integer(psb_lpk_), 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 = 's_create_matrix' + call psb_erractionsave(err_act) + + call psb_info(ctxt, iam, np) + + + if (present(f)) then + f_ => f + else + f_ => s_null_func_3d + end if + + if (present(partition)) then + if ((1<= partition).and.(partition <= 3)) then + partition_ = partition + else + write(*,*) 'Invalid partition choice ',partition,' defaulting to 3' + partition_ = 3 + end if + else + partition_ = 3 + end if + deltah = sone/(idim+2) + sqdeltah = deltah*deltah + deltah2 = 2.0_psb_spk_* deltah + + if (present(partition)) then + if ((1<= partition).and.(partition <= 3)) then + partition_ = partition + else + write(*,*) 'Invalid partition choice ',partition,' defaulting to 3' + partition_ = 3 + end if + else + partition_ = 3 + end if + + ! initialize array descriptor and sparse matrix storage. provide an + ! estimate of the number of non zeroes + + m = (1_psb_lpk_*idim)*idim*idim + n = m + nnz = 7*((n+np-1)/np) + if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n + t0 = psb_wtime() + select case(partition_) + case(1) + ! A BLOCK partition + 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(ctxt,nt) + if (nt /= m) then + write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + + ! + ! First example of use of CDALL: specify for each process a number of + ! contiguous rows + ! + call psb_cdall(ctxt,desc_a,info,nl=nr) + myidx = desc_a%get_global_indices() + nlr = size(myidx) + + case(2) + ! A partition defined by the user through IV + + if (present(iv)) then + if (size(iv) /= m) then + write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + else + write(psb_err_unit,*) iam, 'Initialization error: IV not present' + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + + ! + ! Second example of use of CDALL: specify for each row the + ! process that owns it + ! + call psb_cdall(ctxt,desc_a,info,vg=iv) + myidx = desc_a%get_global_indices() + nlr = size(myidx) + + case(3) + ! A 3-dimensional partition + + ! A nifty MPI function will split the process list + npdims = 0 + call mpi_dims_create(np,3,npdims,info) + npx = npdims(1) + npy = npdims(2) + npz = npdims(3) + + allocate(bndx(0:npx),bndy(0:npy),bndz(0:npz)) + ! We can reuse idx2ijk for process indices as well. + call idx2ijk(iamx,iamy,iamz,iam,npx,npy,npz,base=0) + ! Now let's split the 3D cube in hexahedra + call dist1Didx(bndx,idim,npx) + mynx = bndx(iamx+1)-bndx(iamx) + call dist1Didx(bndy,idim,npy) + myny = bndy(iamy+1)-bndy(iamy) + call dist1Didx(bndz,idim,npz) + mynz = bndz(iamz+1)-bndz(iamz) + + ! How many indices do I own? + nlr = mynx*myny*mynz + allocate(myidx(nlr)) + ! Now, let's generate the list of indices I own + nr = 0 + do i=bndx(iamx),bndx(iamx+1)-1 + do j=bndy(iamy),bndy(iamy+1)-1 + do k=bndz(iamz),bndz(iamz+1)-1 + nr = nr + 1 + call ijk2idx(myidx(nr),i,j,k,idim,idim,idim) + end do + end do + end do + if (nr /= nlr) then + write(psb_err_unit,*) iam,iamx,iamy,iamz, 'Initialization error: NR vs NLR ',& + & nr,nlr,mynx,myny,mynz + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + end if + + ! + ! Third example of use of CDALL: specify for each process + ! the set of global indices it owns. + ! + call psb_cdall(ctxt,desc_a,info,vl=myidx) + + case default + write(psb_err_unit,*) iam, 'Initialization error: should not get here' + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end select + + + 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(ctxt) + 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 + + + ! loop over rows belonging to current process in a block + ! distribution. + + call psb_barrier(ctxt) + 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 + call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim) + ! x, y, z 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 + call ijk2idx(icol(icoeff),ix-1,iy,iz,idim,idim,idim) + 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 + call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim) + 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 + call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + ! term depending on (x,y,z) + val(icoeff)=(2*sone)*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah & + & + c(x,y,z) + call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim) + 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 + call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim) + 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 + call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim) + 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 + call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim) + 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(:)=szero + 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(ctxt) + t1 = psb_wtime() + call psb_cdasb(desc_a,info) + tcdasb = psb_wtime()-t1 + call psb_barrier(ctxt) + 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(ctxt) + 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(ctxt) + ttot = psb_wtime() - t0 + + call psb_amx(ctxt,talc) + call psb_amx(ctxt,tgen) + call psb_amx(ctxt,tasb) + call psb_amx(ctxt,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(ctxt) + return + end if + return + end subroutine amg_s_gen_pde3d + + + + ! + ! subroutine to allocate and fill in the coefficient matrix and + ! the rhs. + ! + subroutine amg_s_gen_pde2d(ctxt,idim,a,bv,xv,desc_a,afmt,& + & a1,a2,b1,b2,c,g,info,f,amold,vmold,partition, nrl,iv) + use psb_base_mod + use psb_util_mod + ! + ! Discretizes the partial differential equation + ! + ! d d(u) d d(u) b1 d(u) b2 d(u) + ! - -- a1 ---- - -- a1 ---- + ----- + ------ + c u = f + ! dx dx dy dy 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_) :: info + type(psb_ctxt_type) :: ctxt + character :: afmt*5 + procedure(s_func_2d), optional :: f + class(psb_s_base_sparse_mat), optional :: amold + class(psb_s_base_vect_type), optional :: vmold + integer(psb_ipk_), optional :: partition, 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,xph,xmh,yph,ymh,zph,zmh + integer(psb_ipk_) :: nnz,nr,nlr,i,j,ii,ib,k, partition_ + integer(psb_lpk_) :: m,n,glob_row,nt + integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner + ! For 2D partition + ! Note: integer control variables going directly into an MPI call + ! must be 4 bytes, i.e. psb_mpk_ + integer(psb_mpk_) :: npdims(2), npp, minfo + integer(psb_ipk_) :: npx,npy,iamx,iamy,mynx,myny + integer(psb_ipk_), allocatable :: bndx(:),bndy(:) + ! Process grid + integer(psb_ipk_) :: np, iam + integer(psb_ipk_) :: icoeff + integer(psb_lpk_), allocatable :: irow(:),icol(:),myidx(:) + real(psb_spk_), allocatable :: val(:) + ! deltah dimension of each grid cell + ! deltat discretization time + real(psb_spk_) :: deltah, sqdeltah, deltah2, dd + real(psb_spk_), parameter :: rhs=0.d0,one=sone,zero=0.d0 + 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(ctxt, iam, np) + + + if (present(f)) then + f_ => f + else + f_ => s_null_func_2d + end if + + deltah = sone/(idim+2) + sqdeltah = deltah*deltah + deltah2 = 2.0_psb_spk_* deltah + + + if (present(partition)) then + if ((1<= partition).and.(partition <= 3)) then + partition_ = partition + else + write(*,*) 'Invalid partition choice ',partition,' defaulting to 3' + partition_ = 3 + end if + else + partition_ = 3 + end if + + ! initialize array descriptor and sparse matrix storage. provide an + ! estimate of the number of non zeroes + + m = (1_psb_lpk_)*idim*idim + n = m + nnz = 7*((n+np-1)/np) + if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n + t0 = psb_wtime() + select case(partition_) + case(1) + ! A BLOCK partition + 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(ctxt,nt) + if (nt /= m) then + write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + + ! + ! First example of use of CDALL: specify for each process a number of + ! contiguous rows + ! + call psb_cdall(ctxt,desc_a,info,nl=nr) + myidx = desc_a%get_global_indices() + nlr = size(myidx) + + case(2) + ! A partition defined by the user through IV + + if (present(iv)) then + if (size(iv) /= m) then + write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + else + write(psb_err_unit,*) iam, 'Initialization error: IV not present' + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + + ! + ! Second example of use of CDALL: specify for each row the + ! process that owns it + ! + call psb_cdall(ctxt,desc_a,info,vg=iv) + myidx = desc_a%get_global_indices() + nlr = size(myidx) + + case(3) + ! A 2-dimensional partition + + ! A nifty MPI function will split the process list + npdims = 0 + call mpi_dims_create(np,2,npdims,info) + npx = npdims(1) + npy = npdims(2) + + allocate(bndx(0:npx),bndy(0:npy)) + ! We can reuse idx2ijk for process indices as well. + call idx2ijk(iamx,iamy,iam,npx,npy,base=0) + ! Now let's split the 2D square in rectangles + call dist1Didx(bndx,idim,npx) + mynx = bndx(iamx+1)-bndx(iamx) + call dist1Didx(bndy,idim,npy) + myny = bndy(iamy+1)-bndy(iamy) + + ! How many indices do I own? + nlr = mynx*myny + allocate(myidx(nlr)) + ! Now, let's generate the list of indices I own + nr = 0 + do i=bndx(iamx),bndx(iamx+1)-1 + do j=bndy(iamy),bndy(iamy+1)-1 + nr = nr + 1 + call ijk2idx(myidx(nr),i,j,idim,idim) + end do + end do + if (nr /= nlr) then + write(psb_err_unit,*) iam,iamx,iamy, 'Initialization error: NR vs NLR ',& + & nr,nlr,mynx,myny + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + end if + + ! + ! Third example of use of CDALL: specify for each process + ! the set of global indices it owns. + ! + call psb_cdall(ctxt,desc_a,info,vl=myidx) + + case default + write(psb_err_unit,*) iam, 'Initialization error: should not get here' + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end select + + + 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(ctxt) + 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 + + + ! loop over rows belonging to current process in a block + ! distribution. + + call psb_barrier(ctxt) + 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 + call idx2ijk(ix,iy,glob_row,idim,idim) + ! x, y coordinates + 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 + call ijk2idx(icol(icoeff),ix-1,iy,idim,idim) + 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 + call ijk2idx(icol(icoeff),ix,iy-1,idim,idim) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + ! term depending on (x,y) + val(icoeff)=(2*sone)*(a1(x,y) + a2(x,y))/sqdeltah + c(x,y) + call ijk2idx(icol(icoeff),ix,iy,idim,idim) + 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 + call ijk2idx(icol(icoeff),ix,iy+1,idim,idim) + 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 + call ijk2idx(icol(icoeff),ix+1,iy,idim,idim) + 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(:)=szero + 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(ctxt) + t1 = psb_wtime() + call psb_cdasb(desc_a,info) + tcdasb = psb_wtime()-t1 + call psb_barrier(ctxt) + 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(ctxt) + 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(ctxt) + ttot = psb_wtime() - t0 + + call psb_amx(ctxt,talc) + call psb_amx(ctxt,tgen) + call psb_amx(ctxt,tasb) + call psb_amx(ctxt,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(ctxt) + return + end if + return + end subroutine amg_s_gen_pde2d +end module amg_s_genpde_mod diff --git a/tests/pdegen/amg_s_pde2d.f90 b/tests/pdegen/amg_s_pde2d.f90 index 18078498..682b9d4b 100644 --- a/tests/pdegen/amg_s_pde2d.f90 +++ b/tests/pdegen/amg_s_pde2d.f90 @@ -1,15 +1,15 @@ -! -! +! +! ! AMG4PSBLAS version 1.0 ! Algebraic Multigrid Package ! based on PSBLAS (Parallel Sparse BLAS version 3.5) -! -! (C) Copyright 2020 -! -! Salvatore Filippone -! Pasqua D'Ambra -! Fabio Durastante -! +! +! (C) Copyright 2020 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -21,7 +21,7 @@ ! 3. The name of the AMG4PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -33,23 +33,23 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! +! +! ! File: amg_s_pde2d.f90 ! ! Program: amg_s_pde2d ! This sample program solves a linear system obtained by discretizing a -! PDE with Dirichlet BCs. -! +! PDE with Dirichlet BCs. +! ! ! The PDE is a general second order equation in 2d ! -! a1 dd(u) a2 dd(u) b1 d(u) b2 d(u) +! a1 dd(u) a2 dd(u) b1 d(u) b2 d(u) ! - ------ - ------ ----- + ------ + c u = f -! dxdx dydy dx dy +! dxdx dydy dx dy ! ! with Dirichlet boundary conditions -! u = g +! u = g ! ! on the unit square 0<=x,y<=1. ! @@ -63,495 +63,25 @@ ! 3. A 2D distribution in which the unit square is partitioned ! into rectangles, each one assigned to a process. ! -module amg_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 amg_gen_pde2d - module procedure amg_s_gen_pde2d - end interface amg_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 - - ! - ! functions parametrizing the differential equation - ! - - ! - ! Note: b1 and b2 are the coefficients of the first - ! derivative of the unknown function. The default - ! we apply here is to have them zero, so that the resulting - ! matrix is symmetric/hermitian and suitable for - ! testing with CG and FCG. - ! When testing methods for non-hermitian matrices you can - ! change the B1/B2 functions to e.g. sone/sqrt((2*sone)) - ! - function b1(x,y) - use psb_base_mod, only : psb_spk_, sone, szero - implicit none - real(psb_spk_) :: b1 - real(psb_spk_), intent(in) :: x,y - b1=szero - end function b1 - function b2(x,y) - use psb_base_mod, only : psb_spk_, sone, szero - implicit none - real(psb_spk_) :: b2 - real(psb_spk_), intent(in) :: x,y - b2=szero - end function b2 - function c(x,y) - use psb_base_mod, only : psb_spk_, sone, szero - implicit none - real(psb_spk_) :: c - real(psb_spk_), intent(in) :: x,y - c=0.d0 - end function c - function a1(x,y) - use psb_base_mod, only : psb_spk_, sone, szero - implicit none - real(psb_spk_) :: a1 - real(psb_spk_), intent(in) :: x,y - a1=sone/80 - end function a1 - function a2(x,y) - use psb_base_mod, only : psb_spk_, sone, szero - implicit none - real(psb_spk_) :: a2 - real(psb_spk_), intent(in) :: x,y - a2=sone/80 - end function a2 - function g(x,y) - use psb_base_mod, only : psb_spk_, sone, szero - implicit none - real(psb_spk_) :: g - real(psb_spk_), intent(in) :: x,y - g = szero - if (x == sone) then - g = sone - else if (x == szero) then - g = exp(-y**2) - end if - end function g - - - ! - ! subroutine to allocate and fill in the coefficient matrix and - ! the rhs. - ! - subroutine amg_s_gen_pde2d(ctxt,idim,a,bv,xv,desc_a,afmt,info,& - & f,amold,vmold,imold,partition,nrl,iv) - use psb_base_mod - use psb_util_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 - integer(psb_ipk_) :: idim - type(psb_sspmat_type) :: a - type(psb_s_vect_type) :: xv,bv - type(psb_desc_type) :: desc_a - type(psb_ctxt_type) :: ctxt - integer(psb_ipk_) :: 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 :: partition, 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_) :: nnz,nr,nlr,i,j,ii,ib,k, partition_ - integer(psb_lpk_) :: m,n,glob_row,nt - integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner - ! For 2D partition - ! Note: integer control variables going directly into an MPI call - ! must be 4 bytes, i.e. psb_mpk_ - integer(psb_mpk_) :: npdims(2), npp, minfo - integer(psb_ipk_) :: npx,npy,iamx,iamy,mynx,myny - integer(psb_ipk_), allocatable :: bndx(:),bndy(:) - ! Process grid - integer(psb_ipk_) :: np, iam - integer(psb_ipk_) :: icoeff - integer(psb_lpk_), 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_2d), pointer :: f_ - character(len=20) :: name, ch_err,tmpfmt - - info = psb_success_ - name = 'create_matrix' - call psb_erractionsave(err_act) - - call psb_info(ctxt, iam, np) - - - if (present(f)) then - f_ => f - else - f_ => s_null_func_2d - end if - - deltah = sone/(idim+1) - sqdeltah = deltah*deltah - deltah2 = (2*sone)* deltah - - if (present(partition)) then - if ((1<= partition).and.(partition <= 3)) then - partition_ = partition - else - write(*,*) 'Invalid partition choice ',partition,' defaulting to 3' - partition_ = 3 - end if - else - partition_ = 3 - end if - - ! initialize array descriptor and sparse matrix storage. provide an - ! estimate of the number of non zeroes - - m = (1_psb_lpk_)*idim*idim - n = m - nnz = 7*((n+np-1)/np) - if (iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n - t0 = psb_wtime() - select case(partition_) - case(1) - ! A BLOCK partition - 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(ctxt,nt) - if (nt /= m) then - write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - return - end if - - ! - ! First example of use of CDALL: specify for each process a number of - ! contiguous rows - ! - call psb_cdall(ctxt,desc_a,info,nl=nr) - myidx = desc_a%get_global_indices() - nlr = size(myidx) - - case(2) - ! A partition defined by the user through IV - - if (present(iv)) then - if (size(iv) /= m) then - write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - return - end if - else - write(psb_err_unit,*) iam, 'Initialization error: IV not present' - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - return - end if - - ! - ! Second example of use of CDALL: specify for each row the - ! process that owns it - ! - call psb_cdall(ctxt,desc_a,info,vg=iv) - myidx = desc_a%get_global_indices() - nlr = size(myidx) - - case(3) - ! A 2-dimensional partition - - ! A nifty MPI function will split the process list - npdims = 0 - call mpi_dims_create(np,2,npdims,info) - npx = npdims(1) - npy = npdims(2) - - allocate(bndx(0:npx),bndy(0:npy)) - ! We can reuse idx2ijk for process indices as well. - call idx2ijk(iamx,iamy,iam,npx,npy,base=0) - ! Now let's split the 2D square in rectangles - call dist1Didx(bndx,idim,npx) - mynx = bndx(iamx+1)-bndx(iamx) - call dist1Didx(bndy,idim,npy) - myny = bndy(iamy+1)-bndy(iamy) - - ! How many indices do I own? - nlr = mynx*myny - allocate(myidx(nlr)) - ! Now, let's generate the list of indices I own - nr = 0 - do i=bndx(iamx),bndx(iamx+1)-1 - do j=bndy(iamy),bndy(iamy+1)-1 - nr = nr + 1 - call ijk2idx(myidx(nr),i,j,idim,idim) - end do - end do - if (nr /= nlr) then - write(psb_err_unit,*) iam,iamx,iamy, 'Initialization error: NR vs NLR ',& - & nr,nlr,mynx,myny - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - end if - - ! - ! Third example of use of CDALL: specify for each process - ! the set of global indices it owns. - ! - call psb_cdall(ctxt,desc_a,info,vl=myidx) - - case default - write(psb_err_unit,*) iam, 'Initialization error: should not get here' - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - return - end select - - - 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(ctxt) - 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 - - - ! loop over rows belonging to current process in a block - ! distribution. - - call psb_barrier(ctxt) - 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 - call idx2ijk(ix,iy,glob_row,idim,idim) - ! x, y coordinates - 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 - call ijk2idx(icol(icoeff),ix-1,iy,idim,idim) - 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 - call ijk2idx(icol(icoeff),ix,iy-1,idim,idim) - irow(icoeff) = glob_row - icoeff = icoeff+1 - endif - - ! term depending on (x,y) - val(icoeff)=(2*sone)*(a1(x,y) + a2(x,y))/sqdeltah + c(x,y) - call ijk2idx(icol(icoeff),ix,iy,idim,idim) - 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 - call ijk2idx(icol(icoeff),ix,iy+1,idim,idim) - 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 - call ijk2idx(icol(icoeff),ix+1,iy,idim,idim) - 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(:)=szero - 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(ctxt) - t1 = psb_wtime() - call psb_cdasb(desc_a,info,mold=imold) - tcdasb = psb_wtime()-t1 - call psb_barrier(ctxt) - 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(ctxt) - 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(ctxt) - ttot = psb_wtime() - t0 - - call psb_amx(ctxt,talc) - call psb_amx(ctxt,tgen) - call psb_amx(ctxt,tasb) - call psb_amx(ctxt,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(ctxt,err_act) - - return - end subroutine amg_s_gen_pde2d - -end module amg_s_pde2d_mod - - program amg_s_pde2d use psb_base_mod use amg_prec_mod use psb_krylov_mod use psb_util_mod use data_input - use amg_s_pde2d_mod + use amg_s_pde2d_base_mod + use amg_s_pde2d_exp_mod + use amg_s_pde2d_box_mod + use amg_s_genpde_mod implicit none ! input parameters character(len=20) :: kmethd, ptype - character(len=5) :: afmt + character(len=5) :: afmt, pdecoeff integer(psb_ipk_) :: idim integer(psb_epk_) :: system_size - ! miscellaneous + ! miscellaneous real(psb_dpk_) :: t1, t2, tprec, thier, tslv ! sparse matrix and preconditioner @@ -651,7 +181,7 @@ program amg_s_pde2d call psb_init(ctxt) call psb_info(ctxt,iam,np) - if (iam < 0) then + if (iam < 0) then ! This should not happen, but just in case call psb_exit(ctxt) stop @@ -662,22 +192,37 @@ program amg_s_pde2d ! ! Hello world ! - if (iam == psb_root_) then - write(*,*) 'Welcome to MLD2P4 version: ',amg_version_string_ + if (iam == psb_root_) then + write(*,*) 'Welcome to AMG4PSBLAS version: ',amg_version_string_ write(*,*) 'This is the ',trim(name),' sample program' end if ! ! get parameters ! - call get_parms(ctxt,afmt,idim,s_choice,p_choice) + call get_parms(ctxt,afmt,idim,s_choice,p_choice,pdecoeff) ! - ! allocate and fill in the coefficient matrix, rhs and initial guess + ! allocate and fill in the coefficient matrix, rhs and initial guess ! call psb_barrier(ctxt) t1 = psb_wtime() - call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,info) + select case(psb_toupper(trim(pdecoeff))) + case("CONST") + call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,& + & a1,a2,b1,b2,c,g,info) + case("EXP") + call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,& + & a1_exp,a2_exp,b1_exp,b2_exp,c_exp,g_exp,info) + case("BOX") + call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,& + & a1_box,a2_box,b1_box,b2_box,c_box,g_box,info) + case default + info=psb_err_from_subroutine_ + ch_err='amg_gen_pdecoeff' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end select call psb_barrier(ctxt) t2 = psb_wtime() - t1 if(info /= psb_success_) then @@ -687,6 +232,8 @@ program amg_s_pde2d goto 9999 end if + if (iam == psb_root_) & + & write(psb_out_unit,'("PDE Coefficients : ",a)')pdecoeff if (iam == psb_root_) & & write(psb_out_unit,'("Overall matrix creation time : ",es12.5)')t2 if (iam == psb_root_) & @@ -702,7 +249,7 @@ program amg_s_pde2d case ('JACOBI','L1-JACOBI','GS','FWGS','FBGS') ! 1-level sweeps from "outer_sweeps" call prec%set('smoother_sweeps', p_choice%jsweeps, info) - + case ('BJAC') call prec%set('smoother_sweeps', p_choice%jsweeps, info) call prec%set('sub_solve', p_choice%solve, info) @@ -717,8 +264,8 @@ program amg_s_pde2d call prec%set('sub_solve', p_choice%solve, info) call prec%set('sub_fillin', p_choice%fill, info) call prec%set('sub_iluthrs', p_choice%thr, info) - - case ('ML') + + case ('ML') ! multilevel preconditioner call prec%set('ml_cycle', p_choice%mlcycle, info) @@ -783,7 +330,7 @@ program amg_s_pde2d call prec%set('coarse_sweeps', p_choice%cjswp, info) end select - + ! build the preconditioner call psb_barrier(ctxt) t1 = psb_wtime() @@ -813,7 +360,7 @@ program amg_s_pde2d end if ! - ! iterative method parameters + ! iterative method parameters ! call psb_barrier(ctxt) t1 = psb_wtime() @@ -853,9 +400,10 @@ program amg_s_pde2d call psb_sum(ctxt,descsize) call psb_sum(ctxt,precsize) call prec%descr(iout=psb_out_unit) - if (iam == psb_root_) then + if (iam == psb_root_) then write(psb_out_unit,'("Computed solution on ",i8," processors")') np write(psb_out_unit,'("Linear system size : ",i12)') system_size + write(psb_out_unit,'("PDE Coefficients : ",a)') trim(pdecoeff) write(psb_out_unit,'("Krylov method : ",a)') trim(s_choice%kmethd) write(psb_out_unit,'("Preconditioner : ",a)') trim(p_choice%descr) write(psb_out_unit,'("Iterations to convergence : ",i12)') iter @@ -877,7 +425,7 @@ program amg_s_pde2d end if - ! + ! ! cleanup storage and exit ! call psb_gefree(b,desc_a,info) @@ -904,7 +452,7 @@ contains ! ! get iteration parameters from standard input ! - subroutine get_parms(ctxt,afmt,idim,solve,prec) + subroutine get_parms(ctxt,afmt,idim,solve,prec,pdecoeff) implicit none @@ -913,6 +461,7 @@ contains character(len=*) :: afmt type(solverdata) :: solve type(precdata) :: prec + character(len=*) :: pdecoeff integer(psb_ipk_) :: iam, nm, np, inp_unit character(len=1024) :: filename @@ -937,6 +486,7 @@ contains ! call read_data(afmt,inp_unit) ! matrix storage format call read_data(idim,inp_unit) ! Discretization grid size + call read_data(pdecoeff,inp_unit) ! PDE Coefficients ! Krylov solver data call read_data(solve%kmethd,inp_unit) ! Krylov solver call read_data(solve%istopc,inp_unit) ! stopping criterion @@ -998,6 +548,7 @@ contains call psb_bcast(ctxt,afmt) call psb_bcast(ctxt,idim) + call psb_bcast(ctxt,pdecoeff) call psb_bcast(ctxt,solve%kmethd) call psb_bcast(ctxt,solve%istopc) @@ -1010,7 +561,7 @@ contains call psb_bcast(ctxt,prec%ptype) ! broadcast first (pre-)smoother / 1-lev prec data - call psb_bcast(ctxt,prec%smther) + call psb_bcast(ctxt,prec%smther) call psb_bcast(ctxt,prec%jsweeps) call psb_bcast(ctxt,prec%novr) call psb_bcast(ctxt,prec%restr) @@ -1018,7 +569,7 @@ contains call psb_bcast(ctxt,prec%solve) call psb_bcast(ctxt,prec%fill) call psb_bcast(ctxt,prec%thr) - ! broadcast second (post-)smoother + ! broadcast second (post-)smoother call psb_bcast(ctxt,prec%smther2) call psb_bcast(ctxt,prec%jsweeps2) call psb_bcast(ctxt,prec%novr2) @@ -1027,12 +578,12 @@ contains call psb_bcast(ctxt,prec%solve2) call psb_bcast(ctxt,prec%fill2) call psb_bcast(ctxt,prec%thr2) - + ! broadcast AMG parameters call psb_bcast(ctxt,prec%mlcycle) call psb_bcast(ctxt,prec%outer_sweeps) call psb_bcast(ctxt,prec%maxlevs) - + call psb_bcast(ctxt,prec%aggr_prol) call psb_bcast(ctxt,prec%par_aggr_alg) call psb_bcast(ctxt,prec%aggr_ord) @@ -1044,7 +595,7 @@ contains call psb_bcast(ctxt,prec%athresv) end if call psb_bcast(ctxt,prec%athres) - + call psb_bcast(ctxt,prec%csize) call psb_bcast(ctxt,prec%cmat) call psb_bcast(ctxt,prec%csolve) diff --git a/tests/pdegen/amg_s_pde2d_base_mod.f90 b/tests/pdegen/amg_s_pde2d_base_mod.f90 new file mode 100644 index 00000000..a7223a8b --- /dev/null +++ b/tests/pdegen/amg_s_pde2d_base_mod.f90 @@ -0,0 +1,53 @@ +module amg_s_pde2d_base_mod + use psb_base_mod, only : psb_spk_, szero, sone + real(psb_spk_), save, private :: epsilon=sone/80 +contains + subroutine pde_set_parm(dat) + real(psb_spk_), intent(in) :: dat + epsilon = dat + end subroutine pde_set_parm + ! + ! functions parametrizing the differential equation + ! + function b1(x,y) + use psb_base_mod, only : psb_spk_, szero, sone + real(psb_spk_) :: b1 + real(psb_spk_), intent(in) :: x,y + b1 = szero/1.414_psb_spk_ + end function b1 + function b2(x,y) + use psb_base_mod, only : psb_spk_, szero, sone + real(psb_spk_) :: b2 + real(psb_spk_), intent(in) :: x,y + b2 = szero/1.414_psb_spk_ + end function b2 + function c(x,y) + use psb_base_mod, only : psb_spk_, szero, sone + real(psb_spk_) :: c + real(psb_spk_), intent(in) :: x,y + c = szero + end function c + function a1(x,y) + use psb_base_mod, only : psb_spk_, szero, sone + real(psb_spk_) :: a1 + real(psb_spk_), intent(in) :: x,y + a1=sone*epsilon + end function a1 + function a2(x,y) + use psb_base_mod, only : psb_spk_, szero, sone + real(psb_spk_) :: a2 + real(psb_spk_), intent(in) :: x,y + a2=sone*epsilon + end function a2 + function g(x,y) + use psb_base_mod, only : psb_spk_, szero, sone + real(psb_spk_) :: g + real(psb_spk_), intent(in) :: x,y + g = szero + if (x == sone) then + g = sone + else if (x == szero) then + g = sone + end if + end function g +end module amg_s_pde2d_base_mod diff --git a/tests/pdegen/amg_s_pde2d_box_mod.f90 b/tests/pdegen/amg_s_pde2d_box_mod.f90 new file mode 100644 index 00000000..96d3f2f0 --- /dev/null +++ b/tests/pdegen/amg_s_pde2d_box_mod.f90 @@ -0,0 +1,53 @@ +module amg_s_pde2d_box_mod + use psb_base_mod, only : psb_spk_, szero, sone + real(psb_spk_), save, private :: epsilon=sone/80 +contains + subroutine pde_set_parm(dat) + real(psb_spk_), intent(in) :: dat + epsilon = dat + end subroutine pde_set_parm + ! + ! functions parametrizing the differential equation + ! + function b1_box(x,y) + use psb_base_mod, only : psb_spk_, szero, sone + real(psb_spk_) :: b1_box + real(psb_spk_), intent(in) :: x,y + b1_box = sone/1.414_psb_spk_ + end function b1_box + function b2_box(x,y) + use psb_base_mod, only : psb_spk_, szero, sone + real(psb_spk_) :: b2_box + real(psb_spk_), intent(in) :: x,y + b2_box = sone/1.414_psb_spk_ + end function b2_box + function c_box(x,y) + use psb_base_mod, only : psb_spk_, szero, sone + real(psb_spk_) :: c_box + real(psb_spk_), intent(in) :: x,y + c_box = szero + end function c_box + function a1_box(x,y) + use psb_base_mod, only : psb_spk_, szero, sone + real(psb_spk_) :: a1_box + real(psb_spk_), intent(in) :: x,y + a1_box=sone*epsilon + end function a1_box + function a2_box(x,y) + use psb_base_mod, only : psb_spk_, szero, sone + real(psb_spk_) :: a2_box + real(psb_spk_), intent(in) :: x,y + a2_box=sone*epsilon + end function a2_box + function g_box(x,y) + use psb_base_mod, only : psb_spk_, szero, sone + real(psb_spk_) :: g_box + real(psb_spk_), intent(in) :: x,y + g_box = szero + if (x == sone) then + g_box = sone + else if (x == szero) then + g_box = sone + end if + end function g_box +end module amg_s_pde2d_box_mod diff --git a/tests/pdegen/amg_s_pde2d_exp_mod.f90 b/tests/pdegen/amg_s_pde2d_exp_mod.f90 new file mode 100644 index 00000000..bcb944f7 --- /dev/null +++ b/tests/pdegen/amg_s_pde2d_exp_mod.f90 @@ -0,0 +1,53 @@ +module amg_s_pde2d_exp_mod + use psb_base_mod, only : psb_spk_, sone, szero + real(psb_spk_), save, private :: epsilon=sone/80 +contains + subroutine pde_set_parm(dat) + real(psb_spk_), intent(in) :: dat + epsilon = dat + end subroutine pde_set_parm + ! + ! functions parametrizing the differential equation + ! + function b1_exp(x,y) + use psb_base_mod, only : psb_spk_, sone, szero + real(psb_spk_) :: b1_exp + real(psb_spk_), intent(in) :: x,y + b1_exp = szero + end function b1_exp + function b2_exp(x,y) + use psb_base_mod, only : psb_spk_, sone, szero + real(psb_spk_) :: b2_exp + real(psb_spk_), intent(in) :: x,y + b2_exp = szero + end function b2_exp + function c_exp(x,y) + use psb_base_mod, only : psb_spk_, sone, szero + real(psb_spk_) :: c_exp + real(psb_spk_), intent(in) :: x,y + c_exp = szero + end function c_exp + function a1_exp(x,y) + use psb_base_mod, only : psb_spk_, sone, szero + real(psb_spk_) :: a1_exp + real(psb_spk_), intent(in) :: x,y + a1=sone*epsilon*exp(-(x+y)) + end function a1_exp + function a2_exp(x,y) + use psb_base_mod, only : psb_spk_, sone, szero + real(psb_spk_) :: a2_exp + real(psb_spk_), intent(in) :: x,y + a2=sone*epsilon*exp(-(x+y)) + end function a2_exp + function g_exp(x,y) + use psb_base_mod, only : psb_spk_, sone, szero + real(psb_spk_) :: g_exp + real(psb_spk_), intent(in) :: x,y + g_exp = szero + if (x == sone) then + g_exp = sone + else if (x == szero) then + g_exp = sone + end if + end function g_exp +end module amg_s_pde2d_exp_mod diff --git a/tests/pdegen/amg_s_pde3d.f90 b/tests/pdegen/amg_s_pde3d.f90 index 28b44ec2..f255c874 100644 --- a/tests/pdegen/amg_s_pde3d.f90 +++ b/tests/pdegen/amg_s_pde3d.f90 @@ -1,15 +1,15 @@ -! -! +! +! ! AMG4PSBLAS version 1.0 ! Algebraic Multigrid Package ! based on PSBLAS (Parallel Sparse BLAS version 3.5) -! -! (C) Copyright 2020 -! -! Salvatore Filippone -! Pasqua D'Ambra -! Fabio Durastante -! +! +! (C) Copyright 2020 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -21,7 +21,7 @@ ! 3. The name of the AMG4PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -33,24 +33,24 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! +! +! ! ! File: amg_s_pde3d.f90 ! ! Program: amg_s_pde3d ! This sample program solves a linear system obtained by discretizing a -! PDE with Dirichlet BCs. -! +! PDE with Dirichlet BCs. +! ! ! The PDE is a general second order equation in 3d ! -! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u) +! 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 +! dxdx dydy dzdz dx dy dz ! ! with Dirichlet boundary conditions -! u = g +! u = g ! ! on the unit cube 0<=x,y,z<=1. ! @@ -64,534 +64,25 @@ ! 3. A 3D distribution in which the unit cube is partitioned ! into subcubes, each one assigned to a process. ! -module amg_s_pde3d_mod - use psb_base_mod, only : psb_spk_, psb_ipk_, psb_lpk_, 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, psb_l_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 amg_gen_pde3d - module procedure amg_s_gen_pde3d - end interface amg_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 - ! - ! functions parametrizing the differential equation - ! - ! - ! Note: b1, b2 and b3 are the coefficients of the first - ! derivative of the unknown function. The default - ! we apply here is to have them zero, so that the resulting - ! matrix is symmetric/hermitian and suitable for - ! testing with CG and FCG. - ! When testing methods for non-hermitian matrices you can - ! change the B1/B2/B3 functions to e.g. sone/sqrt((3*sone)) - ! - function b1(x,y,z) - use psb_base_mod, only : psb_spk_, sone, szero - implicit none - real(psb_spk_) :: b1 - real(psb_spk_), intent(in) :: x,y,z - b1=szero - end function b1 - function b2(x,y,z) - use psb_base_mod, only : psb_spk_, sone, szero - implicit none - real(psb_spk_) :: b2 - real(psb_spk_), intent(in) :: x,y,z - b2=szero - end function b2 - function b3(x,y,z) - use psb_base_mod, only : psb_spk_, sone, szero - implicit none - real(psb_spk_) :: b3 - real(psb_spk_), intent(in) :: x,y,z - - b3=szero - end function b3 - function c(x,y,z) - use psb_base_mod, only : psb_spk_, sone, szero - implicit none - real(psb_spk_) :: c - real(psb_spk_), intent(in) :: x,y,z - c=szero - end function c - function a1(x,y,z) - use psb_base_mod, only : psb_spk_, sone, szero - implicit none - real(psb_spk_) :: a1 - real(psb_spk_), intent(in) :: x,y,z - a1=sone/80 - end function a1 - function a2(x,y,z) - use psb_base_mod, only : psb_spk_, sone, szero - implicit none - real(psb_spk_) :: a2 - real(psb_spk_), intent(in) :: x,y,z - a2=sone/80 - end function a2 - function a3(x,y,z) - use psb_base_mod, only : psb_spk_, sone, szero - implicit none - real(psb_spk_) :: a3 - real(psb_spk_), intent(in) :: x,y,z - a3=sone/80 - end function a3 - function g(x,y,z) - use psb_base_mod, only : psb_spk_, sone, szero - implicit none - 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 - end function g - - - ! - ! subroutine to allocate and fill in the coefficient matrix and - ! the rhs. - ! - subroutine amg_s_gen_pde3d(ctxt,idim,a,bv,xv,desc_a,afmt,info,& - & f,amold,vmold,imold,partition,nrl,iv) - use psb_base_mod - use psb_util_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 - integer(psb_ipk_) :: idim - type(psb_sspmat_type) :: a - type(psb_s_vect_type) :: xv,bv - type(psb_desc_type) :: desc_a - type(psb_ctxt_type) :: ctxt - integer(psb_ipk_) :: 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 :: partition, 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_) :: nnz,nr,nlr,i,j,ii,ib,k, partition_ - integer(psb_lpk_) :: m,n,glob_row,nt - integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner - ! For 3D partition - ! Note: integer control variables going directly into an MPI call - ! must be 4 bytes, i.e. psb_mpk_ - integer(psb_mpk_) :: npdims(3), npp, minfo - integer(psb_ipk_) :: npx,npy,npz, iamx,iamy,iamz,mynx,myny,mynz - integer(psb_ipk_), allocatable :: bndx(:),bndy(:),bndz(:) - ! Process grid - integer(psb_ipk_) :: np, iam - integer(psb_ipk_) :: icoeff - integer(psb_lpk_), 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(ctxt, iam, np) - - - if (present(f)) then - f_ => f - else - f_ => s_null_func_3d - end if - - deltah = sone/(idim+1) - sqdeltah = deltah*deltah - deltah2 = (2*sone)* deltah - - if (present(partition)) then - if ((1<= partition).and.(partition <= 3)) then - partition_ = partition - else - write(*,*) 'Invalid partition choice ',partition,' defaulting to 3' - partition_ = 3 - end if - else - partition_ = 3 - end if - - ! initialize array descriptor and sparse matrix storage. provide an - ! estimate of the number of non zeroes - - m = (1_psb_lpk_*idim)*idim*idim - n = m - nnz = 7*((n+np-1)/np) - if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n - t0 = psb_wtime() - select case(partition_) - case(1) - ! A BLOCK partition - 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(ctxt,nt) - if (nt /= m) then - write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - return - end if - - ! - ! First example of use of CDALL: specify for each process a number of - ! contiguous rows - ! - call psb_cdall(ctxt,desc_a,info,nl=nr) - myidx = desc_a%get_global_indices() - nlr = size(myidx) - - case(2) - ! A partition defined by the user through IV - - if (present(iv)) then - if (size(iv) /= m) then - write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - return - end if - else - write(psb_err_unit,*) iam, 'Initialization error: IV not present' - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - return - end if - - ! - ! Second example of use of CDALL: specify for each row the - ! process that owns it - ! - call psb_cdall(ctxt,desc_a,info,vg=iv) - myidx = desc_a%get_global_indices() - nlr = size(myidx) - - case(3) - ! A 3-dimensional partition - - ! A nifty MPI function will split the process list - npdims = 0 - call mpi_dims_create(np,3,npdims,info) - npx = npdims(1) - npy = npdims(2) - npz = npdims(3) - - allocate(bndx(0:npx),bndy(0:npy),bndz(0:npz)) - ! We can reuse idx2ijk for process indices as well. - call idx2ijk(iamx,iamy,iamz,iam,npx,npy,npz,base=0) - ! Now let's split the 3D cube in hexahedra - call dist1Didx(bndx,idim,npx) - mynx = bndx(iamx+1)-bndx(iamx) - call dist1Didx(bndy,idim,npy) - myny = bndy(iamy+1)-bndy(iamy) - call dist1Didx(bndz,idim,npz) - mynz = bndz(iamz+1)-bndz(iamz) - - ! How many indices do I own? - nlr = mynx*myny*mynz - allocate(myidx(nlr)) - ! Now, let's generate the list of indices I own - nr = 0 - do i=bndx(iamx),bndx(iamx+1)-1 - do j=bndy(iamy),bndy(iamy+1)-1 - do k=bndz(iamz),bndz(iamz+1)-1 - nr = nr + 1 - call ijk2idx(myidx(nr),i,j,k,idim,idim,idim) - end do - end do - end do - if (nr /= nlr) then - write(psb_err_unit,*) iam,iamx,iamy,iamz, 'Initialization error: NR vs NLR ',& - & nr,nlr,mynx,myny,mynz - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - end if - - ! - ! Third example of use of CDALL: specify for each process - ! the set of global indices it owns. - ! - call psb_cdall(ctxt,desc_a,info,vl=myidx) - - case default - write(psb_err_unit,*) iam, 'Initialization error: should not get here' - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - return - end select - - - 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(ctxt) - 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 - - - ! loop over rows belonging to current process in a block - ! distribution. - - call psb_barrier(ctxt) - 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 - call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim) - ! x, y, z 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 - call ijk2idx(icol(icoeff),ix-1,iy,iz,idim,idim,idim) - 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 - call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim) - 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 - call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim) - irow(icoeff) = glob_row - icoeff = icoeff+1 - endif - - ! term depending on (x,y,z) - val(icoeff)=(2*sone)*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah & - & + c(x,y,z) - call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim) - 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 - call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim) - 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 - call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim) - 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 - call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim) - 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(:)=szero - 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(ctxt) - t1 = psb_wtime() - call psb_cdasb(desc_a,info,mold=imold) - tcdasb = psb_wtime()-t1 - call psb_barrier(ctxt) - 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(ctxt) - 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(ctxt) - ttot = psb_wtime() - t0 - - call psb_amx(ctxt,talc) - call psb_amx(ctxt,tgen) - call psb_amx(ctxt,tasb) - call psb_amx(ctxt,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(ctxt,err_act) - - return - end subroutine amg_s_gen_pde3d - -end module amg_s_pde3d_mod - program amg_s_pde3d use psb_base_mod use amg_prec_mod use psb_krylov_mod use psb_util_mod use data_input - use amg_s_pde3d_mod + use amg_s_pde3d_base_mod + use amg_s_pde3d_exp_mod + use amg_s_pde3d_gauss_mod + use amg_s_genpde_mod implicit none ! input parameters character(len=20) :: kmethd, ptype - character(len=5) :: afmt + character(len=5) :: afmt, pdecoeff integer(psb_ipk_) :: idim integer(psb_epk_) :: system_size - ! miscellaneous + ! miscellaneous real(psb_dpk_) :: t1, t2, tprec, thier, tslv ! sparse matrix and preconditioner @@ -691,7 +182,7 @@ program amg_s_pde3d call psb_init(ctxt) call psb_info(ctxt,iam,np) - if (iam < 0) then + if (iam < 0) then ! This should not happen, but just in case call psb_exit(ctxt) stop @@ -702,23 +193,40 @@ program amg_s_pde3d ! ! Hello world ! - if (iam == psb_root_) then - write(*,*) 'Welcome to MLD2P4 version: ',amg_version_string_ + if (iam == psb_root_) then + write(*,*) 'Welcome to AMG4PSBLAS version: ',amg_version_string_ write(*,*) 'This is the ',trim(name),' sample program' end if ! ! get parameters ! - call get_parms(ctxt,afmt,idim,s_choice,p_choice) + call get_parms(ctxt,afmt,idim,s_choice,p_choice,pdecoeff) ! - ! allocate and fill in the coefficient matrix, rhs and initial guess + ! allocate and fill in the coefficient matrix, rhs and initial guess ! call psb_barrier(ctxt) t1 = psb_wtime() - call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,info) + select case(psb_toupper(trim(pdecoeff))) + case("CONST") + call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,& + & a1,a2,a3,b1,b2,b3,c,g,info) + case("EXP") + call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,& + & a1_exp,a2_exp,a3_exp,b1_exp,b2_exp,b3_exp,c_exp,g_exp,info) + case("GAUSS") + call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,& + & a1_gauss,a2_gauss,a3_gauss,b1_gauss,b2_gauss,b3_gauss,c_gauss,g_gauss,info) + case default + info=psb_err_from_subroutine_ + ch_err='amg_gen_pdecoeff' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end select + + call psb_barrier(ctxt) t2 = psb_wtime() - t1 if(info /= psb_success_) then @@ -728,6 +236,8 @@ program amg_s_pde3d goto 9999 end if + if (iam == psb_root_) & + & write(psb_out_unit,'("PDE Coefficients : ",a)')pdecoeff if (iam == psb_root_) & & write(psb_out_unit,'("Overall matrix creation time : ",es12.5)')t2 if (iam == psb_root_) & @@ -743,7 +253,7 @@ program amg_s_pde3d case ('JACOBI','L1-JACOBI','GS','FWGS','FBGS') ! 1-level sweeps from "outer_sweeps" call prec%set('smoother_sweeps', p_choice%jsweeps, info) - + case ('BJAC') call prec%set('smoother_sweeps', p_choice%jsweeps, info) call prec%set('sub_solve', p_choice%solve, info) @@ -758,8 +268,8 @@ program amg_s_pde3d call prec%set('sub_solve', p_choice%solve, info) call prec%set('sub_fillin', p_choice%fill, info) call prec%set('sub_iluthrs', p_choice%thr, info) - - case ('ML') + + case ('ML') ! multilevel preconditioner call prec%set('ml_cycle', p_choice%mlcycle, info) @@ -824,7 +334,7 @@ program amg_s_pde3d call prec%set('coarse_sweeps', p_choice%cjswp, info) end select - + ! build the preconditioner call psb_barrier(ctxt) t1 = psb_wtime() @@ -854,7 +364,7 @@ program amg_s_pde3d end if ! - ! iterative method parameters + ! iterative method parameters ! call psb_barrier(ctxt) t1 = psb_wtime() @@ -894,9 +404,10 @@ program amg_s_pde3d call psb_sum(ctxt,descsize) call psb_sum(ctxt,precsize) call prec%descr(iout=psb_out_unit) - if (iam == psb_root_) then + if (iam == psb_root_) then write(psb_out_unit,'("Computed solution on ",i8," processors")') np write(psb_out_unit,'("Linear system size : ",i12)') system_size + write(psb_out_unit,'("PDE Coefficients : ",a)') trim(pdecoeff) write(psb_out_unit,'("Krylov method : ",a)') trim(s_choice%kmethd) write(psb_out_unit,'("Preconditioner : ",a)') trim(p_choice%descr) write(psb_out_unit,'("Iterations to convergence : ",i12)') iter @@ -918,7 +429,7 @@ program amg_s_pde3d end if - ! + ! ! cleanup storage and exit ! call psb_gefree(b,desc_a,info) @@ -945,7 +456,7 @@ contains ! ! get iteration parameters from standard input ! - subroutine get_parms(ctxt,afmt,idim,solve,prec) + subroutine get_parms(ctxt,afmt,idim,solve,prec,pdecoeff) implicit none @@ -954,6 +465,7 @@ contains character(len=*) :: afmt type(solverdata) :: solve type(precdata) :: prec + character(len=*) :: pdecoeff integer(psb_ipk_) :: iam, nm, np, inp_unit character(len=1024) :: filename @@ -978,6 +490,7 @@ contains ! call read_data(afmt,inp_unit) ! matrix storage format call read_data(idim,inp_unit) ! Discretization grid size + call read_data(pdecoeff,inp_unit) ! PDE Coefficients ! Krylov solver data call read_data(solve%kmethd,inp_unit) ! Krylov solver call read_data(solve%istopc,inp_unit) ! stopping criterion @@ -1039,6 +552,7 @@ contains call psb_bcast(ctxt,afmt) call psb_bcast(ctxt,idim) + call psb_bcast(ctxt,pdecoeff) call psb_bcast(ctxt,solve%kmethd) call psb_bcast(ctxt,solve%istopc) @@ -1051,7 +565,7 @@ contains call psb_bcast(ctxt,prec%ptype) ! broadcast first (pre-)smoother / 1-lev prec data - call psb_bcast(ctxt,prec%smther) + call psb_bcast(ctxt,prec%smther) call psb_bcast(ctxt,prec%jsweeps) call psb_bcast(ctxt,prec%novr) call psb_bcast(ctxt,prec%restr) @@ -1059,7 +573,7 @@ contains call psb_bcast(ctxt,prec%solve) call psb_bcast(ctxt,prec%fill) call psb_bcast(ctxt,prec%thr) - ! broadcast second (post-)smoother + ! broadcast second (post-)smoother call psb_bcast(ctxt,prec%smther2) call psb_bcast(ctxt,prec%jsweeps2) call psb_bcast(ctxt,prec%novr2) @@ -1068,12 +582,12 @@ contains call psb_bcast(ctxt,prec%solve2) call psb_bcast(ctxt,prec%fill2) call psb_bcast(ctxt,prec%thr2) - + ! broadcast AMG parameters call psb_bcast(ctxt,prec%mlcycle) call psb_bcast(ctxt,prec%outer_sweeps) call psb_bcast(ctxt,prec%maxlevs) - + call psb_bcast(ctxt,prec%aggr_prol) call psb_bcast(ctxt,prec%par_aggr_alg) call psb_bcast(ctxt,prec%aggr_ord) @@ -1085,7 +599,7 @@ contains call psb_bcast(ctxt,prec%athresv) end if call psb_bcast(ctxt,prec%athres) - + call psb_bcast(ctxt,prec%csize) call psb_bcast(ctxt,prec%cmat) call psb_bcast(ctxt,prec%csolve) diff --git a/tests/pdegen/amg_s_pde3d_base_mod.f90 b/tests/pdegen/amg_s_pde3d_base_mod.f90 new file mode 100644 index 00000000..a58a10ca --- /dev/null +++ b/tests/pdegen/amg_s_pde3d_base_mod.f90 @@ -0,0 +1,65 @@ +module amg_s_pde3d_base_mod + use psb_base_mod, only : psb_spk_, sone + real(psb_spk_), save, private :: epsilon=sone/80 +contains + subroutine pde_set_parm(dat) + real(psb_spk_), intent(in) :: dat + epsilon = dat + end subroutine pde_set_parm + ! + ! functions parametrizing the differential equation + ! + function b1(x,y,z) + use psb_base_mod, only : psb_spk_, sone + real(psb_spk_) :: b1 + real(psb_spk_), intent(in) :: x,y,z + b1=sone/sqrt(3.0_psb_spk_) + end function b1 + function b2(x,y,z) + use psb_base_mod, only : psb_spk_, sone + real(psb_spk_) :: b2 + real(psb_spk_), intent(in) :: x,y,z + b2=sone/sqrt(3.0_psb_spk_) + end function b2 + function b3(x,y,z) + use psb_base_mod, only : psb_spk_, sone + real(psb_spk_) :: b3 + real(psb_spk_), intent(in) :: x,y,z + b3=sone/sqrt(3.0_psb_spk_) + end function b3 + function c(x,y,z) + use psb_base_mod, only : psb_spk_, sone + real(psb_spk_) :: c + real(psb_spk_), intent(in) :: x,y,z + c=szero + 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=epsilon + 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=epsilon + 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=epsilon + end function a3 + function g(x,y,z) + use psb_base_mod, only : psb_spk_, sone, szero + 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 = sone + end if + end function g +end module amg_s_pde3d_base_mod diff --git a/tests/pdegen/amg_s_pde3d_exp_mod.f90 b/tests/pdegen/amg_s_pde3d_exp_mod.f90 new file mode 100644 index 00000000..ad8f48c7 --- /dev/null +++ b/tests/pdegen/amg_s_pde3d_exp_mod.f90 @@ -0,0 +1,65 @@ +module amg_s_pde3d_exp_mod + use psb_base_mod, only : psb_spk_, sone + real(psb_spk_), save, private :: epsilon=sone/160 +contains + subroutine pde_set_parm(dat) + real(psb_spk_), intent(in) :: dat + epsilon = dat + end subroutine pde_set_parm + ! + ! functions parametrizing the differential equation + ! + function b1_exp(x,y,z) + use psb_base_mod, only : psb_spk_, szero + real(psb_spk_) :: b1_exp + real(psb_spk_), intent(in) :: x,y,z + b1_exp=szero/sqrt(3.0_psb_spk_) + end function b1_exp + function b2_exp(x,y,z) + use psb_base_mod, only : psb_spk_, szero + real(psb_spk_) :: b2_exp + real(psb_spk_), intent(in) :: x,y,z + b2_exp=szero/sqrt(3.0_psb_spk_) + end function b2_exp + function b3_exp(x,y,z) + use psb_base_mod, only : psb_spk_, szero + real(psb_spk_) :: b3_exp + real(psb_spk_), intent(in) :: x,y,z + b3_exp=szero/sqrt(3.0_psb_spk_) + end function b3_exp + function c_exp(x,y,z) + use psb_base_mod, only : psb_spk_, szero + real(psb_spk_) :: c_exp + real(psb_spk_), intent(in) :: x,y,z + c_exp=szero + end function c_exp + function a1_exp(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: a1_exp + real(psb_spk_), intent(in) :: x,y,z + a1_exp=epsilon*exp(-(x+y+z)) + end function a1_exp + function a2_exp(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: a2_exp + real(psb_spk_), intent(in) :: x,y,z + a2_exp=epsilon*exp(-(x+y+z)) + end function a2_exp + function a3_exp(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: a3_exp + real(psb_spk_), intent(in) :: x,y,z + a3_exp=epsilon*exp(-(x+y+z)) + end function a3_exp + function g_exp(x,y,z) + use psb_base_mod, only : psb_spk_, sone, szero + real(psb_spk_) :: g_exp + real(psb_spk_), intent(in) :: x,y,z + g_exp = szero + if (x == sone) then + g_exp = sone + else if (x == szero) then + g_exp = sone + end if + end function g_exp +end module amg_s_pde3d_exp_mod diff --git a/tests/pdegen/amg_s_pde3d_gauss_mod.f90 b/tests/pdegen/amg_s_pde3d_gauss_mod.f90 new file mode 100644 index 00000000..5e62fa09 --- /dev/null +++ b/tests/pdegen/amg_s_pde3d_gauss_mod.f90 @@ -0,0 +1,65 @@ +module amg_s_pde3d_gauss_mod + use psb_base_mod, only : psb_spk_, sone + real(psb_spk_), save, private :: epsilon=sone/80 +contains + subroutine pde_set_parm(dat) + real(psb_spk_), intent(in) :: dat + epsilon = dat + end subroutine pde_set_parm + ! + ! functions parametrizing the differential equation + ! + function b1_gauss(x,y,z) + use psb_base_mod, only : psb_spk_, sone + real(psb_spk_) :: b1_gauss + real(psb_spk_), intent(in) :: x,y,z + b1_gauss=sone/sqrt(3.0_psb_spk_)-2*x*exp(-(x**2+y**2+z**2)) + end function b1_gauss + function b2_gauss(x,y,z) + use psb_base_mod, only : psb_spk_, sone + real(psb_spk_) :: b2_gauss + real(psb_spk_), intent(in) :: x,y,z + b2_gauss=sone/sqrt(3.0_psb_spk_)-2*y*exp(-(x**2+y**2+z**2)) + end function b2_gauss + function b3_gauss(x,y,z) + use psb_base_mod, only : psb_spk_, sone + real(psb_spk_) :: b3_gauss + real(psb_spk_), intent(in) :: x,y,z + b3_gauss=sone/sqrt(3.0_psb_spk_)-2*z*exp(-(x**2+y**2+z**2)) + end function b3_gauss + function c_gauss(x,y,z) + use psb_base_mod, only : psb_spk_, szero + real(psb_spk_) :: c_gauss + real(psb_spk_), intent(in) :: x,y,z + c=szero + end function c_gauss + function a1_gauss(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: a1_gauss + real(psb_spk_), intent(in) :: x,y,z + a1_gauss=epsilon*exp(-(x**2+y**2+z**2)) + end function a1_gauss + function a2_gauss(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: a2_gauss + real(psb_spk_), intent(in) :: x,y,z + a2_gauss=epsilon*exp(-(x**2+y**2+z**2)) + end function a2_gauss + function a3_gauss(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: a3_gauss + real(psb_spk_), intent(in) :: x,y,z + a3_gauss=epsilon*exp(-(x**2+y**2+z**2)) + end function a3_gauss + function g_gauss(x,y,z) + use psb_base_mod, only : psb_spk_, sone, szero + real(psb_spk_) :: g_gauss + real(psb_spk_), intent(in) :: x,y,z + g_gauss = szero + if (x == sone) then + g_gauss = sone + else if (x == szero) then + g_gauss = sone + end if + end function g_gauss +end module amg_s_pde3d_gauss_mod diff --git a/tests/pdegen/runs/amg_pde2d.inp b/tests/pdegen/runs/amg_pde2d.inp index 50f2e9ad..6c315155 100644 --- a/tests/pdegen/runs/amg_pde2d.inp +++ b/tests/pdegen/runs/amg_pde2d.inp @@ -1,20 +1,21 @@ %%%%%%%%%%% General arguments % Lines starting with % are ignored. -CSR ! Storage format CSR COO JAD +CSR ! Storage format CSR COO JAD 0200 ! IDIM; domain size. Linear system size is IDIM**2 +CONST ! PDECOEFF: CONST, EXP, BOX Coefficients of the PDE CG ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES 2 ! ISTOPC 00500 ! ITMAX 1 ! ITRACE -30 ! IRST (restart for RGMRES and BiCGSTABL) +30 ! IRST (restart for RGMRES and BiCGSTABL) 1.d-6 ! EPS %%%%%%%%%%% Main preconditioner choices %%%%%%%%%%%%%%%% -ML-VCYCLE-FBGS-R-UMF ! Longer descriptive name for preconditioner (up to 20 chars) +ML-VCYCLE-FBGS-R-UMF ! Longer descriptive name for preconditioner (up to 20 chars) ML ! Preconditioner type: NONE JACOBI GS FBGS BJAC AS ML %%%%%%%%%%% First smoother (for all levels but coarsest) %%%%%%%%%%%%%%%% -BJAC ! Smoother type JACOBI FBGS GS BWGS BJAC AS. For 1-level, repeats previous. +BJAC ! Smoother type JACOBI FBGS GS BWGS BJAC AS. For 1-level, repeats previous. 1 ! Number of sweeps for smoother 0 ! Number of overlap layers for AS preconditioner -HALO ! AS restriction operator: NONE HALO +HALO ! AS restriction operator: NONE HALO NONE ! AS prolongation operator: NONE SUM AVG ILU ! Subdomain solver for BJAC/AS: JACOBI GS BGS ILU ILUT MILU MUMPS SLU UMF 0 ! Fill level P for ILU(P) and ILU(T,P) @@ -23,11 +24,11 @@ ILU ! Subdomain solver for BJAC/AS: JACOBI GS BGS ILU IL NONE ! 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 +HALO ! AS restriction operator: NONE HALO NONE ! AS prolongation operator: NONE SUM AVG ILU ! Subdomain solver for BJAC/AS: JACOBI GS 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) +1.d-4 ! Threshold T for ILU(T,P) %%%%%%%%%%% Multilevel parameters %%%%%%%%%%%%%%%% VCYCLE ! Type of multilevel CYCLE: VCYCLE WCYCLE KCYCLE MULT ADD 4 ! Number of outer sweeps for ML @@ -39,12 +40,12 @@ NATURAL ! Ordering of aggregation NATURAL DEGREE FILTER ! 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 -0.05 0.025 ! Thresholds +0.05 0.025 ! Thresholds -0.0100d0 ! Smoothed aggregation threshold, ignored if < 0 %%%%%%%%%%% Coarse level solver %%%%%%%%%%%%%%%% -UMF ! Coarsest-level solver: MUMPS UMF SLU SLUDIST JACOBI GS BJAC +UMF ! Coarsest-level solver: MUMPS UMF SLU SLUDIST JACOBI GS BJAC UMF ! Coarsest-level subsolver for BJAC: ILU ILUT MILU UMF MUMPS SLU REPL ! 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) +1.d-4 ! Coarsest-level threshold T for ILU(T,P) 1 ! Number of sweeps for JACOBI/GS/BJAC coarsest-level solver diff --git a/tests/pdegen/runs/amg_pde3d.inp b/tests/pdegen/runs/amg_pde3d.inp index 7d1bcbb7..71ec48b1 100644 --- a/tests/pdegen/runs/amg_pde3d.inp +++ b/tests/pdegen/runs/amg_pde3d.inp @@ -1,19 +1,20 @@ %%%%%%%%%%% General arguments % Lines starting with % are ignored. -CSR ! Storage format CSR COO JAD +CSR ! Storage format CSR COO JAD 0080 ! IDIM; domain size. Linear system size is IDIM**3 +CONST ! PDECOEFF: CONST, EXP, GAUSS Coefficients of the PDE BICGSTAB ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES 2 ! ISTOPC 00500 ! ITMAX 1 ! ITRACE -30 ! IRST (restart for RGMRES and BiCGSTABL) +30 ! IRST (restart for RGMRES and BiCGSTABL) 1.d-6 ! EPS -ML-VCYCLE-FBGS-R-UMF ! Longer descriptive name for preconditioner (up to 20 chars) +ML-VCYCLE-FBGS-R-UMF ! Longer descriptive name for preconditioner (up to 20 chars) ML ! Preconditioner type: NONE JACOBI GS FBGS BJAC AS ML %%%%%%%%%%% First smoother (for all levels but coarsest) %%%%%%%%%%%%%%%% -FBGS ! Smoother type JACOBI FBGS GS BWGS BJAC AS. For 1-level, repeats previous. +FBGS ! Smoother type JACOBI FBGS GS BWGS BJAC AS. For 1-level, repeats previous. 1 ! Number of sweeps for smoother 0 ! Number of overlap layers for AS preconditioner -HALO ! AS restriction operator: NONE HALO +HALO ! AS restriction operator: NONE HALO NONE ! AS prolongation operator: NONE SUM AVG ILU ! Subdomain solver for BJAC/AS: JACOBI GS BGS ILU ILUT MILU MUMPS SLU UMF 0 ! Fill level P for ILU(P) and ILU(T,P) @@ -22,11 +23,11 @@ ILU ! Subdomain solver for BJAC/AS: JACOBI GS BGS ILU IL NONE ! 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 +HALO ! AS restriction operator: NONE HALO NONE ! AS prolongation operator: NONE SUM AVG ILU ! Subdomain solver for BJAC/AS: JACOBI GS 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) +1.d-4 ! Threshold T for ILU(T,P) %%%%%%%%%%% Multilevel parameters %%%%%%%%%%%%%%%% VCYCLE ! Type of multilevel CYCLE: VCYCLE WCYCLE KCYCLE MULT ADD 4 ! Number of outer sweeps for ML @@ -38,12 +39,12 @@ 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 -0.05 0.025 ! Thresholds +0.05 0.025 ! Thresholds -0.0100d0 ! Smoothed aggregation threshold, ignored if < 0 %%%%%%%%%%% Coarse level solver %%%%%%%%%%%%%%%% -UMF ! Coarsest-level solver: MUMPS UMF SLU SLUDIST JACOBI GS BJAC +UMF ! Coarsest-level solver: MUMPS UMF SLU SLUDIST JACOBI GS BJAC UMF ! Coarsest-level subsolver for BJAC: ILU ILUT MILU UMF MUMPS SLU -REPL ! Coarsest-level matrix distribution: DIST REPL +REPL ! 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) +1.d-4 ! Coarsest-level threshold T for ILU(T,P) 1 ! Number of sweeps for JACOBI/GS/BJAC coarsest-level solver