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

Implemented MILU(0).

Browse Source
stopcriterion
Salvatore Filippone 17 years ago
parent 8a7f307da9
commit 6742f1e24b
10 changed files with 66 additions and 40 deletions
Show Stats Download Patch File Download Diff File

4
mlprec/mld_dbjac_aply.f90
Unescape Escape View File

@ -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 if (prec%iprcparm(smooth_sweeps_) == 1) then
select case(prec%iprcparm(sub_solve_)) select case(prec%iprcparm(sub_solve_))
case(ilu_n_,ilu_t_) case(ilu_n_,milu_n_,ilu_t_)
select case(toupper(trans)) select case(toupper(trans))
case('N') case('N')
@ -204,7 +204,7 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
tx = dzero tx = dzero
ty = dzero ty = dzero
select case(prec%iprcparm(sub_solve_)) 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_) do i=1, prec%iprcparm(smooth_sweeps_)
! X(k+1) = M^-1*(b-N*X(k)) ! X(k+1) = M^-1*(b-N*X(k))
ty(1:n_row) = x(1:n_row) ty(1:n_row) = x(1:n_row)

4
mlprec/mld_dbjac_bld.f90
Unescape Escape View File

@ -182,7 +182,7 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
select case(p%iprcparm(sub_solve_)) 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_) call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_)
if (info /= 0) then if (info /= 0) then
@ -268,7 +268,7 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
select case(p%iprcparm(sub_solve_)) 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 if (p%iprcparm(smooth_sweeps_) > 1) then

3
mlprec/mld_dilu_bld.f90
Unescape Escape View File

@ -149,7 +149,8 @@ subroutine mld_dilu_bld(a,desc_a,p,upd,info,blck)
! Ok, factor the matrix. ! Ok, factor the matrix.
! !
t5 = psb_wtime() 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 if(info/=0) then
info=4010 info=4010
ch_err='mld_ilu_fct' ch_err='mld_ilu_fct'

32
mlprec/mld_dilu_fct.f90
Unescape Escape View File

@ -34,7 +34,7 @@
!!$ POSSIBILITY OF SUCH DAMAGE. !!$ 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 ! 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 psb_base_mod
use mld_prec_mod, mld_protect_name => mld_dilu_fct
implicit none implicit none
! .. Scalar Arguments .. ! .. Scalar Arguments ..
integer, intent(out) :: info integer, intent(in) :: ialg
integer, intent(out) :: info
! .. Array Arguments .. ! .. Array Arguments ..
type(psb_dspmat_type),intent(in) :: a type(psb_dspmat_type),intent(in) :: a
type(psb_dspmat_type),intent(inout) :: l,u 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 blck_%m=0
endif endif
!!$ write(0,*) 'ilu_fct: ',size(l%ia2),size(u%ia2),a%m,blck_%m call mld_dilu_fctint(ialg,m,a%m,a,blck_%m,blck_,&
call mld_dilu_fctint(m,a%m,a,blck_%m,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.ne.0) then if(info.ne.0) then
info=4010 info=4010
@ -130,10 +131,11 @@ subroutine mld_dilu_fct(a,l,u,d,info,blck)
return return
contains 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) & d,laspk,lia1,lia2,uaspk,uia1,uia2,l1,l2,info)
implicit none implicit none
integer, intent(in) :: ialg
type(psb_dspmat_type) :: a,b type(psb_dspmat_type) :: a,b
integer :: m,ma,mb,l1,l2,info integer :: m,ma,mb,l1,l2,info
integer, dimension(:) :: lia1,lia2,uia1,uia2 integer, dimension(:) :: lia1,lia2,uia1,uia2
@ -285,9 +287,13 @@ contains
enddo enddo
end if end if
! !
! for milu al=1.; for ilu al=0. ! If we get here we missed the cycle updateloop,
! al = 1.d0 ! which means that this entry does not match; thus
! dia = dia - al*temp*aup(jj) ! we take it out of diagonal for MILU.
!
if (ialg == milu_n_) then
dia = dia - temp*uaspk(jj)
end if
enddo updateloop enddo updateloop
enddo enddo
! !
@ -425,9 +431,13 @@ contains
enddo enddo
end if end if
! !
! for milu al=1.; for ilu al=0. ! If we get here we missed the cycle updateloop,
! al = 1.d0 ! which means that this entry does not match; thus
! dia = dia - al*temp*aup(jj) ! we take it out of diagonal for MILU.
!
if (ialg == milu_n_) then
dia = dia - temp*uaspk(jj)
end if
enddo updateloopb enddo updateloopb
enddo enddo
! !

