mld2p4-smoother-2side:

mlprec/impl/mld_dmlprec_aply.f90 mlprec/impl/mld_dmlprec_bld.f90 mlprec/mld_d_gs_solver.f90 mlprec/mld_d_onelev_mod.f90 tests/pdegen/runs/ppde.inp First steps towards BW-gs as a 2nd smoother.
9 years ago · 6f06a48d25
parent 51fa5ef069
commit 6f06a48d25
5 changed files with 168 additions and 41 deletions
--- a/mlprec/impl/mld_dmlprec_aply.f90
+++ b/mlprec/impl/mld_dmlprec_aply.f90
@ -406,7 +406,7 @@ contains

    ! Arguments
    integer(psb_ipk_)                  :: level 
-    type(mld_dprec_type), intent(inout) :: p
+    type(mld_dprec_type), target, intent(inout) :: p
    type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:)
    character, intent(in)             :: trans
    real(psb_dpk_),target            :: work(:)
@ -539,7 +539,6 @@ contains
          end if

          ! This is one step of post-smoothing 
-
          if (level < nlev) then 
            call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) 
            if (info /= psb_success_) then
@ -572,7 +571,7 @@ contains
            end if

            sweeps = p%precv(level)%parms%sweeps_post 
-            call p%precv(level)%sm%apply(done,&
+            call p%precv(level)%sm2%apply(done,&
                 & mlprec_wrk(level)%x2l,done,mlprec_wrk(level)%y2l,&
                 & p%precv(level)%base_desc, trans,&
                 & sweeps,work,info)
@ -583,8 +582,9 @@ contains
            end if
           
          else
+            ! Here at coarse level
            sweeps = p%precv(level)%parms%sweeps 
-            call p%precv(level)%sm%apply(done,&
+            call p%precv(level)%sm2%apply(done,&
                 & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,&
                 & p%precv(level)%base_desc, trans,&
                 & sweeps,work,info)
@ -624,7 +624,7 @@ contains
          else
            sweeps = p%precv(level)%parms%sweeps
          end if
-          call p%precv(level)%sm%apply(done,&
+          call p%precv(level)%sm2%apply(done,&
               & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,&
               & p%precv(level)%base_desc, trans,&
               & sweeps,work,info)
@ -830,6 +830,13 @@ contains

      case(mld_twoside_smooth_)

+        ! CHECK
+        if (.not.(associated(p%precv(level)%sm2,p%precv(level)%sm2a))) then
+          write(0,*) 'inner_ml_aply: unassociated sm2 at level ',level
+          call psb_errpush(psb_err_internal_error_,name,&
+               & a_err='Error during restriction')
+          goto 9999
+        end if
        nc2l  = p%precv(level)%base_desc%get_local_cols()
        nr2l  = p%precv(level)%base_desc%get_local_rows()
        allocate(mlprec_wrk(level)%ty(nc2l), mlprec_wrk(level)%tx(nc2l), stat=info)
@ -866,10 +873,19 @@ contains
        else
          sweeps = p%precv(level)%parms%sweeps
        end if
-        if (info == psb_success_) call p%precv(level)%sm%apply(done,&
-             & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,&
-             & p%precv(level)%base_desc, trans,&
-             & sweeps,work,info)
+
+        if (trans == 'N') then 
+          if (info == psb_success_) call p%precv(level)%sm%apply(done,&
+               & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,&
+               & p%precv(level)%base_desc, trans,&
+               & sweeps,work,info)
+        else
+          if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
+               & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,&
+               & p%precv(level)%base_desc, trans,&
+               & sweeps,work,info)
+        end if
+          
        if (info /= psb_success_) then
          call psb_errpush(psb_err_internal_error_,name,&
               & a_err='Error during smoother_apply')
@ -930,10 +946,18 @@ contains
          else
            sweeps = p%precv(level)%parms%sweeps_pre
          end if
-          if (info == psb_success_) call p%precv(level)%sm%apply(done,&
-               & mlprec_wrk(level)%tx,done,mlprec_wrk(level)%y2l,&
-               & p%precv(level)%base_desc, trans,&
-               & sweeps,work,info)
+          if (trans == 'N') then 
+            if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
+                 & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,&
+                 & p%precv(level)%base_desc, trans,&
+                 & sweeps,work,info)
+          else
+            if (info == psb_success_) call p%precv(level)%sm%apply(done,&
+                 & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,&
+                 & p%precv(level)%base_desc, trans,&
+                 & sweeps,work,info)
+          end if
+          
          if (info /= psb_success_) then
            call psb_errpush(psb_err_internal_error_,name,&
                 & a_err='Error during smoother_apply')
@ -1043,7 +1067,7 @@ subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
  level = 1

  call psb_geaxpby(done,x,dzero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
-  call mlprec_wrk(level)%vy2l%set(dzero) 
+  call mlprec_wrk(level)%vy2l%zero() 

  call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info)    

@ -1122,7 +1146,9 @@ contains

    nc2l  = p%precv(level)%base_desc%get_local_cols()
    nr2l  = p%precv(level)%base_desc%get_local_rows()
-
+    if(debug_level > 1) then
+      write(debug_unit,*) me,' inner_ml_aply at level ',level
+    end if

    select case(p%precv(level)%parms%ml_type) 

@ -1250,7 +1276,7 @@ contains


            sweeps = p%precv(level)%parms%sweeps_post 
-            call p%precv(level)%sm%apply(done,&
+            call p%precv(level)%sm2%apply(done,&
                 & mlprec_wrk(level)%vx2l,done,mlprec_wrk(level)%vy2l,&
                 & p%precv(level)%base_desc, trans,&
                 & sweeps,work,info)
@ -1262,7 +1288,7 @@ contains
            
          else
            sweeps = p%precv(level)%parms%sweeps 
-            call p%precv(level)%sm%apply(done,&
+            call p%precv(level)%sm2%apply(done,&
                 & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
                 & p%precv(level)%base_desc, trans,&
                 & sweeps,work,info)
@ -1302,7 +1328,7 @@ contains
          else
            sweeps = p%precv(level)%parms%sweeps
          end if
-          call p%precv(level)%sm%apply(done,&
+          call p%precv(level)%sm2%apply(done,&
               & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
               & p%precv(level)%base_desc, trans,&
               & sweeps,work,info)
@ -1536,13 +1562,20 @@ contains
        else
          sweeps = p%precv(level)%parms%sweeps
        end if
-        if (info == psb_success_) call p%precv(level)%sm%apply(done,&
-             & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
-             & p%precv(level)%base_desc, trans,&
-             & sweeps,work,info)
+        if (trans == 'N') then 
+          if (info == psb_success_) call p%precv(level)%sm%apply(done,&
+               & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
+               & p%precv(level)%base_desc, trans,&
+               & sweeps,work,info)
+        else
+          if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
+               & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
+               & p%precv(level)%base_desc, trans,&
+               & sweeps,work,info)
+        end if
        if (info /= psb_success_) then
          call psb_errpush(psb_err_internal_error_,name,&
-               & a_err='Error during smoother_apply')
+               & a_err='Error during 1st smoother_apply')
          goto 9999
        end if
        
@ -1605,13 +1638,20 @@ contains
          else
            sweeps = p%precv(level)%parms%sweeps_pre
          end if
-          if (info == psb_success_) call p%precv(level)%sm%apply(done,&
-               & mlprec_wrk(level)%vtx,done,mlprec_wrk(level)%vy2l,&
-               & p%precv(level)%base_desc, trans,&
-               & sweeps,work,info)
+          if (trans == 'N') then 
+            if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
+                 & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
+                 & p%precv(level)%base_desc, trans,&
+                 & sweeps,work,info)
+          else
+            if (info == psb_success_) call p%precv(level)%sm%apply(done,&
+                 & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
+                 & p%precv(level)%base_desc, trans,&
+                 & sweeps,work,info)
+          end if
          if (info /= psb_success_) then
            call psb_errpush(psb_err_internal_error_,name,&
-                 & a_err='Error during smoother_apply')
+                 & a_err='Error during 2nd smoother_apply')
            goto 9999
          end if

