mld2p4-newset:


			
			
				stopcriterion
			
			
		
Salvatore Filippone 14 years ago
parent 670e7f1ebd
commit d96b3e4b86

@ -53,8 +53,10 @@ module mld_d_as_smoother
!
type(psb_dspmat_type) :: nd
type(psb_desc_type) :: desc_data
integer :: novr, restr, prol
integer :: novr, restr, prol, nd_nnz_tot
contains
procedure, pass(sm) :: check => d_as_smoother_check
procedure, pass(sm) :: dump => d_as_smoother_dmp
procedure, pass(sm) :: build => d_as_smoother_bld
procedure, pass(sm) :: apply => d_as_smoother_apply
procedure, pass(sm) :: free => d_as_smoother_free
@ -63,13 +65,16 @@ module mld_d_as_smoother
procedure, pass(sm) :: setr => d_as_smoother_setr
procedure, pass(sm) :: descr => d_as_smoother_descr
procedure, pass(sm) :: sizeof => d_as_smoother_sizeof
procedure, pass(sm) :: default => d_as_smoother_default
end type mld_d_as_smoother_type
private :: d_as_smoother_bld, d_as_smoother_apply, &
& d_as_smoother_free, d_as_smoother_seti, &
& d_as_smoother_setc, d_as_smoother_setr,&
& d_as_smoother_descr, d_as_smoother_sizeof
& d_as_smoother_descr, d_as_smoother_sizeof, &
& d_as_smoother_check, d_as_smoother_default,&
& d_as_smoother_dmp
character(len=6), parameter, private :: &
& restrict_names(0:4)=(/'none ','halo ',' ',' ',' '/)
@ -79,6 +84,73 @@ module mld_d_as_smoother
contains
subroutine d_as_smoother_default(sm)
use psb_sparse_mod
Implicit None
! Arguments
class(mld_d_as_smoother_type), intent(inout) :: sm
!!$ sm%restr = psb_halo_
!!$ sm%prol = psb_none_
!!$ sm%novr = 1
!!$
!!$
!!$ if (allocated(sm%sv)) then
!!$ call sm%sv%default()
!!$ end if
!!$
return
end subroutine d_as_smoother_default
subroutine d_as_smoother_check(sm,info)
use psb_sparse_mod
Implicit None
! Arguments
class(mld_d_as_smoother_type), intent(inout) :: sm
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='d_as_smoother_check'
call psb_erractionsave(err_act)
info = psb_success_
!!$
!!$ call mld_check_def(sm%restr,&
!!$ & 'Restrictor',psb_halo_,is_legal_restrict)
!!$ call mld_check_def(sm%prol,&
!!$ & 'Prolongator',psb_none_,is_legal_prolong)
!!$ call mld_check_def(sm%novr,&
!!$ & 'Overlap layers ',0,is_legal_n_ovr)
!!$
!!$
!!$ if (allocated(sm%sv)) then
!!$ call sm%sv%check(info)
!!$ else
!!$ info=3111
!!$ call psb_errpush(info,name)
!!$ goto 9999
!!$ end if
!!$
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_as_smoother_check
subroutine d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info)
use psb_sparse_mod
type(psb_desc_type), intent(in) :: desc_data
@ -100,6 +172,8 @@ contains
call psb_erractionsave(err_act)
info = psb_success_
ictxt = psb_cd_get_context(desc_data)
call psb_info (ictxt,me,np)
trans_ = psb_toupper(trans)
select case(trans_)
@ -250,9 +324,11 @@ contains
goto 9999
end select
write(0,*) me,' Entry to inner slver in AS ',tx
call sm%sv%apply(done,tx,dzero,ty,sm%desc_data,trans_,aux,info)
write(0,*) me,' out from inner slver in AS ',ty
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Error in sub_aply Jacobi Sweeps = 1')
goto 9999
@ -401,12 +477,13 @@ contains
! and Y(j) is the approximate solution at sweep j.
!
ww(1:n_row) = tx(1:n_row)
call psb_spmm(-done,sm%nd,tx,done,ww,sm%desc_data,info,work=aux,trans=trans_)
write(0,*) me,' Entry to spmm in AS-ND',ty
call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,work=aux,trans=trans_)
if (info /= psb_success_) exit
write(0,*) me,' Entry to inner slver in AS-ND',ww
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info)
write(0,*) me,' Exit from inner slver in AS-ND ',ty
if (info /= psb_success_) exit
@ -525,7 +602,7 @@ contains
integer, intent(out) :: info
! Local variables
type(psb_dspmat_type) :: blck, atmp
integer :: n_row,n_col, nrow_a, nhalo, novr, data_
integer :: n_row,n_col, nrow_a, nhalo, novr, data_, nzeros
real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
integer :: ictxt,np,me,i, err_act, debug_unit, debug_level
character(len=20) :: name='d_as_smoother_bld', ch_err
@ -631,6 +708,10 @@ contains
call psb_errpush(psb_err_from_subroutine_,name,a_err='clip & psb_spcnv csr 4')
goto 9999
end if
nzeros = sm%nd%get_nzeros()
write(0,*) me,' ND nzeors ',nzeros
call psb_sum(ictxt,nzeros)
sm%nd_nnz_tot = nzeros
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' end'
@ -871,4 +952,50 @@ contains
return
end function d_as_smoother_sizeof
subroutine d_as_smoother_dmp(sm,ictxt,level,info,prefix,head,smoother,solver)
use psb_sparse_mod
implicit none
class(mld_d_as_smoother_type), intent(in) :: sm
integer, intent(in) :: ictxt,level
integer, intent(out) :: info
character(len=*), intent(in), optional :: prefix, head
logical, optional, intent(in) :: smoother, solver
integer :: i, j, il1, iln, lname, lev
integer :: icontxt,iam, np
character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than
logical :: smoother_
! len of prefix_
info = 0
if (present(prefix)) then
prefix_ = trim(prefix(1:min(len(prefix),len(prefix_))))
else
prefix_ = "dump_smth_d"
end if
call psb_info(ictxt,iam,np)
if (present(smoother)) then
smoother_ = smoother
else
smoother_ = .false.
end if
lname = len_trim(prefix_)
fname = trim(prefix_)
write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
lname = lname + 5
if (smoother_) then
write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_nd.mtx'
if (sm%nd%is_asb()) &
& call sm%nd%print(fname,head=head)
end if
! At base level do nothing for the smoother
if (allocated(sm%sv)) &
& call sm%sv%dump(ictxt,level,info,solver=solver)
end subroutine d_as_smoother_dmp
end module mld_d_as_smoother

