Define a SHIFT argument to compute ILU( A+shft I)

maint-3.8.1
Salvatore Filippone 1 year ago
parent 3e29e603d2
commit d3fcd566d9

@ -130,7 +130,7 @@
! greater than 0. If the overlap is 0 or the matrix has been reordered ! greater than 0. If the overlap is 0 or the matrix has been reordered
! (see psb_fact_bld), then blck is empty. ! (see psb_fact_bld), then blck is empty.
! !
subroutine psb_cilu0_fact(ialg,a,l,u,d,info,blck, upd) subroutine psb_cilu0_fact(ialg,a,l,u,d,info,blck, upd,shft)
use psb_base_mod use psb_base_mod
use psb_c_ilu_fact_mod, psb_protect_name => psb_cilu0_fact use psb_c_ilu_fact_mod, psb_protect_name => psb_cilu0_fact
@ -145,11 +145,13 @@ subroutine psb_cilu0_fact(ialg,a,l,u,d,info,blck, upd)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
type(psb_cspmat_type),intent(in), optional, target :: blck type(psb_cspmat_type),intent(in), optional, target :: blck
character, intent(in), optional :: upd character, intent(in), optional :: upd
complex(psb_spk_), intent(in), optional :: shft
! Local variables ! Local variables
integer(psb_ipk_) :: l1, l2, m, err_act integer(psb_ipk_) :: l1, l2, m, err_act
type(psb_cspmat_type), pointer :: blck_ type(psb_cspmat_type), pointer :: blck_
type(psb_c_csr_sparse_mat) :: ll, uu type(psb_c_csr_sparse_mat) :: ll, uu
complex(psb_spk_) :: shft_
character :: upd_ character :: upd_
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
@ -177,6 +179,11 @@ subroutine psb_cilu0_fact(ialg,a,l,u,d,info,blck, upd)
else else
upd_ = 'F' upd_ = 'F'
end if end if
if (present(shft)) then
shft_ = shft
else
shft_ = czero
end if
m = a%get_nrows() + blck_%get_nrows() m = a%get_nrows() + blck_%get_nrows()
if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.& if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.&
@ -193,7 +200,7 @@ subroutine psb_cilu0_fact(ialg,a,l,u,d,info,blck, upd)
! Compute the ILU(0) or the MILU(0) factorization, depending on ialg ! Compute the ILU(0) or the MILU(0) factorization, depending on ialg
! !
call psb_cilu0_factint(ialg,a,blck_,& call psb_cilu0_factint(ialg,a,blck_,&
& d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,upd_,info) & d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,upd_,shft_,info)
if(info.ne.0) then if(info.ne.0) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_cilu0_factint' ch_err='psb_cilu0_factint'
@ -314,7 +321,7 @@ contains
! Error code. ! Error code.
! !
subroutine psb_cilu0_factint(ialg,a,b,& subroutine psb_cilu0_factint(ialg,a,b,&
& d,lval,lja,lirp,uval,uja,uirp,l1,l2,upd,info) & d,lval,lja,lirp,uval,uja,uirp,l1,l2,upd,shft,info)
implicit none implicit none
@ -325,6 +332,7 @@ contains
integer(psb_ipk_), intent(inout) :: lja(:),lirp(:),uja(:),uirp(:) integer(psb_ipk_), intent(inout) :: lja(:),lirp(:),uja(:),uirp(:)
complex(psb_spk_), intent(inout) :: lval(:),uval(:),d(:) complex(psb_spk_), intent(inout) :: lval(:),uval(:),d(:)
character, intent(in) :: upd character, intent(in) :: upd
complex(psb_spk_), intent(in) :: shft
! Local variables ! Local variables
integer(psb_ipk_) :: i,j,k,l,low1,low2,kk,jj,ll, ktrw,err_act, m integer(psb_ipk_) :: i,j,k,l,low1,low2,kk,jj,ll, ktrw,err_act, m
@ -382,14 +390,14 @@ contains
! into lval/d(i)/uval ! into lval/d(i)/uval
! !
call ilu_copyin(i,ma,a,i,ione,m,l1,lja,lval,& call ilu_copyin(i,ma,a,i,ione,m,l1,lja,lval,&
& d(i),l2,uja,uval,ktrw,trw,upd) & d(i),l2,uja,uval,ktrw,trw,upd,shft)
else else
! !
! Copy the i-th local row of the matrix, stored in b ! Copy the i-th local row of the matrix, stored in b
! (as (i-ma)-th row), into lval/d(i)/uval ! (as (i-ma)-th row), into lval/d(i)/uval
! !
call ilu_copyin(i-ma,mb,b,i,ione,m,l1,lja,lval,& call ilu_copyin(i-ma,mb,b,i,ione,m,l1,lja,lval,&
& d(i),l2,uja,uval,ktrw,trw,upd) & d(i),l2,uja,uval,ktrw,trw,upd,shft)
endif endif
lirp(i+1) = l1 + 1 lirp(i+1) = l1 + 1
@ -583,7 +591,7 @@ contains
! every nrb calls to copyin. If A is in CSR format it is unused. ! every nrb calls to copyin. If A is in CSR format it is unused.
! !
subroutine ilu_copyin(i,m,a,jd,jmin,jmax,l1,lja,lval,& subroutine ilu_copyin(i,m,a,jd,jmin,jmax,l1,lja,lval,&
& dia,l2,uja,uval,ktrw,trw,upd) & dia,l2,uja,uval,ktrw,trw,upd,shft)
use psb_base_mod use psb_base_mod
@ -597,6 +605,7 @@ contains
integer(psb_ipk_), intent(inout) :: lja(:), uja(:) integer(psb_ipk_), intent(inout) :: lja(:), uja(:)
complex(psb_spk_), intent(inout) :: lval(:), uval(:), dia complex(psb_spk_), intent(inout) :: lval(:), uval(:), dia
character, intent(in) :: upd character, intent(in) :: upd
complex(psb_spk_), intent(in) :: shft
! Local variables ! Local variables
integer(psb_ipk_) :: k,j,info,irb, nz integer(psb_ipk_) :: k,j,info,irb, nz
integer(psb_ipk_), parameter :: nrb=40 integer(psb_ipk_), parameter :: nrb=40
@ -625,7 +634,7 @@ contains
lval(l1) = aa%val(j) lval(l1) = aa%val(j)
lja(l1) = k lja(l1) = k
else if (k == jd) then else if (k == jd) then
dia = aa%val(j) dia = aa%val(j) + shft
else if ((k > jd).and.(k <= jmax)) then else if ((k > jd).and.(k <= jmax)) then
l2 = l2 + 1 l2 = l2 + 1
uval(l2) = aa%val(j) uval(l2) = aa%val(j)
@ -665,7 +674,7 @@ contains
lval(l1) = trw%val(ktrw) lval(l1) = trw%val(ktrw)
lja(l1) = k lja(l1) = k
else if (k == jd) then else if (k == jd) then
dia = trw%val(ktrw) dia = trw%val(ktrw) + shft
else if ((k > jd).and.(k <= jmax)) then else if ((k > jd).and.(k <= jmax)) then
l2 = l2 + 1 l2 = l2 + 1
uval(l2) = trw%val(ktrw) uval(l2) = trw%val(ktrw)

@ -127,7 +127,7 @@
! greater than 0. If the overlap is 0 or the matrix has been reordered ! greater than 0. If the overlap is 0 or the matrix has been reordered
! (see psb_fact_bld), then blck does not contain any row. ! (see psb_fact_bld), then blck does not contain any row.
! !
subroutine psb_ciluk_fact(fill_in,ialg,a,l,u,d,info,blck) subroutine psb_ciluk_fact(fill_in,ialg,a,l,u,d,info,blck,shft)
use psb_base_mod use psb_base_mod
use psb_c_ilu_fact_mod, psb_protect_name => psb_ciluk_fact use psb_c_ilu_fact_mod, psb_protect_name => psb_ciluk_fact
@ -141,6 +141,7 @@ subroutine psb_ciluk_fact(fill_in,ialg,a,l,u,d,info,blck)
type(psb_cspmat_type),intent(inout) :: l,u type(psb_cspmat_type),intent(inout) :: l,u
type(psb_cspmat_type),intent(in), optional, target :: blck type(psb_cspmat_type),intent(in), optional, target :: blck
complex(psb_spk_), intent(inout) :: d(:) complex(psb_spk_), intent(inout) :: d(:)
complex(psb_spk_), intent(in), optional :: shft
! Local Variables ! Local Variables
integer(psb_ipk_) :: l1, l2, m, err_act integer(psb_ipk_) :: l1, l2, m, err_act
@ -184,7 +185,7 @@ subroutine psb_ciluk_fact(fill_in,ialg,a,l,u,d,info,blck)
! Compute the ILU(k) or the MILU(k) factorization, depending on ialg ! Compute the ILU(k) or the MILU(k) factorization, depending on ialg
! !
call psb_ciluk_factint(fill_in,ialg,a,blck_,& call psb_ciluk_factint(fill_in,ialg,a,blck_,&
& d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,info) & d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,info,shft)
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_ciluk_factint' ch_err='psb_ciluk_factint'
@ -298,7 +299,7 @@ contains
! Error code. ! Error code.
! !
subroutine psb_ciluk_factint(fill_in,ialg,a,b,& subroutine psb_ciluk_factint(fill_in,ialg,a,b,&
& d,lval,lja,lirp,uval,uja,uirp,l1,l2,info) & d,lval,lja,lirp,uval,uja,uirp,l1,l2,info,shft)
use psb_base_mod use psb_base_mod
@ -311,6 +312,7 @@ contains
integer(psb_ipk_), allocatable, intent(inout) :: lja(:),lirp(:),uja(:),uirp(:) integer(psb_ipk_), allocatable, intent(inout) :: lja(:),lirp(:),uja(:),uirp(:)
complex(psb_spk_), allocatable, intent(inout) :: lval(:),uval(:) complex(psb_spk_), allocatable, intent(inout) :: lval(:),uval(:)
complex(psb_spk_), intent(inout) :: d(:) complex(psb_spk_), intent(inout) :: d(:)
complex(psb_spk_), intent(in) :: shft
! Local variables ! Local variables
integer(psb_ipk_) :: ma,mb,i, ktrw,err_act,nidx, m integer(psb_ipk_) :: ma,mb,i, ktrw,err_act,nidx, m
@ -400,13 +402,13 @@ contains
! !
! Copy into trw the i-th local row of the matrix, stored in a ! Copy into trw the i-th local row of the matrix, stored in a
! !
call iluk_copyin(i,ma,a,ione,m,row,rowlevs,heap,ktrw,trw,info) call iluk_copyin(i,ma,a,ione,m,row,rowlevs,heap,ktrw,trw,info,shft)
else else
! !
! Copy into trw the i-th local row of the matrix, stored in b ! Copy into trw the i-th local row of the matrix, stored in b
! (as (i-ma)-th row) ! (as (i-ma)-th row)
! !
call iluk_copyin(i-ma,mb,b,ione,m,row,rowlevs,heap,ktrw,trw,info) call iluk_copyin(i-ma,mb,b,ione,m,row,rowlevs,heap,ktrw,trw,info,shft)
endif endif
! Do an elimination step on the current row. It turns out we only ! Do an elimination step on the current row. It turns out we only
@ -516,7 +518,7 @@ contains
! until we empty the buffer. Thus we will make a call to psb_sp_getblk ! until we empty the buffer. Thus we will make a call to psb_sp_getblk
! every nrb calls to copyin. If A is in CSR format it is unused. ! every nrb calls to copyin. If A is in CSR format it is unused.
! !
subroutine iluk_copyin(i,m,a,jmin,jmax,row,rowlevs,heap,ktrw,trw,info) subroutine iluk_copyin(i,m,a,jmin,jmax,row,rowlevs,heap,ktrw,trw,info,shft)
use psb_base_mod use psb_base_mod
@ -530,6 +532,8 @@ contains
integer(psb_ipk_), intent(inout) :: rowlevs(:) integer(psb_ipk_), intent(inout) :: rowlevs(:)
complex(psb_spk_), intent(inout) :: row(:) complex(psb_spk_), intent(inout) :: row(:)
type(psb_i_heap), intent(inout) :: heap type(psb_i_heap), intent(inout) :: heap
complex(psb_spk_), intent(in) :: shft
! Local variables ! Local variables
integer(psb_ipk_) :: k,j,irb,err_act,nz integer(psb_ipk_) :: k,j,irb,err_act,nz
@ -554,6 +558,7 @@ contains
k = aa%ja(j) k = aa%ja(j)
if ((jmin<=k).and.(k<=jmax)) then if ((jmin<=k).and.(k<=jmax)) then
row(k) = aa%val(j) row(k) = aa%val(j)
if (k==i) row(k) = row(k) + shft
rowlevs(k) = 0 rowlevs(k) = 0
call heap%insert(k,info) call heap%insert(k,info)
end if end if
@ -587,6 +592,7 @@ contains
k = trw%ja(ktrw) k = trw%ja(ktrw)
if ((jmin<=k).and.(k<=jmax)) then if ((jmin<=k).and.(k<=jmax)) then
row(k) = trw%val(ktrw) row(k) = trw%val(ktrw)
if (k==i) row(k) = row(k) + shft
rowlevs(k) = 0 rowlevs(k) = 0
call heap%insert(k,info) call heap%insert(k,info)
end if end if
@ -670,7 +676,8 @@ contains
! Note: this argument is intent(inout) and not only intent(out) ! Note: this argument is intent(inout) and not only intent(out)
! to retain its allocation, done by this routine. ! to retain its allocation, done by this routine.
! !
subroutine iluk_fact(fill_in,i,row,rowlevs,heap,d,uja,uirp,uval,uplevs,nidx,idxs,info) subroutine iluk_fact(fill_in,i,row,rowlevs,heap,d,&
& uja,uirp,uval,uplevs,nidx,idxs,info)
use psb_base_mod use psb_base_mod

