Move newslv to samples and fix it

cmake
sfilippone 1 year ago
parent 94caab8aa9
commit 4a9f8e9da0

@ -0,0 +1,34 @@
AMGDIR=../..
AMGINCDIR=$(AMGDIR)/include
include $(AMGINCDIR)/Make.inc.amg4psblas
AMGMODDIR=$(AMGDIR)/modules
AMGLIBDIR=$(AMGDIR)/lib
AMG_LIBS=-L$(AMGLIBDIR) -lpsb_linsolve -lamg_prec -lpsb_prec
FINCLUDES=$(FMFLAG). $(FMFLAG)$(AMGMODDIR) $(FMFLAG)$(AMGINCDIR) $(PSBLAS_INCLUDES) $(FIFLAG).
PD3DOBJS=amg_pde3d_newslv.o data_input.o amg_d_tlu_solver.o amg_d_tlu_solver_impl.o
PSOBJS=spde.o data_input.o
EXEDIR=./runs
LINKOPT=
all: amg_pde3d_newslv
amg_pde3d_newslv: $(PD3DOBJS)
$(FLINK) $(LINKOPT) $(PD3DOBJS) -o amg_pde3d_newslv $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
/bin/mv amg_pde3d_newslv $(EXEDIR)
amg_pde3d_newslv.o amg_d_tlu_solver_impl.o: data_input.o amg_d_tlu_solver.o
clean:
/bin/rm -f $(PD3DOBJS) $(EXEDIR)/amg_pde3d_newslv
verycleanlib:
(cd ../..; make veryclean)
lib:
(cd ../../; make library)

