Added row-preserving filtering

amgprec/amg_base_prec_type.F90

@@ -297,9 +297,10 @@ module amg_base_prec_type
   !
   ! Legal values for entry: amg_aggr_filter_
   !
-  integer(psb_ipk_), parameter :: amg_no_filter_mat_  = 0
-  integer(psb_ipk_), parameter :: amg_filter_mat_     = 1
-  integer(psb_ipk_), parameter :: amg_max_filter_mat_ = amg_filter_mat_
+  integer(psb_ipk_), parameter :: amg_no_filter_mat_   = 0
+  integer(psb_ipk_), parameter :: amg_filter_mat_      = 1
+  integer(psb_ipk_), parameter :: amg_filter_prow_mat_ = 2
+  integer(psb_ipk_), parameter :: amg_max_filter_mat_  = amg_filter_prow_mat_
   !
   ! Legal values for entry: amg_aggr_ord_
   !
@@ -380,7 +381,8 @@ module amg_base_prec_type
        &  aggr_prols(0:4)=(/'unsmoothed    ','smoothed      ',&
        &  'l1-smoothed   ','min energy    ','bizr. smoothed'/)
   character(len=15), parameter, private :: &
-       &  aggr_filters(0:1)=(/'no filtering  ','filtering     '/)
+       &  aggr_filters(0:2)=(/'no filtering  ','filtering     ',&
+       &  'filtering rsum'/)
   character(len=15), parameter, private :: &
        &  matrix_names(0:1)=(/'distributed   ','replicated    '/)
   character(len=18), parameter, private :: &
@@ -591,6 +593,8 @@ contains
       val = amg_eig_est_
     case('FILTER')
       val = amg_filter_mat_
+    case('FILTERROWSUM')
+      val = amg_filter_prow_mat_
     case('NOFILTER','NO_FILTER')
       val = amg_no_filter_mat_
     case('OUTER_SWEEPS')

amgprec/impl/aggregator/amg_caggrmat_smth_bld.f90

@@ -110,6 +110,7 @@ subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
   use amg_base_prec_type
   use amg_c_inner_mod, amg_protect_name => amg_caggrmat_smth_bld
   use amg_c_base_aggregator_mod
+!  use, intrinsic :: ieee_arithmetic
 
   implicit none
 
@@ -192,7 +193,7 @@ subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
 
   naggrm1 = sum(nlaggr(1:me))
   naggrp1 = sum(nlaggr(1:me+1))
-  filter_mat = (parms%aggr_filter == amg_filter_mat_)
+  filter_mat = (parms%aggr_filter == amg_filter_mat_).or.(parms%aggr_filter == amg_filter_prow_mat_)
 
   !
   ! naggr: number of local aggregates
@@ -220,7 +221,7 @@ subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
     ! \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))
+        adiag(i) = adiag(i) + l1rwsum(i) - abs(adiag(i))
     end do
     !$OMP end parallel do
   end if
@@ -257,8 +258,13 @@ subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
       if (jd == -1) then 
         ! if (.not.do_l1correction) 
         write(0,*) 'Wrong input: we need the diagonal!!!!', i
-      else
+      else if (parms%aggr_filter == amg_filter_mat_) then
+        ! We perform filtering in the standard way assuming that A is an M-matrix
         acsrf%val(jd)=acsrf%val(jd)-tmp
+      else if (parms%aggr_filter == amg_filter_prow_mat_) then
+        ! We are probably doing l1-correction, hence we want to preserve the
+        ! row sum of the matrix: note the change in sign
+        acsrf%val(jd)=acsrf%val(jd)+tmp
       end if
     enddo
     !$OMP end parallel do 
@@ -266,7 +272,6 @@ subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
     call acsrf%clean_zeros(info)
   end if
 
-
   !$OMP parallel do private(i) schedule(static)
   do i=1,size(adiag)
     if (adiag(i) /= czero) then
@@ -278,11 +283,11 @@ subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
   !$OMP end parallel do 
   if (parms%aggr_omega_alg == amg_eig_est_) then 
 
-    !if (do_l1correction) then
-    !  ! For l1-Jacobi this can be estimated with 1
-    !  parms%aggr_omega_val = done
-    !
-    if (parms%aggr_eig == amg_max_norm_) then 
+    if ( (parms%aggr_filter == amg_filter_prow_mat_).and.(do_l1correction) ) then
+      ! For l1-Jacobi this can be estimated with 1:
+      ! this makes sense only if we are preserving the row-sum!
+      parms%aggr_omega_val = done
+    else if (parms%aggr_eig == amg_max_norm_) then 
       allocate(arwsum(nrow))
       call acsr%arwsum(arwsum)      
       anorm = maxval(abs(adiag(1:nrow)*arwsum(1:nrow)))

