fdurastante
/
amg4psblas
mirror of https://github.com/sfilippone/amg4psblas.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

mld2p4:

Browse Source 
mlprec/mld_cmlprec_aply.f90
 mlprec/mld_cprecaply.f90
 mlprec/mld_dmlprec_aply.f90
 mlprec/mld_dprecaply.f90
 mlprec/mld_inner_mod.f90
 mlprec/mld_smlprec_aply.f90
 mlprec/mld_sprecaply.f90
 mlprec/mld_zmlprec_aply.f90
 mlprec/mld_zprecaply.f90

Fixed internal interfaces in mlprec_aply.
stopcriterion
Salvatore Filippone 16 years ago
parent 08dd566600
commit 7681de7a2e
9 changed files with 548 additions and 543 deletions
Show Stats Download Patch File Download Diff File

224
mlprec/mld_cmlprec_aply.f90
Unescape Escape View File

@ -78,41 +78,41 @@
! Arguments: ! Arguments:
! alpha - complex(psb_spk_), input. ! alpha - complex(psb_spk_), input.
! The scalar alpha. ! The scalar alpha.
! precv - type(mld_conelev_type), dimension(:), input. ! p - type(mld_cprec_type), input.
! The array of one-level preconditioner data structures containing the ! The array of one-level preconditioner data structures containing the
! local parts of the preconditioners to be applied at each level. ! local parts of the preconditioners to be applied at each level.
! Note that nlev = size(precv) = number of levels. ! Note that nlev = size(p%precv) = number of levels.
! precv(ilev)%prec - type(psb_cbaseprec_type) ! p%precv(ilev)%prec - type(psb_cbaseprec_type)
! The "base" preconditioner for the current level ! The "base" preconditioner for the current level
! precv(ilev)%ac - type(psb_cspmat_type) ! p%precv(ilev)%ac - type(psb_cspmat_type)
! The local part of the matrix A(ilev). ! The local part of the matrix A(ilev).
! precv(ilev)%desc_ac - type(psb_desc_type). ! p%precv(ilev)%desc_ac - type(psb_desc_type).
! The communication descriptor associated to the sparse ! The communication descriptor associated to the sparse
! matrix A(ilev) ! matrix A(ilev)
! precv(ilev)%map - type(psb_inter_desc_type) ! p%precv(ilev)%map - type(psb_inter_desc_type)
! Stores the linear operators mapping between levels ! Stores the linear operators mapping between levels
! (ilev-1) and (ilev). These are the restriction and ! (ilev-1) and (ilev). These are the restriction and
! prolongation operators described in the sequel. ! prolongation operators described in the sequel.
! precv(ilev)%iprcparm - integer, dimension(:), allocatable. ! p%precv(ilev)%iprcparm - integer, dimension(:), allocatable.
! The integer parameters defining the multilevel ! The integer parameters defining the multilevel
! strategy ! strategy
! precv(ilev)%rprcparm - real(psb_spk_), dimension(:), allocatable. ! p%precv(ilev)%rprcparm - real(psb_spk_), dimension(:), allocatable.
! The real parameters defining the multilevel strategy ! The real parameters defining the multilevel strategy
! precv(ilev)%mlia - integer, dimension(:), allocatable. ! p%precv(ilev)%mlia - integer, dimension(:), allocatable.
! The aggregation map (ilev-1) --> (ilev). ! The aggregation map (ilev-1) --> (ilev).
! In case of non-smoothed aggregation, it is used ! In case of non-smoothed aggregation, it is used
! instead of mld_sm_pr_. ! instead of mld_sm_pr_.
! precv(ilev)%nlaggr - integer, dimension(:), allocatable. ! p%precv(ilev)%nlaggr - integer, dimension(:), allocatable.
! The number of aggregates (rows of A(ilev)) on the ! The number of aggregates (rows of A(ilev)) on the
! various processes. ! various processes.
! precv(ilev)%base_a - type(psb_cspmat_type), pointer. ! p%precv(ilev)%base_a - type(psb_cspmat_type), pointer.
! Pointer (really a pointer!) to the base matrix of ! Pointer (really a pointer!) to the base matrix of
! the current level, i.e. the local part of A(ilev); ! the current level, i.e. the local part of A(ilev);
! so we have a unified treatment of residuals. We ! so we have a unified treatment of residuals. We
! need this to avoid passing explicitly the matrix ! need this to avoid passing explicitly the matrix
! A(ilev) to the routine which applies the ! A(ilev) to the routine which applies the
! preconditioner. ! preconditioner.
! precv(ilev)%base_desc - type(psb_desc_type), pointer. ! p%precv(ilev)%base_desc - type(psb_desc_type), pointer.
! Pointer to the communication descriptor associated ! Pointer to the communication descriptor associated
! to the sparse matrix pointed by base_a. ! to the sparse matrix pointed by base_a.
! !
@ -136,10 +136,10 @@
! Note that when the LU factorization of the matrix A(ilev) is computed instead of ! Note that when the LU factorization of the matrix A(ilev) is computed instead of
! the ILU one, by using UMFPACK or SuperLU, the corresponding L and U factors ! the ILU one, by using UMFPACK or SuperLU, the corresponding L and U factors
! are stored in data structures provided by UMFPACK or SuperLU and pointed by ! are stored in data structures provided by UMFPACK or SuperLU and pointed by
! precv(ilev)%prec%iprcparm(mld_umf_ptr) or precv(ilev)%prec%iprcparm(mld_slu_ptr), ! p%precv(ilev)%prec%iprcparm(mld_umf_ptr) or p%precv(ilev)%prec%iprcparm(mld_slu_ptr),
! respectively. ! respectively.
! !
subroutine mld_cmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info) subroutine mld_cmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
use psb_base_mod use psb_base_mod
use mld_inner_mod, mld_protect_name => mld_cmlprec_aply use mld_inner_mod, mld_protect_name => mld_cmlprec_aply
@ -148,7 +148,7 @@ subroutine mld_cmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
! Arguments ! Arguments
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
type(mld_conelev_type), intent(in) :: precv(:) type(mld_cprec_type), intent(in) :: p
complex(psb_spk_),intent(in) :: alpha,beta complex(psb_spk_),intent(in) :: alpha,beta
complex(psb_spk_),intent(in) :: x(:) complex(psb_spk_),intent(in) :: x(:)
complex(psb_spk_),intent(inout) :: y(:) complex(psb_spk_),intent(inout) :: y(:)
@ -173,11 +173,11 @@ subroutine mld_cmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(precv) & ' Entry ', size(p%precv)
trans_ = psb_toupper(trans) trans_ = psb_toupper(trans)
select case(precv(2)%iprcparm(mld_ml_type_)) select case(p%precv(2)%iprcparm(mld_ml_type_))
case(mld_no_ml_) case(mld_no_ml_)
! !
@ -191,7 +191,7 @@ subroutine mld_cmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
! Additive multilevel ! Additive multilevel
! !
call add_ml_aply(alpha,precv,x,beta,y,desc_data,trans_,work,info) call add_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info)
case(mld_mult_ml_) case(mld_mult_ml_)
! !
@ -202,15 +202,15 @@ subroutine mld_cmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
! Note that the transpose switches pre <-> post. ! Note that the transpose switches pre <-> post.
! !
select case(precv(2)%iprcparm(mld_smoother_pos_)) select case(p%precv(2)%iprcparm(mld_smoother_pos_))
case(mld_post_smooth_) case(mld_post_smooth_)
select case (trans_) select case (trans_)
case('N') case('N')
call mlt_post_ml_aply(alpha,precv,x,beta,y,desc_data,trans_,work,info) call mlt_post_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info)
case('T','C') case('T','C')
call mlt_pre_ml_aply(alpha,precv,x,beta,y,desc_data,trans_,work,info) call mlt_pre_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info)
case default case default
info = 4001 info = 4001
call psb_errpush(info,name,a_err='invalid trans') call psb_errpush(info,name,a_err='invalid trans')
@ -221,9 +221,9 @@ subroutine mld_cmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
select case (trans_) select case (trans_)
case('N') case('N')
call mlt_pre_ml_aply(alpha,precv,x,beta,y,desc_data,trans_,work,info) call mlt_pre_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info)
case('T','C') case('T','C')
call mlt_post_ml_aply(alpha,precv,x,beta,y,desc_data,trans_,work,info) call mlt_post_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info)
case default case default
info = 4001 info = 4001
call psb_errpush(info,name,a_err='invalid trans') call psb_errpush(info,name,a_err='invalid trans')
@ -232,12 +232,12 @@ subroutine mld_cmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
case(mld_twoside_smooth_) case(mld_twoside_smooth_)
call mlt_twoside_ml_aply(alpha,precv,x,beta,y,desc_data,trans_,work,info) call mlt_twoside_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info)
case default case default
info = 4013 info = 4013
call psb_errpush(info,name,a_err='invalid smooth_pos',& call psb_errpush(info,name,a_err='invalid smooth_pos',&
& i_Err=(/precv(2)%iprcparm(mld_smoother_pos_),0,0,0,0/)) & i_Err=(/p%precv(2)%iprcparm(mld_smoother_pos_),0,0,0,0/))
goto 9999 goto 9999
end select end select
@ -245,7 +245,7 @@ subroutine mld_cmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
case default case default
info = 4013 info = 4013
call psb_errpush(info,name,a_err='invalid mltype',& call psb_errpush(info,name,a_err='invalid mltype',&
& i_Err=(/precv(2)%iprcparm(mld_ml_type_),0,0,0,0/)) & i_Err=(/p%precv(2)%iprcparm(mld_ml_type_),0,0,0,0/))
goto 9999 goto 9999
end select end select
@ -272,7 +272,7 @@ contains
! Y = beta*Y + alpha*op(M^(-1))*X, ! Y = beta*Y + alpha*op(M^(-1))*X,
! where ! where
! - M is an additive multilevel domain decomposition (Schwarz) preconditioner ! - M is an additive multilevel domain decomposition (Schwarz) preconditioner
! associated to a certain matrix A and stored in the array precv, ! associated to a certain matrix A and stored in the array p%precv,
! - op(M^(-1)) is M^(-1) or its (conjugate) transpose, according to ! - op(M^(-1)) is M^(-1) or its (conjugate) transpose, according to
! the value of trans, ! the value of trans,
! - X and Y are vectors, ! - X and Y are vectors,
@ -287,7 +287,7 @@ contains
! !
! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners', ! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners',
! each representing the part of the preconditioner associated to a certain level. ! each representing the part of the preconditioner associated to a certain level.
! For each level ilev, the base preconditioner K(ilev) is stored in precv(ilev) ! For each level ilev, the base preconditioner K(ilev) is stored in p%precv(ilev)
! and is associated to a matrix A(ilev), obtained by 'tranferring' the original ! and is associated to a matrix A(ilev), obtained by 'tranferring' the original
! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed ! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed
! aggregation. ! aggregation.
@ -335,13 +335,13 @@ contains
! !
! 4. Yext = beta*Yext + alpha*Y(1) ! 4. Yext = beta*Yext + alpha*Y(1)
! !
subroutine add_ml_aply(alpha,precv,x,beta,y,desc_data,trans,work,info) subroutine add_ml_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
implicit none implicit none
! Arguments ! Arguments
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
type(mld_conelev_type), intent(in) :: precv(:) type(mld_cprec_type), intent(in) :: p
complex(psb_spk_),intent(in) :: alpha,beta complex(psb_spk_),intent(in) :: alpha,beta
complex(psb_spk_),intent(in) :: x(:) complex(psb_spk_),intent(in) :: x(:)
complex(psb_spk_),intent(inout) :: y(:) complex(psb_spk_),intent(inout) :: y(:)
@ -371,9 +371,9 @@ contains
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(precv) & ' Entry ', size(p%precv)
nlev = size(precv) nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info) allocate(mlprec_wrk(nlev),stat=info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
@ -396,8 +396,8 @@ contains
mlprec_wrk(1)%x2l(:) = x(:) mlprec_wrk(1)%x2l(:) = x(:)
mlprec_wrk(1)%y2l(:) = czero mlprec_wrk(1)%y2l(:) = czero
call mld_baseprec_aply(alpha,precv(1)%prec,x,beta,y,& call mld_baseprec_aply(alpha,p%precv(1)%prec,x,beta,y,&
& precv(1)%base_desc,trans,work,info) & p%precv(1)%base_desc,trans,work,info)
if (info /=0) then if (info /=0) then
call psb_errpush(4010,name,a_err='baseprec_aply') call psb_errpush(4010,name,a_err='baseprec_aply')
goto 9999 goto 9999
@ -408,8 +408,8 @@ contains
! For each level except the finest one ... ! For each level except the finest one ...
! !
do ilev = 2, nlev do ilev = 2, nlev
nc2l = psb_cd_get_local_cols(precv(ilev)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(ilev)%base_desc)
nr2l = psb_cd_get_local_rows(precv(ilev)%base_desc) nr2l = psb_cd_get_local_rows(p%precv(ilev)%base_desc)
allocate(mlprec_wrk(ilev)%x2l(nc2l),mlprec_wrk(ilev)%y2l(nc2l),& allocate(mlprec_wrk(ilev)%x2l(nc2l),mlprec_wrk(ilev)%y2l(nc2l),&
& stat=info) & stat=info)
if (info /= 0) then if (info /= 0) then
@ -422,7 +422,7 @@ contains
! Apply prolongator transpose, i.e. restriction ! Apply prolongator transpose, i.e. restriction
call psb_map_X2Y(cone,mlprec_wrk(ilev-1)%x2l,& call psb_map_X2Y(cone,mlprec_wrk(ilev-1)%x2l,&
& czero,mlprec_wrk(ilev)%x2l,& & czero,mlprec_wrk(ilev)%x2l,&
& precv(ilev)%map,info,work=work) & p%precv(ilev)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction') call psb_errpush(4001,name,a_err='Error during restriction')
@ -432,9 +432,9 @@ contains
! !
! Apply the base preconditioner ! Apply the base preconditioner
! !
call mld_baseprec_aply(cone,precv(ilev)%prec,& call mld_baseprec_aply(cone,p%precv(ilev)%prec,&
& mlprec_wrk(ilev)%x2l,czero,mlprec_wrk(ilev)%y2l,& & mlprec_wrk(ilev)%x2l,czero,mlprec_wrk(ilev)%y2l,&
& precv(ilev)%base_desc,trans,work,info) & p%precv(ilev)%base_desc,trans,work,info)
enddo enddo
@ -445,15 +445,15 @@ contains
! !
do ilev =nlev,2,-1 do ilev =nlev,2,-1
nc2l = psb_cd_get_local_cols(precv(ilev)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(ilev)%base_desc)
nr2l = psb_cd_get_local_rows(precv(ilev)%base_desc) nr2l = psb_cd_get_local_rows(p%precv(ilev)%base_desc)
! !
! Apply prolongator ! Apply prolongator
! !
call psb_map_Y2X(cone,mlprec_wrk(ilev)%y2l,& call psb_map_Y2X(cone,mlprec_wrk(ilev)%y2l,&
& cone,mlprec_wrk(ilev-1)%y2l,& & cone,mlprec_wrk(ilev-1)%y2l,&
& precv(ilev)%map,info,work=work) & p%precv(ilev)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during prolongation') call psb_errpush(4001,name,a_err='Error during prolongation')
@ -466,7 +466,7 @@ contains
! !
! Compute the output vector Y ! Compute the output vector Y
! !
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,cone,y,precv(1)%base_desc,info) call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,cone,y,p%precv(1)%base_desc,info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4001,name,a_err='Error on final update') call psb_errpush(4001,name,a_err='Error on final update')
goto 9999 goto 9999
@ -500,7 +500,7 @@ contains
! Y = beta*Y + alpha*op(M^(-1))*X, ! Y = beta*Y + alpha*op(M^(-1))*X,
! where ! where
! - M is a hybrid multilevel domain decomposition (Schwarz) preconditioner ! - M is a hybrid multilevel domain decomposition (Schwarz) preconditioner
! associated to a certain matrix A and stored in the array precv, ! associated to a certain matrix A and stored in the array p%precv,
! - op(M^(-1)) is M^(-1) or its (conjugate) transpose, according to ! - op(M^(-1)) is M^(-1) or its (conjugate) transpose, according to
! the value of trans, ! the value of trans,
! - X and Y are vectors, ! - X and Y are vectors,
@ -515,7 +515,7 @@ contains
! !
! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners', ! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners',
! each representing the part of the preconditioner associated to a certain level. ! each representing the part of the preconditioner associated to a certain level.
! For each level ilev, the base preconditioner K(ilev) is stored in precv(ilev) ! For each level ilev, the base preconditioner K(ilev) is stored in p%precv(ilev)
! and is associated to a matrix A(ilev), obtained by 'tranferring' the original ! and is associated to a matrix A(ilev), obtained by 'tranferring' the original
! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed ! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed
! aggregation. ! aggregation.
@ -571,13 +571,13 @@ contains
! 6. Yext = beta*Yext + alpha*Y(1) ! 6. Yext = beta*Yext + alpha*Y(1)
! !
! !
subroutine mlt_pre_ml_aply(alpha,precv,x,beta,y,desc_data,trans,work,info) subroutine mlt_pre_ml_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
implicit none implicit none
! Arguments ! Arguments
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
type(mld_conelev_type), intent(in) :: precv(:) type(mld_cprec_type), intent(in) :: p
complex(psb_spk_),intent(in) :: alpha,beta complex(psb_spk_),intent(in) :: alpha,beta
complex(psb_spk_),intent(in) :: x(:) complex(psb_spk_),intent(in) :: x(:)
complex(psb_spk_),intent(inout) :: y(:) complex(psb_spk_),intent(inout) :: y(:)
@ -607,9 +607,9 @@ contains
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(precv) & ' Entry ', size(p%precv)
nlev = size(precv) nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info) allocate(mlprec_wrk(nlev),stat=info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
@ -621,7 +621,7 @@ contains
! !
! Copy the input vector X ! Copy the input vector X
! !
nc2l = psb_cd_get_local_cols(precv(1)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(1)%base_desc)
allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), & allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), &
& mlprec_wrk(1)%tx(nc2l), stat=info) & mlprec_wrk(1)%tx(nc2l), stat=info)
@ -638,8 +638,8 @@ contains
! !
! Apply the base preconditioner at the finest level ! Apply the base preconditioner at the finest level
! !
call mld_baseprec_aply(cone,precv(1)%prec,mlprec_wrk(1)%x2l,& call mld_baseprec_aply(cone,p%precv(1)%prec,mlprec_wrk(1)%x2l,&
& czero,mlprec_wrk(1)%y2l,precv(1)%base_desc,& & czero,mlprec_wrk(1)%y2l,p%precv(1)%base_desc,&
& trans,work,info) & trans,work,info)
if (info /=0) then if (info /=0) then
call psb_errpush(4010,name,a_err=' baseprec_aply') call psb_errpush(4010,name,a_err=' baseprec_aply')
@ -653,8 +653,8 @@ contains
! !
mlprec_wrk(1)%tx = mlprec_wrk(1)%x2l mlprec_wrk(1)%tx = mlprec_wrk(1)%x2l
call psb_spmm(-cone,precv(1)%base_a,mlprec_wrk(1)%y2l,& call psb_spmm(-cone,p%precv(1)%base_a,mlprec_wrk(1)%y2l,&
& cone,mlprec_wrk(1)%tx,precv(1)%base_desc,info,& & cone,mlprec_wrk(1)%tx,p%precv(1)%base_desc,info,&
& work=work,trans=trans) & work=work,trans=trans)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err=' fine level residual') call psb_errpush(4001,name,a_err=' fine level residual')
@ -668,8 +668,8 @@ contains
! !
do ilev = 2, nlev do ilev = 2, nlev
nc2l = psb_cd_get_local_cols(precv(ilev)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(ilev)%base_desc)
nr2l = psb_cd_get_local_rows(precv(ilev)%base_desc) nr2l = psb_cd_get_local_rows(p%precv(ilev)%base_desc)
allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%y2l(nc2l),& allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%y2l(nc2l),&
& mlprec_wrk(ilev)%x2l(nc2l), stat=info) & mlprec_wrk(ilev)%x2l(nc2l), stat=info)
@ -683,7 +683,7 @@ contains
! Apply prolongator transpose, i.e. restriction ! Apply prolongator transpose, i.e. restriction
call psb_map_X2Y(cone,mlprec_wrk(ilev-1)%tx,& call psb_map_X2Y(cone,mlprec_wrk(ilev-1)%tx,&
& czero,mlprec_wrk(ilev)%x2l,& & czero,mlprec_wrk(ilev)%x2l,&
& precv(ilev)%map,info,work=work) & p%precv(ilev)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction') call psb_errpush(4001,name,a_err='Error during restriction')
@ -693,17 +693,17 @@ contains
! !
! Apply the base preconditioner ! Apply the base preconditioner
! !
call mld_baseprec_aply(cone,precv(ilev)%prec,mlprec_wrk(ilev)%x2l,& call mld_baseprec_aply(cone,p%precv(ilev)%prec,mlprec_wrk(ilev)%x2l,&
& czero,mlprec_wrk(ilev)%y2l,precv(ilev)%base_desc,trans,work,info) & czero,mlprec_wrk(ilev)%y2l,p%precv(ilev)%base_desc,trans,work,info)
! !
! Compute the residual (at all levels but the coarsest one) ! Compute the residual (at all levels but the coarsest one)
! !
if (ilev < nlev) then if (ilev < nlev) then
mlprec_wrk(ilev)%tx = mlprec_wrk(ilev)%x2l mlprec_wrk(ilev)%tx = mlprec_wrk(ilev)%x2l
if (info == 0) call psb_spmm(-cone,precv(ilev)%base_a,& if (info == 0) call psb_spmm(-cone,p%precv(ilev)%base_a,&
& mlprec_wrk(ilev)%y2l,cone,mlprec_wrk(ilev)%tx,& & mlprec_wrk(ilev)%y2l,cone,mlprec_wrk(ilev)%tx,&
& precv(ilev)%base_desc,info,work=work,trans=trans) & p%precv(ilev)%base_desc,info,work=work,trans=trans)
endif endif
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error on up sweep residual') call psb_errpush(4001,name,a_err='Error on up sweep residual')
@ -722,7 +722,7 @@ contains
! !
call psb_map_Y2X(cone,mlprec_wrk(ilev+1)%y2l,& call psb_map_Y2X(cone,mlprec_wrk(ilev+1)%y2l,&
& cone,mlprec_wrk(ilev)%y2l,& & cone,mlprec_wrk(ilev)%y2l,&
& precv(ilev+1)%map,info,work=work) & p%precv(ilev+1)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during prolongation') call psb_errpush(4001,name,a_err='Error during prolongation')
@ -736,7 +736,7 @@ contains
! Compute the output vector Y ! Compute the output vector Y
! !
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,& call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,&
& precv(1)%base_desc,info) & p%precv(1)%base_desc,info)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error on final update') call psb_errpush(4001,name,a_err='Error on final update')
goto 9999 goto 9999
@ -769,7 +769,7 @@ contains
! Y = beta*Y + alpha*op(M^(-1))*X, ! Y = beta*Y + alpha*op(M^(-1))*X,
! where ! where
! - M is a hybrid multilevel domain decomposition (Schwarz) preconditioner ! - M is a hybrid multilevel domain decomposition (Schwarz) preconditioner
! associated to a certain matrix A and stored in the array precv, ! associated to a certain matrix A and stored in the array p%precv,
! - op(M^(-1)) is M^(-1) or its (conjugate) transpose, according to ! - op(M^(-1)) is M^(-1) or its (conjugate) transpose, according to
! the value of trans, ! the value of trans,
! - X and Y are vectors, ! - X and Y are vectors,
@ -784,7 +784,7 @@ contains
! !
! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners', ! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners',
! each representing the part of the preconditioner associated to a certain level. ! each representing the part of the preconditioner associated to a certain level.
! For each level ilev, the base preconditioner K(ilev) is stored in precv(ilev) ! For each level ilev, the base preconditioner K(ilev) is stored in p%precv(ilev)
! and is associated to a matrix A(ilev), obtained by 'tranferring' the original ! and is associated to a matrix A(ilev), obtained by 'tranferring' the original
! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed ! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed
! aggregation. ! aggregation.
@ -831,13 +831,13 @@ contains
! 5. Yext = beta*Yext + alpha*Y(1) ! 5. Yext = beta*Yext + alpha*Y(1)
! !
! !
subroutine mlt_post_ml_aply(alpha,precv,x,beta,y,desc_data,trans,work,info) subroutine mlt_post_ml_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
implicit none implicit none
! Arguments ! Arguments
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
type(mld_conelev_type), intent(in) :: precv(:) type(mld_cprec_type), intent(in) :: p
complex(psb_spk_),intent(in) :: alpha,beta complex(psb_spk_),intent(in) :: alpha,beta
complex(psb_spk_),intent(in) :: x(:) complex(psb_spk_),intent(in) :: x(:)
complex(psb_spk_),intent(inout) :: y(:) complex(psb_spk_),intent(inout) :: y(:)
@ -867,9 +867,9 @@ contains
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(precv) & ' Entry ', size(p%precv)
nlev = size(precv) nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info) allocate(mlprec_wrk(nlev),stat=info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
@ -885,7 +885,7 @@ contains
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' desc_data status',allocated(desc_data%matrix_data) & ' desc_data status',allocated(desc_data%matrix_data)
nc2l = psb_cd_get_local_cols(precv(1)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(1)%base_desc)
allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), & allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), &
& mlprec_wrk(1)%tx(nc2l), stat=info) & mlprec_wrk(1)%tx(nc2l), stat=info)
@ -897,9 +897,9 @@ contains
end if end if
call psb_geaxpby(cone,x,czero,mlprec_wrk(1)%tx,& call psb_geaxpby(cone,x,czero,mlprec_wrk(1)%tx,&
& precv(1)%base_desc,info) & p%precv(1)%base_desc,info)
call psb_geaxpby(cone,x,czero,mlprec_wrk(1)%x2l,& call psb_geaxpby(cone,x,czero,mlprec_wrk(1)%x2l,&
& precv(1)%base_desc,info) & p%precv(1)%base_desc,info)
! !
! STEP 2 ! STEP 2
@ -908,13 +908,13 @@ contains
! !
do ilev=2, nlev do ilev=2, nlev
nc2l = psb_cd_get_local_cols(precv(ilev)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(ilev)%base_desc)
nr2l = psb_cd_get_local_rows(precv(ilev)%base_desc) nr2l = psb_cd_get_local_rows(p%precv(ilev)%base_desc)
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name), & & write(debug_unit,*) me,' ',trim(name), &
& ' starting up sweep ',& & ' starting up sweep ',&
& ilev,allocated(precv(ilev)%iprcparm),nc2l, nr2l & ilev,allocated(p%precv(ilev)%iprcparm),nc2l, nr2l
allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%y2l(nc2l),& allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%y2l(nc2l),&
& mlprec_wrk(ilev)%x2l(nc2l), stat=info) & mlprec_wrk(ilev)%x2l(nc2l), stat=info)
@ -929,7 +929,7 @@ contains
! Apply prolongator transpose, i.e. restriction ! Apply prolongator transpose, i.e. restriction
call psb_map_X2Y(cone,mlprec_wrk(ilev-1)%x2l,& call psb_map_X2Y(cone,mlprec_wrk(ilev-1)%x2l,&
& czero,mlprec_wrk(ilev)%x2l,& & czero,mlprec_wrk(ilev)%x2l,&
& precv(ilev)%map,info,work=work) & p%precv(ilev)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction') call psb_errpush(4001,name,a_err='Error during restriction')
@ -940,7 +940,7 @@ contains
! update x2l ! update x2l
! !
call psb_geaxpby(cone,mlprec_wrk(ilev)%x2l,czero,mlprec_wrk(ilev)%tx,& call psb_geaxpby(cone,mlprec_wrk(ilev)%x2l,czero,mlprec_wrk(ilev)%tx,&
& precv(ilev)%base_desc,info) & p%precv(ilev)%base_desc,info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4001,name,a_err='Error in update') call psb_errpush(4001,name,a_err='Error in update')
goto 9999 goto 9999
@ -957,8 +957,8 @@ contains
! !
! Apply the base preconditioner at the coarsest level ! Apply the base preconditioner at the coarsest level
! !
call mld_baseprec_aply(cone,precv(nlev)%prec,mlprec_wrk(nlev)%x2l, & call mld_baseprec_aply(cone,p%precv(nlev)%prec,mlprec_wrk(nlev)%x2l, &
& czero, mlprec_wrk(nlev)%y2l,precv(nlev)%base_desc,trans,work,info) & czero, mlprec_wrk(nlev)%y2l,p%precv(nlev)%base_desc,trans,work,info)
if (info /=0) then if (info /=0) then
call psb_errpush(4010,name,a_err='baseprec_aply') call psb_errpush(4010,name,a_err='baseprec_aply')
goto 9999 goto 9999
@ -983,7 +983,7 @@ contains
! !
call psb_map_Y2X(cone,mlprec_wrk(ilev+1)%y2l,& call psb_map_Y2X(cone,mlprec_wrk(ilev+1)%y2l,&
& czero,mlprec_wrk(ilev)%y2l,& & czero,mlprec_wrk(ilev)%y2l,&
& precv(ilev+1)%map,info,work=work) & p%precv(ilev+1)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during prolongation') call psb_errpush(4001,name,a_err='Error during prolongation')
@ -993,15 +993,15 @@ contains
! !
! Compute the residual ! Compute the residual
! !
call psb_spmm(-cone,precv(ilev)%base_a,mlprec_wrk(ilev)%y2l,& call psb_spmm(-cone,p%precv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
& cone,mlprec_wrk(ilev)%tx,precv(ilev)%base_desc,info,& & cone,mlprec_wrk(ilev)%tx,p%precv(ilev)%base_desc,info,&
& work=work,trans=trans) & work=work,trans=trans)
! !
! Apply the base preconditioner ! Apply the base preconditioner
! !
if (info == 0) call mld_baseprec_aply(cone,precv(ilev)%prec,& if (info == 0) call mld_baseprec_aply(cone,p%precv(ilev)%prec,&
& mlprec_wrk(ilev)%tx,cone,mlprec_wrk(ilev)%y2l,precv(ilev)%base_desc,& & mlprec_wrk(ilev)%tx,cone,mlprec_wrk(ilev)%y2l,p%precv(ilev)%base_desc,&
& trans,work,info) & trans,work,info)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err=' spmm/baseprec_aply') call psb_errpush(4001,name,a_err=' spmm/baseprec_aply')
@ -1018,7 +1018,7 @@ contains
! !
! Compute the output vector Y ! Compute the output vector Y
! !
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,precv(1)%base_desc,info) call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,p%precv(1)%base_desc,info)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err=' Final update') call psb_errpush(4001,name,a_err=' Final update')
@ -1055,7 +1055,7 @@ contains
! where ! where
! - M is a symmetrized hybrid multilevel domain decomposition (Schwarz) ! - M is a symmetrized hybrid multilevel domain decomposition (Schwarz)
! preconditioner associated to a certain matrix A and stored in the array ! preconditioner associated to a certain matrix A and stored in the array
! precv, ! p%precv,
! - op(M^(-1)) is M^(-1) or its (conjugate) transpose, according to ! - op(M^(-1)) is M^(-1) or its (conjugate) transpose, according to
! the value of trans, ! the value of trans,
! - X and Y are vectors, ! - X and Y are vectors,
@ -1071,7 +1071,7 @@ contains
! !
! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners', ! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners',
! each representing the part of the preconditioner associated to a certain level. ! each representing the part of the preconditioner associated to a certain level.
! For each level ilev, the base preconditioner K(ilev) is stored in precv(ilev) ! For each level ilev, the base preconditioner K(ilev) is stored in p%precv(ilev)
! and is associated to a matrix A(ilev), obtained by 'tranferring' the original ! and is associated to a matrix A(ilev), obtained by 'tranferring' the original
! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed ! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed
! aggregation. ! aggregation.
@ -1129,13 +1129,13 @@ contains
! !
! 6. Yext = beta*Yext + alpha*Y(1) ! 6. Yext = beta*Yext + alpha*Y(1)
! !
subroutine mlt_twoside_ml_aply(alpha,precv,x,beta,y,desc_data,trans,work,info) subroutine mlt_twoside_ml_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
implicit none implicit none
! Arguments ! Arguments
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
type(mld_conelev_type), intent(in) :: precv(:) type(mld_cprec_type), intent(in) :: p
complex(psb_spk_),intent(in) :: alpha,beta complex(psb_spk_),intent(in) :: alpha,beta
complex(psb_spk_),intent(in) :: x(:) complex(psb_spk_),intent(in) :: x(:)
complex(psb_spk_),intent(inout) :: y(:) complex(psb_spk_),intent(inout) :: y(:)
@ -1165,9 +1165,9 @@ contains
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(precv) & ' Entry ', size(p%precv)
nlev = size(precv) nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info) allocate(mlprec_wrk(nlev),stat=info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
@ -1178,7 +1178,7 @@ contains
! !
! Copy the input vector X ! Copy the input vector X
! !
nc2l = psb_cd_get_local_cols(precv(1)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(1)%base_desc)
allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), & allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), &
& mlprec_wrk(1)%ty(nc2l), mlprec_wrk(1)%tx(nc2l), stat=info) & mlprec_wrk(1)%ty(nc2l), mlprec_wrk(1)%tx(nc2l), stat=info)
@ -1191,17 +1191,17 @@ contains
end if end if
call psb_geaxpby(cone,x,czero,mlprec_wrk(1)%x2l,& call psb_geaxpby(cone,x,czero,mlprec_wrk(1)%x2l,&
& precv(1)%base_desc,info) & p%precv(1)%base_desc,info)
call psb_geaxpby(cone,x,czero,mlprec_wrk(1)%tx,& call psb_geaxpby(cone,x,czero,mlprec_wrk(1)%tx,&
& precv(1)%base_desc,info) & p%precv(1)%base_desc,info)
! !
! STEP 2 ! STEP 2
! !
! Apply the base preconditioner at the finest level ! Apply the base preconditioner at the finest level
! !
call mld_baseprec_aply(cone,precv(1)%prec,mlprec_wrk(1)%x2l,& call mld_baseprec_aply(cone,p%precv(1)%prec,mlprec_wrk(1)%x2l,&
& czero,mlprec_wrk(1)%y2l,precv(1)%base_desc,& & czero,mlprec_wrk(1)%y2l,p%precv(1)%base_desc,&
& trans,work,info) & trans,work,info)
! !
! STEP 3 ! STEP 3
@ -1209,8 +1209,8 @@ contains
! Compute the residual at the finest level ! Compute the residual at the finest level
! !
mlprec_wrk(1)%ty = mlprec_wrk(1)%x2l mlprec_wrk(1)%ty = mlprec_wrk(1)%x2l
if (info == 0) call psb_spmm(-cone,precv(1)%base_a,mlprec_wrk(1)%y2l,& if (info == 0) call psb_spmm(-cone,p%precv(1)%base_a,mlprec_wrk(1)%y2l,&
& cone,mlprec_wrk(1)%ty,precv(1)%base_desc,info,& & cone,mlprec_wrk(1)%ty,p%precv(1)%base_desc,info,&
& work=work,trans=trans) & work=work,trans=trans)
if (info /=0) then if (info /=0) then
call psb_errpush(4010,name,a_err='Fine level baseprec/residual') call psb_errpush(4010,name,a_err='Fine level baseprec/residual')
@ -1224,8 +1224,8 @@ contains
! !
do ilev = 2, nlev do ilev = 2, nlev
nc2l = psb_cd_get_local_cols(precv(ilev)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(ilev)%base_desc)
nr2l = psb_cd_get_local_rows(precv(ilev)%base_desc) nr2l = psb_cd_get_local_rows(p%precv(ilev)%base_desc)
allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%ty(nc2l),& allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%ty(nc2l),&
& mlprec_wrk(ilev)%y2l(nc2l),mlprec_wrk(ilev)%x2l(nc2l), stat=info) & mlprec_wrk(ilev)%y2l(nc2l),mlprec_wrk(ilev)%x2l(nc2l), stat=info)
@ -1240,7 +1240,7 @@ contains
! Apply prolongator transpose, i.e. restriction ! Apply prolongator transpose, i.e. restriction
call psb_map_X2Y(cone,mlprec_wrk(ilev-1)%ty,& call psb_map_X2Y(cone,mlprec_wrk(ilev-1)%ty,&
& czero,mlprec_wrk(ilev)%x2l,& & czero,mlprec_wrk(ilev)%x2l,&
& precv(ilev)%map,info,work=work) & p%precv(ilev)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction') call psb_errpush(4001,name,a_err='Error during restriction')
@ -1248,21 +1248,21 @@ contains
end if end if
call psb_geaxpby(cone,mlprec_wrk(ilev)%x2l,czero,mlprec_wrk(ilev)%tx,& call psb_geaxpby(cone,mlprec_wrk(ilev)%x2l,czero,mlprec_wrk(ilev)%tx,&
& precv(ilev)%base_desc,info) & p%precv(ilev)%base_desc,info)
! !
! Apply the base preconditioner ! Apply the base preconditioner
! !
if (info == 0) call mld_baseprec_aply(cone,precv(ilev)%prec,& if (info == 0) call mld_baseprec_aply(cone,p%precv(ilev)%prec,&
& mlprec_wrk(ilev)%x2l,czero,mlprec_wrk(ilev)%y2l,& & mlprec_wrk(ilev)%x2l,czero,mlprec_wrk(ilev)%y2l,&
&precv(ilev)%base_desc,trans,work,info) &p%precv(ilev)%base_desc,trans,work,info)
! !
! Compute the residual (at all levels but the coarsest one) ! Compute the residual (at all levels but the coarsest one)
! !
if(ilev < nlev) then if(ilev < nlev) then
mlprec_wrk(ilev)%ty = mlprec_wrk(ilev)%x2l mlprec_wrk(ilev)%ty = mlprec_wrk(ilev)%x2l
if (info == 0) call psb_spmm(-cone,precv(ilev)%base_a,& if (info == 0) call psb_spmm(-cone,p%precv(ilev)%base_a,&
& mlprec_wrk(ilev)%y2l,cone,mlprec_wrk(ilev)%ty,& & mlprec_wrk(ilev)%y2l,cone,mlprec_wrk(ilev)%ty,&
& precv(ilev)%base_desc,info,work=work,trans=trans) & p%precv(ilev)%base_desc,info,work=work,trans=trans)
endif endif
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='baseprec_aply/residual') call psb_errpush(4001,name,a_err='baseprec_aply/residual')
@ -1283,7 +1283,7 @@ contains
! !
call psb_map_Y2X(cone,mlprec_wrk(ilev+1)%y2l,& call psb_map_Y2X(cone,mlprec_wrk(ilev+1)%y2l,&
& cone,mlprec_wrk(ilev)%y2l,& & cone,mlprec_wrk(ilev)%y2l,&
& precv(ilev+1)%map,info,work=work) & p%precv(ilev+1)%map,info,work=work)
if (info /=0 ) then if (info /=0 ) then
call psb_errpush(4001,name,a_err='Error during restriction') call psb_errpush(4001,name,a_err='Error during restriction')
@ -1293,14 +1293,14 @@ contains
! !
! Compute the residual ! Compute the residual
! !
call psb_spmm(-cone,precv(ilev)%base_a,mlprec_wrk(ilev)%y2l,& call psb_spmm(-cone,p%precv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
& cone,mlprec_wrk(ilev)%tx,precv(ilev)%base_desc,info,& & cone,mlprec_wrk(ilev)%tx,p%precv(ilev)%base_desc,info,&
& work=work,trans=trans) & work=work,trans=trans)
! !
! Apply the base preconditioner ! Apply the base preconditioner
! !
if (info == 0) call mld_baseprec_aply(cone,precv(ilev)%prec,mlprec_wrk(ilev)%tx,& if (info == 0) call mld_baseprec_aply(cone,p%precv(ilev)%prec,mlprec_wrk(ilev)%tx,&
& cone,mlprec_wrk(ilev)%y2l,precv(ilev)%base_desc, trans, work,info) & cone,mlprec_wrk(ilev)%y2l,p%precv(ilev)%base_desc, trans, work,info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4001,name,a_err='Error: residual/baseprec_aply') call psb_errpush(4001,name,a_err='Error: residual/baseprec_aply')
goto 9999 goto 9999
@ -1313,7 +1313,7 @@ contains
! Compute the output vector Y ! Compute the output vector Y
! !
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,& call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,&
& precv(1)%base_desc,info) & p%precv(1)%base_desc,info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4001,name,a_err='Error final update') call psb_errpush(4001,name,a_err='Error final update')