--- a/mlprec/impl/mld_dmlprec_bld.f90
+++ b/mlprec/impl/mld_dmlprec_bld.f90
@ -495,7 +495,7 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info,amold,vmold,imold)

    call p%precv(i)%sm%build(p%precv(i)%base_a,p%precv(i)%base_desc,&
         & 'F',info,amold=amold,vmold=vmold,imold=imold)
-
+    p%precv(i)%sm2 => p%precv(i)%sm
    if ((info == psb_success_).and.(i>1)) then 
      call p%precv(i)%cnv(info,amold=amold,vmold=vmold,imold=imold)
    end if
--- a/mlprec/mld_d_gs_solver.f90
+++ b/mlprec/mld_d_gs_solver.f90
@ -74,6 +74,16 @@ module mld_d_gs_solver
    procedure, nopass   :: is_iterative => d_gs_solver_is_iterative
  end type mld_d_gs_solver_type

+  type, extends(mld_d_gs_solver_type) :: mld_d_bwgs_solver_type
+  contains
+!!$    procedure, pass(sv) :: build   => mld_d_bwgs_solver_bld
+!!$    procedure, pass(sv) :: apply_v => mld_d_bwgs_solver_apply_vect
+!!$    procedure, pass(sv) :: apply_a => mld_d_bwgs_solver_apply
+    procedure, nopass   :: get_fmt    => d_bwgs_solver_get_fmt
+    procedure, pass(sv) :: descr   => d_bwgs_solver_descr
+  end type mld_d_bwgs_solver_type
+
+  

  private :: d_gs_solver_bld, d_gs_solver_apply, &
       &  d_gs_solver_free,   d_gs_solver_seti, &
