~wbetz/fesslix/flxeigen

« back to all changes in this revision

Viewing changes to src/flxMtx_Eigen.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.h"
 
19
#include "flxMtx_Eigen_ARP.h"
 
20
#include "flxMtx_Eigen_GSL.h"
 
21
 
 
22
 
 
23
 
 
24
void MtxEigenValue(FlxMtx_baseS& Amtx, const tuint M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors, const int Mode)
 
25
{
 
26
  const tuint N = Amtx.nrows();
 
27
  if (N < M) {
 
28
    std::ostringstream ssV;
 
29
    ssV << "Cannot compute more Eigenvalues (" << M << ") than number of DOFs in the system(" <<  N << ").";
 
30
    throw FlxException("MtxEigenValue_G_0.0", ssV.str());
 
31
  }
 
32
  
 
33
  #if FLX_DEBUG
 
34
    if (EigenValues.get_N() != M) {
 
35
      std::ostringstream ssV;
 
36
      ssV << "Problem with size of vector. (" << EigenValues.get_N() << "; " << M << ")";
 
37
      throw FlxException("MtxEigenValue_G_0.1", ssV.str() );
 
38
    }
 
39
    for (tuint i=0;i<M;i++) {
 
40
      if (Eigenvectors[i].get_N()!=N) {
 
41
        std::ostringstream ssV;
 
42
        ssV << "Problem with size of vector.";
 
43
        throw FlxException("MtxEigenValue_G_0.2", ssV.str() );
 
44
      }
 
45
    }
 
46
  #endif
 
47
  
 
48
  if (Mode == 1) {
 
49
    MtxEigenValue_ARP(Amtx,M,EigenValues,Eigenvectors);
 
50
  } else if (Mode == 2) {
 
51
    MtxEigenValue_GSL(Amtx,M,EigenValues,Eigenvectors);
 
52
  } else if (Mode == 3) {
 
53
    MtxEigenValue_GSLstab(Amtx,M,EigenValues,Eigenvectors);
 
54
  } else {
 
55
    throw FlxException_Crude("MtxEigenValue_G_0.3");
 
56
  }
 
57
  
 
58
  EV_orientation(M,Eigenvectors);
 
59
  
 
60
}
 
61
 
 
62
void MtxEigenValue(FlxMtx& Amtx, const tuint M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors, const int Mode)
 
63
{
 
64
  const tuint N = Amtx.nrows();
 
65
  #if FLX_DEBUG
 
66
    if (N!=Amtx.ncols()) {
 
67
      throw FlxException_Crude("MtxEigenValue_G_2.1");
 
68
    }
 
69
  #endif
 
70
  if (N < M) {
 
71
    std::ostringstream ssV;
 
72
    ssV << "Cannot compute more Eigenvalues (" << M << ") than number of DOFs in the system(" <<  N << ").";
 
73
    throw FlxException("MtxEigenValue_G_2.2", ssV.str());
 
74
  }
 
75
  
 
76
  #if FLX_DEBUG
 
77
    if (EigenValues.get_N() != M) {
 
78
      std::ostringstream ssV;
 
79
      ssV << "Problem with size of vector. (" << EigenValues.get_N() << "; " << M << ")";
 
80
      throw FlxException("MtxEigenValue_G_2.1", ssV.str() );
 
81
    }
 
82
    for (tuint i=0;i<M;i++) {
 
83
      if (Eigenvectors[i].get_N()!=N) {
 
84
        std::ostringstream ssV;
 
85
        ssV << "Problem with size of vector.";
 
86
        throw FlxException("MtxEigenValue_G_2.2", ssV.str() );
 
87
      }
 
88
    }
 
89
  #endif
 
90
  
 
91
  if (Mode == 1) {
 
92
    throw FlxException_NotImplemented("MtxEigenValue_G_2.3");
 
93
  } else if (Mode == 2) {
 
94
    MtxEigenValue_GSL(Amtx,M,EigenValues,Eigenvectors);
 
95
  } else {
 
96
    throw FlxException_Crude("MtxEigenValue_G_2.3");
 
97
  }
 
98
  
 
99
  EV_orientation(M,Eigenvectors);
 
100
  
 
101
}
 
102
 
 
103
void MtxEigenValue(FlxMtx_baseS& Amtx, FlxMtx_baseS& Bmtx, const tuint M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors, const int Mode)
 
