Mercurial > octave
changeset 31810:ce32aca67acf
maint: Remove unnecessary indent after OCTAVE_BEGIN_NAMESPACE in liboctave/
* Range.h, oct-norm.cc, lo-utils.cc:
Remove unnecessary indent after OCTAVE_BEGIN_NAMESPACE in liboctave/.
| author | Rik <rik@octave.org> |
|---|---|
| date | Thu, 02 Feb 2023 12:34:03 -0800 |
| parents | 32acdc376a36 |
| children | 8e6a9cf412c1 |
| files | liboctave/array/Range.h liboctave/numeric/oct-norm.cc liboctave/util/lo-utils.cc |
| diffstat | 3 files changed, 1198 insertions(+), 1198 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/array/Range.h Thu Feb 02 12:32:15 2023 -0800 +++ b/liboctave/array/Range.h Thu Feb 02 12:34:03 2023 -0800 @@ -40,359 +40,359 @@ OCTAVE_BEGIN_NAMESPACE(octave) - // For now, only define for floating point types. However, we only - // need range<float> as a temporary local variable in make_float_range - // in ov.cc. +// For now, only define for floating point types. However, we only +// need range<float> as a temporary local variable in make_float_range +// in ov.cc. - template <typename T> - class - range<T, typename std::enable_if<std::is_floating_point<T>::value>::type> - { - public: +template <typename T> +class +range<T, typename std::enable_if<std::is_floating_point<T>::value>::type> +{ +public: - range () - : m_base (0), m_increment (0), m_limit (0), m_final (0), m_numel (0), - m_reverse (false) - { } + range () + : m_base (0), m_increment (0), m_limit (0), m_final (0), m_numel (0), + m_reverse (false) + { } - // LIMIT is an upper limit and may be outside the range of actual - // values. For floating point ranges, we perform a tolerant check - // to attempt to capture limit in the set of values if it is "close" - // to the value of base + a multiple of the increment. + // LIMIT is an upper limit and may be outside the range of actual + // values. For floating point ranges, we perform a tolerant check + // to attempt to capture limit in the set of values if it is "close" + // to the value of base + a multiple of the increment. - range (const T& base, const T& increment, const T& limit, - bool reverse = false) - : m_base (base), m_increment (increment), m_limit (limit), - m_final (), m_numel (), m_reverse (reverse) - { - init (); - } + range (const T& base, const T& increment, const T& limit, + bool reverse = false) + : m_base (base), m_increment (increment), m_limit (limit), + m_final (), m_numel (), m_reverse (reverse) + { + init (); + } - range (const T& base, const T& limit) - : m_base (base), m_increment (1), m_limit (limit), m_final (), m_numel (), - m_reverse (false) - { - init (); - } + range (const T& base, const T& limit) + : m_base (base), m_increment (1), m_limit (limit), m_final (), m_numel (), + m_reverse (false) + { + init (); + } - // Allow conversion from (presumably) properly constructed Range - // objects. The values of base, limit, increment, and numel must be - // consistent. + // Allow conversion from (presumably) properly constructed Range + // objects. The values of base, limit, increment, and numel must be + // consistent. - // FIXME: Actually check that base, limit, increment, and numel are - // consistent? + // FIXME: Actually check that base, limit, increment, and numel are + // consistent? - range (const T& base, const T& increment, const T& limit, - octave_idx_type numel, bool reverse = false) - : m_base (base), m_increment (increment), m_limit (limit), - m_final (limit), m_numel (numel), m_reverse (reverse) - { } + range (const T& base, const T& increment, const T& limit, + octave_idx_type numel, bool reverse = false) + : m_base (base), m_increment (increment), m_limit (limit), + m_final (limit), m_numel (numel), m_reverse (reverse) + { } - // We don't use a constructor for this because it will conflict with - // range<T> (base, limit, increment) when T is octave_idx_type. + // We don't use a constructor for this because it will conflict with + // range<T> (base, limit, increment) when T is octave_idx_type. - static range<T> make_n_element_range (const T& base, const T& increment, - octave_idx_type numel, - bool reverse = false) - { - // We could just make this constructor public, but it allows - // inconsistent ranges to be constructed. And it is probably much - // clearer to see "make_n_element_range" instead of puzzling over the - // purpose of this strange constructor form. + static range<T> make_n_element_range (const T& base, const T& increment, + octave_idx_type numel, + bool reverse = false) + { + // We could just make this constructor public, but it allows + // inconsistent ranges to be constructed. And it is probably much + // clearer to see "make_n_element_range" instead of puzzling over the + // purpose of this strange constructor form. - T final_val = (reverse ? base - (numel - 1) * increment - : base + (numel - 1) * increment); + T final_val = (reverse ? base - (numel - 1) * increment + : base + (numel - 1) * increment); - return range<T> (base, increment, final_val, numel, reverse); - } + return range<T> (base, increment, final_val, numel, reverse); + } - range (const range<T>& r) - : m_base (r.m_base), m_increment (r.m_increment), - m_limit (r.m_limit), m_final (r.m_final), - m_numel (r.m_numel), m_reverse (r.m_reverse) - { } + range (const range<T>& r) + : m_base (r.m_base), m_increment (r.m_increment), + m_limit (r.m_limit), m_final (r.m_final), + m_numel (r.m_numel), m_reverse (r.m_reverse) + { } - range<T>& operator = (const range<T>& r) - { - if (this != &r) - { - m_base = r.m_base; - m_increment = r.m_increment; - m_limit = r.m_limit; - m_final = r.m_final; - m_numel = r.m_numel; - m_reverse = r.m_reverse; - } + range<T>& operator = (const range<T>& r) + { + if (this != &r) + { + m_base = r.m_base; + m_increment = r.m_increment; + m_limit = r.m_limit; + m_final = r.m_final; + m_numel = r.m_numel; + m_reverse = r.m_reverse; + } - return *this; - } + return *this; + } - ~range () = default; + ~range () = default; - T base () const { return m_base; } - T increment () const { return m_increment; } - T limit () const { return m_limit; } - bool reverse () const { return m_reverse; } + T base () const { return m_base; } + T increment () const { return m_increment; } + T limit () const { return m_limit; } + bool reverse () const { return m_reverse; } - T final_value () const { return m_final; } + T final_value () const { return m_final; } - T min () const - { - return (m_numel > 0 - ? ((m_reverse ? m_increment > T (0) - : m_increment > T (0)) ? base () : final_value ()) - : T (0)); - } + T min () const + { + return (m_numel > 0 + ? ((m_reverse ? m_increment > T (0) + : m_increment > T (0)) ? base () : final_value ()) + : T (0)); + } - T max () const - { - return (m_numel > 0 - ? ((m_reverse ? m_increment < T (0) - : m_increment > T (0)) ? final_value () : base ()) - : T (0)); - } + T max () const + { + return (m_numel > 0 + ? ((m_reverse ? m_increment < T (0) + : m_increment > T (0)) ? final_value () : base ()) + : T (0)); + } - octave_idx_type numel () const { return m_numel; } + octave_idx_type numel () const { return m_numel; } - // To support things like "for i = 1:Inf; ...; end" that are - // required for Matlab compatibility, creation of a range object - // like 1:Inf is allowed with m_numel set to - // numeric_limits<octave_idx_type>::max(). However, it is not - // possible to store these ranges. The following function allows - // us to easily distinguish ranges with an infinite number of - // elements. There are specializations for double and float. + // To support things like "for i = 1:Inf; ...; end" that are + // required for Matlab compatibility, creation of a range object + // like 1:Inf is allowed with m_numel set to + // numeric_limits<octave_idx_type>::max(). However, it is not + // possible to store these ranges. The following function allows + // us to easily distinguish ranges with an infinite number of + // elements. There are specializations for double and float. - bool is_storable () const { return true; } + bool is_storable () const { return true; } - dim_vector dims () const { return dim_vector (1, m_numel); } + dim_vector dims () const { return dim_vector (1, m_numel); } - octave_idx_type rows () const { return 1; } + octave_idx_type rows () const { return 1; } - octave_idx_type cols () const { return numel (); } - octave_idx_type columns () const { return numel (); } + octave_idx_type cols () const { return numel (); } + octave_idx_type columns () const { return numel (); } - bool isempty () const { return numel () == 0; } + bool isempty () const { return numel () == 0; } - bool all_elements_are_ints () const { return true; } + bool all_elements_are_ints () const { return true; } - sortmode issorted (sortmode mode = ASCENDING) const - { - if (m_numel > 1 && (m_reverse ? m_increment < T (0) - : m_increment > T (0))) - mode = ((mode == DESCENDING) ? UNSORTED : ASCENDING); - else if (m_numel > 1 && (m_reverse ? m_increment > T (0) - : m_increment < T (0))) - mode = ((mode == ASCENDING) ? UNSORTED : DESCENDING); - else - mode = ((mode == UNSORTED) ? ASCENDING : mode); + sortmode issorted (sortmode mode = ASCENDING) const + { + if (m_numel > 1 && (m_reverse ? m_increment < T (0) + : m_increment > T (0))) + mode = ((mode == DESCENDING) ? UNSORTED : ASCENDING); + else if (m_numel > 1 && (m_reverse ? m_increment > T (0) + : m_increment < T (0))) + mode = ((mode == ASCENDING) ? UNSORTED : DESCENDING); + else + mode = ((mode == UNSORTED) ? ASCENDING : mode); - return mode; - } + return mode; + } - OCTAVE_API octave_idx_type nnz () const; + OCTAVE_API octave_idx_type nnz () const; - // Support for single-index subscripting, without generating matrix cache. + // Support for single-index subscripting, without generating matrix cache. - T checkelem (octave_idx_type i) const - { - if (i < 0 || i >= m_numel) - err_index_out_of_range (2, 2, i+1, m_numel, dims ()); + T checkelem (octave_idx_type i) const + { + if (i < 0 || i >= m_numel) + err_index_out_of_range (2, 2, i+1, m_numel, dims ()); - if (i == 0) - // Required for proper NaN handling. - return (m_numel == 1 ? final_value () : m_base); - else if (i < m_numel - 1) - return (m_reverse ? m_base + T (i) * m_increment - : m_base + T (i) * m_increment); - else - return final_value (); - } + if (i == 0) + // Required for proper NaN handling. + return (m_numel == 1 ? final_value () : m_base); + else if (i < m_numel - 1) + return (m_reverse ? m_base + T (i) * m_increment + : m_base + T (i) * m_increment); + else + return final_value (); + } - T checkelem (octave_idx_type i, octave_idx_type j) const - { - // Ranges are *always* row vectors. - if (i != 0) - err_index_out_of_range (1, 1, i+1, m_numel, dims ()); + T checkelem (octave_idx_type i, octave_idx_type j) const + { + // Ranges are *always* row vectors. + if (i != 0) + err_index_out_of_range (1, 1, i+1, m_numel, dims ()); - return checkelem (j); - } + return checkelem (j); + } - T elem (octave_idx_type i) const - { - if (i == 0) - // Required for proper NaN handling. - return (m_numel == 1 ? final_value () : m_base); - else if (i < m_numel - 1) - return (m_reverse ? m_base - T (i) * m_increment - : m_base + T (i) * m_increment); - else - return final_value (); - } + T elem (octave_idx_type i) const + { + if (i == 0) + // Required for proper NaN handling. + return (m_numel == 1 ? final_value () : m_base); + else if (i < m_numel - 1) + return (m_reverse ? m_base - T (i) * m_increment + : m_base + T (i) * m_increment); + else + return final_value (); + } - T elem (octave_idx_type /* i */, octave_idx_type j) const - { - return elem (j); - } + T elem (octave_idx_type /* i */, octave_idx_type j) const + { + return elem (j); + } - T operator () (octave_idx_type i) const - { - return elem (i); - } + T operator () (octave_idx_type i) const + { + return elem (i); + } - T operator () (octave_idx_type i, octave_idx_type j) const - { - return elem (i, j); - } + T operator () (octave_idx_type i, octave_idx_type j) const + { + return elem (i, j); + } - Array<T> index (const idx_vector& idx) const - { - Array<T> retval; + Array<T> index (const idx_vector& idx) const + { + Array<T> retval; - octave_idx_type n = m_numel; + octave_idx_type n = m_numel; - if (idx.is_colon ()) - { - retval = array_value ().reshape (dim_vector (m_numel, 1)); - } - else - { - if (idx.extent (n) != n) - err_index_out_of_range (1, 1, idx.extent (n), n, dims ()); + if (idx.is_colon ()) + { + retval = array_value ().reshape (dim_vector (m_numel, 1)); + } + else + { + if (idx.extent (n) != n) + err_index_out_of_range (1, 1, idx.extent (n), n, dims ()); - dim_vector idx_dims = idx.orig_dimensions (); - octave_idx_type idx_len = idx.length (n); + dim_vector idx_dims = idx.orig_dimensions (); + octave_idx_type idx_len = idx.length (n); - // taken from Array.cc. - if (n != 1 && idx_dims.isvector ()) - idx_dims = dim_vector (1, idx_len); + // taken from Array.cc. + if (n != 1 && idx_dims.isvector ()) + idx_dims = dim_vector (1, idx_len); - retval.clear (idx_dims); + retval.clear (idx_dims); - // Loop over all values in IDX, executing the lambda - // expression for each index value. + // Loop over all values in IDX, executing the lambda + // expression for each index value. - T *array = retval.fortran_vec (); + T *array = retval.fortran_vec (); - idx.loop (n, [=, &array] (octave_idx_type i) - { - if (i == 0) - // Required for proper NaN handling. - *array++ = (m_numel == 0 ? m_final : m_base); - else if (i < m_numel - 1) - *array++ = (m_reverse ? m_base - T (i) * m_increment - : m_base + T (i) * m_increment); - else - *array++ = m_final; - }); - } + idx.loop (n, [=, &array] (octave_idx_type i) + { + if (i == 0) + // Required for proper NaN handling. + *array++ = (m_numel == 0 ? m_final : m_base); + else if (i < m_numel - 1) + *array++ = (m_reverse ? m_base - T (i) * m_increment + : m_base + T (i) * m_increment); + else + *array++ = m_final; + }); + } - return retval; - } + return retval; + } - Array<T> diag (octave_idx_type k) const - { - return array_value ().diag (k); - } + Array<T> diag (octave_idx_type k) const + { + return array_value ().diag (k); + } - Array<T> array_value () const - { - octave_idx_type nel = numel (); + Array<T> array_value () const + { + octave_idx_type nel = numel (); - Array<T> retval (dim_vector (1, nel)); + Array<T> retval (dim_vector (1, nel)); - if (nel == 1) - // Required for proper NaN handling. - retval(0) = final_value (); - else if (nel > 1) - { - // The first element must always be *exactly* the base. - // E.g, -0 would otherwise become +0 in the loop (-0 + 0*increment). - retval(0) = m_base; + if (nel == 1) + // Required for proper NaN handling. + retval(0) = final_value (); + else if (nel > 1) + { + // The first element must always be *exactly* the base. + // E.g, -0 would otherwise become +0 in the loop (-0 + 0*increment). + retval(0) = m_base; - if (m_reverse) - for (octave_idx_type i = 1; i < nel - 1; i++) - retval.xelem (i) = m_base - i * m_increment; - else - for (octave_idx_type i = 1; i < nel - 1; i++) - retval.xelem (i) = m_base + i * m_increment; + if (m_reverse) + for (octave_idx_type i = 1; i < nel - 1; i++) + retval.xelem (i) = m_base - i * m_increment; + else + for (octave_idx_type i = 1; i < nel - 1; i++) + retval.xelem (i) = m_base + i * m_increment; - retval.xelem (nel - 1) = final_value (); - } + retval.xelem (nel - 1) = final_value (); + } - return retval; - } + return retval; + } - private: +private: - T m_base; - T m_increment; - T m_limit; - T m_final; - octave_idx_type m_numel; - bool m_reverse; + T m_base; + T m_increment; + T m_limit; + T m_final; + octave_idx_type m_numel; + bool m_reverse; - // Setting the number of elements to zero when the increment is zero - // is intentional and matches the behavior of Matlab's colon - // operator. + // Setting the number of elements to zero when the increment is zero + // is intentional and matches the behavior of Matlab's colon + // operator. - // These calculations are appropriate for integer ranges. There are - // specializations for double and float. + // These calculations are appropriate for integer ranges. There are + // specializations for double and float. - void init () - { - if (m_reverse) - { - m_numel = ((m_increment == T (0) - || (m_limit > m_base && m_increment > T (0)) - || (m_limit < m_base && m_increment < T (0))) - ? T (0) - : (m_base - m_limit - m_increment) / m_increment); + void init () + { + if (m_reverse) + { + m_numel = ((m_increment == T (0) + || (m_limit > m_base && m_increment > T (0)) + || (m_limit < m_base && m_increment < T (0))) + ? T (0) + : (m_base - m_limit - m_increment) / m_increment); - m_final = m_base - (m_numel - 1) * m_increment; - } - else - { - m_numel = ((m_increment == T (0) - || (m_limit > m_base && m_increment < T (0)) - || (m_limit < m_base && m_increment > T (0))) - ? T (0) - : (m_limit - m_base + m_increment) / m_increment); + m_final = m_base - (m_numel - 1) * m_increment; + } + else + { + m_numel = ((m_increment == T (0) + || (m_limit > m_base && m_increment < T (0)) + || (m_limit < m_base && m_increment > T (0))) + ? T (0) + : (m_limit - m_base + m_increment) / m_increment); - m_final = m_base + (m_numel - 1) * m_increment; - } - } - }; + m_final = m_base + (m_numel - 1) * m_increment; + } + } +}; - // Specializations defined externally. +// Specializations defined externally. - template <> OCTAVE_API bool range<double>::all_elements_are_ints () const; - template <> OCTAVE_API bool range<float>::all_elements_are_ints () const; +template <> OCTAVE_API bool range<double>::all_elements_are_ints () const; +template <> OCTAVE_API bool range<float>::all_elements_are_ints () const; - template <> OCTAVE_API void range<double>::init (); - template <> OCTAVE_API void range<float>::init (); +template <> OCTAVE_API void range<double>::init (); +template <> OCTAVE_API void range<float>::init (); - // For now, only define for floating point types. However, we only - // need range<float> as a temporary local variable in make_float_range - // in ov.cc. +// For now, only define for floating point types. However, we only +// need range<float> as a temporary local variable in make_float_range +// in ov.cc. #if 0 - template <> OCTAVE_API void range<octave_int8>::init (); - template <> OCTAVE_API void range<octave_int16>::init (); - template <> OCTAVE_API void range<octave_int32>::init (); - template <> OCTAVE_API void range<octave_int64>::init (); - template <> OCTAVE_API void range<octave_uint8>::init (); - template <> OCTAVE_API void range<octave_uint16>::init (); - template <> OCTAVE_API void range<octave_uint32>::init (); - template <> OCTAVE_API void range<octave_uint64>::init (); +template <> OCTAVE_API void range<octave_int8>::init (); +template <> OCTAVE_API void range<octave_int16>::init (); +template <> OCTAVE_API void range<octave_int32>::init (); +template <> OCTAVE_API void range<octave_int64>::init (); +template <> OCTAVE_API void range<octave_uint8>::init (); +template <> OCTAVE_API void range<octave_uint16>::init (); +template <> OCTAVE_API void range<octave_uint32>::init (); +template <> OCTAVE_API void range<octave_uint64>::init (); #endif - template <> OCTAVE_API bool range<double>::is_storable () const; - template <> OCTAVE_API bool range<float>::is_storable () const; +template <> OCTAVE_API bool range<double>::is_storable () const; +template <> OCTAVE_API bool range<float>::is_storable () const; - template <> OCTAVE_API octave_idx_type range<double>::nnz () const; - template <> OCTAVE_API octave_idx_type range<float>::nnz () const; +template <> OCTAVE_API octave_idx_type range<double>::nnz () const; +template <> OCTAVE_API octave_idx_type range<float>::nnz () const; OCTAVE_END_NAMESPACE(octave)
--- a/liboctave/numeric/oct-norm.cc Thu Feb 02 12:32:15 2023 -0800 +++ b/liboctave/numeric/oct-norm.cc Thu Feb 02 12:34:03 2023 -0800 @@ -65,564 +65,564 @@ OCTAVE_BEGIN_NAMESPACE(octave) - // Theory: norm accumulator is an object that has an accum method able - // to handle both real and complex element, and a cast operator - // returning the intermediate norm. Reference: Higham, N. "Estimating - // the Matrix p-Norm." Numer. Math. 62, 539-555, 1992. - - // norm accumulator for the p-norm - template <typename R> - class norm_accumulator_p - { - public: - norm_accumulator_p () { } // we need this one for Array - norm_accumulator_p (R pp) : m_p(pp), m_scl(0), m_sum(1) { } +// Theory: norm accumulator is an object that has an accum method able +// to handle both real and complex element, and a cast operator +// returning the intermediate norm. Reference: Higham, N. "Estimating +// the Matrix p-Norm." Numer. Math. 62, 539-555, 1992. - template <typename U> - void accum (U val) - { - octave_quit (); - R t = std::abs (val); - if (m_scl == t) // we need this to handle Infs properly - m_sum += 1; - else if (m_scl < t) - { - m_sum *= std::pow (m_scl/t, m_p); - m_sum += 1; - m_scl = t; - } - else if (t != 0) - m_sum += std::pow (t/m_scl, m_p); - } - - operator R () { return m_scl * std::pow (m_sum, 1/m_p); } - - private: - R m_p, m_scl, m_sum; - }; - - // norm accumulator for the minus p-pseudonorm - template <typename R> - class norm_accumulator_mp - { - public: - norm_accumulator_mp () { } // we need this one for Array - norm_accumulator_mp (R pp) : m_p(pp), m_scl(0), m_sum(1) { } +// norm accumulator for the p-norm +template <typename R> +class norm_accumulator_p +{ +public: + norm_accumulator_p () { } // we need this one for Array + norm_accumulator_p (R pp) : m_p(pp), m_scl(0), m_sum(1) { } - template <typename U> - void accum (U val) - { - octave_quit (); - R t = 1 / std::abs (val); - if (m_scl == t) + template <typename U> + void accum (U val) + { + octave_quit (); + R t = std::abs (val); + if (m_scl == t) // we need this to handle Infs properly + m_sum += 1; + else if (m_scl < t) + { + m_sum *= std::pow (m_scl/t, m_p); m_sum += 1; - else if (m_scl < t) - { - m_sum *= std::pow (m_scl/t, m_p); - m_sum += 1; - m_scl = t; - } - else if (t != 0) - m_sum += std::pow (t/m_scl, m_p); - } - - operator R () { return m_scl * std::pow (m_sum, -1/m_p); } - - private: - R m_p, m_scl, m_sum; - }; - - // norm accumulator for the 2-norm (euclidean) - template <typename R> - class norm_accumulator_2 - { - public: - norm_accumulator_2 () : m_scl(0), m_sum(1) { } - - void accum (R val) - { - R t = std::abs (val); - if (m_scl == t) - m_sum += 1; - else if (m_scl < t) - { - m_sum *= pow2 (m_scl/t); - m_sum += 1; - m_scl = t; - } - else if (t != 0) - m_sum += pow2 (t/m_scl); - } - - void accum (std::complex<R> val) - { - accum (val.real ()); - accum (val.imag ()); - } + m_scl = t; + } + else if (t != 0) + m_sum += std::pow (t/m_scl, m_p); + } - operator R () { return m_scl * std::sqrt (m_sum); } - - private: - static inline R pow2 (R x) { return x*x; } - - //-------- - - R m_scl, m_sum; - }; + operator R () { return m_scl * std::pow (m_sum, 1/m_p); } - // norm accumulator for the 1-norm (city metric) - template <typename R> - class norm_accumulator_1 - { - public: - norm_accumulator_1 () : m_sum (0) { } - template <typename U> - void accum (U val) - { - m_sum += std::abs (val); - } - - operator R () { return m_sum; } +private: + R m_p, m_scl, m_sum; +}; - private: - R m_sum; - }; - - // norm accumulator for the inf-norm (max metric) - template <typename R> - class norm_accumulator_inf - { - public: - norm_accumulator_inf () : m_max (0) { } - template <typename U> - void accum (U val) - { - if (math::isnan (val)) - m_max = numeric_limits<R>::NaN (); - else - m_max = std::max (m_max, std::abs (val)); - } - - operator R () { return m_max; } - - private: - R m_max; - }; +// norm accumulator for the minus p-pseudonorm +template <typename R> +class norm_accumulator_mp +{ +public: + norm_accumulator_mp () { } // we need this one for Array + norm_accumulator_mp (R pp) : m_p(pp), m_scl(0), m_sum(1) { } - // norm accumulator for the -inf pseudonorm (min abs value) - template <typename R> - class norm_accumulator_minf + template <typename U> + void accum (U val) { - public: - norm_accumulator_minf () : m_min (numeric_limits<R>::Inf ()) { } - template <typename U> - void accum (U val) - { - if (math::isnan (val)) - m_min = numeric_limits<R>::NaN (); - else - m_min = std::min (m_min, std::abs (val)); - } - - operator R () { return m_min; } - - private: - R m_min; - }; - - // norm accumulator for the 0-pseudonorm (hamming distance) - template <typename R> - class norm_accumulator_0 - { - public: - norm_accumulator_0 () : m_num (0) { } - template <typename U> - void accum (U val) - { - if (val != static_cast<U> (0)) ++m_num; - } - - operator R () { return m_num; } - - private: - unsigned int m_num; - }; - - // OK, we're armed :) Now let's go for the fun - - template <typename T, typename R, typename ACC> - inline void vector_norm (const Array<T>& v, R& res, ACC acc) - { - for (octave_idx_type i = 0; i < v.numel (); i++) - acc.accum (v(i)); - - res = acc; + octave_quit (); + R t = 1 / std::abs (val); + if (m_scl == t) + m_sum += 1; + else if (m_scl < t) + { + m_sum *= std::pow (m_scl/t, m_p); + m_sum += 1; + m_scl = t; + } + else if (t != 0) + m_sum += std::pow (t/m_scl, m_p); } - // dense versions - template <typename T, typename R, typename ACC> - void column_norms (const MArray<T>& m, MArray<R>& res, ACC acc) + operator R () { return m_scl * std::pow (m_sum, -1/m_p); } + +private: + R m_p, m_scl, m_sum; +}; + +// norm accumulator for the 2-norm (euclidean) +template <typename R> +class norm_accumulator_2 +{ +public: + norm_accumulator_2 () : m_scl(0), m_sum(1) { } + + void accum (R val) { - res = MArray<R> (dim_vector (1, m.columns ())); - for (octave_idx_type j = 0; j < m.columns (); j++) + R t = std::abs (val); + if (m_scl == t) + m_sum += 1; + else if (m_scl < t) { - ACC accj = acc; - for (octave_idx_type i = 0; i < m.rows (); i++) - accj.accum (m(i, j)); - - res.xelem (j) = accj; + m_sum *= pow2 (m_scl/t); + m_sum += 1; + m_scl = t; } + else if (t != 0) + m_sum += pow2 (t/m_scl); } - template <typename T, typename R, typename ACC> - void row_norms (const MArray<T>& m, MArray<R>& res, ACC acc) + void accum (std::complex<R> val) { - res = MArray<R> (dim_vector (m.rows (), 1)); - std::vector<ACC> acci (m.rows (), acc); - for (octave_idx_type j = 0; j < m.columns (); j++) - { - for (octave_idx_type i = 0; i < m.rows (); i++) - acci[i].accum (m(i, j)); - } + accum (val.real ()); + accum (val.imag ()); + } + + operator R () { return m_scl * std::sqrt (m_sum); } + +private: + static inline R pow2 (R x) { return x*x; } + + //-------- - for (octave_idx_type i = 0; i < m.rows (); i++) - res.xelem (i) = acci[i]; + R m_scl, m_sum; +}; + +// norm accumulator for the 1-norm (city metric) +template <typename R> +class norm_accumulator_1 +{ +public: + norm_accumulator_1 () : m_sum (0) { } + template <typename U> + void accum (U val) + { + m_sum += std::abs (val); } - // sparse versions - template <typename T, typename R, typename ACC> - void column_norms (const MSparse<T>& m, MArray<R>& res, ACC acc) + operator R () { return m_sum; } + +private: + R m_sum; +}; + +// norm accumulator for the inf-norm (max metric) +template <typename R> +class norm_accumulator_inf +{ +public: + norm_accumulator_inf () : m_max (0) { } + template <typename U> + void accum (U val) { - res = MArray<R> (dim_vector (1, m.columns ())); - for (octave_idx_type j = 0; j < m.columns (); j++) - { - ACC accj = acc; - for (octave_idx_type k = m.cidx (j); k < m.cidx (j+1); k++) - accj.accum (m.data (k)); - - res.xelem (j) = accj; - } + if (math::isnan (val)) + m_max = numeric_limits<R>::NaN (); + else + m_max = std::max (m_max, std::abs (val)); } - template <typename T, typename R, typename ACC> - void row_norms (const MSparse<T>& m, MArray<R>& res, ACC acc) + operator R () { return m_max; } + +private: + R m_max; +}; + +// norm accumulator for the -inf pseudonorm (min abs value) +template <typename R> +class norm_accumulator_minf +{ +public: + norm_accumulator_minf () : m_min (numeric_limits<R>::Inf ()) { } + template <typename U> + void accum (U val) { - res = MArray<R> (dim_vector (m.rows (), 1)); - std::vector<ACC> acci (m.rows (), acc); - for (octave_idx_type j = 0; j < m.columns (); j++) - { - for (octave_idx_type k = m.cidx (j); k < m.cidx (j+1); k++) - acci[m.ridx (k)].accum (m.data (k)); - } - - for (octave_idx_type i = 0; i < m.rows (); i++) - res.xelem (i) = acci[i]; + if (math::isnan (val)) + m_min = numeric_limits<R>::NaN (); + else + m_min = std::min (m_min, std::abs (val)); } - // now the dispatchers -#define DEFINE_DISPATCHER(FCN_NAME, ARG_TYPE, RES_TYPE) \ - template <typename T, typename R> \ - RES_TYPE FCN_NAME (const ARG_TYPE& v, R p) \ - { \ - RES_TYPE res; \ - if (p == 2) \ - FCN_NAME (v, res, norm_accumulator_2<R> ()); \ - else if (p == 1) \ - FCN_NAME (v, res, norm_accumulator_1<R> ()); \ - else if (lo_ieee_isinf (p)) \ - { \ - if (p > 0) \ - FCN_NAME (v, res, norm_accumulator_inf<R> ()); \ - else \ - FCN_NAME (v, res, norm_accumulator_minf<R> ()); \ - } \ - else if (p == 0) \ - FCN_NAME (v, res, norm_accumulator_0<R> ()); \ - else if (p > 0) \ - FCN_NAME (v, res, norm_accumulator_p<R> (p)); \ - else \ - FCN_NAME (v, res, norm_accumulator_mp<R> (p)); \ - return res; \ + operator R () { return m_min; } + +private: + R m_min; +}; + +// norm accumulator for the 0-pseudonorm (hamming distance) +template <typename R> +class norm_accumulator_0 +{ +public: + norm_accumulator_0 () : m_num (0) { } + template <typename U> + void accum (U val) + { + if (val != static_cast<U> (0)) ++m_num; } - DEFINE_DISPATCHER (vector_norm, MArray<T>, R) - DEFINE_DISPATCHER (column_norms, MArray<T>, MArray<R>) - DEFINE_DISPATCHER (row_norms, MArray<T>, MArray<R>) - DEFINE_DISPATCHER (column_norms, MSparse<T>, MArray<R>) - DEFINE_DISPATCHER (row_norms, MSparse<T>, MArray<R>) + operator R () { return m_num; } + +private: + unsigned int m_num; +}; + +// OK, we're armed :) Now let's go for the fun + +template <typename T, typename R, typename ACC> +inline void vector_norm (const Array<T>& v, R& res, ACC acc) +{ + for (octave_idx_type i = 0; i < v.numel (); i++) + acc.accum (v(i)); + + res = acc; +} + +// dense versions +template <typename T, typename R, typename ACC> +void column_norms (const MArray<T>& m, MArray<R>& res, ACC acc) +{ + res = MArray<R> (dim_vector (1, m.columns ())); + for (octave_idx_type j = 0; j < m.columns (); j++) + { + ACC accj = acc; + for (octave_idx_type i = 0; i < m.rows (); i++) + accj.accum (m(i, j)); - // The approximate subproblem in Higham's method. Find lambda and mu such - // that norm ([lambda, mu], p) == 1 and norm (y*lambda + col*mu, p) is - // maximized. - // Real version. As in Higham's paper. - template <typename ColVectorT, typename R> - static void - higham_subp (const ColVectorT& y, const ColVectorT& col, - octave_idx_type nsamp, R p, R& lambda, R& mu) - { - R nrm = 0; - for (octave_idx_type i = 0; i < nsamp; i++) - { - octave_quit (); - R fi = i * static_cast<R> (M_PI) / nsamp; - R lambda1 = cos (fi); - R mu1 = sin (fi); - R lmnr = std::pow (std::pow (std::abs (lambda1), p) + - std::pow (std::abs (mu1), p), 1/p); - lambda1 /= lmnr; mu1 /= lmnr; - R nrm1 = vector_norm (lambda1 * y + mu1 * col, p); - if (nrm1 > nrm) - { - lambda = lambda1; - mu = mu1; - nrm = nrm1; - } - } - } + res.xelem (j) = accj; + } +} + +template <typename T, typename R, typename ACC> +void row_norms (const MArray<T>& m, MArray<R>& res, ACC acc) +{ + res = MArray<R> (dim_vector (m.rows (), 1)); + std::vector<ACC> acci (m.rows (), acc); + for (octave_idx_type j = 0; j < m.columns (); j++) + { + for (octave_idx_type i = 0; i < m.rows (); i++) + acci[i].accum (m(i, j)); + } + + for (octave_idx_type i = 0; i < m.rows (); i++) + res.xelem (i) = acci[i]; +} + +// sparse versions +template <typename T, typename R, typename ACC> +void column_norms (const MSparse<T>& m, MArray<R>& res, ACC acc) +{ + res = MArray<R> (dim_vector (1, m.columns ())); + for (octave_idx_type j = 0; j < m.columns (); j++) + { + ACC accj = acc; + for (octave_idx_type k = m.cidx (j); k < m.cidx (j+1); k++) + accj.accum (m.data (k)); + + res.xelem (j) = accj; + } +} - // Complex version. Higham's paper does not deal with complex case, so we - // use a simple extension. First, guess the magnitudes as in real version, - // then try to rotate lambda to improve further. - template <typename ColVectorT, typename R> - static void - higham_subp (const ColVectorT& y, const ColVectorT& col, - octave_idx_type nsamp, R p, - std::complex<R>& lambda, std::complex<R>& mu) - { - typedef std::complex<R> CR; - R nrm = 0; - lambda = 1.0; - CR lamcu = lambda / std::abs (lambda); - // Probe magnitudes - for (octave_idx_type i = 0; i < nsamp; i++) - { - octave_quit (); - R fi = i * static_cast<R> (M_PI) / nsamp; - R lambda1 = cos (fi); - R mu1 = sin (fi); - R lmnr = std::pow (std::pow (std::abs (lambda1), p) + - std::pow (std::abs (mu1), p), 1/p); - lambda1 /= lmnr; mu1 /= lmnr; - R nrm1 = vector_norm (lambda1 * lamcu * y + mu1 * col, p); - if (nrm1 > nrm) - { - lambda = lambda1 * lamcu; - mu = mu1; - nrm = nrm1; - } - } - R lama = std::abs (lambda); - // Probe orientation - for (octave_idx_type i = 0; i < nsamp; i++) - { - octave_quit (); - R fi = i * static_cast<R> (M_PI) / nsamp; - lamcu = CR (cos (fi), sin (fi)); - R nrm1 = vector_norm (lama * lamcu * y + mu * col, p); - if (nrm1 > nrm) - { - lambda = lama * lamcu; - nrm = nrm1; - } - } - } +template <typename T, typename R, typename ACC> +void row_norms (const MSparse<T>& m, MArray<R>& res, ACC acc) +{ + res = MArray<R> (dim_vector (m.rows (), 1)); + std::vector<ACC> acci (m.rows (), acc); + for (octave_idx_type j = 0; j < m.columns (); j++) + { + for (octave_idx_type k = m.cidx (j); k < m.cidx (j+1); k++) + acci[m.ridx (k)].accum (m.data (k)); + } + + for (octave_idx_type i = 0; i < m.rows (); i++) + res.xelem (i) = acci[i]; +} - // the p-dual element (should work for both real and complex) - template <typename T, typename R> - inline T elem_dual_p (T x, R p) - { - return math::signum (x) * std::pow (std::abs (x), p-1); - } +// now the dispatchers +#define DEFINE_DISPATCHER(FCN_NAME, ARG_TYPE, RES_TYPE) \ +template <typename T, typename R> \ +RES_TYPE FCN_NAME (const ARG_TYPE& v, R p) \ +{ \ + RES_TYPE res; \ + if (p == 2) \ + FCN_NAME (v, res, norm_accumulator_2<R> ()); \ + else if (p == 1) \ + FCN_NAME (v, res, norm_accumulator_1<R> ()); \ + else if (lo_ieee_isinf (p)) \ + { \ + if (p > 0) \ + FCN_NAME (v, res, norm_accumulator_inf<R> ()); \ + else \ + FCN_NAME (v, res, norm_accumulator_minf<R> ()); \ + } \ + else if (p == 0) \ + FCN_NAME (v, res, norm_accumulator_0<R> ()); \ + else if (p > 0) \ + FCN_NAME (v, res, norm_accumulator_p<R> (p)); \ + else \ + FCN_NAME (v, res, norm_accumulator_mp<R> (p)); \ + return res; \ +} - // the VectorT is used for vectors, but actually it has to be - // a Matrix type to allow all the operations. For instance SparseMatrix - // does not support multiplication with column/row vectors. - // the dual vector - template <typename VectorT, typename R> - VectorT dual_p (const VectorT& x, R p, R q) - { - VectorT res (x.dims ()); - for (octave_idx_type i = 0; i < x.numel (); i++) - res.xelem (i) = elem_dual_p (x(i), p); - return res / vector_norm (res, q); - } +DEFINE_DISPATCHER (vector_norm, MArray<T>, R) +DEFINE_DISPATCHER (column_norms, MArray<T>, MArray<R>) +DEFINE_DISPATCHER (row_norms, MArray<T>, MArray<R>) +DEFINE_DISPATCHER (column_norms, MSparse<T>, MArray<R>) +DEFINE_DISPATCHER (row_norms, MSparse<T>, MArray<R>) + +// The approximate subproblem in Higham's method. Find lambda and mu such +// that norm ([lambda, mu], p) == 1 and norm (y*lambda + col*mu, p) is +// maximized. +// Real version. As in Higham's paper. +template <typename ColVectorT, typename R> +static void +higham_subp (const ColVectorT& y, const ColVectorT& col, + octave_idx_type nsamp, R p, R& lambda, R& mu) +{ + R nrm = 0; + for (octave_idx_type i = 0; i < nsamp; i++) + { + octave_quit (); + R fi = i * static_cast<R> (M_PI) / nsamp; + R lambda1 = cos (fi); + R mu1 = sin (fi); + R lmnr = std::pow (std::pow (std::abs (lambda1), p) + + std::pow (std::abs (mu1), p), 1/p); + lambda1 /= lmnr; mu1 /= lmnr; + R nrm1 = vector_norm (lambda1 * y + mu1 * col, p); + if (nrm1 > nrm) + { + lambda = lambda1; + mu = mu1; + nrm = nrm1; + } + } +} - // Higham's hybrid method - template <typename MatrixT, typename VectorT, typename R> - R higham (const MatrixT& m, R p, R tol, int maxiter, - VectorT& x) - { - x.resize (m.columns (), 1); - // the OSE part - VectorT y(m.rows (), 1, 0), z(m.rows (), 1); - typedef typename VectorT::element_type RR; - RR lambda = 0; - RR mu = 1; - for (octave_idx_type k = 0; k < m.columns (); k++) - { - octave_quit (); - VectorT col (m.column (k)); - if (k > 0) - higham_subp (y, col, 4*k, p, lambda, mu); - for (octave_idx_type i = 0; i < k; i++) - x(i) *= lambda; - x(k) = mu; - y = lambda * y + mu * col; - } +// Complex version. Higham's paper does not deal with complex case, so we +// use a simple extension. First, guess the magnitudes as in real version, +// then try to rotate lambda to improve further. +template <typename ColVectorT, typename R> +static void +higham_subp (const ColVectorT& y, const ColVectorT& col, + octave_idx_type nsamp, R p, + std::complex<R>& lambda, std::complex<R>& mu) +{ + typedef std::complex<R> CR; + R nrm = 0; + lambda = 1.0; + CR lamcu = lambda / std::abs (lambda); + // Probe magnitudes + for (octave_idx_type i = 0; i < nsamp; i++) + { + octave_quit (); + R fi = i * static_cast<R> (M_PI) / nsamp; + R lambda1 = cos (fi); + R mu1 = sin (fi); + R lmnr = std::pow (std::pow (std::abs (lambda1), p) + + std::pow (std::abs (mu1), p), 1/p); + lambda1 /= lmnr; mu1 /= lmnr; + R nrm1 = vector_norm (lambda1 * lamcu * y + mu1 * col, p); + if (nrm1 > nrm) + { + lambda = lambda1 * lamcu; + mu = mu1; + nrm = nrm1; + } + } + R lama = std::abs (lambda); + // Probe orientation + for (octave_idx_type i = 0; i < nsamp; i++) + { + octave_quit (); + R fi = i * static_cast<R> (M_PI) / nsamp; + lamcu = CR (cos (fi), sin (fi)); + R nrm1 = vector_norm (lama * lamcu * y + mu * col, p); + if (nrm1 > nrm) + { + lambda = lama * lamcu; + nrm = nrm1; + } + } +} - // the PM part - x = x / vector_norm (x, p); - R q = p/(p-1); +// the p-dual element (should work for both real and complex) +template <typename T, typename R> +inline T elem_dual_p (T x, R p) +{ + return math::signum (x) * std::pow (std::abs (x), p-1); +} - R gamma = 0, gamma1; - int iter = 0; - while (iter < maxiter) - { - octave_quit (); - y = m*x; - gamma1 = gamma; - gamma = vector_norm (y, p); - z = dual_p (y, p, q); - z = z.hermitian (); - z = z * m; - - if (iter > 0 && (vector_norm (z, q) <= gamma - || (gamma - gamma1) <= tol*gamma)) - break; - - z = z.hermitian (); - x = dual_p (z, q, p); - iter++; - } - - return gamma; - } - - // derive column vector and SVD types +// the VectorT is used for vectors, but actually it has to be +// a Matrix type to allow all the operations. For instance SparseMatrix +// does not support multiplication with column/row vectors. +// the dual vector +template <typename VectorT, typename R> +VectorT dual_p (const VectorT& x, R p, R q) +{ + VectorT res (x.dims ()); + for (octave_idx_type i = 0; i < x.numel (); i++) + res.xelem (i) = elem_dual_p (x(i), p); + return res / vector_norm (res, q); +} - static const char *p_less1_gripe = "xnorm: p must be >= 1"; +// Higham's hybrid method +template <typename MatrixT, typename VectorT, typename R> +R higham (const MatrixT& m, R p, R tol, int maxiter, + VectorT& x) +{ + x.resize (m.columns (), 1); + // the OSE part + VectorT y(m.rows (), 1, 0), z(m.rows (), 1); + typedef typename VectorT::element_type RR; + RR lambda = 0; + RR mu = 1; + for (octave_idx_type k = 0; k < m.columns (); k++) + { + octave_quit (); + VectorT col (m.column (k)); + if (k > 0) + higham_subp (y, col, 4*k, p, lambda, mu); + for (octave_idx_type i = 0; i < k; i++) + x(i) *= lambda; + x(k) = mu; + y = lambda * y + mu * col; + } - // Static constant to control the maximum number of iterations. 100 seems to - // be a good value. Eventually, we can provide a means to change this - // constant from Octave. - static int max_norm_iter = 100; + // the PM part + x = x / vector_norm (x, p); + R q = p/(p-1); + + R gamma = 0, gamma1; + int iter = 0; + while (iter < maxiter) + { + octave_quit (); + y = m*x; + gamma1 = gamma; + gamma = vector_norm (y, p); + z = dual_p (y, p, q); + z = z.hermitian (); + z = z * m; + + if (iter > 0 && (vector_norm (z, q) <= gamma + || (gamma - gamma1) <= tol*gamma)) + break; + + z = z.hermitian (); + x = dual_p (z, q, p); + iter++; + } + + return gamma; +} + +// derive column vector and SVD types - // version with SVD for dense matrices - template <typename MatrixT, typename VectorT, typename R> - R svd_matrix_norm (const MatrixT& m, R p, VectorT) - { - R res = 0; - if (p == 2) - { - math::svd<MatrixT> fact (m, math::svd<MatrixT>::Type::sigma_only); - res = fact.singular_values () (0, 0); - } - else if (p == 1) - res = xcolnorms (m, static_cast<R> (1)).max (); - else if (lo_ieee_isinf (p) && p > 1) - res = xrownorms (m, static_cast<R> (1)).max (); - else if (p > 1) - { - VectorT x; - const R sqrteps = std::sqrt (std::numeric_limits<R>::epsilon ()); - res = higham (m, p, sqrteps, max_norm_iter, x); - } - else - (*current_liboctave_error_handler) ("%s", p_less1_gripe); +static const char *p_less1_gripe = "xnorm: p must be >= 1"; + +// Static constant to control the maximum number of iterations. 100 seems to +// be a good value. Eventually, we can provide a means to change this +// constant from Octave. +static int max_norm_iter = 100; - return res; - } +// version with SVD for dense matrices +template <typename MatrixT, typename VectorT, typename R> +R svd_matrix_norm (const MatrixT& m, R p, VectorT) +{ + R res = 0; + if (p == 2) + { + math::svd<MatrixT> fact (m, math::svd<MatrixT>::Type::sigma_only); + res = fact.singular_values () (0, 0); + } + else if (p == 1) + res = xcolnorms (m, static_cast<R> (1)).max (); + else if (lo_ieee_isinf (p) && p > 1) + res = xrownorms (m, static_cast<R> (1)).max (); + else if (p > 1) + { + VectorT x; + const R sqrteps = std::sqrt (std::numeric_limits<R>::epsilon ()); + res = higham (m, p, sqrteps, max_norm_iter, x); + } + else + (*current_liboctave_error_handler) ("%s", p_less1_gripe); - // SVD-free version for sparse matrices - template <typename MatrixT, typename VectorT, typename R> - R matrix_norm (const MatrixT& m, R p, VectorT) - { - R res = 0; - if (p == 1) - res = xcolnorms (m, static_cast<R> (1)).max (); - else if (lo_ieee_isinf (p) && p > 1) - res = xrownorms (m, static_cast<R> (1)).max (); - else if (p > 1) - { - VectorT x; - const R sqrteps = std::sqrt (std::numeric_limits<R>::epsilon ()); - res = higham (m, p, sqrteps, max_norm_iter, x); - } - else - (*current_liboctave_error_handler) ("%s", p_less1_gripe); + return res; +} - return res; - } +// SVD-free version for sparse matrices +template <typename MatrixT, typename VectorT, typename R> +R matrix_norm (const MatrixT& m, R p, VectorT) +{ + R res = 0; + if (p == 1) + res = xcolnorms (m, static_cast<R> (1)).max (); + else if (lo_ieee_isinf (p) && p > 1) + res = xrownorms (m, static_cast<R> (1)).max (); + else if (p > 1) + { + VectorT x; + const R sqrteps = std::sqrt (std::numeric_limits<R>::epsilon ()); + res = higham (m, p, sqrteps, max_norm_iter, x); + } + else + (*current_liboctave_error_handler) ("%s", p_less1_gripe); - // and finally, here's what we've promised in the header file + return res; +} + +// and finally, here's what we've promised in the header file #define DEFINE_XNORM_FCNS(PREFIX, RTYPE) \ - RTYPE xnorm (const PREFIX##ColumnVector& x, RTYPE p) \ - { \ - return vector_norm (x, p); \ - } \ - RTYPE xnorm (const PREFIX##RowVector& x, RTYPE p) \ - { \ - return vector_norm (x, p); \ - } \ - RTYPE xnorm (const PREFIX##Matrix& x, RTYPE p) \ - { \ - return svd_matrix_norm (x, p, PREFIX##Matrix ()); \ - } \ - RTYPE xfrobnorm (const PREFIX##Matrix& x) \ - { \ - return vector_norm (x, static_cast<RTYPE> (2)); \ - } +RTYPE xnorm (const PREFIX##ColumnVector& x, RTYPE p) \ +{ \ + return vector_norm (x, p); \ +} \ +RTYPE xnorm (const PREFIX##RowVector& x, RTYPE p) \ +{ \ + return vector_norm (x, p); \ +} \ +RTYPE xnorm (const PREFIX##Matrix& x, RTYPE p) \ +{ \ + return svd_matrix_norm (x, p, PREFIX##Matrix ()); \ +} \ +RTYPE xfrobnorm (const PREFIX##Matrix& x) \ +{ \ + return vector_norm (x, static_cast<RTYPE> (2)); \ +} - DEFINE_XNORM_FCNS(, double) - DEFINE_XNORM_FCNS(Complex, double) - DEFINE_XNORM_FCNS(Float, float) - DEFINE_XNORM_FCNS(FloatComplex, float) +DEFINE_XNORM_FCNS(, double) +DEFINE_XNORM_FCNS(Complex, double) +DEFINE_XNORM_FCNS(Float, float) +DEFINE_XNORM_FCNS(FloatComplex, float) - // this is needed to avoid copying the sparse matrix for xfrobnorm - template <typename T, typename R> - inline void array_norm_2 (const T *v, octave_idx_type n, R& res) - { - norm_accumulator_2<R> acc; - for (octave_idx_type i = 0; i < n; i++) - acc.accum (v[i]); +// this is needed to avoid copying the sparse matrix for xfrobnorm +template <typename T, typename R> +inline void array_norm_2 (const T *v, octave_idx_type n, R& res) +{ + norm_accumulator_2<R> acc; + for (octave_idx_type i = 0; i < n; i++) + acc.accum (v[i]); - res = acc; - } + res = acc; +} #define DEFINE_XNORM_SPARSE_FCNS(PREFIX, RTYPE) \ - RTYPE xnorm (const Sparse##PREFIX##Matrix& x, RTYPE p) \ - { \ - return matrix_norm (x, p, PREFIX##Matrix ()); \ - } \ - RTYPE xfrobnorm (const Sparse##PREFIX##Matrix& x) \ - { \ - RTYPE res; \ - array_norm_2 (x.data (), x.nnz (), res); \ - return res; \ - } +RTYPE xnorm (const Sparse##PREFIX##Matrix& x, RTYPE p) \ +{ \ + return matrix_norm (x, p, PREFIX##Matrix ()); \ +} \ +RTYPE xfrobnorm (const Sparse##PREFIX##Matrix& x) \ +{ \ + RTYPE res; \ + array_norm_2 (x.data (), x.nnz (), res); \ + return res; \ +} - DEFINE_XNORM_SPARSE_FCNS(, double) - DEFINE_XNORM_SPARSE_FCNS(Complex, double) +DEFINE_XNORM_SPARSE_FCNS(, double) +DEFINE_XNORM_SPARSE_FCNS(Complex, double) #define DEFINE_COLROW_NORM_FCNS(PREFIX, RPREFIX, RTYPE) \ - RPREFIX##RowVector \ - xcolnorms (const PREFIX##Matrix& m, RTYPE p) \ - { \ - return column_norms (m, p); \ - } \ - RPREFIX##ColumnVector \ - xrownorms (const PREFIX##Matrix& m, RTYPE p) \ - { \ - return row_norms (m, p); \ - } \ +RPREFIX##RowVector \ +xcolnorms (const PREFIX##Matrix& m, RTYPE p) \ +{ \ + return column_norms (m, p); \ +} \ +RPREFIX##ColumnVector \ +xrownorms (const PREFIX##Matrix& m, RTYPE p) \ +{ \ + return row_norms (m, p); \ +} \ - DEFINE_COLROW_NORM_FCNS(, , double) - DEFINE_COLROW_NORM_FCNS(Complex, , double) - DEFINE_COLROW_NORM_FCNS(Float, Float, float) - DEFINE_COLROW_NORM_FCNS(FloatComplex, Float, float) +DEFINE_COLROW_NORM_FCNS(, , double) +DEFINE_COLROW_NORM_FCNS(Complex, , double) +DEFINE_COLROW_NORM_FCNS(Float, Float, float) +DEFINE_COLROW_NORM_FCNS(FloatComplex, Float, float) - DEFINE_COLROW_NORM_FCNS(Sparse, , double) - DEFINE_COLROW_NORM_FCNS(SparseComplex, , double) +DEFINE_COLROW_NORM_FCNS(Sparse, , double) +DEFINE_COLROW_NORM_FCNS(SparseComplex, , double) OCTAVE_END_NAMESPACE(octave)
--- a/liboctave/util/lo-utils.cc Thu Feb 02 12:32:15 2023 -0800 +++ b/liboctave/util/lo-utils.cc Thu Feb 02 12:34:03 2023 -0800 @@ -47,479 +47,479 @@ OCTAVE_BEGIN_NAMESPACE(octave) - bool is_int_or_inf_or_nan (double x) - { - return math::isnan (x) || math::x_nint (x) == x; - } +bool is_int_or_inf_or_nan (double x) +{ + return math::isnan (x) || math::x_nint (x) == x; +} + +bool too_large_for_float (double x) +{ + return (math::isfinite (x) + && fabs (x) > std::numeric_limits<float>::max ()); +} + +bool too_large_for_float (const Complex& x) +{ + return (too_large_for_float (x.real ()) + || too_large_for_float (x.imag ())); +} - bool too_large_for_float (double x) - { - return (math::isfinite (x) - && fabs (x) > std::numeric_limits<float>::max ()); - } +bool is_int_or_inf_or_nan (float x) +{ + return math::isnan (x) || math::x_nint (x) == x; +} + +// Save a string. - bool too_large_for_float (const Complex& x) - { - return (too_large_for_float (x.real ()) - || too_large_for_float (x.imag ())); - } +char * strsave (const char *s) +{ + if (! s) + return nullptr; + + int len = strlen (s); + char *tmp = new char [len+1]; + tmp = strcpy (tmp, s); + return tmp; +} - bool is_int_or_inf_or_nan (float x) - { - return math::isnan (x) || math::x_nint (x) == x; - } +std::string fgets (FILE *f) +{ + bool eof; + return fgets (f, eof); +} - // Save a string. +std::string fgets (FILE *f, bool& eof) +{ + eof = false; + + std::string retval; - char * strsave (const char *s) - { - if (! s) - return nullptr; + int grow_size = 1024; + int max_size = grow_size; + + char *buf = static_cast<char *> (std::malloc (max_size)); + if (! buf) + (*current_liboctave_error_handler) ("octave_fgets: unable to malloc %d bytes", max_size); - int len = strlen (s); - char *tmp = new char [len+1]; - tmp = strcpy (tmp, s); - return tmp; - } + char *bufptr = buf; + int len = 0; + + do + { + if (std::fgets (bufptr, grow_size, f)) + { + len = strlen (bufptr); - std::string fgets (FILE *f) - { - bool eof; - return fgets (f, eof); - } + if (len == grow_size - 1) + { + int tmp = bufptr - buf + grow_size - 1; + grow_size *= 2; + max_size += grow_size; + auto tmpbuf = static_cast<char *> (std::realloc (buf, max_size)); + if (! tmpbuf) + { + free (buf); + (*current_liboctave_error_handler) ("octave_fgets: unable to realloc %d bytes", max_size); + } + buf = tmpbuf; + bufptr = buf + tmp; - std::string fgets (FILE *f, bool& eof) - { - eof = false; + if (*(bufptr-1) == '\n') + { + *bufptr = '\0'; + retval = buf; + } + } + else if (bufptr[len-1] != '\n') + { + bufptr[len++] = '\n'; + bufptr[len] = '\0'; + retval = buf; + } + else + retval = buf; + } + else + { + if (len == 0) + { + eof = true; - std::string retval; + free (buf); + + buf = nullptr; + } + + break; + } + } + while (retval.empty ()); - int grow_size = 1024; - int max_size = grow_size; + free (buf); + + octave_quit (); + + return retval; +} + +std::string fgetl (FILE *f) +{ + bool eof; + return fgetl (f, eof); +} - char *buf = static_cast<char *> (std::malloc (max_size)); - if (! buf) - (*current_liboctave_error_handler) ("octave_fgets: unable to malloc %d bytes", max_size); +std::string fgetl (FILE *f, bool& eof) +{ + std::string retval = fgets (f, eof); + + if (! retval.empty () && retval.back () == '\n') + retval.pop_back (); - char *bufptr = buf; - int len = 0; + return retval; +} + +template <typename T> +T +read_value (std::istream& is) +{ + T retval; + is >> retval; + return retval; +} - do +template OCTAVE_API bool read_value<bool> (std::istream& is); +template OCTAVE_API octave_int8 read_value<octave_int8> (std::istream& is); +template OCTAVE_API octave_int16 read_value<octave_int16> (std::istream& is); +template OCTAVE_API octave_int32 read_value<octave_int32> (std::istream& is); +template OCTAVE_API octave_int64 read_value<octave_int64> (std::istream& is); +template OCTAVE_API octave_uint8 read_value<octave_uint8> (std::istream& is); +template OCTAVE_API octave_uint16 read_value<octave_uint16> (std::istream& is); +template OCTAVE_API octave_uint32 read_value<octave_uint32> (std::istream& is); +template OCTAVE_API octave_uint64 read_value<octave_uint64> (std::istream& is); + +// Note that the caller is responsible for repositioning the stream on +// failure. + +template <typename T> +T +read_inf_nan_na (std::istream& is, char c0) +{ + T val = 0.0; + + switch (c0) + { + case 'i': case 'I': { - if (std::fgets (bufptr, grow_size, f)) + char c1 = is.get (); + if (c1 == 'n' || c1 == 'N') { - len = strlen (bufptr); - - if (len == grow_size - 1) + char c2 = is.get (); + if (c2 == 'f' || c2 == 'F') { - int tmp = bufptr - buf + grow_size - 1; - grow_size *= 2; - max_size += grow_size; - auto tmpbuf = static_cast<char *> (std::realloc (buf, max_size)); - if (! tmpbuf) - { - free (buf); - (*current_liboctave_error_handler) ("octave_fgets: unable to realloc %d bytes", max_size); - } - buf = tmpbuf; - bufptr = buf + tmp; - - if (*(bufptr-1) == '\n') - { - *bufptr = '\0'; - retval = buf; - } - } - else if (bufptr[len-1] != '\n') - { - bufptr[len++] = '\n'; - bufptr[len] = '\0'; - retval = buf; + val = std::numeric_limits<T>::infinity (); + is.peek (); // Potentially set EOF bit } else - retval = buf; + is.setstate (std::ios::failbit); + } + else + is.setstate (std::ios::failbit); + } + break; + + case 'n': case 'N': + { + char c1 = is.get (); + if (c1 == 'a' || c1 == 'A') + { + char c2 = is.get (); + if (c2 == 'n' || c2 == 'N') + { + val = std::numeric_limits<T>::quiet_NaN (); + is.peek (); // Potentially set EOF bit + } + else + { + val = numeric_limits<T>::NA (); + if (c2 != std::istream::traits_type::eof ()) + is.putback (c2); + else + is.clear (is.rdstate () & ~std::ios::failbit); + } } else + is.setstate (std::ios::failbit); + } + break; + + default: + (*current_liboctave_error_handler) + ("read_inf_nan_na: invalid character '%c'", c0); + } + + return val; +} + +// Read a double value. Discard any sign on NaN and NA. + +template <typename T> +double +read_fp_value (std::istream& is) +{ + T val = 0.0; + + // FIXME: resetting stream position is likely to fail unless we are + // reading from a file. + std::streampos pos = is.tellg (); + + is >> std::ws; // skip through whitespace and advance stream pointer + + bool neg = false; + char c1 = is.get (); + switch (c1) + { + case '-': + neg = true; + OCTAVE_FALLTHROUGH; + + case '+': + { + char c2 = 0; + c2 = is.get (); + if (c2 == 'i' || c2 == 'I' || c2 == 'n' || c2 == 'N') + val = read_inf_nan_na<T> (is, c2); + else if (isspace (c2)) + is.setstate (std::ios::failbit); + else { - if (len == 0) - { - eof = true; - - free (buf); - - buf = nullptr; - } + is.putback (c2); + is >> val; + } - break; - } + if (neg && ! math::isnan (val) && ! is.fail ()) + val = -val; } - while (retval.empty ()); - - free (buf); - - octave_quit (); - - return retval; - } + break; - std::string fgetl (FILE *f) - { - bool eof; - return fgetl (f, eof); - } - - std::string fgetl (FILE *f, bool& eof) - { - std::string retval = fgets (f, eof); + case 'i': case 'I': + case 'n': case 'N': + val = read_inf_nan_na<T> (is, c1); + break; - if (! retval.empty () && retval.back () == '\n') - retval.pop_back (); - - return retval; - } - - template <typename T> - T - read_value (std::istream& is) - { - T retval; - is >> retval; - return retval; - } + default: + is.putback (c1); + is >> val; + break; + } - template OCTAVE_API bool read_value<bool> (std::istream& is); - template OCTAVE_API octave_int8 read_value<octave_int8> (std::istream& is); - template OCTAVE_API octave_int16 read_value<octave_int16> (std::istream& is); - template OCTAVE_API octave_int32 read_value<octave_int32> (std::istream& is); - template OCTAVE_API octave_int64 read_value<octave_int64> (std::istream& is); - template OCTAVE_API octave_uint8 read_value<octave_uint8> (std::istream& is); - template OCTAVE_API octave_uint16 read_value<octave_uint16> (std::istream& is); - template OCTAVE_API octave_uint32 read_value<octave_uint32> (std::istream& is); - template OCTAVE_API octave_uint64 read_value<octave_uint64> (std::istream& is); - - // Note that the caller is responsible for repositioning the stream on - // failure. - - template <typename T> - T - read_inf_nan_na (std::istream& is, char c0) - { - T val = 0.0; - - switch (c0) - { - case 'i': case 'I': + std::ios::iostate status = is.rdstate (); + if (status & std::ios::failbit) + { + // Convert MAX_VAL returned by C++ streams for very large numbers to Inf + if (val == std::numeric_limits<T>::max ()) + { + if (neg) + val = -std::numeric_limits<T>::infinity (); + else + val = std::numeric_limits<T>::infinity (); + is.clear (status & ~std::ios::failbit); + } + else { - char c1 = is.get (); - if (c1 == 'n' || c1 == 'N') - { - char c2 = is.get (); - if (c2 == 'f' || c2 == 'F') - { - val = std::numeric_limits<T>::infinity (); - is.peek (); // Potentially set EOF bit - } - else - is.setstate (std::ios::failbit); - } - else - is.setstate (std::ios::failbit); + // True error. + // Reset stream to original position, clear eof bit, pass status on. + is.clear (); + is.seekg (pos); + is.setstate (status & ~std::ios_base::eofbit); } - break; + } + + return val; +} + +template <typename T> +std::complex<T> +read_cx_fp_value (std::istream& is) +{ + T re = 0.0; + T im = 0.0; - case 'n': case 'N': + std::complex<T> cx = 0.0; + + char ch = ' '; + + while (isspace (ch)) + ch = is.get (); + + if (ch == '(') + { + re = read_value<T> (is); + ch = is.get (); + + if (ch == ',') { - char c1 = is.get (); - if (c1 == 'a' || c1 == 'A') - { - char c2 = is.get (); - if (c2 == 'n' || c2 == 'N') - { - val = std::numeric_limits<T>::quiet_NaN (); - is.peek (); // Potentially set EOF bit - } - else - { - val = numeric_limits<T>::NA (); - if (c2 != std::istream::traits_type::eof ()) - is.putback (c2); - else - is.clear (is.rdstate () & ~std::ios::failbit); - } - } + im = read_value<T> (is); + ch = is.get (); + + if (ch == ')') + cx = std::complex<T> (re, im); else is.setstate (std::ios::failbit); } - break; - - default: - (*current_liboctave_error_handler) - ("read_inf_nan_na: invalid character '%c'", c0); - } - - return val; - } - - // Read a double value. Discard any sign on NaN and NA. - - template <typename T> - double - read_fp_value (std::istream& is) - { - T val = 0.0; - - // FIXME: resetting stream position is likely to fail unless we are - // reading from a file. - std::streampos pos = is.tellg (); - - is >> std::ws; // skip through whitespace and advance stream pointer + else if (ch == ')') + cx = re; + else + is.setstate (std::ios::failbit); + } + else + { + is.putback (ch); + cx = read_value<T> (is); + } - bool neg = false; - char c1 = is.get (); - switch (c1) - { - case '-': - neg = true; - OCTAVE_FALLTHROUGH; + return cx; +} - case '+': - { - char c2 = 0; - c2 = is.get (); - if (c2 == 'i' || c2 == 'I' || c2 == 'n' || c2 == 'N') - val = read_inf_nan_na<T> (is, c2); - else if (isspace (c2)) - is.setstate (std::ios::failbit); - else - { - is.putback (c2); - is >> val; - } +// FIXME: Could we use traits and enable_if to avoid duplication in the +// following specializations? - if (neg && ! math::isnan (val) && ! is.fail ()) - val = -val; - } - break; +template <> OCTAVE_API double read_value (std::istream& is) +{ + return read_fp_value<double> (is); +} - case 'i': case 'I': - case 'n': case 'N': - val = read_inf_nan_na<T> (is, c1); - break; +template <> OCTAVE_API Complex read_value (std::istream& is) +{ + return read_cx_fp_value<double> (is); +} - default: - is.putback (c1); - is >> val; - break; - } +template <> OCTAVE_API float read_value (std::istream& is) +{ + return read_fp_value<float> (is); +} - std::ios::iostate status = is.rdstate (); - if (status & std::ios::failbit) - { - // Convert MAX_VAL returned by C++ streams for very large numbers to Inf - if (val == std::numeric_limits<T>::max ()) - { - if (neg) - val = -std::numeric_limits<T>::infinity (); - else - val = std::numeric_limits<T>::infinity (); - is.clear (status & ~std::ios::failbit); - } - else - { - // True error. - // Reset stream to original position, clear eof bit, pass status on. - is.clear (); - is.seekg (pos); - is.setstate (status & ~std::ios_base::eofbit); - } - } +template <> OCTAVE_API FloatComplex read_value (std::istream& is) +{ + return read_cx_fp_value<float> (is); +} - return val; - } - - template <typename T> - std::complex<T> - read_cx_fp_value (std::istream& is) - { - T re = 0.0; - T im = 0.0; - - std::complex<T> cx = 0.0; - - char ch = ' '; - - while (isspace (ch)) - ch = is.get (); - - if (ch == '(') - { - re = read_value<T> (is); - ch = is.get (); - - if (ch == ',') - { - im = read_value<T> (is); - ch = is.get (); +template <typename T> +void +write_value (std::ostream& os, const T& value) +{ + os << value; +} - if (ch == ')') - cx = std::complex<T> (re, im); - else - is.setstate (std::ios::failbit); - } - else if (ch == ')') - cx = re; - else - is.setstate (std::ios::failbit); - } - else - { - is.putback (ch); - cx = read_value<T> (is); - } - - return cx; - } - - // FIXME: Could we use traits and enable_if to avoid duplication in the - // following specializations? +template OCTAVE_API void +write_value<bool> (std::ostream& os, const bool& value); +template OCTAVE_API void +write_value<octave_int8> (std::ostream& os, const octave_int8& value); +template OCTAVE_API void +write_value<octave_int16> (std::ostream& os, const octave_int16& value); +template OCTAVE_API void +write_value<octave_int32> (std::ostream& os, const octave_int32& value); +template OCTAVE_API void +write_value<octave_int64> (std::ostream& os, const octave_int64& value); +template OCTAVE_API void +write_value<octave_uint8> (std::ostream& os, const octave_uint8& value); +template OCTAVE_API void +write_value<octave_uint16> (std::ostream& os, const octave_uint16& value); +template OCTAVE_API void +write_value<octave_uint32> (std::ostream& os, const octave_uint32& value); +template OCTAVE_API void +write_value<octave_uint64> (std::ostream& os, const octave_uint64& value); - template <> OCTAVE_API double read_value (std::istream& is) - { - return read_fp_value<double> (is); - } - - template <> OCTAVE_API Complex read_value (std::istream& is) - { - return read_cx_fp_value<double> (is); - } +// Note: precision is supposed to be managed outside of this function by +// setting stream parameters. - template <> OCTAVE_API float read_value (std::istream& is) - { - return read_fp_value<float> (is); - } - - template <> OCTAVE_API FloatComplex read_value (std::istream& is) - { - return read_cx_fp_value<float> (is); - } - - template <typename T> - void - write_value (std::ostream& os, const T& value) - { +template <> OCTAVE_API void +write_value (std::ostream& os, const double& value) +{ + if (lo_ieee_is_NA (value)) + os << "NA"; + else if (lo_ieee_isnan (value)) + os << "NaN"; + else if (lo_ieee_isinf (value)) + os << (value < 0 ? "-Inf" : "Inf"); + else os << value; - } +} - template OCTAVE_API void - write_value<bool> (std::ostream& os, const bool& value); - template OCTAVE_API void - write_value<octave_int8> (std::ostream& os, const octave_int8& value); - template OCTAVE_API void - write_value<octave_int16> (std::ostream& os, const octave_int16& value); - template OCTAVE_API void - write_value<octave_int32> (std::ostream& os, const octave_int32& value); - template OCTAVE_API void - write_value<octave_int64> (std::ostream& os, const octave_int64& value); - template OCTAVE_API void - write_value<octave_uint8> (std::ostream& os, const octave_uint8& value); - template OCTAVE_API void - write_value<octave_uint16> (std::ostream& os, const octave_uint16& value); - template OCTAVE_API void - write_value<octave_uint32> (std::ostream& os, const octave_uint32& value); - template OCTAVE_API void - write_value<octave_uint64> (std::ostream& os, const octave_uint64& value); +template <> OCTAVE_API void +write_value (std::ostream& os, const Complex& value) +{ + os << '('; + write_value<double> (os, real (value)); + os << ','; + write_value<double> (os, imag (value)); + os << ')'; +} - // Note: precision is supposed to be managed outside of this function by - // setting stream parameters. - - template <> OCTAVE_API void - write_value (std::ostream& os, const double& value) - { - if (lo_ieee_is_NA (value)) - os << "NA"; - else if (lo_ieee_isnan (value)) - os << "NaN"; - else if (lo_ieee_isinf (value)) - os << (value < 0 ? "-Inf" : "Inf"); - else - os << value; - } +// Note: precision is supposed to be managed outside of this function by +// setting stream parameters. - template <> OCTAVE_API void - write_value (std::ostream& os, const Complex& value) - { - os << '('; - write_value<double> (os, real (value)); - os << ','; - write_value<double> (os, imag (value)); - os << ')'; - } - - // Note: precision is supposed to be managed outside of this function by - // setting stream parameters. +template <> OCTAVE_API void +write_value (std::ostream& os, const float& value) +{ + if (lo_ieee_is_NA (value)) + os << "NA"; + else if (lo_ieee_isnan (value)) + os << "NaN"; + else if (lo_ieee_isinf (value)) + os << (value < 0 ? "-Inf" : "Inf"); + else + os << value; +} - template <> OCTAVE_API void - write_value (std::ostream& os, const float& value) - { - if (lo_ieee_is_NA (value)) - os << "NA"; - else if (lo_ieee_isnan (value)) - os << "NaN"; - else if (lo_ieee_isinf (value)) - os << (value < 0 ? "-Inf" : "Inf"); - else - os << value; - } - - template <> OCTAVE_API void - write_value (std::ostream& os, const FloatComplex& value) - { - os << '('; - write_value<float> (os, real (value)); - os << ','; - write_value<float> (os, imag (value)); - os << ')'; - } +template <> OCTAVE_API void +write_value (std::ostream& os, const FloatComplex& value) +{ + os << '('; + write_value<float> (os, real (value)); + os << ','; + write_value<float> (os, imag (value)); + os << ')'; +} OCTAVE_BEGIN_NAMESPACE(math) - bool int_multiply_overflow (int a, int b, int *r) - { - return octave_i_multiply_overflow_wrapper (a, b, r); - } + bool int_multiply_overflow (int a, int b, int *r) + { + return octave_i_multiply_overflow_wrapper (a, b, r); + } - bool int_multiply_overflow (long int a, long int b, long int *r) - { - return octave_li_multiply_overflow_wrapper (a, b, r); - } + bool int_multiply_overflow (long int a, long int b, long int *r) + { + return octave_li_multiply_overflow_wrapper (a, b, r); + } #if defined (OCTAVE_HAVE_LONG_LONG_INT) - bool int_multiply_overflow (long long int a, long long int b, - long long int *r) - { - return octave_lli_multiply_overflow_wrapper (a, b, r); - } + bool int_multiply_overflow (long long int a, long long int b, + long long int *r) + { + return octave_lli_multiply_overflow_wrapper (a, b, r); + } #endif - bool int_multiply_overflow (unsigned int a, unsigned int b, - unsigned int *r) - { - return octave_ui_multiply_overflow_wrapper (a, b, r); - } + bool int_multiply_overflow (unsigned int a, unsigned int b, + unsigned int *r) + { + return octave_ui_multiply_overflow_wrapper (a, b, r); + } - bool int_multiply_overflow (unsigned long int a, unsigned long int b, - unsigned long int *r) - { - return octave_uli_multiply_overflow_wrapper (a, b, r); - } + bool int_multiply_overflow (unsigned long int a, unsigned long int b, + unsigned long int *r) + { + return octave_uli_multiply_overflow_wrapper (a, b, r); + } #if defined (OCTAVE_HAVE_UNSIGNED_LONG_LONG_INT) - bool int_multiply_overflow (unsigned long long int a, - unsigned long long int b, - unsigned long long int *r) - { - return octave_ulli_multiply_overflow_wrapper (a, b, r); - } + bool int_multiply_overflow (unsigned long long int a, + unsigned long long int b, + unsigned long long int *r) + { + return octave_ulli_multiply_overflow_wrapper (a, b, r); + } #endif OCTAVE_END_NAMESPACE(math)