From 97bc2e3134d65f2f82e3086ffff7ee9d7b1bbb42 Mon Sep 17 00:00:00 2001 From: Weslley da Silva Pereira Date: Thu, 27 Mar 2025 15:42:03 -0600 Subject: [PATCH 1/5] Adds new icamax with exception handling --- BLAS/SRC/CMakeLists.txt | 4 +- BLAS/SRC/{ => DEPRECATED}/icamax.f | 0 BLAS/SRC/icamax.f90 | 188 +++++++++++++++++++++++++++++ 3 files changed, 190 insertions(+), 2 deletions(-) rename BLAS/SRC/{ => DEPRECATED}/icamax.f (100%) create mode 100644 BLAS/SRC/icamax.f90 diff --git a/BLAS/SRC/CMakeLists.txt b/BLAS/SRC/CMakeLists.txt index b9e6f7c4a5..bb1aefc6f8 100644 --- a/BLAS/SRC/CMakeLists.txt +++ b/BLAS/SRC/CMakeLists.txt @@ -32,7 +32,7 @@ set(SBLAS1 isamax.f sasum.f saxpy.f scopy.f sdot.f snrm2.f90 srot.f srotg.f90 sscal.f sswap.f sdsdot.f srotmg.f srotm.f) -set(CBLAS1 scabs1.f scasum.f scnrm2.f90 icamax.f caxpy.f ccopy.f +set(CBLAS1 scabs1.f scasum.f scnrm2.f90 icamax.f90 caxpy.f ccopy.f cdotc.f cdotu.f csscal.f crotg.f90 cscal.f cswap.f csrot.f) set(DBLAS1 idamax.f dasum.f daxpy.f dcopy.f ddot.f dnrm2.f90 @@ -49,7 +49,7 @@ set(CB1AUX sswap.f) set(ZB1AUX - icamax.f idamax.f + icamax.f90 idamax.f cgemm.f cherk.f cscal.f ctrsm.f dasum.f daxpy.f dcopy.f ddot.f dgemm.f dgemv.f dnrm2.f90 drot.f dscal.f dswap.f diff --git a/BLAS/SRC/icamax.f b/BLAS/SRC/DEPRECATED/icamax.f similarity index 100% rename from BLAS/SRC/icamax.f rename to BLAS/SRC/DEPRECATED/icamax.f diff --git a/BLAS/SRC/icamax.f90 b/BLAS/SRC/icamax.f90 new file mode 100644 index 0000000000..8fccc26f37 --- /dev/null +++ b/BLAS/SRC/icamax.f90 @@ -0,0 +1,188 @@ +!> \brief \b ICAMAX +! +! =========== DOCUMENTATION =========== +! +! Online html documentation available at +! http://www.netlib.org/lapack/explore-html/ +! +! Definition: +! =========== +! +! INTEGER FUNCTION ICAMAX(N,CX,INCX) +! +! .. Scalar Arguments .. +! INTEGER INCX,N +! .. +! .. Array Arguments .. +! COMPLEX CX(*) +! .. +! +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> ICAMAX finds the index of the first element having maximum |Re(.)| + |Im(.)| +!> \endverbatim +! +! Arguments: +! ========== +! +!> \param[in] N +!> \verbatim +!> N is INTEGER +!> number of elements in input vector(s) +!> \endverbatim +!> +!> \param[in] CX +!> \verbatim +!> CX is COMPLEX array, dimension ( 1 + ( N - 1 )*abs( INCX ) ) +!> \endverbatim +!> +!> \param[in] INCX +!> \verbatim +!> INCX is INTEGER +!> storage spacing between elements of CX +!> \endverbatim +! +! Authors: +! ======== +! +!> James Demmel, University of California Berkeley, USA +!> Weslley Pereira, University of Colorado Denver, USA +! +!> \ingroup iamax +! +!> \par Further Details: +! ===================== +!> +!> \verbatim +!> +!> James Demmel et al. Proposed Consistent Exception Handling for the BLAS and +!> LAPACK, 2022 (https://arxiv.org/abs/2207.09281). +!> +!> \endverbatim +!> +! ===================================================================== +integer function icamax(n, x, incx) + integer, parameter :: wp = kind(1.e0) +! +! -- Reference BLAS level1 routine -- +! -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +! +! .. Constants .. + real(wp), parameter :: hugeval = huge(0.0_wp) +! +! .. Scalar Arguments .. + integer :: n, incx +! +! .. Array Arguments .. + complex(wp) :: x(*) +! .. +! .. Local Scalars .. + integer :: i, j, ix, jx + real(wp) :: val, smax, scaledsmax +! +! Quick return if possible +! + icamax = 0 + if (n < 1 .or. incx < 1) return +! + icamax = 1 + if (n == 1) return +! + icamax = 0 + scaledsmax = 0 + smax = -1 +! +! scaledsmax = 1 indicates that x(i) finite but +! abs(real(x(i))) + abs(imag(x(i))) is not finite +! + if (incx == 1) then + ! code for increment equal to 1 + do i = 1, n + if (isnan(real(x(i))) .or. isnan(imag(x(i)))) then + ! return when first NaN found + icamax = i + return + elseif (abs(real(x(i))) > hugeval .or. abs(imag(x(i))) > hugeval) then + ! keep looking for first NaN + do j = i+1, n + if (isnan(real(x(j))) .or. isnan(imag(x(j)))) then + ! return when first NaN found + icamax = j + return + endif + enddo + ! record location of first Inf + icamax = i + return + else ! still no Inf found yet + if (scaledsmax == 0) then + ! no abs(real(x(i))) + abs(imag(x(i))) = Inf yet + val = abs(real(x(i))) + abs(imag(x(i))) + if (abs(val) > hugeval) then + scaledsmax = 1 + smax = 0.25*abs(real(x(i))) + 0.25*abs(imag(x(i))) + icamax = i + elseif (val > smax) then ! everything finite so far + smax = val + icamax = i + endif + else ! scaledsmax = 1 + val = 0.25*abs(real(x(i))) + 0.25*abs(imag(x(i))) + if (val > smax) then + smax = val + icamax = i + endif + endif + endif + end do + else + ! code for increment not equal to 1 + ix = 1 + do i = 1, n + if (isnan(real(x(ix))) .or. isnan(imag(x(ix)))) then + ! return when first NaN found + icamax = i + return + elseif (abs(real(x(ix))) > hugeval .or. abs(imag(x(ix))) > hugeval) then + ! keep looking for first NaN + jx = ix + incx + do j = i+1, n + if (isnan(real(x(jx))) .or. isnan(imag(x(jx)))) then + ! return when first NaN found + icamax = j + return + endif + jx = jx + incx + enddo + ! record location of first Inf + icamax = i + return + else ! still no Inf found yet + if (scaledsmax == 0) then + ! no abs(real(x(ix))) + abs(imag(x(ix))) = Inf yet + val = abs(real(x(ix))) + abs(imag(x(ix))) + if (abs(val) > hugeval) then + scaledsmax = 1 + smax = 0.25*abs(real(x(ix))) + 0.25*abs(imag(x(ix))) + icamax = i + elseif (val > smax) then ! everything finite so far + smax = val + icamax = i + endif + else ! scaledsmax = 1 + val = 0.25*abs(real(x(ix))) + 0.25*abs(imag(x(ix))) + if (val > smax) then + smax = val + icamax = i + endif + endif + endif + ix = ix + incx + end do + endif +end From 50dc1bf23804ea7d36d94fbdb4bd0096afcd3d90 Mon Sep 17 00:00:00 2001 From: Weslley da Silva Pereira Date: Thu, 27 Mar 2025 16:12:17 -0600 Subject: [PATCH 2/5] Adds new izamax with exception handling --- BLAS/SRC/CMakeLists.txt | 2 +- BLAS/SRC/{ => DEPRECATED}/izamax.f | 0 BLAS/SRC/icamax.f90 | 10 +- BLAS/SRC/izamax.f90 | 188 +++++++++++++++++++++++++++++ 4 files changed, 194 insertions(+), 6 deletions(-) rename BLAS/SRC/{ => DEPRECATED}/izamax.f (100%) create mode 100644 BLAS/SRC/izamax.f90 diff --git a/BLAS/SRC/CMakeLists.txt b/BLAS/SRC/CMakeLists.txt index bb1aefc6f8..f41bef6d61 100644 --- a/BLAS/SRC/CMakeLists.txt +++ b/BLAS/SRC/CMakeLists.txt @@ -40,7 +40,7 @@ set(DBLAS1 idamax.f dasum.f daxpy.f dcopy.f ddot.f dnrm2.f90 set(DB1AUX sscal.f isamax.f) -set(ZBLAS1 dcabs1.f dzasum.f dznrm2.f90 izamax.f zaxpy.f zcopy.f +set(ZBLAS1 dcabs1.f dzasum.f dznrm2.f90 izamax.f90 zaxpy.f zcopy.f zdotc.f zdotu.f zdscal.f zrotg.f90 zscal.f zswap.f zdrot.f) set(CB1AUX diff --git a/BLAS/SRC/izamax.f b/BLAS/SRC/DEPRECATED/izamax.f similarity index 100% rename from BLAS/SRC/izamax.f rename to BLAS/SRC/DEPRECATED/izamax.f diff --git a/BLAS/SRC/icamax.f90 b/BLAS/SRC/icamax.f90 index 8fccc26f37..192f5ccb55 100644 --- a/BLAS/SRC/icamax.f90 +++ b/BLAS/SRC/icamax.f90 @@ -8,13 +8,13 @@ ! Definition: ! =========== ! -! INTEGER FUNCTION ICAMAX(N,CX,INCX) +! INTEGER FUNCTION ICAMAX(N,X,INCX) ! ! .. Scalar Arguments .. ! INTEGER INCX,N ! .. ! .. Array Arguments .. -! COMPLEX CX(*) +! COMPLEX X(*) ! .. ! ! @@ -35,15 +35,15 @@ !> number of elements in input vector(s) !> \endverbatim !> -!> \param[in] CX +!> \param[in] X !> \verbatim -!> CX is COMPLEX array, dimension ( 1 + ( N - 1 )*abs( INCX ) ) +!> X is COMPLEX array, dimension ( 1 + ( N - 1 )*abs( INCX ) ) !> \endverbatim !> !> \param[in] INCX !> \verbatim !> INCX is INTEGER -!> storage spacing between elements of CX +!> storage spacing between elements of X !> \endverbatim ! ! Authors: diff --git a/BLAS/SRC/izamax.f90 b/BLAS/SRC/izamax.f90 new file mode 100644 index 0000000000..c36e623ea6 --- /dev/null +++ b/BLAS/SRC/izamax.f90 @@ -0,0 +1,188 @@ +!> \brief \b IZAMAX +! +! =========== DOCUMENTATION =========== +! +! Online html documentation available at +! http://www.netlib.org/lapack/explore-html/ +! +! Definition: +! =========== +! +! INTEGER FUNCTION IZAMAX(N,X,INCX) +! +! .. Scalar Arguments .. +! INTEGER INCX,N +! .. +! .. Array Arguments .. +! DOUBLE COMPLEX X(*) +! .. +! +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> IZAMAX finds the index of the first element having maximum |Re(.)| + |Im(.)| +!> \endverbatim +! +! Arguments: +! ========== +! +!> \param[in] N +!> \verbatim +!> N is INTEGER +!> number of elements in input vector(s) +!> \endverbatim +!> +!> \param[in] X +!> \verbatim +!> X is DOUBLE COMPLEX array, dimension ( 1 + ( N - 1 )*abs( INCX ) ) +!> \endverbatim +!> +!> \param[in] INCX +!> \verbatim +!> INCX is INTEGER +!> storage spacing between elements of X +!> \endverbatim +! +! Authors: +! ======== +! +!> James Demmel, University of California Berkeley, USA +!> Weslley Pereira, University of Colorado Denver, USA +! +!> \ingroup iamax +! +!> \par Further Details: +! ===================== +!> +!> \verbatim +!> +!> James Demmel et al. Proposed Consistent Exception Handling for the BLAS and +!> LAPACK, 2022 (https://arxiv.org/abs/2207.09281). +!> +!> \endverbatim +!> +! ===================================================================== +integer function izamax(n, x, incx) + integer, parameter :: wp = kind(1.d0) +! +! -- Reference BLAS level1 routine -- +! -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +! +! .. Constants .. + real(wp), parameter :: hugeval = huge(0.0_wp) +! +! .. Scalar Arguments .. + integer :: n, incx +! +! .. Array Arguments .. + complex(wp) :: x(*) +! .. +! .. Local Scalars .. + integer :: i, j, ix, jx + real(wp) :: val, smax, scaledsmax +! +! Quick return if possible +! + izamax = 0 + if (n < 1 .or. incx < 1) return +! + izamax = 1 + if (n == 1) return +! + izamax = 0 + scaledsmax = 0 + smax = -1 +! +! scaledsmax = 1 indicates that x(i) finite but +! abs(real(x(i))) + abs(imag(x(i))) is not finite +! + if (incx == 1) then + ! code for increment equal to 1 + do i = 1, n + if (isnan(real(x(i))) .or. isnan(imag(x(i)))) then + ! return when first NaN found + izamax = i + return + elseif (abs(real(x(i))) > hugeval .or. abs(imag(x(i))) > hugeval) then + ! keep looking for first NaN + do j = i+1, n + if (isnan(real(x(j))) .or. isnan(imag(x(j)))) then + ! return when first NaN found + izamax = j + return + endif + enddo + ! record location of first Inf + izamax = i + return + else ! still no Inf found yet + if (scaledsmax == 0) then + ! no abs(real(x(i))) + abs(imag(x(i))) = Inf yet + val = abs(real(x(i))) + abs(imag(x(i))) + if (abs(val) > hugeval) then + scaledsmax = 1 + smax = 0.25*abs(real(x(i))) + 0.25*abs(imag(x(i))) + izamax = i + elseif (val > smax) then ! everything finite so far + smax = val + izamax = i + endif + else ! scaledsmax = 1 + val = 0.25*abs(real(x(i))) + 0.25*abs(imag(x(i))) + if (val > smax) then + smax = val + izamax = i + endif + endif + endif + end do + else + ! code for increment not equal to 1 + ix = 1 + do i = 1, n + if (isnan(real(x(ix))) .or. isnan(imag(x(ix)))) then + ! return when first NaN found + izamax = i + return + elseif (abs(real(x(ix))) > hugeval .or. abs(imag(x(ix))) > hugeval) then + ! keep looking for first NaN + jx = ix + incx + do j = i+1, n + if (isnan(real(x(jx))) .or. isnan(imag(x(jx)))) then + ! return when first NaN found + izamax = j + return + endif + jx = jx + incx + enddo + ! record location of first Inf + izamax = i + return + else ! still no Inf found yet + if (scaledsmax == 0) then + ! no abs(real(x(ix))) + abs(imag(x(ix))) = Inf yet + val = abs(real(x(ix))) + abs(imag(x(ix))) + if (abs(val) > hugeval) then + scaledsmax = 1 + smax = 0.25*abs(real(x(ix))) + 0.25*abs(imag(x(ix))) + izamax = i + elseif (val > smax) then ! everything finite so far + smax = val + izamax = i + endif + else ! scaledsmax = 1 + val = 0.25*abs(real(x(ix))) + 0.25*abs(imag(x(ix))) + if (val > smax) then + smax = val + izamax = i + endif + endif + endif + ix = ix + incx + end do + endif +end From e3a48809c816cb5bca2ccf1d74071bdca51704c5 Mon Sep 17 00:00:00 2001 From: Weslley da Silva Pereira Date: Fri, 28 Mar 2025 08:23:58 -0600 Subject: [PATCH 3/5] use `/=` instead of `isnan` --- BLAS/SRC/icamax.f90 | 8 ++++---- BLAS/SRC/izamax.f90 | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/BLAS/SRC/icamax.f90 b/BLAS/SRC/icamax.f90 index 192f5ccb55..1b658f74b5 100644 --- a/BLAS/SRC/icamax.f90 +++ b/BLAS/SRC/icamax.f90 @@ -103,14 +103,14 @@ integer function icamax(n, x, incx) if (incx == 1) then ! code for increment equal to 1 do i = 1, n - if (isnan(real(x(i))) .or. isnan(imag(x(i)))) then + if (x(i) /= x(i)) then ! return when first NaN found icamax = i return elseif (abs(real(x(i))) > hugeval .or. abs(imag(x(i))) > hugeval) then ! keep looking for first NaN do j = i+1, n - if (isnan(real(x(j))) .or. isnan(imag(x(j)))) then + if (x(j) /= x(j)) then ! return when first NaN found icamax = j return @@ -144,7 +144,7 @@ integer function icamax(n, x, incx) ! code for increment not equal to 1 ix = 1 do i = 1, n - if (isnan(real(x(ix))) .or. isnan(imag(x(ix)))) then + if (x(ix) /= x(ix)) then ! return when first NaN found icamax = i return @@ -152,7 +152,7 @@ integer function icamax(n, x, incx) ! keep looking for first NaN jx = ix + incx do j = i+1, n - if (isnan(real(x(jx))) .or. isnan(imag(x(jx)))) then + if (x(jx) /= x(jx)) then ! return when first NaN found icamax = j return diff --git a/BLAS/SRC/izamax.f90 b/BLAS/SRC/izamax.f90 index c36e623ea6..2ce8541172 100644 --- a/BLAS/SRC/izamax.f90 +++ b/BLAS/SRC/izamax.f90 @@ -103,14 +103,14 @@ integer function izamax(n, x, incx) if (incx == 1) then ! code for increment equal to 1 do i = 1, n - if (isnan(real(x(i))) .or. isnan(imag(x(i)))) then + if (x(i) /= x(i)) then ! return when first NaN found izamax = i return elseif (abs(real(x(i))) > hugeval .or. abs(imag(x(i))) > hugeval) then ! keep looking for first NaN do j = i+1, n - if (isnan(real(x(j))) .or. isnan(imag(x(j)))) then + if (x(j) /= x(j)) then ! return when first NaN found izamax = j return @@ -144,7 +144,7 @@ integer function izamax(n, x, incx) ! code for increment not equal to 1 ix = 1 do i = 1, n - if (isnan(real(x(ix))) .or. isnan(imag(x(ix)))) then + if (x(ix) /= x(ix)) then ! return when first NaN found izamax = i return @@ -152,7 +152,7 @@ integer function izamax(n, x, incx) ! keep looking for first NaN jx = ix + incx do j = i+1, n - if (isnan(real(x(jx))) .or. isnan(imag(x(jx)))) then + if (x(jx) /= x(jx)) then ! return when first NaN found izamax = j return From e338bb87e11e08e1462f5a9468e17150bff4a94e Mon Sep 17 00:00:00 2001 From: Weslley da Silva Pereira Date: Fri, 28 Mar 2025 09:24:58 -0600 Subject: [PATCH 4/5] Change my affiliation in icamax and izamax --- BLAS/SRC/icamax.f90 | 2 +- BLAS/SRC/izamax.f90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/BLAS/SRC/icamax.f90 b/BLAS/SRC/icamax.f90 index 1b658f74b5..32c9cbdc8f 100644 --- a/BLAS/SRC/icamax.f90 +++ b/BLAS/SRC/icamax.f90 @@ -50,7 +50,7 @@ ! ======== ! !> James Demmel, University of California Berkeley, USA -!> Weslley Pereira, University of Colorado Denver, USA +!> Weslley Pereira, National Renewable Energy Laboratory, USA ! !> \ingroup iamax ! diff --git a/BLAS/SRC/izamax.f90 b/BLAS/SRC/izamax.f90 index 2ce8541172..9a906698b8 100644 --- a/BLAS/SRC/izamax.f90 +++ b/BLAS/SRC/izamax.f90 @@ -50,7 +50,7 @@ ! ======== ! !> James Demmel, University of California Berkeley, USA -!> Weslley Pereira, University of Colorado Denver, USA +!> Weslley Pereira, National Renewable Energy Laboratory, USA ! !> \ingroup iamax ! From d70c8c040daceefadb7481418c8a0899d04e1f22 Mon Sep 17 00:00:00 2001 From: Weslley da Silva Pereira Date: Mon, 31 Mar 2025 13:32:24 -0600 Subject: [PATCH 5/5] Improvements following @angsch's review --- BLAS/SRC/icamax.f90 | 25 +++++++++++++------------ BLAS/SRC/izamax.f90 | 25 +++++++++++++------------ 2 files changed, 26 insertions(+), 24 deletions(-) diff --git a/BLAS/SRC/icamax.f90 b/BLAS/SRC/icamax.f90 index 32c9cbdc8f..cf19199f96 100644 --- a/BLAS/SRC/icamax.f90 +++ b/BLAS/SRC/icamax.f90 @@ -83,7 +83,8 @@ integer function icamax(n, x, incx) ! .. ! .. Local Scalars .. integer :: i, j, ix, jx - real(wp) :: val, smax, scaledsmax + real(wp) :: val, smax + logical :: scaledsmax ! ! Quick return if possible ! @@ -94,11 +95,11 @@ integer function icamax(n, x, incx) if (n == 1) return ! icamax = 0 - scaledsmax = 0 + scaledsmax = .false. smax = -1 ! -! scaledsmax = 1 indicates that x(i) finite but -! abs(real(x(i))) + abs(imag(x(i))) is not finite +! scaledsmax = .true. indicates that x(icamax) is finite but +! abs(real(x(icamax))) + abs(imag(x(icamax))) overflows ! if (incx == 1) then ! code for increment equal to 1 @@ -120,18 +121,18 @@ integer function icamax(n, x, incx) icamax = i return else ! still no Inf found yet - if (scaledsmax == 0) then + if (.not. scaledsmax) then ! no abs(real(x(i))) + abs(imag(x(i))) = Inf yet val = abs(real(x(i))) + abs(imag(x(i))) - if (abs(val) > hugeval) then - scaledsmax = 1 + if (val > hugeval) then + scaledsmax = .true. smax = 0.25*abs(real(x(i))) + 0.25*abs(imag(x(i))) icamax = i elseif (val > smax) then ! everything finite so far smax = val icamax = i endif - else ! scaledsmax = 1 + else ! scaledsmax val = 0.25*abs(real(x(i))) + 0.25*abs(imag(x(i))) if (val > smax) then smax = val @@ -163,18 +164,18 @@ integer function icamax(n, x, incx) icamax = i return else ! still no Inf found yet - if (scaledsmax == 0) then + if (.not. scaledsmax) then ! no abs(real(x(ix))) + abs(imag(x(ix))) = Inf yet val = abs(real(x(ix))) + abs(imag(x(ix))) - if (abs(val) > hugeval) then - scaledsmax = 1 + if (val > hugeval) then + scaledsmax = .true. smax = 0.25*abs(real(x(ix))) + 0.25*abs(imag(x(ix))) icamax = i elseif (val > smax) then ! everything finite so far smax = val icamax = i endif - else ! scaledsmax = 1 + else ! scaledsmax val = 0.25*abs(real(x(ix))) + 0.25*abs(imag(x(ix))) if (val > smax) then smax = val diff --git a/BLAS/SRC/izamax.f90 b/BLAS/SRC/izamax.f90 index 9a906698b8..6bbb5a1686 100644 --- a/BLAS/SRC/izamax.f90 +++ b/BLAS/SRC/izamax.f90 @@ -83,7 +83,8 @@ integer function izamax(n, x, incx) ! .. ! .. Local Scalars .. integer :: i, j, ix, jx - real(wp) :: val, smax, scaledsmax + real(wp) :: val, smax + logical :: scaledsmax ! ! Quick return if possible ! @@ -94,11 +95,11 @@ integer function izamax(n, x, incx) if (n == 1) return ! izamax = 0 - scaledsmax = 0 + scaledsmax = .false. smax = -1 ! -! scaledsmax = 1 indicates that x(i) finite but -! abs(real(x(i))) + abs(imag(x(i))) is not finite +! scaledsmax = .true. indicates that x(izamax) is finite but +! abs(real(x(izamax))) + abs(imag(x(izamax))) overflows ! if (incx == 1) then ! code for increment equal to 1 @@ -120,18 +121,18 @@ integer function izamax(n, x, incx) izamax = i return else ! still no Inf found yet - if (scaledsmax == 0) then + if (.not. scaledsmax) then ! no abs(real(x(i))) + abs(imag(x(i))) = Inf yet val = abs(real(x(i))) + abs(imag(x(i))) - if (abs(val) > hugeval) then - scaledsmax = 1 + if (val > hugeval) then + scaledsmax = .true. smax = 0.25*abs(real(x(i))) + 0.25*abs(imag(x(i))) izamax = i elseif (val > smax) then ! everything finite so far smax = val izamax = i endif - else ! scaledsmax = 1 + else ! scaledsmax val = 0.25*abs(real(x(i))) + 0.25*abs(imag(x(i))) if (val > smax) then smax = val @@ -163,18 +164,18 @@ integer function izamax(n, x, incx) izamax = i return else ! still no Inf found yet - if (scaledsmax == 0) then + if (.not. scaledsmax) then ! no abs(real(x(ix))) + abs(imag(x(ix))) = Inf yet val = abs(real(x(ix))) + abs(imag(x(ix))) - if (abs(val) > hugeval) then - scaledsmax = 1 + if (val > hugeval) then + scaledsmax = .true. smax = 0.25*abs(real(x(ix))) + 0.25*abs(imag(x(ix))) izamax = i elseif (val > smax) then ! everything finite so far smax = val izamax = i endif - else ! scaledsmax = 1 + else ! scaledsmax val = 0.25*abs(real(x(ix))) + 0.25*abs(imag(x(ix))) if (val > smax) then smax = val