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

psblas3:

Browse Source 
base/modules/psb_c_base_mat_mod.f90
 base/modules/psb_c_base_vect_mod.f90
 base/modules/psb_c_csc_mat_mod.f90
 base/modules/psb_c_csr_mat_mod.f90
 base/modules/psb_c_linmap_mod.f90
 base/modules/psb_c_vect_mod.f90
 base/modules/psb_const_mod.F90
 base/modules/psb_d_base_mat_mod.f90
 base/modules/psb_d_base_vect_mod.f90
 base/modules/psb_d_csc_mat_mod.f90
 base/modules/psb_d_csr_mat_mod.f90
 base/modules/psb_d_linmap_mod.f90
 base/modules/psb_d_vect_mod.f90
 base/modules/psb_s_base_mat_mod.f90
 base/modules/psb_s_base_vect_mod.f90
 base/modules/psb_s_csc_mat_mod.f90
 base/modules/psb_s_csr_mat_mod.f90
 base/modules/psb_s_linmap_mod.f90
 base/modules/psb_s_vect_mod.f90
 base/modules/psb_serial_mod.f90
 base/modules/psb_z_base_mat_mod.f90
 base/modules/psb_z_base_vect_mod.f90
 base/modules/psb_z_csc_mat_mod.f90
 base/modules/psb_z_csr_mat_mod.f90
 base/modules/psb_z_linmap_mod.f90
 base/modules/psb_z_vect_mod.f90
 base/psblas/psb_camax.f90
 base/psblas/psb_casum.f90
 base/psblas/psb_damax.f90
 base/psblas/psb_dasum.f90
 base/psblas/psb_samax.f90
 base/psblas/psb_sasum.f90
 base/psblas/psb_zamax.f90
 base/psblas/psb_zasum.f90
 base/serial/Makefile

Work on long integers: set up a restart point.
psblas3-type-indexed
Salvatore Filippone 13 years ago
parent 235f1de194
commit 62ae695d5f
35 changed files with 1002 additions and 1262 deletions
Show Stats Download Patch File Download Diff File

6
base/modules/psb_c_base_mat_mod.f90
Unescape Escape View File

@ -1170,9 +1170,9 @@ contains
if (allocated(a%ja)) deallocate(a%ja) if (allocated(a%ja)) deallocate(a%ja)
if (allocated(a%val)) deallocate(a%val) if (allocated(a%val)) deallocate(a%val)
call a%set_null() call a%set_null()
call a%set_nrows(0) call a%set_nrows(izero)
call a%set_ncols(0) call a%set_ncols(izero)
call a%set_nzeros(0) call a%set_nzeros(izero)
return return

2
base/modules/psb_c_base_vect_mod.f90
Unescape Escape View File

@ -150,7 +150,7 @@ contains
integer(psb_ipk_) :: info integer(psb_ipk_) :: info
this%v = x this%v = x
call this%asb(size(x),info) call this%asb(size(x,kind=psb_ipk_),info)
end function constructor end function constructor

4
base/modules/psb_c_csc_mat_mod.f90
Unescape Escape View File

@ -527,8 +527,8 @@ contains
if (allocated(a%ia)) deallocate(a%ia) if (allocated(a%ia)) deallocate(a%ia)
if (allocated(a%val)) deallocate(a%val) if (allocated(a%val)) deallocate(a%val)
call a%set_null() call a%set_null()
call a%set_nrows(0) call a%set_nrows(izero)
call a%set_ncols(0) call a%set_ncols(izero)
return return

4
base/modules/psb_c_csr_mat_mod.f90
Unescape Escape View File

@ -527,8 +527,8 @@ contains
if (allocated(a%ja)) deallocate(a%ja) if (allocated(a%ja)) deallocate(a%ja)
if (allocated(a%val)) deallocate(a%val) if (allocated(a%val)) deallocate(a%val)
call a%set_null() call a%set_null()
call a%set_nrows(0) call a%set_nrows(izero)
call a%set_ncols(0) call a%set_ncols(izero)
return return

8
base/modules/psb_c_linmap_mod.f90
Unescape Escape View File

@ -55,7 +55,7 @@ module psb_c_linmap_mod
interface psb_map_X2Y interface psb_map_X2Y
subroutine psb_c_map_X2Y(alpha,x,beta,y,map,info,work) subroutine psb_c_map_X2Y(alpha,x,beta,y,map,info,work)
use psb_const_mod use psb_const_mod
import :: psb_ipk_, psb_clinmap_type import :: psb_clinmap_type
implicit none implicit none
type(psb_clinmap_type), intent(in) :: map type(psb_clinmap_type), intent(in) :: map
complex(psb_spk_), intent(in) :: alpha,beta complex(psb_spk_), intent(in) :: alpha,beta
@ -67,7 +67,7 @@ module psb_c_linmap_mod
subroutine psb_c_map_X2Y_vect(alpha,x,beta,y,map,info,work) subroutine psb_c_map_X2Y_vect(alpha,x,beta,y,map,info,work)
use psb_const_mod use psb_const_mod
use psb_c_vect_mod use psb_c_vect_mod
import :: psb_ipk_, psb_clinmap_type import :: psb_clinmap_type
implicit none implicit none
type(psb_clinmap_type), intent(in) :: map type(psb_clinmap_type), intent(in) :: map
complex(psb_spk_), intent(in) :: alpha,beta complex(psb_spk_), intent(in) :: alpha,beta
@ -80,7 +80,7 @@ module psb_c_linmap_mod
interface psb_map_Y2X interface psb_map_Y2X
subroutine psb_c_map_Y2X(alpha,x,beta,y,map,info,work) subroutine psb_c_map_Y2X(alpha,x,beta,y,map,info,work)
use psb_const_mod use psb_const_mod
import :: psb_ipk_, psb_clinmap_type import :: psb_clinmap_type
implicit none implicit none
type(psb_clinmap_type), intent(in) :: map type(psb_clinmap_type), intent(in) :: map
complex(psb_spk_), intent(in) :: alpha,beta complex(psb_spk_), intent(in) :: alpha,beta
@ -92,7 +92,7 @@ module psb_c_linmap_mod
subroutine psb_c_map_Y2X_vect(alpha,x,beta,y,map,info,work) subroutine psb_c_map_Y2X_vect(alpha,x,beta,y,map,info,work)
use psb_const_mod use psb_const_mod
use psb_c_vect_mod use psb_c_vect_mod
import :: psb_ipk_, psb_clinmap_type import :: psb_clinmap_type
implicit none implicit none
type(psb_clinmap_type), intent(in) :: map type(psb_clinmap_type), intent(in) :: map
complex(psb_spk_), intent(in) :: alpha,beta complex(psb_spk_), intent(in) :: alpha,beta

2
base/modules/psb_c_vect_mod.f90
Unescape Escape View File

@ -163,7 +163,7 @@ contains
if (info == 0) call this%v%bld(x) if (info == 0) call this%v%bld(x)
call this%asb(size(x),info) call this%asb(size(x,kind=psb_ipk_),info)
end function constructor end function constructor

3
base/modules/psb_const_mod.F90
Unescape Escape View File

@ -44,6 +44,9 @@ module psb_const_mod
! This is an 8-byte integer, and normally different from default integer. ! This is an 8-byte integer, and normally different from default integer.
integer, parameter :: longndig=12 integer, parameter :: longndig=12
integer, parameter :: psb_long_int_k_ = selected_int_kind(longndig) integer, parameter :: psb_long_int_k_ = selected_int_kind(longndig)
! This is always a 4-byte integer, for MPI-related stuff
integer, parameter :: mpindig=8
integer, parameter :: psb_mpik_ = selected_int_kind(mpindig)
! !
! These must be the kind parameter corresponding to MPI_DOUBLE_PRECISION ! These must be the kind parameter corresponding to MPI_DOUBLE_PRECISION
! and MPI_REAL ! and MPI_REAL

6
base/modules/psb_d_base_mat_mod.f90
Unescape Escape View File

@ -1170,9 +1170,9 @@ contains
if (allocated(a%ja)) deallocate(a%ja) if (allocated(a%ja)) deallocate(a%ja)
if (allocated(a%val)) deallocate(a%val) if (allocated(a%val)) deallocate(a%val)
call a%set_null() call a%set_null()
call a%set_nrows(0) call a%set_nrows(izero)
call a%set_ncols(0) call a%set_ncols(izero)
call a%set_nzeros(0) call a%set_nzeros(izero)
return return

2
base/modules/psb_d_base_vect_mod.f90
Unescape Escape View File

@ -150,7 +150,7 @@ contains
integer(psb_ipk_) :: info integer(psb_ipk_) :: info
this%v = x this%v = x
call this%asb(size(x),info) call this%asb(size(x,kind=psb_ipk_),info)
end function constructor end function constructor

4
base/modules/psb_d_csc_mat_mod.f90
Unescape Escape View File

@ -527,8 +527,8 @@ contains
if (allocated(a%ia)) deallocate(a%ia) if (allocated(a%ia)) deallocate(a%ia)
if (allocated(a%val)) deallocate(a%val) if (allocated(a%val)) deallocate(a%val)
call a%set_null() call a%set_null()
call a%set_nrows(0) call a%set_nrows(izero)
call a%set_ncols(0) call a%set_ncols(izero)
return return

4
base/modules/psb_d_csr_mat_mod.f90
Unescape Escape View File

@ -527,8 +527,8 @@ contains
if (allocated(a%ja)) deallocate(a%ja) if (allocated(a%ja)) deallocate(a%ja)
if (allocated(a%val)) deallocate(a%val) if (allocated(a%val)) deallocate(a%val)
call a%set_null() call a%set_null()
call a%set_nrows(0) call a%set_nrows(izero)
call a%set_ncols(0) call a%set_ncols(izero)
return return

8
base/modules/psb_d_linmap_mod.f90
Unescape Escape View File

@ -55,7 +55,7 @@ module psb_d_linmap_mod
interface psb_map_X2Y interface psb_map_X2Y
subroutine psb_d_map_X2Y(alpha,x,beta,y,map,info,work) subroutine psb_d_map_X2Y(alpha,x,beta,y,map,info,work)
use psb_const_mod use psb_const_mod
import :: psb_ipk_, psb_dlinmap_type import :: psb_dlinmap_type
implicit none implicit none
type(psb_dlinmap_type), intent(in) :: map type(psb_dlinmap_type), intent(in) :: map
real(psb_dpk_), intent(in) :: alpha,beta real(psb_dpk_), intent(in) :: alpha,beta
@ -67,7 +67,7 @@ module psb_d_linmap_mod
subroutine psb_d_map_X2Y_vect(alpha,x,beta,y,map,info,work) subroutine psb_d_map_X2Y_vect(alpha,x,beta,y,map,info,work)
use psb_const_mod use psb_const_mod
use psb_d_vect_mod use psb_d_vect_mod
import :: psb_ipk_, psb_dlinmap_type import :: psb_dlinmap_type
implicit none implicit none
type(psb_dlinmap_type), intent(in) :: map type(psb_dlinmap_type), intent(in) :: map
real(psb_dpk_), intent(in) :: alpha,beta real(psb_dpk_), intent(in) :: alpha,beta
@ -80,7 +80,7 @@ module psb_d_linmap_mod
interface psb_map_Y2X interface psb_map_Y2X
subroutine psb_d_map_Y2X(alpha,x,beta,y,map,info,work) subroutine psb_d_map_Y2X(alpha,x,beta,y,map,info,work)
use psb_const_mod use psb_const_mod
import :: psb_ipk_, psb_dlinmap_type import :: psb_dlinmap_type
implicit none implicit none
type(psb_dlinmap_type), intent(in) :: map type(psb_dlinmap_type), intent(in) :: map
real(psb_dpk_), intent(in) :: alpha,beta real(psb_dpk_), intent(in) :: alpha,beta
@ -92,7 +92,7 @@ module psb_d_linmap_mod
subroutine psb_d_map_Y2X_vect(alpha,x,beta,y,map,info,work) subroutine psb_d_map_Y2X_vect(alpha,x,beta,y,map,info,work)
use psb_const_mod use psb_const_mod
use psb_d_vect_mod use psb_d_vect_mod
import :: psb_ipk_, psb_dlinmap_type import :: psb_dlinmap_type
implicit none implicit none
type(psb_dlinmap_type), intent(in) :: map type(psb_dlinmap_type), intent(in) :: map
real(psb_dpk_), intent(in) :: alpha,beta real(psb_dpk_), intent(in) :: alpha,beta