@ -1,15 +1,15 @@
!
!
! MLD2P4 version 2.2
! MultiLevel Domain Decomposition Parallel Preconditioners Package
! based on PSBLAS (Parallel Sparse BLAS version 3.5)
!
! (C) Copyright 2008-2018
!
! Salvatore Filippone
! Pasqua D'Ambra
! Daniela di Serafino
!
!
!
! AMG4PSBLAS version 1.0
! Algebraic Multigrid Package
! based on PSBLAS (Parallel Sparse BLAS version 3.7)
!
! (C) Copyright 2021
!
! Salvatore Filippone
! Pasqua D'Ambra
! Fabio Durastante
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
@ -18,14 +18,14 @@
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the MLD2P4 group or the names of its contributors may
! 3. The name of the AMG4PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
@ -33,23 +33,22 @@
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
!
!
! File: mld_d_tlu_solver_mod.f90
!
! Module: mld_d_tlu_solver_mod
! File: amg_d_tlu_solver_mod.f90
!
! Module: amg_d_tlu_solver_mod
!
! This module serves as an example of how to define a new solver and integrate
! it in MLD2P4 via the P%SET(sv,info) method.
! it in AMG4PSBLAS via the P%SET(sv,info) method.
! In this example we are extending the ILU solver by implementing a new factorization algorithm.
! In actual reality, we are just giving a new name to ILU(0), but this should be sufficient to show
! the basics.
!
! The code is divided in two files:
! 1. The interface file (this one)
! 2. The implementation file (mld_d_tlu_solver_impl.f90)
! 2. The implementation file (amg_d_tlu_solver_impl.f90)
!
! The separation between interface and implementation is an essential part of the
! object-oriented design. The most appropriate tool would be to have the implementation
@ -59,12 +58,12 @@
!
!
module mld_d_tlu_solver
module amg_d_tlu_solver
use mld_d_ilu_solver
! use mld_d_ilu_fact_mod
use amg_d_ilu_solver
! use amg_d_ilu_fact_mod
type, extends(mld_d_ilu_solver_type) :: mld_d_tlu_solver_type
type, extends(amg_d_ilu_solver_type) :: amg_d_tlu_solver_type
!
! These are already defined in the ILU solver type; since we
! are supposedly implementing a new factorization strategy, the
@ -82,12 +81,12 @@ module mld_d_tlu_solver
! in common among all possible ILU factorizations
!
!
!procedure, pass(sv) :: dump => mld_d_tlu_solver_dmp
!procedure, pass(sv) :: dump => amg_d_tlu_solver_dmp
!procedure, pass(sv) :: ccheck => d_tlu_solver_check
!procedure, pass(sv) :: clone => mld_d_tlu_solver_clone
!procedure, pass(sv) :: cnv => mld_d_tlu_solver_cnv
!procedure, pass(sv) :: apply_v => mld_d_tlu_solver_apply_vect
!procedure, pass(sv) :: apply_a => mld_d_tlu_solver_apply
!procedure, pass(sv) :: clone => amg_d_tlu_solver_clone
!procedure, pass(sv) :: cnv => amg_d_tlu_solver_cnv
!procedure, pass(sv) :: apply_v => amg_d_tlu_solver_apply_vect
!procedure, pass(sv) :: apply_a => amg_d_tlu_solver_apply
!procedure, pass(sv) :: free => d_tlu_solver_free
!procedure, pass(sv) :: seti => d_tlu_solver_seti
!procedure, pass(sv) :: setc => d_tlu_solver_setc
@ -106,28 +105,28 @@ module mld_d_tlu_solver
!
procedure, pass(sv) :: descr => d_tlu_solver_descr
procedure, pass(sv) :: default => d_tlu_solver_default
procedure, pass(sv) :: build => mld_d_tlu_solver_bld
procedure, pass(sv) :: build => amg_d_tlu_solver_bld
procedure, nopass :: get_fmt => d_tlu_solver_get_fmt
end type mld_d_tlu_solver_type
end type amg_d_tlu_solver_type
private :: d_tlu_solver_get_fmt, d_tlu_solver_descr, d_tlu_solver_default
interface
subroutine mld_d_tlu_solver_bld(a,desc_a,sv,info,b,amold,vmold, imold)
import :: psb_desc_type, mld_d_tlu_solver_type, psb_d_vect_type, psb_dpk_, &
subroutine amg_d_tlu_solver_bld(a,desc_a,sv,info,b,amold,vmold, imold)
import :: psb_desc_type, amg_d_tlu_solver_type, psb_d_vect_type, psb_dpk_, &
& psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,&
& psb_ipk_, psb_i_base_vect_type
implicit none
type(psb_dspmat_type), intent(in), target :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(mld_d_tlu_solver_type), intent(inout) :: sv
Type(psb_desc_type), Intent(inout) :: desc_a
class(amg_d_tlu_solver_type), intent(inout) :: sv
integer(psb_ipk_), intent(out) :: info
type(psb_dspmat_type), intent(in), target, optional :: b
class(psb_d_base_sparse_mat), intent(in), optional :: amold
class(psb_d_base_vect_type), intent(in), optional :: vmold
class(psb_i_base_vect_type), intent(in), optional :: imold
end subroutine mld_d_tlu_solver_bld
end subroutine amg_d_tlu_solver_bld
end interface
contains
@ -141,9 +140,9 @@ contains
Implicit None
! Arguments
class(mld_d_tlu_solver_type), intent(inout) :: sv
class(amg_d_tlu_solver_type), intent(inout) :: sv
sv%fact_type = mld_ilu_n_
sv%fact_type = amg_ilu_n_
sv%fill_in = 0
sv%thresh = dzero
@ -157,20 +156,22 @@ contains
val = "TLU solver"
end function d_tlu_solver_get_fmt
subroutine d_tlu_solver_descr(sv,info,iout,coarse)
subroutine d_tlu_solver_descr(sv,info,iout,coarse,prefix)
Implicit None
! Arguments
class(mld_d_tlu_solver_type), intent(in) :: sv
class(amg_d_tlu_solver_type), intent(in) :: sv
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: iout
logical, intent(in), optional :: coarse
character(len=*), intent(in), optional :: prefix
! Local variables
integer(psb_ipk_) :: err_act
character(len=20), parameter :: name='mld_d_tlu_solver_descr'
character(len=20), parameter :: name='amg_d_tlu_solver_descr'
integer(psb_ipk_) :: iout_
character(1024) :: prefix_
call psb_erractionsave(err_act)
info = psb_success_
@ -179,14 +180,20 @@ contains
else
iout_ = psb_out_unit
endif
if (present(prefix)) then
prefix_ = prefix
else
prefix_ = ""
end if
write(iout_,*) ' Incomplete factorization solver: New Factorization TLU '
write(iout_,*) trim(prefix_), ' Incomplete factorization solver: New Factorization TLU '
select case(sv%fact_type)
case(mld_ilu_n_,mld_milu_n_)
write(iout_,*) ' Fill level:',sv%fill_in
case(mld_ilu_t_)
write(iout_,*) ' Fill level:',sv%fill_in
write(iout_,*) ' Fill threshold :',sv%thresh
case(amg_ilu_n_,amg_milu_n_)
write(iout_,*) trim(prefix_), ' Fill level:',sv%fill_in
case(amg_ilu_t_)
write(iout_,*) trim(prefix_), ' Fill level:',sv%fill_in
write(iout_,*) trim(prefix_), ' Fill threshold :',sv%thresh
end select
call psb_erractionrestore(err_act)
@ -196,4 +203,4 @@ contains
return
end subroutine d_tlu_solver_descr
end module mld_d_tlu_solver
end module amg_d_tlu_solver

