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

Fixes for support to remapping after merging new context. Needs more testing.

Browse Source
implement-ainv
Salvatore Filippone 4 years ago
parent d14bd31b4a
commit 788211c794
32 changed files with 1816 additions and 1076 deletions
Show Stats Download Patch File Download Diff File

84
amgprec/amg_c_onelev_mod.f90
Unescape Escape View File

@ -157,7 +157,7 @@ module amg_c_onelev_mod
type(psb_cspmat_type) :: ac_pre_remap type(psb_cspmat_type) :: ac_pre_remap
type(psb_desc_type) :: desc_ac_pre_remap type(psb_desc_type) :: desc_ac_pre_remap
integer(psb_ipk_) :: idest integer(psb_ipk_) :: idest
integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:) integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:)
contains contains
procedure, pass(rmp) :: clone => c_remap_data_clone procedure, pass(rmp) :: clone => c_remap_data_clone
end type amg_c_remap_data_type end type amg_c_remap_data_type
@ -704,7 +704,17 @@ contains
info = psb_success_ info = psb_success_
nwv = lv%get_wrksz() nwv = lv%get_wrksz()
if (.not.allocated(lv%wrk)) allocate(lv%wrk,stat=info) if (.not.allocated(lv%wrk)) allocate(lv%wrk,stat=info)
if (info == 0) call lv%wrk%alloc(nwv,lv%base_desc,info,vmold=vmold) if (info == 0) then
if (lv%remap_data%desc_ac_pre_remap%is_asb()) then
!
! Need to fix this, we need two different allocations
!
call lv%wrk%alloc(nwv,lv%base_desc,info,vmold=vmold,&
& desc2=lv%remap_data%desc_ac_pre_remap)
else
call lv%wrk%alloc(nwv,lv%base_desc,info,vmold=vmold)
end if
end if
end subroutine c_base_onelev_allocate_wrk end subroutine c_base_onelev_allocate_wrk
@ -724,7 +734,7 @@ contains
end if end if
end subroutine c_base_onelev_free_wrk end subroutine c_base_onelev_free_wrk
subroutine c_wrk_alloc(wk,nwv,desc,info,vmold) subroutine c_wrk_alloc(wk,nwv,desc,info,vmold, desc2)
use psb_base_mod use psb_base_mod
Implicit None Implicit None
@ -735,25 +745,67 @@ contains
type(psb_desc_type), intent(in) :: desc type(psb_desc_type), intent(in) :: desc
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
class(psb_c_base_vect_type), intent(in), optional :: vmold class(psb_c_base_vect_type), intent(in), optional :: vmold
type(psb_desc_type), intent(in), optional :: desc2
! !
integer(psb_ipk_) :: i integer(psb_ipk_) :: i
info = psb_success_ info = psb_success_
call wk%free(info) call wk%free(info)
call psb_geasb(wk%vx2l,desc,info,& if (present(desc2)) then
& scratch=.true.,mold=vmold) !!$ write(0,*) 'Check on wrk_alloc 2',&
call psb_geasb(wk%vy2l,desc,info,& !!$ & desc2%get_local_rows(), desc%get_local_rows(),&
& scratch=.true.,mold=vmold) !!$ & desc2%get_local_cols(),desc%get_local_cols()
call psb_geasb(wk%vtx,desc,info,& !!$ flush(0)
& scratch=.true.,mold=vmold) if (desc2%get_local_cols()>desc%get_local_cols()) then
call psb_geasb(wk%vty,desc,info,& call psb_geasb(wk%vx2l,desc2,info,&
& scratch=.true.,mold=vmold) & scratch=.true.,mold=vmold)
allocate(wk%wv(nwv),stat=info) call psb_geasb(wk%vy2l,desc2,info,&
do i=1,nwv & scratch=.true.,mold=vmold)
call psb_geasb(wk%wv(i),desc,info,& call psb_geasb(wk%vtx,desc2,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(wk%vty,desc2,info,&
& scratch=.true.,mold=vmold)
allocate(wk%wv(nwv),stat=info)
do i=1,nwv
call psb_geasb(wk%wv(i),desc2,info,&
& scratch=.true.,mold=vmold)
end do
else
!!$ write(0,*) 'Check on wrk_alloc 1.5 ',&
!!$ & desc%get_local_rows(),&
!!$ & desc%get_local_cols()
call psb_geasb(wk%vx2l,desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(wk%vy2l,desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(wk%vtx,desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(wk%vty,desc,info,&
& scratch=.true.,mold=vmold)
allocate(wk%wv(nwv),stat=info)
do i=1,nwv
call psb_geasb(wk%wv(i),desc,info,&
& scratch=.true.,mold=vmold)
end do
end if
else
!!$ write(0,*) 'Check on wrk_alloc 1 ',&
!!$ & desc%get_local_rows(),&
!!$ & desc%get_local_cols()
call psb_geasb(wk%vx2l,desc,info,&
& scratch=.true.,mold=vmold) & scratch=.true.,mold=vmold)
end do call psb_geasb(wk%vy2l,desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(wk%vtx,desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(wk%vty,desc,info,&
& scratch=.true.,mold=vmold)
allocate(wk%wv(nwv),stat=info)
do i=1,nwv
call psb_geasb(wk%wv(i),desc,info,&
& scratch=.true.,mold=vmold)
end do
end if
end subroutine c_wrk_alloc end subroutine c_wrk_alloc
subroutine c_wrk_free(wk,info) subroutine c_wrk_free(wk,info)

10
amgprec/amg_c_prec_type.f90
Unescape Escape View File

@ -738,15 +738,15 @@ contains
character(len=*), intent(in), optional :: prefix, head character(len=*), intent(in), optional :: prefix, head
logical, optional, intent(in) :: smoother, solver,ac, rp, tprol, global_num logical, optional, intent(in) :: smoother, solver,ac, rp, tprol, global_num
integer(psb_ipk_) :: i, j, il1, iln, lev integer(psb_ipk_) :: i, j, il1, iln, lev
type(psb_ctxt_type) :: icontxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: iam, np, iproc_ integer(psb_ipk_) :: iam, np, iproc_
character(len=80) :: prefix_ character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than character(len=120) :: fname ! len should be at least 20 more than
! len of prefix_ ! len of prefix_
info = 0 info = 0
icontxt = prec%ctxt ctxt = prec%ctxt
call psb_info(icontxt,iam,np) call psb_info(ctxt,iam,np)
iln = size(prec%precv) iln = size(prec%precv)
if (present(istart)) then if (present(istart)) then
@ -812,13 +812,13 @@ contains
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
! Local vars ! Local vars
integer(psb_ipk_) :: i, j, ln, lev integer(psb_ipk_) :: i, j, ln, lev
type(psb_ctxt_type) :: icontxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: iam, np integer(psb_ipk_) :: iam, np
info = psb_success_ info = psb_success_
select type(pout => precout) select type(pout => precout)
class is (amg_cprec_type) class is (amg_cprec_type)
pout%ctxt = prec%ctxt pout%ctxt = prec%ctxt
pout%ag_data = prec%ag_data pout%ag_data = prec%ag_data
pout%outer_sweeps = prec%outer_sweeps pout%outer_sweeps = prec%outer_sweeps
if (allocated(prec%precv)) then if (allocated(prec%precv)) then

84
amgprec/amg_d_onelev_mod.f90
Unescape Escape View File

@ -157,7 +157,7 @@ module amg_d_onelev_mod
type(psb_dspmat_type) :: ac_pre_remap type(psb_dspmat_type) :: ac_pre_remap
type(psb_desc_type) :: desc_ac_pre_remap type(psb_desc_type) :: desc_ac_pre_remap
integer(psb_ipk_) :: idest integer(psb_ipk_) :: idest
integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:) integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:)
contains contains
procedure, pass(rmp) :: clone => d_remap_data_clone procedure, pass(rmp) :: clone => d_remap_data_clone
end type amg_d_remap_data_type end type amg_d_remap_data_type
@ -704,7 +704,17 @@ contains
info = psb_success_ info = psb_success_
nwv = lv%get_wrksz() nwv = lv%get_wrksz()
if (.not.allocated(lv%wrk)) allocate(lv%wrk,stat=info) if (.not.allocated(lv%wrk)) allocate(lv%wrk,stat=info)
if (info == 0) call lv%wrk%alloc(nwv,lv%base_desc,info,vmold=vmold) if (info == 0) then
if (lv%remap_data%desc_ac_pre_remap%is_asb()) then
!
! Need to fix this, we need two different allocations
!
call lv%wrk%alloc(nwv,lv%base_desc,info,vmold=vmold,&
& desc2=lv%remap_data%desc_ac_pre_remap)
else
call lv%wrk%alloc(nwv,lv%base_desc,info,vmold=vmold)
end if
end if
end subroutine d_base_onelev_allocate_wrk end subroutine d_base_onelev_allocate_wrk
@ -724,7 +734,7 @@ contains
end if end if
end subroutine d_base_onelev_free_wrk end subroutine d_base_onelev_free_wrk
subroutine d_wrk_alloc(wk,nwv,desc,info,vmold) subroutine d_wrk_alloc(wk,nwv,desc,info,vmold, desc2)
use psb_base_mod use psb_base_mod
Implicit None Implicit None
@ -735,25 +745,67 @@ contains
type(psb_desc_type), intent(in) :: desc type(psb_desc_type), intent(in) :: desc
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
class(psb_d_base_vect_type), intent(in), optional :: vmold class(psb_d_base_vect_type), intent(in), optional :: vmold
type(psb_desc_type), intent(in), optional :: desc2
! !
integer(psb_ipk_) :: i integer(psb_ipk_) :: i
info = psb_success_ info = psb_success_
call wk%free(info) call wk%free(info)
call psb_geasb(wk%vx2l,desc,info,& if (present(desc2)) then
& scratch=.true.,mold=vmold) !!$ write(0,*) 'Check on wrk_alloc 2',&
call psb_geasb(wk%vy2l,desc,info,& !!$ & desc2%get_local_rows(), desc%get_local_rows(),&
& scratch=.true.,mold=vmold) !!$ & desc2%get_local_cols(),desc%get_local_cols()
call psb_geasb(wk%vtx,desc,info,& !!$ flush(0)
& scratch=.true.,mold=vmold) if (desc2%get_local_cols()>desc%get_local_cols()) then
call psb_geasb(wk%vty,desc,info,& call psb_geasb(wk%vx2l,desc2,info,&
& scratch=.true.,mold=vmold) & scratch=.true.,mold=vmold)
allocate(wk%wv(nwv),stat=info) call psb_geasb(wk%vy2l,desc2,info,&
do i=1,nwv & scratch=.true.,mold=vmold)
call psb_geasb(wk%wv(i),desc,info,& call psb_geasb(wk%vtx,desc2,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(wk%vty,desc2,info,&
& scratch=.true.,mold=vmold)
allocate(wk%wv(nwv),stat=info)
do i=1,nwv
call psb_geasb(wk%wv(i),desc2,info,&
& scratch=.true.,mold=vmold)
end do
else
!!$ write(0,*) 'Check on wrk_alloc 1.5 ',&
!!$ & desc%get_local_rows(),&
!!$ & desc%get_local_cols()
call psb_geasb(wk%vx2l,desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(wk%vy2l,desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(wk%vtx,desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(wk%vty,desc,info,&
& scratch=.true.,mold=vmold)
allocate(wk%wv(nwv),stat=info)
do i=1,nwv
call psb_geasb(wk%wv(i),desc,info,&
& scratch=.true.,mold=vmold)
end do
end if
else
!!$ write(0,*) 'Check on wrk_alloc 1 ',&
!!$ & desc%get_local_rows(),&
!!$ & desc%get_local_cols()
call psb_geasb(wk%vx2l,desc,info,&
& scratch=.true.,mold=vmold) & scratch=.true.,mold=vmold)
end do call psb_geasb(wk%vy2l,desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(wk%vtx,desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(wk%vty,desc,info,&
& scratch=.true.,mold=vmold)
allocate(wk%wv(nwv),stat=info)
do i=1,nwv
call psb_geasb(wk%wv(i),desc,info,&
& scratch=.true.,mold=vmold)
end do
end if
end subroutine d_wrk_alloc end subroutine d_wrk_alloc
subroutine d_wrk_free(wk,info) subroutine d_wrk_free(wk,info)

10
amgprec/amg_d_prec_type.f90
Unescape Escape View File

@ -738,15 +738,15 @@ contains
character(len=*), intent(in), optional :: prefix, head character(len=*), intent(in), optional :: prefix, head
logical, optional, intent(in) :: smoother, solver,ac, rp, tprol, global_num logical, optional, intent(in) :: smoother, solver,ac, rp, tprol, global_num
integer(psb_ipk_) :: i, j, il1, iln, lev integer(psb_ipk_) :: i, j, il1, iln, lev
type(psb_ctxt_type) :: icontxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: iam, np, iproc_ integer(psb_ipk_) :: iam, np, iproc_
character(len=80) :: prefix_ character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than character(len=120) :: fname ! len should be at least 20 more than
! len of prefix_ ! len of prefix_
info = 0 info = 0
icontxt = prec%ctxt ctxt = prec%ctxt
call psb_info(icontxt,iam,np) call psb_info(ctxt,iam,np)
iln = size(prec%precv) iln = size(prec%precv)
if (present(istart)) then if (present(istart)) then
@ -812,13 +812,13 @@ contains
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
! Local vars ! Local vars
integer(psb_ipk_) :: i, j, ln, lev integer(psb_ipk_) :: i, j, ln, lev
type(psb_ctxt_type) :: icontxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: iam, np integer(psb_ipk_) :: iam, np
info = psb_success_ info = psb_success_
select type(pout => precout) select type(pout => precout)
class is (amg_dprec_type) class is (amg_dprec_type)
pout%ctxt = prec%ctxt pout%ctxt = prec%ctxt
pout%ag_data = prec%ag_data pout%ag_data = prec%ag_data
pout%outer_sweeps = prec%outer_sweeps pout%outer_sweeps = prec%outer_sweeps
if (allocated(prec%precv)) then if (allocated(prec%precv)) then

84
amgprec/amg_s_onelev_mod.f90
Unescape Escape View File

@ -157,7 +157,7 @@ module amg_s_onelev_mod
type(psb_sspmat_type) :: ac_pre_remap type(psb_sspmat_type) :: ac_pre_remap
type(psb_desc_type) :: desc_ac_pre_remap type(psb_desc_type) :: desc_ac_pre_remap
integer(psb_ipk_) :: idest integer(psb_ipk_) :: idest
integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:) integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:)
contains contains
procedure, pass(rmp) :: clone => s_remap_data_clone procedure, pass(rmp) :: clone => s_remap_data_clone
end type amg_s_remap_data_type end type amg_s_remap_data_type
@ -704,7 +704,17 @@ contains
info = psb_success_ info = psb_success_
nwv = lv%get_wrksz() nwv = lv%get_wrksz()
if (.not.allocated(lv%wrk)) allocate(lv%wrk,stat=info) if (.not.allocated(lv%wrk)) allocate(lv%wrk,stat=info)
if (info == 0) call lv%wrk%alloc(nwv,lv%base_desc,info,vmold=vmold) if (info == 0) then
if (lv%remap_data%desc_ac_pre_remap%is_asb()) then
!
! Need to fix this, we need two different allocations
!
call lv%wrk%alloc(nwv,lv%base_desc,info,vmold=vmold,&
& desc2=lv%remap_data%desc_ac_pre_remap)
else
call lv%wrk%alloc(nwv,lv%base_desc,info,vmold=vmold)
end if
end if
end subroutine s_base_onelev_allocate_wrk end subroutine s_base_onelev_allocate_wrk
@ -724,7 +734,7 @@ contains
end if end if
end subroutine s_base_onelev_free_wrk end subroutine s_base_onelev_free_wrk
subroutine s_wrk_alloc(wk,nwv,desc,info,vmold) subroutine s_wrk_alloc(wk,nwv,desc,info,vmold, desc2)
use psb_base_mod use psb_base_mod
Implicit None Implicit None
@ -735,25 +745,67 @@ contains
type(psb_desc_type), intent(in) :: desc type(psb_desc_type), intent(in) :: desc
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
class(psb_s_base_vect_type), intent(in), optional :: vmold class(psb_s_base_vect_type), intent(in), optional :: vmold
type(psb_desc_type), intent(in), optional :: desc2
! !
integer(psb_ipk_) :: i integer(psb_ipk_) :: i
info = psb_success_ info = psb_success_
call wk%free(info) call wk%free(info)
call psb_geasb(wk%vx2l,desc,info,& if (present(desc2)) then
& scratch=.true.,mold=vmold) !!$ write(0,*) 'Check on wrk_alloc 2',&
call psb_geasb(wk%vy2l,desc,info,& !!$ & desc2%get_local_rows(), desc%get_local_rows(),&
& scratch=.true.,mold=vmold) !!$ & desc2%get_local_cols(),desc%get_local_cols()
call psb_geasb(wk%vtx,desc,info,& !!$ flush(0)
& scratch=.true.,mold=vmold) if (desc2%get_local_cols()>desc%get_local_cols()) then
call psb_geasb(wk%vty,desc,info,& call psb_geasb(wk%vx2l,desc2,info,&
& scratch=.true.,mold=vmold) & scratch=.true.,mold=vmold)
allocate(wk%wv(nwv),stat=info) call psb_geasb(wk%vy2l,desc2,info,&
do i=1,nwv & scratch=.true.,mold=vmold)
call psb_geasb(wk%wv(i),desc,info,& call psb_geasb(wk%vtx,desc2,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(wk%vty,desc2,info,&
& scratch=.true.,mold=vmold)
allocate(wk%wv(nwv),stat=info)
do i=1,nwv
call psb_geasb(wk%wv(i),desc2,info,&
& scratch=.true.,mold=vmold)
end do
else
!!$ write(0,*) 'Check on wrk_alloc 1.5 ',&
!!$ & desc%get_local_rows(),&
!!$ & desc%get_local_cols()
call psb_geasb(wk%vx2l,desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(wk%vy2l,desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(wk%vtx,desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(wk%vty,desc,info,&
& scratch=.true.,mold=vmold)
allocate(wk%wv(nwv),stat=info)
do i=1,nwv
call psb_geasb(wk%wv(i),desc,info,&
& scratch=.true.,mold=vmold)
end do
end if
else
!!$ write(0,*) 'Check on wrk_alloc 1 ',&
!!$ & desc%get_local_rows(),&
!!$ & desc%get_local_cols()
call psb_geasb(wk%vx2l,desc,info,&
& scratch=.true.,mold=vmold) & scratch=.true.,mold=vmold)
end do call psb_geasb(wk%vy2l,desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(wk%vtx,desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(wk%vty,desc,info,&
& scratch=.true.,mold=vmold)
allocate(wk%wv(nwv),stat=info)
do i=1,nwv
call psb_geasb(wk%wv(i),desc,info,&
& scratch=.true.,mold=vmold)
end do
end if
end subroutine s_wrk_alloc end subroutine s_wrk_alloc
subroutine s_wrk_free(wk,info) subroutine s_wrk_free(wk,info)

10
amgprec/amg_s_prec_type.f90
Unescape Escape View File

@ -738,15 +738,15 @@ contains
character(len=*), intent(in), optional :: prefix, head character(len=*), intent(in), optional :: prefix, head
logical, optional, intent(in) :: smoother, solver,ac, rp, tprol, global_num logical, optional, intent(in) :: smoother, solver,ac, rp, tprol, global_num
integer(psb_ipk_) :: i, j, il1, iln, lev integer(psb_ipk_) :: i, j, il1, iln, lev
type(psb_ctxt_type) :: icontxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: iam, np, iproc_ integer(psb_ipk_) :: iam, np, iproc_
character(len=80) :: prefix_ character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than character(len=120) :: fname ! len should be at least 20 more than
! len of prefix_ ! len of prefix_
info = 0 info = 0
icontxt = prec%ctxt ctxt = prec%ctxt
call psb_info(icontxt,iam,np) call psb_info(ctxt,iam,np)
iln = size(prec%precv) iln = size(prec%precv)
if (present(istart)) then if (present(istart)) then
@ -812,13 +812,13 @@ contains
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
! Local vars ! Local vars
integer(psb_ipk_) :: i, j, ln, lev integer(psb_ipk_) :: i, j, ln, lev
type(psb_ctxt_type) :: icontxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: iam, np integer(psb_ipk_) :: iam, np
info = psb_success_ info = psb_success_
select type(pout => precout) select type(pout => precout)
class is (amg_sprec_type) class is (amg_sprec_type)
pout%ctxt = prec%ctxt pout%ctxt = prec%ctxt
pout%ag_data = prec%ag_data pout%ag_data = prec%ag_data
pout%outer_sweeps = prec%outer_sweeps pout%outer_sweeps = prec%outer_sweeps
if (allocated(prec%precv)) then if (allocated(prec%precv)) then

84
amgprec/amg_z_onelev_mod.f90
Unescape Escape View File

@ -157,7 +157,7 @@ module amg_z_onelev_mod
type(psb_zspmat_type) :: ac_pre_remap type(psb_zspmat_type) :: ac_pre_remap
type(psb_desc_type) :: desc_ac_pre_remap type(psb_desc_type) :: desc_ac_pre_remap
integer(psb_ipk_) :: idest integer(psb_ipk_) :: idest
integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:) integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:)
contains contains
procedure, pass(rmp) :: clone => z_remap_data_clone procedure, pass(rmp) :: clone => z_remap_data_clone
end type amg_z_remap_data_type end type amg_z_remap_data_type
@ -704,7 +704,17 @@ contains
info = psb_success_ info = psb_success_
nwv = lv%get_wrksz() nwv = lv%get_wrksz()
if (.not.allocated(lv%wrk)) allocate(lv%wrk,stat=info) if (.not.allocated(lv%wrk)) allocate(lv%wrk,stat=info)
if (info == 0) call lv%wrk%alloc(nwv,lv%base_desc,info,vmold=vmold) if (info == 0) then
if (lv%remap_data%desc_ac_pre_remap%is_asb()) then
!
! Need to fix this, we need two different allocations
!
call lv%wrk%alloc(nwv,lv%base_desc,info,vmold=vmold,&
& desc2=lv%remap_data%desc_ac_pre_remap)
else
call lv%wrk%alloc(nwv,lv%base_desc,info,vmold=vmold)
end if
end if
end subroutine z_base_onelev_allocate_wrk end subroutine z_base_onelev_allocate_wrk
@ -724,7 +734,7 @@ contains
end if end if
end subroutine z_base_onelev_free_wrk end subroutine z_base_onelev_free_wrk
subroutine z_wrk_alloc(wk,nwv,desc,info,vmold) subroutine z_wrk_alloc(wk,nwv,desc,info,vmold, desc2)
use psb_base_mod use psb_base_mod
Implicit None Implicit None
@ -735,25 +745,67 @@ contains
type(psb_desc_type), intent(in) :: desc type(psb_desc_type), intent(in) :: desc
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
class(psb_z_base_vect_type), intent(in), optional :: vmold class(psb_z_base_vect_type), intent(in), optional :: vmold
type(psb_desc_type), intent(in), optional :: desc2
! !
integer(psb_ipk_) :: i integer(psb_ipk_) :: i
info = psb_success_ info = psb_success_
call wk%free(info) call wk%free(info)
call psb_geasb(wk%vx2l,desc,info,& if (present(desc2)) then
& scratch=.true.,mold=vmold) !!$ write(0,*) 'Check on wrk_alloc 2',&
call psb_geasb(wk%vy2l,desc,info,& !!$ & desc2%get_local_rows(), desc%get_local_rows(),&
& scratch=.true.,mold=vmold) !!$ & desc2%get_local_cols(),desc%get_local_cols()
call psb_geasb(wk%vtx,desc,info,& !!$ flush(0)
& scratch=.true.,mold=vmold) if (desc2%get_local_cols()>desc%get_local_cols()) then
call psb_geasb(wk%vty,desc,info,& call psb_geasb(wk%vx2l,desc2,info,&
& scratch=.true.,mold=vmold) & scratch=.true.,mold=vmold)
allocate(wk%wv(nwv),stat=info) call psb_geasb(wk%vy2l,desc2,info,&
do i=1,nwv & scratch=.true.,mold=vmold)
call psb_geasb(wk%wv(i),desc,info,& call psb_geasb(wk%vtx,desc2,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(wk%vty,desc2,info,&
& scratch=.true.,mold=vmold)
allocate(wk%wv(nwv),stat=info)
do i=1,nwv
call psb_geasb(wk%wv(i),desc2,info,&
& scratch=.true.,mold=vmold)
end do
else
!!$ write(0,*) 'Check on wrk_alloc 1.5 ',&
!!$ & desc%get_local_rows(),&
!!$ & desc%get_local_cols()
call psb_geasb(wk%vx2l,desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(wk%vy2l,desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(wk%vtx,desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(wk%vty,desc,info,&
& scratch=.true.,mold=vmold)
allocate(wk%wv(nwv),stat=info)
do i=1,nwv
call psb_geasb(wk%wv(i),desc,info,&
& scratch=.true.,mold=vmold)
end do
end if
else
!!$ write(0,*) 'Check on wrk_alloc 1 ',&
!!$ & desc%get_local_rows(),&
!!$ & desc%get_local_cols()
call psb_geasb(wk%vx2l,desc,info,&
& scratch=.true.,mold=vmold) & scratch=.true.,mold=vmold)
end do call psb_geasb(wk%vy2l,desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(wk%vtx,desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(wk%vty,desc,info,&
& scratch=.true.,mold=vmold)
allocate(wk%wv(nwv),stat=info)
do i=1,nwv
call psb_geasb(wk%wv(i),desc,info,&
& scratch=.true.,mold=vmold)
end do
end if
end subroutine z_wrk_alloc end subroutine z_wrk_alloc
subroutine z_wrk_free(wk,info) subroutine z_wrk_free(wk,info)

10
amgprec/amg_z_prec_type.f90
Unescape Escape View File

@ -738,15 +738,15 @@ contains
character(len=*), intent(in), optional :: prefix, head character(len=*), intent(in), optional :: prefix, head
logical, optional, intent(in) :: smoother, solver,ac, rp, tprol, global_num logical, optional, intent(in) :: smoother, solver,ac, rp, tprol, global_num
integer(psb_ipk_) :: i, j, il1, iln, lev integer(psb_ipk_) :: i, j, il1, iln, lev
type(psb_ctxt_type) :: icontxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: iam, np, iproc_ integer(psb_ipk_) :: iam, np, iproc_
character(len=80) :: prefix_ character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than character(len=120) :: fname ! len should be at least 20 more than
! len of prefix_ ! len of prefix_
info = 0 info = 0
icontxt = prec%ctxt ctxt = prec%ctxt
call psb_info(icontxt,iam,np) call psb_info(ctxt,iam,np)
iln = size(prec%precv) iln = size(prec%precv)
if (present(istart)) then if (present(istart)) then
@ -812,13 +812,13 @@ contains
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
! Local vars ! Local vars
integer(psb_ipk_) :: i, j, ln, lev integer(psb_ipk_) :: i, j, ln, lev
type(psb_ctxt_type) :: icontxt type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: iam, np integer(psb_ipk_) :: iam, np
info = psb_success_ info = psb_success_
select type(pout => precout) select type(pout => precout)
class is (amg_zprec_type) class is (amg_zprec_type)
pout%ctxt = prec%ctxt pout%ctxt = prec%ctxt
pout%ag_data = prec%ag_data pout%ag_data = prec%ag_data
pout%outer_sweeps = prec%outer_sweeps pout%outer_sweeps = prec%outer_sweeps
if (allocated(prec%precv)) then if (allocated(prec%precv)) then

25
amgprec/impl/amg_c_hierarchy_bld.f90
Unescape Escape View File

@ -455,6 +455,31 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info)
end do end do
end if end if
!write(0,*) 'Should we remap? '
if (.true..and.(np>=4)) then
write(0,*) 'Going for remapping '
if (.true.) then
associate(lv=>prec%precv(iszv), rmp => prec%precv(iszv)%remap_data)
call lv%desc_ac%clone(rmp%desc_ac_pre_remap,info)
call lv%ac%clone(rmp%ac_pre_remap,info)
if (np >= 8) then
call psb_remap(np/4,rmp%desc_ac_pre_remap,rmp%ac_pre_remap,&
& rmp%idest,rmp%isrc,rmp%nrsrc,rmp%naggr,lv%desc_ac,lv%ac,info)
else
call psb_remap(np/2,rmp%desc_ac_pre_remap,rmp%ac_pre_remap,&
& rmp%idest,rmp%isrc,rmp%nrsrc,rmp%naggr,lv%desc_ac,lv%ac,info)
end if
write(0,*) me,' Out of remapping ',rmp%desc_ac_pre_remap%get_fmt(),' ',&
& lv%desc_ac%get_fmt(),sum(lv%linmap%naggr),sum(rmp%naggr)
lv%linmap%naggr(:) = rmp%naggr(:)
lv%linmap%p_desc_V => rmp%desc_ac_pre_remap
lv%base_a => lv%ac
lv%base_desc => lv%desc_ac
end associate
end if
end if
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Internal hierarchy build' ) & a_err='Internal hierarchy build' )

6
amgprec/impl/amg_c_smoothers_bld.f90
Unescape Escape View File

@ -117,6 +117,10 @@ subroutine amg_c_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
info = psb_success_ info = psb_success_
ctxt = desc_a%get_context() ctxt = desc_a%get_context()
call psb_info(ctxt, me, np) call psb_info(ctxt, me, np)
if (me <0) then
!!$ write(0,*) 'out of CTXT, should not do anything '
goto 9998
end if
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
@ -289,6 +293,7 @@ subroutine amg_c_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
! !
! build the base preconditioner at level i ! build the base preconditioner at level i
! !
!!$ write(0,*) me,' Building at level ',i
call prec%precv(i)%bld(info,amold=amold,vmold=vmold,imold=imold,ilv=i) call prec%precv(i)%bld(info,amold=amold,vmold=vmold,imold=imold,ilv=i)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -304,6 +309,7 @@ subroutine amg_c_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Exiting with',iszv,' levels' & 'Exiting with',iszv,' levels'
9998 continue
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

338
amgprec/impl/amg_cmlprec_aply.f90
Unescape Escape View File

@ -393,38 +393,39 @@ contains
if(debug_level > 1) then if(debug_level > 1) then
write(debug_unit,*) me,' Start inner_ml_aply at level ',level, info write(debug_unit,*) me,' Start inner_ml_aply at level ',level, info
end if end if
if (me >= 0) then
select case(p%precv(level)%parms%ml_cycle)
select case(p%precv(level)%parms%ml_cycle) case(amg_no_ml_)
!
case(amg_no_ml_) ! No preconditioning, should not really get here
! !
! No preconditioning, should not really get here call psb_errpush(psb_err_internal_error_,name,&
! & a_err='amg_no_ml_ in mlprc_aply?')
call psb_errpush(psb_err_internal_error_,name,& goto 9999
& a_err='amg_no_ml_ in mlprc_aply?')
goto 9999
case(amg_add_ml_) case(amg_add_ml_)
call amg_c_inner_add(p, level, trans, work) call amg_c_inner_add(p, level, trans, work)
case(amg_mult_ml_,amg_vcycle_ml_, amg_wcycle_ml_) case(amg_mult_ml_,amg_vcycle_ml_, amg_wcycle_ml_)
call amg_c_inner_mult(p, level, trans, work) call amg_c_inner_mult(p, level, trans, work)
case(amg_kcycle_ml_, amg_kcyclesym_ml_) case(amg_kcycle_ml_, amg_kcyclesym_ml_)
call amg_c_inner_k_cycle(p, level, trans, work) call amg_c_inner_k_cycle(p, level, trans, work)
case default case default
info = psb_err_from_subroutine_ai_ info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='invalid ml_cycle',& call psb_errpush(info,name,a_err='invalid ml_cycle',&
& i_Err=(/p%precv(level)%parms%ml_cycle,izero,izero,izero,izero/)) & i_Err=(/p%precv(level)%parms%ml_cycle,izero,izero,izero,izero/))
goto 9999 goto 9999
end select end select
if(debug_level > 1) then if(debug_level > 1) then
write(debug_unit,*) me,' End inner_ml_aply at level ',level write(debug_unit,*) me,' End inner_ml_aply at level ',level
end if
end if end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -492,28 +493,30 @@ contains
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv) & wv => p%precv(level)%wrk%wv)
if (allocated(p%precv(level)%sm2a)) then if (me >= 0) then
call psb_geaxpby(cone,vx2l,czero,vy2l,base_desc,info) if (allocated(p%precv(level)%sm2a)) then
call psb_geaxpby(cone,vx2l,czero,vy2l,base_desc,info)
sweeps = max(p%precv(level)%parms%sweeps_pre,p%precv(level)%parms%sweeps_post) sweeps = max(p%precv(level)%parms%sweeps_pre,p%precv(level)%parms%sweeps_post)
do k=1, sweeps do k=1, sweeps
call p%precv(level)%sm%apply(cone,& call p%precv(level)%sm%apply(cone,&
& vy2l,czero,vty,& & vy2l,czero,vty,&
& base_desc, trans,& & base_desc, trans,&
& ione,work,wv,info,init='Z') & ione,work,wv,info,init='Z')
call p%precv(level)%sm2a%apply(cone,& call p%precv(level)%sm2a%apply(cone,&
& vty,czero,vy2l,& & vty,czero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& ione,work,wv,info,init='Z') & ione,work,wv,info,init='Z')
end do end do
else else
sweeps = p%precv(level)%parms%sweeps_pre sweeps = p%precv(level)%parms%sweeps_pre
call p%precv(level)%sm%apply(cone,& call p%precv(level)%sm%apply(cone,&
& vx2l,czero,vy2l,& & vx2l,czero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,wv,info,init='Z') & sweeps,work,wv,info,init='Z')
end if
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -602,7 +605,6 @@ contains
end if end if
ctxt = p%precv(level)%base_desc%get_context() ctxt = p%precv(level)%base_desc%get_context()
call psb_info(ctxt, me, np) call psb_info(ctxt, me, np)
if(debug_level > 1) then if(debug_level > 1) then
write(debug_unit,*) me,' inner_mult at level ',level write(debug_unit,*) me,' inner_mult at level ',level
end if end if
@ -623,22 +625,25 @@ contains
! !
if (pre) then if (pre) then
if (trans == 'N') then if (me >=0) then
sweeps = p%precv(level)%parms%sweeps_pre !!$ write(0,*) me,'Applying smoother pre ', level
if (info == psb_success_) call p%precv(level)%sm%apply(cone,& if (trans == 'N') then
& vx2l,czero,vy2l,base_desc, trans,& sweeps = p%precv(level)%parms%sweeps_pre
& sweeps,work,wv,info,init='Z') if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
else & vx2l,czero,vy2l,base_desc, trans,&
sweeps = p%precv(level)%parms%sweeps_post & sweeps,work,wv,info,init='Z')
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,& else
& vx2l,czero,vy2l, base_desc, trans,& sweeps = p%precv(level)%parms%sweeps_post
& sweeps,work,wv,info,init='Z') if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
end if & vx2l,czero,vy2l, base_desc, trans,&
& sweeps,work,wv,info,init='Z')
if (info /= psb_success_) then end if
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during PRE smoother_apply') if (info /= psb_success_) then
goto 9999 call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during PRE smoother_apply')
goto 9999
end if
end if end if
endif endif
! !
@ -696,11 +701,13 @@ contains
if (p%precv(level)%parms%ml_cycle == amg_wcycle_ml_) then if (p%precv(level)%parms%ml_cycle == amg_wcycle_ml_) then
call psb_geaxpby(cone,vx2l, czero,vty,& if (me >=0) then
& base_desc,info) call psb_geaxpby(cone,vx2l, czero,vty,&
if (info == psb_success_) call psb_spmm(-cone,base_a,& & base_desc,info)
& vy2l,cone,vty,& if (info == psb_success_) call psb_spmm(-cone,base_a,&
& base_desc,info,work=work,trans=trans) & vy2l,cone,vty,&
& base_desc,info,work=work,trans=trans)
end if
if (info == psb_success_) & if (info == psb_success_) &
& call p%precv(level+1)%map_rstr(cone,vty,& & call p%precv(level+1)%map_rstr(cone,vty,&
& czero,p%precv(level+1)%wrk%vx2l,info,work=work,& & czero,p%precv(level+1)%wrk%vx2l,info,work=work,&
@ -728,31 +735,33 @@ contains
if (post) then if (post) then
call psb_geaxpby(cone,vx2l,& if (me >=0) then
& czero,vty,& call psb_geaxpby(cone,vx2l,&
& base_desc,info) & czero,vty,&
if (info == psb_success_) call psb_spmm(-cone,base_a,& & base_desc,info)
& vy2l, cone,vty,base_desc,info,& if (info == psb_success_) call psb_spmm(-cone,base_a,&
& work=work,trans=trans) & vy2l, cone,vty,base_desc,info,&
if (info /= psb_success_) then & work=work,trans=trans)
call psb_errpush(psb_err_internal_error_,name,& if (info /= psb_success_) then
& a_err='Error during residue') call psb_errpush(psb_err_internal_error_,name,&
goto 9999 & a_err='Error during residue')
end if goto 9999
end if
!
! Apply the second smoother !
! ! Apply the second smoother
if (trans == 'N') then !
sweeps = p%precv(level)%parms%sweeps_post if (trans == 'N') then
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,& sweeps = p%precv(level)%parms%sweeps_post
& vty,cone,vy2l, base_desc, trans,& if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
& sweeps,work,wv,info,init='Z') & vty,cone,vy2l, base_desc, trans,&
else & sweeps,work,wv,info,init='Z')
sweeps = p%precv(level)%parms%sweeps_pre else
if (info == psb_success_) call p%precv(level)%sm%apply(cone,& sweeps = p%precv(level)%parms%sweeps_pre
& vty,cone,vy2l, base_desc, trans,& if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& sweeps,work,wv,info,init='Z') & vty,cone,vy2l, base_desc, trans,&
& sweeps,work,wv,info,init='Z')
end if
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -764,11 +773,14 @@ contains
endif endif
else if (level == nlev) then else if (level == nlev) then
!!$ write(0,*) me,'Applying smoother at top level ',psb_errstatus_fatal()
sweeps = p%precv(level)%parms%sweeps_pre if (me >=0) then
if (info == psb_success_) call p%precv(level)%sm%apply(cone,& sweeps = p%precv(level)%parms%sweeps_pre
& vx2l,czero,vy2l,base_desc, trans,& if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& sweeps,work,wv,info) & vx2l,czero,vy2l,base_desc, trans,&
& sweeps,work,wv,info)
end if
!!$ write(0,*) me,' Done applying smoother at top level ',psb_errstatus_fatal()
else else
@ -778,7 +790,7 @@ contains
goto 9999 goto 9999
end if end if
end associate end associate
9998 continue
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -847,52 +859,53 @@ contains
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv(8:)) & wv => p%precv(level)%wrk%wv(8:))
if (level == nlev) then
!
! Apply smoother
!
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& vx2l,czero,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
else if (level < nlev) then if (level == nlev) then
if (me >= 0) then
if (trans == 'N') then !
! Apply smoother
!
sweeps = p%precv(level)%parms%sweeps_pre sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,& if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& vx2l,czero,vy2l,base_desc, trans,& & vx2l,czero,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z') & sweeps,work,wv,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
& vx2l,czero,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-PRE smoother_apply')
goto 9999
end if end if
else if (level < nlev) then
if (me >= 0) then
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& vx2l,czero,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
& vx2l,czero,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-PRE smoother_apply')
goto 9999
end if
! !
! Compute the residual and call recursively ! Compute the residual and call recursively
! !
call psb_geaxpby(cone,vx2l,& call psb_geaxpby(cone,vx2l,&
& czero,vty,& & czero,vty,&
& base_desc,info) & base_desc,info)
if (info == psb_success_) call psb_spmm(-cone,base_a,& if (info == psb_success_) call psb_spmm(-cone,base_a,&
& vy2l,cone,vty,base_desc,info,work=work,trans=trans) & vy2l,cone,vty,base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue') & a_err='Error during residue')
goto 9999 goto 9999
end if
end if end if
! Apply the restriction ! Apply the restriction
call p%precv(level + 1)%map_rstr(cone,vty,& call p%precv(level + 1)%map_rstr(cone,vty,&
& czero,p%precv(level + 1)%wrk%vx2l,& & czero,p%precv(level + 1)%wrk%vx2l,&
@ -940,40 +953,42 @@ contains
& a_err='Error during prolongation') & a_err='Error during prolongation')
goto 9999 goto 9999
end if end if
if (me >= 0) then
!
! Compute the residual
!
call psb_geaxpby(cone,vx2l,&
& czero,vty,base_desc,info)
call psb_spmm(-cone,base_a,vy2l,&
& cone,vty,base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!
! Apply the smoother
!
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
& vty,cone,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& vty,cone,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
end if
! if (info /= psb_success_) then
! Compute the residual call psb_errpush(psb_err_internal_error_,name,&
! & a_err='Error during POST smoother_apply')
call psb_geaxpby(cone,vx2l,& goto 9999
& czero,vty,base_desc,info) end if
call psb_spmm(-cone,base_a,vy2l,&
& cone,vty,base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!
! Apply the smoother
!
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
& vty,cone,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& vty,cone,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
end if end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during POST smoother_apply')
goto 9999
end if
else else
info = psb_err_internal_error_ info = psb_err_internal_error_
@ -991,7 +1006,6 @@ contains
end subroutine amg_c_inner_k_cycle end subroutine amg_c_inner_k_cycle
recursive subroutine amg_cinneritkcycle(p, level, trans, work, innersolv) recursive subroutine amg_cinneritkcycle(p, level, trans, work, innersolv)
use psb_base_mod use psb_base_mod
use amg_prec_mod use amg_prec_mod

25
amgprec/impl/amg_d_hierarchy_bld.f90
Unescape Escape View File

@ -455,6 +455,31 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info)
end do end do
end if end if
!write(0,*) 'Should we remap? '
if (.true..and.(np>=4)) then
write(0,*) 'Going for remapping '
if (.true.) then
associate(lv=>prec%precv(iszv), rmp => prec%precv(iszv)%remap_data)
call lv%desc_ac%clone(rmp%desc_ac_pre_remap,info)
call lv%ac%clone(rmp%ac_pre_remap,info)
if (np >= 8) then
call psb_remap(np/4,rmp%desc_ac_pre_remap,rmp%ac_pre_remap,&
& rmp%idest,rmp%isrc,rmp%nrsrc,rmp%naggr,lv%desc_ac,lv%ac,info)
else
call psb_remap(np/2,rmp%desc_ac_pre_remap,rmp%ac_pre_remap,&
& rmp%idest,rmp%isrc,rmp%nrsrc,rmp%naggr,lv%desc_ac,lv%ac,info)
end if
write(0,*) me,' Out of remapping ',rmp%desc_ac_pre_remap%get_fmt(),' ',&
& lv%desc_ac%get_fmt(),sum(lv%linmap%naggr),sum(rmp%naggr)
lv%linmap%naggr(:) = rmp%naggr(:)
lv%linmap%p_desc_V => rmp%desc_ac_pre_remap
lv%base_a => lv%ac
lv%base_desc => lv%desc_ac
end associate
end if
end if
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Internal hierarchy build' ) & a_err='Internal hierarchy build' )

6
amgprec/impl/amg_d_smoothers_bld.f90
Unescape Escape View File

@ -117,6 +117,10 @@ subroutine amg_d_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
info = psb_success_ info = psb_success_
ctxt = desc_a%get_context() ctxt = desc_a%get_context()
call psb_info(ctxt, me, np) call psb_info(ctxt, me, np)
if (me <0) then
!!$ write(0,*) 'out of CTXT, should not do anything '
goto 9998
end if
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
@ -289,6 +293,7 @@ subroutine amg_d_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
! !
! build the base preconditioner at level i ! build the base preconditioner at level i
! !
!!$ write(0,*) me,' Building at level ',i
call prec%precv(i)%bld(info,amold=amold,vmold=vmold,imold=imold,ilv=i) call prec%precv(i)%bld(info,amold=amold,vmold=vmold,imold=imold,ilv=i)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -304,6 +309,7 @@ subroutine amg_d_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Exiting with',iszv,' levels' & 'Exiting with',iszv,' levels'
9998 continue
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

338
amgprec/impl/amg_dmlprec_aply.f90
Unescape Escape View File

@ -393,38 +393,39 @@ contains
if(debug_level > 1) then if(debug_level > 1) then
write(debug_unit,*) me,' Start inner_ml_aply at level ',level, info write(debug_unit,*) me,' Start inner_ml_aply at level ',level, info
end if end if
if (me >= 0) then
select case(p%precv(level)%parms%ml_cycle)
select case(p%precv(level)%parms%ml_cycle) case(amg_no_ml_)
!
case(amg_no_ml_) ! No preconditioning, should not really get here
! !
! No preconditioning, should not really get here call psb_errpush(psb_err_internal_error_,name,&
! & a_err='amg_no_ml_ in mlprc_aply?')
call psb_errpush(psb_err_internal_error_,name,& goto 9999
& a_err='amg_no_ml_ in mlprc_aply?')
goto 9999
case(amg_add_ml_) case(amg_add_ml_)
call amg_d_inner_add(p, level, trans, work) call amg_d_inner_add(p, level, trans, work)
case(amg_mult_ml_,amg_vcycle_ml_, amg_wcycle_ml_) case(amg_mult_ml_,amg_vcycle_ml_, amg_wcycle_ml_)
call amg_d_inner_mult(p, level, trans, work) call amg_d_inner_mult(p, level, trans, work)
case(amg_kcycle_ml_, amg_kcyclesym_ml_) case(amg_kcycle_ml_, amg_kcyclesym_ml_)
call amg_d_inner_k_cycle(p, level, trans, work) call amg_d_inner_k_cycle(p, level, trans, work)
case default case default
info = psb_err_from_subroutine_ai_ info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='invalid ml_cycle',& call psb_errpush(info,name,a_err='invalid ml_cycle',&
& i_Err=(/p%precv(level)%parms%ml_cycle,izero,izero,izero,izero/)) & i_Err=(/p%precv(level)%parms%ml_cycle,izero,izero,izero,izero/))
goto 9999 goto 9999
end select end select
if(debug_level > 1) then if(debug_level > 1) then
write(debug_unit,*) me,' End inner_ml_aply at level ',level write(debug_unit,*) me,' End inner_ml_aply at level ',level
end if
end if end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -492,28 +493,30 @@ contains
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv) & wv => p%precv(level)%wrk%wv)
if (allocated(p%precv(level)%sm2a)) then if (me >= 0) then
call psb_geaxpby(done,vx2l,dzero,vy2l,base_desc,info) if (allocated(p%precv(level)%sm2a)) then
call psb_geaxpby(done,vx2l,dzero,vy2l,base_desc,info)
sweeps = max(p%precv(level)%parms%sweeps_pre,p%precv(level)%parms%sweeps_post) sweeps = max(p%precv(level)%parms%sweeps_pre,p%precv(level)%parms%sweeps_post)
do k=1, sweeps do k=1, sweeps
call p%precv(level)%sm%apply(done,& call p%precv(level)%sm%apply(done,&
& vy2l,dzero,vty,& & vy2l,dzero,vty,&
& base_desc, trans,& & base_desc, trans,&
& ione,work,wv,info,init='Z') & ione,work,wv,info,init='Z')
call p%precv(level)%sm2a%apply(done,& call p%precv(level)%sm2a%apply(done,&
& vty,dzero,vy2l,& & vty,dzero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& ione,work,wv,info,init='Z') & ione,work,wv,info,init='Z')
end do end do
else else
sweeps = p%precv(level)%parms%sweeps_pre sweeps = p%precv(level)%parms%sweeps_pre
call p%precv(level)%sm%apply(done,& call p%precv(level)%sm%apply(done,&
& vx2l,dzero,vy2l,& & vx2l,dzero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,wv,info,init='Z') & sweeps,work,wv,info,init='Z')
end if
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -602,7 +605,6 @@ contains
end if end if
ctxt = p%precv(level)%base_desc%get_context() ctxt = p%precv(level)%base_desc%get_context()
call psb_info(ctxt, me, np) call psb_info(ctxt, me, np)
if(debug_level > 1) then if(debug_level > 1) then
write(debug_unit,*) me,' inner_mult at level ',level write(debug_unit,*) me,' inner_mult at level ',level
end if end if
@ -623,22 +625,25 @@ contains
! !
if (pre) then if (pre) then
if (trans == 'N') then if (me >=0) then
sweeps = p%precv(level)%parms%sweeps_pre !!$ write(0,*) me,'Applying smoother pre ', level
if (info == psb_success_) call p%precv(level)%sm%apply(done,& if (trans == 'N') then
& vx2l,dzero,vy2l,base_desc, trans,& sweeps = p%precv(level)%parms%sweeps_pre
& sweeps,work,wv,info,init='Z') if (info == psb_success_) call p%precv(level)%sm%apply(done,&
else & vx2l,dzero,vy2l,base_desc, trans,&
sweeps = p%precv(level)%parms%sweeps_post & sweeps,work,wv,info,init='Z')
if (info == psb_success_) call p%precv(level)%sm2%apply(done,& else
& vx2l,dzero,vy2l, base_desc, trans,& sweeps = p%precv(level)%parms%sweeps_post
& sweeps,work,wv,info,init='Z') if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
end if & vx2l,dzero,vy2l, base_desc, trans,&
& sweeps,work,wv,info,init='Z')
if (info /= psb_success_) then end if
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during PRE smoother_apply') if (info /= psb_success_) then
goto 9999 call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during PRE smoother_apply')
goto 9999
end if
end if end if
endif endif
! !
@ -696,11 +701,13 @@ contains
if (p%precv(level)%parms%ml_cycle == amg_wcycle_ml_) then if (p%precv(level)%parms%ml_cycle == amg_wcycle_ml_) then
call psb_geaxpby(done,vx2l, dzero,vty,& if (me >=0) then
& base_desc,info) call psb_geaxpby(done,vx2l, dzero,vty,&
if (info == psb_success_) call psb_spmm(-done,base_a,& & base_desc,info)
& vy2l,done,vty,& if (info == psb_success_) call psb_spmm(-done,base_a,&
& base_desc,info,work=work,trans=trans) & vy2l,done,vty,&
& base_desc,info,work=work,trans=trans)
end if
if (info == psb_success_) & if (info == psb_success_) &
& call p%precv(level+1)%map_rstr(done,vty,& & call p%precv(level+1)%map_rstr(done,vty,&
& dzero,p%precv(level+1)%wrk%vx2l,info,work=work,& & dzero,p%precv(level+1)%wrk%vx2l,info,work=work,&
@ -728,31 +735,33 @@ contains
if (post) then if (post) then
call psb_geaxpby(done,vx2l,& if (me >=0) then
& dzero,vty,& call psb_geaxpby(done,vx2l,&
& base_desc,info) & dzero,vty,&
if (info == psb_success_) call psb_spmm(-done,base_a,& & base_desc,info)
& vy2l, done,vty,base_desc,info,& if (info == psb_success_) call psb_spmm(-done,base_a,&
& work=work,trans=trans) & vy2l, done,vty,base_desc,info,&
if (info /= psb_success_) then & work=work,trans=trans)
call psb_errpush(psb_err_internal_error_,name,& if (info /= psb_success_) then
& a_err='Error during residue') call psb_errpush(psb_err_internal_error_,name,&
goto 9999 & a_err='Error during residue')
end if goto 9999
end if
!
! Apply the second smoother !
! ! Apply the second smoother
if (trans == 'N') then !
sweeps = p%precv(level)%parms%sweeps_post if (trans == 'N') then
if (info == psb_success_) call p%precv(level)%sm2%apply(done,& sweeps = p%precv(level)%parms%sweeps_post
& vty,done,vy2l, base_desc, trans,& if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& sweeps,work,wv,info,init='Z') & vty,done,vy2l, base_desc, trans,&
else & sweeps,work,wv,info,init='Z')
sweeps = p%precv(level)%parms%sweeps_pre else
if (info == psb_success_) call p%precv(level)%sm%apply(done,& sweeps = p%precv(level)%parms%sweeps_pre
& vty,done,vy2l, base_desc, trans,& if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& sweeps,work,wv,info,init='Z') & vty,done,vy2l, base_desc, trans,&
& sweeps,work,wv,info,init='Z')
end if
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -764,11 +773,14 @@ contains
endif endif
else if (level == nlev) then else if (level == nlev) then
!!$ write(0,*) me,'Applying smoother at top level ',psb_errstatus_fatal()
sweeps = p%precv(level)%parms%sweeps_pre if (me >=0) then
if (info == psb_success_) call p%precv(level)%sm%apply(done,& sweeps = p%precv(level)%parms%sweeps_pre
& vx2l,dzero,vy2l,base_desc, trans,& if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& sweeps,work,wv,info) & vx2l,dzero,vy2l,base_desc, trans,&
& sweeps,work,wv,info)
end if
!!$ write(0,*) me,' Done applying smoother at top level ',psb_errstatus_fatal()
else else
@ -778,7 +790,7 @@ contains
goto 9999 goto 9999
end if end if
end associate end associate
9998 continue
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -847,52 +859,53 @@ contains
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv(8:)) & wv => p%precv(level)%wrk%wv(8:))
if (level == nlev) then
!
! Apply smoother
!
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& vx2l,dzero,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
else if (level < nlev) then if (level == nlev) then
if (me >= 0) then
if (trans == 'N') then !
! Apply smoother
!
sweeps = p%precv(level)%parms%sweeps_pre sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,& if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& vx2l,dzero,vy2l,base_desc, trans,& & vx2l,dzero,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z') & sweeps,work,wv,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& vx2l,dzero,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-PRE smoother_apply')
goto 9999
end if end if
else if (level < nlev) then
if (me >= 0) then
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& vx2l,dzero,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& vx2l,dzero,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-PRE smoother_apply')
goto 9999
end if
! !
! Compute the residual and call recursively ! Compute the residual and call recursively
! !
call psb_geaxpby(done,vx2l,& call psb_geaxpby(done,vx2l,&
& dzero,vty,& & dzero,vty,&
& base_desc,info) & base_desc,info)
if (info == psb_success_) call psb_spmm(-done,base_a,& if (info == psb_success_) call psb_spmm(-done,base_a,&
& vy2l,done,vty,base_desc,info,work=work,trans=trans) & vy2l,done,vty,base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue') & a_err='Error during residue')
goto 9999 goto 9999
end if
end if end if
! Apply the restriction ! Apply the restriction
call p%precv(level + 1)%map_rstr(done,vty,& call p%precv(level + 1)%map_rstr(done,vty,&
& dzero,p%precv(level + 1)%wrk%vx2l,& & dzero,p%precv(level + 1)%wrk%vx2l,&
@ -940,40 +953,42 @@ contains
& a_err='Error during prolongation') & a_err='Error during prolongation')
goto 9999 goto 9999
end if end if
if (me >= 0) then
!
! Compute the residual
!
call psb_geaxpby(done,vx2l,&
& dzero,vty,base_desc,info)
call psb_spmm(-done,base_a,vy2l,&
& done,vty,base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!
! Apply the smoother
!
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& vty,done,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& vty,done,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
end if
! if (info /= psb_success_) then
! Compute the residual call psb_errpush(psb_err_internal_error_,name,&
! & a_err='Error during POST smoother_apply')
call psb_geaxpby(done,vx2l,& goto 9999
& dzero,vty,base_desc,info) end if
call psb_spmm(-done,base_a,vy2l,&
& done,vty,base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!
! Apply the smoother
!
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& vty,done,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& vty,done,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
end if end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during POST smoother_apply')
goto 9999
end if
else else
info = psb_err_internal_error_ info = psb_err_internal_error_
@ -991,7 +1006,6 @@ contains
end subroutine amg_d_inner_k_cycle end subroutine amg_d_inner_k_cycle
recursive subroutine amg_dinneritkcycle(p, level, trans, work, innersolv) recursive subroutine amg_dinneritkcycle(p, level, trans, work, innersolv)
use psb_base_mod use psb_base_mod
use amg_prec_mod use amg_prec_mod

25
amgprec/impl/amg_s_hierarchy_bld.f90
Unescape Escape View File

@ -455,6 +455,31 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info)
end do end do
end if end if
!write(0,*) 'Should we remap? '
if (.true..and.(np>=4)) then
write(0,*) 'Going for remapping '
if (.true.) then
associate(lv=>prec%precv(iszv), rmp => prec%precv(iszv)%remap_data)
call lv%desc_ac%clone(rmp%desc_ac_pre_remap,info)
call lv%ac%clone(rmp%ac_pre_remap,info)
if (np >= 8) then
call psb_remap(np/4,rmp%desc_ac_pre_remap,rmp%ac_pre_remap,&
& rmp%idest,rmp%isrc,rmp%nrsrc,rmp%naggr,lv%desc_ac,lv%ac,info)
else
call psb_remap(np/2,rmp%desc_ac_pre_remap,rmp%ac_pre_remap,&
& rmp%idest,rmp%isrc,rmp%nrsrc,rmp%naggr,lv%desc_ac,lv%ac,info)
end if
write(0,*) me,' Out of remapping ',rmp%desc_ac_pre_remap%get_fmt(),' ',&
& lv%desc_ac%get_fmt(),sum(lv%linmap%naggr),sum(rmp%naggr)
lv%linmap%naggr(:) = rmp%naggr(:)
lv%linmap%p_desc_V => rmp%desc_ac_pre_remap
lv%base_a => lv%ac
lv%base_desc => lv%desc_ac
end associate
end if
end if
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Internal hierarchy build' ) & a_err='Internal hierarchy build' )

6
amgprec/impl/amg_s_smoothers_bld.f90
Unescape Escape View File

@ -117,6 +117,10 @@ subroutine amg_s_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
info = psb_success_ info = psb_success_
ctxt = desc_a%get_context() ctxt = desc_a%get_context()
call psb_info(ctxt, me, np) call psb_info(ctxt, me, np)
if (me <0) then
!!$ write(0,*) 'out of CTXT, should not do anything '
goto 9998
end if
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
@ -289,6 +293,7 @@ subroutine amg_s_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
! !
! build the base preconditioner at level i ! build the base preconditioner at level i
! !
!!$ write(0,*) me,' Building at level ',i
call prec%precv(i)%bld(info,amold=amold,vmold=vmold,imold=imold,ilv=i) call prec%precv(i)%bld(info,amold=amold,vmold=vmold,imold=imold,ilv=i)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -304,6 +309,7 @@ subroutine amg_s_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Exiting with',iszv,' levels' & 'Exiting with',iszv,' levels'
9998 continue
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

338
amgprec/impl/amg_smlprec_aply.f90
Unescape Escape View File

@ -393,38 +393,39 @@ contains
if(debug_level > 1) then if(debug_level > 1) then
write(debug_unit,*) me,' Start inner_ml_aply at level ',level, info write(debug_unit,*) me,' Start inner_ml_aply at level ',level, info
end if end if
if (me >= 0) then
select case(p%precv(level)%parms%ml_cycle)
select case(p%precv(level)%parms%ml_cycle) case(amg_no_ml_)
!
case(amg_no_ml_) ! No preconditioning, should not really get here
! !
! No preconditioning, should not really get here call psb_errpush(psb_err_internal_error_,name,&
! & a_err='amg_no_ml_ in mlprc_aply?')
call psb_errpush(psb_err_internal_error_,name,& goto 9999
& a_err='amg_no_ml_ in mlprc_aply?')
goto 9999
case(amg_add_ml_) case(amg_add_ml_)
call amg_s_inner_add(p, level, trans, work) call amg_s_inner_add(p, level, trans, work)
case(amg_mult_ml_,amg_vcycle_ml_, amg_wcycle_ml_) case(amg_mult_ml_,amg_vcycle_ml_, amg_wcycle_ml_)
call amg_s_inner_mult(p, level, trans, work) call amg_s_inner_mult(p, level, trans, work)
case(amg_kcycle_ml_, amg_kcyclesym_ml_) case(amg_kcycle_ml_, amg_kcyclesym_ml_)
call amg_s_inner_k_cycle(p, level, trans, work) call amg_s_inner_k_cycle(p, level, trans, work)
case default case default
info = psb_err_from_subroutine_ai_ info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='invalid ml_cycle',& call psb_errpush(info,name,a_err='invalid ml_cycle',&
& i_Err=(/p%precv(level)%parms%ml_cycle,izero,izero,izero,izero/)) & i_Err=(/p%precv(level)%parms%ml_cycle,izero,izero,izero,izero/))
goto 9999 goto 9999
end select end select
if(debug_level > 1) then if(debug_level > 1) then
write(debug_unit,*) me,' End inner_ml_aply at level ',level write(debug_unit,*) me,' End inner_ml_aply at level ',level
end if
end if end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -492,28 +493,30 @@ contains
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv) & wv => p%precv(level)%wrk%wv)
if (allocated(p%precv(level)%sm2a)) then if (me >= 0) then
call psb_geaxpby(sone,vx2l,szero,vy2l,base_desc,info) if (allocated(p%precv(level)%sm2a)) then
call psb_geaxpby(sone,vx2l,szero,vy2l,base_desc,info)
sweeps = max(p%precv(level)%parms%sweeps_pre,p%precv(level)%parms%sweeps_post) sweeps = max(p%precv(level)%parms%sweeps_pre,p%precv(level)%parms%sweeps_post)
do k=1, sweeps do k=1, sweeps
call p%precv(level)%sm%apply(sone,& call p%precv(level)%sm%apply(sone,&
& vy2l,szero,vty,& & vy2l,szero,vty,&
& base_desc, trans,& & base_desc, trans,&
& ione,work,wv,info,init='Z') & ione,work,wv,info,init='Z')
call p%precv(level)%sm2a%apply(sone,& call p%precv(level)%sm2a%apply(sone,&
& vty,szero,vy2l,& & vty,szero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& ione,work,wv,info,init='Z') & ione,work,wv,info,init='Z')
end do end do
else else
sweeps = p%precv(level)%parms%sweeps_pre sweeps = p%precv(level)%parms%sweeps_pre
call p%precv(level)%sm%apply(sone,& call p%precv(level)%sm%apply(sone,&
& vx2l,szero,vy2l,& & vx2l,szero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,wv,info,init='Z') & sweeps,work,wv,info,init='Z')
end if
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -602,7 +605,6 @@ contains
end if end if
ctxt = p%precv(level)%base_desc%get_context() ctxt = p%precv(level)%base_desc%get_context()
call psb_info(ctxt, me, np) call psb_info(ctxt, me, np)
if(debug_level > 1) then if(debug_level > 1) then
write(debug_unit,*) me,' inner_mult at level ',level write(debug_unit,*) me,' inner_mult at level ',level
end if end if
@ -623,22 +625,25 @@ contains
! !
if (pre) then if (pre) then
if (trans == 'N') then if (me >=0) then
sweeps = p%precv(level)%parms%sweeps_pre !!$ write(0,*) me,'Applying smoother pre ', level
if (info == psb_success_) call p%precv(level)%sm%apply(sone,& if (trans == 'N') then
& vx2l,szero,vy2l,base_desc, trans,& sweeps = p%precv(level)%parms%sweeps_pre
& sweeps,work,wv,info,init='Z') if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
else & vx2l,szero,vy2l,base_desc, trans,&
sweeps = p%precv(level)%parms%sweeps_post & sweeps,work,wv,info,init='Z')
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,& else
& vx2l,szero,vy2l, base_desc, trans,& sweeps = p%precv(level)%parms%sweeps_post
& sweeps,work,wv,info,init='Z') if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
end if & vx2l,szero,vy2l, base_desc, trans,&
& sweeps,work,wv,info,init='Z')
if (info /= psb_success_) then end if
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during PRE smoother_apply') if (info /= psb_success_) then
goto 9999 call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during PRE smoother_apply')
goto 9999
end if
end if end if
endif endif
! !
@ -696,11 +701,13 @@ contains
if (p%precv(level)%parms%ml_cycle == amg_wcycle_ml_) then if (p%precv(level)%parms%ml_cycle == amg_wcycle_ml_) then
call psb_geaxpby(sone,vx2l, szero,vty,& if (me >=0) then
& base_desc,info) call psb_geaxpby(sone,vx2l, szero,vty,&
if (info == psb_success_) call psb_spmm(-sone,base_a,& & base_desc,info)
& vy2l,sone,vty,& if (info == psb_success_) call psb_spmm(-sone,base_a,&
& base_desc,info,work=work,trans=trans) & vy2l,sone,vty,&
& base_desc,info,work=work,trans=trans)
end if
if (info == psb_success_) & if (info == psb_success_) &
& call p%precv(level+1)%map_rstr(sone,vty,& & call p%precv(level+1)%map_rstr(sone,vty,&
& szero,p%precv(level+1)%wrk%vx2l,info,work=work,& & szero,p%precv(level+1)%wrk%vx2l,info,work=work,&
@ -728,31 +735,33 @@ contains
if (post) then if (post) then
call psb_geaxpby(sone,vx2l,& if (me >=0) then
& szero,vty,& call psb_geaxpby(sone,vx2l,&
& base_desc,info) & szero,vty,&
if (info == psb_success_) call psb_spmm(-sone,base_a,& & base_desc,info)
& vy2l, sone,vty,base_desc,info,& if (info == psb_success_) call psb_spmm(-sone,base_a,&
& work=work,trans=trans) & vy2l, sone,vty,base_desc,info,&
if (info /= psb_success_) then & work=work,trans=trans)
call psb_errpush(psb_err_internal_error_,name,& if (info /= psb_success_) then
& a_err='Error during residue') call psb_errpush(psb_err_internal_error_,name,&
goto 9999 & a_err='Error during residue')
end if goto 9999
end if
!
! Apply the second smoother !
! ! Apply the second smoother
if (trans == 'N') then !
sweeps = p%precv(level)%parms%sweeps_post if (trans == 'N') then
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,& sweeps = p%precv(level)%parms%sweeps_post
& vty,sone,vy2l, base_desc, trans,& if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
& sweeps,work,wv,info,init='Z') & vty,sone,vy2l, base_desc, trans,&
else & sweeps,work,wv,info,init='Z')
sweeps = p%precv(level)%parms%sweeps_pre else
if (info == psb_success_) call p%precv(level)%sm%apply(sone,& sweeps = p%precv(level)%parms%sweeps_pre
& vty,sone,vy2l, base_desc, trans,& if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& sweeps,work,wv,info,init='Z') & vty,sone,vy2l, base_desc, trans,&
& sweeps,work,wv,info,init='Z')
end if
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -764,11 +773,14 @@ contains
endif endif
else if (level == nlev) then else if (level == nlev) then
!!$ write(0,*) me,'Applying smoother at top level ',psb_errstatus_fatal()
sweeps = p%precv(level)%parms%sweeps_pre if (me >=0) then
if (info == psb_success_) call p%precv(level)%sm%apply(sone,& sweeps = p%precv(level)%parms%sweeps_pre
& vx2l,szero,vy2l,base_desc, trans,& if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& sweeps,work,wv,info) & vx2l,szero,vy2l,base_desc, trans,&
& sweeps,work,wv,info)
end if
!!$ write(0,*) me,' Done applying smoother at top level ',psb_errstatus_fatal()
else else
@ -778,7 +790,7 @@ contains
goto 9999 goto 9999
end if end if
end associate end associate
9998 continue
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -847,52 +859,53 @@ contains
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv(8:)) & wv => p%precv(level)%wrk%wv(8:))
if (level == nlev) then
!
! Apply smoother
!
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& vx2l,szero,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
else if (level < nlev) then if (level == nlev) then
if (me >= 0) then
if (trans == 'N') then !
! Apply smoother
!
sweeps = p%precv(level)%parms%sweeps_pre sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,& if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& vx2l,szero,vy2l,base_desc, trans,& & vx2l,szero,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z') & sweeps,work,wv,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
& vx2l,szero,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-PRE smoother_apply')
goto 9999
end if end if
else if (level < nlev) then
if (me >= 0) then
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& vx2l,szero,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
& vx2l,szero,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-PRE smoother_apply')
goto 9999
end if
! !
! Compute the residual and call recursively ! Compute the residual and call recursively
! !
call psb_geaxpby(sone,vx2l,& call psb_geaxpby(sone,vx2l,&
& szero,vty,& & szero,vty,&
& base_desc,info) & base_desc,info)
if (info == psb_success_) call psb_spmm(-sone,base_a,& if (info == psb_success_) call psb_spmm(-sone,base_a,&
& vy2l,sone,vty,base_desc,info,work=work,trans=trans) & vy2l,sone,vty,base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue') & a_err='Error during residue')
goto 9999 goto 9999
end if
end if end if
! Apply the restriction ! Apply the restriction
call p%precv(level + 1)%map_rstr(sone,vty,& call p%precv(level + 1)%map_rstr(sone,vty,&
& szero,p%precv(level + 1)%wrk%vx2l,& & szero,p%precv(level + 1)%wrk%vx2l,&
@ -940,40 +953,42 @@ contains
& a_err='Error during prolongation') & a_err='Error during prolongation')
goto 9999 goto 9999
end if end if
if (me >= 0) then
!
! Compute the residual
!
call psb_geaxpby(sone,vx2l,&
& szero,vty,base_desc,info)
call psb_spmm(-sone,base_a,vy2l,&
& sone,vty,base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!
! Apply the smoother
!
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
& vty,sone,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& vty,sone,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
end if
! if (info /= psb_success_) then
! Compute the residual call psb_errpush(psb_err_internal_error_,name,&
! & a_err='Error during POST smoother_apply')
call psb_geaxpby(sone,vx2l,& goto 9999
& szero,vty,base_desc,info) end if
call psb_spmm(-sone,base_a,vy2l,&
& sone,vty,base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!
! Apply the smoother
!
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
& vty,sone,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& vty,sone,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
end if end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during POST smoother_apply')
goto 9999
end if
else else
info = psb_err_internal_error_ info = psb_err_internal_error_
@ -991,7 +1006,6 @@ contains
end subroutine amg_s_inner_k_cycle end subroutine amg_s_inner_k_cycle
recursive subroutine amg_sinneritkcycle(p, level, trans, work, innersolv) recursive subroutine amg_sinneritkcycle(p, level, trans, work, innersolv)
use psb_base_mod use psb_base_mod
use amg_prec_mod use amg_prec_mod

25
amgprec/impl/amg_z_hierarchy_bld.f90
Unescape Escape View File

@ -455,6 +455,31 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info)
end do end do
end if end if
!write(0,*) 'Should we remap? '
if (.true..and.(np>=4)) then
write(0,*) 'Going for remapping '
if (.true.) then
associate(lv=>prec%precv(iszv), rmp => prec%precv(iszv)%remap_data)
call lv%desc_ac%clone(rmp%desc_ac_pre_remap,info)
call lv%ac%clone(rmp%ac_pre_remap,info)
if (np >= 8) then
call psb_remap(np/4,rmp%desc_ac_pre_remap,rmp%ac_pre_remap,&
& rmp%idest,rmp%isrc,rmp%nrsrc,rmp%naggr,lv%desc_ac,lv%ac,info)
else
call psb_remap(np/2,rmp%desc_ac_pre_remap,rmp%ac_pre_remap,&
& rmp%idest,rmp%isrc,rmp%nrsrc,rmp%naggr,lv%desc_ac,lv%ac,info)
end if
write(0,*) me,' Out of remapping ',rmp%desc_ac_pre_remap%get_fmt(),' ',&
& lv%desc_ac%get_fmt(),sum(lv%linmap%naggr),sum(rmp%naggr)
lv%linmap%naggr(:) = rmp%naggr(:)
lv%linmap%p_desc_V => rmp%desc_ac_pre_remap
lv%base_a => lv%ac
lv%base_desc => lv%desc_ac
end associate
end if
end if
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Internal hierarchy build' ) & a_err='Internal hierarchy build' )

6
amgprec/impl/amg_z_smoothers_bld.f90
Unescape Escape View File

@ -117,6 +117,10 @@ subroutine amg_z_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
info = psb_success_ info = psb_success_
ctxt = desc_a%get_context() ctxt = desc_a%get_context()
call psb_info(ctxt, me, np) call psb_info(ctxt, me, np)
if (me <0) then
!!$ write(0,*) 'out of CTXT, should not do anything '
goto 9998
end if
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
@ -289,6 +293,7 @@ subroutine amg_z_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
! !
! build the base preconditioner at level i ! build the base preconditioner at level i
! !
!!$ write(0,*) me,' Building at level ',i
call prec%precv(i)%bld(info,amold=amold,vmold=vmold,imold=imold,ilv=i) call prec%precv(i)%bld(info,amold=amold,vmold=vmold,imold=imold,ilv=i)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -304,6 +309,7 @@ subroutine amg_z_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Exiting with',iszv,' levels' & 'Exiting with',iszv,' levels'
9998 continue
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

338
amgprec/impl/amg_zmlprec_aply.f90
Unescape Escape View File

@ -393,38 +393,39 @@ contains
if(debug_level > 1) then if(debug_level > 1) then
write(debug_unit,*) me,' Start inner_ml_aply at level ',level, info write(debug_unit,*) me,' Start inner_ml_aply at level ',level, info
end if end if
if (me >= 0) then
select case(p%precv(level)%parms%ml_cycle)
select case(p%precv(level)%parms%ml_cycle) case(amg_no_ml_)
!
case(amg_no_ml_) ! No preconditioning, should not really get here
! !
! No preconditioning, should not really get here call psb_errpush(psb_err_internal_error_,name,&
! & a_err='amg_no_ml_ in mlprc_aply?')
call psb_errpush(psb_err_internal_error_,name,& goto 9999
& a_err='amg_no_ml_ in mlprc_aply?')
goto 9999
case(amg_add_ml_) case(amg_add_ml_)
call amg_z_inner_add(p, level, trans, work) call amg_z_inner_add(p, level, trans, work)
case(amg_mult_ml_,amg_vcycle_ml_, amg_wcycle_ml_) case(amg_mult_ml_,amg_vcycle_ml_, amg_wcycle_ml_)
call amg_z_inner_mult(p, level, trans, work) call amg_z_inner_mult(p, level, trans, work)
case(amg_kcycle_ml_, amg_kcyclesym_ml_) case(amg_kcycle_ml_, amg_kcyclesym_ml_)
call amg_z_inner_k_cycle(p, level, trans, work) call amg_z_inner_k_cycle(p, level, trans, work)
case default case default
info = psb_err_from_subroutine_ai_ info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='invalid ml_cycle',& call psb_errpush(info,name,a_err='invalid ml_cycle',&
& i_Err=(/p%precv(level)%parms%ml_cycle,izero,izero,izero,izero/)) & i_Err=(/p%precv(level)%parms%ml_cycle,izero,izero,izero,izero/))
goto 9999 goto 9999
end select end select
if(debug_level > 1) then if(debug_level > 1) then
write(debug_unit,*) me,' End inner_ml_aply at level ',level write(debug_unit,*) me,' End inner_ml_aply at level ',level
end if
end if end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -492,28 +493,30 @@ contains
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv) & wv => p%precv(level)%wrk%wv)
if (allocated(p%precv(level)%sm2a)) then if (me >= 0) then
call psb_geaxpby(zone,vx2l,zzero,vy2l,base_desc,info) if (allocated(p%precv(level)%sm2a)) then
call psb_geaxpby(zone,vx2l,zzero,vy2l,base_desc,info)
sweeps = max(p%precv(level)%parms%sweeps_pre,p%precv(level)%parms%sweeps_post) sweeps = max(p%precv(level)%parms%sweeps_pre,p%precv(level)%parms%sweeps_post)
do k=1, sweeps do k=1, sweeps
call p%precv(level)%sm%apply(zone,& call p%precv(level)%sm%apply(zone,&
& vy2l,zzero,vty,& & vy2l,zzero,vty,&
& base_desc, trans,& & base_desc, trans,&
& ione,work,wv,info,init='Z') & ione,work,wv,info,init='Z')
call p%precv(level)%sm2a%apply(zone,& call p%precv(level)%sm2a%apply(zone,&
& vty,zzero,vy2l,& & vty,zzero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& ione,work,wv,info,init='Z') & ione,work,wv,info,init='Z')
end do end do
else else
sweeps = p%precv(level)%parms%sweeps_pre sweeps = p%precv(level)%parms%sweeps_pre
call p%precv(level)%sm%apply(zone,& call p%precv(level)%sm%apply(zone,&
& vx2l,zzero,vy2l,& & vx2l,zzero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,wv,info,init='Z') & sweeps,work,wv,info,init='Z')
end if
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -602,7 +605,6 @@ contains
end if end if
ctxt = p%precv(level)%base_desc%get_context() ctxt = p%precv(level)%base_desc%get_context()
call psb_info(ctxt, me, np) call psb_info(ctxt, me, np)
if(debug_level > 1) then if(debug_level > 1) then
write(debug_unit,*) me,' inner_mult at level ',level write(debug_unit,*) me,' inner_mult at level ',level
end if end if
@ -623,22 +625,25 @@ contains
! !
if (pre) then if (pre) then
if (trans == 'N') then if (me >=0) then
sweeps = p%precv(level)%parms%sweeps_pre !!$ write(0,*) me,'Applying smoother pre ', level
if (info == psb_success_) call p%precv(level)%sm%apply(zone,& if (trans == 'N') then
& vx2l,zzero,vy2l,base_desc, trans,& sweeps = p%precv(level)%parms%sweeps_pre
& sweeps,work,wv,info,init='Z') if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
else & vx2l,zzero,vy2l,base_desc, trans,&
sweeps = p%precv(level)%parms%sweeps_post & sweeps,work,wv,info,init='Z')
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,& else
& vx2l,zzero,vy2l, base_desc, trans,& sweeps = p%precv(level)%parms%sweeps_post
& sweeps,work,wv,info,init='Z') if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
end if & vx2l,zzero,vy2l, base_desc, trans,&
& sweeps,work,wv,info,init='Z')
if (info /= psb_success_) then end if
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during PRE smoother_apply') if (info /= psb_success_) then
goto 9999 call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during PRE smoother_apply')
goto 9999
end if
end if end if
endif endif
! !
@ -696,11 +701,13 @@ contains
if (p%precv(level)%parms%ml_cycle == amg_wcycle_ml_) then if (p%precv(level)%parms%ml_cycle == amg_wcycle_ml_) then
call psb_geaxpby(zone,vx2l, zzero,vty,& if (me >=0) then
& base_desc,info) call psb_geaxpby(zone,vx2l, zzero,vty,&
if (info == psb_success_) call psb_spmm(-zone,base_a,& & base_desc,info)
& vy2l,zone,vty,& if (info == psb_success_) call psb_spmm(-zone,base_a,&
& base_desc,info,work=work,trans=trans) & vy2l,zone,vty,&
& base_desc,info,work=work,trans=trans)
end if
if (info == psb_success_) & if (info == psb_success_) &
& call p%precv(level+1)%map_rstr(zone,vty,& & call p%precv(level+1)%map_rstr(zone,vty,&
& zzero,p%precv(level+1)%wrk%vx2l,info,work=work,& & zzero,p%precv(level+1)%wrk%vx2l,info,work=work,&
@ -728,31 +735,33 @@ contains
if (post) then if (post) then
call psb_geaxpby(zone,vx2l,& if (me >=0) then
& zzero,vty,& call psb_geaxpby(zone,vx2l,&
& base_desc,info) & zzero,vty,&
if (info == psb_success_) call psb_spmm(-zone,base_a,& & base_desc,info)
& vy2l, zone,vty,base_desc,info,& if (info == psb_success_) call psb_spmm(-zone,base_a,&
& work=work,trans=trans) & vy2l, zone,vty,base_desc,info,&
if (info /= psb_success_) then & work=work,trans=trans)
call psb_errpush(psb_err_internal_error_,name,& if (info /= psb_success_) then
& a_err='Error during residue') call psb_errpush(psb_err_internal_error_,name,&
goto 9999 & a_err='Error during residue')
end if goto 9999
end if
!
! Apply the second smoother !
! ! Apply the second smoother
if (trans == 'N') then !
sweeps = p%precv(level)%parms%sweeps_post if (trans == 'N') then
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,& sweeps = p%precv(level)%parms%sweeps_post
& vty,zone,vy2l, base_desc, trans,& if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
& sweeps,work,wv,info,init='Z') & vty,zone,vy2l, base_desc, trans,&
else & sweeps,work,wv,info,init='Z')
sweeps = p%precv(level)%parms%sweeps_pre else
if (info == psb_success_) call p%precv(level)%sm%apply(zone,& sweeps = p%precv(level)%parms%sweeps_pre
& vty,zone,vy2l, base_desc, trans,& if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& sweeps,work,wv,info,init='Z') & vty,zone,vy2l, base_desc, trans,&
& sweeps,work,wv,info,init='Z')
end if
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -764,11 +773,14 @@ contains
endif endif
else if (level == nlev) then else if (level == nlev) then
!!$ write(0,*) me,'Applying smoother at top level ',psb_errstatus_fatal()
sweeps = p%precv(level)%parms%sweeps_pre if (me >=0) then
if (info == psb_success_) call p%precv(level)%sm%apply(zone,& sweeps = p%precv(level)%parms%sweeps_pre
& vx2l,zzero,vy2l,base_desc, trans,& if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& sweeps,work,wv,info) & vx2l,zzero,vy2l,base_desc, trans,&
& sweeps,work,wv,info)
end if
!!$ write(0,*) me,' Done applying smoother at top level ',psb_errstatus_fatal()
else else
@ -778,7 +790,7 @@ contains
goto 9999 goto 9999
end if end if
end associate end associate
9998 continue
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -847,52 +859,53 @@ contains
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv(8:)) & wv => p%precv(level)%wrk%wv(8:))
if (level == nlev) then
!
! Apply smoother
!
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& vx2l,zzero,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
else if (level < nlev) then if (level == nlev) then
if (me >= 0) then
if (trans == 'N') then !
! Apply smoother
!
sweeps = p%precv(level)%parms%sweeps_pre sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,& if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& vx2l,zzero,vy2l,base_desc, trans,& & vx2l,zzero,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z') & sweeps,work,wv,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
& vx2l,zzero,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-PRE smoother_apply')
goto 9999
end if end if
else if (level < nlev) then
if (me >= 0) then
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& vx2l,zzero,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
& vx2l,zzero,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-PRE smoother_apply')
goto 9999
end if
! !
! Compute the residual and call recursively ! Compute the residual and call recursively
! !
call psb_geaxpby(zone,vx2l,& call psb_geaxpby(zone,vx2l,&
& zzero,vty,& & zzero,vty,&
& base_desc,info) & base_desc,info)
if (info == psb_success_) call psb_spmm(-zone,base_a,& if (info == psb_success_) call psb_spmm(-zone,base_a,&
& vy2l,zone,vty,base_desc,info,work=work,trans=trans) & vy2l,zone,vty,base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue') & a_err='Error during residue')
goto 9999 goto 9999
end if
end if end if
! Apply the restriction ! Apply the restriction
call p%precv(level + 1)%map_rstr(zone,vty,& call p%precv(level + 1)%map_rstr(zone,vty,&
& zzero,p%precv(level + 1)%wrk%vx2l,& & zzero,p%precv(level + 1)%wrk%vx2l,&
@ -940,40 +953,42 @@ contains
& a_err='Error during prolongation') & a_err='Error during prolongation')
goto 9999 goto 9999
end if end if
if (me >= 0) then
!
! Compute the residual
!
call psb_geaxpby(zone,vx2l,&
& zzero,vty,base_desc,info)
call psb_spmm(-zone,base_a,vy2l,&
& zone,vty,base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!
! Apply the smoother
!
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
& vty,zone,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& vty,zone,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
end if
! if (info /= psb_success_) then
! Compute the residual call psb_errpush(psb_err_internal_error_,name,&
! & a_err='Error during POST smoother_apply')
call psb_geaxpby(zone,vx2l,& goto 9999
& zzero,vty,base_desc,info) end if
call psb_spmm(-zone,base_a,vy2l,&
& zone,vty,base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!
! Apply the smoother
!
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
& vty,zone,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& vty,zone,vy2l,base_desc, trans,&
& sweeps,work,wv,info,init='Z')
end if end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during POST smoother_apply')
goto 9999
end if
else else
info = psb_err_internal_error_ info = psb_err_internal_error_
@ -991,7 +1006,6 @@ contains
end subroutine amg_z_inner_k_cycle end subroutine amg_z_inner_k_cycle
recursive subroutine amg_zinneritkcycle(p, level, trans, work, innersolv) recursive subroutine amg_zinneritkcycle(p, level, trans, work, innersolv)
use psb_base_mod use psb_base_mod
use amg_prec_mod use amg_prec_mod

129
amgprec/impl/level/amg_c_base_onelev_build.f90
Unescape Escape View File

@ -71,72 +71,81 @@ subroutine amg_c_base_onelev_build(lv,info,amold,vmold,imold,ilv)
ctxt = lv%base_desc%get_ctxt() ctxt = lv%base_desc%get_ctxt()
call psb_info(ctxt,me,np) call psb_info(ctxt,me,np)
if (.not.allocated(lv%sm)) then !
!! Error: should have called amg_dprecinit ! At top level(s) I may be using
info=3111 ! a context with less processes
call psb_errpush(info,name) !
goto 9999 if (me < 0) then
end if !!$ write(0,*) 'onelevbld: I am excluded from this one '
if (.not.allocated(lv%sm%sv)) then else
!! Error: should have called amg_dprecinit !!$ write(0,*) me,' Going to build smoothers at this level '
info=3111 if (.not.allocated(lv%sm)) then
call psb_errpush(info,name) !! Error: should have called amg_dprecinit
goto 9999 info=3111
end if call psb_errpush(info,name)
lv%ac_nz_loc = lv%ac%get_nzeros() goto 9999
lv%ac_nz_tot = lv%ac_nz_loc end if
select case(lv%parms%coarse_mat) if (.not.allocated(lv%sm%sv)) then
case(amg_distr_mat_) !! Error: should have called amg_dprecinit
call psb_sum(ctxt,lv%ac_nz_tot) info=3111
case(amg_repl_mat_) call psb_errpush(info,name)
! Do nothing goto 9999
case default end if
! Should never get here lv%ac_nz_loc = lv%ac%get_nzeros()
info = psb_err_internal_error_ lv%ac_nz_tot = lv%ac_nz_loc
call psb_errpush(info,name,a_err='Wrong lv%parms') select case(lv%parms%coarse_mat)
goto 9999 case(amg_distr_mat_)
end select call psb_sum(ctxt,lv%ac_nz_tot)
case(amg_repl_mat_)
! Do nothing
case default
! Should never get here
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='Wrong lv%parms')
goto 9999
end select
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Calling mlprcbld at level ',i & 'Calling mlprcbld at level ',i
call amg_check_def(lv%parms%sweeps_pre,& call amg_check_def(lv%parms%sweeps_pre,&
& 'Jacobi sweeps',izero,is_int_non_negative) & 'Jacobi sweeps',izero,is_int_non_negative)
call amg_check_def(lv%parms%sweeps_post,& call amg_check_def(lv%parms%sweeps_post,&
& 'Jacobi sweeps',izero,is_int_non_negative) & 'Jacobi sweeps',izero,is_int_non_negative)
call lv%sm%build(lv%base_a,lv%base_desc,info) call lv%sm%build(lv%base_a,lv%base_desc,info)
if (info == 0) then if (info == 0) then
if (allocated(lv%sm2a)) then if (allocated(lv%sm2a)) then
call lv%sm2a%build(lv%base_a,lv%base_desc,info) call lv%sm2a%build(lv%base_a,lv%base_desc,info)
lv%sm2 => lv%sm2a lv%sm2 => lv%sm2a
else else
lv%sm2 => lv%sm lv%sm2 => lv%sm
end if
end if
if (info /=0 ) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Smoother bld error')
goto 9999
end if end if
end if
if (info /=0 ) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Smoother bld error')
goto 9999
end if
if (lv%sm%sv%is_global()) then if (lv%sm%sv%is_global()) then
if ((lv%parms%sweeps_pre>1).or.(lv%parms%sweeps_post>1)) then if ((lv%parms%sweeps_pre>1).or.(lv%parms%sweeps_post>1)) then
lv%parms%sweeps_pre = 1 lv%parms%sweeps_pre = 1
lv%parms%sweeps_post = 1 lv%parms%sweeps_post = 1
if (me == 0) then if (me == 0) then
write(debug_unit,*) write(debug_unit,*)
if (present(ilv)) then if (present(ilv)) then
write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),& write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),&
& '" at level ',ilv & '" at level ',ilv
write(debug_unit,*) ' is configured as a global solver ' write(debug_unit,*) ' is configured as a global solver '
else else
write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),& write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),&
& '" is configured as a global solver ' & '" is configured as a global solver '
end if
write(debug_unit,*) ' Pre and post sweeps at this level reset to 1'
end if end if
write(debug_unit,*) ' Pre and post sweeps at this level reset to 1'
end if end if
end if end if
end if end if