@ -82,7 +92,8 @@ module mld_d_gs_solver
       &  d_gs_solver_default, d_gs_solver_dmp, &
       &  d_gs_solver_apply_vect, d_gs_solver_get_nzeros, &
       &  d_gs_solver_get_fmt, d_gs_solver_check,&
-       &  d_gs_solver_is_iterative
+       &  d_gs_solver_is_iterative, &
+       &  d_bwgs_solver_get_fmt, d_bwgs_solver_descr


  interface 
@ -452,10 +463,10 @@ contains
    endif

    if (sv%eps<=dzero) then 
-      write(iout_,*) '  Gauss-Seidel iterative solver with  ',&
+      write(iout_,*) '  Forward Gauss-Seidel iterative solver with  ',&
           &  sv%sweeps,' sweeps'
    else
-      write(iout_,*) '  Gauss-Seidel iterative solver with  tolerance',&
+      write(iout_,*) '  Forward Gauss-Seidel iterative solver with  tolerance',&
           &  sv%eps,' and maxit', sv%sweeps
    end if
    
@ -500,7 +511,7 @@ contains
    implicit none 
    character(len=32)  :: val

-    val = "Gauss-Seidel solver"
+    val = "Forward Gauss-Seidel solver"
  end function d_gs_solver_get_fmt


@ -515,4 +526,51 @@ contains
    val = .true.
  end function d_gs_solver_is_iterative

+  
+  subroutine d_bwgs_solver_descr(sv,info,iout,coarse)
+
+    Implicit None
+
+    ! Arguments
+    class(mld_d_bwgs_solver_type), intent(in) :: sv
+    integer(psb_ipk_), intent(out)             :: info
+    integer(psb_ipk_), intent(in), optional    :: iout
+    logical, intent(in), optional       :: coarse
+
+    ! Local variables
+    integer(psb_ipk_)      :: err_act
+    character(len=20), parameter :: name='mld_d_bwgs_solver_descr'
+    integer(psb_ipk_) :: iout_
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+    if (present(iout)) then 
+      iout_ = iout 
+    else
+      iout_ = 6
+    endif
+
+    if (sv%eps<=dzero) then 
+      write(iout_,*) '  Backward Gauss-Seidel iterative solver with  ',&
+           &  sv%sweeps,' sweeps'
+    else
+      write(iout_,*) '  Backward Gauss-Seidel iterative solver with  tolerance',&
+           &  sv%eps,' and maxit', sv%sweeps
+    end if
+    
+    call psb_erractionrestore(err_act)
+    return
+
+9999 call psb_error_handler(err_act)
+    return
+  end subroutine d_bwgs_solver_descr
+
+  function d_bwgs_solver_get_fmt() result(val)
+    implicit none 
+    character(len=32)  :: val
+
+    val = "Backward Gauss-Seidel solver"
+  end function d_bwgs_solver_get_fmt
+
+
 end module mld_d_gs_solver
--- a/mlprec/mld_d_onelev_mod.f90
+++ b/mlprec/mld_d_onelev_mod.f90
@ -121,7 +121,8 @@ module mld_d_onelev_mod
  !
  !
  type mld_d_onelev_type
-    class(mld_d_base_smoother_type), allocatable :: sm
+    class(mld_d_base_smoother_type), allocatable :: sm, sm2a
+    class(mld_d_base_smoother_type), pointer     :: sm2
    type(mld_dml_parms)              :: parms 
    type(psb_dspmat_type)            :: ac
    integer(psb_ipk_)                :: ac_nz_loc, ac_nz_tot
@ -331,6 +332,8 @@ contains
    val = 0
    if (allocated(lv%sm)) &
         &  val =  lv%sm%get_nzeros()
+    if (allocated(lv%sm2a)) &
+         &  val =  val + lv%sm2a%get_nzeros()
  end function d_base_onelev_get_nzeros

  function d_base_onelev_sizeof(lv) result(val)
