Defined new variant of polynomial smoother.

amgprec/amg_d_poly_smoother.f90

@@ -63,7 +63,7 @@ module amg_d_poly_smoother
     integer(psb_ipk_)              :: rho_estimate_iterations=10
     type(psb_dspmat_type), pointer :: pa => null()
     real(psb_dpk_), allocatable    :: poly_beta(:)
-    real(psb_dpk_), allocatable    :: poly_a(:)
+    real(psb_dpk_)                 :: cf_a = dzero
     real(psb_dpk_)                 :: rho_ba = -done
   contains
     procedure, pass(sm) :: apply_v => amg_d_poly_smoother_apply_vect

amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90

@@ -63,7 +63,6 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
   integer(psb_ipk_)   :: np, me, i, err_act
   character           :: trans_, init_
   real(psb_dpk_)      :: res, resdenum
-  real(psb_dpk_)      :: cz, cr
   character(len=20)   :: name='d_poly_smoother_apply_v'
 
   call psb_erractionsave(err_act)
@@ -131,62 +130,100 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
 
     select case(sm%variant)
     case(amg_poly_lottes_)
-      !  b == x 
+      block 
-      !  x == tx
+        real(psb_dpk_)      :: cz, cr
-      ! 
+        !  b == x 
-      do i=1, sm%pdegree
+        !  x == tx
-        !   B r_{k-1}
+        ! 
-        call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Z')
+        do i=1, sm%pdegree
-        cz = (2*i*done-3)/(2*i*done+done)
+          !   B r_{k-1}
-        cr = (8*i*done-4)/((2*i*done+done)*sm%rho_ba)
+          call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Z')
-        ! z_k =  cz z_{k-1} + cr ty = cz z_{k-1} + cr Br_{k-1}
+          cz = (2*i*done-3)/(2*i*done+done)
-        call psb_geaxpby(cr,ty,cz,tz,desc_data,info)
+          cr = (8*i*done-4)/((2*i*done+done)*sm%rho_ba)
-        ! r_k =    b-Ax_k  = x -A tx
+          ! z_k =  cz z_{k-1} + cr ty = cz z_{k-1} + cr Br_{k-1}
-        call psb_geaxpby(done,tz,done,tx,desc_data,info)
+          call psb_geaxpby(cr,ty,cz,tz,desc_data,info)
-        if (.false.) then
+          ! r_k =    b-Ax_k  = x -A tx
-          call psb_geaxpby(done,x,dzero,r,desc_data,info)
+          call psb_geaxpby(done,tz,done,tx,desc_data,info)
-          call psb_spmm(-done,sm%pa,tx,done,r,desc_data,info,work=aux,trans=trans_)
+          if (.false.) then
-        else
+            call psb_geaxpby(done,x,dzero,r,desc_data,info)
-          call psb_spmm(-done,sm%pa,tz,done,r,desc_data,info,work=aux,trans=trans_)
+            call psb_spmm(-done,sm%pa,tx,done,r,desc_data,info,work=aux,trans=trans_)
-        end if
+          else
+            call psb_spmm(-done,sm%pa,tz,done,r,desc_data,info,work=aux,trans=trans_)
+          end if
 !!$        res  = psb_genrm2(r,desc_data,info)
 !!$        write(0,*) 'Polynomial smoother ',i,res
-        ! x_k =  x_{k-1} + z_k
+          ! x_k =  x_{k-1} + z_k
-      end do
+        end do
+      end block
 
     case(amg_poly_lottes_beta_)
 
-      !  b == x 
+      block 
-      !  x == tx
+        real(psb_dpk_)      :: cz, cr
-      ! 
+        !  b == x 
-      do i=1, sm%pdegree
+        !  x == tx
-        !   B r_{k-1}
+        !
-        call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Z')
+        do i=1, sm%pdegree
-        cz = (2*i*done-3)/(2*i*done+done)
+          !   B r_{k-1}
-        cr = (8*i*done-4)/((2*i*done+done)*sm%rho_ba)
+          call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Z')
-        ! z_k =  cz z_{k-1} + cr ty = cz z_{k-1} + cr Br_{k-1}
+          cz = (2*i*done-3)/(2*i*done+done)
-        call psb_geaxpby(cr,ty,cz,tz,desc_data,info)
+          cr = (8*i*done-4)/((2*i*done+done)*sm%rho_ba)
-        ! r_k =    b-Ax_k  = x -A tx
+          ! z_k =  cz z_{k-1} + cr ty = cz z_{k-1} + cr Br_{k-1}
-        call psb_geaxpby(sm%poly_beta(i),tz,done,tx,desc_data,info)
+          call psb_geaxpby(cr,ty,cz,tz,desc_data,info)
-        if (.false.) then
+          ! r_k =    b-Ax_k  = x -A tx
-          call psb_geaxpby(done,x,dzero,r,desc_data,info)
+          call psb_geaxpby(sm%poly_beta(i),tz,done,tx,desc_data,info)
-          call psb_spmm(-done,sm%pa,tx,done,r,desc_data,info,work=aux,trans=trans_)
+          if (.false.) then
-        else
+            call psb_geaxpby(done,x,dzero,r,desc_data,info)
-          call psb_spmm(-done,sm%pa,tz,done,r,desc_data,info,work=aux,trans=trans_)
+            call psb_spmm(-done,sm%pa,tx,done,r,desc_data,info,work=aux,trans=trans_)
-        end if
+          else
+            call psb_spmm(-done,sm%pa,tz,done,r,desc_data,info,work=aux,trans=trans_)
+          end if
 !!$        res  = psb_genrm2(r,desc_data,info)
 !!$        write(0,*) 'Polynomial smoother ',i,res
