fdurastante
/
psblas3
mirror of https://github.com/sfilippone/psblas3.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

psblas2-dev:

Browse Source 
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.
psblas3-type-indexed
Salvatore Filippone 17 years ago
parent 7a484fa3a0
commit f0e69ab0da
7 changed files with 225 additions and 96 deletions
Show Stats Download Patch File Download Diff File

42
base/modules/psb_serial_mod.f90
Unescape Escape View File

@ -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

41
base/serial/psb_ctransp.f90
Unescape Escape View File

@ -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

43
base/serial/psb_dtransp.f90
Unescape Escape View File

@ -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

41
base/serial/psb_stransp.f90
Unescape Escape View File

@ -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

41
base/serial/psb_ztransp.f90
Unescape Escape View File

@ -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

33
krylov/psb_drgmres.f90
Unescape Escape View File

@ -341,7 +341,6 @@ 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
@ -353,23 +352,6 @@ subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist
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
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

36
krylov/psb_srgmres.f90
Unescape Escape View File

@ -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,7 +341,7 @@ 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
@ -353,23 +353,6 @@ subroutine psb_srgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist
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
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

Write Preview
Loading…
Cancel
Save