@ -123,7 +123,7 @@
! greater than 0. If the overlap is 0 or the matrix has been reordered ! greater than 0. If the overlap is 0 or the matrix has been reordered
! (see psb_fact_bld), then blck does not contain any row. ! (see psb_fact_bld), then blck does not contain any row.
! !
subroutine psb_cilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) subroutine psb_cilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale,shft)
use psb_base_mod use psb_base_mod
use psb_c_ilu_fact_mod, psb_protect_name => psb_cilut_fact use psb_c_ilu_fact_mod, psb_protect_name => psb_cilut_fact
@ -139,6 +139,7 @@ subroutine psb_cilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale)
complex(psb_spk_), intent(inout) :: d(:) complex(psb_spk_), intent(inout) :: d(:)
type(psb_cspmat_type),intent(in), optional, target :: blck type(psb_cspmat_type),intent(in), optional, target :: blck
integer(psb_ipk_), intent(in), optional :: iscale integer(psb_ipk_), intent(in), optional :: iscale
complex(psb_spk_), intent(in), optional :: shft
! Local Variables ! Local Variables
integer(psb_ipk_) :: l1, l2, m, err_act, iscale_ integer(psb_ipk_) :: l1, l2, m, err_act, iscale_
@ -206,7 +207,7 @@ subroutine psb_cilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale)
! Compute the ILU(k,t) factorization ! Compute the ILU(k,t) factorization
! !
call psb_cilut_factint(fill_in,thres,a,blck_,& call psb_cilut_factint(fill_in,thres,a,blck_,&
& d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,info,scale) & d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,info,scale,shft)
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_cilut_factint' ch_err='psb_cilut_factint'
@ -316,7 +317,7 @@ contains
! Error code. ! Error code.
! !
subroutine psb_cilut_factint(fill_in,thres,a,b,& subroutine psb_cilut_factint(fill_in,thres,a,b,&
& d,lval,lja,lirp,uval,uja,uirp,l1,l2,info,scale) & d,lval,lja,lirp,uval,uja,uirp,l1,l2,info,scale,shft)
use psb_base_mod use psb_base_mod
@ -331,6 +332,7 @@ contains
complex(psb_spk_), allocatable, intent(inout) :: lval(:),uval(:) complex(psb_spk_), allocatable, intent(inout) :: lval(:),uval(:)
complex(psb_spk_), intent(inout) :: d(:) complex(psb_spk_), intent(inout) :: d(:)
real(psb_spk_), intent(in), optional :: scale real(psb_spk_), intent(in), optional :: scale
complex(psb_spk_), intent(in) :: shft
! Local Variables ! Local Variables
integer(psb_ipk_) :: i, ktrw,err_act,nidx,nlw,nup,jmaxup, ma, mb, m integer(psb_ipk_) :: i, ktrw,err_act,nidx,nlw,nup,jmaxup, ma, mb, m
@ -401,10 +403,10 @@ contains
d(i) = czero d(i) = czero
if (i<=ma) then if (i<=ma) then
call ilut_copyin(i,ma,a,i,ione,m,nlw,nup,jmaxup,nrmi,weight,& call ilut_copyin(i,ma,a,i,ione,m,nlw,nup,jmaxup,nrmi,weight,&
& row,heap,ktrw,trw,info) & row,heap,ktrw,trw,info,shft)
else else
call ilut_copyin(i-ma,mb,b,i,ione,m,nlw,nup,jmaxup,nrmi,weight,& call ilut_copyin(i-ma,mb,b,i,ione,m,nlw,nup,jmaxup,nrmi,weight,&
& row,heap,ktrw,trw,info) & row,heap,ktrw,trw,info,shft)
endif endif
! !
@ -540,7 +542,7 @@ contains
! every nrb calls to copyin. If A is in CSR format it is unused. ! every nrb calls to copyin. If A is in CSR format it is unused.
! !
subroutine ilut_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,& subroutine ilut_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,&
& nrmi,weight,row,heap,ktrw,trw,info) & nrmi,weight,row,heap,ktrw,trw,info,shft)
use psb_base_mod use psb_base_mod
implicit none implicit none
type(psb_cspmat_type), intent(in) :: a type(psb_cspmat_type), intent(in) :: a
@ -551,6 +553,7 @@ contains
complex(psb_spk_), intent(inout) :: row(:) complex(psb_spk_), intent(inout) :: row(:)
real(psb_spk_), intent(in) :: weight real(psb_spk_), intent(in) :: weight
type(psb_i_heap), intent(inout) :: heap type(psb_i_heap), intent(inout) :: heap
complex(psb_spk_), intent(in) :: shft
integer(psb_ipk_) :: k,j,irb,kin,nz integer(psb_ipk_) :: k,j,irb,kin,nz
integer(psb_ipk_), parameter :: nrb=40 integer(psb_ipk_), parameter :: nrb=40
@ -597,6 +600,7 @@ contains
call heap%insert(k,info) call heap%insert(k,info)
if (info /= psb_success_) exit if (info /= psb_success_) exit
if (k<jd) nlw = nlw + 1 if (k<jd) nlw = nlw + 1
if (k == jd) row(k) = row(k) + (shft*weight)
if (k>jd) then if (k>jd) then
nup = nup + 1 nup = nup + 1
if (abs(row(k))>dmaxup) then if (abs(row(k))>dmaxup) then
@ -648,6 +652,7 @@ contains
call heap%insert(k,info) call heap%insert(k,info)
if (info /= psb_success_) exit if (info /= psb_success_) exit
if (k<jd) nlw = nlw + 1 if (k<jd) nlw = nlw + 1
if (k == jd) row(k) = row(k) + (shft*weight)
if (k>jd) then if (k>jd) then
nup = nup + 1 nup = nup + 1
if (abs(row(k))>dmaxup) then if (abs(row(k))>dmaxup) then