104
{
 
105
  const tuint N = Amtx.nrows();
 
106
  if (N < M) {
 
107
    std::ostringstream ssV;
 
108
    ssV << "Cannot compute more Eigenvalues (" << M << ") than number of DOFs in the system (" <<  N << ").";
 
109
    throw FlxException("MtxEigenValue_G_1.0", ssV.str());
 
110
  }
 
111
  
 
112
  #if FLX_DEBUG
 
113
    if (EigenValues.get_N() != M) {
 
114
      std::ostringstream ssV;
 
115
      ssV << "Problem with size of vector. (" << EigenValues.get_N() << "; " << M << ")";
 
116
      throw FlxException("MtxEigenValue_G_1.1", ssV.str() );
 
117
    }
 
118
    for (tuint i=0;i<M;i++) {
 
119
      if (Eigenvectors[i].get_N()!=Amtx.ncols()) {
 
120
        std::ostringstream ssV;
 
121
        ssV << "Problem with size of vector.";
 
122
        throw FlxException("MtxEigenValue_G_1.2", ssV.str() );
 
123
      }
 
124
    }
 
125
    if (Amtx.ncols()!=Bmtx.ncols()) {
 
126
      std::ostringstream ssV;
 
127
      ssV << "Matrix sizes do not match";
 
128
      throw FlxException("MtxEigenValue_G_1.3", ssV.str() );
 
129
    }
 
130
  #endif
 
131
  
 
132
  if (Mode == 1) {
 
133
    MtxEigenValue_ARP(Amtx,Bmtx,M,EigenValues,Eigenvectors);
 
134
  } else if (Mode == 2) {
 
135
    MtxEigenValue_GSL(Amtx,Bmtx,M,EigenValues,Eigenvectors,1);
 
136
  } else if (Mode == 3) {
 
137
    MtxEigenValue_GSLstab(Amtx,Bmtx,M,EigenValues,Eigenvectors);
 
138
  } else {
 
139
    throw FlxException_Crude("MtxEigenValue_G_1.4");
 
140
  }
 
141
  
 
142
  EV_orientation(M,Eigenvectors);
 
143
}
 
144
 
 
145
void EV_orientation(const tuint M, std::vector< flxVec >& Eigenvectors)
 
146
{
 
147
  // define orientation of eigenvectors
 
148
  for (tuint i=0;i<M;++i) {
 
149
    flxVec& EV = Eigenvectors[i];
 
150
    const tdouble mN = EV.get_absMean();
 
151
    const tdouble r = mN*0.1;
 
152
    const tuint N = EV.get_N();
 
153
    const tdouble* EVp = EV.get_tmp_vptr_const();
 
154
    #if FLX_DEBUG
 
155
    bool b=false;
 
156
    #endif
 
157
    for (tuint j=0;j<N;++j) {
 
158
      if (fabs(EVp[j])>=r) {
 
159
        if (EVp[j]<0.) {
 
160
          EV *= -1;
 
161
        }
 
162
        #if FLX_DEBUG
 
163
        b=true;
 
164
        #endif
 
165
        break;
 
166
      }
 
167
    }
 
168
    #if FLX_DEBUG
 
169
    if (b==false) {
 
170
      throw FlxException_Crude("EV_orientation");
 
171
    }
 
172
    #endif
 
173
  }
 
174
}
 
175
 
 
176
 
 
177
void flxeigen_logInfo(std::ostream& lout)
 
178
{
 
179
  lout << " FlxEigen: " << std::endl;
 
180
  lout << "   version: " << FLXEIGEN_VERSION_REV << " (" << FLXEIGEN_VERSION_DATE << ")" << std::endl;
 
181
  lout << "   compiled with the options ..." << std::endl;
 
182
  // FLXDEBUG
 
183
    lout << "     FLX_DEBUG                     ";
 
184
    #if FLX_DEBUG
 
185
      lout << "ON";
 
186
    #else 
 
187
      lout << "OFF";
 
188
    #endif
 
189
    lout << std::endl;
 
190
  // FLXDEBUG_COUT
 
191
    lout << "     FLX_DEBUG_COUT                ";
 
192
    #if FLX_DEBUG_COUT
 
193
      lout << "ON";
 
194
    #else 
 
195
      lout << "OFF";
 
196
    #endif
 
197
    lout << std::endl;
 
198
  // FLX_USE_ARPACK
 
199
    lout << "     FLX_USE_ARPACK                ";
 
200
    #if FLX_USE_ARPACK
 
201
      lout << "ON";
 
202
    #else 
 
203
      lout << "OFF";
 
204
    #endif  
 
205
    lout << std::endl;
 
206
  // FLX_USE_GSL
 
207
    lout << "     FLX_USE_GSL                   ";
 
208
    #if FLX_USE_GSL 
 
209
      lout << "ON";
 
210
    #else 
 
211
      lout << "OFF";
 
212
    #endif  
 
213
    lout << std::endl;
 
214
}
 
215
 
 
216
 
 
217