52
amgprec/impl/level/amg_c_base_onelev_map_prol.F90
Unescape Escape View File

@ -47,11 +47,61 @@ subroutine amg_c_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt
complex(psb_spk_), optional :: work(:) complex(psb_spk_), optional :: work(:)
type(psb_c_vect_type), optional, target, intent(inout) :: vtx,vty type(psb_c_vect_type), optional, target, intent(inout) :: vtx,vty
!!$ write(0,*) 'New map_prol',lv%remap_data%ac_pre_remap%is_asb()
if (lv%remap_data%ac_pre_remap%is_asb()) then if (lv%remap_data%ac_pre_remap%is_asb()) then
! !
! Remap has happened, deal with it ! Remap has happened, deal with it
! !
write(0,*) 'Remap handling not implemented yet ' !!$ write(0,*) 'Remap handling '
block
type(psb_ctxt_type) :: ctxt, nctxt
integer(psb_ipk_) :: i,j,ip,idest, nsrc, nrl, nrc, kp
integer(psb_ipk_) :: me, np, rme, rnp
complex(psb_spk_), allocatable :: rsnd(:), rrcv(:)
type(psb_c_vect_type) :: tv
ctxt = lv%remap_data%desc_ac_pre_remap%get_ctxt()
call psb_info(ctxt,me,np)
!!$ write(0,*) 'Old context ',me,np,psb_errstatus_fatal()
nctxt = lv%desc_ac%get_ctxt()
call psb_info(nctxt,rme,rnp)
!!$ write(0,*) 'New context ',rme,rnp,psb_errstatus_fatal()
idest = lv%remap_data%idest
associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc)
!!$ write(0,*) 'Should apply maps, then receive data from ',idest,' to ',me,psb_errstatus_fatal()
nsrc = size(isrc)
nrl = lv%remap_data%desc_ac_pre_remap%get_local_rows()
nrc = lv%remap_data%desc_ac_pre_remap%get_local_cols()
if (rme >=0) then
allocate(rrcv(sum(nrsrc)))
rrcv = vect_v%get_vect()
!!$ write(0,*) me,rme,' Size check ',size(rrcv),lv%desc_ac%get_local_rows(),psb_errstatus_fatal()
kp = 0
do i = 1,size(isrc)
ip = isrc(i)
nrl = nrsrc(i)
!!$ write(0,*) me,' Sending to ',ip,nrl,kp+1,kp+nrl
call psb_snd(ctxt,rrcv(kp+1:kp+nrl),ip)
kp = kp + nrl
end do
end if
nrl = lv%remap_data%desc_ac_pre_remap%get_local_rows()
call psb_geall(tv,lv%remap_data%desc_ac_pre_remap,info)
!!$ write(0,*) me, ' Allocated ',nrl,info,psb_errstatus_fatal()
call psb_geasb(tv,lv%remap_data%desc_ac_pre_remap,info)
!!$ write(0,*) me,' Size of TV ',nrl,tv%get_nrows(),info
!!$ write(0,*) me,' Receiving from ',idest,nrl,psb_errstatus_fatal()
call psb_realloc(nrc,rsnd,info)
call psb_rcv(ctxt,rsnd(1:nrl),idest)
call tv%set_vect(rsnd)
call lv%linmap%map_V2U(alpha,tv,beta,vect_u,info,&
& work=work,vtx=vtx,vty=vty)
end associate
!!$ write(0,*) me, ' Prolongator with remap done '
!!$ flush(0)
!!$ call psb_barrier(ctxt)
end block
else else
! Default transfer ! Default transfer
call lv%linmap%map_V2U(alpha,vect_v,beta,vect_u,info,& call lv%linmap%map_V2U(alpha,vect_v,beta,vect_u,info,&

43
amgprec/impl/level/amg_c_base_onelev_map_rstr.F90
Unescape Escape View File

@ -36,7 +36,8 @@
! !
! !
subroutine amg_c_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,work,vtx,vty) subroutine amg_c_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
& work,vtx,vty)
use psb_base_mod use psb_base_mod
use amg_c_onelev_mod, amg_protect_name => amg_c_base_onelev_map_rstr_v use amg_c_onelev_mod, amg_protect_name => amg_c_base_onelev_map_rstr_v
implicit none implicit none
@ -47,24 +48,52 @@ subroutine amg_c_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,work,vt
complex(psb_spk_), optional :: work(:) complex(psb_spk_), optional :: work(:)
type(psb_c_vect_type), optional, target, intent(inout) :: vtx,vty type(psb_c_vect_type), optional, target, intent(inout) :: vtx,vty
!!$ write(0,*) 'New map_rstr',lv%remap_data%ac_pre_remap%is_asb()
if (lv%remap_data%ac_pre_remap%is_asb()) then if (lv%remap_data%ac_pre_remap%is_asb()) then
! !
! Remap has happened, deal with it ! Remap has happened, deal with it
! !
write(0,*) 'Remap handling not implemented yet ' !!$ write(0,*) 'Remap handling not implemented yet '
block block
integer(psb_ipk_) :: i,j,ip,nctxt,ictxt, idest type(psb_ctxt_type) :: ctxt, nctxt
integer(psb_ipk_) :: i,j,ip, idest, nsrc, nrl, kp
integer(psb_ipk_) :: me, np, rme, rnp integer(psb_ipk_) :: me, np, rme, rnp
complex(psb_spk_), allocatable :: rsnd(:), rrcv(:)
type(psb_c_vect_type) :: tv
ictxt = lv%remap_data%desc_ac_pre_remap%get_ctxt() ctxt = lv%remap_data%desc_ac_pre_remap%get_ctxt()
call psb_info(ictxt,me,np) call psb_info(ctxt,me,np)
nctxt = lv%desc_ac%get_ctxt() nctxt = lv%desc_ac%get_ctxt()
call psb_info(nctxt,rme,rnp) call psb_info(nctxt,rme,rnp)
!!$ write(0,*) 'New context ',rme,rnp
idest = lv%remap_data%idest idest = lv%remap_data%idest
associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc) associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc)
write(0,*) 'Should apply maps, then send data from ',me,' to ',idest !!$ write(0,*) 'Should apply maps, then send data from ',me,' to ',idest
if (rme >= 0) write(0,*) rme, ' Receiving data from ',isrc(:) !!$ if (rme >= 0) write(0,*) rme, ' Receiving data from ',isrc(:)
nsrc = size(isrc)
nrl = lv%remap_data%desc_ac_pre_remap%get_local_rows()
call psb_geall(tv,lv%remap_data%desc_ac_pre_remap,info)
call psb_geasb(tv,lv%remap_data%desc_ac_pre_remap,info)
!!$ write(0,*) me,' Size of TV ',tv%get_nrows()
call lv%linmap%map_U2V(alpha,vect_u,beta,tv,info,&
& work=work,vtx=vtx,vty=vty)
rsnd = tv%get_vect()
call psb_snd(ctxt,rsnd(1:nrl),idest)
if (rme >=0) then
allocate(rrcv(sum(nrsrc)))
!!$ write(0,*) me,rme,' Size check ',size(rrcv)!,lv%desc_ac%get_local_rows()
kp = 0
do i = 1,size(isrc)
ip = isrc(i)
nrl = nrsrc(i)
!!$ write(0,*) me,' Receiving from ',ip,nrl,kp+1,kp+nrl,size(rrcv)
call psb_rcv(ctxt,rrcv(kp+1:kp+nrl),ip)
kp = kp + nrl
end do
call vect_v%set_vect(rrcv)
end if
end associate end associate
!!$ write(0,*) me, ' Restrictor with remap done '
end block end block
else else

