fdurastante
/
amg4psblas
mirror of https://github.com/sfilippone/amg4psblas.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Fixed target attribute on A.

Browse Source
stopcriterion
Salvatore Filippone 18 years ago
parent 0436b0136f
commit 09830f6814
2 changed files with 433 additions and 363 deletions
Show Stats Download Patch File Download Diff File

407
mld_daggrmap_bld.f90
Unescape Escape View File

@ -40,24 +40,26 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
implicit none implicit none
integer, intent(in) :: aggr_type integer, intent(in) :: aggr_type
type(psb_dspmat_type), intent(in) :: a type(psb_dspmat_type), intent(in), target :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, allocatable :: ilaggr(:),nlaggr(:) integer, allocatable :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info integer, intent(out) :: info
! Locals ! Locals
integer, allocatable :: ils(:), neigh(:) integer, allocatable :: ils(:), neigh(:)
integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m
type(psb_dspmat_type), target :: atmp, atrans
type(psb_dspmat_type), pointer :: apnt
logical :: recovery logical :: recovery
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
integer ::ictxt,np,me,err_act integer :: ictxt,np,me,err_act
integer :: nrow, ncol, n_ne integer :: nrow, ncol, n_ne
integer, parameter :: one=1, two=2 integer, parameter :: one=1, two=2
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
if(psb_get_errstatus().ne.0) return if(psb_get_errstatus().ne.0) return
info=0 info=0
name = 'mld_daggrmap_bld' name = 'mld_aggrmap_bld'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
! !
! Note. At the time being we are ignoring aggr_type ! Note. At the time being we are ignoring aggr_type
@ -69,59 +71,143 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
nrow = psb_cd_get_local_rows(desc_a) nrow = psb_cd_get_local_rows(desc_a)
ncol = psb_cd_get_local_cols(desc_a) ncol = psb_cd_get_local_cols(desc_a)
nr = a%m select case (aggr_type)
allocate(ilaggr(nr),neigh(nr),stat=info) case (dec_aggr_,sym_dec_aggr_)
if(info.ne.0) then
info=4025
call psb_errpush(info,name,i_err=(/2*nr,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
do i=1, nr
ilaggr(i) = -(nr+1)
end do
! Note: -(nr+1) Untouched as yet
! -i 1<=i<=nr Adjacent to aggregate i
! i 1<=i<=nr Belonging to aggregate i
! nr = a%m
! Phase one: group nodes together. allocate(ilaggr(nr),neigh(nr),stat=info)
! Very simple minded strategy. if(info.ne.0) then
! info=4025
naggr = 0 call psb_errpush(info,name,i_err=(/2*nr,0,0,0,0/),&
nlp = 0 & a_err='integer')
do goto 9999
icnt = 0 end if
do i=1, nr do i=1, nr
if (ilaggr(i) == -(nr+1)) then ilaggr(i) = -(nr+1)
! end do
! 1. Untouched nodes are marked >0 together if (aggr_type == dec_aggr_) then
! with their neighbours apnt => a
! else
icnt = icnt + 1 call psb_sp_clip(a,atmp,info,imax=nr,jmax=nr,&
naggr = naggr + 1 & rscale=.false.,cscale=.false.)
ilaggr(i) = naggr atmp%m=nr
atmp%k=nr
call psb_transp(atmp,atrans,fmt='COO')
call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.)
atmp%m=nr
atmp%k=nr
call psb_sp_free(atrans,info)
call psb_ipcoo2csr(atmp,info)
apnt => atmp
end if
call psb_neigh(a,i,neigh,n_ne,info,lev=one)
if (info/=0) then ! Note: -(nr+1) Untouched as yet
info=4010 ! -i 1<=i<=nr Adjacent to aggregate i
ch_err='psb_neigh' ! i 1<=i<=nr Belonging to aggregate i
call psb_errpush(info,name,a_err=ch_err)
goto 9999 !
end if ! Phase one: group nodes together.
do k=1, n_ne ! Very simple minded strategy.
j = neigh(k) !
if ((1<=j).and.(j<=nr)) then naggr = 0
ilaggr(j) = naggr nlp = 0
do
icnt = 0
do i=1, nr
if (ilaggr(i) == -(nr+1)) then
!
! 1. Untouched nodes are marked >0 together
! with their neighbours
!
icnt = icnt + 1
naggr = naggr + 1
ilaggr(i) = naggr
call psb_neigh(apnt,i,neigh,n_ne,info,lev=one)
if (info/=0) then
info=4010
ch_err='psb_neigh'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
do k=1, n_ne
j = neigh(k)
if ((1<=j).and.(j<=nr)) then
ilaggr(j) = naggr
!!$ if (ilaggr(j) < 0) ilaggr(j) = naggr !!$ if (ilaggr(j) < 0) ilaggr(j) = naggr
!!$ if (ilaggr(j) == -(nr+1)) ilaggr(j) = naggr !!$ if (ilaggr(j) == -(nr+1)) ilaggr(j) = naggr
endif endif
enddo enddo
!
! 2. Untouched neighbours of these nodes are marked <0.
!
call psb_neigh(apnt,i,neigh,n_ne,info,lev=two)
if (info/=0) then
info=4010
ch_err='psb_neigh'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
do n = 1, n_ne
m = neigh(n)
if ((1<=m).and.(m<=nr)) then
if (ilaggr(m) == -(nr+1)) ilaggr(m) = -naggr
endif
enddo
endif
enddo
nlp = nlp + 1
if (icnt == 0) exit
enddo
if (debug) then
write(0,*) 'Check 1:',count(ilaggr == -(nr+1)),(a%ia1(i),i=a%ia2(1),a%ia2(2)-1)
end if
!
! Phase two: sweep over leftovers.
!
allocate(ils(naggr+10),stat=info)
if(info.ne.0) then
info=4025
call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
do i=1, size(ils)
ils(i) = 0
end do
do i=1, nr
n = ilaggr(i)
if (n>0) then
if (n>naggr) then
write(0,*) 'loc_Aggregate: n > naggr 1 ? ',n,naggr
else
ils(n) = ils(n) + 1
end if
end if
end do
if (debug) then
write(0,*) 'Phase 1: number of aggregates ',naggr
write(0,*) 'Phase 1: nodes aggregated ',sum(ils)
end if
recovery=.false.
do i=1, nr
if (ilaggr(i) < 0) then
! !
! 2. Untouched neighbours of these nodes are marked <0. ! Now some silly rule to break ties:
! Group with smallest adjacent aggregate.
! !
call psb_neigh(a,i,neigh,n_ne,info,lev=two) isz = nr+1
ia = -1
call psb_neigh(apnt,i,neigh,n_ne,info,lev=one)
if (info/=0) then if (info/=0) then
info=4010 info=4010
ch_err='psb_neigh' ch_err='psb_neigh'
@ -129,158 +215,107 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
goto 9999 goto 9999
end if end if
do n = 1, n_ne do j=1, n_ne
m = neigh(n) k = neigh(j)
if ((1<=m).and.(m<=nr)) then if ((1<=k).and.(k<=nr)) then
if (ilaggr(m) == -(nr+1)) ilaggr(m) = -naggr n = ilaggr(k)
if (n>0) then
if (n>naggr) then
write(0,*) 'loc_Aggregate: n > naggr 2? ',n,naggr
end if
if (ils(n) < isz) then
ia = n
isz = ils(n)
endif
endif
endif endif
enddo enddo
endif if (ia == -1) then
enddo if (ilaggr(i) > -(nr+1)) then
nlp = nlp + 1 ilaggr(i) = abs(ilaggr(i))
if (icnt == 0) exit if (ilaggr(I)>naggr) then
enddo write(0,*) 'loc_Aggregate: n > naggr 3? ',ilaggr(i),naggr
if (debug) then end if
write(0,*) 'Check 1:',count(ilaggr == -(nr+1)),(a%ia1(i),i=a%ia2(1),a%ia2(2)-1) ils(ilaggr(i)) = ils(ilaggr(i)) + 1
end if !
! This might happen if the pattern is non symmetric.
! ! Need a better handling.
! Phase two: sweep over leftovers. !
! recovery = .true.
allocate(ils(naggr+10),stat=info) else
if(info.ne.0) then write(0,*) 'Unrecoverable error !!',ilaggr(i), nr
info=4025 endif
call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),& else
& a_err='integer') ilaggr(i) = ia
goto 9999 if (ia>naggr) then
end if write(0,*) 'loc_Aggregate: n > naggr 4? ',ia,naggr
end if
do i=1, size(ils) ils(ia) = ils(ia) + 1
ils(i) = 0 endif
end do
do i=1, nr
n = ilaggr(i)
if (n>0) then
if (n>naggr) then
write(0,*) 'loc_Aggregate: n > naggr 1 ? ',n,naggr
else
ils(n) = ils(n) + 1
end if end if
enddo
if (recovery) then
write(0,*) 'Had to recover from strange situation in loc_aggregate.'
write(0,*) 'Perhaps an unsymmetric pattern?'
endif
if (debug) then
write(0,*) 'Phase 2: number of aggregates ',naggr
write(0,*) 'Phase 2: nodes aggregated ',sum(ils)
do i=1, naggr
write(*,*) 'Size of aggregate ',i,' :',count(ilaggr==i), ils(i)
enddo
write(*,*) maxval(ils(1:naggr))
write(*,*) 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops'
end if end if
end do
if (debug) then
write(0,*) 'Phase 1: number of aggregates ',naggr
write(0,*) 'Phase 1: nodes aggregated ',sum(ils)
end if
recovery=.false. !!$ write(0,*) 'desc_a loc_aggr 4 : ', desc_a%matrix_data(m_)
do i=1, nr if (count(ilaggr<0) >0) then
if (ilaggr(i) < 0) then write(0,*) 'Fatal error: some leftovers!!!'
! endif
! Now some silly rule to break ties:
! Group with smallest adjacent aggregate.
!
isz = nr+1
ia = -1
call psb_neigh(a,i,neigh,n_ne,info,lev=one)
if (info/=0) then
info=4010
ch_err='psb_neigh'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
do j=1, n_ne
k = neigh(j)
if ((1<=k).and.(k<=nr)) then
n = ilaggr(k)
if (n>0) then
if (n>naggr) then
write(0,*) 'loc_Aggregate: n > naggr 2? ',n,naggr
end if
if (ils(n) < isz) then deallocate(ils,neigh,stat=info)
ia = n if (info/=0) then
isz = ils(n) info=4000
endif call psb_errpush(info,name)
endif goto 9999
endif end if
enddo
if (ia == -1) then
if (ilaggr(i) > -(nr+1)) then
ilaggr(i) = abs(ilaggr(i))
if (ilaggr(I)>naggr) then
write(0,*) 'loc_Aggregate: n > naggr 3? ',ilaggr(i),naggr
end if
ils(ilaggr(i)) = ils(ilaggr(i)) + 1
!
! This might happen if the pattern is non symmetric.
! Need a better handling.
!
recovery = .true.
else
write(0,*) 'Unrecoverable error !!',ilaggr(i), nr
endif
else
ilaggr(i) = ia
if (ia>naggr) then
write(0,*) 'loc_Aggregate: n > naggr 4? ',ia,naggr
end if
ils(ia) = ils(ia) + 1 if (nrow /= size(ilaggr)) then
endif write(0,*) 'SOmething wrong ilaggr ',nrow,size(ilaggr)
endif
call psb_realloc(ncol,ilaggr,info)
if (info/=0) then
info=4010
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if end if
enddo
if (recovery) then
write(0,*) 'Had to recover from strange situation in loc_aggregate.'
write(0,*) 'Perhaps an unsymmetric pattern?'
endif
if (debug) then
write(0,*) 'Phase 2: number of aggregates ',naggr
write(0,*) 'Phase 2: nodes aggregated ',sum(ils)
do i=1, naggr
write(*,*) 'Size of aggregate ',i,' :',count(ilaggr==i), ils(i)
enddo
write(*,*) maxval(ils(1:naggr))
write(*,*) 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops'
end if
!!$ write(0,*) 'desc_a loc_aggr 4 : ', desc_a%matrix_data(m_) allocate(nlaggr(np),stat=info)
if (count(ilaggr<0) >0) then if (info/=0) then
write(0,*) 'Fatal error: some leftovers!!!' info=4025
endif call psb_errpush(info,name,i_err=(/np,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
deallocate(ils,neigh,stat=info) nlaggr(:) = 0
if (info/=0) then nlaggr(me+1) = naggr
info=4000 call psb_sum(ictxt,nlaggr(1:np))
call psb_errpush(info,name)
goto 9999
end if
if (nrow /= size(ilaggr)) then if (aggr_type == sym_dec_aggr_) then
write(0,*) 'SOmething wrong ilaggr ',nrow,size(ilaggr) call psb_sp_free(atmp,info)
endif end if
call psb_realloc(ncol,ilaggr,info)
if (info/=0) then
info=4010
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
allocate(nlaggr(np),stat=info) case default
if (info/=0) then
info=4025
call psb_errpush(info,name,i_err=(/np,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
nlaggr(:) = 0 write(0,*) 'Unimplemented aggregation algorithm ',aggr_type
nlaggr(me+1) = naggr info = -1
call psb_sum(ictxt,nlaggr(1:np)) call psb_errpush(30,name,i_err=(/1,aggr_type,0,0,0/))
goto 9999
end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

387
mld_zaggrmap_bld.f90
Unescape Escape View File

@ -40,13 +40,15 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
implicit none implicit none
integer, intent(in) :: aggr_type integer, intent(in) :: aggr_type
type(psb_zspmat_type), intent(in) :: a type(psb_zspmat_type), intent(in), target :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, allocatable :: ilaggr(:),nlaggr(:) integer, allocatable :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info integer, intent(out) :: info
! Locals ! Locals
integer, allocatable :: ils(:), neigh(:) integer, allocatable :: ils(:), neigh(:)
integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m
type(psb_zspmat_type), target :: atmp, atrans
type(psb_zspmat_type), pointer :: apnt
logical :: recovery logical :: recovery
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
@ -57,7 +59,7 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
if(psb_get_errstatus().ne.0) return if(psb_get_errstatus().ne.0) return
info=0 info=0
name = 'mld_zaggrmap_bld' name = 'mld_aggrmap_bld'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
! !
! Note. At the time being we are ignoring aggr_type ! Note. At the time being we are ignoring aggr_type
@ -69,114 +71,135 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
nrow = psb_cd_get_local_rows(desc_a) nrow = psb_cd_get_local_rows(desc_a)
ncol = psb_cd_get_local_cols(desc_a) ncol = psb_cd_get_local_cols(desc_a)
nr = a%m select case (aggr_type)
allocate(ilaggr(nr),neigh(nr),stat=info) case (dec_aggr_,sym_dec_aggr_)
if(info.ne.0) then
info=4025
call psb_errpush(info,name,i_err=(/2*nr,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
do i=1, nr
ilaggr(i) = -(nr+1)
end do
! Note: -(nr+1) Untouched as yet
! -i 1<=i<=nr Adjacent to aggregate i
! i 1<=i<=nr Belonging to aggregate i
! nr = a%m
! Phase one: group nodes together. allocate(ilaggr(nr),neigh(nr),stat=info)
! Very simple minded strategy. if(info.ne.0) then
! info=4025
naggr = 0 call psb_errpush(info,name,i_err=(/2*nr,0,0,0,0/),&
nlp = 0 & a_err='integer')
do goto 9999
icnt = 0 end if
do i=1, nr
do i=1, nr
ilaggr(i) = -(nr+1)
end do
if (aggr_type == dec_aggr_) then
apnt => a
else
call psb_sp_clip(a,atmp,info,imax=nr,jmax=nr,&
& rscale=.false.,cscale=.false.)
atmp%m=nr
atmp%k=nr
call psb_transp(atmp,atrans,fmt='COO')
call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.)
atmp%m=nr
atmp%k=nr
call psb_sp_free(atrans,info)
call psb_ipcoo2csr(atmp,info)
apnt => atmp
end if
! Note: -(nr+1) Untouched as yet
! -i 1<=i<=nr Adjacent to aggregate i
! i 1<=i<=nr Belonging to aggregate i
!
! Phase one: group nodes together.
! Very simple minded strategy.
!
naggr = 0
nlp = 0
do
icnt = 0
do i=1, nr
if (ilaggr(i) == -(nr+1)) then if (ilaggr(i) == -(nr+1)) then
! !
! 1. Untouched nodes are marked >0 together ! 1. Untouched nodes are marked >0 together
! with their neighbours ! with their neighbours
! !
icnt = icnt + 1 icnt = icnt + 1
naggr = naggr + 1 naggr = naggr + 1
ilaggr(i) = naggr ilaggr(i) = naggr
call psb_neigh(a,i,neigh,n_ne,info,lev=one) call psb_neigh(apnt,i,neigh,n_ne,info,lev=one)
if (info/=0) then if (info/=0) then
info=4010 info=4010
ch_err='psb_neigh' ch_err='psb_neigh'
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
do k=1, n_ne do k=1, n_ne
j = neigh(k) j = neigh(k)
if ((1<=j).and.(j<=nr)) then if ((1<=j).and.(j<=nr)) then
ilaggr(j) = naggr ilaggr(j) = naggr
!!$ if (ilaggr(j) < 0) ilaggr(j) = naggr !!$ if (ilaggr(j) < 0) ilaggr(j) = naggr
!!$ if (ilaggr(j) == -(nr+1)) ilaggr(j) = naggr !!$ if (ilaggr(j) == -(nr+1)) ilaggr(j) = naggr
endif endif
enddo enddo
! !
! 2. Untouched neighbours of these nodes are marked <0. ! 2. Untouched neighbours of these nodes are marked <0.
! !
call psb_neigh(a,i,neigh,n_ne,info,lev=two) call psb_neigh(apnt,i,neigh,n_ne,info,lev=two)
if (info/=0) then if (info/=0) then
info=4010 info=4010
ch_err='psb_neigh' ch_err='psb_neigh'
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
do n = 1, n_ne do n = 1, n_ne
m = neigh(n) m = neigh(n)
if ((1<=m).and.(m<=nr)) then if ((1<=m).and.(m<=nr)) then
if (ilaggr(m) == -(nr+1)) ilaggr(m) = -naggr if (ilaggr(m) == -(nr+1)) ilaggr(m) = -naggr
endif endif
enddo enddo
endif endif
enddo enddo
nlp = nlp + 1 nlp = nlp + 1
if (icnt == 0) exit if (icnt == 0) exit
enddo enddo
if (debug) then if (debug) then
write(0,*) 'Check 1:',count(ilaggr == -(nr+1)),(a%ia1(i),i=a%ia2(1),a%ia2(2)-1) write(0,*) 'Check 1:',count(ilaggr == -(nr+1)),(a%ia1(i),i=a%ia2(1),a%ia2(2)-1)
end if end if
! !
! Phase two: sweep over leftovers. ! Phase two: sweep over leftovers.
! !
allocate(ils(naggr+10),stat=info) allocate(ils(naggr+10),stat=info)
if(info.ne.0) then if(info.ne.0) then
info=4025 info=4025
call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),& call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),&
& a_err='integer') & a_err='integer')
goto 9999 goto 9999
end if end if
do i=1, size(ils) do i=1, size(ils)
ils(i) = 0 ils(i) = 0
end do end do
do i=1, nr do i=1, nr
n = ilaggr(i) n = ilaggr(i)
if (n>0) then if (n>0) then
if (n>naggr) then if (n>naggr) then
write(0,*) 'loc_Aggregate: n > naggr 1 ? ',n,naggr write(0,*) 'loc_Aggregate: n > naggr 1 ? ',n,naggr
else else
ils(n) = ils(n) + 1 ils(n) = ils(n) + 1
end if end if
end if end if
end do end do
if (debug) then if (debug) then
write(0,*) 'Phase 1: number of aggregates ',naggr write(0,*) 'Phase 1: number of aggregates ',naggr
write(0,*) 'Phase 1: nodes aggregated ',sum(ils) write(0,*) 'Phase 1: nodes aggregated ',sum(ils)
end if end if
recovery=.false. recovery=.false.
do i=1, nr do i=1, nr
if (ilaggr(i) < 0) then if (ilaggr(i) < 0) then
! !
! Now some silly rule to break ties: ! Now some silly rule to break ties:
! Group with smallest adjacent aggregate. ! Group with smallest adjacent aggregate.
@ -184,103 +207,115 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
isz = nr+1 isz = nr+1
ia = -1 ia = -1
call psb_neigh(a,i,neigh,n_ne,info,lev=one) call psb_neigh(apnt,i,neigh,n_ne,info,lev=one)
if (info/=0) then if (info/=0) then
info=4010 info=4010
ch_err='psb_neigh' ch_err='psb_neigh'
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
do j=1, n_ne do j=1, n_ne
k = neigh(j) k = neigh(j)
if ((1<=k).and.(k<=nr)) then if ((1<=k).and.(k<=nr)) then
n = ilaggr(k) n = ilaggr(k)
if (n>0) then if (n>0) then
if (n>naggr) then if (n>naggr) then
write(0,*) 'loc_Aggregate: n > naggr 2? ',n,naggr write(0,*) 'loc_Aggregate: n > naggr 2? ',n,naggr
end if end if
if (ils(n) < isz) then if (ils(n) < isz) then
ia = n ia = n
isz = ils(n) isz = ils(n)
endif
endif endif
endif endif
endif
enddo enddo
if (ia == -1) then if (ia == -1) then
if (ilaggr(i) > -(nr+1)) then if (ilaggr(i) > -(nr+1)) then
ilaggr(i) = abs(ilaggr(i)) ilaggr(i) = abs(ilaggr(i))
if (ilaggr(I)>naggr) then if (ilaggr(I)>naggr) then
write(0,*) 'loc_Aggregate: n > naggr 3? ',ilaggr(i),naggr write(0,*) 'loc_Aggregate: n > naggr 3? ',ilaggr(i),naggr
end if end if
ils(ilaggr(i)) = ils(ilaggr(i)) + 1 ils(ilaggr(i)) = ils(ilaggr(i)) + 1
! !
! This might happen if the pattern is non symmetric. ! This might happen if the pattern is non symmetric.
! Need a better handling. ! Need a better handling.
! !
recovery = .true. recovery = .true.
else else
write(0,*) 'Unrecoverable error !!',ilaggr(i), nr write(0,*) 'Unrecoverable error !!',ilaggr(i), nr
endif endif
else else
ilaggr(i) = ia ilaggr(i) = ia
if (ia>naggr) then if (ia>naggr) then
write(0,*) 'loc_Aggregate: n > naggr 4? ',ia,naggr write(0,*) 'loc_Aggregate: n > naggr 4? ',ia,naggr
end if end if
ils(ia) = ils(ia) + 1 ils(ia) = ils(ia) + 1
endif endif
end if end if
enddo enddo
if (recovery) then if (recovery) then
write(0,*) 'Had to recover from strange situation in loc_aggregate.' write(0,*) 'Had to recover from strange situation in loc_aggregate.'
write(0,*) 'Perhaps an unsymmetric pattern?' write(0,*) 'Perhaps an unsymmetric pattern?'
endif endif
if (debug) then if (debug) then
write(0,*) 'Phase 2: number of aggregates ',naggr write(0,*) 'Phase 2: number of aggregates ',naggr
write(0,*) 'Phase 2: nodes aggregated ',sum(ils) write(0,*) 'Phase 2: nodes aggregated ',sum(ils)
do i=1, naggr do i=1, naggr
write(*,*) 'Size of aggregate ',i,' :',count(ilaggr==i), ils(i) write(*,*) 'Size of aggregate ',i,' :',count(ilaggr==i), ils(i)
enddo enddo
write(*,*) maxval(ils(1:naggr)) write(*,*) maxval(ils(1:naggr))
write(*,*) 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops' write(*,*) 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops'
end if end if
!!$ write(0,*) 'desc_a loc_aggr 4 : ', desc_a%matrix_data(m_) !!$ write(0,*) 'desc_a loc_aggr 4 : ', desc_a%matrix_data(m_)
if (count(ilaggr<0) >0) then if (count(ilaggr<0) >0) then
write(0,*) 'Fatal error: some leftovers!!!' write(0,*) 'Fatal error: some leftovers!!!'
endif endif
deallocate(ils,neigh,stat=info) deallocate(ils,neigh,stat=info)
if (info/=0) then if (info/=0) then
info=4000 info=4000
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
if (nrow /= size(ilaggr)) then if (nrow /= size(ilaggr)) then
write(0,*) 'SOmething wrong ilaggr ',nrow,size(ilaggr) write(0,*) 'SOmething wrong ilaggr ',nrow,size(ilaggr)
endif endif
call psb_realloc(ncol,ilaggr,info) call psb_realloc(ncol,ilaggr,info)
if (info/=0) then if (info/=0) then
info=4010 info=4010
ch_err='psb_realloc' ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
allocate(nlaggr(np),stat=info) allocate(nlaggr(np),stat=info)
if (info/=0) then if (info/=0) then
info=4025 info=4025
call psb_errpush(info,name,i_err=(/np,0,0,0,0/),& call psb_errpush(info,name,i_err=(/np,0,0,0,0/),&
& a_err='integer') & a_err='integer')
goto 9999 goto 9999
end if end if
nlaggr(:) = 0
nlaggr(me+1) = naggr
call psb_sum(ictxt,nlaggr(1:np))
if (aggr_type == sym_dec_aggr_) then
call psb_sp_free(atmp,info)
end if
case default
nlaggr(:) = 0 write(0,*) 'Unimplemented aggregation algorithm ',aggr_type
nlaggr(me+1) = naggr info = -1
call psb_sum(ictxt,nlaggr(1:np)) call psb_errpush(30,name,i_err=(/1,aggr_type,0,0,0/))
goto 9999
end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -288,8 +323,8 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
9999 continue 9999 continue
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then if (err_act.eq.psb_act_abort_) then
call psb_error() call psb_error()
return return
end if end if
return return

Write Preview
Loading…
Cancel
Save