From f0e69ab0dae40979a52bfe06ed0dc0559bef8f95 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Fri, 30 May 2008 13:08:11 +0000 Subject: [PATCH] psblas2-dev: base/modules/psb_serial_mod.f90 base/serial/psb_ctransp.f90 base/serial/psb_dtransp.f90 base/serial/psb_stransp.f90 base/serial/psb_ztransp.f90 Defined in-place transpose. krylov/psb_drgmres.f90 krylov/psb_srgmres.f90 Fixed for good alternate code path. --- base/modules/psb_serial_mod.f90 | 42 ++++++++++++++++++++---- base/serial/psb_ctransp.f90 | 41 +++++++++++++++++++++++ base/serial/psb_dtransp.f90 | 43 +++++++++++++++++++++++- base/serial/psb_stransp.f90 | 41 +++++++++++++++++++++++ base/serial/psb_ztransp.f90 | 41 +++++++++++++++++++++++ krylov/psb_drgmres.f90 | 55 +++++++------------------------ krylov/psb_srgmres.f90 | 58 ++++++++------------------------- 7 files changed, 225 insertions(+), 96 deletions(-) diff --git a/base/modules/psb_serial_mod.f90 b/base/modules/psb_serial_mod.f90 index 7ab2414c..e59ea422 100644 --- a/base/modules/psb_serial_mod.f90 +++ b/base/modules/psb_serial_mod.f90 @@ -526,40 +526,70 @@ module psb_serial_mod interface psb_transp subroutine psb_stransp(a,b,c,fmt) use psb_spmat_type - type(psb_sspmat_type) :: a,b + type(psb_sspmat_type), intent(in) :: a + type(psb_sspmat_type), intent(out) :: b integer, optional :: c character(len=*), optional :: fmt end subroutine psb_stransp + subroutine psb_stransp1(a,c,fmt) + use psb_spmat_type + type(psb_sspmat_type) :: a + integer, optional :: c + character(len=*), optional :: fmt + end subroutine psb_stransp1 subroutine psb_dtransp(a,b,c,fmt) use psb_spmat_type - type(psb_dspmat_type) :: a,b + type(psb_dspmat_type), intent(in) :: a + type(psb_dspmat_type), intent(out) :: b integer, optional :: c character(len=*), optional :: fmt end subroutine psb_dtransp + subroutine psb_dtransp1(a,c,fmt) + use psb_spmat_type + type(psb_dspmat_type) :: a + integer, optional :: c + character(len=*), optional :: fmt + end subroutine psb_dtransp1 subroutine psb_ctransp(a,b,c,fmt) use psb_spmat_type - type(psb_cspmat_type) :: a,b + type(psb_cspmat_type), intent(in) :: a + type(psb_cspmat_type), intent(out) :: b integer, optional :: c character(len=*), optional :: fmt end subroutine psb_ctransp + subroutine psb_ctransp1(a,c,fmt) + use psb_spmat_type + type(psb_cspmat_type) :: a + integer, optional :: c + character(len=*), optional :: fmt + end subroutine psb_ctransp1 subroutine psb_ztransp(a,b,c,fmt) use psb_spmat_type - type(psb_zspmat_type) :: a,b + type(psb_zspmat_type), intent(in) :: a + type(psb_zspmat_type), intent(out) :: b integer, optional :: c character(len=*), optional :: fmt end subroutine psb_ztransp + subroutine psb_ztransp1(a,c,fmt) + use psb_spmat_type + type(psb_zspmat_type) :: a + integer, optional :: c + character(len=*), optional :: fmt + end subroutine psb_ztransp1 end interface interface psb_transc subroutine psb_ctransc(a,b,c,fmt) use psb_spmat_type - type(psb_cspmat_type) :: a,b + type(psb_cspmat_type), intent(in) :: a + type(psb_cspmat_type), intent(out) :: b integer, optional :: c character(len=*), optional :: fmt end subroutine psb_ctransc subroutine psb_ztransc(a,b,c,fmt) use psb_spmat_type - type(psb_zspmat_type) :: a,b + type(psb_zspmat_type), intent(in) :: a + type(psb_zspmat_type), intent(out) :: b integer, optional :: c character(len=*), optional :: fmt end subroutine psb_ztransc diff --git a/base/serial/psb_ctransp.f90 b/base/serial/psb_ctransp.f90 index 91f1d869..d58f809c 100644 --- a/base/serial/psb_ctransp.f90 +++ b/base/serial/psb_ctransp.f90 @@ -78,3 +78,44 @@ subroutine psb_ctransp(a,b,c,fmt) return end subroutine psb_ctransp + +subroutine psb_ctransp1(a,c,fmt) + use psb_spmat_type + use psb_tools_mod + use psb_string_mod + use psb_serial_mod, psb_protect_name => psb_ctransp1 + implicit none + + type(psb_cspmat_type), intent(inout) :: a + integer, optional :: c + character(len=*), optional :: fmt + + character(len=5) :: fmt_ + integer ::c_, info + integer, allocatable :: itmp(:) + + if (present(c)) then + c_=c + else + c_=1 + endif + if (present(fmt)) then + fmt_ = psb_toupper(fmt) + else + fmt_='CSR' + endif + + call psb_spcnv(a,info,afmt='coo') + + if (info /= 0) then + write(0,*) 'transp: info from CSDP ',info + return + end if + call psb_transfer(a%ia1,itmp,info) + call psb_transfer(a%ia2,a%ia1,info) + call psb_transfer(itmp,a%ia2,info) + + call psb_spcnv(a,info,afmt=fmt_) + + return +end subroutine psb_ctransp1 diff --git a/base/serial/psb_dtransp.f90 b/base/serial/psb_dtransp.f90 index d1f4a3b8..b3cfa028 100644 --- a/base/serial/psb_dtransp.f90 +++ b/base/serial/psb_dtransp.f90 @@ -40,7 +40,7 @@ subroutine psb_dtransp(a,b,c,fmt) use psb_serial_mod, psb_protect_name => psb_dtransp implicit none - type(psb_dspmat_type), intent(inout) :: a + type(psb_dspmat_type), intent(in) :: a type(psb_dspmat_type), intent(out) :: b integer, optional :: c character(len=*), optional :: fmt @@ -78,3 +78,44 @@ subroutine psb_dtransp(a,b,c,fmt) return end subroutine psb_dtransp + +subroutine psb_dtransp1(a,c,fmt) + use psb_spmat_type + use psb_tools_mod + use psb_string_mod + use psb_serial_mod, psb_protect_name => psb_dtransp1 + implicit none + + type(psb_dspmat_type), intent(inout) :: a + integer, optional :: c + character(len=*), optional :: fmt + + character(len=5) :: fmt_ + integer ::c_, info + integer, allocatable :: itmp(:) + + if (present(c)) then + c_=c + else + c_=1 + endif + if (present(fmt)) then + fmt_ = psb_toupper(fmt) + else + fmt_='CSR' + endif + + call psb_spcnv(a,info,afmt='coo') + + if (info /= 0) then + write(0,*) 'transp: info from CSDP ',info + return + end if + call psb_transfer(a%ia1,itmp,info) + call psb_transfer(a%ia2,a%ia1,info) + call psb_transfer(itmp,a%ia2,info) + + call psb_spcnv(a,info,afmt=fmt_) + + return +end subroutine psb_dtransp1 diff --git a/base/serial/psb_stransp.f90 b/base/serial/psb_stransp.f90 index 213993a6..62e44a55 100644 --- a/base/serial/psb_stransp.f90 +++ b/base/serial/psb_stransp.f90 @@ -78,3 +78,44 @@ subroutine psb_stransp(a,b,c,fmt) return end subroutine psb_stransp + +subroutine psb_stransp1(a,c,fmt) + use psb_spmat_type + use psb_tools_mod + use psb_string_mod + use psb_serial_mod, psb_protect_name => psb_stransp1 + implicit none + + type(psb_sspmat_type), intent(inout) :: a + integer, optional :: c + character(len=*), optional :: fmt + + character(len=5) :: fmt_ + integer ::c_, info + integer, allocatable :: itmp(:) + + if (present(c)) then + c_=c + else + c_=1 + endif + if (present(fmt)) then + fmt_ = psb_toupper(fmt) + else + fmt_='CSR' + endif + + call psb_spcnv(a,info,afmt='coo') + + if (info /= 0) then + write(0,*) 'transp: info from CSDP ',info + return + end if + call psb_transfer(a%ia1,itmp,info) + call psb_transfer(a%ia2,a%ia1,info) + call psb_transfer(itmp,a%ia2,info) + + call psb_spcnv(a,info,afmt=fmt_) + + return +end subroutine psb_stransp1 diff --git a/base/serial/psb_ztransp.f90 b/base/serial/psb_ztransp.f90 index e085a604..abb6dd30 100644 --- a/base/serial/psb_ztransp.f90 +++ b/base/serial/psb_ztransp.f90 @@ -78,3 +78,44 @@ subroutine psb_ztransp(a,b,c,fmt) return end subroutine psb_ztransp + +subroutine psb_ztransp1(a,c,fmt) + use psb_spmat_type + use psb_tools_mod + use psb_string_mod + use psb_serial_mod, psb_protect_name => psb_ztransp1 + implicit none + + type(psb_zspmat_type), intent(inout) :: a + integer, optional :: c + character(len=*), optional :: fmt + + character(len=5) :: fmt_ + integer ::c_, info + integer, allocatable :: itmp(:) + + if (present(c)) then + c_=c + else + c_=1 + endif + if (present(fmt)) then + fmt_ = psb_toupper(fmt) + else + fmt_='CSR' + endif + + call psb_spcnv(a,info,afmt='coo') + + if (info /= 0) then + write(0,*) 'transp: info from CSDP ',info + return + end if + call psb_transfer(a%ia1,itmp,info) + call psb_transfer(a%ia2,a%ia1,info) + call psb_transfer(itmp,a%ia2,info) + + call psb_spcnv(a,info,afmt=fmt_) + + return +end subroutine psb_ztransp1 diff --git a/krylov/psb_drgmres.f90 b/krylov/psb_drgmres.f90 index 389f8d24..fc606371 100644 --- a/krylov/psb_drgmres.f90 +++ b/krylov/psb_drgmres.f90 @@ -341,35 +341,17 @@ subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist h(i+1,i) = psb_genrm2(w,desc_a,info) scal=done/h(i+1,i) call psb_geaxpby(scal,w,dzero,v(:,i+1),desc_a,info) - if (use_drot) then - do k=2,i - call drot(1,h(k-1,i),1,h(k,i),1,c(k-1),s(k-1)) - enddo - - rti = h(i,i) - rti1 = h(i+1,i) - call drotg(rti,rti1,c(i),s(i)) - call drot(1,h(i,i),1,h(i+1,i),1,c(i),s(i)) - h(i+1,i) = dzero - call drot(1,rs(i),1,rs(i+1),1,c(i),s(i)) - - else - do k=2,i - dt = h(k-1,i) - h(k-1,i) = c(k-1)*dt + s(k-1)*h(k,i) - h(k,i) = -s(k-1)*dt + c(k-1)*h(k,i) - enddo - gm = safe_dn2(h(i,i),h(i+1,i)) - if (debug_level >= psb_debug_ext_) & - & write(debug_unit,*) me,' ',trim(name),' GM : ',gm - gm = max(gm,epstol) - - c(i) = h(i,i)/gm - s(i) = h(i+1,i)/gm - rs(i+1) = -s(i)*rs(i) - rs(i) = c(i)*rs(i) - h(i,i) = c(i)*h(i,i)+s(i)*h(i+1,i) - endif + do k=2,i + call drot(1,h(k-1,i),1,h(k,i),1,c(k-1),s(k-1)) + enddo + + rti = h(i,i) + rti1 = h(i+1,i) + call drotg(rti,rti1,c(i),s(i)) + call drot(1,h(i,i),1,h(i+1,i),1,c(i),s(i)) + h(i+1,i) = dzero + call drot(1,rs(i),1,rs(i+1),1,c(i),s(i)) + if (istop_ == 1) then ! ! build x and then compute the residual and its infinity norm @@ -481,21 +463,6 @@ subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist end if return -contains - function safe_dn2(a,b) - real(psb_dpk_), intent(in) :: a, b - real(psb_dpk_) :: safe_dn2 - real(psb_dpk_) :: t - - t = max(abs(a),abs(b)) - if (t==0.d0) then - safe_dn2 = 0.d0 - else - safe_dn2 = t * sqrt(abs(a/t)**2 + abs(b/t)**2) - endif - return - end function safe_dn2 - end subroutine psb_drgmres diff --git a/krylov/psb_srgmres.f90 b/krylov/psb_srgmres.f90 index 2e84e7c1..978cc39c 100644 --- a/krylov/psb_srgmres.f90 +++ b/krylov/psb_srgmres.f90 @@ -129,7 +129,7 @@ subroutine psb_srgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist Real(psb_spk_) :: scal, gm, rti, rti1 Integer ::litmax, naux, mglob, it,k, itrace_,& & np,me, n_row, n_col, nl, int_err(5) - Logical, Parameter :: exchange=.True., noexchange=.False., use_srot=.true. + Logical, Parameter :: exchange=.True., noexchange=.False. Integer, Parameter :: irmax = 8 Integer :: itx, i, isvch, ictxt,istop_, err_act integer :: debug_level, debug_unit @@ -341,35 +341,18 @@ subroutine psb_srgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist h(i+1,i) = psb_genrm2(w,desc_a,info) scal=sone/h(i+1,i) call psb_geaxpby(scal,w,szero,v(:,i+1),desc_a,info) - if (use_srot) then - do k=2,i - call srot(1,h(k-1,i),1,h(k,i),1,c(k-1),s(k-1)) - enddo - - rti = h(i,i) - rti1 = h(i+1,i) - call srotg(rti,rti1,c(i),s(i)) - call srot(1,h(i,i),1,h(i+1,i),1,c(i),s(i)) - h(i+1,i) = szero - call srot(1,rs(i),1,rs(i+1),1,c(i),s(i)) - - else - do k=2,i - dt = h(k-1,i) - h(k-1,i) = c(k-1)*dt + s(k-1)*h(k,i) - h(k,i) = -s(k-1)*dt + c(k-1)*h(k,i) - enddo - gm = safe_dn2(h(i,i),h(i+1,i)) - if (debug_level >= psb_debug_ext_) & - & write(debug_unit,*) me,' ',trim(name),' GM : ',gm - gm = max(gm,s_epstol) - - c(i) = h(i,i)/gm - s(i) = h(i+1,i)/gm - rs(i+1) = -s(i)*rs(i) - rs(i) = c(i)*rs(i) - h(i,i) = c(i)*h(i,i)+s(i)*h(i+1,i) - endif + + do k=2,i + call srot(1,h(k-1,i),1,h(k,i),1,c(k-1),s(k-1)) + enddo + + rti = h(i,i) + rti1 = h(i+1,i) + call srotg(rti,rti1,c(i),s(i)) + call srot(1,h(i,i),1,h(i+1,i),1,c(i),s(i)) + h(i+1,i) = szero + call srot(1,rs(i),1,rs(i+1),1,c(i),s(i)) + if (istop_ == 1) then ! ! build x and then compute the residual and its infinity norm @@ -500,21 +483,6 @@ subroutine psb_srgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist end if return -contains - function safe_dn2(a,b) - real(psb_spk_), intent(in) :: a, b - real(psb_spk_) :: safe_dn2 - real(psb_spk_) :: t - - t = max(abs(a),abs(b)) - if (t==0.d0) then - safe_dn2 = 0.d0 - else - safe_dn2 = t * sqrt(abs(a/t)**2 + abs(b/t)**2) - endif - return - end function safe_dn2 - end subroutine psb_srgmres