129
amgprec/impl/level/amg_d_base_onelev_build.f90
Unescape Escape View File

@ -71,72 +71,81 @@ subroutine amg_d_base_onelev_build(lv,info,amold,vmold,imold,ilv)
ctxt = lv%base_desc%get_ctxt() ctxt = lv%base_desc%get_ctxt()
call psb_info(ctxt,me,np) call psb_info(ctxt,me,np)
if (.not.allocated(lv%sm)) then !
!! Error: should have called amg_dprecinit ! At top level(s) I may be using
info=3111 ! a context with less processes
call psb_errpush(info,name) !
goto 9999 if (me < 0) then
end if !!$ write(0,*) 'onelevbld: I am excluded from this one '
if (.not.allocated(lv%sm%sv)) then else
!! Error: should have called amg_dprecinit !!$ write(0,*) me,' Going to build smoothers at this level '
info=3111 if (.not.allocated(lv%sm)) then
call psb_errpush(info,name) !! Error: should have called amg_dprecinit
goto 9999 info=3111
end if call psb_errpush(info,name)
lv%ac_nz_loc = lv%ac%get_nzeros() goto 9999
lv%ac_nz_tot = lv%ac_nz_loc end if
select case(lv%parms%coarse_mat) if (.not.allocated(lv%sm%sv)) then
case(amg_distr_mat_) !! Error: should have called amg_dprecinit
call psb_sum(ctxt,lv%ac_nz_tot) info=3111
case(amg_repl_mat_) call psb_errpush(info,name)
! Do nothing goto 9999
case default end if
! Should never get here lv%ac_nz_loc = lv%ac%get_nzeros()
info = psb_err_internal_error_ lv%ac_nz_tot = lv%ac_nz_loc
call psb_errpush(info,name,a_err='Wrong lv%parms') select case(lv%parms%coarse_mat)
goto 9999 case(amg_distr_mat_)
end select call psb_sum(ctxt,lv%ac_nz_tot)
case(amg_repl_mat_)
! Do nothing
case default
! Should never get here
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='Wrong lv%parms')
goto 9999
end select
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Calling mlprcbld at level ',i & 'Calling mlprcbld at level ',i
call amg_check_def(lv%parms%sweeps_pre,& call amg_check_def(lv%parms%sweeps_pre,&
& 'Jacobi sweeps',izero,is_int_non_negative) & 'Jacobi sweeps',izero,is_int_non_negative)
call amg_check_def(lv%parms%sweeps_post,& call amg_check_def(lv%parms%sweeps_post,&
& 'Jacobi sweeps',izero,is_int_non_negative) & 'Jacobi sweeps',izero,is_int_non_negative)
call lv%sm%build(lv%base_a,lv%base_desc,info) call lv%sm%build(lv%base_a,lv%base_desc,info)
if (info == 0) then if (info == 0) then
if (allocated(lv%sm2a)) then if (allocated(lv%sm2a)) then
call lv%sm2a%build(lv%base_a,lv%base_desc,info) call lv%sm2a%build(lv%base_a,lv%base_desc,info)
lv%sm2 => lv%sm2a lv%sm2 => lv%sm2a
else else
lv%sm2 => lv%sm lv%sm2 => lv%sm
end if
end if
if (info /=0 ) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Smoother bld error')
goto 9999
end if end if
end if
if (info /=0 ) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Smoother bld error')
goto 9999
end if
if (lv%sm%sv%is_global()) then if (lv%sm%sv%is_global()) then
if ((lv%parms%sweeps_pre>1).or.(lv%parms%sweeps_post>1)) then if ((lv%parms%sweeps_pre>1).or.(lv%parms%sweeps_post>1)) then
lv%parms%sweeps_pre = 1 lv%parms%sweeps_pre = 1
lv%parms%sweeps_post = 1 lv%parms%sweeps_post = 1
if (me == 0) then if (me == 0) then
write(debug_unit,*) write(debug_unit,*)
if (present(ilv)) then if (present(ilv)) then
write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),& write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),&
& '" at level ',ilv & '" at level ',ilv
write(debug_unit,*) ' is configured as a global solver ' write(debug_unit,*) ' is configured as a global solver '
else else
write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),& write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),&
& '" is configured as a global solver ' & '" is configured as a global solver '
end if
write(debug_unit,*) ' Pre and post sweeps at this level reset to 1'
end if end if
write(debug_unit,*) ' Pre and post sweeps at this level reset to 1'
end if end if
end if end if
end if end if

