ug4
Loading...
Searching...
No Matches
lapack_interface.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2010-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
47#ifndef __H__UG__MARTIN_ALGEBRA__LAPACK_INTERFACE_H__
48#define __H__UG__MARTIN_ALGEBRA__LAPACK_INTERFACE_H__
49
50#ifdef LAPACK_AVAILABLE
51
52#include "lapack.h"
53#include "../small_algebra.h"
54
55#include "lapack_invert.h"
57
58namespace ug{
59
60extern "C" // solve generalized eigenvalue system
61 void dgegv_(char *cCalcLeftEV, char *cCalcRightEV, lapack_int *n,
62 lapack_double *pColMajorMatrixA, lapack_int *lda, lapack_double *pColMajorMatrixB,
63 lapack_int *ldb, lapack_double *alphar, lapack_double *alphai, lapack_double *beta,
64 lapack_double *vl, lapack_int *ldvl, lapack_double *vr, lapack_int *ldvr, lapack_double *pWork,
65 lapack_int *worksize, lapack_int *info);
66
67inline lapack_int gegv(bool calcLeftEV, bool calcRightEV, lapack_int n,
68 lapack_double *pColMajorMatrixA, lapack_int leading_dim_a, lapack_double *pColMajorMatrixB, lapack_int leading_dim_b,
69 lapack_double *alphar, lapack_double *alphai, lapack_double *beta, lapack_double *vl, lapack_int leading_dim_vl,
70 lapack_double *vr, lapack_int leading_dim_vr, lapack_double *work, lapack_int worksize)
71{
72 char jobvl = calcLeftEV ? 'V' : 'N';
73 char jobvr = calcRightEV ? 'V' : 'N';
74 lapack_int info=0;
75 dgegv_(&jobvl, &jobvr, &n, pColMajorMatrixA, &leading_dim_a, pColMajorMatrixB, &leading_dim_b, alphar, alphai, beta, vl, &leading_dim_vl, vr, &leading_dim_vr, work, &worksize, &info);
76 return info;
77}
78
80
81
83// GeneralizedEigenvalueProblem
84//-------------------------------
96template<typename A_type, typename B_type>
97int GeneralizedEigenvalueProblem(DenseMatrix<A_type> &A, DenseMatrix<A_type> &X,
98 DenseVector<B_type> &lambda, DenseMatrix<A_type> &B, bool bSortEigenvalues=false)
99{
100 size_t N = A.num_rows();
101 UG_ASSERT(N == A.num_cols() && N == B.num_rows() && N == B.num_cols(), "");
102
103 B_type alphar; alphar.resize(N);
104 B_type alphai; alphai.resize(N);
105 B_type beta; beta.resize(N);
106
107
108 int worksize = -1;
109 double dWorksize = 0;
110
111 int info;
112
113 if(A_type::ordering == RowMajor)
114 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);
115 if(A_type::ordering == ColMajor)
116 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);
117 UG_COND_THROW(info != 0, "gegv: failed to detect worksize");
118
119 worksize = (int)dWorksize;
120 double *dwork = new double[worksize];
121 if(A_type::ordering == RowMajor)
122 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);
123 if(A_type::ordering == ColMajor)
124 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);
125 UG_COND_THROW(info != 0, "gegv: failed calculate");
126
127 delete[] dwork;
128
129 for(size_t i=0; i<N; i++)
130 {
131 lambda[i] = alphar[i]/beta[i];
132 UG_COND_THROW(dabs(alphai[i]) > 1e-11, "gegv: complex values");
133 }
134
135 if(bSortEigenvalues)
136 {
137
138 // bubblesort
139 for(int i=N-1; i>0; i--)
140 {
141 for(int j=0; j<i; j++)
142 if(lambda[j] > lambda[j+1])
143 {
144 std::swap(lambda[j], lambda[j+1]);
145 for(size_t r = 0; r<X.num_rows(); r++)
146 std::swap(X(r, j), X(r, j+1));
147 }
148 }
149 }
150
151 return info;
152}
153
154
155}
156
157#endif // LAPACK_AVAILABLE
158#endif
#define UG_ASSERT(expr, msg)
Definition assert.h:70
#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
double dabs(double a)
Definition blocks.h:41
@ ColMajor
Definition storage.h:49
@ RowMajor
Definition storage.h:48
int GeneralizedEigenvalueProblem(A_type &A, A_type &X, B_type &lambda, A_type &B, bool bSortEigenvalues=false)
Definition no_lapack.h:41
int lapack_int
Definition solve_deficit.cpp:39
double lapack_double
Definition solve_deficit.cpp:41
class SparseMatrix B_type