2
base/modules/psb_d_vect_mod.f90
Unescape Escape View File

@ -163,7 +163,7 @@ contains
if (info == 0) call this%v%bld(x) if (info == 0) call this%v%bld(x)
call this%asb(size(x),info) call this%asb(size(x,kind=psb_ipk_),info)
end function constructor end function constructor

6
base/modules/psb_s_base_mat_mod.f90
Unescape Escape View File

@ -1170,9 +1170,9 @@ contains
if (allocated(a%ja)) deallocate(a%ja) if (allocated(a%ja)) deallocate(a%ja)
if (allocated(a%val)) deallocate(a%val) if (allocated(a%val)) deallocate(a%val)
call a%set_null() call a%set_null()
call a%set_nrows(0) call a%set_nrows(izero)
call a%set_ncols(0) call a%set_ncols(izero)
call a%set_nzeros(0) call a%set_nzeros(izero)
return return

2
base/modules/psb_s_base_vect_mod.f90
Unescape Escape View File

@ -150,7 +150,7 @@ contains
integer(psb_ipk_) :: info integer(psb_ipk_) :: info
this%v = x this%v = x
call this%asb(size(x),info) call this%asb(size(x,kind=psb_ipk_),info)
end function constructor end function constructor

4
base/modules/psb_s_csc_mat_mod.f90
Unescape Escape View File

@ -527,8 +527,8 @@ contains
if (allocated(a%ia)) deallocate(a%ia) if (allocated(a%ia)) deallocate(a%ia)
if (allocated(a%val)) deallocate(a%val) if (allocated(a%val)) deallocate(a%val)
call a%set_null() call a%set_null()
call a%set_nrows(0) call a%set_nrows(izero)
call a%set_ncols(0) call a%set_ncols(izero)
return return

4
base/modules/psb_s_csr_mat_mod.f90
Unescape Escape View File

@ -527,8 +527,8 @@ contains
if (allocated(a%ja)) deallocate(a%ja) if (allocated(a%ja)) deallocate(a%ja)
if (allocated(a%val)) deallocate(a%val) if (allocated(a%val)) deallocate(a%val)
call a%set_null() call a%set_null()
call a%set_nrows(0) call a%set_nrows(izero)
call a%set_ncols(0) call a%set_ncols(izero)
return return

8
base/modules/psb_s_linmap_mod.f90
Unescape Escape View File

@ -55,7 +55,7 @@ module psb_s_linmap_mod
interface psb_map_X2Y interface psb_map_X2Y
subroutine psb_s_map_X2Y(alpha,x,beta,y,map,info,work) subroutine psb_s_map_X2Y(alpha,x,beta,y,map,info,work)
use psb_const_mod use psb_const_mod
import :: psb_ipk_, psb_slinmap_type import :: psb_slinmap_type
implicit none implicit none
type(psb_slinmap_type), intent(in) :: map type(psb_slinmap_type), intent(in) :: map
real(psb_spk_), intent(in) :: alpha,beta real(psb_spk_), intent(in) :: alpha,beta
@ -67,7 +67,7 @@ module psb_s_linmap_mod
subroutine psb_s_map_X2Y_vect(alpha,x,beta,y,map,info,work) subroutine psb_s_map_X2Y_vect(alpha,x,beta,y,map,info,work)
use psb_const_mod use psb_const_mod
use psb_s_vect_mod use psb_s_vect_mod
import :: psb_ipk_, psb_slinmap_type import :: psb_slinmap_type
implicit none implicit none
type(psb_slinmap_type), intent(in) :: map type(psb_slinmap_type), intent(in) :: map
real(psb_spk_), intent(in) :: alpha,beta real(psb_spk_), intent(in) :: alpha,beta
@ -80,7 +80,7 @@ module psb_s_linmap_mod
interface psb_map_Y2X interface psb_map_Y2X
subroutine psb_s_map_Y2X(alpha,x,beta,y,map,info,work) subroutine psb_s_map_Y2X(alpha,x,beta,y,map,info,work)
use psb_const_mod use psb_const_mod
import :: psb_ipk_, psb_slinmap_type import :: psb_slinmap_type
implicit none implicit none
type(psb_slinmap_type), intent(in) :: map type(psb_slinmap_type), intent(in) :: map
real(psb_spk_), intent(in) :: alpha,beta real(psb_spk_), intent(in) :: alpha,beta
@ -92,7 +92,7 @@ module psb_s_linmap_mod
subroutine psb_s_map_Y2X_vect(alpha,x,beta,y,map,info,work) subroutine psb_s_map_Y2X_vect(alpha,x,beta,y,map,info,work)
use psb_const_mod use psb_const_mod
use psb_s_vect_mod use psb_s_vect_mod
import :: psb_ipk_, psb_slinmap_type import :: psb_slinmap_type
implicit none implicit none
type(psb_slinmap_type), intent(in) :: map type(psb_slinmap_type), intent(in) :: map
real(psb_spk_), intent(in) :: alpha,beta real(psb_spk_), intent(in) :: alpha,beta

2
base/modules/psb_s_vect_mod.f90
Unescape Escape View File

@ -163,7 +163,7 @@ contains
if (info == 0) call this%v%bld(x) if (info == 0) call this%v%bld(x)
call this%asb(size(x),info) call this%asb(size(x,kind=psb_ipk_),info)
end function constructor end function constructor

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

@ -40,6 +40,63 @@ module psb_serial_mod
& psb_gth => psi_gth,& & psb_gth => psi_gth,&
& psb_sct => psi_sct & psb_sct => psi_sct
interface psb_nrm1
module procedure psb_snrm1, psb_dnrm1, psb_cnrm1, psb_znrm1
end interface psb_nrm1
interface psb_amax
function psb_samax_s(n, x) result(val)
import :: psb_ipk_, psb_spk_, psb_dpk_
integer(psb_ipk_), intent(in) :: n
real(psb_spk_), intent(in) :: x(:)
real(psb_spk_) :: val
end function psb_samax_s
function psb_damax_s(n, x) result(val)
import :: psb_ipk_, psb_spk_, psb_dpk_
integer(psb_ipk_), intent(in) :: n
real(psb_dpk_), intent(in) :: x(:)
real(psb_dpk_) :: val
end function psb_damax_s
function psb_camax_s(n, x) result(val)
import :: psb_ipk_, psb_spk_, psb_dpk_
integer(psb_ipk_), intent(in) :: n
complex(psb_spk_), intent(in) :: x(:)
real(psb_spk_) :: val
end function psb_camax_s
function psb_zamax_s(n, x) result(val)
import :: psb_ipk_, psb_spk_, psb_dpk_
integer(psb_ipk_), intent(in) :: n
complex(psb_dpk_), intent(in) :: x(:)
real(psb_dpk_) :: val
end function psb_zamax_s
end interface psb_amax
interface psb_asum
function psb_sasum_s(n, x) result(val)
import :: psb_ipk_, psb_spk_, psb_dpk_
integer(psb_ipk_), intent(in) :: n
real(psb_spk_), intent(in) :: x(:)
real(psb_spk_) :: val
end function psb_sasum_s
function psb_dasum_s(n, x) result(val)
import :: psb_ipk_, psb_spk_, psb_dpk_
integer(psb_ipk_), intent(in) :: n
real(psb_dpk_), intent(in) :: x(:)
real(psb_dpk_) :: val
end function psb_dasum_s
function psb_casum_s(n, x) result(val)
import :: psb_ipk_, psb_spk_, psb_dpk_
integer(psb_ipk_), intent(in) :: n
complex(psb_spk_), intent(in) :: x(:)
real(psb_spk_) :: val
end function psb_casum_s
function psb_zasum_s(n, x) result(val)
import :: psb_ipk_, psb_spk_, psb_dpk_
integer(psb_ipk_), intent(in) :: n
complex(psb_dpk_), intent(in) :: x(:)
real(psb_dpk_) :: val
end function psb_zasum_s
end interface psb_asum
interface psb_symbmm interface psb_symbmm
subroutine psb_ssymbmm(a,b,c,info) subroutine psb_ssymbmm(a,b,c,info)
@ -480,6 +537,31 @@ module psb_serial_mod
contains contains
elemental function psb_snrm1(x) result(res)
real(psb_spk_), intent(in) :: x
real(psb_spk_) :: val
val = abs( x )
end function psb_snrm1
elemental function psb_dnrm1(x) result(res)
real(psb_dpk_), intent(in) :: x
real(psb_dpk_) :: val
val = abs( x )
end function psb_dnrm1
elemental function psb_cnrm1(x) result(res)
complex(psb_spk_), intent(in) :: x
real(psb_spk_) :: val
val = abs( real( x ) ) + abs( aimag( x ) )
end function psb_cnrm1
elemental function psb_znrm1(x) result(res)
complex(psb_dpk_), intent(in) :: x
real(psb_dpk_) :: val
val = abs( real( x ) ) + abs( aimag( x ) )
end function psb_znrm1
subroutine psb_scsprt(iout,a,iv,head,ivr,ivc) subroutine psb_scsprt(iout,a,iv,head,ivr,ivc)
use psb_mat_mod use psb_mat_mod
integer(psb_ipk_), intent(in) :: iout integer(psb_ipk_), intent(in) :: iout

6
base/modules/psb_z_base_mat_mod.f90
Unescape Escape View File

@ -1170,9 +1170,9 @@ contains
if (allocated(a%ja)) deallocate(a%ja) if (allocated(a%ja)) deallocate(a%ja)
if (allocated(a%val)) deallocate(a%val) if (allocated(a%val)) deallocate(a%val)
call a%set_null() call a%set_null()
call a%set_nrows(0) call a%set_nrows(izero)
call a%set_ncols(0) call a%set_ncols(izero)
call a%set_nzeros(0) call a%set_nzeros(izero)
return return

2
base/modules/psb_z_base_vect_mod.f90
Unescape Escape View File

@ -150,7 +150,7 @@ contains
integer(psb_ipk_) :: info integer(psb_ipk_) :: info
this%v = x this%v = x
call this%asb(size(x),info) call this%asb(size(x,kind=psb_ipk_),info)
end function constructor end function constructor

4
base/modules/psb_z_csc_mat_mod.f90
Unescape Escape View File

@ -527,8 +527,8 @@ contains
if (allocated(a%ia)) deallocate(a%ia) if (allocated(a%ia)) deallocate(a%ia)
if (allocated(a%val)) deallocate(a%val) if (allocated(a%val)) deallocate(a%val)
call a%set_null() call a%set_null()
call a%set_nrows(0) call a%set_nrows(izero)
call a%set_ncols(0) call a%set_ncols(izero)
return return

4
base/modules/psb_z_csr_mat_mod.f90
Unescape Escape View File

@ -527,8 +527,8 @@ contains
if (allocated(a%ja)) deallocate(a%ja) if (allocated(a%ja)) deallocate(a%ja)
if (allocated(a%val)) deallocate(a%val) if (allocated(a%val)) deallocate(a%val)
call a%set_null() call a%set_null()
call a%set_nrows(0) call a%set_nrows(izero)
call a%set_ncols(0) call a%set_ncols(izero)
return return

8
base/modules/psb_z_linmap_mod.f90
Unescape Escape View File