52
amgprec/impl/level/amg_d_base_onelev_map_prol.F90
Unescape Escape View File

@ -47,11 +47,61 @@ subroutine amg_d_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt
real(psb_dpk_), optional :: work(:) real(psb_dpk_), optional :: work(:)
type(psb_d_vect_type), optional, target, intent(inout) :: vtx,vty type(psb_d_vect_type), optional, target, intent(inout) :: vtx,vty
!!$ write(0,*) 'New map_prol',lv%remap_data%ac_pre_remap%is_asb()
if (lv%remap_data%ac_pre_remap%is_asb()) then if (lv%remap_data%ac_pre_remap%is_asb()) then
! !
! Remap has happened, deal with it ! Remap has happened, deal with it
! !
write(0,*) 'Remap handling not implemented yet ' !!$ write(0,*) 'Remap handling '
block
type(psb_ctxt_type) :: ctxt, nctxt
integer(psb_ipk_) :: i,j,ip,idest, nsrc, nrl, nrc, kp
integer(psb_ipk_) :: me, np, rme, rnp
real(psb_dpk_), allocatable :: rsnd(:), rrcv(:)
type(psb_d_vect_type) :: tv
ctxt = lv%remap_data%desc_ac_pre_remap%get_ctxt()
call psb_info(ctxt,me,np)
!!$ write(0,*) 'Old context ',me,np,psb_errstatus_fatal()
nctxt = lv%desc_ac%get_ctxt()
call psb_info(nctxt,rme,rnp)
!!$ write(0,*) 'New context ',rme,rnp,psb_errstatus_fatal()
idest = lv%remap_data%idest
associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc)
!!$ write(0,*) 'Should apply maps, then receive data from ',idest,' to ',me,psb_errstatus_fatal()
nsrc = size(isrc)
nrl = lv%remap_data%desc_ac_pre_remap%get_local_rows()
nrc = lv%remap_data%desc_ac_pre_remap%get_local_cols()
if (rme >=0) then
allocate(rrcv(sum(nrsrc)))
rrcv = vect_v%get_vect()
!!$ write(0,*) me,rme,' Size check ',size(rrcv),lv%desc_ac%get_local_rows(),psb_errstatus_fatal()
kp = 0
do i = 1,size(isrc)
ip = isrc(i)
nrl = nrsrc(i)
!!$ write(0,*) me,' Sending to ',ip,nrl,kp+1,kp+nrl
call psb_snd(ctxt,rrcv(kp+1:kp+nrl),ip)
kp = kp + nrl
end do
end if
nrl = lv%remap_data%desc_ac_pre_remap%get_local_rows()
call psb_geall(tv,lv%remap_data%desc_ac_pre_remap,info)
!!$ write(0,*) me, ' Allocated ',nrl,info,psb_errstatus_fatal()
call psb_geasb(tv,lv%remap_data%desc_ac_pre_remap,info)
!!$ write(0,*) me,' Size of TV ',nrl,tv%get_nrows(),info
!!$ write(0,*) me,' Receiving from ',idest,nrl,psb_errstatus_fatal()
call psb_realloc(nrc,rsnd,info)
call psb_rcv(ctxt,rsnd(1:nrl),idest)
call tv%set_vect(rsnd)
call lv%linmap%map_V2U(alpha,tv,beta,vect_u,info,&
& work=work,vtx=vtx,vty=vty)
end associate
!!$ write(0,*) me, ' Prolongator with remap done '
!!$ flush(0)
!!$ call psb_barrier(ctxt)
end block
else else
! Default transfer ! Default transfer
call lv%linmap%map_V2U(alpha,vect_v,beta,vect_u,info,& call lv%linmap%map_V2U(alpha,vect_v,beta,vect_u,info,&

43
amgprec/impl/level/amg_d_base_onelev_map_rstr.F90
Unescape Escape View File

@ -36,7 +36,8 @@
! !
! !
subroutine amg_d_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,work,vtx,vty) subroutine amg_d_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
& work,vtx,vty)
use psb_base_mod use psb_base_mod
use amg_d_onelev_mod, amg_protect_name => amg_d_base_onelev_map_rstr_v use amg_d_onelev_mod, amg_protect_name => amg_d_base_onelev_map_rstr_v
implicit none implicit none
@ -47,24 +48,52 @@ subroutine amg_d_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,work,vt
real(psb_dpk_), optional :: work(:) real(psb_dpk_), optional :: work(:)
type(psb_d_vect_type), optional, target, intent(inout) :: vtx,vty type(psb_d_vect_type), optional, target, intent(inout) :: vtx,vty
!!$ write(0,*) 'New map_rstr',lv%remap_data%ac_pre_remap%is_asb()
if (lv%remap_data%ac_pre_remap%is_asb()) then if (lv%remap_data%ac_pre_remap%is_asb()) then
! !
! Remap has happened, deal with it ! Remap has happened, deal with it
! !
write(0,*) 'Remap handling not implemented yet ' !!$ write(0,*) 'Remap handling not implemented yet '
block block
integer(psb_ipk_) :: i,j,ip,nctxt,ictxt, idest type(psb_ctxt_type) :: ctxt, nctxt
integer(psb_ipk_) :: i,j,ip, idest, nsrc, nrl, kp
integer(psb_ipk_) :: me, np, rme, rnp integer(psb_ipk_) :: me, np, rme, rnp
real(psb_dpk_), allocatable :: rsnd(:), rrcv(:)
type(psb_d_vect_type) :: tv
ictxt = lv%remap_data%desc_ac_pre_remap%get_ctxt() ctxt = lv%remap_data%desc_ac_pre_remap%get_ctxt()
call psb_info(ictxt,me,np) call psb_info(ctxt,me,np)
nctxt = lv%desc_ac%get_ctxt() nctxt = lv%desc_ac%get_ctxt()
call psb_info(nctxt,rme,rnp) call psb_info(nctxt,rme,rnp)
!!$ write(0,*) 'New context ',rme,rnp
idest = lv%remap_data%idest idest = lv%remap_data%idest
associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc) associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc)
write(0,*) 'Should apply maps, then send data from ',me,' to ',idest !!$ write(0,*) 'Should apply maps, then send data from ',me,' to ',idest
if (rme >= 0) write(0,*) rme, ' Receiving data from ',isrc(:) !!$ if (rme >= 0) write(0,*) rme, ' Receiving data from ',isrc(:)
nsrc = size(isrc)
nrl = lv%remap_data%desc_ac_pre_remap%get_local_rows()
call psb_geall(tv,lv%remap_data%desc_ac_pre_remap,info)
call psb_geasb(tv,lv%remap_data%desc_ac_pre_remap,info)
!!$ write(0,*) me,' Size of TV ',tv%get_nrows()
call lv%linmap%map_U2V(alpha,vect_u,beta,tv,info,&
& work=work,vtx=vtx,vty=vty)
rsnd = tv%get_vect()
call psb_snd(ctxt,rsnd(1:nrl),idest)
if (rme >=0) then
allocate(rrcv(sum(nrsrc)))
!!$ write(0,*) me,rme,' Size check ',size(rrcv)!,lv%desc_ac%get_local_rows()
kp = 0
do i = 1,size(isrc)
ip = isrc(i)
nrl = nrsrc(i)
!!$ write(0,*) me,' Receiving from ',ip,nrl,kp+1,kp+nrl,size(rrcv)
call psb_rcv(ctxt,rrcv(kp+1:kp+nrl),ip)
kp = kp + nrl
end do
call vect_v%set_vect(rrcv)
end if
end associate end associate
!!$ write(0,*) me, ' Restrictor with remap done '
end block end block
else else

