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.h"
19
#include "flxMtx_Eigen_ARP.h"
20
#include "flxMtx_Eigen_GSL.h"
24
void MtxEigenValue(FlxMtx_baseS& Amtx, const tuint M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors, const int Mode)
26
const tuint N = Amtx.nrows();
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());
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() );
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() );
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);
55
throw FlxException_Crude("MtxEigenValue_G_0.3");
58
EV_orientation(M,Eigenvectors);
62
void MtxEigenValue(FlxMtx& Amtx, const tuint M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors, const int Mode)
64
const tuint N = Amtx.nrows();
66
if (N!=Amtx.ncols()) {
67
throw FlxException_Crude("MtxEigenValue_G_2.1");
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());
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() );
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() );
92
throw FlxException_NotImplemented("MtxEigenValue_G_2.3");
93
} else if (Mode == 2) {
94
MtxEigenValue_GSL(Amtx,M,EigenValues,Eigenvectors);
96
throw FlxException_Crude("MtxEigenValue_G_2.3");
99
EV_orientation(M,Eigenvectors);
103
void MtxEigenValue(FlxMtx_baseS& Amtx, FlxMtx_baseS& Bmtx, const tuint M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors, const int Mode)
105
const tuint N = Amtx.nrows();
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());
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() );
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() );
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() );
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);
139
throw FlxException_Crude("MtxEigenValue_G_1.4");
142
EV_orientation(M,Eigenvectors);
145
void EV_orientation(const tuint M, std::vector< flxVec >& Eigenvectors)
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();
157
for (tuint j=0;j<N;++j) {
158
if (fabs(EVp[j])>=r) {
170
throw FlxException_Crude("EV_orientation");
177
void flxeigen_logInfo(std::ostream& lout)
179
lout << " FlxEigen: " << std::endl;
180
lout << " version: " << FLXEIGEN_VERSION_REV << " (" << FLXEIGEN_VERSION_DATE << ")" << std::endl;
181
lout << " compiled with the options ..." << std::endl;
183
lout << " FLX_DEBUG ";
191
lout << " FLX_DEBUG_COUT ";
199
lout << " FLX_USE_ARPACK ";
207
lout << " FLX_USE_GSL ";