@ -55,7 +55,7 @@ module psb_z_linmap_mod
interface psb_map_X2Y interface psb_map_X2Y
subroutine psb_z_map_X2Y(alpha,x,beta,y,map,info,work) subroutine psb_z_map_X2Y(alpha,x,beta,y,map,info,work)
use psb_const_mod use psb_const_mod
import :: psb_ipk_, psb_zlinmap_type import :: psb_zlinmap_type
implicit none implicit none
type(psb_zlinmap_type), intent(in) :: map type(psb_zlinmap_type), intent(in) :: map
complex(psb_dpk_), intent(in) :: alpha,beta complex(psb_dpk_), intent(in) :: alpha,beta
@ -67,7 +67,7 @@ module psb_z_linmap_mod
subroutine psb_z_map_X2Y_vect(alpha,x,beta,y,map,info,work) subroutine psb_z_map_X2Y_vect(alpha,x,beta,y,map,info,work)
use psb_const_mod use psb_const_mod
use psb_z_vect_mod use psb_z_vect_mod
import :: psb_ipk_, psb_zlinmap_type import :: psb_zlinmap_type
implicit none implicit none
type(psb_zlinmap_type), intent(in) :: map type(psb_zlinmap_type), intent(in) :: map
complex(psb_dpk_), intent(in) :: alpha,beta complex(psb_dpk_), intent(in) :: alpha,beta
@ -80,7 +80,7 @@ module psb_z_linmap_mod
interface psb_map_Y2X interface psb_map_Y2X
subroutine psb_z_map_Y2X(alpha,x,beta,y,map,info,work) subroutine psb_z_map_Y2X(alpha,x,beta,y,map,info,work)
use psb_const_mod use psb_const_mod
import :: psb_ipk_, psb_zlinmap_type import :: psb_zlinmap_type
implicit none implicit none
type(psb_zlinmap_type), intent(in) :: map type(psb_zlinmap_type), intent(in) :: map
complex(psb_dpk_), intent(in) :: alpha,beta complex(psb_dpk_), intent(in) :: alpha,beta
@ -92,7 +92,7 @@ module psb_z_linmap_mod
subroutine psb_z_map_Y2X_vect(alpha,x,beta,y,map,info,work) subroutine psb_z_map_Y2X_vect(alpha,x,beta,y,map,info,work)
use psb_const_mod use psb_const_mod
use psb_z_vect_mod use psb_z_vect_mod
import :: psb_ipk_, psb_zlinmap_type import :: psb_zlinmap_type
implicit none implicit none
type(psb_zlinmap_type), intent(in) :: map type(psb_zlinmap_type), intent(in) :: map
complex(psb_dpk_), intent(in) :: alpha,beta complex(psb_dpk_), intent(in) :: alpha,beta

2
base/modules/psb_z_vect_mod.f90
Unescape Escape View File

@ -163,7 +163,7 @@ contains
if (info == 0) call this%v%bld(x) if (info == 0) call this%v%bld(x)
call this%asb(size(x),info) call this%asb(size(x,kind=psb_ipk_),info)
end function constructor end function constructor

125
base/psblas/psb_camax.f90
Unescape Escape View File

@ -44,7 +44,7 @@
! info - integer. Return code ! info - integer. Return code
! jx - integer(optional). The column offset. ! jx - integer(optional). The column offset.
! !
function psb_camax(x,desc_a, info, jx) function psb_camax(x,desc_a, info, jx) result(res)
use psb_base_mod, psb_protect_name => psb_camax use psb_base_mod, psb_protect_name => psb_camax
implicit none implicit none
@ -53,25 +53,20 @@ function psb_camax(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(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: jx integer(psb_ipk_), optional, intent(in) :: jx
real(psb_spk_) :: psb_camax real(psb_spk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ix, ijx, m, imax, icamax & err_act, iix, jjx, ix, ijx, m, ldx
real(psb_spk_) :: amax character(len=20) :: name, ch_err
character(len=20) :: name, ch_err
complex :: cdum
real :: cabs1
cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
name='psb_camax' name='psb_camax'
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)
amax=0.d0
ictxt=desc_a%get_context() ictxt = desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if (np == -1) then if (np == -1) then
@ -88,8 +83,9 @@ function psb_camax(x,desc_a, info, jx)
endif endif
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = size(x,1)
call psb_chkvect(m,1,size(x,1),ix,ijx,desc_a,info,iix,jjx) 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'
@ -105,16 +101,13 @@ function psb_camax(x,desc_a, info, jx)
! compute local max ! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
imax=icamax(desc_a%get_local_rows()-iix+1,x(iix,jjx),1) res = psb_amax(desc_a%get_local_rows()-iix+1,x(:,jjx))
amax=cabs1(x(iix+imax-1,jjx))
else else
amax = szero res = szero
end if end if
! compute global max ! compute global max
call psb_amx(ictxt, amax) call psb_amx(ictxt, res)
psb_camax=amax
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -174,32 +167,27 @@ end function psb_camax
! 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_camaxv (x,desc_a, info) function psb_camaxv (x,desc_a, info) result(res)
use psb_base_mod, psb_protect_name => psb_camaxv use psb_base_mod, psb_protect_name => psb_camaxv
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_camaxv real(psb_spk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, jx, ix, m, imax, icamax & err_act, iix, jjx, jx, ix, m, ldx
real(psb_spk_) :: amax
complex(psb_spk_) :: cmax
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
complex :: cdum
real :: cabs1
cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
name='psb_camaxv' name='psb_camaxv'
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)
amax=0.d0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
@ -214,8 +202,9 @@ function psb_camaxv (x,desc_a, info)
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,1),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'
@ -231,17 +220,13 @@ function psb_camaxv (x,desc_a, info)
! compute local max ! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
imax=icamax(desc_a%get_local_rows()-iix+1,x(iix),1) res = psb_amax(desc_a%get_local_rows()-iix+1,x)
cmax=(x(iix+imax-1))
amax=cabs1(cmax)
else else
amax = szero res = szero
end if end if
! compute global max ! compute global max
call psb_amx(ictxt, amax) call psb_amx(ictxt, res)
psb_camaxv=amax
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -269,12 +254,11 @@ function psb_camax_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, jx, ix, m, imax, isamax & err_act, iix, jjx, jx, ix, m
real(psb_spk_) :: amax
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_camaxv' name='psb_camaxv'
@ -282,7 +266,6 @@ function psb_camax_vect(x, desc_a, info) result(res)
info=psb_success_ info=psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
amax=szero
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
@ -302,8 +285,7 @@ function psb_camax_vect(x, desc_a, info) result(res)
jx = 1 jx = 1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
call psb_chkvect(m,ione,x%get_nrows(),ix,jx,desc_a,info,iix,jjx)
call psb_chkvect(m,1,x%get_nrows(),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'
@ -319,15 +301,13 @@ function psb_camax_vect(x, desc_a, info) result(res)
! compute local max ! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
amax=x%amax(desc_a%get_local_rows()) res = x%amax(desc_a%get_local_rows())
else else
amax = szero res = szero
end if end if
! compute global max ! compute global max
call psb_amx(ictxt, amax) call psb_amx(ictxt, res)
res=amax
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -397,27 +377,21 @@ subroutine psb_camaxvs(res,x,desc_a, info)
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_), intent(out) :: res real(psb_spk_), intent(out) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ix, ijx, m, imax, icamax & err_act, iix, jjx, ix, ijx, m, ldx
real(psb_spk_) :: amax character(len=20) :: name, ch_err
character(len=20) :: name, ch_err
complex(psb_spk_) :: cmax
complex :: cdum
real :: cabs1
cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
name='psb_camaxvs' name='psb_camaxvs'
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)
amax=0.d0
ictxt=desc_a%get_context() ictxt = desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if (np == -1) then if (np == -1) then
@ -430,7 +404,8 @@ subroutine psb_camaxvs(res,x,desc_a, info)
ijx=1 ijx=1
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'
@ -446,17 +421,13 @@ subroutine psb_camaxvs(res,x,desc_a, info)
! compute local max ! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
imax=icamax(desc_a%get_local_rows()-iix+1,x(iix),1) res = psb_amax(desc_a%get_local_rows()-iix+1,x)
cmax=(x(iix+imax-1))
amax=cabs1(cmax)
else else
amax = szero res = szero
end if end if
! compute global max ! compute global max
call psb_amx(ictxt, amax) call psb_amx(ictxt, res)
res = amax
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -528,21 +499,14 @@ subroutine psb_cmamaxs(res,x,desc_a, info,jx)
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ix, ijx, m, imax, i, k, icamax & err_act, iix, jjx, ix, ijx, m, ldx, i, k
real(psb_spk_) :: amax
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
complex(psb_spk_) :: cmax
complex :: cdum
real :: cabs1
cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
name='psb_cmamaxs' name='psb_cmamaxs'
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)
amax=0.d0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
@ -561,8 +525,8 @@ subroutine psb_cmamaxs(res,x,desc_a, info,jx)
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
k = min(size(x,2),size(res,1)) k = min(size(x,2),size(res,1))
ldx = size(x,1)
call psb_chkvect(m,1,size(x,1),ix,ijx,desc_a,info,iix,jjx) 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'
@ -576,15 +540,12 @@ subroutine psb_cmamaxs(res,x,desc_a, info,jx)
goto 9999 goto 9999
end if end if
res(1:k) = szero
! compute local max ! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
do i=1,k do i=1,k
imax=icamax(desc_a%get_local_rows()-iix+1,x(iix,jjx+i-1),1) res(i) = psb_amax(desc_a%get_local_rows()-iix+1,x(:,jjx+i-1))
cmax=(x(iix+imax-1,jjx+i-1)) end do
res(i)=cabs1(cmax)
end do
else
amax = szero
end if end if
! compute global max ! compute global max

183
base/psblas/psb_casum.f90
Unescape Escape View File

@ -44,7 +44,7 @@
! info - integer. Return code ! info - integer. Return code
! jx - integer(optional). The column offset. ! jx - integer(optional). The column offset.
! !
function psb_casum (x,desc_a, info, jx) function psb_casum (x,desc_a, info, jx) result(res)
use psb_base_mod, psb_protect_name => psb_casum use psb_base_mod, psb_protect_name => psb_casum
implicit none implicit none
@ -53,24 +53,18 @@ function psb_casum (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(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: jx integer(psb_ipk_), optional, intent(in) :: jx
real(psb_spk_) :: psb_casum real(psb_spk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me, & integer(psb_ipk_) :: ictxt, np, me, &
& err_act, iix, jjx, ix, ijx, m, i, idx, ndm & err_act, iix, jjx, ix, ijx, m, i, idx, ndm, ldx
real(psb_spk_) :: asum, scasum
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
complex(psb_spk_) :: cmax
complex :: cdum
real :: cabs1
cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
name='psb_casum' name='psb_casum'
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)
asum=0.d0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
@ -89,9 +83,9 @@ function psb_casum (x,desc_a, info, jx)
endif endif
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = size(x,1)
! check vector correctness ! check vector correctness
call psb_chkvect(m,1,size(x,1),ix,ijx,desc_a,info,iix,jjx) 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'
@ -106,29 +100,21 @@ function psb_casum (x,desc_a, info, jx)
end if end if
! compute local max ! compute local max
if ((m /= 0)) then if(desc_a%get_local_rows() > 0) then
if(desc_a%get_local_rows() > 0) then res = psb_asum(desc_a%get_local_rows()-iix+1,x(:,jjx))
asum=scasum(desc_a%get_local_rows()-iix+1,x(iix:,jjx),ione)
! adjust asum because overlapped elements are computed more than once
do i=1,size(desc_a%ovrlap_elem,1)
idx = desc_a%ovrlap_elem(i,1)
ndm = desc_a%ovrlap_elem(i,2)
asum = asum - (real(ndm-1)/real(ndm))*cabs1(x(idx,jjx))
end do
else
asum=0.0
end if
! compute global sum
call psb_sum(ictxt, asum)
else
asum=0.0
end if
! adjust res because overlapped elements are computed more than once
do i=1,size(desc_a%ovrlap_elem,1)
idx = desc_a%ovrlap_elem(i,1)
ndm = desc_a%ovrlap_elem(i,2)
res = res - (real(ndm-1)/real(ndm))*psb_nrm1(x(idx,jjx))
end do
psb_casum=asum else
res = szero
end if
! compute global sum
call psb_sum(ictxt, res)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -137,31 +123,25 @@ function psb_casum (x,desc_a, info, jx)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then if (err_act == psb_act_abort_) then
call psb_error(ictxt) call psb_error(ictxt)
return return
end if end if
return return
end function psb_casum end function psb_casum
function psb_casum_vect(x, desc_a, info) result(res) function psb_casum_vect(x, desc_a, info) result(res)
use psb_penv_mod use psb_base_mod, psb_protect_name => psb_casum_vect
use psb_serial_mod
use psb_descriptor_type
use psb_check_mod
use psb_error_mod
use psb_c_vect_mod
implicit none implicit none
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, jx, ix, m, imax & err_act, iix, jjx, jx, ix, m, imax
real(psb_spk_) :: asum
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_casumv' name='psb_casumv'
@ -169,7 +149,6 @@ function psb_casum_vect(x, desc_a, info) result(res)
info=psb_success_ info=psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
asum=0.d0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
@ -191,8 +170,7 @@ function psb_casum_vect(x, desc_a, info) result(res)
jx = 1 jx = 1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
call psb_chkvect(m,ione,x%get_nrows(),ix,jx,desc_a,info,iix,jjx)
call psb_chkvect(m,1,x%get_nrows(),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'
@ -208,15 +186,13 @@ function psb_casum_vect(x, desc_a, info) result(res)
! compute local max ! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
asum=x%asum(desc_a%get_local_rows()) res = x%asum(desc_a%get_local_rows())
else else
asum = szero res = szero
end if end if
! compute global sum ! compute global sum
call psb_sum(ictxt, asum) call psb_sum(ictxt, res)
res=asum
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -276,33 +252,26 @@ end function psb_casum_vect
! 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_casumv(x,desc_a, info) function psb_casumv(x,desc_a, info) result(res)
use psb_base_mod, psb_protect_name => psb_casumv use psb_base_mod, psb_protect_name => psb_casumv
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_casumv real(psb_spk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, jx, ix, m, i, idx, ndm & err_act, iix, jjx, jx, ix, m, i, idx, ndm, ldx
real(psb_spk_) :: asum, scasum
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
complex(psb_spk_) :: cmax
complex :: cdum
real :: cabs1
cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
name='psb_casumv' name='psb_casumv'
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)
asum=0.d0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
@ -316,9 +285,9 @@ function psb_casumv(x,desc_a, info)
jx=1 jx=1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = size(x,1)
! check vector correctness ! check vector correctness
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'
@ -333,28 +302,22 @@ function psb_casumv(x,desc_a, info)
end if end if
! compute local max ! compute local max
if ((m /= 0)) then if(desc_a%get_local_rows() > 0) then
if(desc_a%get_local_rows() > 0) then res = psb_asum(desc_a%get_local_rows(),x)
asum=scasum(desc_a%get_local_rows(),x,ione)
! adjust asum because overlapped elements are computed more than once
! adjust asum 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) res = res - (real(ndm-1)/real(ndm))*psb_nrm1(x(idx))
asum = asum - (real(ndm-1)/real(ndm))*cabs1(x(idx)) end do
end do
else
asum=0.d0
end if
! compute global sum
call psb_sum(ictxt, asum)
else else
asum=0.d0 res = szero
end if end if
psb_casumv=asum ! compute global sum
call psb_sum(ictxt, res)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -419,28 +382,21 @@ subroutine psb_casumvs(res,x,desc_a, info)
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, ix, jx, m, i, idx, ndm & err_act, iix, jjx, ix, jx, m, i, idx, ndm, ldx
real(psb_spk_) :: asum, scasum
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
complex(psb_spk_) :: cmax
complex :: cdum
real :: cabs1
cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
name='psb_casumvs' name='psb_casumvs'
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)
asum=0.d0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
@ -454,9 +410,9 @@ subroutine psb_casumvs(res,x,desc_a, info)
jx = 1 jx = 1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = size(x,1)
! check vector correctness ! check vector correctness
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'
@ -471,29 +427,22 @@ subroutine psb_casumvs(res,x,desc_a, info)
end if end if
! compute local max ! compute local max
if ((m /= 0)) then if(desc_a%get_local_rows() > 0) then
if(desc_a%get_local_rows() > 0) then res = psb_asum(desc_a%get_local_rows(),x)
asum=scasum(desc_a%get_local_rows(),x,ione)
! adjust asum because overlapped elements are computed more than once
! adjust asum 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) res = res - (real(ndm-1)/real(ndm))*psb_nrm1(x(idx))
asum = asum - (real(ndm-1)/real(ndm))*cabs1(x(idx)) end do
end do
else
asum=0.d0
end if
! compute global sum
call psb_sum(ictxt,asum)
else else
asum=0.d0 res = szero
end if end if
! compute global sum
res = asum call psb_sum(ictxt,res)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