2
mlprec/mld_cprecaply.f90
Unescape Escape View File

@ -127,7 +127,7 @@ subroutine mld_cprecaply(prec,x,y,desc_data,info,trans,work)
goto 9999 goto 9999
end if end if
if (size(prec%precv) >1) then if (size(prec%precv) >1) then
call mld_mlprec_aply(cone,prec%precv,x,czero,y,desc_data,trans_,work_,info) call mld_mlprec_aply(cone,prec,x,czero,y,desc_data,trans_,work_,info)
if(info /= 0) then if(info /= 0) then
call psb_errpush(4010,name,a_err='mld_cmlprec_aply') call psb_errpush(4010,name,a_err='mld_cmlprec_aply')
goto 9999 goto 9999

229
mlprec/mld_dmlprec_aply.f90
Unescape Escape View File

@ -46,7 +46,7 @@
! Y = beta*Y + alpha*op(M^(-1))*X, ! Y = beta*Y + alpha*op(M^(-1))*X,
! where ! where
! - M is a multilevel domain decomposition (Schwarz) preconditioner associated ! - M is a multilevel domain decomposition (Schwarz) preconditioner associated
! to a certain matrix A and stored in the array precv, ! to a certain matrix A and stored in the array p%precv,
! - op(M^(-1)) is M^(-1) or its transpose, according to the value of trans, ! - op(M^(-1)) is M^(-1) or its transpose, according to the value of trans,
! - X and Y are vectors, ! - X and Y are vectors,
! - alpha and beta are scalars. ! - alpha and beta are scalars.
@ -57,7 +57,7 @@
! !
! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners', ! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners',
! each representing the part of the preconditioner associated to a certain level. ! each representing the part of the preconditioner associated to a certain level.
! For each level ilev, the preconditioner K(ilev) is stored in precv(ilev) ! For each level ilev, the preconditioner K(ilev) is stored in p%precv(ilev)
! and is associated to a matrix A(ilev), obtained by 'tranferring' the original ! and is associated to a matrix A(ilev), obtained by 'tranferring' the original
! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed ! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed
! aggregation. ! aggregation.
@ -78,41 +78,41 @@
! Arguments: ! Arguments:
! alpha - real(psb_dpk_), input. ! alpha - real(psb_dpk_), input.
! The scalar alpha. ! The scalar alpha.
! precv - type(mld_donelev_type), dimension(:), input. ! p - type(mld_dprec_type), input.
! The array of one-level preconditioner data structures containing the ! The array of one-level preconditioner data structures containing the
! local parts of the preconditioners to be applied at each level. ! local parts of the preconditioners to be applied at each level.
! Note that nlev = size(precv) = number of levels. ! Note that nlev = size(p%precv) = number of levels.
! precv(ilev)%prec - type(psb_dbaseprec_type) ! p%precv(ilev)%prec - type(psb_dbaseprec_type)
! The "base" preconditioner for the current level ! The "base" preconditioner for the current level
! precv(ilev)%ac - type(psb_dspmat_type) ! p%precv(ilev)%ac - type(psb_dspmat_type)
! The local part of the matrix A(ilev). ! The local part of the matrix A(ilev).
! precv(ilev)%desc_ac - type(psb_desc_type). ! p%precv(ilev)%desc_ac - type(psb_desc_type).
! The communication descriptor associated to the sparse ! The communication descriptor associated to the sparse
! matrix A(ilev) ! matrix A(ilev)
! precv(ilev)%map - type(psb_inter_desc_type) ! p%precv(ilev)%map - type(psb_inter_desc_type)
! Stores the linear operators mapping between levels ! Stores the linear operators mapping between levels
! (ilev-1) and (ilev). These are the restriction and ! (ilev-1) and (ilev). These are the restriction and
! prolongation operators described in the sequel. ! prolongation operators described in the sequel.
! precv(ilev)%iprcparm - integer, dimension(:), allocatable. ! p%precv(ilev)%iprcparm - integer, dimension(:), allocatable.
! The integer parameters defining the multilevel ! The integer parameters defining the multilevel
! strategy ! strategy
! precv(ilev)%rprcparm - real(psb_dpk_), dimension(:), allocatable. ! p%precv(ilev)%rprcparm - real(psb_dpk_), dimension(:), allocatable.
! The real parameters defining the multilevel strategy ! The real parameters defining the multilevel strategy
! precv(ilev)%mlia - integer, dimension(:), allocatable. ! p%precv(ilev)%mlia - integer, dimension(:), allocatable.
! The aggregation map (ilev-1) --> (ilev). ! The aggregation map (ilev-1) --> (ilev).
! In case of non-smoothed aggregation, it is used ! In case of non-smoothed aggregation, it is used
! instead of mld_sm_pr_. ! instead of mld_sm_pr_.
! precv(ilev)%nlaggr - integer, dimension(:), allocatable. ! p%precv(ilev)%nlaggr - integer, dimension(:), allocatable.
! The number of aggregates (rows of A(ilev)) on the ! The number of aggregates (rows of A(ilev)) on the
! various processes. ! various processes.
! precv(ilev)%base_a - type(psb_dspmat_type), pointer. ! p%precv(ilev)%base_a - type(psb_dspmat_type), pointer.
! Pointer (really a pointer!) to the base matrix of ! Pointer (really a pointer!) to the base matrix of
! the current level, i.e. the local part of A(ilev); ! the current level, i.e. the local part of A(ilev);
! so we have a unified treatment of residuals. We ! so we have a unified treatment of residuals. We
! need this to avoid passing explicitly the matrix ! need this to avoid passing explicitly the matrix
! A(ilev) to the routine which applies the ! A(ilev) to the routine which applies the
! preconditioner. ! preconditioner.
! precv(ilev)%base_desc - type(psb_desc_type), pointer. ! p%precv(ilev)%base_desc - type(psb_desc_type), pointer.
! Pointer to the communication descriptor associated ! Pointer to the communication descriptor associated
! to the sparse matrix pointed by base_a. ! to the sparse matrix pointed by base_a.
! !
@ -136,10 +136,10 @@
! Note that when the LU factorization of the matrix A(ilev) is computed instead of ! Note that when the LU factorization of the matrix A(ilev) is computed instead of
! the ILU one, by using UMFPACK or SuperLU, the corresponding L and U factors ! the ILU one, by using UMFPACK or SuperLU, the corresponding L and U factors
! are stored in data structures provided by UMFPACK or SuperLU and pointed by ! are stored in data structures provided by UMFPACK or SuperLU and pointed by
! precv(ilev)%prec%iprcparm(mld_umf_ptr) or precv(ilev)%prec%iprcparm(mld_slu_ptr), ! p%precv(ilev)%prec%iprcparm(mld_umf_ptr) or p%precv(ilev)%prec%iprcparm(mld_slu_ptr),
! respectively. ! respectively.
! !
subroutine mld_dmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info) subroutine mld_dmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
use psb_base_mod use psb_base_mod
use mld_inner_mod, mld_protect_name => mld_dmlprec_aply use mld_inner_mod, mld_protect_name => mld_dmlprec_aply
@ -148,7 +148,8 @@ subroutine mld_dmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
! Arguments ! Arguments
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
type(mld_donelev_type), intent(in) :: precv(:) type(mld_dprec_type), intent(in) :: p
!!$ type(mld_donelev_type), intent(in) :: precv(:)
real(psb_dpk_),intent(in) :: alpha,beta real(psb_dpk_),intent(in) :: alpha,beta
real(psb_dpk_),intent(in) :: x(:) real(psb_dpk_),intent(in) :: x(:)
real(psb_dpk_),intent(inout) :: y(:) real(psb_dpk_),intent(inout) :: y(:)
@ -173,11 +174,11 @@ subroutine mld_dmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(precv) & ' Entry ', size(p%precv)
trans_ = psb_toupper(trans) trans_ = psb_toupper(trans)
select case(precv(2)%iprcparm(mld_ml_type_)) select case(p%precv(2)%iprcparm(mld_ml_type_))
case(mld_no_ml_) case(mld_no_ml_)
! !
@ -191,7 +192,7 @@ subroutine mld_dmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
! Additive multilevel ! Additive multilevel
! !
call add_ml_aply(alpha,precv,x,beta,y,desc_data,trans_,work,info) call add_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info)
case(mld_mult_ml_) case(mld_mult_ml_)
! !
@ -202,15 +203,15 @@ subroutine mld_dmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
! Note that the transpose switches pre <-> post. ! Note that the transpose switches pre <-> post.
! !
select case(precv(2)%iprcparm(mld_smoother_pos_)) select case(p%precv(2)%iprcparm(mld_smoother_pos_))
case(mld_post_smooth_) case(mld_post_smooth_)
select case (trans_) select case (trans_)
case('N') case('N')
call mlt_post_ml_aply(alpha,precv,x,beta,y,desc_data,trans_,work,info) call mlt_post_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info)
case('T','C') case('T','C')
call mlt_pre_ml_aply(alpha,precv,x,beta,y,desc_data,trans_,work,info) call mlt_pre_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info)
case default case default
info = 4001 info = 4001
call psb_errpush(info,name,a_err='invalid trans') call psb_errpush(info,name,a_err='invalid trans')
@ -221,9 +222,9 @@ subroutine mld_dmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
select case (trans_) select case (trans_)
case('N') case('N')
call mlt_pre_ml_aply(alpha,precv,x,beta,y,desc_data,trans_,work,info) call mlt_pre_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info)
case('T','C') case('T','C')
call mlt_post_ml_aply(alpha,precv,x,beta,y,desc_data,trans_,work,info) call mlt_post_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info)
case default case default
info = 4001 info = 4001
call psb_errpush(info,name,a_err='invalid trans') call psb_errpush(info,name,a_err='invalid trans')
@ -232,12 +233,12 @@ subroutine mld_dmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
case(mld_twoside_smooth_) case(mld_twoside_smooth_)
call mlt_twoside_ml_aply(alpha,precv,x,beta,y,desc_data,trans_,work,info) call mlt_twoside_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info)
case default case default
info = 4013 info = 4013
call psb_errpush(info,name,a_err='invalid smooth_pos',& call psb_errpush(info,name,a_err='invalid smooth_pos',&
& i_Err=(/precv(2)%iprcparm(mld_smoother_pos_),0,0,0,0/)) & i_Err=(/p%precv(2)%iprcparm(mld_smoother_pos_),0,0,0,0/))
goto 9999 goto 9999
end select end select
@ -245,7 +246,7 @@ subroutine mld_dmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
case default case default
info = 4013 info = 4013
call psb_errpush(info,name,a_err='invalid mltype',& call psb_errpush(info,name,a_err='invalid mltype',&
& i_Err=(/precv(2)%iprcparm(mld_ml_type_),0,0,0,0/)) & i_Err=(/p%precv(2)%iprcparm(mld_ml_type_),0,0,0,0/))
goto 9999 goto 9999
end select end select
@ -334,13 +335,14 @@ contains
! !
! 4. Yext = beta*Yext + alpha*Y(1) ! 4. Yext = beta*Yext + alpha*Y(1)
! !
subroutine add_ml_aply(alpha,precv,x,beta,y,desc_data,trans,work,info) subroutine add_ml_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
implicit none implicit none
! Arguments ! Arguments
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
type(mld_donelev_type), intent(in) :: precv(:) type(mld_dprec_type), intent(in) :: p
!!$ type(mld_donelev_type), intent(in) :: precv(:)
real(psb_dpk_),intent(in) :: alpha,beta real(psb_dpk_),intent(in) :: alpha,beta
real(psb_dpk_),intent(in) :: x(:) real(psb_dpk_),intent(in) :: x(:)
real(psb_dpk_),intent(inout) :: y(:) real(psb_dpk_),intent(inout) :: y(:)
@ -370,9 +372,9 @@ contains
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(precv) & ' Entry ', size(p%precv)
nlev = size(precv) nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info) allocate(mlprec_wrk(nlev),stat=info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
@ -395,8 +397,8 @@ contains
mlprec_wrk(1)%x2l(:) = x(:) mlprec_wrk(1)%x2l(:) = x(:)
mlprec_wrk(1)%y2l(:) = dzero mlprec_wrk(1)%y2l(:) = dzero
call mld_baseprec_aply(alpha,precv(1)%prec,x,beta,y,& call mld_baseprec_aply(alpha,p%precv(1)%prec,x,beta,y,&
& precv(1)%base_desc,trans,work,info) & p%precv(1)%base_desc,trans,work,info)
if (info /=0) then if (info /=0) then
call psb_errpush(4010,name,a_err='baseprec_aply') call psb_errpush(4010,name,a_err='baseprec_aply')
goto 9999 goto 9999
@ -407,8 +409,8 @@ contains
! For each level except the finest one ... ! For each level except the finest one ...
! !
do ilev = 2, nlev do ilev = 2, nlev
nc2l = psb_cd_get_local_cols(precv(ilev)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(ilev)%base_desc)
nr2l = psb_cd_get_local_rows(precv(ilev)%base_desc) nr2l = psb_cd_get_local_rows(p%precv(ilev)%base_desc)
allocate(mlprec_wrk(ilev)%x2l(nc2l),mlprec_wrk(ilev)%y2l(nc2l),& allocate(mlprec_wrk(ilev)%x2l(nc2l),mlprec_wrk(ilev)%y2l(nc2l),&
& stat=info) & stat=info)
if (info /= 0) then if (info /= 0) then
@ -421,7 +423,7 @@ contains
! Apply prolongator transpose, i.e. restriction ! Apply prolongator transpose, i.e. restriction
call psb_map_X2Y(done,mlprec_wrk(ilev-1)%x2l,& call psb_map_X2Y(done,mlprec_wrk(ilev-1)%x2l,&
& dzero,mlprec_wrk(ilev)%x2l,& & dzero,mlprec_wrk(ilev)%x2l,&
& precv(ilev)%map,info,work=work) & p%precv(ilev)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction') call psb_errpush(4001,name,a_err='Error during restriction')
@ -431,9 +433,9 @@ contains
! !
! Apply the base preconditioner ! Apply the base preconditioner
! !
call mld_baseprec_aply(done,precv(ilev)%prec,& call mld_baseprec_aply(done,p%precv(ilev)%prec,&
& mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%y2l,& & mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%y2l,&
& precv(ilev)%base_desc, trans,work,info) & p%precv(ilev)%base_desc, trans,work,info)
enddo enddo
@ -444,15 +446,15 @@ contains
! !
do ilev =nlev,2,-1 do ilev =nlev,2,-1
nc2l = psb_cd_get_local_cols(precv(ilev)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(ilev)%base_desc)
nr2l = psb_cd_get_local_rows(precv(ilev)%base_desc) nr2l = psb_cd_get_local_rows(p%precv(ilev)%base_desc)
! !
! Apply prolongator ! Apply prolongator
! !
call psb_map_Y2X(done,mlprec_wrk(ilev)%y2l,& call psb_map_Y2X(done,mlprec_wrk(ilev)%y2l,&
& done,mlprec_wrk(ilev-1)%y2l,& & done,mlprec_wrk(ilev-1)%y2l,&
& precv(ilev)%map,info,work=work) & p%precv(ilev)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during prolongation') call psb_errpush(4001,name,a_err='Error during prolongation')
@ -465,7 +467,7 @@ contains
! !
! Compute the output vector Y ! Compute the output vector Y
! !
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,done,y,precv(1)%base_desc,info) call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,done,y,p%precv(1)%base_desc,info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4001,name,a_err='Error on final update') call psb_errpush(4001,name,a_err='Error on final update')
goto 9999 goto 9999
@ -499,7 +501,7 @@ contains
! Y = beta*Y + alpha*op(M^(-1))*X, ! Y = beta*Y + alpha*op(M^(-1))*X,
! where ! where
! - M is a hybrid multilevel domain decomposition (Schwarz) preconditioner ! - M is a hybrid multilevel domain decomposition (Schwarz) preconditioner
! associated to a certain matrix A and stored in the array precv, ! associated to a certain matrix A and stored in the array p%precv,
! - op(M^(-1)) is M^(-1) or its transpose, according to the value of trans, ! - op(M^(-1)) is M^(-1) or its transpose, according to the value of trans,
! - X and Y are vectors, ! - X and Y are vectors,
! - alpha and beta are scalars. ! - alpha and beta are scalars.
@ -513,7 +515,7 @@ contains
! !
! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners', ! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners',
! each representing the part of the preconditioner associated to a certain level. ! each representing the part of the preconditioner associated to a certain level.
! For each level ilev, the base preconditioner K(ilev) is stored in precv(ilev) ! For each level ilev, the base preconditioner K(ilev) is stored in p%precv(ilev)
! and is associated to a matrix A(ilev), obtained by 'tranferring' the original ! and is associated to a matrix A(ilev), obtained by 'tranferring' the original
! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed ! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed
! aggregation. ! aggregation.
@ -569,13 +571,14 @@ contains
! 6. Yext = beta*Yext + alpha*Y(1) ! 6. Yext = beta*Yext + alpha*Y(1)
! !
! !
subroutine mlt_pre_ml_aply(alpha,precv,x,beta,y,desc_data,trans,work,info) subroutine mlt_pre_ml_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
implicit none implicit none
! Arguments ! Arguments
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
type(mld_donelev_type), intent(in) :: precv(:) type(mld_dprec_type), intent(in) :: p
!!$ type(mld_donelev_type), intent(in) :: precv(:)
real(psb_dpk_),intent(in) :: alpha,beta real(psb_dpk_),intent(in) :: alpha,beta
real(psb_dpk_),intent(in) :: x(:) real(psb_dpk_),intent(in) :: x(:)
real(psb_dpk_),intent(inout) :: y(:) real(psb_dpk_),intent(inout) :: y(:)
@ -605,9 +608,9 @@ contains
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(precv) & ' Entry ', size(p%precv)
nlev = size(precv) nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info) allocate(mlprec_wrk(nlev),stat=info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
@ -619,7 +622,7 @@ contains
! !
! Copy the input vector X ! Copy the input vector X
! !
nc2l = psb_cd_get_local_cols(precv(1)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(1)%base_desc)
allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), & allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), &
& mlprec_wrk(1)%tx(nc2l), stat=info) & mlprec_wrk(1)%tx(nc2l), stat=info)
@ -636,8 +639,8 @@ contains
! !
! Apply the base preconditioner at the finest level ! Apply the base preconditioner at the finest level
! !
call mld_baseprec_aply(done,precv(1)%prec,mlprec_wrk(1)%x2l,& call mld_baseprec_aply(done,p%precv(1)%prec,mlprec_wrk(1)%x2l,&
& dzero,mlprec_wrk(1)%y2l,precv(1)%base_desc,& & dzero,mlprec_wrk(1)%y2l,p%precv(1)%base_desc,&
& trans,work,info) & trans,work,info)
if (info /=0) then if (info /=0) then
call psb_errpush(4010,name,a_err=' baseprec_aply') call psb_errpush(4010,name,a_err=' baseprec_aply')
@ -651,8 +654,8 @@ contains
! !
mlprec_wrk(1)%tx = mlprec_wrk(1)%x2l mlprec_wrk(1)%tx = mlprec_wrk(1)%x2l
call psb_spmm(-done,precv(1)%base_a,mlprec_wrk(1)%y2l,& call psb_spmm(-done,p%precv(1)%base_a,mlprec_wrk(1)%y2l,&
& done,mlprec_wrk(1)%tx,precv(1)%base_desc,info,& & done,mlprec_wrk(1)%tx,p%precv(1)%base_desc,info,&
& work=work,trans=trans) & work=work,trans=trans)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err=' fine level residual') call psb_errpush(4001,name,a_err=' fine level residual')
@ -666,8 +669,8 @@ contains
! !
do ilev = 2, nlev do ilev = 2, nlev
nc2l = psb_cd_get_local_cols(precv(ilev)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(ilev)%base_desc)
nr2l = psb_cd_get_local_rows(precv(ilev)%base_desc) nr2l = psb_cd_get_local_rows(p%precv(ilev)%base_desc)
allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%y2l(nc2l),& allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%y2l(nc2l),&
& mlprec_wrk(ilev)%x2l(nc2l), stat=info) & mlprec_wrk(ilev)%x2l(nc2l), stat=info)
@ -681,7 +684,7 @@ contains
! Apply prolongator transpose, i.e. restriction ! Apply prolongator transpose, i.e. restriction
call psb_map_X2Y(done,mlprec_wrk(ilev-1)%tx,& call psb_map_X2Y(done,mlprec_wrk(ilev-1)%tx,&
& dzero,mlprec_wrk(ilev)%x2l,& & dzero,mlprec_wrk(ilev)%x2l,&
& precv(ilev)%map,info,work=work) & p%precv(ilev)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction') call psb_errpush(4001,name,a_err='Error during restriction')
@ -691,17 +694,17 @@ contains
! !
! Apply the base preconditioner ! Apply the base preconditioner
! !
call mld_baseprec_aply(done,precv(ilev)%prec,mlprec_wrk(ilev)%x2l,& call mld_baseprec_aply(done,p%precv(ilev)%prec,mlprec_wrk(ilev)%x2l,&
& dzero,mlprec_wrk(ilev)%y2l,precv(ilev)%base_desc,trans,work,info) & dzero,mlprec_wrk(ilev)%y2l,p%precv(ilev)%base_desc,trans,work,info)
! !
! Compute the residual (at all levels but the coarsest one) ! Compute the residual (at all levels but the coarsest one)
! !
if (ilev < nlev) then if (ilev < nlev) then
mlprec_wrk(ilev)%tx = mlprec_wrk(ilev)%x2l mlprec_wrk(ilev)%tx = mlprec_wrk(ilev)%x2l
if (info == 0) call psb_spmm(-done,precv(ilev)%base_a,& if (info == 0) call psb_spmm(-done,p%precv(ilev)%base_a,&
& mlprec_wrk(ilev)%y2l,done,mlprec_wrk(ilev)%tx,& & mlprec_wrk(ilev)%y2l,done,mlprec_wrk(ilev)%tx,&
& precv(ilev)%base_desc,info,work=work,trans=trans) & p%precv(ilev)%base_desc,info,work=work,trans=trans)
endif endif
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error on up sweep residual') call psb_errpush(4001,name,a_err='Error on up sweep residual')
@ -720,7 +723,7 @@ contains
! !
call psb_map_Y2X(done,mlprec_wrk(ilev+1)%y2l,& call psb_map_Y2X(done,mlprec_wrk(ilev+1)%y2l,&
& done,mlprec_wrk(ilev)%y2l,& & done,mlprec_wrk(ilev)%y2l,&
& precv(ilev+1)%map,info,work=work) & p%precv(ilev+1)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during prolongation') call psb_errpush(4001,name,a_err='Error during prolongation')
@ -734,7 +737,7 @@ contains
! Compute the output vector Y ! Compute the output vector Y
! !
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,& call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,&
& precv(1)%base_desc,info) & p%precv(1)%base_desc,info)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error on final update') call psb_errpush(4001,name,a_err='Error on final update')
goto 9999 goto 9999
@ -767,7 +770,7 @@ contains
! Y = beta*Y + alpha*op(M^(-1))*X, ! Y = beta*Y + alpha*op(M^(-1))*X,
! where ! where
! - M is a hybrid multilevel domain decomposition (Schwarz) preconditioner ! - M is a hybrid multilevel domain decomposition (Schwarz) preconditioner
! associated to a certain matrix A and stored in the array precv, ! associated to a certain matrix A and stored in the array p%precv,
! - op(M^(-1)) is M^(-1) or its transpose, according to the value of trans, ! - op(M^(-1)) is M^(-1) or its transpose, according to the value of trans,
! - X and Y are vectors, ! - X and Y are vectors,
! - alpha and beta are scalars. ! - alpha and beta are scalars.
@ -781,7 +784,7 @@ contains
! !
! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners', ! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners',
! each representing the part of the preconditioner associated to a certain level. ! each representing the part of the preconditioner associated to a certain level.
! For each level ilev, the base preconditioner K(ilev) is stored in precv(ilev) ! For each level ilev, the base preconditioner K(ilev) is stored in p%precv(ilev)
! and is associated to a matrix A(ilev), obtained by 'tranferring' the original ! and is associated to a matrix A(ilev), obtained by 'tranferring' the original
! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed ! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed
! aggregation. ! aggregation.
@ -828,13 +831,14 @@ contains
! 5. Yext = beta*Yext + alpha*Y(1) ! 5. Yext = beta*Yext + alpha*Y(1)
! !
! !
subroutine mlt_post_ml_aply(alpha,precv,x,beta,y,desc_data,trans,work,info) subroutine mlt_post_ml_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
implicit none implicit none
! Arguments ! Arguments
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
type(mld_donelev_type), intent(in) :: precv(:) type(mld_dprec_type), intent(in) :: p
!!$ type(mld_donelev_type), intent(in) :: precv(:)
real(psb_dpk_),intent(in) :: alpha,beta real(psb_dpk_),intent(in) :: alpha,beta
real(psb_dpk_),intent(in) :: x(:) real(psb_dpk_),intent(in) :: x(:)
real(psb_dpk_),intent(inout) :: y(:) real(psb_dpk_),intent(inout) :: y(:)
@ -864,9 +868,9 @@ contains
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(precv) & ' Entry ', size(p%precv)
nlev = size(precv) nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info) allocate(mlprec_wrk(nlev),stat=info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
@ -882,7 +886,7 @@ contains
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' desc_data status',allocated(desc_data%matrix_data) & ' desc_data status',allocated(desc_data%matrix_data)
nc2l = psb_cd_get_local_cols(precv(1)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(1)%base_desc)
allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), & allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), &
& mlprec_wrk(1)%tx(nc2l), stat=info) & mlprec_wrk(1)%tx(nc2l), stat=info)
@ -894,9 +898,9 @@ contains
end if end if
call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%tx,& call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%tx,&
& precv(1)%base_desc,info) & p%precv(1)%base_desc,info)
call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%x2l,& call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%x2l,&
& precv(1)%base_desc,info) & p%precv(1)%base_desc,info)
! !
! STEP 2 ! STEP 2
@ -905,13 +909,13 @@ contains
! !
do ilev=2, nlev do ilev=2, nlev
nc2l = psb_cd_get_local_cols(precv(ilev)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(ilev)%base_desc)
nr2l = psb_cd_get_local_rows(precv(ilev)%base_desc) nr2l = psb_cd_get_local_rows(p%precv(ilev)%base_desc)
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name), & & write(debug_unit,*) me,' ',trim(name), &
& ' starting up sweep ',& & ' starting up sweep ',&
& ilev,allocated(precv(ilev)%iprcparm),nc2l, nr2l & ilev,allocated(p%precv(ilev)%iprcparm),nc2l, nr2l
allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%y2l(nc2l),& allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%y2l(nc2l),&
& mlprec_wrk(ilev)%x2l(nc2l), stat=info) & mlprec_wrk(ilev)%x2l(nc2l), stat=info)
@ -926,7 +930,7 @@ contains
! Apply prolongator transpose, i.e. restriction ! Apply prolongator transpose, i.e. restriction
call psb_map_X2Y(done,mlprec_wrk(ilev-1)%x2l,& call psb_map_X2Y(done,mlprec_wrk(ilev-1)%x2l,&
& dzero,mlprec_wrk(ilev)%x2l,& & dzero,mlprec_wrk(ilev)%x2l,&
& precv(ilev)%map,info,work=work) & p%precv(ilev)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction') call psb_errpush(4001,name,a_err='Error during restriction')
@ -937,7 +941,7 @@ contains
! update x2l ! update x2l
! !
call psb_geaxpby(done,mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%tx,& call psb_geaxpby(done,mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%tx,&
& precv(ilev)%base_desc,info) & p%precv(ilev)%base_desc,info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4001,name,a_err='Error in update') call psb_errpush(4001,name,a_err='Error in update')
goto 9999 goto 9999
@ -954,8 +958,8 @@ contains
! !
! Apply the base preconditioner at the coarsest level ! Apply the base preconditioner at the coarsest level
! !
call mld_baseprec_aply(done,precv(nlev)%prec,mlprec_wrk(nlev)%x2l, & call mld_baseprec_aply(done,p%precv(nlev)%prec,mlprec_wrk(nlev)%x2l, &
& dzero, mlprec_wrk(nlev)%y2l,precv(nlev)%base_desc,trans,work,info) & dzero, mlprec_wrk(nlev)%y2l,p%precv(nlev)%base_desc,trans,work,info)
if (info /=0) then if (info /=0) then
call psb_errpush(4010,name,a_err='baseprec_aply') call psb_errpush(4010,name,a_err='baseprec_aply')
goto 9999 goto 9999
@ -980,7 +984,7 @@ contains
! !
call psb_map_Y2X(done,mlprec_wrk(ilev+1)%y2l,& call psb_map_Y2X(done,mlprec_wrk(ilev+1)%y2l,&
& dzero,mlprec_wrk(ilev)%y2l,& & dzero,mlprec_wrk(ilev)%y2l,&
& precv(ilev+1)%map,info,work=work) & p%precv(ilev+1)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during prolongation') call psb_errpush(4001,name,a_err='Error during prolongation')
@ -990,15 +994,15 @@ contains
! !
! Compute the residual ! Compute the residual
! !
call psb_spmm(-done,precv(ilev)%base_a,mlprec_wrk(ilev)%y2l,& call psb_spmm(-done,p%precv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
& done,mlprec_wrk(ilev)%tx,precv(ilev)%base_desc,info,& & done,mlprec_wrk(ilev)%tx,p%precv(ilev)%base_desc,info,&
& work=work,trans=trans) & work=work,trans=trans)
! !
! Apply the base preconditioner ! Apply the base preconditioner
! !
if (info == 0) call mld_baseprec_aply(done,precv(ilev)%prec,& if (info == 0) call mld_baseprec_aply(done,p%precv(ilev)%prec,&
& mlprec_wrk(ilev)%tx,done,mlprec_wrk(ilev)%y2l,precv(ilev)%base_desc,& & mlprec_wrk(ilev)%tx,done,mlprec_wrk(ilev)%y2l,p%precv(ilev)%base_desc,&
& trans,work,info) & trans,work,info)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err=' spmm/baseprec_aply') call psb_errpush(4001,name,a_err=' spmm/baseprec_aply')
@ -1015,7 +1019,7 @@ contains
! !
! Compute the output vector Y ! Compute the output vector Y
! !
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,precv(1)%base_desc,info) call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,p%precv(1)%base_desc,info)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err=' Final update') call psb_errpush(4001,name,a_err=' Final update')
@ -1052,7 +1056,7 @@ contains
! where ! where
! - M is a symmetrized hybrid multilevel domain decomposition (Schwarz) ! - M is a symmetrized hybrid multilevel domain decomposition (Schwarz)
! preconditioner associated to a certain matrix A and stored in the array ! preconditioner associated to a certain matrix A and stored in the array
! precv, ! p%precv,
! - op(M^(-1)) is M^(-1) or its transpose, according to the value of trans, ! - op(M^(-1)) is M^(-1) or its transpose, according to the value of trans,
! - X and Y are vectors, ! - X and Y are vectors,
! - alpha and beta are scalars. ! - alpha and beta are scalars.
@ -1067,7 +1071,7 @@ contains
! !
! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners', ! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners',
! each representing the part of the preconditioner associated to a certain level. ! each representing the part of the preconditioner associated to a certain level.
! For each level ilev, the base preconditioner K(ilev) is stored in precv(ilev) ! For each level ilev, the base preconditioner K(ilev) is stored in p%precv(ilev)
! and is associated to a matrix A(ilev), obtained by 'tranferring' the original ! and is associated to a matrix A(ilev), obtained by 'tranferring' the original
! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed ! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed
! aggregation. ! aggregation.
@ -1125,13 +1129,14 @@ contains
! !
! 6. Yext = beta*Yext + alpha*Y(1) ! 6. Yext = beta*Yext + alpha*Y(1)
! !
subroutine mlt_twoside_ml_aply(alpha,precv,x,beta,y,desc_data,trans,work,info) subroutine mlt_twoside_ml_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
implicit none implicit none
! Arguments ! Arguments
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
type(mld_donelev_type), intent(in) :: precv(:) type(mld_dprec_type), intent(in) :: p
!!$ type(mld_donelev_type), intent(in) :: p%precv(:)
real(psb_dpk_),intent(in) :: alpha,beta real(psb_dpk_),intent(in) :: alpha,beta
real(psb_dpk_),intent(in) :: x(:) real(psb_dpk_),intent(in) :: x(:)
real(psb_dpk_),intent(inout) :: y(:) real(psb_dpk_),intent(inout) :: y(:)
@ -1161,9 +1166,9 @@ contains
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(precv) & ' Entry ', size(p%precv)
nlev = size(precv) nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info) allocate(mlprec_wrk(nlev),stat=info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
@ -1174,7 +1179,7 @@ contains
! !
! Copy the input vector X ! Copy the input vector X
! !
nc2l = psb_cd_get_local_cols(precv(1)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(1)%base_desc)
allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), & allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), &
& mlprec_wrk(1)%ty(nc2l), mlprec_wrk(1)%tx(nc2l), stat=info) & mlprec_wrk(1)%ty(nc2l), mlprec_wrk(1)%tx(nc2l), stat=info)
@ -1187,17 +1192,17 @@ contains
end if end if
call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%x2l,& call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%x2l,&
& precv(1)%base_desc,info) & p%precv(1)%base_desc,info)
call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%tx,& call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%tx,&
& precv(1)%base_desc,info) & p%precv(1)%base_desc,info)
! !
! STEP 2 ! STEP 2
! !
! Apply the base preconditioner at the finest level ! Apply the base preconditioner at the finest level
! !
call mld_baseprec_aply(done,precv(1)%prec,mlprec_wrk(1)%x2l,& call mld_baseprec_aply(done,p%precv(1)%prec,mlprec_wrk(1)%x2l,&
& dzero,mlprec_wrk(1)%y2l,precv(1)%base_desc,& & dzero,mlprec_wrk(1)%y2l,p%precv(1)%base_desc,&
& trans,work,info) & trans,work,info)
! !
! STEP 3 ! STEP 3
@ -1205,8 +1210,8 @@ contains
! Compute the residual at the finest level ! Compute the residual at the finest level
! !
mlprec_wrk(1)%ty = mlprec_wrk(1)%x2l mlprec_wrk(1)%ty = mlprec_wrk(1)%x2l
if (info == 0) call psb_spmm(-done,precv(1)%base_a,mlprec_wrk(1)%y2l,& if (info == 0) call psb_spmm(-done,p%precv(1)%base_a,mlprec_wrk(1)%y2l,&
& done,mlprec_wrk(1)%ty,precv(1)%base_desc,info,& & done,mlprec_wrk(1)%ty,p%precv(1)%base_desc,info,&
& work=work,trans=trans) & work=work,trans=trans)
if (info /=0) then if (info /=0) then
call psb_errpush(4010,name,a_err='Fine level baseprec/residual') call psb_errpush(4010,name,a_err='Fine level baseprec/residual')
@ -1220,8 +1225,8 @@ contains
! !
do ilev = 2, nlev do ilev = 2, nlev
nc2l = psb_cd_get_local_cols(precv(ilev)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(ilev)%base_desc)
nr2l = psb_cd_get_local_rows(precv(ilev)%base_desc) nr2l = psb_cd_get_local_rows(p%precv(ilev)%base_desc)
allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%ty(nc2l),& allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%ty(nc2l),&
& mlprec_wrk(ilev)%y2l(nc2l),mlprec_wrk(ilev)%x2l(nc2l), stat=info) & mlprec_wrk(ilev)%y2l(nc2l),mlprec_wrk(ilev)%x2l(nc2l), stat=info)
@ -1236,7 +1241,7 @@ contains
! Apply prolongator transpose, i.e. restriction ! Apply prolongator transpose, i.e. restriction
call psb_map_X2Y(done,mlprec_wrk(ilev-1)%ty,& call psb_map_X2Y(done,mlprec_wrk(ilev-1)%ty,&
& dzero,mlprec_wrk(ilev)%x2l,& & dzero,mlprec_wrk(ilev)%x2l,&
& precv(ilev)%map,info,work=work) & p%precv(ilev)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction') call psb_errpush(4001,name,a_err='Error during restriction')
@ -1244,21 +1249,21 @@ contains
end if end if
call psb_geaxpby(done,mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%tx,& call psb_geaxpby(done,mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%tx,&
& precv(ilev)%base_desc,info) & p%precv(ilev)%base_desc,info)
! !
! Apply the base preconditioner ! Apply the base preconditioner
! !
if (info == 0) call mld_baseprec_aply(done,precv(ilev)%prec,& if (info == 0) call mld_baseprec_aply(done,p%precv(ilev)%prec,&
& mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%y2l,& & mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%y2l,&
& precv(ilev)%base_desc,trans,work,info) & p%precv(ilev)%base_desc,trans,work,info)
! !
! Compute the residual (at all levels but the coarsest one) ! Compute the residual (at all levels but the coarsest one)
! !
if(ilev < nlev) then if(ilev < nlev) then
mlprec_wrk(ilev)%ty = mlprec_wrk(ilev)%x2l mlprec_wrk(ilev)%ty = mlprec_wrk(ilev)%x2l
if (info == 0) call psb_spmm(-done,precv(ilev)%base_a,& if (info == 0) call psb_spmm(-done,p%precv(ilev)%base_a,&
& mlprec_wrk(ilev)%y2l,done,mlprec_wrk(ilev)%ty,& & mlprec_wrk(ilev)%y2l,done,mlprec_wrk(ilev)%ty,&
& precv(ilev)%base_desc,info,work=work,trans=trans) & p%precv(ilev)%base_desc,info,work=work,trans=trans)
endif endif
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='baseprec_aply/residual') call psb_errpush(4001,name,a_err='baseprec_aply/residual')
@ -1279,7 +1284,7 @@ contains
! !
call psb_map_Y2X(done,mlprec_wrk(ilev+1)%y2l,& call psb_map_Y2X(done,mlprec_wrk(ilev+1)%y2l,&
& done,mlprec_wrk(ilev)%y2l,& & done,mlprec_wrk(ilev)%y2l,&
& precv(ilev+1)%map,info,work=work) & p%precv(ilev+1)%map,info,work=work)
if (info /=0 ) then if (info /=0 ) then
call psb_errpush(4001,name,a_err='Error during restriction') call psb_errpush(4001,name,a_err='Error during restriction')
@ -1289,14 +1294,14 @@ contains
! !
! Compute the residual ! Compute the residual
! !
call psb_spmm(-done,precv(ilev)%base_a,mlprec_wrk(ilev)%y2l,& call psb_spmm(-done,p%precv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
& done,mlprec_wrk(ilev)%tx,precv(ilev)%base_desc,info,& & done,mlprec_wrk(ilev)%tx,p%precv(ilev)%base_desc,info,&
& work=work,trans=trans) & work=work,trans=trans)
! !
! Apply the base preconditioner ! Apply the base preconditioner
! !
if (info == 0) call mld_baseprec_aply(done,precv(ilev)%prec,mlprec_wrk(ilev)%tx,& if (info == 0) call mld_baseprec_aply(done,p%precv(ilev)%prec,mlprec_wrk(ilev)%tx,&
& done,mlprec_wrk(ilev)%y2l,precv(ilev)%base_desc, trans, work,info) & done,mlprec_wrk(ilev)%y2l,p%precv(ilev)%base_desc, trans, work,info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4001,name,a_err='Error: residual/baseprec_aply') call psb_errpush(4001,name,a_err='Error: residual/baseprec_aply')
goto 9999 goto 9999
@ -1309,7 +1314,7 @@ contains
! Compute the output vector Y ! Compute the output vector Y
! !
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,& call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,&
& precv(1)%base_desc,info) & p%precv(1)%base_desc,info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4001,name,a_err='Error final update') call psb_errpush(4001,name,a_err='Error final update')

