mirror of
https://github.com/sfilippone/amg4psblas.git
synced 2026-10-06 22:55:12 +00:00
mld2p4-2:
config/pac.m4 configure mlprec/Makefile mlprec/mld_base_prec_type.f90 mlprec/mld_c_as_smoother.f90 mlprec/mld_c_id_solver.f90 mlprec/mld_c_ilu_solver.f90 mlprec/mld_c_inner_mod.f90 mlprec/mld_c_jac_smoother.f90 mlprec/mld_c_move_alloc_mod.f90 mlprec/mld_c_prec_mod.f90 mlprec/mld_c_prec_type.f90 mlprec/mld_c_slu_solver.f90 mlprec/mld_caggrmap_bld.f90 mlprec/mld_caggrmat_asb.f90 mlprec/mld_caggrmat_nosmth_asb.F90 mlprec/mld_caggrmat_smth_asb.F90 mlprec/mld_ccoarse_bld.f90 mlprec/mld_cilu0_fact.f90 mlprec/mld_ciluk_fact.f90 mlprec/mld_cilut_fact.f90 mlprec/mld_cmlprec_aply.f90 mlprec/mld_cmlprec_bld.f90 mlprec/mld_cprecaply.f90 mlprec/mld_cprecbld.f90 mlprec/mld_cprecinit.F90 mlprec/mld_cprecset.F90 mlprec/mld_cslu_bld.f90 mlprec/mld_cslu_interface.c mlprec/mld_cslud_bld.f90 mlprec/mld_csp_renum.f90 mlprec/mld_cumf_bld.f90 mlprec/mld_d_as_smoother.f90 mlprec/mld_d_id_solver.f90 mlprec/mld_d_ilu_solver.f90 mlprec/mld_d_inner_mod.f90 mlprec/mld_d_move_alloc_mod.f90 mlprec/mld_d_prec_mod.f90 mlprec/mld_d_prec_type.f90 mlprec/mld_d_slu_solver.f90 mlprec/mld_d_sludist_solver.f90 mlprec/mld_d_umf_solver.f90 mlprec/mld_daggrmap_bld.f90 mlprec/mld_daggrmat_asb.f90 mlprec/mld_daggrmat_minnrg_asb.F90 mlprec/mld_daggrmat_nosmth_asb.F90 mlprec/mld_daggrmat_smth_asb.F90 mlprec/mld_dcoarse_bld.f90 mlprec/mld_dilu0_fact.f90 mlprec/mld_diluk_fact.f90 mlprec/mld_dilut_fact.f90 mlprec/mld_dmlprec_aply.f90 mlprec/mld_dmlprec_bld.f90 mlprec/mld_dprecaply.f90 mlprec/mld_dprecbld.f90 mlprec/mld_dprecinit.F90 mlprec/mld_dprecset.F90 mlprec/mld_dslu_bld.f90 mlprec/mld_dslu_interface.c mlprec/mld_dslud_bld.f90 mlprec/mld_dslud_interface.c mlprec/mld_dsp_renum.f90 mlprec/mld_dumf_interface.c mlprec/mld_inner_mod.f90 mlprec/mld_move_alloc_mod.f90 mlprec/mld_prec_mod.f90 mlprec/mld_s_as_smoother.f90 mlprec/mld_s_id_solver.f90 mlprec/mld_s_ilu_solver.f90 mlprec/mld_s_inner_mod.f90 mlprec/mld_s_jac_smoother.f90 mlprec/mld_s_move_alloc_mod.f90 mlprec/mld_s_prec_mod.f90 mlprec/mld_s_prec_type.f90 mlprec/mld_s_slu_solver.f90 mlprec/mld_saggrmap_bld.f90 mlprec/mld_saggrmat_asb.f90 mlprec/mld_saggrmat_nosmth_asb.F90 mlprec/mld_saggrmat_smth_asb.F90 mlprec/mld_scoarse_bld.f90 mlprec/mld_silu0_fact.f90 mlprec/mld_siluk_fact.f90 mlprec/mld_silut_fact.f90 mlprec/mld_smlprec_aply.f90 mlprec/mld_smlprec_bld.f90 mlprec/mld_sprecaply.f90 mlprec/mld_sprecbld.f90 mlprec/mld_sprecinit.F90 mlprec/mld_sprecset.F90 mlprec/mld_sslu_bld.f90 mlprec/mld_sslu_interface.c mlprec/mld_sslud_bld.f90 mlprec/mld_ssp_renum.f90 mlprec/mld_sumf_bld.f90 mlprec/mld_z_as_smoother.f90 mlprec/mld_z_id_solver.f90 mlprec/mld_z_ilu_solver.f90 mlprec/mld_z_inner_mod.f90 mlprec/mld_z_jac_smoother.f90 mlprec/mld_z_move_alloc_mod.f90 mlprec/mld_z_prec_mod.f90 mlprec/mld_z_prec_type.f90 mlprec/mld_z_slu_solver.f90 mlprec/mld_z_umf_solver.f90 mlprec/mld_zaggrmap_bld.f90 mlprec/mld_zaggrmat_asb.f90 mlprec/mld_zaggrmat_nosmth_asb.F90 mlprec/mld_zaggrmat_smth_asb.F90 mlprec/mld_zas_aply.f90 mlprec/mld_zas_bld.f90 mlprec/mld_zbaseprec_aply.f90 mlprec/mld_zbaseprec_bld.f90 mlprec/mld_zcoarse_bld.f90 mlprec/mld_zdiag_bld.f90 mlprec/mld_zfact_bld.f90 mlprec/mld_zilu0_fact.f90 mlprec/mld_zilu_bld.f90 mlprec/mld_ziluk_fact.f90 mlprec/mld_zilut_fact.f90 mlprec/mld_zmlprec_aply.f90 mlprec/mld_zmlprec_bld.f90 mlprec/mld_zprecaply.f90 mlprec/mld_zprecbld.f90 mlprec/mld_zprecinit.F90 mlprec/mld_zprecset.F90 mlprec/mld_zslu_bld.f90 mlprec/mld_zslu_interface.c mlprec/mld_zslud_bld.f90 mlprec/mld_zsp_renum.f90 mlprec/mld_zumf_bld.f90 tests/newslv tests/newslv/Makefile tests/newslv/data_input.f90 tests/newslv/mld_d_tlu_solver.f90 tests/newslv/ppde.f90 tests/newslv/runs tests/newslv/runs/ppde.inp tests/newslv/spde.f90 tests/pdegen/ppde.f90 tests/pdegen/runs/ppde.inp Merged from newset branch.
This commit is contained in:
Salvatore Filippone
@@ -0,0 +1,39 @@
|
||||
MLDDIR=../..
|
||||
include $(MLDDIR)/Make.inc
|
||||
include $(MLDDIR)/Make.inc
|
||||
PSBLIBDIR=$(PSBLASDIR)/lib/
|
||||
PSBINCDIR=$(PSBLASDIR)/include
|
||||
MLDLIBDIR=$(MLDDIR)/lib
|
||||
MLD_LIB=-L$(MLDLIBDIR) -lpsb_krylov -lmld_prec -lpsb_prec
|
||||
PSBLAS_LIB= -L$(PSBLIBDIR) -lpsb_util -lpsb_base
|
||||
FINCLUDES=$(FMFLAG). $(FMFLAG)$(MLDLIBDIR) $(FMFLAG)$(PSBINCDIR) $(FIFLAG).
|
||||
|
||||
PDOBJS=ppde.o data_input.o mld_d_tlu_solver.o
|
||||
PSOBJS=spde.o data_input.o
|
||||
EXEDIR=./runs
|
||||
|
||||
all: ppde spde
|
||||
|
||||
ppde: $(PDOBJS)
|
||||
$(F90LINK) $(PDOBJS) -o ppde $(MLD_LIB) $(PSBLAS_LIB) $(LDLIBS)
|
||||
/bin/mv ppde $(EXEDIR)
|
||||
|
||||
|
||||
spde: $(PSOBJS)
|
||||
$(F90LINK) $(PSOBJS) -o spde $(MLD_LIB) $(PSBLAS_LIB) $(LDLIBS)
|
||||
/bin/mv spde $(EXEDIR)
|
||||
|
||||
ppde.o spde.o: data_input.o
|
||||
ppde.o: mld_d_tlu_solver.o
|
||||
|
||||
|
||||
clean:
|
||||
/bin/rm -f $(PDOBJS) $(EXEDIR)/ppde $(PSOBJS) $(EXEDIR)/spde
|
||||
|
||||
verycleanlib:
|
||||
(cd ../..; make veryclean)
|
||||
lib:
|
||||
(cd ../../; make library)
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,187 @@
|
||||
!!$
|
||||
!!$
|
||||
!!$ MLD2P4 version 2.0
|
||||
!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package
|
||||
!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0)
|
||||
!!$
|
||||
!!$ (C) Copyright 2008,2009,2010
|
||||
!!$
|
||||
!!$ 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.
|
||||
!!$
|
||||
!!$
|
||||
module data_input
|
||||
|
||||
interface read_data
|
||||
module procedure read_char, read_int,&
|
||||
& read_double, read_single,&
|
||||
& string_read_char, string_read_int,&
|
||||
& string_read_double, string_read_single
|
||||
end interface read_data
|
||||
interface trim_string
|
||||
module procedure trim_string
|
||||
end interface
|
||||
|
||||
character(len=4096), private :: charbuf
|
||||
character, private, parameter :: def_marker="!"
|
||||
|
||||
contains
|
||||
|
||||
subroutine read_char(val,file,marker)
|
||||
character(len=*), intent(out) :: val
|
||||
integer, intent(in) :: file
|
||||
character(len=1), optional, intent(in) :: marker
|
||||
|
||||
read(file,'(a)')charbuf
|
||||
call read_data(val,charbuf,marker)
|
||||
|
||||
end subroutine read_char
|
||||
|
||||
subroutine read_int(val,file,marker)
|
||||
integer, intent(out) :: val
|
||||
integer, intent(in) :: file
|
||||
character(len=1), optional, intent(in) :: marker
|
||||
|
||||
read(file,'(a)')charbuf
|
||||
call read_data(val,charbuf,marker)
|
||||
|
||||
end subroutine read_int
|
||||
subroutine read_single(val,file,marker)
|
||||
use psb_sparse_mod
|
||||
real(psb_spk_), intent(out) :: val
|
||||
integer, intent(in) :: file
|
||||
character(len=1), optional, intent(in) :: marker
|
||||
|
||||
read(file,'(a)')charbuf
|
||||
call read_data(val,charbuf,marker)
|
||||
|
||||
end subroutine read_single
|
||||
subroutine read_double(val,file,marker)
|
||||
use psb_sparse_mod
|
||||
real(psb_dpk_), intent(out) :: val
|
||||
integer, intent(in) :: file
|
||||
character(len=1), optional, intent(in) :: marker
|
||||
|
||||
read(file,'(a)')charbuf
|
||||
call read_data(val,charbuf,marker)
|
||||
|
||||
end subroutine read_double
|
||||
|
||||
subroutine string_read_char(val,file,marker)
|
||||
character(len=*), intent(out) :: val
|
||||
character(len=*), intent(in) :: file
|
||||
character(len=1), optional, intent(in) :: marker
|
||||
character(len=1) :: marker_
|
||||
character(len=1024) :: charbuf
|
||||
integer :: 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),'(a)') val
|
||||
end subroutine string_read_char
|
||||
|
||||
subroutine string_read_int(val,file,marker)
|
||||
integer, intent(out) :: val
|
||||
character(len=*), intent(in) :: file
|
||||
character(len=1), optional, intent(in) :: marker
|
||||
character(len=1) :: marker_
|
||||
character(len=1024) :: charbuf
|
||||
integer :: 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
|
||||
subroutine string_read_single(val,file,marker)
|
||||
use psb_sparse_mod
|
||||
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 :: 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
|
||||
subroutine string_read_double(val,file,marker)
|
||||
use psb_sparse_mod
|
||||
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 :: 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
|
||||
|
||||
function trim_string(string,marker)
|
||||
character(len=*), intent(in) :: string
|
||||
character(len=1), optional, intent(in) :: marker
|
||||
character(len=len(string)) :: trim_string
|
||||
character(len=1) :: marker_
|
||||
integer :: idx
|
||||
if (present(marker)) then
|
||||
marker_ = marker
|
||||
else
|
||||
marker_ = def_marker
|
||||
end if
|
||||
idx=index(string,marker_)
|
||||
trim_string = adjustl(string(idx:))
|
||||
end function trim_string
|
||||
end module data_input
|
||||
|
||||
@@ -0,0 +1,729 @@
|
||||
!!$
|
||||
!!$
|
||||
!!$ MLD2P4 version 2.0
|
||||
!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package
|
||||
!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0)
|
||||
!!$
|
||||
!!$ (C) Copyright 2008,2009,2010, 2010
|
||||
!!$
|
||||
!!$ 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.
|
||||
!!$
|
||||
!!$
|
||||
!
|
||||
!
|
||||
!
|
||||
!
|
||||
!
|
||||
!
|
||||
|
||||
module mld_d_tlu_solver
|
||||
|
||||
use mld_d_prec_type
|
||||
|
||||
type, extends(mld_d_base_solver_type) :: mld_d_tlu_solver_type
|
||||
type(psb_dspmat_type) :: l, u
|
||||
real(psb_dpk_), allocatable :: d(:)
|
||||
integer :: fact_type, fill_in
|
||||
real(psb_dpk_) :: thresh
|
||||
contains
|
||||
procedure, pass(sv) :: dump => d_tlu_solver_dmp
|
||||
procedure, pass(sv) :: build => d_tlu_solver_bld
|
||||
procedure, pass(sv) :: apply => 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
|
||||
procedure, pass(sv) :: setr => d_tlu_solver_setr
|
||||
procedure, pass(sv) :: descr => d_tlu_solver_descr
|
||||
procedure, pass(sv) :: sizeof => d_tlu_solver_sizeof
|
||||
procedure, pass(sv) :: default => d_tlu_solver_default
|
||||
end type mld_d_tlu_solver_type
|
||||
|
||||
|
||||
private :: d_tlu_solver_bld, d_tlu_solver_apply, &
|
||||
& d_tlu_solver_free, d_tlu_solver_seti, &
|
||||
& d_tlu_solver_setc, d_tlu_solver_setr,&
|
||||
& d_tlu_solver_descr, d_tlu_solver_sizeof, &
|
||||
& d_tlu_solver_default, d_tlu_solver_dmp
|
||||
|
||||
|
||||
interface mld_ilu0_fact
|
||||
subroutine mld_dilu0_fact(ialg,a,l,u,d,info,blck,upd)
|
||||
use psb_sparse_mod, only : psb_dspmat_type, psb_dpk_
|
||||
integer, intent(in) :: ialg
|
||||
integer, intent(out) :: info
|
||||
type(psb_dspmat_type),intent(in) :: a
|
||||
type(psb_dspmat_type),intent(inout) :: l,u
|
||||
type(psb_dspmat_type),intent(in), optional, target :: blck
|
||||
character, intent(in), optional :: upd
|
||||
real(psb_dpk_), intent(inout) :: d(:)
|
||||
end subroutine mld_dilu0_fact
|
||||
end interface
|
||||
|
||||
interface mld_iluk_fact
|
||||
subroutine mld_diluk_fact(fill_in,ialg,a,l,u,d,info,blck)
|
||||
use psb_sparse_mod, only : psb_dspmat_type, psb_dpk_
|
||||
integer, intent(in) :: fill_in,ialg
|
||||
integer, intent(out) :: info
|
||||
type(psb_dspmat_type),intent(in) :: a
|
||||
type(psb_dspmat_type),intent(inout) :: l,u
|
||||
type(psb_dspmat_type),intent(in), optional, target :: blck
|
||||
real(psb_dpk_), intent(inout) :: d(:)
|
||||
end subroutine mld_diluk_fact
|
||||
end interface
|
||||
|
||||
interface mld_ilut_fact
|
||||
subroutine mld_dilut_fact(fill_in,thres,a,l,u,d,info,blck)
|
||||
use psb_sparse_mod, only : psb_dspmat_type, psb_dpk_
|
||||
integer, intent(in) :: fill_in
|
||||
real(psb_dpk_), intent(in) :: thres
|
||||
integer, intent(out) :: info
|
||||
type(psb_dspmat_type),intent(in) :: a
|
||||
type(psb_dspmat_type),intent(inout) :: l,u
|
||||
type(psb_dspmat_type),intent(in), optional, target :: blck
|
||||
real(psb_dpk_), intent(inout) :: d(:)
|
||||
end subroutine mld_dilut_fact
|
||||
end interface
|
||||
|
||||
character(len=15), parameter, private :: &
|
||||
& fact_names(0:mld_slv_delta_+4)=(/&
|
||||
& 'none ','none ',&
|
||||
& 'none ','none ',&
|
||||
& 'none ','DIAG ?? ',&
|
||||
& 'ILU(n) ',&
|
||||
& 'MILU(n) ','ILU(t,n) '/)
|
||||
|
||||
|
||||
contains
|
||||
|
||||
subroutine d_tlu_solver_default(sv)
|
||||
|
||||
use psb_sparse_mod
|
||||
|
||||
Implicit None
|
||||
|
||||
! Arguments
|
||||
class(mld_d_tlu_solver_type), intent(inout) :: sv
|
||||
|
||||
sv%fact_type = mld_ilu_n_
|
||||
sv%fill_in = 0
|
||||
sv%thresh = dzero
|
||||
|
||||
return
|
||||
end subroutine d_tlu_solver_default
|
||||
|
||||
subroutine d_tlu_solver_check(sv,info)
|
||||
|
||||
use psb_sparse_mod
|
||||
|
||||
Implicit None
|
||||
|
||||
! Arguments
|
||||
class(mld_d_tlu_solver_type), intent(inout) :: sv
|
||||
integer, intent(out) :: info
|
||||
Integer :: err_act
|
||||
character(len=20) :: name='d_tlu_solver_check'
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
info = psb_success_
|
||||
|
||||
call mld_check_def(sv%fact_type,&
|
||||
& 'Factorization',mld_ilu_n_,is_legal_ilu_fact)
|
||||
|
||||
select case(sv%fact_type)
|
||||
case(mld_ilu_n_,mld_milu_n_)
|
||||
call mld_check_def(sv%fill_in,&
|
||||
& 'Level',0,is_legal_ml_lev)
|
||||
case(mld_ilu_t_)
|
||||
call mld_check_def(sv%thresh,&
|
||||
& 'Eps',dzero,is_legal_fact_thrs)
|
||||
end select
|
||||
|
||||
if (info /= psb_success_) goto 9999
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 continue
|
||||
call psb_erractionrestore(err_act)
|
||||
if (err_act == psb_act_abort_) then
|
||||
call psb_error()
|
||||
return
|
||||
end if
|
||||
return
|
||||
end subroutine d_tlu_solver_check
|
||||
|
||||
|
||||
subroutine d_tlu_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
|
||||
use psb_sparse_mod
|
||||
type(psb_desc_type), intent(in) :: desc_data
|
||||
class(mld_d_tlu_solver_type), intent(in) :: sv
|
||||
real(psb_dpk_),intent(in) :: x(:)
|
||||
real(psb_dpk_),intent(inout) :: y(:)
|
||||
real(psb_dpk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
real(psb_dpk_),target, intent(inout) :: work(:)
|
||||
integer, intent(out) :: info
|
||||
|
||||
integer :: n_row,n_col
|
||||
real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
|
||||
integer :: ictxt,np,me,i, err_act
|
||||
character :: trans_
|
||||
character(len=20) :: name='d_tlu_solver_apply'
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
info = psb_success_
|
||||
|
||||
trans_ = psb_toupper(trans)
|
||||
select case(trans_)
|
||||
case('N')
|
||||
case('T','C')
|
||||
case default
|
||||
call psb_errpush(psb_err_iarg_invalid_i_,name)
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
n_row = psb_cd_get_local_rows(desc_data)
|
||||
n_col = psb_cd_get_local_cols(desc_data)
|
||||
|
||||
if (n_col <= size(work)) then
|
||||
ww => work(1:n_col)
|
||||
if ((4*n_col+n_col) <= size(work)) then
|
||||
aux => work(n_col+1:)
|
||||
else
|
||||
allocate(aux(4*n_col),stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,i_err=(/4*n_col,0,0,0,0/),&
|
||||
& a_err='real(psb_dpk_)')
|
||||
goto 9999
|
||||
end if
|
||||
endif
|
||||
else
|
||||
allocate(ww(n_col),aux(4*n_col),stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,i_err=(/5*n_col,0,0,0,0/),&
|
||||
& a_err='real(psb_dpk_)')
|
||||
goto 9999
|
||||
end if
|
||||
endif
|
||||
|
||||
select case(trans_)
|
||||
case('N')
|
||||
call psb_spsm(done,sv%l,x,dzero,ww,desc_data,info,&
|
||||
& trans=trans_,scale='L',diag=sv%d,choice=psb_none_,work=aux)
|
||||
|
||||
if (info == psb_success_) call psb_spsm(alpha,sv%u,ww,beta,y,desc_data,info,&
|
||||
& trans=trans_,scale='U',choice=psb_none_, work=aux)
|
||||
|
||||
case('T','C')
|
||||
call psb_spsm(done,sv%u,x,dzero,ww,desc_data,info,&
|
||||
& trans=trans_,scale='L',diag=sv%d,choice=psb_none_,work=aux)
|
||||
if (info == psb_success_) call psb_spsm(alpha,sv%l,ww,beta,y,desc_data,info,&
|
||||
& trans=trans_,scale='U',choice=psb_none_,work=aux)
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Invalid TRANS in ILU subsolve')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
|
||||
if (info /= psb_success_) then
|
||||
|
||||
call psb_errpush(psb_err_internal_error_,name,&
|
||||
& a_err='Error in subsolve')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
if (n_col <= size(work)) then
|
||||
if ((4*n_col+n_col) <= size(work)) then
|
||||
else
|
||||
deallocate(aux)
|
||||
endif
|
||||
else
|
||||
deallocate(ww,aux)
|
||||
endif
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 continue
|
||||
call psb_erractionrestore(err_act)
|
||||
if (err_act == psb_act_abort_) then
|
||||
call psb_error()
|
||||
return
|
||||
end if
|
||||
return
|
||||
|
||||
end subroutine d_tlu_solver_apply
|
||||
|
||||
subroutine d_tlu_solver_bld(a,desc_a,sv,upd,info,b)
|
||||
|
||||
use psb_sparse_mod
|
||||
|
||||
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
|
||||
character, intent(in) :: upd
|
||||
integer, intent(out) :: info
|
||||
type(psb_dspmat_type), intent(in), target, optional :: b
|
||||
! Local variables
|
||||
integer :: n_row,n_col, nrow_a, nztota
|
||||
real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
|
||||
integer :: ictxt,np,me,i, err_act, debug_unit, debug_level
|
||||
character(len=20) :: name='d_tlu_solver_bld', ch_err
|
||||
|
||||
info=psb_success_
|
||||
call psb_erractionsave(err_act)
|
||||
debug_unit = psb_get_debug_unit()
|
||||
debug_level = psb_get_debug_level()
|
||||
ictxt = psb_cd_get_context(desc_a)
|
||||
call psb_info(ictxt, me, np)
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
|
||||
|
||||
n_row = psb_cd_get_local_rows(desc_a)
|
||||
|
||||
if (psb_toupper(upd) == 'F') then
|
||||
nrow_a = a%get_nrows()
|
||||
nztota = a%get_nzeros()
|
||||
if (present(b)) then
|
||||
nztota = nztota + b%get_nzeros()
|
||||
end if
|
||||
|
||||
call sv%l%csall(n_row,n_row,info,nztota)
|
||||
if (info == psb_success_) call sv%u%csall(n_row,n_row,info,nztota)
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psb_sp_all'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
if (allocated(sv%d)) then
|
||||
if (size(sv%d) < n_row) then
|
||||
deallocate(sv%d)
|
||||
endif
|
||||
endif
|
||||
if (.not.allocated(sv%d)) then
|
||||
allocate(sv%d(n_row),stat=info)
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
endif
|
||||
|
||||
|
||||
select case(sv%fact_type)
|
||||
|
||||
case (mld_ilu_t_)
|
||||
!
|
||||
! ILU(k,t)
|
||||
!
|
||||
select case(sv%fill_in)
|
||||
|
||||
case(:-1)
|
||||
! Error: fill-in <= -1
|
||||
call psb_errpush(psb_err_input_value_invalid_i_,&
|
||||
& name,i_err=(/3,sv%fill_in,0,0,0/))
|
||||
goto 9999
|
||||
|
||||
case(0:)
|
||||
! Fill-in >= 0
|
||||
call mld_ilut_fact(sv%fill_in,sv%thresh,&
|
||||
& a, sv%l,sv%u,sv%d,info,blck=b)
|
||||
end select
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='mld_ilut_fact'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
case(mld_ilu_n_,mld_milu_n_)
|
||||
!
|
||||
! ILU(k) and MILU(k)
|
||||
!
|
||||
select case(sv%fill_in)
|
||||
case(:-1)
|
||||
! Error: fill-in <= -1
|
||||
call psb_errpush(psb_err_input_value_invalid_i_,&
|
||||
& name,i_err=(/3,sv%fill_in,0,0,0/))
|
||||
goto 9999
|
||||
case(0)
|
||||
! Fill-in 0
|
||||
! Separate implementation of ILU(0) for better performance.
|
||||
! There seems to be a problem with the separate implementation of MILU(0),
|
||||
! contained into mld_ilu0_fact. This must be investigated. For the time being,
|
||||
! resort to the implementation of MILU(k) with k=0.
|
||||
if (sv%fact_type == mld_ilu_n_) then
|
||||
call mld_ilu0_fact(sv%fact_type,a,sv%l,sv%u,&
|
||||
& sv%d,info,blck=b,upd=upd)
|
||||
else
|
||||
call mld_iluk_fact(sv%fill_in,sv%fact_type,&
|
||||
& a,sv%l,sv%u,sv%d,info,blck=b)
|
||||
endif
|
||||
case(1:)
|
||||
! Fill-in >= 1
|
||||
! The same routine implements both ILU(k) and MILU(k)
|
||||
call mld_iluk_fact(sv%fill_in,sv%fact_type,&
|
||||
& a,sv%l,sv%u,sv%d,info,blck=b)
|
||||
end select
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='mld_iluk_fact'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
case default
|
||||
! If we end up here, something was wrong up in the call chain.
|
||||
info = psb_err_input_value_invalid_i_
|
||||
call psb_errpush(psb_err_input_value_invalid_i_,name,&
|
||||
& i_err=(/3,sv%fact_type,0,0,0/))
|
||||
goto 9999
|
||||
|
||||
end select
|
||||
else
|
||||
! Here we should add checks for reuse of L and U.
|
||||
! For the time being just throw an error.
|
||||
info = 31
|
||||
call psb_errpush(info, name, i_err=(/3,0,0,0,0/),a_err=upd)
|
||||
goto 9999
|
||||
|
||||
!
|
||||
! What is an update of a factorization??
|
||||
! A first attempt could be to reuse EXACTLY the existing indices
|
||||
! as if it was an ILU(0) (since, effectively, the sparsity pattern
|
||||
! should not grow beyond what is already there).
|
||||
!
|
||||
call mld_ilu0_fact(sv%fact_type,a,&
|
||||
& sv%l,sv%u,&
|
||||
& sv%d,info,blck=b,upd=upd)
|
||||
|
||||
end if
|
||||
|
||||
call sv%l%set_asb()
|
||||
call sv%l%trim()
|
||||
call sv%u%set_asb()
|
||||
call sv%u%trim()
|
||||
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' end'
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 continue
|
||||
call psb_erractionrestore(err_act)
|
||||
if (err_act == psb_act_abort_) then
|
||||
call psb_error()
|
||||
return
|
||||
end if
|
||||
return
|
||||
end subroutine d_tlu_solver_bld
|
||||
|
||||
|
||||
subroutine d_tlu_solver_seti(sv,what,val,info)
|
||||
|
||||
use psb_sparse_mod
|
||||
|
||||
Implicit None
|
||||
|
||||
! Arguments
|
||||
class(mld_d_tlu_solver_type), intent(inout) :: sv
|
||||
integer, intent(in) :: what
|
||||
integer, intent(in) :: val
|
||||
integer, intent(out) :: info
|
||||
Integer :: err_act
|
||||
character(len=20) :: name='d_tlu_solver_seti'
|
||||
|
||||
info = psb_success_
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
select case(what)
|
||||
case(mld_sub_solve_)
|
||||
sv%fact_type = val
|
||||
case(mld_sub_fillin_)
|
||||
sv%fill_in = val
|
||||
case default
|
||||
!!$ write(0,*) name,': Error: invalid WHAT'
|
||||
!!$ info = -2
|
||||
end select
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 continue
|
||||
call psb_erractionrestore(err_act)
|
||||
if (err_act == psb_act_abort_) then
|
||||
call psb_error()
|
||||
return
|
||||
end if
|
||||
return
|
||||
end subroutine d_tlu_solver_seti
|
||||
|
||||
subroutine d_tlu_solver_setc(sv,what,val,info)
|
||||
|
||||
use psb_sparse_mod
|
||||
|
||||
Implicit None
|
||||
|
||||
! Arguments
|
||||
class(mld_d_tlu_solver_type), intent(inout) :: sv
|
||||
integer, intent(in) :: what
|
||||
character(len=*), intent(in) :: val
|
||||
integer, intent(out) :: info
|
||||
Integer :: err_act, ival
|
||||
character(len=20) :: name='d_tlu_solver_setc'
|
||||
|
||||
info = psb_success_
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
|
||||
call mld_stringval(val,ival,info)
|
||||
if (info == psb_success_) call sv%set(what,ival,info)
|
||||
if (info /= psb_success_) then
|
||||
info = psb_err_from_subroutine_
|
||||
call psb_errpush(info, name)
|
||||
goto 9999
|
||||
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()
|
||||
return
|
||||
end if
|
||||
return
|
||||
end subroutine d_tlu_solver_setc
|
||||
|
||||
subroutine d_tlu_solver_setr(sv,what,val,info)
|
||||
|
||||
use psb_sparse_mod
|
||||
|
||||
Implicit None
|
||||
|
||||
! Arguments
|
||||
class(mld_d_tlu_solver_type), intent(inout) :: sv
|
||||
integer, intent(in) :: what
|
||||
real(psb_dpk_), intent(in) :: val
|
||||
integer, intent(out) :: info
|
||||
Integer :: err_act
|
||||
character(len=20) :: name='d_tlu_solver_setr'
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
info = psb_success_
|
||||
|
||||
select case(what)
|
||||
case(mld_sub_iluthrs_)
|
||||
sv%thresh = val
|
||||
case default
|
||||
!!$ write(0,*) name,': Error: invalid WHAT'
|
||||
!!$ info = -2
|
||||
!!$ goto 9999
|
||||
end select
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 continue
|
||||
call psb_erractionrestore(err_act)
|
||||
if (err_act == psb_act_abort_) then
|
||||
call psb_error()
|
||||
return
|
||||
end if
|
||||
return
|
||||
end subroutine d_tlu_solver_setr
|
||||
|
||||
subroutine d_tlu_solver_free(sv,info)
|
||||
|
||||
use psb_sparse_mod
|
||||
|
||||
Implicit None
|
||||
|
||||
! Arguments
|
||||
class(mld_d_tlu_solver_type), intent(inout) :: sv
|
||||
integer, intent(out) :: info
|
||||
Integer :: err_act
|
||||
character(len=20) :: name='d_tlu_solver_free'
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
info = psb_success_
|
||||
|
||||
if (allocated(sv%d)) then
|
||||
deallocate(sv%d,stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info = psb_err_alloc_dealloc_
|
||||
call psb_errpush(info,name)
|
||||
goto 9999
|
||||
end if
|
||||
end if
|
||||
call sv%l%free()
|
||||
call sv%u%free()
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 continue
|
||||
call psb_erractionrestore(err_act)
|
||||
if (err_act == psb_act_abort_) then
|
||||
call psb_error()
|
||||
return
|
||||
end if
|
||||
return
|
||||
end subroutine d_tlu_solver_free
|
||||
|
||||
subroutine d_tlu_solver_descr(sv,info,iout,coarse)
|
||||
|
||||
use psb_sparse_mod
|
||||
|
||||
Implicit None
|
||||
|
||||
! Arguments
|
||||
class(mld_d_tlu_solver_type), intent(in) :: sv
|
||||
integer, intent(out) :: info
|
||||
integer, intent(in), optional :: iout
|
||||
logical, intent(in), optional :: coarse
|
||||
|
||||
! Local variables
|
||||
integer :: err_act
|
||||
integer :: ictxt, me, np
|
||||
character(len=20), parameter :: name='mld_d_tlu_solver_descr'
|
||||
integer :: iout_
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
info = psb_success_
|
||||
if (present(iout)) then
|
||||
iout_ = iout
|
||||
else
|
||||
iout_ = 6
|
||||
endif
|
||||
|
||||
write(iout_,*) ' TLU: test a new solver kind'
|
||||
write(iout_,*) ' Incomplete factorization solver: ',&
|
||||
& fact_names(sv%fact_type)
|
||||
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
|
||||
end select
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 continue
|
||||
call psb_erractionrestore(err_act)
|
||||
if (err_act == psb_act_abort_) then
|
||||
call psb_error()
|
||||
return
|
||||
end if
|
||||
return
|
||||
end subroutine d_tlu_solver_descr
|
||||
|
||||
function d_tlu_solver_sizeof(sv) result(val)
|
||||
use psb_sparse_mod
|
||||
implicit none
|
||||
! Arguments
|
||||
class(mld_d_tlu_solver_type), intent(in) :: sv
|
||||
integer(psb_long_int_k_) :: val
|
||||
integer :: i
|
||||
|
||||
val = 2*psb_sizeof_int + psb_sizeof_dp
|
||||
if (allocated(sv%d)) val = val + psb_sizeof_dp * size(sv%d)
|
||||
val = val + psb_sizeof(sv%l)
|
||||
val = val + psb_sizeof(sv%u)
|
||||
|
||||
return
|
||||
end function d_tlu_solver_sizeof
|
||||
|
||||
subroutine d_tlu_solver_dmp(sv,ictxt,level,info,prefix,head,solver)
|
||||
use psb_sparse_mod
|
||||
implicit none
|
||||
class(mld_d_tlu_solver_type), intent(in) :: sv
|
||||
integer, intent(in) :: ictxt,level
|
||||
integer, intent(out) :: info
|
||||
character(len=*), intent(in), optional :: prefix, head
|
||||
logical, optional, intent(in) :: solver
|
||||
integer :: i, j, il1, iln, lname, lev
|
||||
integer :: icontxt,iam, np
|
||||
character(len=80) :: prefix_
|
||||
character(len=120) :: fname ! len should be at least 20 more than
|
||||
logical :: solver_
|
||||
! len of prefix_
|
||||
|
||||
info = 0
|
||||
|
||||
if (present(prefix)) then
|
||||
prefix_ = trim(prefix(1:min(len(prefix),len(prefix_))))
|
||||
else
|
||||
prefix_ = "dump_slv_d"
|
||||
end if
|
||||
|
||||
call psb_info(ictxt,iam,np)
|
||||
|
||||
if (present(solver)) then
|
||||
solver_ = solver
|
||||
else
|
||||
solver_ = .false.
|
||||
end if
|
||||
lname = len_trim(prefix_)
|
||||
fname = trim(prefix_)
|
||||
write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
|
||||
lname = lname + 5
|
||||
|
||||
if (solver_) then
|
||||
write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_lower.mtx'
|
||||
if (sv%l%is_asb()) &
|
||||
& call sv%l%print(fname,head=head)
|
||||
write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_diag.mtx'
|
||||
if (allocated(sv%d)) &
|
||||
& call psb_geprt(fname,sv%d,head=head)
|
||||
write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_upper.mtx'
|
||||
if (sv%u%is_asb()) &
|
||||
& call sv%u%print(fname,head=head)
|
||||
|
||||
end if
|
||||
|
||||
end subroutine d_tlu_solver_dmp
|
||||
|
||||
|
||||
end module mld_d_tlu_solver
|
||||
@@ -0,0 +1,751 @@
|
||||
!!$
|
||||
!!$
|
||||
!!$ MLD2P4 version 2.0
|
||||
!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package
|
||||
!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0)
|
||||
!!$
|
||||
!!$ (C) Copyright 2008,2009,2010
|
||||
!!$
|
||||
!!$ 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: ppde.f90
|
||||
!
|
||||
! Program: ppde
|
||||
! This sample program solves a linear system obtained by discretizing a
|
||||
! PDE with Dirichlet BCs.
|
||||
!
|
||||
!
|
||||
! The PDE is a general second order equation in 3d
|
||||
!
|
||||
! b1 dd(u) b2 dd(u) b3 dd(u) a1 d(u) a2 d(u) a3 d(u)
|
||||
! - ------ - ------ - ------ - ----- - ------ - ------ + a4 u = 0
|
||||
! dxdx dydy dzdz dx dy dz
|
||||
!
|
||||
! with Dirichlet boundary conditions, on the unit cube 0<=x,y,z<=1.
|
||||
!
|
||||
! Example taken from:
|
||||
! C.T.Kelley
|
||||
! Iterative Methods for Linear and Nonlinear Equations
|
||||
! SIAM 1995
|
||||
!
|
||||
! In this sample program the index space of the discretized
|
||||
! computational domain is first numbered sequentially in a standard way,
|
||||
! then the corresponding vector is distributed according to a BLOCK
|
||||
! data distribution.
|
||||
!
|
||||
! Boundary conditions are set in a very simple way, by adding
|
||||
! equations of the form
|
||||
!
|
||||
! u(x,y) = exp(-x^2-y^2-z^2)
|
||||
!
|
||||
! Note that if a1=a2=a3=a4=0., the PDE is the well-known Laplace equation.
|
||||
!
|
||||
program ppde
|
||||
use psb_sparse_mod
|
||||
use mld_prec_mod
|
||||
use psb_krylov_mod
|
||||
use psb_util_mod
|
||||
use data_input
|
||||
use mld_d_tlu_solver
|
||||
implicit none
|
||||
|
||||
! input parameters
|
||||
character(len=20) :: kmethd, ptype
|
||||
character(len=5) :: afmt
|
||||
integer :: 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
|
||||
type(mld_d_tlu_solver_type) :: tlusv
|
||||
! descriptor
|
||||
type(psb_desc_type) :: desc_a
|
||||
! dense matrices
|
||||
real(psb_dpk_), allocatable :: b(:), x(:)
|
||||
! 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_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
|
||||
! 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='pde90'
|
||||
call psb_set_errverbosity(2)
|
||||
|
||||
|
||||
!
|
||||
! 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 create_matrix(idim,a,b,x,desc_a,ictxt,afmt,info)
|
||||
call psb_barrier(ictxt)
|
||||
t2 = psb_wtime() - t1
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='create_matrix'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
|
||||
if (iam == psb_root_) write(*,'("Overall matrix creation time : ",es12.5)')t2
|
||||
if (iam == psb_root_) write(*,'(" ")')
|
||||
!
|
||||
! 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)
|
||||
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 mld_inner_precset(prec,tlusv,info,ilev=nlv)
|
||||
|
||||
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 psb_amx(ictxt,tprec)
|
||||
|
||||
if (iam == psb_root_) write(*,'("Preconditioner time : ",es12.5)')tprec
|
||||
if (iam == psb_root_) call mld_precdescr(prec,info)
|
||||
if (iam == psb_root_) write(*,'(" ")')
|
||||
|
||||
!
|
||||
! iterative method parameters
|
||||
!
|
||||
if(iam == psb_root_) write(*,'("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 = psb_sizeof(a)
|
||||
descsize = psb_sizeof(desc_a)
|
||||
precsize = mld_sizeof(prec)
|
||||
call psb_sum(ictxt,amatsize)
|
||||
call psb_sum(ictxt,descsize)
|
||||
call psb_sum(ictxt,precsize)
|
||||
if (iam == psb_root_) then
|
||||
write(*,'(" ")')
|
||||
write(*,'("Time to solve matrix : ",es12.5)')t2
|
||||
write(*,'("Time per iteration : ",es12.5)')t2/iter
|
||||
write(*,'("Number of iterations : ",i0)')iter
|
||||
write(*,'("Convergence indicator on exit : ",es12.5)')err
|
||||
write(*,'("Info on exit : ",i0)')info
|
||||
write(*,'("Total memory occupation for A: ",i12)')amatsize
|
||||
write(*,'("Total memory occupation for DESC_A: ",i12)')descsize
|
||||
write(*,'("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(*,'("Solving matrix : ell1")')
|
||||
write(*,'("Grid dimensions : ",i4,"x",i4,"x",i4)')idim,idim,idim
|
||||
write(*,'("Number of processors : ",i0)') np
|
||||
write(*,'("Data distribution : BLOCK")')
|
||||
write(*,'("Preconditioner : ",a)') prectype%descr
|
||||
write(*,'("Iterative method : ",a)') kmethd
|
||||
write(*,'(" ")')
|
||||
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
|
||||
|
||||
!
|
||||
! 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_sparse_mod
|
||||
implicit none
|
||||
integer :: idim
|
||||
integer, parameter :: nb=20
|
||||
real(psb_dpk_), allocatable :: 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
|
||||
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)
|
||||
|
||||
! 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(*,'("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(0,*) 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 = psb_cd_get_local_rows(desc_a)
|
||||
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)-a1(x,y,z)
|
||||
val(element) = val(element)/(deltah*&
|
||||
& deltah)
|
||||
zt(k) = exp(-y**2-z**2)*(-val(element))
|
||||
else
|
||||
val(element)=-b1(x,y,z)-a1(x,y,z)
|
||||
val(element) = val(element)/(deltah*&
|
||||
& deltah)
|
||||
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)-a2(x,y,z)
|
||||
val(element) = val(element)/(deltah*&
|
||||
& deltah)
|
||||
zt(k) = exp(-y**2-z**2)*exp(-x)*(-val(element))
|
||||
else
|
||||
val(element)=-b2(x,y,z)-a2(x,y,z)
|
||||
val(element) = val(element)/(deltah*deltah)
|
||||
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)-a3(x,y,z)
|
||||
val(element) = val(element)/(deltah*deltah)
|
||||
zt(k) = exp(-y**2-z**2)*exp(-x)*(-val(element))
|
||||
else
|
||||
val(element)=-b3(x,y,z)-a3(x,y,z)
|
||||
val(element) = val(element)/(deltah*deltah)
|
||||
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) + a1(x,y,z)&
|
||||
& + a2(x,y,z) + a3(x,y,z)
|
||||
val(element) = val(element)/(deltah*deltah)
|
||||
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)=-b1(x,y,z)
|
||||
val(element) = val(element)/(deltah*deltah)
|
||||
zt(k) = exp(-y**2-z**2)*exp(-x)*(-val(element))
|
||||
else
|
||||
val(element)=-b1(x,y,z)
|
||||
val(element) = val(element)/(deltah*deltah)
|
||||
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)
|
||||
val(element) = val(element)/(deltah*deltah)
|
||||
zt(k) = exp(-y**2-z**2)*exp(-x)*(-val(element))
|
||||
else
|
||||
val(element)=-b2(x,y,z)
|
||||
val(element) = val(element)/(deltah*deltah)
|
||||
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)=-b3(x,y,z)
|
||||
val(element) = val(element)/(deltah*deltah)
|
||||
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(*,'("The matrix has been generated and assembled in ",a3," format.")')&
|
||||
& ch_err(1:3)
|
||||
write(*,'("-allocation time : ",es12.5)') talc
|
||||
write(*,'("-coeff. gen. time : ",es12.5)') tgen
|
||||
write(*,'("-assembly time : ",es12.5)') tasb
|
||||
write(*,'("-total time : ",es12.5)') ttot
|
||||
|
||||
end if
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 continue
|
||||
call psb_erractionrestore(err_act)
|
||||
if (err_act == psb_act_abort_) then
|
||||
call psb_error(ictxt)
|
||||
return
|
||||
end if
|
||||
return
|
||||
end subroutine create_matrix
|
||||
end program ppde
|
||||
!
|
||||
! functions parametrizing the differential equation
|
||||
!
|
||||
function a1(x,y,z)
|
||||
use psb_sparse_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_sparse_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_sparse_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_sparse_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_sparse_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_sparse_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_sparse_mod, only : psb_dpk_
|
||||
real(psb_dpk_) :: b3
|
||||
real(psb_dpk_) :: x,y,z
|
||||
b3=1.d0
|
||||
end function b3
|
||||
|
||||
|
||||
@@ -0,0 +1,30 @@
|
||||
BICGSTAB ! Iterative method: BiCGSTAB BiCG CGS RGMRES BiCGSTABL CG
|
||||
CSR ! Storage format CSR COO JAD
|
||||
040 ! IDIM; domain size is idim**3
|
||||
2 ! ISTOPC
|
||||
0200 ! ITMAX
|
||||
-1 ! ITRACE
|
||||
30 ! IRST (restart for RGMRES and BiCGSTABL)
|
||||
1.d-6 ! EPS
|
||||
3L-M-RAS-I-D4 ! Longer descriptive name for preconditioner (up to 20 chars)
|
||||
ML ! Preconditioner NONE JACOBI BJAC AS ML
|
||||
0 ! Number of overlap layers for AS preconditioner at finest level
|
||||
HALO ! Restriction operator NONE HALO
|
||||
NONE ! Prolongation operator NONE SUM AVG
|
||||
ILU ! Subdomain solver DSCALE ILU MILU ILUT UMF SLU
|
||||
0 ! Level-set N for ILU(N)
|
||||
1.d-4 ! Threshold T for ILU(T,P)
|
||||
1 ! Smoother/Jacobi sweeps
|
||||
AS ! Smoother type JACOBI BJAC AS ignored for non-ML
|
||||
3 ! Number of levels in a multilevel preconditioner
|
||||
SMOOTHED ! Kind of aggregation: SMOOTHED, NONSMOOTHED, MINENERGY
|
||||
DEC ! Type of aggregation DEC SYMDEC GLB
|
||||
MULT ! Type of multilevel correction: ADD MULT
|
||||
POST ! Side of multiplicative correction PRE POST TWOSIDE (ignored for ADD)
|
||||
DIST ! Coarse level: matrix distribution DIST REPL
|
||||
BJAC ! Coarse level: solver JACOBI BJAC UMF SLU SLUDIST
|
||||
ILU ! Coarse level: subsolver DSCALE ILU UMF SLU SLUDIST
|
||||
1 ! Coarse level: Level-set N for ILU(N)
|
||||
1.d-4 ! Coarse level: Threshold T for ILU(T,P)
|
||||
4 ! Coarse level: Number of Jacobi sweeps
|
||||
0.10d0 ! Smoother Aggregation Threshold: >= 0.0
|
||||
@@ -0,0 +1,745 @@
|
||||
!!$
|
||||
!!$
|
||||
!!$ MLD2P4 version 2.0
|
||||
!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package
|
||||
!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0)
|
||||
!!$
|
||||
!!$ (C) Copyright 2008,2009,2010
|
||||
!!$
|
||||
!!$ 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: ppde.f90
|
||||
!
|
||||
! Program: ppde
|
||||
! This sample program solves a linear system obtained by discretizing a
|
||||
! PDE with Dirichlet BCs.
|
||||
!
|
||||
!
|
||||
! The PDE is a general second order equation in 3d
|
||||
!
|
||||
! b1 dd(u) b2 dd(u) b3 dd(u) a1 d(u) a2 d(u) a3 d(u)
|
||||
! - ------ - ------ - ------ - ----- - ------ - ------ + a4 u = 0
|
||||
! dxdx dydy dzdz dx dy dz
|
||||
!
|
||||
! with Dirichlet boundary conditions, on the unit cube 0<=x,y,z<=1.
|
||||
!
|
||||
! Example taken from:
|
||||
! C.T.Kelley
|
||||
! Iterative Methods for Linear and Nonlinear Equations
|
||||
! SIAM 1995
|
||||
!
|
||||
! In this sample program the index space of the discretized
|
||||
! computational domain is first numbered sequentially in a standard way,
|
||||
! then the corresponding vector is distributed according to a BLOCK
|
||||
! data distribution.
|
||||
!
|
||||
! Boundary conditions are set in a very simple way, by adding
|
||||
! equations of the form
|
||||
!
|
||||
! u(x,y) = exp(-x^2-y^2-z^2)
|
||||
!
|
||||
! Note that if a1=a2=a3=a4=0., the PDE is the well-known Laplace equation.
|
||||
!
|
||||
program spde
|
||||
use psb_sparse_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
|
||||
real(psb_spk_), allocatable :: b(:), x(:)
|
||||
! 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
|
||||
! 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='pde90'
|
||||
call psb_set_errverbosity(2)
|
||||
|
||||
!
|
||||
! 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 create_matrix(idim,a,b,x,desc_a,ictxt,afmt,info)
|
||||
call psb_barrier(ictxt)
|
||||
t2 = psb_wtime() - t1
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='create_matrix'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
if (iam == psb_root_) write(*,'("Overall matrix creation time : ",es12.5)')t2
|
||||
if (iam == psb_root_) write(*,'(" ")')
|
||||
!
|
||||
! 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)
|
||||
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 psb_amx(ictxt,tprec)
|
||||
|
||||
if (iam == psb_root_) write(*,'("Preconditioner time : ",es12.5)')tprec
|
||||
if (iam == psb_root_) call mld_precdescr(prec,info)
|
||||
if (iam == psb_root_) write(*,'(" ")')
|
||||
|
||||
!
|
||||
! iterative method parameters
|
||||
!
|
||||
if(iam == psb_root_) write(*,'("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 = psb_sizeof(a)
|
||||
descsize = psb_sizeof(desc_a)
|
||||
precsize = mld_sizeof(prec)
|
||||
call psb_sum(ictxt,amatsize)
|
||||
call psb_sum(ictxt,descsize)
|
||||
call psb_sum(ictxt,precsize)
|
||||
if (iam == psb_root_) then
|
||||
write(*,'(" ")')
|
||||
write(*,'("Time to solve matrix : ",es12.5)')t2
|
||||
write(*,'("Time per iteration : ",es12.5)')t2/iter
|
||||
write(*,'("Number of iterations : ",i0)')iter
|
||||
write(*,'("Convergence indicator on exit : ",es12.5)')err
|
||||
write(*,'("Info on exit : ",i0)')info
|
||||
write(*,'("Total memory occupation for A: ",i12)')amatsize
|
||||
write(*,'("Total memory occupation for DESC_A: ",i12)')descsize
|
||||
write(*,'("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(*,'("Solving matrix : ell1")')
|
||||
write(*,'("Grid dimensions : ",i4,"x",i4,"x",i4)')idim,idim,idim
|
||||
write(*,'("Number of processors : ",i0)') np
|
||||
write(*,'("Data distribution : BLOCK")')
|
||||
write(*,'("Preconditioner : ",a)') prectype%descr
|
||||
write(*,'("Iterative method : ",a)') kmethd
|
||||
write(*,'(" ")')
|
||||
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
|
||||
|
||||
!
|
||||
! 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_sparse_mod
|
||||
implicit none
|
||||
integer :: idim
|
||||
integer, parameter :: nb=20
|
||||
real(psb_spk_), allocatable :: 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
|
||||
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)
|
||||
|
||||
! 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(0,'("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(0,*) 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 = psb_cd_get_local_rows(desc_a)
|
||||
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)-a1(x,y,z)
|
||||
val(element) = val(element)/(deltah*&
|
||||
& deltah)
|
||||
zt(k) = exp(-y**2-z**2)*(-val(element))
|
||||
else
|
||||
val(element)=-b1(x,y,z)-a1(x,y,z)
|
||||
val(element) = val(element)/(deltah*&
|
||||
& deltah)
|
||||
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)-a2(x,y,z)
|
||||
val(element) = val(element)/(deltah*&
|
||||
& deltah)
|
||||
zt(k) = exp(-y**2-z**2)*exp(-x)*(-val(element))
|
||||
else
|
||||
val(element)=-b2(x,y,z)-a2(x,y,z)
|
||||
val(element) = val(element)/(deltah*deltah)
|
||||
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)-a3(x,y,z)
|
||||
val(element) = val(element)/(deltah*deltah)
|
||||
zt(k) = exp(-y**2-z**2)*exp(-x)*(-val(element))
|
||||
else
|
||||
val(element)=-b3(x,y,z)-a3(x,y,z)
|
||||
val(element) = val(element)/(deltah*deltah)
|
||||
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) + a1(x,y,z)&
|
||||
& + a2(x,y,z) + a3(x,y,z)
|
||||
val(element) = val(element)/(deltah*deltah)
|
||||
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)=-b1(x,y,z)
|
||||
val(element) = val(element)/(deltah*deltah)
|
||||
zt(k) = exp(-y**2-z**2)*exp(-x)*(-val(element))
|
||||
else
|
||||
val(element)=-b1(x,y,z)
|
||||
val(element) = val(element)/(deltah*deltah)
|
||||
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)
|
||||
val(element) = val(element)/(deltah*deltah)
|
||||
zt(k) = exp(-y**2-z**2)*exp(-x)*(-val(element))
|
||||
else
|
||||
val(element)=-b2(x,y,z)
|
||||
val(element) = val(element)/(deltah*deltah)
|
||||
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)=-b3(x,y,z)
|
||||
val(element) = val(element)/(deltah*deltah)
|
||||
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(*,'("The matrix has been generated and assembled in ",a3," format.")')&
|
||||
& ch_err(1:3)
|
||||
write(*,'("-allocation time : ",es12.5)') talc
|
||||
write(*,'("-coeff. gen. time : ",es12.5)') tgen
|
||||
write(*,'("-assembly time : ",es12.5)') tasb
|
||||
write(*,'("-total time : ",es12.5)') ttot
|
||||
|
||||
end if
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 continue
|
||||
call psb_erractionrestore(err_act)
|
||||
if (err_act == psb_act_abort_) then
|
||||
call psb_error(ictxt)
|
||||
return
|
||||
end if
|
||||
return
|
||||
end subroutine create_matrix
|
||||
end program spde
|
||||
!
|
||||
! functions parametrizing the differential equation
|
||||
!
|
||||
function a1(x,y,z)
|
||||
use psb_sparse_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_sparse_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_sparse_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_sparse_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_sparse_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_sparse_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_sparse_mod, only : psb_spk_
|
||||
real(psb_spk_) :: b3
|
||||
real(psb_spk_) :: x,y,z
|
||||
b3=1.e0
|
||||
end function b3
|
||||
|
||||
|
||||
@@ -217,7 +217,7 @@ program ppde
|
||||
end if
|
||||
|
||||
tprec = psb_wtime()-t1
|
||||
call prec%dump(info,prefix='test-ml',ac=.true.,solver=.true.,smoother=.true.)
|
||||
!!$ call prec%dump(info,prefix='test-ml',ac=.true.,solver=.true.,smoother=.true.)
|
||||
|
||||
call psb_amx(ictxt,tprec)
|
||||
|
||||
|
||||
@@ -1,6 +1,6 @@
|
||||
BICGSTAB ! Iterative method: BiCGSTAB BiCG CGS RGMRES BiCGSTABL CG
|
||||
CSR ! Storage format CSR COO JAD
|
||||
003 ! IDIM; domain size is idim**3
|
||||
040 ! IDIM; domain size is idim**3
|
||||
2 ! ISTOPC
|
||||
0002 ! ITMAX
|
||||
01 ! ITRACE
|
||||
@@ -15,16 +15,16 @@ ILU ! Subdomain solver DSCALE ILU MILU ILUT UMF SLU
|
||||
1 ! Level-set N for ILU(N)
|
||||
1.d-4 ! Threshold T for ILU(T,P)
|
||||
1 ! Smoother/Jacobi sweeps
|
||||
AS ! Smoother type JACOBI BJAC AS ignored for non-ML
|
||||
2 ! Number of levels in a multilevel preconditioner
|
||||
SMOOTHED ! Kind of aggregation: SMOOTHED, NONSMOOTHED, MINENERGY
|
||||
AS ! Smoother type JACOBI BJAC AS; ignored for non-ML
|
||||
3 ! Number of levels in a multilevel preconditioner
|
||||
SMOOTHED ! Kind of aggregation: SMOOTHED, NONSMOOTHED
|
||||
DEC ! Type of aggregation DEC SYMDEC GLB
|
||||
MULT ! Type of multilevel correction: ADD MULT
|
||||
POST ! Side of multiplicative correction PRE POST TWOSIDE (ignored for ADD)
|
||||
POST ! Side of correction PRE POST TWOSIDE (ignored for ADD)
|
||||
DIST ! Coarse level: matrix distribution DIST REPL
|
||||
BJAC ! Coarse level: solver JACOBI BJAC UMF SLU SLUDIST
|
||||
ILU ! Coarse level: subsolver DSCALE ILU UMF SLU SLUDIST
|
||||
0 ! Coarse level: Level-set N for ILU(N)
|
||||
UMF ! Coarse level: subsolver DSCALE ILU UMF SLU SLUDIST
|
||||
1 ! Coarse level: Level-set N for ILU(N)
|
||||
1.d-4 ! Coarse level: Threshold T for ILU(T,P)
|
||||
4 ! Coarse level: Number of Jacobi sweeps
|
||||
0.10d0 ! Smoother Aggregation Threshold: >= 0.0
|
||||
|
||||
Reference in New Issue
Block a user