Dot.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2006-2008, 2010 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_DOT_H
#define EIGEN_DOT_H

// IWYU pragma: private
#include "./InternalHeaderCheck.h"

namespace Eigen {

namespace internal {

template <typename Derived, typename Scalar = typename traits<Derived>::Scalar>
struct squared_norm_impl {
  using Real = typename NumTraits<Scalar>::Real;
  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Real run(const Derived& a) {
    Scalar result = a.unaryExpr(squared_norm_functor<Scalar>()).sum();
    return numext::real(result) + numext::imag(result);
  }
};

template <typename Derived>
struct squared_norm_impl<Derived, bool> {
  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool run(const Derived& a) { return a.any(); }
};

}  // end namespace internal

template <typename Derived>
template <typename OtherDerived>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
    typename ScalarBinaryOpTraits<typename internal::traits<Derived>::Scalar,
                                  typename internal::traits<OtherDerived>::Scalar>::ReturnType
    MatrixBase<Derived>::dot(const MatrixBase<OtherDerived>& other) const {
  return internal::dot_impl<Derived, OtherDerived>::run(derived(), other.derived());
}

//---------- implementation of L2 norm and related functions ----------

template <typename Derived>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
MatrixBase<Derived>::squaredNorm() const {
  return internal::squared_norm_impl<Derived>::run(derived());
}

template <typename Derived>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
MatrixBase<Derived>::norm() const {
  return numext::sqrt(squaredNorm());
}

template <typename Derived>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const typename MatrixBase<Derived>::PlainObject MatrixBase<Derived>::normalized()
    const {
  typedef typename internal::nested_eval<Derived, 2>::type Nested_;
  Nested_ n(derived());
  RealScalar z = n.squaredNorm();
  // NOTE: after extensive benchmarking, this conditional does not impact performance, at least on recent x86 CPU
  if (z > RealScalar(0))
    return n / numext::sqrt(z);
  else
    return n;
}

template <typename Derived>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void MatrixBase<Derived>::normalize() {
  RealScalar z = squaredNorm();
  // NOTE: after extensive benchmarking, this conditional does not impact performance, at least on recent x86 CPU
  if (z > RealScalar(0)) derived() /= numext::sqrt(z);
}

template <typename Derived>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const typename MatrixBase<Derived>::PlainObject
MatrixBase<Derived>::stableNormalized() const {
  typedef typename internal::nested_eval<Derived, 3>::type Nested_;
  Nested_ n(derived());
  RealScalar w = n.cwiseAbs().maxCoeff();
  RealScalar z = (n / w).squaredNorm();
  if (z > RealScalar(0))
    return n / (numext::sqrt(z) * w);
  else
    return n;
}

template <typename Derived>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void MatrixBase<Derived>::stableNormalize() {
  RealScalar w = cwiseAbs().maxCoeff();
  RealScalar z = (derived() / w).squaredNorm();
  if (z > RealScalar(0)) derived() /= numext::sqrt(z) * w;
}

//---------- implementation of other norms ----------

namespace internal {

template <typename Derived, int p>
struct lpNorm_selector {
  typedef typename NumTraits<typename traits<Derived>::Scalar>::Real RealScalar;
  EIGEN_DEVICE_FUNC static inline RealScalar run(const MatrixBase<Derived>& m) {
    EIGEN_USING_STD(pow)
    return pow(m.cwiseAbs().array().pow(p).sum(), RealScalar(1) / p);
  }
};

template <typename Derived>
struct lpNorm_selector<Derived, 1> {
  EIGEN_DEVICE_FUNC static inline typename NumTraits<typename traits<Derived>::Scalar>::Real run(
      const MatrixBase<Derived>& m) {
    return m.cwiseAbs().sum();
  }
};

template <typename Derived>
struct lpNorm_selector<Derived, 2> {
  EIGEN_DEVICE_FUNC static inline typename NumTraits<typename traits<Derived>::Scalar>::Real run(
      const MatrixBase<Derived>& m) {
    return m.norm();
  }
};

template <typename Derived>
struct lpNorm_selector<Derived, Infinity> {
  typedef typename NumTraits<typename traits<Derived>::Scalar>::Real RealScalar;
  EIGEN_DEVICE_FUNC static inline RealScalar run(const MatrixBase<Derived>& m) {
    if (Derived::SizeAtCompileTime == 0 || (Derived::SizeAtCompileTime == Dynamic && m.size() == 0))
      return RealScalar(0);
    return m.cwiseAbs().maxCoeff();
  }
};

}  // end namespace internal

template <typename Derived>
template <int p>
#ifndef EIGEN_PARSED_BY_DOXYGEN
EIGEN_DEVICE_FUNC inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
#else
EIGEN_DEVICE_FUNC MatrixBase<Derived>::RealScalar
#endif
MatrixBase<Derived>::lpNorm() const {
  return internal::lpNorm_selector<Derived, p>::run(*this);
}

//---------- implementation of isOrthogonal / isUnitary ----------

template <typename Derived>
template <typename OtherDerived>
bool MatrixBase<Derived>::isOrthogonal(const MatrixBase<OtherDerived>& other, const RealScalar& prec) const {
  typename internal::nested_eval<Derived, 2>::type nested(derived());
  typename internal::nested_eval<OtherDerived, 2>::type otherNested(other.derived());
  return numext::abs2(nested.dot(otherNested)) <= prec * prec * nested.squaredNorm() * otherNested.squaredNorm();
}

template <typename Derived>
bool MatrixBase<Derived>::isUnitary(const RealScalar& prec) const {
  typename internal::nested_eval<Derived, 1>::type self(derived());
  for (Index i = 0; i < cols(); ++i) {
    if (!internal::isApprox(self.col(i).squaredNorm(), static_cast<RealScalar>(1), prec)) return false;
    for (Index j = 0; j < i; ++j)
      if (!internal::isMuchSmallerThan(self.col(i).dot(self.col(j)), static_cast<Scalar>(1), prec)) return false;
  }
  return true;
}

}  // end namespace Eigen

#endif  // EIGEN_DOT_H