diff --git a/mlprec/Makefile b/mlprec/Makefile index 557cc78c..bf21b92e 100644 --- a/mlprec/Makefile +++ b/mlprec/Makefile @@ -33,55 +33,13 @@ INNEROBJS= mld_dcoarse_bld.o mld_dmlprec_bld.o mld_dslu_bld.o \ mld_zilu0_fact.o mld_ziluk_fact.o mld_zilut_fact.o mld_zaggrmap_bld.o \ mld_zmlprec_aply.o mld_zslud_bld.o mld_zaggrmat_asb.o \ $(MPFOBJS) -# -# MPFOBJS=mld_saggrmat_nosmth_asb.o mld_saggrmat_smth_asb.o \ -# mld_daggrmat_nosmth_asb.o mld_daggrmat_smth_asb.o mld_daggrmat_minnrg_asb.o\ -# mld_caggrmat_nosmth_asb.o mld_caggrmat_smth_asb.o \ -# mld_zaggrmat_nosmth_asb.o mld_zaggrmat_smth_asb.o -# MPCOBJS=mld_sslud_interface.o mld_dslud_interface.o mld_cslud_interface.o mld_zslud_interface.o -# INNEROBJS=mld_scoarse_bld.o mld_dcoarse_bld.o \ -# mld_ccoarse_bld.o mld_zcoarse_bld.o \ -# mld_smlprec_bld.o mld_dmlprec_bld.o mld_cmlprec_bld.o mld_zmlprec_bld.o\ -# mld_sas_bld.o mld_sslu_bld.o mld_sumf_bld.o mld_silu0_fact.o\ -# mld_ssp_renum.o mld_sfact_bld.o mld_silu_bld.o \ -# mld_sbaseprec_bld.o mld_sdiag_bld.o mld_saggrmap_bld.o \ -# mld_smlprec_aply.o mld_sslud_bld.o\ -# mld_sbaseprec_aply.o mld_ssub_aply.o mld_ssub_solve.o \ -# mld_sas_aply.o mld_saggrmat_asb.o \ -# mld_das_bld.o mld_dslu_bld.o mld_dumf_bld.o mld_dilu0_fact.o\ -# mld_dsp_renum.o mld_dfact_bld.o mld_dilu_bld.o \ -# mld_dbaseprec_bld.o mld_ddiag_bld.o mld_daggrmap_bld.o \ -# mld_dmlprec_aply.o mld_dslud_bld.o\ -# mld_dbaseprec_aply.o mld_dsub_aply.o mld_dsub_solve.o \ -# mld_das_aply.o mld_daggrmat_asb.o \ -# mld_cas_bld.o mld_cslu_bld.o mld_cumf_bld.o mld_cilu0_fact.o\ -# mld_csp_renum.o mld_cfact_bld.o mld_cilu_bld.o \ -# mld_cbaseprec_bld.o mld_cdiag_bld.o mld_caggrmap_bld.o \ -# mld_cmlprec_aply.o mld_cslud_bld.o\ -# mld_cbaseprec_aply.o mld_csub_aply.o mld_csub_solve.o \ -# mld_cas_aply.o mld_caggrmat_asb.o\ -# mld_zas_bld.o mld_zslu_bld.o mld_zumf_bld.o mld_zilu0_fact.o\ -# mld_zsp_renum.o mld_zfact_bld.o mld_zilu_bld.o \ -# mld_zbaseprec_bld.o mld_zdiag_bld.o mld_zaggrmap_bld.o \ -# mld_zmlprec_aply.o mld_zslud_bld.o\ -# mld_zbaseprec_aply.o mld_zsub_aply.o mld_zsub_solve.o \ -# mld_zas_aply.o mld_zaggrmat_asb.o\ -# mld_siluk_fact.o mld_ciluk_fact.o mld_silut_fact.o mld_cilut_fact.o \ -# mld_diluk_fact.o mld_ziluk_fact.o mld_dilut_fact.o mld_zilut_fact.o \ -# $(MPFOBJS) + OUTEROBJS=mld_dprecbld.o mld_dprecset.o mld_dprecinit.o mld_dprecaply.o \ mld_sprecbld.o mld_sprecset.o mld_sprecinit.o mld_sprecaply.o \ mld_cprecbld.o mld_cprecset.o mld_cprecinit.o mld_cprecaply.o \ mld_zprecbld.o mld_zprecset.o mld_zprecinit.o mld_zprecaply.o -#OUTEROBJS=mld_sprecbld.o mld_sprecset.o mld_sprecinit.o\ -# mld_sprecaply.o \ -# mld_dprecbld.o mld_dprecset.o mld_dprecinit.o\ -# mld_dprecaply.o \ -# mld_cprecbld.o mld_cprecset.o mld_cprecinit.o\ -# mld_cprecaply.o \ -# mld_zprecbld.o mld_zprecset.o mld_zprecinit.o \ -# mld_zprecaply.o + F90OBJS=$(OUTEROBJS) $(INNEROBJS) COBJS= mld_sslu_interface.o mld_sumf_interface.o \ @@ -109,11 +67,11 @@ mld_prec_type.o: mld_s_prec_type.o mld_d_prec_type.o mld_c_prec_type.o mld_z_pre mld_prec_mod.o mld_innner_mod.o: mld_prec_type.o mld_inner_mod.o: mld_move_alloc_mod.o mld_move_alloc_mod.o: mld_prec_type.o -mld_d_diag_solver.o mld_d_ilu_solver.o: mld_d_prec_type.o +mld_d_umf_solver.o mld_d_diag_solver.o mld_d_ilu_solver.o: mld_d_prec_type.o mld_d_as_smoother.o mld_d_jac_smoother.o: mld_d_prec_type.o mld_d_jac_smoother.o: mld_d_diag_solver.o mld_dprecinit.o mld_dprecset.o: mld_d_diag_solver.o mld_d_ilu_solver.o \ - mld_d_as_smoother.o mld_d_jac_smoother.o + mld_d_umf_solver.o mld_d_as_smoother.o mld_d_jac_smoother.o $(MODOBJS): $(PSBINCDIR)/psb_sparse_mod$(.mod) diff --git a/mlprec/mld_csub_aply.f90 b/mlprec/mld_csub_aply.f90 deleted file mode 100644 index b216c5b5..00000000 --- a/mlprec/mld_csub_aply.f90 +++ /dev/null @@ -1,296 +0,0 @@ -!!$ -!!$ -!!$ MLD2P4 version 2.0 -!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package -!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) -!!$ -!!$ (C) Copyright 2008,2009,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. -!!$ -!!$ -! File mld_csub_aply.f90 -! -! Subroutine: mld_csub_aply -! Version: complex -! -! This routine computes -! -! Y = beta*Y + alpha*op(K^(-1))*X, -! -! where -! - K is a suitable matrix, as specified below, -! - op(K^(-1)) is K^(-1) or its transpose, according to the value of the -! argument trans, -! - X and Y are vectors, -! - alpha and beta are scalars. -! -! Depending on K, alpha and beta (and on the communication descriptor desc_data -! - see the arguments below), the above computation may correspond to one of -! the following tasks: -! -! 1. Application of a block-Jacobi preconditioner associated to a matrix A -! distributed among the processes. Here K is the preconditioner, op(K^(-1)) -! = K^(-1), alpha = 1 and beta = 0. -! -! 2. Application of block-Jacobi sweeps to compute an approximate solution of -! a linear system -! A*Y = X, -! -! distributed among the processes (note that a single block-Jacobi sweep, -! with null starting guess, corresponds to the application of a block-Jacobi -! preconditioner). Here K^(-1) denotes the iteration matrix of the -! block-Jacobi solver, op(K^(-1)) = K^(-1), alpha = 1 and beta = 0. -! -! 3. Solution, through the LU factorization, of a linear system -! -! A*Y = X, -! -! distributed among the processes. Here K = L*U = A, op(K^(-1)) = K^(-1), -! alpha = 1 and beta = 0. -! -! 4. (Approximate) solution, through the LU or incomplete LU factorization, of -! a linear system -! A*Y = X, -! -! replicated on the processes. Here K = L*U = A or K = L*U ~ A, op(K^(-1)) = -! K^(-1), alpha = 1 and beta = 0. -! -! The block-Jacobi preconditioner or solver and the L and U factors of the LU -! or ILU factorizations have been built by the routine mld_fact_bld and stored -! into the 'base preconditioner' data structure prec. See mld_fact_bld for more -! details. -! -! This routine is used by mld_as_aply, to apply a 'base' block-Jacobi or -! Additive Schwarz (AS) preconditioner at any level of a multilevel preconditioner, -! or a block-Jacobi or LU or ILU solver at the coarsest level of a multilevel -! preconditioner. -! -! Tasks 1, 3 and 4 may be selected when prec%iprcparm(mld_smoother_sweeps_) = 1, -! while task 2 is selected when prec%iprcparm(mld_smoother_sweeps_) > 1. -! Furthermore, tasks 1, 2 and 3 may be performed when the matrix A is distributed -! among the processes (p%precv(ilev)%iprcparm(mld_coarse_mat_) = mld_distr_mat_, -! where p%precv(ilev) is the one-level data structure associated to the level -! ilev at which mld_sub_aply is called), while task 4 may be performed when A -! is replicated on the processes (p%precv(ilev)%iprcparm(mld_coarse_mat_) = -! mld_repl_mat_). Note that the matrix A is distributed among the processes -! at each level of the multilevel preconditioner, except the coarsest one, where -! it may be either distributed or replicated on the processes. Tasks 2, 3 and 4 -! are performed only at the coarsest level. Note also that this routine manages -! implicitly the fact that the matrix is distributed or replicated, i.e. it does not -! make any explicit reference to the value of p%precv(ilev)%iprcparm(mld_coarse_mat_). -! -! Arguments: -! -! alpha - complex(psb_spk_), input. -! The scalar alpha. -! prec - type(mld_cbaseprec_type), input. -! The 'base preconditioner' data structure containing the local -! part of the preconditioner or solver. -! x - complex(psb_spk_), dimension(:), input. -! The local part of the vector X. -! beta - complex(psb_spk_), input. -! The scalar beta. -! y - complex(psb_spk_), dimension(:), input/output. -! The local part of the vector Y. -! desc_data - type(psb_desc_type), input. -! The communication descriptor associated to the matrix to be -! preconditioned or 'inverted'. -! trans - character(len=1), input. -! If trans='N','n' then op(K^(-1)) = K^(-1); -! if trans='T','t' then op(K^(-1)) = K^(-T) (transpose of K^(-1)). -! if trans='C','c' then op(K^(-1)) = K^(-C) (transpose conjugate of K^(-1)). -! If prec%iprcparm(mld_smoother_sweeps_) > 1, the value of trans provided -! in input is ignored. -! work - complex(psb_spk_), dimension (:), target. -! Workspace. Its size must be at least 4*psb_cd_get_local_cols(desc_data). -! info - integer, output. -! Error code. -! -subroutine mld_csub_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) - - use psb_sparse_mod - use mld_inner_mod, mld_protect_name => mld_csub_aply - - implicit none - - ! Arguments - type(psb_desc_type), intent(in) :: desc_data - type(mld_cbaseprec_type), intent(in) :: prec - complex(psb_spk_),intent(in) :: x(:) - complex(psb_spk_),intent(inout) :: y(:) - complex(psb_spk_),intent(in) :: alpha,beta - character(len=1), intent(in) :: trans - complex(psb_spk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - ! Local variables - integer :: n_row,n_col - complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character(len=20) :: name - character :: trans_ - - name='mld_csub_aply' - info = psb_success_ - call psb_erractionsave(err_act) - - ictxt=psb_cd_get_context(desc_data) - call psb_info(ictxt, me, np) - - 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 - - - n_row = psb_cd_get_local_rows(desc_data) - n_col = psb_cd_get_local_cols(desc_data) - - if (n_col <= size(work)) then - ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then - aux => work(n_col+1:) - else - allocate(aux(4*n_col),stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/4*n_col,0,0,0,0/),& - & a_err='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 (prec%iprcparm(mld_smoother_sweeps_) == 1) then - - call mld_sub_solve(alpha,prec,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 (prec%iprcparm(mld_smoother_sweeps_) > 1) then - ! - ! - ! Apply prec%iprcparm(mld_smoother_sweeps_) sweeps of a block-Jacobi solver - ! to compute an approximate solution of a linear system. - ! - - if (size(prec%av) < mld_ap_nd_) then - info = psb_err_from_subroutine_non_ - goto 9999 - endif - - 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, prec%iprcparm(mld_smoother_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,prec%av(mld_ap_nd_),tx,cone,ty,& - & prec%desc_data,info,work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call mld_sub_solve(cone,prec,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,prec%iprcparm(mld_smoother_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.eq.psb_act_abort_) then - call psb_error() - return - end if - return - -end subroutine mld_csub_aply diff --git a/mlprec/mld_csub_solve.f90 b/mlprec/mld_csub_solve.f90 deleted file mode 100644 index 53590a8d..00000000 --- a/mlprec/mld_csub_solve.f90 +++ /dev/null @@ -1,324 +0,0 @@ -!!$ -!!$ -!!$ MLD2P4 version 2.0 -!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package -!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) -!!$ -!!$ (C) Copyright 2008,2009,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. -!!$ -!!$ -! File mld_csub_solve.f90 -! -! Subroutine: mld_csub_solve -! Version: complex -! -! This routine computes -! -! Y = beta*Y + alpha*op(K^(-1))*X, -! -! where -! - K is a factored matrix, as specified below, -! - op(K^(-1)) is K^(-1) or its transpose, according to the value of the -! argument trans, -! - X and Y are vectors, -! - alpha and beta are scalars. -! -! Depending on K, alpha and beta (and on the communication descriptor desc_data -! - see the arguments below), the above computation may correspond to one of -! the following tasks: -! -! 1. approximate solution of a linear system -! -! A*Y = X, -! -! by using the L and U factors computed with an ILU factorization of A. -! In this case K = L*U ~ A, alpha = 1 and beta = 0. The factors L and U -! (and the matrix A) are either distributed and block-diagonal or replicated. -! -! 2. Solution of a linear system -! -! A*Y = X, -! -! by using the L and U factors computed with a LU factorization of A. In this -! case K = L*U = A, alpha = 1 and beta = 0. The LU factorization is performed -! by one of the following auxiliary pakages: -! a. UMFPACK, -! b. SuperLU, -! c. SuperLU_Dist. -! In the cases a. and b., the factors L and U (and the matrix A) are either -! distributed and block diagonal) or replicated; in the case c., L, U (and A) -! are distributed. -! -! This routine is used by mld_dsub_aply, to apply a 'base' block-Jacobi or -! Additive Schwarz (AS) preconditioner at any level of a multilevel preconditioner, -! or a block-Jacobi or LU or ILU solver at the coarsest level of a multilevel -! preconditioner. -! -! -! Arguments: -! -! alpha - complex(psb_spk_), input. -! The scalar alpha. -! prec - type(mld_cbaseprec_type), input. -! The 'base preconditioner' data structure containing the local -! part of the L and U factors of the matrix A. -! x - complex(psb_spk_), dimension(:), input. -! The local part of the vector X. -! beta - complex(psb_spk_), input. -! The scalar beta. -! y - complex(psb_spk_), dimension(:), input/output. -! The local part of the vector Y. -! desc_data - type(psb_desc_type), input. -! The communication descriptor associated to the matrix to be -! preconditioned or 'inverted'. -! trans - character(len=1), input. -! If trans='N','n' then op(K^(-1)) = K^(-1); -! if trans='T','t' then op(K^(-1)) = K^(-T) (transpose of K^(-1)). -! if trans='C','c' then op(K^(-1)) = K^(-C) (transpose conjugate of K^(-1)). -! If prec%iprcparm(mld_smoother_sweeps_) > 1, the value of trans provided -! in input is ignored. -! work - complex(psb_spk_), dimension (:), target. -! Workspace. Its size must be at least 4*psb_cd_get_local_cols(desc_data). -! info - integer, output. -! Error code. -! -subroutine mld_csub_solve(alpha,prec,x,beta,y,desc_data,trans,work,info) - - use psb_sparse_mod - use mld_inner_mod, mld_protect_name => mld_csub_solve - - implicit none - - ! Arguments - type(psb_desc_type), intent(in) :: desc_data - type(mld_cbaseprec_type), intent(in) :: prec - complex(psb_spk_),intent(in) :: x(:) - complex(psb_spk_),intent(inout) :: y(:) - complex(psb_spk_),intent(in) :: alpha,beta - character(len=1), intent(in) :: trans - complex(psb_spk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - ! Local variables - integer :: n_row,n_col - complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character(len=20) :: name - character :: trans_ - - interface - subroutine mld_cumf_solve(flag,m,x,b,n,ptr,info) - use psb_sparse_mod - integer, intent(in) :: flag,m,n,ptr - integer, intent(out) :: info - complex(psb_spk_), intent(in) :: b(*) - complex(psb_spk_), intent(inout) :: x(*) - end subroutine mld_cumf_solve - end interface - - name='mld_csub_solve' - info = psb_success_ - call psb_erractionsave(err_act) - - ictxt=psb_cd_get_context(desc_data) - call psb_info(ictxt, me, np) - - 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 - - - n_row = psb_cd_get_local_rows(desc_data) - n_col = psb_cd_get_local_cols(desc_data) - - if (n_col <= size(work)) then - ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then - aux => work(n_col+1:) - else - allocate(aux(4*n_col),stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/4*n_col,0,0,0,0/),& - & a_err='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 - - - select case(prec%iprcparm(mld_sub_solve_)) - case(mld_ilu_n_,mld_milu_n_,mld_ilu_t_) - ! - ! Apply a block-Jacobi preconditioner with ILU(k)/MILU(k)/ILU(k,t) - ! factorization of the blocks (distributed matrix) or approximately - ! solve a system through ILU(k)/MILU(k)/ILU(k,t) (replicated matrix). - ! - - select case(trans_) - case('N') - - call psb_spsm(cone,prec%av(mld_l_pr_),x,czero,ww,desc_data,info,& - & trans=trans_,unit='L',diag=prec%d,choice=psb_none_,work=aux) - if (info == psb_success_) call psb_spsm(alpha,prec%av(mld_u_pr_),ww,beta,y,desc_data,info,& - & trans=trans_,unit='U',choice=psb_none_, work=aux) - - case('T') - call psb_spsm(cone,prec%av(mld_u_pr_),x,czero,ww,desc_data,info,& - & trans=trans_,unit='L',diag=prec%d,choice=psb_none_, work=aux) - if(info == psb_success_) call psb_spsm(alpha,prec%av(mld_l_pr_),ww,beta,y,desc_data,info,& - & trans=trans_,unit='U',choice=psb_none_,work=aux) - - case('C') - call psb_spsm(cone,prec%av(mld_u_pr_),x,czero,ww,desc_data,info,& - & trans=trans_,unit='L',diag=conjg(prec%d),choice=psb_none_, work=aux) - if(info == psb_success_) call psb_spsm(alpha,prec%av(mld_l_pr_),ww,beta,y,desc_data,info,& - & trans=trans_,unit='U',choice=psb_none_,work=aux) - case default - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in ILU subsolve') - goto 9999 - end select - - case(mld_slu_) - ! - ! Apply a block-Jacobi preconditioner with LU factorization of the - ! blocks (distributed matrix) or approximately solve a local linear - ! system through LU (replicated matrix). The SuperLU package is used - ! to apply the LU factorization in both cases. - ! - - ww(1:n_row) = x(1:n_row) - - select case(trans_) - case('N') - call mld_cslu_solve(0,n_row,1,ww,n_row,prec%iprcparm(mld_slu_ptr_),info) - case('T') - call mld_cslu_solve(1,n_row,1,ww,n_row,prec%iprcparm(mld_slu_ptr_),info) - case('C') - call mld_cslu_solve(2,n_row,1,ww,n_row,prec%iprcparm(mld_slu_ptr_),info) - case default - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in SLU subsolve') - goto 9999 - end select - - if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info) - - case(mld_sludist_) - ! - ! Solve a distributed linear system with the LU factorization. - ! The SuperLU_DIST package is used. - ! - - ww(1:n_row) = x(1:n_row) - - select case(trans_) - case('N') - call mld_csludist_solve(0,n_row,1,ww,n_row,prec%iprcparm(mld_slud_ptr_),info) - case('T') - call mld_csludist_solve(1,n_row,1,ww,n_row,prec%iprcparm(mld_slud_ptr_),info) - case('C') - call mld_csludist_solve(2,n_row,1,ww,n_row,prec%iprcparm(mld_slud_ptr_),info) - case default - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in SLUDist subsolve') - goto 9999 - end select - - if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info) - - case (mld_umf_) - ! - ! Apply a block-Jacobi preconditioner with LU factorization of the - ! blocks (distributed matrix) or approximately solve a local linear - ! system through LU (replicated matrix). The UMFPACK package is used - ! to apply the LU factorization in both cases. - ! - - select case(trans_) - case('N') - call mld_cumf_solve(0,n_row,ww,x,n_row,prec%iprcparm(mld_umf_numptr_),info) - case('T') - call mld_cumf_solve(1,n_row,ww,x,n_row,prec%iprcparm(mld_umf_numptr_),info) - case('C') - call mld_cumf_solve(2,n_row,ww,x,n_row,prec%iprcparm(mld_umf_numptr_),info) - case default - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in UMF subsolve') - goto 9999 - end select - - if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info) - - case default - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid mld_sub_solve_') - goto 9999 - - end select - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error in subsolve ') - goto 9999 - endif - - if (n_col <= size(work)) then - if ((4*n_col+n_col) <= size(work)) then - else - deallocate(aux) - endif - else - deallocate(ww,aux) - endif - - call psb_erractionrestore(err_act) - return - -9999 continue - - call psb_erractionrestore(err_act) - if (err_act.eq.psb_act_abort_) then - call psb_error() - return - end if - return - -end subroutine mld_csub_solve diff --git a/mlprec/mld_dsub_aply.f90 b/mlprec/mld_dsub_aply.f90 deleted file mode 100644 index 2b7ce3cf..00000000 --- a/mlprec/mld_dsub_aply.f90 +++ /dev/null @@ -1,296 +0,0 @@ -!!$ -!!$ -!!$ MLD2P4 version 2.0 -!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package -!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) -!!$ -!!$ (C) Copyright 2008,2009,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. -!!$ -!!$ -! File mld_dsub_aply.f90 -! -! Subroutine: mld_dsub_aply -! Version: real -! -! This routine computes -! -! Y = beta*Y + alpha*op(K^(-1))*X, -! -! where -! - K is a suitable matrix, as specified below, -! - op(K^(-1)) is K^(-1) or its transpose, according to the value of the -! argument trans, -! - X and Y are vectors, -! - alpha and beta are scalars. -! -! Depending on K, alpha and beta (and on the communication descriptor desc_data -! - see the arguments below), the above computation may correspond to one of -! the following tasks: -! -! 1. Application of a block-Jacobi preconditioner associated to a matrix A -! distributed among the processes. Here K is the preconditioner, op(K^(-1)) -! = K^(-1), alpha = 1 and beta = 0. -! -! 2. Application of block-Jacobi sweeps to compute an approximate solution of -! a linear system -! A*Y = X, -! -! distributed among the processes (note that a single block-Jacobi sweep, -! with null starting guess, corresponds to the application of a block-Jacobi -! preconditioner). Here K^(-1) denotes the iteration matrix of the -! block-Jacobi solver, op(K^(-1)) = K^(-1), alpha = 1 and beta = 0. -! -! 3. Solution, through the LU factorization, of a linear system -! -! A*Y = X, -! -! distributed among the processes. Here K = L*U = A, op(K^(-1)) = K^(-1), -! alpha = 1 and beta = 0. -! -! 4. (Approximate) solution, through the LU or incomplete LU factorization, of -! a linear system -! A*Y = X, -! -! replicated on the processes. Here K = L*U = A or K = L*U ~ A, op(K^(-1)) = -! K^(-1), alpha = 1 and beta = 0. -! -! The block-Jacobi preconditioner or solver and the L and U factors of the LU -! or ILU factorizations have been built by the routine mld_fact_bld and stored -! into the 'base preconditioner' data structure prec. See mld_fact_bld for more -! details. -! -! This routine is used by mld_as_aply, to apply a 'base' block-Jacobi or -! Additive Schwarz (AS) preconditioner at any level of a multilevel preconditioner, -! or a block-Jacobi or LU or ILU solver at the coarsest level of a multilevel -! preconditioner. -! -! Tasks 1, 3 and 4 may be selected when prec%iprcparm(mld_smoother_sweeps_) = 1, -! while task 2 is selected when prec%iprcparm(mld_smoother_sweeps_) > 1. -! Furthermore, tasks 1, 2 and 3 may be performed when the matrix A is distributed -! among the processes (p%precv(ilev)%iprcparm(mld_coarse_mat_) = mld_distr_mat_, -! where p%precv(ilev) is the one-level data structure associated to the level -! ilev at which mld_sub_aply is called), while task 4 may be performed when A -! is replicated on the processes (p%precv(ilev)%iprcparm(mld_coarse_mat_) = -! mld_repl_mat_). Note that the matrix A is distributed among the processes -! at each level of the multilevel preconditioner, except the coarsest one, where -! it may be either distributed or replicated on the processes. Tasks 2, 3 and 4 -! are performed only at the coarsest level. Note also that this routine manages -! implicitly the fact that the matrix is distributed or replicated, i.e. it does not -! make any explicit reference to the value of p%precv(ilev)%iprcparm(mld_coarse_mat_). -! -! Arguments: -! -! alpha - real(psb_dpk_), input. -! The scalar alpha. -! prec - type(mld_dbaseprec_type), input. -! The 'base preconditioner' data structure containing the local -! part of the preconditioner or solver. -! x - real(psb_dpk_), dimension(:), input. -! The local part of the vector X. -! beta - real(psb_dpk_), input. -! The scalar beta. -! y - real(psb_dpk_), dimension(:), input/output. -! The local part of the vector Y. -! desc_data - type(psb_desc_type), input. -! The communication descriptor associated to the matrix to be -! preconditioned or 'inverted'. -! trans - character(len=1), input. -! If trans='N','n' then op(K^(-1)) = K^(-1); -! if trans='T','t' then op(K^(-1)) = K^(-T) (transpose of K^(-1)). -! If prec%iprcparm(mld_smoother_sweeps_) > 1, the value of trans provided -! in input is ignored. -! work - real(psb_dpk_), dimension (:), target. -! Workspace. Its size must be at least 4*psb_cd_get_local_cols(desc_data). -! info - integer, output. -! Error code. -! -subroutine mld_dsub_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) - - use psb_sparse_mod - use mld_inner_mod, mld_protect_name => mld_dsub_aply - - implicit none - - ! Arguments - type(psb_desc_type), intent(in) :: desc_data - type(mld_dbaseprec_type), intent(in) :: prec - real(psb_dpk_),intent(in) :: x(:) - real(psb_dpk_),intent(inout) :: y(:) - real(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - real(psb_dpk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - ! Local variables - integer :: n_row,n_col - real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character(len=20) :: name - character :: trans_ - - name='mld_dsub_aply' - info = psb_success_ - call psb_erractionsave(err_act) - - ictxt=psb_cd_get_context(desc_data) - call psb_info(ictxt, me, np) - - 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 - - - n_row = psb_cd_get_local_rows(desc_data) - n_col = psb_cd_get_local_cols(desc_data) - - if (n_col <= size(work)) then - ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then - aux => work(n_col+1:) - else - allocate(aux(4*n_col),stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/4*n_col,0,0,0,0/),& - & a_err='real(psb_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 (prec%iprcparm(mld_smoother_sweeps_) == 1) then - - call mld_sub_solve(alpha,prec,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 (prec%iprcparm(mld_smoother_sweeps_) > 1) then - ! - ! - ! Apply prec%iprcparm(mld_smoother_sweeps_) sweeps of a block-Jacobi solver - ! to compute an approximate solution of a linear system. - ! - ! - - if (size(prec%av) < mld_ap_nd_) then - info = psb_err_from_subroutine_non_ - goto 9999 - endif - - 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, prec%iprcparm(mld_smoother_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,prec%av(mld_ap_nd_),tx,done,ty,& - & prec%desc_data,info,work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call mld_sub_solve(done,prec,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,prec%iprcparm(mld_smoother_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.eq.psb_act_abort_) then - call psb_error() - return - end if - return - -end subroutine mld_dsub_aply diff --git a/mlprec/mld_dsub_solve.f90 b/mlprec/mld_dsub_solve.f90 deleted file mode 100644 index 778e79c0..00000000 --- a/mlprec/mld_dsub_solve.f90 +++ /dev/null @@ -1,312 +0,0 @@ -!!$ -!!$ -!!$ MLD2P4 version 2.0 -!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package -!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) -!!$ -!!$ (C) Copyright 2008,2009,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. -!!$ -!!$ -! File mld_dsub_solve.f90 -! -! Subroutine: mld_dsub_solve -! Version: real -! -! This routine computes -! -! Y = beta*Y + alpha*op(K^(-1))*X, -! -! where -! - K is a factored matrix, as specified below, -! - op(K^(-1)) is K^(-1) or its transpose, according to the value of the -! argument trans, -! - X and Y are vectors, -! - alpha and beta are scalars. -! -! Depending on K, alpha and beta (and on the communication descriptor desc_data -! - see the arguments below), the above computation may correspond to one of -! the following tasks: -! -! 1. approximate solution of a linear system -! -! A*Y = X, -! -! by using the L and U factors computed with an ILU factorization of A. -! In this case K = L*U ~ A, alpha = 1 and beta = 0. The factors L and U -! (and the matrix A) are either distributed and block-diagonal or replicated. -! -! 2. Solution of a linear system -! -! A*Y = X, -! -! by using the L and U factors computed with a LU factorization of A. In this -! case K = L*U = A, alpha = 1 and beta = 0. The LU factorization is performed -! by one of the following auxiliary pakages: -! a. UMFPACK, -! b. SuperLU, -! c. SuperLU_Dist. -! In the cases a. and b., the factors L and U (and the matrix A) are either -! distributed and block diagonal) or replicated; in the case c., L, U (and A) -! are distributed. -! -! This routine is used by mld_dsub_aply, to apply a 'base' block-Jacobi or -! Additive Schwarz (AS) preconditioner at any level of a multilevel preconditioner, -! or a block-Jacobi or LU or ILU solver at the coarsest level of a multilevel -! preconditioner. -! -! -! Arguments: -! -! alpha - real(psb_dpk_), input. -! The scalar alpha. -! prec - type(mld_dbaseprec_type), input. -! The 'base preconditioner' data structure containing the local -! part of the L and U factors of the matrix A. -! x - real(psb_dpk_), dimension(:), input. -! The local part of the vector X. -! beta - real(psb_dpk_), input. -! The scalar beta. -! y - real(psb_dpk_), dimension(:), input/output. -! The local part of the vector Y. -! desc_data - type(psb_desc_type), input. -! The communication descriptor associated to the matrix to be -! preconditioned or 'inverted'. -! trans - character(len=1), input. -! If trans='N','n' then op(K^(-1)) = K^(-1); -! if trans='T','t' then op(K^(-1)) = K^(-T) (transpose of K^(-1)). -! If prec%iprcparm(mld_smoother_sweeps_) > 1, the value of trans provided -! in input is ignored. -! work - real(psb_dpk_), dimension (:), target. -! Workspace. Its size must be at least 4*psb_cd_get_local_cols(desc_data). -! info - integer, output. -! Error code. -! -subroutine mld_dsub_solve(alpha,prec,x,beta,y,desc_data,trans,work,info) - - use psb_sparse_mod - use mld_inner_mod, mld_protect_name => mld_dsub_solve - - implicit none - - ! Arguments - type(psb_desc_type), intent(in) :: desc_data - type(mld_dbaseprec_type), intent(in) :: prec - real(psb_dpk_),intent(in) :: x(:) - real(psb_dpk_),intent(inout) :: y(:) - real(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - real(psb_dpk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - ! Local variables - integer :: n_row,n_col - real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character(len=20) :: name - character :: trans_ - - interface - subroutine mld_dumf_solve(flag,m,x,b,n,ptr,info) - use psb_sparse_mod - integer, intent(in) :: flag,m,n,ptr - integer, intent(out) :: info - real(psb_dpk_), intent(in) :: b(*) - real(psb_dpk_), intent(inout) :: x(*) - end subroutine mld_dumf_solve - end interface - - name='mld_dsub_solve' - info = psb_success_ - call psb_erractionsave(err_act) - - ictxt=psb_cd_get_context(desc_data) - call psb_info(ictxt, me, np) - - 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 - - - n_row = psb_cd_get_local_rows(desc_data) - n_col = psb_cd_get_local_cols(desc_data) - - if (n_col <= size(work)) then - ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then - aux => work(n_col+1:) - else - allocate(aux(4*n_col),stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/4*n_col,0,0,0,0/),& - & a_err='real(psb_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 - - - select case(prec%iprcparm(mld_sub_solve_)) - case(mld_ilu_n_,mld_milu_n_,mld_ilu_t_) - ! - ! Apply a block-Jacobi preconditioner with ILU(k)/MILU(k)/ILU(k,t) - ! factorization of the blocks (distributed matrix) or approximately - ! solve a system through ILU(k)/MILU(k)/ILU(k,t) (replicated matrix). - ! - - select case(trans_) - case('N') - call psb_spsm(done,prec%av(mld_l_pr_),x,dzero,ww,desc_data,info,& - & trans=trans_,scale='L',diag=prec%d,choice=psb_none_,work=aux) - - if (info == psb_success_) call psb_spsm(alpha,prec%av(mld_u_pr_),ww,beta,y,desc_data,info,& - & trans=trans_,scale='U',choice=psb_none_, work=aux) - - case('T','C') - call psb_spsm(done,prec%av(mld_u_pr_),x,dzero,ww,desc_data,info,& - & trans=trans_,scale='L',diag=prec%d,choice=psb_none_,work=aux) - if (info == psb_success_) call psb_spsm(alpha,prec%av(mld_l_pr_),ww,beta,y,desc_data,info,& - & trans=trans_,scale='U',choice=psb_none_,work=aux) - case default - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in ILU subsolve') - goto 9999 - end select - - case(mld_slu_) - ! - ! Apply a block-Jacobi preconditioner with LU factorization of the - ! blocks (distributed matrix) or approximately solve a local linear - ! system through LU (replicated matrix). The SuperLU package is used - ! to apply the LU factorization in both cases. - ! - - ww(1:n_row) = x(1:n_row) - - select case(trans_) - case('N') - call mld_dslu_solve(0,n_row,1,ww,n_row,prec%iprcparm(mld_slu_ptr_),info) - case('T','C') - call mld_dslu_solve(1,n_row,1,ww,n_row,prec%iprcparm(mld_slu_ptr_),info) - case default - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in SLU subsolve') - goto 9999 - end select - - if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info) - - case(mld_sludist_) - ! - ! Solve a distributed linear system with the LU factorization. - ! The SuperLU_DIST package is used. - ! - - ww(1:n_row) = x(1:n_row) - - select case(trans_) - case('N') - call mld_dsludist_solve(0,n_row,1,ww,n_row,prec%iprcparm(mld_slud_ptr_),info) - case('T','C') - call mld_dsludist_solve(1,n_row,1,ww,n_row,prec%iprcparm(mld_slud_ptr_),info) - case default - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in SLUDist subsolve') - goto 9999 - end select - - if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info) - - case (mld_umf_) - ! - ! Apply a block-Jacobi preconditioner with LU factorization of the - ! blocks (distributed matrix) or approximately solve a local linear - ! system through LU (replicated matrix). The UMFPACK package is used - ! to apply the LU factorization in both cases. - ! - - select case(trans_) - case('N') - call mld_dumf_solve(0,n_row,ww,x,n_row,prec%iprcparm(mld_umf_numptr_),info) - case('T','C') - call mld_dumf_solve(1,n_row,ww,x,n_row,prec%iprcparm(mld_umf_numptr_),info) - case default - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in UMF subsolve') - goto 9999 - end select - - if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info) - - case default - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid mld_sub_solve_') - goto 9999 - - end select - - if (info /= psb_success_) then - - call psb_errpush(psb_err_internal_error_,name,a_err='Error in subsolve') - goto 9999 - endif - - if (n_col <= size(work)) then - if ((4*n_col+n_col) <= size(work)) then - else - deallocate(aux) - endif - else - deallocate(ww,aux) - endif - - call psb_erractionrestore(err_act) - return - -9999 continue - - call psb_erractionrestore(err_act) - if (err_act.eq.psb_act_abort_) then - call psb_error() - return - end if - return - -end subroutine mld_dsub_solve diff --git a/mlprec/mld_inner_mod.f90 b/mlprec/mld_inner_mod.f90 index ad505949..630ded4a 100644 --- a/mlprec/mld_inner_mod.f90 +++ b/mlprec/mld_inner_mod.f90 @@ -245,110 +245,6 @@ module mld_inner_mod end interface - interface mld_sub_aply - subroutine mld_ssub_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) - use psb_sparse_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_ - use mld_prec_type, only : mld_sbaseprec_type - type(psb_desc_type), intent(in) :: desc_data - type(mld_sbaseprec_type), intent(in) :: prec - real(psb_spk_),intent(in) :: x(:) - real(psb_spk_),intent(inout) :: y(:) - real(psb_spk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - real(psb_spk_),target,intent(inout) :: work(:) - integer, intent(out) :: info - end subroutine mld_ssub_aply - subroutine mld_dsub_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) - use psb_sparse_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_ - use mld_prec_type, only : mld_dbaseprec_type - type(psb_desc_type), intent(in) :: desc_data - type(mld_dbaseprec_type), intent(in) :: prec - real(psb_dpk_),intent(in) :: x(:) - real(psb_dpk_),intent(inout) :: y(:) - real(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - real(psb_dpk_),target,intent(inout) :: work(:) - integer, intent(out) :: info - end subroutine mld_dsub_aply - subroutine mld_csub_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) - use psb_sparse_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_ - use mld_prec_type, only : mld_cbaseprec_type - type(psb_desc_type), intent(in) :: desc_data - type(mld_cbaseprec_type), intent(in) :: prec - complex(psb_spk_),intent(in) :: x(:) - complex(psb_spk_),intent(inout) :: y(:) - complex(psb_spk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - complex(psb_spk_),target,intent(inout) :: work(:) - integer, intent(out) :: info - end subroutine mld_csub_aply - subroutine mld_zsub_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) - use psb_sparse_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_ - use mld_prec_type, only : mld_zbaseprec_type - type(psb_desc_type), intent(in) :: desc_data - type(mld_zbaseprec_type), intent(in) :: prec - complex(psb_dpk_),intent(in) :: x(:) - complex(psb_dpk_),intent(inout) :: y(:) - complex(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - complex(psb_dpk_),target,intent(inout) :: work(:) - integer, intent(out) :: info - end subroutine mld_zsub_aply - end interface - - - interface mld_sub_solve - subroutine mld_ssub_solve(alpha,prec,x,beta,y,desc_data,trans,work,info) - use psb_sparse_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_ - use mld_prec_type, only : mld_sbaseprec_type - type(psb_desc_type), intent(in) :: desc_data - type(mld_sbaseprec_type), intent(in) :: prec - real(psb_spk_),intent(in) :: x(:) - real(psb_spk_),intent(inout) :: y(:) - real(psb_spk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - real(psb_spk_),target,intent(inout) :: work(:) - integer, intent(out) :: info - end subroutine mld_ssub_solve - subroutine mld_dsub_solve(alpha,prec,x,beta,y,desc_data,trans,work,info) - use psb_sparse_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_ - use mld_prec_type, only : mld_dbaseprec_type - type(psb_desc_type), intent(in) :: desc_data - type(mld_dbaseprec_type), intent(in) :: prec - real(psb_dpk_),intent(in) :: x(:) - real(psb_dpk_),intent(inout) :: y(:) - real(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - real(psb_dpk_),target,intent(inout) :: work(:) - integer, intent(out) :: info - end subroutine mld_dsub_solve - subroutine mld_csub_solve(alpha,prec,x,beta,y,desc_data,trans,work,info) - use psb_sparse_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_ - use mld_prec_type, only : mld_cbaseprec_type - type(psb_desc_type), intent(in) :: desc_data - type(mld_cbaseprec_type), intent(in) :: prec - complex(psb_spk_),intent(in) :: x(:) - complex(psb_spk_),intent(inout) :: y(:) - complex(psb_spk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - complex(psb_spk_),target,intent(inout) :: work(:) - integer, intent(out) :: info - end subroutine mld_csub_solve - subroutine mld_zsub_solve(alpha,prec,x,beta,y,desc_data,trans,work,info) - use psb_sparse_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_ - use mld_prec_type, only : mld_zbaseprec_type - type(psb_desc_type), intent(in) :: desc_data - type(mld_zbaseprec_type), intent(in) :: prec - complex(psb_dpk_),intent(in) :: x(:) - complex(psb_dpk_),intent(inout) :: y(:) - complex(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - complex(psb_dpk_),target,intent(inout) :: work(:) - integer, intent(out) :: info - end subroutine mld_zsub_solve - end interface - - interface mld_asmat_bld Subroutine mld_sasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) use psb_sparse_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_ diff --git a/mlprec/mld_ssub_aply.f90 b/mlprec/mld_ssub_aply.f90 deleted file mode 100644 index 934ae233..00000000 --- a/mlprec/mld_ssub_aply.f90 +++ /dev/null @@ -1,297 +0,0 @@ -!!$ -!!$ -!!$ MLD2P4 version 2.0 -!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package -!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) -!!$ -!!$ (C) Copyright 2008,2009,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. -!!$ -!!$ -! File mld_ssub_aply.f90 -! -! Subroutine: mld_ssub_aply -! Version: real -! -! This routine computes -! -! Y = beta*Y + alpha*op(K^(-1))*X, -! -! where -! - K is a suitable matrix, as specified below, -! - op(K^(-1)) is K^(-1) or its transpose, according to the value of the -! argument trans, -! - X and Y are vectors, -! - alpha and beta are scalars. -! -! Depending on K, alpha and beta (and on the communication descriptor desc_data -! - see the arguments below), the above computation may correspond to one of -! the following tasks: -! -! 1. Application of a block-Jacobi preconditioner associated to a matrix A -! distributed among the processes. Here K is the preconditioner, op(K^(-1)) -! = K^(-1), alpha = 1 and beta = 0. -! -! 2. Application of block-Jacobi sweeps to compute an approximate solution of -! a linear system -! A*Y = X, -! -! distributed among the processes (note that a single block-Jacobi sweep, -! with null starting guess, corresponds to the application of a block-Jacobi -! preconditioner). Here K^(-1) denotes the iteration matrix of the -! block-Jacobi solver, op(K^(-1)) = K^(-1), alpha = 1 and beta = 0. -! -! 3. Solution, through the LU factorization, of a linear system -! -! A*Y = X, -! -! distributed among the processes. Here K = L*U = A, op(K^(-1)) = K^(-1), -! alpha = 1 and beta = 0. -! -! 4. (Approximate) solution, through the LU or incomplete LU factorization, of -! a linear system -! A*Y = X, -! -! replicated on the processes. Here K = L*U = A or K = L*U ~ A, op(K^(-1)) = -! K^(-1), alpha = 1 and beta = 0. -! -! The block-Jacobi preconditioner or solver and the L and U factors of the LU -! or ILU factorizations have been built by the routine mld_fact_bld and stored -! into the 'base preconditioner' data structure prec. See mld_fact_bld for more -! details. -! -! This routine is used by mld_as_aply, to apply a 'base' block-Jacobi or -! Additive Schwarz (AS) preconditioner at any level of a multilevel preconditioner, -! or a block-Jacobi or LU or ILU solver at the coarsest level of a multilevel -! preconditioner. -! -! Tasks 1, 3 and 4 may be selected when prec%iprcparm(mld_smoother_sweeps_) = 1, -! while task 2 is selected when prec%iprcparm(mld_smoother_sweeps_) > 1. -! Furthermore, tasks 1, 2 and 3 may be performed when the matrix A is distributed -! among the processes (p%precv(ilev)%iprcparm(mld_coarse_mat_) = mld_distr_mat_, -! where p%precv(ilev) is the one-level data structure associated to the level -! ilev at which mld_sub_aply is called), while task 4 may be performed when A -! is replicated on the processes (p%precv(ilev)%iprcparm(mld_coarse_mat_) = -! mld_repl_mat_). Note that the matrix A is distributed among the processes -! at each level of the multilevel preconditioner, except the coarsest one, where -! it may be either distributed or replicated on the processes. Tasks 2, 3 and 4 -! are performed only at the coarsest level. Note also that this routine manages -! implicitly the fact that the matrix is distributed or replicated, i.e. it does not -! make any explicit reference to the value of p%precv(ilev)%iprcparm(mld_coarse_mat_). -! -! Arguments: -! -! alpha - real(psb_spk_), input. -! The scalar alpha. -! prec - type(mld_sbaseprec_type), input. -! The 'base preconditioner' data structure containing the local -! part of the preconditioner or solver. -! x - real(psb_spk_), dimension(:), input. -! The local part of the vector X. -! beta - real(psb_spk_), input. -! The scalar beta. -! y - real(psb_spk_), dimension(:), input/output. -! The local part of the vector Y. -! desc_data - type(psb_desc_type), input. -! The communication descriptor associated to the matrix to be -! preconditioned or 'inverted'. -! trans - character(len=1), input. -! If trans='N','n' then op(K^(-1)) = K^(-1); -! if trans='T','t' then op(K^(-1)) = K^(-T) (transpose of K^(-1)). -! If prec%iprcparm(mld_smoother_sweeps_) > 1, the value of trans provided -! in input is ignored. -! work - real(psb_spk_), dimension (:), target. -! Workspace. Its size must be at least 4*psb_cd_get_local_cols(desc_data). -! info - integer, output. -! Error code. -! -subroutine mld_ssub_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) - - use psb_sparse_mod - use mld_inner_mod, mld_protect_name => mld_ssub_aply - - implicit none - - ! Arguments - type(psb_desc_type), intent(in) :: desc_data - type(mld_sbaseprec_type), intent(in) :: prec - real(psb_spk_),intent(in) :: x(:) - real(psb_spk_),intent(inout) :: y(:) - real(psb_spk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - real(psb_spk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - ! Local variables - integer :: n_row,n_col - real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character(len=20) :: name - character :: trans_ - - name='mld_ssub_aply' - info = psb_success_ - call psb_erractionsave(err_act) - - ictxt=psb_cd_get_context(desc_data) - call psb_info(ictxt, me, np) - - 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 - - - n_row = psb_cd_get_local_rows(desc_data) - n_col = psb_cd_get_local_cols(desc_data) - - if (n_col <= size(work)) then - ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then - aux => work(n_col+1:) - else - allocate(aux(4*n_col),stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/4*n_col,0,0,0,0/),& - & a_err='real(psb_spk_)') - goto 9999 - end if - endif - else - allocate(ww(n_col),aux(4*n_col),stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/5*n_col,0,0,0,0/),& - & a_err='real(psb_spk_)') - goto 9999 - end if - endif - - if (prec%iprcparm(mld_smoother_sweeps_) == 1) then - - call mld_sub_solve(alpha,prec,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 (prec%iprcparm(mld_smoother_sweeps_) > 1) then - ! - ! - ! Apply prec%iprcparm(mld_smoother_sweeps_) sweeps of a block-Jacobi solver - ! to compute an approximate solution of a linear system. - ! - ! - - if (size(prec%av) < mld_ap_nd_) then - info = psb_err_from_subroutine_non_ - goto 9999 - endif - - 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, prec%iprcparm(mld_smoother_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,prec%av(mld_ap_nd_),tx,sone,ty,& - & prec%desc_data,info,work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call mld_sub_solve(sone,prec,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,prec%iprcparm(mld_smoother_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.eq.psb_act_abort_) then - call psb_error() - return - end if - return - -end subroutine mld_ssub_aply - diff --git a/mlprec/mld_ssub_solve.f90 b/mlprec/mld_ssub_solve.f90 deleted file mode 100644 index ae6f2c0a..00000000 --- a/mlprec/mld_ssub_solve.f90 +++ /dev/null @@ -1,312 +0,0 @@ -!!$ -!!$ -!!$ MLD2P4 version 2.0 -!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package -!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) -!!$ -!!$ (C) Copyright 2008,2009,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. -!!$ -!!$ -! File mld_ssub_solve.f90 -! -! Subroutine: mld_ssub_solve -! Version: real -! -! This routine computes -! -! Y = beta*Y + alpha*op(K^(-1))*X, -! -! where -! - K is a factored matrix, as specified below, -! - op(K^(-1)) is K^(-1) or its transpose, according to the value of the -! argument trans, -! - X and Y are vectors, -! - alpha and beta are scalars. -! -! Depending on K, alpha and beta (and on the communication descriptor desc_data -! - see the arguments below), the above computation may correspond to one of -! the following tasks: -! -! 1. approximate solution of a linear system -! -! A*Y = X, -! -! by using the L and U factors computed with an ILU factorization of A. -! In this case K = L*U ~ A, alpha = 1 and beta = 0. The factors L and U -! (and the matrix A) are either distributed and block-diagonal or replicated. -! -! 2. Solution of a linear system -! -! A*Y = X, -! -! by using the L and U factors computed with a LU factorization of A. In this -! case K = L*U = A, alpha = 1 and beta = 0. The LU factorization is performed -! by one of the following auxiliary pakages: -! a. UMFPACK, -! b. SuperLU, -! c. SuperLU_Dist. -! In the cases a. and b., the factors L and U (and the matrix A) are either -! distributed and block diagonal) or replicated; in the case c., L, U (and A) -! are distributed. -! -! This routine is used by mld_ssub_aply, to apply a 'base' block-Jacobi or -! Additive Schwarz (AS) preconditioner at any level of a multilevel preconditioner, -! or a block-Jacobi or LU or ILU solver at the coarsest level of a multilevel -! preconditioner. -! -! -! Arguments: -! -! alpha - real(psb_spk_), input. -! The scalar alpha. -! prec - type(mld_sbaseprec_type), input. -! The 'base preconditioner' data structure containing the local -! part of the L and U factors of the matrix A. -! x - real(psb_spk_), dimension(:), input. -! The local part of the vector X. -! beta - real(psb_spk_), input. -! The scalar beta. -! y - real(psb_spk_), dimension(:), input/output. -! The local part of the vector Y. -! desc_data - type(psb_desc_type), input. -! The communication descriptor associated to the matrix to be -! preconditioned or 'inverted'. -! trans - character(len=1), input. -! If trans='N','n' then op(K^(-1)) = K^(-1); -! if trans='T','t' then op(K^(-1)) = K^(-T) (transpose of K^(-1)). -! If prec%iprcparm(mld_smoother_sweeps_) > 1, the value of trans provided -! in input is ignored. -! work - real(psb_spk_), dimension (:), target. -! Workspace. Its size must be at least 4*psb_cd_get_local_cols(desc_data). -! info - integer, output. -! Error code. -! -subroutine mld_ssub_solve(alpha,prec,x,beta,y,desc_data,trans,work,info) - - use psb_sparse_mod - use mld_inner_mod, mld_protect_name => mld_ssub_solve - - implicit none - - ! Arguments - type(psb_desc_type), intent(in) :: desc_data - type(mld_sbaseprec_type), intent(in) :: prec - real(psb_spk_),intent(in) :: x(:) - real(psb_spk_),intent(inout) :: y(:) - real(psb_spk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - real(psb_spk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - ! Local variables - integer :: n_row,n_col - real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character(len=20) :: name - character :: trans_ - - interface - subroutine mld_sumf_solve(flag,m,x,b,n,ptr,info) - use psb_sparse_mod - integer, intent(in) :: flag,m,n,ptr - integer, intent(out) :: info - real(psb_spk_), intent(in) :: b(*) - real(psb_spk_), intent(inout) :: x(*) - end subroutine mld_sumf_solve - end interface - - name='mld_ssub_solve' - info = psb_success_ - call psb_erractionsave(err_act) - - ictxt=psb_cd_get_context(desc_data) - call psb_info(ictxt, me, np) - - 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 - - - n_row = psb_cd_get_local_rows(desc_data) - n_col = psb_cd_get_local_cols(desc_data) - - if (n_col <= size(work)) then - ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then - aux => work(n_col+1:) - else - allocate(aux(4*n_col),stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/4*n_col,0,0,0,0/),& - & a_err='real(psb_spk_)') - goto 9999 - end if - endif - else - allocate(ww(n_col),aux(4*n_col),stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/5*n_col,0,0,0,0/),& - & a_err='real(psb_spk_)') - goto 9999 - end if - endif - - - select case(prec%iprcparm(mld_sub_solve_)) - case(mld_ilu_n_,mld_milu_n_,mld_ilu_t_) - ! - ! Apply a block-Jacobi preconditioner with ILU(k)/MILU(k)/ILU(k,t) - ! factorization of the blocks (distributed matrix) or approximately - ! solve a system through ILU(k)/MILU(k)/ILU(k,t) (replicated matrix). - ! - - select case(trans_) - case('N') - - call psb_spsm(sone,prec%av(mld_l_pr_),x,szero,ww,desc_data,info,& - & trans=trans_,unit='L',diag=prec%d,choice=psb_none_,work=aux) - if (info == psb_success_) call psb_spsm(alpha,prec%av(mld_u_pr_),ww,beta,y,desc_data,info,& - & trans=trans_,unit='U',choice=psb_none_, work=aux) - - case('T','C') - call psb_spsm(sone,prec%av(mld_u_pr_),x,szero,ww,desc_data,info,& - & trans=trans_,unit='L',diag=prec%d,choice=psb_none_,work=aux) - if (info == psb_success_) call psb_spsm(alpha,prec%av(mld_l_pr_),ww,beta,y,desc_data,info,& - & trans=trans_,unit='U',choice=psb_none_,work=aux) - case default - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in ILU subsolve') - goto 9999 - end select - - case(mld_slu_) - ! - ! Apply a block-Jacobi preconditioner with LU factorization of the - ! blocks (distributed matrix) or approximately solve a local linear - ! system through LU (replicated matrix). The SuperLU package is used - ! to apply the LU factorization in both cases. - ! - - ww(1:n_row) = x(1:n_row) - - select case(trans_) - case('N') - call mld_sslu_solve(0,n_row,1,ww,n_row,prec%iprcparm(mld_slu_ptr_),info) - case('T','C') - call mld_sslu_solve(1,n_row,1,ww,n_row,prec%iprcparm(mld_slu_ptr_),info) - case default - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in SLU subsolve') - goto 9999 - end select - - if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info) - - case(mld_sludist_) - ! - ! Solve a distributed linear system with the LU factorization. - ! The SuperLU_DIST package is used. - ! - - ww(1:n_row) = x(1:n_row) - - select case(trans_) - case('N') - call mld_ssludist_solve(0,n_row,1,ww,n_row,prec%iprcparm(mld_slud_ptr_),info) - case('T','C') - call mld_ssludist_solve(1,n_row,1,ww,n_row,prec%iprcparm(mld_slud_ptr_),info) - case default - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in SLUDist subsolve') - goto 9999 - end select - - if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info) - - case (mld_umf_) - ! - ! Apply a block-Jacobi preconditioner with LU factorization of the - ! blocks (distributed matrix) or approximately solve a local linear - ! system through LU (replicated matrix). The UMFPACK package is used - ! to apply the LU factorization in both cases. - ! - - select case(trans_) - case('N') - call mld_sumf_solve(0,n_row,ww,x,n_row,prec%iprcparm(mld_umf_numptr_),info) - case('T','C') - call mld_sumf_solve(1,n_row,ww,x,n_row,prec%iprcparm(mld_umf_numptr_),info) - case default - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in UMF subsolve') - goto 9999 - end select - - if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info) - - case default - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid mld_sub_solve_') - goto 9999 - - end select - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error in subsolve') - goto 9999 - endif - - if (n_col <= size(work)) then - if ((4*n_col+n_col) <= size(work)) then - else - deallocate(aux) - endif - else - deallocate(ww,aux) - endif - - call psb_erractionrestore(err_act) - return - -9999 continue - - call psb_erractionrestore(err_act) - if (err_act.eq.psb_act_abort_) then - call psb_error() - return - end if - return - -end subroutine mld_ssub_solve - diff --git a/mlprec/mld_zsub_aply.f90 b/mlprec/mld_zsub_aply.f90 deleted file mode 100644 index 2629dc5e..00000000 --- a/mlprec/mld_zsub_aply.f90 +++ /dev/null @@ -1,297 +0,0 @@ -!!$ -!!$ -!!$ MLD2P4 version 2.0 -!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package -!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) -!!$ -!!$ (C) Copyright 2008,2009,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. -!!$ -!!$ -! File mld_zsub_aply.f90 -! -! Subroutine: mld_zsub_aply -! Version: complex -! -! This routine computes -! -! Y = beta*Y + alpha*op(K^(-1))*X, -! -! where -! - K is a suitable matrix, as specified below, -! - op(K^(-1)) is K^(-1) or its transpose, according to the value of the -! argument trans, -! - X and Y are vectors, -! - alpha and beta are scalars. -! -! Depending on K, alpha and beta (and on the communication descriptor desc_data -! - see the arguments below), the above computation may correspond to one of -! the following tasks: -! -! 1. Application of a block-Jacobi preconditioner associated to a matrix A -! distributed among the processes. Here K is the preconditioner, op(K^(-1)) -! = K^(-1), alpha = 1 and beta = 0. -! -! 2. Application of block-Jacobi sweeps to compute an approximate solution of -! a linear system -! A*Y = X, -! -! distributed among the processes (note that a single block-Jacobi sweep, -! with null starting guess, corresponds to the application of a block-Jacobi -! preconditioner). Here K^(-1) denotes the iteration matrix of the -! block-Jacobi solver, op(K^(-1)) = K^(-1), alpha = 1 and beta = 0. -! -! 3. Solution, through the LU factorization, of a linear system -! -! A*Y = X, -! -! distributed among the processes. Here K = L*U = A, op(K^(-1)) = K^(-1), -! alpha = 1 and beta = 0. -! -! 4. (Approximate) solution, through the LU or incomplete LU factorization, of -! a linear system -! A*Y = X, -! -! replicated on the processes. Here K = L*U = A or K = L*U ~ A, op(K^(-1)) = -! K^(-1), alpha = 1 and beta = 0. -! -! The block-Jacobi preconditioner or solver and the L and U factors of the LU -! or ILU factorizations have been built by the routine mld_fact_bld and stored -! into the 'base preconditioner' data structure prec. See mld_fact_bld for more -! details. -! -! This routine is used by mld_as_aply, to apply a 'base' block-Jacobi or -! Additive Schwarz (AS) preconditioner at any level of a multilevel preconditioner, -! or a block-Jacobi or LU or ILU solver at the coarsest level of a multilevel -! preconditioner. -! -! Tasks 1, 3 and 4 may be selected when prec%iprcparm(mld_smoother_sweeps_) = 1, -! while task 2 is selected when prec%iprcparm(mld_smoother_sweeps_) > 1. -! Furthermore, tasks 1, 2 and 3 may be performed when the matrix A is distributed -! among the processes (p%precv(ilev)%iprcparm(mld_coarse_mat_) = mld_distr_mat_, -! where p%precv(ilev) is the one-level data structure associated to the level -! ilev at which mld_sub_aply is called), while task 4 may be performed when A -! is replicated on the processes (p%precv(ilev)%iprcparm(mld_coarse_mat_) = -! mld_repl_mat_). Note that the matrix A is distributed among the processes -! at each level of the multilevel preconditioner, except the coarsest one, where -! it may be either distributed or replicated on the processes. Tasks 2, 3 and 4 -! are performed only at the coarsest level. Note also that this routine manages -! implicitly the fact that the matrix is distributed or replicated, i.e. it does not -! make any explicit reference to the value of p%precv(ilev)%iprcparm(mld_coarse_mat_). -! -! Arguments: -! -! alpha - complex(psb_dpk_), input. -! The scalar alpha. -! prec - type(mld_zbaseprec_type), input. -! The 'base preconditioner' data structure containing the local -! part of the preconditioner or solver. -! x - complex(psb_dpk_), dimension(:), input. -! The local part of the vector X. -! beta - complex(psb_dpk_), input. -! The scalar beta. -! y - complex(psb_dpk_), dimension(:), input/output. -! The local part of the vector Y. -! desc_data - type(psb_desc_type), input. -! The communication descriptor associated to the matrix to be -! preconditioned or 'inverted'. -! trans - character(len=1), input. -! If trans='N','n' then op(K^(-1)) = K^(-1); -! if trans='T','t' then op(K^(-1)) = K^(-T) (transpose of K^(-1)). -! if trans='C','c' then op(K^(-1)) = K^(-C) (transpose conjugate of K^(-1)). -! If prec%iprcparm(mld_smoother_sweeps_) > 1, the value of trans provided -! in input is ignored. -! work - complex(psb_dpk_), dimension (:), target. -! Workspace. Its size must be at least 4*psb_cd_get_local_cols(desc_data). -! info - integer, output. -! Error code. -! -subroutine mld_zsub_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) - - use psb_sparse_mod - use mld_inner_mod, mld_protect_name => mld_zsub_aply - - implicit none - - ! Arguments - type(psb_desc_type), intent(in) :: desc_data - type(mld_zbaseprec_type), intent(in) :: prec - complex(psb_dpk_),intent(in) :: x(:) - complex(psb_dpk_),intent(inout) :: y(:) - complex(psb_dpk_),intent(in) :: alpha,beta - character(len=1), intent(in) :: trans - complex(psb_dpk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - ! Local variables - integer :: n_row,n_col - complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character(len=20) :: name - character :: trans_ - - name='mld_zsub_aply' - info = psb_success_ - call psb_erractionsave(err_act) - - ictxt=psb_cd_get_context(desc_data) - call psb_info(ictxt, me, np) - - 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 - - - n_row = psb_cd_get_local_rows(desc_data) - n_col = psb_cd_get_local_cols(desc_data) - - if (n_col <= size(work)) then - ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then - aux => work(n_col+1:) - else - allocate(aux(4*n_col),stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/4*n_col,0,0,0,0/),& - & a_err='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 (prec%iprcparm(mld_smoother_sweeps_) == 1) then - - call mld_sub_solve(alpha,prec,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 (prec%iprcparm(mld_smoother_sweeps_) > 1) then - ! - ! - ! Apply prec%iprcparm(mld_smoother_sweeps_) sweeps of a block-Jacobi solver - ! to compute an approximate solution of a linear system. - ! - - if (size(prec%av) < mld_ap_nd_) then - info = psb_err_from_subroutine_non_ - goto 9999 - endif - - 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, prec%iprcparm(mld_smoother_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,prec%av(mld_ap_nd_),tx,zone,ty,& - & prec%desc_data,info,work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call mld_sub_solve(zone,prec,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,prec%iprcparm(mld_smoother_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.eq.psb_act_abort_) then - call psb_error() - return - end if - return - -end subroutine mld_zsub_aply - diff --git a/mlprec/mld_zsub_solve.f90 b/mlprec/mld_zsub_solve.f90 deleted file mode 100644 index 8f7d7372..00000000 --- a/mlprec/mld_zsub_solve.f90 +++ /dev/null @@ -1,325 +0,0 @@ -!!$ -!!$ -!!$ MLD2P4 version 2.0 -!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package -!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) -!!$ -!!$ (C) Copyright 2008,2009,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. -!!$ -!!$ -! File mld_zsub_solve.f90 -! -! Subroutine: mld_zsub_solve -! Version: complex -! -! This routine computes -! -! Y = beta*Y + alpha*op(K^(-1))*X, -! -! where -! - K is a factored matrix, as specified below, -! - op(K^(-1)) is K^(-1) or its transpose, according to the value of the -! argument trans, -! - X and Y are vectors, -! - alpha and beta are scalars. -! -! Depending on K, alpha and beta (and on the communication descriptor desc_data -! - see the arguments below), the above computation may correspond to one of -! the following tasks: -! -! 1. approximate solution of a linear system -! -! A*Y = X, -! -! by using the L and U factors computed with an ILU factorization of A. -! In this case K = L*U ~ A, alpha = 1 and beta = 0. The factors L and U -! (and the matrix A) are either distributed and block-diagonal or replicated. -! -! 2. Solution of a linear system -! -! A*Y = X, -! -! by using the L and U factors computed with a LU factorization of A. In this -! case K = L*U = A, alpha = 1 and beta = 0. The LU factorization is performed -! by one of the following auxiliary pakages: -! a. UMFPACK, -! b. SuperLU, -! c. SuperLU_Dist. -! In the cases a. and b., the factors L and U (and the matrix A) are either -! distributed and block diagonal) or replicated; in the case c., L, U (and A) -! are distributed. -! -! This routine is used by mld_dsub_aply, to apply a 'base' block-Jacobi or -! Additive Schwarz (AS) preconditioner at any level of a multilevel preconditioner, -! or a block-Jacobi or LU or ILU solver at the coarsest level of a multilevel -! preconditioner. -! -! -! Arguments: -! -! alpha - complex(psb_dpk_), input. -! The scalar alpha. -! prec - type(mld_zbaseprec_type), input. -! The 'base preconditioner' data structure containing the local -! part of the L and U factors of the matrix A. -! x - complex(psb_dpk_), dimension(:), input. -! The local part of the vector X. -! beta - complex(psb_dpk_), input. -! The scalar beta. -! y - complex(psb_dpk_), dimension(:), input/output. -! The local part of the vector Y. -! desc_data - type(psb_desc_type), input. -! The communication descriptor associated to the matrix to be -! preconditioned or 'inverted'. -! trans - character(len=1), input. -! If trans='N','n' then op(K^(-1)) = K^(-1); -! if trans='T','t' then op(K^(-1)) = K^(-T) (transpose of K^(-1)). -! if trans='C','c' then op(K^(-1)) = K^(-C) (transpose conjugate of K^(-1)). -! If prec%iprcparm(mld_smoother_sweeps_) > 1, the value of trans provided -! in input is ignored. -! work - complex(psb_dpk_), dimension (:), target. -! Workspace. Its size must be at least 4*psb_cd_get_local_cols(desc_data). -! info - integer, output. -! Error code. -! -subroutine mld_zsub_solve(alpha,prec,x,beta,y,desc_data,trans,work,info) - - use psb_sparse_mod - use mld_inner_mod, mld_protect_name => mld_zsub_solve - - implicit none - - ! Arguments - type(psb_desc_type), intent(in) :: desc_data - type(mld_zbaseprec_type), intent(in) :: prec - complex(psb_dpk_),intent(in) :: x(:) - complex(psb_dpk_),intent(inout) :: y(:) - complex(psb_dpk_),intent(in) :: alpha,beta - character(len=1), intent(in) :: trans - complex(psb_dpk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - ! Local variables - integer :: n_row,n_col - complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character(len=20) :: name - character :: trans_ - - interface - subroutine mld_zumf_solve(flag,m,x,b,n,ptr,info) - use psb_sparse_mod - integer, intent(in) :: flag,m,n,ptr - integer, intent(out) :: info - complex(psb_dpk_), intent(in) :: b(*) - complex(psb_dpk_), intent(inout) :: x(*) - end subroutine mld_zumf_solve - end interface - - name='mld_zsub_solve' - info = psb_success_ - call psb_erractionsave(err_act) - - ictxt=psb_cd_get_context(desc_data) - call psb_info(ictxt, me, np) - - 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 - - - n_row = psb_cd_get_local_rows(desc_data) - n_col = psb_cd_get_local_cols(desc_data) - - if (n_col <= size(work)) then - ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then - aux => work(n_col+1:) - else - allocate(aux(4*n_col),stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/4*n_col,0,0,0,0/),& - & a_err='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 - - - select case(prec%iprcparm(mld_sub_solve_)) - case(mld_ilu_n_,mld_milu_n_,mld_ilu_t_) - ! - ! Apply a block-Jacobi preconditioner with ILU(k)/MILU(k)/ILU(k,t) - ! factorization of the blocks (distributed matrix) or approximately - ! solve a system through ILU(k)/MILU(k)/ILU(k,t) (replicated matrix). - ! - - select case(trans_) - case('N') - - call psb_spsm(zone,prec%av(mld_l_pr_),x,zzero,ww,desc_data,info,& - & trans=trans_,unit='L',diag=prec%d,choice=psb_none_,work=aux) - if (info == psb_success_) call psb_spsm(alpha,prec%av(mld_u_pr_),ww,beta,y,desc_data,info,& - & trans=trans_,unit='U',choice=psb_none_, work=aux) - - case('T') - call psb_spsm(zone,prec%av(mld_u_pr_),x,zzero,ww,desc_data,info,& - & trans=trans_,unit='L',diag=prec%d,choice=psb_none_, work=aux) - if(info == psb_success_) call psb_spsm(alpha,prec%av(mld_l_pr_),ww,beta,y,desc_data,info,& - & trans=trans_,unit='U',choice=psb_none_,work=aux) - - case('C') - call psb_spsm(zone,prec%av(mld_u_pr_),x,zzero,ww,desc_data,info,& - & trans=trans_,unit='L',diag=conjg(prec%d),choice=psb_none_, work=aux) - if(info == psb_success_) call psb_spsm(alpha,prec%av(mld_l_pr_),ww,beta,y,desc_data,info,& - & trans=trans_,unit='U',choice=psb_none_,work=aux) - case default - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in ILU subsolve') - goto 9999 - end select - - case(mld_slu_) - ! - ! Apply a block-Jacobi preconditioner with LU factorization of the - ! blocks (distributed matrix) or approximately solve a local linear - ! system through LU (replicated matrix). The SuperLU package is used - ! to apply the LU factorization in both cases. - ! - - ww(1:n_row) = x(1:n_row) - - select case(trans_) - case('N') - call mld_zslu_solve(0,n_row,1,ww,n_row,prec%iprcparm(mld_slu_ptr_),info) - case('T') - call mld_zslu_solve(1,n_row,1,ww,n_row,prec%iprcparm(mld_slu_ptr_),info) - case('C') - call mld_zslu_solve(2,n_row,1,ww,n_row,prec%iprcparm(mld_slu_ptr_),info) - case default - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in SLU subsolve') - goto 9999 - end select - - if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info) - - case(mld_sludist_) - ! - ! Solve a distributed linear system with the LU factorization. - ! The SuperLU_DIST package is used. - ! - - ww(1:n_row) = x(1:n_row) - - select case(trans_) - case('N') - call mld_zsludist_solve(0,n_row,1,ww,n_row,prec%iprcparm(mld_slud_ptr_),info) - case('T') - call mld_zsludist_solve(1,n_row,1,ww,n_row,prec%iprcparm(mld_slud_ptr_),info) - case('C') - call mld_zsludist_solve(2,n_row,1,ww,n_row,prec%iprcparm(mld_slud_ptr_),info) - case default - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in SLUDist subsolve') - goto 9999 - end select - - if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info) - - case (mld_umf_) - ! - ! Apply a block-Jacobi preconditioner with LU factorization of the - ! blocks (distributed matrix) or approximately solve a local linear - ! system through LU (replicated matrix). The UMFPACK package is used - ! to apply the LU factorization in both cases. - ! - - select case(trans_) - case('N') - call mld_zumf_solve(0,n_row,ww,x,n_row,prec%iprcparm(mld_umf_numptr_),info) - case('T') - call mld_zumf_solve(1,n_row,ww,x,n_row,prec%iprcparm(mld_umf_numptr_),info) - case('C') - call mld_zumf_solve(2,n_row,ww,x,n_row,prec%iprcparm(mld_umf_numptr_),info) - case default - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in UMF subsolve') - goto 9999 - end select - - if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info) - - case default - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid mld_sub_solve_') - goto 9999 - - end select - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Error in subsolve ') - goto 9999 - endif - - if (n_col <= size(work)) then - if ((4*n_col+n_col) <= size(work)) then - else - deallocate(aux) - endif - else - deallocate(ww,aux) - endif - - call psb_erractionrestore(err_act) - return - -9999 continue - - call psb_erractionrestore(err_act) - if (err_act.eq.psb_act_abort_) then - call psb_error() - return - end if - return - -end subroutine mld_zsub_solve -