From da6bde21307cd6aa4e98759e0d3510140a3ef8b6 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Fri, 5 Aug 2016 10:31:17 +0000 Subject: [PATCH] mld2p4-2: mlprec/impl/smoother/mld_c_jac_smoother_apply.f90 mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90 mlprec/impl/smoother/mld_d_jac_smoother_apply.f90 mlprec/impl/smoother/mld_d_jac_smoother_apply_vect.f90 mlprec/impl/smoother/mld_s_jac_smoother_apply.f90 mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90 mlprec/impl/smoother/mld_z_jac_smoother_apply.f90 mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90 Fixed application for various cases SWEEPS >=1. --- .../smoother/mld_c_jac_smoother_apply.f90 | 92 +++++++++--------- .../mld_c_jac_smoother_apply_vect.f90 | 96 +++++++++---------- .../smoother/mld_d_jac_smoother_apply.f90 | 92 +++++++++--------- .../mld_d_jac_smoother_apply_vect.f90 | 96 +++++++++---------- .../smoother/mld_s_jac_smoother_apply.f90 | 92 +++++++++--------- .../mld_s_jac_smoother_apply_vect.f90 | 96 +++++++++---------- .../smoother/mld_z_jac_smoother_apply.f90 | 92 +++++++++--------- .../mld_z_jac_smoother_apply_vect.f90 | 96 +++++++++---------- 8 files changed, 364 insertions(+), 388 deletions(-) diff --git a/mlprec/impl/smoother/mld_c_jac_smoother_apply.f90 b/mlprec/impl/smoother/mld_c_jac_smoother_apply.f90 index 3d842f29..04a49eec 100644 --- a/mlprec/impl/smoother/mld_c_jac_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_c_jac_smoother_apply.f90 @@ -55,7 +55,7 @@ subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,& integer(psb_ipk_) :: n_row,n_col complex(psb_spk_), allocatable :: tx(:),ty(:) - complex(psb_spk_), pointer :: ww(:), aux(:) + complex(psb_spk_), pointer :: aux(:) integer(psb_ipk_) :: ictxt,np,me,i, err_act character :: trans_, init_ character(len=20) :: name='c_jac_smoother_apply' @@ -89,71 +89,73 @@ subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,& n_row = desc_data%get_local_rows() n_col = desc_data%get_local_cols() - if (n_col <= size(work)) then - ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then - aux => work(n_col+1:) - else - allocate(aux(4*n_col),stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,& - & i_err=(/4*n_col,izero,izero,izero,izero/),& - & a_err='complex(psb_spk_)') - goto 9999 - end if - endif + if (4*n_col <= size(work)) then + aux => work(:) else - allocate(ww(n_col),aux(4*n_col),stat=info) + allocate(aux(4*n_col),stat=info) if (info /= psb_success_) then info=psb_err_alloc_request_ call psb_errpush(info,name,& - & i_err=(/5*n_col,izero,izero,izero,izero/),& + & i_err=(/4*n_col,izero,izero,izero,izero/),& & a_err='complex(psb_spk_)') goto 9999 end if endif - -!!$ if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_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,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 >= 1) then + + if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_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,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 ! ! ! Apply multiple sweeps of a block-Jacobi solver ! to compute an approximate solution of a linear system. ! ! - call psb_geasb(tx,desc_data,info) call psb_geasb(ty,desc_data,info) + ! + ! 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') - ty(:) = czero + + call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,info,init='Z') + case('Y') + call psb_geaxpby(cone,x,czero,tx,desc_data,info) call psb_geaxpby(cone,y,czero,ty,desc_data,info) + call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,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(cone,x,czero,tx,desc_data,info) call psb_geaxpby(cone,initu,czero,ty,desc_data,info) + call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,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 + + 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 @@ -186,23 +188,17 @@ subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,& goto 9999 end if -!!$ else -!!$ -!!$ info = psb_err_iarg_neg_ -!!$ call psb_errpush(info,name,& -!!$ & i_err=(/itwo,sweeps,izero,izero,izero/)) -!!$ goto 9999 -!!$ -!!$ endif + else + + info = psb_err_iarg_neg_ + call psb_errpush(info,name,& + & i_err=(/itwo,sweeps,izero,izero,izero/)) + goto 9999 + endif - if (n_col <= size(work)) then - if ((4*n_col+n_col) <= size(work)) then - else - deallocate(aux) - endif - else - deallocate(ww,aux) + if (.not.(4*n_col <= size(work))) then + deallocate(aux) endif call psb_erractionrestore(err_act) diff --git a/mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90 index 885a16df..d95e61c4 100644 --- a/mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90 @@ -56,10 +56,10 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& integer(psb_ipk_) :: n_row,n_col type(psb_c_vect_type) :: tx, ty - complex(psb_spk_), pointer :: ww(:), aux(:) + complex(psb_spk_), pointer :: aux(:) integer(psb_ipk_) :: ictxt,np,me,i, err_act character :: trans_, init_ - character(len=20) :: name='c_jac_smoother_apply' + character(len=20) :: name='c_jac_smoother_apply_v' call psb_erractionsave(err_act) @@ -90,71 +90,74 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& n_row = desc_data%get_local_rows() n_col = desc_data%get_local_cols() - if (n_col <= size(work)) then - ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then - aux => work(n_col+1:) - else - allocate(aux(4*n_col),stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,& - & i_err=(/4*n_col,izero,izero,izero,izero/),& - & a_err='complex(psb_spk_)') - goto 9999 - end if - endif + if (4*n_col <= size(work)) then + aux => work(:) else - allocate(ww(n_col),aux(4*n_col),stat=info) + allocate(aux(4*n_col),stat=info) if (info /= psb_success_) then info=psb_err_alloc_request_ call psb_errpush(info,name,& - & i_err=(/5*n_col,izero,izero,izero,izero/),& + & i_err=(/4*n_col,izero,izero,izero,izero/),& & a_err='complex(psb_spk_)') goto 9999 end if endif - -!!$ if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_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,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 + + if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_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,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 ! ! ! Apply multiple sweeps of a block-Jacobi solver ! to compute an approximate solution of a linear system. ! ! - call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.) call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.) + + ! + ! 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 ty%zero() + + call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,info,init='Z') + case('Y') + call psb_geaxpby(cone,x,czero,tx,desc_data,info) call psb_geaxpby(cone,y,czero,ty,desc_data,info) + call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,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(cone,x,czero,tx,desc_data,info) call psb_geaxpby(cone,initu,czero,ty,desc_data,info) + call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,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 + + 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 @@ -188,22 +191,17 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& goto 9999 end if -!!$ else -!!$ -!!$ info = psb_err_iarg_neg_ -!!$ call psb_errpush(info,name,& -!!$ & i_err=(/itwo,sweeps,izero,izero,izero/)) -!!$ goto 9999 -!!$ -!!$ endif - - if (n_col <= size(work)) then - if ((4*n_col+n_col) <= size(work)) then - else - deallocate(aux) - endif else - deallocate(ww,aux) + + 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 call psb_erractionrestore(err_act) diff --git a/mlprec/impl/smoother/mld_d_jac_smoother_apply.f90 b/mlprec/impl/smoother/mld_d_jac_smoother_apply.f90 index 2d4cbd4a..eaf81d75 100644 --- a/mlprec/impl/smoother/mld_d_jac_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_d_jac_smoother_apply.f90 @@ -55,7 +55,7 @@ subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,& integer(psb_ipk_) :: n_row,n_col real(psb_dpk_), allocatable :: tx(:),ty(:) - real(psb_dpk_), pointer :: ww(:), aux(:) + real(psb_dpk_), pointer :: aux(:) integer(psb_ipk_) :: ictxt,np,me,i, err_act character :: trans_, init_ character(len=20) :: name='d_jac_smoother_apply' @@ -89,71 +89,73 @@ subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,& n_row = desc_data%get_local_rows() n_col = desc_data%get_local_cols() - if (n_col <= size(work)) then - ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then - aux => work(n_col+1:) - else - allocate(aux(4*n_col),stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,& - & i_err=(/4*n_col,izero,izero,izero,izero/),& - & a_err='real(psb_dpk_)') - goto 9999 - end if - endif + if (4*n_col <= size(work)) then + aux => work(:) else - allocate(ww(n_col),aux(4*n_col),stat=info) + allocate(aux(4*n_col),stat=info) if (info /= psb_success_) then info=psb_err_alloc_request_ call psb_errpush(info,name,& - & i_err=(/5*n_col,izero,izero,izero,izero/),& + & i_err=(/4*n_col,izero,izero,izero,izero/),& & a_err='real(psb_dpk_)') goto 9999 end if endif - -!!$ if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_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,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 >= 1) then + + if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_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,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 ! ! ! Apply multiple sweeps of a block-Jacobi solver ! to compute an approximate solution of a linear system. ! ! - call psb_geasb(tx,desc_data,info) call psb_geasb(ty,desc_data,info) + ! + ! 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') - ty(:) = dzero + + call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,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,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,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 + + 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 @@ -186,23 +188,17 @@ subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,& goto 9999 end if -!!$ else -!!$ -!!$ info = psb_err_iarg_neg_ -!!$ call psb_errpush(info,name,& -!!$ & i_err=(/itwo,sweeps,izero,izero,izero/)) -!!$ goto 9999 -!!$ -!!$ endif + else + + info = psb_err_iarg_neg_ + call psb_errpush(info,name,& + & i_err=(/itwo,sweeps,izero,izero,izero/)) + goto 9999 + endif - if (n_col <= size(work)) then - if ((4*n_col+n_col) <= size(work)) then - else - deallocate(aux) - endif - else - deallocate(ww,aux) + if (.not.(4*n_col <= size(work))) then + deallocate(aux) endif call psb_erractionrestore(err_act) diff --git a/mlprec/impl/smoother/mld_d_jac_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_d_jac_smoother_apply_vect.f90 index f138f929..fb796a47 100644 --- a/mlprec/impl/smoother/mld_d_jac_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_d_jac_smoother_apply_vect.f90 @@ -56,10 +56,10 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& integer(psb_ipk_) :: n_row,n_col type(psb_d_vect_type) :: tx, ty - real(psb_dpk_), pointer :: ww(:), aux(:) + real(psb_dpk_), pointer :: aux(:) integer(psb_ipk_) :: ictxt,np,me,i, err_act character :: trans_, init_ - character(len=20) :: name='d_jac_smoother_apply' + character(len=20) :: name='d_jac_smoother_apply_v' call psb_erractionsave(err_act) @@ -90,71 +90,74 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& n_row = desc_data%get_local_rows() n_col = desc_data%get_local_cols() - if (n_col <= size(work)) then - ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then - aux => work(n_col+1:) - else - allocate(aux(4*n_col),stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,& - & i_err=(/4*n_col,izero,izero,izero,izero/),& - & a_err='real(psb_dpk_)') - goto 9999 - end if - endif + if (4*n_col <= size(work)) then + aux => work(:) else - allocate(ww(n_col),aux(4*n_col),stat=info) + allocate(aux(4*n_col),stat=info) if (info /= psb_success_) then info=psb_err_alloc_request_ call psb_errpush(info,name,& - & i_err=(/5*n_col,izero,izero,izero,izero/),& + & i_err=(/4*n_col,izero,izero,izero,izero/),& & a_err='real(psb_dpk_)') goto 9999 end if endif - -!!$ if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_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,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 + + if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_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,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 ! ! ! Apply multiple sweeps of a block-Jacobi solver ! to compute an approximate solution of a linear system. ! ! - call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.) call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.) + + ! + ! 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 ty%zero() + + call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,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,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,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 + + 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 @@ -188,22 +191,17 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& goto 9999 end if -!!$ else -!!$ -!!$ info = psb_err_iarg_neg_ -!!$ call psb_errpush(info,name,& -!!$ & i_err=(/itwo,sweeps,izero,izero,izero/)) -!!$ goto 9999 -!!$ -!!$ endif - - if (n_col <= size(work)) then - if ((4*n_col+n_col) <= size(work)) then - else - deallocate(aux) - endif else - deallocate(ww,aux) + + 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 call psb_erractionrestore(err_act) diff --git a/mlprec/impl/smoother/mld_s_jac_smoother_apply.f90 b/mlprec/impl/smoother/mld_s_jac_smoother_apply.f90 index 3f2f10e8..fbb9bcf8 100644 --- a/mlprec/impl/smoother/mld_s_jac_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_s_jac_smoother_apply.f90 @@ -55,7 +55,7 @@ subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,& integer(psb_ipk_) :: n_row,n_col real(psb_spk_), allocatable :: tx(:),ty(:) - real(psb_spk_), pointer :: ww(:), aux(:) + real(psb_spk_), pointer :: aux(:) integer(psb_ipk_) :: ictxt,np,me,i, err_act character :: trans_, init_ character(len=20) :: name='s_jac_smoother_apply' @@ -89,71 +89,73 @@ subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,& n_row = desc_data%get_local_rows() n_col = desc_data%get_local_cols() - if (n_col <= size(work)) then - ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then - aux => work(n_col+1:) - else - allocate(aux(4*n_col),stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,& - & i_err=(/4*n_col,izero,izero,izero,izero/),& - & a_err='real(psb_spk_)') - goto 9999 - end if - endif + if (4*n_col <= size(work)) then + aux => work(:) else - allocate(ww(n_col),aux(4*n_col),stat=info) + allocate(aux(4*n_col),stat=info) if (info /= psb_success_) then info=psb_err_alloc_request_ call psb_errpush(info,name,& - & i_err=(/5*n_col,izero,izero,izero,izero/),& + & i_err=(/4*n_col,izero,izero,izero,izero/),& & a_err='real(psb_spk_)') goto 9999 end if endif - -!!$ if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_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,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 >= 1) then + + if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_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,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 ! ! ! Apply multiple sweeps of a block-Jacobi solver ! to compute an approximate solution of a linear system. ! ! - call psb_geasb(tx,desc_data,info) call psb_geasb(ty,desc_data,info) + ! + ! 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') - ty(:) = szero + + call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,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,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,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 + + 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 @@ -186,23 +188,17 @@ subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,& goto 9999 end if -!!$ else -!!$ -!!$ info = psb_err_iarg_neg_ -!!$ call psb_errpush(info,name,& -!!$ & i_err=(/itwo,sweeps,izero,izero,izero/)) -!!$ goto 9999 -!!$ -!!$ endif + else + + info = psb_err_iarg_neg_ + call psb_errpush(info,name,& + & i_err=(/itwo,sweeps,izero,izero,izero/)) + goto 9999 + endif - if (n_col <= size(work)) then - if ((4*n_col+n_col) <= size(work)) then - else - deallocate(aux) - endif - else - deallocate(ww,aux) + if (.not.(4*n_col <= size(work))) then + deallocate(aux) endif call psb_erractionrestore(err_act) diff --git a/mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90 index 27e5fd6a..96b96e57 100644 --- a/mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90 @@ -56,10 +56,10 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& integer(psb_ipk_) :: n_row,n_col type(psb_s_vect_type) :: tx, ty - real(psb_spk_), pointer :: ww(:), aux(:) + real(psb_spk_), pointer :: aux(:) integer(psb_ipk_) :: ictxt,np,me,i, err_act character :: trans_, init_ - character(len=20) :: name='s_jac_smoother_apply' + character(len=20) :: name='s_jac_smoother_apply_v' call psb_erractionsave(err_act) @@ -90,71 +90,74 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& n_row = desc_data%get_local_rows() n_col = desc_data%get_local_cols() - if (n_col <= size(work)) then - ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then - aux => work(n_col+1:) - else - allocate(aux(4*n_col),stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,& - & i_err=(/4*n_col,izero,izero,izero,izero/),& - & a_err='real(psb_spk_)') - goto 9999 - end if - endif + if (4*n_col <= size(work)) then + aux => work(:) else - allocate(ww(n_col),aux(4*n_col),stat=info) + allocate(aux(4*n_col),stat=info) if (info /= psb_success_) then info=psb_err_alloc_request_ call psb_errpush(info,name,& - & i_err=(/5*n_col,izero,izero,izero,izero/),& + & i_err=(/4*n_col,izero,izero,izero,izero/),& & a_err='real(psb_spk_)') goto 9999 end if endif - -!!$ if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_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,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 + + if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_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,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 ! ! ! Apply multiple sweeps of a block-Jacobi solver ! to compute an approximate solution of a linear system. ! ! - call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.) call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.) + + ! + ! 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 ty%zero() + + call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,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,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,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 + + 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 @@ -188,22 +191,17 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& goto 9999 end if -!!$ else -!!$ -!!$ info = psb_err_iarg_neg_ -!!$ call psb_errpush(info,name,& -!!$ & i_err=(/itwo,sweeps,izero,izero,izero/)) -!!$ goto 9999 -!!$ -!!$ endif - - if (n_col <= size(work)) then - if ((4*n_col+n_col) <= size(work)) then - else - deallocate(aux) - endif else - deallocate(ww,aux) + + 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 call psb_erractionrestore(err_act) diff --git a/mlprec/impl/smoother/mld_z_jac_smoother_apply.f90 b/mlprec/impl/smoother/mld_z_jac_smoother_apply.f90 index 7204e2c0..7c651894 100644 --- a/mlprec/impl/smoother/mld_z_jac_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_z_jac_smoother_apply.f90 @@ -55,7 +55,7 @@ subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,& integer(psb_ipk_) :: n_row,n_col complex(psb_dpk_), allocatable :: tx(:),ty(:) - complex(psb_dpk_), pointer :: ww(:), aux(:) + complex(psb_dpk_), pointer :: aux(:) integer(psb_ipk_) :: ictxt,np,me,i, err_act character :: trans_, init_ character(len=20) :: name='z_jac_smoother_apply' @@ -89,71 +89,73 @@ subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,& n_row = desc_data%get_local_rows() n_col = desc_data%get_local_cols() - if (n_col <= size(work)) then - ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then - aux => work(n_col+1:) - else - allocate(aux(4*n_col),stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,& - & i_err=(/4*n_col,izero,izero,izero,izero/),& - & a_err='complex(psb_dpk_)') - goto 9999 - end if - endif + if (4*n_col <= size(work)) then + aux => work(:) else - allocate(ww(n_col),aux(4*n_col),stat=info) + allocate(aux(4*n_col),stat=info) if (info /= psb_success_) then info=psb_err_alloc_request_ call psb_errpush(info,name,& - & i_err=(/5*n_col,izero,izero,izero,izero/),& + & i_err=(/4*n_col,izero,izero,izero,izero/),& & a_err='complex(psb_dpk_)') goto 9999 end if endif - -!!$ if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_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,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 >= 1) then + + if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_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,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 ! ! ! Apply multiple sweeps of a block-Jacobi solver ! to compute an approximate solution of a linear system. ! ! - call psb_geasb(tx,desc_data,info) call psb_geasb(ty,desc_data,info) + ! + ! 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') - ty(:) = zzero + + call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,info,init='Z') + case('Y') + call psb_geaxpby(zone,x,zzero,tx,desc_data,info) call psb_geaxpby(zone,y,zzero,ty,desc_data,info) + call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,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(zone,x,zzero,tx,desc_data,info) call psb_geaxpby(zone,initu,zzero,ty,desc_data,info) + call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,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 + + 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 @@ -186,23 +188,17 @@ subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,& goto 9999 end if -!!$ else -!!$ -!!$ info = psb_err_iarg_neg_ -!!$ call psb_errpush(info,name,& -!!$ & i_err=(/itwo,sweeps,izero,izero,izero/)) -!!$ goto 9999 -!!$ -!!$ endif + else + + info = psb_err_iarg_neg_ + call psb_errpush(info,name,& + & i_err=(/itwo,sweeps,izero,izero,izero/)) + goto 9999 + endif - if (n_col <= size(work)) then - if ((4*n_col+n_col) <= size(work)) then - else - deallocate(aux) - endif - else - deallocate(ww,aux) + if (.not.(4*n_col <= size(work))) then + deallocate(aux) endif call psb_erractionrestore(err_act) diff --git a/mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90 index 726b1687..020f4cc2 100644 --- a/mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90 @@ -56,10 +56,10 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& integer(psb_ipk_) :: n_row,n_col type(psb_z_vect_type) :: tx, ty - complex(psb_dpk_), pointer :: ww(:), aux(:) + complex(psb_dpk_), pointer :: aux(:) integer(psb_ipk_) :: ictxt,np,me,i, err_act character :: trans_, init_ - character(len=20) :: name='z_jac_smoother_apply' + character(len=20) :: name='z_jac_smoother_apply_v' call psb_erractionsave(err_act) @@ -90,71 +90,74 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& n_row = desc_data%get_local_rows() n_col = desc_data%get_local_cols() - if (n_col <= size(work)) then - ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then - aux => work(n_col+1:) - else - allocate(aux(4*n_col),stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,& - & i_err=(/4*n_col,izero,izero,izero,izero/),& - & a_err='complex(psb_dpk_)') - goto 9999 - end if - endif + if (4*n_col <= size(work)) then + aux => work(:) else - allocate(ww(n_col),aux(4*n_col),stat=info) + allocate(aux(4*n_col),stat=info) if (info /= psb_success_) then info=psb_err_alloc_request_ call psb_errpush(info,name,& - & i_err=(/5*n_col,izero,izero,izero,izero/),& + & i_err=(/4*n_col,izero,izero,izero,izero/),& & a_err='complex(psb_dpk_)') goto 9999 end if endif - -!!$ if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_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,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 + + if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_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,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 ! ! ! Apply multiple sweeps of a block-Jacobi solver ! to compute an approximate solution of a linear system. ! ! - call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.) call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.) + + ! + ! 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 ty%zero() + + call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,info,init='Z') + case('Y') + call psb_geaxpby(zone,x,zzero,tx,desc_data,info) call psb_geaxpby(zone,y,zzero,ty,desc_data,info) + call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,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(zone,x,zzero,tx,desc_data,info) call psb_geaxpby(zone,initu,zzero,ty,desc_data,info) + call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,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 + + 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 @@ -188,22 +191,17 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& goto 9999 end if -!!$ else -!!$ -!!$ info = psb_err_iarg_neg_ -!!$ call psb_errpush(info,name,& -!!$ & i_err=(/itwo,sweeps,izero,izero,izero/)) -!!$ goto 9999 -!!$ -!!$ endif - - if (n_col <= size(work)) then - if ((4*n_col+n_col) <= size(work)) then - else - deallocate(aux) - endif else - deallocate(ww,aux) + + 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 call psb_erractionrestore(err_act)