From e4135c4b442580d961b6511cd00ec18edc738a52 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Fri, 4 Jan 2008 10:22:03 +0000 Subject: [PATCH] Fixed transpose (hopefully!) --- mlprec/mld_dbjac_aply.f90 | 58 +++-- mlprec/mld_dmlprec_aply.f90 | 100 ++++---- mlprec/mld_zbjac_aply.f90 | 71 ++++-- mlprec/mld_zmlprec_aply.f90 | 496 ++++++++++++++++++------------------ 4 files changed, 391 insertions(+), 334 deletions(-) diff --git a/mlprec/mld_dbjac_aply.f90 b/mlprec/mld_dbjac_aply.f90 index 1f2df5c9..bbe4da39 100644 --- a/mlprec/mld_dbjac_aply.f90 +++ b/mlprec/mld_dbjac_aply.f90 @@ -152,6 +152,7 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) real(kind(1.d0)), pointer :: ww(:), aux(:), tx(:),ty(:) integer :: ictxt,np,me,i, err_act character(len=20) :: name + character :: trans_ interface subroutine mld_dumf_solve(flag,m,x,b,n,ptr,info) @@ -169,7 +170,8 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ictxt=psb_cd_get_context(desc_data) call psb_info(ictxt, me, np) - select case(toupper(trans)) + trans_ = toupper(trans) + select case(trans_) case('N') case('T','C') case default @@ -217,19 +219,19 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! solve a system through ILU(k)/MILU(k)/ILU(k,t) (replicated matrix). ! - select case(toupper(trans)) + select case(trans_) case('N') call psb_spsm(done,prec%av(mld_l_pr_),x,dzero,ww,desc_data,info,& - & trans='N',unit='L',diag=prec%d,choice=psb_none_,work=aux) + & trans=trans_,unit='L',diag=prec%d,choice=psb_none_,work=aux) if (info == 0) call psb_spsm(alpha,prec%av(mld_u_pr_),ww,beta,y,desc_data,info,& - & trans='N',unit='U',choice=psb_none_, work=aux) + & trans=trans_,unit='U',choice=psb_none_, work=aux) case('T','C') call psb_spsm(done,prec%av(mld_u_pr_),x,dzero,ww,desc_data,info,& - & trans=toupper(trans),unit='L',diag=prec%d,choice=psb_none_,work=aux) + & trans=trans_,unit='L',diag=prec%d,choice=psb_none_,work=aux) if (info == 0) call psb_spsm(alpha,prec%av(mld_l_pr_),ww,beta,y,desc_data,info,& - & trans=toupper(trans),unit='U',choice=psb_none_,work=aux) + & trans=trans_,unit='U',choice=psb_none_,work=aux) case default call psb_errpush(4001,name,a_err='Invalid TRANS in ILU subsolve') goto 9999 @@ -245,7 +247,7 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ww(1:n_row) = x(1:n_row) - select case(toupper(trans)) + select case(trans_) case('N') call mld_dslu_solve(0,n_row,1,ww,n_row,prec%iprcparm(mld_slu_ptr_),info) case('T','C') @@ -265,7 +267,7 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ww(1:n_row) = x(1:n_row) - select case(toupper(trans)) + select case(trans_) case('N') call mld_dsludist_solve(0,n_row,1,ww,n_row,prec%iprcparm(mld_slud_ptr_),info) case('T','C') @@ -285,7 +287,7 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! to apply the LU factorization in both cases. ! - select case(toupper(trans)) + select case(trans_) case('N') call mld_dumf_solve(0,n_row,ww,x,n_row,prec%iprcparm(mld_umf_numptr_),info) case('T','C') @@ -331,16 +333,17 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) goto 9999 end if - tx = dzero - ty = dzero select case(prec%iprcparm(mld_sub_solve_)) case(mld_ilu_n_,mld_milu_n_,mld_ilu_t_) ! ! Use ILU(k)/MILU(k)/ILU(k,t) on the blocks. ! - select case(toupper(trans)) + select case(trans_) case('N') + + tx = dzero + ty = dzero do i=1, prec%iprcparm(mld_smooth_sweeps_) ! ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the @@ -353,14 +356,17 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) if (info /=0) exit call psb_spsm(done,prec%av(mld_l_pr_),ty,dzero,ww,& & prec%desc_data,info,& - & trans='N',unit='L',diag=prec%d,choice=psb_none_,work=aux) + & trans=trans_,unit='L',diag=prec%d,choice=psb_none_,work=aux) if (info /=0) exit call psb_spsm(done,prec%av(mld_u_pr_),ww,dzero,tx,& & prec%desc_data,info,& - & trans='N',unit='U',choice=psb_none_,work=aux) + & trans=trans_,unit='U',choice=psb_none_,work=aux) if (info /=0) exit end do + case('T','C') + tx = dzero + ty = dzero do i=1, prec%iprcparm(mld_smooth_sweeps_) ! ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the @@ -369,15 +375,15 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! ty(1:n_row) = x(1:n_row) call psb_spmm(-done,prec%av(mld_ap_nd_),tx,done,ty,& - & prec%desc_data,info,work=aux,trans=toupper(trans)) + & prec%desc_data,info,work=aux,trans=trans_) if (info /=0) exit call psb_spsm(done,prec%av(mld_u_pr_),ty,dzero,ww,& & prec%desc_data,info,& - & trans=toupper(trans),unit='L',diag=prec%d,choice=psb_none_,work=aux) + & trans=trans_,unit='L',diag=prec%d,choice=psb_none_,work=aux) if (info /=0) exit call psb_spsm(done,prec%av(mld_l_pr_),ww,dzero,tx,& & prec%desc_data,info,& - & trans=toupper(trans),unit='U',choice=psb_none_,work=aux) + & trans=trans_,unit='U',choice=psb_none_,work=aux) if (info /=0) exit end do @@ -400,8 +406,10 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! Use the LU factorization from SuperLU. ! - select case(toupper(trans)) + select case(trans_) case('N') + tx = dzero + ty = dzero do i=1, prec%iprcparm(mld_smooth_sweeps_) ! ! Compute Y(k+1) = D^(-1)*(X-ND*Y(k)), where D and ND are the @@ -417,7 +425,10 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) if(info /= 0) exit tx(1:n_row) = ty(1:n_row) end do + case('T','C') + tx = dzero + ty = dzero do i=1, prec%iprcparm(mld_smooth_sweeps_) ! ! Compute Y(k+1) = D^(-1)*(X-ND*Y(k)), where D and ND are the @@ -426,7 +437,7 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! ty(1:n_row) = x(1:n_row) call psb_spmm(-done,prec%av(mld_ap_nd_),tx,done,ty,& - & prec%desc_data,info,work=aux,trans=toupper(trans)) + & prec%desc_data,info,work=aux,trans=trans_) if(info /= 0) exit call mld_dslu_solve(1,n_row,1,ty,n_row,prec%iprcparm(mld_slu_ptr_),info) if(info /= 0) exit @@ -444,8 +455,10 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! Use the LU factorization from UMFPACK. ! - select case(toupper(trans)) + select case(trans_) case('N') + tx = dzero + ty = dzero do i=1, prec%iprcparm(mld_smooth_sweeps_) ! ! Compute Y(k+1) = D^(-1)*(X-ND*Y(k)), where D and ND are the @@ -462,7 +475,10 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) if (info /= 0) exit tx(1:n_row) = ww(1:n_row) end do + case('T','C') + tx = dzero + ty = dzero do i=1, prec%iprcparm(mld_smooth_sweeps_) ! ! Compute Y(k+1) = D^(-1)*(X-ND*Y(k)), where D and ND are the @@ -471,7 +487,7 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! ty(1:n_row) = x(1:n_row) call psb_spmm(-done,prec%av(mld_ap_nd_),tx,done,ty,& - & prec%desc_data,info,work=aux,trans=toupper(trans)) + & prec%desc_data,info,work=aux,trans=trans_) if (info /= 0) exit call mld_dumf_solve(1,n_row,ww,ty,n_row,& diff --git a/mlprec/mld_dmlprec_aply.f90 b/mlprec/mld_dmlprec_aply.f90 index 0d10e455..9d5d7c91 100644 --- a/mlprec/mld_dmlprec_aply.f90 +++ b/mlprec/mld_dmlprec_aply.f90 @@ -175,7 +175,7 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) real(kind(0.d0)),intent(in) :: alpha,beta real(kind(0.d0)),intent(in) :: x(:) real(kind(0.d0)),intent(inout) :: y(:) - character :: trans + character, intent(in) :: trans real(kind(0.d0)),target :: work(:) integer, intent(out) :: info @@ -185,6 +185,7 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) integer :: debug_level, debug_unit integer :: ismth, nlev, ilev, icm character(len=20) :: name + character :: trans_ name='mld_dmlprec_aply' info = 0 @@ -199,6 +200,8 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) & write(debug_unit,*) me,' ',trim(name),& & ' Entry ', size(baseprecv) + trans_ = toupper(trans) + select case(baseprecv(2)%iprcparm(mld_ml_type_)) case(mld_no_ml_) @@ -211,7 +214,7 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) case(mld_add_ml_) - call add_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) + call add_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans_,work,info) case(mld_mult_ml_) @@ -225,16 +228,35 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) select case(baseprecv(2)%iprcparm(mld_smooth_pos_)) case(mld_post_smooth_) + + select case (trans_) + case('N') + call mlt_post_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans_,work,info) + case('T','C') + call mlt_pre_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans_,work,info) + case default + info = 4001 + call psb_errpush(info,name,a_err='invalid trans') + goto 9999 + end select - call mlt_post_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) - + case(mld_pre_smooth_) - call mlt_pre_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) + select case (trans_) + case('N') + call mlt_pre_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans_,work,info) + case('T','C') + call mlt_post_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans_,work,info) + case default + info = 4001 + call psb_errpush(info,name,a_err='invalid trans') + goto 9999 + end select case(mld_twoside_smooth_) - call mlt_twoside_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) + call mlt_twoside_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans_,work,info) case default info = 4013 @@ -273,16 +295,17 @@ contains real(kind(0.d0)),intent(in) :: alpha,beta real(kind(0.d0)),intent(in) :: x(:) real(kind(0.d0)),intent(inout) :: y(:) - character :: trans + character, intent(in) :: trans real(kind(0.d0)),target :: work(:) integer, intent(out) :: info ! Local variables - integer :: n_row,n_col - integer :: ictxt,np,me,i, nr2l,nc2l,err_act - integer :: debug_level, debug_unit - integer :: ismth, nlev, ilev, icm - character(len=20) :: name + integer :: n_row,n_col + integer :: ictxt,np,me,i, nr2l,nc2l,err_act + integer :: debug_level, debug_unit + integer :: ismth, nlev, ilev, icm + character(len=20) :: name + type psb_mlprec_wrk_type real(kind(1.d0)), allocatable :: tx(:), ty(:), x2l(:), y2l(:) end type psb_mlprec_wrk_type @@ -483,6 +506,7 @@ contains call psb_errpush(4001,name,a_err='Error on final update') goto 9999 end if + deallocate(mlprec_wrk,stat=info) if (info /= 0) then call psb_errpush(4000,name) @@ -511,16 +535,17 @@ contains real(kind(0.d0)),intent(in) :: alpha,beta real(kind(0.d0)),intent(in) :: x(:) real(kind(0.d0)),intent(inout) :: y(:) - character :: trans + character, intent(in) :: trans real(kind(0.d0)),target :: work(:) integer, intent(out) :: info ! Local variables - integer :: n_row,n_col - integer :: ictxt,np,me,i, nr2l,nc2l,err_act - integer :: debug_level, debug_unit - integer :: ismth, nlev, ilev, icm - character(len=20) :: name + integer :: n_row,n_col + integer :: ictxt,np,me,i, nr2l,nc2l,err_act + integer :: debug_level, debug_unit + integer :: ismth, nlev, ilev, icm + character(len=20) :: name + type psb_mlprec_wrk_type real(kind(1.d0)), allocatable :: tx(:), ty(:), x2l(:), y2l(:) end type psb_mlprec_wrk_type @@ -759,11 +784,6 @@ contains goto 9999 end if - - - - - deallocate(mlprec_wrk,stat=info) if (info /= 0) then call psb_errpush(4000,name) @@ -780,11 +800,9 @@ contains return end if return - end subroutine mlt_pre_ml_aply - subroutine mlt_post_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) implicit none ! Arguments @@ -793,16 +811,17 @@ contains real(kind(0.d0)),intent(in) :: alpha,beta real(kind(0.d0)),intent(in) :: x(:) real(kind(0.d0)),intent(inout) :: y(:) - character :: trans + character, intent(in) :: trans real(kind(0.d0)),target :: work(:) integer, intent(out) :: info ! Local variables - integer :: n_row,n_col - integer :: ictxt,np,me,i, nr2l,nc2l,err_act - integer :: debug_level, debug_unit - integer :: ismth, nlev, ilev, icm - character(len=20) :: name + integer :: n_row,n_col + integer :: ictxt,np,me,i, nr2l,nc2l,err_act + integer :: debug_level, debug_unit + integer :: ismth, nlev, ilev, icm + character(len=20) :: name + type psb_mlprec_wrk_type real(kind(1.d0)), allocatable :: tx(:), ty(:), x2l(:), y2l(:) end type psb_mlprec_wrk_type @@ -1071,7 +1090,6 @@ contains return end if return - end subroutine mlt_post_ml_aply @@ -1083,16 +1101,17 @@ contains real(kind(0.d0)),intent(in) :: alpha,beta real(kind(0.d0)),intent(in) :: x(:) real(kind(0.d0)),intent(inout) :: y(:) - character :: trans + character, intent(in) :: trans real(kind(0.d0)),target :: work(:) integer, intent(out) :: info ! Local variables - integer :: n_row,n_col - integer :: ictxt,np,me,i, nr2l,nc2l,err_act - integer :: debug_level, debug_unit - integer :: ismth, nlev, ilev, icm - character(len=20) :: name + integer :: n_row,n_col + integer :: ictxt,np,me,i, nr2l,nc2l,err_act + integer :: debug_level, debug_unit + integer :: ismth, nlev, ilev, icm + character(len=20) :: name + type psb_mlprec_wrk_type real(kind(1.d0)), allocatable :: tx(:), ty(:), x2l(:), y2l(:) end type psb_mlprec_wrk_type @@ -1118,7 +1137,6 @@ contains goto 9999 end if - ! ! Pre- and post-smoothing (symmetrized) ! @@ -1374,11 +1392,7 @@ contains return end if return - end subroutine mlt_twoside_ml_aply - - - end subroutine mld_dmlprec_aply diff --git a/mlprec/mld_zbjac_aply.f90 b/mlprec/mld_zbjac_aply.f90 index b07f88c9..65d099e8 100644 --- a/mlprec/mld_zbjac_aply.f90 +++ b/mlprec/mld_zbjac_aply.f90 @@ -152,6 +152,7 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) complex(kind(1.d0)), pointer :: ww(:), aux(:), tx(:),ty(:) integer :: ictxt,np,me,i, err_act character(len=20) :: name + character :: trans_ interface subroutine mld_zumf_solve(flag,m,x,b,n,ptr,info) @@ -169,7 +170,8 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ictxt=psb_cd_get_context(desc_data) call psb_info(ictxt, me, np) - select case(toupper(trans)) + trans_ = toupper(trans) + select case(trans_) case('N') case('T','C') case default @@ -217,19 +219,19 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! solve a system through ILU(k)/MILU(k)/ILU(k,t) (replicated matrix). ! - select case(toupper(trans)) + select case(trans_) case('N') call psb_spsm(zone,prec%av(mld_l_pr_),x,zzero,ww,desc_data,info,& - & trans='N',unit='L',diag=prec%d,choice=psb_none_,work=aux) + & trans=trans_,unit='L',diag=prec%d,choice=psb_none_,work=aux) if (info == 0) call psb_spsm(alpha,prec%av(mld_u_pr_),ww,beta,y,desc_data,info,& - & trans='N',unit='U',choice=psb_none_, work=aux) + & trans=trans_,unit='U',choice=psb_none_, work=aux) case('T','C') call psb_spsm(zone,prec%av(mld_u_pr_),x,zzero,ww,desc_data,info,& - & trans=trans,unit='L',diag=prec%d,choice=psb_none_, work=aux) + & trans=trans_,unit='L',diag=prec%d,choice=psb_none_, work=aux) if(info ==0) call psb_spsm(alpha,prec%av(mld_l_pr_),ww,beta,y,desc_data,info,& - & trans=trans,unit='U',choice=psb_none_,work=aux) + & trans=trans_,unit='U',choice=psb_none_,work=aux) case default call psb_errpush(4001,name,a_err='Invalid TRANS in ILU subsolve') goto 9999 @@ -245,7 +247,7 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ww(1:n_row) = x(1:n_row) - select case(toupper(trans)) + select case(trans_) case('N') call mld_zslu_solve(0,n_row,1,ww,n_row,prec%iprcparm(mld_slu_ptr_),info) case('T') @@ -267,7 +269,7 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ww(1:n_row) = x(1:n_row) - select case(toupper(trans)) + select case(trans_) case('N') call mld_zsludist_solve(0,n_row,1,ww,n_row,prec%iprcparm(mld_slud_ptr_),info) case('T') @@ -289,7 +291,7 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! to apply the LU factorization in both cases. ! - select case(toupper(trans)) + select case(trans_) case('N') call mld_zumf_solve(0,n_row,ww,x,n_row,prec%iprcparm(mld_umf_numptr_),info) case('T') @@ -337,17 +339,16 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) goto 9999 end if - tx = zzero - ty = zzero - select case(prec%iprcparm(mld_sub_solve_)) case(mld_ilu_n_,mld_milu_n_,mld_ilu_t_) ! ! Use ILU(k)/MILU(k)/ILU(k,t) on the blocks. ! - select case(toupper(trans)) - case('N') + select case(trans_) + case('N') + tx = zzero + ty = zzero do i=1, prec%iprcparm(mld_smooth_sweeps_) ! ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the @@ -360,14 +361,17 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) if (info /=0) exit call psb_spsm(zone,prec%av(mld_l_pr_),ty,zzero,ww,& & prec%desc_data,info,& - & trans='N',unit='L',diag=prec%d,choice=psb_none_,work=aux) + & trans=trans_,unit='L',diag=prec%d,choice=psb_none_,work=aux) if (info /=0) exit call psb_spsm(zone,prec%av(mld_u_pr_),ww,zzero,tx,& & prec%desc_data,info,& - & trans='N',unit='U',choice=psb_none_,work=aux) + & trans=trans_,unit='U',choice=psb_none_,work=aux) if (info /=0) exit end do + case('T','C') + tx = zzero + ty = zzero do i=1, prec%iprcparm(mld_smooth_sweeps_) ! ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the @@ -376,17 +380,18 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! ty(1:n_row) = x(1:n_row) call psb_spmm(-zone,prec%av(mld_ap_nd_),tx,zone,ty,& - & prec%desc_data,info,work=aux,trans=toupper(trans)) + & prec%desc_data,info,work=aux,trans=trans_) if (info /=0) exit call psb_spsm(zone,prec%av(mld_u_pr_),ty,zzero,ww,& & prec%desc_data,info,& - & trans=toupper(trans),unit='L',diag=prec%d,choice=psb_none_,work=aux) + & trans=trans_,unit='L',diag=prec%d,choice=psb_none_,work=aux) if (info /=0) exit call psb_spsm(zone,prec%av(mld_l_pr_),ww,zzero,tx,& & prec%desc_data,info,& - & trans=toupper(trans),unit='U',choice=psb_none_,work=aux) + & trans=trans_,unit='U',choice=psb_none_,work=aux) if (info /=0) exit end do + case default call psb_errpush(4001,name,a_err='Invalid TRANS in ILU subsolve') goto 9999 @@ -406,8 +411,10 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! Use the LU factorization from SuperLU. ! - select case(toupper(trans)) + select case(trans_) case('N') + tx = zzero + ty = zzero do i=1, prec%iprcparm(mld_smooth_sweeps_) ! ! Compute Y(k+1) = D^(-1)*(X-ND*Y(k)), where D and ND are the @@ -423,7 +430,10 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) if (info /= 0) exit tx(1:n_row) = ty(1:n_row) end do + case('T') + tx = zzero + ty = zzero do i=1, prec%iprcparm(mld_smooth_sweeps_) ! ! Compute Y(k+1) = D^(-1)*(X-ND*Y(k)), where D and ND are the @@ -432,14 +442,17 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! ty(1:n_row) = x(1:n_row) call psb_spmm(-zone,prec%av(mld_ap_nd_),tx,zone,ty,& - & prec%desc_data,info,work=aux,trans=toupper(trans)) + & prec%desc_data,info,work=aux,trans=trans_) if (info /= 0) exit call mld_zslu_solve(1,n_row,1,ty,n_row,prec%iprcparm(mld_slu_ptr_),info) if (info /= 0) exit tx(1:n_row) = ty(1:n_row) end do + case('C') + tx = zzero + ty = zzero do i=1, prec%iprcparm(mld_smooth_sweeps_) ! ! Compute Y(k+1) = D^(-1)*(X-ND*Y(k)), where D and ND are the @@ -448,7 +461,7 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! ty(1:n_row) = x(1:n_row) call psb_spmm(-zone,prec%av(mld_ap_nd_),tx,zone,ty,& - & prec%desc_data,info,work=aux,trans=toupper(trans)) + & prec%desc_data,info,work=aux,trans=trans_) if (info /= 0) exit call mld_zslu_solve(2,n_row,1,ty,n_row,prec%iprcparm(mld_slu_ptr_),info) @@ -466,8 +479,10 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! Use the LU factorization from UMFPACK. ! - select case(toupper(trans)) + select case(trans_) case('N') + tx = zzero + ty = zzero do i=1, prec%iprcparm(mld_smooth_sweeps_) ! ! Compute Y(k+1) = D^(-1)*(X-ND*Y(k)), where D and ND are the @@ -484,7 +499,10 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) if (info /= 0) exit tx(1:n_row) = ww(1:n_row) end do + case('T') + tx = zzero + ty = zzero do i=1, prec%iprcparm(mld_smooth_sweeps_) ! ! Compute Y(k+1) = D^(-1)*(X-ND*Y(k)), where D and ND are the @@ -493,7 +511,7 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! ty(1:n_row) = x(1:n_row) call psb_spmm(-zone,prec%av(mld_ap_nd_),tx,zone,ty,& - & prec%desc_data,info,work=aux,trans=toupper(trans)) + & prec%desc_data,info,work=aux,trans=trans_) if (info /= 0) exit call mld_zumf_solve(1,n_row,ww,ty,n_row,& @@ -501,7 +519,10 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) if (info /= 0) exit tx(1:n_row) = ww(1:n_row) end do + case('C') + tx = zzero + ty = zzero do i=1, prec%iprcparm(mld_smooth_sweeps_) ! ! Compute Y(k+1) = D^(-1)*(X-ND*Y(k)), where D and ND are the @@ -510,7 +531,7 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! ty(1:n_row) = x(1:n_row) call psb_spmm(-zone,prec%av(mld_ap_nd_),tx,zone,ty,& - & prec%desc_data,info,work=aux,trans=toupper(trans)) + & prec%desc_data,info,work=aux,trans=trans_) if (info /= 0) exit call mld_zumf_solve(2,n_row,ww,ty,n_row,& diff --git a/mlprec/mld_zmlprec_aply.f90 b/mlprec/mld_zmlprec_aply.f90 index 6241a8aa..18787bf8 100644 --- a/mlprec/mld_zmlprec_aply.f90 +++ b/mlprec/mld_zmlprec_aply.f90 @@ -175,17 +175,17 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) complex(kind(1.d0)),intent(in) :: alpha,beta complex(kind(1.d0)),intent(in) :: x(:) complex(kind(1.d0)),intent(inout) :: y(:) - character :: trans + character, intent(in) :: trans complex(kind(1.d0)),target :: work(:) integer, intent(out) :: info ! Local variables - integer :: n_row,n_col - integer :: ictxt,np,me,i, nr2l,nc2l,err_act + integer :: n_row,n_col + integer :: ictxt,np,me,i, nr2l,nc2l,err_act integer :: debug_level, debug_unit integer :: ismth, nlev, ilev, icm character(len=20) :: name - + character :: trans_ name = 'mld_zmlprec_aply' info = 0 @@ -200,6 +200,7 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) & write(debug_unit,*) me,' ',trim(name),& & ' Entry ', size(baseprecv) + trans_ = toupper(trans) select case(baseprecv(2)%iprcparm(mld_ml_type_)) @@ -213,7 +214,7 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) case(mld_add_ml_) - call add_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) + call add_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans_,work,info) case(mld_mult_ml_) @@ -227,16 +228,35 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) select case(baseprecv(2)%iprcparm(mld_smooth_pos_)) case(mld_post_smooth_) + + select case (trans_) + case('N') + call mlt_post_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans_,work,info) + case('T','C') + call mlt_pre_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans_,work,info) + case default + info = 4001 + call psb_errpush(info,name,a_err='invalid trans') + goto 9999 + end select - call mlt_post_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) - + case(mld_pre_smooth_) - call mlt_pre_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) + select case (trans_) + case('N') + call mlt_pre_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans_,work,info) + case('T','C') + call mlt_post_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans_,work,info) + case default + info = 4001 + call psb_errpush(info,name,a_err='invalid trans') + goto 9999 + end select case(mld_twoside_smooth_) - call mlt_twoside_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) + call mlt_twoside_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans_,work,info) case default info = 4013 @@ -265,20 +285,17 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) end if return - contains subroutine add_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) - - implicit none - + implicit none ! Arguments type(psb_desc_type),intent(in) :: desc_data type(mld_zbaseprc_type), intent(in) :: baseprecv(:) complex(kind(1.d0)),intent(in) :: alpha,beta complex(kind(1.d0)),intent(in) :: x(:) complex(kind(1.d0)),intent(inout) :: y(:) - character :: trans + character, intent(in) :: trans complex(kind(1.d0)),target :: work(:) integer, intent(out) :: info @@ -314,7 +331,6 @@ contains goto 9999 end if - ! ! Additive multilevel ! @@ -492,7 +508,6 @@ contains goto 9999 end if - deallocate(mlprec_wrk,stat=info) if (info /= 0) then call psb_errpush(4000,name) @@ -509,20 +524,19 @@ contains return end if return - end subroutine add_ml_aply - - subroutine mlt_post_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) + end subroutine add_ml_aply - implicit none + subroutine mlt_pre_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) + implicit none ! Arguments type(psb_desc_type),intent(in) :: desc_data type(mld_zbaseprc_type), intent(in) :: baseprecv(:) complex(kind(1.d0)),intent(in) :: alpha,beta complex(kind(1.d0)),intent(in) :: x(:) complex(kind(1.d0)),intent(inout) :: y(:) - character :: trans + character, intent(in) :: trans complex(kind(1.d0)),target :: work(:) integer, intent(out) :: info @@ -538,7 +552,7 @@ contains end type psb_mlprec_wrk_type type(psb_mlprec_wrk_type), allocatable :: mlprec_wrk(:) - name = 'mlt_post_ml_aply' + name = 'mlt_pre_ml_aply' info = 0 call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() @@ -558,35 +572,42 @@ contains goto 9999 end if - ! - ! Post-smoothing - ! - ! 1. X(1) = Xext + ! Pre-smoothing ! - ! 2. DO ilev=2, nlev + ! 1. X(1) = Xext ! - ! ! Transfer X(ilev-1) to the next coarser level. - ! X(ilev) = AV(ilev; sm_pr_t_)*X(ilev-1) + ! 2. ! Apply the base preconditioner at the finest level. + ! Y(1) = (K(1)^(-1))*X(1) ! - ! ENDDO - ! - ! 3.! Apply the preconditioner at the coarsest level. - ! Y(nlev) = (K(nlev)^(-1))*X(nlev) + ! 3. ! Compute the residual at the finest level. + ! TX(1) = X(1) - A(1)*Y(1) ! - ! 4. DO ilev=nlev-1,1,-1 + ! 4. DO ilev=2, nlev ! - ! ! Transfer Y(ilev+1) to the next finer level. - ! Y(ilev) = AV(ilev+1; sm_pr_)*Y(ilev+1) + ! ! Transfer the residual to the current (coarser) level. + ! X(ilev) = AV(ilev; sm_pr_t_)*TX(ilev-1) ! - ! ! Compute the residual at the current level and apply to it the - ! ! base preconditioner. The sum over the subdomains is carried out - ! ! in the application of K(ilev). - ! Y(ilev) = Y(ilev) + (K(ilev)^(-1))*(X(ilev)-A(ilev)*Y(ilev)) + ! ! Apply the base preconditioner at the current level. + ! ! The sum over the subdomains is carried out in the + ! ! application of K(ilev). + ! Y(ilev) = (K(ilev)^(-1))*X(ilev) ! - ! ENDDO + ! ! Compute the residual at the current level (except at + ! ! the coarsest level). + ! IF (ilev < nlev) + ! TX(ilev) = (X(ilev)-A(ilev)*Y(ilev)) ! - ! 5. Yext = beta*Yext + alpha*Y(1) + ! ENDDO + ! + ! 5. DO ilev=nlev-1,1,-1 + ! + ! ! Transfer Y(ilev+1) to the next finer level + ! Y(ilev) = Y(ilev) + AV(ilev+1; sm_pr_)*Y(ilev+1) + ! + ! ENDDO + ! + ! 6. Yext = beta*Yext + alpha*Y(1) ! ! @@ -594,30 +615,55 @@ contains ! ! Copy the input vector X ! - if (debug_level >= psb_debug_inner_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' desc_data status',allocated(desc_data%matrix_data) - n_col = psb_cd_get_local_cols(desc_data) nc2l = psb_cd_get_local_cols(baseprecv(1)%desc_data) allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), & & mlprec_wrk(1)%tx(nc2l), stat=info) - mlprec_wrk(1)%x2l(:) = zzero - mlprec_wrk(1)%y2l(:) = zzero - mlprec_wrk(1)%tx(:) = zzero + if (info /= 0) then + info=4025 + call psb_errpush(info,name,i_err=(/4*nc2l,0,0,0,0/),& + & a_err='real(kind(1.d0))') + goto 9999 + end if - call psb_geaxpby(zone,x,zzero,mlprec_wrk(1)%tx,& - & baseprecv(1)%base_desc,info) - call psb_geaxpby(zone,x,zzero,mlprec_wrk(1)%x2l,& - & baseprecv(1)%base_desc,info) + mlprec_wrk(1)%y2l(:) = zzero + mlprec_wrk(1)%x2l(:) = x ! ! STEP 2 ! + ! Apply the base preconditioner at the finest level + ! + call mld_baseprec_aply(zone,baseprecv(1),mlprec_wrk(1)%x2l,& + & zzero,mlprec_wrk(1)%y2l,& + & baseprecv(1)%base_desc,& + & trans,work,info) + if (info /=0) then + call psb_errpush(4010,name,a_err=' baseprec_aply') + goto 9999 + end if + + ! + ! STEP 3 + ! + ! Compute the residual at the finest level + ! + mlprec_wrk(1)%tx = mlprec_wrk(1)%x2l + + call psb_spmm(-zone,baseprecv(1)%base_a,mlprec_wrk(1)%y2l,& + & zone,mlprec_wrk(1)%tx,baseprecv(1)%base_desc,info,work=work) + if (info /=0) then + call psb_errpush(4001,name,a_err=' fine level residual') + goto 9999 + end if + + ! + ! STEP 4 + ! ! For each level but the finest one ... ! - do ilev=2, nlev + do ilev = 2, nlev n_row = psb_cd_get_local_rows(baseprecv(ilev-1)%base_desc) n_col = psb_cd_get_local_cols(baseprecv(ilev-1)%desc_data) @@ -626,15 +672,8 @@ contains ismth = baseprecv(ilev)%iprcparm(mld_smooth_kind_) icm = baseprecv(ilev)%iprcparm(mld_coarse_mat_) - if (debug_level >= psb_debug_inner_) & - & write(debug_unit,*) me,' ',trim(name), & - & ' starting up sweep ',& - & ilev,allocated(baseprecv(ilev)%iprcparm),n_row,n_col,& - & nc2l, nr2l,ismth - allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%y2l(nc2l),& & mlprec_wrk(ilev)%x2l(nc2l), stat=info) - if (info /= 0) then info=4025 call psb_errpush(info,name,i_err=(/4*nc2l,0,0,0,0/),& @@ -645,36 +684,34 @@ contains mlprec_wrk(ilev)%x2l(:) = zzero mlprec_wrk(ilev)%y2l(:) = zzero mlprec_wrk(ilev)%tx(:) = zzero + if (ismth /= mld_no_smooth_) then ! ! Apply the smoothed prolongator transpose ! - if (debug_level >= psb_debug_inner_) & - & write(debug_unit,*) me,' ',trim(name), ' up sweep ', ilev - - call psb_halo(mlprec_wrk(ilev-1)%x2l,& - & baseprecv(ilev-1)%base_desc,info,work=work) + call psb_halo(mlprec_wrk(ilev-1)%tx,baseprecv(ilev-1)%base_desc,& + & info,work=work) if (info == 0) call psb_csmm(zone,baseprecv(ilev)%av(mld_sm_pr_t_),& - & mlprec_wrk(ilev-1)%x2l,zzero,mlprec_wrk(ilev)%x2l,info) + & mlprec_wrk(ilev-1)%tx,zzero,mlprec_wrk(ilev)%x2l,info) else ! ! Apply the raw aggregation map transpose (take a shortcut) ! + mlprec_wrk(ilev)%x2l = zzero do i=1,n_row mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) = & & mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) + & - & mlprec_wrk(ilev-1)%x2l(i) + & mlprec_wrk(ilev-1)%tx(i) end do - end if if (info /=0) then call psb_errpush(4001,name,a_err='Error during restriction') goto 9999 end if - if (icm == mld_repl_mat_) Then + if (icm ==mld_repl_mat_) then call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nr2l)) - else if (icm /= mld_distr_mat_) Then + else if (icm /= mld_distr_mat_) then info = 4013 call psb_errpush(info,name,a_err='invalid mld_coarse_mat_',& & i_Err=(/icm,0,0,0,0/)) @@ -682,64 +719,49 @@ contains endif ! - ! update x2l + ! Apply the base preconditioner ! - call psb_geaxpby(zone,mlprec_wrk(ilev)%x2l,zzero,mlprec_wrk(ilev)%tx,& - & baseprecv(ilev)%base_desc,info) - if (info /= 0) then - call psb_errpush(4001,name,a_err='Error in update') + call mld_baseprec_aply(zone,baseprecv(ilev),mlprec_wrk(ilev)%x2l,& + & zzero,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%desc_data, 'N',work,info) + + ! + ! Compute the residual (at all levels but the coarsest one) + ! + if (ilev < nlev) then + mlprec_wrk(ilev)%tx = mlprec_wrk(ilev)%x2l + if (info == 0) call psb_spmm(-zone,baseprecv(ilev)%base_a,& + & mlprec_wrk(ilev)%y2l,zone,mlprec_wrk(ilev)%tx,& + & baseprecv(ilev)%base_desc,info,work=work) + endif + if (info /=0) then + call psb_errpush(4001,name,a_err='Error on up sweep residual') goto 9999 end if - - if (debug_level >= psb_debug_inner_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' done up sweep ', ilev - enddo ! - ! STEP 3 - ! - ! Apply the base preconditioner at the coarsest level - ! - call mld_baseprec_aply(zone,baseprecv(nlev),mlprec_wrk(nlev)%x2l, & - & zzero, mlprec_wrk(nlev)%y2l,baseprecv(nlev)%desc_data,'N',work,info) - if (info /=0) then - call psb_errpush(4010,name,a_err='baseprec_aply') - goto 9999 - end if - - if (debug_level >= psb_debug_inner_) write(debug_unit,*) & - & me,' ',trim(name), ' done baseprec_aply ', nlev - - ! - ! STEP 4 + ! STEP 5 ! ! For each level but the coarsest one ... ! - do ilev=nlev-1, 1, -1 - - if (debug_level >= psb_debug_inner_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' starting down sweep',ilev + do ilev = nlev-1, 1, -1 ismth = baseprecv(ilev+1)%iprcparm(mld_smooth_kind_) n_row = psb_cd_get_local_rows(baseprecv(ilev)%base_desc) - if (ismth /= mld_no_smooth_) then + if (ismth /= mld_no_smooth_) then ! ! Apply the smoothed prolongator - ! + ! if (ismth == mld_smooth_prol_) & - & call psb_halo(mlprec_wrk(ilev+1)%y2l,baseprecv(ilev+1)%desc_data,& - & info,work=work) + & call psb_halo(mlprec_wrk(ilev+1)%y2l,& + & baseprecv(ilev+1)%desc_data,info,work=work) if (info == 0) call psb_csmm(zone,baseprecv(ilev+1)%av(mld_sm_pr_),& - & mlprec_wrk(ilev+1)%y2l, zzero,mlprec_wrk(ilev)%y2l,info) + & mlprec_wrk(ilev+1)%y2l,zone,mlprec_wrk(ilev)%y2l,info) else ! ! Apply the raw aggregation map (take a shortcut) ! - mlprec_wrk(ilev)%y2l(:) = zzero do i=1, n_row mlprec_wrk(ilev)%y2l(i) = mlprec_wrk(ilev)%y2l(i) + & & mlprec_wrk(ilev+1)%y2l(baseprecv(ilev+1)%mlia(i)) @@ -749,42 +771,20 @@ contains call psb_errpush(4001,name,a_err='Error during prolongation') goto 9999 end if - - ! - ! Compute the residual - ! - call psb_spmm(-zone,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,& - & zone,mlprec_wrk(ilev)%tx,baseprecv(ilev)%base_desc,info,work=work) - - ! - ! Apply the base preconditioner - ! - if (info == 0) call mld_baseprec_aply(zone,baseprecv(ilev),mlprec_wrk(ilev)%tx,& - & zone,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%base_desc, trans, work,info) - if (info /=0) then - call psb_errpush(4001,name,a_err=' spmm/baseprec_aply') - goto 9999 - end if - - if (debug_level >= psb_debug_inner_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' done down sweep',ilev enddo ! - ! STEP 5 + ! STEP 6 ! ! Compute the output vector Y ! - call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,baseprecv(1)%base_desc,info) - + call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,& + & baseprecv(1)%base_desc,info) if (info /=0) then - call psb_errpush(4001,name,a_err=' Final update') + call psb_errpush(4001,name,a_err='Error on final update') goto 9999 end if - - deallocate(mlprec_wrk,stat=info) if (info /= 0) then call psb_errpush(4000,name) @@ -801,21 +801,18 @@ contains return end if return - end subroutine mlt_post_ml_aply - - - - subroutine mlt_pre_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) + end subroutine mlt_pre_ml_aply - implicit none + subroutine mlt_post_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) + implicit none ! Arguments type(psb_desc_type),intent(in) :: desc_data type(mld_zbaseprc_type), intent(in) :: baseprecv(:) complex(kind(1.d0)),intent(in) :: alpha,beta complex(kind(1.d0)),intent(in) :: x(:) complex(kind(1.d0)),intent(inout) :: y(:) - character :: trans + character, intent(in) :: trans complex(kind(1.d0)),target :: work(:) integer, intent(out) :: info @@ -831,7 +828,7 @@ contains end type psb_mlprec_wrk_type type(psb_mlprec_wrk_type), allocatable :: mlprec_wrk(:) - name = 'mlt_pre_ml_aply' + name = 'mlt_post_ml_aply' info = 0 call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() @@ -852,41 +849,33 @@ contains end if ! - ! Pre-smoothing - ! - ! 1. X(1) = Xext - ! - ! 2. ! Apply the base preconditioner at the finest level. - ! Y(1) = (K(1)^(-1))*X(1) - ! - ! 3. ! Compute the residual at the finest level. - ! TX(1) = X(1) - A(1)*Y(1) + ! Post-smoothing + ! + ! 1. X(1) = Xext ! - ! 4. DO ilev=2, nlev + ! 2. DO ilev=2, nlev ! - ! ! Transfer the residual to the current (coarser) level. - ! X(ilev) = AV(ilev; sm_pr_t_)*TX(ilev-1) + ! ! Transfer X(ilev-1) to the next coarser level. + ! X(ilev) = AV(ilev; sm_pr_t_)*X(ilev-1) ! - ! ! Apply the base preconditioner at the current level. - ! ! The sum over the subdomains is carried out in the - ! ! application of K(ilev). - ! Y(ilev) = (K(ilev)^(-1))*X(ilev) + ! ENDDO + ! + ! 3.! Apply the preconditioner at the coarsest level. + ! Y(nlev) = (K(nlev)^(-1))*X(nlev) ! - ! ! Compute the residual at the current level (except at - ! ! the coarsest level). - ! IF (ilev < nlev) - ! TX(ilev) = (X(ilev)-A(ilev)*Y(ilev)) + ! 4. DO ilev=nlev-1,1,-1 ! - ! ENDDO + ! ! Transfer Y(ilev+1) to the next finer level. + ! Y(ilev) = AV(ilev+1; sm_pr_)*Y(ilev+1) ! - ! 5. DO ilev=nlev-1,1,-1 + ! ! Compute the residual at the current level and apply to it the + ! ! base preconditioner. The sum over the subdomains is carried out + ! ! in the application of K(ilev). + ! Y(ilev) = Y(ilev) + (K(ilev)^(-1))*(X(ilev)-A(ilev)*Y(ilev)) ! - ! ! Transfer Y(ilev+1) to the next finer level - ! Y(ilev) = Y(ilev) + AV(ilev+1; sm_pr_)*Y(ilev+1) + ! ENDDO ! - ! ENDDO - ! - ! 6. Yext = beta*Yext + alpha*Y(1) + ! 5. Yext = beta*Yext + alpha*Y(1) ! ! @@ -894,55 +883,30 @@ contains ! ! Copy the input vector X ! + if (debug_level >= psb_debug_inner_) & + & write(debug_unit,*) me,' ',trim(name),& + & ' desc_data status',allocated(desc_data%matrix_data) + n_col = psb_cd_get_local_cols(desc_data) nc2l = psb_cd_get_local_cols(baseprecv(1)%desc_data) allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), & & mlprec_wrk(1)%tx(nc2l), stat=info) - if (info /= 0) then - info=4025 - call psb_errpush(info,name,i_err=(/4*nc2l,0,0,0,0/),& - & a_err='real(kind(1.d0))') - goto 9999 - end if - + mlprec_wrk(1)%x2l(:) = zzero mlprec_wrk(1)%y2l(:) = zzero - mlprec_wrk(1)%x2l(:) = x - - ! - ! STEP 2 - ! - ! Apply the base preconditioner at the finest level - ! - call mld_baseprec_aply(zone,baseprecv(1),mlprec_wrk(1)%x2l,& - & zzero,mlprec_wrk(1)%y2l,& - & baseprecv(1)%base_desc,& - & trans,work,info) - if (info /=0) then - call psb_errpush(4010,name,a_err=' baseprec_aply') - goto 9999 - end if - - ! - ! STEP 3 - ! - ! Compute the residual at the finest level - ! - mlprec_wrk(1)%tx = mlprec_wrk(1)%x2l + mlprec_wrk(1)%tx(:) = zzero - call psb_spmm(-zone,baseprecv(1)%base_a,mlprec_wrk(1)%y2l,& - & zone,mlprec_wrk(1)%tx,baseprecv(1)%base_desc,info,work=work) - if (info /=0) then - call psb_errpush(4001,name,a_err=' fine level residual') - goto 9999 - end if + call psb_geaxpby(zone,x,zzero,mlprec_wrk(1)%tx,& + & baseprecv(1)%base_desc,info) + call psb_geaxpby(zone,x,zzero,mlprec_wrk(1)%x2l,& + & baseprecv(1)%base_desc,info) ! - ! STEP 4 + ! STEP 2 ! ! For each level but the finest one ... ! - do ilev = 2, nlev + do ilev=2, nlev n_row = psb_cd_get_local_rows(baseprecv(ilev-1)%base_desc) n_col = psb_cd_get_local_cols(baseprecv(ilev-1)%desc_data) @@ -951,8 +915,15 @@ contains ismth = baseprecv(ilev)%iprcparm(mld_smooth_kind_) icm = baseprecv(ilev)%iprcparm(mld_coarse_mat_) + if (debug_level >= psb_debug_inner_) & + & write(debug_unit,*) me,' ',trim(name), & + & ' starting up sweep ',& + & ilev,allocated(baseprecv(ilev)%iprcparm),n_row,n_col,& + & nc2l, nr2l,ismth + allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%y2l(nc2l),& & mlprec_wrk(ilev)%x2l(nc2l), stat=info) + if (info /= 0) then info=4025 call psb_errpush(info,name,i_err=(/4*nc2l,0,0,0,0/),& @@ -963,34 +934,36 @@ contains mlprec_wrk(ilev)%x2l(:) = zzero mlprec_wrk(ilev)%y2l(:) = zzero mlprec_wrk(ilev)%tx(:) = zzero - if (ismth /= mld_no_smooth_) then ! ! Apply the smoothed prolongator transpose ! - call psb_halo(mlprec_wrk(ilev-1)%tx,baseprecv(ilev-1)%base_desc,& - & info,work=work) + if (debug_level >= psb_debug_inner_) & + & write(debug_unit,*) me,' ',trim(name), ' up sweep ', ilev + + call psb_halo(mlprec_wrk(ilev-1)%x2l,& + & baseprecv(ilev-1)%base_desc,info,work=work) if (info == 0) call psb_csmm(zone,baseprecv(ilev)%av(mld_sm_pr_t_),& - & mlprec_wrk(ilev-1)%tx,zzero,mlprec_wrk(ilev)%x2l,info) + & mlprec_wrk(ilev-1)%x2l,zzero,mlprec_wrk(ilev)%x2l,info) else ! ! Apply the raw aggregation map transpose (take a shortcut) ! - mlprec_wrk(ilev)%x2l = zzero do i=1,n_row mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) = & & mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) + & - & mlprec_wrk(ilev-1)%tx(i) + & mlprec_wrk(ilev-1)%x2l(i) end do + end if if (info /=0) then call psb_errpush(4001,name,a_err='Error during restriction') goto 9999 end if - if (icm ==mld_repl_mat_) then + if (icm == mld_repl_mat_) Then call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nr2l)) - else if (icm /= mld_distr_mat_) then + else if (icm /= mld_distr_mat_) Then info = 4013 call psb_errpush(info,name,a_err='invalid mld_coarse_mat_',& & i_Err=(/icm,0,0,0,0/)) @@ -998,49 +971,64 @@ contains endif ! - ! Apply the base preconditioner - ! - call mld_baseprec_aply(zone,baseprecv(ilev),mlprec_wrk(ilev)%x2l,& - & zzero,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%desc_data, 'N',work,info) - - ! - ! Compute the residual (at all levels but the coarsest one) + ! update x2l ! - if (ilev < nlev) then - mlprec_wrk(ilev)%tx = mlprec_wrk(ilev)%x2l - if (info == 0) call psb_spmm(-zone,baseprecv(ilev)%base_a,& - & mlprec_wrk(ilev)%y2l,zone,mlprec_wrk(ilev)%tx,& - & baseprecv(ilev)%base_desc,info,work=work) - endif - if (info /=0) then - call psb_errpush(4001,name,a_err='Error on up sweep residual') + call psb_geaxpby(zone,mlprec_wrk(ilev)%x2l,zzero,mlprec_wrk(ilev)%tx,& + & baseprecv(ilev)%base_desc,info) + if (info /= 0) then + call psb_errpush(4001,name,a_err='Error in update') goto 9999 end if + + if (debug_level >= psb_debug_inner_) & + & write(debug_unit,*) me,' ',trim(name),& + & ' done up sweep ', ilev + enddo ! - ! STEP 5 + ! STEP 3 + ! + ! Apply the base preconditioner at the coarsest level + ! + call mld_baseprec_aply(zone,baseprecv(nlev),mlprec_wrk(nlev)%x2l, & + & zzero, mlprec_wrk(nlev)%y2l,baseprecv(nlev)%desc_data,'N',work,info) + if (info /=0) then + call psb_errpush(4010,name,a_err='baseprec_aply') + goto 9999 + end if + + if (debug_level >= psb_debug_inner_) write(debug_unit,*) & + & me,' ',trim(name), ' done baseprec_aply ', nlev + + ! + ! STEP 4 ! ! For each level but the coarsest one ... ! - do ilev = nlev-1, 1, -1 + do ilev=nlev-1, 1, -1 + + if (debug_level >= psb_debug_inner_) & + & write(debug_unit,*) me,' ',trim(name),& + & ' starting down sweep',ilev ismth = baseprecv(ilev+1)%iprcparm(mld_smooth_kind_) n_row = psb_cd_get_local_rows(baseprecv(ilev)%base_desc) - if (ismth /= mld_no_smooth_) then + if (ismth /= mld_no_smooth_) then ! ! Apply the smoothed prolongator - ! + ! if (ismth == mld_smooth_prol_) & - & call psb_halo(mlprec_wrk(ilev+1)%y2l,& - & baseprecv(ilev+1)%desc_data,info,work=work) + & call psb_halo(mlprec_wrk(ilev+1)%y2l,baseprecv(ilev+1)%desc_data,& + & info,work=work) if (info == 0) call psb_csmm(zone,baseprecv(ilev+1)%av(mld_sm_pr_),& - & mlprec_wrk(ilev+1)%y2l,zone,mlprec_wrk(ilev)%y2l,info) + & mlprec_wrk(ilev+1)%y2l, zzero,mlprec_wrk(ilev)%y2l,info) else ! ! Apply the raw aggregation map (take a shortcut) ! + mlprec_wrk(ilev)%y2l(:) = zzero do i=1, n_row mlprec_wrk(ilev)%y2l(i) = mlprec_wrk(ilev)%y2l(i) + & & mlprec_wrk(ilev+1)%y2l(baseprecv(ilev+1)%mlia(i)) @@ -1050,17 +1038,37 @@ contains call psb_errpush(4001,name,a_err='Error during prolongation') goto 9999 end if + + ! + ! Compute the residual + ! + call psb_spmm(-zone,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,& + & zone,mlprec_wrk(ilev)%tx,baseprecv(ilev)%base_desc,info,work=work) + + ! + ! Apply the base preconditioner + ! + if (info == 0) call mld_baseprec_aply(zone,baseprecv(ilev),mlprec_wrk(ilev)%tx,& + & zone,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%base_desc, trans, work,info) + if (info /=0) then + call psb_errpush(4001,name,a_err=' spmm/baseprec_aply') + goto 9999 + end if + + if (debug_level >= psb_debug_inner_) & + & write(debug_unit,*) me,' ',trim(name),& + & ' done down sweep',ilev enddo ! - ! STEP 6 + ! STEP 5 ! ! Compute the output vector Y ! - call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,& - & baseprecv(1)%base_desc,info) + call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,baseprecv(1)%base_desc,info) + if (info /=0) then - call psb_errpush(4001,name,a_err='Error on final update') + call psb_errpush(4001,name,a_err=' Final update') goto 9999 end if @@ -1082,20 +1090,18 @@ contains return end if return - end subroutine mlt_pre_ml_aply + end subroutine mlt_post_ml_aply subroutine mlt_twoside_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) - - implicit none - + implicit none ! Arguments type(psb_desc_type),intent(in) :: desc_data type(mld_zbaseprc_type), intent(in) :: baseprecv(:) complex(kind(1.d0)),intent(in) :: alpha,beta complex(kind(1.d0)),intent(in) :: x(:) complex(kind(1.d0)),intent(inout) :: y(:) - character :: trans + character, intent(in) :: trans complex(kind(1.d0)),target :: work(:) integer, intent(out) :: info