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

psblas3-integer8:

Browse Source 
base/psblas/psb_caxpby.f90
 base/psblas/psb_cnrm2.f90
 base/psblas/psb_cnrmi.f90
 base/psblas/psb_cspmm.f90
 base/psblas/psb_cspsm.f90
 base/psblas/psb_daxpby.f90
 base/psblas/psb_dnrm2.f90
 base/psblas/psb_dnrmi.f90
 base/psblas/psb_dspmm.f90
 base/psblas/psb_dspnrm1.f90
 base/psblas/psb_dspsm.f90
 base/psblas/psb_saxpby.f90
 base/psblas/psb_snrm2.f90
 base/psblas/psb_snrmi.f90
 base/psblas/psb_sspmm.f90
 base/psblas/psb_sspsm.f90
 base/psblas/psb_zaxpby.f90
 base/psblas/psb_znrm2.f90
 base/psblas/psb_znrmi.f90
 base/psblas/psb_zspmm.f90
 base/psblas/psb_zspsm.f90


subdir psblas should work fine now.
psblas3-type-indexed
Salvatore Filippone 13 years ago
parent 934759b828
commit de5a399d55
21 changed files with 714 additions and 824 deletions
Show Stats Download Patch File Download Diff File

36
base/psblas/psb_caxpby.f90
Unescape Escape View File

