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

mld2p4-extaggr:

Browse Source 
mlprec/impl/Makefile
 mlprec/impl/mld_c_dec_map_bld.f90
 mlprec/impl/mld_c_hierarchy_bld.f90
 mlprec/impl/mld_c_lev_aggrmap_bld.f90
 mlprec/impl/mld_c_lev_aggrmat_asb.f90
 mlprec/impl/mld_caggrmap_bld.f90
 mlprec/impl/mld_caggrmat_asb.f90
 mlprec/impl/mld_caggrmat_biz_asb.f90
 mlprec/impl/mld_caggrmat_minnrg_asb.f90
 mlprec/impl/mld_caggrmat_nosmth_asb.f90
 mlprec/impl/mld_caggrmat_smth_asb.f90
 mlprec/impl/mld_d_bld_mlhier_array.f90
 mlprec/impl/mld_d_dec_map_bld.f90.new
 mlprec/impl/mld_d_dec_map_bld.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_smth_asb.f90
 mlprec/impl/mld_dcoarse_bld.f90
 mlprec/impl/mld_dprecinit.F90
 mlprec/impl/mld_s_dec_map_bld.f90
 mlprec/impl/mld_s_hierarchy_bld.f90
 mlprec/impl/mld_s_lev_aggrmap_bld.f90
 mlprec/impl/mld_s_lev_aggrmat_asb.f90
 mlprec/impl/mld_saggrmap_bld.f90
 mlprec/impl/mld_saggrmat_asb.f90
 mlprec/impl/mld_saggrmat_biz_asb.f90
 mlprec/impl/mld_saggrmat_minnrg_asb.f90
 mlprec/impl/mld_saggrmat_nosmth_asb.f90
 mlprec/impl/mld_saggrmat_smth_asb.f90
 mlprec/impl/mld_z_dec_map_bld.f90
 mlprec/impl/mld_z_hierarchy_bld.f90
 mlprec/impl/mld_z_lev_aggrmap_bld.f90
 mlprec/impl/mld_z_lev_aggrmat_asb.f90
 mlprec/impl/mld_zaggrmap_bld.f90
 mlprec/impl/mld_zaggrmat_asb.f90
 mlprec/impl/mld_zaggrmat_biz_asb.f90
 mlprec/impl/mld_zaggrmat_minnrg_asb.f90
 mlprec/impl/mld_zaggrmat_nosmth_asb.f90
 mlprec/impl/mld_zaggrmat_smth_asb.f90
 mlprec/mld_c_inner_mod.f90
 mlprec/mld_d_inner_mod.f90
 mlprec/mld_s_inner_mod.f90
 mlprec/mld_z_inner_mod.f90