129
amgprec/impl/level/amg_s_base_onelev_build.f90
Unescape Escape View File

@ -71,72 +71,81 @@ subroutine amg_s_base_onelev_build(lv,info,amold,vmold,imold,ilv)
ctxt = lv%base_desc%get_ctxt() ctxt = lv%base_desc%get_ctxt()
call psb_info(ctxt,me,np) call psb_info(ctxt,me,np)
if (.not.allocated(lv%sm)) then !
!! Error: should have called amg_dprecinit ! At top level(s) I may be using
info=3111 ! a context with less processes
call psb_errpush(info,name) !
goto 9999 if (me < 0) then
end if !!$ write(0,*) 'onelevbld: I am excluded from this one '
if (.not.allocated(lv%sm%sv)) then else
!! Error: should have called amg_dprecinit !!$ write(0,*) me,' Going to build smoothers at this level '
info=3111 if (.not.allocated(lv%sm)) then
call psb_errpush(info,name) !! Error: should have called amg_dprecinit
goto 9999 info=3111
end if call psb_errpush(info,name)
lv%ac_nz_loc = lv%ac%get_nzeros() goto 9999
lv%ac_nz_tot = lv%ac_nz_loc end if
select case(lv%parms%coarse_mat) if (.not.allocated(lv%sm%sv)) then
case(amg_distr_mat_) !! Error: should have called amg_dprecinit
call psb_sum(ctxt,lv%ac_nz_tot) info=3111
case(amg_repl_mat_) call psb_errpush(info,name)
! Do nothing goto 9999
case default end if
! Should never get here lv%ac_nz_loc = lv%ac%get_nzeros()
info = psb_err_internal_error_ lv%ac_nz_tot = lv%ac_nz_loc
call psb_errpush(info,name,a_err='Wrong lv%parms') select case(lv%parms%coarse_mat)
goto 9999 case(amg_distr_mat_)
end select call psb_sum(ctxt,lv%ac_nz_tot)
case(amg_repl_mat_)
! Do nothing
case default
! Should never get here
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='Wrong lv%parms')
goto 9999
end select
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Calling mlprcbld at level ',i & 'Calling mlprcbld at level ',i
call amg_check_def(lv%parms%sweeps_pre,& call amg_check_def(lv%parms%sweeps_pre,&
& 'Jacobi sweeps',izero,is_int_non_negative) & 'Jacobi sweeps',izero,is_int_non_negative)
call amg_check_def(lv%parms%sweeps_post,& call amg_check_def(lv%parms%sweeps_post,&
& 'Jacobi sweeps',izero,is_int_non_negative) & 'Jacobi sweeps',izero,is_int_non_negative)
call lv%sm%build(lv%base_a,lv%base_desc,info) call lv%sm%build(lv%base_a,lv%base_desc,info)
if (info == 0) then if (info == 0) then
if (allocated(lv%sm2a)) then if (allocated(lv%sm2a)) then
call lv%sm2a%build(lv%base_a,lv%base_desc,info) call lv%sm2a%build(lv%base_a,lv%base_desc,info)
lv%sm2 => lv%sm2a lv%sm2 => lv%sm2a
else else
lv%sm2 => lv%sm lv%sm2 => lv%sm
end if
end if
if (info /=0 ) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Smoother bld error')
goto 9999
end if end if
end if
if (info /=0 ) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Smoother bld error')
goto 9999
end if
if (lv%sm%sv%is_global()) then if (lv%sm%sv%is_global()) then
if ((lv%parms%sweeps_pre>1).or.(lv%parms%sweeps_post>1)) then if ((lv%parms%sweeps_pre>1).or.(lv%parms%sweeps_post>1)) then
lv%parms%sweeps_pre = 1 lv%parms%sweeps_pre = 1
lv%parms%sweeps_post = 1 lv%parms%sweeps_post = 1
if (me == 0) then if (me == 0) then
write(debug_unit,*) write(debug_unit,*)
if (present(ilv)) then if (present(ilv)) then
write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),& write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),&
& '" at level ',ilv & '" at level ',ilv
write(debug_unit,*) ' is configured as a global solver ' write(debug_unit,*) ' is configured as a global solver '
else else
write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),& write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),&
& '" is configured as a global solver ' & '" is configured as a global solver '
end if
write(debug_unit,*) ' Pre and post sweeps at this level reset to 1'
end if end if
write(debug_unit,*) ' Pre and post sweeps at this level reset to 1'
end if end if
end if end if
end if end if