249
base/psblas/psb_damax.f90
Unescape Escape View File

@ -39,39 +39,34 @@
! where sub( X ) denotes X(1:N,JX:). ! where sub( X ) denotes X(1:N,JX:).
! !
! Arguments: ! Arguments:
! x(:,:) - real The input vector. ! x(:,:) - real The input vector.
! 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. ! jx - integer(optional). The column offset.
! !
function psb_damax (x,desc_a, info, jx) function psb_damax(x,desc_a, info, jx) result(res)
use psb_penv_mod use psb_base_mod, psb_protect_name => psb_damax
use psb_serial_mod
use psb_descriptor_type
use psb_check_mod
use psb_error_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
integer(psb_ipk_), optional, intent(in) :: jx integer(psb_ipk_), optional, intent(in) :: jx
real(psb_dpk_) :: psb_damax real(psb_dpk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ix, ijx, m, imax, idamax & err_act, iix, jjx, ix, ijx, m, ldx
real(psb_dpk_) :: amax character(len=20) :: name, ch_err
character(len=20) :: name, ch_err
name='psb_damax' name='psb_damax'
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)
amax=0.d0
ictxt=desc_a%get_context() ictxt = desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if (np == -1) then if (np == -1) then
@ -82,39 +77,37 @@ function psb_damax (x,desc_a, info, jx)
ix = 1 ix = 1
if (present(jx)) then if (present(jx)) then
ijx = jx ijx = jx
else else
ijx = 1 ijx = 1
endif endif
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = size(x,1)
call psb_chkvect(m,1,size(x,1),ix,ijx,desc_a,info,iix,jjx) 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'
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
if (iix /= 1) then if (iix /= 1) then
info=psb_err_ix_n1_iy_n1_unsupported_ info=psb_err_ix_n1_iy_n1_unsupported_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
! compute local max ! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
imax=idamax(desc_a%get_local_rows()-iix+1,x(iix,jjx),1) res = psb_amax(desc_a%get_local_rows()-iix+1,x(:,jjx))
amax=abs(x(iix+imax-1,jjx))
else else
amax = dzero res = dzero
end if end if
! compute global max ! compute global max
call psb_amx(ictxt, amax) call psb_amx(ictxt, res)
psb_damax=amax
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -123,8 +116,8 @@ function psb_damax (x,desc_a, info, jx)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then if (err_act == psb_act_abort_) then
call psb_error(ictxt) call psb_error(ictxt)
return return
end if end if
return return
end function psb_damax end function psb_damax
@ -170,27 +163,24 @@ end function psb_damax
! normi := max(abs(X(i)) ! normi := max(abs(X(i))
! !
! Arguments: ! Arguments:
! x(:) - real The input vector. ! x(:) - real The input vector.
! 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_damaxv (x,desc_a, info) function psb_damaxv (x,desc_a, info) result(res)
use psb_penv_mod use psb_base_mod, psb_protect_name => psb_damaxv
use psb_serial_mod
use psb_descriptor_type
use psb_check_mod
use psb_error_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_damaxv real(psb_dpk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, jx, ix, m, imax, idamax & err_act, iix, jjx, jx, ix, m, ldx
real(psb_dpk_) :: amax
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_damaxv' name='psb_damaxv'
@ -198,7 +188,6 @@ function psb_damaxv (x,desc_a, info)
info=psb_success_ info=psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
amax=0.d0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
@ -213,33 +202,31 @@ function psb_damaxv (x,desc_a, info)
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,1),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'
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
if (iix /= 1) then if (iix /= 1) then
info=psb_err_ix_n1_iy_n1_unsupported_ info=psb_err_ix_n1_iy_n1_unsupported_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
! compute local max ! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
imax=idamax(desc_a%get_local_rows()-iix+1,x(iix),1) res = psb_amax(desc_a%get_local_rows()-iix+1,x)
amax=abs(x(iix+imax-1))
else else
amax = dzero res = dzero
end if end if
! compute global max ! compute global max
call psb_amx(ictxt, amax) call psb_amx(ictxt, res)
psb_damaxv=amax
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -248,12 +235,13 @@ function psb_damaxv (x,desc_a, info)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then if (err_act == psb_act_abort_) then
call psb_error(ictxt) call psb_error(ictxt)
return return
end if end if
return return
end function psb_damaxv end function psb_damaxv
function psb_damax_vect(x, desc_a, info) result(res) function psb_damax_vect(x, desc_a, info) result(res)
use psb_penv_mod use psb_penv_mod
use psb_serial_mod use psb_serial_mod
@ -266,20 +254,18 @@ function psb_damax_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, jx, ix, m, imax, idamax & err_act, iix, jjx, jx, ix, m
real(psb_dpk_) :: amax character(len=20) :: name, ch_err
character(len=20) :: name, ch_err
name='psb_damaxv' name='psb_damaxv'
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)
amax=0.d0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
@ -299,8 +285,7 @@ function psb_damax_vect(x, desc_a, info) result(res)
jx = 1 jx = 1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
call psb_chkvect(m,ione,x%get_nrows(),ix,jx,desc_a,info,iix,jjx)
call psb_chkvect(m,1,x%get_nrows(),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'
@ -316,15 +301,13 @@ function psb_damax_vect(x, desc_a, info) result(res)
! compute local max ! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
amax=x%amax(desc_a%get_local_rows()) res = x%amax(desc_a%get_local_rows())
else else
amax = dzero res = dzero
end if end if
! compute global max ! compute global max
call psb_amx(ictxt, amax) call psb_amx(ictxt, res)
res=amax
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -381,35 +364,34 @@ end function psb_damax_vect
! where sub( X ) denotes X(1:N,JX:). ! where sub( X ) denotes X(1:N,JX:).
! !
! Arguments: ! Arguments:
! res - real. The result. ! res - real The result.
! x(:,:) - real The input vector. ! x(:,:) - real The input vector.
! 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. ! jx - integer(optional). The column offset.
! !
subroutine psb_damaxvs (res,x,desc_a, info) subroutine psb_damaxvs(res,x,desc_a, info)
use psb_base_mod, psb_protect_name => psb_damaxvs use psb_base_mod, psb_protect_name => psb_damaxvs
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_), intent(out) :: res real(psb_dpk_), intent(out) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ix, ijx, m, imax, idamax & err_act, iix, jjx, ix, ijx, m, ldx
real(psb_dpk_) :: amax character(len=20) :: name, ch_err
character(len=20) :: name, ch_err
name='psb_damaxvs' name='psb_damaxvs'
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)
amax=0.d0
ictxt=desc_a%get_context() ictxt = desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if (np == -1) then if (np == -1) then
@ -422,33 +404,30 @@ subroutine psb_damaxvs (res,x,desc_a, info)
ijx=1 ijx=1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx=size(x,1)
call psb_chkvect(m,1,size(x,1),ix,ijx,desc_a,info,iix,jjx) 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'
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
if (iix /= 1) then if (iix /= 1) then
info=psb_err_ix_n1_iy_n1_unsupported_ info=psb_err_ix_n1_iy_n1_unsupported_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
! compute local max ! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
imax=idamax(desc_a%get_local_rows()-iix+1,x(iix),1) res = psb_amax(desc_a%get_local_rows()-iix+1,x)
amax=abs(x(iix+imax-1))
else else
amax = dzero res = dzero
end if end if
! compute global max ! compute global max
call psb_amx(ictxt, amax) call psb_amx(ictxt, res)
res = amax
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -457,8 +436,8 @@ subroutine psb_damaxvs (res,x,desc_a, info)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then if (err_act == psb_act_abort_) then
call psb_error(ictxt) call psb_error(ictxt)
return return
end if end if
return return
end subroutine psb_damaxvs end subroutine psb_damaxvs
@ -502,34 +481,32 @@ end subroutine psb_damaxvs
! normi := max(abs(X(i)) ! normi := max(abs(X(i))
! !
! Arguments: ! Arguments:
! res(:) - real The result. ! res(:) - real. The result.
! x(:,:) - real The input vector. ! x(:,:) - real The input vector.
! 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_dmamaxs (res,x,desc_a, info,jx) subroutine psb_dmamaxs(res,x,desc_a, info,jx)
use psb_base_mod, psb_protect_name => psb_dmamaxs use psb_base_mod, psb_protect_name => psb_dmamaxs
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
integer(psb_ipk_), optional, intent(in) :: jx integer(psb_ipk_), optional, intent(in) :: jx
real(psb_dpk_), intent(out) :: res(:) real(psb_dpk_), intent(out) :: res(:)
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ix, ijx, m, imax, i, k, idamax & err_act, iix, jjx, ix, ijx, m, ldx, i, k
real(psb_dpk_) :: amax
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_dmamaxs' name='psb_dmamaxs'
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)
amax=0.d0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
@ -548,8 +525,8 @@ subroutine psb_dmamaxs (res,x,desc_a, info,jx)
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
k = min(size(x,2),size(res,1)) k = min(size(x,2),size(res,1))
ldx = size(x,1)
call psb_chkvect(m,1,size(x,1),ix,ijx,desc_a,info,iix,jjx) 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'
@ -563,14 +540,12 @@ subroutine psb_dmamaxs (res,x,desc_a, info,jx)
goto 9999 goto 9999
end if end if
res(1:k) = dzero
! compute local max ! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
do i=1,k do i=1,k
imax=idamax(desc_a%get_local_rows()-iix+1,x(iix,jjx+i-1),1) res(i) = psb_amax(desc_a%get_local_rows()-iix+1,x(:,jjx+i-1))
res(i)=abs(x(iix+imax-1,jjx+i-1)) end do
end do
else
amax = dzero
end if end if
! compute global max ! compute global max

355
base/psblas/psb_dasum.f90
Unescape Escape View File

@ -39,38 +39,32 @@
! where sub( X ) denotes X(1:N,JX:). ! where sub( X ) denotes X(1:N,JX:).
! !
! Arguments: ! Arguments:
! x(:,:) - real The input vector. ! x(:,:) - real The input vector.
! 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. ! jx - integer(optional). The column offset.
! !
function psb_dasum (x,desc_a, info, jx) function psb_dasum (x,desc_a, info, jx) result(res)
use psb_base_mod, psb_protect_name => psb_dasum
use psb_serial_mod
use psb_descriptor_type
use psb_check_mod
use psb_error_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
integer(psb_ipk_), optional, intent(in) :: jx integer(psb_ipk_), optional, intent(in) :: jx
real(psb_dpk_) :: psb_dasum real(psb_dpk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me, err_act, & integer(psb_ipk_) :: ictxt, np, me, &
& iix, jjx, ix, ijx, m, i, idx, ndm & err_act, iix, jjx, ix, ijx, m, i, idx, ndm, ldx
real(psb_dpk_) :: asum, dasum character(len=20) :: name, ch_err
character(len=20) :: name, ch_err
name='psb_dasum' name='psb_dasum'
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)
asum=0.d0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
@ -89,9 +83,9 @@ function psb_dasum (x,desc_a, info, jx)
endif endif
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = size(x,1)
! check vector correctness ! check vector correctness
call psb_chkvect(m,1,size(x,1),ix,ijx,desc_a,info,iix,jjx) 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'
@ -106,31 +100,21 @@ function psb_dasum (x,desc_a, info, jx)
end if end if
! compute local max ! compute local max
if ((m /= 0)) then if(desc_a%get_local_rows() > 0) then
if(desc_a%get_local_rows() > 0) then res = psb_asum(desc_a%get_local_rows()-iix+1,x(:,jjx))
asum=dasum(desc_a%get_local_rows()-iix+1,x(iix:,jjx),ione)
! adjust asum because overlapped elements are computed more than once
do i=1,size(desc_a%ovrlap_elem,1)
idx = desc_a%ovrlap_elem(i,1)
ndm = desc_a%ovrlap_elem(i,2)
asum = asum - (real(ndm-1)/real(ndm))*abs(x(idx,jjx))
end do
! compute global sum
call psb_sum(ictxt, asum)
else
asum=0.d0
! compute global sum
call psb_sum(ictxt, asum)
end if
else
asum=0.d0
end if
! adjust res because overlapped elements are computed more than once
do i=1,size(desc_a%ovrlap_elem,1)
idx = desc_a%ovrlap_elem(i,1)
ndm = desc_a%ovrlap_elem(i,2)
res = res - (real(ndm-1)/real(ndm))*psb_nrm1(x(idx,jjx))
end do
psb_dasum=asum else
res = dzero
end if
! compute global sum
call psb_sum(ictxt, res)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -146,73 +130,25 @@ function psb_dasum (x,desc_a, info, jx)
end function psb_dasum end function psb_dasum
!!$ function psb_dasum_vect(x, desc_a, info) result(res)
!!$ Parallel Sparse BLAS version 3.0 use psb_base_mod, psb_protect_name => psb_dasum_vect
!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
!
! Function: psb_dasumv
! Computes norm1 of X
!
! norm1 := sum(X(i))
!
! Arguments:
! x(:) - real The input vector.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
!
function psb_dasumv (x,desc_a, info)
use psb_serial_mod
use psb_descriptor_type
use psb_check_mod
use psb_error_mod
use psb_penv_mod
implicit none implicit none
real(psb_dpk_), intent(in) :: x(:) real(psb_dpk_) :: res
type(psb_desc_type), intent(in) :: desc_a type(psb_d_vect_type), intent (inout) :: x
integer(psb_ipk_), intent(out) :: info type(psb_desc_type), intent (in) :: desc_a
real(psb_dpk_) :: psb_dasumv integer(psb_ipk_), intent(out) :: info
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me, err_act, iix, jjx, jx, ix, m, i, idx, ndm integer(psb_ipk_) :: ictxt, np, me,&
real(psb_dpk_) :: asum, dasum & err_act, iix, jjx, jx, ix, m, imax
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_dasumv' name='psb_dasumv'
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)
asum=0.d0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
@ -223,13 +159,18 @@ function psb_dasumv (x,desc_a, info)
goto 9999 goto 9999
endif endif
if (.not.allocated(x%v)) then
info = psb_err_invalid_vect_state_
call psb_errpush(info,name)
goto 9999
endif
ix = 1 ix = 1
jx=1 jx = 1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
call psb_chkvect(m,ione,x%get_nrows(),ix,jx,desc_a,info,iix,jjx)
! check vector correctness
call psb_chkvect(m,1,size(x),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'
@ -244,29 +185,14 @@ function psb_dasumv (x,desc_a, info)
end if end if
! compute local max ! compute local max
if ((m /= 0)) then if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
if(desc_a%get_local_rows() > 0) then res = x%asum(desc_a%get_local_rows())
asum=dasum(desc_a%get_local_rows(),x,ione)
! adjust asum because overlapped elements are computed more than once
do i=1,size(desc_a%ovrlap_elem,1)
idx = desc_a%ovrlap_elem(i,1)
ndm = desc_a%ovrlap_elem(i,2)
asum = asum - (real(ndm-1)/real(ndm))*abs(x(idx))
end do
! compute global sum
call psb_sum(ictxt, asum)
else
asum=0.d0
! compute global sum
call psb_sum(ictxt, asum)
end if
else else
asum=0.d0 res = dzero
end if end if
psb_dasumv=asum ! compute global sum
call psb_sum(ictxt, res)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -279,36 +205,73 @@ function psb_dasumv (x,desc_a, info)
return return
end if end if
return return
end function psb_dasumv
end function psb_dasum_vect
!!$
!!$ Parallel Sparse BLAS version 3.0
!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
!
! Function: psb_dasumv
! Computes norm1 of X
!
! norm1 := sum(X(i))
!
! Arguments:
! x(:) - real The input vector.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
!
function psb_dasumv(x,desc_a, info) result(res)
use psb_base_mod, psb_protect_name => psb_dasumv
function psb_dasum_vect(x, desc_a, info) result(res)
use psb_penv_mod
use psb_serial_mod
use psb_descriptor_type
use psb_check_mod
use psb_error_mod
use psb_d_vect_mod
implicit none implicit none
real(psb_dpk_) :: res real(psb_dpk_), intent(in) :: 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 real(psb_dpk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, jx, ix, m, imax & err_act, iix, jjx, jx, ix, m, i, idx, ndm, ldx
real(psb_dpk_) :: asum
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_dasumv' name='psb_dasumv'
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)
asum=0.d0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
@ -318,19 +281,13 @@ function psb_dasum_vect(x, desc_a, info) result(res)
goto 9999 goto 9999
endif endif
if (.not.allocated(x%v)) then
info = psb_err_invalid_vect_state_
call psb_errpush(info,name)
goto 9999
endif
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,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) ! check vector correctness
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'
@ -345,16 +302,22 @@ function psb_dasum_vect(x, desc_a, info) result(res)
end if end if
! compute local max ! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then if(desc_a%get_local_rows() > 0) then
asum=x%asum(desc_a%get_local_rows()) res = psb_asum(desc_a%get_local_rows(),x)
! adjust asum because overlapped elements are computed more than once
do i=1,size(desc_a%ovrlap_elem,1)
idx = desc_a%ovrlap_elem(i,1)
ndm = desc_a%ovrlap_elem(i,2)
res = res - (real(ndm-1)/real(ndm))*psb_nrm1(x(idx))
end do
else else
asum = dzero res = dzero
end if end if
! compute global sum ! compute global sum
call psb_sum(ictxt, asum) call psb_sum(ictxt, res)
res=asum
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -367,9 +330,7 @@ function psb_dasum_vect(x, desc_a, info) result(res)
return return
end if end if
return return
end function psb_dasumv
end function psb_dasum_vect
!!$ !!$
@ -410,33 +371,32 @@ end function psb_dasum_vect
! norm1 := sum(X(i)) ! norm1 := sum(X(i))
! !
! Arguments: ! Arguments:
! res - real The result. ! res - real. The result.
! x(:) - real The input vector. ! x(:) - real The input vector.
! 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. ! jx - integer(optional). The column offset.
! !
subroutine psb_dasumvs(res,x,desc_a, info) subroutine psb_dasumvs(res,x,desc_a, info)
use psb_base_mod, psb_protect_name => psb_dasumvs use psb_base_mod, psb_protect_name => psb_dasumvs
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, err_act, iix, jjx, ix, jx, m, i, idx, ndm integer(psb_ipk_) :: ictxt, np, me,&
real(psb_dpk_) :: asum, dasum & err_act, iix, jjx, ix, jx, m, i, idx, ndm, ldx
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_dasumvs' name='psb_dasumvs'
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)
asum=0.d0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
@ -450,9 +410,9 @@ subroutine psb_dasumvs(res,x,desc_a, info)
jx = 1 jx = 1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = size(x,1)
! check vector correctness ! check vector correctness
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'
@ -467,31 +427,22 @@ subroutine psb_dasumvs(res,x,desc_a, info)
end if end if
! compute local max ! compute local max
if ((m /= 0)) then if(desc_a%get_local_rows() > 0) then
if(desc_a%get_local_rows() > 0) then res = psb_asum(desc_a%get_local_rows(),x)
asum=dasum(desc_a%get_local_rows(),x,ione)
! adjust asum because overlapped elements are computed more than once
! adjust asum 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) res = res - (real(ndm-1)/real(ndm))*psb_nrm1(x(idx))
asum = asum - (real(ndm-1)/real(ndm))*abs(x(idx)) end do
end do
! compute global sum
call psb_sum(ictxt,asum)
else
asum=0.d0
! compute global sum
call psb_sum(ictxt, asum)
end if
else else
asum=0.d0 res = dzero
end if end if
! compute global sum
res = asum call psb_sum(ictxt,res)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -505,7 +456,3 @@ subroutine psb_dasumvs(res,x,desc_a, info)
end if end if
return return
end subroutine psb_dasumvs end subroutine psb_dasumvs