@ -1,15 +1,15 @@
!
!
! MLD2P4 version 2.2
! MultiLevel Domain Decomposition Parallel Preconditioners Package
! based on PSBLAS (Parallel Sparse BLAS version 3.5)
!
! (C) Copyright 2008-2018
!
! Salvatore Filippone
! Pasqua D'Ambra
! Daniela di Serafino
!
!
!
! AMG4PSBLAS version 1.0
! Algebraic Multigrid Package
! based on PSBLAS (Parallel Sparse BLAS version 3.7)
!
! (C) Copyright 2021
!
! Salvatore Filippone
! Pasqua D'Ambra
! Fabio Durastante
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
@ -18,14 +18,14 @@
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the MLD2P4 group or the names of its contributors may
! 3. The name of the AMG4PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
@ -33,36 +33,36 @@
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
!
! File: mld_d_tlu_solver_impl.f90
! File: amg_d_tlu_solver_impl.f90
!
! This is the implementation file corresponding to mld_d_tlu_solver_mod.
! This is the implementation file corresponding to amg_d_tlu_solver_mod.
!
! In this example we are extending the ILU solver; we pretend to have a new
! factorization method, but since we are only interested in the interfacing,
! we are simply giving a new name to ILU(0).
!
!
subroutine mld_d_tlu_solver_bld(a,desc_a,sv,info,b,amold,vmold)
subroutine amg_d_tlu_solver_bld(a,desc_a,sv,info,b,amold,vmold)
use psb_base_mod
use mld_d_tlu_solver, mld_protect_name => mld_d_tlu_solver_bld
use mld_d_ilu_fact_mod
use amg_d_tlu_solver, amg_protect_name => amg_d_tlu_solver_bld
Implicit None
! Arguments
type(psb_dspmat_type), intent(in), target :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(mld_d_tlu_solver_type), intent(inout) :: sv
Type(psb_desc_type), Intent(inout) :: desc_a
class(amg_d_tlu_solver_type), intent(inout) :: sv
integer, intent(out) :: info
type(psb_dspmat_type), intent(in), target, optional :: b
class(psb_d_base_sparse_mat), intent(in), optional :: amold
class(psb_d_base_vect_type), intent(in), optional :: vmold
! Local variables
integer :: n_row,n_col, nrow_a, nztota
integer :: ctxt,np,me,i, err_act, debug_unit, debug_level
integer :: np,me,i, err_act, debug_unit, debug_level
type(psb_ctxt_type) :: ctxt
character(len=20) :: name='d_tlu_solver_bld', ch_err
info=psb_success_
@ -105,7 +105,7 @@ subroutine mld_d_tlu_solver_bld(a,desc_a,sv,info,b,amold,vmold)
! Fill-in 0, simple implementation.
call mld_ilu0_fact(sv%fact_type,a,sv%l,sv%u,&
call psb_ilu0_fact(sv%fact_type,a,sv%l,sv%u,&
& sv%d,info,blck=b)
call sv%l%set_asb()
@ -128,5 +128,5 @@ subroutine mld_d_tlu_solver_bld(a,desc_a,sv,info,b,amold,vmold)
9999 call psb_error_handler(err_act)
return
end subroutine mld_d_tlu_solver_bld
end subroutine amg_d_tlu_solver_bld

File diff suppressed because it is too large Load Diff

@ -1,14 +1,14 @@
!
!
! MLD2P4 version 2.2
! MultiLevel Domain Decomposition Parallel Preconditioners Package
! based on PSBLAS (Parallel Sparse BLAS version 3.5)
! AMG4PSBLAS version 1.0
! Algebraic Multigrid Package
! based on PSBLAS (Parallel Sparse BLAS version 3.7)
!
! (C) Copyright 2008-2018
! (C) Copyright 2021
!
! Salvatore Filippone
! Pasqua D'Ambra
! Daniela di Serafino
! Fabio Durastante
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
@ -18,14 +18,14 @@
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the MLD2P4 group or the names of its contributors may
! 3. The name of the AMG4PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
@ -39,10 +39,12 @@ module data_input
use psb_base_mod, only : psb_spk_, psb_dpk_, psb_ipk_
interface read_data
module procedure read_char, read_int,&
& read_double, read_single, read_logical,&
& string_read_char, string_read_int,&
& string_read_double, string_read_single, &
module procedure read_char, read_int, read_int_array,&
& read_double, read_double_array, &
& read_single, read_single_array, read_logical,&
& string_read_char, string_read_int, string_read_int_array,&
& string_read_double, string_read_double_array,&
& string_read_single, string_read_single_array, &
& string_read_logical
end interface read_data
interface trim_string
@ -51,15 +53,33 @@ module data_input
character(len=4096), private :: charbuf
character, private, parameter :: def_marker="!"
character, private, parameter :: cmt_marker="%"
contains
subroutine get_buffer(file,buffer)
integer(psb_ipk_), intent(in) :: file
character(len=*), intent(inout) :: buffer
integer :: idx
do
read(file,'(a)',end=999) buffer
buffer = adjustl(buffer)
idx=index(charbuf,cmt_marker)
if (idx == 1 ) then
cycle
else
exit
end if
end do
999 continue
return
end subroutine get_buffer
subroutine read_logical(val,file,marker)
logical, intent(out) :: val
integer(psb_ipk_), intent(in) :: file
character(len=1), optional, intent(in) :: marker
read(file,'(a)')charbuf
call get_buffer(file,charbuf)
call read_data(val,charbuf,marker)
end subroutine read_logical
@ -69,7 +89,7 @@ contains
integer(psb_ipk_), intent(in) :: file
character(len=1), optional, intent(in) :: marker
read(file,'(a)')charbuf
call get_buffer(file,charbuf)
call read_data(val,charbuf,marker)
end subroutine read_char
@ -79,29 +99,61 @@ contains
integer(psb_ipk_), intent(in) :: file
character(len=1), optional, intent(in) :: marker
read(file,'(a)')charbuf
call get_buffer(file,charbuf)
call read_data(val,charbuf,marker)
end subroutine read_int
subroutine read_int_array(val,file,marker)
integer(psb_ipk_), intent(out) :: val(:)
integer(psb_ipk_), intent(in) :: file
character(len=1), optional, intent(in) :: marker
call get_buffer(file,charbuf)
call read_data(val,charbuf,marker)
end subroutine read_int_array
subroutine read_single(val,file,marker)
real(psb_spk_), intent(out) :: val
integer(psb_ipk_), intent(in) :: file
character(len=1), optional, intent(in) :: marker
read(file,'(a)')charbuf
call get_buffer(file,charbuf)
call read_data(val,charbuf,marker)
end subroutine read_single
subroutine read_single_array(val,file,marker)
real(psb_spk_), intent(out) :: val(:)
integer(psb_ipk_), intent(in) :: file
character(len=1), optional, intent(in) :: marker
call get_buffer(file,charbuf)
call read_data(val,charbuf,marker)
end subroutine read_single_array
subroutine read_double(val,file,marker)
real(psb_dpk_), intent(out) :: val
integer(psb_ipk_), intent(in) :: file
character(len=1), optional, intent(in) :: marker
read(file,'(a)')charbuf
call get_buffer(file,charbuf)
call read_data(val,charbuf,marker)
end subroutine read_double
subroutine read_double_array(val,file,marker)
real(psb_dpk_), intent(out) :: val(:)
integer(psb_ipk_), intent(in) :: file
character(len=1), optional, intent(in) :: marker
call get_buffer(file,charbuf)
call read_data(val,charbuf,marker)
end subroutine read_double_array
subroutine string_read_char(val,file,marker)
character(len=*), intent(out) :: val
character(len=*), intent(in) :: file
@ -140,6 +192,25 @@ contains
read(charbuf(1:idx-1),*) val
end subroutine string_read_int
subroutine string_read_int_array(val,file,marker)
integer(psb_ipk_), intent(out) :: val(:)
character(len=*), intent(in) :: file
character(len=1), optional, intent(in) :: marker
character(len=1) :: marker_
character(len=1024) :: charbuf
integer(psb_ipk_) :: idx
if (present(marker)) then
marker_ = marker
else
marker_ = def_marker
end if
read(file,'(a)')charbuf
charbuf = adjustl(charbuf)
idx=index(charbuf,marker_)
if (idx == 0) idx = len(charbuf)+1
read(charbuf(1:idx-1),*) val(:)
end subroutine string_read_int_array
subroutine string_read_single(val,file,marker)
real(psb_spk_), intent(out) :: val
character(len=*), intent(in) :: file
@ -159,6 +230,25 @@ contains
read(charbuf(1:idx-1),*) val
end subroutine string_read_single
subroutine string_read_single_array(val,file,marker)
real(psb_spk_), intent(out) :: val(:)
character(len=*), intent(in) :: file
character(len=1), optional, intent(in) :: marker
character(len=1) :: marker_
character(len=1024) :: charbuf
integer(psb_ipk_) :: idx
if (present(marker)) then
marker_ = marker
else
marker_ = def_marker
end if
read(file,'(a)')charbuf
charbuf = adjustl(charbuf)
idx=index(charbuf,marker_)
if (idx == 0) idx = len(charbuf)+1
read(charbuf(1:idx-1),*) val(:)
end subroutine string_read_single_array
subroutine string_read_double(val,file,marker)
real(psb_dpk_), intent(out) :: val
character(len=*), intent(in) :: file
@ -178,6 +268,25 @@ contains
read(charbuf(1:idx-1),*) val
end subroutine string_read_double
subroutine string_read_double_array(val,file,marker)
real(psb_dpk_), intent(out) :: val(:)
character(len=*), intent(in) :: file
character(len=1), optional, intent(in) :: marker
character(len=1) :: marker_
character(len=1024) :: charbuf
integer(psb_ipk_) :: idx
if (present(marker)) then
marker_ = marker
else
marker_ = def_marker
end if
read(file,'(a)')charbuf
charbuf = adjustl(charbuf)
idx=index(charbuf,marker_)
if (idx == 0) idx = len(charbuf)+1
read(charbuf(1:idx-1),*) val(:)
end subroutine string_read_double_array
subroutine string_read_logical(val,file,marker)
logical, intent(out) :: val
character(len=*), intent(in) :: file

@ -1,34 +0,0 @@
MLDDIR=../..
MLDINCDIR=$(MLDDIR)/include
include $(MLDINCDIR)/Make.inc.amg4psblas
MLDMODDIR=$(MLDDIR)/modules
MLDLIBDIR=$(MLDDIR)/lib
MLD_LIBS=-L$(MLDLIBDIR) -lpsb_linsolve -lmld_prec -lpsb_prec
FINCLUDES=$(FMFLAG). $(FMFLAG)$(MLDMODDIR) $(FMFLAG)$(MLDINCDIR) $(PSBLAS_INCLUDES) $(FIFLAG).
PD3DOBJS=mld_pde3d_newslv.o data_input.o mld_d_tlu_solver.o mld_d_tlu_solver_impl.o
PSOBJS=spde.o data_input.o
EXEDIR=./runs
LINKOPT=
all: mld_pde3d_newslv
mld_pde3d_newslv: $(PD3DOBJS)
$(FLINK) $(LINKOPT) $(PD3DOBJS) -o mld_pde3d_newslv $(MLD_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
/bin/mv mld_pde3d_newslv $(EXEDIR)
mld_pde3d_newslv.o mld_d_tlu_solver_impl.o: data_input.o mld_d_tlu_solver.o
clean:
/bin/rm -f $(PD3DOBJS) $(EXEDIR)/mld_pde3d_newslv
verycleanlib:
(cd ../..; make veryclean)
lib:
(cd ../../; make library)

@ -1,385 +0,0 @@
!
!
! MLD2P4 version 2.2
! MultiLevel Domain Decomposition Parallel Preconditioners Package
! based on PSBLAS (Parallel Sparse BLAS version 3.5)
!
! (C) Copyright 2008-2018
!
! Salvatore Filippone
! Pasqua D'Ambra
! Daniela di Serafino
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the MLD2P4 group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
!
! File: mld_d_pde3d.f90
!
! Program: mld_d_pde3d
! This sample program solves a linear system obtained by discretizing a
! PDE with Dirichlet BCs.
!
!
! The PDE is a general second order equation in 3d
!
! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u)
! - ------ - ------ - ------ + ----- + ------ + ------ + c u = f
! dxdx dydy dzdz dx dy dz
!
! with Dirichlet boundary conditions
! u = g
!
! on the unit cube 0<=x,y,z<=1.
!
!
! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation.
!
! In this sample program the index space of the discretized
! computational domain is first numbered sequentially in a standard way,
! then the corresponding vector is distributed according to a BLOCK
! data distribution.
!
module mld_d_pde3d_mod
contains
!
! functions parametrizing the differential equation
!
function b1(x,y,z)
use psb_base_mod, only : psb_dpk_,done,dzero
real(psb_dpk_) :: b1
real(psb_dpk_), intent(in) :: x,y,z
b1=dzero/sqrt((3*done))
end function b1
function b2(x,y,z)
use psb_base_mod, only : psb_dpk_,done,dzero
real(psb_dpk_) :: b2
real(psb_dpk_), intent(in) :: x,y,z
b2=dzero/sqrt((3*done))
end function b2
function b3(x,y,z)
use psb_base_mod, only : psb_dpk_,done,dzero
real(psb_dpk_) :: b3
real(psb_dpk_), intent(in) :: x,y,z
b3=dzero/sqrt((3*done))
end function b3
function c(x,y,z)
use psb_base_mod, only : psb_dpk_,done,dzero
real(psb_dpk_) :: c
real(psb_dpk_), intent(in) :: x,y,z
c=dzero
end function c
function a1(x,y,z)
use psb_base_mod, only : psb_dpk_,done,dzero
real(psb_dpk_) :: a1
real(psb_dpk_), intent(in) :: x,y,z
a1=done!/80
end function a1
function a2(x,y,z)
use psb_base_mod, only : psb_dpk_,done,dzero
real(psb_dpk_) :: a2
real(psb_dpk_), intent(in) :: x,y,z
a2=done!/80
end function a2
function a3(x,y,z)
use psb_base_mod, only : psb_dpk_,done,dzero
real(psb_dpk_) :: a3
real(psb_dpk_), intent(in) :: x,y,z
a3=done!/80
end function a3
function g(x,y,z)
use psb_base_mod, only : psb_dpk_,done,dzero
real(psb_dpk_) :: g
real(psb_dpk_), intent(in) :: x,y,z
g = dzero
if (x == done) then
g = done
else if (x == dzero) then
g = exp(y**2-z**2)
end if
end function g
end module mld_d_pde3d_mod
program mld_d_pde3d
use psb_base_mod
use mld_prec_mod
use psb_krylov_mod
use psb_util_mod
use data_input
use mld_d_pde3d_mod
use mld_d_tlu_solver
implicit none
! input parameters
character(len=20) :: kmethd, ptype
character(len=5) :: afmt
integer(psb_ipk_) :: idim
! miscellaneous
real(psb_dpk_) :: t1, t2, tprec, thier, tslv
! sparse matrix and preconditioner
type(psb_dspmat_type) :: a
type(mld_dprec_type) :: prec
type(mld_d_tlu_solver_type) :: tlusv
! descriptor
type(psb_desc_type) :: desc_a
! dense vectors
type(psb_d_vect_type) :: x,b
! parallel environment
integer(psb_ipk_) :: ctxt, iam, np
! solver parameters
integer(psb_ipk_) :: iter, itmax,itrace, istopc, irst, nlv
integer(psb_epk_) :: amatsize, precsize, descsize
real(psb_dpk_) :: err, eps
! other variables
integer(psb_ipk_) :: info, i
character(len=20) :: name,ch_err
info=psb_success_
call psb_init(ctxt)
call psb_info(ctxt,iam,np)
if (iam < 0) then
! This should not happen, but just in case
call psb_exit(ctxt)
stop
endif
if(psb_get_errstatus() /= 0) goto 9999
name='mld_d_pde3d'
call psb_set_errverbosity(itwo)
!
! Hello world
!
if (iam == psb_root_) then
write(*,*) 'Welcome to MLD2P4 version: ',mld_version_string_
write(*,*) 'This is the ',trim(name),' sample program'
end if
!
! get parameters
!
call get_parms(ctxt,kmethd,afmt,idim,istopc,itmax,itrace,irst,eps)
!
! allocate and fill in the coefficient matrix, rhs and initial guess
!
call psb_barrier(ctxt)
t1 = psb_wtime()
call psb_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,&
& a1,a2,a3,b1,b2,b3,c,g,info)
call psb_barrier(ctxt)
t2 = psb_wtime() - t1
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='create_matrix'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (iam == psb_root_) &
& write(psb_out_unit,'("Overall matrix creation time : ",es12.5)')t2
if (iam == psb_root_) &
& write(psb_out_unit,'(" ")')
!
! prepare the preconditioner: an ML with defaults, but with TLU solver at
! intermediate levels. All other parameters are at default values.
!
call prec%init('ML', info)
call psb_barrier(ctxt)
t1 = psb_wtime()
call prec%hierarchy_build(a,desc_a,info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='prec%hierarchy_bld'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
thier = psb_wtime()-t1
nlv = prec%get_nlevs()
call prec%set(tlusv, info,ilev=1,ilmax=max(1,nlv-1))
call psb_barrier(ctxt)
t1 = psb_wtime()
call prec%smoothers_build(a,desc_a,info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='prec%smoothers_build'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
tprec = psb_wtime()-t1
call psb_amx(ctxt,thier)
call psb_amx(ctxt,tprec)
if (iam == psb_root_) &
& write(psb_out_unit,'("Preconditioner time : ",es12.5)') tprec+thier
if (iam == psb_root_) call prec%descr(info)
if (iam == psb_root_) &
& write(psb_out_unit,'(" ")')
!
! iterative method parameters
!
if(iam == psb_root_) &
& write(psb_out_unit,'("Calling iterative method ",a)')kmethd
call psb_barrier(ctxt)
t1 = psb_wtime()
call psb_krylov(kmethd,a,prec,b,x,eps,desc_a,info,&
& itmax=itmax,iter=iter,err=err,itrace=itrace,istop=istopc,irst=irst)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='solver routine'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call psb_barrier(ctxt)
tslv = psb_wtime() - t1
call psb_amx(ctxt,tslv)
amatsize = a%sizeof()
descsize = desc_a%sizeof()
precsize = prec%sizeof()
call psb_sum(ctxt,amatsize)
call psb_sum(ctxt,descsize)
call psb_sum(ctxt,precsize)
if (iam == psb_root_) then
write(psb_out_unit,'(" ")')
write(psb_out_unit,'("Numer of levels of aggr. hierarchy: ",i12)') prec%get_nlevs()
write(psb_out_unit,'("Time to build aggr. hierarchy : ",es12.5)') thier
write(psb_out_unit,'("Time to build smoothers : ",es12.5)') tprec
write(psb_out_unit,'("Total preconditioner time : ",es12.5)') tprec+thier
write(psb_out_unit,'("Time to solve system : ",es12.5)') tslv
write(psb_out_unit,'("Time per iteration : ",es12.5)') tslv/iter
write(psb_out_unit,'("Number of iterations : ",i0)') iter
write(psb_out_unit,'("Convergence indicator on exit : ",es12.5)') err
write(psb_out_unit,'("Info on exit : ",i0)') info
write(psb_out_unit,'("Total memory occupation for A: ",i12)') amatsize
write(psb_out_unit,'("Storage format for A: ",a)') trim(a%get_fmt())
write(psb_out_unit,'("Total memory occupation for DESC_A: ",i12)') descsize
write(psb_out_unit,'("Storage format for DESC_A: ",a)') trim(desc_a%get_fmt())
write(psb_out_unit,'("Total memory occupation for PREC: ",i12)') precsize
end if
!
! cleanup storage and exit
!
call psb_gefree(b,desc_a,info)
call psb_gefree(x,desc_a,info)
call psb_spfree(a,desc_a,info)
call prec%free(info)
call psb_cdfree(desc_a,info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='free routine'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call psb_exit(ctxt)
stop
9999 continue
call psb_error(ctxt)
contains
!
! get iteration parameters from standard input
!
subroutine get_parms(ctxt,kmethd,afmt,idim,istopc,itmax,itrace,irst,eps)
integer(psb_ipk_) :: ctxt
character(len=*) :: kmethd, afmt
integer(psb_ipk_) :: idim, istopc,itmax,itrace,irst
integer(psb_ipk_) :: np, iam, info
real(psb_dpk_) :: eps
character(len=20) :: buffer
call psb_info(ctxt, iam, np)
if (iam == psb_root_) then
call read_data(kmethd,psb_inp_unit)
call read_data(afmt,psb_inp_unit)
call read_data(idim,psb_inp_unit)
call read_data(istopc,psb_inp_unit)
call read_data(itmax,psb_inp_unit)
call read_data(itrace,psb_inp_unit)
call read_data(irst,psb_inp_unit)
call read_data(eps,psb_inp_unit)
end if
! broadcast parameters to all processors
call psb_bcast(ctxt,kmethd)
call psb_bcast(ctxt,afmt)
call psb_bcast(ctxt,idim)
call psb_bcast(ctxt,istopc)
call psb_bcast(ctxt,itmax)
call psb_bcast(ctxt,itrace)
call psb_bcast(ctxt,irst)
call psb_bcast(ctxt,eps)
if (iam == psb_root_) then
write(psb_out_unit,'("Solving matrix : ell1")')
write(psb_out_unit,'("Grid dimensions : ",i4,"x",i4,"x",i4)')idim,idim,idim
write(psb_out_unit,'("Number of processors : ",i0)') np
write(psb_out_unit,'("Data distribution : BLOCK")')
write(psb_out_unit,'("Preconditioner : ",a)') 'ML-TLU'
write(psb_out_unit,'("Iterative method : ",a)') kmethd
write(psb_out_unit,'(" ")')
endif
return
end subroutine get_parms
!
! print an error message
!
subroutine pr_usage(iout)
integer(psb_ipk_) :: iout
write(iout,*)'incorrect parameter(s) found'
write(iout,*)' usage: mld_d_pde3d methd prec dim &
&[istop itmax itrace]'
write(iout,*)' where:'
write(iout,*)' methd: cgstab cgs rgmres bicgstabl'
write(iout,*)' prec : bjac diag none'
write(iout,*)' dim number of points along each axis'
write(iout,*)' the size of the resulting linear '
write(iout,*)' system is dim**3'
write(iout,*)' istop stopping criterion 1, 2 '
write(iout,*)' itmax maximum number of iterations [500] '
write(iout,*)' itrace <=0 (no tracing, default) or '
write(iout,*)' >= 1 do tracing every itrace'
write(iout,*)' iterations '
end subroutine pr_usage
end program mld_d_pde3d
Loading…
Cancel
Save