diff --git a/amgprec/amg_c_onelev_mod.f90 b/amgprec/amg_c_onelev_mod.f90 index 11f283ea..66e3632f 100644 --- a/amgprec/amg_c_onelev_mod.f90 +++ b/amgprec/amg_c_onelev_mod.f90 @@ -157,7 +157,7 @@ module amg_c_onelev_mod type(psb_cspmat_type) :: ac_pre_remap type(psb_desc_type) :: desc_ac_pre_remap integer(psb_ipk_) :: idest - integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:) + integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:) contains procedure, pass(rmp) :: clone => c_remap_data_clone end type amg_c_remap_data_type @@ -704,7 +704,17 @@ contains info = psb_success_ nwv = lv%get_wrksz() 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 @@ -724,7 +734,7 @@ contains end if 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 Implicit None @@ -735,25 +745,67 @@ contains type(psb_desc_type), intent(in) :: desc integer(psb_ipk_), intent(out) :: info class(psb_c_base_vect_type), intent(in), optional :: vmold + type(psb_desc_type), intent(in), optional :: desc2 ! integer(psb_ipk_) :: i info = psb_success_ call wk%free(info) - 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,& + if (present(desc2)) then +!!$ write(0,*) 'Check on wrk_alloc 2',& +!!$ & desc2%get_local_rows(), desc%get_local_rows(),& +!!$ & desc2%get_local_cols(),desc%get_local_cols() +!!$ flush(0) + if (desc2%get_local_cols()>desc%get_local_cols()) then + call psb_geasb(wk%vx2l,desc2,info,& + & scratch=.true.,mold=vmold) + call psb_geasb(wk%vy2l,desc2,info,& + & scratch=.true.,mold=vmold) + 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) - 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 subroutine c_wrk_free(wk,info) diff --git a/amgprec/amg_c_prec_type.f90 b/amgprec/amg_c_prec_type.f90 index de049810..db885e54 100644 --- a/amgprec/amg_c_prec_type.f90 +++ b/amgprec/amg_c_prec_type.f90 @@ -738,15 +738,15 @@ contains character(len=*), intent(in), optional :: prefix, head logical, optional, intent(in) :: smoother, solver,ac, rp, tprol, global_num integer(psb_ipk_) :: i, j, il1, iln, lev - type(psb_ctxt_type) :: icontxt + type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: iam, np, iproc_ character(len=80) :: prefix_ character(len=120) :: fname ! len should be at least 20 more than ! len of prefix_ info = 0 - icontxt = prec%ctxt - call psb_info(icontxt,iam,np) + ctxt = prec%ctxt + call psb_info(ctxt,iam,np) iln = size(prec%precv) if (present(istart)) then @@ -812,13 +812,13 @@ contains integer(psb_ipk_), intent(out) :: info ! Local vars integer(psb_ipk_) :: i, j, ln, lev - type(psb_ctxt_type) :: icontxt + type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: iam, np info = psb_success_ select type(pout => precout) class is (amg_cprec_type) - pout%ctxt = prec%ctxt + pout%ctxt = prec%ctxt pout%ag_data = prec%ag_data pout%outer_sweeps = prec%outer_sweeps if (allocated(prec%precv)) then diff --git a/amgprec/amg_d_onelev_mod.f90 b/amgprec/amg_d_onelev_mod.f90 index 511c3986..b0a92419 100644 --- a/amgprec/amg_d_onelev_mod.f90 +++ b/amgprec/amg_d_onelev_mod.f90 @@ -157,7 +157,7 @@ module amg_d_onelev_mod type(psb_dspmat_type) :: ac_pre_remap type(psb_desc_type) :: desc_ac_pre_remap integer(psb_ipk_) :: idest - integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:) + integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:) contains procedure, pass(rmp) :: clone => d_remap_data_clone end type amg_d_remap_data_type @@ -704,7 +704,17 @@ contains info = psb_success_ nwv = lv%get_wrksz() 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 @@ -724,7 +734,7 @@ contains end if 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 Implicit None @@ -735,25 +745,67 @@ contains type(psb_desc_type), intent(in) :: desc integer(psb_ipk_), intent(out) :: info class(psb_d_base_vect_type), intent(in), optional :: vmold + type(psb_desc_type), intent(in), optional :: desc2 ! integer(psb_ipk_) :: i info = psb_success_ call wk%free(info) - 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,& + if (present(desc2)) then +!!$ write(0,*) 'Check on wrk_alloc 2',& +!!$ & desc2%get_local_rows(), desc%get_local_rows(),& +!!$ & desc2%get_local_cols(),desc%get_local_cols() +!!$ flush(0) + if (desc2%get_local_cols()>desc%get_local_cols()) then + call psb_geasb(wk%vx2l,desc2,info,& + & scratch=.true.,mold=vmold) + call psb_geasb(wk%vy2l,desc2,info,& + & scratch=.true.,mold=vmold) + 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) - 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 subroutine d_wrk_free(wk,info) diff --git a/amgprec/amg_d_prec_type.f90 b/amgprec/amg_d_prec_type.f90 index cca6a067..9d383f70 100644 --- a/amgprec/amg_d_prec_type.f90 +++ b/amgprec/amg_d_prec_type.f90 @@ -738,15 +738,15 @@ contains character(len=*), intent(in), optional :: prefix, head logical, optional, intent(in) :: smoother, solver,ac, rp, tprol, global_num integer(psb_ipk_) :: i, j, il1, iln, lev - type(psb_ctxt_type) :: icontxt + type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: iam, np, iproc_ character(len=80) :: prefix_ character(len=120) :: fname ! len should be at least 20 more than ! len of prefix_ info = 0 - icontxt = prec%ctxt - call psb_info(icontxt,iam,np) + ctxt = prec%ctxt + call psb_info(ctxt,iam,np) iln = size(prec%precv) if (present(istart)) then @@ -812,13 +812,13 @@ contains integer(psb_ipk_), intent(out) :: info ! Local vars integer(psb_ipk_) :: i, j, ln, lev - type(psb_ctxt_type) :: icontxt + type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: iam, np info = psb_success_ select type(pout => precout) class is (amg_dprec_type) - pout%ctxt = prec%ctxt + pout%ctxt = prec%ctxt pout%ag_data = prec%ag_data pout%outer_sweeps = prec%outer_sweeps if (allocated(prec%precv)) then diff --git a/amgprec/amg_s_onelev_mod.f90 b/amgprec/amg_s_onelev_mod.f90 index 06d13ef5..8c606fca 100644 --- a/amgprec/amg_s_onelev_mod.f90 +++ b/amgprec/amg_s_onelev_mod.f90 @@ -157,7 +157,7 @@ module amg_s_onelev_mod type(psb_sspmat_type) :: ac_pre_remap type(psb_desc_type) :: desc_ac_pre_remap integer(psb_ipk_) :: idest - integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:) + integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:) contains procedure, pass(rmp) :: clone => s_remap_data_clone end type amg_s_remap_data_type @@ -704,7 +704,17 @@ contains info = psb_success_ nwv = lv%get_wrksz() 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 @@ -724,7 +734,7 @@ contains end if 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 Implicit None @@ -735,25 +745,67 @@ contains type(psb_desc_type), intent(in) :: desc integer(psb_ipk_), intent(out) :: info class(psb_s_base_vect_type), intent(in), optional :: vmold + type(psb_desc_type), intent(in), optional :: desc2 ! integer(psb_ipk_) :: i info = psb_success_ call wk%free(info) - 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,& + if (present(desc2)) then +!!$ write(0,*) 'Check on wrk_alloc 2',& +!!$ & desc2%get_local_rows(), desc%get_local_rows(),& +!!$ & desc2%get_local_cols(),desc%get_local_cols() +!!$ flush(0) + if (desc2%get_local_cols()>desc%get_local_cols()) then + call psb_geasb(wk%vx2l,desc2,info,& + & scratch=.true.,mold=vmold) + call psb_geasb(wk%vy2l,desc2,info,& + & scratch=.true.,mold=vmold) + 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) - 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 subroutine s_wrk_free(wk,info) diff --git a/amgprec/amg_s_prec_type.f90 b/amgprec/amg_s_prec_type.f90 index 05b29fd0..f449f958 100644 --- a/amgprec/amg_s_prec_type.f90 +++ b/amgprec/amg_s_prec_type.f90 @@ -738,15 +738,15 @@ contains character(len=*), intent(in), optional :: prefix, head logical, optional, intent(in) :: smoother, solver,ac, rp, tprol, global_num integer(psb_ipk_) :: i, j, il1, iln, lev - type(psb_ctxt_type) :: icontxt + type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: iam, np, iproc_ character(len=80) :: prefix_ character(len=120) :: fname ! len should be at least 20 more than ! len of prefix_ info = 0 - icontxt = prec%ctxt - call psb_info(icontxt,iam,np) + ctxt = prec%ctxt + call psb_info(ctxt,iam,np) iln = size(prec%precv) if (present(istart)) then @@ -812,13 +812,13 @@ contains integer(psb_ipk_), intent(out) :: info ! Local vars integer(psb_ipk_) :: i, j, ln, lev - type(psb_ctxt_type) :: icontxt + type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: iam, np info = psb_success_ select type(pout => precout) class is (amg_sprec_type) - pout%ctxt = prec%ctxt + pout%ctxt = prec%ctxt pout%ag_data = prec%ag_data pout%outer_sweeps = prec%outer_sweeps if (allocated(prec%precv)) then diff --git a/amgprec/amg_z_onelev_mod.f90 b/amgprec/amg_z_onelev_mod.f90 index 6328218c..811c757d 100644 --- a/amgprec/amg_z_onelev_mod.f90 +++ b/amgprec/amg_z_onelev_mod.f90 @@ -157,7 +157,7 @@ module amg_z_onelev_mod type(psb_zspmat_type) :: ac_pre_remap type(psb_desc_type) :: desc_ac_pre_remap integer(psb_ipk_) :: idest - integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:) + integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:) contains procedure, pass(rmp) :: clone => z_remap_data_clone end type amg_z_remap_data_type @@ -704,7 +704,17 @@ contains info = psb_success_ nwv = lv%get_wrksz() 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 @@ -724,7 +734,7 @@ contains end if 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 Implicit None @@ -735,25 +745,67 @@ contains type(psb_desc_type), intent(in) :: desc integer(psb_ipk_), intent(out) :: info class(psb_z_base_vect_type), intent(in), optional :: vmold + type(psb_desc_type), intent(in), optional :: desc2 ! integer(psb_ipk_) :: i info = psb_success_ call wk%free(info) - 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,& + if (present(desc2)) then +!!$ write(0,*) 'Check on wrk_alloc 2',& +!!$ & desc2%get_local_rows(), desc%get_local_rows(),& +!!$ & desc2%get_local_cols(),desc%get_local_cols() +!!$ flush(0) + if (desc2%get_local_cols()>desc%get_local_cols()) then + call psb_geasb(wk%vx2l,desc2,info,& + & scratch=.true.,mold=vmold) + call psb_geasb(wk%vy2l,desc2,info,& + & scratch=.true.,mold=vmold) + 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) - 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 subroutine z_wrk_free(wk,info) diff --git a/amgprec/amg_z_prec_type.f90 b/amgprec/amg_z_prec_type.f90 index b16cea5b..533c3ec3 100644 --- a/amgprec/amg_z_prec_type.f90 +++ b/amgprec/amg_z_prec_type.f90 @@ -738,15 +738,15 @@ contains character(len=*), intent(in), optional :: prefix, head logical, optional, intent(in) :: smoother, solver,ac, rp, tprol, global_num integer(psb_ipk_) :: i, j, il1, iln, lev - type(psb_ctxt_type) :: icontxt + type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: iam, np, iproc_ character(len=80) :: prefix_ character(len=120) :: fname ! len should be at least 20 more than ! len of prefix_ info = 0 - icontxt = prec%ctxt - call psb_info(icontxt,iam,np) + ctxt = prec%ctxt + call psb_info(ctxt,iam,np) iln = size(prec%precv) if (present(istart)) then @@ -812,13 +812,13 @@ contains integer(psb_ipk_), intent(out) :: info ! Local vars integer(psb_ipk_) :: i, j, ln, lev - type(psb_ctxt_type) :: icontxt + type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: iam, np info = psb_success_ select type(pout => precout) class is (amg_zprec_type) - pout%ctxt = prec%ctxt + pout%ctxt = prec%ctxt pout%ag_data = prec%ag_data pout%outer_sweeps = prec%outer_sweeps if (allocated(prec%precv)) then diff --git a/amgprec/impl/amg_c_hierarchy_bld.f90 b/amgprec/impl/amg_c_hierarchy_bld.f90 index 4b59c562..cfb7e53d 100644 --- a/amgprec/impl/amg_c_hierarchy_bld.f90 +++ b/amgprec/impl/amg_c_hierarchy_bld.f90 @@ -455,6 +455,31 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info) end do 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 call psb_errpush(psb_err_internal_error_,name,& & a_err='Internal hierarchy build' ) diff --git a/amgprec/impl/amg_c_smoothers_bld.f90 b/amgprec/impl/amg_c_smoothers_bld.f90 index 5dffeec6..bb5907a1 100644 --- a/amgprec/impl/amg_c_smoothers_bld.f90 +++ b/amgprec/impl/amg_c_smoothers_bld.f90 @@ -117,6 +117,10 @@ subroutine amg_c_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold) info = psb_success_ ctxt = desc_a%get_context() 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_) & & 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 ! +!!$ write(0,*) me,' Building at level ',i call prec%precv(i)%bld(info,amold=amold,vmold=vmold,imold=imold,ilv=i) 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),& & 'Exiting with',iszv,' levels' +9998 continue call psb_erractionrestore(err_act) return diff --git a/amgprec/impl/amg_cmlprec_aply.f90 b/amgprec/impl/amg_cmlprec_aply.f90 index f669caeb..32c525d5 100644 --- a/amgprec/impl/amg_cmlprec_aply.f90 +++ b/amgprec/impl/amg_cmlprec_aply.f90 @@ -393,38 +393,39 @@ contains if(debug_level > 1) then write(debug_unit,*) me,' Start inner_ml_aply at level ',level, info end if - - select case(p%precv(level)%parms%ml_cycle) - - case(amg_no_ml_) - ! - ! No preconditioning, should not really get here - ! - call psb_errpush(psb_err_internal_error_,name,& - & a_err='amg_no_ml_ in mlprc_aply?') - goto 9999 - - case(amg_add_ml_) - - call amg_c_inner_add(p, level, trans, work) - - case(amg_mult_ml_,amg_vcycle_ml_, amg_wcycle_ml_) - - call amg_c_inner_mult(p, level, trans, work) - - case(amg_kcycle_ml_, amg_kcyclesym_ml_) - - call amg_c_inner_k_cycle(p, level, trans, work) - - case default - info = psb_err_from_subroutine_ai_ - call psb_errpush(info,name,a_err='invalid ml_cycle',& - & i_Err=(/p%precv(level)%parms%ml_cycle,izero,izero,izero,izero/)) - goto 9999 - - end select - if(debug_level > 1) then - write(debug_unit,*) me,' End inner_ml_aply at level ',level + if (me >= 0) then + select case(p%precv(level)%parms%ml_cycle) + + case(amg_no_ml_) + ! + ! No preconditioning, should not really get here + ! + call psb_errpush(psb_err_internal_error_,name,& + & a_err='amg_no_ml_ in mlprc_aply?') + goto 9999 + + case(amg_add_ml_) + + call amg_c_inner_add(p, level, trans, work) + + case(amg_mult_ml_,amg_vcycle_ml_, amg_wcycle_ml_) + + call amg_c_inner_mult(p, level, trans, work) + + case(amg_kcycle_ml_, amg_kcyclesym_ml_) + + call amg_c_inner_k_cycle(p, level, trans, work) + + case default + info = psb_err_from_subroutine_ai_ + call psb_errpush(info,name,a_err='invalid ml_cycle',& + & i_Err=(/p%precv(level)%parms%ml_cycle,izero,izero,izero,izero/)) + goto 9999 + + end select + if(debug_level > 1) then + write(debug_unit,*) me,' End inner_ml_aply at level ',level + end if end if call psb_erractionrestore(err_act) @@ -492,28 +493,30 @@ contains & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & wv => p%precv(level)%wrk%wv) - 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) - do k=1, sweeps + if (me >= 0) then + 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) + do k=1, sweeps + call p%precv(level)%sm%apply(cone,& + & vy2l,czero,vty,& + & base_desc, trans,& + & ione,work,wv,info,init='Z') + + call p%precv(level)%sm2a%apply(cone,& + & vty,czero,vy2l,& + & base_desc, trans,& + & ione,work,wv,info,init='Z') + end do + + else + sweeps = p%precv(level)%parms%sweeps_pre call p%precv(level)%sm%apply(cone,& - & vy2l,czero,vty,& - & base_desc, trans,& - & ione,work,wv,info,init='Z') - - call p%precv(level)%sm2a%apply(cone,& - & vty,czero,vy2l,& + & vx2l,czero,vy2l,& & base_desc, trans,& - & ione,work,wv,info,init='Z') - end do - - else - sweeps = p%precv(level)%parms%sweeps_pre - call p%precv(level)%sm%apply(cone,& - & vx2l,czero,vy2l,& - & base_desc, trans,& - & sweeps,work,wv,info,init='Z') + & sweeps,work,wv,info,init='Z') + end if end if if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -602,7 +605,6 @@ contains end if ctxt = p%precv(level)%base_desc%get_context() call psb_info(ctxt, me, np) - if(debug_level > 1) then write(debug_unit,*) me,' inner_mult at level ',level end if @@ -623,22 +625,25 @@ contains ! if (pre) 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 PRE smoother_apply') - goto 9999 + if (me >=0) then +!!$ write(0,*) me,'Applying smoother pre ', level + 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 PRE smoother_apply') + goto 9999 + end if end if endif ! @@ -696,11 +701,13 @@ contains if (p%precv(level)%parms%ml_cycle == amg_wcycle_ml_) then - call psb_geaxpby(cone,vx2l, czero,vty,& - & base_desc,info) - if (info == psb_success_) call psb_spmm(-cone,base_a,& - & vy2l,cone,vty,& - & base_desc,info,work=work,trans=trans) + if (me >=0) then + call psb_geaxpby(cone,vx2l, czero,vty,& + & base_desc,info) + if (info == psb_success_) call psb_spmm(-cone,base_a,& + & vy2l,cone,vty,& + & base_desc,info,work=work,trans=trans) + end if if (info == psb_success_) & & call p%precv(level+1)%map_rstr(cone,vty,& & czero,p%precv(level+1)%wrk%vx2l,info,work=work,& @@ -728,31 +735,33 @@ contains if (post) then - call psb_geaxpby(cone,vx2l,& - & czero,vty,& - & base_desc,info) - if (info == psb_success_) 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 second 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') + if (me >=0) then + call psb_geaxpby(cone,vx2l,& + & czero,vty,& + & base_desc,info) + if (info == psb_success_) 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 second 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 @@ -764,11 +773,14 @@ contains endif else if (level == nlev) 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) +!!$ write(0,*) me,'Applying smoother at top level ',psb_errstatus_fatal() + if (me >=0) 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) + end if +!!$ write(0,*) me,' Done applying smoother at top level ',psb_errstatus_fatal() else @@ -778,7 +790,7 @@ contains goto 9999 end if end associate - +9998 continue call psb_erractionrestore(err_act) return @@ -829,70 +841,71 @@ contains end if ctxt = p%precv(level)%base_desc%get_context() call psb_info(ctxt, me, np) - + if(debug_level > 1) then write(debug_unit,*) me,name,' start at level ',level end if - + if ((level<1).or.(level>nlev)) then info = psb_err_internal_error_ call psb_errpush(info,name,& & a_err='Invalid LEVEL>NLEV') goto 9999 end if - + !K cycle associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& & 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,& & 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 (trans == 'N') then + if (level == nlev) then + if (me >= 0) 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 - 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 + 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,& - & czero,vty,& - & base_desc,info) + call psb_geaxpby(cone,vx2l,& + & czero,vty,& + & base_desc,info) - if (info == psb_success_) 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 + if (info == psb_success_) 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 end if - ! Apply the restriction call p%precv(level + 1)%map_rstr(cone,vty,& & czero,p%precv(level + 1)%wrk%vx2l,& @@ -940,40 +953,42 @@ contains & a_err='Error during prolongation') goto 9999 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 - ! - ! 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') + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during POST smoother_apply') + goto 9999 + 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 info = psb_err_internal_error_ @@ -990,8 +1005,7 @@ contains return end subroutine amg_c_inner_k_cycle - - + recursive subroutine amg_cinneritkcycle(p, level, trans, work, innersolv) use psb_base_mod use amg_prec_mod diff --git a/amgprec/impl/amg_d_hierarchy_bld.f90 b/amgprec/impl/amg_d_hierarchy_bld.f90 index f6183fc8..c0140f92 100644 --- a/amgprec/impl/amg_d_hierarchy_bld.f90 +++ b/amgprec/impl/amg_d_hierarchy_bld.f90 @@ -455,6 +455,31 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info) end do 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 call psb_errpush(psb_err_internal_error_,name,& & a_err='Internal hierarchy build' ) diff --git a/amgprec/impl/amg_d_smoothers_bld.f90 b/amgprec/impl/amg_d_smoothers_bld.f90 index e9dd0d93..a71b6d85 100644 --- a/amgprec/impl/amg_d_smoothers_bld.f90 +++ b/amgprec/impl/amg_d_smoothers_bld.f90 @@ -117,6 +117,10 @@ subroutine amg_d_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold) info = psb_success_ ctxt = desc_a%get_context() 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_) & & 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 ! +!!$ write(0,*) me,' Building at level ',i call prec%precv(i)%bld(info,amold=amold,vmold=vmold,imold=imold,ilv=i) 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),& & 'Exiting with',iszv,' levels' +9998 continue call psb_erractionrestore(err_act) return diff --git a/amgprec/impl/amg_dmlprec_aply.f90 b/amgprec/impl/amg_dmlprec_aply.f90 index 4606643b..592ce2d4 100644 --- a/amgprec/impl/amg_dmlprec_aply.f90 +++ b/amgprec/impl/amg_dmlprec_aply.f90 @@ -393,38 +393,39 @@ contains if(debug_level > 1) then write(debug_unit,*) me,' Start inner_ml_aply at level ',level, info end if - - select case(p%precv(level)%parms%ml_cycle) - - case(amg_no_ml_) - ! - ! No preconditioning, should not really get here - ! - call psb_errpush(psb_err_internal_error_,name,& - & a_err='amg_no_ml_ in mlprc_aply?') - goto 9999 - - case(amg_add_ml_) - - call amg_d_inner_add(p, level, trans, work) - - case(amg_mult_ml_,amg_vcycle_ml_, amg_wcycle_ml_) - - call amg_d_inner_mult(p, level, trans, work) - - case(amg_kcycle_ml_, amg_kcyclesym_ml_) - - call amg_d_inner_k_cycle(p, level, trans, work) - - case default - info = psb_err_from_subroutine_ai_ - call psb_errpush(info,name,a_err='invalid ml_cycle',& - & i_Err=(/p%precv(level)%parms%ml_cycle,izero,izero,izero,izero/)) - goto 9999 - - end select - if(debug_level > 1) then - write(debug_unit,*) me,' End inner_ml_aply at level ',level + if (me >= 0) then + select case(p%precv(level)%parms%ml_cycle) + + case(amg_no_ml_) + ! + ! No preconditioning, should not really get here + ! + call psb_errpush(psb_err_internal_error_,name,& + & a_err='amg_no_ml_ in mlprc_aply?') + goto 9999 + + case(amg_add_ml_) + + call amg_d_inner_add(p, level, trans, work) + + case(amg_mult_ml_,amg_vcycle_ml_, amg_wcycle_ml_) + + call amg_d_inner_mult(p, level, trans, work) + + case(amg_kcycle_ml_, amg_kcyclesym_ml_) + + call amg_d_inner_k_cycle(p, level, trans, work) + + case default + info = psb_err_from_subroutine_ai_ + call psb_errpush(info,name,a_err='invalid ml_cycle',& + & i_Err=(/p%precv(level)%parms%ml_cycle,izero,izero,izero,izero/)) + goto 9999 + + end select + if(debug_level > 1) then + write(debug_unit,*) me,' End inner_ml_aply at level ',level + end if end if call psb_erractionrestore(err_act) @@ -492,28 +493,30 @@ contains & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & wv => p%precv(level)%wrk%wv) - 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) - do k=1, sweeps + if (me >= 0) then + 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) + do k=1, sweeps + call p%precv(level)%sm%apply(done,& + & vy2l,dzero,vty,& + & base_desc, trans,& + & ione,work,wv,info,init='Z') + + call p%precv(level)%sm2a%apply(done,& + & vty,dzero,vy2l,& + & base_desc, trans,& + & ione,work,wv,info,init='Z') + end do + + else + sweeps = p%precv(level)%parms%sweeps_pre call p%precv(level)%sm%apply(done,& - & vy2l,dzero,vty,& - & base_desc, trans,& - & ione,work,wv,info,init='Z') - - call p%precv(level)%sm2a%apply(done,& - & vty,dzero,vy2l,& + & vx2l,dzero,vy2l,& & base_desc, trans,& - & ione,work,wv,info,init='Z') - end do - - else - sweeps = p%precv(level)%parms%sweeps_pre - call p%precv(level)%sm%apply(done,& - & vx2l,dzero,vy2l,& - & base_desc, trans,& - & sweeps,work,wv,info,init='Z') + & sweeps,work,wv,info,init='Z') + end if end if if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -602,7 +605,6 @@ contains end if ctxt = p%precv(level)%base_desc%get_context() call psb_info(ctxt, me, np) - if(debug_level > 1) then write(debug_unit,*) me,' inner_mult at level ',level end if @@ -623,22 +625,25 @@ contains ! if (pre) 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 PRE smoother_apply') - goto 9999 + if (me >=0) then +!!$ write(0,*) me,'Applying smoother pre ', level + 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 PRE smoother_apply') + goto 9999 + end if end if endif ! @@ -696,11 +701,13 @@ contains if (p%precv(level)%parms%ml_cycle == amg_wcycle_ml_) then - call psb_geaxpby(done,vx2l, dzero,vty,& - & base_desc,info) - if (info == psb_success_) call psb_spmm(-done,base_a,& - & vy2l,done,vty,& - & base_desc,info,work=work,trans=trans) + if (me >=0) then + call psb_geaxpby(done,vx2l, dzero,vty,& + & base_desc,info) + if (info == psb_success_) call psb_spmm(-done,base_a,& + & vy2l,done,vty,& + & base_desc,info,work=work,trans=trans) + end if if (info == psb_success_) & & call p%precv(level+1)%map_rstr(done,vty,& & dzero,p%precv(level+1)%wrk%vx2l,info,work=work,& @@ -728,31 +735,33 @@ contains if (post) then - call psb_geaxpby(done,vx2l,& - & dzero,vty,& - & base_desc,info) - if (info == psb_success_) 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 second 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') + if (me >=0) then + call psb_geaxpby(done,vx2l,& + & dzero,vty,& + & base_desc,info) + if (info == psb_success_) 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 second 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 @@ -764,11 +773,14 @@ contains endif else if (level == nlev) 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) +!!$ write(0,*) me,'Applying smoother at top level ',psb_errstatus_fatal() + if (me >=0) 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) + end if +!!$ write(0,*) me,' Done applying smoother at top level ',psb_errstatus_fatal() else @@ -778,7 +790,7 @@ contains goto 9999 end if end associate - +9998 continue call psb_erractionrestore(err_act) return @@ -829,70 +841,71 @@ contains end if ctxt = p%precv(level)%base_desc%get_context() call psb_info(ctxt, me, np) - + if(debug_level > 1) then write(debug_unit,*) me,name,' start at level ',level end if - + if ((level<1).or.(level>nlev)) then info = psb_err_internal_error_ call psb_errpush(info,name,& & a_err='Invalid LEVEL>NLEV') goto 9999 end if - + !K cycle associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& & 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,& & 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 (trans == 'N') then + if (level == nlev) then + if (me >= 0) 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 - 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 + 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,& - & dzero,vty,& - & base_desc,info) + call psb_geaxpby(done,vx2l,& + & dzero,vty,& + & base_desc,info) - if (info == psb_success_) 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 + if (info == psb_success_) 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 end if - ! Apply the restriction call p%precv(level + 1)%map_rstr(done,vty,& & dzero,p%precv(level + 1)%wrk%vx2l,& @@ -940,40 +953,42 @@ contains & a_err='Error during prolongation') goto 9999 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 - ! - ! 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') + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during POST smoother_apply') + goto 9999 + 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 info = psb_err_internal_error_ @@ -990,8 +1005,7 @@ contains return end subroutine amg_d_inner_k_cycle - - + recursive subroutine amg_dinneritkcycle(p, level, trans, work, innersolv) use psb_base_mod use amg_prec_mod diff --git a/amgprec/impl/amg_s_hierarchy_bld.f90 b/amgprec/impl/amg_s_hierarchy_bld.f90 index edfe6add..93cad383 100644 --- a/amgprec/impl/amg_s_hierarchy_bld.f90 +++ b/amgprec/impl/amg_s_hierarchy_bld.f90 @@ -455,6 +455,31 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info) end do 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 call psb_errpush(psb_err_internal_error_,name,& & a_err='Internal hierarchy build' ) diff --git a/amgprec/impl/amg_s_smoothers_bld.f90 b/amgprec/impl/amg_s_smoothers_bld.f90 index ba54fd19..a6fcb596 100644 --- a/amgprec/impl/amg_s_smoothers_bld.f90 +++ b/amgprec/impl/amg_s_smoothers_bld.f90 @@ -117,6 +117,10 @@ subroutine amg_s_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold) info = psb_success_ ctxt = desc_a%get_context() 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_) & & 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 ! +!!$ write(0,*) me,' Building at level ',i call prec%precv(i)%bld(info,amold=amold,vmold=vmold,imold=imold,ilv=i) 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),& & 'Exiting with',iszv,' levels' +9998 continue call psb_erractionrestore(err_act) return diff --git a/amgprec/impl/amg_smlprec_aply.f90 b/amgprec/impl/amg_smlprec_aply.f90 index fd3f55af..d2ab3fed 100644 --- a/amgprec/impl/amg_smlprec_aply.f90 +++ b/amgprec/impl/amg_smlprec_aply.f90 @@ -393,38 +393,39 @@ contains if(debug_level > 1) then write(debug_unit,*) me,' Start inner_ml_aply at level ',level, info end if - - select case(p%precv(level)%parms%ml_cycle) - - case(amg_no_ml_) - ! - ! No preconditioning, should not really get here - ! - call psb_errpush(psb_err_internal_error_,name,& - & a_err='amg_no_ml_ in mlprc_aply?') - goto 9999 - - case(amg_add_ml_) - - call amg_s_inner_add(p, level, trans, work) - - case(amg_mult_ml_,amg_vcycle_ml_, amg_wcycle_ml_) - - call amg_s_inner_mult(p, level, trans, work) - - case(amg_kcycle_ml_, amg_kcyclesym_ml_) - - call amg_s_inner_k_cycle(p, level, trans, work) - - case default - info = psb_err_from_subroutine_ai_ - call psb_errpush(info,name,a_err='invalid ml_cycle',& - & i_Err=(/p%precv(level)%parms%ml_cycle,izero,izero,izero,izero/)) - goto 9999 - - end select - if(debug_level > 1) then - write(debug_unit,*) me,' End inner_ml_aply at level ',level + if (me >= 0) then + select case(p%precv(level)%parms%ml_cycle) + + case(amg_no_ml_) + ! + ! No preconditioning, should not really get here + ! + call psb_errpush(psb_err_internal_error_,name,& + & a_err='amg_no_ml_ in mlprc_aply?') + goto 9999 + + case(amg_add_ml_) + + call amg_s_inner_add(p, level, trans, work) + + case(amg_mult_ml_,amg_vcycle_ml_, amg_wcycle_ml_) + + call amg_s_inner_mult(p, level, trans, work) + + case(amg_kcycle_ml_, amg_kcyclesym_ml_) + + call amg_s_inner_k_cycle(p, level, trans, work) + + case default + info = psb_err_from_subroutine_ai_ + call psb_errpush(info,name,a_err='invalid ml_cycle',& + & i_Err=(/p%precv(level)%parms%ml_cycle,izero,izero,izero,izero/)) + goto 9999 + + end select + if(debug_level > 1) then + write(debug_unit,*) me,' End inner_ml_aply at level ',level + end if end if call psb_erractionrestore(err_act) @@ -492,28 +493,30 @@ contains & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & wv => p%precv(level)%wrk%wv) - 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) - do k=1, sweeps + if (me >= 0) then + 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) + do k=1, sweeps + call p%precv(level)%sm%apply(sone,& + & vy2l,szero,vty,& + & base_desc, trans,& + & ione,work,wv,info,init='Z') + + call p%precv(level)%sm2a%apply(sone,& + & vty,szero,vy2l,& + & base_desc, trans,& + & ione,work,wv,info,init='Z') + end do + + else + sweeps = p%precv(level)%parms%sweeps_pre call p%precv(level)%sm%apply(sone,& - & vy2l,szero,vty,& - & base_desc, trans,& - & ione,work,wv,info,init='Z') - - call p%precv(level)%sm2a%apply(sone,& - & vty,szero,vy2l,& + & vx2l,szero,vy2l,& & base_desc, trans,& - & ione,work,wv,info,init='Z') - end do - - else - sweeps = p%precv(level)%parms%sweeps_pre - call p%precv(level)%sm%apply(sone,& - & vx2l,szero,vy2l,& - & base_desc, trans,& - & sweeps,work,wv,info,init='Z') + & sweeps,work,wv,info,init='Z') + end if end if if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -602,7 +605,6 @@ contains end if ctxt = p%precv(level)%base_desc%get_context() call psb_info(ctxt, me, np) - if(debug_level > 1) then write(debug_unit,*) me,' inner_mult at level ',level end if @@ -623,22 +625,25 @@ contains ! if (pre) 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 PRE smoother_apply') - goto 9999 + if (me >=0) then +!!$ write(0,*) me,'Applying smoother pre ', level + 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 PRE smoother_apply') + goto 9999 + end if end if endif ! @@ -696,11 +701,13 @@ contains if (p%precv(level)%parms%ml_cycle == amg_wcycle_ml_) then - call psb_geaxpby(sone,vx2l, szero,vty,& - & base_desc,info) - if (info == psb_success_) call psb_spmm(-sone,base_a,& - & vy2l,sone,vty,& - & base_desc,info,work=work,trans=trans) + if (me >=0) then + call psb_geaxpby(sone,vx2l, szero,vty,& + & base_desc,info) + if (info == psb_success_) call psb_spmm(-sone,base_a,& + & vy2l,sone,vty,& + & base_desc,info,work=work,trans=trans) + end if if (info == psb_success_) & & call p%precv(level+1)%map_rstr(sone,vty,& & szero,p%precv(level+1)%wrk%vx2l,info,work=work,& @@ -728,31 +735,33 @@ contains if (post) then - call psb_geaxpby(sone,vx2l,& - & szero,vty,& - & base_desc,info) - if (info == psb_success_) 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 second 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') + if (me >=0) then + call psb_geaxpby(sone,vx2l,& + & szero,vty,& + & base_desc,info) + if (info == psb_success_) 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 second 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 @@ -764,11 +773,14 @@ contains endif else if (level == nlev) 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) +!!$ write(0,*) me,'Applying smoother at top level ',psb_errstatus_fatal() + if (me >=0) 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) + end if +!!$ write(0,*) me,' Done applying smoother at top level ',psb_errstatus_fatal() else @@ -778,7 +790,7 @@ contains goto 9999 end if end associate - +9998 continue call psb_erractionrestore(err_act) return @@ -829,70 +841,71 @@ contains end if ctxt = p%precv(level)%base_desc%get_context() call psb_info(ctxt, me, np) - + if(debug_level > 1) then write(debug_unit,*) me,name,' start at level ',level end if - + if ((level<1).or.(level>nlev)) then info = psb_err_internal_error_ call psb_errpush(info,name,& & a_err='Invalid LEVEL>NLEV') goto 9999 end if - + !K cycle associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& & 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,& & 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 (trans == 'N') then + if (level == nlev) then + if (me >= 0) 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 - 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 + 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,& - & szero,vty,& - & base_desc,info) + call psb_geaxpby(sone,vx2l,& + & szero,vty,& + & base_desc,info) - if (info == psb_success_) 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 + if (info == psb_success_) 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 end if - ! Apply the restriction call p%precv(level + 1)%map_rstr(sone,vty,& & szero,p%precv(level + 1)%wrk%vx2l,& @@ -940,40 +953,42 @@ contains & a_err='Error during prolongation') goto 9999 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 - ! - ! 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') + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during POST smoother_apply') + goto 9999 + 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 info = psb_err_internal_error_ @@ -990,8 +1005,7 @@ contains return end subroutine amg_s_inner_k_cycle - - + recursive subroutine amg_sinneritkcycle(p, level, trans, work, innersolv) use psb_base_mod use amg_prec_mod diff --git a/amgprec/impl/amg_z_hierarchy_bld.f90 b/amgprec/impl/amg_z_hierarchy_bld.f90 index 15f4d48a..173c515e 100644 --- a/amgprec/impl/amg_z_hierarchy_bld.f90 +++ b/amgprec/impl/amg_z_hierarchy_bld.f90 @@ -455,6 +455,31 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info) end do 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 call psb_errpush(psb_err_internal_error_,name,& & a_err='Internal hierarchy build' ) diff --git a/amgprec/impl/amg_z_smoothers_bld.f90 b/amgprec/impl/amg_z_smoothers_bld.f90 index d2cdd7d9..73836341 100644 --- a/amgprec/impl/amg_z_smoothers_bld.f90 +++ b/amgprec/impl/amg_z_smoothers_bld.f90 @@ -117,6 +117,10 @@ subroutine amg_z_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold) info = psb_success_ ctxt = desc_a%get_context() 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_) & & 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 ! +!!$ write(0,*) me,' Building at level ',i call prec%precv(i)%bld(info,amold=amold,vmold=vmold,imold=imold,ilv=i) 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),& & 'Exiting with',iszv,' levels' +9998 continue call psb_erractionrestore(err_act) return diff --git a/amgprec/impl/amg_zmlprec_aply.f90 b/amgprec/impl/amg_zmlprec_aply.f90 index 46dba71a..eca31ebd 100644 --- a/amgprec/impl/amg_zmlprec_aply.f90 +++ b/amgprec/impl/amg_zmlprec_aply.f90 @@ -393,38 +393,39 @@ contains if(debug_level > 1) then write(debug_unit,*) me,' Start inner_ml_aply at level ',level, info end if - - select case(p%precv(level)%parms%ml_cycle) - - case(amg_no_ml_) - ! - ! No preconditioning, should not really get here - ! - call psb_errpush(psb_err_internal_error_,name,& - & a_err='amg_no_ml_ in mlprc_aply?') - goto 9999 - - case(amg_add_ml_) - - call amg_z_inner_add(p, level, trans, work) - - case(amg_mult_ml_,amg_vcycle_ml_, amg_wcycle_ml_) - - call amg_z_inner_mult(p, level, trans, work) - - case(amg_kcycle_ml_, amg_kcyclesym_ml_) - - call amg_z_inner_k_cycle(p, level, trans, work) - - case default - info = psb_err_from_subroutine_ai_ - call psb_errpush(info,name,a_err='invalid ml_cycle',& - & i_Err=(/p%precv(level)%parms%ml_cycle,izero,izero,izero,izero/)) - goto 9999 - - end select - if(debug_level > 1) then - write(debug_unit,*) me,' End inner_ml_aply at level ',level + if (me >= 0) then + select case(p%precv(level)%parms%ml_cycle) + + case(amg_no_ml_) + ! + ! No preconditioning, should not really get here + ! + call psb_errpush(psb_err_internal_error_,name,& + & a_err='amg_no_ml_ in mlprc_aply?') + goto 9999 + + case(amg_add_ml_) + + call amg_z_inner_add(p, level, trans, work) + + case(amg_mult_ml_,amg_vcycle_ml_, amg_wcycle_ml_) + + call amg_z_inner_mult(p, level, trans, work) + + case(amg_kcycle_ml_, amg_kcyclesym_ml_) + + call amg_z_inner_k_cycle(p, level, trans, work) + + case default + info = psb_err_from_subroutine_ai_ + call psb_errpush(info,name,a_err='invalid ml_cycle',& + & i_Err=(/p%precv(level)%parms%ml_cycle,izero,izero,izero,izero/)) + goto 9999 + + end select + if(debug_level > 1) then + write(debug_unit,*) me,' End inner_ml_aply at level ',level + end if end if call psb_erractionrestore(err_act) @@ -492,28 +493,30 @@ contains & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & wv => p%precv(level)%wrk%wv) - 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) - do k=1, sweeps + if (me >= 0) then + 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) + do k=1, sweeps + call p%precv(level)%sm%apply(zone,& + & vy2l,zzero,vty,& + & base_desc, trans,& + & ione,work,wv,info,init='Z') + + call p%precv(level)%sm2a%apply(zone,& + & vty,zzero,vy2l,& + & base_desc, trans,& + & ione,work,wv,info,init='Z') + end do + + else + sweeps = p%precv(level)%parms%sweeps_pre call p%precv(level)%sm%apply(zone,& - & vy2l,zzero,vty,& - & base_desc, trans,& - & ione,work,wv,info,init='Z') - - call p%precv(level)%sm2a%apply(zone,& - & vty,zzero,vy2l,& + & vx2l,zzero,vy2l,& & base_desc, trans,& - & ione,work,wv,info,init='Z') - end do - - else - sweeps = p%precv(level)%parms%sweeps_pre - call p%precv(level)%sm%apply(zone,& - & vx2l,zzero,vy2l,& - & base_desc, trans,& - & sweeps,work,wv,info,init='Z') + & sweeps,work,wv,info,init='Z') + end if end if if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -602,7 +605,6 @@ contains end if ctxt = p%precv(level)%base_desc%get_context() call psb_info(ctxt, me, np) - if(debug_level > 1) then write(debug_unit,*) me,' inner_mult at level ',level end if @@ -623,22 +625,25 @@ contains ! if (pre) 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 PRE smoother_apply') - goto 9999 + if (me >=0) then +!!$ write(0,*) me,'Applying smoother pre ', level + 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 PRE smoother_apply') + goto 9999 + end if end if endif ! @@ -696,11 +701,13 @@ contains if (p%precv(level)%parms%ml_cycle == amg_wcycle_ml_) then - call psb_geaxpby(zone,vx2l, zzero,vty,& - & base_desc,info) - if (info == psb_success_) call psb_spmm(-zone,base_a,& - & vy2l,zone,vty,& - & base_desc,info,work=work,trans=trans) + if (me >=0) then + call psb_geaxpby(zone,vx2l, zzero,vty,& + & base_desc,info) + if (info == psb_success_) call psb_spmm(-zone,base_a,& + & vy2l,zone,vty,& + & base_desc,info,work=work,trans=trans) + end if if (info == psb_success_) & & call p%precv(level+1)%map_rstr(zone,vty,& & zzero,p%precv(level+1)%wrk%vx2l,info,work=work,& @@ -728,31 +735,33 @@ contains if (post) then - call psb_geaxpby(zone,vx2l,& - & zzero,vty,& - & base_desc,info) - if (info == psb_success_) 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 second 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') + if (me >=0) then + call psb_geaxpby(zone,vx2l,& + & zzero,vty,& + & base_desc,info) + if (info == psb_success_) 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 second 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 @@ -764,11 +773,14 @@ contains endif else if (level == nlev) 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) +!!$ write(0,*) me,'Applying smoother at top level ',psb_errstatus_fatal() + if (me >=0) 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) + end if +!!$ write(0,*) me,' Done applying smoother at top level ',psb_errstatus_fatal() else @@ -778,7 +790,7 @@ contains goto 9999 end if end associate - +9998 continue call psb_erractionrestore(err_act) return @@ -829,70 +841,71 @@ contains end if ctxt = p%precv(level)%base_desc%get_context() call psb_info(ctxt, me, np) - + if(debug_level > 1) then write(debug_unit,*) me,name,' start at level ',level end if - + if ((level<1).or.(level>nlev)) then info = psb_err_internal_error_ call psb_errpush(info,name,& & a_err='Invalid LEVEL>NLEV') goto 9999 end if - + !K cycle associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& & 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,& & 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 (trans == 'N') then + if (level == nlev) then + if (me >= 0) 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 - 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 + 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,& - & zzero,vty,& - & base_desc,info) + call psb_geaxpby(zone,vx2l,& + & zzero,vty,& + & base_desc,info) - if (info == psb_success_) 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 + if (info == psb_success_) 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 end if - ! Apply the restriction call p%precv(level + 1)%map_rstr(zone,vty,& & zzero,p%precv(level + 1)%wrk%vx2l,& @@ -940,40 +953,42 @@ contains & a_err='Error during prolongation') goto 9999 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 - ! - ! 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') + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during POST smoother_apply') + goto 9999 + 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 info = psb_err_internal_error_ @@ -990,8 +1005,7 @@ contains return end subroutine amg_z_inner_k_cycle - - + recursive subroutine amg_zinneritkcycle(p, level, trans, work, innersolv) use psb_base_mod use amg_prec_mod diff --git a/amgprec/impl/level/amg_c_base_onelev_build.f90 b/amgprec/impl/level/amg_c_base_onelev_build.f90 index 9af46a8f..444bcdd0 100644 --- a/amgprec/impl/level/amg_c_base_onelev_build.f90 +++ b/amgprec/impl/level/amg_c_base_onelev_build.f90 @@ -71,72 +71,81 @@ subroutine amg_c_base_onelev_build(lv,info,amold,vmold,imold,ilv) ctxt = lv%base_desc%get_ctxt() call psb_info(ctxt,me,np) - if (.not.allocated(lv%sm)) then - !! Error: should have called amg_dprecinit - info=3111 - call psb_errpush(info,name) - goto 9999 - end if - if (.not.allocated(lv%sm%sv)) then - !! Error: should have called amg_dprecinit - info=3111 - call psb_errpush(info,name) - goto 9999 - end if - lv%ac_nz_loc = lv%ac%get_nzeros() - lv%ac_nz_tot = lv%ac_nz_loc - select case(lv%parms%coarse_mat) - case(amg_distr_mat_) - 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_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Calling mlprcbld at level ',i - call amg_check_def(lv%parms%sweeps_pre,& - & 'Jacobi sweeps',izero,is_int_non_negative) - call amg_check_def(lv%parms%sweeps_post,& - & 'Jacobi sweeps',izero,is_int_non_negative) - - call lv%sm%build(lv%base_a,lv%base_desc,info) - if (info == 0) then - if (allocated(lv%sm2a)) then - call lv%sm2a%build(lv%base_a,lv%base_desc,info) - lv%sm2 => lv%sm2a - else - lv%sm2 => lv%sm + ! + ! At top level(s) I may be using + ! a context with less processes + ! + if (me < 0) then +!!$ write(0,*) 'onelevbld: I am excluded from this one ' + else +!!$ write(0,*) me,' Going to build smoothers at this level ' + if (.not.allocated(lv%sm)) then + !! Error: should have called amg_dprecinit + info=3111 + call psb_errpush(info,name) + goto 9999 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%parms%sweeps_pre>1).or.(lv%parms%sweeps_post>1)) then - lv%parms%sweeps_pre = 1 - lv%parms%sweeps_post = 1 - if (me == 0) then - write(debug_unit,*) - if (present(ilv)) then - write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),& - & '" at level ',ilv - write(debug_unit,*) ' is configured as a global solver ' - else - write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),& - & '" is configured as a global solver ' + if (.not.allocated(lv%sm%sv)) then + !! Error: should have called amg_dprecinit + info=3111 + call psb_errpush(info,name) + goto 9999 + end if + lv%ac_nz_loc = lv%ac%get_nzeros() + lv%ac_nz_tot = lv%ac_nz_loc + select case(lv%parms%coarse_mat) + case(amg_distr_mat_) + 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_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Calling mlprcbld at level ',i + call amg_check_def(lv%parms%sweeps_pre,& + & 'Jacobi sweeps',izero,is_int_non_negative) + call amg_check_def(lv%parms%sweeps_post,& + & 'Jacobi sweeps',izero,is_int_non_negative) + + call lv%sm%build(lv%base_a,lv%base_desc,info) + if (info == 0) then + if (allocated(lv%sm2a)) then + call lv%sm2a%build(lv%base_a,lv%base_desc,info) + lv%sm2 => lv%sm2a + else + 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 + + if (lv%sm%sv%is_global()) then + if ((lv%parms%sweeps_pre>1).or.(lv%parms%sweeps_post>1)) then + lv%parms%sweeps_pre = 1 + lv%parms%sweeps_post = 1 + if (me == 0) then + write(debug_unit,*) + if (present(ilv)) then + write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),& + & '" at level ',ilv + write(debug_unit,*) ' is configured as a global solver ' + else + write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),& + & '" is configured as a global solver ' + end if + write(debug_unit,*) ' Pre and post sweeps at this level reset to 1' end if - write(debug_unit,*) ' Pre and post sweeps at this level reset to 1' end if end if end if diff --git a/amgprec/impl/level/amg_c_base_onelev_map_prol.F90 b/amgprec/impl/level/amg_c_base_onelev_map_prol.F90 index 0ec5dda1..82b4beb7 100644 --- a/amgprec/impl/level/amg_c_base_onelev_map_prol.F90 +++ b/amgprec/impl/level/amg_c_base_onelev_map_prol.F90 @@ -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(:) 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 ! ! 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 ! Default transfer call lv%linmap%map_V2U(alpha,vect_v,beta,vect_u,info,& diff --git a/amgprec/impl/level/amg_c_base_onelev_map_rstr.F90 b/amgprec/impl/level/amg_c_base_onelev_map_rstr.F90 index 3f714e57..6d3c0598 100644 --- a/amgprec/impl/level/amg_c_base_onelev_map_rstr.F90 +++ b/amgprec/impl/level/amg_c_base_onelev_map_rstr.F90 @@ -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 amg_c_onelev_mod, amg_protect_name => amg_c_base_onelev_map_rstr_v 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(:) 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 ! ! Remap has happened, deal with it ! - write(0,*) 'Remap handling not implemented yet ' +!!$ write(0,*) 'Remap handling not implemented yet ' 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 + complex(psb_spk_), allocatable :: rsnd(:), rrcv(:) + type(psb_c_vect_type) :: tv - ictxt = lv%remap_data%desc_ac_pre_remap%get_ctxt() - call psb_info(ictxt,me,np) + ctxt = lv%remap_data%desc_ac_pre_remap%get_ctxt() + call psb_info(ctxt,me,np) nctxt = lv%desc_ac%get_ctxt() call psb_info(nctxt,rme,rnp) +!!$ write(0,*) 'New context ',rme,rnp idest = lv%remap_data%idest associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc) - write(0,*) 'Should apply maps, then send data from ',me,' to ',idest - if (rme >= 0) write(0,*) rme, ' Receiving data from ',isrc(:) +!!$ write(0,*) 'Should apply maps, then send data from ',me,' to ',idest +!!$ 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 +!!$ write(0,*) me, ' Restrictor with remap done ' end block else diff --git a/amgprec/impl/level/amg_d_base_onelev_build.f90 b/amgprec/impl/level/amg_d_base_onelev_build.f90 index 345fbfd4..24fcff89 100644 --- a/amgprec/impl/level/amg_d_base_onelev_build.f90 +++ b/amgprec/impl/level/amg_d_base_onelev_build.f90 @@ -71,72 +71,81 @@ subroutine amg_d_base_onelev_build(lv,info,amold,vmold,imold,ilv) ctxt = lv%base_desc%get_ctxt() call psb_info(ctxt,me,np) - if (.not.allocated(lv%sm)) then - !! Error: should have called amg_dprecinit - info=3111 - call psb_errpush(info,name) - goto 9999 - end if - if (.not.allocated(lv%sm%sv)) then - !! Error: should have called amg_dprecinit - info=3111 - call psb_errpush(info,name) - goto 9999 - end if - lv%ac_nz_loc = lv%ac%get_nzeros() - lv%ac_nz_tot = lv%ac_nz_loc - select case(lv%parms%coarse_mat) - case(amg_distr_mat_) - 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_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Calling mlprcbld at level ',i - call amg_check_def(lv%parms%sweeps_pre,& - & 'Jacobi sweeps',izero,is_int_non_negative) - call amg_check_def(lv%parms%sweeps_post,& - & 'Jacobi sweeps',izero,is_int_non_negative) - - call lv%sm%build(lv%base_a,lv%base_desc,info) - if (info == 0) then - if (allocated(lv%sm2a)) then - call lv%sm2a%build(lv%base_a,lv%base_desc,info) - lv%sm2 => lv%sm2a - else - lv%sm2 => lv%sm + ! + ! At top level(s) I may be using + ! a context with less processes + ! + if (me < 0) then +!!$ write(0,*) 'onelevbld: I am excluded from this one ' + else +!!$ write(0,*) me,' Going to build smoothers at this level ' + if (.not.allocated(lv%sm)) then + !! Error: should have called amg_dprecinit + info=3111 + call psb_errpush(info,name) + goto 9999 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%parms%sweeps_pre>1).or.(lv%parms%sweeps_post>1)) then - lv%parms%sweeps_pre = 1 - lv%parms%sweeps_post = 1 - if (me == 0) then - write(debug_unit,*) - if (present(ilv)) then - write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),& - & '" at level ',ilv - write(debug_unit,*) ' is configured as a global solver ' - else - write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),& - & '" is configured as a global solver ' + if (.not.allocated(lv%sm%sv)) then + !! Error: should have called amg_dprecinit + info=3111 + call psb_errpush(info,name) + goto 9999 + end if + lv%ac_nz_loc = lv%ac%get_nzeros() + lv%ac_nz_tot = lv%ac_nz_loc + select case(lv%parms%coarse_mat) + case(amg_distr_mat_) + 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_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Calling mlprcbld at level ',i + call amg_check_def(lv%parms%sweeps_pre,& + & 'Jacobi sweeps',izero,is_int_non_negative) + call amg_check_def(lv%parms%sweeps_post,& + & 'Jacobi sweeps',izero,is_int_non_negative) + + call lv%sm%build(lv%base_a,lv%base_desc,info) + if (info == 0) then + if (allocated(lv%sm2a)) then + call lv%sm2a%build(lv%base_a,lv%base_desc,info) + lv%sm2 => lv%sm2a + else + 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 + + if (lv%sm%sv%is_global()) then + if ((lv%parms%sweeps_pre>1).or.(lv%parms%sweeps_post>1)) then + lv%parms%sweeps_pre = 1 + lv%parms%sweeps_post = 1 + if (me == 0) then + write(debug_unit,*) + if (present(ilv)) then + write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),& + & '" at level ',ilv + write(debug_unit,*) ' is configured as a global solver ' + else + write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),& + & '" is configured as a global solver ' + end if + write(debug_unit,*) ' Pre and post sweeps at this level reset to 1' end if - write(debug_unit,*) ' Pre and post sweeps at this level reset to 1' end if end if end if diff --git a/amgprec/impl/level/amg_d_base_onelev_map_prol.F90 b/amgprec/impl/level/amg_d_base_onelev_map_prol.F90 index 4d7e27e0..ed33a38e 100644 --- a/amgprec/impl/level/amg_d_base_onelev_map_prol.F90 +++ b/amgprec/impl/level/amg_d_base_onelev_map_prol.F90 @@ -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(:) 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 ! ! 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 ! Default transfer call lv%linmap%map_V2U(alpha,vect_v,beta,vect_u,info,& diff --git a/amgprec/impl/level/amg_d_base_onelev_map_rstr.F90 b/amgprec/impl/level/amg_d_base_onelev_map_rstr.F90 index 455bdac0..9318fc86 100644 --- a/amgprec/impl/level/amg_d_base_onelev_map_rstr.F90 +++ b/amgprec/impl/level/amg_d_base_onelev_map_rstr.F90 @@ -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 amg_d_onelev_mod, amg_protect_name => amg_d_base_onelev_map_rstr_v 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(:) 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 ! ! Remap has happened, deal with it ! - write(0,*) 'Remap handling not implemented yet ' +!!$ write(0,*) 'Remap handling not implemented yet ' 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 + real(psb_dpk_), allocatable :: rsnd(:), rrcv(:) + type(psb_d_vect_type) :: tv - ictxt = lv%remap_data%desc_ac_pre_remap%get_ctxt() - call psb_info(ictxt,me,np) + ctxt = lv%remap_data%desc_ac_pre_remap%get_ctxt() + call psb_info(ctxt,me,np) nctxt = lv%desc_ac%get_ctxt() call psb_info(nctxt,rme,rnp) +!!$ write(0,*) 'New context ',rme,rnp idest = lv%remap_data%idest associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc) - write(0,*) 'Should apply maps, then send data from ',me,' to ',idest - if (rme >= 0) write(0,*) rme, ' Receiving data from ',isrc(:) +!!$ write(0,*) 'Should apply maps, then send data from ',me,' to ',idest +!!$ 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 +!!$ write(0,*) me, ' Restrictor with remap done ' end block else diff --git a/amgprec/impl/level/amg_s_base_onelev_build.f90 b/amgprec/impl/level/amg_s_base_onelev_build.f90 index 4a2e207a..4b80fb5d 100644 --- a/amgprec/impl/level/amg_s_base_onelev_build.f90 +++ b/amgprec/impl/level/amg_s_base_onelev_build.f90 @@ -71,72 +71,81 @@ subroutine amg_s_base_onelev_build(lv,info,amold,vmold,imold,ilv) ctxt = lv%base_desc%get_ctxt() call psb_info(ctxt,me,np) - if (.not.allocated(lv%sm)) then - !! Error: should have called amg_dprecinit - info=3111 - call psb_errpush(info,name) - goto 9999 - end if - if (.not.allocated(lv%sm%sv)) then - !! Error: should have called amg_dprecinit - info=3111 - call psb_errpush(info,name) - goto 9999 - end if - lv%ac_nz_loc = lv%ac%get_nzeros() - lv%ac_nz_tot = lv%ac_nz_loc - select case(lv%parms%coarse_mat) - case(amg_distr_mat_) - 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_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Calling mlprcbld at level ',i - call amg_check_def(lv%parms%sweeps_pre,& - & 'Jacobi sweeps',izero,is_int_non_negative) - call amg_check_def(lv%parms%sweeps_post,& - & 'Jacobi sweeps',izero,is_int_non_negative) - - call lv%sm%build(lv%base_a,lv%base_desc,info) - if (info == 0) then - if (allocated(lv%sm2a)) then - call lv%sm2a%build(lv%base_a,lv%base_desc,info) - lv%sm2 => lv%sm2a - else - lv%sm2 => lv%sm + ! + ! At top level(s) I may be using + ! a context with less processes + ! + if (me < 0) then +!!$ write(0,*) 'onelevbld: I am excluded from this one ' + else +!!$ write(0,*) me,' Going to build smoothers at this level ' + if (.not.allocated(lv%sm)) then + !! Error: should have called amg_dprecinit + info=3111 + call psb_errpush(info,name) + goto 9999 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%parms%sweeps_pre>1).or.(lv%parms%sweeps_post>1)) then - lv%parms%sweeps_pre = 1 - lv%parms%sweeps_post = 1 - if (me == 0) then - write(debug_unit,*) - if (present(ilv)) then - write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),& - & '" at level ',ilv - write(debug_unit,*) ' is configured as a global solver ' - else - write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),& - & '" is configured as a global solver ' + if (.not.allocated(lv%sm%sv)) then + !! Error: should have called amg_dprecinit + info=3111 + call psb_errpush(info,name) + goto 9999 + end if + lv%ac_nz_loc = lv%ac%get_nzeros() + lv%ac_nz_tot = lv%ac_nz_loc + select case(lv%parms%coarse_mat) + case(amg_distr_mat_) + 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_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Calling mlprcbld at level ',i + call amg_check_def(lv%parms%sweeps_pre,& + & 'Jacobi sweeps',izero,is_int_non_negative) + call amg_check_def(lv%parms%sweeps_post,& + & 'Jacobi sweeps',izero,is_int_non_negative) + + call lv%sm%build(lv%base_a,lv%base_desc,info) + if (info == 0) then + if (allocated(lv%sm2a)) then + call lv%sm2a%build(lv%base_a,lv%base_desc,info) + lv%sm2 => lv%sm2a + else + 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 + + if (lv%sm%sv%is_global()) then + if ((lv%parms%sweeps_pre>1).or.(lv%parms%sweeps_post>1)) then + lv%parms%sweeps_pre = 1 + lv%parms%sweeps_post = 1 + if (me == 0) then + write(debug_unit,*) + if (present(ilv)) then + write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),& + & '" at level ',ilv + write(debug_unit,*) ' is configured as a global solver ' + else + write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),& + & '" is configured as a global solver ' + end if + write(debug_unit,*) ' Pre and post sweeps at this level reset to 1' end if - write(debug_unit,*) ' Pre and post sweeps at this level reset to 1' end if end if end if diff --git a/amgprec/impl/level/amg_s_base_onelev_map_prol.F90 b/amgprec/impl/level/amg_s_base_onelev_map_prol.F90 index 1b1777b4..81a6334f 100644 --- a/amgprec/impl/level/amg_s_base_onelev_map_prol.F90 +++ b/amgprec/impl/level/amg_s_base_onelev_map_prol.F90 @@ -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(:) 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 ! ! 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 ! Default transfer call lv%linmap%map_V2U(alpha,vect_v,beta,vect_u,info,& diff --git a/amgprec/impl/level/amg_s_base_onelev_map_rstr.F90 b/amgprec/impl/level/amg_s_base_onelev_map_rstr.F90 index fad8ceac..0ccd31d7 100644 --- a/amgprec/impl/level/amg_s_base_onelev_map_rstr.F90 +++ b/amgprec/impl/level/amg_s_base_onelev_map_rstr.F90 @@ -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 amg_s_onelev_mod, amg_protect_name => amg_s_base_onelev_map_rstr_v 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(:) 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 ! ! Remap has happened, deal with it ! - write(0,*) 'Remap handling not implemented yet ' +!!$ write(0,*) 'Remap handling not implemented yet ' 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 + real(psb_spk_), allocatable :: rsnd(:), rrcv(:) + type(psb_s_vect_type) :: tv - ictxt = lv%remap_data%desc_ac_pre_remap%get_ctxt() - call psb_info(ictxt,me,np) + ctxt = lv%remap_data%desc_ac_pre_remap%get_ctxt() + call psb_info(ctxt,me,np) nctxt = lv%desc_ac%get_ctxt() call psb_info(nctxt,rme,rnp) +!!$ write(0,*) 'New context ',rme,rnp idest = lv%remap_data%idest associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc) - write(0,*) 'Should apply maps, then send data from ',me,' to ',idest - if (rme >= 0) write(0,*) rme, ' Receiving data from ',isrc(:) +!!$ write(0,*) 'Should apply maps, then send data from ',me,' to ',idest +!!$ 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 +!!$ write(0,*) me, ' Restrictor with remap done ' end block else diff --git a/amgprec/impl/level/amg_z_base_onelev_build.f90 b/amgprec/impl/level/amg_z_base_onelev_build.f90 index 57020edf..b497ef1c 100644 --- a/amgprec/impl/level/amg_z_base_onelev_build.f90 +++ b/amgprec/impl/level/amg_z_base_onelev_build.f90 @@ -71,72 +71,81 @@ subroutine amg_z_base_onelev_build(lv,info,amold,vmold,imold,ilv) ctxt = lv%base_desc%get_ctxt() call psb_info(ctxt,me,np) - if (.not.allocated(lv%sm)) then - !! Error: should have called amg_dprecinit - info=3111 - call psb_errpush(info,name) - goto 9999 - end if - if (.not.allocated(lv%sm%sv)) then - !! Error: should have called amg_dprecinit - info=3111 - call psb_errpush(info,name) - goto 9999 - end if - lv%ac_nz_loc = lv%ac%get_nzeros() - lv%ac_nz_tot = lv%ac_nz_loc - select case(lv%parms%coarse_mat) - case(amg_distr_mat_) - 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_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Calling mlprcbld at level ',i - call amg_check_def(lv%parms%sweeps_pre,& - & 'Jacobi sweeps',izero,is_int_non_negative) - call amg_check_def(lv%parms%sweeps_post,& - & 'Jacobi sweeps',izero,is_int_non_negative) - - call lv%sm%build(lv%base_a,lv%base_desc,info) - if (info == 0) then - if (allocated(lv%sm2a)) then - call lv%sm2a%build(lv%base_a,lv%base_desc,info) - lv%sm2 => lv%sm2a - else - lv%sm2 => lv%sm + ! + ! At top level(s) I may be using + ! a context with less processes + ! + if (me < 0) then +!!$ write(0,*) 'onelevbld: I am excluded from this one ' + else +!!$ write(0,*) me,' Going to build smoothers at this level ' + if (.not.allocated(lv%sm)) then + !! Error: should have called amg_dprecinit + info=3111 + call psb_errpush(info,name) + goto 9999 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%parms%sweeps_pre>1).or.(lv%parms%sweeps_post>1)) then - lv%parms%sweeps_pre = 1 - lv%parms%sweeps_post = 1 - if (me == 0) then - write(debug_unit,*) - if (present(ilv)) then - write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),& - & '" at level ',ilv - write(debug_unit,*) ' is configured as a global solver ' - else - write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),& - & '" is configured as a global solver ' + if (.not.allocated(lv%sm%sv)) then + !! Error: should have called amg_dprecinit + info=3111 + call psb_errpush(info,name) + goto 9999 + end if + lv%ac_nz_loc = lv%ac%get_nzeros() + lv%ac_nz_tot = lv%ac_nz_loc + select case(lv%parms%coarse_mat) + case(amg_distr_mat_) + 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_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Calling mlprcbld at level ',i + call amg_check_def(lv%parms%sweeps_pre,& + & 'Jacobi sweeps',izero,is_int_non_negative) + call amg_check_def(lv%parms%sweeps_post,& + & 'Jacobi sweeps',izero,is_int_non_negative) + + call lv%sm%build(lv%base_a,lv%base_desc,info) + if (info == 0) then + if (allocated(lv%sm2a)) then + call lv%sm2a%build(lv%base_a,lv%base_desc,info) + lv%sm2 => lv%sm2a + else + 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 + + if (lv%sm%sv%is_global()) then + if ((lv%parms%sweeps_pre>1).or.(lv%parms%sweeps_post>1)) then + lv%parms%sweeps_pre = 1 + lv%parms%sweeps_post = 1 + if (me == 0) then + write(debug_unit,*) + if (present(ilv)) then + write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),& + & '" at level ',ilv + write(debug_unit,*) ' is configured as a global solver ' + else + write(debug_unit,*) 'Warning: the solver "',trim(lv%sm%sv%get_fmt()),& + & '" is configured as a global solver ' + end if + write(debug_unit,*) ' Pre and post sweeps at this level reset to 1' end if - write(debug_unit,*) ' Pre and post sweeps at this level reset to 1' end if end if end if diff --git a/amgprec/impl/level/amg_z_base_onelev_map_prol.F90 b/amgprec/impl/level/amg_z_base_onelev_map_prol.F90 index 16c982cc..1fe59b44 100644 --- a/amgprec/impl/level/amg_z_base_onelev_map_prol.F90 +++ b/amgprec/impl/level/amg_z_base_onelev_map_prol.F90 @@ -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(:) 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 ! ! 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 ! Default transfer call lv%linmap%map_V2U(alpha,vect_v,beta,vect_u,info,& diff --git a/amgprec/impl/level/amg_z_base_onelev_map_rstr.F90 b/amgprec/impl/level/amg_z_base_onelev_map_rstr.F90 index d857e5b2..a71c5a54 100644 --- a/amgprec/impl/level/amg_z_base_onelev_map_rstr.F90 +++ b/amgprec/impl/level/amg_z_base_onelev_map_rstr.F90 @@ -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 amg_z_onelev_mod, amg_protect_name => amg_z_base_onelev_map_rstr_v 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(:) 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 ! ! Remap has happened, deal with it ! - write(0,*) 'Remap handling not implemented yet ' +!!$ write(0,*) 'Remap handling not implemented yet ' 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 + complex(psb_dpk_), allocatable :: rsnd(:), rrcv(:) + type(psb_z_vect_type) :: tv - ictxt = lv%remap_data%desc_ac_pre_remap%get_ctxt() - call psb_info(ictxt,me,np) + ctxt = lv%remap_data%desc_ac_pre_remap%get_ctxt() + call psb_info(ctxt,me,np) nctxt = lv%desc_ac%get_ctxt() call psb_info(nctxt,rme,rnp) +!!$ write(0,*) 'New context ',rme,rnp idest = lv%remap_data%idest associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc) - write(0,*) 'Should apply maps, then send data from ',me,' to ',idest - if (rme >= 0) write(0,*) rme, ' Receiving data from ',isrc(:) +!!$ write(0,*) 'Should apply maps, then send data from ',me,' to ',idest +!!$ 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 +!!$ write(0,*) me, ' Restrictor with remap done ' end block else