diff --git a/base/modules/comm/psb_c_linmap_mod.f90 b/base/modules/comm/psb_c_linmap_mod.f90 index 69e3359c..5786ecad 100644 --- a/base/modules/comm/psb_c_linmap_mod.f90 +++ b/base/modules/comm/psb_c_linmap_mod.f90 @@ -65,7 +65,7 @@ module psb_c_linmap_mod integer(psb_ipk_), intent(out) :: info complex(psb_spk_), optional :: work(:) end subroutine psb_c_map_X2Y - subroutine psb_c_map_X2Y_vect(alpha,x,beta,y,map,info,work) + subroutine psb_c_map_X2Y_vect(alpha,x,beta,y,map,info,work,vtx,vty) use psb_c_vect_mod, only : psb_c_vect_type import :: psb_ipk_, psb_spk_, psb_clinmap_type implicit none @@ -74,6 +74,7 @@ module psb_c_linmap_mod type(psb_c_vect_type), intent(inout) :: x,y integer(psb_ipk_), intent(out) :: info complex(psb_spk_), optional :: work(:) + type(psb_c_vect_type), optional, target, intent(inout) :: vtx,vty end subroutine psb_c_map_X2Y_vect end interface @@ -88,7 +89,7 @@ module psb_c_linmap_mod integer(psb_ipk_), intent(out) :: info complex(psb_spk_), optional :: work(:) end subroutine psb_c_map_Y2X - subroutine psb_c_map_Y2X_vect(alpha,x,beta,y,map,info,work) + subroutine psb_c_map_Y2X_vect(alpha,x,beta,y,map,info,work,vtx,vty) use psb_c_vect_mod, only : psb_c_vect_type import :: psb_ipk_, psb_spk_, psb_clinmap_type implicit none @@ -97,6 +98,7 @@ module psb_c_linmap_mod type(psb_c_vect_type), intent(inout) :: x,y integer(psb_ipk_), intent(out) :: info complex(psb_spk_), optional :: work(:) + type(psb_c_vect_type), optional, target, intent(inout) :: vtx,vty end subroutine psb_c_map_Y2X_vect end interface @@ -170,8 +172,9 @@ contains class(psb_c_base_sparse_mat), intent(in), optional :: mold class(psb_i_base_vect_type), intent(in), optional :: imold - call map%map_X2Y%cscnv(info,type=type,mold=mold) - if (info == psb_success_)& + if (map%map_X2Y%is_asb())& + & call map%map_X2Y%cscnv(info,type=type,mold=mold) + if (info == psb_success_ .and.map%map_Y2X%is_asb())& & call map%map_Y2X%cscnv(info,type=type,mold=mold) if (present(imold)) then call map%desc_X%cnv(mold=imold) diff --git a/base/modules/comm/psb_d_linmap_mod.f90 b/base/modules/comm/psb_d_linmap_mod.f90 index c513bfa1..1494446c 100644 --- a/base/modules/comm/psb_d_linmap_mod.f90 +++ b/base/modules/comm/psb_d_linmap_mod.f90 @@ -65,7 +65,7 @@ module psb_d_linmap_mod integer(psb_ipk_), intent(out) :: info real(psb_dpk_), optional :: work(:) end subroutine psb_d_map_X2Y - subroutine psb_d_map_X2Y_vect(alpha,x,beta,y,map,info,work) + subroutine psb_d_map_X2Y_vect(alpha,x,beta,y,map,info,work,vtx,vty) use psb_d_vect_mod, only : psb_d_vect_type import :: psb_ipk_, psb_dpk_, psb_dlinmap_type implicit none @@ -74,6 +74,7 @@ module psb_d_linmap_mod type(psb_d_vect_type), intent(inout) :: x,y integer(psb_ipk_), intent(out) :: info real(psb_dpk_), optional :: work(:) + type(psb_d_vect_type), optional, target, intent(inout) :: vtx,vty end subroutine psb_d_map_X2Y_vect end interface @@ -88,7 +89,7 @@ module psb_d_linmap_mod integer(psb_ipk_), intent(out) :: info real(psb_dpk_), optional :: work(:) end subroutine psb_d_map_Y2X - subroutine psb_d_map_Y2X_vect(alpha,x,beta,y,map,info,work) + subroutine psb_d_map_Y2X_vect(alpha,x,beta,y,map,info,work,vtx,vty) use psb_d_vect_mod, only : psb_d_vect_type import :: psb_ipk_, psb_dpk_, psb_dlinmap_type implicit none @@ -97,6 +98,7 @@ module psb_d_linmap_mod type(psb_d_vect_type), intent(inout) :: x,y integer(psb_ipk_), intent(out) :: info real(psb_dpk_), optional :: work(:) + type(psb_d_vect_type), optional, target, intent(inout) :: vtx,vty end subroutine psb_d_map_Y2X_vect end interface @@ -170,8 +172,9 @@ contains class(psb_d_base_sparse_mat), intent(in), optional :: mold class(psb_i_base_vect_type), intent(in), optional :: imold - call map%map_X2Y%cscnv(info,type=type,mold=mold) - if (info == psb_success_)& + if (map%map_X2Y%is_asb())& + & call map%map_X2Y%cscnv(info,type=type,mold=mold) + if (info == psb_success_ .and.map%map_Y2X%is_asb())& & call map%map_Y2X%cscnv(info,type=type,mold=mold) if (present(imold)) then call map%desc_X%cnv(mold=imold) diff --git a/base/modules/comm/psb_s_linmap_mod.f90 b/base/modules/comm/psb_s_linmap_mod.f90 index eb94decb..a2ef105b 100644 --- a/base/modules/comm/psb_s_linmap_mod.f90 +++ b/base/modules/comm/psb_s_linmap_mod.f90 @@ -65,7 +65,7 @@ module psb_s_linmap_mod integer(psb_ipk_), intent(out) :: info real(psb_spk_), optional :: work(:) end subroutine psb_s_map_X2Y - subroutine psb_s_map_X2Y_vect(alpha,x,beta,y,map,info,work) + subroutine psb_s_map_X2Y_vect(alpha,x,beta,y,map,info,work,vtx,vty) use psb_s_vect_mod, only : psb_s_vect_type import :: psb_ipk_, psb_spk_, psb_slinmap_type implicit none @@ -74,6 +74,7 @@ module psb_s_linmap_mod type(psb_s_vect_type), intent(inout) :: x,y integer(psb_ipk_), intent(out) :: info real(psb_spk_), optional :: work(:) + type(psb_s_vect_type), optional, target, intent(inout) :: vtx,vty end subroutine psb_s_map_X2Y_vect end interface @@ -88,7 +89,7 @@ module psb_s_linmap_mod integer(psb_ipk_), intent(out) :: info real(psb_spk_), optional :: work(:) end subroutine psb_s_map_Y2X - subroutine psb_s_map_Y2X_vect(alpha,x,beta,y,map,info,work) + subroutine psb_s_map_Y2X_vect(alpha,x,beta,y,map,info,work,vtx,vty) use psb_s_vect_mod, only : psb_s_vect_type import :: psb_ipk_, psb_spk_, psb_slinmap_type implicit none @@ -97,6 +98,7 @@ module psb_s_linmap_mod type(psb_s_vect_type), intent(inout) :: x,y integer(psb_ipk_), intent(out) :: info real(psb_spk_), optional :: work(:) + type(psb_s_vect_type), optional, target, intent(inout) :: vtx,vty end subroutine psb_s_map_Y2X_vect end interface @@ -170,8 +172,9 @@ contains class(psb_s_base_sparse_mat), intent(in), optional :: mold class(psb_i_base_vect_type), intent(in), optional :: imold - call map%map_X2Y%cscnv(info,type=type,mold=mold) - if (info == psb_success_)& + if (map%map_X2Y%is_asb())& + & call map%map_X2Y%cscnv(info,type=type,mold=mold) + if (info == psb_success_ .and.map%map_Y2X%is_asb())& & call map%map_Y2X%cscnv(info,type=type,mold=mold) if (present(imold)) then call map%desc_X%cnv(mold=imold) diff --git a/base/modules/comm/psb_z_linmap_mod.f90 b/base/modules/comm/psb_z_linmap_mod.f90 index 93a626de..d4a815d3 100644 --- a/base/modules/comm/psb_z_linmap_mod.f90 +++ b/base/modules/comm/psb_z_linmap_mod.f90 @@ -65,7 +65,7 @@ module psb_z_linmap_mod integer(psb_ipk_), intent(out) :: info complex(psb_dpk_), optional :: work(:) end subroutine psb_z_map_X2Y - subroutine psb_z_map_X2Y_vect(alpha,x,beta,y,map,info,work) + subroutine psb_z_map_X2Y_vect(alpha,x,beta,y,map,info,work,vtx,vty) use psb_z_vect_mod, only : psb_z_vect_type import :: psb_ipk_, psb_dpk_, psb_zlinmap_type implicit none @@ -74,6 +74,7 @@ module psb_z_linmap_mod type(psb_z_vect_type), intent(inout) :: x,y integer(psb_ipk_), intent(out) :: info complex(psb_dpk_), optional :: work(:) + type(psb_z_vect_type), optional, target, intent(inout) :: vtx,vty end subroutine psb_z_map_X2Y_vect end interface @@ -88,7 +89,7 @@ module psb_z_linmap_mod integer(psb_ipk_), intent(out) :: info complex(psb_dpk_), optional :: work(:) end subroutine psb_z_map_Y2X - subroutine psb_z_map_Y2X_vect(alpha,x,beta,y,map,info,work) + subroutine psb_z_map_Y2X_vect(alpha,x,beta,y,map,info,work,vtx,vty) use psb_z_vect_mod, only : psb_z_vect_type import :: psb_ipk_, psb_dpk_, psb_zlinmap_type implicit none @@ -97,6 +98,7 @@ module psb_z_linmap_mod type(psb_z_vect_type), intent(inout) :: x,y integer(psb_ipk_), intent(out) :: info complex(psb_dpk_), optional :: work(:) + type(psb_z_vect_type), optional, target, intent(inout) :: vtx,vty end subroutine psb_z_map_Y2X_vect end interface @@ -170,8 +172,9 @@ contains class(psb_z_base_sparse_mat), intent(in), optional :: mold class(psb_i_base_vect_type), intent(in), optional :: imold - call map%map_X2Y%cscnv(info,type=type,mold=mold) - if (info == psb_success_)& + if (map%map_X2Y%is_asb())& + & call map%map_X2Y%cscnv(info,type=type,mold=mold) + if (info == psb_success_ .and.map%map_Y2X%is_asb())& & call map%map_Y2X%cscnv(info,type=type,mold=mold) if (present(imold)) then call map%desc_X%cnv(mold=imold) diff --git a/base/tools/psb_c_map.f90 b/base/tools/psb_c_map.f90 index 7270e077..fbf7f057 100644 --- a/base/tools/psb_c_map.f90 +++ b/base/tools/psb_c_map.f90 @@ -111,7 +111,7 @@ subroutine psb_c_map_X2Y(alpha,x,beta,y,map,info,work) end subroutine psb_c_map_X2Y -subroutine psb_c_map_X2Y_vect(alpha,x,beta,y,map,info,work) +subroutine psb_c_map_X2Y_vect(alpha,x,beta,y,map,info,work,vtx,vty) use psb_base_mod, psb_protect_name => psb_c_map_X2Y_vect implicit none type(psb_clinmap_type), intent(in) :: map @@ -119,11 +119,13 @@ subroutine psb_c_map_X2Y_vect(alpha,x,beta,y,map,info,work) type(psb_c_vect_type), intent(inout) :: x,y integer(psb_ipk_), intent(out) :: info complex(psb_spk_), optional :: work(:) + type(psb_c_vect_type), optional, target, intent(inout) :: vtx,vty ! Local - type(psb_c_vect_type) :: xt, yt + type(psb_c_vect_type), target :: xt, yt + type(psb_c_vect_type),pointer :: ptx, pty complex(psb_spk_), allocatable :: xta(:), yta(:) integer(psb_ipk_) :: i, j, nr1, nc1,nr2, nc2 ,& - & map_kind, nr, ictxt + & map_kind, nr, ictxt, iam, np character(len=20), parameter :: name='psb_map_X2Yv' info = psb_success_ @@ -139,54 +141,66 @@ subroutine psb_c_map_X2Y_vect(alpha,x,beta,y,map,info,work) case(psb_map_aggr_) ictxt = map%p_desc_Y%get_context() + call psb_info(ictxt,iam,np) nr2 = map%p_desc_Y%get_global_rows() - nc2 = map%p_desc_Y%get_local_cols() - call yt%bld(nc2,mold=x%v) + nc2 = map%p_desc_Y%get_local_cols() + if (present(vty)) then + pty => vty + else + call yt%bld(nc2,mold=x%v) + pty => yt + end if if (info == psb_success_) call psb_halo(x,map%p_desc_X,info,work=work) - if (info == psb_success_) call psb_csmm(cone,map%map_X2Y,x,czero,yt,info) - if ((info == psb_success_) .and. map%p_desc_Y%is_repl()) then - yta = yt%get_vect() + if (info == psb_success_) call psb_csmm(cone,map%map_X2Y,x,czero,pty,info) + if ((info == psb_success_) .and. map%p_desc_Y%is_repl().and.(np>1)) then + yta = pty%get_vect() call psb_sum(ictxt,yta(1:nr2)) - call yt%set(yta) + call pty%set(yta) end if - if (info == psb_success_) call psb_geaxpby(alpha,yt,beta,y,map%p_desc_Y,info) + if (info == psb_success_) call psb_geaxpby(alpha,pty,beta,y,map%p_desc_Y,info) if (info /= psb_success_) then write(psb_err_unit,*) trim(name),' Error from inner routines',info info = -1 else - call yt%free(info) + if (.not.present(vty)) call yt%free(info) end if case(psb_map_gen_linear_) ictxt = map%desc_Y%get_context() + call psb_info(ictxt,iam,np) nr1 = map%desc_X%get_local_rows() nc1 = map%desc_X%get_local_cols() nr2 = map%desc_Y%get_global_rows() nc2 = map%desc_Y%get_local_cols() - call xt%bld(nc1,mold=x%v) - call yt%bld(nc2,mold=y%v) - - xta = x%get_vect() - call xt%set(xta(1:nr1)) - if (info == psb_success_) call psb_halo(xt,map%desc_X,info,work=work) - if (info == psb_success_) call psb_csmm(cone,map%map_X2Y,xt,czero,yt,info) - if ((info == psb_success_) .and. map%desc_Y%is_repl()) then - yta = yt%get_vect() + if (present(vtx).and.present(vty)) then + ptx => vtx + pty => vty + else + call xt%bld(nc1,mold=x%v) + call yt%bld(nc2,mold=y%v) + ptx => xt + pty => yt + end if + + call psb_geaxpby(cone,x,@XZERO,ptx,map%desc_X,info) + if (info == psb_success_) call psb_halo(ptx,map%desc_X,info,work=work) + if (info == psb_success_) call psb_csmm(cone,map%map_X2Y,ptx,czero,pty,info) + if ((info == psb_success_) .and. map%desc_Y%is_repl().and.(np>1)) then + yta = pty%get_vect() call psb_sum(ictxt,yta(1:nr2)) - call yt%set(yta) + call pty%set(yta) end if - if (info == psb_success_) call psb_geaxpby(alpha,yt,beta,y,map%desc_Y,info) + if (info == psb_success_) call psb_geaxpby(alpha,pty,beta,y,map%desc_Y,info) if (info /= psb_success_) then write(psb_err_unit,*) trim(name),' Error from inner routines',info info = -1 else - call xt%free(info) - call yt%free(info) + if (.not.present(vtx)) call xt%free(info) + if (.not.present(vty)) call yt%free(info) end if - case default write(psb_err_unit,*) trim(name),' Invalid descriptor input', & & map_kind, psb_map_aggr_, psb_map_gen_linear_ diff --git a/base/tools/psb_d_map.f90 b/base/tools/psb_d_map.f90 index e6165855..467d6b71 100644 --- a/base/tools/psb_d_map.f90 +++ b/base/tools/psb_d_map.f90 @@ -111,7 +111,7 @@ subroutine psb_d_map_X2Y(alpha,x,beta,y,map,info,work) end subroutine psb_d_map_X2Y -subroutine psb_d_map_X2Y_vect(alpha,x,beta,y,map,info,work) +subroutine psb_d_map_X2Y_vect(alpha,x,beta,y,map,info,work,vtx,vty) use psb_base_mod, psb_protect_name => psb_d_map_X2Y_vect implicit none type(psb_dlinmap_type), intent(in) :: map @@ -119,11 +119,13 @@ subroutine psb_d_map_X2Y_vect(alpha,x,beta,y,map,info,work) type(psb_d_vect_type), intent(inout) :: x,y integer(psb_ipk_), intent(out) :: info real(psb_dpk_), optional :: work(:) + type(psb_d_vect_type), optional, target, intent(inout) :: vtx,vty ! Local - type(psb_d_vect_type) :: xt, yt + type(psb_d_vect_type), target :: xt, yt + type(psb_d_vect_type),pointer :: ptx, pty real(psb_dpk_), allocatable :: xta(:), yta(:) integer(psb_ipk_) :: i, j, nr1, nc1,nr2, nc2 ,& - & map_kind, nr, ictxt + & map_kind, nr, ictxt, iam, np character(len=20), parameter :: name='psb_map_X2Yv' info = psb_success_ @@ -139,54 +141,66 @@ subroutine psb_d_map_X2Y_vect(alpha,x,beta,y,map,info,work) case(psb_map_aggr_) ictxt = map%p_desc_Y%get_context() + call psb_info(ictxt,iam,np) nr2 = map%p_desc_Y%get_global_rows() - nc2 = map%p_desc_Y%get_local_cols() - call yt%bld(nc2,mold=x%v) + nc2 = map%p_desc_Y%get_local_cols() + if (present(vty)) then + pty => vty + else + call yt%bld(nc2,mold=x%v) + pty => yt + end if if (info == psb_success_) call psb_halo(x,map%p_desc_X,info,work=work) - if (info == psb_success_) call psb_csmm(done,map%map_X2Y,x,dzero,yt,info) - if ((info == psb_success_) .and. map%p_desc_Y%is_repl()) then - yta = yt%get_vect() + if (info == psb_success_) call psb_csmm(done,map%map_X2Y,x,dzero,pty,info) + if ((info == psb_success_) .and. map%p_desc_Y%is_repl().and.(np>1)) then + yta = pty%get_vect() call psb_sum(ictxt,yta(1:nr2)) - call yt%set(yta) + call pty%set(yta) end if - if (info == psb_success_) call psb_geaxpby(alpha,yt,beta,y,map%p_desc_Y,info) + if (info == psb_success_) call psb_geaxpby(alpha,pty,beta,y,map%p_desc_Y,info) if (info /= psb_success_) then write(psb_err_unit,*) trim(name),' Error from inner routines',info info = -1 else - call yt%free(info) + if (.not.present(vty)) call yt%free(info) end if case(psb_map_gen_linear_) ictxt = map%desc_Y%get_context() + call psb_info(ictxt,iam,np) nr1 = map%desc_X%get_local_rows() nc1 = map%desc_X%get_local_cols() nr2 = map%desc_Y%get_global_rows() nc2 = map%desc_Y%get_local_cols() - call xt%bld(nc1,mold=x%v) - call yt%bld(nc2,mold=y%v) - - xta = x%get_vect() - call xt%set(xta(1:nr1)) - if (info == psb_success_) call psb_halo(xt,map%desc_X,info,work=work) - if (info == psb_success_) call psb_csmm(done,map%map_X2Y,xt,dzero,yt,info) - if ((info == psb_success_) .and. map%desc_Y%is_repl()) then - yta = yt%get_vect() + if (present(vtx).and.present(vty)) then + ptx => vtx + pty => vty + else + call xt%bld(nc1,mold=x%v) + call yt%bld(nc2,mold=y%v) + ptx => xt + pty => yt + end if + + call psb_geaxpby(done,x,@XZERO,ptx,map%desc_X,info) + if (info == psb_success_) call psb_halo(ptx,map%desc_X,info,work=work) + if (info == psb_success_) call psb_csmm(done,map%map_X2Y,ptx,dzero,pty,info) + if ((info == psb_success_) .and. map%desc_Y%is_repl().and.(np>1)) then + yta = pty%get_vect() call psb_sum(ictxt,yta(1:nr2)) - call yt%set(yta) + call pty%set(yta) end if - if (info == psb_success_) call psb_geaxpby(alpha,yt,beta,y,map%desc_Y,info) + if (info == psb_success_) call psb_geaxpby(alpha,pty,beta,y,map%desc_Y,info) if (info /= psb_success_) then write(psb_err_unit,*) trim(name),' Error from inner routines',info info = -1 else - call xt%free(info) - call yt%free(info) + if (.not.present(vtx)) call xt%free(info) + if (.not.present(vty)) call yt%free(info) end if - case default write(psb_err_unit,*) trim(name),' Invalid descriptor input', & & map_kind, psb_map_aggr_, psb_map_gen_linear_ diff --git a/base/tools/psb_s_map.f90 b/base/tools/psb_s_map.f90 index cb89af84..21343957 100644 --- a/base/tools/psb_s_map.f90 +++ b/base/tools/psb_s_map.f90 @@ -111,7 +111,7 @@ subroutine psb_s_map_X2Y(alpha,x,beta,y,map,info,work) end subroutine psb_s_map_X2Y -subroutine psb_s_map_X2Y_vect(alpha,x,beta,y,map,info,work) +subroutine psb_s_map_X2Y_vect(alpha,x,beta,y,map,info,work,vtx,vty) use psb_base_mod, psb_protect_name => psb_s_map_X2Y_vect implicit none type(psb_slinmap_type), intent(in) :: map @@ -119,11 +119,13 @@ subroutine psb_s_map_X2Y_vect(alpha,x,beta,y,map,info,work) type(psb_s_vect_type), intent(inout) :: x,y integer(psb_ipk_), intent(out) :: info real(psb_spk_), optional :: work(:) + type(psb_s_vect_type), optional, target, intent(inout) :: vtx,vty ! Local - type(psb_s_vect_type) :: xt, yt + type(psb_s_vect_type), target :: xt, yt + type(psb_s_vect_type),pointer :: ptx, pty real(psb_spk_), allocatable :: xta(:), yta(:) integer(psb_ipk_) :: i, j, nr1, nc1,nr2, nc2 ,& - & map_kind, nr, ictxt + & map_kind, nr, ictxt, iam, np character(len=20), parameter :: name='psb_map_X2Yv' info = psb_success_ @@ -139,54 +141,66 @@ subroutine psb_s_map_X2Y_vect(alpha,x,beta,y,map,info,work) case(psb_map_aggr_) ictxt = map%p_desc_Y%get_context() + call psb_info(ictxt,iam,np) nr2 = map%p_desc_Y%get_global_rows() - nc2 = map%p_desc_Y%get_local_cols() - call yt%bld(nc2,mold=x%v) + nc2 = map%p_desc_Y%get_local_cols() + if (present(vty)) then + pty => vty + else + call yt%bld(nc2,mold=x%v) + pty => yt + end if if (info == psb_success_) call psb_halo(x,map%p_desc_X,info,work=work) - if (info == psb_success_) call psb_csmm(sone,map%map_X2Y,x,szero,yt,info) - if ((info == psb_success_) .and. map%p_desc_Y%is_repl()) then - yta = yt%get_vect() + if (info == psb_success_) call psb_csmm(sone,map%map_X2Y,x,szero,pty,info) + if ((info == psb_success_) .and. map%p_desc_Y%is_repl().and.(np>1)) then + yta = pty%get_vect() call psb_sum(ictxt,yta(1:nr2)) - call yt%set(yta) + call pty%set(yta) end if - if (info == psb_success_) call psb_geaxpby(alpha,yt,beta,y,map%p_desc_Y,info) + if (info == psb_success_) call psb_geaxpby(alpha,pty,beta,y,map%p_desc_Y,info) if (info /= psb_success_) then write(psb_err_unit,*) trim(name),' Error from inner routines',info info = -1 else - call yt%free(info) + if (.not.present(vty)) call yt%free(info) end if case(psb_map_gen_linear_) ictxt = map%desc_Y%get_context() + call psb_info(ictxt,iam,np) nr1 = map%desc_X%get_local_rows() nc1 = map%desc_X%get_local_cols() nr2 = map%desc_Y%get_global_rows() nc2 = map%desc_Y%get_local_cols() - call xt%bld(nc1,mold=x%v) - call yt%bld(nc2,mold=y%v) - - xta = x%get_vect() - call xt%set(xta(1:nr1)) - if (info == psb_success_) call psb_halo(xt,map%desc_X,info,work=work) - if (info == psb_success_) call psb_csmm(sone,map%map_X2Y,xt,szero,yt,info) - if ((info == psb_success_) .and. map%desc_Y%is_repl()) then - yta = yt%get_vect() + if (present(vtx).and.present(vty)) then + ptx => vtx + pty => vty + else + call xt%bld(nc1,mold=x%v) + call yt%bld(nc2,mold=y%v) + ptx => xt + pty => yt + end if + + call psb_geaxpby(sone,x,@XZERO,ptx,map%desc_X,info) + if (info == psb_success_) call psb_halo(ptx,map%desc_X,info,work=work) + if (info == psb_success_) call psb_csmm(sone,map%map_X2Y,ptx,szero,pty,info) + if ((info == psb_success_) .and. map%desc_Y%is_repl().and.(np>1)) then + yta = pty%get_vect() call psb_sum(ictxt,yta(1:nr2)) - call yt%set(yta) + call pty%set(yta) end if - if (info == psb_success_) call psb_geaxpby(alpha,yt,beta,y,map%desc_Y,info) + if (info == psb_success_) call psb_geaxpby(alpha,pty,beta,y,map%desc_Y,info) if (info /= psb_success_) then write(psb_err_unit,*) trim(name),' Error from inner routines',info info = -1 else - call xt%free(info) - call yt%free(info) + if (.not.present(vtx)) call xt%free(info) + if (.not.present(vty)) call yt%free(info) end if - case default write(psb_err_unit,*) trim(name),' Invalid descriptor input', & & map_kind, psb_map_aggr_, psb_map_gen_linear_ diff --git a/base/tools/psb_z_map.f90 b/base/tools/psb_z_map.f90 index b7f11982..dc92ca73 100644 --- a/base/tools/psb_z_map.f90 +++ b/base/tools/psb_z_map.f90 @@ -111,7 +111,7 @@ subroutine psb_z_map_X2Y(alpha,x,beta,y,map,info,work) end subroutine psb_z_map_X2Y -subroutine psb_z_map_X2Y_vect(alpha,x,beta,y,map,info,work) +subroutine psb_z_map_X2Y_vect(alpha,x,beta,y,map,info,work,vtx,vty) use psb_base_mod, psb_protect_name => psb_z_map_X2Y_vect implicit none type(psb_zlinmap_type), intent(in) :: map @@ -119,11 +119,13 @@ subroutine psb_z_map_X2Y_vect(alpha,x,beta,y,map,info,work) type(psb_z_vect_type), intent(inout) :: x,y integer(psb_ipk_), intent(out) :: info complex(psb_dpk_), optional :: work(:) + type(psb_z_vect_type), optional, target, intent(inout) :: vtx,vty ! Local - type(psb_z_vect_type) :: xt, yt + type(psb_z_vect_type), target :: xt, yt + type(psb_z_vect_type),pointer :: ptx, pty complex(psb_dpk_), allocatable :: xta(:), yta(:) integer(psb_ipk_) :: i, j, nr1, nc1,nr2, nc2 ,& - & map_kind, nr, ictxt + & map_kind, nr, ictxt, iam, np character(len=20), parameter :: name='psb_map_X2Yv' info = psb_success_ @@ -139,54 +141,66 @@ subroutine psb_z_map_X2Y_vect(alpha,x,beta,y,map,info,work) case(psb_map_aggr_) ictxt = map%p_desc_Y%get_context() + call psb_info(ictxt,iam,np) nr2 = map%p_desc_Y%get_global_rows() - nc2 = map%p_desc_Y%get_local_cols() - call yt%bld(nc2,mold=x%v) + nc2 = map%p_desc_Y%get_local_cols() + if (present(vty)) then + pty => vty + else + call yt%bld(nc2,mold=x%v) + pty => yt + end if if (info == psb_success_) call psb_halo(x,map%p_desc_X,info,work=work) - if (info == psb_success_) call psb_csmm(zone,map%map_X2Y,x,zzero,yt,info) - if ((info == psb_success_) .and. map%p_desc_Y%is_repl()) then - yta = yt%get_vect() + if (info == psb_success_) call psb_csmm(zone,map%map_X2Y,x,zzero,pty,info) + if ((info == psb_success_) .and. map%p_desc_Y%is_repl().and.(np>1)) then + yta = pty%get_vect() call psb_sum(ictxt,yta(1:nr2)) - call yt%set(yta) + call pty%set(yta) end if - if (info == psb_success_) call psb_geaxpby(alpha,yt,beta,y,map%p_desc_Y,info) + if (info == psb_success_) call psb_geaxpby(alpha,pty,beta,y,map%p_desc_Y,info) if (info /= psb_success_) then write(psb_err_unit,*) trim(name),' Error from inner routines',info info = -1 else - call yt%free(info) + if (.not.present(vty)) call yt%free(info) end if case(psb_map_gen_linear_) ictxt = map%desc_Y%get_context() + call psb_info(ictxt,iam,np) nr1 = map%desc_X%get_local_rows() nc1 = map%desc_X%get_local_cols() nr2 = map%desc_Y%get_global_rows() nc2 = map%desc_Y%get_local_cols() - call xt%bld(nc1,mold=x%v) - call yt%bld(nc2,mold=y%v) - - xta = x%get_vect() - call xt%set(xta(1:nr1)) - if (info == psb_success_) call psb_halo(xt,map%desc_X,info,work=work) - if (info == psb_success_) call psb_csmm(zone,map%map_X2Y,xt,zzero,yt,info) - if ((info == psb_success_) .and. map%desc_Y%is_repl()) then - yta = yt%get_vect() + if (present(vtx).and.present(vty)) then + ptx => vtx + pty => vty + else + call xt%bld(nc1,mold=x%v) + call yt%bld(nc2,mold=y%v) + ptx => xt + pty => yt + end if + + call psb_geaxpby(zone,x,@XZERO,ptx,map%desc_X,info) + if (info == psb_success_) call psb_halo(ptx,map%desc_X,info,work=work) + if (info == psb_success_) call psb_csmm(zone,map%map_X2Y,ptx,zzero,pty,info) + if ((info == psb_success_) .and. map%desc_Y%is_repl().and.(np>1)) then + yta = pty%get_vect() call psb_sum(ictxt,yta(1:nr2)) - call yt%set(yta) + call pty%set(yta) end if - if (info == psb_success_) call psb_geaxpby(alpha,yt,beta,y,map%desc_Y,info) + if (info == psb_success_) call psb_geaxpby(alpha,pty,beta,y,map%desc_Y,info) if (info /= psb_success_) then write(psb_err_unit,*) trim(name),' Error from inner routines',info info = -1 else - call xt%free(info) - call yt%free(info) + if (.not.present(vtx)) call xt%free(info) + if (.not.present(vty)) call yt%free(info) end if - case default write(psb_err_unit,*) trim(name),' Invalid descriptor input', & & map_kind, psb_map_aggr_, psb_map_gen_linear_