Merge branch 'MixedI8' into mergemixed. Fixed configure for MUMPS vs LPK8

stopcriterion
Salvatore Filippone 6 years ago
commit b8900beee3

@ -414,6 +414,7 @@ LDFLAGS=" $PSBLAS_LIBS $save_LDFLAGS"
dnl ac_compile='${MPIFC-$FC} -c -o conftest${ac_objext} $FMFLAG$PSBLAS_DIR/include $FMFLAG$PSBLAS_DIR/lib conftest.$ac_ext 1>&5'
dnl ac_link='${MPIFC-$FC} -o conftest${ac_exeext} $FCFLAGS $LDFLAGS conftest.$ac_ext $FMFLAG$PSBLAS_DIR/include -L$PSBLAS_DIR/lib -lpsb_base $LIBS 1>&5'
ac_link='${MPIFC-$FC} -o conftest${ac_exeext} $FCFLAGS conftest.$ac_ext $LDFLAGS $LIBS 1>&5'
dnl Warning : square brackets are EVIL!
AC_LINK_IFELSE([
@ -439,6 +440,20 @@ AC_LINK_IFELSE([
end program test],
[pac_cv_psblas_patchlevel=`./conftest${ac_exeext} | sed 's/^ *//'`],
[pac_cv_psblas_patchlevel="unknown"])
AC_LINK_IFELSE([
program test
use psb_base_mod, only : psb_sizeof_ip
print *,psb_sizeof_ip
end program test],
[pac_cv_psblas_sizeof_ip=`./conftest${ac_exeext} | sed 's/^ *//'`],
[pac_cv_psblas_sizeof_ip="unknown"])
AC_LINK_IFELSE([
program test
use psb_base_mod, only : psb_sizeof_lp
print *,psb_sizeof_lp
end program test],
[pac_cv_psblas_sizeof_lp=`./conftest${ac_exeext} | sed 's/^ *//'`],
[pac_cv_psblas_sizeof_lp="unknown"])
LDFLAGS="$save_LDFLAGS";
FCFLAGS="$save_FCFLAGS";
@ -635,7 +650,7 @@ if test "x$pac_slu_header_ok" == "xyes" ; then
LIBS="$SLU_LIBS -lm $save_LIBS";
AC_TRY_LINK_FUNC(superlu_malloc,
[mld2p4_cv_have_superlu=yes;pac_slu_lib_ok=yes;],
[mld2p4_cv_have_superlu=no;pac_slu_lib_ok=no; SLU_LIBS=""; ])
[mld2p4_cv_have_superlu=no;pac_slu_lib_ok=no; SLU_LIBS="";])
fi
if test "x$pac_slu_lib_ok" == "xno" ; then
dnl Maybe lib64?

9671
configure vendored

File diff suppressed because it is too large Load Diff

@ -495,6 +495,7 @@ if (( $pac_cv_psblas_major < 3 )) ||
else
AC_MSG_NOTICE([Am configuring with PSBLAS version $pac_cv_psblas_major.$pac_cv_psblas_minor.$pac_cv_psblas_patchlevel.])
fi
AC_MSG_NOTICE([PSBLAS size of LPK "$pac_cv_psblas_sizeof_lp".])
###############################################################################
# Parachute rules for ar and ranlib ... (could cause problems)
###############################################################################
@ -629,7 +630,12 @@ AC_LANG([C])
PAC_CHECK_MUMPS
if test "x$mld2p4_cv_have_mumps" == "xyes" ; then
if test "x$pac_cv_psblas_sizeof_lp" == "x8" ; then
AC_MSG_NOTICE([PSBLAS defines PSB_LPK_ as $pac_cv_psblas_sizeof_lp. We cannot interface MUMPS. ])
MUMPS_FLAGS=""
MUMPS_LIBS=""
mld2p4_cv_have_mumps="no";
elif test "x$mld2p4_cv_have_mumps" == "xyes" ; then
MUMPS_FLAGS="-DHave_MUMPS_ $MUMPS_INCLUDES"
if test "x$pac_mumps_fmods_ok" == "xyes" ; then
FDEFINES="$mld_cv_define_prepend-DHAVE_MUMPS_ $mld_cv_define_prepend-DHAVE_MUMPS_MODULES_ $MUMPS_MODULES $FDEFINES"

@ -81,7 +81,7 @@ program mld_cexample_1lev
! other variables
integer :: i,info,j,m_problem
integer(psb_long_int_k_) :: amatsize, precsize, descsize
integer(psb_epk_) :: amatsize, precsize, descsize
integer :: ierr, ircode
real(psb_spk_) :: resmx, resmxp
real(psb_dpk_) :: t1, t2, tprec

@ -96,7 +96,7 @@ program mld_cexample_ml
! other variables
integer :: choice
integer :: i,info,j,m_problem
integer(psb_long_int_k_) :: amatsize, precsize, descsize
integer(psb_epk_) :: amatsize, precsize, descsize
integer :: ierr, ircode
real(psb_spk_) :: resmx, resmxp
real(psb_dpk_) :: t1, t2, tprec

@ -81,7 +81,7 @@ program mld_dexample_1lev
! other variables
integer :: i,info,j,m_problem
integer(psb_long_int_k_) :: amatsize, precsize, descsize
integer(psb_epk_) :: amatsize, precsize, descsize
integer :: ierr, ircode
real(psb_dpk_) :: resmx, resmxp
real(psb_dpk_) :: t1, t2, tprec

@ -96,7 +96,7 @@ program mld_dexample_ml
! other variables
integer :: choice
integer :: i,info,j,m_problem
integer(psb_long_int_k_) :: amatsize, precsize, descsize
integer(psb_epk_) :: amatsize, precsize, descsize
integer :: ierr, ircode
real(psb_dpk_) :: resmx, resmxp
real(psb_dpk_) :: t1, t2, tprec

@ -81,7 +81,7 @@ program mld_sexample_1lev
! other variables
integer :: i,info,j,m_problem
integer(psb_long_int_k_) :: amatsize, precsize, descsize
integer(psb_epk_) :: amatsize, precsize, descsize
integer :: ierr, ircode
real(psb_spk_) :: resmx, resmxp
real(psb_dpk_) :: t1, t2, tprec

@ -96,7 +96,7 @@ program mld_sexample_ml
! other variables
integer :: choice
integer :: i,info,j,m_problem
integer(psb_long_int_k_) :: amatsize, precsize, descsize
integer(psb_epk_) :: amatsize, precsize, descsize
integer :: ierr, ircode
real(psb_spk_) :: resmx, resmxp
real(psb_dpk_) :: t1, t2, tprec

@ -81,7 +81,7 @@ program mld_zexample_1lev
! other variables
integer :: i,info,j,m_problem
integer(psb_long_int_k_) :: amatsize, precsize, descsize
integer(psb_epk_) :: amatsize, precsize, descsize
integer :: ierr, ircode
real(psb_dpk_) :: resmx, resmxp
real(psb_dpk_) :: t1, t2, tprec

@ -96,7 +96,7 @@ program mld_zexample_ml
! other variables
integer :: choice
integer :: i,info,j,m_problem
integer(psb_long_int_k_) :: amatsize, precsize, descsize
integer(psb_epk_) :: amatsize, precsize, descsize
integer :: ierr, ircode
real(psb_dpk_) :: resmx, resmxp
real(psb_dpk_) :: t1, t2, tprec

@ -89,7 +89,7 @@ program mld_dexample_1lev
! other variables
integer :: i,info,j
integer(psb_long_int_k_) :: amatsize, precsize, descsize
integer(psb_epk_) :: amatsize, precsize, descsize
integer :: idim, nlev, ierr, ircode
real(psb_dpk_) :: resmx, resmxp
real(psb_dpk_) :: t1, t2, tprec

@ -108,7 +108,7 @@ program mld_dexample_ml
! other variables
integer :: choice
integer :: i,info,j
integer(psb_long_int_k_) :: amatsize, precsize, descsize
integer(psb_epk_) :: amatsize, precsize, descsize
integer :: idim, ierr, ircode
real(psb_dpk_) :: resmx, resmxp
real(psb_dpk_) :: t1, t2, tprec

@ -89,7 +89,7 @@ program mld_sexample_1lev
! other variables
integer :: i,info,j
integer(psb_long_int_k_) :: amatsize, precsize, descsize
integer(psb_epk_) :: amatsize, precsize, descsize
integer :: idim, nlev, ierr, ircode
real(psb_spk_) :: resmx, resmxp
real(psb_dpk_) :: t1, t2, tprec

@ -108,7 +108,7 @@ program mld_sexample_ml
! other variables
integer :: choice
integer :: i,info,j
integer(psb_long_int_k_) :: amatsize, precsize, descsize
integer(psb_epk_) :: amatsize, precsize, descsize
integer :: idim, ierr, ircode
real(psb_spk_) :: resmx, resmxp
real(psb_dpk_) :: t1, t2, tprec

@ -143,23 +143,25 @@ subroutine mld_c_dec_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,ac,op_p
type(mld_sml_parms), intent(inout) :: parms
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
type(psb_cspmat_type), intent(out) :: ac,op_restr
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(psb_lcspmat_type), intent(inout) :: op_prol
type(psb_lcspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
character(len=20) :: name
integer(psb_mpik_) :: ictxt, np, me
type(psb_c_coo_sparse_mat) :: acoo, bcoo
type(psb_c_csr_sparse_mat) :: acsr1
integer(psb_ipk_) :: nzl,ntaggr
integer(psb_ipk_) :: err_act
character(len=20) :: name
integer(psb_mpk_) :: ictxt, np, me
type(psb_lc_coo_sparse_mat) :: acoo, bcoo
type(psb_lc_csr_sparse_mat) :: acsr1
integer(psb_lpk_) :: nzl,ntaggr
integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: debug_level, debug_unit
name='mld_c_dec_aggregator_mat_asb'
if (psb_get_errstatus().ne.0) return
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
info = psb_success_

@ -81,20 +81,22 @@ subroutine mld_c_dec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_
type(mld_sml_parms), intent(inout) :: parms
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_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
type(psb_lcspmat_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_) :: ntaggr
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name
integer(psb_mpk_) :: ictxt, np, me
integer(psb_ipk_) :: err_act
integer(psb_lpk_) :: ntaggr
integer(psb_ipk_) :: debug_level, debug_unit
name='mld_c_dec_aggregator_tprol'
if (psb_get_errstatus().ne.0) return
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
info = psb_success_

@ -87,22 +87,24 @@ subroutine mld_c_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info)
! Arguments
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), allocatable, intent(inout) :: ilaggr(:),nlaggr(:)
type(psb_cspmat_type), intent(out) :: op_prol
integer(psb_lpk_), allocatable, intent(inout) :: ilaggr(:),nlaggr(:)
type(psb_lcspmat_type), intent(out) :: op_prol
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m,naggrm1, naggrp1, ntaggr
type(psb_c_coo_sparse_mat) :: tmpcoo
integer(psb_lpk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m,naggrm1, naggrp1, ntaggr
type(psb_lc_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
integer(psb_lpk_) :: nrow, ncol, n_ne
character(len=20) :: name, ch_err
if(psb_get_errstatus() /= 0) return
info=psb_success_
name = 'mld_map_to_tprol'
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!

@ -80,12 +80,13 @@ subroutine mld_c_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
real(psb_spk_), intent(in) :: theta
integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),&
& ideg(:), idxs(:), tmpaggr(:)
integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),&
& ideg(:), idxs(:)
integer(psb_lpk_), allocatable :: tmpaggr(:)
complex(psb_spk_), allocatable :: val(:), diag(:)
integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, nc, naggr,i,j,m, nz, ilg, ii, ip
type(psb_c_csr_sparse_mat) :: acsr
@ -96,10 +97,12 @@ subroutine mld_c_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
integer(psb_ipk_) :: nrow, ncol, n_ne
character(len=20) :: name, ch_err
if (psb_get_errstatus() /= 0) return
info=psb_success_
name = 'mld_soc1_map_bld'
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!

@ -79,12 +79,13 @@ subroutine mld_c_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
real(psb_spk_), intent(in) :: theta
integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),&
& ideg(:), idxs(:), tmpaggr(:)
& ideg(:), idxs(:)
integer(psb_lpk_), allocatable :: tmpaggr(:)
complex(psb_spk_), allocatable :: val(:), diag(:)
integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr,nc,naggr,i,j,m, nz, ilg, ii, ip, ip1,nzcnt
type(psb_c_csr_sparse_mat) :: acsr, muij, s_neigh
@ -96,10 +97,12 @@ subroutine mld_c_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
integer(psb_ipk_) :: nrow, ncol, n_ne
character(len=20) :: name, ch_err
if (psb_get_errstatus() /= 0) return
info=psb_success_
name = 'mld_soc2_map_bld'
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!

@ -85,21 +85,24 @@ subroutine mld_c_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,
type(mld_sml_parms), intent(inout) :: parms
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_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
type(psb_lcspmat_type), intent(out) :: op_prol
integer(psb_ipk_), intent(out) :: info
! Local variables
type(psb_cspmat_type) :: atmp, atrans
character(len=20) :: name
integer(psb_mpik_) :: ictxt, np, me
integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: ntaggr, nr
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name
integer(psb_mpk_) :: ictxt, np, me
integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: nr
integer(psb_lpk_) :: ntaggr
integer(psb_ipk_) :: debug_level, debug_unit
name='mld_c_symdec_aggregator_tprol'
if (psb_get_errstatus().ne.0) return
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
info = psb_success_

