diff --git a/base/modules/psb_c_base_mat_mod.f90 b/base/modules/psb_c_base_mat_mod.f90 index 64f77e87..2334ef30 100644 --- a/base/modules/psb_c_base_mat_mod.f90 +++ b/base/modules/psb_c_base_mat_mod.f90 @@ -1170,9 +1170,9 @@ contains if (allocated(a%ja)) deallocate(a%ja) if (allocated(a%val)) deallocate(a%val) call a%set_null() - call a%set_nrows(0) - call a%set_ncols(0) - call a%set_nzeros(0) + call a%set_nrows(izero) + call a%set_ncols(izero) + call a%set_nzeros(izero) return diff --git a/base/modules/psb_c_base_vect_mod.f90 b/base/modules/psb_c_base_vect_mod.f90 index 384d49b3..34aa4007 100644 --- a/base/modules/psb_c_base_vect_mod.f90 +++ b/base/modules/psb_c_base_vect_mod.f90 @@ -150,7 +150,7 @@ contains integer(psb_ipk_) :: info this%v = x - call this%asb(size(x),info) + call this%asb(size(x,kind=psb_ipk_),info) end function constructor diff --git a/base/modules/psb_c_csc_mat_mod.f90 b/base/modules/psb_c_csc_mat_mod.f90 index 60034d64..e23a73ec 100644 --- a/base/modules/psb_c_csc_mat_mod.f90 +++ b/base/modules/psb_c_csc_mat_mod.f90 @@ -527,8 +527,8 @@ contains if (allocated(a%ia)) deallocate(a%ia) if (allocated(a%val)) deallocate(a%val) call a%set_null() - call a%set_nrows(0) - call a%set_ncols(0) + call a%set_nrows(izero) + call a%set_ncols(izero) return diff --git a/base/modules/psb_c_csr_mat_mod.f90 b/base/modules/psb_c_csr_mat_mod.f90 index fb990144..88c7930f 100644 --- a/base/modules/psb_c_csr_mat_mod.f90 +++ b/base/modules/psb_c_csr_mat_mod.f90 @@ -527,8 +527,8 @@ contains if (allocated(a%ja)) deallocate(a%ja) if (allocated(a%val)) deallocate(a%val) call a%set_null() - call a%set_nrows(0) - call a%set_ncols(0) + call a%set_nrows(izero) + call a%set_ncols(izero) return diff --git a/base/modules/psb_c_linmap_mod.f90 b/base/modules/psb_c_linmap_mod.f90 index 7f49477d..d3efe916 100644 --- a/base/modules/psb_c_linmap_mod.f90 +++ b/base/modules/psb_c_linmap_mod.f90 @@ -55,7 +55,7 @@ module psb_c_linmap_mod interface psb_map_X2Y subroutine psb_c_map_X2Y(alpha,x,beta,y,map,info,work) use psb_const_mod - import :: psb_ipk_, psb_clinmap_type + import :: psb_clinmap_type implicit none type(psb_clinmap_type), intent(in) :: map 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) use psb_const_mod use psb_c_vect_mod - import :: psb_ipk_, psb_clinmap_type + import :: psb_clinmap_type implicit none type(psb_clinmap_type), intent(in) :: map complex(psb_spk_), intent(in) :: alpha,beta @@ -80,7 +80,7 @@ module psb_c_linmap_mod interface psb_map_Y2X subroutine psb_c_map_Y2X(alpha,x,beta,y,map,info,work) use psb_const_mod - import :: psb_ipk_, psb_clinmap_type + import :: psb_clinmap_type implicit none type(psb_clinmap_type), intent(in) :: map 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) use psb_const_mod use psb_c_vect_mod - import :: psb_ipk_, psb_clinmap_type + import :: psb_clinmap_type implicit none type(psb_clinmap_type), intent(in) :: map complex(psb_spk_), intent(in) :: alpha,beta diff --git a/base/modules/psb_c_vect_mod.f90 b/base/modules/psb_c_vect_mod.f90 index 29d809bf..c910cd16 100644 --- a/base/modules/psb_c_vect_mod.f90 +++ b/base/modules/psb_c_vect_mod.f90 @@ -163,7 +163,7 @@ contains 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 diff --git a/base/modules/psb_const_mod.F90 b/base/modules/psb_const_mod.F90 index 3fc12f94..c71da304 100644 --- a/base/modules/psb_const_mod.F90 +++ b/base/modules/psb_const_mod.F90 @@ -44,6 +44,9 @@ module psb_const_mod ! This is an 8-byte integer, and normally different from default integer. integer, parameter :: longndig=12 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 ! and MPI_REAL diff --git a/base/modules/psb_d_base_mat_mod.f90 b/base/modules/psb_d_base_mat_mod.f90 index ea682074..8e22c7e9 100644 --- a/base/modules/psb_d_base_mat_mod.f90 +++ b/base/modules/psb_d_base_mat_mod.f90 @@ -1170,9 +1170,9 @@ contains if (allocated(a%ja)) deallocate(a%ja) if (allocated(a%val)) deallocate(a%val) call a%set_null() - call a%set_nrows(0) - call a%set_ncols(0) - call a%set_nzeros(0) + call a%set_nrows(izero) + call a%set_ncols(izero) + call a%set_nzeros(izero) return diff --git a/base/modules/psb_d_base_vect_mod.f90 b/base/modules/psb_d_base_vect_mod.f90 index d7c00e64..3ee19d4d 100644 --- a/base/modules/psb_d_base_vect_mod.f90 +++ b/base/modules/psb_d_base_vect_mod.f90 @@ -150,7 +150,7 @@ contains integer(psb_ipk_) :: info this%v = x - call this%asb(size(x),info) + call this%asb(size(x,kind=psb_ipk_),info) end function constructor diff --git a/base/modules/psb_d_csc_mat_mod.f90 b/base/modules/psb_d_csc_mat_mod.f90 index d99edf39..cbacf911 100644 --- a/base/modules/psb_d_csc_mat_mod.f90 +++ b/base/modules/psb_d_csc_mat_mod.f90 @@ -527,8 +527,8 @@ contains if (allocated(a%ia)) deallocate(a%ia) if (allocated(a%val)) deallocate(a%val) call a%set_null() - call a%set_nrows(0) - call a%set_ncols(0) + call a%set_nrows(izero) + call a%set_ncols(izero) return diff --git a/base/modules/psb_d_csr_mat_mod.f90 b/base/modules/psb_d_csr_mat_mod.f90 index e698f315..5633ac8d 100644 --- a/base/modules/psb_d_csr_mat_mod.f90 +++ b/base/modules/psb_d_csr_mat_mod.f90 @@ -527,8 +527,8 @@ contains if (allocated(a%ja)) deallocate(a%ja) if (allocated(a%val)) deallocate(a%val) call a%set_null() - call a%set_nrows(0) - call a%set_ncols(0) + call a%set_nrows(izero) + call a%set_ncols(izero) return diff --git a/base/modules/psb_d_linmap_mod.f90 b/base/modules/psb_d_linmap_mod.f90 index 15273488..ac8598fd 100644 --- a/base/modules/psb_d_linmap_mod.f90 +++ b/base/modules/psb_d_linmap_mod.f90 @@ -55,7 +55,7 @@ module psb_d_linmap_mod interface psb_map_X2Y subroutine psb_d_map_X2Y(alpha,x,beta,y,map,info,work) use psb_const_mod - import :: psb_ipk_, psb_dlinmap_type + import :: psb_dlinmap_type implicit none type(psb_dlinmap_type), intent(in) :: map 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) use psb_const_mod use psb_d_vect_mod - import :: psb_ipk_, psb_dlinmap_type + import :: psb_dlinmap_type implicit none type(psb_dlinmap_type), intent(in) :: map real(psb_dpk_), intent(in) :: alpha,beta @@ -80,7 +80,7 @@ module psb_d_linmap_mod interface psb_map_Y2X subroutine psb_d_map_Y2X(alpha,x,beta,y,map,info,work) use psb_const_mod - import :: psb_ipk_, psb_dlinmap_type + import :: psb_dlinmap_type implicit none type(psb_dlinmap_type), intent(in) :: map 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) use psb_const_mod use psb_d_vect_mod - import :: psb_ipk_, psb_dlinmap_type + import :: psb_dlinmap_type implicit none type(psb_dlinmap_type), intent(in) :: map real(psb_dpk_), intent(in) :: alpha,beta diff --git a/base/modules/psb_d_vect_mod.f90 b/base/modules/psb_d_vect_mod.f90 index d3e10562..98f60e66 100644 --- a/base/modules/psb_d_vect_mod.f90 +++ b/base/modules/psb_d_vect_mod.f90 @@ -163,7 +163,7 @@ contains 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 diff --git a/base/modules/psb_s_base_mat_mod.f90 b/base/modules/psb_s_base_mat_mod.f90 index 8898526b..f5c0a757 100644 --- a/base/modules/psb_s_base_mat_mod.f90 +++ b/base/modules/psb_s_base_mat_mod.f90 @@ -1170,9 +1170,9 @@ contains if (allocated(a%ja)) deallocate(a%ja) if (allocated(a%val)) deallocate(a%val) call a%set_null() - call a%set_nrows(0) - call a%set_ncols(0) - call a%set_nzeros(0) + call a%set_nrows(izero) + call a%set_ncols(izero) + call a%set_nzeros(izero) return diff --git a/base/modules/psb_s_base_vect_mod.f90 b/base/modules/psb_s_base_vect_mod.f90 index ab337d9f..a3f38bd9 100644 --- a/base/modules/psb_s_base_vect_mod.f90 +++ b/base/modules/psb_s_base_vect_mod.f90 @@ -150,7 +150,7 @@ contains integer(psb_ipk_) :: info this%v = x - call this%asb(size(x),info) + call this%asb(size(x,kind=psb_ipk_),info) end function constructor diff --git a/base/modules/psb_s_csc_mat_mod.f90 b/base/modules/psb_s_csc_mat_mod.f90 index 8c5a2e30..99fe2600 100644 --- a/base/modules/psb_s_csc_mat_mod.f90 +++ b/base/modules/psb_s_csc_mat_mod.f90 @@ -527,8 +527,8 @@ contains if (allocated(a%ia)) deallocate(a%ia) if (allocated(a%val)) deallocate(a%val) call a%set_null() - call a%set_nrows(0) - call a%set_ncols(0) + call a%set_nrows(izero) + call a%set_ncols(izero) return diff --git a/base/modules/psb_s_csr_mat_mod.f90 b/base/modules/psb_s_csr_mat_mod.f90 index f302d612..0ed5d338 100644 --- a/base/modules/psb_s_csr_mat_mod.f90 +++ b/base/modules/psb_s_csr_mat_mod.f90 @@ -527,8 +527,8 @@ contains if (allocated(a%ja)) deallocate(a%ja) if (allocated(a%val)) deallocate(a%val) call a%set_null() - call a%set_nrows(0) - call a%set_ncols(0) + call a%set_nrows(izero) + call a%set_ncols(izero) return diff --git a/base/modules/psb_s_linmap_mod.f90 b/base/modules/psb_s_linmap_mod.f90 index 86b2ad70..577c3fc6 100644 --- a/base/modules/psb_s_linmap_mod.f90 +++ b/base/modules/psb_s_linmap_mod.f90 @@ -55,7 +55,7 @@ module psb_s_linmap_mod interface psb_map_X2Y subroutine psb_s_map_X2Y(alpha,x,beta,y,map,info,work) use psb_const_mod - import :: psb_ipk_, psb_slinmap_type + import :: psb_slinmap_type implicit none type(psb_slinmap_type), intent(in) :: map 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) use psb_const_mod use psb_s_vect_mod - import :: psb_ipk_, psb_slinmap_type + import :: psb_slinmap_type implicit none type(psb_slinmap_type), intent(in) :: map real(psb_spk_), intent(in) :: alpha,beta @@ -80,7 +80,7 @@ module psb_s_linmap_mod interface psb_map_Y2X subroutine psb_s_map_Y2X(alpha,x,beta,y,map,info,work) use psb_const_mod - import :: psb_ipk_, psb_slinmap_type + import :: psb_slinmap_type implicit none type(psb_slinmap_type), intent(in) :: map 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) use psb_const_mod use psb_s_vect_mod - import :: psb_ipk_, psb_slinmap_type + import :: psb_slinmap_type implicit none type(psb_slinmap_type), intent(in) :: map real(psb_spk_), intent(in) :: alpha,beta diff --git a/base/modules/psb_s_vect_mod.f90 b/base/modules/psb_s_vect_mod.f90 index baf22d60..c89392c6 100644 --- a/base/modules/psb_s_vect_mod.f90 +++ b/base/modules/psb_s_vect_mod.f90 @@ -163,7 +163,7 @@ contains 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 diff --git a/base/modules/psb_serial_mod.f90 b/base/modules/psb_serial_mod.f90 index 046fac81..331dc33f 100644 --- a/base/modules/psb_serial_mod.f90 +++ b/base/modules/psb_serial_mod.f90 @@ -40,6 +40,63 @@ module psb_serial_mod & psb_gth => psi_gth,& & 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 subroutine psb_ssymbmm(a,b,c,info) @@ -480,6 +537,31 @@ module psb_serial_mod 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) use psb_mat_mod integer(psb_ipk_), intent(in) :: iout diff --git a/base/modules/psb_z_base_mat_mod.f90 b/base/modules/psb_z_base_mat_mod.f90 index 76668dd0..0404abc1 100644 --- a/base/modules/psb_z_base_mat_mod.f90 +++ b/base/modules/psb_z_base_mat_mod.f90 @@ -1170,9 +1170,9 @@ contains if (allocated(a%ja)) deallocate(a%ja) if (allocated(a%val)) deallocate(a%val) call a%set_null() - call a%set_nrows(0) - call a%set_ncols(0) - call a%set_nzeros(0) + call a%set_nrows(izero) + call a%set_ncols(izero) + call a%set_nzeros(izero) return diff --git a/base/modules/psb_z_base_vect_mod.f90 b/base/modules/psb_z_base_vect_mod.f90 index 04b66c68..bca0a227 100644 --- a/base/modules/psb_z_base_vect_mod.f90 +++ b/base/modules/psb_z_base_vect_mod.f90 @@ -150,7 +150,7 @@ contains integer(psb_ipk_) :: info this%v = x - call this%asb(size(x),info) + call this%asb(size(x,kind=psb_ipk_),info) end function constructor diff --git a/base/modules/psb_z_csc_mat_mod.f90 b/base/modules/psb_z_csc_mat_mod.f90 index cf7c1032..46be0d36 100644 --- a/base/modules/psb_z_csc_mat_mod.f90 +++ b/base/modules/psb_z_csc_mat_mod.f90 @@ -527,8 +527,8 @@ contains if (allocated(a%ia)) deallocate(a%ia) if (allocated(a%val)) deallocate(a%val) call a%set_null() - call a%set_nrows(0) - call a%set_ncols(0) + call a%set_nrows(izero) + call a%set_ncols(izero) return diff --git a/base/modules/psb_z_csr_mat_mod.f90 b/base/modules/psb_z_csr_mat_mod.f90 index 1cbb1d33..b3efebb7 100644 --- a/base/modules/psb_z_csr_mat_mod.f90 +++ b/base/modules/psb_z_csr_mat_mod.f90 @@ -527,8 +527,8 @@ contains if (allocated(a%ja)) deallocate(a%ja) if (allocated(a%val)) deallocate(a%val) call a%set_null() - call a%set_nrows(0) - call a%set_ncols(0) + call a%set_nrows(izero) + call a%set_ncols(izero) return diff --git a/base/modules/psb_z_linmap_mod.f90 b/base/modules/psb_z_linmap_mod.f90 index 74fd6d71..234958cf 100644 --- a/base/modules/psb_z_linmap_mod.f90 +++ b/base/modules/psb_z_linmap_mod.f90 @@ -55,7 +55,7 @@ module psb_z_linmap_mod interface psb_map_X2Y subroutine psb_z_map_X2Y(alpha,x,beta,y,map,info,work) use psb_const_mod - import :: psb_ipk_, psb_zlinmap_type + import :: psb_zlinmap_type implicit none type(psb_zlinmap_type), intent(in) :: map 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) use psb_const_mod use psb_z_vect_mod - import :: psb_ipk_, psb_zlinmap_type + import :: psb_zlinmap_type implicit none type(psb_zlinmap_type), intent(in) :: map complex(psb_dpk_), intent(in) :: alpha,beta @@ -80,7 +80,7 @@ module psb_z_linmap_mod interface psb_map_Y2X subroutine psb_z_map_Y2X(alpha,x,beta,y,map,info,work) use psb_const_mod - import :: psb_ipk_, psb_zlinmap_type + import :: psb_zlinmap_type implicit none type(psb_zlinmap_type), intent(in) :: map 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) use psb_const_mod use psb_z_vect_mod - import :: psb_ipk_, psb_zlinmap_type + import :: psb_zlinmap_type implicit none type(psb_zlinmap_type), intent(in) :: map complex(psb_dpk_), intent(in) :: alpha,beta diff --git a/base/modules/psb_z_vect_mod.f90 b/base/modules/psb_z_vect_mod.f90 index 4a9da1a4..546aab04 100644 --- a/base/modules/psb_z_vect_mod.f90 +++ b/base/modules/psb_z_vect_mod.f90 @@ -163,7 +163,7 @@ contains 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 diff --git a/base/psblas/psb_camax.f90 b/base/psblas/psb_camax.f90 index 88b67047..e7c42592 100644 --- a/base/psblas/psb_camax.f90 +++ b/base/psblas/psb_camax.f90 @@ -44,7 +44,7 @@ ! info - integer. Return code ! 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 implicit none @@ -53,25 +53,20 @@ function psb_camax(x,desc_a, info, jx) type(psb_desc_type), intent(in) :: desc_a integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: jx - real(psb_spk_) :: psb_camax + real(psb_spk_) :: res ! locals integer(psb_ipk_) :: ictxt, np, me,& - & err_act, iix, jjx, ix, ijx, m, imax, icamax - real(psb_spk_) :: amax - character(len=20) :: name, ch_err - complex :: cdum - real :: cabs1 - cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) ) + & err_act, iix, jjx, ix, ijx, m, ldx + character(len=20) :: name, ch_err name='psb_camax' if(psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - amax=0.d0 - ictxt=desc_a%get_context() + ictxt = desc_a%get_context() call psb_info(ictxt, me, np) if (np == -1) then @@ -88,8 +83,9 @@ function psb_camax(x,desc_a, info, jx) endif 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 info=psb_err_from_subroutine_ ch_err='psb_chkvect' @@ -105,16 +101,13 @@ function psb_camax(x,desc_a, info, jx) ! compute local max if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then - imax=icamax(desc_a%get_local_rows()-iix+1,x(iix,jjx),1) - amax=cabs1(x(iix+imax-1,jjx)) + res = psb_amax(desc_a%get_local_rows()-iix+1,x(:,jjx)) else - amax = szero + res = szero end if ! compute global max - call psb_amx(ictxt, amax) - - psb_camax=amax + call psb_amx(ictxt, res) call psb_erractionrestore(err_act) return @@ -174,32 +167,27 @@ end function psb_camax ! desc_a - type(psb_desc_type). The communication descriptor. ! 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 implicit none complex(psb_spk_), intent(in) :: x(:) type(psb_desc_type), intent(in) :: desc_a - integer(psb_ipk_), intent(out) :: info - real(psb_spk_) :: psb_camaxv + integer(psb_ipk_), intent(out) :: info + real(psb_spk_) :: res ! locals integer(psb_ipk_) :: ictxt, np, me,& - & err_act, iix, jjx, jx, ix, m, imax, icamax - real(psb_spk_) :: amax - complex(psb_spk_) :: cmax + & err_act, iix, jjx, jx, ix, m, ldx + character(len=20) :: name, ch_err - complex :: cdum - real :: cabs1 - cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) ) name='psb_camaxv' if(psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - amax=0.d0 ictxt=desc_a%get_context() @@ -214,8 +202,9 @@ function psb_camaxv (x,desc_a, info) jx = 1 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 info=psb_err_from_subroutine_ ch_err='psb_chkvect' @@ -231,17 +220,13 @@ function psb_camaxv (x,desc_a, info) ! compute local max if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then - imax=icamax(desc_a%get_local_rows()-iix+1,x(iix),1) - cmax=(x(iix+imax-1)) - amax=cabs1(cmax) + res = psb_amax(desc_a%get_local_rows()-iix+1,x) else - amax = szero + res = szero end if ! compute global max - call psb_amx(ictxt, amax) - - psb_camaxv=amax + call psb_amx(ictxt, res) call psb_erractionrestore(err_act) return @@ -269,12 +254,11 @@ function psb_camax_vect(x, desc_a, info) result(res) real(psb_spk_) :: res type(psb_c_vect_type), intent (inout) :: x type(psb_desc_type), intent (in) :: desc_a - integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(out) :: info ! locals integer(psb_ipk_) :: ictxt, np, me,& - & err_act, iix, jjx, jx, ix, m, imax, isamax - real(psb_spk_) :: amax + & err_act, iix, jjx, jx, ix, m character(len=20) :: name, ch_err name='psb_camaxv' @@ -282,7 +266,6 @@ function psb_camax_vect(x, desc_a, info) result(res) info=psb_success_ call psb_erractionsave(err_act) - amax=szero ictxt=desc_a%get_context() call psb_info(ictxt, me, np) @@ -302,8 +285,7 @@ function psb_camax_vect(x, desc_a, info) result(res) jx = 1 m = desc_a%get_global_rows() - - call psb_chkvect(m,1,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) + 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' @@ -319,15 +301,13 @@ function psb_camax_vect(x, desc_a, info) result(res) ! compute local max 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 - amax = szero + res = szero end if ! compute global max - call psb_amx(ictxt, amax) - - res=amax + call psb_amx(ictxt, res) call psb_erractionrestore(err_act) return @@ -397,27 +377,21 @@ subroutine psb_camaxvs(res,x,desc_a, info) complex(psb_spk_), intent(in) :: x(:) type(psb_desc_type), intent(in) :: desc_a - integer(psb_ipk_), intent(out) :: info - real(psb_spk_), intent(out) :: res + integer(psb_ipk_), intent(out) :: info + real(psb_spk_), intent(out) :: res ! locals integer(psb_ipk_) :: ictxt, np, me,& - & err_act, iix, jjx, ix, ijx, m, imax, icamax - real(psb_spk_) :: amax - character(len=20) :: name, ch_err - complex(psb_spk_) :: cmax - complex :: cdum - real :: cabs1 - cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) ) + & err_act, iix, jjx, ix, ijx, m, ldx + character(len=20) :: name, ch_err name='psb_camaxvs' if(psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - amax=0.d0 - ictxt=desc_a%get_context() + ictxt = desc_a%get_context() call psb_info(ictxt, me, np) if (np == -1) then @@ -430,7 +404,8 @@ subroutine psb_camaxvs(res,x,desc_a, info) ijx=1 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 info=psb_err_from_subroutine_ ch_err='psb_chkvect' @@ -446,17 +421,13 @@ subroutine psb_camaxvs(res,x,desc_a, info) ! compute local max if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then - imax=icamax(desc_a%get_local_rows()-iix+1,x(iix),1) - cmax=(x(iix+imax-1)) - amax=cabs1(cmax) + res = psb_amax(desc_a%get_local_rows()-iix+1,x) else - amax = szero + res = szero end if ! compute global max - call psb_amx(ictxt, amax) - - res = amax + call psb_amx(ictxt, res) call psb_erractionrestore(err_act) return @@ -528,21 +499,14 @@ subroutine psb_cmamaxs(res,x,desc_a, info,jx) ! locals integer(psb_ipk_) :: ictxt, np, me,& - & err_act, iix, jjx, ix, ijx, m, imax, i, k, icamax - real(psb_spk_) :: amax + & err_act, iix, jjx, ix, ijx, m, ldx, i, k 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' if (psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - amax=0.d0 - ictxt=desc_a%get_context() 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() k = min(size(x,2),size(res,1)) - - 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 info=psb_err_from_subroutine_ ch_err='psb_chkvect' @@ -576,20 +540,17 @@ subroutine psb_cmamaxs(res,x,desc_a, info,jx) goto 9999 end if + res(1:k) = szero ! compute local max if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then - do i=1,k - imax=icamax(desc_a%get_local_rows()-iix+1,x(iix,jjx+i-1),1) - cmax=(x(iix+imax-1,jjx+i-1)) - res(i)=cabs1(cmax) - end do - else - amax = szero + do i=1,k + res(i) = psb_amax(desc_a%get_local_rows()-iix+1,x(:,jjx+i-1)) + end do end if ! compute global max call psb_amx(ictxt, res(1:k)) - + call psb_erractionrestore(err_act) return diff --git a/base/psblas/psb_casum.f90 b/base/psblas/psb_casum.f90 index 0af133c3..96d166d3 100644 --- a/base/psblas/psb_casum.f90 +++ b/base/psblas/psb_casum.f90 @@ -44,7 +44,7 @@ ! info - integer. Return code ! 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 implicit none @@ -53,24 +53,18 @@ function psb_casum (x,desc_a, info, jx) type(psb_desc_type), intent(in) :: desc_a integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: jx - real(psb_spk_) :: psb_casum + real(psb_spk_) :: res ! locals integer(psb_ipk_) :: ictxt, np, me, & - & err_act, iix, jjx, ix, ijx, m, i, idx, ndm - real(psb_spk_) :: asum, scasum + & err_act, iix, jjx, ix, ijx, m, i, idx, ndm, ldx 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' if(psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - asum=0.d0 ictxt=desc_a%get_context() @@ -89,9 +83,9 @@ function psb_casum (x,desc_a, info, jx) endif m = desc_a%get_global_rows() - + ldx = size(x,1) ! 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 info=psb_err_from_subroutine_ ch_err='psb_chkvect' @@ -106,29 +100,21 @@ function psb_casum (x,desc_a, info, jx) end if ! compute local max - if ((m /= 0)) then - if(desc_a%get_local_rows() > 0) then - 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) + if(desc_a%get_local_rows() > 0) then + res = psb_asum(desc_a%get_local_rows()-iix+1,x(:,jjx)) + + ! 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 + else - asum=0.0 + res = szero end if - - - psb_casum=asum + ! compute global sum + call psb_sum(ictxt, res) call psb_erractionrestore(err_act) return @@ -137,31 +123,25 @@ function psb_casum (x,desc_a, info, jx) call psb_erractionrestore(err_act) if (err_act == psb_act_abort_) then - call psb_error(ictxt) - return + call psb_error(ictxt) + return end if return end function psb_casum function psb_casum_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_c_vect_mod + use psb_base_mod, psb_protect_name => psb_casum_vect implicit none real(psb_spk_) :: res type(psb_c_vect_type), intent (inout) :: x type(psb_desc_type), intent (in) :: desc_a - integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(out) :: info ! locals integer(psb_ipk_) :: ictxt, np, me,& & err_act, iix, jjx, jx, ix, m, imax - real(psb_spk_) :: asum character(len=20) :: name, ch_err name='psb_casumv' @@ -169,7 +149,6 @@ function psb_casum_vect(x, desc_a, info) result(res) info=psb_success_ call psb_erractionsave(err_act) - asum=0.d0 ictxt=desc_a%get_context() @@ -191,8 +170,7 @@ function psb_casum_vect(x, desc_a, info) result(res) jx = 1 m = desc_a%get_global_rows() - - call psb_chkvect(m,1,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) + 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' @@ -208,15 +186,13 @@ function psb_casum_vect(x, desc_a, info) result(res) ! compute local max 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 - asum = szero + res = szero end if ! compute global sum - call psb_sum(ictxt, asum) - - res=asum + call psb_sum(ictxt, res) call psb_erractionrestore(err_act) return @@ -276,33 +252,26 @@ end function psb_casum_vect ! desc_a - type(psb_desc_type). The communication descriptor. ! 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 implicit none complex(psb_spk_), intent(in) :: x(:) - type(psb_desc_type), intent(in) :: desc_a - integer(psb_ipk_), intent(out) :: info - real(psb_spk_) :: psb_casumv + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + real(psb_spk_) :: res ! locals integer(psb_ipk_) :: ictxt, np, me,& - & err_act, iix, jjx, jx, ix, m, i, idx, ndm - real(psb_spk_) :: asum, scasum + & err_act, iix, jjx, jx, ix, m, i, idx, ndm, ldx 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' if(psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - asum=0.d0 - ictxt=desc_a%get_context() call psb_info(ictxt, me, np) @@ -316,9 +285,9 @@ function psb_casumv(x,desc_a, info) jx=1 m = desc_a%get_global_rows() - + ldx = size(x,1) ! 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 info=psb_err_from_subroutine_ ch_err='psb_chkvect' @@ -333,28 +302,22 @@ function psb_casumv(x,desc_a, info) end if ! compute local max - if ((m /= 0)) then - if(desc_a%get_local_rows() > 0) then - asum=scasum(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))*cabs1(x(idx)) - end do - - else - asum=0.d0 - end if - - ! compute global sum - call psb_sum(ictxt, asum) + if(desc_a%get_local_rows() > 0) then + 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 - asum=0.d0 + res = szero end if - psb_casumv=asum + ! compute global sum + call psb_sum(ictxt, res) call psb_erractionrestore(err_act) return @@ -419,28 +382,21 @@ subroutine psb_casumvs(res,x,desc_a, info) implicit none - complex(psb_spk_), intent(in) :: x(:) - real(psb_spk_), intent(out) :: res - type(psb_desc_type), intent(in) :: desc_a - integer(psb_ipk_), intent(out) :: info + complex(psb_spk_), intent(in) :: x(:) + real(psb_spk_), intent(out) :: res + 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, ix, jx, m, i, idx, ndm - real(psb_spk_) :: asum, scasum + & err_act, iix, jjx, ix, jx, m, i, idx, ndm, ldx 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' if(psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - asum=0.d0 - ictxt=desc_a%get_context() call psb_info(ictxt, me, np) @@ -454,9 +410,9 @@ subroutine psb_casumvs(res,x,desc_a, info) jx = 1 m = desc_a%get_global_rows() - + ldx = size(x,1) ! 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 info=psb_err_from_subroutine_ ch_err='psb_chkvect' @@ -471,29 +427,22 @@ subroutine psb_casumvs(res,x,desc_a, info) end if ! compute local max - if ((m /= 0)) then - if(desc_a%get_local_rows() > 0) then - asum=scasum(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))*cabs1(x(idx)) - end do - - else - asum=0.d0 - end if - - ! compute global sum - call psb_sum(ictxt,asum) + if(desc_a%get_local_rows() > 0) then + 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 - asum=0.d0 + res = szero end if - - res = asum + ! compute global sum + call psb_sum(ictxt,res) call psb_erractionrestore(err_act) return diff --git a/base/psblas/psb_damax.f90 b/base/psblas/psb_damax.f90 index 220c672a..8acafde0 100644 --- a/base/psblas/psb_damax.f90 +++ b/base/psblas/psb_damax.f90 @@ -39,39 +39,34 @@ ! where sub( X ) denotes X(1:N,JX:). ! ! Arguments: -! x(:,:) - real The input vector. -! desc_a - type(psb_desc_type). The communication descriptor. -! info - integer. Return code -! jx - integer(optional). The column offset. +! x(:,:) - real The input vector. +! desc_a - type(psb_desc_type). The communication descriptor. +! info - integer. Return code +! jx - integer(optional). The column offset. ! -function psb_damax (x,desc_a, info, jx) - use psb_penv_mod - use psb_serial_mod - use psb_descriptor_type - use psb_check_mod - use psb_error_mod +function psb_damax(x,desc_a, info, jx) result(res) + use psb_base_mod, psb_protect_name => psb_damax + implicit none - real(psb_dpk_), intent(in) :: x(:,:) - type(psb_desc_type), intent(in) :: desc_a - integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_), optional, intent(in) :: jx - real(psb_dpk_) :: psb_damax + real(psb_dpk_), intent(in) :: x(:,:) + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: jx + real(psb_dpk_) :: res ! locals integer(psb_ipk_) :: ictxt, np, me,& - & err_act, iix, jjx, ix, ijx, m, imax, idamax - real(psb_dpk_) :: amax - character(len=20) :: name, ch_err + & err_act, iix, jjx, ix, ijx, m, ldx + character(len=20) :: name, ch_err name='psb_damax' if(psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - amax=0.d0 - ictxt=desc_a%get_context() + ictxt = desc_a%get_context() call psb_info(ictxt, me, np) if (np == -1) then @@ -79,42 +74,40 @@ function psb_damax (x,desc_a, info, jx) call psb_errpush(info,name) goto 9999 endif - + ix = 1 if (present(jx)) then - ijx = jx + ijx = jx else - ijx = 1 + ijx = 1 endif 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 - info=psb_err_from_subroutine_ - ch_err='psb_chkvect' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 + 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 + 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 - imax=idamax(desc_a%get_local_rows()-iix+1,x(iix,jjx),1) - amax=abs(x(iix+imax-1,jjx)) + res = psb_amax(desc_a%get_local_rows()-iix+1,x(:,jjx)) else - amax = dzero + res = dzero end if - + ! compute global max - call psb_amx(ictxt, amax) - - psb_damax=amax + call psb_amx(ictxt, res) call psb_erractionrestore(err_act) return @@ -123,8 +116,8 @@ function psb_damax (x,desc_a, info, jx) call psb_erractionrestore(err_act) if (err_act == psb_act_abort_) then - call psb_error(ictxt) - return + call psb_error(ictxt) + return end if return end function psb_damax @@ -162,7 +155,7 @@ end function psb_damax !!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE !!$ POSSIBILITY OF SUCH DAMAGE. !!$ -!!$ +!!$ ! ! Function: psb_damaxv ! Searches the absolute max of X. @@ -170,27 +163,24 @@ end function psb_damax ! normi := max(abs(X(i)) ! ! Arguments: -! x(:) - real The input vector. -! desc_a - type(psb_desc_type). The communication descriptor. -! info - integer. Return code +! x(:) - real The input vector. +! desc_a - type(psb_desc_type). The communication descriptor. +! info - integer. Return code ! -function psb_damaxv (x,desc_a, info) - use psb_penv_mod - use psb_serial_mod - use psb_descriptor_type - use psb_check_mod - use psb_error_mod +function psb_damaxv (x,desc_a, info) result(res) + use psb_base_mod, psb_protect_name => psb_damaxv + implicit none - real(psb_dpk_), intent(in) :: x(:) - type(psb_desc_type), intent(in) :: desc_a - integer(psb_ipk_), intent(out) :: info - real(psb_dpk_) :: psb_damaxv + real(psb_dpk_), intent(in) :: x(:) + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + real(psb_dpk_) :: res ! locals integer(psb_ipk_) :: ictxt, np, me,& - & err_act, iix, jjx, jx, ix, m, imax, idamax - real(psb_dpk_) :: amax + & err_act, iix, jjx, jx, ix, m, ldx + character(len=20) :: name, ch_err name='psb_damaxv' @@ -198,7 +188,6 @@ function psb_damaxv (x,desc_a, info) info=psb_success_ call psb_erractionsave(err_act) - amax=0.d0 ictxt=desc_a%get_context() @@ -208,38 +197,36 @@ function psb_damaxv (x,desc_a, info) call psb_errpush(info,name) goto 9999 endif - + ix = 1 jx = 1 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 - info=psb_err_from_subroutine_ - ch_err='psb_chkvect' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 + 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 + 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 - imax=idamax(desc_a%get_local_rows()-iix+1,x(iix),1) - amax=abs(x(iix+imax-1)) + res = psb_amax(desc_a%get_local_rows()-iix+1,x) else - amax = dzero + res = dzero end if - + ! compute global max - call psb_amx(ictxt, amax) - - psb_damaxv=amax + call psb_amx(ictxt, res) call psb_erractionrestore(err_act) return @@ -248,12 +235,13 @@ function psb_damaxv (x,desc_a, info) call psb_erractionrestore(err_act) if (err_act == psb_act_abort_) then - call psb_error(ictxt) - return + call psb_error(ictxt) + return end if return end function psb_damaxv + function psb_damax_vect(x, desc_a, info) result(res) use psb_penv_mod use psb_serial_mod @@ -266,20 +254,18 @@ function psb_damax_vect(x, desc_a, info) result(res) real(psb_dpk_) :: res type(psb_d_vect_type), intent (inout) :: x type(psb_desc_type), intent (in) :: desc_a - integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(out) :: info ! locals integer(psb_ipk_) :: ictxt, np, me,& - & err_act, iix, jjx, jx, ix, m, imax, idamax - real(psb_dpk_) :: amax - character(len=20) :: name, ch_err + & err_act, iix, jjx, jx, ix, m + character(len=20) :: name, ch_err name='psb_damaxv' if(psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - amax=0.d0 ictxt=desc_a%get_context() call psb_info(ictxt, me, np) @@ -299,8 +285,7 @@ function psb_damax_vect(x, desc_a, info) result(res) jx = 1 m = desc_a%get_global_rows() - - call psb_chkvect(m,1,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) + 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' @@ -316,15 +301,13 @@ function psb_damax_vect(x, desc_a, info) result(res) ! compute local max 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 - amax = dzero + res = dzero end if ! compute global max - call psb_amx(ictxt, amax) - - res=amax + call psb_amx(ictxt, res) call psb_erractionrestore(err_act) return @@ -381,35 +364,34 @@ end function psb_damax_vect ! where sub( X ) denotes X(1:N,JX:). ! ! Arguments: -! res - real. The result. -! x(:,:) - real The input vector. -! desc_a - type(psb_desc_type). The communication descriptor. -! info - integer. Return code -! jx - integer(optional). The column offset. +! res - real The result. +! x(:,:) - real The input vector. +! desc_a - type(psb_desc_type). The communication descriptor. +! info - integer. Return code +! jx - integer(optional). The column offset. ! -subroutine psb_damaxvs (res,x,desc_a, info) - use psb_base_mod, psb_protect_name => psb_damaxvs +subroutine psb_damaxvs(res,x,desc_a, info) + use psb_base_mod, psb_protect_name => psb_damaxvs + implicit none - real(psb_dpk_), intent(in) :: x(:) - type(psb_desc_type), intent(in) :: desc_a - integer(psb_ipk_), intent(out) :: info - real(psb_dpk_), intent(out) :: res + real(psb_dpk_), intent(in) :: x(:) + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + real(psb_dpk_), intent(out) :: res ! locals integer(psb_ipk_) :: ictxt, np, me,& - & err_act, iix, jjx, ix, ijx, m, imax, idamax - real(psb_dpk_) :: amax - character(len=20) :: name, ch_err + & err_act, iix, jjx, ix, ijx, m, ldx + character(len=20) :: name, ch_err name='psb_damaxvs' if(psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - amax=0.d0 - ictxt=desc_a%get_context() + ictxt = desc_a%get_context() call psb_info(ictxt, me, np) if (np == -1) then @@ -422,33 +404,30 @@ subroutine psb_damaxvs (res,x,desc_a, info) ijx=1 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 - info=psb_err_from_subroutine_ - ch_err='psb_chkvect' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 + 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 + 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 - imax=idamax(desc_a%get_local_rows()-iix+1,x(iix),1) - amax=abs(x(iix+imax-1)) + res = psb_amax(desc_a%get_local_rows()-iix+1,x) else - amax = dzero + res = dzero end if - + ! compute global max - call psb_amx(ictxt, amax) - - res = amax + call psb_amx(ictxt, res) call psb_erractionrestore(err_act) return @@ -457,8 +436,8 @@ subroutine psb_damaxvs (res,x,desc_a, info) call psb_erractionrestore(err_act) if (err_act == psb_act_abort_) then - call psb_error(ictxt) - return + call psb_error(ictxt) + return end if return end subroutine psb_damaxvs @@ -502,34 +481,32 @@ end subroutine psb_damaxvs ! normi := max(abs(X(i)) ! ! Arguments: -! res(:) - real The result. -! x(:,:) - real The input vector. -! desc_a - type(psb_desc_type). The communication descriptor. -! info - integer. Return code +! res(:) - real. The result. +! x(:,:) - real The input vector. +! desc_a - type(psb_desc_type). The communication descriptor. +! info - integer. Return code ! -subroutine psb_dmamaxs (res,x,desc_a, info,jx) - use psb_base_mod, psb_protect_name => psb_dmamaxs +subroutine psb_dmamaxs(res,x,desc_a, info,jx) + use psb_base_mod, psb_protect_name => psb_dmamaxs + implicit none - real(psb_dpk_), intent(in) :: x(:,:) - type(psb_desc_type), intent(in) :: desc_a - integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_), optional, intent(in) :: jx - real(psb_dpk_), intent(out) :: res(:) + real(psb_dpk_), intent(in) :: x(:,:) + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: jx + real(psb_dpk_), intent(out) :: res(:) ! locals integer(psb_ipk_) :: ictxt, np, me,& - & err_act, iix, jjx, ix, ijx, m, imax, i, k, idamax - real(psb_dpk_) :: amax + & err_act, iix, jjx, ix, ijx, m, ldx, i, k character(len=20) :: name, ch_err name='psb_dmamaxs' - if (psb_errstatus_fatal()) return + if (psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - amax=0.d0 - ictxt=desc_a%get_context() 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() k = min(size(x,2),size(res,1)) - - 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 info=psb_err_from_subroutine_ ch_err='psb_chkvect' @@ -563,19 +540,17 @@ subroutine psb_dmamaxs (res,x,desc_a, info,jx) goto 9999 end if + res(1:k) = dzero ! compute local max if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then - do i=1,k - imax=idamax(desc_a%get_local_rows()-iix+1,x(iix,jjx+i-1),1) - res(i)=abs(x(iix+imax-1,jjx+i-1)) - end do - else - amax = dzero + do i=1,k + res(i) = psb_amax(desc_a%get_local_rows()-iix+1,x(:,jjx+i-1)) + end do end if ! compute global max call psb_amx(ictxt, res(1:k)) - + call psb_erractionrestore(err_act) return diff --git a/base/psblas/psb_dasum.f90 b/base/psblas/psb_dasum.f90 index d43afa25..17510eee 100644 --- a/base/psblas/psb_dasum.f90 +++ b/base/psblas/psb_dasum.f90 @@ -39,38 +39,32 @@ ! where sub( X ) denotes X(1:N,JX:). ! ! Arguments: -! x(:,:) - real The input vector. -! desc_a - type(psb_desc_type). The communication descriptor. -! info - integer. Return code -! jx - integer(optional). The column offset. +! x(:,:) - real The input vector. +! desc_a - type(psb_desc_type). The communication descriptor. +! info - integer. Return code +! 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 - real(psb_dpk_), intent(in) :: x(:,:) + real(psb_dpk_), intent(in) :: x(:,:) type(psb_desc_type), intent(in) :: desc_a integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: jx - real(psb_dpk_) :: psb_dasum + real(psb_dpk_) :: res ! locals - integer(psb_ipk_) :: ictxt, np, me, err_act, & - & iix, jjx, ix, ijx, m, i, idx, ndm - real(psb_dpk_) :: asum, dasum - character(len=20) :: name, ch_err + integer(psb_ipk_) :: ictxt, np, me, & + & err_act, iix, jjx, ix, ijx, m, i, idx, ndm, ldx + character(len=20) :: name, ch_err name='psb_dasum' - if (psb_get_errstatus() /= 0) return + if(psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - asum=0.d0 ictxt=desc_a%get_context() @@ -89,9 +83,9 @@ function psb_dasum (x,desc_a, info, jx) endif m = desc_a%get_global_rows() - + ldx = size(x,1) ! 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 info=psb_err_from_subroutine_ ch_err='psb_chkvect' @@ -106,31 +100,21 @@ function psb_dasum (x,desc_a, info, jx) end if ! compute local max - if ((m /= 0)) then - if(desc_a%get_local_rows() > 0) then - 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 + if(desc_a%get_local_rows() > 0) then + res = psb_asum(desc_a%get_local_rows()-iix+1,x(:,jjx)) + ! 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) return @@ -146,73 +130,25 @@ function psb_dasum (x,desc_a, info, jx) end function psb_dasum -!!$ -!!$ 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) - - use psb_serial_mod - use psb_descriptor_type - use psb_check_mod - use psb_error_mod - use psb_penv_mod +function psb_dasum_vect(x, desc_a, info) result(res) + use psb_base_mod, psb_protect_name => psb_dasum_vect implicit none - real(psb_dpk_), intent(in) :: x(:) - type(psb_desc_type), intent(in) :: desc_a - integer(psb_ipk_), intent(out) :: info - real(psb_dpk_) :: psb_dasumv + real(psb_dpk_) :: res + type(psb_d_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, i, idx, ndm - real(psb_dpk_) :: asum, dasum - character(len=20) :: name, ch_err + integer(psb_ipk_) :: ictxt, np, me,& + & err_act, iix, jjx, jx, ix, m, imax + character(len=20) :: name, ch_err name='psb_dasumv' - if(psb_get_errstatus() /= 0) return + if (psb_errstatus_fatal()) return info=psb_success_ call psb_erractionsave(err_act) - asum=0.d0 ictxt=desc_a%get_context() @@ -223,13 +159,18 @@ function psb_dasumv (x,desc_a, info) 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 + jx = 1 m = desc_a%get_global_rows() - - ! check vector correctness - call psb_chkvect(m,1,size(x),ix,jx,desc_a,info,iix,jjx) + 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' @@ -244,29 +185,14 @@ function psb_dasumv (x,desc_a, info) end if ! compute local max - if ((m /= 0)) then - if(desc_a%get_local_rows() > 0) then - 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 - asum=0.d0 + if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then + res = x%asum(desc_a%get_local_rows()) + else + res = dzero end if - psb_dasumv=asum + ! compute global sum + call psb_sum(ictxt, res) call psb_erractionrestore(err_act) return @@ -279,36 +205,73 @@ function psb_dasumv (x,desc_a, info) return end if 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 - real(psb_dpk_) :: res - type(psb_d_vect_type), intent (inout) :: x - type(psb_desc_type), intent (in) :: desc_a - integer(psb_ipk_), intent(out) :: info + real(psb_dpk_), intent(in) :: x(:) + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + real(psb_dpk_) :: res ! locals integer(psb_ipk_) :: ictxt, np, me,& - & err_act, iix, jjx, jx, ix, m, imax - real(psb_dpk_) :: asum + & err_act, iix, jjx, jx, ix, m, i, idx, ndm, ldx character(len=20) :: name, ch_err name='psb_dasumv' - if (psb_errstatus_fatal()) return + if(psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - asum=0.d0 - ictxt=desc_a%get_context() call psb_info(ictxt, me, np) @@ -318,19 +281,13 @@ function psb_dasum_vect(x, desc_a, info) result(res) 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 + jx=1 m = desc_a%get_global_rows() - - call psb_chkvect(m,1,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) + ldx = size(x,1) + ! check vector correctness + call psb_chkvect(m,ione,ldx,ix,jx,desc_a,info,iix,jjx) if(info /= psb_success_) then info=psb_err_from_subroutine_ ch_err='psb_chkvect' @@ -345,16 +302,22 @@ function psb_dasum_vect(x, desc_a, info) result(res) end if ! compute local max - if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then - asum=x%asum(desc_a%get_local_rows()) - else - asum = dzero + if(desc_a%get_local_rows() > 0) then + 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 + res = dzero end if ! compute global sum - call psb_sum(ictxt, asum) - - res=asum + call psb_sum(ictxt, res) call psb_erractionrestore(err_act) return @@ -367,9 +330,7 @@ function psb_dasum_vect(x, desc_a, info) result(res) return end if return - -end function psb_dasum_vect - +end function psb_dasumv !!$ @@ -402,7 +363,7 @@ end function psb_dasum_vect !!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE !!$ POSSIBILITY OF SUCH DAMAGE. !!$ -!!$ +!!$ ! ! Subroutine: psb_dasumvs ! Computes norm1 of X @@ -410,33 +371,32 @@ end function psb_dasum_vect ! norm1 := sum(X(i)) ! ! Arguments: -! res - real The result. -! x(:) - real The input vector. -! desc_a - type(psb_desc_type). The communication descriptor. -! info - integer. Return code -! jx - integer(optional). The column offset. +! res - real. The result. +! x(:) - real The input vector. +! desc_a - type(psb_desc_type). The communication descriptor. +! info - integer. Return code +! jx - integer(optional). The column offset. ! subroutine psb_dasumvs(res,x,desc_a, info) use psb_base_mod, psb_protect_name => psb_dasumvs + implicit none - real(psb_dpk_), intent(in) :: x(:) - real(psb_dpk_), intent(out) :: res - type(psb_desc_type), intent(in) :: desc_a - integer(psb_ipk_), intent(out) :: info + real(psb_dpk_), intent(in) :: x(:) + real(psb_dpk_), intent(out) :: res + 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, ix, jx, m, i, idx, ndm - real(psb_dpk_) :: asum, dasum - character(len=20) :: name, ch_err + integer(psb_ipk_) :: ictxt, np, me,& + & err_act, iix, jjx, ix, jx, m, i, idx, ndm, ldx + character(len=20) :: name, ch_err name='psb_dasumvs' - if (psb_errstatus_fatal()) return + if(psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - asum=0.d0 - ictxt=desc_a%get_context() call psb_info(ictxt, me, np) @@ -450,9 +410,9 @@ subroutine psb_dasumvs(res,x,desc_a, info) jx = 1 m = desc_a%get_global_rows() - + ldx = size(x,1) ! 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 info=psb_err_from_subroutine_ ch_err='psb_chkvect' @@ -467,31 +427,22 @@ subroutine psb_dasumvs(res,x,desc_a, info) end if ! compute local max - if ((m /= 0)) then - if(desc_a%get_local_rows() > 0) then - 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 + if(desc_a%get_local_rows() > 0) then + 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 - asum=0.d0 + res = dzero end if - - res = asum + ! compute global sum + call psb_sum(ictxt,res) call psb_erractionrestore(err_act) return @@ -505,7 +456,3 @@ subroutine psb_dasumvs(res,x,desc_a, info) end if return end subroutine psb_dasumvs - - - - diff --git a/base/psblas/psb_samax.f90 b/base/psblas/psb_samax.f90 index 9388fb39..aa87e0af 100644 --- a/base/psblas/psb_samax.f90 +++ b/base/psblas/psb_samax.f90 @@ -39,39 +39,34 @@ ! where sub( X ) denotes X(1:N,JX:). ! ! Arguments: -! x(:,:) - real The input vector. -! desc_a - type(psb_desc_type). The communication descriptor. -! info - integer. Return code -! jx - integer(optional). The column offset. +! x(:,:) - real The input vector. +! desc_a - type(psb_desc_type). The communication descriptor. +! info - integer. Return code +! jx - integer(optional). The column offset. ! -function psb_samax (x,desc_a, info, jx) - use psb_penv_mod - use psb_serial_mod - use psb_descriptor_type - use psb_check_mod - use psb_error_mod +function psb_samax(x,desc_a, info, jx) result(res) + use psb_base_mod, psb_protect_name => psb_samax + implicit none - real(psb_spk_), intent(in) :: x(:,:) - type(psb_desc_type), intent(in) :: desc_a - integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_), optional, intent(in) :: jx - real(psb_spk_) :: psb_samax + real(psb_spk_), intent(in) :: x(:,:) + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: jx + real(psb_spk_) :: res ! locals integer(psb_ipk_) :: ictxt, np, me,& - & err_act, iix, jjx, ix, ijx, m, imax, isamax - real(psb_spk_) :: amax - character(len=20) :: name, ch_err + & err_act, iix, jjx, ix, ijx, m, ldx + character(len=20) :: name, ch_err name='psb_samax' if(psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - amax=0.d0 - ictxt=desc_a%get_context() + ictxt = desc_a%get_context() call psb_info(ictxt, me, np) if (np == -1) then @@ -79,42 +74,40 @@ function psb_samax (x,desc_a, info, jx) call psb_errpush(info,name) goto 9999 endif - + ix = 1 if (present(jx)) then - ijx = jx + ijx = jx else - ijx = 1 + ijx = 1 endif 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 - info=psb_err_from_subroutine_ - ch_err='psb_chkvect' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 + 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 + 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 - imax=isamax(desc_a%get_local_rows()-iix+1,x(iix,jjx),1) - amax=abs(x(iix+imax-1,jjx)) + res = psb_amax(desc_a%get_local_rows()-iix+1,x(:,jjx)) else - amax = szero + res = szero end if - + ! compute global max - call psb_amx(ictxt, amax) - - psb_samax=amax + call psb_amx(ictxt, res) call psb_erractionrestore(err_act) return @@ -123,8 +116,8 @@ function psb_samax (x,desc_a, info, jx) call psb_erractionrestore(err_act) if (err_act == psb_act_abort_) then - call psb_error(ictxt) - return + call psb_error(ictxt) + return end if return end function psb_samax @@ -162,7 +155,7 @@ end function psb_samax !!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE !!$ POSSIBILITY OF SUCH DAMAGE. !!$ -!!$ +!!$ ! ! Function: psb_samaxv ! Searches the absolute max of X. @@ -170,27 +163,24 @@ end function psb_samax ! normi := max(abs(X(i)) ! ! Arguments: -! x(:) - real The input vector. -! desc_a - type(psb_desc_type). The communication descriptor. -! info - integer. Return code +! x(:) - real The input vector. +! desc_a - type(psb_desc_type). The communication descriptor. +! info - integer. Return code ! -function psb_samaxv (x,desc_a, info) - use psb_penv_mod - use psb_serial_mod - use psb_descriptor_type - use psb_check_mod - use psb_error_mod +function psb_samaxv (x,desc_a, info) result(res) + use psb_base_mod, psb_protect_name => psb_samaxv + implicit none - real(psb_spk_), intent(in) :: x(:) - type(psb_desc_type), intent(in) :: desc_a - integer(psb_ipk_), intent(out) :: info - real(psb_spk_) :: psb_samaxv + real(psb_spk_), intent(in) :: x(:) + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + real(psb_spk_) :: res ! locals integer(psb_ipk_) :: ictxt, np, me,& - & err_act, iix, jjx, jx, ix, m, imax, isamax - real(psb_spk_) :: amax + & err_act, iix, jjx, jx, ix, m, ldx + character(len=20) :: name, ch_err name='psb_samaxv' @@ -198,7 +188,6 @@ function psb_samaxv (x,desc_a, info) info=psb_success_ call psb_erractionsave(err_act) - amax=0.d0 ictxt=desc_a%get_context() @@ -208,38 +197,36 @@ function psb_samaxv (x,desc_a, info) call psb_errpush(info,name) goto 9999 endif - + ix = 1 jx = 1 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 - info=psb_err_from_subroutine_ - ch_err='psb_chkvect' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 + 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 + 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 - imax=isamax(desc_a%get_local_rows()-iix+1,x(iix),1) - amax=abs(x(iix+imax-1)) + res = psb_amax(desc_a%get_local_rows()-iix+1,x) else - amax = szero + res = szero end if - + ! compute global max - call psb_amx(ictxt, amax) - - psb_samaxv=amax + call psb_amx(ictxt, res) call psb_erractionrestore(err_act) return @@ -248,8 +235,8 @@ function psb_samaxv (x,desc_a, info) call psb_erractionrestore(err_act) if (err_act == psb_act_abort_) then - call psb_error(ictxt) - return + call psb_error(ictxt) + return end if return end function psb_samaxv @@ -267,20 +254,18 @@ function psb_samax_vect(x, desc_a, info) result(res) real(psb_spk_) :: res type(psb_s_vect_type), intent (inout) :: x type(psb_desc_type), intent (in) :: desc_a - integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(out) :: info ! locals integer(psb_ipk_) :: ictxt, np, me,& - & err_act, iix, jjx, jx, ix, m, imax, isamax - real(psb_spk_) :: amax - character(len=20) :: name, ch_err + & err_act, iix, jjx, jx, ix, m + character(len=20) :: name, ch_err name='psb_samaxv' if(psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - amax=0.d0 ictxt=desc_a%get_context() call psb_info(ictxt, me, np) @@ -300,8 +285,7 @@ function psb_samax_vect(x, desc_a, info) result(res) jx = 1 m = desc_a%get_global_rows() - - call psb_chkvect(m,1,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) + 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' @@ -317,15 +301,13 @@ function psb_samax_vect(x, desc_a, info) result(res) ! compute local max 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 - amax = szero + res = szero end if ! compute global max - call psb_amx(ictxt, amax) - - res=amax + call psb_amx(ictxt, res) call psb_erractionrestore(err_act) return @@ -382,35 +364,34 @@ end function psb_samax_vect ! where sub( X ) denotes X(1:N,JX:). ! ! Arguments: -! res - real. The result. -! x(:,:) - real The input vector. -! desc_a - type(psb_desc_type). The communication descriptor. -! info - integer. Return code -! jx - integer(optional). The column offset. +! res - real The result. +! x(:,:) - real The input vector. +! desc_a - type(psb_desc_type). The communication descriptor. +! info - integer. Return code +! jx - integer(optional). The column offset. ! -subroutine psb_samaxvs (res,x,desc_a, info) - use psb_base_mod, psb_protect_name => psb_samaxvs +subroutine psb_samaxvs(res,x,desc_a, info) + use psb_base_mod, psb_protect_name => psb_samaxvs + implicit none - real(psb_spk_), intent(in) :: x(:) - type(psb_desc_type), intent(in) :: desc_a - integer(psb_ipk_), intent(out) :: info - real(psb_spk_), intent(out) :: res + real(psb_spk_), intent(in) :: x(:) + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + real(psb_spk_), intent(out) :: res ! locals integer(psb_ipk_) :: ictxt, np, me,& - & err_act, iix, jjx, ix, ijx, m, imax, isamax - real(psb_spk_) :: amax - character(len=20) :: name, ch_err + & err_act, iix, jjx, ix, ijx, m, ldx + character(len=20) :: name, ch_err name='psb_samaxvs' if(psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - amax=0.d0 - ictxt=desc_a%get_context() + ictxt = desc_a%get_context() call psb_info(ictxt, me, np) if (np == -1) then @@ -423,33 +404,30 @@ subroutine psb_samaxvs (res,x,desc_a, info) ijx=1 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 - info=psb_err_from_subroutine_ - ch_err='psb_chkvect' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 + 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 + 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 - imax=isamax(desc_a%get_local_rows()-iix+1,x(iix),1) - amax=abs(x(iix+imax-1)) + res = psb_amax(desc_a%get_local_rows()-iix+1,x) else - amax = szero + res = szero end if - + ! compute global max - call psb_amx(ictxt, amax) - - res = amax + call psb_amx(ictxt, res) call psb_erractionrestore(err_act) return @@ -458,8 +436,8 @@ subroutine psb_samaxvs (res,x,desc_a, info) call psb_erractionrestore(err_act) if (err_act == psb_act_abort_) then - call psb_error(ictxt) - return + call psb_error(ictxt) + return end if return end subroutine psb_samaxvs @@ -503,25 +481,25 @@ end subroutine psb_samaxvs ! normi := max(abs(X(i)) ! ! Arguments: -! res(:) - real The result. -! x(:,:) - real The input vector. -! desc_a - type(psb_desc_type). The communication descriptor. -! info - integer. Return code +! res(:) - real. The result. +! x(:,:) - real The input vector. +! desc_a - type(psb_desc_type). The communication descriptor. +! info - integer. Return code ! -subroutine psb_smamaxs (res,x,desc_a, info,jx) - use psb_base_mod, psb_protect_name => psb_smamaxs +subroutine psb_smamaxs(res,x,desc_a, info,jx) + use psb_base_mod, psb_protect_name => psb_smamaxs + implicit none - real(psb_spk_), intent(in) :: x(:,:) - type(psb_desc_type), intent(in) :: desc_a - integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_), optional, intent(in) :: jx - real(psb_spk_), intent(out) :: res(:) + real(psb_spk_), intent(in) :: x(:,:) + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: jx + real(psb_spk_), intent(out) :: res(:) ! locals integer(psb_ipk_) :: ictxt, np, me,& - & err_act, iix, jjx, ix, ijx, m, imax, i, k, isamax - real(psb_spk_) :: amax + & err_act, iix, jjx, ix, ijx, m, ldx, i, k character(len=20) :: name, ch_err name='psb_smamaxs' @@ -529,8 +507,6 @@ subroutine psb_smamaxs (res,x,desc_a, info,jx) info=psb_success_ call psb_erractionsave(err_act) - amax=0.d0 - ictxt=desc_a%get_context() 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() k = min(size(x,2),size(res,1)) - - 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 info=psb_err_from_subroutine_ ch_err='psb_chkvect' @@ -564,19 +540,17 @@ subroutine psb_smamaxs (res,x,desc_a, info,jx) goto 9999 end if + res(1:k) = szero ! compute local max if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then - do i=1,k - imax=isamax(desc_a%get_local_rows()-iix+1,x(iix,jjx+i-1),1) - res(i)=abs(x(iix+imax-1,jjx+i-1)) - end do - else - amax = szero + do i=1,k + res(i) = psb_amax(desc_a%get_local_rows()-iix+1,x(:,jjx+i-1)) + end do end if ! compute global max call psb_amx(ictxt, res(1:k)) - + call psb_erractionrestore(err_act) return diff --git a/base/psblas/psb_sasum.f90 b/base/psblas/psb_sasum.f90 index 3469beaf..b8ce8653 100644 --- a/base/psblas/psb_sasum.f90 +++ b/base/psblas/psb_sasum.f90 @@ -39,38 +39,32 @@ ! where sub( X ) denotes X(1:N,JX:). ! ! Arguments: -! x(:,:) - real The input vector. -! desc_a - type(psb_desc_type). The communication descriptor. -! info - integer. Return code -! jx - integer(optional). The column offset. +! x(:,:) - real The input vector. +! desc_a - type(psb_desc_type). The communication descriptor. +! info - integer. Return code +! 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 - real(psb_spk_), intent(in) :: x(:,:) + real(psb_spk_), intent(in) :: x(:,:) type(psb_desc_type), intent(in) :: desc_a integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: jx - real(psb_spk_) :: psb_sasum + real(psb_spk_) :: res ! locals - integer(psb_ipk_) :: ictxt, np, me, err_act, & - & iix, jjx, ix, ijx, m, i, idx, ndm - real(psb_spk_) :: asum, sasum - character(len=20) :: name, ch_err + integer(psb_ipk_) :: ictxt, np, me, & + & err_act, iix, jjx, ix, ijx, m, i, idx, ndm, ldx + character(len=20) :: name, ch_err name='psb_sasum' if(psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - asum=0.0 ictxt=desc_a%get_context() @@ -89,9 +83,9 @@ function psb_sasum (x,desc_a, info, jx) endif m = desc_a%get_global_rows() - + ldx = size(x,1) ! 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 info=psb_err_from_subroutine_ ch_err='psb_chkvect' @@ -106,31 +100,21 @@ function psb_sasum (x,desc_a, info, jx) end if ! compute local max - if ((m /= 0)) then - if(desc_a%get_local_rows() > 0) then - 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 + if(desc_a%get_local_rows() > 0) then + res = psb_asum(desc_a%get_local_rows()-iix+1,x(:,jjx)) + ! 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) return @@ -146,73 +130,25 @@ function psb_sasum (x,desc_a, info, jx) end function psb_sasum -!!$ -!!$ 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) - - use psb_serial_mod - use psb_descriptor_type - use psb_check_mod - use psb_error_mod - use psb_penv_mod +function psb_sasum_vect(x, desc_a, info) result(res) + use psb_base_mod, psb_protect_name => psb_sasum_vect implicit none - real(psb_spk_), intent(in) :: x(:) - type(psb_desc_type), intent(in) :: desc_a - integer(psb_ipk_), intent(out) :: info - real(psb_spk_) :: psb_sasumv + real(psb_spk_) :: res + type(psb_s_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, i, idx, ndm - real(psb_spk_) :: asum, sasum - character(len=20) :: name, ch_err + integer(psb_ipk_) :: ictxt, np, me,& + & err_act, iix, jjx, jx, ix, m, imax + character(len=20) :: name, ch_err name='psb_sasumv' - if(psb_get_errstatus() /= 0) return + if (psb_errstatus_fatal()) return info=psb_success_ call psb_erractionsave(err_act) - asum=0.0 ictxt=desc_a%get_context() @@ -223,13 +159,18 @@ function psb_sasumv (x,desc_a, info) 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 + jx = 1 m = desc_a%get_global_rows() - - ! check vector correctness - call psb_chkvect(m,1,size(x),ix,jx,desc_a,info,iix,jjx) + 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' @@ -244,29 +185,14 @@ function psb_sasumv (x,desc_a, info) end if ! compute local max - if ((m /= 0)) then - if(desc_a%get_local_rows() > 0) then - 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 - asum=0.0 + if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then + res = x%asum(desc_a%get_local_rows()) + else + res = szero end if - psb_sasumv=asum + ! compute global sum + call psb_sum(ictxt, res) call psb_erractionrestore(err_act) return @@ -279,36 +205,73 @@ function psb_sasumv (x,desc_a, info) return end if 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 - real(psb_spk_) :: res - type(psb_s_vect_type), intent (inout) :: x - type(psb_desc_type), intent (in) :: desc_a - integer(psb_ipk_), intent(out) :: info + real(psb_spk_), intent(in) :: x(:) + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + real(psb_spk_) :: res ! locals integer(psb_ipk_) :: ictxt, np, me,& - & err_act, iix, jjx, jx, ix, m, imax - real(psb_spk_) :: asum + & err_act, iix, jjx, jx, ix, m, i, idx, ndm, ldx character(len=20) :: name, ch_err name='psb_sasumv' - if (psb_errstatus_fatal()) return + if(psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - asum=0.d0 - ictxt=desc_a%get_context() call psb_info(ictxt, me, np) @@ -318,19 +281,13 @@ function psb_sasum_vect(x, desc_a, info) result(res) 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 + jx=1 m = desc_a%get_global_rows() - - call psb_chkvect(m,1,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) + ldx = size(x,1) + ! check vector correctness + call psb_chkvect(m,ione,ldx,ix,jx,desc_a,info,iix,jjx) if(info /= psb_success_) then info=psb_err_from_subroutine_ ch_err='psb_chkvect' @@ -345,16 +302,22 @@ function psb_sasum_vect(x, desc_a, info) result(res) end if ! compute local max - if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then - asum=x%asum(desc_a%get_local_rows()) - else - asum = szero + if(desc_a%get_local_rows() > 0) then + 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 + res = szero end if ! compute global sum - call psb_sum(ictxt, asum) - - res=asum + call psb_sum(ictxt, res) call psb_erractionrestore(err_act) return @@ -367,9 +330,7 @@ function psb_sasum_vect(x, desc_a, info) result(res) return end if return - -end function psb_sasum_vect - +end function psb_sasumv !!$ @@ -402,7 +363,7 @@ end function psb_sasum_vect !!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE !!$ POSSIBILITY OF SUCH DAMAGE. !!$ -!!$ +!!$ ! ! Subroutine: psb_sasumvs ! Computes norm1 of X @@ -410,33 +371,32 @@ end function psb_sasum_vect ! norm1 := sum(X(i)) ! ! Arguments: -! res - real The result. -! x(:) - real The input vector. -! desc_a - type(psb_desc_type). The communication descriptor. -! info - integer. Return code -! jx - integer(optional). The column offset. +! res - real. The result. +! x(:) - real The input vector. +! desc_a - type(psb_desc_type). The communication descriptor. +! info - integer. Return code +! jx - integer(optional). The column offset. ! subroutine psb_sasumvs(res,x,desc_a, info) use psb_base_mod, psb_protect_name => psb_sasumvs + implicit none - real(psb_spk_), intent(in) :: x(:) - real(psb_spk_), intent(out) :: res - type(psb_desc_type), intent(in) :: desc_a - integer(psb_ipk_), intent(out) :: info + real(psb_spk_), intent(in) :: x(:) + real(psb_spk_), intent(out) :: res + 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, ix, jx, m, i, idx, ndm - real(psb_spk_) :: asum, sasum - character(len=20) :: name, ch_err + integer(psb_ipk_) :: ictxt, np, me,& + & err_act, iix, jjx, ix, jx, m, i, idx, ndm, ldx + character(len=20) :: name, ch_err name='psb_sasumvs' if(psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - asum=0.0 - ictxt=desc_a%get_context() call psb_info(ictxt, me, np) @@ -450,9 +410,9 @@ subroutine psb_sasumvs(res,x,desc_a, info) jx = 1 m = desc_a%get_global_rows() - + ldx = size(x,1) ! 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 info=psb_err_from_subroutine_ ch_err='psb_chkvect' @@ -467,31 +427,22 @@ subroutine psb_sasumvs(res,x,desc_a, info) end if ! compute local max - if ((m /= 0)) then - if(desc_a%get_local_rows() > 0) then - 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 + if(desc_a%get_local_rows() > 0) then + 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 - asum=0.0 + res = szero end if - - res = asum + ! compute global sum + call psb_sum(ictxt,res) call psb_erractionrestore(err_act) return diff --git a/base/psblas/psb_zamax.f90 b/base/psblas/psb_zamax.f90 index 5bf7cf8c..d564cda0 100644 --- a/base/psblas/psb_zamax.f90 +++ b/base/psblas/psb_zamax.f90 @@ -44,37 +44,29 @@ ! info - integer. Return code ! jx - integer(optional). The column offset. ! -function psb_zamax (x,desc_a, info, jx) - use psb_penv_mod - use psb_serial_mod - use psb_descriptor_type - use psb_check_mod - use psb_error_mod +function psb_zamax(x,desc_a, info, jx) result(res) + use psb_base_mod, psb_protect_name => psb_zamax + implicit none - complex(psb_dpk_), intent(in) :: x(:,:) - type(psb_desc_type), intent(in) :: desc_a - integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_), optional, intent(in) :: jx - real(psb_dpk_) :: psb_zamax + complex(psb_dpk_), intent(in) :: x(:,:) + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: jx + real(psb_dpk_) :: res ! locals integer(psb_ipk_) :: ictxt, np, me,& - & err_act, iix, jjx, ix, ijx, m, imax, izamax - real(psb_dpk_) :: amax - character(len=20) :: name, ch_err - double complex :: zdum - double precision :: cabs1 - cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) ) + & err_act, iix, jjx, ix, ijx, m, ldx + character(len=20) :: name, ch_err name='psb_zamax' if(psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - amax=0.d0 - ictxt=desc_a%get_context() + ictxt = desc_a%get_context() call psb_info(ictxt, me, np) if (np == -1) then @@ -91,8 +83,9 @@ function psb_zamax (x,desc_a, info, jx) endif 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 info=psb_err_from_subroutine_ ch_err='psb_chkvect' @@ -108,16 +101,13 @@ function psb_zamax (x,desc_a, info, jx) ! compute local max if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then - imax=izamax(desc_a%get_local_rows()-iix+1,x(iix,jjx),1) - amax=cabs1(x(iix+imax-1,jjx)) + res = psb_amax(desc_a%get_local_rows()-iix+1,x(:,jjx)) else - amax = dzero + res = dzero end if ! compute global max - call psb_amx(ictxt, amax) - - psb_zamax=amax + call psb_amx(ictxt, res) call psb_erractionrestore(err_act) return @@ -133,93 +123,6 @@ function psb_zamax (x,desc_a, info, jx) 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. ! info - integer. Return code ! -function psb_zamaxv (x,desc_a, info) - use psb_penv_mod - use psb_serial_mod - use psb_descriptor_type - use psb_check_mod - use psb_error_mod +function psb_zamaxv (x,desc_a, info) result(res) + use psb_base_mod, psb_protect_name => psb_zamaxv + implicit none complex(psb_dpk_), intent(in) :: x(:) - type(psb_desc_type), intent(in) :: desc_a - integer(psb_ipk_), intent(out) :: info - real(psb_dpk_) :: psb_zamaxv + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + real(psb_dpk_) :: res ! locals integer(psb_ipk_) :: ictxt, np, me,& - & err_act, iix, jjx, jx, ix, m, imax, izamax - real(psb_dpk_) :: amax - complex(psb_dpk_) :: cmax + & err_act, iix, jjx, jx, ix, m, ldx + character(len=20) :: name, ch_err - double complex :: zdum - double precision :: cabs1 - cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) ) name='psb_zamaxv' if(psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - amax=0.d0 ictxt=desc_a%get_context() @@ -307,8 +202,9 @@ function psb_zamaxv (x,desc_a, info) jx = 1 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 info=psb_err_from_subroutine_ ch_err='psb_chkvect' @@ -324,17 +220,13 @@ function psb_zamaxv (x,desc_a, info) ! compute local max if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then - imax=izamax(desc_a%get_local_rows()-iix+1,x(iix),1) - cmax=(x(iix+imax-1)) - amax=cabs1(cmax) + res = psb_amax(desc_a%get_local_rows()-iix+1,x) else - amax = dzero + res = dzero end if ! compute global max - call psb_amx(ictxt, amax) - - psb_zamaxv=amax + call psb_amx(ictxt, res) call psb_erractionrestore(err_act) return @@ -349,6 +241,89 @@ function psb_zamaxv (x,desc_a, info) return 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 !!$ (C) Copyright 2006, 2007, 2008, 2009, 2010 @@ -397,31 +372,26 @@ end function psb_zamaxv ! subroutine psb_zamaxvs(res,x,desc_a, info) use psb_base_mod, psb_protect_name => psb_zamaxvs + implicit none complex(psb_dpk_), intent(in) :: x(:) - type(psb_desc_type), intent(in) :: desc_a - integer(psb_ipk_), intent(out) :: info - real(psb_dpk_), intent(out) :: res + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + real(psb_dpk_), intent(out) :: res ! locals integer(psb_ipk_) :: ictxt, np, me,& - & err_act, iix, jjx, ix, ijx, m, imax, izamax - real(psb_dpk_) :: amax - character(len=20) :: name, ch_err - complex(psb_dpk_) :: cmax - double complex :: zdum - double precision :: cabs1 - cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) ) + & err_act, iix, jjx, ix, ijx, m, ldx + character(len=20) :: name, ch_err name='psb_zamaxvs' if(psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - amax=0.d0 - ictxt=desc_a%get_context() + ictxt = desc_a%get_context() call psb_info(ictxt, me, np) if (np == -1) then @@ -434,7 +404,8 @@ subroutine psb_zamaxvs(res,x,desc_a, info) ijx=1 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 info=psb_err_from_subroutine_ ch_err='psb_chkvect' @@ -450,17 +421,13 @@ subroutine psb_zamaxvs(res,x,desc_a, info) ! compute local max if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then - imax=izamax(desc_a%get_local_rows()-iix+1,x(iix),1) - cmax=(x(iix+imax-1)) - amax=cabs1(cmax) + res = psb_amax(desc_a%get_local_rows()-iix+1,x) else - amax = dzero + res = dzero end if ! compute global max - call psb_amx(ictxt, amax) - - res = amax + call psb_amx(ictxt, res) call psb_erractionrestore(err_act) return @@ -521,31 +488,25 @@ end subroutine psb_zamaxvs ! subroutine psb_zmamaxs(res,x,desc_a, info,jx) use psb_base_mod, psb_protect_name => psb_zmamaxs + implicit none complex(psb_dpk_), intent(in) :: x(:,:) - type(psb_desc_type), intent(in) :: desc_a - integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_), optional, intent(in) :: jx + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: jx real(psb_dpk_), intent(out) :: res(:) ! locals integer(psb_ipk_) :: ictxt, np, me,& - & err_act, iix, jjx, ix, ijx, m, imax, i, k, izamax - real(psb_dpk_) :: amax + & err_act, iix, jjx, ix, ijx, m, ldx, i, k 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' if (psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - amax=0.d0 - ictxt=desc_a%get_context() 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() k = min(size(x,2),size(res,1)) - - 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 info=psb_err_from_subroutine_ ch_err='psb_chkvect' @@ -579,20 +540,17 @@ subroutine psb_zmamaxs(res,x,desc_a, info,jx) goto 9999 end if + res(1:k) = dzero ! compute local max if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then - do i=1,k - imax=izamax(desc_a%get_local_rows()-iix+1,x(iix,jjx+i-1),1) - cmax=(x(iix+imax-1,jjx+i-1)) - res(i)=cabs1(cmax) - end do - else - amax = dzero + do i=1,k + res(i) = psb_amax(desc_a%get_local_rows()-iix+1,x(:,jjx+i-1)) + end do end if ! compute global max call psb_amx(ictxt, res(1:k)) - + call psb_erractionrestore(err_act) return diff --git a/base/psblas/psb_zasum.f90 b/base/psblas/psb_zasum.f90 index 777d9a03..24c2b666 100644 --- a/base/psblas/psb_zasum.f90 +++ b/base/psblas/psb_zasum.f90 @@ -44,37 +44,27 @@ ! info - integer. Return code ! 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 complex(psb_dpk_), intent(in) :: x(:,:) type(psb_desc_type), intent(in) :: desc_a integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: jx - real(psb_dpk_) :: psb_zasum + real(psb_dpk_) :: res ! locals integer(psb_ipk_) :: ictxt, np, me, & - & err_act, iix, jjx, ix, ijx, m, i, idx, ndm - real(psb_dpk_) :: asum, dzasum + & err_act, iix, jjx, ix, ijx, m, i, idx, ndm, ldx 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' if(psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - asum=0.d0 ictxt=desc_a%get_context() @@ -93,9 +83,9 @@ function psb_zasum (x,desc_a, info, jx) endif m = desc_a%get_global_rows() - + ldx = size(x,1) ! 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 info=psb_err_from_subroutine_ ch_err='psb_chkvect' @@ -110,29 +100,21 @@ function psb_zasum (x,desc_a, info, jx) end if ! compute local max - if ((m /= 0)) then - if(desc_a%get_local_rows() > 0) then - 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) + if(desc_a%get_local_rows() > 0) then + res = psb_asum(desc_a%get_local_rows()-iix+1,x(:,jjx)) + + ! 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 + else - asum=0.d0 + res = dzero end if - - - psb_zasum=asum + ! compute global sum + call psb_sum(ictxt, res) call psb_erractionrestore(err_act) return @@ -141,31 +123,25 @@ function psb_zasum (x,desc_a, info, jx) call psb_erractionrestore(err_act) if (err_act == psb_act_abort_) then - call psb_error(ictxt) - return + call psb_error(ictxt) + return end if return end function psb_zasum function psb_zasum_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 + use psb_base_mod, psb_protect_name => psb_zasum_vect 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 + integer(psb_ipk_), intent(out) :: info ! locals integer(psb_ipk_) :: ictxt, np, me,& & err_act, iix, jjx, jx, ix, m, imax - real(psb_dpk_) :: asum character(len=20) :: name, ch_err name='psb_zasumv' @@ -173,7 +149,6 @@ function psb_zasum_vect(x, desc_a, info) result(res) info=psb_success_ call psb_erractionsave(err_act) - asum=0.d0 ictxt=desc_a%get_context() @@ -195,8 +170,7 @@ function psb_zasum_vect(x, desc_a, info) result(res) jx = 1 m = desc_a%get_global_rows() - - call psb_chkvect(m,1,x%get_nrows(),ix,jx,desc_a,info,iix,jjx) + 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' @@ -212,15 +186,13 @@ function psb_zasum_vect(x, desc_a, info) result(res) ! compute local max 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 - asum = dzero + res = dzero end if ! compute global sum - call psb_sum(ictxt, asum) - - res=asum + call psb_sum(ictxt, res) call psb_erractionrestore(err_act) return @@ -233,7 +205,7 @@ function psb_zasum_vect(x, desc_a, info) result(res) return end if return - + end function psb_zasum_vect @@ -280,37 +252,26 @@ end function psb_zasum_vect ! desc_a - type(psb_desc_type). The communication descriptor. ! 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 complex(psb_dpk_), intent(in) :: x(:) - type(psb_desc_type), intent(in) :: desc_a - integer(psb_ipk_), intent(out) :: info - real(psb_dpk_) :: psb_zasumv + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + real(psb_dpk_) :: res ! locals integer(psb_ipk_) :: ictxt, np, me,& - & err_act, iix, jjx, jx, ix, m, i, idx, ndm - real(psb_dpk_) :: asum, dzasum + & err_act, iix, jjx, jx, ix, m, i, idx, ndm, ldx 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' if(psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - asum=0.d0 - ictxt=desc_a%get_context() call psb_info(ictxt, me, np) @@ -324,9 +285,9 @@ function psb_zasumv(x,desc_a, info) jx=1 m = desc_a%get_global_rows() - + ldx = size(x,1) ! 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 info=psb_err_from_subroutine_ ch_err='psb_chkvect' @@ -341,30 +302,22 @@ function psb_zasumv(x,desc_a, info) end if ! compute local max - if ((m /= 0)) then - if(desc_a%get_local_rows() > 0) then - asum=dzasum(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))*cabs1(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 + if(desc_a%get_local_rows() > 0) then + 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 - asum=0.d0 + res = dzero end if - psb_zasumv=asum + ! compute global sum + call psb_sum(ictxt, res) call psb_erractionrestore(err_act) return @@ -426,30 +379,24 @@ end function psb_zasumv ! subroutine psb_zasumvs(res,x,desc_a, info) use psb_base_mod, psb_protect_name => psb_zasumvs + implicit none - complex(psb_dpk_), intent(in) :: x(:) - real(psb_dpk_), intent(out) :: res - type(psb_desc_type), intent(in) :: desc_a - integer(psb_ipk_), intent(out) :: info + complex(psb_dpk_), intent(in) :: x(:) + real(psb_dpk_), intent(out) :: res + 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, ix, jx, m, i, idx, ndm - real(psb_dpk_) :: asum, dzasum + & err_act, iix, jjx, ix, jx, m, i, idx, ndm, ldx 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' if(psb_get_errstatus() /= 0) return info=psb_success_ call psb_erractionsave(err_act) - asum=0.d0 - ictxt=desc_a%get_context() call psb_info(ictxt, me, np) @@ -463,9 +410,9 @@ subroutine psb_zasumvs(res,x,desc_a, info) jx = 1 m = desc_a%get_global_rows() - + ldx = size(x,1) ! 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 info=psb_err_from_subroutine_ ch_err='psb_chkvect' @@ -480,31 +427,22 @@ subroutine psb_zasumvs(res,x,desc_a, info) end if ! compute local max - if ((m /= 0)) then - if(desc_a%get_local_rows() > 0) then - asum=dzasum(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))*cabs1(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 + if(desc_a%get_local_rows() > 0) then + 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 - asum=0.d0 + res = dzero end if - - res = asum + ! compute global sum + call psb_sum(ictxt,res) call psb_erractionrestore(err_act) return diff --git a/base/serial/Makefile b/base/serial/Makefile index 3600a49f..764d481c 100644 --- a/base/serial/Makefile +++ b/base/serial/Makefile @@ -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_sgeprt.o psb_dgeprt.o psb_cgeprt.o psb_zgeprt.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=.. MODDIR=../modules