diff --git a/mlprec/impl/solver/Makefile b/mlprec/impl/solver/Makefile index 047aff65..4fb01f25 100644 --- a/mlprec/impl/solver/Makefile +++ b/mlprec/impl/solver/Makefile @@ -66,6 +66,9 @@ mld_d_diag_solver_apply_vect.o \ mld_d_diag_solver_bld.o \ mld_d_diag_solver_clone.o \ mld_d_diag_solver_cnv.o \ +mld_d_bwgs_solver_bld.o \ +mld_d_bwgs_solver_apply.o \ +mld_d_bwgs_solver_apply_vect.o \ mld_d_gs_solver_bld.o \ mld_d_gs_solver_clone.o \ mld_d_gs_solver_cnv.o \ diff --git a/mlprec/impl/solver/mld_d_bwgs_solver_apply.f90 b/mlprec/impl/solver/mld_d_bwgs_solver_apply.f90 new file mode 100644 index 00000000..075e4f52 --- /dev/null +++ b/mlprec/impl/solver/mld_d_bwgs_solver_apply.f90 @@ -0,0 +1,196 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ 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. +!!$ +!!$ +subroutine mld_d_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + + use psb_base_mod + use mld_d_gs_solver, mld_protect_name => mld_d_bwgs_solver_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_bwgs_solver_type), intent(inout) :: sv + real(psb_dpk_),intent(inout) :: x(:) + real(psb_dpk_),intent(inout) :: y(:) + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + + integer(psb_ipk_) :: n_row,n_col, itx + real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) + real(psb_dpk_), allocatable :: temp(:),wv(:),xit(:) + integer(psb_ipk_) :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='d_bwgs_solver_apply' + + call psb_erractionsave(err_act) + ictxt = desc_data%get_ctxt() + call psb_info(ictxt,me,np) + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') +!!$ case('T') +!!$ case('C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + n_row = desc_data%get_local_rows() + n_col = desc_data%get_local_cols() + + 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,izero,izero,izero,izero/),& + & a_err='real(psb_dpk_)') + 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,izero,izero,izero,izero/),& + & a_err='real(psb_dpk_)') + goto 9999 + end if + endif + + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,& + & i_err=(/5*n_col,izero,izero,izero,izero/),& + & a_err='real(psb_dpk_)') + goto 9999 + end if + + call psb_geasb(wv,desc_data,info) + call psb_geasb(xit,desc_data,info) + + select case(trans_) + case('N') + if (sv%eps <=dzero) then + ! + ! Fixed number of iterations + ! + ! + ! WARNING: this is not completely satisfactory. We are assuming here Y + ! as the initial guess, but this is only working if we are called from the + ! current JAC smoother loop. A good solution would be to have a separate + ! input argument as the initial guess + ! +!!$ write(0,*) 'GS Iteration with ',sv%sweeps + call psb_geaxpby(done,y,dzero,xit,desc_data,info) + do itx=1,sv%sweeps + call psb_geaxpby(done,x,dzero,wv,desc_data,info) + ! Update with U. The off-diagonal block is taken care + ! from the Jacobi smoother, hence this is purely local. + call psb_spmm(-done,sv%u,xit,done,wv,desc_data,info,doswap=.false.) + call psb_spsm(done,sv%l,wv,dzero,xit,desc_data,info) +!!$ temp = xit%get_vect() +!!$ write(0,*) me,'GS Iteration ',itx,':',temp(1:n_row) + end do + + call psb_geaxpby(alpha,xit,beta,y,desc_data,info) + + else + ! + ! Iterations to convergence, not implemented right now. + ! + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='EPS>0 not implemented in GS subsolve') + goto 9999 + + end if +!!$ case('T') +!!$ call psb_spsm(done,sv%u,x,dzero,wv,desc_data,info,& +!!$ & trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux) +!!$ if (info == psb_success_) call psb_spsm(alpha,sv%l,wv,beta,y,desc_data,info,& +!!$ & trans=trans_,scale='U',choice=psb_none_,work=aux) +!!$ +!!$ case('C') +!!$ +!!$ call psb_spsm(done,sv%u,x,dzero,wv,desc_data,info,& +!!$ & trans=trans_,scale='U',choice=psb_none_,work=aux) +!!$ +!!$ call wv1%mlt(done,sv%dv,wv,dzero,info,conjgx=trans_) +!!$ +!!$ if (info == psb_success_) call psb_spsm(alpha,sv%l,wv1,beta,y,desc_data,info,& +!!$ & trans=trans_,scale='U',choice=psb_none_,work=aux) + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid TRANS in GS subsolve') + goto 9999 + end select + + + if (info /= psb_success_) then + + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in subsolve') + 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 call psb_error_handler(err_act) + + return + +end subroutine mld_d_bwgs_solver_apply diff --git a/mlprec/impl/solver/mld_d_bwgs_solver_apply_vect.f90 b/mlprec/impl/solver/mld_d_bwgs_solver_apply_vect.f90 new file mode 100644 index 00000000..f0ab056f --- /dev/null +++ b/mlprec/impl/solver/mld_d_bwgs_solver_apply_vect.f90 @@ -0,0 +1,200 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ 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. +!!$ +!!$ +subroutine mld_d_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + + use psb_base_mod + use mld_d_gs_solver, mld_protect_name => mld_d_bwgs_solver_apply_vect + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_bwgs_solver_type), intent(inout) :: sv + type(psb_d_vect_type),intent(inout) :: x + type(psb_d_vect_type),intent(inout) :: y + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + + integer(psb_ipk_) :: n_row,n_col, itx + type(psb_d_vect_type) :: wv, xit + real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) + real(psb_dpk_), allocatable :: temp(:) + integer(psb_ipk_) :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='d_bwgs_solver_apply' + + call psb_erractionsave(err_act) + ictxt = desc_data%get_ctxt() + call psb_info(ictxt,me,np) + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') +!!$ case('T') +!!$ case('C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + n_row = desc_data%get_local_rows() + n_col = desc_data%get_local_cols() + + + if (x%get_nrows() < n_row) then + info = 36 + call psb_errpush(info,name,& + & i_err=(/itwo,n_row,izero,izero,izero/)) + goto 9999 + end if + if (y%get_nrows() < n_row) then + info = 36 + call psb_errpush(info,name,& + & i_err=(/ithree,n_row,izero,izero,izero/)) + goto 9999 + end if + + + + 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) + endif + else + allocate(ww(n_col),aux(4*n_col),stat=info) + endif + + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,& + & i_err=(/5*n_col,izero,izero,izero,izero/),& + & a_err='real(psb_dpk_)') + goto 9999 + end if + + call psb_geasb(wv,desc_data,info,mold=x%v,scratch=.true.) + call psb_geasb(xit,desc_data,info,mold=x%v,scratch=.true.) + + select case(trans_) + case('N') + if (sv%eps <=dzero) then + ! + ! Fixed number of iterations + ! + ! + ! WARNING: this is not completely satisfactory. We are assuming here Y + ! as the initial guess, but this is only working if we are called from the + ! current JAC smoother loop. A good solution would be to have a separate + ! input argument as the initial guess + ! +!!$ write(0,*) 'GS Iteration with ',sv%sweeps + call psb_geaxpby(done,y,dzero,xit,desc_data,info) + do itx=1,sv%sweeps + call psb_geaxpby(done,x,dzero,wv,desc_data,info) + ! Update with U. The off-diagonal block is taken care + ! from the Jacobi smoother, hence this is purely local. + call psb_spmm(-done,sv%u,xit,done,wv,desc_data,info,doswap=.false.) + call psb_spsm(done,sv%l,wv,dzero,xit,desc_data,info) +!!$ temp = xit%get_vect() +!!$ write(0,*) me,'GS Iteration ',itx,':',temp(1:n_row) + end do + + call psb_geaxpby(alpha,xit,beta,y,desc_data,info) + + else + ! + ! Iterations to convergence, not implemented right now. + ! + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='EPS>0 not implemented in GS subsolve') + goto 9999 + + end if +!!$ case('T') +!!$ call psb_spsm(done,sv%u,x,dzero,wv,desc_data,info,& +!!$ & trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux) +!!$ if (info == psb_success_) call psb_spsm(alpha,sv%l,wv,beta,y,desc_data,info,& +!!$ & trans=trans_,scale='U',choice=psb_none_,work=aux) +!!$ +!!$ case('C') +!!$ +!!$ call psb_spsm(done,sv%u,x,dzero,wv,desc_data,info,& +!!$ & trans=trans_,scale='U',choice=psb_none_,work=aux) +!!$ +!!$ call wv1%mlt(done,sv%dv,wv,dzero,info,conjgx=trans_) +!!$ +!!$ if (info == psb_success_) call psb_spsm(alpha,sv%l,wv1,beta,y,desc_data,info,& +!!$ & trans=trans_,scale='U',choice=psb_none_,work=aux) + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid TRANS in GS subsolve') + goto 9999 + end select + + + if (info /= psb_success_) then + + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in subsolve') + goto 9999 + endif + call wv%free(info) + call xit%free(info) + 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 call psb_error_handler(err_act) + + return + +end subroutine mld_d_bwgs_solver_apply_vect diff --git a/mlprec/impl/solver/mld_d_bwgs_solver_bld.f90 b/mlprec/impl/solver/mld_d_bwgs_solver_bld.f90 new file mode 100644 index 00000000..7547cd38 --- /dev/null +++ b/mlprec/impl/solver/mld_d_bwgs_solver_bld.f90 @@ -0,0 +1,117 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ 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. +!!$ +!!$ +subroutine mld_d_bwgs_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold,imold) + + use psb_base_mod + use mld_d_gs_solver, mld_protect_name => mld_d_bwgs_solver_bld + + Implicit None + + ! Arguments + type(psb_dspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_d_bwgs_solver_type), intent(inout) :: sv + character, intent(in) :: upd + integer(psb_ipk_), intent(out) :: info + type(psb_dspmat_type), intent(in), target, optional :: b + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + ! Local variables + integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota + integer(psb_ipk_) :: ictxt,np,me,i, err_act, debug_unit, debug_level + character(len=20) :: name='d_bwgs_solver_bld', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = desc_a%get_context() + call psb_info(ictxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + + n_row = desc_a%get_local_rows() + + if (psb_toupper(upd) == 'F') then + nrow_a = a%get_nrows() + nztota = a%get_nzeros() +!!$ if (present(b)) then +!!$ nztota = nztota + b%get_nzeros() +!!$ end if + if (sv%eps <= dzero) then + ! + ! This cuts out the off-diagonal part, because it's supposed to + ! be handled by the outer Jacobi smoother. + ! + call a%tril(sv%l,info) + call a%triu(sv%u,info,diag=1,jmax=nrow_a) + + else + + info = psb_err_missing_override_method_ + call psb_errpush(info,name) + goto 9999 + end if + + + + call sv%l%set_asb() + call sv%l%trim() + call sv%u%set_asb() + call sv%u%trim() + + if (present(amold)) then + call sv%l%cscnv(info,mold=amold) + call sv%u%cscnv(info,mold=amold) + end if + + end if + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine mld_d_bwgs_solver_bld diff --git a/mlprec/mld_d_gs_solver.f90 b/mlprec/mld_d_gs_solver.f90 index 230c5d89..feed37c5 100644 --- a/mlprec/mld_d_gs_solver.f90 +++ b/mlprec/mld_d_gs_solver.f90 @@ -76,9 +76,9 @@ module mld_d_gs_solver type, extends(mld_d_gs_solver_type) :: mld_d_bwgs_solver_type contains -!!$ procedure, pass(sv) :: build => mld_d_bwgs_solver_bld -!!$ procedure, pass(sv) :: apply_v => mld_d_bwgs_solver_apply_vect -!!$ procedure, pass(sv) :: apply_a => mld_d_bwgs_solver_apply + procedure, pass(sv) :: build => mld_d_bwgs_solver_bld + procedure, pass(sv) :: apply_v => mld_d_bwgs_solver_apply_vect + procedure, pass(sv) :: apply_a => mld_d_bwgs_solver_apply procedure, nopass :: get_fmt => d_bwgs_solver_get_fmt procedure, pass(sv) :: descr => d_bwgs_solver_descr end type mld_d_bwgs_solver_type @@ -110,6 +110,19 @@ module mld_d_gs_solver real(psb_dpk_),target, intent(inout) :: work(:) integer(psb_ipk_), intent(out) :: info end subroutine mld_d_gs_solver_apply_vect + subroutine mld_d_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, mld_d_bwgs_solver_type, psb_d_vect_type, psb_dpk_, & + & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type, psb_ipk_ + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_bwgs_solver_type), intent(inout) :: sv + type(psb_d_vect_type),intent(inout) :: x + type(psb_d_vect_type),intent(inout) :: y + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + end subroutine mld_d_bwgs_solver_apply_vect end interface interface @@ -126,6 +139,19 @@ module mld_d_gs_solver real(psb_dpk_),target, intent(inout) :: work(:) integer(psb_ipk_), intent(out) :: info end subroutine mld_d_gs_solver_apply + subroutine mld_d_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, mld_d_bwgs_solver_type, psb_d_vect_type, psb_dpk_, & + & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type, psb_ipk_ + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_bwgs_solver_type), intent(inout) :: sv + real(psb_dpk_),intent(inout) :: x(:) + real(psb_dpk_),intent(inout) :: y(:) + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + end subroutine mld_d_bwgs_solver_apply end interface interface @@ -144,6 +170,21 @@ module mld_d_gs_solver class(psb_d_base_vect_type), intent(in), optional :: vmold class(psb_i_base_vect_type), intent(in), optional :: imold end subroutine mld_d_gs_solver_bld + subroutine mld_d_bwgs_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold,imold) + import :: psb_desc_type, mld_d_bwgs_solver_type, psb_d_vect_type, psb_dpk_, & + & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,& + & psb_ipk_, psb_i_base_vect_type + implicit none + type(psb_dspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_d_bwgs_solver_type), intent(inout) :: sv + character, intent(in) :: upd + integer(psb_ipk_), intent(out) :: info + type(psb_dspmat_type), intent(in), target, optional :: b + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + end subroutine mld_d_bwgs_solver_bld end interface interface diff --git a/tests/pdegen/ppde3d-gs.f90 b/tests/pdegen/ppde3d-gs.f90 index 8cf6a9d4..1374c4b3 100644 --- a/tests/pdegen/ppde3d-gs.f90 +++ b/tests/pdegen/ppde3d-gs.f90 @@ -171,6 +171,7 @@ program ppde3d integer(psb_ipk_) :: nlev ! Number of levels in multilevel prec. character(len=16) :: aggrkind ! smoothed/raw aggregatin character(len=16) :: aggr_alg ! local or global aggregation + character(len=16) :: aggr_ord ! Ordering for aggregation character(len=16) :: mltype ! additive or multiplicative 2nd level prec character(len=16) :: smthpos ! side: pre, post, both smoothing integer(psb_ipk_) :: csize ! aggregation size at which to stop. @@ -255,6 +256,7 @@ program ppde3d call mld_precset(prec,'sub_iluthrs', prectype%thr1, info) call mld_precset(prec,'aggr_kind', prectype%aggrkind,info) call mld_precset(prec,'aggr_alg', prectype%aggr_alg,info) + call mld_precset(prec,'aggr_ord', prectype%aggr_ord,info) call mld_precset(prec,'ml_type', prectype%mltype, info) call mld_precset(prec,'smoother_pos', prectype%smthpos, info) if (prectype%athres >= dzero) & @@ -400,6 +402,7 @@ contains call read_data(prectype%nlev,psb_inp_unit) ! Number of levels in multilevel prec. call read_data(prectype%aggrkind,psb_inp_unit) ! smoothed/raw aggregatin call read_data(prectype%aggr_alg,psb_inp_unit) ! local or global aggregation + call read_data(prectype%aggr_ord,psb_inp_unit) ! aggregation ordering call read_data(prectype%mltype,psb_inp_unit) ! additive or multiplicative 2nd level prec call read_data(prectype%smthpos,psb_inp_unit) ! side: pre, post, both smoothing call read_data(prectype%cmat,psb_inp_unit) ! coarse mat @@ -439,6 +442,7 @@ contains call psb_bcast(ictxt,prectype%nlev) ! Number of levels in multilevel prec. call psb_bcast(ictxt,prectype%aggrkind) ! smoothed/raw aggregatin call psb_bcast(ictxt,prectype%aggr_alg) ! local or global aggregation + call psb_bcast(ictxt,prectype%aggr_ord) ! aggregation ordering call psb_bcast(ictxt,prectype%mltype) ! additive or multiplicative 2nd level prec call psb_bcast(ictxt,prectype%smthpos) ! side: pre, post, both smoothing call psb_bcast(ictxt,prectype%cmat) ! coarse mat