Merged internal docs and html headers.

stopcriterion
Salvatore Filippone 17 years ago
parent d8a4ceb5a5
commit 2272f944be

@ -1,6 +1,8 @@
Changelog. A lot less detailed than usual, at least for past
history.
2007/12/21: Merge version with prologues and internal docs.
2007/11/15: Created pargen example.
2007/11/14: Fix INTENT(IN) on X vector in preconditioner routines.

@ -41,7 +41,7 @@
!
! This routine builds a mapping from the row indices of the fine-level matrix
! to the row indices of the coarse-level matrix, according to a decoupled
! aggregation algorithm. This mapping will be used by mld_daggrmat_asb to
! aggregation algorithm. This mapping will be used by mld_aggrmat_asb to
! build the coarse-level matrix.
!
! The aggregation algorithm is a parallel version of that described in
@ -79,29 +79,30 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
implicit none
! Arguments
! Arguments
integer, intent(in) :: aggr_type
type(psb_dspmat_type), intent(in), target :: a
type(psb_desc_type), intent(in) :: desc_a
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
! Local variables
! Local variables
integer, allocatable :: ils(:), neigh(:)
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, parameter :: debug=.false.
integer :: debug_level, debug_unit
integer :: ictxt,np,me,err_act
integer :: nrow, ncol, n_ne
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() /= 0) return
info=0
name = 'mld_aggrmap_bld'
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!
! Note. At the time being we are ignoring aggr_type so
! that we only have decoupled aggregation. This might
@ -115,10 +116,9 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
select case (aggr_type)
case (mld_dec_aggr_,mld_sym_dec_aggr_)
nr = a%m
allocate(ilaggr(nr),neigh(nr),stat=info)
if(info.ne.0) then
if(info /= 0) then
info=4025
call psb_errpush(info,name,i_err=(/2*nr,0,0,0,0/),&
& a_err='integer')
@ -135,18 +135,23 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
& 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.)
if (info == 0) call psb_transp(atmp,atrans,fmt='COO')
if (info == 0) 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)
if (info == 0) call psb_sp_free(atrans,info)
if (info == 0) call psb_ipcoo2csr(atmp,info)
apnt => atmp
if (info/=0) then
info=4001
call psb_errpush(info,name,a_err='init apnt')
goto 9999
end if
end if
!
! Meaning of variables in the loops beloc
! -(nr+1) Untouched as yet
! Note: -(nr+1) Untouched as yet
! -i 1<=i<=nr Adjacent to aggregate i
! i 1<=i<=nr Belonging to aggregate i
@ -168,11 +173,10 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
naggr = naggr + 1
ilaggr(i) = naggr
call psb_neigh(apnt,i,neigh,n_ne,info,lev=one)
call psb_neigh(apnt,i,neigh,n_ne,info,lev=1)
if (info/=0) then
info=4010
ch_err='psb_neigh'
call psb_errpush(info,name,a_err=ch_err)
call psb_errpush(info,name,a_err='psb_neigh')
goto 9999
end if
do k=1, n_ne
@ -184,11 +188,10 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
!
! 2. Untouched neighbours of these nodes are marked <0.
!
call psb_neigh(apnt,i,neigh,n_ne,info,lev=two)
call psb_neigh(apnt,i,neigh,n_ne,info,lev=2)
if (info/=0) then
info=4010
ch_err='psb_neigh'
call psb_errpush(info,name,a_err=ch_err)
call psb_errpush(info,name,a_err='psb_neigh')
goto 9999
end if
@ -203,15 +206,17 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
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)
if (debug_level >= psb_debug_outer_) then
write(debug_unit,*) me,' ',trim(name),&
& ' 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
if(info /= 0) then
info=4025
call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),&
& a_err='integer')
@ -225,16 +230,20 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
n = ilaggr(i)
if (n>0) then
if (n>naggr) then
write(0,*) 'loc_Aggregate: n > naggr 1 ? ',n,naggr
info=4001
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 1 ?')
goto 9999
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)
if (debug_level >= psb_debug_outer_) then
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 1: number of aggregates ',naggr
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 1: nodes aggregated ',sum(ils)
end if
recovery=.false.
@ -247,11 +256,10 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
isz = nr+1
ia = -1
call psb_neigh(apnt,i,neigh,n_ne,info,lev=one)
call psb_neigh(apnt,i,neigh,n_ne,info,lev=1)
if (info/=0) then
info=4010
ch_err='psb_neigh'
call psb_errpush(info,name,a_err=ch_err)
call psb_errpush(info,name,a_err='psb_neigh')
goto 9999
end if
@ -261,7 +269,9 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
n = ilaggr(k)
if (n>0) then
if (n>naggr) then
write(0,*) 'loc_Aggregate: n > naggr 2? ',n,naggr
info=4001
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?')
goto 9999
end if
if (ils(n) < isz) then
@ -275,7 +285,9 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
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
info=4001
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 3 ?')
goto 9999
end if
ils(ilaggr(i)) = ils(ilaggr(i)) + 1
!
@ -284,35 +296,44 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
!
recovery = .true.
else
write(0,*) 'Unrecoverable error !!',ilaggr(i), nr
info=4001
call psb_errpush(info,name,a_err='Unrecoverable error !!')
goto 9999
endif
else
ilaggr(i) = ia
if (ia>naggr) then
write(0,*) 'loc_Aggregate: n > naggr 4? ',ia,naggr
info=4001
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 4? ')
goto 9999
end if
ils(ia) = ils(ia) + 1
endif
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)
if (debug_level >= psb_debug_outer_) then
if (recovery) then
write(debug_unit,*) me,' ',trim(name),&
& 'Had to recover from strange situation in loc_aggregate.'
write(debug_unit,*) me,' ',trim(name),&
& 'Perhaps an unsymmetric pattern?'
endif
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 2: number of aggregates ',naggr,sum(ils)
do i=1, naggr
write(*,*) 'Size of aggregate ',i,' :',count(ilaggr==i), ils(i)
write(debug_unit,*) me,' ',trim(name),&
& 'Size of aggregate ',i,' :',count(ilaggr==i), ils(i)
enddo
write(*,*) maxval(ils(1:naggr))
write(*,*) 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops'
write(debug_unit,*) me,' ',trim(name),&
& maxval(ils(1:naggr))
write(debug_unit,*) me,' ',trim(name),&
& 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops'
end if
!!$ write(0,*) 'desc_a loc_aggr 4 : ', desc_a%matrix_data(m_)
if (count(ilaggr<0) >0) then
write(0,*) 'Fatal error: some leftovers!!!'
info=4001
call psb_errpush(info,name,a_err='Fatal error: some leftovers')
goto 9999
endif
deallocate(ils,neigh,stat=info)
@ -322,9 +343,6 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
goto 9999
end if
if (nrow /= size(ilaggr)) then
write(0,*) 'SOmething wrong ilaggr ',nrow,size(ilaggr)
endif
call psb_realloc(ncol,ilaggr,info)
if (info/=0) then
info=4010
@ -351,7 +369,6 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
case default
write(0,*) 'Unimplemented aggregation algorithm ',aggr_type
info = -1
call psb_errpush(30,name,i_err=(/1,aggr_type,0,0,0/))
goto 9999

@ -107,7 +107,6 @@ subroutine mld_daggrmat_asb(a,desc_a,ac,desc_ac,p,info)
integer, intent(out) :: info
! Local variables
logical, parameter :: aggr_dump=.false.
integer ::ictxt,np,me, err_act, icomm
character(len=20) :: name
@ -125,24 +124,22 @@ subroutine mld_daggrmat_asb(a,desc_a,ac,desc_ac,p,info)
case (mld_no_smooth_)
call mld_aggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='raw_aggregate')
call psb_errpush(4010,name,a_err='mld_aggrmat_raw_asb')
goto 9999
end if
if (aggr_dump) call psb_csprt(90+me,ac,head='% Raw aggregate.')
case(mld_smooth_prol_,mld_biz_prol_)
if (aggr_dump) call psb_csprt(70+me,a,head='% Input matrix')
call mld_aggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
call mld_aggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='smooth_aggregate')
call psb_errpush(4010,name,a_err='mld_aggrmat_smth_asb')
goto 9999
end if
if (aggr_dump) call psb_csprt(90+me,ac,head='% Smooth aggregate.')
case default
call psb_errpush(4010,name,a_err=name)
call psb_errpush(4001,name,a_err='Invalid aggr kind')
goto 9999
end select