2
mlprec/mld_dprecaply.f90
Unescape Escape View File

@ -127,7 +127,7 @@ subroutine mld_dprecaply(prec,x,y,desc_data,info,trans,work)
goto 9999 goto 9999
end if end if
if (size(prec%precv) >1) then if (size(prec%precv) >1) then
call mld_mlprec_aply(done,prec%precv,x,dzero,y,desc_data,trans_,work_,info) call mld_mlprec_aply(done,prec,x,dzero,y,desc_data,trans_,work_,info)
if(info /= 0) then if(info /= 0) then
call psb_errpush(4010,name,a_err='mld_dmlprec_aply') call psb_errpush(4010,name,a_err='mld_dmlprec_aply')
goto 9999 goto 9999

24
mlprec/mld_inner_mod.f90
Unescape Escape View File

@ -186,11 +186,11 @@ module mld_inner_mod
end interface end interface
interface mld_mlprec_aply interface mld_mlprec_aply
subroutine mld_smlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info) subroutine mld_smlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_ use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_
use mld_prec_type, only : mld_sbaseprec_type, mld_sonelev_type use mld_prec_type, only : mld_sbaseprec_type, mld_sprec_type
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
type(mld_sonelev_type), intent(in) :: precv(:) type(mld_sprec_type), intent(in) :: p
real(psb_spk_),intent(in) :: alpha,beta real(psb_spk_),intent(in) :: alpha,beta
real(psb_spk_),intent(in) :: x(:) real(psb_spk_),intent(in) :: x(:)
real(psb_spk_),intent(inout) :: y(:) real(psb_spk_),intent(inout) :: y(:)
@ -198,11 +198,11 @@ module mld_inner_mod
real(psb_spk_),target :: work(:) real(psb_spk_),target :: work(:)
integer, intent(out) :: info integer, intent(out) :: info
end subroutine mld_smlprec_aply end subroutine mld_smlprec_aply
subroutine mld_dmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info) subroutine mld_dmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_ use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_
use mld_prec_type, only : mld_dbaseprec_type, mld_donelev_type use mld_prec_type, only : mld_dbaseprec_type, mld_dprec_type
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
type(mld_donelev_type), intent(in) :: precv(:) type(mld_dprec_type), intent(in) :: p
real(psb_dpk_),intent(in) :: alpha,beta real(psb_dpk_),intent(in) :: alpha,beta
real(psb_dpk_),intent(in) :: x(:) real(psb_dpk_),intent(in) :: x(:)
real(psb_dpk_),intent(inout) :: y(:) real(psb_dpk_),intent(inout) :: y(:)
@ -210,11 +210,11 @@ module mld_inner_mod
real(psb_dpk_),target :: work(:) real(psb_dpk_),target :: work(:)
integer, intent(out) :: info integer, intent(out) :: info
end subroutine mld_dmlprec_aply end subroutine mld_dmlprec_aply
subroutine mld_cmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) subroutine mld_cmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_ use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_
use mld_prec_type, only : mld_cbaseprec_type, mld_conelev_type use mld_prec_type, only : mld_cbaseprec_type, mld_cprec_type
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
type(mld_conelev_type), intent(in) :: baseprecv(:) type(mld_cprec_type), intent(in) :: p
complex(psb_spk_),intent(in) :: alpha,beta complex(psb_spk_),intent(in) :: alpha,beta
complex(psb_spk_),intent(in) :: x(:) complex(psb_spk_),intent(in) :: x(:)
complex(psb_spk_),intent(inout) :: y(:) complex(psb_spk_),intent(inout) :: y(:)
@ -222,11 +222,11 @@ module mld_inner_mod
complex(psb_spk_),target :: work(:) complex(psb_spk_),target :: work(:)
integer, intent(out) :: info integer, intent(out) :: info
end subroutine mld_cmlprec_aply end subroutine mld_cmlprec_aply
subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) subroutine mld_zmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_ use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_
use mld_prec_type, only : mld_zbaseprec_type, mld_zonelev_type use mld_prec_type, only : mld_zbaseprec_type, mld_zprec_type
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
type(mld_zonelev_type), intent(in) :: baseprecv(:) type(mld_zprec_type), intent(in) :: p
complex(psb_dpk_),intent(in) :: alpha,beta complex(psb_dpk_),intent(in) :: alpha,beta
complex(psb_dpk_),intent(in) :: x(:) complex(psb_dpk_),intent(in) :: x(:)
complex(psb_dpk_),intent(inout) :: y(:) complex(psb_dpk_),intent(inout) :: y(:)

