diff --git a/mlprec/Makefile b/mlprec/Makefile index a29dbbdc..ac269caa 100644 --- a/mlprec/Makefile +++ b/mlprec/Makefile @@ -29,6 +29,7 @@ INNEROBJS= mld_dcoarse_bld.o mld_dmlprec_bld.o mld_dslu_bld.o mld_dumf_bld.o \ mld_cmlprec_aply.o mld_cslud_bld.o mld_caggrmat_asb.o \ mld_zcoarse_bld.o mld_zmlprec_bld.o mld_zslu_bld.o mld_zumf_bld.o \ mld_zilu0_fact.o mld_ziluk_fact.o mld_zilut_fact.o mld_zaggrmap_bld.o \ + mld_zmlprec_aply.o mld_zslud_bld.o mld_zaggrmat_asb.o \ $(MPFOBJS) # diff --git a/mlprec/mld_cmlprec_aply.f90 b/mlprec/mld_cmlprec_aply.f90 index 09a2bedc..ad9665ab 100644 --- a/mlprec/mld_cmlprec_aply.f90 +++ b/mlprec/mld_cmlprec_aply.f90 @@ -357,7 +357,8 @@ subroutine mld_cmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info) nlev = size(p%precv) allocate(mlprec_wrk(nlev),stat=info) if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='Allocate') goto 9999 end if level = 1 @@ -379,7 +380,8 @@ subroutine mld_cmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info) call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Inner prec aply') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Inner prec aply') goto 9999 end if @@ -387,7 +389,8 @@ subroutine mld_cmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info) & p%precv(level)%base_desc,info) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error final update') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error final update') goto 9999 end if @@ -432,7 +435,8 @@ contains nlev = size(p%precv) if ((level < 1) .or. (level > nlev)) then - call psb_errpush(psb_err_internal_error_,name,a_err='wrong call level to inner_ml') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong call level to inner_ml') goto 9999 end if ictxt = psb_cd_get_context(p%precv(level)%base_desc) @@ -459,7 +463,8 @@ contains ! No preconditioning, should not really get here ! write(0,*) 'MLD_NO_ML_ in inner_ml ',level - call psb_errpush(psb_err_internal_error_,name,a_err='mld_no_ml_ in mlprc_aply?') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='mld_no_ml_ in mlprc_aply?') goto 9999 case(mld_add_ml_) @@ -474,7 +479,8 @@ contains & p%precv(level)%map,info,work=work) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') goto 9999 end if @@ -523,7 +529,8 @@ contains & p%precv(level)%map,info,work=work) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') goto 9999 end if end if @@ -574,7 +581,8 @@ contains & p%precv(level)%map,info,work=work) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') goto 9999 end if @@ -634,7 +642,8 @@ contains & p%precv(level)%map,info,work=work) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') goto 9999 end if @@ -686,7 +695,8 @@ contains & p%precv(level)%map,info,work=work) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') goto 9999 end if @@ -748,7 +758,8 @@ contains & p%precv(level)%map,info,work=work) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') goto 9999 end if end if @@ -791,7 +802,8 @@ contains & p%precv(level+1)%map,info,work=work) if (info /= psb_success_ ) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') goto 9999 end if @@ -814,7 +826,8 @@ contains & p%precv(level)%base_desc, trans,& & sweeps,work,info) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error: residual/baseprec_aply') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error: residual/baseprec_aply') goto 9999 end if diff --git a/mlprec/mld_dmlprec_aply.f90 b/mlprec/mld_dmlprec_aply.f90 index 32f5a7cf..82ba5158 100644 --- a/mlprec/mld_dmlprec_aply.f90 +++ b/mlprec/mld_dmlprec_aply.f90 @@ -379,7 +379,8 @@ subroutine mld_dmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info) call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Inner prec aply') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Inner prec aply') goto 9999 end if @@ -387,7 +388,8 @@ subroutine mld_dmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info) & p%precv(level)%base_desc,info) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error final update') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error final update') goto 9999 end if @@ -432,7 +434,8 @@ contains nlev = size(p%precv) if ((level < 1) .or. (level > nlev)) then - call psb_errpush(psb_err_internal_error_,name,a_err='wrong call level to inner_ml') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong call level to inner_ml') goto 9999 end if ictxt = psb_cd_get_context(p%precv(level)%base_desc) @@ -459,7 +462,8 @@ contains ! No preconditioning, should not really get here ! write(0,*) 'MLD_NO_ML_ in inner_ml ',level - call psb_errpush(psb_err_internal_error_,name,a_err='mld_no_ml_ in mlprc_aply?') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='mld_no_ml_ in mlprc_aply?') goto 9999 case(mld_add_ml_) @@ -474,7 +478,8 @@ contains & p%precv(level)%map,info,work=work) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') goto 9999 end if @@ -524,7 +529,8 @@ contains & p%precv(level)%map,info,work=work) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') goto 9999 end if end if @@ -575,7 +581,8 @@ contains & p%precv(level)%map,info,work=work) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') goto 9999 end if @@ -635,7 +642,8 @@ contains & p%precv(level)%map,info,work=work) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') goto 9999 end if @@ -687,7 +695,8 @@ contains & p%precv(level)%map,info,work=work) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') goto 9999 end if @@ -749,7 +758,8 @@ contains & p%precv(level)%map,info,work=work) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') goto 9999 end if end if @@ -792,7 +802,8 @@ contains & p%precv(level+1)%map,info,work=work) if (info /= psb_success_ ) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') goto 9999 end if @@ -815,7 +826,8 @@ contains & p%precv(level)%base_desc, trans,& & sweeps,work,info) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error: residual/baseprec_aply') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error: residual/baseprec_aply') goto 9999 end if diff --git a/mlprec/mld_smlprec_aply.f90 b/mlprec/mld_smlprec_aply.f90 index 35de2390..c335352b 100644 --- a/mlprec/mld_smlprec_aply.f90 +++ b/mlprec/mld_smlprec_aply.f90 @@ -379,7 +379,8 @@ subroutine mld_smlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info) call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Inner prec aply') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Inner prec aply') goto 9999 end if @@ -387,7 +388,8 @@ subroutine mld_smlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info) & p%precv(level)%base_desc,info) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error final update') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error final update') goto 9999 end if @@ -432,7 +434,8 @@ contains nlev = size(p%precv) if ((level < 1) .or. (level > nlev)) then - call psb_errpush(psb_err_internal_error_,name,a_err='wrong call level to inner_ml') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong call level to inner_ml') goto 9999 end if ictxt = psb_cd_get_context(p%precv(level)%base_desc) @@ -459,7 +462,8 @@ contains ! No preconditioning, should not really get here ! write(0,*) 'MLD_NO_ML_ in inner_ml ',level - call psb_errpush(psb_err_internal_error_,name,a_err='mld_no_ml_ in mlprc_aply?') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='mld_no_ml_ in mlprc_aply?') goto 9999 case(mld_add_ml_) @@ -474,7 +478,8 @@ contains & p%precv(level)%map,info,work=work) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') goto 9999 end if @@ -524,7 +529,8 @@ contains & p%precv(level)%map,info,work=work) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') goto 9999 end if end if @@ -575,7 +581,8 @@ contains & p%precv(level)%map,info,work=work) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') goto 9999 end if @@ -635,7 +642,8 @@ contains & p%precv(level)%map,info,work=work) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') goto 9999 end if @@ -687,7 +695,8 @@ contains & p%precv(level)%map,info,work=work) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') goto 9999 end if @@ -749,7 +758,8 @@ contains & p%precv(level)%map,info,work=work) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') goto 9999 end if end if @@ -792,7 +802,8 @@ contains & p%precv(level+1)%map,info,work=work) if (info /= psb_success_ ) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') goto 9999 end if @@ -815,7 +826,8 @@ contains & p%precv(level)%base_desc, trans,& & sweeps,work,info) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error: residual/baseprec_aply') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error: residual/baseprec_aply') goto 9999 end if diff --git a/mlprec/mld_zmlprec_aply.f90 b/mlprec/mld_zmlprec_aply.f90 index 89075b0d..eb72ae31 100644 --- a/mlprec/mld_zmlprec_aply.f90 +++ b/mlprec/mld_zmlprec_aply.f90 @@ -138,6 +138,179 @@ ! 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 only +! +! 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 only +! +! +! 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 +! +! For details on the hybrid multiplicative multilevel Schwarz preconditioner +! with pre- and post-smoothing (symmetrized multiplicative multilevel), 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_zmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info) @@ -158,9 +331,13 @@ subroutine mld_zmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info) ! Local variables integer :: ictxt, np, me, err_act - integer :: debug_level, debug_unit + integer :: debug_level, debug_unit, nlev,nc2l,nr2l,level character(len=20) :: name character :: trans_ + type psb_mlprec_wrk_type + complex(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) + end type psb_mlprec_wrk_type + type(psb_mlprec_wrk_type), allocatable :: mlprec_wrk(:) name = 'mld_zmlprec_aply' info = psb_success_ @@ -177,78 +354,43 @@ subroutine mld_zmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info) trans_ = psb_toupper(trans) - select case(p%precv(2)%iprcparm(mld_ml_type_)) - - case(mld_no_ml_) - ! - ! No preconditioning, should not really get here - ! - call psb_errpush(psb_err_internal_error_,name,a_err='mld_no_ml_ in mlprc_aply?') + nlev = size(p%precv) + allocate(mlprec_wrk(nlev),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,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 /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/size(x)+size(y),0,0,0,0/),& + & a_err='real(psb_dpk_)') + goto 9999 + end if - 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_) - - 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 = psb_err_internal_error_ - call psb_errpush(info,name,a_err='invalid trans') - goto 9999 - end select - - case(mld_pre_smooth_) - - 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 = psb_err_internal_error_ - call psb_errpush(info,name,a_err='invalid trans') - goto 9999 - end select + mlprec_wrk(level)%x2l(:) = x(:) + mlprec_wrk(level)%y2l(:) = zzero - case(mld_twoside_smooth_) + call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info) - call mlt_twoside_ml_aply(alpha,p,x,beta,y,desc_data,trans_,work,info) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,a_err='Inner prec aply') + goto 9999 + end if - case default - info = psb_err_from_subroutine_ai_ - 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 + call psb_geaxpby(alpha,mlprec_wrk(level)%y2l,beta,y,& + & p%precv(level)%base_desc,info) - end select - - case default - info = psb_err_from_subroutine_ai_ - call psb_errpush(info,name,a_err='invalid mltype',& - & i_Err=(/p%precv(2)%iprcparm(mld_ml_type_),0,0,0,0/)) - goto 9999 + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,a_err='Error final update') + goto 9999 + end if - end select call psb_erractionrestore(err_act) return @@ -262,157 +404,46 @@ subroutine mld_zmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info) return contains - ! - ! Subroutine: add_ml_aply - ! Version: complex - ! 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 (conjugate) 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_zmlprec_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) = 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) + + recursive subroutine inner_ml_aply(level,p,mlprec_wrk,trans,work,info) implicit none ! Arguments - type(psb_desc_type),intent(in) :: desc_data - type(mld_zprec_type), intent(in) :: p - complex(psb_dpk_),intent(in) :: alpha,beta - complex(psb_dpk_),intent(in) :: x(:) - complex(psb_dpk_),intent(inout) :: y(:) - character, intent(in) :: trans - complex(psb_dpk_),target :: work(:) - integer, intent(out) :: info + integer :: level + type(mld_zprec_type), intent(in) :: p + type(psb_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:) + character, intent(in) :: trans + complex(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 + integer :: nlev, ilev, sweeps - type psb_mlprec_wrk_type - complex(psb_dpk_), allocatable :: tx(:),ty(:),x2l(:),y2l(:) - end type psb_mlprec_wrk_type - type(psb_mlprec_wrk_type), allocatable :: mlprec_wrk(:) + character(len=20) :: name - name = 'add_ml_aply' + name = 'inner_ml_aply' info = psb_success_ 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 /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,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 /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/size(x)+size(y),0,0,0,0/),& - & a_err='complex(psb_dpk_)') + if ((level < 1) .or. (level > nlev)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong call level to inner_ml') goto 9999 end if + ictxt = psb_cd_get_context(p%precv(level)%base_desc) + call psb_info(ictxt, me, np) - mlprec_wrk(1)%x2l(:) = x(:) - mlprec_wrk(1)%y2l(:) = zzero - call mld_baseprec_aply(alpha,p%precv(1)%prec,x,beta,y,& - & p%precv(1)%base_desc,trans,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='baseprec_aply') - goto 9999 - end if - ! - ! STEP 2 - ! - ! For each level except 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)%x2l(nc2l),mlprec_wrk(ilev)%y2l(nc2l),& + if (level > 1) then + 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 /= psb_success_) then info=psb_err_alloc_request_ @@ -420,924 +451,400 @@ contains & a_err='complex(psb_dpk_)') goto 9999 end if - - ! Apply the prolongator transpose, i.e. the restriction - call psb_map_X2Y(zone,mlprec_wrk(ilev-1)%x2l,& - & zzero,mlprec_wrk(ilev)%x2l,& - & p%precv(ilev)%map,info,work=work) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') - goto 9999 - end if - - ! - ! Apply the base preconditioner - ! - call mld_baseprec_aply(zone,p%precv(ilev)%prec,& - & mlprec_wrk(ilev)%x2l,zzero,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(zone,mlprec_wrk(ilev)%y2l,& - & zone,mlprec_wrk(ilev-1)%y2l,& - & p%precv(ilev)%map,info,work=work) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,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,zone,y,p%precv(1)%base_desc,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error on final update') - goto 9999 - end if - - deallocate(mlprec_wrk,stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_dealloc_,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: complex - ! Note: internal subroutine of mld_zmlprec_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 (conjugate) 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_zmlprec_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_zprec_type), intent(in) :: p - complex(psb_dpk_),intent(in) :: alpha,beta - complex(psb_dpk_),intent(in) :: x(:) - complex(psb_dpk_),intent(inout) :: y(:) - character, intent(in) :: trans - complex(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 - complex(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 = psb_success_ - 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 /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,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 /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/4*nc2l,0,0,0,0/),& - & a_err='complex(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(zone,p%precv(1)%prec,mlprec_wrk(1)%x2l,& - & zzero,mlprec_wrk(1)%y2l,p%precv(1)%base_desc,& - & trans,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err=' baseprec_aply') - goto 9999 - end if - - ! - ! STEP 3 - ! - ! Compute the residual at the finest level - ! - mlprec_wrk(1)%tx = mlprec_wrk(1)%x2l - - call psb_spmm(-zone,p%precv(1)%base_a,mlprec_wrk(1)%y2l,& - & zone,mlprec_wrk(1)%tx,p%precv(1)%base_desc,info,& - & work=work,trans=trans) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,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 /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/4*nc2l,0,0,0,0/),& - & a_err='complex(psb_dpk_)') - goto 9999 - end if - - ! Apply the prolongator transpose, i.e. the restriction - call psb_map_X2Y(zone,mlprec_wrk(ilev-1)%tx,& - & zzero,mlprec_wrk(ilev)%x2l,& - & p%precv(ilev)%map,info,work=work) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') - goto 9999 - end if - - ! - ! Apply the base preconditioner - ! - call mld_baseprec_aply(zone,p%precv(ilev)%prec,mlprec_wrk(ilev)%x2l,& - & zzero,mlprec_wrk(ilev)%y2l,p%precv(ilev)%base_desc,trans,work,info) + select case(p%precv(level)%iprcparm(mld_ml_type_)) + case(mld_no_ml_) ! - ! Compute the residual (at all levels but the coarsest one) - ! - if (ilev < nlev) then - mlprec_wrk(ilev)%tx = mlprec_wrk(ilev)%x2l - if (info == psb_success_) call psb_spmm(-zone,p%precv(ilev)%base_a,& - & mlprec_wrk(ilev)%y2l,zone,mlprec_wrk(ilev)%tx,& - & p%precv(ilev)%base_desc,info,work=work,trans=trans) - endif - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,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(zone,mlprec_wrk(ilev+1)%y2l,& - & zone,mlprec_wrk(ilev)%y2l,& - & p%precv(ilev+1)%map,info,work=work) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,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 /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error on final update') - goto 9999 - end if - - deallocate(mlprec_wrk,stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_dealloc_,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: complex - ! Note: internal subroutine of mld_zmlprec_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 (conjugate) 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_zmlprec_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_zprec_type), intent(in) :: p - complex(psb_dpk_),intent(in) :: alpha,beta - complex(psb_dpk_),intent(in) :: x(:) - complex(psb_dpk_),intent(inout) :: y(:) - character, intent(in) :: trans - complex(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 - complex(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 = psb_success_ - 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 /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + ! No preconditioning, should not really get here + ! + write(0,*) 'MLD_NO_ML_ in inner_ml ',level + call psb_errpush(psb_err_internal_error_,name,& + & a_err='mld_no_ml_ in mlprc_aply?') 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 /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/4*nc2l,0,0,0,0/),& - & a_err='complex(psb_dpk_)') - goto 9999 - end if - - call psb_geaxpby(zone,x,zzero,mlprec_wrk(1)%tx,& - & p%precv(1)%base_desc,info) - call psb_geaxpby(zone,x,zzero,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 /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/4*nc2l,0,0,0,0/),& - & a_err='complex(psb_dpk_)') - goto 9999 - end if - - ! Apply the prolongator transpose, i.e. the restriction - call psb_map_X2Y(zone,mlprec_wrk(ilev-1)%x2l,& - & zzero,mlprec_wrk(ilev)%x2l,& - & p%precv(ilev)%map,info,work=work) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') - goto 9999 - end if + case(mld_add_ml_) ! - ! update x2l + ! Additive multilevel ! - call psb_geaxpby(zone,mlprec_wrk(ilev)%x2l,zzero,mlprec_wrk(ilev)%tx,& - & p%precv(ilev)%base_desc,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,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(zone,p%precv(nlev)%prec,mlprec_wrk(nlev)%x2l, & - & zzero, mlprec_wrk(nlev)%y2l,p%precv(nlev)%base_desc,trans,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,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 + if (level > 1) then + ! Apply the restriction + call psb_map_X2Y(zone,mlprec_wrk(level-1)%x2l,& + & zzero,mlprec_wrk(level)%x2l,& + & p%precv(level)%map,info,work=work) - ! - ! Apply the prolongator - ! - call psb_map_Y2X(zone,mlprec_wrk(ilev+1)%y2l,& - & zzero,mlprec_wrk(ilev)%y2l,& - & p%precv(ilev+1)%map,info,work=work) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during prolongation') - goto 9999 - end if - - ! - ! Compute the residual - ! - call psb_spmm(-zone,p%precv(ilev)%base_a,mlprec_wrk(ilev)%y2l,& - & zone,mlprec_wrk(ilev)%tx,p%precv(ilev)%base_desc,info,& - & work=work,trans=trans) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') + goto 9999 + end if - ! - ! Apply the base preconditioner - ! - if (info == psb_success_) call mld_baseprec_aply(zone,p%precv(ilev)%prec,& - & mlprec_wrk(ilev)%tx,zone,mlprec_wrk(ilev)%y2l,p%precv(ilev)%base_desc,& - &trans,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,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 /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err=' Final update') - goto 9999 - end if - + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_) + call p%precv(level)%sm%apply(zone,& + & mlprec_wrk(level)%x2l,zzero,mlprec_wrk(level)%y2l,& + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) + + if (info /= psb_success_) goto 9999 + if (level < nlev) then + call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) + if (info /= psb_success_) goto 9999 + ! + ! Apply the prolongator + ! + call psb_map_Y2X(zone,mlprec_wrk(level+1)%y2l,& + & zone,mlprec_wrk(level)%y2l,& + & p%precv(level+1)%map,info,work=work) + if (info /= psb_success_) goto 9999 - - deallocate(mlprec_wrk,stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_dealloc_,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: complex - ! Note: internal subroutine of mld_zmlprec_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 (conjugate) 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_zmlprec_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_zprec_type), intent(in) :: p - complex(psb_dpk_),intent(in) :: alpha,beta - complex(psb_dpk_),intent(in) :: x(:) - complex(psb_dpk_),intent(inout) :: y(:) - character, intent(in) :: trans - complex(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 - complex(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 = psb_success_ - 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 /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,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 /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/4*nc2l,0,0,0,0/),& - & a_err='complex(psb_dpk_)') - goto 9999 - end if - - call psb_geaxpby(zone,x,zzero,mlprec_wrk(1)%x2l,& - & p%precv(1)%base_desc,info) - call psb_geaxpby(zone,x,zzero,mlprec_wrk(1)%tx,& - & p%precv(1)%base_desc,info) - - ! - ! STEP 2 - ! - ! Apply the base preconditioner at the finest level - ! - call mld_baseprec_aply(zone,p%precv(1)%prec,mlprec_wrk(1)%x2l,& - & zzero,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 == psb_success_) call psb_spmm(-zone,p%precv(1)%base_a,mlprec_wrk(1)%y2l,& - & zone,mlprec_wrk(1)%ty,p%precv(1)%base_desc,info,& - & work=work,trans=trans) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,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 /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/4*nc2l,0,0,0,0/),& - & a_err='complex(psb_dpk_)') - goto 9999 - end if - - ! Apply the prolongator transpose, i.e. the restriction - call psb_map_X2Y(zone,mlprec_wrk(ilev-1)%ty,& - & zzero,mlprec_wrk(ilev)%x2l,& - & p%precv(ilev)%map,info,work=work) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') - goto 9999 end if - call psb_geaxpby(zone,mlprec_wrk(ilev)%x2l,zzero,mlprec_wrk(ilev)%tx,& - & p%precv(ilev)%base_desc,info) + case(mld_mult_ml_) + ! + ! Multiplicative multilevel (multiplicative among the levels, additive inside + ! each level) ! - ! Apply the base preconditioner + ! Pre/post-smoothing versions. + ! Note that the transpose switches pre <-> post. ! - if (info == psb_success_) call mld_baseprec_aply(zone,p%precv(ilev)%prec,& - & mlprec_wrk(ilev)%x2l,zzero,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 == psb_success_) call psb_spmm(-zone,p%precv(ilev)%base_a,& - & mlprec_wrk(ilev)%y2l,zone,mlprec_wrk(ilev)%ty,& - & p%precv(ilev)%base_desc,info,work=work,trans=trans) - endif - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,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 + select case(p%precv(level)%iprcparm(mld_smoother_pos_)) + + case(mld_post_smooth_) + + select case (trans_) + case('N') + + if (level > 1) then + ! Apply the restriction + call psb_map_X2Y(zone,mlprec_wrk(level-1)%x2l,& + & zzero,mlprec_wrk(level)%x2l,& + & p%precv(level)%map,info,work=work) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') + goto 9999 + end if + end if + + ! This is one step of post-smoothing + + if (level < nlev) then + call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) + if (info /= psb_success_) goto 9999 + ! + ! Apply the prolongator + ! + call psb_map_Y2X(zone,mlprec_wrk(level+1)%y2l,& + & zzero,mlprec_wrk(level)%y2l,& + & p%precv(level+1)%map,info,work=work) + if (info /= psb_success_) goto 9999 + ! + ! Compute the residual + ! + call psb_spmm(-zone,p%precv(level)%base_a,mlprec_wrk(level)%y2l,& + & zone,mlprec_wrk(level)%x2l,p%precv(level)%base_desc,info,& + & work=work,trans=trans) + if (info /= psb_success_) goto 9999 + + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_post_) + call p%precv(level)%sm%apply(zone,& + & mlprec_wrk(level)%x2l,zone,mlprec_wrk(level)%y2l,& + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) + else + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_) + call p%precv(level)%sm%apply(zone,& + & mlprec_wrk(level)%x2l,zzero,mlprec_wrk(level)%y2l,& + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) + + end if + + case('T','C') + + ! Post-smoothing transpose is pre-smoothing + + + if (level > 1) then + ! Apply the restriction + call psb_map_X2Y(zone,mlprec_wrk(level-1)%x2l,& + & zzero,mlprec_wrk(level)%x2l,& + & p%precv(level)%map,info,work=work) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') + goto 9999 + end if + + + end if + + ! + ! Apply the base preconditioner + ! + if (level < nlev) then + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_post_) + else + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_) + end if + call p%precv(level)%sm%apply(zone,& + & mlprec_wrk(level)%x2l,zzero,mlprec_wrk(level)%y2l,& + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) + + if (info /= psb_success_) goto 9999 + + ! + ! Compute the residual (at all levels but the coarsest one) + ! + if (level < nlev) then + call psb_spmm(-zone,p%precv(level)%base_a,& + & mlprec_wrk(level)%y2l,zone,mlprec_wrk(level)%x2l,& + & p%precv(level)%base_desc,info,work=work,trans=trans) + if (info /= psb_success_) goto 9999 + call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) + if (info /= psb_success_) goto 9999 + + call psb_map_Y2X(zone,mlprec_wrk(level+1)%y2l,& + & zone,mlprec_wrk(level)%y2l,& + & p%precv(level+1)%map,info,work=work) + if (info /= psb_success_) goto 9999 + + end if + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='invalid trans') + goto 9999 + end select + + case(mld_pre_smooth_) + + select case (trans_) + case('N') + ! One step of pre-smoothing + + + if (level > 1) then + ! Apply the restriction + call psb_map_X2Y(zone,mlprec_wrk(level-1)%x2l,& + & zzero,mlprec_wrk(level)%x2l,& + & p%precv(level)%map,info,work=work) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') + goto 9999 + end if + + end if + + ! + ! Apply the base preconditioner + ! + if (level < nlev) then + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_pre_) + else + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_) + end if + call p%precv(level)%sm%apply(zone,& + & mlprec_wrk(level)%x2l,zzero,mlprec_wrk(level)%y2l,& + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) + + if (info /= psb_success_) goto 9999 + + ! + ! Compute the residual (at all levels but the coarsest one) + ! + if (level < nlev) then + call psb_spmm(-zone,p%precv(level)%base_a,& + & mlprec_wrk(level)%y2l,zone,mlprec_wrk(level)%x2l,& + & p%precv(level)%base_desc,info,work=work,trans=trans) + if (info /= psb_success_) goto 9999 + call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) + if (info /= psb_success_) goto 9999 + + call psb_map_Y2X(zone,mlprec_wrk(level+1)%y2l,& + & zone,mlprec_wrk(level)%y2l,& + & p%precv(level+1)%map,info,work=work) + if (info /= psb_success_) goto 9999 + + end if + + + case('T','C') + + ! pre-smooth transpose is post-smoothing + + + if (level > 1) then + ! Apply the restriction + call psb_map_X2Y(zone,mlprec_wrk(level-1)%x2l,& + & zzero,mlprec_wrk(level)%x2l,& + & p%precv(level)%map,info,work=work) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') + goto 9999 + end if + + end if + + if (level < nlev) then + call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) + if (info /= psb_success_) goto 9999 + ! + ! Apply the prolongator + ! + call psb_map_Y2X(zone,mlprec_wrk(level+1)%y2l,& + & zzero,mlprec_wrk(level)%y2l,& + & p%precv(level+1)%map,info,work=work) + if (info /= psb_success_) goto 9999 + ! + ! Compute the residual + ! + call psb_spmm(-zone,p%precv(level)%base_a,mlprec_wrk(level)%y2l,& + & zone,mlprec_wrk(level)%x2l,p%precv(level)%base_desc,info,& + & work=work,trans=trans) + if (info /= psb_success_) goto 9999 + + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_pre_) + call p%precv(level)%sm%apply(zone,& + & mlprec_wrk(level)%x2l,zone,mlprec_wrk(level)%y2l,& + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) + else + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_) + call p%precv(level)%sm%apply(zone,& + & mlprec_wrk(level)%x2l,zzero,mlprec_wrk(level)%y2l,& + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) + end if + + case default + info = psb_err_internal_error_ + 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) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/2*nc2l,0,0,0,0/),& + & a_err='real(psb_dpk_)') + goto 9999 + end if + + if (level > 1) then + ! Apply the restriction + call psb_map_X2Y(zone,mlprec_wrk(level-1)%ty,& + & zzero,mlprec_wrk(level)%x2l,& + & p%precv(level)%map,info,work=work) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') + goto 9999 + end if + end if + call psb_geaxpby(zone,mlprec_wrk(level)%x2l,zzero,mlprec_wrk(level)%tx,& + & p%precv(level)%base_desc,info) + ! + ! Apply the base preconditioner + ! + if (level < nlev) then + if (trans == 'N') then + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_pre_) + else + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_post_) + end if + else + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_) + end if + if (info == psb_success_) call p%precv(level)%sm%apply(zone,& + & mlprec_wrk(level)%x2l,zzero,mlprec_wrk(level)%y2l,& + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) + ! + ! Compute the residual (at all levels but the coarsest one) + ! and call recursively + ! + if(level < nlev) then + mlprec_wrk(level)%ty = mlprec_wrk(level)%x2l + if (info == psb_success_) call psb_spmm(-zone,p%precv(level)%base_a,& + & mlprec_wrk(level)%y2l,zone,mlprec_wrk(level)%ty,& + & p%precv(level)%base_desc,info,work=work,trans=trans) + + call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) + + + ! + ! Apply the prolongator + ! + call psb_map_Y2X(zone,mlprec_wrk(level+1)%y2l,& + & zone,mlprec_wrk(level)%y2l,& + & p%precv(level+1)%map,info,work=work) + + if (info /= psb_success_ ) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') + goto 9999 + end if + + ! + ! Compute the residual + ! + call psb_spmm(-zone,p%precv(level)%base_a,mlprec_wrk(level)%y2l,& + & zone,mlprec_wrk(level)%tx,p%precv(level)%base_desc,info,& + & work=work,trans=trans) + ! + ! Apply the base preconditioner + ! + if (trans == 'N') then + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_post_) + else + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_pre_) + end if + if (info == psb_success_) call p%precv(level)%sm%apply(zone,& + & mlprec_wrk(level)%tx,zone,mlprec_wrk(level)%y2l,& + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error: residual/baseprec_aply') + goto 9999 + end if + + endif - ! - ! Apply the prolongator - ! - call psb_map_Y2X(zone,mlprec_wrk(ilev+1)%y2l,& - & zone,mlprec_wrk(ilev)%y2l,& - & p%precv(ilev+1)%map,info,work=work) - - if (info /= psb_success_ ) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error during restriction') - goto 9999 - end if - - ! - ! Compute the residual - ! - call psb_spmm(-zone,p%precv(ilev)%base_a,mlprec_wrk(ilev)%y2l,& - & zone,mlprec_wrk(ilev)%tx,p%precv(ilev)%base_desc,info,& - & work=work,trans=trans) - ! - ! Apply the base preconditioner - ! - if (info == psb_success_) call mld_baseprec_aply(zone,p%precv(ilev)%prec,mlprec_wrk(ilev)%tx,& - & zone,mlprec_wrk(ilev)%y2l,p%precv(ilev)%base_desc, trans, work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,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 /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error final update') - goto 9999 - end if + case default + info = psb_err_from_subroutine_ai_ + call psb_errpush(info,name,a_err='invalid smooth_pos',& + & i_Err=(/p%precv(level)%iprcparm(mld_smoother_pos_),0,0,0,0/)) + goto 9999 + end select + case default + info = psb_err_from_subroutine_ai_ + call psb_errpush(info,name,a_err='invalid mltype',& + & i_Err=(/p%precv(level)%iprcparm(mld_ml_type_),0,0,0,0/)) + goto 9999 - deallocate(mlprec_wrk,stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_dealloc_,name) - goto 9999 - end if + end select call psb_erractionrestore(err_act) return @@ -1349,7 +856,8 @@ contains return end if return - end subroutine mlt_twoside_ml_aply + + end subroutine inner_ml_aply end subroutine mld_zmlprec_aply diff --git a/mlprec/mld_zslud_bld.f90 b/mlprec/mld_zslud_bld.f90 index ebce9562..2a2ffb08 100644 --- a/mlprec/mld_zslud_bld.f90 +++ b/mlprec/mld_zslud_bld.f90 @@ -84,7 +84,7 @@ subroutine mld_zsludist_bld(a,desc_a,p,info) & mglob,ifrst,ibcheck,nrow,ncol,npr,npc character(len=20) :: name, ch_err - if(psb_get_errstatus().ne.0) return + if (psb_get_errstatus().ne.0) return info=psb_success_ name='mld_zslud_bld' call psb_erractionsave(err_act) @@ -93,47 +93,56 @@ subroutine mld_zsludist_bld(a,desc_a,p,info) call psb_info(ictxt, me, np) - if (psb_toupper(a%fida) /= 'CSR') then + select type(aa=>a%a) + type is (psb_z_csr_sparse_mat) + + ! + ! WARN: we need to check for a BLOCK distribution (this is the + ! distribution required by SuperLU_DIST) + ! + nrow = psb_cd_get_local_rows(desc_a) + ncol = psb_cd_get_local_cols(desc_a) + call psb_loc_to_glob(1,ifrst,desc_a,info) + call psb_loc_to_glob(nrow,ibcheck,desc_a,info) + ibcheck = ibcheck - ifrst + 1 + ibcheck = ibcheck - nrow + call psb_amx(ictxt,ibcheck) + if (ibcheck > 0) then + write(0,*) 'Warning: does not look like a BLOCK distribution' + info=psb_err_unsupported_format_ + ch_err = aa%get_fmt() + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + + mglob = psb_cd_get_global_rows(desc_a) + nzt = aa%get_nzeros() + + npr = np + npc = 1 + call psb_loc_to_glob(aa%ja(1:nzt),desc_a,info,iact='I') + + ! + ! Compute the LU factorization + ! + call mld_zsludist_fact(mglob,nrow,nzt,ifrst,& + & aa%val,aa%irp,aa%ja,p%iprcparm(mld_slud_ptr_),& + & npr, npc, info) + if (info /= psb_success_) then + ch_err='psb_sludist_fact' + call psb_errpush(4110,name,a_err=ch_err,i_err=(/info,0,0,0,0/)) + goto 9999 + end if + + call psb_glob_to_loc(aa%ja(1:nzt),desc_a,info,iact='I') + + class default info=psb_err_unsupported_format_ - call psb_errpush(info,name,a_err=a%fida) - goto 9999 - endif - - ! - ! WARN: we need to check for a BLOCK distribution (this is the - ! distribution required by SuperLU_DIST) - ! - nrow = psb_cd_get_local_rows(desc_a) - ncol = psb_cd_get_local_cols(desc_a) - call psb_loc_to_glob(1,ifrst,desc_a,info) - call psb_loc_to_glob(nrow,ibcheck,desc_a,info) - ibcheck = ibcheck - ifrst + 1 - ibcheck = ibcheck - nrow - call psb_amx(ictxt,ibcheck) - if (ibcheck > 0) then - write(0,*) 'Warning: does not look like a BLOCK distribution' - endif - - mglob = psb_cd_get_global_rows(desc_a) - nzt = psb_sp_get_nnzeros(a) - - npr = np - npc = 1 - call psb_loc_to_glob(a%ia1(1:nzt),desc_a,info,iact='I') - - ! - ! Compute the LU factorization - ! - call mld_zsludist_fact(mglob,nrow,nzt,ifrst,& - & a%aspk,a%ia2,a%ia1,p%iprcparm(mld_slud_ptr_),& - & npr, npc, info) - if (info /= psb_success_) then - ch_err='psb_sludist_fact' - call psb_errpush(4110,name,a_err=ch_err,i_err=(/info,0,0,0,0/)) + ch_err = aa%get_fmt() + call psb_errpush(info,name,a_err=ch_err) goto 9999 - end if - - call psb_glob_to_loc(a%ia1(1:nzt),desc_a,info,iact='I') + end select + call psb_erractionrestore(err_act) return