You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
530 lines
18 KiB
Fortran
530 lines
18 KiB
Fortran
!
|
|
!
|
|
! AMG4PSBLAS version 1.0
|
|
! Algebraic Multigrid Package
|
|
! based on PSBLAS (Parallel Sparse BLAS version 3.7)
|
|
!
|
|
! (C) Copyright 2021
|
|
!
|
|
! Salvatore Filippone
|
|
! Pasqua D'Ambra
|
|
! Fabio Durastante
|
|
!
|
|
! Redistribution and use in source and binary forms, with or without
|
|
! modification, are permitted provided that the following conditions
|
|
! are met:
|
|
! 1. Redistributions of source code must retain the above copyright
|
|
! notice, this list of conditions and the following disclaimer.
|
|
! 2. Redistributions in binary form must reproduce the above copyright
|
|
! notice, this list of conditions, and the following disclaimer in the
|
|
! documentation and/or other materials provided with the distribution.
|
|
! 3. The name of the AMG4PSBLAS group or the names of its contributors may
|
|
! not be used to endorse or promote products derived from this
|
|
! software without specific written permission.
|
|
!
|
|
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
|
|
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
|
|
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
|
|
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
|
|
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
|
|
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
|
|
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
|
|
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
|
|
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
|
|
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
|
|
! POSSIBILITY OF SUCH DAMAGE.
|
|
!
|
|
module amg_d_decmatch_mod
|
|
|
|
use iso_c_binding
|
|
use psb_base_cbind_mod
|
|
|
|
interface new_Match_If
|
|
function dnew_Match_If(ipar,matching,lambda,nr, irp, ja, val, diag, w, mate) &
|
|
& bind(c,name="dnew_Match_If") result(res)
|
|
use iso_c_binding
|
|
import :: psb_c_ipk_, psb_c_lpk_, psb_c_mpk_, psb_c_epk_
|
|
implicit none
|
|
|
|
integer(psb_c_ipk_) :: res
|
|
integer(psb_c_ipk_), value :: nr,ipar,matching
|
|
real(c_double), value :: lambda
|
|
type(c_ptr), value :: irp, ja, mate
|
|
type(c_ptr), value :: val, diag, w
|
|
end function dnew_Match_If
|
|
end interface new_Match_If
|
|
|
|
interface amg_build_decmatch
|
|
module procedure amg_dbuild_decmatch
|
|
end interface amg_build_decmatch
|
|
|
|
logical, parameter, private :: print_statistics=.false.
|
|
contains
|
|
|
|
subroutine amg_ddecmatch_build_prol(w,a,desc_a,ilaggr,nlaggr,prol,info,&
|
|
& symmetrize,reproducible,display_inp, display_out, print_out, &
|
|
& parallel, matching,lambda)
|
|
use psb_base_mod
|
|
use psb_util_mod
|
|
use iso_c_binding
|
|
implicit none
|
|
real(psb_dpk_), allocatable, intent(inout) :: w(:)
|
|
type(psb_dspmat_type), intent(inout) :: a
|
|
type(psb_desc_type) :: desc_a
|
|
integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:)
|
|
integer(psb_lpk_), allocatable, intent(out) :: nlaggr(:)
|
|
type(psb_ldspmat_type), intent(out) :: prol
|
|
integer(psb_ipk_), intent(out) :: info
|
|
logical, optional, intent(in) :: display_inp, display_out, reproducible
|
|
logical, optional, intent(in) :: symmetrize, print_out, parallel
|
|
integer(psb_ipk_), optional, intent(in) :: matching
|
|
real(psb_dpk_), optional, intent(in) :: lambda
|
|
|
|
!
|
|
!
|
|
type(psb_ctxt_type) :: ictxt
|
|
integer(psb_ipk_) :: iam, np, iown
|
|
integer(psb_ipk_) :: nr, nc, sweep, nzl, ncsave, nct, idx
|
|
integer(psb_lpk_) :: i, k, kg, idxg, ntaggr, naggrm1, naggrp1, &
|
|
& ip, nlpairs, nlsingl, nunmatched, lnr
|
|
real(psb_dpk_) :: wk, widx, wmax, nrmagg
|
|
real(psb_dpk_), allocatable :: wtemp(:)
|
|
integer(psb_ipk_), allocatable :: mate(:)
|
|
integer(psb_lpk_), allocatable :: ilv(:)
|
|
integer(psb_ipk_), save :: cnt=1
|
|
character(len=256) :: aname
|
|
type(psb_ld_coo_sparse_mat) :: tmpcoo
|
|
logical :: display_out_, print_out_, reproducible_, parallel_
|
|
integer(psb_ipk_) :: matching_
|
|
real(psb_dpk_) :: lambda_
|
|
logical, parameter :: dump=.false., debug=.false., dump_mate=.false., &
|
|
& debug_ilaggr=.false., debug_sync=.false.
|
|
integer(psb_ipk_), save :: idx_bldmtc=-1, idx_phase1=-1, idx_phase2=-1, idx_phase3=-1
|
|
logical, parameter :: do_timings=.true.
|
|
|
|
ictxt = desc_a%get_ctxt()
|
|
call psb_info(ictxt,iam,np)
|
|
|
|
if ((do_timings).and.(idx_phase1==-1)) &
|
|
& idx_phase1 = psb_get_timer_idx("MBP_BLDP: phase1 ")
|
|
if ((do_timings).and.(idx_bldmtc==-1)) &
|
|
& idx_bldmtc = psb_get_timer_idx("MBP_BLDP: buil_matching")
|
|
if ((do_timings).and.(idx_phase2==-1)) &
|
|
& idx_phase2 = psb_get_timer_idx("MBP_BLDP: phase2 ")
|
|
if ((do_timings).and.(idx_phase3==-1)) &
|
|
& idx_phase3 = psb_get_timer_idx("MBP_BLDP: phase3 ")
|
|
|
|
if (do_timings) call psb_tic(idx_phase1)
|
|
|
|
if (present(display_out)) then
|
|
display_out_ = display_out
|
|
else
|
|
display_out_ = .false.
|
|
end if
|
|
if (present(print_out)) then
|
|
print_out_ = print_out
|
|
else
|
|
print_out_ = .false.
|
|
end if
|
|
if (present(reproducible)) then
|
|
reproducible_ = reproducible
|
|
else
|
|
reproducible_ = .false.
|
|
end if
|
|
|
|
if (present(parallel)) then
|
|
parallel_ = parallel
|
|
else
|
|
parallel_ = .true.
|
|
end if
|
|
|
|
if (present(matching)) then
|
|
matching_ = matching
|
|
else
|
|
matching_ = 2
|
|
end if
|
|
|
|
if (present(lambda)) then
|
|
lambda_ = lambda
|
|
else
|
|
lambda_ = 2.0
|
|
end if
|
|
|
|
allocate(nlaggr(0:np-1),stat=info)
|
|
if (info /= 0) then
|
|
return
|
|
end if
|
|
|
|
nlaggr = 0
|
|
ilv = [(i,i=1,desc_a%get_local_cols())]
|
|
call desc_a%l2gip(ilv,info,owned=.false.)
|
|
|
|
call psb_geall(ilaggr,desc_a,info)
|
|
ilaggr = -1
|
|
call psb_geasb(ilaggr,desc_a,info)
|
|
nr = a%get_nrows()
|
|
nc = a%get_ncols()
|
|
if (size(w) < nc) then
|
|
call psb_realloc(nc,w,info)
|
|
end if
|
|
call psb_halo(w,desc_a,info)
|
|
|
|
if (debug) write(0,*) iam,' buildprol into buildmatching:',&
|
|
& nr, nc
|
|
if (debug_sync) then
|
|
call psb_barrier(ictxt)
|
|
if (iam == 0) write(0,*)' buildprol into buildmatching:',&
|
|
& nr, nc
|
|
end if
|
|
if (do_timings) call psb_toc(idx_phase1)
|
|
if (do_timings) call psb_tic(idx_bldmtc)
|
|
call amg_dbuild_decmatch(parallel_,matching_,lambda_,w,a,desc_a,mate,info)
|
|
if (do_timings) call psb_toc(idx_bldmtc)
|
|
if (debug) write(0,*) iam,' buildprol from buildmatching:',&
|
|
& info
|
|
if (debug_sync) then
|
|
call psb_barrier(ictxt)
|
|
if (iam == 0) write(0,*)' out from buildmatching:', info
|
|
end if
|
|
|
|
if (info == 0) then
|
|
if (do_timings) call psb_tic(idx_phase2)
|
|
if (debug_sync) then
|
|
call psb_barrier(ictxt)
|
|
if (iam == 0) write(0,*)' Into building the tentative prol:'
|
|
end if
|
|
|
|
call psb_geall(wtemp,desc_a,info)
|
|
wtemp = dzero
|
|
call psb_geasb(wtemp,desc_a,info)
|
|
|
|
nlaggr(iam) = 0
|
|
nlpairs = 0
|
|
nlsingl = 0
|
|
nunmatched = 0
|
|
!
|
|
! First sweep
|
|
! On return from build_matching, mate has been converted to local numbering,
|
|
! so assigning to idx is OK.
|
|
!
|
|
do k=1, nr
|
|
idx = mate(k)
|
|
!
|
|
! Figure out an allocation of aggregates to processes
|
|
!
|
|
if (idx < 0) then
|
|
!
|
|
! Unmatched vertex, potential singleton.
|
|
!
|
|
nunmatched = nunmatched + 1
|
|
if (abs(w(k))<epsilon(nrmagg)) then
|
|
! Keep it unaggregated
|
|
wtemp(k) = dzero
|
|
else
|
|
! Create a singleton aggregate
|
|
nlaggr(iam) = nlaggr(iam) + 1
|
|
ilaggr(k) = nlaggr(iam)
|
|
wtemp(k) = w(k)/abs(w(k))
|
|
nlsingl = nlsingl + 1
|
|
end if
|
|
!!$ write(0,*) k,mate(k),ilaggr(k),' negative match ',abs(w(k)), epsilon(nrmagg)
|
|
else if (idx > nc) then
|
|
write(0,*) 'Impossible: mate(k) > nc'
|
|
cycle
|
|
else
|
|
|
|
if (ilaggr(k) == -1) then
|
|
|
|
wk = w(k)
|
|
widx = w(idx)
|
|
wmax = max(abs(wk),abs(widx))
|
|
nrmagg = wmax*sqrt((wk/wmax)**2+(widx/wmax)**2)
|
|
if (nrmagg > epsilon(nrmagg)) then
|
|
if (idx <= nr) then
|
|
if (ilaggr(idx) == -1) then
|
|
! Now, if both vertices are local, the aggregate is local
|
|
! (kinda obvious).
|
|
nlaggr(iam) = nlaggr(iam) + 1
|
|
ilaggr(k) = nlaggr(iam)
|
|
ilaggr(idx) = nlaggr(iam)
|
|
wtemp(k) = w(k)/nrmagg
|
|
wtemp(idx) = w(idx)/nrmagg
|
|
end if
|
|
nlpairs = nlpairs+1
|
|
else
|
|
write(0,*) 'Really? mate(k) > nr? ',mate(k),nr
|
|
end if
|
|
else
|
|
if (abs(w(k))<epsilon(nrmagg)) then
|
|
! Keep it unaggregated
|
|
wtemp(k) = dzero
|
|
else
|
|
! Create a singleton aggregate
|
|
nlaggr(iam) = nlaggr(iam) + 1
|
|
ilaggr(k) = nlaggr(iam)
|
|
wtemp(k) = w(k)/abs(w(k))
|
|
nlsingl = nlsingl + 1
|
|
end if
|
|
end if
|
|
end if
|
|
end if
|
|
end do
|
|
if (do_timings) call psb_toc(idx_phase2)
|
|
if (do_timings) call psb_tic(idx_phase3)
|
|
|
|
! Ok, now compute offsets, gather halo and fix non-local
|
|
! aggregates (those where ilaggr == -2)
|
|
call psb_sum(ictxt,nlaggr)
|
|
ntaggr = sum(nlaggr(0:np-1))
|
|
naggrm1 = sum(nlaggr(0:iam-1))
|
|
naggrp1 = sum(nlaggr(0:iam))
|
|
!
|
|
! Shift all indices already assigned (i.e. >0)
|
|
!
|
|
do k=1,nr
|
|
if (ilaggr(k) > 0) then
|
|
ilaggr(k) = ilaggr(k) + naggrm1
|
|
!!$ else
|
|
!!$ write(0,*) 'Leftover ILAGGR',k,ilaggr(k),mate(k),abs(w(k)),epsilon(nrmagg)
|
|
end if
|
|
end do
|
|
call psb_halo(ilaggr,desc_a,info)
|
|
call psb_halo(wtemp,desc_a,info)
|
|
! Cleanup as yet unmarked entries
|
|
do k=1,nr
|
|
if (ilaggr(k) == -2) then
|
|
idx = mate(k)
|
|
if (idx > nr) then
|
|
i = ilaggr(idx)
|
|
if (i > 0) then
|
|
ilaggr(k) = i
|
|
else
|
|
write(0,*) 'Error : unresolved (paired) index ',k,idx,i,nr,nc, ilv(k),ilv(idx)
|
|
end if
|
|
else
|
|
write(0,*) 'Error : unresolved (paired) index ',k,idx,i,nr,nc, ilv(k),ilv(idx)
|
|
end if
|
|
end if
|
|
if (ilaggr(k) <0) then
|
|
write(0,*) 'Decmatch: Funny number: ',k,ilv(k),ilaggr(k),wtemp(k)
|
|
end if
|
|
end do
|
|
if (debug_sync) then
|
|
call psb_barrier(ictxt)
|
|
if (iam == 0) write(0,*)' Done building the tentative prol:'
|
|
end if
|
|
|
|
|
|
if (dump_mate) then
|
|
block
|
|
integer(psb_lpk_), allocatable :: glaggr(:)
|
|
write(aname,'(a,i3.3,a,i3.3,a)') 'mateg-',cnt,'-p',iam,'.mtx'
|
|
open(20,file=aname)
|
|
write(20,'(a,I3,a)') '% sparse vector on process ',iam,' '
|
|
do k=1, nr
|
|
write(20,'(3(I8,1X))') ilv(k),ilv(mate(k))
|
|
end do
|
|
close(20)
|
|
write(aname,'(a,i3.3,a,i3.3,a)') 'nloc-',cnt,'-p',iam,'.mtx'
|
|
open(20,file=aname)
|
|
write(20,'(a,I3,a)') '% sparse vector on process ',iam,' '
|
|
write(20,'(a,I12,a)') 'nlpairs ',nlpairs
|
|
write(20,'(a,I12,a)') 'nlsingl ',nlsingl
|
|
write(20,'(a,I12,a)') 'nlaggr(iam) ',nlaggr(iam)
|
|
close(20)
|
|
write(aname,'(a,i3.3,a,i3.3,a,i3.3,a)') 'ilaggr-',cnt,'-i',iam,'-p',np,'.mtx'
|
|
open(20,file=aname)
|
|
write(20,'(a,I3,a)') '% sparse vector on process ',iam,' '
|
|
do k=1, nr
|
|
write(20,'(3(I8,1X))') ilv(k),ilaggr(k)
|
|
end do
|
|
close(20)
|
|
|
|
write(aname,'(a,i3.3,a,i3.3,a)') 'glaggr-',cnt,'-p',np,'.mtx'
|
|
call psb_gather(glaggr,ilaggr,desc_a,info,root=izero)
|
|
if (iam==0) call mm_array_write(glaggr,'Aggregates ',info,filename=aname)
|
|
|
|
cnt=cnt+1
|
|
end block
|
|
end if
|
|
block
|
|
integer(psb_lpk_) :: v(3)
|
|
v(1) = nunmatched
|
|
v(2) = nlsingl
|
|
v(3) = nlpairs
|
|
call psb_sum(ictxt,v)
|
|
nunmatched = v(1)
|
|
nlsingl = v(2)
|
|
nlpairs = v(3)
|
|
|
|
end block
|
|
if (print_statistics) then
|
|
if (iam == 0) then
|
|
write(0,*) 'Matching statistics: Unmatched nodes ',&
|
|
& nunmatched,' Singletons:',nlsingl,' Pairs:',nlpairs
|
|
end if
|
|
end if
|
|
|
|
if (display_out_) then
|
|
block
|
|
integer(psb_ipk_) :: idx
|
|
!
|
|
! And finally print out
|
|
!
|
|
do i=0,np-1
|
|
call psb_barrier(ictxt)
|
|
if (iam == i) then
|
|
write(0,*) 'Process ', iam,' hosts aggregates: (',naggrm1+1,' : ',naggrp1-1,')'
|
|
do k=1, nr
|
|
idx = mate(k)
|
|
kg = ilv(k)
|
|
if (idx >0) then
|
|
idxg = ilv(idx)
|
|
else
|
|
idxg = -1
|
|
end if
|
|
if (idx < 0) then
|
|
write(0,*) kg,': singleton (',kg,' ( Proc',iam,') ) into aggregate => ', ilaggr(k)
|
|
else if (idx <= nr) then
|
|
write(0,*) kg,': paired with (',idxg,' ( Proc',iam,') ) into aggregate => ', ilaggr(k)
|
|
else
|
|
call desc_a%indxmap%qry_halo_owner(idx,iown,info)
|
|
write(0,*) kg,': paired with (',idxg,' ( Proc',iown,') ) into aggregate => ', ilaggr(k)
|
|
end if
|
|
end do
|
|
flush(0)
|
|
end if
|
|
end do
|
|
end block
|
|
end if
|
|
|
|
! Dirty trick: allocate tmpcoo with local
|
|
! number of aggregates, then change to ntaggr.
|
|
! Just to make sure the allocation is not global
|
|
lnr = nr
|
|
call tmpcoo%allocate(lnr,nlaggr(iam),lnr)
|
|
k = 0
|
|
do i=1,nr
|
|
!
|
|
! Note: at this point, a value ilaggr(i)<=0
|
|
! tags an unaggregated row, and it has to be
|
|
! left alone (i.e.: it should stay at fine level only)
|
|
!
|
|
if (ilaggr(i)>0) then
|
|
k = k + 1
|
|
tmpcoo%val(k) = wtemp(i)
|
|
tmpcoo%ia(k) = i
|
|
tmpcoo%ja(k) = ilaggr(i)
|
|
end if
|
|
end do
|
|
call tmpcoo%set_nzeros(k)
|
|
call tmpcoo%set_dupl(psb_dupl_add_)
|
|
call tmpcoo%set_sorted() ! This is now in row-major
|
|
|
|
if (display_out_) then
|
|
call psb_barrier(ictxt)
|
|
flush(0)
|
|
|
|
if (iam == 0) write(0,*) 'Prolongator: '
|
|
flush(0)
|
|
do i=0,np-1
|
|
call psb_barrier(ictxt)
|
|
if (iam == i) then
|
|
do k=1, nr
|
|
write(0,*) ilv(tmpcoo%ia(k)),tmpcoo%ja(k), tmpcoo%val(k)
|
|
end do
|
|
flush(0)
|
|
end if
|
|
end do
|
|
end if
|
|
|
|
call prol%mv_from(tmpcoo)
|
|
if (do_timings) call psb_toc(idx_phase3)
|
|
|
|
if (print_out_) then
|
|
write(aname,'(a,i3.3,a)') 'prol-g-',iam,'.mtx'
|
|
call prol%print(fname=aname,head='Test ',ivr=ilv)
|
|
write(aname,'(a,i3.3,a)') 'prol-',iam,'.mtx'
|
|
call prol%print(fname=aname,head='Test ')
|
|
end if
|
|
|
|
else
|
|
write(0,*) iam,' : error from Matching: ',info
|
|
end if
|
|
|
|
end subroutine amg_ddecmatch_build_prol
|
|
|
|
subroutine amg_dbuild_decmatch(parallel,matching,lambda,w,a,desc_a,mate,info)
|
|
use psb_base_mod
|
|
use psb_util_mod
|
|
use iso_c_binding
|
|
implicit none
|
|
logical, intent(in) :: parallel
|
|
integer(psb_ipk_), intent(in) :: matching
|
|
real(psb_dpk_), intent(in) :: lambda
|
|
real(psb_dpk_), target :: w(:)
|
|
type(psb_dspmat_type), intent(inout) :: a
|
|
type(psb_desc_type) :: desc_a
|
|
integer(psb_ipk_), allocatable, intent(out), target :: mate(:)
|
|
integer(psb_ipk_), intent(out) :: info
|
|
|
|
type(psb_d_csr_sparse_mat), target :: tcsr
|
|
real(psb_dpk_), allocatable, target :: diag(:)
|
|
real(psb_dpk_) :: ph0t,ph1t,ph2t
|
|
type(psb_ctxt_type) :: ictxt
|
|
integer(psb_ipk_) :: iam, np
|
|
integer(psb_ipk_) :: nr, nc, nz, i, nunmatch, ipar
|
|
integer(psb_ipk_), save :: cnt=2
|
|
logical, parameter :: debug=.false., dump_ahat=.false., debug_sync=.false.
|
|
logical, parameter :: old_style=.false., sort_minp=.true.
|
|
character(len=40) :: name='build_matching', fname
|
|
integer(psb_ipk_), save :: idx_cmboxp=-1, idx_bldahat=-1, idx_phase2=-1, idx_phase3=-1
|
|
logical, parameter :: do_timings=.true.
|
|
|
|
ictxt = desc_a%get_ctxt()
|
|
call psb_info(ictxt,iam,np)
|
|
|
|
nr = a%get_nrows()
|
|
call a%cp_to(tcsr)
|
|
call psb_realloc(nr,mate,info)
|
|
diag = a%get_diag(info)
|
|
if (parallel) then
|
|
ipar = 2
|
|
else
|
|
ipar = 1
|
|
end if
|
|
!
|
|
! Now call matching!
|
|
!
|
|
if (debug) write(0,*) iam,' buildmatching into NewMatch:'
|
|
if (do_timings) call psb_tic(idx_cmboxp)
|
|
info = dnew_Match_If(ipar,matching,lambda,nr,c_loc(tcsr%irp),c_loc(tcsr%ja),&
|
|
& c_loc(tcsr%val),c_loc(diag),c_loc(w),c_loc(mate))
|
|
if (do_timings) call psb_toc(idx_cmboxp)
|
|
if (debug) write(0,*) iam,' buildmatching from NewMatch:', info
|
|
if (debug_sync) then
|
|
call psb_max(ictxt,info)
|
|
if (iam == 0) write(0,*)' done NewMatch', info
|
|
end if
|
|
|
|
if (do_timings) call psb_tic(idx_phase3)
|
|
nunmatch = count(mate(1:nr)<=0)
|
|
! call psb_sum(ictxt,nunmatch)
|
|
!if (nunmatch /= 0) write(0,*) iam,' Unmatched nodes local imbalance ',nunmatch
|
|
! if (count(mate(1:nr)<0) /= nunmatch) write(0,*) iam,' Matching results ?',&
|
|
! & nunmatch, count(mate(1:nr)<0)
|
|
if (debug_sync) then
|
|
call psb_barrier(ictxt)
|
|
if (iam == 0) write(0,*)' done build_matching '
|
|
end if
|
|
|
|
if (do_timings) call psb_toc(idx_phase3)
|
|
return
|
|
|
|
9999 continue
|
|
call psb_error(ictxt)
|
|
|
|
end subroutine amg_dbuild_decmatch
|
|
|
|
end module amg_d_decmatch_mod
|