diff --git a/mlprec/mld_dbjac_bld.f90 b/mlprec/mld_dbjac_bld.f90 index 85b78f31..9399b473 100644 --- a/mlprec/mld_dbjac_bld.f90 +++ b/mlprec/mld_dbjac_bld.f90 @@ -120,7 +120,6 @@ subroutine mld_dbjac_bld(a,p,upd,info,data) integer :: err_act, n_row, nrow_a,n_col, data_ integer :: ictxt,np,me character(len=20) :: name, ch_err - character(len=5), parameter :: coofmt='COO', csrfmt='CSR' if(psb_get_errstatus().ne.0) return info=0 diff --git a/mlprec/mld_dilu_fct.f90 b/mlprec/mld_dilu_fct.f90 index 839b0901..5df3f678 100644 --- a/mlprec/mld_dilu_fct.f90 +++ b/mlprec/mld_dilu_fct.f90 @@ -279,11 +279,12 @@ contains implicit none ! Arguments - integer, intent(in) :: ialg - type(psb_dspmat_type) :: a,b - integer :: m,ma,mb,l1,l2,info - integer, dimension(:) :: lia1,lia2,uia1,uia2 - real(kind(1.d0)), dimension(:) :: laspk,uaspk,d + integer, intent(in) :: ialg + type(psb_dspmat_type),intent(in) :: a,b + integer,intent(inout) :: m,l1,l2,info + integer, intent(in) :: ma,mb + integer, dimension(:), intent(inout) :: lia1,lia2,uia1,uia2 + real(kind(1.d0)), dimension(:),intent(inout) :: laspk,uaspk,d ! Local variables integer :: i,j,k,l,low1,low2,kk,jj,ll, ktrw,err_act @@ -547,15 +548,17 @@ contains implicit none ! Arguments - type(psb_dspmat_type) :: a,trw - integer :: i,m,ktrw,jd,jmin,jmax,l1,l2 - integer :: lia1(:), uia1(:) - real(kind(1.d0)) :: laspk(:), uaspk(:), dia + type(psb_dspmat_type), intent(in) :: a + type(psb_dspmat_type), intent(inout) :: trw + integer, intent(in) :: i,m,jd,jmin,jmax + integer, intent(inout) :: ktrw,l1,l2 + integer, intent(inout) :: lia1(:), uia1(:) + real(kind(1.d0)), intent(inout) :: laspk(:), uaspk(:), dia ! Local variables integer :: k,j,info,irb integer, parameter :: nrb=16 - character(len=20), parameter :: name='mld_dilu_fctint' + character(len=20), parameter :: name='ilu_copyin' character(len=20) :: ch_err if (psb_get_errstatus() /= 0) return diff --git a/mlprec/mld_diluk_fct.f90 b/mlprec/mld_diluk_fct.f90 index 7cbdd753..5cbfe632 100644 --- a/mlprec/mld_diluk_fct.f90 +++ b/mlprec/mld_diluk_fct.f90 @@ -142,7 +142,7 @@ subroutine mld_diluk_fct(fill_in,ialg,a,l,u,d,info,blck) ! ! Compute the ILU(k) or the MILU(k) factorization, depending on ialg ! - call mld_diluk_fctint(fill_in,ialg,m,a%m,a,blck_%m,blck_,& + call mld_diluk_fctint(fill_in,ialg,m,a,blck_,& & d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info) if (info /= 0) then info=4010 @@ -221,8 +221,6 @@ contains ! m - integer, output. ! The total number of rows of the local matrix to be factorized, ! i.e. ma+mb. - ! ma - integer, input. - ! The number of rows of the local submatrix stored into a. ! a - type(psb_dspmat_type), input. ! The sparse matrix structure containing the local matrix to be ! factorized. Note that, if the 'base' Additive Schwarz preconditioner @@ -230,8 +228,6 @@ contains ! (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 ! 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. ! The sparse matrix structure containing the remote rows of the ! matrix to be factorized, that have been retrieved by mld_asmat_bld @@ -265,7 +261,7 @@ contains ! info - integer, output. ! Error code. ! - subroutine mld_diluk_fctint(fill_in,ialg,m,ma,a,mb,b,& + subroutine mld_diluk_fctint(fill_in,ialg,m,a,b,& & d,laspk,lia1,lia2,uaspk,uia1,uia2,l1,l2,info) use psb_base_mod @@ -273,15 +269,15 @@ contains implicit none ! Arguments - integer, intent(in) :: fill_in, ialg - type(psb_dspmat_type) :: a,b - integer :: m,ma,mb,l1,l2,info - integer, dimension(:), allocatable :: lia1,lia2,uia1,uia2 - real(kind(1.d0)), dimension(:), allocatable :: laspk,uaspk - real(kind(1.d0)), dimension(:) :: d + integer, intent(in) :: fill_in, ialg + type(psb_dspmat_type), intent(in) :: a,b + integer, intent(inout) :: m,l1,l2,info + integer, allocatable, intent(inout) :: lia1(:),lia2(:),uia1(:),uia2(:) + real(kind(1.d0)), allocatable, intent(inout) :: laspk(:),uaspk(:) + real(kind(1.d0)), intent(inout) :: d(:) ! Local variables - integer :: i, ktrw,err_act,nidx + integer :: ma,mb,i, ktrw,err_act,nidx integer, allocatable :: uplevs(:), rowlevs(:),idxs(:) real(kind(1.d0)), allocatable :: row(:) type(psb_int_heap) :: heap @@ -293,7 +289,23 @@ contains info=0 call psb_erractionsave(err_act) - m = ma+mb + select case(ialg) + case(mld_ilu_n_,mld_milu_n_) + ! Ok + case default + info=35 + call psb_errpush(info,name,i_err=(/2,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 + + ma = a%m + mb = b%m + m = ma+mb ! ! Allocate a temporary buffer for the iluk_copyin function @@ -346,27 +358,31 @@ contains ! ! Copy into trw the i-th local row of the matrix, stored in a ! - call iluk_copyin(i,ma,a,1,m,row,rowlevs,heap,ktrw,trw) + call iluk_copyin(i,ma,a,1,m,row,rowlevs,heap,ktrw,trw,info) 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,1,m,row,rowlevs,heap,ktrw,trw) + call iluk_copyin(i-ma,mb,b,1,m,row,rowlevs,heap,ktrw,trw,info) endif - + ! Do an elimination step on the current row. It turns out we only ! need to keep track of fill levels for the upper triangle, hence we ! do not have a lowlevs variable. ! - call iluk_fact(fill_in,i,m,row,rowlevs,heap,& - & d,uia1,uia2,uaspk,uplevs,nidx,idxs) + if (info == 0) call iluk_fact(fill_in,i,m,row,rowlevs,heap,& + & d,uia1,uia2,uaspk,uplevs,nidx,idxs,info) ! ! Copy the row into laspk/d(i)/uaspk ! - call iluk_copyout(fill_in,ialg,i,m,row,rowlevs,nidx,idxs,& - & l1,l2,lia1,lia2,laspk,d,uia1,uia2,uaspk,uplevs) - + if (info == 0) call iluk_copyout(fill_in,ialg,i,m,row,rowlevs,nidx,idxs,& + & l1,l2,lia1,lia2,laspk,d,uia1,uia2,uaspk,uplevs,info) + if (info /= 0) then + info=4001 + call psb_errpush(info,name,a_err='Copy/factor loop') + goto 9999 + end if end do ! @@ -462,22 +478,25 @@ 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) + subroutine iluk_copyin(i,m,a,jmin,jmax,row,rowlevs,heap,ktrw,trw,info) use psb_base_mod implicit none ! Arguments - type(psb_dspmat_type) :: a,trw - integer :: i, rowlevs(:),m,ktrw,jmin,jmax - real(kind(1.d0)) :: row(:) - type(psb_int_heap) :: heap + type(psb_dspmat_type), intent(in) :: a + type(psb_dspmat_type), intent(inout) :: trw + integer, intent(in) :: i,m,jmin,jmax + integer, intent(inout) :: ktrw,info + integer, intent(inout) :: rowlevs(:) + real(kind(1.d0)), intent(inout) :: row(:) + type(psb_int_heap), intent(inout) :: heap ! Local variables - integer :: k,j,info,irb + integer :: k,j,irb,err_act integer, parameter :: nrb=16 - character(len=20), parameter :: name='mld_diluk_fctint' + character(len=20), parameter :: name='iluk_copyin' character(len=20) :: ch_err if (psb_get_errstatus() /= 0) return @@ -619,25 +638,29 @@ 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,m,row,rowlevs,heap,d,uia1,uia2,uaspk,uplevs,nidx,idxs) + subroutine iluk_fact(fill_in,i,m,row,rowlevs,heap,d,uia1,uia2,uaspk,uplevs,nidx,idxs,info) use psb_base_mod implicit none ! Arguments - type(psb_int_heap) :: heap - integer :: i,m, rowlevs(:),fill_in,nidx - integer, allocatable :: idxs(:) - integer :: uia1(:),uia2(:),uplevs(:) - real(kind(1.d0)) :: row(:), uaspk(:),d(:) + type(psb_int_heap), intent(inout) :: heap + integer, intent(in) :: i,m, fill_in + integer, intent(inout) :: nidx,info + integer, intent(inout) :: rowlevs(:) + integer, allocatable, intent(inout) :: idxs(:) + integer, intent(inout) :: uia1(:),uia2(:),uplevs(:) + real(kind(1.d0)), intent(inout) :: row(:), uaspk(:),d(:) ! Local variables - integer :: k,j,lrwk,jj,info, lastk + integer :: k,j,lrwk,jj,lastk, iret real(kind(1.d0)) :: rwk + info = 0 if (.not.allocated(idxs)) then allocate(idxs(200),stat=info) + if (info /= 0) return endif nidx = 0 lastk = -1 @@ -646,9 +669,9 @@ contains ! Do while there are indices to be processed ! do - - call psb_heap_get_first(k,heap,info) - if (info < 0) exit + ! Beware: (iret < 0) means that the heap is empty, not an error. + call psb_heap_get_first(k,heap,iret) + if (iret < 0) return ! ! Just in case an index has been put on the heap more than once. @@ -659,6 +682,7 @@ contains nidx = nidx + 1 if (nidx>size(idxs)) then call psb_realloc(nidx+psb_heap_resize,idxs,info) + if (info /= 0) return end if idxs(nidx) = k @@ -674,7 +698,8 @@ contains do jj=uia2(k),uia2(k+1)-1 j = uia1(jj) if (j<=k) then - write(0,*) 'Error in accessing upper mat???',j,k,jj + info = -i + return endif ! ! Insert the index into the heap for further processing. @@ -685,6 +710,7 @@ contains ! if (rowlevs(j)<0) then call psb_insert_heap(j,heap,info) + if (info /= 0) return rowlevs(j) = abs(rowlevs(j)) end if ! @@ -788,25 +814,27 @@ contains ! U factor above the current row. ! subroutine iluk_copyout(fill_in,ialg,i,m,row,rowlevs,nidx,idxs,& - & l1,l2,lia1,lia2,laspk,d,uia1,uia2,uaspk,uplevs) + & l1,l2,lia1,lia2,laspk,d,uia1,uia2,uaspk,uplevs,info) use psb_base_mod implicit none ! Arguments - integer :: fill_in, ialg, i,rowlevs(:), l1, l2, m, nidx, idxs(:) - integer, allocatable :: uia1(:), uia2(:), lia1(:), lia2(:),uplevs(:) - real(kind(1.d0)),allocatable :: uaspk(:), laspk(:) - real(kind(1.d0)) :: row(:), d(:) + integer, intent(in) :: fill_in, ialg, i, m, nidx + integer, intent(inout) :: l1, l2, info + integer, intent(inout) :: rowlevs(:), idxs(:) + integer, allocatable, intent(inout) :: uia1(:), uia2(:), lia1(:), lia2(:),uplevs(:) + real(kind(1.d0)), allocatable, intent(inout) :: uaspk(:), laspk(:) + real(kind(1.d0)), intent(inout) :: row(:), d(:) ! Local variables - integer :: j,isz,info,err_act,int_err(5),idxp + integer :: j,isz,err_act,int_err(5),idxp character(len=20), parameter :: name='mld_diluk_fctint' character(len=20) :: ch_err if (psb_get_errstatus() /= 0) return - info=0 + info = 0 call psb_erractionsave(err_act) d(i) = dzero diff --git a/mlprec/mld_zbjac_bld.f90 b/mlprec/mld_zbjac_bld.f90 index 93e15e2c..68de775d 100644 --- a/mlprec/mld_zbjac_bld.f90 +++ b/mlprec/mld_zbjac_bld.f90 @@ -121,7 +121,6 @@ subroutine mld_zbjac_bld(a,p,upd,info,data) integer :: err_act, n_row, nrow_a,n_col, data_ integer :: ictxt,np,me character(len=20) :: name, ch_err - character(len=5), parameter :: coofmt='COO', csrfmt='CSR' if(psb_get_errstatus().ne.0) return info=0 diff --git a/mlprec/mld_zilu_fct.f90 b/mlprec/mld_zilu_fct.f90 index 52a3a8e4..3a78ba37 100644 --- a/mlprec/mld_zilu_fct.f90 +++ b/mlprec/mld_zilu_fct.f90 @@ -278,11 +278,12 @@ contains implicit none ! Arguments - integer, intent(in) :: ialg - type(psb_zspmat_type) :: a,b - integer :: m,ma,mb,l1,l2,info - integer, dimension(:) :: lia1,lia2,uia1,uia2 - complex(kind(1.d0)), dimension(:) :: laspk,uaspk,d + integer, intent(in) :: ialg + type(psb_zspmat_type),intent(in) :: a,b + integer,intent(inout) :: m,l1,l2,info + integer, intent(in) :: ma,mb + integer, dimension(:), intent(inout) :: lia1,lia2,uia1,uia2 + complex(kind(1.d0)), dimension(:), intent(inout) :: laspk,uaspk,d ! Local variables integer :: i,j,k,l,low1,low2,kk,jj,ll, ktrw,err_act @@ -546,15 +547,17 @@ contains implicit none ! Arguments - type(psb_zspmat_type) :: a,trw - integer :: i,m,ktrw,jd,jmin,jmax,l1,l2 - integer :: lia1(:),uia1(:) - complex(kind(1.d0)) :: laspk(:), uaspk(:), dia + type(psb_zspmat_type), intent(in) :: a + type(psb_zspmat_type), intent(inout) :: trw + integer, intent(in) :: i,m,jd,jmin,jmax + integer, intent(inout) :: ktrw,l1,l2 + integer, intent(inout) :: lia1(:), uia1(:) + complex(kind(1.d0)), intent(inout) :: laspk(:), uaspk(:), dia ! Local variables integer :: k,j,info,irb integer, parameter :: nrb=16 - character(len=20), parameter :: name='mld_dilu_fctint' + character(len=20), parameter :: name='ilu_copyin' character(len=20) :: ch_err if (psb_get_errstatus() /= 0) return diff --git a/mlprec/mld_ziluk_fct.f90 b/mlprec/mld_ziluk_fct.f90 index ddda262d..71ab7f34 100644 --- a/mlprec/mld_ziluk_fct.f90 +++ b/mlprec/mld_ziluk_fct.f90 @@ -141,7 +141,7 @@ subroutine mld_ziluk_fct(fill_in,ialg,a,l,u,d,info,blck) ! ! Compute the ILU(k) or the MILU(k) factorization, depending on ialg ! - call mld_ziluk_fctint(fill_in,ialg,m,a%m,a,blck_%m,blck_,& + call mld_ziluk_fctint(fill_in,ialg,m,a,blck_,& & d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info) if (info /= 0) then info=4010 @@ -220,8 +220,6 @@ contains ! m - integer, output. ! The total number of rows of the local matrix to be factorized, ! i.e. ma+mb. - ! ma - integer, input. - ! The number of rows of the local submatrix stored into a. ! a - type(psb_zspmat_type), input. ! The sparse matrix structure containing the local matrix to be ! factorized. Note that, if the 'base' Additive Schwarz preconditioner @@ -229,8 +227,6 @@ contains ! (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 ! 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. ! The sparse matrix structure containing the remote rows of the ! matrix to be factorized, that have been retrieved by mld_asmat_bld @@ -264,7 +260,7 @@ contains ! info - integer, output. ! Error code. ! - subroutine mld_ziluk_fctint(fill_in,ialg,m,ma,a,mb,b,& + subroutine mld_ziluk_fctint(fill_in,ialg,m,a,b,& & d,laspk,lia1,lia2,uaspk,uia1,uia2,l1,l2,info) use psb_base_mod @@ -272,15 +268,15 @@ contains implicit none ! Arguments - integer, intent(in) :: fill_in, ialg - type(psb_zspmat_type) :: a,b - integer :: m,ma,mb,l1,l2,info - integer, dimension(:), allocatable :: lia1,lia2,uia1,uia2 - complex(kind(1.d0)), dimension(:), allocatable :: laspk,uaspk - complex(kind(1.d0)), dimension(:) :: d + integer, intent(in) :: fill_in, ialg + type(psb_zspmat_type), intent(in) :: a,b + integer, intent(inout) :: m,l1,l2,info + integer, allocatable, intent(inout) :: lia1(:),lia2(:),uia1(:),uia2(:) + complex(kind(1.d0)), allocatable, intent(inout) :: laspk(:),uaspk(:) + complex(kind(1.d0)), intent(inout) :: d(:) ! Local variables - integer :: i, ktrw,err_act, nidx + integer :: ma,mb,i, ktrw,err_act,nidx integer, allocatable :: uplevs(:), rowlevs(:),idxs(:) complex(kind(1.d0)), allocatable :: row(:) type(psb_int_heap) :: heap @@ -293,7 +289,23 @@ contains info=0 call psb_erractionsave(err_act) - m = ma+mb + select case(ialg) + case(mld_ilu_n_,mld_milu_n_) + ! Ok + case default + info=35 + call psb_errpush(info,name,i_err=(/2,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 + + ma = a%m + mb = b%m + m = ma+mb ! ! Allocate a temporary buffer for the iluk_copyin function @@ -346,27 +358,31 @@ contains ! ! Copy into trw the i-th local row of the matrix, stored in a ! - call iluk_copyin(i,ma,a,1,m,row,rowlevs,heap,ktrw,trw) + call iluk_copyin(i,ma,a,1,m,row,rowlevs,heap,ktrw,trw,info) 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,1,m,row,rowlevs,heap,ktrw,trw) + call iluk_copyin(i-ma,mb,b,1,m,row,rowlevs,heap,ktrw,trw,info) endif ! Do an elimination step on the current row. It turns out we only ! need to keep track of fill levels for the upper triangle, hence we ! do not have a lowlevs variable. ! - call iluk_fact(fill_in,i,m,row,rowlevs,heap,& - & d,uia1,uia2,uaspk,uplevs,nidx,idxs) + if (info == 0) call iluk_fact(fill_in,i,m,row,rowlevs,heap,& + & d,uia1,uia2,uaspk,uplevs,nidx,idxs,info) ! ! Copy the row into laspk/d(i)/uaspk ! - call iluk_copyout(fill_in,ialg,i,m,row,rowlevs,nidx,idxs,& - & l1,l2,lia1,lia2,laspk,d,uia1,uia2,uaspk,uplevs) - + if (info == 0) call iluk_copyout(fill_in,ialg,i,m,row,rowlevs,nidx,idxs,& + & l1,l2,lia1,lia2,laspk,d,uia1,uia2,uaspk,uplevs,info) + if (info /= 0) then + info=4001 + call psb_errpush(info,name,a_err='Copy/factor loop') + goto 9999 + end if end do ! @@ -462,22 +478,25 @@ 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) + subroutine iluk_copyin(i,m,a,jmin,jmax,row,rowlevs,heap,ktrw,trw,info) use psb_base_mod implicit none ! Arguments - type(psb_zspmat_type) :: a,trw - integer :: i, rowlevs(:),m,ktrw,jmin,jmax - complex(kind(1.d0)) :: row(:) - type(psb_int_heap) :: heap + type(psb_zspmat_type), intent(in) :: a + type(psb_zspmat_type), intent(inout) :: trw + integer, intent(in) :: i,m,jmin,jmax + integer, intent(inout) :: ktrw,info + integer, intent(inout) :: rowlevs(:) + complex(kind(1.d0)), intent(inout) :: row(:) + type(psb_int_heap), intent(inout) :: heap ! Local variables - integer :: k,j,info,irb + integer :: k,j,irb,err_act integer, parameter :: nrb=16 - character(len=20), parameter :: name='mld_ziluk_fctint' + character(len=20), parameter :: name='iluk_copyin' character(len=20) :: ch_err if (psb_get_errstatus() /= 0) return @@ -619,25 +638,29 @@ 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,m,row,rowlevs,heap,d,uia1,uia2,uaspk,uplevs,nidx,idxs) + subroutine iluk_fact(fill_in,i,m,row,rowlevs,heap,d,uia1,uia2,uaspk,uplevs,nidx,idxs,info) use psb_base_mod implicit none ! Arguments - type(psb_int_heap) :: heap - integer :: i,m, rowlevs(:),fill_in,nidx - integer, allocatable :: idxs(:) - integer :: uia1(:),uia2(:),uplevs(:) - complex(kind(1.d0)) :: row(:), uaspk(:),d(:) + type(psb_int_heap), intent(inout) :: heap + integer, intent(in) :: i,m, fill_in + integer, intent(inout) :: nidx,info + integer, intent(inout) :: rowlevs(:) + integer, allocatable, intent(inout) :: idxs(:) + integer, intent(inout) :: uia1(:),uia2(:),uplevs(:) + complex(kind(1.d0)), intent(inout) :: row(:), uaspk(:),d(:) ! Local variables - integer :: k,j,lrwk,jj,info, lastk + integer :: k,j,lrwk,jj,lastk, iret complex(kind(1.d0)) :: rwk + info = 0 if (.not.allocated(idxs)) then allocate(idxs(200),stat=info) + if (info /= 0) return endif nidx = 0 lastk = -1 @@ -646,9 +669,9 @@ contains ! Do while there are indices to be processed ! do - - call psb_heap_get_first(k,heap,info) - if (info < 0) exit + ! Beware: (iret < 0) means that the heap is empty, not an error. + call psb_heap_get_first(k,heap,iret) + if (iret < 0) return ! ! Just in case an index has been put on the heap more than once. @@ -659,6 +682,7 @@ contains nidx = nidx + 1 if (nidx>size(idxs)) then call psb_realloc(nidx+psb_heap_resize,idxs,info) + if (info /= 0) return end if idxs(nidx) = k if ((row(k) /= zzero).and.(rowlevs(k) <= fill_in).and.(k