228
mlprec/mld_smlprec_aply.f90
Unescape Escape View File

@ -46,7 +46,7 @@
! Y = beta*Y + alpha*op(M^(-1))*X, ! Y = beta*Y + alpha*op(M^(-1))*X,
! where ! where
! - M is a multilevel domain decomposition (Schwarz) preconditioner associated ! - M is a multilevel domain decomposition (Schwarz) preconditioner associated
! to a certain matrix A and stored in the array precv, ! to a certain matrix A and stored in the array p%precv,
! - op(M^(-1)) is M^(-1) or its transpose, according to the value of trans, ! - op(M^(-1)) is M^(-1) or its transpose, according to the value of trans,
! - X and Y are vectors, ! - X and Y are vectors,
! - alpha and beta are scalars. ! - alpha and beta are scalars.
@ -57,7 +57,7 @@
! !
! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners', ! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners',
! each representing the part of the preconditioner associated to a certain level. ! each representing the part of the preconditioner associated to a certain level.
! For each level ilev, the preconditioner K(ilev) is stored in precv(ilev) ! For each level ilev, the preconditioner K(ilev) is stored in p%precv(ilev)
! and is associated to a matrix A(ilev), obtained by 'tranferring' the original ! and is associated to a matrix A(ilev), obtained by 'tranferring' the original
! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed ! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed
! aggregation. ! aggregation.
@ -78,41 +78,41 @@
! Arguments: ! Arguments:
! alpha - real(psb_spk_), input. ! alpha - real(psb_spk_), input.
! The scalar alpha. ! The scalar alpha.
! precv - type(mld_sonelev_type), dimension(:), input. ! p - type(mld_sprec_type), input.
! The array of one-level preconditioner data structures containing the ! The array of one-level preconditioner data structures containing the
! local parts of the preconditioners to be applied at each level. ! local parts of the preconditioners to be applied at each level.
! Note that nlev = size(precv) = number of levels. ! Note that nlev = size(p%precv) = number of levels.
! precv(ilev)%prec - type(psb_sbaseprec_type) ! p%precv(ilev)%prec - type(psb_sbaseprec_type)
! The "base" preconditioner for the current level ! The "base" preconditioner for the current level
! precv(ilev)%ac - type(psb_sspmat_type) ! p%precv(ilev)%ac - type(psb_sspmat_type)
! The local part of the matrix A(ilev). ! The local part of the matrix A(ilev).
! precv(ilev)%desc_ac - type(psb_desc_type). ! p%precv(ilev)%desc_ac - type(psb_desc_type).
! The communication descriptor associated to the sparse ! The communication descriptor associated to the sparse
! matrix A(ilev) ! matrix A(ilev)
! precv(ilev)%map - type(psb_inter_desc_type) ! p%precv(ilev)%map - type(psb_inter_desc_type)
! Stores the linear operators mapping between levels ! Stores the linear operators mapping between levels
! (ilev-1) and (ilev). These are the restriction and ! (ilev-1) and (ilev). These are the restriction and
! prolongation operators described in the sequel. ! prolongation operators described in the sequel.
! precv(ilev)%iprcparm - integer, dimension(:), allocatable. ! p%precv(ilev)%iprcparm - integer, dimension(:), allocatable.
! The integer parameters defining the multilevel ! The integer parameters defining the multilevel
! strategy ! strategy
! precv(ilev)%rprcparm - real(psb_spk_), dimension(:), allocatable. ! p%precv(ilev)%rprcparm - real(psb_spk_), dimension(:), allocatable.
! The real parameters defining the multilevel strategy ! The real parameters defining the multilevel strategy
! precv(ilev)%mlia - integer, dimension(:), allocatable. ! p%precv(ilev)%mlia - integer, dimension(:), allocatable.
! The aggregation map (ilev-1) --> (ilev). ! The aggregation map (ilev-1) --> (ilev).
! In case of non-smoothed aggregation, it is used ! In case of non-smoothed aggregation, it is used
! instead of mld_sm_pr_. ! instead of mld_sm_pr_.
! precv(ilev)%nlaggr - integer, dimension(:), allocatable. ! p%precv(ilev)%nlaggr - integer, dimension(:), allocatable.
! The number of aggregates (rows of A(ilev)) on the ! The number of aggregates (rows of A(ilev)) on the
! various processes. ! various processes.
! precv(ilev)%base_a - type(psb_sspmat_type), pointer. ! p%precv(ilev)%base_a - type(psb_sspmat_type), pointer.
! Pointer (really a pointer!) to the base matrix of ! Pointer (really a pointer!) to the base matrix of
! the current level, i.e. the local part of A(ilev); ! the current level, i.e. the local part of A(ilev);
! so we have a unified treatment of residuals. We ! so we have a unified treatment of residuals. We
! need this to avoid passing explicitly the matrix ! need this to avoid passing explicitly the matrix
! A(ilev) to the routine which applies the ! A(ilev) to the routine which applies the
! preconditioner. ! preconditioner.
! precv(ilev)%base_desc - type(psb_desc_type), pointer. ! p%precv(ilev)%base_desc - type(psb_desc_type), pointer.
! Pointer to the communication descriptor associated ! Pointer to the communication descriptor associated
! to the sparse matrix pointed by base_a. ! to the sparse matrix pointed by base_a.
! !
@ -136,10 +136,10 @@
! Note that when the LU factorization of the matrix A(ilev) is computed instead of ! Note that when the LU factorization of the matrix A(ilev) is computed instead of
! the ILU one, by using UMFPACK or SuperLU, the corresponding L and U factors ! the ILU one, by using UMFPACK or SuperLU, the corresponding L and U factors
! are stored in data structures provided by UMFPACK or SuperLU and pointed by ! are stored in data structures provided by UMFPACK or SuperLU and pointed by
! precv(ilev)%prec%iprcparm(mld_umf_ptr) or precv(ilev)%prec%iprcparm(mld_slu_ptr), ! p%precv(ilev)%prec%iprcparm(mld_umf_ptr) or p%precv(ilev)%prec%iprcparm(mld_slu_ptr),
! respectively. ! respectively.
! !
subroutine mld_smlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info) subroutine mld_smlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
use psb_base_mod use psb_base_mod
use mld_inner_mod, mld_protect_name => mld_smlprec_aply use mld_inner_mod, mld_protect_name => mld_smlprec_aply
@ -148,7 +148,7 @@ subroutine mld_smlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
! Arguments ! Arguments
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
type(mld_sonelev_type), intent(in) :: precv(:) type(mld_sprec_type), intent(in) :: p
real(psb_spk_),intent(in) :: alpha,beta real(psb_spk_),intent(in) :: alpha,beta
real(psb_spk_),intent(in) :: x(:) real(psb_spk_),intent(in) :: x(:)
real(psb_spk_),intent(inout) :: y(:) real(psb_spk_),intent(inout) :: y(:)
@ -173,11 +173,11 @@ subroutine mld_smlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(precv) & ' Entry ', size(p%precv)
trans_ = psb_toupper(trans) trans_ = psb_toupper(trans)
select case(precv(2)%iprcparm(mld_ml_type_)) select case(p%precv(2)%iprcparm(mld_ml_type_))
case(mld_no_ml_) case(mld_no_ml_)
! !
@ -191,7 +191,7 @@ subroutine mld_smlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
! Additive multilevel ! Additive multilevel
! !
call add_ml_aply(alpha,precv,x,beta,y,desc_data,trans_,work,info) call add_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info)
case(mld_mult_ml_) case(mld_mult_ml_)
! !
@ -202,15 +202,15 @@ subroutine mld_smlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
! Note that the transpose switches pre <-> post. ! Note that the transpose switches pre <-> post.
! !
select case(precv(2)%iprcparm(mld_smoother_pos_)) select case(p%precv(2)%iprcparm(mld_smoother_pos_))
case(mld_post_smooth_) case(mld_post_smooth_)
select case (trans_) select case (trans_)
case('N') case('N')
call mlt_post_ml_aply(alpha,precv,x,beta,y,desc_data,trans_,work,info) call mlt_post_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info)
case('T','C') case('T','C')
call mlt_pre_ml_aply(alpha,precv,x,beta,y,desc_data,trans_,work,info) call mlt_pre_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info)
case default case default
info = 4001 info = 4001
call psb_errpush(info,name,a_err='invalid trans') call psb_errpush(info,name,a_err='invalid trans')
@ -221,9 +221,9 @@ subroutine mld_smlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
select case (trans_) select case (trans_)
case('N') case('N')
call mlt_pre_ml_aply(alpha,precv,x,beta,y,desc_data,trans_,work,info) call mlt_pre_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info)
case('T','C') case('T','C')
call mlt_post_ml_aply(alpha,precv,x,beta,y,desc_data,trans_,work,info) call mlt_post_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info)
case default case default
info = 4001 info = 4001
call psb_errpush(info,name,a_err='invalid trans') call psb_errpush(info,name,a_err='invalid trans')
@ -232,12 +232,12 @@ subroutine mld_smlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
case(mld_twoside_smooth_) case(mld_twoside_smooth_)
call mlt_twoside_ml_aply(alpha,precv,x,beta,y,desc_data,trans_,work,info) call mlt_twoside_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info)
case default case default
info = 4013 info = 4013
call psb_errpush(info,name,a_err='invalid smooth_pos',& call psb_errpush(info,name,a_err='invalid smooth_pos',&
& i_Err=(/precv(2)%iprcparm(mld_smoother_pos_),0,0,0,0/)) & i_Err=(/p%precv(2)%iprcparm(mld_smoother_pos_),0,0,0,0/))
goto 9999 goto 9999
end select end select
@ -245,7 +245,7 @@ subroutine mld_smlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
case default case default
info = 4013 info = 4013
call psb_errpush(info,name,a_err='invalid mltype',& call psb_errpush(info,name,a_err='invalid mltype',&
& i_Err=(/precv(2)%iprcparm(mld_ml_type_),0,0,0,0/)) & i_Err=(/p%precv(2)%iprcparm(mld_ml_type_),0,0,0,0/))
goto 9999 goto 9999
end select end select
@ -272,7 +272,7 @@ contains
! Y = beta*Y + alpha*op(M^(-1))*X, ! Y = beta*Y + alpha*op(M^(-1))*X,
! where ! where
! - M is an additive multilevel domain decomposition (Schwarz) preconditioner ! - M is an additive multilevel domain decomposition (Schwarz) preconditioner
! associated to a certain matrix A and stored in the array precv, ! associated to a certain matrix A and stored in the array p%precv,
! - op(M^(-1)) is M^(-1) or its transpose, according to the value of trans, ! - op(M^(-1)) is M^(-1) or its transpose, according to the value of trans,
! - X and Y are vectors, ! - X and Y are vectors,
! - alpha and beta are scalars. ! - alpha and beta are scalars.
@ -286,7 +286,7 @@ contains
! !
! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners', ! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners',
! each representing the part of the preconditioner associated to a certain level. ! each representing the part of the preconditioner associated to a certain level.
! For each level ilev, the base preconditioner K(ilev) is stored in precv(ilev) ! For each level ilev, the base preconditioner K(ilev) is stored in p%precv(ilev)
! and is associated to a matrix A(ilev), obtained by 'tranferring' the original ! and is associated to a matrix A(ilev), obtained by 'tranferring' the original
! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed ! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed
! aggregation. ! aggregation.
@ -334,13 +334,13 @@ contains
! !
! 4. Yext = beta*Yext + alpha*Y(1) ! 4. Yext = beta*Yext + alpha*Y(1)
! !
subroutine add_ml_aply(alpha,precv,x,beta,y,desc_data,trans,work,info) subroutine add_ml_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
implicit none implicit none
! Arguments ! Arguments
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
type(mld_sonelev_type), intent(in) :: precv(:) type(mld_sprec_type), intent(in) :: p
real(psb_spk_),intent(in) :: alpha,beta real(psb_spk_),intent(in) :: alpha,beta
real(psb_spk_),intent(in) :: x(:) real(psb_spk_),intent(in) :: x(:)
real(psb_spk_),intent(inout) :: y(:) real(psb_spk_),intent(inout) :: y(:)
@ -370,9 +370,9 @@ contains
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(precv) & ' Entry ', size(p%precv)
nlev = size(precv) nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info) allocate(mlprec_wrk(nlev),stat=info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
@ -395,8 +395,8 @@ contains
mlprec_wrk(1)%x2l(:) = x(:) mlprec_wrk(1)%x2l(:) = x(:)
mlprec_wrk(1)%y2l(:) = szero mlprec_wrk(1)%y2l(:) = szero
call mld_baseprec_aply(alpha,precv(1)%prec,x,beta,y,& call mld_baseprec_aply(alpha,p%precv(1)%prec,x,beta,y,&
& precv(1)%base_desc,trans,work,info) & p%precv(1)%base_desc,trans,work,info)
if (info /=0) then if (info /=0) then
call psb_errpush(4010,name,a_err='baseprec_aply') call psb_errpush(4010,name,a_err='baseprec_aply')
goto 9999 goto 9999
@ -407,8 +407,8 @@ contains
! For each level except the finest one ... ! For each level except the finest one ...
! !
do ilev = 2, nlev do ilev = 2, nlev
nc2l = psb_cd_get_local_cols(precv(ilev)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(ilev)%base_desc)
nr2l = psb_cd_get_local_rows(precv(ilev)%base_desc) nr2l = psb_cd_get_local_rows(p%precv(ilev)%base_desc)
allocate(mlprec_wrk(ilev)%x2l(nc2l),mlprec_wrk(ilev)%y2l(nc2l),& allocate(mlprec_wrk(ilev)%x2l(nc2l),mlprec_wrk(ilev)%y2l(nc2l),&
& stat=info) & stat=info)
if (info /= 0) then if (info /= 0) then
@ -421,7 +421,7 @@ contains
! Apply prolongator transpose, i.e. restriction ! Apply prolongator transpose, i.e. restriction
call psb_map_X2Y(sone,mlprec_wrk(ilev-1)%x2l,& call psb_map_X2Y(sone,mlprec_wrk(ilev-1)%x2l,&
& szero,mlprec_wrk(ilev)%x2l,& & szero,mlprec_wrk(ilev)%x2l,&
& precv(ilev)%map,info,work=work) & p%precv(ilev)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction') call psb_errpush(4001,name,a_err='Error during restriction')
@ -431,9 +431,9 @@ contains
! !
! Apply the base preconditioner ! Apply the base preconditioner
! !
call mld_baseprec_aply(sone,precv(ilev)%prec,& call mld_baseprec_aply(sone,p%precv(ilev)%prec,&
& mlprec_wrk(ilev)%x2l,szero,mlprec_wrk(ilev)%y2l,& & mlprec_wrk(ilev)%x2l,szero,mlprec_wrk(ilev)%y2l,&
& precv(ilev)%base_desc, trans,work,info) & p%precv(ilev)%base_desc, trans,work,info)
enddo enddo
@ -444,15 +444,15 @@ contains
! !
do ilev =nlev,2,-1 do ilev =nlev,2,-1
nc2l = psb_cd_get_local_cols(precv(ilev)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(ilev)%base_desc)
nr2l = psb_cd_get_local_rows(precv(ilev)%base_desc) nr2l = psb_cd_get_local_rows(p%precv(ilev)%base_desc)
! !
! Apply prolongator ! Apply prolongator
! !
call psb_map_Y2X(sone,mlprec_wrk(ilev)%y2l,& call psb_map_Y2X(sone,mlprec_wrk(ilev)%y2l,&
& sone,mlprec_wrk(ilev-1)%y2l,& & sone,mlprec_wrk(ilev-1)%y2l,&
& precv(ilev)%map,info,work=work) & p%precv(ilev)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during prolongation') call psb_errpush(4001,name,a_err='Error during prolongation')
@ -465,7 +465,7 @@ contains
! !
! Compute the output vector Y ! Compute the output vector Y
! !
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,sone,y,precv(1)%base_desc,info) call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,sone,y,p%precv(1)%base_desc,info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4001,name,a_err='Error on final update') call psb_errpush(4001,name,a_err='Error on final update')
goto 9999 goto 9999
@ -499,7 +499,7 @@ contains
! Y = beta*Y + alpha*op(M^(-1))*X, ! Y = beta*Y + alpha*op(M^(-1))*X,
! where ! where
! - M is a hybrid multilevel domain decomposition (Schwarz) preconditioner ! - M is a hybrid multilevel domain decomposition (Schwarz) preconditioner
! associated to a certain matrix A and stored in the array precv, ! associated to a certain matrix A and stored in the array p%precv,
! - op(M^(-1)) is M^(-1) or its transpose, according to the value of trans, ! - op(M^(-1)) is M^(-1) or its transpose, according to the value of trans,
! - X and Y are vectors, ! - X and Y are vectors,
! - alpha and beta are scalars. ! - alpha and beta are scalars.
@ -513,7 +513,7 @@ contains
! !
! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners', ! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners',
! each representing the part of the preconditioner associated to a certain level. ! each representing the part of the preconditioner associated to a certain level.
! For each level ilev, the base preconditioner K(ilev) is stored in precv(ilev) ! For each level ilev, the base preconditioner K(ilev) is stored in p%precv(ilev)
! and is associated to a matrix A(ilev), obtained by 'tranferring' the original ! and is associated to a matrix A(ilev), obtained by 'tranferring' the original
! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed ! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed
! aggregation. ! aggregation.
@ -569,13 +569,13 @@ contains
! 6. Yext = beta*Yext + alpha*Y(1) ! 6. Yext = beta*Yext + alpha*Y(1)
! !
! !
subroutine mlt_pre_ml_aply(alpha,precv,x,beta,y,desc_data,trans,work,info) subroutine mlt_pre_ml_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
implicit none implicit none
! Arguments ! Arguments
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
type(mld_sonelev_type), intent(in) :: precv(:) type(mld_sprec_type), intent(in) :: p
real(psb_spk_),intent(in) :: alpha,beta real(psb_spk_),intent(in) :: alpha,beta
real(psb_spk_),intent(in) :: x(:) real(psb_spk_),intent(in) :: x(:)
real(psb_spk_),intent(inout) :: y(:) real(psb_spk_),intent(inout) :: y(:)
@ -605,9 +605,9 @@ contains
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(precv) & ' Entry ', size(p%precv)
nlev = size(precv) nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info) allocate(mlprec_wrk(nlev),stat=info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
@ -619,7 +619,7 @@ contains
! !
! Copy the input vector X ! Copy the input vector X
! !
nc2l = psb_cd_get_local_cols(precv(1)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(1)%base_desc)
allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), & allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), &
& mlprec_wrk(1)%tx(nc2l), stat=info) & mlprec_wrk(1)%tx(nc2l), stat=info)
@ -636,8 +636,8 @@ contains
! !
! Apply the base preconditioner at the finest level ! Apply the base preconditioner at the finest level
! !
call mld_baseprec_aply(sone,precv(1)%prec,mlprec_wrk(1)%x2l,& call mld_baseprec_aply(sone,p%precv(1)%prec,mlprec_wrk(1)%x2l,&
& szero,mlprec_wrk(1)%y2l,precv(1)%base_desc,& & szero,mlprec_wrk(1)%y2l,p%precv(1)%base_desc,&
& trans,work,info) & trans,work,info)
if (info /=0) then if (info /=0) then
call psb_errpush(4010,name,a_err=' baseprec_aply') call psb_errpush(4010,name,a_err=' baseprec_aply')
@ -651,8 +651,8 @@ contains
! !
mlprec_wrk(1)%tx = mlprec_wrk(1)%x2l mlprec_wrk(1)%tx = mlprec_wrk(1)%x2l
call psb_spmm(-sone,precv(1)%base_a,mlprec_wrk(1)%y2l,& call psb_spmm(-sone,p%precv(1)%base_a,mlprec_wrk(1)%y2l,&
& sone,mlprec_wrk(1)%tx,precv(1)%base_desc,info,& & sone,mlprec_wrk(1)%tx,p%precv(1)%base_desc,info,&
& work=work,trans=trans) & work=work,trans=trans)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err=' fine level residual') call psb_errpush(4001,name,a_err=' fine level residual')
@ -666,8 +666,8 @@ contains
! !
do ilev = 2, nlev do ilev = 2, nlev
nc2l = psb_cd_get_local_cols(precv(ilev)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(ilev)%base_desc)
nr2l = psb_cd_get_local_rows(precv(ilev)%base_desc) nr2l = psb_cd_get_local_rows(p%precv(ilev)%base_desc)
allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%y2l(nc2l),& allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%y2l(nc2l),&
& mlprec_wrk(ilev)%x2l(nc2l), stat=info) & mlprec_wrk(ilev)%x2l(nc2l), stat=info)
@ -681,7 +681,7 @@ contains
! Apply prolongator transpose, i.e. restriction ! Apply prolongator transpose, i.e. restriction
call psb_map_X2Y(sone,mlprec_wrk(ilev-1)%tx,& call psb_map_X2Y(sone,mlprec_wrk(ilev-1)%tx,&
& szero,mlprec_wrk(ilev)%x2l,& & szero,mlprec_wrk(ilev)%x2l,&
& precv(ilev)%map,info,work=work) & p%precv(ilev)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction') call psb_errpush(4001,name,a_err='Error during restriction')
@ -691,17 +691,17 @@ contains
! !
! Apply the base preconditioner ! Apply the base preconditioner
! !
call mld_baseprec_aply(sone,precv(ilev)%prec,mlprec_wrk(ilev)%x2l,& call mld_baseprec_aply(sone,p%precv(ilev)%prec,mlprec_wrk(ilev)%x2l,&
& szero,mlprec_wrk(ilev)%y2l,precv(ilev)%base_desc,trans,work,info) & szero,mlprec_wrk(ilev)%y2l,p%precv(ilev)%base_desc,trans,work,info)
! !
! Compute the residual (at all levels but the coarsest one) ! Compute the residual (at all levels but the coarsest one)
! !
if (ilev < nlev) then if (ilev < nlev) then
mlprec_wrk(ilev)%tx = mlprec_wrk(ilev)%x2l mlprec_wrk(ilev)%tx = mlprec_wrk(ilev)%x2l
if (info == 0) call psb_spmm(-sone,precv(ilev)%base_a,& if (info == 0) call psb_spmm(-sone,p%precv(ilev)%base_a,&
& mlprec_wrk(ilev)%y2l,sone,mlprec_wrk(ilev)%tx,& & mlprec_wrk(ilev)%y2l,sone,mlprec_wrk(ilev)%tx,&
& precv(ilev)%base_desc,info,work=work,trans=trans) & p%precv(ilev)%base_desc,info,work=work,trans=trans)
endif endif
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error on up sweep residual') call psb_errpush(4001,name,a_err='Error on up sweep residual')
@ -720,7 +720,7 @@ contains
! !
call psb_map_Y2X(sone,mlprec_wrk(ilev+1)%y2l,& call psb_map_Y2X(sone,mlprec_wrk(ilev+1)%y2l,&
& sone,mlprec_wrk(ilev)%y2l,& & sone,mlprec_wrk(ilev)%y2l,&
& precv(ilev+1)%map,info,work=work) & p%precv(ilev+1)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during prolongation') call psb_errpush(4001,name,a_err='Error during prolongation')
@ -734,7 +734,7 @@ contains
! Compute the output vector Y ! Compute the output vector Y
! !
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,& call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,&
& precv(1)%base_desc,info) & p%precv(1)%base_desc,info)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error on final update') call psb_errpush(4001,name,a_err='Error on final update')
goto 9999 goto 9999
@ -767,7 +767,7 @@ contains
! Y = beta*Y + alpha*op(M^(-1))*X, ! Y = beta*Y + alpha*op(M^(-1))*X,
! where ! where
! - M is a hybrid multilevel domain decomposition (Schwarz) preconditioner ! - M is a hybrid multilevel domain decomposition (Schwarz) preconditioner
! associated to a certain matrix A and stored in the array precv, ! associated to a certain matrix A and stored in the array p%precv,
! - op(M^(-1)) is M^(-1) or its transpose, according to the value of trans, ! - op(M^(-1)) is M^(-1) or its transpose, according to the value of trans,
! - X and Y are vectors, ! - X and Y are vectors,
! - alpha and beta are scalars. ! - alpha and beta are scalars.
@ -781,7 +781,7 @@ contains
! !
! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners', ! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners',
! each representing the part of the preconditioner associated to a certain level. ! each representing the part of the preconditioner associated to a certain level.
! For each level ilev, the base preconditioner K(ilev) is stored in precv(ilev) ! For each level ilev, the base preconditioner K(ilev) is stored in p%precv(ilev)
! and is associated to a matrix A(ilev), obtained by 'tranferring' the original ! and is associated to a matrix A(ilev), obtained by 'tranferring' the original
! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed ! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed
! aggregation. ! aggregation.
@ -828,13 +828,13 @@ contains
! 5. Yext = beta*Yext + alpha*Y(1) ! 5. Yext = beta*Yext + alpha*Y(1)
! !
! !
subroutine mlt_post_ml_aply(alpha,precv,x,beta,y,desc_data,trans,work,info) subroutine mlt_post_ml_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
implicit none implicit none
! Arguments ! Arguments
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
type(mld_sonelev_type), intent(in) :: precv(:) type(mld_sprec_type), intent(in) :: p
real(psb_spk_),intent(in) :: alpha,beta real(psb_spk_),intent(in) :: alpha,beta
real(psb_spk_),intent(in) :: x(:) real(psb_spk_),intent(in) :: x(:)
real(psb_spk_),intent(inout) :: y(:) real(psb_spk_),intent(inout) :: y(:)
@ -864,9 +864,9 @@ contains
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(precv) & ' Entry ', size(p%precv)
nlev = size(precv) nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info) allocate(mlprec_wrk(nlev),stat=info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
@ -882,7 +882,7 @@ contains
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' desc_data status',allocated(desc_data%matrix_data) & ' desc_data status',allocated(desc_data%matrix_data)
nc2l = psb_cd_get_local_cols(precv(1)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(1)%base_desc)
allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), & allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), &
& mlprec_wrk(1)%tx(nc2l), stat=info) & mlprec_wrk(1)%tx(nc2l), stat=info)
@ -894,9 +894,9 @@ contains
end if end if
call psb_geaxpby(sone,x,szero,mlprec_wrk(1)%tx,& call psb_geaxpby(sone,x,szero,mlprec_wrk(1)%tx,&
& precv(1)%base_desc,info) & p%precv(1)%base_desc,info)
call psb_geaxpby(sone,x,szero,mlprec_wrk(1)%x2l,& call psb_geaxpby(sone,x,szero,mlprec_wrk(1)%x2l,&
& precv(1)%base_desc,info) & p%precv(1)%base_desc,info)
! !
! STEP 2 ! STEP 2
@ -905,13 +905,13 @@ contains
! !
do ilev=2, nlev do ilev=2, nlev
nc2l = psb_cd_get_local_cols(precv(ilev)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(ilev)%base_desc)
nr2l = psb_cd_get_local_rows(precv(ilev)%base_desc) nr2l = psb_cd_get_local_rows(p%precv(ilev)%base_desc)
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name), & & write(debug_unit,*) me,' ',trim(name), &
& ' starting up sweep ',& & ' starting up sweep ',&
& ilev,allocated(precv(ilev)%iprcparm),nc2l, nr2l & ilev,allocated(p%precv(ilev)%iprcparm),nc2l, nr2l
allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%y2l(nc2l),& allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%y2l(nc2l),&
& mlprec_wrk(ilev)%x2l(nc2l), stat=info) & mlprec_wrk(ilev)%x2l(nc2l), stat=info)
@ -926,7 +926,7 @@ contains
! Apply prolongator transpose, i.e. restriction ! Apply prolongator transpose, i.e. restriction
call psb_map_X2Y(sone,mlprec_wrk(ilev-1)%x2l,& call psb_map_X2Y(sone,mlprec_wrk(ilev-1)%x2l,&
& szero,mlprec_wrk(ilev)%x2l,& & szero,mlprec_wrk(ilev)%x2l,&
& precv(ilev)%map,info,work=work) & p%precv(ilev)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction') call psb_errpush(4001,name,a_err='Error during restriction')
@ -937,7 +937,7 @@ contains
! update x2l ! update x2l
! !
call psb_geaxpby(sone,mlprec_wrk(ilev)%x2l,szero,mlprec_wrk(ilev)%tx,& call psb_geaxpby(sone,mlprec_wrk(ilev)%x2l,szero,mlprec_wrk(ilev)%tx,&
& precv(ilev)%base_desc,info) & p%precv(ilev)%base_desc,info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4001,name,a_err='Error in update') call psb_errpush(4001,name,a_err='Error in update')
goto 9999 goto 9999
@ -954,8 +954,8 @@ contains
! !
! Apply the base preconditioner at the coarsest level ! Apply the base preconditioner at the coarsest level
! !
call mld_baseprec_aply(sone,precv(nlev)%prec,mlprec_wrk(nlev)%x2l, & call mld_baseprec_aply(sone,p%precv(nlev)%prec,mlprec_wrk(nlev)%x2l, &
& szero, mlprec_wrk(nlev)%y2l,precv(nlev)%base_desc,trans,work,info) & szero, mlprec_wrk(nlev)%y2l,p%precv(nlev)%base_desc,trans,work,info)
if (info /=0) then if (info /=0) then
call psb_errpush(4010,name,a_err='baseprec_aply') call psb_errpush(4010,name,a_err='baseprec_aply')
goto 9999 goto 9999
@ -980,7 +980,7 @@ contains
! !
call psb_map_Y2X(sone,mlprec_wrk(ilev+1)%y2l,& call psb_map_Y2X(sone,mlprec_wrk(ilev+1)%y2l,&
& szero,mlprec_wrk(ilev)%y2l,& & szero,mlprec_wrk(ilev)%y2l,&
& precv(ilev+1)%map,info,work=work) & p%precv(ilev+1)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during prolongation') call psb_errpush(4001,name,a_err='Error during prolongation')
@ -990,15 +990,15 @@ contains
! !
! Compute the residual ! Compute the residual
! !
call psb_spmm(-sone,precv(ilev)%base_a,mlprec_wrk(ilev)%y2l,& call psb_spmm(-sone,p%precv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
& sone,mlprec_wrk(ilev)%tx,precv(ilev)%base_desc,info,& & sone,mlprec_wrk(ilev)%tx,p%precv(ilev)%base_desc,info,&
& work=work,trans=trans) & work=work,trans=trans)
! !
! Apply the base preconditioner ! Apply the base preconditioner
! !
if (info == 0) call mld_baseprec_aply(sone,precv(ilev)%prec,& if (info == 0) call mld_baseprec_aply(sone,p%precv(ilev)%prec,&
& mlprec_wrk(ilev)%tx,sone,mlprec_wrk(ilev)%y2l,precv(ilev)%base_desc,& & mlprec_wrk(ilev)%tx,sone,mlprec_wrk(ilev)%y2l,p%precv(ilev)%base_desc,&
& trans,work,info) & trans,work,info)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err=' spmm/baseprec_aply') call psb_errpush(4001,name,a_err=' spmm/baseprec_aply')
@ -1015,7 +1015,7 @@ contains
! !
! Compute the output vector Y ! Compute the output vector Y
! !
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,precv(1)%base_desc,info) call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,p%precv(1)%base_desc,info)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err=' Final update') call psb_errpush(4001,name,a_err=' Final update')
@ -1052,7 +1052,7 @@ contains
! where ! where
! - M is a symmetrized hybrid multilevel domain decomposition (Schwarz) ! - M is a symmetrized hybrid multilevel domain decomposition (Schwarz)
! preconditioner associated to a certain matrix A and stored in the array ! preconditioner associated to a certain matrix A and stored in the array
! precv, ! p%precv,
! - op(M^(-1)) is M^(-1) or its transpose, according to the value of trans, ! - op(M^(-1)) is M^(-1) or its transpose, according to the value of trans,
! - X and Y are vectors, ! - X and Y are vectors,
! - alpha and beta are scalars. ! - alpha and beta are scalars.
@ -1067,7 +1067,7 @@ contains
! !
! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners', ! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners',
! each representing the part of the preconditioner associated to a certain level. ! each representing the part of the preconditioner associated to a certain level.
! For each level ilev, the base preconditioner K(ilev) is stored in precv(ilev) ! For each level ilev, the base preconditioner K(ilev) is stored in p%precv(ilev)
! and is associated to a matrix A(ilev), obtained by 'tranferring' the original ! and is associated to a matrix A(ilev), obtained by 'tranferring' the original
! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed ! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed
! aggregation. ! aggregation.
@ -1125,13 +1125,13 @@ contains
! !
! 6. Yext = beta*Yext + alpha*Y(1) ! 6. Yext = beta*Yext + alpha*Y(1)
! !
subroutine mlt_twoside_ml_aply(alpha,precv,x,beta,y,desc_data,trans,work,info) subroutine mlt_twoside_ml_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
implicit none implicit none
! Arguments ! Arguments
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
type(mld_sonelev_type), intent(in) :: precv(:) type(mld_sprec_type), intent(in) :: p
real(psb_spk_),intent(in) :: alpha,beta real(psb_spk_),intent(in) :: alpha,beta
real(psb_spk_),intent(in) :: x(:) real(psb_spk_),intent(in) :: x(:)
real(psb_spk_),intent(inout) :: y(:) real(psb_spk_),intent(inout) :: y(:)
@ -1161,9 +1161,9 @@ contains
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(precv) & ' Entry ', size(p%precv)
nlev = size(precv) nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info) allocate(mlprec_wrk(nlev),stat=info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
@ -1174,7 +1174,7 @@ contains
! !
! Copy the input vector X ! Copy the input vector X
! !
nc2l = psb_cd_get_local_cols(precv(1)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(1)%base_desc)
allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), & allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), &
& mlprec_wrk(1)%ty(nc2l), mlprec_wrk(1)%tx(nc2l), stat=info) & mlprec_wrk(1)%ty(nc2l), mlprec_wrk(1)%tx(nc2l), stat=info)
@ -1187,17 +1187,17 @@ contains
end if end if
call psb_geaxpby(sone,x,szero,mlprec_wrk(1)%x2l,& call psb_geaxpby(sone,x,szero,mlprec_wrk(1)%x2l,&
& precv(1)%base_desc,info) & p%precv(1)%base_desc,info)
call psb_geaxpby(sone,x,szero,mlprec_wrk(1)%tx,& call psb_geaxpby(sone,x,szero,mlprec_wrk(1)%tx,&
& precv(1)%base_desc,info) & p%precv(1)%base_desc,info)
! !
! STEP 2 ! STEP 2
! !
! Apply the base preconditioner at the finest level ! Apply the base preconditioner at the finest level
! !
call mld_baseprec_aply(sone,precv(1)%prec,mlprec_wrk(1)%x2l,& call mld_baseprec_aply(sone,p%precv(1)%prec,mlprec_wrk(1)%x2l,&
& szero,mlprec_wrk(1)%y2l,precv(1)%base_desc,& & szero,mlprec_wrk(1)%y2l,p%precv(1)%base_desc,&
& trans,work,info) & trans,work,info)
! !
! STEP 3 ! STEP 3
@ -1205,8 +1205,8 @@ contains
! Compute the residual at the finest level ! Compute the residual at the finest level
! !
mlprec_wrk(1)%ty = mlprec_wrk(1)%x2l mlprec_wrk(1)%ty = mlprec_wrk(1)%x2l
if (info == 0) call psb_spmm(-sone,precv(1)%base_a,mlprec_wrk(1)%y2l,& if (info == 0) call psb_spmm(-sone,p%precv(1)%base_a,mlprec_wrk(1)%y2l,&
& sone,mlprec_wrk(1)%ty,precv(1)%base_desc,info,& & sone,mlprec_wrk(1)%ty,p%precv(1)%base_desc,info,&
& work=work,trans=trans) & work=work,trans=trans)
if (info /=0) then if (info /=0) then
call psb_errpush(4010,name,a_err='Fine level baseprec/residual') call psb_errpush(4010,name,a_err='Fine level baseprec/residual')
@ -1220,8 +1220,8 @@ contains
! !
do ilev = 2, nlev do ilev = 2, nlev
nc2l = psb_cd_get_local_cols(precv(ilev)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(ilev)%base_desc)
nr2l = psb_cd_get_local_rows(precv(ilev)%base_desc) nr2l = psb_cd_get_local_rows(p%precv(ilev)%base_desc)
allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%ty(nc2l),& allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%ty(nc2l),&
& mlprec_wrk(ilev)%y2l(nc2l),mlprec_wrk(ilev)%x2l(nc2l), stat=info) & mlprec_wrk(ilev)%y2l(nc2l),mlprec_wrk(ilev)%x2l(nc2l), stat=info)
@ -1236,7 +1236,7 @@ contains
! Apply prolongator transpose, i.e. restriction ! Apply prolongator transpose, i.e. restriction
call psb_map_X2Y(sone,mlprec_wrk(ilev-1)%ty,& call psb_map_X2Y(sone,mlprec_wrk(ilev-1)%ty,&
& szero,mlprec_wrk(ilev)%x2l,& & szero,mlprec_wrk(ilev)%x2l,&
& precv(ilev)%map,info,work=work) & p%precv(ilev)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction') call psb_errpush(4001,name,a_err='Error during restriction')
@ -1244,21 +1244,21 @@ contains
end if end if
call psb_geaxpby(sone,mlprec_wrk(ilev)%x2l,szero,mlprec_wrk(ilev)%tx,& call psb_geaxpby(sone,mlprec_wrk(ilev)%x2l,szero,mlprec_wrk(ilev)%tx,&
& precv(ilev)%base_desc,info) & p%precv(ilev)%base_desc,info)
! !
! Apply the base preconditioner ! Apply the base preconditioner
! !
if (info == 0) call mld_baseprec_aply(sone,precv(ilev)%prec,& if (info == 0) call mld_baseprec_aply(sone,p%precv(ilev)%prec,&
& mlprec_wrk(ilev)%x2l,szero,mlprec_wrk(ilev)%y2l,& & mlprec_wrk(ilev)%x2l,szero,mlprec_wrk(ilev)%y2l,&
& precv(ilev)%base_desc,trans,work,info) & p%precv(ilev)%base_desc,trans,work,info)
! !
! Compute the residual (at all levels but the coarsest one) ! Compute the residual (at all levels but the coarsest one)
! !
if(ilev < nlev) then if(ilev < nlev) then
mlprec_wrk(ilev)%ty = mlprec_wrk(ilev)%x2l mlprec_wrk(ilev)%ty = mlprec_wrk(ilev)%x2l
if (info == 0) call psb_spmm(-sone,precv(ilev)%base_a,& if (info == 0) call psb_spmm(-sone,p%precv(ilev)%base_a,&
& mlprec_wrk(ilev)%y2l,sone,mlprec_wrk(ilev)%ty,& & mlprec_wrk(ilev)%y2l,sone,mlprec_wrk(ilev)%ty,&
& precv(ilev)%base_desc,info,work=work,trans=trans) & p%precv(ilev)%base_desc,info,work=work,trans=trans)
endif endif
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='baseprec_aply/residual') call psb_errpush(4001,name,a_err='baseprec_aply/residual')
@ -1279,7 +1279,7 @@ contains
! !
call psb_map_Y2X(sone,mlprec_wrk(ilev+1)%y2l,& call psb_map_Y2X(sone,mlprec_wrk(ilev+1)%y2l,&
& sone,mlprec_wrk(ilev)%y2l,& & sone,mlprec_wrk(ilev)%y2l,&
& precv(ilev+1)%map,info,work=work) & p%precv(ilev+1)%map,info,work=work)
if (info /=0 ) then if (info /=0 ) then
call psb_errpush(4001,name,a_err='Error during restriction') call psb_errpush(4001,name,a_err='Error during restriction')
@ -1289,14 +1289,14 @@ contains
! !
! Compute the residual ! Compute the residual
! !
call psb_spmm(-sone,precv(ilev)%base_a,mlprec_wrk(ilev)%y2l,& call psb_spmm(-sone,p%precv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
& sone,mlprec_wrk(ilev)%tx,precv(ilev)%base_desc,info,& & sone,mlprec_wrk(ilev)%tx,p%precv(ilev)%base_desc,info,&
& work=work,trans=trans) & work=work,trans=trans)
! !
! Apply the base preconditioner ! Apply the base preconditioner
! !
if (info == 0) call mld_baseprec_aply(sone,precv(ilev)%prec,mlprec_wrk(ilev)%tx,& if (info == 0) call mld_baseprec_aply(sone,p%precv(ilev)%prec,mlprec_wrk(ilev)%tx,&
& sone,mlprec_wrk(ilev)%y2l,precv(ilev)%base_desc, trans, work,info) & sone,mlprec_wrk(ilev)%y2l,p%precv(ilev)%base_desc, trans, work,info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4001,name,a_err='Error: residual/baseprec_aply') call psb_errpush(4001,name,a_err='Error: residual/baseprec_aply')
goto 9999 goto 9999
@ -1309,7 +1309,7 @@ contains
! Compute the output vector Y ! Compute the output vector Y
! !
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,& call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,&
& precv(1)%base_desc,info) & p%precv(1)%base_desc,info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4001,name,a_err='Error final update') call psb_errpush(4001,name,a_err='Error final update')

