1
/* Fesslix - Stochastic Analysis
2
* Copyright (C) 2010-2016 Wolfgang Betz
4
* Fesslix is free software: you can redistribute it and/or modify
5
* it under the terms of the GNU General Public License as published by
6
* the Free Software Foundation, either version 3 of the License, or
7
* (at your option) any later version.
9
* Fesslix is distributed in the hope that it will be useful,
10
* but WITHOUT ANY WARRANTY; without even the implied warranty of
11
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12
* GNU General Public License for more details.
14
* You should have received a copy of the GNU General Public License
15
* along with Fesslix. If not, see <http://www.gnu.org/licenses/>.
18
#include "flxMtx_Eigen_GSL.h"
22
#include <gsl/gsl_math.h>
23
#include <gsl/gsl_eigen.h>
24
#include <gsl/gsl_complex.h>
25
#include <gsl/gsl_complex_math.h>
27
void MtxEigenValue_GSL(FlxMtx_baseS& Amtx, const int M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors)
29
const std::size_t n = Amtx.nrows();
30
// get the matrix to solve
31
tdouble* data = new tdouble[n*n];
32
Amtx.get_VecPointer_full(data);
33
gsl_matrix_view m = gsl_matrix_view_array (data, n, n);
34
// allocate some memory
35
gsl_vector *eval = gsl_vector_alloc (n);
36
gsl_matrix *evec = gsl_matrix_alloc (n, n);
37
gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc (n);
39
gsl_eigen_symmv (&m.matrix, eval, evec, w);
40
gsl_eigen_symmv_free (w);
41
gsl_eigen_symmv_sort (eval, evec, GSL_EIGEN_SORT_VAL_DESC);
44
for (int i = 0; i < M; ++i)
46
EigenValues[i] = gsl_vector_get (eval, i);
47
gsl_vector_view evec_i = gsl_matrix_column (evec, i);
48
for (tuint j=0;j<n;++j) {
49
Eigenvectors[i][j] = gsl_vector_get(&evec_i.vector,j);
54
gsl_vector_free (eval);
55
gsl_matrix_free (evec);
60
const bool ensure_complex_is_real_GSL(gsl_complex cn)
62
const tdouble ca = gsl_complex_abs(cn);
63
const tdouble ia = fabs(GSL_IMAG(cn));
64
return (ia/ca <= GlobalVar.TOL());
68
void MtxEigenValue_GSL(FlxMtx& Amtx, const int M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors)
70
const std::size_t n = Amtx.nrows();
71
// get the matrix to solve
72
tdouble* data = &(Amtx.operator()(0,0));
73
gsl_matrix_view m = gsl_matrix_view_array (data, n, n);
74
// allocate some memory
75
gsl_vector_complex *eval = gsl_vector_complex_alloc (n);
76
gsl_matrix_complex *evec = gsl_matrix_complex_alloc (n, n);
77
gsl_eigen_nonsymmv_workspace * w = gsl_eigen_nonsymmv_alloc (n);
79
gsl_eigen_nonsymmv (&m.matrix, eval, evec, w);
80
gsl_eigen_nonsymmv_free (w);
81
gsl_eigen_nonsymmv_sort (eval, evec, GSL_EIGEN_SORT_ABS_DESC);
84
for (int i = 0; i < M; ++i)
86
gsl_complex eval_i = gsl_vector_complex_get (eval, i);
87
gsl_vector_complex_view evec_i = gsl_matrix_complex_column (evec, i);
89
// make sure the eigenvalues are real
90
if (!ensure_complex_is_real_GSL(eval_i)) {
91
throw FlxException_Crude("MtxEigenValue_GSL_1");
93
EigenValues[i] = GSL_REAL(eval_i);
95
// make sure the eigenvectors are real
96
for (tuint j=0;j<n;++j) {
97
gsl_complex z = gsl_vector_complex_get(&evec_i.vector, j);
98
if (!ensure_complex_is_real_GSL(z)) {
99
throw FlxException_Crude("MtxEigenValue_GSL_2");
101
Eigenvectors[i][j] = GSL_REAL(z);
106
gsl_vector_complex_free(eval);
107
gsl_matrix_complex_free(evec);
111
void MtxEigenValue_GSL(FlxMtx_baseS& Amtx, FlxMtx_baseS& Bmtx, const int M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors, const tuint sortID)
113
const size_t n = Amtx.nrows();
114
// get the matrix to solve
115
tdouble* data = new tdouble[n*n];
116
tdouble* dataB = new tdouble[n*n];
117
Amtx.get_VecPointer_full(data);
118
Bmtx.get_VecPointer_full(dataB);
119
gsl_matrix_view m = gsl_matrix_view_array (data, n, n);
120
gsl_matrix_view b = gsl_matrix_view_array (dataB, n, n);
121
// allocate some memory
122
gsl_vector *eval = gsl_vector_alloc (n);
123
gsl_matrix *evec = gsl_matrix_alloc (n, n);
124
gsl_eigen_gensymmv_workspace * w = gsl_eigen_gensymmv_alloc(n);
126
gsl_eigen_gensymmv (&m.matrix,&b.matrix, eval, evec, w);
127
gsl_eigen_gensymmv_free (w);
129
gsl_eigen_gensymmv_sort (eval, evec, GSL_EIGEN_SORT_VAL_DESC);
131
gsl_eigen_gensymmv_sort (eval, evec, GSL_EIGEN_SORT_VAL_ASC);
134
// assign the results
135
for (int i = 0; i < M; ++i)
137
EigenValues[i] = gsl_vector_get (eval, i);
138
gsl_vector_view evec_i = gsl_matrix_column (evec, i);
139
for (tuint j=0;j<n;++j) {
140
Eigenvectors[i][j] = gsl_vector_get(&evec_i.vector,j);
145
gsl_vector_free (eval);
146
gsl_matrix_free (evec);
151
void MtxEigenValue_GSLstab(FlxMtx_baseS& Amtx, const int M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors)
153
throw FlxException_NotImplemented("MtxEigenValue_GSLstab");
156
void MtxEigenValue_GSLstab(FlxMtx_baseS& Amtx, FlxMtx_baseS& Bmtx, const int M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors)
158
GlobalVar.slog(4) << std::endl;
159
GlobalVar.slog(4) << "Solving eigenvalue problem - stabilized version" << std::endl;
160
GlobalVar.slog(4) << "-----------------------------------------------" << std::endl;
161
GlobalVar.slog(4) << " Problem: Bx = lMx" << std::endl << std::endl;
162
const size_t n = Amtx.nrows();
163
// get the matrix to solve
164
tdouble* data = new tdouble[n*n];
165
tdouble* dataB = new tdouble[n*n];
166
Amtx.get_VecPointer_full(data);
167
Bmtx.get_VecPointer_full(dataB);
168
gsl_matrix_view m = gsl_matrix_view_array (data, n, n);
169
gsl_matrix_view b = gsl_matrix_view_array (dataB, n, n);
170
// allocate some memory
171
gsl_vector *eval = gsl_vector_alloc (n);
172
gsl_matrix *evec = gsl_matrix_alloc (n, n);
173
// stabilize matrix B
174
const tdouble CONDMIN = 1e-7;
175
gsl_eigen_symmv_workspace * wB = gsl_eigen_symmv_alloc(n);
176
gsl_eigen_symmv(&b.matrix,eval, evec, wB);
177
gsl_eigen_symmv_free(wB);
178
gsl_eigen_symmv_sort (eval,evec,GSL_EIGEN_SORT_VAL_DESC);
179
tdouble maxev = gsl_vector_get (eval, 0);
180
tdouble minev = maxev;
183
for (tuint i = 1; i < n; ++i) {
184
const tdouble cev = gsl_vector_get (eval, i);
186
gsl_vector_set(eval,i,CONDMIN*maxev);
191
if (cev/maxev < CONDMIN) {
192
gsl_vector_set(eval,i,CONDMIN*maxev);
197
GlobalVar.slog(4) << " Matrix M" << std::endl;
198
GlobalVar.slog(4) << " --------" << std::endl;
199
GlobalVar.slog(4) << " condition number: " << GlobalVar.Double2String(maxev/minev) << std::endl;
200
GlobalVar.slog(4) << " corrected eigenvalues: " << corn << std::endl;
201
GlobalVar.slog(4) << " largest positive ev.: " << maxev << std::endl;
202
GlobalVar.slog(4) << " smallest positive ev.: " << minev << std::endl;
204
GlobalVar.slog(4) << " smallest negative ev.: " << negev << std::endl;
207
GlobalVar.slog(4) << " new condition number: " << GlobalVar.Double2String(1.0/CONDMIN) << std::endl;
209
for (tuint k=0;k<n;++k) {
210
for (tuint j=0;j<n;++j) {
212
for (tuint i=0;i<n;++i) {
213
t += gsl_matrix_get(evec,k,n-i-1)*gsl_matrix_get(evec,j,n-i-1)*gsl_vector_get (eval,n-i-1);
215
gsl_matrix_set(&b.matrix,k,j,t);
220
gsl_eigen_gensymmv_workspace * w = gsl_eigen_gensymmv_alloc(n);
221
gsl_eigen_gensymmv (&m.matrix,&b.matrix, eval, evec, w);
222
gsl_eigen_gensymmv_free (w);
223
gsl_eigen_gensymmv_sort (eval, evec, GSL_EIGEN_SORT_VAL_DESC);
225
// assign the results
226
for (int i = 0; i < M; ++i)
228
EigenValues[i] = gsl_vector_get (eval, i);
229
gsl_vector_view evec_i = gsl_matrix_column (evec, i);
230
for (tuint j=0;j<n;++j) {
231
Eigenvectors[i][j] = gsl_vector_get(&evec_i.vector,j);
236
gsl_vector_free (eval);
237
gsl_matrix_free (evec);
240
GlobalVar.slog(4) << std::endl;
245
void MtxEigenValue_GSL(FlxMtx& Amtx, const int M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors)
247
std::ostringstream ssV;
248
ssV << "Fesslix has been compiled without GSL-support.";
249
throw FlxException("MtxEigenValue_GSL_76", ssV.str() );
252
void MtxEigenValue_GSL(FlxMtx_baseS& Amtx, const int M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors)
254
std::ostringstream ssV;
255
ssV << "Fesslix has been compiled without GSL-support.";
256
throw FlxException("MtxEigenValue_GSL_77", ssV.str() );
259
void MtxEigenValue_GSL(FlxMtx_baseS& Amtx, FlxMtx_baseS& Bmtx, const int M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors, const tuint sortID)
261
std::ostringstream ssV;
262
ssV << "Fesslix has been compiled without GSL-support.";
263
throw FlxException("MtxEigenValue_GSL_78", ssV.str() );
266
void MtxEigenValue_GSLstab(FlxMtx_baseS& Amtx, const int M, flxVec& EigenValues, std::vector<flxVec>& Eigenvectors)
268
std::ostringstream ssV;
269
ssV << "Fesslix has been compiled without GSL-support.";
270
throw FlxException("MtxEigenValue_GSLstab_77", ssV.str() );
273
void MtxEigenValue_GSLstab(FlxMtx_baseS& Amtx, FlxMtx_baseS& Bmtx, const int M, flxVec& EigenValues, std::vector<flxVec>& Eigenvectors)
275
std::ostringstream ssV;
276
ssV << "Fesslix has been compiled without GSL-support.";
277
throw FlxException("MtxEigenValue_GSLstab_78", ssV.str() );