@ -53,6 +53,7 @@ module mld_d_ilu_solver
integer :: fact_type, fill_in
real(psb_dpk_) :: thresh
contains
procedure, pass(sv) :: dump => d_ilu_solver_dmp
procedure, pass(sv) :: build => d_ilu_solver_bld
procedure, pass(sv) :: apply => d_ilu_solver_apply
procedure, pass(sv) :: free => d_ilu_solver_free
@ -67,7 +68,8 @@ module mld_d_ilu_solver
private :: d_ilu_solver_bld, d_ilu_solver_apply, &
& d_ilu_solver_free, d_ilu_solver_seti, &
& d_ilu_solver_setc, d_ilu_solver_setr,&
& d_ilu_solver_descr, d_ilu_solver_sizeof
& d_ilu_solver_descr, d_ilu_solver_sizeof, &
& d_ilu_solver_dmp
interface mld_ilu0_fact
@ -605,4 +607,55 @@ contains
return
end function d_ilu_solver_sizeof
subroutine d_ilu_solver_dmp(sv,ictxt,level,info,prefix,head,solver)
use psb_sparse_mod
implicit none
class(mld_d_ilu_solver_type), intent(in) :: sv
integer, intent(in) :: ictxt,level
integer, intent(out) :: info
character(len=*), intent(in), optional :: prefix, head
logical, optional, intent(in) :: solver
integer :: i, j, il1, iln, lname, lev
integer :: icontxt,iam, np
character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than
logical :: solver_
! len of prefix_
info = 0
if (present(prefix)) then
prefix_ = trim(prefix(1:min(len(prefix),len(prefix_))))
else
prefix_ = "dump_slv_d"
end if
call psb_info(ictxt,iam,np)
if (present(solver)) then
solver_ = solver
else
solver_ = .false.
end if
lname = len_trim(prefix_)
fname = trim(prefix_)
write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
lname = lname + 5
if (solver_) then
write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_lower.mtx'
if (sv%l%is_asb()) &
& call sv%l%print(fname,head=head)
write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_diag.mtx'
if (allocated(sv%d)) &
& call psb_geprt(fname,sv%d,head=head)
write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_upper.mtx'
if (sv%u%is_asb()) &
& call sv%u%print(fname,head=head)
end if
end subroutine d_ilu_solver_dmp
end module mld_d_ilu_solver