@ -130,7 +130,7 @@
! greater than 0. If the overlap is 0 or the matrix has been reordered ! greater than 0. If the overlap is 0 or the matrix has been reordered
! (see psb_fact_bld), then blck is empty. ! (see psb_fact_bld), then blck is empty.
! !
subroutine psb_dilu0_fact(ialg,a,l,u,d,info,blck, upd) subroutine psb_dilu0_fact(ialg,a,l,u,d,info,blck, upd,shft)
use psb_base_mod use psb_base_mod
use psb_d_ilu_fact_mod, psb_protect_name => psb_dilu0_fact use psb_d_ilu_fact_mod, psb_protect_name => psb_dilu0_fact
@ -145,11 +145,13 @@ subroutine psb_dilu0_fact(ialg,a,l,u,d,info,blck, upd)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
type(psb_dspmat_type),intent(in), optional, target :: blck type(psb_dspmat_type),intent(in), optional, target :: blck
character, intent(in), optional :: upd character, intent(in), optional :: upd
real(psb_dpk_), intent(in), optional :: shft
! Local variables ! Local variables
integer(psb_ipk_) :: l1, l2, m, err_act integer(psb_ipk_) :: l1, l2, m, err_act
type(psb_dspmat_type), pointer :: blck_ type(psb_dspmat_type), pointer :: blck_
type(psb_d_csr_sparse_mat) :: ll, uu type(psb_d_csr_sparse_mat) :: ll, uu
real(psb_dpk_) :: shft_
character :: upd_ character :: upd_
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
@ -177,6 +179,11 @@ subroutine psb_dilu0_fact(ialg,a,l,u,d,info,blck, upd)
else else
upd_ = 'F' upd_ = 'F'
end if end if
if (present(shft)) then
shft_ = shft
else
shft_ = dzero
end if
m = a%get_nrows() + blck_%get_nrows() m = a%get_nrows() + blck_%get_nrows()
if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.& if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.&
@ -193,7 +200,7 @@ subroutine psb_dilu0_fact(ialg,a,l,u,d,info,blck, upd)
! Compute the ILU(0) or the MILU(0) factorization, depending on ialg ! Compute the ILU(0) or the MILU(0) factorization, depending on ialg
! !
call psb_dilu0_factint(ialg,a,blck_,& call psb_dilu0_factint(ialg,a,blck_,&
& d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,upd_,info) & d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,upd_,shft_,info)
if(info.ne.0) then if(info.ne.0) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_dilu0_factint' ch_err='psb_dilu0_factint'
@ -314,7 +321,7 @@ contains
! Error code. ! Error code.
! !
subroutine psb_dilu0_factint(ialg,a,b,& subroutine psb_dilu0_factint(ialg,a,b,&
& d,lval,lja,lirp,uval,uja,uirp,l1,l2,upd,info) & d,lval,lja,lirp,uval,uja,uirp,l1,l2,upd,shft,info)
implicit none implicit none
@ -325,6 +332,7 @@ contains
integer(psb_ipk_), intent(inout) :: lja(:),lirp(:),uja(:),uirp(:) integer(psb_ipk_), intent(inout) :: lja(:),lirp(:),uja(:),uirp(:)
real(psb_dpk_), intent(inout) :: lval(:),uval(:),d(:) real(psb_dpk_), intent(inout) :: lval(:),uval(:),d(:)
character, intent(in) :: upd character, intent(in) :: upd
real(psb_dpk_), intent(in) :: shft
! Local variables ! Local variables
integer(psb_ipk_) :: i,j,k,l,low1,low2,kk,jj,ll, ktrw,err_act, m integer(psb_ipk_) :: i,j,k,l,low1,low2,kk,jj,ll, ktrw,err_act, m
@ -382,14 +390,14 @@ contains
! into lval/d(i)/uval ! into lval/d(i)/uval
! !
call ilu_copyin(i,ma,a,i,ione,m,l1,lja,lval,& call ilu_copyin(i,ma,a,i,ione,m,l1,lja,lval,&
& d(i),l2,uja,uval,ktrw,trw,upd) & d(i),l2,uja,uval,ktrw,trw,upd,shft)
else else
! !
! Copy the i-th local row of the matrix, stored in b ! Copy the i-th local row of the matrix, stored in b
! (as (i-ma)-th row), into lval/d(i)/uval ! (as (i-ma)-th row), into lval/d(i)/uval
! !
call ilu_copyin(i-ma,mb,b,i,ione,m,l1,lja,lval,& call ilu_copyin(i-ma,mb,b,i,ione,m,l1,lja,lval,&
& d(i),l2,uja,uval,ktrw,trw,upd) & d(i),l2,uja,uval,ktrw,trw,upd,shft)
endif endif
lirp(i+1) = l1 + 1 lirp(i+1) = l1 + 1
@ -583,7 +591,7 @@ contains
! every nrb calls to copyin. If A is in CSR format it is unused. ! every nrb calls to copyin. If A is in CSR format it is unused.
! !
subroutine ilu_copyin(i,m,a,jd,jmin,jmax,l1,lja,lval,& subroutine ilu_copyin(i,m,a,jd,jmin,jmax,l1,lja,lval,&
& dia,l2,uja,uval,ktrw,trw,upd) & dia,l2,uja,uval,ktrw,trw,upd,shft)
use psb_base_mod use psb_base_mod
@ -597,6 +605,7 @@ contains
integer(psb_ipk_), intent(inout) :: lja(:), uja(:) integer(psb_ipk_), intent(inout) :: lja(:), uja(:)
real(psb_dpk_), intent(inout) :: lval(:), uval(:), dia real(psb_dpk_), intent(inout) :: lval(:), uval(:), dia
character, intent(in) :: upd character, intent(in) :: upd
real(psb_dpk_), intent(in) :: shft
! Local variables ! Local variables
integer(psb_ipk_) :: k,j,info,irb, nz integer(psb_ipk_) :: k,j,info,irb, nz
integer(psb_ipk_), parameter :: nrb=40 integer(psb_ipk_), parameter :: nrb=40
@ -625,7 +634,7 @@ contains
lval(l1) = aa%val(j) lval(l1) = aa%val(j)
lja(l1) = k lja(l1) = k
else if (k == jd) then else if (k == jd) then
dia = aa%val(j) dia = aa%val(j) + shft
else if ((k > jd).and.(k <= jmax)) then else if ((k > jd).and.(k <= jmax)) then
l2 = l2 + 1 l2 = l2 + 1
uval(l2) = aa%val(j) uval(l2) = aa%val(j)
@ -665,7 +674,7 @@ contains
lval(l1) = trw%val(ktrw) lval(l1) = trw%val(ktrw)
lja(l1) = k lja(l1) = k
else if (k == jd) then else if (k == jd) then
dia = trw%val(ktrw) dia = trw%val(ktrw) + shft
else if ((k > jd).and.(k <= jmax)) then else if ((k > jd).and.(k <= jmax)) then
l2 = l2 + 1 l2 = l2 + 1
uval(l2) = trw%val(ktrw) uval(l2) = trw%val(ktrw)

@ -127,7 +127,7 @@
! greater than 0. If the overlap is 0 or the matrix has been reordered ! greater than 0. If the overlap is 0 or the matrix has been reordered
! (see psb_fact_bld), then blck does not contain any row. ! (see psb_fact_bld), then blck does not contain any row.
! !
subroutine psb_diluk_fact(fill_in,ialg,a,l,u,d,info,blck) subroutine psb_diluk_fact(fill_in,ialg,a,l,u,d,info,blck,shft)
use psb_base_mod use psb_base_mod
use psb_d_ilu_fact_mod, psb_protect_name => psb_diluk_fact use psb_d_ilu_fact_mod, psb_protect_name => psb_diluk_fact
@ -141,6 +141,7 @@ subroutine psb_diluk_fact(fill_in,ialg,a,l,u,d,info,blck)
type(psb_dspmat_type),intent(inout) :: l,u type(psb_dspmat_type),intent(inout) :: l,u
type(psb_dspmat_type),intent(in), optional, target :: blck type(psb_dspmat_type),intent(in), optional, target :: blck
real(psb_dpk_), intent(inout) :: d(:) real(psb_dpk_), intent(inout) :: d(:)
real(psb_dpk_), intent(in), optional :: shft
! Local Variables ! Local Variables
integer(psb_ipk_) :: l1, l2, m, err_act integer(psb_ipk_) :: l1, l2, m, err_act
@ -184,7 +185,7 @@ subroutine psb_diluk_fact(fill_in,ialg,a,l,u,d,info,blck)
! Compute the ILU(k) or the MILU(k) factorization, depending on ialg ! Compute the ILU(k) or the MILU(k) factorization, depending on ialg
! !
call psb_diluk_factint(fill_in,ialg,a,blck_,& call psb_diluk_factint(fill_in,ialg,a,blck_,&
& d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,info) & d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,info,shft)
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_diluk_factint' ch_err='psb_diluk_factint'
@ -298,7 +299,7 @@ contains
! Error code. ! Error code.
! !
subroutine psb_diluk_factint(fill_in,ialg,a,b,& subroutine psb_diluk_factint(fill_in,ialg,a,b,&
& d,lval,lja,lirp,uval,uja,uirp,l1,l2,info) & d,lval,lja,lirp,uval,uja,uirp,l1,l2,info,shft)
use psb_base_mod use psb_base_mod
@ -311,6 +312,7 @@ contains
integer(psb_ipk_), allocatable, intent(inout) :: lja(:),lirp(:),uja(:),uirp(:) integer(psb_ipk_), allocatable, intent(inout) :: lja(:),lirp(:),uja(:),uirp(:)
real(psb_dpk_), allocatable, intent(inout) :: lval(:),uval(:) real(psb_dpk_), allocatable, intent(inout) :: lval(:),uval(:)
real(psb_dpk_), intent(inout) :: d(:) real(psb_dpk_), intent(inout) :: d(:)
real(psb_dpk_), intent(in) :: shft
! Local variables ! Local variables
integer(psb_ipk_) :: ma,mb,i, ktrw,err_act,nidx, m integer(psb_ipk_) :: ma,mb,i, ktrw,err_act,nidx, m
@ -400,13 +402,13 @@ contains
! !
! Copy into trw the i-th local row of the matrix, stored in a ! Copy into trw the i-th local row of the matrix, stored in a
! !
call iluk_copyin(i,ma,a,ione,m,row,rowlevs,heap,ktrw,trw,info) call iluk_copyin(i,ma,a,ione,m,row,rowlevs,heap,ktrw,trw,info,shft)
else else
! !
! Copy into trw the i-th local row of the matrix, stored in b ! Copy into trw the i-th local row of the matrix, stored in b
! (as (i-ma)-th row) ! (as (i-ma)-th row)
! !
call iluk_copyin(i-ma,mb,b,ione,m,row,rowlevs,heap,ktrw,trw,info) call iluk_copyin(i-ma,mb,b,ione,m,row,rowlevs,heap,ktrw,trw,info,shft)
endif endif
! Do an elimination step on the current row. It turns out we only ! Do an elimination step on the current row. It turns out we only
@ -516,7 +518,7 @@ contains
! until we empty the buffer. Thus we will make a call to psb_sp_getblk ! until we empty the buffer. Thus we will make a call to psb_sp_getblk
! every nrb calls to copyin. If A is in CSR format it is unused. ! every nrb calls to copyin. If A is in CSR format it is unused.
! !
subroutine iluk_copyin(i,m,a,jmin,jmax,row,rowlevs,heap,ktrw,trw,info) subroutine iluk_copyin(i,m,a,jmin,jmax,row,rowlevs,heap,ktrw,trw,info,shft)
use psb_base_mod use psb_base_mod
@ -530,6 +532,8 @@ contains
integer(psb_ipk_), intent(inout) :: rowlevs(:) integer(psb_ipk_), intent(inout) :: rowlevs(:)
real(psb_dpk_), intent(inout) :: row(:) real(psb_dpk_), intent(inout) :: row(:)
type(psb_i_heap), intent(inout) :: heap type(psb_i_heap), intent(inout) :: heap
real(psb_dpk_), intent(in) :: shft
! Local variables ! Local variables
integer(psb_ipk_) :: k,j,irb,err_act,nz integer(psb_ipk_) :: k,j,irb,err_act,nz
@ -554,6 +558,7 @@ contains
k = aa%ja(j) k = aa%ja(j)
if ((jmin<=k).and.(k<=jmax)) then if ((jmin<=k).and.(k<=jmax)) then
row(k) = aa%val(j) row(k) = aa%val(j)
if (k==i) row(k) = row(k) + shft
rowlevs(k) = 0 rowlevs(k) = 0
call heap%insert(k,info) call heap%insert(k,info)
end if end if
@ -587,6 +592,7 @@ contains
k = trw%ja(ktrw) k = trw%ja(ktrw)
if ((jmin<=k).and.(k<=jmax)) then if ((jmin<=k).and.(k<=jmax)) then
row(k) = trw%val(ktrw) row(k) = trw%val(ktrw)
if (k==i) row(k) = row(k) + shft
rowlevs(k) = 0 rowlevs(k) = 0
call heap%insert(k,info) call heap%insert(k,info)
end if end if
@ -670,7 +676,8 @@ contains
! Note: this argument is intent(inout) and not only intent(out) ! Note: this argument is intent(inout) and not only intent(out)
! to retain its allocation, done by this routine. ! to retain its allocation, done by this routine.
! !
subroutine iluk_fact(fill_in,i,row,rowlevs,heap,d,uja,uirp,uval,uplevs,nidx,idxs,info) subroutine iluk_fact(fill_in,i,row,rowlevs,heap,d,&
& uja,uirp,uval,uplevs,nidx,idxs,info)
use psb_base_mod use psb_base_mod

