// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_TECTONIC_UTILS_MAXEIGENVALUES_HH
#define DUNE_TECTONIC_UTILS_MAXEIGENVALUES_HH

/** \file
 * \brief max eigenvalue computations for the FieldMatrix class
 */

#include <iostream>
#include <cmath>

#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fmatrixev.hh>

template <typename K>
static K maxEigenvalue(const Dune::FieldMatrix<K, 1, 1>& matrix) {
    return matrix[0][0];
}

/** \brief calculates the eigenvalues of a symmetric field matrix
    \param[in]  matrix matrix eigenvalues are calculated for
    \param[out] max eigenvalue
*/
template <typename K>
static K maxEigenvalue(const Dune::FieldMatrix<K, 2, 2>& matrix) {
    const K p = 0.5 * (matrix[0][0] + matrix [1][1]);
    const K p2 = p - matrix[1][1];
    K q = p2 * p2 + matrix[1][0] * matrix[0][1];
    if( q < 0 && q > -1e-14 ) q = 0;
    if (q < 0) {
        std::cout << matrix << std::endl;
        // Complex eigenvalues are either caused by non-symmetric matrices or by round-off errors
        DUNE_THROW(Dune::MathError, "Complex eigenvalue detected (which this implementation cannot handle).");
    }

    // get square root
    q = std::sqrt(q);

    return p + q;
}

/** \brief Calculates the eigenvalues of a symmetric 3x3 field matrix
    \param[in]  matrix eigenvalues are calculated for
    \param[out] max eigenvalue

    \note If the input matrix is not symmetric the behavior of this method is undefined.
 */
template <typename K>
static K maxEigenvalue(const Dune::FieldMatrix<K, 3, 3>& matrix) {
    Dune::FieldVector<K, 3> eigenvalues;
    Dune::FMatrixHelp::eigenValues(matrix, eigenvalues);
    return eigenvalues[0];
}

  /** @} end documentation */

#endif