fdurastante
/
amg4psblas
mirror of https://github.com/sfilippone/amg4psblas.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Fixed interfaces of internal routines for ILUT.

Browse Source
stopcriterion
Salvatore Filippone 17 years ago
parent 581cd9b6ac
commit 55d8db4c62
2 changed files with 242 additions and 118 deletions
Show Stats Download Patch File Download Diff File

166
mlprec/mld_dilut_fct.f90
Unescape Escape View File

@ -116,6 +116,20 @@ subroutine mld_dilut_fct(fill_in,thres,ialg,a,l,u,d,info,blck)
info = 0 info = 0
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
select case(ialg)
case(mld_ilu_n_,mld_milu_n_)
! Ok
case default
info=35
call psb_errpush(info,name,i_err=(/3,ialg,0,0,0/))
goto 9999
end select
if (fill_in < 0) then
info=35
call psb_errpush(info,name,i_err=(/1,fill_in,0,0,0/))
goto 9999
end if
! !
! Point to / allocate memory for the incomplete factorization ! Point to / allocate memory for the incomplete factorization
! !
@ -141,7 +155,7 @@ subroutine mld_dilut_fct(fill_in,thres,ialg,a,l,u,d,info,blck)
! !
! Compute the ILU(k,t) factorization ! Compute the ILU(k,t) factorization
! !
call mld_dilut_fctint(fill_in,thres,ialg,m,a%m,a,blck_%m,blck_,& call mld_dilut_fctint(fill_in,thres,m,a,blck_,&
& d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info) & d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info)
if (info /= 0) then if (info /= 0) then
info=4010 info=4010
@ -218,8 +232,6 @@ contains
! m - integer, output. ! m - integer, output.
! The total number of rows of the local matrix to be factorized, ! The total number of rows of the local matrix to be factorized,
! i.e. ma+mb. ! i.e. ma+mb.
! ma - integer, input.
! The number of rows of the local submatrix stored into a.
! a - type(psb_dspmat_type), input. ! a - type(psb_dspmat_type), input.
! The sparse matrix structure containing the local matrix to be ! The sparse matrix structure containing the local matrix to be
! factorized. Note that, if the 'base' Additive Schwarz preconditioner ! factorized. Note that, if the 'base' Additive Schwarz preconditioner
@ -227,8 +239,6 @@ contains
! (see mld_bjac_bld), then a contains only the 'original' local part ! (see mld_bjac_bld), then a contains only the 'original' local part
! of the matrix to be factorized, i.e. the rows of the matrix held ! of the matrix to be factorized, i.e. the rows of the matrix held
! by the calling process according to the initial data distribution. ! by the calling process according to the initial data distribution.
! mb - integer, input.
! The number of rows of the local submatrix stored into b.
! b - type(psb_dspmat_type), input. ! b - type(psb_dspmat_type), input.
! The sparse matrix structure containing the remote rows of the ! The sparse matrix structure containing the remote rows of the
! matrix to be factorized, that have been retrieved by mld_asmat_bld ! matrix to be factorized, that have been retrieved by mld_asmat_bld
@ -262,7 +272,7 @@ contains
! info - integer, output. ! info - integer, output.
! Error code. ! Error code.
! !
subroutine mld_dilut_fctint(fill_in,thres,ialg,m,ma,a,mb,b,& subroutine mld_dilut_fctint(fill_in,thres,m,a,b,&
& d,laspk,lia1,lia2,uaspk,uia1,uia2,l1,l2,info) & d,laspk,lia1,lia2,uaspk,uia1,uia2,l1,l2,info)
use psb_base_mod use psb_base_mod
@ -270,16 +280,16 @@ contains
implicit none implicit none
! Arguments ! Arguments
integer, intent(in) :: fill_in, ialg integer, intent(in) :: fill_in
real(kind(1.d0)), intent(in) :: thres real(kind(1.d0)), intent(in) :: thres
type(psb_dspmat_type) :: a,b type(psb_dspmat_type), intent(in) :: a,b
integer :: m,ma,mb,l1,l2,info integer, intent(inout) :: m,l1,l2,info
integer, dimension(:), allocatable :: lia1,lia2,uia1,uia2 integer, allocatable, intent(inout) :: lia1(:),lia2(:),uia1(:),uia2(:)
real(kind(1.d0)), dimension(:), allocatable :: laspk,uaspk real(kind(1.d0)), allocatable, intent(inout) :: laspk(:),uaspk(:)
real(kind(1.d0)), dimension(:) :: d real(kind(1.d0)), intent(inout) :: d(:)
! Local Variables ! Local Variables
integer :: i, ktrw,err_act,nidx,nlw,nup,jmaxup integer :: i, ktrw,err_act,nidx,nlw,nup,jmaxup, ma, mb
real(kind(1.d0)) :: nrmi real(kind(1.d0)) :: nrmi
integer, allocatable :: idxs(:) integer, allocatable :: idxs(:)
real(kind(1.d0)), allocatable :: row(:) real(kind(1.d0)), allocatable :: row(:)
@ -292,6 +302,9 @@ contains
info = 0 info = 0
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
ma = a%m
mb = b%m
m = ma+mb m = ma+mb
! !
@ -339,21 +352,27 @@ contains
! !
d(i) = dzero d(i) = dzero
if (i<=ma) then if (i<=ma) then
call ilut_copyin(i,ma,a,i,1,m,nlw,nup,jmaxup,nrmi,row,heap,ktrw,trw) call ilut_copyin(i,ma,a,i,1,m,nlw,nup,jmaxup,nrmi,row,heap,ktrw,trw,info)
else else
call ilut_copyin(i-ma,mb,b,i,1,m,nlw,nup,jmaxup,nrmi,row,heap,ktrw,trw) call ilut_copyin(i-ma,mb,b,i,1,m,nlw,nup,jmaxup,nrmi,row,heap,ktrw,trw,info)
endif endif
! !
! Do an elimination step on current row ! Do an elimination step on current row
! !
call ilut_fact(fill_in,thres,i,m,nrmi,row,heap,& if (info == 0) call ilut_fact(fill_in,thres,i,m,nrmi,row,heap,&
& d,uia1,uia2,uaspk,nidx,idxs) & d,uia1,uia2,uaspk,nidx,idxs,info)
! !
! Copy the row into laspk/d(i)/uaspk ! Copy the row into laspk/d(i)/uaspk
! !
call ilut_copyout(fill_in,thres,i,m,nlw,nup,jmaxup,nrmi,row,nidx,idxs,& if (info == 0) call ilut_copyout(fill_in,thres,i,m,nlw,nup,jmaxup,nrmi,row,nidx,idxs,&
& l1,l2,lia1,lia2,laspk,d,uia1,uia2,uaspk) & l1,l2,lia1,lia2,laspk,d,uia1,uia2,uaspk,info)
if (info /= 0) then
info=4001
call psb_errpush(info,name,a_err='Copy/factor loop')
goto 9999
end if
end do end do
@ -467,15 +486,17 @@ 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 ilut_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,nrmi,row,heap,ktrw,trw) subroutine ilut_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,nrmi,row,heap,ktrw,trw,info)
use psb_base_mod use psb_base_mod
implicit none implicit none
type(psb_dspmat_type) :: a,trw type(psb_dspmat_type), intent(in) :: a
integer :: i, m,ktrw,jmin,jmax,jd,nlw,nup,jmaxup type(psb_dspmat_type), intent(inout) :: trw
real(kind(1.d0)) :: nrmi,row(:) integer, intent(in) :: i, m,jmin,jmax,jd
type(psb_int_heap) :: heap integer, intent(inout) :: ktrw,nlw,nup,jmaxup,info
real(kind(1.d0)), intent(inout) :: nrmi,row(:)
type(psb_int_heap), intent(inout) :: heap
integer :: k,j,info,irb,kin,nz integer :: k,j,irb,kin,nz
integer, parameter :: nrb=16 integer, parameter :: nrb=16
real(kind(1.d0)) :: dmaxup real(kind(1.d0)) :: dmaxup
real(kind(1.d0)), external :: dnrm2 real(kind(1.d0)), external :: dnrm2
@ -487,6 +508,11 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
call psb_init_heap(heap,info) call psb_init_heap(heap,info)
if (info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_init_heap')
goto 9999
end if
! !
! nrmi is the norm of the current sparse row (for the time being, ! nrmi is the norm of the current sparse row (for the time being,
@ -512,6 +538,11 @@ contains
if ((jmin<=k).and.(k<=jmax)) then if ((jmin<=k).and.(k<=jmax)) then
row(k) = a%aspk(j) row(k) = a%aspk(j)
call psb_insert_heap(k,heap,info) call psb_insert_heap(k,heap,info)
if (info /= 0) then
info=4010
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
@ -539,8 +570,7 @@ contains
call psb_sp_getblk(i,a,trw,info,lrw=i+irb-1) call psb_sp_getblk(i,a,trw,info,lrw=i+irb-1)
if (info /= 0) then if (info /= 0) then
info=4010 info=4010
ch_err='psb_sp_getblk' call psb_errpush(info,name,a_err='psb_sp_getblk')
call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
ktrw=1 ktrw=1
@ -554,6 +584,11 @@ contains
if ((jmin<=k).and.(k<=jmax)) then if ((jmin<=k).and.(k<=jmax)) then
row(k) = trw%aspk(ktrw) row(k) = trw%aspk(ktrw)
call psb_insert_heap(k,heap,info) call psb_insert_heap(k,heap,info)
if (info /= 0) then
info=4010
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
@ -646,26 +681,28 @@ contains
! to retain its allocation, done by this routine. ! to retain its allocation, done by this routine.
! !
subroutine ilut_fact(fill_in,thres,i,m,nrmi,row,heap,& subroutine ilut_fact(fill_in,thres,i,m,nrmi,row,heap,&
& d,uia1,uia2,uaspk,nidx,idxs) & d,uia1,uia2,uaspk,nidx,idxs,info)
use psb_base_mod use psb_base_mod
implicit none implicit none
! Arguments ! Arguments
type(psb_int_heap) :: heap type(psb_int_heap), intent(inout) :: heap
integer :: i,m,fill_in,nidx integer, intent(in) :: i,m,fill_in
integer, intent(inout) :: nidx,info
real(kind(1.d0)), intent(in) :: thres,nrmi real(kind(1.d0)), intent(in) :: thres,nrmi
integer, allocatable :: idxs(:) integer, allocatable, intent(inout) :: idxs(:)
integer :: uia1(:),uia2(:) integer, intent(inout) :: uia1(:),uia2(:)
real(kind(1.d0)) :: row(:), uaspk(:),d(:) real(kind(1.d0)), intent(inout) :: row(:), uaspk(:),d(:)
! Local Variables ! Local Variables
integer :: k,j,jj,info, lastk integer :: k,j,jj,lastk,iret
real(kind(1.d0)) :: rwk real(kind(1.d0)) :: rwk
info = 0
call psb_ensure_size(200,idxs,info) call psb_ensure_size(200,idxs,info)
if (info /= 0) return
nidx = 0 nidx = 0
lastk = -1 lastk = -1
! !
@ -673,8 +710,8 @@ contains
! !
do do
call psb_heap_get_first(k,heap,info) call psb_heap_get_first(k,heap,iret)
if (info < 0) exit if (iret < 0) exit
! !
! An index may have been put on the heap more than once. ! An index may have been put on the heap more than once.
@ -704,7 +741,8 @@ contains
do jj=uia2(k),uia2(k+1)-1 do jj=uia2(k),uia2(k+1)-1
j = uia1(jj) j = uia1(jj)
if (j<=k) then if (j<=k) then
write(0,*) 'Error in accessing upper mat???',j,k,jj info = -i
return
endif endif
! !
! Update row(j) and, if it is not to be discarded, insert ! Update row(j) and, if it is not to be discarded, insert
@ -721,6 +759,7 @@ contains
! Do the insertion. ! Do the insertion.
! !
call psb_insert_heap(j,heap,info) call psb_insert_heap(j,heap,info)
if (info /= 0) return
endif endif
end do end do
end if end if
@ -731,6 +770,7 @@ contains
! !
nidx = nidx + 1 nidx = nidx + 1
call psb_ensure_size(nidx,idxs,info,addsz=psb_heap_resize) call psb_ensure_size(nidx,idxs,info,addsz=psb_heap_resize)
if (info /= 0) return
idxs(nidx) = k idxs(nidx) = k
end do end do
@ -827,28 +867,29 @@ contains
! 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) & nidx,idxs,l1,l2,lia1,lia2,laspk,d,uia1,uia2,uaspk,info)
use psb_base_mod use psb_base_mod
implicit none implicit none
! Arguments ! Arguments
integer :: fill_in,i,l1,l2,m,nidx,idxs(:) integer, intent(in) :: fill_in,i,m,nidx,nlw,nup,jmaxup
integer :: nlw,nup,jmaxup integer, intent(in) :: idxs(:)
integer, allocatable :: uia1(:),uia2(:), lia1(:),lia2(:) integer, intent(inout) :: l1,l2, info
integer, allocatable, intent(inout) :: uia1(:),uia2(:), lia1(:),lia2(:)
real(kind(1.d0)), intent(in) :: thres,nrmi real(kind(1.d0)), intent(in) :: thres,nrmi
real(kind(1.d0)),allocatable :: uaspk(:), laspk(:) real(kind(1.d0)),allocatable, intent(inout) :: uaspk(:), laspk(:)
real(kind(1.d0)) :: row(:), d(:) real(kind(1.d0)), intent(inout) :: row(:), d(:)
! Local variables ! Local variables
real(kind(1.d0)),allocatable :: xw(:) real(kind(1.d0)),allocatable :: xw(:)
integer, allocatable :: xwid(:), indx(:) integer, allocatable :: xwid(:), indx(:)
real(kind(1.d0)) :: witem real(kind(1.d0)) :: witem
integer :: widx integer :: widx
integer :: k,isz,info,err_act,int_err(5),idxp, nz integer :: k,isz,err_act,int_err(5),idxp, nz
type(psb_double_idx_heap) :: heap type(psb_double_idx_heap) :: heap
character(len=20), parameter :: name='mld_dilut_fctint' character(len=20), parameter :: name='ilut_copyout'
character(len=20) :: ch_err character(len=20) :: ch_err
logical :: fndmaxup logical :: fndmaxup
@ -898,7 +939,11 @@ contains
xw(nz) = witem xw(nz) = witem
xwid(nz) = widx xwid(nz) = widx
call psb_insert_heap(witem,widx,heap,info) call psb_insert_heap(witem,widx,heap,info)
if (info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_insert_heap')
goto 9999
end if
end do end do
! !
@ -912,6 +957,12 @@ contains
nz = nlw+fill_in nz = nlw+fill_in
do k=1,nz do k=1,nz
call psb_heap_get_first(witem,widx,heap,info) call psb_heap_get_first(witem,widx,heap,info)
if (info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_heap_get_first')
goto 9999
end if
xw(k) = witem xw(k) = witem
xwid(k) = widx xwid(k) = widx
end do end do
@ -957,12 +1008,12 @@ contains
end if end if
end if end if
if (idxp > size(idxs)) then if (idxp > size(idxs)) then
write(0,*) 'Warning: missing diagonal element in the row ' !!$ write(0,*) 'Warning: missing diagonal element in the row '
else else
if (idxs(idxp) > i) then if (idxs(idxp) > i) then
write(0,*) 'Warning: missing diagonal element in the row ' !!$ write(0,*) 'Warning: missing diagonal element in the row '
else if (idxs(idxp) /= i) then else if (idxs(idxp) /= i) then
write(0,*) 'Warning: impossible error: diagonal has vanished' !!$ write(0,*) 'Warning: impossible error: diagonal has vanished'
else else
! !
! Copy the diagonal entry ! Copy the diagonal entry
@ -993,6 +1044,12 @@ contains
! !
call psb_init_heap(heap,info,dir=psb_asort_down_) call psb_init_heap(heap,info,dir=psb_asort_down_)
if (info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_init_heap')
goto 9999
end if
nz = 0 nz = 0
do do
@ -1000,11 +1057,11 @@ contains
if (idxp > nidx) exit if (idxp > nidx) exit
widx = idxs(idxp) widx = idxs(idxp)
if (widx <= i) then if (widx <= i) then
write(0,*) 'Warning: lower triangle in upper copy',widx,i,idxp,idxs(idxp) !!$ write(0,*) 'Warning: lower triangle in upper copy',widx,i,idxp,idxs(idxp)
cycle cycle
end if end if
if (widx > m) then if (widx > m) then
write(0,*) 'Warning: impossible value',widx,i,idxp,idxs(idxp) !!$ write(0,*) 'Warning: impossible value',widx,i,idxp,idxs(idxp)
cycle cycle
end if end if
witem = row(widx) witem = row(widx)
@ -1019,6 +1076,11 @@ contains
xw(nz) = witem xw(nz) = witem
xwid(nz) = widx xwid(nz) = widx
call psb_insert_heap(witem,widx,heap,info) call psb_insert_heap(witem,widx,heap,info)
if (info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_insert_heap')
goto 9999
end if
end do end do

170
mlprec/mld_zilut_fct.f90
Unescape Escape View File

@ -115,6 +115,20 @@ subroutine mld_zilut_fct(fill_in,thres,ialg,a,l,u,d,info,blck)
info = 0 info = 0
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
select case(ialg)
case(mld_ilu_n_,mld_milu_n_)
! Ok
case default
info=35
call psb_errpush(info,name,i_err=(/3,ialg,0,0,0/))
goto 9999
end select
if (fill_in < 0) then
info=35
call psb_errpush(info,name,i_err=(/1,fill_in,0,0,0/))
goto 9999
end if
! !
! Point to / allocate memory for the incomplete factorization ! Point to / allocate memory for the incomplete factorization
! !
@ -140,7 +154,7 @@ subroutine mld_zilut_fct(fill_in,thres,ialg,a,l,u,d,info,blck)
! !
! Compute the ILU(k,t) factorization ! Compute the ILU(k,t) factorization
! !
call mld_zilut_fctint(fill_in,thres,ialg,m,a%m,a,blck_%m,blck_,& call mld_zilut_fctint(fill_in,thres,m,a,blck_,&
& d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info) & d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info)
if (info /= 0) then if (info /= 0) then
info=4010 info=4010
@ -217,8 +231,6 @@ contains
! m - integer, output. ! m - integer, output.
! The total number of rows of the local matrix to be factorized, ! The total number of rows of the local matrix to be factorized,
! i.e. ma+mb. ! i.e. ma+mb.
! ma - integer, input.
! The number of rows of the local submatrix stored into a.
! a - type(psb_zspmat_type), input. ! a - type(psb_zspmat_type), input.
! The sparse matrix structure containing the local matrix to be ! The sparse matrix structure containing the local matrix to be
! factorized. Note that, if the 'base' Additive Schwarz preconditioner ! factorized. Note that, if the 'base' Additive Schwarz preconditioner
@ -226,8 +238,6 @@ contains
! (see mld_bjac_bld), then a contains only the 'original' local part ! (see mld_bjac_bld), then a contains only the 'original' local part
! of the matrix to be factorized, i.e. the rows of the matrix held ! of the matrix to be factorized, i.e. the rows of the matrix held
! by the calling process according to the initial data distribution. ! by the calling process according to the initial data distribution.
! mb - integer, input.
! The number of rows of the local submatrix stored into b.
! b - type(psb_zspmat_type), input. ! b - type(psb_zspmat_type), input.
! The sparse matrix structure containing the remote rows of the ! The sparse matrix structure containing the remote rows of the
! matrix to be factorized, that have been retrieved by mld_asmat_bld ! matrix to be factorized, that have been retrieved by mld_asmat_bld
@ -261,7 +271,7 @@ contains
! info - integer, output. ! info - integer, output.
! Error code. ! Error code.
! !
subroutine mld_zilut_fctint(fill_in,thres,ialg,m,ma,a,mb,b,& subroutine mld_zilut_fctint(fill_in,thres,m,a,b,&
& d,laspk,lia1,lia2,uaspk,uia1,uia2,l1,l2,info) & d,laspk,lia1,lia2,uaspk,uia1,uia2,l1,l2,info)
use psb_base_mod use psb_base_mod
@ -269,16 +279,16 @@ contains
implicit none implicit none
! Arguments ! Arguments
integer, intent(in) :: fill_in, ialg integer, intent(in) :: fill_in
real(kind(1.d0)), intent(in) :: thres real(kind(1.d0)), intent(in) :: thres
type(psb_zspmat_type) :: a,b type(psb_zspmat_type), intent(in) :: a,b
integer :: m,ma,mb,l1,l2,info integer, intent(inout) :: m,l1,l2,info
integer, dimension(:), allocatable :: lia1,lia2,uia1,uia2 integer, allocatable, intent(inout) :: lia1(:),lia2(:),uia1(:),uia2(:)
complex(kind(1.d0)), dimension(:), allocatable :: laspk,uaspk complex(kind(1.d0)), allocatable, intent(inout) :: laspk(:),uaspk(:)
complex(kind(1.d0)), dimension(:) :: d complex(kind(1.d0)), intent(inout) :: d(:)
! Local Variables ! Local Variables
integer :: i, ktrw,err_act,nidx,nlw,nup,jmaxup integer :: i, ktrw,err_act,nidx,nlw,nup,jmaxup, ma, mb
real(kind(1.d0)) :: nrmi real(kind(1.d0)) :: nrmi
integer, allocatable :: idxs(:) integer, allocatable :: idxs(:)
complex(kind(1.d0)), allocatable :: row(:) complex(kind(1.d0)), allocatable :: row(:)
@ -291,6 +301,9 @@ contains
info = 0 info = 0
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
ma = a%m
mb = b%m
m = ma+mb m = ma+mb
! !
@ -338,21 +351,27 @@ contains
! !
d(i) = zzero d(i) = zzero
if (i<=ma) then if (i<=ma) then
call ilut_copyin(i,ma,a,i,1,m,nlw,nup,jmaxup,nrmi,row,heap,ktrw,trw) call ilut_copyin(i,ma,a,i,1,m,nlw,nup,jmaxup,nrmi,row,heap,ktrw,trw,info)
else else
call ilut_copyin(i-ma,mb,b,i,1,m,nlw,nup,jmaxup,nrmi,row,heap,ktrw,trw) call ilut_copyin(i-ma,mb,b,i,1,m,nlw,nup,jmaxup,nrmi,row,heap,ktrw,trw,info)
endif endif
! !
! Do an elimination step on current row ! Do an elimination step on current row
! !
call ilut_fact(fill_in,thres,i,m,nrmi,row,heap,& if (info == 0) call ilut_fact(fill_in,thres,i,m,nrmi,row,heap,&
& d,uia1,uia2,uaspk,nidx,idxs) & d,uia1,uia2,uaspk,nidx,idxs,info)
! !
! Copy the row into laspk/d(i)/uaspk ! Copy the row into laspk/d(i)/uaspk
! !
call ilut_copyout(fill_in,thres,i,m,nlw,nup,jmaxup,nrmi,row,nidx,idxs,& if (info == 0) call ilut_copyout(fill_in,thres,i,m,nlw,nup,jmaxup,nrmi,row,nidx,idxs,&
& l1,l2,lia1,lia2,laspk,d,uia1,uia2,uaspk) & l1,l2,lia1,lia2,laspk,d,uia1,uia2,uaspk,info)
if (info /= 0) then
info=4001
call psb_errpush(info,name,a_err='Copy/factor loop')
goto 9999
end if
end do end do
@ -466,16 +485,18 @@ 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 ilut_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,nrmi,row,heap,ktrw,trw) subroutine ilut_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,nrmi,row,heap,ktrw,trw,info)
use psb_base_mod use psb_base_mod
implicit none implicit none
type(psb_zspmat_type) :: a,trw type(psb_zspmat_type), intent(in) :: a
integer :: i, m,ktrw,jmin,jmax,jd,nlw,nup,jmaxup type(psb_zspmat_type), intent(inout) :: trw
real(kind(1.d0)) :: nrmi integer, intent(in) :: i, m,jmin,jmax,jd
complex(kind(1.d0)) :: row(:) integer, intent(inout) :: ktrw,nlw,nup,jmaxup,info
type(psb_int_heap) :: heap real(kind(1.d0)), intent(inout) :: nrmi
complex(kind(1.d0)), intent(inout) :: row(:)
integer :: k,j,info,irb,kin,nz type(psb_int_heap), intent(inout) :: heap
integer :: k,j,irb,kin,nz
integer, parameter :: nrb=16 integer, parameter :: nrb=16
real(kind(1.d0)) :: dmaxup real(kind(1.d0)) :: dmaxup
real(kind(1.d0)), external :: dznrm2 real(kind(1.d0)), external :: dznrm2
@ -487,6 +508,11 @@ contains
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
call psb_init_heap(heap,info) call psb_init_heap(heap,info)
if (info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_init_heap')
goto 9999
end if
! !
! nrmi is the norm of the current sparse row (for the time being, ! nrmi is the norm of the current sparse row (for the time being,
@ -512,6 +538,11 @@ contains
if ((jmin<=k).and.(k<=jmax)) then if ((jmin<=k).and.(k<=jmax)) then
row(k) = a%aspk(j) row(k) = a%aspk(j)
call psb_insert_heap(k,heap,info) call psb_insert_heap(k,heap,info)
if (info /= 0) then
info=4010
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
@ -539,8 +570,7 @@ contains
call psb_sp_getblk(i,a,trw,info,lrw=i+irb-1) call psb_sp_getblk(i,a,trw,info,lrw=i+irb-1)
if (info /= 0) then if (info /= 0) then
info=4010 info=4010
ch_err='psb_sp_getblk' call psb_errpush(info,name,a_err='psb_sp_getblk')
call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
ktrw=1 ktrw=1
@ -554,6 +584,11 @@ contains
if ((jmin<=k).and.(k<=jmax)) then if ((jmin<=k).and.(k<=jmax)) then
row(k) = trw%aspk(ktrw) row(k) = trw%aspk(ktrw)
call psb_insert_heap(k,heap,info) call psb_insert_heap(k,heap,info)
if (info /= 0) then
info=4010
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
@ -646,26 +681,28 @@ contains
! to retain its allocation, done by this routine. ! to retain its allocation, done by this routine.
! !
subroutine ilut_fact(fill_in,thres,i,m,nrmi,row,heap,& subroutine ilut_fact(fill_in,thres,i,m,nrmi,row,heap,&
& d,uia1,uia2,uaspk,nidx,idxs) & d,uia1,uia2,uaspk,nidx,idxs,info)
use psb_base_mod use psb_base_mod
implicit none implicit none
! Arguments ! Arguments
type(psb_int_heap) :: heap type(psb_int_heap), intent(inout) :: heap
integer :: i,m,fill_in,nidx integer, intent(in) :: i,m,fill_in
integer, intent(inout) :: nidx,info
real(kind(1.d0)), intent(in) :: thres,nrmi real(kind(1.d0)), intent(in) :: thres,nrmi
integer, allocatable :: idxs(:) integer, allocatable, intent(inout) :: idxs(:)
integer :: uia1(:),uia2(:) integer, intent(inout) :: uia1(:),uia2(:)
complex(kind(1.d0)) :: row(:), uaspk(:),d(:) complex(kind(1.d0)), intent(inout) :: row(:), uaspk(:),d(:)
! Local Variables ! Local Variables
integer :: k,j,jj,info, lastk integer :: k,j,jj,lastk, iret
complex(kind(1.d0)) :: rwk complex(kind(1.d0)) :: rwk
info = 0
call psb_ensure_size(200,idxs,info) call psb_ensure_size(200,idxs,info)
if (info /= 0) return
nidx = 0 nidx = 0
lastk = -1 lastk = -1
! !
@ -673,8 +710,8 @@ contains
! !
do do
call psb_heap_get_first(k,heap,info) call psb_heap_get_first(k,heap,iret)
if (info < 0) exit if (iret < 0) exit
! !
! An index may have been put on the heap more than once. ! An index may have been put on the heap more than once.
@ -704,7 +741,8 @@ contains
do jj=uia2(k),uia2(k+1)-1 do jj=uia2(k),uia2(k+1)-1
j = uia1(jj) j = uia1(jj)
if (j<=k) then if (j<=k) then
write(0,*) 'Error in accessing upper mat???',j,k,jj info = -i
return
endif endif
! !
! Update row(j) and, if it is not to be discarded, insert ! Update row(j) and, if it is not to be discarded, insert
@ -721,6 +759,7 @@ contains
! Do the insertion. ! Do the insertion.
! !
call psb_insert_heap(j,heap,info) call psb_insert_heap(j,heap,info)
if (info /= 0) return
endif endif
end do end do
end if end if
@ -731,6 +770,7 @@ contains
! !
nidx = nidx + 1 nidx = nidx + 1
call psb_ensure_size(nidx,idxs,info,addsz=psb_heap_resize) call psb_ensure_size(nidx,idxs,info,addsz=psb_heap_resize)
if (info /= 0) return
idxs(nidx) = k idxs(nidx) = k
end do end do
@ -827,28 +867,29 @@ contains
! 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) & nidx,idxs,l1,l2,lia1,lia2,laspk,d,uia1,uia2,uaspk,info)
use psb_base_mod use psb_base_mod
implicit none implicit none
! Arguments ! Arguments
integer :: fill_in,i,l1,l2,m,nidx,idxs(:) integer, intent(in) :: fill_in,i,m,nidx,nlw,nup,jmaxup
integer :: nlw,nup,jmaxup integer, intent(in) :: idxs(:)
integer, allocatable :: uia1(:),uia2(:), lia1(:),lia2(:) integer, intent(inout) :: l1,l2, info
integer, allocatable, intent(inout) :: uia1(:),uia2(:), lia1(:),lia2(:)
real(kind(1.d0)), intent(in) :: thres,nrmi real(kind(1.d0)), intent(in) :: thres,nrmi
complex(kind(1.d0)),allocatable :: uaspk(:), laspk(:) complex(kind(1.d0)),allocatable, intent(inout) :: uaspk(:), laspk(:)
complex(kind(1.d0)) :: row(:), d(:) complex(kind(1.d0)), intent(inout) :: row(:), d(:)
! Local variables ! Local variables
complex(kind(1.d0)),allocatable :: xw(:) complex(kind(1.d0)),allocatable :: xw(:)
integer, allocatable :: xwid(:), indx(:) integer, allocatable :: xwid(:), indx(:)
complex(kind(1.d0)) :: witem complex(kind(1.d0)) :: witem
integer :: widx integer :: widx
integer :: k,isz,info,err_act,int_err(5),idxp, nz integer :: k,isz,err_act,int_err(5),idxp, nz
type(psb_dcomplex_idx_heap) :: heap type(psb_dcomplex_idx_heap) :: heap
character(len=20), parameter :: name='mld_zilut_fctint' character(len=20), parameter :: name='ilut_copyout'
character(len=20) :: ch_err character(len=20) :: ch_err
logical :: fndmaxup logical :: fndmaxup
@ -898,7 +939,11 @@ contains
xw(nz) = witem xw(nz) = witem
xwid(nz) = widx xwid(nz) = widx
call psb_insert_heap(witem,widx,heap,info) call psb_insert_heap(witem,widx,heap,info)
if (info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_insert_heap')
goto 9999
end if
end do end do
! !
@ -912,6 +957,12 @@ contains
nz = nlw+fill_in nz = nlw+fill_in
do k=1,nz do k=1,nz
call psb_heap_get_first(witem,widx,heap,info) call psb_heap_get_first(witem,widx,heap,info)
if (info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_heap_get_first')
goto 9999
end if
xw(k) = witem xw(k) = witem
xwid(k) = widx xwid(k) = widx
end do end do
@ -957,12 +1008,12 @@ contains
end if end if
end if end if
if (idxp > size(idxs)) then if (idxp > size(idxs)) then
write(0,*) 'Warning: missing diagonal element in the row ' !!$ write(0,*) 'Warning: missing diagonal element in the row '
else else
if (idxs(idxp) > i) then if (idxs(idxp) > i) then
write(0,*) 'Warning: missing diagonal element in the row ' !!$ write(0,*) 'Warning: missing diagonal element in the row '
else if (idxs(idxp) /= i) then else if (idxs(idxp) /= i) then
write(0,*) 'Warning: impossible error: diagonal has vanished' !!$ write(0,*) 'Warning: impossible error: diagonal has vanished'
else else
! !
! Copy the diagonal entry ! Copy the diagonal entry
@ -993,6 +1044,12 @@ contains
! !
call psb_init_heap(heap,info,dir=psb_asort_down_) call psb_init_heap(heap,info,dir=psb_asort_down_)
if (info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_init_heap')
goto 9999
end if
nz = 0 nz = 0
do do
@ -1000,11 +1057,11 @@ contains
if (idxp > nidx) exit if (idxp > nidx) exit
widx = idxs(idxp) widx = idxs(idxp)
if (widx <= i) then if (widx <= i) then
write(0,*) 'Warning: lower triangle in upper copy',widx,i,idxp,idxs(idxp) !!$ write(0,*) 'Warning: lower triangle in upper copy',widx,i,idxp,idxs(idxp)
cycle cycle
end if end if
if (widx > m) then if (widx > m) then
write(0,*) 'Warning: impossible value',widx,i,idxp,idxs(idxp) !!$ write(0,*) 'Warning: impossible value',widx,i,idxp,idxs(idxp)
cycle cycle
end if end if
witem = row(widx) witem = row(widx)
@ -1019,6 +1076,11 @@ contains
xw(nz) = witem xw(nz) = witem
xwid(nz) = widx xwid(nz) = widx
call psb_insert_heap(witem,widx,heap,info) call psb_insert_heap(witem,widx,heap,info)
if (info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_insert_heap')
goto 9999
end if
end do end do

Write Preview
Loading…
Cancel
Save