amgprec/impl/aggregator/amg_daggrmat_smth_bld.f90

@@ -110,6 +110,7 @@ subroutine amg_daggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
   use amg_base_prec_type
   use amg_d_inner_mod, amg_protect_name => amg_daggrmat_smth_bld
   use amg_d_base_aggregator_mod
+!  use, intrinsic :: ieee_arithmetic
 
   implicit none
 
@@ -192,7 +193,7 @@ subroutine amg_daggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
 
   naggrm1 = sum(nlaggr(1:me))
   naggrp1 = sum(nlaggr(1:me+1))
-  filter_mat = (parms%aggr_filter == amg_filter_mat_)
+  filter_mat = (parms%aggr_filter == amg_filter_mat_).or.(parms%aggr_filter == amg_filter_prow_mat_)
 
   !
   ! naggr: number of local aggregates
@@ -220,7 +221,7 @@ subroutine amg_daggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
     ! \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))
+        adiag(i) = adiag(i) + l1rwsum(i) - abs(adiag(i))
     end do
     !$OMP end parallel do
   end if
@@ -257,8 +258,13 @@ subroutine amg_daggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
       if (jd == -1) then 
         ! if (.not.do_l1correction) 
         write(0,*) 'Wrong input: we need the diagonal!!!!', i
-      else
+      else if (parms%aggr_filter == amg_filter_mat_) then
+        ! We perform filtering in the standard way assuming that A is an M-matrix
         acsrf%val(jd)=acsrf%val(jd)-tmp
+      else if (parms%aggr_filter == amg_filter_prow_mat_) then
+        ! We are probably doing l1-correction, hence we want to preserve the
+        ! row sum of the matrix: note the change in sign
+        acsrf%val(jd)=acsrf%val(jd)+tmp
       end if
     enddo
     !$OMP end parallel do 
@@ -266,7 +272,6 @@ subroutine amg_daggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
     call acsrf%clean_zeros(info)
   end if
 
-
   !$OMP parallel do private(i) schedule(static)
   do i=1,size(adiag)
     if (adiag(i) /= dzero) then
@@ -278,11 +283,11 @@ subroutine amg_daggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
   !$OMP end parallel do 
   if (parms%aggr_omega_alg == amg_eig_est_) then 
 
-    !if (do_l1correction) then
-    !  ! For l1-Jacobi this can be estimated with 1
-    !  parms%aggr_omega_val = done
-    !
-    if (parms%aggr_eig == amg_max_norm_) then 
+    if ( (parms%aggr_filter == amg_filter_prow_mat_).and.(do_l1correction) ) then
+      ! For l1-Jacobi this can be estimated with 1:
+      ! this makes sense only if we are preserving the row-sum!
+      parms%aggr_omega_val = done
+    else if (parms%aggr_eig == amg_max_norm_) then 
       allocate(arwsum(nrow))
       call acsr%arwsum(arwsum)      
       anorm = maxval(abs(adiag(1:nrow)*arwsum(1:nrow)))

amgprec/impl/aggregator/amg_saggrmat_smth_bld.f90

@@ -110,6 +110,7 @@ subroutine amg_saggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
   use amg_base_prec_type
   use amg_s_inner_mod, amg_protect_name => amg_saggrmat_smth_bld
   use amg_s_base_aggregator_mod
+!  use, intrinsic :: ieee_arithmetic
 
   implicit none
 
@@ -192,7 +193,7 @@ subroutine amg_saggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
 
   naggrm1 = sum(nlaggr(1:me))
   naggrp1 = sum(nlaggr(1:me+1))
-  filter_mat = (parms%aggr_filter == amg_filter_mat_)
+  filter_mat = (parms%aggr_filter == amg_filter_mat_).or.(parms%aggr_filter == amg_filter_prow_mat_)
 
   !
   ! naggr: number of local aggregates
@@ -220,7 +221,7 @@ subroutine amg_saggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
     ! \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))
+        adiag(i) = adiag(i) + l1rwsum(i) - abs(adiag(i))
     end do
     !$OMP end parallel do
   end if
@@ -257,8 +258,13 @@ subroutine amg_saggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
       if (jd == -1) then 
         ! if (.not.do_l1correction) 
         write(0,*) 'Wrong input: we need the diagonal!!!!', i
-      else
+      else if (parms%aggr_filter == amg_filter_mat_) then
+        ! We perform filtering in the standard way assuming that A is an M-matrix
         acsrf%val(jd)=acsrf%val(jd)-tmp
