From 3e22e9e963f572d70d991726d489de1fae87b917 Mon Sep 17 00:00:00 2001 From: Cirdans-Home Date: Thu, 19 Nov 2020 16:02:15 +0100 Subject: [PATCH] implemented biconjugation --- prec/Makefile | 12 +- prec/impl/Makefile | 2 + prec/impl/psb_c_ainv_bld.f90 | 224 ++++++++++++++++++ prec/impl/psb_c_bjacprec_impl.f90 | 189 ++++++++++++++- prec/impl/psb_c_prec_type_impl.f90 | 19 ++ prec/impl/psb_c_sparsify.f90 | 8 +- prec/impl/psb_crwclip.f90 | 4 +- prec/impl/psb_csparse_biconjg_llk.F90 | 3 +- prec/impl/psb_csparse_biconjg_llk_noth.F90 | 2 +- prec/impl/psb_csparse_biconjg_mlk.F90 | 2 +- prec/impl/psb_csparse_biconjg_s_ft_llk.F90 | 8 +- prec/impl/psb_csparse_biconjg_s_llk.F90 | 2 +- prec/impl/psb_d_ainv_bld.f90 | 224 ++++++++++++++++++ prec/impl/psb_d_bjacprec_impl.f90 | 189 ++++++++++++++- prec/impl/psb_d_prec_type_impl.f90 | 19 ++ prec/impl/psb_d_sparsify.f90 | 8 +- prec/impl/psb_drwclip.f90 | 4 +- prec/impl/psb_dsparse_biconjg_llk.F90 | 3 +- prec/impl/psb_dsparse_biconjg_llk_noth.F90 | 2 +- prec/impl/psb_dsparse_biconjg_mlk.F90 | 2 +- prec/impl/psb_dsparse_biconjg_s_ft_llk.F90 | 2 +- prec/impl/psb_dsparse_biconjg_s_llk.F90 | 2 +- prec/impl/psb_s_ainv_bld.f90 | 224 ++++++++++++++++++ prec/impl/psb_s_bjacprec_impl.f90 | 189 ++++++++++++++- prec/impl/psb_s_prec_type_impl.f90 | 19 ++ prec/impl/psb_s_sparsify.f90 | 8 +- prec/impl/psb_srwclip.f90 | 4 +- prec/impl/psb_ssparse_biconjg_llk.F90 | 3 +- prec/impl/psb_ssparse_biconjg_llk_noth.F90 | 2 +- prec/impl/psb_ssparse_biconjg_mlk.F90 | 2 +- prec/impl/psb_ssparse_biconjg_s_ft_llk.F90 | 8 +- prec/impl/psb_ssparse_biconjg_s_llk.F90 | 2 +- prec/impl/psb_z_ainv_bld.f90 | 224 ++++++++++++++++++ prec/impl/psb_z_bjacprec_impl.f90 | 189 ++++++++++++++- prec/impl/psb_z_prec_type_impl.f90 | 19 ++ prec/impl/psb_z_sparsify.f90 | 8 +- prec/impl/psb_zrwclip.f90 | 4 +- prec/impl/psb_zsparse_biconjg_llk.F90 | 3 +- prec/impl/psb_zsparse_biconjg_llk_noth.F90 | 2 +- prec/impl/psb_zsparse_biconjg_mlk.F90 | 2 +- prec/impl/psb_zsparse_biconjg_s_ft_llk.F90 | 8 +- prec/impl/psb_zsparse_biconjg_s_llk.F90 | 2 +- ...e_ainv_mod.F90 => psb_c_ainv_fact_mod.f90} | 70 +++--- prec/psb_c_biconjg_mod.F90 | 7 +- prec/psb_c_bjacprec.f90 | 5 +- prec/psb_d_ainv_fact_mod.f90 | 98 ++++++++ prec/psb_d_biconjg_mod.F90 | 7 +- prec/psb_d_bjacprec.f90 | 5 +- prec/psb_prec_const_mod.f90 | 18 +- prec/psb_s_ainv_fact_mod.f90 | 98 ++++++++ prec/psb_s_biconjg_mod.F90 | 7 +- prec/psb_s_bjacprec.f90 | 5 +- prec/psb_z_ainv_fact_mod.f90 | 98 ++++++++ prec/psb_z_biconjg_mod.F90 | 7 +- prec/psb_z_bjacprec.f90 | 5 +- test/pargen/psb_d_pde3d.f90 | 16 +- test/pargen/psb_s_pde3d.f90 | 16 +- 57 files changed, 2140 insertions(+), 174 deletions(-) create mode 100644 prec/impl/psb_c_ainv_bld.f90 create mode 100644 prec/impl/psb_d_ainv_bld.f90 create mode 100644 prec/impl/psb_s_ainv_bld.f90 create mode 100644 prec/impl/psb_z_ainv_bld.f90 rename prec/{psb_base_ainv_mod.F90 => psb_c_ainv_fact_mod.f90} (64%) create mode 100644 prec/psb_d_ainv_fact_mod.f90 create mode 100644 prec/psb_s_ainv_fact_mod.f90 create mode 100644 prec/psb_z_ainv_fact_mod.f90 diff --git a/prec/Makefile b/prec/Makefile index b44172ad..5b887780 100644 --- a/prec/Makefile +++ b/prec/Makefile @@ -52,12 +52,14 @@ psb_s_bjacprec.o psb_s_diagprec.o psb_s_nullprec.o: psb_prec_mod.o psb_s_base_pr psb_d_bjacprec.o psb_d_diagprec.o psb_d_nullprec.o: psb_prec_mod.o psb_d_base_prec_mod.o psb_c_bjacprec.o psb_c_diagprec.o psb_c_nullprec.o: psb_prec_mod.o psb_c_base_prec_mod.o psb_z_bjacprec.o psb_z_diagprec.o psb_z_nullprec.o: psb_prec_mod.o psb_z_base_prec_mod.o -psb_s_bjacprec.o: psb_s_ilu_fact_mod.o -psb_d_bjacprec.o: psb_d_ilu_fact_mod.o -psb_c_bjacprec.o: psb_c_ilu_fact_mod.o -psb_z_bjacprec.o: psb_z_ilu_fact_mod.o +psb_s_bjacprec.o: psb_s_ilu_fact_mod.o psb_s_ainv_fact_mod.o +psb_d_bjacprec.o: psb_d_ilu_fact_mod.o psb_d_ainv_fact_mod.o +psb_c_bjacprec.o: psb_c_ilu_fact_mod.o psb_c_ainv_fact_mod.o +psb_z_bjacprec.o: psb_z_ilu_fact_mod.o psb_z_ainv_fact_mod.o +psb_c_ainv_fact_mod.o: psb_prec_const_mod.o psb_ainv_tools_mod.o psb_s_ainv_fact_mod.o psb_d_ainv_fact_mod.o psb_c_ainv_fact_mod.o psb_z_ainv_fact_mod.o psb_ainv_tools_mod.o: psb_c_ainv_tools_mod.o psb_d_ainv_tools_mod.o psb_s_ainv_tools_mod.o psb_z_ainv_tools_mod.o -psb_biconjg_mod.o: psb_base_ainv_mod.o psb_c_biconjg_mod.o psb_d_biconjg_mod.o psb_s_biconjg_mod.o psb_z_biconjg_mod.o +psb_biconjg_mod.o: psb_prec_const_mod.o psb_c_biconjg_mod.o \ + psb_d_biconjg_mod.o psb_s_biconjg_mod.o psb_z_biconjg_mod.o veryclean: clean /bin/rm -f $(LIBNAME) *$(.mod) diff --git a/prec/impl/Makefile b/prec/impl/Makefile index 111c8d30..284eb6e0 100644 --- a/prec/impl/Makefile +++ b/prec/impl/Makefile @@ -33,6 +33,8 @@ OBJS=psb_s_prec_type_impl.o psb_d_prec_type_impl.o \ psb_ssparse_biconjg_llk_noth.o psb_ssparse_biconjg_llk.o \ psb_ssparse_biconjg_mlk.o psb_ssparse_biconjg_s_ft_llk.o \ psb_ssparse_biconjg_s_llk.o \ + psb_d_ainv_bld.o psb_c_ainv_bld.o psb_s_ainv_bld.o \ + psb_z_ainv_bld.o LIBNAME=$(PRECLIBNAME) COBJS= diff --git a/prec/impl/psb_c_ainv_bld.f90 b/prec/impl/psb_c_ainv_bld.f90 new file mode 100644 index 00000000..48e6993d --- /dev/null +++ b/prec/impl/psb_c_ainv_bld.f90 @@ -0,0 +1,224 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! 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 PSBLAS 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 PSBLAS 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. +! +! Moved here from AMG-AINV, original copyright below. +! +! +! AMG-AINV: Approximate Inverse plugin for +! AMG4PSBLAS version 1.0 +! +! (C) Copyright 2020 +! +! Salvatore Filippone University of Rome Tor Vergata +! +! 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 AMG4PSBLAS 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 AMG4PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +subroutine psb_c_ainv_bld(a,alg,fillin,thresh,wmat,d,zmat,desc,info,blck,iscale) + + use psb_base_mod + use psb_prec_const_mod + use psb_c_biconjg_mod + + implicit none + + ! Arguments + type(psb_cspmat_type), intent(in), target :: a + integer(psb_ipk_), intent(in) :: fillin,alg + real(psb_spk_), intent(in) :: thresh + type(psb_cspmat_type), intent(inout) :: wmat, zmat + complex(psb_spk_), allocatable :: d(:) + Type(psb_desc_type), Intent(in) :: desc + integer(psb_ipk_), intent(out) :: info + type(psb_cspmat_type), intent(in), optional :: blck + integer(psb_ipk_), intent(in), optional :: iscale + ! + integer(psb_ipk_) :: i, nztota, err_act, n_row, nrow_a + type(psb_c_coo_sparse_mat) :: acoo + type(psb_c_csr_sparse_mat) :: acsr + type(psb_cspmat_type) :: atmp + real(psb_spk_), allocatable :: arws(:), acls(:) + complex(psb_spk_), allocatable :: pq(:), ad(:) + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: nzrmax, iscale_ + real(psb_spk_) :: sp_thresh + complex(psb_spk_) :: weight + character(len=20) :: name, ch_err + + + info = psb_success_ + name = 'psb_cainv_bld' + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = psb_cd_get_context(desc) + call psb_info(ictxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + iscale_ = psb_ilu_scale_none_ + if (present(iscale)) iscale_ = iscale + weight = cone + ! + ! Check the memory available to hold the W and Z factors + ! and allocate it if needed + ! + nrow_a = a%get_nrows() + nztota = a%get_nzeros() + + if (present(blck)) then + nztota = nztota + blck%get_nzeros() + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & ': out get_nnzeros',nrow_a,nztota,& + & a%get_nrows(),a%get_ncols(),a%get_nzeros() + + + n_row = desc%get_local_rows() + allocate(pq(n_row),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + nzrmax = fillin + sp_thresh = thresh + + ! + ! Ok, let's start first with Z (i.e. Upper) + ! + call a%csclip(acoo,info,imax=n_row,jmax=n_row) + call acsr%mv_from_coo(acoo,info) + select case(iscale_) + case(psb_ilu_scale_none_) + ! Ok, do nothing. + + case(psb_ilu_scale_maxval_) + weight = acsr%maxval() + weight = cone/weight + call acsr%scal(weight,info) + + case(psb_ilu_scale_arcsum_) + allocate(arws(n_row),acls(n_row),ad(n_row),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + call acsr%arwsum(arws) + call acsr%aclsum(acls) + ad(1:n_row) = sqrt(sqrt(arws(1:n_row)*acls(1:n_row))) + ad(1:n_row) = cone/ad(1:n_row) + call acsr%scal(ad,info,side='L') + call acsr%scal(ad,info,side='R') + case default + call psb_errpush(psb_err_from_subroutine_,name,a_err='wrong iscale') + goto 9999 + end select + + ! + ! Here for the actual workhorses. + ! Only biconjg is surviving for now.... + ! + call psb_sparse_biconjg(alg,n_row,acsr,pq,& + & zmat,wmat,nzrmax,sp_thresh,info) + ! Done. Hopefully.... + + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='sparse_biconjg') + goto 9999 + end if + call atmp%mv_from(acsr) + + ! + ! Is this right??? + ! + do i=1, n_row + if (abs(pq(i)) < d_epstol) then + pq(i) = cone + else + pq(i) = cone/pq(i) + end if + end do + + select case(iscale_) + case(psb_ilu_scale_none_) + ! Ok, do nothing. + case(psb_ilu_scale_maxval_) + pq(:) = pq(:)*weight + + case(psb_ilu_scale_arcsum_) + call zmat%scal(ad,info,side='L') + call wmat%scal(ad,info,side='R') + + case default + call psb_errpush(psb_err_from_subroutine_,name,a_err='wrong iscale') + goto 9999 + end select + + call psb_move_alloc(pq,d,info) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_c_ainv_bld diff --git a/prec/impl/psb_c_bjacprec_impl.f90 b/prec/impl/psb_c_bjacprec_impl.f90 index 02e23da3..3bfbff64 100644 --- a/prec/impl/psb_c_bjacprec_impl.f90 +++ b/prec/impl/psb_c_bjacprec_impl.f90 @@ -195,6 +195,46 @@ subroutine psb_c_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) goto 9999 end if + case(psb_f_ainv_) + + select case(trans_) + case('N') + call psb_spmm(cone,prec%av(psb_l_pr_),x,czero,wv,desc_data,info,& + & trans=trans_, work=aux) + + call wv1%mlt(cone,prec%dv,wv,czero,info) + + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,& + & beta,y,desc_data,info,& + & trans=trans_, work=aux) + + case('T') + call psb_spmm(cone,prec%av(psb_u_pr_),x,czero,wv,desc_data,info,& + & trans=trans_, work=aux) + + call wv1%mlt(cone,prec%dv,wv,czero,info) + + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,& + & beta,y,desc_data,info,& + & trans=trans_,work=aux) + + case('C') + + call psb_spmm(cone,prec%av(psb_u_pr_),x,czero,wv,desc_data,info,& + & trans=trans_,work=aux) + + call wv1%mlt(cone,prec%dv,wv,czero,info,conjgx=trans_) + + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,& + & beta,y,desc_data,info,& + & trans=trans_,work=aux) + + end select + if (info /= psb_success_) then + ch_err="psb_spsm" + goto 9999 + end if + case default info = psb_err_internal_error_ call psb_errpush(info,name,a_err='Invalid factorization') @@ -242,6 +282,7 @@ subroutine psb_c_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) ! Local variables integer(psb_ipk_) :: n_row,n_col complex(psb_spk_), pointer :: ww(:), aux(:) + type(psb_d_vect_type) :: tx,ty integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act, ierr(5) integer(psb_ipk_) :: debug_level, debug_unit @@ -346,6 +387,29 @@ subroutine psb_c_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) goto 9999 end if + case(psb_f_ainv_) + ! Application of approximate inverse preconditioner, just some spmm + ! call psb_spmm(alpha, a, x, beta, y,desc_a, info, & + ! & trans, work) + + + select case(trans_) + case('N','T') + call psb_spmm(cone,prec%av(psb_l_pr_),x,czero,ww,desc_data,info,& + & trans=trans_, work=aux) + ww(1:n_row) = ww(1:n_row)*prec%dv%v%v(1:n_row) + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),ww,& + & beta,y,desc_data,info, trans=trans_, work=aux) + + case('C') + call psb_spmm(cone,prec%av(psb_l_pr_),x,czero,ww,desc_data,info,& + & trans=trans_, work=aux) + ww(1:n_row) = ww(1:n_row)*conjg(prec%dv%v%v(1:n_row)) + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),ww,& + & beta,y,desc_data,info, trans=trans_, work=aux) + + end select + case default info = psb_err_internal_error_ @@ -402,9 +466,13 @@ subroutine psb_c_bjac_precinit(prec,info) prec%iprcparm(psb_ilu_fill_in_) = 0 prec%iprcparm(psb_ilu_ialg_) = psb_ilu_n_ prec%iprcparm(psb_ilu_scale_) = psb_ilu_scale_none_ + prec%iprcparm(psb_inv_fillin_) = 0 + prec%iprcparm(psb_ainv_alg_) = psb_ainv_llk_ + prec%rprcparm(:) = 0 prec%rprcparm(psb_fact_eps_) = 1E-1_psb_spk_ + prec%rprcparm(psb_inv_thresh_) = 1E-1_psb_spk_ call psb_erractionrestore(err_act) @@ -420,6 +488,7 @@ subroutine psb_c_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) use psb_base_mod use psb_c_ilu_fact_mod + use psb_c_ainv_fact_mod use psb_c_bjacprec, psb_protect_name => psb_c_bjac_precbld Implicit None @@ -432,12 +501,12 @@ subroutine psb_c_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) class(psb_i_base_vect_type), intent(in), optional :: imold ! .. Local Scalars .. - integer(psb_ipk_) :: i, m, ialg, fill_in, iscale + integer(psb_ipk_) :: i, m, ialg, fill_in, iscale, inv_fill, iinvalg integer(psb_ipk_) :: ierr(5) character :: trans, unitd type(psb_cspmat_type), allocatable :: lf, uf complex(psb_spk_), allocatable :: dd(:) - real(psb_spk_) :: fact_eps + real(psb_spk_) :: fact_eps, inv_thresh integer(psb_ipk_) :: nztota, err_act, n_row, nrow_a,n_col, nhalo integer(psb_ipk_) :: ictxt,np,me character(len=20) :: name='c_bjac_precbld' @@ -468,8 +537,13 @@ subroutine psb_c_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) ! We check if all the information contained in the preconditioner structure ! are meaningful, otherwise we give an error and get out of the build ! procedure - ialg = prec%iprcparm(psb_ilu_ialg_) - if ((ialg == psb_ilu_n_).or.(ialg == psb_milu_n_).or.(ialg == psb_ilu_t_)) then + ialg = prec%iprcparm(psb_ilu_ialg_) ! Integer for ILU type algorithm + iinvalg = prec%iprcparm(psb_ainv_alg_) ! Integer for AINV type algorithm + if ((ialg == psb_ilu_n_).or.(ialg == psb_milu_n_).or.& + & (ialg == psb_ilu_t_).or.(iinvalg == psb_ainv_llk_).or.& + & (iinvalg == psb_ainv_s_llk_).or.(iinvalg == psb_ainv_s_ft_llk_).or.& + & (iinvalg == psb_ainv_llk_noth_).or.(iinvalg == psb_ainv_mlk_).or.& + & (iinvalg == psb_ainv_lmx_)) then ! Do nothing: admissible request else info=psb_err_from_subroutine_ @@ -492,24 +566,32 @@ subroutine psb_c_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) goto 9999 end if fact_eps = prec%rprcparm(psb_fact_eps_) - if(fact_eps > 1 ) then + if( (fact_eps > 1).and.(prec%iprcparm(psb_f_type_) == psb_f_ainv_) ) then info=psb_err_from_subroutine_ ch_err='psb_fact_eps_' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if - fill_in = prec%iprcparm(psb_ilu_fill_in_) - if(fill_in < 0) then + inv_thresh = prec%rprcparm(psb_inv_thresh_) + if( (inv_thresh > 1) ) then info=psb_err_from_subroutine_ - ch_err='psb_ilu_fill_in_' + ch_err='psb_fact_eps_' call psb_errpush(info,name,a_err=ch_err) goto 9999 - else if (fill_in == 0) then - ! If the requested level of fill is equal to zero, we default to the - ! specialized ILU(0) routine - prec%iprcparm(psb_f_type_) = psb_f_ilu_n_ end if - + fill_in = prec%iprcparm(psb_ilu_fill_in_) + if (prec%iprcparm(psb_f_type_) == psb_f_ilu_n_) then + if(fill_in < 0) then + info=psb_err_from_subroutine_ + ch_err='psb_ilu_fill_in_' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + else if (fill_in == 0) then + ! If the requested level of fill is equal to zero, we default to the + ! specialized ILU(0) routine + prec%iprcparm(psb_f_type_) = psb_f_ilu_n_ + end if + end if ! Select on the type of factorization to be used select case(prec%iprcparm(psb_f_type_)) @@ -732,6 +814,78 @@ subroutine psb_c_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) goto 9999 end if + case(psb_f_ainv_) + ! Approximate Inverse Factorizations based on variants of the incomplete + ! biconjugation algorithms + if (allocated(prec%av)) then + if (size(prec%av) < psb_bp_ilu_avsz) then + do i=1,size(prec%av) + call prec%av(i)%free() + enddo + deallocate(prec%av,stat=info) + endif + end if + if (.not.allocated(prec%av)) then + allocate(prec%av(psb_max_avsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + end if + endif + + nrow_a = desc_a%get_local_rows() + nztota = a%get_nzeros() + + n_col = desc_a%get_local_cols() + nhalo = n_col-nrow_a + n_row = nrow_a + + allocate(lf,uf,stat=info) + if (info == psb_success_) call lf%allocate(n_row,n_row,nztota) + if (info == psb_success_) call uf%allocate(n_row,n_row,nztota) + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(dd(n_row),stat=info) + if (info == psb_success_) then + allocate(prec%dv, stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(prec%dv%v,mold=vmold,stat=info) + else + allocate(psb_c_base_vect_type :: prec%dv%v,stat=info) + end if + end if + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + endif + ! Computin the factorization + call psb_ainv_fact(a,iinvalg,inv_fill,inv_thresh,lf,dd,uf,desc_a,info,iscale=iscale) + + if(info == psb_success_) then + call prec%av(psb_l_pr_)%mv_from(lf%a) + call prec%av(psb_u_pr_)%mv_from(uf%a) + call prec%av(psb_l_pr_)%set_asb() + call prec%av(psb_u_pr_)%set_asb() + call prec%av(psb_l_pr_)%trim() + call prec%av(psb_u_pr_)%trim() + call prec%dv%bld(dd) + ! call move_alloc(dd,prec%d) + else + info=psb_err_from_subroutine_ + ch_err='psb_ilut_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + case(psb_f_none_) info=psb_err_from_subroutine_ ch_err='Inconsistent prec psb_f_none_' @@ -792,6 +946,12 @@ subroutine psb_c_bjac_precseti(prec,what,val,info) case (psb_ilu_scale_) prec%iprcparm(psb_ilu_scale_) = val + case (psb_ainv_alg_) + prec%iprcparm(psb_ainv_alg_) = val + + case(psb_inv_fillin_) + prec%iprcparm(psb_inv_fillin_) = val + case default write(psb_err_unit,'(i0," is invalid, ignoring user specification")') what @@ -833,6 +993,9 @@ subroutine psb_c_bjac_precsetr(prec,what,val,info) case (psb_fact_eps_) prec%rprcparm(psb_fact_eps_) = val + case (psb_inv_thresh_) + prec%rprcparm(psb_inv_thresh_) = val + case default write(psb_err_unit,'(i0," is invalid, ignoring user specification")') what diff --git a/prec/impl/psb_c_prec_type_impl.f90 b/prec/impl/psb_c_prec_type_impl.f90 index 1d99a85d..6e941bb1 100644 --- a/prec/impl/psb_c_prec_type_impl.f90 +++ b/prec/impl/psb_c_prec_type_impl.f90 @@ -356,6 +356,8 @@ subroutine psb_ccprecseti(prec,what,val,info,ilev,ilmax,pos,idx) select case (psb_toupper(what)) case ("SUB_FILLIN") call prec%prec%precset(psb_ilu_fill_in_,val,info) + case('INV_FILLIN') + call prec%prec%precset(psb_inv_fillin_,val,info) case default info = psb_err_invalid_args_combination_ write(psb_err_unit,*) name,& @@ -390,6 +392,8 @@ subroutine psb_ccprecsetr(prec,what,val,info,ilev,ilmax,pos,idx) select case (psb_toupper(what)) case('SUB_ILUTHRS') call prec%prec%precset(psb_fact_eps_,val,info) + case('INV_THRESH') + call prec%prec%precset(psb_inv_thresh_,val,info) case default info = psb_err_invalid_args_combination_ write(psb_err_unit,*) name,& @@ -431,6 +435,8 @@ subroutine psb_ccprecsetc(prec,what,string,info,ilev,ilmax,pos,idx) case("ILUT") call prec%prec%precset(psb_f_type_,psb_f_ilu_t_,info) call prec%prec%precset(psb_ilu_ialg_,psb_ilu_t_,info) + case("AINV") + call prec%prec%precset(psb_f_type_,psb_f_ainv_,info) case default ! Default to ILU(0) factorization call prec%prec%precset(psb_f_type_,psb_f_ilu_n_,info) @@ -450,6 +456,19 @@ subroutine psb_ccprecsetc(prec,what,string,info,ilev,ilmax,pos,idx) case default call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_none_,info) end select + case ("AINV_ALG") + select case (psb_toupper(string)) + case("LLK") + call prec%prec%precset(psb_ainv_alg_,psb_ainv_llk_,info) + case("SYM-LLK") + call prec%prec%precset(psb_ainv_alg_,psb_ainv_s_llk_,info) + case("STAB-LLK") + call prec%prec%precset(psb_ainv_alg_,psb_ainv_s_ft_llk_,info) + case("MLK","LMX") + call prec%prec%precset(psb_ainv_alg_,psb_ainv_mlk_,info) + case default + call prec%prec%precset(psb_ainv_alg_,psb_ainv_llk_,info) + end select case default end select diff --git a/prec/impl/psb_c_sparsify.f90 b/prec/impl/psb_c_sparsify.f90 index b89ce3d0..a215117e 100644 --- a/prec/impl/psb_c_sparsify.f90 +++ b/prec/impl/psb_c_sparsify.f90 @@ -63,7 +63,7 @@ ! POSSIBILITY OF SUCH DAMAGE. ! ! -subroutine amg_c_sparsify(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,info,istart,iheap,ikr) +subroutine psb_c_sparsify(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,info,istart,iheap,ikr) use psb_base_mod implicit none @@ -177,10 +177,10 @@ subroutine amg_c_sparsify(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,info,istart,ihe return -end subroutine amg_c_sparsify +end subroutine psb_c_sparsify -subroutine amg_c_sparsify_list(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,lhead,listv,ikr,info) +subroutine psb_c_sparsify_list(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,lhead,listv,ikr,info) use psb_base_mod implicit none @@ -258,4 +258,4 @@ subroutine amg_c_sparsify_list(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,lhead,list return -end subroutine amg_c_sparsify_list +end subroutine psb_c_sparsify_list diff --git a/prec/impl/psb_crwclip.f90 b/prec/impl/psb_crwclip.f90 index 941725d2..ade1171f 100644 --- a/prec/impl/psb_crwclip.f90 +++ b/prec/impl/psb_crwclip.f90 @@ -63,7 +63,7 @@ ! POSSIBILITY OF SUCH DAMAGE. ! ! -subroutine amg_c_rwclip(nz,ia,ja,val,imin,imax,jmin,jmax) +subroutine psb_c_rwclip(nz,ia,ja,val,imin,imax,jmin,jmax) use psb_base_mod implicit none @@ -87,4 +87,4 @@ subroutine amg_c_rwclip(nz,ia,ja,val,imin,imax,jmin,jmax) end if end do nz = j -end subroutine amg_c_rwclip +end subroutine psb_c_rwclip diff --git a/prec/impl/psb_csparse_biconjg_llk.F90 b/prec/impl/psb_csparse_biconjg_llk.F90 index e7ba35a7..9f4c629f 100644 --- a/prec/impl/psb_csparse_biconjg_llk.F90 +++ b/prec/impl/psb_csparse_biconjg_llk.F90 @@ -34,7 +34,7 @@ ! subroutine psb_csparse_biconjg_llk(n,a,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod - use psb_ainv_tools_mod + use psb_c_ainv_tools_mod use psb_c_biconjg_mod, psb_protect_name => psb_csparse_biconjg_llk ! ! Left-looking variant @@ -224,6 +224,7 @@ subroutine psb_csparse_biconjg_llk(n,a,p,z,w,nzrmax,sp_thresh,info) ! ! Sparsify current ZVAL and put into ZMAT ! + write(psb_out_unit,'("z(1) = ",f16.14)') zval(1) call sparsify(i,nzrmax,sp_thresh,n,zval,nzrz,ia,val,info,iheap=heap,ikr=izkr) if (info /= psb_success_) then info = psb_err_internal_error_ diff --git a/prec/impl/psb_csparse_biconjg_llk_noth.F90 b/prec/impl/psb_csparse_biconjg_llk_noth.F90 index ed1c19ff..5fe472eb 100644 --- a/prec/impl/psb_csparse_biconjg_llk_noth.F90 +++ b/prec/impl/psb_csparse_biconjg_llk_noth.F90 @@ -34,7 +34,7 @@ ! subroutine psb_csparse_biconjg_llk_noth(n,a,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod - use psb_ainv_tools_mod + use psb_c_ainv_tools_mod use psb_c_biconjg_mod, psb_protect_name => psb_csparse_biconjg_llk_noth ! diff --git a/prec/impl/psb_csparse_biconjg_mlk.F90 b/prec/impl/psb_csparse_biconjg_mlk.F90 index 8e8b7a74..58838c82 100644 --- a/prec/impl/psb_csparse_biconjg_mlk.F90 +++ b/prec/impl/psb_csparse_biconjg_mlk.F90 @@ -34,7 +34,7 @@ ! subroutine psb_csparse_biconjg_mlk(n,a,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod - use psb_ainv_tools_mod + use psb_c_ainv_tools_mod use psb_c_biconjg_mod, psb_protect_name => psb_csparse_biconjg_mlk ! ! Left-looking variant diff --git a/prec/impl/psb_csparse_biconjg_s_ft_llk.F90 b/prec/impl/psb_csparse_biconjg_s_ft_llk.F90 index 910b74b9..64af2ed5 100644 --- a/prec/impl/psb_csparse_biconjg_s_ft_llk.F90 +++ b/prec/impl/psb_csparse_biconjg_s_ft_llk.F90 @@ -34,7 +34,7 @@ ! subroutine psb_csparse_biconjg_s_ft_llk(n,a,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod - use psb_ainv_tools_mod + use psb_c_ainv_tools_mod use psb_c_biconjg_mod, psb_protect_name => psb_csparse_biconjg_s_ft_llk ! @@ -164,7 +164,7 @@ subroutine psb_csparse_biconjg_s_ft_llk(n,a,p,z,w,nzrmax,sp_thresh,info) ip2 = w%icp(j+1) - 1 nzra = max(0,ip2 - ip1 + 1) nzww = 0 - call psb_d_spvspm(cone,a,nzra,w%ia(ip1:ip2),w%val(ip1:ip2),& + call psb_c_spvspm(cone,a,nzra,w%ia(ip1:ip2),w%val(ip1:ip2),& & czero,nzww,iww,ww,info) p(i) = psb_spge_dot(nzww,iww,ww,zval) @@ -299,7 +299,7 @@ subroutine psb_csparse_biconjg_s_ft_llk(n,a,p,z,w,nzrmax,sp_thresh,info) ip2 = z%icp(j+1) - 1 nzra = max(0,ip2 - ip1 + 1) nzww = 0 - call psb_d_spmspv(cone,ac,nzra,z%ia(ip1:ip2),z%val(ip1:ip2),& + call psb_c_spmspv(cone,ac,nzra,z%ia(ip1:ip2),z%val(ip1:ip2),& & czero,nzww,iww,ww,info) q(i) = psb_spge_dot(nzww,iww,ww,zval) @@ -384,7 +384,7 @@ subroutine psb_csparse_biconjg_s_ft_llk(n,a,p,z,w,nzrmax,sp_thresh,info) nzww = 0 nzrz = z%icp(i+1)-z%icp(i) ipz1 = z%icp(i) - call psb_d_spmspv(cone,ac,& + call psb_c_spmspv(cone,ac,& & nzrz,z%ia(ipz1:ipz1+nzrz-1),z%val(ipz1:ipz1+nzrz-1),& & czero,nzww,iww,ww,info) tmpq = psb_spdot_srtd(nzww,iww,ww,nzrw,ia,val) diff --git a/prec/impl/psb_csparse_biconjg_s_llk.F90 b/prec/impl/psb_csparse_biconjg_s_llk.F90 index 25d90f1d..b38b2a0c 100644 --- a/prec/impl/psb_csparse_biconjg_s_llk.F90 +++ b/prec/impl/psb_csparse_biconjg_s_llk.F90 @@ -34,7 +34,7 @@ ! subroutine psb_csparse_biconjg_s_llk(n,a,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod - use psb_ainv_tools_mod + use psb_c_ainv_tools_mod use psb_c_biconjg_mod, psb_protect_name => psb_csparse_biconjg_s_llk ! diff --git a/prec/impl/psb_d_ainv_bld.f90 b/prec/impl/psb_d_ainv_bld.f90 new file mode 100644 index 00000000..67099730 --- /dev/null +++ b/prec/impl/psb_d_ainv_bld.f90 @@ -0,0 +1,224 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! 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 PSBLAS 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 PSBLAS 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. +! +! Moved here from AMG-AINV, original copyright below. +! +! +! AMG-AINV: Approximate Inverse plugin for +! AMG4PSBLAS version 1.0 +! +! (C) Copyright 2020 +! +! Salvatore Filippone University of Rome Tor Vergata +! +! 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 AMG4PSBLAS 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 AMG4PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +subroutine psb_d_ainv_bld(a,alg,fillin,thresh,wmat,d,zmat,desc,info,blck,iscale) + + use psb_base_mod + use psb_prec_const_mod + use psb_d_biconjg_mod + + implicit none + + ! Arguments + type(psb_dspmat_type), intent(in), target :: a + integer(psb_ipk_), intent(in) :: fillin,alg + real(psb_dpk_), intent(in) :: thresh + type(psb_dspmat_type), intent(inout) :: wmat, zmat + real(psb_dpk_), allocatable :: d(:) + Type(psb_desc_type), Intent(in) :: desc + integer(psb_ipk_), intent(out) :: info + type(psb_dspmat_type), intent(in), optional :: blck + integer(psb_ipk_), intent(in), optional :: iscale + ! + integer(psb_ipk_) :: i, nztota, err_act, n_row, nrow_a + type(psb_d_coo_sparse_mat) :: acoo + type(psb_d_csr_sparse_mat) :: acsr + type(psb_dspmat_type) :: atmp + real(psb_dpk_), allocatable :: arws(:), acls(:) + real(psb_dpk_), allocatable :: pq(:), ad(:) + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: nzrmax, iscale_ + real(psb_dpk_) :: sp_thresh + real(psb_dpk_) :: weight + character(len=20) :: name, ch_err + + + info = psb_success_ + name = 'psb_dainv_bld' + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = psb_cd_get_context(desc) + call psb_info(ictxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + iscale_ = psb_ilu_scale_none_ + if (present(iscale)) iscale_ = iscale + weight = done + ! + ! Check the memory available to hold the W and Z factors + ! and allocate it if needed + ! + nrow_a = a%get_nrows() + nztota = a%get_nzeros() + + if (present(blck)) then + nztota = nztota + blck%get_nzeros() + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & ': out get_nnzeros',nrow_a,nztota,& + & a%get_nrows(),a%get_ncols(),a%get_nzeros() + + + n_row = desc%get_local_rows() + allocate(pq(n_row),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + nzrmax = fillin + sp_thresh = thresh + + ! + ! Ok, let's start first with Z (i.e. Upper) + ! + call a%csclip(acoo,info,imax=n_row,jmax=n_row) + call acsr%mv_from_coo(acoo,info) + select case(iscale_) + case(psb_ilu_scale_none_) + ! Ok, do nothing. + + case(psb_ilu_scale_maxval_) + weight = acsr%maxval() + weight = done/weight + call acsr%scal(weight,info) + + case(psb_ilu_scale_arcsum_) + allocate(arws(n_row),acls(n_row),ad(n_row),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + call acsr%arwsum(arws) + call acsr%aclsum(acls) + ad(1:n_row) = sqrt(sqrt(arws(1:n_row)*acls(1:n_row))) + ad(1:n_row) = done/ad(1:n_row) + call acsr%scal(ad,info,side='L') + call acsr%scal(ad,info,side='R') + case default + call psb_errpush(psb_err_from_subroutine_,name,a_err='wrong iscale') + goto 9999 + end select + + ! + ! Here for the actual workhorses. + ! Only biconjg is surviving for now.... + ! + call psb_sparse_biconjg(alg,n_row,acsr,pq,& + & zmat,wmat,nzrmax,sp_thresh,info) + ! Done. Hopefully.... + + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='sparse_biconjg') + goto 9999 + end if + call atmp%mv_from(acsr) + + ! + ! Is this right??? + ! + do i=1, n_row + if (abs(pq(i)) < d_epstol) then + pq(i) = done + else + pq(i) = done/pq(i) + end if + end do + + select case(iscale_) + case(psb_ilu_scale_none_) + ! Ok, do nothing. + case(psb_ilu_scale_maxval_) + pq(:) = pq(:)*weight + + case(psb_ilu_scale_arcsum_) + call zmat%scal(ad,info,side='L') + call wmat%scal(ad,info,side='R') + + case default + call psb_errpush(psb_err_from_subroutine_,name,a_err='wrong iscale') + goto 9999 + end select + + call psb_move_alloc(pq,d,info) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_d_ainv_bld diff --git a/prec/impl/psb_d_bjacprec_impl.f90 b/prec/impl/psb_d_bjacprec_impl.f90 index 898fd224..32f0f9e3 100644 --- a/prec/impl/psb_d_bjacprec_impl.f90 +++ b/prec/impl/psb_d_bjacprec_impl.f90 @@ -195,6 +195,46 @@ subroutine psb_d_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) goto 9999 end if + case(psb_f_ainv_) + + select case(trans_) + case('N') + call psb_spmm(done,prec%av(psb_l_pr_),x,dzero,wv,desc_data,info,& + & trans=trans_, work=aux) + + call wv1%mlt(done,prec%dv,wv,dzero,info) + + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,& + & beta,y,desc_data,info,& + & trans=trans_, work=aux) + + case('T') + call psb_spmm(done,prec%av(psb_u_pr_),x,dzero,wv,desc_data,info,& + & trans=trans_, work=aux) + + call wv1%mlt(done,prec%dv,wv,dzero,info) + + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,& + & beta,y,desc_data,info,& + & trans=trans_,work=aux) + + case('C') + + call psb_spmm(done,prec%av(psb_u_pr_),x,dzero,wv,desc_data,info,& + & trans=trans_,work=aux) + + call wv1%mlt(done,prec%dv,wv,dzero,info,conjgx=trans_) + + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,& + & beta,y,desc_data,info,& + & trans=trans_,work=aux) + + end select + if (info /= psb_success_) then + ch_err="psb_spsm" + goto 9999 + end if + case default info = psb_err_internal_error_ call psb_errpush(info,name,a_err='Invalid factorization') @@ -242,6 +282,7 @@ subroutine psb_d_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) ! Local variables integer(psb_ipk_) :: n_row,n_col real(psb_dpk_), pointer :: ww(:), aux(:) + type(psb_d_vect_type) :: tx,ty integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act, ierr(5) integer(psb_ipk_) :: debug_level, debug_unit @@ -346,6 +387,29 @@ subroutine psb_d_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) goto 9999 end if + case(psb_f_ainv_) + ! Application of approximate inverse preconditioner, just some spmm + ! call psb_spmm(alpha, a, x, beta, y,desc_a, info, & + ! & trans, work) + + + select case(trans_) + case('N','T') + call psb_spmm(done,prec%av(psb_l_pr_),x,dzero,ww,desc_data,info,& + & trans=trans_, work=aux) + ww(1:n_row) = ww(1:n_row)*prec%dv%v%v(1:n_row) + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),ww,& + & beta,y,desc_data,info, trans=trans_, work=aux) + + case('C') + call psb_spmm(done,prec%av(psb_l_pr_),x,dzero,ww,desc_data,info,& + & trans=trans_, work=aux) + ww(1:n_row) = ww(1:n_row)*(prec%dv%v%v(1:n_row)) + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),ww,& + & beta,y,desc_data,info, trans=trans_, work=aux) + + end select + case default info = psb_err_internal_error_ @@ -402,9 +466,13 @@ subroutine psb_d_bjac_precinit(prec,info) prec%iprcparm(psb_ilu_fill_in_) = 0 prec%iprcparm(psb_ilu_ialg_) = psb_ilu_n_ prec%iprcparm(psb_ilu_scale_) = psb_ilu_scale_none_ + prec%iprcparm(psb_inv_fillin_) = 0 + prec%iprcparm(psb_ainv_alg_) = psb_ainv_llk_ + prec%rprcparm(:) = 0 prec%rprcparm(psb_fact_eps_) = 1E-1_psb_dpk_ + prec%rprcparm(psb_inv_thresh_) = 1E-1_psb_dpk_ call psb_erractionrestore(err_act) @@ -420,6 +488,7 @@ subroutine psb_d_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) use psb_base_mod use psb_d_ilu_fact_mod + use psb_d_ainv_fact_mod use psb_d_bjacprec, psb_protect_name => psb_d_bjac_precbld Implicit None @@ -432,12 +501,12 @@ subroutine psb_d_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) class(psb_i_base_vect_type), intent(in), optional :: imold ! .. Local Scalars .. - integer(psb_ipk_) :: i, m, ialg, fill_in, iscale + integer(psb_ipk_) :: i, m, ialg, fill_in, iscale, inv_fill, iinvalg integer(psb_ipk_) :: ierr(5) character :: trans, unitd type(psb_dspmat_type), allocatable :: lf, uf real(psb_dpk_), allocatable :: dd(:) - real(psb_dpk_) :: fact_eps + real(psb_dpk_) :: fact_eps, inv_thresh integer(psb_ipk_) :: nztota, err_act, n_row, nrow_a,n_col, nhalo integer(psb_ipk_) :: ictxt,np,me character(len=20) :: name='d_bjac_precbld' @@ -468,8 +537,13 @@ subroutine psb_d_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) ! We check if all the information contained in the preconditioner structure ! are meaningful, otherwise we give an error and get out of the build ! procedure - ialg = prec%iprcparm(psb_ilu_ialg_) - if ((ialg == psb_ilu_n_).or.(ialg == psb_milu_n_).or.(ialg == psb_ilu_t_)) then + ialg = prec%iprcparm(psb_ilu_ialg_) ! Integer for ILU type algorithm + iinvalg = prec%iprcparm(psb_ainv_alg_) ! Integer for AINV type algorithm + if ((ialg == psb_ilu_n_).or.(ialg == psb_milu_n_).or.& + & (ialg == psb_ilu_t_).or.(iinvalg == psb_ainv_llk_).or.& + & (iinvalg == psb_ainv_s_llk_).or.(iinvalg == psb_ainv_s_ft_llk_).or.& + & (iinvalg == psb_ainv_llk_noth_).or.(iinvalg == psb_ainv_mlk_).or.& + & (iinvalg == psb_ainv_lmx_)) then ! Do nothing: admissible request else info=psb_err_from_subroutine_ @@ -492,24 +566,32 @@ subroutine psb_d_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) goto 9999 end if fact_eps = prec%rprcparm(psb_fact_eps_) - if(fact_eps > 1 ) then + if( (fact_eps > 1).and.(prec%iprcparm(psb_f_type_) == psb_f_ainv_) ) then info=psb_err_from_subroutine_ ch_err='psb_fact_eps_' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if - fill_in = prec%iprcparm(psb_ilu_fill_in_) - if(fill_in < 0) then + inv_thresh = prec%rprcparm(psb_inv_thresh_) + if( (inv_thresh > 1) ) then info=psb_err_from_subroutine_ - ch_err='psb_ilu_fill_in_' + ch_err='psb_fact_eps_' call psb_errpush(info,name,a_err=ch_err) goto 9999 - else if (fill_in == 0) then - ! If the requested level of fill is equal to zero, we default to the - ! specialized ILU(0) routine - prec%iprcparm(psb_f_type_) = psb_f_ilu_n_ end if - + fill_in = prec%iprcparm(psb_ilu_fill_in_) + if (prec%iprcparm(psb_f_type_) == psb_f_ilu_n_) then + if(fill_in < 0) then + info=psb_err_from_subroutine_ + ch_err='psb_ilu_fill_in_' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + else if (fill_in == 0) then + ! If the requested level of fill is equal to zero, we default to the + ! specialized ILU(0) routine + prec%iprcparm(psb_f_type_) = psb_f_ilu_n_ + end if + end if ! Select on the type of factorization to be used select case(prec%iprcparm(psb_f_type_)) @@ -732,6 +814,78 @@ subroutine psb_d_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) goto 9999 end if + case(psb_f_ainv_) + ! Approximate Inverse Factorizations based on variants of the incomplete + ! biconjugation algorithms + if (allocated(prec%av)) then + if (size(prec%av) < psb_bp_ilu_avsz) then + do i=1,size(prec%av) + call prec%av(i)%free() + enddo + deallocate(prec%av,stat=info) + endif + end if + if (.not.allocated(prec%av)) then + allocate(prec%av(psb_max_avsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + end if + endif + + nrow_a = desc_a%get_local_rows() + nztota = a%get_nzeros() + + n_col = desc_a%get_local_cols() + nhalo = n_col-nrow_a + n_row = nrow_a + + allocate(lf,uf,stat=info) + if (info == psb_success_) call lf%allocate(n_row,n_row,nztota) + if (info == psb_success_) call uf%allocate(n_row,n_row,nztota) + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(dd(n_row),stat=info) + if (info == psb_success_) then + allocate(prec%dv, stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(prec%dv%v,mold=vmold,stat=info) + else + allocate(psb_d_base_vect_type :: prec%dv%v,stat=info) + end if + end if + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + endif + ! Computin the factorization + call psb_ainv_fact(a,iinvalg,inv_fill,inv_thresh,lf,dd,uf,desc_a,info,iscale=iscale) + + if(info == psb_success_) then + call prec%av(psb_l_pr_)%mv_from(lf%a) + call prec%av(psb_u_pr_)%mv_from(uf%a) + call prec%av(psb_l_pr_)%set_asb() + call prec%av(psb_u_pr_)%set_asb() + call prec%av(psb_l_pr_)%trim() + call prec%av(psb_u_pr_)%trim() + call prec%dv%bld(dd) + ! call move_alloc(dd,prec%d) + else + info=psb_err_from_subroutine_ + ch_err='psb_ilut_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + case(psb_f_none_) info=psb_err_from_subroutine_ ch_err='Inconsistent prec psb_f_none_' @@ -792,6 +946,12 @@ subroutine psb_d_bjac_precseti(prec,what,val,info) case (psb_ilu_scale_) prec%iprcparm(psb_ilu_scale_) = val + case (psb_ainv_alg_) + prec%iprcparm(psb_ainv_alg_) = val + + case(psb_inv_fillin_) + prec%iprcparm(psb_inv_fillin_) = val + case default write(psb_err_unit,'(i0," is invalid, ignoring user specification")') what @@ -833,6 +993,9 @@ subroutine psb_d_bjac_precsetr(prec,what,val,info) case (psb_fact_eps_) prec%rprcparm(psb_fact_eps_) = val + case (psb_inv_thresh_) + prec%rprcparm(psb_inv_thresh_) = val + case default write(psb_err_unit,'(i0," is invalid, ignoring user specification")') what diff --git a/prec/impl/psb_d_prec_type_impl.f90 b/prec/impl/psb_d_prec_type_impl.f90 index 5c07bebd..c0376c6c 100644 --- a/prec/impl/psb_d_prec_type_impl.f90 +++ b/prec/impl/psb_d_prec_type_impl.f90 @@ -356,6 +356,8 @@ subroutine psb_dcprecseti(prec,what,val,info,ilev,ilmax,pos,idx) select case (psb_toupper(what)) case ("SUB_FILLIN") call prec%prec%precset(psb_ilu_fill_in_,val,info) + case('INV_FILLIN') + call prec%prec%precset(psb_inv_fillin_,val,info) case default info = psb_err_invalid_args_combination_ write(psb_err_unit,*) name,& @@ -390,6 +392,8 @@ subroutine psb_dcprecsetr(prec,what,val,info,ilev,ilmax,pos,idx) select case (psb_toupper(what)) case('SUB_ILUTHRS') call prec%prec%precset(psb_fact_eps_,val,info) + case('INV_THRESH') + call prec%prec%precset(psb_inv_thresh_,val,info) case default info = psb_err_invalid_args_combination_ write(psb_err_unit,*) name,& @@ -431,6 +435,8 @@ subroutine psb_dcprecsetc(prec,what,string,info,ilev,ilmax,pos,idx) case("ILUT") call prec%prec%precset(psb_f_type_,psb_f_ilu_t_,info) call prec%prec%precset(psb_ilu_ialg_,psb_ilu_t_,info) + case("AINV") + call prec%prec%precset(psb_f_type_,psb_f_ainv_,info) case default ! Default to ILU(0) factorization call prec%prec%precset(psb_f_type_,psb_f_ilu_n_,info) @@ -450,6 +456,19 @@ subroutine psb_dcprecsetc(prec,what,string,info,ilev,ilmax,pos,idx) case default call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_none_,info) end select + case ("AINV_ALG") + select case (psb_toupper(string)) + case("LLK") + call prec%prec%precset(psb_ainv_alg_,psb_ainv_llk_,info) + case("SYM-LLK") + call prec%prec%precset(psb_ainv_alg_,psb_ainv_s_llk_,info) + case("STAB-LLK") + call prec%prec%precset(psb_ainv_alg_,psb_ainv_s_ft_llk_,info) + case("MLK","LMX") + call prec%prec%precset(psb_ainv_alg_,psb_ainv_mlk_,info) + case default + call prec%prec%precset(psb_ainv_alg_,psb_ainv_llk_,info) + end select case default end select diff --git a/prec/impl/psb_d_sparsify.f90 b/prec/impl/psb_d_sparsify.f90 index 264f9157..f3760b93 100644 --- a/prec/impl/psb_d_sparsify.f90 +++ b/prec/impl/psb_d_sparsify.f90 @@ -63,7 +63,7 @@ ! POSSIBILITY OF SUCH DAMAGE. ! ! -subroutine amg_d_sparsify(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,info,istart,iheap,ikr) +subroutine psb_d_sparsify(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,info,istart,iheap,ikr) use psb_base_mod implicit none @@ -177,10 +177,10 @@ subroutine amg_d_sparsify(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,info,istart,ihe return -end subroutine amg_d_sparsify +end subroutine psb_d_sparsify -subroutine amg_d_sparsify_list(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,lhead,listv,ikr,info) +subroutine psb_d_sparsify_list(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,lhead,listv,ikr,info) use psb_base_mod implicit none @@ -258,4 +258,4 @@ subroutine amg_d_sparsify_list(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,lhead,list return -end subroutine amg_d_sparsify_list +end subroutine psb_d_sparsify_list diff --git a/prec/impl/psb_drwclip.f90 b/prec/impl/psb_drwclip.f90 index 528bde71..97aea428 100644 --- a/prec/impl/psb_drwclip.f90 +++ b/prec/impl/psb_drwclip.f90 @@ -63,7 +63,7 @@ ! POSSIBILITY OF SUCH DAMAGE. ! ! -subroutine amg_d_rwclip(nz,ia,ja,val,imin,imax,jmin,jmax) +subroutine psb_d_rwclip(nz,ia,ja,val,imin,imax,jmin,jmax) use psb_base_mod implicit none @@ -87,4 +87,4 @@ subroutine amg_d_rwclip(nz,ia,ja,val,imin,imax,jmin,jmax) end if end do nz = j -end subroutine amg_d_rwclip +end subroutine psb_d_rwclip diff --git a/prec/impl/psb_dsparse_biconjg_llk.F90 b/prec/impl/psb_dsparse_biconjg_llk.F90 index 63e4fa49..b0fa62a9 100644 --- a/prec/impl/psb_dsparse_biconjg_llk.F90 +++ b/prec/impl/psb_dsparse_biconjg_llk.F90 @@ -34,7 +34,7 @@ ! subroutine psb_dsparse_biconjg_llk(n,a,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod - use psb_ainv_tools_mod + use psb_d_ainv_tools_mod use psb_d_biconjg_mod, psb_protect_name => psb_dsparse_biconjg_llk ! ! Left-looking variant @@ -224,6 +224,7 @@ subroutine psb_dsparse_biconjg_llk(n,a,p,z,w,nzrmax,sp_thresh,info) ! ! Sparsify current ZVAL and put into ZMAT ! + write(psb_out_unit,'("z(1) = ",f16.14)') zval(1) call sparsify(i,nzrmax,sp_thresh,n,zval,nzrz,ia,val,info,iheap=heap,ikr=izkr) if (info /= psb_success_) then info = psb_err_internal_error_ diff --git a/prec/impl/psb_dsparse_biconjg_llk_noth.F90 b/prec/impl/psb_dsparse_biconjg_llk_noth.F90 index 2a946bdf..447cb68f 100644 --- a/prec/impl/psb_dsparse_biconjg_llk_noth.F90 +++ b/prec/impl/psb_dsparse_biconjg_llk_noth.F90 @@ -34,7 +34,7 @@ ! subroutine psb_dsparse_biconjg_llk_noth(n,a,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod - use psb_ainv_tools_mod + use psb_d_ainv_tools_mod use psb_d_biconjg_mod, psb_protect_name => psb_dsparse_biconjg_llk_noth ! diff --git a/prec/impl/psb_dsparse_biconjg_mlk.F90 b/prec/impl/psb_dsparse_biconjg_mlk.F90 index c52b2224..aae56a5a 100644 --- a/prec/impl/psb_dsparse_biconjg_mlk.F90 +++ b/prec/impl/psb_dsparse_biconjg_mlk.F90 @@ -34,7 +34,7 @@ ! subroutine psb_dsparse_biconjg_mlk(n,a,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod - use psb_ainv_tools_mod + use psb_d_ainv_tools_mod use psb_d_biconjg_mod, psb_protect_name => psb_dsparse_biconjg_mlk ! ! Left-looking variant diff --git a/prec/impl/psb_dsparse_biconjg_s_ft_llk.F90 b/prec/impl/psb_dsparse_biconjg_s_ft_llk.F90 index 6318afdc..fbf4bc02 100644 --- a/prec/impl/psb_dsparse_biconjg_s_ft_llk.F90 +++ b/prec/impl/psb_dsparse_biconjg_s_ft_llk.F90 @@ -34,7 +34,7 @@ ! subroutine psb_dsparse_biconjg_s_ft_llk(n,a,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod - use psb_ainv_tools_mod + use psb_d_ainv_tools_mod use psb_d_biconjg_mod, psb_protect_name => psb_dsparse_biconjg_s_ft_llk ! diff --git a/prec/impl/psb_dsparse_biconjg_s_llk.F90 b/prec/impl/psb_dsparse_biconjg_s_llk.F90 index fbc8891b..72257d44 100644 --- a/prec/impl/psb_dsparse_biconjg_s_llk.F90 +++ b/prec/impl/psb_dsparse_biconjg_s_llk.F90 @@ -34,7 +34,7 @@ ! subroutine psb_dsparse_biconjg_s_llk(n,a,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod - use psb_ainv_tools_mod + use psb_d_ainv_tools_mod use psb_d_biconjg_mod, psb_protect_name => psb_dsparse_biconjg_s_llk ! diff --git a/prec/impl/psb_s_ainv_bld.f90 b/prec/impl/psb_s_ainv_bld.f90 new file mode 100644 index 00000000..fe844ce3 --- /dev/null +++ b/prec/impl/psb_s_ainv_bld.f90 @@ -0,0 +1,224 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! 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 PSBLAS 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 PSBLAS 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. +! +! Moved here from AMG-AINV, original copyright below. +! +! +! AMG-AINV: Approximate Inverse plugin for +! AMG4PSBLAS version 1.0 +! +! (C) Copyright 2020 +! +! Salvatore Filippone University of Rome Tor Vergata +! +! 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 AMG4PSBLAS 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 AMG4PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +subroutine psb_s_ainv_bld(a,alg,fillin,thresh,wmat,d,zmat,desc,info,blck,iscale) + + use psb_base_mod + use psb_prec_const_mod + use psb_s_biconjg_mod + + implicit none + + ! Arguments + type(psb_sspmat_type), intent(in), target :: a + integer(psb_ipk_), intent(in) :: fillin,alg + real(psb_spk_), intent(in) :: thresh + type(psb_sspmat_type), intent(inout) :: wmat, zmat + real(psb_spk_), allocatable :: d(:) + Type(psb_desc_type), Intent(in) :: desc + integer(psb_ipk_), intent(out) :: info + type(psb_sspmat_type), intent(in), optional :: blck + integer(psb_ipk_), intent(in), optional :: iscale + ! + integer(psb_ipk_) :: i, nztota, err_act, n_row, nrow_a + type(psb_s_coo_sparse_mat) :: acoo + type(psb_s_csr_sparse_mat) :: acsr + type(psb_sspmat_type) :: atmp + real(psb_spk_), allocatable :: arws(:), acls(:) + real(psb_spk_), allocatable :: pq(:), ad(:) + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: nzrmax, iscale_ + real(psb_spk_) :: sp_thresh + real(psb_spk_) :: weight + character(len=20) :: name, ch_err + + + info = psb_success_ + name = 'psb_sainv_bld' + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = psb_cd_get_context(desc) + call psb_info(ictxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + iscale_ = psb_ilu_scale_none_ + if (present(iscale)) iscale_ = iscale + weight = sone + ! + ! Check the memory available to hold the W and Z factors + ! and allocate it if needed + ! + nrow_a = a%get_nrows() + nztota = a%get_nzeros() + + if (present(blck)) then + nztota = nztota + blck%get_nzeros() + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & ': out get_nnzeros',nrow_a,nztota,& + & a%get_nrows(),a%get_ncols(),a%get_nzeros() + + + n_row = desc%get_local_rows() + allocate(pq(n_row),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + nzrmax = fillin + sp_thresh = thresh + + ! + ! Ok, let's start first with Z (i.e. Upper) + ! + call a%csclip(acoo,info,imax=n_row,jmax=n_row) + call acsr%mv_from_coo(acoo,info) + select case(iscale_) + case(psb_ilu_scale_none_) + ! Ok, do nothing. + + case(psb_ilu_scale_maxval_) + weight = acsr%maxval() + weight = sone/weight + call acsr%scal(weight,info) + + case(psb_ilu_scale_arcsum_) + allocate(arws(n_row),acls(n_row),ad(n_row),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + call acsr%arwsum(arws) + call acsr%aclsum(acls) + ad(1:n_row) = sqrt(sqrt(arws(1:n_row)*acls(1:n_row))) + ad(1:n_row) = sone/ad(1:n_row) + call acsr%scal(ad,info,side='L') + call acsr%scal(ad,info,side='R') + case default + call psb_errpush(psb_err_from_subroutine_,name,a_err='wrong iscale') + goto 9999 + end select + + ! + ! Here for the actual workhorses. + ! Only biconjg is surviving for now.... + ! + call psb_sparse_biconjg(alg,n_row,acsr,pq,& + & zmat,wmat,nzrmax,sp_thresh,info) + ! Done. Hopefully.... + + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='sparse_biconjg') + goto 9999 + end if + call atmp%mv_from(acsr) + + ! + ! Is this right??? + ! + do i=1, n_row + if (abs(pq(i)) < d_epstol) then + pq(i) = sone + else + pq(i) = sone/pq(i) + end if + end do + + select case(iscale_) + case(psb_ilu_scale_none_) + ! Ok, do nothing. + case(psb_ilu_scale_maxval_) + pq(:) = pq(:)*weight + + case(psb_ilu_scale_arcsum_) + call zmat%scal(ad,info,side='L') + call wmat%scal(ad,info,side='R') + + case default + call psb_errpush(psb_err_from_subroutine_,name,a_err='wrong iscale') + goto 9999 + end select + + call psb_move_alloc(pq,d,info) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_s_ainv_bld diff --git a/prec/impl/psb_s_bjacprec_impl.f90 b/prec/impl/psb_s_bjacprec_impl.f90 index b2545890..cd33b26c 100644 --- a/prec/impl/psb_s_bjacprec_impl.f90 +++ b/prec/impl/psb_s_bjacprec_impl.f90 @@ -195,6 +195,46 @@ subroutine psb_s_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) goto 9999 end if + case(psb_f_ainv_) + + select case(trans_) + case('N') + call psb_spmm(sone,prec%av(psb_l_pr_),x,szero,wv,desc_data,info,& + & trans=trans_, work=aux) + + call wv1%mlt(sone,prec%dv,wv,szero,info) + + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,& + & beta,y,desc_data,info,& + & trans=trans_, work=aux) + + case('T') + call psb_spmm(sone,prec%av(psb_u_pr_),x,szero,wv,desc_data,info,& + & trans=trans_, work=aux) + + call wv1%mlt(sone,prec%dv,wv,szero,info) + + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,& + & beta,y,desc_data,info,& + & trans=trans_,work=aux) + + case('C') + + call psb_spmm(sone,prec%av(psb_u_pr_),x,szero,wv,desc_data,info,& + & trans=trans_,work=aux) + + call wv1%mlt(sone,prec%dv,wv,szero,info,conjgx=trans_) + + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,& + & beta,y,desc_data,info,& + & trans=trans_,work=aux) + + end select + if (info /= psb_success_) then + ch_err="psb_spsm" + goto 9999 + end if + case default info = psb_err_internal_error_ call psb_errpush(info,name,a_err='Invalid factorization') @@ -242,6 +282,7 @@ subroutine psb_s_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) ! Local variables integer(psb_ipk_) :: n_row,n_col real(psb_spk_), pointer :: ww(:), aux(:) + type(psb_d_vect_type) :: tx,ty integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act, ierr(5) integer(psb_ipk_) :: debug_level, debug_unit @@ -346,6 +387,29 @@ subroutine psb_s_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) goto 9999 end if + case(psb_f_ainv_) + ! Application of approximate inverse preconditioner, just some spmm + ! call psb_spmm(alpha, a, x, beta, y,desc_a, info, & + ! & trans, work) + + + select case(trans_) + case('N','T') + call psb_spmm(sone,prec%av(psb_l_pr_),x,szero,ww,desc_data,info,& + & trans=trans_, work=aux) + ww(1:n_row) = ww(1:n_row)*prec%dv%v%v(1:n_row) + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),ww,& + & beta,y,desc_data,info, trans=trans_, work=aux) + + case('C') + call psb_spmm(sone,prec%av(psb_l_pr_),x,szero,ww,desc_data,info,& + & trans=trans_, work=aux) + ww(1:n_row) = ww(1:n_row)*(prec%dv%v%v(1:n_row)) + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),ww,& + & beta,y,desc_data,info, trans=trans_, work=aux) + + end select + case default info = psb_err_internal_error_ @@ -402,9 +466,13 @@ subroutine psb_s_bjac_precinit(prec,info) prec%iprcparm(psb_ilu_fill_in_) = 0 prec%iprcparm(psb_ilu_ialg_) = psb_ilu_n_ prec%iprcparm(psb_ilu_scale_) = psb_ilu_scale_none_ + prec%iprcparm(psb_inv_fillin_) = 0 + prec%iprcparm(psb_ainv_alg_) = psb_ainv_llk_ + prec%rprcparm(:) = 0 prec%rprcparm(psb_fact_eps_) = 1E-1_psb_spk_ + prec%rprcparm(psb_inv_thresh_) = 1E-1_psb_spk_ call psb_erractionrestore(err_act) @@ -420,6 +488,7 @@ subroutine psb_s_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) use psb_base_mod use psb_s_ilu_fact_mod + use psb_s_ainv_fact_mod use psb_s_bjacprec, psb_protect_name => psb_s_bjac_precbld Implicit None @@ -432,12 +501,12 @@ subroutine psb_s_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) class(psb_i_base_vect_type), intent(in), optional :: imold ! .. Local Scalars .. - integer(psb_ipk_) :: i, m, ialg, fill_in, iscale + integer(psb_ipk_) :: i, m, ialg, fill_in, iscale, inv_fill, iinvalg integer(psb_ipk_) :: ierr(5) character :: trans, unitd type(psb_sspmat_type), allocatable :: lf, uf real(psb_spk_), allocatable :: dd(:) - real(psb_spk_) :: fact_eps + real(psb_spk_) :: fact_eps, inv_thresh integer(psb_ipk_) :: nztota, err_act, n_row, nrow_a,n_col, nhalo integer(psb_ipk_) :: ictxt,np,me character(len=20) :: name='s_bjac_precbld' @@ -468,8 +537,13 @@ subroutine psb_s_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) ! We check if all the information contained in the preconditioner structure ! are meaningful, otherwise we give an error and get out of the build ! procedure - ialg = prec%iprcparm(psb_ilu_ialg_) - if ((ialg == psb_ilu_n_).or.(ialg == psb_milu_n_).or.(ialg == psb_ilu_t_)) then + ialg = prec%iprcparm(psb_ilu_ialg_) ! Integer for ILU type algorithm + iinvalg = prec%iprcparm(psb_ainv_alg_) ! Integer for AINV type algorithm + if ((ialg == psb_ilu_n_).or.(ialg == psb_milu_n_).or.& + & (ialg == psb_ilu_t_).or.(iinvalg == psb_ainv_llk_).or.& + & (iinvalg == psb_ainv_s_llk_).or.(iinvalg == psb_ainv_s_ft_llk_).or.& + & (iinvalg == psb_ainv_llk_noth_).or.(iinvalg == psb_ainv_mlk_).or.& + & (iinvalg == psb_ainv_lmx_)) then ! Do nothing: admissible request else info=psb_err_from_subroutine_ @@ -492,24 +566,32 @@ subroutine psb_s_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) goto 9999 end if fact_eps = prec%rprcparm(psb_fact_eps_) - if(fact_eps > 1 ) then + if( (fact_eps > 1).and.(prec%iprcparm(psb_f_type_) == psb_f_ainv_) ) then info=psb_err_from_subroutine_ ch_err='psb_fact_eps_' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if - fill_in = prec%iprcparm(psb_ilu_fill_in_) - if(fill_in < 0) then + inv_thresh = prec%rprcparm(psb_inv_thresh_) + if( (inv_thresh > 1) ) then info=psb_err_from_subroutine_ - ch_err='psb_ilu_fill_in_' + ch_err='psb_fact_eps_' call psb_errpush(info,name,a_err=ch_err) goto 9999 - else if (fill_in == 0) then - ! If the requested level of fill is equal to zero, we default to the - ! specialized ILU(0) routine - prec%iprcparm(psb_f_type_) = psb_f_ilu_n_ end if - + fill_in = prec%iprcparm(psb_ilu_fill_in_) + if (prec%iprcparm(psb_f_type_) == psb_f_ilu_n_) then + if(fill_in < 0) then + info=psb_err_from_subroutine_ + ch_err='psb_ilu_fill_in_' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + else if (fill_in == 0) then + ! If the requested level of fill is equal to zero, we default to the + ! specialized ILU(0) routine + prec%iprcparm(psb_f_type_) = psb_f_ilu_n_ + end if + end if ! Select on the type of factorization to be used select case(prec%iprcparm(psb_f_type_)) @@ -732,6 +814,78 @@ subroutine psb_s_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) goto 9999 end if + case(psb_f_ainv_) + ! Approximate Inverse Factorizations based on variants of the incomplete + ! biconjugation algorithms + if (allocated(prec%av)) then + if (size(prec%av) < psb_bp_ilu_avsz) then + do i=1,size(prec%av) + call prec%av(i)%free() + enddo + deallocate(prec%av,stat=info) + endif + end if + if (.not.allocated(prec%av)) then + allocate(prec%av(psb_max_avsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + end if + endif + + nrow_a = desc_a%get_local_rows() + nztota = a%get_nzeros() + + n_col = desc_a%get_local_cols() + nhalo = n_col-nrow_a + n_row = nrow_a + + allocate(lf,uf,stat=info) + if (info == psb_success_) call lf%allocate(n_row,n_row,nztota) + if (info == psb_success_) call uf%allocate(n_row,n_row,nztota) + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(dd(n_row),stat=info) + if (info == psb_success_) then + allocate(prec%dv, stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(prec%dv%v,mold=vmold,stat=info) + else + allocate(psb_s_base_vect_type :: prec%dv%v,stat=info) + end if + end if + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + endif + ! Computin the factorization + call psb_ainv_fact(a,iinvalg,inv_fill,inv_thresh,lf,dd,uf,desc_a,info,iscale=iscale) + + if(info == psb_success_) then + call prec%av(psb_l_pr_)%mv_from(lf%a) + call prec%av(psb_u_pr_)%mv_from(uf%a) + call prec%av(psb_l_pr_)%set_asb() + call prec%av(psb_u_pr_)%set_asb() + call prec%av(psb_l_pr_)%trim() + call prec%av(psb_u_pr_)%trim() + call prec%dv%bld(dd) + ! call move_alloc(dd,prec%d) + else + info=psb_err_from_subroutine_ + ch_err='psb_ilut_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + case(psb_f_none_) info=psb_err_from_subroutine_ ch_err='Inconsistent prec psb_f_none_' @@ -792,6 +946,12 @@ subroutine psb_s_bjac_precseti(prec,what,val,info) case (psb_ilu_scale_) prec%iprcparm(psb_ilu_scale_) = val + case (psb_ainv_alg_) + prec%iprcparm(psb_ainv_alg_) = val + + case(psb_inv_fillin_) + prec%iprcparm(psb_inv_fillin_) = val + case default write(psb_err_unit,'(i0," is invalid, ignoring user specification")') what @@ -833,6 +993,9 @@ subroutine psb_s_bjac_precsetr(prec,what,val,info) case (psb_fact_eps_) prec%rprcparm(psb_fact_eps_) = val + case (psb_inv_thresh_) + prec%rprcparm(psb_inv_thresh_) = val + case default write(psb_err_unit,'(i0," is invalid, ignoring user specification")') what diff --git a/prec/impl/psb_s_prec_type_impl.f90 b/prec/impl/psb_s_prec_type_impl.f90 index 837c3106..b06f51ff 100644 --- a/prec/impl/psb_s_prec_type_impl.f90 +++ b/prec/impl/psb_s_prec_type_impl.f90 @@ -356,6 +356,8 @@ subroutine psb_scprecseti(prec,what,val,info,ilev,ilmax,pos,idx) select case (psb_toupper(what)) case ("SUB_FILLIN") call prec%prec%precset(psb_ilu_fill_in_,val,info) + case('INV_FILLIN') + call prec%prec%precset(psb_inv_fillin_,val,info) case default info = psb_err_invalid_args_combination_ write(psb_err_unit,*) name,& @@ -390,6 +392,8 @@ subroutine psb_scprecsetr(prec,what,val,info,ilev,ilmax,pos,idx) select case (psb_toupper(what)) case('SUB_ILUTHRS') call prec%prec%precset(psb_fact_eps_,val,info) + case('INV_THRESH') + call prec%prec%precset(psb_inv_thresh_,val,info) case default info = psb_err_invalid_args_combination_ write(psb_err_unit,*) name,& @@ -431,6 +435,8 @@ subroutine psb_scprecsetc(prec,what,string,info,ilev,ilmax,pos,idx) case("ILUT") call prec%prec%precset(psb_f_type_,psb_f_ilu_t_,info) call prec%prec%precset(psb_ilu_ialg_,psb_ilu_t_,info) + case("AINV") + call prec%prec%precset(psb_f_type_,psb_f_ainv_,info) case default ! Default to ILU(0) factorization call prec%prec%precset(psb_f_type_,psb_f_ilu_n_,info) @@ -450,6 +456,19 @@ subroutine psb_scprecsetc(prec,what,string,info,ilev,ilmax,pos,idx) case default call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_none_,info) end select + case ("AINV_ALG") + select case (psb_toupper(string)) + case("LLK") + call prec%prec%precset(psb_ainv_alg_,psb_ainv_llk_,info) + case("SYM-LLK") + call prec%prec%precset(psb_ainv_alg_,psb_ainv_s_llk_,info) + case("STAB-LLK") + call prec%prec%precset(psb_ainv_alg_,psb_ainv_s_ft_llk_,info) + case("MLK","LMX") + call prec%prec%precset(psb_ainv_alg_,psb_ainv_mlk_,info) + case default + call prec%prec%precset(psb_ainv_alg_,psb_ainv_llk_,info) + end select case default end select diff --git a/prec/impl/psb_s_sparsify.f90 b/prec/impl/psb_s_sparsify.f90 index 191b1da5..b271d8a8 100644 --- a/prec/impl/psb_s_sparsify.f90 +++ b/prec/impl/psb_s_sparsify.f90 @@ -63,7 +63,7 @@ ! POSSIBILITY OF SUCH DAMAGE. ! ! -subroutine amg_s_sparsify(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,info,istart,iheap,ikr) +subroutine psb_s_sparsify(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,info,istart,iheap,ikr) use psb_base_mod implicit none @@ -177,10 +177,10 @@ subroutine amg_s_sparsify(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,info,istart,ihe return -end subroutine amg_s_sparsify +end subroutine psb_s_sparsify -subroutine amg_s_sparsify_list(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,lhead,listv,ikr,info) +subroutine psb_s_sparsify_list(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,lhead,listv,ikr,info) use psb_base_mod implicit none @@ -258,4 +258,4 @@ subroutine amg_s_sparsify_list(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,lhead,list return -end subroutine amg_s_sparsify_list +end subroutine psb_s_sparsify_list diff --git a/prec/impl/psb_srwclip.f90 b/prec/impl/psb_srwclip.f90 index d9c303dd..f57207d7 100644 --- a/prec/impl/psb_srwclip.f90 +++ b/prec/impl/psb_srwclip.f90 @@ -63,7 +63,7 @@ ! POSSIBILITY OF SUCH DAMAGE. ! ! -subroutine amg_s_rwclip(nz,ia,ja,val,imin,imax,jmin,jmax) +subroutine psb_s_rwclip(nz,ia,ja,val,imin,imax,jmin,jmax) use psb_base_mod implicit none @@ -87,4 +87,4 @@ subroutine amg_s_rwclip(nz,ia,ja,val,imin,imax,jmin,jmax) end if end do nz = j -end subroutine amg_s_rwclip +end subroutine psb_s_rwclip diff --git a/prec/impl/psb_ssparse_biconjg_llk.F90 b/prec/impl/psb_ssparse_biconjg_llk.F90 index 96274016..324ebe9b 100644 --- a/prec/impl/psb_ssparse_biconjg_llk.F90 +++ b/prec/impl/psb_ssparse_biconjg_llk.F90 @@ -34,7 +34,7 @@ ! subroutine psb_ssparse_biconjg_llk(n,a,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod - use psb_ainv_tools_mod + use psb_s_ainv_tools_mod use psb_s_biconjg_mod, psb_protect_name => psb_ssparse_biconjg_llk ! ! Left-looking variant @@ -224,6 +224,7 @@ subroutine psb_ssparse_biconjg_llk(n,a,p,z,w,nzrmax,sp_thresh,info) ! ! Sparsify current ZVAL and put into ZMAT ! + write(psb_out_unit,'("z(1) = ",f16.14)') zval(1) call sparsify(i,nzrmax,sp_thresh,n,zval,nzrz,ia,val,info,iheap=heap,ikr=izkr) if (info /= psb_success_) then info = psb_err_internal_error_ diff --git a/prec/impl/psb_ssparse_biconjg_llk_noth.F90 b/prec/impl/psb_ssparse_biconjg_llk_noth.F90 index d0b73d19..0683750a 100644 --- a/prec/impl/psb_ssparse_biconjg_llk_noth.F90 +++ b/prec/impl/psb_ssparse_biconjg_llk_noth.F90 @@ -34,7 +34,7 @@ ! subroutine psb_ssparse_biconjg_llk_noth(n,a,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod - use psb_ainv_tools_mod + use psb_s_ainv_tools_mod use psb_s_biconjg_mod, psb_protect_name => psb_ssparse_biconjg_llk_noth ! diff --git a/prec/impl/psb_ssparse_biconjg_mlk.F90 b/prec/impl/psb_ssparse_biconjg_mlk.F90 index 8d5db1b2..7fc2db48 100644 --- a/prec/impl/psb_ssparse_biconjg_mlk.F90 +++ b/prec/impl/psb_ssparse_biconjg_mlk.F90 @@ -34,7 +34,7 @@ ! subroutine psb_ssparse_biconjg_mlk(n,a,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod - use psb_ainv_tools_mod + use psb_s_ainv_tools_mod use psb_s_biconjg_mod, psb_protect_name => psb_ssparse_biconjg_mlk ! ! Left-looking variant diff --git a/prec/impl/psb_ssparse_biconjg_s_ft_llk.F90 b/prec/impl/psb_ssparse_biconjg_s_ft_llk.F90 index 581edfc9..e8287e84 100644 --- a/prec/impl/psb_ssparse_biconjg_s_ft_llk.F90 +++ b/prec/impl/psb_ssparse_biconjg_s_ft_llk.F90 @@ -34,7 +34,7 @@ ! subroutine psb_ssparse_biconjg_s_ft_llk(n,a,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod - use psb_ainv_tools_mod + use psb_s_ainv_tools_mod use psb_s_biconjg_mod, psb_protect_name => psb_ssparse_biconjg_s_ft_llk ! @@ -164,7 +164,7 @@ subroutine psb_ssparse_biconjg_s_ft_llk(n,a,p,z,w,nzrmax,sp_thresh,info) ip2 = w%icp(j+1) - 1 nzra = max(0,ip2 - ip1 + 1) nzww = 0 - call psb_d_spvspm(sone,a,nzra,w%ia(ip1:ip2),w%val(ip1:ip2),& + call psb_s_spvspm(sone,a,nzra,w%ia(ip1:ip2),w%val(ip1:ip2),& & szero,nzww,iww,ww,info) p(i) = psb_spge_dot(nzww,iww,ww,zval) @@ -299,7 +299,7 @@ subroutine psb_ssparse_biconjg_s_ft_llk(n,a,p,z,w,nzrmax,sp_thresh,info) ip2 = z%icp(j+1) - 1 nzra = max(0,ip2 - ip1 + 1) nzww = 0 - call psb_d_spmspv(sone,ac,nzra,z%ia(ip1:ip2),z%val(ip1:ip2),& + call psb_s_spmspv(sone,ac,nzra,z%ia(ip1:ip2),z%val(ip1:ip2),& & szero,nzww,iww,ww,info) q(i) = psb_spge_dot(nzww,iww,ww,zval) @@ -384,7 +384,7 @@ subroutine psb_ssparse_biconjg_s_ft_llk(n,a,p,z,w,nzrmax,sp_thresh,info) nzww = 0 nzrz = z%icp(i+1)-z%icp(i) ipz1 = z%icp(i) - call psb_d_spmspv(sone,ac,& + call psb_s_spmspv(sone,ac,& & nzrz,z%ia(ipz1:ipz1+nzrz-1),z%val(ipz1:ipz1+nzrz-1),& & szero,nzww,iww,ww,info) tmpq = psb_spdot_srtd(nzww,iww,ww,nzrw,ia,val) diff --git a/prec/impl/psb_ssparse_biconjg_s_llk.F90 b/prec/impl/psb_ssparse_biconjg_s_llk.F90 index 9d971d8d..42e71089 100644 --- a/prec/impl/psb_ssparse_biconjg_s_llk.F90 +++ b/prec/impl/psb_ssparse_biconjg_s_llk.F90 @@ -34,7 +34,7 @@ ! subroutine psb_ssparse_biconjg_s_llk(n,a,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod - use psb_ainv_tools_mod + use psb_s_ainv_tools_mod use psb_s_biconjg_mod, psb_protect_name => psb_ssparse_biconjg_s_llk ! diff --git a/prec/impl/psb_z_ainv_bld.f90 b/prec/impl/psb_z_ainv_bld.f90 new file mode 100644 index 00000000..b16cee6b --- /dev/null +++ b/prec/impl/psb_z_ainv_bld.f90 @@ -0,0 +1,224 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! 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 PSBLAS 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 PSBLAS 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. +! +! Moved here from AMG-AINV, original copyright below. +! +! +! AMG-AINV: Approximate Inverse plugin for +! AMG4PSBLAS version 1.0 +! +! (C) Copyright 2020 +! +! Salvatore Filippone University of Rome Tor Vergata +! +! 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 AMG4PSBLAS 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 AMG4PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +subroutine psb_z_ainv_bld(a,alg,fillin,thresh,wmat,d,zmat,desc,info,blck,iscale) + + use psb_base_mod + use psb_prec_const_mod + use psb_z_biconjg_mod + + implicit none + + ! Arguments + type(psb_zspmat_type), intent(in), target :: a + integer(psb_ipk_), intent(in) :: fillin,alg + real(psb_dpk_), intent(in) :: thresh + type(psb_zspmat_type), intent(inout) :: wmat, zmat + complex(psb_dpk_), allocatable :: d(:) + Type(psb_desc_type), Intent(in) :: desc + integer(psb_ipk_), intent(out) :: info + type(psb_zspmat_type), intent(in), optional :: blck + integer(psb_ipk_), intent(in), optional :: iscale + ! + integer(psb_ipk_) :: i, nztota, err_act, n_row, nrow_a + type(psb_z_coo_sparse_mat) :: acoo + type(psb_z_csr_sparse_mat) :: acsr + type(psb_zspmat_type) :: atmp + real(psb_dpk_), allocatable :: arws(:), acls(:) + complex(psb_dpk_), allocatable :: pq(:), ad(:) + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: nzrmax, iscale_ + real(psb_dpk_) :: sp_thresh + complex(psb_dpk_) :: weight + character(len=20) :: name, ch_err + + + info = psb_success_ + name = 'psb_zainv_bld' + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = psb_cd_get_context(desc) + call psb_info(ictxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + iscale_ = psb_ilu_scale_none_ + if (present(iscale)) iscale_ = iscale + weight = zone + ! + ! Check the memory available to hold the W and Z factors + ! and allocate it if needed + ! + nrow_a = a%get_nrows() + nztota = a%get_nzeros() + + if (present(blck)) then + nztota = nztota + blck%get_nzeros() + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & ': out get_nnzeros',nrow_a,nztota,& + & a%get_nrows(),a%get_ncols(),a%get_nzeros() + + + n_row = desc%get_local_rows() + allocate(pq(n_row),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + nzrmax = fillin + sp_thresh = thresh + + ! + ! Ok, let's start first with Z (i.e. Upper) + ! + call a%csclip(acoo,info,imax=n_row,jmax=n_row) + call acsr%mv_from_coo(acoo,info) + select case(iscale_) + case(psb_ilu_scale_none_) + ! Ok, do nothing. + + case(psb_ilu_scale_maxval_) + weight = acsr%maxval() + weight = zone/weight + call acsr%scal(weight,info) + + case(psb_ilu_scale_arcsum_) + allocate(arws(n_row),acls(n_row),ad(n_row),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + call acsr%arwsum(arws) + call acsr%aclsum(acls) + ad(1:n_row) = sqrt(sqrt(arws(1:n_row)*acls(1:n_row))) + ad(1:n_row) = zone/ad(1:n_row) + call acsr%scal(ad,info,side='L') + call acsr%scal(ad,info,side='R') + case default + call psb_errpush(psb_err_from_subroutine_,name,a_err='wrong iscale') + goto 9999 + end select + + ! + ! Here for the actual workhorses. + ! Only biconjg is surviving for now.... + ! + call psb_sparse_biconjg(alg,n_row,acsr,pq,& + & zmat,wmat,nzrmax,sp_thresh,info) + ! Done. Hopefully.... + + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='sparse_biconjg') + goto 9999 + end if + call atmp%mv_from(acsr) + + ! + ! Is this right??? + ! + do i=1, n_row + if (abs(pq(i)) < d_epstol) then + pq(i) = zone + else + pq(i) = zone/pq(i) + end if + end do + + select case(iscale_) + case(psb_ilu_scale_none_) + ! Ok, do nothing. + case(psb_ilu_scale_maxval_) + pq(:) = pq(:)*weight + + case(psb_ilu_scale_arcsum_) + call zmat%scal(ad,info,side='L') + call wmat%scal(ad,info,side='R') + + case default + call psb_errpush(psb_err_from_subroutine_,name,a_err='wrong iscale') + goto 9999 + end select + + call psb_move_alloc(pq,d,info) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_z_ainv_bld diff --git a/prec/impl/psb_z_bjacprec_impl.f90 b/prec/impl/psb_z_bjacprec_impl.f90 index 2a8960ea..b0650401 100644 --- a/prec/impl/psb_z_bjacprec_impl.f90 +++ b/prec/impl/psb_z_bjacprec_impl.f90 @@ -195,6 +195,46 @@ subroutine psb_z_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) goto 9999 end if + case(psb_f_ainv_) + + select case(trans_) + case('N') + call psb_spmm(zone,prec%av(psb_l_pr_),x,zzero,wv,desc_data,info,& + & trans=trans_, work=aux) + + call wv1%mlt(zone,prec%dv,wv,zzero,info) + + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,& + & beta,y,desc_data,info,& + & trans=trans_, work=aux) + + case('T') + call psb_spmm(zone,prec%av(psb_u_pr_),x,zzero,wv,desc_data,info,& + & trans=trans_, work=aux) + + call wv1%mlt(zone,prec%dv,wv,zzero,info) + + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,& + & beta,y,desc_data,info,& + & trans=trans_,work=aux) + + case('C') + + call psb_spmm(zone,prec%av(psb_u_pr_),x,zzero,wv,desc_data,info,& + & trans=trans_,work=aux) + + call wv1%mlt(zone,prec%dv,wv,zzero,info,conjgx=trans_) + + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,& + & beta,y,desc_data,info,& + & trans=trans_,work=aux) + + end select + if (info /= psb_success_) then + ch_err="psb_spsm" + goto 9999 + end if + case default info = psb_err_internal_error_ call psb_errpush(info,name,a_err='Invalid factorization') @@ -242,6 +282,7 @@ subroutine psb_z_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) ! Local variables integer(psb_ipk_) :: n_row,n_col complex(psb_dpk_), pointer :: ww(:), aux(:) + type(psb_d_vect_type) :: tx,ty integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act, ierr(5) integer(psb_ipk_) :: debug_level, debug_unit @@ -346,6 +387,29 @@ subroutine psb_z_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) goto 9999 end if + case(psb_f_ainv_) + ! Application of approximate inverse preconditioner, just some spmm + ! call psb_spmm(alpha, a, x, beta, y,desc_a, info, & + ! & trans, work) + + + select case(trans_) + case('N','T') + call psb_spmm(zone,prec%av(psb_l_pr_),x,zzero,ww,desc_data,info,& + & trans=trans_, work=aux) + ww(1:n_row) = ww(1:n_row)*prec%dv%v%v(1:n_row) + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),ww,& + & beta,y,desc_data,info, trans=trans_, work=aux) + + case('C') + call psb_spmm(zone,prec%av(psb_l_pr_),x,zzero,ww,desc_data,info,& + & trans=trans_, work=aux) + ww(1:n_row) = ww(1:n_row)*conjg(prec%dv%v%v(1:n_row)) + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),ww,& + & beta,y,desc_data,info, trans=trans_, work=aux) + + end select + case default info = psb_err_internal_error_ @@ -402,9 +466,13 @@ subroutine psb_z_bjac_precinit(prec,info) prec%iprcparm(psb_ilu_fill_in_) = 0 prec%iprcparm(psb_ilu_ialg_) = psb_ilu_n_ prec%iprcparm(psb_ilu_scale_) = psb_ilu_scale_none_ + prec%iprcparm(psb_inv_fillin_) = 0 + prec%iprcparm(psb_ainv_alg_) = psb_ainv_llk_ + prec%rprcparm(:) = 0 prec%rprcparm(psb_fact_eps_) = 1E-1_psb_dpk_ + prec%rprcparm(psb_inv_thresh_) = 1E-1_psb_dpk_ call psb_erractionrestore(err_act) @@ -420,6 +488,7 @@ subroutine psb_z_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) use psb_base_mod use psb_z_ilu_fact_mod + use psb_z_ainv_fact_mod use psb_z_bjacprec, psb_protect_name => psb_z_bjac_precbld Implicit None @@ -432,12 +501,12 @@ subroutine psb_z_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) class(psb_i_base_vect_type), intent(in), optional :: imold ! .. Local Scalars .. - integer(psb_ipk_) :: i, m, ialg, fill_in, iscale + integer(psb_ipk_) :: i, m, ialg, fill_in, iscale, inv_fill, iinvalg integer(psb_ipk_) :: ierr(5) character :: trans, unitd type(psb_zspmat_type), allocatable :: lf, uf complex(psb_dpk_), allocatable :: dd(:) - real(psb_dpk_) :: fact_eps + real(psb_dpk_) :: fact_eps, inv_thresh integer(psb_ipk_) :: nztota, err_act, n_row, nrow_a,n_col, nhalo integer(psb_ipk_) :: ictxt,np,me character(len=20) :: name='z_bjac_precbld' @@ -468,8 +537,13 @@ subroutine psb_z_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) ! We check if all the information contained in the preconditioner structure ! are meaningful, otherwise we give an error and get out of the build ! procedure - ialg = prec%iprcparm(psb_ilu_ialg_) - if ((ialg == psb_ilu_n_).or.(ialg == psb_milu_n_).or.(ialg == psb_ilu_t_)) then + ialg = prec%iprcparm(psb_ilu_ialg_) ! Integer for ILU type algorithm + iinvalg = prec%iprcparm(psb_ainv_alg_) ! Integer for AINV type algorithm + if ((ialg == psb_ilu_n_).or.(ialg == psb_milu_n_).or.& + & (ialg == psb_ilu_t_).or.(iinvalg == psb_ainv_llk_).or.& + & (iinvalg == psb_ainv_s_llk_).or.(iinvalg == psb_ainv_s_ft_llk_).or.& + & (iinvalg == psb_ainv_llk_noth_).or.(iinvalg == psb_ainv_mlk_).or.& + & (iinvalg == psb_ainv_lmx_)) then ! Do nothing: admissible request else info=psb_err_from_subroutine_ @@ -492,24 +566,32 @@ subroutine psb_z_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) goto 9999 end if fact_eps = prec%rprcparm(psb_fact_eps_) - if(fact_eps > 1 ) then + if( (fact_eps > 1).and.(prec%iprcparm(psb_f_type_) == psb_f_ainv_) ) then info=psb_err_from_subroutine_ ch_err='psb_fact_eps_' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if - fill_in = prec%iprcparm(psb_ilu_fill_in_) - if(fill_in < 0) then + inv_thresh = prec%rprcparm(psb_inv_thresh_) + if( (inv_thresh > 1) ) then info=psb_err_from_subroutine_ - ch_err='psb_ilu_fill_in_' + ch_err='psb_fact_eps_' call psb_errpush(info,name,a_err=ch_err) goto 9999 - else if (fill_in == 0) then - ! If the requested level of fill is equal to zero, we default to the - ! specialized ILU(0) routine - prec%iprcparm(psb_f_type_) = psb_f_ilu_n_ end if - + fill_in = prec%iprcparm(psb_ilu_fill_in_) + if (prec%iprcparm(psb_f_type_) == psb_f_ilu_n_) then + if(fill_in < 0) then + info=psb_err_from_subroutine_ + ch_err='psb_ilu_fill_in_' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + else if (fill_in == 0) then + ! If the requested level of fill is equal to zero, we default to the + ! specialized ILU(0) routine + prec%iprcparm(psb_f_type_) = psb_f_ilu_n_ + end if + end if ! Select on the type of factorization to be used select case(prec%iprcparm(psb_f_type_)) @@ -732,6 +814,78 @@ subroutine psb_z_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) goto 9999 end if + case(psb_f_ainv_) + ! Approximate Inverse Factorizations based on variants of the incomplete + ! biconjugation algorithms + if (allocated(prec%av)) then + if (size(prec%av) < psb_bp_ilu_avsz) then + do i=1,size(prec%av) + call prec%av(i)%free() + enddo + deallocate(prec%av,stat=info) + endif + end if + if (.not.allocated(prec%av)) then + allocate(prec%av(psb_max_avsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + end if + endif + + nrow_a = desc_a%get_local_rows() + nztota = a%get_nzeros() + + n_col = desc_a%get_local_cols() + nhalo = n_col-nrow_a + n_row = nrow_a + + allocate(lf,uf,stat=info) + if (info == psb_success_) call lf%allocate(n_row,n_row,nztota) + if (info == psb_success_) call uf%allocate(n_row,n_row,nztota) + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(dd(n_row),stat=info) + if (info == psb_success_) then + allocate(prec%dv, stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(prec%dv%v,mold=vmold,stat=info) + else + allocate(psb_z_base_vect_type :: prec%dv%v,stat=info) + end if + end if + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + endif + ! Computin the factorization + call psb_ainv_fact(a,iinvalg,inv_fill,inv_thresh,lf,dd,uf,desc_a,info,iscale=iscale) + + if(info == psb_success_) then + call prec%av(psb_l_pr_)%mv_from(lf%a) + call prec%av(psb_u_pr_)%mv_from(uf%a) + call prec%av(psb_l_pr_)%set_asb() + call prec%av(psb_u_pr_)%set_asb() + call prec%av(psb_l_pr_)%trim() + call prec%av(psb_u_pr_)%trim() + call prec%dv%bld(dd) + ! call move_alloc(dd,prec%d) + else + info=psb_err_from_subroutine_ + ch_err='psb_ilut_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + case(psb_f_none_) info=psb_err_from_subroutine_ ch_err='Inconsistent prec psb_f_none_' @@ -792,6 +946,12 @@ subroutine psb_z_bjac_precseti(prec,what,val,info) case (psb_ilu_scale_) prec%iprcparm(psb_ilu_scale_) = val + case (psb_ainv_alg_) + prec%iprcparm(psb_ainv_alg_) = val + + case(psb_inv_fillin_) + prec%iprcparm(psb_inv_fillin_) = val + case default write(psb_err_unit,'(i0," is invalid, ignoring user specification")') what @@ -833,6 +993,9 @@ subroutine psb_z_bjac_precsetr(prec,what,val,info) case (psb_fact_eps_) prec%rprcparm(psb_fact_eps_) = val + case (psb_inv_thresh_) + prec%rprcparm(psb_inv_thresh_) = val + case default write(psb_err_unit,'(i0," is invalid, ignoring user specification")') what diff --git a/prec/impl/psb_z_prec_type_impl.f90 b/prec/impl/psb_z_prec_type_impl.f90 index 3f4afe66..19d38302 100644 --- a/prec/impl/psb_z_prec_type_impl.f90 +++ b/prec/impl/psb_z_prec_type_impl.f90 @@ -356,6 +356,8 @@ subroutine psb_zcprecseti(prec,what,val,info,ilev,ilmax,pos,idx) select case (psb_toupper(what)) case ("SUB_FILLIN") call prec%prec%precset(psb_ilu_fill_in_,val,info) + case('INV_FILLIN') + call prec%prec%precset(psb_inv_fillin_,val,info) case default info = psb_err_invalid_args_combination_ write(psb_err_unit,*) name,& @@ -390,6 +392,8 @@ subroutine psb_zcprecsetr(prec,what,val,info,ilev,ilmax,pos,idx) select case (psb_toupper(what)) case('SUB_ILUTHRS') call prec%prec%precset(psb_fact_eps_,val,info) + case('INV_THRESH') + call prec%prec%precset(psb_inv_thresh_,val,info) case default info = psb_err_invalid_args_combination_ write(psb_err_unit,*) name,& @@ -431,6 +435,8 @@ subroutine psb_zcprecsetc(prec,what,string,info,ilev,ilmax,pos,idx) case("ILUT") call prec%prec%precset(psb_f_type_,psb_f_ilu_t_,info) call prec%prec%precset(psb_ilu_ialg_,psb_ilu_t_,info) + case("AINV") + call prec%prec%precset(psb_f_type_,psb_f_ainv_,info) case default ! Default to ILU(0) factorization call prec%prec%precset(psb_f_type_,psb_f_ilu_n_,info) @@ -450,6 +456,19 @@ subroutine psb_zcprecsetc(prec,what,string,info,ilev,ilmax,pos,idx) case default call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_none_,info) end select + case ("AINV_ALG") + select case (psb_toupper(string)) + case("LLK") + call prec%prec%precset(psb_ainv_alg_,psb_ainv_llk_,info) + case("SYM-LLK") + call prec%prec%precset(psb_ainv_alg_,psb_ainv_s_llk_,info) + case("STAB-LLK") + call prec%prec%precset(psb_ainv_alg_,psb_ainv_s_ft_llk_,info) + case("MLK","LMX") + call prec%prec%precset(psb_ainv_alg_,psb_ainv_mlk_,info) + case default + call prec%prec%precset(psb_ainv_alg_,psb_ainv_llk_,info) + end select case default end select diff --git a/prec/impl/psb_z_sparsify.f90 b/prec/impl/psb_z_sparsify.f90 index e19ef19e..9bfb2428 100644 --- a/prec/impl/psb_z_sparsify.f90 +++ b/prec/impl/psb_z_sparsify.f90 @@ -63,7 +63,7 @@ ! POSSIBILITY OF SUCH DAMAGE. ! ! -subroutine amg_z_sparsify(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,info,istart,iheap,ikr) +subroutine psb_z_sparsify(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,info,istart,iheap,ikr) use psb_base_mod implicit none @@ -177,10 +177,10 @@ subroutine amg_z_sparsify(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,info,istart,ihe return -end subroutine amg_z_sparsify +end subroutine psb_z_sparsify -subroutine amg_z_sparsify_list(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,lhead,listv,ikr,info) +subroutine psb_z_sparsify_list(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,lhead,listv,ikr,info) use psb_base_mod implicit none @@ -258,4 +258,4 @@ subroutine amg_z_sparsify_list(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,lhead,list return -end subroutine amg_z_sparsify_list +end subroutine psb_z_sparsify_list diff --git a/prec/impl/psb_zrwclip.f90 b/prec/impl/psb_zrwclip.f90 index 41de9603..574ebcf8 100644 --- a/prec/impl/psb_zrwclip.f90 +++ b/prec/impl/psb_zrwclip.f90 @@ -63,7 +63,7 @@ ! POSSIBILITY OF SUCH DAMAGE. ! ! -subroutine amg_z_rwclip(nz,ia,ja,val,imin,imax,jmin,jmax) +subroutine psb_z_rwclip(nz,ia,ja,val,imin,imax,jmin,jmax) use psb_base_mod implicit none @@ -87,4 +87,4 @@ subroutine amg_z_rwclip(nz,ia,ja,val,imin,imax,jmin,jmax) end if end do nz = j -end subroutine amg_z_rwclip +end subroutine psb_z_rwclip diff --git a/prec/impl/psb_zsparse_biconjg_llk.F90 b/prec/impl/psb_zsparse_biconjg_llk.F90 index 0d7c9a54..7913ab85 100644 --- a/prec/impl/psb_zsparse_biconjg_llk.F90 +++ b/prec/impl/psb_zsparse_biconjg_llk.F90 @@ -34,7 +34,7 @@ ! subroutine psb_zsparse_biconjg_llk(n,a,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod - use psb_ainv_tools_mod + use psb_z_ainv_tools_mod use psb_z_biconjg_mod, psb_protect_name => psb_zsparse_biconjg_llk ! ! Left-looking variant @@ -224,6 +224,7 @@ subroutine psb_zsparse_biconjg_llk(n,a,p,z,w,nzrmax,sp_thresh,info) ! ! Sparsify current ZVAL and put into ZMAT ! + write(psb_out_unit,'("z(1) = ",f16.14)') zval(1) call sparsify(i,nzrmax,sp_thresh,n,zval,nzrz,ia,val,info,iheap=heap,ikr=izkr) if (info /= psb_success_) then info = psb_err_internal_error_ diff --git a/prec/impl/psb_zsparse_biconjg_llk_noth.F90 b/prec/impl/psb_zsparse_biconjg_llk_noth.F90 index 91fe8659..65975a24 100644 --- a/prec/impl/psb_zsparse_biconjg_llk_noth.F90 +++ b/prec/impl/psb_zsparse_biconjg_llk_noth.F90 @@ -34,7 +34,7 @@ ! subroutine psb_zsparse_biconjg_llk_noth(n,a,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod - use psb_ainv_tools_mod + use psb_z_ainv_tools_mod use psb_z_biconjg_mod, psb_protect_name => psb_zsparse_biconjg_llk_noth ! diff --git a/prec/impl/psb_zsparse_biconjg_mlk.F90 b/prec/impl/psb_zsparse_biconjg_mlk.F90 index 6e8e705f..52c30c2a 100644 --- a/prec/impl/psb_zsparse_biconjg_mlk.F90 +++ b/prec/impl/psb_zsparse_biconjg_mlk.F90 @@ -34,7 +34,7 @@ ! subroutine psb_zsparse_biconjg_mlk(n,a,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod - use psb_ainv_tools_mod + use psb_z_ainv_tools_mod use psb_z_biconjg_mod, psb_protect_name => psb_zsparse_biconjg_mlk ! ! Left-looking variant diff --git a/prec/impl/psb_zsparse_biconjg_s_ft_llk.F90 b/prec/impl/psb_zsparse_biconjg_s_ft_llk.F90 index a8f545be..541a755c 100644 --- a/prec/impl/psb_zsparse_biconjg_s_ft_llk.F90 +++ b/prec/impl/psb_zsparse_biconjg_s_ft_llk.F90 @@ -34,7 +34,7 @@ ! subroutine psb_zsparse_biconjg_s_ft_llk(n,a,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod - use psb_ainv_tools_mod + use psb_z_ainv_tools_mod use psb_z_biconjg_mod, psb_protect_name => psb_zsparse_biconjg_s_ft_llk ! @@ -164,7 +164,7 @@ subroutine psb_zsparse_biconjg_s_ft_llk(n,a,p,z,w,nzrmax,sp_thresh,info) ip2 = w%icp(j+1) - 1 nzra = max(0,ip2 - ip1 + 1) nzww = 0 - call psb_d_spvspm(zone,a,nzra,w%ia(ip1:ip2),w%val(ip1:ip2),& + call psb_z_spvspm(zone,a,nzra,w%ia(ip1:ip2),w%val(ip1:ip2),& & zzero,nzww,iww,ww,info) p(i) = psb_spge_dot(nzww,iww,ww,zval) @@ -299,7 +299,7 @@ subroutine psb_zsparse_biconjg_s_ft_llk(n,a,p,z,w,nzrmax,sp_thresh,info) ip2 = z%icp(j+1) - 1 nzra = max(0,ip2 - ip1 + 1) nzww = 0 - call psb_d_spmspv(zone,ac,nzra,z%ia(ip1:ip2),z%val(ip1:ip2),& + call psb_z_spmspv(zone,ac,nzra,z%ia(ip1:ip2),z%val(ip1:ip2),& & zzero,nzww,iww,ww,info) q(i) = psb_spge_dot(nzww,iww,ww,zval) @@ -384,7 +384,7 @@ subroutine psb_zsparse_biconjg_s_ft_llk(n,a,p,z,w,nzrmax,sp_thresh,info) nzww = 0 nzrz = z%icp(i+1)-z%icp(i) ipz1 = z%icp(i) - call psb_d_spmspv(zone,ac,& + call psb_z_spmspv(zone,ac,& & nzrz,z%ia(ipz1:ipz1+nzrz-1),z%val(ipz1:ipz1+nzrz-1),& & zzero,nzww,iww,ww,info) tmpq = psb_spdot_srtd(nzww,iww,ww,nzrw,ia,val) diff --git a/prec/impl/psb_zsparse_biconjg_s_llk.F90 b/prec/impl/psb_zsparse_biconjg_s_llk.F90 index 51114753..e4d6624a 100644 --- a/prec/impl/psb_zsparse_biconjg_s_llk.F90 +++ b/prec/impl/psb_zsparse_biconjg_s_llk.F90 @@ -34,7 +34,7 @@ ! subroutine psb_zsparse_biconjg_s_llk(n,a,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod - use psb_ainv_tools_mod + use psb_z_ainv_tools_mod use psb_z_biconjg_mod, psb_protect_name => psb_zsparse_biconjg_s_llk ! diff --git a/prec/psb_base_ainv_mod.F90 b/prec/psb_c_ainv_fact_mod.f90 similarity index 64% rename from prec/psb_base_ainv_mod.F90 rename to prec/psb_c_ainv_fact_mod.f90 index 71bd5197..200ff561 100644 --- a/prec/psb_base_ainv_mod.F90 +++ b/prec/psb_c_ainv_fact_mod.f90 @@ -28,15 +28,18 @@ ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. ! -! Moved here from AMG-AINV, original copyright below. +! Moved here from MLD2P4, original copyright below. ! ! -! AMG-AINV: Approximate Inverse plugin for -! AMG4PSBLAS version 1.0 +! MLD2P4 version 2.2 +! MultiLevel Domain Decomposition Parallel Preconditioners Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) ! -! (C) Copyright 2020 +! (C) Copyright 2008-2018 ! -! Salvatore Filippone University of Rome Tor Vergata +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions @@ -46,14 +49,14 @@ ! 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 AMG4PSBLAS group or the names of its contributors may +! 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 AMG4PSBLAS GROUP OR ITS CONTRIBUTORS +! 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 @@ -64,31 +67,32 @@ ! ! ! -module psb_base_ainv_mod - - use psb_prec_mod - - integer, parameter :: psb_inv_fillin_ = 100 ! To check for compatibility - integer, parameter :: psb_ainv_alg_ = psb_inv_fillin_ + 1 - integer, parameter :: psb_inv_thresh_ = 102 ! To check for compatibility -#if 0 - integer, parameter :: psb_ainv_orth1_ = psb_inv_thresh_ + 1 - integer, parameter :: psb_ainv_orth2_ = psb_ainv_orth1_ + 1 - integer, parameter :: psb_ainv_orth3_ = psb_ainv_orth2_ + 1 - integer, parameter :: psb_ainv_orth4_ = psb_ainv_orth3_ + 1 - integer, parameter :: psb_ainv_llk_ = psb_ainv_orth4_ + 1 -#else - integer, parameter :: psb_ainv_llk_ = psb_inv_thresh_ + 1 -#endif - integer, parameter :: psb_ainv_s_llk_ = psb_ainv_llk_ + 1 - integer, parameter :: psb_ainv_s_ft_llk_ = psb_ainv_s_llk_ + 1 - integer, parameter :: psb_ainv_llk_noth_ = psb_ainv_s_ft_llk_ + 1 - integer, parameter :: psb_ainv_mlk_ = psb_ainv_llk_noth_ + 1 - integer, parameter :: psb_ainv_lmx_ = psb_ainv_mlk_ -#if defined(HAVE_TUMA_SAINV) - integer, parameter :: psb_ainv_s_tuma_ = psb_ainv_lmx_ + 1 - integer, parameter :: psb_ainv_l_tuma_ = psb_ainv_s_tuma_ + 1 -#endif +! +! File: psb_c_ainv_fact_mod.f90 +! +! Module: psb_c_ainv_fact_mod +! +! This module defines some interfaces used internally by the implementation of +! psb_c_ainv_solver, but not visible to the end user. +! +! +module psb_c_ainv_fact_mod + use psb_base_mod + use psb_prec_const_mod + interface psb_ainv_fact + subroutine psb_c_ainv_bld(a,alg,fillin,thresh,wmat,d,zmat,desc,info,blck,iscale) + import psb_cspmat_type, psb_spk_, psb_ipk_, psb_desc_type + type(psb_cspmat_type), intent(in), target :: a + integer(psb_ipk_), intent(in) :: fillin,alg + real(psb_spk_), intent(in) :: thresh + type(psb_cspmat_type), intent(inout) :: wmat, zmat + complex(psb_spk_), allocatable :: d(:) + Type(psb_desc_type), Intent(in) :: desc + integer(psb_ipk_), intent(out) :: info + type(psb_cspmat_type), intent(in), optional :: blck + integer(psb_ipk_), intent(in), optional :: iscale + end subroutine psb_c_ainv_bld + end interface -end module psb_base_ainv_mod +end module psb_c_ainv_fact_mod diff --git a/prec/psb_c_biconjg_mod.F90 b/prec/psb_c_biconjg_mod.F90 index f3292639..44fabb34 100644 --- a/prec/psb_c_biconjg_mod.F90 +++ b/prec/psb_c_biconjg_mod.F90 @@ -100,11 +100,13 @@ ! module psb_c_biconjg_mod + use psb_base_mod + use psb_prec_const_mod + interface psb_sparse_biconjg module procedure psb_csparse_biconjg end interface - abstract interface subroutine psb_csparse_biconjg_variant(n,a,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod, only : psb_c_csr_sparse_mat, psb_c_csc_sparse_mat, & @@ -131,12 +133,11 @@ module psb_c_biconjg_mod & psb_csparse_tuma_lainv #endif - contains subroutine psb_csparse_biconjg(alg,n,acsr,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod - use psb_base_ainv_mod + use psb_prec_const_mod integer(psb_ipk_), intent(in) :: alg,n type(psb_c_csr_sparse_mat), intent(in) :: acsr type(psb_cspmat_type), intent(out) :: z, w diff --git a/prec/psb_c_bjacprec.f90 b/prec/psb_c_bjacprec.f90 index ca4cbcb8..fb1d5429 100644 --- a/prec/psb_c_bjacprec.f90 +++ b/prec/psb_c_bjacprec.f90 @@ -33,6 +33,7 @@ module psb_c_bjacprec use psb_c_base_prec_mod use psb_c_ilu_fact_mod + use psb_c_ainv_fact_mod type, extends(psb_c_base_prec_type) :: psb_c_bjac_prec_type integer(psb_ipk_), allocatable :: iprcparm(:) @@ -58,8 +59,8 @@ module psb_c_bjacprec end type psb_c_bjac_prec_type character(len=15), parameter, private :: & - & fact_names(0:3)=(/'None ','ILU(0) ',& - & 'ILU(n) ','ILU(eps) '/) + & fact_names(0:4)=(/'None ','ILU(0) ',& + & 'ILU(n) ','ILU(eps) ','AINV(eps) '/) private :: psb_c_bjac_sizeof, psb_c_bjac_precdescr, psb_c_bjac_get_nzeros diff --git a/prec/psb_d_ainv_fact_mod.f90 b/prec/psb_d_ainv_fact_mod.f90 new file mode 100644 index 00000000..8eb6fbc8 --- /dev/null +++ b/prec/psb_d_ainv_fact_mod.f90 @@ -0,0 +1,98 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! 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 PSBLAS 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 PSBLAS 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. +! +! Moved here from MLD2P4, original copyright below. +! +! +! MLD2P4 version 2.2 +! MultiLevel Domain Decomposition Parallel Preconditioners Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2008-2018 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! +! 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: psb_d_ainv_fact_mod.f90 +! +! Module: psb_d_ainv_fact_mod +! +! This module defines some interfaces used internally by the implementation of +! psb_d_ainv_solver, but not visible to the end user. +! +! +module psb_d_ainv_fact_mod + use psb_base_mod + use psb_prec_const_mod + + interface psb_ainv_fact + subroutine psb_d_ainv_bld(a,alg,fillin,thresh,wmat,d,zmat,desc,info,blck,iscale) + import psb_dspmat_type, psb_dpk_, psb_ipk_, psb_desc_type + type(psb_dspmat_type), intent(in), target :: a + integer(psb_ipk_), intent(in) :: fillin,alg + real(psb_dpk_), intent(in) :: thresh + type(psb_dspmat_type), intent(inout) :: wmat, zmat + real(psb_dpk_), allocatable :: d(:) + Type(psb_desc_type), Intent(in) :: desc + integer(psb_ipk_), intent(out) :: info + type(psb_dspmat_type), intent(in), optional :: blck + integer(psb_ipk_), intent(in), optional :: iscale + end subroutine psb_d_ainv_bld + end interface + +end module psb_d_ainv_fact_mod diff --git a/prec/psb_d_biconjg_mod.F90 b/prec/psb_d_biconjg_mod.F90 index 2c68437a..0542b2ed 100644 --- a/prec/psb_d_biconjg_mod.F90 +++ b/prec/psb_d_biconjg_mod.F90 @@ -100,11 +100,13 @@ ! module psb_d_biconjg_mod + use psb_base_mod + use psb_prec_const_mod + interface psb_sparse_biconjg module procedure psb_dsparse_biconjg end interface - abstract interface subroutine psb_dsparse_biconjg_variant(n,a,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod, only : psb_d_csr_sparse_mat, psb_d_csc_sparse_mat, & @@ -131,12 +133,11 @@ module psb_d_biconjg_mod & psb_dsparse_tuma_lainv #endif - contains subroutine psb_dsparse_biconjg(alg,n,acsr,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod - use psb_base_ainv_mod + use psb_prec_const_mod integer(psb_ipk_), intent(in) :: alg,n type(psb_d_csr_sparse_mat), intent(in) :: acsr type(psb_dspmat_type), intent(out) :: z, w diff --git a/prec/psb_d_bjacprec.f90 b/prec/psb_d_bjacprec.f90 index 0452d148..3fb3e2a0 100644 --- a/prec/psb_d_bjacprec.f90 +++ b/prec/psb_d_bjacprec.f90 @@ -33,6 +33,7 @@ module psb_d_bjacprec use psb_d_base_prec_mod use psb_d_ilu_fact_mod + use psb_d_ainv_fact_mod type, extends(psb_d_base_prec_type) :: psb_d_bjac_prec_type integer(psb_ipk_), allocatable :: iprcparm(:) @@ -58,8 +59,8 @@ module psb_d_bjacprec end type psb_d_bjac_prec_type character(len=15), parameter, private :: & - & fact_names(0:3)=(/'None ','ILU(0) ',& - & 'ILU(n) ','ILU(eps) '/) + & fact_names(0:4)=(/'None ','ILU(0) ',& + & 'ILU(n) ','ILU(eps) ','AINV(eps) '/) private :: psb_d_bjac_sizeof, psb_d_bjac_precdescr, psb_d_bjac_get_nzeros diff --git a/prec/psb_prec_const_mod.f90 b/prec/psb_prec_const_mod.f90 index 0e0e019b..b21ae9da 100644 --- a/prec/psb_prec_const_mod.f90 +++ b/prec/psb_prec_const_mod.f90 @@ -54,9 +54,11 @@ module psb_prec_const_mod ! Entries in rprcparm: ILU(E) epsilon, smoother omega integer(psb_ipk_), parameter :: psb_ilu_scale_=7 integer(psb_ipk_), parameter :: psb_fact_eps_=1 - integer(psb_ipk_), parameter :: psb_rfpsz=4 + integer(psb_ipk_), parameter :: psb_rfpsz=8 ! Factorization types: none, ILU(0), ILU(N), ILU(N,E) integer(psb_ipk_), parameter :: psb_f_none_=0,psb_f_ilu_n_=1,psb_f_ilu_k_=2,psb_f_ilu_t_=3 + ! Approximate Inverse factorization type: AINV + integer(psb_ipk_), parameter :: psb_f_ainv_=4 ! Fields for sparse matrices ensembles: integer(psb_ipk_), parameter :: psb_l_pr_=1, psb_u_pr_=2, psb_bp_ilu_avsz=2 integer(psb_ipk_), parameter :: psb_max_avsz=psb_bp_ilu_avsz @@ -71,6 +73,20 @@ module psb_prec_const_mod integer(psb_ipk_), parameter :: psb_ilu_scale_aclsum_ = 4 integer(psb_ipk_), parameter :: psb_ilu_scale_arcsum_ = 5 + ! Numerical parameters relative to Approximate Inverse Preconditioners + integer, parameter :: psb_inv_fillin_ = 3 + integer, parameter :: psb_ainv_alg_ = psb_inv_fillin_ + 1 + integer, parameter :: psb_inv_thresh_ = 3 + integer, parameter :: psb_ainv_llk_ = psb_inv_thresh_ + 1 + integer, parameter :: psb_ainv_s_llk_ = psb_ainv_llk_ + 1 + integer, parameter :: psb_ainv_s_ft_llk_ = psb_ainv_s_llk_ + 1 + integer, parameter :: psb_ainv_llk_noth_ = psb_ainv_s_ft_llk_ + 1 + integer, parameter :: psb_ainv_mlk_ = psb_ainv_llk_noth_ + 1 + integer, parameter :: psb_ainv_lmx_ = psb_ainv_mlk_ +#if defined(HAVE_TUMA_SAINV) + integer, parameter :: psb_ainv_s_tuma_ = psb_ainv_lmx_ + 1 + integer, parameter :: psb_ainv_l_tuma_ = psb_ainv_s_tuma_ + 1 +#endif interface psb_check_def diff --git a/prec/psb_s_ainv_fact_mod.f90 b/prec/psb_s_ainv_fact_mod.f90 new file mode 100644 index 00000000..bc7f1d12 --- /dev/null +++ b/prec/psb_s_ainv_fact_mod.f90 @@ -0,0 +1,98 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! 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 PSBLAS 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 PSBLAS 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. +! +! Moved here from MLD2P4, original copyright below. +! +! +! MLD2P4 version 2.2 +! MultiLevel Domain Decomposition Parallel Preconditioners Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2008-2018 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! +! 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: psb_s_ainv_fact_mod.f90 +! +! Module: psb_s_ainv_fact_mod +! +! This module defines some interfaces used internally by the implementation of +! psb_s_ainv_solver, but not visible to the end user. +! +! +module psb_s_ainv_fact_mod + use psb_base_mod + use psb_prec_const_mod + + interface psb_ainv_fact + subroutine psb_s_ainv_bld(a,alg,fillin,thresh,wmat,d,zmat,desc,info,blck,iscale) + import psb_sspmat_type, psb_spk_, psb_ipk_, psb_desc_type + type(psb_sspmat_type), intent(in), target :: a + integer(psb_ipk_), intent(in) :: fillin,alg + real(psb_spk_), intent(in) :: thresh + type(psb_sspmat_type), intent(inout) :: wmat, zmat + real(psb_spk_), allocatable :: d(:) + Type(psb_desc_type), Intent(in) :: desc + integer(psb_ipk_), intent(out) :: info + type(psb_sspmat_type), intent(in), optional :: blck + integer(psb_ipk_), intent(in), optional :: iscale + end subroutine psb_s_ainv_bld + end interface + +end module psb_s_ainv_fact_mod diff --git a/prec/psb_s_biconjg_mod.F90 b/prec/psb_s_biconjg_mod.F90 index ec5f14b8..54022344 100644 --- a/prec/psb_s_biconjg_mod.F90 +++ b/prec/psb_s_biconjg_mod.F90 @@ -100,11 +100,13 @@ ! module psb_s_biconjg_mod + use psb_base_mod + use psb_prec_const_mod + interface psb_sparse_biconjg module procedure psb_ssparse_biconjg end interface - abstract interface subroutine psb_ssparse_biconjg_variant(n,a,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod, only : psb_s_csr_sparse_mat, psb_s_csc_sparse_mat, & @@ -131,12 +133,11 @@ module psb_s_biconjg_mod & psb_ssparse_tuma_lainv #endif - contains subroutine psb_ssparse_biconjg(alg,n,acsr,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod - use psb_base_ainv_mod + use psb_prec_const_mod integer(psb_ipk_), intent(in) :: alg,n type(psb_s_csr_sparse_mat), intent(in) :: acsr type(psb_sspmat_type), intent(out) :: z, w diff --git a/prec/psb_s_bjacprec.f90 b/prec/psb_s_bjacprec.f90 index 2cdd184e..6c7e9c9e 100644 --- a/prec/psb_s_bjacprec.f90 +++ b/prec/psb_s_bjacprec.f90 @@ -33,6 +33,7 @@ module psb_s_bjacprec use psb_s_base_prec_mod use psb_s_ilu_fact_mod + use psb_s_ainv_fact_mod type, extends(psb_s_base_prec_type) :: psb_s_bjac_prec_type integer(psb_ipk_), allocatable :: iprcparm(:) @@ -58,8 +59,8 @@ module psb_s_bjacprec end type psb_s_bjac_prec_type character(len=15), parameter, private :: & - & fact_names(0:3)=(/'None ','ILU(0) ',& - & 'ILU(n) ','ILU(eps) '/) + & fact_names(0:4)=(/'None ','ILU(0) ',& + & 'ILU(n) ','ILU(eps) ','AINV(eps) '/) private :: psb_s_bjac_sizeof, psb_s_bjac_precdescr, psb_s_bjac_get_nzeros diff --git a/prec/psb_z_ainv_fact_mod.f90 b/prec/psb_z_ainv_fact_mod.f90 new file mode 100644 index 00000000..490fe132 --- /dev/null +++ b/prec/psb_z_ainv_fact_mod.f90 @@ -0,0 +1,98 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! 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 PSBLAS 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 PSBLAS 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. +! +! Moved here from MLD2P4, original copyright below. +! +! +! MLD2P4 version 2.2 +! MultiLevel Domain Decomposition Parallel Preconditioners Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2008-2018 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! +! 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: psb_z_ainv_fact_mod.f90 +! +! Module: psb_z_ainv_fact_mod +! +! This module defines some interfaces used internally by the implementation of +! psb_z_ainv_solver, but not visible to the end user. +! +! +module psb_z_ainv_fact_mod + use psb_base_mod + use psb_prec_const_mod + + interface psb_ainv_fact + subroutine psb_z_ainv_bld(a,alg,fillin,thresh,wmat,d,zmat,desc,info,blck,iscale) + import psb_zspmat_type, psb_dpk_, psb_ipk_, psb_desc_type + type(psb_zspmat_type), intent(in), target :: a + integer(psb_ipk_), intent(in) :: fillin,alg + real(psb_dpk_), intent(in) :: thresh + type(psb_zspmat_type), intent(inout) :: wmat, zmat + complex(psb_dpk_), allocatable :: d(:) + Type(psb_desc_type), Intent(in) :: desc + integer(psb_ipk_), intent(out) :: info + type(psb_zspmat_type), intent(in), optional :: blck + integer(psb_ipk_), intent(in), optional :: iscale + end subroutine psb_z_ainv_bld + end interface + +end module psb_z_ainv_fact_mod diff --git a/prec/psb_z_biconjg_mod.F90 b/prec/psb_z_biconjg_mod.F90 index ccaa3aa8..4467a2f3 100644 --- a/prec/psb_z_biconjg_mod.F90 +++ b/prec/psb_z_biconjg_mod.F90 @@ -100,11 +100,13 @@ ! module psb_z_biconjg_mod + use psb_base_mod + use psb_prec_const_mod + interface psb_sparse_biconjg module procedure psb_zsparse_biconjg end interface - abstract interface subroutine psb_zsparse_biconjg_variant(n,a,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod, only : psb_z_csr_sparse_mat, psb_z_csc_sparse_mat, & @@ -131,12 +133,11 @@ module psb_z_biconjg_mod & psb_zsparse_tuma_lainv #endif - contains subroutine psb_zsparse_biconjg(alg,n,acsr,p,z,w,nzrmax,sp_thresh,info) use psb_base_mod - use psb_base_ainv_mod + use psb_prec_const_mod integer(psb_ipk_), intent(in) :: alg,n type(psb_z_csr_sparse_mat), intent(in) :: acsr type(psb_zspmat_type), intent(out) :: z, w diff --git a/prec/psb_z_bjacprec.f90 b/prec/psb_z_bjacprec.f90 index 13bb4bcd..23e826b1 100644 --- a/prec/psb_z_bjacprec.f90 +++ b/prec/psb_z_bjacprec.f90 @@ -33,6 +33,7 @@ module psb_z_bjacprec use psb_z_base_prec_mod use psb_z_ilu_fact_mod + use psb_z_ainv_fact_mod type, extends(psb_z_base_prec_type) :: psb_z_bjac_prec_type integer(psb_ipk_), allocatable :: iprcparm(:) @@ -58,8 +59,8 @@ module psb_z_bjacprec end type psb_z_bjac_prec_type character(len=15), parameter, private :: & - & fact_names(0:3)=(/'None ','ILU(0) ',& - & 'ILU(n) ','ILU(eps) '/) + & fact_names(0:4)=(/'None ','ILU(0) ',& + & 'ILU(n) ','ILU(eps) ','AINV(eps) '/) private :: psb_z_bjac_sizeof, psb_z_bjac_precdescr, psb_z_bjac_get_nzeros diff --git a/test/pargen/psb_d_pde3d.f90 b/test/pargen/psb_d_pde3d.f90 index 4bac87a4..c79f23b9 100644 --- a/test/pargen/psb_d_pde3d.f90 +++ b/test/pargen/psb_d_pde3d.f90 @@ -671,15 +671,20 @@ program psb_d_pde3d call prec%set('sub_solve', parms%alg, info) select case (psb_toupper(parms%alg)) case ("ILU") - call prec%set('sub_fillin', parms%fill, info) - call prec%set('ilu_alg', parms%ilu_alg, info) + call prec%set('sub_fillin', parms%fill, info) + call prec%set('ilu_alg', parms%ilu_alg, info) case ("ILUT") - call prec%set('sub_fillin', parms%fill, info) - call prec%set('sub_iluthrs', parms%thresh, info) + call prec%set('sub_fillin', parms%fill, info) + call prec%set('sub_iluthrs', parms%thresh, info) call prec%set('ilut_scale', parms%ilut_scale, info) + case ("AINV") + call prec%set('inv_thresh', parms%inv_thresh, info) case default ! Do nothing, use default setting in the init routine end select + select case (psb_toupper(parms%orth_alg)) + + end select else ! nothing to set for NONE or DIAG preconditioner end if @@ -881,8 +886,9 @@ contains write(psb_out_unit,'("Threshold : ",es12.5)') parms%thresh write(psb_out_unit,'("Invese Fill in : ",i0)') parms%inv_fill write(psb_out_unit,'("Inverse Threshold : ",es12.5)') parms%inv_thresh - case ('AINVT','AORTH') + case ('AINV','AORTH') write(psb_out_unit,'("Inverse Threshold : ",es12.5)') parms%inv_thresh + write(psb_out_unit,'("Orthogonalization : ",a)') parms%orth_alg case default write(psb_out_unit,'("Unknown diagonal solver")') end select diff --git a/test/pargen/psb_s_pde3d.f90 b/test/pargen/psb_s_pde3d.f90 index 00f9037d..ff1721ab 100644 --- a/test/pargen/psb_s_pde3d.f90 +++ b/test/pargen/psb_s_pde3d.f90 @@ -671,15 +671,20 @@ program psb_s_pde3d call prec%set('sub_solve', parms%alg, info) select case (psb_toupper(parms%alg)) case ("ILU") - call prec%set('sub_fillin', parms%fill, info) - call prec%set('ilu_alg', parms%ilu_alg, info) + call prec%set('sub_fillin', parms%fill, info) + call prec%set('ilu_alg', parms%ilu_alg, info) case ("ILUT") - call prec%set('sub_fillin', parms%fill, info) - call prec%set('sub_iluthrs', parms%thresh, info) + call prec%set('sub_fillin', parms%fill, info) + call prec%set('sub_iluthrs', parms%thresh, info) call prec%set('ilut_scale', parms%ilut_scale, info) + case ("AINV") + call prec%set('inv_thresh', parms%inv_thresh, info) case default ! Do nothing, use default setting in the init routine end select + select case (psb_toupper(parms%orth_alg)) + + end select else ! nothing to set for NONE or DIAG preconditioner end if @@ -881,8 +886,9 @@ contains write(psb_out_unit,'("Threshold : ",es12.5)') parms%thresh write(psb_out_unit,'("Invese Fill in : ",i0)') parms%inv_fill write(psb_out_unit,'("Inverse Threshold : ",es12.5)') parms%inv_thresh - case ('AINVT','AORTH') + case ('AINV','AORTH') write(psb_out_unit,'("Inverse Threshold : ",es12.5)') parms%inv_thresh + write(psb_out_unit,'("Orthogonalization : ",a)') parms%orth_alg case default write(psb_out_unit,'("Unknown diagonal solver")') end select