diff --git a/samples/newslv/Makefile b/samples/newslv/Makefile new file mode 100644 index 00000000..72b9fc22 --- /dev/null +++ b/samples/newslv/Makefile @@ -0,0 +1,34 @@ +AMGDIR=../.. +AMGINCDIR=$(AMGDIR)/include +include $(AMGINCDIR)/Make.inc.amg4psblas +AMGMODDIR=$(AMGDIR)/modules +AMGLIBDIR=$(AMGDIR)/lib +AMG_LIBS=-L$(AMGLIBDIR) -lpsb_linsolve -lamg_prec -lpsb_prec +FINCLUDES=$(FMFLAG). $(FMFLAG)$(AMGMODDIR) $(FMFLAG)$(AMGINCDIR) $(PSBLAS_INCLUDES) $(FIFLAG). + +PD3DOBJS=amg_pde3d_newslv.o data_input.o amg_d_tlu_solver.o amg_d_tlu_solver_impl.o +PSOBJS=spde.o data_input.o +EXEDIR=./runs +LINKOPT= + +all: amg_pde3d_newslv + + +amg_pde3d_newslv: $(PD3DOBJS) + $(FLINK) $(LINKOPT) $(PD3DOBJS) -o amg_pde3d_newslv $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS) + /bin/mv amg_pde3d_newslv $(EXEDIR) + + +amg_pde3d_newslv.o amg_d_tlu_solver_impl.o: data_input.o amg_d_tlu_solver.o + + +clean: + /bin/rm -f $(PD3DOBJS) $(EXEDIR)/amg_pde3d_newslv + +verycleanlib: + (cd ../..; make veryclean) +lib: + (cd ../../; make library) + + + diff --git a/tests/newslv/amg_d_tlu_solver.f90 b/samples/newslv/amg_d_tlu_solver.f90 similarity index 71% rename from tests/newslv/amg_d_tlu_solver.f90 rename to samples/newslv/amg_d_tlu_solver.f90 index af6fb6bd..c357cdeb 100644 --- a/tests/newslv/amg_d_tlu_solver.f90 +++ b/samples/newslv/amg_d_tlu_solver.f90 @@ -1,15 +1,15 @@ -! -! -! MLD2P4 version 2.2 -! MultiLevel Domain Decomposition Parallel Preconditioners Package -! based on PSBLAS (Parallel Sparse BLAS version 3.5) -! -! (C) Copyright 2008-2018 -! -! Salvatore Filippone -! Pasqua D'Ambra -! Daniela di Serafino -! +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -18,14 +18,14 @@ ! 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 +! 3. The name of the AMG4PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS 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 @@ -33,23 +33,22 @@ ! 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: mld_d_tlu_solver_mod.f90 ! -! Module: mld_d_tlu_solver_mod +! File: amg_d_tlu_solver_mod.f90 +! +! Module: amg_d_tlu_solver_mod ! ! This module serves as an example of how to define a new solver and integrate -! it in MLD2P4 via the P%SET(sv,info) method. +! it in AMG4PSBLAS via the P%SET(sv,info) method. ! In this example we are extending the ILU solver by implementing a new factorization algorithm. ! In actual reality, we are just giving a new name to ILU(0), but this should be sufficient to show ! the basics. ! ! The code is divided in two files: ! 1. The interface file (this one) -! 2. The implementation file (mld_d_tlu_solver_impl.f90) +! 2. The implementation file (amg_d_tlu_solver_impl.f90) ! ! The separation between interface and implementation is an essential part of the ! object-oriented design. The most appropriate tool would be to have the implementation @@ -59,12 +58,12 @@ ! ! -module mld_d_tlu_solver +module amg_d_tlu_solver - use mld_d_ilu_solver - ! use mld_d_ilu_fact_mod + use amg_d_ilu_solver + ! use amg_d_ilu_fact_mod - type, extends(mld_d_ilu_solver_type) :: mld_d_tlu_solver_type + type, extends(amg_d_ilu_solver_type) :: amg_d_tlu_solver_type ! ! These are already defined in the ILU solver type; since we ! are supposedly implementing a new factorization strategy, the @@ -82,12 +81,12 @@ module mld_d_tlu_solver ! in common among all possible ILU factorizations ! ! - !procedure, pass(sv) :: dump => mld_d_tlu_solver_dmp + !procedure, pass(sv) :: dump => amg_d_tlu_solver_dmp !procedure, pass(sv) :: ccheck => d_tlu_solver_check - !procedure, pass(sv) :: clone => mld_d_tlu_solver_clone - !procedure, pass(sv) :: cnv => mld_d_tlu_solver_cnv - !procedure, pass(sv) :: apply_v => mld_d_tlu_solver_apply_vect - !procedure, pass(sv) :: apply_a => mld_d_tlu_solver_apply + !procedure, pass(sv) :: clone => amg_d_tlu_solver_clone + !procedure, pass(sv) :: cnv => amg_d_tlu_solver_cnv + !procedure, pass(sv) :: apply_v => amg_d_tlu_solver_apply_vect + !procedure, pass(sv) :: apply_a => amg_d_tlu_solver_apply !procedure, pass(sv) :: free => d_tlu_solver_free !procedure, pass(sv) :: seti => d_tlu_solver_seti !procedure, pass(sv) :: setc => d_tlu_solver_setc @@ -106,28 +105,28 @@ module mld_d_tlu_solver ! procedure, pass(sv) :: descr => d_tlu_solver_descr procedure, pass(sv) :: default => d_tlu_solver_default - procedure, pass(sv) :: build => mld_d_tlu_solver_bld + procedure, pass(sv) :: build => amg_d_tlu_solver_bld procedure, nopass :: get_fmt => d_tlu_solver_get_fmt - end type mld_d_tlu_solver_type + end type amg_d_tlu_solver_type private :: d_tlu_solver_get_fmt, d_tlu_solver_descr, d_tlu_solver_default interface - subroutine mld_d_tlu_solver_bld(a,desc_a,sv,info,b,amold,vmold, imold) - import :: psb_desc_type, mld_d_tlu_solver_type, psb_d_vect_type, psb_dpk_, & + subroutine amg_d_tlu_solver_bld(a,desc_a,sv,info,b,amold,vmold, imold) + import :: psb_desc_type, amg_d_tlu_solver_type, psb_d_vect_type, psb_dpk_, & & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,& & psb_ipk_, psb_i_base_vect_type implicit none type(psb_dspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(mld_d_tlu_solver_type), intent(inout) :: sv + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_d_tlu_solver_type), intent(inout) :: sv integer(psb_ipk_), intent(out) :: info type(psb_dspmat_type), intent(in), target, optional :: b class(psb_d_base_sparse_mat), intent(in), optional :: amold class(psb_d_base_vect_type), intent(in), optional :: vmold class(psb_i_base_vect_type), intent(in), optional :: imold - end subroutine mld_d_tlu_solver_bld + end subroutine amg_d_tlu_solver_bld end interface contains @@ -141,9 +140,9 @@ contains Implicit None ! Arguments - class(mld_d_tlu_solver_type), intent(inout) :: sv + class(amg_d_tlu_solver_type), intent(inout) :: sv - sv%fact_type = mld_ilu_n_ + sv%fact_type = amg_ilu_n_ sv%fill_in = 0 sv%thresh = dzero @@ -157,20 +156,22 @@ contains val = "TLU solver" end function d_tlu_solver_get_fmt - subroutine d_tlu_solver_descr(sv,info,iout,coarse) + subroutine d_tlu_solver_descr(sv,info,iout,coarse,prefix) Implicit None ! Arguments - class(mld_d_tlu_solver_type), intent(in) :: sv + class(amg_d_tlu_solver_type), intent(in) :: sv integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: iout logical, intent(in), optional :: coarse + character(len=*), intent(in), optional :: prefix ! Local variables integer(psb_ipk_) :: err_act - character(len=20), parameter :: name='mld_d_tlu_solver_descr' + character(len=20), parameter :: name='amg_d_tlu_solver_descr' integer(psb_ipk_) :: iout_ + character(1024) :: prefix_ call psb_erractionsave(err_act) info = psb_success_ @@ -179,14 +180,20 @@ contains else iout_ = psb_out_unit endif + if (present(prefix)) then + prefix_ = prefix + else + prefix_ = "" + end if + - write(iout_,*) ' Incomplete factorization solver: New Factorization TLU ' + write(iout_,*) trim(prefix_), ' Incomplete factorization solver: New Factorization TLU ' select case(sv%fact_type) - case(mld_ilu_n_,mld_milu_n_) - write(iout_,*) ' Fill level:',sv%fill_in - case(mld_ilu_t_) - write(iout_,*) ' Fill level:',sv%fill_in - write(iout_,*) ' Fill threshold :',sv%thresh + case(amg_ilu_n_,amg_milu_n_) + write(iout_,*) trim(prefix_), ' Fill level:',sv%fill_in + case(amg_ilu_t_) + write(iout_,*) trim(prefix_), ' Fill level:',sv%fill_in + write(iout_,*) trim(prefix_), ' Fill threshold :',sv%thresh end select call psb_erractionrestore(err_act) @@ -196,4 +203,4 @@ contains return end subroutine d_tlu_solver_descr -end module mld_d_tlu_solver +end module amg_d_tlu_solver diff --git a/tests/newslv/amg_d_tlu_solver_impl.f90 b/samples/newslv/amg_d_tlu_solver_impl.f90 similarity index 79% rename from tests/newslv/amg_d_tlu_solver_impl.f90 rename to samples/newslv/amg_d_tlu_solver_impl.f90 index 209e85e4..6593ce9a 100644 --- a/tests/newslv/amg_d_tlu_solver_impl.f90 +++ b/samples/newslv/amg_d_tlu_solver_impl.f90 @@ -1,15 +1,15 @@ -! -! -! MLD2P4 version 2.2 -! MultiLevel Domain Decomposition Parallel Preconditioners Package -! based on PSBLAS (Parallel Sparse BLAS version 3.5) -! -! (C) Copyright 2008-2018 -! -! Salvatore Filippone -! Pasqua D'Ambra -! Daniela di Serafino -! +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -18,14 +18,14 @@ ! 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 +! 3. The name of the AMG4PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS 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 @@ -33,36 +33,36 @@ ! 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: mld_d_tlu_solver_impl.f90 +! File: amg_d_tlu_solver_impl.f90 ! -! This is the implementation file corresponding to mld_d_tlu_solver_mod. +! This is the implementation file corresponding to amg_d_tlu_solver_mod. ! ! In this example we are extending the ILU solver; we pretend to have a new ! factorization method, but since we are only interested in the interfacing, ! we are simply giving a new name to ILU(0). ! ! -subroutine mld_d_tlu_solver_bld(a,desc_a,sv,info,b,amold,vmold) +subroutine amg_d_tlu_solver_bld(a,desc_a,sv,info,b,amold,vmold) use psb_base_mod - use mld_d_tlu_solver, mld_protect_name => mld_d_tlu_solver_bld - use mld_d_ilu_fact_mod + use amg_d_tlu_solver, amg_protect_name => amg_d_tlu_solver_bld Implicit None ! Arguments type(psb_dspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(mld_d_tlu_solver_type), intent(inout) :: sv + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_d_tlu_solver_type), intent(inout) :: sv integer, intent(out) :: info type(psb_dspmat_type), intent(in), target, optional :: b class(psb_d_base_sparse_mat), intent(in), optional :: amold class(psb_d_base_vect_type), intent(in), optional :: vmold ! Local variables integer :: n_row,n_col, nrow_a, nztota - integer :: ctxt,np,me,i, err_act, debug_unit, debug_level + integer :: np,me,i, err_act, debug_unit, debug_level + type(psb_ctxt_type) :: ctxt character(len=20) :: name='d_tlu_solver_bld', ch_err info=psb_success_ @@ -105,7 +105,7 @@ subroutine mld_d_tlu_solver_bld(a,desc_a,sv,info,b,amold,vmold) ! Fill-in 0, simple implementation. - call mld_ilu0_fact(sv%fact_type,a,sv%l,sv%u,& + call psb_ilu0_fact(sv%fact_type,a,sv%l,sv%u,& & sv%d,info,blck=b) call sv%l%set_asb() @@ -128,5 +128,5 @@ subroutine mld_d_tlu_solver_bld(a,desc_a,sv,info,b,amold,vmold) 9999 call psb_error_handler(err_act) return -end subroutine mld_d_tlu_solver_bld +end subroutine amg_d_tlu_solver_bld diff --git a/samples/newslv/amg_pde3d_newslv.f90 b/samples/newslv/amg_pde3d_newslv.f90 new file mode 100644 index 00000000..52902a13 --- /dev/null +++ b/samples/newslv/amg_pde3d_newslv.f90 @@ -0,0 +1,1324 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 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 AMG4PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS 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. +! +module amg_d_pde3d_poisson_mod + use psb_base_mod, only : psb_dpk_, done, dzero + real(psb_dpk_), save, private :: epsilon=done/80 +contains + subroutine pde_set_parm3d_poisson(dat) + real(psb_dpk_), intent(in) :: dat + epsilon = dat + end subroutine pde_set_parm3d_poisson + ! + ! functions parametrizing the differential equation + ! + function b1_poisson(x,y,z) + implicit none + real(psb_dpk_) :: b1_poisson + real(psb_dpk_), intent(in) :: x,y,z + b1_poisson=dzero + end function b1_poisson + function b2_poisson(x,y,z) + implicit none + real(psb_dpk_) :: b2_poisson + real(psb_dpk_), intent(in) :: x,y,z + b2_poisson=dzero + end function b2_poisson + function b3_poisson(x,y,z) + implicit none + real(psb_dpk_) :: b3_poisson + real(psb_dpk_), intent(in) :: x,y,z + b3_poisson=dzero + end function b3_poisson + function c_poisson(x,y,z) + implicit none + real(psb_dpk_) :: c_poisson + real(psb_dpk_), intent(in) :: x,y,z + c_poisson=dzero + end function c_poisson + function a1_poisson(x,y,z) + implicit none + real(psb_dpk_) :: a1_poisson + real(psb_dpk_), intent(in) :: x,y,z + a1_poisson=epsilon + end function a1_poisson + function a2_poisson(x,y,z) + implicit none + real(psb_dpk_) :: a2_poisson + real(psb_dpk_), intent(in) :: x,y,z + a2_poisson=epsilon + end function a2_poisson + function a3_poisson(x,y,z) + implicit none + real(psb_dpk_) :: a3_poisson + real(psb_dpk_), intent(in) :: x,y,z + a3_poisson=epsilon + end function a3_poisson + function g_poisson(x,y,z) + implicit none + real(psb_dpk_) :: g_poisson + real(psb_dpk_), intent(in) :: x,y,z + g_poisson = dzero + if (x == done) then + g_poisson = done + else if (x == dzero) then + g_poisson = done + end if + end function g_poisson +end module amg_d_pde3d_poisson_mod +module amg_d_genpde_mod + + use psb_base_mod, only : psb_dpk_, psb_ipk_, psb_desc_type,& + & psb_dspmat_type, psb_d_vect_type, dzero, done,& + & psb_d_base_sparse_mat, psb_d_base_vect_type, psb_i_base_vect_type + + interface + function d_func_3d(x,y,z) result(val) + import :: psb_dpk_ + real(psb_dpk_), intent(in) :: x,y,z + real(psb_dpk_) :: val + end function d_func_3d + end interface + + interface amg_gen_pde3d + module procedure amg_d_gen_pde3d + end interface amg_gen_pde3d + + interface + function d_func_2d(x,y) result(val) + import :: psb_dpk_ + real(psb_dpk_), intent(in) :: x,y + real(psb_dpk_) :: val + end function d_func_2d + end interface + + interface amg_gen_pde2d + module procedure amg_d_gen_pde2d + end interface amg_gen_pde2d + +contains + + function d_null_func_2d(x,y) result(val) + + real(psb_dpk_), intent(in) :: x,y + real(psb_dpk_) :: val + + val = dzero + + end function d_null_func_2d + + function d_null_func_3d(x,y,z) result(val) + + real(psb_dpk_), intent(in) :: x,y,z + real(psb_dpk_) :: val + + val = dzero + + end function d_null_func_3d + + ! + ! subroutine to allocate and fill in the coefficient matrix and + ! the rhs. + ! + subroutine amg_d_gen_pde3d(ctxt,idim,a,bv,xv,desc_a,afmt,& + & a1,a2,a3,b1,b2,b3,c,g,info,f,amold,vmold,partition, nrl,iv) + use psb_base_mod + use psb_util_mod + ! + ! Discretizes the partial differential equation + ! + ! d a1 d(u) d a1 d(u) d a1 d(u) b1 d(u) b2 d(u) b3 d(u) + ! - ------ - ------ - ------ + ----- + ------ + ------ + c u = f + ! dx dx dy dy dz dz dx dy dz + ! + ! with Dirichlet boundary conditions + ! u = g + ! + ! on the unit cube 0<=x,y,z<=1. + ! + ! + ! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation. + ! + implicit none + procedure(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_) :: info + type(psb_ctxt_type) :: ctxt + character :: afmt*5 + procedure(d_func_3d), optional :: f + class(psb_d_base_sparse_mat), optional :: amold + class(psb_d_base_vect_type), optional :: vmold + integer(psb_ipk_), optional :: partition, nrl,iv(:) + + ! Local variables. + + integer(psb_ipk_), parameter :: nb=20 + type(psb_d_csc_sparse_mat) :: acsc + type(psb_d_coo_sparse_mat) :: acoo + type(psb_d_csr_sparse_mat) :: acsr + integer(psb_ipk_) :: nnz,nr,nlr,i,j,ii,ib,k, partition_ + integer(psb_lpk_) :: m,n,glob_row,nt + integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner + ! For 3D partition + ! Note: integer control variables going directly into an MPI call + ! must be 4 bytes, i.e. psb_mpk_ + integer(psb_mpk_) :: npdims(3), npp, minfo + integer(psb_ipk_) :: npx,npy,npz, iamx,iamy,iamz,mynx,myny,mynz + integer(psb_ipk_), allocatable :: bndx(:),bndy(:),bndz(:) + ! Process grid + integer(psb_ipk_) :: np, iam + integer(psb_ipk_) :: icoeff + integer(psb_lpk_), allocatable :: myidx(:) + ! deltah dimension of each grid cell + ! deltat discretization time + real(psb_dpk_) :: deltah, sqdeltah, deltah2 + real(psb_dpk_), parameter :: rhs=dzero,one=done,zero=dzero + real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb + integer(psb_ipk_) :: err_act + procedure(d_func_3d), pointer :: f_ + character(len=20) :: name, ch_err,tmpfmt + + info = psb_success_ + name = 'd_create_matrix' + call psb_erractionsave(err_act) + + call psb_info(ctxt, iam, np) + + + if (present(f)) then + f_ => f + else + f_ => d_null_func_3d + end if + + if (present(partition)) then + if ((1<= partition).and.(partition <= 3)) then + partition_ = partition + else + write(*,*) 'Invalid partition choice ',partition,' defaulting to 3' + partition_ = 3 + end if + else + partition_ = 3 + end if + deltah = done/(idim+2) + sqdeltah = deltah*deltah + deltah2 = 2.0_psb_dpk_* deltah + + if (present(partition)) then + if ((1<= partition).and.(partition <= 3)) then + partition_ = partition + else + write(*,*) 'Invalid partition choice ',partition,' defaulting to 3' + partition_ = 3 + end if + else + partition_ = 3 + end if + + ! initialize array descriptor and sparse matrix storage. provide an + ! estimate of the number of non zeroes + + m = (1_psb_lpk_*idim)*idim*idim + n = m + nnz = 7*((n+np-1)/np) + if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n + t0 = psb_wtime() + select case(partition_) + case(1) + ! A BLOCK partition + if (present(nrl)) then + nr = nrl + else + ! + ! Using a simple BLOCK distribution. + ! + nt = (m+np-1)/np + nr = max(0,min(nt,m-(iam*nt))) + end if + + nt = nr + call psb_sum(ctxt,nt) + if (nt /= m) then + write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + + ! + ! First example of use of CDALL: specify for each process a number of + ! contiguous rows + ! + call psb_cdall(ctxt,desc_a,info,nl=nr) + myidx = desc_a%get_global_indices() + nlr = size(myidx) + + case(2) + ! A partition defined by the user through IV + + if (present(iv)) then + if (size(iv) /= m) then + write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + else + write(psb_err_unit,*) iam, 'Initialization error: IV not present' + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + + ! + ! Second example of use of CDALL: specify for each row the + ! process that owns it + ! + call psb_cdall(ctxt,desc_a,info,vg=iv) + myidx = desc_a%get_global_indices() + nlr = size(myidx) + + case(3) + ! A 3-dimensional partition + + ! A nifty MPI function will split the process list + npdims = 0 + call mpi_dims_create(np,3,npdims,info) + npx = npdims(1) + npy = npdims(2) + npz = npdims(3) + + allocate(bndx(0:npx),bndy(0:npy),bndz(0:npz)) + ! We can reuse idx2ijk for process indices as well. + call idx2ijk(iamx,iamy,iamz,iam,npx,npy,npz,base=0) + ! Now let's split the 3D cube in hexahedra + call dist1Didx(bndx,idim,npx) + mynx = bndx(iamx+1)-bndx(iamx) + call dist1Didx(bndy,idim,npy) + myny = bndy(iamy+1)-bndy(iamy) + call dist1Didx(bndz,idim,npz) + mynz = bndz(iamz+1)-bndz(iamz) + + ! How many indices do I own? + nlr = mynx*myny*mynz + allocate(myidx(nlr)) + ! Now, let's generate the list of indices I own + nr = 0 + do i=bndx(iamx),bndx(iamx+1)-1 + do j=bndy(iamy),bndy(iamy+1)-1 + do k=bndz(iamz),bndz(iamz+1)-1 + nr = nr + 1 + call ijk2idx(myidx(nr),i,j,k,idim,idim,idim) + end do + end do + end do + if (nr /= nlr) then + write(psb_err_unit,*) iam,iamx,iamy,iamz, 'Initialization error: NR vs NLR ',& + & nr,nlr,mynx,myny,mynz + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + end if + + ! + ! Third example of use of CDALL: specify for each process + ! the set of global indices it owns. + ! + call psb_cdall(ctxt,desc_a,info,vl=myidx) + + ! + ! Specify process topology + ! + block + ! + ! Use adjcncy methods + ! + integer(psb_mpk_), allocatable :: neighbours(:) + integer(psb_mpk_) :: cnt + logical, parameter :: debug_adj=.true. + if (debug_adj.and.(np > 1)) then + cnt = 0 + allocate(neighbours(np)) + if (iamx < npx-1) then + cnt = cnt + 1 + call ijk2idx(neighbours(cnt),iamx+1,iamy,iamz,npx,npy,npz,base=0) + end if + if (iamy < npy-1) then + cnt = cnt + 1 + call ijk2idx(neighbours(cnt),iamx,iamy+1,iamz,npx,npy,npz,base=0) + end if + if (iamz < npz-1) then + cnt = cnt + 1 + call ijk2idx(neighbours(cnt),iamx,iamy,iamz+1,npx,npy,npz,base=0) + end if + if (iamx >0) then + cnt = cnt + 1 + call ijk2idx(neighbours(cnt),iamx-1,iamy,iamz,npx,npy,npz,base=0) + end if + if (iamy >0) then + cnt = cnt + 1 + call ijk2idx(neighbours(cnt),iamx,iamy-1,iamz,npx,npy,npz,base=0) + end if + if (iamz >0) then + cnt = cnt + 1 + call ijk2idx(neighbours(cnt),iamx,iamy,iamz-1,npx,npy,npz,base=0) + end if + call psb_realloc(cnt, neighbours,info) + call desc_a%set_p_adjcncy(neighbours) + !write(0,*) iam,' Check on neighbours: ',desc_a%get_p_adjcncy() + end if + end block + + case default + write(psb_err_unit,*) iam, 'Initialization error: should not get here' + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end select + + + if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz) + ! define rhs from boundary conditions; also build initial guess + if (info == psb_success_) call psb_geall(xv,desc_a,info) + if (info == psb_success_) call psb_geall(bv,desc_a,info) + + call psb_barrier(ctxt) + talc = psb_wtime()-t0 + + call psb_barrier(ctxt) + t1 = psb_wtime() + !$omp parallel shared(deltah,myidx,a,desc_a) + ! + block + integer(psb_ipk_) :: i,j,k,ii,ib,icoeff, ix,iy,iz, ith,nth + integer(psb_lpk_) :: glob_row + integer(psb_lpk_), allocatable :: irow(:),icol(:) + real(psb_dpk_), allocatable :: val(:) + real(psb_dpk_) :: x,y,z, zt(nb) + nth = 1 + ith = 0 + allocate(val(20*nb),irow(20*nb),& + &icol(20*nb),stat=info) + if (info /= psb_success_ ) then + info=psb_err_alloc_dealloc_ + call psb_errpush(info,name) + !goto 9999 + endif + + !$omp do schedule(dynamic) + ! + do ii=1, nlr, nb + if (info /= psb_success_) cycle + ib = min(nb,nlr-ii+1) + icoeff = 1 + do k=1,ib + i=ii+k-1 + ! local matrix pointer + glob_row=myidx(i) + ! compute gridpoint coordinates + call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim) + ! x, y, z coordinates + x = (ix-1)*deltah + y = (iy-1)*deltah + z = (iz-1)*deltah + zt(k) = f_(x,y,z) + ! internal point: build discretization + ! + ! term depending on (x-1,y,z) + ! + val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2 + if (ix == 1) then + zt(k) = g(dzero,y,z)*(-val(icoeff)) + zt(k) + else + call ijk2idx(icol(icoeff),ix-1,iy,iz,idim,idim,idim) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y-1,z) + val(icoeff) = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2 + if (iy == 1) then + zt(k) = g(x,dzero,z)*(-val(icoeff)) + zt(k) + else + call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y,z-1) + val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2 + if (iz == 1) then + zt(k) = g(x,y,dzero)*(-val(icoeff)) + zt(k) + else + call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + ! term depending on (x,y,z) + val(icoeff)=(2*done)*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah & + & + c(x,y,z) + call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim) + irow(icoeff) = glob_row + icoeff = icoeff+1 + ! term depending on (x,y,z+1) + val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2 + if (iz == idim) then + zt(k) = g(x,y,done)*(-val(icoeff)) + zt(k) + else + call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y+1,z) + val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2 + if (iy == idim) then + zt(k) = g(x,done,z)*(-val(icoeff)) + zt(k) + else + call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x+1,y,z) + val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2 + if (ix==idim) then + zt(k) = g(done,y,z)*(-val(icoeff)) + zt(k) + else + call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + end do + !write(0,*) ' Outer in_parallel ',omp_in_parallel() + call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info) + if(info /= psb_success_) cycle + call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info) + if(info /= psb_success_) cycle + zt(:)=dzero + call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info) + if(info /= psb_success_) cycle + end do + !$omp end do + + deallocate(val,irow,icol) + end block + !$omp end parallel + + 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 + + + call psb_barrier(ctxt) + t1 = psb_wtime() + call psb_cdasb(desc_a,info) + tcdasb = psb_wtime()-t1 + call psb_barrier(ctxt) + t1 = psb_wtime() + if (info == psb_success_) then + if (present(amold)) then + call psb_spasb(a,desc_a,info,mold=amold) + else + call psb_spasb(a,desc_a,info,afmt=afmt) + end if + end if + call psb_barrier(ctxt) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='asb rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + if (info == psb_success_) call psb_geasb(xv,desc_a,info,mold=vmold) + if (info == psb_success_) call psb_geasb(bv,desc_a,info,mold=vmold) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='asb rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + tasb = psb_wtime()-t1 + call psb_barrier(ctxt) + ttot = psb_wtime() - t0 + + call psb_amx(ctxt,talc) + call psb_amx(ctxt,tgen) + call psb_amx(ctxt,tasb) + call psb_amx(ctxt,ttot) + if(iam == psb_root_) then + tmpfmt = a%get_fmt() + write(psb_out_unit,'("The matrix has been generated and assembled in ",a3," format.")')& + & tmpfmt + write(psb_out_unit,'("-allocation time : ",es12.5)') talc + write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen + write(psb_out_unit,'("-desc asbly time : ",es12.5)') tcdasb + write(psb_out_unit,'("- mat asbly time : ",es12.5)') tasb + write(psb_out_unit,'("-total time : ",es12.5)') ttot + + end if + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ctxt,err_act) + + return + end subroutine amg_d_gen_pde3d + + + + ! + ! subroutine to allocate and fill in the coefficient matrix and + ! the rhs. + ! + subroutine amg_d_gen_pde2d(ctxt,idim,a,bv,xv,desc_a,afmt,& + & a1,a2,b1,b2,c,g,info,f,amold,vmold,partition, nrl,iv) + use psb_base_mod + use psb_util_mod + ! + ! Discretizes the partial differential equation + ! + ! d d(u) d d(u) b1 d(u) b2 d(u) + ! - -- a1 ---- - -- a1 ---- + ----- + ------ + c u = f + ! dx dx dy dy dx dy + ! + ! with Dirichlet boundary conditions + ! u = g + ! + ! on the unit square 0<=x,y<=1. + ! + ! + ! Note that if b1=b2=c=0., the PDE is the Laplace equation. + ! + implicit none + procedure(d_func_2d) :: b1,b2,c,a1,a2,g + integer(psb_ipk_) :: idim + type(psb_dspmat_type) :: a + type(psb_d_vect_type) :: xv,bv + type(psb_desc_type) :: desc_a + integer(psb_ipk_) :: info + type(psb_ctxt_type) :: ctxt + character :: afmt*5 + procedure(d_func_2d), optional :: f + class(psb_d_base_sparse_mat), optional :: amold + class(psb_d_base_vect_type), optional :: vmold + integer(psb_ipk_), optional :: partition, nrl,iv(:) + ! Local variables. + + integer(psb_ipk_), parameter :: nb=20 + type(psb_d_csc_sparse_mat) :: acsc + type(psb_d_coo_sparse_mat) :: acoo + type(psb_d_csr_sparse_mat) :: acsr + integer(psb_ipk_) :: nnz,nr,nlr,i,j,ii,ib,k, partition_ + integer(psb_lpk_) :: m,n,glob_row,nt + integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner + ! For 2D partition + ! Note: integer control variables going directly into an MPI call + ! must be 4 bytes, i.e. psb_mpk_ + integer(psb_mpk_) :: npdims(2), npp, minfo + integer(psb_ipk_) :: npx,npy,iamx,iamy,mynx,myny + integer(psb_ipk_), allocatable :: bndx(:),bndy(:) + ! Process grid + integer(psb_ipk_) :: np, iam + integer(psb_ipk_) :: icoeff + integer(psb_lpk_), allocatable :: myidx(:) + ! deltah dimension of each grid cell + ! deltat discretization time + real(psb_dpk_) :: deltah, sqdeltah, deltah2, dd + real(psb_dpk_), parameter :: rhs=0.d0,one=done,zero=0.d0 + real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb + integer(psb_ipk_) :: err_act + procedure(d_func_2d), pointer :: f_ + character(len=20) :: name, ch_err,tmpfmt + + info = psb_success_ + name = 'create_matrix' + call psb_erractionsave(err_act) + + call psb_info(ctxt, iam, np) + + + if (present(f)) then + f_ => f + else + f_ => d_null_func_2d + end if + + deltah = done/(idim+2) + sqdeltah = deltah*deltah + deltah2 = 2.0_psb_dpk_* deltah + + + if (present(partition)) then + if ((1<= partition).and.(partition <= 3)) then + partition_ = partition + else + write(*,*) 'Invalid partition choice ',partition,' defaulting to 3' + partition_ = 3 + end if + else + partition_ = 3 + end if + + ! initialize array descriptor and sparse matrix storage. provide an + ! estimate of the number of non zeroes + + m = (1_psb_lpk_)*idim*idim + n = m + nnz = 7*((n+np-1)/np) + if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n + t0 = psb_wtime() + select case(partition_) + case(1) + ! A BLOCK partition + if (present(nrl)) then + nr = nrl + else + ! + ! Using a simple BLOCK distribution. + ! + nt = (m+np-1)/np + nr = max(0,min(nt,m-(iam*nt))) + end if + + nt = nr + call psb_sum(ctxt,nt) + if (nt /= m) then + write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + + ! + ! First example of use of CDALL: specify for each process a number of + ! contiguous rows + ! + call psb_cdall(ctxt,desc_a,info,nl=nr) + myidx = desc_a%get_global_indices() + nlr = size(myidx) + + case(2) + ! A partition defined by the user through IV + + if (present(iv)) then + if (size(iv) /= m) then + write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + else + write(psb_err_unit,*) iam, 'Initialization error: IV not present' + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + + ! + ! Second example of use of CDALL: specify for each row the + ! process that owns it + ! + call psb_cdall(ctxt,desc_a,info,vg=iv) + myidx = desc_a%get_global_indices() + nlr = size(myidx) + + case(3) + ! A 2-dimensional partition + + ! A nifty MPI function will split the process list + npdims = 0 + call mpi_dims_create(np,2,npdims,info) + npx = npdims(1) + npy = npdims(2) + + allocate(bndx(0:npx),bndy(0:npy)) + ! We can reuse idx2ijk for process indices as well. + call idx2ijk(iamx,iamy,iam,npx,npy,base=0) + ! Now let's split the 2D square in rectangles + call dist1Didx(bndx,idim,npx) + mynx = bndx(iamx+1)-bndx(iamx) + call dist1Didx(bndy,idim,npy) + myny = bndy(iamy+1)-bndy(iamy) + + ! How many indices do I own? + nlr = mynx*myny + allocate(myidx(nlr)) + ! Now, let's generate the list of indices I own + nr = 0 + do i=bndx(iamx),bndx(iamx+1)-1 + do j=bndy(iamy),bndy(iamy+1)-1 + nr = nr + 1 + call ijk2idx(myidx(nr),i,j,idim,idim) + end do + end do + if (nr /= nlr) then + write(psb_err_unit,*) iam,iamx,iamy, 'Initialization error: NR vs NLR ',& + & nr,nlr,mynx,myny + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + end if + + ! + ! Third example of use of CDALL: specify for each process + ! the set of global indices it owns. + ! + call psb_cdall(ctxt,desc_a,info,vl=myidx) + + ! + ! Specify process topology + ! + block + ! + ! Use adjcncy methods + ! + integer(psb_mpk_), allocatable :: neighbours(:) + integer(psb_mpk_) :: cnt + logical, parameter :: debug_adj=.true. + if (debug_adj.and.(np > 1)) then + cnt = 0 + allocate(neighbours(np)) + if (iamx < npx-1) then + cnt = cnt + 1 + call ijk2idx(neighbours(cnt),iamx+1,iamy,npx,npy,base=0) + end if + if (iamy < npy-1) then + cnt = cnt + 1 + call ijk2idx(neighbours(cnt),iamx,iamy+1,npx,npy,base=0) + end if + if (iamx >0) then + cnt = cnt + 1 + call ijk2idx(neighbours(cnt),iamx-1,iamy,npx,npy,base=0) + end if + if (iamy >0) then + cnt = cnt + 1 + call ijk2idx(neighbours(cnt),iamx,iamy-1,npx,npy,base=0) + end if + call psb_realloc(cnt, neighbours,info) + call desc_a%set_p_adjcncy(neighbours) + !write(0,*) iam,' Check on neighbours: ',desc_a%get_p_adjcncy() + end if + end block + + case default + write(psb_err_unit,*) iam, 'Initialization error: should not get here' + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end select + + + if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz) + ! define rhs from boundary conditions; also build initial guess + if (info == psb_success_) call psb_geall(xv,desc_a,info) + if (info == psb_success_) call psb_geall(bv,desc_a,info) + + call psb_barrier(ctxt) + talc = psb_wtime()-t0 + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='allocation rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + call psb_barrier(ctxt) + t1 = psb_wtime() + !$omp parallel shared(deltah,myidx,a,desc_a) + ! + block + integer(psb_ipk_) :: i,j,k,ii,ib,icoeff, ix,iy,iz, ith,nth + integer(psb_lpk_) :: glob_row + integer(psb_lpk_), allocatable :: irow(:),icol(:) + real(psb_dpk_), allocatable :: val(:) + real(psb_dpk_) :: x,y,z, zt(nb) + nth = 1 + ith = 0 + allocate(val(20*nb),irow(20*nb),& + &icol(20*nb),stat=info) + if (info /= psb_success_ ) then + info=psb_err_alloc_dealloc_ + call psb_errpush(info,name) + !goto 9999 + endif + + ! loop over rows belonging to current process in a block + ! distribution. + !$omp do schedule(dynamic) + ! + do ii=1, nlr,nb + ib = min(nb,nlr-ii+1) + icoeff = 1 + do k=1,ib + i=ii+k-1 + ! local matrix pointer + glob_row=myidx(i) + ! compute gridpoint coordinates + call idx2ijk(ix,iy,glob_row,idim,idim) + ! x, y coordinates + x = (ix-1)*deltah + y = (iy-1)*deltah + + zt(k) = f_(x,y) + ! internal point: build discretization + ! + ! term depending on (x-1,y) + ! + val(icoeff) = -a1(x,y)/sqdeltah-b1(x,y)/deltah2 + if (ix == 1) then + zt(k) = g(dzero,y)*(-val(icoeff)) + zt(k) + else + call ijk2idx(icol(icoeff),ix-1,iy,idim,idim) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x,y-1) + val(icoeff) = -a2(x,y)/sqdeltah-b2(x,y)/deltah2 + if (iy == 1) then + zt(k) = g(x,dzero)*(-val(icoeff)) + zt(k) + else + call ijk2idx(icol(icoeff),ix,iy-1,idim,idim) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + ! term depending on (x,y) + val(icoeff)=(2*done)*(a1(x,y) + a2(x,y))/sqdeltah + c(x,y) + call ijk2idx(icol(icoeff),ix,iy,idim,idim) + irow(icoeff) = glob_row + icoeff = icoeff+1 + ! term depending on (x,y+1) + val(icoeff)=-a2(x,y)/sqdeltah+b2(x,y)/deltah2 + if (iy == idim) then + zt(k) = g(x,done)*(-val(icoeff)) + zt(k) + else + call ijk2idx(icol(icoeff),ix,iy+1,idim,idim) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + ! term depending on (x+1,y) + val(icoeff)=-a1(x,y)/sqdeltah+b1(x,y)/deltah2 + if (ix==idim) then + zt(k) = g(done,y)*(-val(icoeff)) + zt(k) + else + call ijk2idx(icol(icoeff),ix+1,iy,idim,idim) + irow(icoeff) = glob_row + icoeff = icoeff+1 + endif + + end do + call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info) + if(info /= psb_success_) cycle + call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info) + if(info /= psb_success_) cycle + zt(:)=dzero + call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info) + if(info /= psb_success_) cycle + end do + !$omp end do + + deallocate(val,irow,icol) + end block + !$omp end parallel + + 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 + + call psb_barrier(ctxt) + t1 = psb_wtime() + call psb_cdasb(desc_a,info) + tcdasb = psb_wtime()-t1 + call psb_barrier(ctxt) + t1 = psb_wtime() + if (info == psb_success_) then + if (present(amold)) then + call psb_spasb(a,desc_a,info,mold=amold) + else + call psb_spasb(a,desc_a,info,afmt=afmt) + end if + end if + call psb_barrier(ctxt) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='asb rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + if (info == psb_success_) call psb_geasb(xv,desc_a,info,mold=vmold) + if (info == psb_success_) call psb_geasb(bv,desc_a,info,mold=vmold) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='asb rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + tasb = psb_wtime()-t1 + call psb_barrier(ctxt) + ttot = psb_wtime() - t0 + + call psb_amx(ctxt,talc) + call psb_amx(ctxt,tgen) + call psb_amx(ctxt,tasb) + call psb_amx(ctxt,ttot) + if(iam == psb_root_) then + tmpfmt = a%get_fmt() + write(psb_out_unit,'("The matrix has been generated and assembled in ",a3," format.")')& + & tmpfmt + write(psb_out_unit,'("-allocation time : ",es12.5)') talc + write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen + write(psb_out_unit,'("-desc asbly time : ",es12.5)') tcdasb + write(psb_out_unit,'("- mat asbly time : ",es12.5)') tasb + write(psb_out_unit,'("-total time : ",es12.5)') ttot + + end if + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ctxt) + return + end if + return + end subroutine amg_d_gen_pde2d +end module amg_d_genpde_mod + +program amg_d_pde3d + use psb_base_mod + use amg_prec_mod + use psb_linsolve_mod + use psb_util_mod + use data_input + use amg_d_genpde_mod + use amg_d_pde3d_poisson_mod + use amg_d_tlu_solver + implicit none + + ! input parameters + character(len=20) :: kmethd, ptype + character(len=5) :: afmt + integer(psb_ipk_) :: idim + + ! miscellaneous + real(psb_dpk_) :: t1, t2, tprec, thier, tslv + + ! sparse matrix and preconditioner + type(psb_dspmat_type) :: a + type(amg_dprec_type) :: prec + type(amg_d_tlu_solver_type) :: tlusv + ! descriptor + type(psb_desc_type) :: desc_a + ! dense vectors + type(psb_d_vect_type) :: x,b + ! parallel environment + type(psb_ctxt_type) :: ctxt + integer(psb_ipk_) :: iam, np + + ! solver parameters + integer(psb_ipk_) :: iter, itmax,itrace, istopc, irst, nlv + integer(psb_epk_) :: amatsize, precsize, descsize + real(psb_dpk_) :: err, eps + + ! other variables + integer(psb_ipk_) :: info, i + character(len=20) :: name,ch_err + + info=psb_success_ + + + call psb_init(ctxt) + call psb_info(ctxt,iam,np) + + if (iam < 0) then + ! This should not happen, but just in case + call psb_exit(ctxt) + stop + endif + if(psb_get_errstatus() /= 0) goto 9999 + name='amg_d_pde3d' + call psb_set_errverbosity(itwo) + ! + ! Hello world + ! + if (iam == psb_root_) then + write(*,*) 'Welcome to AMG4PSBLAS version: ',amg_version_string_ + write(*,*) 'This is the ',trim(name),' sample program' + end if + + ! + ! get parameters + ! + call get_parms(ctxt,kmethd,afmt,idim,istopc,itmax,itrace,irst,eps) + + ! + ! allocate and fill in the coefficient matrix, rhs and initial guess + ! + + call psb_barrier(ctxt) + t1 = psb_wtime() + call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,& + & a1_poisson,a2_poisson,a3_poisson,& + & b1_poisson,b2_poisson,b3_poisson,c_poisson,g_poisson,info) + call psb_barrier(ctxt) + 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,'(" ")') + ! + ! prepare the preconditioner: an ML with defaults, but with TLU solver at + ! intermediate levels. All other parameters are at default values. + ! + call prec%init(ctxt,'ML', info) + + call psb_barrier(ctxt) + t1 = psb_wtime() + call prec%hierarchy_build(a,desc_a,info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='prec%hierarchy_bld' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + thier = psb_wtime()-t1 + nlv = prec%get_nlevs() + call prec%set(tlusv, info,ilev=1,ilmax=max(1,nlv-1)) + + call psb_barrier(ctxt) + t1 = psb_wtime() + call prec%smoothers_build(a,desc_a,info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='prec%smoothers_build' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + tprec = psb_wtime()-t1 + + call psb_amx(ctxt,thier) + call psb_amx(ctxt,tprec) + + if (iam == psb_root_) & + & write(psb_out_unit,'("Preconditioner time : ",es12.5)') tprec+thier + if (iam == psb_root_) call prec%descr(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(ctxt) + 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(ctxt) + tslv = psb_wtime() - t1 + call psb_amx(ctxt,tslv) + + amatsize = a%sizeof() + descsize = desc_a%sizeof() + precsize = prec%sizeof() + call psb_sum(ctxt,amatsize) + call psb_sum(ctxt,descsize) + call psb_sum(ctxt,precsize) + if (iam == psb_root_) then + write(psb_out_unit,'(" ")') + write(psb_out_unit,'("Numer of levels of aggr. hierarchy: ",i12)') prec%get_nlevs() + write(psb_out_unit,'("Time to build aggr. hierarchy : ",es12.5)') thier + write(psb_out_unit,'("Time to build smoothers : ",es12.5)') tprec + write(psb_out_unit,'("Total preconditioner time : ",es12.5)') tprec+thier + write(psb_out_unit,'("Time to solve system : ",es12.5)') tslv + write(psb_out_unit,'("Time per iteration : ",es12.5)') tslv/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,'("Storage format for A: ",a)') trim(a%get_fmt()) + write(psb_out_unit,'("Total memory occupation for DESC_A: ",i12)') descsize + write(psb_out_unit,'("Storage format for DESC_A: ",a)') trim(desc_a%get_fmt()) + 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 prec%free(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 + call psb_exit(ctxt) + stop + +9999 continue + call psb_error(ctxt) + +contains + ! + ! get iteration parameters from standard input + ! + subroutine get_parms(ctxt,kmethd,afmt,idim,istopc,itmax,itrace,irst,eps) + + type(psb_ctxt_type) :: ctxt + character(len=*) :: kmethd, afmt + integer(psb_ipk_) :: idim, istopc,itmax,itrace,irst + integer(psb_ipk_) :: np, iam, info, inp_unit + real(psb_dpk_) :: eps + character(len=20) :: buffer + character(len=1024) :: filename + + call psb_info(ctxt, iam, np) + + + if (iam == psb_root_) then + if (command_argument_count()>0) then + call get_command_argument(1,filename) + inp_unit = 30 + open(inp_unit,file=filename,action='read',iostat=info) + if (info /= 0) then + write(psb_err_unit,*) 'Could not open file ',filename,' for input' + call psb_abort(ctxt) + stop + else + write(psb_err_unit,*) 'Opened file ',trim(filename),' for input' + end if + else + inp_unit=psb_inp_unit + end if + call read_data(kmethd,inp_unit) + call read_data(afmt,inp_unit) + call read_data(idim,inp_unit) + call read_data(istopc,inp_unit) + call read_data(itmax,inp_unit) + call read_data(itrace,inp_unit) + call read_data(irst,inp_unit) + call read_data(eps,inp_unit) + + if (inp_unit /= psb_inp_unit) then + close(inp_unit) + end if + end if + + ! broadcast parameters to all processors + call psb_bcast(ctxt,kmethd) + call psb_bcast(ctxt,afmt) + call psb_bcast(ctxt,idim) + call psb_bcast(ctxt,istopc) + call psb_bcast(ctxt,itmax) + call psb_bcast(ctxt,itrace) + call psb_bcast(ctxt,irst) + call psb_bcast(ctxt,eps) + + 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)') 'ML-TLU' + 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(psb_ipk_) :: iout + write(iout,*)'incorrect parameter(s) found' + write(iout,*)' usage: amg_d_pde3d 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 + +end program amg_d_pde3d diff --git a/tests/newslv/data_input.f90 b/samples/newslv/data_input.f90 similarity index 62% rename from tests/newslv/data_input.f90 rename to samples/newslv/data_input.f90 index e6a8582b..6b961352 100644 --- a/tests/newslv/data_input.f90 +++ b/samples/newslv/data_input.f90 @@ -1,14 +1,14 @@ ! ! -! MLD2P4 version 2.2 -! MultiLevel Domain Decomposition Parallel Preconditioners Package -! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) ! -! (C) Copyright 2008-2018 +! (C) Copyright 2021 ! ! Salvatore Filippone ! Pasqua D'Ambra -! Daniela di Serafino +! Fabio Durastante ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions @@ -18,14 +18,14 @@ ! 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 +! 3. The name of the AMG4PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. ! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS 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 @@ -39,10 +39,12 @@ module data_input use psb_base_mod, only : psb_spk_, psb_dpk_, psb_ipk_ interface read_data - module procedure read_char, read_int,& - & read_double, read_single, read_logical,& - & string_read_char, string_read_int,& - & string_read_double, string_read_single, & + module procedure read_char, read_int, read_int_array,& + & read_double, read_double_array, & + & read_single, read_single_array, read_logical,& + & string_read_char, string_read_int, string_read_int_array,& + & string_read_double, string_read_double_array,& + & string_read_single, string_read_single_array, & & string_read_logical end interface read_data interface trim_string @@ -51,15 +53,33 @@ module data_input character(len=4096), private :: charbuf character, private, parameter :: def_marker="!" + character, private, parameter :: cmt_marker="%" contains - + subroutine get_buffer(file,buffer) + integer(psb_ipk_), intent(in) :: file + character(len=*), intent(inout) :: buffer + integer :: idx + do + read(file,'(a)',end=999) buffer + buffer = adjustl(buffer) + idx=index(charbuf,cmt_marker) + if (idx == 1 ) then + cycle + else + exit + end if + end do +999 continue + return + end subroutine get_buffer + subroutine read_logical(val,file,marker) logical, intent(out) :: val integer(psb_ipk_), intent(in) :: file character(len=1), optional, intent(in) :: marker - read(file,'(a)')charbuf + call get_buffer(file,charbuf) call read_data(val,charbuf,marker) end subroutine read_logical @@ -69,7 +89,7 @@ contains integer(psb_ipk_), intent(in) :: file character(len=1), optional, intent(in) :: marker - read(file,'(a)')charbuf + call get_buffer(file,charbuf) call read_data(val,charbuf,marker) end subroutine read_char @@ -79,29 +99,61 @@ contains integer(psb_ipk_), intent(in) :: file character(len=1), optional, intent(in) :: marker - read(file,'(a)')charbuf + call get_buffer(file,charbuf) call read_data(val,charbuf,marker) end subroutine read_int + + subroutine read_int_array(val,file,marker) + integer(psb_ipk_), intent(out) :: val(:) + integer(psb_ipk_), intent(in) :: file + character(len=1), optional, intent(in) :: marker + + call get_buffer(file,charbuf) + call read_data(val,charbuf,marker) + + end subroutine read_int_array + subroutine read_single(val,file,marker) real(psb_spk_), intent(out) :: val integer(psb_ipk_), intent(in) :: file character(len=1), optional, intent(in) :: marker - read(file,'(a)')charbuf + call get_buffer(file,charbuf) call read_data(val,charbuf,marker) end subroutine read_single + + subroutine read_single_array(val,file,marker) + real(psb_spk_), intent(out) :: val(:) + integer(psb_ipk_), intent(in) :: file + character(len=1), optional, intent(in) :: marker + + call get_buffer(file,charbuf) + call read_data(val,charbuf,marker) + + end subroutine read_single_array + subroutine read_double(val,file,marker) real(psb_dpk_), intent(out) :: val integer(psb_ipk_), intent(in) :: file character(len=1), optional, intent(in) :: marker - read(file,'(a)')charbuf + call get_buffer(file,charbuf) call read_data(val,charbuf,marker) end subroutine read_double + subroutine read_double_array(val,file,marker) + real(psb_dpk_), intent(out) :: val(:) + integer(psb_ipk_), intent(in) :: file + character(len=1), optional, intent(in) :: marker + + call get_buffer(file,charbuf) + call read_data(val,charbuf,marker) + + end subroutine read_double_array + subroutine string_read_char(val,file,marker) character(len=*), intent(out) :: val character(len=*), intent(in) :: file @@ -140,6 +192,25 @@ contains read(charbuf(1:idx-1),*) val end subroutine string_read_int + subroutine string_read_int_array(val,file,marker) + integer(psb_ipk_), intent(out) :: val(:) + character(len=*), intent(in) :: file + character(len=1), optional, intent(in) :: marker + character(len=1) :: marker_ + character(len=1024) :: charbuf + integer(psb_ipk_) :: idx + if (present(marker)) then + marker_ = marker + else + marker_ = def_marker + end if + read(file,'(a)')charbuf + charbuf = adjustl(charbuf) + idx=index(charbuf,marker_) + if (idx == 0) idx = len(charbuf)+1 + read(charbuf(1:idx-1),*) val(:) + end subroutine string_read_int_array + subroutine string_read_single(val,file,marker) real(psb_spk_), intent(out) :: val character(len=*), intent(in) :: file @@ -159,6 +230,25 @@ contains read(charbuf(1:idx-1),*) val end subroutine string_read_single + subroutine string_read_single_array(val,file,marker) + real(psb_spk_), intent(out) :: val(:) + character(len=*), intent(in) :: file + character(len=1), optional, intent(in) :: marker + character(len=1) :: marker_ + character(len=1024) :: charbuf + integer(psb_ipk_) :: idx + if (present(marker)) then + marker_ = marker + else + marker_ = def_marker + end if + read(file,'(a)')charbuf + charbuf = adjustl(charbuf) + idx=index(charbuf,marker_) + if (idx == 0) idx = len(charbuf)+1 + read(charbuf(1:idx-1),*) val(:) + end subroutine string_read_single_array + subroutine string_read_double(val,file,marker) real(psb_dpk_), intent(out) :: val character(len=*), intent(in) :: file @@ -178,6 +268,25 @@ contains read(charbuf(1:idx-1),*) val end subroutine string_read_double + subroutine string_read_double_array(val,file,marker) + real(psb_dpk_), intent(out) :: val(:) + character(len=*), intent(in) :: file + character(len=1), optional, intent(in) :: marker + character(len=1) :: marker_ + character(len=1024) :: charbuf + integer(psb_ipk_) :: idx + if (present(marker)) then + marker_ = marker + else + marker_ = def_marker + end if + read(file,'(a)')charbuf + charbuf = adjustl(charbuf) + idx=index(charbuf,marker_) + if (idx == 0) idx = len(charbuf)+1 + read(charbuf(1:idx-1),*) val(:) + end subroutine string_read_double_array + subroutine string_read_logical(val,file,marker) logical, intent(out) :: val character(len=*), intent(in) :: file diff --git a/tests/newslv/runs/mld_pde.inp b/samples/newslv/runs/amg_pde.inp similarity index 100% rename from tests/newslv/runs/mld_pde.inp rename to samples/newslv/runs/amg_pde.inp diff --git a/tests/newslv/Makefile b/tests/newslv/Makefile deleted file mode 100644 index 76fd30c0..00000000 --- a/tests/newslv/Makefile +++ /dev/null @@ -1,34 +0,0 @@ -MLDDIR=../.. -MLDINCDIR=$(MLDDIR)/include -include $(MLDINCDIR)/Make.inc.amg4psblas -MLDMODDIR=$(MLDDIR)/modules -MLDLIBDIR=$(MLDDIR)/lib -MLD_LIBS=-L$(MLDLIBDIR) -lpsb_linsolve -lmld_prec -lpsb_prec -FINCLUDES=$(FMFLAG). $(FMFLAG)$(MLDMODDIR) $(FMFLAG)$(MLDINCDIR) $(PSBLAS_INCLUDES) $(FIFLAG). - -PD3DOBJS=mld_pde3d_newslv.o data_input.o mld_d_tlu_solver.o mld_d_tlu_solver_impl.o -PSOBJS=spde.o data_input.o -EXEDIR=./runs -LINKOPT= - -all: mld_pde3d_newslv - - -mld_pde3d_newslv: $(PD3DOBJS) - $(FLINK) $(LINKOPT) $(PD3DOBJS) -o mld_pde3d_newslv $(MLD_LIBS) $(PSBLAS_LIBS) $(LDLIBS) - /bin/mv mld_pde3d_newslv $(EXEDIR) - - -mld_pde3d_newslv.o mld_d_tlu_solver_impl.o: data_input.o mld_d_tlu_solver.o - - -clean: - /bin/rm -f $(PD3DOBJS) $(EXEDIR)/mld_pde3d_newslv - -verycleanlib: - (cd ../..; make veryclean) -lib: - (cd ../../; make library) - - - diff --git a/tests/newslv/amg_pde3d_newslv.f90 b/tests/newslv/amg_pde3d_newslv.f90 deleted file mode 100644 index 7367408d..00000000 --- a/tests/newslv/amg_pde3d_newslv.f90 +++ /dev/null @@ -1,385 +0,0 @@ -! -! -! MLD2P4 version 2.2 -! MultiLevel Domain Decomposition Parallel Preconditioners Package -! based on PSBLAS (Parallel Sparse BLAS version 3.5) -! -! (C) Copyright 2008-2018 -! -! Salvatore Filippone -! Pasqua D'Ambra -! Daniela di Serafino -! -! 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: mld_d_pde3d.f90 -! -! Program: mld_d_pde3d -! 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. -! -module mld_d_pde3d_mod -contains - ! - ! functions parametrizing the differential equation - ! - function b1(x,y,z) - use psb_base_mod, only : psb_dpk_,done,dzero - real(psb_dpk_) :: b1 - real(psb_dpk_), intent(in) :: x,y,z - b1=dzero/sqrt((3*done)) - end function b1 - function b2(x,y,z) - use psb_base_mod, only : psb_dpk_,done,dzero - real(psb_dpk_) :: b2 - real(psb_dpk_), intent(in) :: x,y,z - b2=dzero/sqrt((3*done)) - end function b2 - function b3(x,y,z) - use psb_base_mod, only : psb_dpk_,done,dzero - real(psb_dpk_) :: b3 - real(psb_dpk_), intent(in) :: x,y,z - b3=dzero/sqrt((3*done)) - end function b3 - function c(x,y,z) - use psb_base_mod, only : psb_dpk_,done,dzero - real(psb_dpk_) :: c - real(psb_dpk_), intent(in) :: x,y,z - c=dzero - end function c - function a1(x,y,z) - use psb_base_mod, only : psb_dpk_,done,dzero - real(psb_dpk_) :: a1 - real(psb_dpk_), intent(in) :: x,y,z - a1=done!/80 - end function a1 - function a2(x,y,z) - use psb_base_mod, only : psb_dpk_,done,dzero - real(psb_dpk_) :: a2 - real(psb_dpk_), intent(in) :: x,y,z - a2=done!/80 - end function a2 - function a3(x,y,z) - use psb_base_mod, only : psb_dpk_,done,dzero - real(psb_dpk_) :: a3 - real(psb_dpk_), intent(in) :: x,y,z - a3=done!/80 - end function a3 - function g(x,y,z) - use psb_base_mod, only : psb_dpk_,done,dzero - real(psb_dpk_) :: g - real(psb_dpk_), intent(in) :: x,y,z - g = dzero - if (x == done) then - g = done - else if (x == dzero) then - g = exp(y**2-z**2) - end if - end function g -end module mld_d_pde3d_mod - -program mld_d_pde3d - use psb_base_mod - use mld_prec_mod - use psb_krylov_mod - use psb_util_mod - use data_input - use mld_d_pde3d_mod - use mld_d_tlu_solver - implicit none - - ! input parameters - character(len=20) :: kmethd, ptype - character(len=5) :: afmt - integer(psb_ipk_) :: idim - - ! miscellaneous - real(psb_dpk_) :: t1, t2, tprec, thier, tslv - - ! sparse matrix and preconditioner - type(psb_dspmat_type) :: a - type(mld_dprec_type) :: prec - type(mld_d_tlu_solver_type) :: tlusv - ! descriptor - type(psb_desc_type) :: desc_a - ! dense vectors - type(psb_d_vect_type) :: x,b - ! parallel environment - integer(psb_ipk_) :: ctxt, iam, np - - ! solver parameters - integer(psb_ipk_) :: iter, itmax,itrace, istopc, irst, nlv - integer(psb_epk_) :: amatsize, precsize, descsize - real(psb_dpk_) :: err, eps - - ! other variables - integer(psb_ipk_) :: info, i - character(len=20) :: name,ch_err - - info=psb_success_ - - - call psb_init(ctxt) - call psb_info(ctxt,iam,np) - - if (iam < 0) then - ! This should not happen, but just in case - call psb_exit(ctxt) - stop - endif - if(psb_get_errstatus() /= 0) goto 9999 - name='mld_d_pde3d' - call psb_set_errverbosity(itwo) - ! - ! 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(ctxt,kmethd,afmt,idim,istopc,itmax,itrace,irst,eps) - - ! - ! allocate and fill in the coefficient matrix, rhs and initial guess - ! - - call psb_barrier(ctxt) - t1 = psb_wtime() - call psb_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,& - & a1,a2,a3,b1,b2,b3,c,g,info) - call psb_barrier(ctxt) - 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,'(" ")') - ! - ! prepare the preconditioner: an ML with defaults, but with TLU solver at - ! intermediate levels. All other parameters are at default values. - ! - call prec%init('ML', info) - - call psb_barrier(ctxt) - t1 = psb_wtime() - call prec%hierarchy_build(a,desc_a,info) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='prec%hierarchy_bld' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - thier = psb_wtime()-t1 - nlv = prec%get_nlevs() - call prec%set(tlusv, info,ilev=1,ilmax=max(1,nlv-1)) - - call psb_barrier(ctxt) - t1 = psb_wtime() - call prec%smoothers_build(a,desc_a,info) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='prec%smoothers_build' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - tprec = psb_wtime()-t1 - - call psb_amx(ctxt,thier) - call psb_amx(ctxt,tprec) - - if (iam == psb_root_) & - & write(psb_out_unit,'("Preconditioner time : ",es12.5)') tprec+thier - if (iam == psb_root_) call prec%descr(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(ctxt) - 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(ctxt) - tslv = psb_wtime() - t1 - call psb_amx(ctxt,tslv) - - amatsize = a%sizeof() - descsize = desc_a%sizeof() - precsize = prec%sizeof() - call psb_sum(ctxt,amatsize) - call psb_sum(ctxt,descsize) - call psb_sum(ctxt,precsize) - if (iam == psb_root_) then - write(psb_out_unit,'(" ")') - write(psb_out_unit,'("Numer of levels of aggr. hierarchy: ",i12)') prec%get_nlevs() - write(psb_out_unit,'("Time to build aggr. hierarchy : ",es12.5)') thier - write(psb_out_unit,'("Time to build smoothers : ",es12.5)') tprec - write(psb_out_unit,'("Total preconditioner time : ",es12.5)') tprec+thier - write(psb_out_unit,'("Time to solve system : ",es12.5)') tslv - write(psb_out_unit,'("Time per iteration : ",es12.5)') tslv/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,'("Storage format for A: ",a)') trim(a%get_fmt()) - write(psb_out_unit,'("Total memory occupation for DESC_A: ",i12)') descsize - write(psb_out_unit,'("Storage format for DESC_A: ",a)') trim(desc_a%get_fmt()) - 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 prec%free(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 - call psb_exit(ctxt) - stop - -9999 continue - call psb_error(ctxt) - -contains - ! - ! get iteration parameters from standard input - ! - subroutine get_parms(ctxt,kmethd,afmt,idim,istopc,itmax,itrace,irst,eps) - - integer(psb_ipk_) :: ctxt - character(len=*) :: kmethd, afmt - integer(psb_ipk_) :: idim, istopc,itmax,itrace,irst - integer(psb_ipk_) :: np, iam, info - real(psb_dpk_) :: eps - character(len=20) :: buffer - - call psb_info(ctxt, iam, np) - - if (iam == psb_root_) then - call read_data(kmethd,psb_inp_unit) - call read_data(afmt,psb_inp_unit) - call read_data(idim,psb_inp_unit) - call read_data(istopc,psb_inp_unit) - call read_data(itmax,psb_inp_unit) - call read_data(itrace,psb_inp_unit) - call read_data(irst,psb_inp_unit) - call read_data(eps,psb_inp_unit) - end if - - ! broadcast parameters to all processors - call psb_bcast(ctxt,kmethd) - call psb_bcast(ctxt,afmt) - call psb_bcast(ctxt,idim) - call psb_bcast(ctxt,istopc) - call psb_bcast(ctxt,itmax) - call psb_bcast(ctxt,itrace) - call psb_bcast(ctxt,irst) - call psb_bcast(ctxt,eps) - - 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)') 'ML-TLU' - 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(psb_ipk_) :: iout - write(iout,*)'incorrect parameter(s) found' - write(iout,*)' usage: mld_d_pde3d 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 - -end program mld_d_pde3d