Adds new icamax with exception handling

BLAS/SRC/CMakeLists.txt

@@ -32,15 +32,15 @@
 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
 	drot.f drotg.f90 dscal.f dsdot.f dswap.f drotmg.f drotm.f)
 
 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
@@ -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

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,X,INCX)
+!
+!       .. Scalar Arguments ..
+!       INTEGER INCX,N
+!       ..
+!       .. Array Arguments ..
+!       COMPLEX X(*)
+!       ..
+!
+!
+!> \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] X
+!> \verbatim
+!>          X is 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, National Renewable Energy Laboratory, 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 (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 (x(j) /= 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 (x(ix) /= 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 (x(jx) /= 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