From a8b63aea62bcfc92f9776b3c502f8ba4b048e9d8 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Tue, 2 Nov 2010 16:20:05 +0000 Subject: [PATCH] mld2p4: Makefile mld_dmlprec_aply.f90 mld_s_jac_smoother.f03 mld_smlprec_aply.f90 Second step of SINGLE PRECISION implementation. --- mlprec/Makefile | 3 +- mlprec/mld_dmlprec_aply.f90 | 3 - mlprec/mld_s_jac_smoother.f03 | 508 ++++++++++ mlprec/mld_smlprec_aply.f90 | 1684 ++++++++++++--------------------- 4 files changed, 1104 insertions(+), 1094 deletions(-) create mode 100644 mlprec/mld_s_jac_smoother.f03 diff --git a/mlprec/Makefile b/mlprec/Makefile index 07d1c647..77795d7c 100644 --- a/mlprec/Makefile +++ b/mlprec/Makefile @@ -10,7 +10,7 @@ MODOBJS=mld_base_prec_type.o \ mld_s_prec_type.o mld_d_prec_type.o mld_c_prec_type.o mld_z_prec_type.o \ mld_prec_type.o mld_prec_mod.o mld_inner_mod.o mld_move_alloc_mod.o\ mld_d_ilu_solver.o mld_d_diag_solver.o mld_d_jac_smoother.o mld_d_as_smoother.o \ - mld_s_ilu_solver.o mld_s_diag_solver.o mld_s_as_smoother.o + mld_s_ilu_solver.o mld_s_diag_solver.o mld_s_jac_smoother.o mld_s_as_smoother.o MPFOBJS=mld_daggrmat_nosmth_asb.o mld_daggrmat_smth_asb.o mld_daggrmat_minnrg_asb.o \ mld_saggrmat_nosmth_asb.o mld_saggrmat_smth_asb.o MPCOBJS=mld_sslud_interface.o mld_dslud_interface.o mld_cslud_interface.o mld_zslud_interface.o @@ -19,6 +19,7 @@ INNEROBJS= mld_dcoarse_bld.o mld_dmlprec_bld.o mld_dslu_bld.o mld_dumf_bld.o \ mld_dmlprec_aply.o mld_dslud_bld.o mld_daggrmat_asb.o \ mld_scoarse_bld.o mld_saggrmap_bld.o \ mld_silu0_fact.o mld_siluk_fact.o mld_silut_fact.o \ + mld_smlprec_aply.o \ $(MPFOBJS) # diff --git a/mlprec/mld_dmlprec_aply.f90 b/mlprec/mld_dmlprec_aply.f90 index 348ffae7..32f5a7cf 100644 --- a/mlprec/mld_dmlprec_aply.f90 +++ b/mlprec/mld_dmlprec_aply.f90 @@ -481,9 +481,6 @@ contains end if sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_) -!!$ 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,& diff --git a/mlprec/mld_s_jac_smoother.f03 b/mlprec/mld_s_jac_smoother.f03 new file mode 100644 index 00000000..5f528761 --- /dev/null +++ b/mlprec/mld_s_jac_smoother.f03 @@ -0,0 +1,508 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) +!!$ +!!$ (C) Copyright 2008,2009,2010, 2010 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +! +! +! +! +! +! +module mld_s_jac_smoother + + use mld_s_prec_type + + type, extends(mld_s_base_smoother_type) :: mld_s_jac_smoother_type + ! The local solver component is inherited from the + ! parent type. + ! class(mld_s_base_solver_type), allocatable :: sv + ! + type(psb_sspmat_type) :: nd + contains + procedure, pass(sm) :: build => s_jac_smoother_bld + procedure, pass(sm) :: apply => s_jac_smoother_apply + procedure, pass(sm) :: free => s_jac_smoother_free + procedure, pass(sm) :: seti => s_jac_smoother_seti + procedure, pass(sm) :: setc => s_jac_smoother_setc + procedure, pass(sm) :: setr => s_jac_smoother_setr + procedure, pass(sm) :: descr => s_jac_smoother_descr + procedure, pass(sm) :: sizeof => s_jac_smoother_sizeof + end type mld_s_jac_smoother_type + + + private :: s_jac_smoother_bld, s_jac_smoother_apply, & + & s_jac_smoother_free, s_jac_smoother_seti, & + & s_jac_smoother_setc, s_jac_smoother_setr,& + & s_jac_smoother_descr, s_jac_smoother_sizeof + + + +contains + + subroutine s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) + use psb_sparse_mod + type(psb_desc_type), intent(in) :: desc_data + class(mld_s_jac_smoother_type), intent(in) :: sm + real(psb_spk_),intent(in) :: x(:) + real(psb_spk_),intent(inout) :: y(:) + real(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + integer, intent(in) :: sweeps + real(psb_spk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='d_jac_smoother_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T','C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + if (.not.allocated(sm%sv)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + n_row = psb_cd_get_local_rows(desc_data) + n_col = psb_cd_get_local_cols(desc_data) + + if (n_col <= size(work)) then + ww => work(1:n_col) + if ((4*n_col+n_col) <= size(work)) then + aux => work(n_col+1:) + else + allocate(aux(4*n_col),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/4*n_col,0,0,0,0/),& + & a_err='real(psb_spk_)') + goto 9999 + end if + endif + else + allocate(ww(n_col),aux(4*n_col),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/5*n_col,0,0,0,0/),& + & a_err='real(psb_spk_)') + goto 9999 + end if + endif + + if (sweeps == 1) then + + call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,a_err='Error in sub_aply Jacobi Sweeps = 1') + goto 9999 + endif + + else if (sweeps > 1) then + + ! + ! + ! Apply multiple sweeps of a block-Jacobi solver + ! to compute an approximate solution of a linear system. + ! + ! + allocate(tx(n_col),ty(n_col),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/2*n_col,0,0,0,0/),& + & a_err='real(psb_spk_)') + goto 9999 + end if + + tx = szero + ty = szero + do i=1, sweeps + ! + ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the + ! block diagonal part and the remaining part of the local matrix + ! and Y(j) is the approximate solution at sweep j. + ! + ty(1:n_row) = x(1:n_row) + call psb_spmm(-sone,sm%nd,tx,sone,ty,desc_data,info,work=aux,trans=trans_) + + if (info /= psb_success_) exit + + call sm%sv%apply(sone,ty,szero,tx,desc_data,trans_,aux,info) + + if (info /= psb_success_) exit + end do + + if (info == psb_success_) call psb_geaxpby(alpha,tx,beta,y,desc_data,info) + + if (info /= psb_success_) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='subsolve with Jacobi sweeps > 1') + goto 9999 + end if + + deallocate(tx,ty,stat=info) + if (info /= psb_success_) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='final cleanup with Jacobi sweeps > 1') + goto 9999 + end if + + else + + info = psb_err_iarg_neg_ + call psb_errpush(info,name,& + & i_err=(/2,sweeps,0,0,0/)) + goto 9999 + + endif + + + if (n_col <= size(work)) then + if ((4*n_col+n_col) <= size(work)) then + else + deallocate(aux) + endif + else + deallocate(ww,aux) + endif + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine s_jac_smoother_apply + + subroutine s_jac_smoother_bld(a,desc_a,sm,upd,info) + + use psb_sparse_mod + use mld_s_diag_solver + Implicit None + + ! Arguments + type(psb_sspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_s_jac_smoother_type), intent(inout) :: sm + character, intent(in) :: upd + integer, intent(out) :: info + ! Local variables + integer :: n_row,n_col, nrow_a, nztota + real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me,i, err_act, debug_unit, debug_level + character(len=20) :: name='d_jac_smoother_bld', ch_err + + 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_a) + call psb_info(ictxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + + n_row = psb_cd_get_local_rows(desc_a) + + nrow_a = a%get_nrows() + nztota = a%get_nzeros() + select type (smsv => sm%sv) + type is (mld_s_diag_solver_type) + call a%clip_diag(sm%nd,info) + class default + call a%csclip(sm%nd,info,& + & jmin=nrow_a+1,rscale=.false.,cscale=.false.) + end select + if (info == psb_success_) call sm%nd%cscnv(info,& + & type='csr',dupl=psb_dupl_add_) + if (info == psb_success_) & + & call sm%sv%build(a,desc_a,upd,info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='clip & psb_spcnv csr 4') + goto 9999 + end if + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine s_jac_smoother_bld + + + subroutine s_jac_smoother_seti(sm,what,val,info) + + use psb_sparse_mod + + Implicit None + + ! Arguments + class(mld_s_jac_smoother_type), intent(inout) :: sm + integer, intent(in) :: what + integer, intent(in) :: val + integer, intent(out) :: info + Integer :: err_act + character(len=20) :: name='d_jac_smoother_seti' + + info = psb_success_ + call psb_erractionsave(err_act) + + select case(what) +!!$ case(mld_smoother_sweeps_) +!!$ sm%sweeps = val + case default + if (allocated(sm%sv)) then + call sm%sv%set(what,val,info) + else +!!$ write(0,*) trim(name),' Missing component, not setting!' +!!$ info = 1121 + end if + end select + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine s_jac_smoother_seti + + subroutine s_jac_smoother_setc(sm,what,val,info) + + use psb_sparse_mod + + Implicit None + + ! Arguments + class(mld_s_jac_smoother_type), intent(inout) :: sm + integer, intent(in) :: what + character(len=*), intent(in) :: val + integer, intent(out) :: info + Integer :: err_act, ival + character(len=20) :: name='d_jac_smoother_setc' + + info = psb_success_ + call psb_erractionsave(err_act) + + + call mld_stringval(val,ival,info) + if (info == psb_success_) call sm%set(what,ival,info) + if (info /= psb_success_) then + info = psb_err_from_subroutine_ + call psb_errpush(info, name) + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine s_jac_smoother_setc + + subroutine s_jac_smoother_setr(sm,what,val,info) + + use psb_sparse_mod + + Implicit None + + ! Arguments + class(mld_s_jac_smoother_type), intent(inout) :: sm + integer, intent(in) :: what + real(psb_spk_), intent(in) :: val + integer, intent(out) :: info + Integer :: err_act + character(len=20) :: name='d_jac_smoother_setr' + + call psb_erractionsave(err_act) + info = psb_success_ + + + if (allocated(sm%sv)) then + call sm%sv%set(what,val,info) + else +!!$ write(0,*) trim(name),' Missing component, not setting!' +!!$ info = 1121 + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine s_jac_smoother_setr + + subroutine s_jac_smoother_free(sm,info) + + use psb_sparse_mod + + Implicit None + + ! Arguments + class(mld_s_jac_smoother_type), intent(inout) :: sm + integer, intent(out) :: info + Integer :: err_act + character(len=20) :: name='d_jac_smoother_free' + + call psb_erractionsave(err_act) + info = psb_success_ + + + + if (allocated(sm%sv)) then + call sm%sv%free(info) + if (info == psb_success_) deallocate(sm%sv,stat=info) + if (info /= psb_success_) then + info = psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + end if + end if + call sm%nd%free() + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine s_jac_smoother_free + + subroutine s_jac_smoother_descr(sm,info,iout) + + use psb_sparse_mod + + Implicit None + + ! Arguments + class(mld_s_jac_smoother_type), intent(in) :: sm + integer, intent(out) :: info + integer, intent(in), optional :: iout + + ! Local variables + integer :: err_act + integer :: ictxt, me, np + character(len=20), parameter :: name='mld_s_jac_smoother_descr' + integer :: iout_ + + call psb_erractionsave(err_act) + info = psb_success_ + if (present(iout)) then + iout_ = iout + else + iout_ = 6 + endif + + write(iout_,*) ' Block Jacobi smoother ' + write(iout_,*) ' Local solver:' + if (allocated(sm%sv)) then + call sm%sv%descr(info,iout_) + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine s_jac_smoother_descr + + function s_jac_smoother_sizeof(sm) result(val) + use psb_sparse_mod + implicit none + ! Arguments + class(mld_s_jac_smoother_type), intent(in) :: sm + integer(psb_long_int_k_) :: val + integer :: i + + val = psb_sizeof_int + if (allocated(sm%sv)) val = val + sm%sv%sizeof() + val = val + sm%nd%sizeof() + + return + end function s_jac_smoother_sizeof + +end module mld_s_jac_smoother diff --git a/mlprec/mld_smlprec_aply.f90 b/mlprec/mld_smlprec_aply.f90 index 2d0f4c2a..35de2390 100644 --- a/mlprec/mld_smlprec_aply.f90 +++ b/mlprec/mld_smlprec_aply.f90 @@ -58,7 +58,8 @@ ! A 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 +! 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. @@ -137,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_smlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info) use psb_sparse_mod @@ -157,9 +331,13 @@ subroutine mld_smlprec_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 + real(psb_spk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) + end type psb_mlprec_wrk_type + type(psb_mlprec_wrk_type), allocatable :: mlprec_wrk(:) name='mld_smlprec_aply' info = psb_success_ @@ -176,78 +354,43 @@ subroutine mld_smlprec_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_spk_)') + 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(:) = szero - 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 @@ -261,156 +404,45 @@ subroutine mld_smlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info) return contains - ! - ! Subroutine: add_ml_aply - ! Version: real - ! Note: internal subroutine of mld_smlprec_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_smlprec_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_sprec_type), intent(in) :: p - real(psb_spk_),intent(in) :: alpha,beta - real(psb_spk_),intent(in) :: x(:) - real(psb_spk_),intent(inout) :: y(:) - character, intent(in) :: trans - real(psb_spk_),target :: work(:) - integer, intent(out) :: info + integer :: level + type(mld_sprec_type), intent(in) :: p + type(psb_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:) + character, intent(in) :: trans + real(psb_spk_),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 - real(psb_spk_), 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='real(psb_spk_)') + 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(:) = szero - 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_ @@ -418,921 +450,392 @@ contains & a_err='real(psb_spk_)') goto 9999 end if - - ! Apply the prolongator transpose, i.e. the restriction - call psb_map_X2Y(sone,mlprec_wrk(ilev-1)%x2l,& - & szero,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(sone,p%precv(ilev)%prec,& - & mlprec_wrk(ilev)%x2l,szero,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(sone,mlprec_wrk(ilev)%y2l,& - & sone,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,sone,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: real - ! Note: internal subroutine of mld_smlprec_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_smlprec_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_sprec_type), intent(in) :: p - real(psb_spk_),intent(in) :: alpha,beta - real(psb_spk_),intent(in) :: x(:) - real(psb_spk_),intent(inout) :: y(:) - character, intent(in) :: trans - real(psb_spk_),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_spk_), 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='real(psb_spk_)') - goto 9999 - end if - - mlprec_wrk(1)%x2l(:) = x - ! - ! STEP 2 - ! - ! Apply the base preconditioner at the finest level - ! - call mld_baseprec_aply(sone,p%precv(1)%prec,mlprec_wrk(1)%x2l,& - & szero,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(-sone,p%precv(1)%base_a,mlprec_wrk(1)%y2l,& - & sone,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='real(psb_spk_)') - goto 9999 - end if - - ! Apply the prolongator transpose, i.e. the restriction - call psb_map_X2Y(sone,mlprec_wrk(ilev-1)%tx,& - & szero,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(sone,p%precv(ilev)%prec,mlprec_wrk(ilev)%x2l,& - & szero,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(-sone,p%precv(ilev)%base_a,& - & mlprec_wrk(ilev)%y2l,sone,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(sone,mlprec_wrk(ilev+1)%y2l,& - & sone,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: real - ! Note: internal subroutine of mld_smlprec_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_smlprec_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_sprec_type), intent(in) :: p - real(psb_spk_),intent(in) :: alpha,beta - real(psb_spk_),intent(in) :: x(:) - real(psb_spk_),intent(inout) :: y(:) - character, intent(in) :: trans - real(psb_spk_),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_spk_), 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='real(psb_spk_)') - goto 9999 - end if - - call psb_geaxpby(sone,x,szero,mlprec_wrk(1)%tx,& - & p%precv(1)%base_desc,info) - call psb_geaxpby(sone,x,szero,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='real(psb_spk_)') - goto 9999 - end if - - ! Apply the prolongator transpose, i.e. the restriction - call psb_map_X2Y(sone,mlprec_wrk(ilev-1)%x2l,& - & szero,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(sone,mlprec_wrk(ilev)%x2l,szero,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(sone,p%precv(nlev)%prec,mlprec_wrk(nlev)%x2l, & - & szero, 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(sone,mlprec_wrk(level-1)%x2l,& + & szero,mlprec_wrk(level)%x2l,& + & p%precv(level)%map,info,work=work) - ! - ! Apply the prolongator - ! - call psb_map_Y2X(sone,mlprec_wrk(ilev+1)%y2l,& - & szero,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(-sone,p%precv(ilev)%base_a,mlprec_wrk(ilev)%y2l,& - & sone,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(sone,p%precv(ilev)%prec,& - & mlprec_wrk(ilev)%tx,sone,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(sone,& + & mlprec_wrk(level)%x2l,szero,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(sone,mlprec_wrk(level+1)%y2l,& + & sone,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: real - ! Note: internal subroutine of mld_smlprec_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_smlprec_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_sprec_type), intent(in) :: p - real(psb_spk_),intent(in) :: alpha,beta - real(psb_spk_),intent(in) :: x(:) - real(psb_spk_),intent(inout) :: y(:) - character, intent(in) :: trans - real(psb_spk_),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_spk_), 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='real(psb_spk_)') - goto 9999 - end if - - call psb_geaxpby(sone,x,szero,mlprec_wrk(1)%x2l,& - & p%precv(1)%base_desc,info) - call psb_geaxpby(sone,x,szero,mlprec_wrk(1)%tx,& - & p%precv(1)%base_desc,info) - - ! - ! STEP 2 - ! - ! Apply the base preconditioner at the finest level - ! - call mld_baseprec_aply(sone,p%precv(1)%prec,mlprec_wrk(1)%x2l,& - & szero,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(-sone,p%precv(1)%base_a,mlprec_wrk(1)%y2l,& - & sone,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='real(psb_spk_)') - goto 9999 - end if - - ! Apply the prolongator transpose, i.e. the restriction - call psb_map_X2Y(sone,mlprec_wrk(ilev-1)%ty,& - & szero,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(sone,mlprec_wrk(ilev)%x2l,szero,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(sone,p%precv(ilev)%prec,& - & mlprec_wrk(ilev)%x2l,szero,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(-sone,p%precv(ilev)%base_a,& - & mlprec_wrk(ilev)%y2l,sone,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(sone,mlprec_wrk(level-1)%x2l,& + & szero,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(sone,mlprec_wrk(level+1)%y2l,& + & szero,mlprec_wrk(level)%y2l,& + & p%precv(level+1)%map,info,work=work) + if (info /= psb_success_) goto 9999 + ! + ! Compute the residual + ! + call psb_spmm(-sone,p%precv(level)%base_a,mlprec_wrk(level)%y2l,& + & sone,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(sone,& + & mlprec_wrk(level)%x2l,sone,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(sone,& + & mlprec_wrk(level)%x2l,szero,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(sone,mlprec_wrk(level-1)%x2l,& + & szero,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(sone,& + & mlprec_wrk(level)%x2l,szero,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(-sone,p%precv(level)%base_a,& + & mlprec_wrk(level)%y2l,sone,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(sone,mlprec_wrk(level+1)%y2l,& + & sone,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(sone,mlprec_wrk(level-1)%x2l,& + & szero,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(sone,& + & mlprec_wrk(level)%x2l,szero,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(-sone,p%precv(level)%base_a,& + & mlprec_wrk(level)%y2l,sone,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(sone,mlprec_wrk(level+1)%y2l,& + & sone,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(sone,mlprec_wrk(level-1)%x2l,& + & szero,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(sone,mlprec_wrk(level+1)%y2l,& + & szero,mlprec_wrk(level)%y2l,& + & p%precv(level+1)%map,info,work=work) + if (info /= psb_success_) goto 9999 + ! + ! Compute the residual + ! + call psb_spmm(-sone,p%precv(level)%base_a,mlprec_wrk(level)%y2l,& + & sone,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(sone,& + & mlprec_wrk(level)%x2l,sone,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(sone,& + & mlprec_wrk(level)%x2l,szero,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_spk_)') + goto 9999 + end if + + if (level > 1) then + ! Apply the restriction + call psb_map_X2Y(sone,mlprec_wrk(level-1)%ty,& + & szero,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(sone,mlprec_wrk(level)%x2l,szero,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(sone,& + & mlprec_wrk(level)%x2l,szero,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(-sone,p%precv(level)%base_a,& + & mlprec_wrk(level)%y2l,sone,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(sone,mlprec_wrk(level+1)%y2l,& + & sone,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(-sone,p%precv(level)%base_a,mlprec_wrk(level)%y2l,& + & sone,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(sone,& + & mlprec_wrk(level)%tx,sone,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(sone,mlprec_wrk(ilev+1)%y2l,& - & sone,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(-sone,p%precv(ilev)%base_a,mlprec_wrk(ilev)%y2l,& - & sone,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(sone,p%precv(ilev)%prec,mlprec_wrk(ilev)%tx,& - & sone,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 @@ -1344,7 +847,8 @@ contains return end if return - end subroutine mlt_twoside_ml_aply + + end subroutine inner_ml_aply end subroutine mld_smlprec_aply