diff --git a/prec/impl/psb_c_ilu0_fact.f90 b/prec/impl/psb_c_ilu0_fact.f90 index 1a3e1046..9e039b19 100644 --- a/prec/impl/psb_c_ilu0_fact.f90 +++ b/prec/impl/psb_c_ilu0_fact.f90 @@ -130,7 +130,7 @@ ! greater than 0. If the overlap is 0 or the matrix has been reordered ! (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_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 type(psb_cspmat_type),intent(in), optional, target :: blck character, intent(in), optional :: upd + complex(psb_spk_), intent(in), optional :: shft ! Local variables integer(psb_ipk_) :: l1, l2, m, err_act type(psb_cspmat_type), pointer :: blck_ type(psb_c_csr_sparse_mat) :: ll, uu + complex(psb_spk_) :: shft_ character :: upd_ character(len=20) :: name, ch_err @@ -177,7 +179,12 @@ subroutine psb_cilu0_fact(ialg,a,l,u,d,info,blck, upd) else upd_ = 'F' end if - + if (present(shft)) then + shft_ = shft + else + shft_ = czero + end if + m = a%get_nrows() + blck_%get_nrows() if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.& & (m > size(d)) ) then @@ -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 ! 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 info=psb_err_from_subroutine_ ch_err='psb_cilu0_factint' @@ -314,7 +321,7 @@ contains ! Error code. ! 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 @@ -325,6 +332,7 @@ contains integer(psb_ipk_), intent(inout) :: lja(:),lirp(:),uja(:),uirp(:) complex(psb_spk_), intent(inout) :: lval(:),uval(:),d(:) character, intent(in) :: upd + complex(psb_spk_), intent(in) :: shft ! Local variables 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 ! 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 ! ! Copy the i-th local row of the matrix, stored in b ! (as (i-ma)-th row), into lval/d(i)/uval ! 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 lirp(i+1) = l1 + 1 @@ -583,7 +591,7 @@ contains ! 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,& - & dia,l2,uja,uval,ktrw,trw,upd) + & dia,l2,uja,uval,ktrw,trw,upd,shft) use psb_base_mod @@ -597,6 +605,7 @@ contains integer(psb_ipk_), intent(inout) :: lja(:), uja(:) complex(psb_spk_), intent(inout) :: lval(:), uval(:), dia character, intent(in) :: upd + complex(psb_spk_), intent(in) :: shft ! Local variables integer(psb_ipk_) :: k,j,info,irb, nz integer(psb_ipk_), parameter :: nrb=40 @@ -625,7 +634,7 @@ contains lval(l1) = aa%val(j) lja(l1) = k else if (k == jd) then - dia = aa%val(j) + dia = aa%val(j) + shft else if ((k > jd).and.(k <= jmax)) then l2 = l2 + 1 uval(l2) = aa%val(j) @@ -665,7 +674,7 @@ contains lval(l1) = trw%val(ktrw) lja(l1) = k else if (k == jd) then - dia = trw%val(ktrw) + dia = trw%val(ktrw) + shft else if ((k > jd).and.(k <= jmax)) then l2 = l2 + 1 uval(l2) = trw%val(ktrw) diff --git a/prec/impl/psb_c_iluk_fact.f90 b/prec/impl/psb_c_iluk_fact.f90 index c4ebc678..ddda6d20 100644 --- a/prec/impl/psb_c_iluk_fact.f90 +++ b/prec/impl/psb_c_iluk_fact.f90 @@ -127,7 +127,7 @@ ! 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. ! -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_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(in), optional, target :: blck complex(psb_spk_), intent(inout) :: d(:) + complex(psb_spk_), intent(in), optional :: shft ! Local Variables 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 ! 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 info=psb_err_from_subroutine_ ch_err='psb_ciluk_factint' @@ -298,7 +299,7 @@ contains ! Error code. ! 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 @@ -311,6 +312,7 @@ contains integer(psb_ipk_), allocatable, intent(inout) :: lja(:),lirp(:),uja(:),uirp(:) complex(psb_spk_), allocatable, intent(inout) :: lval(:),uval(:) complex(psb_spk_), intent(inout) :: d(:) + complex(psb_spk_), intent(in) :: shft ! Local variables 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 ! - 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 ! ! Copy into trw the i-th local row of the matrix, stored in b ! (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 ! 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 ! 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 @@ -530,6 +532,8 @@ contains integer(psb_ipk_), intent(inout) :: rowlevs(:) complex(psb_spk_), intent(inout) :: row(:) type(psb_i_heap), intent(inout) :: heap + complex(psb_spk_), intent(in) :: shft + ! Local variables integer(psb_ipk_) :: k,j,irb,err_act,nz @@ -554,6 +558,7 @@ contains k = aa%ja(j) if ((jmin<=k).and.(k<=jmax)) then row(k) = aa%val(j) + if (k==i) row(k) = row(k) + shft rowlevs(k) = 0 call heap%insert(k,info) end if @@ -587,6 +592,7 @@ contains k = trw%ja(ktrw) if ((jmin<=k).and.(k<=jmax)) then row(k) = trw%val(ktrw) + if (k==i) row(k) = row(k) + shft rowlevs(k) = 0 call heap%insert(k,info) end if @@ -670,7 +676,8 @@ contains ! Note: this argument is intent(inout) and not only intent(out) ! 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 diff --git a/prec/impl/psb_c_ilut_fact.f90 b/prec/impl/psb_c_ilut_fact.f90 index 633899de..997b3b84 100644 --- a/prec/impl/psb_c_ilut_fact.f90 +++ b/prec/impl/psb_c_ilut_fact.f90 @@ -123,7 +123,7 @@ ! 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. ! -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_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(:) type(psb_cspmat_type),intent(in), optional, target :: blck integer(psb_ipk_), intent(in), optional :: iscale + complex(psb_spk_), intent(in), optional :: shft ! Local Variables 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 ! 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 info=psb_err_from_subroutine_ ch_err='psb_cilut_factint' @@ -316,7 +317,7 @@ contains ! Error code. ! 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 @@ -331,6 +332,7 @@ contains complex(psb_spk_), allocatable, intent(inout) :: lval(:),uval(:) complex(psb_spk_), intent(inout) :: d(:) real(psb_spk_), intent(in), optional :: scale + complex(psb_spk_), intent(in) :: shft ! Local Variables integer(psb_ipk_) :: i, ktrw,err_act,nidx,nlw,nup,jmaxup, ma, mb, m @@ -401,10 +403,10 @@ contains d(i) = czero if (i<=ma) then 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 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 ! @@ -540,7 +542,7 @@ contains ! 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,& - & nrmi,weight,row,heap,ktrw,trw,info) + & nrmi,weight,row,heap,ktrw,trw,info,shft) use psb_base_mod implicit none type(psb_cspmat_type), intent(in) :: a @@ -551,6 +553,7 @@ contains complex(psb_spk_), intent(inout) :: row(:) real(psb_spk_), intent(in) :: weight type(psb_i_heap), intent(inout) :: heap + complex(psb_spk_), intent(in) :: shft integer(psb_ipk_) :: k,j,irb,kin,nz integer(psb_ipk_), parameter :: nrb=40 @@ -597,6 +600,7 @@ contains call heap%insert(k,info) if (info /= psb_success_) exit if (kjd) then nup = nup + 1 if (abs(row(k))>dmaxup) then @@ -648,6 +652,7 @@ contains call heap%insert(k,info) if (info /= psb_success_) exit if (kjd) then nup = nup + 1 if (abs(row(k))>dmaxup) then diff --git a/prec/impl/psb_d_ilu0_fact.f90 b/prec/impl/psb_d_ilu0_fact.f90 index 478eedfa..29968e0c 100644 --- a/prec/impl/psb_d_ilu0_fact.f90 +++ b/prec/impl/psb_d_ilu0_fact.f90 @@ -130,7 +130,7 @@ ! greater than 0. If the overlap is 0 or the matrix has been reordered ! (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_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 type(psb_dspmat_type),intent(in), optional, target :: blck character, intent(in), optional :: upd + real(psb_dpk_), intent(in), optional :: shft ! Local variables integer(psb_ipk_) :: l1, l2, m, err_act type(psb_dspmat_type), pointer :: blck_ type(psb_d_csr_sparse_mat) :: ll, uu + real(psb_dpk_) :: shft_ character :: upd_ character(len=20) :: name, ch_err @@ -177,7 +179,12 @@ subroutine psb_dilu0_fact(ialg,a,l,u,d,info,blck, upd) else upd_ = 'F' end if - + if (present(shft)) then + shft_ = shft + else + shft_ = dzero + end if + m = a%get_nrows() + blck_%get_nrows() if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.& & (m > size(d)) ) then @@ -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 ! 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 info=psb_err_from_subroutine_ ch_err='psb_dilu0_factint' @@ -314,7 +321,7 @@ contains ! Error code. ! 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 @@ -325,6 +332,7 @@ contains integer(psb_ipk_), intent(inout) :: lja(:),lirp(:),uja(:),uirp(:) real(psb_dpk_), intent(inout) :: lval(:),uval(:),d(:) character, intent(in) :: upd + real(psb_dpk_), intent(in) :: shft ! Local variables 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 ! 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 ! ! Copy the i-th local row of the matrix, stored in b ! (as (i-ma)-th row), into lval/d(i)/uval ! 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 lirp(i+1) = l1 + 1 @@ -583,7 +591,7 @@ contains ! 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,& - & dia,l2,uja,uval,ktrw,trw,upd) + & dia,l2,uja,uval,ktrw,trw,upd,shft) use psb_base_mod @@ -597,6 +605,7 @@ contains integer(psb_ipk_), intent(inout) :: lja(:), uja(:) real(psb_dpk_), intent(inout) :: lval(:), uval(:), dia character, intent(in) :: upd + real(psb_dpk_), intent(in) :: shft ! Local variables integer(psb_ipk_) :: k,j,info,irb, nz integer(psb_ipk_), parameter :: nrb=40 @@ -625,7 +634,7 @@ contains lval(l1) = aa%val(j) lja(l1) = k else if (k == jd) then - dia = aa%val(j) + dia = aa%val(j) + shft else if ((k > jd).and.(k <= jmax)) then l2 = l2 + 1 uval(l2) = aa%val(j) @@ -665,7 +674,7 @@ contains lval(l1) = trw%val(ktrw) lja(l1) = k else if (k == jd) then - dia = trw%val(ktrw) + dia = trw%val(ktrw) + shft else if ((k > jd).and.(k <= jmax)) then l2 = l2 + 1 uval(l2) = trw%val(ktrw) diff --git a/prec/impl/psb_d_iluk_fact.f90 b/prec/impl/psb_d_iluk_fact.f90 index 544ec987..8f1e6c7b 100644 --- a/prec/impl/psb_d_iluk_fact.f90 +++ b/prec/impl/psb_d_iluk_fact.f90 @@ -127,7 +127,7 @@ ! 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. ! -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_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(in), optional, target :: blck real(psb_dpk_), intent(inout) :: d(:) + real(psb_dpk_), intent(in), optional :: shft ! Local Variables 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 ! 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 info=psb_err_from_subroutine_ ch_err='psb_diluk_factint' @@ -298,7 +299,7 @@ contains ! Error code. ! 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 @@ -311,6 +312,7 @@ contains integer(psb_ipk_), allocatable, intent(inout) :: lja(:),lirp(:),uja(:),uirp(:) real(psb_dpk_), allocatable, intent(inout) :: lval(:),uval(:) real(psb_dpk_), intent(inout) :: d(:) + real(psb_dpk_), intent(in) :: shft ! Local variables 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 ! - 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 ! ! Copy into trw the i-th local row of the matrix, stored in b ! (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 ! 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 ! 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 @@ -530,6 +532,8 @@ contains integer(psb_ipk_), intent(inout) :: rowlevs(:) real(psb_dpk_), intent(inout) :: row(:) type(psb_i_heap), intent(inout) :: heap + real(psb_dpk_), intent(in) :: shft + ! Local variables integer(psb_ipk_) :: k,j,irb,err_act,nz @@ -554,6 +558,7 @@ contains k = aa%ja(j) if ((jmin<=k).and.(k<=jmax)) then row(k) = aa%val(j) + if (k==i) row(k) = row(k) + shft rowlevs(k) = 0 call heap%insert(k,info) end if @@ -587,6 +592,7 @@ contains k = trw%ja(ktrw) if ((jmin<=k).and.(k<=jmax)) then row(k) = trw%val(ktrw) + if (k==i) row(k) = row(k) + shft rowlevs(k) = 0 call heap%insert(k,info) end if @@ -670,7 +676,8 @@ contains ! Note: this argument is intent(inout) and not only intent(out) ! 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 diff --git a/prec/impl/psb_d_ilut_fact.f90 b/prec/impl/psb_d_ilut_fact.f90 index 6c2dc698..c6079b55 100644 --- a/prec/impl/psb_d_ilut_fact.f90 +++ b/prec/impl/psb_d_ilut_fact.f90 @@ -123,7 +123,7 @@ ! 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. ! -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_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(:) type(psb_dspmat_type),intent(in), optional, target :: blck integer(psb_ipk_), intent(in), optional :: iscale + real(psb_dpk_), intent(in), optional :: shft ! Local Variables 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 ! 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 info=psb_err_from_subroutine_ ch_err='psb_dilut_factint' @@ -316,7 +317,7 @@ contains ! Error code. ! 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 @@ -331,6 +332,7 @@ contains real(psb_dpk_), allocatable, intent(inout) :: lval(:),uval(:) real(psb_dpk_), intent(inout) :: d(:) real(psb_dpk_), intent(in), optional :: scale + real(psb_dpk_), intent(in) :: shft ! Local Variables integer(psb_ipk_) :: i, ktrw,err_act,nidx,nlw,nup,jmaxup, ma, mb, m @@ -401,10 +403,10 @@ contains d(i) = czero if (i<=ma) then 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 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 ! @@ -540,7 +542,7 @@ contains ! 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,& - & nrmi,weight,row,heap,ktrw,trw,info) + & nrmi,weight,row,heap,ktrw,trw,info,shft) use psb_base_mod implicit none type(psb_dspmat_type), intent(in) :: a @@ -551,6 +553,7 @@ contains real(psb_dpk_), intent(inout) :: row(:) real(psb_dpk_), intent(in) :: weight type(psb_i_heap), intent(inout) :: heap + real(psb_dpk_), intent(in) :: shft integer(psb_ipk_) :: k,j,irb,kin,nz integer(psb_ipk_), parameter :: nrb=40 @@ -597,6 +600,7 @@ contains call heap%insert(k,info) if (info /= psb_success_) exit if (kjd) then nup = nup + 1 if (abs(row(k))>dmaxup) then @@ -648,6 +652,7 @@ contains call heap%insert(k,info) if (info /= psb_success_) exit if (kjd) then nup = nup + 1 if (abs(row(k))>dmaxup) then diff --git a/prec/impl/psb_s_ilu0_fact.f90 b/prec/impl/psb_s_ilu0_fact.f90 index b6f442e9..faa21fda 100644 --- a/prec/impl/psb_s_ilu0_fact.f90 +++ b/prec/impl/psb_s_ilu0_fact.f90 @@ -130,7 +130,7 @@ ! greater than 0. If the overlap is 0 or the matrix has been reordered ! (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_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 type(psb_sspmat_type),intent(in), optional, target :: blck character, intent(in), optional :: upd + real(psb_spk_), intent(in), optional :: shft ! Local variables integer(psb_ipk_) :: l1, l2, m, err_act type(psb_sspmat_type), pointer :: blck_ type(psb_s_csr_sparse_mat) :: ll, uu + real(psb_spk_) :: shft_ character :: upd_ character(len=20) :: name, ch_err @@ -177,7 +179,12 @@ subroutine psb_silu0_fact(ialg,a,l,u,d,info,blck, upd) else upd_ = 'F' end if - + if (present(shft)) then + shft_ = shft + else + shft_ = szero + end if + m = a%get_nrows() + blck_%get_nrows() if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.& & (m > size(d)) ) then @@ -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 ! 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 info=psb_err_from_subroutine_ ch_err='psb_silu0_factint' @@ -314,7 +321,7 @@ contains ! Error code. ! 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 @@ -325,6 +332,7 @@ contains integer(psb_ipk_), intent(inout) :: lja(:),lirp(:),uja(:),uirp(:) real(psb_spk_), intent(inout) :: lval(:),uval(:),d(:) character, intent(in) :: upd + real(psb_spk_), intent(in) :: shft ! Local variables 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 ! 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 ! ! Copy the i-th local row of the matrix, stored in b ! (as (i-ma)-th row), into lval/d(i)/uval ! 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 lirp(i+1) = l1 + 1 @@ -583,7 +591,7 @@ contains ! 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,& - & dia,l2,uja,uval,ktrw,trw,upd) + & dia,l2,uja,uval,ktrw,trw,upd,shft) use psb_base_mod @@ -597,6 +605,7 @@ contains integer(psb_ipk_), intent(inout) :: lja(:), uja(:) real(psb_spk_), intent(inout) :: lval(:), uval(:), dia character, intent(in) :: upd + real(psb_spk_), intent(in) :: shft ! Local variables integer(psb_ipk_) :: k,j,info,irb, nz integer(psb_ipk_), parameter :: nrb=40 @@ -625,7 +634,7 @@ contains lval(l1) = aa%val(j) lja(l1) = k else if (k == jd) then - dia = aa%val(j) + dia = aa%val(j) + shft else if ((k > jd).and.(k <= jmax)) then l2 = l2 + 1 uval(l2) = aa%val(j) @@ -665,7 +674,7 @@ contains lval(l1) = trw%val(ktrw) lja(l1) = k else if (k == jd) then - dia = trw%val(ktrw) + dia = trw%val(ktrw) + shft else if ((k > jd).and.(k <= jmax)) then l2 = l2 + 1 uval(l2) = trw%val(ktrw) diff --git a/prec/impl/psb_s_iluk_fact.f90 b/prec/impl/psb_s_iluk_fact.f90 index 6129663b..17ef4eef 100644 --- a/prec/impl/psb_s_iluk_fact.f90 +++ b/prec/impl/psb_s_iluk_fact.f90 @@ -127,7 +127,7 @@ ! 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. ! -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_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(in), optional, target :: blck real(psb_spk_), intent(inout) :: d(:) + real(psb_spk_), intent(in), optional :: shft ! Local Variables 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 ! 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 info=psb_err_from_subroutine_ ch_err='psb_siluk_factint' @@ -298,7 +299,7 @@ contains ! Error code. ! 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 @@ -311,6 +312,7 @@ contains integer(psb_ipk_), allocatable, intent(inout) :: lja(:),lirp(:),uja(:),uirp(:) real(psb_spk_), allocatable, intent(inout) :: lval(:),uval(:) real(psb_spk_), intent(inout) :: d(:) + real(psb_spk_), intent(in) :: shft ! Local variables 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 ! - 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 ! ! Copy into trw the i-th local row of the matrix, stored in b ! (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 ! 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 ! 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 @@ -530,6 +532,8 @@ contains integer(psb_ipk_), intent(inout) :: rowlevs(:) real(psb_spk_), intent(inout) :: row(:) type(psb_i_heap), intent(inout) :: heap + real(psb_spk_), intent(in) :: shft + ! Local variables integer(psb_ipk_) :: k,j,irb,err_act,nz @@ -554,6 +558,7 @@ contains k = aa%ja(j) if ((jmin<=k).and.(k<=jmax)) then row(k) = aa%val(j) + if (k==i) row(k) = row(k) + shft rowlevs(k) = 0 call heap%insert(k,info) end if @@ -587,6 +592,7 @@ contains k = trw%ja(ktrw) if ((jmin<=k).and.(k<=jmax)) then row(k) = trw%val(ktrw) + if (k==i) row(k) = row(k) + shft rowlevs(k) = 0 call heap%insert(k,info) end if @@ -670,7 +676,8 @@ contains ! Note: this argument is intent(inout) and not only intent(out) ! 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 diff --git a/prec/impl/psb_s_ilut_fact.f90 b/prec/impl/psb_s_ilut_fact.f90 index 43cacf41..76b514a3 100644 --- a/prec/impl/psb_s_ilut_fact.f90 +++ b/prec/impl/psb_s_ilut_fact.f90 @@ -123,7 +123,7 @@ ! 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. ! -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_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(:) type(psb_sspmat_type),intent(in), optional, target :: blck integer(psb_ipk_), intent(in), optional :: iscale + real(psb_spk_), intent(in), optional :: shft ! Local Variables 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 ! 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 info=psb_err_from_subroutine_ ch_err='psb_silut_factint' @@ -316,7 +317,7 @@ contains ! Error code. ! 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 @@ -331,6 +332,7 @@ contains real(psb_spk_), allocatable, intent(inout) :: lval(:),uval(:) real(psb_spk_), intent(inout) :: d(:) real(psb_spk_), intent(in), optional :: scale + real(psb_spk_), intent(in) :: shft ! Local Variables integer(psb_ipk_) :: i, ktrw,err_act,nidx,nlw,nup,jmaxup, ma, mb, m @@ -401,10 +403,10 @@ contains d(i) = czero if (i<=ma) then 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 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 ! @@ -540,7 +542,7 @@ contains ! 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,& - & nrmi,weight,row,heap,ktrw,trw,info) + & nrmi,weight,row,heap,ktrw,trw,info,shft) use psb_base_mod implicit none type(psb_sspmat_type), intent(in) :: a @@ -551,6 +553,7 @@ contains real(psb_spk_), intent(inout) :: row(:) real(psb_spk_), intent(in) :: weight type(psb_i_heap), intent(inout) :: heap + real(psb_spk_), intent(in) :: shft integer(psb_ipk_) :: k,j,irb,kin,nz integer(psb_ipk_), parameter :: nrb=40 @@ -597,6 +600,7 @@ contains call heap%insert(k,info) if (info /= psb_success_) exit if (kjd) then nup = nup + 1 if (abs(row(k))>dmaxup) then @@ -648,6 +652,7 @@ contains call heap%insert(k,info) if (info /= psb_success_) exit if (kjd) then nup = nup + 1 if (abs(row(k))>dmaxup) then diff --git a/prec/impl/psb_z_ilu0_fact.f90 b/prec/impl/psb_z_ilu0_fact.f90 index 26322e95..c6c0fc55 100644 --- a/prec/impl/psb_z_ilu0_fact.f90 +++ b/prec/impl/psb_z_ilu0_fact.f90 @@ -130,7 +130,7 @@ ! greater than 0. If the overlap is 0 or the matrix has been reordered ! (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_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 type(psb_zspmat_type),intent(in), optional, target :: blck character, intent(in), optional :: upd + complex(psb_dpk_), intent(in), optional :: shft ! Local variables integer(psb_ipk_) :: l1, l2, m, err_act type(psb_zspmat_type), pointer :: blck_ type(psb_z_csr_sparse_mat) :: ll, uu + complex(psb_dpk_) :: shft_ character :: upd_ character(len=20) :: name, ch_err @@ -177,7 +179,12 @@ subroutine psb_zilu0_fact(ialg,a,l,u,d,info,blck, upd) else upd_ = 'F' end if - + if (present(shft)) then + shft_ = shft + else + shft_ = zzero + end if + m = a%get_nrows() + blck_%get_nrows() if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.& & (m > size(d)) ) then @@ -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 ! 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 info=psb_err_from_subroutine_ ch_err='psb_zilu0_factint' @@ -314,7 +321,7 @@ contains ! Error code. ! 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 @@ -325,6 +332,7 @@ contains integer(psb_ipk_), intent(inout) :: lja(:),lirp(:),uja(:),uirp(:) complex(psb_dpk_), intent(inout) :: lval(:),uval(:),d(:) character, intent(in) :: upd + complex(psb_dpk_), intent(in) :: shft ! Local variables 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 ! 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 ! ! Copy the i-th local row of the matrix, stored in b ! (as (i-ma)-th row), into lval/d(i)/uval ! 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 lirp(i+1) = l1 + 1 @@ -583,7 +591,7 @@ contains ! 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,& - & dia,l2,uja,uval,ktrw,trw,upd) + & dia,l2,uja,uval,ktrw,trw,upd,shft) use psb_base_mod @@ -597,6 +605,7 @@ contains integer(psb_ipk_), intent(inout) :: lja(:), uja(:) complex(psb_dpk_), intent(inout) :: lval(:), uval(:), dia character, intent(in) :: upd + complex(psb_dpk_), intent(in) :: shft ! Local variables integer(psb_ipk_) :: k,j,info,irb, nz integer(psb_ipk_), parameter :: nrb=40 @@ -625,7 +634,7 @@ contains lval(l1) = aa%val(j) lja(l1) = k else if (k == jd) then - dia = aa%val(j) + dia = aa%val(j) + shft else if ((k > jd).and.(k <= jmax)) then l2 = l2 + 1 uval(l2) = aa%val(j) @@ -665,7 +674,7 @@ contains lval(l1) = trw%val(ktrw) lja(l1) = k else if (k == jd) then - dia = trw%val(ktrw) + dia = trw%val(ktrw) + shft else if ((k > jd).and.(k <= jmax)) then l2 = l2 + 1 uval(l2) = trw%val(ktrw) diff --git a/prec/impl/psb_z_iluk_fact.f90 b/prec/impl/psb_z_iluk_fact.f90 index 1a398cda..7675226a 100644 --- a/prec/impl/psb_z_iluk_fact.f90 +++ b/prec/impl/psb_z_iluk_fact.f90 @@ -127,7 +127,7 @@ ! 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. ! -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_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(in), optional, target :: blck complex(psb_dpk_), intent(inout) :: d(:) + complex(psb_dpk_), intent(in), optional :: shft ! Local Variables 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 ! 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 info=psb_err_from_subroutine_ ch_err='psb_ziluk_factint' @@ -298,7 +299,7 @@ contains ! Error code. ! 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 @@ -311,6 +312,7 @@ contains integer(psb_ipk_), allocatable, intent(inout) :: lja(:),lirp(:),uja(:),uirp(:) complex(psb_dpk_), allocatable, intent(inout) :: lval(:),uval(:) complex(psb_dpk_), intent(inout) :: d(:) + complex(psb_dpk_), intent(in) :: shft ! Local variables 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 ! - 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 ! ! Copy into trw the i-th local row of the matrix, stored in b ! (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 ! 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 ! 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 @@ -530,6 +532,8 @@ contains integer(psb_ipk_), intent(inout) :: rowlevs(:) complex(psb_dpk_), intent(inout) :: row(:) type(psb_i_heap), intent(inout) :: heap + complex(psb_dpk_), intent(in) :: shft + ! Local variables integer(psb_ipk_) :: k,j,irb,err_act,nz @@ -554,6 +558,7 @@ contains k = aa%ja(j) if ((jmin<=k).and.(k<=jmax)) then row(k) = aa%val(j) + if (k==i) row(k) = row(k) + shft rowlevs(k) = 0 call heap%insert(k,info) end if @@ -587,6 +592,7 @@ contains k = trw%ja(ktrw) if ((jmin<=k).and.(k<=jmax)) then row(k) = trw%val(ktrw) + if (k==i) row(k) = row(k) + shft rowlevs(k) = 0 call heap%insert(k,info) end if @@ -670,7 +676,8 @@ contains ! Note: this argument is intent(inout) and not only intent(out) ! 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 diff --git a/prec/impl/psb_z_ilut_fact.f90 b/prec/impl/psb_z_ilut_fact.f90 index 291dc778..2f004725 100644 --- a/prec/impl/psb_z_ilut_fact.f90 +++ b/prec/impl/psb_z_ilut_fact.f90 @@ -123,7 +123,7 @@ ! 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. ! -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_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(:) type(psb_zspmat_type),intent(in), optional, target :: blck integer(psb_ipk_), intent(in), optional :: iscale + complex(psb_dpk_), intent(in), optional :: shft ! Local Variables 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 ! 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 info=psb_err_from_subroutine_ ch_err='psb_zilut_factint' @@ -316,7 +317,7 @@ contains ! Error code. ! 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 @@ -331,6 +332,7 @@ contains complex(psb_dpk_), allocatable, intent(inout) :: lval(:),uval(:) complex(psb_dpk_), intent(inout) :: d(:) real(psb_dpk_), intent(in), optional :: scale + complex(psb_dpk_), intent(in) :: shft ! Local Variables integer(psb_ipk_) :: i, ktrw,err_act,nidx,nlw,nup,jmaxup, ma, mb, m @@ -401,10 +403,10 @@ contains d(i) = czero if (i<=ma) then 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 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 ! @@ -540,7 +542,7 @@ contains ! 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,& - & nrmi,weight,row,heap,ktrw,trw,info) + & nrmi,weight,row,heap,ktrw,trw,info,shft) use psb_base_mod implicit none type(psb_zspmat_type), intent(in) :: a @@ -551,6 +553,7 @@ contains complex(psb_dpk_), intent(inout) :: row(:) real(psb_dpk_), intent(in) :: weight type(psb_i_heap), intent(inout) :: heap + complex(psb_dpk_), intent(in) :: shft integer(psb_ipk_) :: k,j,irb,kin,nz integer(psb_ipk_), parameter :: nrb=40 @@ -597,6 +600,7 @@ contains call heap%insert(k,info) if (info /= psb_success_) exit if (kjd) then nup = nup + 1 if (abs(row(k))>dmaxup) then @@ -648,6 +652,7 @@ contains call heap%insert(k,info) if (info /= psb_success_) exit if (kjd) then nup = nup + 1 if (abs(row(k))>dmaxup) then diff --git a/prec/psb_c_ilu_fact_mod.f90 b/prec/psb_c_ilu_fact_mod.f90 index 45d06211..0fae1fc5 100644 --- a/prec/psb_c_ilu_fact_mod.f90 +++ b/prec/psb_c_ilu_fact_mod.f90 @@ -80,7 +80,7 @@ module psb_c_ilu_fact_mod use psb_base_mod use psb_prec_const_mod 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_ integer(psb_ipk_), intent(in) :: ialg 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 character, intent(in), optional :: upd complex(psb_spk_), intent(inout) :: d(:) + complex(psb_spk_), intent(in), optional :: shft end subroutine psb_cilu0_fact end interface 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_ integer(psb_ipk_), intent(in) :: fill_in,ialg 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(in), optional, target :: blck complex(psb_spk_), intent(inout) :: d(:) + complex(psb_spk_), intent(in), optional :: shft end subroutine psb_ciluk_fact end interface 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_ integer(psb_ipk_), intent(in) :: fill_in real(psb_spk_), intent(in) :: thres @@ -115,6 +117,7 @@ module psb_c_ilu_fact_mod complex(psb_spk_), intent(inout) :: d(:) type(psb_cspmat_type),intent(in), optional, target :: blck integer(psb_ipk_), intent(in), optional :: iscale + complex(psb_spk_), intent(in), optional :: shft end subroutine psb_cilut_fact end interface diff --git a/prec/psb_d_ilu_fact_mod.f90 b/prec/psb_d_ilu_fact_mod.f90 index 02753a4c..6354573d 100644 --- a/prec/psb_d_ilu_fact_mod.f90 +++ b/prec/psb_d_ilu_fact_mod.f90 @@ -80,7 +80,7 @@ module psb_d_ilu_fact_mod use psb_base_mod use psb_prec_const_mod 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_ integer(psb_ipk_), intent(in) :: ialg 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 character, intent(in), optional :: upd real(psb_dpk_), intent(inout) :: d(:) + real(psb_dpk_), intent(in), optional :: shft end subroutine psb_dilu0_fact end interface 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_ integer(psb_ipk_), intent(in) :: fill_in,ialg 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(in), optional, target :: blck real(psb_dpk_), intent(inout) :: d(:) + real(psb_dpk_), intent(in), optional :: shft end subroutine psb_diluk_fact end interface 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_ integer(psb_ipk_), intent(in) :: fill_in real(psb_dpk_), intent(in) :: thres @@ -115,6 +117,7 @@ module psb_d_ilu_fact_mod real(psb_dpk_), intent(inout) :: d(:) type(psb_dspmat_type),intent(in), optional, target :: blck integer(psb_ipk_), intent(in), optional :: iscale + real(psb_dpk_), intent(in), optional :: shft end subroutine psb_dilut_fact end interface diff --git a/prec/psb_s_ilu_fact_mod.f90 b/prec/psb_s_ilu_fact_mod.f90 index 6334df15..4021adc9 100644 --- a/prec/psb_s_ilu_fact_mod.f90 +++ b/prec/psb_s_ilu_fact_mod.f90 @@ -80,7 +80,7 @@ module psb_s_ilu_fact_mod use psb_base_mod use psb_prec_const_mod 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_ integer(psb_ipk_), intent(in) :: ialg 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 character, intent(in), optional :: upd real(psb_spk_), intent(inout) :: d(:) + real(psb_spk_), intent(in), optional :: shft end subroutine psb_silu0_fact end interface 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_ integer(psb_ipk_), intent(in) :: fill_in,ialg 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(in), optional, target :: blck real(psb_spk_), intent(inout) :: d(:) + real(psb_spk_), intent(in), optional :: shft end subroutine psb_siluk_fact end interface 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_ integer(psb_ipk_), intent(in) :: fill_in real(psb_spk_), intent(in) :: thres @@ -115,6 +117,7 @@ module psb_s_ilu_fact_mod real(psb_spk_), intent(inout) :: d(:) type(psb_sspmat_type),intent(in), optional, target :: blck integer(psb_ipk_), intent(in), optional :: iscale + real(psb_spk_), intent(in), optional :: shft end subroutine psb_silut_fact end interface diff --git a/prec/psb_z_ilu_fact_mod.f90 b/prec/psb_z_ilu_fact_mod.f90 index 220d673f..4793b43b 100644 --- a/prec/psb_z_ilu_fact_mod.f90 +++ b/prec/psb_z_ilu_fact_mod.f90 @@ -80,7 +80,7 @@ module psb_z_ilu_fact_mod use psb_base_mod use psb_prec_const_mod 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_ integer(psb_ipk_), intent(in) :: ialg 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 character, intent(in), optional :: upd complex(psb_dpk_), intent(inout) :: d(:) + complex(psb_dpk_), intent(in), optional :: shft end subroutine psb_zilu0_fact end interface 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_ integer(psb_ipk_), intent(in) :: fill_in,ialg 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(in), optional, target :: blck complex(psb_dpk_), intent(inout) :: d(:) + complex(psb_dpk_), intent(in), optional :: shft end subroutine psb_ziluk_fact end interface 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_ integer(psb_ipk_), intent(in) :: fill_in real(psb_dpk_), intent(in) :: thres @@ -115,6 +117,7 @@ module psb_z_ilu_fact_mod complex(psb_dpk_), intent(inout) :: d(:) type(psb_zspmat_type),intent(in), optional, target :: blck integer(psb_ipk_), intent(in), optional :: iscale + complex(psb_dpk_), intent(in), optional :: shft end subroutine psb_zilut_fact end interface