DUNE-FEM (unstable)

dynmatrixev.hh
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
4// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5#ifndef DUNE_DYNMATRIXEIGENVALUES_HH
6#define DUNE_DYNMATRIXEIGENVALUES_HH
7
8#include <algorithm>
9#include <memory>
10
11#include <dune-common-config.hh> // HAVE_LAPACK
12
13#include "dynmatrix.hh"
14#include "fmatrixev.hh"
15
24namespace Dune {
25
26 namespace DynamicMatrixHelp {
27
28#if HAVE_LAPACK
29 using Dune::FMatrixHelp::eigenValuesNonsymLapackCall;
30#endif
31
40 template <typename K, class C>
41 static void eigenValuesNonSym(const DynamicMatrix<K>& matrix,
42 DynamicVector<C>& eigenValues,
43 std::vector<DynamicVector<K>>* eigenVectors = nullptr
44 )
45 {
46
47#if HAVE_LAPACK
48 {
49 const long int N = matrix.rows();
50 const char jobvl = 'n';
51 const char jobvr = eigenVectors ? 'v' : 'n';
52
53
54 // matrix to put into dgeev
55 auto matrixVector = std::make_unique<double[]>(N*N);
56
57 // copy matrix
58 int row = 0;
59 for(int i=0; i<N; ++i)
60 {
61 for(int j=0; j<N; ++j, ++row)
62 {
63 matrixVector[ row ] = matrix[ i ][ j ];
64 }
65 }
66
67 // working memory
68 auto eigenR = std::make_unique<double[]>(N);
69 auto eigenI = std::make_unique<double[]>(N);
70
71 const long int lwork = eigenVectors ? 4*N : 3*N;
72 auto work = std::make_unique<double[]>(lwork);
73 auto vr = eigenVectors ? std::make_unique<double[]>(N*N) : std::unique_ptr<double[]>{};
74
75 // return value information
76 long int info = 0;
77
78 // call LAPACK routine (see fmatrixev_ext.cc)
79 eigenValuesNonsymLapackCall(&jobvl, &jobvr, &N, matrixVector.get(), &N,
80 eigenR.get(), eigenI.get(), nullptr, &N, vr.get(), &N, work.get(),
81 &lwork, &info);
82
83 if( info != 0 )
84 {
85 std::cerr << "For matrix " << matrix << " eigenvalue calculation failed! " << std::endl;
86 DUNE_THROW(InvalidStateException,"eigenValues: Eigenvalue calculation failed!");
87 }
88
89 eigenValues.resize(N);
90 for (int i=0; i<N; ++i)
91 eigenValues[i] = std::complex<double>(eigenR[i], eigenI[i]);
92
93 if (eigenVectors) {
94 eigenVectors->resize(N);
95 for (int i = 0; i < N; ++i) {
96 auto& v = (*eigenVectors)[i];
97 v.resize(N);
98 std::copy(vr.get() + N*i, vr.get() + N*(i+1), &v[0]);
99 }
100 }
101 }
102#else // #if HAVE_LAPACK
103 DUNE_THROW(NotImplemented,"LAPACK not found!");
104#endif
105 }
106 }
107
108}
110#endif
constexpr size_type rows() const
number of rows
Definition: densematrix.hh:709
Construct a matrix with a dynamic size.
Definition: dynmatrix.hh:61
Construct a vector with a dynamic size.
Definition: dynvector.hh:59
Default exception if a function was called while the object is not in a valid state for that function...
Definition: exceptions.hh:281
Default exception for dummy implementations.
Definition: exceptions.hh:263
This file implements a dense matrix with dynamic numbers of rows and columns.
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
Eigenvalue computations for the FieldMatrix class.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
Dune namespace.
Definition: alignedallocator.hh:13
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)