@ -269,12 +269,15 @@ contains
end select
if (info == psb_success_) call sm%nd%cscnv(info,&
& type='csr',dupl=psb_dupl_add_)
if (info == psb_success_) &
& call sm%sv%build(a,desc_a,upd,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='clip & psb_spcnv csr 4')
goto 9999
end if
call sm%sv%build(a,desc_a,upd,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='solver_build')
goto 9999
end if
nzeros = sm%nd%get_nzeros()
call psb_sum(ictxt,nzeros)
sm%nnz_nd_tot = nzeros

@ -177,6 +177,7 @@ module mld_d_prec_type
type mld_d_base_solver_type
contains
procedure, pass(sv) :: dump => d_base_solver_dmp
procedure, pass(sv) :: build => d_base_solver_bld
procedure, pass(sv) :: apply => d_base_solver_apply
procedure, pass(sv) :: free => d_base_solver_free
@ -192,6 +193,7 @@ module mld_d_prec_type
type mld_d_base_smoother_type
class(mld_d_base_solver_type), allocatable :: sv
contains
procedure, pass(sm) :: dump => d_base_smoother_dmp
procedure, pass(sm) :: build => d_base_smoother_bld
procedure, pass(sm) :: apply => d_base_smoother_apply
procedure, pass(sm) :: free => d_base_smoother_free
@ -221,6 +223,7 @@ module mld_d_prec_type
type(psb_desc_type), pointer :: base_desc => null()
type(psb_dlinmap_type) :: map
contains
procedure, pass(lv) :: dump => d_base_onelev_dump
procedure, pass(lv) :: seti => d_base_onelev_seti
procedure, pass(lv) :: setr => d_base_onelev_setr
procedure, pass(lv) :: setc => d_base_onelev_setc
@ -233,18 +236,21 @@ module mld_d_prec_type
contains
procedure, pass(prec) :: d_apply2v => mld_d_apply2v
procedure, pass(prec) :: d_apply1v => mld_d_apply1v
procedure, pass(prec) :: dump => mld_d_dump
end type mld_dprec_type
private :: d_base_solver_bld, d_base_solver_apply, &
& d_base_solver_free, d_base_solver_seti, &
& d_base_solver_setc, d_base_solver_setr, &
& d_base_solver_descr, d_base_solver_sizeof, &
& d_base_solver_default, &
& d_base_solver_default, d_base_solver_dmp, &
& d_base_smoother_bld, d_base_smoother_apply, &
& d_base_smoother_free, d_base_smoother_seti, &
& d_base_smoother_setc, d_base_smoother_setr,&
& d_base_smoother_descr, d_base_smoother_sizeof, &
& d_base_smoother_default
& d_base_smoother_default, d_base_smoother_dmp, &
& d_base_onelev_dump, d_base_onelev_seti, &
& d_base_onelev_setr, d_base_onelev_setc
!
@ -1405,5 +1411,224 @@ contains
return
end subroutine d_base_onelev_setr
subroutine mld_d_dump(prec,info,istart,iend,prefix,head,ac,smoother,solver)
use psb_sparse_mod
implicit none
class(mld_dprec_type), intent(in) :: prec
integer, intent(out) :: info
integer, intent(in), optional :: istart, iend
character(len=*), intent(in), optional :: prefix, head
logical, optional, intent(in) :: smoother, solver,ac
integer :: i, j, il1, iln, lname, lev
integer :: icontxt,iam, np
character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than
! len of prefix_
info = 0
iln = size(prec%precv)
if (present(istart)) then
il1 = max(1,istart)
else
il1 = 2
end if
if (present(iend)) then
iln = min(iln, iend)
end if
do lev=il1, iln
call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,&
& ac=ac,smoother=smoother,solver=solver)
end do
end subroutine mld_d_dump
subroutine d_base_onelev_dump(lv,level,info,prefix,head,ac,smoother,solver)
use psb_sparse_mod
implicit none
class(mld_donelev_type), intent(in) :: lv
integer, intent(in) :: level
integer, intent(out) :: info
character(len=*), intent(in), optional :: prefix, head
logical, optional, intent(in) :: ac, smoother, solver
integer :: i, j, il1, iln, lname, lev
integer :: icontxt,iam, np
character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than
logical :: ac_
! len of prefix_
info = 0
if (present(prefix)) then
prefix_ = trim(prefix(1:min(len(prefix),len(prefix_))))
else
prefix_ = "dump_lev_d"
end if
if (associated(lv%base_desc)) then
icontxt = psb_cd_get_context(lv%base_desc)
call psb_info(icontxt,iam,np)
else
icontxt = -1
iam = -1
end if
if (present(ac)) then
ac_ = ac
else
ac_ = .false.
end if
lname = len_trim(prefix_)
fname = trim(prefix_)
write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
lname = lname + 5
if (level >= 2) then
write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_ac.mtx'
write(0,*) 'Filename ',fname
if (ac_) call lv%ac%print(fname,head=head)
end if
if (allocated(lv%sm)) &
& call lv%sm%dump(icontxt,level,info,smoother=smoother,solver=solver)
end subroutine d_base_onelev_dump
subroutine d_base_smoother_dmp(sm,ictxt,level,info,prefix,head,smoother,solver)
use psb_sparse_mod
implicit none
class(mld_d_base_smoother_type), intent(in) :: sm
integer, intent(in) :: ictxt,level
integer, intent(out) :: info
character(len=*), intent(in), optional :: prefix, head
logical, optional, intent(in) :: smoother, solver
integer :: i, j, il1, iln, lname, lev
integer :: icontxt,iam, np
character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than
logical :: smoother_
! len of prefix_
info = 0
if (present(prefix)) then
prefix_ = trim(prefix(1:min(len(prefix),len(prefix_))))
else
prefix_ = "dump_smth_d"
end if
call psb_info(ictxt,iam,np)
if (present(smoother)) then
smoother_ = smoother
else
smoother_ = .false.
end if
lname = len_trim(prefix_)
fname = trim(prefix_)
write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
lname = lname + 5
! At base level do nothing for the smoother
if (allocated(sm%sv)) &
& call sm%sv%dump(ictxt,level,info,solver=solver)
end subroutine d_base_smoother_dmp
subroutine d_base_solver_dmp(sv,ictxt,level,info,prefix,head,solver)
use psb_sparse_mod
implicit none
class(mld_d_base_solver_type), intent(in) :: sv
integer, intent(in) :: ictxt,level
integer, intent(out) :: info
character(len=*), intent(in), optional :: prefix, head
logical, optional, intent(in) :: solver
integer :: i, j, il1, iln, lname, lev
integer :: icontxt,iam, np
character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than
logical :: solver_
! len of prefix_
info = 0
if (present(prefix)) then
prefix_ = trim(prefix(1:min(len(prefix),len(prefix_))))
else
prefix_ = "dump_slv_d"
end if
call psb_info(ictxt,iam,np)
if (present(solver)) then
solver_ = solver
else
solver_ = .false.
end if
lname = len_trim(prefix_)
fname = trim(prefix_)
write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
lname = lname + 5
! At base level do nothing for the solver
end subroutine d_base_solver_dmp
!!$
!!$
!!$ subroutine mld_d_precdump_fact(prec,info,istart,iend,prefix,head)
!!$ use psb_base_mod
!!$ implicit none
!!$ type(mld_dprec_type), intent(in) :: prec
!!$ integer, intent(out) :: info
!!$ integer, intent(in), optional :: istart, iend
!!$ character(len=*), intent(in), optional :: prefix,head
!!$ integer :: i, j, il1, iln, lname, lev
!!$ integer :: icontxt,iam, np
!!$ character(len=80) :: prefix_
!!$ character(len=120) :: fname ! len should be at least 20 more than
!!$ ! len of prefix_
!!$
!!$ info = 0
!!$
!!$ if (.not.mld_is_asb(prec)) then
!!$ info = -1
!!$ write(psb_err_unit,*) 'Trying to dump a non-built preconditioner'
!!$ return
!!$ end if
!!$
!!$ il1 = 1
!!$ iln = size(prec%precv)
!!$ if (present(istart)) then
!!$ il1 = max(1,istart)
!!$ end if
!!$ if (present(iend)) then
!!$ iln = min(iln, iend)
!!$ end if
!!$ if (present(prefix)) then
!!$ prefix_ = trim(prefix(1:min(len(prefix),len(prefix_))))
!!$ else
!!$ prefix_ = "dump_fact_d"
!!$ end if
!!$
!!$ icontxt = psb_cd_get_context(prec%precv(1)%prec%desc_data)
!!$ call psb_info(icontxt,iam,np)
!!$ lname = len_trim(prefix_)
!!$ fname = trim(prefix_)
!!$ write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
!!$ lname = lname + 5
!!$ do lev=il1, iln
!!$ write(fname(lname+1:),'(a,i3.3,a)')'_l',lev,'_lower.mtx'
!!$ if (psb_is_asb(prec%precv(lev)%prec%av(mld_l_pr_))) &
!!$ & call psb_csprt(fname,prec%precv(lev)%prec%av(mld_l_pr_),head=head)
!!$ write(fname(lname+1:),'(a,i3.3,a)')'_l',lev,'_diag.mtx'
!!$ if (allocated(prec%precv(lev)%prec%d)) &
!!$ & call psb_geprt(fname,prec%precv(lev)%prec%d,head=head)
!!$ write(fname(lname+1:),'(a,i3.3,a)')'_l',lev,'_upper.mtx'
!!$ if (psb_is_asb(prec%precv(lev)%prec%av(mld_u_pr_))) &
!!$ & call psb_csprt(fname,prec%precv(lev)%prec%av(mld_u_pr_),head=head)
!!$ end do
!!$
!!$ end subroutine mld_d_precdump_fact
end module mld_d_prec_type

