iamax.hpp

Go to the documentation of this file.

 
//
// Copyright (c) 2017-2021, University of Tennessee. All rights reserved.
// Copyright (c) 2021-2023, University of Colorado Denver. All rights reserved.
//
// This file is part of <T>LAPACK.
// <T>LAPACK 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 TLAPACK_BLAS_IAMAX_HH
#define TLAPACK_BLAS_IAMAX_HH
 
#include "tlapack/base/utils.hpp"
 
namespace tlapack {
 
template <class abs_f>
struct IamaxOpts : public EcOpts {
    constexpr IamaxOpts(abs_f absf, const EcOpts& opts = {})
        : EcOpts(opts), absf(absf){};
 
    abs_f absf;  
};
 
template <TLAPACK_VECTOR vector_t, class abs_f>
size_type<vector_t> iamax_ec(const vector_t& x, abs_f absf)
{
    // data traits
    using idx_t = size_type<vector_t>;
    using T = type_t<vector_t>;
    using real_t = real_type<T>;
 
    // constants
    const real_t oneFourth(0.25);
    const idx_t n = size(x);
 
    // quick return
    if (n <= 0) return 0;
 
    bool scaledsmax = false;  // indicates whether |x_i| = Inf
    real_t smax(-1);
    idx_t index = -1;
    idx_t i = 0;
    for (; i < n; ++i) {
        if (isnan(x[i])) {
            // return when first NaN found
            return i;
        }
        else if (isinf(x[i])) {
            // keep looking for first NaN
            for (idx_t k = i + 1; k < n; ++k) {
                if (isnan(x[k])) {
                    // return when first NaN found
                    return k;
                }
            }
 
            // return the position of the first Inf
            return i;
        }
        else {  // still no Inf found yet
            if (is_real<T>) {
                real_t a = absf(x[i]);
                if (a > smax) {
                    smax = a;
                    index = i;
                }
            }
            else if (!scaledsmax) {  // no |x_i| = Inf  yet
                real_t a = absf(x[i]);
                if (isinf(a)) {
                    scaledsmax = true;
                    smax = absf(oneFourth * x[i]);
                    index = i;
                }
                else if (a > smax) {
                    smax = a;
                    index = i;
                }
            }
            else {  // scaledsmax = true
                real_t a = absf(oneFourth * x[i]);
                if (a > smax) {
                    smax = a;
                    index = i;
                }
            }
        }
    }
 
    return index;
}
 
template <TLAPACK_VECTOR vector_t, class abs_f>
size_type<vector_t> iamax_nc(const vector_t& x, abs_f absf)
{
    // data traits
    using idx_t = size_type<vector_t>;
    using T = type_t<vector_t>;
    using real_t = real_type<T>;
 
    // constants
    const real_t oneFourth(0.25);
    const idx_t n = size(x);
 
    // quick return
    if (n <= 0) return 0;
 
    bool scaledsmax = false;  // indicates whether |x_i| = Inf
    real_t smax(-1);
    idx_t index = -1;
    idx_t i = 0;
    for (; i < n; ++i) {
        if (isinf(x[i])) {
            // return the position of the first Inf
            return i;
        }
        else {  // still no Inf found yet
            if (is_real<T>) {
                real_t a = absf(x[i]);
                if (a > smax) {
                    smax = a;
                    index = i;
                }
            }
            else if (!scaledsmax) {  // no |x_i| = Inf  yet
                real_t a = absf(x[i]);
                if (isinf(a)) {
                    scaledsmax = true;
                    smax = absf(oneFourth * x[i]);
                    index = i;
                }
                else if (a > smax) {
                    smax = a;
                    index = i;
                }
            }
            else {  // scaledsmax = true
                real_t a = absf(oneFourth * x[i]);
                if (a > smax) {
                    smax = a;
                    index = i;
                }
            }
        }
    }
 
    return (index != idx_t(-1)) ? index : 0;
}
 
template <TLAPACK_VECTOR vector_t, class abs_f>
size_type<vector_t> iamax(const vector_t& x, const IamaxOpts<abs_f>& opts)
{
    return (opts.ec.nan == true) ? iamax_ec(x, opts.absf)
                                 : iamax_nc(x, opts.absf);
}
 
template <TLAPACK_VECTOR vector_t, disable_if_allow_optblas_t<vector_t> = 0>
size_type<vector_t> iamax(const vector_t& x)
{
    using T = type_t<vector_t>;
    using real_t = real_type<T>;
 
    return iamax(x, IamaxOpts([](const T& x) -> real_t { return abs1(x); }));
}
 
#ifdef TLAPACK_USE_LAPACKPP
 
template <TLAPACK_LEGACY_VECTOR vector_t,
          enable_if_allow_optblas_t<vector_t> = 0>
size_type<vector_t> iamax(vector_t const& x)
{
    // Legacy objects
    auto x_ = legacy_vector(x);
 
    // Constants to forward
    const auto& n = x_.n;
 
    return ::blas::iamax(n, x_.ptr, x_.inc);
}
 
#endif
 
}  // namespace tlapack
 
#endif  //  #ifndef TLAPACK_BLAS_IAMAX_HH