diff --git a/prec/Makefile b/prec/Makefile index 35420bc4..ec5892fe 100644 --- a/prec/Makefile +++ b/prec/Makefile @@ -56,10 +56,10 @@ psb_s_bjacprec.o psb_s_diagprec.o psb_s_nullprec.o: psb_prec_mod.o psb_s_base_pr psb_d_bjacprec.o psb_d_diagprec.o psb_d_nullprec.o: psb_prec_mod.o psb_d_base_prec_mod.o psb_c_bjacprec.o psb_c_diagprec.o psb_c_nullprec.o: psb_prec_mod.o psb_c_base_prec_mod.o psb_z_bjacprec.o psb_z_diagprec.o psb_z_nullprec.o: psb_prec_mod.o psb_z_base_prec_mod.o -psb_s_bjacprec.o: psb_s_ilu_fact_mod.o psb_s_ainv_fact_mod.o psb_s_invk_fact_mod psb_s_invt_fact_mod -psb_d_bjacprec.o: psb_d_ilu_fact_mod.o psb_d_ainv_fact_mod.o psb_d_invk_fact_mod psb_d_invt_fact_mod -psb_c_bjacprec.o: psb_c_ilu_fact_mod.o psb_c_ainv_fact_mod.o psb_c_invk_fact_mod psb_c_invt_fact_mod -psb_z_bjacprec.o: psb_z_ilu_fact_mod.o psb_z_ainv_fact_mod.o psb_z_invk_fact_mod psb_z_invt_fact_mod +psb_s_bjacprec.o: psb_s_ilu_fact_mod.o psb_s_ainv_fact_mod.o psb_s_invk_fact_mod.o psb_s_invt_fact_mod.o +psb_d_bjacprec.o: psb_d_ilu_fact_mod.o psb_d_ainv_fact_mod.o psb_d_invk_fact_mod.o psb_d_invt_fact_mod.o +psb_c_bjacprec.o: psb_c_ilu_fact_mod.o psb_c_ainv_fact_mod.o psb_c_invk_fact_mod.o psb_c_invt_fact_mod.o +psb_z_bjacprec.o: psb_z_ilu_fact_mod.o psb_z_ainv_fact_mod.o psb_z_invk_fact_mod.o psb_z_invt_fact_mod.o psb_d_ainv_fact_mod.o: psb_prec_const_mod.o psb_ainv_tools_mod.o psb_s_ainv_fact_mod.o: psb_prec_const_mod.o psb_ainv_tools_mod.o psb_c_ainv_fact_mod.o: psb_prec_const_mod.o psb_ainv_tools_mod.o diff --git a/prec/impl/psb_c_invt_fact.f90 b/prec/impl/psb_c_invt_fact.f90 index 272cb4ab..f0a5bd85 100644 --- a/prec/impl/psb_c_invt_fact.f90 +++ b/prec/impl/psb_c_invt_fact.f90 @@ -629,3 +629,110 @@ subroutine psb_c_invt_copyout(fill_in,thres,i,m,nlw,nup,jmaxup,nrmi,row, & 9999 call psb_error_handler(err_act) return end subroutine psb_c_invt_copyout + +subroutine psb_c_invt_inv(thres,i,nrmi,row,heap,irwt,ja,irp,val,nidx,idxs,info) + + use psb_base_mod + use psb_c_invt_fact_mod, psb_protect_name => psb_c_invt_inv + + implicit none + + ! Arguments + type(psb_i_heap), intent(inout) :: heap + integer(psb_ipk_), intent(in) :: i + integer(psb_ipk_), intent(inout) :: nidx,info + integer(psb_ipk_), intent(inout) :: irwt(:) + real(psb_spk_), intent(in) :: thres,nrmi + integer(psb_ipk_), allocatable, intent(inout) :: idxs(:) + integer(psb_ipk_), intent(in) :: ja(:),irp(:) + complex(psb_spk_), intent(in) :: val(:) + complex(psb_spk_), intent(inout) :: row(:) + + ! Local Variables + integer(psb_ipk_) :: k,j,jj,lastk,iret + real(psb_dpk_) :: rwk, alpha + + info = psb_success_ + + call psb_ensure_size(200, idxs, info) + if (info /= psb_success_) return + nidx = 1 + idxs(1) = i + lastk = i + irwt(i) = 1 +!!$ write(0,*) 'Drop Threshold ',thres*nrmi + ! + ! Do while there are indices to be processed + ! + do + + call heap%get_first(k,iret) + if (iret < 0) exit + + ! + ! An index may have been put on the heap more than once. + ! Should not happen, but just in case. + ! + if (k == lastk) cycle + lastk = k + + ! + ! Dropping rule based on the threshold: compare the absolute + ! value of each updated entry of row with thres * 2-norm of row. + ! + rwk = row(k) + + if (abs(rwk) < thres*nrmi) then + ! + ! Drop the entry. + ! + row(k) = dzero + irwt(k) = 0 + cycle + else + ! + ! Note: since U is scaled while copying it out (see ilut_copyout), + ! we can use rwk in the update below. + ! + do jj=irp(k),irp(k+1)-1 + j = ja(jj) + if (j<=k) then + info = -i + return + endif + ! + ! Update row(j) and, if it is not to be discarded, insert + ! its index into the heap for further processing. + ! + row(j) = row(j) - rwk * val(jj) + if (irwt(j) == 0) then + if (abs(row(j)) < thres*nrmi) then + ! + ! Drop the entry. + ! + row(j) = dzero + else + ! + ! Do the insertion. + ! + call heap%insert(j,info) + if (info /= psb_success_) return + irwt(j) = 1 + end if + end if + end do + end if + + ! + ! If we get here it is an index we need to keep on copyout. + ! + + nidx = nidx + 1 + call psb_ensure_size(nidx,idxs,info,addsz=psb_heap_resize) + if (info /= psb_success_) return + idxs(nidx) = k + irwt(k) = 0 + end do + + irwt(i) = 0 +end subroutine psb_c_invt_inv diff --git a/prec/impl/psb_d_invt_fact.f90 b/prec/impl/psb_d_invt_fact.f90 index 0312733c..21b11949 100644 --- a/prec/impl/psb_d_invt_fact.f90 +++ b/prec/impl/psb_d_invt_fact.f90 @@ -629,3 +629,110 @@ subroutine psb_d_invt_copyout(fill_in,thres,i,m,nlw,nup,jmaxup,nrmi,row, & 9999 call psb_error_handler(err_act) return end subroutine psb_d_invt_copyout + +subroutine psb_d_invt_inv(thres,i,nrmi,row,heap,irwt,ja,irp,val,nidx,idxs,info) + + use psb_base_mod + use psb_d_invt_fact_mod, psb_protect_name => psb_d_invt_inv + + implicit none + + ! Arguments + type(psb_i_heap), intent(inout) :: heap + integer(psb_ipk_), intent(in) :: i + integer(psb_ipk_), intent(inout) :: nidx,info + integer(psb_ipk_), intent(inout) :: irwt(:) + real(psb_dpk_), intent(in) :: thres,nrmi + integer(psb_ipk_), allocatable, intent(inout) :: idxs(:) + integer(psb_ipk_), intent(in) :: ja(:),irp(:) + real(psb_dpk_), intent(in) :: val(:) + real(psb_dpk_), intent(inout) :: row(:) + + ! Local Variables + integer(psb_ipk_) :: k,j,jj,lastk,iret + real(psb_dpk_) :: rwk, alpha + + info = psb_success_ + + call psb_ensure_size(200, idxs, info) + if (info /= psb_success_) return + nidx = 1 + idxs(1) = i + lastk = i + irwt(i) = 1 +!!$ write(0,*) 'Drop Threshold ',thres*nrmi + ! + ! Do while there are indices to be processed + ! + do + + call heap%get_first(k,iret) + if (iret < 0) exit + + ! + ! An index may have been put on the heap more than once. + ! Should not happen, but just in case. + ! + if (k == lastk) cycle + lastk = k + + ! + ! Dropping rule based on the threshold: compare the absolute + ! value of each updated entry of row with thres * 2-norm of row. + ! + rwk = row(k) + + if (abs(rwk) < thres*nrmi) then + ! + ! Drop the entry. + ! + row(k) = dzero + irwt(k) = 0 + cycle + else + ! + ! Note: since U is scaled while copying it out (see ilut_copyout), + ! we can use rwk in the update below. + ! + do jj=irp(k),irp(k+1)-1 + j = ja(jj) + if (j<=k) then + info = -i + return + endif + ! + ! Update row(j) and, if it is not to be discarded, insert + ! its index into the heap for further processing. + ! + row(j) = row(j) - rwk * val(jj) + if (irwt(j) == 0) then + if (abs(row(j)) < thres*nrmi) then + ! + ! Drop the entry. + ! + row(j) = dzero + else + ! + ! Do the insertion. + ! + call heap%insert(j,info) + if (info /= psb_success_) return + irwt(j) = 1 + end if + end if + end do + end if + + ! + ! If we get here it is an index we need to keep on copyout. + ! + + nidx = nidx + 1 + call psb_ensure_size(nidx,idxs,info,addsz=psb_heap_resize) + if (info /= psb_success_) return + idxs(nidx) = k + irwt(k) = 0 + end do + + irwt(i) = 0 +end subroutine psb_d_invt_inv diff --git a/prec/impl/psb_s_invt_fact.f90 b/prec/impl/psb_s_invt_fact.f90 index 3a261ca6..95bcd42e 100644 --- a/prec/impl/psb_s_invt_fact.f90 +++ b/prec/impl/psb_s_invt_fact.f90 @@ -629,3 +629,110 @@ subroutine psb_s_invt_copyout(fill_in,thres,i,m,nlw,nup,jmaxup,nrmi,row, & 9999 call psb_error_handler(err_act) return end subroutine psb_s_invt_copyout + +subroutine psb_s_invt_inv(thres,i,nrmi,row,heap,irwt,ja,irp,val,nidx,idxs,info) + + use psb_base_mod + use psb_s_invt_fact_mod, psb_protect_name => psb_s_invt_inv + + implicit none + + ! Arguments + type(psb_i_heap), intent(inout) :: heap + integer(psb_ipk_), intent(in) :: i + integer(psb_ipk_), intent(inout) :: nidx,info + integer(psb_ipk_), intent(inout) :: irwt(:) + real(psb_spk_), intent(in) :: thres,nrmi + integer(psb_ipk_), allocatable, intent(inout) :: idxs(:) + integer(psb_ipk_), intent(in) :: ja(:),irp(:) + real(psb_spk_), intent(in) :: val(:) + real(psb_spk_), intent(inout) :: row(:) + + ! Local Variables + integer(psb_ipk_) :: k,j,jj,lastk,iret + real(psb_dpk_) :: rwk, alpha + + info = psb_success_ + + call psb_ensure_size(200, idxs, info) + if (info /= psb_success_) return + nidx = 1 + idxs(1) = i + lastk = i + irwt(i) = 1 +!!$ write(0,*) 'Drop Threshold ',thres*nrmi + ! + ! Do while there are indices to be processed + ! + do + + call heap%get_first(k,iret) + if (iret < 0) exit + + ! + ! An index may have been put on the heap more than once. + ! Should not happen, but just in case. + ! + if (k == lastk) cycle + lastk = k + + ! + ! Dropping rule based on the threshold: compare the absolute + ! value of each updated entry of row with thres * 2-norm of row. + ! + rwk = row(k) + + if (abs(rwk) < thres*nrmi) then + ! + ! Drop the entry. + ! + row(k) = dzero + irwt(k) = 0 + cycle + else + ! + ! Note: since U is scaled while copying it out (see ilut_copyout), + ! we can use rwk in the update below. + ! + do jj=irp(k),irp(k+1)-1 + j = ja(jj) + if (j<=k) then + info = -i + return + endif + ! + ! Update row(j) and, if it is not to be discarded, insert + ! its index into the heap for further processing. + ! + row(j) = row(j) - rwk * val(jj) + if (irwt(j) == 0) then + if (abs(row(j)) < thres*nrmi) then + ! + ! Drop the entry. + ! + row(j) = dzero + else + ! + ! Do the insertion. + ! + call heap%insert(j,info) + if (info /= psb_success_) return + irwt(j) = 1 + end if + end if + end do + end if + + ! + ! If we get here it is an index we need to keep on copyout. + ! + + nidx = nidx + 1 + call psb_ensure_size(nidx,idxs,info,addsz=psb_heap_resize) + if (info /= psb_success_) return + idxs(nidx) = k + irwt(k) = 0 + end do + + irwt(i) = 0 +end subroutine psb_s_invt_inv diff --git a/prec/impl/psb_z_invt_fact.f90 b/prec/impl/psb_z_invt_fact.f90 index 17e38c67..74b557cd 100644 --- a/prec/impl/psb_z_invt_fact.f90 +++ b/prec/impl/psb_z_invt_fact.f90 @@ -629,3 +629,110 @@ subroutine psb_z_invt_copyout(fill_in,thres,i,m,nlw,nup,jmaxup,nrmi,row, & 9999 call psb_error_handler(err_act) return end subroutine psb_z_invt_copyout + +subroutine psb_z_invt_inv(thres,i,nrmi,row,heap,irwt,ja,irp,val,nidx,idxs,info) + + use psb_base_mod + use psb_z_invt_fact_mod, psb_protect_name => psb_z_invt_inv + + implicit none + + ! Arguments + type(psb_i_heap), intent(inout) :: heap + integer(psb_ipk_), intent(in) :: i + integer(psb_ipk_), intent(inout) :: nidx,info + integer(psb_ipk_), intent(inout) :: irwt(:) + real(psb_dpk_), intent(in) :: thres,nrmi + integer(psb_ipk_), allocatable, intent(inout) :: idxs(:) + integer(psb_ipk_), intent(in) :: ja(:),irp(:) + complex(psb_dpk_), intent(in) :: val(:) + complex(psb_dpk_), intent(inout) :: row(:) + + ! Local Variables + integer(psb_ipk_) :: k,j,jj,lastk,iret + real(psb_dpk_) :: rwk, alpha + + info = psb_success_ + + call psb_ensure_size(200, idxs, info) + if (info /= psb_success_) return + nidx = 1 + idxs(1) = i + lastk = i + irwt(i) = 1 +!!$ write(0,*) 'Drop Threshold ',thres*nrmi + ! + ! Do while there are indices to be processed + ! + do + + call heap%get_first(k,iret) + if (iret < 0) exit + + ! + ! An index may have been put on the heap more than once. + ! Should not happen, but just in case. + ! + if (k == lastk) cycle + lastk = k + + ! + ! Dropping rule based on the threshold: compare the absolute + ! value of each updated entry of row with thres * 2-norm of row. + ! + rwk = row(k) + + if (abs(rwk) < thres*nrmi) then + ! + ! Drop the entry. + ! + row(k) = dzero + irwt(k) = 0 + cycle + else + ! + ! Note: since U is scaled while copying it out (see ilut_copyout), + ! we can use rwk in the update below. + ! + do jj=irp(k),irp(k+1)-1 + j = ja(jj) + if (j<=k) then + info = -i + return + endif + ! + ! Update row(j) and, if it is not to be discarded, insert + ! its index into the heap for further processing. + ! + row(j) = row(j) - rwk * val(jj) + if (irwt(j) == 0) then + if (abs(row(j)) < thres*nrmi) then + ! + ! Drop the entry. + ! + row(j) = dzero + else + ! + ! Do the insertion. + ! + call heap%insert(j,info) + if (info /= psb_success_) return + irwt(j) = 1 + end if + end if + end do + end if + + ! + ! If we get here it is an index we need to keep on copyout. + ! + + nidx = nidx + 1 + call psb_ensure_size(nidx,idxs,info,addsz=psb_heap_resize) + if (info /= psb_success_) return + idxs(nidx) = k + irwt(k) = 0 + end do + + irwt(i) = 0 +end subroutine psb_z_invt_inv diff --git a/prec/psb_c_invt_fact_mod.f90 b/prec/psb_c_invt_fact_mod.f90 index d7bf2f14..841c39b1 100644 --- a/prec/psb_c_invt_fact_mod.f90 +++ b/prec/psb_c_invt_fact_mod.f90 @@ -145,6 +145,24 @@ module psb_c_invt_fact_mod end subroutine psb_c_invt_copyout end interface + interface psb_invt_inv + subroutine psb_c_invt_inv(thres,i,nrmi,row,heap,irwt,ja,irp,val, & + & nidx,idxs,info) + ! import + import psb_spk_, psb_i_heap, psb_ipk_ + ! Arguments + type(psb_i_heap), intent(inout) :: heap + integer(psb_ipk_), intent(in) :: i + integer(psb_ipk_), intent(inout) :: nidx,info + integer(psb_ipk_), intent(inout) :: irwt(:) + real(psb_spk_), intent(in) :: thres,nrmi + integer(psb_ipk_), allocatable, intent(inout) :: idxs(:) + integer(psb_ipk_), intent(in) :: ja(:),irp(:) + complex(psb_spk_), intent(in) :: val(:) + complex(psb_spk_), intent(inout) :: row(:) + end subroutine + end interface + contains end module psb_c_invt_fact_mod diff --git a/prec/psb_d_invt_fact_mod.f90 b/prec/psb_d_invt_fact_mod.f90 index e954d938..f38c1c2b 100644 --- a/prec/psb_d_invt_fact_mod.f90 +++ b/prec/psb_d_invt_fact_mod.f90 @@ -145,6 +145,24 @@ module psb_d_invt_fact_mod end subroutine psb_d_invt_copyout end interface + interface psb_invt_inv + subroutine psb_d_invt_inv(thres,i,nrmi,row,heap,irwt,ja,irp,val, & + & nidx,idxs,info) + ! import + import psb_dpk_, psb_i_heap, psb_ipk_ + ! Arguments + type(psb_i_heap), intent(inout) :: heap + integer(psb_ipk_), intent(in) :: i + integer(psb_ipk_), intent(inout) :: nidx,info + integer(psb_ipk_), intent(inout) :: irwt(:) + real(psb_dpk_), intent(in) :: thres,nrmi + integer(psb_ipk_), allocatable, intent(inout) :: idxs(:) + integer(psb_ipk_), intent(in) :: ja(:),irp(:) + real(psb_dpk_), intent(in) :: val(:) + real(psb_dpk_), intent(inout) :: row(:) + end subroutine + end interface + contains end module psb_d_invt_fact_mod diff --git a/prec/psb_s_invt_fact_mod.f90 b/prec/psb_s_invt_fact_mod.f90 index 4c147a9e..2c9ce38c 100644 --- a/prec/psb_s_invt_fact_mod.f90 +++ b/prec/psb_s_invt_fact_mod.f90 @@ -145,6 +145,24 @@ module psb_s_invt_fact_mod end subroutine psb_s_invt_copyout end interface + interface psb_invt_inv + subroutine psb_s_invt_inv(thres,i,nrmi,row,heap,irwt,ja,irp,val, & + & nidx,idxs,info) + ! import + import psb_spk_, psb_i_heap, psb_ipk_ + ! Arguments + type(psb_i_heap), intent(inout) :: heap + integer(psb_ipk_), intent(in) :: i + integer(psb_ipk_), intent(inout) :: nidx,info + integer(psb_ipk_), intent(inout) :: irwt(:) + real(psb_spk_), intent(in) :: thres,nrmi + integer(psb_ipk_), allocatable, intent(inout) :: idxs(:) + integer(psb_ipk_), intent(in) :: ja(:),irp(:) + real(psb_spk_), intent(in) :: val(:) + real(psb_spk_), intent(inout) :: row(:) + end subroutine + end interface + contains end module psb_s_invt_fact_mod diff --git a/prec/psb_z_invt_fact_mod.f90 b/prec/psb_z_invt_fact_mod.f90 index c0b3d4a3..1cdf32f4 100644 --- a/prec/psb_z_invt_fact_mod.f90 +++ b/prec/psb_z_invt_fact_mod.f90 @@ -145,6 +145,24 @@ module psb_z_invt_fact_mod end subroutine psb_z_invt_copyout end interface + interface psb_invt_inv + subroutine psb_z_invt_inv(thres,i,nrmi,row,heap,irwt,ja,irp,val, & + & nidx,idxs,info) + ! import + import psb_dpk_, psb_i_heap, psb_ipk_ + ! Arguments + type(psb_i_heap), intent(inout) :: heap + integer(psb_ipk_), intent(in) :: i + integer(psb_ipk_), intent(inout) :: nidx,info + integer(psb_ipk_), intent(inout) :: irwt(:) + real(psb_dpk_), intent(in) :: thres,nrmi + integer(psb_ipk_), allocatable, intent(inout) :: idxs(:) + integer(psb_ipk_), intent(in) :: ja(:),irp(:) + complex(psb_dpk_), intent(in) :: val(:) + complex(psb_dpk_), intent(inout) :: row(:) + end subroutine + end interface + contains end module psb_z_invt_fact_mod