246
base/psblas/psb_samax.f90
Unescape Escape View File

@ -39,39 +39,34 @@
! where sub( X ) denotes X(1:N,JX:). ! where sub( X ) denotes X(1:N,JX:).
! !
! Arguments: ! Arguments:
! x(:,:) - real The input vector. ! x(:,:) - real The input vector.
! 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. ! jx - integer(optional). The column offset.
! !
function psb_samax (x,desc_a, info, jx) function psb_samax(x,desc_a, info, jx) result(res)
use psb_penv_mod use psb_base_mod, psb_protect_name => psb_samax
use psb_serial_mod
use psb_descriptor_type
use psb_check_mod
use psb_error_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
integer(psb_ipk_), optional, intent(in) :: jx integer(psb_ipk_), optional, intent(in) :: jx
real(psb_spk_) :: psb_samax real(psb_spk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ix, ijx, m, imax, isamax & err_act, iix, jjx, ix, ijx, m, ldx
real(psb_spk_) :: amax character(len=20) :: name, ch_err
character(len=20) :: name, ch_err
name='psb_samax' name='psb_samax'
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)
amax=0.d0
ictxt=desc_a%get_context() ictxt = desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if (np == -1) then if (np == -1) then
@ -82,39 +77,37 @@ function psb_samax (x,desc_a, info, jx)
ix = 1 ix = 1
if (present(jx)) then if (present(jx)) then
ijx = jx ijx = jx
else else
ijx = 1 ijx = 1
endif endif
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = size(x,1)
call psb_chkvect(m,1,size(x,1),ix,ijx,desc_a,info,iix,jjx) 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'
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
if (iix /= 1) then if (iix /= 1) then
info=psb_err_ix_n1_iy_n1_unsupported_ info=psb_err_ix_n1_iy_n1_unsupported_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
! compute local max ! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
imax=isamax(desc_a%get_local_rows()-iix+1,x(iix,jjx),1) res = psb_amax(desc_a%get_local_rows()-iix+1,x(:,jjx))
amax=abs(x(iix+imax-1,jjx))
else else
amax = szero res = szero
end if end if
! compute global max ! compute global max
call psb_amx(ictxt, amax) call psb_amx(ictxt, res)
psb_samax=amax
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -123,8 +116,8 @@ function psb_samax (x,desc_a, info, jx)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then if (err_act == psb_act_abort_) then
call psb_error(ictxt) call psb_error(ictxt)
return return
end if end if
return return
end function psb_samax end function psb_samax
@ -170,27 +163,24 @@ end function psb_samax
! normi := max(abs(X(i)) ! normi := max(abs(X(i))
! !
! Arguments: ! Arguments:
! x(:) - real The input vector. ! x(:) - real The input vector.
! 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_samaxv (x,desc_a, info) function psb_samaxv (x,desc_a, info) result(res)
use psb_penv_mod use psb_base_mod, psb_protect_name => psb_samaxv
use psb_serial_mod
use psb_descriptor_type
use psb_check_mod
use psb_error_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_samaxv real(psb_spk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, jx, ix, m, imax, isamax & err_act, iix, jjx, jx, ix, m, ldx
real(psb_spk_) :: amax
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_samaxv' name='psb_samaxv'
@ -198,7 +188,6 @@ function psb_samaxv (x,desc_a, info)
info=psb_success_ info=psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
amax=0.d0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
@ -213,33 +202,31 @@ function psb_samaxv (x,desc_a, info)
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,1),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'
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
if (iix /= 1) then if (iix /= 1) then
info=psb_err_ix_n1_iy_n1_unsupported_ info=psb_err_ix_n1_iy_n1_unsupported_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
! compute local max ! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
imax=isamax(desc_a%get_local_rows()-iix+1,x(iix),1) res = psb_amax(desc_a%get_local_rows()-iix+1,x)
amax=abs(x(iix+imax-1))
else else
amax = szero res = szero
end if end if
! compute global max ! compute global max
call psb_amx(ictxt, amax) call psb_amx(ictxt, res)
psb_samaxv=amax
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -248,8 +235,8 @@ function psb_samaxv (x,desc_a, info)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then if (err_act == psb_act_abort_) then
call psb_error(ictxt) call psb_error(ictxt)
return return
end if end if
return return
end function psb_samaxv end function psb_samaxv
@ -267,20 +254,18 @@ function psb_samax_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, jx, ix, m, imax, isamax & err_act, iix, jjx, jx, ix, m
real(psb_spk_) :: amax character(len=20) :: name, ch_err
character(len=20) :: name, ch_err
name='psb_samaxv' name='psb_samaxv'
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)
amax=0.d0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
@ -300,8 +285,7 @@ function psb_samax_vect(x, desc_a, info) result(res)
jx = 1 jx = 1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
call psb_chkvect(m,ione,x%get_nrows(),ix,jx,desc_a,info,iix,jjx)
call psb_chkvect(m,1,x%get_nrows(),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'
@ -317,15 +301,13 @@ function psb_samax_vect(x, desc_a, info) result(res)
! compute local max ! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
amax=x%amax(desc_a%get_local_rows()) res = x%amax(desc_a%get_local_rows())
else else
amax = szero res = szero
end if end if
! compute global max ! compute global max
call psb_amx(ictxt, amax) call psb_amx(ictxt, res)
res=amax
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -382,35 +364,34 @@ end function psb_samax_vect
! where sub( X ) denotes X(1:N,JX:). ! where sub( X ) denotes X(1:N,JX:).
! !
! Arguments: ! Arguments:
! res - real. The result. ! res - real The result.
! x(:,:) - real The input vector. ! x(:,:) - real The input vector.
! 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. ! jx - integer(optional). The column offset.
! !
subroutine psb_samaxvs (res,x,desc_a, info) subroutine psb_samaxvs(res,x,desc_a, info)
use psb_base_mod, psb_protect_name => psb_samaxvs use psb_base_mod, psb_protect_name => psb_samaxvs
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_), intent(out) :: res real(psb_spk_), intent(out) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ix, ijx, m, imax, isamax & err_act, iix, jjx, ix, ijx, m, ldx
real(psb_spk_) :: amax character(len=20) :: name, ch_err
character(len=20) :: name, ch_err
name='psb_samaxvs' name='psb_samaxvs'
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)
amax=0.d0
ictxt=desc_a%get_context() ictxt = desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if (np == -1) then if (np == -1) then
@ -423,33 +404,30 @@ subroutine psb_samaxvs (res,x,desc_a, info)
ijx=1 ijx=1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx=size(x,1)
call psb_chkvect(m,1,size(x,1),ix,ijx,desc_a,info,iix,jjx) 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'
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
if (iix /= 1) then if (iix /= 1) then
info=psb_err_ix_n1_iy_n1_unsupported_ info=psb_err_ix_n1_iy_n1_unsupported_
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
! compute local max ! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
imax=isamax(desc_a%get_local_rows()-iix+1,x(iix),1) res = psb_amax(desc_a%get_local_rows()-iix+1,x)
amax=abs(x(iix+imax-1))
else else
amax = szero res = szero
end if end if
! compute global max ! compute global max
call psb_amx(ictxt, amax) call psb_amx(ictxt, res)
res = amax
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -458,8 +436,8 @@ subroutine psb_samaxvs (res,x,desc_a, info)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then if (err_act == psb_act_abort_) then
call psb_error(ictxt) call psb_error(ictxt)
return return
end if end if
return return
end subroutine psb_samaxvs end subroutine psb_samaxvs
@ -503,25 +481,25 @@ end subroutine psb_samaxvs
! normi := max(abs(X(i)) ! normi := max(abs(X(i))
! !
! Arguments: ! Arguments:
! res(:) - real The result. ! res(:) - real. The result.
! x(:,:) - real The input vector. ! x(:,:) - real The input vector.
! 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_smamaxs (res,x,desc_a, info,jx) subroutine psb_smamaxs(res,x,desc_a, info,jx)
use psb_base_mod, psb_protect_name => psb_smamaxs use psb_base_mod, psb_protect_name => psb_smamaxs
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
integer(psb_ipk_), optional, intent(in) :: jx integer(psb_ipk_), optional, intent(in) :: jx
real(psb_spk_), intent(out) :: res(:) real(psb_spk_), intent(out) :: res(:)
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ix, ijx, m, imax, i, k, isamax & err_act, iix, jjx, ix, ijx, m, ldx, i, k
real(psb_spk_) :: amax
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_smamaxs' name='psb_smamaxs'
@ -529,8 +507,6 @@ subroutine psb_smamaxs (res,x,desc_a, info,jx)
info=psb_success_ info=psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
amax=0.d0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
@ -549,8 +525,8 @@ subroutine psb_smamaxs (res,x,desc_a, info,jx)
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
k = min(size(x,2),size(res,1)) k = min(size(x,2),size(res,1))
ldx = size(x,1)
call psb_chkvect(m,1,size(x,1),ix,ijx,desc_a,info,iix,jjx) 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'
@ -564,14 +540,12 @@ subroutine psb_smamaxs (res,x,desc_a, info,jx)
goto 9999 goto 9999
end if end if
res(1:k) = szero
! compute local max ! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
do i=1,k do i=1,k
imax=isamax(desc_a%get_local_rows()-iix+1,x(iix,jjx+i-1),1) res(i) = psb_amax(desc_a%get_local_rows()-iix+1,x(:,jjx+i-1))
res(i)=abs(x(iix+imax-1,jjx+i-1)) end do
end do
else
amax = szero
end if end if
! compute global max ! compute global max

347
base/psblas/psb_sasum.f90
Unescape Escape View File

@ -39,38 +39,32 @@
! where sub( X ) denotes X(1:N,JX:). ! where sub( X ) denotes X(1:N,JX:).
! !
! Arguments: ! Arguments:
! x(:,:) - real The input vector. ! x(:,:) - real The input vector.
! 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. ! jx - integer(optional). The column offset.
! !
function psb_sasum (x,desc_a, info, jx) function psb_sasum (x,desc_a, info, jx) result(res)
use psb_base_mod, psb_protect_name => psb_sasum
use psb_serial_mod
use psb_descriptor_type
use psb_check_mod
use psb_error_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
integer(psb_ipk_), optional, intent(in) :: jx integer(psb_ipk_), optional, intent(in) :: jx
real(psb_spk_) :: psb_sasum real(psb_spk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me, err_act, & integer(psb_ipk_) :: ictxt, np, me, &
& iix, jjx, ix, ijx, m, i, idx, ndm & err_act, iix, jjx, ix, ijx, m, i, idx, ndm, ldx
real(psb_spk_) :: asum, sasum character(len=20) :: name, ch_err
character(len=20) :: name, ch_err
name='psb_sasum' name='psb_sasum'
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)
asum=0.0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
@ -89,9 +83,9 @@ function psb_sasum (x,desc_a, info, jx)
endif endif
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = size(x,1)
! check vector correctness ! check vector correctness
call psb_chkvect(m,1,size(x,1),ix,ijx,desc_a,info,iix,jjx) 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'
@ -106,31 +100,21 @@ function psb_sasum (x,desc_a, info, jx)
end if end if
! compute local max ! compute local max
if ((m /= 0)) then if(desc_a%get_local_rows() > 0) then
if(desc_a%get_local_rows() > 0) then res = psb_asum(desc_a%get_local_rows()-iix+1,x(:,jjx))
asum=sasum(desc_a%get_local_rows()-iix+1,x(iix:,jjx),ione)
! adjust asum because overlapped elements are computed more than once
do i=1,size(desc_a%ovrlap_elem,1)
idx = desc_a%ovrlap_elem(i,1)
ndm = desc_a%ovrlap_elem(i,2)
asum = asum - (real(ndm-1)/real(ndm))*abs(x(idx,jjx))
end do
! compute global sum
call psb_sum(ictxt, asum)
else
asum=0.0
! compute global sum
call psb_sum(ictxt, asum)
end if
else
asum=0.0
end if
! adjust res because overlapped elements are computed more than once
do i=1,size(desc_a%ovrlap_elem,1)
idx = desc_a%ovrlap_elem(i,1)
ndm = desc_a%ovrlap_elem(i,2)
res = res - (real(ndm-1)/real(ndm))*psb_nrm1(x(idx,jjx))
end do
psb_sasum=asum else
res = szero
end if
! compute global sum
call psb_sum(ictxt, res)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -146,73 +130,25 @@ function psb_sasum (x,desc_a, info, jx)
end function psb_sasum end function psb_sasum
!!$ function psb_sasum_vect(x, desc_a, info) result(res)
!!$ Parallel Sparse BLAS version 3.0 use psb_base_mod, psb_protect_name => psb_sasum_vect
!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
!
! Function: psb_sasumv
! Computes norm1 of X
!
! norm1 := sum(X(i))
!
! Arguments:
! x(:) - real The input vector.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
!
function psb_sasumv (x,desc_a, info)
use psb_serial_mod
use psb_descriptor_type
use psb_check_mod
use psb_error_mod
use psb_penv_mod
implicit none implicit none
real(psb_spk_), intent(in) :: x(:) real(psb_spk_) :: res
type(psb_desc_type), intent(in) :: desc_a type(psb_s_vect_type), intent (inout) :: x
integer(psb_ipk_), intent(out) :: info type(psb_desc_type), intent (in) :: desc_a
real(psb_spk_) :: psb_sasumv integer(psb_ipk_), intent(out) :: info
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me, err_act, iix, jjx, jx, ix, m, i, idx, ndm integer(psb_ipk_) :: ictxt, np, me,&
real(psb_spk_) :: asum, sasum & err_act, iix, jjx, jx, ix, m, imax
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_sasumv' name='psb_sasumv'
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)
asum=0.0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
@ -223,13 +159,18 @@ function psb_sasumv (x,desc_a, info)
goto 9999 goto 9999
endif endif
if (.not.allocated(x%v)) then
info = psb_err_invalid_vect_state_
call psb_errpush(info,name)
goto 9999
endif
ix = 1 ix = 1
jx=1 jx = 1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
call psb_chkvect(m,ione,x%get_nrows(),ix,jx,desc_a,info,iix,jjx)
! check vector correctness
call psb_chkvect(m,1,size(x),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'
@ -244,29 +185,14 @@ function psb_sasumv (x,desc_a, info)
end if end if
! compute local max ! compute local max
if ((m /= 0)) then if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
if(desc_a%get_local_rows() > 0) then res = x%asum(desc_a%get_local_rows())
asum=sasum(desc_a%get_local_rows(),x,ione)
! adjust asum because overlapped elements are computed more than once
do i=1,size(desc_a%ovrlap_elem,1)
idx = desc_a%ovrlap_elem(i,1)
ndm = desc_a%ovrlap_elem(i,2)
asum = asum - (real(ndm-1)/real(ndm))*abs(x(idx))
end do
! compute global sum
call psb_sum(ictxt, asum)
else
asum=0.0
! compute global sum
call psb_sum(ictxt, asum)
end if
else else
asum=0.0 res = szero
end if end if
psb_sasumv=asum ! compute global sum
call psb_sum(ictxt, res)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -279,36 +205,73 @@ function psb_sasumv (x,desc_a, info)
return return
end if end if
return return
end function psb_sasumv
end function psb_sasum_vect
!!$
!!$ Parallel Sparse BLAS version 3.0
!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
!
! Function: psb_sasumv
! Computes norm1 of X
!
! norm1 := sum(X(i))
!
! Arguments:
! x(:) - real The input vector.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
!
function psb_sasumv(x,desc_a, info) result(res)
use psb_base_mod, psb_protect_name => psb_sasumv
function psb_sasum_vect(x, desc_a, info) result(res)
use psb_penv_mod
use psb_serial_mod
use psb_descriptor_type
use psb_check_mod
use psb_error_mod
use psb_s_vect_mod
implicit none implicit none
real(psb_spk_) :: res real(psb_spk_), intent(in) :: 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 real(psb_spk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, jx, ix, m, imax & err_act, iix, jjx, jx, ix, m, i, idx, ndm, ldx
real(psb_spk_) :: asum
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_sasumv' name='psb_sasumv'
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)
asum=0.d0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
@ -318,19 +281,13 @@ function psb_sasum_vect(x, desc_a, info) result(res)
goto 9999 goto 9999
endif endif
if (.not.allocated(x%v)) then
info = psb_err_invalid_vect_state_
call psb_errpush(info,name)
goto 9999
endif
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,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) ! check vector correctness
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'
@ -345,16 +302,22 @@ function psb_sasum_vect(x, desc_a, info) result(res)
end if end if
! compute local max ! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then if(desc_a%get_local_rows() > 0) then
asum=x%asum(desc_a%get_local_rows()) res = psb_asum(desc_a%get_local_rows(),x)
! adjust asum because overlapped elements are computed more than once
do i=1,size(desc_a%ovrlap_elem,1)
idx = desc_a%ovrlap_elem(i,1)
ndm = desc_a%ovrlap_elem(i,2)
res = res - (real(ndm-1)/real(ndm))*psb_nrm1(x(idx))
end do
else else
asum = szero res = szero
end if end if
! compute global sum ! compute global sum
call psb_sum(ictxt, asum) call psb_sum(ictxt, res)
res=asum
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -367,9 +330,7 @@ function psb_sasum_vect(x, desc_a, info) result(res)
return return
end if end if
return return
end function psb_sasumv
end function psb_sasum_vect
!!$ !!$
@ -410,33 +371,32 @@ end function psb_sasum_vect
! norm1 := sum(X(i)) ! norm1 := sum(X(i))
! !
! Arguments: ! Arguments:
! res - real The result. ! res - real. The result.
! x(:) - real The input vector. ! x(:) - real The input vector.
! 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. ! jx - integer(optional). The column offset.
! !
subroutine psb_sasumvs(res,x,desc_a, info) subroutine psb_sasumvs(res,x,desc_a, info)
use psb_base_mod, psb_protect_name => psb_sasumvs use psb_base_mod, psb_protect_name => psb_sasumvs
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, err_act, iix, jjx, ix, jx, m, i, idx, ndm integer(psb_ipk_) :: ictxt, np, me,&
real(psb_spk_) :: asum, sasum & err_act, iix, jjx, ix, jx, m, i, idx, ndm, ldx
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_sasumvs' name='psb_sasumvs'
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)
asum=0.0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
@ -450,9 +410,9 @@ subroutine psb_sasumvs(res,x,desc_a, info)
jx = 1 jx = 1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = size(x,1)
! check vector correctness ! check vector correctness
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'
@ -467,31 +427,22 @@ subroutine psb_sasumvs(res,x,desc_a, info)
end if end if
! compute local max ! compute local max
if ((m /= 0)) then if(desc_a%get_local_rows() > 0) then
if(desc_a%get_local_rows() > 0) then res = psb_asum(desc_a%get_local_rows(),x)
asum=sasum(desc_a%get_local_rows(),x,ione)
! adjust asum because overlapped elements are computed more than once
! adjust asum 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) res = res - (real(ndm-1)/real(ndm))*psb_nrm1(x(idx))
asum = asum - (real(ndm-1)/real(ndm))*abs(x(idx)) end do
end do
! compute global sum
call psb_sum(ictxt,asum)
else
asum=0.0
! compute global sum
call psb_sum(ictxt, asum)
end if
else else
asum=0.0 res = szero
end if end if
! compute global sum
res = asum call psb_sum(ictxt,res)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