@ -123,7 +123,7 @@
! greater than 0. If the overlap is 0 or the matrix has been reordered ! greater than 0. If the overlap is 0 or the matrix has been reordered
! (see psb_fact_bld), then blck does not contain any row. ! (see psb_fact_bld), then blck does not contain any row.
! !
subroutine psb_dilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) subroutine psb_dilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale,shft)
use psb_base_mod use psb_base_mod
use psb_d_ilu_fact_mod, psb_protect_name => psb_dilut_fact use psb_d_ilu_fact_mod, psb_protect_name => psb_dilut_fact
@ -139,6 +139,7 @@ subroutine psb_dilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale)
real(psb_dpk_), intent(inout) :: d(:) real(psb_dpk_), intent(inout) :: d(:)
type(psb_dspmat_type),intent(in), optional, target :: blck type(psb_dspmat_type),intent(in), optional, target :: blck
integer(psb_ipk_), intent(in), optional :: iscale integer(psb_ipk_), intent(in), optional :: iscale
real(psb_dpk_), intent(in), optional :: shft
! Local Variables ! Local Variables
integer(psb_ipk_) :: l1, l2, m, err_act, iscale_ integer(psb_ipk_) :: l1, l2, m, err_act, iscale_
@ -206,7 +207,7 @@ subroutine psb_dilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale)
! Compute the ILU(k,t) factorization ! Compute the ILU(k,t) factorization
! !
call psb_dilut_factint(fill_in,thres,a,blck_,& call psb_dilut_factint(fill_in,thres,a,blck_,&
& d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,info,scale) & d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,info,scale,shft)
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_dilut_factint' ch_err='psb_dilut_factint'
@ -316,7 +317,7 @@ contains
! Error code. ! Error code.
! !
subroutine psb_dilut_factint(fill_in,thres,a,b,& subroutine psb_dilut_factint(fill_in,thres,a,b,&
& d,lval,lja,lirp,uval,uja,uirp,l1,l2,info,scale) & d,lval,lja,lirp,uval,uja,uirp,l1,l2,info,scale,shft)
use psb_base_mod use psb_base_mod
@ -331,6 +332,7 @@ contains
real(psb_dpk_), allocatable, intent(inout) :: lval(:),uval(:) real(psb_dpk_), allocatable, intent(inout) :: lval(:),uval(:)
real(psb_dpk_), intent(inout) :: d(:) real(psb_dpk_), intent(inout) :: d(:)
real(psb_dpk_), intent(in), optional :: scale real(psb_dpk_), intent(in), optional :: scale
real(psb_dpk_), intent(in) :: shft
! Local Variables ! Local Variables
integer(psb_ipk_) :: i, ktrw,err_act,nidx,nlw,nup,jmaxup, ma, mb, m integer(psb_ipk_) :: i, ktrw,err_act,nidx,nlw,nup,jmaxup, ma, mb, m
@ -401,10 +403,10 @@ contains
d(i) = czero d(i) = czero
if (i<=ma) then if (i<=ma) then
call ilut_copyin(i,ma,a,i,ione,m,nlw,nup,jmaxup,nrmi,weight,& call ilut_copyin(i,ma,a,i,ione,m,nlw,nup,jmaxup,nrmi,weight,&
& row,heap,ktrw,trw,info) & row,heap,ktrw,trw,info,shft)
else else
call ilut_copyin(i-ma,mb,b,i,ione,m,nlw,nup,jmaxup,nrmi,weight,& call ilut_copyin(i-ma,mb,b,i,ione,m,nlw,nup,jmaxup,nrmi,weight,&
& row,heap,ktrw,trw,info) & row,heap,ktrw,trw,info,shft)
endif endif
! !
@ -540,7 +542,7 @@ contains
! every nrb calls to copyin. If A is in CSR format it is unused. ! every nrb calls to copyin. If A is in CSR format it is unused.
! !
subroutine ilut_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,& subroutine ilut_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,&
& nrmi,weight,row,heap,ktrw,trw,info) & nrmi,weight,row,heap,ktrw,trw,info,shft)
use psb_base_mod use psb_base_mod
implicit none implicit none
type(psb_dspmat_type), intent(in) :: a type(psb_dspmat_type), intent(in) :: a
@ -551,6 +553,7 @@ contains
real(psb_dpk_), intent(inout) :: row(:) real(psb_dpk_), intent(inout) :: row(:)
real(psb_dpk_), intent(in) :: weight real(psb_dpk_), intent(in) :: weight
type(psb_i_heap), intent(inout) :: heap type(psb_i_heap), intent(inout) :: heap
real(psb_dpk_), intent(in) :: shft
integer(psb_ipk_) :: k,j,irb,kin,nz integer(psb_ipk_) :: k,j,irb,kin,nz
integer(psb_ipk_), parameter :: nrb=40 integer(psb_ipk_), parameter :: nrb=40
@ -597,6 +600,7 @@ contains
call heap%insert(k,info) call heap%insert(k,info)
if (info /= psb_success_) exit if (info /= psb_success_) exit
if (k<jd) nlw = nlw + 1 if (k<jd) nlw = nlw + 1
if (k == jd) row(k) = row(k) + (shft*weight)
if (k>jd) then if (k>jd) then
nup = nup + 1 nup = nup + 1
if (abs(row(k))>dmaxup) then if (abs(row(k))>dmaxup) then
@ -648,6 +652,7 @@ contains
call heap%insert(k,info) call heap%insert(k,info)
if (info /= psb_success_) exit if (info /= psb_success_) exit
if (k<jd) nlw = nlw + 1 if (k<jd) nlw = nlw + 1
if (k == jd) row(k) = row(k) + (shft*weight)
if (k>jd) then if (k>jd) then
nup = nup + 1 nup = nup + 1
if (abs(row(k))>dmaxup) then if (abs(row(k))>dmaxup) then

@ -130,7 +130,7 @@
! greater than 0. If the overlap is 0 or the matrix has been reordered ! greater than 0. If the overlap is 0 or the matrix has been reordered
! (see psb_fact_bld), then blck is empty. ! (see psb_fact_bld), then blck is empty.
! !
subroutine psb_silu0_fact(ialg,a,l,u,d,info,blck, upd) subroutine psb_silu0_fact(ialg,a,l,u,d,info,blck, upd,shft)
use psb_base_mod use psb_base_mod
use psb_s_ilu_fact_mod, psb_protect_name => psb_silu0_fact use psb_s_ilu_fact_mod, psb_protect_name => psb_silu0_fact
@ -145,11 +145,13 @@ subroutine psb_silu0_fact(ialg,a,l,u,d,info,blck, upd)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
type(psb_sspmat_type),intent(in), optional, target :: blck type(psb_sspmat_type),intent(in), optional, target :: blck
character, intent(in), optional :: upd character, intent(in), optional :: upd
real(psb_spk_), intent(in), optional :: shft
! Local variables ! Local variables
integer(psb_ipk_) :: l1, l2, m, err_act integer(psb_ipk_) :: l1, l2, m, err_act
type(psb_sspmat_type), pointer :: blck_ type(psb_sspmat_type), pointer :: blck_
type(psb_s_csr_sparse_mat) :: ll, uu type(psb_s_csr_sparse_mat) :: ll, uu
real(psb_spk_) :: shft_
character :: upd_ character :: upd_
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
@ -177,6 +179,11 @@ subroutine psb_silu0_fact(ialg,a,l,u,d,info,blck, upd)
else else
upd_ = 'F' upd_ = 'F'
end if end if
if (present(shft)) then
shft_ = shft
else
shft_ = szero
end if
m = a%get_nrows() + blck_%get_nrows() m = a%get_nrows() + blck_%get_nrows()
if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.& if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.&
@ -193,7 +200,7 @@ subroutine psb_silu0_fact(ialg,a,l,u,d,info,blck, upd)
! Compute the ILU(0) or the MILU(0) factorization, depending on ialg ! Compute the ILU(0) or the MILU(0) factorization, depending on ialg
! !
call psb_silu0_factint(ialg,a,blck_,& call psb_silu0_factint(ialg,a,blck_,&
& d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,upd_,info) & d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,upd_,shft_,info)
if(info.ne.0) then if(info.ne.0) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_silu0_factint' ch_err='psb_silu0_factint'
@ -314,7 +321,7 @@ contains
! Error code. ! Error code.
! !
subroutine psb_silu0_factint(ialg,a,b,& subroutine psb_silu0_factint(ialg,a,b,&
& d,lval,lja,lirp,uval,uja,uirp,l1,l2,upd,info) & d,lval,lja,lirp,uval,uja,uirp,l1,l2,upd,shft,info)
implicit none implicit none
@ -325,6 +332,7 @@ contains
integer(psb_ipk_), intent(inout) :: lja(:),lirp(:),uja(:),uirp(:) integer(psb_ipk_), intent(inout) :: lja(:),lirp(:),uja(:),uirp(:)
real(psb_spk_), intent(inout) :: lval(:),uval(:),d(:) real(psb_spk_), intent(inout) :: lval(:),uval(:),d(:)
character, intent(in) :: upd character, intent(in) :: upd
real(psb_spk_), intent(in) :: shft
! Local variables ! Local variables
integer(psb_ipk_) :: i,j,k,l,low1,low2,kk,jj,ll, ktrw,err_act, m integer(psb_ipk_) :: i,j,k,l,low1,low2,kk,jj,ll, ktrw,err_act, m
@ -382,14 +390,14 @@ contains
! into lval/d(i)/uval ! into lval/d(i)/uval
! !
call ilu_copyin(i,ma,a,i,ione,m,l1,lja,lval,& call ilu_copyin(i,ma,a,i,ione,m,l1,lja,lval,&
& d(i),l2,uja,uval,ktrw,trw,upd) & d(i),l2,uja,uval,ktrw,trw,upd,shft)
else else
! !
! Copy the i-th local row of the matrix, stored in b ! Copy the i-th local row of the matrix, stored in b
! (as (i-ma)-th row), into lval/d(i)/uval ! (as (i-ma)-th row), into lval/d(i)/uval
! !
call ilu_copyin(i-ma,mb,b,i,ione,m,l1,lja,lval,& call ilu_copyin(i-ma,mb,b,i,ione,m,l1,lja,lval,&
& d(i),l2,uja,uval,ktrw,trw,upd) & d(i),l2,uja,uval,ktrw,trw,upd,shft)
endif endif
lirp(i+1) = l1 + 1 lirp(i+1) = l1 + 1
@ -583,7 +591,7 @@ contains
! every nrb calls to copyin. If A is in CSR format it is unused. ! every nrb calls to copyin. If A is in CSR format it is unused.
! !
subroutine ilu_copyin(i,m,a,jd,jmin,jmax,l1,lja,lval,& subroutine ilu_copyin(i,m,a,jd,jmin,jmax,l1,lja,lval,&
& dia,l2,uja,uval,ktrw,trw,upd) & dia,l2,uja,uval,ktrw,trw,upd,shft)
use psb_base_mod use psb_base_mod
@ -597,6 +605,7 @@ contains
integer(psb_ipk_), intent(inout) :: lja(:), uja(:) integer(psb_ipk_), intent(inout) :: lja(:), uja(:)
real(psb_spk_), intent(inout) :: lval(:), uval(:), dia real(psb_spk_), intent(inout) :: lval(:), uval(:), dia
character, intent(in) :: upd character, intent(in) :: upd
real(psb_spk_), intent(in) :: shft
! Local variables ! Local variables
integer(psb_ipk_) :: k,j,info,irb, nz integer(psb_ipk_) :: k,j,info,irb, nz
integer(psb_ipk_), parameter :: nrb=40 integer(psb_ipk_), parameter :: nrb=40
@ -625,7 +634,7 @@ contains
lval(l1) = aa%val(j) lval(l1) = aa%val(j)
lja(l1) = k lja(l1) = k
else if (k == jd) then else if (k == jd) then
dia = aa%val(j) dia = aa%val(j) + shft
else if ((k > jd).and.(k <= jmax)) then else if ((k > jd).and.(k <= jmax)) then
l2 = l2 + 1 l2 = l2 + 1
uval(l2) = aa%val(j) uval(l2) = aa%val(j)
@ -665,7 +674,7 @@ contains
lval(l1) = trw%val(ktrw) lval(l1) = trw%val(ktrw)
lja(l1) = k lja(l1) = k
else if (k == jd) then else if (k == jd) then
dia = trw%val(ktrw) dia = trw%val(ktrw) + shft
else if ((k > jd).and.(k <= jmax)) then else if ((k > jd).and.(k <= jmax)) then
l2 = l2 + 1 l2 = l2 + 1
uval(l2) = trw%val(ktrw) uval(l2) = trw%val(ktrw)

