~wbetz/fesslix/flxeigen

« back to all changes in this revision

Viewing changes to src/flxMtx_Eigen_GSL.cpp

  • Committer: Wolfgang Betz
  • Date: 2016-06-17 14:28:23 UTC
  • Revision ID: wolfgang.betz@fesslix.org-20160617142823-ao5x99dj066qkckg
Initial import

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
/* Fesslix - Stochastic Analysis
 
2
 * Copyright (C) 2010-2016 Wolfgang Betz
 
3
 *
 
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.
 
8
 *
 
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.
 
13
 *
 
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/>. 
 
16
 */
 
17
 
 
18
#include "flxMtx_Eigen_GSL.h"
 
19
 
 
20
#if FLX_USE_GSL
 
21
 
 
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>
 
26
 
 
27
void MtxEigenValue_GSL(FlxMtx_baseS& Amtx, const int M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors)
 
28
{
 
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);
 
38
  // solve the problem
 
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);
 
42
  
 
43
  // assign the results
 
44
    for (int i = 0; i < M; ++i)
 
45
      {
 
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);
 
50
        }
 
51
      }
 
52
 
 
53
  // free the memory
 
54
    gsl_vector_free (eval);
 
55
    gsl_matrix_free (evec);
 
56
    delete [] data;
 
57
 
 
58
}
 
59
 
 
60
const bool ensure_complex_is_real_GSL(gsl_complex cn) 
 
61
{  
 
62
  const tdouble ca = gsl_complex_abs(cn);
 
63
  const tdouble ia = fabs(GSL_IMAG(cn));
 
64
  return (ia/ca <= GlobalVar.TOL());    
 
65
}
 
66
 
 
67
 
 
68
void MtxEigenValue_GSL(FlxMtx& Amtx, const int M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors)
 
69
{
 
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);
 
78
  // solve the problem
 
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);
 
82
  
 
83
  // assign the results
 
84
    for (int i = 0; i < M; ++i)
 
85
      {
 
86
        gsl_complex eval_i = gsl_vector_complex_get (eval, i);
 
87
        gsl_vector_complex_view evec_i = gsl_matrix_complex_column (evec, i);
 
88
        
 
89
        // make sure the eigenvalues are real
 
90
        if (!ensure_complex_is_real_GSL(eval_i)) {
 
91
          throw FlxException_Crude("MtxEigenValue_GSL_1");
 
92
        }
 
93
        EigenValues[i] = GSL_REAL(eval_i);
 
94
        
 
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");
 
100
            }
 
101
            Eigenvectors[i][j] = GSL_REAL(z);
 
102
          }
 
103
      }
 
104
 
 
105
  // free the memory
 
106
    gsl_vector_complex_free(eval);
 
107
    gsl_matrix_complex_free(evec);
 
108
 
 
109
}
 
110
 
 
111
void MtxEigenValue_GSL(FlxMtx_baseS& Amtx, FlxMtx_baseS& Bmtx, const int M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors, const tuint sortID)
 
112
{
 
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);
 
125
  // solve the problem
 
126
    gsl_eigen_gensymmv (&m.matrix,&b.matrix, eval, evec, w);
 
127
    gsl_eigen_gensymmv_free (w);
 
128
    if (sortID==1) {
 
129
      gsl_eigen_gensymmv_sort (eval, evec, GSL_EIGEN_SORT_VAL_DESC);
 
130
    } else {
 
131
      gsl_eigen_gensymmv_sort (eval, evec, GSL_EIGEN_SORT_VAL_ASC);
 
132
    }
 
133
  
 
134
  // assign the results
 
135
    for (int i = 0; i < M; ++i)
 
136
      {
 
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);
 
141
        }
 
142
      }
 
143
 
 
144
  // free the memory
 
145
    gsl_vector_free (eval);
 
146
    gsl_matrix_free (evec);
 
147
    delete [] data;
 
148
    delete [] dataB;
 
149
}
 
150
 
 
151
void MtxEigenValue_GSLstab(FlxMtx_baseS& Amtx, const int M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors)
 
152
{
 
153
  throw FlxException_NotImplemented("MtxEigenValue_GSLstab");
 
154
}
 
