diff --git a/prec/Makefile b/prec/Makefile index 15284a93..5fa10d4d 100644 --- a/prec/Makefile +++ b/prec/Makefile @@ -17,7 +17,11 @@ MODOBJS=psb_prec_const_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 \ psb_biconjg_mod.o psb_c_biconjg_mod.o psb_d_biconjg_mod.o psb_s_biconjg_mod.o \ - psb_z_biconjg_mod.o + psb_z_biconjg_mod.o \ + psb_c_invt_fact_mod.o psb_d_invt_fact_mod.o psb_s_invt_fact_mod.o \ + psb_z_invt_fact_mod.o\ + psb_c_invk_fact_mod.o psb_d_invk_fact_mod.o psb_s_invk_fact_mod.o \ + psb_z_invk_fact_mod.o LIBNAME=$(PRECLIBNAME) COBJS= @@ -63,6 +67,14 @@ psb_z_ainv_fact_mod.o: psb_prec_const_mod.o psb_ainv_tools_mod.o psb_ainv_tools_mod.o: psb_c_ainv_tools_mod.o psb_d_ainv_tools_mod.o psb_s_ainv_tools_mod.o psb_z_ainv_tools_mod.o psb_biconjg_mod.o: psb_prec_const_mod.o psb_c_biconjg_mod.o \ psb_d_biconjg_mod.o psb_s_biconjg_mod.o psb_z_biconjg_mod.o +psb_c_invt_fact_mod.o: psb_prec_const_mod.o psb_c_ilu_fact_mod.o +psb_d_invt_fact_mod.o: psb_prec_const_mod.o psb_d_ilu_fact_mod.o +psb_s_invt_fact_mod.o: psb_prec_const_mod.o psb_s_ilu_fact_mod.o +psb_z_invt_fact_mod.o: psb_prec_const_mod.o psb_z_ilu_fact_mod.o +psb_c_invk_fact_mod.o: psb_prec_const_mod.o psb_c_ilu_fact_mod.o +psb_d_invk_fact_mod.o: psb_prec_const_mod.o psb_d_ilu_fact_mod.o +psb_s_invk_fact_mod.o: psb_prec_const_mod.o psb_s_ilu_fact_mod.o +psb_z_invk_fact_mod.o: psb_prec_const_mod.o psb_z_ilu_fact_mod.o veryclean: clean /bin/rm -f $(LIBNAME) *$(.mod) diff --git a/prec/impl/Makefile b/prec/impl/Makefile index 284eb6e0..57d7c304 100644 --- a/prec/impl/Makefile +++ b/prec/impl/Makefile @@ -34,7 +34,9 @@ OBJS=psb_s_prec_type_impl.o psb_d_prec_type_impl.o \ psb_ssparse_biconjg_mlk.o psb_ssparse_biconjg_s_ft_llk.o \ psb_ssparse_biconjg_s_llk.o \ psb_d_ainv_bld.o psb_c_ainv_bld.o psb_s_ainv_bld.o \ - psb_z_ainv_bld.o + psb_z_ainv_bld.o \ + psb_c_invt_fact.o psb_d_invt_fact.o psb_s_invt_fact.o psb_z_invt_fact.o\ + psb_c_invk_fact.o psb_d_invk_fact.o psb_s_invk_fact.o psb_z_invk_fact.o LIBNAME=$(PRECLIBNAME) COBJS= diff --git a/prec/impl/psb_c_invk_fact.f90 b/prec/impl/psb_c_invk_fact.f90 new file mode 100644 index 00000000..c7c98970 --- /dev/null +++ b/prec/impl/psb_c_invk_fact.f90 @@ -0,0 +1,494 @@ +! +! +! 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_invk_bld(a,fill1, fill2,lmat,d,umat,desc,info,blck) + + use psb_base_mod + use psb_c_invk_fact_mod, psb_protect_name => psb_c_invk_bld + use psb_c_ilu_fact_mod + implicit none + + ! Arguments + type(psb_cspmat_type), intent(in), target :: a + integer(psb_ipk_), intent(in) :: fill1, fill2 + type(psb_cspmat_type), intent(inout) :: lmat, umat + complex(psb_spk_), allocatable :: d(:) + Type(psb_desc_type), Intent(inout) :: desc + integer(psb_ipk_), intent(out) :: info + type(psb_cspmat_type), intent(in), optional :: blck + ! + integer(psb_ipk_) :: i, nztota, err_act, n_row, nrow_a, n_col + type(psb_cspmat_type) :: atmp + complex(psb_spk_), allocatable :: pq(:), pd(:) + integer(psb_ipk_), allocatable :: uplevs(:) + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: nzrmax + character(len=20) :: name, ch_err + + + info = psb_success_ + name='psb_cinvk_bld' + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = psb_cd_get_context(desc) + call psb_info(ictxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + ! + ! Check the memory available to hold the incomplete L and U factors + ! and allocate it if needed + ! + nrow_a = a%get_nrows() + nztota = a%get_nzeros() + + if (present(blck)) then + nztota = nztota + blck%get_nzeros() + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & ': out get_nnzeros',nrow_a,nztota,& + & a%get_nrows(),a%get_ncols(),a%get_nzeros() + + + n_row = psb_cd_get_local_rows(desc) + n_col = psb_cd_get_local_cols(desc) + allocate(pd(n_row),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + + call lmat%allocate(n_row,n_row,info,nz=nztota) + if (info == psb_success_) call umat%allocate(n_row,n_row,info,nz=nztota) + + + call psb_iluk_fact(fill1,psb_ilu_n_,a,lmat,umat,pd,info,blck=blck) + !,uplevs=uplevs) + !call psb_ciluk_fact(fill1,psb_ilu_n_,a,lmat,umat,pd,info,blck=blck) + + if (info == psb_success_) call atmp%allocate(n_row,n_row,info,nz=nztota) + if(info/=0) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + + ! + ! Compute the aprox U^-1 and L^-1 + ! + call psb_sparse_invk(n_row,umat,atmp,fill2,info) + if (info == psb_success_) call psb_move_alloc(atmp,umat,info) + if (info == psb_success_) call lmat%transp() + if (info == psb_success_) call psb_sparse_invk(n_row,lmat,atmp,fill2,info) + if (info == psb_success_) call psb_move_alloc(atmp,lmat,info) + if (info == psb_success_) call lmat%transp() + ! Done. Hopefully.... + + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='invt') + goto 9999 + end if + + call psb_move_alloc(pd,d,info) + call lmat%set_asb() + call lmat%trim() + call umat%set_asb() + call umat%trim() + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_c_invk_bld + +subroutine psb_csparse_invk(n,a,z,fill_in,info,inlevs) + + use psb_base_mod + use psb_c_invk_fact_mod, psb_protect_name => psb_csparse_invk + + integer(psb_ipk_), intent(in) :: n + type(psb_cspmat_type), intent(in) :: a + type(psb_cspmat_type), intent(inout) :: z + integer(psb_ipk_), intent(in) :: fill_in + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: inlevs(:) + ! + integer(psb_ipk_) :: i,j,k, err_act, nz, nzra, nzrz, ipz1,ipz2, nzz, ip1, ip2, l2 + integer(psb_ipk_), allocatable :: ia(:), ja(:), iz(:), jz(:) + complex(psb_spk_), allocatable :: zw(:), val(:), valz(:) + integer(psb_ipk_), allocatable :: uplevs(:), rowlevs(:), idxs(:) + complex(psb_spk_), allocatable :: row(:) + type(psb_c_coo_sparse_mat) :: trw + type(psb_c_csr_sparse_mat) :: acsr, zcsr + integer(psb_ipk_) :: ktrw, nidx + type(psb_i_heap) :: heap + + real(psb_dpk_) :: alpha + character(len=20) :: name='psb_sp_invk' + + info = psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + + if (.not.(a%is_triangle().and.a%is_unit().and.a%is_upper())) then + write(psb_err_unit,*) 'Wrong A ' + info = psb_err_internal_error_ + call psb_errpush(psb_err_internal_error_,name,a_err='wrong A') + goto 9999 + end if + call a%cp_to(acsr) + call trw%allocate(izero,izero,ione) + if (info == psb_success_) allocate(zw(n),iz(n),valz(n),& + & row(n),rowlevs(n),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + allocate(uplevs(acsr%get_nzeros()),stat=info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='Allocate') + goto 9999 + end if + uplevs(:) = 0 + row(:) = czero + rowlevs(:) = -(n+1) + + call zcsr%allocate(n,n,n*fill_in) + call zcsr%set_triangle() + call zcsr%set_unit(.false.) + call zcsr%set_upper() + call psb_ensure_size(n+1, idxs, info) + + + ! + ! + zcsr%irp(1) = 1 + nzz = 0 + + l2 = 0 + outer: do i = 1, n-1 + ! ZW = e_i + call psb_invk_copyin(i,n,acsr,ione,n,row,rowlevs,heap,ktrw,trw,info,& + & sign=-sone,inlevs=inlevs) + row(i) = cone + rowlevs(i) = 0 + + ! Update loop + call psb_invk_inv(fill_in,i,row,rowlevs,heap,& + & acsr%ja,acsr%irp,acsr%val,uplevs,nidx,idxs,info) + + call psb_invk_copyout(fill_in,i,n,row,rowlevs,nidx,idxs,& + & l2,zcsr%ja,zcsr%irp,zcsr%val,info) + + nzz = l2 + end do outer + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='mainloop') + goto 9999 + end if + ipz1 = nzz+1 + call psb_ensure_size(ipz1,zcsr%val,info) + call psb_ensure_size(ipz1,zcsr%ja,info) + zcsr%val(ipz1) = cone + zcsr%ja(ipz1) = n + zcsr%irp(n+1) = ipz1+1 + call zcsr%set_sorted() + call z%mv_from(zcsr) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_csparse_invk + +subroutine psb_c_invk_copyin(i,m,a,jmin,jmax,row,rowlevs,heap,ktrw,trw,info,sign,inlevs) + + use psb_base_mod + use psb_c_invk_fact_mod, psb_protect_name => psb_c_invk_copyin + + implicit none + + ! Arguments + type(psb_c_csr_sparse_mat), intent(in) :: a + type(psb_c_coo_sparse_mat), intent(inout) :: trw + integer(psb_ipk_), intent(in) :: i,m,jmin,jmax + integer(psb_ipk_), intent(inout) :: ktrw,info + integer(psb_ipk_), intent(inout) :: rowlevs(:) + complex(psb_spk_), intent(inout) :: row(:) + type(psb_i_heap), intent(inout) :: heap + real(psb_spk_), optional, intent(in) :: sign + integer(psb_ipk_), intent(in), optional :: inlevs(:) + + ! Local variables + integer(psb_ipk_) :: k,j,irb,err_act, nz + integer(psb_ipk_), parameter :: nrb=16 + real(psb_spk_) :: sign_ + character(len=20), parameter :: name='invk_copyin' + character(len=20) :: ch_err + + info = psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + + call heap%init(info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_init_heap') + goto 9999 + end if + + if (present(sign)) then + sign_ = sign + else + sign_ = sone + end if + + + ! + ! Take a fast shortcut if the matrix is stored in CSR format + ! + if (present(inlevs)) then + do j = a%irp(i), a%irp(i+1) - 1 + k = a%ja(j) + if ((jmin<=k).and.(k<=jmax)) then + row(k) = sign_ * a%val(j) + rowlevs(k) = inlevs(j) + call heap%insert(k,info) + end if + end do + else + do j = a%irp(i), a%irp(i+1) - 1 + k = a%ja(j) + if ((jmin<=k).and.(k<=jmax)) then + row(k) = sign_ * a%val(j) + rowlevs(k) = 0 + call heap%insert(k,info) + end if + end do + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_c_invk_copyin + + +subroutine psb_c_invk_copyout(fill_in,i,m,row,rowlevs,nidx,idxs,& + & l2,uia1,uia2,uaspk,info) + + use psb_base_mod + use psb_c_invk_fact_mod, psb_protect_name => psb_c_invk_copyout + + implicit none + + ! Arguments + integer(psb_ipk_), intent(in) :: fill_in, i, m, nidx + integer(psb_ipk_), intent(inout) :: l2, info + integer(psb_ipk_), intent(inout) :: rowlevs(:), idxs(:) + integer(psb_ipk_), allocatable, intent(inout) :: uia1(:), uia2(:) + complex(psb_spk_), allocatable, intent(inout) :: uaspk(:) + complex(psb_spk_), intent(inout) :: row(:) + + ! Local variables + integer(psb_ipk_) :: j,isz,err_act,int_err(5),idxp + character(len=20), parameter :: name='psb_ciluk_factint' + character(len=20) :: ch_err + + info = psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + + do idxp=1,nidx + + j = idxs(idxp) + + if (j>=i) then + ! + ! Copy the upper part of the row + ! + if (rowlevs(j) <= fill_in) then + l2 = l2 + 1 + if (size(uaspk) < l2) then + ! + ! Figure out a good reallocation size! + ! + isz = max(int(1.2*l2),l2+100) + call psb_realloc(isz,uaspk,info) + if (info == psb_success_) call psb_realloc(isz,uia1,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='Allocate') + goto 9999 + end if + end if + uia1(l2) = j + uaspk(l2) = row(j) + end if + ! + ! Re-initialize row(j) and rowlevs(j) + ! + row(j) = czero + rowlevs(j) = -(m+1) + end if + end do + + uia2(i+1) = l2 + 1 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_c_invk_copyout + +subroutine psb_cinvk_inv(fill_in,i,row,rowlevs,heap,ja,irp,val,uplevs,& + & nidx,idxs,info) + + use psb_base_mod + use psb_c_invk_fact_mod, psb_protect_name => psb_cinvk_inv + + implicit none + + ! Arguments + type(psb_i_heap), intent(inout) :: heap + integer(psb_ipk_), intent(in) :: i, fill_in + integer(psb_ipk_), intent(inout) :: nidx,info + integer(psb_ipk_), intent(inout) :: rowlevs(:) + integer(psb_ipk_), allocatable, intent(inout) :: idxs(:) + integer(psb_ipk_), intent(in) :: ja(:),irp(:),uplevs(:) + complex(psb_spk_), intent(in) :: val(:) + complex(psb_spk_), intent(inout) :: row(:) + + ! Local variables + integer(psb_ipk_) :: k,j,lrwk,jj,lastk, iret + real(psb_dpk_) :: rwk + + + info = psb_success_ + + call psb_ensure_size(200, idxs, info) + if (info /= psb_success_) return + nidx = 1 + idxs(1) = i + lastk = i + + ! + ! Do while there are indices to be processed + ! + do + ! Beware: (iret < 0) means that the heap is empty, not an error. + call heap%get_first(k,iret) + if (iret < 0) then +!!$ write(psb_err_unit,*) 'IINVK: ',i,' returning at ',lastk + return + end if + + ! + ! Just in case an index has been put on the heap more than once. + ! + if (k == lastk) cycle + + lastk = k + nidx = nidx + 1 + if (nidx>size(idxs)) then + call psb_realloc(nidx+psb_heap_resize,idxs,info) + if (info /= psb_success_) return + end if + idxs(nidx) = k + + if ((row(k) /= czero).and.(rowlevs(k) <= fill_in)) then + ! + ! Note: since U is scaled while copying it out (see iluk_copyout), + ! we can use rwk in the update below + ! + rwk = row(k) + lrwk = rowlevs(k) + + do jj=irp(k),irp(k+1)-1 + j = ja(jj) + if (j<=k) then + info = -i + return + endif + ! + ! Insert the index into the heap for further processing. + ! The fill levels are initialized to a negative value. If we find + ! one, it means that it is an as yet untouched index, so we need + ! to insert it; otherwise it is already on the heap, there is no + ! need to insert it more than once. + ! + if (rowlevs(j)<0) then + call heap%insert(j,info) + if (info /= psb_success_) return + rowlevs(j) = abs(rowlevs(j)) + end if + ! + ! Update row(j) and the corresponding fill level + ! + row(j) = row(j) - rwk * val(jj) + rowlevs(j) = min(rowlevs(j),lrwk+uplevs(jj)+1) + end do + + end if + end do + +end subroutine psb_cinvk_inv diff --git a/prec/impl/psb_c_invt_fact.f90 b/prec/impl/psb_c_invt_fact.f90 new file mode 100644 index 00000000..272cb4ab --- /dev/null +++ b/prec/impl/psb_c_invt_fact.f90 @@ -0,0 +1,631 @@ +! +! 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 AMG4PSBLAS, 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_invt_bld(a,fillin,invfill,thresh,invthresh,& + & lmat,d,umat,desc,info,blck) + + use psb_base_mod + use psb_c_invt_fact_mod, psb_protect_name => psb_c_invt_bld + use psb_c_ilu_fact_mod + implicit none + + ! Arguments + type(psb_cspmat_type), intent(in), target :: a + integer(psb_ipk_), intent(in) :: fillin,invfill + real(psb_spk_), intent(in) :: thresh + real(psb_spk_), intent(in) :: invthresh + type(psb_cspmat_type), intent(inout) :: lmat, umat + complex(psb_spk_), allocatable :: d(:) + Type(psb_desc_type), Intent(inout) :: desc + integer(psb_ipk_), intent(out) :: info + type(psb_cspmat_type), intent(in), optional :: blck + ! + integer(psb_ipk_) :: i, nztota, err_act, n_row, nrow_a, n_col + type(psb_cspmat_type) :: atmp + complex(psb_spk_), allocatable :: pq(:), pd(:), w(:) + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: nzrmax + real(psb_spk_) :: sp_thresh + character(len=20) :: name, ch_err, fname + + + info = psb_success_ + name='psb_cinvt_fact' + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = psb_cd_get_context(desc) + call psb_info(ictxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + ! + ! Check the memory available to hold the incomplete L and U factors + ! and allocate it if needed + ! + nrow_a = a%get_nrows() + nztota = a%get_nzeros() + + if (present(blck)) then + nztota = nztota + blck%get_nzeros() + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & ': out get_nnzeros',nrow_a,nztota,& + & a%get_nrows(),a%get_ncols(),a%get_nzeros() + + + n_row = psb_cd_get_local_rows(desc) + n_col = psb_cd_get_local_cols(desc) + allocate(pd(n_row),w(n_row),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + nzrmax = fillin + sp_thresh = thresh + + call lmat%allocate(n_row,n_row,info,nz=nztota) + if (info == psb_success_) call umat%allocate(n_row,n_row,info,nz=nztota) + + if (info == 0) call psb_ilut_fact(nzrmax,sp_thresh,& + & a,lmat,umat,pd,info,blck=blck,iscale=psb_ilu_scale_maxval_) + + if (info == psb_success_) call atmp%allocate(n_row,n_row,info,nz=nztota) + if(info/=0) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if (.false.) then +!!$ if (debug_level >= psb_debug_inner_) then + write(fname,'(a,i0,a)') 'invt-lo-',me,'.mtx' + call lmat%print(fname,head="INVTLOW") + write(fname,'(a,i0,a)') 'invt-up-',me,'.mtx' + call umat%print(fname,head="INVTUPP") + end if + + ! + ! Compute the approx U^-1 and L^-1 + ! + nzrmax = invfill + call psb_csparse_invt(n_row,umat,atmp,nzrmax,invthresh,info) + if (info == psb_success_) call psb_move_alloc(atmp,umat,info) + if (info == psb_success_) call lmat%transp() + if (info == psb_success_) call psb_csparse_invt(n_row,lmat,atmp,nzrmax,invthresh,info) + if (info == psb_success_) call psb_move_alloc(atmp,lmat,info) + if (info == psb_success_) call lmat%transp() + ! Done. Hopefully.... + + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='invt') + goto 9999 + end if + + call psb_move_alloc(pd,d,info) + call lmat%set_asb() + call lmat%trim() + call umat%set_asb() + call umat%trim() + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_c_invt_bld + +subroutine psb_csparse_invt(n,a,z,nzrmax,sp_thresh,info) + + use psb_base_mod + use psb_c_invt_fact_mod, psb_protect_name => psb_csparse_invt + + implicit none + integer(psb_ipk_), intent(in) :: n + type(psb_cspmat_type), intent(in) :: a + type(psb_cspmat_type), intent(inout) :: z + integer(psb_ipk_), intent(in) :: nzrmax + real(psb_spk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(out) :: info + ! + integer(psb_ipk_) :: i,j,k, err_act, nz, nzra, nzrz, ipz1,ipz2, nzz, ip1, ip2, l2 + integer(psb_ipk_), allocatable :: ia(:), ja(:), iz(:),jz(:) + complex(psb_spk_), allocatable :: zw(:), val(:), valz(:) + integer(psb_ipk_), allocatable :: uplevs(:), rowlevs(:),idxs(:) + complex(psb_spk_), allocatable :: row(:) + type(psb_c_coo_sparse_mat) :: trw + type(psb_c_csr_sparse_mat) :: acsr, zcsr + integer(psb_ipk_) :: ktrw, nidx, nlw,nup,jmaxup + type(psb_i_heap) :: heap + real(psb_spk_) :: alpha, nrmi + character(len=20) :: name='psb_sp_invt' + + info = psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + + if (.not.(a%is_triangle().and.a%is_unit().and.a%is_upper())) then + write(psb_err_unit,*) 'Wrong A ' + info = psb_err_internal_error_ + call psb_errpush(psb_err_internal_error_,name,a_err='wrong A') + goto 9999 + end if + call a%cp_to(acsr) + call trw%allocate(izero,izero,ione) + if (info == psb_success_) allocate(zw(n),iz(n),valz(n),& + & row(n),rowlevs(n),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + call zcsr%allocate(n,n,n*nzrmax) + call zcsr%set_triangle() + call zcsr%set_unit(.false.) + call zcsr%set_upper() + ! + ! + nzz = 0 + row(:) = czero + rowlevs(:) = 0 + l2 = 0 + zcsr%irp(1) = 1 + + outer: do i = 1, n-1 + ! ZW = e_i + call psb_c_invt_copyin(i,n,acsr,i,ione,n,nlw,nup,jmaxup,nrmi,row,& + & heap,rowlevs,ktrw,trw,info,sign=-sone) + if (info /= 0) exit + row(i) = cone + ! Adjust norm + if (nrmi < sone) then + nrmi = sqrt(sone + nrmi**2) + else + nrmi = nrmi*sqrt(cone+cone/(nrmi**2)) + end if + + call psb_invt_inv(sp_thresh,i,nrmi,row,heap,rowlevs,& + & acsr%ja,acsr%irp,acsr%val,nidx,idxs,info) + if (info /= 0) exit +!!$ write(0,*) 'Calling copyout ',nzrmax,nlw,nup,nidx,l2 + call psb_c_invt_copyout(nzrmax,sp_thresh,i,n,nlw,nup,jmaxup,nrmi,row,& + & nidx,idxs,l2,zcsr%ja,zcsr%irp,zcsr%val,info) + if (info /= 0) exit + nzz = l2 + end do outer + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='mainloop') + goto 9999 + end if + + ipz1 = nzz+1 + call psb_ensure_size(ipz1,zcsr%val,info) + call psb_ensure_size(ipz1,zcsr%ja,info) + zcsr%val(ipz1) = cone + zcsr%ja(ipz1) = n + zcsr%irp(n+1) = ipz1+1 + + call z%mv_from(zcsr) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_csparse_invt + +subroutine psb_c_invt_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,nrmi,row,heap,& + & irwt,ktrw,trw,info,sign) + use psb_base_mod + use psb_d_invt_fact_mod, psb_protect_name => psb_d_invt_copyin + implicit none + type(psb_c_csr_sparse_mat), intent(in) :: a + type(psb_c_coo_sparse_mat), intent(inout) :: trw + integer(psb_ipk_), intent(in) :: i, m,jmin,jmax,jd + integer(psb_ipk_), intent(inout) :: ktrw,nlw,nup,jmaxup,info + integer(psb_ipk_), intent(inout) :: irwt(:) + real(psb_spk_), intent(inout) :: nrmi + complex(psb_spk_), intent(inout) :: row(:) + type(psb_i_heap), intent(inout) :: heap + real(psb_spk_), intent(in), optional :: sign + ! + integer(psb_ipk_) :: k,j,irb,kin,nz, err_act + integer(psb_ipk_), parameter :: nrb=16 + real(psb_dpk_) :: dmaxup, sign_ + real(psb_dpk_), external :: dnrm2 + character(len=20), parameter :: name='invt_copyin' + + info = psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + + call heap%init(info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_init_heap') + goto 9999 + end if + sign_ = sone + if (present(sign)) sign_ = sign + ! + ! nrmi is the norm of the current sparse row (for the time being, + ! we use the 2-norm). + ! NOTE: the 2-norm below includes also elements that are outside + ! [jmin:jmax] strictly. Is this really important? TO BE CHECKED. + ! + + nlw = 0 + nup = 0 + jmaxup = 0 + dmaxup = czero + nrmi = czero + + do j = a%irp(i), a%irp(i+1) - 1 + k = a%ja(j) + if ((jmin<=k).and.(k<=jmax)) then + row(k) = sign_ * a%val(j) + call heap%insert(k,info) + irwt(k) = 1 + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_insert_heap') + goto 9999 + end if + end if + if (kjd) then + nup = nup + 1 + if (abs(row(k))>dmaxup) then + jmaxup = k + dmaxup = abs(row(k)) + end if + end if + end do + nz = a%irp(i+1) - a%irp(i) + nrmi = dnrm2(nz,a%val(a%irp(i):),ione) + + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_c_invt_copyin + +subroutine psb_c_invt_copyout(fill_in,thres,i,m,nlw,nup,jmaxup,nrmi,row, & + & nidx,idxs,l2,ja,irp,val,info) + + use psb_base_mod + use psb_c_invt_fact_mod, psb_protect_name => psb_c_invt_copyout + + implicit none + + ! Arguments + integer(psb_ipk_), intent(in) :: fill_in,i,m,nidx,nlw,nup,jmaxup + integer(psb_ipk_), intent(in) :: idxs(:) + integer(psb_ipk_), intent(inout) :: l2, info + integer(psb_ipk_), allocatable, intent(inout) :: ja(:),irp(:) + real(psb_spk_), intent(in) :: thres,nrmi + complex(psb_spk_),allocatable, intent(inout) :: val(:) + complex(psb_spk_), intent(inout) :: row(:) + + ! Local variables + real(psb_dpk_),allocatable :: xw(:) + integer(psb_ipk_), allocatable :: xwid(:), indx(:) + real(psb_dpk_) :: witem, wmin + integer(psb_ipk_) :: widx + integer(psb_ipk_) :: k,isz,err_act,int_err(5),idxp, nz + type(psb_d_idx_heap) :: heap + character(len=20), parameter :: name='invt_copyout' + character(len=20) :: ch_err + logical :: fndmaxup + + info = psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + + ! + ! Here we need to apply also the dropping rule base on the fill-in. + ! We do it by putting into a heap the elements that are not dropped + ! by using the 2-norm rule, and then copying them out. + ! + ! The heap goes down on the entry absolute value, so the first item + ! is the largest absolute value. + ! +!!$ write(0,*) 'invt_copyout ',nidx,nup+fill_in + call heap%init(info,dir=psb_asort_down_) + + if (info == psb_success_) allocate(xwid(nidx),xw(nidx),indx(nidx),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/3*nidx/),& + & a_err='real(psb_dpk_)') + goto 9999 + end if + + ! + ! First the lower part + ! + + nz = 0 + idxp = 0 + + do + + idxp = idxp + 1 + if (idxp > nidx) exit + if (idxs(idxp) >= i) exit + widx = idxs(idxp) + witem = row(widx) + ! + ! Dropping rule based on the 2-norm + ! + if (abs(witem) < thres*nrmi) cycle + + nz = nz + 1 + xw(nz) = witem + xwid(nz) = widx + call heap%insert(witem,widx,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_insert_heap') + goto 9999 + end if + end do + + if (nz > 1) then + write(psb_err_unit,*) 'Warning: lower triangle from invt???? ' + end if + + + if (idxp <= size(idxs)) then + if (idxs(idxp) < i) then + do + idxp = idxp + 1 + if (idxp > nidx) exit + if (idxs(idxp) >= i) exit + end do + end if + end if + idxp = idxp - 1 + nz = 0 + wmin=HUGE(wmin) + if (.false.) then + do + + idxp = idxp + 1 + if (idxp > nidx) exit + widx = idxs(idxp) + if (widx < i) then + write(psb_err_unit,*) 'Warning: lower triangle in upper copy',widx,i,idxp,idxs(idxp) + cycle + end if + if (widx > m) then + cycle + end if + witem = row(widx) + ! + ! Dropping rule based on the 2-norm. But keep the jmaxup-th entry anyway. + ! + if ((widx /= jmaxup) .and. (widx /= i) .and. (abs(witem) < thres*nrmi)) then + cycle + end if + if ((widx/=jmaxup).and.(nz > nup+fill_in)) then + if (abs(witem) < wmin) cycle + endif + wmin = min(abs(witem),wmin) + nz = nz + 1 + xw(nz) = witem + xwid(nz) = widx + call heap%insert(witem,widx,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_insert_heap') + goto 9999 + end if + + end do + + ! + ! Now we have to take out the first nup-fill_in entries. But make sure + ! we include entry jmaxup. + ! + if (nz <= nup+fill_in) then + ! + ! Just copy everything from xw + ! + fndmaxup=.true. + else + fndmaxup = .false. + nz = nup+fill_in + do k=1,nz + call heap%get_first(witem,widx,info) + xw(k) = witem + xwid(k) = widx + if (widx == jmaxup) fndmaxup=.true. + end do + end if + if ((i nidx) exit + widx = idxs(idxp) + if (widx < i) then + write(psb_err_unit,*) 'Warning: lower triangle in upper copy',widx,i,idxp,idxs(idxp) + cycle + end if + if (widx > m) then + cycle + end if + witem = row(widx) + ! + ! Dropping rule based on the 2-norm. But keep the jmaxup-th entry anyway. + ! + if ((widx /= i) .and. (abs(witem) < thres*nrmi)) then + cycle + end if + if (nz > nup+fill_in) then + if (abs(witem) < wmin) cycle + endif + wmin = min(abs(witem),wmin) + nz = nz + 1 + xw(nz) = witem + xwid(nz) = widx + call heap%insert(witem,widx,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_insert_heap') + goto 9999 + end if + + end do + + ! + ! Now we have to take out the first nup-fill_in entries. But make sure + ! we include entry jmaxup. + ! + if (nz > nup+fill_in) then + nz = nup+fill_in + do k=1,nz + call heap%get_first(witem,widx,info) + xw(k) = witem + xwid(k) = widx + end do + end if + end if + + ! + ! Now we put things back into ascending column order + ! + call psb_msort(xwid(1:nz),indx(1:nz),dir=psb_sort_up_) + + ! + ! Copy out the upper part of the row + ! + do k=1,nz + l2 = l2 + 1 + if (size(val) < l2) then + ! + ! Figure out a good reallocation size! + ! + isz = max(int(1.2*l2),l2+100) + call psb_realloc(isz,val,info) + if (info == psb_success_) call psb_realloc(isz,ja,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='Allocate') + goto 9999 + end if + end if + ja(l2) = xwid(k) + val(l2) = xw(indx(k)) + end do + + ! + ! Set row to zero + ! + do idxp=1,nidx + row(idxs(idxp)) = czero + end do + + irp(i+1) = l2 + 1 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_c_invt_copyout diff --git a/prec/impl/psb_d_invk_fact.f90 b/prec/impl/psb_d_invk_fact.f90 new file mode 100644 index 00000000..e0dc598c --- /dev/null +++ b/prec/impl/psb_d_invk_fact.f90 @@ -0,0 +1,494 @@ +! +! +! 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_invk_bld(a,fill1, fill2,lmat,d,umat,desc,info,blck) + + use psb_base_mod + use psb_d_invk_fact_mod, psb_protect_name => psb_d_invk_bld + use psb_d_ilu_fact_mod + implicit none + + ! Arguments + type(psb_dspmat_type), intent(in), target :: a + integer(psb_ipk_), intent(in) :: fill1, fill2 + type(psb_dspmat_type), intent(inout) :: lmat, umat + real(psb_dpk_), allocatable :: d(:) + Type(psb_desc_type), Intent(inout) :: desc + integer(psb_ipk_), intent(out) :: info + type(psb_dspmat_type), intent(in), optional :: blck + ! + integer(psb_ipk_) :: i, nztota, err_act, n_row, nrow_a, n_col + type(psb_dspmat_type) :: atmp + real(psb_dpk_), allocatable :: pq(:), pd(:) + integer(psb_ipk_), allocatable :: uplevs(:) + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: nzrmax + character(len=20) :: name, ch_err + + + info = psb_success_ + name='psb_dinvk_bld' + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = psb_cd_get_context(desc) + call psb_info(ictxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + ! + ! Check the memory available to hold the incomplete L and U factors + ! and allocate it if needed + ! + nrow_a = a%get_nrows() + nztota = a%get_nzeros() + + if (present(blck)) then + nztota = nztota + blck%get_nzeros() + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & ': out get_nnzeros',nrow_a,nztota,& + & a%get_nrows(),a%get_ncols(),a%get_nzeros() + + + n_row = psb_cd_get_local_rows(desc) + n_col = psb_cd_get_local_cols(desc) + allocate(pd(n_row),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + + call lmat%allocate(n_row,n_row,info,nz=nztota) + if (info == psb_success_) call umat%allocate(n_row,n_row,info,nz=nztota) + + + call psb_iluk_fact(fill1,psb_ilu_n_,a,lmat,umat,pd,info,blck=blck) + !,uplevs=uplevs) + !call psb_diluk_fact(fill1,psb_ilu_n_,a,lmat,umat,pd,info,blck=blck) + + if (info == psb_success_) call atmp%allocate(n_row,n_row,info,nz=nztota) + if(info/=0) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + + ! + ! Compute the aprox U^-1 and L^-1 + ! + call psb_sparse_invk(n_row,umat,atmp,fill2,info) + if (info == psb_success_) call psb_move_alloc(atmp,umat,info) + if (info == psb_success_) call lmat%transp() + if (info == psb_success_) call psb_sparse_invk(n_row,lmat,atmp,fill2,info) + if (info == psb_success_) call psb_move_alloc(atmp,lmat,info) + if (info == psb_success_) call lmat%transp() + ! Done. Hopefully.... + + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='invt') + goto 9999 + end if + + call psb_move_alloc(pd,d,info) + call lmat%set_asb() + call lmat%trim() + call umat%set_asb() + call umat%trim() + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_d_invk_bld + +subroutine psb_dsparse_invk(n,a,z,fill_in,info,inlevs) + + use psb_base_mod + use psb_d_invk_fact_mod, psb_protect_name => psb_dsparse_invk + + integer(psb_ipk_), intent(in) :: n + type(psb_dspmat_type), intent(in) :: a + type(psb_dspmat_type), intent(inout) :: z + integer(psb_ipk_), intent(in) :: fill_in + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: inlevs(:) + ! + integer(psb_ipk_) :: i,j,k, err_act, nz, nzra, nzrz, ipz1,ipz2, nzz, ip1, ip2, l2 + integer(psb_ipk_), allocatable :: ia(:), ja(:), iz(:), jz(:) + real(psb_dpk_), allocatable :: zw(:), val(:), valz(:) + integer(psb_ipk_), allocatable :: uplevs(:), rowlevs(:), idxs(:) + real(psb_dpk_), allocatable :: row(:) + type(psb_d_coo_sparse_mat) :: trw + type(psb_d_csr_sparse_mat) :: acsr, zcsr + integer(psb_ipk_) :: ktrw, nidx + type(psb_i_heap) :: heap + + real(psb_dpk_) :: alpha + character(len=20) :: name='psb_sp_invk' + + info = psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + + if (.not.(a%is_triangle().and.a%is_unit().and.a%is_upper())) then + write(psb_err_unit,*) 'Wrong A ' + info = psb_err_internal_error_ + call psb_errpush(psb_err_internal_error_,name,a_err='wrong A') + goto 9999 + end if + call a%cp_to(acsr) + call trw%allocate(izero,izero,ione) + if (info == psb_success_) allocate(zw(n),iz(n),valz(n),& + & row(n),rowlevs(n),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + allocate(uplevs(acsr%get_nzeros()),stat=info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='Allocate') + goto 9999 + end if + uplevs(:) = 0 + row(:) = dzero + rowlevs(:) = -(n+1) + + call zcsr%allocate(n,n,n*fill_in) + call zcsr%set_triangle() + call zcsr%set_unit(.false.) + call zcsr%set_upper() + call psb_ensure_size(n+1, idxs, info) + + + ! + ! + zcsr%irp(1) = 1 + nzz = 0 + + l2 = 0 + outer: do i = 1, n-1 + ! ZW = e_i + call psb_invk_copyin(i,n,acsr,ione,n,row,rowlevs,heap,ktrw,trw,info,& + & sign=-done,inlevs=inlevs) + row(i) = done + rowlevs(i) = 0 + + ! Update loop + call psb_invk_inv(fill_in,i,row,rowlevs,heap,& + & acsr%ja,acsr%irp,acsr%val,uplevs,nidx,idxs,info) + + call psb_invk_copyout(fill_in,i,n,row,rowlevs,nidx,idxs,& + & l2,zcsr%ja,zcsr%irp,zcsr%val,info) + + nzz = l2 + end do outer + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='mainloop') + goto 9999 + end if + ipz1 = nzz+1 + call psb_ensure_size(ipz1,zcsr%val,info) + call psb_ensure_size(ipz1,zcsr%ja,info) + zcsr%val(ipz1) = done + zcsr%ja(ipz1) = n + zcsr%irp(n+1) = ipz1+1 + call zcsr%set_sorted() + call z%mv_from(zcsr) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_dsparse_invk + +subroutine psb_d_invk_copyin(i,m,a,jmin,jmax,row,rowlevs,heap,ktrw,trw,info,sign,inlevs) + + use psb_base_mod + use psb_d_invk_fact_mod, psb_protect_name => psb_d_invk_copyin + + implicit none + + ! Arguments + type(psb_d_csr_sparse_mat), intent(in) :: a + type(psb_d_coo_sparse_mat), intent(inout) :: trw + integer(psb_ipk_), intent(in) :: i,m,jmin,jmax + integer(psb_ipk_), intent(inout) :: ktrw,info + integer(psb_ipk_), intent(inout) :: rowlevs(:) + real(psb_dpk_), intent(inout) :: row(:) + type(psb_i_heap), intent(inout) :: heap + real(psb_dpk_), optional, intent(in) :: sign + integer(psb_ipk_), intent(in), optional :: inlevs(:) + + ! Local variables + integer(psb_ipk_) :: k,j,irb,err_act, nz + integer(psb_ipk_), parameter :: nrb=16 + real(psb_dpk_) :: sign_ + character(len=20), parameter :: name='invk_copyin' + character(len=20) :: ch_err + + info = psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + + call heap%init(info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_init_heap') + goto 9999 + end if + + if (present(sign)) then + sign_ = sign + else + sign_ = done + end if + + + ! + ! Take a fast shortcut if the matrix is stored in CSR format + ! + if (present(inlevs)) then + do j = a%irp(i), a%irp(i+1) - 1 + k = a%ja(j) + if ((jmin<=k).and.(k<=jmax)) then + row(k) = sign_ * a%val(j) + rowlevs(k) = inlevs(j) + call heap%insert(k,info) + end if + end do + else + do j = a%irp(i), a%irp(i+1) - 1 + k = a%ja(j) + if ((jmin<=k).and.(k<=jmax)) then + row(k) = sign_ * a%val(j) + rowlevs(k) = 0 + call heap%insert(k,info) + end if + end do + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_d_invk_copyin + + +subroutine psb_d_invk_copyout(fill_in,i,m,row,rowlevs,nidx,idxs,& + & l2,uia1,uia2,uaspk,info) + + use psb_base_mod + use psb_d_invk_fact_mod, psb_protect_name => psb_d_invk_copyout + + implicit none + + ! Arguments + integer(psb_ipk_), intent(in) :: fill_in, i, m, nidx + integer(psb_ipk_), intent(inout) :: l2, info + integer(psb_ipk_), intent(inout) :: rowlevs(:), idxs(:) + integer(psb_ipk_), allocatable, intent(inout) :: uia1(:), uia2(:) + real(psb_dpk_), allocatable, intent(inout) :: uaspk(:) + real(psb_dpk_), intent(inout) :: row(:) + + ! Local variables + integer(psb_ipk_) :: j,isz,err_act,int_err(5),idxp + character(len=20), parameter :: name='psb_diluk_factint' + character(len=20) :: ch_err + + info = psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + + do idxp=1,nidx + + j = idxs(idxp) + + if (j>=i) then + ! + ! Copy the upper part of the row + ! + if (rowlevs(j) <= fill_in) then + l2 = l2 + 1 + if (size(uaspk) < l2) then + ! + ! Figure out a good reallocation size! + ! + isz = max(int(1.2*l2),l2+100) + call psb_realloc(isz,uaspk,info) + if (info == psb_success_) call psb_realloc(isz,uia1,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='Allocate') + goto 9999 + end if + end if + uia1(l2) = j + uaspk(l2) = row(j) + end if + ! + ! Re-initialize row(j) and rowlevs(j) + ! + row(j) = dzero + rowlevs(j) = -(m+1) + end if + end do + + uia2(i+1) = l2 + 1 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_d_invk_copyout + +subroutine psb_dinvk_inv(fill_in,i,row,rowlevs,heap,ja,irp,val,uplevs,& + & nidx,idxs,info) + + use psb_base_mod + use psb_d_invk_fact_mod, psb_protect_name => psb_dinvk_inv + + implicit none + + ! Arguments + type(psb_i_heap), intent(inout) :: heap + integer(psb_ipk_), intent(in) :: i, fill_in + integer(psb_ipk_), intent(inout) :: nidx,info + integer(psb_ipk_), intent(inout) :: rowlevs(:) + integer(psb_ipk_), allocatable, intent(inout) :: idxs(:) + integer(psb_ipk_), intent(in) :: ja(:),irp(:),uplevs(:) + real(psb_dpk_), intent(in) :: val(:) + real(psb_dpk_), intent(inout) :: row(:) + + ! Local variables + integer(psb_ipk_) :: k,j,lrwk,jj,lastk, iret + real(psb_dpk_) :: rwk + + + info = psb_success_ + + call psb_ensure_size(200, idxs, info) + if (info /= psb_success_) return + nidx = 1 + idxs(1) = i + lastk = i + + ! + ! Do while there are indices to be processed + ! + do + ! Beware: (iret < 0) means that the heap is empty, not an error. + call heap%get_first(k,iret) + if (iret < 0) then +!!$ write(psb_err_unit,*) 'IINVK: ',i,' returning at ',lastk + return + end if + + ! + ! Just in case an index has been put on the heap more than once. + ! + if (k == lastk) cycle + + lastk = k + nidx = nidx + 1 + if (nidx>size(idxs)) then + call psb_realloc(nidx+psb_heap_resize,idxs,info) + if (info /= psb_success_) return + end if + idxs(nidx) = k + + if ((row(k) /= dzero).and.(rowlevs(k) <= fill_in)) then + ! + ! Note: since U is scaled while copying it out (see iluk_copyout), + ! we can use rwk in the update below + ! + rwk = row(k) + lrwk = rowlevs(k) + + do jj=irp(k),irp(k+1)-1 + j = ja(jj) + if (j<=k) then + info = -i + return + endif + ! + ! Insert the index into the heap for further processing. + ! The fill levels are initialized to a negative value. If we find + ! one, it means that it is an as yet untouched index, so we need + ! to insert it; otherwise it is already on the heap, there is no + ! need to insert it more than once. + ! + if (rowlevs(j)<0) then + call heap%insert(j,info) + if (info /= psb_success_) return + rowlevs(j) = abs(rowlevs(j)) + end if + ! + ! Update row(j) and the corresponding fill level + ! + row(j) = row(j) - rwk * val(jj) + rowlevs(j) = min(rowlevs(j),lrwk+uplevs(jj)+1) + end do + + end if + end do + +end subroutine psb_dinvk_inv diff --git a/prec/impl/psb_d_invt_fact.f90 b/prec/impl/psb_d_invt_fact.f90 new file mode 100644 index 00000000..0312733c --- /dev/null +++ b/prec/impl/psb_d_invt_fact.f90 @@ -0,0 +1,631 @@ +! +! 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 AMG4PSBLAS, 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_invt_bld(a,fillin,invfill,thresh,invthresh,& + & lmat,d,umat,desc,info,blck) + + use psb_base_mod + use psb_d_invt_fact_mod, psb_protect_name => psb_d_invt_bld + use psb_d_ilu_fact_mod + implicit none + + ! Arguments + type(psb_dspmat_type), intent(in), target :: a + integer(psb_ipk_), intent(in) :: fillin,invfill + real(psb_dpk_), intent(in) :: thresh + real(psb_dpk_), intent(in) :: invthresh + type(psb_dspmat_type), intent(inout) :: lmat, umat + real(psb_dpk_), allocatable :: d(:) + Type(psb_desc_type), Intent(inout) :: desc + integer(psb_ipk_), intent(out) :: info + type(psb_dspmat_type), intent(in), optional :: blck + ! + integer(psb_ipk_) :: i, nztota, err_act, n_row, nrow_a, n_col + type(psb_dspmat_type) :: atmp + real(psb_dpk_), allocatable :: pq(:), pd(:), w(:) + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: nzrmax + real(psb_dpk_) :: sp_thresh + character(len=20) :: name, ch_err, fname + + + info = psb_success_ + name='psb_dinvt_fact' + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = psb_cd_get_context(desc) + call psb_info(ictxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + ! + ! Check the memory available to hold the incomplete L and U factors + ! and allocate it if needed + ! + nrow_a = a%get_nrows() + nztota = a%get_nzeros() + + if (present(blck)) then + nztota = nztota + blck%get_nzeros() + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & ': out get_nnzeros',nrow_a,nztota,& + & a%get_nrows(),a%get_ncols(),a%get_nzeros() + + + n_row = psb_cd_get_local_rows(desc) + n_col = psb_cd_get_local_cols(desc) + allocate(pd(n_row),w(n_row),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + nzrmax = fillin + sp_thresh = thresh + + call lmat%allocate(n_row,n_row,info,nz=nztota) + if (info == psb_success_) call umat%allocate(n_row,n_row,info,nz=nztota) + + if (info == 0) call psb_ilut_fact(nzrmax,sp_thresh,& + & a,lmat,umat,pd,info,blck=blck,iscale=psb_ilu_scale_maxval_) + + if (info == psb_success_) call atmp%allocate(n_row,n_row,info,nz=nztota) + if(info/=0) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if (.false.) then +!!$ if (debug_level >= psb_debug_inner_) then + write(fname,'(a,i0,a)') 'invt-lo-',me,'.mtx' + call lmat%print(fname,head="INVTLOW") + write(fname,'(a,i0,a)') 'invt-up-',me,'.mtx' + call umat%print(fname,head="INVTUPP") + end if + + ! + ! Compute the approx U^-1 and L^-1 + ! + nzrmax = invfill + call psb_dsparse_invt(n_row,umat,atmp,nzrmax,invthresh,info) + if (info == psb_success_) call psb_move_alloc(atmp,umat,info) + if (info == psb_success_) call lmat%transp() + if (info == psb_success_) call psb_dsparse_invt(n_row,lmat,atmp,nzrmax,invthresh,info) + if (info == psb_success_) call psb_move_alloc(atmp,lmat,info) + if (info == psb_success_) call lmat%transp() + ! Done. Hopefully.... + + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='invt') + goto 9999 + end if + + call psb_move_alloc(pd,d,info) + call lmat%set_asb() + call lmat%trim() + call umat%set_asb() + call umat%trim() + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_d_invt_bld + +subroutine psb_dsparse_invt(n,a,z,nzrmax,sp_thresh,info) + + use psb_base_mod + use psb_d_invt_fact_mod, psb_protect_name => psb_dsparse_invt + + implicit none + integer(psb_ipk_), intent(in) :: n + type(psb_dspmat_type), intent(in) :: a + type(psb_dspmat_type), intent(inout) :: z + integer(psb_ipk_), intent(in) :: nzrmax + real(psb_dpk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(out) :: info + ! + integer(psb_ipk_) :: i,j,k, err_act, nz, nzra, nzrz, ipz1,ipz2, nzz, ip1, ip2, l2 + integer(psb_ipk_), allocatable :: ia(:), ja(:), iz(:),jz(:) + real(psb_dpk_), allocatable :: zw(:), val(:), valz(:) + integer(psb_ipk_), allocatable :: uplevs(:), rowlevs(:),idxs(:) + real(psb_dpk_), allocatable :: row(:) + type(psb_d_coo_sparse_mat) :: trw + type(psb_d_csr_sparse_mat) :: acsr, zcsr + integer(psb_ipk_) :: ktrw, nidx, nlw,nup,jmaxup + type(psb_i_heap) :: heap + real(psb_dpk_) :: alpha, nrmi + character(len=20) :: name='psb_sp_invt' + + info = psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + + if (.not.(a%is_triangle().and.a%is_unit().and.a%is_upper())) then + write(psb_err_unit,*) 'Wrong A ' + info = psb_err_internal_error_ + call psb_errpush(psb_err_internal_error_,name,a_err='wrong A') + goto 9999 + end if + call a%cp_to(acsr) + call trw%allocate(izero,izero,ione) + if (info == psb_success_) allocate(zw(n),iz(n),valz(n),& + & row(n),rowlevs(n),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + call zcsr%allocate(n,n,n*nzrmax) + call zcsr%set_triangle() + call zcsr%set_unit(.false.) + call zcsr%set_upper() + ! + ! + nzz = 0 + row(:) = dzero + rowlevs(:) = 0 + l2 = 0 + zcsr%irp(1) = 1 + + outer: do i = 1, n-1 + ! ZW = e_i + call psb_d_invt_copyin(i,n,acsr,i,ione,n,nlw,nup,jmaxup,nrmi,row,& + & heap,rowlevs,ktrw,trw,info,sign=-done) + if (info /= 0) exit + row(i) = done + ! Adjust norm + if (nrmi < done) then + nrmi = sqrt(done + nrmi**2) + else + nrmi = nrmi*sqrt(done+done/(nrmi**2)) + end if + + call psb_invt_inv(sp_thresh,i,nrmi,row,heap,rowlevs,& + & acsr%ja,acsr%irp,acsr%val,nidx,idxs,info) + if (info /= 0) exit +!!$ write(0,*) 'Calling copyout ',nzrmax,nlw,nup,nidx,l2 + call psb_d_invt_copyout(nzrmax,sp_thresh,i,n,nlw,nup,jmaxup,nrmi,row,& + & nidx,idxs,l2,zcsr%ja,zcsr%irp,zcsr%val,info) + if (info /= 0) exit + nzz = l2 + end do outer + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='mainloop') + goto 9999 + end if + + ipz1 = nzz+1 + call psb_ensure_size(ipz1,zcsr%val,info) + call psb_ensure_size(ipz1,zcsr%ja,info) + zcsr%val(ipz1) = done + zcsr%ja(ipz1) = n + zcsr%irp(n+1) = ipz1+1 + + call z%mv_from(zcsr) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_dsparse_invt + +subroutine psb_d_invt_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,nrmi,row,heap,& + & irwt,ktrw,trw,info,sign) + use psb_base_mod + use psb_d_invt_fact_mod, psb_protect_name => psb_d_invt_copyin + implicit none + type(psb_d_csr_sparse_mat), intent(in) :: a + type(psb_d_coo_sparse_mat), intent(inout) :: trw + integer(psb_ipk_), intent(in) :: i, m,jmin,jmax,jd + integer(psb_ipk_), intent(inout) :: ktrw,nlw,nup,jmaxup,info + integer(psb_ipk_), intent(inout) :: irwt(:) + real(psb_dpk_), intent(inout) :: nrmi + real(psb_dpk_), intent(inout) :: row(:) + type(psb_i_heap), intent(inout) :: heap + real(psb_dpk_), intent(in), optional :: sign + ! + integer(psb_ipk_) :: k,j,irb,kin,nz, err_act + integer(psb_ipk_), parameter :: nrb=16 + real(psb_dpk_) :: dmaxup, sign_ + real(psb_dpk_), external :: dnrm2 + character(len=20), parameter :: name='invt_copyin' + + info = psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + + call heap%init(info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_init_heap') + goto 9999 + end if + sign_ = done + if (present(sign)) sign_ = sign + ! + ! nrmi is the norm of the current sparse row (for the time being, + ! we use the 2-norm). + ! NOTE: the 2-norm below includes also elements that are outside + ! [jmin:jmax] strictly. Is this really important? TO BE CHECKED. + ! + + nlw = 0 + nup = 0 + jmaxup = 0 + dmaxup = dzero + nrmi = dzero + + do j = a%irp(i), a%irp(i+1) - 1 + k = a%ja(j) + if ((jmin<=k).and.(k<=jmax)) then + row(k) = sign_ * a%val(j) + call heap%insert(k,info) + irwt(k) = 1 + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_insert_heap') + goto 9999 + end if + end if + if (kjd) then + nup = nup + 1 + if (abs(row(k))>dmaxup) then + jmaxup = k + dmaxup = abs(row(k)) + end if + end if + end do + nz = a%irp(i+1) - a%irp(i) + nrmi = dnrm2(nz,a%val(a%irp(i):),ione) + + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_d_invt_copyin + +subroutine psb_d_invt_copyout(fill_in,thres,i,m,nlw,nup,jmaxup,nrmi,row, & + & nidx,idxs,l2,ja,irp,val,info) + + use psb_base_mod + use psb_d_invt_fact_mod, psb_protect_name => psb_d_invt_copyout + + implicit none + + ! Arguments + integer(psb_ipk_), intent(in) :: fill_in,i,m,nidx,nlw,nup,jmaxup + integer(psb_ipk_), intent(in) :: idxs(:) + integer(psb_ipk_), intent(inout) :: l2, info + integer(psb_ipk_), allocatable, intent(inout) :: ja(:),irp(:) + real(psb_dpk_), intent(in) :: thres,nrmi + real(psb_dpk_),allocatable, intent(inout) :: val(:) + real(psb_dpk_), intent(inout) :: row(:) + + ! Local variables + real(psb_dpk_),allocatable :: xw(:) + integer(psb_ipk_), allocatable :: xwid(:), indx(:) + real(psb_dpk_) :: witem, wmin + integer(psb_ipk_) :: widx + integer(psb_ipk_) :: k,isz,err_act,int_err(5),idxp, nz + type(psb_d_idx_heap) :: heap + character(len=20), parameter :: name='invt_copyout' + character(len=20) :: ch_err + logical :: fndmaxup + + info = psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + + ! + ! Here we need to apply also the dropping rule base on the fill-in. + ! We do it by putting into a heap the elements that are not dropped + ! by using the 2-norm rule, and then copying them out. + ! + ! The heap goes down on the entry absolute value, so the first item + ! is the largest absolute value. + ! +!!$ write(0,*) 'invt_copyout ',nidx,nup+fill_in + call heap%init(info,dir=psb_asort_down_) + + if (info == psb_success_) allocate(xwid(nidx),xw(nidx),indx(nidx),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/3*nidx/),& + & a_err='real(psb_dpk_)') + goto 9999 + end if + + ! + ! First the lower part + ! + + nz = 0 + idxp = 0 + + do + + idxp = idxp + 1 + if (idxp > nidx) exit + if (idxs(idxp) >= i) exit + widx = idxs(idxp) + witem = row(widx) + ! + ! Dropping rule based on the 2-norm + ! + if (abs(witem) < thres*nrmi) cycle + + nz = nz + 1 + xw(nz) = witem + xwid(nz) = widx + call heap%insert(witem,widx,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_insert_heap') + goto 9999 + end if + end do + + if (nz > 1) then + write(psb_err_unit,*) 'Warning: lower triangle from invt???? ' + end if + + + if (idxp <= size(idxs)) then + if (idxs(idxp) < i) then + do + idxp = idxp + 1 + if (idxp > nidx) exit + if (idxs(idxp) >= i) exit + end do + end if + end if + idxp = idxp - 1 + nz = 0 + wmin=HUGE(wmin) + if (.false.) then + do + + idxp = idxp + 1 + if (idxp > nidx) exit + widx = idxs(idxp) + if (widx < i) then + write(psb_err_unit,*) 'Warning: lower triangle in upper copy',widx,i,idxp,idxs(idxp) + cycle + end if + if (widx > m) then + cycle + end if + witem = row(widx) + ! + ! Dropping rule based on the 2-norm. But keep the jmaxup-th entry anyway. + ! + if ((widx /= jmaxup) .and. (widx /= i) .and. (abs(witem) < thres*nrmi)) then + cycle + end if + if ((widx/=jmaxup).and.(nz > nup+fill_in)) then + if (abs(witem) < wmin) cycle + endif + wmin = min(abs(witem),wmin) + nz = nz + 1 + xw(nz) = witem + xwid(nz) = widx + call heap%insert(witem,widx,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_insert_heap') + goto 9999 + end if + + end do + + ! + ! Now we have to take out the first nup-fill_in entries. But make sure + ! we include entry jmaxup. + ! + if (nz <= nup+fill_in) then + ! + ! Just copy everything from xw + ! + fndmaxup=.true. + else + fndmaxup = .false. + nz = nup+fill_in + do k=1,nz + call heap%get_first(witem,widx,info) + xw(k) = witem + xwid(k) = widx + if (widx == jmaxup) fndmaxup=.true. + end do + end if + if ((i nidx) exit + widx = idxs(idxp) + if (widx < i) then + write(psb_err_unit,*) 'Warning: lower triangle in upper copy',widx,i,idxp,idxs(idxp) + cycle + end if + if (widx > m) then + cycle + end if + witem = row(widx) + ! + ! Dropping rule based on the 2-norm. But keep the jmaxup-th entry anyway. + ! + if ((widx /= i) .and. (abs(witem) < thres*nrmi)) then + cycle + end if + if (nz > nup+fill_in) then + if (abs(witem) < wmin) cycle + endif + wmin = min(abs(witem),wmin) + nz = nz + 1 + xw(nz) = witem + xwid(nz) = widx + call heap%insert(witem,widx,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_insert_heap') + goto 9999 + end if + + end do + + ! + ! Now we have to take out the first nup-fill_in entries. But make sure + ! we include entry jmaxup. + ! + if (nz > nup+fill_in) then + nz = nup+fill_in + do k=1,nz + call heap%get_first(witem,widx,info) + xw(k) = witem + xwid(k) = widx + end do + end if + end if + + ! + ! Now we put things back into ascending column order + ! + call psb_msort(xwid(1:nz),indx(1:nz),dir=psb_sort_up_) + + ! + ! Copy out the upper part of the row + ! + do k=1,nz + l2 = l2 + 1 + if (size(val) < l2) then + ! + ! Figure out a good reallocation size! + ! + isz = max(int(1.2*l2),l2+100) + call psb_realloc(isz,val,info) + if (info == psb_success_) call psb_realloc(isz,ja,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='Allocate') + goto 9999 + end if + end if + ja(l2) = xwid(k) + val(l2) = xw(indx(k)) + end do + + ! + ! Set row to zero + ! + do idxp=1,nidx + row(idxs(idxp)) = dzero + end do + + irp(i+1) = l2 + 1 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_d_invt_copyout diff --git a/prec/impl/psb_s_invk_fact.f90 b/prec/impl/psb_s_invk_fact.f90 new file mode 100644 index 00000000..1f30c978 --- /dev/null +++ b/prec/impl/psb_s_invk_fact.f90 @@ -0,0 +1,494 @@ +! +! +! 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_invk_bld(a,fill1, fill2,lmat,d,umat,desc,info,blck) + + use psb_base_mod + use psb_s_invk_fact_mod, psb_protect_name => psb_s_invk_bld + use psb_s_ilu_fact_mod + implicit none + + ! Arguments + type(psb_sspmat_type), intent(in), target :: a + integer(psb_ipk_), intent(in) :: fill1, fill2 + type(psb_sspmat_type), intent(inout) :: lmat, umat + real(psb_spk_), allocatable :: d(:) + Type(psb_desc_type), Intent(inout) :: desc + integer(psb_ipk_), intent(out) :: info + type(psb_sspmat_type), intent(in), optional :: blck + ! + integer(psb_ipk_) :: i, nztota, err_act, n_row, nrow_a, n_col + type(psb_sspmat_type) :: atmp + real(psb_spk_), allocatable :: pq(:), pd(:) + integer(psb_ipk_), allocatable :: uplevs(:) + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: nzrmax + character(len=20) :: name, ch_err + + + info = psb_success_ + name='psb_sinvk_bld' + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = psb_cd_get_context(desc) + call psb_info(ictxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + ! + ! Check the memory available to hold the incomplete L and U factors + ! and allocate it if needed + ! + nrow_a = a%get_nrows() + nztota = a%get_nzeros() + + if (present(blck)) then + nztota = nztota + blck%get_nzeros() + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & ': out get_nnzeros',nrow_a,nztota,& + & a%get_nrows(),a%get_ncols(),a%get_nzeros() + + + n_row = psb_cd_get_local_rows(desc) + n_col = psb_cd_get_local_cols(desc) + allocate(pd(n_row),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + + call lmat%allocate(n_row,n_row,info,nz=nztota) + if (info == psb_success_) call umat%allocate(n_row,n_row,info,nz=nztota) + + + call psb_iluk_fact(fill1,psb_ilu_n_,a,lmat,umat,pd,info,blck=blck) + !,uplevs=uplevs) + !call psb_siluk_fact(fill1,psb_ilu_n_,a,lmat,umat,pd,info,blck=blck) + + if (info == psb_success_) call atmp%allocate(n_row,n_row,info,nz=nztota) + if(info/=0) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + + ! + ! Compute the aprox U^-1 and L^-1 + ! + call psb_sparse_invk(n_row,umat,atmp,fill2,info) + if (info == psb_success_) call psb_move_alloc(atmp,umat,info) + if (info == psb_success_) call lmat%transp() + if (info == psb_success_) call psb_sparse_invk(n_row,lmat,atmp,fill2,info) + if (info == psb_success_) call psb_move_alloc(atmp,lmat,info) + if (info == psb_success_) call lmat%transp() + ! Done. Hopefully.... + + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='invt') + goto 9999 + end if + + call psb_move_alloc(pd,d,info) + call lmat%set_asb() + call lmat%trim() + call umat%set_asb() + call umat%trim() + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_s_invk_bld + +subroutine psb_ssparse_invk(n,a,z,fill_in,info,inlevs) + + use psb_base_mod + use psb_s_invk_fact_mod, psb_protect_name => psb_ssparse_invk + + integer(psb_ipk_), intent(in) :: n + type(psb_sspmat_type), intent(in) :: a + type(psb_sspmat_type), intent(inout) :: z + integer(psb_ipk_), intent(in) :: fill_in + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: inlevs(:) + ! + integer(psb_ipk_) :: i,j,k, err_act, nz, nzra, nzrz, ipz1,ipz2, nzz, ip1, ip2, l2 + integer(psb_ipk_), allocatable :: ia(:), ja(:), iz(:), jz(:) + real(psb_spk_), allocatable :: zw(:), val(:), valz(:) + integer(psb_ipk_), allocatable :: uplevs(:), rowlevs(:), idxs(:) + real(psb_spk_), allocatable :: row(:) + type(psb_s_coo_sparse_mat) :: trw + type(psb_s_csr_sparse_mat) :: acsr, zcsr + integer(psb_ipk_) :: ktrw, nidx + type(psb_i_heap) :: heap + + real(psb_dpk_) :: alpha + character(len=20) :: name='psb_sp_invk' + + info = psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + + if (.not.(a%is_triangle().and.a%is_unit().and.a%is_upper())) then + write(psb_err_unit,*) 'Wrong A ' + info = psb_err_internal_error_ + call psb_errpush(psb_err_internal_error_,name,a_err='wrong A') + goto 9999 + end if + call a%cp_to(acsr) + call trw%allocate(izero,izero,ione) + if (info == psb_success_) allocate(zw(n),iz(n),valz(n),& + & row(n),rowlevs(n),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + allocate(uplevs(acsr%get_nzeros()),stat=info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='Allocate') + goto 9999 + end if + uplevs(:) = 0 + row(:) = szero + rowlevs(:) = -(n+1) + + call zcsr%allocate(n,n,n*fill_in) + call zcsr%set_triangle() + call zcsr%set_unit(.false.) + call zcsr%set_upper() + call psb_ensure_size(n+1, idxs, info) + + + ! + ! + zcsr%irp(1) = 1 + nzz = 0 + + l2 = 0 + outer: do i = 1, n-1 + ! ZW = e_i + call psb_invk_copyin(i,n,acsr,ione,n,row,rowlevs,heap,ktrw,trw,info,& + & sign=-sone,inlevs=inlevs) + row(i) = sone + rowlevs(i) = 0 + + ! Update loop + call psb_invk_inv(fill_in,i,row,rowlevs,heap,& + & acsr%ja,acsr%irp,acsr%val,uplevs,nidx,idxs,info) + + call psb_invk_copyout(fill_in,i,n,row,rowlevs,nidx,idxs,& + & l2,zcsr%ja,zcsr%irp,zcsr%val,info) + + nzz = l2 + end do outer + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='mainloop') + goto 9999 + end if + ipz1 = nzz+1 + call psb_ensure_size(ipz1,zcsr%val,info) + call psb_ensure_size(ipz1,zcsr%ja,info) + zcsr%val(ipz1) = sone + zcsr%ja(ipz1) = n + zcsr%irp(n+1) = ipz1+1 + call zcsr%set_sorted() + call z%mv_from(zcsr) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_ssparse_invk + +subroutine psb_s_invk_copyin(i,m,a,jmin,jmax,row,rowlevs,heap,ktrw,trw,info,sign,inlevs) + + use psb_base_mod + use psb_s_invk_fact_mod, psb_protect_name => psb_s_invk_copyin + + implicit none + + ! Arguments + type(psb_s_csr_sparse_mat), intent(in) :: a + type(psb_s_coo_sparse_mat), intent(inout) :: trw + integer(psb_ipk_), intent(in) :: i,m,jmin,jmax + integer(psb_ipk_), intent(inout) :: ktrw,info + integer(psb_ipk_), intent(inout) :: rowlevs(:) + real(psb_spk_), intent(inout) :: row(:) + type(psb_i_heap), intent(inout) :: heap + real(psb_spk_), optional, intent(in) :: sign + integer(psb_ipk_), intent(in), optional :: inlevs(:) + + ! Local variables + integer(psb_ipk_) :: k,j,irb,err_act, nz + integer(psb_ipk_), parameter :: nrb=16 + real(psb_spk_) :: sign_ + character(len=20), parameter :: name='invk_copyin' + character(len=20) :: ch_err + + info = psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + + call heap%init(info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_init_heap') + goto 9999 + end if + + if (present(sign)) then + sign_ = sign + else + sign_ = sone + end if + + + ! + ! Take a fast shortcut if the matrix is stored in CSR format + ! + if (present(inlevs)) then + do j = a%irp(i), a%irp(i+1) - 1 + k = a%ja(j) + if ((jmin<=k).and.(k<=jmax)) then + row(k) = sign_ * a%val(j) + rowlevs(k) = inlevs(j) + call heap%insert(k,info) + end if + end do + else + do j = a%irp(i), a%irp(i+1) - 1 + k = a%ja(j) + if ((jmin<=k).and.(k<=jmax)) then + row(k) = sign_ * a%val(j) + rowlevs(k) = 0 + call heap%insert(k,info) + end if + end do + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_s_invk_copyin + + +subroutine psb_s_invk_copyout(fill_in,i,m,row,rowlevs,nidx,idxs,& + & l2,uia1,uia2,uaspk,info) + + use psb_base_mod + use psb_s_invk_fact_mod, psb_protect_name => psb_s_invk_copyout + + implicit none + + ! Arguments + integer(psb_ipk_), intent(in) :: fill_in, i, m, nidx + integer(psb_ipk_), intent(inout) :: l2, info + integer(psb_ipk_), intent(inout) :: rowlevs(:), idxs(:) + integer(psb_ipk_), allocatable, intent(inout) :: uia1(:), uia2(:) + real(psb_spk_), allocatable, intent(inout) :: uaspk(:) + real(psb_spk_), intent(inout) :: row(:) + + ! Local variables + integer(psb_ipk_) :: j,isz,err_act,int_err(5),idxp + character(len=20), parameter :: name='psb_siluk_factint' + character(len=20) :: ch_err + + info = psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + + do idxp=1,nidx + + j = idxs(idxp) + + if (j>=i) then + ! + ! Copy the upper part of the row + ! + if (rowlevs(j) <= fill_in) then + l2 = l2 + 1 + if (size(uaspk) < l2) then + ! + ! Figure out a good reallocation size! + ! + isz = max(int(1.2*l2),l2+100) + call psb_realloc(isz,uaspk,info) + if (info == psb_success_) call psb_realloc(isz,uia1,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='Allocate') + goto 9999 + end if + end if + uia1(l2) = j + uaspk(l2) = row(j) + end if + ! + ! Re-initialize row(j) and rowlevs(j) + ! + row(j) = szero + rowlevs(j) = -(m+1) + end if + end do + + uia2(i+1) = l2 + 1 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_s_invk_copyout + +subroutine psb_sinvk_inv(fill_in,i,row,rowlevs,heap,ja,irp,val,uplevs,& + & nidx,idxs,info) + + use psb_base_mod + use psb_s_invk_fact_mod, psb_protect_name => psb_sinvk_inv + + implicit none + + ! Arguments + type(psb_i_heap), intent(inout) :: heap + integer(psb_ipk_), intent(in) :: i, fill_in + integer(psb_ipk_), intent(inout) :: nidx,info + integer(psb_ipk_), intent(inout) :: rowlevs(:) + integer(psb_ipk_), allocatable, intent(inout) :: idxs(:) + integer(psb_ipk_), intent(in) :: ja(:),irp(:),uplevs(:) + real(psb_spk_), intent(in) :: val(:) + real(psb_spk_), intent(inout) :: row(:) + + ! Local variables + integer(psb_ipk_) :: k,j,lrwk,jj,lastk, iret + real(psb_dpk_) :: rwk + + + info = psb_success_ + + call psb_ensure_size(200, idxs, info) + if (info /= psb_success_) return + nidx = 1 + idxs(1) = i + lastk = i + + ! + ! Do while there are indices to be processed + ! + do + ! Beware: (iret < 0) means that the heap is empty, not an error. + call heap%get_first(k,iret) + if (iret < 0) then +!!$ write(psb_err_unit,*) 'IINVK: ',i,' returning at ',lastk + return + end if + + ! + ! Just in case an index has been put on the heap more than once. + ! + if (k == lastk) cycle + + lastk = k + nidx = nidx + 1 + if (nidx>size(idxs)) then + call psb_realloc(nidx+psb_heap_resize,idxs,info) + if (info /= psb_success_) return + end if + idxs(nidx) = k + + if ((row(k) /= szero).and.(rowlevs(k) <= fill_in)) then + ! + ! Note: since U is scaled while copying it out (see iluk_copyout), + ! we can use rwk in the update below + ! + rwk = row(k) + lrwk = rowlevs(k) + + do jj=irp(k),irp(k+1)-1 + j = ja(jj) + if (j<=k) then + info = -i + return + endif + ! + ! Insert the index into the heap for further processing. + ! The fill levels are initialized to a negative value. If we find + ! one, it means that it is an as yet untouched index, so we need + ! to insert it; otherwise it is already on the heap, there is no + ! need to insert it more than once. + ! + if (rowlevs(j)<0) then + call heap%insert(j,info) + if (info /= psb_success_) return + rowlevs(j) = abs(rowlevs(j)) + end if + ! + ! Update row(j) and the corresponding fill level + ! + row(j) = row(j) - rwk * val(jj) + rowlevs(j) = min(rowlevs(j),lrwk+uplevs(jj)+1) + end do + + end if + end do + +end subroutine psb_sinvk_inv diff --git a/prec/impl/psb_s_invt_fact.f90 b/prec/impl/psb_s_invt_fact.f90 new file mode 100644 index 00000000..3a261ca6 --- /dev/null +++ b/prec/impl/psb_s_invt_fact.f90 @@ -0,0 +1,631 @@ +! +! 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 AMG4PSBLAS, 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_invt_bld(a,fillin,invfill,thresh,invthresh,& + & lmat,d,umat,desc,info,blck) + + use psb_base_mod + use psb_s_invt_fact_mod, psb_protect_name => psb_s_invt_bld + use psb_s_ilu_fact_mod + implicit none + + ! Arguments + type(psb_sspmat_type), intent(in), target :: a + integer(psb_ipk_), intent(in) :: fillin,invfill + real(psb_spk_), intent(in) :: thresh + real(psb_spk_), intent(in) :: invthresh + type(psb_sspmat_type), intent(inout) :: lmat, umat + real(psb_spk_), allocatable :: d(:) + Type(psb_desc_type), Intent(inout) :: desc + integer(psb_ipk_), intent(out) :: info + type(psb_sspmat_type), intent(in), optional :: blck + ! + integer(psb_ipk_) :: i, nztota, err_act, n_row, nrow_a, n_col + type(psb_sspmat_type) :: atmp + real(psb_spk_), allocatable :: pq(:), pd(:), w(:) + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: nzrmax + real(psb_spk_) :: sp_thresh + character(len=20) :: name, ch_err, fname + + + info = psb_success_ + name='psb_sinvt_fact' + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = psb_cd_get_context(desc) + call psb_info(ictxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + ! + ! Check the memory available to hold the incomplete L and U factors + ! and allocate it if needed + ! + nrow_a = a%get_nrows() + nztota = a%get_nzeros() + + if (present(blck)) then + nztota = nztota + blck%get_nzeros() + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & ': out get_nnzeros',nrow_a,nztota,& + & a%get_nrows(),a%get_ncols(),a%get_nzeros() + + + n_row = psb_cd_get_local_rows(desc) + n_col = psb_cd_get_local_cols(desc) + allocate(pd(n_row),w(n_row),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + nzrmax = fillin + sp_thresh = thresh + + call lmat%allocate(n_row,n_row,info,nz=nztota) + if (info == psb_success_) call umat%allocate(n_row,n_row,info,nz=nztota) + + if (info == 0) call psb_ilut_fact(nzrmax,sp_thresh,& + & a,lmat,umat,pd,info,blck=blck,iscale=psb_ilu_scale_maxval_) + + if (info == psb_success_) call atmp%allocate(n_row,n_row,info,nz=nztota) + if(info/=0) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if (.false.) then +!!$ if (debug_level >= psb_debug_inner_) then + write(fname,'(a,i0,a)') 'invt-lo-',me,'.mtx' + call lmat%print(fname,head="INVTLOW") + write(fname,'(a,i0,a)') 'invt-up-',me,'.mtx' + call umat%print(fname,head="INVTUPP") + end if + + ! + ! Compute the approx U^-1 and L^-1 + ! + nzrmax = invfill + call psb_ssparse_invt(n_row,umat,atmp,nzrmax,invthresh,info) + if (info == psb_success_) call psb_move_alloc(atmp,umat,info) + if (info == psb_success_) call lmat%transp() + if (info == psb_success_) call psb_ssparse_invt(n_row,lmat,atmp,nzrmax,invthresh,info) + if (info == psb_success_) call psb_move_alloc(atmp,lmat,info) + if (info == psb_success_) call lmat%transp() + ! Done. Hopefully.... + + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='invt') + goto 9999 + end if + + call psb_move_alloc(pd,d,info) + call lmat%set_asb() + call lmat%trim() + call umat%set_asb() + call umat%trim() + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_s_invt_bld + +subroutine psb_ssparse_invt(n,a,z,nzrmax,sp_thresh,info) + + use psb_base_mod + use psb_s_invt_fact_mod, psb_protect_name => psb_ssparse_invt + + implicit none + integer(psb_ipk_), intent(in) :: n + type(psb_sspmat_type), intent(in) :: a + type(psb_sspmat_type), intent(inout) :: z + integer(psb_ipk_), intent(in) :: nzrmax + real(psb_spk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(out) :: info + ! + integer(psb_ipk_) :: i,j,k, err_act, nz, nzra, nzrz, ipz1,ipz2, nzz, ip1, ip2, l2 + integer(psb_ipk_), allocatable :: ia(:), ja(:), iz(:),jz(:) + real(psb_spk_), allocatable :: zw(:), val(:), valz(:) + integer(psb_ipk_), allocatable :: uplevs(:), rowlevs(:),idxs(:) + real(psb_spk_), allocatable :: row(:) + type(psb_s_coo_sparse_mat) :: trw + type(psb_s_csr_sparse_mat) :: acsr, zcsr + integer(psb_ipk_) :: ktrw, nidx, nlw,nup,jmaxup + type(psb_i_heap) :: heap + real(psb_spk_) :: alpha, nrmi + character(len=20) :: name='psb_sp_invt' + + info = psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + + if (.not.(a%is_triangle().and.a%is_unit().and.a%is_upper())) then + write(psb_err_unit,*) 'Wrong A ' + info = psb_err_internal_error_ + call psb_errpush(psb_err_internal_error_,name,a_err='wrong A') + goto 9999 + end if + call a%cp_to(acsr) + call trw%allocate(izero,izero,ione) + if (info == psb_success_) allocate(zw(n),iz(n),valz(n),& + & row(n),rowlevs(n),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + call zcsr%allocate(n,n,n*nzrmax) + call zcsr%set_triangle() + call zcsr%set_unit(.false.) + call zcsr%set_upper() + ! + ! + nzz = 0 + row(:) = szero + rowlevs(:) = 0 + l2 = 0 + zcsr%irp(1) = 1 + + outer: do i = 1, n-1 + ! ZW = e_i + call psb_s_invt_copyin(i,n,acsr,i,ione,n,nlw,nup,jmaxup,nrmi,row,& + & heap,rowlevs,ktrw,trw,info,sign=-sone) + if (info /= 0) exit + row(i) = sone + ! Adjust norm + if (nrmi < sone) then + nrmi = sqrt(sone + nrmi**2) + else + nrmi = nrmi*sqrt(sone+sone/(nrmi**2)) + end if + + call psb_invt_inv(sp_thresh,i,nrmi,row,heap,rowlevs,& + & acsr%ja,acsr%irp,acsr%val,nidx,idxs,info) + if (info /= 0) exit +!!$ write(0,*) 'Calling copyout ',nzrmax,nlw,nup,nidx,l2 + call psb_s_invt_copyout(nzrmax,sp_thresh,i,n,nlw,nup,jmaxup,nrmi,row,& + & nidx,idxs,l2,zcsr%ja,zcsr%irp,zcsr%val,info) + if (info /= 0) exit + nzz = l2 + end do outer + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='mainloop') + goto 9999 + end if + + ipz1 = nzz+1 + call psb_ensure_size(ipz1,zcsr%val,info) + call psb_ensure_size(ipz1,zcsr%ja,info) + zcsr%val(ipz1) = sone + zcsr%ja(ipz1) = n + zcsr%irp(n+1) = ipz1+1 + + call z%mv_from(zcsr) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_ssparse_invt + +subroutine psb_s_invt_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,nrmi,row,heap,& + & irwt,ktrw,trw,info,sign) + use psb_base_mod + use psb_d_invt_fact_mod, psb_protect_name => psb_d_invt_copyin + implicit none + type(psb_s_csr_sparse_mat), intent(in) :: a + type(psb_s_coo_sparse_mat), intent(inout) :: trw + integer(psb_ipk_), intent(in) :: i, m,jmin,jmax,jd + integer(psb_ipk_), intent(inout) :: ktrw,nlw,nup,jmaxup,info + integer(psb_ipk_), intent(inout) :: irwt(:) + real(psb_spk_), intent(inout) :: nrmi + real(psb_spk_), intent(inout) :: row(:) + type(psb_i_heap), intent(inout) :: heap + real(psb_spk_), intent(in), optional :: sign + ! + integer(psb_ipk_) :: k,j,irb,kin,nz, err_act + integer(psb_ipk_), parameter :: nrb=16 + real(psb_dpk_) :: dmaxup, sign_ + real(psb_dpk_), external :: dnrm2 + character(len=20), parameter :: name='invt_copyin' + + info = psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + + call heap%init(info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_init_heap') + goto 9999 + end if + sign_ = sone + if (present(sign)) sign_ = sign + ! + ! nrmi is the norm of the current sparse row (for the time being, + ! we use the 2-norm). + ! NOTE: the 2-norm below includes also elements that are outside + ! [jmin:jmax] strictly. Is this really important? TO BE CHECKED. + ! + + nlw = 0 + nup = 0 + jmaxup = 0 + dmaxup = szero + nrmi = szero + + do j = a%irp(i), a%irp(i+1) - 1 + k = a%ja(j) + if ((jmin<=k).and.(k<=jmax)) then + row(k) = sign_ * a%val(j) + call heap%insert(k,info) + irwt(k) = 1 + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_insert_heap') + goto 9999 + end if + end if + if (kjd) then + nup = nup + 1 + if (abs(row(k))>dmaxup) then + jmaxup = k + dmaxup = abs(row(k)) + end if + end if + end do + nz = a%irp(i+1) - a%irp(i) + nrmi = dnrm2(nz,a%val(a%irp(i):),ione) + + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_s_invt_copyin + +subroutine psb_s_invt_copyout(fill_in,thres,i,m,nlw,nup,jmaxup,nrmi,row, & + & nidx,idxs,l2,ja,irp,val,info) + + use psb_base_mod + use psb_s_invt_fact_mod, psb_protect_name => psb_s_invt_copyout + + implicit none + + ! Arguments + integer(psb_ipk_), intent(in) :: fill_in,i,m,nidx,nlw,nup,jmaxup + integer(psb_ipk_), intent(in) :: idxs(:) + integer(psb_ipk_), intent(inout) :: l2, info + integer(psb_ipk_), allocatable, intent(inout) :: ja(:),irp(:) + real(psb_spk_), intent(in) :: thres,nrmi + real(psb_spk_),allocatable, intent(inout) :: val(:) + real(psb_spk_), intent(inout) :: row(:) + + ! Local variables + real(psb_dpk_),allocatable :: xw(:) + integer(psb_ipk_), allocatable :: xwid(:), indx(:) + real(psb_dpk_) :: witem, wmin + integer(psb_ipk_) :: widx + integer(psb_ipk_) :: k,isz,err_act,int_err(5),idxp, nz + type(psb_d_idx_heap) :: heap + character(len=20), parameter :: name='invt_copyout' + character(len=20) :: ch_err + logical :: fndmaxup + + info = psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + + ! + ! Here we need to apply also the dropping rule base on the fill-in. + ! We do it by putting into a heap the elements that are not dropped + ! by using the 2-norm rule, and then copying them out. + ! + ! The heap goes down on the entry absolute value, so the first item + ! is the largest absolute value. + ! +!!$ write(0,*) 'invt_copyout ',nidx,nup+fill_in + call heap%init(info,dir=psb_asort_down_) + + if (info == psb_success_) allocate(xwid(nidx),xw(nidx),indx(nidx),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/3*nidx/),& + & a_err='real(psb_dpk_)') + goto 9999 + end if + + ! + ! First the lower part + ! + + nz = 0 + idxp = 0 + + do + + idxp = idxp + 1 + if (idxp > nidx) exit + if (idxs(idxp) >= i) exit + widx = idxs(idxp) + witem = row(widx) + ! + ! Dropping rule based on the 2-norm + ! + if (abs(witem) < thres*nrmi) cycle + + nz = nz + 1 + xw(nz) = witem + xwid(nz) = widx + call heap%insert(witem,widx,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_insert_heap') + goto 9999 + end if + end do + + if (nz > 1) then + write(psb_err_unit,*) 'Warning: lower triangle from invt???? ' + end if + + + if (idxp <= size(idxs)) then + if (idxs(idxp) < i) then + do + idxp = idxp + 1 + if (idxp > nidx) exit + if (idxs(idxp) >= i) exit + end do + end if + end if + idxp = idxp - 1 + nz = 0 + wmin=HUGE(wmin) + if (.false.) then + do + + idxp = idxp + 1 + if (idxp > nidx) exit + widx = idxs(idxp) + if (widx < i) then + write(psb_err_unit,*) 'Warning: lower triangle in upper copy',widx,i,idxp,idxs(idxp) + cycle + end if + if (widx > m) then + cycle + end if + witem = row(widx) + ! + ! Dropping rule based on the 2-norm. But keep the jmaxup-th entry anyway. + ! + if ((widx /= jmaxup) .and. (widx /= i) .and. (abs(witem) < thres*nrmi)) then + cycle + end if + if ((widx/=jmaxup).and.(nz > nup+fill_in)) then + if (abs(witem) < wmin) cycle + endif + wmin = min(abs(witem),wmin) + nz = nz + 1 + xw(nz) = witem + xwid(nz) = widx + call heap%insert(witem,widx,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_insert_heap') + goto 9999 + end if + + end do + + ! + ! Now we have to take out the first nup-fill_in entries. But make sure + ! we include entry jmaxup. + ! + if (nz <= nup+fill_in) then + ! + ! Just copy everything from xw + ! + fndmaxup=.true. + else + fndmaxup = .false. + nz = nup+fill_in + do k=1,nz + call heap%get_first(witem,widx,info) + xw(k) = witem + xwid(k) = widx + if (widx == jmaxup) fndmaxup=.true. + end do + end if + if ((i nidx) exit + widx = idxs(idxp) + if (widx < i) then + write(psb_err_unit,*) 'Warning: lower triangle in upper copy',widx,i,idxp,idxs(idxp) + cycle + end if + if (widx > m) then + cycle + end if + witem = row(widx) + ! + ! Dropping rule based on the 2-norm. But keep the jmaxup-th entry anyway. + ! + if ((widx /= i) .and. (abs(witem) < thres*nrmi)) then + cycle + end if + if (nz > nup+fill_in) then + if (abs(witem) < wmin) cycle + endif + wmin = min(abs(witem),wmin) + nz = nz + 1 + xw(nz) = witem + xwid(nz) = widx + call heap%insert(witem,widx,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_insert_heap') + goto 9999 + end if + + end do + + ! + ! Now we have to take out the first nup-fill_in entries. But make sure + ! we include entry jmaxup. + ! + if (nz > nup+fill_in) then + nz = nup+fill_in + do k=1,nz + call heap%get_first(witem,widx,info) + xw(k) = witem + xwid(k) = widx + end do + end if + end if + + ! + ! Now we put things back into ascending column order + ! + call psb_msort(xwid(1:nz),indx(1:nz),dir=psb_sort_up_) + + ! + ! Copy out the upper part of the row + ! + do k=1,nz + l2 = l2 + 1 + if (size(val) < l2) then + ! + ! Figure out a good reallocation size! + ! + isz = max(int(1.2*l2),l2+100) + call psb_realloc(isz,val,info) + if (info == psb_success_) call psb_realloc(isz,ja,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='Allocate') + goto 9999 + end if + end if + ja(l2) = xwid(k) + val(l2) = xw(indx(k)) + end do + + ! + ! Set row to zero + ! + do idxp=1,nidx + row(idxs(idxp)) = szero + end do + + irp(i+1) = l2 + 1 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_s_invt_copyout diff --git a/prec/impl/psb_z_invk_fact.f90 b/prec/impl/psb_z_invk_fact.f90 new file mode 100644 index 00000000..4f109701 --- /dev/null +++ b/prec/impl/psb_z_invk_fact.f90 @@ -0,0 +1,494 @@ +! +! +! 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_invk_bld(a,fill1, fill2,lmat,d,umat,desc,info,blck) + + use psb_base_mod + use psb_z_invk_fact_mod, psb_protect_name => psb_z_invk_bld + use psb_z_ilu_fact_mod + implicit none + + ! Arguments + type(psb_zspmat_type), intent(in), target :: a + integer(psb_ipk_), intent(in) :: fill1, fill2 + type(psb_zspmat_type), intent(inout) :: lmat, umat + complex(psb_dpk_), allocatable :: d(:) + Type(psb_desc_type), Intent(inout) :: desc + integer(psb_ipk_), intent(out) :: info + type(psb_zspmat_type), intent(in), optional :: blck + ! + integer(psb_ipk_) :: i, nztota, err_act, n_row, nrow_a, n_col + type(psb_zspmat_type) :: atmp + complex(psb_dpk_), allocatable :: pq(:), pd(:) + integer(psb_ipk_), allocatable :: uplevs(:) + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: nzrmax + character(len=20) :: name, ch_err + + + info = psb_success_ + name='psb_zinvk_bld' + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = psb_cd_get_context(desc) + call psb_info(ictxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + ! + ! Check the memory available to hold the incomplete L and U factors + ! and allocate it if needed + ! + nrow_a = a%get_nrows() + nztota = a%get_nzeros() + + if (present(blck)) then + nztota = nztota + blck%get_nzeros() + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & ': out get_nnzeros',nrow_a,nztota,& + & a%get_nrows(),a%get_ncols(),a%get_nzeros() + + + n_row = psb_cd_get_local_rows(desc) + n_col = psb_cd_get_local_cols(desc) + allocate(pd(n_row),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + + call lmat%allocate(n_row,n_row,info,nz=nztota) + if (info == psb_success_) call umat%allocate(n_row,n_row,info,nz=nztota) + + + call psb_iluk_fact(fill1,psb_ilu_n_,a,lmat,umat,pd,info,blck=blck) + !,uplevs=uplevs) + !call psb_ziluk_fact(fill1,psb_ilu_n_,a,lmat,umat,pd,info,blck=blck) + + if (info == psb_success_) call atmp%allocate(n_row,n_row,info,nz=nztota) + if(info/=0) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + + ! + ! Compute the aprox U^-1 and L^-1 + ! + call psb_sparse_invk(n_row,umat,atmp,fill2,info) + if (info == psb_success_) call psb_move_alloc(atmp,umat,info) + if (info == psb_success_) call lmat%transp() + if (info == psb_success_) call psb_sparse_invk(n_row,lmat,atmp,fill2,info) + if (info == psb_success_) call psb_move_alloc(atmp,lmat,info) + if (info == psb_success_) call lmat%transp() + ! Done. Hopefully.... + + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='invt') + goto 9999 + end if + + call psb_move_alloc(pd,d,info) + call lmat%set_asb() + call lmat%trim() + call umat%set_asb() + call umat%trim() + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_z_invk_bld + +subroutine psb_zsparse_invk(n,a,z,fill_in,info,inlevs) + + use psb_base_mod + use psb_z_invk_fact_mod, psb_protect_name => psb_zsparse_invk + + integer(psb_ipk_), intent(in) :: n + type(psb_zspmat_type), intent(in) :: a + type(psb_zspmat_type), intent(inout) :: z + integer(psb_ipk_), intent(in) :: fill_in + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: inlevs(:) + ! + integer(psb_ipk_) :: i,j,k, err_act, nz, nzra, nzrz, ipz1,ipz2, nzz, ip1, ip2, l2 + integer(psb_ipk_), allocatable :: ia(:), ja(:), iz(:), jz(:) + complex(psb_dpk_), allocatable :: zw(:), val(:), valz(:) + integer(psb_ipk_), allocatable :: uplevs(:), rowlevs(:), idxs(:) + complex(psb_dpk_), allocatable :: row(:) + type(psb_z_coo_sparse_mat) :: trw + type(psb_z_csr_sparse_mat) :: acsr, zcsr + integer(psb_ipk_) :: ktrw, nidx + type(psb_i_heap) :: heap + + real(psb_dpk_) :: alpha + character(len=20) :: name='psb_sp_invk' + + info = psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + + if (.not.(a%is_triangle().and.a%is_unit().and.a%is_upper())) then + write(psb_err_unit,*) 'Wrong A ' + info = psb_err_internal_error_ + call psb_errpush(psb_err_internal_error_,name,a_err='wrong A') + goto 9999 + end if + call a%cp_to(acsr) + call trw%allocate(izero,izero,ione) + if (info == psb_success_) allocate(zw(n),iz(n),valz(n),& + & row(n),rowlevs(n),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + allocate(uplevs(acsr%get_nzeros()),stat=info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='Allocate') + goto 9999 + end if + uplevs(:) = 0 + row(:) = zzero + rowlevs(:) = -(n+1) + + call zcsr%allocate(n,n,n*fill_in) + call zcsr%set_triangle() + call zcsr%set_unit(.false.) + call zcsr%set_upper() + call psb_ensure_size(n+1, idxs, info) + + + ! + ! + zcsr%irp(1) = 1 + nzz = 0 + + l2 = 0 + outer: do i = 1, n-1 + ! ZW = e_i + call psb_invk_copyin(i,n,acsr,ione,n,row,rowlevs,heap,ktrw,trw,info,& + & sign=-done,inlevs=inlevs) + row(i) = zone + rowlevs(i) = 0 + + ! Update loop + call psb_invk_inv(fill_in,i,row,rowlevs,heap,& + & acsr%ja,acsr%irp,acsr%val,uplevs,nidx,idxs,info) + + call psb_invk_copyout(fill_in,i,n,row,rowlevs,nidx,idxs,& + & l2,zcsr%ja,zcsr%irp,zcsr%val,info) + + nzz = l2 + end do outer + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='mainloop') + goto 9999 + end if + ipz1 = nzz+1 + call psb_ensure_size(ipz1,zcsr%val,info) + call psb_ensure_size(ipz1,zcsr%ja,info) + zcsr%val(ipz1) = zone + zcsr%ja(ipz1) = n + zcsr%irp(n+1) = ipz1+1 + call zcsr%set_sorted() + call z%mv_from(zcsr) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_zsparse_invk + +subroutine psb_z_invk_copyin(i,m,a,jmin,jmax,row,rowlevs,heap,ktrw,trw,info,sign,inlevs) + + use psb_base_mod + use psb_z_invk_fact_mod, psb_protect_name => psb_z_invk_copyin + + implicit none + + ! Arguments + type(psb_z_csr_sparse_mat), intent(in) :: a + type(psb_z_coo_sparse_mat), intent(inout) :: trw + integer(psb_ipk_), intent(in) :: i,m,jmin,jmax + integer(psb_ipk_), intent(inout) :: ktrw,info + integer(psb_ipk_), intent(inout) :: rowlevs(:) + complex(psb_dpk_), intent(inout) :: row(:) + type(psb_i_heap), intent(inout) :: heap + real(psb_dpk_), optional, intent(in) :: sign + integer(psb_ipk_), intent(in), optional :: inlevs(:) + + ! Local variables + integer(psb_ipk_) :: k,j,irb,err_act, nz + integer(psb_ipk_), parameter :: nrb=16 + real(psb_dpk_) :: sign_ + character(len=20), parameter :: name='invk_copyin' + character(len=20) :: ch_err + + info = psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + + call heap%init(info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_init_heap') + goto 9999 + end if + + if (present(sign)) then + sign_ = sign + else + sign_ = done + end if + + + ! + ! Take a fast shortcut if the matrix is stored in CSR format + ! + if (present(inlevs)) then + do j = a%irp(i), a%irp(i+1) - 1 + k = a%ja(j) + if ((jmin<=k).and.(k<=jmax)) then + row(k) = sign_ * a%val(j) + rowlevs(k) = inlevs(j) + call heap%insert(k,info) + end if + end do + else + do j = a%irp(i), a%irp(i+1) - 1 + k = a%ja(j) + if ((jmin<=k).and.(k<=jmax)) then + row(k) = sign_ * a%val(j) + rowlevs(k) = 0 + call heap%insert(k,info) + end if + end do + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_z_invk_copyin + + +subroutine psb_z_invk_copyout(fill_in,i,m,row,rowlevs,nidx,idxs,& + & l2,uia1,uia2,uaspk,info) + + use psb_base_mod + use psb_z_invk_fact_mod, psb_protect_name => psb_z_invk_copyout + + implicit none + + ! Arguments + integer(psb_ipk_), intent(in) :: fill_in, i, m, nidx + integer(psb_ipk_), intent(inout) :: l2, info + integer(psb_ipk_), intent(inout) :: rowlevs(:), idxs(:) + integer(psb_ipk_), allocatable, intent(inout) :: uia1(:), uia2(:) + complex(psb_dpk_), allocatable, intent(inout) :: uaspk(:) + complex(psb_dpk_), intent(inout) :: row(:) + + ! Local variables + integer(psb_ipk_) :: j,isz,err_act,int_err(5),idxp + character(len=20), parameter :: name='psb_ziluk_factint' + character(len=20) :: ch_err + + info = psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + + do idxp=1,nidx + + j = idxs(idxp) + + if (j>=i) then + ! + ! Copy the upper part of the row + ! + if (rowlevs(j) <= fill_in) then + l2 = l2 + 1 + if (size(uaspk) < l2) then + ! + ! Figure out a good reallocation size! + ! + isz = max(int(1.2*l2),l2+100) + call psb_realloc(isz,uaspk,info) + if (info == psb_success_) call psb_realloc(isz,uia1,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='Allocate') + goto 9999 + end if + end if + uia1(l2) = j + uaspk(l2) = row(j) + end if + ! + ! Re-initialize row(j) and rowlevs(j) + ! + row(j) = zzero + rowlevs(j) = -(m+1) + end if + end do + + uia2(i+1) = l2 + 1 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_z_invk_copyout + +subroutine psb_zinvk_inv(fill_in,i,row,rowlevs,heap,ja,irp,val,uplevs,& + & nidx,idxs,info) + + use psb_base_mod + use psb_z_invk_fact_mod, psb_protect_name => psb_zinvk_inv + + implicit none + + ! Arguments + type(psb_i_heap), intent(inout) :: heap + integer(psb_ipk_), intent(in) :: i, fill_in + integer(psb_ipk_), intent(inout) :: nidx,info + integer(psb_ipk_), intent(inout) :: rowlevs(:) + integer(psb_ipk_), allocatable, intent(inout) :: idxs(:) + integer(psb_ipk_), intent(in) :: ja(:),irp(:),uplevs(:) + complex(psb_dpk_), intent(in) :: val(:) + complex(psb_dpk_), intent(inout) :: row(:) + + ! Local variables + integer(psb_ipk_) :: k,j,lrwk,jj,lastk, iret + real(psb_dpk_) :: rwk + + + info = psb_success_ + + call psb_ensure_size(200, idxs, info) + if (info /= psb_success_) return + nidx = 1 + idxs(1) = i + lastk = i + + ! + ! Do while there are indices to be processed + ! + do + ! Beware: (iret < 0) means that the heap is empty, not an error. + call heap%get_first(k,iret) + if (iret < 0) then +!!$ write(psb_err_unit,*) 'IINVK: ',i,' returning at ',lastk + return + end if + + ! + ! Just in case an index has been put on the heap more than once. + ! + if (k == lastk) cycle + + lastk = k + nidx = nidx + 1 + if (nidx>size(idxs)) then + call psb_realloc(nidx+psb_heap_resize,idxs,info) + if (info /= psb_success_) return + end if + idxs(nidx) = k + + if ((row(k) /= zzero).and.(rowlevs(k) <= fill_in)) then + ! + ! Note: since U is scaled while copying it out (see iluk_copyout), + ! we can use rwk in the update below + ! + rwk = row(k) + lrwk = rowlevs(k) + + do jj=irp(k),irp(k+1)-1 + j = ja(jj) + if (j<=k) then + info = -i + return + endif + ! + ! Insert the index into the heap for further processing. + ! The fill levels are initialized to a negative value. If we find + ! one, it means that it is an as yet untouched index, so we need + ! to insert it; otherwise it is already on the heap, there is no + ! need to insert it more than once. + ! + if (rowlevs(j)<0) then + call heap%insert(j,info) + if (info /= psb_success_) return + rowlevs(j) = abs(rowlevs(j)) + end if + ! + ! Update row(j) and the corresponding fill level + ! + row(j) = row(j) - rwk * val(jj) + rowlevs(j) = min(rowlevs(j),lrwk+uplevs(jj)+1) + end do + + end if + end do + +end subroutine psb_zinvk_inv diff --git a/prec/impl/psb_z_invt_fact.f90 b/prec/impl/psb_z_invt_fact.f90 new file mode 100644 index 00000000..17e38c67 --- /dev/null +++ b/prec/impl/psb_z_invt_fact.f90 @@ -0,0 +1,631 @@ +! +! 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 AMG4PSBLAS, 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_invt_bld(a,fillin,invfill,thresh,invthresh,& + & lmat,d,umat,desc,info,blck) + + use psb_base_mod + use psb_z_invt_fact_mod, psb_protect_name => psb_z_invt_bld + use psb_z_ilu_fact_mod + implicit none + + ! Arguments + type(psb_zspmat_type), intent(in), target :: a + integer(psb_ipk_), intent(in) :: fillin,invfill + real(psb_dpk_), intent(in) :: thresh + real(psb_dpk_), intent(in) :: invthresh + type(psb_zspmat_type), intent(inout) :: lmat, umat + complex(psb_dpk_), allocatable :: d(:) + Type(psb_desc_type), Intent(inout) :: desc + integer(psb_ipk_), intent(out) :: info + type(psb_zspmat_type), intent(in), optional :: blck + ! + integer(psb_ipk_) :: i, nztota, err_act, n_row, nrow_a, n_col + type(psb_zspmat_type) :: atmp + complex(psb_dpk_), allocatable :: pq(:), pd(:), w(:) + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: nzrmax + real(psb_dpk_) :: sp_thresh + character(len=20) :: name, ch_err, fname + + + info = psb_success_ + name='psb_zinvt_fact' + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = psb_cd_get_context(desc) + call psb_info(ictxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + ! + ! Check the memory available to hold the incomplete L and U factors + ! and allocate it if needed + ! + nrow_a = a%get_nrows() + nztota = a%get_nzeros() + + if (present(blck)) then + nztota = nztota + blck%get_nzeros() + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & ': out get_nnzeros',nrow_a,nztota,& + & a%get_nrows(),a%get_ncols(),a%get_nzeros() + + + n_row = psb_cd_get_local_rows(desc) + n_col = psb_cd_get_local_cols(desc) + allocate(pd(n_row),w(n_row),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + nzrmax = fillin + sp_thresh = thresh + + call lmat%allocate(n_row,n_row,info,nz=nztota) + if (info == psb_success_) call umat%allocate(n_row,n_row,info,nz=nztota) + + if (info == 0) call psb_ilut_fact(nzrmax,sp_thresh,& + & a,lmat,umat,pd,info,blck=blck,iscale=psb_ilu_scale_maxval_) + + if (info == psb_success_) call atmp%allocate(n_row,n_row,info,nz=nztota) + if(info/=0) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if (.false.) then +!!$ if (debug_level >= psb_debug_inner_) then + write(fname,'(a,i0,a)') 'invt-lo-',me,'.mtx' + call lmat%print(fname,head="INVTLOW") + write(fname,'(a,i0,a)') 'invt-up-',me,'.mtx' + call umat%print(fname,head="INVTUPP") + end if + + ! + ! Compute the approx U^-1 and L^-1 + ! + nzrmax = invfill + call psb_zsparse_invt(n_row,umat,atmp,nzrmax,invthresh,info) + if (info == psb_success_) call psb_move_alloc(atmp,umat,info) + if (info == psb_success_) call lmat%transp() + if (info == psb_success_) call psb_zsparse_invt(n_row,lmat,atmp,nzrmax,invthresh,info) + if (info == psb_success_) call psb_move_alloc(atmp,lmat,info) + if (info == psb_success_) call lmat%transp() + ! Done. Hopefully.... + + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='invt') + goto 9999 + end if + + call psb_move_alloc(pd,d,info) + call lmat%set_asb() + call lmat%trim() + call umat%set_asb() + call umat%trim() + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_z_invt_bld + +subroutine psb_zsparse_invt(n,a,z,nzrmax,sp_thresh,info) + + use psb_base_mod + use psb_z_invt_fact_mod, psb_protect_name => psb_zsparse_invt + + implicit none + integer(psb_ipk_), intent(in) :: n + type(psb_zspmat_type), intent(in) :: a + type(psb_zspmat_type), intent(inout) :: z + integer(psb_ipk_), intent(in) :: nzrmax + real(psb_dpk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(out) :: info + ! + integer(psb_ipk_) :: i,j,k, err_act, nz, nzra, nzrz, ipz1,ipz2, nzz, ip1, ip2, l2 + integer(psb_ipk_), allocatable :: ia(:), ja(:), iz(:),jz(:) + complex(psb_dpk_), allocatable :: zw(:), val(:), valz(:) + integer(psb_ipk_), allocatable :: uplevs(:), rowlevs(:),idxs(:) + complex(psb_dpk_), allocatable :: row(:) + type(psb_z_coo_sparse_mat) :: trw + type(psb_z_csr_sparse_mat) :: acsr, zcsr + integer(psb_ipk_) :: ktrw, nidx, nlw,nup,jmaxup + type(psb_i_heap) :: heap + real(psb_dpk_) :: alpha, nrmi + character(len=20) :: name='psb_sp_invt' + + info = psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + + if (.not.(a%is_triangle().and.a%is_unit().and.a%is_upper())) then + write(psb_err_unit,*) 'Wrong A ' + info = psb_err_internal_error_ + call psb_errpush(psb_err_internal_error_,name,a_err='wrong A') + goto 9999 + end if + call a%cp_to(acsr) + call trw%allocate(izero,izero,ione) + if (info == psb_success_) allocate(zw(n),iz(n),valz(n),& + & row(n),rowlevs(n),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + call zcsr%allocate(n,n,n*nzrmax) + call zcsr%set_triangle() + call zcsr%set_unit(.false.) + call zcsr%set_upper() + ! + ! + nzz = 0 + row(:) = zzero + rowlevs(:) = 0 + l2 = 0 + zcsr%irp(1) = 1 + + outer: do i = 1, n-1 + ! ZW = e_i + call psb_z_invt_copyin(i,n,acsr,i,ione,n,nlw,nup,jmaxup,nrmi,row,& + & heap,rowlevs,ktrw,trw,info,sign=-done) + if (info /= 0) exit + row(i) = zone + ! Adjust norm + if (nrmi < done) then + nrmi = sqrt(done + nrmi**2) + else + nrmi = nrmi*sqrt(zone+zone/(nrmi**2)) + end if + + call psb_invt_inv(sp_thresh,i,nrmi,row,heap,rowlevs,& + & acsr%ja,acsr%irp,acsr%val,nidx,idxs,info) + if (info /= 0) exit +!!$ write(0,*) 'Calling copyout ',nzrmax,nlw,nup,nidx,l2 + call psb_z_invt_copyout(nzrmax,sp_thresh,i,n,nlw,nup,jmaxup,nrmi,row,& + & nidx,idxs,l2,zcsr%ja,zcsr%irp,zcsr%val,info) + if (info /= 0) exit + nzz = l2 + end do outer + if (info /= psb_success_) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='mainloop') + goto 9999 + end if + + ipz1 = nzz+1 + call psb_ensure_size(ipz1,zcsr%val,info) + call psb_ensure_size(ipz1,zcsr%ja,info) + zcsr%val(ipz1) = zone + zcsr%ja(ipz1) = n + zcsr%irp(n+1) = ipz1+1 + + call z%mv_from(zcsr) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_zsparse_invt + +subroutine psb_z_invt_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,nrmi,row,heap,& + & irwt,ktrw,trw,info,sign) + use psb_base_mod + use psb_d_invt_fact_mod, psb_protect_name => psb_d_invt_copyin + implicit none + type(psb_z_csr_sparse_mat), intent(in) :: a + type(psb_z_coo_sparse_mat), intent(inout) :: trw + integer(psb_ipk_), intent(in) :: i, m,jmin,jmax,jd + integer(psb_ipk_), intent(inout) :: ktrw,nlw,nup,jmaxup,info + integer(psb_ipk_), intent(inout) :: irwt(:) + real(psb_dpk_), intent(inout) :: nrmi + complex(psb_dpk_), intent(inout) :: row(:) + type(psb_i_heap), intent(inout) :: heap + real(psb_dpk_), intent(in), optional :: sign + ! + integer(psb_ipk_) :: k,j,irb,kin,nz, err_act + integer(psb_ipk_), parameter :: nrb=16 + real(psb_dpk_) :: dmaxup, sign_ + real(psb_dpk_), external :: dnrm2 + character(len=20), parameter :: name='invt_copyin' + + info = psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + + call heap%init(info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_init_heap') + goto 9999 + end if + sign_ = done + if (present(sign)) sign_ = sign + ! + ! nrmi is the norm of the current sparse row (for the time being, + ! we use the 2-norm). + ! NOTE: the 2-norm below includes also elements that are outside + ! [jmin:jmax] strictly. Is this really important? TO BE CHECKED. + ! + + nlw = 0 + nup = 0 + jmaxup = 0 + dmaxup = zzero + nrmi = zzero + + do j = a%irp(i), a%irp(i+1) - 1 + k = a%ja(j) + if ((jmin<=k).and.(k<=jmax)) then + row(k) = sign_ * a%val(j) + call heap%insert(k,info) + irwt(k) = 1 + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_insert_heap') + goto 9999 + end if + end if + if (kjd) then + nup = nup + 1 + if (abs(row(k))>dmaxup) then + jmaxup = k + dmaxup = abs(row(k)) + end if + end if + end do + nz = a%irp(i+1) - a%irp(i) + nrmi = dnrm2(nz,a%val(a%irp(i):),ione) + + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_z_invt_copyin + +subroutine psb_z_invt_copyout(fill_in,thres,i,m,nlw,nup,jmaxup,nrmi,row, & + & nidx,idxs,l2,ja,irp,val,info) + + use psb_base_mod + use psb_z_invt_fact_mod, psb_protect_name => psb_z_invt_copyout + + implicit none + + ! Arguments + integer(psb_ipk_), intent(in) :: fill_in,i,m,nidx,nlw,nup,jmaxup + integer(psb_ipk_), intent(in) :: idxs(:) + integer(psb_ipk_), intent(inout) :: l2, info + integer(psb_ipk_), allocatable, intent(inout) :: ja(:),irp(:) + real(psb_dpk_), intent(in) :: thres,nrmi + complex(psb_dpk_),allocatable, intent(inout) :: val(:) + complex(psb_dpk_), intent(inout) :: row(:) + + ! Local variables + real(psb_dpk_),allocatable :: xw(:) + integer(psb_ipk_), allocatable :: xwid(:), indx(:) + real(psb_dpk_) :: witem, wmin + integer(psb_ipk_) :: widx + integer(psb_ipk_) :: k,isz,err_act,int_err(5),idxp, nz + type(psb_d_idx_heap) :: heap + character(len=20), parameter :: name='invt_copyout' + character(len=20) :: ch_err + logical :: fndmaxup + + info = psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_; goto 9999 + end if + + ! + ! Here we need to apply also the dropping rule base on the fill-in. + ! We do it by putting into a heap the elements that are not dropped + ! by using the 2-norm rule, and then copying them out. + ! + ! The heap goes down on the entry absolute value, so the first item + ! is the largest absolute value. + ! +!!$ write(0,*) 'invt_copyout ',nidx,nup+fill_in + call heap%init(info,dir=psb_asort_down_) + + if (info == psb_success_) allocate(xwid(nidx),xw(nidx),indx(nidx),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/3*nidx/),& + & a_err='real(psb_dpk_)') + goto 9999 + end if + + ! + ! First the lower part + ! + + nz = 0 + idxp = 0 + + do + + idxp = idxp + 1 + if (idxp > nidx) exit + if (idxs(idxp) >= i) exit + widx = idxs(idxp) + witem = row(widx) + ! + ! Dropping rule based on the 2-norm + ! + if (abs(witem) < thres*nrmi) cycle + + nz = nz + 1 + xw(nz) = witem + xwid(nz) = widx + call heap%insert(witem,widx,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_insert_heap') + goto 9999 + end if + end do + + if (nz > 1) then + write(psb_err_unit,*) 'Warning: lower triangle from invt???? ' + end if + + + if (idxp <= size(idxs)) then + if (idxs(idxp) < i) then + do + idxp = idxp + 1 + if (idxp > nidx) exit + if (idxs(idxp) >= i) exit + end do + end if + end if + idxp = idxp - 1 + nz = 0 + wmin=HUGE(wmin) + if (.false.) then + do + + idxp = idxp + 1 + if (idxp > nidx) exit + widx = idxs(idxp) + if (widx < i) then + write(psb_err_unit,*) 'Warning: lower triangle in upper copy',widx,i,idxp,idxs(idxp) + cycle + end if + if (widx > m) then + cycle + end if + witem = row(widx) + ! + ! Dropping rule based on the 2-norm. But keep the jmaxup-th entry anyway. + ! + if ((widx /= jmaxup) .and. (widx /= i) .and. (abs(witem) < thres*nrmi)) then + cycle + end if + if ((widx/=jmaxup).and.(nz > nup+fill_in)) then + if (abs(witem) < wmin) cycle + endif + wmin = min(abs(witem),wmin) + nz = nz + 1 + xw(nz) = witem + xwid(nz) = widx + call heap%insert(witem,widx,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_insert_heap') + goto 9999 + end if + + end do + + ! + ! Now we have to take out the first nup-fill_in entries. But make sure + ! we include entry jmaxup. + ! + if (nz <= nup+fill_in) then + ! + ! Just copy everything from xw + ! + fndmaxup=.true. + else + fndmaxup = .false. + nz = nup+fill_in + do k=1,nz + call heap%get_first(witem,widx,info) + xw(k) = witem + xwid(k) = widx + if (widx == jmaxup) fndmaxup=.true. + end do + end if + if ((i nidx) exit + widx = idxs(idxp) + if (widx < i) then + write(psb_err_unit,*) 'Warning: lower triangle in upper copy',widx,i,idxp,idxs(idxp) + cycle + end if + if (widx > m) then + cycle + end if + witem = row(widx) + ! + ! Dropping rule based on the 2-norm. But keep the jmaxup-th entry anyway. + ! + if ((widx /= i) .and. (abs(witem) < thres*nrmi)) then + cycle + end if + if (nz > nup+fill_in) then + if (abs(witem) < wmin) cycle + endif + wmin = min(abs(witem),wmin) + nz = nz + 1 + xw(nz) = witem + xwid(nz) = widx + call heap%insert(witem,widx,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_insert_heap') + goto 9999 + end if + + end do + + ! + ! Now we have to take out the first nup-fill_in entries. But make sure + ! we include entry jmaxup. + ! + if (nz > nup+fill_in) then + nz = nup+fill_in + do k=1,nz + call heap%get_first(witem,widx,info) + xw(k) = witem + xwid(k) = widx + end do + end if + end if + + ! + ! Now we put things back into ascending column order + ! + call psb_msort(xwid(1:nz),indx(1:nz),dir=psb_sort_up_) + + ! + ! Copy out the upper part of the row + ! + do k=1,nz + l2 = l2 + 1 + if (size(val) < l2) then + ! + ! Figure out a good reallocation size! + ! + isz = max(int(1.2*l2),l2+100) + call psb_realloc(isz,val,info) + if (info == psb_success_) call psb_realloc(isz,ja,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='Allocate') + goto 9999 + end if + end if + ja(l2) = xwid(k) + val(l2) = xw(indx(k)) + end do + + ! + ! Set row to zero + ! + do idxp=1,nidx + row(idxs(idxp)) = zzero + end do + + irp(i+1) = l2 + 1 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_z_invt_copyout diff --git a/prec/psb_c_bjacprec.f90 b/prec/psb_c_bjacprec.f90 index fb1d5429..847baf00 100644 --- a/prec/psb_c_bjacprec.f90 +++ b/prec/psb_c_bjacprec.f90 @@ -59,8 +59,9 @@ module psb_c_bjacprec end type psb_c_bjac_prec_type character(len=15), parameter, private :: & - & fact_names(0:4)=(/'None ','ILU(0) ',& - & 'ILU(n) ','ILU(eps) ','AINV(eps) '/) + & fact_names(0:5)=(/'None ','ILU(0) ',& + & 'ILU(n) ','ILU(eps) ','AINV(eps) ',& + & 'INVT '/) private :: psb_c_bjac_sizeof, psb_c_bjac_precdescr, psb_c_bjac_get_nzeros diff --git a/prec/psb_c_invk_fact_mod.f90 b/prec/psb_c_invk_fact_mod.f90 new file mode 100644 index 00000000..620a8adf --- /dev/null +++ b/prec/psb_c_invk_fact_mod.f90 @@ -0,0 +1,165 @@ +! +! 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 AMG4PSBLAS, 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. +! +! +! File: psb_c_invk_fact_mod.f90 +! +! Module: psb_c_invk_fact_mod +! +! This module defines some interfaces used internally by the implementation of +! psb_c_invk_solver, but not visible to the end user. +! +! +module psb_c_invk_fact_mod + + use psb_base_mod + use psb_prec_const_mod + + interface psb_invk_fact + subroutine psb_c_invk_bld(a,fill1, fill2,lmat,d,umat,desc,info,blck) + ! import + import psb_cspmat_type, psb_ipk_, psb_spk_, psb_desc_type + ! Arguments + type(psb_cspmat_type), intent(in), target :: a + integer(psb_ipk_), intent(in) :: fill1, fill2 + type(psb_cspmat_type), intent(inout) :: lmat, umat + complex(psb_spk_), allocatable :: d(:) + Type(psb_desc_type), Intent(inout) :: desc + integer(psb_ipk_), intent(out) :: info + type(psb_cspmat_type), intent(in), optional :: blck + end subroutine + end interface + + ! Inner workings + interface psb_sparse_invk + subroutine psb_csparse_invk(n,a,z,fill_in,info,inlevs) + ! Import + import psb_ipk_, psb_cspmat_type + ! Arguments + integer(psb_ipk_), intent(in) :: n + type(psb_cspmat_type), intent(in) :: a + type(psb_cspmat_type), intent(inout) :: z + integer(psb_ipk_), intent(in) :: fill_in + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: inlevs(:) + end subroutine psb_csparse_invk + end interface + + interface psb_invk_inv + subroutine psb_cinvk_inv(fill_in,i,row,rowlevs,heap,uia1,uia2,uaspk,uplevs,& + & nidx,idxs,info) + + use psb_base_mod, only : psb_cspmat_type, psb_spk_, psb_i_heap, psb_ipk_ + implicit none + + ! Arguments + type(psb_i_heap), intent(inout) :: heap + integer(psb_ipk_), intent(in) :: i, fill_in + integer(psb_ipk_), intent(inout) :: nidx,info + integer(psb_ipk_), intent(inout) :: rowlevs(:) + integer(psb_ipk_), allocatable, intent(inout) :: idxs(:) + integer(psb_ipk_), intent(in) :: uia1(:),uia2(:),uplevs(:) + complex(psb_spk_), intent(in) :: uaspk(:) + complex(psb_spk_), intent(inout) :: row(:) + + + end subroutine psb_cinvk_inv + end interface + + interface psb_invk_copyin + subroutine psb_c_invk_copyin(i,m,a,jmin,jmax,row,rowlevs,heap,& + & ktrw,trw,info,sign,inlevs) + ! Import + use psb_base_mod, only : psb_c_csr_sparse_mat, psb_c_coo_sparse_mat,& + & psb_spk_, psb_i_heap, psb_ipk_ + implicit none + ! Arguments + type(psb_c_csr_sparse_mat), intent(in) :: a + type(psb_c_coo_sparse_mat), intent(inout) :: trw + integer(psb_ipk_), intent(in) :: i,m,jmin,jmax + integer(psb_ipk_), intent(inout) :: ktrw,info + integer(psb_ipk_), intent(inout) :: rowlevs(:) + complex(psb_spk_), intent(inout) :: row(:) + type(psb_i_heap), intent(inout) :: heap + real(psb_spk_), optional, intent(in) :: sign + integer(psb_ipk_), intent(in), optional :: inlevs(:) + end subroutine psb_c_invk_copyin + end interface + + interface psb_invk_copyout + subroutine psb_c_invk_copyout(fill_in,i,m,row,rowlevs,nidx,idxs,& + & l2,uia1,uia2,uaspk,info) + use psb_base_mod, only : psb_cspmat_type, psb_spk_, psb_i_heap, psb_ipk_ + implicit none + ! Arguments + integer(psb_ipk_), intent(in) :: fill_in, i, m, nidx + integer(psb_ipk_), intent(inout) :: l2, info + integer(psb_ipk_), intent(inout) :: rowlevs(:), idxs(:) + integer(psb_ipk_), allocatable, intent(inout) :: uia1(:), uia2(:) + complex(psb_spk_), allocatable, intent(inout) :: uaspk(:) + complex(psb_spk_), intent(inout) :: row(:) + end subroutine psb_c_invk_copyout + end interface + +end module diff --git a/prec/psb_c_invt_fact_mod.f90 b/prec/psb_c_invt_fact_mod.f90 new file mode 100644 index 00000000..d7bf2f14 --- /dev/null +++ b/prec/psb_c_invt_fact_mod.f90 @@ -0,0 +1,150 @@ +! +! 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 AMG4PSBLAS, 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. +! +! +! File: psb_c_invt_fact_mod.f90 +! +! Module: psb_c_invt_fact_mod +! +! This module defines some interfaces used internally by the implementation of +! psb_c_invt_solver, but not visible to the end user. +! +! +module psb_c_invt_fact_mod + + use psb_base_mod + use psb_prec_const_mod + + interface psb_invt_fact + subroutine psb_c_invt_bld(a,fillin,invfill,thresh,invthresh,& + & lmat,d,umat,desc,info,blck) + ! Import + import psb_cspmat_type, psb_spk_, psb_ipk_, psb_desc_type + ! Arguments + type(psb_cspmat_type), intent(in), target :: a + integer(psb_ipk_), intent(in) :: fillin,invfill + real(psb_spk_), intent(in) :: thresh + real(psb_spk_), intent(in) :: invthresh + type(psb_cspmat_type), intent(inout) :: lmat, umat + complex(psb_spk_), allocatable :: d(:) + Type(psb_desc_type), Intent(inout) :: desc + integer(psb_ipk_), intent(out) :: info + type(psb_cspmat_type), intent(in), optional :: blck + end subroutine psb_c_invt_bld + end interface + + ! Interfaces for the inner workings + + interface + subroutine psb_csparse_invt(n,a,z,nzrmax,sp_thresh,info) + ! Import + import psb_cspmat_type, psb_spk_, psb_ipk_ + ! Arguments + integer(psb_ipk_), intent(in) :: n + type(psb_cspmat_type), intent(in) :: a + type(psb_cspmat_type), intent(inout) :: z + integer(psb_ipk_), intent(in) :: nzrmax + real(psb_spk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(out) :: info + end subroutine psb_csparse_invt + end interface + + interface + subroutine psb_c_invt_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,nrmi,row,heap,& + & irwt,ktrw,trw,info,sign) + ! Import + import psb_c_csr_sparse_mat, psb_c_coo_sparse_mat, psb_spk_, & + & psb_ipk_, psb_i_heap + ! Arguments + type(psb_c_csr_sparse_mat), intent(in) :: a + type(psb_c_coo_sparse_mat), intent(inout) :: trw + integer(psb_ipk_), intent(in) :: i, m,jmin,jmax,jd + integer(psb_ipk_), intent(inout) :: ktrw,nlw,nup,jmaxup,info + integer(psb_ipk_), intent(inout) :: irwt(:) + real(psb_spk_), intent(inout) :: nrmi + complex(psb_spk_), intent(inout) :: row(:) + type(psb_i_heap), intent(inout) :: heap + real(psb_spk_), intent(in), optional :: sign + end subroutine psb_c_invt_copyin + end interface + + interface + subroutine psb_c_invt_copyout(fill_in,thres,i,m,nlw,nup,jmaxup,nrmi,row, & + & nidx,idxs,l2,ja,irp,val,info) + ! Import + import psb_spk_, psb_ipk_ + ! Arguments + integer(psb_ipk_), intent(in) :: fill_in,i,m,nidx,nlw,nup,jmaxup + integer(psb_ipk_), intent(in) :: idxs(:) + integer(psb_ipk_), intent(inout) :: l2, info + integer(psb_ipk_), allocatable, intent(inout) :: ja(:),irp(:) + real(psb_spk_), intent(in) :: thres,nrmi + complex(psb_spk_),allocatable, intent(inout) :: val(:) + complex(psb_spk_), intent(inout) :: row(:) + end subroutine psb_c_invt_copyout + end interface + +contains + +end module psb_c_invt_fact_mod diff --git a/prec/psb_d_bjacprec.f90 b/prec/psb_d_bjacprec.f90 index 3fb3e2a0..1cb90e05 100644 --- a/prec/psb_d_bjacprec.f90 +++ b/prec/psb_d_bjacprec.f90 @@ -59,8 +59,9 @@ module psb_d_bjacprec end type psb_d_bjac_prec_type character(len=15), parameter, private :: & - & fact_names(0:4)=(/'None ','ILU(0) ',& - & 'ILU(n) ','ILU(eps) ','AINV(eps) '/) + & fact_names(0:5)=(/'None ','ILU(0) ',& + & 'ILU(n) ','ILU(eps) ','AINV(eps) ',& + & 'INVT '/) private :: psb_d_bjac_sizeof, psb_d_bjac_precdescr, psb_d_bjac_get_nzeros diff --git a/prec/psb_d_invk_fact_mod.f90 b/prec/psb_d_invk_fact_mod.f90 new file mode 100644 index 00000000..2bd97198 --- /dev/null +++ b/prec/psb_d_invk_fact_mod.f90 @@ -0,0 +1,165 @@ +! +! 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 AMG4PSBLAS, 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. +! +! +! File: psb_d_invk_fact_mod.f90 +! +! Module: psb_d_invk_fact_mod +! +! This module defines some interfaces used internally by the implementation of +! psb_d_invk_solver, but not visible to the end user. +! +! +module psb_d_invk_fact_mod + + use psb_base_mod + use psb_prec_const_mod + + interface psb_invk_fact + subroutine psb_d_invk_bld(a,fill1, fill2,lmat,d,umat,desc,info,blck) + ! import + import psb_dspmat_type, psb_ipk_, psb_dpk_, psb_desc_type + ! Arguments + type(psb_dspmat_type), intent(in), target :: a + integer(psb_ipk_), intent(in) :: fill1, fill2 + type(psb_dspmat_type), intent(inout) :: lmat, umat + real(psb_dpk_), allocatable :: d(:) + Type(psb_desc_type), Intent(inout) :: desc + integer(psb_ipk_), intent(out) :: info + type(psb_dspmat_type), intent(in), optional :: blck + end subroutine + end interface + + ! Inner workings + interface psb_sparse_invk + subroutine psb_dsparse_invk(n,a,z,fill_in,info,inlevs) + ! Import + import psb_ipk_, psb_dspmat_type + ! Arguments + integer(psb_ipk_), intent(in) :: n + type(psb_dspmat_type), intent(in) :: a + type(psb_dspmat_type), intent(inout) :: z + integer(psb_ipk_), intent(in) :: fill_in + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: inlevs(:) + end subroutine psb_dsparse_invk + end interface + + interface psb_invk_inv + subroutine psb_dinvk_inv(fill_in,i,row,rowlevs,heap,uia1,uia2,uaspk,uplevs,& + & nidx,idxs,info) + + use psb_base_mod, only : psb_dspmat_type, psb_dpk_, psb_i_heap, psb_ipk_ + implicit none + + ! Arguments + type(psb_i_heap), intent(inout) :: heap + integer(psb_ipk_), intent(in) :: i, fill_in + integer(psb_ipk_), intent(inout) :: nidx,info + integer(psb_ipk_), intent(inout) :: rowlevs(:) + integer(psb_ipk_), allocatable, intent(inout) :: idxs(:) + integer(psb_ipk_), intent(in) :: uia1(:),uia2(:),uplevs(:) + real(psb_dpk_), intent(in) :: uaspk(:) + real(psb_dpk_), intent(inout) :: row(:) + + + end subroutine psb_dinvk_inv + end interface + + interface psb_invk_copyin + subroutine psb_d_invk_copyin(i,m,a,jmin,jmax,row,rowlevs,heap,& + & ktrw,trw,info,sign,inlevs) + ! Import + use psb_base_mod, only : psb_d_csr_sparse_mat, psb_d_coo_sparse_mat,& + & psb_dpk_, psb_i_heap, psb_ipk_ + implicit none + ! Arguments + type(psb_d_csr_sparse_mat), intent(in) :: a + type(psb_d_coo_sparse_mat), intent(inout) :: trw + integer(psb_ipk_), intent(in) :: i,m,jmin,jmax + integer(psb_ipk_), intent(inout) :: ktrw,info + integer(psb_ipk_), intent(inout) :: rowlevs(:) + real(psb_dpk_), intent(inout) :: row(:) + type(psb_i_heap), intent(inout) :: heap + real(psb_dpk_), optional, intent(in) :: sign + integer(psb_ipk_), intent(in), optional :: inlevs(:) + end subroutine psb_d_invk_copyin + end interface + + interface psb_invk_copyout + subroutine psb_d_invk_copyout(fill_in,i,m,row,rowlevs,nidx,idxs,& + & l2,uia1,uia2,uaspk,info) + use psb_base_mod, only : psb_dspmat_type, psb_dpk_, psb_i_heap, psb_ipk_ + implicit none + ! Arguments + integer(psb_ipk_), intent(in) :: fill_in, i, m, nidx + integer(psb_ipk_), intent(inout) :: l2, info + integer(psb_ipk_), intent(inout) :: rowlevs(:), idxs(:) + integer(psb_ipk_), allocatable, intent(inout) :: uia1(:), uia2(:) + real(psb_dpk_), allocatable, intent(inout) :: uaspk(:) + real(psb_dpk_), intent(inout) :: row(:) + end subroutine psb_d_invk_copyout + end interface + +end module diff --git a/prec/psb_d_invt_fact_mod.f90 b/prec/psb_d_invt_fact_mod.f90 new file mode 100644 index 00000000..e954d938 --- /dev/null +++ b/prec/psb_d_invt_fact_mod.f90 @@ -0,0 +1,150 @@ +! +! 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 AMG4PSBLAS, 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. +! +! +! File: psb_d_invt_fact_mod.f90 +! +! Module: psb_d_invt_fact_mod +! +! This module defines some interfaces used internally by the implementation of +! psb_d_invt_solver, but not visible to the end user. +! +! +module psb_d_invt_fact_mod + + use psb_base_mod + use psb_prec_const_mod + + interface psb_invt_fact + subroutine psb_d_invt_bld(a,fillin,invfill,thresh,invthresh,& + & lmat,d,umat,desc,info,blck) + ! Import + import psb_dspmat_type, psb_dpk_, psb_ipk_, psb_desc_type + ! Arguments + type(psb_dspmat_type), intent(in), target :: a + integer(psb_ipk_), intent(in) :: fillin,invfill + real(psb_dpk_), intent(in) :: thresh + real(psb_dpk_), intent(in) :: invthresh + type(psb_dspmat_type), intent(inout) :: lmat, umat + real(psb_dpk_), allocatable :: d(:) + Type(psb_desc_type), Intent(inout) :: desc + integer(psb_ipk_), intent(out) :: info + type(psb_dspmat_type), intent(in), optional :: blck + end subroutine psb_d_invt_bld + end interface + + ! Interfaces for the inner workings + + interface + subroutine psb_dsparse_invt(n,a,z,nzrmax,sp_thresh,info) + ! Import + import psb_dspmat_type, psb_dpk_, psb_ipk_ + ! Arguments + integer(psb_ipk_), intent(in) :: n + type(psb_dspmat_type), intent(in) :: a + type(psb_dspmat_type), intent(inout) :: z + integer(psb_ipk_), intent(in) :: nzrmax + real(psb_dpk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(out) :: info + end subroutine psb_dsparse_invt + end interface + + interface + subroutine psb_d_invt_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,nrmi,row,heap,& + & irwt,ktrw,trw,info,sign) + ! Import + import psb_d_csr_sparse_mat, psb_d_coo_sparse_mat, psb_dpk_, & + & psb_ipk_, psb_i_heap + ! Arguments + type(psb_d_csr_sparse_mat), intent(in) :: a + type(psb_d_coo_sparse_mat), intent(inout) :: trw + integer(psb_ipk_), intent(in) :: i, m,jmin,jmax,jd + integer(psb_ipk_), intent(inout) :: ktrw,nlw,nup,jmaxup,info + integer(psb_ipk_), intent(inout) :: irwt(:) + real(psb_dpk_), intent(inout) :: nrmi + real(psb_dpk_), intent(inout) :: row(:) + type(psb_i_heap), intent(inout) :: heap + real(psb_dpk_), intent(in), optional :: sign + end subroutine psb_d_invt_copyin + end interface + + interface + subroutine psb_d_invt_copyout(fill_in,thres,i,m,nlw,nup,jmaxup,nrmi,row, & + & nidx,idxs,l2,ja,irp,val,info) + ! Import + import psb_dpk_, psb_ipk_ + ! Arguments + integer(psb_ipk_), intent(in) :: fill_in,i,m,nidx,nlw,nup,jmaxup + integer(psb_ipk_), intent(in) :: idxs(:) + integer(psb_ipk_), intent(inout) :: l2, info + integer(psb_ipk_), allocatable, intent(inout) :: ja(:),irp(:) + real(psb_dpk_), intent(in) :: thres,nrmi + real(psb_dpk_),allocatable, intent(inout) :: val(:) + real(psb_dpk_), intent(inout) :: row(:) + end subroutine psb_d_invt_copyout + end interface + +contains + +end module psb_d_invt_fact_mod diff --git a/prec/psb_prec_const_mod.f90 b/prec/psb_prec_const_mod.f90 index b21ae9da..d8d9a43f 100644 --- a/prec/psb_prec_const_mod.f90 +++ b/prec/psb_prec_const_mod.f90 @@ -58,7 +58,7 @@ module psb_prec_const_mod ! Factorization types: none, ILU(0), ILU(N), ILU(N,E) integer(psb_ipk_), parameter :: psb_f_none_=0,psb_f_ilu_n_=1,psb_f_ilu_k_=2,psb_f_ilu_t_=3 ! Approximate Inverse factorization type: AINV - integer(psb_ipk_), parameter :: psb_f_ainv_=4 + integer(psb_ipk_), parameter :: psb_f_ainv_=4, psb_f_invt_=5 ! Fields for sparse matrices ensembles: integer(psb_ipk_), parameter :: psb_l_pr_=1, psb_u_pr_=2, psb_bp_ilu_avsz=2 integer(psb_ipk_), parameter :: psb_max_avsz=psb_bp_ilu_avsz diff --git a/prec/psb_s_bjacprec.f90 b/prec/psb_s_bjacprec.f90 index 6c7e9c9e..028aabe8 100644 --- a/prec/psb_s_bjacprec.f90 +++ b/prec/psb_s_bjacprec.f90 @@ -59,8 +59,9 @@ module psb_s_bjacprec end type psb_s_bjac_prec_type character(len=15), parameter, private :: & - & fact_names(0:4)=(/'None ','ILU(0) ',& - & 'ILU(n) ','ILU(eps) ','AINV(eps) '/) + & fact_names(0:5)=(/'None ','ILU(0) ',& + & 'ILU(n) ','ILU(eps) ','AINV(eps) ',& + & 'INVT '/) private :: psb_s_bjac_sizeof, psb_s_bjac_precdescr, psb_s_bjac_get_nzeros diff --git a/prec/psb_s_invk_fact_mod.f90 b/prec/psb_s_invk_fact_mod.f90 new file mode 100644 index 00000000..6b0d3553 --- /dev/null +++ b/prec/psb_s_invk_fact_mod.f90 @@ -0,0 +1,165 @@ +! +! 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 AMG4PSBLAS, 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. +! +! +! File: psb_s_invk_fact_mod.f90 +! +! Module: psb_s_invk_fact_mod +! +! This module defines some interfaces used internally by the implementation of +! psb_s_invk_solver, but not visible to the end user. +! +! +module psb_s_invk_fact_mod + + use psb_base_mod + use psb_prec_const_mod + + interface psb_invk_fact + subroutine psb_s_invk_bld(a,fill1, fill2,lmat,d,umat,desc,info,blck) + ! import + import psb_sspmat_type, psb_ipk_, psb_spk_, psb_desc_type + ! Arguments + type(psb_sspmat_type), intent(in), target :: a + integer(psb_ipk_), intent(in) :: fill1, fill2 + type(psb_sspmat_type), intent(inout) :: lmat, umat + real(psb_spk_), allocatable :: d(:) + Type(psb_desc_type), Intent(inout) :: desc + integer(psb_ipk_), intent(out) :: info + type(psb_sspmat_type), intent(in), optional :: blck + end subroutine + end interface + + ! Inner workings + interface psb_sparse_invk + subroutine psb_ssparse_invk(n,a,z,fill_in,info,inlevs) + ! Import + import psb_ipk_, psb_sspmat_type + ! Arguments + integer(psb_ipk_), intent(in) :: n + type(psb_sspmat_type), intent(in) :: a + type(psb_sspmat_type), intent(inout) :: z + integer(psb_ipk_), intent(in) :: fill_in + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: inlevs(:) + end subroutine psb_ssparse_invk + end interface + + interface psb_invk_inv + subroutine psb_sinvk_inv(fill_in,i,row,rowlevs,heap,uia1,uia2,uaspk,uplevs,& + & nidx,idxs,info) + + use psb_base_mod, only : psb_sspmat_type, psb_spk_, psb_i_heap, psb_ipk_ + implicit none + + ! Arguments + type(psb_i_heap), intent(inout) :: heap + integer(psb_ipk_), intent(in) :: i, fill_in + integer(psb_ipk_), intent(inout) :: nidx,info + integer(psb_ipk_), intent(inout) :: rowlevs(:) + integer(psb_ipk_), allocatable, intent(inout) :: idxs(:) + integer(psb_ipk_), intent(in) :: uia1(:),uia2(:),uplevs(:) + real(psb_spk_), intent(in) :: uaspk(:) + real(psb_spk_), intent(inout) :: row(:) + + + end subroutine psb_sinvk_inv + end interface + + interface psb_invk_copyin + subroutine psb_s_invk_copyin(i,m,a,jmin,jmax,row,rowlevs,heap,& + & ktrw,trw,info,sign,inlevs) + ! Import + use psb_base_mod, only : psb_s_csr_sparse_mat, psb_s_coo_sparse_mat,& + & psb_spk_, psb_i_heap, psb_ipk_ + implicit none + ! Arguments + type(psb_s_csr_sparse_mat), intent(in) :: a + type(psb_s_coo_sparse_mat), intent(inout) :: trw + integer(psb_ipk_), intent(in) :: i,m,jmin,jmax + integer(psb_ipk_), intent(inout) :: ktrw,info + integer(psb_ipk_), intent(inout) :: rowlevs(:) + real(psb_spk_), intent(inout) :: row(:) + type(psb_i_heap), intent(inout) :: heap + real(psb_spk_), optional, intent(in) :: sign + integer(psb_ipk_), intent(in), optional :: inlevs(:) + end subroutine psb_s_invk_copyin + end interface + + interface psb_invk_copyout + subroutine psb_s_invk_copyout(fill_in,i,m,row,rowlevs,nidx,idxs,& + & l2,uia1,uia2,uaspk,info) + use psb_base_mod, only : psb_sspmat_type, psb_spk_, psb_i_heap, psb_ipk_ + implicit none + ! Arguments + integer(psb_ipk_), intent(in) :: fill_in, i, m, nidx + integer(psb_ipk_), intent(inout) :: l2, info + integer(psb_ipk_), intent(inout) :: rowlevs(:), idxs(:) + integer(psb_ipk_), allocatable, intent(inout) :: uia1(:), uia2(:) + real(psb_spk_), allocatable, intent(inout) :: uaspk(:) + real(psb_spk_), intent(inout) :: row(:) + end subroutine psb_s_invk_copyout + end interface + +end module diff --git a/prec/psb_s_invt_fact_mod.f90 b/prec/psb_s_invt_fact_mod.f90 new file mode 100644 index 00000000..4c147a9e --- /dev/null +++ b/prec/psb_s_invt_fact_mod.f90 @@ -0,0 +1,150 @@ +! +! 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 AMG4PSBLAS, 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. +! +! +! File: psb_s_invt_fact_mod.f90 +! +! Module: psb_s_invt_fact_mod +! +! This module defines some interfaces used internally by the implementation of +! psb_s_invt_solver, but not visible to the end user. +! +! +module psb_s_invt_fact_mod + + use psb_base_mod + use psb_prec_const_mod + + interface psb_invt_fact + subroutine psb_s_invt_bld(a,fillin,invfill,thresh,invthresh,& + & lmat,d,umat,desc,info,blck) + ! Import + import psb_sspmat_type, psb_spk_, psb_ipk_, psb_desc_type + ! Arguments + type(psb_sspmat_type), intent(in), target :: a + integer(psb_ipk_), intent(in) :: fillin,invfill + real(psb_spk_), intent(in) :: thresh + real(psb_spk_), intent(in) :: invthresh + type(psb_sspmat_type), intent(inout) :: lmat, umat + real(psb_spk_), allocatable :: d(:) + Type(psb_desc_type), Intent(inout) :: desc + integer(psb_ipk_), intent(out) :: info + type(psb_sspmat_type), intent(in), optional :: blck + end subroutine psb_s_invt_bld + end interface + + ! Interfaces for the inner workings + + interface + subroutine psb_ssparse_invt(n,a,z,nzrmax,sp_thresh,info) + ! Import + import psb_sspmat_type, psb_spk_, psb_ipk_ + ! Arguments + integer(psb_ipk_), intent(in) :: n + type(psb_sspmat_type), intent(in) :: a + type(psb_sspmat_type), intent(inout) :: z + integer(psb_ipk_), intent(in) :: nzrmax + real(psb_spk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(out) :: info + end subroutine psb_ssparse_invt + end interface + + interface + subroutine psb_s_invt_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,nrmi,row,heap,& + & irwt,ktrw,trw,info,sign) + ! Import + import psb_s_csr_sparse_mat, psb_s_coo_sparse_mat, psb_spk_, & + & psb_ipk_, psb_i_heap + ! Arguments + type(psb_s_csr_sparse_mat), intent(in) :: a + type(psb_s_coo_sparse_mat), intent(inout) :: trw + integer(psb_ipk_), intent(in) :: i, m,jmin,jmax,jd + integer(psb_ipk_), intent(inout) :: ktrw,nlw,nup,jmaxup,info + integer(psb_ipk_), intent(inout) :: irwt(:) + real(psb_spk_), intent(inout) :: nrmi + real(psb_spk_), intent(inout) :: row(:) + type(psb_i_heap), intent(inout) :: heap + real(psb_spk_), intent(in), optional :: sign + end subroutine psb_s_invt_copyin + end interface + + interface + subroutine psb_s_invt_copyout(fill_in,thres,i,m,nlw,nup,jmaxup,nrmi,row, & + & nidx,idxs,l2,ja,irp,val,info) + ! Import + import psb_spk_, psb_ipk_ + ! Arguments + integer(psb_ipk_), intent(in) :: fill_in,i,m,nidx,nlw,nup,jmaxup + integer(psb_ipk_), intent(in) :: idxs(:) + integer(psb_ipk_), intent(inout) :: l2, info + integer(psb_ipk_), allocatable, intent(inout) :: ja(:),irp(:) + real(psb_spk_), intent(in) :: thres,nrmi + real(psb_spk_),allocatable, intent(inout) :: val(:) + real(psb_spk_), intent(inout) :: row(:) + end subroutine psb_s_invt_copyout + end interface + +contains + +end module psb_s_invt_fact_mod diff --git a/prec/psb_z_bjacprec.f90 b/prec/psb_z_bjacprec.f90 index 23e826b1..4ca879db 100644 --- a/prec/psb_z_bjacprec.f90 +++ b/prec/psb_z_bjacprec.f90 @@ -59,8 +59,9 @@ module psb_z_bjacprec end type psb_z_bjac_prec_type character(len=15), parameter, private :: & - & fact_names(0:4)=(/'None ','ILU(0) ',& - & 'ILU(n) ','ILU(eps) ','AINV(eps) '/) + & fact_names(0:5)=(/'None ','ILU(0) ',& + & 'ILU(n) ','ILU(eps) ','AINV(eps) ',& + & 'INVT '/) private :: psb_z_bjac_sizeof, psb_z_bjac_precdescr, psb_z_bjac_get_nzeros diff --git a/prec/psb_z_invk_fact_mod.f90 b/prec/psb_z_invk_fact_mod.f90 new file mode 100644 index 00000000..0a1e5faf --- /dev/null +++ b/prec/psb_z_invk_fact_mod.f90 @@ -0,0 +1,165 @@ +! +! 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 AMG4PSBLAS, 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. +! +! +! File: psb_z_invk_fact_mod.f90 +! +! Module: psb_z_invk_fact_mod +! +! This module defines some interfaces used internally by the implementation of +! psb_z_invk_solver, but not visible to the end user. +! +! +module psb_z_invk_fact_mod + + use psb_base_mod + use psb_prec_const_mod + + interface psb_invk_fact + subroutine psb_z_invk_bld(a,fill1, fill2,lmat,d,umat,desc,info,blck) + ! import + import psb_zspmat_type, psb_ipk_, psb_dpk_, psb_desc_type + ! Arguments + type(psb_zspmat_type), intent(in), target :: a + integer(psb_ipk_), intent(in) :: fill1, fill2 + type(psb_zspmat_type), intent(inout) :: lmat, umat + complex(psb_dpk_), allocatable :: d(:) + Type(psb_desc_type), Intent(inout) :: desc + integer(psb_ipk_), intent(out) :: info + type(psb_zspmat_type), intent(in), optional :: blck + end subroutine + end interface + + ! Inner workings + interface psb_sparse_invk + subroutine psb_zsparse_invk(n,a,z,fill_in,info,inlevs) + ! Import + import psb_ipk_, psb_zspmat_type + ! Arguments + integer(psb_ipk_), intent(in) :: n + type(psb_zspmat_type), intent(in) :: a + type(psb_zspmat_type), intent(inout) :: z + integer(psb_ipk_), intent(in) :: fill_in + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: inlevs(:) + end subroutine psb_zsparse_invk + end interface + + interface psb_invk_inv + subroutine psb_zinvk_inv(fill_in,i,row,rowlevs,heap,uia1,uia2,uaspk,uplevs,& + & nidx,idxs,info) + + use psb_base_mod, only : psb_zspmat_type, psb_dpk_, psb_i_heap, psb_ipk_ + implicit none + + ! Arguments + type(psb_i_heap), intent(inout) :: heap + integer(psb_ipk_), intent(in) :: i, fill_in + integer(psb_ipk_), intent(inout) :: nidx,info + integer(psb_ipk_), intent(inout) :: rowlevs(:) + integer(psb_ipk_), allocatable, intent(inout) :: idxs(:) + integer(psb_ipk_), intent(in) :: uia1(:),uia2(:),uplevs(:) + complex(psb_dpk_), intent(in) :: uaspk(:) + complex(psb_dpk_), intent(inout) :: row(:) + + + end subroutine psb_zinvk_inv + end interface + + interface psb_invk_copyin + subroutine psb_z_invk_copyin(i,m,a,jmin,jmax,row,rowlevs,heap,& + & ktrw,trw,info,sign,inlevs) + ! Import + use psb_base_mod, only : psb_z_csr_sparse_mat, psb_z_coo_sparse_mat,& + & psb_dpk_, psb_i_heap, psb_ipk_ + implicit none + ! Arguments + type(psb_z_csr_sparse_mat), intent(in) :: a + type(psb_z_coo_sparse_mat), intent(inout) :: trw + integer(psb_ipk_), intent(in) :: i,m,jmin,jmax + integer(psb_ipk_), intent(inout) :: ktrw,info + integer(psb_ipk_), intent(inout) :: rowlevs(:) + complex(psb_dpk_), intent(inout) :: row(:) + type(psb_i_heap), intent(inout) :: heap + real(psb_dpk_), optional, intent(in) :: sign + integer(psb_ipk_), intent(in), optional :: inlevs(:) + end subroutine psb_z_invk_copyin + end interface + + interface psb_invk_copyout + subroutine psb_z_invk_copyout(fill_in,i,m,row,rowlevs,nidx,idxs,& + & l2,uia1,uia2,uaspk,info) + use psb_base_mod, only : psb_zspmat_type, psb_dpk_, psb_i_heap, psb_ipk_ + implicit none + ! Arguments + integer(psb_ipk_), intent(in) :: fill_in, i, m, nidx + integer(psb_ipk_), intent(inout) :: l2, info + integer(psb_ipk_), intent(inout) :: rowlevs(:), idxs(:) + integer(psb_ipk_), allocatable, intent(inout) :: uia1(:), uia2(:) + complex(psb_dpk_), allocatable, intent(inout) :: uaspk(:) + complex(psb_dpk_), intent(inout) :: row(:) + end subroutine psb_z_invk_copyout + end interface + +end module diff --git a/prec/psb_z_invt_fact_mod.f90 b/prec/psb_z_invt_fact_mod.f90 new file mode 100644 index 00000000..c0b3d4a3 --- /dev/null +++ b/prec/psb_z_invt_fact_mod.f90 @@ -0,0 +1,150 @@ +! +! 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 AMG4PSBLAS, 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. +! +! +! File: psb_z_invt_fact_mod.f90 +! +! Module: psb_z_invt_fact_mod +! +! This module defines some interfaces used internally by the implementation of +! psb_z_invt_solver, but not visible to the end user. +! +! +module psb_z_invt_fact_mod + + use psb_base_mod + use psb_prec_const_mod + + interface psb_invt_fact + subroutine psb_z_invt_bld(a,fillin,invfill,thresh,invthresh,& + & lmat,d,umat,desc,info,blck) + ! Import + import psb_zspmat_type, psb_dpk_, psb_ipk_, psb_desc_type + ! Arguments + type(psb_zspmat_type), intent(in), target :: a + integer(psb_ipk_), intent(in) :: fillin,invfill + real(psb_dpk_), intent(in) :: thresh + real(psb_dpk_), intent(in) :: invthresh + type(psb_zspmat_type), intent(inout) :: lmat, umat + complex(psb_dpk_), allocatable :: d(:) + Type(psb_desc_type), Intent(inout) :: desc + integer(psb_ipk_), intent(out) :: info + type(psb_zspmat_type), intent(in), optional :: blck + end subroutine psb_z_invt_bld + end interface + + ! Interfaces for the inner workings + + interface + subroutine psb_zsparse_invt(n,a,z,nzrmax,sp_thresh,info) + ! Import + import psb_zspmat_type, psb_dpk_, psb_ipk_ + ! Arguments + integer(psb_ipk_), intent(in) :: n + type(psb_zspmat_type), intent(in) :: a + type(psb_zspmat_type), intent(inout) :: z + integer(psb_ipk_), intent(in) :: nzrmax + real(psb_dpk_), intent(in) :: sp_thresh + integer(psb_ipk_), intent(out) :: info + end subroutine psb_zsparse_invt + end interface + + interface + subroutine psb_z_invt_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,nrmi,row,heap,& + & irwt,ktrw,trw,info,sign) + ! Import + import psb_z_csr_sparse_mat, psb_z_coo_sparse_mat, psb_dpk_, & + & psb_ipk_, psb_i_heap + ! Arguments + type(psb_z_csr_sparse_mat), intent(in) :: a + type(psb_z_coo_sparse_mat), intent(inout) :: trw + integer(psb_ipk_), intent(in) :: i, m,jmin,jmax,jd + integer(psb_ipk_), intent(inout) :: ktrw,nlw,nup,jmaxup,info + integer(psb_ipk_), intent(inout) :: irwt(:) + real(psb_dpk_), intent(inout) :: nrmi + complex(psb_dpk_), intent(inout) :: row(:) + type(psb_i_heap), intent(inout) :: heap + real(psb_dpk_), intent(in), optional :: sign + end subroutine psb_z_invt_copyin + end interface + + interface + subroutine psb_z_invt_copyout(fill_in,thres,i,m,nlw,nup,jmaxup,nrmi,row, & + & nidx,idxs,l2,ja,irp,val,info) + ! Import + import psb_dpk_, psb_ipk_ + ! Arguments + integer(psb_ipk_), intent(in) :: fill_in,i,m,nidx,nlw,nup,jmaxup + integer(psb_ipk_), intent(in) :: idxs(:) + integer(psb_ipk_), intent(inout) :: l2, info + integer(psb_ipk_), allocatable, intent(inout) :: ja(:),irp(:) + real(psb_dpk_), intent(in) :: thres,nrmi + complex(psb_dpk_),allocatable, intent(inout) :: val(:) + complex(psb_dpk_), intent(inout) :: row(:) + end subroutine psb_z_invt_copyout + end interface + +contains + +end module psb_z_invt_fact_mod