@ -127,7 +127,7 @@
! greater than 0. If the overlap is 0 or the matrix has been reordered ! greater than 0. If the overlap is 0 or the matrix has been reordered
! (see psb_fact_bld), then blck does not contain any row. ! (see psb_fact_bld), then blck does not contain any row.
! !
subroutine psb_siluk_fact(fill_in,ialg,a,l,u,d,info,blck) subroutine psb_siluk_fact(fill_in,ialg,a,l,u,d,info,blck,shft)
use psb_base_mod use psb_base_mod
use psb_s_ilu_fact_mod, psb_protect_name => psb_siluk_fact use psb_s_ilu_fact_mod, psb_protect_name => psb_siluk_fact
@ -141,6 +141,7 @@ subroutine psb_siluk_fact(fill_in,ialg,a,l,u,d,info,blck)
type(psb_sspmat_type),intent(inout) :: l,u type(psb_sspmat_type),intent(inout) :: l,u
type(psb_sspmat_type),intent(in), optional, target :: blck type(psb_sspmat_type),intent(in), optional, target :: blck
real(psb_spk_), intent(inout) :: d(:) real(psb_spk_), intent(inout) :: d(:)
real(psb_spk_), intent(in), optional :: shft
! Local Variables ! Local Variables
integer(psb_ipk_) :: l1, l2, m, err_act integer(psb_ipk_) :: l1, l2, m, err_act
@ -184,7 +185,7 @@ subroutine psb_siluk_fact(fill_in,ialg,a,l,u,d,info,blck)
! Compute the ILU(k) or the MILU(k) factorization, depending on ialg ! Compute the ILU(k) or the MILU(k) factorization, depending on ialg
! !
call psb_siluk_factint(fill_in,ialg,a,blck_,& call psb_siluk_factint(fill_in,ialg,a,blck_,&
& d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,info) & d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,info,shft)
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_siluk_factint' ch_err='psb_siluk_factint'
@ -298,7 +299,7 @@ contains
! Error code. ! Error code.
! !
subroutine psb_siluk_factint(fill_in,ialg,a,b,& subroutine psb_siluk_factint(fill_in,ialg,a,b,&
& d,lval,lja,lirp,uval,uja,uirp,l1,l2,info) & d,lval,lja,lirp,uval,uja,uirp,l1,l2,info,shft)
use psb_base_mod use psb_base_mod
@ -311,6 +312,7 @@ contains
integer(psb_ipk_), allocatable, intent(inout) :: lja(:),lirp(:),uja(:),uirp(:) integer(psb_ipk_), allocatable, intent(inout) :: lja(:),lirp(:),uja(:),uirp(:)
real(psb_spk_), allocatable, intent(inout) :: lval(:),uval(:) real(psb_spk_), allocatable, intent(inout) :: lval(:),uval(:)
real(psb_spk_), intent(inout) :: d(:) real(psb_spk_), intent(inout) :: d(:)
real(psb_spk_), intent(in) :: shft
! Local variables ! Local variables
integer(psb_ipk_) :: ma,mb,i, ktrw,err_act,nidx, m integer(psb_ipk_) :: ma,mb,i, ktrw,err_act,nidx, m
@ -400,13 +402,13 @@ contains
! !
! Copy into trw the i-th local row of the matrix, stored in a ! Copy into trw the i-th local row of the matrix, stored in a
! !
call iluk_copyin(i,ma,a,ione,m,row,rowlevs,heap,ktrw,trw,info) call iluk_copyin(i,ma,a,ione,m,row,rowlevs,heap,ktrw,trw,info,shft)
else else
! !
! Copy into trw the i-th local row of the matrix, stored in b ! Copy into trw the i-th local row of the matrix, stored in b
! (as (i-ma)-th row) ! (as (i-ma)-th row)
! !
call iluk_copyin(i-ma,mb,b,ione,m,row,rowlevs,heap,ktrw,trw,info) call iluk_copyin(i-ma,mb,b,ione,m,row,rowlevs,heap,ktrw,trw,info,shft)
endif endif
! Do an elimination step on the current row. It turns out we only ! Do an elimination step on the current row. It turns out we only
@ -516,7 +518,7 @@ contains
! until we empty the buffer. Thus we will make a call to psb_sp_getblk ! until we empty the buffer. Thus we will make a call to psb_sp_getblk
! every nrb calls to copyin. If A is in CSR format it is unused. ! every nrb calls to copyin. If A is in CSR format it is unused.
! !
subroutine iluk_copyin(i,m,a,jmin,jmax,row,rowlevs,heap,ktrw,trw,info) subroutine iluk_copyin(i,m,a,jmin,jmax,row,rowlevs,heap,ktrw,trw,info,shft)
use psb_base_mod use psb_base_mod
@ -530,6 +532,8 @@ contains
integer(psb_ipk_), intent(inout) :: rowlevs(:) integer(psb_ipk_), intent(inout) :: rowlevs(:)
real(psb_spk_), intent(inout) :: row(:) real(psb_spk_), intent(inout) :: row(:)
type(psb_i_heap), intent(inout) :: heap type(psb_i_heap), intent(inout) :: heap
real(psb_spk_), intent(in) :: shft
! Local variables ! Local variables
integer(psb_ipk_) :: k,j,irb,err_act,nz integer(psb_ipk_) :: k,j,irb,err_act,nz
@ -554,6 +558,7 @@ contains
k = aa%ja(j) k = aa%ja(j)
if ((jmin<=k).and.(k<=jmax)) then if ((jmin<=k).and.(k<=jmax)) then
row(k) = aa%val(j) row(k) = aa%val(j)
if (k==i) row(k) = row(k) + shft
rowlevs(k) = 0 rowlevs(k) = 0
call heap%insert(k,info) call heap%insert(k,info)
end if end if
@ -587,6 +592,7 @@ contains
k = trw%ja(ktrw) k = trw%ja(ktrw)
if ((jmin<=k).and.(k<=jmax)) then if ((jmin<=k).and.(k<=jmax)) then
row(k) = trw%val(ktrw) row(k) = trw%val(ktrw)
if (k==i) row(k) = row(k) + shft
rowlevs(k) = 0 rowlevs(k) = 0
call heap%insert(k,info) call heap%insert(k,info)
end if end if
@ -670,7 +676,8 @@ contains
! Note: this argument is intent(inout) and not only intent(out) ! Note: this argument is intent(inout) and not only intent(out)
! to retain its allocation, done by this routine. ! to retain its allocation, done by this routine.
! !
subroutine iluk_fact(fill_in,i,row,rowlevs,heap,d,uja,uirp,uval,uplevs,nidx,idxs,info) subroutine iluk_fact(fill_in,i,row,rowlevs,heap,d,&
& uja,uirp,uval,uplevs,nidx,idxs,info)
use psb_base_mod use psb_base_mod

@ -123,7 +123,7 @@
! greater than 0. If the overlap is 0 or the matrix has been reordered ! greater than 0. If the overlap is 0 or the matrix has been reordered
! (see psb_fact_bld), then blck does not contain any row. ! (see psb_fact_bld), then blck does not contain any row.
! !
subroutine psb_silut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) subroutine psb_silut_fact(fill_in,thres,a,l,u,d,info,blck,iscale,shft)
use psb_base_mod use psb_base_mod
use psb_s_ilu_fact_mod, psb_protect_name => psb_silut_fact use psb_s_ilu_fact_mod, psb_protect_name => psb_silut_fact
@ -139,6 +139,7 @@ subroutine psb_silut_fact(fill_in,thres,a,l,u,d,info,blck,iscale)
real(psb_spk_), intent(inout) :: d(:) real(psb_spk_), intent(inout) :: d(:)
type(psb_sspmat_type),intent(in), optional, target :: blck type(psb_sspmat_type),intent(in), optional, target :: blck
integer(psb_ipk_), intent(in), optional :: iscale integer(psb_ipk_), intent(in), optional :: iscale
real(psb_spk_), intent(in), optional :: shft
! Local Variables ! Local Variables
integer(psb_ipk_) :: l1, l2, m, err_act, iscale_ integer(psb_ipk_) :: l1, l2, m, err_act, iscale_
@ -206,7 +207,7 @@ subroutine psb_silut_fact(fill_in,thres,a,l,u,d,info,blck,iscale)
! Compute the ILU(k,t) factorization ! Compute the ILU(k,t) factorization
! !
call psb_silut_factint(fill_in,thres,a,blck_,& call psb_silut_factint(fill_in,thres,a,blck_,&
& d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,info,scale) & d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,info,scale,shft)
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_silut_factint' ch_err='psb_silut_factint'
@ -316,7 +317,7 @@ contains
! Error code. ! Error code.
! !
subroutine psb_silut_factint(fill_in,thres,a,b,& subroutine psb_silut_factint(fill_in,thres,a,b,&
& d,lval,lja,lirp,uval,uja,uirp,l1,l2,info,scale) & d,lval,lja,lirp,uval,uja,uirp,l1,l2,info,scale,shft)
use psb_base_mod use psb_base_mod
@ -331,6 +332,7 @@ contains
real(psb_spk_), allocatable, intent(inout) :: lval(:),uval(:) real(psb_spk_), allocatable, intent(inout) :: lval(:),uval(:)
real(psb_spk_), intent(inout) :: d(:) real(psb_spk_), intent(inout) :: d(:)
real(psb_spk_), intent(in), optional :: scale real(psb_spk_), intent(in), optional :: scale
real(psb_spk_), intent(in) :: shft
! Local Variables ! Local Variables
integer(psb_ipk_) :: i, ktrw,err_act,nidx,nlw,nup,jmaxup, ma, mb, m integer(psb_ipk_) :: i, ktrw,err_act,nidx,nlw,nup,jmaxup, ma, mb, m
@ -401,10 +403,10 @@ contains
d(i) = czero d(i) = czero
if (i<=ma) then if (i<=ma) then
call ilut_copyin(i,ma,a,i,ione,m,nlw,nup,jmaxup,nrmi,weight,& call ilut_copyin(i,ma,a,i,ione,m,nlw,nup,jmaxup,nrmi,weight,&
& row,heap,ktrw,trw,info) & row,heap,ktrw,trw,info,shft)
else else
call ilut_copyin(i-ma,mb,b,i,ione,m,nlw,nup,jmaxup,nrmi,weight,& call ilut_copyin(i-ma,mb,b,i,ione,m,nlw,nup,jmaxup,nrmi,weight,&
& row,heap,ktrw,trw,info) & row,heap,ktrw,trw,info,shft)
endif endif
! !
@ -540,7 +542,7 @@ contains
! every nrb calls to copyin. If A is in CSR format it is unused. ! every nrb calls to copyin. If A is in CSR format it is unused.
! !
subroutine ilut_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,& subroutine ilut_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,&
& nrmi,weight,row,heap,ktrw,trw,info) & nrmi,weight,row,heap,ktrw,trw,info,shft)
use psb_base_mod use psb_base_mod
implicit none implicit none
type(psb_sspmat_type), intent(in) :: a type(psb_sspmat_type), intent(in) :: a
@ -551,6 +553,7 @@ contains
real(psb_spk_), intent(inout) :: row(:) real(psb_spk_), intent(inout) :: row(:)
real(psb_spk_), intent(in) :: weight real(psb_spk_), intent(in) :: weight
type(psb_i_heap), intent(inout) :: heap type(psb_i_heap), intent(inout) :: heap
real(psb_spk_), intent(in) :: shft
integer(psb_ipk_) :: k,j,irb,kin,nz integer(psb_ipk_) :: k,j,irb,kin,nz
integer(psb_ipk_), parameter :: nrb=40 integer(psb_ipk_), parameter :: nrb=40
@ -597,6 +600,7 @@ contains
call heap%insert(k,info) call heap%insert(k,info)
if (info /= psb_success_) exit if (info /= psb_success_) exit
if (k<jd) nlw = nlw + 1 if (k<jd) nlw = nlw + 1
if (k == jd) row(k) = row(k) + (shft*weight)
if (k>jd) then if (k>jd) then
nup = nup + 1 nup = nup + 1
if (abs(row(k))>dmaxup) then if (abs(row(k))>dmaxup) then
@ -648,6 +652,7 @@ contains
call heap%insert(k,info) call heap%insert(k,info)
if (info /= psb_success_) exit if (info /= psb_success_) exit
if (k<jd) nlw = nlw + 1 if (k<jd) nlw = nlw + 1
if (k == jd) row(k) = row(k) + (shft*weight)
if (k>jd) then if (k>jd) then
nup = nup + 1 nup = nup + 1
if (abs(row(k))>dmaxup) then if (abs(row(k))>dmaxup) then