155
 
 
156
void MtxEigenValue_GSLstab(FlxMtx_baseS& Amtx, FlxMtx_baseS& Bmtx, const int M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors)
 
157
{
 
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;
 
181
    tdouble negev = 0.0;
 
182
    tuint corn = 0;
 
183
    for (tuint i = 1; i < n; ++i) {
 
184
      const tdouble cev = gsl_vector_get (eval, i);
 
185
      if (cev<0.0) {
 
186
        gsl_vector_set(eval,i,CONDMIN*maxev);
 
187
        ++corn;
 
188
        negev = cev;
 
189
      } else {
 
190
        minev = cev;
 
191
        if (cev/maxev < CONDMIN) {
 
192
          gsl_vector_set(eval,i,CONDMIN*maxev);
 
193
          ++corn;
 
194
        }
 
195
      }
 
196
    }
 
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;
 
203
    if (negev<0.0) {
 
204
    GlobalVar.slog(4) << "  smallest negative ev.: " << negev << std::endl;
 
205
    }
 
206
    if (corn>0) {
 
207
    GlobalVar.slog(4) << "  new condition number:  " << GlobalVar.Double2String(1.0/CONDMIN) << std::endl;
 
208
    }
 
209
    for (tuint k=0;k<n;++k) {
 
210
      for (tuint j=0;j<n;++j) {
 
211
        tdouble t = 0.0;
 
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);
 
214
        }
 
215
        gsl_matrix_set(&b.matrix,k,j,t);
 
216
      }
 
217
    }
 
218
 
 
219
  // solve the problem
 
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);
 
224
  
 
225
  // assign the results
 
226
    for (int i = 0; i < M; ++i)
 
227
      {
 
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);
 
232
        }
 
233
      }
 
234
 
 
235
  // free the memory
 
236
    gsl_vector_free (eval);
 
237
    gsl_matrix_free (evec);
 
238
    delete [] data;
 
239
    delete [] dataB;
 
240
  GlobalVar.slog(4) << std::endl;
 
241
}
 
242
 
 
243
#else
 
244
 
 
245
void MtxEigenValue_GSL(FlxMtx& Amtx, const int M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors)
 
246
{
 
247
  std::ostringstream ssV;
 
248
  ssV << "Fesslix has been compiled without GSL-support.";
 
249
  throw FlxException("MtxEigenValue_GSL_76", ssV.str() );
 
250
}
 
251
 
 
252
void MtxEigenValue_GSL(FlxMtx_baseS& Amtx, const int M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors)
 
253
{
 
254
  std::ostringstream ssV;
 
255
  ssV << "Fesslix has been compiled without GSL-support.";
 
256
  throw FlxException("MtxEigenValue_GSL_77", ssV.str() );
 
257
}
 
258
 
 
259
void MtxEigenValue_GSL(FlxMtx_baseS& Amtx, FlxMtx_baseS& Bmtx, const int M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors, const tuint sortID)
 
260
{
 
261
  std::ostringstream ssV;
 
262
  ssV << "Fesslix has been compiled without GSL-support.";
 
263
  throw FlxException("MtxEigenValue_GSL_78", ssV.str() );
 
264
}
 
265
 
 
266
void MtxEigenValue_GSLstab(FlxMtx_baseS& Amtx, const int M, flxVec& EigenValues, std::vector<flxVec>& Eigenvectors)
 
267
{
 
268
  std::ostringstream ssV;
 
269
  ssV << "Fesslix has been compiled without GSL-support.";
 
270
  throw FlxException("MtxEigenValue_GSLstab_77", ssV.str() );
 
271
}
 
272
 
 
273
void MtxEigenValue_GSLstab(FlxMtx_baseS& Amtx, FlxMtx_baseS& Bmtx, const int M, flxVec& EigenValues, std::vector<flxVec>& Eigenvectors)
 
274
{
 
275
  std::ostringstream ssV;
 
276
  ssV << "Fesslix has been compiled without GSL-support.";
 
277
  throw FlxException("MtxEigenValue_GSLstab_78", ssV.str() );
 
278
}
 
279
 
 
280
#endif
 
281