Loading...
Searching...
No Matches
dynmatrixev.hh
Go to the documentation of this file.
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
#ifndef DUNE_DYNMATRIXEIGENVALUES_HH
#define DUNE_DYNMATRIXEIGENVALUES_HH
#include <algorithm>
#include <memory>
#include <dune-common-config.hh> // HAVE_LAPACK
#include "dynmatrix.hh"
#include "fmatrixev.hh"
namespace Dune {
namespace DynamicMatrixHelp {
#if HAVE_LAPACK
using Dune::FMatrixHelp::eigenValuesNonsymLapackCall;
#endif
template <typename K, class C>
static void eigenValuesNonSym(const DynamicMatrix<K>& matrix,
DynamicVector<C>& eigenValues,
std::vector<DynamicVector<K>>* eigenVectors = nullptr
)
{
#if HAVE_LAPACK
{
const long int N = matrix.rows();
const char jobvl = 'n';
const char jobvr = eigenVectors ? 'v' : 'n';
// matrix to put into dgeev
auto matrixVector = std::make_unique<double[]>(N*N);
// copy matrix
int row = 0;
for(int i=0; i<N; ++i)
{
for(int j=0; j<N; ++j, ++row)
{
matrixVector[ row ] = matrix[ i ][ j ];
}
}
// working memory
auto eigenR = std::make_unique<double[]>(N);
auto eigenI = std::make_unique<double[]>(N);
const long int lwork = eigenVectors ? 4*N : 3*N;
auto work = std::make_unique<double[]>(lwork);
auto vr = eigenVectors ? std::make_unique<double[]>(N*N) : std::unique_ptr<double[]>{};
// return value information
long int info = 0;
// call LAPACK routine (see fmatrixev_ext.cc)
eigenValuesNonsymLapackCall(&jobvl, &jobvr, &N, matrixVector.get(), &N,
eigenR.get(), eigenI.get(), nullptr, &N, vr.get(), &N, work.get(),
&lwork, &info);
if( info != 0 )
{
std::cerr << "For matrix " << matrix << " eigenvalue calculation failed! " << std::endl;
DUNE_THROW(InvalidStateException,"eigenValues: Eigenvalue calculation failed!");
}
eigenValues.resize(N);
for (int i=0; i<N; ++i)
eigenValues[i] = std::complex<double>(eigenR[i], eigenI[i]);
if (eigenVectors) {
eigenVectors->resize(N);
for (int i = 0; i < N; ++i) {
auto& v = (*eigenVectors)[i];
v.resize(N);
std::copy(vr.get() + N*i, vr.get() + N*(i+1), &v[0]);
}
}
}
#else // #if HAVE_LAPACK
DUNE_THROW(NotImplemented,"LAPACK not found!");
#endif
}
}
}
#endif
dynmatrix.hh
This file implements a dense matrix with dynamic numbers of rows and columns.
fmatrixev.hh
Eigenvalue computations for the FieldMatrix class.
DUNE_THROW
#define DUNE_THROW(E,...)
Definition exceptions.hh:314
Dune
Dune namespace
Definition alignedallocator.hh:13
Dune::DynamicMatrixHelp
Definition dynmatrixev.hh:26
Dune::DynamicMatrixHelp::eigenValuesNonSym
static void eigenValuesNonSym(const DynamicMatrix< K > &matrix, DynamicVector< C > &eigenValues, std::vector< DynamicVector< K > > *eigenVectors=nullptr)
calculates the eigenvalues of a symmetric field matrix
Definition dynmatrixev.hh:41
Dune::DenseMatrix::rows
constexpr size_type rows() const
number of rows
Definition densematrix.hh:714
Dune::DynamicMatrix
Construct a matrix with a dynamic size.
Definition dynmatrix.hh:61
Dune::DynamicVector
Construct a vector with a dynamic size.
Definition dynvector.hh:59
Dune::DynamicVector::resize
void resize(size_type n, value_type c=value_type())
Definition dynvector.hh:141
Dune::NotImplemented
Default exception for dummy implementations.
Definition exceptions.hh:357
Dune::InvalidStateException
Default exception if a function was called while the object is not in a valid state for that function...
Definition exceptions.hh:375