@ -130,7 +130,7 @@
! greater than 0. If the overlap is 0 or the matrix has been reordered ! greater than 0. If the overlap is 0 or the matrix has been reordered
! (see psb_fact_bld), then blck is empty. ! (see psb_fact_bld), then blck is empty.
! !
subroutine psb_zilu0_fact(ialg,a,l,u,d,info,blck, upd) subroutine psb_zilu0_fact(ialg,a,l,u,d,info,blck, upd,shft)
use psb_base_mod use psb_base_mod
use psb_z_ilu_fact_mod, psb_protect_name => psb_zilu0_fact use psb_z_ilu_fact_mod, psb_protect_name => psb_zilu0_fact
@ -145,11 +145,13 @@ subroutine psb_zilu0_fact(ialg,a,l,u,d,info,blck, upd)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
type(psb_zspmat_type),intent(in), optional, target :: blck type(psb_zspmat_type),intent(in), optional, target :: blck
character, intent(in), optional :: upd character, intent(in), optional :: upd
complex(psb_dpk_), intent(in), optional :: shft
! Local variables ! Local variables
integer(psb_ipk_) :: l1, l2, m, err_act integer(psb_ipk_) :: l1, l2, m, err_act
type(psb_zspmat_type), pointer :: blck_ type(psb_zspmat_type), pointer :: blck_
type(psb_z_csr_sparse_mat) :: ll, uu type(psb_z_csr_sparse_mat) :: ll, uu
complex(psb_dpk_) :: shft_
character :: upd_ character :: upd_
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
@ -177,6 +179,11 @@ subroutine psb_zilu0_fact(ialg,a,l,u,d,info,blck, upd)
else else
upd_ = 'F' upd_ = 'F'
end if end if
if (present(shft)) then
shft_ = shft
else
shft_ = zzero
end if
m = a%get_nrows() + blck_%get_nrows() m = a%get_nrows() + blck_%get_nrows()
if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.& if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.&
@ -193,7 +200,7 @@ subroutine psb_zilu0_fact(ialg,a,l,u,d,info,blck, upd)
! Compute the ILU(0) or the MILU(0) factorization, depending on ialg ! Compute the ILU(0) or the MILU(0) factorization, depending on ialg
! !
call psb_zilu0_factint(ialg,a,blck_,& call psb_zilu0_factint(ialg,a,blck_,&
& d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,upd_,info) & d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,upd_,shft_,info)
if(info.ne.0) then if(info.ne.0) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_zilu0_factint' ch_err='psb_zilu0_factint'
@ -314,7 +321,7 @@ contains
! Error code. ! Error code.
! !
subroutine psb_zilu0_factint(ialg,a,b,& subroutine psb_zilu0_factint(ialg,a,b,&
& d,lval,lja,lirp,uval,uja,uirp,l1,l2,upd,info) & d,lval,lja,lirp,uval,uja,uirp,l1,l2,upd,shft,info)
implicit none implicit none
@ -325,6 +332,7 @@ contains
integer(psb_ipk_), intent(inout) :: lja(:),lirp(:),uja(:),uirp(:) integer(psb_ipk_), intent(inout) :: lja(:),lirp(:),uja(:),uirp(:)
complex(psb_dpk_), intent(inout) :: lval(:),uval(:),d(:) complex(psb_dpk_), intent(inout) :: lval(:),uval(:),d(:)
character, intent(in) :: upd character, intent(in) :: upd
complex(psb_dpk_), intent(in) :: shft
! Local variables ! Local variables
integer(psb_ipk_) :: i,j,k,l,low1,low2,kk,jj,ll, ktrw,err_act, m integer(psb_ipk_) :: i,j,k,l,low1,low2,kk,jj,ll, ktrw,err_act, m
@ -382,14 +390,14 @@ contains
! into lval/d(i)/uval ! into lval/d(i)/uval
! !
call ilu_copyin(i,ma,a,i,ione,m,l1,lja,lval,& call ilu_copyin(i,ma,a,i,ione,m,l1,lja,lval,&
& d(i),l2,uja,uval,ktrw,trw,upd) & d(i),l2,uja,uval,ktrw,trw,upd,shft)
else else
! !
! Copy the i-th local row of the matrix, stored in b ! Copy the i-th local row of the matrix, stored in b
! (as (i-ma)-th row), into lval/d(i)/uval ! (as (i-ma)-th row), into lval/d(i)/uval
! !
call ilu_copyin(i-ma,mb,b,i,ione,m,l1,lja,lval,& call ilu_copyin(i-ma,mb,b,i,ione,m,l1,lja,lval,&
& d(i),l2,uja,uval,ktrw,trw,upd) & d(i),l2,uja,uval,ktrw,trw,upd,shft)
endif endif
lirp(i+1) = l1 + 1 lirp(i+1) = l1 + 1
@ -583,7 +591,7 @@ contains
! every nrb calls to copyin. If A is in CSR format it is unused. ! every nrb calls to copyin. If A is in CSR format it is unused.
! !
subroutine ilu_copyin(i,m,a,jd,jmin,jmax,l1,lja,lval,& subroutine ilu_copyin(i,m,a,jd,jmin,jmax,l1,lja,lval,&
& dia,l2,uja,uval,ktrw,trw,upd) & dia,l2,uja,uval,ktrw,trw,upd,shft)
use psb_base_mod use psb_base_mod
@ -597,6 +605,7 @@ contains
integer(psb_ipk_), intent(inout) :: lja(:), uja(:) integer(psb_ipk_), intent(inout) :: lja(:), uja(:)
complex(psb_dpk_), intent(inout) :: lval(:), uval(:), dia complex(psb_dpk_), intent(inout) :: lval(:), uval(:), dia
character, intent(in) :: upd character, intent(in) :: upd
complex(psb_dpk_), intent(in) :: shft
! Local variables ! Local variables
integer(psb_ipk_) :: k,j,info,irb, nz integer(psb_ipk_) :: k,j,info,irb, nz
integer(psb_ipk_), parameter :: nrb=40 integer(psb_ipk_), parameter :: nrb=40
@ -625,7 +634,7 @@ contains
lval(l1) = aa%val(j) lval(l1) = aa%val(j)
lja(l1) = k lja(l1) = k
else if (k == jd) then else if (k == jd) then
dia = aa%val(j) dia = aa%val(j) + shft
else if ((k > jd).and.(k <= jmax)) then else if ((k > jd).and.(k <= jmax)) then
l2 = l2 + 1 l2 = l2 + 1
uval(l2) = aa%val(j) uval(l2) = aa%val(j)
@ -665,7 +674,7 @@ contains
lval(l1) = trw%val(ktrw) lval(l1) = trw%val(ktrw)
lja(l1) = k lja(l1) = k
else if (k == jd) then else if (k == jd) then
dia = trw%val(ktrw) dia = trw%val(ktrw) + shft
else if ((k > jd).and.(k <= jmax)) then else if ((k > jd).and.(k <= jmax)) then
l2 = l2 + 1 l2 = l2 + 1
uval(l2) = trw%val(ktrw) uval(l2) = trw%val(ktrw)

