diff --git a/mlprec/impl/Makefile b/mlprec/impl/Makefile index 49fce95a..87465261 100644 --- a/mlprec/impl/Makefile +++ b/mlprec/impl/Makefile @@ -25,25 +25,25 @@ DINNEROBJS= mld_dcoarse_bld.o mld_dmlprec_bld.o \ mld_dilu0_fact.o mld_diluk_fact.o mld_dilut_fact.o mld_daggrmap_bld.o \ mld_d_dec_map_bld.o mld_dmlprec_aply.o mld_daggrmat_asb.o \ $(DMPFOBJS) mld_d_base_solver_impl.o mld_d_base_smoother_impl.o mld_d_onelev_impl.o\ - mld_d_as_smoother_impl.o + mld_d_as_smoother_impl.o mld_d_jac_smoother_impl.o SINNEROBJS= mld_scoarse_bld.o mld_smlprec_bld.o \ mld_silu0_fact.o mld_siluk_fact.o mld_silut_fact.o mld_saggrmap_bld.o \ mld_s_dec_map_bld.o mld_smlprec_aply.o mld_saggrmat_asb.o \ $(SMPFOBJS) mld_s_base_solver_impl.o mld_s_base_smoother_impl.o mld_s_onelev_impl.o\ - mld_s_as_smoother_impl.o + mld_s_as_smoother_impl.o mld_s_jac_smoother_impl.o ZINNEROBJS= mld_zcoarse_bld.o mld_zmlprec_bld.o \ mld_zilu0_fact.o mld_ziluk_fact.o mld_zilut_fact.o mld_zaggrmap_bld.o \ mld_z_dec_map_bld.o mld_zmlprec_aply.o mld_zaggrmat_asb.o \ $(ZMPFOBJS) mld_z_base_solver_impl.o mld_z_base_smoother_impl.o mld_z_onelev_impl.o\ - mld_z_as_smoother_impl.o + mld_z_as_smoother_impl.o mld_z_jac_smoother_impl.o CINNEROBJS= mld_ccoarse_bld.o mld_cmlprec_bld.o \ mld_cilu0_fact.o mld_ciluk_fact.o mld_cilut_fact.o mld_caggrmap_bld.o \ mld_c_dec_map_bld.o mld_cmlprec_aply.o mld_caggrmat_asb.o \ $(CMPFOBJS) mld_c_base_solver_impl.o mld_c_base_smoother_impl.o mld_c_onelev_impl.o\ - mld_c_as_smoother_impl.o + mld_c_as_smoother_impl.o mld_c_jac_smoother_impl.o INNEROBJS= $(SINNEROBJS) $(DINNEROBJS) $(CINNEROBJS) $(ZINNEROBJS) diff --git a/mlprec/impl/mld_c_jac_smoother_impl.f90 b/mlprec/impl/mld_c_jac_smoother_impl.f90 new file mode 100644 index 00000000..fee018a7 --- /dev/null +++ b/mlprec/impl/mld_c_jac_smoother_impl.f90 @@ -0,0 +1,433 @@ +!!$ +!!$ +!!$ 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. +!!$ +!!$ +! +! +! +! +subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) + use psb_base_mod + use mld_c_jac_smoother, mld_protect_name => mld_c_jac_smoother_apply_vect + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_c_jac_smoother_type), intent(inout) :: sm + type(psb_c_vect_type),intent(inout) :: x + type(psb_c_vect_type),intent(inout) :: y + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + integer, intent(in) :: sweeps + complex(psb_spk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + type(psb_c_vect_type) :: tx, ty + complex(psb_spk_), pointer :: ww(:), aux(:) + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='c_jac_smoother_apply' + + call psb_erractionsave(err_act) + + 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 + + if (.not.allocated(sm%sv)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + 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,0,0,0,0/),& + & a_err='complex(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='complex(psb_spk_)') + goto 9999 + end if + endif + + if ((sweeps == 1).or.(sm%nnz_nd_tot==0)) 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. + ! + ! + call tx%bld(x%get_nrows(),mold=x%v) + call tx%set(czero) + call ty%bld(x%get_nrows(),mold=x%v) + + 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. + ! + call psb_geaxpby(cone,x,czero,ty,desc_data,info) + call psb_spmm(-cone,sm%nd,tx,cone,ty,desc_data,info,work=aux,trans=trans_) + + if (info /= psb_success_) exit + + call sm%sv%apply(cone,ty,czero,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 + + call tx%free(info) + if (info == psb_success_) call ty%free(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 mld_c_jac_smoother_apply_vect + +subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) + use psb_base_mod + use mld_c_jac_smoother, mld_protect_name => mld_c_jac_smoother_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_c_jac_smoother_type), intent(in) :: sm + complex(psb_spk_),intent(inout) :: x(:) + complex(psb_spk_),intent(inout) :: y(:) + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + integer, intent(in) :: sweeps + complex(psb_spk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='c_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 = 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,0,0,0,0/),& + & a_err='complex(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='complex(psb_spk_)') + goto 9999 + end if + endif + + if ((sweeps == 1).or.(sm%nnz_nd_tot==0)) 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='complex(psb_spk_)') + goto 9999 + end if + + tx = czero + ty = czero + 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(-cone,sm%nd,tx,cone,ty,desc_data,info,work=aux,trans=trans_) + + if (info /= psb_success_) exit + + call sm%sv%apply(cone,ty,czero,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 mld_c_jac_smoother_apply + +subroutine mld_c_jac_smoother_bld(a,desc_a,sm,upd,info,amold,vmold) + + use psb_base_mod + use mld_c_diag_solver + use mld_c_jac_smoother, mld_protect_name => mld_c_jac_smoother_bld + Implicit None + + ! Arguments + type(psb_cspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_c_jac_smoother_type), intent(inout) :: sm + character, intent(in) :: upd + integer, intent(out) :: info + class(psb_c_base_sparse_mat), intent(in), optional :: amold + class(psb_c_base_vect_type), intent(in), optional :: vmold + ! Local variables + integer :: n_row,n_col, nrow_a, nztota, nzeros + complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me,i, err_act, debug_unit, debug_level + character(len=20) :: name='c_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 = 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() + + nrow_a = a%get_nrows() + nztota = a%get_nzeros() + select type (smsv => sm%sv) + type is (mld_c_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_) then + if (present(amold)) then + call sm%nd%cscnv(info,& + & mold=amold,dupl=psb_dupl_add_) + else + call sm%nd%cscnv(info,& + & type='csr',dupl=psb_dupl_add_) + endif + end if + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='clip & psb_spcnv csr 4') + goto 9999 + end if + + call sm%sv%build(a,desc_a,upd,info,amold=amold,vmold=vmold) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='solver build') + goto 9999 + end if + nzeros = sm%nd%get_nzeros() + call psb_sum(ictxt,nzeros) + sm%nnz_nd_tot = nzeros + 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 mld_c_jac_smoother_bld + diff --git a/mlprec/impl/mld_d_jac_smoother_impl.f90 b/mlprec/impl/mld_d_jac_smoother_impl.f90 new file mode 100644 index 00000000..ddff4c02 --- /dev/null +++ b/mlprec/impl/mld_d_jac_smoother_impl.f90 @@ -0,0 +1,433 @@ +!!$ +!!$ +!!$ 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. +!!$ +!!$ +! +! +! +! +subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) + use psb_base_mod + use mld_d_jac_smoother, mld_protect_name => mld_d_jac_smoother_apply_vect + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_jac_smoother_type), intent(inout) :: sm + 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 + integer, intent(in) :: sweeps + real(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + type(psb_d_vect_type) :: tx, ty + real(psb_dpk_), pointer :: ww(:), aux(:) + 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') + case('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 = 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,0,0,0,0/),& + & 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,0,0,0,0/),& + & a_err='real(psb_dpk_)') + goto 9999 + end if + endif + + if ((sweeps == 1).or.(sm%nnz_nd_tot==0)) 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. + ! + ! + call tx%bld(x%get_nrows(),mold=x%v) + call tx%set(dzero) + call ty%bld(x%get_nrows(),mold=x%v) + + 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. + ! + call psb_geaxpby(done,x,dzero,ty,desc_data,info) + call psb_spmm(-done,sm%nd,tx,done,ty,desc_data,info,work=aux,trans=trans_) + + if (info /= psb_success_) exit + + call sm%sv%apply(done,ty,dzero,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 + + call tx%free(info) + if (info == psb_success_) call ty%free(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 mld_d_jac_smoother_apply_vect + +subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) + use psb_base_mod + use mld_d_jac_smoother, mld_protect_name => mld_d_jac_smoother_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_jac_smoother_type), intent(in) :: sm + 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 + integer, intent(in) :: sweeps + real(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + real(psb_dpk_), 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 = 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,0,0,0,0/),& + & 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,0,0,0,0/),& + & a_err='real(psb_dpk_)') + goto 9999 + end if + endif + + if ((sweeps == 1).or.(sm%nnz_nd_tot==0)) 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_dpk_)') + goto 9999 + end if + + tx = dzero + ty = dzero + 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(-done,sm%nd,tx,done,ty,desc_data,info,work=aux,trans=trans_) + + if (info /= psb_success_) exit + + call sm%sv%apply(done,ty,dzero,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 mld_d_jac_smoother_apply + +subroutine mld_d_jac_smoother_bld(a,desc_a,sm,upd,info,amold,vmold) + + use psb_base_mod + use mld_d_diag_solver + use mld_d_jac_smoother, mld_protect_name => mld_d_jac_smoother_bld + Implicit None + + ! Arguments + type(psb_dspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_d_jac_smoother_type), intent(inout) :: sm + character, intent(in) :: upd + integer, intent(out) :: info + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold + ! Local variables + integer :: n_row,n_col, nrow_a, nztota, nzeros + real(psb_dpk_), 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 = 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() + + nrow_a = a%get_nrows() + nztota = a%get_nzeros() + select type (smsv => sm%sv) + type is (mld_d_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_) then + if (present(amold)) then + call sm%nd%cscnv(info,& + & mold=amold,dupl=psb_dupl_add_) + else + call sm%nd%cscnv(info,& + & type='csr',dupl=psb_dupl_add_) + endif + end if + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='clip & psb_spcnv csr 4') + goto 9999 + end if + + call sm%sv%build(a,desc_a,upd,info,amold=amold,vmold=vmold) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='solver build') + goto 9999 + end if + nzeros = sm%nd%get_nzeros() + call psb_sum(ictxt,nzeros) + sm%nnz_nd_tot = nzeros + 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 mld_d_jac_smoother_bld + diff --git a/mlprec/impl/mld_daggrmat_minnrg_asb.F90 b/mlprec/impl/mld_daggrmat_minnrg_asb.F90 index af068f9c..5a23e2a0 100644 --- a/mlprec/impl/mld_daggrmat_minnrg_asb.F90 +++ b/mlprec/impl/mld_daggrmat_minnrg_asb.F90 @@ -249,10 +249,14 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) call dap%clone(atmp,info) - call psb_sphalo(atmp,desc_a,am4,info,& + if (info == psb_success_) call psb_sphalo(atmp,desc_a,am4,info,& & colcnv=.false.,rowscale=.true.,outfmt='CSR ') if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=am4) if (info == psb_success_) call am4%free() + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='sphalo 2') + goto 9999 + end if call psb_symbmm(da,atmp,dadap,info) call psb_numbmm(da,atmp,dadap) diff --git a/mlprec/impl/mld_s_jac_smoother_impl.f90 b/mlprec/impl/mld_s_jac_smoother_impl.f90 new file mode 100644 index 00000000..c8115c6b --- /dev/null +++ b/mlprec/impl/mld_s_jac_smoother_impl.f90 @@ -0,0 +1,433 @@ +!!$ +!!$ +!!$ 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. +!!$ +!!$ +! +! +! +! +subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) + use psb_base_mod + use mld_s_jac_smoother, mld_protect_name => mld_s_jac_smoother_apply_vect + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_s_jac_smoother_type), intent(inout) :: sm + type(psb_s_vect_type),intent(inout) :: x + type(psb_s_vect_type),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 + type(psb_s_vect_type) :: tx, ty + real(psb_spk_), pointer :: ww(:), aux(:) + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='s_jac_smoother_apply' + + call psb_erractionsave(err_act) + + 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 + + if (.not.allocated(sm%sv)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + 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,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).or.(sm%nnz_nd_tot==0)) 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. + ! + ! + call tx%bld(x%get_nrows(),mold=x%v) + call tx%set(szero) + call ty%bld(x%get_nrows(),mold=x%v) + + 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. + ! + call psb_geaxpby(sone,x,szero,ty,desc_data,info) + 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 + + call tx%free(info) + if (info == psb_success_) call ty%free(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 mld_s_jac_smoother_apply_vect + +subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) + use psb_base_mod + use mld_s_jac_smoother, mld_protect_name => mld_s_jac_smoother_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_s_jac_smoother_type), intent(in) :: sm + real(psb_spk_),intent(inout) :: 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='s_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 = 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,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).or.(sm%nnz_nd_tot==0)) 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 mld_s_jac_smoother_apply + +subroutine mld_s_jac_smoother_bld(a,desc_a,sm,upd,info,amold,vmold) + + use psb_base_mod + use mld_s_diag_solver + use mld_s_jac_smoother, mld_protect_name => mld_s_jac_smoother_bld + 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 + class(psb_s_base_sparse_mat), intent(in), optional :: amold + class(psb_s_base_vect_type), intent(in), optional :: vmold + ! Local variables + integer :: n_row,n_col, nrow_a, nztota, nzeros + real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me,i, err_act, debug_unit, debug_level + character(len=20) :: name='s_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 = 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() + + 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_) then + if (present(amold)) then + call sm%nd%cscnv(info,& + & mold=amold,dupl=psb_dupl_add_) + else + call sm%nd%cscnv(info,& + & type='csr',dupl=psb_dupl_add_) + endif + end if + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='clip & psb_spcnv csr 4') + goto 9999 + end if + + call sm%sv%build(a,desc_a,upd,info,amold=amold,vmold=vmold) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='solver build') + goto 9999 + end if + nzeros = sm%nd%get_nzeros() + call psb_sum(ictxt,nzeros) + sm%nnz_nd_tot = nzeros + 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 mld_s_jac_smoother_bld + diff --git a/mlprec/impl/mld_z_jac_smoother_impl.f90 b/mlprec/impl/mld_z_jac_smoother_impl.f90 new file mode 100644 index 00000000..0dbd1ee1 --- /dev/null +++ b/mlprec/impl/mld_z_jac_smoother_impl.f90 @@ -0,0 +1,433 @@ +!!$ +!!$ +!!$ 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. +!!$ +!!$ +! +! +! +! +subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) + use psb_base_mod + use mld_z_jac_smoother, mld_protect_name => mld_z_jac_smoother_apply_vect + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_z_jac_smoother_type), intent(inout) :: sm + type(psb_z_vect_type),intent(inout) :: x + type(psb_z_vect_type),intent(inout) :: y + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + integer, intent(in) :: sweeps + complex(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + type(psb_z_vect_type) :: tx, ty + complex(psb_dpk_), pointer :: ww(:), aux(:) + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='z_jac_smoother_apply' + + call psb_erractionsave(err_act) + + 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 + + if (.not.allocated(sm%sv)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + 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,0,0,0,0/),& + & a_err='complex(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,0,0,0,0/),& + & a_err='complex(psb_dpk_)') + goto 9999 + end if + endif + + if ((sweeps == 1).or.(sm%nnz_nd_tot==0)) 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. + ! + ! + call tx%bld(x%get_nrows(),mold=x%v) + call tx%set(zzero) + call ty%bld(x%get_nrows(),mold=x%v) + + 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. + ! + call psb_geaxpby(zone,x,zzero,ty,desc_data,info) + call psb_spmm(-zone,sm%nd,tx,zone,ty,desc_data,info,work=aux,trans=trans_) + + if (info /= psb_success_) exit + + call sm%sv%apply(zone,ty,zzero,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 + + call tx%free(info) + if (info == psb_success_) call ty%free(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 mld_z_jac_smoother_apply_vect + +subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) + use psb_base_mod + use mld_z_jac_smoother, mld_protect_name => mld_z_jac_smoother_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_z_jac_smoother_type), intent(in) :: sm + complex(psb_dpk_),intent(inout) :: x(:) + complex(psb_dpk_),intent(inout) :: y(:) + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + integer, intent(in) :: sweeps + complex(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='z_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 = 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,0,0,0,0/),& + & a_err='complex(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,0,0,0,0/),& + & a_err='complex(psb_dpk_)') + goto 9999 + end if + endif + + if ((sweeps == 1).or.(sm%nnz_nd_tot==0)) 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='complex(psb_dpk_)') + goto 9999 + end if + + tx = zzero + ty = zzero + 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(-zone,sm%nd,tx,zone,ty,desc_data,info,work=aux,trans=trans_) + + if (info /= psb_success_) exit + + call sm%sv%apply(zone,ty,zzero,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 mld_z_jac_smoother_apply + +subroutine mld_z_jac_smoother_bld(a,desc_a,sm,upd,info,amold,vmold) + + use psb_base_mod + use mld_z_diag_solver + use mld_z_jac_smoother, mld_protect_name => mld_z_jac_smoother_bld + Implicit None + + ! Arguments + type(psb_zspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_z_jac_smoother_type), intent(inout) :: sm + character, intent(in) :: upd + integer, intent(out) :: info + class(psb_z_base_sparse_mat), intent(in), optional :: amold + class(psb_z_base_vect_type), intent(in), optional :: vmold + ! Local variables + integer :: n_row,n_col, nrow_a, nztota, nzeros + complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me,i, err_act, debug_unit, debug_level + character(len=20) :: name='z_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 = 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() + + nrow_a = a%get_nrows() + nztota = a%get_nzeros() + select type (smsv => sm%sv) + type is (mld_z_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_) then + if (present(amold)) then + call sm%nd%cscnv(info,& + & mold=amold,dupl=psb_dupl_add_) + else + call sm%nd%cscnv(info,& + & type='csr',dupl=psb_dupl_add_) + endif + end if + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='clip & psb_spcnv csr 4') + goto 9999 + end if + + call sm%sv%build(a,desc_a,upd,info,amold=amold,vmold=vmold) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='solver build') + goto 9999 + end if + nzeros = sm%nd%get_nzeros() + call psb_sum(ictxt,nzeros) + sm%nnz_nd_tot = nzeros + 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 mld_z_jac_smoother_bld + diff --git a/mlprec/mld_c_as_smoother.f90 b/mlprec/mld_c_as_smoother.f90 index fca84f42..bc1b07e3 100644 --- a/mlprec/mld_c_as_smoother.f90 +++ b/mlprec/mld_c_as_smoother.f90 @@ -243,8 +243,6 @@ contains subroutine c_as_smoother_descr(sm,info,iout,coarse) - use psb_base_mod - Implicit None ! Arguments diff --git a/mlprec/mld_c_base_solver_mod.f90 b/mlprec/mld_c_base_solver_mod.f90 index 9b0057de..ed21cf64 100644 --- a/mlprec/mld_c_base_solver_mod.f90 +++ b/mlprec/mld_c_base_solver_mod.f90 @@ -55,6 +55,8 @@ module mld_c_base_solver_mod use mld_base_prec_type use psb_base_mod, only : psb_desc_type, psb_cspmat_type, psb_long_int_k_, & + & psb_sizeof, psb_free, psb_cdfree, psb_errpush, psb_act_abort_,& + & psb_erractionsave, psb_erractionrestore, psb_error, psb_get_errstatus, psb_success_,& & psb_c_vect_type, psb_c_base_vect_type, psb_c_base_sparse_mat, psb_spk_ ! ! diff --git a/mlprec/mld_c_jac_smoother.f90 b/mlprec/mld_c_jac_smoother.f90 index ed46bdaf..a3374d9c 100644 --- a/mlprec/mld_c_jac_smoother.f90 +++ b/mlprec/mld_c_jac_smoother.f90 @@ -54,9 +54,9 @@ module mld_c_jac_smoother type(psb_cspmat_type) :: nd integer :: nnz_nd_tot contains - procedure, pass(sm) :: build => c_jac_smoother_bld - procedure, pass(sm) :: apply_v => c_jac_smoother_apply_vect - procedure, pass(sm) :: apply_a => c_jac_smoother_apply + procedure, pass(sm) :: build => mld_c_jac_smoother_bld + procedure, pass(sm) :: apply_v => mld_c_jac_smoother_apply_vect + procedure, pass(sm) :: apply_a => mld_c_jac_smoother_apply procedure, pass(sm) :: free => c_jac_smoother_free procedure, pass(sm) :: seti => c_jac_smoother_seti procedure, pass(sm) :: setc => c_jac_smoother_setc @@ -67,408 +67,63 @@ module mld_c_jac_smoother end type mld_c_jac_smoother_type - private :: c_jac_smoother_bld, c_jac_smoother_apply, & - & c_jac_smoother_free, c_jac_smoother_seti, & + private :: c_jac_smoother_free, c_jac_smoother_seti, & & c_jac_smoother_setc, c_jac_smoother_setr,& & c_jac_smoother_descr, c_jac_smoother_sizeof, & - & c_jac_smoother_apply_vect, c_jac_smoother_get_nzeros - - - + & c_jac_smoother_get_nzeros + + + interface mld_c_jac_smoother_apply_vect + subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) + import :: psb_desc_type, mld_c_jac_smoother_type, psb_c_vect_type, psb_spk_, & + & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type + + type(psb_desc_type), intent(in) :: desc_data + class(mld_c_jac_smoother_type), intent(inout) :: sm + type(psb_c_vect_type),intent(inout) :: x + type(psb_c_vect_type),intent(inout) :: y + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + integer, intent(in) :: sweeps + complex(psb_spk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + end subroutine mld_c_jac_smoother_apply_vect + end interface mld_c_jac_smoother_apply_vect + + interface mld_c_jac_smoother_apply + subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) + import :: psb_desc_type, mld_c_jac_smoother_type, psb_c_vect_type, psb_spk_, & + & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type + type(psb_desc_type), intent(in) :: desc_data + class(mld_c_jac_smoother_type), intent(in) :: sm + complex(psb_spk_),intent(inout) :: x(:) + complex(psb_spk_),intent(inout) :: y(:) + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + integer, intent(in) :: sweeps + complex(psb_spk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + end subroutine mld_c_jac_smoother_apply + end interface mld_c_jac_smoother_apply + + interface mld_c_jac_smoother_bld + subroutine mld_c_jac_smoother_bld(a,desc_a,sm,upd,info,amold,vmold) + import :: psb_desc_type, mld_c_jac_smoother_type, psb_c_vect_type, psb_spk_, & + & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type + type(psb_cspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_c_jac_smoother_type), intent(inout) :: sm + character, intent(in) :: upd + integer, intent(out) :: info + class(psb_c_base_sparse_mat), intent(in), optional :: amold + class(psb_c_base_vect_type), intent(in), optional :: vmold + end subroutine mld_c_jac_smoother_bld + end interface mld_c_jac_smoother_bld + contains - - subroutine c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) - use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_c_jac_smoother_type), intent(inout) :: sm - type(psb_c_vect_type),intent(inout) :: x - type(psb_c_vect_type),intent(inout) :: y - complex(psb_spk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - integer, intent(in) :: sweeps - complex(psb_spk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - integer :: n_row,n_col - type(psb_c_vect_type) :: tx, ty - complex(psb_spk_), pointer :: ww(:), aux(:) - integer :: ictxt,np,me,i, err_act - character :: trans_ - character(len=20) :: name='c_jac_smoother_apply' - - call psb_erractionsave(err_act) - - 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 - - if (.not.allocated(sm%sv)) then - info = 1121 - call psb_errpush(info,name) - goto 9999 - end if - - 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,0,0,0,0/),& - & a_err='complex(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='complex(psb_spk_)') - goto 9999 - end if - endif - - if ((sweeps == 1).or.(sm%nnz_nd_tot==0)) 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. - ! - ! - call tx%bld(x%get_nrows(),mold=x%v) - call tx%set(czero) - call ty%bld(x%get_nrows(),mold=x%v) - - 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. - ! - call psb_geaxpby(cone,x,czero,ty,desc_data,info) - call psb_spmm(-cone,sm%nd,tx,cone,ty,desc_data,info,work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(cone,ty,czero,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 - - call tx%free(info) - if (info == psb_success_) call ty%free(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 c_jac_smoother_apply_vect - - subroutine c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) - use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_c_jac_smoother_type), intent(in) :: sm - complex(psb_spk_),intent(inout) :: x(:) - complex(psb_spk_),intent(inout) :: y(:) - complex(psb_spk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - integer, intent(in) :: sweeps - complex(psb_spk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - integer :: n_row,n_col - complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character :: trans_ - character(len=20) :: name='c_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 = 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,0,0,0,0/),& - & a_err='complex(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='complex(psb_spk_)') - goto 9999 - end if - endif - - if ((sweeps == 1).or.(sm%nnz_nd_tot==0)) 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='complex(psb_spk_)') - goto 9999 - end if - - tx = czero - ty = czero - 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(-cone,sm%nd,tx,cone,ty,desc_data,info,work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(cone,ty,czero,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 c_jac_smoother_apply - - subroutine c_jac_smoother_bld(a,desc_a,sm,upd,info,amold,vmold) - - use psb_base_mod - use mld_c_diag_solver - Implicit None - - ! Arguments - type(psb_cspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(mld_c_jac_smoother_type), intent(inout) :: sm - character, intent(in) :: upd - integer, intent(out) :: info - class(psb_c_base_sparse_mat), intent(in), optional :: amold - class(psb_c_base_vect_type), intent(in), optional :: vmold - ! Local variables - integer :: n_row,n_col, nrow_a, nztota, nzeros - complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act, debug_unit, debug_level - character(len=20) :: name='c_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 = 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() - - nrow_a = a%get_nrows() - nztota = a%get_nzeros() - select type (smsv => sm%sv) - type is (mld_c_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_) then - if (present(amold)) then - call sm%nd%cscnv(info,& - & mold=amold,dupl=psb_dupl_add_) - else - call sm%nd%cscnv(info,& - & type='csr',dupl=psb_dupl_add_) - endif - end if - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='clip & psb_spcnv csr 4') - goto 9999 - end if - - call sm%sv%build(a,desc_a,upd,info,amold=amold,vmold=vmold) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='solver build') - goto 9999 - end if - nzeros = sm%nd%get_nzeros() - call psb_sum(ictxt,nzeros) - sm%nnz_nd_tot = nzeros - 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 c_jac_smoother_bld - - subroutine c_jac_smoother_seti(sm,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -483,14 +138,11 @@ contains call psb_erractionsave(err_act) select case(what) -!!$ case(mld_smoother_sweeps_) -!!$ sm%sweeps = val +! !$ 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 @@ -508,8 +160,6 @@ contains subroutine c_jac_smoother_setc(sm,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -546,7 +196,6 @@ contains subroutine c_jac_smoother_setr(sm,what,val,info) - use psb_base_mod Implicit None @@ -564,9 +213,6 @@ contains 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) @@ -583,7 +229,6 @@ contains subroutine c_jac_smoother_free(sm,info) - use psb_base_mod Implicit None @@ -623,7 +268,6 @@ contains subroutine c_jac_smoother_descr(sm,info,iout,coarse) - use psb_base_mod Implicit None @@ -674,7 +318,7 @@ contains end subroutine c_jac_smoother_descr function c_jac_smoother_sizeof(sm) result(val) - use psb_base_mod + implicit none ! Arguments class(mld_c_jac_smoother_type), intent(in) :: sm @@ -689,7 +333,7 @@ contains end function c_jac_smoother_sizeof function c_jac_smoother_get_nzeros(sm) result(val) - use psb_base_mod + implicit none ! Arguments class(mld_c_jac_smoother_type), intent(in) :: sm diff --git a/mlprec/mld_c_prec_type.f90 b/mlprec/mld_c_prec_type.f90 index a6ad8144..5fbcf07e 100644 --- a/mlprec/mld_c_prec_type.f90 +++ b/mlprec/mld_c_prec_type.f90 @@ -117,9 +117,8 @@ module mld_c_prec_type interface mld_precaply subroutine mld_cprecaply_vect(prec,x,y,desc_data,info,trans,work) - use psb_base_mod, only : psb_cspmat_type, psb_desc_type, & - & psb_spk_, psb_c_vect_type - import mld_cprec_type + import :: psb_cspmat_type, psb_desc_type, & + & psb_spk_, psb_c_vect_type, mld_cprec_type type(psb_desc_type),intent(in) :: desc_data type(mld_cprec_type), intent(inout) :: prec type(psb_c_vect_type),intent(inout) :: x @@ -129,8 +128,7 @@ module mld_c_prec_type complex(psb_spk_),intent(inout), optional, target :: work(:) end subroutine mld_cprecaply_vect subroutine mld_cprecaply(prec,x,y,desc_data,info,trans,work) - use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_ - import mld_cprec_type + import :: psb_cspmat_type, psb_desc_type, psb_spk_, mld_cprec_type type(psb_desc_type),intent(in) :: desc_data type(mld_cprec_type), intent(in) :: prec complex(psb_spk_),intent(inout) :: x(:) @@ -140,8 +138,7 @@ module mld_c_prec_type complex(psb_spk_),intent(inout), optional, target :: work(:) end subroutine mld_cprecaply subroutine mld_cprecaply1(prec,x,desc_data,info,trans) - use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_ - import mld_cprec_type + import :: psb_cspmat_type, psb_desc_type, psb_spk_, mld_cprec_type type(psb_desc_type),intent(in) :: desc_data type(mld_cprec_type), intent(in) :: prec complex(psb_spk_),intent(inout) :: x(:) diff --git a/mlprec/mld_d_as_smoother.f90 b/mlprec/mld_d_as_smoother.f90 index b06add65..fede0318 100644 --- a/mlprec/mld_d_as_smoother.f90 +++ b/mlprec/mld_d_as_smoother.f90 @@ -243,8 +243,6 @@ contains subroutine d_as_smoother_descr(sm,info,iout,coarse) - use psb_base_mod - Implicit None ! Arguments diff --git a/mlprec/mld_d_base_solver_mod.f90 b/mlprec/mld_d_base_solver_mod.f90 index 173a3b58..6c9be740 100644 --- a/mlprec/mld_d_base_solver_mod.f90 +++ b/mlprec/mld_d_base_solver_mod.f90 @@ -55,6 +55,8 @@ module mld_d_base_solver_mod use mld_base_prec_type use psb_base_mod, only : psb_desc_type, psb_dspmat_type, psb_long_int_k_, & + & psb_sizeof, psb_free, psb_cdfree, psb_errpush, psb_act_abort_,& + & psb_erractionsave, psb_erractionrestore, psb_error, psb_get_errstatus, psb_success_,& & psb_d_vect_type, psb_d_base_vect_type, psb_d_base_sparse_mat, psb_dpk_ ! ! diff --git a/mlprec/mld_d_jac_smoother.f90 b/mlprec/mld_d_jac_smoother.f90 index 3e865828..cc8e1f15 100644 --- a/mlprec/mld_d_jac_smoother.f90 +++ b/mlprec/mld_d_jac_smoother.f90 @@ -54,9 +54,9 @@ module mld_d_jac_smoother type(psb_dspmat_type) :: nd integer :: nnz_nd_tot contains - procedure, pass(sm) :: build => d_jac_smoother_bld - procedure, pass(sm) :: apply_v => d_jac_smoother_apply_vect - procedure, pass(sm) :: apply_a => d_jac_smoother_apply + procedure, pass(sm) :: build => mld_d_jac_smoother_bld + procedure, pass(sm) :: apply_v => mld_d_jac_smoother_apply_vect + procedure, pass(sm) :: apply_a => mld_d_jac_smoother_apply procedure, pass(sm) :: free => d_jac_smoother_free procedure, pass(sm) :: seti => d_jac_smoother_seti procedure, pass(sm) :: setc => d_jac_smoother_setc @@ -67,408 +67,63 @@ module mld_d_jac_smoother end type mld_d_jac_smoother_type - private :: d_jac_smoother_bld, d_jac_smoother_apply, & - & d_jac_smoother_free, d_jac_smoother_seti, & + private :: d_jac_smoother_free, d_jac_smoother_seti, & & d_jac_smoother_setc, d_jac_smoother_setr,& & d_jac_smoother_descr, d_jac_smoother_sizeof, & - & d_jac_smoother_apply_vect, d_jac_smoother_get_nzeros - - - + & d_jac_smoother_get_nzeros + + + interface mld_d_jac_smoother_apply_vect + subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) + import :: psb_desc_type, mld_d_jac_smoother_type, psb_d_vect_type, psb_dpk_, & + & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type + + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_jac_smoother_type), intent(inout) :: sm + 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 + integer, intent(in) :: sweeps + real(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + end subroutine mld_d_jac_smoother_apply_vect + end interface mld_d_jac_smoother_apply_vect + + interface mld_d_jac_smoother_apply + subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) + import :: psb_desc_type, mld_d_jac_smoother_type, psb_d_vect_type, psb_dpk_, & + & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_jac_smoother_type), intent(in) :: sm + 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 + integer, intent(in) :: sweeps + real(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + end subroutine mld_d_jac_smoother_apply + end interface mld_d_jac_smoother_apply + + interface mld_d_jac_smoother_bld + subroutine mld_d_jac_smoother_bld(a,desc_a,sm,upd,info,amold,vmold) + import :: psb_desc_type, mld_d_jac_smoother_type, psb_d_vect_type, psb_dpk_, & + & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type + type(psb_dspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_d_jac_smoother_type), intent(inout) :: sm + character, intent(in) :: upd + integer, intent(out) :: info + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold + end subroutine mld_d_jac_smoother_bld + end interface mld_d_jac_smoother_bld + contains - - subroutine d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) - use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_d_jac_smoother_type), intent(inout) :: sm - 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 - integer, intent(in) :: sweeps - real(psb_dpk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - integer :: n_row,n_col - type(psb_d_vect_type) :: tx, ty - real(psb_dpk_), pointer :: ww(:), aux(:) - 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') - case('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 = 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,0,0,0,0/),& - & 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,0,0,0,0/),& - & a_err='real(psb_dpk_)') - goto 9999 - end if - endif - - if ((sweeps == 1).or.(sm%nnz_nd_tot==0)) 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. - ! - ! - call tx%bld(x%get_nrows(),mold=x%v) - call tx%set(dzero) - call ty%bld(x%get_nrows(),mold=x%v) - - 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. - ! - call psb_geaxpby(done,x,dzero,ty,desc_data,info) - call psb_spmm(-done,sm%nd,tx,done,ty,desc_data,info,work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(done,ty,dzero,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 - - call tx%free(info) - if (info == psb_success_) call ty%free(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 d_jac_smoother_apply_vect - - subroutine d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) - use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_d_jac_smoother_type), intent(in) :: sm - 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 - integer, intent(in) :: sweeps - real(psb_dpk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - integer :: n_row,n_col - real(psb_dpk_), 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 = 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,0,0,0,0/),& - & 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,0,0,0,0/),& - & a_err='real(psb_dpk_)') - goto 9999 - end if - endif - - if ((sweeps == 1).or.(sm%nnz_nd_tot==0)) 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_dpk_)') - goto 9999 - end if - - tx = dzero - ty = dzero - 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(-done,sm%nd,tx,done,ty,desc_data,info,work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(done,ty,dzero,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 d_jac_smoother_apply - - subroutine d_jac_smoother_bld(a,desc_a,sm,upd,info,amold,vmold) - - use psb_base_mod - use mld_d_diag_solver - Implicit None - - ! Arguments - type(psb_dspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(mld_d_jac_smoother_type), intent(inout) :: sm - character, intent(in) :: upd - integer, intent(out) :: info - class(psb_d_base_sparse_mat), intent(in), optional :: amold - class(psb_d_base_vect_type), intent(in), optional :: vmold - ! Local variables - integer :: n_row,n_col, nrow_a, nztota, nzeros - real(psb_dpk_), 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 = 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() - - nrow_a = a%get_nrows() - nztota = a%get_nzeros() - select type (smsv => sm%sv) - type is (mld_d_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_) then - if (present(amold)) then - call sm%nd%cscnv(info,& - & mold=amold,dupl=psb_dupl_add_) - else - call sm%nd%cscnv(info,& - & type='csr',dupl=psb_dupl_add_) - endif - end if - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='clip & psb_spcnv csr 4') - goto 9999 - end if - - call sm%sv%build(a,desc_a,upd,info,amold=amold,vmold=vmold) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='solver build') - goto 9999 - end if - nzeros = sm%nd%get_nzeros() - call psb_sum(ictxt,nzeros) - sm%nnz_nd_tot = nzeros - 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 d_jac_smoother_bld - - subroutine d_jac_smoother_seti(sm,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -483,14 +138,11 @@ contains call psb_erractionsave(err_act) select case(what) -!!$ case(mld_smoother_sweeps_) -!!$ sm%sweeps = val +! !$ 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 @@ -508,8 +160,6 @@ contains subroutine d_jac_smoother_setc(sm,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -546,7 +196,6 @@ contains subroutine d_jac_smoother_setr(sm,what,val,info) - use psb_base_mod Implicit None @@ -564,9 +213,6 @@ contains 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) @@ -583,7 +229,6 @@ contains subroutine d_jac_smoother_free(sm,info) - use psb_base_mod Implicit None @@ -623,7 +268,6 @@ contains subroutine d_jac_smoother_descr(sm,info,iout,coarse) - use psb_base_mod Implicit None @@ -674,7 +318,7 @@ contains end subroutine d_jac_smoother_descr function d_jac_smoother_sizeof(sm) result(val) - use psb_base_mod + implicit none ! Arguments class(mld_d_jac_smoother_type), intent(in) :: sm @@ -689,7 +333,7 @@ contains end function d_jac_smoother_sizeof function d_jac_smoother_get_nzeros(sm) result(val) - use psb_base_mod + implicit none ! Arguments class(mld_d_jac_smoother_type), intent(in) :: sm diff --git a/mlprec/mld_d_prec_type.f90 b/mlprec/mld_d_prec_type.f90 index 2c05be91..33ff7333 100644 --- a/mlprec/mld_d_prec_type.f90 +++ b/mlprec/mld_d_prec_type.f90 @@ -117,9 +117,8 @@ module mld_d_prec_type interface mld_precaply subroutine mld_dprecaply_vect(prec,x,y,desc_data,info,trans,work) - use psb_base_mod, only : psb_dspmat_type, psb_desc_type, & - & psb_dpk_, psb_d_vect_type - import mld_dprec_type + import :: psb_dspmat_type, psb_desc_type, & + & psb_dpk_, psb_d_vect_type, mld_dprec_type type(psb_desc_type),intent(in) :: desc_data type(mld_dprec_type), intent(inout) :: prec type(psb_d_vect_type),intent(inout) :: x @@ -129,8 +128,7 @@ module mld_d_prec_type real(psb_dpk_),intent(inout), optional, target :: work(:) end subroutine mld_dprecaply_vect subroutine mld_dprecaply(prec,x,y,desc_data,info,trans,work) - use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_ - import mld_dprec_type + import :: psb_dspmat_type, psb_desc_type, psb_dpk_, mld_dprec_type type(psb_desc_type),intent(in) :: desc_data type(mld_dprec_type), intent(in) :: prec real(psb_dpk_),intent(inout) :: x(:) @@ -140,8 +138,7 @@ module mld_d_prec_type real(psb_dpk_),intent(inout), optional, target :: work(:) end subroutine mld_dprecaply subroutine mld_dprecaply1(prec,x,desc_data,info,trans) - use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_ - import mld_dprec_type + import :: psb_dspmat_type, psb_desc_type, psb_dpk_, mld_dprec_type type(psb_desc_type),intent(in) :: desc_data type(mld_dprec_type), intent(in) :: prec real(psb_dpk_),intent(inout) :: x(:) diff --git a/mlprec/mld_s_as_smoother.f90 b/mlprec/mld_s_as_smoother.f90 index 7bbe3310..b023a31e 100644 --- a/mlprec/mld_s_as_smoother.f90 +++ b/mlprec/mld_s_as_smoother.f90 @@ -243,8 +243,6 @@ contains subroutine s_as_smoother_descr(sm,info,iout,coarse) - use psb_base_mod - Implicit None ! Arguments diff --git a/mlprec/mld_s_base_solver_mod.f90 b/mlprec/mld_s_base_solver_mod.f90 index 094c1e93..9efc12f1 100644 --- a/mlprec/mld_s_base_solver_mod.f90 +++ b/mlprec/mld_s_base_solver_mod.f90 @@ -55,6 +55,8 @@ module mld_s_base_solver_mod use mld_base_prec_type use psb_base_mod, only : psb_desc_type, psb_sspmat_type, psb_long_int_k_, & + & psb_sizeof, psb_free, psb_cdfree, psb_errpush, psb_act_abort_,& + & psb_erractionsave, psb_erractionrestore, psb_error, psb_get_errstatus, psb_success_,& & psb_s_vect_type, psb_s_base_vect_type, psb_s_base_sparse_mat, psb_spk_ ! ! diff --git a/mlprec/mld_s_jac_smoother.f90 b/mlprec/mld_s_jac_smoother.f90 index 62b52f20..594d49e1 100644 --- a/mlprec/mld_s_jac_smoother.f90 +++ b/mlprec/mld_s_jac_smoother.f90 @@ -54,9 +54,9 @@ module mld_s_jac_smoother type(psb_sspmat_type) :: nd integer :: nnz_nd_tot contains - procedure, pass(sm) :: build => s_jac_smoother_bld - procedure, pass(sm) :: apply_v => s_jac_smoother_apply_vect - procedure, pass(sm) :: apply_a => s_jac_smoother_apply + procedure, pass(sm) :: build => mld_s_jac_smoother_bld + procedure, pass(sm) :: apply_v => mld_s_jac_smoother_apply_vect + procedure, pass(sm) :: apply_a => mld_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 @@ -67,408 +67,63 @@ module mld_s_jac_smoother end type mld_s_jac_smoother_type - private :: s_jac_smoother_bld, s_jac_smoother_apply, & - & s_jac_smoother_free, s_jac_smoother_seti, & + private :: s_jac_smoother_free, s_jac_smoother_seti, & & s_jac_smoother_setc, s_jac_smoother_setr,& & s_jac_smoother_descr, s_jac_smoother_sizeof, & - & s_jac_smoother_apply_vect, s_jac_smoother_get_nzeros - - - + & s_jac_smoother_get_nzeros + + + interface mld_s_jac_smoother_apply_vect + subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) + import :: psb_desc_type, mld_s_jac_smoother_type, psb_s_vect_type, psb_spk_, & + & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type + + type(psb_desc_type), intent(in) :: desc_data + class(mld_s_jac_smoother_type), intent(inout) :: sm + type(psb_s_vect_type),intent(inout) :: x + type(psb_s_vect_type),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 + end subroutine mld_s_jac_smoother_apply_vect + end interface mld_s_jac_smoother_apply_vect + + interface mld_s_jac_smoother_apply + subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) + import :: psb_desc_type, mld_s_jac_smoother_type, psb_s_vect_type, psb_spk_, & + & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type + type(psb_desc_type), intent(in) :: desc_data + class(mld_s_jac_smoother_type), intent(in) :: sm + real(psb_spk_),intent(inout) :: 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 + end subroutine mld_s_jac_smoother_apply + end interface mld_s_jac_smoother_apply + + interface mld_s_jac_smoother_bld + subroutine mld_s_jac_smoother_bld(a,desc_a,sm,upd,info,amold,vmold) + import :: psb_desc_type, mld_s_jac_smoother_type, psb_s_vect_type, psb_spk_, & + & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type + 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 + class(psb_s_base_sparse_mat), intent(in), optional :: amold + class(psb_s_base_vect_type), intent(in), optional :: vmold + end subroutine mld_s_jac_smoother_bld + end interface mld_s_jac_smoother_bld + contains - - subroutine s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) - use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_s_jac_smoother_type), intent(inout) :: sm - type(psb_s_vect_type),intent(inout) :: x - type(psb_s_vect_type),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 - type(psb_s_vect_type) :: tx, ty - real(psb_spk_), pointer :: ww(:), aux(:) - integer :: ictxt,np,me,i, err_act - character :: trans_ - character(len=20) :: name='s_jac_smoother_apply' - - call psb_erractionsave(err_act) - - 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 - - if (.not.allocated(sm%sv)) then - info = 1121 - call psb_errpush(info,name) - goto 9999 - end if - - 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,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).or.(sm%nnz_nd_tot==0)) 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. - ! - ! - call tx%bld(x%get_nrows(),mold=x%v) - call tx%set(szero) - call ty%bld(x%get_nrows(),mold=x%v) - - 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. - ! - call psb_geaxpby(sone,x,szero,ty,desc_data,info) - 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 - - call tx%free(info) - if (info == psb_success_) call ty%free(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_vect - - subroutine s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) - use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_s_jac_smoother_type), intent(in) :: sm - real(psb_spk_),intent(inout) :: 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='s_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 = 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,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).or.(sm%nnz_nd_tot==0)) 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,amold,vmold) - - use psb_base_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 - class(psb_s_base_sparse_mat), intent(in), optional :: amold - class(psb_s_base_vect_type), intent(in), optional :: vmold - ! Local variables - integer :: n_row,n_col, nrow_a, nztota, nzeros - real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act, debug_unit, debug_level - character(len=20) :: name='s_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 = 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() - - 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_) then - if (present(amold)) then - call sm%nd%cscnv(info,& - & mold=amold,dupl=psb_dupl_add_) - else - call sm%nd%cscnv(info,& - & type='csr',dupl=psb_dupl_add_) - endif - end if - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='clip & psb_spcnv csr 4') - goto 9999 - end if - - call sm%sv%build(a,desc_a,upd,info,amold=amold,vmold=vmold) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='solver build') - goto 9999 - end if - nzeros = sm%nd%get_nzeros() - call psb_sum(ictxt,nzeros) - sm%nnz_nd_tot = nzeros - 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_base_mod - Implicit None ! Arguments @@ -483,14 +138,11 @@ contains call psb_erractionsave(err_act) select case(what) -!!$ case(mld_smoother_sweeps_) -!!$ sm%sweeps = val +! !$ 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 @@ -508,8 +160,6 @@ contains subroutine s_jac_smoother_setc(sm,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -546,7 +196,6 @@ contains subroutine s_jac_smoother_setr(sm,what,val,info) - use psb_base_mod Implicit None @@ -564,9 +213,6 @@ contains 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) @@ -583,7 +229,6 @@ contains subroutine s_jac_smoother_free(sm,info) - use psb_base_mod Implicit None @@ -623,7 +268,6 @@ contains subroutine s_jac_smoother_descr(sm,info,iout,coarse) - use psb_base_mod Implicit None @@ -674,7 +318,7 @@ contains end subroutine s_jac_smoother_descr function s_jac_smoother_sizeof(sm) result(val) - use psb_base_mod + implicit none ! Arguments class(mld_s_jac_smoother_type), intent(in) :: sm @@ -689,7 +333,7 @@ contains end function s_jac_smoother_sizeof function s_jac_smoother_get_nzeros(sm) result(val) - use psb_base_mod + implicit none ! Arguments class(mld_s_jac_smoother_type), intent(in) :: sm diff --git a/mlprec/mld_s_prec_type.f90 b/mlprec/mld_s_prec_type.f90 index b462a128..6ed1e80f 100644 --- a/mlprec/mld_s_prec_type.f90 +++ b/mlprec/mld_s_prec_type.f90 @@ -117,9 +117,8 @@ module mld_s_prec_type interface mld_precaply subroutine mld_sprecaply_vect(prec,x,y,desc_data,info,trans,work) - use psb_base_mod, only : psb_sspmat_type, psb_desc_type, & - & psb_spk_, psb_s_vect_type - import mld_sprec_type + import :: psb_sspmat_type, psb_desc_type, & + & psb_spk_, psb_s_vect_type, mld_sprec_type type(psb_desc_type),intent(in) :: desc_data type(mld_sprec_type), intent(inout) :: prec type(psb_s_vect_type),intent(inout) :: x @@ -129,8 +128,7 @@ module mld_s_prec_type real(psb_spk_),intent(inout), optional, target :: work(:) end subroutine mld_sprecaply_vect subroutine mld_sprecaply(prec,x,y,desc_data,info,trans,work) - use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_ - import mld_sprec_type + import :: psb_sspmat_type, psb_desc_type, psb_spk_, mld_sprec_type type(psb_desc_type),intent(in) :: desc_data type(mld_sprec_type), intent(in) :: prec real(psb_spk_),intent(inout) :: x(:) @@ -140,8 +138,7 @@ module mld_s_prec_type real(psb_spk_),intent(inout), optional, target :: work(:) end subroutine mld_sprecaply subroutine mld_sprecaply1(prec,x,desc_data,info,trans) - use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_ - import mld_sprec_type + import :: psb_sspmat_type, psb_desc_type, psb_spk_, mld_sprec_type type(psb_desc_type),intent(in) :: desc_data type(mld_sprec_type), intent(in) :: prec real(psb_spk_),intent(inout) :: x(:) diff --git a/mlprec/mld_z_as_smoother.f90 b/mlprec/mld_z_as_smoother.f90 index a0bc2f3d..9b7455e6 100644 --- a/mlprec/mld_z_as_smoother.f90 +++ b/mlprec/mld_z_as_smoother.f90 @@ -243,8 +243,6 @@ contains subroutine z_as_smoother_descr(sm,info,iout,coarse) - use psb_base_mod - Implicit None ! Arguments diff --git a/mlprec/mld_z_base_solver_mod.f90 b/mlprec/mld_z_base_solver_mod.f90 index e743b5b0..3244d328 100644 --- a/mlprec/mld_z_base_solver_mod.f90 +++ b/mlprec/mld_z_base_solver_mod.f90 @@ -55,6 +55,8 @@ module mld_z_base_solver_mod use mld_base_prec_type use psb_base_mod, only : psb_desc_type, psb_zspmat_type, psb_long_int_k_, & + & psb_sizeof, psb_free, psb_cdfree, psb_errpush, psb_act_abort_,& + & psb_erractionsave, psb_erractionrestore, psb_error, psb_get_errstatus, psb_success_,& & psb_z_vect_type, psb_z_base_vect_type, psb_z_base_sparse_mat, psb_dpk_ ! ! diff --git a/mlprec/mld_z_jac_smoother.f90 b/mlprec/mld_z_jac_smoother.f90 index f1ccef09..499a0a42 100644 --- a/mlprec/mld_z_jac_smoother.f90 +++ b/mlprec/mld_z_jac_smoother.f90 @@ -54,9 +54,9 @@ module mld_z_jac_smoother type(psb_zspmat_type) :: nd integer :: nnz_nd_tot contains - procedure, pass(sm) :: build => z_jac_smoother_bld - procedure, pass(sm) :: apply_v => z_jac_smoother_apply_vect - procedure, pass(sm) :: apply_a => z_jac_smoother_apply + procedure, pass(sm) :: build => mld_z_jac_smoother_bld + procedure, pass(sm) :: apply_v => mld_z_jac_smoother_apply_vect + procedure, pass(sm) :: apply_a => mld_z_jac_smoother_apply procedure, pass(sm) :: free => z_jac_smoother_free procedure, pass(sm) :: seti => z_jac_smoother_seti procedure, pass(sm) :: setc => z_jac_smoother_setc @@ -67,408 +67,63 @@ module mld_z_jac_smoother end type mld_z_jac_smoother_type - private :: z_jac_smoother_bld, z_jac_smoother_apply, & - & z_jac_smoother_free, z_jac_smoother_seti, & + private :: z_jac_smoother_free, z_jac_smoother_seti, & & z_jac_smoother_setc, z_jac_smoother_setr,& & z_jac_smoother_descr, z_jac_smoother_sizeof, & - & z_jac_smoother_apply_vect, z_jac_smoother_get_nzeros - - - + & z_jac_smoother_get_nzeros + + + interface mld_z_jac_smoother_apply_vect + subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) + import :: psb_desc_type, mld_z_jac_smoother_type, psb_z_vect_type, psb_dpk_, & + & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type + + type(psb_desc_type), intent(in) :: desc_data + class(mld_z_jac_smoother_type), intent(inout) :: sm + type(psb_z_vect_type),intent(inout) :: x + type(psb_z_vect_type),intent(inout) :: y + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + integer, intent(in) :: sweeps + complex(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + end subroutine mld_z_jac_smoother_apply_vect + end interface mld_z_jac_smoother_apply_vect + + interface mld_z_jac_smoother_apply + subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) + import :: psb_desc_type, mld_z_jac_smoother_type, psb_z_vect_type, psb_dpk_, & + & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type + type(psb_desc_type), intent(in) :: desc_data + class(mld_z_jac_smoother_type), intent(in) :: sm + complex(psb_dpk_),intent(inout) :: x(:) + complex(psb_dpk_),intent(inout) :: y(:) + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + integer, intent(in) :: sweeps + complex(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + end subroutine mld_z_jac_smoother_apply + end interface mld_z_jac_smoother_apply + + interface mld_z_jac_smoother_bld + subroutine mld_z_jac_smoother_bld(a,desc_a,sm,upd,info,amold,vmold) + import :: psb_desc_type, mld_z_jac_smoother_type, psb_z_vect_type, psb_dpk_, & + & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type + type(psb_zspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_z_jac_smoother_type), intent(inout) :: sm + character, intent(in) :: upd + integer, intent(out) :: info + class(psb_z_base_sparse_mat), intent(in), optional :: amold + class(psb_z_base_vect_type), intent(in), optional :: vmold + end subroutine mld_z_jac_smoother_bld + end interface mld_z_jac_smoother_bld + contains - - subroutine z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) - use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_z_jac_smoother_type), intent(inout) :: sm - type(psb_z_vect_type),intent(inout) :: x - type(psb_z_vect_type),intent(inout) :: y - complex(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - integer, intent(in) :: sweeps - complex(psb_dpk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - integer :: n_row,n_col - type(psb_z_vect_type) :: tx, ty - complex(psb_dpk_), pointer :: ww(:), aux(:) - integer :: ictxt,np,me,i, err_act - character :: trans_ - character(len=20) :: name='z_jac_smoother_apply' - - call psb_erractionsave(err_act) - - 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 - - if (.not.allocated(sm%sv)) then - info = 1121 - call psb_errpush(info,name) - goto 9999 - end if - - 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,0,0,0,0/),& - & a_err='complex(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,0,0,0,0/),& - & a_err='complex(psb_dpk_)') - goto 9999 - end if - endif - - if ((sweeps == 1).or.(sm%nnz_nd_tot==0)) 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. - ! - ! - call tx%bld(x%get_nrows(),mold=x%v) - call tx%set(zzero) - call ty%bld(x%get_nrows(),mold=x%v) - - 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. - ! - call psb_geaxpby(zone,x,zzero,ty,desc_data,info) - call psb_spmm(-zone,sm%nd,tx,zone,ty,desc_data,info,work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(zone,ty,zzero,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 - - call tx%free(info) - if (info == psb_success_) call ty%free(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 z_jac_smoother_apply_vect - - subroutine z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) - use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_z_jac_smoother_type), intent(in) :: sm - complex(psb_dpk_),intent(inout) :: x(:) - complex(psb_dpk_),intent(inout) :: y(:) - complex(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - integer, intent(in) :: sweeps - complex(psb_dpk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - integer :: n_row,n_col - complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character :: trans_ - character(len=20) :: name='z_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 = 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,0,0,0,0/),& - & a_err='complex(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,0,0,0,0/),& - & a_err='complex(psb_dpk_)') - goto 9999 - end if - endif - - if ((sweeps == 1).or.(sm%nnz_nd_tot==0)) 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='complex(psb_dpk_)') - goto 9999 - end if - - tx = zzero - ty = zzero - 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(-zone,sm%nd,tx,zone,ty,desc_data,info,work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(zone,ty,zzero,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 z_jac_smoother_apply - - subroutine z_jac_smoother_bld(a,desc_a,sm,upd,info,amold,vmold) - - use psb_base_mod - use mld_z_diag_solver - Implicit None - - ! Arguments - type(psb_zspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(mld_z_jac_smoother_type), intent(inout) :: sm - character, intent(in) :: upd - integer, intent(out) :: info - class(psb_z_base_sparse_mat), intent(in), optional :: amold - class(psb_z_base_vect_type), intent(in), optional :: vmold - ! Local variables - integer :: n_row,n_col, nrow_a, nztota, nzeros - complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act, debug_unit, debug_level - character(len=20) :: name='z_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 = 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() - - nrow_a = a%get_nrows() - nztota = a%get_nzeros() - select type (smsv => sm%sv) - type is (mld_z_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_) then - if (present(amold)) then - call sm%nd%cscnv(info,& - & mold=amold,dupl=psb_dupl_add_) - else - call sm%nd%cscnv(info,& - & type='csr',dupl=psb_dupl_add_) - endif - end if - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='clip & psb_spcnv csr 4') - goto 9999 - end if - - call sm%sv%build(a,desc_a,upd,info,amold=amold,vmold=vmold) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='solver build') - goto 9999 - end if - nzeros = sm%nd%get_nzeros() - call psb_sum(ictxt,nzeros) - sm%nnz_nd_tot = nzeros - 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 z_jac_smoother_bld - - subroutine z_jac_smoother_seti(sm,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -483,14 +138,11 @@ contains call psb_erractionsave(err_act) select case(what) -!!$ case(mld_smoother_sweeps_) -!!$ sm%sweeps = val +! !$ 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 @@ -508,8 +160,6 @@ contains subroutine z_jac_smoother_setc(sm,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -546,7 +196,6 @@ contains subroutine z_jac_smoother_setr(sm,what,val,info) - use psb_base_mod Implicit None @@ -564,9 +213,6 @@ contains 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) @@ -583,7 +229,6 @@ contains subroutine z_jac_smoother_free(sm,info) - use psb_base_mod Implicit None @@ -623,7 +268,6 @@ contains subroutine z_jac_smoother_descr(sm,info,iout,coarse) - use psb_base_mod Implicit None @@ -674,7 +318,7 @@ contains end subroutine z_jac_smoother_descr function z_jac_smoother_sizeof(sm) result(val) - use psb_base_mod + implicit none ! Arguments class(mld_z_jac_smoother_type), intent(in) :: sm @@ -689,7 +333,7 @@ contains end function z_jac_smoother_sizeof function z_jac_smoother_get_nzeros(sm) result(val) - use psb_base_mod + implicit none ! Arguments class(mld_z_jac_smoother_type), intent(in) :: sm diff --git a/mlprec/mld_z_prec_type.f90 b/mlprec/mld_z_prec_type.f90 index d9ceeaa4..fb0d6416 100644 --- a/mlprec/mld_z_prec_type.f90 +++ b/mlprec/mld_z_prec_type.f90 @@ -117,9 +117,8 @@ module mld_z_prec_type interface mld_precaply subroutine mld_zprecaply_vect(prec,x,y,desc_data,info,trans,work) - use psb_base_mod, only : psb_zspmat_type, psb_desc_type, & - & psb_dpk_, psb_z_vect_type - import mld_zprec_type + import :: psb_zspmat_type, psb_desc_type, & + & psb_dpk_, psb_z_vect_type, mld_zprec_type type(psb_desc_type),intent(in) :: desc_data type(mld_zprec_type), intent(inout) :: prec type(psb_z_vect_type),intent(inout) :: x @@ -129,8 +128,7 @@ module mld_z_prec_type complex(psb_dpk_),intent(inout), optional, target :: work(:) end subroutine mld_zprecaply_vect subroutine mld_zprecaply(prec,x,y,desc_data,info,trans,work) - use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_ - import mld_zprec_type + import :: psb_zspmat_type, psb_desc_type, psb_dpk_, mld_zprec_type type(psb_desc_type),intent(in) :: desc_data type(mld_zprec_type), intent(in) :: prec complex(psb_dpk_),intent(inout) :: x(:) @@ -140,8 +138,7 @@ module mld_z_prec_type complex(psb_dpk_),intent(inout), optional, target :: work(:) end subroutine mld_zprecaply subroutine mld_zprecaply1(prec,x,desc_data,info,trans) - use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_ - import mld_zprec_type + import :: psb_zspmat_type, psb_desc_type, psb_dpk_, mld_zprec_type type(psb_desc_type),intent(in) :: desc_data type(mld_zprec_type), intent(in) :: prec complex(psb_dpk_),intent(inout) :: x(:) diff --git a/tests/pdegen/runs/ppde.inp b/tests/pdegen/runs/ppde.inp index 39f0ee41..46fe8d3d 100644 --- a/tests/pdegen/runs/ppde.inp +++ b/tests/pdegen/runs/ppde.inp @@ -1,6 +1,6 @@ BICGSTAB ! Iterative method: BiCGSTAB BiCG CGS RGMRES BiCGSTABL CG CSR ! Storage format CSR COO JAD -080 ! IDIM; domain size is idim**3 +040 ! IDIM; domain size is idim**3 2 ! ISTOPC 0100 ! ITMAX -1 ! ITRACE