@ -343,6 +346,7 @@ contains
    val = val + lv%ac%sizeof()
    val = val + lv%map%sizeof() 
    if (allocated(lv%sm))  val = val + lv%sm%sizeof()
+    if (allocated(lv%sm2a))  val = val + lv%sm2a%sizeof()
  end function d_base_onelev_sizeof


@ -353,7 +357,7 @@ contains

    nullify(lv%base_a) 
    nullify(lv%base_desc) 
-
+    nullify(lv%sm2)
  end subroutine d_base_onelev_nullify

  !
@ -371,7 +375,7 @@ contains
    Implicit None

    ! Arguments
-    class(mld_d_onelev_type), intent(inout) :: lv 
+    class(mld_d_onelev_type), target, intent(inout) :: lv 

    lv%parms%sweeps          = 1
    lv%parms%sweeps_pre      = 1
@ -388,6 +392,12 @@ contains
    lv%parms%aggr_thresh     = dzero

    if (allocated(lv%sm)) call lv%sm%default()
+    if (allocated(lv%sm2a)) then
+      call lv%sm2a%default()
+      lv%sm2 => lv%sm2a
+    else
+      lv%sm2 => lv%sm
+    end if

    return

@ -401,7 +411,7 @@ contains

    ! Arguments
    class(mld_d_onelev_type), target, intent(inout) :: lv 
-    class(mld_d_onelev_type), intent(inout)         :: lvout
+    class(mld_d_onelev_type), target, intent(inout) :: lvout
    integer(psb_ipk_), intent(out)                  :: info
 
    info = psb_success_
@ -413,6 +423,16 @@ contains
        if (info==psb_success_) deallocate(lvout%sm,stat=info)
      end if
    end if
+    if (allocated(lv%sm2a)) then 
+      call lv%sm%clone(lvout%sm2a,info)
+      lvout%sm2 => lvout%sm2a
+    else 
+      if (allocated(lvout%sm2a)) then 
+        call lvout%sm2a%free(info)
+        if (info==psb_success_) deallocate(lvout%sm2a,stat=info)
+      end if
+      lvout%sm2 => lvout%sm
+    end if
    if (info == psb_success_) call lv%parms%clone(lvout%parms,info)
    if (info == psb_success_) call lv%ac%clone(lvout%ac,info)
    if (info == psb_success_) call lv%desc_ac%clone(lvout%desc_ac,info)
@ -428,12 +448,21 @@ contains
  subroutine mld_d_onelev_move_alloc(a, b,info)
    use psb_base_mod
    implicit none
-    type(mld_d_onelev_type), intent(inout) :: a, b
+    type(mld_d_onelev_type), target, intent(inout) :: a, b
    integer(psb_ipk_), intent(out) :: info 

    call b%free(info)
    b%parms  = a%parms
-    call move_alloc(a%sm,b%sm)
+    if (associated(a%sm2,a%sm2a)) then 
+      call move_alloc(a%sm,b%sm)
+      call move_alloc(a%sm2a,b%sm2a)
+      b%sm2 =>b%sm2a
+    else
+      call move_alloc(a%sm,b%sm)
+      call move_alloc(a%sm2a,b%sm2a)
+      b%sm2 =>b%sm
+    end if
+    
    if (info == psb_success_) call psb_move_alloc(a%ac,b%ac,info) 
    if (info == psb_success_) call psb_move_alloc(a%desc_ac,b%desc_ac,info) 
    if (info == psb_success_) call psb_move_alloc(a%map,b%map,info) 
--- a/tests/pdegen/runs/ppde.inp
+++ b/tests/pdegen/runs/ppde.inp
@ -3,7 +3,7 @@ CSR                         ! Storage format CSR COO JAD
 0100                        ! IDIM; domain size is idim**3
 2                           ! ISTOPC
 0100                        ! ITMAX
-1                           ! ITRACE
+10                          ! ITRACE
 30                          ! IRST     (restart for RGMRES and BiCGSTABL) 
 1.d-6                       ! EPS
 3L-MUL-RAS-BJAC4-ILU        ! Descriptive name for preconditioner (up to 40 chars) 
@ -12,7 +12,7 @@ ML                        ! Preconditioner NONE  JACOBI  BJAC  AS ML
 HALO                        ! Restriction operator NONE HALO 
 NONE                        ! Prolongation operator NONE SUM AVG
 GS                         ! Subdomain solver DSCALE ILU MILU ILUT UMF SLU 
-1                           ! sweeps for GS
+2                           ! sweeps for GS
 0                           ! Level-set N for ILU(N), and P for ILUT
 1.d-4                       ! Threshold T for ILU(T,P) 
 1                           ! Smoother/Jacobi  sweeps