fdurastante/amg4psblas
1
0
Fork 0
mirror of https://github.com/sfilippone/amg4psblas.git synced 2026-10-20 21:02:35 +00:00
Code Issues Packages Projects Releases Wiki Activity

mld2p4-2:

Browse Source 
config/pac.m4
 configure
 mlprec/impl/Makefile
 mlprec/impl/level/Makefile
 mlprec/impl/mld_d_bld_mlhier_aggsize.f90
 mlprec/impl/mld_d_bld_mlhier_array.f90
 mlprec/impl/mld_dmlprec_bld.f90
 mlprec/mld_d_inner_mod.f90
 mlprec/mld_d_onelev_mod.f90
 tests/pdegen/runs/ppde.inp

Bump minimim GNU compiler version to 4.8.
Unpack mlprec_bld, step 1.
This commit is contained in:
Salvatore Filippone
2016-05-20 13:42:21 +00:00
parent 81a3e58358
commit aa34612da0
10 changed files with 504 additions and 362 deletions
Download Patch File Download Diff File Expand all files Collapse all files
config/pac.m4
+1 -1
View File
@@ -129,7 +129,7 @@ AC_DEFUN(PAC_HAVE_MODERN_GFORTRAN,
ac_fc=${MPIFC-$FC};
AC_COMPILE_IFELSE([
program main
#if ( __GNUC__ >= 4 && __GNUC_MINOR__ >= 6 ) || ( __GNUC__ > 4 )
#if ( __GNUC__ >= 4 && __GNUC_MINOR__ >= 8 ) || ( __GNUC__ > 4 )
print *, "ok"
#else
this program will fail
configure
Vendored
+3 -3
View File
@@ -1393,8 +1393,8 @@ Optional Packages:
--with-extra-libs List additional link flags here. For example,
--with-extra-libs=-lspecial_system_lib or
--with-extra-libs=-L/path/to/libs
--with-mumps=LIBNAME Specify the libname for MUMPS. Default: "-lsmumps
-ldmumps -lcmumps -lzmumps -lmumps_common -lpord"
--with-mumps=LIBNAME Specify the libname for MUMPS. Default: autodetect
with minimum "-lmumps_common -lpord"
--with-mumpsdir=DIR Specify the directory for MUMPS library and
includes. Note: you will need to add auxiliary
libraries with --extra-libs; this depends on how
@@ -5116,7 +5116,7 @@ fi
if test "${with_mumps+set}" = set; then
withval=$with_mumps; mld2p4_cv_mumps=$withval
else
mld2p4_cv_mumps='-lsmumps -ldmumps -lcmumps -lzmumps -lmumps_common -lpord'
mld2p4_cv_mumps='-lmumps_common -lpord'
fi
mlprec/impl/Makefile
+1 -1
View File
@@ -22,7 +22,7 @@ MPFOBJS=$(SMPFOBJS) $(DMPFOBJS) $(CMPFOBJS) $(ZMPFOBJS)
MPCOBJS=mld_sslud_interface.o mld_dslud_interface.o mld_cslud_interface.o mld_zslud_interface.o
DINNEROBJS= mld_dcoarse_bld.o mld_dmlprec_bld.o \
DINNEROBJS= mld_dcoarse_bld.o mld_dmlprec_bld.o mld_d_bld_mlhier_aggsize.o mld_d_bld_mlhier_array.o \
mld_dilu0_fact.o mld_diluk_fact.o mld_dilut_fact.o mld_daggrmap_bld.o \
mld_d_dec_map_bld.o mld_dmlprec_aply.o mld_daggrmat_asb.o \
$(DMPFOBJS)
mlprec/impl/level/Makefile
+1
View File
@@ -21,6 +21,7 @@ mld_c_base_onelev_seti.o \
mld_c_base_onelev_setr.o \
mld_c_base_onelev_setsm.o \
mld_c_base_onelev_setsv.o \
mld_d_base_onelev_build.o \
mld_d_base_onelev_check.o \
mld_d_base_onelev_cnv.o \
mld_d_base_onelev_csetc.o \
mlprec/impl/mld_d_bld_mlhier_aggsize.f90
+235
View File
@@ -0,0 +1,235 @@
!!$
!!$
!!$ 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.
!!$
!!$
!
! Build an aggregation hierarchy with a target aggregation size
!
!
subroutine mld_d_bld_mlhier_aggsize(casize,a,desc_a,iszv,precv,info)
use psb_base_mod
use mld_d_inner_mod, mld_protect_name => mld_d_bld_mlhier_aggsize
use mld_d_prec_mod
implicit none
integer(psb_ipk_), intent(in) :: casize
integer(psb_ipk_), intent(inout) :: iszv
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(:)
integer(psb_ipk_), intent(out) :: info
! Local
integer(psb_ipk_) :: ictxt, me,np
integer(psb_ipk_) :: err,i,k, err_act, newsz
integer(psb_ipk_) :: ipv(mld_ifpsz_), val
integer(psb_ipk_) :: int_err(5)
character :: upd_
class(mld_d_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm
type(mld_dml_parms) :: baseparms, medparms, coarseparms
type(mld_d_onelev_node), pointer :: head, tail, newnode, current
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
name = 'mld_bld_mlhier_aggsize'
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_a%get_ctxt()
call psb_info(ictxt,me,np)
!
! New strategy to build according to coarse size.
!
coarseparms = precv(iszv)%parms
baseparms = precv(1)%parms
medparms = precv(2)%parms
allocate(coarse_sm, source=precv(iszv)%sm,stat=info)
if (info == psb_success_) &
& allocate(med_sm, source=precv(2)%sm,stat=info)
if (info == psb_success_) &
& allocate(base_sm, source=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
!
! Replicated matrix should only ever happen at coarse level.
!
call mld_check_def(baseparms%coarse_mat,'Coarse matrix',&
& mld_distr_mat_,is_distr_ml_coarse_mat)
!
! Now build a doubly linked list
!
allocate(newnode,stat=info)
if (info /= psb_success_) then
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='List start'); goto 9999
end if
head => newnode
tail => newnode
newnode%item%base_a => a
newnode%item%base_desc => desc_a
newnode%item%parms = baseparms
newsz = 1
current => head
list_build_loop: do
allocate(newnode,stat=info)
if (info /= psb_success_) then
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='List start'); goto 9999
end if
current%next => newnode
newnode%prev => current
newsz = newsz + 1
newnode%item%parms = medparms
newnode%item%parms%aggr_thresh = current%item%parms%aggr_thresh/2
call mld_coarse_bld(current%item%base_a, current%item%base_desc, &
& newnode%item,info)
if (info /= psb_success_) then
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='build next level'); goto 9999
end if
if (newsz>2) then
if (all(current%item%map%naggr == newnode%item%map%naggr)) then
!
! We are not gaining anything
!
newsz = newsz-1
current%next => null()
call newnode%item%free(info)
if (info == psb_success_) deallocate(newnode,stat=info)
if (info /= psb_success_) then
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='Deallocate at list end'); goto 9999
end if
exit list_build_loop
end if
end if
current => current%next
tail => current
if (sum(current%item%map%naggr) <= casize) then
!
! Target reached; but we may need to rebuild.
!
exit list_build_loop
end if
end do list_build_loop
!
! At this point, we are at the list tail,
! and it needs to be rebuilt in case the parms were
! different.
!
! But the threshold has to be fixed before rebuliding
coarseparms%aggr_thresh = current%item%parms%aggr_thresh
current%item%parms = coarseparms
call mld_coarse_bld(current%prev%item%base_a,&
& current%prev%item%base_desc, &
& current%item,info)
if (info /= psb_success_) then
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='build next level'); goto 9999
end if
!
! Ok, now allocate the output vector and fix items.
!
do i=1,iszv
if (info == psb_success_) call precv(i)%free(info)
end do
if (info == psb_success_) deallocate(precv,stat=info)
if (info == psb_success_) allocate(precv(newsz),stat=info)
if (info /= psb_success_) then
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='Reallocate precv'); goto 9999
end if
newnode => head
do i=1, newsz
current => newnode
! First do a move_alloc.
! This handles the AC, DESC_AC and MAP fields
if (info == psb_success_) &
& call current%item%move_alloc(precv(i),info)
! Now set the smoother/solver parts.
if (info == psb_success_) then
if (i ==1) then
! This is a workaround for a bug in gfortran 4.7.2
allocate(precv(i)%sm,source=base_sm,stat=info)
else if (i < newsz) then
allocate(precv(i)%sm,source=med_sm,stat=info)
else
allocate(precv(i)%sm,source=coarse_sm,stat=info)
end if
end if
if (info /= psb_success_) then
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='list cpy'); goto 9999
end if
if (i == 1) then
precv(i)%base_a => a
precv(i)%base_desc => desc_a
else
precv(i)%base_a => precv(i)%ac
precv(i)%base_desc => precv(i)%desc_ac
precv(i)%map%p_desc_X => precv(i-1)%base_desc
precv(i)%map%p_desc_Y => precv(i)%base_desc
end if
newnode => current%next
deallocate(current)
end do
call base_sm%free(info)
if (info == psb_success_) call med_sm%free(info)
if (info == psb_success_) call coarse_sm%free(info)
if (info == psb_success_) deallocate(coarse_sm,med_sm,base_sm,stat=info)
if (info /= psb_success_) then
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='final cleanup'); goto 9999
end if
iszv = newsz
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_d_bld_mlhier_aggsize
mlprec/impl/mld_d_bld_mlhier_array.f90
+188
View File
@@ -0,0 +1,188 @@
!!$
!!$
!!$ 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_bld_mlhier_array(a,desc_a,iszv,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) :: iszv
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
integer(psb_ipk_) :: ipv(mld_ifpsz_), val
type(mld_d_onelev_type), allocatable :: tprecv(:)
integer(psb_ipk_) :: int_err(5)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
name = 'mld_bld_array_hierarchy'
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_a%get_ctxt()
call psb_info(ictxt,me,np)
!
!
! Build the matrix and the transfer operators corresponding
! to the remaining levels
!
!
! Check on the iprcparm contents: they should be the same
! on all processes.
!
call psb_bcast(ictxt,precv(1)%parms)
!
! Finest level first; remember to fix base_a and base_desc
!
precv(1)%base_a => a
precv(1)%base_desc => desc_a
iszv = size(precv)
array_build_loop: do i=2, iszv
!
! Check on the iprcparm contents: they should be the same
! on all processes.
!
call psb_bcast(ictxt,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(precv(i)%parms%coarse_mat,'Coarse matrix',&
& mld_distr_mat_,is_distr_ml_coarse_mat)
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Calling mlprcbld at level ',i
!
! Build the mapping between levels i-1 and i and the matrix
! at level i
!
if (info == psb_success_) call mld_coarse_bld(precv(i-1)%base_a,&
& precv(i-1)%base_desc, precv(i),info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Init upper level preconditioner')
goto 9999
endif
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Return from ',i,' call to mlprcbld ',info
if (i>2) then
if (all(precv(i)%map%naggr == precv(i-1)%map%naggr)) then
newsz=i-1
end if
call psb_bcast(ictxt,newsz)
if (newsz > 0) exit array_build_loop
end if
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,&
& a_err='prec reallocation')
goto 9999
endif
do i=1,newsz-1
call precv(i)%move_alloc(tprecv(i),info)
end do
call precv(iszv)%move_alloc(tprecv(newsz),info)
do i=newsz+1, iszv
call precv(i)%free(info)
end do
call move_alloc(tprecv,precv)
! Ignore errors from transfer
info = psb_success_
!
! Restart
iszv = newsz
! Fix the pointers, but the level 1 should
! be already OK
do i=2, iszv - 1
precv(i)%base_a => precv(i)%ac
precv(i)%base_desc => precv(i)%desc_ac
precv(i)%map%p_desc_X => precv(i-1)%base_desc
precv(i)%map%p_desc_Y => precv(i)%base_desc
end do
i = iszv
if (info == psb_success_) call mld_coarse_bld(precv(i-1)%base_a,&
& precv(i-1)%base_desc, precv(i),info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='coarse rebuild')
goto 9999
endif
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_d_bld_mlhier_array
mlprec/impl/mld_dmlprec_bld.f90
+29 -346
View File
@@ -99,9 +99,6 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info,amold,vmold,imold)
integer(psb_ipk_) :: ipv(mld_ifpsz_), val
integer(psb_ipk_) :: int_err(5)
character :: upd_
class(mld_d_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm
type(mld_dml_parms) :: baseparms, medparms, coarseparms
type(mld_d_onelev_node), pointer :: head, tail, newnode, current
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
@@ -124,18 +121,18 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info,amold,vmold,imold)
!
! For the time being we are commenting out the UPDATE argument
! we plan to resurrect it later.
!!$ if (present(upd)) then
!!$ if (debug_level >= psb_debug_outer_) &
!!$ & write(debug_unit,*) me,' ',trim(name),'UPD ', upd
!!$
!!$ if ((psb_toupper(upd).eq.'F').or.(psb_toupper(upd).eq.'T')) then
!!$ upd_=psb_toupper(upd)
!!$ else
!!$ upd_='F'
!!$ endif
!!$ else
!!$ upd_='F'
!!$ endif
! !$ if (present(upd)) then
! !$ if (debug_level >= psb_debug_outer_) &
! !$ & write(debug_unit,*) me,' ',trim(name),'UPD ', upd
! !$
! !$ if ((psb_toupper(upd).eq.'F').or.(psb_toupper(upd).eq.'T')) then
! !$ upd_=psb_toupper(upd)
! !$ else
! !$ upd_='F'
! !$ endif
! !$ else
! !$ upd_='F'
! !$ endif
upd_ = 'F'
if (.not.allocated(p%precv)) then
@@ -148,9 +145,9 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info,amold,vmold,imold)
!
! Check to ensure all procs have the same
!
newsz = -1
newsz = -1
casize = p%coarse_aggr_size
iszv = size(p%precv)
iszv = size(p%precv)
call psb_bcast(ictxt,iszv)
call psb_bcast(ictxt,casize)
if (casize > 0) then
@@ -177,292 +174,20 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info,amold,vmold,imold)
if (casize>0) then
!
! New strategy to build according to coarse size.
if (casize>0) then
!
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
! This should become the default strategy, we specify a target aggregation size.
!
! Replicated matrix should only ever happen at coarse level.
!
call mld_check_def(baseparms%coarse_mat,'Coarse matrix',&
& mld_distr_mat_,is_distr_ml_coarse_mat)
!
! Now build a doubly linked list
!
allocate(newnode,stat=info)
if (info /= psb_success_) then
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='List start'); goto 9999
end if
head => newnode
tail => newnode
newnode%item%base_a => a
newnode%item%base_desc => desc_a
newnode%item%parms = baseparms
newsz = 1
current => head
list_build_loop: do
allocate(newnode,stat=info)
if (info /= psb_success_) then
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='List start'); goto 9999
end if
current%next => newnode
newnode%prev => current
newsz = newsz + 1
newnode%item%parms = medparms
newnode%item%parms%aggr_thresh = current%item%parms%aggr_thresh/2
call mld_coarse_bld(current%item%base_a, current%item%base_desc, &
& newnode%item,info)
if (info /= psb_success_) then
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='build next level'); goto 9999
end if
if (newsz>2) then
if (all(current%item%map%naggr == newnode%item%map%naggr)) then
!
! We are not gaining anything
!
newsz = newsz-1
current%next => null()
call newnode%item%free(info)
if (info == psb_success_) deallocate(newnode,stat=info)
if (info /= psb_success_) then
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='Deallocate at list end'); goto 9999
end if
exit list_build_loop
end if
end if
current => current%next
tail => current
if (sum(newnode%item%map%naggr) <= casize) then
!
! Target reached; but we may need to rebuild.
!
exit list_build_loop
end if
end do list_build_loop
!
! At this point, we are at the list tail,
! and it needs to be rebuilt in case the parms were
! different.
!
! But the threshold has to be fixed before rebuliding
coarseparms%aggr_thresh = current%item%parms%aggr_thresh
current%item%parms = coarseparms
call mld_coarse_bld(current%prev%item%base_a,&
& current%prev%item%base_desc, &
& current%item,info)
if (info /= psb_success_) then
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='build next level'); goto 9999
end if
!
! Ok, now allocate the output vector and fix items.
!
do i=1,iszv
if (info == psb_success_) call p%precv(i)%free(info)
end do
if (info == psb_success_) deallocate(p%precv,stat=info)
if (info == psb_success_) allocate(p%precv(newsz),stat=info)
if (info /= psb_success_) then
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='Reallocate precv'); goto 9999
end if
newnode => head
do i=1, newsz
current => newnode
! First do a move_alloc.
! This handles the AC, DESC_AC and MAP fields
if (info == psb_success_) &
& call current%item%move_alloc(p%precv(i),info)
! Now set the smoother/solver parts.
if (info == psb_success_) then
if (i ==1) then
! This is a workaround for a bug in gfortran 4.7.2
call doallc(i,p%precv,base_sm,info)
! !$ allocate(p%precv(i)%sm,source=base_sm,stat=info)
else if (i < newsz) then
call doallc(i,p%precv,med_sm,info)
! !$ allocate(p%precv(i)%sm,source=med_sm,stat=info)
else
call doallc(i,p%precv,coarse_sm,info)
! !$ allocate(p%precv(i)%sm,source=coarse_sm,stat=info)
end if
end if
if (info /= psb_success_) then
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='list cpy'); goto 9999
end if
if (i == 1) then
p%precv(i)%base_a => a
p%precv(i)%base_desc => desc_a
else
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 if
newnode => current%next
deallocate(current)
end do
call base_sm%free(info)
if (info == psb_success_) call med_sm%free(info)
if (info == psb_success_) call coarse_sm%free(info)
if (info == psb_success_) deallocate(coarse_sm,med_sm,base_sm,stat=info)
if (info /= psb_success_) then
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='final cleanup'); goto 9999
end if
iszv = newsz
call mld_bld_mlhier_aggsize(casize,a,desc_a,iszv,p%precv,info)
else
!
! Default, oldstyle
!
! Oldstyle with fixed number of levels.
!
! Build the matrix and the transfer operators corresponding
! to the remaining levels
!
!
! Check on the iprcparm contents: they should be the same
! on all processes.
!
call psb_bcast(ictxt,p%precv(1)%parms)
!
! Finest level first; remember to fix base_a and base_desc
!
p%precv(1)%base_a => a
p%precv(1)%base_desc => desc_a
do i=2, iszv
!
! Check on the iprcparm contents: they should be the same
! on all processes.
!
call psb_bcast(ictxt,p%precv(i)%parms)
!
! Sanity checks on the parameters
!
if (i<iszv) then
!
! A replicated matrix only makes sense at the coarsest level
!
call mld_check_def(p%precv(i)%parms%coarse_mat,'Coarse matrix',&
& mld_distr_mat_,is_distr_ml_coarse_mat)
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Calling mlprcbld at level ',i
!
! Build the mapping between levels i-1 and i and the matrix
! at level i
!
if (info == psb_success_) call mld_coarse_bld(p%precv(i-1)%base_a,&
& p%precv(i-1)%base_desc, p%precv(i),info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Init upper level preconditioner')
goto 9999
endif
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Return from ',i,' call to mlprcbld ',info
if (i>2) then
if (all(p%precv(i)%map%naggr == p%precv(i-1)%map%naggr)) then
newsz=i-1
end if
call psb_bcast(ictxt,newsz)
if (newsz > 0) exit
end if
end do
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(t_prec%precv(newsz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='prec reallocation')
goto 9999
endif
do i=1,newsz-1
call p%precv(i)%move_alloc(t_prec%precv(i),info)
end do
call p%precv(iszv)%move_alloc(t_prec%precv(newsz),info)
do i=newsz+1, iszv
call p%precv(i)%free(info)
end do
call t_prec%move_alloc(p,info)
! Ignore errors from transfer
info = psb_success_
!
! Restart
iszv = newsz
! Fix the pointers, but the level 1 should
! be already OK
do i=2, iszv - 1
p%precv(i)%base_a => p%precv(i)%ac
p%precv(i)%base_desc => p%precv(i)%desc_ac
p%precv(i)%map%p_desc_X => p%precv(i-1)%base_desc
p%precv(i)%map%p_desc_Y => p%precv(i)%base_desc
end do
i = iszv
if (info == psb_success_) call mld_coarse_bld(p%precv(i-1)%base_a,&
& p%precv(i-1)%base_desc, p%precv(i),info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='coarse rebuild')
goto 9999
endif
end if
call mld_bld_mlhier_array(a,desc_a,p%precv,info)
end if
! Fix nzeros
do i=1,iszv
p%precv(i)%ac_nz_loc = p%precv(i)%ac%get_nzeros()
p%precv(i)%ac_nz_tot = p%precv(i)%ac_nz_loc
call psb_sum(ictxt,p%precv(i)%ac_nz_tot)
end do
!
! The coarse space hierarchy has been built.
!
! Now do the preconditioner build.
!
@@ -471,51 +196,22 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info,amold,vmold,imold)
!
! build the base preconditioner at level i
!
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Calling mlprcbld at level ',i
call mld_check_def(p%precv(i)%parms%sweeps,&
& 'Jacobi sweeps',ione,is_int_positive)
call mld_check_def(p%precv(i)%parms%sweeps_pre,&
& 'Jacobi sweeps',ione,is_int_positive)
call mld_check_def(p%precv(i)%parms%sweeps_post,&
& 'Jacobi sweeps',ione,is_int_positive)
if (.not.allocated(p%precv(i)%sm)) then
!! Error: should have called mld_dprecinit
info=3111
call psb_errpush(info,name)
goto 9999
end if
if (.not.allocated(p%precv(i)%sm%sv)) then
!! Error: should have called mld_dprecinit
info=3111
call psb_errpush(info,name)
goto 9999
end if
call p%precv(i)%sm%build(p%precv(i)%base_a,p%precv(i)%base_desc,&
& 'F',info,amold=amold,vmold=vmold,imold=imold)
if (info == 0) then
if (allocated(p%precv(i)%sm2a)) then
call p%precv(i)%sm2a%build(a,desc_a,upd_,info,&
& amold=amold,vmold=vmold,imold=imold)
p%precv(i)%sm2 => p%precv(i)%sm2a
else
p%precv(i)%sm2 => p%precv(i)%sm
end if
end if
call p%precv(i)%bld(info,amold=amold,vmold=vmold,imold=imold)
if (info /= psb_success_) then
write(ch_err,'(a,i7)') 'Error @ level',i
call psb_errpush(psb_err_internal_error_,name,&
& a_err='One level preconditioner build.')
& a_err=ch_err)
goto 9999
endif
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Return from ',i,' call to mlprcbld ',info
end do
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Exiting with',iszv,' levels'
call psb_erractionrestore(err_act)
return
@@ -523,18 +219,5 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info,amold,vmold,imold)
9999 call psb_error_handler(err_act)
return
contains
subroutine doallc(i,v,src,info)
type(mld_d_onelev_type), intent(inout) :: v(:)
integer(psb_ipk_), intent(in) :: i
class(mld_d_base_smoother_type), intent(in) :: src
integer(psb_ipk_), intent(out) :: info
if ((lbound(v,1)<=i).and.(i<=ubound(v,1))) then
allocate(v(i)%sm,source=src,stat=info)
end if
end subroutine doallc
end subroutine mld_dmlprec_bld
mlprec/mld_d_inner_mod.f90
+26 -4
View File
@@ -64,7 +64,6 @@ module mld_d_inner_mod
end subroutine mld_dmlprec_bld
end interface mld_mlprec_bld
interface mld_mlprec_aply
subroutine mld_dmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_
@@ -95,7 +94,6 @@ module mld_d_inner_mod
end subroutine mld_dmlprec_aply_vect
end interface mld_mlprec_aply
interface mld_coarse_bld
subroutine mld_dcoarse_bld(a,desc_a,p,info)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_
@@ -108,6 +106,32 @@ module mld_d_inner_mod
end subroutine mld_dcoarse_bld
end interface mld_coarse_bld
interface mld_bld_mlhier_aggsize
subroutine mld_d_bld_mlhier_aggsize(casize,a,desc_a,iszv,precv,info)
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(in) :: casize
integer(psb_ipk_), intent(inout) :: iszv
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(:)
integer(psb_ipk_), intent(out) :: info
end subroutine mld_d_bld_mlhier_aggsize
end interface mld_bld_mlhier_aggsize
interface mld_bld_mlhier_array
subroutine mld_d_bld_mlhier_array(a,desc_a,precv,info)
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
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(:)
integer(psb_ipk_), intent(out) :: info
end subroutine mld_d_bld_mlhier_array
end interface mld_bld_mlhier_array
interface mld_aggrmap_bld
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_
@@ -122,7 +146,6 @@ module mld_d_inner_mod
end subroutine mld_daggrmap_bld
end interface mld_aggrmap_bld
interface mld_dec_map_bld
subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_
@@ -166,7 +189,6 @@ module mld_d_inner_mod
end subroutine mld_daggrmat_var_asb
end interface
procedure(mld_daggrmat_var_asb) :: mld_daggrmat_nosmth_asb, &
& mld_daggrmat_smth_asb, mld_daggrmat_minnrg_asb, &
& mld_daggrmat_biz_asb
mlprec/mld_d_onelev_mod.f90
+18 -5
View File
@@ -56,8 +56,8 @@ module mld_d_onelev_mod
use mld_base_prec_type
use mld_d_base_smoother_mod
use psb_base_mod, only : psb_dspmat_type, psb_d_vect_type, &
& psb_d_base_vect_type, psb_dlinmap_type, psb_dpk_, &
use psb_base_mod, only : psb_dspmat_type, psb_d_base_sparse_mat,&
& psb_d_vect_type, psb_d_base_vect_type, psb_dlinmap_type, psb_dpk_, &
& psb_ipk_, psb_long_int_k_, psb_desc_type, psb_i_base_vect_type, &
& psb_erractionsave, psb_error_handler
!
@@ -131,6 +131,7 @@ module mld_d_onelev_mod
type(psb_desc_type), pointer :: base_desc => null()
type(psb_dlinmap_type) :: map
contains
procedure, pass(lv) :: bld => mld_d_base_onelev_build
procedure, pass(lv) :: clone => d_base_onelev_clone
procedure, pass(lv) :: cnv => mld_d_base_onelev_cnv
procedure, pass(lv) :: descr => mld_d_base_onelev_descr
@@ -164,8 +165,20 @@ module mld_d_onelev_mod
& d_base_onelev_nullify, d_base_onelev_get_nzeros, &
& d_base_onelev_clone, d_base_onelev_move_alloc
interface
subroutine mld_d_base_onelev_build(lv,info,amold,vmold,imold)
import :: psb_d_base_sparse_mat, psb_d_base_vect_type, &
& psb_i_base_vect_type, psb_dpk_, mld_d_onelev_type, &
& psb_ipk_, psb_long_int_k_, psb_desc_type
implicit none
class(mld_d_onelev_type), target, intent(inout) :: lv
integer(psb_ipk_), intent(out) :: info
class(psb_d_base_sparse_mat), intent(in), optional :: amold
class(psb_d_base_vect_type), intent(in), optional :: vmold
class(psb_i_base_vect_type), intent(in), optional :: imold
end subroutine mld_d_base_onelev_build
end interface
interface
subroutine mld_d_base_onelev_descr(lv,il,nl,ilmin,info,iout)
import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, &
@@ -482,7 +495,7 @@ contains
subroutine d_base_onelev_move_alloc(lv, b,info)
use psb_base_mod
use psb_base_mod, only : psb_move_alloc
implicit none
class(mld_d_onelev_type), target, intent(inout) :: lv, b
integer(psb_ipk_), intent(out) :: info
tests/pdegen/runs/ppde.inp
+2 -2
View File
@@ -1,6 +1,6 @@
BICGSTAB ! Iterative method: BiCGSTAB BiCG CGS RGMRES BiCGSTABL CG
CSR ! Storage format CSR COO JAD
0040 ! IDIM; domain size is idim**3
0100 ! IDIM; domain size is idim**3
2 ! ISTOPC
2000 ! ITMAX
10 ! ITRACE
@@ -30,4 +30,4 @@ ILU ! Coarse level: subsolver DSCALE GS ILU UMF SLU SLU
1.d-4 ! Coarse level: Threshold T for ILU(T,P)
4 ! Coarse level: Number of Jacobi sweeps
-0.10d0 ! Smoother Aggregation Threshold: >= 0.0 default if <0
-10 ! Coarse size limit to determine levels. If <0, then use NL
1000 ! Coarse size limit to determine levels. If <0, then use NL