First hardcoded implementation of l1 smooth aggregation

amgprec/impl/aggregator/amg_daggrmat_smth_bld.f90

@@ -112,11 +112,11 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   implicit none
 
   ! Arguments
-  type(psb_dspmat_type), intent(in)      :: a
+  type(psb_dspmat_type), intent(in)        :: a
   type(psb_desc_type), intent(inout)       :: desc_a
   integer(psb_lpk_), intent(inout)         :: ilaggr(:), nlaggr(:)
-  type(amg_dml_parms), intent(inout)      :: parms 
+  type(amg_dml_parms), intent(inout)       :: parms 
-  type(psb_dspmat_type), intent(inout)      :: op_prol,ac,op_restr
+  type(psb_dspmat_type), intent(inout)     :: op_prol,ac,op_restr
   type(psb_ldspmat_type), intent(inout)    :: t_prol
   type(psb_desc_type), intent(inout)       :: desc_ac
   integer(psb_ipk_), intent(out)           :: info
@@ -132,7 +132,7 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   type(psb_d_coo_sparse_mat) :: coo_prol, coo_restr
   type(psb_d_csr_sparse_mat) :: acsr1,  acsrf, csr_prol, acsr
   real(psb_dpk_), allocatable :: adiag(:)
-  real(psb_dpk_), allocatable :: arwsum(:)
+  real(psb_dpk_), allocatable :: arwsum(:),l1rwsum(:)
   integer(psb_ipk_)  :: ierr(5)
   logical            :: filter_mat
   integer(psb_ipk_)            :: debug_level, debug_unit, err_act
@@ -141,6 +141,7 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   logical, parameter :: debug_new=.false.
   character(len=80) :: filename
   logical, parameter :: do_timings=.false.
+  logical, parameter :: do_l1correction=.true.
   integer(psb_ipk_), save :: idx_spspmm=-1, idx_phase1=-1, idx_gtrans=-1, idx_phase2=-1, idx_refine=-1
   integer(psb_ipk_), save :: idx_phase3=-1, idx_cdasb=-1, idx_ptap=-1
 
@@ -200,6 +201,21 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   if (info == psb_success_) &
        & call psb_halo(adiag,desc_a,info)
   if (info == psb_success_) call a%cp_to(acsr)
+  ! Get the l1-diagonal of D
+  if (do_l1correction) then
+    allocate(l1rwsum(nrow))
+    call acsr%arwsum(l1rwsum) 
+    if (info == psb_success_) &
+      & call psb_realloc(ncol,l1rwsum,info)
+    if (info == psb_success_) &
+      & call psb_halo(l1rwsum,desc_a,info)
+    ! \tilde{D}_{i,i} = \sum_{j \ne i} |a_{i,j}|
+    !$OMP parallel do private(i) schedule(static)
+    do i=1,size(adiag)
+      adiag(i) = adiag(i) + l1rwsum(i) - abs(adiag(i))
+    end do
+    !$OMP end parallel do 
+  end if
 
   if(info /= psb_success_) then
     call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
@@ -230,7 +246,7 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
 
       enddo
       if (jd == -1) then 
-        write(0,*) 'Wrong input: we need the diagonal!!!!', i
+        if (.not.do_l1correction) write(0,*) 'Wrong input: we need the diagonal!!!!', i
       else
         acsrf%val(jd)=acsrf%val(jd)-tmp
       end if
@@ -240,7 +256,6 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
     call acsrf%clean_zeros(info)
   end if
 
-
   !$OMP parallel do private(i) schedule(static)
   do i=1,size(adiag)
     if (adiag(i) /= dzero) then
@@ -249,7 +264,8 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
       adiag(i) = done
     end if
   end do
-  !$OMP end parallel do 
+  !$OMP end parallel do
+
   if (parms%aggr_omega_alg == amg_eig_est_) then 
 
     if (parms%aggr_eig == amg_max_norm_) then 
@@ -259,7 +275,9 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
       call psb_amx(ctxt,anorm)
       omega = 4.d0/(3.d0*anorm)
       parms%aggr_omega_val = omega 
-
+    else if (do_l1correction) then
+      ! For l1-Jacobi this can be estimated with 1
+      parms%aggr_omega_val = done
     else 
       info = psb_err_internal_error_
       call psb_errpush(info,name,a_err='invalid amg_aggr_eig_')
@@ -323,6 +341,9 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
        & write(debug_unit,*) me,' ',trim(name),&
        & 'Done smooth_aggregate '
   call psb_erractionrestore(err_act)
+
+  if (allocated(l1rwsum)) deallocate(l1rwsum)
+  if (allocated(arwsum)) deallocate(arwsum)
   return
 
 9999 continue