From 5a1dace4e32c24eb30099a405ea8d767e2e3862b Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Thu, 21 Mar 2019 18:08:14 +0000 Subject: [PATCH] Fixed singleton marking in SOC1/SOC2. Fixed singleton handling for unsmoothed aggregation and minenergy. --- mlprec/impl/aggregator/mld_c_map_to_tprol.f90 | 17 +++++++--- mlprec/impl/aggregator/mld_c_soc1_map_bld.f90 | 31 ++++++++++++------- mlprec/impl/aggregator/mld_c_soc2_map_bld.f90 | 31 +++++++++++++------ .../aggregator/mld_caggrmat_minnrg_asb.f90 | 12 +++++-- .../aggregator/mld_caggrmat_nosmth_asb.f90 | 4 +++ mlprec/impl/aggregator/mld_d_map_to_tprol.f90 | 17 +++++++--- mlprec/impl/aggregator/mld_d_soc1_map_bld.f90 | 31 ++++++++++++------- mlprec/impl/aggregator/mld_d_soc2_map_bld.f90 | 31 +++++++++++++------ .../aggregator/mld_daggrmat_minnrg_asb.f90 | 12 +++++-- .../aggregator/mld_daggrmat_nosmth_asb.f90 | 4 +++ mlprec/impl/aggregator/mld_s_map_to_tprol.f90 | 17 +++++++--- mlprec/impl/aggregator/mld_s_soc1_map_bld.f90 | 31 ++++++++++++------- mlprec/impl/aggregator/mld_s_soc2_map_bld.f90 | 31 +++++++++++++------ .../aggregator/mld_saggrmat_minnrg_asb.f90 | 12 +++++-- .../aggregator/mld_saggrmat_nosmth_asb.f90 | 4 +++ mlprec/impl/aggregator/mld_z_map_to_tprol.f90 | 17 +++++++--- mlprec/impl/aggregator/mld_z_soc1_map_bld.f90 | 31 ++++++++++++------- mlprec/impl/aggregator/mld_z_soc2_map_bld.f90 | 31 +++++++++++++------ .../aggregator/mld_zaggrmat_minnrg_asb.f90 | 12 +++++-- .../aggregator/mld_zaggrmat_nosmth_asb.f90 | 4 +++ 20 files changed, 276 insertions(+), 104 deletions(-) diff --git a/mlprec/impl/aggregator/mld_c_map_to_tprol.f90 b/mlprec/impl/aggregator/mld_c_map_to_tprol.f90 index e7464638..1c2c5516 100644 --- a/mlprec/impl/aggregator/mld_c_map_to_tprol.f90 +++ b/mlprec/impl/aggregator/mld_c_map_to_tprol.f90 @@ -125,12 +125,21 @@ subroutine mld_c_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) end if call tmpcoo%allocate(ncol,ntaggr,ncol) + k = 0 do i=1,ncol - tmpcoo%val(i) = cone - tmpcoo%ia(i) = i - tmpcoo%ja(i) = ilaggr(i) + ! + ! Note: at this point, a value ilaggr(i)<=0 + ! tags a "singleton" row, and it has to be + ! left alone. + ! + if (ilaggr(i)>0) then + k = k + 1 + tmpcoo%val(k) = cone + tmpcoo%ia(k) = i + tmpcoo%ja(k) = ilaggr(i) + end if end do - call tmpcoo%set_nzeros(ncol) + call tmpcoo%set_nzeros(k) call tmpcoo%set_dupl(psb_dupl_add_) call tmpcoo%set_sorted() ! At this point this is in row-major call op_prol%mv_from(tmpcoo) diff --git a/mlprec/impl/aggregator/mld_c_soc1_map_bld.f90 b/mlprec/impl/aggregator/mld_c_soc1_map_bld.f90 index 23a13ceb..53bebfe0 100644 --- a/mlprec/impl/aggregator/mld_c_soc1_map_bld.f90 +++ b/mlprec/impl/aggregator/mld_c_soc1_map_bld.f90 @@ -95,6 +95,7 @@ subroutine mld_c_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) 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_) :: nrglob character(len=20) :: name, ch_err info=psb_success_ @@ -110,6 +111,7 @@ subroutine mld_c_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) call psb_info(ictxt,me,np) nrow = desc_a%get_local_rows() ncol = desc_a%get_local_cols() + nrglob = desc_a%get_global_rows() nr = a%get_nrows() nc = a%get_ncols() @@ -261,12 +263,6 @@ subroutine mld_c_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) nz = (acsr%irp(i+1)-acsr%irp(i)) icol(1:nz) = acsr%ja(acsr%irp(i):acsr%irp(i+1)-1) val(1:nz) = acsr%val(acsr%irp(i):acsr%irp(i+1)-1) -!!$ call a%csget(i,i,nz,irow,icol,val,info,chksz=.false.) -!!$ if (info /= psb_success_) then -!!$ info=psb_err_from_subroutine_ -!!$ call psb_errpush(info,name,a_err='psb_sp_getrow') -!!$ goto 9999 -!!$ end if ! ! Find its strongly connected neighbourhood not ! already aggregated, and make it into a new aggregate. @@ -302,11 +298,24 @@ subroutine mld_c_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) end do step3 - if (count(ilaggr<0) >0) then - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='Fatal error: some leftovers') - goto 9999 - endif + ! Any leftovers? + do i=1, nr + if (ilaggr(i) < 0) then + nz = (acsr%irp(i+1)-acsr%irp(i)) + if (nz == 1) then + ! Mark explicitly as a singleton so that + ! it will be ignored in map_to_tprol. + ! Need to use -(nrglob+nr) to make sure + ! it's still negative when shifted and combined with + ! other processes. + ilaggr(i) = -(nrglob+nr) + else + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: non-singleton leftovers') + goto 9999 + endif + end if + end do if (naggr > ncol) then !write(0,*) name,'Error : naggr > ncol',naggr,ncol diff --git a/mlprec/impl/aggregator/mld_c_soc2_map_bld.f90 b/mlprec/impl/aggregator/mld_c_soc2_map_bld.f90 index d3185a68..6b9bdc6e 100644 --- a/mlprec/impl/aggregator/mld_c_soc2_map_bld.f90 +++ b/mlprec/impl/aggregator/mld_c_soc2_map_bld.f90 @@ -88,6 +88,7 @@ subroutine mld_c_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) 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 + integer(psb_lpk_) :: nrglob type(psb_c_csr_sparse_mat) :: acsr, muij, s_neigh type(psb_c_coo_sparse_mat) :: s_neigh_coo real(psb_spk_) :: cpling, tcl @@ -108,8 +109,9 @@ subroutine mld_c_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! ictxt=desc_a%get_context() call psb_info(ictxt,me,np) - nrow = desc_a%get_local_rows() - ncol = desc_a%get_local_cols() + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + nrglob = desc_a%get_global_rows() nr = a%get_nrows() nc = a%get_ncols() @@ -288,15 +290,26 @@ subroutine mld_c_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) end if end do step3 - - if (count(ilaggr<0) >0) then - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='Fatal error: some leftovers') - goto 9999 - endif + ! Any leftovers? + do i=1, nr + if (ilaggr(i) <= 0) then + nz = (s_neigh%irp(i+1)-s_neigh%irp(i)) + if (nz <= 1) then + ! Mark explicitly as a singleton so that + ! it will be ignored in map_to_tprol. + ! Need to use -(nrglob+nr) to make sure + ! it's still negative when shifted and combined with + ! other processes. + ilaggr(i) = -(nrglob+nr) + else + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: non-singleton leftovers') + goto 9999 + endif + end if + end do if (naggr > ncol) then - write(0,*) name,'Error : naggr > ncol',naggr,ncol info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Fatal error: naggr>ncol') goto 9999 diff --git a/mlprec/impl/aggregator/mld_caggrmat_minnrg_asb.f90 b/mlprec/impl/aggregator/mld_caggrmat_minnrg_asb.f90 index 09e8298b..12eb9379 100644 --- a/mlprec/impl/aggregator/mld_caggrmat_minnrg_asb.f90 +++ b/mlprec/impl/aggregator/mld_caggrmat_minnrg_asb.f90 @@ -260,7 +260,11 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re ! Compute omega_int ommx = czero do i=1, ncol - omi(i) = omp(ilaggr(i)) + if (ilaggr(i) >0) then + omi(i) = omp(ilaggr(i)) + else + omi(i) = czero + end if if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) end do ! Compute omega_fine @@ -409,7 +413,11 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re ! Compute omega_int ommx = czero do i=1, ncol - omi(i) = omp(ilaggr(i)) + if (ilaggr(i) >0) then + omi(i) = omp(ilaggr(i)) + else + omi(i) = czero + end if if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) end do ! Compute omega_fine diff --git a/mlprec/impl/aggregator/mld_caggrmat_nosmth_asb.f90 b/mlprec/impl/aggregator/mld_caggrmat_nosmth_asb.f90 index 4224e209..7eb18dcc 100644 --- a/mlprec/impl/aggregator/mld_caggrmat_nosmth_asb.f90 +++ b/mlprec/impl/aggregator/mld_caggrmat_nosmth_asb.f90 @@ -178,6 +178,10 @@ subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re ac_coo%ia(k) = ilaggr(ac_coo%ia(i)) ac_coo%ja(k) = ilaggr(ac_coo%ja(i)) ac_coo%val(k) = ac_coo%val(i) + ! At this point, there may be negative entries, + ! because that's how ILAGGR marks singletons + ! If this is the case, roll back K + if ((ac_coo%ia(k)<=0).or.(ac_coo%ja(k)<=0)) k = k-1 enddo call ac_coo%set_nrows(naggr) call ac_coo%set_ncols(naggr) diff --git a/mlprec/impl/aggregator/mld_d_map_to_tprol.f90 b/mlprec/impl/aggregator/mld_d_map_to_tprol.f90 index dd8b1447..b8117d94 100644 --- a/mlprec/impl/aggregator/mld_d_map_to_tprol.f90 +++ b/mlprec/impl/aggregator/mld_d_map_to_tprol.f90 @@ -125,12 +125,21 @@ subroutine mld_d_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) end if call tmpcoo%allocate(ncol,ntaggr,ncol) + k = 0 do i=1,ncol - tmpcoo%val(i) = done - tmpcoo%ia(i) = i - tmpcoo%ja(i) = ilaggr(i) + ! + ! Note: at this point, a value ilaggr(i)<=0 + ! tags a "singleton" row, and it has to be + ! left alone. + ! + if (ilaggr(i)>0) then + k = k + 1 + tmpcoo%val(k) = done + tmpcoo%ia(k) = i + tmpcoo%ja(k) = ilaggr(i) + end if end do - call tmpcoo%set_nzeros(ncol) + call tmpcoo%set_nzeros(k) call tmpcoo%set_dupl(psb_dupl_add_) call tmpcoo%set_sorted() ! At this point this is in row-major call op_prol%mv_from(tmpcoo) diff --git a/mlprec/impl/aggregator/mld_d_soc1_map_bld.f90 b/mlprec/impl/aggregator/mld_d_soc1_map_bld.f90 index ab0d570d..24b83650 100644 --- a/mlprec/impl/aggregator/mld_d_soc1_map_bld.f90 +++ b/mlprec/impl/aggregator/mld_d_soc1_map_bld.f90 @@ -95,6 +95,7 @@ subroutine mld_d_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) 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_) :: nrglob character(len=20) :: name, ch_err info=psb_success_ @@ -110,6 +111,7 @@ subroutine mld_d_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) call psb_info(ictxt,me,np) nrow = desc_a%get_local_rows() ncol = desc_a%get_local_cols() + nrglob = desc_a%get_global_rows() nr = a%get_nrows() nc = a%get_ncols() @@ -261,12 +263,6 @@ subroutine mld_d_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) nz = (acsr%irp(i+1)-acsr%irp(i)) icol(1:nz) = acsr%ja(acsr%irp(i):acsr%irp(i+1)-1) val(1:nz) = acsr%val(acsr%irp(i):acsr%irp(i+1)-1) -!!$ call a%csget(i,i,nz,irow,icol,val,info,chksz=.false.) -!!$ if (info /= psb_success_) then -!!$ info=psb_err_from_subroutine_ -!!$ call psb_errpush(info,name,a_err='psb_sp_getrow') -!!$ goto 9999 -!!$ end if ! ! Find its strongly connected neighbourhood not ! already aggregated, and make it into a new aggregate. @@ -302,11 +298,24 @@ subroutine mld_d_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) end do step3 - if (count(ilaggr<0) >0) then - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='Fatal error: some leftovers') - goto 9999 - endif + ! Any leftovers? + do i=1, nr + if (ilaggr(i) < 0) then + nz = (acsr%irp(i+1)-acsr%irp(i)) + if (nz == 1) then + ! Mark explicitly as a singleton so that + ! it will be ignored in map_to_tprol. + ! Need to use -(nrglob+nr) to make sure + ! it's still negative when shifted and combined with + ! other processes. + ilaggr(i) = -(nrglob+nr) + else + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: non-singleton leftovers') + goto 9999 + endif + end if + end do if (naggr > ncol) then !write(0,*) name,'Error : naggr > ncol',naggr,ncol diff --git a/mlprec/impl/aggregator/mld_d_soc2_map_bld.f90 b/mlprec/impl/aggregator/mld_d_soc2_map_bld.f90 index ae9c5346..0307ca73 100644 --- a/mlprec/impl/aggregator/mld_d_soc2_map_bld.f90 +++ b/mlprec/impl/aggregator/mld_d_soc2_map_bld.f90 @@ -88,6 +88,7 @@ subroutine mld_d_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) 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 + integer(psb_lpk_) :: nrglob type(psb_d_csr_sparse_mat) :: acsr, muij, s_neigh type(psb_d_coo_sparse_mat) :: s_neigh_coo real(psb_dpk_) :: cpling, tcl @@ -108,8 +109,9 @@ subroutine mld_d_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! ictxt=desc_a%get_context() call psb_info(ictxt,me,np) - nrow = desc_a%get_local_rows() - ncol = desc_a%get_local_cols() + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + nrglob = desc_a%get_global_rows() nr = a%get_nrows() nc = a%get_ncols() @@ -288,15 +290,26 @@ subroutine mld_d_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) end if end do step3 - - if (count(ilaggr<0) >0) then - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='Fatal error: some leftovers') - goto 9999 - endif + ! Any leftovers? + do i=1, nr + if (ilaggr(i) <= 0) then + nz = (s_neigh%irp(i+1)-s_neigh%irp(i)) + if (nz <= 1) then + ! Mark explicitly as a singleton so that + ! it will be ignored in map_to_tprol. + ! Need to use -(nrglob+nr) to make sure + ! it's still negative when shifted and combined with + ! other processes. + ilaggr(i) = -(nrglob+nr) + else + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: non-singleton leftovers') + goto 9999 + endif + end if + end do if (naggr > ncol) then - write(0,*) name,'Error : naggr > ncol',naggr,ncol info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Fatal error: naggr>ncol') goto 9999 diff --git a/mlprec/impl/aggregator/mld_daggrmat_minnrg_asb.f90 b/mlprec/impl/aggregator/mld_daggrmat_minnrg_asb.f90 index 1b33ae3b..e6bf12af 100644 --- a/mlprec/impl/aggregator/mld_daggrmat_minnrg_asb.f90 +++ b/mlprec/impl/aggregator/mld_daggrmat_minnrg_asb.f90 @@ -260,7 +260,11 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re ! Compute omega_int ommx = dzero do i=1, ncol - omi(i) = omp(ilaggr(i)) + if (ilaggr(i) >0) then + omi(i) = omp(ilaggr(i)) + else + omi(i) = dzero + end if if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) end do ! Compute omega_fine @@ -409,7 +413,11 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re ! Compute omega_int ommx = dzero do i=1, ncol - omi(i) = omp(ilaggr(i)) + if (ilaggr(i) >0) then + omi(i) = omp(ilaggr(i)) + else + omi(i) = dzero + end if if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) end do ! Compute omega_fine diff --git a/mlprec/impl/aggregator/mld_daggrmat_nosmth_asb.f90 b/mlprec/impl/aggregator/mld_daggrmat_nosmth_asb.f90 index 15bbb30a..15c8cb4b 100644 --- a/mlprec/impl/aggregator/mld_daggrmat_nosmth_asb.f90 +++ b/mlprec/impl/aggregator/mld_daggrmat_nosmth_asb.f90 @@ -178,6 +178,10 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re ac_coo%ia(k) = ilaggr(ac_coo%ia(i)) ac_coo%ja(k) = ilaggr(ac_coo%ja(i)) ac_coo%val(k) = ac_coo%val(i) + ! At this point, there may be negative entries, + ! because that's how ILAGGR marks singletons + ! If this is the case, roll back K + if ((ac_coo%ia(k)<=0).or.(ac_coo%ja(k)<=0)) k = k-1 enddo call ac_coo%set_nrows(naggr) call ac_coo%set_ncols(naggr) diff --git a/mlprec/impl/aggregator/mld_s_map_to_tprol.f90 b/mlprec/impl/aggregator/mld_s_map_to_tprol.f90 index ee1bc70e..d0d68654 100644 --- a/mlprec/impl/aggregator/mld_s_map_to_tprol.f90 +++ b/mlprec/impl/aggregator/mld_s_map_to_tprol.f90 @@ -125,12 +125,21 @@ subroutine mld_s_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) end if call tmpcoo%allocate(ncol,ntaggr,ncol) + k = 0 do i=1,ncol - tmpcoo%val(i) = sone - tmpcoo%ia(i) = i - tmpcoo%ja(i) = ilaggr(i) + ! + ! Note: at this point, a value ilaggr(i)<=0 + ! tags a "singleton" row, and it has to be + ! left alone. + ! + if (ilaggr(i)>0) then + k = k + 1 + tmpcoo%val(k) = sone + tmpcoo%ia(k) = i + tmpcoo%ja(k) = ilaggr(i) + end if end do - call tmpcoo%set_nzeros(ncol) + call tmpcoo%set_nzeros(k) call tmpcoo%set_dupl(psb_dupl_add_) call tmpcoo%set_sorted() ! At this point this is in row-major call op_prol%mv_from(tmpcoo) diff --git a/mlprec/impl/aggregator/mld_s_soc1_map_bld.f90 b/mlprec/impl/aggregator/mld_s_soc1_map_bld.f90 index ef0f0045..4c197a1a 100644 --- a/mlprec/impl/aggregator/mld_s_soc1_map_bld.f90 +++ b/mlprec/impl/aggregator/mld_s_soc1_map_bld.f90 @@ -95,6 +95,7 @@ subroutine mld_s_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) 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_) :: nrglob character(len=20) :: name, ch_err info=psb_success_ @@ -110,6 +111,7 @@ subroutine mld_s_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) call psb_info(ictxt,me,np) nrow = desc_a%get_local_rows() ncol = desc_a%get_local_cols() + nrglob = desc_a%get_global_rows() nr = a%get_nrows() nc = a%get_ncols() @@ -261,12 +263,6 @@ subroutine mld_s_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) nz = (acsr%irp(i+1)-acsr%irp(i)) icol(1:nz) = acsr%ja(acsr%irp(i):acsr%irp(i+1)-1) val(1:nz) = acsr%val(acsr%irp(i):acsr%irp(i+1)-1) -!!$ call a%csget(i,i,nz,irow,icol,val,info,chksz=.false.) -!!$ if (info /= psb_success_) then -!!$ info=psb_err_from_subroutine_ -!!$ call psb_errpush(info,name,a_err='psb_sp_getrow') -!!$ goto 9999 -!!$ end if ! ! Find its strongly connected neighbourhood not ! already aggregated, and make it into a new aggregate. @@ -302,11 +298,24 @@ subroutine mld_s_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) end do step3 - if (count(ilaggr<0) >0) then - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='Fatal error: some leftovers') - goto 9999 - endif + ! Any leftovers? + do i=1, nr + if (ilaggr(i) < 0) then + nz = (acsr%irp(i+1)-acsr%irp(i)) + if (nz == 1) then + ! Mark explicitly as a singleton so that + ! it will be ignored in map_to_tprol. + ! Need to use -(nrglob+nr) to make sure + ! it's still negative when shifted and combined with + ! other processes. + ilaggr(i) = -(nrglob+nr) + else + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: non-singleton leftovers') + goto 9999 + endif + end if + end do if (naggr > ncol) then !write(0,*) name,'Error : naggr > ncol',naggr,ncol diff --git a/mlprec/impl/aggregator/mld_s_soc2_map_bld.f90 b/mlprec/impl/aggregator/mld_s_soc2_map_bld.f90 index 682fb25c..770c51ac 100644 --- a/mlprec/impl/aggregator/mld_s_soc2_map_bld.f90 +++ b/mlprec/impl/aggregator/mld_s_soc2_map_bld.f90 @@ -88,6 +88,7 @@ subroutine mld_s_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) 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 + integer(psb_lpk_) :: nrglob type(psb_s_csr_sparse_mat) :: acsr, muij, s_neigh type(psb_s_coo_sparse_mat) :: s_neigh_coo real(psb_spk_) :: cpling, tcl @@ -108,8 +109,9 @@ subroutine mld_s_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! ictxt=desc_a%get_context() call psb_info(ictxt,me,np) - nrow = desc_a%get_local_rows() - ncol = desc_a%get_local_cols() + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + nrglob = desc_a%get_global_rows() nr = a%get_nrows() nc = a%get_ncols() @@ -288,15 +290,26 @@ subroutine mld_s_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) end if end do step3 - - if (count(ilaggr<0) >0) then - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='Fatal error: some leftovers') - goto 9999 - endif + ! Any leftovers? + do i=1, nr + if (ilaggr(i) <= 0) then + nz = (s_neigh%irp(i+1)-s_neigh%irp(i)) + if (nz <= 1) then + ! Mark explicitly as a singleton so that + ! it will be ignored in map_to_tprol. + ! Need to use -(nrglob+nr) to make sure + ! it's still negative when shifted and combined with + ! other processes. + ilaggr(i) = -(nrglob+nr) + else + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: non-singleton leftovers') + goto 9999 + endif + end if + end do if (naggr > ncol) then - write(0,*) name,'Error : naggr > ncol',naggr,ncol info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Fatal error: naggr>ncol') goto 9999 diff --git a/mlprec/impl/aggregator/mld_saggrmat_minnrg_asb.f90 b/mlprec/impl/aggregator/mld_saggrmat_minnrg_asb.f90 index 0f293e00..56da4362 100644 --- a/mlprec/impl/aggregator/mld_saggrmat_minnrg_asb.f90 +++ b/mlprec/impl/aggregator/mld_saggrmat_minnrg_asb.f90 @@ -260,7 +260,11 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re ! Compute omega_int ommx = szero do i=1, ncol - omi(i) = omp(ilaggr(i)) + if (ilaggr(i) >0) then + omi(i) = omp(ilaggr(i)) + else + omi(i) = szero + end if if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) end do ! Compute omega_fine @@ -409,7 +413,11 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re ! Compute omega_int ommx = szero do i=1, ncol - omi(i) = omp(ilaggr(i)) + if (ilaggr(i) >0) then + omi(i) = omp(ilaggr(i)) + else + omi(i) = szero + end if if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) end do ! Compute omega_fine diff --git a/mlprec/impl/aggregator/mld_saggrmat_nosmth_asb.f90 b/mlprec/impl/aggregator/mld_saggrmat_nosmth_asb.f90 index e588982c..9b8d99bf 100644 --- a/mlprec/impl/aggregator/mld_saggrmat_nosmth_asb.f90 +++ b/mlprec/impl/aggregator/mld_saggrmat_nosmth_asb.f90 @@ -178,6 +178,10 @@ subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re ac_coo%ia(k) = ilaggr(ac_coo%ia(i)) ac_coo%ja(k) = ilaggr(ac_coo%ja(i)) ac_coo%val(k) = ac_coo%val(i) + ! At this point, there may be negative entries, + ! because that's how ILAGGR marks singletons + ! If this is the case, roll back K + if ((ac_coo%ia(k)<=0).or.(ac_coo%ja(k)<=0)) k = k-1 enddo call ac_coo%set_nrows(naggr) call ac_coo%set_ncols(naggr) diff --git a/mlprec/impl/aggregator/mld_z_map_to_tprol.f90 b/mlprec/impl/aggregator/mld_z_map_to_tprol.f90 index 0a757ea1..b006b976 100644 --- a/mlprec/impl/aggregator/mld_z_map_to_tprol.f90 +++ b/mlprec/impl/aggregator/mld_z_map_to_tprol.f90 @@ -125,12 +125,21 @@ subroutine mld_z_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) end if call tmpcoo%allocate(ncol,ntaggr,ncol) + k = 0 do i=1,ncol - tmpcoo%val(i) = zone - tmpcoo%ia(i) = i - tmpcoo%ja(i) = ilaggr(i) + ! + ! Note: at this point, a value ilaggr(i)<=0 + ! tags a "singleton" row, and it has to be + ! left alone. + ! + if (ilaggr(i)>0) then + k = k + 1 + tmpcoo%val(k) = zone + tmpcoo%ia(k) = i + tmpcoo%ja(k) = ilaggr(i) + end if end do - call tmpcoo%set_nzeros(ncol) + call tmpcoo%set_nzeros(k) call tmpcoo%set_dupl(psb_dupl_add_) call tmpcoo%set_sorted() ! At this point this is in row-major call op_prol%mv_from(tmpcoo) diff --git a/mlprec/impl/aggregator/mld_z_soc1_map_bld.f90 b/mlprec/impl/aggregator/mld_z_soc1_map_bld.f90 index 4496683b..249a5ec0 100644 --- a/mlprec/impl/aggregator/mld_z_soc1_map_bld.f90 +++ b/mlprec/impl/aggregator/mld_z_soc1_map_bld.f90 @@ -95,6 +95,7 @@ subroutine mld_z_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) 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_) :: nrglob character(len=20) :: name, ch_err info=psb_success_ @@ -110,6 +111,7 @@ subroutine mld_z_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) call psb_info(ictxt,me,np) nrow = desc_a%get_local_rows() ncol = desc_a%get_local_cols() + nrglob = desc_a%get_global_rows() nr = a%get_nrows() nc = a%get_ncols() @@ -261,12 +263,6 @@ subroutine mld_z_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) nz = (acsr%irp(i+1)-acsr%irp(i)) icol(1:nz) = acsr%ja(acsr%irp(i):acsr%irp(i+1)-1) val(1:nz) = acsr%val(acsr%irp(i):acsr%irp(i+1)-1) -!!$ call a%csget(i,i,nz,irow,icol,val,info,chksz=.false.) -!!$ if (info /= psb_success_) then -!!$ info=psb_err_from_subroutine_ -!!$ call psb_errpush(info,name,a_err='psb_sp_getrow') -!!$ goto 9999 -!!$ end if ! ! Find its strongly connected neighbourhood not ! already aggregated, and make it into a new aggregate. @@ -302,11 +298,24 @@ subroutine mld_z_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) end do step3 - if (count(ilaggr<0) >0) then - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='Fatal error: some leftovers') - goto 9999 - endif + ! Any leftovers? + do i=1, nr + if (ilaggr(i) < 0) then + nz = (acsr%irp(i+1)-acsr%irp(i)) + if (nz == 1) then + ! Mark explicitly as a singleton so that + ! it will be ignored in map_to_tprol. + ! Need to use -(nrglob+nr) to make sure + ! it's still negative when shifted and combined with + ! other processes. + ilaggr(i) = -(nrglob+nr) + else + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: non-singleton leftovers') + goto 9999 + endif + end if + end do if (naggr > ncol) then !write(0,*) name,'Error : naggr > ncol',naggr,ncol diff --git a/mlprec/impl/aggregator/mld_z_soc2_map_bld.f90 b/mlprec/impl/aggregator/mld_z_soc2_map_bld.f90 index d7a78662..967ba393 100644 --- a/mlprec/impl/aggregator/mld_z_soc2_map_bld.f90 +++ b/mlprec/impl/aggregator/mld_z_soc2_map_bld.f90 @@ -88,6 +88,7 @@ subroutine mld_z_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) 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 + integer(psb_lpk_) :: nrglob type(psb_z_csr_sparse_mat) :: acsr, muij, s_neigh type(psb_z_coo_sparse_mat) :: s_neigh_coo real(psb_dpk_) :: cpling, tcl @@ -108,8 +109,9 @@ subroutine mld_z_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! ictxt=desc_a%get_context() call psb_info(ictxt,me,np) - nrow = desc_a%get_local_rows() - ncol = desc_a%get_local_cols() + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + nrglob = desc_a%get_global_rows() nr = a%get_nrows() nc = a%get_ncols() @@ -288,15 +290,26 @@ subroutine mld_z_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) end if end do step3 - - if (count(ilaggr<0) >0) then - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='Fatal error: some leftovers') - goto 9999 - endif + ! Any leftovers? + do i=1, nr + if (ilaggr(i) <= 0) then + nz = (s_neigh%irp(i+1)-s_neigh%irp(i)) + if (nz <= 1) then + ! Mark explicitly as a singleton so that + ! it will be ignored in map_to_tprol. + ! Need to use -(nrglob+nr) to make sure + ! it's still negative when shifted and combined with + ! other processes. + ilaggr(i) = -(nrglob+nr) + else + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: non-singleton leftovers') + goto 9999 + endif + end if + end do if (naggr > ncol) then - write(0,*) name,'Error : naggr > ncol',naggr,ncol info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Fatal error: naggr>ncol') goto 9999 diff --git a/mlprec/impl/aggregator/mld_zaggrmat_minnrg_asb.f90 b/mlprec/impl/aggregator/mld_zaggrmat_minnrg_asb.f90 index 8f121c6a..043855aa 100644 --- a/mlprec/impl/aggregator/mld_zaggrmat_minnrg_asb.f90 +++ b/mlprec/impl/aggregator/mld_zaggrmat_minnrg_asb.f90 @@ -260,7 +260,11 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re ! Compute omega_int ommx = zzero do i=1, ncol - omi(i) = omp(ilaggr(i)) + if (ilaggr(i) >0) then + omi(i) = omp(ilaggr(i)) + else + omi(i) = zzero + end if if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) end do ! Compute omega_fine @@ -409,7 +413,11 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re ! Compute omega_int ommx = zzero do i=1, ncol - omi(i) = omp(ilaggr(i)) + if (ilaggr(i) >0) then + omi(i) = omp(ilaggr(i)) + else + omi(i) = zzero + end if if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) end do ! Compute omega_fine diff --git a/mlprec/impl/aggregator/mld_zaggrmat_nosmth_asb.f90 b/mlprec/impl/aggregator/mld_zaggrmat_nosmth_asb.f90 index 30a8502b..152b2676 100644 --- a/mlprec/impl/aggregator/mld_zaggrmat_nosmth_asb.f90 +++ b/mlprec/impl/aggregator/mld_zaggrmat_nosmth_asb.f90 @@ -178,6 +178,10 @@ subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re ac_coo%ia(k) = ilaggr(ac_coo%ia(i)) ac_coo%ja(k) = ilaggr(ac_coo%ja(i)) ac_coo%val(k) = ac_coo%val(i) + ! At this point, there may be negative entries, + ! because that's how ILAGGR marks singletons + ! If this is the case, roll back K + if ((ac_coo%ia(k)<=0).or.(ac_coo%ja(k)<=0)) k = k-1 enddo call ac_coo%set_nrows(naggr) call ac_coo%set_ncols(naggr)