1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
|
// Copyright (c) 2017-2023, University of Tennessee. All rights reserved.
// SPDX-License-Identifier: BSD-3-Clause
// This program is free software: you can redistribute it and/or modify it under
// the terms of the BSD 3-Clause license. See the accompanying LICENSE file.
#ifndef BLAS_IAMAX_HH
#define BLAS_IAMAX_HH
#include "blas/util.hh"
#include <limits>
namespace blas {
// =============================================================================
/// @return Index of infinity-norm of vector, $|| x ||_{inf}$,
/// $\text{argmax}_{i=0}^{n-1} |Re(x_i)| + |Im(x_i)|$.
/// Returns -1 if n = 0.
///
/// Generic implementation for arbitrary data types.
///
/// @param[in] n
/// Number of elements in x. n >= 0.
///
/// @param[in] x
/// The n-element vector x, in an array of length (n-1)*incx + 1.
///
/// @param[in] incx
/// Stride between elements of x. incx > 0.
///
/// @ingroup iamax
template <typename T>
int64_t iamax(
int64_t n,
T const *x, int64_t incx )
{
typedef real_type<T> real_t;
// check arguments
blas_error_if( n < 0 ); // standard BLAS returns, doesn't fail
blas_error_if( incx <= 0 ); // standard BLAS returns, doesn't fail
// todo: check NAN
real_t result = -1;
int64_t index = -1;
if (incx == 1) {
// unit stride
for (int64_t i = 0; i < n; ++i) {
real_t tmp = abs1( x[i] );
if (tmp > result) {
result = tmp;
index = i;
}
}
}
else {
// non-unit stride
int64_t ix = 0;
for (int64_t i = 0; i < n; ++i) {
real_t tmp = abs1( x[ix] );
if (tmp > result) {
result = tmp;
index = i;
}
ix += incx;
}
}
return index;
}
} // namespace blas
#endif // #ifndef BLAS_IAMAX_HH
|
