From 85da2dbd27d145b94c126cf169937ed5fc69793d Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Wed, 15 Mar 2017 16:55:20 +0000 Subject: [PATCH] psblas3: base/tools/psb_c_map.f90 base/tools/psb_d_map.f90 base/tools/psb_s_map.f90 base/tools/psb_z_map.f90 Fixed X2Y and Y2X pointer magic. --- base/tools/psb_c_map.f90 | 74 +++++++++++++++++++++++++--------------- base/tools/psb_d_map.f90 | 74 +++++++++++++++++++++++++--------------- base/tools/psb_s_map.f90 | 74 +++++++++++++++++++++++++--------------- base/tools/psb_z_map.f90 | 74 +++++++++++++++++++++++++--------------- 4 files changed, 184 insertions(+), 112 deletions(-) diff --git a/base/tools/psb_c_map.f90 b/base/tools/psb_c_map.f90 index fbf7f057..18a6072d 100644 --- a/base/tools/psb_c_map.f90 +++ b/base/tools/psb_c_map.f90 @@ -147,7 +147,7 @@ subroutine psb_c_map_X2Y_vect(alpha,x,beta,y,map,info,work,vtx,vty) if (present(vty)) then pty => vty else - call yt%bld(nc2,mold=x%v) + call psb_geasb(yt,map%p_desc_Y,info,scratch=.true.,mold=x%v) pty => yt end if if (info == psb_success_) call psb_halo(x,map%p_desc_X,info,work=work) @@ -178,13 +178,13 @@ subroutine psb_c_map_X2Y_vect(alpha,x,beta,y,map,info,work,vtx,vty) ptx => vtx pty => vty else - call xt%bld(nc1,mold=x%v) - call yt%bld(nc2,mold=y%v) + call psb_geasb(xt,map%desc_X,info,scratch=.true.,mold=x%v) + call psb_geasb(yt,map%desc_Y,info,scratch=.true.,mold=x%v) ptx => xt pty => yt end if - call psb_geaxpby(cone,x,@XZERO,ptx,map%desc_X,info) + call psb_geaxpby(cone,x,czero,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 @@ -197,8 +197,10 @@ subroutine psb_c_map_X2Y_vect(alpha,x,beta,y,map,info,work,vtx,vty) write(psb_err_unit,*) trim(name),' Error from inner routines',info info = -1 else - if (.not.present(vtx)) call xt%free(info) - if (.not.present(vty)) call yt%free(info) + if (.not.(present(vtx).and.present(vty) )) then + call xt%free(info) + call yt%free(info) + end if end if case default @@ -290,7 +292,7 @@ subroutine psb_c_map_Y2X(alpha,x,beta,y,map,info,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_base_mod, psb_protect_name => psb_c_map_Y2X_vect implicit none type(psb_clinmap_type), intent(in) :: map @@ -298,8 +300,10 @@ subroutine psb_c_map_Y2X_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) :: xt, yt + type(psb_c_vect_type), optional, target, intent(inout) :: vtx,vty + ! Local + 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 @@ -320,20 +324,25 @@ subroutine psb_c_map_Y2X_vect(alpha,x,beta,y,map,info,work) ictxt = map%p_desc_X%get_context() nr2 = map%p_desc_X%get_global_rows() nc2 = map%p_desc_X%get_local_cols() - call yt%bld(nc2,mold=y%v) + if (present(vty)) then + pty => vty + else + call psb_geasb(yt,map%p_desc_X,info,scratch=.true.,mold=x%v) + pty => yt + end if if (info == psb_success_) call psb_halo(x,map%p_desc_Y,info,work=work) - if (info == psb_success_) call psb_csmm(cone,map%map_Y2X,x,czero,yt,info) - if ((info == psb_success_) .and. map%p_desc_X%is_repl()) then - yta = yt%get_vect() + if (info == psb_success_) call psb_csmm(cone,map%map_Y2X,x,czero,pty,info) + if ((info == psb_success_) .and. map%p_desc_X%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_X,info) + if (info == psb_success_) call psb_geaxpby(alpha,pty,beta,y,map%p_desc_X,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_) @@ -343,25 +352,34 @@ subroutine psb_c_map_Y2X_vect(alpha,x,beta,y,map,info,work) nc1 = map%desc_Y%get_local_cols() nr2 = map%desc_X%get_global_rows() nc2 = map%desc_X%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 (present(vtx).and.present(vty)) then + ptx => vtx + pty => vty + else + call psb_geasb(xt,map%desc_Y,info,scratch=.true.,mold=x%v) + call psb_geasb(yt,map%desc_X,info,scratch=.true.,mold=x%v) + ptx => xt + pty => yt + end if - if (info == psb_success_) call psb_halo(xt,map%desc_Y,info,work=work) - if (info == psb_success_) call psb_csmm(cone,map%map_Y2X,xt,czero,yt,info) - if ((info == psb_success_) .and. map%desc_X%is_repl()) then - yta = yt%get_vect() + call psb_geaxpby(cone,x,czero,ptx,map%desc_X,info) + + if (info == psb_success_) call psb_halo(ptx,map%desc_Y,info,work=work) + if (info == psb_success_) call psb_csmm(cone,map%map_Y2X,ptx,czero,pty,info) + if ((info == psb_success_) .and. map%desc_X%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_X,info) + if (info == psb_success_) call psb_geaxpby(alpha,pty,beta,y,map%desc_X,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).and.present(vty) )) then + call xt%free(info) + call yt%free(info) + end if end if case default diff --git a/base/tools/psb_d_map.f90 b/base/tools/psb_d_map.f90 index 467d6b71..0e827e91 100644 --- a/base/tools/psb_d_map.f90 +++ b/base/tools/psb_d_map.f90 @@ -147,7 +147,7 @@ subroutine psb_d_map_X2Y_vect(alpha,x,beta,y,map,info,work,vtx,vty) if (present(vty)) then pty => vty else - call yt%bld(nc2,mold=x%v) + call psb_geasb(yt,map%p_desc_Y,info,scratch=.true.,mold=x%v) pty => yt end if if (info == psb_success_) call psb_halo(x,map%p_desc_X,info,work=work) @@ -178,13 +178,13 @@ subroutine psb_d_map_X2Y_vect(alpha,x,beta,y,map,info,work,vtx,vty) ptx => vtx pty => vty else - call xt%bld(nc1,mold=x%v) - call yt%bld(nc2,mold=y%v) + call psb_geasb(xt,map%desc_X,info,scratch=.true.,mold=x%v) + call psb_geasb(yt,map%desc_Y,info,scratch=.true.,mold=x%v) ptx => xt pty => yt end if - call psb_geaxpby(done,x,@XZERO,ptx,map%desc_X,info) + call psb_geaxpby(done,x,dzero,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 @@ -197,8 +197,10 @@ subroutine psb_d_map_X2Y_vect(alpha,x,beta,y,map,info,work,vtx,vty) write(psb_err_unit,*) trim(name),' Error from inner routines',info info = -1 else - if (.not.present(vtx)) call xt%free(info) - if (.not.present(vty)) call yt%free(info) + if (.not.(present(vtx).and.present(vty) )) then + call xt%free(info) + call yt%free(info) + end if end if case default @@ -290,7 +292,7 @@ subroutine psb_d_map_Y2X(alpha,x,beta,y,map,info,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_base_mod, psb_protect_name => psb_d_map_Y2X_vect implicit none type(psb_dlinmap_type), intent(in) :: map @@ -298,8 +300,10 @@ subroutine psb_d_map_Y2X_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) :: xt, yt + type(psb_d_vect_type), optional, target, intent(inout) :: vtx,vty + ! Local + 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 @@ -320,20 +324,25 @@ subroutine psb_d_map_Y2X_vect(alpha,x,beta,y,map,info,work) ictxt = map%p_desc_X%get_context() nr2 = map%p_desc_X%get_global_rows() nc2 = map%p_desc_X%get_local_cols() - call yt%bld(nc2,mold=y%v) + if (present(vty)) then + pty => vty + else + call psb_geasb(yt,map%p_desc_X,info,scratch=.true.,mold=x%v) + pty => yt + end if if (info == psb_success_) call psb_halo(x,map%p_desc_Y,info,work=work) - if (info == psb_success_) call psb_csmm(done,map%map_Y2X,x,dzero,yt,info) - if ((info == psb_success_) .and. map%p_desc_X%is_repl()) then - yta = yt%get_vect() + if (info == psb_success_) call psb_csmm(done,map%map_Y2X,x,dzero,pty,info) + if ((info == psb_success_) .and. map%p_desc_X%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_X,info) + if (info == psb_success_) call psb_geaxpby(alpha,pty,beta,y,map%p_desc_X,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_) @@ -343,25 +352,34 @@ subroutine psb_d_map_Y2X_vect(alpha,x,beta,y,map,info,work) nc1 = map%desc_Y%get_local_cols() nr2 = map%desc_X%get_global_rows() nc2 = map%desc_X%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 (present(vtx).and.present(vty)) then + ptx => vtx + pty => vty + else + call psb_geasb(xt,map%desc_Y,info,scratch=.true.,mold=x%v) + call psb_geasb(yt,map%desc_X,info,scratch=.true.,mold=x%v) + ptx => xt + pty => yt + end if - if (info == psb_success_) call psb_halo(xt,map%desc_Y,info,work=work) - if (info == psb_success_) call psb_csmm(done,map%map_Y2X,xt,dzero,yt,info) - if ((info == psb_success_) .and. map%desc_X%is_repl()) then - yta = yt%get_vect() + call psb_geaxpby(done,x,dzero,ptx,map%desc_X,info) + + if (info == psb_success_) call psb_halo(ptx,map%desc_Y,info,work=work) + if (info == psb_success_) call psb_csmm(done,map%map_Y2X,ptx,dzero,pty,info) + if ((info == psb_success_) .and. map%desc_X%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_X,info) + if (info == psb_success_) call psb_geaxpby(alpha,pty,beta,y,map%desc_X,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).and.present(vty) )) then + call xt%free(info) + call yt%free(info) + end if end if case default diff --git a/base/tools/psb_s_map.f90 b/base/tools/psb_s_map.f90 index 21343957..9d0e8801 100644 --- a/base/tools/psb_s_map.f90 +++ b/base/tools/psb_s_map.f90 @@ -147,7 +147,7 @@ subroutine psb_s_map_X2Y_vect(alpha,x,beta,y,map,info,work,vtx,vty) if (present(vty)) then pty => vty else - call yt%bld(nc2,mold=x%v) + call psb_geasb(yt,map%p_desc_Y,info,scratch=.true.,mold=x%v) pty => yt end if if (info == psb_success_) call psb_halo(x,map%p_desc_X,info,work=work) @@ -178,13 +178,13 @@ subroutine psb_s_map_X2Y_vect(alpha,x,beta,y,map,info,work,vtx,vty) ptx => vtx pty => vty else - call xt%bld(nc1,mold=x%v) - call yt%bld(nc2,mold=y%v) + call psb_geasb(xt,map%desc_X,info,scratch=.true.,mold=x%v) + call psb_geasb(yt,map%desc_Y,info,scratch=.true.,mold=x%v) ptx => xt pty => yt end if - call psb_geaxpby(sone,x,@XZERO,ptx,map%desc_X,info) + call psb_geaxpby(sone,x,szero,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 @@ -197,8 +197,10 @@ subroutine psb_s_map_X2Y_vect(alpha,x,beta,y,map,info,work,vtx,vty) write(psb_err_unit,*) trim(name),' Error from inner routines',info info = -1 else - if (.not.present(vtx)) call xt%free(info) - if (.not.present(vty)) call yt%free(info) + if (.not.(present(vtx).and.present(vty) )) then + call xt%free(info) + call yt%free(info) + end if end if case default @@ -290,7 +292,7 @@ subroutine psb_s_map_Y2X(alpha,x,beta,y,map,info,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_base_mod, psb_protect_name => psb_s_map_Y2X_vect implicit none type(psb_slinmap_type), intent(in) :: map @@ -298,8 +300,10 @@ subroutine psb_s_map_Y2X_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) :: xt, yt + type(psb_s_vect_type), optional, target, intent(inout) :: vtx,vty + ! Local + 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 @@ -320,20 +324,25 @@ subroutine psb_s_map_Y2X_vect(alpha,x,beta,y,map,info,work) ictxt = map%p_desc_X%get_context() nr2 = map%p_desc_X%get_global_rows() nc2 = map%p_desc_X%get_local_cols() - call yt%bld(nc2,mold=y%v) + if (present(vty)) then + pty => vty + else + call psb_geasb(yt,map%p_desc_X,info,scratch=.true.,mold=x%v) + pty => yt + end if if (info == psb_success_) call psb_halo(x,map%p_desc_Y,info,work=work) - if (info == psb_success_) call psb_csmm(sone,map%map_Y2X,x,szero,yt,info) - if ((info == psb_success_) .and. map%p_desc_X%is_repl()) then - yta = yt%get_vect() + if (info == psb_success_) call psb_csmm(sone,map%map_Y2X,x,szero,pty,info) + if ((info == psb_success_) .and. map%p_desc_X%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_X,info) + if (info == psb_success_) call psb_geaxpby(alpha,pty,beta,y,map%p_desc_X,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_) @@ -343,25 +352,34 @@ subroutine psb_s_map_Y2X_vect(alpha,x,beta,y,map,info,work) nc1 = map%desc_Y%get_local_cols() nr2 = map%desc_X%get_global_rows() nc2 = map%desc_X%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 (present(vtx).and.present(vty)) then + ptx => vtx + pty => vty + else + call psb_geasb(xt,map%desc_Y,info,scratch=.true.,mold=x%v) + call psb_geasb(yt,map%desc_X,info,scratch=.true.,mold=x%v) + ptx => xt + pty => yt + end if - if (info == psb_success_) call psb_halo(xt,map%desc_Y,info,work=work) - if (info == psb_success_) call psb_csmm(sone,map%map_Y2X,xt,szero,yt,info) - if ((info == psb_success_) .and. map%desc_X%is_repl()) then - yta = yt%get_vect() + call psb_geaxpby(sone,x,szero,ptx,map%desc_X,info) + + if (info == psb_success_) call psb_halo(ptx,map%desc_Y,info,work=work) + if (info == psb_success_) call psb_csmm(sone,map%map_Y2X,ptx,szero,pty,info) + if ((info == psb_success_) .and. map%desc_X%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_X,info) + if (info == psb_success_) call psb_geaxpby(alpha,pty,beta,y,map%desc_X,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).and.present(vty) )) then + call xt%free(info) + call yt%free(info) + end if end if case default diff --git a/base/tools/psb_z_map.f90 b/base/tools/psb_z_map.f90 index dc92ca73..c65be176 100644 --- a/base/tools/psb_z_map.f90 +++ b/base/tools/psb_z_map.f90 @@ -147,7 +147,7 @@ subroutine psb_z_map_X2Y_vect(alpha,x,beta,y,map,info,work,vtx,vty) if (present(vty)) then pty => vty else - call yt%bld(nc2,mold=x%v) + call psb_geasb(yt,map%p_desc_Y,info,scratch=.true.,mold=x%v) pty => yt end if if (info == psb_success_) call psb_halo(x,map%p_desc_X,info,work=work) @@ -178,13 +178,13 @@ subroutine psb_z_map_X2Y_vect(alpha,x,beta,y,map,info,work,vtx,vty) ptx => vtx pty => vty else - call xt%bld(nc1,mold=x%v) - call yt%bld(nc2,mold=y%v) + call psb_geasb(xt,map%desc_X,info,scratch=.true.,mold=x%v) + call psb_geasb(yt,map%desc_Y,info,scratch=.true.,mold=x%v) ptx => xt pty => yt end if - call psb_geaxpby(zone,x,@XZERO,ptx,map%desc_X,info) + call psb_geaxpby(zone,x,zzero,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 @@ -197,8 +197,10 @@ subroutine psb_z_map_X2Y_vect(alpha,x,beta,y,map,info,work,vtx,vty) write(psb_err_unit,*) trim(name),' Error from inner routines',info info = -1 else - if (.not.present(vtx)) call xt%free(info) - if (.not.present(vty)) call yt%free(info) + if (.not.(present(vtx).and.present(vty) )) then + call xt%free(info) + call yt%free(info) + end if end if case default @@ -290,7 +292,7 @@ subroutine psb_z_map_Y2X(alpha,x,beta,y,map,info,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_base_mod, psb_protect_name => psb_z_map_Y2X_vect implicit none type(psb_zlinmap_type), intent(in) :: map @@ -298,8 +300,10 @@ subroutine psb_z_map_Y2X_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) :: xt, yt + type(psb_z_vect_type), optional, target, intent(inout) :: vtx,vty + ! Local + 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 @@ -320,20 +324,25 @@ subroutine psb_z_map_Y2X_vect(alpha,x,beta,y,map,info,work) ictxt = map%p_desc_X%get_context() nr2 = map%p_desc_X%get_global_rows() nc2 = map%p_desc_X%get_local_cols() - call yt%bld(nc2,mold=y%v) + if (present(vty)) then + pty => vty + else + call psb_geasb(yt,map%p_desc_X,info,scratch=.true.,mold=x%v) + pty => yt + end if if (info == psb_success_) call psb_halo(x,map%p_desc_Y,info,work=work) - if (info == psb_success_) call psb_csmm(zone,map%map_Y2X,x,zzero,yt,info) - if ((info == psb_success_) .and. map%p_desc_X%is_repl()) then - yta = yt%get_vect() + if (info == psb_success_) call psb_csmm(zone,map%map_Y2X,x,zzero,pty,info) + if ((info == psb_success_) .and. map%p_desc_X%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_X,info) + if (info == psb_success_) call psb_geaxpby(alpha,pty,beta,y,map%p_desc_X,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_) @@ -343,25 +352,34 @@ subroutine psb_z_map_Y2X_vect(alpha,x,beta,y,map,info,work) nc1 = map%desc_Y%get_local_cols() nr2 = map%desc_X%get_global_rows() nc2 = map%desc_X%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 (present(vtx).and.present(vty)) then + ptx => vtx + pty => vty + else + call psb_geasb(xt,map%desc_Y,info,scratch=.true.,mold=x%v) + call psb_geasb(yt,map%desc_X,info,scratch=.true.,mold=x%v) + ptx => xt + pty => yt + end if - if (info == psb_success_) call psb_halo(xt,map%desc_Y,info,work=work) - if (info == psb_success_) call psb_csmm(zone,map%map_Y2X,xt,zzero,yt,info) - if ((info == psb_success_) .and. map%desc_X%is_repl()) then - yta = yt%get_vect() + call psb_geaxpby(zone,x,zzero,ptx,map%desc_X,info) + + if (info == psb_success_) call psb_halo(ptx,map%desc_Y,info,work=work) + if (info == psb_success_) call psb_csmm(zone,map%map_Y2X,ptx,zzero,pty,info) + if ((info == psb_success_) .and. map%desc_X%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_X,info) + if (info == psb_success_) call psb_geaxpby(alpha,pty,beta,y,map%desc_X,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).and.present(vty) )) then + call xt%free(info) + call yt%free(info) + end if end if case default