diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90 index 358636d8..662f1c1d 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90 @@ -49,7 +49,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& type(psb_d_vect_type),intent(inout) :: y real(psb_dpk_),intent(in) :: alpha,beta character(len=1),intent(in) :: trans - integer(psb_ipk_), intent(in) :: sweeps ! this is ignored here, the polynomial degree dictates the value + integer(psb_ipk_), intent(in) :: sweeps! this is ignored here, the polynomial degree dictates the value real(psb_dpk_),target, intent(inout) :: work(:) type(psb_d_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info @@ -121,13 +121,13 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& call psb_geaxpby(done,x,dzero,r,desc_data,info) call tx%zero() call ty%zero() - call tz%zero() + call tz%zero() select case(sm%variant) case(amg_poly_lottes_) - block + block real(psb_dpk_) :: cz, cr - ! b == x + ! b == x ! x == tx ! do i=1, sm%pdegree @@ -135,27 +135,31 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Z') cz = (2*i*done-3)/(2*i*done+done) cr = (8*i*done-4)/((2*i*done+done)*sm%rho_ba) - ! z_k = cz z_{k-1} + cr ty = cz z_{k-1} + cr Br_{k-1} - call psb_geaxpby(cr,ty,cz,tz,desc_data,info) - ! r_k = b-Ax_k = x -A tx - call psb_geaxpby(done,tz,done,tx,desc_data,info) + if (.false.) then + ! z_k = cz z_{k-1} + cr ty = cz z_{k-1} + cr Br_{k-1} + call psb_geaxpby(cr,ty,cz,tz,desc_data,info) + ! r_k = b-Ax_k = x -A tx + call psb_geaxpby(done,tz,done,tx,desc_data,info) + else + call psb_abgdxyz(cr,cz,done,done,ty,tz,tx,desc_data,info) + end if if (.false.) then call psb_geaxpby(done,x,dzero,r,desc_data,info) call psb_spmm(-done,sm%pa,tx,done,r,desc_data,info,work=aux,trans=trans_) 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 LOTTES ',i,res +!!$ res = psb_genrm2(r,desc_data,info) +!!$ write(0,*) 'Polynomial smoother LOTTES',i,res ! x_k = x_{k-1} + z_k end do end block case(amg_poly_lottes_beta_) - block + block real(psb_dpk_) :: cz, cr - ! b == x + ! b == x ! x == tx ! if (allocated(sm%poly_beta)) then @@ -171,29 +175,33 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Z') cz = (2*i*done-3)/(2*i*done+done) cr = (8*i*done-4)/((2*i*done+done)*sm%rho_ba) - ! z_k = cz z_{k-1} + cr ty = cz z_{k-1} + cr Br_{k-1} - call psb_geaxpby(cr,ty,cz,tz,desc_data,info) - ! r_k = b-Ax_k = x -A tx - call psb_geaxpby(sm%poly_beta(i),tz,done,tx,desc_data,info) + if (.false.) then + ! z_k = cz z_{k-1} + cr ty = cz z_{k-1} + cr Br_{k-1} + call psb_geaxpby(cr,ty,cz,tz,desc_data,info) + ! r_k = b-Ax_k = x -A tx + call psb_geaxpby(sm%poly_beta(i),tz,done,tx,desc_data,info) + else + call psb_abgdxyz(cr,cz,sm%poly_beta(i),done,ty,tz,tx,desc_data,info) + end if if (.false.) then call psb_geaxpby(done,x,dzero,r,desc_data,info) call psb_spmm(-done,sm%pa,tx,done,r,desc_data,info,work=aux,trans=trans_) 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 LOTTES_BETA ',i,res +!!$ res = psb_genrm2(r,desc_data,info) +!!$ write(0,*) 'Polynomial smoother LOTTES_BETA',i,res ! x_k = x_{k-1} + z_k end do end block case(amg_poly_new_) - block + block real(psb_dpk_) :: sigma, theta, delta, rho_old, rho - ! b == x + ! b == x ! x == tx ! - + sm%cf_a = amg_d_poly_a_vect(sm%pdegree) theta = (done+sm%cf_a)/2 delta = (done-sm%cf_a)/2 @@ -201,12 +209,17 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& 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) + if (.false.) then + call psb_geaxpby((done/theta),r,dzero,tz,desc_data,info) + call psb_geaxpby(done,tz,done,tx,desc_data,info) + else + call psb_abgdxyz((done/theta),dzero,done,done,r,tz,tx,desc_data,info) + end if + ! 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') @@ -214,9 +227,14 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& ! ! 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) + if (.false.) then + call psb_geaxpby((2*rho/delta),r,(rho*rho_old),tz,desc_data,info) + call psb_geaxpby(done,tz,done,tx,desc_data,info) + else + call psb_abgdxyz((2*rho/delta),(rho*rho_old),done,done,r,tz,tx,desc_data,info) + end if !!$ res = psb_genrm2(r,desc_data,info) -!!$ write(0,*) 'Polynomial smoother NEW ',i,res +!!$ write(0,*) 'Polynomial smoother ',i,res ! x_k = x_{k-1} + z_k end do end block @@ -239,10 +257,205 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& end if end associate + + + +!!$ if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then +!!$ ! if .not.sv%is_iterative, there's no need to pass init +!!$ call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info) +!!$ +!!$ if (info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,& +!!$ & name,a_err='Error in sub_aply Jacobi Sweeps = 1') +!!$ goto 9999 +!!$ endif +!!$ +!!$ else if (sweeps >= 0) then +!!$ select type (smsv => sm%sv) +!!$ class is (amg_d_diag_solver_type) +!!$ ! +!!$ ! This means we are dealing with a pure Jacobi smoother/solver. +!!$ ! +!!$ associate(tx => wv(1), ty => wv(2)) +!!$ select case (init_) +!!$ case('Z') +!!$ +!!$ call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z') +!!$ +!!$ case('Y') +!!$ call psb_geaxpby(done,x,dzero,tx,desc_data,info) +!!$ call psb_geaxpby(done,y,dzero,ty,desc_data,info) +!!$ call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_) +!!$ call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') +!!$ +!!$ case('U') +!!$ if (.not.present(initu)) then +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ & a_err='missing initu to smoother_apply') +!!$ goto 9999 +!!$ end if +!!$ call psb_geaxpby(done,x,dzero,tx,desc_data,info) +!!$ call psb_geaxpby(done,initu,dzero,ty,desc_data,info) +!!$ call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_) +!!$ call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') +!!$ +!!$ case default +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ & a_err='wrong init to smoother_apply') +!!$ goto 9999 +!!$ end select +!!$ +!!$ do i=1, sweeps-1 +!!$ ! +!!$ ! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)), +!!$ ! where is the diagonal and A the matrix. +!!$ ! +!!$ call psb_geaxpby(done,x,dzero,tx,desc_data,info) +!!$ call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_) +!!$ +!!$ if (info /= psb_success_) exit +!!$ +!!$ call sm%sv%apply(done,tx,done,ty,desc_data,trans_,aux,wv(3:),info,init='Y') +!!$ +!!$ if (info /= psb_success_) exit +!!$ +!!$ if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then +!!$ call psb_geaxpby(done,x,dzero,r,r,desc_data,info) +!!$ call psb_spmm(-done,sm%pa,ty,done,r,desc_data,info) +!!$ res = psb_genrm2(r,desc_data,info) +!!$ if( sm%printres ) then +!!$ call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol) +!!$ end if +!!$ if ( res < sm%tol*resdenum ) then +!!$ if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) & +!!$ & call log_conv("BJAC",me,i,1,res,resdenum,sm%tol) +!!$ exit +!!$ end if +!!$ end if +!!$ +!!$ end do +!!$ +!!$ +!!$ if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) +!!$ +!!$ if (info /= psb_success_) then +!!$ info=psb_err_internal_error_ +!!$ call psb_errpush(info,name,& +!!$ & a_err='subsolve with Jacobi sweeps > 1') +!!$ goto 9999 +!!$ end if +!!$ +!!$ +!!$ end associate +!!$ +!!$ class default +!!$ ! +!!$ ! +!!$ ! Apply multiple sweeps of a block-Jacobi solver +!!$ ! to compute an approximate solution of a linear system. +!!$ ! +!!$ ! +!!$ if (size(wv) < 2) then +!!$ info = psb_err_internal_error_ +!!$ call psb_errpush(info,name,& +!!$ & a_err='invalid wv size in smoother_apply') +!!$ goto 9999 +!!$ end if +!!$ associate(tx => wv(1), ty => wv(2)) +!!$ +!!$ ! +!!$ ! Unroll the first iteration and fold it inside SELECT CASE +!!$ ! this will save one AXPBY and one SPMM when INIT=Z, and will be +!!$ ! significant when sweeps=1 (a common case) +!!$ ! +!!$ select case (init_) +!!$ case('Z') +!!$ +!!$ call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z') +!!$ +!!$ case('Y') +!!$ call psb_geaxpby(done,x,dzero,tx,desc_data,info) +!!$ call psb_geaxpby(done,y,dzero,ty,desc_data,info) +!!$ call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_) +!!$ call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') +!!$ +!!$ case('U') +!!$ if (.not.present(initu)) then +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ & a_err='missing initu to smoother_apply') +!!$ goto 9999 +!!$ end if +!!$ call psb_geaxpby(done,x,dzero,tx,desc_data,info) +!!$ call psb_geaxpby(done,initu,dzero,ty,desc_data,info) +!!$ call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_) +!!$ call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') +!!$ +!!$ case default +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ & a_err='wrong init to smoother_apply') +!!$ goto 9999 +!!$ end select +!!$ +!!$ do i=1, sweeps-1 +!!$ ! +!!$ ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the +!!$ ! block diagonal part and the remaining part of the local matrix +!!$ ! and Y(j) is the approximate solution at sweep j. +!!$ ! +!!$ call psb_geaxpby(done,x,dzero,tx,desc_data,info) +!!$ call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_) +!!$ +!!$ if (info /= psb_success_) exit +!!$ +!!$ call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') +!!$ +!!$ if (info /= psb_success_) exit +!!$ +!!$ if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then +!!$ call psb_geaxpby(done,x,dzero,r,r,desc_data,info) +!!$ call psb_spmm(-done,sm%pa,ty,done,r,desc_data,info) +!!$ res = psb_genrm2(r,desc_data,info) +!!$ if( sm%printres ) then +!!$ call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol) +!!$ end if +!!$ if (res < sm%tol*resdenum ) then +!!$ if( (sm%printres).and.( mod(sm%printiter,sm%checkiter) /=0 ) ) & +!!$ & call log_conv("BJAC",me,i,1,res,resdenum,sm%tol) +!!$ exit +!!$ end if +!!$ end if +!!$ +!!$ end do +!!$ +!!$ if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) +!!$ +!!$ if (info /= psb_success_) then +!!$ info=psb_err_internal_error_ +!!$ call psb_errpush(info,name,& +!!$ & a_err='subsolve with Jacobi sweeps > 1') +!!$ goto 9999 +!!$ end if +!!$ +!!$ end associate +!!$ end select +!!$ +!!$ else +!!$ +!!$ info = psb_err_iarg_neg_ +!!$ call psb_errpush(info,name,& +!!$ & i_err=(/itwo,sweeps,izero,izero,izero/)) +!!$ goto 9999 +!!$ +!!$ endif +!!$ if (.not.(4*n_col <= size(work))) then deallocate(aux) endif +!!$ if(sm%checkres) then +!!$ call psb_gefree(r,desc_data,info) +!!$ end if + call psb_erractionrestore(err_act) return diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 index 8408beff..ac3f1010 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 @@ -56,7 +56,7 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) type(psb_dspmat_type) :: tmpa integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, nzeros type(psb_ctxt_type) :: ctxt - real(psb_dpk_), allocatable :: da(:), dsv(:) + real(psb_dpk_), allocatable :: da(:), dsv(:) integer(psb_ipk_) :: np, me, i, err_act, debug_unit, debug_level character(len=20) :: name='d_poly_smoother_bld', ch_err @@ -75,39 +75,39 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) nrow_a = a%get_nrows() nztota = a%get_nzeros() select case(sm%variant) - case(amg_poly_lottes_) - ! do nothing - case(amg_poly_lottes_beta_) - if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then - call psb_realloc(sm%pdegree,sm%poly_beta,info) - sm%poly_beta(1:sm%pdegree) = amg_d_poly_beta_mat(1:sm%pdegree,sm%pdegree) - else - info = psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='invalid sm%degree for poly_beta') - goto 9999 - end if - case(amg_poly_new_) - write(psb_out_unit,*) "Nella fase di build sm%pdegree = ",sm%pdegree - if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then - !Ok - 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(amg_poly_lottes_) + ! do nothing + case(amg_poly_lottes_beta_) + if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then + call psb_realloc(sm%pdegree,sm%poly_beta,info) + sm%poly_beta(1:sm%pdegree) = amg_d_poly_beta_mat(1:sm%pdegree,sm%pdegree) + else + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='invalid sm%degree for poly_beta') + goto 9999 + end if + case(amg_poly_new_) - case default + if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then + !Ok + 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%variant') + & a_err='invalid sm%degree for poly_a') goto 9999 - end select + end if + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='invalid sm%variant') + goto 9999 + end select sm%pa => a - if (.not.allocated(sm%sv)) then + if (.not.allocated(sm%sv)) then info = psb_err_internal_error_ call psb_errpush(info,name,& & a_err='unallocated sm%sv') @@ -120,7 +120,7 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) goto 9999 end if -!!$ if (.false.) then +!!$ if (.false.) then !!$ select type(ssv => sm%sv) !!$ class is(amg_d_l1_diag_solver_type) !!$ da = a%arwsum(info) @@ -128,7 +128,7 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) !!$ sm%rho_ba = maxval(da(1:n_row)*dsv(1:n_row)) !!$ class default !!$ write(0,*) 'PolySmoother BUILD: only L1-Jacobi/L1-DIAG for now ',ssv%get_fmt() -!!$ sm%rho_ba = done +!!$ sm%rho_ba = done !!$ end select !!$ else if (sm%rho_ba <= dzero) then @@ -152,9 +152,9 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) call sm%sv%apply_v(done,tt,dzero,tz,desc_a,'NoTrans',work,wv,info) ! z_{k+1} = BA q_k do i=1,sm%rho_estimate_iterations znrm = psb_genrm2(tz,desc_a,info) ! znrm = |z_k|_2 - call psb_geaxpby((done/znrm),tz,dzero,tq,desc_a,info) ! q_k = z_k/znrm + call psb_geaxpby((done/znrm),tz,dzero,tq,desc_a,info) ! q_k = z_k/znrm call psb_spmm(done,a,tq,dzero,tt,desc_a,info) ! t_{k+1} = BA q_k - call sm%sv%apply_v(done,tt,dzero,tz,desc_a,'NoTrans',work,wv,info) ! z_{k+1} = B t_{k+1} + call sm%sv%apply_v(done,tt,dzero,tz,desc_a,'NoTrans',work,wv,info) ! z_{k+1} = B t_{k+1} lambda = psb_gedot(tq,tz,desc_a,info) ! lambda = q_k^T z_{k+1} = q_k^T BA q_k !write(0,*) 'BLD: lambda estimate ',i,lambda end do diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_descr.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_descr.f90 index 469e1bb6..9e4803d4 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_descr.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_descr.f90 @@ -1,15 +1,15 @@ -! -! +! +! ! AMG4PSBLAS version 1.0 ! Algebraic Multigrid Package ! based on PSBLAS (Parallel Sparse BLAS version 3.7) -! -! (C) Copyright 2021 -! -! Salvatore Filippone -! Pasqua D'Ambra -! Fabio Durastante -! +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -21,7 +21,7 @@ ! 3. The name of the AMG4PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -33,8 +33,8 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! +! +! subroutine amg_d_poly_smoother_descr(sm,info,iout,coarse,prefix) use psb_base_mod @@ -59,13 +59,13 @@ subroutine amg_d_poly_smoother_descr(sm,info,iout,coarse,prefix) call psb_erractionsave(err_act) info = psb_success_ - if (present(coarse)) then + if (present(coarse)) then coarse_ = coarse else coarse_ = .false. end if - if (present(iout)) then - iout_ = iout + if (present(iout)) then + iout_ = iout else iout_ = psb_out_unit endif diff --git a/amgprec/impl/smoother/amg_s_poly_smoother_apply_vect.f90 b/amgprec/impl/smoother/amg_s_poly_smoother_apply_vect.f90 index b2b7b7df..c76621a9 100644 --- a/amgprec/impl/smoother/amg_s_poly_smoother_apply_vect.f90 +++ b/amgprec/impl/smoother/amg_s_poly_smoother_apply_vect.f90 @@ -49,7 +49,7 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& type(psb_s_vect_type),intent(inout) :: y real(psb_spk_),intent(in) :: alpha,beta character(len=1),intent(in) :: trans - integer(psb_ipk_), intent(in) :: sweeps ! this is ignored here, the polynomial degree dictates the value + integer(psb_ipk_), intent(in) :: sweeps! this is ignored here, the polynomial degree dictates the value real(psb_spk_),target, intent(inout) :: work(:) type(psb_s_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info @@ -121,13 +121,13 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& call psb_geaxpby(sone,x,szero,r,desc_data,info) call tx%zero() call ty%zero() - call tz%zero() + call tz%zero() select case(sm%variant) case(amg_poly_lottes_) - block + block real(psb_spk_) :: cz, cr - ! b == x + ! b == x ! x == tx ! do i=1, sm%pdegree @@ -135,27 +135,31 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& call sm%sv%apply(sone,r,szero,ty,desc_data,trans_,aux,wv(5:),info,init='Z') cz = (2*i*sone-3)/(2*i*sone+sone) cr = (8*i*sone-4)/((2*i*sone+sone)*sm%rho_ba) - ! z_k = cz z_{k-1} + cr ty = cz z_{k-1} + cr Br_{k-1} - call psb_geaxpby(cr,ty,cz,tz,desc_data,info) - ! r_k = b-Ax_k = x -A tx - call psb_geaxpby(sone,tz,sone,tx,desc_data,info) + if (.false.) then + ! z_k = cz z_{k-1} + cr ty = cz z_{k-1} + cr Br_{k-1} + call psb_geaxpby(cr,ty,cz,tz,desc_data,info) + ! r_k = b-Ax_k = x -A tx + call psb_geaxpby(sone,tz,sone,tx,desc_data,info) + else + call psb_abgdxyz(cr,cz,sone,sone,ty,tz,tx,desc_data,info) + end if if (.false.) then call psb_geaxpby(sone,x,szero,r,desc_data,info) call psb_spmm(-sone,sm%pa,tx,sone,r,desc_data,info,work=aux,trans=trans_) else call psb_spmm(-sone,sm%pa,tz,sone,r,desc_data,info,work=aux,trans=trans_) end if -!!$ res = psb_genrm2(r,desc_data,info) -!!$ write(0,*) 'Polynomial smoother LOTTES ',i,res +!!$ res = psb_genrm2(r,desc_data,info) +!!$ write(0,*) 'Polynomial smoother LOTTES',i,res ! x_k = x_{k-1} + z_k end do end block case(amg_poly_lottes_beta_) - block + block real(psb_spk_) :: cz, cr - ! b == x + ! b == x ! x == tx ! if (allocated(sm%poly_beta)) then @@ -171,29 +175,33 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& call sm%sv%apply(sone,r,szero,ty,desc_data,trans_,aux,wv(5:),info,init='Z') cz = (2*i*sone-3)/(2*i*sone+sone) cr = (8*i*sone-4)/((2*i*sone+sone)*sm%rho_ba) - ! z_k = cz z_{k-1} + cr ty = cz z_{k-1} + cr Br_{k-1} - call psb_geaxpby(cr,ty,cz,tz,desc_data,info) - ! r_k = b-Ax_k = x -A tx - call psb_geaxpby(sm%poly_beta(i),tz,sone,tx,desc_data,info) + if (.false.) then + ! z_k = cz z_{k-1} + cr ty = cz z_{k-1} + cr Br_{k-1} + call psb_geaxpby(cr,ty,cz,tz,desc_data,info) + ! r_k = b-Ax_k = x -A tx + call psb_geaxpby(sm%poly_beta(i),tz,sone,tx,desc_data,info) + else + call psb_abgdxyz(cr,cz,sm%poly_beta(i),sone,ty,tz,tx,desc_data,info) + end if if (.false.) then call psb_geaxpby(sone,x,szero,r,desc_data,info) call psb_spmm(-sone,sm%pa,tx,sone,r,desc_data,info,work=aux,trans=trans_) else call psb_spmm(-sone,sm%pa,tz,sone,r,desc_data,info,work=aux,trans=trans_) end if -!!$ res = psb_genrm2(r,desc_data,info) -!!$ write(0,*) 'Polynomial smoother LOTTES_BETA ',i,res +!!$ res = psb_genrm2(r,desc_data,info) +!!$ write(0,*) 'Polynomial smoother LOTTES_BETA',i,res ! x_k = x_{k-1} + z_k end do end block case(amg_poly_new_) - block + block real(psb_spk_) :: sigma, theta, delta, rho_old, rho - ! b == x + ! b == x ! x == tx ! - + sm%cf_a = amg_d_poly_a_vect(sm%pdegree) theta = (sone+sm%cf_a)/2 delta = (sone-sm%cf_a)/2 @@ -201,12 +209,17 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& rho_old = sone/sigma call sm%sv%apply(sone,r,szero,ty,desc_data,trans_,aux,wv(5:),info,init='Z') call psb_geaxpby((sone/sm%rho_ba),ty,szero,r,desc_data,info) - call psb_geaxpby((sone/theta),r,szero,tz,desc_data,info) + if (.false.) then + call psb_geaxpby((sone/theta),r,szero,tz,desc_data,info) + call psb_geaxpby(sone,tz,sone,tx,desc_data,info) + else + call psb_abgdxyz((sone/theta),szero,sone,sone,r,tz,tx,desc_data,info) + end if + ! tz == d do i=1, sm%pdegree - ! x_{k+1} = x_k + d_k - call psb_geaxpby(sone,tz,sone,tx,desc_data,info) ! + ! ! r_{k-1} = r_k - (1/rho(BA)) B A d_k call psb_spmm(sone,sm%pa,tz,szero,ty,desc_data,info,work=aux,trans=trans_) call sm%sv%apply(-sone,ty,sone,r,desc_data,trans_,aux,wv(5:),info,init='Z') @@ -214,9 +227,14 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& ! ! d_{k+1} = (rho rho_old) d_k + 2(rho/delta) r_{k+1} rho = sone/(2*sigma - rho_old) - call psb_geaxpby((2*rho/delta),r,(rho*rho_old),tz,desc_data,info) + if (.false.) then + call psb_geaxpby((2*rho/delta),r,(rho*rho_old),tz,desc_data,info) + call psb_geaxpby(sone,tz,sone,tx,desc_data,info) + else + call psb_abgdxyz((2*rho/delta),(rho*rho_old),sone,sone,r,tz,tx,desc_data,info) + end if !!$ res = psb_genrm2(r,desc_data,info) -!!$ write(0,*) 'Polynomial smoother NEW ',i,res +!!$ write(0,*) 'Polynomial smoother ',i,res ! x_k = x_{k-1} + z_k end do end block @@ -239,10 +257,205 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& end if end associate + + + +!!$ if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then +!!$ ! if .not.sv%is_iterative, there's no need to pass init +!!$ call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info) +!!$ +!!$ if (info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,& +!!$ & name,a_err='Error in sub_aply Jacobi Sweeps = 1') +!!$ goto 9999 +!!$ endif +!!$ +!!$ else if (sweeps >= 0) then +!!$ select type (smsv => sm%sv) +!!$ class is (amg_s_diag_solver_type) +!!$ ! +!!$ ! This means we are dealing with a pure Jacobi smoother/solver. +!!$ ! +!!$ associate(tx => wv(1), ty => wv(2)) +!!$ select case (init_) +!!$ case('Z') +!!$ +!!$ call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Z') +!!$ +!!$ case('Y') +!!$ call psb_geaxpby(sone,x,szero,tx,desc_data,info) +!!$ call psb_geaxpby(sone,y,szero,ty,desc_data,info) +!!$ call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_) +!!$ call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') +!!$ +!!$ case('U') +!!$ if (.not.present(initu)) then +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ & a_err='missing initu to smoother_apply') +!!$ goto 9999 +!!$ end if +!!$ call psb_geaxpby(sone,x,szero,tx,desc_data,info) +!!$ call psb_geaxpby(sone,initu,szero,ty,desc_data,info) +!!$ call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_) +!!$ call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') +!!$ +!!$ case default +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ & a_err='wrong init to smoother_apply') +!!$ goto 9999 +!!$ end select +!!$ +!!$ do i=1, sweeps-1 +!!$ ! +!!$ ! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)), +!!$ ! where is the diagonal and A the matrix. +!!$ ! +!!$ call psb_geaxpby(sone,x,szero,tx,desc_data,info) +!!$ call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_) +!!$ +!!$ if (info /= psb_success_) exit +!!$ +!!$ call sm%sv%apply(sone,tx,sone,ty,desc_data,trans_,aux,wv(3:),info,init='Y') +!!$ +!!$ if (info /= psb_success_) exit +!!$ +!!$ if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then +!!$ call psb_geaxpby(sone,x,szero,r,r,desc_data,info) +!!$ call psb_spmm(-sone,sm%pa,ty,sone,r,desc_data,info) +!!$ res = psb_genrm2(r,desc_data,info) +!!$ if( sm%printres ) then +!!$ call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol) +!!$ end if +!!$ if ( res < sm%tol*resdenum ) then +!!$ if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) & +!!$ & call log_conv("BJAC",me,i,1,res,resdenum,sm%tol) +!!$ exit +!!$ end if +!!$ end if +!!$ +!!$ end do +!!$ +!!$ +!!$ if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) +!!$ +!!$ if (info /= psb_success_) then +!!$ info=psb_err_internal_error_ +!!$ call psb_errpush(info,name,& +!!$ & a_err='subsolve with Jacobi sweeps > 1') +!!$ goto 9999 +!!$ end if +!!$ +!!$ +!!$ end associate +!!$ +!!$ class default +!!$ ! +!!$ ! +!!$ ! Apply multiple sweeps of a block-Jacobi solver +!!$ ! to compute an approximate solution of a linear system. +!!$ ! +!!$ ! +!!$ if (size(wv) < 2) then +!!$ info = psb_err_internal_error_ +!!$ call psb_errpush(info,name,& +!!$ & a_err='invalid wv size in smoother_apply') +!!$ goto 9999 +!!$ end if +!!$ associate(tx => wv(1), ty => wv(2)) +!!$ +!!$ ! +!!$ ! Unroll the first iteration and fold it inside SELECT CASE +!!$ ! this will save one AXPBY and one SPMM when INIT=Z, and will be +!!$ ! significant when sweeps=1 (a common case) +!!$ ! +!!$ select case (init_) +!!$ case('Z') +!!$ +!!$ call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Z') +!!$ +!!$ case('Y') +!!$ call psb_geaxpby(sone,x,szero,tx,desc_data,info) +!!$ call psb_geaxpby(sone,y,szero,ty,desc_data,info) +!!$ call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_) +!!$ call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') +!!$ +!!$ case('U') +!!$ if (.not.present(initu)) then +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ & a_err='missing initu to smoother_apply') +!!$ goto 9999 +!!$ end if +!!$ call psb_geaxpby(sone,x,szero,tx,desc_data,info) +!!$ call psb_geaxpby(sone,initu,szero,ty,desc_data,info) +!!$ call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_) +!!$ call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') +!!$ +!!$ case default +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ & a_err='wrong init to smoother_apply') +!!$ goto 9999 +!!$ end select +!!$ +!!$ do i=1, sweeps-1 +!!$ ! +!!$ ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the +!!$ ! block diagonal part and the remaining part of the local matrix +!!$ ! and Y(j) is the approximate solution at sweep j. +!!$ ! +!!$ call psb_geaxpby(sone,x,szero,tx,desc_data,info) +!!$ call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_) +!!$ +!!$ if (info /= psb_success_) exit +!!$ +!!$ call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') +!!$ +!!$ if (info /= psb_success_) exit +!!$ +!!$ if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then +!!$ call psb_geaxpby(sone,x,szero,r,r,desc_data,info) +!!$ call psb_spmm(-sone,sm%pa,ty,sone,r,desc_data,info) +!!$ res = psb_genrm2(r,desc_data,info) +!!$ if( sm%printres ) then +!!$ call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol) +!!$ end if +!!$ if (res < sm%tol*resdenum ) then +!!$ if( (sm%printres).and.( mod(sm%printiter,sm%checkiter) /=0 ) ) & +!!$ & call log_conv("BJAC",me,i,1,res,resdenum,sm%tol) +!!$ exit +!!$ end if +!!$ end if +!!$ +!!$ end do +!!$ +!!$ if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) +!!$ +!!$ if (info /= psb_success_) then +!!$ info=psb_err_internal_error_ +!!$ call psb_errpush(info,name,& +!!$ & a_err='subsolve with Jacobi sweeps > 1') +!!$ goto 9999 +!!$ end if +!!$ +!!$ end associate +!!$ end select +!!$ +!!$ else +!!$ +!!$ info = psb_err_iarg_neg_ +!!$ call psb_errpush(info,name,& +!!$ & i_err=(/itwo,sweeps,izero,izero,izero/)) +!!$ goto 9999 +!!$ +!!$ endif +!!$ if (.not.(4*n_col <= size(work))) then deallocate(aux) endif +!!$ if(sm%checkres) then +!!$ call psb_gefree(r,desc_data,info) +!!$ end if + call psb_erractionrestore(err_act) return diff --git a/amgprec/impl/smoother/amg_s_poly_smoother_bld.f90 b/amgprec/impl/smoother/amg_s_poly_smoother_bld.f90 index 3b0f5846..44df090f 100644 --- a/amgprec/impl/smoother/amg_s_poly_smoother_bld.f90 +++ b/amgprec/impl/smoother/amg_s_poly_smoother_bld.f90 @@ -56,7 +56,7 @@ subroutine amg_s_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) type(psb_sspmat_type) :: tmpa integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, nzeros type(psb_ctxt_type) :: ctxt - real(psb_spk_), allocatable :: da(:), dsv(:) + real(psb_spk_), allocatable :: da(:), dsv(:) integer(psb_ipk_) :: np, me, i, err_act, debug_unit, debug_level character(len=20) :: name='d_poly_smoother_bld', ch_err @@ -75,39 +75,39 @@ subroutine amg_s_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) nrow_a = a%get_nrows() nztota = a%get_nzeros() select case(sm%variant) - case(amg_poly_lottes_) - ! do nothing - case(amg_poly_lottes_beta_) - if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then - call psb_realloc(sm%pdegree,sm%poly_beta,info) - sm%poly_beta(1:sm%pdegree) = amg_d_poly_beta_mat(1:sm%pdegree,sm%pdegree) - else - info = psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='invalid sm%degree for poly_beta') - goto 9999 - end if - case(amg_poly_new_) - write(psb_out_unit,*) "Nella fase di build sm%pdegree = ",sm%pdegree - if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then - !Ok - 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(amg_poly_lottes_) + ! do nothing + case(amg_poly_lottes_beta_) + if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then + call psb_realloc(sm%pdegree,sm%poly_beta,info) + sm%poly_beta(1:sm%pdegree) = amg_d_poly_beta_mat(1:sm%pdegree,sm%pdegree) + else + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='invalid sm%degree for poly_beta') + goto 9999 + end if + case(amg_poly_new_) - case default + if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then + !Ok + 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%variant') + & a_err='invalid sm%degree for poly_a') goto 9999 - end select + end if + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='invalid sm%variant') + goto 9999 + end select sm%pa => a - if (.not.allocated(sm%sv)) then + if (.not.allocated(sm%sv)) then info = psb_err_internal_error_ call psb_errpush(info,name,& & a_err='unallocated sm%sv') @@ -120,7 +120,7 @@ subroutine amg_s_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) goto 9999 end if -!!$ if (.false.) then +!!$ if (.false.) then !!$ select type(ssv => sm%sv) !!$ class is(amg_s_l1_diag_solver_type) !!$ da = a%arwsum(info) @@ -128,7 +128,7 @@ subroutine amg_s_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) !!$ sm%rho_ba = maxval(da(1:n_row)*dsv(1:n_row)) !!$ class default !!$ write(0,*) 'PolySmoother BUILD: only L1-Jacobi/L1-DIAG for now ',ssv%get_fmt() -!!$ sm%rho_ba = sone +!!$ sm%rho_ba = sone !!$ end select !!$ else if (sm%rho_ba <= szero) then @@ -152,9 +152,9 @@ subroutine amg_s_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) call sm%sv%apply_v(sone,tt,szero,tz,desc_a,'NoTrans',work,wv,info) ! z_{k+1} = BA q_k do i=1,sm%rho_estimate_iterations znrm = psb_genrm2(tz,desc_a,info) ! znrm = |z_k|_2 - call psb_geaxpby((sone/znrm),tz,szero,tq,desc_a,info) ! q_k = z_k/znrm + call psb_geaxpby((sone/znrm),tz,szero,tq,desc_a,info) ! q_k = z_k/znrm call psb_spmm(sone,a,tq,szero,tt,desc_a,info) ! t_{k+1} = BA q_k - call sm%sv%apply_v(sone,tt,szero,tz,desc_a,'NoTrans',work,wv,info) ! z_{k+1} = B t_{k+1} + call sm%sv%apply_v(sone,tt,szero,tz,desc_a,'NoTrans',work,wv,info) ! z_{k+1} = B t_{k+1} lambda = psb_gedot(tq,tz,desc_a,info) ! lambda = q_k^T z_{k+1} = q_k^T BA q_k !write(0,*) 'BLD: lambda estimate ',i,lambda end do diff --git a/amgprec/impl/smoother/amg_s_poly_smoother_descr.f90 b/amgprec/impl/smoother/amg_s_poly_smoother_descr.f90 index 89ba79ab..86152eb9 100644 --- a/amgprec/impl/smoother/amg_s_poly_smoother_descr.f90 +++ b/amgprec/impl/smoother/amg_s_poly_smoother_descr.f90 @@ -1,15 +1,15 @@ -! -! +! +! ! AMG4PSBLAS version 1.0 ! Algebraic Multigrid Package ! based on PSBLAS (Parallel Sparse BLAS version 3.7) -! -! (C) Copyright 2021 -! -! Salvatore Filippone -! Pasqua D'Ambra -! Fabio Durastante -! +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -21,7 +21,7 @@ ! 3. The name of the AMG4PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -33,8 +33,8 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! +! +! subroutine amg_s_poly_smoother_descr(sm,info,iout,coarse,prefix) use psb_base_mod @@ -59,13 +59,13 @@ subroutine amg_s_poly_smoother_descr(sm,info,iout,coarse,prefix) call psb_erractionsave(err_act) info = psb_success_ - if (present(coarse)) then + if (present(coarse)) then coarse_ = coarse else coarse_ = .false. end if - if (present(iout)) then - iout_ = iout + if (present(iout)) then + iout_ = iout else iout_ = psb_out_unit endif