10
mlprec/mld_prec_mod.f90
Unescape Escape View File

@ -420,17 +420,19 @@ module mld_prec_mod
interface mld_ilu_fct 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 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(in) :: a
type(psb_dspmat_type),intent(inout) :: l,u type(psb_dspmat_type),intent(inout) :: l,u
type(psb_dspmat_type),intent(in), optional, target :: blck type(psb_dspmat_type),intent(in), optional, target :: blck
real(kind(1.d0)), intent(inout) :: d(:) real(kind(1.d0)), intent(inout) :: d(:)
end subroutine mld_dilu_fct 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 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(in) :: a
type(psb_zspmat_type),intent(inout) :: l,u type(psb_zspmat_type),intent(inout) :: l,u
type(psb_zspmat_type),intent(in), optional, target :: blck type(psb_zspmat_type),intent(in), optional, target :: blck

11
mlprec/mld_prec_type.f90
Unescape Escape View File

@ -167,8 +167,8 @@ module mld_prec_type
integer, parameter :: pre_smooth_=1, post_smooth_=2, twoside_smooth_=3,& integer, parameter :: pre_smooth_=1, post_smooth_=2, twoside_smooth_=3,&
& max_smooth_=twoside_smooth_ & max_smooth_=twoside_smooth_
! Legal values for entry: sub_solve_ ! Legal values for entry: sub_solve_
integer, parameter :: f_none_=0,ilu_n_=1,ilu_t_=2,slu_=3 integer, parameter :: f_none_=0,ilu_n_=1,milu_n_=2, ilu_t_=3
integer, parameter :: umf_=4, sludist_=5 integer, parameter :: slu_=4, umf_=5, sludist_=6
! Legal values for entry: aggr_alg_ ! Legal values for entry: aggr_alg_
integer, parameter :: dec_aggr_=0, sym_dec_aggr_=1 integer, parameter :: dec_aggr_=0, sym_dec_aggr_=1
integer, parameter :: glb_aggr_=2, new_dec_aggr_=3 integer, parameter :: glb_aggr_=2, new_dec_aggr_=3
@ -184,7 +184,7 @@ module mld_prec_type
! Legal values for entry: sub_ren_ ! Legal values for entry: sub_ren_
integer, parameter :: renum_none_=0, renum_glb_=1, renum_gps_=2 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 :: fact_eps_=1
integer, parameter :: aggr_damp_=2 integer, parameter :: aggr_damp_=2
integer, parameter :: aggr_thresh_=3 integer, parameter :: aggr_thresh_=3
@ -217,8 +217,9 @@ module mld_prec_type
& ml_names(0:3)=(/'None ','Additive ','Multiplicative',& & ml_names(0:3)=(/'None ','Additive ','Multiplicative',&
& 'New ML '/) & 'New ML '/)
character(len=15), parameter, private :: & character(len=15), parameter, private :: &
& fact_names(0:5)=(/'None ','ILU(n) ',& & fact_names(0:6)=(/'None ','ILU(n) ',&
& 'ILU(T) ','Sparse SuperLU','UMFPACK Sp. LU',& & 'MILU(n) ','ILU(T) ',&
& 'Sparse SuperLU','UMFPACK Sp. LU',&
& 'SuperLU_Dist '/) & 'SuperLU_Dist '/)
interface mld_base_precfree interface mld_base_precfree

4
mlprec/mld_zbjac_aply.f90
Unescape Escape View File

