fdurastante/amg4psblas
1
0
Fork 0
mirror of https://github.com/sfilippone/amg4psblas.git synced 2026-10-06 14:44:55 +00:00
Code Issues Packages Projects Releases Wiki Activity

mld2p4-extaggr:

Browse Source 
Make.inc.in
 Makefile
 mlprec/impl/Makefile
 mlprec/impl/mld_d_bld_mlhier_array.f90
 mlprec/impl/mld_d_dec_map_bld.f90
 mlprec/impl/mld_d_hierarchy_bld.f90
 mlprec/impl/mld_d_lev_aggrmap_bld.f90
 mlprec/impl/mld_d_lev_aggrmat_asb.f90
 mlprec/impl/mld_daggrmap_bld.f90
 mlprec/impl/mld_daggrmat_asb.f90
 mlprec/impl/mld_daggrmat_biz_asb.f90
 mlprec/impl/mld_daggrmat_minnrg_asb.f90
 mlprec/impl/mld_daggrmat_nosmth_asb.f90
 mlprec/impl/mld_daggrmat_smth_asb.f90
 mlprec/impl/mld_dcoarse_bld.f90
 mlprec/impl/mld_dprecinit.F90
 mlprec/mld_base_prec_type.F90
 mlprec/mld_d_inner_mod.f90
 tests/pdegen/mld_d_pde3d.f90
 tests/pdegen/runs/ppde.inp