52
amgprec/impl/level/amg_s_base_onelev_map_prol.F90
Unescape Escape View File

@ -47,11 +47,61 @@ subroutine amg_s_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt
real(psb_spk_), optional :: work(:) real(psb_spk_), optional :: work(:)
type(psb_s_vect_type), optional, target, intent(inout) :: vtx,vty type(psb_s_vect_type), optional, target, intent(inout) :: vtx,vty
!!$ write(0,*) 'New map_prol',lv%remap_data%ac_pre_remap%is_asb()
if (lv%remap_data%ac_pre_remap%is_asb()) then if (lv%remap_data%ac_pre_remap%is_asb()) then
! !
! Remap has happened, deal with it ! Remap has happened, deal with it
! !
write(0,*) 'Remap handling not implemented yet ' !!$ write(0,*) 'Remap handling '
block
type(psb_ctxt_type) :: ctxt, nctxt
integer(psb_ipk_) :: i,j,ip,idest, nsrc, nrl, nrc, kp
integer(psb_ipk_) :: me, np, rme, rnp
real(psb_spk_), allocatable :: rsnd(:), rrcv(:)
type(psb_s_vect_type) :: tv
ctxt = lv%remap_data%desc_ac_pre_remap%get_ctxt()
call psb_info(ctxt,me,np)
!!$ write(0,*) 'Old context ',me,np,psb_errstatus_fatal()
nctxt = lv%desc_ac%get_ctxt()
call psb_info(nctxt,rme,rnp)
!!$ write(0,*) 'New context ',rme,rnp,psb_errstatus_fatal()
idest = lv%remap_data%idest
associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc)
!!$ write(0,*) 'Should apply maps, then receive data from ',idest,' to ',me,psb_errstatus_fatal()
nsrc = size(isrc)
nrl = lv%remap_data%desc_ac_pre_remap%get_local_rows()
nrc = lv%remap_data%desc_ac_pre_remap%get_local_cols()
if (rme >=0) then
allocate(rrcv(sum(nrsrc)))
rrcv = vect_v%get_vect()
!!$ write(0,*) me,rme,' Size check ',size(rrcv),lv%desc_ac%get_local_rows(),psb_errstatus_fatal()
kp = 0
do i = 1,size(isrc)
ip = isrc(i)
nrl = nrsrc(i)
!!$ write(0,*) me,' Sending to ',ip,nrl,kp+1,kp+nrl
call psb_snd(ctxt,rrcv(kp+1:kp+nrl),ip)
kp = kp + nrl
end do
end if
nrl = lv%remap_data%desc_ac_pre_remap%get_local_rows()
call psb_geall(tv,lv%remap_data%desc_ac_pre_remap,info)
!!$ write(0,*) me, ' Allocated ',nrl,info,psb_errstatus_fatal()
call psb_geasb(tv,lv%remap_data%desc_ac_pre_remap,info)
!!$ write(0,*) me,' Size of TV ',nrl,tv%get_nrows(),info
!!$ write(0,*) me,' Receiving from ',idest,nrl,psb_errstatus_fatal()
call psb_realloc(nrc,rsnd,info)
call psb_rcv(ctxt,rsnd(1:nrl),idest)
call tv%set_vect(rsnd)
call lv%linmap%map_V2U(alpha,tv,beta,vect_u,info,&
& work=work,vtx=vtx,vty=vty)
end associate
!!$ write(0,*) me, ' Prolongator with remap done '
!!$ flush(0)
!!$ call psb_barrier(ctxt)
end block
else else
! Default transfer ! Default transfer
call lv%linmap%map_V2U(alpha,vect_v,beta,vect_u,info,& call lv%linmap%map_V2U(alpha,vect_v,beta,vect_u,info,&

43
amgprec/impl/level/amg_s_base_onelev_map_rstr.F90
Unescape Escape View File

@ -36,7 +36,8 @@
! !
! !
subroutine amg_s_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,work,vtx,vty) subroutine amg_s_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
& work,vtx,vty)
use psb_base_mod use psb_base_mod
use amg_s_onelev_mod, amg_protect_name => amg_s_base_onelev_map_rstr_v use amg_s_onelev_mod, amg_protect_name => amg_s_base_onelev_map_rstr_v
implicit none implicit none
@ -47,24 +48,52 @@ subroutine amg_s_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,work,vt
real(psb_spk_), optional :: work(:) real(psb_spk_), optional :: work(:)
type(psb_s_vect_type), optional, target, intent(inout) :: vtx,vty type(psb_s_vect_type), optional, target, intent(inout) :: vtx,vty
!!$ write(0,*) 'New map_rstr',lv%remap_data%ac_pre_remap%is_asb()
if (lv%remap_data%ac_pre_remap%is_asb()) then if (lv%remap_data%ac_pre_remap%is_asb()) then
! !
! Remap has happened, deal with it ! Remap has happened, deal with it
! !
write(0,*) 'Remap handling not implemented yet ' !!$ write(0,*) 'Remap handling not implemented yet '
block block
integer(psb_ipk_) :: i,j,ip,nctxt,ictxt, idest type(psb_ctxt_type) :: ctxt, nctxt
integer(psb_ipk_) :: i,j,ip, idest, nsrc, nrl, kp
integer(psb_ipk_) :: me, np, rme, rnp integer(psb_ipk_) :: me, np, rme, rnp
real(psb_spk_), allocatable :: rsnd(:), rrcv(:)
type(psb_s_vect_type) :: tv
ictxt = lv%remap_data%desc_ac_pre_remap%get_ctxt() ctxt = lv%remap_data%desc_ac_pre_remap%get_ctxt()
call psb_info(ictxt,me,np) call psb_info(ctxt,me,np)
nctxt = lv%desc_ac%get_ctxt() nctxt = lv%desc_ac%get_ctxt()
call psb_info(nctxt,rme,rnp) call psb_info(nctxt,rme,rnp)
!!$ write(0,*) 'New context ',rme,rnp
idest = lv%remap_data%idest idest = lv%remap_data%idest
associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc) associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc)
write(0,*) 'Should apply maps, then send data from ',me,' to ',idest !!$ write(0,*) 'Should apply maps, then send data from ',me,' to ',idest
if (rme >= 0) write(0,*) rme, ' Receiving data from ',isrc(:) !!$ if (rme >= 0) write(0,*) rme, ' Receiving data from ',isrc(:)
nsrc = size(isrc)
nrl = lv%remap_data%desc_ac_pre_remap%get_local_rows()
call psb_geall(tv,lv%remap_data%desc_ac_pre_remap,info)
call psb_geasb(tv,lv%remap_data%desc_ac_pre_remap,info)
!!$ write(0,*) me,' Size of TV ',tv%get_nrows()
call lv%linmap%map_U2V(alpha,vect_u,beta,tv,info,&
& work=work,vtx=vtx,vty=vty)
rsnd = tv%get_vect()
call psb_snd(ctxt,rsnd(1:nrl),idest)
if (rme >=0) then
allocate(rrcv(sum(nrsrc)))
!!$ write(0,*) me,rme,' Size check ',size(rrcv)!,lv%desc_ac%get_local_rows()
kp = 0
do i = 1,size(isrc)
ip = isrc(i)
nrl = nrsrc(i)
!!$ write(0,*) me,' Receiving from ',ip,nrl,kp+1,kp+nrl,size(rrcv)
call psb_rcv(ctxt,rrcv(kp+1:kp+nrl),ip)
kp = kp + nrl
end do
call vect_v%set_vect(rrcv)
end if
end associate end associate
!!$ write(0,*) me, ' Restrictor with remap done '
end block end block
else else

