|
|
@ -95,7 +95,7 @@
|
|
|
|
subroutine mld_zilut_fact(fill_in,thres,a,l,u,d,info,blck)
|
|
|
|
subroutine mld_zilut_fact(fill_in,thres,a,l,u,d,info,blck)
|
|
|
|
|
|
|
|
|
|
|
|
use psb_sparse_mod
|
|
|
|
use psb_sparse_mod
|
|
|
|
use mld_inner_mod, mld_protect_name => mld_zilut_fact
|
|
|
|
use mld_inner_mod!, mld_protect_name => mld_zilut_fact
|
|
|
|
|
|
|
|
|
|
|
|
implicit none
|
|
|
|
implicit none
|
|
|
|
|
|
|
|
|
|
|
@ -112,6 +112,7 @@ subroutine mld_zilut_fact(fill_in,thres,a,l,u,d,info,blck)
|
|
|
|
integer :: l1, l2, m, err_act
|
|
|
|
integer :: 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
|
|
|
|
character(len=20) :: name, ch_err
|
|
|
|
character(len=20) :: name, ch_err
|
|
|
|
|
|
|
|
|
|
|
|
name='mld_zilut_fact'
|
|
|
|
name='mld_zilut_fact'
|
|
|
@ -130,26 +131,32 @@ subroutine mld_zilut_fact(fill_in,thres,a,l,u,d,info,blck)
|
|
|
|
blck_ => blck
|
|
|
|
blck_ => blck
|
|
|
|
else
|
|
|
|
else
|
|
|
|
allocate(blck_,stat=info)
|
|
|
|
allocate(blck_,stat=info)
|
|
|
|
if (info /= psb_success_) then
|
|
|
|
if (info == psb_success_) call blck_%csall(0,0,info,1)
|
|
|
|
call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
|
|
|
|
|
|
|
|
goto 9999
|
|
|
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
call psb_sp_all(0,0,blck_,1,info)
|
|
|
|
|
|
|
|
if (info /= psb_success_) then
|
|
|
|
if (info /= psb_success_) then
|
|
|
|
info=psb_err_from_subroutine_
|
|
|
|
info=psb_err_from_subroutine_
|
|
|
|
ch_err='psb_sp_all'
|
|
|
|
ch_err='csall'
|
|
|
|
call psb_errpush(info,name,a_err=ch_err)
|
|
|
|
call psb_errpush(info,name,a_err=ch_err)
|
|
|
|
goto 9999
|
|
|
|
goto 9999
|
|
|
|
end if
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
|
|
endif
|
|
|
|
endif
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
m = a%get_nrows() + blck_%get_nrows()
|
|
|
|
|
|
|
|
if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.&
|
|
|
|
|
|
|
|
& (m > size(d)) ) then
|
|
|
|
|
|
|
|
write(0,*) 'Wrong allocation status for L,D,U? ',&
|
|
|
|
|
|
|
|
& l%get_nrows(),size(d),u%get_nrows()
|
|
|
|
|
|
|
|
info = -1
|
|
|
|
|
|
|
|
return
|
|
|
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
call l%mv_to(ll)
|
|
|
|
|
|
|
|
call u%mv_to(uu)
|
|
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! Compute the ILU(k,t) factorization
|
|
|
|
! Compute the ILU(k,t) factorization
|
|
|
|
!
|
|
|
|
!
|
|
|
|
call mld_zilut_factint(fill_in,thres,m,a,blck_,&
|
|
|
|
call mld_zilut_factint(fill_in,thres,a,blck_,&
|
|
|
|
& d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info)
|
|
|
|
& d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,info)
|
|
|
|
if (info /= psb_success_) then
|
|
|
|
if (info /= psb_success_) then
|
|
|
|
info=psb_err_from_subroutine_
|
|
|
|
info=psb_err_from_subroutine_
|
|
|
|
ch_err='mld_zilut_factint'
|
|
|
|
ch_err='mld_zilut_factint'
|
|
|
@ -160,31 +167,29 @@ subroutine mld_zilut_fact(fill_in,thres,a,l,u,d,info,blck)
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! Store information on the L and U sparse matrices
|
|
|
|
! Store information on the L and U sparse matrices
|
|
|
|
!
|
|
|
|
!
|
|
|
|
l%infoa(1) = l1
|
|
|
|
call l%mv_from(ll)
|
|
|
|
l%fida = 'CSR'
|
|
|
|
call l%set_triangle()
|
|
|
|
l%descra = 'TLU'
|
|
|
|
call l%set_unit()
|
|
|
|
u%infoa(1) = l2
|
|
|
|
call l%set_lower()
|
|
|
|
u%fida = 'CSR'
|
|
|
|
call u%mv_from(uu)
|
|
|
|
u%descra = 'TUU'
|
|
|
|
call u%set_triangle()
|
|
|
|
l%m = m
|
|
|
|
call u%set_unit()
|
|
|
|
l%k = m
|
|
|
|
call u%set_upper()
|
|
|
|
u%m = m
|
|
|
|
|
|
|
|
u%k = m
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! Nullify the pointer / deallocate the memory
|
|
|
|
! Nullify pointer / deallocate memory
|
|
|
|
!
|
|
|
|
!
|
|
|
|
if (present(blck)) then
|
|
|
|
if (present(blck)) then
|
|
|
|
blck_ => null()
|
|
|
|
blck_ => null()
|
|
|
|
else
|
|
|
|
else
|
|
|
|
call psb_sp_free(blck_,info)
|
|
|
|
call blck_%free()
|
|
|
|
if (info /= psb_success_) then
|
|
|
|
deallocate(blck_,stat=info)
|
|
|
|
info=psb_err_from_subroutine_
|
|
|
|
if(info.ne.0) then
|
|
|
|
ch_err='psb_sp_free'
|
|
|
|
info=psb_err_from_subroutine_
|
|
|
|
call psb_errpush(info,name,a_err=ch_err)
|
|
|
|
ch_err='psb_sp_free'
|
|
|
|
goto 9999
|
|
|
|
call psb_errpush(info,name,a_err=ch_err)
|
|
|
|
|
|
|
|
goto 9999
|
|
|
|
end if
|
|
|
|
end if
|
|
|
|
deallocate(blck_)
|
|
|
|
|
|
|
|
endif
|
|
|
|
endif
|
|
|
|
|
|
|
|
|
|
|
|
call psb_erractionrestore(err_act)
|
|
|
|
call psb_erractionrestore(err_act)
|
|
|
@ -241,53 +246,53 @@ contains
|
|
|
|
! d - complex(psb_dpk_), dimension(:), output.
|
|
|
|
! d - complex(psb_dpk_), dimension(:), output.
|
|
|
|
! The inverse of the diagonal entries of the U factor in the incomplete
|
|
|
|
! The inverse of the diagonal entries of the U factor in the incomplete
|
|
|
|
! factorization.
|
|
|
|
! factorization.
|
|
|
|
! laspk - complex(psb_dpk_), dimension(:), input/output.
|
|
|
|
! lval - complex(psb_dpk_), dimension(:), input/output.
|
|
|
|
! The L factor in the incomplete factorization.
|
|
|
|
! The L factor in the incomplete factorization.
|
|
|
|
! lia1 - integer, dimension(:), input/output.
|
|
|
|
! lia1 - integer, dimension(:), input/output.
|
|
|
|
! The column indices of the nonzero entries of the L factor,
|
|
|
|
! The column indices of the nonzero entries of the L factor,
|
|
|
|
! according to the CSR storage format.
|
|
|
|
! according to the CSR storage format.
|
|
|
|
! lia2 - integer, dimension(:), input/output.
|
|
|
|
! lirp - integer, dimension(:), input/output.
|
|
|
|
! The indices identifying the first nonzero entry of each row
|
|
|
|
! The indices identifying the first nonzero entry of each row
|
|
|
|
! of the L factor in laspk, according to the CSR storage format.
|
|
|
|
! of the L factor in lval, according to the CSR storage format.
|
|
|
|
! uaspk - complex(psb_dpk_), dimension(:), input/output.
|
|
|
|
! uval - complex(psb_dpk_), dimension(:), input/output.
|
|
|
|
! The U factor in the incomplete factorization.
|
|
|
|
! The U factor in the incomplete factorization.
|
|
|
|
! The entries of U are stored according to the CSR format.
|
|
|
|
! The entries of U are stored according to the CSR format.
|
|
|
|
! uia1 - integer, dimension(:), input/output.
|
|
|
|
! uja - integer, dimension(:), input/output.
|
|
|
|
! The column indices of the nonzero entries of the U factor,
|
|
|
|
! The column indices of the nonzero entries of the U factor,
|
|
|
|
! according to the CSR storage format.
|
|
|
|
! according to the CSR storage format.
|
|
|
|
! uia2 - integer, dimension(:), input/output.
|
|
|
|
! uirp - integer, dimension(:), input/output.
|
|
|
|
! The indices identifying the first nonzero entry of each row
|
|
|
|
! The indices identifying the first nonzero entry of each row
|
|
|
|
! of the U factor in uaspk, according to the CSR storage format.
|
|
|
|
! of the U factor in uval, according to the CSR storage format.
|
|
|
|
! l1 - integer, output
|
|
|
|
! l1 - integer, output
|
|
|
|
! The number of nonzero entries in laspk.
|
|
|
|
! The number of nonzero entries in lval.
|
|
|
|
! l2 - integer, output
|
|
|
|
! l2 - integer, output
|
|
|
|
! The number of nonzero entries in uaspk.
|
|
|
|
! The number of nonzero entries in uval.
|
|
|
|
! info - integer, output.
|
|
|
|
! info - integer, output.
|
|
|
|
! Error code.
|
|
|
|
! Error code.
|
|
|
|
!
|
|
|
|
!
|
|
|
|
subroutine mld_zilut_factint(fill_in,thres,m,a,b,&
|
|
|
|
subroutine mld_zilut_factint(fill_in,thres,a,b,&
|
|
|
|
& d,laspk,lia1,lia2,uaspk,uia1,uia2,l1,l2,info)
|
|
|
|
& d,lval,lja,lirp,uval,uja,uirp,l1,l2,info)
|
|
|
|
|
|
|
|
|
|
|
|
use psb_sparse_mod
|
|
|
|
use psb_sparse_mod
|
|
|
|
|
|
|
|
|
|
|
|
implicit none
|
|
|
|
implicit none
|
|
|
|
|
|
|
|
|
|
|
|
! Arguments
|
|
|
|
! Arguments
|
|
|
|
integer, intent(in) :: fill_in
|
|
|
|
integer, intent(in) :: fill_in
|
|
|
|
real(psb_dpk_), intent(in) :: thres
|
|
|
|
real(psb_dpk_), intent(in) :: thres
|
|
|
|
type(psb_zspmat_type), intent(in) :: a,b
|
|
|
|
type(psb_zspmat_type),intent(in) :: a,b
|
|
|
|
integer, intent(inout) :: m,l1,l2,info
|
|
|
|
integer,intent(inout) :: l1,l2,info
|
|
|
|
integer, allocatable, intent(inout) :: lia1(:),lia2(:),uia1(:),uia2(:)
|
|
|
|
integer, allocatable, intent(inout) :: lja(:),lirp(:),uja(:),uirp(:)
|
|
|
|
complex(psb_dpk_), allocatable, intent(inout) :: laspk(:),uaspk(:)
|
|
|
|
complex(psb_dpk_), allocatable, intent(inout) :: lval(:),uval(:)
|
|
|
|
complex(psb_dpk_), intent(inout) :: d(:)
|
|
|
|
complex(psb_dpk_), intent(inout) :: d(:)
|
|
|
|
|
|
|
|
|
|
|
|
! Local Variables
|
|
|
|
! Local Variables
|
|
|
|
integer :: i, ktrw,err_act,nidx,nlw,nup,jmaxup, ma, mb
|
|
|
|
integer :: i, ktrw,err_act,nidx,nlw,nup,jmaxup, ma, mb, m
|
|
|
|
real(psb_dpk_) :: nrmi
|
|
|
|
real(psb_dpk_) :: nrmi
|
|
|
|
integer, allocatable :: idxs(:)
|
|
|
|
integer, allocatable :: idxs(:)
|
|
|
|
complex(psb_dpk_), allocatable :: row(:)
|
|
|
|
complex(psb_dpk_), allocatable :: row(:)
|
|
|
|
type(psb_int_heap) :: heap
|
|
|
|
type(psb_int_heap) :: heap
|
|
|
|
type(psb_zspmat_type) :: trw
|
|
|
|
type(psb_z_coo_sparse_mat) :: trw
|
|
|
|
character(len=20), parameter :: name='mld_zilut_factint'
|
|
|
|
character(len=20), parameter :: name='mld_zilut_factint'
|
|
|
|
character(len=20) :: ch_err
|
|
|
|
character(len=20) :: ch_err
|
|
|
|
|
|
|
|
|
|
|
@ -296,16 +301,16 @@ contains
|
|
|
|
call psb_erractionsave(err_act)
|
|
|
|
call psb_erractionsave(err_act)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
ma = a%m
|
|
|
|
ma = a%get_nrows()
|
|
|
|
mb = b%m
|
|
|
|
mb = b%get_nrows()
|
|
|
|
m = ma+mb
|
|
|
|
m = ma+mb
|
|
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! Allocate a temporary buffer for the ilut_copyin function
|
|
|
|
! Allocate a temporary buffer for the ilut_copyin function
|
|
|
|
!
|
|
|
|
!
|
|
|
|
call psb_sp_all(0,0,trw,1,info)
|
|
|
|
call trw%allocate(0,0,1)
|
|
|
|
if (info == psb_success_) call psb_ensure_size(m+1,lia2,info)
|
|
|
|
if (info == psb_success_) call psb_ensure_size(m+1,lirp,info)
|
|
|
|
if (info == psb_success_) call psb_ensure_size(m+1,uia2,info)
|
|
|
|
if (info == psb_success_) call psb_ensure_size(m+1,uirp,info)
|
|
|
|
|
|
|
|
|
|
|
|
if (info /= psb_success_) then
|
|
|
|
if (info /= psb_success_) then
|
|
|
|
info=psb_err_from_subroutine_
|
|
|
|
info=psb_err_from_subroutine_
|
|
|
@ -315,8 +320,8 @@ contains
|
|
|
|
|
|
|
|
|
|
|
|
l1=0
|
|
|
|
l1=0
|
|
|
|
l2=0
|
|
|
|
l2=0
|
|
|
|
lia2(1) = 1
|
|
|
|
lirp(1) = 1
|
|
|
|
uia2(1) = 1
|
|
|
|
uirp(1) = 1
|
|
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! Allocate memory to hold the entries of a row
|
|
|
|
! Allocate memory to hold the entries of a row
|
|
|
@ -354,12 +359,12 @@ contains
|
|
|
|
! Do an elimination step on current row
|
|
|
|
! Do an elimination step on current row
|
|
|
|
!
|
|
|
|
!
|
|
|
|
if (info == psb_success_) call ilut_fact(thres,i,nrmi,row,heap,&
|
|
|
|
if (info == psb_success_) call ilut_fact(thres,i,nrmi,row,heap,&
|
|
|
|
& d,uia1,uia2,uaspk,nidx,idxs,info)
|
|
|
|
& d,uja,uirp,uval,nidx,idxs,info)
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! Copy the row into laspk/d(i)/uaspk
|
|
|
|
! Copy the row into lval/d(i)/uval
|
|
|
|
!
|
|
|
|
!
|
|
|
|
if (info == psb_success_) call ilut_copyout(fill_in,thres,i,m,nlw,nup,jmaxup,nrmi,row,nidx,idxs,&
|
|
|
|
if (info == psb_success_) call ilut_copyout(fill_in,thres,i,m,nlw,nup,jmaxup,nrmi,row,nidx,idxs,&
|
|
|
|
& l1,l2,lia1,lia2,laspk,d,uia1,uia2,uaspk,info)
|
|
|
|
& l1,l2,lja,lirp,lval,d,uja,uirp,uval,info)
|
|
|
|
|
|
|
|
|
|
|
|
if (info /= psb_success_) then
|
|
|
|
if (info /= psb_success_) then
|
|
|
|
info=psb_err_internal_error_
|
|
|
|
info=psb_err_internal_error_
|
|
|
@ -378,7 +383,7 @@ contains
|
|
|
|
call psb_errpush(info,name,a_err='Deallocate')
|
|
|
|
call psb_errpush(info,name,a_err='Deallocate')
|
|
|
|
goto 9999
|
|
|
|
goto 9999
|
|
|
|
end if
|
|
|
|
end if
|
|
|
|
if (info == psb_success_) call psb_sp_free(trw,info)
|
|
|
|
if (info == psb_success_) call trw%free()
|
|
|
|
if (info /= psb_success_) then
|
|
|
|
if (info /= psb_success_) then
|
|
|
|
info=psb_err_from_subroutine_
|
|
|
|
info=psb_err_from_subroutine_
|
|
|
|
ch_err='psb_sp_free'
|
|
|
|
ch_err='psb_sp_free'
|
|
|
@ -482,17 +487,17 @@ contains
|
|
|
|
subroutine ilut_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,nrmi,row,heap,ktrw,trw,info)
|
|
|
|
subroutine ilut_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,nrmi,row,heap,ktrw,trw,info)
|
|
|
|
use psb_sparse_mod
|
|
|
|
use psb_sparse_mod
|
|
|
|
implicit none
|
|
|
|
implicit none
|
|
|
|
type(psb_zspmat_type), intent(in) :: a
|
|
|
|
type(psb_zspmat_type), intent(in) :: a
|
|
|
|
type(psb_zspmat_type), intent(inout) :: trw
|
|
|
|
type(psb_z_coo_sparse_mat), intent(inout) :: trw
|
|
|
|
integer, intent(in) :: i, m,jmin,jmax,jd
|
|
|
|
integer, intent(in) :: i, m,jmin,jmax,jd
|
|
|
|
integer, intent(inout) :: ktrw,nlw,nup,jmaxup,info
|
|
|
|
integer, intent(inout) :: ktrw,nlw,nup,jmaxup,info
|
|
|
|
real(psb_dpk_), intent(inout) :: nrmi
|
|
|
|
real(psb_dpk_), intent(inout) :: nrmi
|
|
|
|
complex(psb_dpk_), intent(inout) :: row(:)
|
|
|
|
complex(psb_dpk_), intent(inout) :: row(:)
|
|
|
|
type(psb_int_heap), intent(inout) :: heap
|
|
|
|
type(psb_int_heap), intent(inout) :: heap
|
|
|
|
|
|
|
|
|
|
|
|
integer :: k,j,irb,kin,nz
|
|
|
|
integer :: k,j,irb,kin,nz
|
|
|
|
integer, parameter :: nrb=16
|
|
|
|
integer, parameter :: nrb=40
|
|
|
|
real(psb_dpk_) :: dmaxup
|
|
|
|
real(psb_dpk_) :: dmaxup
|
|
|
|
real(psb_dpk_), external :: dznrm2
|
|
|
|
real(psb_dpk_), external :: dznrm2
|
|
|
|
character(len=20), parameter :: name='mld_zilut_factint'
|
|
|
|
character(len=20), parameter :: name='mld_zilut_factint'
|
|
|
|
|
|
|
|
|
|
|
@ -517,25 +522,21 @@ contains
|
|
|
|
nlw = 0
|
|
|
|
nlw = 0
|
|
|
|
nup = 0
|
|
|
|
nup = 0
|
|
|
|
jmaxup = 0
|
|
|
|
jmaxup = 0
|
|
|
|
dmaxup = dzero
|
|
|
|
dmaxup = szero
|
|
|
|
nrmi = dzero
|
|
|
|
nrmi = szero
|
|
|
|
|
|
|
|
|
|
|
|
if (psb_toupper(a%fida) == 'CSR') then
|
|
|
|
select type (aa=> a%a)
|
|
|
|
|
|
|
|
type is (psb_z_csr_sparse_mat)
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! Take a fast shortcut if the matrix is stored in CSR format
|
|
|
|
! Take a fast shortcut if the matrix is stored in CSR format
|
|
|
|
!
|
|
|
|
!
|
|
|
|
|
|
|
|
|
|
|
|
do j = a%ia2(i), a%ia2(i+1) - 1
|
|
|
|
do j = aa%irp(i), aa%irp(i+1) - 1
|
|
|
|
k = a%ia1(j)
|
|
|
|
k = aa%ja(j)
|
|
|
|
if ((jmin<=k).and.(k<=jmax)) then
|
|
|
|
if ((jmin<=k).and.(k<=jmax)) then
|
|
|
|
row(k) = a%aspk(j)
|
|
|
|
row(k) = aa%val(j)
|
|
|
|
call psb_insert_heap(k,heap,info)
|
|
|
|
call psb_insert_heap(k,heap,info)
|
|
|
|
if (info /= psb_success_) then
|
|
|
|
if (info /= psb_success_) exit
|
|
|
|
info=psb_err_from_subroutine_
|
|
|
|
|
|
|
|
call psb_errpush(info,name,a_err='psb_insert_heap')
|
|
|
|
|
|
|
|
goto 9999
|
|
|
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
end if
|
|
|
|
end if
|
|
|
|
if (k<jd) nlw = nlw + 1
|
|
|
|
if (k<jd) nlw = nlw + 1
|
|
|
|
if (k>jd) then
|
|
|
|
if (k>jd) then
|
|
|
@ -546,9 +547,17 @@ contains
|
|
|
|
end if
|
|
|
|
end if
|
|
|
|
end if
|
|
|
|
end if
|
|
|
|
end do
|
|
|
|
end do
|
|
|
|
nz = a%ia2(i+1) - a%ia2(i)
|
|
|
|
if (info /= psb_success_) then
|
|
|
|
nrmi = dznrm2(nz,a%aspk(a%ia2(i)),ione)
|
|
|
|
info=psb_err_from_subroutine_
|
|
|
|
else
|
|
|
|
call psb_errpush(info,name,a_err='psb_insert_heap')
|
|
|
|
|
|
|
|
goto 9999
|
|
|
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
nz = aa%irp(i+1) - aa%irp(i)
|
|
|
|
|
|
|
|
nrmi = dznrm2(nz,aa%val(aa%irp(i)),ione)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
class default
|
|
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! Otherwise use psb_sp_getblk, slower but able (in principle) of
|
|
|
|
! Otherwise use psb_sp_getblk, slower but able (in principle) of
|
|
|
@ -560,7 +569,7 @@ contains
|
|
|
|
|
|
|
|
|
|
|
|
if ((mod(i,nrb) == 1).or.(nrb == 1)) then
|
|
|
|
if ((mod(i,nrb) == 1).or.(nrb == 1)) then
|
|
|
|
irb = min(m-i+1,nrb)
|
|
|
|
irb = min(m-i+1,nrb)
|
|
|
|
call psb_sp_getblk(i,a,trw,info,lrw=i+irb-1)
|
|
|
|
call aa%csget(i,i+irb-1,trw,info)
|
|
|
|
if (info /= psb_success_) then
|
|
|
|
if (info /= psb_success_) then
|
|
|
|
info=psb_err_from_subroutine_
|
|
|
|
info=psb_err_from_subroutine_
|
|
|
|
call psb_errpush(info,name,a_err='psb_sp_getblk')
|
|
|
|
call psb_errpush(info,name,a_err='psb_sp_getblk')
|
|
|
@ -570,18 +579,16 @@ contains
|
|
|
|
end if
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
|
|
kin = ktrw
|
|
|
|
kin = ktrw
|
|
|
|
|
|
|
|
nz = trw%get_nzeros()
|
|
|
|
do
|
|
|
|
do
|
|
|
|
if (ktrw > trw%infoa(psb_nnz_)) exit
|
|
|
|
if (ktrw > nz) exit
|
|
|
|
if (trw%ia1(ktrw) > i) exit
|
|
|
|
if (trw%ia(ktrw) > i) exit
|
|
|
|
k = trw%ia2(ktrw)
|
|
|
|
k = trw%ja(ktrw)
|
|
|
|
if ((jmin<=k).and.(k<=jmax)) then
|
|
|
|
if ((jmin<=k).and.(k<=jmax)) then
|
|
|
|
row(k) = trw%aspk(ktrw)
|
|
|
|
row(k) = trw%val(ktrw)
|
|
|
|
call psb_insert_heap(k,heap,info)
|
|
|
|
call psb_insert_heap(k,heap,info)
|
|
|
|
if (info /= psb_success_) then
|
|
|
|
if (info /= psb_success_) exit
|
|
|
|
info=psb_err_from_subroutine_
|
|
|
|
|
|
|
|
call psb_errpush(info,name,a_err='psb_insert_heap')
|
|
|
|
|
|
|
|
goto 9999
|
|
|
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
end if
|
|
|
|
end if
|
|
|
|
if (k<jd) nlw = nlw + 1
|
|
|
|
if (k<jd) nlw = nlw + 1
|
|
|
|
if (k>jd) then
|
|
|
|
if (k>jd) then
|
|
|
@ -594,8 +601,9 @@ contains
|
|
|
|
ktrw = ktrw + 1
|
|
|
|
ktrw = ktrw + 1
|
|
|
|
enddo
|
|
|
|
enddo
|
|
|
|
nz = ktrw - kin
|
|
|
|
nz = ktrw - kin
|
|
|
|
nrmi = dznrm2(nz,trw%aspk(kin),ione)
|
|
|
|
nrmi = dznrm2(nz,trw%val(kin),ione)
|
|
|
|
end if
|
|
|
|
end select
|
|
|
|
|
|
|
|
|
|
|
|
call psb_erractionrestore(err_act)
|
|
|
|
call psb_erractionrestore(err_act)
|
|
|
|
return
|
|
|
|
return
|
|
|
|
|
|
|
|
|
|
|
@ -645,17 +653,17 @@ contains
|
|
|
|
! d - complex(psb_dpk_), input.
|
|
|
|
! d - complex(psb_dpk_), input.
|
|
|
|
! The inverse of the diagonal entries of the part of the U factor
|
|
|
|
! The inverse of the diagonal entries of the part of the U factor
|
|
|
|
! above the current row (see ilut_copyout).
|
|
|
|
! above the current row (see ilut_copyout).
|
|
|
|
! uia1 - integer, dimension(:), input.
|
|
|
|
! uja - integer, dimension(:), input.
|
|
|
|
! The column indices of the nonzero entries of the part of the U
|
|
|
|
! The column indices of the nonzero entries of the part of the U
|
|
|
|
! factor above the current row, stored in uaspk row by row (see
|
|
|
|
! factor above the current row, stored in uval row by row (see
|
|
|
|
! ilut_copyout, called by mld_zilut_factint), according to the CSR
|
|
|
|
! ilut_copyout, called by mld_zilut_factint), according to the CSR
|
|
|
|
! storage format.
|
|
|
|
! storage format.
|
|
|
|
! uia2 - integer, dimension(:), input.
|
|
|
|
! uirp - integer, dimension(:), input.
|
|
|
|
! The indices identifying the first nonzero entry of each row of
|
|
|
|
! The indices identifying the first nonzero entry of each row of
|
|
|
|
! the U factor above the current row, stored in uaspk row by row
|
|
|
|
! the U factor above the current row, stored in uval row by row
|
|
|
|
! (see ilut_copyout, called by mld_zilut_factint), according to
|
|
|
|
! (see ilut_copyout, called by mld_zilut_factint), according to
|
|
|
|
! the CSR storage format.
|
|
|
|
! the CSR storage format.
|
|
|
|
! uaspk - complex(psb_dpk_), dimension(:), input.
|
|
|
|
! uval - complex(psb_dpk_), dimension(:), input.
|
|
|
|
! The entries of the U factor above the current row (except the
|
|
|
|
! The entries of the U factor above the current row (except the
|
|
|
|
! diagonal ones), stored according to the CSR format.
|
|
|
|
! diagonal ones), stored according to the CSR format.
|
|
|
|
! nidx - integer, output.
|
|
|
|
! nidx - integer, output.
|
|
|
@ -669,7 +677,7 @@ 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 ilut_fact(thres,i,nrmi,row,heap,d,uia1,uia2,uaspk,nidx,idxs,info)
|
|
|
|
subroutine ilut_fact(thres,i,nrmi,row,heap,d,uja,uirp,uval,nidx,idxs,info)
|
|
|
|
|
|
|
|
|
|
|
|
use psb_sparse_mod
|
|
|
|
use psb_sparse_mod
|
|
|
|
|
|
|
|
|
|
|
@ -681,8 +689,8 @@ contains
|
|
|
|
integer, intent(inout) :: nidx,info
|
|
|
|
integer, intent(inout) :: nidx,info
|
|
|
|
real(psb_dpk_), intent(in) :: thres,nrmi
|
|
|
|
real(psb_dpk_), intent(in) :: thres,nrmi
|
|
|
|
integer, allocatable, intent(inout) :: idxs(:)
|
|
|
|
integer, allocatable, intent(inout) :: idxs(:)
|
|
|
|
integer, intent(inout) :: uia1(:),uia2(:)
|
|
|
|
integer, intent(inout) :: uja(:),uirp(:)
|
|
|
|
complex(psb_dpk_), intent(inout) :: row(:), uaspk(:),d(:)
|
|
|
|
complex(psb_dpk_), intent(inout) :: row(:), uval(:),d(:)
|
|
|
|
|
|
|
|
|
|
|
|
! Local Variables
|
|
|
|
! Local Variables
|
|
|
|
integer :: k,j,jj,lastk, iret
|
|
|
|
integer :: k,j,jj,lastk, iret
|
|
|
@ -726,8 +734,8 @@ contains
|
|
|
|
! Note: since U is scaled while copying it out (see ilut_copyout),
|
|
|
|
! Note: since U is scaled while copying it out (see ilut_copyout),
|
|
|
|
! we can use rwk in the update below.
|
|
|
|
! we can use rwk in the update below.
|
|
|
|
!
|
|
|
|
!
|
|
|
|
do jj=uia2(k),uia2(k+1)-1
|
|
|
|
do jj=uirp(k),uirp(k+1)-1
|
|
|
|
j = uia1(jj)
|
|
|
|
j = uja(jj)
|
|
|
|
if (j<=k) then
|
|
|
|
if (j<=k) then
|
|
|
|
info = -i
|
|
|
|
info = -i
|
|
|
|
return
|
|
|
|
return
|
|
|
@ -736,7 +744,7 @@ contains
|
|
|
|
! Update row(j) and, if it is not to be discarded, insert
|
|
|
|
! Update row(j) and, if it is not to be discarded, insert
|
|
|
|
! its index into the heap for further processing.
|
|
|
|
! its index into the heap for further processing.
|
|
|
|
!
|
|
|
|
!
|
|
|
|
row(j) = row(j) - rwk * uaspk(jj)
|
|
|
|
row(j) = row(j) - rwk * uval(jj)
|
|
|
|
if (abs(row(j)) < thres*nrmi) then
|
|
|
|
if (abs(row(j)) < thres*nrmi) then
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! Drop the entry.
|
|
|
|
! Drop the entry.
|
|
|
@ -771,8 +779,8 @@ contains
|
|
|
|
! Note: internal subroutine of mld_zilut_fact
|
|
|
|
! Note: internal subroutine of mld_zilut_fact
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! This routine copies a matrix row, computed by ilut_fact by applying an
|
|
|
|
! This routine copies a matrix row, computed by ilut_fact by applying an
|
|
|
|
! elimination step of the ILU(k,t) factorization, into the arrays laspk,
|
|
|
|
! elimination step of the ILU(k,t) factorization, into the arrays lval,
|
|
|
|
! uaspk, d, corresponding to the L factor, the U factor and the diagonal
|
|
|
|
! uval, d, corresponding to the L factor, the U factor and the diagonal
|
|
|
|
! of U, respectively.
|
|
|
|
! of U, respectively.
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! Note that
|
|
|
|
! Note that
|
|
|
@ -781,11 +789,11 @@ contains
|
|
|
|
! the 'lower part' of the row, and the nup+k ones in the 'upper part';
|
|
|
|
! the 'lower part' of the row, and the nup+k ones in the 'upper part';
|
|
|
|
! - the entry in the upper part of the row which has maximum absolute value
|
|
|
|
! - the entry in the upper part of the row which has maximum absolute value
|
|
|
|
! in the original matrix is included in the above nup+k entries anyway;
|
|
|
|
! in the original matrix is included in the above nup+k entries anyway;
|
|
|
|
! - the part of the row stored into uaspk is scaled by the corresponding
|
|
|
|
! - the part of the row stored into uval is scaled by the corresponding
|
|
|
|
! diagonal entry, according to the LDU form of the incomplete factorization;
|
|
|
|
! diagonal entry, according to the LDU form of the incomplete factorization;
|
|
|
|
! - the inverse of the diagonal entries of U is actually stored into d; this
|
|
|
|
! - the inverse of the diagonal entries of U is actually stored into d; this
|
|
|
|
! is then managed in the solve stage associated to the ILU(k,t) factorization;
|
|
|
|
! is then managed in the solve stage associated to the ILU(k,t) factorization;
|
|
|
|
! - the row entries are stored in laspk and uaspk according to the CSR format;
|
|
|
|
! - the row entries are stored in lval and uval according to the CSR format;
|
|
|
|
! - the array row is re-initialized for future use in mld_ilut_fact(see also
|
|
|
|
! - the array row is re-initialized for future use in mld_ilut_fact(see also
|
|
|
|
! ilut_copyin and ilut_fact).
|
|
|
|
! ilut_copyin and ilut_fact).
|
|
|
|
!
|
|
|
|
!
|
|
|
@ -825,49 +833,49 @@ contains
|
|
|
|
! examined during the elimination step carried out by the routine
|
|
|
|
! examined during the elimination step carried out by the routine
|
|
|
|
! ilut_fact.
|
|
|
|
! ilut_fact.
|
|
|
|
! l1 - integer, input/output.
|
|
|
|
! l1 - integer, input/output.
|
|
|
|
! Pointer to the last occupied entry of laspk.
|
|
|
|
! Pointer to the last occupied entry of lval.
|
|
|
|
! l2 - integer, input/output.
|
|
|
|
! l2 - integer, input/output.
|
|
|
|
! Pointer to the last occupied entry of uaspk.
|
|
|
|
! Pointer to the last occupied entry of uval.
|
|
|
|
! lia1 - integer, dimension(:), input/output.
|
|
|
|
! lja - integer, dimension(:), input/output.
|
|
|
|
! The column indices of the nonzero entries of the L factor,
|
|
|
|
! The column indices of the nonzero entries of the L factor,
|
|
|
|
! copied in laspk row by row (see mld_zilut_factint), according
|
|
|
|
! copied in lval row by row (see mld_zilut_factint), according
|
|
|
|
! to the CSR storage format.
|
|
|
|
! to the CSR storage format.
|
|
|
|
! lia2 - integer, dimension(:), input/output.
|
|
|
|
! lirp - integer, dimension(:), input/output.
|
|
|
|
! The indices identifying the first nonzero entry of each row
|
|
|
|
! The indices identifying the first nonzero entry of each row
|
|
|
|
! of the L factor, copied in laspk row by row (see
|
|
|
|
! of the L factor, copied in lval row by row (see
|
|
|
|
! mld_zilut_factint), according to the CSR storage format.
|
|
|
|
! mld_zilut_factint), according to the CSR storage format.
|
|
|
|
! laspk - complex(psb_dpk_), dimension(:), input/output.
|
|
|
|
! lval - complex(psb_dpk_), dimension(:), input/output.
|
|
|
|
! The array where the entries of the row corresponding to the
|
|
|
|
! The array where the entries of the row corresponding to the
|
|
|
|
! L factor are copied.
|
|
|
|
! L factor are copied.
|
|
|
|
! d - complex(psb_dpk_), dimension(:), input/output.
|
|
|
|
! d - complex(psb_dpk_), dimension(:), input/output.
|
|
|
|
! The array where the inverse of the diagonal entry of the
|
|
|
|
! The array where the inverse of the diagonal entry of the
|
|
|
|
! row is copied (only d(i) is used by the routine).
|
|
|
|
! row is copied (only d(i) is used by the routine).
|
|
|
|
! uia1 - integer, dimension(:), input/output.
|
|
|
|
! uja - integer, dimension(:), input/output.
|
|
|
|
! The column indices of the nonzero entries of the U factor
|
|
|
|
! The column indices of the nonzero entries of the U factor
|
|
|
|
! copied in uaspk row by row (see mld_zilut_factint), according
|
|
|
|
! copied in uval row by row (see mld_zilut_factint), according
|
|
|
|
! to the CSR storage format.
|
|
|
|
! to the CSR storage format.
|
|
|
|
! uia2 - integer, dimension(:), input/output.
|
|
|
|
! uirp - integer, dimension(:), input/output.
|
|
|
|
! The indices identifying the first nonzero entry of each row
|
|
|
|
! The indices identifying the first nonzero entry of each row
|
|
|
|
! of the U factor copied in uaspk row by row (see
|
|
|
|
! of the U factor copied in uval row by row (see
|
|
|
|
! mld_zilu_fctint), according to the CSR storage format.
|
|
|
|
! mld_zilu_fctint), according to the CSR storage format.
|
|
|
|
! uaspk - complex(psb_dpk_), dimension(:), input/output.
|
|
|
|
! uval - complex(psb_dpk_), dimension(:), input/output.
|
|
|
|
! The array where the entries of the row corresponding to the
|
|
|
|
! The array where the entries of the row corresponding to the
|
|
|
|
! U factor are copied.
|
|
|
|
! U factor are copied.
|
|
|
|
!
|
|
|
|
!
|
|
|
|
subroutine ilut_copyout(fill_in,thres,i,m,nlw,nup,jmaxup,nrmi,row, &
|
|
|
|
subroutine ilut_copyout(fill_in,thres,i,m,nlw,nup,jmaxup,nrmi,row, &
|
|
|
|
& nidx,idxs,l1,l2,lia1,lia2,laspk,d,uia1,uia2,uaspk,info)
|
|
|
|
& nidx,idxs,l1,l2,lja,lirp,lval,d,uja,uirp,uval,info)
|
|
|
|
|
|
|
|
|
|
|
|
use psb_sparse_mod
|
|
|
|
use psb_sparse_mod
|
|
|
|
|
|
|
|
|
|
|
|
implicit none
|
|
|
|
implicit none
|
|
|
|
|
|
|
|
|
|
|
|
! Arguments
|
|
|
|
! Arguments
|
|
|
|
integer, intent(in) :: fill_in,i,m,nidx,nlw,nup,jmaxup
|
|
|
|
integer, intent(in) :: fill_in,i,m,nidx,nlw,nup,jmaxup
|
|
|
|
integer, intent(in) :: idxs(:)
|
|
|
|
integer, intent(in) :: idxs(:)
|
|
|
|
integer, intent(inout) :: l1,l2, info
|
|
|
|
integer, intent(inout) :: l1,l2, info
|
|
|
|
integer, allocatable, intent(inout) :: uia1(:),uia2(:), lia1(:),lia2(:)
|
|
|
|
integer, allocatable, intent(inout) :: uja(:),uirp(:), lja(:),lirp(:)
|
|
|
|
real(psb_dpk_), intent(in) :: thres,nrmi
|
|
|
|
real(psb_dpk_), intent(in) :: thres,nrmi
|
|
|
|
complex(psb_dpk_),allocatable, intent(inout) :: uaspk(:), laspk(:)
|
|
|
|
complex(psb_dpk_),allocatable, intent(inout) :: uval(:), lval(:)
|
|
|
|
complex(psb_dpk_), intent(inout) :: row(:), d(:)
|
|
|
|
complex(psb_dpk_), intent(inout) :: row(:), d(:)
|
|
|
|
|
|
|
|
|
|
|
|
! Local variables
|
|
|
|
! Local variables
|
|
|
@ -966,21 +974,21 @@ contains
|
|
|
|
!
|
|
|
|
!
|
|
|
|
do k=1,nz
|
|
|
|
do k=1,nz
|
|
|
|
l1 = l1 + 1
|
|
|
|
l1 = l1 + 1
|
|
|
|
if (size(laspk) < l1) then
|
|
|
|
if (size(lval) < l1) then
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! Figure out a good reallocation size!
|
|
|
|
! Figure out a good reallocation size!
|
|
|
|
!
|
|
|
|
!
|
|
|
|
isz = (max((l1/i)*m,int(1.2*l1),l1+100))
|
|
|
|
isz = (max((l1/i)*m,int(1.2*l1),l1+100))
|
|
|
|
call psb_realloc(isz,laspk,info)
|
|
|
|
call psb_realloc(isz,lval,info)
|
|
|
|
if (info == psb_success_) call psb_realloc(isz,lia1,info)
|
|
|
|
if (info == psb_success_) call psb_realloc(isz,lja,info)
|
|
|
|
if (info /= psb_success_) then
|
|
|
|
if (info /= psb_success_) then
|
|
|
|
info=psb_err_from_subroutine_
|
|
|
|
info=psb_err_from_subroutine_
|
|
|
|
call psb_errpush(info,name,a_err='Allocate')
|
|
|
|
call psb_errpush(info,name,a_err='Allocate')
|
|
|
|
goto 9999
|
|
|
|
goto 9999
|
|
|
|
end if
|
|
|
|
end if
|
|
|
|
end if
|
|
|
|
end if
|
|
|
|
lia1(l1) = xwid(k)
|
|
|
|
lja(l1) = xwid(k)
|
|
|
|
laspk(l1) = xw(indx(k))
|
|
|
|
lval(l1) = xw(indx(k))
|
|
|
|
end do
|
|
|
|
end do
|
|
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
!
|
|
|
@ -1022,7 +1030,7 @@ contains
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! Compute 1/pivot
|
|
|
|
! Compute 1/pivot
|
|
|
|
!
|
|
|
|
!
|
|
|
|
d(i) = done/d(i)
|
|
|
|
d(i) = zone/d(i)
|
|
|
|
end if
|
|
|
|
end if
|
|
|
|
end if
|
|
|
|
end if
|
|
|
|
end if
|
|
|
|
end if
|
|
|
@ -1112,21 +1120,21 @@ contains
|
|
|
|
!
|
|
|
|
!
|
|
|
|
do k=1,nz
|
|
|
|
do k=1,nz
|
|
|
|
l2 = l2 + 1
|
|
|
|
l2 = l2 + 1
|
|
|
|
if (size(uaspk) < l2) then
|
|
|
|
if (size(uval) < l2) then
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! Figure out a good reallocation size!
|
|
|
|
! Figure out a good reallocation size!
|
|
|
|
!
|
|
|
|
!
|
|
|
|
isz = max((l2/i)*m,int(1.2*l2),l2+100)
|
|
|
|
isz = max((l2/i)*m,int(1.2*l2),l2+100)
|
|
|
|
call psb_realloc(isz,uaspk,info)
|
|
|
|
call psb_realloc(isz,uval,info)
|
|
|
|
if (info == psb_success_) call psb_realloc(isz,uia1,info)
|
|
|
|
if (info == psb_success_) call psb_realloc(isz,uja,info)
|
|
|
|
if (info /= psb_success_) then
|
|
|
|
if (info /= psb_success_) then
|
|
|
|
info=psb_err_from_subroutine_
|
|
|
|
info=psb_err_from_subroutine_
|
|
|
|
call psb_errpush(info,name,a_err='Allocate')
|
|
|
|
call psb_errpush(info,name,a_err='Allocate')
|
|
|
|
goto 9999
|
|
|
|
goto 9999
|
|
|
|
end if
|
|
|
|
end if
|
|
|
|
end if
|
|
|
|
end if
|
|
|
|
uia1(l2) = xwid(k)
|
|
|
|
uja(l2) = xwid(k)
|
|
|
|
uaspk(l2) = d(i)*xw(indx(k))
|
|
|
|
uval(l2) = d(i)*xw(indx(k))
|
|
|
|
end do
|
|
|
|
end do
|
|
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
!
|
|
|
@ -1138,10 +1146,10 @@ contains
|
|
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! Store the pointers to the first non occupied entry of in
|
|
|
|
! Store the pointers to the first non occupied entry of in
|
|
|
|
! laspk and uaspk
|
|
|
|
! lval and uval
|
|
|
|
!
|
|
|
|
!
|
|
|
|
lia2(i+1) = l1 + 1
|
|
|
|
lirp(i+1) = l1 + 1
|
|
|
|
uia2(i+1) = l2 + 1
|
|
|
|
uirp(i+1) = l2 + 1
|
|
|
|
|
|
|
|
|
|
|
|
call psb_erractionrestore(err_act)
|
|
|
|
call psb_erractionrestore(err_act)
|
|
|
|
return
|
|
|
|
return
|
|
|
|