Refactored map bld/mat asb routines. Now there is no need to call
coarse_bld twice. To be cleaned up yet.
This commit is contained in:
Salvatore Filippone
2016-09-30 13:59:41 +00:00
parent ac521deb6b
commit e1c05d0e22
20 changed files with 1103 additions and 434 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Make.inc.in
+1
View File
@@ -15,6 +15,7 @@ INSTALL_DIR=@INSTALL_DIR@
INSTALL_LIBDIR=@INSTALL_LIBDIR@
INSTALL_INCLUDEDIR=@INSTALL_INCLUDEDIR@
INSTALL_DOCSDIR=@INSTALL_DOCSDIR@
INSTALL_SAMPLESDIR=@INSTALL_SAMPLESDIR@
##########################################################
Makefile
+4
View File
@@ -29,6 +29,10 @@ install: all
/bin/cp -fr docs/*pdf docs/html $(INSTALL_DOCSDIR))
(./mkdir.sh $(INSTALL_DOCSDIR) && \
$(INSTALL_DATA) README LICENSE $(INSTALL_DOCSDIR))
(./mkdir.sh $(INSTALL_SAMPLESDIR) && ./mkdir.sh $(INSTALL_SAMPLESDIR)/simple &&\
./mkdir.sh $(INSTALL_SAMPLESDIR)/advanced && \
(cd examples; /bin/cp -fr pdegen fileread $(INSTALL_SAMPLESDIR)/simple ) && \
(cd tests; /bin/cp -fr pdegen fileread $(INSTALL_SAMPLESDIR)/advanced ))
cleanlib:
(cd lib; /bin/rm -f *.a *$(.mod) *$(.fh))
(cd include; /bin/rm -f *.a *$(.mod) *$(.fh))
mlprec/impl/Makefile
+1 -1
View File
@@ -26,7 +26,7 @@ DINNEROBJS= mld_dcoarse_bld.o mld_dmlprec_bld.o mld_d_bld_mlhier_aggsize.o mld
mld_d_ml_prec_bld.o mld_d_hierarchy_bld.o \
mld_dilu0_fact.o mld_diluk_fact.o mld_dilut_fact.o mld_daggrmap_bld.o \
mld_d_dec_map_bld.o mld_dmlprec_aply.o mld_daggrmat_asb.o \
$(DMPFOBJS) mld_d_extprol_bld.o
$(DMPFOBJS) mld_d_extprol_bld.o mld_d_lev_aggrmap_bld.o mld_d_lev_aggrmat_asb.o
SINNEROBJS= mld_scoarse_bld.o mld_smlprec_bld.o mld_s_bld_mlhier_aggsize.o mld_s_bld_mlhier_array.o \
mld_s_ml_prec_bld.o mld_s_hierarchy_bld.o \
mlprec/impl/mld_d_bld_mlhier_array.f90
+37 -17
View File
@@ -37,23 +37,25 @@
!!$
!!$
subroutine mld_d_bld_mlhier_array(nplevs,a,desc_a,precv,info)
subroutine mld_d_bld_mlhier_array(nplevs,casize,mnaggratio,a,desc_a,precv,info)
use psb_base_mod
use mld_d_inner_mod, mld_protect_name => mld_d_bld_mlhier_array
use mld_d_prec_mod
implicit none
integer(psb_ipk_), intent(inout) :: nplevs
integer(psb_ipk_), intent(inout) :: nplevs, casize
real(psb_dpk_) :: mnaggratio
type(psb_dspmat_type),intent(in), target :: a
type(psb_desc_type), intent(inout), target :: desc_a
type(mld_d_onelev_type),intent(inout), allocatable, target :: precv(:)
integer(psb_ipk_), intent(out) :: info
! Local
integer(psb_ipk_) :: ictxt, me,np
integer(psb_ipk_) :: err,i,k, err_act, newsz, iszv
integer(psb_ipk_) :: err,i,k, err_act, newsz, iszv, iaggsize
integer(psb_ipk_) :: ipv(mld_ifpsz_), val
class(mld_d_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm
type(mld_dml_parms) :: baseparms, medparms, coarseparms
type(mld_d_onelev_type), allocatable :: tprecv(:)
real(psb_dpk_) :: sizeratio
integer(psb_ipk_) :: int_err(5)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
@@ -178,27 +180,45 @@ subroutine mld_d_bld_mlhier_array(nplevs,a,desc_a,precv,info)
& write(debug_unit,*) me,' ',trim(name),&
& 'Return from ',i,' call to mlprcbld ',info
if (i>2) then
iaggsize = sum(precv(i)%map%naggr)
if (iaggsize <= casize) then
newsz = i
end if
if (i>2) then
sizeratio = iaggsize
sizeratio = sum(precv(i-1)%map%naggr)/sizeratio
if (sizeratio < mnaggratio) then
if (sizeratio > 1) then
newsz = i
else
!
! We are not gaining
!
newsz = newsz-1
end if
end if
if (all(precv(i)%map%naggr == precv(i-1)%map%naggr)) then
newsz=i-1
if (me == 0) then
write(debug_unit,*) trim(name),&
&': Warning: aggregates from level ',&
& newsz
write(debug_unit,*) trim(name),&
&': to level ',&
& iszv,' coincide.'
write(debug_unit,*) trim(name),&
&': Number of levels actually used :',newsz
write(debug_unit,*)
end if
end if
call psb_bcast(ictxt,newsz)
if (newsz > 0) exit array_build_loop
end if
call psb_bcast(ictxt,newsz)
if (newsz > 0) exit array_build_loop
end do array_build_loop
if (newsz > 0) then
if (me == 0) then
write(debug_unit,*) trim(name),&
&': Warning: aggregates from level ',&
& newsz
write(debug_unit,*) trim(name),&
&': to level ',&
& iszv,' coincide.'
write(debug_unit,*) trim(name),&
&': Number of levels actually used :',newsz
write(debug_unit,*)
end if
allocate(tprecv(newsz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,&
mlprec/impl/mld_d_dec_map_bld.f90
+30 -23
View File
@@ -54,12 +54,12 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
! Local variables
integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),&
& ideg(:), idxs(:)
& ideg(:), idxs(:), tmpaggr(:)
real(psb_dpk_), allocatable :: val(:), diag(:)
integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii, ip
integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii, ip, nisolate,nthird
type(psb_d_csr_sparse_mat) :: acsr
real(psb_dpk_) :: cpling, tcl
logical :: recovery, candidate
logical :: recovery, disjoint
integer(psb_ipk_) :: debug_level, debug_unit,err_act
integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: nrow, ncol, n_ne
@@ -114,6 +114,7 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
!
naggr = 0
icnt = 0
nisolate = 0
step1: do ii=1, nr
i = idxs(ii)
@@ -138,20 +139,20 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
end if
end if
enddo
if (ip < 1) then
write(0,*) "Should at least contain the node itself ! "
cycle step1
end if
candidate = .true.
do k=1, ip
candidate = candidate .and. (ilaggr(icol(k)) == -(nr+1))
end do
if (candidate) then
!
! If the whole strongly coupled neighborhood of I is
! as yet unconnected, turn it into the next aggregate
!
!!$ if (ip <= 1) then
!!$ nisolate = nisolate + 1
!!$ cycle step1
!!$ end if
!
! If the whole strongly coupled neighborhood of I is
! as yet unconnected, turn it into the next aggregate.
! Same if ip==0 (in which case, neighborhood only
! contains I even if it does not look from matrix)
!
disjoint = all(ilaggr(icol(1:ip)) == -(nr+1)).or.(ip==0)
if (disjoint) then
icnt = icnt + 1
naggr = naggr + 1
do k=1, ip
@@ -161,6 +162,7 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
end if
endif
enddo step1
if (debug_level >= psb_debug_outer_) then
write(debug_unit,*) me,' ',trim(name),&
& ' Check 1:',count(ilaggr == -(nr+1))
@@ -168,7 +170,8 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
!
! Phase two: join the neighbours
!
!
tmpaggr = ilaggr
step2: do ii=1,nr
i = idxs(ii)
@@ -188,8 +191,9 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
do k=1, nz
j = icol(k)
if ((1<=j).and.(j<=nr)) then
if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j))))&
& .and. (ilaggr(j) > 0).and. (abs(val(k)) > cpling)) then
!!$ if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j))))&
!!$ & .and. (tmpaggr(j) > 0).and. (abs(val(k)) > cpling)) then
if ((tmpaggr(j) > 0).and. (abs(val(k)) > cpling)) then
ip = k
cpling = abs(val(k))
end if
@@ -205,10 +209,12 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
!
! Phase three: sweep over leftovers, if any
!
nthird = 0
step3: do ii=1,nr
i = idxs(ii)
if (ilaggr(i) == -(nr+1)) then
if (ilaggr(i) < 0) then
nthird = nthird + 1
call a%csget(i,i,nz,irow,icol,val,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
@@ -216,8 +222,8 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
goto 9999
end if
!
! Find the most strongly connected neighbour that is
! already aggregated, if any, and join its aggregate
! Find its strongly connected neighbourhood not
! already aggregated, and make it into a new aggregate.
!
cpling = dzero
ip = 0
@@ -282,7 +288,8 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
nlaggr(:) = 0
nlaggr(me+1) = naggr
call psb_sum(ictxt,nlaggr(1:np))
!write(*,*) me,'Info from dec_map_bld: ',naggr,nisolate,nthird
call psb_erractionrestore(err_act)
return
mlprec/impl/mld_d_hierarchy_bld.f90
+338 -15
View File
@@ -94,9 +94,13 @@ subroutine mld_d_hierarchy_bld(a,desc_a,p,info,amold,vmold,imold)
! Local Variables
integer(psb_ipk_) :: ictxt, me,np
integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize, nplevs, mxplevs
real(psb_dpk_) :: mnaggratio
integer(psb_ipk_) :: ipv(mld_ifpsz_), val
integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize, nplevs, mxplevs, iaggsize
real(psb_dpk_) :: mnaggratio, sizeratio
class(mld_d_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm
type(mld_dml_parms) :: baseparms, medparms, coarseparms
integer(psb_ipk_), allocatable :: ilaggr(:), nlaggr(:)
type(psb_dspmat_type) :: op_prol
type(mld_d_onelev_type), allocatable :: tprecv(:)
integer(psb_ipk_) :: int_err(5)
character :: upd_
integer(psb_ipk_) :: debug_level, debug_unit
@@ -191,10 +195,23 @@ subroutine mld_d_hierarchy_bld(a,desc_a,p,info,amold,vmold,imold)
goto 9999
endif
!
! The strategy:
! 1. The maximum number of levels should be already encoded in the
! size of the array;
! 2. If the user did not specify anything, then a default coarse size
! is generated, and the number of levels is set to the maximum;
! 3. If the number of levels has been specified, make sure it's capped
! at the maximum;
! 4. If the size of the array is different from target number of levels,
! reallocate;
! 5. Build the matrix hierarchy, stopping early if either the target
! coarse size is hit, or the gain fall below the mmin_aggr_ratio
! threshold.
!
if (nplevs <= 0) then
!
! This should become the default strategy, we specify a target aggregation size.
!
if (casize <=0) then
!
! Default to the cubic root of the size at base level.
@@ -204,15 +221,321 @@ subroutine mld_d_hierarchy_bld(a,desc_a,p,info,amold,vmold,imold)
casize = max(casize,ione)
casize = casize*40_psb_ipk_
end if
call mld_bld_mlhier_aggsize(casize,mxplevs,mnaggratio,a,desc_a,p%precv,info)
else
!
! Oldstyle with fixed number of levels.
!
nplevs = max(itwo,min(nplevs,mxplevs))
call mld_bld_mlhier_array(nplevs,a,desc_a,p%precv,info)
nplevs = mxplevs
end if
nplevs = max(itwo,min(nplevs,mxplevs))
coarseparms = p%precv(iszv)%parms
baseparms = p%precv(1)%parms
medparms = p%precv(2)%parms
allocate(coarse_sm, source=p%precv(iszv)%sm,stat=info)
if (info == psb_success_) &
& allocate(med_sm, source=p%precv(2)%sm,stat=info)
if (info == psb_success_) &
& allocate(base_sm, source=p%precv(1)%sm,stat=info)
if (info /= psb_success_) then
write(0,*) 'Error in saving smoothers',info
call psb_errpush(psb_err_internal_error_,name,a_err='Base level precbuild.')
goto 9999
end if
!
! First set desired number of levels
!
if (iszv /= nplevs) then
allocate(tprecv(nplevs),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='prec reallocation')
goto 9999
endif
tprecv(1)%parms = baseparms
allocate(tprecv(1)%sm,source=base_sm,stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='prec reallocation')
goto 9999
endif
do i=2,nplevs-1
tprecv(i)%parms = medparms
allocate(tprecv(i)%sm,source=med_sm,stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='prec reallocation')
goto 9999
endif
end do
tprecv(nplevs)%parms = coarseparms
allocate(tprecv(nplevs)%sm,source=coarse_sm,stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='prec reallocation')
goto 9999
endif
do i=1,iszv
call p%precv(i)%free(info)
end do
call move_alloc(tprecv,p%precv)
iszv = size(p%precv)
end if
!
! Finest level first; remember to fix base_a and base_desc
!
p%precv(1)%base_a => a
p%precv(1)%base_desc => desc_a
newsz = 0
if (.false.) then
old_array_build_loop: do i=2, iszv
!
! Check on the iprcparm contents: they should be the same
! on all processes.
!
call psb_bcast(ictxt,p%precv(i)%parms)
!
! Sanity checks on the parameters
!
if (i<iszv) then
!
! A replicated matrix only makes sense at the coarsest level
!
call mld_check_def(p%precv(i)%parms%coarse_mat,'Coarse matrix',&
& mld_distr_mat_,is_distr_ml_coarse_mat)
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Calling mlprcbld at level ',i
!
! Build the mapping between levels i-1 and i and the matrix
! at level i
!
if (info == psb_success_) call mld_coarse_bld(p%precv(i-1)%base_a,&
& p%precv(i-1)%base_desc, p%precv(i),info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Init upper level preconditioner')
goto 9999
endif
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Return from ',i,' call to mlprcbld ',info
iaggsize = sum(p%precv(i)%map%naggr)
if (iaggsize <= casize) then
newsz = i
end if
if (i>2) then
sizeratio = iaggsize
sizeratio = sum(p%precv(i-1)%map%naggr)/sizeratio
if (sizeratio < mnaggratio) then
if (sizeratio > 1) then
newsz = i
else
!
! We are not gaining
!
newsz = newsz-1
end if
end if
if (all(p%precv(i)%map%naggr == p%precv(i-1)%map%naggr)) then
newsz=i-1
if (me == 0) then
write(debug_unit,*) trim(name),&
&': Warning: aggregates from level ',&
& newsz
write(debug_unit,*) trim(name),&
&': to level ',&
& iszv,' coincide.'
write(debug_unit,*) trim(name),&
&': Number of levels actually used :',newsz
write(debug_unit,*)
end if
end if
end if
call psb_bcast(ictxt,newsz)
if (newsz > 0) exit old_array_build_loop
end do old_array_build_loop
if (newsz > 0) then
!
! We exited early from the build loop, need to fix
! the size.
!
allocate(tprecv(newsz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='prec reallocation')
goto 9999
endif
do i=1,newsz-1
call p%precv(i)%move_alloc(tprecv(i),info)
end do
call p%precv(iszv)%move_alloc(tprecv(newsz),info)
do i=newsz+1, iszv
call p%precv(i)%free(info)
end do
call move_alloc(tprecv,p%precv)
! Ignore errors from transfer
info = psb_success_
!
! Restart
iszv = newsz
! Fix the pointers, but the level 1 should
! be already OK
do i=2, iszv - 1
p%precv(i)%base_a => p%precv(i)%ac
p%precv(i)%base_desc => p%precv(i)%desc_ac
p%precv(i)%map%p_desc_X => p%precv(i-1)%base_desc
p%precv(i)%map%p_desc_Y => p%precv(i)%base_desc
end do
i = iszv
if (info == psb_success_) call mld_coarse_bld(p%precv(i-1)%base_a,&
& p%precv(i-1)%base_desc, p%precv(i),info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='coarse rebuild')
goto 9999
endif
end if
else
array_build_loop: do i=2, iszv
!
! Check on the iprcparm contents: they should be the same
! on all processes.
!
call psb_bcast(ictxt,p%precv(i)%parms)
!
! Sanity checks on the parameters
!
if (i<iszv) then
!
! A replicated matrix only makes sense at the coarsest level
!
call mld_check_def(p%precv(i)%parms%coarse_mat,'Coarse matrix',&
& mld_distr_mat_,is_distr_ml_coarse_mat)
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Calling mlprcbld at level ',i
!
! Build the mapping between levels i-1 and i and the matrix
! at level i
!
if (info == psb_success_) call mld_aggrmap_bld(p%precv(i),&
& p%precv(i-1)%base_a,p%precv(i-1)%base_desc,&
& ilaggr,nlaggr,op_prol,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Map build')
goto 9999
endif
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Return from ',i,' call to mlprcbld ',info
!
! Check for early termination of aggregation loop.
!
iaggsize = sum(nlaggr)
if (iaggsize <= casize) then
newsz = i
end if
if (i>2) then
sizeratio = iaggsize
sizeratio = sum(p%precv(i-1)%map%naggr)/sizeratio
if (sizeratio < mnaggratio) then
if (sizeratio > 1) then
newsz = i
else
!
! We are not gaining
!
newsz = newsz-1
end if
end if
if (all(nlaggr == p%precv(i-1)%map%naggr)) then
newsz=i-1
if (me == 0) then
write(debug_unit,*) trim(name),&
&': Warning: aggregates from level ',&
& newsz
write(debug_unit,*) trim(name),&
&': to level ',&
& iszv,' coincide.'
write(debug_unit,*) trim(name),&
&': Number of levels actually used :',newsz
write(debug_unit,*)
end if
end if
end if
call psb_bcast(ictxt,newsz)
if (newsz > 0) &
& call p%precv(i)%parms%get_coarse(p%precv(iszv)%parms)
if (info == psb_success_) call mld_aggrmat_asb(p%precv(i),&
& p%precv(i-1)%base_a,p%precv(i-1)%base_desc,&
& ilaggr,nlaggr,op_prol,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Map build')
goto 9999
endif
if (newsz > 0) exit array_build_loop
end do array_build_loop
if (newsz > 0) then
!
! We exited early from the build loop, need to fix
! the size.
!
allocate(tprecv(newsz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='prec reallocation')
goto 9999
endif
do i=1,newsz
call p%precv(i)%move_alloc(tprecv(i),info)
end do
do i=newsz+1, iszv
call p%precv(i)%free(info)
end do
call move_alloc(tprecv,p%precv)
! Ignore errors from transfer
info = psb_success_
!
! Restart
iszv = newsz
! Fix the pointers, but the level 1 should
! be already OK
do i=2, iszv
p%precv(i)%base_a => p%precv(i)%ac
p%precv(i)%base_desc => p%precv(i)%desc_ac
p%precv(i)%map%p_desc_X => p%precv(i-1)%base_desc
p%precv(i)%map%p_desc_Y => p%precv(i)%base_desc
end do
end if
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Internal hierarchy build' )
@@ -220,7 +543,7 @@ subroutine mld_d_hierarchy_bld(a,desc_a,p,info,amold,vmold,imold)
endif
iszv = size(p%precv)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Exiting with',iszv,' levels'
mlprec/impl/mld_d_lev_aggrmap_bld.f90
+148
View File
@@ -0,0 +1,148 @@
!!$
!!$
!!$ MLD2P4 version 2.0
!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package
!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3)
!!$
!!$ (C) Copyright 2008, 2010, 2012, 2015
!!$
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$ Pasqua D'Ambra ICAR-CNR, Naples
!!$ Daniela di Serafino Second University of Naples
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the MLD2P4 group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
! File: mld_dcoarse_bld.f90
!
! Subroutine: mld_dcoarse_bld
! Version: real
!
! This routine builds the matrix associated to the current level of the
! multilevel preconditioner from the matrix associated to the previous level,
! by using a smoothed aggregation technique (therefore, it also builds the
! prolongation and restriction operators mapping the current level to the
! previous one and vice versa). Then the routine builds the base preconditioner
! at the current level.
! The current level is regarded as the coarse one, while the previous as
! the fine one. This is in agreement with the fact that the routine is called,
! by mld_mlprec_bld, only on levels >=2.
!
!
! Arguments:
! a - type(psb_dspmat_type).
! The sparse matrix structure containing the local part of the
! fine-level matrix.
! desc_a - type(psb_desc_type), input.
! The communication descriptor of a.
! p - type(mld_d_onelev_type), input/output.
! The 'one-level' data structure containing the local part
! of the base preconditioner to be built as well as
! information concerning the prolongator and its transpose.
! info - integer, output.
! Error code.
!
subroutine mld_d_lev_aggrmap_bld(p,a,desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod
use mld_d_inner_mod, mld_protect_name => mld_d_lev_aggrmap_bld
implicit none
! Arguments
type(mld_d_onelev_type), intent(inout), target :: p
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
type(psb_dspmat_type), intent(out) :: op_prol
integer(psb_ipk_), intent(out) :: info
! Local variables
character(len=20) :: name
integer(psb_mpik_) :: ictxt, np, me
integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: nzl, ntaggr
integer(psb_ipk_) :: debug_level, debug_unit
name='mld_d_lev_aggrmap_bld'
if (psb_get_errstatus().ne.0) return
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
info = psb_success_
ictxt = desc_a%get_context()
call psb_info(ictxt,me,np)
call mld_check_def(p%parms%ml_type,'Multilevel type',&
& mld_mult_ml_,is_legal_ml_type)
call mld_check_def(p%parms%aggr_alg,'Aggregation',&
& mld_dec_aggr_,is_legal_ml_aggr_alg)
call mld_check_def(p%parms%aggr_ord,'Ordering',&
& mld_aggr_ord_nat_,is_legal_ml_aggr_ord)
call mld_check_def(p%parms%aggr_thresh,'Aggr_Thresh',dzero,is_legal_d_aggr_thrs)
select case(p%parms%aggr_alg)
case (mld_dec_aggr_, mld_sym_dec_aggr_)
!
! Build a mapping between the row indices of the fine-level matrix
! and the row indices of the coarse-level matrix, according to a decoupled
! aggregation algorithm. This also defines a tentative prolongator from
! the coarse to the fine level.
!
call mld_aggrmap_bld(p%parms%aggr_alg,p%parms%aggr_ord,p%parms%aggr_thresh,&
& a,desc_a,ilaggr,nlaggr,op_prol,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmap_bld')
goto 9999
end if
case (mld_bcmatch_aggr_)
write(0,*) 'Matching is not implemented yet '
info = -1111
call psb_errpush(psb_err_input_value_invalid_i_,name,&
& i_err=(/ione,p%parms%aggr_alg,izero,izero,izero/))
goto 9999
case default
info = -1
call psb_errpush(psb_err_input_value_invalid_i_,name,&
& i_err=(/ione,p%parms%aggr_alg,izero,izero,izero/))
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_d_lev_aggrmap_bld
mlprec/impl/mld_d_lev_aggrmat_asb.f90
+252
View File
@@ -0,0 +1,252 @@
!!$
!!$
!!$ MLD2P4 version 2.0
!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package
!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3)
!!$
!!$ (C) Copyright 2008, 2010, 2012, 2015
!!$
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$ Pasqua D'Ambra ICAR-CNR, Naples
!!$ Daniela di Serafino Second University of Naples
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the MLD2P4 group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
! File: mld_dcoarse_bld.f90
!
! Subroutine: mld_dcoarse_bld
! Version: real
!
! This routine builds the matrix associated to the current level of the
! multilevel preconditioner from the matrix associated to the previous level,
! by using a smoothed aggregation technique (therefore, it also builds the
! prolongation and restriction operators mapping the current level to the
! previous one and vice versa). Then the routine builds the base preconditioner
! at the current level.
! The current level is regarded as the coarse one, while the previous as
! the fine one. This is in agreement with the fact that the routine is called,
! by mld_mlprec_bld, only on levels >=2.
!
!
! Arguments:
! a - type(psb_dspmat_type).
! The sparse matrix structure containing the local part of the
! fine-level matrix.
! desc_a - type(psb_desc_type), input.
! The communication descriptor of a.
! p - type(mld_d_onelev_type), input/output.
! The 'one-level' data structure containing the local part
! of the base preconditioner to be built as well as
! information concerning the prolongator and its transpose.
! info - integer, output.
! Error code.
!
subroutine mld_d_lev_aggrmat_asb(p,a,desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod
use mld_d_inner_mod, mld_protect_name => mld_d_lev_aggrmat_asb
implicit none
! Arguments
type(mld_d_onelev_type), intent(inout), target :: p
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:),nlaggr(:)
type(psb_dspmat_type), intent(inout) :: op_prol
integer(psb_ipk_), intent(out) :: info
! Local variables
character(len=20) :: name
integer(psb_mpik_) :: ictxt, np, me
integer(psb_ipk_) :: err_act
type(psb_dspmat_type) :: ac, op_restr
type(psb_d_coo_sparse_mat) :: acoo, bcoo
type(psb_d_csr_sparse_mat) :: acsr1
integer(psb_ipk_) :: nzl, ntaggr
integer(psb_ipk_) :: debug_level, debug_unit
name='mld_dcoarse_bld'
if (psb_get_errstatus().ne.0) return
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
info = psb_success_
ictxt = desc_a%get_context()
call psb_info(ictxt,me,np)
call mld_check_def(p%parms%aggr_kind,'Smoother',&
& mld_smooth_prol_,is_legal_ml_aggr_kind)
call mld_check_def(p%parms%coarse_mat,'Coarse matrix',&
& mld_distr_mat_,is_legal_ml_coarse_mat)
call mld_check_def(p%parms%aggr_filter,'Use filtered matrix',&
& mld_no_filter_mat_,is_legal_aggr_filter)
call mld_check_def(p%parms%smoother_pos,'smooth_pos',&
& mld_pre_smooth_,is_legal_ml_smooth_pos)
call mld_check_def(p%parms%aggr_omega_alg,'Omega Alg.',&
& mld_eig_est_,is_legal_ml_aggr_omega_alg)
call mld_check_def(p%parms%aggr_eig,'Eigenvalue estimate',&
& mld_max_norm_,is_legal_ml_aggr_eig)
call mld_check_def(p%parms%aggr_omega_val,'Omega',dzero,is_legal_d_omega)
!
! Build the coarse-level matrix from the fine-level one, starting from
! the mapping defined by mld_aggrmap_bld and applying the aggregation
! algorithm specified by p%iprcparm(mld_aggr_kind_)
!
call mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p%parms,ac,op_prol,op_restr,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_asb')
goto 9999
end if
! Common code refactored here.
ntaggr = sum(nlaggr)
select case(p%parms%coarse_mat)
case(mld_distr_mat_)
call ac%mv_to(bcoo)
if (p%parms%clean_zeros) call bcoo%clean_zeros(info)
nzl = bcoo%get_nzeros()
if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1))
if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%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.'
call p%ac%mv_from(bcoo)
call p%ac%set_nrows(p%desc_ac%get_local_rows())
call p%ac%set_ncols(p%desc_ac%get_local_cols())
call p%ac%set_asb()
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free')
goto 9999
end if
if (np>1) then
call op_prol%mv_to(acsr1)
nzl = acsr1%get_nzeros()
call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I')
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc')
goto 9999
end if
call op_prol%mv_from(acsr1)
endif
call op_prol%set_ncols(p%desc_ac%get_local_cols())
if (np>1) then
call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_)
call op_restr%mv_to(acoo)
nzl = acoo%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I')
call acoo%set_dupl(psb_dupl_add_)
if (info == psb_success_) call op_restr%mv_from(acoo)
if (info == psb_success_) call op_restr%cscnv(info,type='csr')
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Converting op_restr to local')
goto 9999
end if
end if
!
! Clip to local rows.
!
call op_restr%set_nrows(p%desc_ac%get_local_rows())
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done ac '
case(mld_repl_mat_)
!
!
call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if ((info == psb_success_).and.p%parms%clean_zeros) call ac%clean_zeros(info)
if (info == psb_success_) &
& call psb_gather(p%ac,ac,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if (info /= psb_success_) goto 9999
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
goto 9999
end if
!
! Copy the prolongation/restriction matrices into the descriptor map.
! op_restr => PR^T i.e. restriction operator
! op_prol => PR i.e. prolongation operator
!
p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999
end if
!
! Fix the base_a and base_desc pointers for handling of residuals.
! This is correct because this routine is only called at levels >=2.
!
p%base_a => p%ac
p%base_desc => p%desc_ac
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_d_lev_aggrmat_asb
mlprec/impl/mld_daggrmap_bld.f90
+28 -3
View File
@@ -79,7 +79,7 @@
! info - integer, output.
! Error code.
!
subroutine mld_daggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info)
subroutine mld_daggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod
use mld_d_inner_mod, mld_protect_name => mld_daggrmap_bld
@@ -92,14 +92,15 @@ subroutine mld_daggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info)
real(psb_dpk_), intent(in) :: theta
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
type(psb_dspmat_type), intent(out) :: op_prol
integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_), allocatable :: ils(:), neigh(:)
integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m
integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m,naggrm1, naggrp1, ntaggr
type(psb_dspmat_type) :: atmp, atrans
logical :: recovery
type(psb_d_coo_sparse_mat) :: tmpcoo
integer(psb_ipk_) :: debug_level, debug_unit,err_act
integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: nrow, ncol, n_ne
@@ -151,6 +152,30 @@ subroutine mld_daggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info)
goto 9999
end if
naggr = nlaggr(me+1)
ntaggr = sum(nlaggr)
naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1))
ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1
call psb_halo(ilaggr,desc_a,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo')
goto 9999
end if
call tmpcoo%allocate(ncol,ntaggr,ncol)
do i=1,ncol
tmpcoo%val(i) = done
tmpcoo%ia(i) = i
tmpcoo%ja(i) = ilaggr(i)
end do
call tmpcoo%set_nzeros(ncol)
call tmpcoo%set_dupl(psb_dupl_add_)
call tmpcoo%set_sorted() ! At this point this is in row-major
call op_prol%mv_from(tmpcoo)
call psb_erractionrestore(err_act)
return
mlprec/impl/mld_daggrmat_asb.f90
-110
View File
@@ -166,116 +166,6 @@ subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,inf
goto 9999
end if
!!$
!!$ ntaggr = sum(nlaggr)
!!$
!!$ select case(p%parms%coarse_mat)
!!$
!!$ case(mld_distr_mat_)
!!$
!!$ call ac%mv_to(bcoo)
!!$ if (p%parms%clean_zeros) call bcoo%clean_zeros(info)
!!$ nzl = bcoo%get_nzeros()
!!$
!!$ if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1))
!!$ if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
!!$ if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
!!$ if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
!!$ if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
!!$ if (info /= psb_success_) then
!!$ call psb_errpush(psb_err_internal_error_,name,&
!!$ & a_err='Creating p%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.'
!!$ call p%ac%mv_from(bcoo)
!!$
!!$ call p%ac%set_nrows(p%desc_ac%get_local_rows())
!!$ call p%ac%set_ncols(p%desc_ac%get_local_cols())
!!$ call p%ac%set_asb()
!!$
!!$ if (info /= psb_success_) then
!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free')
!!$ goto 9999
!!$ end if
!!$
!!$ if (np>1) then
!!$ call op_prol%mv_to(acsr1)
!!$ nzl = acsr1%get_nzeros()
!!$ call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I')
!!$ if(info /= psb_success_) then
!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc')
!!$ goto 9999
!!$ end if
!!$ call op_prol%mv_from(acsr1)
!!$ endif
!!$ call op_prol%set_ncols(p%desc_ac%get_local_cols())
!!$
!!$ if (np>1) then
!!$ call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_)
!!$ call op_restr%mv_to(acoo)
!!$ nzl = acoo%get_nzeros()
!!$ if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I')
!!$ call acoo%set_dupl(psb_dupl_add_)
!!$ if (info == psb_success_) call op_restr%mv_from(acoo)
!!$ if (info == psb_success_) call op_restr%cscnv(info,type='csr')
!!$ if(info /= psb_success_) then
!!$ call psb_errpush(psb_err_internal_error_,name,&
!!$ & a_err='Converting op_restr to local')
!!$ goto 9999
!!$ end if
!!$ end if
!!$ !
!!$ ! Clip to local rows.
!!$ !
!!$ call op_restr%set_nrows(p%desc_ac%get_local_rows())
!!$
!!$ if (debug_level >= psb_debug_outer_) &
!!$ & write(debug_unit,*) me,' ',trim(name),&
!!$ & 'Done ac '
!!$
!!$ case(mld_repl_mat_)
!!$ !
!!$ !
!!$ call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
!!$ if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
!!$ if ((info == psb_success_).and.p%parms%clean_zeros) call ac%clean_zeros(info)
!!$ if (info == psb_success_) &
!!$ & call psb_gather(p%ac,ac,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
!!$
!!$ if (info /= psb_success_) goto 9999
!!$
!!$ case default
!!$ info = psb_err_internal_error_
!!$ call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
!!$ goto 9999
!!$ end select
!!$
!!$ call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
!!$ if(info /= psb_success_) then
!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
!!$ goto 9999
!!$ end if
!!$
!!$ !
!!$ ! Copy the prolongation/restriction matrices into the descriptor map.
!!$ ! op_restr => PR^T i.e. restriction operator
!!$ ! op_prol => PR i.e. prolongation operator
!!$ !
!!$
!!$ p%map = psb_linmap(psb_map_aggr_,desc_a,&
!!$ & p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
!!$ if (info == psb_success_) call op_prol%free()
!!$ if (info == psb_success_) call op_restr%free()
!!$ if(info /= psb_success_) then
!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
!!$ goto 9999
!!$ end if
call psb_erractionrestore(err_act)
return
mlprec/impl/mld_daggrmat_biz_asb.f90
+4 -21
View File
@@ -89,7 +89,8 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_dml_parms), intent(inout) :: parms
type(psb_dspmat_type), intent(out) :: ac,op_prol,op_restr
type(psb_dspmat_type), intent(inout) :: op_prol
type(psb_dspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
@@ -128,19 +129,8 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
naggr = nlaggr(me+1)
ntaggr = sum(nlaggr)
naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1))
filter_mat = (parms%aggr_filter == mld_filter_mat_)
ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1
call psb_halo(ilaggr,desc_a,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo')
goto 9999
end if
! naggr: number of local aggregates
! nrow: local rows.
!
@@ -157,17 +147,10 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
end if
! 1. Allocate Ptilde in sparse matrix form
call tmpcoo%allocate(ncol,naggr,ncol)
do i=1,nrow
tmpcoo%val(i) = done
tmpcoo%ia(i) = i
tmpcoo%ja(i) = ilaggr(i)
end do
call tmpcoo%set_nzeros(nrow)
call tmpcoo%set_dupl(psb_dupl_add_)
call op_prol%mv_to(tmpcoo)
call ptilde%mv_from_coo(tmpcoo,info)
if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (info /= psb_success_) goto 9999
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
mlprec/impl/mld_daggrmat_minnrg_asb.f90
+3 -20
View File
@@ -109,7 +109,8 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_dml_parms), intent(inout) :: parms
type(psb_dspmat_type), intent(out) :: ac,op_prol,op_restr
type(psb_dspmat_type), intent(inout) :: op_prol
type(psb_dspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
@@ -167,14 +168,6 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
filter_mat = (parms%aggr_filter == mld_filter_mat_)
ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1
call psb_halo(ilaggr,desc_a,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo')
goto 9999
end if
! naggr: number of local aggregates
! nrow: local rows.
!
@@ -209,20 +202,10 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
! 1. Allocate Ptilde in sparse matrix form
call tmpcoo%allocate(ncol,ntaggr,ncol)
do i=1,ncol
tmpcoo%val(i) = done
tmpcoo%ia(i) = i
tmpcoo%ja(i) = ilaggr(i)
end do
call tmpcoo%set_nzeros(ncol)
call tmpcoo%set_dupl(psb_dupl_add_)
call tmpcoo%set_asb()
call op_prol%mv_to(tmpcoo)
call ptilde%mv_from(tmpcoo)
call ptilde%cscnv(info,type='csr')
!!$ call local_dump(me,ptilde,'csr-ptilde','Ptilde-1')
if (info == psb_success_) call a%cscnv(am3,info,type='csr',dupl=psb_dupl_add_)
if (info == psb_success_) call a%cscnv(da,info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) then
mlprec/impl/mld_daggrmat_nosmth_asb.f90
+5 -26
View File
@@ -92,7 +92,8 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_dml_parms), intent(inout) :: parms
type(psb_dspmat_type), intent(out) :: ac,op_prol,op_restr
type(psb_dspmat_type), intent(inout) :: op_prol
type(psb_dspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
@@ -124,34 +125,12 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
ntaggr = sum(nlaggr)
naggrm1=sum(nlaggr(1:me))
do i=1, nrow
ilaggr(i) = ilaggr(i) + naggrm1
end do
call psb_halo(ilaggr,desc_a,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo')
goto 9999
end if
call acoo%allocate(ncol,ntaggr,ncol)
do i=1,nrow
acoo%val(i) = done
acoo%ia(i) = i
acoo%ja(i) = ilaggr(i)
end do
call acoo%set_dupl(psb_dupl_add_)
call acoo%set_nzeros(nrow)
call acoo%set_asb()
call acoo%fix(info)
call op_prol%mv_from(acoo)
call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info == psb_success_) call op_prol%transp(op_restr)
if (info /= psb_success_) goto 9999
call op_prol%transp(op_restr)
call a%cp_to(ac_coo)
nzt = ac_coo%get_nzeros()
mlprec/impl/mld_daggrmat_smth_asb.f90
+4 -17
View File
@@ -104,7 +104,8 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_dml_parms), intent(inout) :: parms
type(psb_dspmat_type), intent(out) :: ac,op_prol,op_restr
type(psb_dspmat_type), intent(inout) :: op_prol
type(psb_dspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
@@ -147,13 +148,6 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
naggrp1 = sum(nlaggr(1:me+1))
filter_mat = (parms%aggr_filter == mld_filter_mat_)
ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1
call psb_halo(ilaggr,desc_a,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo')
goto 9999
end if
! naggr: number of local aggregates
! nrow: local rows.
@@ -172,17 +166,10 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
end if
! 1. Allocate Ptilde in sparse matrix form
call tmpcoo%allocate(ncol,ntaggr,ncol)
do i=1,ncol
tmpcoo%val(i) = done
tmpcoo%ia(i) = i
tmpcoo%ja(i) = ilaggr(i)
end do
call tmpcoo%set_nzeros(ncol)
call tmpcoo%set_dupl(psb_dupl_add_)
call tmpcoo%set_sorted() ! At this point this is in row-major
call op_prol%mv_to(tmpcoo)
call ptilde%mv_from_coo(tmpcoo,info)
if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (info /= psb_success_) goto 9999
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
mlprec/impl/mld_dcoarse_bld.f90
+173 -147
View File
@@ -98,187 +98,213 @@ subroutine mld_dcoarse_bld(a,desc_a,p,info)
ictxt = desc_a%get_context()
call psb_info(ictxt,me,np)
call mld_check_def(p%parms%ml_type,'Multilevel type',&
& mld_mult_ml_,is_legal_ml_type)
call mld_check_def(p%parms%aggr_alg,'Aggregation',&
& mld_dec_aggr_,is_legal_ml_aggr_alg)
call mld_check_def(p%parms%aggr_ord,'Ordering',&
& mld_aggr_ord_nat_,is_legal_ml_aggr_ord)
call mld_check_def(p%parms%aggr_kind,'Smoother',&
& mld_smooth_prol_,is_legal_ml_aggr_kind)
call mld_check_def(p%parms%coarse_mat,'Coarse matrix',&
& mld_distr_mat_,is_legal_ml_coarse_mat)
call mld_check_def(p%parms%aggr_filter,'Use filtered matrix',&
& mld_no_filter_mat_,is_legal_aggr_filter)
call mld_check_def(p%parms%smoother_pos,'smooth_pos',&
& mld_pre_smooth_,is_legal_ml_smooth_pos)
call mld_check_def(p%parms%aggr_omega_alg,'Omega Alg.',&
& mld_eig_est_,is_legal_ml_aggr_omega_alg)
call mld_check_def(p%parms%aggr_eig,'Eigenvalue estimate',&
& mld_max_norm_,is_legal_ml_aggr_eig)
call mld_check_def(p%parms%aggr_omega_val,'Omega',dzero,is_legal_d_omega)
call mld_check_def(p%parms%aggr_thresh,'Aggr_Thresh',dzero,is_legal_d_aggr_thrs)
select case(p%parms%aggr_alg)
case (mld_dec_aggr_, mld_sym_dec_aggr_)
!
! Build a mapping between the row indices of the fine-level matrix
! and the row indices of the coarse-level matrix, according to a decoupled
! aggregation algorithm. This also defines a tentative prolongator from
! the coarse to the fine level.
!
call mld_aggrmap_bld(p%parms%aggr_alg,p%parms%aggr_ord,p%parms%aggr_thresh,&
& a,desc_a,ilaggr,nlaggr,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmap_bld')
if (.true.) then
call mld_aggrmap_bld(p,&
& a,desc_a,ilaggr,nlaggr,op_prol,info)
else
call mld_check_def(p%parms%ml_type,'Multilevel type',&
& mld_mult_ml_,is_legal_ml_type)
call mld_check_def(p%parms%aggr_alg,'Aggregation',&
& mld_dec_aggr_,is_legal_ml_aggr_alg)
call mld_check_def(p%parms%aggr_ord,'Ordering',&
& mld_aggr_ord_nat_,is_legal_ml_aggr_ord)
call mld_check_def(p%parms%aggr_thresh,'Aggr_Thresh',dzero,is_legal_d_aggr_thrs)
select case(p%parms%aggr_alg)
case (mld_dec_aggr_, mld_sym_dec_aggr_)
!
! Build a mapping between the row indices of the fine-level matrix
! and the row indices of the coarse-level matrix, according to a decoupled
! aggregation algorithm. This also defines a tentative prolongator from
! the coarse to the fine level.
!
call mld_aggrmap_bld(p%parms%aggr_alg,p%parms%aggr_ord,p%parms%aggr_thresh,&
& a,desc_a,ilaggr,nlaggr,op_prol,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmap_bld')
goto 9999
end if
case (mld_bcmatch_aggr_)
write(0,*) 'Matching is not implemented yet '
info = -1111
call psb_errpush(psb_err_input_value_invalid_i_,name,&
& i_err=(/ione,p%parms%aggr_alg,izero,izero,izero/))
goto 9999
end if
case default
info = -1
call psb_errpush(psb_err_input_value_invalid_i_,name,&
& i_err=(/ione,p%parms%aggr_alg,izero,izero,izero/))
goto 9999
end select
end if
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmap_bld')
goto 9999
end if
if (.true.) then
call mld_aggrmat_asb(p,&
& a,desc_a,ilaggr,nlaggr,op_prol,info)
else
call mld_check_def(p%parms%aggr_kind,'Smoother',&
& mld_smooth_prol_,is_legal_ml_aggr_kind)
call mld_check_def(p%parms%coarse_mat,'Coarse matrix',&
& mld_distr_mat_,is_legal_ml_coarse_mat)
call mld_check_def(p%parms%aggr_filter,'Use filtered matrix',&
& mld_no_filter_mat_,is_legal_aggr_filter)
call mld_check_def(p%parms%smoother_pos,'smooth_pos',&
& mld_pre_smooth_,is_legal_ml_smooth_pos)
call mld_check_def(p%parms%aggr_omega_alg,'Omega Alg.',&
& mld_eig_est_,is_legal_ml_aggr_omega_alg)
call mld_check_def(p%parms%aggr_eig,'Eigenvalue estimate',&
& mld_max_norm_,is_legal_ml_aggr_eig)
call mld_check_def(p%parms%aggr_omega_val,'Omega',dzero,is_legal_d_omega)
!
! Build the coarse-level matrix from the fine-level one, starting from
! the mapping defined by mld_aggrmap_bld and applying the aggregation
! algorithm specified by p%iprcparm(mld_aggr_kind_)
!
call mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p%parms,ac,op_prol,op_restr,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_asb')
goto 9999
end if
case (mld_bcmatch_aggr_)
write(0,*) 'Matching is not implemented yet '
info = -1111
call psb_errpush(psb_err_input_value_invalid_i_,name,&
& i_err=(/ione,p%parms%aggr_alg,izero,izero,izero/))
goto 9999
case default
info = -1
call psb_errpush(psb_err_input_value_invalid_i_,name,&
& i_err=(/ione,p%parms%aggr_alg,izero,izero,izero/))
goto 9999
! Common code refactored here.
end select
ntaggr = sum(nlaggr)
! Common code refactored here.
ntaggr = sum(nlaggr)
select case(p%parms%coarse_mat)
select case(p%parms%coarse_mat)
case(mld_distr_mat_)
case(mld_distr_mat_)
call ac%mv_to(bcoo)
if (p%parms%clean_zeros) call bcoo%clean_zeros(info)
nzl = bcoo%get_nzeros()
call ac%mv_to(bcoo)
if (p%parms%clean_zeros) call bcoo%clean_zeros(info)
nzl = bcoo%get_nzeros()
if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1))
if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%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.'
call p%ac%mv_from(bcoo)
if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1))
if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%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.'
call p%ac%mv_from(bcoo)
call p%ac%set_nrows(p%desc_ac%get_local_rows())
call p%ac%set_ncols(p%desc_ac%get_local_cols())
call p%ac%set_asb()
call p%ac%set_nrows(p%desc_ac%get_local_rows())
call p%ac%set_ncols(p%desc_ac%get_local_cols())
call p%ac%set_asb()
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free')
goto 9999
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free')
if (np>1) then
call op_prol%mv_to(acsr1)
nzl = acsr1%get_nzeros()
call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I')
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc')
goto 9999
end if
call op_prol%mv_from(acsr1)
endif
call op_prol%set_ncols(p%desc_ac%get_local_cols())
if (np>1) then
call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_)
call op_restr%mv_to(acoo)
nzl = acoo%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I')
call acoo%set_dupl(psb_dupl_add_)
if (info == psb_success_) call op_restr%mv_from(acoo)
if (info == psb_success_) call op_restr%cscnv(info,type='csr')
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Converting op_restr to local')
goto 9999
end if
end if
!
! Clip to local rows.
!
call op_restr%set_nrows(p%desc_ac%get_local_rows())
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done ac '
case(mld_repl_mat_)
!
!
call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if ((info == psb_success_).and.p%parms%clean_zeros) call ac%clean_zeros(info)
if (info == psb_success_) &
& call psb_gather(p%ac,ac,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if (info /= psb_success_) goto 9999
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
goto 9999
end if
if (np>1) then
call op_prol%mv_to(acsr1)
nzl = acsr1%get_nzeros()
call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I')
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc')
goto 9999
end if
call op_prol%mv_from(acsr1)
endif
call op_prol%set_ncols(p%desc_ac%get_local_cols())
!
! Copy the prolongation/restriction matrices into the descriptor map.
! op_restr => PR^T i.e. restriction operator
! op_prol => PR i.e. prolongation operator
!
if (np>1) then
call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_)
call op_restr%mv_to(acoo)
nzl = acoo%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I')
call acoo%set_dupl(psb_dupl_add_)
if (info == psb_success_) call op_restr%mv_from(acoo)
if (info == psb_success_) call op_restr%cscnv(info,type='csr')
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Converting op_restr to local')
goto 9999
end if
p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999
end if
!
! Clip to local rows.
! Fix the base_a and base_desc pointers for handling of residuals.
! This is correct because this routine is only called at levels >=2.
!
call op_restr%set_nrows(p%desc_ac%get_local_rows())
p%base_a => p%ac
p%base_desc => p%desc_ac
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done ac '
case(mld_repl_mat_)
!
!
call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if ((info == psb_success_).and.p%parms%clean_zeros) call ac%clean_zeros(info)
if (info == psb_success_) &
& call psb_gather(p%ac,ac,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if (info /= psb_success_) goto 9999
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_asb')
goto 9999
end if
!
! Copy the prolongation/restriction matrices into the descriptor map.
! op_restr => PR^T i.e. restriction operator
! op_prol => PR i.e. prolongation operator
!
p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999
end if
!
! Fix the base_a and base_desc pointers for handling of residuals.
! This is correct because this routine is only called at levels >=2.
!
p%base_a => p%ac
p%base_desc => p%desc_ac
call psb_erractionrestore(err_act)
return
mlprec/impl/mld_dprecinit.F90
+1 -1
View File
@@ -171,7 +171,7 @@ subroutine mld_dprecinit(p,ptype,info,nlev)
nlev_ = max(1,nlev)
p%n_prec_levs = nlev_
else
nlev_ = 3
nlev_ = p%max_prec_levs
p%n_prec_levs = -ione
end if
ilev_ = 1
mlprec/mld_base_prec_type.F90
+13 -4
View File
@@ -103,11 +103,12 @@ module mld_base_prec_type
integer(psb_ipk_) :: coarse_mat, coarse_solve
logical :: clean_zeros=.true.
contains
procedure, pass(pm) :: clone => ml_parms_clone
procedure, pass(pm) :: descr => ml_parms_descr
procedure, pass(pm) :: mldescr => ml_parms_mldescr
procedure, pass(pm) :: get_coarse => ml_parms_get_coarse
procedure, pass(pm) :: clone => ml_parms_clone
procedure, pass(pm) :: descr => ml_parms_descr
procedure, pass(pm) :: mldescr => ml_parms_mldescr
procedure, pass(pm) :: coarsedescr => ml_parms_coarsedescr
procedure, pass(pm) :: printout => ml_parms_printout
procedure, pass(pm) :: printout => ml_parms_printout
end type mld_ml_parms
@@ -506,6 +507,14 @@ contains
end select
end function mld_stringval
subroutine ml_parms_get_coarse(pm,pmin)
implicit none
class(mld_ml_parms), intent(inout) :: pm
class(mld_ml_parms), intent(in) :: pmin
pm%coarse_mat = pmin%coarse_mat
pm%coarse_solve = pmin%coarse_solve
end subroutine ml_parms_get_coarse
subroutine ml_parms_printout(pm,iout)
mlprec/mld_d_inner_mod.f90
+32 -17
View File
@@ -124,11 +124,12 @@ module mld_d_inner_mod
end interface mld_bld_mlhier_aggsize
interface mld_bld_mlhier_array
subroutine mld_d_bld_mlhier_array(nplevs,a,desc_a,precv,info)
use psb_base_mod, only : psb_ipk_, psb_dspmat_type, psb_desc_type
subroutine mld_d_bld_mlhier_array(nplevs,casize,mnaggratio,a,desc_a,precv,info)
use psb_base_mod, only : psb_ipk_, psb_dspmat_type, psb_desc_type, psb_dpk_
use mld_d_prec_type, only : mld_d_onelev_type
implicit none
integer(psb_ipk_), intent(inout) :: nplevs
integer(psb_ipk_), intent(inout) :: nplevs, casize
real(psb_dpk_) :: mnaggratio
type(psb_dspmat_type),intent(in), target :: a
type(psb_desc_type), intent(inout), target :: desc_a
type(mld_d_onelev_type), allocatable, target, intent(inout) :: precv(:)
@@ -137,7 +138,18 @@ module mld_d_inner_mod
end interface mld_bld_mlhier_array
interface mld_aggrmap_bld
subroutine mld_daggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info)
subroutine mld_d_lev_aggrmap_bld(p,a,desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_
use mld_d_prec_type, only : mld_d_onelev_type
implicit none
type(mld_d_onelev_type), intent(inout), target :: p
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
type(psb_dspmat_type), intent(out) :: op_prol
integer(psb_ipk_), intent(out) :: info
end subroutine mld_d_lev_aggrmap_bld
subroutine mld_daggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_
implicit none
integer(psb_ipk_), intent(in) :: iorder
@@ -146,6 +158,7 @@ module mld_d_inner_mod
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
type(psb_dspmat_type), intent(out) :: op_prol
integer(psb_ipk_), intent(out) :: info
end subroutine mld_daggrmap_bld
end interface mld_aggrmap_bld
@@ -165,18 +178,19 @@ module mld_d_inner_mod
end interface mld_dec_map_bld
!!$ interface mld_aggrmat_asb
!!$ subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
!!$ use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_
!!$ use mld_d_prec_type, only : mld_d_onelev_type
!!$ implicit none
!!$ type(psb_dspmat_type), intent(in) :: a
!!$ type(psb_desc_type), intent(in) :: desc_a
!!$ integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:)
!!$ type(mld_d_onelev_type), intent(inout), target :: p
!!$ integer(psb_ipk_), intent(out) :: info
!!$ end subroutine mld_daggrmat_asb
!!$ end interface mld_aggrmat_asb
interface mld_aggrmat_asb
subroutine mld_d_lev_aggrmat_asb(p,a,desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_
use mld_d_prec_type, only : mld_d_onelev_type
implicit none
type(mld_d_onelev_type), intent(inout), target :: p
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:),nlaggr(:)
type(psb_dspmat_type), intent(inout) :: op_prol
integer(psb_ipk_), intent(out) :: info
end subroutine mld_d_lev_aggrmat_asb
end interface mld_aggrmat_asb
@@ -189,7 +203,8 @@ module mld_d_inner_mod
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_dml_parms), intent(inout) :: parms
type(psb_dspmat_type), intent(out) :: ac,op_prol,op_restr
type(psb_dspmat_type), intent(inout) :: op_prol
type(psb_dspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
end subroutine mld_daggrmat_var_asb
end interface
tests/pdegen/mld_d_pde3d.f90
+17 -2
View File
@@ -186,6 +186,8 @@ program mld_d_pde3d
type(precdata) :: prectype
type(psb_d_coo_sparse_mat) :: acoo
! other variables
character(len=20) :: dump_prefix
logical :: dump_sol=.false., dump_prec=.false.
integer(psb_ipk_) :: info, i
character(len=20) :: name,ch_err
@@ -214,7 +216,8 @@ program mld_d_pde3d
!
! get parameters
!
call get_parms(ictxt,kmethd,prectype,afmt,idim,istopc,itmax,itrace,irst,eps)
call get_parms(ictxt,kmethd,prectype,afmt,idim,istopc,itmax,itrace,irst,eps,&
&dump_prec,dump_prefix)
!
! allocate and fill in the coefficient matrix, rhs and initial guess
@@ -380,6 +383,10 @@ program mld_d_pde3d
write(psb_out_unit,'("Total memory occupation for PREC: ",i12)') precsize
end if
if (dump_prec) call prec%dump(info,prefix=trim(dump_prefix),&
& ac=.true.,solver=.true.,smoother=.true.,rp=.true.,global_num=.true.)
!
! cleanup storage and exit
!
@@ -404,13 +411,17 @@ contains
!
! get iteration parameters from standard input
!
subroutine get_parms(ictxt,kmethd,prectype,afmt,idim,istopc,itmax,itrace,irst,eps)
subroutine get_parms(ictxt,kmethd,prectype,afmt,idim,istopc,itmax,itrace,irst,eps,&
& dump_prec,dump_prefix)
integer(psb_ipk_) :: ictxt
type(precdata) :: prectype
character(len=*) :: kmethd, afmt
integer(psb_ipk_) :: idim, istopc,itmax,itrace,irst
integer(psb_ipk_) :: np, iam, info
real(psb_dpk_) :: eps
logical :: dump_prec
character(len=*) :: dump_prefix
character(len=20) :: buffer
call psb_info(ictxt, iam, np)
@@ -424,6 +435,8 @@ contains
call read_data(itrace,psb_inp_unit)
call read_data(irst,psb_inp_unit)
call read_data(eps,psb_inp_unit)
call read_data(dump_prec,psb_inp_unit)
call read_data(dump_prefix,psb_inp_unit)
call read_data(prectype%descr,psb_inp_unit) ! verbose description of the prec
call read_data(prectype%prec,psb_inp_unit) ! overall prectype
call read_data(prectype%nlevs,psb_inp_unit) ! Prescribed number of levels
@@ -462,6 +475,8 @@ contains
call psb_bcast(ictxt,itrace)
call psb_bcast(ictxt,irst)
call psb_bcast(ictxt,eps)
call psb_bcast(ictxt,dump_prec)
call psb_bcast(ictxt,dump_prefix)
call psb_bcast(ictxt,prectype%descr) ! verbose description of the prec
call psb_bcast(ictxt,prectype%prec) ! overall prectype
call psb_bcast(ictxt,prectype%nlevs) ! Prescribed number of levels
tests/pdegen/runs/ppde.inp
+12 -10
View File
@@ -6,30 +6,32 @@ CSR ! Storage format CSR COO JAD
10 ! ITRACE
30 ! IRST (restart for RGMRES and BiCGSTABL)
1.d-6 ! EPS
F ! Dump preconditioner on file T F
test-ml-unsm-our ! File prefix for preconditioner dump
ML-MUL-RAS-ILU ! Descriptive name for preconditioner (up to 40 chars)
ML ! Preconditioner NONE JACOBI BJAC AS ML
-1 ! If ML: Prescribed number of levels; if <0, ignore it and use coarse size.
-010 ! If ML: Target coarse size. If <0, then use library default
-4 ! If ML: Prescribed number of levels; if <0, ignore it and use coarse size.
-8000 ! If ML: Target coarse size. If <0, then use library default
-1.5d0 ! If ML: Minimum aggregation ratio; if <0 use library default
-0.10d0 ! If ML: Smoother Aggregation Threshold: >= 0.0 default if <0
-20 ! If ML: Maximum acceptable number of levels; if <0 use library default
SMOOTHED ! Type of aggregation: SMOOTHED, NONSMOOTHED, MINENERGY
SYMDEC ! Type of aggregation: DEC SYMDEC
SMOOTHED ! Type of aggregation: SMOOTHED, UNSMOOTHED, MINENERGY
DEC ! Type of aggregation: DEC SYMDEC
NATURAL ! Ordering of aggregation: NATURAL DEGREE
MULT ! Type of multilevel correction: ADD MULT KCYCLE VCYCLE WCYCLE KCYCLESYM
TWOSIDE ! Side of correction: PRE POST TWOSIDE (ignored for ADD)
4 ! Smoother sweeps
2 ! Smoother sweeps
BJAC ! Smoother type JACOBI BJAC AS; ignored for non-ML
0 ! Number of overlap layers for AS preconditioner (at finest level)
HALO ! AS Restriction operator NONE HALO
NONE ! AS Prolongation operator NONE SUM AVG
GS ! Subdomain solver DSCALE ILU MILU ILUT FWGS BWGS MUMPS UMF SLU
4 ! Solver sweeps for GS
FWGS ! Subdomain solver DSCALE ILU MILU ILUT FWGS BWGS MUMPS UMF SLU
1 ! Solver sweeps for GS
0 ! Level-set N for ILU(N), and P for ILUT
1.d-4 ! Threshold T for ILU(T,P)
DIST ! Coarse level: matrix distribution DIST REPL
BJAC ! Coarse level: solver JACOBI BJAC UMF SLU SLUDIST MUMPS
ILU ! Coarse level: subsolver DSCALE GS BWGS ILU UMF SLU SLUDIST MUMPS
1 ! Coarse level: Level-set N for ILU(N)
FWGS ! Coarse level: subsolver DSCALE GS BWGS ILU UMF SLU SLUDIST MUMPS
0 ! Coarse level: Level-set N for ILU(N)
1.d-4 ! Coarse level: Threshold T for ILU(T,P)
4 ! Coarse level: Number of Jacobi sweeps
2 ! Coarse level: Number of Jacobi sweeps