@ -95,7 +95,6 @@ subroutine mld_daggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
integer, intent(out) :: info
! Local variables
logical, parameter :: aggr_dump=.false.
integer ::ictxt,np,me, err_act, icomm
character(len=20) :: name
type(psb_dspmat_type) :: b
@ -202,7 +201,7 @@ subroutine mld_daggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
ac%k = ntaggr
ac%infoa(psb_nnz_) = nzac
ac%fida='COO'
ac%descra='G'
ac%descra='GUN'
call psb_spcnv(ac,info,afmt='coo',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='sp_free')
@ -212,19 +211,10 @@ subroutine mld_daggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
else if (p%iprcparm(mld_coarse_mat_) == mld_distr_mat_) then
call psb_cdall(ictxt,desc_ac,info,nl=naggr)
if (info == 0) call psb_cdasb(desc_ac,info)
if (info == 0) call psb_sp_clone(b,ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdall')
goto 9999
end if
call psb_cdasb(desc_ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdasb')
goto 9999
end if
call psb_sp_clone(b,ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spclone')
call psb_errpush(4001,name,a_err='Build ac, desc_ac')
goto 9999
end if
call psb_sp_free(b,info)
@ -234,8 +224,9 @@ subroutine mld_daggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
end if
else
write(0,*) 'Unknown p%iprcparm(coarse_mat) in aggregate_sp',p%iprcparm(mld_coarse_mat_)
info = 4001
call psb_errpush(4001,name,a_err='invalid mld_coarse_mat_')
goto 9999
end if
deallocate(nzbr,idisp)

@ -117,14 +117,15 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
integer, pointer :: nzbr(:), idisp(:)
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k
logical, parameter :: aggr_dump=.false.
integer ::ictxt,np,me, err_act, icomm
character(len=20) :: name
type(psb_dspmat_type), pointer :: am1,am2
type(psb_dspmat_type) :: am3,am4
logical :: ml_global_nmb
logical, parameter :: test_dump=.false.,debug=.false.
integer :: nz
integer, allocatable :: ia(:), ja(:)
real(kind(1.d0)), allocatable :: val(:)
integer :: debug_level, debug_unit
integer, parameter :: ncmax=16
real(kind(1.d0)) :: omega, anorm, tmp, dg
@ -132,10 +133,12 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
if(psb_get_errstatus().ne.0) return
info=0
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = psb_cd_get_context(desc_a)
icomm = psb_cd_get_mpic(desc_a)
ictxt=psb_cd_get_context(desc_a)
ictxt = psb_cd_get_context(desc_a)
call psb_info(ictxt, me, np)
@ -164,38 +167,27 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
goto 9999
end if
naggrm1 = sum(p%nlaggr(1:me))
naggrp1 = sum(p%nlaggr(1:me+1))
ml_global_nmb = ( (p%iprcparm(mld_aggr_kind_) == mld_smooth_prol_).or.&
& ( (p%iprcparm(mld_aggr_kind_) == mld_biz_prol_).and.&
& (p%iprcparm(mld_coarse_mat_) == mld_repl_mat_)) )
if (ml_global_nmb) then
p%mlia(1:nrow) = p%mlia(1:nrow) + naggrm1
call psb_halo(p%mlia,desc_a,info)
if(info /= 0) then
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_halo')
goto 9999
end if
end if
if (aggr_dump) then
open(30+me)
write(30+me,*) '% Aggregation map'
do i=1,ncol
write(30+me,*) i,p%mlia(i)
end do
close(30+me)
end if
! naggr: number of local aggregates
! nrow: local rows.
!
allocate(p%dorig(nrow),stat=info)
if (info /= 0) then
info=4025
call psb_errpush(info,name,i_err=(/nrow,0,0,0,0/),&
@ -218,13 +210,6 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
end if
end do
! where (p%dorig /= dzero)
! p%dorig = done / p%dorig
! elsewhere
! p%dorig = done
! end where
! 1. Allocate Ptilde in sparse matrix form
am4%fida='COO'
am4%m=ncol
@ -235,7 +220,8 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
am4%k=naggr
call psb_sp_all(ncol,naggr,am4,ncol,info)
endif
if(info /= 0) then
if (info /= 0) then
call psb_errpush(4010,name,a_err='spall')
goto 9999
end if
@ -258,16 +244,14 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
call psb_spcnv(am4,info,afmt='csr',dupl=psb_dupl_add_)
if(info /= 0) then
if (info==0) call psb_spcnv(a,am3,info,afmt='csr')
if (info /= 0) then
call psb_errpush(4010,name,a_err='spcnv')
goto 9999
end if
call psb_sp_clone(a,am3,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spclone')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Initial copies done.'
!
! WARNING: the cycles below assume that AM3 does have
@ -275,7 +259,7 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
! Should we switch to something safer?
!
call psb_sp_scal(am3,p%dorig,info)
if(info /= 0) goto 9999
if (info /= 0) goto 9999
if (p%iprcparm(mld_aggr_eig_) == mld_max_norm_) then
@ -284,25 +268,33 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
!
! This only works with CSR.
!
anorm = dzero
dg = done
do i=1,am3%m
tmp = dzero
do j=am3%ia2(i),am3%ia2(i+1)-1
if (am3%ia1(j) <= am3%m) then
tmp = tmp + dabs(am3%aspk(j))
endif
if (am3%ia1(j) == i ) then
dg = dabs(am3%aspk(j))
end if
end do
anorm = max(anorm,tmp/dg)
enddo
call psb_amx(ictxt,anorm)
if (toupper(am3%fida)=='CSR') then
anorm = dzero
dg = done
do i=1,am3%m
tmp = dzero
do j=am3%ia2(i),am3%ia2(i+1)-1
if (am3%ia1(j) <= am3%m) then
tmp = tmp + dabs(am3%aspk(j))
endif
if (am3%ia1(j) == i ) then
dg = dabs(am3%aspk(j))
end if
end do
anorm = max(anorm,tmp/dg)
enddo
call psb_amx(ictxt,anorm)
else
info = 4001
endif
else
anorm = psb_spnrmi(am3,desc_a,info)
endif
if (info /= 0) then
call psb_errpush(4001,name,a_err='Invalid AM3 storage format')
goto 9999
end if
omega = 4.d0/(3.d0*anorm)
p%dprcparm(mld_aggr_damp_) = omega
@ -311,8 +303,9 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
omega = p%dprcparm(mld_aggr_damp_)
else if (p%iprcparm(mld_aggr_eig_) /= mld_user_choice_) then
write(0,*) me,'Error: invalid choice for OMEGA in blaggrmat?? ',&
& p%iprcparm(mld_aggr_eig_)
info = 4001
call psb_errpush(info,name,a_err='invalid mld_aggr_eig_')
goto 9999
end if
@ -326,41 +319,14 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
end if
end do
end do
else if (toupper(am3%fida)=='COO') then
do j=1,am3%infoa(psb_nnz_)
if (am3%ia1(j) /= am3%ia2(j)) then
am3%aspk(j) = - omega*am3%aspk(j)
else
am3%aspk(j) = done - omega*am3%aspk(j)
endif
end do
call psb_spcnv(am3,info,afmt='csr',dupl=psb_dupl_add_)
if (info /=0) then
call psb_errpush(4010,name,a_err='spcnv am3')
goto 9999
end if
else
write(0,*) 'Missing implementation of I sum'
call psb_errpush(4010,name)
call psb_errpush(4001,name,a_err='Invalid AM3 storage format')
goto 9999
end if
if (test_dump) then
open(30+me)
write(30+me,*) 'OMEGA: ',omega
do i=1,size(p%dorig)
write(30+me,*) p%dorig(i)
end do
close(30+me)
end if
if (test_dump) call &
& psb_csprt(20+me,am4,head='% Operator Ptilde.',ivr=desc_a%loc_to_glob)
if (test_dump) call psb_csprt(40+me,am3,head='% (I-wDA)',ivr=desc_a%loc_to_glob,&
& ivc=desc_a%loc_to_glob)
if (debug) write(0,*) me,'Done gather, going for SYMBMM 1'
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done gather, going for SYMBMM 1'
!
! Symbmm90 does the allocation for its result.
!
@ -376,7 +342,9 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
call psb_numbmm(am3,am4,am1)
if (debug) write(0,*) me,'Done NUMBMM 1'
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
call psb_sp_free(am4,info)
if(info /= 0) then
@ -391,35 +359,15 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
!
call psb_sphalo(am1,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_sphalo')
goto 9999
end if
call psb_rwextd(ncol,am1,info,b=am4)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_rwextd')
goto 9999
end if
call psb_sp_free(am4,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
goto 9999
end if
if (info == 0) call psb_rwextd(ncol,am1,info,b=am4)
if (info == 0) call psb_sp_free(am4,info)
else
call psb_rwextd(ncol,am1,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='rwextd')
goto 9999
end if
endif
if (test_dump) &
& call psb_csprt(60+me,am1,head='% (I-wDA)Pt',ivr=desc_a%loc_to_glob)
if(info /= 0) then
call psb_errpush(4001,name,a_err='Halo of am1')
goto 9999
end if
call psb_symbmm(a,am1,am3,info)
if(info /= 0) then
@ -428,7 +376,9 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
end if
call psb_numbmm(a,am1,am3)
if (debug) write(0,*) me,'Done NUMBMM 2'
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2'
if (p%iprcparm(mld_aggr_kind_) == mld_smooth_prol_) then
call psb_transp(am1,am2,fmt='COO')
@ -446,7 +396,6 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
am2%ia2(i) = am2%ia2(k)
end if
end do
am2%infoa(psb_nnz_) = i
call psb_spcnv(am2,info,afmt='csr',dupl=psb_dupl_add_)
if (info /=0) then
@ -456,63 +405,38 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
else
call psb_transp(am1,am2)
endif
if (debug) write(0,*) me,'starting sphalo/ rwxtd'
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd'
if (p%iprcparm(mld_aggr_kind_) == mld_smooth_prol_) then
if (p%iprcparm(mld_aggr_kind_) == mld_smooth_prol_) then
! am2 = ((i-wDA)Ptilde)^T
call psb_sphalo(am3,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_sphalo')
goto 9999
end if
call psb_rwextd(ncol,am3,info,b=am4)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_rwextd')
goto 9999
end if
call psb_sp_free(am4,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
goto 9999
end if
if (info == 0) call psb_rwextd(ncol,am3,info,b=am4)
if (info == 0) call psb_sp_free(am4,info)
else if (p%iprcparm(mld_aggr_kind_) == mld_biz_prol_) then
call psb_rwextd(ncol,am3,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_rwextd')
goto 9999
end if
endif
if (debug) write(0,*) me,'starting symbmm 3'
call psb_symbmm(am2,am3,b,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='symbmm 3')
call psb_errpush(4001,name,a_err='Extend am3')
goto 9999
end if
if (debug) write(0,*) me,'starting numbmm 3'
call psb_numbmm(am2,am3,b)
if (debug) write(0,*) me,'Done NUMBMM 3'
!!$ if (aggr_dump) call csprt(50+me,am1,head='% Operator PTrans.')
call psb_sp_free(am3,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
goto 9999
end if
call psb_spcnv(b,info,afmt='coo',dupl=psb_dupl_add_)
if (info /=0) then
call psb_errpush(4010,name,a_err='spcnv b')
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting symbmm 3'
call psb_symbmm(am2,am3,b,info)
if (info == 0) call psb_numbmm(am2,am3,b)
if (info == 0) call psb_sp_free(am3,info)
if (info == 0) call psb_spcnv(b,info,afmt='coo',dupl=psb_dupl_add_)
if (info /= 0) then
call psb_errpush(4001,name,a_err='Build b = am2 x am3')
goto 9999
end if
if (test_dump) call psb_csprt(80+me,b,head='% Smoothed aggregate AC.')
select case(p%iprcparm(mld_aggr_kind_))
@ -523,55 +447,30 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
case(mld_distr_mat_)
call psb_sp_clone(b,ac,info)
if(info /= 0) goto 9999
nzac = ac%infoa(psb_nnz_)
nzl = ac%infoa(psb_nnz_)
call psb_cdall(ictxt,desc_ac,info,nl=p%nlaggr(me+1))
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdall')
goto 9999
end if
call psb_cdins(nzl,ac%ia1,ac%ia2,desc_ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdins')
goto 9999
end if
if (debug) write(0,*) me,'Created aux descr. distr.'
call psb_cdasb(desc_ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdasb')
goto 9999
end if
if (debug) write(0,*) me,'Asmbld aux descr. distr.'
call psb_glob_to_loc(ac%ia1(1:nzl),desc_ac,info,iact='I')
if(info /= 0) then
call psb_errpush(4010,name,a_err='psglob_to_loc')
goto 9999
end if
call psb_glob_to_loc(ac%ia2(1:nzl),desc_ac,info,iact='I')
if(info /= 0) then
call psb_errpush(4010,name,a_err='psglob_to_loc')
if (info == 0) call psb_cdall(ictxt,desc_ac,info,nl=p%nlaggr(me+1))
if (info == 0) call psb_cdins(nzl,ac%ia1,ac%ia2,desc_ac,info)
if (info == 0) call psb_cdasb(desc_ac,info)
if (info == 0) call psb_glob_to_loc(ac%ia1(1:nzl),desc_ac,info,iact='I')
if (info == 0) call psb_glob_to_loc(ac%ia2(1:nzl),desc_ac,info,iact='I')
if (info /= 0) then
call psb_errpush(4001,name,a_err='Creating desc_ac and converting ac')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
ac%m=desc_ac%matrix_data(psb_n_row_)
ac%k=desc_ac%matrix_data(psb_n_col_)
ac%fida='COO'
ac%descra='G'
ac%descra='GUN'
call psb_sp_free(b,info)
if (info == 0) deallocate(nzbr,idisp,stat=info)
if(info /= 0) then
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
goto 9999
end if
@ -579,7 +478,6 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
if (np>1) then
nzl = psb_sp_get_nnzeros(am1)
call psb_glob_to_loc(am1%ia1(1:nzl),desc_ac,info,'I')
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_glob_to_loc')
goto 9999
@ -589,44 +487,31 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
if (np>1) then
call psb_spcnv(am2,info,afmt='coo',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spcnv')
goto 9999
end if
nzl = am2%infoa(psb_nnz_)
call psb_glob_to_loc(am2%ia1(1:nzl),desc_ac,info,'I')
if (info == 0) call psb_glob_to_loc(am2%ia1(1:nzl),desc_ac,info,'I')
if (info == 0) call psb_spcnv(am2,info,afmt='csr',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_glob_to_loc')
goto 9999
end if
call psb_spcnv(am2,info,afmt='csr',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spcnv')
call psb_errpush(4001,name,a_err='Converting am2 to local')
goto 9999
end if
end if
am2%m=desc_ac%matrix_data(psb_n_col_)
if (debug) write(0,*) me,'Done ac '
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done ac '
case(mld_repl_mat_)
!
!
call psb_cdall(ictxt,desc_ac,info,mg=ntaggr,repl=.true.)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdall')
goto 9999
end if
nzbr(:) = 0
nzbr(me+1) = b%infoa(psb_nnz_)
call psb_sum(ictxt,nzbr(1:np))
nzac = sum(nzbr)
call psb_sp_all(ntaggr,ntaggr,ac,nzac,info)
if(info /= 0) goto 9999
if (info == 0) call psb_sp_all(ntaggr,ntaggr,ac,nzac,info)
if (info /= 0) goto 9999
do ip=1,np
idisp(ip) = sum(nzbr(1:ip-1))
@ -635,30 +520,36 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
call mpi_allgatherv(b%aspk,ndx,mpi_double_precision,ac%aspk,nzbr,idisp,&
& mpi_double_precision,icomm,info)
call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,&
if (info == 0) call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,&
& mpi_integer,icomm,info)
call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,&
if (info == 0) call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,&
& mpi_integer,icomm,info)
if(info /= 0) goto 9999
if (info /= 0) then
call psb_errpush(4001,name,a_err=' from mpi_allgatherv')
goto 9999
end if
ac%m = ntaggr
ac%k = ntaggr
ac%infoa(psb_nnz_) = nzac
ac%fida='COO'
ac%descra='G'
ac%descra='GUN'
call psb_spcnv(ac,info,afmt='coo',dupl=psb_dupl_add_)
if(info /= 0) goto 9999
call psb_sp_free(b,info)
if(info /= 0) goto 9999
if (me==0) then
if (test_dump) call psb_csprt(80+me,ac,head='% Smoothed aggregate AC.')
endif
deallocate(nzbr,idisp)
deallocate(nzbr,idisp,stat=info)
if (info /= 0) then
info = 4000
call psb_errpush(info,name)
goto 9999
end if
case default
write(0,*) 'Inconsistent input in smooth_new_aggregate'
info = 4001
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
@ -669,25 +560,11 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
case(mld_distr_mat_)
call psb_sp_clone(b,ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spclone')
goto 9999
end if
call psb_cdall(ictxt,desc_ac,info,nl=naggr)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdall')
goto 9999
end if
call psb_cdasb(desc_ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdasb')
goto 9999
end if
call psb_sp_free(b,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='sp_free')
if (info == 0) call psb_cdall(ictxt,desc_ac,info,nl=naggr)
if (info == 0) call psb_cdasb(desc_ac,info)
if (info == 0) call psb_sp_free(b,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Build desc_ac, ac')
goto 9999
end if
@ -718,13 +595,12 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
call mpi_allgatherv(b%aspk,ndx,mpi_double_precision,ac%aspk,nzbr,idisp,&
& mpi_double_precision,icomm,info)
call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,&
if (info == 0) call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,&
& mpi_integer,icomm,info)
call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,&
if (info == 0) call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,&
& mpi_integer,icomm,info)
if(info /= 0) then
info=-1
call psb_errpush(info,name)
if (info /= 0) then
call psb_errpush(4001,name,a_err=' from mpi_allgatherv')
goto 9999
end if
@ -733,7 +609,7 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
ac%k = ntaggr
ac%infoa(psb_nnz_) = nzac
ac%fida='COO'
ac%descra='G'
ac%descra='GUN'
call psb_spcnv(ac,info,afmt='coo',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spcnv')
@ -745,8 +621,23 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
goto 9999
end if
case default
info = 4001
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
deallocate(nzbr,idisp)
deallocate(nzbr,idisp,stat=info)
if (info /= 0) then
info = 4000
call psb_errpush(info,name)
goto 9999
end if
case default
info = 4001
call psb_errpush(info,name,a_err='invalid mld_smooth_prol_')
goto 9999
end select
@ -756,7 +647,9 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
goto 9999
end if
if (debug) write(0,*) me,'Done smooth_aggregate '
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done smooth_aggregate '
call psb_erractionrestore(err_act)
return

@ -78,7 +78,7 @@ subroutine mld_dasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
Implicit None
! Arguments
! Arguments
integer, intent(in) :: ptype,novr
Type(psb_dspmat_type), Intent(in) :: a
Type(psb_dspmat_type), Intent(inout) :: blk
@ -88,21 +88,24 @@ subroutine mld_dasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
Character, Intent(in) :: upd
character(len=5), optional :: outfmt
! Local variables
real(kind(1.d0)) :: t1,t2,t3
! Local variables
integer icomm
Integer :: np,me,nnzero,&
& ictxt, n_col,int_err(5),&
& tot_recv, n_row,nhalo, nrow_a,err_act
Logical,Parameter :: debug=.false.
integer :: debug_level, debug_unit
character(len=20) :: name, ch_err
name='mld_dasmat_bld'
if(psb_get_errstatus().ne.0) return
if(psb_get_errstatus() /= 0) return
info=0
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
If(debug) Write(0,*)'IN DASMATBLD ', upd
If (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' start ', upd
ictxt = psb_cd_get_context(desc_data)
icomm = psb_cd_get_mpic(desc_data)
@ -121,7 +124,6 @@ subroutine mld_dasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
! Block-Jacobi preconditioner. Copy the descriptor, just in case
! we want to renumber the rows and columns of the matrix.
!
If(debug) Write(0,*)' asmatbld calling allocate '
call psb_sp_all(0,0,blk,1,info)
if(info /= 0) then
info=4010
@ -131,10 +133,11 @@ subroutine mld_dasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
end if
blk%fida = 'COO'
blk%infoa(psb_nnz_) = 0
If(debug) Write(0,*)' asmatbld done spallocate'
If (upd == 'F') Then
call psb_cdcpy(desc_data,desc_p,info)
If(debug) Write(0,*)' asmatbld done cdcpy'
If(debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' done cdcpy'
if(info /= 0) then
info=4010
ch_err='psb_cdcpy'
@ -144,12 +147,9 @@ subroutine mld_dasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
endif
case(mld_as_)
!
! Additive Schwarz
!
if (novr < 0) then
info=3
int_err(1)=novr
@ -161,7 +161,9 @@ subroutine mld_dasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
!
! Actually, this is just block Jacobi
!
If(debug) Write(0,*)' asmatbld calling allocate novr=0'
If(debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' calling allocate novr=0'
call psb_sp_all(0,0,blk,1,info)
if(info /= 0) then
info=4010
@ -171,25 +173,28 @@ subroutine mld_dasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
end if
blk%fida = 'COO'
blk%infoa(psb_nnz_) = 0
if (debug) write(0,*) 'Calling desccpy'
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Calling desccpy'
if (upd == 'F') then
call psb_cdcpy(desc_data,desc_p,info)
If(debug) Write(0,*)' asmatbld done cdcpy'
If(debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' done cdcpy'
if(info /= 0) then
info=4010
ch_err='psb_cdcpy'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (debug) write(0,*) 'Early return from asmatbld: P>=3 N_OVR=0'
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Early return: P>=3 N_OVR=0'
endif
return
endif
If(debug)Write(0,*)'BEGIN dasmatbld',me,upd,novr
t1 = psb_wtime()
If (upd == 'F') Then
!
! Build the auxiliary descriptor desc_p%matrix_data(psb_n_row_).
@ -200,7 +205,7 @@ subroutine mld_dasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
! a descriptor for an extended stencil in a PDE solver.
!
call psb_cdbldext(a,desc_data,novr,desc_p,info,extype=psb_ovt_asov_)
if(info /= 0) then
if (info /= 0) then
info=4010
ch_err='psb_cdbldext'
call psb_errpush(info,name,a_err=ch_err)
@ -208,7 +213,9 @@ subroutine mld_dasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
end if
Endif
if(debug) write(0,*) me,' From cdbldext _:',desc_p%matrix_data(psb_n_row_),&
if(debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' From cdbldext _:',desc_p%matrix_data(psb_n_row_),&
& desc_p%matrix_data(psb_n_col_)
!
@ -216,30 +223,35 @@ subroutine mld_dasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
!
n_row = desc_p%matrix_data(psb_n_row_)
t2 = psb_wtime()
if (debug) write(0,*) 'Before sphalo ',blk%fida,blk%m,psb_nnz_,blk%infoa(psb_nnz_)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Before sphalo ',blk%fida,blk%m,psb_nnz_,blk%infoa(psb_nnz_)
Call psb_sphalo(a,desc_p,blk,info,&
& outfmt=outfmt,data=psb_comm_ext_,rowscale=.true.)
if(info /= 0) then
if (info /= 0) then
info=4010
ch_err='psb_sphalo'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (debug) write(0,*) 'After psb_sphalo ',&
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'After psb_sphalo ',&
& blk%fida,blk%m,psb_nnz_,blk%infoa(psb_nnz_)
case default
if(info /= 0) then
info=4000
info=4001
ch_err='Invalid ptype'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
End select
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),'Done'
call psb_erractionrestore(err_act)
return

@ -98,7 +98,6 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
real(kind(1.d0)), pointer :: ww(:), aux(:), tx(:),ty(:)
character ::diagl, diagu
integer :: ictxt,np,me,isz, err_act
logical, parameter :: debug=.false., debugprt=.false.
character(len=20) :: name, ch_err
name='mld_dbaseprec_aply'
@ -164,7 +163,7 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
!
call mld_bjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='mld_bjac_aply'
goto 9999
@ -180,7 +179,7 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
! shortcut: this fixes performance for RAS(0) == BJA
!
call mld_bjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='psb_bjacaply'
goto 9999
@ -230,9 +229,6 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
endif
if (debugprt) write(0,*)' vdiag: ',prec%d(:)
if (debug) write(0,*) 'Bi-CGSTAB with Additive Schwarz prec'
tx(1:nrow_d) = x(1:nrow_d)
tx(nrow_d+1:isz) = dzero
@ -247,8 +243,8 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
goto 9999
end if
else if (prec%iprcparm(mld_sub_restr_) /= psb_none_) then
write(0,*) 'Problem in PREC_APLY: Unknown value for restriction ',&
&prec%iprcparm(mld_sub_restr_)
call psb_errpush(4001,name,a_err='Invalid mld_sub_restr_')
goto 9999
end if
!
@ -270,7 +266,7 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
! preconditioner). The resulting vector is ty.
!
call mld_bjac_aply(done,prec,tx,dzero,ty,prec%desc_data,trans,aux,info)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='mld_bjac_aply'
goto 9999
@ -281,7 +277,7 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
!
if (prec%iprcparm(mld_sub_ren_)>0) then
call psb_gelp('n',prec%invperm,ty,info)
if(info /=0) then
if(info /= 0) then
info=4010
ch_err='psb_gelp'
goto 9999
@ -309,8 +305,8 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
end if
case default
write(0,*) 'Problem in PREC_APLY: Unknown value for prolongation ',&
& prec%iprcparm(mld_sub_prol_)
call psb_errpush(4001,name,a_err='Invalid mld_sub_prol_')
goto 9999
end select
!
@ -330,9 +326,8 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
end if
case default
write(0,*) 'Invalid PRE%PREC ',prec%iprcparm(mld_prec_type_),':',&
& mld_min_prec_,mld_noprec_,mld_diag_,mld_bjac_,mld_as_
call psb_errpush(4001,name,a_err='Invalid mld_prec_type_')
goto 9999
end select
call psb_erractionrestore(err_act)
@ -341,7 +336,7 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
9999 continue
call psb_errpush(info,name,i_err=int_err,a_err=ch_err)
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if

@ -77,34 +77,33 @@ subroutine mld_dbaseprc_bld(a,desc_a,p,info,upd)
! Local variables
Integer :: err, n_row, n_col,ictxt, me,np,mglob, err_act
integer :: int_err(5)
character :: iupd
logical, parameter :: debug=.false.
integer,parameter :: iroot=0,iout=60,ilout=40
integer :: debug_level, debug_unit
character(len=20) :: name, ch_err
if(psb_get_errstatus().ne.0) return
if (psb_get_errstatus() /= 0) return
name = 'mld_dbaseprc_bld'
info=0
err=0
call psb_erractionsave(err_act)
name = 'mld_dbaseprc_bld'
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
if (debug) write(0,*) 'Entering baseprc_bld'
info = 0
int_err(1) = 0
ictxt = psb_cd_get_context(desc_a)
n_row = psb_cd_get_local_rows(desc_a)
n_col = psb_cd_get_local_cols(desc_a)
mglob = psb_cd_get_global_rows(desc_a)
if (debug) write(0,*) 'Preconditioner Blacs_gridinfo'
call psb_info(ictxt, me, np)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' start'
if (present(upd)) then
if (debug) write(0,*) 'UPD ', upd
if ((UPD.eq.'F').or.(UPD.eq.'T')) then
IUPD=UPD
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),'UPD ', upd
if ((toupper(UPD) == 'F').or.(toupper(UPD) == 'T')) then
IUPD=toupper(UPD)
else
IUPD='F'
endif
@ -140,7 +139,9 @@ subroutine mld_dbaseprc_bld(a,desc_a,p,info,upd)
! Diagonal preconditioner
call mld_diag_bld(a,desc_a,p,iupd,info)
if(debug) write(0,*)me,': out of mld_diag_bld'
if(debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ': out of mld_diag_bld'
if(info /= 0) then
info=4010
ch_err='mld_diag_bld'
@ -168,8 +169,9 @@ subroutine mld_dbaseprc_bld(a,desc_a,p,info,upd)
p%iprcparm(mld_smooth_sweeps_) = 1
end if
if (debug) write(0,*)me, ': Calling mld_bjac_bld'
if (debug) call psb_barrier(ictxt)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ': Calling mld_bjac_bld'
! Build the local part of the base preconditioner
call mld_bjac_bld(a,desc_a,p,iupd,info)
@ -180,7 +182,7 @@ subroutine mld_dbaseprc_bld(a,desc_a,p,info,upd)
end if
case default
info=4010
info=4001
ch_err='Unknown mld_prec_type_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
@ -190,6 +192,8 @@ subroutine mld_dbaseprc_bld(a,desc_a,p,info,upd)
p%base_a => a
p%base_desc => desc_a
p%iprcparm(mld_prec_status_) = mld_prec_built_
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),': Done'
call psb_erractionrestore(err_act)
return

@ -150,8 +150,7 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
! Local variables
integer :: n_row,n_col
real(kind(1.d0)), pointer :: ww(:), aux(:), tx(:),ty(:)
integer :: ictxt,np,me,i, err_act, int_err(5)
logical,parameter :: debug=.false., debugprt=.false.
integer :: ictxt,np,me,i, err_act
character(len=20) :: name
interface
@ -206,10 +205,6 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
end if
endif
if (debug) then
write(0,*) me,' mld_bjac_APLY: ',prec%iprcparm(mld_sub_solve_),prec%iprcparm(mld_smooth_sweeps_)
end if
if (prec%iprcparm(mld_smooth_sweeps_) == 1) then
!
! TASKS 1, 3 and 4
@ -228,19 +223,17 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
call psb_spsm(done,prec%av(mld_l_pr_),x,dzero,ww,desc_data,info,&
& trans='N',unit='L',diag=prec%d,choice=psb_none_,work=aux)
if(info /=0) goto 9999
call psb_spsm(alpha,prec%av(mld_u_pr_),ww,beta,y,desc_data,info,&
if (info == 0) call psb_spsm(alpha,prec%av(mld_u_pr_),ww,beta,y,desc_data,info,&
& trans='N',unit='U',choice=psb_none_, work=aux)
if(info /=0) goto 9999
case('T','C')
call psb_spsm(done,prec%av(mld_u_pr_),x,dzero,ww,desc_data,info,&
& trans=trans,unit='L',diag=prec%d,choice=psb_none_,work=aux)
if(info /=0) goto 9999
call psb_spsm(alpha,prec%av(mld_l_pr_),ww,beta,y,desc_data,info,&
if (info == 0) call psb_spsm(alpha,prec%av(mld_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans,unit='U',choice=psb_none_,work=aux)
if(info /=0) goto 9999
case default
call psb_errpush(4001,name,a_err='Invalid TRANS in ILU subsolve')
goto 9999
end select
case(mld_slu_)
@ -258,10 +251,12 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
call mld_dslu_solve(0,n_row,1,ww,n_row,prec%iprcparm(mld_slu_ptr_),info)
case('T','C')
call mld_dslu_solve(1,n_row,1,ww,n_row,prec%iprcparm(mld_slu_ptr_),info)
case default
call psb_errpush(4001,name,a_err='Invalid TRANS in SLU subsolve')
goto 9999
end select
if(info /=0) goto 9999
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
if (info ==0) call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
case(mld_sludist_)
!
@ -276,10 +271,12 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
call mld_dsludist_solve(0,n_row,1,ww,n_row,prec%iprcparm(mld_slud_ptr_),info)
case('T','C')
call mld_dsludist_solve(1,n_row,1,ww,n_row,prec%iprcparm(mld_slud_ptr_),info)
case default
call psb_errpush(4001,name,a_err='Invalid TRANS in SLUDist subsolve')
goto 9999
end select
if(info /=0) goto 9999
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
if (info == 0) call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
case (mld_umf_)
!
@ -295,16 +292,22 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
call mld_dumf_solve(0,n_row,ww,x,n_row,prec%iprcparm(mld_umf_numptr_),info)
case('T','C')
call mld_dumf_solve(1,n_row,ww,x,n_row,prec%iprcparm(mld_umf_numptr_),info)
case default
call psb_errpush(4001,name,a_err='Invalid TRANS in UMF subsolve')
goto 9999
end select
if(info /=0) goto 9999
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
if (info == 0) call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
case default
write(0,*) 'Unknown factorization type in mld_bjac_aply',prec%iprcparm(mld_sub_solve_)
call psb_errpush(4001,name,a_err='Invalid mld_sub_solve_')
goto 9999
end select
if (debugprt) write(0,*)' Y: ',y(:)
if (info /= 0) then
call psb_errpush(4001,name,a_err='Error in subsolve Jacobi Sweeps = 1')
goto 9999
endif
else if (prec%iprcparm(mld_smooth_sweeps_) > 1) then
@ -346,23 +349,23 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
ty(1:n_row) = x(1:n_row)
call psb_spmm(-done,prec%av(mld_ap_nd_),tx,done,ty,&
& prec%desc_data,info,work=aux)
if(info /=0) goto 9999
if (info /=0) exit
call psb_spsm(done,prec%av(mld_l_pr_),ty,dzero,ww,&
& prec%desc_data,info,&
& trans='N',unit='L',diag=prec%d,choice=psb_none_,work=aux)
if(info /=0) goto 9999
if (info /=0) exit
call psb_spsm(done,prec%av(mld_u_pr_),ww,dzero,tx,&
& prec%desc_data,info,&
& trans='N',unit='U',choice=psb_none_,work=aux)
if(info /=0) goto 9999
if (info /=0) exit
end do
case(mld_sludist_)
!
! Wrong choice: SuperLU_DIST
!
write(0,*) 'No sense in having SuperLU_DIST with multiple Jacobi sweeps'
info=4010
info = 4001
call psb_errpush(4001,name,a_err='Invalid SuperLU_DIST with Jacobi sweeps >1')
goto 9999
case(mld_slu_)
@ -379,10 +382,10 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
ty(1:n_row) = x(1:n_row)
call psb_spmm(-done,prec%av(mld_ap_nd_),tx,done,ty,&
& prec%desc_data,info,work=aux)
if(info /=0) goto 9999
if(info /=0) exit
call mld_dslu_solve(0,n_row,1,ty,n_row,prec%iprcparm(mld_slu_ptr_),info)
if(info /=0) goto 9999
if(info /=0) exit
tx(1:n_row) = ty(1:n_row)
end do
@ -400,23 +403,34 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
ty(1:n_row) = x(1:n_row)
call psb_spmm(-done,prec%av(mld_ap_nd_),tx,done,ty,&
& prec%desc_data,info,work=aux)
if(info /=0) goto 9999
if(info /=0) exit
call mld_dumf_solve(0,n_row,ww,ty,n_row,&
& prec%iprcparm(mld_umf_numptr_),info)
if(info /=0) goto 9999
if(info /=0) exit
tx(1:n_row) = ww(1:n_row)
end do
case default
call psb_errpush(4001,name,a_err='Invalid mld_sub_solve_')
goto 9999
end select
if (info /= 0) then
info=4001
call psb_errpush(info,name,a_err='subsolve with Jacobi sweeps > 1')
goto 9999
end if
!
! Put the result into the output vector Y.
!
call psb_geaxpby(alpha,tx,beta,y,desc_data,info)
deallocate(tx,ty)
deallocate(tx,ty,stat=info)
if (info /= 0) then
info=4001
call psb_errpush(info,name,a_err='final cleanup with Jacobi sweeps > 1')
goto 9999
end if
else
@ -440,7 +454,6 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
return
9999 continue
call psb_errpush(info,name,i_err=int_err)
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()

@ -59,19 +59,14 @@
! 2. setup of block-Jacobi sweeps to compute an approximate solution of a
! linear system
! A*Y = X,
!
! distributed among the processes (allowed only at the coarsest level);
!
! 3. LU factorization of a linear system
!
! A*Y = X,
!
! distributed among the processes (allowed only at the coarsest level);
!
! 4. LU or incomplete LU factorization of a linear system
!
! A*Y = X,
!
! replicated on the processes (allowed only at the coarsest level).
!
! The following factorizations are available:
@ -116,7 +111,7 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
integer :: int_err(5)
character :: trans, unitd
type(psb_dspmat_type) :: blck, atmp
logical, parameter :: debugprt=.false., debug=.false., aggr_dump=.false.
integer :: debug_level, debug_unit
integer :: err_act, n_row, nrow_a,n_col
integer :: ictxt,np,me
character(len=20) :: name
@ -126,8 +121,9 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
info=0
name='mld_dbjac_bld'
call psb_erractionsave(err_act)
ictxt=psb_cd_get_context(desc_a)
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)
m = a%m
@ -152,9 +148,9 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
call psb_nullify_sp(atmp)
if(debug) write(0,*)me,': calling mld_asmat_bld',&
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),': Start',&
& p%iprcparm(mld_prec_type_),p%iprcparm(mld_n_ovr_)
if (debug) call psb_barrier(ictxt)
!
! Build the communication descriptor for the Additive Schwarz
@ -166,22 +162,14 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
call mld_asmat_bld(p%iprcparm(mld_prec_type_),p%iprcparm(mld_n_ovr_),a,&
& blck,desc_a,upd,p%desc_data,info,outfmt=csrfmt)
if (debugprt) then
open(60+me)
call psb_csprt(60+me,a,head='% A')
close(60+me)
open(70+me)
call psb_csprt(70+me,blck,head='% BLCK')
close(70+me)
endif
if(info/=0) then
if (info/=0) then
call psb_errpush(4010,name,a_err='mld_asmat_bld')
goto 9999
end if
if (debug) write(0,*)me,': out of mld_asmat_bld'
if (debug) call psb_barrier(ictxt)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ': out of mld_asmat_bld'
!
! Treat separately the case the local matrix has to be reordered
@ -200,7 +188,6 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
! matrix is stored into atmp, using the COO format.
!
call mld_sp_renum(a,desc_a,blck,p,atmp,info)
if (info/=0) then
call psb_errpush(4010,name,a_err='mld_sp_renum')
goto 9999
@ -212,10 +199,10 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
!
call psb_sp_clip(atmp,p%av(mld_ap_nd_),info,&
& jmin=atmp%m+1,rscale=.false.,cscale=.false.)
call psb_spcnv(p%av(mld_ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv csr 1')
if (info == 0) call psb_spcnv(p%av(mld_ap_nd_),info,&
& afmt='csr',dupl=psb_dupl_add_)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv')
goto 9999
end if
@ -232,14 +219,9 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
end if
if (debugprt) then
call psb_barrier(ictxt)
open(40+me)
call psb_csprt(40+me,atmp,head='% Local matrix')
close(40+me)
endif
if (debug) write(0,*) me,' Factoring rows ',&
&atmp%m,a%m,blck%m,atmp%ia2(atmp%m+1)-1
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' Factoring rows ',&
& atmp%m,a%m,blck%m,atmp%ia2(atmp%m+1)-1
!
! Compute a factorization of the diagonal block of the local matrix,
@ -248,90 +230,56 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
select case(p%iprcparm(mld_sub_solve_))
case(mld_ilu_n_,mld_milu_n_,mld_ilu_t_)
!
! ILU(k)/MILU(k)/ILU(k,t) factorization.
!
!
! ILU(k)/MILU(k)/ILU(k,t) factorization.
!
call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv csr 2')
goto 9999
end if
call mld_ilu_bld(atmp,p%desc_data,p,upd,info)
if (info == 0) call mld_ilu_bld(atmp,p%desc_data,p,upd,info)
if (info/=0) then
call psb_errpush(4010,name,a_err='mld_ilu_bld')
goto 9999
end if
if (debugprt) then
open(80+me)
call psb_csprt(80+me,p%av(mld_l_pr_),head='% Local L factor')
write(80+me,*) '% Diagonal: ',p%av(mld_l_pr_)%m
do i=1,p%av(mld_l_pr_)%m
write(80+me,*) i,i,p%d(i)
enddo
call psb_csprt(80+me,p%av(mld_u_pr_),head='% Local U factor')
close(80+me)
endif
case(mld_slu_)
!
! LU factorization through the SuperLU package.
!
!
! LU factorization through the SuperLU package.
!
call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_)
if (info == 0) call mld_slu_bld(atmp,p%desc_data,p,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv csr 3')
goto 9999
end if
call mld_slu_bld(atmp,p%desc_data,p,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='mld_slu_bld')
goto 9999
end if
case(mld_umf_)
!
! LU factorization through the UMFPACK package.
!
!
! LU factorization through the UMFPACK package.
!
call psb_spcnv(atmp,info,afmt='csc',dupl=psb_dupl_add_)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv csc')
goto 9999
end if
call mld_umf_bld(atmp,p%desc_data,p,info)
if(debug) write(0,*)me,': Done mld_umf_bld ',info
if (info == 0) call mld_umf_bld(atmp,p%desc_data,p,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='mld_umf_bld')
goto 9999
end if
case(mld_f_none_)
!
! Error: no factorization required.
!
info=4010
!
! Error: no factorization required.
!
info=4001
call psb_errpush(info,name,a_err='Inconsistent prec mld_f_none_')
goto 9999
case default
info=4010
info=4001
call psb_errpush(info,name,a_err='Unknown mld_sub_solve_')
goto 9999
end select
call psb_sp_free(atmp,info)
if(info/=0) then
if (info/=0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
goto 9999
end if
@ -347,11 +295,10 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
!
select case(p%iprcparm(mld_sub_solve_))
case(mld_ilu_n_,mld_milu_n_,mld_ilu_t_)
!
! ILU(k)/MILU(k)/ILU(k,t) factorization.
!
!
! ILU(k)/MILU(k)/ILU(k,t) factorization.
!
!
! In case of multiple block-Jacobi sweeps, clip into p%av(ap_nd_)
@ -367,13 +314,13 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
! given that the output from CLIP is in COO.
call psb_sp_clip(a,p%av(mld_ap_nd_),info,&
& jmin=nrow_a+1,rscale=.false.,cscale=.false.)
call psb_sp_clip(blck,atmp,info,&
if (info == 0) call psb_sp_clip(blck,atmp,info,&
& jmin=nrow_a+1,rscale=.false.,cscale=.false.)
call psb_rwextd(n_row,p%av(mld_ap_nd_),info,b=atmp)
call psb_spcnv(p%av(mld_ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv csr 4')
if (info == 0) call psb_rwextd(n_row,p%av(mld_ap_nd_),info,b=atmp)
if (info == 0) call psb_spcnv(p%av(mld_ap_nd_),info,&
& afmt='csr',dupl=psb_dupl_add_)
if (info /= 0) then
call psb_errpush(4010,name,a_err='clip & psb_spcnv csr 4')
goto 9999
end if
@ -390,45 +337,23 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
end if
call psb_sp_free(atmp,info)
end if
!
! Compute the incomplete LU factorization.
!
call mld_ilu_bld(a,desc_a,p,upd,info,blck=blck)
if(info/=0) then
if (info/=0) then
call psb_errpush(4010,name,a_err='mld_ilu_bld')
goto 9999
end if
if (debugprt) then
open(80+me)
call psb_csprt(80+me,p%av(mld_l_pr_),head='% Local L factor')
write(80+me,*) '% Diagonal: ',p%av(mld_l_pr_)%m
do i=1,p%av(mld_l_pr_)%m
write(80+me,*) i,i,p%d(i)
enddo
call psb_csprt(80+me,p%av(mld_u_pr_),head='% Local U factor')
close(80+me)
endif
case(mld_slu_)
!
! LU factorization through the SuperLU package.
!
call psb_spcnv(a,atmp,info,afmt='coo')
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv')
goto 9999
end if
!
! LU factorization through the SuperLU package.
!
n_row = psb_cd_get_local_rows(p%desc_data)
n_col = psb_cd_get_local_cols(p%desc_data)
call psb_rwextd(n_row,atmp,info,b=blck)
call psb_spcnv(a,atmp,info,afmt='coo')
if (info == 0) call psb_rwextd(n_row,atmp,info,b=blck)
!
! In case of multiple block-Jacobi sweeps, clip into p%av(ap_nd_)
@ -437,11 +362,12 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
!
if (p%iprcparm(mld_smooth_sweeps_) > 1) then
call psb_sp_clip(atmp,p%av(mld_ap_nd_),info,&
if (info == 0) call psb_sp_clip(atmp,p%av(mld_ap_nd_),info,&
& jmin=atmp%m+1,rscale=.false.,cscale=.false.)
call psb_spcnv(p%av(mld_ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
if(info /= 0) then
if (info == 0) call psb_spcnv(p%av(mld_ap_nd_),info,&
& afmt='csr',dupl=psb_dupl_add_)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv csr 6')
goto 9999
end if
@ -458,19 +384,18 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
p%iprcparm(mld_smooth_sweeps_) = 1
end if
endif
!
! Compute the LU factorization.
!
if (info == 0) call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_)
if (info == 0) call mld_slu_bld(atmp,p%desc_data,p,info)
if(info /= 0) then
if (info /= 0) then
call psb_errpush(4010,name,a_err='mld_slu_bld')
goto 9999
end if
call psb_sp_free(atmp,info)
if(info/=0) then
if (info/=0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
goto 9999
end if
@ -481,23 +406,15 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
! when the matrix is distributed among the processes.
! NOTE: Should have NO overlap here!!!!
!
call psb_spcnv(a,atmp,info,afmt='csr')
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv')
goto 9999
end if
n_row = psb_cd_get_local_rows(p%desc_data)
n_col = psb_cd_get_local_cols(p%desc_data)
if (info == 0) call mld_sludist_bld(atmp,p%desc_data,p,info)
if(info /= 0) then
if (info /= 0) then
call psb_errpush(4010,name,a_err='mld_slu_bld')
goto 9999
end if
call psb_sp_free(atmp,info)
if(info/=0) then
if (info/=0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
goto 9999
end if
@ -527,13 +444,12 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
call psb_sp_clip(atmp,p%av(mld_ap_nd_),info,&
& jmin=atmp%m+1,rscale=.false.,cscale=.false.)
call psb_spcnv(p%av(mld_ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
if(info /= 0) then
if (info == 0) call psb_spcnv(p%av(mld_ap_nd_),info,&
& afmt='csr',dupl=psb_dupl_add_)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv csr 8')
goto 9999
end if
k = psb_sp_get_nnzeros(p%av(mld_ap_nd_))
call psb_sum(ictxt,k)
@ -551,19 +467,17 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
! Compute the LU factorization.
!
if (info == 0) call psb_ipcoo2csc(atmp,info,clshr=.true.)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_ipcoo2csc')
goto 9999
end if
call mld_umf_bld(atmp,p%desc_data,p,info)
if(debug) write(0,*)me,': Done mld_umf_bld ',info
if (info == 0) call mld_umf_bld(atmp,p%desc_data,p,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ': Done mld_umf_bld ',info
if (info /= 0) then
call psb_errpush(4010,name,a_err='mld_umf_bld')
goto 9999
end if
call psb_sp_free(atmp,info)
if(info/=0) then
if (info/=0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
goto 9999
end if
@ -573,31 +487,30 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
!
! Error: no factorization required.
!
info=4010
info=4001
call psb_errpush(info,name,a_err='Inconsistent prec mld_f_none_')
goto 9999
case default
info=4010
info=4001
call psb_errpush(info,name,a_err='Unknown mld_sub_solve_')
goto 9999
end select
case default
info=4010
info=4001
call psb_errpush(info,name,a_err='Invalid renum_')
goto 9999
end select
call psb_sp_free(blck,info)
if(info/=0) then
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
goto 9999
end if
if (debug) write(0,*) me,'End of ilu_bld'
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),'End '
call psb_erractionrestore(err_act)

@ -72,30 +72,25 @@ subroutine mld_ddiag_bld(a,desc_a,p,upd,info)
! Local variables
Integer :: err, n_row, n_col,I,j,k,ictxt,&
& me,np,mglob,lw, err_act
integer :: int_err(5)
logical, parameter :: debug=.false.
integer,parameter :: iroot=0,iout=60,ilout=40
integer :: debug_level, debug_unit
character(len=20) :: name, ch_err
if(psb_get_errstatus().ne.0) return
info=0
err=0
call psb_erractionsave(err_act)
name = 'mld_ddiag_bld'
if (debug) write(0,*) 'Entering diagsc_bld'
info = 0
int_err(1) = 0
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
name = 'mld_ddiag_bld'
info = 0
ictxt = psb_cd_get_context(desc_a)
n_row = psb_cd_get_local_rows(desc_a)
n_col = psb_cd_get_local_cols(desc_a)
mglob = psb_cd_get_global_rows(desc_a)
if (debug) write(0,*) 'Preconditioner Blacs_gridinfo'
call psb_info(ictxt, me, np)
if (debug) write(0,*) 'Precond: Diagonal'
if (debug_level >= psb_debug_outer_)&
& write(debug_unit,*) me,' ',trim(name),' Enter'
call psb_realloc(n_col,p%d,info)
if (info /= 0) then
@ -122,7 +117,6 @@ subroutine mld_ddiag_bld(a,desc_a,p,upd,info)
goto 9999
end if
if (debug) write(ilout+me,*) 'VDIAG ',n_row
!
! The i-th diagonal entry of the preconditioner is set to one if the
! corresponding entry a_ii of the sparse matrix A is zero; otherwise
@ -134,8 +128,6 @@ subroutine mld_ddiag_bld(a,desc_a,p,upd,info)
else
p%d(i) = done/p%d(i)
endif
if (debug) write(ilout+me,*) i,desc_a%loc_to_glob(i), p%d(i)
end do
if (a%pl(1) /= 0) then
@ -151,8 +143,8 @@ subroutine mld_ddiag_bld(a,desc_a,p,upd,info)
end if
endif
if (debug) write(*,*) 'Preconditioner DIAG computed OK'
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),'Done'
call psb_erractionrestore(err_act)
return

@ -106,18 +106,20 @@ subroutine mld_dilu_bld(a,desc_a,p,upd,info,blck)
! Local Variables
integer :: i, nztota, err_act, n_row, nrow_a
character :: trans, unitd
logical, parameter :: debugprt=.false., debug=.false., aggr_dump=.false.
integer :: ictxt,np,me
integer :: debug_level, debug_unit
integer :: ictxt,np,me
character(len=20) :: name, ch_err
if(psb_get_errstatus().ne.0) return
info=0
name='mld_dilu_bld'
call psb_erractionsave(err_act)
ictxt=psb_cd_get_context(desc_a)
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'
trans = 'N'
unitd = 'U'
@ -145,15 +147,15 @@ subroutine mld_dilu_bld(a,desc_a,p,upd,info,blck)
goto 9999
end if
endif
!!$ call psb_csprt(50+me,a,head='% (A)')
nrow_a = psb_cd_get_local_rows(desc_a)
nztota = psb_sp_get_nnzeros(a)
if (present(blck)) then
nztota = nztota + psb_sp_get_nnzeros(blck)
end if
if (debug) write(0,*)me,': out get_nnzeros',nztota,a%m,a%k
if (debug) call psb_barrier(ictxt)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ': out get_nnzeros',nztota,a%m,a%k
n_row = p%desc_data%matrix_data(psb_n_row_)
p%av(mld_l_pr_)%m = n_row
@ -253,22 +255,6 @@ subroutine mld_dilu_bld(a,desc_a,p,upd,info,blck)
end select
if (debugprt) then
!
! Print out the factors on file.
!
open(80+me)
call psb_csprt(80+me,p%av(mld_l_pr_),head='% Local L factor')
write(80+me,*) '% Diagonal: ',p%av(mld_l_pr_)%m
do i=1,p%av(mld_l_pr_)%m
write(80+me,*) i,i,p%d(i)
enddo
call psb_csprt(80+me,p%av(mld_u_pr_),head='% Local U factor')
close(80+me)
endif
if (psb_sp_getifld(psb_upd_,p%av(mld_u_pr_),info) /= psb_upd_perm_) then
call psb_sp_trim(p%av(mld_u_pr_),info)
endif
@ -277,7 +263,8 @@ subroutine mld_dilu_bld(a,desc_a,p,upd,info,blck)
call psb_sp_trim(p%av(mld_l_pr_),info)
endif
if (debug) write(0,*) me,'End of ilu_bld'
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),'End'
call psb_erractionrestore(err_act)
return

@ -116,7 +116,6 @@ subroutine mld_dilu_fct(ialg,a,l,u,d,info,blck)
integer :: l1, l2,m,err_act
type(psb_dspmat_type), pointer :: blck_
character(len=20) :: name, ch_err
logical, parameter :: debug=.false.
name='mld_dilu_fct'
info = 0
@ -290,7 +289,6 @@ contains
integer :: i,j,k,l,low1,low2,kk,jj,ll, ktrw,err_act
real(kind(1.d0)) :: dia,temp
integer, parameter :: nrb=16
logical,parameter :: debug=.false.
type(psb_dspmat_type) :: trw
integer :: int_err(5)
character(len=20) :: name, ch_err
@ -302,7 +300,7 @@ contains
call psb_nullify_sp(trw)
trw%m=0
trw%k=0
if(debug) write(0,*)'LUINT Allocating TRW'
call psb_sp_all(trw,1,info)
if(info.ne.0) then
info=4010
@ -310,20 +308,18 @@ contains
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if(debug) write(0,*)'LUINT Done Allocating TRW'
lia2(1) = 1
uia2(1) = 1
l1 = 0
l2 = 0
m = ma+mb
if(debug) write(0,*)'In DCSRLU Begin cycle',m,ma,mb
!
! Cycle over the matrix rows
!
do i = 1, m
if(debug) write(0,*)'LUINT: Loop index ',i,ma
d(i) = dzero
if (i <= ma) then
@ -448,7 +444,6 @@ contains
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if(debug) write(0,*)'Leaving ilu_fct'
call psb_erractionrestore(err_act)
return
@ -465,7 +460,7 @@ contains
!
! Subroutine: ilu_copyin
! Version: real
! Note: internal subroutine of mld_dilu_fct.
! Note: internal subroutine of mld_dilu_fct
!
! This routine copies a row of a sparse matrix A, stored in the psb_dspmat_type
! data structure a, into the arrays laspk and uaspk and into the scalar variable

@ -112,14 +112,11 @@ subroutine mld_diluk_fct(fill_in,ialg,a,l,u,d,info,blck)
type(psb_dspmat_type), pointer :: blck_
character(len=20) :: name, ch_err
logical, parameter :: debug=.false.
name='mld_diluk_fct'
info = 0
call psb_erractionsave(err_act)
if (debug) write(0,*) 'mld_diluk_fct: start'
!
! Point to / allocate memory for the incomplete factorization
!
@ -145,7 +142,6 @@ subroutine mld_diluk_fct(fill_in,ialg,a,l,u,d,info,blck)
!
! Compute the ILU(k) or the MILU(k) factorization, depending on ialg
!
if (debug) write(0,*) 'mld_diluk_fct: calling fctint'
call mld_diluk_fctint(fill_in,ialg,m,a%m,a,blck_%m,blck_,&
& d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info)
if (info /= 0) then
@ -289,7 +285,6 @@ contains
integer, allocatable :: uplevs(:), rowlevs(:),idxs(:)
real(kind(1.d0)), allocatable :: row(:)
type(psb_int_heap) :: heap
logical,parameter :: debug=.false.
type(psb_dspmat_type) :: trw
character(len=20), parameter :: name='mld_diluk_fctint'
character(len=20) :: ch_err
@ -303,7 +298,6 @@ contains
!
! Allocate a temporary buffer for the iluk_copyin function
!
if (debug) write(0,*)'LUINT Allocating TRW'
call psb_sp_all(0,0,trw,1,info)
if (info==0) call psb_ensure_size(m+1,lia2,info)
if (info==0) call psb_ensure_size(m+1,uia2,info)
@ -313,15 +307,12 @@ contains
call psb_errpush(info,name,a_err='psb_sp_all')
goto 9999
end if
if (debug) write(0,*)'LUINT Done Allocating TRW'
l1=0
l2=0
lia2(1) = 1
uia2(1) = 1
if (debug) write(0,*)'In DCSRLU Begin cycle',m,ma,mb
!
! Allocate memory to hold the entries of a row and the corresponding
! fill levels
@ -342,8 +333,6 @@ contains
!
do i = 1, m
if (debug.and.(mod(i,500)==1)) write(0,*)'LUINT: Loop index ',i,ma
!
! At each iteration of the loop we keep in a heap the column indices
! affected by the factorization. The heap is initialized and filled
@ -366,8 +355,6 @@ contains
call iluk_copyin(i-ma,mb,b,1,m,row,rowlevs,heap,ktrw,trw)
endif
if (debug) write(0,*)'LUINT: input Copy done'
! Do an elimination step on the current row. It turns out we only
! need to keep track of fill levels for the upper triangle, hence we
! do not have a lowlevs variable.
@ -398,7 +385,6 @@ contains
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (debug) write(0,*)'Leaving ilu_fct'
call psb_erractionrestore(err_act)
return

@ -111,14 +111,11 @@ subroutine mld_dilut_fct(fill_in,thres,ialg,a,l,u,d,info,blck)
type(psb_dspmat_type), pointer :: blck_
character(len=20) :: name, ch_err
logical, parameter :: debug=.false.
name='mld_dilut_fct'
info = 0
call psb_erractionsave(err_act)
if (debug) write(0,*) 'mld_dilut_fct: start'
!
! Point to / allocate memory for the incomplete factorization
!
@ -144,7 +141,6 @@ subroutine mld_dilut_fct(fill_in,thres,ialg,a,l,u,d,info,blck)
!
! Compute the ILU(k,t) factorization
!
if (debug) write(0,*) 'mld_dilut_fct: calling fctint'
call mld_dilut_fctint(fill_in,thres,ialg,m,a%m,a,blck_%m,blck_,&
& d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info)
if (info /= 0) then
@ -288,7 +284,6 @@ contains
integer, allocatable :: idxs(:)
real(kind(1.d0)), allocatable :: row(:)
type(psb_int_heap) :: heap
logical,parameter :: debug=.false.
type(psb_dspmat_type) :: trw
character(len=20), parameter :: name='mld_dilut_fctint'
character(len=20) :: ch_err
@ -302,7 +297,6 @@ contains
!
! Allocate a temporary buffer for the ilut_copyin function
!
if (debug) write(0,*)'LUINT Allocating TRW'
call psb_sp_all(0,0,trw,1,info)
if (info==0) call psb_ensure_size(m+1,lia2,info)
if (info==0) call psb_ensure_size(m+1,uia2,info)
@ -312,15 +306,12 @@ contains
call psb_errpush(info,name,a_err='psb_sp_all')
goto 9999
end if
if (debug) write(0,*)'LUINT Done Allocating TRW'
l1=0
l2=0
lia2(1) = 1
uia2(1) = 1
if (debug) write(0,*)'In DCSRLU Begin cycle',m,ma,mb
!
! Allocate memory to hold the entries of a row
!
@ -338,7 +329,6 @@ contains
!
do i = 1, m
if (debug) write(0,*)'LUINT: Loop index ',i
!
! At each iteration of the loop we keep in a heap the column indices
! affected by the factorization. The heap is initialized and filled
@ -354,7 +344,6 @@ contains
call ilut_copyin(i-ma,mb,b,i,1,m,nlw,nup,jmaxup,nrmi,row,heap,ktrw,trw)
endif
if (debug) write(0,*)'LUINT: input Copy done'
!
! Do an elimination step on current row
!
@ -384,7 +373,6 @@ contains
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (debug) write(0,*)'Leaving ilu_fct'
call psb_erractionrestore(err_act)
return
@ -675,7 +663,6 @@ contains
! Local Variables
integer :: k,j,jj,info, lastk
real(kind(1.d0)) :: rwk
logical, parameter :: debug=.false.
call psb_ensure_size(200,idxs,info)
@ -748,12 +735,6 @@ contains
end do
if (debug) then
write(0,*) 'At end of factint: ',i,nidx
write(0,*) idxs(1:nidx)
write(0,*) row(idxs(1:nidx))
end if
end subroutine ilut_fact
!
@ -870,7 +851,6 @@ contains
character(len=20), parameter :: name='mld_dilut_fctint'
character(len=20) :: ch_err
logical :: fndmaxup
logical, parameter :: debug=.false.
if (psb_get_errstatus() /= 0) return
info=0
@ -909,10 +889,6 @@ contains
if (idxs(idxp) >= i) exit
widx = idxs(idxp)
witem = row(widx)
if (debug) then
write(0,*) 'Lower: Deciding on drop of item ',witem,widx,thres,nrmi,thres*nrmi
end if
!
! Dropping rule based on the 2-norm
!
@ -1032,11 +1008,6 @@ contains
cycle
end if
witem = row(widx)
if (debug) then
write(0,*) 'Upper: Deciding on drop of item ',witem,widx,&
& jmaxup,thres,nrmi,thres*nrmi
end if
!
! Dropping rule based on the 2-norm. But keep the jmaxup-th entry anyway.
!
@ -1051,14 +1022,6 @@ contains
end do
if (debug) then
write(0,*) 'Row ',i,' copyout: after first round at upper:',nz,jmaxup
write(0,*) xwid(1:nz)
write(0,*) xw(1:nz)
write(0,*) 'Dumping heap'
call psb_dump_heap(0,heap,info)
end if
!
! Now we have to take out the first nup-fill_in entries. But make sure
! we include entry jmaxup.
@ -1093,11 +1056,6 @@ contains
! Now we put things back into ascending column order
!
call psb_msort(xwid(1:nz),indx(1:nz),dir=psb_sort_up_)
if (debug) then
write(0,*) 'Row ',i,' copyout: after sort at upper:',nz,jmaxup
write(0,*) xwid(1:nz)
write(0,*) xw(indx(1:nz))
end if
!
! Copy out the upper part of the row

@ -64,11 +64,11 @@
! level 1 is the finest level and A(1) is the matrix A.
!
! For a general description of (parallel) multilevel preconditioners see
! 1. B.F. Smith, P.E. Bjorstad & W.D. Gropp,
! - B.F. Smith, P.E. Bjorstad & W.D. Gropp,
! Domain decomposition: parallel multilevel methods for elliptic partial
! differential equations,
! Cambridge University Press, 1996.
! 2. K. Stuben,
! - K. Stuben,
! Algebraic Multigrid (AMG): An Introduction with Applications,
! GMD Report N. 70, 1999.
!
@ -113,11 +113,11 @@
! The real parameters defining the base preconditioner
! K(ilev).
! baseprecv(ilev)%perm - integer, dimension(:), allocatable.
! The row and column permutations applied to the local
! The row and column permutations applied to the local
! part of A(ilev) (defined only if baseprecv(ilev)%
! iprcparm(mld_sub_ren_)>0).
! baseprecv(ilev)%invperm - integer, dimension(:), allocatable.
! The inverse of the permutation stored in
! The inverse of the permutation stored in
! baseprecv(ilev)%perm.
! baseprecv(ilev)%mlia - integer, dimension(:), allocatable.
! The aggregation map (ilev-1) --> (ilev).
@ -182,7 +182,7 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
! Local variables
integer :: n_row,n_col
integer :: ictxt,np,me,i, nr2l,nc2l,err_act
logical, parameter :: debug=.false., debugprt=.false.
integer :: debug_level, debug_unit
integer :: ismth, nlev, ilev, icm
character(len=20) :: name
@ -194,12 +194,15 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
name='mld_dmlprec_aply'
info = 0
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = psb_cd_get_context(desc_data)
call psb_info(ictxt, me, np)
if (debug) write(0,*) me,'Entry to mlprec_aply ',&
& size(baseprecv)
if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(baseprecv)
nlev = size(baseprecv)
allocate(mlprec_wrk(nlev),stat=info)
@ -215,7 +218,7 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
! No preconditioning, should not really get here
!
call psb_errpush(4010,name,a_err='mld_no_ml_ in mlprc_aply?')
call psb_errpush(4001,name,a_err='mld_no_ml_ in mlprc_aply?')
goto 9999
@ -260,7 +263,10 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
call mld_baseprec_aply(alpha,baseprecv(1),x,beta,y,&
& baseprecv(1)%base_desc,trans,work,info)
if(info /=0) goto 9999
if (info /=0) then
call psb_errpush(4010,name,a_err='baseprec_aply')
goto 9999
end if
allocate(mlprec_wrk(1)%x2l(size(x)),mlprec_wrk(1)%y2l(size(y)), stat=info)
if (info /= 0) then
@ -308,11 +314,8 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
call psb_halo(mlprec_wrk(ilev-1)%x2l,baseprecv(ilev-1)%base_desc,&
& info,work=work)
if(info /=0) goto 9999
call psb_csmm(done,baseprecv(ilev)%av(mld_sm_pr_t_),mlprec_wrk(ilev-1)%x2l,&
& dzero,mlprec_wrk(ilev)%x2l,info)
if(info /=0) goto 9999
if (info == 0) call psb_csmm(done,baseprecv(ilev)%av(mld_sm_pr_t_),&
& mlprec_wrk(ilev-1)%x2l,dzero,mlprec_wrk(ilev)%x2l,info)
else
!
! Apply the raw aggregation map transpose (take a shortcut)
@ -324,11 +327,19 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
end do
end if
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction')
goto 9999
end if
if (icm == mld_repl_mat_) then
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nr2l))
else if (icm /= mld_distr_mat_) then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(mld_coarse_mat_) ',icm
info = 4013
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_',&
& i_Err=(/icm,0,0,0,0/))
goto 9999
endif
!
@ -361,8 +372,6 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
call psb_csmm(done,baseprecv(ilev)%av(mld_sm_pr_),mlprec_wrk(ilev)%y2l,&
& done,mlprec_wrk(ilev-1)%y2l,info)
if(info /=0) goto 9999
else
!
! Apply the raw aggregation map (take a shortcut)
@ -373,6 +382,10 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
enddo
end if
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during prolognation')
goto 9999
end if
end do
!
@ -381,8 +394,10 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
! Compute the output vector Y
!
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,done,y,baseprecv(1)%base_desc,info)
if(info /=0) goto 9999
if (info /= 0) then
call psb_errpush(4001,name,a_err='Error on final update')
goto 9999
end if
case(mld_mult_ml_)
@ -432,8 +447,9 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
! Copy the input vector X
!
if (debug) write(0,*) me, 'mlprec_aply desc_data',&
& allocated(desc_data%matrix_data)
if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' desc_data status',allocated(desc_data%matrix_data)
n_col = psb_cd_get_local_cols(desc_data)
nc2l = psb_cd_get_local_cols(baseprecv(1)%desc_data)
@ -442,7 +458,7 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
& mlprec_wrk(1)%tx(nc2l), stat=info)
mlprec_wrk(1)%x2l(:) = dzero
mlprec_wrk(1)%y2l(:) = dzero
mlprec_wrk(1)%tx(:) = dzero
mlprec_wrk(1)%tx(:) = dzero
call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%tx,&
& baseprecv(1)%base_desc,info)
@ -463,7 +479,9 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
ismth = baseprecv(ilev)%iprcparm(mld_aggr_kind_)
icm = baseprecv(ilev)%iprcparm(mld_coarse_mat_)
if (debug) write(0,*) me, 'mlprec_aply starting up sweep ',&
if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name), &
& ' starting up sweep ',&
& ilev,allocated(baseprecv(ilev)%iprcparm),n_row,n_col,&
& nc2l, nr2l,ismth
@ -479,21 +497,19 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
mlprec_wrk(ilev)%x2l(:) = dzero
mlprec_wrk(ilev)%y2l(:) = dzero
mlprec_wrk(ilev)%tx(:) = dzero
mlprec_wrk(ilev)%tx(:) = dzero
if (ismth /= mld_no_smooth_) then
!
! Apply the smoothed prolongator transpose
!
if (debug) write(0,*) me, 'mlprec_aply halo in up sweep ', ilev
if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name), ' up sweep ', ilev
call psb_halo(mlprec_wrk(ilev-1)%x2l,&
& baseprecv(ilev-1)%base_desc,info,work=work)
if(info /=0) goto 9999
if (debug) write(0,*) me, 'mlprec_aply csmm in up sweep ', ilev
call psb_csmm(done,baseprecv(ilev)%av(mld_sm_pr_t_),mlprec_wrk(ilev-1)%x2l, &
& dzero,mlprec_wrk(ilev)%x2l,info)
if(info /=0) goto 9999
if (info == 0) call psb_csmm(done,baseprecv(ilev)%av(mld_sm_pr_t_),&
& mlprec_wrk(ilev-1)%x2l,dzero,mlprec_wrk(ilev)%x2l,info)
else
!
! Apply the raw aggregation map transpose (take a shortcut)
@ -505,20 +521,18 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
end do
end if
if (debug) write(0,*) me, 'mlprec_aply possible sum in up sweep ', &
& ilev,icm,associated(baseprecv(ilev)%base_desc),mld_repl_mat_
if (debug) write(0,*) me, 'mlprec_aply geaxpby in up sweep X', &
& ilev,associated(baseprecv(ilev)%base_desc),&
& baseprecv(ilev)%base_desc%matrix_data(psb_n_row_),&
& baseprecv(ilev)%base_desc%matrix_data(psb_n_col_),&
& size(mlprec_wrk(ilev)%tx),size(mlprec_wrk(ilev)%x2l)
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction')
goto 9999
end if
if (icm == mld_repl_mat_) Then
if (debug) write(0,*) 'Entering psb_sum ',nr2l
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nr2l))
else if (icm /= mld_distr_mat_) Then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(mld_coarse_mat_) ', icm
info = 4013
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_',&
& i_Err=(/icm,0,0,0,0/))
goto 9999
endif
!
@ -526,9 +540,14 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
call psb_geaxpby(done,mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%tx,&
& baseprecv(ilev)%base_desc,info)
if(info /=0) goto 9999
if (debug) write(0,*) me, 'mlprec_aply done up sweep ',&
& ilev
if (info /= 0) then
call psb_errpush(4001,name,a_err='Error in update')
goto 9999
end if
if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' done up sweep ', ilev
enddo
@ -539,10 +558,13 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
call mld_baseprec_aply(done,baseprecv(nlev),mlprec_wrk(nlev)%x2l, &
& dzero, mlprec_wrk(nlev)%y2l,baseprecv(nlev)%desc_data,'N',work,info)
if (info /=0) then
call psb_errpush(4010,name,a_err='baseprec_aply')
goto 9999
end if
if(info /=0) goto 9999
if (debug) write(0,*) me, 'mlprec_aply done prc_apl ',&
& nlev
if (debug_level >= psb_debug_inner_) write(debug_unit,*) &
& me,' ',trim(name), ' done baseprec_aply ', nlev
!
! STEP 4
@ -551,7 +573,10 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
do ilev=nlev-1, 1, -1
if (debug) write(0,*) me, 'mlprec_aply starting down sweep',ilev
if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' starting down sweep',ilev
ismth = baseprecv(ilev+1)%iprcparm(mld_aggr_kind_)
n_row = psb_cd_get_local_rows(baseprecv(ilev)%base_desc)
@ -562,10 +587,8 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
if (ismth == mld_smooth_prol_) &
& call psb_halo(mlprec_wrk(ilev+1)%y2l,baseprecv(ilev+1)%desc_data,&
& info,work=work)
call psb_csmm(done,baseprecv(ilev+1)%av(mld_sm_pr_),mlprec_wrk(ilev+1)%y2l,&
& dzero,mlprec_wrk(ilev)%y2l,info)
if(info /=0) goto 9999
if (info == 0) call psb_csmm(done,baseprecv(ilev+1)%av(mld_sm_pr_),&
& mlprec_wrk(ilev+1)%y2l, dzero,mlprec_wrk(ilev)%y2l,info)
else
!
! Apply the raw aggregation map (take a shortcut)
@ -575,7 +598,10 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
mlprec_wrk(ilev)%y2l(i) = mlprec_wrk(ilev)%y2l(i) + &
& mlprec_wrk(ilev+1)%y2l(baseprecv(ilev+1)%mlia(i))
enddo
end if
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during prolongation')
goto 9999
end if
!
@ -584,16 +610,19 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
call psb_spmm(-done,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
& done,mlprec_wrk(ilev)%tx,baseprecv(ilev)%base_desc,info,work=work)
if(info /=0) goto 9999
!
! Apply the base preconditioner
!
call mld_baseprec_aply(done,baseprecv(ilev),mlprec_wrk(ilev)%tx,&
if (info == 0) call mld_baseprec_aply(done,baseprecv(ilev),mlprec_wrk(ilev)%tx,&
& done,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%base_desc, trans, work,info)
if (info /=0) then
call psb_errpush(4001,name,a_err=' spmm/baseprec_aply')
goto 9999
end if
if(info /=0) goto 9999
if (debug) write(0,*) me, 'mlprec_aply done down sweep',ilev
if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' done down sweep',ilev
enddo
!
@ -603,8 +632,10 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,baseprecv(1)%base_desc,info)
if(info /=0) goto 9999
if (info /=0) then
call psb_errpush(4001,name,a_err=' Final update')
goto 9999
end if
case(mld_pre_smooth_)
@ -675,9 +706,10 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
& dzero,mlprec_wrk(1)%y2l,&
& baseprecv(1)%base_desc,&
& trans,work,info)
if(info /=0) goto 9999
if (info /=0) then
call psb_errpush(4010,name,a_err=' baseprec_aply')
goto 9999
end if
!
! STEP 3
@ -688,7 +720,10 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
call psb_spmm(-done,baseprecv(1)%base_a,mlprec_wrk(1)%y2l,&
& done,mlprec_wrk(1)%tx,baseprecv(1)%base_desc,info,work=work)
if(info /=0) goto 9999
if (info /=0) then
call psb_errpush(4001,name,a_err=' fine level residual')
goto 9999
end if
!
! STEP 4
@ -715,7 +750,7 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
mlprec_wrk(ilev)%x2l(:) = dzero
mlprec_wrk(ilev)%y2l(:) = dzero
mlprec_wrk(ilev)%tx(:) = dzero
mlprec_wrk(ilev)%tx(:) = dzero
if (ismth /= mld_no_smooth_) then
!
@ -723,12 +758,8 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
call psb_halo(mlprec_wrk(ilev-1)%tx,baseprecv(ilev-1)%base_desc,&
& info,work=work)
if(info /=0) goto 9999
call psb_csmm(done,baseprecv(ilev)%av(mld_sm_pr_t_),mlprec_wrk(ilev-1)%tx,dzero,&
& mlprec_wrk(ilev)%x2l,info)
if(info /=0) goto 9999
if (info == 0) call psb_csmm(done,baseprecv(ilev)%av(mld_sm_pr_t_),&
& mlprec_wrk(ilev-1)%tx,dzero,mlprec_wrk(ilev)%x2l,info)
else
!
! Apply the raw aggregation map transpose (take a shortcut)
@ -740,11 +771,18 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
& mlprec_wrk(ilev-1)%tx(i)
end do
end if
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction')
goto 9999
end if
if (icm ==mld_repl_mat_) then
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nr2l))
else if (icm /= mld_distr_mat_) then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(mld_coarse_mat_) ', icm
info = 4013
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_',&
& i_Err=(/icm,0,0,0,0/))
goto 9999
endif
!
@ -753,18 +791,19 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
call mld_baseprec_aply(done,baseprecv(ilev),mlprec_wrk(ilev)%x2l,&
& dzero,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%desc_data, 'N',work,info)
if(info /= 0) goto 9999
!
! Compute the residual (at all levels but the coarsest one)
!
if (ilev < nlev) then
mlprec_wrk(ilev)%tx = mlprec_wrk(ilev)%x2l
call psb_spmm(-done,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
& done,mlprec_wrk(ilev)%tx,baseprecv(ilev)%base_desc,info,work=work)
if(info /=0) goto 9999
if (info == 0) call psb_spmm(-done,baseprecv(ilev)%base_a,&
& mlprec_wrk(ilev)%y2l,done,mlprec_wrk(ilev)%tx,&
& baseprecv(ilev)%base_desc,info,work=work)
endif
if (info /=0) then
call psb_errpush(4001,name,a_err='Error on up sweep residual')
goto 9999
end if
enddo
!
@ -784,10 +823,8 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
if (ismth == mld_smooth_prol_) &
& call psb_halo(mlprec_wrk(ilev+1)%y2l,&
& baseprecv(ilev+1)%desc_data,info,work=work)
call psb_csmm(done,baseprecv(ilev+1)%av(mld_sm_pr_),mlprec_wrk(ilev+1)%y2l,&
& done,mlprec_wrk(ilev)%y2l,info)
if(info /=0) goto 9999
if (info == 0) call psb_csmm(done,baseprecv(ilev+1)%av(mld_sm_pr_),&
& mlprec_wrk(ilev+1)%y2l,done,mlprec_wrk(ilev)%y2l,info)
else
!
! Apply the raw aggregation map (take a shortcut)
@ -796,9 +833,11 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
mlprec_wrk(ilev)%y2l(i) = mlprec_wrk(ilev)%y2l(i) + &
& mlprec_wrk(ilev+1)%y2l(baseprecv(ilev+1)%mlia(i))
enddo
end if
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during prolongation')
goto 9999
end if
enddo
!
@ -808,8 +847,11 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,&
& baseprecv(1)%base_desc,info)
if (info /=0) then
call psb_errpush(4001,name,a_err='Error on final update')
goto 9999
end if
if(info /=0) goto 9999
case(mld_twoside_smooth_)
@ -891,18 +933,18 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
& dzero,mlprec_wrk(1)%y2l,&
& baseprecv(1)%base_desc,&
& trans,work,info)
if(info /=0) goto 9999
!
! STEP 3
!
! Compute the residual at the finest level
!
mlprec_wrk(1)%ty = mlprec_wrk(1)%x2l
call psb_spmm(-done,baseprecv(1)%base_a,mlprec_wrk(1)%y2l,&
if (info == 0) call psb_spmm(-done,baseprecv(1)%base_a,mlprec_wrk(1)%y2l,&
& done,mlprec_wrk(1)%ty,baseprecv(1)%base_desc,info,work=work)
if(info /=0) goto 9999
if (info /=0) then
call psb_errpush(4010,name,a_err='Fine level baseprec/residual')
goto 9999
end if
!
! STEP 4
@ -938,11 +980,8 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
call psb_halo(mlprec_wrk(ilev-1)%ty,baseprecv(ilev-1)%base_desc,&
& info,work=work)
if(info /=0) goto 9999
call psb_csmm(done,baseprecv(ilev)%av(mld_sm_pr_t_),mlprec_wrk(ilev-1)%ty,dzero,&
& mlprec_wrk(ilev)%x2l,info)
if(info /=0) goto 9999
if (info == 0) call psb_csmm(done,baseprecv(ilev)%av(mld_sm_pr_t_),&
& mlprec_wrk(ilev-1)%ty,dzero,mlprec_wrk(ilev)%x2l,info)
else
!
! Apply the raw aggregation map transpose (take a shortcut)
@ -954,34 +993,41 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
& mlprec_wrk(ilev-1)%ty(i)
end do
end if
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction')
goto 9999
end if
if (icm == mld_repl_mat_) then
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nr2l))
else if (icm /= mld_distr_mat_) then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(mld_coarse_mat_) ', icm
info = 4013
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_',&
& i_Err=(/icm,0,0,0,0/))
goto 9999
endif
call psb_geaxpby(done,mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%tx,&
& baseprecv(ilev)%base_desc,info)
if(info /=0) goto 9999
!
! Apply the base preconditioner
!
call mld_baseprec_aply(done,baseprecv(ilev),mlprec_wrk(ilev)%x2l,&
& dzero,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%desc_data, 'N',work,info)
if(info /=0) goto 9999
if (info == 0) call mld_baseprec_aply(done,baseprecv(ilev),&
& mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%y2l,&
& baseprecv(ilev)%desc_data, 'N',work,info)
!
! Compute the residual (at all levels but the coarsest one)
!
if(ilev < nlev) then
mlprec_wrk(ilev)%ty = mlprec_wrk(ilev)%x2l
call psb_spmm(-done,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
& done,mlprec_wrk(ilev)%ty,baseprecv(ilev)%base_desc,info,work=work)
if(info /=0) goto 9999
if (info == 0) call psb_spmm(-done,baseprecv(ilev)%base_a,&
& mlprec_wrk(ilev)%y2l,done,mlprec_wrk(ilev)%ty,&
& baseprecv(ilev)%base_desc,info,work=work)
endif
if (info /=0) then
call psb_errpush(4001,name,a_err='baseprec_aply/residual')
goto 9999
end if
enddo
@ -1002,10 +1048,8 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
if (ismth == mld_smooth_prol_) &
& call psb_halo(mlprec_wrk(ilev+1)%y2l,baseprecv(ilev+1)%desc_data,&
& info,work=work)
call psb_csmm(done,baseprecv(ilev+1)%av(mld_sm_pr_),mlprec_wrk(ilev+1)%y2l,&
& done,mlprec_wrk(ilev)%y2l,info)
if(info /=0) goto 9999
if (info == 0) call psb_csmm(done,baseprecv(ilev+1)%av(mld_sm_pr_),&
& mlprec_wrk(ilev+1)%y2l, done,mlprec_wrk(ilev)%y2l,info)
else
!
! Apply the raw aggregation map (take a shortcut)
@ -1014,7 +1058,10 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
mlprec_wrk(ilev)%y2l(i) = mlprec_wrk(ilev)%y2l(i) + &
& mlprec_wrk(ilev+1)%y2l(baseprecv(ilev+1)%mlia(i))
enddo
end if
if (info /=0 ) then
call psb_errpush(4001,name,a_err='Error during restriction')
goto 9999
end if
!
@ -1022,17 +1069,15 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
call psb_spmm(-done,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
& done,mlprec_wrk(ilev)%tx,baseprecv(ilev)%base_desc,info,work=work)
if(info /=0) goto 9999
!
! Apply the base preconditioner
!
call mld_baseprec_aply(done,baseprecv(ilev),mlprec_wrk(ilev)%tx,&
if (info == 0) call mld_baseprec_aply(done,baseprecv(ilev),mlprec_wrk(ilev)%tx,&
& done,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%base_desc, trans, work,info)
if(info /=0) goto 9999
if (info /= 0) then
call psb_errpush(4001,name,a_err='Error: residual/baseprec_aply')
goto 9999
end if
enddo
!
@ -1043,30 +1088,37 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,&
& baseprecv(1)%base_desc,info)
if(info /=0) goto 9999
if (info /= 0) then
call psb_errpush(4001,name,a_err='Error final update')
goto 9999
end if
case default
call psb_errpush(4013,name,a_err='wrong smooth_pos',&
info = 4013
call psb_errpush(info,name,a_err='invalid smooth_pos',&
& i_Err=(/baseprecv(2)%iprcparm(mld_smooth_pos_),0,0,0,0/))
goto 9999
end select
case default
call psb_errpush(4013,name,a_err='wrong mltype',&
info = 4013
call psb_errpush(info,name,a_err='invalid mltype',&
& i_Err=(/baseprecv(2)%iprcparm(mld_ml_type_),0,0,0,0/))
goto 9999
end select
deallocate(mlprec_wrk)
deallocate(mlprec_wrk,stat=info)
if (info /= 0) then
call psb_errpush(4000,name)
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_errpush(info,name)
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()

@ -73,16 +73,15 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info)
integer :: err_act
character(len=20) :: name, ch_err
logical, parameter :: debug=.false.
type(psb_dspmat_type) :: ac
integer :: ictxt, np, me
name='psb_dmlprec_bld'
if (psb_get_errstatus().ne.0) return
call psb_erractionsave(err_act)
info = 0
ictxt = psb_cd_get_context(desc_a)
call psb_info(ictxt,me,np)
call psb_erractionsave(err_act)
if (.not.allocated(p%iprcparm)) then
info = 2222
@ -120,14 +119,10 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info)
!
call mld_aggrmap_bld(p%iprcparm(mld_aggr_alg_),a,desc_a,p%nlaggr,p%mlia,info)
if(info /= 0) then
info=4010
ch_err='mld_aggrmap_bld'
call psb_errpush(info,name,a_err=ch_err)
call psb_errpush(4010,name,a_err='mld_aggrmap_bld')
goto 9999
end if
if (debug) write(0,*) 'Out from genaggrmap',p%nlaggr
!
! Build the coarse-level matrix from the fine level one, starting from
! the mapping defined by mld_aggrmap_bld and applying the aggregation
@ -137,22 +132,16 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info)
call psb_nullify_desc(desc_ac)
call mld_aggrmat_asb(a,desc_a,ac,desc_ac,p,info)
if(info /= 0) then
info=4010
ch_err='mld_aggrmat_asb'
call psb_errpush(info,name,a_err=ch_err)
call psb_errpush(4010,name,a_err='mld_aggrmat_asb')
goto 9999
end if
if (debug) write(0,*) 'Out from bldaggrmat',desc_ac%matrix_data(:)
!
! Build the 'base preconditioner' corresponding to the coarse level
!
call mld_baseprc_bld(ac,desc_ac,p,info)
if (debug) write(0,*) 'Out from baseprcbld',info
if(info /= 0) then
info=4010
ch_err='mld_baseprc_bld'
call psb_errpush(info,name,a_err=ch_err)
if (info /= 0) then
call psb_errpush(4010,name,a_err='mld_baseprc_bld')
goto 9999
end if
@ -165,12 +154,10 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info)
!
call psb_sp_transfer(ac,p%av(mld_ac_),info)
p%base_a => p%av(mld_ac_)
call psb_cdtransfer(desc_ac,p%desc_ac,info)
if (info==0) call psb_cdtransfer(desc_ac,p%desc_ac,info)
if (info /= 0) then
info=4010
ch_err='psb_cdtransfer'
call psb_errpush(info,name,a_err=ch_err)
call psb_errpush(4010,name,a_err='psb_cdtransfer')
goto 9999
end if
p%base_desc => p%desc_ac

@ -50,7 +50,7 @@
!
!
! Arguments:
! prec - type(mld_dprec_type), input.
! prec - type(mld_dprec_type), input.
! The preconditioner data structure containing the local part
! of the preconditioner to be applied.
! x - real(kind(0.d0)), dimension(:), input.
@ -66,7 +66,8 @@
! If trans='N','n' then op(M^(-1)) = M^(-1);
! if trans='T','t' then op(M^(-1)) = M^(-T) (transpose of M^(-1)).
! work - real(kind(0.d0)), dimension (:), optional, target.
! Workspace. Its size must be at least 4*psb_cd_get_local_cols(desc_data).
! Workspace. Its size must be at
! least 4*psb_cd_get_local_cols(desc_data).
!
subroutine mld_dprec_aply(prec,x,y,desc_data,info,trans,work)
@ -88,7 +89,6 @@ subroutine mld_dprec_aply(prec,x,y,desc_data,info,trans,work)
character :: trans_
real(kind(1.d0)), pointer :: work_(:)
integer :: ictxt,np,me,err_act,iwsz
logical,parameter :: debug=.false., debugprt=.false.
character(len=20) :: name
name='mld_dprec_aply'
@ -118,10 +118,12 @@ subroutine mld_dprec_aply(prec,x,y,desc_data,info,trans,work)
end if
if (.not.(allocated(prec%baseprecv))) then
write(0,*) 'Inconsistent preconditioner: neither ML nor BASE?'
!! Error 1: should call mld_dprecbld
info=3112
call psb_errpush(info,name)
goto 9999
end if
if (size(prec%baseprecv) >1) then
if (debug) write(0,*) 'Into mlprec_aply',size(x),size(y)
call mld_mlprec_aply(done,prec%baseprecv,x,dzero,y,desc_data,trans_,work_,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='mld_dmlprec_aply')
@ -131,7 +133,10 @@ subroutine mld_dprec_aply(prec,x,y,desc_data,info,trans,work)
else if (size(prec%baseprecv) == 1) then
call mld_baseprec_aply(done,prec%baseprecv(1),x,dzero,y,desc_data,trans_, work_,info)
else
write(0,*) 'Inconsistent preconditioner: size of baseprecv???'
info = 4013
call psb_errpush(info,name,a_err='Invalid size of baseprecv',&
& i_Err=(/size(prec%baseprecv),0,0,0,0/))
goto 9999
endif
if (present(work)) then
@ -171,7 +176,7 @@ end subroutine mld_dprec_aply
!
!
! Arguments:
! prec - type(mld_dprec_type), input.
! prec - type(mld_dprec_type), input.
! The preconditioner data structure containing the local part
! of the preconditioner to be applied.
! x - real(kind(0.d0)), dimension(:), input/output.
@ -200,9 +205,8 @@ subroutine mld_dprec_aply1(prec,x,desc_data,info,trans)
character(len=1), optional :: trans
! Local variables
logical,parameter :: debug=.false., debugprt=.false.
character :: trans_
integer :: ictxt,np,me, err_act
integer :: ictxt,np,me, err_act
real(kind(1.d0)), pointer :: WW(:), w1(:)
character(len=20) :: name
@ -226,17 +230,25 @@ subroutine mld_dprec_aply1(prec,x,desc_data,info,trans)
& a_err='real(kind(1.d0))')
goto 9999
end if
if (debug) write(0,*) 'prec_aply1 size(x) ',size(x), size(ww),size(w1)
call mld_dprec_aply(prec,x,ww,desc_data,info,trans_,work=w1)
if(info /=0) goto 9999
call mld_precaply(prec,x,ww,desc_data,info,trans_,work=w1)
if (info /= 0) then
call psb_errpush(4010,name,a_err='mld_precaply')
goto 9999
end if
x(:) = ww(:)
deallocate(ww,W1)
deallocate(ww,W1,stat=info)
if (info /= 0) then
info = 4000
call psb_errpush(info,name)
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_errpush(info,name)
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()

@ -41,7 +41,7 @@
! Contains: subroutine init_baseprc_av
!
! This routine builds the preconditioner according to the requirements made by
! the user trough the subroutines mld_dprecinit and mld_dprecset.
! the user trough the subroutines mld_precinit and mld_precset.
!
! A multilevel preconditioner is regarded as an array of 'base preconditioners',
! each representing the part of the preconditioner associated to a certain level.
@ -76,32 +76,36 @@ subroutine mld_dprecbld(a,desc_a,p,info,upd)
character, intent(in), optional :: upd
! Local Variables
Integer :: err,i,k,ictxt, me,np, err_act
Integer :: err,i,k,ictxt, me,np, err_act, iszv
integer :: int_err(5)
character :: iupd
logical, parameter :: debug=.false.
integer,parameter :: iroot=0,iout=60,ilout=40
character(len=20) :: name, ch_err
integer :: debug_level, debug_unit
character(len=20) :: name, ch_err
if(psb_get_errstatus().ne.0) return
if (psb_get_errstatus().ne.0) return
info=0
err=0
call psb_erractionsave(err_act)
name = 'mld_dprecbld'
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
if (debug) write(0,*) 'Entering precbld',desc_a%matrix_data(:)
name = 'mld_dprecbld'
info = 0
int_err(1) = 0
ictxt = psb_cd_get_context(desc_a)
if (debug) write(0,*) 'Preconditioner psb_info'
call psb_info(ictxt, me, np)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Entering ',desc_a%matrix_data(:)
if (present(upd)) then
if (debug) write(0,*) 'UPD ', upd
if ((upd.eq.'F').or.(upd.eq.'T')) then
iupd=upd
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),'UPD ', upd
if ((toupper(upd).eq.'F').or.(toupper(upd).eq.'T')) then
iupd=toupper(upd)
else
iupd='F'
endif
@ -110,86 +114,81 @@ subroutine mld_dprecbld(a,desc_a,p,info,upd)
endif
if (.not.allocated(p%baseprecv)) then
!! Error 1: should call mld_dprecset
info=4010
ch_err='unallocated bpv'
call psb_errpush(info,name,a_err=ch_err)
!! Error: should have called mld_dprecinit
info=3111
call psb_errpush(info,name)
goto 9999
end if
!
! Should add check to ensure all procs have the same ...
! Check to ensure all procs have the same
!
iszv = size(p%baseprecv)
call psb_bcast(ictxt,iszv)
if (iszv /= size(p%baseprecv)) then
info=4001
call psb_errpush(info,name,a_err='Inconsistent size of baseprecv')
goto 9999
end if
if (size(p%baseprecv) >= 1) then
if (iszv >= 1) then
!
! Allocate the av component of the preconditioner data type
! at the finest level
! Allocate and build the fine level preconditioner
!
call init_baseprc_av(p%baseprecv(1),info)
if (info == 0) call mld_baseprc_bld(a,desc_a,p%baseprecv(1),info,iupd)
if (info /= 0) then
info=4010
ch_err='allocate'
call psb_errpush(info,name,a_err=ch_err)
call psb_errpush(4001,name,a_err='Base level precbuild.')
goto 9999
endif
!
! Build the base preconditioner corresponding to the finest
! level
!
call mld_baseprc_bld(a,desc_a,p%baseprecv(1),info,iupd)
end if
else
info=4010
ch_err='size bpv'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
if (size(p%baseprecv) > 1) then
if (iszv > 1) then
!
! Build the base preconditioners corresponding to the remaining
! levels
!
do i=2, size(p%baseprecv)
do i=2, iszv
!
! Allocate the av component of the preconditioner data type
! at level i
!
call init_baseprc_av(p%baseprecv(i),info)
if (info /= 0) then
info=4010
ch_err='allocate'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
if (i<size(p%baseprecv)) then
if (i<iszv) then
!
! A replicated matrix only makes sense at the coarsest level
!
call mld_check_def(p%baseprecv(i)%iprcparm(mld_coarse_mat_),'Coarse matrix',&
& mld_distr_mat_,is_distr_ml_coarse_mat)
end if
call init_baseprc_av(p%baseprecv(i),info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Calling mlprcbld at level ',i
!
! Build the base preconditioner corresponding to level i
!
call mld_mlprec_bld(p%baseprecv(i-1)%base_a,p%baseprecv(i-1)%base_desc,&
& p%baseprecv(i),info)
if (info == 0) call mld_mlprec_bld(p%baseprecv(i-1)%base_a,&
& p%baseprecv(i-1)%base_desc, p%baseprecv(i),info)
if (info /= 0) then
info=4010
call psb_errpush(info,name)
call psb_errpush(4001,name,a_err='Init & build upper level preconditioner')
goto 9999
endif
if (debug) then
write(0,*) 'Return from ',i-1,' call to mlprcbld ',info
endif
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Return from ',i,' call to mlprcbld ',info
end do
endif
@ -211,12 +210,15 @@ contains
type(mld_dbaseprc_type), intent(inout) :: p
integer :: info
if (allocated(p%av)) then
!
! We have not yet decided what to do
!
if (size(p%av) /= mld_max_avsz_) then
deallocate(p%av,stat=info)
if (info /= 0) return
endif
end if
if (.not.(allocated(p%av))) then
allocate(p%av(mld_max_avsz_),stat=info)
if (info /= 0) return
end if
allocate(p%av(mld_max_avsz_),stat=info)
!!$ if (info /= 0) return
do k=1,size(p%av)
call psb_nullify_sp(p%av(k))
end do

@ -182,9 +182,6 @@ subroutine mld_dprecinit(p,ptype,info,nlev)
else
nlev_ = 2
end if
if (nlev_ == 1) then
write(0,*) 'Warning: requested ML preconditioner with NLEV=1'
endif
ilev_ = 1
allocate(p%baseprecv(nlev_),stat=info)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)

@ -87,8 +87,7 @@ subroutine mld_dprecseti(p,what,val,info,ilev)
info = 0
if (.not.allocated(p%baseprecv)) then
write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
info = -1
info = 3111
return
endif
nlev_ = size(p%baseprecv)
@ -100,13 +99,11 @@ subroutine mld_dprecseti(p,what,val,info,ilev)
end if
if ((ilev_<1).or.(ilev_ > nlev_)) then
write(0,*) 'PRECSET ERRROR: ilev out of bounds'
info = -1
return
endif
if (.not.allocated(p%baseprecv(ilev_)%iprcparm)) then
write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
info = -1
info = 3111
return
endif
@ -251,8 +248,7 @@ subroutine mld_dprecsetc(p,what,string,info,ilev)
info = 0
if (.not.allocated(p%baseprecv)) then
write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
info = -1
info = 3111
return
endif
nlev_ = size(p%baseprecv)
@ -269,8 +265,7 @@ subroutine mld_dprecsetc(p,what,string,info,ilev)
return
endif
if (.not.allocated(p%baseprecv(ilev_)%iprcparm)) then
write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
info = -1
info = 3111
return
endif
@ -445,8 +440,7 @@ subroutine mld_dprecsetd(p,what,val,info,ilev)
end if
if (.not.allocated(p%baseprecv)) then
write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
info = -1
info = 3111
return
endif
nlev_ = size(p%baseprecv)
@ -457,8 +451,7 @@ subroutine mld_dprecsetd(p,what,val,info,ilev)
return
endif
if (.not.allocated(p%baseprecv(ilev_)%dprcparm)) then
write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
info = -1
info = 3111
return
endif

@ -82,7 +82,6 @@ subroutine mld_dslu_bld(a,desc_a,p,info)
! Local variables
integer :: nzt,ictxt,me,np,err_act
logical, parameter :: debug=.false.
character(len=20) :: name, ch_err
if(psb_get_errstatus().ne.0) return
@ -95,19 +94,12 @@ subroutine mld_dslu_bld(a,desc_a,p,info)
call psb_info(ictxt, me, np)
if (toupper(a%fida) /= 'CSR') then
write(0,*) 'Unimplemented input to mld_slu_BLD'
info=135
call psb_errpush(info,name,a_err=a%fida)
goto 9999
endif
nzt = psb_sp_get_nnzeros(a)
if (Debug) then
write(0,*) me,'Calling mld_slu_factor ',nzt,a%m,&
& a%k,p%desc_data%matrix_data(psb_n_row_)
call psb_barrier(ictxt)
endif
!
! Compute the LU factorization
!
@ -120,11 +112,6 @@ subroutine mld_dslu_bld(a,desc_a,p,info)
goto 9999
end if
if (Debug) then
write(0,*) me, 'SPLUBLD: Done mld_slu_Factor',info,p%iprcparm(mld_slu_ptr_)
call psb_barrier(ictxt)
endif
call psb_erractionrestore(err_act)
return

@ -80,7 +80,6 @@ subroutine mld_dsludist_bld(a,desc_a,p,info)
! Local variables
integer :: nzt,ictxt,me,np,err_act,&
& mglob,ifrst,ibcheck,nrow,ncol,npr,npc
logical, parameter :: debug=.false.
character(len=20) :: name, ch_err
if (psb_get_errstatus().ne.0) return
@ -93,7 +92,8 @@ subroutine mld_dsludist_bld(a,desc_a,p,info)
call psb_info(ictxt, me, np)
if (toupper(a%fida) /= 'CSR') then
write(0,*) 'Unimplemented input to mld_slu_BLD'
info=135
call psb_errpush(info,name,a_err=a%fida)
goto 9999
endif

@ -87,9 +87,8 @@ subroutine mld_dumf_bld(a,desc_a,p,info)
integer, intent(out) :: info
! Local variables
integer :: nzt,ictxt,me,np,err_act
integer :: i_err(5)
logical, parameter :: debug=.false.
integer :: nzt,ictxt,me,np,err_act
integer :: i_err(5)
character(len=20) :: name
info=0
@ -104,18 +103,8 @@ subroutine mld_dumf_bld(a,desc_a,p,info)
goto 9999
endif
nzt = psb_sp_get_nnzeros(a)
if (Debug) then
write(0,*) me,'Calling mld_umf_factor ',nzt,a%m,&
& a%k,p%desc_data%matrix_data(psb_n_row_)
open(80+me)
call psb_csprt(80+me,a)
close(80+me)
call psb_barrier(ictxt)
endif
!
! Compute the LU factorization
!
@ -130,11 +119,6 @@ subroutine mld_dumf_bld(a,desc_a,p,info)
goto 9999
end if
if (Debug) then
write(0,*) me, 'UMFBLD: Done mld_umf_Factor',info,p%iprcparm(mld_umf_numptr_)
call psb_barrier(ictxt)
endif
call psb_erractionrestore(err_act)
return

@ -41,7 +41,7 @@
!
! This routine builds a mapping from the row indices of the fine-level matrix
! to the row indices of the coarse-level matrix, according to a decoupled
! aggregation algorithm. This mapping will be used by mld_daggrmat_asb to
! aggregation algorithm. This mapping will be used by mld_aggrmat_asb to
! build the coarse-level matrix.
!
! The aggregation algorithm is a parallel version of that described in
@ -93,16 +93,17 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
type(psb_zspmat_type), pointer :: apnt
logical :: recovery
logical, parameter :: debug=.false.
integer ::ictxt,np,me,err_act
integer :: debug_level, debug_unit
integer :: ictxt,np,me,err_act
integer :: nrow, ncol, n_ne
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() /= 0) return
info=0
name = 'mld_aggrmap_bld'
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!
! Note. At the time being we are ignoring aggr_type so
! that we only have decoupled aggregation. This might
@ -116,10 +117,9 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
select case (aggr_type)
case (mld_dec_aggr_,mld_sym_dec_aggr_)
nr = a%m
allocate(ilaggr(nr),neigh(nr),stat=info)
if(info.ne.0) then
if(info /= 0) then
info=4025
call psb_errpush(info,name,i_err=(/2*nr,0,0,0,0/),&
& a_err='integer')
@ -136,13 +136,19 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
& 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.)
if (info == 0) call psb_transp(atmp,atrans,fmt='COO')
if (info == 0) 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)
if (info == 0) call psb_sp_free(atrans,info)
if (info == 0) call psb_ipcoo2csr(atmp,info)
apnt => atmp
if (info/=0) then
info=4001
call psb_errpush(info,name,a_err='init apnt')
goto 9999
end if
end if
@ -168,11 +174,10 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
naggr = naggr + 1
ilaggr(i) = naggr
call psb_neigh(apnt,i,neigh,n_ne,info,lev=one)
call psb_neigh(apnt,i,neigh,n_ne,info,lev=1)
if (info/=0) then
info=4010
ch_err='psb_neigh'
call psb_errpush(info,name,a_err=ch_err)
call psb_errpush(info,name,a_err='psb_neigh')
goto 9999
end if
do k=1, n_ne
@ -184,11 +189,10 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
!
! 2. Untouched neighbours of these nodes are marked <0.
!
call psb_neigh(apnt,i,neigh,n_ne,info,lev=two)
call psb_neigh(apnt,i,neigh,n_ne,info,lev=2)
if (info/=0) then
info=4010
ch_err='psb_neigh'
call psb_errpush(info,name,a_err=ch_err)
call psb_errpush(info,name,a_err='psb_neigh')
goto 9999
end if
@ -203,15 +207,17 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
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)
if (debug_level >= psb_debug_outer_) then
write(debug_unit,*) me,' ',trim(name),&
& ' 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
if(info /= 0) then
info=4025
call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),&
& a_err='integer')
@ -225,16 +231,20 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
n = ilaggr(i)
if (n>0) then
if (n>naggr) then
write(0,*) 'loc_Aggregate: n > naggr 1 ? ',n,naggr
info=4001
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 1 ?')
goto 9999
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)
if (debug_level >= psb_debug_outer_) then
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 1: number of aggregates ',naggr
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 1: nodes aggregated ',sum(ils)
end if
recovery=.false.
@ -247,11 +257,10 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
isz = nr+1
ia = -1
call psb_neigh(apnt,i,neigh,n_ne,info,lev=one)
call psb_neigh(apnt,i,neigh,n_ne,info,lev=1)
if (info/=0) then
info=4010
ch_err='psb_neigh'
call psb_errpush(info,name,a_err=ch_err)
call psb_errpush(info,name,a_err='psb_neigh')
goto 9999
end if
@ -261,7 +270,9 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
n = ilaggr(k)
if (n>0) then
if (n>naggr) then
write(0,*) 'loc_Aggregate: n > naggr 2? ',n,naggr
info=4001
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?')
goto 9999
end if
if (ils(n) < isz) then
@ -275,7 +286,9 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
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
info=4001
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 3 ?')
goto 9999
end if
ils(ilaggr(i)) = ils(ilaggr(i)) + 1
!
@ -284,35 +297,44 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
!
recovery = .true.
else
write(0,*) 'Unrecoverable error !!',ilaggr(i), nr
info=4001
call psb_errpush(info,name,a_err='Unrecoverable error !!')
goto 9999
endif
else
ilaggr(i) = ia
if (ia>naggr) then
write(0,*) 'loc_Aggregate: n > naggr 4? ',ia,naggr
info=4001
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 4? ')
goto 9999
end if
ils(ia) = ils(ia) + 1
endif
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)
if (debug_level >= psb_debug_outer_) then
if (recovery) then
write(debug_unit,*) me,' ',trim(name),&
& 'Had to recover from strange situation in loc_aggregate.'
write(debug_unit,*) me,' ',trim(name),&
& 'Perhaps an unsymmetric pattern?'
endif
write(debug_unit,*) me,' ',trim(name),&
& 'Phase 2: number of aggregates ',naggr,sum(ils)
do i=1, naggr
write(*,*) 'Size of aggregate ',i,' :',count(ilaggr==i), ils(i)
write(debug_unit,*) me,' ',trim(name),&
& 'Size of aggregate ',i,' :',count(ilaggr==i), ils(i)
enddo
write(*,*) maxval(ils(1:naggr))
write(*,*) 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops'
write(debug_unit,*) me,' ',trim(name),&
& maxval(ils(1:naggr))
write(debug_unit,*) me,' ',trim(name),&
& 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops'
end if
!!$ write(0,*) 'desc_a loc_aggr 4 : ', desc_a%matrix_data(m_)
if (count(ilaggr<0) >0) then
write(0,*) 'Fatal error: some leftovers!!!'
info=4001
call psb_errpush(info,name,a_err='Fatal error: some leftovers')
goto 9999
endif
deallocate(ils,neigh,stat=info)
@ -322,9 +344,6 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
goto 9999
end if
if (nrow /= size(ilaggr)) then
write(0,*) 'SOmething wrong ilaggr ',nrow,size(ilaggr)
endif
call psb_realloc(ncol,ilaggr,info)
if (info/=0) then
info=4010
@ -351,7 +370,6 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
case default
write(0,*) 'Unimplemented aggregation algorithm ',aggr_type
info = -1
call psb_errpush(30,name,i_err=(/1,aggr_type,0,0,0/))
goto 9999