@ -514,20 +514,21 @@ contains
! Pre/post-smoothing versions.
! Note that the transpose switches pre <-> post.
!
write(0,*) me,' inner_ml: mult ', level
select case(p%precv(level)%iprcparm(mld_smoother_pos_))
case(mld_post_smooth_)
select case (trans_)
case('N')
write(0,*) me,' inner_ml: post',level
if (level > 1) then
! Apply the restriction
call psb_map_X2Y(done,mlprec_wrk(level-1)%x2l,&
& dzero,mlprec_wrk(level)%x2l,&
& p%precv(level)%map,info,work=work)
write(0,*) me,' inner_ml: entry x2l:',mlprec_wrk(level)%x2l
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
@ -556,12 +557,14 @@ contains
if (info /= psb_success_) goto 9999
sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_post_)
write(0,*) me,' inner_ml: apply ',level,' sweeps: ',sweeps
call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%x2l,done,mlprec_wrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_)
write(0,*) me,' inner_ml: apply ',level,' sweeps: ',sweeps
call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
@ -569,6 +572,8 @@ contains
end if
write(0,*) me,' inner_ml: exit y2l:',mlprec_wrk(level)%y2l
case('T','C')
! Post-smoothing transpose is pre-smoothing

@ -217,6 +217,7 @@ program ppde
end if
tprec = psb_wtime()-t1
call prec%dump(info,prefix='test-ml',ac=.true.,solver=.true.,smoother=.true.)
call psb_amx(ictxt,tprec)

