diff --git a/mlprec/mld_dmlprec_aply.f90 b/mlprec/mld_dmlprec_aply.f90 index 2d75dec4..75a742a7 100644 --- a/mlprec/mld_dmlprec_aply.f90 +++ b/mlprec/mld_dmlprec_aply.f90 @@ -138,7 +138,180 @@ ! are stored in data structures provided by UMFPACK or SuperLU and pointed by ! p%precv(ilev)%prec%iprcparm(mld_umf_ptr) or p%precv(ilev)%prec%iprcparm(mld_slu_ptr), ! respectively. +! +! This routine is formulated in a recursive way, so it is very compact. +! In the original code the recursive formulation was explicitly unrolled. +! The description of the various alternatives is given below in the explicit +! formulation, hopefully it will be clear enough when related to the +! recursive formulation. +! +! This routine computes +! Y = beta*Y + alpha*op(M^(-1))*X, +! where +! - M is a multilevel domain decomposition (Schwarz) preconditioner +! associated to a certain matrix A and stored in p, +! - op(M^(-1)) is M^(-1) or its transpose, according to the value of trans, +! - X and Y are vectors, +! - alpha and beta are scalars. +! +! For each level we have as many submatrices as processes (except for the coarsest +! level where we might have a replicated index space) and each process takes care +! of one submatrix. +! +! The multilevel preconditioner is regarded as an array of 'one-level' data structures, +! each containing the part of the preconditioner associated to a certain level +! (for more details see the description of mld_Tonelev_type in mld_prec_type.f90). +! For each level ilev, the 'base preconditioner' K(ilev) is stored in +! p%precv(ilev)%prec +! 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 +! aggregation. +! The levels are numbered in increasing order starting from the finest one, i.e. +! level 1 is the finest level and A(1) is the matrix A. +! +! +! Additive multilevel +! This is additive both within the levels and among levels. +! +! For details on the additive multilevel Schwarz preconditioner see the +! Algorithm 3.1.1 in the book: +! B.F. Smith, P.E. Bjorstad & W.D. Gropp, +! Domain decomposition: parallel multilevel methods for elliptic partial +! differential equations, Cambridge University Press, 1996. +! +! (P(ilev) denotes the smoothed prolongator from level ilev to level +! ilev-1, while PT(ilev) denotes its transpose, i.e. the corresponding +! restriction operator from level ilev-1 to level ilev). +! +! 0. Transfer the outer vector Xest to X(1) (inner X at level 1). +! +! 1. If ilev >1 Transfer X(ilev-1) to the current level. +! X(ilev) = PT(ilev)*X(ilev-1) +! +! 2. Apply the base preconditioner at the current level. +! ! The sum over the subdomains is carried out in the +! ! application of K(ilev). +! Y(ilev) = (K(ilev)^(-1))*X(ilev) +! +! 3. If ilev < nlevel +! a. Call recursively itself +! b. Transfer Y(ilev+1) to the current level. +! Y(ilev) = Y(ilev) + P(ilev+1)*Y(ilev+1) +! +! 4. if ilev == 1 Transfer the inner Y to the external +! Yext = beta*Yext + alpha*Y(1) +! +! +! +! Hybrid multiplicative---pre-smoothing ! +! The preconditioner M is hybrid in the sense that it is multiplicative through the +! levels and additive inside a level. +! +! For details on the pre-smoothed hybrid multiplicative multilevel Schwarz +! preconditioner, see the Algorithm 3.2.1 in the book: +! B.F. Smith, P.E. Bjorstad & W.D. Gropp, +! Domain decomposition: parallel multilevel methods for elliptic partial +! differential equations, Cambridge University Press, 1996. +! +! +! 0. Transfer the outer vector Xest to X(1) (inner X at level 1). +! +! 1. If ilev >1 Transfer X(ilev-1) to the current level. +! X(ilev) = PT(ilev)*X(ilev-1) +! +! 2. Apply the base preconditioner at the current level. +! ! The sum over the subdomains is carried out in the +! ! application of K(ilev). +! Y(ilev) = (K(ilev)^(-1))*X(ilev) +! +! 3. If ilev < nlevel +! a. Compute the residula +! r(ilev) = y(ilev) - A(ilev)*x(ilev) +! b. Call recursively itself passing +! r(ilev) for transfer to the next level +! +! c. Transfer Y(ilev+1) to the current level. +! Y(ilev) = Y(ilev) + P(ilev+1)*Y(ilev+1) +! +! 4. if ilev == 1 Transfer the inner Y to the external +! Yext = beta*Yext + alpha*Y(1) +! +! +! +! Hybrid multiplicative, post-smoothing variant +! +! +! +! 0. Transfer the outer vector Xest to X(1) (inner X at level 1). +! +! 1. If ilev >1 Transfer X(ilev-1) to the current level. +! X(ilev) = PT(ilev)*X(ilev-1) +! +! 2. If ilev < nlev +! +! a. Call recursively itself passing +! r(ilev) for transfer to the next level +! b. Transfer Y(ilev+1) to the current level. +! Y(ilev) = Y(ilev) + P(ilev+1)*Y(ilev+1) +! +! c. Compute the residual +! r(ilev) = y(ilev) - A(ilev)*x(ilev) +! +! 3. Apply the base preconditioner at the current level to the residual +! ! The sum over the subdomains is carried out in the +! ! application of K(ilev). +! Y(ilev) = y(ilev) + (K(ilev)^(-1))*r(ilev) +! +! +! 4. if ilev == 1 Transfer the inner Y to the external +! Yext = beta*Yext + alpha*Y(1) +! +! +! +! Hybrid multiplicative, pre- and post-smoothing (two-side) variant +! +! +! For details on the symmetrized hybrid multiplicative multilevel Schwarz +! preconditioner, see the Algorithm 3.2.2 of the book: +! B.F. Smith, P.E. Bjorstad & W.D. Gropp, +! Domain decomposition: parallel multilevel methods for elliptic partial +! differential equations, Cambridge University Press, 1996. +! +! +! 0. Transfer the outer vector Xest to X(1) (inner X at level 1). +! +! 1. If ilev >1 Transfer X(ilev-1) to the current level. +! X(ilev) = PT(ilev)*X(ilev-1) +! +! 2. Apply the base preconditioner at the current level. +! ! The sum over the subdomains is carried out in the +! ! application of K(ilev). +! Y(ilev) = (K(ilev)^(-1))*X(ilev) +! +! 3. If ilev < nlevel +! a. Compute the residula +! r(ilev) = y(ilev) - A(ilev)*x(ilev) +! b. Call recursively itself passing +! r(ilev) for transfer to the next level +! +! c. Transfer Y(ilev+1) to the current level. +! Y(ilev) = Y(ilev) + P(ilev+1)*Y(ilev+1) +! +! d. Compute the residual +! r(ilev) = y(ilev) - A(ilev)*x(ilev) +! +! 4. Apply the base preconditioner at the current level to the residual +! ! The sum over the subdomains is carried out in the +! ! application of K(ilev). +! Y(ilev) = y(ilev) + (K(ilev)^(-1))*r(ilev) +! +! +! 5. if ilev == 1 Transfer the inner Y to the external +! Yext = beta*Yext + alpha*Y(1) +! +! +! subroutine mld_dmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info) use psb_base_mod @@ -181,120 +354,44 @@ subroutine mld_dmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info) trans_ = psb_toupper(trans) - if (.false.) then - nlev = size(p%precv) - allocate(mlprec_wrk(nlev),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - level = 1 - - nc2l = psb_cd_get_local_cols(p%precv(level)%base_desc) - nr2l = psb_cd_get_local_rows(p%precv(level)%base_desc) - allocate(mlprec_wrk(level)%x2l(nc2l),mlprec_wrk(level)%y2l(nc2l),& - & stat=info) - if (info /= 0) then - info=4025 - call psb_errpush(info,name,i_err=(/size(x)+size(y),0,0,0,0/),& - & a_err='real(psb_dpk_)') - goto 9999 - end if - - mlprec_wrk(level)%x2l(:) = x(:) - mlprec_wrk(level)%y2l(:) = dzero - - call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info) - - if (info /= 0) then - call psb_errpush(4001,name,a_err='Inner prec aply') - goto 9999 - end if - - call psb_geaxpby(alpha,mlprec_wrk(level)%y2l,beta,y,& - & p%precv(level)%base_desc,info) - - if (info /= 0) then - call psb_errpush(4001,name,a_err='Error final update') - goto 9999 - end if - - - else - select case(p%precv(2)%iprcparm(mld_ml_type_)) - - case(mld_no_ml_) - ! - ! No preconditioning, should not really get here - ! - call psb_errpush(4001,name,a_err='mld_no_ml_ in mlprc_aply?') - goto 9999 - - case(mld_add_ml_) - ! - ! Additive multilevel - ! - - call add_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info) - - case(mld_mult_ml_) - ! - ! Multiplicative multilevel (multiplicative among the levels, additive inside - ! each level) - ! - ! Pre/post-smoothing versions. - ! Note that the transpose switches pre <-> post. - ! - - select case(p%precv(2)%iprcparm(mld_smoother_pos_)) - - case(mld_post_smooth_) + nlev = size(p%precv) + allocate(mlprec_wrk(nlev),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + level = 1 + + nc2l = psb_cd_get_local_cols(p%precv(level)%base_desc) + nr2l = psb_cd_get_local_rows(p%precv(level)%base_desc) + allocate(mlprec_wrk(level)%x2l(nc2l),mlprec_wrk(level)%y2l(nc2l),& + & stat=info) + if (info /= 0) then + info=4025 + call psb_errpush(info,name,i_err=(/size(x)+size(y),0,0,0,0/),& + & a_err='real(psb_dpk_)') + goto 9999 + end if - select case (trans_) - case('N') - call mlt_post_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info) - - case('T','C') - call mlt_pre_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info) - case default - info = 4001 - call psb_errpush(info,name,a_err='invalid trans') - goto 9999 - end select + mlprec_wrk(level)%x2l(:) = x(:) + mlprec_wrk(level)%y2l(:) = dzero - case(mld_pre_smooth_) + call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info) - select case (trans_) - case('N') - call mlt_pre_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info) - case('T','C') - call mlt_post_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info) - case default - info = 4001 - call psb_errpush(info,name,a_err='invalid trans') - goto 9999 - end select - - case(mld_twoside_smooth_) + if (info /= 0) then + call psb_errpush(4001,name,a_err='Inner prec aply') + goto 9999 + end if - call mlt_twoside_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info) + call psb_geaxpby(alpha,mlprec_wrk(level)%y2l,beta,y,& + & p%precv(level)%base_desc,info) - case default - info = 4013 - call psb_errpush(info,name,a_err='invalid smooth_pos',& - & i_Err=(/p%precv(2)%iprcparm(mld_smoother_pos_),0,0,0,0/)) - goto 9999 - - end select + if (info /= 0) then + call psb_errpush(4001,name,a_err='Error final update') + goto 9999 + end if - case default - info = 4013 - call psb_errpush(info,name,a_err='invalid mltype',& - & i_Err=(/p%precv(2)%iprcparm(mld_ml_type_),0,0,0,0/)) - goto 9999 - end select - endif call psb_erractionrestore(err_act) return @@ -383,9 +480,14 @@ contains end if - call mld_baseprec_aply(done,p%precv(level)%prec,& +!!$ call mld_baseprec_aply(done,p%precv(level)%prec,& +!!$ & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,& +!!$ & p%precv(level)%base_desc, trans,work,info) + call p%precv(level)%sm%apply(done,& & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,& & p%precv(level)%base_desc, trans,work,info) + + if (info /= 0) goto 9999 if (level < nlev) then call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) @@ -448,20 +550,17 @@ contains & work=work,trans=trans) if (info /= 0) goto 9999 - - call mld_baseprec_aply(done,p%precv(level)%prec,& + call p%precv(level)%sm%apply(done,& & mlprec_wrk(level)%x2l,done,mlprec_wrk(level)%y2l,& - & p%precv(level)%base_desc,& - & trans,work,info) + & p%precv(level)%base_desc, trans,work,info) else - call mld_baseprec_aply(done,p%precv(level)%prec,& + + call p%precv(level)%sm%apply(done,& & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,& - & p%precv(level)%base_desc,& - & trans,work,info) + & p%precv(level)%base_desc, trans,work,info) end if - case('T','C') ! Post-smoothing transpose is pre-smoothing @@ -484,18 +583,20 @@ contains ! ! Apply the base preconditioner ! - call mld_baseprec_aply(done,p%precv(ilev)%prec,mlprec_wrk(ilev)%x2l,& - & dzero,mlprec_wrk(ilev)%y2l,p%precv(ilev)%base_desc,trans,work,info) + call p%precv(level)%sm%apply(done,& + & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,& + & p%precv(level)%base_desc, trans,work,info) + if (info /= 0) goto 9999 ! ! Compute the residual (at all levels but the coarsest one) ! - if (ilev < nlev) then - call psb_spmm(-done,p%precv(ilev)%base_a,& - & mlprec_wrk(ilev)%y2l,done,mlprec_wrk(ilev)%x2l,& - & p%precv(ilev)%base_desc,info,work=work,trans=trans) + if (level < nlev) then + call psb_spmm(-done,p%precv(level)%base_a,& + & mlprec_wrk(level)%y2l,done,mlprec_wrk(level)%x2l,& + & p%precv(level)%base_desc,info,work=work,trans=trans) if (info /= 0) goto 9999 call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) if (info /= 0) goto 9999 @@ -536,10 +637,10 @@ contains ! ! Apply the base preconditioner ! - call mld_baseprec_aply(done,p%precv(level)%prec,& - & mlprec_wrk(level)%x2l,& - & dzero,mlprec_wrk(level)%y2l,& - & p%precv(level)%base_desc,trans,work,info) + call p%precv(level)%sm%apply(done,& + & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,& + & p%precv(level)%base_desc, trans,work,info) + if (info /= 0) goto 9999 @@ -598,17 +699,15 @@ contains & work=work,trans=trans) if (info /= 0) goto 9999 - - call mld_baseprec_aply(done,p%precv(level)%prec,& + call p%precv(level)%sm%apply(done,& & mlprec_wrk(level)%x2l,done,mlprec_wrk(level)%y2l,& - & p%precv(level)%base_desc,& - & trans,work,info) + & p%precv(level)%base_desc, trans,work,info) else - call mld_baseprec_aply(done,p%precv(level)%prec,& + call p%precv(level)%sm%apply(done,& & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,& - & p%precv(level)%base_desc,& - & trans,work,info) + & p%precv(level)%base_desc, trans,work,info) + end if case default @@ -616,9 +715,9 @@ contains call psb_errpush(info,name,a_err='invalid trans') goto 9999 end select - + case(mld_twoside_smooth_) - + nc2l = psb_cd_get_local_cols(p%precv(level)%base_desc) nr2l = psb_cd_get_local_rows(p%precv(level)%base_desc) allocate(mlprec_wrk(level)%ty(nc2l), mlprec_wrk(level)%tx(nc2l), stat=info) @@ -645,9 +744,10 @@ contains ! ! Apply the base preconditioner ! - if (info == 0) call mld_baseprec_aply(done,p%precv(level)%prec,& + if (info == 0) call p%precv(level)%sm%apply(done,& & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,& - & p%precv(level)%base_desc,trans,work,info) + & p%precv(level)%base_desc, trans,work,info) + ! ! Compute the residual (at all levels but the coarsest one) ! and call recursively @@ -682,10 +782,11 @@ contains ! ! Apply the base preconditioner ! - if (info == 0) call mld_baseprec_aply(done,p%precv(level)%prec,& - & mlprec_wrk(level)%tx,& - & done,mlprec_wrk(level)%y2l,p%precv(level)%base_desc,& - & trans, work,info) + if (info == 0) call p%precv(level)%sm%apply(done,& + & mlprec_wrk(level)%tx,done,mlprec_wrk(level)%y2l,& + & p%precv(level)%base_desc, trans,work,info) + + if (info /= 0) then call psb_errpush(4001,name,a_err='Error: residual/baseprec_aply') goto 9999 @@ -723,1085 +824,5 @@ contains end subroutine inner_ml_aply - ! - ! Subroutine: add_ml_aply - ! Version: real - ! Note: internal subroutine of mld_dmlprec_aply. - ! - ! This routine computes - ! - ! Y = beta*Y + alpha*op(M^(-1))*X, - ! where - ! - M is an additive multilevel domain decomposition (Schwarz) preconditioner - ! associated to a certain matrix A and stored in p, - ! - op(M^(-1)) is M^(-1) or its transpose, according to the value of trans, - ! - X and Y are vectors, - ! - alpha and beta are scalars. - ! - ! The preconditioner M is additive both through the levels and inside each - ! level. - ! - ! For each level we have as many submatrices as processes (except for the coarsest - ! level where we might have a replicated index space) and each process takes care - ! of one submatrix. - ! - ! The multilevel preconditioner is regarded as an array of 'one-level' data structures, - ! each containing the part of the preconditioner associated to a certain level - ! (for more details see the description of mld_Tonelev_type in mld_prec_type.f90). - ! For each level ilev, the 'base preconditioner' K(ilev) is stored in - ! p%precv(ilev)%prec - ! 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 - ! aggregation. - ! - ! The levels are numbered in increasing order starting from the finest one, i.e. - ! level 1 is the finest level and A(1) is the matrix A. - ! - ! For details on the additive multilevel Schwarz preconditioner see the - ! Algorithm 3.1.1 in the book: - ! B.F. Smith, P.E. Bjorstad & W.D. Gropp, - ! Domain decomposition: parallel multilevel methods for elliptic partial - ! differential equations, Cambridge University Press, 1996. - ! - ! For a description of the arguments see mld_dmlprec_aply. - ! - ! A sketch of the algorithm implemented in this routine is provided below. - ! (P(ilev) denotes the smoothed prolongator from level ilev to level - ! ilev-1, while PT(ilev) denotes its transpose, i.e. the corresponding - ! restriction operator from level ilev-1 to level ilev). - ! - ! 1. ! Apply the base preconditioner at level 1. - ! ! The sum over the subdomains is carried out in the - ! ! application of K(1). - ! X(1) = Xest - ! Y(1) = (K(1)^(-1))*X(1) - ! - ! 2. DO ilev=2,nlev - ! - ! ! Transfer X(ilev-1) to the next coarser level. - ! X(ilev) = PT(ilev)*X(ilev-1) - ! - ! ! Apply the base preconditioner at the current level. - ! ! The sum over the subdomains is carried out in the - ! ! application of K(ilev). - ! Y(ilev) = (K(ilev)^(-1))*X(ilev) - ! - ! ENDDO - ! - ! 3. DO ilev=nlev-1,1,-1 - ! - ! ! Transfer Y(ilev+1) to the next finer level. - ! Y(ilev) = Y(ilev) + P(ilev+1)*Y(ilev+1) - ! - ! ENDDO - ! - ! 4. Yext = beta*Yext + alpha*Y(1) - ! - subroutine add_ml_aply(alpha,p,x,beta,y,desc_data,trans,work,info) - - implicit none - - ! Arguments - type(psb_desc_type),intent(in) :: desc_data - type(mld_dprec_type), intent(in) :: p - real(psb_dpk_),intent(in) :: alpha,beta - real(psb_dpk_),intent(in) :: x(:) - real(psb_dpk_),intent(inout) :: y(:) - character, intent(in) :: trans - real(psb_dpk_),target :: work(:) - integer, intent(out) :: info - - ! Local variables - integer :: ictxt,np,me,i, nr2l,nc2l,err_act - integer :: debug_level, debug_unit - integer :: nlev, ilev - character(len=20) :: name - -!!$ type psb_mlprec_wrk_type -!!$ real(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) -!!$ end type psb_mlprec_wrk_type - type(psb_mlprec_wrk_type), allocatable :: mlprec_wrk(:) - - name='add_ml_aply' - info = 0 - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - - ictxt = psb_cd_get_context(desc_data) - call psb_info(ictxt, me, np) - - if (debug_level >= psb_debug_inner_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' Entry ', size(p%precv) - - nlev = size(p%precv) - allocate(mlprec_wrk(nlev),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - ! - ! STEP 1 - ! - ! Apply the base preconditioner at the finest level - ! - allocate(mlprec_wrk(1)%x2l(size(x)),mlprec_wrk(1)%y2l(size(y)), stat=info) - if (info /= 0) then - info=4025 - call psb_errpush(info,name,i_err=(/size(x)+size(y),0,0,0,0/),& - & a_err='real(psb_dpk_)') - goto 9999 - end if - - mlprec_wrk(1)%x2l(:) = x(:) - mlprec_wrk(1)%y2l(:) = dzero - ! - ! STEP 2 - ! - ! For each level except the finest one ... - ! - do ilev = 1, nlev - if (ilev > 1) then - nc2l = psb_cd_get_local_cols(p%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),& - & stat=info) - if (info /= 0) then - info=4025 - call psb_errpush(info,name,i_err=(/2*nc2l,0,0,0,0/),& - & a_err='real(psb_dpk_)') - goto 9999 - end if - - ! Apply the prolongator transpose, i.e. the restriction - call psb_map_X2Y(done,mlprec_wrk(ilev-1)%x2l,& - & dzero,mlprec_wrk(ilev)%x2l,& - & p%precv(ilev)%map,info,work=work) - - if (info /=0) then - call psb_errpush(4001,name,a_err='Error during restriction') - goto 9999 - end if - end if - - ! - ! Apply the base preconditioner - ! - call mld_baseprec_aply(done,p%precv(ilev)%prec,& - & mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%y2l,& - & p%precv(ilev)%base_desc, trans,work,info) - - enddo - - ! - ! STEP 3 - ! - ! For each level except the finest one ... - ! - do ilev =nlev,2,-1 - - nc2l = psb_cd_get_local_cols(p%precv(ilev)%base_desc) - nr2l = psb_cd_get_local_rows(p%precv(ilev)%base_desc) - - ! - ! Apply the prolongator - ! - call psb_map_Y2X(done,mlprec_wrk(ilev)%y2l,& - & done,mlprec_wrk(ilev-1)%y2l,& - & p%precv(ilev)%map,info,work=work) - - if (info /=0) then - call psb_errpush(4001,name,a_err='Error during prolongation') - goto 9999 - end if - end do - - ! - ! STEP 4 - ! - ! Compute the output vector Y - ! - call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,p%precv(1)%base_desc,info) - if (info /= 0) then - call psb_errpush(4001,name,a_err='Error on final update') - goto 9999 - end if - - deallocate(mlprec_wrk,stat=info) - if (info /= 0) then - call psb_errpush(4000,name) - goto 9999 - end if - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act.eq.psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine add_ml_aply - ! - ! Subroutine: mlt_pre_ml_aply - ! Version: real - ! Note: internal subroutine of mld_dmlprec_aply. - ! - ! This routine computes - ! - ! Y = beta*Y + alpha*op(M^(-1))*X, - ! where - ! - M is a hybrid multilevel domain decomposition (Schwarz) preconditioner - ! 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, - ! - X and Y are vectors, - ! - alpha and beta are scalars. - ! - ! The preconditioner M is hybrid in the sense that it is multiplicative through the - ! levels and additive inside a level; pre-smoothing only is applied at each level. - ! - ! For each level we have as many submatrices as processes (except for the coarsest - ! level where we might have a replicated index space) and each process takes care - ! of one submatrix. - ! - ! The multilevel preconditioner is regarded as an array of 'one-level' data structures, - ! each containing the part of the preconditioner associated to a certain level - ! (for more details see the description of mld_Tonelev_type in mld_prec_type.f90). - ! For each level ilev, the 'base preconditioner' K(ilev) is stored in - ! p%precv(ilev)%prec - ! 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 - ! aggregation. - ! - ! The levels are numbered in increasing order starting from the finest one, i.e. - ! level 1 is the finest level and A(1) is the matrix A. - ! - ! For details on the pre-smoothed hybrid multiplicative multilevel Schwarz - ! preconditioner, see the Algorithm 3.2.1 in the book: - ! B.F. Smith, P.E. Bjorstad & W.D. Gropp, - ! Domain decomposition: parallel multilevel methods for elliptic partial - ! differential equations, Cambridge University Press, 1996. - ! - ! For a description of the arguments see mld_dmlprec_aply. - ! - ! A sketch of the algorithm implemented in this routine is provided below. - ! (P(ilev) denotes the smoothed prolongator from level ilev to level - ! ilev-1, while PT(ilev) denotes its transpose, i.e. the corresponding - ! restriction operator from level ilev-1 to level ilev). - ! - ! 1. X(1) = Xext - ! - ! 2. ! Apply the base preconditioner at the finest level. - ! Y(1) = (K(1)^(-1))*X(1) - ! - ! 3. ! Compute the residual at the finest level. - ! TX(1) = X(1) - A(1)*Y(1) - ! - ! 4. DO ilev=2, nlev - ! - ! ! Transfer the residual to the current (coarser) level. - ! X(ilev) = PT(ilev)*TX(ilev-1) - ! - ! ! Apply the base preconditioner at the current level. - ! ! The sum over the subdomains is carried out in the - ! ! application of K(ilev). - ! Y(ilev) = (K(ilev)^(-1))*X(ilev) - ! - ! ! Compute the residual at the current level (except at - ! ! the coarsest level). - ! IF (ilev < nlev) - ! TX(ilev) = (X(ilev)-A(ilev)*Y(ilev)) - ! - ! ENDDO - ! - ! 5. DO ilev=nlev-1,1,-1 - ! - ! ! Transfer Y(ilev+1) to the next finer level - ! Y(ilev) = Y(ilev) + P(ilev+1)*Y(ilev+1) - ! - ! ENDDO - ! - ! 6. Yext = beta*Yext + alpha*Y(1) - ! - ! - subroutine mlt_pre_ml_aply(alpha,p,x,beta,y,desc_data,trans,work,info) - - implicit none - - ! Arguments - type(psb_desc_type),intent(in) :: desc_data - type(mld_dprec_type), intent(in) :: p - real(psb_dpk_),intent(in) :: alpha,beta - real(psb_dpk_),intent(in) :: x(:) - real(psb_dpk_),intent(inout) :: y(:) - character, intent(in) :: trans - real(psb_dpk_),target :: work(:) - integer, intent(out) :: info - - ! Local variables - integer :: ictxt,np,me,i, nr2l,nc2l,err_act - integer :: debug_level, debug_unit - integer :: nlev, ilev - character(len=20) :: name - -!!$ type psb_mlprec_wrk_type -!!$ real(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) -!!$ end type psb_mlprec_wrk_type - type(psb_mlprec_wrk_type), allocatable :: mlprec_wrk(:) - - name='mlt_pre_ml_aply' - info = 0 - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - - ictxt = psb_cd_get_context(desc_data) - call psb_info(ictxt, me, np) - - if (debug_level >= psb_debug_inner_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' Entry ', size(p%precv) - - nlev = size(p%precv) - allocate(mlprec_wrk(nlev),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - ! - ! STEP 1 - ! - ! Copy the input vector X - ! - nc2l = psb_cd_get_local_cols(p%precv(1)%base_desc) - - allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), & - & mlprec_wrk(1)%tx(nc2l), stat=info) - if (info /= 0) then - info=4025 - call psb_errpush(info,name,i_err=(/4*nc2l,0,0,0,0/),& - & a_err='real(psb_dpk_)') - goto 9999 - end if - - mlprec_wrk(1)%x2l(:) = x - ! - ! STEP 2 - ! - ! Apply the base preconditioner at the finest level - ! - call mld_baseprec_aply(done,p%precv(1)%prec,mlprec_wrk(1)%x2l,& - & dzero,mlprec_wrk(1)%y2l,p%precv(1)%base_desc,& - & trans,work,info) - if (info /=0) then - call psb_errpush(4010,name,a_err=' baseprec_aply') - goto 9999 - end if - - ! - ! STEP 3 - ! - ! Compute the residual at the finest level - ! - mlprec_wrk(1)%tx = mlprec_wrk(1)%x2l - - call psb_spmm(-done,p%precv(1)%base_a,mlprec_wrk(1)%y2l,& - & done,mlprec_wrk(1)%tx,p%precv(1)%base_desc,info,& - & work=work,trans=trans) - if (info /=0) then - call psb_errpush(4001,name,a_err=' fine level residual') - goto 9999 - end if - - ! - ! STEP 4 - ! - ! For each level but the finest one ... - ! - do ilev = 2, nlev - - nc2l = psb_cd_get_local_cols(p%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),& - & mlprec_wrk(ilev)%x2l(nc2l), stat=info) - if (info /= 0) then - info=4025 - call psb_errpush(info,name,i_err=(/4*nc2l,0,0,0,0/),& - & a_err='real(psb_dpk_)') - goto 9999 - end if - - ! Apply the prolongator transpose, i.e. the restriction - call psb_map_X2Y(done,mlprec_wrk(ilev-1)%tx,& - & dzero,mlprec_wrk(ilev)%x2l,& - & p%precv(ilev)%map,info,work=work) - - if (info /=0) then - call psb_errpush(4001,name,a_err='Error during restriction') - goto 9999 - end if - - ! - ! Apply the base preconditioner - ! - call mld_baseprec_aply(done,p%precv(ilev)%prec,mlprec_wrk(ilev)%x2l,& - & dzero,mlprec_wrk(ilev)%y2l,p%precv(ilev)%base_desc,trans,work,info) - - ! - ! Compute the residual (at all levels but the coarsest one) - ! - if (ilev < nlev) then - mlprec_wrk(ilev)%tx = mlprec_wrk(ilev)%x2l - if (info == 0) call psb_spmm(-done,p%precv(ilev)%base_a,& - & mlprec_wrk(ilev)%y2l,done,mlprec_wrk(ilev)%tx,& - & p%precv(ilev)%base_desc,info,work=work,trans=trans) - endif - if (info /=0) then - call psb_errpush(4001,name,a_err='Error on up sweep residual') - goto 9999 - end if - enddo - - ! - ! STEP 5 - ! - ! For each level but the coarsest one ... - ! - do ilev = nlev-1, 1, -1 - ! - ! Apply the prolongator - ! - call psb_map_Y2X(done,mlprec_wrk(ilev+1)%y2l,& - & done,mlprec_wrk(ilev)%y2l,& - & p%precv(ilev+1)%map,info,work=work) - - if (info /=0) then - call psb_errpush(4001,name,a_err='Error during prolongation') - goto 9999 - end if - enddo - - ! - ! STEP 6 - ! - ! Compute the output vector Y - ! - call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,& - & p%precv(1)%base_desc,info) - if (info /=0) then - call psb_errpush(4001,name,a_err='Error on final update') - goto 9999 - end if - - deallocate(mlprec_wrk,stat=info) - if (info /= 0) then - call psb_errpush(4000,name) - goto 9999 - end if - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act.eq.psb_act_abort_) then - call psb_error() - return - end if - return - end subroutine mlt_pre_ml_aply - ! - ! Subroutine: mlt_post_ml_aply - ! Version: real - ! Note: internal subroutine of mld_dmlprec_aply. - ! - ! This routine computes - ! - ! Y = beta*Y + alpha*op(M^(-1))*X, - ! where - ! - M is a hybrid multilevel domain decomposition (Schwarz) preconditioner - ! 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, - ! - X and Y are vectors, - ! - alpha and beta are scalars. - ! - ! The preconditioner M is hybrid in the sense that it is multiplicative through the - ! levels and additive inside a level; post-smoothing only is applied at each level. - ! - ! For each level we have as many submatrices as processes (except for the coarsest - ! level where we might have a replicated index space) and each process takes care - ! of one submatrix. - ! - ! The multilevel preconditioner is regarded as an array of 'one-level' data structures, - ! each containing the part of the preconditioner associated to a certain level - ! (for more details see the description of mld_Tonelev_type in mld_prec_type.f90). - ! For each level ilev, the 'base preconditioner' K(ilev) is stored in - ! p%precv(ilev)%prec - ! 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 - ! aggregation. - ! - ! The levels are numbered in increasing order starting from the finest one, i.e. - ! level 1 is the finest level and A(1) is the matrix A. - ! - ! For details on hybrid multiplicative multilevel Schwarz preconditioners, see - ! B.F. Smith, P.E. Bjorstad & W.D. Gropp, - ! Domain decomposition: parallel multilevel methods for elliptic partial - ! differential equations, Cambridge University Press, 1996. - ! - ! For a description of the arguments see mld_dmlprec_aply. - ! - ! A sketch of the algorithm implemented in this routine is provided below. - ! (P(ilev) denotes the smoothed prolongator from level ilev to level - ! ilev-1, while PT(ilev) denotes its transpose, i.e. the corresponding - ! restriction operator from level ilev-1 to level ilev). - ! - ! 1. X(1) = Xext - ! - ! 2. DO ilev=2, nlev - ! - ! ! Transfer X(ilev-1) to the next coarser level. - ! X(ilev) = PT(ilev)*X(ilev-1) - ! - ! ENDDO - ! - ! 3.! Apply the preconditioner at the coarsest level. - ! Y(nlev) = (K(nlev)^(-1))*X(nlev) - ! - ! 4. DO ilev=nlev-1,1,-1 - ! - ! ! Transfer Y(ilev+1) to the next finer level. - ! Y(ilev) = P(ilev+1)*Y(ilev+1) - ! - ! ! Compute the residual at the current level and apply to it the - ! ! base preconditioner. The sum over the subdomains is carried out - ! ! in the application of K(ilev). - ! Y(ilev) = Y(ilev) + (K(ilev)^(-1))*(X(ilev)-A(ilev)*Y(ilev)) - ! - ! ENDDO - ! - ! 5. Yext = beta*Yext + alpha*Y(1) - ! - ! - subroutine mlt_post_ml_aply(alpha,p,x,beta,y,desc_data,trans,work,info) - - implicit none - - ! Arguments - type(psb_desc_type),intent(in) :: desc_data - type(mld_dprec_type), intent(in) :: p - real(psb_dpk_),intent(in) :: alpha,beta - real(psb_dpk_),intent(in) :: x(:) - real(psb_dpk_),intent(inout) :: y(:) - character, intent(in) :: trans - real(psb_dpk_),target :: work(:) - integer, intent(out) :: info - - ! Local variables - integer :: ictxt,np,me,i, nr2l,nc2l,err_act - integer :: debug_level, debug_unit - integer :: nlev, ilev - character(len=20) :: name - -!!$ type psb_mlprec_wrk_type -!!$ real(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) -!!$ end type psb_mlprec_wrk_type - type(psb_mlprec_wrk_type), allocatable :: mlprec_wrk(:) - - name='mlt_post_ml_aply' - info = 0 - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - - ictxt = psb_cd_get_context(desc_data) - call psb_info(ictxt, me, np) - - if (debug_level >= psb_debug_inner_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' Entry ', size(p%precv) - - nlev = size(p%precv) - allocate(mlprec_wrk(nlev),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - ! - ! STEP 1 - ! - ! Copy the input vector X - ! - if (debug_level >= psb_debug_inner_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' desc_data status',allocated(desc_data%matrix_data) - - nc2l = psb_cd_get_local_cols(p%precv(1)%base_desc) - - allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), & - & mlprec_wrk(1)%tx(nc2l), stat=info) - if (info /= 0) then - info=4025 - call psb_errpush(info,name,i_err=(/4*nc2l,0,0,0,0/),& - & a_err='real(psb_dpk_)') - goto 9999 - end if - - call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%tx,& - & p%precv(1)%base_desc,info) - call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%x2l,& - & p%precv(1)%base_desc,info) - - ! - ! STEP 2 - ! - ! For each level but the finest one ... - ! - do ilev=2, nlev - - nc2l = psb_cd_get_local_cols(p%precv(ilev)%base_desc) - nr2l = psb_cd_get_local_rows(p%precv(ilev)%base_desc) - - if (debug_level >= psb_debug_inner_) & - & write(debug_unit,*) me,' ',trim(name), & - & ' starting up sweep ',& - & ilev,allocated(p%precv(ilev)%iprcparm),nc2l, nr2l - - allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%y2l(nc2l),& - & mlprec_wrk(ilev)%x2l(nc2l), stat=info) - - if (info /= 0) then - info=4025 - call psb_errpush(info,name,i_err=(/4*nc2l,0,0,0,0/),& - & a_err='real(psb_dpk_)') - goto 9999 - end if - - ! Apply the prolongator transpose, i.e. the restriction - call psb_map_X2Y(done,mlprec_wrk(ilev-1)%x2l,& - & dzero,mlprec_wrk(ilev)%x2l,& - & p%precv(ilev)%map,info,work=work) - - if (info /=0) then - call psb_errpush(4001,name,a_err='Error during restriction') - goto 9999 - end if - - ! - ! update x2l - ! - call psb_geaxpby(done,mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%tx,& - & p%precv(ilev)%base_desc,info) - if (info /= 0) then - call psb_errpush(4001,name,a_err='Error in update') - goto 9999 - end if - - if (debug_level >= psb_debug_inner_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' done up sweep ', ilev - - enddo - - ! - ! STEP 3 - ! - ! Apply the base preconditioner at the coarsest level - ! - call mld_baseprec_aply(done,p%precv(nlev)%prec,mlprec_wrk(nlev)%x2l, & - & dzero, mlprec_wrk(nlev)%y2l,p%precv(nlev)%base_desc,trans,work,info) - if (info /=0) then - call psb_errpush(4010,name,a_err='baseprec_aply') - goto 9999 - end if - - if (debug_level >= psb_debug_inner_) write(debug_unit,*) & - & me,' ',trim(name), ' done baseprec_aply ', nlev - - ! - ! STEP 4 - ! - ! For each level but the coarsest one ... - ! - do ilev=nlev-1, 1, -1 - - if (debug_level >= psb_debug_inner_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' starting down sweep',ilev - - ! - ! Apply the prolongator - ! - call psb_map_Y2X(done,mlprec_wrk(ilev+1)%y2l,& - & dzero,mlprec_wrk(ilev)%y2l,& - & p%precv(ilev+1)%map,info,work=work) - - if (info /=0) then - call psb_errpush(4001,name,a_err='Error during prolongation') - goto 9999 - end if - - ! - ! Compute the residual - ! - call psb_spmm(-done,p%precv(ilev)%base_a,mlprec_wrk(ilev)%y2l,& - & done,mlprec_wrk(ilev)%tx,p%precv(ilev)%base_desc,info,& - & work=work,trans=trans) - - ! - ! Apply the base preconditioner - ! - if (info == 0) call mld_baseprec_aply(done,p%precv(ilev)%prec,& - & mlprec_wrk(ilev)%tx,done,mlprec_wrk(ilev)%y2l,p%precv(ilev)%base_desc,& - & trans,work,info) - if (info /=0) then - call psb_errpush(4001,name,a_err=' spmm/baseprec_aply') - goto 9999 - end if - - if (debug_level >= psb_debug_inner_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' done down sweep',ilev - enddo - - ! - ! STEP 5 - ! - ! Compute the output vector Y - ! - call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,p%precv(1)%base_desc,info) - - if (info /=0) then - call psb_errpush(4001,name,a_err=' Final update') - goto 9999 - end if - - - - deallocate(mlprec_wrk,stat=info) - if (info /= 0) then - call psb_errpush(4000,name) - goto 9999 - end if - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act.eq.psb_act_abort_) then - call psb_error() - return - end if - return - end subroutine mlt_post_ml_aply - ! - ! Subroutine: mlt_twoside_ml_aply - ! Version: real - ! Note: internal subroutine of mld_dmlprec_aply. - ! - ! This routine computes - ! - ! Y = beta*Y + alpha*op(M^(-1))*X, - ! where - ! - M is a symmetrized hybrid multilevel domain decomposition (Schwarz) - ! preconditioner 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, - ! - X and Y are vectors, - ! - alpha and beta are scalars. - ! - ! The preconditioner M is hybrid in the sense that it is multiplicative through - ! the levels and additive inside a level; it is symmetrized since pre-smoothing - ! and post-smoothing are applied at each level. - ! - ! For each level we have as many submatrices as processes (except for the coarsest - ! level where we might have a replicated index space) and each process takes care - ! of one submatrix. - ! - ! The multilevel preconditioner is regarded as an array of 'one-level' data structures, - ! each containing the part of the preconditioner associated to a certain level - ! (for more details see the description of mld_Tonelev_type in mld_prec_type.f90). - ! For each level ilev, the 'base preconditioner' K(ilev) is stored in - ! p%precv(ilev)%prec - ! 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 - ! aggregation. - ! - ! The levels are numbered in increasing order starting from the finest one, i.e. - ! level 1 is the finest level and A(1) is the matrix A. - ! - ! For details on the symmetrized hybrid multiplicative multilevel Schwarz - ! preconditioner, see the Algorithm 3.2.2 of the book: - ! B.F. Smith, P.E. Bjorstad & W.D. Gropp, - ! Domain decomposition: parallel multilevel methods for elliptic partial - ! differential equations, Cambridge University Press, 1996. - ! - ! For a description of the arguments see mld_dmlprec_aply. - ! - ! A sketch of the algorithm implemented in this routine is provided below. - ! (P(ilev) denotes the smoothed prolongator from level ilev to level - ! ilev-1, while PT(ilev) denotes its transpose, i.e. the corresponding - ! restriction operator from level ilev-1 to level ilev). - ! - ! 1. X(1) = Xext - ! - ! 2. ! Apply the base peconditioner at the finest level - ! Y(1) = (K(1)^(-1))*X(1) - ! - ! 3. ! Compute the residual at the finest level - ! TX(1) = X(1) - A(1)*Y(1) - ! - ! 4. DO ilev=2, nlev - ! - ! ! Transfer the residual to the current (coarser) level - ! X(ilev) = PT(ilev)*TX(ilev-1) - ! - ! ! Apply the base preconditioner at the current level. - ! ! The sum over the subdomains is carried out in the - ! ! application of K(ilev) - ! Y(ilev) = (K(ilev)^(-1))*X(ilev) - ! - ! ! Compute the residual at the current level - ! ! (except for ilev=nlev) - ! if(ilev < nlev)then - ! TX(ilev) = (X(ilev)-A(ilev)*Y(ilev)) - ! endif - ! - ! ENDDO - ! - ! 5. DO ilev=NLEV-1,1,-1 - ! - ! ! Transfer Y(ilev+1) to the next finer level - ! Y(ilev) = Y(ilev) + P(ilev+1)*Y(ilev+1) - ! - ! ! Compute the residual at the current level and apply to it the - ! ! base preconditioner. The sum over the subdomains is carried out - ! ! in the application of K(ilev) - ! Y(ilev) = Y(ilev) + (K(ilev)**(-1))*(X(ilev)-A(ilev)*Y(ilev)) - ! - ! ENDDO - ! - ! 6. Yext = beta*Yext + alpha*Y(1) - ! - subroutine mlt_twoside_ml_aply(alpha,p,x,beta,y,desc_data,trans,work,info) - - implicit none - - ! Arguments - type(psb_desc_type),intent(in) :: desc_data - type(mld_dprec_type), intent(in) :: p - real(psb_dpk_),intent(in) :: alpha,beta - real(psb_dpk_),intent(in) :: x(:) - real(psb_dpk_),intent(inout) :: y(:) - character, intent(in) :: trans - real(psb_dpk_),target :: work(:) - integer, intent(out) :: info - - ! Local variables - integer :: ictxt,np,me,i, nr2l,nc2l,err_act - integer :: debug_level, debug_unit - integer :: nlev, ilev - character(len=20) :: name - -!!$ type psb_mlprec_wrk_type -!!$ real(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) -!!$ end type psb_mlprec_wrk_type - type(psb_mlprec_wrk_type), allocatable :: mlprec_wrk(:) - - name='mlt_twoside_ml_aply' - info = 0 - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - - ictxt = psb_cd_get_context(desc_data) - call psb_info(ictxt, me, np) - - if (debug_level >= psb_debug_inner_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' Entry ', size(p%precv) - - nlev = size(p%precv) - allocate(mlprec_wrk(nlev),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - ! STEP 1 - ! - ! Copy the input vector X - ! - nc2l = psb_cd_get_local_cols(p%precv(1)%base_desc) - - allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), & - & mlprec_wrk(1)%ty(nc2l), mlprec_wrk(1)%tx(nc2l), stat=info) - - if (info /= 0) then - info=4025 - call psb_errpush(info,name,i_err=(/4*nc2l,0,0,0,0/),& - & a_err='real(psb_dpk_)') - goto 9999 - end if - - call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%x2l,& - & p%precv(1)%base_desc,info) - call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%tx,& - & p%precv(1)%base_desc,info) - - ! - ! STEP 2 - ! - ! Apply the base preconditioner at the finest level - ! - call mld_baseprec_aply(done,p%precv(1)%prec,mlprec_wrk(1)%x2l,& - & dzero,mlprec_wrk(1)%y2l,p%precv(1)%base_desc,& - & trans,work,info) - ! - ! STEP 3 - ! - ! Compute the residual at the finest level - ! - mlprec_wrk(1)%ty = mlprec_wrk(1)%x2l - if (info == 0) call psb_spmm(-done,p%precv(1)%base_a,mlprec_wrk(1)%y2l,& - & done,mlprec_wrk(1)%ty,p%precv(1)%base_desc,info,& - & work=work,trans=trans) - if (info /=0) then - call psb_errpush(4010,name,a_err='Fine level baseprec/residual') - goto 9999 - end if - - ! - ! STEP 4 - ! - ! For each level but the finest one ... - ! - do ilev = 2, nlev - - nc2l = psb_cd_get_local_cols(p%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),& - & mlprec_wrk(ilev)%y2l(nc2l),mlprec_wrk(ilev)%x2l(nc2l), stat=info) - - if (info /= 0) then - info=4025 - call psb_errpush(info,name,i_err=(/4*nc2l,0,0,0,0/),& - & a_err='real(psb_dpk_)') - goto 9999 - end if - - ! Apply the prolongator transpose, i.e. the restriction - call psb_map_X2Y(done,mlprec_wrk(ilev-1)%ty,& - & dzero,mlprec_wrk(ilev)%x2l,& - & p%precv(ilev)%map,info,work=work) - - if (info /=0) then - call psb_errpush(4001,name,a_err='Error during restriction') - goto 9999 - end if - - call psb_geaxpby(done,mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%tx,& - & p%precv(ilev)%base_desc,info) - ! - ! Apply the base preconditioner - ! - if (info == 0) call mld_baseprec_aply(done,p%precv(ilev)%prec,& - & mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%y2l,& - & p%precv(ilev)%base_desc,trans,work,info) - ! - ! Compute the residual (at all levels but the coarsest one) - ! - if(ilev < nlev) then - mlprec_wrk(ilev)%ty = mlprec_wrk(ilev)%x2l - if (info == 0) call psb_spmm(-done,p%precv(ilev)%base_a,& - & mlprec_wrk(ilev)%y2l,done,mlprec_wrk(ilev)%ty,& - & p%precv(ilev)%base_desc,info,work=work,trans=trans) - endif - if (info /=0) then - call psb_errpush(4001,name,a_err='baseprec_aply/residual') - goto 9999 - end if - - enddo - - ! - ! STEP 5 - ! - ! For each level but the coarsest one ... - ! - do ilev=nlev-1, 1, -1 - - ! - ! Apply the prolongator - ! - call psb_map_Y2X(done,mlprec_wrk(ilev+1)%y2l,& - & done,mlprec_wrk(ilev)%y2l,& - & p%precv(ilev+1)%map,info,work=work) - - if (info /=0 ) then - call psb_errpush(4001,name,a_err='Error during restriction') - goto 9999 - end if - - ! - ! Compute the residual - ! - call psb_spmm(-done,p%precv(ilev)%base_a,mlprec_wrk(ilev)%y2l,& - & done,mlprec_wrk(ilev)%tx,p%precv(ilev)%base_desc,info,& - & work=work,trans=trans) - ! - ! Apply the base preconditioner - ! - if (info == 0) call mld_baseprec_aply(done,p%precv(ilev)%prec,mlprec_wrk(ilev)%tx,& - & done,mlprec_wrk(ilev)%y2l,p%precv(ilev)%base_desc, trans, work,info) - if (info /= 0) then - call psb_errpush(4001,name,a_err='Error: residual/baseprec_aply') - goto 9999 - end if - enddo - - ! - ! STEP 6 - ! - ! Compute the output vector Y - ! - call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,& - & p%precv(1)%base_desc,info) - - if (info /= 0) then - call psb_errpush(4001,name,a_err='Error final update') - goto 9999 - end if - - - - deallocate(mlprec_wrk,stat=info) - if (info /= 0) then - call psb_errpush(4000,name) - goto 9999 - end if - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act.eq.psb_act_abort_) then - call psb_error() - return - end if - return - end subroutine mlt_twoside_ml_aply - end subroutine mld_dmlprec_aply diff --git a/mlprec/mld_dmlprec_bld.f90 b/mlprec/mld_dmlprec_bld.f90 index 2f802cac..300d1b33 100644 --- a/mlprec/mld_dmlprec_bld.f90 +++ b/mlprec/mld_dmlprec_bld.f90 @@ -69,6 +69,10 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info) use psb_base_mod use mld_inner_mod, mld_protect_name => mld_dmlprec_bld use mld_prec_mod + use mld_d_jac_smoother + use mld_d_as_smoother + use mld_d_diag_solver + use mld_d_ilu_solver Implicit None @@ -306,8 +310,60 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info) call mld_check_def(p%precv(i)%prec%iprcparm(mld_smoother_sweeps_),& & 'Jacobi sweeps',1,is_legal_jac_sweeps) - call mld_baseprec_bld(p%precv(i)%base_a,p%precv(i)%base_desc,& - & p%precv(i)%prec,info) + ! + ! Test version for beginning of OO stuff. + ! + if (allocated(p%precv(i)%sm)) then + call p%precv(i)%sm%free(info) + if (info ==0) deallocate(p%precv(i)%sm,stat=info) + if (info /= 0) then + call psb_errpush(4000,name,a_err='One level preconditioner build.') + goto 9999 + endif + end if + select case (p%precv(i)%prec%iprcparm(mld_smoother_type_)) + case(mld_diag_, mld_bjac_, mld_pjac_) + allocate(mld_d_jac_smoother_type :: p%precv(i)%sm, stat=info) + case(mld_as_) + allocate(mld_d_as_smoother_type :: p%precv(i)%sm, stat=info) + case default + info = -1 + end select + if (info /= 0) then + write(0,*) ' Smoother allocation error',info,& + & p%precv(i)%prec%iprcparm(mld_smoother_type_) + call psb_errpush(4001,name,a_err='One level preconditioner build.') + goto 9999 + endif + call p%precv(i)%sm%set(mld_sub_restr_,p%precv(i)%prec%iprcparm(mld_sub_restr_),info) + call p%precv(i)%sm%set(mld_sub_prol_,p%precv(i)%prec%iprcparm(mld_sub_prol_),info) + call p%precv(i)%sm%set(mld_sub_ovr_,p%precv(i)%prec%iprcparm(mld_sub_ovr_),info) + call p%precv(i)%sm%set(mld_smoother_sweeps_,& + & p%precv(i)%prec%iprcparm(mld_smoother_sweeps_),info) + + select case (p%precv(i)%prec%iprcparm(mld_sub_solve_)) + case(mld_ilu_n_,mld_milu_n_,mld_ilu_t_) + allocate(mld_d_ilu_solver_type :: p%precv(i)%sm%sv, stat=info) + if (info == 0) call p%precv(i)%sm%sv%set(mld_sub_solve_,& + & p%precv(i)%prec%iprcparm(mld_sub_solve_),info) + if (info == 0) call p%precv(i)%sm%sv%set(mld_sub_fillin_,& + & p%precv(i)%prec%iprcparm(mld_sub_fillin_),info) + if (info == 0) call p%precv(i)%sm%sv%set(mld_sub_iluthrs_,& + & p%precv(i)%prec%rprcparm(mld_sub_iluthrs_),info) + case(mld_diag_scale_) + allocate(mld_d_diag_solver_type :: p%precv(i)%sm%sv, stat=info) + case default + info = -1 + end select + + if (info /= 0) then + write(0,*) ' Solver allocation error',info,& + & p%precv(i)%prec%iprcparm(mld_sub_solve_) + call psb_errpush(4001,name,a_err='One level preconditioner build.') + goto 9999 + endif + + call p%precv(i)%sm%build(p%precv(i)%base_a,p%precv(i)%base_desc,'F',info) if (info /= 0) then call psb_errpush(4001,name,a_err='One level preconditioner build.') diff --git a/mlprec/mld_dprecset.f90 b/mlprec/mld_dprecset.f90 index e70b5684..38ba5b5d 100644 --- a/mlprec/mld_dprecset.f90 +++ b/mlprec/mld_dprecset.f90 @@ -140,7 +140,10 @@ subroutine mld_dprecseti(p,what,val,info,ilev) ! Rules for fine level are slightly different. ! select case(what) - case(mld_smoother_type_,mld_sub_solve_,mld_sub_restr_,mld_sub_prol_,& + case(mld_smoother_type_) + p%precv(ilev_)%iprcparm(what) = val + p%precv(ilev_)%prec%iprcparm(what) = val + case(mld_sub_solve_,mld_sub_restr_,mld_sub_prol_,& & mld_sub_ren_,mld_sub_ovr_,mld_sub_fillin_,mld_smoother_sweeps_) p%precv(ilev_)%prec%iprcparm(what) = val case(mld_ml_type_,mld_aggr_alg_,mld_aggr_kind_,& @@ -153,7 +156,10 @@ subroutine mld_dprecseti(p,what,val,info,ilev) else if (ilev_ > 1) then select case(what) - case(mld_smoother_type_,mld_sub_solve_,mld_sub_restr_,mld_sub_prol_,& + case(mld_smoother_type_) + p%precv(ilev_)%iprcparm(what) = val + p%precv(ilev_)%prec%iprcparm(what) = val + case(mld_sub_solve_,mld_sub_restr_,mld_sub_prol_,& & mld_sub_ren_,mld_sub_ovr_,mld_sub_fillin_,& & mld_smoother_sweeps_) p%precv(ilev_)%prec%iprcparm(what) = val @@ -220,21 +226,32 @@ subroutine mld_dprecseti(p,what,val,info,ilev) ! levels ! select case(what) - case(mld_smoother_type_,mld_sub_solve_,mld_sub_restr_,mld_sub_prol_,& + case(mld_sub_solve_,mld_sub_restr_,mld_sub_prol_,& & mld_sub_ren_,mld_sub_ovr_,mld_sub_fillin_,& & mld_smoother_sweeps_) do ilev_=1,max(1,nlev_-1) if (.not.allocated(p%precv(ilev_)%iprcparm)) then write(0,*) name,& - &': Error: uninitialized preconditioner component,',& - &' should call MLD_PRECINIT' + &': Error: uninitialized preconditioner component, should call MLD_PRECINIT' + info = -1 + return + endif + p%precv(ilev_)%iprcparm(what) = val + p%precv(ilev_)%prec%iprcparm(what) = val + end do + case(mld_smoother_type_) + do ilev_=1,nlev_ + if (.not.allocated(p%precv(ilev_)%iprcparm)) then + write(0,*) name,& + &': Error: uninitialized preconditioner component, should call MLD_PRECINIT' info = -1 return endif + p%precv(ilev_)%iprcparm(what) = val p%precv(ilev_)%prec%iprcparm(what) = val end do case(mld_ml_type_,mld_aggr_alg_,mld_aggr_kind_,& - & mld_smoother_pos_,mld_aggr_omega_alg_,mld_aggr_eig_) + & mld_smoother_pos_,mld_aggr_omega_alg_,mld_aggr_eig_,mld_aggr_filter_) do ilev_=1,nlev_ if (.not.allocated(p%precv(ilev_)%iprcparm)) then write(0,*) name,& @@ -265,15 +282,18 @@ subroutine mld_dprecseti(p,what,val,info,ilev) endif if (nlev_ > 1) then - p%precv(nlev_)%iprcparm(mld_coarse_solve_) = val + p%precv(nlev_)%iprcparm(mld_coarse_solve_) = val p%precv(nlev_)%prec%iprcparm(mld_smoother_type_) = mld_bjac_ - p%precv(nlev_)%iprcparm(mld_coarse_mat_) = mld_distr_mat_ + p%precv(nlev_)%iprcparm(mld_coarse_mat_) = mld_distr_mat_ select case (val) case(mld_umf_, mld_slu_) - p%precv(nlev_)%iprcparm(mld_coarse_mat_) = mld_repl_mat_ + p%precv(nlev_)%iprcparm(mld_coarse_mat_) = mld_repl_mat_ p%precv(nlev_)%prec%iprcparm(mld_sub_solve_) = val case(mld_sludist_) p%precv(nlev_)%prec%iprcparm(mld_sub_solve_) = val +!!$ case(mld_jac_) +!!$ p%precv(nlev_)%prec%iprcparm(mld_smoother_type_) = mld_jac_ +!!$ p%precv(nlev_)%prec%iprcparm(mld_sub_solve_) = mld_diag_scale_ end select endif case(mld_coarse_subsolve_) diff --git a/mlprec/mld_move_alloc.f90 b/mlprec/mld_move_alloc_mod.f90 similarity index 100% rename from mlprec/mld_move_alloc.f90 rename to mlprec/mld_move_alloc_mod.f90 diff --git a/tests/pdegen/ppde.f90 b/tests/pdegen/ppde.f90 index 235b3065..e0ccc108 100644 --- a/tests/pdegen/ppde.f90 +++ b/tests/pdegen/ppde.f90 @@ -104,11 +104,13 @@ program ppde character(len=20) :: descr ! verbose description of the prec character(len=10) :: prec ! overall prectype integer :: novr ! number of overlap layers + integer :: jsweeps ! Jacobi/smoother sweeps character(len=16) :: restr ! restriction over application of as character(len=16) :: prol ! prolongation over application of as - character(len=16) :: solve ! Factorization type: ILU, SuperLU, UMFPACK. + character(len=16) :: solve ! Solver type: ILU, SuperLU, UMFPACK. integer :: fill1 ! Fill-in for factorization 1 real(psb_dpk_) :: thr1 ! Threshold for fact. 1 ILU(T) + character(len=16) :: smther ! Smoother integer :: nlev ! Number of levels in multilevel prec. character(len=16) :: aggrkind ! smoothed/raw aggregatin character(len=16) :: aggr_alg ! local or global aggregation @@ -171,17 +173,15 @@ program ppde if (psb_toupper(prectype%prec) =='ML') then nlv = prectype%nlev - else - nlv = 1 - end if - call mld_precinit(prec,prectype%prec,info,nlev=nlv) - call mld_precset(prec,mld_sub_ovr_,prectype%novr,info) - call mld_precset(prec,mld_sub_restr_,prectype%restr,info) - call mld_precset(prec,mld_sub_prol_,prectype%prol,info) - call mld_precset(prec,mld_sub_solve_,prectype%solve,info) - call mld_precset(prec,mld_sub_fillin_,prectype%fill1,info) - call mld_precset(prec,mld_sub_iluthrs_,prectype%thr1,info) - if (psb_toupper(prectype%prec) =='ML') then + call mld_precinit(prec,prectype%prec, info, nlev=nlv) + call mld_precset(prec,mld_smoother_type_, prectype%smther, info) + call mld_precset(prec,mld_smoother_sweeps_, prectype%jsweeps, info) + call mld_precset(prec,mld_sub_ovr_, prectype%novr, info) + call mld_precset(prec,mld_sub_restr_, prectype%restr, info) + call mld_precset(prec,mld_sub_prol_, prectype%prol, info) + call mld_precset(prec,mld_sub_solve_, prectype%solve, info) + call mld_precset(prec,mld_sub_fillin_, prectype%fill1, info) + call mld_precset(prec,mld_sub_iluthrs_, prectype%thr1, info) call mld_precset(prec,mld_aggr_kind_, prectype%aggrkind,info) call mld_precset(prec,mld_aggr_alg_, prectype%aggr_alg,info) call mld_precset(prec,mld_ml_type_, prectype%mltype, info) @@ -193,8 +193,17 @@ program ppde call mld_precset(prec,mld_coarse_fillin_, prectype%cfill, info) call mld_precset(prec,mld_coarse_iluthrs_, prectype%cthres, info) call mld_precset(prec,mld_coarse_sweeps_, prectype%cjswp, info) - end if - + else + nlv = 1 + call mld_precinit(prec,prectype%prec, info, nlev=nlv) + call mld_precset(prec,mld_smoother_sweeps_, prectype%jsweeps, info) + call mld_precset(prec,mld_sub_ovr_, prectype%novr, info) + call mld_precset(prec,mld_sub_restr_, prectype%restr, info) + call mld_precset(prec,mld_sub_prol_, prectype%prol, info) + call mld_precset(prec,mld_sub_solve_, prectype%solve, info) + call mld_precset(prec,mld_sub_fillin_, prectype%fill1, info) + call mld_precset(prec,mld_sub_iluthrs_, prectype%thr1, info) + end if call psb_barrier(ictxt) t1 = psb_wtime() call mld_precbld(a,desc_a,prec,info) @@ -305,7 +314,9 @@ contains call read_data(prectype%solve,5) ! Factorization type: ILU, SuperLU, UMFPACK. call read_data(prectype%fill1,5) ! Fill-in for factorization 1 call read_data(prectype%thr1,5) ! Threshold for fact. 1 ILU(T) + call read_data(prectype%jsweeps,5) ! Jacobi sweeps for PJAC if (psb_toupper(prectype%prec) == 'ML') then + call read_data(prectype%smther,5) ! Smoother type. call read_data(prectype%nlev,5) ! Number of levels in multilevel prec. call read_data(prectype%aggrkind,5) ! smoothed/raw aggregatin call read_data(prectype%aggr_alg,5) ! local or global aggregation @@ -340,7 +351,9 @@ contains call psb_bcast(ictxt,prectype%solve) ! Factorization type: ILU, SuperLU, UMFPACK. call psb_bcast(ictxt,prectype%fill1) ! Fill-in for factorization 1 call psb_bcast(ictxt,prectype%thr1) ! Threshold for fact. 1 ILU(T) + call psb_bcast(ictxt,prectype%jsweeps) ! Threshold for fact. 1 ILU(T) if (psb_toupper(prectype%prec) == 'ML') then + call psb_bcast(ictxt,prectype%smther) ! Smoother type. call psb_bcast(ictxt,prectype%nlev) ! Number of levels in multilevel prec. call psb_bcast(ictxt,prectype%aggrkind) ! smoothed/raw aggregatin call psb_bcast(ictxt,prectype%aggr_alg) ! local or global aggregation @@ -352,7 +365,7 @@ contains call psb_bcast(ictxt,prectype%cfill) ! Fill-in for factorization 1 call psb_bcast(ictxt,prectype%cthres) ! Threshold for fact. 1 ILU(T) call psb_bcast(ictxt,prectype%cjswp) ! Jacobi sweeps - call psb_bcast(ictxt,prectype%athres) ! smoother aggr thresh + call psb_bcast(ictxt,prectype%athres) ! smoother aggr thresh end if if (iam==psb_root_) then diff --git a/tests/pdegen/runs/ppde.inp b/tests/pdegen/runs/ppde.inp index 91998620..b86a2da6 100644 --- a/tests/pdegen/runs/ppde.inp +++ b/tests/pdegen/runs/ppde.inp @@ -2,26 +2,28 @@ BICGSTAB ! Iterative method: BiCGSTAB BiCG CGS RGMRES BiCGSTA CSR ! Storage format CSR COO JAD 40 ! IDIM; domain size is idim**3 2 ! ISTOPC -00100 ! ITMAX +0010 ! ITMAX 01 ! ITRACE 30 ! IRST (restart for RGMRES and BiCGSTABL) -1.d-7 ! EPS -RAS ! Longer descriptive name for preconditioner (up to 20 chars) -AS ! Preconditioner NONE DIAG BJAC AS ML -0 ! Number of overlap layers for AS preconditioner at finest level +1.d-6 ! EPS +3L-M-RAS-I-D4 ! Longer descriptive name for preconditioner (up to 20 chars) +ML ! Preconditioner NONE JACOBI BJAC AS ML +1 ! Number of overlap layers for AS preconditioner at finest level HALO ! Restriction operator NONE HALO NONE ! Prolongation operator NONE SUM AVG -ILU ! Subdomain solver ILU MILU ILUT UMF SLU -1 ! Level-set N for ILU(N) +ILU ! Subdomain solver DSCALE ILU MILU ILUT UMF SLU +0 ! Level-set N for ILU(N) 1.d-4 ! Threshold T for ILU(T,P) +1 ! Smoother/Jacobi sweeps +AS ! Smoother type JACOBI BJAC AS ignored for non-ML 3 ! Number of levels in a multilevel preconditioner SMOOTHED ! Kind of aggregation: SMOOTHED, NONSMOOTHED, MINENERGY DEC ! Type of aggregation DEC SYMDEC GLB MULT ! Type of multilevel correction: ADD MULT -TWOSIDE ! Side of mult. correction PRE POST TWOSIDE (ignored for ADD) -REPL ! Coarse level: matrix distribution DIST REPL -BJAC ! Coarse level: solver BJAC UMF SLU SLUDIST -ILUT ! Coarse level: subsolver ILU UMF SLU SLUDIST +POST ! Side of multiplicative correction PRE POST TWOSIDE (ignored for ADD) +DIST ! Coarse level: matrix distribution DIST REPL +BJAC ! Coarse level: solver JACOBI BJAC UMF SLU SLUDIST +ILU ! Coarse level: subsolver DSCALE ILU UMF SLU SLUDIST 1 ! Coarse level: Level-set N for ILU(N) 1.d-4 ! Coarse level: Threshold T for ILU(T,P) 4 ! Coarse level: Number of Jacobi sweeps