diff --git a/mlprec/mld_dbjac_aply.f90 b/mlprec/mld_dbjac_aply.f90 index 4179d455..683e2612 100644 --- a/mlprec/mld_dbjac_aply.f90 +++ b/mlprec/mld_dbjac_aply.f90 @@ -114,7 +114,7 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) if (prec%iprcparm(smooth_sweeps_) == 1) then select case(prec%iprcparm(sub_solve_)) - case(ilu_n_,ilu_t_) + case(ilu_n_,milu_n_,ilu_t_) select case(toupper(trans)) case('N') @@ -204,7 +204,7 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) tx = dzero ty = dzero select case(prec%iprcparm(sub_solve_)) - case(ilu_n_,ilu_t_) + case(ilu_n_,milu_n_,ilu_t_) do i=1, prec%iprcparm(smooth_sweeps_) ! X(k+1) = M^-1*(b-N*X(k)) ty(1:n_row) = x(1:n_row) diff --git a/mlprec/mld_dbjac_bld.f90 b/mlprec/mld_dbjac_bld.f90 index e588e7ea..768ad176 100644 --- a/mlprec/mld_dbjac_bld.f90 +++ b/mlprec/mld_dbjac_bld.f90 @@ -182,7 +182,7 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info) select case(p%iprcparm(sub_solve_)) - case(ilu_n_,ilu_t_) + case(ilu_n_,milu_n_,ilu_t_) call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_) if (info /= 0) then @@ -268,7 +268,7 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info) select case(p%iprcparm(sub_solve_)) - case(ilu_n_,ilu_t_) + case(ilu_n_,milu_n_,ilu_t_) if (p%iprcparm(smooth_sweeps_) > 1) then diff --git a/mlprec/mld_dilu_bld.f90 b/mlprec/mld_dilu_bld.f90 index 37afddbe..d0889ad1 100644 --- a/mlprec/mld_dilu_bld.f90 +++ b/mlprec/mld_dilu_bld.f90 @@ -149,7 +149,8 @@ subroutine mld_dilu_bld(a,desc_a,p,upd,info,blck) ! Ok, factor the matrix. ! t5 = psb_wtime() - call mld_ilu_fct(a,p%av(l_pr_),p%av(u_pr_),p%d,info,blck=blck) + call mld_ilu_fct(p%iprcparm(sub_solve_),a,p%av(l_pr_),p%av(u_pr_),& + & p%d,info,blck=blck) if(info/=0) then info=4010 ch_err='mld_ilu_fct' diff --git a/mlprec/mld_dilu_fct.f90 b/mlprec/mld_dilu_fct.f90 index 303f0b2d..1cf29c2c 100644 --- a/mlprec/mld_dilu_fct.f90 +++ b/mlprec/mld_dilu_fct.f90 @@ -34,7 +34,7 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -subroutine mld_dilu_fct(a,l,u,d,info,blck) +subroutine mld_dilu_fct(ialg,a,l,u,d,info,blck) ! ! This routine copies and factors "on the fly" from A and BLCK @@ -42,9 +42,11 @@ subroutine mld_dilu_fct(a,l,u,d,info,blck) ! ! use psb_base_mod + use mld_prec_mod, mld_protect_name => mld_dilu_fct implicit none ! .. Scalar Arguments .. - integer, intent(out) :: info + integer, intent(in) :: ialg + integer, intent(out) :: info ! .. Array Arguments .. type(psb_dspmat_type),intent(in) :: a type(psb_dspmat_type),intent(inout) :: l,u @@ -85,8 +87,7 @@ subroutine mld_dilu_fct(a,l,u,d,info,blck) blck_%m=0 endif -!!$ write(0,*) 'ilu_fct: ',size(l%ia2),size(u%ia2),a%m,blck_%m - call mld_dilu_fctint(m,a%m,a,blck_%m,blck_,& + call mld_dilu_fctint(ialg,m,a%m,a,blck_%m,blck_,& & d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info) if(info.ne.0) then info=4010 @@ -130,10 +131,11 @@ subroutine mld_dilu_fct(a,l,u,d,info,blck) return contains - subroutine mld_dilu_fctint(m,ma,a,mb,b,& + subroutine mld_dilu_fctint(ialg,m,ma,a,mb,b,& & d,laspk,lia1,lia2,uaspk,uia1,uia2,l1,l2,info) implicit none + integer, intent(in) :: ialg type(psb_dspmat_type) :: a,b integer :: m,ma,mb,l1,l2,info integer, dimension(:) :: lia1,lia2,uia1,uia2 @@ -285,9 +287,13 @@ contains enddo end if ! - ! for milu al=1.; for ilu al=0. - ! al = 1.d0 - ! dia = dia - al*temp*aup(jj) + ! If we get here we missed the cycle updateloop, + ! which means that this entry does not match; thus + ! we take it out of diagonal for MILU. + ! + if (ialg == milu_n_) then + dia = dia - temp*uaspk(jj) + end if enddo updateloop enddo ! @@ -425,9 +431,13 @@ contains enddo end if ! - ! for milu al=1.; for ilu al=0. - ! al = 1.d0 - ! dia = dia - al*temp*aup(jj) + ! If we get here we missed the cycle updateloop, + ! which means that this entry does not match; thus + ! we take it out of diagonal for MILU. + ! + if (ialg == milu_n_) then + dia = dia - temp*uaspk(jj) + end if enddo updateloopb enddo ! diff --git a/mlprec/mld_prec_mod.f90 b/mlprec/mld_prec_mod.f90 index 7bf1fe31..11e0aec7 100644 --- a/mlprec/mld_prec_mod.f90 +++ b/mlprec/mld_prec_mod.f90 @@ -420,17 +420,19 @@ module mld_prec_mod interface mld_ilu_fct - subroutine mld_dilu_fct(a,l,u,d,info,blck) + subroutine mld_dilu_fct(ialg,a,l,u,d,info,blck) use psb_base_mod - integer, intent(out) :: info + integer, intent(in) :: ialg + integer, intent(out) :: info type(psb_dspmat_type),intent(in) :: a type(psb_dspmat_type),intent(inout) :: l,u type(psb_dspmat_type),intent(in), optional, target :: blck real(kind(1.d0)), intent(inout) :: d(:) end subroutine mld_dilu_fct - subroutine mld_zilu_fct(a,l,u,d,info,blck) + subroutine mld_zilu_fct(ialg,a,l,u,d,info,blck) use psb_base_mod - integer, intent(out) :: info + integer, intent(in) :: ialg + integer, intent(out) :: info type(psb_zspmat_type),intent(in) :: a type(psb_zspmat_type),intent(inout) :: l,u type(psb_zspmat_type),intent(in), optional, target :: blck diff --git a/mlprec/mld_prec_type.f90 b/mlprec/mld_prec_type.f90 index d45f1914..907425f0 100644 --- a/mlprec/mld_prec_type.f90 +++ b/mlprec/mld_prec_type.f90 @@ -167,8 +167,8 @@ module mld_prec_type integer, parameter :: pre_smooth_=1, post_smooth_=2, twoside_smooth_=3,& & max_smooth_=twoside_smooth_ ! Legal values for entry: sub_solve_ - integer, parameter :: f_none_=0,ilu_n_=1,ilu_t_=2,slu_=3 - integer, parameter :: umf_=4, sludist_=5 + integer, parameter :: f_none_=0,ilu_n_=1,milu_n_=2, ilu_t_=3 + integer, parameter :: slu_=4, umf_=5, sludist_=6 ! Legal values for entry: aggr_alg_ integer, parameter :: dec_aggr_=0, sym_dec_aggr_=1 integer, parameter :: glb_aggr_=2, new_dec_aggr_=3 @@ -184,7 +184,7 @@ module mld_prec_type ! Legal values for entry: sub_ren_ integer, parameter :: renum_none_=0, renum_glb_=1, renum_gps_=2 - ! Entries in dprcparm: ILU(E) epsilon, smoother omega + ! Entries in dprcparm: ILU(T) epsilon, smoother omega integer, parameter :: fact_eps_=1 integer, parameter :: aggr_damp_=2 integer, parameter :: aggr_thresh_=3 @@ -217,8 +217,9 @@ module mld_prec_type & ml_names(0:3)=(/'None ','Additive ','Multiplicative',& & 'New ML '/) character(len=15), parameter, private :: & - & fact_names(0:5)=(/'None ','ILU(n) ',& - & 'ILU(T) ','Sparse SuperLU','UMFPACK Sp. LU',& + & fact_names(0:6)=(/'None ','ILU(n) ',& + & 'MILU(n) ','ILU(T) ',& + & 'Sparse SuperLU','UMFPACK Sp. LU',& & 'SuperLU_Dist '/) interface mld_base_precfree diff --git a/mlprec/mld_zbjac_aply.f90 b/mlprec/mld_zbjac_aply.f90 index 09bd72de..61b3e5f0 100644 --- a/mlprec/mld_zbjac_aply.f90 +++ b/mlprec/mld_zbjac_aply.f90 @@ -112,7 +112,7 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) select case(prec%iprcparm(sub_solve_)) - case(ilu_n_,ilu_t_) + case(ilu_n_,milu_n_,ilu_t_) select case(toupper(trans)) case('N') @@ -208,7 +208,7 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) tx = zzero ty = zzero select case(prec%iprcparm(sub_solve_)) - case(ilu_n_,ilu_t_) + case(ilu_n_,milu_n_,ilu_t_) do i=1, prec%iprcparm(smooth_sweeps_) ! X(k+1) = M^-1*(b-N*X(k)) ty(1:n_row) = x(1:n_row) diff --git a/mlprec/mld_zbjac_bld.f90 b/mlprec/mld_zbjac_bld.f90 index 4ad1d370..c133267e 100644 --- a/mlprec/mld_zbjac_bld.f90 +++ b/mlprec/mld_zbjac_bld.f90 @@ -183,7 +183,7 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info) select case(p%iprcparm(sub_solve_)) - case(ilu_n_,ilu_t_) + case(ilu_n_,milu_n_,ilu_t_) call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_) if (info /= 0) then @@ -269,7 +269,7 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info) select case(p%iprcparm(sub_solve_)) - case(ilu_n_,ilu_t_) + case(ilu_n_,milu_n_,ilu_t_) if (p%iprcparm(smooth_sweeps_) > 1) then diff --git a/mlprec/mld_zilu_bld.f90 b/mlprec/mld_zilu_bld.f90 index afd45347..77240e4e 100644 --- a/mlprec/mld_zilu_bld.f90 +++ b/mlprec/mld_zilu_bld.f90 @@ -148,7 +148,8 @@ subroutine mld_zilu_bld(a,desc_a,p,upd,info,blck) ! Ok, factor the matrix. ! t5 = psb_wtime() - call mld_ilu_fct(a,p%av(l_pr_),p%av(u_pr_),p%d,info,blck=blck) + call mld_ilu_fct(p%iprcparm(sub_solve_),a,p%av(l_pr_),p%av(u_pr_),& + & p%d,info,blck=blck) if(info/=0) then info=4010 ch_err='mld_ilu_fct' diff --git a/mlprec/mld_zilu_fct.f90 b/mlprec/mld_zilu_fct.f90 index 468f284d..f2e9e7fc 100644 --- a/mlprec/mld_zilu_fct.f90 +++ b/mlprec/mld_zilu_fct.f90 @@ -34,7 +34,7 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -subroutine mld_zilu_fct(a,l,u,d,info,blck) +subroutine mld_zilu_fct(ialg,a,l,u,d,info,blck) ! ! This routine copies and factors "on the fly" from A and BLCK @@ -42,9 +42,11 @@ subroutine mld_zilu_fct(a,l,u,d,info,blck) ! ! use psb_base_mod + use mld_prec_mod, mld_protect_name => mld_zilu_fct implicit none ! .. Scalar Arguments .. - integer, intent(out) :: info + integer, intent(in) :: ialg + integer, intent(out) :: info ! .. Array Arguments .. type(psb_zspmat_type),intent(in) :: a type(psb_zspmat_type),intent(inout) :: l,u @@ -82,7 +84,7 @@ subroutine mld_zilu_fct(a,l,u,d,info,blck) blck_%m=0 endif - call mld_zilu_fctint(m,a%m,a,blck_%m,blck_,& + call mld_zilu_fctint(ialg,m,a%m,a,blck_%m,blck_,& & d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info) if(info.ne.0) then info=4010 @@ -126,10 +128,11 @@ subroutine mld_zilu_fct(a,l,u,d,info,blck) return contains - subroutine mld_zilu_fctint(m,ma,a,mb,b,& + subroutine mld_zilu_fctint(ialg,m,ma,a,mb,b,& & d,laspk,lia1,lia2,uaspk,uia1,uia2,l1,l2,info) implicit none + integer, intent(in) :: ialg type(psb_zspmat_type) :: a,b integer :: m,ma,mb,l1,l2,info integer, dimension(:) :: lia1,lia2,uia1,uia2 @@ -281,9 +284,13 @@ contains enddo end if ! - ! for milu al=1.; for ilu al=0. - ! al = 1.d0 - ! dia = dia - al*temp*aup(jj) + ! If we get here we missed the cycle updateloop, + ! which means that this entry does not match; thus + ! we take it out of diagonal for MILU. + ! + if (ialg == milu_n_) then + dia = dia - temp*uaspk(jj) + end if enddo updateloop enddo ! @@ -417,9 +424,13 @@ contains enddo end if ! - ! for milu al=1.; for ilu al=0. - ! al = 1.d0 - ! dia = dia - al*temp*aup(jj) + ! If we get here we missed the cycle updateloop, + ! which means that this entry does not match; thus + ! we take it out of diagonal for MILU. + ! + if (ialg == milu_n_) then + dia = dia - temp*uaspk(jj) + end if enddo updateloopb enddo !