312
base/psblas/psb_zamax.f90
Unescape Escape View File

@ -44,37 +44,29 @@
! info - integer. Return code ! info - integer. Return code
! jx - integer(optional). The column offset. ! jx - integer(optional). The column offset.
! !
function psb_zamax (x,desc_a, info, jx) function psb_zamax(x,desc_a, info, jx) result(res)
use psb_penv_mod use psb_base_mod, psb_protect_name => psb_zamax
use psb_serial_mod
use psb_descriptor_type
use psb_check_mod
use psb_error_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
integer(psb_ipk_), optional, intent(in) :: jx integer(psb_ipk_), optional, intent(in) :: jx
real(psb_dpk_) :: psb_zamax real(psb_dpk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ix, ijx, m, imax, izamax & err_act, iix, jjx, ix, ijx, m, ldx
real(psb_dpk_) :: amax character(len=20) :: name, ch_err
character(len=20) :: name, ch_err
double complex :: zdum
double precision :: cabs1
cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
name='psb_zamax' name='psb_zamax'
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)
amax=0.d0
ictxt=desc_a%get_context() ictxt = desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if (np == -1) then if (np == -1) then
@ -91,8 +83,9 @@ function psb_zamax (x,desc_a, info, jx)
endif endif
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = size(x,1)
call psb_chkvect(m,1,size(x,1),ix,ijx,desc_a,info,iix,jjx) 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'
@ -108,16 +101,13 @@ function psb_zamax (x,desc_a, info, jx)
! compute local max ! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
imax=izamax(desc_a%get_local_rows()-iix+1,x(iix,jjx),1) res = psb_amax(desc_a%get_local_rows()-iix+1,x(:,jjx))
amax=cabs1(x(iix+imax-1,jjx))
else else
amax = dzero res = dzero
end if end if
! compute global max ! compute global max
call psb_amx(ictxt, amax) call psb_amx(ictxt, res)
psb_zamax=amax
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -133,93 +123,6 @@ function psb_zamax (x,desc_a, info, jx)
end function psb_zamax end function psb_zamax
function psb_zamax_vect(x, desc_a, info) result(res)
use psb_penv_mod
use psb_serial_mod
use psb_descriptor_type
use psb_check_mod
use psb_error_mod
use psb_z_vect_mod
implicit none
real(psb_dpk_) :: res
type(psb_z_vect_type), intent (inout) :: x
type(psb_desc_type), intent (in) :: desc_a
integer(psb_ipk_), intent(out) :: info
! locals
integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, jx, ix, m, imax, isamax
real(psb_dpk_) :: amax
character(len=20) :: name, ch_err
name='psb_zamaxv'
if(psb_get_errstatus() /= 0) return
info=psb_success_
call psb_erractionsave(err_act)
amax=dzero
ictxt=desc_a%get_context()
call psb_info(ictxt, me, np)
if (np == -1) then
info = psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
if (.not.allocated(x%v)) then
info = psb_err_invalid_vect_state_
call psb_errpush(info,name)
goto 9999
endif
ix = 1
jx = 1
m = desc_a%get_global_rows()
call psb_chkvect(m,1,x%get_nrows(),ix,jx,desc_a,info,iix,jjx)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_chkvect'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (iix /= 1) then
info=psb_err_ix_n1_iy_n1_unsupported_
call psb_errpush(info,name)
goto 9999
end if
! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
amax=x%amax(desc_a%get_local_rows())
else
amax = dzero
end if
! compute global max
call psb_amx(ictxt, amax)
res=amax
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error(ictxt)
return
end if
return
end function psb_zamax_vect
!!$ !!$
@ -264,35 +167,27 @@ end function psb_zamax_vect
! 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_zamaxv (x,desc_a, info) function psb_zamaxv (x,desc_a, info) result(res)
use psb_penv_mod use psb_base_mod, psb_protect_name => psb_zamaxv
use psb_serial_mod
use psb_descriptor_type
use psb_check_mod
use psb_error_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_zamaxv real(psb_dpk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, jx, ix, m, imax, izamax & err_act, iix, jjx, jx, ix, m, ldx
real(psb_dpk_) :: amax
complex(psb_dpk_) :: cmax
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
double complex :: zdum
double precision :: cabs1
cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
name='psb_zamaxv' name='psb_zamaxv'
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)
amax=0.d0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
@ -307,8 +202,9 @@ function psb_zamaxv (x,desc_a, info)
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,1),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'
@ -324,17 +220,13 @@ function psb_zamaxv (x,desc_a, info)
! compute local max ! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
imax=izamax(desc_a%get_local_rows()-iix+1,x(iix),1) res = psb_amax(desc_a%get_local_rows()-iix+1,x)
cmax=(x(iix+imax-1))
amax=cabs1(cmax)
else else
amax = dzero res = dzero
end if end if
! compute global max ! compute global max
call psb_amx(ictxt, amax) call psb_amx(ictxt, res)
psb_zamaxv=amax
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -349,6 +241,89 @@ function psb_zamaxv (x,desc_a, info)
return return
end function psb_zamaxv end function psb_zamaxv
function psb_zamax_vect(x, desc_a, info) result(res)
use psb_penv_mod
use psb_serial_mod
use psb_descriptor_type
use psb_check_mod
use psb_error_mod
use psb_z_vect_mod
implicit none
real(psb_dpk_) :: res
type(psb_z_vect_type), intent (inout) :: x
type(psb_desc_type), intent (in) :: desc_a
integer(psb_ipk_), intent(out) :: info
! locals
integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, jx, ix, m
character(len=20) :: name, ch_err
name='psb_zamaxv'
if(psb_get_errstatus() /= 0) return
info=psb_success_
call psb_erractionsave(err_act)
ictxt=desc_a%get_context()
call psb_info(ictxt, me, np)
if (np == -1) then
info = psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
if (.not.allocated(x%v)) then
info = psb_err_invalid_vect_state_
call psb_errpush(info,name)
goto 9999
endif
ix = 1
jx = 1
m = desc_a%get_global_rows()
call psb_chkvect(m,ione,x%get_nrows(),ix,jx,desc_a,info,iix,jjx)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_chkvect'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (iix /= 1) then
info=psb_err_ix_n1_iy_n1_unsupported_
call psb_errpush(info,name)
goto 9999
end if
! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
res = x%amax(desc_a%get_local_rows())
else
res = dzero
end if
! compute global max
call psb_amx(ictxt, res)
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error(ictxt)
return
end if
return
end function psb_zamax_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
@ -397,31 +372,26 @@ end function psb_zamaxv
! !
subroutine psb_zamaxvs(res,x,desc_a, info) subroutine psb_zamaxvs(res,x,desc_a, info)
use psb_base_mod, psb_protect_name => psb_zamaxvs use psb_base_mod, psb_protect_name => psb_zamaxvs
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_), intent(out) :: res real(psb_dpk_), intent(out) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ix, ijx, m, imax, izamax & err_act, iix, jjx, ix, ijx, m, ldx
real(psb_dpk_) :: amax character(len=20) :: name, ch_err
character(len=20) :: name, ch_err
complex(psb_dpk_) :: cmax
double complex :: zdum
double precision :: cabs1
cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
name='psb_zamaxvs' name='psb_zamaxvs'
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)
amax=0.d0
ictxt=desc_a%get_context() ictxt = desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if (np == -1) then if (np == -1) then
@ -434,7 +404,8 @@ subroutine psb_zamaxvs(res,x,desc_a, info)
ijx=1 ijx=1
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'
@ -450,17 +421,13 @@ subroutine psb_zamaxvs(res,x,desc_a, info)
! compute local max ! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
imax=izamax(desc_a%get_local_rows()-iix+1,x(iix),1) res = psb_amax(desc_a%get_local_rows()-iix+1,x)
cmax=(x(iix+imax-1))
amax=cabs1(cmax)
else else
amax = dzero res = dzero
end if end if
! compute global max ! compute global max
call psb_amx(ictxt, amax) call psb_amx(ictxt, res)
res = amax
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -521,31 +488,25 @@ end subroutine psb_zamaxvs
! !
subroutine psb_zmamaxs(res,x,desc_a, info,jx) subroutine psb_zmamaxs(res,x,desc_a, info,jx)
use psb_base_mod, psb_protect_name => psb_zmamaxs use psb_base_mod, psb_protect_name => psb_zmamaxs
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
integer(psb_ipk_), optional, intent(in) :: jx integer(psb_ipk_), optional, intent(in) :: jx
real(psb_dpk_), intent(out) :: res(:) real(psb_dpk_), intent(out) :: res(:)
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, ix, ijx, m, imax, i, k, izamax & err_act, iix, jjx, ix, ijx, m, ldx, i, k
real(psb_dpk_) :: amax
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
complex(psb_dpk_) :: cmax
double complex :: zdum
double precision :: cabs1
cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
name='psb_zmamaxs' name='psb_zmamaxs'
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)
amax=0.d0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
@ -564,8 +525,8 @@ subroutine psb_zmamaxs(res,x,desc_a, info,jx)
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
k = min(size(x,2),size(res,1)) k = min(size(x,2),size(res,1))
ldx = size(x,1)
call psb_chkvect(m,1,size(x,1),ix,ijx,desc_a,info,iix,jjx) 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'
@ -579,15 +540,12 @@ subroutine psb_zmamaxs(res,x,desc_a, info,jx)
goto 9999 goto 9999
end if end if
res(1:k) = dzero
! compute local max ! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
do i=1,k do i=1,k
imax=izamax(desc_a%get_local_rows()-iix+1,x(iix,jjx+i-1),1) res(i) = psb_amax(desc_a%get_local_rows()-iix+1,x(:,jjx+i-1))
cmax=(x(iix+imax-1,jjx+i-1)) end do
res(i)=cabs1(cmax)
end do
else
amax = dzero
end if end if
! compute global max ! compute global max