2
mlprec/mld_sprecaply.f90
Unescape Escape View File

@ -127,7 +127,7 @@ subroutine mld_sprecaply(prec,x,y,desc_data,info,trans,work)
goto 9999 goto 9999
end if end if
if (size(prec%precv) >1) then if (size(prec%precv) >1) then
call mld_mlprec_aply(sone,prec%precv,x,szero,y,desc_data,trans_,work_,info) call mld_mlprec_aply(sone,prec,x,szero,y,desc_data,trans_,work_,info)
if(info /= 0) then if(info /= 0) then
call psb_errpush(4010,name,a_err='mld_smlprec_aply') call psb_errpush(4010,name,a_err='mld_smlprec_aply')
goto 9999 goto 9999

224
mlprec/mld_zmlprec_aply.f90
Unescape Escape View File

@ -78,41 +78,41 @@
! Arguments: ! Arguments:
! alpha - complex(psb_dpk_), input. ! alpha - complex(psb_dpk_), input.
! The scalar alpha. ! The scalar alpha.
! precv - type(mld_zonelev_type), dimension(:), input. ! p - type(mld_zprec_type), input.
! The array of one-level preconditioner data structures containing the ! The array of one-level preconditioner data structures containing the
! local parts of the preconditioners to be applied at each level. ! local parts of the preconditioners to be applied at each level.
! Note that nlev = size(precv) = number of levels. ! Note that nlev = size(p%precv) = number of levels.
! precv(ilev)%prec - type(psb_zbaseprec_type) ! p%p%precv(ilev)%prec - type(psb_zbaseprec_type)
! The "base" preconditioner for the current level ! The "base" preconditioner for the current level
! precv(ilev)%ac - type(psb_zspmat_type) ! p%precv(ilev)%ac - type(psb_zspmat_type)
! The local part of the matrix A(ilev). ! The local part of the matrix A(ilev).
! precv(ilev)%desc_ac - type(psb_desc_type). ! p%precv(ilev)%desc_ac - type(psb_desc_type).
! The communication descriptor associated to the sparse ! The communication descriptor associated to the sparse
! matrix A(ilev) ! matrix A(ilev)
! precv(ilev)%map - type(psb_inter_desc_type) ! p%precv(ilev)%map - type(psb_inter_desc_type)
! Stores the linear operators mapping between levels ! Stores the linear operators mapping between levels
! (ilev-1) and (ilev). These are the restriction and ! (ilev-1) and (ilev). These are the restriction and
! prolongation operators described in the sequel. ! prolongation operators described in the sequel.
! precv(ilev)%iprcparm - integer, dimension(:), allocatable. ! p%precv(ilev)%iprcparm - integer, dimension(:), allocatable.
! The integer parameters defining the multilevel ! The integer parameters defining the multilevel
! strategy ! strategy
! precv(ilev)%rprcparm - real(psb_dpk_), dimension(:), allocatable. ! p%precv(ilev)%rprcparm - real(psb_dpk_), dimension(:), allocatable.
! The real parameters defining the multilevel strategy ! The real parameters defining the multilevel strategy
! precv(ilev)%mlia - integer, dimension(:), allocatable. ! p%precv(ilev)%mlia - integer, dimension(:), allocatable.
! The aggregation map (ilev-1) --> (ilev). ! The aggregation map (ilev-1) --> (ilev).
! In case of non-smoothed aggregation, it is used ! In case of non-smoothed aggregation, it is used
! instead of mld_sm_pr_. ! instead of mld_sm_pr_.
! precv(ilev)%nlaggr - integer, dimension(:), allocatable. ! p%precv(ilev)%nlaggr - integer, dimension(:), allocatable.
! The number of aggregates (rows of A(ilev)) on the ! The number of aggregates (rows of A(ilev)) on the
! various processes. ! various processes.
! precv(ilev)%base_a - type(psb_zspmat_type), pointer. ! p%precv(ilev)%base_a - type(psb_zspmat_type), pointer.
! Pointer (really a pointer!) to the base matrix of ! Pointer (really a pointer!) to the base matrix of
! the current level, i.e. the local part of A(ilev); ! the current level, i.e. the local part of A(ilev);
! so we have a unified treatment of residuals. We ! so we have a unified treatment of residuals. We
! need this to avoid passing explicitly the matrix ! need this to avoid passing explicitly the matrix
! A(ilev) to the routine which applies the ! A(ilev) to the routine which applies the
! preconditioner. ! preconditioner.
! precv(ilev)%base_desc - type(psb_desc_type), pointer. ! p%precv(ilev)%base_desc - type(psb_desc_type), pointer.
! Pointer to the communication descriptor associated ! Pointer to the communication descriptor associated
! to the sparse matrix pointed by base_a. ! to the sparse matrix pointed by base_a.
! !
@ -136,10 +136,10 @@
! Note that when the LU factorization of the matrix A(ilev) is computed instead of ! Note that when the LU factorization of the matrix A(ilev) is computed instead of
! the ILU one, by using UMFPACK or SuperLU, the corresponding L and U factors ! the ILU one, by using UMFPACK or SuperLU, the corresponding L and U factors
! are stored in data structures provided by UMFPACK or SuperLU and pointed by ! are stored in data structures provided by UMFPACK or SuperLU and pointed by
! precv(ilev)%prec%iprcparm(mld_umf_ptr) or precv(ilev)%prec%iprcparm(mld_slu_ptr), ! p%precv(ilev)%prec%iprcparm(mld_umf_ptr) or p%precv(ilev)%prec%iprcparm(mld_slu_ptr),
! respectively. ! respectively.
! !
subroutine mld_zmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info) subroutine mld_zmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
use psb_base_mod use psb_base_mod
use mld_inner_mod, mld_protect_name => mld_zmlprec_aply use mld_inner_mod, mld_protect_name => mld_zmlprec_aply
@ -148,7 +148,7 @@ subroutine mld_zmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
! Arguments ! Arguments
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
type(mld_zonelev_type), intent(in) :: precv(:) type(mld_zprec_type), intent(in) :: p
complex(psb_dpk_),intent(in) :: alpha,beta complex(psb_dpk_),intent(in) :: alpha,beta
complex(psb_dpk_),intent(in) :: x(:) complex(psb_dpk_),intent(in) :: x(:)
complex(psb_dpk_),intent(inout) :: y(:) complex(psb_dpk_),intent(inout) :: y(:)
@ -173,11 +173,11 @@ subroutine mld_zmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(precv) & ' Entry ', size(p%precv)
trans_ = psb_toupper(trans) trans_ = psb_toupper(trans)
select case(precv(2)%iprcparm(mld_ml_type_)) select case(p%precv(2)%iprcparm(mld_ml_type_))
case(mld_no_ml_) case(mld_no_ml_)
! !
@ -191,7 +191,7 @@ subroutine mld_zmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
! Additive multilevel ! Additive multilevel
! !
call add_ml_aply(alpha,precv,x,beta,y,desc_data,trans_,work,info) call add_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info)
case(mld_mult_ml_) case(mld_mult_ml_)
! !
@ -202,15 +202,15 @@ subroutine mld_zmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
! Note that the transpose switches pre <-> post. ! Note that the transpose switches pre <-> post.
! !
select case(precv(2)%iprcparm(mld_smoother_pos_)) select case(p%precv(2)%iprcparm(mld_smoother_pos_))
case(mld_post_smooth_) case(mld_post_smooth_)
select case (trans_) select case (trans_)
case('N') case('N')
call mlt_post_ml_aply(alpha,precv,x,beta,y,desc_data,trans_,work,info) call mlt_post_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info)
case('T','C') case('T','C')
call mlt_pre_ml_aply(alpha,precv,x,beta,y,desc_data,trans_,work,info) call mlt_pre_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info)
case default case default
info = 4001 info = 4001
call psb_errpush(info,name,a_err='invalid trans') call psb_errpush(info,name,a_err='invalid trans')
@ -221,9 +221,9 @@ subroutine mld_zmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
select case (trans_) select case (trans_)
case('N') case('N')
call mlt_pre_ml_aply(alpha,precv,x,beta,y,desc_data,trans_,work,info) call mlt_pre_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info)
case('T','C') case('T','C')
call mlt_post_ml_aply(alpha,precv,x,beta,y,desc_data,trans_,work,info) call mlt_post_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info)
case default case default
info = 4001 info = 4001
call psb_errpush(info,name,a_err='invalid trans') call psb_errpush(info,name,a_err='invalid trans')
@ -232,12 +232,12 @@ subroutine mld_zmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
case(mld_twoside_smooth_) case(mld_twoside_smooth_)
call mlt_twoside_ml_aply(alpha,precv,x,beta,y,desc_data,trans_,work,info) call mlt_twoside_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info)
case default case default
info = 4013 info = 4013
call psb_errpush(info,name,a_err='invalid smooth_pos',& call psb_errpush(info,name,a_err='invalid smooth_pos',&
& i_Err=(/precv(2)%iprcparm(mld_smoother_pos_),0,0,0,0/)) & i_Err=(/p%precv(2)%iprcparm(mld_smoother_pos_),0,0,0,0/))
goto 9999 goto 9999
end select end select
@ -245,7 +245,7 @@ subroutine mld_zmlprec_aply(alpha,precv,x,beta,y,desc_data,trans,work,info)
case default case default
info = 4013 info = 4013
call psb_errpush(info,name,a_err='invalid mltype',& call psb_errpush(info,name,a_err='invalid mltype',&
& i_Err=(/precv(2)%iprcparm(mld_ml_type_),0,0,0,0/)) & i_Err=(/p%precv(2)%iprcparm(mld_ml_type_),0,0,0,0/))
goto 9999 goto 9999
end select end select
@ -272,7 +272,7 @@ contains
! Y = beta*Y + alpha*op(M^(-1))*X, ! Y = beta*Y + alpha*op(M^(-1))*X,
! where ! where
! - M is an additive multilevel domain decomposition (Schwarz) preconditioner ! - M is an additive multilevel domain decomposition (Schwarz) preconditioner
! associated to a certain matrix A and stored in the array precv, ! associated to a certain matrix A and stored in the array p%precv,
! - op(M^(-1)) is M^(-1) or its (conjugate) transpose, according to ! - op(M^(-1)) is M^(-1) or its (conjugate) transpose, according to
! the value of trans, ! the value of trans,
! - X and Y are vectors, ! - X and Y are vectors,
@ -287,7 +287,7 @@ contains
! !
! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners', ! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners',
! each representing the part of the preconditioner associated to a certain level. ! each representing the part of the preconditioner associated to a certain level.
! For each level ilev, the base preconditioner K(ilev) is stored in precv(ilev) ! For each level ilev, the base preconditioner K(ilev) is stored in p%precv(ilev)
! and is associated to a matrix A(ilev), obtained by 'tranferring' the original ! and is associated to a matrix A(ilev), obtained by 'tranferring' the original
! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed ! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed
! aggregation. ! aggregation.
@ -335,13 +335,13 @@ contains
! !
! 4. Yext = beta*Yext + alpha*Y(1) ! 4. Yext = beta*Yext + alpha*Y(1)
! !
subroutine add_ml_aply(alpha,precv,x,beta,y,desc_data,trans,work,info) subroutine add_ml_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
implicit none implicit none
! Arguments ! Arguments
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
type(mld_zonelev_type), intent(in) :: precv(:) type(mld_zprec_type), intent(in) :: p
complex(psb_dpk_),intent(in) :: alpha,beta complex(psb_dpk_),intent(in) :: alpha,beta
complex(psb_dpk_),intent(in) :: x(:) complex(psb_dpk_),intent(in) :: x(:)
complex(psb_dpk_),intent(inout) :: y(:) complex(psb_dpk_),intent(inout) :: y(:)
@ -371,9 +371,9 @@ contains
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(precv) & ' Entry ', size(p%precv)
nlev = size(precv) nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info) allocate(mlprec_wrk(nlev),stat=info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
@ -396,8 +396,8 @@ contains
mlprec_wrk(1)%x2l(:) = x(:) mlprec_wrk(1)%x2l(:) = x(:)
mlprec_wrk(1)%y2l(:) = zzero mlprec_wrk(1)%y2l(:) = zzero
call mld_baseprec_aply(alpha,precv(1)%prec,x,beta,y,& call mld_baseprec_aply(alpha,p%precv(1)%prec,x,beta,y,&
& precv(1)%base_desc,trans,work,info) & p%precv(1)%base_desc,trans,work,info)
if (info /=0) then if (info /=0) then
call psb_errpush(4010,name,a_err='baseprec_aply') call psb_errpush(4010,name,a_err='baseprec_aply')
goto 9999 goto 9999
@ -408,8 +408,8 @@ contains
! For each level except the finest one ... ! For each level except the finest one ...
! !
do ilev = 2, nlev do ilev = 2, nlev
nc2l = psb_cd_get_local_cols(precv(ilev)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(ilev)%base_desc)
nr2l = psb_cd_get_local_rows(precv(ilev)%base_desc) nr2l = psb_cd_get_local_rows(p%precv(ilev)%base_desc)
allocate(mlprec_wrk(ilev)%x2l(nc2l),mlprec_wrk(ilev)%y2l(nc2l),& allocate(mlprec_wrk(ilev)%x2l(nc2l),mlprec_wrk(ilev)%y2l(nc2l),&
& stat=info) & stat=info)
if (info /= 0) then if (info /= 0) then
@ -422,7 +422,7 @@ contains
! Apply prolongator transpose, i.e. restriction ! Apply prolongator transpose, i.e. restriction
call psb_map_X2Y(zone,mlprec_wrk(ilev-1)%x2l,& call psb_map_X2Y(zone,mlprec_wrk(ilev-1)%x2l,&
& zzero,mlprec_wrk(ilev)%x2l,& & zzero,mlprec_wrk(ilev)%x2l,&
& precv(ilev)%map,info,work=work) & p%precv(ilev)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction') call psb_errpush(4001,name,a_err='Error during restriction')
@ -432,9 +432,9 @@ contains
! !
! Apply the base preconditioner ! Apply the base preconditioner
! !
call mld_baseprec_aply(zone,precv(ilev)%prec,& call mld_baseprec_aply(zone,p%precv(ilev)%prec,&
& mlprec_wrk(ilev)%x2l,zzero,mlprec_wrk(ilev)%y2l,& & mlprec_wrk(ilev)%x2l,zzero,mlprec_wrk(ilev)%y2l,&
& precv(ilev)%base_desc,trans,work,info) & p%precv(ilev)%base_desc,trans,work,info)
enddo enddo
@ -445,15 +445,15 @@ contains
! !
do ilev =nlev,2,-1 do ilev =nlev,2,-1
nc2l = psb_cd_get_local_cols(precv(ilev)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(ilev)%base_desc)
nr2l = psb_cd_get_local_rows(precv(ilev)%base_desc) nr2l = psb_cd_get_local_rows(p%precv(ilev)%base_desc)
! !
! Apply prolongator ! Apply prolongator
! !
call psb_map_Y2X(zone,mlprec_wrk(ilev)%y2l,& call psb_map_Y2X(zone,mlprec_wrk(ilev)%y2l,&
& zone,mlprec_wrk(ilev-1)%y2l,& & zone,mlprec_wrk(ilev-1)%y2l,&
& precv(ilev)%map,info,work=work) & p%precv(ilev)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during prolongation') call psb_errpush(4001,name,a_err='Error during prolongation')
@ -466,7 +466,7 @@ contains
! !
! Compute the output vector Y ! Compute the output vector Y
! !
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,zone,y,precv(1)%base_desc,info) call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,zone,y,p%precv(1)%base_desc,info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4001,name,a_err='Error on final update') call psb_errpush(4001,name,a_err='Error on final update')
goto 9999 goto 9999
@ -500,7 +500,7 @@ contains
! Y = beta*Y + alpha*op(M^(-1))*X, ! Y = beta*Y + alpha*op(M^(-1))*X,
! where ! where
! - M is a hybrid multilevel domain decomposition (Schwarz) preconditioner ! - M is a hybrid multilevel domain decomposition (Schwarz) preconditioner
! associated to a certain matrix A and stored in the array precv, ! associated to a certain matrix A and stored in the array p%precv,
! - op(M^(-1)) is M^(-1) or its (conjugate) transpose, according to ! - op(M^(-1)) is M^(-1) or its (conjugate) transpose, according to
! the value of trans, ! the value of trans,
! - X and Y are vectors, ! - X and Y are vectors,
@ -515,7 +515,7 @@ contains
! !
! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners', ! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners',
! each representing the part of the preconditioner associated to a certain level. ! each representing the part of the preconditioner associated to a certain level.
! For each level ilev, the base preconditioner K(ilev) is stored in precv(ilev) ! For each level ilev, the base preconditioner K(ilev) is stored in p%precv(ilev)
! and is associated to a matrix A(ilev), obtained by 'tranferring' the original ! and is associated to a matrix A(ilev), obtained by 'tranferring' the original
! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed ! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed
! aggregation. ! aggregation.
@ -571,13 +571,13 @@ contains
! 6. Yext = beta*Yext + alpha*Y(1) ! 6. Yext = beta*Yext + alpha*Y(1)
! !
! !
subroutine mlt_pre_ml_aply(alpha,precv,x,beta,y,desc_data,trans,work,info) subroutine mlt_pre_ml_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
implicit none implicit none
! Arguments ! Arguments
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
type(mld_zonelev_type), intent(in) :: precv(:) type(mld_zprec_type), intent(in) :: p
complex(psb_dpk_),intent(in) :: alpha,beta complex(psb_dpk_),intent(in) :: alpha,beta
complex(psb_dpk_),intent(in) :: x(:) complex(psb_dpk_),intent(in) :: x(:)
complex(psb_dpk_),intent(inout) :: y(:) complex(psb_dpk_),intent(inout) :: y(:)
@ -607,9 +607,9 @@ contains
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(precv) & ' Entry ', size(p%precv)
nlev = size(precv) nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info) allocate(mlprec_wrk(nlev),stat=info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
@ -621,7 +621,7 @@ contains
! !
! Copy the input vector X ! Copy the input vector X
! !
nc2l = psb_cd_get_local_cols(precv(1)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(1)%base_desc)
allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), & allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), &
& mlprec_wrk(1)%tx(nc2l), stat=info) & mlprec_wrk(1)%tx(nc2l), stat=info)
@ -638,8 +638,8 @@ contains
! !
! Apply the base preconditioner at the finest level ! Apply the base preconditioner at the finest level
! !
call mld_baseprec_aply(zone,precv(1)%prec,mlprec_wrk(1)%x2l,& call mld_baseprec_aply(zone,p%precv(1)%prec,mlprec_wrk(1)%x2l,&
& zzero,mlprec_wrk(1)%y2l,precv(1)%base_desc,& & zzero,mlprec_wrk(1)%y2l,p%precv(1)%base_desc,&
& trans,work,info) & trans,work,info)
if (info /=0) then if (info /=0) then
call psb_errpush(4010,name,a_err=' baseprec_aply') call psb_errpush(4010,name,a_err=' baseprec_aply')
@ -653,8 +653,8 @@ contains
! !
mlprec_wrk(1)%tx = mlprec_wrk(1)%x2l mlprec_wrk(1)%tx = mlprec_wrk(1)%x2l
call psb_spmm(-zone,precv(1)%base_a,mlprec_wrk(1)%y2l,& call psb_spmm(-zone,p%precv(1)%base_a,mlprec_wrk(1)%y2l,&
& zone,mlprec_wrk(1)%tx,precv(1)%base_desc,info,& & zone,mlprec_wrk(1)%tx,p%precv(1)%base_desc,info,&
& work=work,trans=trans) & work=work,trans=trans)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err=' fine level residual') call psb_errpush(4001,name,a_err=' fine level residual')
@ -668,8 +668,8 @@ contains
! !
do ilev = 2, nlev do ilev = 2, nlev
nc2l = psb_cd_get_local_cols(precv(ilev)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(ilev)%base_desc)
nr2l = psb_cd_get_local_rows(precv(ilev)%base_desc) nr2l = psb_cd_get_local_rows(p%precv(ilev)%base_desc)
allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%y2l(nc2l),& allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%y2l(nc2l),&
& mlprec_wrk(ilev)%x2l(nc2l), stat=info) & mlprec_wrk(ilev)%x2l(nc2l), stat=info)
@ -683,7 +683,7 @@ contains
! Apply prolongator transpose, i.e. restriction ! Apply prolongator transpose, i.e. restriction
call psb_map_X2Y(zone,mlprec_wrk(ilev-1)%tx,& call psb_map_X2Y(zone,mlprec_wrk(ilev-1)%tx,&
& zzero,mlprec_wrk(ilev)%x2l,& & zzero,mlprec_wrk(ilev)%x2l,&
& precv(ilev)%map,info,work=work) & p%precv(ilev)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction') call psb_errpush(4001,name,a_err='Error during restriction')
@ -693,17 +693,17 @@ contains
! !
! Apply the base preconditioner ! Apply the base preconditioner
! !
call mld_baseprec_aply(zone,precv(ilev)%prec,mlprec_wrk(ilev)%x2l,& call mld_baseprec_aply(zone,p%precv(ilev)%prec,mlprec_wrk(ilev)%x2l,&
& zzero,mlprec_wrk(ilev)%y2l,precv(ilev)%base_desc,trans,work,info) & zzero,mlprec_wrk(ilev)%y2l,p%precv(ilev)%base_desc,trans,work,info)
! !
! Compute the residual (at all levels but the coarsest one) ! Compute the residual (at all levels but the coarsest one)
! !
if (ilev < nlev) then if (ilev < nlev) then
mlprec_wrk(ilev)%tx = mlprec_wrk(ilev)%x2l mlprec_wrk(ilev)%tx = mlprec_wrk(ilev)%x2l
if (info == 0) call psb_spmm(-zone,precv(ilev)%base_a,& if (info == 0) call psb_spmm(-zone,p%precv(ilev)%base_a,&
& mlprec_wrk(ilev)%y2l,zone,mlprec_wrk(ilev)%tx,& & mlprec_wrk(ilev)%y2l,zone,mlprec_wrk(ilev)%tx,&
& precv(ilev)%base_desc,info,work=work,trans=trans) & p%precv(ilev)%base_desc,info,work=work,trans=trans)
endif endif
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error on up sweep residual') call psb_errpush(4001,name,a_err='Error on up sweep residual')
@ -722,7 +722,7 @@ contains
! !
call psb_map_Y2X(zone,mlprec_wrk(ilev+1)%y2l,& call psb_map_Y2X(zone,mlprec_wrk(ilev+1)%y2l,&
& zone,mlprec_wrk(ilev)%y2l,& & zone,mlprec_wrk(ilev)%y2l,&
& precv(ilev+1)%map,info,work=work) & p%precv(ilev+1)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during prolongation') call psb_errpush(4001,name,a_err='Error during prolongation')
@ -736,7 +736,7 @@ contains
! Compute the output vector Y ! Compute the output vector Y
! !
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,& call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,&
& precv(1)%base_desc,info) & p%precv(1)%base_desc,info)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error on final update') call psb_errpush(4001,name,a_err='Error on final update')
goto 9999 goto 9999
@ -769,7 +769,7 @@ contains
! Y = beta*Y + alpha*op(M^(-1))*X, ! Y = beta*Y + alpha*op(M^(-1))*X,
! where ! where
! - M is a hybrid multilevel domain decomposition (Schwarz) preconditioner ! - M is a hybrid multilevel domain decomposition (Schwarz) preconditioner
! associated to a certain matrix A and stored in the array precv, ! associated to a certain matrix A and stored in the array p%precv,
! - op(M^(-1)) is M^(-1) or its (conjugate) transpose, according to ! - op(M^(-1)) is M^(-1) or its (conjugate) transpose, according to
! the value of trans, ! the value of trans,
! - X and Y are vectors, ! - X and Y are vectors,
@ -784,7 +784,7 @@ contains
! !
! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners', ! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners',
! each representing the part of the preconditioner associated to a certain level. ! each representing the part of the preconditioner associated to a certain level.
! For each level ilev, the base preconditioner K(ilev) is stored in precv(ilev) ! For each level ilev, the base preconditioner K(ilev) is stored in p%precv(ilev)
! and is associated to a matrix A(ilev), obtained by 'tranferring' the original ! and is associated to a matrix A(ilev), obtained by 'tranferring' the original
! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed ! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed
! aggregation. ! aggregation.
@ -831,13 +831,13 @@ contains
! 5. Yext = beta*Yext + alpha*Y(1) ! 5. Yext = beta*Yext + alpha*Y(1)
! !
! !
subroutine mlt_post_ml_aply(alpha,precv,x,beta,y,desc_data,trans,work,info) subroutine mlt_post_ml_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
implicit none implicit none
! Arguments ! Arguments
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
type(mld_zonelev_type), intent(in) :: precv(:) type(mld_zprec_type), intent(in) :: p
complex(psb_dpk_),intent(in) :: alpha,beta complex(psb_dpk_),intent(in) :: alpha,beta
complex(psb_dpk_),intent(in) :: x(:) complex(psb_dpk_),intent(in) :: x(:)
complex(psb_dpk_),intent(inout) :: y(:) complex(psb_dpk_),intent(inout) :: y(:)
@ -867,9 +867,9 @@ contains
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(precv) & ' Entry ', size(p%precv)
nlev = size(precv) nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info) allocate(mlprec_wrk(nlev),stat=info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
@ -885,7 +885,7 @@ contains
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' desc_data status',allocated(desc_data%matrix_data) & ' desc_data status',allocated(desc_data%matrix_data)
nc2l = psb_cd_get_local_cols(precv(1)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(1)%base_desc)
allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), & allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), &
& mlprec_wrk(1)%tx(nc2l), stat=info) & mlprec_wrk(1)%tx(nc2l), stat=info)
@ -897,9 +897,9 @@ contains
end if end if
call psb_geaxpby(zone,x,zzero,mlprec_wrk(1)%tx,& call psb_geaxpby(zone,x,zzero,mlprec_wrk(1)%tx,&
& precv(1)%base_desc,info) & p%precv(1)%base_desc,info)
call psb_geaxpby(zone,x,zzero,mlprec_wrk(1)%x2l,& call psb_geaxpby(zone,x,zzero,mlprec_wrk(1)%x2l,&
& precv(1)%base_desc,info) & p%precv(1)%base_desc,info)
! !
! STEP 2 ! STEP 2
@ -908,13 +908,13 @@ contains
! !
do ilev=2, nlev do ilev=2, nlev
nc2l = psb_cd_get_local_cols(precv(ilev)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(ilev)%base_desc)
nr2l = psb_cd_get_local_rows(precv(ilev)%base_desc) nr2l = psb_cd_get_local_rows(p%precv(ilev)%base_desc)
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name), & & write(debug_unit,*) me,' ',trim(name), &
& ' starting up sweep ',& & ' starting up sweep ',&
& ilev,allocated(precv(ilev)%iprcparm),nc2l, nr2l & ilev,allocated(p%precv(ilev)%iprcparm),nc2l, nr2l
allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%y2l(nc2l),& allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%y2l(nc2l),&
& mlprec_wrk(ilev)%x2l(nc2l), stat=info) & mlprec_wrk(ilev)%x2l(nc2l), stat=info)
@ -929,7 +929,7 @@ contains
! Apply prolongator transpose, i.e. restriction ! Apply prolongator transpose, i.e. restriction
call psb_map_X2Y(zone,mlprec_wrk(ilev-1)%x2l,& call psb_map_X2Y(zone,mlprec_wrk(ilev-1)%x2l,&
& zzero,mlprec_wrk(ilev)%x2l,& & zzero,mlprec_wrk(ilev)%x2l,&
& precv(ilev)%map,info,work=work) & p%precv(ilev)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction') call psb_errpush(4001,name,a_err='Error during restriction')
@ -940,7 +940,7 @@ contains
! update x2l ! update x2l
! !
call psb_geaxpby(zone,mlprec_wrk(ilev)%x2l,zzero,mlprec_wrk(ilev)%tx,& call psb_geaxpby(zone,mlprec_wrk(ilev)%x2l,zzero,mlprec_wrk(ilev)%tx,&
& precv(ilev)%base_desc,info) & p%precv(ilev)%base_desc,info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4001,name,a_err='Error in update') call psb_errpush(4001,name,a_err='Error in update')
goto 9999 goto 9999
@ -957,8 +957,8 @@ contains
! !
! Apply the base preconditioner at the coarsest level ! Apply the base preconditioner at the coarsest level
! !
call mld_baseprec_aply(zone,precv(nlev)%prec,mlprec_wrk(nlev)%x2l, & call mld_baseprec_aply(zone,p%precv(nlev)%prec,mlprec_wrk(nlev)%x2l, &
& zzero, mlprec_wrk(nlev)%y2l,precv(nlev)%base_desc,trans,work,info) & zzero, mlprec_wrk(nlev)%y2l,p%precv(nlev)%base_desc,trans,work,info)
if (info /=0) then if (info /=0) then
call psb_errpush(4010,name,a_err='baseprec_aply') call psb_errpush(4010,name,a_err='baseprec_aply')
goto 9999 goto 9999
@ -983,7 +983,7 @@ contains
! !
call psb_map_Y2X(zone,mlprec_wrk(ilev+1)%y2l,& call psb_map_Y2X(zone,mlprec_wrk(ilev+1)%y2l,&
& zzero,mlprec_wrk(ilev)%y2l,& & zzero,mlprec_wrk(ilev)%y2l,&
& precv(ilev+1)%map,info,work=work) & p%precv(ilev+1)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during prolongation') call psb_errpush(4001,name,a_err='Error during prolongation')
@ -993,15 +993,15 @@ contains
! !
! Compute the residual ! Compute the residual
! !
call psb_spmm(-zone,precv(ilev)%base_a,mlprec_wrk(ilev)%y2l,& call psb_spmm(-zone,p%precv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
& zone,mlprec_wrk(ilev)%tx,precv(ilev)%base_desc,info,& & zone,mlprec_wrk(ilev)%tx,p%precv(ilev)%base_desc,info,&
& work=work,trans=trans) & work=work,trans=trans)
! !
! Apply the base preconditioner ! Apply the base preconditioner
! !
if (info == 0) call mld_baseprec_aply(zone,precv(ilev)%prec,& if (info == 0) call mld_baseprec_aply(zone,p%precv(ilev)%prec,&
& mlprec_wrk(ilev)%tx,zone,mlprec_wrk(ilev)%y2l,precv(ilev)%base_desc,& & mlprec_wrk(ilev)%tx,zone,mlprec_wrk(ilev)%y2l,p%precv(ilev)%base_desc,&
&trans,work,info) &trans,work,info)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err=' spmm/baseprec_aply') call psb_errpush(4001,name,a_err=' spmm/baseprec_aply')
@ -1018,7 +1018,7 @@ contains
! !
! Compute the output vector Y ! Compute the output vector Y
! !
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,precv(1)%base_desc,info) call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,p%precv(1)%base_desc,info)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err=' Final update') call psb_errpush(4001,name,a_err=' Final update')
@ -1055,7 +1055,7 @@ contains
! where ! where
! - M is a symmetrized hybrid multilevel domain decomposition (Schwarz) ! - M is a symmetrized hybrid multilevel domain decomposition (Schwarz)
! preconditioner associated to a certain matrix A and stored in the array ! preconditioner associated to a certain matrix A and stored in the array
! precv, ! p%precv,
! - op(M^(-1)) is M^(-1) or its (conjugate) transpose, according to ! - op(M^(-1)) is M^(-1) or its (conjugate) transpose, according to
! the value of trans, ! the value of trans,
! - X and Y are vectors, ! - X and Y are vectors,
@ -1071,7 +1071,7 @@ contains
! !
! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners', ! The multilevel preconditioner M is regarded as an array of 'one-level preconditioners',
! each representing the part of the preconditioner associated to a certain level. ! each representing the part of the preconditioner associated to a certain level.
! For each level ilev, the base preconditioner K(ilev) is stored in precv(ilev) ! For each level ilev, the base preconditioner K(ilev) is stored in p%precv(ilev)
! and is associated to a matrix A(ilev), obtained by 'tranferring' the original ! and is associated to a matrix A(ilev), obtained by 'tranferring' the original
! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed ! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed
! aggregation. ! aggregation.
@ -1129,13 +1129,13 @@ contains
! !
! 6. Yext = beta*Yext + alpha*Y(1) ! 6. Yext = beta*Yext + alpha*Y(1)
! !
subroutine mlt_twoside_ml_aply(alpha,precv,x,beta,y,desc_data,trans,work,info) subroutine mlt_twoside_ml_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
implicit none implicit none
! Arguments ! Arguments
type(psb_desc_type),intent(in) :: desc_data type(psb_desc_type),intent(in) :: desc_data
type(mld_zonelev_type), intent(in) :: precv(:) type(mld_zprec_type), intent(in) :: p
complex(psb_dpk_),intent(in) :: alpha,beta complex(psb_dpk_),intent(in) :: alpha,beta
complex(psb_dpk_),intent(in) :: x(:) complex(psb_dpk_),intent(in) :: x(:)
complex(psb_dpk_),intent(inout) :: y(:) complex(psb_dpk_),intent(inout) :: y(:)
@ -1165,9 +1165,9 @@ contains
if (debug_level >= psb_debug_inner_) & if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(precv) & ' Entry ', size(p%precv)
nlev = size(precv) nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info) allocate(mlprec_wrk(nlev),stat=info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
@ -1178,7 +1178,7 @@ contains
! !
! Copy the input vector X ! Copy the input vector X
! !
nc2l = psb_cd_get_local_cols(precv(1)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(1)%base_desc)
allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), & allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), &
& mlprec_wrk(1)%ty(nc2l), mlprec_wrk(1)%tx(nc2l), stat=info) & mlprec_wrk(1)%ty(nc2l), mlprec_wrk(1)%tx(nc2l), stat=info)
@ -1191,17 +1191,17 @@ contains
end if end if
call psb_geaxpby(zone,x,zzero,mlprec_wrk(1)%x2l,& call psb_geaxpby(zone,x,zzero,mlprec_wrk(1)%x2l,&
& precv(1)%base_desc,info) & p%precv(1)%base_desc,info)
call psb_geaxpby(zone,x,zzero,mlprec_wrk(1)%tx,& call psb_geaxpby(zone,x,zzero,mlprec_wrk(1)%tx,&
& precv(1)%base_desc,info) & p%precv(1)%base_desc,info)
! !
! STEP 2 ! STEP 2
! !
! Apply the base preconditioner at the finest level ! Apply the base preconditioner at the finest level
! !
call mld_baseprec_aply(zone,precv(1)%prec,mlprec_wrk(1)%x2l,& call mld_baseprec_aply(zone,p%precv(1)%prec,mlprec_wrk(1)%x2l,&
& zzero,mlprec_wrk(1)%y2l,precv(1)%base_desc,& & zzero,mlprec_wrk(1)%y2l,p%precv(1)%base_desc,&
& trans,work,info) & trans,work,info)
! !
! STEP 3 ! STEP 3
@ -1209,8 +1209,8 @@ contains
! Compute the residual at the finest level ! Compute the residual at the finest level
! !
mlprec_wrk(1)%ty = mlprec_wrk(1)%x2l mlprec_wrk(1)%ty = mlprec_wrk(1)%x2l
if (info == 0) call psb_spmm(-zone,precv(1)%base_a,mlprec_wrk(1)%y2l,& if (info == 0) call psb_spmm(-zone,p%precv(1)%base_a,mlprec_wrk(1)%y2l,&
& zone,mlprec_wrk(1)%ty,precv(1)%base_desc,info,& & zone,mlprec_wrk(1)%ty,p%precv(1)%base_desc,info,&
& work=work,trans=trans) & work=work,trans=trans)
if (info /=0) then if (info /=0) then
call psb_errpush(4010,name,a_err='Fine level baseprec/residual') call psb_errpush(4010,name,a_err='Fine level baseprec/residual')
@ -1224,8 +1224,8 @@ contains
! !
do ilev = 2, nlev do ilev = 2, nlev
nc2l = psb_cd_get_local_cols(precv(ilev)%base_desc) nc2l = psb_cd_get_local_cols(p%precv(ilev)%base_desc)
nr2l = psb_cd_get_local_rows(precv(ilev)%base_desc) nr2l = psb_cd_get_local_rows(p%precv(ilev)%base_desc)
allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%ty(nc2l),& allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%ty(nc2l),&
& mlprec_wrk(ilev)%y2l(nc2l),mlprec_wrk(ilev)%x2l(nc2l), stat=info) & mlprec_wrk(ilev)%y2l(nc2l),mlprec_wrk(ilev)%x2l(nc2l), stat=info)
@ -1240,7 +1240,7 @@ contains
! Apply prolongator transpose, i.e. restriction ! Apply prolongator transpose, i.e. restriction
call psb_map_X2Y(zone,mlprec_wrk(ilev-1)%ty,& call psb_map_X2Y(zone,mlprec_wrk(ilev-1)%ty,&
& zzero,mlprec_wrk(ilev)%x2l,& & zzero,mlprec_wrk(ilev)%x2l,&
& precv(ilev)%map,info,work=work) & p%precv(ilev)%map,info,work=work)
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction') call psb_errpush(4001,name,a_err='Error during restriction')
@ -1248,21 +1248,21 @@ contains
end if end if
call psb_geaxpby(zone,mlprec_wrk(ilev)%x2l,zzero,mlprec_wrk(ilev)%tx,& call psb_geaxpby(zone,mlprec_wrk(ilev)%x2l,zzero,mlprec_wrk(ilev)%tx,&
& precv(ilev)%base_desc,info) & p%precv(ilev)%base_desc,info)
! !
! Apply the base preconditioner ! Apply the base preconditioner
! !
if (info == 0) call mld_baseprec_aply(zone,precv(ilev)%prec,& if (info == 0) call mld_baseprec_aply(zone,p%precv(ilev)%prec,&
& mlprec_wrk(ilev)%x2l,zzero,mlprec_wrk(ilev)%y2l,& & mlprec_wrk(ilev)%x2l,zzero,mlprec_wrk(ilev)%y2l,&
&precv(ilev)%base_desc,trans,work,info) &p%precv(ilev)%base_desc,trans,work,info)
! !
! Compute the residual (at all levels but the coarsest one) ! Compute the residual (at all levels but the coarsest one)
! !
if(ilev < nlev) then if(ilev < nlev) then
mlprec_wrk(ilev)%ty = mlprec_wrk(ilev)%x2l mlprec_wrk(ilev)%ty = mlprec_wrk(ilev)%x2l
if (info == 0) call psb_spmm(-zone,precv(ilev)%base_a,& if (info == 0) call psb_spmm(-zone,p%precv(ilev)%base_a,&
& mlprec_wrk(ilev)%y2l,zone,mlprec_wrk(ilev)%ty,& & mlprec_wrk(ilev)%y2l,zone,mlprec_wrk(ilev)%ty,&
& precv(ilev)%base_desc,info,work=work,trans=trans) & p%precv(ilev)%base_desc,info,work=work,trans=trans)
endif endif
if (info /=0) then if (info /=0) then
call psb_errpush(4001,name,a_err='baseprec_aply/residual') call psb_errpush(4001,name,a_err='baseprec_aply/residual')
@ -1283,7 +1283,7 @@ contains
! !
call psb_map_Y2X(zone,mlprec_wrk(ilev+1)%y2l,& call psb_map_Y2X(zone,mlprec_wrk(ilev+1)%y2l,&
& zone,mlprec_wrk(ilev)%y2l,& & zone,mlprec_wrk(ilev)%y2l,&
& precv(ilev+1)%map,info,work=work) & p%precv(ilev+1)%map,info,work=work)
if (info /=0 ) then if (info /=0 ) then
call psb_errpush(4001,name,a_err='Error during restriction') call psb_errpush(4001,name,a_err='Error during restriction')
@ -1293,14 +1293,14 @@ contains
! !
! Compute the residual ! Compute the residual
! !
call psb_spmm(-zone,precv(ilev)%base_a,mlprec_wrk(ilev)%y2l,& call psb_spmm(-zone,p%precv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
& zone,mlprec_wrk(ilev)%tx,precv(ilev)%base_desc,info,& & zone,mlprec_wrk(ilev)%tx,p%precv(ilev)%base_desc,info,&
& work=work,trans=trans) & work=work,trans=trans)
! !
! Apply the base preconditioner ! Apply the base preconditioner
! !
if (info == 0) call mld_baseprec_aply(zone,precv(ilev)%prec,mlprec_wrk(ilev)%tx,& if (info == 0) call mld_baseprec_aply(zone,p%precv(ilev)%prec,mlprec_wrk(ilev)%tx,&
& zone,mlprec_wrk(ilev)%y2l,precv(ilev)%base_desc, trans, work,info) & zone,mlprec_wrk(ilev)%y2l,p%precv(ilev)%base_desc, trans, work,info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4001,name,a_err='Error: residual/baseprec_aply') call psb_errpush(4001,name,a_err='Error: residual/baseprec_aply')
goto 9999 goto 9999
@ -1313,7 +1313,7 @@ contains
! Compute the output vector Y ! Compute the output vector Y
! !
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,& call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,&
& precv(1)%base_desc,info) & p%precv(1)%base_desc,info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4001,name,a_err='Error final update') call psb_errpush(4001,name,a_err='Error final update')

2
mlprec/mld_zprecaply.f90
Unescape Escape View File

@ -127,7 +127,7 @@ subroutine mld_zprecaply(prec,x,y,desc_data,info,trans,work)
goto 9999 goto 9999
end if end if
if (size(prec%precv) >1) then if (size(prec%precv) >1) then
call mld_mlprec_aply(zone,prec%precv,x,zzero,y,desc_data,trans_,work_,info) call mld_mlprec_aply(zone,prec,x,zzero,y,desc_data,trans_,work_,info)
if(info /= 0) then if(info /= 0) then
call psb_errpush(4010,name,a_err='mld_zmlprec_aply') call psb_errpush(4010,name,a_err='mld_zmlprec_aply')
goto 9999 goto 9999

Write Preview
Loading…
Cancel
Save