New organization of aggregation routines.
This commit is contained in:
Salvatore Filippone
2016-09-30 19:19:23 +00:00
parent 42b7a3770a
commit 74ab1d540f
45 changed files with 2581 additions and 1005 deletions
Download Patch File Download Diff File Expand all files Collapse all files
mlprec/impl/Makefile
+3 -3
View File
@@ -32,19 +32,19 @@ SINNEROBJS= mld_scoarse_bld.o mld_smlprec_bld.o mld_s_bld_mlhier_aggsize.o mld
mld_s_ml_prec_bld.o mld_s_hierarchy_bld.o \
mld_silu0_fact.o mld_siluk_fact.o mld_silut_fact.o mld_saggrmap_bld.o \
mld_s_dec_map_bld.o mld_smlprec_aply.o mld_saggrmat_asb.o \
$(SMPFOBJS) mld_s_extprol_bld.o
$(SMPFOBJS) mld_s_extprol_bld.o mld_s_lev_aggrmap_bld.o mld_s_lev_aggrmat_asb.o
ZINNEROBJS= mld_zcoarse_bld.o mld_zmlprec_bld.o mld_z_bld_mlhier_aggsize.o mld_z_bld_mlhier_array.o \
mld_z_ml_prec_bld.o mld_z_hierarchy_bld.o \
mld_zilu0_fact.o mld_ziluk_fact.o mld_zilut_fact.o mld_zaggrmap_bld.o \
mld_z_dec_map_bld.o mld_zmlprec_aply.o mld_zaggrmat_asb.o \
$(ZMPFOBJS) mld_z_extprol_bld.o
$(ZMPFOBJS) mld_z_extprol_bld.o mld_z_lev_aggrmap_bld.o mld_z_lev_aggrmat_asb.o
CINNEROBJS= mld_ccoarse_bld.o mld_cmlprec_bld.o mld_c_bld_mlhier_aggsize.o mld_c_bld_mlhier_array.o \
mld_c_ml_prec_bld.o mld_c_hierarchy_bld.o \
mld_cilu0_fact.o mld_ciluk_fact.o mld_cilut_fact.o mld_caggrmap_bld.o \
mld_c_dec_map_bld.o mld_cmlprec_aply.o mld_caggrmat_asb.o \
$(CMPFOBJS) mld_c_extprol_bld.o
$(CMPFOBJS) mld_c_extprol_bld.o mld_c_lev_aggrmap_bld.o mld_c_lev_aggrmat_asb.o
INNEROBJS= $(SINNEROBJS) $(DINNEROBJS) $(CINNEROBJS) $(ZINNEROBJS)
mlprec/impl/mld_c_dec_map_bld.f90
+18 -21
View File
@@ -54,12 +54,12 @@ subroutine mld_c_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(:)
complex(psb_spk_), allocatable :: val(:), diag(:)
integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii, ip
type(psb_c_csr_sparse_mat) :: acsr
real(psb_spk_) :: cpling, tcl
logical :: recovery, candidate
logical :: disjoint
integer(psb_ipk_) :: debug_level, debug_unit,err_act
integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: nrow, ncol, n_ne
@@ -138,20 +138,15 @@ subroutine mld_c_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 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 +156,7 @@ subroutine mld_c_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 +164,8 @@ subroutine mld_c_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)
@@ -189,7 +186,7 @@ subroutine mld_c_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
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
& .and. (tmpaggr(j) > 0).and. (abs(val(k)) > cpling)) then
ip = k
cpling = abs(val(k))
end if
@@ -208,7 +205,7 @@ subroutine mld_c_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
step3: do ii=1,nr
i = idxs(ii)
if (ilaggr(i) == -(nr+1)) then
if (ilaggr(i) < 0) then
call a%csget(i,i,nz,irow,icol,val,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
@@ -216,10 +213,10 @@ subroutine mld_c_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 = szero
cpling = dzero
ip = 0
do k=1, nz
j = icol(k)
mlprec/impl/mld_c_hierarchy_bld.f90
+217 -15
View File
@@ -94,9 +94,13 @@ subroutine mld_c_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_spk_) :: mnaggratio
integer(psb_ipk_) :: ipv(mld_ifpsz_), val
integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize, nplevs, mxplevs, iaggsize
real(psb_spk_) :: mnaggratio, sizeratio
class(mld_c_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm
type(mld_sml_parms) :: baseparms, medparms, coarseparms
integer(psb_ipk_), allocatable :: ilaggr(:), nlaggr(:)
type(psb_cspmat_type) :: op_prol
type(mld_c_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_c_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,200 @@ subroutine mld_c_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
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_lev_mat_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
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Internal hierarchy build' )
@@ -220,7 +422,7 @@ subroutine mld_c_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_c_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_ccoarse_bld.f90
!
! Subroutine: mld_ccoarse_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_cspmat_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_c_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_c_lev_aggrmap_bld(p,a,desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod
use mld_c_inner_mod, mld_protect_name => mld_c_lev_aggrmap_bld
implicit none
! Arguments
type(mld_c_onelev_type), intent(inout), target :: p
type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
type(psb_cspmat_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_c_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',szero,is_legal_s_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_c_lev_aggrmap_bld
mlprec/impl/mld_c_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_ccoarse_bld.f90
!
! Subroutine: mld_ccoarse_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_cspmat_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_c_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_c_lev_aggrmat_asb(p,a,desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod
use mld_c_inner_mod, mld_protect_name => mld_c_lev_aggrmat_asb
implicit none
! Arguments
type(mld_c_onelev_type), intent(inout), target :: p
type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:),nlaggr(:)
type(psb_cspmat_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_cspmat_type) :: ac, op_restr
type(psb_c_coo_sparse_mat) :: acoo, bcoo
type(psb_c_csr_sparse_mat) :: acsr1
integer(psb_ipk_) :: nzl, ntaggr
integer(psb_ipk_) :: debug_level, debug_unit
name='mld_ccoarse_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',szero,is_legal_s_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_caggrmat_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_c_lev_aggrmat_asb
mlprec/impl/mld_caggrmap_bld.f90
+26 -3
View File
@@ -79,7 +79,7 @@
! info - integer, output.
! Error code.
!
subroutine mld_caggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info)
subroutine mld_caggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod
use mld_c_inner_mod, mld_protect_name => mld_caggrmap_bld
@@ -93,13 +93,14 @@ subroutine mld_caggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info)
type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
type(psb_cspmat_type), intent(out) :: op_prol
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_cspmat_type) :: atmp, atrans
logical :: recovery
type(psb_c_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,28 @@ subroutine mld_caggrmap_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) = cone
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_caggrmat_asb.f90
+3 -112
View File
@@ -98,7 +98,7 @@
! info - integer, output.
! Error code.
!
subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info)
use psb_base_mod
use mld_c_inner_mod, mld_protect_name => mld_caggrmat_asb
@@ -109,11 +109,12 @@ subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_sml_parms), intent(inout) :: parms
type(mld_c_onelev_type), intent(inout), target :: p
type(psb_cspmat_type), intent(inout) :: ac, op_prol,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
type(psb_cspmat_type) :: ac, op_prol,op_restr
type(psb_c_coo_sparse_mat) :: acoo, bcoo
type(psb_c_csr_sparse_mat) :: acsr1
integer(psb_ipk_) :: nzl,ntaggr, err_act
@@ -165,116 +166,6 @@ subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
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_caggrmat_biz_asb.f90
+4 -21
View File
@@ -89,7 +89,8 @@ subroutine mld_caggrmat_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_sml_parms), intent(inout) :: parms
type(psb_cspmat_type), intent(out) :: ac,op_prol,op_restr
type(psb_cspmat_type), intent(inout) :: op_prol
type(psb_cspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
@@ -128,19 +129,8 @@ subroutine mld_caggrmat_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_caggrmat_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) = cone
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_caggrmat_minnrg_asb.f90
+6 -23
View File
@@ -108,8 +108,9 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_sml_parms), intent(inout) :: parms
type(psb_cspmat_type), intent(out) :: ac,op_prol,op_restr
type(mld_sml_parms), intent(inout) :: parms
type(psb_cspmat_type), intent(inout) :: op_prol
type(psb_cspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
@@ -167,14 +168,6 @@ subroutine mld_caggrmat_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_caggrmat_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) = cone
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
@@ -276,7 +259,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
call am3%mv_to(acsr3)
! Compute omega_int
ommx = cmplx(szero,szero)
ommx = czero
do i=1, ncol
omi(i) = omp(ilaggr(i))
if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
@@ -454,7 +437,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
omp = omp/oden
! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10))
! Compute omega_int
ommx = cmplx(szero,szero)
ommx = czero
do i=1, ncol
omi(i) = omp(ilaggr(i))
if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
mlprec/impl/mld_caggrmat_nosmth_asb.f90
+5 -26
View File
@@ -92,7 +92,8 @@ subroutine mld_caggrmat_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_sml_parms), intent(inout) :: parms
type(psb_cspmat_type), intent(out) :: ac,op_prol,op_restr
type(psb_cspmat_type), intent(inout) :: op_prol
type(psb_cspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
@@ -124,34 +125,12 @@ subroutine mld_caggrmat_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) = cone
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_caggrmat_smth_asb.f90
+5 -18
View File
@@ -104,7 +104,8 @@ subroutine mld_caggrmat_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_sml_parms), intent(inout) :: parms
type(psb_cspmat_type), intent(out) :: ac,op_prol,op_restr
type(psb_cspmat_type), intent(inout) :: op_prol
type(psb_cspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
@@ -147,14 +148,7 @@ subroutine mld_caggrmat_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_caggrmat_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) = cone
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_d_bld_mlhier_array.f90
+17 -37
View File
@@ -37,25 +37,23 @@
!!$
!!$
subroutine mld_d_bld_mlhier_array(nplevs,casize,mnaggratio,a,desc_a,precv,info)
subroutine mld_d_bld_mlhier_array(nplevs,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, casize
real(psb_dpk_) :: mnaggratio
integer(psb_ipk_), intent(inout) :: nplevs
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, iaggsize
integer(psb_ipk_) :: err,i,k, err_act, newsz, iszv
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
@@ -180,45 +178,27 @@ subroutine mld_d_bld_mlhier_array(nplevs,casize,mnaggratio,a,desc_a,precv,info)
& write(debug_unit,*) me,' ',trim(name),&
& 'Return from ',i,' call to mlprcbld ',info
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 (i>2) then
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
+3 -13
View File
@@ -56,10 +56,10 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),&
& 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, nisolate,nthird
integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii, ip
type(psb_d_csr_sparse_mat) :: acsr
real(psb_dpk_) :: cpling, tcl
logical :: recovery, disjoint
logical :: disjoint
integer(psb_ipk_) :: debug_level, debug_unit,err_act
integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: nrow, ncol, n_ne
@@ -114,7 +114,6 @@ 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)
@@ -140,11 +139,6 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
end if
enddo
!!$ 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.
@@ -193,7 +187,6 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
if ((1<=j).and.(j<=nr)) 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
@@ -209,12 +202,10 @@ 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) < 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_
@@ -288,8 +279,7 @@ 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_dec_map_bld.f90.new
+343
View File
@@ -0,0 +1,343 @@
!!$
!!$
!!$ 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.
!!$
!!$
subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod
use mld_d_inner_mod, mld_protect_name => mld_d_dec_map_bld
implicit none
! Arguments
integer(psb_ipk_), intent(in) :: iorder
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
real(psb_dpk_), intent(in) :: theta
integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),&
& ideg(:), idxs(:)
real(psb_dpk_), allocatable :: val(:), diag(:)
integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii
type(psb_d_csr_sparse_mat) :: acsr
real(psb_dpk_) :: cpling, tcl
logical :: recovery
integer(psb_ipk_) :: debug_level, debug_unit,err_act
integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: nrow, ncol, n_ne
character(len=20) :: name, ch_err
if (psb_get_errstatus() /= 0) return
info=psb_success_
name = 'mld_dec_map_bld'
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!
ictxt=desc_a%get_context()
call psb_info(ictxt,me,np)
nrow = desc_a%get_local_rows()
ncol = desc_a%get_local_cols()
nr = a%get_nrows()
allocate(ilaggr(nr),neigh(nr),ideg(nr),idxs(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nr,izero,izero,izero,izero/),&
& a_err='integer')
goto 9999
end if
diag = a%get_diag(info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_sp_getdiag')
goto 9999
end if
if (iorder == mld_aggr_ord_nat_) then
do i=1, nr
ilaggr(i) = -(nr+1)
idxs(i) = i
end do
else
call a%cp_to(acsr)
do i=1, nr
ilaggr(i) = -(nr+1)
ideg(i) = acsr%irp(i+1) - acsr%irp(i)
end do
call acsr%free()
call psb_msort(ideg,ix=idxs,dir=psb_sort_down_)
end if
! Note: -(nr+1) Untouched as yet
! -i 1<=i<=nr Adjacent to aggregate i
! i 1<=i<=nr Belonging to aggregate i
!
! Phase one: group nodes together.
! Very simple minded strategy.
!
naggr = 0
nlp = 0
do
icnt = 0
do ii=1, nr
i = idxs(ii)
if (ilaggr(i) == -(nr+1)) then
!
! 1. Untouched nodes are marked >0 together
! with their neighbours
!
icnt = icnt + 1
naggr = naggr + 1
ilaggr(i) = naggr
call a%csget(i,i,nz,irow,icol,val,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='csget')
goto 9999
end if
do k=1, nz
j = icol(k)
if ((1<=j).and.(j<=nr)) then
ilg = ilaggr(j)
if ((ilg<0).and.(i /= j)) then
if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
ilaggr(j) = naggr
else
ilaggr(j) = -naggr
endif
end if
end if
enddo
!
! 2. Untouched neighbours of these nodes are marked <0.
!
call a%get_neigh(i,neigh,n_ne,info,lev=2_psb_ipk_)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_neigh')
goto 9999
end if
do n = 1, n_ne
m = neigh(n)
if ((1<=m).and.(m<=nr)) then
if (ilaggr(m) == -(nr+1)) ilaggr(m) = -naggr
endif
enddo
endif
enddo
nlp = nlp + 1
if (icnt == 0) exit
enddo
if (debug_level >= psb_debug_outer_) then
write(debug_unit,*) me,' ',trim(name),&
& ' Check 1:',count(ilaggr == -(nr+1))
end if
!
! Phase two: sweep over leftovers.
!
allocate(ils(nr),stat=info)
if(info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/naggr+10,izero,izero,izero,izero/),&
& a_err='integer')
goto 9999
end if
do i=1, size(ils)
ils(i) = 0
end do
do i=1, nr
n = ilaggr(i)
if (n>0) then
if (n>naggr) then
info=psb_err_internal_error_
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_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.
do i=1, nr
if (ilaggr(i) < 0) then
!
! Now some silly rule to break ties:
! Group with adjacent aggregate.
!
isz = nr+1
ia = -1
cpling = dzero
call a%csget(i,i,nz,irow,icol,val,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_sp_getrow')
goto 9999
end if
do j=1, nz
k = icol(j)
if ((1<=k).and.(k<=nr).and.(k /= i)) then
tcl = abs(val(j)) / sqrt(abs(diag(i)*diag(k)))
if (abs(val(j)) > theta*sqrt(abs(diag(i)*diag(k)))) then
!!$ if (tcl > theta) then
n = ilaggr(k)
if (n>0) then
if (n>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?')
goto 9999
end if
if ((abs(val(j))>cpling) .or. &
& ((abs(val(j)) == cpling).and. (ils(n) < isz))) then
!!$ if ((tcl > cpling) .or. ((tcl == cpling).and. (ils(n) < isz))) then
ia = n
isz = ils(n)
cpling = abs(val(j))
!!$ cpling = tcl
endif
endif
endif
end if
enddo
if (ia == -1) then
! At this point, the easiest thing is to start a new aggregate
naggr = naggr + 1
ilaggr(i) = naggr
ils(ilaggr(i)) = 1
else
ilaggr(i) = ia
if (ia>naggr) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr ? ')
goto 9999
end if
ils(ia) = ils(ia) + 1
endif
end if
end do
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(debug_unit,*) me,' ',trim(name),&
& 'Size of aggregate ',i,' :',count(ilaggr == i), ils(i)
enddo
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
if (count(ilaggr<0) >0) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Fatal error: some leftovers')
goto 9999
endif
deallocate(ils,neigh,stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
if (naggr > ncol) then
write(0,*) name,'Error : naggr > ncol',naggr,ncol
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Fatal error: naggr>ncol')
goto 9999
end if
call psb_realloc(ncol,ilaggr,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
allocate(nlaggr(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/np,izero,izero,izero,izero/),&
& a_err='integer')
goto 9999
end if
nlaggr(:) = 0
nlaggr(me+1) = naggr
call psb_sum(ictxt,nlaggr(1:np))
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_d_dec_map_bld
mlprec/impl/mld_daggrmap_bld.f90
+3 -5
View File
@@ -92,8 +92,8 @@ subroutine mld_daggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,op_pro
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(:)
type(psb_dspmat_type), intent(out) :: op_prol
integer(psb_ipk_), intent(out) :: info
! Local variables
@@ -152,14 +152,12 @@ subroutine mld_daggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,op_pro
goto 9999
end if
naggr = nlaggr(me+1)
ntaggr = sum(nlaggr)
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
mlprec/impl/mld_daggrmat_asb.f90
+7 -7
View File
@@ -110,9 +110,9 @@ subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,inf
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_dml_parms), intent(inout) :: parms
!!$ type(mld_d_onelev_type), intent(inout), target :: p
integer(psb_ipk_), intent(out) :: info
type(mld_d_onelev_type), intent(inout), target :: p
type(psb_dspmat_type), intent(inout) :: ac, op_prol,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
type(psb_d_coo_sparse_mat) :: acoo, bcoo
@@ -134,26 +134,26 @@ subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,inf
call psb_info(ictxt, me, np)
select case (parms%aggr_kind)
select case (p%parms%aggr_kind)
case (mld_no_smooth_)
call mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,&
& parms,ac,op_prol,op_restr,info)
& p%parms,ac,op_prol,op_restr,info)
case(mld_smooth_prol_)
call mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr, &
& parms,ac,op_prol,op_restr,info)
& p%parms,ac,op_prol,op_restr,info)
case(mld_biz_prol_)
call mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr, &
& parms,ac,op_prol,op_restr,info)
& p%parms,ac,op_prol,op_restr,info)
case(mld_min_energy_)
call mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr, &
& parms,ac,op_prol,op_restr,info)
& p%parms,ac,op_prol,op_restr,info)
case default
info = psb_err_internal_error_
mlprec/impl/mld_daggrmat_biz_asb.f90
+1 -1
View File
@@ -89,7 +89,7 @@ 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(inout) :: op_prol
type(psb_dspmat_type), intent(inout) :: op_prol
type(psb_dspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
mlprec/impl/mld_daggrmat_minnrg_asb.f90
+3 -3
View File
@@ -108,7 +108,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_dml_parms), intent(inout) :: parms
type(mld_dml_parms), intent(inout) :: parms
type(psb_dspmat_type), intent(inout) :: op_prol
type(psb_dspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
@@ -259,7 +259,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
call am3%mv_to(acsr3)
! Compute omega_int
ommx = cmplx(dzero,dzero)
ommx = dzero
do i=1, ncol
omi(i) = omp(ilaggr(i))
if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
@@ -437,7 +437,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
omp = omp/oden
! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10))
! Compute omega_int
ommx = cmplx(dzero,dzero)
ommx = dzero
do i=1, ncol
omi(i) = omp(ilaggr(i))
if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
mlprec/impl/mld_daggrmat_smth_asb.f90
+1 -1
View File
@@ -148,7 +148,7 @@ 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_)
!
! naggr: number of local aggregates
! nrow: local rows.
!
mlprec/impl/mld_dcoarse_bld.f90
+43 -203
View File
@@ -83,229 +83,69 @@ subroutine mld_dcoarse_bld(a,desc_a,p,info)
integer(psb_mpik_) :: ictxt, np, me
integer(psb_ipk_) :: err_act
integer(psb_ipk_), allocatable :: ilaggr(:), nlaggr(:)
type(psb_dspmat_type) :: ac, op_prol,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_
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_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)
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,info)
!
! 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
end if
if(info /= psb_success_) then
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
! 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
end if
!
! 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_aggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_asb')
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_ = p%max_prec_levs
nlev_ = 3
p%n_prec_levs = -ione
end if
ilev_ = 1
mlprec/impl/mld_s_dec_map_bld.f90
+18 -21
View File
@@ -54,12 +54,12 @@ subroutine mld_s_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_spk_), allocatable :: val(:), diag(:)
integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii, ip
type(psb_s_csr_sparse_mat) :: acsr
real(psb_spk_) :: cpling, tcl
logical :: recovery, candidate
logical :: disjoint
integer(psb_ipk_) :: debug_level, debug_unit,err_act
integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: nrow, ncol, n_ne
@@ -138,20 +138,15 @@ subroutine mld_s_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 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 +156,7 @@ subroutine mld_s_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 +164,8 @@ subroutine mld_s_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)
@@ -189,7 +186,7 @@ subroutine mld_s_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
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
& .and. (tmpaggr(j) > 0).and. (abs(val(k)) > cpling)) then
ip = k
cpling = abs(val(k))
end if
@@ -208,7 +205,7 @@ subroutine mld_s_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
step3: do ii=1,nr
i = idxs(ii)
if (ilaggr(i) == -(nr+1)) then
if (ilaggr(i) < 0) then
call a%csget(i,i,nz,irow,icol,val,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
@@ -216,10 +213,10 @@ subroutine mld_s_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 = szero
cpling = dzero
ip = 0
do k=1, nz
j = icol(k)
mlprec/impl/mld_s_hierarchy_bld.f90
+217 -15
View File
@@ -94,9 +94,13 @@ subroutine mld_s_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_spk_) :: mnaggratio
integer(psb_ipk_) :: ipv(mld_ifpsz_), val
integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize, nplevs, mxplevs, iaggsize
real(psb_spk_) :: mnaggratio, sizeratio
class(mld_s_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm
type(mld_sml_parms) :: baseparms, medparms, coarseparms
integer(psb_ipk_), allocatable :: ilaggr(:), nlaggr(:)
type(psb_sspmat_type) :: op_prol
type(mld_s_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_s_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,200 @@ subroutine mld_s_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
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_lev_mat_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
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Internal hierarchy build' )
@@ -220,7 +422,7 @@ subroutine mld_s_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_s_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_scoarse_bld.f90
!
! Subroutine: mld_scoarse_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_sspmat_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_s_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_s_lev_aggrmap_bld(p,a,desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod
use mld_s_inner_mod, mld_protect_name => mld_s_lev_aggrmap_bld
implicit none
! Arguments
type(mld_s_onelev_type), intent(inout), target :: p
type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
type(psb_sspmat_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_s_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',szero,is_legal_s_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_s_lev_aggrmap_bld
mlprec/impl/mld_s_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_scoarse_bld.f90
!
! Subroutine: mld_scoarse_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_sspmat_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_s_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_s_lev_aggrmat_asb(p,a,desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod
use mld_s_inner_mod, mld_protect_name => mld_s_lev_aggrmat_asb
implicit none
! Arguments
type(mld_s_onelev_type), intent(inout), target :: p
type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:),nlaggr(:)
type(psb_sspmat_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_sspmat_type) :: ac, op_restr
type(psb_s_coo_sparse_mat) :: acoo, bcoo
type(psb_s_csr_sparse_mat) :: acsr1
integer(psb_ipk_) :: nzl, ntaggr
integer(psb_ipk_) :: debug_level, debug_unit
name='mld_scoarse_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',szero,is_legal_s_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_saggrmat_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_s_lev_aggrmat_asb
mlprec/impl/mld_saggrmap_bld.f90
+26 -3
View File
@@ -79,7 +79,7 @@
! info - integer, output.
! Error code.
!
subroutine mld_saggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info)
subroutine mld_saggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod
use mld_s_inner_mod, mld_protect_name => mld_saggrmap_bld
@@ -93,13 +93,14 @@ subroutine mld_saggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info)
type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
type(psb_sspmat_type), intent(out) :: op_prol
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_sspmat_type) :: atmp, atrans
logical :: recovery
type(psb_s_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,28 @@ subroutine mld_saggrmap_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) = sone
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_saggrmat_asb.f90
+3 -112
View File
@@ -98,7 +98,7 @@
! info - integer, output.
! Error code.
!
subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info)
use psb_base_mod
use mld_s_inner_mod, mld_protect_name => mld_saggrmat_asb
@@ -109,11 +109,12 @@ subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_sml_parms), intent(inout) :: parms
type(mld_s_onelev_type), intent(inout), target :: p
type(psb_sspmat_type), intent(inout) :: ac, op_prol,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
type(psb_sspmat_type) :: ac, op_prol,op_restr
type(psb_s_coo_sparse_mat) :: acoo, bcoo
type(psb_s_csr_sparse_mat) :: acsr1
integer(psb_ipk_) :: nzl,ntaggr, err_act
@@ -165,116 +166,6 @@ subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
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_saggrmat_biz_asb.f90
+4 -21
View File
@@ -89,7 +89,8 @@ subroutine mld_saggrmat_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_sml_parms), intent(inout) :: parms
type(psb_sspmat_type), intent(out) :: ac,op_prol,op_restr
type(psb_sspmat_type), intent(inout) :: op_prol
type(psb_sspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
@@ -128,19 +129,8 @@ subroutine mld_saggrmat_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_saggrmat_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) = sone
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_saggrmat_minnrg_asb.f90
+6 -23
View File
@@ -108,8 +108,9 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_sml_parms), intent(inout) :: parms
type(psb_sspmat_type), intent(out) :: ac,op_prol,op_restr
type(mld_sml_parms), intent(inout) :: parms
type(psb_sspmat_type), intent(inout) :: op_prol
type(psb_sspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
@@ -167,14 +168,6 @@ subroutine mld_saggrmat_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_saggrmat_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) = sone
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
@@ -276,7 +259,7 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
call am3%mv_to(acsr3)
! Compute omega_int
ommx = cmplx(szero,szero)
ommx = szero
do i=1, ncol
omi(i) = omp(ilaggr(i))
if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
@@ -454,7 +437,7 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
omp = omp/oden
! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10))
! Compute omega_int
ommx = cmplx(szero,szero)
ommx = szero
do i=1, ncol
omi(i) = omp(ilaggr(i))
if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
mlprec/impl/mld_saggrmat_nosmth_asb.f90
+5 -26
View File
@@ -92,7 +92,8 @@ subroutine mld_saggrmat_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_sml_parms), intent(inout) :: parms
type(psb_sspmat_type), intent(out) :: ac,op_prol,op_restr
type(psb_sspmat_type), intent(inout) :: op_prol
type(psb_sspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
@@ -124,34 +125,12 @@ subroutine mld_saggrmat_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) = sone
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_saggrmat_smth_asb.f90
+5 -18
View File
@@ -104,7 +104,8 @@ subroutine mld_saggrmat_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_sml_parms), intent(inout) :: parms
type(psb_sspmat_type), intent(out) :: ac,op_prol,op_restr
type(psb_sspmat_type), intent(inout) :: op_prol
type(psb_sspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
@@ -147,14 +148,7 @@ subroutine mld_saggrmat_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_saggrmat_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) = sone
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_z_dec_map_bld.f90
+17 -20
View File
@@ -54,12 +54,12 @@ subroutine mld_z_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(:)
complex(psb_dpk_), allocatable :: val(:), diag(:)
integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii, ip
type(psb_z_csr_sparse_mat) :: acsr
real(psb_dpk_) :: cpling, tcl
logical :: recovery, candidate
logical :: disjoint
integer(psb_ipk_) :: debug_level, debug_unit,err_act
integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: nrow, ncol, n_ne
@@ -138,20 +138,15 @@ subroutine mld_z_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 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 +156,7 @@ subroutine mld_z_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 +164,8 @@ subroutine mld_z_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)
@@ -189,7 +186,7 @@ subroutine mld_z_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
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
& .and. (tmpaggr(j) > 0).and. (abs(val(k)) > cpling)) then
ip = k
cpling = abs(val(k))
end if
@@ -208,7 +205,7 @@ subroutine mld_z_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
step3: do ii=1,nr
i = idxs(ii)
if (ilaggr(i) == -(nr+1)) then
if (ilaggr(i) < 0) then
call a%csget(i,i,nz,irow,icol,val,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
@@ -216,8 +213,8 @@ subroutine mld_z_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
mlprec/impl/mld_z_hierarchy_bld.f90
+217 -15
View File
@@ -94,9 +94,13 @@ subroutine mld_z_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_z_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_zspmat_type) :: op_prol
type(mld_z_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_z_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,200 @@ subroutine mld_z_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
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_lev_mat_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
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Internal hierarchy build' )
@@ -220,7 +422,7 @@ subroutine mld_z_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_z_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_zcoarse_bld.f90
!
! Subroutine: mld_zcoarse_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_zspmat_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_z_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_z_lev_aggrmap_bld(p,a,desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod
use mld_z_inner_mod, mld_protect_name => mld_z_lev_aggrmap_bld
implicit none
! Arguments
type(mld_z_onelev_type), intent(inout), target :: p
type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
type(psb_zspmat_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_z_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_z_lev_aggrmap_bld
mlprec/impl/mld_z_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_zcoarse_bld.f90
!
! Subroutine: mld_zcoarse_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_zspmat_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_z_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_z_lev_aggrmat_asb(p,a,desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod
use mld_z_inner_mod, mld_protect_name => mld_z_lev_aggrmat_asb
implicit none
! Arguments
type(mld_z_onelev_type), intent(inout), target :: p
type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:),nlaggr(:)
type(psb_zspmat_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_zspmat_type) :: ac, op_restr
type(psb_z_coo_sparse_mat) :: acoo, bcoo
type(psb_z_csr_sparse_mat) :: acsr1
integer(psb_ipk_) :: nzl, ntaggr
integer(psb_ipk_) :: debug_level, debug_unit
name='mld_zcoarse_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_zaggrmat_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_z_lev_aggrmat_asb
mlprec/impl/mld_zaggrmap_bld.f90
+26 -3
View File
@@ -79,7 +79,7 @@
! info - integer, output.
! Error code.
!
subroutine mld_zaggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info)
subroutine mld_zaggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod
use mld_z_inner_mod, mld_protect_name => mld_zaggrmap_bld
@@ -93,13 +93,14 @@ subroutine mld_zaggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info)
type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
type(psb_zspmat_type), intent(out) :: op_prol
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_zspmat_type) :: atmp, atrans
logical :: recovery
type(psb_z_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,28 @@ subroutine mld_zaggrmap_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) = zone
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_zaggrmat_asb.f90
+3 -112
View File
@@ -98,7 +98,7 @@
! info - integer, output.
! Error code.
!
subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info)
use psb_base_mod
use mld_z_inner_mod, mld_protect_name => mld_zaggrmat_asb
@@ -109,11 +109,12 @@ subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_dml_parms), intent(inout) :: parms
type(mld_z_onelev_type), intent(inout), target :: p
type(psb_zspmat_type), intent(inout) :: ac, op_prol,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
type(psb_zspmat_type) :: ac, op_prol,op_restr
type(psb_z_coo_sparse_mat) :: acoo, bcoo
type(psb_z_csr_sparse_mat) :: acsr1
integer(psb_ipk_) :: nzl,ntaggr, err_act
@@ -165,116 +166,6 @@ subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
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_zaggrmat_biz_asb.f90
+4 -21
View File
@@ -89,7 +89,8 @@ subroutine mld_zaggrmat_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_zspmat_type), intent(out) :: ac,op_prol,op_restr
type(psb_zspmat_type), intent(inout) :: op_prol
type(psb_zspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
@@ -128,19 +129,8 @@ subroutine mld_zaggrmat_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_zaggrmat_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) = zone
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_zaggrmat_minnrg_asb.f90
+6 -23
View File
@@ -108,8 +108,9 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_dml_parms), intent(inout) :: parms
type(psb_zspmat_type), intent(out) :: ac,op_prol,op_restr
type(mld_dml_parms), intent(inout) :: parms
type(psb_zspmat_type), intent(inout) :: op_prol
type(psb_zspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
@@ -167,14 +168,6 @@ subroutine mld_zaggrmat_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_zaggrmat_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) = zone
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
@@ -276,7 +259,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
call am3%mv_to(acsr3)
! Compute omega_int
ommx = cmplx(dzero,dzero)
ommx = zzero
do i=1, ncol
omi(i) = omp(ilaggr(i))
if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
@@ -454,7 +437,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
omp = omp/oden
! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10))
! Compute omega_int
ommx = cmplx(dzero,dzero)
ommx = zzero
do i=1, ncol
omi(i) = omp(ilaggr(i))
if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
mlprec/impl/mld_zaggrmat_nosmth_asb.f90
+5 -26
View File
@@ -92,7 +92,8 @@ subroutine mld_zaggrmat_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_zspmat_type), intent(out) :: ac,op_prol,op_restr
type(psb_zspmat_type), intent(inout) :: op_prol
type(psb_zspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
@@ -124,34 +125,12 @@ subroutine mld_zaggrmat_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) = zone
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_zaggrmat_smth_asb.f90
+5 -18
View File
@@ -104,7 +104,8 @@ subroutine mld_zaggrmat_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_zspmat_type), intent(out) :: ac,op_prol,op_restr
type(psb_zspmat_type), intent(inout) :: op_prol
type(psb_zspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
@@ -147,14 +148,7 @@ subroutine mld_zaggrmat_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_zaggrmat_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) = zone
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/mld_c_inner_mod.f90
+28 -2
View File
@@ -124,11 +124,12 @@ module mld_c_inner_mod
end interface mld_bld_mlhier_aggsize
interface mld_bld_mlhier_array
subroutine mld_c_bld_mlhier_array(nplevs,a,desc_a,precv,info)
subroutine mld_c_bld_mlhier_array(nplevs,casize,mnaggratio,a,desc_a,precv,info)
use psb_base_mod, only : psb_ipk_, psb_cspmat_type, psb_desc_type
use mld_c_prec_type, only : mld_c_onelev_type
implicit none
integer(psb_ipk_), intent(inout) :: nplevs
integer(psb_ipk_), intent(inout) :: nplevs, casize
real(psb_spk_) :: mnaggratio
type(psb_cspmat_type),intent(in), target :: a
type(psb_desc_type), intent(inout), target :: desc_a
type(mld_c_onelev_type), allocatable, target, intent(inout) :: precv(:)
@@ -137,6 +138,17 @@ module mld_c_inner_mod
end interface mld_bld_mlhier_array
interface mld_aggrmap_bld
subroutine mld_c_lev_aggrmap_bld(p,a,desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_, psb_ipk_
use mld_c_prec_type, only : mld_c_onelev_type
implicit none
type(mld_c_onelev_type), intent(inout), target :: p
type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: ilaggr(:),nlaggr(:)
type(psb_cspmat_type), intent(out) :: op_prol
integer(psb_ipk_), intent(out) :: info
end subroutine mld_c_lev_aggrmap_bld
subroutine mld_caggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info)
use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_, psb_ipk_
implicit none
@@ -165,6 +177,20 @@ module mld_c_inner_mod
end interface mld_dec_map_bld
interface mld_lev_mat_asb
subroutine mld_c_lev_aggrmat_asb(p,a,desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_, psb_ipk_
use mld_c_prec_type, only : mld_c_onelev_type
implicit none
type(mld_c_onelev_type), intent(inout), target :: p
type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:),nlaggr(:)
type(psb_cspmat_type), intent(inout) :: op_prol
integer(psb_ipk_), intent(out) :: info
end subroutine mld_c_lev_aggrmat_asb
end interface mld_lev_mat_asb
interface mld_aggrmat_asb
subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_, psb_ipk_
mlprec/mld_d_inner_mod.f90
+21 -9
View File
@@ -125,11 +125,11 @@ module mld_d_inner_mod
interface mld_bld_mlhier_array
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 psb_base_mod, only : psb_ipk_, psb_dspmat_type, psb_desc_type
use mld_d_prec_type, only : mld_d_onelev_type
implicit none
integer(psb_ipk_), intent(inout) :: nplevs, casize
real(psb_dpk_) :: mnaggratio
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(:)
@@ -145,11 +145,11 @@ module mld_d_inner_mod
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(:)
integer(psb_ipk_), 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)
subroutine mld_daggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_
implicit none
integer(psb_ipk_), intent(in) :: iorder
@@ -158,7 +158,6 @@ 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
@@ -192,6 +191,20 @@ module mld_d_inner_mod
end subroutine mld_d_lev_aggrmat_asb
end interface mld_lev_mat_asb
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
abstract interface
subroutine mld_daggrmat_var_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info)
@@ -202,8 +215,7 @@ 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(inout) :: op_prol
type(psb_dspmat_type), intent(out) :: ac,op_restr
type(psb_dspmat_type), intent(out) :: ac,op_prol,op_restr
integer(psb_ipk_), intent(out) :: info
end subroutine mld_daggrmat_var_asb
end interface
@@ -211,7 +223,7 @@ module mld_d_inner_mod
procedure(mld_daggrmat_var_asb) :: mld_daggrmat_nosmth_asb, &
& mld_daggrmat_smth_asb, mld_daggrmat_minnrg_asb, &
& mld_daggrmat_biz_asb, mld_daggrmat_asb
& mld_daggrmat_biz_asb
end module mld_d_inner_mod
mlprec/mld_s_inner_mod.f90
+28 -2
View File
@@ -124,11 +124,12 @@ module mld_s_inner_mod
end interface mld_bld_mlhier_aggsize
interface mld_bld_mlhier_array
subroutine mld_s_bld_mlhier_array(nplevs,a,desc_a,precv,info)
subroutine mld_s_bld_mlhier_array(nplevs,casize,mnaggratio,a,desc_a,precv,info)
use psb_base_mod, only : psb_ipk_, psb_sspmat_type, psb_desc_type
use mld_s_prec_type, only : mld_s_onelev_type
implicit none
integer(psb_ipk_), intent(inout) :: nplevs
integer(psb_ipk_), intent(inout) :: nplevs, casize
real(psb_spk_) :: mnaggratio
type(psb_sspmat_type),intent(in), target :: a
type(psb_desc_type), intent(inout), target :: desc_a
type(mld_s_onelev_type), allocatable, target, intent(inout) :: precv(:)
@@ -137,6 +138,17 @@ module mld_s_inner_mod
end interface mld_bld_mlhier_array
interface mld_aggrmap_bld
subroutine mld_s_lev_aggrmap_bld(p,a,desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_, psb_ipk_
use mld_s_prec_type, only : mld_s_onelev_type
implicit none
type(mld_s_onelev_type), intent(inout), target :: p
type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: ilaggr(:),nlaggr(:)
type(psb_sspmat_type), intent(out) :: op_prol
integer(psb_ipk_), intent(out) :: info
end subroutine mld_s_lev_aggrmap_bld
subroutine mld_saggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info)
use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_, psb_ipk_
implicit none
@@ -165,6 +177,20 @@ module mld_s_inner_mod
end interface mld_dec_map_bld
interface mld_lev_mat_asb
subroutine mld_s_lev_aggrmat_asb(p,a,desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_, psb_ipk_
use mld_s_prec_type, only : mld_s_onelev_type
implicit none
type(mld_s_onelev_type), intent(inout), target :: p
type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:),nlaggr(:)
type(psb_sspmat_type), intent(inout) :: op_prol
integer(psb_ipk_), intent(out) :: info
end subroutine mld_s_lev_aggrmat_asb
end interface mld_lev_mat_asb
interface mld_aggrmat_asb
subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_, psb_ipk_
mlprec/mld_z_inner_mod.f90
+28 -2
View File
@@ -124,11 +124,12 @@ module mld_z_inner_mod
end interface mld_bld_mlhier_aggsize
interface mld_bld_mlhier_array
subroutine mld_z_bld_mlhier_array(nplevs,a,desc_a,precv,info)
subroutine mld_z_bld_mlhier_array(nplevs,casize,mnaggratio,a,desc_a,precv,info)
use psb_base_mod, only : psb_ipk_, psb_zspmat_type, psb_desc_type
use mld_z_prec_type, only : mld_z_onelev_type
implicit none
integer(psb_ipk_), intent(inout) :: nplevs
integer(psb_ipk_), intent(inout) :: nplevs, casize
real(psb_dpk_) :: mnaggratio
type(psb_zspmat_type),intent(in), target :: a
type(psb_desc_type), intent(inout), target :: desc_a
type(mld_z_onelev_type), allocatable, target, intent(inout) :: precv(:)
@@ -137,6 +138,17 @@ module mld_z_inner_mod
end interface mld_bld_mlhier_array
interface mld_aggrmap_bld
subroutine mld_z_lev_aggrmap_bld(p,a,desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_, psb_ipk_
use mld_z_prec_type, only : mld_z_onelev_type
implicit none
type(mld_z_onelev_type), intent(inout), target :: p
type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: ilaggr(:),nlaggr(:)
type(psb_zspmat_type), intent(out) :: op_prol
integer(psb_ipk_), intent(out) :: info
end subroutine mld_z_lev_aggrmap_bld
subroutine mld_zaggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info)
use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_, psb_ipk_
implicit none
@@ -165,6 +177,20 @@ module mld_z_inner_mod
end interface mld_dec_map_bld
interface mld_lev_mat_asb
subroutine mld_z_lev_aggrmat_asb(p,a,desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_, psb_ipk_
use mld_z_prec_type, only : mld_z_onelev_type
implicit none
type(mld_z_onelev_type), intent(inout), target :: p
type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:),nlaggr(:)
type(psb_zspmat_type), intent(inout) :: op_prol
integer(psb_ipk_), intent(out) :: info
end subroutine mld_z_lev_aggrmat_asb
end interface mld_lev_mat_asb
interface mld_aggrmat_asb
subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_, psb_ipk_