diff --git a/prec/Makefile b/prec/Makefile index e3b727b7..495f5603 100644 --- a/prec/Makefile +++ b/prec/Makefile @@ -12,12 +12,13 @@ MODOBJS=psb_prec_const_mod.o\ psb_d_diagprec.o psb_d_nullprec.o psb_d_bjacprec.o psb_s_ilu_fact_mod.o \ psb_s_diagprec.o psb_s_nullprec.o psb_s_bjacprec.o psb_d_ilu_fact_mod.o \ psb_c_diagprec.o psb_c_nullprec.o psb_c_bjacprec.o psb_c_ilu_fact_mod.o \ - psb_z_diagprec.o psb_z_nullprec.o psb_z_bjacprec.o psb_z_ilu_fact_mod.o - + psb_z_diagprec.o psb_z_nullprec.o psb_z_bjacprec.o psb_z_ilu_fact_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_ainv_tools_mod.o LIBNAME=$(PRECLIBNAME) COBJS= -FINCLUDES=$(FMFLAG). $(FMFLAG)$(MODDIR) +FINCLUDES=$(FMFLAG). $(FMFLAG)$(MODDIR) OBJS=$(F90OBJS) $(COBJS) $(MPFOBJS) $(MODOBJS) lib: $(OBJS) impld @@ -43,15 +44,16 @@ psb_d_prec_mod.o: psb_prec_type.o psb_c_prec_mod.o: psb_prec_type.o psb_z_prec_mod.o: psb_prec_type.o psb_prec_type.o: psb_s_prec_type.o psb_d_prec_type.o psb_c_prec_type.o psb_z_prec_type.o -psb_prec_mod.o: psb_s_prec_mod.o psb_d_prec_mod.o psb_c_prec_mod.o psb_z_prec_mod.o -psb_s_bjacprec.o psb_s_diagprec.o psb_s_nullprec.o: psb_prec_mod.o psb_s_base_prec_mod.o -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_prec_mod.o: psb_s_prec_mod.o psb_d_prec_mod.o psb_c_prec_mod.o psb_z_prec_mod.o +psb_s_bjacprec.o psb_s_diagprec.o psb_s_nullprec.o: psb_prec_mod.o psb_s_base_prec_mod.o +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_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_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 veryclean: clean /bin/rm -f $(LIBNAME) *$(.mod) @@ -60,4 +62,3 @@ iclean: cd impl && $(MAKE) clean clean: iclean /bin/rm -f $(OBJS) $(LOCAL_MODS) - diff --git a/prec/impl/Makefile b/prec/impl/Makefile index 80e87a54..d2074701 100644 --- a/prec/impl/Makefile +++ b/prec/impl/Makefile @@ -17,13 +17,16 @@ OBJS=psb_s_prec_type_impl.o psb_d_prec_type_impl.o \ psb_cprecbld.o psb_cprecset.o psb_cprecinit.o \ psb_z_diagprec_impl.o psb_z_bjacprec_impl.o psb_z_nullprec_impl.o \ psb_zilu_fct.o psb_z_ilu0_fact.o psb_z_iluk_fact.o psb_z_ilut_fact.o \ - psb_zprecbld.o psb_zprecset.o psb_zprecinit.o + psb_zprecbld.o psb_zprecset.o psb_zprecinit.o \ + psb_c_sparsify.o psb_d_sparsify.o psb_s_sparsify.o psb_z_sparsify.o \ + psb_crwclip.o psb_drwclip.o psb_srwclip.o psb_zrwclip.o \ + psb_c_sp_drop.o psb_d_sp_drop.o psb_s_sp_drop.o psb_z_sp_drop.o LIBNAME=$(PRECLIBNAME) COBJS= -FINCLUDES=$(FMFLAG).. $(FMFLAG)$(MODDIR) +FINCLUDES=$(FMFLAG).. $(FMFLAG)$(MODDIR) -lib: $(OBJS) +lib: $(OBJS) $(AR) $(HERE)/$(LIBNAME) $(OBJS) $(RANLIB) $(HERE)/$(LIBNAME) @@ -31,4 +34,3 @@ veryclean: clean clean: /bin/rm -f $(OBJS) $(LOCAL_MODS) - diff --git a/prec/impl/psb_c_sp_drop.f90 b/prec/impl/psb_c_sp_drop.f90 new file mode 100644 index 00000000..fda59cc7 --- /dev/null +++ b/prec/impl/psb_c_sp_drop.f90 @@ -0,0 +1,168 @@ +! +! 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_sp_drop(idiag,nzrmax,sp_thresh,nz,iz,valz,info) + use psb_base_mod + implicit none + real(psb_spk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(in) :: idiag, nzrmax + integer(psb_ipk_), intent(inout) :: nz + integer(psb_ipk_), intent(inout) :: iz(:) + complex(psb_spk_), intent(inout) :: valz(:) + integer(psb_ipk_), intent(out) :: info + ! + integer(psb_ipk_) :: i, j, idf, nw + complex(psb_spk_) :: witem + integer(psb_ipk_) :: widx + complex(psb_spk_), allocatable :: xw(:) + integer(psb_ipk_), allocatable :: xwid(:), indx(:) + + + info = psb_success_ + + if (nz > min(size(iz),size(valz))) then + write(0,*) 'Serious size problem ',nz,size(iz),size(valz) + info = -2 + return + end if + allocate(xw(nz),xwid(nz),indx(nz),stat=info) + if (info /= psb_success_) then + write(psb_err_unit,*) ' Memory allocation failure in sp_drop',nz,info + return + endif + + ! Always keep the diagonal element + idf = -1 + do i=1, nz + if (iz(i) == idiag) then + idf = i + witem = valz(i) + widx = iz(i) + valz(i) = valz(1) + iz(i) = iz(1) + valz(1) = witem + iz(1) = widx + exit + end if + end do + + if (idf == -1) then + + xw(1:nz) = valz(1:nz) + call psb_qsort(xw(1:nz),indx(1:nz),dir=psb_asort_down_) + do i=1, nz + xwid(i) = iz(indx(i)) + end do + nw = min(nz,nzrmax) + do + if (nw <= 1) exit + if (abs(xw(nw)) < sp_thresh) then + nw = nw - 1 + else + exit + end if + end do + nw = max(nw, 1) + + else + + nw = nz-1 + + xw(1:nw) = valz(2:nz) + + call psb_qsort(xw(1:nw),indx(1:nw),dir=psb_asort_down_) + nw = min(nw,nzrmax-1) + do + if (nw <= 1) exit + if (abs(xw(nw)) < sp_thresh) then + nw = nw - 1 + else + exit + end if + end do + + do i=1, nw + xwid(i) = iz(1+indx(i)) + end do + nw = nw + 1 + xw(nw) = valz(1) + xwid(nw) = iz(1) + end if + + call psb_msort(xwid(1:nw),indx(1:nw),dir=psb_sort_up_) + + do i=1, nw + valz(i) = xw(indx(i)) + iz(i) = xwid(i) + end do + nz = nw + if (nz>nzrmax) write(0,*) 'in sp_drop: ',nw,nzrmax,nz + deallocate(xw,xwid,indx,stat=info) + if (info /= psb_success_) then + write(psb_err_unit,*) ' Memory deallocation failure in sp_drop',info + return + endif + return +end subroutine psb_c_sp_drop diff --git a/prec/impl/psb_c_sparsify.f90 b/prec/impl/psb_c_sparsify.f90 new file mode 100644 index 00000000..b89ce3d0 --- /dev/null +++ b/prec/impl/psb_c_sparsify.f90 @@ -0,0 +1,261 @@ +! +! 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 amg_c_sparsify(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,info,istart,iheap,ikr) + use psb_base_mod + implicit none + + real(psb_spk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(in) :: idiag, n, nzrmax + complex(psb_spk_), intent(inout) :: zw(:) + integer(psb_ipk_), intent(out) :: nz + integer(psb_ipk_), intent(out) :: iz(:) + complex(psb_spk_), intent(out) :: valz(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: istart + type(psb_i_heap), optional :: iheap + integer(psb_ipk_), optional :: ikr(:) + ! + integer(psb_ipk_) :: i, istart_, last_i, iret,k + complex(psb_spk_) :: witem + integer(psb_ipk_) :: widx + complex(psb_spk_), allocatable :: xw(:) + integer(psb_ipk_), allocatable :: xwid(:), indx(:) + type(psb_c_idx_heap) :: heap + + + info = psb_success_ + istart_ = 1 + if (present(istart)) istart_ = max(1,istart) + if (.false.) then + nz = 0 + do i=istart_, n + if ((i == idiag).or.(abs(zw(i)) >= sp_thresh)) then + nz = nz + 1 + iz(nz) = i + valz(nz) = zw(i) + end if + end do + + else + + allocate(xw(nzrmax),xwid(nzrmax),indx(nzrmax),stat=info) + if (info /= psb_success_) then + return + end if + + call heap%init(info,dir=psb_asort_down_) + + ! Keep at least the diagonal + nz = 0 + + if (present(iheap)) then + if (.not.(present(ikr))) then + write(psb_err_unit,*) 'Error: if IHEAP then also IKR' + info = -1 + return + end if + last_i = -1 + do + call iheap%get_first(i,iret) + if (iret < 0) exit + ! An index may have been put on the heap more than once. + if (i == last_i) cycle + last_i = i + if (i == idiag) then + xw(1) = zw(i) + xwid(1) = i + else if (abs(zw(i)) >= sp_thresh) then + call heap%insert(zw(i),i,info) + end if + zw(i) = dzero + ikr(i) = 0 + end do + + else + + do i=istart_, n + if (i == idiag) then + xw(1) = zw(i) + xwid(1) = i + else if (abs(zw(i)) >= sp_thresh) then + call heap%insert(zw(i),i,info) + end if + zw(i) = dzero + end do + end if + + k = 1 + do + if (k == nzrmax) exit + call heap%get_first(witem,widx,info) + if (info == -1) then + info = psb_success_ + exit + endif + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + return + end if + k = k + 1 + xw(k) = witem + xwid(k) = widx + end do + call heap%free(info) + nz = k + call psb_msort(xwid(1:nz),indx(1:nz),dir=psb_sort_up_) +!!$ write(0,*) 'sparsify output for idiag ',idiag,' :',nz,sp_thresh + do i=1, nz + valz(i) = xw(indx(i)) + iz(i) = xwid(i) +!!$ write(0,*) ' ',iz(i),valz(i) + end do + + end if + + return + +end subroutine amg_c_sparsify + + +subroutine amg_c_sparsify_list(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,lhead,listv,ikr,info) + use psb_base_mod + implicit none + + real(psb_spk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(in) :: idiag, n, nzrmax + complex(psb_spk_), intent(inout) :: zw(:) + integer(psb_ipk_), intent(out) :: nz + integer(psb_ipk_), intent(out) :: iz(:) + complex(psb_spk_), intent(out) :: valz(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(inout) :: lhead, listv(:) + integer(psb_ipk_) :: ikr(:) + ! + integer(psb_ipk_) :: i, istart_, last_i, iret,k,current, next + complex(psb_spk_) :: witem + integer(psb_ipk_) :: widx + complex(psb_spk_), allocatable :: xw(:) + integer(psb_ipk_), allocatable :: xwid(:), indx(:) + + + info = psb_success_ + istart_ = 1 + allocate(xw(n),xwid(n),indx(n),stat=info) + + current = lhead + lhead = -1 + i = 0 + do while (current >0) + i = i + 1 + xw(i) = zw(current) + xwid(i) = current + + if (current == idiag) then + ! Bring the diagona into first position + witem = xw(1) + widx = xwid(1) + xw(1) = xw(i) + xwid(1) = xwid(i) + xw(i) = witem + xwid(i) = widx + end if + + zw(current) = dzero + ikr(current) = 0 + + next = listv(current) + listv(current) = -1 + current = next + end do + nz = i + if (nz > 2) call psb_hsort(xw(2:nz),ix=xwid(2:nz),& + & dir=psb_asort_down_,flag=psb_sort_keep_idx_) +!!$ write(0,*) 'Done first msort ' +!!$ write(0,*) ' after first msort for idiag ',idiag,' :',nz,sp_thresh +!!$ do i=1, nz +!!$ write(0,*) ' ',xwid(i),xw(i) +!!$ end do + + i = 2 + do while (i<=nz) + if (abs(xw(i)) < sp_thresh) exit + i = i + 1 + end do +!!$ write(0,*) 'NZ ',nz, i, nzrmax + nz = max(1,min(i-1,nzrmax)) + call psb_msort(xwid(1:nz),ix=indx(1:nz),dir=psb_sort_up_) +!!$ write(0,*) 'Done second msort ' + +!!$ write(0,*) 'sparsify output for idiag ',idiag,' :',nz,i,sp_thresh + do i=1, nz + valz(i) = xw(indx(i)) + iz(i) = xwid(i) +!!$ write(0,*) ' ',iz(i),valz(i) + end do + + return + +end subroutine amg_c_sparsify_list diff --git a/prec/impl/psb_crwclip.f90 b/prec/impl/psb_crwclip.f90 new file mode 100644 index 00000000..941725d2 --- /dev/null +++ b/prec/impl/psb_crwclip.f90 @@ -0,0 +1,90 @@ +! +! 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 amg_c_rwclip(nz,ia,ja,val,imin,imax,jmin,jmax) + use psb_base_mod + + implicit none + integer(psb_ipk_), intent(inout) :: nz + integer(psb_ipk_), intent(inout) :: ia(*), ja(*) + complex(psb_spk_), intent(inout) :: val(*) + integer(psb_ipk_), intent(in) :: imin,imax,jmin,jmax + + integer(psb_ipk_) :: i,j + + j = 0 + do i=1, nz + if ((imin <= ia(i)).and.& + & (ia(i) <= imax).and.& + & (jmin <= ja(i)).and.& + & (ja(i) <= jmax) ) then + j = j + 1 + ia(j) = ia(i) + ja(j) = ja(i) + val(j) = val(i) + end if + end do + nz = j +end subroutine amg_c_rwclip diff --git a/prec/impl/psb_d_sp_drop.f90 b/prec/impl/psb_d_sp_drop.f90 new file mode 100644 index 00000000..67c49b6f --- /dev/null +++ b/prec/impl/psb_d_sp_drop.f90 @@ -0,0 +1,168 @@ +! +! 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_sp_drop(idiag,nzrmax,sp_thresh,nz,iz,valz,info) + use psb_base_mod + implicit none + real(psb_dpk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(in) :: idiag, nzrmax + integer(psb_ipk_), intent(inout) :: nz + integer(psb_ipk_), intent(inout) :: iz(:) + real(psb_dpk_), intent(inout) :: valz(:) + integer(psb_ipk_), intent(out) :: info + ! + integer(psb_ipk_) :: i, j, idf, nw + real(psb_dpk_) :: witem + integer(psb_ipk_) :: widx + real(psb_dpk_), allocatable :: xw(:) + integer(psb_ipk_), allocatable :: xwid(:), indx(:) + + + info = psb_success_ + + if (nz > min(size(iz),size(valz))) then + write(0,*) 'Serious size problem ',nz,size(iz),size(valz) + info = -2 + return + end if + allocate(xw(nz),xwid(nz),indx(nz),stat=info) + if (info /= psb_success_) then + write(psb_err_unit,*) ' Memory allocation failure in sp_drop',nz,info + return + endif + + ! Always keep the diagonal element + idf = -1 + do i=1, nz + if (iz(i) == idiag) then + idf = i + witem = valz(i) + widx = iz(i) + valz(i) = valz(1) + iz(i) = iz(1) + valz(1) = witem + iz(1) = widx + exit + end if + end do + + if (idf == -1) then + + xw(1:nz) = valz(1:nz) + call psb_qsort(xw(1:nz),indx(1:nz),dir=psb_asort_down_) + do i=1, nz + xwid(i) = iz(indx(i)) + end do + nw = min(nz,nzrmax) + do + if (nw <= 1) exit + if (abs(xw(nw)) < sp_thresh) then + nw = nw - 1 + else + exit + end if + end do + nw = max(nw, 1) + + else + + nw = nz-1 + + xw(1:nw) = valz(2:nz) + + call psb_qsort(xw(1:nw),indx(1:nw),dir=psb_asort_down_) + nw = min(nw,nzrmax-1) + do + if (nw <= 1) exit + if (abs(xw(nw)) < sp_thresh) then + nw = nw - 1 + else + exit + end if + end do + + do i=1, nw + xwid(i) = iz(1+indx(i)) + end do + nw = nw + 1 + xw(nw) = valz(1) + xwid(nw) = iz(1) + end if + + call psb_msort(xwid(1:nw),indx(1:nw),dir=psb_sort_up_) + + do i=1, nw + valz(i) = xw(indx(i)) + iz(i) = xwid(i) + end do + nz = nw + if (nz>nzrmax) write(0,*) 'in sp_drop: ',nw,nzrmax,nz + deallocate(xw,xwid,indx,stat=info) + if (info /= psb_success_) then + write(psb_err_unit,*) ' Memory deallocation failure in sp_drop',info + return + endif + return +end subroutine psb_d_sp_drop diff --git a/prec/impl/psb_d_sparsify.f90 b/prec/impl/psb_d_sparsify.f90 new file mode 100644 index 00000000..264f9157 --- /dev/null +++ b/prec/impl/psb_d_sparsify.f90 @@ -0,0 +1,261 @@ +! +! 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 amg_d_sparsify(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,info,istart,iheap,ikr) + use psb_base_mod + implicit none + + real(psb_dpk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(in) :: idiag, n, nzrmax + real(psb_dpk_), intent(inout) :: zw(:) + integer(psb_ipk_), intent(out) :: nz + integer(psb_ipk_), intent(out) :: iz(:) + real(psb_dpk_), intent(out) :: valz(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: istart + type(psb_i_heap), optional :: iheap + integer(psb_ipk_), optional :: ikr(:) + ! + integer(psb_ipk_) :: i, istart_, last_i, iret,k + real(psb_dpk_) :: witem + integer(psb_ipk_) :: widx + real(psb_dpk_), allocatable :: xw(:) + integer(psb_ipk_), allocatable :: xwid(:), indx(:) + type(psb_d_idx_heap) :: heap + + + info = psb_success_ + istart_ = 1 + if (present(istart)) istart_ = max(1,istart) + if (.false.) then + nz = 0 + do i=istart_, n + if ((i == idiag).or.(abs(zw(i)) >= sp_thresh)) then + nz = nz + 1 + iz(nz) = i + valz(nz) = zw(i) + end if + end do + + else + + allocate(xw(nzrmax),xwid(nzrmax),indx(nzrmax),stat=info) + if (info /= psb_success_) then + return + end if + + call heap%init(info,dir=psb_asort_down_) + + ! Keep at least the diagonal + nz = 0 + + if (present(iheap)) then + if (.not.(present(ikr))) then + write(psb_err_unit,*) 'Error: if IHEAP then also IKR' + info = -1 + return + end if + last_i = -1 + do + call iheap%get_first(i,iret) + if (iret < 0) exit + ! An index may have been put on the heap more than once. + if (i == last_i) cycle + last_i = i + if (i == idiag) then + xw(1) = zw(i) + xwid(1) = i + else if (abs(zw(i)) >= sp_thresh) then + call heap%insert(zw(i),i,info) + end if + zw(i) = dzero + ikr(i) = 0 + end do + + else + + do i=istart_, n + if (i == idiag) then + xw(1) = zw(i) + xwid(1) = i + else if (abs(zw(i)) >= sp_thresh) then + call heap%insert(zw(i),i,info) + end if + zw(i) = dzero + end do + end if + + k = 1 + do + if (k == nzrmax) exit + call heap%get_first(witem,widx,info) + if (info == -1) then + info = psb_success_ + exit + endif + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + return + end if + k = k + 1 + xw(k) = witem + xwid(k) = widx + end do + call heap%free(info) + nz = k + call psb_msort(xwid(1:nz),indx(1:nz),dir=psb_sort_up_) +!!$ write(0,*) 'sparsify output for idiag ',idiag,' :',nz,sp_thresh + do i=1, nz + valz(i) = xw(indx(i)) + iz(i) = xwid(i) +!!$ write(0,*) ' ',iz(i),valz(i) + end do + + end if + + return + +end subroutine amg_d_sparsify + + +subroutine amg_d_sparsify_list(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,lhead,listv,ikr,info) + use psb_base_mod + implicit none + + real(psb_dpk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(in) :: idiag, n, nzrmax + real(psb_dpk_), intent(inout) :: zw(:) + integer(psb_ipk_), intent(out) :: nz + integer(psb_ipk_), intent(out) :: iz(:) + real(psb_dpk_), intent(out) :: valz(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(inout) :: lhead, listv(:) + integer(psb_ipk_) :: ikr(:) + ! + integer(psb_ipk_) :: i, istart_, last_i, iret,k,current, next + real(psb_dpk_) :: witem + integer(psb_ipk_) :: widx + real(psb_dpk_), allocatable :: xw(:) + integer(psb_ipk_), allocatable :: xwid(:), indx(:) + + + info = psb_success_ + istart_ = 1 + allocate(xw(n),xwid(n),indx(n),stat=info) + + current = lhead + lhead = -1 + i = 0 + do while (current >0) + i = i + 1 + xw(i) = zw(current) + xwid(i) = current + + if (current == idiag) then + ! Bring the diagona into first position + witem = xw(1) + widx = xwid(1) + xw(1) = xw(i) + xwid(1) = xwid(i) + xw(i) = witem + xwid(i) = widx + end if + + zw(current) = dzero + ikr(current) = 0 + + next = listv(current) + listv(current) = -1 + current = next + end do + nz = i + if (nz > 2) call psb_hsort(xw(2:nz),ix=xwid(2:nz),& + & dir=psb_asort_down_,flag=psb_sort_keep_idx_) +!!$ write(0,*) 'Done first msort ' +!!$ write(0,*) ' after first msort for idiag ',idiag,' :',nz,sp_thresh +!!$ do i=1, nz +!!$ write(0,*) ' ',xwid(i),xw(i) +!!$ end do + + i = 2 + do while (i<=nz) + if (abs(xw(i)) < sp_thresh) exit + i = i + 1 + end do +!!$ write(0,*) 'NZ ',nz, i, nzrmax + nz = max(1,min(i-1,nzrmax)) + call psb_msort(xwid(1:nz),ix=indx(1:nz),dir=psb_sort_up_) +!!$ write(0,*) 'Done second msort ' + +!!$ write(0,*) 'sparsify output for idiag ',idiag,' :',nz,i,sp_thresh + do i=1, nz + valz(i) = xw(indx(i)) + iz(i) = xwid(i) +!!$ write(0,*) ' ',iz(i),valz(i) + end do + + return + +end subroutine amg_d_sparsify_list diff --git a/prec/impl/psb_drwclip.f90 b/prec/impl/psb_drwclip.f90 new file mode 100644 index 00000000..528bde71 --- /dev/null +++ b/prec/impl/psb_drwclip.f90 @@ -0,0 +1,90 @@ +! +! 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 amg_d_rwclip(nz,ia,ja,val,imin,imax,jmin,jmax) + use psb_base_mod + + implicit none + integer(psb_ipk_), intent(inout) :: nz + integer(psb_ipk_), intent(inout) :: ia(*), ja(*) + real(psb_dpk_), intent(inout) :: val(*) + integer(psb_ipk_), intent(in) :: imin,imax,jmin,jmax + + integer(psb_ipk_) :: i,j + + j = 0 + do i=1, nz + if ((imin <= ia(i)).and.& + & (ia(i) <= imax).and.& + & (jmin <= ja(i)).and.& + & (ja(i) <= jmax) ) then + j = j + 1 + ia(j) = ia(i) + ja(j) = ja(i) + val(j) = val(i) + end if + end do + nz = j +end subroutine amg_d_rwclip diff --git a/prec/impl/psb_s_sp_drop.f90 b/prec/impl/psb_s_sp_drop.f90 new file mode 100644 index 00000000..bc297d08 --- /dev/null +++ b/prec/impl/psb_s_sp_drop.f90 @@ -0,0 +1,168 @@ +! +! 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_sp_drop(idiag,nzrmax,sp_thresh,nz,iz,valz,info) + use psb_base_mod + implicit none + real(psb_spk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(in) :: idiag, nzrmax + integer(psb_ipk_), intent(inout) :: nz + integer(psb_ipk_), intent(inout) :: iz(:) + real(psb_spk_), intent(inout) :: valz(:) + integer(psb_ipk_), intent(out) :: info + ! + integer(psb_ipk_) :: i, j, idf, nw + real(psb_spk_) :: witem + integer(psb_ipk_) :: widx + real(psb_spk_), allocatable :: xw(:) + integer(psb_ipk_), allocatable :: xwid(:), indx(:) + + + info = psb_success_ + + if (nz > min(size(iz),size(valz))) then + write(0,*) 'Serious size problem ',nz,size(iz),size(valz) + info = -2 + return + end if + allocate(xw(nz),xwid(nz),indx(nz),stat=info) + if (info /= psb_success_) then + write(psb_err_unit,*) ' Memory allocation failure in sp_drop',nz,info + return + endif + + ! Always keep the diagonal element + idf = -1 + do i=1, nz + if (iz(i) == idiag) then + idf = i + witem = valz(i) + widx = iz(i) + valz(i) = valz(1) + iz(i) = iz(1) + valz(1) = witem + iz(1) = widx + exit + end if + end do + + if (idf == -1) then + + xw(1:nz) = valz(1:nz) + call psb_qsort(xw(1:nz),indx(1:nz),dir=psb_asort_down_) + do i=1, nz + xwid(i) = iz(indx(i)) + end do + nw = min(nz,nzrmax) + do + if (nw <= 1) exit + if (abs(xw(nw)) < sp_thresh) then + nw = nw - 1 + else + exit + end if + end do + nw = max(nw, 1) + + else + + nw = nz-1 + + xw(1:nw) = valz(2:nz) + + call psb_qsort(xw(1:nw),indx(1:nw),dir=psb_asort_down_) + nw = min(nw,nzrmax-1) + do + if (nw <= 1) exit + if (abs(xw(nw)) < sp_thresh) then + nw = nw - 1 + else + exit + end if + end do + + do i=1, nw + xwid(i) = iz(1+indx(i)) + end do + nw = nw + 1 + xw(nw) = valz(1) + xwid(nw) = iz(1) + end if + + call psb_msort(xwid(1:nw),indx(1:nw),dir=psb_sort_up_) + + do i=1, nw + valz(i) = xw(indx(i)) + iz(i) = xwid(i) + end do + nz = nw + if (nz>nzrmax) write(0,*) 'in sp_drop: ',nw,nzrmax,nz + deallocate(xw,xwid,indx,stat=info) + if (info /= psb_success_) then + write(psb_err_unit,*) ' Memory deallocation failure in sp_drop',info + return + endif + return +end subroutine psb_s_sp_drop diff --git a/prec/impl/psb_s_sparsify.f90 b/prec/impl/psb_s_sparsify.f90 new file mode 100644 index 00000000..191b1da5 --- /dev/null +++ b/prec/impl/psb_s_sparsify.f90 @@ -0,0 +1,261 @@ +! +! 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 amg_s_sparsify(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,info,istart,iheap,ikr) + use psb_base_mod + implicit none + + real(psb_spk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(in) :: idiag, n, nzrmax + real(psb_spk_), intent(inout) :: zw(:) + integer(psb_ipk_), intent(out) :: nz + integer(psb_ipk_), intent(out) :: iz(:) + real(psb_spk_), intent(out) :: valz(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: istart + type(psb_i_heap), optional :: iheap + integer(psb_ipk_), optional :: ikr(:) + ! + integer(psb_ipk_) :: i, istart_, last_i, iret,k + real(psb_spk_) :: witem + integer(psb_ipk_) :: widx + real(psb_spk_), allocatable :: xw(:) + integer(psb_ipk_), allocatable :: xwid(:), indx(:) + type(psb_s_idx_heap) :: heap + + + info = psb_success_ + istart_ = 1 + if (present(istart)) istart_ = max(1,istart) + if (.false.) then + nz = 0 + do i=istart_, n + if ((i == idiag).or.(abs(zw(i)) >= sp_thresh)) then + nz = nz + 1 + iz(nz) = i + valz(nz) = zw(i) + end if + end do + + else + + allocate(xw(nzrmax),xwid(nzrmax),indx(nzrmax),stat=info) + if (info /= psb_success_) then + return + end if + + call heap%init(info,dir=psb_asort_down_) + + ! Keep at least the diagonal + nz = 0 + + if (present(iheap)) then + if (.not.(present(ikr))) then + write(psb_err_unit,*) 'Error: if IHEAP then also IKR' + info = -1 + return + end if + last_i = -1 + do + call iheap%get_first(i,iret) + if (iret < 0) exit + ! An index may have been put on the heap more than once. + if (i == last_i) cycle + last_i = i + if (i == idiag) then + xw(1) = zw(i) + xwid(1) = i + else if (abs(zw(i)) >= sp_thresh) then + call heap%insert(zw(i),i,info) + end if + zw(i) = dzero + ikr(i) = 0 + end do + + else + + do i=istart_, n + if (i == idiag) then + xw(1) = zw(i) + xwid(1) = i + else if (abs(zw(i)) >= sp_thresh) then + call heap%insert(zw(i),i,info) + end if + zw(i) = dzero + end do + end if + + k = 1 + do + if (k == nzrmax) exit + call heap%get_first(witem,widx,info) + if (info == -1) then + info = psb_success_ + exit + endif + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + return + end if + k = k + 1 + xw(k) = witem + xwid(k) = widx + end do + call heap%free(info) + nz = k + call psb_msort(xwid(1:nz),indx(1:nz),dir=psb_sort_up_) +!!$ write(0,*) 'sparsify output for idiag ',idiag,' :',nz,sp_thresh + do i=1, nz + valz(i) = xw(indx(i)) + iz(i) = xwid(i) +!!$ write(0,*) ' ',iz(i),valz(i) + end do + + end if + + return + +end subroutine amg_s_sparsify + + +subroutine amg_s_sparsify_list(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,lhead,listv,ikr,info) + use psb_base_mod + implicit none + + real(psb_spk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(in) :: idiag, n, nzrmax + real(psb_spk_), intent(inout) :: zw(:) + integer(psb_ipk_), intent(out) :: nz + integer(psb_ipk_), intent(out) :: iz(:) + real(psb_spk_), intent(out) :: valz(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(inout) :: lhead, listv(:) + integer(psb_ipk_) :: ikr(:) + ! + integer(psb_ipk_) :: i, istart_, last_i, iret,k,current, next + real(psb_spk_) :: witem + integer(psb_ipk_) :: widx + real(psb_spk_), allocatable :: xw(:) + integer(psb_ipk_), allocatable :: xwid(:), indx(:) + + + info = psb_success_ + istart_ = 1 + allocate(xw(n),xwid(n),indx(n),stat=info) + + current = lhead + lhead = -1 + i = 0 + do while (current >0) + i = i + 1 + xw(i) = zw(current) + xwid(i) = current + + if (current == idiag) then + ! Bring the diagona into first position + witem = xw(1) + widx = xwid(1) + xw(1) = xw(i) + xwid(1) = xwid(i) + xw(i) = witem + xwid(i) = widx + end if + + zw(current) = dzero + ikr(current) = 0 + + next = listv(current) + listv(current) = -1 + current = next + end do + nz = i + if (nz > 2) call psb_hsort(xw(2:nz),ix=xwid(2:nz),& + & dir=psb_asort_down_,flag=psb_sort_keep_idx_) +!!$ write(0,*) 'Done first msort ' +!!$ write(0,*) ' after first msort for idiag ',idiag,' :',nz,sp_thresh +!!$ do i=1, nz +!!$ write(0,*) ' ',xwid(i),xw(i) +!!$ end do + + i = 2 + do while (i<=nz) + if (abs(xw(i)) < sp_thresh) exit + i = i + 1 + end do +!!$ write(0,*) 'NZ ',nz, i, nzrmax + nz = max(1,min(i-1,nzrmax)) + call psb_msort(xwid(1:nz),ix=indx(1:nz),dir=psb_sort_up_) +!!$ write(0,*) 'Done second msort ' + +!!$ write(0,*) 'sparsify output for idiag ',idiag,' :',nz,i,sp_thresh + do i=1, nz + valz(i) = xw(indx(i)) + iz(i) = xwid(i) +!!$ write(0,*) ' ',iz(i),valz(i) + end do + + return + +end subroutine amg_s_sparsify_list diff --git a/prec/impl/psb_srwclip.f90 b/prec/impl/psb_srwclip.f90 new file mode 100644 index 00000000..d9c303dd --- /dev/null +++ b/prec/impl/psb_srwclip.f90 @@ -0,0 +1,90 @@ +! +! 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 amg_s_rwclip(nz,ia,ja,val,imin,imax,jmin,jmax) + use psb_base_mod + + implicit none + integer(psb_ipk_), intent(inout) :: nz + integer(psb_ipk_), intent(inout) :: ia(*), ja(*) + real(psb_spk_), intent(inout) :: val(*) + integer(psb_ipk_), intent(in) :: imin,imax,jmin,jmax + + integer(psb_ipk_) :: i,j + + j = 0 + do i=1, nz + if ((imin <= ia(i)).and.& + & (ia(i) <= imax).and.& + & (jmin <= ja(i)).and.& + & (ja(i) <= jmax) ) then + j = j + 1 + ia(j) = ia(i) + ja(j) = ja(i) + val(j) = val(i) + end if + end do + nz = j +end subroutine amg_s_rwclip diff --git a/prec/impl/psb_z_sp_drop.f90 b/prec/impl/psb_z_sp_drop.f90 new file mode 100644 index 00000000..754c76cc --- /dev/null +++ b/prec/impl/psb_z_sp_drop.f90 @@ -0,0 +1,168 @@ +! +! 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_sp_drop(idiag,nzrmax,sp_thresh,nz,iz,valz,info) + use psb_base_mod + implicit none + real(psb_dpk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(in) :: idiag, nzrmax + integer(psb_ipk_), intent(inout) :: nz + integer(psb_ipk_), intent(inout) :: iz(:) + complex(psb_dpk_), intent(inout) :: valz(:) + integer(psb_ipk_), intent(out) :: info + ! + integer(psb_ipk_) :: i, j, idf, nw + complex(psb_dpk_) :: witem + integer(psb_ipk_) :: widx + complex(psb_dpk_), allocatable :: xw(:) + integer(psb_ipk_), allocatable :: xwid(:), indx(:) + + + info = psb_success_ + + if (nz > min(size(iz),size(valz))) then + write(0,*) 'Serious size problem ',nz,size(iz),size(valz) + info = -2 + return + end if + allocate(xw(nz),xwid(nz),indx(nz),stat=info) + if (info /= psb_success_) then + write(psb_err_unit,*) ' Memory allocation failure in sp_drop',nz,info + return + endif + + ! Always keep the diagonal element + idf = -1 + do i=1, nz + if (iz(i) == idiag) then + idf = i + witem = valz(i) + widx = iz(i) + valz(i) = valz(1) + iz(i) = iz(1) + valz(1) = witem + iz(1) = widx + exit + end if + end do + + if (idf == -1) then + + xw(1:nz) = valz(1:nz) + call psb_qsort(xw(1:nz),indx(1:nz),dir=psb_asort_down_) + do i=1, nz + xwid(i) = iz(indx(i)) + end do + nw = min(nz,nzrmax) + do + if (nw <= 1) exit + if (abs(xw(nw)) < sp_thresh) then + nw = nw - 1 + else + exit + end if + end do + nw = max(nw, 1) + + else + + nw = nz-1 + + xw(1:nw) = valz(2:nz) + + call psb_qsort(xw(1:nw),indx(1:nw),dir=psb_asort_down_) + nw = min(nw,nzrmax-1) + do + if (nw <= 1) exit + if (abs(xw(nw)) < sp_thresh) then + nw = nw - 1 + else + exit + end if + end do + + do i=1, nw + xwid(i) = iz(1+indx(i)) + end do + nw = nw + 1 + xw(nw) = valz(1) + xwid(nw) = iz(1) + end if + + call psb_msort(xwid(1:nw),indx(1:nw),dir=psb_sort_up_) + + do i=1, nw + valz(i) = xw(indx(i)) + iz(i) = xwid(i) + end do + nz = nw + if (nz>nzrmax) write(0,*) 'in sp_drop: ',nw,nzrmax,nz + deallocate(xw,xwid,indx,stat=info) + if (info /= psb_success_) then + write(psb_err_unit,*) ' Memory deallocation failure in sp_drop',info + return + endif + return +end subroutine psb_z_sp_drop diff --git a/prec/impl/psb_z_sparsify.f90 b/prec/impl/psb_z_sparsify.f90 new file mode 100644 index 00000000..e19ef19e --- /dev/null +++ b/prec/impl/psb_z_sparsify.f90 @@ -0,0 +1,261 @@ +! +! 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 amg_z_sparsify(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,info,istart,iheap,ikr) + use psb_base_mod + implicit none + + real(psb_dpk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(in) :: idiag, n, nzrmax + complex(psb_dpk_), intent(inout) :: zw(:) + integer(psb_ipk_), intent(out) :: nz + integer(psb_ipk_), intent(out) :: iz(:) + complex(psb_dpk_), intent(out) :: valz(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: istart + type(psb_i_heap), optional :: iheap + integer(psb_ipk_), optional :: ikr(:) + ! + integer(psb_ipk_) :: i, istart_, last_i, iret,k + complex(psb_dpk_) :: witem + integer(psb_ipk_) :: widx + complex(psb_dpk_), allocatable :: xw(:) + integer(psb_ipk_), allocatable :: xwid(:), indx(:) + type(psb_z_idx_heap) :: heap + + + info = psb_success_ + istart_ = 1 + if (present(istart)) istart_ = max(1,istart) + if (.false.) then + nz = 0 + do i=istart_, n + if ((i == idiag).or.(abs(zw(i)) >= sp_thresh)) then + nz = nz + 1 + iz(nz) = i + valz(nz) = zw(i) + end if + end do + + else + + allocate(xw(nzrmax),xwid(nzrmax),indx(nzrmax),stat=info) + if (info /= psb_success_) then + return + end if + + call heap%init(info,dir=psb_asort_down_) + + ! Keep at least the diagonal + nz = 0 + + if (present(iheap)) then + if (.not.(present(ikr))) then + write(psb_err_unit,*) 'Error: if IHEAP then also IKR' + info = -1 + return + end if + last_i = -1 + do + call iheap%get_first(i,iret) + if (iret < 0) exit + ! An index may have been put on the heap more than once. + if (i == last_i) cycle + last_i = i + if (i == idiag) then + xw(1) = zw(i) + xwid(1) = i + else if (abs(zw(i)) >= sp_thresh) then + call heap%insert(zw(i),i,info) + end if + zw(i) = dzero + ikr(i) = 0 + end do + + else + + do i=istart_, n + if (i == idiag) then + xw(1) = zw(i) + xwid(1) = i + else if (abs(zw(i)) >= sp_thresh) then + call heap%insert(zw(i),i,info) + end if + zw(i) = dzero + end do + end if + + k = 1 + do + if (k == nzrmax) exit + call heap%get_first(witem,widx,info) + if (info == -1) then + info = psb_success_ + exit + endif + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + return + end if + k = k + 1 + xw(k) = witem + xwid(k) = widx + end do + call heap%free(info) + nz = k + call psb_msort(xwid(1:nz),indx(1:nz),dir=psb_sort_up_) +!!$ write(0,*) 'sparsify output for idiag ',idiag,' :',nz,sp_thresh + do i=1, nz + valz(i) = xw(indx(i)) + iz(i) = xwid(i) +!!$ write(0,*) ' ',iz(i),valz(i) + end do + + end if + + return + +end subroutine amg_z_sparsify + + +subroutine amg_z_sparsify_list(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,lhead,listv,ikr,info) + use psb_base_mod + implicit none + + real(psb_dpk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(in) :: idiag, n, nzrmax + complex(psb_dpk_), intent(inout) :: zw(:) + integer(psb_ipk_), intent(out) :: nz + integer(psb_ipk_), intent(out) :: iz(:) + complex(psb_dpk_), intent(out) :: valz(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(inout) :: lhead, listv(:) + integer(psb_ipk_) :: ikr(:) + ! + integer(psb_ipk_) :: i, istart_, last_i, iret,k,current, next + complex(psb_dpk_) :: witem + integer(psb_ipk_) :: widx + complex(psb_dpk_), allocatable :: xw(:) + integer(psb_ipk_), allocatable :: xwid(:), indx(:) + + + info = psb_success_ + istart_ = 1 + allocate(xw(n),xwid(n),indx(n),stat=info) + + current = lhead + lhead = -1 + i = 0 + do while (current >0) + i = i + 1 + xw(i) = zw(current) + xwid(i) = current + + if (current == idiag) then + ! Bring the diagona into first position + witem = xw(1) + widx = xwid(1) + xw(1) = xw(i) + xwid(1) = xwid(i) + xw(i) = witem + xwid(i) = widx + end if + + zw(current) = dzero + ikr(current) = 0 + + next = listv(current) + listv(current) = -1 + current = next + end do + nz = i + if (nz > 2) call psb_hsort(xw(2:nz),ix=xwid(2:nz),& + & dir=psb_asort_down_,flag=psb_sort_keep_idx_) +!!$ write(0,*) 'Done first msort ' +!!$ write(0,*) ' after first msort for idiag ',idiag,' :',nz,sp_thresh +!!$ do i=1, nz +!!$ write(0,*) ' ',xwid(i),xw(i) +!!$ end do + + i = 2 + do while (i<=nz) + if (abs(xw(i)) < sp_thresh) exit + i = i + 1 + end do +!!$ write(0,*) 'NZ ',nz, i, nzrmax + nz = max(1,min(i-1,nzrmax)) + call psb_msort(xwid(1:nz),ix=indx(1:nz),dir=psb_sort_up_) +!!$ write(0,*) 'Done second msort ' + +!!$ write(0,*) 'sparsify output for idiag ',idiag,' :',nz,i,sp_thresh + do i=1, nz + valz(i) = xw(indx(i)) + iz(i) = xwid(i) +!!$ write(0,*) ' ',iz(i),valz(i) + end do + + return + +end subroutine amg_z_sparsify_list diff --git a/prec/impl/psb_zrwclip.f90 b/prec/impl/psb_zrwclip.f90 new file mode 100644 index 00000000..41de9603 --- /dev/null +++ b/prec/impl/psb_zrwclip.f90 @@ -0,0 +1,90 @@ +! +! 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 amg_z_rwclip(nz,ia,ja,val,imin,imax,jmin,jmax) + use psb_base_mod + + implicit none + integer(psb_ipk_), intent(inout) :: nz + integer(psb_ipk_), intent(inout) :: ia(*), ja(*) + complex(psb_dpk_), intent(inout) :: val(*) + integer(psb_ipk_), intent(in) :: imin,imax,jmin,jmax + + integer(psb_ipk_) :: i,j + + j = 0 + do i=1, nz + if ((imin <= ia(i)).and.& + & (ia(i) <= imax).and.& + & (jmin <= ja(i)).and.& + & (ja(i) <= jmax) ) then + j = j + 1 + ia(j) = ia(i) + ja(j) = ja(i) + val(j) = val(i) + end if + end do + nz = j +end subroutine amg_z_rwclip diff --git a/prec/psb_ainv_tools_mod.f90 b/prec/psb_ainv_tools_mod.f90 new file mode 100644 index 00000000..561f0f79 --- /dev/null +++ b/prec/psb_ainv_tools_mod.f90 @@ -0,0 +1,6 @@ +module psb_ainv_tools_mod + use psb_c_ainv_tools_mod + use psb_d_ainv_tools_mod + use psb_s_ainv_tools_mod + use psb_z_ainv_tools_mod +end module psb_ainv_tools_mod diff --git a/prec/psb_base_ainv_mod.F90 b/prec/psb_base_ainv_mod.F90 new file mode 100644 index 00000000..71bd5197 --- /dev/null +++ b/prec/psb_base_ainv_mod.F90 @@ -0,0 +1,94 @@ +! +! 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. +! +! +! +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 + + +end module psb_base_ainv_mod diff --git a/prec/psb_c_ainv_tools_mod.f90 b/prec/psb_c_ainv_tools_mod.f90 new file mode 100644 index 00000000..a68d6d69 --- /dev/null +++ b/prec/psb_c_ainv_tools_mod.f90 @@ -0,0 +1,132 @@ +! +! 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. +! +! +! +! +! +! +! +module psb_c_ainv_tools_mod + + interface sp_drop + subroutine psb_c_sp_drop(idiag,nzrmax,sp_thresh,nz,iz,valz,info) + use psb_base_mod, only : psb_spk_, psb_ipk_ + implicit none + real(psb_spk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(in) :: idiag, nzrmax + integer(psb_ipk_), intent(inout) :: nz + integer(psb_ipk_), intent(inout) :: iz(:) + complex(psb_spk_), intent(inout) :: valz(:) + integer(psb_ipk_), intent(out) :: info + end subroutine psb_c_sp_drop + end interface + + interface rwclip + subroutine psb_c_rwclip(nz,ia,ja,val,imin,imax,jmin,jmax) + use psb_base_mod, only : psb_spk_, psb_ipk_ + + implicit none + integer(psb_ipk_), intent(inout) :: nz + integer(psb_ipk_), intent(inout) :: ia(*), ja(*) + complex(psb_spk_), intent(inout) :: val(*) + integer(psb_ipk_), intent(in) :: imin,imax,jmin,jmax + end subroutine psb_c_rwclip + end interface + + interface sparsify + subroutine psb_c_sparsify(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,info, & + & istart,iheap,ikr) + use psb_base_mod, only : psb_spk_, psb_ipk_, psb_i_heap + implicit none + + real(psb_spk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(in) :: idiag, n, nzrmax + complex(psb_spk_), intent(inout) :: zw(:) + integer(psb_ipk_), intent(out) :: nz + integer(psb_ipk_), intent(out) :: iz(:) + complex(psb_spk_), intent(out) :: valz(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: istart + type(psb_i_heap), optional :: iheap + integer(psb_ipk_), optional :: ikr(:) + end subroutine psb_c_sparsify + subroutine psb_c_sparsify_list(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,lhead,listv,ikr,info) + use psb_base_mod, only : psb_spk_, psb_ipk_ + implicit none + + real(psb_spk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(in) :: idiag, n, nzrmax + complex(psb_spk_), intent(inout) :: zw(:) + integer(psb_ipk_), intent(out) :: nz + integer(psb_ipk_), intent(out) :: iz(:) + complex(psb_spk_), intent(out) :: valz(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(inout) :: lhead, listv(:) + integer(psb_ipk_) :: ikr(:) + end subroutine psb_c_sparsify_list + + end interface + +end module psb_c_ainv_tools_mod diff --git a/prec/psb_d_ainv_tools_mod.f90 b/prec/psb_d_ainv_tools_mod.f90 new file mode 100644 index 00000000..7329533b --- /dev/null +++ b/prec/psb_d_ainv_tools_mod.f90 @@ -0,0 +1,132 @@ +! +! 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. +! +! +! +! +! +! +! +module psb_d_ainv_tools_mod + + interface sp_drop + subroutine psb_d_sp_drop(idiag,nzrmax,sp_thresh,nz,iz,valz,info) + use psb_base_mod, only : psb_dpk_, psb_ipk_ + implicit none + real(psb_dpk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(in) :: idiag, nzrmax + integer(psb_ipk_), intent(inout) :: nz + integer(psb_ipk_), intent(inout) :: iz(:) + real(psb_dpk_), intent(inout) :: valz(:) + integer(psb_ipk_), intent(out) :: info + end subroutine psb_d_sp_drop + end interface + + interface rwclip + subroutine psb_d_rwclip(nz,ia,ja,val,imin,imax,jmin,jmax) + use psb_base_mod, only : psb_dpk_, psb_ipk_ + + implicit none + integer(psb_ipk_), intent(inout) :: nz + integer(psb_ipk_), intent(inout) :: ia(*), ja(*) + real(psb_dpk_), intent(inout) :: val(*) + integer(psb_ipk_), intent(in) :: imin,imax,jmin,jmax + end subroutine psb_d_rwclip + end interface + + interface sparsify + subroutine psb_d_sparsify(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,info, & + & istart,iheap,ikr) + use psb_base_mod, only : psb_dpk_, psb_ipk_, psb_i_heap + implicit none + + real(psb_dpk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(in) :: idiag, n, nzrmax + real(psb_dpk_), intent(inout) :: zw(:) + integer(psb_ipk_), intent(out) :: nz + integer(psb_ipk_), intent(out) :: iz(:) + real(psb_dpk_), intent(out) :: valz(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: istart + type(psb_i_heap), optional :: iheap + integer(psb_ipk_), optional :: ikr(:) + end subroutine psb_d_sparsify + subroutine psb_d_sparsify_list(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,lhead,listv,ikr,info) + use psb_base_mod, only : psb_dpk_, psb_ipk_ + implicit none + + real(psb_dpk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(in) :: idiag, n, nzrmax + real(psb_dpk_), intent(inout) :: zw(:) + integer(psb_ipk_), intent(out) :: nz + integer(psb_ipk_), intent(out) :: iz(:) + real(psb_dpk_), intent(out) :: valz(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(inout) :: lhead, listv(:) + integer(psb_ipk_) :: ikr(:) + end subroutine psb_d_sparsify_list + + end interface + +end module psb_d_ainv_tools_mod diff --git a/prec/psb_s_ainv_tools_mod.f90 b/prec/psb_s_ainv_tools_mod.f90 new file mode 100644 index 00000000..caa50164 --- /dev/null +++ b/prec/psb_s_ainv_tools_mod.f90 @@ -0,0 +1,132 @@ +! +! 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. +! +! +! +! +! +! +! +module psb_s_ainv_tools_mod + + interface sp_drop + subroutine psb_s_sp_drop(idiag,nzrmax,sp_thresh,nz,iz,valz,info) + use psb_base_mod, only : psb_spk_, psb_ipk_ + implicit none + real(psb_spk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(in) :: idiag, nzrmax + integer(psb_ipk_), intent(inout) :: nz + integer(psb_ipk_), intent(inout) :: iz(:) + real(psb_spk_), intent(inout) :: valz(:) + integer(psb_ipk_), intent(out) :: info + end subroutine psb_s_sp_drop + end interface + + interface rwclip + subroutine psb_s_rwclip(nz,ia,ja,val,imin,imax,jmin,jmax) + use psb_base_mod, only : psb_spk_, psb_ipk_ + + implicit none + integer(psb_ipk_), intent(inout) :: nz + integer(psb_ipk_), intent(inout) :: ia(*), ja(*) + real(psb_spk_), intent(inout) :: val(*) + integer(psb_ipk_), intent(in) :: imin,imax,jmin,jmax + end subroutine psb_s_rwclip + end interface + + interface sparsify + subroutine psb_s_sparsify(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,info, & + & istart,iheap,ikr) + use psb_base_mod, only : psb_spk_, psb_ipk_, psb_i_heap + implicit none + + real(psb_spk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(in) :: idiag, n, nzrmax + real(psb_spk_), intent(inout) :: zw(:) + integer(psb_ipk_), intent(out) :: nz + integer(psb_ipk_), intent(out) :: iz(:) + real(psb_spk_), intent(out) :: valz(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: istart + type(psb_i_heap), optional :: iheap + integer(psb_ipk_), optional :: ikr(:) + end subroutine psb_s_sparsify + subroutine psb_s_sparsify_list(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,lhead,listv,ikr,info) + use psb_base_mod, only : psb_spk_, psb_ipk_ + implicit none + + real(psb_spk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(in) :: idiag, n, nzrmax + real(psb_spk_), intent(inout) :: zw(:) + integer(psb_ipk_), intent(out) :: nz + integer(psb_ipk_), intent(out) :: iz(:) + real(psb_spk_), intent(out) :: valz(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(inout) :: lhead, listv(:) + integer(psb_ipk_) :: ikr(:) + end subroutine psb_s_sparsify_list + + end interface + +end module psb_s_ainv_tools_mod diff --git a/prec/psb_z_ainv_tools_mod.f90 b/prec/psb_z_ainv_tools_mod.f90 new file mode 100644 index 00000000..f611c2a7 --- /dev/null +++ b/prec/psb_z_ainv_tools_mod.f90 @@ -0,0 +1,132 @@ +! +! 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. +! +! +! +! +! +! +! +module psb_z_ainv_tools_mod + + interface sp_drop + subroutine psb_z_sp_drop(idiag,nzrmax,sp_thresh,nz,iz,valz,info) + use psb_base_mod, only : psb_dpk_, psb_ipk_ + implicit none + real(psb_dpk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(in) :: idiag, nzrmax + integer(psb_ipk_), intent(inout) :: nz + integer(psb_ipk_), intent(inout) :: iz(:) + complex(psb_dpk_), intent(inout) :: valz(:) + integer(psb_ipk_), intent(out) :: info + end subroutine psb_z_sp_drop + end interface + + interface rwclip + subroutine psb_z_rwclip(nz,ia,ja,val,imin,imax,jmin,jmax) + use psb_base_mod, only : psb_dpk_, psb_ipk_ + + implicit none + integer(psb_ipk_), intent(inout) :: nz + integer(psb_ipk_), intent(inout) :: ia(*), ja(*) + complex(psb_dpk_), intent(inout) :: val(*) + integer(psb_ipk_), intent(in) :: imin,imax,jmin,jmax + end subroutine psb_z_rwclip + end interface + + interface sparsify + subroutine psb_z_sparsify(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,info, & + & istart,iheap,ikr) + use psb_base_mod, only : psb_dpk_, psb_ipk_, psb_i_heap + implicit none + + real(psb_dpk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(in) :: idiag, n, nzrmax + complex(psb_dpk_), intent(inout) :: zw(:) + integer(psb_ipk_), intent(out) :: nz + integer(psb_ipk_), intent(out) :: iz(:) + complex(psb_dpk_), intent(out) :: valz(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: istart + type(psb_i_heap), optional :: iheap + integer(psb_ipk_), optional :: ikr(:) + end subroutine psb_z_sparsify + subroutine psb_z_sparsify_list(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,lhead,listv,ikr,info) + use psb_base_mod, only : psb_dpk_, psb_ipk_ + implicit none + + real(psb_dpk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(in) :: idiag, n, nzrmax + complex(psb_dpk_), intent(inout) :: zw(:) + integer(psb_ipk_), intent(out) :: nz + integer(psb_ipk_), intent(out) :: iz(:) + complex(psb_dpk_), intent(out) :: valz(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(inout) :: lhead, listv(:) + integer(psb_ipk_) :: ikr(:) + end subroutine psb_z_sparsify_list + + end interface + +end module psb_z_ainv_tools_mod