From 0aa1bd1c24d2aa78a2c3cf9db233c17dd97a9a61 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Fri, 17 Feb 2012 20:33:44 +0000 Subject: [PATCH] psblas3: base/modules/psb_sort_mod.f90 test/pargen/Makefile test/pargen/ppde.f90 test/pargen/ppde3d.f90 test/pargen/runs/ppde.inp test/pargen/spde.f90 test/pargen/spde3d.f90 util/Makefile util/psb_d_genmat_impl.f90 util/psb_d_genmat_mod.f90 util/psb_d_genpde_impl.f90 util/psb_d_genpde_mod.f90 util/psb_s_genpde_impl.f90 util/psb_s_genpde_mod.f90 util/psb_util_mod.f90 Factored PDE generation code. Defined 2D code. --- base/modules/psb_sort_mod.f90 | 2 +- test/pargen/Makefile | 27 +- test/pargen/{ppde.f90 => ppde3d.f90} | 344 +------------ test/pargen/runs/ppde.inp | 2 +- test/pargen/spde.f90 | 692 --------------------------- test/pargen/spde3d.f90 | 388 +++++++++++++++ util/Makefile | 4 +- util/psb_d_genmat_impl.f90 | 284 ----------- util/psb_d_genmat_mod.f90 | 41 -- util/psb_d_genpde_impl.f90 | 563 ++++++++++++++++++++++ util/psb_d_genpde_mod.f90 | 105 ++++ util/psb_s_genpde_impl.f90 | 563 ++++++++++++++++++++++ util/psb_s_genpde_mod.f90 | 105 ++++ util/psb_util_mod.f90 | 3 +- 14 files changed, 1772 insertions(+), 1351 deletions(-) rename test/pargen/{ppde.f90 => ppde3d.f90} (51%) delete mode 100644 test/pargen/spde.f90 create mode 100644 test/pargen/spde3d.f90 delete mode 100644 util/psb_d_genmat_impl.f90 delete mode 100644 util/psb_d_genmat_mod.f90 create mode 100644 util/psb_d_genpde_impl.f90 create mode 100644 util/psb_d_genpde_mod.f90 create mode 100644 util/psb_s_genpde_impl.f90 create mode 100644 util/psb_s_genpde_mod.f90 diff --git a/base/modules/psb_sort_mod.f90 b/base/modules/psb_sort_mod.f90 index 4d913cca..f50f12c1 100644 --- a/base/modules/psb_sort_mod.f90 +++ b/base/modules/psb_sort_mod.f90 @@ -114,7 +114,7 @@ module psb_sort_mod interface psb_isaperm logical function psb_isaperm(n,eip) import :: psb_ipk_ - integer(psb_ipk_), intent(in) :: n + integer(psb_ipk_), intent(in) :: n integer(psb_ipk_), intent(in) :: eip(n) end function psb_isaperm end interface diff --git a/test/pargen/Makefile b/test/pargen/Makefile index 6b2b027d..8e670ac7 100644 --- a/test/pargen/Makefile +++ b/test/pargen/Makefile @@ -14,20 +14,31 @@ FINCLUDES=$(FMFLAG)$(LIBDIR) $(FMFLAG). EXEDIR=./runs -all: ppde spde +all: ppde3d spde3d ppde2d spde2d -ppde: ppde.o - $(F90LINK) ppde.o -o ppde $(PSBLAS_LIB) $(LDLIBS) - /bin/mv ppde $(EXEDIR) +ppde3d: ppde3d.o + $(F90LINK) ppde3d.o -o ppde3d $(PSBLAS_LIB) $(LDLIBS) + /bin/mv ppde3d $(EXEDIR) -spde: spde.o - $(F90LINK) spde.o -o spde $(PSBLAS_LIB) $(LDLIBS) - /bin/mv spde $(EXEDIR) + +spde3d: spde3d.o + $(F90LINK) spde3d.o -o spde3d $(PSBLAS_LIB) $(LDLIBS) + /bin/mv spde3d $(EXEDIR) + +ppde2d: ppde2d.o + $(F90LINK) ppde2d.o -o ppde2d $(PSBLAS_LIB) $(LDLIBS) + /bin/mv ppde2d $(EXEDIR) + + +spde2d: spde2d.o + $(F90LINK) spde2d.o -o spde2d $(PSBLAS_LIB) $(LDLIBS) + /bin/mv spde2d $(EXEDIR) clean: - /bin/rm -f ppde.o spde.o $(EXEDIR)/ppde + /bin/rm -f ppde3d.o spde3d.o ppde2d.o spde2d.o \ + $(EXEDIR)/ppde3d $(EXEDIR)/spde3d $(EXEDIR)/ppde2d $(EXEDIR)/spde2d verycleanlib: (cd ../..; make veryclean) lib: diff --git a/test/pargen/ppde.f90 b/test/pargen/ppde3d.f90 similarity index 51% rename from test/pargen/ppde.f90 rename to test/pargen/ppde3d.f90 index 87ebbe0c..b710c527 100644 --- a/test/pargen/ppde.f90 +++ b/test/pargen/ppde3d.f90 @@ -29,39 +29,33 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -! File: ppde.f90 +! File: ppde3d.f90 ! -! Program: ppde +! Program: ppde3d ! This sample program solves a linear system obtained by discretizing a ! PDE with Dirichlet BCs. ! ! ! The PDE is a general second order equation in 3d ! -! b1 dd(u) b2 dd(u) b3 dd(u) a1 d(u) a2 d(u) a3 d(u) -! - ------ - ------ - ------ + ----- + ------ + ------ + a4 u = 0 -! dxdx dydy dzdz dx dy dz +! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u) +! - ------ - ------ - ------ + ----- + ------ + ------ + c u = f +! dxdx dydy dzdz dx dy dz ! -! with Dirichlet boundary conditions, on the unit cube 0<=x,y,z<=1. +! with Dirichlet boundary conditions +! u = g ! -! Example taken from: -! C.T.Kelley -! Iterative Methods for Linear and Nonlinear Equations -! SIAM 1995 +! on the unit cube 0<=x,y,z<=1. +! +! +! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation. ! ! In this sample program the index space of the discretized ! computational domain is first numbered sequentially in a standard way, ! then the corresponding vector is distributed according to a BLOCK ! data distribution. ! -! Boundary conditions are set in a very simple way, by adding -! equations of the form -! -! u(x,y) = exp(-x^2-y^2-z^2) -! -! Note that if a1=a2=a3=a4=0., the PDE is the well-known Laplace equation. -! -program ppde +program ppde3d use psb_base_mod use psb_prec_mod use psb_krylov_mod @@ -109,7 +103,7 @@ program ppde stop endif if(psb_get_errstatus() /= 0) goto 9999 - name='pde90' + name='pde3d90' call psb_set_errverbosity(2) ! ! Hello world @@ -128,14 +122,13 @@ program ppde ! call psb_barrier(ictxt) t1 = psb_wtime() - call gen_prob3d(ictxt,idim,a,bv,xxv,desc_a,afmt,& + call psb_gen_prob3d(ictxt,idim,a,bv,xxv,desc_a,afmt,& & a1,a2,a3,b1,b2,b3,c,g,info) -!!$ call create_matrix(idim,a,bv,xxv,desc_a,ictxt,afmt,info) call psb_barrier(ictxt) t2 = psb_wtime() - t1 if(info /= psb_success_) then info=psb_err_from_subroutine_ - ch_err='create_matrix' + ch_err='psb_gen_prob3d' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if @@ -283,7 +276,9 @@ contains call psb_bcast(ictxt,intbuf(1:5)) write(psb_out_unit,'("Solving matrix : ell1")') - write(psb_out_unit,'("Grid dimensions : ",i4,"x",i4,"x",i4)')idim,idim,idim + write(psb_out_unit,& + & '("Grid dimensions : ",i4," x ",i4," x ",i4)') & + & idim,idim,idim write(psb_out_unit,'("Number of processors : ",i0)')np write(psb_out_unit,'("Data distribution : BLOCK")') write(psb_out_unit,'("Preconditioner : ",a)') ptype @@ -315,7 +310,7 @@ contains subroutine pr_usage(iout) integer(psb_ipk_) :: iout write(iout,*)'incorrect parameter(s) found' - write(iout,*)' usage: pde90 methd prec dim & + write(iout,*)' usage: pde3d90 methd prec dim & &[istop itmax itrace]' write(iout,*)' where:' write(iout,*)' methd: cgstab cgs rgmres bicgstabl' @@ -330,299 +325,6 @@ contains write(iout,*)' iterations ' end subroutine pr_usage - ! - ! subroutine to allocate and fill in the coefficient matrix and - ! the rhs. - ! - subroutine create_matrix(idim,a,bv,xxv,desc_a,ictxt,afmt,info) - ! - ! discretize the partial diferential equation - ! - ! b1 dd(u) b2 dd(u) b3 dd(u) a1 d(u) a2 d(u) a3 d(u) - ! - ------ - ------ - ------ + ----- + ------ + ------ + a4 u - ! dxdx dydy dzdz dx dy dz - ! - ! with Dirichlet boundary conditions, on the unit cube 0<=x,y,z<=1. - ! - ! - ! Note that if a1=a2=a3=a4=0., the PDE is the well-known Laplace equation. - ! - use psb_base_mod - implicit none - integer(psb_ipk_) :: idim - integer(psb_ipk_), parameter :: nb=20 - type(psb_d_vect_type) :: xxv,bv - type(psb_desc_type) :: desc_a - integer(psb_ipk_) :: ictxt, info - character :: afmt*5 - type(psb_dspmat_type) :: a - type(psb_d_csc_sparse_mat) :: acsc - type(psb_d_coo_sparse_mat) :: acoo - type(psb_d_csr_sparse_mat) :: acsr - real(psb_dpk_) :: zt(nb),x,y,z - integer(psb_ipk_) :: m,n,nnz,glob_row,nlr,i,ii,ib,k - integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner - integer(psb_ipk_) :: np, iam, nr, nt - integer(psb_ipk_) :: element - integer(psb_ipk_), allocatable :: irow(:),icol(:),myidx(:) - real(psb_dpk_), allocatable :: val(:) - ! deltah dimension of each grid cell - ! deltat discretization time - real(psb_dpk_) :: deltah, sqdeltah, deltah2 - real(psb_dpk_), parameter :: rhs=0.d0,one=1.d0,zero=0.d0 - real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb - real(psb_dpk_) :: a1, a2, a3, a4, b1, b2, b3 - external :: a1, a2, a3, a4, b1, b2, b3 - integer(psb_ipk_) :: err_act - - character(len=20) :: name, ch_err,tmpfmt - - info = psb_success_ - name = 'create_matrix' - call psb_erractionsave(err_act) - - call psb_info(ictxt, iam, np) - - deltah = 1.d0/(idim+2) - sqdeltah = deltah*deltah - deltah2 = 2.d0* deltah - - ! initialize array descriptor and sparse matrix storage. provide an - ! estimate of the number of non zeroes - - m = idim*idim*idim - n = m - nnz = ((n*9)/(np)) - if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n - - ! - ! Using a simple BLOCK distribution. - ! - nt = (m+np-1)/np - nr = max(0,min(nt,m-(iam*nt))) - - nt = nr - call psb_sum(ictxt,nt) - if (nt /= m) write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m - call psb_barrier(ictxt) - t0 = psb_wtime() - call psb_cdall(ictxt,desc_a,info,nl=nr) - if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz) - ! define rhs from boundary conditions; also build initial guess - if (info == psb_success_) call psb_geall(xxv,desc_a,info) - if (info == psb_success_) call psb_geall(bv,desc_a,info) - nlr = desc_a%get_local_rows() - call psb_barrier(ictxt) - talc = psb_wtime()-t0 - - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='allocation rout.' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - ! we build an auxiliary matrix consisting of one row at a - ! time; just a small matrix. might be extended to generate - ! a bunch of rows per call. - ! - allocate(val(20*nb),irow(20*nb),& - &icol(20*nb),myidx(nlr),stat=info) - if (info /= psb_success_ ) then - info=psb_err_alloc_dealloc_ - call psb_errpush(info,name) - goto 9999 - endif - - do i=1,nlr - myidx(i) = i - end do - - - call psb_loc_to_glob(myidx,desc_a,info) - - ! loop over rows belonging to current process in a block - ! distribution. - - call psb_barrier(ictxt) - t1 = psb_wtime() - do ii=1, nlr,nb - ib = min(nb,nlr-ii+1) - element = 1 - do k=1,ib - i=ii+k-1 - ! local matrix pointer - glob_row=myidx(i) - ! compute gridpoint coordinates - if (mod(glob_row,(idim*idim)) == 0) then - ix = glob_row/(idim*idim) - else - ix = glob_row/(idim*idim)+1 - endif - if (mod((glob_row-(ix-1)*idim*idim),idim) == 0) then - iy = (glob_row-(ix-1)*idim*idim)/idim - else - iy = (glob_row-(ix-1)*idim*idim)/idim+1 - endif - iz = glob_row-(ix-1)*idim*idim-(iy-1)*idim - ! x, y, x coordinates - x = ix*deltah - y = iy*deltah - z = iz*deltah - if (glob_row == 1) then - write(0,*) 'Starting from ',ix,iy,iz,x,y,z,deltah - end if - if (glob_row == nt) then - write(0,*) 'Ending at ',ix,iy,iz,x,y,z,deltah - end if - if (i == nlr) then - write(0,*) 'Ending at ',ix,iy,iz,x,y,z,deltah - end if - ! check on boundary points - zt(k) = 0.d0 - ! internal point: build discretization - ! - ! term depending on (x-1,y,z) - ! - if (ix == 1) then - val(element) = -b1(x,y,z)/sqdeltah-a1(x,y,z)/deltah2 - zt(k) = exp(-x**2-y**2-z**2)*(-val(element)) - else - val(element) = -b1(x,y,z)/sqdeltah-a1(x,y,z)/deltah2 - icol(element) = (ix-2)*idim*idim+(iy-1)*idim+(iz) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x,y-1,z) - if (iy == 1) then - val(element) = -b2(x,y,z)/sqdeltah-a2(x,y,z)/deltah2 - zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element) = -b2(x,y,z)/sqdeltah-a2(x,y,z)/deltah2 - icol(element) = (ix-1)*idim*idim+(iy-2)*idim+(iz) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x,y,z-1) - if (iz == 1) then - val(element)=-b3(x,y,z)/sqdeltah-a3(x,y,z)/deltah2 - zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element)=-b3(x,y,z)/sqdeltah-a3(x,y,z)/deltah2 - icol(element) = (ix-1)*idim*idim+(iy-1)*idim+(iz-1) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x,y,z) - val(element)=(2*b1(x,y,z) + 2*b2(x,y,z) + 2*b3(x,y,z))/sqdeltah& - & +a4(x,y,z) - icol(element) = (ix-1)*idim*idim+(iy-1)*idim+(iz) - irow(element) = glob_row - element = element+1 - ! term depending on (x,y,z+1) - if (iz == idim) then - val(element)=-b3(x,y,z)/sqdeltah+a3(x,y,z)/deltah2 - zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element)=-b3(x,y,z)/sqdeltah+a3(x,y,z)/deltah2 - icol(element) = (ix-1)*idim*idim+(iy-1)*idim+(iz+1) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x,y+1,z) - if (iy == idim) then - val(element)=-b2(x,y,z)/sqdeltah+a2(x,y,z)/deltah2 - zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element)=-b2(x,y,z)/sqdeltah+a2(x,y,z)/deltah2 - icol(element) = (ix-1)*idim*idim+(iy)*idim+(iz) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x+1,y,z) - if (ix==idim) then - val(element)=-b1(x,y,z)/sqdeltah+a1(x,y,z)/deltah2 - zt(k) = exp(-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element)=-b1(x,y,z)/sqdeltah+a1(x,y,z)/deltah2 - icol(element) = (ix)*idim*idim+(iy-1)*idim+(iz) - irow(element) = glob_row - element = element+1 - endif - - end do - call psb_spins(element-1,irow,icol,val,a,desc_a,info) - if(info /= psb_success_) exit - call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info) - if(info /= psb_success_) exit - zt(:)=0.d0 - call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xxv,desc_a,info) - if(info /= psb_success_) exit - end do - - tgen = psb_wtime()-t1 - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='insert rout.' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - deallocate(val,irow,icol) - - call psb_barrier(ictxt) - t1 = psb_wtime() - call psb_cdasb(desc_a,info) - tcdasb = psb_wtime()-t1 - call psb_barrier(ictxt) - t1 = psb_wtime() - if (info == psb_success_) & - & call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) - call psb_barrier(ictxt) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='asb rout.' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - if (info == psb_success_) call psb_geasb(xxv,desc_a,info) - if (info == psb_success_) call psb_geasb(bv,desc_a,info) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='asb rout.' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - tasb = psb_wtime()-t1 - call psb_barrier(ictxt) - ttot = psb_wtime() - t0 - - call psb_amx(ictxt,talc) - call psb_amx(ictxt,tgen) - call psb_amx(ictxt,tasb) - call psb_amx(ictxt,ttot) - if(iam == psb_root_) then - tmpfmt = a%get_fmt() - write(psb_out_unit,'("The matrix has been generated and assembled in ",a3," format.")')& - & tmpfmt - write(psb_out_unit,'("-allocation time : ",es12.5)') talc - write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen - write(psb_out_unit,'("-desc asbly time : ",es12.5)') tcdasb - write(psb_out_unit,'("- mat asbly time : ",es12.5)') tasb - write(psb_out_unit,'("-total time : ",es12.5)') ttot - - end if - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error(ictxt) - return - end if - return - end subroutine create_matrix ! ! functions parametrizing the differential equation ! @@ -630,19 +332,19 @@ contains use psb_base_mod, only : psb_dpk_ real(psb_dpk_) :: b1 real(psb_dpk_), intent(in) :: x,y,z - b1=1.414d0 + b1=1.d0/sqrt(3.d0) end function b1 function b2(x,y,z) use psb_base_mod, only : psb_dpk_ real(psb_dpk_) :: b2 real(psb_dpk_), intent(in) :: x,y,z - b2=1.414d0 + b2=1.d0/sqrt(3.d0) end function b2 function b3(x,y,z) use psb_base_mod, only : psb_dpk_ real(psb_dpk_) :: b3 real(psb_dpk_), intent(in) :: x,y,z - b3=1.414d0 + b3=1.d0/sqrt(3.d0) end function b3 function c(x,y,z) use psb_base_mod, only : psb_dpk_ @@ -680,6 +382,6 @@ contains end if end function g -end program ppde +end program ppde3d diff --git a/test/pargen/runs/ppde.inp b/test/pargen/runs/ppde.inp index 60315ead..550d52c6 100644 --- a/test/pargen/runs/ppde.inp +++ b/test/pargen/runs/ppde.inp @@ -2,7 +2,7 @@ BICGSTAB Iterative method BICGSTAB CGS BICG BICGSTABL RGMRES BJAC Preconditioner NONE DIAG BJAC CSR Storage format for matrix A: CSR COO JAD -080 Domain size (acutal system is this**3) +060 Domain size (acutal system is this**3) 2 Stopping criterion 1000 MAXIT -2 ITRACE diff --git a/test/pargen/spde.f90 b/test/pargen/spde.f90 deleted file mode 100644 index 63963fa0..00000000 --- a/test/pargen/spde.f90 +++ /dev/null @@ -1,692 +0,0 @@ -!!$ -!!$ Parallel Sparse BLAS version 2.3.1 -!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010 -!!$ Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari CNRS-IRIT, Toulouse -!!$ -!!$ Redistribution and use in source and binary forms, with or without -!!$ modification, are permitted provided that the following conditions -!!$ are met: -!!$ 1. Redistributions of source code must retain the above copyright -!!$ notice, this list of conditions and the following disclaimer. -!!$ 2. Redistributions in binary form must reproduce the above copyright -!!$ notice, this list of conditions, and the following disclaimer in the -!!$ documentation and/or other materials provided with the distribution. -!!$ 3. The name of the PSBLAS 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 -!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS -!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -!!$ 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: ppde.f90 -! -! Program: ppde -! This sample program solves a linear system obtained by discretizing a -! PDE with Dirichlet BCs. -! -! -! The PDE is a general second order equation in 3d -! -! b1 dd(u) b2 dd(u) b3 dd(u) a1 d(u) a2 d(u) a3 d(u) -! - ------ - ------ - ------ + ----- + ------ + ------ + a4 u = 0 -! dxdx dydy dzdz dx dy dz -! -! with Dirichlet boundary conditions, on the unit cube 0<=x,y,z<=1. -! -! Example taken from: -! C.T.Kelley -! Iterative Methods for Linear and Nonlinear Equations -! SIAM 1995 -! -! In this sample program the index space of the discretized -! computational domain is first numbered sequentially in a standard way, -! then the corresponding vector is distributed according to a BLOCK -! data distribution. -! -! Boundary conditions are set in a very simple way, by adding -! equations of the form -! -! u(x,y) = exp(-x^2-y^2-z^2) -! -! Note that if a1=a2=a3=a4=0., the PDE is the well-known Laplace equation. -! -program ppde - use psb_base_mod - use psb_prec_mod - use psb_krylov_mod - use psb_util_mod - implicit none - - ! input parameters - character(len=20) :: kmethd, ptype - character(len=5) :: afmt - integer(psb_ipk_) :: idim - - ! miscellaneous - real(psb_spk_), parameter :: one = 1.0 - real(psb_dpk_) :: t1, t2, tprec - - ! sparse matrix and preconditioner - type(psb_sspmat_type) :: a - type(psb_sprec_type) :: prec - ! descriptor - type(psb_desc_type) :: desc_a, desc_b - ! dense matrices - type(psb_s_vect_type) :: xxv,bv, vtst - real(psb_spk_), allocatable :: tst(:) - ! blacs parameters - integer(psb_ipk_) :: ictxt, iam, np - - ! solver parameters - integer(psb_ipk_) :: iter, itmax,itrace, istopc, irst - integer(psb_long_int_k_) :: amatsize, precsize, descsize, d2size - real(psb_spk_) :: err, eps - - ! other variables - integer(psb_ipk_) :: info, i - character(len=20) :: name,ch_err - character(len=40) :: fname - - info=psb_success_ - - - call psb_init(ictxt) - call psb_info(ictxt,iam,np) - - if (iam < 0) then - ! This should not happen, but just in case - call psb_exit(ictxt) - stop - endif - if(psb_get_errstatus() /= 0) goto 9999 - name='pde90' - call psb_set_errverbosity(2) - ! - ! Hello world - ! - if (iam == psb_root_) then - write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_ - write(*,*) 'This is the ',trim(name),' sample program' - end if - ! - ! get parameters - ! - call get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst) - - ! - ! allocate and fill in the coefficient matrix, rhs and initial guess - ! - call psb_barrier(ictxt) - t1 = psb_wtime() - call create_matrix(idim,a,bv,xxv,desc_a,ictxt,afmt,info) - call psb_barrier(ictxt) - t2 = psb_wtime() - t1 - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='create_matrix' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - if (iam == psb_root_) write(psb_out_unit,'("Overall matrix creation time : ",es12.5)')t2 - if (iam == psb_root_) write(psb_out_unit,'(" ")') -!!$ write(fname,'(a,i0,a)') 'pde-',idim,'.hb' -!!$ call hb_write(a,info,filename=fname,rhs=b,key='PDEGEN',mtitle='MLD2P4 pdegen Test matrix ') -!!$ write(fname,'(a,i2.2,a,i2.2,a)') 'amat-',iam,'-',np,'.mtx' -!!$ call a%print(fname) -!!$ call psb_cdprt(20+iam,desc_a,short=.false.) -!!$ call psb_cdcpy(desc_a,desc_b,info) -!!$ call psb_set_debug_level(9999) - - call psb_cdbldext(a,desc_a,2,desc_b,info,extype=psb_ovt_asov_) - if (info /= 0) then - write(0,*) 'Error from bldext' - call psb_abort(ictxt) - end if - ! - ! prepare the preconditioner. - ! - if(iam == psb_root_) write(psb_out_unit,'("Setting preconditioner to : ",a)')ptype - call psb_precinit(prec,ptype,info) - - call psb_barrier(ictxt) - t1 = psb_wtime() - call psb_precbld(a,desc_a,prec,info) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_precbld' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - tprec = psb_wtime()-t1 - - call psb_amx(ictxt,tprec) - - if (iam == psb_root_) write(psb_out_unit,'("Preconditioner time : ",es12.5)')tprec - if (iam == psb_root_) write(psb_out_unit,'(" ")') - ! - ! iterative method parameters - ! - if(iam == psb_root_) write(psb_out_unit,'("Calling iterative method ",a)')kmethd - call psb_barrier(ictxt) - t1 = psb_wtime() - eps = 1.d-9 - call psb_krylov(kmethd,a,prec,bv,xxv,eps,desc_a,info,& - & itmax=itmax,iter=iter,err=err,itrace=itrace,istop=istopc,irst=irst) - - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='solver routine' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - call psb_barrier(ictxt) - t2 = psb_wtime() - t1 - call psb_amx(ictxt,t2) - amatsize = a%sizeof() - descsize = desc_a%sizeof() - precsize = prec%sizeof() - call psb_sum(ictxt,amatsize) - call psb_sum(ictxt,descsize) - call psb_sum(ictxt,precsize) - - if (iam == psb_root_) then - write(psb_out_unit,'(" ")') - write(psb_out_unit,'("Time to solve matrix : ",es12.5)')t2 - write(psb_out_unit,'("Time per iteration : ",es12.5)')t2/iter - write(psb_out_unit,'("Number of iterations : ",i0)')iter - write(psb_out_unit,'("Convergence indicator on exit : ",es12.5)')err - write(psb_out_unit,'("Info on exit : ",i0)')info - write(psb_out_unit,'("Total memory occupation for A: ",i12)')amatsize - write(psb_out_unit,'("Total memory occupation for PREC: ",i12)')precsize - write(psb_out_unit,'("Total memory occupation for DESC_A: ",i12)')descsize - write(psb_out_unit,'("Storage type for DESC_A: ",a)') desc_a%indxmap%get_fmt() - write(psb_out_unit,'("Storage type for DESC_B: ",a)') desc_b%indxmap%get_fmt() - end if - - ! - if (.false.) then - call psb_geall(tst,desc_b, info) - call psb_geall(vtst,desc_b, info) - vtst%v%v = iam+1 - call psb_geasb(vtst,desc_b,info) - tst = vtst%get_vect() - call psb_geasb(tst,desc_b,info) - call psb_ovrl(vtst,desc_b,info,update=psb_avg_) - call psb_ovrl(tst,desc_b,info,update=psb_avg_) - write(0,*) iam,' After ovrl:',vtst%v%v - write(0,*) iam,' After ovrl:',tst - end if - - ! - ! cleanup storage and exit - ! - call psb_gefree(bv,desc_a,info) - call psb_gefree(xxv,desc_a,info) - call psb_spfree(a,desc_a,info) - call psb_precfree(prec,info) - call psb_cdfree(desc_a,info) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='free routine' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - -9999 continue - if(info /= psb_success_) then - call psb_error(ictxt) - end if - call psb_exit(ictxt) - stop - -contains - ! - ! get iteration parameters from standard input - ! - subroutine get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst) - integer(psb_ipk_) :: ictxt - character(len=*) :: kmethd, ptype, afmt - integer(psb_ipk_) :: idim, istopc,itmax,itrace,irst - integer(psb_ipk_) :: np, iam - integer(psb_ipk_) :: intbuf(10), ip - - call psb_info(ictxt, iam, np) - - if (iam == 0) then - read(psb_inp_unit,*) ip - if (ip >= 3) then - read(psb_inp_unit,*) kmethd - read(psb_inp_unit,*) ptype - read(psb_inp_unit,*) afmt - - ! broadcast parameters to all processors - call psb_bcast(ictxt,kmethd) - call psb_bcast(ictxt,afmt) - call psb_bcast(ictxt,ptype) - - - read(psb_inp_unit,*) idim - if (ip >= 4) then - read(psb_inp_unit,*) istopc - else - istopc=1 - endif - if (ip >= 5) then - read(psb_inp_unit,*) itmax - else - itmax=500 - endif - if (ip >= 6) then - read(psb_inp_unit,*) itrace - else - itrace=-1 - endif - if (ip >= 7) then - read(psb_inp_unit,*) irst - else - irst=1 - endif - ! broadcast parameters to all processors - - intbuf(1) = idim - intbuf(2) = istopc - intbuf(3) = itmax - intbuf(4) = itrace - intbuf(5) = irst - call psb_bcast(ictxt,intbuf(1:5)) - - write(psb_out_unit,'("Solving matrix : ell1")') - write(psb_out_unit,'("Grid dimensions : ",i4,"x",i4,"x",i4)')idim,idim,idim - write(psb_out_unit,'("Number of processors : ",i0)')np - write(psb_out_unit,'("Data distribution : BLOCK")') - write(psb_out_unit,'("Preconditioner : ",a)') ptype - write(psb_out_unit,'("Iterative method : ",a)') kmethd - write(psb_out_unit,'(" ")') - else - ! wrong number of parameter, print an error message and exit - call pr_usage(0) - call psb_abort(ictxt) - stop 1 - endif - else - call psb_bcast(ictxt,kmethd) - call psb_bcast(ictxt,afmt) - call psb_bcast(ictxt,ptype) - call psb_bcast(ictxt,intbuf(1:5)) - idim = intbuf(1) - istopc = intbuf(2) - itmax = intbuf(3) - itrace = intbuf(4) - irst = intbuf(5) - end if - return - - end subroutine get_parms - ! - ! print an error message - ! - subroutine pr_usage(iout) - integer(psb_ipk_) :: iout - write(iout,*)'incorrect parameter(s) found' - write(iout,*)' usage: pde90 methd prec dim & - &[istop itmax itrace]' - write(iout,*)' where:' - write(iout,*)' methd: cgstab cgs rgmres bicgstabl' - write(iout,*)' prec : bjac diag none' - write(iout,*)' dim number of points along each axis' - write(iout,*)' the size of the resulting linear ' - write(iout,*)' system is dim**3' - write(iout,*)' istop stopping criterion 1, 2 ' - write(iout,*)' itmax maximum number of iterations [500] ' - write(iout,*)' itrace <=0 (no tracing, default) or ' - write(iout,*)' >= 1 do tracing every itrace' - write(iout,*)' iterations ' - end subroutine pr_usage - - ! - ! subroutine to allocate and fill in the coefficient matrix and - ! the rhs. - ! - subroutine create_matrix(idim,a,bv,xxv,desc_a,ictxt,afmt,info) - ! - ! discretize the partial diferential equation - ! - ! b1 dd(u) b2 dd(u) b3 dd(u) a1 d(u) a2 d(u) a3 d(u) - ! - ------ - ------ - ------ + ----- + ------ + ------ + a4 u - ! dxdx dydy dzdz dx dy dz - ! - ! with Dirichlet boundary conditions, on the unit cube 0<=x,y,z<=1. - ! - ! Boundary conditions are set in a very simple way, by adding - ! equations of the form - ! - ! u(x,y) = exp(-x^2-y^2-z^2) - ! - ! Note that if a1=a2=a3=a4=0., the PDE is the well-known Laplace equation. - ! - use psb_base_mod - use psb_mat_mod - implicit none - integer(psb_ipk_) :: idim - integer(psb_ipk_), parameter :: nb=20 - type(psb_s_vect_type) :: xxv,bv - type(psb_desc_type) :: desc_a - integer(psb_ipk_) :: ictxt, info - character :: afmt*5 - type(psb_sspmat_type) :: a - type(psb_s_csc_sparse_mat) :: acsc - type(psb_s_coo_sparse_mat) :: acoo - type(psb_s_csr_sparse_mat) :: acsr - real(psb_spk_) :: zt(nb),x,y,z - integer(psb_ipk_) :: m,n,nnz,glob_row,nlr,i,ii,ib,k - integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner - integer(psb_ipk_) :: np, iam, nr, nt - integer(psb_ipk_) :: element - integer(psb_ipk_), allocatable :: irow(:),icol(:),myidx(:) - real(psb_spk_), allocatable :: val(:) - ! deltah dimension of each grid cell - ! deltat discretization time - real(psb_spk_) :: deltah, sqdeltah, deltah2 - real(psb_spk_),parameter :: rhs=0.0,one=1.0,zero=0.0 - real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen - real(psb_spk_) :: a1, a2, a3, a4, b1, b2, b3 - external :: a1, a2, a3, a4, b1, b2, b3 - integer(psb_ipk_) :: err_act - - character(len=20) :: name, ch_err,tmpfmt - - info = psb_success_ - name = 'create_matrix' - call psb_erractionsave(err_act) - - call psb_info(ictxt, iam, np) - - deltah = 1.0/(idim-1) - sqdeltah = deltah*deltah - deltah2 = 2.0* deltah - - ! initialize array descriptor and sparse matrix storage. provide an - ! estimate of the number of non zeroes - - m = idim*idim*idim - n = m - nnz = ((n*9)/(np)) - if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n - - ! - ! Using a simple BLOCK distribution. - ! - nt = (m+np-1)/np - nr = max(0,min(nt,m-(iam*nt))) - - nt = nr - call psb_sum(ictxt,nt) - if (nt /= m) write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m - call psb_barrier(ictxt) - t0 = psb_wtime() - call psb_cdall(ictxt,desc_a,info,nl=nr) - if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz) - ! define rhs from boundary conditions; also build initial guess - if (info == psb_success_) call psb_geall(xxv,desc_a,info) - if (info == psb_success_) call psb_geall(bv,desc_a,info) - nlr = desc_a%get_local_rows() - call psb_barrier(ictxt) - talc = psb_wtime()-t0 - - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='allocation rout.' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - ! we build an auxiliary matrix consisting of one row at a - ! time; just a small matrix. might be extended to generate - ! a bunch of rows per call. - ! - allocate(val(20*nb),irow(20*nb),& - &icol(20*nb),myidx(nlr),stat=info) - if (info /= psb_success_ ) then - info=psb_err_alloc_dealloc_ - call psb_errpush(info,name) - goto 9999 - endif - - do i=1,nlr - myidx(i) = i - end do - - - call psb_loc_to_glob(myidx,desc_a,info) - - ! loop over rows belonging to current process in a block - ! distribution. - - call psb_barrier(ictxt) - t1 = psb_wtime() - do ii=1, nlr,nb - ib = min(nb,nlr-ii+1) - element = 1 - do k=1,ib - i=ii+k-1 - ! local matrix pointer - glob_row=myidx(i) - ! compute gridpoint coordinates - if (mod(glob_row,(idim*idim)) == 0) then - ix = glob_row/(idim*idim) - else - ix = glob_row/(idim*idim)+1 - endif - if (mod((glob_row-(ix-1)*idim*idim),idim) == 0) then - iy = (glob_row-(ix-1)*idim*idim)/idim - else - iy = (glob_row-(ix-1)*idim*idim)/idim+1 - endif - iz = glob_row-(ix-1)*idim*idim-(iy-1)*idim - ! x, y, x coordinates - x = ix*deltah - y = iy*deltah - z = iz*deltah - - ! check on boundary points - zt(k) = 0.d0 - ! internal point: build discretization - ! - ! term depending on (x-1,y,z) - ! - if (ix == 1) then - val(element) = -b1(x,y,z)/sqdeltah-a1(x,y,z)/deltah2 - zt(k) = exp(-x**2-y**2-z**2)*(-val(element)) - else - val(element) = -b1(x,y,z)/sqdeltah-a1(x,y,z)/deltah2 - icol(element) = (ix-2)*idim*idim+(iy-1)*idim+(iz) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x,y-1,z) - if (iy == 1) then - val(element) = -b2(x,y,z)/sqdeltah-a2(x,y,z)/deltah2 - zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element) = -b2(x,y,z)/sqdeltah-a2(x,y,z)/deltah2 - icol(element) = (ix-1)*idim*idim+(iy-2)*idim+(iz) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x,y,z-1) - if (iz == 1) then - val(element)=-b3(x,y,z)/sqdeltah-a3(x,y,z)/deltah2 - zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element)=-b3(x,y,z)/sqdeltah-a3(x,y,z)/deltah2 - icol(element) = (ix-1)*idim*idim+(iy-1)*idim+(iz-1) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x,y,z) - val(element)=(2*b1(x,y,z) + 2*b2(x,y,z) + 2*b3(x,y,z))/sqdeltah& - & +a4(x,y,z) - icol(element) = (ix-1)*idim*idim+(iy-1)*idim+(iz) - irow(element) = glob_row - element = element+1 - ! term depending on (x,y,z+1) - if (iz == idim) then - val(element)=-b3(x,y,z)/sqdeltah+a3(x,y,z)/deltah2 - zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element)=-b3(x,y,z)/sqdeltah+a3(x,y,z)/deltah2 - icol(element) = (ix-1)*idim*idim+(iy-1)*idim+(iz+1) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x,y+1,z) - if (iy == idim) then - val(element)=-b2(x,y,z)/sqdeltah+a2(x,y,z)/deltah2 - zt(k) = exp(-x**2-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element)=-b2(x,y,z)/sqdeltah+a2(x,y,z)/deltah2 - icol(element) = (ix-1)*idim*idim+(iy)*idim+(iz) - irow(element) = glob_row - element = element+1 - endif - ! term depending on (x+1,y,z) - if (ix==idim) then - val(element)=-b1(x,y,z)/sqdeltah+a1(x,y,z)/deltah2 - zt(k) = exp(-y**2-z**2)*exp(-x)*(-val(element)) - else - val(element)=-b1(x,y,z)/sqdeltah+a1(x,y,z)/deltah2 - icol(element) = (ix)*idim*idim+(iy-1)*idim+(iz) - irow(element) = glob_row - element = element+1 - endif - - end do - call psb_spins(element-1,irow,icol,val,a,desc_a,info) - if(info /= psb_success_) exit - call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info) - if(info /= psb_success_) exit - zt(:)=0.d0 - call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xxv,desc_a,info) - if(info /= psb_success_) exit - end do - - tgen = psb_wtime()-t1 - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='insert rout.' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - deallocate(val,irow,icol) - - call psb_barrier(ictxt) - t1 = psb_wtime() - call psb_cdasb(desc_a,info) - if (info == psb_success_) & - & call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) - call psb_barrier(ictxt) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='asb rout.' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - if (info == psb_success_) call psb_geasb(xxv,desc_a,info) - if (info == psb_success_) call psb_geasb(bv,desc_a,info) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='asb rout.' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - tasb = psb_wtime()-t1 - call psb_barrier(ictxt) - ttot = psb_wtime() - t0 - - call psb_amx(ictxt,talc) - call psb_amx(ictxt,tgen) - call psb_amx(ictxt,tasb) - call psb_amx(ictxt,ttot) - if(iam == psb_root_) then - tmpfmt = a%get_fmt() - write(psb_out_unit,'("The matrix has been generated and assembled in ",a3," format.")')& - & tmpfmt - write(psb_out_unit,'("-allocation time : ",es12.5)') talc - write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen - write(psb_out_unit,'("-assembly time : ",es12.5)') tasb - write(psb_out_unit,'("-total time : ",es12.5)') ttot - - end if - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error(ictxt) - return - end if - return - end subroutine create_matrix -end program ppde -! -! functions parametrizing the differential equation -! -function a1(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: a1 - real(psb_spk_) :: x,y,z - a1=1.414e0 -end function a1 -function a2(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: a2 - real(psb_spk_) :: x,y,z - a2=1.414e0 -end function a2 -function a3(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: a3 - real(psb_spk_) :: x,y,z - a3=1.414e0 -end function a3 -function a4(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: a4 - real(psb_spk_) :: x,y,z - a4=0.e0 -end function a4 -function b1(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: b1 - real(psb_spk_) :: x,y,z - b1=1.e0/80 -end function b1 -function b2(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: b2 - real(psb_spk_) :: x,y,z - b2=1.e0/80 -end function b2 -function b3(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: b3 - real(psb_spk_) :: x,y,z - b3=1.e0/80 -end function b3 - - diff --git a/test/pargen/spde3d.f90 b/test/pargen/spde3d.f90 new file mode 100644 index 00000000..b3673f22 --- /dev/null +++ b/test/pargen/spde3d.f90 @@ -0,0 +1,388 @@ +!!$ +!!$ Parallel Sparse BLAS version 2.3.1 +!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the PSBLAS 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 +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ 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: spde3d.f90 +! +! Program: spde3d +! This sample program solves a linear system obtained by discretizing a +! 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) +! - ------ - ------ - ------ + ----- + ------ + ------ + 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. +! +! In this sample program the index space of the discretized +! computational domain is first numbered sequentially in a standard way, +! then the corresponding vector is distributed according to a BLOCK +! data distribution. +! +! +program spde3d + use psb_base_mod + use psb_prec_mod + use psb_krylov_mod + use psb_util_mod + implicit none + + ! input parameters + character(len=20) :: kmethd, ptype + character(len=5) :: afmt + integer(psb_ipk_) :: idim + + ! miscellaneous + real(psb_spk_), parameter :: one = 1.e0 + real(psb_dpk_) :: t1, t2, tprec + + ! sparse matrix and preconditioner + type(psb_sspmat_type) :: a + type(psb_sprec_type) :: prec + ! descriptor + type(psb_desc_type) :: desc_a + ! dense vectors + type(psb_s_vect_type) :: xxv,bv + ! parallel environment + integer(psb_ipk_) :: ictxt, iam, np + + ! solver parameters + integer(psb_ipk_) :: iter, itmax,itrace, istopc, irst + integer(psb_long_int_k_) :: amatsize, precsize, descsize, d2size + real(psb_spk_) :: err, eps + + ! other variables + integer(psb_ipk_) :: info, i + character(len=20) :: name,ch_err + character(len=40) :: fname + + info=psb_success_ + + + call psb_init(ictxt) + call psb_info(ictxt,iam,np) + + if (iam < 0) then + ! This should not happen, but just in case + call psb_exit(ictxt) + stop + endif + if(psb_get_errstatus() /= 0) goto 9999 + name='pde3d90' + call psb_set_errverbosity(2) + ! + ! Hello world + ! + if (iam == psb_root_) then + write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_ + write(*,*) 'This is the ',trim(name),' sample program' + end if + ! + ! get parameters + ! + call get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst) + + ! + ! allocate and fill in the coefficient matrix, rhs and initial guess + ! + call psb_barrier(ictxt) + t1 = psb_wtime() + call psb_gen_prob3d(ictxt,idim,a,bv,xxv,desc_a,afmt,& + & a1,a2,a3,b1,b2,b3,c,g,info) + call psb_barrier(ictxt) + t2 = psb_wtime() - t1 + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_gen_prob3d' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + if (iam == psb_root_) write(psb_out_unit,'("Overall matrix creation time : ",es12.5)')t2 + if (iam == psb_root_) write(psb_out_unit,'(" ")') + ! + ! prepare the preconditioner. + ! + if(iam == psb_root_) write(psb_out_unit,'("Setting preconditioner to : ",a)')ptype + call psb_precinit(prec,ptype,info) + + call psb_barrier(ictxt) + t1 = psb_wtime() + call psb_precbld(a,desc_a,prec,info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_precbld' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + tprec = psb_wtime()-t1 + + call psb_amx(ictxt,tprec) + + if (iam == psb_root_) write(psb_out_unit,'("Preconditioner time : ",es12.5)')tprec + if (iam == psb_root_) write(psb_out_unit,'(" ")') + ! + ! iterative method parameters + ! + if(iam == psb_root_) write(psb_out_unit,'("Calling iterative method ",a)')kmethd + call psb_barrier(ictxt) + t1 = psb_wtime() + eps = 1.d-9 + call psb_krylov(kmethd,a,prec,bv,xxv,eps,desc_a,info,& + & itmax=itmax,iter=iter,err=err,itrace=itrace,istop=istopc,irst=irst) + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='solver routine' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + call psb_barrier(ictxt) + t2 = psb_wtime() - t1 + call psb_amx(ictxt,t2) + amatsize = a%sizeof() + descsize = desc_a%sizeof() + precsize = prec%sizeof() + call psb_sum(ictxt,amatsize) + call psb_sum(ictxt,descsize) + call psb_sum(ictxt,precsize) + + if (iam == psb_root_) then + write(psb_out_unit,'(" ")') + write(psb_out_unit,'("Time to solve matrix : ",es12.5)')t2 + write(psb_out_unit,'("Time per iteration : ",es12.5)')t2/iter + write(psb_out_unit,'("Number of iterations : ",i0)')iter + write(psb_out_unit,'("Convergence indicator on exit : ",es12.5)')err + write(psb_out_unit,'("Info on exit : ",i0)')info + write(psb_out_unit,'("Total memory occupation for A: ",i12)')amatsize + write(psb_out_unit,'("Total memory occupation for PREC: ",i12)')precsize + write(psb_out_unit,'("Total memory occupation for DESC_A: ",i12)')descsize + write(psb_out_unit,'("Storage type for DESC_A: ",a)') desc_a%indxmap%get_fmt() + end if + + + ! + ! cleanup storage and exit + ! + call psb_gefree(bv,desc_a,info) + call psb_gefree(xxv,desc_a,info) + call psb_spfree(a,desc_a,info) + call psb_precfree(prec,info) + call psb_cdfree(desc_a,info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='free routine' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + +9999 continue + if(info /= psb_success_) then + call psb_error(ictxt) + end if + call psb_exit(ictxt) + stop + +contains + ! + ! get iteration parameters from standard input + ! + subroutine get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst) + integer(psb_ipk_) :: ictxt + character(len=*) :: kmethd, ptype, afmt + integer(psb_ipk_) :: idim, istopc,itmax,itrace,irst + integer(psb_ipk_) :: np, iam + integer(psb_ipk_) :: intbuf(10), ip + + call psb_info(ictxt, iam, np) + + if (iam == 0) then + read(psb_inp_unit,*) ip + if (ip >= 3) then + read(psb_inp_unit,*) kmethd + read(psb_inp_unit,*) ptype + read(psb_inp_unit,*) afmt + + ! broadcast parameters to all processors + call psb_bcast(ictxt,kmethd) + call psb_bcast(ictxt,afmt) + call psb_bcast(ictxt,ptype) + + + read(psb_inp_unit,*) idim + if (ip >= 4) then + read(psb_inp_unit,*) istopc + else + istopc=1 + endif + if (ip >= 5) then + read(psb_inp_unit,*) itmax + else + itmax=500 + endif + if (ip >= 6) then + read(psb_inp_unit,*) itrace + else + itrace=-1 + endif + if (ip >= 7) then + read(psb_inp_unit,*) irst + else + irst=1 + endif + ! broadcast parameters to all processors + + intbuf(1) = idim + intbuf(2) = istopc + intbuf(3) = itmax + intbuf(4) = itrace + intbuf(5) = irst + call psb_bcast(ictxt,intbuf(1:5)) + + write(psb_out_unit,'("Solving matrix : ell1")') + write(psb_out_unit,& + & '("Grid dimensions : ",i4," x ",i4," x ",i4)') & + & idim,idim,idim + write(psb_out_unit,'("Number of processors : ",i0)')np + write(psb_out_unit,'("Data distribution : BLOCK")') + write(psb_out_unit,'("Preconditioner : ",a)') ptype + write(psb_out_unit,'("Iterative method : ",a)') kmethd + write(psb_out_unit,'(" ")') + else + ! wrong number of parameter, print an error message and exit + call pr_usage(0) + call psb_abort(ictxt) + stop 1 + endif + else + call psb_bcast(ictxt,kmethd) + call psb_bcast(ictxt,afmt) + call psb_bcast(ictxt,ptype) + call psb_bcast(ictxt,intbuf(1:5)) + idim = intbuf(1) + istopc = intbuf(2) + itmax = intbuf(3) + itrace = intbuf(4) + irst = intbuf(5) + end if + return + + end subroutine get_parms + ! + ! print an error message + ! + subroutine pr_usage(iout) + integer(psb_ipk_) :: iout + write(iout,*)'incorrect parameter(s) found' + write(iout,*)' usage: pde3d90 methd prec dim & + &[istop itmax itrace]' + write(iout,*)' where:' + write(iout,*)' methd: cgstab cgs rgmres bicgstabl' + write(iout,*)' prec : bjac diag none' + write(iout,*)' dim number of points along each axis' + write(iout,*)' the size of the resulting linear ' + write(iout,*)' system is dim**3' + write(iout,*)' istop stopping criterion 1, 2 ' + write(iout,*)' itmax maximum number of iterations [500] ' + write(iout,*)' itrace <=0 (no tracing, default) or ' + write(iout,*)' >= 1 do tracing every itrace' + write(iout,*)' iterations ' + end subroutine pr_usage + + ! + ! functions parametrizing the differential equation + ! + function b1(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: b1 + real(psb_spk_), intent(in) :: x,y,z + b1=1.e0/sqrt(3.e0) + end function b1 + function b2(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: b2 + real(psb_spk_), intent(in) :: x,y,z + b2=1.e0/sqrt(3.e0) + end function b2 + function b3(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: b3 + real(psb_spk_), intent(in) :: x,y,z + b3=1.e0/sqrt(3.e0) + end function b3 + function c(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: c + real(psb_spk_), intent(in) :: x,y,z + c=0.e0 + end function c + function a1(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: a1 + real(psb_spk_), intent(in) :: x,y,z + a1=1.e0/80 + end function a1 + function a2(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: a2 + real(psb_spk_), intent(in) :: x,y,z + a2=1.e0/80 + end function a2 + function a3(x,y,z) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: a3 + real(psb_spk_), intent(in) :: x,y,z + a3=1.e0/80 + end function a3 + function g(x,y,z) + use psb_base_mod, only : psb_spk_, sone + real(psb_spk_) :: g + real(psb_spk_), intent(in) :: x,y,z + g = szero + if (x == sone) then + g = sone + else if (x == szero) then + g = exp(y**2-z**2) + end if + end function g + +end program spde3d + + diff --git a/util/Makefile b/util/Makefile index 864ece2c..f5405ef1 100644 --- a/util/Makefile +++ b/util/Makefile @@ -7,7 +7,7 @@ HERE=. BASEOBJS= psb_blockpart_mod.o psb_metispart_mod.o \ psb_hbio_mod.o psb_mmio_mod.o psb_mat_dist_mod.o \ - psb_renum_mod.o psb_gps_mod.o psb_d_genmat_mod.o + psb_renum_mod.o psb_gps_mod.o psb_d_genpde_mod.o psb_s_genpde_mod.o IMPLOBJS= psb_s_hbio_impl.o psb_d_hbio_impl.o \ psb_c_hbio_impl.o psb_z_hbio_impl.o \ psb_s_mmio_impl.o psb_d_mmio_impl.o \ @@ -16,7 +16,7 @@ IMPLOBJS= psb_s_hbio_impl.o psb_d_hbio_impl.o \ psb_c_mat_dist_impl.o psb_z_mat_dist_impl.o \ psb_s_renum_impl.o psb_d_renum_impl.o \ psb_c_renum_impl.o psb_z_renum_impl.o \ - psb_d_genmat_impl.o + psb_d_genpde_impl.o psb_s_genpde_impl.o MODOBJS=psb_util_mod.o $(BASEOBJS) COBJS=psb_amd_order.o OBJS=$(MODOBJS) $(IMPLOBJS) $(COBJS) diff --git a/util/psb_d_genmat_impl.f90 b/util/psb_d_genmat_impl.f90 deleted file mode 100644 index 21322bf5..00000000 --- a/util/psb_d_genmat_impl.f90 +++ /dev/null @@ -1,284 +0,0 @@ -! -! subroutine to allocate and fill in the coefficient matrix and -! the rhs. -! -subroutine gen_prob3d(ictxt,idim,a,bv,xv,desc_a,afmt,a1,a2,a3,b1,b2,b3,c,g,info) - use psb_base_mod - use psb_d_genmat_mod, psb_protect_name => gen_prob3d - ! - ! Discretizes the partial differential equation - ! - ! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u) - ! - ------ - ------ - ------ + ----- + ------ + ------ + c u = 0 - ! dxdx dydy dzdz dx dy dz - ! - ! with Dirichlet boundary conditions - ! u = g - ! - ! on the unit cube 0<=x,y,z<=1. - ! - ! - ! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation. - ! - implicit none - procedure(d_func_3d) :: b1,b2,b3,c,a1,a2,a3,g - integer(psb_ipk_) :: idim - type(psb_dspmat_type) :: a - type(psb_d_vect_type) :: xv,bv - type(psb_desc_type) :: desc_a - integer(psb_ipk_) :: ictxt, info - character :: afmt*5 - - ! Local variables. - - integer(psb_ipk_), parameter :: nb=20 - type(psb_d_csc_sparse_mat) :: acsc - type(psb_d_coo_sparse_mat) :: acoo - type(psb_d_csr_sparse_mat) :: acsr - real(psb_dpk_) :: zt(nb),x,y,z - integer(psb_ipk_) :: m,n,nnz,glob_row,nlr,i,ii,ib,k - integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner - integer(psb_ipk_) :: np, iam, nr, nt - integer(psb_ipk_) :: icoeff - integer(psb_ipk_), allocatable :: irow(:),icol(:),myidx(:) - real(psb_dpk_), allocatable :: val(:) - ! deltah dimension of each grid cell - ! deltat discretization time - real(psb_dpk_) :: deltah, sqdeltah, deltah2 - real(psb_dpk_), parameter :: rhs=0.d0,one=1.d0,zero=0.d0 - real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb - integer(psb_ipk_) :: err_act - - character(len=20) :: name, ch_err,tmpfmt - - info = psb_success_ - name = 'create_matrix' - call psb_erractionsave(err_act) - - call psb_info(ictxt, iam, np) - - deltah = 1.d0/(idim+2) - sqdeltah = deltah*deltah - deltah2 = 2.d0* deltah - - ! initialize array descriptor and sparse matrix storage. provide an - ! estimate of the number of non zeroes - - m = idim*idim*idim - n = m - nnz = ((n*9)/(np)) - if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n - - ! - ! Using a simple BLOCK distribution. - ! - nt = (m+np-1)/np - nr = max(0,min(nt,m-(iam*nt))) - - nt = nr - call psb_sum(ictxt,nt) - if (nt /= m) write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m - call psb_barrier(ictxt) - t0 = psb_wtime() - call psb_cdall(ictxt,desc_a,info,nl=nr) - if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz) - ! define rhs from boundary conditions; also build initial guess - if (info == psb_success_) call psb_geall(xv,desc_a,info) - if (info == psb_success_) call psb_geall(bv,desc_a,info) - nlr = desc_a%get_local_rows() - call psb_barrier(ictxt) - talc = psb_wtime()-t0 - - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='allocation rout.' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - ! we build an auxiliary matrix consisting of one row at a - ! time; just a small matrix. might be extended to generate - ! a bunch of rows per call. - ! - allocate(val(20*nb),irow(20*nb),& - &icol(20*nb),myidx(nlr),stat=info) - if (info /= psb_success_ ) then - info=psb_err_alloc_dealloc_ - call psb_errpush(info,name) - goto 9999 - endif - - do i=1,nlr - myidx(i) = i - end do - - - call psb_loc_to_glob(myidx,desc_a,info) - - ! loop over rows belonging to current process in a block - ! distribution. - - call psb_barrier(ictxt) - t1 = psb_wtime() - do ii=1, nlr,nb - ib = min(nb,nlr-ii+1) - icoeff = 1 - do k=1,ib - i=ii+k-1 - ! local matrix pointer - glob_row=myidx(i) - ! compute gridpoint coordinates - if (mod(glob_row,(idim*idim)) == 0) then - ix = glob_row/(idim*idim) - else - ix = glob_row/(idim*idim)+1 - endif - if (mod((glob_row-(ix-1)*idim*idim),idim) == 0) then - iy = (glob_row-(ix-1)*idim*idim)/idim - else - iy = (glob_row-(ix-1)*idim*idim)/idim+1 - endif - iz = glob_row-(ix-1)*idim*idim-(iy-1)*idim - ! x, y, x coordinates - x = ix*deltah - y = iy*deltah - z = iz*deltah - zt(k) = 0.d0 - ! internal point: build discretization - ! - ! term depending on (x-1,y,z) - ! - val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2 - if (ix == 1) then - zt(k) = g(dzero,y,z)*(-val(icoeff)) - else - icol(icoeff) = (ix-2)*idim*idim+(iy-1)*idim+(iz) - irow(icoeff) = glob_row - icoeff = icoeff+1 - endif - ! term depending on (x,y-1,z) - val(icoeff) = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2 - if (iy == 1) then - zt(k) = g(x,dzero,z)*(-val(icoeff)) - else - icol(icoeff) = (ix-1)*idim*idim+(iy-2)*idim+(iz) - irow(icoeff) = glob_row - icoeff = icoeff+1 - endif - ! term depending on (x,y,z-1) - val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2 - if (iz == 1) then - zt(k) = g(x,y,dzero)*(-val(icoeff)) - else - icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz-1) - irow(icoeff) = glob_row - icoeff = icoeff+1 - endif - - ! term depending on (x,y,z) - val(icoeff)=(2*a1(x,y,z) + 2*a2(x,y,z) + 2*a3(x,y,z))/sqdeltah& - & +c(x,y,z) - icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz) - irow(icoeff) = glob_row - icoeff = icoeff+1 - ! term depending on (x,y,z+1) - val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2 - if (iz == idim) then - zt(k) = g(x,y,done)*(-val(icoeff)) - else - icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz+1) - irow(icoeff) = glob_row - icoeff = icoeff+1 - endif - ! term depending on (x,y+1,z) - val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2 - if (iy == idim) then - zt(k) = g(x,done,z)*(-val(icoeff)) - else - icol(icoeff) = (ix-1)*idim*idim+(iy)*idim+(iz) - irow(icoeff) = glob_row - icoeff = icoeff+1 - endif - ! term depending on (x+1,y,z) - val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2 - if (ix==idim) then - zt(k) = g(done,y,z)*(-val(icoeff)) - else - icol(icoeff) = (ix)*idim*idim+(iy-1)*idim+(iz) - irow(icoeff) = glob_row - icoeff = icoeff+1 - endif - - end do - call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info) - if(info /= psb_success_) exit - call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info) - if(info /= psb_success_) exit - zt(:)=0.d0 - call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info) - if(info /= psb_success_) exit - end do - - tgen = psb_wtime()-t1 - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='insert rout.' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - deallocate(val,irow,icol) - - call psb_barrier(ictxt) - t1 = psb_wtime() - call psb_cdasb(desc_a,info) - tcdasb = psb_wtime()-t1 - call psb_barrier(ictxt) - t1 = psb_wtime() - if (info == psb_success_) & - & call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) - call psb_barrier(ictxt) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='asb rout.' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - if (info == psb_success_) call psb_geasb(xv,desc_a,info) - if (info == psb_success_) call psb_geasb(bv,desc_a,info) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='asb rout.' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - tasb = psb_wtime()-t1 - call psb_barrier(ictxt) - ttot = psb_wtime() - t0 - - call psb_amx(ictxt,talc) - call psb_amx(ictxt,tgen) - call psb_amx(ictxt,tasb) - call psb_amx(ictxt,ttot) - if(iam == psb_root_) then - tmpfmt = a%get_fmt() - write(psb_out_unit,'("The matrix has been generated and assembled in ",a3," format.")')& - & tmpfmt - write(psb_out_unit,'("-allocation time : ",es12.5)') talc - write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen - write(psb_out_unit,'("-desc asbly time : ",es12.5)') tcdasb - write(psb_out_unit,'("- mat asbly time : ",es12.5)') tasb - write(psb_out_unit,'("-total time : ",es12.5)') ttot - - end if - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error(ictxt) - return - end if - return -end subroutine gen_prob3d diff --git a/util/psb_d_genmat_mod.f90 b/util/psb_d_genmat_mod.f90 deleted file mode 100644 index e66bde9f..00000000 --- a/util/psb_d_genmat_mod.f90 +++ /dev/null @@ -1,41 +0,0 @@ -module psb_d_genmat_mod - use psb_base_mod - interface - function d_func_3d(x,y,z) result(val) - import :: psb_dpk_ - real(psb_dpk_), intent(in) :: x,y,z - real(psb_dpk_) :: val - end function d_func_3d - end interface - - interface - subroutine gen_prob3d(ictxt,idim,a,bv,xv,desc_a,afmt,a1,a2,a3,b1,b2,b3,c,g,info) - ! - ! Discretizes the partial differential equation - ! - ! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u) - ! - ------ - ------ - ------ + ----- + ------ + ------ + c u = 0 - ! dxdx dydy dzdz dx dy dz - ! - ! with Dirichlet boundary conditions - ! u = g - ! - ! on the unit cube 0<=x,y,z<=1. - ! - ! - ! Note that if a1=a2=a3=c=0., the PDE is the Laplace equation. - ! - import :: psb_ipk_, psb_desc_type, psb_dspmat_type, psb_d_vect_type, d_func_3d - implicit none - procedure(d_func_3d) :: a1,a2,a3,c,b1,b2,b3,g - integer(psb_ipk_) :: idim - type(psb_dspmat_type) :: a - type(psb_d_vect_type) :: xv,bv - type(psb_desc_type) :: desc_a - integer(psb_ipk_) :: ictxt, info - character :: afmt*5 - end subroutine gen_prob3d - end interface - - -end module psb_d_genmat_mod diff --git a/util/psb_d_genpde_impl.f90 b/util/psb_d_genpde_impl.f90 new file mode 100644 index 00000000..bac64daf --- /dev/null +++ b/util/psb_d_genpde_impl.f90 @@ -0,0 +1,563 @@ +! +! subroutine to allocate and fill in the coefficient matrix and +! the rhs. +! +subroutine psb_d_gen_prob3d(ictxt,idim,a,bv,xv,desc_a,afmt,a1,a2,a3,b1,b2,b3,c,g,info,f) + use psb_base_mod + use psb_d_genpde_mod, psb_protect_name => psb_d_gen_prob3d + ! + ! Discretizes the partial differential equation + ! + ! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u) + ! - ------ - ------ - ------ + ----- + ------ + ------ + c u = f + ! dxdx dydy dzdz dx dy dz + ! + ! with Dirichlet boundary conditions + ! u = g + ! + ! on the unit cube 0<=x,y,z<=1. + ! + ! + ! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation. + ! + implicit none + procedure(d_func_3d) :: b1,b2,b3,c,a1,a2,a3,g + integer(psb_ipk_) :: idim + type(psb_dspmat_type) :: a + type(psb_d_vect_type) :: xv,bv + type(psb_desc_type) :: desc_a + integer(psb_ipk_) :: ictxt, info + character :: afmt*5 + procedure(d_func_3d), optional :: f + + ! Local variables. + + integer(psb_ipk_), parameter :: nb=20 + type(psb_d_csc_sparse_mat) :: acsc + type(psb_d_coo_sparse_mat) :: acoo + type(psb_d_csr_sparse_mat) :: acsr + real(psb_dpk_) :: zt(nb),x,y,z + integer(psb_ipk_) :: m,n,nnz,glob_row,nlr,i,ii,ib,k + integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner + integer(psb_ipk_) :: np, iam, nr, nt + integer(psb_ipk_) :: icoeff + integer(psb_ipk_), allocatable :: irow(:),icol(:),myidx(:) + real(psb_dpk_), allocatable :: val(:) + ! deltah dimension of each grid cell + ! deltat discretization time + real(psb_dpk_) :: deltah, sqdeltah, deltah2 + real(psb_dpk_), parameter :: rhs=0.d0,one=1.d0,zero=0.d0 + real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb + integer(psb_ipk_) :: err_act + procedure(d_func_3d), pointer :: f_ + character(len=20) :: name, ch_err,tmpfmt + + info = psb_success_ + name = 'create_matrix' + call psb_erractionsave(err_act) + + call psb_info(ictxt, iam, np) + + + if (present(f)) then + f_ => f + else + f_ => d_null_func_3d + end if + + deltah = 1.d0/(idim+2) + sqdeltah = deltah*deltah + deltah2 = 2.d0* deltah + + ! initialize array descriptor and sparse matrix storage. provide an + ! estimate of the number of non zeroes + + m = idim*idim*idim + n = m + nnz = ((n*9)/(np)) + if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n + + ! + ! Using a simple BLOCK distribution. + ! + nt = (m+np-1)/np + nr = max(0,min(nt,m-(iam*nt))) + + nt = nr + call psb_sum(ictxt,nt) + if (nt /= m) write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m + call psb_barrier(ictxt) + t0 = psb_wtime() + call psb_cdall(ictxt,desc_a,info,nl=nr) + if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz) + ! define rhs from boundary conditions; also build initial guess + if (info == psb_success_) call psb_geall(xv,desc_a,info) + if (info == psb_success_) call psb_geall(bv,desc_a,info) + nlr = desc_a%get_local_rows() + call psb_barrier(ictxt) + talc = psb_wtime()-t0 + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='allocation rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + ! we build an auxiliary matrix consisting of one row at a + ! time; just a small matrix. might be extended to generate + ! a bunch of rows per call. + ! + allocate(val(20*nb),irow(20*nb),& + &icol(20*nb),myidx(nlr),stat=info) + if (info /= psb_success_ ) then + info=psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + endif + + do i=1,nlr + myidx(i) = i + end do + + + call psb_loc_to_glob(myidx,desc_a,info) + + ! loop over rows belonging to current process in a block + ! distribution. + + call psb_barrier(ictxt) + t1 = psb_wtime() + do ii=1, nlr,nb + ib = min(nb,nlr-ii+1) + icoeff = 1 + do k=1,ib + i=ii+k-1 + ! local matrix pointer + glob_row=myidx(i) + ! compute gridpoint coordinates + if (mod(glob_row,(idim*idim)) == 0) then + ix = glob_row/(idim*idim) + else + ix = glob_row/(idim*idim)+1 + endif + if (mod((glob_row-(ix-1)*idim*idim),idim) == 0) then + iy = (glob_row-(ix-1)*idim*idim)/idim + else + iy = (glob_row-(ix-1)*idim*idim)/idim+1 + endif + iz = glob_row-(ix-1)*idim*idim-(iy-1)*idim + ! x, y, x coordinates + x = ix*deltah + y = iy*deltah + z = iz*deltah + zt(k) = f_(x,y,z) + ! internal point: build discretization + ! + ! term depending on (x-1,y,z) + ! + val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2 + if (ix == 1) then + zt(k) = g(dzero,y,z)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-2)*idim*idim+(iy-1)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y-1,z) + val(icoeff) = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2 + if (iy == 1) then + zt(k) = g(x,dzero,z)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim*idim+(iy-2)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y,z-1) + val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2 + if (iz == 1) then + zt(k) = g(x,y,dzero)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz-1) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + ! term depending on (x,y,z) + val(icoeff)=2.d0*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah & + & + c(x,y,z) + icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + ! term depending on (x,y,z+1) + val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2 + if (iz == idim) then + zt(k) = g(x,y,done)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz+1) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y+1,z) + val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2 + if (iy == idim) then + zt(k) = g(x,done,z)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim*idim+(iy)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x+1,y,z) + val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2 + if (ix==idim) then + zt(k) = g(done,y,z)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix)*idim*idim+(iy-1)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + end do + call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info) + if(info /= psb_success_) exit + call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info) + if(info /= psb_success_) exit + zt(:)=0.d0 + call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info) + if(info /= psb_success_) exit + end do + + tgen = psb_wtime()-t1 + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='insert rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + deallocate(val,irow,icol) + + call psb_barrier(ictxt) + t1 = psb_wtime() + call psb_cdasb(desc_a,info) + tcdasb = psb_wtime()-t1 + call psb_barrier(ictxt) + t1 = psb_wtime() + if (info == psb_success_) & + & call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) + call psb_barrier(ictxt) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='asb rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + if (info == psb_success_) call psb_geasb(xv,desc_a,info) + if (info == psb_success_) call psb_geasb(bv,desc_a,info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='asb rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + tasb = psb_wtime()-t1 + call psb_barrier(ictxt) + ttot = psb_wtime() - t0 + + call psb_amx(ictxt,talc) + call psb_amx(ictxt,tgen) + call psb_amx(ictxt,tasb) + call psb_amx(ictxt,ttot) + if(iam == psb_root_) then + tmpfmt = a%get_fmt() + write(psb_out_unit,'("The matrix has been generated and assembled in ",a3," format.")')& + & tmpfmt + write(psb_out_unit,'("-allocation time : ",es12.5)') talc + write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen + write(psb_out_unit,'("-desc asbly time : ",es12.5)') tcdasb + write(psb_out_unit,'("- mat asbly time : ",es12.5)') tasb + write(psb_out_unit,'("-total time : ",es12.5)') ttot + + end if + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + return +end subroutine psb_d_gen_prob3d + + + +! +! subroutine to allocate and fill in the coefficient matrix and +! the rhs. +! +subroutine psb_d_gen_prob2d(ictxt,idim,a,bv,xv,desc_a,afmt,a1,a2,b1,b2,c,g,info,f) + use psb_base_mod + use psb_d_genpde_mod, psb_protect_name => psb_d_gen_prob2d + ! + ! Discretizes the partial differential equation + ! + ! a1 dd(u) a2 dd(u) b1 d(u) b2 d(u) + ! - ------ - ------ + ----- + ------ + c u = f + ! dxdx dydy dx dy + ! + ! with Dirichlet boundary conditions + ! u = g + ! + ! on the unit square 0<=x,y<=1. + ! + ! + ! Note that if b1=b2=c=0., the PDE is the Laplace equation. + ! + implicit none + procedure(d_func_2d) :: b1,b2,c,a1,a2,g + integer(psb_ipk_) :: idim + type(psb_dspmat_type) :: a + type(psb_d_vect_type) :: xv,bv + type(psb_desc_type) :: desc_a + integer(psb_ipk_) :: ictxt, info + character :: afmt*5 + procedure(d_func_2d), optional :: f + + ! Local variables. + + integer(psb_ipk_), parameter :: nb=20 + type(psb_d_csc_sparse_mat) :: acsc + type(psb_d_coo_sparse_mat) :: acoo + type(psb_d_csr_sparse_mat) :: acsr + real(psb_dpk_) :: zt(nb),x,y,z + integer(psb_ipk_) :: m,n,nnz,glob_row,nlr,i,ii,ib,k + integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner + integer(psb_ipk_) :: np, iam, nr, nt + integer(psb_ipk_) :: icoeff + integer(psb_ipk_), allocatable :: irow(:),icol(:),myidx(:) + real(psb_dpk_), allocatable :: val(:) + ! deltah dimension of each grid cell + ! deltat discretization time + real(psb_dpk_) :: deltah, sqdeltah, deltah2 + real(psb_dpk_), parameter :: rhs=0.d0,one=1.d0,zero=0.d0 + real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb + integer(psb_ipk_) :: err_act + procedure(d_func_2d), pointer :: f_ + character(len=20) :: name, ch_err,tmpfmt + + info = psb_success_ + name = 'create_matrix' + call psb_erractionsave(err_act) + + call psb_info(ictxt, iam, np) + + + if (present(f)) then + f_ => f + else + f_ => d_null_func_2d + end if + + deltah = 1.d0/(idim+2) + sqdeltah = deltah*deltah + deltah2 = 2.d0* deltah + + ! initialize array descriptor and sparse matrix storage. provide an + ! estimate of the number of non zeroes + + m = idim*idim + n = m + nnz = ((n*7)/(np)) + if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n + + ! + ! Using a simple BLOCK distribution. + ! + nt = (m+np-1)/np + nr = max(0,min(nt,m-(iam*nt))) + + nt = nr + call psb_sum(ictxt,nt) + if (nt /= m) write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m + call psb_barrier(ictxt) + t0 = psb_wtime() + call psb_cdall(ictxt,desc_a,info,nl=nr) + if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz) + ! define rhs from boundary conditions; also build initial guess + if (info == psb_success_) call psb_geall(xv,desc_a,info) + if (info == psb_success_) call psb_geall(bv,desc_a,info) + nlr = desc_a%get_local_rows() + call psb_barrier(ictxt) + talc = psb_wtime()-t0 + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='allocation rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + ! we build an auxiliary matrix consisting of one row at a + ! time; just a small matrix. might be extended to generate + ! a bunch of rows per call. + ! + allocate(val(20*nb),irow(20*nb),& + &icol(20*nb),myidx(nlr),stat=info) + if (info /= psb_success_ ) then + info=psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + endif + + do i=1,nlr + myidx(i) = i + end do + + + call psb_loc_to_glob(myidx,desc_a,info) + + ! loop over rows belonging to current process in a block + ! distribution. + + call psb_barrier(ictxt) + t1 = psb_wtime() + do ii=1, nlr,nb + ib = min(nb,nlr-ii+1) + icoeff = 1 + do k=1,ib + i=ii+k-1 + ! local matrix pointer + glob_row=myidx(i) + ! compute gridpoint coordinates + if (mod(glob_row,(idim)) == 0) then + ix = glob_row/(idim) + else + ix = glob_row/(idim)+1 + endif + iy = (glob_row-(ix-1)*idim) + ! x, y + x = ix*deltah + y = iy*deltah + + zt(k) = f_(x,y) + ! internal point: build discretization + ! + ! term depending on (x-1,y) + ! + val(icoeff) = -a1(x,y)/sqdeltah-b1(x,y)/deltah2 + if (ix == 1) then + zt(k) = g(dzero,y)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-2)*idim+iy + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y-1) + val(icoeff) = -a2(x,y)/sqdeltah-b2(x,y)/deltah2 + if (iy == 1) then + zt(k) = g(x,dzero)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim+(iy-1) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + ! term depending on (x,y) + val(icoeff)=2.d0*(a1(x,y) + a2(x,y))/sqdeltah + c(x,y) + icol(icoeff) = (ix-1)*idim+iy + irow(icoeff) = glob_row + icoeff = icoeff+1 + ! term depending on (x,y+1) + val(icoeff)=-a2(x,y)/sqdeltah+b2(x,y)/deltah2 + if (iy == idim) then + zt(k) = g(x,done)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim+(iy+1) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x+1,y) + val(icoeff)=-a1(x,y)/sqdeltah+b1(x,y)/deltah2 + if (ix==idim) then + zt(k) = g(done,y)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix)*idim+(iy) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + end do + call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info) + if(info /= psb_success_) exit + call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info) + if(info /= psb_success_) exit + zt(:)=0.d0 + call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info) + if(info /= psb_success_) exit + end do + + tgen = psb_wtime()-t1 + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='insert rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + deallocate(val,irow,icol) + + call psb_barrier(ictxt) + t1 = psb_wtime() + call psb_cdasb(desc_a,info) + tcdasb = psb_wtime()-t1 + call psb_barrier(ictxt) + t1 = psb_wtime() + if (info == psb_success_) & + & call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) + call psb_barrier(ictxt) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='asb rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + if (info == psb_success_) call psb_geasb(xv,desc_a,info) + if (info == psb_success_) call psb_geasb(bv,desc_a,info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='asb rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + tasb = psb_wtime()-t1 + call psb_barrier(ictxt) + ttot = psb_wtime() - t0 + + call psb_amx(ictxt,talc) + call psb_amx(ictxt,tgen) + call psb_amx(ictxt,tasb) + call psb_amx(ictxt,ttot) + if(iam == psb_root_) then + tmpfmt = a%get_fmt() + write(psb_out_unit,'("The matrix has been generated and assembled in ",a3," format.")')& + & tmpfmt + write(psb_out_unit,'("-allocation time : ",es12.5)') talc + write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen + write(psb_out_unit,'("-desc asbly time : ",es12.5)') tcdasb + write(psb_out_unit,'("- mat asbly time : ",es12.5)') tasb + write(psb_out_unit,'("-total time : ",es12.5)') ttot + + end if + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + return +end subroutine psb_d_gen_prob2d diff --git a/util/psb_d_genpde_mod.f90 b/util/psb_d_genpde_mod.f90 new file mode 100644 index 00000000..9935c46a --- /dev/null +++ b/util/psb_d_genpde_mod.f90 @@ -0,0 +1,105 @@ +module psb_d_genpde_mod + + use psb_base_mod, only : psb_dpk_, psb_ipk_, psb_desc_type,& + & psb_dspmat_type, psb_d_vect_type, dzero + + 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 psb_gen_prob3d + subroutine psb_d_gen_prob3d(ictxt,idim,a,bv,xv,desc_a,afmt, & + & a1,a2,a3,b1,b2,b3,c,g,info,f) + ! + ! 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. + ! + import :: psb_ipk_, psb_desc_type, psb_dspmat_type, psb_d_vect_type, d_func_3d + implicit none + procedure(d_func_3d) :: a1,a2,a3,c,b1,b2,b3,g + integer(psb_ipk_) :: idim + type(psb_dspmat_type) :: a + type(psb_d_vect_type) :: xv,bv + type(psb_desc_type) :: desc_a + integer(psb_ipk_) :: ictxt, info + character :: afmt*5 + procedure(d_func_3d), optional :: f + end subroutine psb_d_gen_prob3d + end interface + + + 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 psb_gen_prob2d + subroutine psb_d_gen_prob2d(ictxt,idim,a,bv,xv,desc_a,afmt,a1,a2,b1,b2,c,g,info,f) + ! + ! 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. + ! + import :: psb_ipk_, psb_desc_type, psb_dspmat_type, psb_d_vect_type, d_func_2d + implicit none + procedure(d_func_2d) :: a1,a2,c,b1,b2,g + integer(psb_ipk_) :: idim + type(psb_dspmat_type) :: a + type(psb_d_vect_type) :: xv,bv + type(psb_desc_type) :: desc_a + integer(psb_ipk_) :: ictxt, info + character :: afmt*5 + procedure(d_func_2d), optional :: f + end subroutine psb_d_gen_prob2d + end interface + +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 + + 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 + + +end module psb_d_genpde_mod diff --git a/util/psb_s_genpde_impl.f90 b/util/psb_s_genpde_impl.f90 new file mode 100644 index 00000000..fe362135 --- /dev/null +++ b/util/psb_s_genpde_impl.f90 @@ -0,0 +1,563 @@ +! +! subroutine to allocate and fill in the coefficient matrix and +! the rhs. +! +subroutine psb_s_gen_prob3d(ictxt,idim,a,bv,xv,desc_a,afmt,a1,a2,a3,b1,b2,b3,c,g,info,f) + use psb_base_mod + use psb_s_genpde_mod, psb_protect_name => psb_s_gen_prob3d + ! + ! Discretizes the partial differential equation + ! + ! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u) + ! - ------ - ------ - ------ + ----- + ------ + ------ + c u = f + ! dxdx dydy dzdz dx dy dz + ! + ! with Dirichlet boundary conditions + ! u = g + ! + ! on the unit cube 0<=x,y,z<=1. + ! + ! + ! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation. + ! + implicit none + procedure(d_func_3d) :: b1,b2,b3,c,a1,a2,a3,g + integer(psb_ipk_) :: idim + type(psb_sspmat_type) :: a + type(psb_s_vect_type) :: xv,bv + type(psb_desc_type) :: desc_a + integer(psb_ipk_) :: ictxt, info + character :: afmt*5 + procedure(d_func_3d), optional :: f + + ! Local variables. + + integer(psb_ipk_), parameter :: nb=20 + type(psb_s_csc_sparse_mat) :: acsc + type(psb_s_coo_sparse_mat) :: acoo + type(psb_s_csr_sparse_mat) :: acsr + real(psb_spk_) :: zt(nb),x,y,z + integer(psb_ipk_) :: m,n,nnz,glob_row,nlr,i,ii,ib,k + integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner + integer(psb_ipk_) :: np, iam, nr, nt + integer(psb_ipk_) :: icoeff + integer(psb_ipk_), allocatable :: irow(:),icol(:),myidx(:) + real(psb_spk_), allocatable :: val(:) + ! deltah dimension of each grid cell + ! deltat discretization time + real(psb_spk_) :: deltah, sqdeltah, deltah2 + real(psb_spk_), parameter :: rhs=0.e0,one=1.e0,zero=0.e0 + real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb + integer(psb_ipk_) :: err_act + procedure(d_func_3d), pointer :: f_ + character(len=20) :: name, ch_err,tmpfmt + + info = psb_success_ + name = 'create_matrix' + call psb_erractionsave(err_act) + + call psb_info(ictxt, iam, np) + + + if (present(f)) then + f_ => f + else + f_ => d_null_func_3d + end if + + deltah = 1.e0/(idim+2) + sqdeltah = deltah*deltah + deltah2 = 2.e0* deltah + + ! initialize array descriptor and sparse matrix storage. provide an + ! estimate of the number of non zeroes + + m = idim*idim*idim + n = m + nnz = ((n*9)/(np)) + if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n + + ! + ! Using a simple BLOCK distribution. + ! + nt = (m+np-1)/np + nr = max(0,min(nt,m-(iam*nt))) + + nt = nr + call psb_sum(ictxt,nt) + if (nt /= m) write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m + call psb_barrier(ictxt) + t0 = psb_wtime() + call psb_cdall(ictxt,desc_a,info,nl=nr) + if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz) + ! define rhs from boundary conditions; also build initial guess + if (info == psb_success_) call psb_geall(xv,desc_a,info) + if (info == psb_success_) call psb_geall(bv,desc_a,info) + nlr = desc_a%get_local_rows() + call psb_barrier(ictxt) + talc = psb_wtime()-t0 + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='allocation rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + ! we build an auxiliary matrix consisting of one row at a + ! time; just a small matrix. might be extended to generate + ! a bunch of rows per call. + ! + allocate(val(20*nb),irow(20*nb),& + &icol(20*nb),myidx(nlr),stat=info) + if (info /= psb_success_ ) then + info=psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + endif + + do i=1,nlr + myidx(i) = i + end do + + + call psb_loc_to_glob(myidx,desc_a,info) + + ! loop over rows belonging to current process in a block + ! distribution. + + call psb_barrier(ictxt) + t1 = psb_wtime() + do ii=1, nlr,nb + ib = min(nb,nlr-ii+1) + icoeff = 1 + do k=1,ib + i=ii+k-1 + ! local matrix pointer + glob_row=myidx(i) + ! compute gridpoint coordinates + if (mod(glob_row,(idim*idim)) == 0) then + ix = glob_row/(idim*idim) + else + ix = glob_row/(idim*idim)+1 + endif + if (mod((glob_row-(ix-1)*idim*idim),idim) == 0) then + iy = (glob_row-(ix-1)*idim*idim)/idim + else + iy = (glob_row-(ix-1)*idim*idim)/idim+1 + endif + iz = glob_row-(ix-1)*idim*idim-(iy-1)*idim + ! x, y, x coordinates + x = ix*deltah + y = iy*deltah + z = iz*deltah + zt(k) = f_(x,y,z) + ! internal point: build discretization + ! + ! term depending on (x-1,y,z) + ! + val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2 + if (ix == 1) then + zt(k) = g(szero,y,z)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-2)*idim*idim+(iy-1)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y-1,z) + val(icoeff) = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2 + if (iy == 1) then + zt(k) = g(x,szero,z)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim*idim+(iy-2)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y,z-1) + val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2 + if (iz == 1) then + zt(k) = g(x,y,szero)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz-1) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + ! term depending on (x,y,z) + val(icoeff)=2.e0*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah & + & + c(x,y,z) + icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + ! term depending on (x,y,z+1) + val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2 + if (iz == idim) then + zt(k) = g(x,y,sone)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz+1) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y+1,z) + val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2 + if (iy == idim) then + zt(k) = g(x,sone,z)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim*idim+(iy)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x+1,y,z) + val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2 + if (ix==idim) then + zt(k) = g(sone,y,z)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix)*idim*idim+(iy-1)*idim+(iz) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + end do + call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info) + if(info /= psb_success_) exit + call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info) + if(info /= psb_success_) exit + zt(:)=0.e0 + call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info) + if(info /= psb_success_) exit + end do + + tgen = psb_wtime()-t1 + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='insert rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + deallocate(val,irow,icol) + + call psb_barrier(ictxt) + t1 = psb_wtime() + call psb_cdasb(desc_a,info) + tcdasb = psb_wtime()-t1 + call psb_barrier(ictxt) + t1 = psb_wtime() + if (info == psb_success_) & + & call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) + call psb_barrier(ictxt) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='asb rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + if (info == psb_success_) call psb_geasb(xv,desc_a,info) + if (info == psb_success_) call psb_geasb(bv,desc_a,info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='asb rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + tasb = psb_wtime()-t1 + call psb_barrier(ictxt) + ttot = psb_wtime() - t0 + + call psb_amx(ictxt,talc) + call psb_amx(ictxt,tgen) + call psb_amx(ictxt,tasb) + call psb_amx(ictxt,ttot) + if(iam == psb_root_) then + tmpfmt = a%get_fmt() + write(psb_out_unit,'("The matrix has been generated and assembled in ",a3," format.")')& + & tmpfmt + write(psb_out_unit,'("-allocation time : ",es12.5)') talc + write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen + write(psb_out_unit,'("-desc asbly time : ",es12.5)') tcdasb + write(psb_out_unit,'("- mat asbly time : ",es12.5)') tasb + write(psb_out_unit,'("-total time : ",es12.5)') ttot + + end if + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + return +end subroutine psb_s_gen_prob3d + + + +! +! subroutine to allocate and fill in the coefficient matrix and +! the rhs. +! +subroutine psb_s_gen_prob2d(ictxt,idim,a,bv,xv,desc_a,afmt,a1,a2,b1,b2,c,g,info,f) + use psb_base_mod + use psb_s_genpde_mod, psb_protect_name => psb_s_gen_prob2d + ! + ! Discretizes the partial differential equation + ! + ! a1 dd(u) a2 dd(u) b1 d(u) b2 d(u) + ! - ------ - ------ + ----- + ------ + c u = f + ! dxdx dydy dx dy + ! + ! with Dirichlet boundary conditions + ! u = g + ! + ! on the unit square 0<=x,y<=1. + ! + ! + ! Note that if b1=b2=c=0., the PDE is the Laplace equation. + ! + implicit none + procedure(d_func_2d) :: b1,b2,c,a1,a2,g + integer(psb_ipk_) :: idim + type(psb_sspmat_type) :: a + type(psb_s_vect_type) :: xv,bv + type(psb_desc_type) :: desc_a + integer(psb_ipk_) :: ictxt, info + character :: afmt*5 + procedure(d_func_2d), optional :: f + + ! Local variables. + + integer(psb_ipk_), parameter :: nb=20 + type(psb_s_csc_sparse_mat) :: acsc + type(psb_s_coo_sparse_mat) :: acoo + type(psb_s_csr_sparse_mat) :: acsr + real(psb_spk_) :: zt(nb),x,y,z + integer(psb_ipk_) :: m,n,nnz,glob_row,nlr,i,ii,ib,k + integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner + integer(psb_ipk_) :: np, iam, nr, nt + integer(psb_ipk_) :: icoeff + integer(psb_ipk_), allocatable :: irow(:),icol(:),myidx(:) + real(psb_spk_), allocatable :: val(:) + ! deltah dimension of each grid cell + ! deltat discretization time + real(psb_spk_) :: deltah, sqdeltah, deltah2 + real(psb_spk_), parameter :: rhs=0.e0,one=1.e0,zero=0.e0 + real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb + integer(psb_ipk_) :: err_act + procedure(d_func_2d), pointer :: f_ + character(len=20) :: name, ch_err,tmpfmt + + info = psb_success_ + name = 'create_matrix' + call psb_erractionsave(err_act) + + call psb_info(ictxt, iam, np) + + + if (present(f)) then + f_ => f + else + f_ => d_null_func_2d + end if + + deltah = 1.e0/(idim+2) + sqdeltah = deltah*deltah + deltah2 = 2.e0* deltah + + ! initialize array descriptor and sparse matrix storage. provide an + ! estimate of the number of non zeroes + + m = idim*idim + n = m + nnz = ((n*7)/(np)) + if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n + + ! + ! Using a simple BLOCK distribution. + ! + nt = (m+np-1)/np + nr = max(0,min(nt,m-(iam*nt))) + + nt = nr + call psb_sum(ictxt,nt) + if (nt /= m) write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m + call psb_barrier(ictxt) + t0 = psb_wtime() + call psb_cdall(ictxt,desc_a,info,nl=nr) + if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz) + ! define rhs from boundary conditions; also build initial guess + if (info == psb_success_) call psb_geall(xv,desc_a,info) + if (info == psb_success_) call psb_geall(bv,desc_a,info) + nlr = desc_a%get_local_rows() + call psb_barrier(ictxt) + talc = psb_wtime()-t0 + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='allocation rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + ! we build an auxiliary matrix consisting of one row at a + ! time; just a small matrix. might be extended to generate + ! a bunch of rows per call. + ! + allocate(val(20*nb),irow(20*nb),& + &icol(20*nb),myidx(nlr),stat=info) + if (info /= psb_success_ ) then + info=psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + endif + + do i=1,nlr + myidx(i) = i + end do + + + call psb_loc_to_glob(myidx,desc_a,info) + + ! loop over rows belonging to current process in a block + ! distribution. + + call psb_barrier(ictxt) + t1 = psb_wtime() + do ii=1, nlr,nb + ib = min(nb,nlr-ii+1) + icoeff = 1 + do k=1,ib + i=ii+k-1 + ! local matrix pointer + glob_row=myidx(i) + ! compute gridpoint coordinates + if (mod(glob_row,(idim)) == 0) then + ix = glob_row/(idim) + else + ix = glob_row/(idim)+1 + endif + iy = (glob_row-(ix-1)*idim) + ! x, y + x = ix*deltah + y = iy*deltah + + zt(k) = f_(x,y) + ! internal point: build discretization + ! + ! term depending on (x-1,y) + ! + val(icoeff) = -a1(x,y)/sqdeltah-b1(x,y)/deltah2 + if (ix == 1) then + zt(k) = g(szero,y)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-2)*idim+iy + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y-1) + val(icoeff) = -a2(x,y)/sqdeltah-b2(x,y)/deltah2 + if (iy == 1) then + zt(k) = g(x,szero)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim+(iy-1) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + ! term depending on (x,y) + val(icoeff)=2.e0*(a1(x,y) + a2(x,y))/sqdeltah + c(x,y) + icol(icoeff) = (ix-1)*idim+iy + irow(icoeff) = glob_row + icoeff = icoeff+1 + ! term depending on (x,y+1) + val(icoeff)=-a2(x,y)/sqdeltah+b2(x,y)/deltah2 + if (iy == idim) then + zt(k) = g(x,sone)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix-1)*idim+(iy+1) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x+1,y) + val(icoeff)=-a1(x,y)/sqdeltah+b1(x,y)/deltah2 + if (ix==idim) then + zt(k) = g(sone,y)*(-val(icoeff)) + zt(k) + else + icol(icoeff) = (ix)*idim+(iy) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + end do + call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info) + if(info /= psb_success_) exit + call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info) + if(info /= psb_success_) exit + zt(:)=0.e0 + call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info) + if(info /= psb_success_) exit + end do + + tgen = psb_wtime()-t1 + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='insert rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + deallocate(val,irow,icol) + + call psb_barrier(ictxt) + t1 = psb_wtime() + call psb_cdasb(desc_a,info) + tcdasb = psb_wtime()-t1 + call psb_barrier(ictxt) + t1 = psb_wtime() + if (info == psb_success_) & + & call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) + call psb_barrier(ictxt) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='asb rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + if (info == psb_success_) call psb_geasb(xv,desc_a,info) + if (info == psb_success_) call psb_geasb(bv,desc_a,info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='asb rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + tasb = psb_wtime()-t1 + call psb_barrier(ictxt) + ttot = psb_wtime() - t0 + + call psb_amx(ictxt,talc) + call psb_amx(ictxt,tgen) + call psb_amx(ictxt,tasb) + call psb_amx(ictxt,ttot) + if(iam == psb_root_) then + tmpfmt = a%get_fmt() + write(psb_out_unit,'("The matrix has been generated and assembled in ",a3," format.")')& + & tmpfmt + write(psb_out_unit,'("-allocation time : ",es12.5)') talc + write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen + write(psb_out_unit,'("-desc asbly time : ",es12.5)') tcdasb + write(psb_out_unit,'("- mat asbly time : ",es12.5)') tasb + write(psb_out_unit,'("-total time : ",es12.5)') ttot + + end if + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + return +end subroutine psb_s_gen_prob2d diff --git a/util/psb_s_genpde_mod.f90 b/util/psb_s_genpde_mod.f90 new file mode 100644 index 00000000..8087d457 --- /dev/null +++ b/util/psb_s_genpde_mod.f90 @@ -0,0 +1,105 @@ +module psb_s_genpde_mod + + use psb_base_mod, only : psb_spk_, psb_ipk_, psb_desc_type,& + & psb_sspmat_type, psb_s_vect_type, szero + + interface + function d_func_3d(x,y,z) result(val) + import :: psb_spk_ + real(psb_spk_), intent(in) :: x,y,z + real(psb_spk_) :: val + end function d_func_3d + end interface + + interface psb_gen_prob3d + subroutine psb_s_gen_prob3d(ictxt,idim,a,bv,xv,desc_a,afmt,& + & a1,a2,a3,b1,b2,b3,c,g,info,f) + ! + ! 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. + ! + import :: psb_ipk_, psb_desc_type, psb_sspmat_type, psb_s_vect_type, d_func_3d + implicit none + procedure(d_func_3d) :: a1,a2,a3,c,b1,b2,b3,g + integer(psb_ipk_) :: idim + type(psb_sspmat_type) :: a + type(psb_s_vect_type) :: xv,bv + type(psb_desc_type) :: desc_a + integer(psb_ipk_) :: ictxt, info + character :: afmt*5 + procedure(d_func_3d), optional :: f + end subroutine psb_s_gen_prob3d + end interface + + + interface + function d_func_2d(x,y) result(val) + import :: psb_spk_ + real(psb_spk_), intent(in) :: x,y + real(psb_spk_) :: val + end function d_func_2d + end interface + + interface psb_gen_prob2d + subroutine psb_s_gen_prob2d(ictxt,idim,a,bv,xv,desc_a,afmt,a1,a2,b1,b2,c,g,info,f) + ! + ! 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. + ! + import :: psb_ipk_, psb_desc_type, psb_sspmat_type, psb_s_vect_type, d_func_2d + implicit none + procedure(d_func_2d) :: a1,a2,c,b1,b2,g + integer(psb_ipk_) :: idim + type(psb_sspmat_type) :: a + type(psb_s_vect_type) :: xv,bv + type(psb_desc_type) :: desc_a + integer(psb_ipk_) :: ictxt, info + character :: afmt*5 + procedure(d_func_2d), optional :: f + end subroutine psb_s_gen_prob2d + end interface + +contains + + function d_null_func_3d(x,y,z) result(val) + + real(psb_spk_), intent(in) :: x,y,z + real(psb_spk_) :: val + + val = szero + + end function d_null_func_3d + + function d_null_func_2d(x,y) result(val) + + real(psb_spk_), intent(in) :: x,y + real(psb_spk_) :: val + + val = szero + + end function d_null_func_2d + + +end module psb_s_genpde_mod diff --git a/util/psb_util_mod.f90 b/util/psb_util_mod.f90 index 4b78019c..f38ceac6 100644 --- a/util/psb_util_mod.f90 +++ b/util/psb_util_mod.f90 @@ -38,6 +38,7 @@ module psb_util_mod use psb_mmio_mod use psb_mat_dist_mod use psb_renum_mod - use psb_d_genmat_mod + use psb_d_genpde_mod + use psb_s_genpde_mod end module psb_util_mod