@ -107,9 +107,8 @@ subroutine mld_zaggrmat_asb(a,desc_a,ac,desc_ac,p,info)
integer, intent(out) :: info
! Local variables
logical, parameter :: aggr_dump=.false.
integer :: ictxt,np,me, err_act,icomm
character(len=20) :: name
integer ::ictxt,np,me, err_act, icomm
character(len=20) :: name
name='mld_aggrmat_asb'
if(psb_get_errstatus().ne.0) return
@ -125,24 +124,22 @@ subroutine mld_zaggrmat_asb(a,desc_a,ac,desc_ac,p,info)
case (mld_no_smooth_)
call mld_aggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='raw_aggregate')
call psb_errpush(4010,name,a_err='mld_aggrmat_raw_asb')
goto 9999
end if
if (aggr_dump) call psb_csprt(90+me,ac,head='% Raw aggregate.')
case(mld_smooth_prol_,mld_biz_prol_)
if (aggr_dump) call psb_csprt(70+me,a,head='% Input matrix')
call mld_aggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
call mld_aggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='smooth_aggregate')
call psb_errpush(4010,name,a_err='mld_aggrmat_smth_asb')
goto 9999
end if
if (aggr_dump) call psb_csprt(90+me,ac,head='% Smooth aggregate.')
case default
call psb_errpush(4010,name,a_err=name)
call psb_errpush(4001,name,a_err='Invalid aggr kind')
goto 9999
end select