@ -88,33 +88,35 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
implicit none
! Arguments
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(inout) :: op_prol
type(psb_cspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_sml_parms), intent(inout) :: parms
type(psb_lcspmat_type), intent(inout) :: op_prol
type(psb_lcspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_) :: nrow, nglob, ncol, ntaggr, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw, err_act
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, ip, &
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw
integer(psb_ipk_) ::ictxt, np, me
character(len=20) :: name
type(psb_cspmat_type) :: am3, am4,tmp_prol
type(psb_c_coo_sparse_mat) :: tmpcoo
type(psb_c_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
type(psb_lcspmat_type) :: am3, am4,tmp_prol, la
type(psb_lc_coo_sparse_mat) :: tmpcoo
type(psb_lc_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
complex(psb_spk_), allocatable :: adiag(:)
integer(psb_ipk_) :: ierr(5)
logical :: filter_mat
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: debug_level, debug_unit, err_act
integer(psb_ipk_), parameter :: ncmax=16
real(psb_spk_) :: anorm, omega, tmp, dg, theta
name='mld_aggrmat_biz_asb'
if(psb_get_errstatus().ne.0) return
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
@ -128,10 +130,8 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
ncol = desc_a%get_local_cols()
theta = parms%aggr_thresh
naggr = nlaggr(me+1)
ntaggr = sum(nlaggr)
filter_mat = (parms%aggr_filter == mld_filter_mat_)
! naggr: number of local aggregates
@ -143,8 +143,8 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
& call psb_realloc(ncol,adiag,info)
if (info == psb_success_) &
& call psb_halo(adiag,desc_a,info)
if(info /= psb_success_) then
if (info == psb_success_) call a%cp_to_l(la)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
goto 9999
end if
@ -152,7 +152,7 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
! 1. Allocate Ptilde in sparse matrix form
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_) call la%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (info /= psb_success_) goto 9999
if (debug_level >= psb_debug_outer_) &
@ -163,9 +163,9 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
!
! Build the filtered matrix Af from A
!
if (info == psb_success_) call a%cscnv(acsrf,info,dupl=psb_dupl_add_)
if (info == psb_success_) call acsr3%cp_to_fmt(acsrf,info)
do i=1,nrow
do i=1, nrow
tmp = czero
jd = -1
do j=acsrf%irp(i),acsrf%irp(i+1)-1
@ -271,17 +271,15 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
! Doing it this way means to consider diag(Af_i)
!
!
call psb_symbmm(acsrf,ptilde,acsr1,info)
call psb_spspmm(acsrf,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
goto 9999
end if
call psb_numbmm(acsrf,ptilde,acsr1)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
& 'Done SPSPMM 1'
else
!
@ -307,17 +305,15 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
! Doing it this way means to consider diag(A_i)
!
!
call psb_symbmm(acsr3,ptilde,acsr1,info)
call psb_spspmm(acsr3,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
goto 9999
end if
call psb_numbmm(acsr3,ptilde,acsr1)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
& 'Done SPSPMM 1'
end if
call ptilde%free()
@ -331,16 +327,15 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
goto 9999
end if
call psb_symbmm(a,tmp_prol,am3,info)
call psb_spspmm(la,tmp_prol,am3,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2')
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 2')
goto 9999
end if
call psb_numbmm(a,tmp_prol,am3)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2',parms%aggr_prol, mld_smooth_prol_
& 'Done SPSPMM 2',parms%aggr_prol, mld_smooth_prol_
call tmp_prol%transp(op_restr)
if (debug_level >= psb_debug_outer_) &
@ -353,23 +348,18 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting symbmm 3'
call psb_symbmm(op_restr,am3,ac,info)
if (info == psb_success_) call psb_numbmm(op_restr,am3,ac)
& 'starting spspmm 3'
call psb_spspmm(op_restr,am3,ac,info)
if (info == psb_success_) call am3%free()
if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_)
if (info == psb_success_) call ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Build b = op_restr x am3')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done smooth_aggregate '

@ -114,37 +114,38 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
! Arguments
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(inout) :: op_prol
type(psb_cspmat_type), intent(out) :: ac,op_restr
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_sml_parms), intent(inout) :: parms
type(psb_lcspmat_type), intent(inout) :: op_prol
type(psb_lcspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_), allocatable :: nzbr(:), idisp(:)
integer(psb_ipk_) :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrt, err_act
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrt
integer(psb_ipk_) :: ictxt,np,me, icomm
character(len=20) :: name
type(psb_cspmat_type) :: af, ptilde, rtilde, atran, atp, atdatp
type(psb_cspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da
type(psb_cspmat_type) :: dat, datp, datdatp, atmp3, tmp_prol
type(psb_c_coo_sparse_mat) :: tmpcoo
type(psb_c_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf
type(psb_c_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc
type(psb_lcspmat_type) :: la, af, ptilde, rtilde, atran, atp, atdatp
type(psb_lcspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da
type(psb_lcspmat_type) :: dat, datp, datdatp, atmp3, tmp_prol
type(psb_lc_coo_sparse_mat) :: tmpcoo
type(psb_lc_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf
type(psb_lc_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc
complex(psb_spk_), allocatable :: adiag(:), adinv(:)
complex(psb_spk_), allocatable :: omf(:), omp(:), omi(:), oden(:)
logical :: filter_mat
integer(psb_ipk_) :: ierr(5)
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: debug_level, debug_unit, err_act
integer(psb_ipk_), parameter :: ncmax=16
real(psb_spk_) :: anorm, theta
complex(psb_spk_) :: tmp, alpha, beta, ommx
name='mld_aggrmat_minnrg'
if(psb_get_errstatus().ne.0) return
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
@ -163,13 +164,6 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
naggr = nlaggr(me+1)
ntaggr = sum(nlaggr)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_; ierr(1)=2*np;
call psb_errpush(info,name,i_err=ierr,a_err='integer')
goto 9999
end if
naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1))
@ -193,6 +187,11 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
& call psb_realloc(ncol,adiag,info)
if (info == psb_success_) &
& call psb_halo(adiag,desc_a,info)
if (info == psb_success_) call a%cp_to_l(la)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
goto 9999
end if
do i=1,size(adiag)
if (adiag(i) /= czero) then
@ -202,10 +201,6 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
end if
end do
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
goto 9999
end if
! 1. Allocate Ptilde in sparse matrix form
@ -213,8 +208,8 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
call ptilde%mv_from(tmpcoo)
call ptilde%cscnv(info,type='csr')
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_) call la%cscnv(am3,info,type='csr',dupl=psb_dupl_add_)
if (info == psb_success_) call la%cscnv(da,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
@ -225,11 +220,10 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
call da%scal(adinv,info)
call psb_symbmm(da,ptilde,dap,info)
if (info == psb_success_) call psb_numbmm(da,ptilde,dap)
call psb_spspmm(da,ptilde,dap,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
goto 9999
end if
@ -240,8 +234,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=am4)
if (info == psb_success_) call am4%free()
call psb_symbmm(da,atmp,dadap,info)
call psb_numbmm(da,atmp,dadap)
call psb_spspmm(da,atmp,dadap,info)
call atmp%free()
! !$ write(0,*) 'Columns of AP',psb_sp_get_ncols(ap)
@ -249,7 +242,6 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
call dap%mv_to(csc_dap)
call dadap%mv_to(csc_dadap)
call csc_mat_col_prod(csc_dap,csc_dadap,omp,info)
call csc_mat_col_prod(csc_dadap,csc_dadap,oden,info)
call psb_sum(ictxt,omp)
@ -287,7 +279,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
!
! Build the filtered matrix Af from A
!
call a%cscnv(acsrf,info,dupl=psb_dupl_add_)
call la%cscnv(acsrf,info,dupl=psb_dupl_add_)
do i=1,nrow
tmp = czero
@ -327,23 +319,14 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
call af%mv_from(acsrf)
!
! Symbmm90 does the allocation for its result.
!
! op_prol = (I-w*D*Af)Ptilde
! Doing it this way means to consider diag(Af_i)
!
!
call psb_symbmm(af,ptilde,op_prol,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999
end if
call psb_numbmm(af,ptilde,op_prol)
call psb_spspmm(af,ptilde,op_prol,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
& 'Done SPSPMM 1'
else
!
! Build the smoothed prolongator using the original matrix
@ -363,19 +346,11 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
& write(debug_unit,*) me,' ',trim(name),&
& 'Done gather, going for SYMBMM 1'
!
! Symbmm90 does the allocation for its result.
!
! op_prol = (I-w*D*A)Ptilde
!
!
call psb_symbmm(am3,ptilde,op_prol,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999
end if
call psb_numbmm(am3,ptilde,op_prol)
call psb_spspmm(am3,ptilde,op_prol,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
@ -387,11 +362,11 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
! Ok, let's start over with the restrictor
!
call ptilde%transc(rtilde)
call a%cscnv(atmp,info,type='csr')
call la%cscnv(atmp,info,type='csr')
call psb_sphalo(atmp,desc_a,am4,info,&
& colcnv=.true.,rowscale=.true.)
nrt = am4%get_nrows()
call am4%csclip(atmp2,info,ione,nrt,ione,ncol)
call am4%csclip(atmp2,info,lone,nrt,lone,ncol)
call atmp2%cscnv(info,type='CSR')
if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=atmp2)
call am4%free()
@ -400,13 +375,12 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
! This is to compute the transpose. It ONLY works if the
! original A has a symmetric pattern.
call atmp%transc(atmp2)
call atmp2%csclip(dat,info,ione,nrow,ione,ncol)
call atmp2%csclip(dat,info,lone,nrow,lone,ncol)
call dat%cscnv(info,type='csr')
call dat%scal(adinv,info)
! Now for the product.
call psb_symbmm(dat,ptilde,datp,info)
if (info == psb_success_) call psb_numbmm(dat,ptilde,datp)
call psb_spspmm(dat,ptilde,datp,info)
call datp%clone(atmp2,info)
call psb_sphalo(atmp2,desc_a,am4,info,&
@ -484,8 +458,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
call rtilde%mv_from(tmpcoo)
call rtilde%cscnv(info,type='csr')
call psb_symbmm(rtilde,atmp,op_restr,info)
call psb_numbmm(rtilde,atmp,op_restr)
call psb_spspmm(rtilde,atmp,op_restr,info)
!
! Now we have to gather the halo of op_prol, and add it to itself
@ -527,16 +500,10 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
& write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd'
call psb_symbmm(a,tmp_prol,am3,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='symbmm 2')
goto 9999
end if
call psb_numbmm(a,tmp_prol,am3)
call psb_spspmm(la,tmp_prol,am3,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2'
& 'Done SPSPMM 2'
call psb_sphalo(am3,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.)
@ -552,8 +519,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
& write(debug_unit,*) me,' ',trim(name),&
& 'Done sphalo/ rwxtd'
call psb_symbmm(op_restr,am3,ac,info)
if (info == psb_success_) call psb_numbmm(op_restr,am3,ac)
call psb_spspmm(op_restr,am3,ac,info)
if (info == psb_success_) call am3%free()
if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_)
@ -581,11 +547,11 @@ contains
subroutine csc_mat_col_prod(a,b,v,info)
implicit none
type(psb_c_csc_sparse_mat), intent(in) :: a, b
type(psb_lc_csc_sparse_mat), intent(in) :: a, b
complex(psb_spk_), intent(out) :: v(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i,j,k, nr, nc,iap,nra,ibp,nrb
integer(psb_lpk_) :: i,j,k, nr, nc,iap,nra,ibp,nrb
info = psb_success_
nc = a%get_ncols()
@ -609,11 +575,11 @@ contains
subroutine csr_mat_row_prod(a,b,v,info)
implicit none
type(psb_c_csr_sparse_mat), intent(in) :: a, b
type(psb_lc_csr_sparse_mat), intent(in) :: a, b
complex(psb_spk_), intent(out) :: v(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i,j,k, nr, nc,iap,nca,ibp,ncb
integer(psb_lpk_) :: i,j,k, nr, nc,iap,nca,ibp,ncb
info = psb_success_
nr = a%get_nrows()
@ -637,12 +603,12 @@ contains
function sparse_srtd_dot(nv1,iv1,v1,nv2,iv2,v2) result(dot)
implicit none
integer(psb_ipk_), intent(in) :: nv1,nv2
integer(psb_ipk_), intent(in) :: iv1(:), iv2(:)
integer(psb_lpk_), intent(in) :: nv1,nv2
integer(psb_lpk_), intent(in) :: iv1(:), iv2(:)
complex(psb_spk_), intent(in) :: v1(:),v2(:)
complex(psb_spk_) :: dot
integer(psb_ipk_) :: i,j,k, ip1, ip2
integer(psb_lpk_) :: i,j,k, ip1, ip2
dot = czero
ip1 = 1
@ -665,7 +631,7 @@ contains
end function sparse_srtd_dot
subroutine local_dump(me,mat,name,header)
type(psb_cspmat_type), intent(in) :: mat
type(psb_lcspmat_type), intent(in) :: mat
integer(psb_ipk_), intent(in) :: me
character(len=*), intent(in) :: name
character(len=*), intent(in) :: header

@ -106,28 +106,30 @@ subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
! Arguments
type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:)
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_sml_parms), intent(inout) :: parms
type(psb_cspmat_type), intent(inout) :: op_prol
type(psb_cspmat_type), intent(out) :: ac,op_restr
type(psb_lcspmat_type), intent(inout) :: op_prol
type(psb_lcspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo
integer(psb_ipk_) :: ictxt, np, me, icomm, minfo
character(len=20) :: name
integer(psb_ipk_) :: ierr(5)
type(psb_c_coo_sparse_mat) :: ac_coo, tmpcoo
type(psb_c_csr_sparse_mat) :: acsr1, acsr2
integer(psb_ipk_) :: ierr(5)
type(psb_lcspmat_type) :: la
type(psb_lc_coo_sparse_mat) :: ac_coo, tmpcoo
type(psb_lc_csr_sparse_mat) :: acsr1, acsr2
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nrow, nglob, ncol, ntaggr, nzl, ip, &
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nzl, ip, &
& naggr, nzt, naggrm1, naggrp1, i, k
name='mld_aggrmat_nosmth_asb'
if(psb_get_errstatus().ne.0) return
info=psb_success_
name = 'mld_aggrmat_nosmth_asb'
info = psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
ictxt = desc_a%get_context()
icomm = desc_a%get_mpic()
@ -136,24 +138,23 @@ subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
nrow = desc_a%get_local_rows()
ncol = desc_a%get_local_cols()
naggr = nlaggr(me+1)
ntaggr = sum(nlaggr)
naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1))
call op_prol%cp_to(tmpcoo)
call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_)
call tmpcoo%transp()
nzl = tmpcoo%get_nzeros()
i=0
!
! Now we have to fix this. The only rows of tmpcoo/op_restr that are correct
! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
!
do k=1, nzl
nzl = tmpcoo%get_nzeros()
i = 0
do k = 1, nzl
if ((naggrm1 < tmpcoo%ia(k)) .and.(tmpcoo%ia(k) <= naggrp1)) then
i = i+1
tmpcoo%val(i) = tmpcoo%val(k)
@ -168,11 +169,11 @@ subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
if (info /= psb_success_) goto 9999
call a%cp_to(ac_coo)
call a%cp_to_l(la)
call la%mv_to(ac_coo)
nzt = ac_coo%get_nzeros()
k = 0
do i=1, nzt
k = 0
do i = 1, nzt
k = k + 1
ac_coo%ia(k) = ilaggr(ac_coo%ia(i))
ac_coo%ja(k) = ilaggr(ac_coo%ja(i))

@ -65,14 +65,7 @@
! specified by the user through mld_cprecinit and mld_zprecset.
! On output from this routine the entries of AC, op_prol, op_restr
! are still in "global numbering" mode; this is fixed in the calling routine
! mld_c_lev_aggrmat_asb.
!
! For more details see
! M. Brezina and P. Vanek, A black-box iterative solver based on a
! two-level Schwarz method, Computing, 63 (1999), 233-263.
! P. D'Ambra, D. di Serafino and S. Filippone, On the development of
! PSBLAS-based parallel two-level Schwarz preconditioners, Appl. Num. Math.
! 57 (2007), 1181-1196.
! aggregator%mat_asb.
!
!
! Arguments:
@ -117,37 +110,38 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
implicit none
! Arguments
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(inout) :: op_prol
type(psb_cspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_sml_parms), intent(inout) :: parms
type(psb_lcspmat_type), intent(inout) :: op_prol
type(psb_lcspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_) :: nrow, nglob, ncol, ntaggr, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw, err_act
integer(psb_ipk_) ::ictxt, np, me
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, ip, &
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw
integer(psb_ipk_) :: ictxt, np, me
character(len=20) :: name
type(psb_cspmat_type) :: am3, am4, tmp_prol
type(psb_c_coo_sparse_mat) :: tmpcoo
type(psb_c_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
type(psb_lcspmat_type) :: la, am3, am4, tmp_prol
type(psb_lc_coo_sparse_mat) :: tmpcoo
type(psb_lc_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
complex(psb_spk_), allocatable :: adiag(:)
integer(psb_ipk_) :: ierr(5)
logical :: filter_mat
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: debug_level, debug_unit, err_act
integer(psb_ipk_), parameter :: ncmax=16
real(psb_spk_) :: anorm, omega, tmp, dg, theta
name='mld_aggrmat_smth_asb'
if(psb_get_errstatus().ne.0) return
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context()
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
@ -176,6 +170,7 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
& call psb_realloc(ncol,adiag,info)
if (info == psb_success_) &
& call psb_halo(adiag,desc_a,info)
if (info == psb_success_) call a%cp_to_l(la)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
@ -185,7 +180,7 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
! 1. Allocate Ptilde in sparse matrix form
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_) call la%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (info /= psb_success_) goto 9999
if (debug_level >= psb_debug_outer_) &
@ -198,7 +193,7 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
!
if (info == psb_success_) call acsr3%cp_to_fmt(acsrf,info)
do i=1,nrow
do i=1, nrow
tmp = czero
jd = -1
do j=acsrf%irp(i),acsrf%irp(i+1)-1
@ -342,7 +337,7 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
goto 9999
end if
call psb_spspmm(a,tmp_prol,am3,info)
call psb_spspmm(la,tmp_prol,am3,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 2')
goto 9999

@ -143,23 +143,25 @@ subroutine mld_d_dec_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,ac,op_p
type(mld_dml_parms), intent(inout) :: parms
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(psb_dspmat_type), intent(inout) :: op_prol
type(psb_dspmat_type), intent(out) :: ac,op_restr
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(psb_ldspmat_type), intent(inout) :: op_prol
type(psb_ldspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
character(len=20) :: name
integer(psb_mpik_) :: ictxt, np, me
type(psb_d_coo_sparse_mat) :: acoo, bcoo
type(psb_d_csr_sparse_mat) :: acsr1
integer(psb_ipk_) :: nzl,ntaggr
integer(psb_ipk_) :: err_act
character(len=20) :: name
integer(psb_mpk_) :: ictxt, np, me
type(psb_ld_coo_sparse_mat) :: acoo, bcoo
type(psb_ld_csr_sparse_mat) :: acsr1
integer(psb_lpk_) :: nzl,ntaggr
integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: debug_level, debug_unit
name='mld_d_dec_aggregator_mat_asb'
if (psb_get_errstatus().ne.0) return
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
info = psb_success_

@ -81,20 +81,22 @@ subroutine mld_d_dec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_
type(mld_dml_parms), intent(inout) :: parms
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_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
type(psb_ldspmat_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_) :: ntaggr
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name
integer(psb_mpk_) :: ictxt, np, me
integer(psb_ipk_) :: err_act
integer(psb_lpk_) :: ntaggr
integer(psb_ipk_) :: debug_level, debug_unit
name='mld_d_dec_aggregator_tprol'
if (psb_get_errstatus().ne.0) return
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
info = psb_success_

@ -87,22 +87,24 @@ subroutine mld_d_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info)
! Arguments
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), allocatable, intent(inout) :: ilaggr(:),nlaggr(:)
type(psb_dspmat_type), intent(out) :: op_prol
integer(psb_lpk_), allocatable, intent(inout) :: ilaggr(:),nlaggr(:)
type(psb_ldspmat_type), intent(out) :: op_prol
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m,naggrm1, naggrp1, ntaggr
type(psb_d_coo_sparse_mat) :: tmpcoo
integer(psb_lpk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m,naggrm1, naggrp1, ntaggr
type(psb_ld_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
integer(psb_lpk_) :: nrow, ncol, n_ne
character(len=20) :: name, ch_err
if(psb_get_errstatus() /= 0) return
info=psb_success_
name = 'mld_map_to_tprol'
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!

@ -80,12 +80,13 @@ subroutine mld_d_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
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_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),&
& ideg(:), idxs(:), tmpaggr(:)
integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),&
& ideg(:), idxs(:)
integer(psb_lpk_), allocatable :: tmpaggr(:)
real(psb_dpk_), allocatable :: val(:), diag(:)
integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, nc, naggr,i,j,m, nz, ilg, ii, ip
type(psb_d_csr_sparse_mat) :: acsr
@ -96,10 +97,12 @@ subroutine mld_d_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
integer(psb_ipk_) :: nrow, ncol, n_ne
character(len=20) :: name, ch_err
if (psb_get_errstatus() /= 0) return
info=psb_success_
name = 'mld_soc1_map_bld'
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!

@ -79,12 +79,13 @@ subroutine mld_d_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
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_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),&
& ideg(:), idxs(:), tmpaggr(:)
& ideg(:), idxs(:)
integer(psb_lpk_), allocatable :: tmpaggr(:)
real(psb_dpk_), allocatable :: val(:), diag(:)
integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr,nc,naggr,i,j,m, nz, ilg, ii, ip, ip1,nzcnt
type(psb_d_csr_sparse_mat) :: acsr, muij, s_neigh
@ -96,10 +97,12 @@ subroutine mld_d_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
integer(psb_ipk_) :: nrow, ncol, n_ne
character(len=20) :: name, ch_err
if (psb_get_errstatus() /= 0) return
info=psb_success_
name = 'mld_soc2_map_bld'
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!

@ -85,21 +85,24 @@ subroutine mld_d_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,
type(mld_dml_parms), intent(inout) :: parms
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_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
type(psb_ldspmat_type), intent(out) :: op_prol
integer(psb_ipk_), intent(out) :: info
! Local variables
type(psb_dspmat_type) :: atmp, atrans
character(len=20) :: name
integer(psb_mpik_) :: ictxt, np, me
integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: ntaggr, nr
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name
integer(psb_mpk_) :: ictxt, np, me
integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: nr
integer(psb_lpk_) :: ntaggr
integer(psb_ipk_) :: debug_level, debug_unit
name='mld_d_symdec_aggregator_tprol'
if (psb_get_errstatus().ne.0) return
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
info = psb_success_

@ -88,33 +88,35 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
implicit none
! Arguments
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(psb_dspmat_type), intent(inout) :: op_prol
type(psb_dspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_dml_parms), intent(inout) :: parms
type(psb_ldspmat_type), intent(inout) :: op_prol
type(psb_ldspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_) :: nrow, nglob, ncol, ntaggr, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw, err_act
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, ip, &
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw
integer(psb_ipk_) ::ictxt, np, me
character(len=20) :: name
type(psb_dspmat_type) :: am3, am4,tmp_prol
type(psb_d_coo_sparse_mat) :: tmpcoo
type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
type(psb_ldspmat_type) :: am3, am4,tmp_prol, la
type(psb_ld_coo_sparse_mat) :: tmpcoo
type(psb_ld_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
real(psb_dpk_), allocatable :: adiag(:)
integer(psb_ipk_) :: ierr(5)
logical :: filter_mat
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: debug_level, debug_unit, err_act
integer(psb_ipk_), parameter :: ncmax=16
real(psb_dpk_) :: anorm, omega, tmp, dg, theta
name='mld_aggrmat_biz_asb'
if(psb_get_errstatus().ne.0) return
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
@ -128,10 +130,8 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
ncol = desc_a%get_local_cols()
theta = parms%aggr_thresh
naggr = nlaggr(me+1)
ntaggr = sum(nlaggr)
filter_mat = (parms%aggr_filter == mld_filter_mat_)
! naggr: number of local aggregates
@ -143,8 +143,8 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
& call psb_realloc(ncol,adiag,info)
if (info == psb_success_) &
& call psb_halo(adiag,desc_a,info)
if(info /= psb_success_) then
if (info == psb_success_) call a%cp_to_l(la)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
goto 9999
end if
@ -152,7 +152,7 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
! 1. Allocate Ptilde in sparse matrix form
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_) call la%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (info /= psb_success_) goto 9999
if (debug_level >= psb_debug_outer_) &
@ -163,9 +163,9 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
!
! Build the filtered matrix Af from A
!
if (info == psb_success_) call a%cscnv(acsrf,info,dupl=psb_dupl_add_)
if (info == psb_success_) call acsr3%cp_to_fmt(acsrf,info)
do i=1,nrow
do i=1, nrow
tmp = dzero
jd = -1
do j=acsrf%irp(i),acsrf%irp(i+1)-1
@ -271,17 +271,15 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
! Doing it this way means to consider diag(Af_i)
!
!
call psb_symbmm(acsrf,ptilde,acsr1,info)
call psb_spspmm(acsrf,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
goto 9999
end if
call psb_numbmm(acsrf,ptilde,acsr1)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
& 'Done SPSPMM 1'
else
!
@ -307,17 +305,15 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
! Doing it this way means to consider diag(A_i)
!
!
call psb_symbmm(acsr3,ptilde,acsr1,info)
call psb_spspmm(acsr3,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
goto 9999
end if
call psb_numbmm(acsr3,ptilde,acsr1)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
& 'Done SPSPMM 1'
end if
call ptilde%free()
@ -331,16 +327,15 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
goto 9999
end if
call psb_symbmm(a,tmp_prol,am3,info)
call psb_spspmm(la,tmp_prol,am3,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2')
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 2')
goto 9999
end if
call psb_numbmm(a,tmp_prol,am3)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2',parms%aggr_prol, mld_smooth_prol_
& 'Done SPSPMM 2',parms%aggr_prol, mld_smooth_prol_
call tmp_prol%transp(op_restr)
if (debug_level >= psb_debug_outer_) &
@ -353,23 +348,18 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting symbmm 3'
call psb_symbmm(op_restr,am3,ac,info)
if (info == psb_success_) call psb_numbmm(op_restr,am3,ac)
& 'starting spspmm 3'
call psb_spspmm(op_restr,am3,ac,info)
if (info == psb_success_) call am3%free()
if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_)
if (info == psb_success_) call ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Build b = op_restr x am3')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done smooth_aggregate '

@ -114,37 +114,38 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
! Arguments
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(psb_dspmat_type), intent(inout) :: op_prol
type(psb_dspmat_type), intent(out) :: ac,op_restr
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_dml_parms), intent(inout) :: parms
type(psb_ldspmat_type), intent(inout) :: op_prol
type(psb_ldspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_), allocatable :: nzbr(:), idisp(:)
integer(psb_ipk_) :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrt, err_act
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrt
integer(psb_ipk_) :: ictxt,np,me, icomm
character(len=20) :: name
type(psb_dspmat_type) :: af, ptilde, rtilde, atran, atp, atdatp
type(psb_dspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da
type(psb_dspmat_type) :: dat, datp, datdatp, atmp3, tmp_prol
type(psb_d_coo_sparse_mat) :: tmpcoo
type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf
type(psb_d_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc
type(psb_ldspmat_type) :: la, af, ptilde, rtilde, atran, atp, atdatp
type(psb_ldspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da
type(psb_ldspmat_type) :: dat, datp, datdatp, atmp3, tmp_prol
type(psb_ld_coo_sparse_mat) :: tmpcoo
type(psb_ld_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf
type(psb_ld_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc
real(psb_dpk_), allocatable :: adiag(:), adinv(:)
real(psb_dpk_), allocatable :: omf(:), omp(:), omi(:), oden(:)
logical :: filter_mat
integer(psb_ipk_) :: ierr(5)
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: debug_level, debug_unit, err_act
integer(psb_ipk_), parameter :: ncmax=16
real(psb_dpk_) :: anorm, theta
real(psb_dpk_) :: tmp, alpha, beta, ommx
name='mld_aggrmat_minnrg'
if(psb_get_errstatus().ne.0) return
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
@ -163,13 +164,6 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
naggr = nlaggr(me+1)
ntaggr = sum(nlaggr)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_; ierr(1)=2*np;
call psb_errpush(info,name,i_err=ierr,a_err='integer')
goto 9999
end if
naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1))
@ -193,6 +187,11 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
& call psb_realloc(ncol,adiag,info)
if (info == psb_success_) &
& call psb_halo(adiag,desc_a,info)
if (info == psb_success_) call a%cp_to_l(la)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
goto 9999
end if
do i=1,size(adiag)
if (adiag(i) /= dzero) then
@ -202,10 +201,6 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
end if
end do
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
goto 9999
end if
! 1. Allocate Ptilde in sparse matrix form
@ -213,8 +208,8 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
call ptilde%mv_from(tmpcoo)
call ptilde%cscnv(info,type='csr')
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_) call la%cscnv(am3,info,type='csr',dupl=psb_dupl_add_)
if (info == psb_success_) call la%cscnv(da,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
@ -225,11 +220,10 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
call da%scal(adinv,info)
call psb_symbmm(da,ptilde,dap,info)
if (info == psb_success_) call psb_numbmm(da,ptilde,dap)
call psb_spspmm(da,ptilde,dap,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
goto 9999
end if
@ -240,8 +234,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=am4)
if (info == psb_success_) call am4%free()
call psb_symbmm(da,atmp,dadap,info)
call psb_numbmm(da,atmp,dadap)
call psb_spspmm(da,atmp,dadap,info)
call atmp%free()
! !$ write(0,*) 'Columns of AP',psb_sp_get_ncols(ap)
@ -249,7 +242,6 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
call dap%mv_to(csc_dap)
call dadap%mv_to(csc_dadap)
call csc_mat_col_prod(csc_dap,csc_dadap,omp,info)
call csc_mat_col_prod(csc_dadap,csc_dadap,oden,info)
call psb_sum(ictxt,omp)
@ -287,7 +279,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
!
! Build the filtered matrix Af from A
!
call a%cscnv(acsrf,info,dupl=psb_dupl_add_)
call la%cscnv(acsrf,info,dupl=psb_dupl_add_)
do i=1,nrow
tmp = dzero
@ -327,23 +319,14 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
call af%mv_from(acsrf)
!
! Symbmm90 does the allocation for its result.
!
! op_prol = (I-w*D*Af)Ptilde
! Doing it this way means to consider diag(Af_i)
!
!
call psb_symbmm(af,ptilde,op_prol,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999
end if
call psb_numbmm(af,ptilde,op_prol)
call psb_spspmm(af,ptilde,op_prol,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
& 'Done SPSPMM 1'
else
!
! Build the smoothed prolongator using the original matrix
@ -363,19 +346,11 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
& write(debug_unit,*) me,' ',trim(name),&
& 'Done gather, going for SYMBMM 1'
!
! Symbmm90 does the allocation for its result.
!
! op_prol = (I-w*D*A)Ptilde
!
!
call psb_symbmm(am3,ptilde,op_prol,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999
end if
call psb_numbmm(am3,ptilde,op_prol)
call psb_spspmm(am3,ptilde,op_prol,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
@ -387,11 +362,11 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
! Ok, let's start over with the restrictor
!
call ptilde%transc(rtilde)
call a%cscnv(atmp,info,type='csr')
call la%cscnv(atmp,info,type='csr')
call psb_sphalo(atmp,desc_a,am4,info,&
& colcnv=.true.,rowscale=.true.)
nrt = am4%get_nrows()
call am4%csclip(atmp2,info,ione,nrt,ione,ncol)
call am4%csclip(atmp2,info,lone,nrt,lone,ncol)
call atmp2%cscnv(info,type='CSR')
if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=atmp2)
call am4%free()
@ -400,13 +375,12 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
! This is to compute the transpose. It ONLY works if the
! original A has a symmetric pattern.
call atmp%transc(atmp2)
call atmp2%csclip(dat,info,ione,nrow,ione,ncol)
call atmp2%csclip(dat,info,lone,nrow,lone,ncol)
call dat%cscnv(info,type='csr')
call dat%scal(adinv,info)
! Now for the product.
call psb_symbmm(dat,ptilde,datp,info)
if (info == psb_success_) call psb_numbmm(dat,ptilde,datp)
call psb_spspmm(dat,ptilde,datp,info)
call datp%clone(atmp2,info)
call psb_sphalo(atmp2,desc_a,am4,info,&
@ -484,8 +458,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
call rtilde%mv_from(tmpcoo)
call rtilde%cscnv(info,type='csr')
call psb_symbmm(rtilde,atmp,op_restr,info)
call psb_numbmm(rtilde,atmp,op_restr)
call psb_spspmm(rtilde,atmp,op_restr,info)
!
! Now we have to gather the halo of op_prol, and add it to itself
@ -527,16 +500,10 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
& write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd'
call psb_symbmm(a,tmp_prol,am3,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='symbmm 2')
goto 9999
end if
call psb_numbmm(a,tmp_prol,am3)
call psb_spspmm(la,tmp_prol,am3,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2'
& 'Done SPSPMM 2'
call psb_sphalo(am3,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.)
@ -552,8 +519,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
& write(debug_unit,*) me,' ',trim(name),&
& 'Done sphalo/ rwxtd'
call psb_symbmm(op_restr,am3,ac,info)
if (info == psb_success_) call psb_numbmm(op_restr,am3,ac)
call psb_spspmm(op_restr,am3,ac,info)
if (info == psb_success_) call am3%free()
if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_)
@ -581,11 +547,11 @@ contains
subroutine csc_mat_col_prod(a,b,v,info)
implicit none
type(psb_d_csc_sparse_mat), intent(in) :: a, b
type(psb_ld_csc_sparse_mat), intent(in) :: a, b
real(psb_dpk_), intent(out) :: v(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i,j,k, nr, nc,iap,nra,ibp,nrb
integer(psb_lpk_) :: i,j,k, nr, nc,iap,nra,ibp,nrb
info = psb_success_
nc = a%get_ncols()
@ -609,11 +575,11 @@ contains
subroutine csr_mat_row_prod(a,b,v,info)
implicit none
type(psb_d_csr_sparse_mat), intent(in) :: a, b
type(psb_ld_csr_sparse_mat), intent(in) :: a, b
real(psb_dpk_), intent(out) :: v(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i,j,k, nr, nc,iap,nca,ibp,ncb
integer(psb_lpk_) :: i,j,k, nr, nc,iap,nca,ibp,ncb
info = psb_success_
nr = a%get_nrows()
@ -637,12 +603,12 @@ contains
function sparse_srtd_dot(nv1,iv1,v1,nv2,iv2,v2) result(dot)
implicit none
integer(psb_ipk_), intent(in) :: nv1,nv2
integer(psb_ipk_), intent(in) :: iv1(:), iv2(:)
integer(psb_lpk_), intent(in) :: nv1,nv2
integer(psb_lpk_), intent(in) :: iv1(:), iv2(:)
real(psb_dpk_), intent(in) :: v1(:),v2(:)
real(psb_dpk_) :: dot
integer(psb_ipk_) :: i,j,k, ip1, ip2
integer(psb_lpk_) :: i,j,k, ip1, ip2
dot = dzero
ip1 = 1
@ -665,7 +631,7 @@ contains
end function sparse_srtd_dot
subroutine local_dump(me,mat,name,header)
type(psb_dspmat_type), intent(in) :: mat
type(psb_ldspmat_type), intent(in) :: mat
integer(psb_ipk_), intent(in) :: me
character(len=*), intent(in) :: name
character(len=*), intent(in) :: header

@ -106,28 +106,30 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
! Arguments
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:)
integer(psb_lpk_), 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_ldspmat_type), intent(inout) :: op_prol
type(psb_ldspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo
integer(psb_ipk_) :: ictxt, np, me, icomm, minfo
character(len=20) :: name
integer(psb_ipk_) :: ierr(5)
type(psb_d_coo_sparse_mat) :: ac_coo, tmpcoo
type(psb_d_csr_sparse_mat) :: acsr1, acsr2
integer(psb_ipk_) :: ierr(5)
type(psb_ldspmat_type) :: la
type(psb_ld_coo_sparse_mat) :: ac_coo, tmpcoo
type(psb_ld_csr_sparse_mat) :: acsr1, acsr2
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nrow, nglob, ncol, ntaggr, nzl, ip, &
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nzl, ip, &
& naggr, nzt, naggrm1, naggrp1, i, k
name='mld_aggrmat_nosmth_asb'
if(psb_get_errstatus().ne.0) return
info=psb_success_
name = 'mld_aggrmat_nosmth_asb'
info = psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
ictxt = desc_a%get_context()
icomm = desc_a%get_mpic()
@ -136,24 +138,23 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
nrow = desc_a%get_local_rows()
ncol = desc_a%get_local_cols()
naggr = nlaggr(me+1)
ntaggr = sum(nlaggr)
naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1))
call op_prol%cp_to(tmpcoo)
call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_)
call tmpcoo%transp()
nzl = tmpcoo%get_nzeros()
i=0
!
! Now we have to fix this. The only rows of tmpcoo/op_restr that are correct
! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
!
do k=1, nzl
nzl = tmpcoo%get_nzeros()
i = 0
do k = 1, nzl
if ((naggrm1 < tmpcoo%ia(k)) .and.(tmpcoo%ia(k) <= naggrp1)) then
i = i+1
tmpcoo%val(i) = tmpcoo%val(k)
@ -168,11 +169,11 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
if (info /= psb_success_) goto 9999
call a%cp_to(ac_coo)
call a%cp_to_l(la)
call la%mv_to(ac_coo)
nzt = ac_coo%get_nzeros()
k = 0
do i=1, nzt
k = 0
do i = 1, nzt
k = k + 1
ac_coo%ia(k) = ilaggr(ac_coo%ia(i))
ac_coo%ja(k) = ilaggr(ac_coo%ja(i))

@ -65,14 +65,7 @@
! specified by the user through mld_dprecinit and mld_zprecset.
! On output from this routine the entries of AC, op_prol, op_restr
! are still in "global numbering" mode; this is fixed in the calling routine
! mld_d_lev_aggrmat_asb.
!
! For more details see
! M. Brezina and P. Vanek, A black-box iterative solver based on a
! two-level Schwarz method, Computing, 63 (1999), 233-263.
! P. D'Ambra, D. di Serafino and S. Filippone, On the development of
! PSBLAS-based parallel two-level Schwarz preconditioners, Appl. Num. Math.
! 57 (2007), 1181-1196.
! aggregator%mat_asb.
!
!
! Arguments:
@ -117,37 +110,38 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
implicit none
! Arguments
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(psb_dspmat_type), intent(inout) :: op_prol
type(psb_dspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_dml_parms), intent(inout) :: parms
type(psb_ldspmat_type), intent(inout) :: op_prol
type(psb_ldspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_) :: nrow, nglob, ncol, ntaggr, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw, err_act
integer(psb_ipk_) ::ictxt, np, me
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, ip, &
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw
integer(psb_ipk_) :: ictxt, np, me
character(len=20) :: name
type(psb_dspmat_type) :: am3, am4, tmp_prol
type(psb_d_coo_sparse_mat) :: tmpcoo
type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
type(psb_ldspmat_type) :: la, am3, am4, tmp_prol
type(psb_ld_coo_sparse_mat) :: tmpcoo
type(psb_ld_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
real(psb_dpk_), allocatable :: adiag(:)
integer(psb_ipk_) :: ierr(5)
logical :: filter_mat
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: debug_level, debug_unit, err_act
integer(psb_ipk_), parameter :: ncmax=16
real(psb_dpk_) :: anorm, omega, tmp, dg, theta
name='mld_aggrmat_smth_asb'
if(psb_get_errstatus().ne.0) return
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context()
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
@ -176,6 +170,7 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
& call psb_realloc(ncol,adiag,info)
if (info == psb_success_) &
& call psb_halo(adiag,desc_a,info)
if (info == psb_success_) call a%cp_to_l(la)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
@ -185,7 +180,7 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
! 1. Allocate Ptilde in sparse matrix form
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_) call la%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (info /= psb_success_) goto 9999
if (debug_level >= psb_debug_outer_) &
@ -198,7 +193,7 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
!
if (info == psb_success_) call acsr3%cp_to_fmt(acsrf,info)
do i=1,nrow
do i=1, nrow
tmp = dzero
jd = -1
do j=acsrf%irp(i),acsrf%irp(i+1)-1
@ -342,7 +337,7 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
goto 9999
end if
call psb_spspmm(a,tmp_prol,am3,info)
call psb_spspmm(la,tmp_prol,am3,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 2')
goto 9999

@ -143,23 +143,25 @@ subroutine mld_s_dec_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,ac,op_p
type(mld_sml_parms), intent(inout) :: parms
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
type(psb_sspmat_type), intent(out) :: ac,op_restr
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(psb_lsspmat_type), intent(inout) :: op_prol
type(psb_lsspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
character(len=20) :: name
integer(psb_mpik_) :: ictxt, np, me
type(psb_s_coo_sparse_mat) :: acoo, bcoo
type(psb_s_csr_sparse_mat) :: acsr1
integer(psb_ipk_) :: nzl,ntaggr
integer(psb_ipk_) :: err_act
character(len=20) :: name
integer(psb_mpk_) :: ictxt, np, me
type(psb_ls_coo_sparse_mat) :: acoo, bcoo
type(psb_ls_csr_sparse_mat) :: acsr1
integer(psb_lpk_) :: nzl,ntaggr
integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: debug_level, debug_unit
name='mld_s_dec_aggregator_mat_asb'
if (psb_get_errstatus().ne.0) return
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
info = psb_success_

@ -81,20 +81,22 @@ subroutine mld_s_dec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_
type(mld_sml_parms), intent(inout) :: parms
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_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
type(psb_lsspmat_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_) :: ntaggr
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name
integer(psb_mpk_) :: ictxt, np, me
integer(psb_ipk_) :: err_act
integer(psb_lpk_) :: ntaggr
integer(psb_ipk_) :: debug_level, debug_unit
name='mld_s_dec_aggregator_tprol'
if (psb_get_errstatus().ne.0) return
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
info = psb_success_

@ -87,22 +87,24 @@ subroutine mld_s_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info)
! Arguments
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), allocatable, intent(inout) :: ilaggr(:),nlaggr(:)
type(psb_sspmat_type), intent(out) :: op_prol
integer(psb_lpk_), allocatable, intent(inout) :: ilaggr(:),nlaggr(:)
type(psb_lsspmat_type), intent(out) :: op_prol
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m,naggrm1, naggrp1, ntaggr
type(psb_s_coo_sparse_mat) :: tmpcoo
integer(psb_lpk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m,naggrm1, naggrp1, ntaggr
type(psb_ls_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
integer(psb_lpk_) :: nrow, ncol, n_ne
character(len=20) :: name, ch_err
if(psb_get_errstatus() /= 0) return
info=psb_success_
name = 'mld_map_to_tprol'
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!

@ -80,12 +80,13 @@ subroutine mld_s_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
real(psb_spk_), intent(in) :: theta
integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),&
& ideg(:), idxs(:), tmpaggr(:)
integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),&
& ideg(:), idxs(:)
integer(psb_lpk_), allocatable :: tmpaggr(:)
real(psb_spk_), allocatable :: val(:), diag(:)
integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, nc, naggr,i,j,m, nz, ilg, ii, ip
type(psb_s_csr_sparse_mat) :: acsr
@ -96,10 +97,12 @@ subroutine mld_s_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
integer(psb_ipk_) :: nrow, ncol, n_ne
character(len=20) :: name, ch_err
if (psb_get_errstatus() /= 0) return
info=psb_success_
name = 'mld_soc1_map_bld'
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!

@ -79,12 +79,13 @@ subroutine mld_s_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
real(psb_spk_), intent(in) :: theta
integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),&
& ideg(:), idxs(:), tmpaggr(:)
& ideg(:), idxs(:)
integer(psb_lpk_), allocatable :: tmpaggr(:)
real(psb_spk_), allocatable :: val(:), diag(:)
integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr,nc,naggr,i,j,m, nz, ilg, ii, ip, ip1,nzcnt
type(psb_s_csr_sparse_mat) :: acsr, muij, s_neigh
@ -96,10 +97,12 @@ subroutine mld_s_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
integer(psb_ipk_) :: nrow, ncol, n_ne
character(len=20) :: name, ch_err
if (psb_get_errstatus() /= 0) return
info=psb_success_
name = 'mld_soc2_map_bld'
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!

@ -85,21 +85,24 @@ subroutine mld_s_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,
type(mld_sml_parms), intent(inout) :: parms
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_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
type(psb_lsspmat_type), intent(out) :: op_prol
integer(psb_ipk_), intent(out) :: info
! Local variables
type(psb_sspmat_type) :: atmp, atrans
character(len=20) :: name
integer(psb_mpik_) :: ictxt, np, me
integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: ntaggr, nr
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name
integer(psb_mpk_) :: ictxt, np, me
integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: nr
integer(psb_lpk_) :: ntaggr
integer(psb_ipk_) :: debug_level, debug_unit
name='mld_s_symdec_aggregator_tprol'
if (psb_get_errstatus().ne.0) return
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
info = psb_success_

@ -88,33 +88,35 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
implicit none
! Arguments
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(inout) :: op_prol
type(psb_sspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_sml_parms), intent(inout) :: parms
type(psb_lsspmat_type), intent(inout) :: op_prol
type(psb_lsspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_) :: nrow, nglob, ncol, ntaggr, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw, err_act
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, ip, &
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw
integer(psb_ipk_) ::ictxt, np, me
character(len=20) :: name
type(psb_sspmat_type) :: am3, am4,tmp_prol
type(psb_s_coo_sparse_mat) :: tmpcoo
type(psb_s_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
type(psb_lsspmat_type) :: am3, am4,tmp_prol, la
type(psb_ls_coo_sparse_mat) :: tmpcoo
type(psb_ls_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
real(psb_spk_), allocatable :: adiag(:)
integer(psb_ipk_) :: ierr(5)
logical :: filter_mat
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: debug_level, debug_unit, err_act
integer(psb_ipk_), parameter :: ncmax=16
real(psb_spk_) :: anorm, omega, tmp, dg, theta
name='mld_aggrmat_biz_asb'
if(psb_get_errstatus().ne.0) return
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
@ -128,10 +130,8 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
ncol = desc_a%get_local_cols()
theta = parms%aggr_thresh
naggr = nlaggr(me+1)
ntaggr = sum(nlaggr)
filter_mat = (parms%aggr_filter == mld_filter_mat_)
! naggr: number of local aggregates
@ -143,8 +143,8 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
& call psb_realloc(ncol,adiag,info)
if (info == psb_success_) &
& call psb_halo(adiag,desc_a,info)
if(info /= psb_success_) then
if (info == psb_success_) call a%cp_to_l(la)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
goto 9999
end if
@ -152,7 +152,7 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
! 1. Allocate Ptilde in sparse matrix form
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_) call la%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (info /= psb_success_) goto 9999
if (debug_level >= psb_debug_outer_) &
@ -163,9 +163,9 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
!
! Build the filtered matrix Af from A
!
if (info == psb_success_) call a%cscnv(acsrf,info,dupl=psb_dupl_add_)
if (info == psb_success_) call acsr3%cp_to_fmt(acsrf,info)
do i=1,nrow
do i=1, nrow
tmp = szero
jd = -1
do j=acsrf%irp(i),acsrf%irp(i+1)-1
@ -271,17 +271,15 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
! Doing it this way means to consider diag(Af_i)
!
!
call psb_symbmm(acsrf,ptilde,acsr1,info)
call psb_spspmm(acsrf,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
goto 9999
end if
call psb_numbmm(acsrf,ptilde,acsr1)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
& 'Done SPSPMM 1'
else
!
@ -307,17 +305,15 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
! Doing it this way means to consider diag(A_i)
!
!
call psb_symbmm(acsr3,ptilde,acsr1,info)
call psb_spspmm(acsr3,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
goto 9999
end if
call psb_numbmm(acsr3,ptilde,acsr1)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
& 'Done SPSPMM 1'
end if
call ptilde%free()
@ -331,16 +327,15 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
goto 9999
end if
call psb_symbmm(a,tmp_prol,am3,info)
call psb_spspmm(la,tmp_prol,am3,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2')
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 2')
goto 9999
end if
call psb_numbmm(a,tmp_prol,am3)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2',parms%aggr_prol, mld_smooth_prol_
& 'Done SPSPMM 2',parms%aggr_prol, mld_smooth_prol_
call tmp_prol%transp(op_restr)
if (debug_level >= psb_debug_outer_) &
@ -353,23 +348,18 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting symbmm 3'
call psb_symbmm(op_restr,am3,ac,info)
if (info == psb_success_) call psb_numbmm(op_restr,am3,ac)
& 'starting spspmm 3'
call psb_spspmm(op_restr,am3,ac,info)
if (info == psb_success_) call am3%free()
if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_)
if (info == psb_success_) call ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Build b = op_restr x am3')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done smooth_aggregate '

@ -114,37 +114,38 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
! Arguments
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(inout) :: op_prol
type(psb_sspmat_type), intent(out) :: ac,op_restr
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_sml_parms), intent(inout) :: parms
type(psb_lsspmat_type), intent(inout) :: op_prol
type(psb_lsspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_), allocatable :: nzbr(:), idisp(:)
integer(psb_ipk_) :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrt, err_act
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrt
integer(psb_ipk_) :: ictxt,np,me, icomm
character(len=20) :: name
type(psb_sspmat_type) :: af, ptilde, rtilde, atran, atp, atdatp
type(psb_sspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da
type(psb_sspmat_type) :: dat, datp, datdatp, atmp3, tmp_prol
type(psb_s_coo_sparse_mat) :: tmpcoo
type(psb_s_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf
type(psb_s_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc
type(psb_lsspmat_type) :: la, af, ptilde, rtilde, atran, atp, atdatp
type(psb_lsspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da
type(psb_lsspmat_type) :: dat, datp, datdatp, atmp3, tmp_prol
type(psb_ls_coo_sparse_mat) :: tmpcoo
type(psb_ls_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf
type(psb_ls_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc
real(psb_spk_), allocatable :: adiag(:), adinv(:)
real(psb_spk_), allocatable :: omf(:), omp(:), omi(:), oden(:)
logical :: filter_mat
integer(psb_ipk_) :: ierr(5)
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: debug_level, debug_unit, err_act
integer(psb_ipk_), parameter :: ncmax=16
real(psb_spk_) :: anorm, theta
real(psb_spk_) :: tmp, alpha, beta, ommx
name='mld_aggrmat_minnrg'
if(psb_get_errstatus().ne.0) return
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
@ -163,13 +164,6 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
naggr = nlaggr(me+1)
ntaggr = sum(nlaggr)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_; ierr(1)=2*np;
call psb_errpush(info,name,i_err=ierr,a_err='integer')
goto 9999
end if
naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1))
@ -193,6 +187,11 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
& call psb_realloc(ncol,adiag,info)
if (info == psb_success_) &
& call psb_halo(adiag,desc_a,info)
if (info == psb_success_) call a%cp_to_l(la)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
goto 9999
end if
do i=1,size(adiag)
if (adiag(i) /= szero) then
@ -202,10 +201,6 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
end if
end do
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
goto 9999
end if
! 1. Allocate Ptilde in sparse matrix form
@ -213,8 +208,8 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
call ptilde%mv_from(tmpcoo)
call ptilde%cscnv(info,type='csr')
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_) call la%cscnv(am3,info,type='csr',dupl=psb_dupl_add_)
if (info == psb_success_) call la%cscnv(da,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
@ -225,11 +220,10 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
call da%scal(adinv,info)
call psb_symbmm(da,ptilde,dap,info)
if (info == psb_success_) call psb_numbmm(da,ptilde,dap)
call psb_spspmm(da,ptilde,dap,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
goto 9999
end if
@ -240,8 +234,7 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=am4)
if (info == psb_success_) call am4%free()
call psb_symbmm(da,atmp,dadap,info)
call psb_numbmm(da,atmp,dadap)
call psb_spspmm(da,atmp,dadap,info)
call atmp%free()
! !$ write(0,*) 'Columns of AP',psb_sp_get_ncols(ap)
@ -249,7 +242,6 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
call dap%mv_to(csc_dap)
call dadap%mv_to(csc_dadap)
call csc_mat_col_prod(csc_dap,csc_dadap,omp,info)
call csc_mat_col_prod(csc_dadap,csc_dadap,oden,info)
call psb_sum(ictxt,omp)
@ -287,7 +279,7 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
!
! Build the filtered matrix Af from A
!
call a%cscnv(acsrf,info,dupl=psb_dupl_add_)
call la%cscnv(acsrf,info,dupl=psb_dupl_add_)
do i=1,nrow
tmp = szero
@ -327,23 +319,14 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
call af%mv_from(acsrf)
!
! Symbmm90 does the allocation for its result.
!
! op_prol = (I-w*D*Af)Ptilde
! Doing it this way means to consider diag(Af_i)
!
!
call psb_symbmm(af,ptilde,op_prol,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999
end if
call psb_numbmm(af,ptilde,op_prol)
call psb_spspmm(af,ptilde,op_prol,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
& 'Done SPSPMM 1'
else
!
! Build the smoothed prolongator using the original matrix
@ -363,19 +346,11 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
& write(debug_unit,*) me,' ',trim(name),&
& 'Done gather, going for SYMBMM 1'
!
! Symbmm90 does the allocation for its result.
!
! op_prol = (I-w*D*A)Ptilde
!
!
call psb_symbmm(am3,ptilde,op_prol,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999
end if
call psb_numbmm(am3,ptilde,op_prol)
call psb_spspmm(am3,ptilde,op_prol,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
@ -387,11 +362,11 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
! Ok, let's start over with the restrictor
!
call ptilde%transc(rtilde)
call a%cscnv(atmp,info,type='csr')
call la%cscnv(atmp,info,type='csr')
call psb_sphalo(atmp,desc_a,am4,info,&
& colcnv=.true.,rowscale=.true.)
nrt = am4%get_nrows()
call am4%csclip(atmp2,info,ione,nrt,ione,ncol)
call am4%csclip(atmp2,info,lone,nrt,lone,ncol)
call atmp2%cscnv(info,type='CSR')
if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=atmp2)
call am4%free()
@ -400,13 +375,12 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
! This is to compute the transpose. It ONLY works if the
! original A has a symmetric pattern.
call atmp%transc(atmp2)
call atmp2%csclip(dat,info,ione,nrow,ione,ncol)
call atmp2%csclip(dat,info,lone,nrow,lone,ncol)
call dat%cscnv(info,type='csr')
call dat%scal(adinv,info)
! Now for the product.
call psb_symbmm(dat,ptilde,datp,info)
if (info == psb_success_) call psb_numbmm(dat,ptilde,datp)
call psb_spspmm(dat,ptilde,datp,info)
call datp%clone(atmp2,info)
call psb_sphalo(atmp2,desc_a,am4,info,&
@ -484,8 +458,7 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
call rtilde%mv_from(tmpcoo)
call rtilde%cscnv(info,type='csr')
call psb_symbmm(rtilde,atmp,op_restr,info)
call psb_numbmm(rtilde,atmp,op_restr)
call psb_spspmm(rtilde,atmp,op_restr,info)
!
! Now we have to gather the halo of op_prol, and add it to itself
@ -527,16 +500,10 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
& write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd'
call psb_symbmm(a,tmp_prol,am3,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='symbmm 2')
goto 9999
end if
call psb_numbmm(a,tmp_prol,am3)
call psb_spspmm(la,tmp_prol,am3,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2'
& 'Done SPSPMM 2'
call psb_sphalo(am3,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.)
@ -552,8 +519,7 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
& write(debug_unit,*) me,' ',trim(name),&
& 'Done sphalo/ rwxtd'
call psb_symbmm(op_restr,am3,ac,info)
if (info == psb_success_) call psb_numbmm(op_restr,am3,ac)
call psb_spspmm(op_restr,am3,ac,info)
if (info == psb_success_) call am3%free()
if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_)
@ -581,11 +547,11 @@ contains
subroutine csc_mat_col_prod(a,b,v,info)
implicit none
type(psb_s_csc_sparse_mat), intent(in) :: a, b
type(psb_ls_csc_sparse_mat), intent(in) :: a, b
real(psb_spk_), intent(out) :: v(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i,j,k, nr, nc,iap,nra,ibp,nrb
integer(psb_lpk_) :: i,j,k, nr, nc,iap,nra,ibp,nrb
info = psb_success_
nc = a%get_ncols()
@ -609,11 +575,11 @@ contains
subroutine csr_mat_row_prod(a,b,v,info)
implicit none
type(psb_s_csr_sparse_mat), intent(in) :: a, b
type(psb_ls_csr_sparse_mat), intent(in) :: a, b
real(psb_spk_), intent(out) :: v(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i,j,k, nr, nc,iap,nca,ibp,ncb
integer(psb_lpk_) :: i,j,k, nr, nc,iap,nca,ibp,ncb
info = psb_success_
nr = a%get_nrows()
@ -637,12 +603,12 @@ contains
function sparse_srtd_dot(nv1,iv1,v1,nv2,iv2,v2) result(dot)
implicit none
integer(psb_ipk_), intent(in) :: nv1,nv2
integer(psb_ipk_), intent(in) :: iv1(:), iv2(:)
integer(psb_lpk_), intent(in) :: nv1,nv2
integer(psb_lpk_), intent(in) :: iv1(:), iv2(:)
real(psb_spk_), intent(in) :: v1(:),v2(:)
real(psb_spk_) :: dot
integer(psb_ipk_) :: i,j,k, ip1, ip2
integer(psb_lpk_) :: i,j,k, ip1, ip2
dot = szero
ip1 = 1
@ -665,7 +631,7 @@ contains
end function sparse_srtd_dot
subroutine local_dump(me,mat,name,header)
type(psb_sspmat_type), intent(in) :: mat
type(psb_lsspmat_type), intent(in) :: mat
integer(psb_ipk_), intent(in) :: me
character(len=*), intent(in) :: name
character(len=*), intent(in) :: header

@ -106,28 +106,30 @@ subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
! Arguments
type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:)
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_sml_parms), intent(inout) :: parms
type(psb_sspmat_type), intent(inout) :: op_prol
type(psb_sspmat_type), intent(out) :: ac,op_restr
type(psb_lsspmat_type), intent(inout) :: op_prol
type(psb_lsspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo
integer(psb_ipk_) :: ictxt, np, me, icomm, minfo
character(len=20) :: name
integer(psb_ipk_) :: ierr(5)
type(psb_s_coo_sparse_mat) :: ac_coo, tmpcoo
type(psb_s_csr_sparse_mat) :: acsr1, acsr2
integer(psb_ipk_) :: ierr(5)
type(psb_lsspmat_type) :: la
type(psb_ls_coo_sparse_mat) :: ac_coo, tmpcoo
type(psb_ls_csr_sparse_mat) :: acsr1, acsr2
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nrow, nglob, ncol, ntaggr, nzl, ip, &
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nzl, ip, &
& naggr, nzt, naggrm1, naggrp1, i, k
name='mld_aggrmat_nosmth_asb'
if(psb_get_errstatus().ne.0) return
info=psb_success_
name = 'mld_aggrmat_nosmth_asb'
info = psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
ictxt = desc_a%get_context()
icomm = desc_a%get_mpic()
@ -136,24 +138,23 @@ subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
nrow = desc_a%get_local_rows()
ncol = desc_a%get_local_cols()
naggr = nlaggr(me+1)
ntaggr = sum(nlaggr)
naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1))
call op_prol%cp_to(tmpcoo)
call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_)
call tmpcoo%transp()
nzl = tmpcoo%get_nzeros()
i=0
!
! Now we have to fix this. The only rows of tmpcoo/op_restr that are correct
! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
!
do k=1, nzl
nzl = tmpcoo%get_nzeros()
i = 0
do k = 1, nzl
if ((naggrm1 < tmpcoo%ia(k)) .and.(tmpcoo%ia(k) <= naggrp1)) then
i = i+1
tmpcoo%val(i) = tmpcoo%val(k)
@ -168,11 +169,11 @@ subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
if (info /= psb_success_) goto 9999
call a%cp_to(ac_coo)
call a%cp_to_l(la)
call la%mv_to(ac_coo)
nzt = ac_coo%get_nzeros()
k = 0
do i=1, nzt
k = 0
do i = 1, nzt
k = k + 1
ac_coo%ia(k) = ilaggr(ac_coo%ia(i))
ac_coo%ja(k) = ilaggr(ac_coo%ja(i))

@ -65,14 +65,7 @@
! specified by the user through mld_sprecinit and mld_zprecset.
! On output from this routine the entries of AC, op_prol, op_restr
! are still in "global numbering" mode; this is fixed in the calling routine
! mld_s_lev_aggrmat_asb.
!
! For more details see
! M. Brezina and P. Vanek, A black-box iterative solver based on a
! two-level Schwarz method, Computing, 63 (1999), 233-263.
! P. D'Ambra, D. di Serafino and S. Filippone, On the development of
! PSBLAS-based parallel two-level Schwarz preconditioners, Appl. Num. Math.
! 57 (2007), 1181-1196.
! aggregator%mat_asb.
!
!
! Arguments:
@ -117,37 +110,38 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
implicit none
! Arguments
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(inout) :: op_prol
type(psb_sspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_sml_parms), intent(inout) :: parms
type(psb_lsspmat_type), intent(inout) :: op_prol
type(psb_lsspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_) :: nrow, nglob, ncol, ntaggr, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw, err_act
integer(psb_ipk_) ::ictxt, np, me
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, ip, &
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw
integer(psb_ipk_) :: ictxt, np, me
character(len=20) :: name
type(psb_sspmat_type) :: am3, am4, tmp_prol
type(psb_s_coo_sparse_mat) :: tmpcoo
type(psb_s_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
type(psb_lsspmat_type) :: la, am3, am4, tmp_prol
type(psb_ls_coo_sparse_mat) :: tmpcoo
type(psb_ls_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
real(psb_spk_), allocatable :: adiag(:)
integer(psb_ipk_) :: ierr(5)
logical :: filter_mat
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: debug_level, debug_unit, err_act
integer(psb_ipk_), parameter :: ncmax=16
real(psb_spk_) :: anorm, omega, tmp, dg, theta
name='mld_aggrmat_smth_asb'
if(psb_get_errstatus().ne.0) return
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context()
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
@ -176,6 +170,7 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
& call psb_realloc(ncol,adiag,info)
if (info == psb_success_) &
& call psb_halo(adiag,desc_a,info)
if (info == psb_success_) call a%cp_to_l(la)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
@ -185,7 +180,7 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
! 1. Allocate Ptilde in sparse matrix form
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_) call la%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (info /= psb_success_) goto 9999
if (debug_level >= psb_debug_outer_) &
@ -198,7 +193,7 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
!
if (info == psb_success_) call acsr3%cp_to_fmt(acsrf,info)
do i=1,nrow
do i=1, nrow
tmp = szero
jd = -1
do j=acsrf%irp(i),acsrf%irp(i+1)-1
@ -342,7 +337,7 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
goto 9999
end if
call psb_spspmm(a,tmp_prol,am3,info)
call psb_spspmm(la,tmp_prol,am3,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 2')
goto 9999

@ -143,23 +143,25 @@ subroutine mld_z_dec_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,ac,op_p
type(mld_dml_parms), intent(inout) :: parms
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
type(psb_zspmat_type), intent(out) :: ac,op_restr
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(psb_lzspmat_type), intent(inout) :: op_prol
type(psb_lzspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
character(len=20) :: name
integer(psb_mpik_) :: ictxt, np, me
type(psb_z_coo_sparse_mat) :: acoo, bcoo
type(psb_z_csr_sparse_mat) :: acsr1
integer(psb_ipk_) :: nzl,ntaggr
integer(psb_ipk_) :: err_act
character(len=20) :: name
integer(psb_mpk_) :: ictxt, np, me
type(psb_lz_coo_sparse_mat) :: acoo, bcoo
type(psb_lz_csr_sparse_mat) :: acsr1
integer(psb_lpk_) :: nzl,ntaggr
integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: debug_level, debug_unit
name='mld_z_dec_aggregator_mat_asb'
if (psb_get_errstatus().ne.0) return
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
info = psb_success_

@ -81,20 +81,22 @@ subroutine mld_z_dec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_
type(mld_dml_parms), intent(inout) :: parms
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_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
type(psb_lzspmat_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_) :: ntaggr
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name
integer(psb_mpk_) :: ictxt, np, me
integer(psb_ipk_) :: err_act
integer(psb_lpk_) :: ntaggr
integer(psb_ipk_) :: debug_level, debug_unit
name='mld_z_dec_aggregator_tprol'
if (psb_get_errstatus().ne.0) return
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
info = psb_success_

@ -87,22 +87,24 @@ subroutine mld_z_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info)
! Arguments
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), allocatable, intent(inout) :: ilaggr(:),nlaggr(:)
type(psb_zspmat_type), intent(out) :: op_prol
integer(psb_lpk_), allocatable, intent(inout) :: ilaggr(:),nlaggr(:)
type(psb_lzspmat_type), intent(out) :: op_prol
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m,naggrm1, naggrp1, ntaggr
type(psb_z_coo_sparse_mat) :: tmpcoo
integer(psb_lpk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m,naggrm1, naggrp1, ntaggr
type(psb_lz_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
integer(psb_lpk_) :: nrow, ncol, n_ne
character(len=20) :: name, ch_err
if(psb_get_errstatus() /= 0) return
info=psb_success_
name = 'mld_map_to_tprol'
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!

@ -80,12 +80,13 @@ subroutine mld_z_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
type(psb_zspmat_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_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),&
& ideg(:), idxs(:), tmpaggr(:)
integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),&
& ideg(:), idxs(:)
integer(psb_lpk_), allocatable :: tmpaggr(:)
complex(psb_dpk_), allocatable :: val(:), diag(:)
integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, nc, naggr,i,j,m, nz, ilg, ii, ip
type(psb_z_csr_sparse_mat) :: acsr
@ -96,10 +97,12 @@ subroutine mld_z_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
integer(psb_ipk_) :: nrow, ncol, n_ne
character(len=20) :: name, ch_err
if (psb_get_errstatus() /= 0) return
info=psb_success_
name = 'mld_soc1_map_bld'
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!

@ -79,12 +79,13 @@ subroutine mld_z_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
type(psb_zspmat_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_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),&
& ideg(:), idxs(:), tmpaggr(:)
& ideg(:), idxs(:)
integer(psb_lpk_), allocatable :: tmpaggr(:)
complex(psb_dpk_), allocatable :: val(:), diag(:)
integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr,nc,naggr,i,j,m, nz, ilg, ii, ip, ip1,nzcnt
type(psb_z_csr_sparse_mat) :: acsr, muij, s_neigh
@ -96,10 +97,12 @@ subroutine mld_z_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
integer(psb_ipk_) :: nrow, ncol, n_ne
character(len=20) :: name, ch_err
if (psb_get_errstatus() /= 0) return
info=psb_success_
name = 'mld_soc2_map_bld'
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
!

@ -85,21 +85,24 @@ subroutine mld_z_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,
type(mld_dml_parms), intent(inout) :: parms
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_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
type(psb_lzspmat_type), intent(out) :: op_prol
integer(psb_ipk_), intent(out) :: info
! Local variables
type(psb_zspmat_type) :: atmp, atrans
character(len=20) :: name
integer(psb_mpik_) :: ictxt, np, me
integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: ntaggr, nr
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name
integer(psb_mpk_) :: ictxt, np, me
integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: nr
integer(psb_lpk_) :: ntaggr
integer(psb_ipk_) :: debug_level, debug_unit
name='mld_z_symdec_aggregator_tprol'
if (psb_get_errstatus().ne.0) return
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
info = psb_success_

@ -88,33 +88,35 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
implicit none
! Arguments
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(inout) :: op_prol
type(psb_zspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_dml_parms), intent(inout) :: parms
type(psb_lzspmat_type), intent(inout) :: op_prol
type(psb_lzspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_) :: nrow, nglob, ncol, ntaggr, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw, err_act
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, ip, &
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw
integer(psb_ipk_) ::ictxt, np, me
character(len=20) :: name
type(psb_zspmat_type) :: am3, am4,tmp_prol
type(psb_z_coo_sparse_mat) :: tmpcoo
type(psb_z_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
type(psb_lzspmat_type) :: am3, am4,tmp_prol, la
type(psb_lz_coo_sparse_mat) :: tmpcoo
type(psb_lz_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
complex(psb_dpk_), allocatable :: adiag(:)
integer(psb_ipk_) :: ierr(5)
logical :: filter_mat
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: debug_level, debug_unit, err_act
integer(psb_ipk_), parameter :: ncmax=16
real(psb_dpk_) :: anorm, omega, tmp, dg, theta
name='mld_aggrmat_biz_asb'
if(psb_get_errstatus().ne.0) return
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
@ -128,10 +130,8 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
ncol = desc_a%get_local_cols()
theta = parms%aggr_thresh
naggr = nlaggr(me+1)
ntaggr = sum(nlaggr)
filter_mat = (parms%aggr_filter == mld_filter_mat_)
! naggr: number of local aggregates
@ -143,8 +143,8 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
& call psb_realloc(ncol,adiag,info)
if (info == psb_success_) &
& call psb_halo(adiag,desc_a,info)
if(info /= psb_success_) then
if (info == psb_success_) call a%cp_to_l(la)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
goto 9999
end if
@ -152,7 +152,7 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
! 1. Allocate Ptilde in sparse matrix form
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_) call la%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (info /= psb_success_) goto 9999
if (debug_level >= psb_debug_outer_) &
@ -163,9 +163,9 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
!
! Build the filtered matrix Af from A
!
if (info == psb_success_) call a%cscnv(acsrf,info,dupl=psb_dupl_add_)
if (info == psb_success_) call acsr3%cp_to_fmt(acsrf,info)
do i=1,nrow
do i=1, nrow
tmp = zzero
jd = -1
do j=acsrf%irp(i),acsrf%irp(i+1)-1
@ -271,17 +271,15 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
! Doing it this way means to consider diag(Af_i)
!
!
call psb_symbmm(acsrf,ptilde,acsr1,info)
call psb_spspmm(acsrf,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
goto 9999
end if
call psb_numbmm(acsrf,ptilde,acsr1)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
& 'Done SPSPMM 1'
else
!
@ -307,17 +305,15 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
! Doing it this way means to consider diag(A_i)
!
!
call psb_symbmm(acsr3,ptilde,acsr1,info)
call psb_spspmm(acsr3,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
goto 9999
end if
call psb_numbmm(acsr3,ptilde,acsr1)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
& 'Done SPSPMM 1'
end if
call ptilde%free()
@ -331,16 +327,15 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
goto 9999
end if
call psb_symbmm(a,tmp_prol,am3,info)
call psb_spspmm(la,tmp_prol,am3,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2')
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 2')
goto 9999
end if
call psb_numbmm(a,tmp_prol,am3)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2',parms%aggr_prol, mld_smooth_prol_
& 'Done SPSPMM 2',parms%aggr_prol, mld_smooth_prol_
call tmp_prol%transp(op_restr)
if (debug_level >= psb_debug_outer_) &
@ -353,23 +348,18 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting symbmm 3'
call psb_symbmm(op_restr,am3,ac,info)
if (info == psb_success_) call psb_numbmm(op_restr,am3,ac)
& 'starting spspmm 3'
call psb_spspmm(op_restr,am3,ac,info)
if (info == psb_success_) call am3%free()
if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_)
if (info == psb_success_) call ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Build b = op_restr x am3')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done smooth_aggregate '

@ -114,37 +114,38 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
! Arguments
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(inout) :: op_prol
type(psb_zspmat_type), intent(out) :: ac,op_restr
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_dml_parms), intent(inout) :: parms
type(psb_lzspmat_type), intent(inout) :: op_prol
type(psb_lzspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_), allocatable :: nzbr(:), idisp(:)
integer(psb_ipk_) :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrt, err_act
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrt
integer(psb_ipk_) :: ictxt,np,me, icomm
character(len=20) :: name
type(psb_zspmat_type) :: af, ptilde, rtilde, atran, atp, atdatp
type(psb_zspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da
type(psb_zspmat_type) :: dat, datp, datdatp, atmp3, tmp_prol
type(psb_z_coo_sparse_mat) :: tmpcoo
type(psb_z_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf
type(psb_z_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc
type(psb_lzspmat_type) :: la, af, ptilde, rtilde, atran, atp, atdatp
type(psb_lzspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da
type(psb_lzspmat_type) :: dat, datp, datdatp, atmp3, tmp_prol
type(psb_lz_coo_sparse_mat) :: tmpcoo
type(psb_lz_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf
type(psb_lz_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc
complex(psb_dpk_), allocatable :: adiag(:), adinv(:)
complex(psb_dpk_), allocatable :: omf(:), omp(:), omi(:), oden(:)
logical :: filter_mat
integer(psb_ipk_) :: ierr(5)
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: debug_level, debug_unit, err_act
integer(psb_ipk_), parameter :: ncmax=16
real(psb_dpk_) :: anorm, theta
complex(psb_dpk_) :: tmp, alpha, beta, ommx
name='mld_aggrmat_minnrg'
if(psb_get_errstatus().ne.0) return
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
@ -163,13 +164,6 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
naggr = nlaggr(me+1)
ntaggr = sum(nlaggr)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_; ierr(1)=2*np;
call psb_errpush(info,name,i_err=ierr,a_err='integer')
goto 9999
end if
naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1))
@ -193,6 +187,11 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
& call psb_realloc(ncol,adiag,info)
if (info == psb_success_) &
& call psb_halo(adiag,desc_a,info)
if (info == psb_success_) call a%cp_to_l(la)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
goto 9999
end if
do i=1,size(adiag)
if (adiag(i) /= zzero) then
@ -202,10 +201,6 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
end if
end do
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
goto 9999
end if
! 1. Allocate Ptilde in sparse matrix form
@ -213,8 +208,8 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
call ptilde%mv_from(tmpcoo)
call ptilde%cscnv(info,type='csr')
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_) call la%cscnv(am3,info,type='csr',dupl=psb_dupl_add_)
if (info == psb_success_) call la%cscnv(da,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
@ -225,11 +220,10 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
call da%scal(adinv,info)
call psb_symbmm(da,ptilde,dap,info)
if (info == psb_success_) call psb_numbmm(da,ptilde,dap)
call psb_spspmm(da,ptilde,dap,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
goto 9999
end if
@ -240,8 +234,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=am4)
if (info == psb_success_) call am4%free()
call psb_symbmm(da,atmp,dadap,info)
call psb_numbmm(da,atmp,dadap)
call psb_spspmm(da,atmp,dadap,info)
call atmp%free()
! !$ write(0,*) 'Columns of AP',psb_sp_get_ncols(ap)
@ -249,7 +242,6 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
call dap%mv_to(csc_dap)
call dadap%mv_to(csc_dadap)
call csc_mat_col_prod(csc_dap,csc_dadap,omp,info)
call csc_mat_col_prod(csc_dadap,csc_dadap,oden,info)
call psb_sum(ictxt,omp)
@ -287,7 +279,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
!
! Build the filtered matrix Af from A
!
call a%cscnv(acsrf,info,dupl=psb_dupl_add_)
call la%cscnv(acsrf,info,dupl=psb_dupl_add_)
do i=1,nrow
tmp = zzero
@ -327,23 +319,14 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
call af%mv_from(acsrf)
!
! Symbmm90 does the allocation for its result.
!
! op_prol = (I-w*D*Af)Ptilde
! Doing it this way means to consider diag(Af_i)
!
!
call psb_symbmm(af,ptilde,op_prol,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999
end if
call psb_numbmm(af,ptilde,op_prol)
call psb_spspmm(af,ptilde,op_prol,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
& 'Done SPSPMM 1'
else
!
! Build the smoothed prolongator using the original matrix
@ -363,19 +346,11 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
& write(debug_unit,*) me,' ',trim(name),&
& 'Done gather, going for SYMBMM 1'
!
! Symbmm90 does the allocation for its result.
!
! op_prol = (I-w*D*A)Ptilde
!
!
call psb_symbmm(am3,ptilde,op_prol,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999
end if
call psb_numbmm(am3,ptilde,op_prol)
call psb_spspmm(am3,ptilde,op_prol,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
@ -387,11 +362,11 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
! Ok, let's start over with the restrictor
!
call ptilde%transc(rtilde)
call a%cscnv(atmp,info,type='csr')
call la%cscnv(atmp,info,type='csr')
call psb_sphalo(atmp,desc_a,am4,info,&
& colcnv=.true.,rowscale=.true.)
nrt = am4%get_nrows()
call am4%csclip(atmp2,info,ione,nrt,ione,ncol)
call am4%csclip(atmp2,info,lone,nrt,lone,ncol)
call atmp2%cscnv(info,type='CSR')
if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=atmp2)
call am4%free()
@ -400,13 +375,12 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
! This is to compute the transpose. It ONLY works if the
! original A has a symmetric pattern.
call atmp%transc(atmp2)
call atmp2%csclip(dat,info,ione,nrow,ione,ncol)
call atmp2%csclip(dat,info,lone,nrow,lone,ncol)
call dat%cscnv(info,type='csr')
call dat%scal(adinv,info)
! Now for the product.
call psb_symbmm(dat,ptilde,datp,info)
if (info == psb_success_) call psb_numbmm(dat,ptilde,datp)
call psb_spspmm(dat,ptilde,datp,info)
call datp%clone(atmp2,info)
call psb_sphalo(atmp2,desc_a,am4,info,&
@ -484,8 +458,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
call rtilde%mv_from(tmpcoo)
call rtilde%cscnv(info,type='csr')
call psb_symbmm(rtilde,atmp,op_restr,info)
call psb_numbmm(rtilde,atmp,op_restr)
call psb_spspmm(rtilde,atmp,op_restr,info)
!
! Now we have to gather the halo of op_prol, and add it to itself
@ -527,16 +500,10 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
& write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd'
call psb_symbmm(a,tmp_prol,am3,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='symbmm 2')
goto 9999
end if
call psb_numbmm(a,tmp_prol,am3)
call psb_spspmm(la,tmp_prol,am3,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2'
& 'Done SPSPMM 2'
call psb_sphalo(am3,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.)
@ -552,8 +519,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
& write(debug_unit,*) me,' ',trim(name),&
& 'Done sphalo/ rwxtd'
call psb_symbmm(op_restr,am3,ac,info)
if (info == psb_success_) call psb_numbmm(op_restr,am3,ac)
call psb_spspmm(op_restr,am3,ac,info)
if (info == psb_success_) call am3%free()
if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_)
@ -581,11 +547,11 @@ contains
subroutine csc_mat_col_prod(a,b,v,info)
implicit none
type(psb_z_csc_sparse_mat), intent(in) :: a, b
type(psb_lz_csc_sparse_mat), intent(in) :: a, b
complex(psb_dpk_), intent(out) :: v(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i,j,k, nr, nc,iap,nra,ibp,nrb
integer(psb_lpk_) :: i,j,k, nr, nc,iap,nra,ibp,nrb
info = psb_success_
nc = a%get_ncols()
@ -609,11 +575,11 @@ contains
subroutine csr_mat_row_prod(a,b,v,info)
implicit none
type(psb_z_csr_sparse_mat), intent(in) :: a, b
type(psb_lz_csr_sparse_mat), intent(in) :: a, b
complex(psb_dpk_), intent(out) :: v(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i,j,k, nr, nc,iap,nca,ibp,ncb
integer(psb_lpk_) :: i,j,k, nr, nc,iap,nca,ibp,ncb
info = psb_success_
nr = a%get_nrows()
@ -637,12 +603,12 @@ contains
function sparse_srtd_dot(nv1,iv1,v1,nv2,iv2,v2) result(dot)
implicit none
integer(psb_ipk_), intent(in) :: nv1,nv2
integer(psb_ipk_), intent(in) :: iv1(:), iv2(:)
integer(psb_lpk_), intent(in) :: nv1,nv2
integer(psb_lpk_), intent(in) :: iv1(:), iv2(:)
complex(psb_dpk_), intent(in) :: v1(:),v2(:)
complex(psb_dpk_) :: dot
integer(psb_ipk_) :: i,j,k, ip1, ip2
integer(psb_lpk_) :: i,j,k, ip1, ip2
dot = zzero
ip1 = 1
@ -665,7 +631,7 @@ contains
end function sparse_srtd_dot
subroutine local_dump(me,mat,name,header)
type(psb_zspmat_type), intent(in) :: mat
type(psb_lzspmat_type), intent(in) :: mat
integer(psb_ipk_), intent(in) :: me
character(len=*), intent(in) :: name
character(len=*), intent(in) :: header

@ -106,28 +106,30 @@ subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
! Arguments
type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:)
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_dml_parms), intent(inout) :: parms
type(psb_zspmat_type), intent(inout) :: op_prol
type(psb_zspmat_type), intent(out) :: ac,op_restr
type(psb_lzspmat_type), intent(inout) :: op_prol
type(psb_lzspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo
integer(psb_ipk_) :: ictxt, np, me, icomm, minfo
character(len=20) :: name
integer(psb_ipk_) :: ierr(5)
type(psb_z_coo_sparse_mat) :: ac_coo, tmpcoo
type(psb_z_csr_sparse_mat) :: acsr1, acsr2
integer(psb_ipk_) :: ierr(5)
type(psb_lzspmat_type) :: la
type(psb_lz_coo_sparse_mat) :: ac_coo, tmpcoo
type(psb_lz_csr_sparse_mat) :: acsr1, acsr2
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nrow, nglob, ncol, ntaggr, nzl, ip, &
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nzl, ip, &
& naggr, nzt, naggrm1, naggrp1, i, k
name='mld_aggrmat_nosmth_asb'
if(psb_get_errstatus().ne.0) return
info=psb_success_
name = 'mld_aggrmat_nosmth_asb'
info = psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
ictxt = desc_a%get_context()
icomm = desc_a%get_mpic()
@ -136,24 +138,23 @@ subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
nrow = desc_a%get_local_rows()
ncol = desc_a%get_local_cols()
naggr = nlaggr(me+1)
ntaggr = sum(nlaggr)
naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1))
call op_prol%cp_to(tmpcoo)
call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_)
call tmpcoo%transp()
nzl = tmpcoo%get_nzeros()
i=0
!
! Now we have to fix this. The only rows of tmpcoo/op_restr that are correct
! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
!
do k=1, nzl
nzl = tmpcoo%get_nzeros()
i = 0
do k = 1, nzl
if ((naggrm1 < tmpcoo%ia(k)) .and.(tmpcoo%ia(k) <= naggrp1)) then
i = i+1
tmpcoo%val(i) = tmpcoo%val(k)
@ -168,11 +169,11 @@ subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
if (info /= psb_success_) goto 9999
call a%cp_to(ac_coo)
call a%cp_to_l(la)
call la%mv_to(ac_coo)
nzt = ac_coo%get_nzeros()
k = 0
do i=1, nzt
k = 0
do i = 1, nzt
k = k + 1
ac_coo%ia(k) = ilaggr(ac_coo%ia(i))
ac_coo%ja(k) = ilaggr(ac_coo%ja(i))

@ -65,14 +65,7 @@
! specified by the user through mld_zprecinit and mld_zprecset.
! On output from this routine the entries of AC, op_prol, op_restr
! are still in "global numbering" mode; this is fixed in the calling routine
! mld_z_lev_aggrmat_asb.
!
! For more details see
! M. Brezina and P. Vanek, A black-box iterative solver based on a
! two-level Schwarz method, Computing, 63 (1999), 233-263.
! P. D'Ambra, D. di Serafino and S. Filippone, On the development of
! PSBLAS-based parallel two-level Schwarz preconditioners, Appl. Num. Math.
! 57 (2007), 1181-1196.
! aggregator%mat_asb.
!
!
! Arguments:
@ -117,37 +110,38 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
implicit none
! Arguments
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(inout) :: op_prol
type(psb_zspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_dml_parms), intent(inout) :: parms
type(psb_lzspmat_type), intent(inout) :: op_prol
type(psb_lzspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_) :: nrow, nglob, ncol, ntaggr, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw, err_act
integer(psb_ipk_) ::ictxt, np, me
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, ip, &
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw
integer(psb_ipk_) :: ictxt, np, me
character(len=20) :: name
type(psb_zspmat_type) :: am3, am4, tmp_prol
type(psb_z_coo_sparse_mat) :: tmpcoo
type(psb_z_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
type(psb_lzspmat_type) :: la, am3, am4, tmp_prol
type(psb_lz_coo_sparse_mat) :: tmpcoo
type(psb_lz_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
complex(psb_dpk_), allocatable :: adiag(:)
integer(psb_ipk_) :: ierr(5)
logical :: filter_mat
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: debug_level, debug_unit, err_act
integer(psb_ipk_), parameter :: ncmax=16
real(psb_dpk_) :: anorm, omega, tmp, dg, theta
name='mld_aggrmat_smth_asb'
if(psb_get_errstatus().ne.0) return
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context()
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
@ -176,6 +170,7 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
& call psb_realloc(ncol,adiag,info)
if (info == psb_success_) &
& call psb_halo(adiag,desc_a,info)
if (info == psb_success_) call a%cp_to_l(la)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
@ -185,7 +180,7 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
! 1. Allocate Ptilde in sparse matrix form
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_) call la%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (info /= psb_success_) goto 9999
if (debug_level >= psb_debug_outer_) &
@ -198,7 +193,7 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
!
if (info == psb_success_) call acsr3%cp_to_fmt(acsrf,info)
do i=1,nrow
do i=1, nrow
tmp = zzero
jd = -1
do j=acsrf%irp(i),acsrf%irp(i+1)-1
@ -342,7 +337,7 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
goto 9999
end if
call psb_spspmm(a,tmp_prol,am3,info)
call psb_spspmm(la,tmp_prol,am3,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 2')
goto 9999

@ -53,10 +53,12 @@ subroutine mld_c_base_onelev_build(lv,info,amold,vmold,imold)
character(len=20) :: name, ch_err
name = 'mld_onelev_build'
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
if (.not.associated(lv%base_desc)) then

@ -95,24 +95,29 @@ subroutine mld_c_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
class(mld_c_onelev_type), intent(inout), target :: lv
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_lpk_), intent(inout) :: nlaggr(:)
integer(psb_lpk_), intent(inout) :: ilaggr(:)
type(psb_lcspmat_type), intent(inout) :: op_prol
integer(psb_ipk_), intent(out) :: info
! Local variables
character(len=24) :: name
integer(psb_mpik_) :: ictxt, np, me
integer(psb_ipk_) :: 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
type(psb_lcspmat_type) :: lac, lac1, op_restr
type(psb_cspmat_type) :: ac, iop_restr, iop_prol
type(psb_lc_coo_sparse_mat) :: acoo, bcoo
type(psb_lc_csr_sparse_mat) :: acsr1
integer(psb_lpk_) :: ntaggr, nr, nc
integer(psb_ipk_) :: nzl, inl
integer(psb_ipk_) :: debug_level, debug_unit
name='mld_c_onelev_mat_asb'
if (psb_get_errstatus().ne.0) return
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
info = psb_success_
@ -137,7 +142,7 @@ subroutine mld_c_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
! the mapping defined by mld_aggrmap_bld and applying the aggregation
! algorithm specified by lv%iprcparm(mld_aggr_prol_)
!
call lv%aggr%mat_asb(lv%parms,a,desc_a,ilaggr,nlaggr,ac,op_prol,op_restr,info)
call lv%aggr%mat_asb(lv%parms,a,desc_a,ilaggr,nlaggr,lac,op_prol,op_restr,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_asb')
@ -153,10 +158,16 @@ subroutine mld_c_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
case(mld_distr_mat_)
call ac%mv_to(bcoo)
call lac%mv_to(bcoo)
nzl = bcoo%get_nzeros()
if (info == psb_success_) call psb_cdall(ictxt,lv%desc_ac,info,nl=nlaggr(me+1))
inl = nlaggr(me+1)
if (inl < nlaggr(me+1)) then
info = psb_err_bad_int_cnv_
call psb_errpush(info,name,&
& l_err=(/nlaggr(me+1),inl*1_psb_lpk_/))
goto 9999
end if
if (info == psb_success_) call psb_cdall(ictxt,lv%desc_ac,info,nl=inl)
if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,lv%desc_ac,info)
if (info == psb_success_) call psb_cdasb(lv%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),lv%desc_ac,info,iact='I')
@ -169,7 +180,8 @@ subroutine mld_c_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call lv%ac%mv_from(bcoo)
call lac%mv_from(bcoo)
call lv%ac%mv_from_l(lac)
call lv%ac%set_nrows(lv%desc_ac%get_local_rows())
call lv%ac%set_ncols(lv%desc_ac%get_local_cols())
@ -190,7 +202,8 @@ subroutine mld_c_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
end if
call op_prol%mv_from(acsr1)
endif
call op_prol%set_ncols(lv%desc_ac%get_local_cols())
nc = lv%desc_ac%get_local_cols()
call op_prol%set_ncols(nc)
if (np>1) then
call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_)
@ -209,7 +222,8 @@ subroutine mld_c_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
!
! Clip to local rows.
!
call op_restr%set_nrows(lv%desc_ac%get_local_rows())
nr = lv%desc_ac%get_local_rows()
call op_restr%set_nrows(nr)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
@ -218,11 +232,12 @@ subroutine mld_c_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
case(mld_repl_mat_)
!
!
call psb_cdall(ictxt,lv%desc_ac,info,mg=ntaggr,repl=.true.)
if (info == psb_success_) call psb_cdasb(lv%desc_ac,info)
if (info == psb_success_) &
& call psb_gather(lv%ac,ac,lv%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
& call psb_gather(lac1,lac,lv%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
call lv%ac%mv_from_l(lac1)
if (info /= psb_success_) goto 9999
case default
@ -242,11 +257,12 @@ subroutine mld_c_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
! op_restr => PR^T i.e. restriction operator
! op_prol => PR i.e. prolongation operator
!
call iop_restr%mv_from_l(op_restr)
call iop_prol%mv_from_l(op_prol)
lv%map = psb_linmap(psb_map_aggr_,desc_a,&
& lv%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
& lv%desc_ac,iop_restr,iop_prol,ilaggr,nlaggr)
if (info == psb_success_) call iop_prol%free()
if (info == psb_success_) call iop_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999

@ -53,10 +53,12 @@ subroutine mld_d_base_onelev_build(lv,info,amold,vmold,imold)
character(len=20) :: name, ch_err
name = 'mld_onelev_build'
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
if (.not.associated(lv%base_desc)) then

@ -95,24 +95,29 @@ subroutine mld_d_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
class(mld_d_onelev_type), intent(inout), target :: lv
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(inout) :: ilaggr(:),nlaggr(:)
type(psb_dspmat_type), intent(inout) :: op_prol
integer(psb_lpk_), intent(inout) :: nlaggr(:)
integer(psb_lpk_), intent(inout) :: ilaggr(:)
type(psb_ldspmat_type), intent(inout) :: op_prol
integer(psb_ipk_), intent(out) :: info
! Local variables
character(len=24) :: name
integer(psb_mpik_) :: ictxt, np, me
integer(psb_ipk_) :: ictxt, np, me
integer(psb_ipk_) :: err_act
type(psb_dspmat_type) :: ac, op_restr
type(psb_d_coo_sparse_mat) :: acoo, bcoo
type(psb_d_csr_sparse_mat) :: acsr1
integer(psb_ipk_) :: nzl, ntaggr
type(psb_ldspmat_type) :: lac, lac1, op_restr
type(psb_dspmat_type) :: ac, iop_restr, iop_prol
type(psb_ld_coo_sparse_mat) :: acoo, bcoo
type(psb_ld_csr_sparse_mat) :: acsr1
integer(psb_lpk_) :: ntaggr, nr, nc
integer(psb_ipk_) :: nzl, inl
integer(psb_ipk_) :: debug_level, debug_unit
name='mld_d_onelev_mat_asb'
if (psb_get_errstatus().ne.0) return
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
info = psb_success_
@ -137,7 +142,7 @@ subroutine mld_d_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
! the mapping defined by mld_aggrmap_bld and applying the aggregation
! algorithm specified by lv%iprcparm(mld_aggr_prol_)
!
call lv%aggr%mat_asb(lv%parms,a,desc_a,ilaggr,nlaggr,ac,op_prol,op_restr,info)
call lv%aggr%mat_asb(lv%parms,a,desc_a,ilaggr,nlaggr,lac,op_prol,op_restr,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_asb')
@ -153,10 +158,16 @@ subroutine mld_d_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
case(mld_distr_mat_)
call ac%mv_to(bcoo)
call lac%mv_to(bcoo)
nzl = bcoo%get_nzeros()
if (info == psb_success_) call psb_cdall(ictxt,lv%desc_ac,info,nl=nlaggr(me+1))
inl = nlaggr(me+1)
if (inl < nlaggr(me+1)) then
info = psb_err_bad_int_cnv_
call psb_errpush(info,name,&
& l_err=(/nlaggr(me+1),inl*1_psb_lpk_/))
goto 9999
end if
if (info == psb_success_) call psb_cdall(ictxt,lv%desc_ac,info,nl=inl)
if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,lv%desc_ac,info)
if (info == psb_success_) call psb_cdasb(lv%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),lv%desc_ac,info,iact='I')
@ -169,7 +180,8 @@ subroutine mld_d_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call lv%ac%mv_from(bcoo)
call lac%mv_from(bcoo)
call lv%ac%mv_from_l(lac)
call lv%ac%set_nrows(lv%desc_ac%get_local_rows())
call lv%ac%set_ncols(lv%desc_ac%get_local_cols())
@ -190,7 +202,8 @@ subroutine mld_d_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
end if
call op_prol%mv_from(acsr1)
endif
call op_prol%set_ncols(lv%desc_ac%get_local_cols())
nc = lv%desc_ac%get_local_cols()
call op_prol%set_ncols(nc)
if (np>1) then
call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_)
@ -209,7 +222,8 @@ subroutine mld_d_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
!
! Clip to local rows.
!
call op_restr%set_nrows(lv%desc_ac%get_local_rows())
nr = lv%desc_ac%get_local_rows()
call op_restr%set_nrows(nr)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
@ -218,11 +232,12 @@ subroutine mld_d_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
case(mld_repl_mat_)
!
!
call psb_cdall(ictxt,lv%desc_ac,info,mg=ntaggr,repl=.true.)
if (info == psb_success_) call psb_cdasb(lv%desc_ac,info)
if (info == psb_success_) &
& call psb_gather(lv%ac,ac,lv%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
& call psb_gather(lac1,lac,lv%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
call lv%ac%mv_from_l(lac1)
if (info /= psb_success_) goto 9999
case default
@ -242,11 +257,12 @@ subroutine mld_d_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
! op_restr => PR^T i.e. restriction operator
! op_prol => PR i.e. prolongation operator
!
call iop_restr%mv_from_l(op_restr)
call iop_prol%mv_from_l(op_prol)
lv%map = psb_linmap(psb_map_aggr_,desc_a,&
& lv%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
& lv%desc_ac,iop_restr,iop_prol,ilaggr,nlaggr)
if (info == psb_success_) call iop_prol%free()
if (info == psb_success_) call iop_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999

@ -53,10 +53,12 @@ subroutine mld_s_base_onelev_build(lv,info,amold,vmold,imold)
character(len=20) :: name, ch_err
name = 'mld_onelev_build'
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
if (.not.associated(lv%base_desc)) then

@ -95,24 +95,29 @@ subroutine mld_s_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
class(mld_s_onelev_type), intent(inout), target :: lv
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_lpk_), intent(inout) :: nlaggr(:)
integer(psb_lpk_), intent(inout) :: ilaggr(:)
type(psb_lsspmat_type), intent(inout) :: op_prol
integer(psb_ipk_), intent(out) :: info
! Local variables
character(len=24) :: name
integer(psb_mpik_) :: ictxt, np, me
integer(psb_ipk_) :: 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
type(psb_lsspmat_type) :: lac, lac1, op_restr
type(psb_sspmat_type) :: ac, iop_restr, iop_prol
type(psb_ls_coo_sparse_mat) :: acoo, bcoo
type(psb_ls_csr_sparse_mat) :: acsr1
integer(psb_lpk_) :: ntaggr, nr, nc
integer(psb_ipk_) :: nzl, inl
integer(psb_ipk_) :: debug_level, debug_unit
name='mld_s_onelev_mat_asb'
if (psb_get_errstatus().ne.0) return
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
info = psb_success_
@ -137,7 +142,7 @@ subroutine mld_s_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
! the mapping defined by mld_aggrmap_bld and applying the aggregation
! algorithm specified by lv%iprcparm(mld_aggr_prol_)
!
call lv%aggr%mat_asb(lv%parms,a,desc_a,ilaggr,nlaggr,ac,op_prol,op_restr,info)
call lv%aggr%mat_asb(lv%parms,a,desc_a,ilaggr,nlaggr,lac,op_prol,op_restr,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_asb')
@ -153,10 +158,16 @@ subroutine mld_s_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
case(mld_distr_mat_)
call ac%mv_to(bcoo)
call lac%mv_to(bcoo)
nzl = bcoo%get_nzeros()
if (info == psb_success_) call psb_cdall(ictxt,lv%desc_ac,info,nl=nlaggr(me+1))
inl = nlaggr(me+1)
if (inl < nlaggr(me+1)) then
info = psb_err_bad_int_cnv_
call psb_errpush(info,name,&
& l_err=(/nlaggr(me+1),inl*1_psb_lpk_/))
goto 9999
end if
if (info == psb_success_) call psb_cdall(ictxt,lv%desc_ac,info,nl=inl)
if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,lv%desc_ac,info)
if (info == psb_success_) call psb_cdasb(lv%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),lv%desc_ac,info,iact='I')
@ -169,7 +180,8 @@ subroutine mld_s_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call lv%ac%mv_from(bcoo)
call lac%mv_from(bcoo)
call lv%ac%mv_from_l(lac)
call lv%ac%set_nrows(lv%desc_ac%get_local_rows())
call lv%ac%set_ncols(lv%desc_ac%get_local_cols())
@ -190,7 +202,8 @@ subroutine mld_s_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
end if
call op_prol%mv_from(acsr1)
endif
call op_prol%set_ncols(lv%desc_ac%get_local_cols())
nc = lv%desc_ac%get_local_cols()
call op_prol%set_ncols(nc)
if (np>1) then
call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_)
@ -209,7 +222,8 @@ subroutine mld_s_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
!
! Clip to local rows.
!
call op_restr%set_nrows(lv%desc_ac%get_local_rows())
nr = lv%desc_ac%get_local_rows()
call op_restr%set_nrows(nr)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
@ -218,11 +232,12 @@ subroutine mld_s_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
case(mld_repl_mat_)
!
!
call psb_cdall(ictxt,lv%desc_ac,info,mg=ntaggr,repl=.true.)
if (info == psb_success_) call psb_cdasb(lv%desc_ac,info)
if (info == psb_success_) &
& call psb_gather(lv%ac,ac,lv%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
& call psb_gather(lac1,lac,lv%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
call lv%ac%mv_from_l(lac1)
if (info /= psb_success_) goto 9999
case default
@ -242,11 +257,12 @@ subroutine mld_s_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
! op_restr => PR^T i.e. restriction operator
! op_prol => PR i.e. prolongation operator
!
call iop_restr%mv_from_l(op_restr)
call iop_prol%mv_from_l(op_prol)
lv%map = psb_linmap(psb_map_aggr_,desc_a,&
& lv%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
& lv%desc_ac,iop_restr,iop_prol,ilaggr,nlaggr)
if (info == psb_success_) call iop_prol%free()
if (info == psb_success_) call iop_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999

@ -53,10 +53,12 @@ subroutine mld_z_base_onelev_build(lv,info,amold,vmold,imold)
character(len=20) :: name, ch_err
name = 'mld_onelev_build'
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
if (.not.associated(lv%base_desc)) then

@ -95,24 +95,29 @@ subroutine mld_z_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
class(mld_z_onelev_type), intent(inout), target :: lv
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_lpk_), intent(inout) :: nlaggr(:)
integer(psb_lpk_), intent(inout) :: ilaggr(:)
type(psb_lzspmat_type), intent(inout) :: op_prol
integer(psb_ipk_), intent(out) :: info
! Local variables
character(len=24) :: name
integer(psb_mpik_) :: ictxt, np, me
integer(psb_ipk_) :: 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
type(psb_lzspmat_type) :: lac, lac1, op_restr
type(psb_zspmat_type) :: ac, iop_restr, iop_prol
type(psb_lz_coo_sparse_mat) :: acoo, bcoo
type(psb_lz_csr_sparse_mat) :: acsr1
integer(psb_lpk_) :: ntaggr, nr, nc
integer(psb_ipk_) :: nzl, inl
integer(psb_ipk_) :: debug_level, debug_unit
name='mld_z_onelev_mat_asb'
if (psb_get_errstatus().ne.0) return
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
info = psb_success_
@ -137,7 +142,7 @@ subroutine mld_z_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
! the mapping defined by mld_aggrmap_bld and applying the aggregation
! algorithm specified by lv%iprcparm(mld_aggr_prol_)
!
call lv%aggr%mat_asb(lv%parms,a,desc_a,ilaggr,nlaggr,ac,op_prol,op_restr,info)
call lv%aggr%mat_asb(lv%parms,a,desc_a,ilaggr,nlaggr,lac,op_prol,op_restr,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_asb')
@ -153,10 +158,16 @@ subroutine mld_z_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
case(mld_distr_mat_)
call ac%mv_to(bcoo)
call lac%mv_to(bcoo)
nzl = bcoo%get_nzeros()
if (info == psb_success_) call psb_cdall(ictxt,lv%desc_ac,info,nl=nlaggr(me+1))
inl = nlaggr(me+1)
if (inl < nlaggr(me+1)) then
info = psb_err_bad_int_cnv_
call psb_errpush(info,name,&
& l_err=(/nlaggr(me+1),inl*1_psb_lpk_/))
goto 9999
end if
if (info == psb_success_) call psb_cdall(ictxt,lv%desc_ac,info,nl=inl)
if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,lv%desc_ac,info)
if (info == psb_success_) call psb_cdasb(lv%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),lv%desc_ac,info,iact='I')
@ -169,7 +180,8 @@ subroutine mld_z_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call lv%ac%mv_from(bcoo)
call lac%mv_from(bcoo)
call lv%ac%mv_from_l(lac)
call lv%ac%set_nrows(lv%desc_ac%get_local_rows())
call lv%ac%set_ncols(lv%desc_ac%get_local_cols())
@ -190,7 +202,8 @@ subroutine mld_z_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
end if
call op_prol%mv_from(acsr1)
endif
call op_prol%set_ncols(lv%desc_ac%get_local_cols())
nc = lv%desc_ac%get_local_cols()
call op_prol%set_ncols(nc)
if (np>1) then
call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_)
@ -209,7 +222,8 @@ subroutine mld_z_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
!
! Clip to local rows.
!
call op_restr%set_nrows(lv%desc_ac%get_local_rows())
nr = lv%desc_ac%get_local_rows()
call op_restr%set_nrows(nr)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
@ -218,11 +232,12 @@ subroutine mld_z_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
case(mld_repl_mat_)
!
!
call psb_cdall(ictxt,lv%desc_ac,info,mg=ntaggr,repl=.true.)
if (info == psb_success_) call psb_cdasb(lv%desc_ac,info)
if (info == psb_success_) &
& call psb_gather(lv%ac,ac,lv%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
& call psb_gather(lac1,lac,lv%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
call lv%ac%mv_from_l(lac1)
if (info /= psb_success_) goto 9999
case default
@ -242,11 +257,12 @@ subroutine mld_z_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
! op_restr => PR^T i.e. restriction operator
! op_prol => PR i.e. prolongation operator
!
call iop_restr%mv_from_l(op_restr)
call iop_prol%mv_from_l(op_prol)
lv%map = psb_linmap(psb_map_aggr_,desc_a,&
& lv%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
& lv%desc_ac,iop_restr,iop_prol,ilaggr,nlaggr)
if (info == psb_success_) call iop_prol%free()
if (info == psb_success_) call iop_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999

@ -108,10 +108,12 @@ subroutine mld_c_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold)
character(len=20) :: name, ch_err
logical, parameter :: debug=.false.
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
@ -124,6 +126,12 @@ subroutine mld_c_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Entering '
#if defined(LPK8)
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Need fix for LPK8')
goto 9999
#else
!
! For the time being we are commenting out the UPDATE argument
! we plan to resurrect it later.
@ -321,6 +329,7 @@ subroutine mld_c_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Exiting with',iszv,' levels'
#endif
call psb_erractionrestore(err_act)
return
@ -346,7 +355,7 @@ contains
! Local variables
character(len=20) :: name
integer(psb_mpik_) :: ictxt, np, me, ncol
integer(psb_mpk_) :: ictxt, np, me, ncol
integer(psb_ipk_) :: err_act,ntaggr,nzl
integer(psb_ipk_), allocatable :: ilaggr(:), nlaggr(:)
type(psb_cspmat_type) :: ac, am2, am3, am4
@ -355,11 +364,18 @@ contains
logical, parameter :: debug=.false.
name='mld_c_extaggr_bld'
if (psb_get_errstatus().ne.0) return
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
info = psb_success_
ictxt = desc_a%get_context()
call psb_info(ictxt,me,np)
#if defined(LPK8)
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Need fix for LPK8')
goto 9999
#else
allocate(nlaggr(np),ilaggr(1))
nlaggr = 0
ilaggr = 0
@ -506,7 +522,7 @@ contains
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999
end if
#endif
call psb_erractionrestore(err_act)
return

@ -86,17 +86,19 @@ subroutine mld_c_hierarchy_bld(a,desc_a,prec,info)
& base_sm2, med_sm2, coarse_sm2
class(mld_c_base_aggregator_type), allocatable :: tmp_aggr
type(mld_sml_parms) :: baseparms, medparms, coarseparms
integer(psb_ipk_), allocatable :: ilaggr(:), nlaggr(:)
type(psb_cspmat_type) :: op_prol
integer(psb_lpk_), allocatable :: ilaggr(:), nlaggr(:)
type(psb_lcspmat_type) :: op_prol
type(mld_c_onelev_type), allocatable :: tprecv(:)
integer(psb_ipk_) :: int_err(5)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()

@ -104,10 +104,12 @@ subroutine mld_c_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()

@ -155,7 +155,10 @@
! local part of the preconditioner to be applied.
! Note that nlev = size(p%precv) = number of levels.
! p%precv(lev)%sm - type(psb_cbaseprec_type)
! The 'smoother' for the current level
! The pre-'smoother' for the current level
! p%precv(lev)%sm2 - type(psb_cbaseprec_type)
! The post-'smoother' for the current level
! may be the same or different from %sm
! p%precv(lev)%ac - type(psb_cspmat_type)
! The local part of the matrix A(lev).
! p%precv(lev)%parms - type(psb_sml_parms)

@ -102,10 +102,12 @@ subroutine mld_cmlprec_bld(a,desc_a,p,info,amold,vmold,imold)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()

@ -416,7 +416,7 @@ subroutine mld_cprecaply2_vect(prec,x,y,desc_data,info,trans,work)
& nswps,work_,wv,info)
end if
end associate
if (psb_get_errstatus() /=0) info = psb_err_internal_error_
if (psb_errstatus_fatal()) info = psb_err_internal_error_
if (info /= 0) then
info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='Smoother application',&
@ -563,7 +563,7 @@ subroutine mld_cprecaply1_vect(prec,x,desc_data,info,trans,work)
if (info == 0) call psb_geaxpby(cone,ww,czero,x,desc_data,info)
end if
if (psb_get_errstatus() /=0) info = psb_err_internal_error_
if (psb_errstatus_fatal()) info = psb_err_internal_error_
if (info /=0) then
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='Smoother application',&

@ -83,10 +83,12 @@ subroutine mld_cprecbld(a,desc_a,prec,info,amold,vmold,imold)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()

@ -108,10 +108,12 @@ subroutine mld_d_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold)
character(len=20) :: name, ch_err
logical, parameter :: debug=.false.
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
@ -124,6 +126,12 @@ subroutine mld_d_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Entering '
#if defined(LPK8)
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Need fix for LPK8')
goto 9999
#else
!
! For the time being we are commenting out the UPDATE argument
! we plan to resurrect it later.
@ -321,6 +329,7 @@ subroutine mld_d_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Exiting with',iszv,' levels'
#endif
call psb_erractionrestore(err_act)
return
@ -346,7 +355,7 @@ contains
! Local variables
character(len=20) :: name
integer(psb_mpik_) :: ictxt, np, me, ncol
integer(psb_mpk_) :: ictxt, np, me, ncol
integer(psb_ipk_) :: err_act,ntaggr,nzl
integer(psb_ipk_), allocatable :: ilaggr(:), nlaggr(:)
type(psb_dspmat_type) :: ac, am2, am3, am4
@ -355,11 +364,18 @@ contains
logical, parameter :: debug=.false.
name='mld_d_extaggr_bld'
if (psb_get_errstatus().ne.0) return
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
info = psb_success_
ictxt = desc_a%get_context()
call psb_info(ictxt,me,np)
#if defined(LPK8)
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Need fix for LPK8')
goto 9999
#else
allocate(nlaggr(np),ilaggr(1))
nlaggr = 0
ilaggr = 0
@ -506,7 +522,7 @@ contains
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999
end if
#endif
call psb_erractionrestore(err_act)
return

@ -86,17 +86,19 @@ subroutine mld_d_hierarchy_bld(a,desc_a,prec,info)
& base_sm2, med_sm2, coarse_sm2
class(mld_d_base_aggregator_type), allocatable :: tmp_aggr
type(mld_dml_parms) :: baseparms, medparms, coarseparms
integer(psb_ipk_), allocatable :: ilaggr(:), nlaggr(:)
type(psb_dspmat_type) :: op_prol
integer(psb_lpk_), allocatable :: ilaggr(:), nlaggr(:)
type(psb_ldspmat_type) :: op_prol
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
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()

@ -104,10 +104,12 @@ subroutine mld_d_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()

@ -155,7 +155,10 @@
! local part of the preconditioner to be applied.
! Note that nlev = size(p%precv) = number of levels.
! p%precv(lev)%sm - type(psb_dbaseprec_type)
! The 'smoother' for the current level
! The pre-'smoother' for the current level
! p%precv(lev)%sm2 - type(psb_dbaseprec_type)
! The post-'smoother' for the current level
! may be the same or different from %sm
! p%precv(lev)%ac - type(psb_dspmat_type)
! The local part of the matrix A(lev).
! p%precv(lev)%parms - type(psb_dml_parms)

@ -102,10 +102,12 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info,amold,vmold,imold)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()

@ -416,7 +416,7 @@ subroutine mld_dprecaply2_vect(prec,x,y,desc_data,info,trans,work)
& nswps,work_,wv,info)
end if
end associate
if (psb_get_errstatus() /=0) info = psb_err_internal_error_
if (psb_errstatus_fatal()) info = psb_err_internal_error_
if (info /= 0) then
info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='Smoother application',&
@ -563,7 +563,7 @@ subroutine mld_dprecaply1_vect(prec,x,desc_data,info,trans,work)
if (info == 0) call psb_geaxpby(done,ww,dzero,x,desc_data,info)
end if
if (psb_get_errstatus() /=0) info = psb_err_internal_error_
if (psb_errstatus_fatal()) info = psb_err_internal_error_
if (info /=0) then
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='Smoother application',&

@ -83,10 +83,12 @@ subroutine mld_dprecbld(a,desc_a,prec,info,amold,vmold,imold)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()

@ -108,10 +108,12 @@ subroutine mld_s_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold)
character(len=20) :: name, ch_err
logical, parameter :: debug=.false.
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
@ -124,6 +126,12 @@ subroutine mld_s_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Entering '
#if defined(LPK8)
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Need fix for LPK8')
goto 9999
#else
!
! For the time being we are commenting out the UPDATE argument
! we plan to resurrect it later.
@ -321,6 +329,7 @@ subroutine mld_s_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Exiting with',iszv,' levels'
#endif
call psb_erractionrestore(err_act)
return
@ -346,7 +355,7 @@ contains
! Local variables
character(len=20) :: name
integer(psb_mpik_) :: ictxt, np, me, ncol
integer(psb_mpk_) :: ictxt, np, me, ncol
integer(psb_ipk_) :: err_act,ntaggr,nzl
integer(psb_ipk_), allocatable :: ilaggr(:), nlaggr(:)
type(psb_sspmat_type) :: ac, am2, am3, am4
@ -355,11 +364,18 @@ contains
logical, parameter :: debug=.false.
name='mld_s_extaggr_bld'
if (psb_get_errstatus().ne.0) return
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
info = psb_success_
ictxt = desc_a%get_context()
call psb_info(ictxt,me,np)
#if defined(LPK8)
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Need fix for LPK8')
goto 9999
#else
allocate(nlaggr(np),ilaggr(1))
nlaggr = 0
ilaggr = 0
@ -506,7 +522,7 @@ contains
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999
end if
#endif
call psb_erractionrestore(err_act)
return

@ -86,17 +86,19 @@ subroutine mld_s_hierarchy_bld(a,desc_a,prec,info)
& base_sm2, med_sm2, coarse_sm2
class(mld_s_base_aggregator_type), allocatable :: tmp_aggr
type(mld_sml_parms) :: baseparms, medparms, coarseparms
integer(psb_ipk_), allocatable :: ilaggr(:), nlaggr(:)
type(psb_sspmat_type) :: op_prol
integer(psb_lpk_), allocatable :: ilaggr(:), nlaggr(:)
type(psb_lsspmat_type) :: op_prol
type(mld_s_onelev_type), allocatable :: tprecv(:)
integer(psb_ipk_) :: int_err(5)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()

@ -104,10 +104,12 @@ subroutine mld_s_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()

@ -155,7 +155,10 @@
! local part of the preconditioner to be applied.
! Note that nlev = size(p%precv) = number of levels.
! p%precv(lev)%sm - type(psb_sbaseprec_type)
! The 'smoother' for the current level
! The pre-'smoother' for the current level
! p%precv(lev)%sm2 - type(psb_sbaseprec_type)
! The post-'smoother' for the current level
! may be the same or different from %sm
! p%precv(lev)%ac - type(psb_sspmat_type)
! The local part of the matrix A(lev).
! p%precv(lev)%parms - type(psb_sml_parms)

@ -102,10 +102,12 @@ subroutine mld_smlprec_bld(a,desc_a,p,info,amold,vmold,imold)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()

@ -416,7 +416,7 @@ subroutine mld_sprecaply2_vect(prec,x,y,desc_data,info,trans,work)
& nswps,work_,wv,info)
end if
end associate
if (psb_get_errstatus() /=0) info = psb_err_internal_error_
if (psb_errstatus_fatal()) info = psb_err_internal_error_
if (info /= 0) then
info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='Smoother application',&
@ -563,7 +563,7 @@ subroutine mld_sprecaply1_vect(prec,x,desc_data,info,trans,work)
if (info == 0) call psb_geaxpby(sone,ww,szero,x,desc_data,info)
end if
if (psb_get_errstatus() /=0) info = psb_err_internal_error_
if (psb_errstatus_fatal()) info = psb_err_internal_error_
if (info /=0) then
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='Smoother application',&

@ -83,10 +83,12 @@ subroutine mld_sprecbld(a,desc_a,prec,info,amold,vmold,imold)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()

@ -108,10 +108,12 @@ subroutine mld_z_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold)
character(len=20) :: name, ch_err
logical, parameter :: debug=.false.
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
@ -124,6 +126,12 @@ subroutine mld_z_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Entering '
#if defined(LPK8)
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Need fix for LPK8')
goto 9999
#else
!
! For the time being we are commenting out the UPDATE argument
! we plan to resurrect it later.
@ -321,6 +329,7 @@ subroutine mld_z_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Exiting with',iszv,' levels'
#endif
call psb_erractionrestore(err_act)
return
@ -346,7 +355,7 @@ contains
! Local variables
character(len=20) :: name
integer(psb_mpik_) :: ictxt, np, me, ncol
integer(psb_mpk_) :: ictxt, np, me, ncol
integer(psb_ipk_) :: err_act,ntaggr,nzl
integer(psb_ipk_), allocatable :: ilaggr(:), nlaggr(:)
type(psb_zspmat_type) :: ac, am2, am3, am4
@ -355,11 +364,18 @@ contains
logical, parameter :: debug=.false.
name='mld_z_extaggr_bld'
if (psb_get_errstatus().ne.0) return
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
info = psb_success_
ictxt = desc_a%get_context()
call psb_info(ictxt,me,np)
#if defined(LPK8)
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='Need fix for LPK8')
goto 9999
#else
allocate(nlaggr(np),ilaggr(1))
nlaggr = 0
ilaggr = 0
@ -506,7 +522,7 @@ contains
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999
end if
#endif
call psb_erractionrestore(err_act)
return

@ -86,17 +86,19 @@ subroutine mld_z_hierarchy_bld(a,desc_a,prec,info)
& base_sm2, med_sm2, coarse_sm2
class(mld_z_base_aggregator_type), allocatable :: tmp_aggr
type(mld_dml_parms) :: baseparms, medparms, coarseparms
integer(psb_ipk_), allocatable :: ilaggr(:), nlaggr(:)
type(psb_zspmat_type) :: op_prol
integer(psb_lpk_), allocatable :: ilaggr(:), nlaggr(:)
type(psb_lzspmat_type) :: op_prol
type(mld_z_onelev_type), allocatable :: tprecv(:)
integer(psb_ipk_) :: int_err(5)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()

@ -104,10 +104,12 @@ subroutine mld_z_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()

@ -155,7 +155,10 @@
! local part of the preconditioner to be applied.
! Note that nlev = size(p%precv) = number of levels.
! p%precv(lev)%sm - type(psb_zbaseprec_type)
! The 'smoother' for the current level
! The pre-'smoother' for the current level
! p%precv(lev)%sm2 - type(psb_zbaseprec_type)
! The post-'smoother' for the current level
! may be the same or different from %sm
! p%precv(lev)%ac - type(psb_zspmat_type)
! The local part of the matrix A(lev).
! p%precv(lev)%parms - type(psb_dml_parms)

@ -102,10 +102,12 @@ subroutine mld_zmlprec_bld(a,desc_a,p,info,amold,vmold,imold)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()

@ -416,7 +416,7 @@ subroutine mld_zprecaply2_vect(prec,x,y,desc_data,info,trans,work)
& nswps,work_,wv,info)
end if
end associate
if (psb_get_errstatus() /=0) info = psb_err_internal_error_
if (psb_errstatus_fatal()) info = psb_err_internal_error_
if (info /= 0) then
info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='Smoother application',&
@ -563,7 +563,7 @@ subroutine mld_zprecaply1_vect(prec,x,desc_data,info,trans,work)
if (info == 0) call psb_geaxpby(zone,ww,zzero,x,desc_data,info)
end if
if (psb_get_errstatus() /=0) info = psb_err_internal_error_
if (psb_errstatus_fatal()) info = psb_err_internal_error_
if (info /=0) then
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='Smoother application',&

@ -83,10 +83,12 @@ subroutine mld_zprecbld(a,desc_a,prec,info,amold,vmold,imold)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
if (psb_get_errstatus().ne.0) return
info=psb_success_
err=0
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()

@ -78,7 +78,7 @@ subroutine mld_c_bwgs_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
! This cuts out the off-diagonal part, because it's supposed to
! be handled by the outer Jacobi smoother.
!
call a%tril(sv%l,info,diag=-1,jmax=nrow_a,u=sv%u)
call a%tril(sv%l,info,diag=-ione,jmax=nrow_a,u=sv%u)
else

@ -76,7 +76,7 @@ subroutine mld_c_gs_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
! This cuts out the off-diagonal part, because it's supposed to
! be handled by the outer Jacobi smoother.
!
call a%tril(sv%l,info,diag=0,jmax=nrow_a,u=sv%u)
call a%tril(sv%l,info,diag=izero,jmax=nrow_a,u=sv%u)
else

@ -64,7 +64,7 @@ subroutine c_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
call psb_erractionsave(err_act)
#if defined(HAVE_MUMPS_)
#if defined(HAVE_MUMPS_) && !defined(LPK8)
info = psb_success_
trans_ = psb_toupper(trans)
select case(trans_)

@ -59,7 +59,7 @@ subroutine c_mumps_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
integer(psb_ipk_) :: err_act
character(len=20) :: name='c_mumps_solver_apply_vect'
#if defined(HAVE_MUMPS_)
#if defined(HAVE_MUMPS_) && !defined(LPK8)
call psb_erractionsave(err_act)

@ -63,7 +63,7 @@
integer(psb_ipk_) :: ictxt, ictxt1, icomm, np, iam, me, i, err_act, debug_unit, debug_level
character(len=20) :: name='c_mumps_solver_bld', ch_err
#if defined(HAVE_MUMPS_)
#if defined(HAVE_MUMPS_) && !defined(LPK8)
info=psb_success_

@ -304,9 +304,11 @@ contains
character(len=20) :: name, ch_err
name='mld_cilu0_factint'
if(psb_get_errstatus().ne.0) return
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
ma = a%get_nrows()
mb = b%get_nrows()
@ -569,9 +571,11 @@ contains
character(len=20), parameter :: name='ilu_copyin'
character(len=20) :: ch_err
if (psb_get_errstatus() /= 0) return
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
if (psb_toupper(upd) == 'F') then
select type(aa => a%a)

@ -289,9 +289,11 @@ contains
character(len=20), parameter :: name='mld_ciluk_factint'
character(len=20) :: ch_err
if (psb_get_errstatus() /= 0) return
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
select case(ialg)
@ -503,9 +505,11 @@ contains
character(len=20), parameter :: name='iluk_copyin'
character(len=20) :: ch_err
if (psb_get_errstatus() /= 0) return
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
call heap%init(info)
select type (aa=> a%a)
@ -829,9 +833,11 @@ contains
character(len=20), parameter :: name='mld_ciluk_factint'
character(len=20) :: ch_err
if (psb_get_errstatus() /= 0) return
info = psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
d(i) = czero

@ -311,9 +311,11 @@ contains
character(len=20), parameter :: name='mld_cilut_factint'
character(len=20) :: ch_err
if (psb_get_errstatus() /= 0) return
info = psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
ma = a%get_nrows()
@ -524,9 +526,11 @@ contains
real(psb_spk_), external :: dnrm2
character(len=20), parameter :: name='mld_cilut_factint'
if (psb_get_errstatus() /= 0) return
info = psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
call heap%init(info)
if (info /= psb_success_) then
@ -907,9 +911,11 @@ contains
character(len=20) :: ch_err
logical :: fndmaxup
if (psb_get_errstatus() /= 0) return
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_; goto 9999
end if
!
! Here we need to apply also the dropping rule base on the fill-in.

@ -78,7 +78,7 @@ subroutine mld_d_bwgs_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
! This cuts out the off-diagonal part, because it's supposed to
! be handled by the outer Jacobi smoother.
!
call a%tril(sv%l,info,diag=-1,jmax=nrow_a,u=sv%u)
call a%tril(sv%l,info,diag=-ione,jmax=nrow_a,u=sv%u)
else

Some files were not shown because too many files have changed in this diff Show More

Loading…
Cancel
Save