-        ! x_k =  x_{k-1} + z_k
+          ! x_k =  x_{k-1} + z_k
-      end do
+        end do
-
+      end block
       
+    case(amg_poly_new_)
+      block 
+        real(psb_dpk_)      :: sigma, theta, delta, rho_old, rho
+        !  b == x 
+        !  x == tx
+        !
+        theta = (done+sm%cf_a)/2
+        delta = (done-sm%cf_a)/2
+        sigma = theta/delta
+        rho_old = done/sigma
+        call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Z')
+        call psb_geaxpby((done/sm%rho_ba),ty,dzero,r,desc_data,info)
+        call psb_geaxpby((done/theta),r,dzero,tz,desc_data,info)
+        ! tz == d
+        do i=1, sm%pdegree
+          ! x_{k+1} = x_k + d_k
+          call psb_geaxpby(done,tz,done,tx,desc_data,info)
+          ! 
+          ! r_{k-1} = r_k - (1/rho(BA)) B A d_k
+          call psb_spmm(done,sm%pa,tz,dzero,ty,desc_data,info,work=aux,trans=trans_)
+          call sm%sv%apply(-done,ty,done,r,desc_data,trans_,aux,wv(5:),info,init='Z')
+
+          !
+          ! d_{k+1} = (rho rho_old) d_k + 2(rho/delta) r_{k+1}
+          rho = done/(2*sigma - rho_old)
+          call psb_geaxpby((2*rho/delta),r,(rho*rho_old),tz,desc_data,info)
+!!$          res  = psb_genrm2(r,desc_data,info)
+!!$          write(0,*) 'Polynomial smoother ',i,res
+          ! x_k =  x_{k-1} + z_k
+        end do
+      end block
+
+
     case default
       info=psb_err_internal_error_
       call psb_errpush(info,name,&
            & a_err='wrong polynomial variant')
       goto 9999
     end select
-    
+
     if (info == psb_success_) call psb_geaxpby(alpha,tx,beta,y,desc_data,info)
 
     if (info /= psb_success_) then

amgprec/impl/smoother/amg_d_poly_smoother_bld.f90

@@ -90,15 +90,16 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
     end if
   case(amg_poly_new_)
 
-    if ((1<=sm%pdegree).and.(sm%pdegree<=6)) then
+    if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then
-      call psb_realloc(sm%pdegree,sm%poly_a,info)
+      !Ok
-      sm%poly_a(1:sm%pdegree) = amg_d_poly_a_vect(1:sm%pdegree)
+      sm%cf_a = amg_d_poly_a_vect(sm%pdegree)
     else
       info = psb_err_internal_error_
       call psb_errpush(info,name,&
            & a_err='invalid sm%degree for poly_a')
       goto 9999
     end if
+    
   case default  
     info = psb_err_internal_error_
     call psb_errpush(info,name,&

amgprec/impl/smoother/amg_d_poly_smoother_descr.f90

@@ -88,6 +88,11 @@ subroutine amg_d_poly_smoother_descr(sm,info,iout,coarse,prefix)
     write(iout_,*) trim(prefix_), '         Degree:        ',sm%pdegree
     write(iout_,*) trim(prefix_), '         rho_ba:        ',sm%rho_ba
     if (allocated(sm%poly_beta)) write(iout_,*) trim(prefix_), '   Coefficients:        ',sm%poly_beta(1:sm%pdegree)
+  case(amg_poly_new_)
+    write(iout_,*) trim(prefix_), '        variant:        ','POLY_NEW'
+    write(iout_,*) trim(prefix_), '         Degree:        ',sm%pdegree
+    write(iout_,*) trim(prefix_), '         rho_ba:        ',sm%rho_ba
+    write(iout_,*) trim(prefix_), '   Coefficient:        ',sm%cf_a
   case default
     write(iout_,*) trim(prefix_), '        variant:        ','UNKNOWN???'
   end select