@ -129,8 +129,8 @@ end subroutine psb_caxpby_vect
! Arguments: ! Arguments:
! alpha - complex,input The scalar used to multiply each component of X ! alpha - complex,input The scalar used to multiply each component of X
! x(:,:) - complex,input The input vector containing the entries of X ! x(:,:) - complex,input The input vector containing the entries of X
! beta - real,input The scalar used to multiply each component of Y ! beta - complex,input The scalar used to multiply each component of Y
! y(:,:) - real,inout The input vector Y ! y(:,:) - complex,inout The input vector Y
! desc_a - type(psb_desc_type) The communication descriptor. ! desc_a - type(psb_desc_type) The communication descriptor.
! info - integer Return code ! info - integer Return code
! jx - integer(optional) The column offset for X ! jx - integer(optional) The column offset for X
@ -150,7 +150,8 @@ subroutine psb_caxpby(alpha, x, beta,y,desc_a,info, n, jx, jy)
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ix, iy, ijx, ijy, m, iiy, in, jjy & err_act, iix, jjx, ix, iy, ijx, ijy, m, iiy, in, jjy, &
& lldx, lldy
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_geaxpby' name='psb_geaxpby'
@ -198,11 +199,12 @@ subroutine psb_caxpby(alpha, x, beta,y,desc_a,info, n, jx, jy)
end if end if
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
lldx = size(x,1)
lldy = size(y,1)
! check vector correctness ! check vector correctness
call psb_chkvect(m,ione,size(x,1),ix,ijx,desc_a,info,iix,jjx) call psb_chkvect(m,ione,lldx,ix,ijx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ione,size(y,1),iy,ijy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ione,lldy,iy,ijy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -219,8 +221,8 @@ subroutine psb_caxpby(alpha, x, beta,y,desc_a,info, n, jx, jy)
if ((in /= 0)) then if ((in /= 0)) then
if(desc_a%get_local_rows() > 0) then if(desc_a%get_local_rows() > 0) then
call caxpby(desc_a%get_local_cols(),in,& call caxpby(desc_a%get_local_cols(),in,&
& alpha,x(iix:,jjx),size(x,1),beta,& & alpha,x(iix:,jjx),lldx,beta,&
& y(iiy:,jjy),size(y,1),info) & y(iiy:,jjy),lldy,info)
end if end if
end if end if
@ -281,8 +283,8 @@ end subroutine psb_caxpby
! Arguments: ! Arguments:
! alpha - complex,input The scalar used to multiply each component of X ! alpha - complex,input The scalar used to multiply each component of X
! x(:) - complex,input The input vector containing the entries of X ! x(:) - complex,input The input vector containing the entries of X
! beta - real,input The scalar used to multiply each component of Y ! beta - complex,input The scalar used to multiply each component of Y
! y(:) - real,inout The input vector Y ! y(:) - complex,inout The input vector Y
! desc_a - type(psb_desc_type) The communication descriptor. ! desc_a - type(psb_desc_type) The communication descriptor.
! info - integer Return code ! info - integer Return code
! !
@ -299,7 +301,8 @@ subroutine psb_caxpbyv(alpha, x, beta,y,desc_a,info)
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ix, iy, m, iiy, jjy & err_act, iix, jjx, ix, iy, m, iiy, jjy, &
& lldx, lldy
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
@ -321,16 +324,17 @@ subroutine psb_caxpbyv(alpha, x, beta,y,desc_a,info)
iy = ione iy = ione
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
lldx = size(x,1)
lldy = size(y,1)
! check vector correctness ! check vector correctness
call psb_chkvect(m,ione,size(x),ix,ione,desc_a,info,iix,jjx) call psb_chkvect(m,ione,lldx,ix,ione,desc_a,info,iix,jjx)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect 1' ch_err='psb_chkvect 1'
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
call psb_chkvect(m,ione,size(y),iy,ione,desc_a,info,iiy,jjy) call psb_chkvect(m,ione,lldy,iy,ione,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect 2' ch_err='psb_chkvect 2'
@ -345,8 +349,8 @@ subroutine psb_caxpbyv(alpha, x, beta,y,desc_a,info)
if(desc_a%get_local_rows() > 0) then if(desc_a%get_local_rows() > 0) then
call caxpby(desc_a%get_local_cols(),ione,& call caxpby(desc_a%get_local_cols(),ione,&
& alpha,x,size(x),beta,& & alpha,x,lldx,beta,&
& y,size(y),info) & y,lldy,info)
end if end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)

169
base/psblas/psb_cnrm2.f90
Unescape Escape View File

@ -44,7 +44,7 @@
! info - integer. Return code ! info - integer. Return code
! jx - integer(optional). The column offset for sub( X ). ! jx - integer(optional). The column offset for sub( X ).
! !
function psb_cnrm2(x, desc_a, info, jx) function psb_cnrm2(x, desc_a, info, jx) result(res)
use psb_descriptor_type use psb_descriptor_type
use psb_check_mod use psb_check_mod
use psb_error_mod use psb_error_mod
@ -55,15 +55,13 @@ function psb_cnrm2(x, desc_a, info, jx)
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(in), optional :: jx integer(psb_ipk_), intent(in), optional :: jx
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
real(psb_spk_) :: psb_cnrm2 real(psb_spk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ndim, ix, ijx, i, m, id, idx, ndm & err_act, iix, jjx, ndim, ix, ijx, i, m, id, idx, ndm, ldx
real(psb_spk_) :: nrm2, scnrm2, dd real(psb_spk_) :: scnrm2, dd
character(len=20) :: name, ch_err
!!$ external scombnrm2
character(len=20) :: name, ch_err
name='psb_cnrm2' name='psb_cnrm2'
if(psb_get_errstatus() /= 0) return if(psb_get_errstatus() /= 0) return
@ -87,7 +85,8 @@ function psb_cnrm2(x, desc_a, info, jx)
endif endif
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
call psb_chkvect(m,1,size(x,1),ix,ijx,desc_a,info,iix,jjx) ldx = size(x,1)
call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -100,29 +99,22 @@ function psb_cnrm2(x, desc_a, info, jx)
goto 9999 goto 9999
end if end if
if(m /= 0) then if (desc_a%get_local_rows() > 0) then
if (desc_a%get_local_rows() > 0) then ndim = desc_a%get_local_rows()
ndim = desc_a%get_local_rows() res = scnrm2( int(ndim,kind=psb_mpik_), x(iix:,jjx), int(ione,kind=psb_mpik_) )
nrm2 = scnrm2( ndim, x(iix:,jjx), ione )
! adjust because overlapped elements are computed more than once
! adjust because overlapped elements are computed more than once do i=1,size(desc_a%ovrlap_elem,1)
do i=1,size(desc_a%ovrlap_elem,1) idx = desc_a%ovrlap_elem(i,1)
idx = desc_a%ovrlap_elem(i,1) ndm = desc_a%ovrlap_elem(i,2)
ndm = desc_a%ovrlap_elem(i,2) dd = real(ndm-1)/real(ndm)
dd = real(ndm-1)/real(ndm) res = res * sqrt(sone - dd*(abs(x(idx,jjx))/res)**2)
nrm2 = nrm2 * sqrt(sone - dd*(abs(x(idx,jjx))/nrm2)**2) end do
end do
else
nrm2 = dzero
end if
else else
nrm2 = dzero res = szero
end if end if
!!$ call pstreecomb(ictxt,'All',1,nrm2,-1,-1,scombnrm2) call psb_nrm2(ictxt,res)
call psb_nrm2(ictxt,nrm2)
psb_cnrm2 = nrm2
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -181,22 +173,22 @@ end function psb_cnrm2
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code ! info - integer. Return code
! !
function psb_cnrm2v(x, desc_a, info) function psb_cnrm2v(x, desc_a, info) result(res)
use psb_descriptor_type use psb_descriptor_type
use psb_check_mod use psb_check_mod
use psb_error_mod use psb_error_mod
use psb_penv_mod use psb_penv_mod
implicit none implicit none
complex(psb_spk_), intent(in) :: x(:) complex(psb_spk_), intent(in) :: x(:)
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
real(psb_spk_) :: psb_cnrm2v real(psb_spk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm & err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm, ldx
real(psb_spk_) :: nrm2, scnrm2, dd real(psb_spk_) :: scnrm2, dd
!!$ external scombnrm2 !!$ external scombnrm2
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
@ -218,8 +210,8 @@ function psb_cnrm2v(x, desc_a, info)
ix = 1 ix = 1
jx=1 jx=1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = size(x,1)
call psb_chkvect(m,1,size(x),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(m,ione,ldx,ix,jx,desc_a,info,iix,jjx)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -232,28 +224,22 @@ function psb_cnrm2v(x, desc_a, info)
goto 9999 goto 9999
end if end if
if(m /= 0) then if (desc_a%get_local_rows() > 0) then
if (desc_a%get_local_rows() > 0) then ndim = desc_a%get_local_rows()
ndim = desc_a%get_local_rows() res = scnrm2( int(ndim,kind=psb_mpik_), x, int(ione,kind=psb_mpik_) )
nrm2 = scnrm2( ndim, x, ione ) ! adjust because overlapped elements are computed more than once
! adjust because overlapped elements are computed more than once do i=1,size(desc_a%ovrlap_elem,1)
do i=1,size(desc_a%ovrlap_elem,1) idx = desc_a%ovrlap_elem(i,1)
idx = desc_a%ovrlap_elem(i,1) ndm = desc_a%ovrlap_elem(i,2)
ndm = desc_a%ovrlap_elem(i,2) dd = real(ndm-1)/real(ndm)
dd = real(ndm-1)/real(ndm) res = res * sqrt(sone - dd*(abs(x(idx))/res)**2)
nrm2 = nrm2 * sqrt(sone - dd*(abs(x(idx))/nrm2)**2) end do
end do
else
nrm2 = dzero
end if
else else
nrm2 = dzero res = szero
end if end if
!!$ call pstreecomb(ictxt,'All',1,nrm2,-1,-1,scombnrm2) call psb_nrm2(ictxt,res)
call psb_nrm2(ictxt,nrm2)
psb_cnrm2v = nrm2
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -281,12 +267,12 @@ function psb_cnrm2_vect(x, desc_a, info) result(res)
real(psb_spk_) :: res real(psb_spk_) :: res
type(psb_c_vect_type), intent (inout) :: x type(psb_c_vect_type), intent (inout) :: x
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm & err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm, ldx
real(psb_spk_) :: nrm2, snrm2, dd real(psb_spk_) :: snrm2, dd
!!$ external dcombnrm2 !!$ external dcombnrm2
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
@ -314,8 +300,8 @@ function psb_cnrm2_vect(x, desc_a, info) result(res)
ix = 1 ix = 1
jx=1 jx=1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = x%get_nrows()
call psb_chkvect(m,1,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(m,ione,ldx,ix,jx,desc_a,info,iix,jjx)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -328,10 +314,9 @@ function psb_cnrm2_vect(x, desc_a, info) result(res)
goto 9999 goto 9999
end if end if
if(m /= 0) then if (desc_a%get_local_rows() > 0) then
if (desc_a%get_local_rows() > 0) then ndim = desc_a%get_local_rows()
ndim = desc_a%get_local_rows() res = x%nrm2(ndim)
nrm2 = x%nrm2(ndim)
!!$ ! adjust because overlapped elements are computed more than once !!$ ! adjust because overlapped elements are computed more than once
!!$ do i=1,size(desc_a%ovrlap_elem,1) !!$ do i=1,size(desc_a%ovrlap_elem,1)
!!$ idx = desc_a%ovrlap_elem(i,1) !!$ idx = desc_a%ovrlap_elem(i,1)
@ -339,17 +324,11 @@ function psb_cnrm2_vect(x, desc_a, info) result(res)
!!$ dd = dble(ndm-1)/dble(ndm) !!$ dd = dble(ndm-1)/dble(ndm)
!!$ nrm2 = nrm2 * sqrt(done - dd*(abs(x(idx))/nrm2)**2) !!$ nrm2 = nrm2 * sqrt(done - dd*(abs(x(idx))/nrm2)**2)
!!$ end do !!$ end do
else
nrm2 = szero
end if
else else
nrm2 = szero res = szero
end if end if
!!$ call pdtreecomb(ictxt,'All',1,nrm2,-1,-1,dcombnrm2) call psb_nrm2(ictxt,res)
call psb_nrm2(ictxt,nrm2)
res = nrm2
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -416,14 +395,14 @@ subroutine psb_cnrm2vs(res, x, desc_a, info)
use psb_penv_mod use psb_penv_mod
implicit none implicit none
complex(psb_spk_), intent(in) :: x(:) complex(psb_spk_), intent(in) :: x(:)
real(psb_spk_), intent(out) :: res real(psb_spk_), intent(out) :: res
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm & err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm, ldx
real(psb_spk_) :: nrm2, scnrm2, dd real(psb_spk_) :: nrm2, scnrm2, dd
!!$ external scombnrm2 !!$ external scombnrm2
@ -446,8 +425,8 @@ subroutine psb_cnrm2vs(res, x, desc_a, info)
ix = 1 ix = 1
jx = 1 jx = 1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = size(x,1)
call psb_chkvect(m,1,size(x),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(m,ione,ldx,ix,jx,desc_a,info,iix,jjx)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -460,27 +439,23 @@ subroutine psb_cnrm2vs(res, x, desc_a, info)
goto 9999 goto 9999
end if end if
if(m /= 0) then if (desc_a%get_local_rows() > 0) then
if (desc_a%get_local_rows() > 0) then ndim = desc_a%get_local_rows()
ndim = desc_a%get_local_rows() res = scnrm2( int(ndim,kind=psb_mpik_), x, int(ione,kind=psb_mpik_) )
nrm2 = scnrm2( ndim, x, ione )
! adjust because overlapped elements are computed more than once ! adjust because overlapped elements are computed more than once
do i=1,size(desc_a%ovrlap_elem,1) do i=1,size(desc_a%ovrlap_elem,1)
idx = desc_a%ovrlap_elem(i,1) idx = desc_a%ovrlap_elem(i,1)
ndm = desc_a%ovrlap_elem(i,2) ndm = desc_a%ovrlap_elem(i,2)
dd = real(ndm-1)/real(ndm) dd = real(ndm-1)/real(ndm)
nrm2 = nrm2 * sqrt(sone - dd*(abs(x(idx))/nrm2)**2) res = res * sqrt(sone - dd*(abs(x(idx))/res)**2)
end do end do
else
nrm2 = dzero
end if
else else
nrm2 = dzero res = szero
end if end if
!!$ call pstreecomb(ictxt,'All',1,nrm2,-1,-1,scombnrm2) call psb_nrm2(ictxt,res)
call psb_nrm2(ictxt,nrm2)
res = nrm2
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

27
base/psblas/psb_cnrmi.f90
Unescape Escape View File

@ -32,28 +32,27 @@
! File: psb_cnrmi.f90 ! File: psb_cnrmi.f90
! !
! Function: psb_cnrmi ! Function: psb_cnrmi
! Forms the approximated norm of a sparse matrix, ! Forms the infinity norm of a sparse matrix,
! !
! normi := max(abs(sum(A(i,j)))) ! nrmi := max_i(abs(sum(A(i,:))))
! !
! Arguments: ! Arguments:
! a - type(psb_dspmat_type). The sparse matrix containing A. ! a - type(psb_cspmat_type). The sparse matrix containing A.
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code ! info - integer. Return code
! !
function psb_cnrmi(a,desc_a,info) function psb_cnrmi(a,desc_a,info) result(res)
use psb_base_mod, psb_protect_name => psb_cnrmi use psb_base_mod, psb_protect_name => psb_cnrmi
implicit none implicit none
type(psb_cspmat_type), intent(in) :: a type(psb_cspmat_type), intent(in) :: a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
real(psb_spk_) :: psb_cnrmi real(psb_spk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, n, iia, jja, ia, ja, mdim, ndim, m & err_act, n, iia, jja, ia, ja, mdim, ndim, m
real(psb_spk_) :: nrmi
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_cnrmi' name='psb_cnrmi'
@ -90,21 +89,19 @@ function psb_cnrmi(a,desc_a,info)
end if end if
if ((m /= 0).and.(n /= 0)) then if ((m /= 0).and.(n /= 0)) then
nrmi = a%csnmi() res = a%csnmi()
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_csnmi' ch_err='psb_cnrmi'
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
else else
nrmi = 0.0 res = szero
end if end if
! compute global max ! compute global max
call psb_amx(ictxt, nrmi) call psb_amx(ictxt, res)
psb_cnrmi = nrmi
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

30
base/psblas/psb_cspmm.f90
Unescape Escape View File

@ -205,9 +205,9 @@ subroutine psb_cspmm(alpha,a,x,beta,y,desc_a,info,&
end if end if
! checking for vectors correctness ! checking for vectors correctness
call psb_chkvect(n,ik,size(x,1),ix,ijx,desc_a,info,iix,jjx) call psb_chkvect(n,ik,lldx,ix,ijx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,size(y,1),iy,ijy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ik,lldy,iy,ijy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -277,9 +277,9 @@ subroutine psb_cspmm(alpha,a,x,beta,y,desc_a,info,&
! checking for vectors correctness ! checking for vectors correctness
call psb_chkvect(m,ik,size(x,1),ix,ijx,desc_a,info,iix,jjx) call psb_chkvect(m,ik,lldx,ix,ijx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(n,ik,size(y,1),iy,ijy,desc_a,info,iiy,jjy) & call psb_chkvect(n,ik,lldy,iy,ijy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -304,7 +304,8 @@ subroutine psb_cspmm(alpha,a,x,beta,y,desc_a,info,&
if (info == psb_success_) call psi_ovrl_upd(x,desc_a,psb_avg_,info) if (info == psb_success_) call psi_ovrl_upd(x,desc_a,psb_avg_,info)
y(nrow+1:ncol,1:ik) = czero y(nrow+1:ncol,1:ik) = czero
if (info == psb_success_) call psb_csmm(alpha,a,x(:,1:ik),beta,y(:,1:ik),info,trans=trans_) if (info == psb_success_) &
& call psb_csmm(alpha,a,x(:,1:ik),beta,y(:,1:ik),info,trans=trans_)
if (debug_level >= psb_debug_comp_) & if (debug_level >= psb_debug_comp_) &
& write(debug_unit,*) me,' ',trim(name),' csmm ', info & write(debug_unit,*) me,' ',trim(name),' csmm ', info
if (info /= psb_success_) then if (info /= psb_success_) then
@ -542,9 +543,9 @@ subroutine psb_cspmv(alpha,a,x,beta,y,desc_a,info,&
end if end if
! checking for vectors correctness ! checking for vectors correctness
call psb_chkvect(n,ik,size(x),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(n,ik,lldx,ix,jx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,size(y),iy,jy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ik,lldy,iy,jy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -582,9 +583,9 @@ subroutine psb_cspmv(alpha,a,x,beta,y,desc_a,info,&
end if end if
! checking for vectors correctness ! checking for vectors correctness
call psb_chkvect(m,ik,size(x),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(m,ik,lldx,ix,jx,desc_a,info,iix,jjx)
if (info == psb_success_)& if (info == psb_success_)&
& call psb_chkvect(n,ik,size(y),iy,jy,desc_a,info,iiy,jjy) & call psb_chkvect(n,ik,lldy,iy,jy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -815,9 +816,9 @@ subroutine psb_cspmv_vect(alpha,a,x,beta,y,desc_a,info,&
end if end if
! checking for vectors correctness ! checking for vectors correctness
call psb_chkvect(n,ik,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(n,ik,lldx,ix,jx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,y%get_nrows(),iy,jy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ik,lldy,iy,jy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -855,9 +856,9 @@ subroutine psb_cspmv_vect(alpha,a,x,beta,y,desc_a,info,&
end if end if
! checking for vectors correctness ! checking for vectors correctness
call psb_chkvect(m,ik,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(m,ik,lldx,ix,jx,desc_a,info,iix,jjx)
if (info == psb_success_)& if (info == psb_success_)&
& call psb_chkvect(n,ik,y%get_nrows(),iy,jy,desc_a,info,iiy,jjy) & call psb_chkvect(n,ik,lldy,iy,jy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -881,8 +882,9 @@ subroutine psb_cspmv_vect(alpha,a,x,beta,y,desc_a,info,&
! !
call psi_ovrl_save(x%v,xvsave,desc_a,info) call psi_ovrl_save(x%v,xvsave,desc_a,info)
if (info == psb_success_) call psi_ovrl_upd(x%v,desc_a,psb_avg_,info) if (info == psb_success_) call psi_ovrl_upd(x%v,desc_a,psb_avg_,info)
!!! THIS SHOULD BE FIXED !!! But beta is almost never /= 0 !!! THIS SHOULD BE FIXED !!! But beta is almost never /= 0
!!$ yp(nrow+1:ncol) = szero !!$ yp(nrow+1:ncol) = czero
! local Matrix-vector product ! local Matrix-vector product
if (info == psb_success_) call psb_csmm(alpha,a,x,beta,y,info,trans=trans_) if (info == psb_success_) call psb_csmm(alpha,a,x,beta,y,info,trans=trans_)

14
base/psblas/psb_cspsm.f90
Unescape Escape View File

@ -98,7 +98,7 @@ subroutine psb_cspsm(alpha,a,x,beta,y,desc_a,info,&
& m, nrow, ncol, liwork, llwork, iiy, jjy, idx, ndm & m, nrow, ncol, liwork, llwork, iiy, jjy, idx, ndm
character :: lscale character :: lscale
integer(psb_ipk_), parameter :: nb=4 integer(psb_ipk_), parameter :: nb=4
complex(psb_spk_),pointer :: iwork(:), xp(:,:), yp(:,:), id(:) complex(psb_spk_),pointer :: iwork(:), xp(:,:), yp(:,:), id(:)
character :: itrans character :: itrans
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
@ -220,9 +220,9 @@ subroutine psb_cspsm(alpha,a,x,beta,y,desc_a,info,&
call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja) call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja)
! checking for vectors correctness ! checking for vectors correctness
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,size(x,1),ix,ijx,desc_a,info,iix,jjx) & call psb_chkvect(m,ik,lldx,ix,ijx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,size(y,1),iy,ijy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ik,lldy,iy,ijy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect/mat' ch_err='psb_chkvect/mat'
@ -483,9 +483,9 @@ subroutine psb_cspsv(alpha,a,x,beta,y,desc_a,info,&
call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja) call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja)
! checking for vectors correctness ! checking for vectors correctness
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,size(x),ix,jx,desc_a,info,iix,jjx) & call psb_chkvect(m,ik,lldx,ix,jx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,size(y),iy,jy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ik,lldy,iy,jy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect/mat' ch_err='psb_chkvect/mat'
@ -684,9 +684,9 @@ subroutine psb_cspsv_vect(alpha,a,x,beta,y,desc_a,info,&
call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja) call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja)
! checking for vectors correctness ! checking for vectors correctness
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) & call psb_chkvect(m,ik,lldx,ix,jx,desc_a,info,iix,jjx)
if (info == psb_success_)& if (info == psb_success_)&
& call psb_chkvect(m,ik,y%get_nrows(),iy,jy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ik,lldy,iy,jy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect/mat' ch_err='psb_chkvect/mat'

74
base/psblas/psb_daxpby.f90
Unescape Escape View File

@ -31,14 +31,13 @@
!!$ !!$
! File: psb_daxpby.f90 ! File: psb_daxpby.f90
subroutine psb_daxpby_vect(alpha, x, beta, y,& subroutine psb_daxpby_vect(alpha, x, beta, y,&
& desc_a, info) & desc_a, info)
use psb_base_mod, psb_protect_name => psb_daxpby_vect use psb_base_mod, psb_protect_name => psb_daxpby_vect
implicit none implicit none
type(psb_d_vect_type), intent (inout) :: x type(psb_d_vect_type), intent (inout) :: x
type(psb_d_vect_type), intent (inout) :: y type(psb_d_vect_type), intent (inout) :: y
real(psb_dpk_), intent (in) :: alpha, beta real(psb_dpk_), intent (in) :: alpha, beta
type(psb_desc_type), intent (in) :: desc_a type(psb_desc_type), intent (in) :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -128,17 +127,18 @@ end subroutine psb_daxpby_vect
! sub( Y ) denotes Y(:,JY). ! sub( Y ) denotes Y(:,JY).
! !
! Arguments: ! Arguments:
! alpha - real The scalar alpha ! alpha - real,input The scalar used to multiply each component of X
! x(:,:) - real The input vector containing the entries of ( X ). ! x(:,:) - real,input The input vector containing the entries of X
! beta - real The scalar beta ! beta - real,input The scalar used to multiply each component of Y
! y(:,:) - real The input vector containing the entries of ( Y ). ! y(:,:) - real,inout The input vector Y
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type) The communication descriptor.
! info - integer. Return code ! info - integer Return code
! jx - integer(optional). The column offset for sub( X ). ! jx - integer(optional) The column offset for X
! jy - integer(optional). The column offset for sub( Y ). ! jy - integer(optional) The column offset for Y
! !
subroutine psb_daxpby(alpha, x, beta,y,desc_a,info, n, jx, jy) subroutine psb_daxpby(alpha, x, beta,y,desc_a,info, n, jx, jy)
use psb_base_mod, psb_protect_name => psb_daxpby use psb_base_mod, psb_protect_name => psb_daxpby
implicit none implicit none
integer(psb_ipk_), intent(in), optional :: n, jx, jy integer(psb_ipk_), intent(in), optional :: n, jx, jy
@ -150,16 +150,16 @@ subroutine psb_daxpby(alpha, x, beta,y,desc_a,info, n, jx, jy)
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ix, iy, ijx, ijy, m, iiy, in, jjy & err_act, iix, jjx, ix, iy, ijx, ijy, m, iiy, in, jjy, &
& lldx, lldy
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_dgeaxpby' name='psb_geaxpby'
if (psb_errstatus_fatal()) return if(psb_get_errstatus() /= 0) return
info=psb_success_ info=psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if (np == -ione) then if (np == -ione) then
info = psb_err_context_error_ info = psb_err_context_error_
@ -199,11 +199,12 @@ subroutine psb_daxpby(alpha, x, beta,y,desc_a,info, n, jx, jy)
end if end if
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
lldx = size(x,1)
lldy = size(y,1)
! check vector correctness ! check vector correctness
call psb_chkvect(m,ione,size(x,1),ix,ijx,desc_a,info,iix,jjx) call psb_chkvect(m,ione,lldx,ix,ijx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ione,size(y,1),iy,ijy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ione,lldy,iy,ijy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -219,9 +220,9 @@ subroutine psb_daxpby(alpha, x, beta,y,desc_a,info, n, jx, jy)
if ((in /= 0)) then if ((in /= 0)) then
if(desc_a%get_local_rows() > 0) then if(desc_a%get_local_rows() > 0) then
call daxpby(desc_a%get_local_rows(),in,& call caxpby(desc_a%get_local_cols(),in,&
& alpha,x(iix:,jjx),size(x,1),beta,& & alpha,x(iix:,jjx),lldx,beta,&
& y(iiy:,jjy),size(y,1),info) & y(iiy:,jjy),lldy,info)
end if end if
end if end if
@ -280,12 +281,12 @@ end subroutine psb_daxpby
! Y := beta * Y + alpha * X ! Y := beta * Y + alpha * X
! !
! Arguments: ! Arguments:
! alpha - real The scalar alpha ! alpha - real,input The scalar used to multiply each component of X
! x(:) - real The input vector containing the entries of ( X ). ! x(:) - real,input The input vector containing the entries of X
! beta - real The scalar beta ! beta - real,input The scalar used to multiply each component of Y
! y(:) - real The input vector containing the entries of ( Y ). ! y(:) - real,inout The input vector Y
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type) The communication descriptor.
! info - integer. Return code ! info - integer Return code
! !
! !
subroutine psb_daxpbyv(alpha, x, beta,y,desc_a,info) subroutine psb_daxpbyv(alpha, x, beta,y,desc_a,info)
@ -300,11 +301,13 @@ subroutine psb_daxpbyv(alpha, x, beta,y,desc_a,info)
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ix, iy, m, iiy, jjy & err_act, iix, jjx, ix, iy, m, iiy, jjy, &
& lldx, lldy
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
logical, parameter :: debug=.false.
name='psb_dgeaxpby' name='psb_geaxpby'
if (psb_errstatus_fatal()) return if(psb_get_errstatus() /= 0) return
info=psb_success_ info=psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -321,16 +324,17 @@ subroutine psb_daxpbyv(alpha, x, beta,y,desc_a,info)
iy = ione iy = ione
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
lldx = size(x,1)
lldy = size(y,1)
! check vector correctness ! check vector correctness
call psb_chkvect(m,ione,size(x),ix,ione,desc_a,info,iix,jjx) call psb_chkvect(m,ione,lldx,ix,ione,desc_a,info,iix,jjx)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect 1' ch_err='psb_chkvect 1'
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
call psb_chkvect(m,ione,size(y),iy,ione,desc_a,info,iiy,jjy) call psb_chkvect(m,ione,lldy,iy,ione,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect 2' ch_err='psb_chkvect 2'
@ -344,9 +348,9 @@ subroutine psb_daxpbyv(alpha, x, beta,y,desc_a,info)
end if end if
if(desc_a%get_local_rows() > 0) then if(desc_a%get_local_rows() > 0) then
call daxpby(desc_a%get_local_rows(),ione,& call caxpby(desc_a%get_local_cols(),ione,&
& alpha,x,size(x),beta,& & alpha,x,lldx,beta,&
& y,size(y),info) & y,lldy,info)
end if end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)

208
base/psblas/psb_dnrm2.f90
Unescape Escape View File

@ -39,12 +39,12 @@
! where sub( X ) denotes X(:,JX). ! where sub( X ) denotes X(:,JX).
! !
! Arguments: ! Arguments:
! x - real,dimension(:,:). The input vector containing the entries of X. ! x(:,:) - real The input vector containing the entries of sub( X ).
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code ! info - integer. Return code
! jx - integer(optional). The column offset for X . ! jx - integer(optional). The column offset for sub( X ).
! !
function psb_dnrm2(x, desc_a, info, jx) function psb_dnrm2(x, desc_a, info, jx) result(res)
use psb_descriptor_type use psb_descriptor_type
use psb_check_mod use psb_check_mod
use psb_error_mod use psb_error_mod
@ -55,17 +55,16 @@ function psb_dnrm2(x, desc_a, info, jx)
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(in), optional :: jx integer(psb_ipk_), intent(in), optional :: jx
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
real(psb_dpk_) :: psb_dnrm2 real(psb_dpk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ndim, ix, ijx, i, m, id, idx, ndm & err_act, iix, jjx, ndim, ix, ijx, i, m, id, idx, ndm, ldx
real(psb_dpk_) :: nrm2, dnrm2, dd real(psb_dpk_) :: dnrm2, dd
!!$ external dcombnrm2 character(len=20) :: name, ch_err
character(len=20) :: name, ch_err
name='psb_dnrm2' name='psb_dnrm2'
if (psb_errstatus_fatal()) return if(psb_get_errstatus() /= 0) return
info=psb_success_ info=psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -86,7 +85,8 @@ function psb_dnrm2(x, desc_a, info, jx)
endif endif
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
call psb_chkvect(m,1,size(x,1),ix,ijx,desc_a,info,iix,jjx) ldx = size(x,1)
call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -99,29 +99,22 @@ function psb_dnrm2(x, desc_a, info, jx)
goto 9999 goto 9999
end if end if
if(m /= 0) then if (desc_a%get_local_rows() > 0) then
if (desc_a%get_local_rows() > 0) then ndim = desc_a%get_local_rows()
ndim = desc_a%get_local_rows() res = dnrm2( int(ndim,kind=psb_mpik_), x(iix:,jjx), int(ione,kind=psb_mpik_) )
nrm2 = dnrm2( ndim, x(iix:,jjx), ione )
! adjust because overlapped elements are computed more than once
! adjust because overlapped elements are computed more than once do i=1,size(desc_a%ovrlap_elem,1)
do i=1,size(desc_a%ovrlap_elem,1) idx = desc_a%ovrlap_elem(i,1)
idx = desc_a%ovrlap_elem(i,1) ndm = desc_a%ovrlap_elem(i,2)
ndm = desc_a%ovrlap_elem(i,2) dd = real(ndm-1)/real(ndm)
dd = dble(ndm-1)/dble(ndm) res = res * sqrt(sone - dd*(abs(x(idx,jjx))/res)**2)
nrm2 = nrm2 * sqrt(done - dd*(abs(x(idx,jjx))/nrm2)**2) end do
end do
else
nrm2 = dzero
end if
else else
nrm2 = dzero res = dzero
end if end if
!!$ call pdtreecomb(ictxt,'All',1,nrm2,-1,-1,dcombnrm2) call psb_nrm2(ictxt,res)
call psb_nrm2(ictxt,nrm2)
psb_dnrm2 = nrm2
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -176,31 +169,32 @@ end function psb_dnrm2
! norm2 := sqrt ( X**T * X) ! norm2 := sqrt ( X**T * X)
! !
! Arguments: ! Arguments:
! x - real,dimension(:). The input vector containing the entries of X. ! x(:) - real The input vector containing the entries of X.
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code ! info - integer. Return code
! !
function psb_dnrm2v(x, desc_a, info) function psb_dnrm2v(x, desc_a, info) result(res)
use psb_descriptor_type use psb_descriptor_type
use psb_check_mod use psb_check_mod
use psb_error_mod use psb_error_mod
use psb_penv_mod use psb_penv_mod
implicit none implicit none
real(psb_dpk_), intent(in) :: x(:) real(psb_dpk_), intent(in) :: x(:)
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
real(psb_dpk_) :: psb_dnrm2v real(psb_dpk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm & err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm, ldx
real(psb_dpk_) :: nrm2, dnrm2, dd real(psb_dpk_) :: dnrm2, dd
!!$ external dcombnrm2
!!$ external scombnrm2
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_dnrm2v' name='psb_dnrm2v'
if (psb_errstatus_fatal()) return if(psb_get_errstatus() /= 0) return
info=psb_success_ info=psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -216,8 +210,8 @@ function psb_dnrm2v(x, desc_a, info)
ix = 1 ix = 1
jx=1 jx=1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = size(x,1)
call psb_chkvect(m,1,size(x),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(m,ione,ldx,ix,jx,desc_a,info,iix,jjx)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -230,28 +224,22 @@ function psb_dnrm2v(x, desc_a, info)
goto 9999 goto 9999
end if end if
if(m /= 0) then if (desc_a%get_local_rows() > 0) then
if (desc_a%get_local_rows() > 0) then ndim = desc_a%get_local_rows()
ndim = desc_a%get_local_rows() res = dnrm2( int(ndim,kind=psb_mpik_), x, int(ione,kind=psb_mpik_) )
nrm2 = dnrm2( ndim, x, ione ) ! adjust because overlapped elements are computed more than once
! adjust because overlapped elements are computed more than once do i=1,size(desc_a%ovrlap_elem,1)
do i=1,size(desc_a%ovrlap_elem,1) idx = desc_a%ovrlap_elem(i,1)
idx = desc_a%ovrlap_elem(i,1) ndm = desc_a%ovrlap_elem(i,2)
ndm = desc_a%ovrlap_elem(i,2) dd = real(ndm-1)/real(ndm)
dd = dble(ndm-1)/dble(ndm) res = res * sqrt(sone - dd*(abs(x(idx))/res)**2)
nrm2 = nrm2 * sqrt(done - dd*(abs(x(idx))/nrm2)**2) end do
end do
else
nrm2 = dzero
end if
else else
nrm2 = dzero res = dzero
end if end if
!!$ call pdtreecomb(ictxt,'All',1,nrm2,-1,-1,dcombnrm2) call psb_nrm2(ictxt,res)
call psb_nrm2(ictxt,nrm2)
psb_dnrm2v = nrm2
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -267,6 +255,7 @@ function psb_dnrm2v(x, desc_a, info)
end function psb_dnrm2v end function psb_dnrm2v
function psb_dnrm2_vect(x, desc_a, info) result(res) function psb_dnrm2_vect(x, desc_a, info) result(res)
use psb_descriptor_type use psb_descriptor_type
use psb_check_mod use psb_check_mod
@ -278,14 +267,14 @@ function psb_dnrm2_vect(x, desc_a, info) result(res)
real(psb_dpk_) :: res real(psb_dpk_) :: res
type(psb_d_vect_type), intent (inout) :: x type(psb_d_vect_type), intent (inout) :: x
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm & err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm, ldx
real(psb_dpk_) :: nrm2, dnrm2, dd real(psb_dpk_) :: snrm2, dd
!!$ external dcombnrm2 !!$ external dcombnrm2
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_dnrm2v' name='psb_dnrm2v'
if (psb_errstatus_fatal()) return if (psb_errstatus_fatal()) return
@ -311,8 +300,8 @@ function psb_dnrm2_vect(x, desc_a, info) result(res)
ix = 1 ix = 1
jx=1 jx=1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = x%get_nrows()
call psb_chkvect(m,1,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(m,ione,ldx,ix,jx,desc_a,info,iix,jjx)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -325,10 +314,9 @@ function psb_dnrm2_vect(x, desc_a, info) result(res)
goto 9999 goto 9999
end if end if
if(m /= 0) then if (desc_a%get_local_rows() > 0) then
if (desc_a%get_local_rows() > 0) then ndim = desc_a%get_local_rows()
ndim = desc_a%get_local_rows() res = x%nrm2(ndim)
nrm2 = x%nrm2(ndim)
!!$ ! adjust because overlapped elements are computed more than once !!$ ! adjust because overlapped elements are computed more than once
!!$ do i=1,size(desc_a%ovrlap_elem,1) !!$ do i=1,size(desc_a%ovrlap_elem,1)
!!$ idx = desc_a%ovrlap_elem(i,1) !!$ idx = desc_a%ovrlap_elem(i,1)
@ -336,17 +324,11 @@ function psb_dnrm2_vect(x, desc_a, info) result(res)
!!$ dd = dble(ndm-1)/dble(ndm) !!$ dd = dble(ndm-1)/dble(ndm)
!!$ nrm2 = nrm2 * sqrt(done - dd*(abs(x(idx))/nrm2)**2) !!$ nrm2 = nrm2 * sqrt(done - dd*(abs(x(idx))/nrm2)**2)
!!$ end do !!$ end do
else
nrm2 = dzero
end if
else else
nrm2 = dzero res = dzero
end if end if
!!$ call pdtreecomb(ictxt,'All',1,nrm2,-1,-1,dcombnrm2) call psb_nrm2(ictxt,res)
call psb_nrm2(ictxt,nrm2)
res = nrm2
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -363,8 +345,6 @@ end function psb_dnrm2_vect
!!$ !!$
!!$ Parallel Sparse BLAS version 3.0 !!$ Parallel Sparse BLAS version 3.0
!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010 !!$ (C) Copyright 2006, 2007, 2008, 2009, 2010
@ -400,13 +380,13 @@ end function psb_dnrm2_vect
! Subroutine: psb_dnrm2vs ! Subroutine: psb_dnrm2vs
! Forms the norm2 of a distributed vector, ! Forms the norm2 of a distributed vector,
! !
! res := sqrt ( X**T * X) ! norm2 := sqrt ( X**T * X)
! !
! Arguments: ! Arguments:
! res - real. The result. ! res - real The result.
! x - real,dimension(:). The input vector containing the entries of X. ! x(:) - real The input vector containing the entries of X.
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code ! info - integer. Return code
! !
subroutine psb_dnrm2vs(res, x, desc_a, info) subroutine psb_dnrm2vs(res, x, desc_a, info)
use psb_descriptor_type use psb_descriptor_type
@ -415,20 +395,21 @@ subroutine psb_dnrm2vs(res, x, desc_a, info)
use psb_penv_mod use psb_penv_mod
implicit none implicit none
real(psb_dpk_), intent(in) :: x(:) real(psb_dpk_), intent(in) :: x(:)
real(psb_dpk_), intent(out) :: res real(psb_dpk_), intent(out) :: res
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm & err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm, ldx
real(psb_dpk_) :: nrm2, dnrm2, dd real(psb_dpk_) :: nrm2, dnrm2, dd
!!$ external dcombnrm2
!!$ external scombnrm2
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_dnrm2' name='psb_dnrm2'
if (psb_errstatus_fatal()) return if(psb_get_errstatus() /= 0) return
info=psb_success_ info=psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -444,8 +425,8 @@ subroutine psb_dnrm2vs(res, x, desc_a, info)
ix = 1 ix = 1
jx = 1 jx = 1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = size(x,1)
call psb_chkvect(m,1,size(x),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(m,ione,ldx,ix,jx,desc_a,info,iix,jjx)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -458,28 +439,23 @@ subroutine psb_dnrm2vs(res, x, desc_a, info)
goto 9999 goto 9999
end if end if
if(m /= 0) then if (desc_a%get_local_rows() > 0) then
if (desc_a%get_local_rows() > 0) then ndim = desc_a%get_local_rows()
ndim = desc_a%get_local_rows() res = dnrm2( int(ndim,kind=psb_mpik_), x, int(ione,kind=psb_mpik_) )
nrm2 = dnrm2( ndim, x, ione )
! adjust because overlapped elements are computed more than once ! adjust because overlapped elements are computed more than once
do i=1,size(desc_a%ovrlap_elem,1) do i=1,size(desc_a%ovrlap_elem,1)
idx = desc_a%ovrlap_elem(i,1) idx = desc_a%ovrlap_elem(i,1)
ndm = desc_a%ovrlap_elem(i,2) ndm = desc_a%ovrlap_elem(i,2)
dd = dble(ndm-1)/dble(ndm) dd = real(ndm-1)/real(ndm)
nrm2 = nrm2 * sqrt(done - dd*(abs(x(idx))/nrm2)**2) res = res * sqrt(sone - dd*(abs(x(idx))/res)**2)
end do end do
else
nrm2 = dzero
end if
else else
nrm2 = dzero res = dzero
end if end if
!!$ call pdtreecomb(ictxt,'All',1,nrm2,-1,-1,dcombnrm2) call psb_nrm2(ictxt,res)
call psb_nrm2(ictxt,nrm2)
res = nrm2
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

36
base/psblas/psb_dnrmi.f90
Unescape Escape View File

@ -32,38 +32,31 @@
! File: psb_dnrmi.f90 ! File: psb_dnrmi.f90
! !
! Function: psb_dnrmi ! Function: psb_dnrmi
! Forms the approximated norm of a sparse matrix, ! Forms the infinity norm of a sparse matrix,
! !
! normi := max(sum(abs(A(i,:)))) ! nrmi := max_i(abs(sum(A(i,:))))
! !
! Arguments: ! Arguments:
! a - type(psb_dspmat_type). The sparse matrix containing A. ! a - type(psb_dspmat_type). The sparse matrix containing A.
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code ! info - integer. Return code
! !
function psb_dnrmi(a,desc_a,info) function psb_dnrmi(a,desc_a,info) result(res)
use psb_descriptor_type use psb_base_mod, psb_protect_name => psb_dnrmi
use psb_serial_mod
use psb_check_mod
use psb_error_mod
use psb_penv_mod
use psb_mat_mod
implicit none implicit none
type(psb_dspmat_type), intent(in) :: a type(psb_dspmat_type), intent(in) :: a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
real(psb_dpk_) :: psb_dnrmi real(psb_dpk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, n, iia, jja, ia, ja, mdim, ndim, m & err_act, n, iia, jja, ia, ja, mdim, ndim, m
real(psb_dpk_) :: nrmi character(len=20) :: name, ch_err
character(len=20) :: name, ch_err
name='psb_dnrmi' name='psb_dnrmi'
psb_dnrmi = dzero if(psb_get_errstatus() /= 0) return
if (psb_errstatus_fatal()) return
info=psb_success_ info=psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -96,20 +89,19 @@ function psb_dnrmi(a,desc_a,info)
end if end if
if ((m /= 0).and.(n /= 0)) then if ((m /= 0).and.(n /= 0)) then
nrmi = a%csnmi() res = a%csnmi()
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_csnmi' ch_err='psb_dnrmi'
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
else else
nrmi = 0.d0 res = dzero
end if end if
! compute global max
call psb_amx(ictxt, nrmi)
psb_dnrmi = nrmi ! compute global max
call psb_amx(ictxt, res)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

83
base/psblas/psb_dspmm.f90
Unescape Escape View File

@ -48,18 +48,18 @@
! vectors and A is a M-by-N distributed matrix. ! vectors and A is a M-by-N distributed matrix.
! !
! Arguments: ! Arguments:
! alpha - real The scalar alpha. ! alpha - real The scalar alpha.
! a - type(psb_zspmat_type). The sparse matrix containing A. ! a - type(psb_dspmat_type). The sparse matrix containing A.
! x(:,:) - real The input vector containing the entries of ( X ). ! x(:,:) - real The input vector containing the entries of ( X ).
! beta - real The scalar beta. ! beta - real The scalar beta.
! y(:,:) - real The input vector containing the entries of ( Y ). ! y(:,:) - real The input vector containing the entries of ( Y ).
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code ! info - integer. Return code
! trans - character(optional). Whether A or A'. Default: 'N' ! trans - character(optional). Whether A or A'. Default: 'N'
! k - integer(optional). The number of right-hand sides. ! k - integer(optional). The number of right-hand sides.
! jx - integer(optional). The column offset for ( X ). Default: 1 ! jx - integer(optional). The column offset for ( X ). Default: 1
! jy - integer(optional). The column offset for ( Y ). Default: 1 ! jy - integer(optional). The column offset for ( Y ). Default: 1
! work(:) - real,(optional). Working area. ! work(:) - real,(optional). Working area.
! doswap - logical(optional). Whether to performe halo updates. ! doswap - logical(optional). Whether to performe halo updates.
! !
subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,& subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
@ -84,16 +84,16 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
& err_act, n, iix, jjx, ia, ja, iia, jja, ix, iy, ik, ijx, ijy,& & err_act, n, iix, jjx, ia, ja, iia, jja, ix, iy, ik, ijx, ijy,&
& m, nrow, ncol, lldx, lldy, liwork, iiy, jjy,& & m, nrow, ncol, lldx, lldy, liwork, iiy, jjy,&
& i, ib, ib1, ip, idx & i, ib, ib1, ip, idx
integer(psb_ipk_), parameter :: nb=4 integer(psb_ipk_), parameter :: nb=4
real(psb_dpk_), pointer :: xp(:,:), yp(:,:), iwork(:) real(psb_dpk_), pointer :: xp(:,:), yp(:,:), iwork(:)
real(psb_dpk_), allocatable :: xvsave(:,:) real(psb_dpk_), allocatable :: xvsave(:,:)
character :: trans_ character :: trans_
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
logical :: aliw, doswap_ logical :: aliw, doswap_
integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: debug_level, debug_unit
name='psb_dspmm' name='psb_dspmm'
if (psb_errstatus_fatal()) return if(psb_get_errstatus() /= 0) return
info=psb_success_ info=psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit() debug_unit = psb_get_debug_unit()
@ -205,9 +205,9 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
end if end if
! checking for vectors correctness ! checking for vectors correctness
call psb_chkvect(n,ik,size(x,1),ix,ijx,desc_a,info,iix,jjx) call psb_chkvect(n,ik,lldx,ix,ijx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,size(y,1),iy,ijy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ik,lldy,iy,ijy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -277,9 +277,9 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
! checking for vectors correctness ! checking for vectors correctness
call psb_chkvect(m,ik,size(x,1),ix,ijx,desc_a,info,iix,jjx) call psb_chkvect(m,ik,lldx,ix,ijx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(n,ik,size(y,1),iy,ijy,desc_a,info,iiy,jjy) & call psb_chkvect(n,ik,lldy,iy,ijy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -304,7 +304,7 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
if (info == psb_success_) call psi_ovrl_upd(x,desc_a,psb_avg_,info) if (info == psb_success_) call psi_ovrl_upd(x,desc_a,psb_avg_,info)
y(nrow+1:ncol,1:ik) = dzero y(nrow+1:ncol,1:ik) = dzero
if (info == psb_success_)& if (info == psb_success_) &
& call psb_csmm(alpha,a,x(:,1:ik),beta,y(:,1:ik),info,trans=trans_) & call psb_csmm(alpha,a,x(:,1:ik),beta,y(:,1:ik),info,trans=trans_)
if (debug_level >= psb_debug_comp_) & if (debug_level >= psb_debug_comp_) &
& write(debug_unit,*) me,' ',trim(name),' csmm ', info & write(debug_unit,*) me,' ',trim(name),' csmm ', info
@ -404,15 +404,15 @@ end subroutine psb_dspmm
! vectors and A is a M-by-N distributed matrix. ! vectors and A is a M-by-N distributed matrix.
! !
! Arguments: ! Arguments:
! alpha - real The scalar alpha. ! alpha - real The scalar alpha.
! a - type(psb_zspmat_type). The sparse matrix containing A. ! a - type(psb_dspmat_type). The sparse matrix containing A.
! x(:) - real The input vector containing the entries of ( X ). ! x(:) - real The input vector containing the entries of ( X ).
! beta - real The scalar beta. ! beta - real The scalar beta.
! y(:) - real The input vector containing the entries of ( Y ). ! y(:) - real The input vector containing the entries of ( Y ).
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code ! info - integer. Return code
! trans - character(optional). Whether A or A'. Default: 'N' ! trans - character(optional). Whether A or A'. Default: 'N'
! work(:) - real,(optional). Working area. ! work(:) - real,(optional). Working area.
! doswap - logical(optional). Whether to performe halo updates. ! doswap - logical(optional). Whether to performe halo updates.
! !
subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,& subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
@ -427,7 +427,7 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
type(psb_dspmat_type), intent(in) :: a type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
real(psb_dpk_), optional, target, intent(inout) :: work(:) real(psb_dpk_), optional, target, intent(inout) :: work(:)
character, intent(in), optional :: trans character, intent(in), optional :: trans
logical, intent(in), optional :: doswap logical, intent(in), optional :: doswap
@ -436,16 +436,16 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
& err_act, n, iix, jjx, ia, ja, iia, jja, ix, iy, ik, & & err_act, n, iix, jjx, ia, ja, iia, jja, ix, iy, ik, &
& m, nrow, ncol, lldx, lldy, liwork, jx, jy, iiy, jjy,& & m, nrow, ncol, lldx, lldy, liwork, jx, jy, iiy, jjy,&
& ib, ip, idx & ib, ip, idx
integer(psb_ipk_), parameter :: nb=4 integer(psb_ipk_), parameter :: nb=4
real(psb_dpk_), pointer :: iwork(:), xp(:), yp(:) real(psb_dpk_), pointer :: iwork(:), xp(:), yp(:)
real(psb_dpk_), allocatable :: xvsave(:) real(psb_dpk_), allocatable :: xvsave(:)
character :: trans_ character :: trans_
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
logical :: aliw, doswap_ logical :: aliw, doswap_
integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: debug_level, debug_unit
name='psb_dspmv' name='psb_dspmv'
if (psb_errstatus_fatal()) return if(psb_get_errstatus() /= 0) return
info=psb_success_ info=psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit() debug_unit = psb_get_debug_unit()
@ -543,9 +543,9 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
end if end if
! checking for vectors correctness ! checking for vectors correctness
call psb_chkvect(n,ik,size(x),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(n,ik,lldx,ix,jx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,size(y),iy,jy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ik,lldy,iy,jy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -583,9 +583,9 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
end if end if
! checking for vectors correctness ! checking for vectors correctness
call psb_chkvect(m,ik,size(x),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(m,ik,lldx,ix,jx,desc_a,info,iix,jjx)
if (info == psb_success_)& if (info == psb_success_)&
& call psb_chkvect(n,ik,size(y),iy,jy,desc_a,info,iiy,jjy) & call psb_chkvect(n,ik,lldy,iy,jy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -682,13 +682,13 @@ subroutine psb_dspmv_vect(alpha,a,x,beta,y,desc_a,info,&
use psi_mod use psi_mod
implicit none implicit none
real(psb_dpk_), intent(in) :: alpha, beta real(psb_dpk_), intent(in) :: alpha, beta
type(psb_d_vect_type), intent(inout) :: x type(psb_d_vect_type), intent(inout) :: x
type(psb_d_vect_type), intent(inout) :: y type(psb_d_vect_type), intent(inout) :: y
type(psb_dspmat_type), intent(in) :: a type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
real(psb_dpk_), optional, target, intent(inout) :: work(:) real(psb_dpk_), optional, target, intent(inout) :: work(:)
character, intent(in), optional :: trans character, intent(in), optional :: trans
logical, intent(in), optional :: doswap logical, intent(in), optional :: doswap
@ -816,9 +816,9 @@ subroutine psb_dspmv_vect(alpha,a,x,beta,y,desc_a,info,&
end if end if
! checking for vectors correctness ! checking for vectors correctness
call psb_chkvect(n,ik,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(n,ik,lldx,ix,jx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,y%get_nrows(),iy,jy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ik,lldy,iy,jy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -856,9 +856,9 @@ subroutine psb_dspmv_vect(alpha,a,x,beta,y,desc_a,info,&
end if end if
! checking for vectors correctness ! checking for vectors correctness
call psb_chkvect(m,ik,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(m,ik,lldx,ix,jx,desc_a,info,iix,jjx)
if (info == psb_success_)& if (info == psb_success_)&
& call psb_chkvect(n,ik,y%get_nrows(),iy,jy,desc_a,info,iiy,jjy) & call psb_chkvect(n,ik,lldy,iy,jy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -882,6 +882,7 @@ subroutine psb_dspmv_vect(alpha,a,x,beta,y,desc_a,info,&
! !
call psi_ovrl_save(x%v,xvsave,desc_a,info) call psi_ovrl_save(x%v,xvsave,desc_a,info)
if (info == psb_success_) call psi_ovrl_upd(x%v,desc_a,psb_avg_,info) if (info == psb_success_) call psi_ovrl_upd(x%v,desc_a,psb_avg_,info)
!!! THIS SHOULD BE FIXED !!! But beta is almost never /= 0 !!! THIS SHOULD BE FIXED !!! But beta is almost never /= 0
!!$ yp(nrow+1:ncol) = dzero !!$ yp(nrow+1:ncol) = dzero
@ -909,7 +910,7 @@ subroutine psb_dspmv_vect(alpha,a,x,beta,y,desc_a,info,&
& write(debug_unit,*) me,' ',trim(name),' swaptran ', info & write(debug_unit,*) me,' ',trim(name),' swaptran ', info
if(info /= psb_success_) then if(info /= psb_success_) then
info = psb_err_from_subroutine_ info = psb_err_from_subroutine_
ch_err='PSI_dSwapTran' ch_err='PSI_SwapTran'
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if

34
base/psblas/psb_dspnrm1.f90
Unescape Escape View File

@ -29,42 +29,34 @@
!!$ POSSIBILITY OF SUCH DAMAGE. !!$ POSSIBILITY OF SUCH DAMAGE.
!!$ !!$
!!$ !!$
! File: psb_dnrmi.f90 ! File: psb_dnrm1.f90
! !
! Function: psb_dnrmi ! Function: psb_dnrm1
! Forms the approximated norm of a sparse matrix, ! Forms the norm1 of a sparse matrix,
! !
! norm1 := max(sum(abs(A(:,j)))) ! norm1 := max_j(sum(abs(A(:,j))))
! !
! Arguments: ! Arguments:
! a - type(psb_dspmat_type). The sparse matrix containing A. ! a - type(psb_dspmat_type). The sparse matrix containing A.
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code ! info - integer. Return code
! !
function psb_dspnrm1(a,desc_a,info) function psb_dspnrm1(a,desc_a,info) result(res)
!!$ use psb_descriptor_type
!!$ use psb_serial_mod
!!$ use psb_check_mod
!!$ use psb_error_mod
!!$ use psb_penv_mod
!!$ use psb_mat_mod
!!$ use psb_tools_mod
use psb_base_mod, psb_protect_name => psb_dspnrm1 use psb_base_mod, psb_protect_name => psb_dspnrm1
implicit none implicit none
type(psb_dspmat_type), intent(in) :: a type(psb_dspmat_type), intent(in) :: a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
real(psb_dpk_) :: psb_dspnrm1 real(psb_dpk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me, nr,nc,& integer(psb_ipk_) :: ictxt, np, me, nr,nc,&
& err_act, n, iia, jja, ia, ja, mdim, ndim, m & err_act, n, iia, jja, ia, ja, mdim, ndim, m
real(psb_dpk_) :: nrm1
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
real(psb_dpk_), allocatable :: v(:) real(psb_dpk_), allocatable :: v(:)
name='psb_dnrm1' name='psb_dspnrm1'
if (psb_errstatus_fatal()) return if (psb_errstatus_fatal()) return
info=psb_success_ info=psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -120,14 +112,12 @@ function psb_dspnrm1(a,desc_a,info)
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
nrm1 = maxval(v(1:nr)) res = maxval(v(1:nr))
else else
nrm1 = 0.d0 res = dzero
end if end if
! compute global max ! compute global max
call psb_amx(ictxt, nrm1) call psb_amx(ictxt, res)
psb_dspnrm1 = nrm1
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

75
base/psblas/psb_dspsm.f90
Unescape Escape View File

@ -56,23 +56,22 @@
! vector and T is a M-by-M distributed triangular matrix. ! vector and T is a M-by-M distributed triangular matrix.
! !
! Arguments: ! Arguments:
! alpha - real The scalar alpha. ! alpha - real. The scalar alpha.
! a - type(psb_zspmat_type). The sparse matrix containing A. ! a - type(psb_dspmat_type). The sparse matrix containing A.
! x(:,:) - real The input vector containing the entries of ( X ). ! x(:,:) - real The input vector containing the entries of ( X ).
! beta - real The scalar beta. ! beta - real The scalar beta.
! y(:,:) - real The input vector containing the entries of ( Y ). ! y(:,:) - real The input vector containing the entries of ( Y ).
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code ! info - integer. Return code
! trans - character(optional). Whether A or A'. If not present 'N' is assumed. ! trans - character(optional). Whether A or A'. If not present 'N' is assumed.
! scale - character(optional). Specify some type of operation with ! scale - character(optional). Specify some type of operation with
! the diagonal matrix D. ! the diagonal matrix D.
! choice - integer(optional). The kind of update to perform on overlap elements. ! choice - integer(optional). The kind of update to perform on overlap elements.
! d(:) - real , optional Matrix for diagonal scaling. ! d(:) - real, optional Matrix for diagonal scaling.
! k - integer(optional). The number of right-hand sides. ! k - integer(optional). The number of right-hand sides.
! jx - integer(optional). The column offset for ( X ). Default: 1 ! jx - integer(optional). The column offset for ( X ). Default: 1
! jy - integer(optional). The column offset for ( Y ). Default: 1 ! jy - integer(optional). The column offset for ( Y ). Default: 1
! work(:) - real , optional Working area. ! work(:) - real, optional Working area.
!
! !
subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,& subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,&
& trans, scale, choice, diag, k, jx, jy, work) & trans, scale, choice, diag, k, jx, jy, work)
@ -83,11 +82,11 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,&
real(psb_dpk_), intent(in) :: alpha, beta real(psb_dpk_), intent(in) :: alpha, beta
real(psb_dpk_), intent(in), target :: x(:,:) real(psb_dpk_), intent(in), target :: x(:,:)
real(psb_dpk_), intent(inout), target :: y(:,:) real(psb_dpk_), intent(inout), target :: y(:,:)
type(psb_dspmat_type), intent(in) :: a type (psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
real(psb_dpk_), intent(in), optional, target :: diag(:) real(psb_dpk_), intent(in), optional, target :: diag(:)
real(psb_dpk_), optional, target, intent(inout) :: work(:) real(psb_dpk_), optional, target, intent(inout) :: work(:)
character, intent(in), optional :: trans, scale character, intent(in), optional :: trans, scale
integer(psb_ipk_), intent(in), optional :: choice integer(psb_ipk_), intent(in), optional :: choice
integer(psb_ipk_), intent(in), optional :: k, jx, jy integer(psb_ipk_), intent(in), optional :: k, jx, jy
@ -99,14 +98,14 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,&
& m, nrow, ncol, liwork, llwork, iiy, jjy, idx, ndm & m, nrow, ncol, liwork, llwork, iiy, jjy, idx, ndm
character :: lscale character :: lscale
integer(psb_ipk_), parameter :: nb=4 integer(psb_ipk_), parameter :: nb=4
real(psb_dpk_),pointer :: iwork(:), xp(:,:), yp(:,:), id(:) real(psb_dpk_),pointer :: iwork(:), xp(:,:), yp(:,:), id(:)
character :: itrans character :: itrans
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
logical :: aliw logical :: aliw
name='psb_dspsm' name='psb_dspsm'
if (psb_errstatus_fatal()) return if(psb_get_errstatus() /= 0) return
info=psb_success_ info=psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -221,9 +220,9 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,&
call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja) call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja)
! checking for vectors correctness ! checking for vectors correctness
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,size(x,1),ix,ijx,desc_a,info,iix,jjx) & call psb_chkvect(m,ik,lldx,ix,ijx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,size(y,1),iy,ijy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ik,lldy,iy,ijy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect/mat' ch_err='psb_chkvect/mat'
@ -339,19 +338,19 @@ end subroutine psb_dspsm
! !
! !
! Arguments: ! Arguments:
! alpha - real The scalar alpha. ! alpha - real. The scalar alpha.
! a - type(psb_zspmat_type). The sparse matrix containing A. ! a - type(psb_dspmat_type). The sparse matrix containing A.
! x(:) - real The input vector containing the entries of ( X ). ! x(:) - real The input vector containing the entries of ( X ).
! beta - real The scalar beta. ! beta - real The scalar beta.
! y(:) - real The input vector containing the entries of ( Y ). ! y(:) - real The input vector containing the entries of ( Y ).
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code ! info - integer. Return code
! trans - character(optional). Whether A or A'. If not present 'N' is assumed. ! trans - character(optional). Whether A or A'. If not present 'N' is assumed.
! scale - character(optional). Specify some type of operation with ! scale - character(optional). Specify some type of operation with
! the diagonal matrix D. ! the diagonal matrix D.
! choice - integer(optional). The kind of update to perform on overlap elements. ! choice - integer(optional). The kind of update to perform on overlap elements.
! d(:) - real , optional Matrix for diagonal scaling. ! d(:) - real, optional Matrix for diagonal scaling.
! work(:) - real , optional Working area. ! work(:) - real, optional Working area.
! !
subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,& subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,&
& trans, scale, choice, diag, work) & trans, scale, choice, diag, work)
@ -366,7 +365,7 @@ subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,&
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
real(psb_dpk_), intent(in), optional, target :: diag(:) real(psb_dpk_), intent(in), optional, target :: diag(:)
real(psb_dpk_), optional, target, intent(inout) :: work(:) real(psb_dpk_), optional, target, intent(inout) :: work(:)
character, intent(in), optional :: trans, scale character, intent(in), optional :: trans, scale
integer(psb_ipk_), intent(in), optional :: choice integer(psb_ipk_), intent(in), optional :: choice
@ -384,7 +383,7 @@ subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,&
logical :: aliw logical :: aliw
name='psb_dspsv' name='psb_dspsv'
if (psb_errstatus_fatal()) return if(psb_get_errstatus() /= 0) return
info=psb_success_ info=psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -471,7 +470,7 @@ subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,&
iwork(1)=0.d0 iwork(1)=0.d0
if (present(diag)) then if(present(diag)) then
lld = size(diag) lld = size(diag)
id => diag id => diag
else else
@ -484,9 +483,9 @@ subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,&
call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja) call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja)
! checking for vectors correctness ! checking for vectors correctness
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,size(x),ix,jx,desc_a,info,iix,jjx) & call psb_chkvect(m,ik,lldx,ix,jx,desc_a,info,iix,jjx)
if (info == psb_success_)& if (info == psb_success_) &
& call psb_chkvect(m,ik,size(y),iy,jy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ik,lldy,iy,jy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect/mat' ch_err='psb_chkvect/mat'
@ -522,7 +521,7 @@ subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,&
end if end if
! update overlap elements ! update overlap elements
if (choice_ > 0) then if(choice_ > 0) then
call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),& call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),&
& done,yp,desc_a,iwork,info,data=psb_comm_ovr_) & done,yp,desc_a,iwork,info,data=psb_comm_ovr_)
@ -550,22 +549,23 @@ subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,&
return return
end subroutine psb_dspsv end subroutine psb_dspsv
subroutine psb_dspsv_vect(alpha,a,x,beta,y,desc_a,info,& subroutine psb_dspsv_vect(alpha,a,x,beta,y,desc_a,info,&
& trans, scale, choice, diag, work) & trans, scale, choice, diag, work)
use psb_base_mod, psb_protect_name => psb_dspsv_vect use psb_base_mod, psb_protect_name => psb_dspsv_vect
use psi_mod use psi_mod
implicit none implicit none
real(psb_dpk_), intent(in) :: alpha, beta real(psb_dpk_), intent(in) :: alpha, beta
type(psb_d_vect_type), intent(inout) :: x type(psb_d_vect_type), intent(inout) :: x
type(psb_d_vect_type), intent(inout) :: y type(psb_d_vect_type), intent(inout) :: y
type(psb_dspmat_type), intent(inout) :: a type(psb_dspmat_type), intent(inout) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
type(psb_d_vect_type), intent(inout), optional :: diag type(psb_d_vect_type), intent(inout), optional :: diag
real(psb_dpk_), optional, target, intent(inout) :: work(:) real(psb_dpk_), optional, target, intent(inout) :: work(:)
character, intent(in), optional :: trans, scale character, intent(in), optional :: trans, scale
integer(psb_ipk_), intent(in), optional :: choice integer(psb_ipk_), intent(in), optional :: choice
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me, & integer(psb_ipk_) :: ictxt, np, me, &
@ -580,7 +580,7 @@ subroutine psb_dspsv_vect(alpha,a,x,beta,y,desc_a,info,&
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
logical :: aliw logical :: aliw
name='psb_dspsv' name='psb_sspsv'
if (psb_errstatus_fatal()) return if (psb_errstatus_fatal()) return
info=psb_success_ info=psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -684,9 +684,9 @@ subroutine psb_dspsv_vect(alpha,a,x,beta,y,desc_a,info,&
call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja) call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja)
! checking for vectors correctness ! checking for vectors correctness
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) & call psb_chkvect(m,ik,lldx,ix,jx,desc_a,info,iix,jjx)
if (info == psb_success_)& if (info == psb_success_)&
& call psb_chkvect(m,ik,y%get_nrows(),iy,jy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ik,lldy,iy,jy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect/mat' ch_err='psb_chkvect/mat'
@ -710,11 +710,6 @@ subroutine psb_dspsv_vect(alpha,a,x,beta,y,desc_a,info,&
end if end if
! Perform local triangular system solve ! Perform local triangular system solve
!!$ if (present(diag)) then
!!$ call psb_cssm(alpha,a,x,beta,y,info,scale=scale,d=diag,trans=trans)
!!$ else
!!$ call psb_cssm(alpha,a,x,beta,y,info,scale=scale,trans=trans)
!!$ end if
if (present(diag)) then if (present(diag)) then
call a%cssm(alpha,x,beta,y,info,scale=scale,d=diag,trans=trans) call a%cssm(alpha,x,beta,y,info,scale=scale,d=diag,trans=trans)
else else

70
base/psblas/psb_saxpby.f90
Unescape Escape View File

@ -37,7 +37,7 @@ subroutine psb_saxpby_vect(alpha, x, beta, y,&
implicit none implicit none
type(psb_s_vect_type), intent (inout) :: x type(psb_s_vect_type), intent (inout) :: x
type(psb_s_vect_type), intent (inout) :: y type(psb_s_vect_type), intent (inout) :: y
real(psb_spk_), intent (in) :: alpha, beta real(psb_spk_), intent (in) :: alpha, beta
type(psb_desc_type), intent (in) :: desc_a type(psb_desc_type), intent (in) :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -115,6 +115,7 @@ subroutine psb_saxpby_vect(alpha, x, beta, y,&
return return
end subroutine psb_saxpby_vect end subroutine psb_saxpby_vect
! !
! Subroutine: psb_saxpby ! Subroutine: psb_saxpby
! Adds one distributed matrix to another, ! Adds one distributed matrix to another,
@ -126,17 +127,18 @@ end subroutine psb_saxpby_vect
! sub( Y ) denotes Y(:,JY). ! sub( Y ) denotes Y(:,JY).
! !
! Arguments: ! Arguments:
! alpha - real The scalar alpha ! alpha - real,input The scalar used to multiply each component of X
! x(:,:) - real The input vector containing the entries of ( X ). ! x(:,:) - real,input The input vector containing the entries of X
! beta - real The scalar beta ! beta - real,input The scalar used to multiply each component of Y
! y(:,:) - real The input vector containing the entries of ( Y ). ! y(:,:) - real,inout The input vector Y
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type) The communication descriptor.
! info - integer. Return code ! info - integer Return code
! jx - integer(optional). The column offset for sub( X ). ! jx - integer(optional) The column offset for X
! jy - integer(optional). The column offset for sub( Y ). ! jy - integer(optional) The column offset for Y
! !
subroutine psb_saxpby(alpha, x, beta,y,desc_a,info, n, jx, jy) subroutine psb_saxpby(alpha, x, beta,y,desc_a,info, n, jx, jy)
use psb_base_mod, psb_protect_name => psb_saxpby use psb_base_mod, psb_protect_name => psb_saxpby
implicit none implicit none
integer(psb_ipk_), intent(in), optional :: n, jx, jy integer(psb_ipk_), intent(in), optional :: n, jx, jy
@ -148,16 +150,16 @@ subroutine psb_saxpby(alpha, x, beta,y,desc_a,info, n, jx, jy)
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ix, iy, ijx, ijy, m, iiy, in, jjy & err_act, iix, jjx, ix, iy, ijx, ijy, m, iiy, in, jjy, &
& lldx, lldy
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_sgeaxpby' name='psb_geaxpby'
if(psb_get_errstatus() /= 0) return if(psb_get_errstatus() /= 0) return
info=psb_success_ info=psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if (np == -ione) then if (np == -ione) then
info = psb_err_context_error_ info = psb_err_context_error_
@ -197,11 +199,12 @@ subroutine psb_saxpby(alpha, x, beta,y,desc_a,info, n, jx, jy)
end if end if
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
lldx = size(x,1)
lldy = size(y,1)
! check vector correctness ! check vector correctness
call psb_chkvect(m,ione,size(x,1),ix,ijx,desc_a,info,iix,jjx) call psb_chkvect(m,ione,lldx,ix,ijx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ione,size(y,1),iy,ijy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ione,lldy,iy,ijy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -217,9 +220,9 @@ subroutine psb_saxpby(alpha, x, beta,y,desc_a,info, n, jx, jy)
if ((in /= 0)) then if ((in /= 0)) then
if(desc_a%get_local_rows() > 0) then if(desc_a%get_local_rows() > 0) then
call saxpby(desc_a%get_local_rows(),in,& call caxpby(desc_a%get_local_cols(),in,&
& alpha,x(iix:,jjx),size(x,1),beta,& & alpha,x(iix:,jjx),lldx,beta,&
& y(iiy:,jjy),size(y,1),info) & y(iiy:,jjy),lldy,info)
end if end if
end if end if
@ -278,12 +281,12 @@ end subroutine psb_saxpby
! Y := beta * Y + alpha * X ! Y := beta * Y + alpha * X
! !
! Arguments: ! Arguments:
! alpha - real The scalar alpha ! alpha - real,input The scalar used to multiply each component of X
! x(:) - real The input vector containing the entries of ( X ). ! x(:) - real,input The input vector containing the entries of X
! beta - real The scalar beta ! beta - real,input The scalar used to multiply each component of Y
! y(:) - real The input vector containing the entries of ( Y ). ! y(:) - real,inout The input vector Y
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type) The communication descriptor.
! info - integer. Return code ! info - integer Return code
! !
! !
subroutine psb_saxpbyv(alpha, x, beta,y,desc_a,info) subroutine psb_saxpbyv(alpha, x, beta,y,desc_a,info)
@ -298,10 +301,12 @@ subroutine psb_saxpbyv(alpha, x, beta,y,desc_a,info)
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ix, iy, m, iiy, jjy & err_act, iix, jjx, ix, iy, m, iiy, jjy, &
& lldx, lldy
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
logical, parameter :: debug=.false.
name='psb_sgeaxpby' name='psb_geaxpby'
if(psb_get_errstatus() /= 0) return if(psb_get_errstatus() /= 0) return
info=psb_success_ info=psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -319,16 +324,17 @@ subroutine psb_saxpbyv(alpha, x, beta,y,desc_a,info)
iy = ione iy = ione
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
lldx = size(x,1)
lldy = size(y,1)
! check vector correctness ! check vector correctness
call psb_chkvect(m,ione,size(x),ix,ione,desc_a,info,iix,jjx) call psb_chkvect(m,ione,lldx,ix,ione,desc_a,info,iix,jjx)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect 1' ch_err='psb_chkvect 1'
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
call psb_chkvect(m,ione,size(y),iy,ione,desc_a,info,iiy,jjy) call psb_chkvect(m,ione,lldy,iy,ione,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect 2' ch_err='psb_chkvect 2'
@ -342,9 +348,9 @@ subroutine psb_saxpbyv(alpha, x, beta,y,desc_a,info)
end if end if
if(desc_a%get_local_rows() > 0) then if(desc_a%get_local_rows() > 0) then
call saxpby(desc_a%get_local_rows(),ione,& call caxpby(desc_a%get_local_cols(),ione,&
& alpha,x,size(x),beta,& & alpha,x,lldx,beta,&
& y,size(y),info) & y,lldy,info)
end if end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)

197
base/psblas/psb_snrm2.f90
Unescape Escape View File

@ -39,12 +39,12 @@
! where sub( X ) denotes X(:,JX). ! where sub( X ) denotes X(:,JX).
! !
! Arguments: ! Arguments:
! x - real,dimension(:,:). The input vector containing the entries of X. ! x(:,:) - real The input vector containing the entries of sub( X ).
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code ! info - integer. Return code
! jx - integer(optional). The column offset for X . ! jx - integer(optional). The column offset for sub( X ).
! !
function psb_snrm2(x, desc_a, info, jx) function psb_snrm2(x, desc_a, info, jx) result(res)
use psb_descriptor_type use psb_descriptor_type
use psb_check_mod use psb_check_mod
use psb_error_mod use psb_error_mod
@ -55,14 +55,13 @@ function psb_snrm2(x, desc_a, info, jx)
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(in), optional :: jx integer(psb_ipk_), intent(in), optional :: jx
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
real(psb_spk_) :: psb_snrm2 real(psb_spk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ndim, ix, ijx, i, m, id, idx, ndm & err_act, iix, jjx, ndim, ix, ijx, i, m, id, idx, ndm, ldx
real(psb_spk_) :: nrm2, snrm2, dd real(psb_spk_) :: snrm2, dd
!!$ external scombnrm2 character(len=20) :: name, ch_err
character(len=20) :: name, ch_err
name='psb_snrm2' name='psb_snrm2'
if(psb_get_errstatus() /= 0) return if(psb_get_errstatus() /= 0) return
@ -86,7 +85,8 @@ function psb_snrm2(x, desc_a, info, jx)
endif endif
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
call psb_chkvect(m,1,size(x,1),ix,ijx,desc_a,info,iix,jjx) ldx = size(x,1)
call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -99,29 +99,22 @@ function psb_snrm2(x, desc_a, info, jx)
goto 9999 goto 9999
end if end if
if(m /= 0) then if (desc_a%get_local_rows() > 0) then
if (desc_a%get_local_rows() > 0) then ndim = desc_a%get_local_rows()
ndim = desc_a%get_local_rows() res = snrm2( int(ndim,kind=psb_mpik_), x(iix:,jjx), int(ione,kind=psb_mpik_) )
nrm2 = snrm2( ndim, x(iix:,jjx), ione )
! adjust because overlapped elements are computed more than once
! adjust because overlapped elements are computed more than once do i=1,size(desc_a%ovrlap_elem,1)
do i=1,size(desc_a%ovrlap_elem,1) idx = desc_a%ovrlap_elem(i,1)
idx = desc_a%ovrlap_elem(i,1) ndm = desc_a%ovrlap_elem(i,2)
ndm = desc_a%ovrlap_elem(i,2) dd = real(ndm-1)/real(ndm)
dd = real(ndm-1)/real(ndm) res = res * sqrt(sone - dd*(abs(x(idx,jjx))/res)**2)
nrm2 = nrm2 * sqrt(sone - dd*(abs(x(idx,jjx))/nrm2)**2) end do
end do
else
nrm2 = szero
end if
else else
nrm2 = szero res = szero
end if end if
!!$ call pstreecomb(ictxt,'All',1,nrm2,-1,-1,scombnrm2) call psb_nrm2(ictxt,res)
call psb_nrm2(ictxt,nrm2)
psb_snrm2 = nrm2
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -176,28 +169,29 @@ end function psb_snrm2
! norm2 := sqrt ( X**T * X) ! norm2 := sqrt ( X**T * X)
! !
! Arguments: ! Arguments:
! x - real,dimension(:). The input vector containing the entries of X. ! x(:) - real The input vector containing the entries of X.
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code ! info - integer. Return code
! !
function psb_snrm2v(x, desc_a, info) function psb_snrm2v(x, desc_a, info) result(res)
use psb_descriptor_type use psb_descriptor_type
use psb_check_mod use psb_check_mod
use psb_error_mod use psb_error_mod
use psb_penv_mod use psb_penv_mod
implicit none implicit none
real(psb_spk_), intent(in) :: x(:) real(psb_spk_), intent(in) :: x(:)
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
real(psb_spk_) :: psb_snrm2v real(psb_spk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm & err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm, ldx
real(psb_spk_) :: nrm2, snrm2, dd real(psb_spk_) :: snrm2, dd
!!$ external scombnrm2 !!$ external scombnrm2
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_snrm2v' name='psb_snrm2v'
if(psb_get_errstatus() /= 0) return if(psb_get_errstatus() /= 0) return
@ -216,8 +210,8 @@ function psb_snrm2v(x, desc_a, info)
ix = 1 ix = 1
jx=1 jx=1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = size(x,1)
call psb_chkvect(m,1,size(x),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(m,ione,ldx,ix,jx,desc_a,info,iix,jjx)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -230,28 +224,22 @@ function psb_snrm2v(x, desc_a, info)
goto 9999 goto 9999
end if end if
if(m /= 0) then if (desc_a%get_local_rows() > 0) then
if (desc_a%get_local_rows() > 0) then ndim = desc_a%get_local_rows()
ndim = desc_a%get_local_rows() res = snrm2( int(ndim,kind=psb_mpik_), x, int(ione,kind=psb_mpik_) )
nrm2 = snrm2( ndim, x, ione ) ! adjust because overlapped elements are computed more than once
! adjust because overlapped elements are computed more than once do i=1,size(desc_a%ovrlap_elem,1)
do i=1,size(desc_a%ovrlap_elem,1) idx = desc_a%ovrlap_elem(i,1)
idx = desc_a%ovrlap_elem(i,1) ndm = desc_a%ovrlap_elem(i,2)
ndm = desc_a%ovrlap_elem(i,2) dd = real(ndm-1)/real(ndm)
dd = real(ndm-1)/real(ndm) res = res * sqrt(sone - dd*(abs(x(idx))/res)**2)
nrm2 = nrm2 * sqrt(sone - dd*(abs(x(idx))/nrm2)**2) end do
end do
else
nrm2 = szero
end if
else else
nrm2 = szero res = szero
end if end if
!!$ call pstreecomb(ictxt,'All',1,nrm2,-1,-1,scombnrm2) call psb_nrm2(ictxt,res)
call psb_nrm2(ictxt,nrm2)
psb_snrm2v = nrm2
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -267,6 +255,7 @@ function psb_snrm2v(x, desc_a, info)
end function psb_snrm2v end function psb_snrm2v
function psb_snrm2_vect(x, desc_a, info) result(res) function psb_snrm2_vect(x, desc_a, info) result(res)
use psb_descriptor_type use psb_descriptor_type
use psb_check_mod use psb_check_mod
@ -278,12 +267,12 @@ function psb_snrm2_vect(x, desc_a, info) result(res)
real(psb_spk_) :: res real(psb_spk_) :: res
type(psb_s_vect_type), intent (inout) :: x type(psb_s_vect_type), intent (inout) :: x
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm & err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm, ldx
real(psb_spk_) :: nrm2, snrm2, dd real(psb_spk_) :: snrm2, dd
!!$ external dcombnrm2 !!$ external dcombnrm2
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
@ -311,8 +300,8 @@ function psb_snrm2_vect(x, desc_a, info) result(res)
ix = 1 ix = 1
jx=1 jx=1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = x%get_nrows()
call psb_chkvect(m,1,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(m,ione,ldx,ix,jx,desc_a,info,iix,jjx)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -325,10 +314,9 @@ function psb_snrm2_vect(x, desc_a, info) result(res)
goto 9999 goto 9999
end if end if
if(m /= 0) then if (desc_a%get_local_rows() > 0) then
if (desc_a%get_local_rows() > 0) then ndim = desc_a%get_local_rows()
ndim = desc_a%get_local_rows() res = x%nrm2(ndim)
nrm2 = x%nrm2(ndim)
!!$ ! adjust because overlapped elements are computed more than once !!$ ! adjust because overlapped elements are computed more than once
!!$ do i=1,size(desc_a%ovrlap_elem,1) !!$ do i=1,size(desc_a%ovrlap_elem,1)
!!$ idx = desc_a%ovrlap_elem(i,1) !!$ idx = desc_a%ovrlap_elem(i,1)
@ -336,17 +324,11 @@ function psb_snrm2_vect(x, desc_a, info) result(res)
!!$ dd = dble(ndm-1)/dble(ndm) !!$ dd = dble(ndm-1)/dble(ndm)
!!$ nrm2 = nrm2 * sqrt(done - dd*(abs(x(idx))/nrm2)**2) !!$ nrm2 = nrm2 * sqrt(done - dd*(abs(x(idx))/nrm2)**2)
!!$ end do !!$ end do
else
nrm2 = szero
end if
else else
nrm2 = szero res = szero
end if end if
!!$ call pdtreecomb(ictxt,'All',1,nrm2,-1,-1,dcombnrm2) call psb_nrm2(ictxt,res)
call psb_nrm2(ictxt,nrm2)
res = nrm2
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -363,7 +345,6 @@ end function psb_snrm2_vect
!!$ !!$
!!$ Parallel Sparse BLAS version 3.0 !!$ Parallel Sparse BLAS version 3.0
!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010 !!$ (C) Copyright 2006, 2007, 2008, 2009, 2010
@ -399,13 +380,13 @@ end function psb_snrm2_vect
! Subroutine: psb_snrm2vs ! Subroutine: psb_snrm2vs
! Forms the norm2 of a distributed vector, ! Forms the norm2 of a distributed vector,
! !
! res := sqrt ( X**T * X) ! norm2 := sqrt ( X**T * X)
! !
! Arguments: ! Arguments:
! res - real. The result. ! res - real The result.
! x - real,dimension(:). The input vector containing the entries of X. ! x(:) - real The input vector containing the entries of X.
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code ! info - integer. Return code
! !
subroutine psb_snrm2vs(res, x, desc_a, info) subroutine psb_snrm2vs(res, x, desc_a, info)
use psb_descriptor_type use psb_descriptor_type
@ -414,15 +395,16 @@ subroutine psb_snrm2vs(res, x, desc_a, info)
use psb_penv_mod use psb_penv_mod
implicit none implicit none
real(psb_spk_), intent(in) :: x(:) real(psb_spk_), intent(in) :: x(:)
real(psb_spk_), intent(out) :: res real(psb_spk_), intent(out) :: res
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm & err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm, ldx
real(psb_spk_) :: nrm2, snrm2, dd real(psb_spk_) :: nrm2, snrm2, dd
!!$ external scombnrm2 !!$ external scombnrm2
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
@ -443,8 +425,8 @@ subroutine psb_snrm2vs(res, x, desc_a, info)
ix = 1 ix = 1
jx = 1 jx = 1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = size(x,1)
call psb_chkvect(m,1,size(x),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(m,ione,ldx,ix,jx,desc_a,info,iix,jjx)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -457,28 +439,23 @@ subroutine psb_snrm2vs(res, x, desc_a, info)
goto 9999 goto 9999
end if end if
if(m /= 0) then if (desc_a%get_local_rows() > 0) then
if (desc_a%get_local_rows() > 0) then ndim = desc_a%get_local_rows()
ndim = desc_a%get_local_rows() res = snrm2( int(ndim,kind=psb_mpik_), x, int(ione,kind=psb_mpik_) )
nrm2 = snrm2( ndim, x, ione )
! adjust because overlapped elements are computed more than once ! adjust because overlapped elements are computed more than once
do i=1,size(desc_a%ovrlap_elem,1) do i=1,size(desc_a%ovrlap_elem,1)
idx = desc_a%ovrlap_elem(i,1) idx = desc_a%ovrlap_elem(i,1)
ndm = desc_a%ovrlap_elem(i,2) ndm = desc_a%ovrlap_elem(i,2)
dd = real(ndm-1)/real(ndm) dd = real(ndm-1)/real(ndm)
nrm2 = nrm2 * sqrt(sone - dd*(abs(x(idx))/nrm2)**2) res = res * sqrt(sone - dd*(abs(x(idx))/res)**2)
end do end do
else
nrm2 = szero
end if
else else
nrm2 = szero res = szero
end if end if
!!$ call pstreecomb(ictxt,'All',1,nrm2,-1,-1,scombnrm2) call psb_nrm2(ictxt,res)
call psb_nrm2(ictxt,nrm2)
res = nrm2
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

36
base/psblas/psb_snrmi.f90
Unescape Escape View File

@ -32,34 +32,28 @@
! File: psb_snrmi.f90 ! File: psb_snrmi.f90
! !
! Function: psb_snrmi ! Function: psb_snrmi
! Forms the approximated norm of a sparse matrix, ! Forms the infinity norm of a sparse matrix,
! !
! normi := max(abs(sum(A(i,j)))) ! nrmi := max_i(abs(sum(A(i,:))))
! !
! Arguments: ! Arguments:
! a - type(psb_sspmat_type). The sparse matrix containing A. ! a - type(psb_sspmat_type). The sparse matrix containing A.
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code ! info - integer. Return code
! !
function psb_snrmi(a,desc_a,info) function psb_snrmi(a,desc_a,info) result(res)
use psb_descriptor_type use psb_base_mod, psb_protect_name => psb_snrmi
use psb_serial_mod
use psb_check_mod
use psb_error_mod
use psb_penv_mod
use psb_mat_mod
implicit none implicit none
type(psb_sspmat_type), intent(in) :: a type(psb_sspmat_type), intent(in) :: a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
real(psb_spk_) :: psb_snrmi real(psb_spk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, n, iia, jja, ia, ja, mdim, ndim, m & err_act, n, iia, jja, ia, ja, mdim, ndim, m
real(psb_spk_) :: nrmi character(len=20) :: name, ch_err
character(len=20) :: name, ch_err
name='psb_snrmi' name='psb_snrmi'
if(psb_get_errstatus() /= 0) return if(psb_get_errstatus() /= 0) return
@ -95,21 +89,19 @@ function psb_snrmi(a,desc_a,info)
end if end if
if ((m /= 0).and.(n /= 0)) then if ((m /= 0).and.(n /= 0)) then
nrmi = a%csnmi() res = a%csnmi()
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_csnmi' ch_err='psb_snrmi'
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
else else
nrmi = 0.0 res = szero
end if end if
! compute global max
call psb_amx(ictxt, nrmi)
psb_snrmi = nrmi ! compute global max
call psb_amx(ictxt, res)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

81
base/psblas/psb_sspmm.f90
Unescape Escape View File

@ -48,18 +48,18 @@
! vectors and A is a M-by-N distributed matrix. ! vectors and A is a M-by-N distributed matrix.
! !
! Arguments: ! Arguments:
! alpha - real The scalar alpha. ! alpha - real The scalar alpha.
! a - type(psb_zspmat_type). The sparse matrix containing A. ! a - type(psb_sspmat_type). The sparse matrix containing A.
! x(:,:) - real The input vector containing the entries of ( X ). ! x(:,:) - real The input vector containing the entries of ( X ).
! beta - real The scalar beta. ! beta - real The scalar beta.
! y(:,:) - real The input vector containing the entries of ( Y ). ! y(:,:) - real The input vector containing the entries of ( Y ).
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code ! info - integer. Return code
! trans - character(optional). Whether A or A'. Default: 'N' ! trans - character(optional). Whether A or A'. Default: 'N'
! k - integer(optional). The number of right-hand sides. ! k - integer(optional). The number of right-hand sides.
! jx - integer(optional). The column offset for ( X ). Default: 1 ! jx - integer(optional). The column offset for ( X ). Default: 1
! jy - integer(optional). The column offset for ( Y ). Default: 1 ! jy - integer(optional). The column offset for ( Y ). Default: 1
! work(:) - real,(optional). Working area. ! work(:) - real,(optional). Working area.
! doswap - logical(optional). Whether to performe halo updates. ! doswap - logical(optional). Whether to performe halo updates.
! !
subroutine psb_sspmm(alpha,a,x,beta,y,desc_a,info,& subroutine psb_sspmm(alpha,a,x,beta,y,desc_a,info,&
@ -84,12 +84,12 @@ subroutine psb_sspmm(alpha,a,x,beta,y,desc_a,info,&
& err_act, n, iix, jjx, ia, ja, iia, jja, ix, iy, ik, ijx, ijy,& & err_act, n, iix, jjx, ia, ja, iia, jja, ix, iy, ik, ijx, ijy,&
& m, nrow, ncol, lldx, lldy, liwork, iiy, jjy,& & m, nrow, ncol, lldx, lldy, liwork, iiy, jjy,&
& i, ib, ib1, ip, idx & i, ib, ib1, ip, idx
integer(psb_ipk_), parameter :: nb=4 integer(psb_ipk_), parameter :: nb=4
real(psb_spk_), pointer :: xp(:,:), yp(:,:), iwork(:) real(psb_spk_), pointer :: xp(:,:), yp(:,:), iwork(:)
real(psb_spk_), allocatable :: xvsave(:,:) real(psb_spk_), allocatable :: xvsave(:,:)
character :: trans_ character :: trans_
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
logical :: aliw, doswap_ logical :: aliw, doswap_
integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: debug_level, debug_unit
name='psb_sspmm' name='psb_sspmm'
@ -205,9 +205,9 @@ subroutine psb_sspmm(alpha,a,x,beta,y,desc_a,info,&
end if end if
! checking for vectors correctness ! checking for vectors correctness
call psb_chkvect(n,ik,size(x,1),ix,ijx,desc_a,info,iix,jjx) call psb_chkvect(n,ik,lldx,ix,ijx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,size(y,1),iy,ijy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ik,lldy,iy,ijy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -277,9 +277,9 @@ subroutine psb_sspmm(alpha,a,x,beta,y,desc_a,info,&
! checking for vectors correctness ! checking for vectors correctness
call psb_chkvect(m,ik,size(x,1),ix,ijx,desc_a,info,iix,jjx) call psb_chkvect(m,ik,lldx,ix,ijx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(n,ik,size(y,1),iy,ijy,desc_a,info,iiy,jjy) & call psb_chkvect(n,ik,lldy,iy,ijy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -304,7 +304,8 @@ subroutine psb_sspmm(alpha,a,x,beta,y,desc_a,info,&
if (info == psb_success_) call psi_ovrl_upd(x,desc_a,psb_avg_,info) if (info == psb_success_) call psi_ovrl_upd(x,desc_a,psb_avg_,info)
y(nrow+1:ncol,1:ik) = szero y(nrow+1:ncol,1:ik) = szero
if (info == psb_success_) call psb_csmm(alpha,a,x(:,1:ik),beta,y(:,1:ik),info,trans=trans_) if (info == psb_success_) &
& call psb_csmm(alpha,a,x(:,1:ik),beta,y(:,1:ik),info,trans=trans_)
if (debug_level >= psb_debug_comp_) & if (debug_level >= psb_debug_comp_) &
& write(debug_unit,*) me,' ',trim(name),' csmm ', info & write(debug_unit,*) me,' ',trim(name),' csmm ', info
if (info /= psb_success_) then if (info /= psb_success_) then
@ -403,15 +404,15 @@ end subroutine psb_sspmm
! vectors and A is a M-by-N distributed matrix. ! vectors and A is a M-by-N distributed matrix.
! !
! Arguments: ! Arguments:
! alpha - real The scalar alpha. ! alpha - real The scalar alpha.
! a - type(psb_zspmat_type). The sparse matrix containing A. ! a - type(psb_sspmat_type). The sparse matrix containing A.
! x(:) - real The input vector containing the entries of ( X ). ! x(:) - real The input vector containing the entries of ( X ).
! beta - real The scalar beta. ! beta - real The scalar beta.
! y(:) - real The input vector containing the entries of ( Y ). ! y(:) - real The input vector containing the entries of ( Y ).
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code ! info - integer. Return code
! trans - character(optional). Whether A or A'. Default: 'N' ! trans - character(optional). Whether A or A'. Default: 'N'
! work(:) - real,(optional). Working area. ! work(:) - real,(optional). Working area.
! doswap - logical(optional). Whether to performe halo updates. ! doswap - logical(optional). Whether to performe halo updates.
! !
subroutine psb_sspmv(alpha,a,x,beta,y,desc_a,info,& subroutine psb_sspmv(alpha,a,x,beta,y,desc_a,info,&
@ -423,10 +424,10 @@ subroutine psb_sspmv(alpha,a,x,beta,y,desc_a,info,&
real(psb_spk_), intent(in) :: alpha, beta real(psb_spk_), intent(in) :: alpha, beta
real(psb_spk_), intent(inout), target :: x(:) real(psb_spk_), intent(inout), target :: x(:)
real(psb_spk_), intent(inout), target :: y(:) real(psb_spk_), intent(inout), target :: y(:)
type(psb_sspmat_type), intent(in) :: a type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
real(psb_spk_), optional, target, intent(inout) :: work(:) real(psb_spk_), optional, target, intent(inout) :: work(:)
character, intent(in), optional :: trans character, intent(in), optional :: trans
logical, intent(in), optional :: doswap logical, intent(in), optional :: doswap
@ -435,12 +436,12 @@ subroutine psb_sspmv(alpha,a,x,beta,y,desc_a,info,&
& err_act, n, iix, jjx, ia, ja, iia, jja, ix, iy, ik, & & err_act, n, iix, jjx, ia, ja, iia, jja, ix, iy, ik, &
& m, nrow, ncol, lldx, lldy, liwork, jx, jy, iiy, jjy,& & m, nrow, ncol, lldx, lldy, liwork, jx, jy, iiy, jjy,&
& ib, ip, idx & ib, ip, idx
integer(psb_ipk_), parameter :: nb=4 integer(psb_ipk_), parameter :: nb=4
real(psb_spk_), pointer :: iwork(:), xp(:), yp(:) real(psb_spk_), pointer :: iwork(:), xp(:), yp(:)
real(psb_spk_), allocatable :: xvsave(:) real(psb_spk_), allocatable :: xvsave(:)
character :: trans_ character :: trans_
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
logical :: aliw, doswap_ logical :: aliw, doswap_
integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: debug_level, debug_unit
name='psb_sspmv' name='psb_sspmv'
@ -542,9 +543,9 @@ subroutine psb_sspmv(alpha,a,x,beta,y,desc_a,info,&
end if end if
! checking for vectors correctness ! checking for vectors correctness
call psb_chkvect(n,ik,size(x),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(n,ik,lldx,ix,jx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,size(y),iy,jy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ik,lldy,iy,jy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -582,9 +583,9 @@ subroutine psb_sspmv(alpha,a,x,beta,y,desc_a,info,&
end if end if
! checking for vectors correctness ! checking for vectors correctness
call psb_chkvect(m,ik,size(x),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(m,ik,lldx,ix,jx,desc_a,info,iix,jjx)
if (info == psb_success_)& if (info == psb_success_)&
& call psb_chkvect(n,ik,size(y),iy,jy,desc_a,info,iiy,jjy) & call psb_chkvect(n,ik,lldy,iy,jy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -674,19 +675,20 @@ subroutine psb_sspmv(alpha,a,x,beta,y,desc_a,info,&
end subroutine psb_sspmv end subroutine psb_sspmv
subroutine psb_sspmv_vect(alpha,a,x,beta,y,desc_a,info,& subroutine psb_sspmv_vect(alpha,a,x,beta,y,desc_a,info,&
& trans, work, doswap) & trans, work, doswap)
use psb_base_mod, psb_protect_name => psb_sspmv_vect use psb_base_mod, psb_protect_name => psb_sspmv_vect
use psi_mod use psi_mod
implicit none implicit none
real(psb_spk_), intent(in) :: alpha, beta real(psb_spk_), intent(in) :: alpha, beta
type(psb_s_vect_type), intent(inout) :: x type(psb_s_vect_type), intent(inout) :: x
type(psb_s_vect_type), intent(inout) :: y type(psb_s_vect_type), intent(inout) :: y
type(psb_sspmat_type), intent(in) :: a type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
real(psb_spk_), optional, target, intent(inout) :: work(:) real(psb_spk_), optional, target, intent(inout) :: work(:)
character, intent(in), optional :: trans character, intent(in), optional :: trans
logical, intent(in), optional :: doswap logical, intent(in), optional :: doswap
@ -814,9 +816,9 @@ subroutine psb_sspmv_vect(alpha,a,x,beta,y,desc_a,info,&
end if end if
! checking for vectors correctness ! checking for vectors correctness
call psb_chkvect(n,ik,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(n,ik,lldx,ix,jx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,y%get_nrows(),iy,jy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ik,lldy,iy,jy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -854,9 +856,9 @@ subroutine psb_sspmv_vect(alpha,a,x,beta,y,desc_a,info,&
end if end if
! checking for vectors correctness ! checking for vectors correctness
call psb_chkvect(m,ik,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(m,ik,lldx,ix,jx,desc_a,info,iix,jjx)
if (info == psb_success_)& if (info == psb_success_)&
& call psb_chkvect(n,ik,y%get_nrows(),iy,jy,desc_a,info,iiy,jjy) & call psb_chkvect(n,ik,lldy,iy,jy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -880,6 +882,7 @@ subroutine psb_sspmv_vect(alpha,a,x,beta,y,desc_a,info,&
! !
call psi_ovrl_save(x%v,xvsave,desc_a,info) call psi_ovrl_save(x%v,xvsave,desc_a,info)
if (info == psb_success_) call psi_ovrl_upd(x%v,desc_a,psb_avg_,info) if (info == psb_success_) call psi_ovrl_upd(x%v,desc_a,psb_avg_,info)
!!! THIS SHOULD BE FIXED !!! But beta is almost never /= 0 !!! THIS SHOULD BE FIXED !!! But beta is almost never /= 0
!!$ yp(nrow+1:ncol) = szero !!$ yp(nrow+1:ncol) = szero
@ -907,7 +910,7 @@ subroutine psb_sspmv_vect(alpha,a,x,beta,y,desc_a,info,&
& write(debug_unit,*) me,' ',trim(name),' swaptran ', info & write(debug_unit,*) me,' ',trim(name),' swaptran ', info
if(info /= psb_success_) then if(info /= psb_success_) then
info = psb_err_from_subroutine_ info = psb_err_from_subroutine_
ch_err='PSI_dSwapTran' ch_err='PSI_SwapTran'
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if

61
base/psblas/psb_sspsm.f90
Unescape Escape View File

@ -56,23 +56,22 @@
! vector and T is a M-by-M distributed triangular matrix. ! vector and T is a M-by-M distributed triangular matrix.
! !
! Arguments: ! Arguments:
! alpha - real The scalar alpha. ! alpha - real. The scalar alpha.
! a - type(psb_zspmat_type). The sparse matrix containing A. ! a - type(psb_sspmat_type). The sparse matrix containing A.
! x(:,:) - real The input vector containing the entries of ( X ). ! x(:,:) - real The input vector containing the entries of ( X ).
! beta - real The scalar beta. ! beta - real The scalar beta.
! y(:,:) - real The input vector containing the entries of ( Y ). ! y(:,:) - real The input vector containing the entries of ( Y ).
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code ! info - integer. Return code
! trans - character(optional). Whether A or A'. If not present 'N' is assumed. ! trans - character(optional). Whether A or A'. If not present 'N' is assumed.
! unitd - character(optional). Specify some type of operation with ! scale - character(optional). Specify some type of operation with
! the diagonal matrix D. ! the diagonal matrix D.
! choice - integer(optional). The kind of update to perform on overlap elements. ! choice - integer(optional). The kind of update to perform on overlap elements.
! d(:) - real , optional Matrix for diagonal scaling. ! d(:) - real, optional Matrix for diagonal scaling.
! k - integer(optional). The number of right-hand sides. ! k - integer(optional). The number of right-hand sides.
! jx - integer(optional). The column offset for ( X ). Default: 1 ! jx - integer(optional). The column offset for ( X ). Default: 1
! jy - integer(optional). The column offset for ( Y ). Default: 1 ! jy - integer(optional). The column offset for ( Y ). Default: 1
! work(:) - real , optional Working area. ! work(:) - real, optional Working area.
!
! !
subroutine psb_sspsm(alpha,a,x,beta,y,desc_a,info,& subroutine psb_sspsm(alpha,a,x,beta,y,desc_a,info,&
& trans, scale, choice, diag, k, jx, jy, work) & trans, scale, choice, diag, k, jx, jy, work)
@ -99,7 +98,7 @@ subroutine psb_sspsm(alpha,a,x,beta,y,desc_a,info,&
& m, nrow, ncol, liwork, llwork, iiy, jjy, idx, ndm & m, nrow, ncol, liwork, llwork, iiy, jjy, idx, ndm
character :: lscale character :: lscale
integer(psb_ipk_), parameter :: nb=4 integer(psb_ipk_), parameter :: nb=4
real(psb_spk_),pointer :: iwork(:), xp(:,:), yp(:,:), id(:) real(psb_spk_),pointer :: iwork(:), xp(:,:), yp(:,:), id(:)
character :: itrans character :: itrans
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
@ -221,9 +220,9 @@ subroutine psb_sspsm(alpha,a,x,beta,y,desc_a,info,&
call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja) call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja)
! checking for vectors correctness ! checking for vectors correctness
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,size(x,1),ix,ijx,desc_a,info,iix,jjx) & call psb_chkvect(m,ik,lldx,ix,ijx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,size(y,1),iy,ijy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ik,lldy,iy,ijy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect/mat' ch_err='psb_chkvect/mat'
@ -339,19 +338,19 @@ end subroutine psb_sspsm
! !
! !
! Arguments: ! Arguments:
! alpha - real The scalar alpha. ! alpha - real. The scalar alpha.
! a - type(psb_zspmat_type). The sparse matrix containing A. ! a - type(psb_sspmat_type). The sparse matrix containing A.
! x(:) - real The input vector containing the entries of ( X ). ! x(:) - real The input vector containing the entries of ( X ).
! beta - real The scalar beta. ! beta - real The scalar beta.
! y(:) - real The input vector containing the entries of ( Y ). ! y(:) - real The input vector containing the entries of ( Y ).
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code ! info - integer. Return code
! trans - character(optional). Whether A or A'. If not present 'N' is assumed. ! trans - character(optional). Whether A or A'. If not present 'N' is assumed.
! scale - character(optional). Specify some type of operation with ! scale - character(optional). Specify some type of operation with
! the diagonal matrix D. ! the diagonal matrix D.
! choice - integer(optional). The kind of update to perform on overlap elements. ! choice - integer(optional). The kind of update to perform on overlap elements.
! d(:) - real , optional Matrix for diagonal scaling. ! d(:) - real, optional Matrix for diagonal scaling.
! work(:) - real , optional Working area. ! work(:) - real, optional Working area.
! !
subroutine psb_sspsv(alpha,a,x,beta,y,desc_a,info,& subroutine psb_sspsv(alpha,a,x,beta,y,desc_a,info,&
& trans, scale, choice, diag, work) & trans, scale, choice, diag, work)
@ -471,7 +470,7 @@ subroutine psb_sspsv(alpha,a,x,beta,y,desc_a,info,&
iwork(1)=0.d0 iwork(1)=0.d0
if (present(diag)) then if(present(diag)) then
lld = size(diag) lld = size(diag)
id => diag id => diag
else else
@ -484,9 +483,9 @@ subroutine psb_sspsv(alpha,a,x,beta,y,desc_a,info,&
call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja) call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja)
! checking for vectors correctness ! checking for vectors correctness
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,size(x),ix,jx,desc_a,info,iix,jjx) & call psb_chkvect(m,ik,lldx,ix,jx,desc_a,info,iix,jjx)
if (info == psb_success_)& if (info == psb_success_) &
& call psb_chkvect(m,ik,size(y),iy,jy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ik,lldy,iy,jy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect/mat' ch_err='psb_chkvect/mat'
@ -522,7 +521,7 @@ subroutine psb_sspsv(alpha,a,x,beta,y,desc_a,info,&
end if end if
! update overlap elements ! update overlap elements
if (choice_ > 0) then if(choice_ > 0) then
call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),& call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),&
& sone,yp,desc_a,iwork,info,data=psb_comm_ovr_) & sone,yp,desc_a,iwork,info,data=psb_comm_ovr_)
@ -550,22 +549,23 @@ subroutine psb_sspsv(alpha,a,x,beta,y,desc_a,info,&
return return
end subroutine psb_sspsv end subroutine psb_sspsv
subroutine psb_sspsv_vect(alpha,a,x,beta,y,desc_a,info,& subroutine psb_sspsv_vect(alpha,a,x,beta,y,desc_a,info,&
& trans, scale, choice, diag, work) & trans, scale, choice, diag, work)
use psb_base_mod, psb_protect_name => psb_sspsv_vect use psb_base_mod, psb_protect_name => psb_sspsv_vect
use psi_mod use psi_mod
implicit none implicit none
real(psb_spk_), intent(in) :: alpha, beta real(psb_spk_), intent(in) :: alpha, beta
type(psb_s_vect_type), intent(inout) :: x type(psb_s_vect_type), intent(inout) :: x
type(psb_s_vect_type), intent(inout) :: y type(psb_s_vect_type), intent(inout) :: y
type(psb_sspmat_type), intent(inout) :: a type(psb_sspmat_type), intent(inout) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
type(psb_s_vect_type), intent(inout), optional :: diag type(psb_s_vect_type), intent(inout), optional :: diag
real(psb_spk_), optional, target, intent(inout) :: work(:) real(psb_spk_), optional, target, intent(inout) :: work(:)
character, intent(in), optional :: trans, scale character, intent(in), optional :: trans, scale
integer(psb_ipk_), intent(in), optional :: choice integer(psb_ipk_), intent(in), optional :: choice
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me, & integer(psb_ipk_) :: ictxt, np, me, &
@ -684,9 +684,9 @@ subroutine psb_sspsv_vect(alpha,a,x,beta,y,desc_a,info,&
call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja) call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja)
! checking for vectors correctness ! checking for vectors correctness
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) & call psb_chkvect(m,ik,lldx,ix,jx,desc_a,info,iix,jjx)
if (info == psb_success_)& if (info == psb_success_)&
& call psb_chkvect(m,ik,y%get_nrows(),iy,jy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ik,lldy,iy,jy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect/mat' ch_err='psb_chkvect/mat'
@ -750,4 +750,3 @@ subroutine psb_sspsv_vect(alpha,a,x,beta,y,desc_a,info,&
return return
end subroutine psb_sspsv_vect end subroutine psb_sspsv_vect

41
base/psblas/psb_zaxpby.f90
Unescape Escape View File

@ -129,8 +129,8 @@ end subroutine psb_zaxpby_vect
! Arguments: ! Arguments:
! alpha - complex,input The scalar used to multiply each component of X ! alpha - complex,input The scalar used to multiply each component of X
! x(:,:) - complex,input The input vector containing the entries of X ! x(:,:) - complex,input The input vector containing the entries of X
! beta - real,input The scalar used to multiply each component of Y ! beta - complex,input The scalar used to multiply each component of Y
! y(:,:) - real,inout The input vector Y ! y(:,:) - complex,inout The input vector Y
! desc_a - type(psb_desc_type) The communication descriptor. ! desc_a - type(psb_desc_type) The communication descriptor.
! info - integer Return code ! info - integer Return code
! jx - integer(optional) The column offset for X ! jx - integer(optional) The column offset for X
@ -138,6 +138,7 @@ end subroutine psb_zaxpby_vect
! !
subroutine psb_zaxpby(alpha, x, beta,y,desc_a,info, n, jx, jy) subroutine psb_zaxpby(alpha, x, beta,y,desc_a,info, n, jx, jy)
use psb_base_mod, psb_protect_name => psb_zaxpby use psb_base_mod, psb_protect_name => psb_zaxpby
implicit none implicit none
integer(psb_ipk_), intent(in), optional :: n, jx, jy integer(psb_ipk_), intent(in), optional :: n, jx, jy
@ -149,7 +150,8 @@ subroutine psb_zaxpby(alpha, x, beta,y,desc_a,info, n, jx, jy)
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ix, iy, ijx, ijy, m, iiy, in, jjy & err_act, iix, jjx, ix, iy, ijx, ijy, m, iiy, in, jjy, &
& lldx, lldy
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_geaxpby' name='psb_geaxpby'
@ -197,11 +199,12 @@ subroutine psb_zaxpby(alpha, x, beta,y,desc_a,info, n, jx, jy)
end if end if
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
lldx = size(x,1)
lldy = size(y,1)
! check vector correctness ! check vector correctness
call psb_chkvect(m,ione,size(x,1),ix,ijx,desc_a,info,iix,jjx) call psb_chkvect(m,ione,lldx,ix,ijx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ione,size(y,1),iy,ijy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ione,lldy,iy,ijy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -217,9 +220,9 @@ subroutine psb_zaxpby(alpha, x, beta,y,desc_a,info, n, jx, jy)
if ((in /= 0)) then if ((in /= 0)) then
if(desc_a%get_local_rows() > 0) then if(desc_a%get_local_rows() > 0) then
call zaxpby(desc_a%get_local_cols(),in,& call caxpby(desc_a%get_local_cols(),in,&
& alpha,x(iix:,jjx),size(x,1),beta,& & alpha,x(iix:,jjx),lldx,beta,&
& y(iiy:,jjy),size(y,1),info) & y(iiy:,jjy),lldy,info)
end if end if
end if end if
@ -280,8 +283,8 @@ end subroutine psb_zaxpby
! Arguments: ! Arguments:
! alpha - complex,input The scalar used to multiply each component of X ! alpha - complex,input The scalar used to multiply each component of X
! x(:) - complex,input The input vector containing the entries of X ! x(:) - complex,input The input vector containing the entries of X
! beta - real,input The scalar used to multiply each component of Y ! beta - complex,input The scalar used to multiply each component of Y
! y(:) - real,inout The input vector Y ! y(:) - complex,inout The input vector Y
! desc_a - type(psb_desc_type) The communication descriptor. ! desc_a - type(psb_desc_type) The communication descriptor.
! info - integer Return code ! info - integer Return code
! !
@ -298,7 +301,8 @@ subroutine psb_zaxpbyv(alpha, x, beta,y,desc_a,info)
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ix, iy, m, iiy, jjy & err_act, iix, jjx, ix, iy, m, iiy, jjy, &
& lldx, lldy
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
@ -320,16 +324,17 @@ subroutine psb_zaxpbyv(alpha, x, beta,y,desc_a,info)
iy = ione iy = ione
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
lldx = size(x,1)
lldy = size(y,1)
! check vector correctness ! check vector correctness
call psb_chkvect(m,ione,size(x),ix,ione,desc_a,info,iix,jjx) call psb_chkvect(m,ione,lldx,ix,ione,desc_a,info,iix,jjx)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect 1' ch_err='psb_chkvect 1'
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
call psb_chkvect(m,ione,size(y),iy,ione,desc_a,info,iiy,jjy) call psb_chkvect(m,ione,lldy,iy,ione,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect 2' ch_err='psb_chkvect 2'
@ -343,9 +348,9 @@ subroutine psb_zaxpbyv(alpha, x, beta,y,desc_a,info)
end if end if
if(desc_a%get_local_rows() > 0) then if(desc_a%get_local_rows() > 0) then
call zaxpby(desc_a%get_local_cols(),ione,& call caxpby(desc_a%get_local_cols(),ione,&
& alpha,x,size(x),beta,& & alpha,x,lldx,beta,&
& y,size(y),info) & y,lldy,info)
end if end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)

174
base/psblas/psb_znrm2.f90
Unescape Escape View File

@ -44,7 +44,7 @@
! info - integer. Return code ! info - integer. Return code
! jx - integer(optional). The column offset for sub( X ). ! jx - integer(optional). The column offset for sub( X ).
! !
function psb_znrm2(x, desc_a, info, jx) function psb_znrm2(x, desc_a, info, jx) result(res)
use psb_descriptor_type use psb_descriptor_type
use psb_check_mod use psb_check_mod
use psb_error_mod use psb_error_mod
@ -55,15 +55,13 @@ function psb_znrm2(x, desc_a, info, jx)
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(in), optional :: jx integer(psb_ipk_), intent(in), optional :: jx
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
real(psb_dpk_) :: psb_znrm2 real(psb_dpk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ndim, ix, ijx, i, m, id, idx, ndm & err_act, iix, jjx, ndim, ix, ijx, i, m, id, idx, ndm, ldx
real(psb_dpk_) :: nrm2, dznrm2, dd real(psb_dpk_) :: dznrm2, dd
character(len=20) :: name, ch_err
!!$ external dcombnrm2
character(len=20) :: name, ch_err
name='psb_znrm2' name='psb_znrm2'
if(psb_get_errstatus() /= 0) return if(psb_get_errstatus() /= 0) return
@ -87,7 +85,8 @@ function psb_znrm2(x, desc_a, info, jx)
endif endif
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
call psb_chkvect(m,1,size(x,1),ix,ijx,desc_a,info,iix,jjx) ldx = size(x,1)
call psb_chkvect(m,ione,ldx,ix,ijx,desc_a,info,iix,jjx)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -100,29 +99,22 @@ function psb_znrm2(x, desc_a, info, jx)
goto 9999 goto 9999
end if end if
if(m /= 0) then if (desc_a%get_local_rows() > 0) then
if (desc_a%get_local_rows() > 0) then ndim = desc_a%get_local_rows()
ndim = desc_a%get_local_rows() res = dznrm2( int(ndim,kind=psb_mpik_), x(iix:,jjx), int(ione,kind=psb_mpik_) )
nrm2 = dznrm2( ndim, x(iix:,jjx), ione )
! adjust because overlapped elements are computed more than once
! adjust because overlapped elements are computed more than once do i=1,size(desc_a%ovrlap_elem,1)
do i=1,size(desc_a%ovrlap_elem,1) idx = desc_a%ovrlap_elem(i,1)
idx = desc_a%ovrlap_elem(i,1) ndm = desc_a%ovrlap_elem(i,2)
ndm = desc_a%ovrlap_elem(i,2) dd = real(ndm-1)/real(ndm)
dd = dble(ndm-1)/dble(ndm) res = res * sqrt(sone - dd*(abs(x(idx,jjx))/res)**2)
nrm2 = nrm2 * sqrt(done - dd*(abs(x(idx,jjx))/nrm2)**2) end do
end do
else
nrm2 = dzero
end if
else else
nrm2 = dzero res = dzero
end if end if
!!$ call pdtreecomb(ictxt,'All',1,nrm2,-1,-1,dcombnrm2) call psb_nrm2(ictxt,res)
call psb_nrm2(ictxt,nrm2)
psb_znrm2 = nrm2
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -181,24 +173,24 @@ end function psb_znrm2
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code ! info - integer. Return code
! !
function psb_znrm2v(x, desc_a, info) function psb_znrm2v(x, desc_a, info) result(res)
use psb_descriptor_type use psb_descriptor_type
use psb_check_mod use psb_check_mod
use psb_error_mod use psb_error_mod
use psb_penv_mod use psb_penv_mod
implicit none implicit none
complex(psb_dpk_), intent(in) :: x(:) complex(psb_dpk_), intent(in) :: x(:)
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
real(psb_dpk_) :: psb_znrm2v real(psb_dpk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm & err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm, ldx
real(psb_dpk_) :: nrm2, dznrm2, dd real(psb_dpk_) :: dznrm2, dd
!!$ external dcombnrm2 !!$ external scombnrm2
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_znrm2v' name='psb_znrm2v'
@ -218,8 +210,8 @@ function psb_znrm2v(x, desc_a, info)
ix = 1 ix = 1
jx=1 jx=1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = size(x,1)
call psb_chkvect(m,1,size(x),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(m,ione,ldx,ix,jx,desc_a,info,iix,jjx)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -232,28 +224,22 @@ function psb_znrm2v(x, desc_a, info)
goto 9999 goto 9999
end if end if
if(m /= 0) then if (desc_a%get_local_rows() > 0) then
if (desc_a%get_local_rows() > 0) then ndim = desc_a%get_local_rows()
ndim = desc_a%get_local_rows() res = dznrm2( int(ndim,kind=psb_mpik_), x, int(ione,kind=psb_mpik_) )
nrm2 = dznrm2( ndim, x, ione ) ! adjust because overlapped elements are computed more than once
! adjust because overlapped elements are computed more than once do i=1,size(desc_a%ovrlap_elem,1)
do i=1,size(desc_a%ovrlap_elem,1) idx = desc_a%ovrlap_elem(i,1)
idx = desc_a%ovrlap_elem(i,1) ndm = desc_a%ovrlap_elem(i,2)
ndm = desc_a%ovrlap_elem(i,2) dd = real(ndm-1)/real(ndm)
dd = dble(ndm-1)/dble(ndm) res = res * sqrt(sone - dd*(abs(x(idx))/res)**2)
nrm2 = nrm2 * sqrt(done - dd*(abs(x(idx))/nrm2)**2) end do
end do
else
nrm2 = dzero
end if
else else
nrm2 = dzero res = dzero
end if end if
!!$ call pdtreecomb(ictxt,'All',1,nrm2,-1,-1,dcombnrm2) call psb_nrm2(ictxt,res)
call psb_nrm2(ictxt,nrm2)
psb_znrm2v = nrm2
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -269,6 +255,7 @@ function psb_znrm2v(x, desc_a, info)
end function psb_znrm2v end function psb_znrm2v
function psb_znrm2_vect(x, desc_a, info) result(res) function psb_znrm2_vect(x, desc_a, info) result(res)
use psb_descriptor_type use psb_descriptor_type
use psb_check_mod use psb_check_mod
@ -280,12 +267,13 @@ function psb_znrm2_vect(x, desc_a, info) result(res)
real(psb_dpk_) :: res real(psb_dpk_) :: res
type(psb_z_vect_type), intent (inout) :: x type(psb_z_vect_type), intent (inout) :: x
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm & err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm, ldx
real(psb_dpk_) :: nrm2 real(psb_dpk_) :: snrm2, dd
!!$ external dcombnrm2
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_znrm2v' name='psb_znrm2v'
@ -312,8 +300,8 @@ function psb_znrm2_vect(x, desc_a, info) result(res)
ix = 1 ix = 1
jx=1 jx=1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = x%get_nrows()
call psb_chkvect(m,1,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(m,ione,ldx,ix,jx,desc_a,info,iix,jjx)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -326,10 +314,9 @@ function psb_znrm2_vect(x, desc_a, info) result(res)
goto 9999 goto 9999
end if end if
if(m /= 0) then if (desc_a%get_local_rows() > 0) then
if (desc_a%get_local_rows() > 0) then ndim = desc_a%get_local_rows()
ndim = desc_a%get_local_rows() res = x%nrm2(ndim)
nrm2 = x%nrm2(ndim)
!!$ ! adjust because overlapped elements are computed more than once !!$ ! adjust because overlapped elements are computed more than once
!!$ do i=1,size(desc_a%ovrlap_elem,1) !!$ do i=1,size(desc_a%ovrlap_elem,1)
!!$ idx = desc_a%ovrlap_elem(i,1) !!$ idx = desc_a%ovrlap_elem(i,1)
@ -337,16 +324,11 @@ function psb_znrm2_vect(x, desc_a, info) result(res)
!!$ dd = dble(ndm-1)/dble(ndm) !!$ dd = dble(ndm-1)/dble(ndm)
!!$ nrm2 = nrm2 * sqrt(done - dd*(abs(x(idx))/nrm2)**2) !!$ nrm2 = nrm2 * sqrt(done - dd*(abs(x(idx))/nrm2)**2)
!!$ end do !!$ end do
else
nrm2 = dzero
end if
else else
nrm2 = dzero res = dzero
end if end if
call psb_nrm2(ictxt,nrm2) call psb_nrm2(ictxt,res)
res = nrm2
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -363,7 +345,6 @@ end function psb_znrm2_vect
!!$ !!$
!!$ Parallel Sparse BLAS version 3.0 !!$ Parallel Sparse BLAS version 3.0
!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010 !!$ (C) Copyright 2006, 2007, 2008, 2009, 2010
@ -414,17 +395,17 @@ subroutine psb_znrm2vs(res, x, desc_a, info)
use psb_penv_mod use psb_penv_mod
implicit none implicit none
complex(psb_dpk_), intent(in) :: x(:) complex(psb_dpk_), intent(in) :: x(:)
real(psb_dpk_), intent(out) :: res real(psb_dpk_), intent(out) :: res
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm & err_act, iix, jjx, ndim, ix, jx, i, m, id, idx, ndm, ldx
real(psb_dpk_) :: nrm2, dznrm2, dd real(psb_dpk_) :: nrm2, dznrm2, dd
!!$ external dcombnrm2 !!$ external scombnrm2
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_znrm2' name='psb_znrm2'
@ -444,8 +425,8 @@ subroutine psb_znrm2vs(res, x, desc_a, info)
ix = 1 ix = 1
jx = 1 jx = 1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = size(x,1)
call psb_chkvect(m,1,size(x),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(m,ione,ldx,ix,jx,desc_a,info,iix,jjx)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -458,28 +439,23 @@ subroutine psb_znrm2vs(res, x, desc_a, info)
goto 9999 goto 9999
end if end if
if(m /= 0) then if (desc_a%get_local_rows() > 0) then
if (desc_a%get_local_rows() > 0) then ndim = desc_a%get_local_rows()
ndim = desc_a%get_local_rows() res = dznrm2( int(ndim,kind=psb_mpik_), x, int(ione,kind=psb_mpik_) )
nrm2 = dznrm2( ndim, x, ione )
! adjust because overlapped elements are computed more than once ! adjust because overlapped elements are computed more than once
do i=1,size(desc_a%ovrlap_elem,1) do i=1,size(desc_a%ovrlap_elem,1)
idx = desc_a%ovrlap_elem(i,1) idx = desc_a%ovrlap_elem(i,1)
ndm = desc_a%ovrlap_elem(i,2) ndm = desc_a%ovrlap_elem(i,2)
dd = dble(ndm-1)/dble(ndm) dd = real(ndm-1)/real(ndm)
nrm2 = nrm2 * sqrt(done - dd*(abs(x(idx))/nrm2)**2) res = res * sqrt(sone - dd*(abs(x(idx))/res)**2)
end do end do
else
nrm2 = dzero
end if
else else
nrm2 = dzero res = dzero
end if end if
!!$ call pdtreecomb(ictxt,'All',1,nrm2,-1,-1,dcombnrm2) call psb_nrm2(ictxt,res)
call psb_nrm2(ictxt,nrm2)
res = nrm2
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

35
base/psblas/psb_znrmi.f90
Unescape Escape View File

@ -32,33 +32,27 @@
! File: psb_znrmi.f90 ! File: psb_znrmi.f90
! !
! Function: psb_znrmi ! Function: psb_znrmi
! Forms the approximated norm of a sparse matrix, ! Forms the infinity norm of a sparse matrix,
! !
! normi := max(abs(sum(A(i,j)))) ! nrmi := max_i(abs(sum(A(i,:))))
! !
! Arguments: ! Arguments:
! a - type(psb_dspmat_type). The sparse matrix containing A. ! a - type(psb_zspmat_type). The sparse matrix containing A.
! desc_a - type(psb_desc_type). The communication descriptor. ! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code ! info - integer. Return code
! !
function psb_znrmi(a,desc_a,info) function psb_znrmi(a,desc_a,info) result(res)
use psb_descriptor_type use psb_base_mod, psb_protect_name => psb_znrmi
use psb_serial_mod
use psb_check_mod
use psb_error_mod
use psb_penv_mod
use psb_mat_mod
implicit none implicit none
type(psb_zspmat_type), intent(in) :: a type(psb_zspmat_type), intent(in) :: a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
real(psb_dpk_) :: psb_znrmi real(psb_dpk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, n, iia, jja, ia, ja, mdim, ndim, m & err_act, n, iia, jja, ia, ja, mdim, ndim, m
real(psb_dpk_) :: nrmi
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_znrmi' name='psb_znrmi'
@ -95,20 +89,19 @@ function psb_znrmi(a,desc_a,info)
end if end if
if ((m /= 0).and.(n /= 0)) then if ((m /= 0).and.(n /= 0)) then
nrmi = a%csnmi() res = a%csnmi()
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_csnmi' ch_err='psb_znrmi'
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
else else
nrmi = 0.d0 res = dzero
end if end if
! compute global max
call psb_amx(ictxt, nrmi)
psb_znrmi = nrmi ! compute global max
call psb_amx(ictxt, res)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

31
base/psblas/psb_zspmm.f90
Unescape Escape View File

@ -205,9 +205,9 @@ subroutine psb_zspmm(alpha,a,x,beta,y,desc_a,info,&
end if end if
! checking for vectors correctness ! checking for vectors correctness
call psb_chkvect(n,ik,size(x,1),ix,ijx,desc_a,info,iix,jjx) call psb_chkvect(n,ik,lldx,ix,ijx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,size(y,1),iy,ijy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ik,lldy,iy,ijy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -277,9 +277,9 @@ subroutine psb_zspmm(alpha,a,x,beta,y,desc_a,info,&
! checking for vectors correctness ! checking for vectors correctness
call psb_chkvect(m,ik,size(x,1),ix,ijx,desc_a,info,iix,jjx) call psb_chkvect(m,ik,lldx,ix,ijx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(n,ik,size(y,1),iy,ijy,desc_a,info,iiy,jjy) & call psb_chkvect(n,ik,lldy,iy,ijy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -304,7 +304,8 @@ subroutine psb_zspmm(alpha,a,x,beta,y,desc_a,info,&
if (info == psb_success_) call psi_ovrl_upd(x,desc_a,psb_avg_,info) if (info == psb_success_) call psi_ovrl_upd(x,desc_a,psb_avg_,info)
y(nrow+1:ncol,1:ik) = zzero y(nrow+1:ncol,1:ik) = zzero
if (info == psb_success_) call psb_csmm(alpha,a,x(:,1:ik),beta,y(:,1:ik),info,trans=trans_) if (info == psb_success_) &
& call psb_csmm(alpha,a,x(:,1:ik),beta,y(:,1:ik),info,trans=trans_)
if (debug_level >= psb_debug_comp_) & if (debug_level >= psb_debug_comp_) &
& write(debug_unit,*) me,' ',trim(name),' csmm ', info & write(debug_unit,*) me,' ',trim(name),' csmm ', info
if (info /= psb_success_) then if (info /= psb_success_) then
@ -542,9 +543,9 @@ subroutine psb_zspmv(alpha,a,x,beta,y,desc_a,info,&
end if end if
! checking for vectors correctness ! checking for vectors correctness
call psb_chkvect(n,ik,size(x),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(n,ik,lldx,ix,jx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,size(y),iy,jy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ik,lldy,iy,jy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -582,9 +583,9 @@ subroutine psb_zspmv(alpha,a,x,beta,y,desc_a,info,&
end if end if
! checking for vectors correctness ! checking for vectors correctness
call psb_chkvect(m,ik,size(x),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(m,ik,lldx,ix,jx,desc_a,info,iix,jjx)
if (info == psb_success_)& if (info == psb_success_)&
& call psb_chkvect(n,ik,size(y),iy,jy,desc_a,info,iiy,jjy) & call psb_chkvect(n,ik,lldy,iy,jy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -674,6 +675,7 @@ subroutine psb_zspmv(alpha,a,x,beta,y,desc_a,info,&
end subroutine psb_zspmv end subroutine psb_zspmv
subroutine psb_zspmv_vect(alpha,a,x,beta,y,desc_a,info,& subroutine psb_zspmv_vect(alpha,a,x,beta,y,desc_a,info,&
& trans, work, doswap) & trans, work, doswap)
use psb_base_mod, psb_protect_name => psb_zspmv_vect use psb_base_mod, psb_protect_name => psb_zspmv_vect
@ -814,9 +816,9 @@ subroutine psb_zspmv_vect(alpha,a,x,beta,y,desc_a,info,&
end if end if
! checking for vectors correctness ! checking for vectors correctness
call psb_chkvect(n,ik,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(n,ik,lldx,ix,jx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,y%get_nrows(),iy,jy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ik,lldy,iy,jy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -854,9 +856,9 @@ subroutine psb_zspmv_vect(alpha,a,x,beta,y,desc_a,info,&
end if end if
! checking for vectors correctness ! checking for vectors correctness
call psb_chkvect(m,ik,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) call psb_chkvect(m,ik,lldx,ix,jx,desc_a,info,iix,jjx)
if (info == psb_success_)& if (info == psb_success_)&
& call psb_chkvect(n,ik,y%get_nrows(),iy,jy,desc_a,info,iiy,jjy) & call psb_chkvect(n,ik,lldy,iy,jy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect' ch_err='psb_chkvect'
@ -880,8 +882,9 @@ subroutine psb_zspmv_vect(alpha,a,x,beta,y,desc_a,info,&
! !
call psi_ovrl_save(x%v,xvsave,desc_a,info) call psi_ovrl_save(x%v,xvsave,desc_a,info)
if (info == psb_success_) call psi_ovrl_upd(x%v,desc_a,psb_avg_,info) if (info == psb_success_) call psi_ovrl_upd(x%v,desc_a,psb_avg_,info)
!!! THIS SHOULD BE FIXED !!! But beta is almost never /= 0 !!! THIS SHOULD BE FIXED !!! But beta is almost never /= 0
!!$ yp(nrow+1:ncol) = szero !!$ yp(nrow+1:ncol) = zzero
! local Matrix-vector product ! local Matrix-vector product
if (info == psb_success_) call psb_csmm(alpha,a,x,beta,y,info,trans=trans_) if (info == psb_success_) call psb_csmm(alpha,a,x,beta,y,info,trans=trans_)

14
base/psblas/psb_zspsm.f90
Unescape Escape View File

@ -98,7 +98,7 @@ subroutine psb_zspsm(alpha,a,x,beta,y,desc_a,info,&
& m, nrow, ncol, liwork, llwork, iiy, jjy, idx, ndm & m, nrow, ncol, liwork, llwork, iiy, jjy, idx, ndm
character :: lscale character :: lscale
integer(psb_ipk_), parameter :: nb=4 integer(psb_ipk_), parameter :: nb=4
complex(psb_dpk_),pointer :: iwork(:), xp(:,:), yp(:,:), id(:) complex(psb_dpk_),pointer :: iwork(:), xp(:,:), yp(:,:), id(:)
character :: itrans character :: itrans
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
@ -220,9 +220,9 @@ subroutine psb_zspsm(alpha,a,x,beta,y,desc_a,info,&
call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja) call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja)
! checking for vectors correctness ! checking for vectors correctness
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,size(x,1),ix,ijx,desc_a,info,iix,jjx) & call psb_chkvect(m,ik,lldx,ix,ijx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,size(y,1),iy,ijy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ik,lldy,iy,ijy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect/mat' ch_err='psb_chkvect/mat'
@ -483,9 +483,9 @@ subroutine psb_zspsv(alpha,a,x,beta,y,desc_a,info,&
call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja) call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja)
! checking for vectors correctness ! checking for vectors correctness
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,size(x),ix,jx,desc_a,info,iix,jjx) & call psb_chkvect(m,ik,lldx,ix,jx,desc_a,info,iix,jjx)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,size(y),iy,jy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ik,lldy,iy,jy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect/mat' ch_err='psb_chkvect/mat'
@ -684,9 +684,9 @@ subroutine psb_zspsv_vect(alpha,a,x,beta,y,desc_a,info,&
call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja) call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja)
! checking for vectors correctness ! checking for vectors correctness
if (info == psb_success_) & if (info == psb_success_) &
& call psb_chkvect(m,ik,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) & call psb_chkvect(m,ik,lldx,ix,jx,desc_a,info,iix,jjx)
if (info == psb_success_)& if (info == psb_success_)&
& call psb_chkvect(m,ik,y%get_nrows(),iy,jy,desc_a,info,iiy,jjy) & call psb_chkvect(m,ik,lldy,iy,jy,desc_a,info,iiy,jjy)
if(info /= psb_success_) then if(info /= psb_success_) then
info=psb_err_from_subroutine_ info=psb_err_from_subroutine_
ch_err='psb_chkvect/mat' ch_err='psb_chkvect/mat'

Write Preview
Loading…
Cancel
Save