From 0c3f67ae3581d0142b8157edaa27cad8315264b4 Mon Sep 17 00:00:00 2001 From: sfilippone Date: Mon, 8 Jun 2026 15:29:00 +0200 Subject: [PATCH] Refactor rgmres step 2 --- linsolve/impl/psb_crgmres.f90 | 21 +++++++++++++-------- linsolve/impl/psb_drgmres.f90 | 21 +++++++++++++-------- linsolve/impl/psb_srgmres.f90 | 21 +++++++++++++-------- linsolve/impl/psb_zrgmres.f90 | 21 +++++++++++++-------- 4 files changed, 52 insertions(+), 32 deletions(-) diff --git a/linsolve/impl/psb_crgmres.f90 b/linsolve/impl/psb_crgmres.f90 index cf6e8dc67..c4bb98605 100644 --- a/linsolve/impl/psb_crgmres.f90 +++ b/linsolve/impl/psb_crgmres.f90 @@ -285,6 +285,7 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& end if call psb_spmm(-cone,a,x,cone,v(1),desc_a,info,work=aux) + if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info) if (info /= psb_success_) then info=psb_err_from_subroutine_non_ call psb_errpush(info,name) @@ -401,10 +402,11 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& ! rst = rs call psb_geaxpby(cone,x,czero,xt,desc_a,info) - call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info) + call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info,s2=s2) call psb_geaxpby(cone,b,czero,w1,desc_a,info) call psb_spmm(-cone,a,xt,cone,w1,desc_a,info,work=aux) + if (present(s1)) call psb_gemlt(s1,w,desc_a,info) rni = psb_geamax(w1,desc_a,info) xni = psb_geamax(xt,desc_a,info) errnum = rni @@ -442,7 +444,8 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& call psb_geaxpby(cone,xt,czero,x,desc_a,info) ! = x = xt case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_) - call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info) ! + call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2) + ! end select @@ -462,7 +465,7 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,& call psb_geaxpby(cone,xt,czero,x,desc_a,info)! x = xt case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_) - call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info) ! + call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2) ! end select if (itx >= itmax_) then @@ -533,11 +536,12 @@ contains ! Rebuild solution X from the space V using the factor ! stored in R ! - subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info) + subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2) complex(psb_spk_) :: c(:), s(:), rs(:), h(:,:) type(psb_c_vect_type) :: v(:), w, w1, x type(psb_desc_type) :: desc_a class(psb_cprec_type) :: prec + type(psb_c_vect_type), intent(inout), optional :: s2 integer(psb_ipk_) :: info integer(psb_ipk_) :: k,n @@ -549,12 +553,13 @@ contains if (debug_level >= psb_debug_ext_) & & write(debug_unit,*) me,' ',trim(name),& & ' Rebuild x-> RS:',rs(1:n) - call w1%zero() + call w%zero() do k=1, n - call psb_geaxpby(rs(k),v(k),cone,w1,desc_a,info) + call psb_geaxpby(rs(k),v(k),cone,w,desc_a,info) end do - call prec%apply(w1,w,desc_a,info) - call psb_geaxpby(cone,w,cone,x,desc_a,info) + if (present(s2)) call psb_gediv(s2,w,desc_a,info) + call prec%apply(w,w1,desc_a,info) + call psb_geaxpby(cone,w1,cone,x,desc_a,info) end subroutine rebuildx end subroutine psb_crgmres_vect diff --git a/linsolve/impl/psb_drgmres.f90 b/linsolve/impl/psb_drgmres.f90 index 7cd4bb47b..a2306cede 100644 --- a/linsolve/impl/psb_drgmres.f90 +++ b/linsolve/impl/psb_drgmres.f90 @@ -285,6 +285,7 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& end if call psb_spmm(-done,a,x,done,v(1),desc_a,info,work=aux) + if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info) if (info /= psb_success_) then info=psb_err_from_subroutine_non_ call psb_errpush(info,name) @@ -401,10 +402,11 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& ! rst = rs call psb_geaxpby(done,x,dzero,xt,desc_a,info) - call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info) + call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info,s2=s2) call psb_geaxpby(done,b,dzero,w1,desc_a,info) call psb_spmm(-done,a,xt,done,w1,desc_a,info,work=aux) + if (present(s1)) call psb_gemlt(s1,w,desc_a,info) rni = psb_geamax(w1,desc_a,info) xni = psb_geamax(xt,desc_a,info) errnum = rni @@ -442,7 +444,8 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& call psb_geaxpby(done,xt,dzero,x,desc_a,info) ! = x = xt case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_) - call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info) ! + call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2) + ! end select @@ -462,7 +465,7 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,& call psb_geaxpby(done,xt,dzero,x,desc_a,info)! x = xt case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_) - call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info) ! + call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2) ! end select if (itx >= itmax_) then @@ -533,11 +536,12 @@ contains ! Rebuild solution X from the space V using the factor ! stored in R ! - subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info) + subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2) real(psb_dpk_) :: c(:), s(:), rs(:), h(:,:) type(psb_d_vect_type) :: v(:), w, w1, x type(psb_desc_type) :: desc_a class(psb_dprec_type) :: prec + type(psb_d_vect_type), intent(inout), optional :: s2 integer(psb_ipk_) :: info integer(psb_ipk_) :: k,n @@ -549,12 +553,13 @@ contains if (debug_level >= psb_debug_ext_) & & write(debug_unit,*) me,' ',trim(name),& & ' Rebuild x-> RS:',rs(1:n) - call w1%zero() + call w%zero() do k=1, n - call psb_geaxpby(rs(k),v(k),done,w1,desc_a,info) + call psb_geaxpby(rs(k),v(k),done,w,desc_a,info) end do - call prec%apply(w1,w,desc_a,info) - call psb_geaxpby(done,w,done,x,desc_a,info) + if (present(s2)) call psb_gediv(s2,w,desc_a,info) + call prec%apply(w,w1,desc_a,info) + call psb_geaxpby(done,w1,done,x,desc_a,info) end subroutine rebuildx end subroutine psb_drgmres_vect diff --git a/linsolve/impl/psb_srgmres.f90 b/linsolve/impl/psb_srgmres.f90 index fc2dfca11..18cd85402 100644 --- a/linsolve/impl/psb_srgmres.f90 +++ b/linsolve/impl/psb_srgmres.f90 @@ -285,6 +285,7 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& end if call psb_spmm(-sone,a,x,sone,v(1),desc_a,info,work=aux) + if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info) if (info /= psb_success_) then info=psb_err_from_subroutine_non_ call psb_errpush(info,name) @@ -401,10 +402,11 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& ! rst = rs call psb_geaxpby(sone,x,szero,xt,desc_a,info) - call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info) + call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info,s2=s2) call psb_geaxpby(sone,b,szero,w1,desc_a,info) call psb_spmm(-sone,a,xt,sone,w1,desc_a,info,work=aux) + if (present(s1)) call psb_gemlt(s1,w,desc_a,info) rni = psb_geamax(w1,desc_a,info) xni = psb_geamax(xt,desc_a,info) errnum = rni @@ -442,7 +444,8 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& call psb_geaxpby(sone,xt,szero,x,desc_a,info) ! = x = xt case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_) - call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info) ! + call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2) + ! end select @@ -462,7 +465,7 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,& call psb_geaxpby(sone,xt,szero,x,desc_a,info)! x = xt case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_) - call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info) ! + call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2) ! end select if (itx >= itmax_) then @@ -533,11 +536,12 @@ contains ! Rebuild solution X from the space V using the factor ! stored in R ! - subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info) + subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2) real(psb_spk_) :: c(:), s(:), rs(:), h(:,:) type(psb_s_vect_type) :: v(:), w, w1, x type(psb_desc_type) :: desc_a class(psb_sprec_type) :: prec + type(psb_s_vect_type), intent(inout), optional :: s2 integer(psb_ipk_) :: info integer(psb_ipk_) :: k,n @@ -549,12 +553,13 @@ contains if (debug_level >= psb_debug_ext_) & & write(debug_unit,*) me,' ',trim(name),& & ' Rebuild x-> RS:',rs(1:n) - call w1%zero() + call w%zero() do k=1, n - call psb_geaxpby(rs(k),v(k),sone,w1,desc_a,info) + call psb_geaxpby(rs(k),v(k),sone,w,desc_a,info) end do - call prec%apply(w1,w,desc_a,info) - call psb_geaxpby(sone,w,sone,x,desc_a,info) + if (present(s2)) call psb_gediv(s2,w,desc_a,info) + call prec%apply(w,w1,desc_a,info) + call psb_geaxpby(sone,w1,sone,x,desc_a,info) end subroutine rebuildx end subroutine psb_srgmres_vect diff --git a/linsolve/impl/psb_zrgmres.f90 b/linsolve/impl/psb_zrgmres.f90 index 63ca823d1..6e683dd4f 100644 --- a/linsolve/impl/psb_zrgmres.f90 +++ b/linsolve/impl/psb_zrgmres.f90 @@ -285,6 +285,7 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& end if call psb_spmm(-zone,a,x,zone,v(1),desc_a,info,work=aux) + if (present(s1)) call psb_gemlt(s1,v(1),desc_a,info) if (info /= psb_success_) then info=psb_err_from_subroutine_non_ call psb_errpush(info,name) @@ -401,10 +402,11 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& ! rst = rs call psb_geaxpby(zone,x,zzero,xt,desc_a,info) - call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info) + call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info,s2=s2) call psb_geaxpby(zone,b,zzero,w1,desc_a,info) call psb_spmm(-zone,a,xt,zone,w1,desc_a,info,work=aux) + if (present(s1)) call psb_gemlt(s1,w,desc_a,info) rni = psb_geamax(w1,desc_a,info) xni = psb_geamax(xt,desc_a,info) errnum = rni @@ -442,7 +444,8 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& call psb_geaxpby(zone,xt,zzero,x,desc_a,info) ! = x = xt case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_) - call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info) ! + call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2) + ! end select @@ -462,7 +465,7 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,& call psb_geaxpby(zone,xt,zzero,x,desc_a,info)! x = xt case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_) - call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info) ! + call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2=s2) ! end select if (itx >= itmax_) then @@ -533,11 +536,12 @@ contains ! Rebuild solution X from the space V using the factor ! stored in R ! - subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info) + subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info,s2) complex(psb_dpk_) :: c(:), s(:), rs(:), h(:,:) type(psb_z_vect_type) :: v(:), w, w1, x type(psb_desc_type) :: desc_a class(psb_zprec_type) :: prec + type(psb_z_vect_type), intent(inout), optional :: s2 integer(psb_ipk_) :: info integer(psb_ipk_) :: k,n @@ -549,12 +553,13 @@ contains if (debug_level >= psb_debug_ext_) & & write(debug_unit,*) me,' ',trim(name),& & ' Rebuild x-> RS:',rs(1:n) - call w1%zero() + call w%zero() do k=1, n - call psb_geaxpby(rs(k),v(k),zone,w1,desc_a,info) + call psb_geaxpby(rs(k),v(k),zone,w,desc_a,info) end do - call prec%apply(w1,w,desc_a,info) - call psb_geaxpby(zone,w,zone,x,desc_a,info) + if (present(s2)) call psb_gediv(s2,w,desc_a,info) + call prec%apply(w,w1,desc_a,info) + call psb_geaxpby(zone,w1,zone,x,desc_a,info) end subroutine rebuildx end subroutine psb_zrgmres_vect