200
base/psblas/psb_zasum.f90
Unescape Escape View File

@ -44,37 +44,27 @@
! info - integer. Return code ! info - integer. Return code
! jx - integer(optional). The column offset. ! jx - integer(optional). The column offset.
! !
function psb_zasum (x,desc_a, info, jx) function psb_zasum (x,desc_a, info, jx) result(res)
use psb_base_mod, psb_protect_name => psb_zasum
use psb_serial_mod
use psb_descriptor_type
use psb_check_mod
use psb_error_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
integer(psb_ipk_), optional, intent(in) :: jx integer(psb_ipk_), optional, intent(in) :: jx
real(psb_dpk_) :: psb_zasum real(psb_dpk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me, & integer(psb_ipk_) :: ictxt, np, me, &
& err_act, iix, jjx, ix, ijx, m, i, idx, ndm & err_act, iix, jjx, ix, ijx, m, i, idx, ndm, ldx
real(psb_dpk_) :: asum, dzasum
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
complex(psb_dpk_) :: cmax
double complex :: zdum
double precision :: cabs1
cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
name='psb_zasum' name='psb_zasum'
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)
asum=0.d0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
@ -93,9 +83,9 @@ function psb_zasum (x,desc_a, info, jx)
endif endif
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = size(x,1)
! check vector correctness ! check vector correctness
call psb_chkvect(m,1,size(x,1),ix,ijx,desc_a,info,iix,jjx) 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'
@ -110,29 +100,21 @@ function psb_zasum (x,desc_a, info, jx)
end if end if
! compute local max ! compute local max
if ((m /= 0)) then if(desc_a%get_local_rows() > 0) then
if(desc_a%get_local_rows() > 0) then res = psb_asum(desc_a%get_local_rows()-iix+1,x(:,jjx))
asum=dzasum(desc_a%get_local_rows()-iix+1,x(iix:,jjx),ione)
! adjust asum because overlapped elements are computed more than once
do i=1,size(desc_a%ovrlap_elem,1)
idx = desc_a%ovrlap_elem(i,1)
ndm = desc_a%ovrlap_elem(i,2)
asum = asum - (real(ndm-1)/real(ndm))*cabs1(x(idx,jjx))
end do
else
asum=0.d0
end if
! compute global sum
call psb_sum(ictxt, asum)
else
asum=0.d0
end if
! adjust res because overlapped elements are computed more than once
do i=1,size(desc_a%ovrlap_elem,1)
idx = desc_a%ovrlap_elem(i,1)
ndm = desc_a%ovrlap_elem(i,2)
res = res - (real(ndm-1)/real(ndm))*psb_nrm1(x(idx,jjx))
end do
psb_zasum=asum else
res = dzero
end if
! compute global sum
call psb_sum(ictxt, res)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -141,31 +123,25 @@ function psb_zasum (x,desc_a, info, jx)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then if (err_act == psb_act_abort_) then
call psb_error(ictxt) call psb_error(ictxt)
return return
end if end if
return return
end function psb_zasum end function psb_zasum
function psb_zasum_vect(x, desc_a, info) result(res) function psb_zasum_vect(x, desc_a, info) result(res)
use psb_penv_mod use psb_base_mod, psb_protect_name => psb_zasum_vect
use psb_serial_mod
use psb_descriptor_type
use psb_check_mod
use psb_error_mod
use psb_z_vect_mod
implicit none implicit none
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, jx, ix, m, imax & err_act, iix, jjx, jx, ix, m, imax
real(psb_dpk_) :: asum
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_zasumv' name='psb_zasumv'
@ -173,7 +149,6 @@ function psb_zasum_vect(x, desc_a, info) result(res)
info=psb_success_ info=psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
asum=0.d0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
@ -195,8 +170,7 @@ function psb_zasum_vect(x, desc_a, info) result(res)
jx = 1 jx = 1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
call psb_chkvect(m,ione,x%get_nrows(),ix,jx,desc_a,info,iix,jjx)
call psb_chkvect(m,1,x%get_nrows(),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'
@ -212,15 +186,13 @@ function psb_zasum_vect(x, desc_a, info) result(res)
! compute local max ! compute local max
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
asum=x%asum(desc_a%get_local_rows()) res = x%asum(desc_a%get_local_rows())
else else
asum = dzero res = dzero
end if end if
! compute global sum ! compute global sum
call psb_sum(ictxt, asum) call psb_sum(ictxt, res)
res=asum
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -280,37 +252,26 @@ end function psb_zasum_vect
! 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_zasumv(x,desc_a, info) function psb_zasumv(x,desc_a, info) result(res)
use psb_base_mod, psb_protect_name => psb_zasumv
use psb_serial_mod
use psb_descriptor_type
use psb_check_mod
use psb_error_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_zasumv real(psb_dpk_) :: res
! locals ! locals
integer(psb_ipk_) :: ictxt, np, me,& integer(psb_ipk_) :: ictxt, np, me,&
& err_act, iix, jjx, jx, ix, m, i, idx, ndm & err_act, iix, jjx, jx, ix, m, i, idx, ndm, ldx
real(psb_dpk_) :: asum, dzasum
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
complex(psb_dpk_) :: cmax
double complex :: zdum
double precision :: cabs1
cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
name='psb_zasumv' name='psb_zasumv'
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)
asum=0.d0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
@ -324,9 +285,9 @@ function psb_zasumv(x,desc_a, info)
jx=1 jx=1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = size(x,1)
! check vector correctness ! check vector correctness
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'
@ -341,30 +302,22 @@ function psb_zasumv(x,desc_a, info)
end if end if
! compute local max ! compute local max
if ((m /= 0)) then if(desc_a%get_local_rows() > 0) then
if(desc_a%get_local_rows() > 0) then res = psb_asum(desc_a%get_local_rows(),x)
asum=dzasum(desc_a%get_local_rows(),x,ione)
! adjust asum because overlapped elements are computed more than once
! adjust asum 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) res = res - (real(ndm-1)/real(ndm))*psb_nrm1(x(idx))
asum = asum - (real(ndm-1)/real(ndm))*cabs1(x(idx)) end do
end do
! compute global sum
call psb_sum(ictxt, asum)
else
asum=0.d0
! compute global sum
call psb_sum(ictxt, asum)
end if
else else
asum=0.d0 res = dzero
end if end if
psb_zasumv=asum ! compute global sum
call psb_sum(ictxt, res)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -426,30 +379,24 @@ end function psb_zasumv
! !
subroutine psb_zasumvs(res,x,desc_a, info) subroutine psb_zasumvs(res,x,desc_a, info)
use psb_base_mod, psb_protect_name => psb_zasumvs use psb_base_mod, psb_protect_name => psb_zasumvs
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, ix, jx, m, i, idx, ndm & err_act, iix, jjx, ix, jx, m, i, idx, ndm, ldx
real(psb_dpk_) :: asum, dzasum
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
complex(psb_dpk_) :: cmax
double complex :: zdum
double precision :: cabs1
cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
name='psb_zasumvs' name='psb_zasumvs'
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)
asum=0.d0
ictxt=desc_a%get_context() ictxt=desc_a%get_context()
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
@ -463,9 +410,9 @@ subroutine psb_zasumvs(res,x,desc_a, info)
jx = 1 jx = 1
m = desc_a%get_global_rows() m = desc_a%get_global_rows()
ldx = size(x,1)
! check vector correctness ! check vector correctness
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'
@ -480,31 +427,22 @@ subroutine psb_zasumvs(res,x,desc_a, info)
end if end if
! compute local max ! compute local max
if ((m /= 0)) then if(desc_a%get_local_rows() > 0) then
if(desc_a%get_local_rows() > 0) then res = psb_asum(desc_a%get_local_rows(),x)
asum=dzasum(desc_a%get_local_rows(),x,ione)
! adjust asum because overlapped elements are computed more than once
! adjust asum 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) res = res - (real(ndm-1)/real(ndm))*psb_nrm1(x(idx))
asum = asum - (real(ndm-1)/real(ndm))*cabs1(x(idx)) end do
end do
! compute global sum
call psb_sum(ictxt,asum)
else
asum=0.d0
! compute global sum
call psb_sum(ictxt, asum)
end if
else else
asum=0.d0 res = dzero
end if end if
! compute global sum
res = asum call psb_sum(ictxt,res)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