+      else if (parms%aggr_filter == amg_filter_prow_mat_) then
+        ! We are probably doing l1-correction, hence we want to preserve the
+        ! row sum of the matrix: note the change in sign
+        acsrf%val(jd)=acsrf%val(jd)+tmp
       end if
     enddo
     !$OMP end parallel do 
@@ -266,7 +272,6 @@ subroutine amg_saggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
     call acsrf%clean_zeros(info)
   end if
 
-
   !$OMP parallel do private(i) schedule(static)
   do i=1,size(adiag)
     if (adiag(i) /= szero) then
@@ -278,11 +283,11 @@ subroutine amg_saggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
   !$OMP end parallel do 
   if (parms%aggr_omega_alg == amg_eig_est_) then 
 
-    !if (do_l1correction) then
-    !  ! For l1-Jacobi this can be estimated with 1
-    !  parms%aggr_omega_val = done
-    !
-    if (parms%aggr_eig == amg_max_norm_) then 
+    if ( (parms%aggr_filter == amg_filter_prow_mat_).and.(do_l1correction) ) then
+      ! For l1-Jacobi this can be estimated with 1:
+      ! this makes sense only if we are preserving the row-sum!
+      parms%aggr_omega_val = done
+    else if (parms%aggr_eig == amg_max_norm_) then 
       allocate(arwsum(nrow))
       call acsr%arwsum(arwsum)      
       anorm = maxval(abs(adiag(1:nrow)*arwsum(1:nrow)))

amgprec/impl/aggregator/amg_zaggrmat_smth_bld.f90

@@ -110,6 +110,7 @@ subroutine amg_zaggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
   use amg_base_prec_type
   use amg_z_inner_mod, amg_protect_name => amg_zaggrmat_smth_bld
   use amg_z_base_aggregator_mod
+!  use, intrinsic :: ieee_arithmetic
 
   implicit none
 
@@ -192,7 +193,7 @@ subroutine amg_zaggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
 
   naggrm1 = sum(nlaggr(1:me))
   naggrp1 = sum(nlaggr(1:me+1))
-  filter_mat = (parms%aggr_filter == amg_filter_mat_)
+  filter_mat = (parms%aggr_filter == amg_filter_mat_).or.(parms%aggr_filter == amg_filter_prow_mat_)
 
   !
   ! naggr: number of local aggregates
@@ -220,7 +221,7 @@ subroutine amg_zaggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
     ! \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))
+        adiag(i) = adiag(i) + l1rwsum(i) - abs(adiag(i))
     end do
     !$OMP end parallel do
   end if
@@ -257,8 +258,13 @@ subroutine amg_zaggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
       if (jd == -1) then 
         ! if (.not.do_l1correction) 
         write(0,*) 'Wrong input: we need the diagonal!!!!', i
-      else
+      else if (parms%aggr_filter == amg_filter_mat_) then
+        ! We perform filtering in the standard way assuming that A is an M-matrix
         acsrf%val(jd)=acsrf%val(jd)-tmp
+      else if (parms%aggr_filter == amg_filter_prow_mat_) then
+        ! We are probably doing l1-correction, hence we want to preserve the
+        ! row sum of the matrix: note the change in sign
+        acsrf%val(jd)=acsrf%val(jd)+tmp
       end if
     enddo
     !$OMP end parallel do 
@@ -266,7 +272,6 @@ subroutine amg_zaggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
     call acsrf%clean_zeros(info)
   end if
 
-
   !$OMP parallel do private(i) schedule(static)
   do i=1,size(adiag)
     if (adiag(i) /= zzero) then
@@ -278,11 +283,11 @@ subroutine amg_zaggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
   !$OMP end parallel do 
   if (parms%aggr_omega_alg == amg_eig_est_) then 
 
-    !if (do_l1correction) then
-    !  ! For l1-Jacobi this can be estimated with 1
-    !  parms%aggr_omega_val = done
-    !
-    if (parms%aggr_eig == amg_max_norm_) then 
+    if ( (parms%aggr_filter == amg_filter_prow_mat_).and.(do_l1correction) ) then
+      ! For l1-Jacobi this can be estimated with 1:
+      ! this makes sense only if we are preserving the row-sum!
+      parms%aggr_omega_val = done
+    else if (parms%aggr_eig == amg_max_norm_) then 
       allocate(arwsum(nrow))
       call acsr%arwsum(arwsum)      
       anorm = maxval(abs(adiag(1:nrow)*arwsum(1:nrow)))