From ad88b4298d1d434994c749c8a4b787dce16b1948 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Sat, 18 Feb 2012 13:17:01 +0000 Subject: [PATCH] mld2p4-2: tests/pdegen/Makefile tests/pdegen/ppde.f90 tests/pdegen/ppde2d.f90 tests/pdegen/ppde3d.f90 tests/pdegen/spde.f90 tests/pdegen/spde2d.f90 tests/pdegen/spde3d.f90 Fixed test cases to use new support in psblas. --- tests/pdegen/Makefile | 29 +- tests/pdegen/ppde2d.f90 | 458 ++++++++++++++++++++++++++ tests/pdegen/{ppde.f90 => ppde3d.f90} | 422 ++++-------------------- tests/pdegen/spde2d.f90 | 457 +++++++++++++++++++++++++ tests/pdegen/{spde.f90 => spde3d.f90} | 423 ++++-------------------- 5 files changed, 1077 insertions(+), 712 deletions(-) create mode 100644 tests/pdegen/ppde2d.f90 rename tests/pdegen/{ppde.f90 => ppde3d.f90} (59%) create mode 100644 tests/pdegen/spde2d.f90 rename tests/pdegen/{spde.f90 => spde3d.f90} (59%) diff --git a/tests/pdegen/Makefile b/tests/pdegen/Makefile index 1453d965..fbc8eb42 100644 --- a/tests/pdegen/Makefile +++ b/tests/pdegen/Makefile @@ -11,25 +11,32 @@ FINCLUDES=$(FMFLAG). $(FMFLAG)$(MLDLIBDIR) $(FMFLAG)$(PSBINCDIR) $(FIFLAG). EXEDIR=./runs -all: spde ppde +all: spde3d ppde3d spde2d ppde2d -check: all - cd $(EXEDIR) ; ./ppde < ppde.inp +ppde3d: ppde3d.o data_input.o + $(F90LINK) ppde3d.o data_input.o -o ppde3d $(MLD_LIB) $(PSBLAS_LIB) $(LDLIBS) + /bin/mv ppde3d $(EXEDIR) -ppde: ppde.o data_input.o - $(F90LINK) ppde.o data_input.o -o ppde $(MLD_LIB) $(PSBLAS_LIB) $(LDLIBS) - /bin/mv ppde $(EXEDIR) +spde3d: spde3d.o data_input.o + $(F90LINK) spde3d.o data_input.o -o spde3d $(MLD_LIB) $(PSBLAS_LIB) $(LDLIBS) + /bin/mv spde3d $(EXEDIR) -spde: spde.o data_input.o - $(F90LINK) spde.o data_input.o -o spde $(MLD_LIB) $(PSBLAS_LIB) $(LDLIBS) - /bin/mv spde $(EXEDIR) +ppde2d: ppde2d.o data_input.o + $(F90LINK) ppde2d.o data_input.o -o ppde2d $(MLD_LIB) $(PSBLAS_LIB) $(LDLIBS) + /bin/mv ppde2d $(EXEDIR) -ppde.o spde.o: data_input.o +spde2d: spde2d.o data_input.o + $(F90LINK) spde2d.o data_input.o -o spde2d $(MLD_LIB) $(PSBLAS_LIB) $(LDLIBS) + /bin/mv spde2d $(EXEDIR) + + +ppde3d.o spde3d.o ppde2d.o spde2d.o: data_input.o clean: - /bin/rm -f data_input.o ppde.o $(EXEDIR)/ppde spde.o $(EXEDIR)/spde + /bin/rm -f data_input.o ppde3d.o spde3d.o ppde2d.o spde2d.o \ + $(EXEDIR)/ppde3d $(EXEDIR)/spde3d $(EXEDIR)/ppde2d $(EXEDIR)/spde2d verycleanlib: (cd ../..; make veryclean) diff --git a/tests/pdegen/ppde2d.f90 b/tests/pdegen/ppde2d.f90 new file mode 100644 index 00000000..69fc77e6 --- /dev/null +++ b/tests/pdegen/ppde2d.f90 @@ -0,0 +1,458 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) +!!$ +!!$ (C) Copyright 2008,2009,2010,2012 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ 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 MLD2P4 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 MLD2P4 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: ppde2d.f90 +! +! Program: ppde2d +! This sample program solves a linear system obtained by discretizing a +! PDE with Dirichlet BCs. +! +! +! The PDE is a general second order equation in 2d +! +! a1 dd(u) a2 dd(u) b1 d(u) b2 d(u) +! - ------ - ------ ----- + ------ + 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. +! +! 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 ppde2d + use psb_base_mod + use mld_prec_mod + use psb_krylov_mod + use psb_util_mod + use data_input + implicit none + + ! input parameters + character(len=20) :: kmethd, ptype + character(len=5) :: afmt + integer(psb_ipk_) :: idim + + ! miscellaneous + real(psb_dpk_), parameter :: one = 1.d0 + real(psb_dpk_) :: t1, t2, tprec + + ! sparse matrix and preconditioner + type(psb_dspmat_type) :: a + type(mld_dprec_type) :: prec + ! descriptor + type(psb_desc_type) :: desc_a + ! dense vectors + type(psb_d_vect_type) :: x,b + ! parallel environment + integer(psb_ipk_) :: ictxt, iam, np + + ! solver parameters + integer :: iter, itmax,itrace, istopc, irst, nlv + integer(psb_long_int_k_) :: amatsize, precsize, descsize + real(psb_dpk_) :: err, eps + + type precdata + character(len=20) :: descr ! verbose description of the prec + character(len=10) :: prec ! overall prectype + integer :: novr ! number of overlap layers + integer :: jsweeps ! Jacobi/smoother sweeps + character(len=16) :: restr ! restriction over application of as + character(len=16) :: prol ! prolongation over application of as + character(len=16) :: solve ! Solver type: ILU, SuperLU, UMFPACK. + integer :: fill1 ! Fill-in for factorization 1 + real(psb_dpk_) :: thr1 ! Threshold for fact. 1 ILU(T) + character(len=16) :: smther ! Smoother + integer :: nlev ! Number of levels in multilevel prec. + character(len=16) :: aggrkind ! smoothed/raw aggregatin + character(len=16) :: aggr_alg ! local or global aggregation + character(len=16) :: mltype ! additive or multiplicative 2nd level prec + character(len=16) :: smthpos ! side: pre, post, both smoothing + character(len=16) :: cmat ! coarse mat + character(len=16) :: csolve ! Coarse solver: bjac, umf, slu, sludist + character(len=16) :: csbsolve ! Coarse subsolver: ILU, ILU(T), SuperLU, UMFPACK. + integer :: cfill ! Fill-in for factorization 1 + real(psb_dpk_) :: cthres ! Threshold for fact. 1 ILU(T) + integer :: cjswp ! Jacobi sweeps + real(psb_dpk_) :: athres ! smoother aggregation threshold + end type precdata + type(precdata) :: prectype + type(psb_d_coo_sparse_mat) :: acoo + ! other variables + integer :: info, i + character(len=20) :: name,ch_err + + 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='pde2d90' + call psb_set_errverbosity(2) + ! + ! Hello world + ! + if (iam == psb_root_) then + write(*,*) 'Welcome to MLD2P4 version: ',mld_version_string_ + write(*,*) 'This is the ',trim(name),' sample program' + end if + ! + ! get parameters + ! + call get_parms(ictxt,kmethd,prectype,afmt,idim,istopc,itmax,itrace,irst,eps) + + ! + ! allocate and fill in the coefficient matrix, rhs and initial guess + ! + call psb_barrier(ictxt) + t1 = psb_wtime() + call psb_gen_pde2d(ictxt,idim,a,b,x,desc_a,afmt,a1,a2,b1,b2,c,g,info) + call psb_barrier(ictxt) + t2 = psb_wtime() - t1 + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_gen_pde2d' + 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 (psb_toupper(prectype%prec) == 'ML') then + nlv = prectype%nlev + call mld_precinit(prec,prectype%prec, info, nlev=nlv) + call mld_precset(prec,mld_smoother_type_, prectype%smther, info) + call mld_precset(prec,mld_smoother_sweeps_, prectype%jsweeps, info) + call mld_precset(prec,mld_sub_ovr_, prectype%novr, info) + call mld_precset(prec,mld_sub_restr_, prectype%restr, info) + call mld_precset(prec,mld_sub_prol_, prectype%prol, info) + call mld_precset(prec,mld_sub_solve_, prectype%solve, info) + call mld_precset(prec,mld_sub_fillin_, prectype%fill1, info) + call mld_precset(prec,mld_sub_iluthrs_, prectype%thr1, info) + call mld_precset(prec,mld_aggr_kind_, prectype%aggrkind,info) + call mld_precset(prec,mld_aggr_alg_, prectype%aggr_alg,info) + call mld_precset(prec,mld_ml_type_, prectype%mltype, info) + call mld_precset(prec,mld_smoother_pos_, prectype%smthpos, info) + if (prectype%athres >= dzero) & + & call mld_precset(prec,mld_aggr_thresh_, prectype%athres, info) + call mld_precset(prec,mld_coarse_solve_, prectype%csolve, info) + call mld_precset(prec,mld_coarse_subsolve_, prectype%csbsolve,info) + call mld_precset(prec,mld_coarse_mat_, prectype%cmat, info) + call mld_precset(prec,mld_coarse_fillin_, prectype%cfill, info) + call mld_precset(prec,mld_coarse_iluthrs_, prectype%cthres, info) + call mld_precset(prec,mld_coarse_sweeps_, prectype%cjswp, info) + else + nlv = 1 + call mld_precinit(prec,prectype%prec, info, nlev=nlv) + call mld_precset(prec,mld_smoother_sweeps_, prectype%jsweeps, info) + call mld_precset(prec,mld_sub_ovr_, prectype%novr, info) + call mld_precset(prec,mld_sub_restr_, prectype%restr, info) + call mld_precset(prec,mld_sub_prol_, prectype%prol, info) + call mld_precset(prec,mld_sub_solve_, prectype%solve, info) + call mld_precset(prec,mld_sub_fillin_, prectype%fill1, info) + call mld_precset(prec,mld_sub_iluthrs_, prectype%thr1, info) + end if + call psb_barrier(ictxt) + t1 = psb_wtime() + call mld_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 prec%dump(info,prefix='test-ml',ac=.true.,solver=.true.,smoother=.true.) + + call psb_amx(ictxt,tprec) + + if (iam == psb_root_) & + & write(psb_out_unit,'("Preconditioner time : ",es12.5)')tprec + if (iam == psb_root_) call mld_precdescr(prec,info) + 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() + call psb_krylov(kmethd,a,prec,b,x,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 DESC_A: ",i12)') descsize + write(psb_out_unit,'("Total memory occupation for PREC: ",i12)') precsize + end if + + ! + ! cleanup storage and exit + ! + call psb_gefree(b,desc_a,info) + call psb_gefree(x,desc_a,info) + call psb_spfree(a,desc_a,info) + call mld_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,prectype,afmt,idim,istopc,itmax,itrace,irst,eps) + integer :: ictxt + type(precdata) :: prectype + character(len=*) :: kmethd, afmt + integer :: idim, istopc,itmax,itrace,irst + integer :: np, iam, info + real(psb_dpk_) :: eps + character(len=20) :: buffer + + call psb_info(ictxt, iam, np) + + if (iam == psb_root_) then + call read_data(kmethd,5) + call read_data(afmt,5) + call read_data(idim,5) + call read_data(istopc,5) + call read_data(itmax,5) + call read_data(itrace,5) + call read_data(irst,5) + call read_data(eps,5) + call read_data(prectype%descr,5) ! verbose description of the prec + call read_data(prectype%prec,5) ! overall prectype + call read_data(prectype%novr,5) ! number of overlap layers + call read_data(prectype%restr,5) ! restriction over application of as + call read_data(prectype%prol,5) ! prolongation over application of as + call read_data(prectype%solve,5) ! Factorization type: ILU, SuperLU, UMFPACK. + call read_data(prectype%fill1,5) ! Fill-in for factorization 1 + call read_data(prectype%thr1,5) ! Threshold for fact. 1 ILU(T) + call read_data(prectype%jsweeps,5) ! Jacobi sweeps for PJAC + if (psb_toupper(prectype%prec) == 'ML') then + call read_data(prectype%smther,5) ! Smoother type. + call read_data(prectype%nlev,5) ! Number of levels in multilevel prec. + call read_data(prectype%aggrkind,5) ! smoothed/raw aggregatin + call read_data(prectype%aggr_alg,5) ! local or global aggregation + call read_data(prectype%mltype,5) ! additive or multiplicative 2nd level prec + call read_data(prectype%smthpos,5) ! side: pre, post, both smoothing + call read_data(prectype%cmat,5) ! coarse mat + call read_data(prectype%csolve,5) ! Factorization type: ILU, SuperLU, UMFPACK. + call read_data(prectype%csbsolve,5) ! Factorization type: ILU, SuperLU, UMFPACK. + call read_data(prectype%cfill,5) ! Fill-in for factorization 1 + call read_data(prectype%cthres,5) ! Threshold for fact. 1 ILU(T) + call read_data(prectype%cjswp,5) ! Jacobi sweeps + call read_data(prectype%athres,5) ! smoother aggr thresh + end if + end if + + ! broadcast parameters to all processors + call psb_bcast(ictxt,kmethd) + call psb_bcast(ictxt,afmt) + call psb_bcast(ictxt,idim) + call psb_bcast(ictxt,istopc) + call psb_bcast(ictxt,itmax) + call psb_bcast(ictxt,itrace) + call psb_bcast(ictxt,irst) + call psb_bcast(ictxt,eps) + + + call psb_bcast(ictxt,prectype%descr) ! verbose description of the prec + call psb_bcast(ictxt,prectype%prec) ! overall prectype + call psb_bcast(ictxt,prectype%novr) ! number of overlap layers + call psb_bcast(ictxt,prectype%restr) ! restriction over application of as + call psb_bcast(ictxt,prectype%prol) ! prolongation over application of as + call psb_bcast(ictxt,prectype%solve) ! Factorization type: ILU, SuperLU, UMFPACK. + call psb_bcast(ictxt,prectype%fill1) ! Fill-in for factorization 1 + call psb_bcast(ictxt,prectype%thr1) ! Threshold for fact. 1 ILU(T) + call psb_bcast(ictxt,prectype%jsweeps) ! Jacobi sweeps + if (psb_toupper(prectype%prec) == 'ML') then + call psb_bcast(ictxt,prectype%smther) ! Smoother type. + call psb_bcast(ictxt,prectype%nlev) ! Number of levels in multilevel prec. + call psb_bcast(ictxt,prectype%aggrkind) ! smoothed/raw aggregatin + call psb_bcast(ictxt,prectype%aggr_alg) ! local or global aggregation + call psb_bcast(ictxt,prectype%mltype) ! additive or multiplicative 2nd level prec + call psb_bcast(ictxt,prectype%smthpos) ! side: pre, post, both smoothing + call psb_bcast(ictxt,prectype%cmat) ! coarse mat + call psb_bcast(ictxt,prectype%csolve) ! Factorization type: ILU, SuperLU, UMFPACK. + call psb_bcast(ictxt,prectype%csbsolve) ! Factorization type: ILU, SuperLU, UMFPACK. + call psb_bcast(ictxt,prectype%cfill) ! Fill-in for factorization 1 + call psb_bcast(ictxt,prectype%cthres) ! Threshold for fact. 1 ILU(T) + call psb_bcast(ictxt,prectype%cjswp) ! Jacobi sweeps + call psb_bcast(ictxt,prectype%athres) ! smoother aggr thresh + end if + + if (iam == psb_root_) then + 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)') prectype%descr + write(psb_out_unit,'("Iterative method : ",a)') kmethd + write(psb_out_unit,'(" ")') + endif + + return + + end subroutine get_parms + ! + ! print an error message + ! + subroutine pr_usage(iout) + integer :: 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 + + + ! + ! functions parametrizing the differential equation + ! + function b1(x,y) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: b1 + real(psb_dpk_), intent(in) :: x,y + b1=1.d0/sqrt(2.d0) + end function b1 + function b2(x,y) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: b2 + real(psb_dpk_), intent(in) :: x,y + b2=1.d0/sqrt(2.d0) + end function b2 + function c(x,y) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: c + real(psb_dpk_), intent(in) :: x,y + c=0.d0 + end function c + function a1(x,y) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: a1 + real(psb_dpk_), intent(in) :: x,y + a1=1.d0/80 + end function a1 + function a2(x,y) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: a2 + real(psb_dpk_), intent(in) :: x,y + a2=1.d0/80 + end function a2 + function g(x,y) + use psb_base_mod, only : psb_dpk_, done + real(psb_dpk_) :: g + real(psb_dpk_), intent(in) :: x,y + g = dzero + if (x == done) then + g = done + else if (x == dzero) then + g = exp(-y**2) + end if + end function g +end program ppde2d + diff --git a/tests/pdegen/ppde.f90 b/tests/pdegen/ppde3d.f90 similarity index 59% rename from tests/pdegen/ppde.f90 rename to tests/pdegen/ppde3d.f90 index bebddeab..db96119c 100644 --- a/tests/pdegen/ppde.f90 +++ b/tests/pdegen/ppde3d.f90 @@ -4,7 +4,7 @@ !!$ MultiLevel Domain Decomposition Parallel Preconditioners Package !!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) !!$ -!!$ (C) Copyright 2008,2009,2010 +!!$ (C) Copyright 2008,2009,2010,2012 !!$ !!$ Salvatore Filippone University of Rome Tor Vergata !!$ Alfredo Buttari CNRS-IRIT, Toulouse @@ -36,39 +36,34 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -! File: ppde.f90 ! -! Program: ppde +! File: ppde3d.f90 +! +! 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 +! u = g ! -! with Dirichlet boundary conditions, on the unit cube 0<=x,y,z<=1. +! on the unit cube 0<=x,y,z<=1. ! -! Example taken from: -! C.T.Kelley -! Iterative Methods for Linear and Nonlinear Equations -! SIAM 1995 +! +! 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 mld_prec_mod use psb_krylov_mod @@ -91,8 +86,7 @@ program ppde ! descriptor type(psb_desc_type) :: desc_a ! dense matrices - type(psb_d_vect_type) :: x,b, vtst - real(psb_dpk_), allocatable :: tst(:) + type(psb_d_vect_type) :: x,b ! blacs parameters integer :: ictxt, iam, np @@ -164,7 +158,8 @@ program ppde call psb_barrier(ictxt) t1 = psb_wtime() - call create_matrix(idim,a,b,x,desc_a,ictxt,afmt,info) + call psb_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,& + & a1,a2,a3,b1,b2,b3,c,g,info) call psb_barrier(ictxt) t2 = psb_wtime() - t1 if(info /= psb_success_) then @@ -419,333 +414,60 @@ contains end subroutine pr_usage ! - ! subroutine to allocate and fill in the coefficient matrix and - ! the rhs. - ! - subroutine create_matrix(idim,a,b,xv,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 - implicit none - integer :: idim - integer, parameter :: nb=20 - type(psb_d_vect_type) :: b,xv - type(psb_desc_type) :: desc_a - integer :: ictxt, info - character :: afmt*5 - type(psb_dspmat_type) :: a - real(psb_dpk_) :: zt(nb),x,y,z - integer :: m,n,nnz,glob_row,nlr,i,ii,ib,k - integer :: ix,iy,iz,ia,indx_owner - integer :: np, iam, nr, nt - integer :: element - integer, allocatable :: irow(:),icol(:),myidx(:) - real(psb_dpk_), allocatable :: val(:) - ! deltah dimension of each grid cell - ! deltat discretization time - real(psb_dpk_) :: deltah, sqdeltah, deltah2 - real(psb_dpk_), parameter :: rhs=0.d0,one=1.d0,zero=0.d0 - real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen - real(psb_dpk_) :: a1, a2, a3, a4, b1, b2, b3 - external :: a1, a2, a3, a4, b1, b2, b3 - integer :: err_act - - character(len=20) :: name, ch_err - - info = psb_success_ - name = 'create_matrix' - call psb_erractionsave(err_act) - - call psb_info(ictxt, iam, np) - - deltah = 1.d0/(idim-1) - 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(b,desc_a,info) - if (info == psb_success_) call psb_geall(xv,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),b,desc_a,info) - if(info /= psb_success_) exit - zt(:)=0.d0 - call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info) - if(info /= psb_success_) exit - end do - - tgen = psb_wtime()-t1 - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - 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 - call psb_geasb(b,desc_a,info) - call psb_geasb(xv,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 - ch_err = a%get_fmt() - write(psb_out_unit,& - & '("The matrix has been generated and assembled in ",a3," format.")')& - & ch_err(1:3) - 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 + ! functions parametrizing the differential equation + ! + function b1(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: b1 + real(psb_dpk_), intent(in) :: x,y,z + b1=1.d0/sqrt(3.d0) + end function b1 + function b2(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: b2 + real(psb_dpk_), intent(in) :: x,y,z + b2=1.d0/sqrt(3.d0) + end function b2 + function b3(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: b3 + real(psb_dpk_), intent(in) :: x,y,z + b3=1.d0/sqrt(3.d0) + end function b3 + function c(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: c + real(psb_dpk_), intent(in) :: x,y,z + c=0.d0 + end function c + function a1(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: a1 + real(psb_dpk_), intent(in) :: x,y,z + a1=1.d0/80 + end function a1 + function a2(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: a2 + real(psb_dpk_), intent(in) :: x,y,z + a2=1.d0/80 + end function a2 + function a3(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: a3 + real(psb_dpk_), intent(in) :: x,y,z + a3=1.d0/80 + end function a3 + function g(x,y,z) + use psb_base_mod, only : psb_dpk_, done + real(psb_dpk_) :: g + real(psb_dpk_), intent(in) :: x,y,z + g = dzero + if (x == done) then + g = done + else if (x == dzero) then + g = exp(y**2-z**2) end if - return - end subroutine create_matrix -end program ppde -! -! functions parametrizing the differential equation -! -function a1(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: a1 - real(psb_dpk_) :: x,y,z - a1=1.d0 -end function a1 -function a2(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: a2 - real(psb_dpk_) :: x,y,z - a2=2.d1*y -end function a2 -function a3(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: a3 - real(psb_dpk_) :: x,y,z - a3=1.d0 -end function a3 -function a4(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: a4 - real(psb_dpk_) :: x,y,z - a4=1.d0 -end function a4 -function b1(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: b1 - real(psb_dpk_) :: x,y,z - b1=1.d0 -end function b1 -function b2(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: b2 - real(psb_dpk_) :: x,y,z - b2=1.d0 -end function b2 -function b3(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: b3 - real(psb_dpk_) :: x,y,z - b3=1.d0 -end function b3 - + end function g +end program ppde3d diff --git a/tests/pdegen/spde2d.f90 b/tests/pdegen/spde2d.f90 new file mode 100644 index 00000000..aca14652 --- /dev/null +++ b/tests/pdegen/spde2d.f90 @@ -0,0 +1,457 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) +!!$ +!!$ (C) Copyright 2008,2009,2010,2012 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ 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 MLD2P4 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 MLD2P4 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: spde2d.f90 +! +! Program: spde2d +! This sample program solves a linear system obtained by discretizing a +! PDE with Dirichlet BCs. +! +! +! The PDE is a general second order equation in 2d +! +! a1 dd(u) a2 dd(u) b1 d(u) b2 d(u) +! - ------ - ------ ----- + ------ + 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. +! +! 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 spde2d + use psb_base_mod + use mld_prec_mod + use psb_krylov_mod + use psb_util_mod + use data_input + implicit none + + ! input parameters + character(len=20) :: kmethd, ptype + character(len=5) :: afmt + integer :: 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(mld_sprec_type) :: prec + ! descriptor + type(psb_desc_type) :: desc_a + ! dense matrices + type(psb_s_vect_type) :: x,b + ! blacs parameters + integer :: ictxt, iam, np + + ! solver parameters + integer :: iter, itmax,itrace, istopc, irst, nlv + integer(psb_long_int_k_) :: amatsize, precsize, descsize + real(psb_spk_) :: err, eps + + type precdata + character(len=20) :: descr ! verbose description of the prec + character(len=10) :: prec ! overall prectype + integer :: novr ! number of overlap layers + integer :: jsweeps ! Jacobi/smoother sweeps + character(len=16) :: restr ! restriction over application of as + character(len=16) :: prol ! prolongation over application of as + character(len=16) :: solve ! Solver type: ILU, SuperLU, UMFPACK. + integer :: fill1 ! Fill-in for factorization 1 + real(psb_spk_) :: thr1 ! Threshold for fact. 1 ILU(T) + character(len=16) :: smther ! Smoother + integer :: nlev ! Number of levels in multilevel prec. + character(len=16) :: aggrkind ! smoothed/raw aggregatin + character(len=16) :: aggr_alg ! local or global aggregation + character(len=16) :: mltype ! additive or multiplicative 2nd level prec + character(len=16) :: smthpos ! side: pre, post, both smoothing + character(len=16) :: cmat ! coarse mat + character(len=16) :: csolve ! Coarse solver: bjac, umf, slu, sludist + character(len=16) :: csbsolve ! Coarse subsolver: ILU, ILU(T), SuperLU, UMFPACK. + integer :: cfill ! Fill-in for factorization 1 + real(psb_spk_) :: cthres ! Threshold for fact. 1 ILU(T) + integer :: cjswp ! Jacobi sweeps + real(psb_spk_) :: athres ! smoother aggregation threshold + end type precdata + type(precdata) :: prectype + type(psb_s_coo_sparse_mat) :: acoo + ! other variables + integer :: info + character(len=20) :: name,ch_err + + 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='pde2d90' + call psb_set_errverbosity(2) + ! + ! Hello world + ! + if (iam == psb_root_) then + write(*,*) 'Welcome to MLD2P4 version: ',mld_version_string_ + write(*,*) 'This is the ',trim(name),' sample program' + end if + ! + ! get parameters + ! + call get_parms(ictxt,kmethd,prectype,afmt,idim,istopc,itmax,itrace,irst,eps) + + ! + ! allocate and fill in the coefficient matrix, rhs and initial guess + ! + call psb_barrier(ictxt) + t1 = psb_wtime() + call psb_gen_pde2d(ictxt,idim,a,b,x,desc_a,afmt,a1,a2,b1,b2,c,g,info) + call psb_barrier(ictxt) + t2 = psb_wtime() - t1 + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_gen_pde2d' + 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 (psb_toupper(prectype%prec) == 'ML') then + nlv = prectype%nlev + call mld_precinit(prec,prectype%prec, info, nlev=nlv) + call mld_precset(prec,mld_smoother_type_, prectype%smther, info) + call mld_precset(prec,mld_smoother_sweeps_, prectype%jsweeps, info) + call mld_precset(prec,mld_sub_ovr_, prectype%novr, info) + call mld_precset(prec,mld_sub_restr_, prectype%restr, info) + call mld_precset(prec,mld_sub_prol_, prectype%prol, info) + call mld_precset(prec,mld_sub_solve_, prectype%solve, info) + call mld_precset(prec,mld_sub_fillin_, prectype%fill1, info) + call mld_precset(prec,mld_sub_iluthrs_, prectype%thr1, info) + call mld_precset(prec,mld_aggr_kind_, prectype%aggrkind,info) + call mld_precset(prec,mld_aggr_alg_, prectype%aggr_alg,info) + call mld_precset(prec,mld_ml_type_, prectype%mltype, info) + call mld_precset(prec,mld_smoother_pos_, prectype%smthpos, info) + if (prectype%athres >= szero) & + & call mld_precset(prec,mld_aggr_thresh_, prectype%athres, info) + call mld_precset(prec,mld_coarse_solve_, prectype%csolve, info) + call mld_precset(prec,mld_coarse_subsolve_, prectype%csbsolve,info) + call mld_precset(prec,mld_coarse_mat_, prectype%cmat, info) + call mld_precset(prec,mld_coarse_fillin_, prectype%cfill, info) + call mld_precset(prec,mld_coarse_iluthrs_, prectype%cthres, info) + call mld_precset(prec,mld_coarse_sweeps_, prectype%cjswp, info) + else + nlv = 1 + call mld_precinit(prec,prectype%prec, info, nlev=nlv) + call mld_precset(prec,mld_smoother_sweeps_, prectype%jsweeps, info) + call mld_precset(prec,mld_sub_ovr_, prectype%novr, info) + call mld_precset(prec,mld_sub_restr_, prectype%restr, info) + call mld_precset(prec,mld_sub_prol_, prectype%prol, info) + call mld_precset(prec,mld_sub_solve_, prectype%solve, info) + call mld_precset(prec,mld_sub_fillin_, prectype%fill1, info) + call mld_precset(prec,mld_sub_iluthrs_, prectype%thr1, info) + end if + call psb_barrier(ictxt) + t1 = psb_wtime() + call mld_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 prec%dump(info,prefix='test-ml',ac=.true.,solver=.true.,smoother=.true.) + + call psb_amx(ictxt,tprec) + + if (iam == psb_root_) & + & write(psb_out_unit,'("Preconditioner time : ",es12.5)')tprec + if (iam == psb_root_) call mld_precdescr(prec,info) + 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() + call psb_krylov(kmethd,a,prec,b,x,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 DESC_A: ",i12)') descsize + write(psb_out_unit,'("Total memory occupation for PREC: ",i12)') precsize + end if + + ! + ! cleanup storage and exit + ! + call psb_gefree(b,desc_a,info) + call psb_gefree(x,desc_a,info) + call psb_spfree(a,desc_a,info) + call mld_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,prectype,afmt,idim,istopc,itmax,itrace,irst,eps) + integer :: ictxt + type(precdata) :: prectype + character(len=*) :: kmethd, afmt + integer :: idim, istopc,itmax,itrace,irst + integer :: np, iam, info + real(psb_spk_) :: eps + character(len=20) :: buffer + + call psb_info(ictxt, iam, np) + + if (iam == psb_root_) then + call read_data(kmethd,5) + call read_data(afmt,5) + call read_data(idim,5) + call read_data(istopc,5) + call read_data(itmax,5) + call read_data(itrace,5) + call read_data(irst,5) + call read_data(eps,5) + call read_data(prectype%descr,5) ! verbose description of the prec + call read_data(prectype%prec,5) ! overall prectype + call read_data(prectype%novr,5) ! number of overlap layers + call read_data(prectype%restr,5) ! restriction over application of as + call read_data(prectype%prol,5) ! prolongation over application of as + call read_data(prectype%solve,5) ! Factorization type: ILU, SuperLU, UMFPACK. + call read_data(prectype%fill1,5) ! Fill-in for factorization 1 + call read_data(prectype%thr1,5) ! Threshold for fact. 1 ILU(T) + call read_data(prectype%jsweeps,5) ! Jacobi sweeps for PJAC + if (psb_toupper(prectype%prec) == 'ML') then + call read_data(prectype%smther,5) ! Smoother type. + call read_data(prectype%nlev,5) ! Number of levels in multilevel prec. + call read_data(prectype%aggrkind,5) ! smoothed/raw aggregatin + call read_data(prectype%aggr_alg,5) ! local or global aggregation + call read_data(prectype%mltype,5) ! additive or multiplicative 2nd level prec + call read_data(prectype%smthpos,5) ! side: pre, post, both smoothing + call read_data(prectype%cmat,5) ! coarse mat + call read_data(prectype%csolve,5) ! Factorization type: ILU, SuperLU, UMFPACK. + call read_data(prectype%csbsolve,5) ! Factorization type: ILU, SuperLU, UMFPACK. + call read_data(prectype%cfill,5) ! Fill-in for factorization 1 + call read_data(prectype%cthres,5) ! Threshold for fact. 1 ILU(T) + call read_data(prectype%cjswp,5) ! Jacobi sweeps + call read_data(prectype%athres,5) ! smoother aggr thresh + end if + end if + + ! broadcast parameters to all processors + call psb_bcast(ictxt,kmethd) + call psb_bcast(ictxt,afmt) + call psb_bcast(ictxt,idim) + call psb_bcast(ictxt,istopc) + call psb_bcast(ictxt,itmax) + call psb_bcast(ictxt,itrace) + call psb_bcast(ictxt,irst) + call psb_bcast(ictxt,eps) + + + call psb_bcast(ictxt,prectype%descr) ! verbose description of the prec + call psb_bcast(ictxt,prectype%prec) ! overall prectype + call psb_bcast(ictxt,prectype%novr) ! number of overlap layers + call psb_bcast(ictxt,prectype%restr) ! restriction over application of as + call psb_bcast(ictxt,prectype%prol) ! prolongation over application of as + call psb_bcast(ictxt,prectype%solve) ! Factorization type: ILU, SuperLU, UMFPACK. + call psb_bcast(ictxt,prectype%fill1) ! Fill-in for factorization 1 + call psb_bcast(ictxt,prectype%thr1) ! Threshold for fact. 1 ILU(T) + call psb_bcast(ictxt,prectype%jsweeps) ! Jacobi sweeps + if (psb_toupper(prectype%prec) == 'ML') then + call psb_bcast(ictxt,prectype%smther) ! Smoother type. + call psb_bcast(ictxt,prectype%nlev) ! Number of levels in multilevel prec. + call psb_bcast(ictxt,prectype%aggrkind) ! smoothed/raw aggregatin + call psb_bcast(ictxt,prectype%aggr_alg) ! local or global aggregation + call psb_bcast(ictxt,prectype%mltype) ! additive or multiplicative 2nd level prec + call psb_bcast(ictxt,prectype%smthpos) ! side: pre, post, both smoothing + call psb_bcast(ictxt,prectype%cmat) ! coarse mat + call psb_bcast(ictxt,prectype%csolve) ! Factorization type: ILU, SuperLU, UMFPACK. + call psb_bcast(ictxt,prectype%csbsolve) ! Factorization type: ILU, SuperLU, UMFPACK. + call psb_bcast(ictxt,prectype%cfill) ! Fill-in for factorization 1 + call psb_bcast(ictxt,prectype%cthres) ! Threshold for fact. 1 ILU(T) + call psb_bcast(ictxt,prectype%cjswp) ! Jacobi sweeps + call psb_bcast(ictxt,prectype%athres) ! smoother aggr thresh + end if + + if (iam == psb_root_) then + 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)') prectype%descr + write(psb_out_unit,'("Iterative method : ",a)') kmethd + write(psb_out_unit,'(" ")') + endif + + return + + end subroutine get_parms + ! + ! print an error message + ! + subroutine pr_usage(iout) + integer :: 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 + + ! + ! functions parametrizing the differential equation + ! + function b1(x,y) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: b1 + real(psb_spk_), intent(in) :: x,y + b1=1.e0/sqrt(2.e0) + end function b1 + function b2(x,y) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: b2 + real(psb_spk_), intent(in) :: x,y + b2=1.e0/sqrt(2.e0) + end function b2 + function c(x,y) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: c + real(psb_spk_), intent(in) :: x,y + c=0.e0 + end function c + function a1(x,y) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: a1 + real(psb_spk_), intent(in) :: x,y + a1=1.e0/80 + end function a1 + function a2(x,y) + use psb_base_mod, only : psb_spk_ + real(psb_spk_) :: a2 + real(psb_spk_), intent(in) :: x,y + a2=1.e0/80 + end function a2 + function g(x,y) + use psb_base_mod, only : psb_spk_, sone + real(psb_spk_) :: g + real(psb_spk_), intent(in) :: x,y + g = szero + if (x == sone) then + g = sone + else if (x == szero) then + g = exp(-y**2) + end if + end function g + +end program spde2d diff --git a/tests/pdegen/spde.f90 b/tests/pdegen/spde3d.f90 similarity index 59% rename from tests/pdegen/spde.f90 rename to tests/pdegen/spde3d.f90 index 3a886a64..32b12f88 100644 --- a/tests/pdegen/spde.f90 +++ b/tests/pdegen/spde3d.f90 @@ -4,7 +4,7 @@ !!$ MultiLevel Domain Decomposition Parallel Preconditioners Package !!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) !!$ -!!$ (C) Copyright 2008,2009,2010 +!!$ (C) Copyright 2008,2009,2010,2012 !!$ !!$ Salvatore Filippone University of Rome Tor Vergata !!$ Alfredo Buttari CNRS-IRIT, Toulouse @@ -36,39 +36,34 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -! File: ppde.f90 ! -! Program: ppde +! 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 ! -! 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 +! u = g ! -! with Dirichlet boundary conditions, on the unit cube 0<=x,y,z<=1. +! on the unit cube 0<=x,y,z<=1. ! -! Example taken from: -! C.T.Kelley -! Iterative Methods for Linear and Nonlinear Equations -! SIAM 1995 +! +! 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 spde +program spde3d use psb_base_mod use mld_prec_mod use psb_krylov_mod @@ -91,8 +86,7 @@ program spde ! descriptor type(psb_desc_type) :: desc_a ! dense matrices - type(psb_s_vect_type) :: x,b, vtst - real(psb_spk_), allocatable :: tst(:) + type(psb_s_vect_type) :: x,b ! blacs parameters integer :: ictxt, iam, np @@ -164,7 +158,8 @@ program spde call psb_barrier(ictxt) t1 = psb_wtime() - call create_matrix(idim,a,b,x,desc_a,ictxt,afmt,info) + call psb_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,& + & a1,a2,a3,b1,b2,b3,c,g,info) call psb_barrier(ictxt) t2 = psb_wtime() - t1 if(info /= psb_success_) then @@ -417,335 +412,61 @@ contains 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,b,xv,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 - implicit none - integer :: idim - integer, parameter :: nb=20 - type(psb_s_vect_type) :: b,xv - type(psb_desc_type) :: desc_a - integer :: ictxt, info - character :: afmt*5 - type(psb_sspmat_type) :: a - real(psb_spk_) :: zt(nb),x,y,z - integer :: m,n,nnz,glob_row,nlr,i,ii,ib,k - integer :: ix,iy,iz,ia,indx_owner - integer :: np, iam, nr, nt - integer :: element - integer, allocatable :: irow(:),icol(:),myidx(:) - real(psb_spk_), allocatable :: val(:) - ! deltah dimension of each grid cell - ! deltat discretization time - real(psb_spk_) :: deltah, 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 :: err_act - - character(len=20) :: name, ch_err - - info = psb_success_ - name = 'create_matrix' - call psb_erractionsave(err_act) - - call psb_info(ictxt, iam, np) - - deltah = 1.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(b,desc_a,info) - if (info == psb_success_) call psb_geall(xv,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),b,desc_a,info) - if(info /= psb_success_) exit - zt(:)=0.d0 - call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info) - if(info /= psb_success_) exit - end do - - tgen = psb_wtime()-t1 - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - 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 - call psb_geasb(b,desc_a,info) - call psb_geasb(xv,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 - ch_err = a%get_fmt() - write(psb_out_unit,& - & '("The matrix has been generated and assembled in ",a3," format.")')& - & ch_err(1:3) - 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 + ! 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 - return - end subroutine create_matrix -end program spde -! -! functions parametrizing the differential equation -! -function a1(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: a1 - real(psb_spk_) :: x,y,z - a1=1.e0 -end function a1 -function a2(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: a2 - real(psb_spk_) :: x,y,z - a2=2.e1*y -end function a2 -function a3(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: a3 - real(psb_spk_) :: x,y,z - a3=1.e0 -end function a3 -function a4(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: a4 - real(psb_spk_) :: x,y,z - a4=1.e0 -end function a4 -function b1(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: b1 - real(psb_spk_) :: x,y,z - b1=1.e0 -end function b1 -function b2(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: b2 - real(psb_spk_) :: x,y,z - b2=1.e0 -end function b2 -function b3(x,y,z) - use psb_base_mod, only : psb_spk_ - real(psb_spk_) :: b3 - real(psb_spk_) :: x,y,z - b3=1.e0 -end function b3 - + end function g +end program spde3d