@ -1,9 +1,9 @@
BICGSTAB ! Iterative method: BiCGSTAB BiCG CGS RGMRES BiCGSTABL CG
CSR ! Storage format CSR COO JAD
040 ! IDIM; domain size is idim**3
003 ! IDIM; domain size is idim**3
2 ! ISTOPC
0200 ! ITMAX
-1 ! ITRACE
0002 ! ITMAX
01 ! ITRACE
30 ! IRST (restart for RGMRES and BiCGSTABL)
1.d-6 ! EPS
3L-M-RAS-I-D4 ! Longer descriptive name for preconditioner (up to 20 chars)
@ -12,7 +12,7 @@ ML ! Preconditioner NONE JACOBI BJAC AS ML
HALO ! Restriction operator NONE HALO
NONE ! Prolongation operator NONE SUM AVG
ILU ! Subdomain solver DSCALE ILU MILU ILUT UMF SLU
0 ! Level-set N for ILU(N)
1 ! Level-set N for ILU(N)
1.d-4 ! Threshold T for ILU(T,P)
1 ! Smoother/Jacobi sweeps
AS ! Smoother type JACOBI BJAC AS ignored for non-ML
@ -23,8 +23,8 @@ MULT ! Type of multilevel correction: ADD MULT
POST ! Side of multiplicative correction PRE POST TWOSIDE (ignored for ADD)
DIST ! Coarse level: matrix distribution DIST REPL
BJAC ! Coarse level: solver JACOBI BJAC UMF SLU SLUDIST
UMF ! Coarse level: subsolver DSCALE ILU UMF SLU SLUDIST
1 ! Coarse level: Level-set N for ILU(N)
ILU ! Coarse level: subsolver DSCALE ILU UMF SLU SLUDIST
0 ! Coarse level: Level-set N for ILU(N)
1.d-4 ! Coarse level: Threshold T for ILU(T,P)
4 ! Coarse level: Number of Jacobi sweeps
0.10d0 ! Smoother Aggregation Threshold: >= 0.0

Loading…
Cancel
Save