129
amgprec/impl/level/amg_z_base_onelev_build.f90
Unescape Escape View File

@ -71,72 +71,81 @@ subroutine amg_z_base_onelev_build(lv,info,amold,vmold,imold,ilv)
ctxt = lv%base_desc%get_ctxt() ctxt = lv%base_desc%get_ctxt()
call psb_info(ctxt,me,np) call psb_info(ctxt,me,np)
if (.not.allocated(lv%sm)) then !
!! Error: should have called amg_dprecinit ! At top level(s) I may be using
info=3111 ! a context with less processes
call psb_errpush(info,name) !
goto 9999 if (me < 0) then
end if !!$ write(0,*) 'onelevbld: I am excluded from this one '
if (.not.allocated(lv%sm%sv)) then else
!! Error: should have called amg_dprecinit !!$ write(0,*) me,' Going to build smoothers at this level '
info=3111 if (.not.allocated(lv%sm)) then
call psb_errpush(info,name) !! Error: should have called amg_dprecinit
goto 9999 info=3111
end if call psb_errpush(info,name)
lv%ac_nz_loc = lv%ac%get_nzeros() goto 9999
lv%ac_nz_tot = lv%ac_nz_loc end if
select case(lv%parms%coarse_mat) if (.not.allocated(lv%sm%sv)) then
case(amg_distr_mat_) !! Error: should have called amg_dprecinit
call psb_sum(ctxt,lv%ac_nz_tot) info=3111
case(amg_repl_mat_) call psb_errpush(info,name)
! Do nothing goto 9999
case default end if
! Should never get here lv%ac_nz_loc = lv%ac%get_nzeros()
info = psb_err_internal_error_ lv%ac_nz_tot = lv%ac_nz_loc
call psb_errpush(info,name,a_err='Wrong lv%parms') select case(lv%parms%coarse_mat)
goto 9999 case(amg_distr_mat_)
end select call psb_sum(ctxt,lv%ac_nz_tot)
case(amg_repl_mat_)
! Do nothing
case default
! Should never get here
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='Wrong lv%parms')
goto 9999
end select
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Calling mlprcbld at level ',i & 'Calling mlprcbld at level ',i
call amg_check_def(lv%parms%sweeps_pre,& call amg_check_def(lv%parms%sweeps_pre,&
& 'Jacobi sweeps',izero,is_int_non_negative) & 'Jacobi sweeps',izero,is_int_non_negative)
call amg_check_def(lv%parms%sweeps_post,& call amg_check_def(lv%parms%sweeps_post,&
& 'Jacobi sweeps',izero,is_int_non_negative) & 'Jacobi sweeps',izero,is_int_non_negative)
call lv%sm%build(lv%base_a,lv%base_desc,info) call lv%sm%build(lv%base_a,lv%base_desc,info)
if (info == 0) then if (info == 0) then
if (allocated(lv%sm2a)) then if (allocated(lv%sm2a)) then
call lv%sm2a%build(lv%base_a,lv%base_desc,info) call lv%sm2a%build(lv%base_a,lv%base_desc,info)
lv%sm2 => lv%sm2a lv%sm2 => lv%sm2a
else else
lv%sm2 => lv%sm lv%sm2 => lv%sm
end if
end if
if (info /=0 ) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Smoother bld error')
goto 9999
end if end if
end if
if (info /=0 ) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Smoother bld error')
goto 9999
end if
if (lv%sm%sv%is_global()) then if (lv%sm%sv%is_global()) then
if ((lv%parms%sweeps_pre>1).or.(lv%parms%sweeps_post>1)) then if ((lv%parms%sweeps_pre>1).or.(lv%parms%sweeps_post>1)) then
lv%parms%sweeps_pre = 1 lv%parms%sweeps_pre = 1
lv%parms%sweeps_post = 1 lv%parms%sweeps_post = 1
if (me == 0) then if (me == 0) then
write(debug_unit,*) write(debug_unit,*)
if (present(ilv)) then if (present(ilv)) then
write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),& write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),&
& '" at level ',ilv & '" at level ',ilv
write(debug_unit,*) ' is configured as a global solver ' write(debug_unit,*) ' is configured as a global solver '
else else
write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),& write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),&
& '" is configured as a global solver ' & '" is configured as a global solver '
end if
write(debug_unit,*) ' Pre and post sweeps at this level reset to 1'
end if end if
write(debug_unit,*) ' Pre and post sweeps at this level reset to 1'
end if end if
end if end if
end if end if