@ -127,7 +127,7 @@
! greater than 0. If the overlap is 0 or the matrix has been reordered ! greater than 0. If the overlap is 0 or the matrix has been reordered
! (see psb_fact_bld), then blck does not contain any row. ! (see psb_fact_bld), then blck does not contain any row.
! !
subroutine psb_ziluk_fact(fill_in,ialg,a,l,u,d,info,blck) subroutine psb_ziluk_fact(fill_in,ialg,a,l,u,d,info,blck,shft)
use psb_base_mod use psb_base_mod
use psb_z_ilu_fact_mod, psb_protect_name => psb_ziluk_fact use psb_z_ilu_fact_mod, psb_protect_name => psb_ziluk_fact
@ -141,6 +141,7 @@ subroutine psb_ziluk_fact(fill_in,ialg,a,l,u,d,info,blck)
type(psb_zspmat_type),intent(inout) :: l,u type(psb_zspmat_type),intent(inout) :: l,u
type(psb_zspmat_type),intent(in), optional, target :: blck type(psb_zspmat_type),intent(in), optional, target :: blck
complex(psb_dpk_), intent(inout) :: d(:) complex(psb_dpk_), intent(inout) :: d(:)
complex(psb_dpk_), intent(in), optional :: shft
! Local Variables ! Local Variables
integer(psb_ipk_) :: l1, l2, m, err_act integer(psb_ipk_) :: l1, l2, m, err_act
@ -184,7 +185,7 @@ subroutine psb_ziluk_fact(fill_in,ialg,a,l,u,d,info,blck)
! Compute the ILU(k) or the MILU(k) factorization, depending on ialg ! Compute the ILU(k) or the MILU(k) factorization, depending on ialg
! !
call psb_ziluk_factint(fill_in,ialg,a,blck_,& call psb_ziluk_factint(fill_in,ialg,a,blck_,&
& d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,info) & d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,info,shft)
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_ziluk_factint' ch_err='psb_ziluk_factint'
@ -298,7 +299,7 @@ contains
! Error code. ! Error code.
! !
subroutine psb_ziluk_factint(fill_in,ialg,a,b,& subroutine psb_ziluk_factint(fill_in,ialg,a,b,&
& d,lval,lja,lirp,uval,uja,uirp,l1,l2,info) & d,lval,lja,lirp,uval,uja,uirp,l1,l2,info,shft)
use psb_base_mod use psb_base_mod
@ -311,6 +312,7 @@ contains
integer(psb_ipk_), allocatable, intent(inout) :: lja(:),lirp(:),uja(:),uirp(:) integer(psb_ipk_), allocatable, intent(inout) :: lja(:),lirp(:),uja(:),uirp(:)
complex(psb_dpk_), allocatable, intent(inout) :: lval(:),uval(:) complex(psb_dpk_), allocatable, intent(inout) :: lval(:),uval(:)
complex(psb_dpk_), intent(inout) :: d(:) complex(psb_dpk_), intent(inout) :: d(:)
complex(psb_dpk_), intent(in) :: shft
! Local variables ! Local variables
integer(psb_ipk_) :: ma,mb,i, ktrw,err_act,nidx, m integer(psb_ipk_) :: ma,mb,i, ktrw,err_act,nidx, m
@ -400,13 +402,13 @@ contains
! !
! Copy into trw the i-th local row of the matrix, stored in a ! Copy into trw the i-th local row of the matrix, stored in a
! !
call iluk_copyin(i,ma,a,ione,m,row,rowlevs,heap,ktrw,trw,info) call iluk_copyin(i,ma,a,ione,m,row,rowlevs,heap,ktrw,trw,info,shft)
else else
! !
! Copy into trw the i-th local row of the matrix, stored in b ! Copy into trw the i-th local row of the matrix, stored in b
! (as (i-ma)-th row) ! (as (i-ma)-th row)
! !
call iluk_copyin(i-ma,mb,b,ione,m,row,rowlevs,heap,ktrw,trw,info) call iluk_copyin(i-ma,mb,b,ione,m,row,rowlevs,heap,ktrw,trw,info,shft)
endif endif
! Do an elimination step on the current row. It turns out we only ! Do an elimination step on the current row. It turns out we only
@ -516,7 +518,7 @@ contains
! until we empty the buffer. Thus we will make a call to psb_sp_getblk ! until we empty the buffer. Thus we will make a call to psb_sp_getblk
! every nrb calls to copyin. If A is in CSR format it is unused. ! every nrb calls to copyin. If A is in CSR format it is unused.
! !
subroutine iluk_copyin(i,m,a,jmin,jmax,row,rowlevs,heap,ktrw,trw,info) subroutine iluk_copyin(i,m,a,jmin,jmax,row,rowlevs,heap,ktrw,trw,info,shft)
use psb_base_mod use psb_base_mod
@ -530,6 +532,8 @@ contains
integer(psb_ipk_), intent(inout) :: rowlevs(:) integer(psb_ipk_), intent(inout) :: rowlevs(:)
complex(psb_dpk_), intent(inout) :: row(:) complex(psb_dpk_), intent(inout) :: row(:)
type(psb_i_heap), intent(inout) :: heap type(psb_i_heap), intent(inout) :: heap
complex(psb_dpk_), intent(in) :: shft
! Local variables ! Local variables
integer(psb_ipk_) :: k,j,irb,err_act,nz integer(psb_ipk_) :: k,j,irb,err_act,nz
@ -554,6 +558,7 @@ contains
k = aa%ja(j) k = aa%ja(j)
if ((jmin<=k).and.(k<=jmax)) then if ((jmin<=k).and.(k<=jmax)) then
row(k) = aa%val(j) row(k) = aa%val(j)
if (k==i) row(k) = row(k) + shft
rowlevs(k) = 0 rowlevs(k) = 0
call heap%insert(k,info) call heap%insert(k,info)
end if end if
@ -587,6 +592,7 @@ contains
k = trw%ja(ktrw) k = trw%ja(ktrw)
if ((jmin<=k).and.(k<=jmax)) then if ((jmin<=k).and.(k<=jmax)) then
row(k) = trw%val(ktrw) row(k) = trw%val(ktrw)
if (k==i) row(k) = row(k) + shft
rowlevs(k) = 0 rowlevs(k) = 0
call heap%insert(k,info) call heap%insert(k,info)
end if end if
@ -670,7 +676,8 @@ contains
! Note: this argument is intent(inout) and not only intent(out) ! Note: this argument is intent(inout) and not only intent(out)
! to retain its allocation, done by this routine. ! to retain its allocation, done by this routine.
! !
subroutine iluk_fact(fill_in,i,row,rowlevs,heap,d,uja,uirp,uval,uplevs,nidx,idxs,info) subroutine iluk_fact(fill_in,i,row,rowlevs,heap,d,&
& uja,uirp,uval,uplevs,nidx,idxs,info)
use psb_base_mod use psb_base_mod

@ -123,7 +123,7 @@
! greater than 0. If the overlap is 0 or the matrix has been reordered ! greater than 0. If the overlap is 0 or the matrix has been reordered
! (see psb_fact_bld), then blck does not contain any row. ! (see psb_fact_bld), then blck does not contain any row.
! !
subroutine psb_zilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) subroutine psb_zilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale,shft)
use psb_base_mod use psb_base_mod
use psb_z_ilu_fact_mod, psb_protect_name => psb_zilut_fact use psb_z_ilu_fact_mod, psb_protect_name => psb_zilut_fact
@ -139,6 +139,7 @@ subroutine psb_zilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale)
complex(psb_dpk_), intent(inout) :: d(:) complex(psb_dpk_), intent(inout) :: d(:)
type(psb_zspmat_type),intent(in), optional, target :: blck type(psb_zspmat_type),intent(in), optional, target :: blck
integer(psb_ipk_), intent(in), optional :: iscale integer(psb_ipk_), intent(in), optional :: iscale
complex(psb_dpk_), intent(in), optional :: shft
! Local Variables ! Local Variables
integer(psb_ipk_) :: l1, l2, m, err_act, iscale_ integer(psb_ipk_) :: l1, l2, m, err_act, iscale_
@ -206,7 +207,7 @@ subroutine psb_zilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale)
! Compute the ILU(k,t) factorization ! Compute the ILU(k,t) factorization
! !
call psb_zilut_factint(fill_in,thres,a,blck_,& call psb_zilut_factint(fill_in,thres,a,blck_,&
& d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,info,scale) & d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,info,scale,shft)
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_zilut_factint' ch_err='psb_zilut_factint'
@ -316,7 +317,7 @@ contains
! Error code. ! Error code.
! !
subroutine psb_zilut_factint(fill_in,thres,a,b,& subroutine psb_zilut_factint(fill_in,thres,a,b,&
& d,lval,lja,lirp,uval,uja,uirp,l1,l2,info,scale) & d,lval,lja,lirp,uval,uja,uirp,l1,l2,info,scale,shft)
use psb_base_mod use psb_base_mod
@ -331,6 +332,7 @@ contains
complex(psb_dpk_), allocatable, intent(inout) :: lval(:),uval(:) complex(psb_dpk_), allocatable, intent(inout) :: lval(:),uval(:)
complex(psb_dpk_), intent(inout) :: d(:) complex(psb_dpk_), intent(inout) :: d(:)
real(psb_dpk_), intent(in), optional :: scale real(psb_dpk_), intent(in), optional :: scale
complex(psb_dpk_), intent(in) :: shft
! Local Variables ! Local Variables
integer(psb_ipk_) :: i, ktrw,err_act,nidx,nlw,nup,jmaxup, ma, mb, m integer(psb_ipk_) :: i, ktrw,err_act,nidx,nlw,nup,jmaxup, ma, mb, m
@ -401,10 +403,10 @@ contains
d(i) = czero d(i) = czero
if (i<=ma) then if (i<=ma) then
call ilut_copyin(i,ma,a,i,ione,m,nlw,nup,jmaxup,nrmi,weight,& call ilut_copyin(i,ma,a,i,ione,m,nlw,nup,jmaxup,nrmi,weight,&
& row,heap,ktrw,trw,info) & row,heap,ktrw,trw,info,shft)
else else
call ilut_copyin(i-ma,mb,b,i,ione,m,nlw,nup,jmaxup,nrmi,weight,& call ilut_copyin(i-ma,mb,b,i,ione,m,nlw,nup,jmaxup,nrmi,weight,&
& row,heap,ktrw,trw,info) & row,heap,ktrw,trw,info,shft)
endif endif
! !
@ -540,7 +542,7 @@ contains
! every nrb calls to copyin. If A is in CSR format it is unused. ! every nrb calls to copyin. If A is in CSR format it is unused.
! !
subroutine ilut_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,& subroutine ilut_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,&
& nrmi,weight,row,heap,ktrw,trw,info) & nrmi,weight,row,heap,ktrw,trw,info,shft)
use psb_base_mod use psb_base_mod
implicit none implicit none
type(psb_zspmat_type), intent(in) :: a type(psb_zspmat_type), intent(in) :: a
@ -551,6 +553,7 @@ contains
complex(psb_dpk_), intent(inout) :: row(:) complex(psb_dpk_), intent(inout) :: row(:)
real(psb_dpk_), intent(in) :: weight real(psb_dpk_), intent(in) :: weight
type(psb_i_heap), intent(inout) :: heap type(psb_i_heap), intent(inout) :: heap
complex(psb_dpk_), intent(in) :: shft
integer(psb_ipk_) :: k,j,irb,kin,nz integer(psb_ipk_) :: k,j,irb,kin,nz
integer(psb_ipk_), parameter :: nrb=40 integer(psb_ipk_), parameter :: nrb=40
@ -597,6 +600,7 @@ contains
call heap%insert(k,info) call heap%insert(k,info)
if (info /= psb_success_) exit if (info /= psb_success_) exit
if (k<jd) nlw = nlw + 1 if (k<jd) nlw = nlw + 1
if (k == jd) row(k) = row(k) + (shft*weight)
if (k>jd) then if (k>jd) then
nup = nup + 1 nup = nup + 1
if (abs(row(k))>dmaxup) then if (abs(row(k))>dmaxup) then
@ -648,6 +652,7 @@ contains
call heap%insert(k,info) call heap%insert(k,info)
if (info /= psb_success_) exit if (info /= psb_success_) exit
if (k<jd) nlw = nlw + 1 if (k<jd) nlw = nlw + 1
if (k == jd) row(k) = row(k) + (shft*weight)
if (k>jd) then if (k>jd) then
nup = nup + 1 nup = nup + 1
if (abs(row(k))>dmaxup) then if (abs(row(k))>dmaxup) then

@ -80,7 +80,7 @@ module psb_c_ilu_fact_mod
use psb_base_mod use psb_base_mod
use psb_prec_const_mod use psb_prec_const_mod
interface psb_ilu0_fact interface psb_ilu0_fact
subroutine psb_cilu0_fact(ialg,a,l,u,d,info,blck,upd) subroutine psb_cilu0_fact(ialg,a,l,u,d,info,blck,upd,shft)
import psb_cspmat_type, psb_spk_, psb_ipk_ import psb_cspmat_type, psb_spk_, psb_ipk_
integer(psb_ipk_), intent(in) :: ialg integer(psb_ipk_), intent(in) :: ialg
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -89,11 +89,12 @@ module psb_c_ilu_fact_mod
type(psb_cspmat_type),intent(in), optional, target :: blck type(psb_cspmat_type),intent(in), optional, target :: blck
character, intent(in), optional :: upd character, intent(in), optional :: upd
complex(psb_spk_), intent(inout) :: d(:) complex(psb_spk_), intent(inout) :: d(:)
complex(psb_spk_), intent(in), optional :: shft
end subroutine psb_cilu0_fact end subroutine psb_cilu0_fact
end interface end interface
interface psb_iluk_fact interface psb_iluk_fact
subroutine psb_ciluk_fact(fill_in,ialg,a,l,u,d,info,blck) subroutine psb_ciluk_fact(fill_in,ialg,a,l,u,d,info,blck,shft)
import psb_cspmat_type, psb_spk_, psb_ipk_ import psb_cspmat_type, psb_spk_, psb_ipk_
integer(psb_ipk_), intent(in) :: fill_in,ialg integer(psb_ipk_), intent(in) :: fill_in,ialg
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -101,11 +102,12 @@ module psb_c_ilu_fact_mod
type(psb_cspmat_type),intent(inout) :: l,u type(psb_cspmat_type),intent(inout) :: l,u
type(psb_cspmat_type),intent(in), optional, target :: blck type(psb_cspmat_type),intent(in), optional, target :: blck
complex(psb_spk_), intent(inout) :: d(:) complex(psb_spk_), intent(inout) :: d(:)
complex(psb_spk_), intent(in), optional :: shft
end subroutine psb_ciluk_fact end subroutine psb_ciluk_fact
end interface end interface
interface psb_ilut_fact interface psb_ilut_fact
subroutine psb_cilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) subroutine psb_cilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale,shft)
import psb_cspmat_type, psb_spk_, psb_ipk_ import psb_cspmat_type, psb_spk_, psb_ipk_
integer(psb_ipk_), intent(in) :: fill_in integer(psb_ipk_), intent(in) :: fill_in
real(psb_spk_), intent(in) :: thres real(psb_spk_), intent(in) :: thres
@ -115,6 +117,7 @@ module psb_c_ilu_fact_mod
complex(psb_spk_), intent(inout) :: d(:) complex(psb_spk_), intent(inout) :: d(:)
type(psb_cspmat_type),intent(in), optional, target :: blck type(psb_cspmat_type),intent(in), optional, target :: blck
integer(psb_ipk_), intent(in), optional :: iscale integer(psb_ipk_), intent(in), optional :: iscale
complex(psb_spk_), intent(in), optional :: shft
end subroutine psb_cilut_fact end subroutine psb_cilut_fact
end interface end interface