4
base/serial/Makefile
Unescape Escape View File

@ -7,7 +7,9 @@ FOBJS = psb_lsame.o psi_serial_impl.o psb_sort_impl.o \
psb_snumbmm.o psb_dnumbmm.o psb_cnumbmm.o psb_znumbmm.o \ psb_snumbmm.o psb_dnumbmm.o psb_cnumbmm.o psb_znumbmm.o \
psb_sgeprt.o psb_dgeprt.o psb_cgeprt.o psb_zgeprt.o\ psb_sgeprt.o psb_dgeprt.o psb_cgeprt.o psb_zgeprt.o\
psb_spdot_srtd.o psb_aspxpby.o psb_spge_dot.o\ psb_spdot_srtd.o psb_aspxpby.o psb_spge_dot.o\
psb_sgelp.o psb_dgelp.o psb_cgelp.o psb_zgelp.o psb_sgelp.o psb_dgelp.o psb_cgelp.o psb_zgelp.o \
psb_samax_s.o psb_damax_s.o psb_camax_s.o psb_zamax_s.o \
psb_sasum_s.o psb_dasum_s.o psb_casum_s.o psb_zasum_s.o
LIBDIR=.. LIBDIR=..
MODDIR=../modules MODDIR=../modules

Write Preview
Loading…
Cancel
Save