52
amgprec/impl/level/amg_z_base_onelev_map_prol.F90
Unescape Escape View File

@ -47,11 +47,61 @@ subroutine amg_z_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt
complex(psb_dpk_), optional :: work(:) complex(psb_dpk_), optional :: work(:)
type(psb_z_vect_type), optional, target, intent(inout) :: vtx,vty type(psb_z_vect_type), optional, target, intent(inout) :: vtx,vty
!!$ write(0,*) 'New map_prol',lv%remap_data%ac_pre_remap%is_asb()
if (lv%remap_data%ac_pre_remap%is_asb()) then if (lv%remap_data%ac_pre_remap%is_asb()) then
! !
! Remap has happened, deal with it ! Remap has happened, deal with it
! !
write(0,*) 'Remap handling not implemented yet ' !!$ write(0,*) 'Remap handling '
block
type(psb_ctxt_type) :: ctxt, nctxt
integer(psb_ipk_) :: i,j,ip,idest, nsrc, nrl, nrc, kp
integer(psb_ipk_) :: me, np, rme, rnp
complex(psb_dpk_), allocatable :: rsnd(:), rrcv(:)
type(psb_z_vect_type) :: tv
ctxt = lv%remap_data%desc_ac_pre_remap%get_ctxt()
call psb_info(ctxt,me,np)
!!$ write(0,*) 'Old context ',me,np,psb_errstatus_fatal()
nctxt = lv%desc_ac%get_ctxt()
call psb_info(nctxt,rme,rnp)
!!$ write(0,*) 'New context ',rme,rnp,psb_errstatus_fatal()
idest = lv%remap_data%idest
associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc)
!!$ write(0,*) 'Should apply maps, then receive data from ',idest,' to ',me,psb_errstatus_fatal()
nsrc = size(isrc)
nrl = lv%remap_data%desc_ac_pre_remap%get_local_rows()
nrc = lv%remap_data%desc_ac_pre_remap%get_local_cols()
if (rme >=0) then
allocate(rrcv(sum(nrsrc)))
rrcv = vect_v%get_vect()
!!$ write(0,*) me,rme,' Size check ',size(rrcv),lv%desc_ac%get_local_rows(),psb_errstatus_fatal()
kp = 0
do i = 1,size(isrc)
ip = isrc(i)
nrl = nrsrc(i)
!!$ write(0,*) me,' Sending to ',ip,nrl,kp+1,kp+nrl
call psb_snd(ctxt,rrcv(kp+1:kp+nrl),ip)
kp = kp + nrl
end do
end if
nrl = lv%remap_data%desc_ac_pre_remap%get_local_rows()
call psb_geall(tv,lv%remap_data%desc_ac_pre_remap,info)
!!$ write(0,*) me, ' Allocated ',nrl,info,psb_errstatus_fatal()
call psb_geasb(tv,lv%remap_data%desc_ac_pre_remap,info)
!!$ write(0,*) me,' Size of TV ',nrl,tv%get_nrows(),info
!!$ write(0,*) me,' Receiving from ',idest,nrl,psb_errstatus_fatal()
call psb_realloc(nrc,rsnd,info)
call psb_rcv(ctxt,rsnd(1:nrl),idest)
call tv%set_vect(rsnd)
call lv%linmap%map_V2U(alpha,tv,beta,vect_u,info,&
& work=work,vtx=vtx,vty=vty)
end associate
!!$ write(0,*) me, ' Prolongator with remap done '
!!$ flush(0)
!!$ call psb_barrier(ctxt)
end block
else else
! Default transfer ! Default transfer
call lv%linmap%map_V2U(alpha,vect_v,beta,vect_u,info,& call lv%linmap%map_V2U(alpha,vect_v,beta,vect_u,info,&

43
amgprec/impl/level/amg_z_base_onelev_map_rstr.F90
Unescape Escape View File

@ -36,7 +36,8 @@
! !
! !
subroutine amg_z_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,work,vtx,vty) subroutine amg_z_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,&
& work,vtx,vty)
use psb_base_mod use psb_base_mod
use amg_z_onelev_mod, amg_protect_name => amg_z_base_onelev_map_rstr_v use amg_z_onelev_mod, amg_protect_name => amg_z_base_onelev_map_rstr_v
implicit none implicit none
@ -47,24 +48,52 @@ subroutine amg_z_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,work,vt
complex(psb_dpk_), optional :: work(:) complex(psb_dpk_), optional :: work(:)
type(psb_z_vect_type), optional, target, intent(inout) :: vtx,vty type(psb_z_vect_type), optional, target, intent(inout) :: vtx,vty
!!$ write(0,*) 'New map_rstr',lv%remap_data%ac_pre_remap%is_asb()
if (lv%remap_data%ac_pre_remap%is_asb()) then if (lv%remap_data%ac_pre_remap%is_asb()) then
! !
! Remap has happened, deal with it ! Remap has happened, deal with it
! !
write(0,*) 'Remap handling not implemented yet ' !!$ write(0,*) 'Remap handling not implemented yet '
block block
integer(psb_ipk_) :: i,j,ip,nctxt,ictxt, idest type(psb_ctxt_type) :: ctxt, nctxt
integer(psb_ipk_) :: i,j,ip, idest, nsrc, nrl, kp
integer(psb_ipk_) :: me, np, rme, rnp integer(psb_ipk_) :: me, np, rme, rnp
complex(psb_dpk_), allocatable :: rsnd(:), rrcv(:)
type(psb_z_vect_type) :: tv
ictxt = lv%remap_data%desc_ac_pre_remap%get_ctxt() ctxt = lv%remap_data%desc_ac_pre_remap%get_ctxt()
call psb_info(ictxt,me,np) call psb_info(ctxt,me,np)
nctxt = lv%desc_ac%get_ctxt() nctxt = lv%desc_ac%get_ctxt()
call psb_info(nctxt,rme,rnp) call psb_info(nctxt,rme,rnp)
!!$ write(0,*) 'New context ',rme,rnp
idest = lv%remap_data%idest idest = lv%remap_data%idest
associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc) associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc)
write(0,*) 'Should apply maps, then send data from ',me,' to ',idest !!$ write(0,*) 'Should apply maps, then send data from ',me,' to ',idest
if (rme >= 0) write(0,*) rme, ' Receiving data from ',isrc(:) !!$ if (rme >= 0) write(0,*) rme, ' Receiving data from ',isrc(:)
nsrc = size(isrc)
nrl = lv%remap_data%desc_ac_pre_remap%get_local_rows()
call psb_geall(tv,lv%remap_data%desc_ac_pre_remap,info)
call psb_geasb(tv,lv%remap_data%desc_ac_pre_remap,info)
!!$ write(0,*) me,' Size of TV ',tv%get_nrows()
call lv%linmap%map_U2V(alpha,vect_u,beta,tv,info,&
& work=work,vtx=vtx,vty=vty)
rsnd = tv%get_vect()
call psb_snd(ctxt,rsnd(1:nrl),idest)
if (rme >=0) then
allocate(rrcv(sum(nrsrc)))
!!$ write(0,*) me,rme,' Size check ',size(rrcv)!,lv%desc_ac%get_local_rows()
kp = 0
do i = 1,size(isrc)
ip = isrc(i)
nrl = nrsrc(i)
!!$ write(0,*) me,' Receiving from ',ip,nrl,kp+1,kp+nrl,size(rrcv)
call psb_rcv(ctxt,rrcv(kp+1:kp+nrl),ip)
kp = kp + nrl
end do
call vect_v%set_vect(rrcv)
end if
end associate end associate
!!$ write(0,*) me, ' Restrictor with remap done '
end block end block
else else

Write Preview
Loading…
Cancel
Save