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 662f1c1d..9d348bc5 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90 @@ -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 @@ -157,9 +157,9 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& case(amg_poly_lottes_beta_) - block + block real(psb_dpk_) :: cz, cr - ! b == x + ! b == x ! x == tx ! if (allocated(sm%poly_beta)) then @@ -196,9 +196,9 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& 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) @@ -234,7 +234,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& 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 ',i,res +!!$ write(0,*) 'Polynomial smoother NEW ',i,res ! x_k = x_{k-1} + z_k end do end block @@ -257,205 +257,10 @@ 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 ac3f1010..dd156912 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 @@ -107,7 +107,7 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) 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 9e4803d4..469e1bb6 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 c76621a9..53c2cb43 100644 --- a/amgprec/impl/smoother/amg_s_poly_smoother_apply_vect.f90 +++ b/amgprec/impl/smoother/amg_s_poly_smoother_apply_vect.f90 @@ -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 @@ -157,9 +157,9 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& case(amg_poly_lottes_beta_) - block + block real(psb_spk_) :: cz, cr - ! b == x + ! b == x ! x == tx ! if (allocated(sm%poly_beta)) then @@ -196,9 +196,9 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& 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) @@ -234,7 +234,7 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& 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 ',i,res +!!$ write(0,*) 'Polynomial smoother NEW ',i,res ! x_k = x_{k-1} + z_k end do end block @@ -257,205 +257,10 @@ 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 44df090f..09b01248 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 @@ -107,7 +107,7 @@ subroutine amg_s_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) 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 86152eb9..89ba79ab 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