@ -95,8 +95,7 @@ subroutine mld_zaggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
integer, intent(out) :: info
! Local variables
logical, parameter :: aggr_dump=.false.
integer ::ictxt,np,me, err_act,icomm
integer ::ictxt,np,me, err_act, icomm
character(len=20) :: name
type(psb_zspmat_type) :: b
integer, pointer :: nzbr(:), idisp(:)
@ -202,7 +201,7 @@ subroutine mld_zaggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
ac%k = ntaggr
ac%infoa(psb_nnz_) = nzac
ac%fida='COO'
ac%descra='G'
ac%descra='GUN'
call psb_spcnv(ac,info,afmt='coo',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='sp_free')
@ -212,19 +211,10 @@ subroutine mld_zaggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
else if (p%iprcparm(mld_coarse_mat_) == mld_distr_mat_) then
call psb_cdall(ictxt,desc_ac,info,nl=naggr)
if (info == 0) call psb_cdasb(desc_ac,info)
if (info == 0) call psb_sp_clone(b,ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdall')
goto 9999
end if
call psb_cdasb(desc_ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdasb')
goto 9999
end if
call psb_sp_clone(b,ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spclone')
call psb_errpush(4001,name,a_err='Build ac, desc_ac')
goto 9999
end if
call psb_sp_free(b,info)
@ -234,8 +224,9 @@ subroutine mld_zaggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
end if
else
write(0,*) 'Unknown p%iprcparm(coarse_mat) in aggregate_sp',p%iprcparm(mld_coarse_mat_)
info = 4001
call psb_errpush(4001,name,a_err='invalid mld_coarse_mat_')
goto 9999
end if
deallocate(nzbr,idisp)

@ -117,14 +117,15 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
integer, pointer :: nzbr(:), idisp(:)
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k
logical, parameter :: aggr_dump=.false.
integer ::ictxt,np,me, err_act, icomm
character(len=20) :: name
type(psb_zspmat_type), pointer :: am1,am2
type(psb_zspmat_type) :: am3,am4
logical :: ml_global_nmb
logical, parameter :: test_dump=.false., debug=.false.
integer :: nz
integer, allocatable :: ia(:), ja(:)
complex(kind(1.d0)), allocatable :: val(:)
integer :: debug_level, debug_unit
integer, parameter :: ncmax=16
real(kind(1.d0)) :: omega, anorm, tmp, dg
@ -132,10 +133,12 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
if(psb_get_errstatus().ne.0) return
info=0
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = psb_cd_get_context(desc_a)
icomm = psb_cd_get_mpic(desc_a)
ictxt=psb_cd_get_context(desc_a)
ictxt = psb_cd_get_context(desc_a)
call psb_info(ictxt, me, np)
@ -164,38 +167,27 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
goto 9999
end if
naggrm1 = sum(p%nlaggr(1:me))
naggrp1 = sum(p%nlaggr(1:me+1))
ml_global_nmb = ( (p%iprcparm(mld_aggr_kind_) == mld_smooth_prol_).or.&
& ( (p%iprcparm(mld_aggr_kind_) == mld_biz_prol_).and.&
& (p%iprcparm(mld_coarse_mat_) == mld_repl_mat_)) )
if (ml_global_nmb) then
p%mlia(1:nrow) = p%mlia(1:nrow) + naggrm1
call psb_halo(p%mlia,desc_a,info)
if(info /= 0) then
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_halo')
goto 9999
end if
end if
if (aggr_dump) then
open(30+me)
write(30+me,*) '% Aggregation map'
do i=1,ncol
write(30+me,*) i,p%mlia(i)
end do
close(30+me)
end if
! naggr: number of local aggregates
! nrow: local rows.
!
allocate(p%dorig(nrow),stat=info)
if (info /= 0) then
info=4025
call psb_errpush(info,name,i_err=(/nrow,0,0,0,0/),&
@ -218,13 +210,6 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
end if
end do
! where (p%dorig /= zzero)
! p%dorig = zone / p%dorig
! elsewhere
! p%dorig = zone
! end where
! 1. Allocate Ptilde in sparse matrix form
am4%fida='COO'
am4%m=ncol
@ -235,7 +220,8 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
am4%k=naggr
call psb_sp_all(ncol,naggr,am4,ncol,info)
endif
if(info /= 0) then
if (info /= 0) then
call psb_errpush(4010,name,a_err='spall')
goto 9999
end if
@ -258,24 +244,19 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
call psb_spcnv(am4,info,afmt='csr',dupl=psb_dupl_add_)
if(info /= 0) then
if (info==0) call psb_spcnv(a,am3,info,afmt='csr',dupl=psb_dupl_add_)
if (info /= 0) then
call psb_errpush(4010,name,a_err='spcnv')
goto 9999
end if
call psb_sp_clone(a,am3,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spclone')
goto 9999
end if
!
! WARNING: the cycles below assume that AM3 does have
! its diagonal elements stored explicitly!!!
! Should we switch to something safer?
!
call psb_sp_scal(am3,p%dorig,info)
if(info /= 0) goto 9999
if (info /= 0) goto 9999
if (p%iprcparm(mld_aggr_eig_) == mld_max_norm_) then
@ -284,25 +265,33 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
!
! This only works with CSR.
!
anorm = dzero
dg = done
do i=1,am3%m
tmp = dzero
do j=am3%ia2(i),am3%ia2(i+1)-1
if (am3%ia1(j) <= am3%m) then
tmp = tmp + abs(am3%aspk(j))
endif
if (am3%ia1(j) == i ) then
dg = abs(am3%aspk(j))
end if
end do
anorm = max(anorm,tmp/dg)
enddo
call psb_amx(ictxt,anorm)
if (toupper(am3%fida)=='CSR') then
anorm = dzero
dg = done
do i=1,am3%m
tmp = dzero
do j=am3%ia2(i),am3%ia2(i+1)-1
if (am3%ia1(j) <= am3%m) then
tmp = tmp + abs(am3%aspk(j))
endif
if (am3%ia1(j) == i ) then
dg = abs(am3%aspk(j))
end if
end do
anorm = max(anorm,tmp/dg)
enddo
call psb_amx(ictxt,anorm)
else
info = 4001
endif
else
anorm = psb_spnrmi(am3,desc_a,info)
endif
if (info /= 0) then
call psb_errpush(4001,name,a_err='Invalid AM3 storage format')
goto 9999
end if
omega = 4.d0/(3.d0*anorm)
p%dprcparm(mld_aggr_damp_) = omega
@ -311,8 +300,9 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
omega = p%dprcparm(mld_aggr_damp_)
else if (p%iprcparm(mld_aggr_eig_) /= mld_user_choice_) then
write(0,*) me,'Error: invalid choice for OMEGA in blaggrmat?? ',&
& p%iprcparm(mld_aggr_eig_)
info = 4001
call psb_errpush(info,name,a_err='invalid mld_aggr_eig_')
goto 9999
end if
@ -326,41 +316,14 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
end if
end do
end do
else if (toupper(am3%fida)=='COO') then
do j=1,am3%infoa(psb_nnz_)
if (am3%ia1(j) /= am3%ia2(j)) then
am3%aspk(j) = - omega*am3%aspk(j)
else
am3%aspk(j) = zone - omega*am3%aspk(j)
endif
end do
call psb_spcnv(am3,info,afmt='csr',dupl=psb_dupl_add_)
if (info /=0) then
call psb_errpush(4010,name,a_err='spcnv am3')
goto 9999
end if
else
write(0,*) 'Missing implementation of I sum'
call psb_errpush(4010,name)
call psb_errpush(4001,name,a_err='Invalid AM3 storage format')
goto 9999
end if
if (test_dump) then
open(30+me)
write(30+me,*) 'OMEGA: ',omega
do i=1,size(p%dorig)
write(30+me,*) p%dorig(i)
end do
close(30+me)
end if
if (test_dump) call &
& psb_csprt(20+me,am4,head='% Operator Ptilde.',ivr=desc_a%loc_to_glob)
if (test_dump) call psb_csprt(40+me,am3,head='% (I-wDA)',ivr=desc_a%loc_to_glob,&
& ivc=desc_a%loc_to_glob)
if (debug) write(0,*) me,'Done gather, going for SYMBMM 1'
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done gather, going for SYMBMM 1'
!
! Symbmm90 does the allocation for its result.
!
@ -376,7 +339,9 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
call psb_numbmm(am3,am4,am1)
if (debug) write(0,*) me,'Done NUMBMM 1'
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
call psb_sp_free(am4,info)
if(info /= 0) then
@ -391,35 +356,15 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
!
call psb_sphalo(am1,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_sphalo')
goto 9999
end if
call psb_rwextd(ncol,am1,info,b=am4)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_rwextd')
goto 9999
end if
call psb_sp_free(am4,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
goto 9999
end if
if (info == 0) call psb_rwextd(ncol,am1,info,b=am4)
if (info == 0) call psb_sp_free(am4,info)
else
call psb_rwextd(ncol,am1,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='rwextd')
goto 9999
end if
endif
if (test_dump) &
& call psb_csprt(60+me,am1,head='% (I-wDA)Pt',ivr=desc_a%loc_to_glob)
if(info /= 0) then
call psb_errpush(4001,name,a_err='Halo of am1')
goto 9999
end if
call psb_symbmm(a,am1,am3,info)
if(info /= 0) then
@ -428,7 +373,9 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
end if
call psb_numbmm(a,am1,am3)
if (debug) write(0,*) me,'Done NUMBMM 2'
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2'
if (p%iprcparm(mld_aggr_kind_) == mld_smooth_prol_) then
call psb_transp(am1,am2,fmt='COO')
@ -446,7 +393,6 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
am2%ia2(i) = am2%ia2(k)
end if
end do
am2%infoa(psb_nnz_) = i
call psb_spcnv(am2,info,afmt='csr',dupl=psb_dupl_add_)
if (info /=0) then
@ -456,63 +402,38 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
else
call psb_transp(am1,am2)
endif
if (debug) write(0,*) me,'starting sphalo/ rwxtd'
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd'
if (p%iprcparm(mld_aggr_kind_) == mld_smooth_prol_) then
if (p%iprcparm(mld_aggr_kind_) == mld_smooth_prol_) then
! am2 = ((i-wDA)Ptilde)^T
call psb_sphalo(am3,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_sphalo')
goto 9999
end if
call psb_rwextd(ncol,am3,info,b=am4)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_rwextd')
goto 9999
end if
call psb_sp_free(am4,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
goto 9999
end if
if (info == 0) call psb_rwextd(ncol,am3,info,b=am4)
if (info == 0) call psb_sp_free(am4,info)
else if (p%iprcparm(mld_aggr_kind_) == mld_biz_prol_) then
call psb_rwextd(ncol,am3,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_rwextd')
goto 9999
end if
endif
if (debug) write(0,*) me,'starting symbmm 3'
call psb_symbmm(am2,am3,b,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='symbmm 3')
call psb_errpush(4001,name,a_err='Extend am3')
goto 9999
end if
if (debug) write(0,*) me,'starting numbmm 3'
call psb_numbmm(am2,am3,b)
if (debug) write(0,*) me,'Done NUMBMM 3'
!!$ if (aggr_dump) call csprt(50+me,am1,head='% Operator PTrans.')
call psb_sp_free(am3,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
goto 9999
end if
call psb_spcnv(b,info,afmt='coo',dupl=psb_dupl_add_)
if (info /=0) then
call psb_errpush(4010,name,a_err='spcnv b')
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting symbmm 3'
call psb_symbmm(am2,am3,b,info)
if (info == 0) call psb_numbmm(am2,am3,b)
if (info == 0) call psb_sp_free(am3,info)
if (info == 0) call psb_spcnv(b,info,afmt='coo',dupl=psb_dupl_add_)
if (info /= 0) then
call psb_errpush(4001,name,a_err='Build b = am2 x am3')
goto 9999
end if
if (test_dump) call psb_csprt(80+me,b,head='% Smoothed aggregate AC.')
select case(p%iprcparm(mld_aggr_kind_))
@ -523,55 +444,30 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
case(mld_distr_mat_)
call psb_sp_clone(b,ac,info)
if(info /= 0) goto 9999
nzac = ac%infoa(psb_nnz_)
nzl = ac%infoa(psb_nnz_)
call psb_cdall(ictxt,desc_ac,info,nl=p%nlaggr(me+1))
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdall')
goto 9999
end if
call psb_cdins(nzl,ac%ia1,ac%ia2,desc_ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdins')
goto 9999
end if
if (debug) write(0,*) me,'Created aux descr. distr.'
call psb_cdasb(desc_ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdasb')
goto 9999
end if
if (debug) write(0,*) me,'Asmbld aux descr. distr.'
call psb_glob_to_loc(ac%ia1(1:nzl),desc_ac,info,iact='I')
if(info /= 0) then
call psb_errpush(4010,name,a_err='psglob_to_loc')
goto 9999
end if
call psb_glob_to_loc(ac%ia2(1:nzl),desc_ac,info,iact='I')
if(info /= 0) then
call psb_errpush(4010,name,a_err='psglob_to_loc')
if (info == 0) call psb_cdall(ictxt,desc_ac,info,nl=p%nlaggr(me+1))
if (info == 0) call psb_cdins(nzl,ac%ia1,ac%ia2,desc_ac,info)
if (info == 0) call psb_cdasb(desc_ac,info)
if (info == 0) call psb_glob_to_loc(ac%ia1(1:nzl),desc_ac,info,iact='I')
if (info == 0) call psb_glob_to_loc(ac%ia2(1:nzl),desc_ac,info,iact='I')
if (info /= 0) then
call psb_errpush(4001,name,a_err='Creating desc_ac and converting ac')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
ac%m=desc_ac%matrix_data(psb_n_row_)
ac%k=desc_ac%matrix_data(psb_n_col_)
ac%fida='COO'
ac%descra='G'
ac%descra='GUN'
call psb_sp_free(b,info)
if (info == 0) deallocate(nzbr,idisp,stat=info)
if(info /= 0) then
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
goto 9999
end if
@ -579,7 +475,6 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
if (np>1) then
nzl = psb_sp_get_nnzeros(am1)
call psb_glob_to_loc(am1%ia1(1:nzl),desc_ac,info,'I')
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_glob_to_loc')
goto 9999
@ -589,44 +484,31 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
if (np>1) then
call psb_spcnv(am2,info,afmt='coo',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spcnv')
goto 9999
end if
nzl = am2%infoa(psb_nnz_)
call psb_glob_to_loc(am2%ia1(1:nzl),desc_ac,info,'I')
if (info == 0) call psb_glob_to_loc(am2%ia1(1:nzl),desc_ac,info,'I')
if (info == 0) call psb_spcnv(am2,info,afmt='csr',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_glob_to_loc')
goto 9999
end if
call psb_spcnv(am2,info,afmt='csr',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spcnv')
call psb_errpush(4001,name,a_err='Converting am2 to local')
goto 9999
end if
end if
am2%m=desc_ac%matrix_data(psb_n_col_)
if (debug) write(0,*) me,'Done ac '
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done ac '
case(mld_repl_mat_)
!
!
call psb_cdall(ictxt,desc_ac,info,mg=ntaggr,repl=.true.)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdall')
goto 9999
end if
nzbr(:) = 0
nzbr(me+1) = b%infoa(psb_nnz_)
call psb_sum(ictxt,nzbr(1:np))
nzac = sum(nzbr)
call psb_sp_all(ntaggr,ntaggr,ac,nzac,info)
if(info /= 0) goto 9999
if (info == 0) call psb_sp_all(ntaggr,ntaggr,ac,nzac,info)
if (info /= 0) goto 9999
do ip=1,np
idisp(ip) = sum(nzbr(1:ip-1))
@ -635,30 +517,36 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
call mpi_allgatherv(b%aspk,ndx,mpi_double_complex,ac%aspk,nzbr,idisp,&
& mpi_double_complex,icomm,info)
call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,&
if (info == 0) call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,&
& mpi_integer,icomm,info)
call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,&
if (info == 0) call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,&
& mpi_integer,icomm,info)
if(info /= 0) goto 9999
if (info /= 0) then
call psb_errpush(4001,name,a_err=' from mpi_allgatherv')
goto 9999
end if
ac%m = ntaggr
ac%k = ntaggr
ac%infoa(psb_nnz_) = nzac
ac%fida='COO'
ac%descra='G'
ac%descra='GUN'
call psb_spcnv(ac,info,afmt='coo',dupl=psb_dupl_add_)
if(info /= 0) goto 9999
call psb_sp_free(b,info)
if(info /= 0) goto 9999
if (me==0) then
if (test_dump) call psb_csprt(80+me,ac,head='% Smoothed aggregate AC.')
endif
deallocate(nzbr,idisp)
deallocate(nzbr,idisp,stat=info)
if (info /= 0) then
info = 4000
call psb_errpush(info,name)
goto 9999
end if
case default
write(0,*) 'Inconsistent input in smooth_new_aggregate'
info = 4001
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
@ -669,25 +557,11 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
case(mld_distr_mat_)
call psb_sp_clone(b,ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spclone')
goto 9999
end if
call psb_cdall(ictxt,desc_ac,info,nl=naggr)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdall')
goto 9999
end if
call psb_cdasb(desc_ac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdasb')
goto 9999
end if
call psb_sp_free(b,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='sp_free')
if (info == 0) call psb_cdall(ictxt,desc_ac,info,nl=naggr)
if (info == 0) call psb_cdasb(desc_ac,info)
if (info == 0) call psb_sp_free(b,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Build desc_ac, ac')
goto 9999
end if
@ -718,13 +592,12 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
call mpi_allgatherv(b%aspk,ndx,mpi_double_complex,ac%aspk,nzbr,idisp,&
& mpi_double_complex,icomm,info)
call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,&
if (info == 0) call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,&
& mpi_integer,icomm,info)
call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,&
if (info == 0) call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,&
& mpi_integer,icomm,info)
if(info /= 0) then
info=-1
call psb_errpush(info,name)
if (info /= 0) then
call psb_errpush(4001,name,a_err=' from mpi_allgatherv')
goto 9999
end if
@ -733,7 +606,7 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
ac%k = ntaggr
ac%infoa(psb_nnz_) = nzac
ac%fida='COO'
ac%descra='G'
ac%descra='GUN'
call psb_spcnv(ac,info,afmt='coo',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spcnv')
@ -745,8 +618,23 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
goto 9999
end if
case default
info = 4001
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
deallocate(nzbr,idisp)
deallocate(nzbr,idisp,stat=info)
if (info /= 0) then
info = 4000
call psb_errpush(info,name)
goto 9999
end if
case default
info = 4001
call psb_errpush(info,name,a_err='invalid mld_smooth_prol_')
goto 9999
end select
@ -756,7 +644,9 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
goto 9999
end if
if (debug) write(0,*) me,'Done smooth_aggregate '
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done smooth_aggregate '
call psb_erractionrestore(err_act)
return

@ -78,7 +78,7 @@ subroutine mld_zasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
Implicit None
! Arguments
! Arguments
integer, intent(in) :: ptype,novr
Type(psb_zspmat_type), Intent(in) :: a
Type(psb_zspmat_type), Intent(inout) :: blk
@ -88,21 +88,24 @@ subroutine mld_zasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
Character, Intent(in) :: upd
character(len=5), optional :: outfmt
! Local variables
real(kind(1.d0)) :: t1,t2,t3
! Local variables
integer icomm
Integer :: np,me,nnzero,&
& ictxt, n_col,int_err(5),&
& tot_recv, n_row,nhalo, nrow_a,err_act
Logical,Parameter :: debug=.false., debugprt=.false.
integer :: debug_level, debug_unit
character(len=20) :: name, ch_err
name='mld_zasmat_bld'
if(psb_get_errstatus().ne.0) return
if(psb_get_errstatus() /= 0) return
info=0
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
If(debug) Write(0,*)'IN DASMATBLD ', upd
If (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' start ', upd
ictxt = psb_cd_get_context(desc_data)
icomm = psb_cd_get_mpic(desc_data)
@ -121,7 +124,6 @@ subroutine mld_zasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
! Block-Jacobi preconditioner. Copy the descriptor, just in case
! we want to renumber the rows and columns of the matrix.
!
If(debug) Write(0,*)' asmatbld calling allocate '
call psb_sp_all(0,0,blk,1,info)
if(info /= 0) then
info=4010
@ -131,10 +133,11 @@ subroutine mld_zasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
end if
blk%fida = 'COO'
blk%infoa(psb_nnz_) = 0
If(debug) Write(0,*)' asmatbld done spallocate'
If (upd == 'F') Then
call psb_cdcpy(desc_data,desc_p,info)
If(debug) Write(0,*)' asmatbld done cdcpy'
If(debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' done cdcpy'
if(info /= 0) then
info=4010
ch_err='psb_cdcpy'
@ -144,12 +147,9 @@ subroutine mld_zasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
endif
case(mld_as_)
!
! Additive Schwarz
!
if (novr < 0) then
info=3
int_err(1)=novr
@ -161,7 +161,9 @@ subroutine mld_zasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
!
! Actually, this is just block Jacobi
!
If(debug) Write(0,*)' asmatbld calling allocate novr=0'
If(debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' calling allocate novr=0'
call psb_sp_all(0,0,blk,1,info)
if(info /= 0) then
info=4010
@ -171,25 +173,28 @@ subroutine mld_zasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
end if
blk%fida = 'COO'
blk%infoa(psb_nnz_) = 0
if (debug) write(0,*) 'Calling desccpy'
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Calling desccpy'
if (upd == 'F') then
call psb_cdcpy(desc_data,desc_p,info)
If(debug) Write(0,*)' asmatbld done cdcpy'
If(debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' done cdcpy'
if(info /= 0) then
info=4010
ch_err='psb_cdcpy'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (debug) write(0,*) 'Early return from asmatbld: P>=3 N_OVR=0'
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Early return: P>=3 N_OVR=0'
endif
return
endif
If(debug)Write(0,*)'BEGIN dasmatbld',me,upd,novr
t1 = psb_wtime()
If (upd == 'F') Then
!
! Build the auxiliary descriptor desc_p%matrix_data(psb_n_row_).
@ -200,7 +205,7 @@ subroutine mld_zasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
! a descriptor for an extended stencil in a PDE solver.
!
call psb_cdbldext(a,desc_data,novr,desc_p,info,extype=psb_ovt_asov_)
if(info /= 0) then
if (info /= 0) then
info=4010
ch_err='psb_cdbldext'
call psb_errpush(info,name,a_err=ch_err)
@ -208,7 +213,9 @@ subroutine mld_zasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
end if
Endif
if(debug) write(0,*) me,' From cdbldext _:',desc_p%matrix_data(psb_n_row_),&
if(debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' From cdbldext _:',desc_p%matrix_data(psb_n_row_),&
& desc_p%matrix_data(psb_n_col_)
!
@ -216,30 +223,35 @@ subroutine mld_zasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
!
n_row = desc_p%matrix_data(psb_n_row_)
t2 = psb_wtime()
if (debug) write(0,*) 'Before sphalo ',blk%fida,blk%m,psb_nnz_,blk%infoa(psb_nnz_)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Before sphalo ',blk%fida,blk%m,psb_nnz_,blk%infoa(psb_nnz_)
Call psb_sphalo(a,desc_p,blk,info,&
& outfmt=outfmt,data=psb_comm_ext_,rowscale=.true.)
if(info /= 0) then
if (info /= 0) then
info=4010
ch_err='psb_sphalo'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (debug) write(0,*) 'After psb_sphalo ',&
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'After psb_sphalo ',&
& blk%fida,blk%m,psb_nnz_,blk%infoa(psb_nnz_)
case default
if(info /= 0) then
info=4000
info=4001
ch_err='Invalid ptype'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
End select
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),'Done'
call psb_erractionrestore(err_act)
return

@ -97,8 +97,7 @@ subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
integer :: n_row,n_col, int_err(5), nrow_d
complex(kind(1.d0)), pointer :: ww(:), aux(:), tx(:),ty(:)
character ::diagl, diagu
integer :: ictxt,np,me, isz, err_act
logical,parameter :: debug=.false., debugprt=.false.
integer :: ictxt,np,me,isz, err_act
character(len=20) :: name, ch_err
name='mld_zbaseprec_aply'
@ -116,10 +115,10 @@ subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
case('N','n')
case('T','t','C','c')
case default
info=40
int_err(1)=6
ch_err(2:2)=trans
goto 9999
info=40
int_err(1)=6
ch_err(2:2)=trans
goto 9999
end select
select case(prec%iprcparm(mld_prec_type_))
@ -164,7 +163,7 @@ subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
!
call mld_bjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='mld_bjac_aply'
goto 9999
@ -180,7 +179,7 @@ subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
! shortcut: this fixes performance for RAS(0) == BJA
!
call mld_bjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='psb_bjacaply'
goto 9999
@ -230,9 +229,6 @@ subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
endif
if (debugprt) write(0,*)' vdiag: ',prec%d(:)
if (debug) write(0,*) 'Bi-CGSTAB with Additive Schwarz prec'
tx(1:nrow_d) = x(1:nrow_d)
tx(nrow_d+1:isz) = zzero
@ -247,8 +243,8 @@ subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
goto 9999
end if
else if (prec%iprcparm(mld_sub_restr_) /= psb_none_) then
write(0,*) 'Problem in PREC_APLY: Unknown value for restriction ',&
&prec%iprcparm(mld_sub_restr_)
call psb_errpush(4001,name,a_err='Invalid mld_sub_restr_')
goto 9999
end if
!
@ -270,7 +266,7 @@ subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
! preconditioner). The resulting vector is ty.
!
call mld_bjac_aply(zone,prec,tx,zzero,ty,prec%desc_data,trans,aux,info)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='mld_bjac_aply'
goto 9999
@ -309,8 +305,8 @@ subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
end if
case default
write(0,*) 'Problem in PREC_APLY: Unknown value for prolongation ',&
& prec%iprcparm(mld_sub_prol_)
call psb_errpush(4001,name,a_err='Invalid mld_sub_prol_')
goto 9999
end select
!
@ -330,9 +326,8 @@ subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
end if
case default
write(0,*) 'Invalid PRE%PREC ',prec%iprcparm(mld_prec_type_),':',&
& mld_min_prec_,mld_noprec_,mld_diag_,mld_bjac_,mld_as_
call psb_errpush(4001,name,a_err='Invalid mld_prec_type_')
goto 9999
end select
call psb_erractionrestore(err_act)
@ -341,7 +336,7 @@ subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
9999 continue
call psb_errpush(info,name,i_err=int_err,a_err=ch_err)
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if

@ -77,34 +77,33 @@ subroutine mld_zbaseprc_bld(a,desc_a,p,info,upd)
! Local variables
Integer :: err, n_row, n_col,ictxt, me,np,mglob, err_act
integer :: int_err(5)
character :: iupd
logical, parameter :: debug=.false.
integer,parameter :: iroot=0,iout=60,ilout=40
integer :: debug_level, debug_unit
character(len=20) :: name, ch_err
if(psb_get_errstatus().ne.0) return
if (psb_get_errstatus() /= 0) return
name = 'mld_zbaseprc_bld'
info=0
err=0
call psb_erractionsave(err_act)
name = 'mld_zbaseprc_bld'
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
if (debug) write(0,*) 'Entering baseprc_bld'
info = 0
int_err(1) = 0
ictxt = psb_cd_get_context(desc_a)
n_row = psb_cd_get_local_rows(desc_a)
n_col = psb_cd_get_local_cols(desc_a)
mglob = psb_cd_get_global_rows(desc_a)
if (debug) write(0,*) 'Preconditioner Blacs_gridinfo'
call psb_info(ictxt, me, np)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' start'
if (present(upd)) then
if (debug) write(0,*) 'UPD ', upd
if ((UPD.eq.'F').or.(UPD.eq.'T')) then
IUPD=UPD
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),'UPD ', upd
if ((toupper(UPD) == 'F').or.(toupper(UPD) == 'T')) then
IUPD=toupper(UPD)
else
IUPD='F'
endif
@ -140,7 +139,9 @@ subroutine mld_zbaseprc_bld(a,desc_a,p,info,upd)
! Diagonal preconditioner
call mld_diag_bld(a,desc_a,p,iupd,info)
if(debug) write(0,*)me,': out of mld_diag_bld'
if(debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ': out of mld_diag_bld'
if(info /= 0) then
info=4010
ch_err='mld_diag_bld'
@ -168,8 +169,9 @@ subroutine mld_zbaseprc_bld(a,desc_a,p,info,upd)
p%iprcparm(mld_smooth_sweeps_) = 1
end if
if (debug) write(0,*)me, ': Calling mld_bjac_bld'
if (debug) call psb_barrier(ictxt)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ': Calling mld_bjac_bld'
! Build the local part of the base preconditioner
call mld_bjac_bld(a,desc_a,p,iupd,info)
@ -180,7 +182,7 @@ subroutine mld_zbaseprc_bld(a,desc_a,p,info,upd)
end if
case default
info=4010
info=4001
ch_err='Unknown mld_prec_type_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
@ -190,6 +192,8 @@ subroutine mld_zbaseprc_bld(a,desc_a,p,info,upd)
p%base_a => a
p%base_desc => desc_a
p%iprcparm(mld_prec_status_) = mld_prec_built_
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),': Done'
call psb_erractionrestore(err_act)
return

@ -150,8 +150,7 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
! Local variables
integer :: n_row,n_col
complex(kind(1.d0)), pointer :: ww(:), aux(:), tx(:),ty(:)
integer :: ictxt,np,me,i, err_act, int_err(5)
logical,parameter :: debug=.false., debugprt=.false.
integer :: ictxt,np,me,i, err_act
character(len=20) :: name
interface
@ -206,10 +205,6 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
end if
endif
if (debug) then
write(0,*) me,' mld_bjac_APLY: ',prec%iprcparm(mld_sub_solve_),prec%iprcparm(mld_smooth_sweeps_)
end if
if (prec%iprcparm(mld_smooth_sweeps_) == 1) then
!
! TASKS 1, 3 and 4
@ -228,19 +223,17 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
call psb_spsm(zone,prec%av(mld_l_pr_),x,zzero,ww,desc_data,info,&
& trans='N',unit='L',diag=prec%d,choice=psb_none_,work=aux)
if(info /=0) goto 9999
call psb_spsm(alpha,prec%av(mld_u_pr_),ww,beta,y,desc_data,info,&
if (info == 0) call psb_spsm(alpha,prec%av(mld_u_pr_),ww,beta,y,desc_data,info,&
& trans='N',unit='U',choice=psb_none_, work=aux)
if(info /=0) goto 9999
case('T','C')
call psb_spsm(zone,prec%av(mld_u_pr_),x,zzero,ww,desc_data,info,&
& trans=trans,unit='L',diag=prec%d,choice=psb_none_, work=aux)
if(info /=0) goto 9999
call psb_spsm(alpha,prec%av(mld_l_pr_),ww,beta,y,desc_data,info,&
if(info ==0) call psb_spsm(alpha,prec%av(mld_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans,unit='U',choice=psb_none_,work=aux)
if(info /=0) goto 9999
case default
call psb_errpush(4001,name,a_err='Invalid TRANS in ILU subsolve')
goto 9999
end select
case(mld_slu_)
@ -260,10 +253,12 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
call mld_zslu_solve(1,n_row,1,ww,n_row,prec%iprcparm(mld_slu_ptr_),info)
case('C')
call mld_zslu_solve(2,n_row,1,ww,n_row,prec%iprcparm(mld_slu_ptr_),info)
case default
call psb_errpush(4001,name,a_err='Invalid TRANS in SLU subsolve')
goto 9999
end select
if(info /=0) goto 9999
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
if (info ==0) call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
case(mld_sludist_)
!
@ -280,10 +275,12 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
call mld_zsludist_solve(1,n_row,1,ww,n_row,prec%iprcparm(mld_slud_ptr_),info)
case('C')
call mld_zsludist_solve(2,n_row,1,ww,n_row,prec%iprcparm(mld_slud_ptr_),info)
case default
call psb_errpush(4001,name,a_err='Invalid TRANS in SLUDist subsolve')
goto 9999
end select
if(info /=0) goto 9999
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
if (info == 0) call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
case (mld_umf_)
!
@ -301,16 +298,22 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
call mld_zumf_solve(1,n_row,ww,x,n_row,prec%iprcparm(mld_umf_numptr_),info)
case('C')
call mld_zumf_solve(2,n_row,ww,x,n_row,prec%iprcparm(mld_umf_numptr_),info)
case default
call psb_errpush(4001,name,a_err='Invalid TRANS in UMF subsolve')
goto 9999
end select
if(info /=0) goto 9999
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
if (info == 0) call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
case default
write(0,*) 'Unknown factorization type in mld_bjac_aply',prec%iprcparm(mld_sub_solve_)
call psb_errpush(4001,name,a_err='Invalid mld_sub_solve_')
goto 9999
end select
if (debugprt) write(0,*)' Y: ',y(:)
if (info /= 0) then
call psb_errpush(4001,name,a_err='Error in subsolve Jacobi Sweeps = 1')
goto 9999
endif
else if (prec%iprcparm(mld_smooth_sweeps_) > 1) then
@ -352,15 +355,15 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
ty(1:n_row) = x(1:n_row)
call psb_spmm(-zone,prec%av(mld_ap_nd_),tx,zone,ty,&
& prec%desc_data,info,work=aux)
if(info /=0) goto 9999
if (info /=0) exit
call psb_spsm(zone,prec%av(mld_l_pr_),ty,zzero,ww,&
& prec%desc_data,info,&
& trans='N',unit='L',diag=prec%d,choice=psb_none_,work=aux)
if(info /=0) goto 9999
if (info /=0) exit
call psb_spsm(zone,prec%av(mld_u_pr_),ww,zzero,tx,&
& prec%desc_data,info,&
& trans='N',unit='U',choice=psb_none_,work=aux)
if(info /=0) goto 9999
if (info /=0) exit
end do
case(mld_sludist_)
@ -385,10 +388,10 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
ty(1:n_row) = x(1:n_row)
call psb_spmm(-zone,prec%av(mld_ap_nd_),tx,zone,ty,&
& prec%desc_data,info,work=aux)
if(info /=0) goto 9999
if (info /=0) exit
call mld_zslu_solve(0,n_row,1,ty,n_row,prec%iprcparm(mld_slu_ptr_),info)
if(info /=0) goto 9999
if (info /=0) exit
tx(1:n_row) = ty(1:n_row)
end do
@ -406,23 +409,34 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
ty(1:n_row) = x(1:n_row)
call psb_spmm(-zone,prec%av(mld_ap_nd_),tx,zone,ty,&
& prec%desc_data,info,work=aux)
if(info /=0) goto 9999
if (info /=0) exit
call mld_zumf_solve(0,n_row,ww,ty,n_row,&
& prec%iprcparm(mld_umf_numptr_),info)
if(info /=0) goto 9999
if (info /=0) exit
tx(1:n_row) = ww(1:n_row)
end do
case default
call psb_errpush(4001,name,a_err='Invalid mld_sub_solve_')
goto 9999
end select
if (info /= 0) then
info=4001
call psb_errpush(info,name,a_err='subsolve with Jacobi sweeps > 1')
goto 9999
end if
!
! Put the result into the output vector Y.
!
call psb_geaxpby(alpha,tx,beta,y,desc_data,info)
deallocate(tx,ty)
deallocate(tx,ty,stat=info)
if (info /= 0) then
info=4001
call psb_errpush(info,name,a_err='final cleanup with Jacobi sweeps > 1')
goto 9999
end if
else
@ -446,7 +460,6 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
return
9999 continue
call psb_errpush(info,name,i_err=int_err)
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()

@ -59,19 +59,14 @@
! 2. setup of block-Jacobi sweeps to compute an approximate solution of a
! linear system
! A*Y = X,
!
! distributed among the processes (allowed only at the coarsest level);
!
! 3. LU factorization of a linear system
!
! A*Y = X,
!
! distributed among the processes (allowed only at the coarsest level);
!
! 4. LU or incomplete LU factorization of a linear system
!
! A*Y = X,
!
! replicated on the processes (allowed only at the coarsest level).
!
! The following factorizations are available:
@ -116,7 +111,7 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
integer :: int_err(5)
character :: trans, unitd
type(psb_zspmat_type) :: blck, atmp
logical, parameter :: debugprt=.false., debug=.false., aggr_dump=.false.
integer :: debug_level, debug_unit
integer :: err_act, n_row, nrow_a,n_col
integer :: ictxt,np,me
character(len=20) :: name
@ -126,8 +121,9 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
info=0
name='mld_zbjac_bld'
call psb_erractionsave(err_act)
ictxt=psb_cd_get_context(desc_a)
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)
m = a%m
@ -152,9 +148,9 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
call psb_nullify_sp(atmp)
if(debug) write(0,*)me,': calling mld_asmat_bld',&
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),': Start',&
& p%iprcparm(mld_prec_type_),p%iprcparm(mld_n_ovr_)
if (debug) call psb_barrier(ictxt)
!
! Build the communication descriptor for the Additive Schwarz
@ -166,22 +162,14 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
call mld_asmat_bld(p%iprcparm(mld_prec_type_),p%iprcparm(mld_n_ovr_),a,&
& blck,desc_a,upd,p%desc_data,info,outfmt=csrfmt)
if (debugprt) then
open(60+me)
call psb_csprt(60+me,a,head='% A')
close(60+me)
open(70+me)
call psb_csprt(70+me,blck,head='% BLCK')
close(70+me)
endif
if(info/=0) then
if (info/=0) then
call psb_errpush(4010,name,a_err='mld_asmat_bld')
goto 9999
end if
if (debug) write(0,*)me,': out of mld_asmat_bld'
if (debug) call psb_barrier(ictxt)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ': out of mld_asmat_bld'
!
! Treat separately the case the local matrix has to be reordered
@ -200,7 +188,6 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
! matrix is stored into atmp, using the COO format.
!
call mld_sp_renum(a,desc_a,blck,p,atmp,info)
if (info/=0) then
call psb_errpush(4010,name,a_err='mld_sp_renum')
goto 9999
@ -212,10 +199,10 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
!
call psb_sp_clip(atmp,p%av(mld_ap_nd_),info,&
& jmin=atmp%m+1,rscale=.false.,cscale=.false.)
call psb_spcnv(p%av(mld_ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv csr 1')
if (info == 0) call psb_spcnv(p%av(mld_ap_nd_),info,&
& afmt='csr',dupl=psb_dupl_add_)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv')
goto 9999
end if
@ -232,14 +219,9 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
end if
if (debugprt) then
call psb_barrier(ictxt)
open(40+me)
call psb_csprt(40+me,atmp,head='% Local matrix')
close(40+me)
endif
if (debug) write(0,*) me,' Factoring rows ',&
&atmp%m,a%m,blck%m,atmp%ia2(atmp%m+1)-1
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' Factoring rows ',&
& atmp%m,a%m,blck%m,atmp%ia2(atmp%m+1)-1
!
! Compute a factorization of the diagonal block of the local matrix,
@ -248,90 +230,56 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
select case(p%iprcparm(mld_sub_solve_))
case(mld_ilu_n_,mld_milu_n_,mld_ilu_t_)
!
! ILU(k)/MILU(k)/ILU(k,t) factorization.
!
!
! ILU(k)/MILU(k)/ILU(k,t) factorization.
!
call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv csr 2')
goto 9999
end if
call mld_ilu_bld(atmp,p%desc_data,p,upd,info)
if (info == 0) call mld_ilu_bld(atmp,p%desc_data,p,upd,info)
if (info/=0) then
call psb_errpush(4010,name,a_err='mld_ilu_bld')
goto 9999
end if
if (debugprt) then
open(80+me)
call psb_csprt(80+me,p%av(mld_l_pr_),head='% Local L factor')
write(80+me,*) '% Diagonal: ',p%av(mld_l_pr_)%m
do i=1,p%av(mld_l_pr_)%m
write(80+me,*) i,i,p%d(i)
enddo
call psb_csprt(80+me,p%av(mld_u_pr_),head='% Local U factor')
close(80+me)
endif
case(mld_slu_)
!
! LU factorization through the SuperLU package.
!
!
! LU factorization through the SuperLU package.
!
call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_)
if (info == 0) call mld_slu_bld(atmp,p%desc_data,p,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv csr 3')
goto 9999
end if
call mld_slu_bld(atmp,p%desc_data,p,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='mld_slu_bld')
goto 9999
end if
case(mld_umf_)
!
! LU factorization through the UMFPACK package.
!
!
! LU factorization through the UMFPACK package.
!
call psb_spcnv(atmp,info,afmt='csc',dupl=psb_dupl_add_)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv csc')
goto 9999
end if
call mld_umf_bld(atmp,p%desc_data,p,info)
if(debug) write(0,*)me,': Done mld_umf_bld ',info
if (info == 0) call mld_umf_bld(atmp,p%desc_data,p,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='mld_umf_bld')
goto 9999
end if
case(mld_f_none_)
!
! Error: no factorization required.
!
info=4010
!
! Error: no factorization required.
!
info=4001
call psb_errpush(info,name,a_err='Inconsistent prec mld_f_none_')
goto 9999
case default
info=4010
info=4001
call psb_errpush(info,name,a_err='Unknown mld_sub_solve_')
goto 9999
end select
call psb_sp_free(atmp,info)
if(info/=0) then
if (info/=0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
goto 9999
end if
@ -347,11 +295,10 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
!
select case(p%iprcparm(mld_sub_solve_))
case(mld_ilu_n_,mld_milu_n_,mld_ilu_t_)
!
! ILU(k)/MILU(k)/ILU(k,t) factorization.
!
!
! ILU(k)/MILU(k)/ILU(k,t) factorization.
!
!
! In case of multiple block-Jacobi sweeps, clip into p%av(ap_nd_)
@ -367,13 +314,13 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
! given that the output from CLIP is in COO.
call psb_sp_clip(a,p%av(mld_ap_nd_),info,&
& jmin=nrow_a+1,rscale=.false.,cscale=.false.)
call psb_sp_clip(blck,atmp,info,&
if (info == 0) call psb_sp_clip(blck,atmp,info,&
& jmin=nrow_a+1,rscale=.false.,cscale=.false.)
call psb_rwextd(n_row,p%av(mld_ap_nd_),info,b=atmp)
call psb_spcnv(p%av(mld_ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv csr 4')
if (info == 0) call psb_rwextd(n_row,p%av(mld_ap_nd_),info,b=atmp)
if (info == 0) call psb_spcnv(p%av(mld_ap_nd_),info,&
& afmt='csr',dupl=psb_dupl_add_)
if (info /= 0) then
call psb_errpush(4010,name,a_err='clip & psb_spcnv csr 4')
goto 9999
end if
@ -390,45 +337,23 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
end if
call psb_sp_free(atmp,info)
end if
!
! Compute the incomplete LU factorization.
!
call mld_ilu_bld(a,desc_a,p,upd,info,blck=blck)
if(info/=0) then
if (info/=0) then
call psb_errpush(4010,name,a_err='mld_ilu_bld')
goto 9999
end if
if (debugprt) then
open(80+me)
call psb_csprt(80+me,p%av(mld_l_pr_),head='% Local L factor')
write(80+me,*) '% Diagonal: ',p%av(mld_l_pr_)%m
do i=1,p%av(mld_l_pr_)%m
write(80+me,*) i,i,p%d(i)
enddo
call psb_csprt(80+me,p%av(mld_u_pr_),head='% Local U factor')
close(80+me)
endif
case(mld_slu_)
!
! LU factorization through the SuperLU package.
!
call psb_spcnv(a,atmp,info,afmt='coo')
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv')
goto 9999
end if
!
! LU factorization through the SuperLU package.
!
n_row = psb_cd_get_local_rows(p%desc_data)
n_col = psb_cd_get_local_cols(p%desc_data)
call psb_rwextd(n_row,atmp,info,b=blck)
call psb_spcnv(a,atmp,info,afmt='coo')
if (info == 0) call psb_rwextd(n_row,atmp,info,b=blck)
!
! In case of multiple block-Jacobi sweeps, clip into p%av(ap_nd_)
@ -437,11 +362,12 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
!
if (p%iprcparm(mld_smooth_sweeps_) > 1) then
call psb_sp_clip(atmp,p%av(mld_ap_nd_),info,&
if (info == 0) call psb_sp_clip(atmp,p%av(mld_ap_nd_),info,&
& jmin=atmp%m+1,rscale=.false.,cscale=.false.)
call psb_spcnv(p%av(mld_ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
if(info /= 0) then
if (info == 0) call psb_spcnv(p%av(mld_ap_nd_),info,&
& afmt='csr',dupl=psb_dupl_add_)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv csr 6')
goto 9999
end if
@ -458,19 +384,18 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
p%iprcparm(mld_smooth_sweeps_) = 1
end if
endif
!
! Compute the LU factorization.
!
if (info == 0) call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_)
if (info == 0) call mld_slu_bld(atmp,p%desc_data,p,info)
if(info /= 0) then
if (info /= 0) then
call psb_errpush(4010,name,a_err='mld_slu_bld')
goto 9999
end if
call psb_sp_free(atmp,info)
if(info/=0) then
if (info/=0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
goto 9999
end if
@ -481,23 +406,15 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
! when the matrix is distributed among the processes.
! NOTE: Should have NO overlap here!!!!
!
call psb_spcnv(a,atmp,info,afmt='csr')
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv')
goto 9999
end if
n_row = psb_cd_get_local_rows(p%desc_data)
n_col = psb_cd_get_local_cols(p%desc_data)
if (info == 0) call mld_sludist_bld(atmp,p%desc_data,p,info)
if(info /= 0) then
if (info /= 0) then
call psb_errpush(4010,name,a_err='mld_slu_bld')
goto 9999
end if
call psb_sp_free(atmp,info)
if(info/=0) then
if (info/=0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
goto 9999
end if
@ -527,13 +444,12 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
call psb_sp_clip(atmp,p%av(mld_ap_nd_),info,&
& jmin=atmp%m+1,rscale=.false.,cscale=.false.)
call psb_spcnv(p%av(mld_ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
if(info /= 0) then
if (info == 0) call psb_spcnv(p%av(mld_ap_nd_),info,&
& afmt='csr',dupl=psb_dupl_add_)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv csr 8')
goto 9999
end if
k = psb_sp_get_nnzeros(p%av(mld_ap_nd_))
call psb_sum(ictxt,k)
@ -551,19 +467,17 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
! Compute the LU factorization.
!
if (info == 0) call psb_ipcoo2csc(atmp,info,clshr=.true.)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_ipcoo2csc')
goto 9999
end if
call mld_umf_bld(atmp,p%desc_data,p,info)
if(debug) write(0,*)me,': Done mld_umf_bld ',info
if (info == 0) call mld_umf_bld(atmp,p%desc_data,p,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ': Done mld_umf_bld ',info
if (info /= 0) then
call psb_errpush(4010,name,a_err='mld_umf_bld')
goto 9999
end if
call psb_sp_free(atmp,info)
if(info/=0) then
if (info/=0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
goto 9999
end if
@ -573,31 +487,30 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
!
! Error: no factorization required.
!
info=4010
info=4001
call psb_errpush(info,name,a_err='Inconsistent prec mld_f_none_')
goto 9999
case default
info=4010
info=4001
call psb_errpush(info,name,a_err='Unknown mld_sub_solve_')
goto 9999
end select
case default
info=4010
info=4001
call psb_errpush(info,name,a_err='Invalid renum_')
goto 9999
end select
call psb_sp_free(blck,info)
if(info/=0) then
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
goto 9999
end if
if (debug) write(0,*) me,'End of ilu_bld'
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),'End '
call psb_erractionrestore(err_act)

@ -72,30 +72,25 @@ subroutine mld_zdiag_bld(a,desc_a,p,upd,info)
! Local variables
Integer :: err, n_row, n_col,I,j,k,ictxt,&
& me,np,mglob,lw, err_act
integer :: int_err(5)
logical, parameter :: debug=.false.
integer,parameter :: iroot=0,iout=60,ilout=40
integer :: debug_level, debug_unit
character(len=20) :: name, ch_err
if(psb_get_errstatus().ne.0) return
info=0
err=0
call psb_erractionsave(err_act)
name = 'mld_zdiag_bld'
if (debug) write(0,*) 'Entering diagsc_bld'
info = 0
int_err(1) = 0
ictxt = psb_cd_get_context(desc_a)
n_row = psb_cd_get_local_rows(desc_a)
n_col = psb_cd_get_local_cols(desc_a)
mglob = psb_cd_get_global_rows(desc_a)
if (debug) write(0,*) 'Preconditioner Blacs_gridinfo'
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
name = 'mld_zdiag_bld'
info = 0
ictxt = psb_cd_get_context(desc_a)
n_row = psb_cd_get_local_rows(desc_a)
n_col = psb_cd_get_local_cols(desc_a)
mglob = psb_cd_get_global_rows(desc_a)
call psb_info(ictxt, me, np)
if (debug) write(0,*) 'Precond: Diagonal'
if (debug_level >= psb_debug_outer_)&
& write(debug_unit,*) me,' ',trim(name),' Enter'
call psb_realloc(n_col,p%d,info)
if (info /= 0) then
@ -122,7 +117,6 @@ subroutine mld_zdiag_bld(a,desc_a,p,upd,info)
goto 9999
end if
if (debug) write(ilout+me,*) 'VDIAG ',n_row
!
! The i-th diagonal entry of the preconditioner is set to one if the
! corresponding entry a_ii of the sparse matrix A is zero; otherwise
@ -134,8 +128,6 @@ subroutine mld_zdiag_bld(a,desc_a,p,upd,info)
else
p%d(i) = zone/p%d(i)
endif
if (debug) write(ilout+me,*) i,desc_a%loc_to_glob(i), p%d(i)
end do
if (a%pl(1) /= 0) then
@ -151,8 +143,8 @@ subroutine mld_zdiag_bld(a,desc_a,p,upd,info)
end if
endif
if (debug) write(*,*) 'Preconditioner DIAG computed OK'
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),'Done'
call psb_erractionrestore(err_act)
return

@ -106,18 +106,20 @@ subroutine mld_zilu_bld(a,desc_a,p,upd,info,blck)
! Local Variables
integer :: i, nztota, err_act, n_row, nrow_a
character :: trans, unitd
logical, parameter :: debugprt=.false., debug=.false., aggr_dump=.false.
integer :: ictxt,np,me
character(len=20) :: name, ch_err
integer :: debug_level, debug_unit
integer :: ictxt,np,me
character(len=20) :: name, ch_err
if(psb_get_errstatus().ne.0) return
info=0
name='mld_zilu_bld'
call psb_erractionsave(err_act)
ictxt=psb_cd_get_context(desc_a)
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'
trans = 'N'
unitd = 'U'
@ -145,15 +147,15 @@ subroutine mld_zilu_bld(a,desc_a,p,upd,info,blck)
goto 9999
end if
endif
!!$ call psb_csprt(50+me,a,head='% (A)')
nrow_a = psb_cd_get_local_rows(desc_a)
nztota = psb_sp_get_nnzeros(a)
if (present(blck)) then
nztota = nztota + psb_sp_get_nnzeros(blck)
end if
if (debug) write(0,*)me,': out get_nnzeros',nztota,a%m,a%k
if (debug) call psb_barrier(ictxt)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ': out get_nnzeros',nztota,a%m,a%k
n_row = p%desc_data%matrix_data(psb_n_row_)
p%av(mld_l_pr_)%m = n_row
@ -253,22 +255,6 @@ subroutine mld_zilu_bld(a,desc_a,p,upd,info,blck)
end select
if (debugprt) then
!
! Print out the factors on file.
!
open(80+me)
call psb_csprt(80+me,p%av(mld_l_pr_),head='% Local L factor')
write(80+me,*) '% Diagonal: ',p%av(mld_l_pr_)%m
do i=1,p%av(mld_l_pr_)%m
write(80+me,*) i,i,p%d(i)
enddo
call psb_csprt(80+me,p%av(mld_u_pr_),head='% Local U factor')
close(80+me)
endif
if (psb_sp_getifld(psb_upd_,p%av(mld_u_pr_),info) /= psb_upd_perm_) then
call psb_sp_trim(p%av(mld_u_pr_),info)
endif
@ -277,7 +263,8 @@ subroutine mld_zilu_bld(a,desc_a,p,upd,info,blck)
call psb_sp_trim(p%av(mld_l_pr_),info)
endif
if (debug) write(0,*) me,'End of ilu_bld'
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),'End'
call psb_erractionrestore(err_act)
return

@ -115,7 +115,6 @@ subroutine mld_zilu_fct(ialg,a,l,u,d,info,blck)
integer :: l1, l2,m,err_act
type(psb_zspmat_type), pointer :: blck_
character(len=20) :: name, ch_err
logical, parameter :: debug=.false.
name='mld_zilu_fct'
info = 0
@ -289,7 +288,6 @@ contains
integer :: i,j,k,l,low1,low2,kk,jj,ll, ktrw,err_act
complex(kind(1.d0)) :: dia,temp
integer, parameter :: nrb=16
logical,parameter :: debug=.false.
type(psb_zspmat_type) :: trw
integer :: int_err(5)
character(len=20) :: name, ch_err
@ -301,7 +299,7 @@ contains
call psb_nullify_sp(trw)
trw%m=0
trw%k=0
if(debug) write(0,*)'LUINT Allocating TRW'
call psb_sp_all(trw,1,info)
if(info.ne.0) then
info=4010
@ -309,20 +307,18 @@ contains
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if(debug) write(0,*)'LUINT Done Allocating TRW'
lia2(1) = 1
uia2(1) = 1
l1 = 0
l2 = 0
m = ma+mb
if(debug) write(0,*)'In DCSRLU Begin cycle',m,ma,mb
!
! Cycle over the matrix rows
!
do i = 1, m
if(debug) write(0,*)'LUINT: Loop index ',i,ma
d(i) = zzero
if (i <= ma) then
@ -447,7 +443,6 @@ contains
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if(debug) write(0,*)'Leaving ilu_fct'
call psb_erractionrestore(err_act)
return

@ -111,14 +111,11 @@ subroutine mld_ziluk_fct(fill_in,ialg,a,l,u,d,info,blck)
type(psb_zspmat_type), pointer :: blck_
character(len=20) :: name, ch_err
logical, parameter :: debug=.false.
name='mld_ziluk_fct'
info = 0
call psb_erractionsave(err_act)
if (debug) write(0,*) 'mld_diluk_fct: start'
!
! Point to / allocate memory for the incomplete factorization
!
@ -144,7 +141,6 @@ subroutine mld_ziluk_fct(fill_in,ialg,a,l,u,d,info,blck)
!
! Compute the ILU(k) or the MILU(k) factorization, depending on ialg
!
if (debug) write(0,*) 'mld_ziluk_fct: calling fctint'
call mld_ziluk_fctint(fill_in,ialg,m,a%m,a,blck_%m,blck_,&
& d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info)
if (info /= 0) then
@ -302,7 +298,6 @@ contains
!
! Allocate a temporary buffer for the iluk_copyin function
!
if (debug) write(0,*)'LUINT Allocating TRW'
call psb_sp_all(0,0,trw,1,info)
if (info==0) call psb_ensure_size(m+1,lia2,info)
if (info==0) call psb_ensure_size(m+1,uia2,info)
@ -312,15 +307,12 @@ contains
call psb_errpush(info,name,a_err='psb_sp_all')
goto 9999
end if
if (debug) write(0,*)'LUINT Done Allocating TRW'
l1=0
l2=0
lia2(1) = 1
uia2(1) = 1
if (debug) write(0,*)'In DCSRLU Begin cycle',m,ma,mb
!
! Allocate memory to hold the entries of a row and the corresponding
! fill levels
@ -341,8 +333,6 @@ contains
!
do i = 1, m
if (debug.and.(mod(i,500)==1)) write(0,*)'LUINT: Loop index ',i,ma
!
! At each iteration of the loop we keep in a heap the column indices
! affected by the factorization. The heap is initialized and filled
@ -365,8 +355,6 @@ contains
call iluk_copyin(i-ma,mb,b,1,m,row,rowlevs,heap,ktrw,trw)
endif
if (debug) write(0,*)'LUINT: input Copy done'
! Do an elimination step on the current row. It turns out we only
! need to keep track of fill levels for the upper triangle, hence we
! do not have a lowlevs variable.
@ -397,7 +385,6 @@ contains
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (debug) write(0,*)'Leaving ilu_fct'
call psb_erractionrestore(err_act)
return

@ -110,14 +110,11 @@ subroutine mld_zilut_fct(fill_in,thres,ialg,a,l,u,d,info,blck)
type(psb_zspmat_type), pointer :: blck_
character(len=20) :: name, ch_err
logical, parameter :: debug=.false.
name='mld_zilut_fct'
info = 0
call psb_erractionsave(err_act)
if (debug) write(0,*) 'mld_zilut_fct: start'
!
! Point to / allocate memory for the incomplete factorization
!
@ -143,7 +140,6 @@ subroutine mld_zilut_fct(fill_in,thres,ialg,a,l,u,d,info,blck)
!
! Compute the ILU(k,t) factorization
!
if (debug) write(0,*) 'mld_zilut_fct: calling fctint'
call mld_zilut_fctint(fill_in,thres,ialg,m,a%m,a,blck_%m,blck_,&
& d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info)
if (info /= 0) then
@ -287,7 +283,6 @@ contains
integer, allocatable :: idxs(:)
complex(kind(1.d0)), allocatable :: row(:)
type(psb_int_heap) :: heap
logical,parameter :: debug=.false.
type(psb_zspmat_type) :: trw
character(len=20), parameter :: name='mld_zilut_fctint'
character(len=20) :: ch_err
@ -301,7 +296,6 @@ contains
!
! Allocate a temporary buffer for the ilut_copyin function
!
if (debug) write(0,*)'LUINT Allocating TRW'
call psb_sp_all(0,0,trw,1,info)
if (info==0) call psb_ensure_size(m+1,lia2,info)
if (info==0) call psb_ensure_size(m+1,uia2,info)
@ -311,15 +305,12 @@ contains
call psb_errpush(info,name,a_err='psb_sp_all')
goto 9999
end if
if (debug) write(0,*)'LUINT Done Allocating TRW'
l1=0
l2=0
lia2(1) = 1
uia2(1) = 1
if (debug) write(0,*)'In DCSRLU Begin cycle',m,ma,mb
!
! Allocate memory to hold the entries of a row
!
@ -337,7 +328,6 @@ contains
!
do i = 1, m
if (debug) write(0,*)'LUINT: Loop index ',i
!
! At each iteration of the loop we keep in a heap the column indices
! affected by the factorization. The heap is initialized and filled
@ -353,7 +343,6 @@ contains
call ilut_copyin(i-ma,mb,b,i,1,m,nlw,nup,jmaxup,nrmi,row,heap,ktrw,trw)
endif
if (debug) write(0,*)'LUINT: input Copy done'
!
! Do an elimination step on current row
!
@ -383,7 +372,6 @@ contains
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (debug) write(0,*)'Leaving ilu_fct'
call psb_erractionrestore(err_act)
return
@ -675,7 +663,6 @@ contains
! Local Variables
integer :: k,j,jj,info, lastk
complex(kind(1.d0)) :: rwk
logical, parameter :: debug=.false.
call psb_ensure_size(200,idxs,info)
@ -748,12 +735,6 @@ contains
end do
if (debug) then
write(0,*) 'At end of factint: ',i,nidx
write(0,*) idxs(1:nidx)
write(0,*) row(idxs(1:nidx))
end if
end subroutine ilut_fact
!
@ -870,7 +851,6 @@ contains
character(len=20), parameter :: name='mld_zilut_fctint'
character(len=20) :: ch_err
logical :: fndmaxup
logical, parameter :: debug=.false.
if (psb_get_errstatus() /= 0) return
info=0
@ -909,10 +889,6 @@ contains
if (idxs(idxp) >= i) exit
widx = idxs(idxp)
witem = row(widx)
if (debug) then
write(0,*) 'Lower: Deciding on drop of item ',witem,widx,thres,nrmi,thres*nrmi
end if
!
! Dropping rule based on the 2-norm
!
@ -1032,11 +1008,6 @@ contains
cycle
end if
witem = row(widx)
if (debug) then
write(0,*) 'Upper: Deciding on drop of item ',witem,widx,&
& jmaxup,thres,nrmi,thres*nrmi
end if
!
! Dropping rule based on the 2-norm. But keep the jmaxup-th entry anyway.
!
@ -1051,14 +1022,6 @@ contains
end do
if (debug) then
write(0,*) 'Row ',i,' copyout: after first round at upper:',nz,jmaxup
write(0,*) xwid(1:nz)
write(0,*) xw(1:nz)
write(0,*) 'Dumping heap'
call psb_dump_heap(0,heap,info)
end if
!
! Now we have to take out the first nup-fill_in entries. But make sure
! we include entry jmaxup.
@ -1093,11 +1056,6 @@ contains
! Now we put things back into ascending column order
!
call psb_msort(xwid(1:nz),indx(1:nz),dir=psb_sort_up_)
if (debug) then
write(0,*) 'Row ',i,' copyout: after sort at upper:',nz,jmaxup
write(0,*) xwid(1:nz)
write(0,*) xw(indx(1:nz))
end if
!
! Copy out the upper part of the row

@ -64,11 +64,11 @@
! level 1 is the finest level and A(1) is the matrix A.
!
! For a general description of (parallel) multilevel preconditioners see
! 1. B.F. Smith, P.E. Bjorstad & W.D. Gropp,
! - B.F. Smith, P.E. Bjorstad & W.D. Gropp,
! Domain decomposition: parallel multilevel methods for elliptic partial
! differential equations,
! Cambridge University Press, 1996.
! 2. K. Stuben,
! - K. Stuben,
! Algebraic Multigrid (AMG): An Introduction with Applications,
! GMD Report N. 70, 1999.
!
@ -181,10 +181,10 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
! Local variables
integer :: n_row,n_col
integer :: ictxt,np,me,i, nr2l,nc2l,err_act
logical, parameter :: debug=.false., debugprt=.false.
integer :: ismth, nlev, ilev, icm
character(len=20) :: name
integer :: ictxt,np,me,i, nr2l,nc2l,err_act
integer :: debug_level, debug_unit
integer :: ismth, nlev, ilev, icm
character(len=20) :: name
type psb_mlprec_wrk_type
complex(kind(1.d0)), allocatable :: tx(:),ty(:),x2l(:),y2l(:)
@ -194,12 +194,15 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
name='mld_zmlprec_aply'
info = 0
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = psb_cd_get_context(desc_data)
call psb_info(ictxt, me, np)
if (debug) write(0,*) me,'Entry to mlprec_aply ',&
& size(baseprecv)
if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(baseprecv)
nlev = size(baseprecv)
allocate(mlprec_wrk(nlev),stat=info)
@ -215,7 +218,7 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
! No preconditioning, should not really get here
!
call psb_errpush(4010,name,a_err='mld_no_ml_ in mlprc_aply?')
call psb_errpush(4001,name,a_err='mld_no_ml_ in mlprc_aply?')
goto 9999
@ -260,7 +263,10 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
call mld_baseprec_aply(alpha,baseprecv(1),x,beta,y,&
& baseprecv(1)%base_desc,trans,work,info)
if(info /=0) goto 9999
if (info /=0) then
call psb_errpush(4010,name,a_err='baseprec_aply')
goto 9999
end if
allocate(mlprec_wrk(1)%x2l(size(x)),mlprec_wrk(1)%y2l(size(y)), stat=info)
if (info /= 0) then
@ -309,12 +315,8 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
call psb_halo(mlprec_wrk(ilev-1)%x2l,baseprecv(ilev-1)%base_desc,&
& info,work=work)
if(info /=0) goto 9999
call psb_csmm(zone,baseprecv(ilev)%av(mld_sm_pr_t_),mlprec_wrk(ilev-1)%x2l,&
& zzero,mlprec_wrk(ilev)%x2l,info)
if(info /=0) goto 9999
if (info == 0) call psb_csmm(zone,baseprecv(ilev)%av(mld_sm_pr_t_),&
& mlprec_wrk(ilev-1)%x2l,zzero,mlprec_wrk(ilev)%x2l,info)
else
!
! Apply the raw aggregation map transpose (take a shortcut)
@ -326,11 +328,19 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
end do
end if
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction')
goto 9999
end if
if (icm == mld_repl_mat_) then
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nr2l))
else if (icm /= mld_distr_mat_) then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(mld_coarse_mat_) ',icm
info = 4013
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_',&
& i_Err=(/icm,0,0,0,0/))
goto 9999
endif
!
@ -363,8 +373,6 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
call psb_csmm(zone,baseprecv(ilev)%av(mld_sm_pr_),mlprec_wrk(ilev)%y2l,&
& zone,mlprec_wrk(ilev-1)%y2l,info)
if(info /=0) goto 9999
else
!
! Apply the raw aggregation map (take a shortcut)
@ -375,6 +383,10 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
enddo
end if
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during prolognation')
goto 9999
end if
end do
!
@ -383,8 +395,10 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
! Compute the output vector Y
!
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,zone,y,baseprecv(1)%base_desc,info)
if(info /=0) goto 9999
if (info /= 0) then
call psb_errpush(4001,name,a_err='Error on final update')
goto 9999
end if
case(mld_mult_ml_)
@ -434,8 +448,9 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
! Copy the input vector X
!
if (debug) write(0,*) me, 'mlprec_aply desc_data',&
& allocated(desc_data%matrix_data)
if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' desc_data status',allocated(desc_data%matrix_data)
n_col = psb_cd_get_local_cols(desc_data)
nc2l = psb_cd_get_local_cols(baseprecv(1)%desc_data)
@ -465,7 +480,9 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
ismth = baseprecv(ilev)%iprcparm(mld_aggr_kind_)
icm = baseprecv(ilev)%iprcparm(mld_coarse_mat_)
if (debug) write(0,*) me, 'mlprec_aply starting up sweep ',&
if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name), &
& ' starting up sweep ',&
& ilev,allocated(baseprecv(ilev)%iprcparm),n_row,n_col,&
& nc2l, nr2l,ismth
@ -481,20 +498,18 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
mlprec_wrk(ilev)%x2l(:) = zzero
mlprec_wrk(ilev)%y2l(:) = zzero
mlprec_wrk(ilev)%tx(:) = zzero
mlprec_wrk(ilev)%tx(:) = zzero
if (ismth /= mld_no_smooth_) then
!
! Apply the smoothed prolongator transpose
!
if (debug) write(0,*) me, 'mlprec_aply halo in up sweep ', ilev
if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name), ' up sweep ', ilev
call psb_halo(mlprec_wrk(ilev-1)%x2l,&
& baseprecv(ilev-1)%base_desc,info,work=work)
if(info /=0) goto 9999
if (debug) write(0,*) me, 'mlprec_aply csmm in up sweep ', ilev
call psb_csmm(zone,baseprecv(ilev)%av(mld_sm_pr_t_),mlprec_wrk(ilev-1)%x2l, &
& zzero,mlprec_wrk(ilev)%x2l,info)
if(info /=0) goto 9999
if (info == 0) call psb_csmm(zone,baseprecv(ilev)%av(mld_sm_pr_t_),&
& mlprec_wrk(ilev-1)%x2l,zzero,mlprec_wrk(ilev)%x2l,info)
else
!
! Apply the raw aggregation map transpose (take a shortcut)
@ -506,20 +521,18 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
end do
end if
if (debug) write(0,*) me, 'mlprec_aply possible sum in up sweep ', &
& ilev,icm,associated(baseprecv(ilev)%base_desc),mld_repl_mat_
if (debug) write(0,*) me, 'mlprec_aply geaxpby in up sweep X', &
& ilev,associated(baseprecv(ilev)%base_desc),&
& baseprecv(ilev)%base_desc%matrix_data(psb_n_row_),&
& baseprecv(ilev)%base_desc%matrix_data(psb_n_col_),&
& size(mlprec_wrk(ilev)%tx),size(mlprec_wrk(ilev)%x2l)
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction')
goto 9999
end if
if (icm == mld_repl_mat_) Then
if (debug) write(0,*) 'Entering psb_sum ',nr2l
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nr2l))
else if (icm /= mld_distr_mat_) Then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(mld_coarse_mat_) ', icm
info = 4013
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_',&
& i_Err=(/icm,0,0,0,0/))
goto 9999
endif
!
@ -527,9 +540,14 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
call psb_geaxpby(zone,mlprec_wrk(ilev)%x2l,zzero,mlprec_wrk(ilev)%tx,&
& baseprecv(ilev)%base_desc,info)
if(info /=0) goto 9999
if (debug) write(0,*) me, 'mlprec_aply done up sweep ',&
& ilev
if (info /= 0) then
call psb_errpush(4001,name,a_err='Error in update')
goto 9999
end if
if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' done up sweep ', ilev
enddo
@ -540,10 +558,13 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
call mld_baseprec_aply(zone,baseprecv(nlev),mlprec_wrk(nlev)%x2l, &
& zzero, mlprec_wrk(nlev)%y2l,baseprecv(nlev)%desc_data,'N',work,info)
if (info /=0) then
call psb_errpush(4010,name,a_err='baseprec_aply')
goto 9999
end if
if(info /=0) goto 9999
if (debug) write(0,*) me, 'mlprec_aply done prc_apl ',&
& nlev
if (debug_level >= psb_debug_inner_) write(debug_unit,*) &
& me,' ',trim(name), ' done baseprec_aply ', nlev
!
! STEP 4
@ -552,7 +573,10 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
do ilev=nlev-1, 1, -1
if (debug) write(0,*) me, 'mlprec_aply starting down sweep',ilev
if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' starting down sweep',ilev
ismth = baseprecv(ilev+1)%iprcparm(mld_aggr_kind_)
n_row = psb_cd_get_local_rows(baseprecv(ilev)%base_desc)
@ -563,10 +587,8 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
if (ismth == mld_smooth_prol_) &
& call psb_halo(mlprec_wrk(ilev+1)%y2l,baseprecv(ilev+1)%desc_data,&
& info,work=work)
call psb_csmm(zone,baseprecv(ilev+1)%av(mld_sm_pr_),mlprec_wrk(ilev+1)%y2l,&
& zzero,mlprec_wrk(ilev)%y2l,info)
if(info /=0) goto 9999
if (info == 0) call psb_csmm(zone,baseprecv(ilev+1)%av(mld_sm_pr_),&
& mlprec_wrk(ilev+1)%y2l, zzero,mlprec_wrk(ilev)%y2l,info)
else
!
! Apply the raw aggregation map (take a shortcut)
@ -576,7 +598,10 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
mlprec_wrk(ilev)%y2l(i) = mlprec_wrk(ilev)%y2l(i) + &
& mlprec_wrk(ilev+1)%y2l(baseprecv(ilev+1)%mlia(i))
enddo
end if
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during prolongation')
goto 9999
end if
!
@ -585,16 +610,19 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
call psb_spmm(-zone,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
& zone,mlprec_wrk(ilev)%tx,baseprecv(ilev)%base_desc,info,work=work)
if(info /=0) goto 9999
!
! Apply the base preconditioner
!
call mld_baseprec_aply(zone,baseprecv(ilev),mlprec_wrk(ilev)%tx,&
if (info == 0) call mld_baseprec_aply(zone,baseprecv(ilev),mlprec_wrk(ilev)%tx,&
& zone,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%base_desc, trans, work,info)
if (info /=0) then
call psb_errpush(4001,name,a_err=' spmm/baseprec_aply')
goto 9999
end if
if(info /=0) goto 9999
if (debug) write(0,*) me, 'mlprec_aply done down sweep',ilev
if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' done down sweep',ilev
enddo
!
@ -604,8 +632,10 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,baseprecv(1)%base_desc,info)
if(info /=0) goto 9999
if (info /=0) then
call psb_errpush(4001,name,a_err=' Final update')
goto 9999
end if
case(mld_pre_smooth_)
@ -676,9 +706,10 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
& zzero,mlprec_wrk(1)%y2l,&
& baseprecv(1)%base_desc,&
& trans,work,info)
if(info /=0) goto 9999
if (info /=0) then
call psb_errpush(4010,name,a_err=' baseprec_aply')
goto 9999
end if
!
! STEP 3
@ -689,7 +720,10 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
call psb_spmm(-zone,baseprecv(1)%base_a,mlprec_wrk(1)%y2l,&
& zone,mlprec_wrk(1)%tx,baseprecv(1)%base_desc,info,work=work)
if(info /=0) goto 9999
if (info /=0) then
call psb_errpush(4001,name,a_err=' fine level residual')
goto 9999
end if
!
! STEP 4
@ -716,8 +750,7 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
mlprec_wrk(ilev)%x2l(:) = zzero
mlprec_wrk(ilev)%y2l(:) = zzero
mlprec_wrk(ilev)%tx(:) = zzero
mlprec_wrk(ilev)%tx(:) = zzero
if (ismth /= mld_no_smooth_) then
!
@ -725,12 +758,8 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
call psb_halo(mlprec_wrk(ilev-1)%tx,baseprecv(ilev-1)%base_desc,&
& info,work=work)
if(info /=0) goto 9999
call psb_csmm(zone,baseprecv(ilev)%av(mld_sm_pr_t_),mlprec_wrk(ilev-1)%tx,zzero,&
& mlprec_wrk(ilev)%x2l,info)
if(info /=0) goto 9999
if (info == 0) call psb_csmm(zone,baseprecv(ilev)%av(mld_sm_pr_t_),&
& mlprec_wrk(ilev-1)%tx,zzero,mlprec_wrk(ilev)%x2l,info)
else
!
! Apply the raw aggregation map transpose (take a shortcut)
@ -742,11 +771,18 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
& mlprec_wrk(ilev-1)%tx(i)
end do
end if
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction')
goto 9999
end if
if (icm ==mld_repl_mat_) then
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nr2l))
else if (icm /= mld_distr_mat_) then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(mld_coarse_mat_) ', icm
info = 4013
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_',&
& i_Err=(/icm,0,0,0,0/))
goto 9999
endif
!
@ -755,18 +791,19 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
call mld_baseprec_aply(zone,baseprecv(ilev),mlprec_wrk(ilev)%x2l,&
& zzero,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%desc_data, 'N',work,info)
if(info /=0) goto 9999
!
! Compute the residual (at all levels but the coarsest one)
!
if (ilev < nlev) then
mlprec_wrk(ilev)%tx = mlprec_wrk(ilev)%x2l
call psb_spmm(-zone,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
& zone,mlprec_wrk(ilev)%tx,baseprecv(ilev)%base_desc,info,work=work)
if(info /=0) goto 9999
if (info == 0) call psb_spmm(-zone,baseprecv(ilev)%base_a,&
& mlprec_wrk(ilev)%y2l,zone,mlprec_wrk(ilev)%tx,&
& baseprecv(ilev)%base_desc,info,work=work)
endif
if (info /=0) then
call psb_errpush(4001,name,a_err='Error on up sweep residual')
goto 9999
end if
enddo
!
@ -786,11 +823,8 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
if (ismth == mld_smooth_prol_) &
& call psb_halo(mlprec_wrk(ilev+1)%y2l,&
& baseprecv(ilev+1)%desc_data,info,work=work)
call psb_csmm(zone,baseprecv(ilev+1)%av(mld_sm_pr_),mlprec_wrk(ilev+1)%y2l,&
& zone,mlprec_wrk(ilev)%y2l,info)
if(info /=0) goto 9999
if (info == 0) call psb_csmm(zone,baseprecv(ilev+1)%av(mld_sm_pr_),&
& mlprec_wrk(ilev+1)%y2l,zone,mlprec_wrk(ilev)%y2l,info)
else
!
! Apply the raw aggregation map (take a shortcut)
@ -799,9 +833,11 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
mlprec_wrk(ilev)%y2l(i) = mlprec_wrk(ilev)%y2l(i) + &
& mlprec_wrk(ilev+1)%y2l(baseprecv(ilev+1)%mlia(i))
enddo
end if
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during prolongation')
goto 9999
end if
enddo
!
@ -811,8 +847,11 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,&
& baseprecv(1)%base_desc,info)
if (info /=0) then
call psb_errpush(4001,name,a_err='Error on final update')
goto 9999
end if
if(info /=0) goto 9999
case(mld_twoside_smooth_)
@ -895,18 +934,18 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
& zzero,mlprec_wrk(1)%y2l,&
& baseprecv(1)%base_desc,&
& trans,work,info)
if(info /=0) goto 9999
!
! STEP 3
!
! Compute the residual at the finest level
!
mlprec_wrk(1)%ty = mlprec_wrk(1)%x2l
call psb_spmm(-zone,baseprecv(1)%base_a,mlprec_wrk(1)%y2l,&
if (info == 0) call psb_spmm(-zone,baseprecv(1)%base_a,mlprec_wrk(1)%y2l,&
& zone,mlprec_wrk(1)%ty,baseprecv(1)%base_desc,info,work=work)
if(info /=0) goto 9999
if (info /=0) then
call psb_errpush(4010,name,a_err='Fine level baseprec/residual')
goto 9999
end if
!
! STEP 4
@ -943,11 +982,8 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
call psb_halo(mlprec_wrk(ilev-1)%ty,baseprecv(ilev-1)%base_desc,&
& info,work=work)
if(info /=0) goto 9999
call psb_csmm(zone,baseprecv(ilev)%av(mld_sm_pr_t_),mlprec_wrk(ilev-1)%ty,zzero,&
& mlprec_wrk(ilev)%x2l,info)
if(info /=0) goto 9999
if (info == 0) call psb_csmm(zone,baseprecv(ilev)%av(mld_sm_pr_t_),&
& mlprec_wrk(ilev-1)%ty,zzero,mlprec_wrk(ilev)%x2l,info)
else
!
! Apply the raw aggregation map transpose (take a shortcut)
@ -959,34 +995,41 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
& mlprec_wrk(ilev-1)%ty(i)
end do
end if
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction')
goto 9999
end if
if (icm == mld_repl_mat_) then
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nr2l))
else if (icm /= mld_distr_mat_) then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(mld_coarse_mat_) ', icm
info = 4013
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_',&
& i_Err=(/icm,0,0,0,0/))
goto 9999
endif
call psb_geaxpby(zone,mlprec_wrk(ilev)%x2l,zzero,mlprec_wrk(ilev)%tx,&
& baseprecv(ilev)%base_desc,info)
if(info /=0) goto 9999
!
! Apply the base preconditioner
!
call mld_baseprec_aply(zone,baseprecv(ilev),mlprec_wrk(ilev)%x2l,&
& zzero,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%desc_data, 'N',work,info)
if(info /=0) goto 9999
if (info == 0) call mld_baseprec_aply(zone,baseprecv(ilev),&
& mlprec_wrk(ilev)%x2l,zzero,mlprec_wrk(ilev)%y2l,&
&baseprecv(ilev)%desc_data, 'N',work,info)
!
! Compute the residual (at all levels but the coarsest one)
!
if(ilev < nlev) then
mlprec_wrk(ilev)%ty = mlprec_wrk(ilev)%x2l
call psb_spmm(-zone,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
& zone,mlprec_wrk(ilev)%ty,baseprecv(ilev)%base_desc,info,work=work)
if(info /=0) goto 9999
if (info == 0) call psb_spmm(-zone,baseprecv(ilev)%base_a,&
& mlprec_wrk(ilev)%y2l,zone,mlprec_wrk(ilev)%ty,&
& baseprecv(ilev)%base_desc,info,work=work)
endif
if (info /=0) then
call psb_errpush(4001,name,a_err='baseprec_aply/residual')
goto 9999
end if
enddo
@ -1007,10 +1050,8 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
if (ismth == mld_smooth_prol_) &
& call psb_halo(mlprec_wrk(ilev+1)%y2l,baseprecv(ilev+1)%desc_data,&
& info,work=work)
call psb_csmm(zone,baseprecv(ilev+1)%av(mld_sm_pr_),mlprec_wrk(ilev+1)%y2l,&
& zone,mlprec_wrk(ilev)%y2l,info)
if(info /=0) goto 9999
if (info == 0) call psb_csmm(zone,baseprecv(ilev+1)%av(mld_sm_pr_),&
& mlprec_wrk(ilev+1)%y2l,zone,mlprec_wrk(ilev)%y2l,info)
else
!
! Apply the raw aggregation map (take a shortcut)
@ -1019,7 +1060,10 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
mlprec_wrk(ilev)%y2l(i) = mlprec_wrk(ilev)%y2l(i) + &
& mlprec_wrk(ilev+1)%y2l(baseprecv(ilev+1)%mlia(i))
enddo
end if
if (info /=0 ) then
call psb_errpush(4001,name,a_err='Error during restriction')
goto 9999
end if
!
@ -1027,17 +1071,15 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
call psb_spmm(-zone,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
& zone,mlprec_wrk(ilev)%tx,baseprecv(ilev)%base_desc,info,work=work)
if(info /=0) goto 9999
!
! Apply the base preconditioner
!
call mld_baseprec_aply(zone,baseprecv(ilev),mlprec_wrk(ilev)%tx,&
if (info == 0) call mld_baseprec_aply(zone,baseprecv(ilev),mlprec_wrk(ilev)%tx,&
& zone,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%base_desc, trans, work,info)
if(info /=0) goto 9999
if (info /= 0) then
call psb_errpush(4001,name,a_err='Error: residual/baseprec_aply')
goto 9999
end if
enddo
!
@ -1048,30 +1090,37 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,&
& baseprecv(1)%base_desc,info)
if(info /=0) goto 9999
if (info /= 0) then
call psb_errpush(4001,name,a_err='Error final update')
goto 9999
end if
case default
call psb_errpush(4013,name,a_err='wrong smooth_pos',&
info = 4013
call psb_errpush(info,name,a_err='invalid smooth_pos',&
& i_Err=(/baseprecv(2)%iprcparm(mld_smooth_pos_),0,0,0,0/))
goto 9999
end select
case default
call psb_errpush(4013,name,a_err='wrong mltype',&
info = 4013
call psb_errpush(info,name,a_err='invalid mltype',&
& i_Err=(/baseprecv(2)%iprcparm(mld_ml_type_),0,0,0,0/))
goto 9999
end select
deallocate(mlprec_wrk)
deallocate(mlprec_wrk,stat=info)
if (info /= 0) then
call psb_errpush(4000,name)
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_errpush(info,name)
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()

@ -73,16 +73,15 @@ subroutine mld_zmlprec_bld(a,desc_a,p,info)
integer :: err_act
character(len=20) :: name, ch_err
logical, parameter :: debug=.false.
type(psb_zspmat_type) :: ac
integer :: ictxt, np, me
name='psb_zmlprec_bld'
if (psb_get_errstatus().ne.0) return
call psb_erractionsave(err_act)
info = 0
ictxt = psb_cd_get_context(desc_a)
call psb_info(ictxt,me,np)
call psb_erractionsave(err_act)
if (.not.allocated(p%iprcparm)) then
info = 2222
@ -120,14 +119,10 @@ subroutine mld_zmlprec_bld(a,desc_a,p,info)
!
call mld_aggrmap_bld(p%iprcparm(mld_aggr_alg_),a,desc_a,p%nlaggr,p%mlia,info)
if(info /= 0) then
info=4010
ch_err='mld_aggrmap_bld'
call psb_errpush(info,name,a_err=ch_err)
call psb_errpush(4010,name,a_err='mld_aggrmap_bld')
goto 9999
end if
if (debug) write(0,*) 'Out from genaggrmap',p%nlaggr
!
! Build the coarse-level matrix from the fine level one, starting from
! the mapping defined by mld_aggrmap_bld and applying the aggregation
@ -137,22 +132,16 @@ subroutine mld_zmlprec_bld(a,desc_a,p,info)
call psb_nullify_desc(desc_ac)
call mld_aggrmat_asb(a,desc_a,ac,desc_ac,p,info)
if(info /= 0) then
info=4010
ch_err='mld_aggrmat_asb'
call psb_errpush(info,name,a_err=ch_err)
call psb_errpush(4010,name,a_err='mld_aggrmat_asb')
goto 9999
end if
if (debug) write(0,*) 'Out from bldaggrmat',desc_ac%matrix_data(:)
!
! Build the 'base preconditioner' corresponding to the coarse level
!
call mld_baseprc_bld(ac,desc_ac,p,info)
if (debug) write(0,*) 'Out from baseprcbld',info
if(info /= 0) then
info=4010
ch_err='mld_baseprc_bld'
call psb_errpush(info,name,a_err=ch_err)
if (info /= 0) then
call psb_errpush(4010,name,a_err='mld_baseprc_bld')
goto 9999
end if
@ -165,12 +154,10 @@ subroutine mld_zmlprec_bld(a,desc_a,p,info)
!
call psb_sp_transfer(ac,p%av(mld_ac_),info)
p%base_a => p%av(mld_ac_)
call psb_cdtransfer(desc_ac,p%desc_ac,info)
if (info==0) call psb_cdtransfer(desc_ac,p%desc_ac,info)
if (info /= 0) then
info=4010
ch_err='psb_cdtransfer'
call psb_errpush(info,name,a_err=ch_err)
call psb_errpush(4010,name,a_err='psb_cdtransfer')
goto 9999
end if
p%base_desc => p%desc_ac

@ -66,7 +66,8 @@
! If trans='N','n' then op(M^(-1)) = M^(-1);
! if trans='T','t' then op(M^(-1)) = M^(-T) (transpose of M^(-1)).
! work - complex(kind(0.d0)), dimension (:), optional, target.
! Workspace. Its size must be at least 4*psb_cd_get_local_cols(desc_data).
! Workspace. Its size must be at
! least 4*psb_cd_get_local_cols(desc_data).
!
subroutine mld_zprec_aply(prec,x,y,desc_data,info,trans,work)
@ -88,7 +89,6 @@ subroutine mld_zprec_aply(prec,x,y,desc_data,info,trans,work)
character :: trans_
complex(kind(1.d0)), pointer :: work_(:)
integer :: ictxt,np,me,err_act,iwsz
logical,parameter :: debug=.false., debugprt=.false.
character(len=20) :: name
name='mld_zprec_aply'
@ -110,8 +110,7 @@ subroutine mld_zprec_aply(prec,x,y,desc_data,info,trans,work)
iwsz = max(1,4*psb_cd_get_local_cols(desc_data))
allocate(work_(iwsz),stat=info)
if (info /= 0) then
info=4025
call psb_errpush(info,name,i_err=(/iwsz,0,0,0,0/),&
call psb_errpush(4025,name,i_err=(/iwsz,0,0,0,0/),&
& a_err='complex(kind(1.d0))')
goto 9999
end if
@ -119,10 +118,12 @@ subroutine mld_zprec_aply(prec,x,y,desc_data,info,trans,work)
end if
if (.not.(allocated(prec%baseprecv))) then
write(0,*) 'Inconsistent preconditioner: neither ML nor BASE?'
!! Error 1: should call mld_dprecbld
info=3112
call psb_errpush(info,name)
goto 9999
end if
if (size(prec%baseprecv) >1) then
if (debug) write(0,*) 'Into mlprec_aply',size(x),size(y)
call mld_mlprec_aply(zone,prec%baseprecv,x,zzero,y,desc_data,trans_,work_,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='mld_zmlprec_aply')
@ -132,7 +133,10 @@ subroutine mld_zprec_aply(prec,x,y,desc_data,info,trans,work)
else if (size(prec%baseprecv) == 1) then
call mld_baseprec_aply(zone,prec%baseprecv(1),x,zzero,y,desc_data,trans_, work_,info)
else
write(0,*) 'Inconsistent preconditioner: size of baseprecv???'
info = 4013
call psb_errpush(info,name,a_err='Invalid size of baseprecv',&
& i_Err=(/size(prec%baseprecv),0,0,0,0/))
goto 9999
endif
if (present(work)) then
@ -202,11 +206,11 @@ subroutine mld_zprec_aply1(prec,x,desc_data,info,trans)
character(len=1), optional :: trans
! Local variables
logical,parameter :: debug=.false., debugprt=.false.
character :: trans_
integer :: ictxt,np,me, err_act
complex(kind(1.d0)), pointer :: WW(:), w1(:)
character(len=20) :: name
name='mld_zprec_aply1'
info = 0
call psb_erractionsave(err_act)
@ -227,17 +231,25 @@ subroutine mld_zprec_aply1(prec,x,desc_data,info,trans)
& a_err='complex(kind(1.d0))')
goto 9999
end if
if (debug) write(0,*) 'prec_aply1 Size(x) ',size(x), size(ww),size(w1)
call mld_zprec_aply(prec,x,ww,desc_data,info,trans_,work=w1)
if(info /=0) goto 9999
call mld_precaply(prec,x,ww,desc_data,info,trans_,work=w1)
if (info /= 0) then
call psb_errpush(4010,name,a_err='mld_precaply')
goto 9999
end if
x(:) = ww(:)
deallocate(ww,W1)
deallocate(ww,W1,stat=info)
if (info /= 0) then
info = 4000
call psb_errpush(info,name)
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_errpush(info,name)
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()

@ -41,7 +41,7 @@
! Contains: subroutine init_baseprc_av
!
! This routine builds the preconditioner according to the requirements made by
! the user trough the subroutines mld_zprecinit and mld_zprecset.
! the user trough the subroutines mld_precinit and mld_precset.
!
! A multilevel preconditioner is regarded as an array of 'base preconditioners',
! each representing the part of the preconditioner associated to a certain level.
@ -76,32 +76,36 @@ subroutine mld_zprecbld(a,desc_a,p,info,upd)
! Local Variables
Integer :: err,i,k,ictxt, me,np, err_act
Integer :: err,i,k,ictxt, me,np, err_act, iszv
integer :: int_err(5)
character :: iupd
logical, parameter :: debug=.false.
integer,parameter :: iroot=0,iout=60,ilout=40
character(len=20) :: name, ch_err
integer :: debug_level, debug_unit
character(len=20) :: name, ch_err
if(psb_get_errstatus().ne.0) return
info=0
err=0
call psb_erractionsave(err_act)
name = 'mld_zprecbld'
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
if (debug) write(0,*) 'Entering precbld',desc_a%matrix_data(:)
name = 'mld_zprecbld'
info = 0
int_err(1) = 0
ictxt = psb_cd_get_context(desc_a)
if (debug) write(0,*) 'Preconditioner psb_info'
call psb_info(ictxt, me, np)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Entering ',desc_a%matrix_data(:)
if (present(upd)) then
if (debug) write(0,*) 'UPD ', upd
if ((upd.eq.'F').or.(upd.eq.'T')) then
iupd=upd
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),'UPD ', upd
if ((toupper(upd).eq.'F').or.(toupper(upd).eq.'T')) then
iupd=toupper(upd)
else
iupd='F'
endif
@ -110,86 +114,81 @@ subroutine mld_zprecbld(a,desc_a,p,info,upd)
endif
if (.not.allocated(p%baseprecv)) then
!! Error 1: should call mld_dprecset
info=4010
ch_err='unallocated bpv'
call psb_errpush(info,name,a_err=ch_err)
!! Error: should have called mld_dprecinit
info=3111
call psb_errpush(info,name)
goto 9999
end if
!
! Should add check to ensure all procs have the same ...
! Check to ensure all procs have the same
!
iszv = size(p%baseprecv)
call psb_bcast(ictxt,iszv)
if (iszv /= size(p%baseprecv)) then
info=4001
call psb_errpush(info,name,a_err='Inconsistent size of baseprecv')
goto 9999
end if
if (size(p%baseprecv) >= 1) then
if (iszv >= 1) then
!
! Allocate the av component of the preconditioner data type
! at the finest level
! Allocate and build the fine level preconditioner
!
call init_baseprc_av(p%baseprecv(1),info)
if (info == 0) call mld_baseprc_bld(a,desc_a,p%baseprecv(1),info,iupd)
if (info /= 0) then
info=4010
ch_err='allocate'
call psb_errpush(info,name,a_err=ch_err)
call psb_errpush(4001,name,a_err='Base level precbuild.')
goto 9999
endif
!
! Build the base preconditioner corresponding to the finest
! level
!
call mld_baseprc_bld(a,desc_a,p%baseprecv(1),info,iupd)
end if
else
info=4010
ch_err='size bpv'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
if (size(p%baseprecv) > 1) then
if (iszv > 1) then
!
! Build the base preconditioners corresponding to the remaining
! levels
!
do i=2, size(p%baseprecv)
do i=2, iszv
!
! Allocate the av component of the preconditioner data type
! at level i
!
call init_baseprc_av(p%baseprecv(i),info)
if (info /= 0) then
info=4010
ch_err='allocate'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
if (i<size(p%baseprecv)) then
if (i<iszv) then
!
! A replicated matrix only makes sense at the coarsest level
!
call mld_check_def(p%baseprecv(i)%iprcparm(mld_coarse_mat_),'Coarse matrix',&
& mld_distr_mat_,is_distr_ml_coarse_mat)
end if
call init_baseprc_av(p%baseprecv(i),info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Calling mlprcbld at level ',i
!
! Build the base preconditioner corresponding to level i
!
call mld_mlprec_bld(p%baseprecv(i-1)%base_a,p%baseprecv(i-1)%base_desc,&
& p%baseprecv(i),info)
if (info == 0) call mld_mlprec_bld(p%baseprecv(i-1)%base_a,&
& p%baseprecv(i-1)%base_desc, p%baseprecv(i),info)
if (info /= 0) then
info=4010
call psb_errpush(info,name)
call psb_errpush(4001,name,a_err='Init & build upper level preconditioner')
goto 9999
endif
if (debug) then
write(0,*) 'Return from ',i-1,' call to mlprcbld ',info
endif
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Return from ',i,' call to mlprcbld ',info
end do
endif
@ -211,12 +210,15 @@ contains
type(mld_zbaseprc_type), intent(inout) :: p
integer :: info
if (allocated(p%av)) then
!
! We have not yet decided what to do
!
if (size(p%av) /= mld_max_avsz_) then
deallocate(p%av,stat=info)
if (info /= 0) return
endif
end if
if (.not.(allocated(p%av))) then
allocate(p%av(mld_max_avsz_),stat=info)
if (info /= 0) return
end if
allocate(p%av(mld_max_avsz_),stat=info)
!!$ if (info /= 0) return
do k=1,size(p%av)
call psb_nullify_sp(p%av(k))
end do

@ -182,9 +182,6 @@ subroutine mld_zprecinit(p,ptype,info,nlev)
else
nlev_ = 2
end if
if (nlev_ == 1) then
write(0,*) 'Warning: requested ML preconditioner with NLEV=1'
endif
ilev_ = 1
allocate(p%baseprecv(nlev_),stat=info)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)

@ -87,8 +87,7 @@ subroutine mld_zprecseti(p,what,val,info,ilev)
info = 0
if (.not.allocated(p%baseprecv)) then
write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
info = -1
info = 3111
return
endif
nlev_ = size(p%baseprecv)
@ -100,13 +99,11 @@ subroutine mld_zprecseti(p,what,val,info,ilev)
end if
if ((ilev_<1).or.(ilev_ > nlev_)) then
write(0,*) 'PRECSET ERRROR: ilev out of bounds'
info = -1
return
endif
if (.not.allocated(p%baseprecv(ilev_)%iprcparm)) then
write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
info = -1
info = 3111
return
endif
@ -251,8 +248,7 @@ subroutine mld_zprecsetc(p,what,string,info,ilev)
info = 0
if (.not.allocated(p%baseprecv)) then
write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
info = -1
info = 3111
return
endif
nlev_ = size(p%baseprecv)
@ -269,8 +265,7 @@ subroutine mld_zprecsetc(p,what,string,info,ilev)
return
endif
if (.not.allocated(p%baseprecv(ilev_)%iprcparm)) then
write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
info = -1
info = 3111
return
endif
@ -445,8 +440,7 @@ subroutine mld_zprecsetd(p,what,val,info,ilev)
end if
if (.not.allocated(p%baseprecv)) then
write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
info = -1
info = 3111
return
endif
nlev_ = size(p%baseprecv)
@ -457,8 +451,7 @@ subroutine mld_zprecsetd(p,what,val,info,ilev)
return
endif
if (.not.allocated(p%baseprecv(ilev_)%dprcparm)) then
write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
info = -1
info = 3111
return
endif

@ -82,7 +82,6 @@ subroutine mld_zslu_bld(a,desc_a,p,info)
! Local variables
integer :: nzt,ictxt,me,np,err_act
logical, parameter :: debug=.false.
character(len=20) :: name, ch_err
if(psb_get_errstatus().ne.0) return
@ -95,19 +94,12 @@ subroutine mld_zslu_bld(a,desc_a,p,info)
call psb_info(ictxt, me, np)
if (toupper(a%fida) /= 'CSR') then
write(0,*) 'Unimplemented input to mld_slu_BLD'
info=135
call psb_errpush(info,name,a_err=a%fida)
goto 9999
endif
nzt = psb_sp_get_nnzeros(a)
if (Debug) then
write(0,*) me,'Calling mld_slu_factor ',nzt,a%m,&
& a%k,p%desc_data%matrix_data(psb_n_row_)
call psb_barrier(ictxt)
endif
!
! Compute the LU factorization
!
@ -120,11 +112,6 @@ subroutine mld_zslu_bld(a,desc_a,p,info)
goto 9999
end if
if (Debug) then
write(0,*) me, 'SPLUBLD: Done mld_slu_Factor',info,p%iprcparm(mld_slu_ptr_)
call psb_barrier(ictxt)
endif
call psb_erractionrestore(err_act)
return

@ -80,8 +80,7 @@ subroutine mld_zsludist_bld(a,desc_a,p,info)
! Local variables
integer :: nzt,ictxt,me,np,err_act,&
& mglob,ifrst,ibcheck,nrow,ncol,npr,npc
logical, parameter :: debug=.false.
character(len=20) :: name, ch_err
character(len=20) :: name, ch_err
if(psb_get_errstatus().ne.0) return
info=0
@ -93,7 +92,8 @@ subroutine mld_zsludist_bld(a,desc_a,p,info)
call psb_info(ictxt, me, np)
if (toupper(a%fida) /= 'CSR') then
write(0,*) 'Unimplemented input to mld_slu_BLD'
info=135
call psb_errpush(info,name,a_err=a%fida)
goto 9999
endif

@ -87,10 +87,9 @@ subroutine mld_zumf_bld(a,desc_a,p,info)
integer, intent(out) :: info
! Local variables
integer :: nzt,ictxt,me,np,err_act
integer :: i_err(5)
logical, parameter :: debug=.false.
character(len=20) :: name
integer :: nzt,ictxt,me,np,err_act
integer :: i_err(5)
character(len=20) :: name
info=0
name='mld_zumf_bld'
@ -104,18 +103,8 @@ subroutine mld_zumf_bld(a,desc_a,p,info)
goto 9999
endif
nzt = psb_sp_get_nnzeros(a)
if (Debug) then
write(0,*) me,'Calling mld_umf_factor ',nzt,a%m,&
& a%k,p%desc_data%matrix_data(psb_n_row_)
open(80+me)
call psb_csprt(80+me,a)
close(80+me)
call psb_barrier(ictxt)
endif
!
! Compute the LU factorization
!
@ -130,11 +119,6 @@ subroutine mld_zumf_bld(a,desc_a,p,info)
goto 9999
end if
if (Debug) then
write(0,*) me, 'UMFBLD: Done mld_umf_Factor',info,p%iprcparm(mld_umf_numptr_)
call psb_barrier(ictxt)
endif
call psb_erractionrestore(err_act)
return

Loading…
Cancel
Save