@ -112,7 +112,7 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
select case(prec%iprcparm(sub_solve_)) select case(prec%iprcparm(sub_solve_))
case(ilu_n_,ilu_t_) case(ilu_n_,milu_n_,ilu_t_)
select case(toupper(trans)) select case(toupper(trans))
case('N') case('N')
@ -208,7 +208,7 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
tx = zzero tx = zzero
ty = zzero ty = zzero
select case(prec%iprcparm(sub_solve_)) 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_) do i=1, prec%iprcparm(smooth_sweeps_)
! X(k+1) = M^-1*(b-N*X(k)) ! X(k+1) = M^-1*(b-N*X(k))
ty(1:n_row) = x(1:n_row) ty(1:n_row) = x(1:n_row)

4
mlprec/mld_zbjac_bld.f90
Unescape Escape View File

@ -183,7 +183,7 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
select case(p%iprcparm(sub_solve_)) 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_) call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_)
if (info /= 0) then if (info /= 0) then
@ -269,7 +269,7 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
select case(p%iprcparm(sub_solve_)) 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 if (p%iprcparm(smooth_sweeps_) > 1) then

3
mlprec/mld_zilu_bld.f90
Unescape Escape View File

@ -148,7 +148,8 @@ subroutine mld_zilu_bld(a,desc_a,p,upd,info,blck)
! Ok, factor the matrix. ! Ok, factor the matrix.
! !
t5 = psb_wtime() 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 if(info/=0) then
info=4010 info=4010
ch_err='mld_ilu_fct' ch_err='mld_ilu_fct'

31
mlprec/mld_zilu_fct.f90
Unescape Escape View File

@ -34,7 +34,7 @@
!!$ POSSIBILITY OF SUCH DAMAGE. !!$ 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 ! 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 psb_base_mod
use mld_prec_mod, mld_protect_name => mld_zilu_fct
implicit none implicit none
! .. Scalar Arguments .. ! .. Scalar Arguments ..
integer, intent(out) :: info integer, intent(in) :: ialg
integer, intent(out) :: info
! .. Array Arguments .. ! .. Array Arguments ..
type(psb_zspmat_type),intent(in) :: a type(psb_zspmat_type),intent(in) :: a
type(psb_zspmat_type),intent(inout) :: l,u 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 blck_%m=0
endif 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) & d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info)
if(info.ne.0) then if(info.ne.0) then
info=4010 info=4010
@ -126,10 +128,11 @@ subroutine mld_zilu_fct(a,l,u,d,info,blck)
return return
contains 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) & d,laspk,lia1,lia2,uaspk,uia1,uia2,l1,l2,info)
implicit none implicit none
integer, intent(in) :: ialg
type(psb_zspmat_type) :: a,b type(psb_zspmat_type) :: a,b
integer :: m,ma,mb,l1,l2,info integer :: m,ma,mb,l1,l2,info
integer, dimension(:) :: lia1,lia2,uia1,uia2 integer, dimension(:) :: lia1,lia2,uia1,uia2
@ -281,9 +284,13 @@ contains
enddo enddo
end if end if
! !
! for milu al=1.; for ilu al=0. ! If we get here we missed the cycle updateloop,
! al = 1.d0 ! which means that this entry does not match; thus
! dia = dia - al*temp*aup(jj) ! we take it out of diagonal for MILU.
!
if (ialg == milu_n_) then
dia = dia - temp*uaspk(jj)
end if
enddo updateloop enddo updateloop
enddo enddo
! !
@ -417,9 +424,13 @@ contains
enddo enddo
end if end if
! !
! for milu al=1.; for ilu al=0. ! If we get here we missed the cycle updateloop,
! al = 1.d0 ! which means that this entry does not match; thus
! dia = dia - al*temp*aup(jj) ! we take it out of diagonal for MILU.
!
if (ialg == milu_n_) then
dia = dia - temp*uaspk(jj)
end if
enddo updateloopb enddo updateloopb
enddo enddo
! !

Write Preview
Loading…
Cancel
Save