ug4
Loading...
Searching...
No Matches
eigenvalue2.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2013-2014: G-CSC, Goethe University Frankfurt
3 * Author: Martin Rupp
4 *
5 * This file is part of UG4.
6 *
7 * UG4 is free software: you can redistribute it and/or modify it under the
8 * terms of the GNU Lesser General Public License version 3 (as published by the
9 * Free Software Foundation) with the following additional attribution
10 * requirements (according to LGPL/GPL v3 §7):
11 *
12 * (1) The following notice must be displayed in the Appropriate Legal Notices
13 * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 *
15 * (2) The following notice must be displayed at a prominent place in the
16 * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 *
18 * (3) The following bibliography is recommended for citation and must be
19 * preserved in all covered files:
20 * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 * parallel geometric multigrid solver on hierarchically distributed grids.
22 * Computing and visualization in science 16, 4 (2013), 151-164"
23 * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 * flexible software system for simulating pde based models on high performance
25 * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 *
27 * This program is distributed in the hope that it will be useful,
28 * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 * GNU Lesser General Public License for more details.
31 */
32
33
34#ifndef __H__UG__EIGENVALUE__
35#define __H__UG__EIGENVALUE__
36
37#ifdef LAPACK_AVAILABLE
38
39#include "lapack.h"
40#include "../small_algebra.h"
41
42#include "lapack_invert.h"
44#include "lapack_interface.h"
45#include <complex>
46#include "../small_matrix/densematrix.h"
47#include "common/error.h"
48
49namespace ug{
50
52
53
55// GeneralizedEigenvalueProblem
56//-------------------------------
68template<typename A_type, typename TLambdaVectorType>
69int GeneralizedEigenvalueProblemComplex(DenseMatrix<A_type> &A, DenseMatrix<A_type> &X,
70 TLambdaVectorType &lambda, DenseMatrix<A_type> &B, bool bSortEigenvalues=false)
71{
72 size_t N = A.num_rows();
73 THROW_IF_NOT_EQUAL_4(N, A.num_cols(), B.num_rows(), B.num_cols());
74
75 std::vector<double> alphar; alphar.resize(N);
76 std::vector<double> alphai; alphai.resize(N);
77 std::vector<double> beta; beta.resize(N);
78
79
80 int worksize = -1;
81 double dWorksize = 0;
82
83 int info;
84
85 if(A_type::ordering == RowMajor)
86 info = gegv(true, false, N, &A(0,0), N, &B(0,0), N, &alphar[0], &alphai[0], &beta[0], &X(0,0), N, NULL, 0, &dWorksize, worksize);
87 if(A_type::ordering == ColMajor)
88 info = gegv(false, true, N, &A(0,0), N, &B(0,0), N, &alphar[0], &alphai[0], &beta[0], NULL, N, &X(0,0), N, &dWorksize, worksize);
89 UG_COND_THROW(info != 0, "gegv: failed to detect worksize");
90
91 worksize = (int)dWorksize;
92 double *dwork = new double[worksize];
93 if(A_type::ordering == RowMajor)
94 info = gegv(true, false, N, &A(0,0), N, &B(0,0), N, &alphar[0], &alphai[0], &beta[0], &X(0,0), N, NULL, 0, dwork, worksize);
95 if(A_type::ordering == ColMajor)
96 info = gegv(false, true, N, &A(0,0), N, &B(0,0), N, &alphar[0], &alphai[0], &beta[0], NULL, N, &X(0,0), N, dwork, worksize);
97 UG_COND_THROW(info != 0, "gegv: failed calculate");
98
99 delete[] dwork;
100
101 for(size_t i=0; i<N; i++)
102 lambda[i] = std::complex<double>(alphar[i]/beta[i], alphai[i]/beta[i]);
103
104 if(bSortEigenvalues)
105 {
106
107 // bubblesort
108 for(int i=N-1; i>0; i--)
109 {
110 for(int j=0; j<i; j++)
111 if(std::abs(lambda[j]) > std::abs(lambda[j+1]))
112 {
113 std::swap(lambda[j], lambda[j+1]);
114 for(size_t r = 0; r<X.num_rows(); r++)
115 std::swap(X(r, j), X(r, j+1));
116 }
117 }
118 }
119
120 return info;
121}
122
123
124}
125#else // LAPACK_AVAILABLE
126namespace ug{
127template<typename A_type, typename TLambdaVectorType>
129 TLambdaVectorType &lambda, DenseMatrix<A_type> &B, bool bSortEigenvalues=false)
130{
131 UG_ASSERT(0, "GeneralizedEigenvalueProblemComplex is only implemented for LAPACK at the moment");
132 return 0;
133}
134}
135#endif // LAPACK_AVAILABLE
136#endif // __H__UG__EIGENVALUE__
Definition densematrix.h:57
#define UG_ASSERT(expr, msg)
Definition assert.h:70
#define THROW_IF_NOT_EQUAL_4(s1, s2, s3, s4)
Definition error.h:183
#define UG_COND_THROW(cond, msg)
UG_COND_THROW(cond, msg) : performs a UG_THROW(msg) if cond == true.
Definition error.h:61
the ug namespace
@ ColMajor
Definition storage.h:49
@ RowMajor
Definition storage.h:48
int GeneralizedEigenvalueProblemComplex(DenseMatrix< A_type > &A, DenseMatrix< A_type > &X, TLambdaVectorType &lambda, DenseMatrix< A_type > &B, bool bSortEigenvalues=false)
Definition eigenvalue2.h:128