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_ARP.h"
22
#include <arpack++/ardsmat.h>
23
#include <arpack++/ardssym.h>
24
#include <arpack++/ardgsym.h>
25
#include <arpack++/blas1c.h>
26
#include <arpack++/lapackc.h>
29
template<class ARMATRIX, class ARFLOAT>
30
void Solution(ARMATRIX &A, ARluSymStdEig<ARFLOAT> &Prob)
32
Prints eigenvalues and eigenvectors of symmetric eigen-problems
33
on standard "cout" stream.
37
int i, n, nconv, mode;
42
nconv = Prob.ConvergedEigenvalues();
43
mode = Prob.GetMode();
45
*(GlobalVar.get_cout()) << std::endl << std::endl << "Testing ARPACK++ class ARluSymStdEig \n";
46
*(GlobalVar.get_cout()) << "Real symmetric eigenvalue problem: A*x - lambda*x" << std::endl;
49
*(GlobalVar.get_cout()) << "Regular mode" << std::endl;
52
*(GlobalVar.get_cout()) << "Shift and invert mode" << std::endl;
54
*(GlobalVar.get_cout()) << std::endl;
56
*(GlobalVar.get_cout()) << "Dimension of the system : " << n << std::endl;
57
*(GlobalVar.get_cout()) << "Number of 'requested' eigenvalues : " << Prob.GetNev() << std::endl;
58
*(GlobalVar.get_cout()) << "Number of 'converged' eigenvalues : " << nconv << std::endl;
59
*(GlobalVar.get_cout()) << "Number of Arnoldi vectors generated: " << Prob.GetNcv() << std::endl;
60
*(GlobalVar.get_cout()) << "Number of iterations taken : " << Prob.GetIter() << std::endl;
61
*(GlobalVar.get_cout()) << std::endl;
63
if (Prob.EigenvaluesFound()) {
65
// Printing eigenvalues.
67
*(GlobalVar.get_cout()) << "Eigenvalues:" << std::endl;
68
for (i=0; i<nconv; i++) {
69
*(GlobalVar.get_cout()) << " lambda[" << (i+1) << "]: " << Prob.Eigenvalue(i) << std::endl;
71
*(GlobalVar.get_cout()) << std::endl;
74
if (Prob.EigenvectorsFound()) {
76
// Printing the residual norm || A*x - lambda*x ||
77
// for the nconv accurately computed eigenvectors.
80
ResNorm = new ARFLOAT[nconv+1];
82
for (i=0; i<nconv; i++) {
83
A.MultMv(Prob.RawEigenvector(i), Ax);
84
axpy(n, -Prob.Eigenvalue(i), Prob.RawEigenvector(i), 1, Ax, 1);
85
ResNorm[i] = nrm2(n, Ax, 1)/fabs(Prob.Eigenvalue(i));
88
for (i=0; i<nconv; i++) {
89
*(GlobalVar.get_cout()) << "||A*x(" << (i+1) << ") - lambda(" << (i+1);
90
*(GlobalVar.get_cout()) << ")*x(" << (i+1) << ")||: " << ResNorm[i] << "\n";
92
*(GlobalVar.get_cout()) << "\n";
102
void MtxEigenValue_ARP(FlxMtx_baseS& Amtx, const int M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors)
104
// Defining variables;
105
const int n = Amtx.ncols(); // Dimension of the problem.
106
tdouble* A= Amtx.get_VecPointer(); // Pointer to an array that stores the lower triangular elements of A.
109
ARdsSymMatrix<tdouble> matrix(n, A,'U');
110
ARluSymStdEig<tdouble> dprob(M, matrix);
112
// Finding eigenvalues and eigenvectors.
113
dprob.FindEigenvectors();
115
if (M != dprob.ConvergedEigenvalues() ) {
116
std::ostringstream ssV;
117
ssV << "Not all eigenvalues converged.";
118
throw FlxException("MtxEigenValue_ARP_2", ssV.str() );
120
if (!dprob.EigenvaluesFound()) {
121
std::ostringstream ssV;
122
ssV << "Eigenvalues could not be obtained.";
123
throw FlxException("MtxEigenValue_ARP_3", ssV.str() );
125
if (!dprob.EigenvectorsFound()) {
126
std::ostringstream ssV;
127
ssV << "Eigenvectors could not be obtained.";
128
throw FlxException("MtxEigenValue_ARP_4", ssV.str() );
131
for (int i=0; i<M; i++) {
132
EigenValues[M-1-i] = dprob.Eigenvalue(i);
135
for (int i=0; i<M; i++) {
136
tdouble* dp = dprob.RawEigenvector(i);
137
for (int j=0;j<n;j++) {
138
Eigenvectors[M-1-i][j]=dp[j];
142
// *(GlobalVar.get_cout()) << "Dimension of the system : " << n << std::endl;
143
// *(GlobalVar.get_cout())<< "Number of 'requested' eigenvalues : " << dprob.GetNev() << std::endl;
144
// *(GlobalVar.get_cout()) << "Number of 'converged' eigenvalues : " << dprob.ConvergedEigenvalues() << std::endl;
145
// *(GlobalVar.get_cout()) << "Number of Arnoldi vectors generated: " << dprob.GetNcv() << std::endl;
146
// *(GlobalVar.get_cout()) << "Number of iterations taken : " << dprob.GetIter() << std::endl;
147
// *(GlobalVar.get_cout()) << std::endl;
149
// Printing solution.
150
// Solution(matrix, dprob);
154
void MtxEigenValue_ARP(FlxMtx_baseS& Amtx, FlxMtx_baseS& Bmtx, const int M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors)
156
// Defining variables;
157
const int n = Amtx.ncols(); // Dimension of the problem.
158
tdouble* A= Amtx.get_VecPointer(); // Pointer to an array that stores the lower triangular elements of A.
159
tdouble* B=Bmtx.get_VecPointer();
162
ARdsSymMatrix<tdouble> matrixA(n, A,uplo);
163
ARdsSymMatrix<tdouble> matrixB(n, B,uplo);
164
ARluSymGenEig<tdouble> dprob(M, matrixA, matrixB);
166
// Finding eigenvalues and eigenvectors.
167
dprob.FindEigenvectors();
169
if (M != dprob.ConvergedEigenvalues() ) {
170
std::ostringstream ssV;
171
ssV << "Not all eigenvalues converged.";
172
throw FlxException("MtxEigenValue_ARP_302", ssV.str() );
174
if (!dprob.EigenvaluesFound()) {
175
std::ostringstream ssV;
176
ssV << "Eigenvalues could not be obtained.";
177
throw FlxException("MtxEigenValue_ARP_303", ssV.str() );
179
if (!dprob.EigenvectorsFound()) {
180
std::ostringstream ssV;
181
ssV << "Eigenvectors could not be obtained.";
182
throw FlxException("MtxEigenValue_ARP_304", ssV.str() );
185
for (int i=0; i<M; i++) {
186
EigenValues[M-1-i] = dprob.Eigenvalue(i);
189
for (int i=0; i<M; i++) {
190
tdouble* dp = dprob.RawEigenvector(i);
191
for (int j=0;j<n;j++) {
192
Eigenvectors[M-1-i][j]=dp[j];
199
void MtxEigenValue_ARP(FlxMtx_baseS& Amtx, const int M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors)
201
std::ostringstream ssV;
202
ssV << "Fesslix has been compiled without ARPACK-support.";
203
throw FlxException("MtxEigenValue_ARP_77", ssV.str() );
206
void MtxEigenValue_ARP(FlxMtx_baseS& Amtx, FlxMtx_baseS& Bmtx, const int M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors)
208
std::ostringstream ssV;
209
ssV << "Fesslix has been compiled without ARPACK-support.";
210
throw FlxException("MtxEigenValue_ARP_78", ssV.str() );