Format: Bloom subset of https://www.debian.org/doc/packaging-manuals/copyright-format/1.0/
Upstream-Name: gtsam

Files: See file headers in repository for details
Copyright: See package copyright in source code for details
License: BSD-3-Clause
 See repository for full license text

Files: See file headers in repository for details
Copyright: See package copyright in source code for details
License: BSD-3-Clause
 CCOLAMD: constrained column approximate minimum degree ordering
 Copyright (C) 2005-2016, Univ. of Florida.  Authors: Timothy A. Davis,
 Sivasankaran Rajamanickam, and Stefan Larimore.  Closely based on COLAMD by
 Davis, Stefan Larimore, in collaboration with Esmond Ng, and John Gilbert.
 http://www.suitesparse.com
 .
 --------------------------------------------------------------------------------
 .
 CCOLAMD license: BSD 3-clause:
 .
     Redistribution and use in source and binary forms, with or without
     modification, are permitted provided that the following conditions are met:
         * Redistributions of source code must retain the above copyright
           notice, this list of conditions and the following disclaimer.
         * Redistributions in binary form must reproduce the above copyright
           notice, this list of conditions and the following disclaimer in the
           documentation and/or other materials provided with the distribution.
         * Neither the name of the organizations to which the authors are
           affiliated, nor the names of its contributors may be used to endorse
           or promote products derived from this software without specific prior
           written permission.
 .
     THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
     AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
     IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
     ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY
     DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
     (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
     SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
     CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
     LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
     OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
     DAMAGE.

Files: See file headers in repository for details
Copyright: See package copyright in source code for details
License: BSD-3-Clause
 // Ceres Solver - A fast non-linear least squares minimizer
 // Copyright 2010, 2011, 2012 Google Inc. All rights reserved.
 // http://code.google.com/p/ceres-solver/
 //
 // Redistribution and use in source and binary forms, with or without
 // modification, are permitted provided that the following conditions are met:
 //
 // * Redistributions of source code must retain the above copyright notice,
 //   this list of conditions and the following disclaimer.
 // * Redistributions in binary form must reproduce the above copyright notice,
 //   this list of conditions and the following disclaimer in the documentation
 //   and/or other materials provided with the distribution.
 // * Neither the name of Google Inc. nor the names of its contributors may be
 //   used to endorse or promote products derived from this software without
 //   specific prior written permission.
 //
 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
 // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 // POSSIBILITY OF SUCH DAMAGE.
 //
 // Author: keir@google.com (Keir Mierle)
 //
 // A simple implementation of N-dimensional dual numbers, for automatically
 // computing exact derivatives of functions.
 //
 // While a complete treatment of the mechanics of automatic differentation is
 // beyond the scope of this header (see
 // http://en.wikipedia.org/wiki/Automatic_differentiation for details), the
 // basic idea is to extend normal arithmetic with an extra element, "e," often
 // denoted with the greek symbol epsilon, such that e != 0 but e^2 = 0. Dual
 // numbers are extensions of the real numbers analogous to complex numbers:
 // whereas complex numbers augment the reals by introducing an imaginary unit i
 // such that i^2 = -1, dual numbers introduce an "infinitesimal" unit e such
 // that e^2 = 0. Dual numbers have two components: the "real" component and the
 // "infinitesimal" component, generally written as x + y*e. Surprisingly, this
 // leads to a convenient method for computing exact derivatives without needing
 // to manipulate complicated symbolic expressions.
 //
 // For example, consider the function
 //
 //   f(x) = x^2 ,
 //
 // evaluated at 10. Using normal arithmetic, f(10) = 100, and df/dx(10) = 20.
 // Next, augument 10 with an infinitesimal to get:
 //
 //   f(10 + e) = (10 + e)^2
 //             = 100 + 2 * 10 * e + e^2
 //             = 100 + 20 * e       -+-
 //                     --            |
 //                     |             +--- This is zero, since e^2 = 0
 //                     |
 //                     +----------------- This is df/dx!
 //
 // Note that the derivative of f with respect to x is simply the infinitesimal
 // component of the value of f(x + e). So, in order to take the derivative of
 // any function, it is only necessary to replace the numeric "object" used in
 // the function with one extended with infinitesimals. The class Jet, defined in
 // this header, is one such example of this, where substitution is done with
 // templates.
 //
 // To handle derivatives of functions taking multiple arguments, different
 // infinitesimals are used, one for each variable to take the derivative of. For
 // example, consider a scalar function of two scalar parameters x and y:
 //
 //   f(x, y) = x^2 + x * y
 //
 // Following the technique above, to compute the derivatives df/dx and df/dy for
 // f(1, 3) involves doing two evaluations of f, the first time replacing x with
 // x + e, the second time replacing y with y + e.
 //
 // For df/dx:
 //
 //   f(1 + e, y) = (1 + e)^2 + (1 + e) * 3
 //               = 1 + 2 * e + 3 + 3 * e
 //               = 4 + 5 * e
 //
 //               --> df/dx = 5
 //
 // For df/dy:
 //
 //   f(1, 3 + e) = 1^2 + 1 * (3 + e)
 //               = 1 + 3 + e
 //               = 4 + e
 //
 //               --> df/dy = 1
 //
 // To take the gradient of f with the implementation of dual numbers ("jets") in
 // this file, it is necessary to create a single jet type which has components
 // for the derivative in x and y, and passing them to a templated version of f:
 //
 //   template<typename T>
 //   T f(const T &x, const T &y) {
 //     return x * x + x * y;
 //   }
 //
 //   // The "2" means there should be 2 dual number components.
 //   Jet<double, 2> x(0);  // Pick the 0th dual number for x.
 //   Jet<double, 2> y(1);  // Pick the 1st dual number for y.
 //   Jet<double, 2> z = f(x, y);
 //
 //   LOG(INFO) << "df/dx = " << z.a[0]
 //             << "df/dy = " << z.a[1];
 //
 // Most users should not use Jet objects directly; a wrapper around Jet objects,
 // which makes computing the derivative, gradient, or jacobian of templated
 // functors simple, is in autodiff.h. Even autodiff.h should not be used
 // directly; instead autodiff_cost_function.h is typically the file of interest.
 //
 // For the more mathematically inclined, this file implements first-order
 // "jets". A 1st order jet is an element of the ring
 //
 //   T[N] = T[t_1, ..., t_N] / (t_1, ..., t_N)^2
 //
 // which essentially means that each jet consists of a "scalar" value 'a' from T
 // and a 1st order perturbation vector 'v' of length N:
 //
 //   x = a + \sum_i v[i] t_i
 //
 // A shorthand is to write an element as x = a + u, where u is the pertubation.
 // Then, the main point about the arithmetic of jets is that the product of
 // perturbations is zero:
 //
 //   (a + u) * (b + v) = ab + av + bu + uv
 //                     = ab + (av + bu) + 0
 //
 // which is what operator* implements below. Addition is simpler:
 //
 //   (a + u) + (b + v) = (a + b) + (u + v).
 //
 // The only remaining question is how to evaluate the function of a jet, for
 // which we use the chain rule:
 //
 //   f(a + u) = f(a) + f'(a) u
 //
 // where f'(a) is the (scalar) derivative of f at a.
 //
 // By pushing these things through sufficiently and suitably templated
 // functions, we can do automatic differentiation. Just be sure to turn on
 // function inlining and common-subexpression elimination, or it will be very
 // slow!
 //
 // WARNING: Most Ceres users should not directly include this file or know the
 // details of how jets work. Instead the suggested method for automatic
 // derivatives is to use autodiff_cost_function.h, which is a wrapper around
 // both jets.h and autodiff.h to make taking derivatives of cost functions for
 // use in Ceres easier.
 .
 #ifndef CERES_PUBLIC_JET_H_
 #define CERES_PUBLIC_JET_H_
 .
 #include <cmath>
 #include <iosfwd>
 #include <iostream>  // NOLINT
 #include <limits>
 #include <string>
 .
 #include <gtsam/3rdparty/ceres/eigen.h>
 #include <gtsam/3rdparty/ceres/fpclassify.h>
 .
 namespace ceres {
 .
 template <typename T, int N>
 struct Jet {
   enum { DIMENSION = N };
 .
   // Default-construct "a" because otherwise this can lead to false errors about
   // uninitialized uses when other classes relying on default constructed T
   // (where T is a Jet<T, N>). This usually only happens in opt mode. Note that
   // the C++ standard mandates that e.g. default constructed doubles are
   // initialized to 0.0; see sections 8.5 of the C++03 standard.
   Jet() : a() {
     v.setZero();
   }
 .
   // Constructor from scalar: a + 0.
   explicit Jet(const T& value) {
     a = value;
     v.setZero();
   }
 .
   // Constructor from scalar plus variable: a + t_i.
   Jet(const T& value, int k) {
     a = value;
     v.setZero();
     v[k] = T(1.0);
   }
 .
   // Constructor from scalar and vector part
   // The use of Eigen::DenseBase allows Eigen expressions
   // to be passed in without being fully evaluated until
   // they are assigned to v
   template<typename Derived>
   EIGEN_STRONG_INLINE Jet(const T& a, const Eigen::DenseBase<Derived> &v)
       : a(a), v(v) {
   }
 .
   // Compound operators
   Jet<T, N>& operator+=(const Jet<T, N> &y) {
     *this = *this + y;
     return *this;
   }
 .
   Jet<T, N>& operator-=(const Jet<T, N> &y) {
     *this = *this - y;
     return *this;
   }
 .
   Jet<T, N>& operator*=(const Jet<T, N> &y) {
     *this = *this * y;
     return *this;
   }
 .
   Jet<T, N>& operator/=(const Jet<T, N> &y) {
     *this = *this / y;
     return *this;
   }
 .
   // The scalar part.
   T a;
 .
   // The infinitesimal part.
   //
   // Note the Eigen::DontAlign bit is needed here because this object
   // gets allocated on the stack and as part of other arrays and
   // structs. Forcing the right alignment there is the source of much
   // pain and suffering. Even if that works, passing Jets around to
   // functions by value has problems because the C++ ABI does not
   // guarantee alignment for function arguments.
   //
   // Setting the DontAlign bit prevents Eigen from using SSE for the
   // various operations on Jets. This is a small performance penalty
   // since the AutoDiff code will still expose much of the code as
   // statically sized loops to the compiler. But given the subtle
   // issues that arise due to alignment, especially when dealing with
   // multiple platforms, it seems to be a trade off worth making.
   Eigen::Matrix<T, N, 1, Eigen::DontAlign> v;
 };
 .
 // Unary +
 template<typename T, int N> inline
 Jet<T, N> const& operator+(const Jet<T, N>& f) {
   return f;
 }
 .
 // TODO(keir): Try adding __attribute__((always_inline)) to these functions to
 // see if it causes a performance increase.
 .
 // Unary -
 template<typename T, int N> inline
 Jet<T, N> operator-(const Jet<T, N>&f) {
   return Jet<T, N>(-f.a, -f.v);
 }
 .
 // Binary +
 template<typename T, int N> inline
 Jet<T, N> operator+(const Jet<T, N>& f,
                     const Jet<T, N>& g) {
   return Jet<T, N>(f.a + g.a, f.v + g.v);
 }
 .
 // Binary + with a scalar: x + s
 template<typename T, int N> inline
 Jet<T, N> operator+(const Jet<T, N>& f, T s) {
   return Jet<T, N>(f.a + s, f.v);
 }
 .
 // Binary + with a scalar: s + x
 template<typename T, int N> inline
 Jet<T, N> operator+(T s, const Jet<T, N>& f) {
   return Jet<T, N>(f.a + s, f.v);
 }
 .
 // Binary -
 template<typename T, int N> inline
 Jet<T, N> operator-(const Jet<T, N>& f,
                     const Jet<T, N>& g) {
   return Jet<T, N>(f.a - g.a, f.v - g.v);
 }
 .
 // Binary - with a scalar: x - s
 template<typename T, int N> inline
 Jet<T, N> operator-(const Jet<T, N>& f, T s) {
   return Jet<T, N>(f.a - s, f.v);
 }
 .
 // Binary - with a scalar: s - x
 template<typename T, int N> inline
 Jet<T, N> operator-(T s, const Jet<T, N>& f) {
   return Jet<T, N>(s - f.a, -f.v);
 }
 .
 // Binary *
 template<typename T, int N> inline
 Jet<T, N> operator*(const Jet<T, N>& f,
                     const Jet<T, N>& g) {
   return Jet<T, N>(f.a * g.a, f.a * g.v + f.v * g.a);
 }
 .
 // Binary * with a scalar: x * s
 template<typename T, int N> inline
 Jet<T, N> operator*(const Jet<T, N>& f, T s) {
   return Jet<T, N>(f.a * s, f.v * s);
 }
 .
 // Binary * with a scalar: s * x
 template<typename T, int N> inline
 Jet<T, N> operator*(T s, const Jet<T, N>& f) {
   return Jet<T, N>(f.a * s, f.v * s);
 }
 .
 // Binary /
 template<typename T, int N> inline
 Jet<T, N> operator/(const Jet<T, N>& f,
                     const Jet<T, N>& g) {
   // This uses:
   //
   //   a + u   (a + u)(b - v)   (a + u)(b - v)
   //   ----- = -------------- = --------------
   //   b + v   (b + v)(b - v)        b^2
   //
   // which holds because v*v = 0.
   const T g_a_inverse = T(1.0) / g.a;
   const T f_a_by_g_a = f.a * g_a_inverse;
   return Jet<T, N>(f.a * g_a_inverse, (f.v - f_a_by_g_a * g.v) * g_a_inverse);
 }
 .
 // Binary / with a scalar: s / x
 template<typename T, int N> inline
 Jet<T, N> operator/(T s, const Jet<T, N>& g) {
   const T minus_s_g_a_inverse2 = -s / (g.a * g.a);
   return Jet<T, N>(s / g.a, g.v * minus_s_g_a_inverse2);
 }
 .
 // Binary / with a scalar: x / s
 template<typename T, int N> inline
 Jet<T, N> operator/(const Jet<T, N>& f, T s) {
   const T s_inverse = 1.0 / s;
   return Jet<T, N>(f.a * s_inverse, f.v * s_inverse);
 }
 .
 // Binary comparison operators for both scalars and jets.
 #define CERES_DEFINE_JET_COMPARISON_OPERATOR(op) \
 template<typename T, int N> inline \
 bool operator op(const Jet<T, N>& f, const Jet<T, N>& g) { \
   return f.a op g.a; \
 } \
 template<typename T, int N> inline \
 bool operator op(const T& s, const Jet<T, N>& g) { \
   return s op g.a; \
 } \
 template<typename T, int N> inline \
 bool operator op(const Jet<T, N>& f, const T& s) { \
   return f.a op s; \
 }
 CERES_DEFINE_JET_COMPARISON_OPERATOR( <  )  // NOLINT
 CERES_DEFINE_JET_COMPARISON_OPERATOR( <= )  // NOLINT
 CERES_DEFINE_JET_COMPARISON_OPERATOR( >  )  // NOLINT
 CERES_DEFINE_JET_COMPARISON_OPERATOR( >= )  // NOLINT
 CERES_DEFINE_JET_COMPARISON_OPERATOR( == )  // NOLINT
 CERES_DEFINE_JET_COMPARISON_OPERATOR( != )  // NOLINT
 #undef CERES_DEFINE_JET_COMPARISON_OPERATOR
 .
 // Pull some functions from namespace std.
 //
 // This is necessary because we want to use the same name (e.g. 'sqrt') for
 // double-valued and Jet-valued functions, but we are not allowed to put
 // Jet-valued functions inside namespace std.
 //
 // TODO(keir): Switch to "using".
 inline double abs     (double x) { return std::abs(x);      }
 inline double log     (double x) { return std::log(x);      }
 inline double exp     (double x) { return std::exp(x);      }
 inline double sqrt    (double x) { return std::sqrt(x);     }
 inline double cos     (double x) { return std::cos(x);      }
 inline double acos    (double x) { return std::acos(x);     }
 inline double sin     (double x) { return std::sin(x);      }
 inline double asin    (double x) { return std::asin(x);     }
 inline double tan     (double x) { return std::tan(x);      }
 inline double atan    (double x) { return std::atan(x);     }
 inline double sinh    (double x) { return std::sinh(x);     }
 inline double cosh    (double x) { return std::cosh(x);     }
 inline double tanh    (double x) { return std::tanh(x);     }
 inline double pow  (double x, double y) { return std::pow(x, y);   }
 inline double atan2(double y, double x) { return std::atan2(y, x); }
 .
 // In general, f(a + h) ~= f(a) + f'(a) h, via the chain rule.
 .
 // abs(x + h) ~= x + h or -(x + h)
 template <typename T, int N> inline
 Jet<T, N> abs(const Jet<T, N>& f) {
   return f.a < T(0.0) ? -f : f;
 }
 .
 // log(a + h) ~= log(a) + h / a
 template <typename T, int N> inline
 Jet<T, N> log(const Jet<T, N>& f) {
   const T a_inverse = T(1.0) / f.a;
   return Jet<T, N>(log(f.a), f.v * a_inverse);
 }
 .
 // exp(a + h) ~= exp(a) + exp(a) h
 template <typename T, int N> inline
 Jet<T, N> exp(const Jet<T, N>& f) {
   const T tmp = exp(f.a);
   return Jet<T, N>(tmp, tmp * f.v);
 }
 .
 // sqrt(a + h) ~= sqrt(a) + h / (2 sqrt(a))
 template <typename T, int N> inline
 Jet<T, N> sqrt(const Jet<T, N>& f) {
   const T tmp = sqrt(f.a);
   const T two_a_inverse = T(1.0) / (T(2.0) * tmp);
   return Jet<T, N>(tmp, f.v * two_a_inverse);
 }
 .
 // cos(a + h) ~= cos(a) - sin(a) h
 template <typename T, int N> inline
 Jet<T, N> cos(const Jet<T, N>& f) {
   return Jet<T, N>(cos(f.a), - sin(f.a) * f.v);
 }
 .
 // acos(a + h) ~= acos(a) - 1 / sqrt(1 - a^2) h
 template <typename T, int N> inline
 Jet<T, N> acos(const Jet<T, N>& f) {
   const T tmp = - T(1.0) / sqrt(T(1.0) - f.a * f.a);
   return Jet<T, N>(acos(f.a), tmp * f.v);
 }
 .
 // sin(a + h) ~= sin(a) + cos(a) h
 template <typename T, int N> inline
 Jet<T, N> sin(const Jet<T, N>& f) {
   return Jet<T, N>(sin(f.a), cos(f.a) * f.v);
 }
 .
 // asin(a + h) ~= asin(a) + 1 / sqrt(1 - a^2) h
 template <typename T, int N> inline
 Jet<T, N> asin(const Jet<T, N>& f) {
   const T tmp = T(1.0) / sqrt(T(1.0) - f.a * f.a);
   return Jet<T, N>(asin(f.a), tmp * f.v);
 }
 .
 // tan(a + h) ~= tan(a) + (1 + tan(a)^2) h
 template <typename T, int N> inline
 Jet<T, N> tan(const Jet<T, N>& f) {
   const T tan_a = tan(f.a);
   const T tmp = T(1.0) + tan_a * tan_a;
   return Jet<T, N>(tan_a, tmp * f.v);
 }
 .
 // atan(a + h) ~= atan(a) + 1 / (1 + a^2) h
 template <typename T, int N> inline
 Jet<T, N> atan(const Jet<T, N>& f) {
   const T tmp = T(1.0) / (T(1.0) + f.a * f.a);
   return Jet<T, N>(atan(f.a), tmp * f.v);
 }
 .
 // sinh(a + h) ~= sinh(a) + cosh(a) h
 template <typename T, int N> inline
 Jet<T, N> sinh(const Jet<T, N>& f) {
   return Jet<T, N>(sinh(f.a), cosh(f.a) * f.v);
 }
 .
 // cosh(a + h) ~= cosh(a) + sinh(a) h
 template <typename T, int N> inline
 Jet<T, N> cosh(const Jet<T, N>& f) {
   return Jet<T, N>(cosh(f.a), sinh(f.a) * f.v);
 }
 .
 // tanh(a + h) ~= tanh(a) + (1 - tanh(a)^2) h
 template <typename T, int N> inline
 Jet<T, N> tanh(const Jet<T, N>& f) {
   const T tanh_a = tanh(f.a);
   const T tmp = T(1.0) - tanh_a * tanh_a;
   return Jet<T, N>(tanh_a, tmp * f.v);
 }
 .
 // Jet Classification. It is not clear what the appropriate semantics are for
 // these classifications. This picks that IsFinite and isnormal are "all"
 // operations, i.e. all elements of the jet must be finite for the jet itself
 // to be finite (or normal). For IsNaN and IsInfinite, the answer is less
 // clear. This takes a "any" approach for IsNaN and IsInfinite such that if any
 // part of a jet is nan or inf, then the entire jet is nan or inf. This leads
 // to strange situations like a jet can be both IsInfinite and IsNaN, but in
 // practice the "any" semantics are the most useful for e.g. checking that
 // derivatives are sane.
 .
 // The jet is finite if all parts of the jet are finite.
 template <typename T, int N> inline
 bool IsFinite(const Jet<T, N>& f) {
   if (!IsFinite(f.a)) {
     return false;
   }
   for (int i = 0; i < N; ++i) {
     if (!IsFinite(f.v[i])) {
       return false;
     }
   }
   return true;
 }
 .
 // The jet is infinite if any part of the jet is infinite.
 template <typename T, int N> inline
 bool IsInfinite(const Jet<T, N>& f) {
   if (IsInfinite(f.a)) {
     return true;
   }
   for (int i = 0; i < N; i++) {
     if (IsInfinite(f.v[i])) {
       return true;
     }
   }
   return false;
 }
 .
 // The jet is NaN if any part of the jet is NaN.
 template <typename T, int N> inline
 bool IsNaN(const Jet<T, N>& f) {
   if (IsNaN(f.a)) {
     return true;
   }
   for (int i = 0; i < N; ++i) {
     if (IsNaN(f.v[i])) {
       return true;
     }
   }
   return false;
 }
 .
 // The jet is normal if all parts of the jet are normal.
 template <typename T, int N> inline
 bool IsNormal(const Jet<T, N>& f) {
   if (!IsNormal(f.a)) {
     return false;
   }
   for (int i = 0; i < N; ++i) {
     if (!IsNormal(f.v[i])) {
       return false;
     }
   }
   return true;
 }
 .
 // atan2(b + db, a + da) ~= atan2(b, a) + (- b da + a db) / (a^2 + b^2)
 //
 // In words: the rate of change of theta is 1/r times the rate of
 // change of (x, y) in the positive angular direction.
 template <typename T, int N> inline
 Jet<T, N> atan2(const Jet<T, N>& g, const Jet<T, N>& f) {
   // Note order of arguments:
   //
   //   f = a + da
   //   g = b + db
 .
   T const tmp = T(1.0) / (f.a * f.a + g.a * g.a);
   return Jet<T, N>(atan2(g.a, f.a), tmp * (- g.a * f.v + f.a * g.v));
 }
 .
 .
 // pow -- base is a differentiable function, exponent is a constant.
 // (a+da)^p ~= a^p + p*a^(p-1) da
 template <typename T, int N> inline
 Jet<T, N> pow(const Jet<T, N>& f, double g) {
   T const tmp = g * pow(f.a, g - T(1.0));
   return Jet<T, N>(pow(f.a, g), tmp * f.v);
 }
 .
 // pow -- base is a constant, exponent is a differentiable function.
 // (a)^(p+dp) ~= a^p + a^p log(a) dp
 template <typename T, int N> inline
 Jet<T, N> pow(double f, const Jet<T, N>& g) {
   T const tmp = pow(f, g.a);
   return Jet<T, N>(tmp, log(f) * tmp * g.v);
 }
 .
 .
 // pow -- both base and exponent are differentiable functions.
 // (a+da)^(b+db) ~= a^b + b * a^(b-1) da + a^b log(a) * db
 template <typename T, int N> inline
 Jet<T, N> pow(const Jet<T, N>& f, const Jet<T, N>& g) {
   T const tmp1 = pow(f.a, g.a);
   T const tmp2 = g.a * pow(f.a, g.a - T(1.0));
   T const tmp3 = tmp1 * log(f.a);
 .
   return Jet<T, N>(tmp1, tmp2 * f.v + tmp3 * g.v);
 }
 .
 // Define the helper functions Eigen needs to embed Jet types.
 //
 // NOTE(keir): machine_epsilon() and precision() are missing, because they don't
 // work with nested template types (e.g. where the scalar is itself templated).
 // Among other things, this means that decompositions of Jet's does not work,
 // for example
 //
 //   Matrix<Jet<T, N> ... > A, x, b;
 //   ...
 //   A.solve(b, &x)
 //
 // does not work and will fail with a strange compiler error.
 //
 // TODO(keir): This is an Eigen 2.0 limitation that is lifted in 3.0. When we
 // switch to 3.0, also add the rest of the specialization functionality.
 template<typename T, int N> inline const Jet<T, N>& ei_conj(const Jet<T, N>& x) { return x;              }  // NOLINT
 template<typename T, int N> inline const Jet<T, N>& ei_real(const Jet<T, N>& x) { return x;              }  // NOLINT
 template<typename T, int N> inline       Jet<T, N>  ei_imag(const Jet<T, N>&  ) { return Jet<T, N>(0.0); }  // NOLINT
 template<typename T, int N> inline       Jet<T, N>  ei_abs (const Jet<T, N>& x) { return fabs(x);        }  // NOLINT
 template<typename T, int N> inline       Jet<T, N>  ei_abs2(const Jet<T, N>& x) { return x * x;          }  // NOLINT
 template<typename T, int N> inline       Jet<T, N>  ei_sqrt(const Jet<T, N>& x) { return sqrt(x);        }  // NOLINT
 template<typename T, int N> inline       Jet<T, N>  ei_exp (const Jet<T, N>& x) { return exp(x);         }  // NOLINT
 template<typename T, int N> inline       Jet<T, N>  ei_log (const Jet<T, N>& x) { return log(x);         }  // NOLINT
 template<typename T, int N> inline       Jet<T, N>  ei_sin (const Jet<T, N>& x) { return sin(x);         }  // NOLINT
 template<typename T, int N> inline       Jet<T, N>  ei_cos (const Jet<T, N>& x) { return cos(x);         }  // NOLINT
 template<typename T, int N> inline       Jet<T, N>  ei_tan (const Jet<T, N>& x) { return tan(x);         }  // NOLINT
 template<typename T, int N> inline       Jet<T, N>  ei_atan(const Jet<T, N>& x) { return atan(x);        }  // NOLINT
 template<typename T, int N> inline       Jet<T, N>  ei_sinh(const Jet<T, N>& x) { return sinh(x);        }  // NOLINT
 template<typename T, int N> inline       Jet<T, N>  ei_cosh(const Jet<T, N>& x) { return cosh(x);        }  // NOLINT
 template<typename T, int N> inline       Jet<T, N>  ei_tanh(const Jet<T, N>& x) { return tanh(x);        }  // NOLINT
 template<typename T, int N> inline       Jet<T, N>  ei_pow (const Jet<T, N>& x, Jet<T, N> y) { return pow(x, y); }  // NOLINT
 .
 // Note: This has to be in the ceres namespace for argument dependent lookup to
 // function correctly. Otherwise statements like CHECK_LE(x, 2.0) fail with
 // strange compile errors.
 template <typename T, int N>
 inline std::ostream &operator<<(std::ostream &s, const Jet<T, N>& z) {
   return s << "[" << z.a << " ; " << z.v.transpose() << "]";
 }
 .
 }  // namespace ceres
 .
 namespace Eigen {
 .
 // Creating a specialization of NumTraits enables placing Jet objects inside
 // Eigen arrays, getting all the goodness of Eigen combined with autodiff.
 template<typename T, int N>
 struct NumTraits<ceres::Jet<T, N> > {
   typedef ceres::Jet<T, N> Real;
   typedef ceres::Jet<T, N> NonInteger;
   typedef ceres::Jet<T, N> Nested;
 .
   static typename ceres::Jet<T, N> dummy_precision() {
     return ceres::Jet<T, N>(1e-12);
   }
 .
   static inline Real epsilon() {
     return Real(std::numeric_limits<T>::epsilon());
   }
 .
   enum {
     IsComplex = 0,
     IsInteger = 0,
     IsSigned,
     ReadCost = 1,
     AddCost = 1,
     // For Jet types, multiplication is more expensive than addition.
     MulCost = 3,
     HasFloatingPoint = 1,
     RequireInitialization = 1
   };
 };
 .
 }  // namespace Eigen
 .
 #endif  // CERES_PUBLIC_JET_H_

Files: See file headers in repository for details
Copyright: See package copyright in source code for details
License: MPL-2.0
 Eigen is primarily MPL2 licensed. See COPYING.MPL2 and these links:
   http://www.mozilla.org/MPL/2.0/
   http://www.mozilla.org/MPL/2.0/FAQ.html
 .
 Some files contain third-party code under BSD or LGPL licenses, whence the other
 COPYING.* files here.
 .
 All the LGPL code is either LGPL 2.1-only, or LGPL 2.1-or-later.
 For this reason, the COPYING.LGPL file contains the LGPL 2.1 text.
 .
 If you want to guarantee that the Eigen code that you are #including is licensed
 under the MPL2 and possibly more permissive licenses (like BSD), #define this
 preprocessor symbol:
   EIGEN_MPL2_ONLY
 For example, with most compilers, you could add this to your project CXXFLAGS:
   -DEIGEN_MPL2_ONLY
 This will cause a compilation error to be generated if you #include any code that is
 LGPL licensed.

Files: See file headers in repository for details
Copyright: See package copyright in source code for details
License: MIT
 The MIT License (MIT); this license applies to GeographicLib,
 versions 1.12 and later.
 .
 Copyright (c) 2008-2017, Charles Karney
 .
 Permission is hereby granted, free of charge, to any person
 obtaining a copy of this software and associated documentation
 files (the "Software"), to deal in the Software without
 restriction, including without limitation the rights to use, copy,
 modify, merge, publish, distribute, sublicense, and/or sell copies
 of the Software, and to permit persons to whom the Software is
 furnished to do so, subject to the following conditions:
 .
 The above copyright notice and this permission notice shall be
 included in all copies or substantial portions of the Software.
 .
 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
 NONINFRINGEMENT.  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
 HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
 WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 DEALINGS IN THE SOFTWARE.

Files: See file headers in repository for details
Copyright: See package copyright in source code for details
License: Apache-2.0
 .
 Copyright & License Notice
 ---------------------------
 .
 Copyright 1995-2013, Regents of the University of Minnesota
 .
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at
 .
 http://www.apache.org/licenses/LICENSE-2.0
 .
 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or 
 implied. See the License for the specific language governing 
 permissions and limitations under the License.

Files: See file headers in repository for details
Copyright: See package copyright in source code for details
License: MPL-2.0
 // Copyright (C) 2016-2025 Yixuan Qiu <yixuan.qiu@cos.name>
 //
 // 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 https://mozilla.org/MPL/2.0/.
 .
 #ifndef SPECTRA_GEN_EIGS_SOLVER_H
 #define SPECTRA_GEN_EIGS_SOLVER_H
 .
 #include <Eigen/Core>
 .
 #include "GenEigsBase.h"
 #include "Util/SelectionRule.h"
 #include "MatOp/DenseGenMatProd.h"
 .
 namespace Spectra {
 .
 ///
 /// \ingroup EigenSolver
 ///
 /// This class implements the eigen solver for general real matrices, i.e.,
 /// to solve \f$Ax=\lambda x\f$ for a possibly non-symmetric \f$A\f$ matrix.
 ///
 /// Most of the background information documented in the SymEigsSolver class
 /// also applies to the GenEigsSolver class here, except that the eigenvalues
 /// and eigenvectors of a general matrix can now be complex-valued.
 ///
 /// \tparam OpType  The name of the matrix operation class. Users could either
 ///                 use the wrapper classes such as DenseGenMatProd and
 ///                 SparseGenMatProd, or define their own that implements the type
 ///                 definition `Scalar` and all the public member functions as in
 ///                 DenseGenMatProd.
 ///
 /// An example that illustrates the usage of GenEigsSolver is give below:
 ///
 /// \code{.cpp}
 /// #include <Eigen/Core>
 /// #include <Spectra/GenEigsSolver.h>
 /// // <Spectra/MatOp/DenseGenMatProd.h> is implicitly included
 /// #include <iostream>
 ///
 /// using namespace Spectra;
 ///
 /// int main()
 /// {
 ///     // We are going to calculate the eigenvalues of M
 ///     Eigen::MatrixXd M = Eigen::MatrixXd::Random(10, 10);
 ///
 ///     // Construct matrix operation object using the wrapper class
 ///     DenseGenMatProd<double> op(M);
 ///
 ///     // Construct eigen solver object, requesting the largest
 ///     // (in magnitude, or norm) three eigenvalues
 ///     GenEigsSolver<DenseGenMatProd<double>> eigs(op, 3, 6);
 ///
 ///     // Initialize and compute
 ///     eigs.init();
 ///     int nconv = eigs.compute(SortRule::LargestMagn);
 ///
 ///     // Retrieve results
 ///     Eigen::VectorXcd evalues;
 ///     if (eigs.info() == CompInfo::Successful)
 ///         evalues = eigs.eigenvalues();
 ///
 ///     std::cout << "Eigenvalues found:\n" << evalues << std::endl;
 ///
 ///     return 0;
 /// }
 /// \endcode
 ///
 /// And also an example for sparse matrices:
 ///
 /// \code{.cpp}
 /// #include <Eigen/Core>
 /// #include <Eigen/SparseCore>
 /// #include <Spectra/GenEigsSolver.h>
 /// #include <Spectra/MatOp/SparseGenMatProd.h>
 /// #include <iostream>
 ///
 /// using namespace Spectra;
 ///
 /// int main()
 /// {
 ///     // A band matrix with 1 on the main diagonal, 2 on the below-main subdiagonal,
 ///     // and 3 on the above-main subdiagonal
 ///     const int n = 10;
 ///     Eigen::SparseMatrix<double> M(n, n);
 ///     M.reserve(Eigen::VectorXi::Constant(n, 3));
 ///     for (int i = 0; i < n; i++)
 ///     {
 ///         M.insert(i, i) = 1.0;
 ///         if (i > 0)
 ///             M.insert(i - 1, i) = 3.0;
 ///         if (i < n - 1)
 ///             M.insert(i + 1, i) = 2.0;
 ///     }
 ///
 ///     // Construct matrix operation object using the wrapper class SparseGenMatProd
 ///     SparseGenMatProd<double> op(M);
 ///
 ///     // Construct eigen solver object, requesting the largest three eigenvalues
 ///     GenEigsSolver<SparseGenMatProd<double>> eigs(op, 3, 6);
 ///
 ///     // Initialize and compute
 ///     eigs.init();
 ///     int nconv = eigs.compute(SortRule::LargestMagn);
 ///
 ///     // Retrieve results
 ///     Eigen::VectorXcd evalues;
 ///     if (eigs.info() == CompInfo::Successful)
 ///         evalues = eigs.eigenvalues();
 ///
 ///     std::cout << "Eigenvalues found:\n" << evalues << std::endl;
 ///
 ///     return 0;
 /// }
 /// \endcode
 template <typename OpType = DenseGenMatProd<double>>
 class GenEigsSolver : public GenEigsBase<OpType, IdentityBOp>
 {
 private:
     using Index = Eigen::Index;
 .
 public:
     ///
     /// Constructor to create a solver object.
     ///
     /// \param op   The matrix operation object that implements
     ///             the matrix-vector multiplication operation of \f$A\f$:
     ///             calculating \f$Av\f$ for any vector \f$v\f$. Users could either
     ///             create the object from the wrapper class such as DenseGenMatProd, or
     ///             define their own that implements all the public members
     ///             as in DenseGenMatProd.
     /// \param nev  Number of eigenvalues requested. This should satisfy \f$1\le nev \le n-2\f$,
     ///             where \f$n\f$ is the size of matrix.
     /// \param ncv  Parameter that controls the convergence speed of the algorithm.
     ///             Typically a larger `ncv` means faster convergence, but it may
     ///             also result in greater memory use and more matrix operations
     ///             in each iteration. This parameter must satisfy \f$nev+2 \le ncv \le n\f$,
     ///             and is advised to take \f$ncv \ge 2\cdot nev + 1\f$.
     ///
     GenEigsSolver(OpType& op, Index nev, Index ncv) :
         GenEigsBase<OpType, IdentityBOp>(op, IdentityBOp(), nev, ncv)
     {}
 };
 .
 }  // namespace Spectra
 .
 #endif  // SPECTRA_GEN_EIGS_SOLVER_H