@ -80,7 +80,7 @@ module psb_d_ilu_fact_mod
use psb_base_mod use psb_base_mod
use psb_prec_const_mod use psb_prec_const_mod
interface psb_ilu0_fact interface psb_ilu0_fact
subroutine psb_dilu0_fact(ialg,a,l,u,d,info,blck,upd) subroutine psb_dilu0_fact(ialg,a,l,u,d,info,blck,upd,shft)
import psb_dspmat_type, psb_dpk_, psb_ipk_ import psb_dspmat_type, psb_dpk_, psb_ipk_
integer(psb_ipk_), intent(in) :: ialg integer(psb_ipk_), intent(in) :: ialg
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -89,11 +89,12 @@ module psb_d_ilu_fact_mod
type(psb_dspmat_type),intent(in), optional, target :: blck type(psb_dspmat_type),intent(in), optional, target :: blck
character, intent(in), optional :: upd character, intent(in), optional :: upd
real(psb_dpk_), intent(inout) :: d(:) real(psb_dpk_), intent(inout) :: d(:)
real(psb_dpk_), intent(in), optional :: shft
end subroutine psb_dilu0_fact end subroutine psb_dilu0_fact
end interface end interface
interface psb_iluk_fact interface psb_iluk_fact
subroutine psb_diluk_fact(fill_in,ialg,a,l,u,d,info,blck) subroutine psb_diluk_fact(fill_in,ialg,a,l,u,d,info,blck,shft)
import psb_dspmat_type, psb_dpk_, psb_ipk_ import psb_dspmat_type, psb_dpk_, psb_ipk_
integer(psb_ipk_), intent(in) :: fill_in,ialg integer(psb_ipk_), intent(in) :: fill_in,ialg
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -101,11 +102,12 @@ module psb_d_ilu_fact_mod
type(psb_dspmat_type),intent(inout) :: l,u type(psb_dspmat_type),intent(inout) :: l,u
type(psb_dspmat_type),intent(in), optional, target :: blck type(psb_dspmat_type),intent(in), optional, target :: blck
real(psb_dpk_), intent(inout) :: d(:) real(psb_dpk_), intent(inout) :: d(:)
real(psb_dpk_), intent(in), optional :: shft
end subroutine psb_diluk_fact end subroutine psb_diluk_fact
end interface end interface
interface psb_ilut_fact interface psb_ilut_fact
subroutine psb_dilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) subroutine psb_dilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale,shft)
import psb_dspmat_type, psb_dpk_, psb_ipk_ import psb_dspmat_type, psb_dpk_, psb_ipk_
integer(psb_ipk_), intent(in) :: fill_in integer(psb_ipk_), intent(in) :: fill_in
real(psb_dpk_), intent(in) :: thres real(psb_dpk_), intent(in) :: thres
@ -115,6 +117,7 @@ module psb_d_ilu_fact_mod
real(psb_dpk_), intent(inout) :: d(:) real(psb_dpk_), intent(inout) :: d(:)
type(psb_dspmat_type),intent(in), optional, target :: blck type(psb_dspmat_type),intent(in), optional, target :: blck
integer(psb_ipk_), intent(in), optional :: iscale integer(psb_ipk_), intent(in), optional :: iscale
real(psb_dpk_), intent(in), optional :: shft
end subroutine psb_dilut_fact end subroutine psb_dilut_fact
end interface end interface

@ -80,7 +80,7 @@ module psb_s_ilu_fact_mod
use psb_base_mod use psb_base_mod
use psb_prec_const_mod use psb_prec_const_mod
interface psb_ilu0_fact interface psb_ilu0_fact
subroutine psb_silu0_fact(ialg,a,l,u,d,info,blck,upd) subroutine psb_silu0_fact(ialg,a,l,u,d,info,blck,upd,shft)
import psb_sspmat_type, psb_spk_, psb_ipk_ import psb_sspmat_type, psb_spk_, psb_ipk_
integer(psb_ipk_), intent(in) :: ialg integer(psb_ipk_), intent(in) :: ialg
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -89,11 +89,12 @@ module psb_s_ilu_fact_mod
type(psb_sspmat_type),intent(in), optional, target :: blck type(psb_sspmat_type),intent(in), optional, target :: blck
character, intent(in), optional :: upd character, intent(in), optional :: upd
real(psb_spk_), intent(inout) :: d(:) real(psb_spk_), intent(inout) :: d(:)
real(psb_spk_), intent(in), optional :: shft
end subroutine psb_silu0_fact end subroutine psb_silu0_fact
end interface end interface
interface psb_iluk_fact interface psb_iluk_fact
subroutine psb_siluk_fact(fill_in,ialg,a,l,u,d,info,blck) subroutine psb_siluk_fact(fill_in,ialg,a,l,u,d,info,blck,shft)
import psb_sspmat_type, psb_spk_, psb_ipk_ import psb_sspmat_type, psb_spk_, psb_ipk_
integer(psb_ipk_), intent(in) :: fill_in,ialg integer(psb_ipk_), intent(in) :: fill_in,ialg
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -101,11 +102,12 @@ module psb_s_ilu_fact_mod
type(psb_sspmat_type),intent(inout) :: l,u type(psb_sspmat_type),intent(inout) :: l,u
type(psb_sspmat_type),intent(in), optional, target :: blck type(psb_sspmat_type),intent(in), optional, target :: blck
real(psb_spk_), intent(inout) :: d(:) real(psb_spk_), intent(inout) :: d(:)
real(psb_spk_), intent(in), optional :: shft
end subroutine psb_siluk_fact end subroutine psb_siluk_fact
end interface end interface
interface psb_ilut_fact interface psb_ilut_fact
subroutine psb_silut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) subroutine psb_silut_fact(fill_in,thres,a,l,u,d,info,blck,iscale,shft)
import psb_sspmat_type, psb_spk_, psb_ipk_ import psb_sspmat_type, psb_spk_, psb_ipk_
integer(psb_ipk_), intent(in) :: fill_in integer(psb_ipk_), intent(in) :: fill_in
real(psb_spk_), intent(in) :: thres real(psb_spk_), intent(in) :: thres
@ -115,6 +117,7 @@ module psb_s_ilu_fact_mod
real(psb_spk_), intent(inout) :: d(:) real(psb_spk_), intent(inout) :: d(:)
type(psb_sspmat_type),intent(in), optional, target :: blck type(psb_sspmat_type),intent(in), optional, target :: blck
integer(psb_ipk_), intent(in), optional :: iscale integer(psb_ipk_), intent(in), optional :: iscale
real(psb_spk_), intent(in), optional :: shft
end subroutine psb_silut_fact end subroutine psb_silut_fact
end interface end interface

@ -80,7 +80,7 @@ module psb_z_ilu_fact_mod
use psb_base_mod use psb_base_mod
use psb_prec_const_mod use psb_prec_const_mod
interface psb_ilu0_fact interface psb_ilu0_fact
subroutine psb_zilu0_fact(ialg,a,l,u,d,info,blck,upd) subroutine psb_zilu0_fact(ialg,a,l,u,d,info,blck,upd,shft)
import psb_zspmat_type, psb_dpk_, psb_ipk_ import psb_zspmat_type, psb_dpk_, psb_ipk_
integer(psb_ipk_), intent(in) :: ialg integer(psb_ipk_), intent(in) :: ialg
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -89,11 +89,12 @@ module psb_z_ilu_fact_mod
type(psb_zspmat_type),intent(in), optional, target :: blck type(psb_zspmat_type),intent(in), optional, target :: blck
character, intent(in), optional :: upd character, intent(in), optional :: upd
complex(psb_dpk_), intent(inout) :: d(:) complex(psb_dpk_), intent(inout) :: d(:)
complex(psb_dpk_), intent(in), optional :: shft
end subroutine psb_zilu0_fact end subroutine psb_zilu0_fact
end interface end interface
interface psb_iluk_fact interface psb_iluk_fact
subroutine psb_ziluk_fact(fill_in,ialg,a,l,u,d,info,blck) subroutine psb_ziluk_fact(fill_in,ialg,a,l,u,d,info,blck,shft)
import psb_zspmat_type, psb_dpk_, psb_ipk_ import psb_zspmat_type, psb_dpk_, psb_ipk_
integer(psb_ipk_), intent(in) :: fill_in,ialg integer(psb_ipk_), intent(in) :: fill_in,ialg
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -101,11 +102,12 @@ module psb_z_ilu_fact_mod
type(psb_zspmat_type),intent(inout) :: l,u type(psb_zspmat_type),intent(inout) :: l,u
type(psb_zspmat_type),intent(in), optional, target :: blck type(psb_zspmat_type),intent(in), optional, target :: blck
complex(psb_dpk_), intent(inout) :: d(:) complex(psb_dpk_), intent(inout) :: d(:)
complex(psb_dpk_), intent(in), optional :: shft
end subroutine psb_ziluk_fact end subroutine psb_ziluk_fact
end interface end interface
interface psb_ilut_fact interface psb_ilut_fact
subroutine psb_zilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) subroutine psb_zilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale,shft)
import psb_zspmat_type, psb_dpk_, psb_ipk_ import psb_zspmat_type, psb_dpk_, psb_ipk_
integer(psb_ipk_), intent(in) :: fill_in integer(psb_ipk_), intent(in) :: fill_in
real(psb_dpk_), intent(in) :: thres real(psb_dpk_), intent(in) :: thres
@ -115,6 +117,7 @@ module psb_z_ilu_fact_mod
complex(psb_dpk_), intent(inout) :: d(:) complex(psb_dpk_), intent(inout) :: d(:)
type(psb_zspmat_type),intent(in), optional, target :: blck type(psb_zspmat_type),intent(in), optional, target :: blck
integer(psb_ipk_), intent(in), optional :: iscale integer(psb_ipk_), intent(in), optional :: iscale
complex(psb_dpk_), intent(in), optional :: shft
end subroutine psb_zilut_fact end subroutine psb_zilut_fact
end interface end interface

Loading…
Cancel
Save