Skip to content

Commit 97bc2e3

Browse files
author
Weslley da Silva Pereira
committed
Adds new icamax with exception handling
1 parent ca82e5e commit 97bc2e3

File tree

3 files changed

+190
-2
lines changed

3 files changed

+190
-2
lines changed

‎BLAS/SRC/CMakeLists.txt

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -32,7 +32,7 @@
3232
set(SBLAS1 isamax.f sasum.f saxpy.f scopy.f sdot.f snrm2.f90
3333
srot.f srotg.f90 sscal.f sswap.f sdsdot.f srotmg.f srotm.f)
3434

35-
set(CBLAS1 scabs1.f scasum.f scnrm2.f90 icamax.f caxpy.f ccopy.f
35+
set(CBLAS1 scabs1.f scasum.f scnrm2.f90 icamax.f90 caxpy.f ccopy.f
3636
cdotc.f cdotu.f csscal.f crotg.f90 cscal.f cswap.f csrot.f)
3737

3838
set(DBLAS1 idamax.f dasum.f daxpy.f dcopy.f ddot.f dnrm2.f90
@@ -49,7 +49,7 @@ set(CB1AUX
4949
sswap.f)
5050

5151
set(ZB1AUX
52-
icamax.f idamax.f
52+
icamax.f90 idamax.f
5353
cgemm.f cherk.f cscal.f ctrsm.f
5454
dasum.f daxpy.f dcopy.f ddot.f dgemm.f dgemv.f dnrm2.f90 drot.f dscal.f
5555
dswap.f

‎BLAS/SRC/icamax.f renamed to ‎BLAS/SRC/DEPRECATED/icamax.f

File renamed without changes.

‎BLAS/SRC/icamax.f90

Lines changed: 188 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,188 @@
1+
!> \brief \b ICAMAX
2+
!
3+
! =========== DOCUMENTATION ===========
4+
!
5+
! Online html documentation available at
6+
! http://www.netlib.org/lapack/explore-html/
7+
!
8+
! Definition:
9+
! ===========
10+
!
11+
! INTEGER FUNCTION ICAMAX(N,CX,INCX)
12+
!
13+
! .. Scalar Arguments ..
14+
! INTEGER INCX,N
15+
! ..
16+
! .. Array Arguments ..
17+
! COMPLEX CX(*)
18+
! ..
19+
!
20+
!
21+
!> \par Purpose:
22+
! =============
23+
!>
24+
!> \verbatim
25+
!>
26+
!> ICAMAX finds the index of the first element having maximum |Re(.)| + |Im(.)|
27+
!> \endverbatim
28+
!
29+
! Arguments:
30+
! ==========
31+
!
32+
!> \param[in] N
33+
!> \verbatim
34+
!> N is INTEGER
35+
!> number of elements in input vector(s)
36+
!> \endverbatim
37+
!>
38+
!> \param[in] CX
39+
!> \verbatim
40+
!> CX is COMPLEX array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
41+
!> \endverbatim
42+
!>
43+
!> \param[in] INCX
44+
!> \verbatim
45+
!> INCX is INTEGER
46+
!> storage spacing between elements of CX
47+
!> \endverbatim
48+
!
49+
! Authors:
50+
! ========
51+
!
52+
!> James Demmel, University of California Berkeley, USA
53+
!> Weslley Pereira, University of Colorado Denver, USA
54+
!
55+
!> \ingroup iamax
56+
!
57+
!> \par Further Details:
58+
! =====================
59+
!>
60+
!> \verbatim
61+
!>
62+
!> James Demmel et al. Proposed Consistent Exception Handling for the BLAS and
63+
!> LAPACK, 2022 (https://arxiv.org/abs/2207.09281).
64+
!>
65+
!> \endverbatim
66+
!>
67+
! =====================================================================
68+
integer function icamax(n, x, incx)
69+
integer, parameter :: wp = kind(1.e0)
70+
!
71+
! -- Reference BLAS level1 routine --
72+
! -- Reference BLAS is a software package provided by Univ. of Tennessee, --
73+
! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
74+
!
75+
! .. Constants ..
76+
real(wp), parameter :: hugeval = huge(0.0_wp)
77+
!
78+
! .. Scalar Arguments ..
79+
integer :: n, incx
80+
!
81+
! .. Array Arguments ..
82+
complex(wp) :: x(*)
83+
! ..
84+
! .. Local Scalars ..
85+
integer :: i, j, ix, jx
86+
real(wp) :: val, smax, scaledsmax
87+
!
88+
! Quick return if possible
89+
!
90+
icamax = 0
91+
if (n < 1 .or. incx < 1) return
92+
!
93+
icamax = 1
94+
if (n == 1) return
95+
!
96+
icamax = 0
97+
scaledsmax = 0
98+
smax = -1
99+
!
100+
! scaledsmax = 1 indicates that x(i) finite but
101+
! abs(real(x(i))) + abs(imag(x(i))) is not finite
102+
!
103+
if (incx == 1) then
104+
! code for increment equal to 1
105+
do i = 1, n
106+
if (isnan(real(x(i))) .or. isnan(imag(x(i)))) then
107+
! return when first NaN found
108+
icamax = i
109+
return
110+
elseif (abs(real(x(i))) > hugeval .or. abs(imag(x(i))) > hugeval) then
111+
! keep looking for first NaN
112+
do j = i+1, n
113+
if (isnan(real(x(j))) .or. isnan(imag(x(j)))) then
114+
! return when first NaN found
115+
icamax = j
116+
return
117+
endif
118+
enddo
119+
! record location of first Inf
120+
icamax = i
121+
return
122+
else ! still no Inf found yet
123+
if (scaledsmax == 0) then
124+
! no abs(real(x(i))) + abs(imag(x(i))) = Inf yet
125+
val = abs(real(x(i))) + abs(imag(x(i)))
126+
if (abs(val) > hugeval) then
127+
scaledsmax = 1
128+
smax = 0.25*abs(real(x(i))) + 0.25*abs(imag(x(i)))
129+
icamax = i
130+
elseif (val > smax) then ! everything finite so far
131+
smax = val
132+
icamax = i
133+
endif
134+
else ! scaledsmax = 1
135+
val = 0.25*abs(real(x(i))) + 0.25*abs(imag(x(i)))
136+
if (val > smax) then
137+
smax = val
138+
icamax = i
139+
endif
140+
endif
141+
endif
142+
end do
143+
else
144+
! code for increment not equal to 1
145+
ix = 1
146+
do i = 1, n
147+
if (isnan(real(x(ix))) .or. isnan(imag(x(ix)))) then
148+
! return when first NaN found
149+
icamax = i
150+
return
151+
elseif (abs(real(x(ix))) > hugeval .or. abs(imag(x(ix))) > hugeval) then
152+
! keep looking for first NaN
153+
jx = ix + incx
154+
do j = i+1, n
155+
if (isnan(real(x(jx))) .or. isnan(imag(x(jx)))) then
156+
! return when first NaN found
157+
icamax = j
158+
return
159+
endif
160+
jx = jx + incx
161+
enddo
162+
! record location of first Inf
163+
icamax = i
164+
return
165+
else ! still no Inf found yet
166+
if (scaledsmax == 0) then
167+
! no abs(real(x(ix))) + abs(imag(x(ix))) = Inf yet
168+
val = abs(real(x(ix))) + abs(imag(x(ix)))
169+
if (abs(val) > hugeval) then
170+
scaledsmax = 1
171+
smax = 0.25*abs(real(x(ix))) + 0.25*abs(imag(x(ix)))
172+
icamax = i
173+
elseif (val > smax) then ! everything finite so far
174+
smax = val
175+
icamax = i
176+
endif
177+
else ! scaledsmax = 1
178+
val = 0.25*abs(real(x(ix))) + 0.25*abs(imag(x(ix)))
179+
if (val > smax) then
180+
smax = val
181+
icamax = i
182+
endif
183+
endif
184+
endif
185+
ix = ix + incx
186+
end do
187+
endif
188+
end

0 commit comments

Comments
 (0)