~wbetz/fesslix/flxeigen

« back to all changes in this revision

Viewing changes to src/flxMtx_Eigen_ARP.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_ARP.h"
 
19
 
 
20
#if FLX_USE_ARPACK
 
21
 
 
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>
 
27
 
 
28
 
 
29
template<class ARMATRIX, class ARFLOAT>
 
30
void Solution(ARMATRIX &A, ARluSymStdEig<ARFLOAT> &Prob)
 
31
/*
 
32
  Prints eigenvalues and eigenvectors of symmetric eigen-problems
 
33
  on standard "cout" stream.
 
34
*/
 
35
{
 
36
 
 
37
  int     i, n, nconv, mode;
 
38
  ARFLOAT *Ax;
 
39
  ARFLOAT *ResNorm;
 
40
 
 
41
  n     = Prob.GetN();
 
42
  nconv = Prob.ConvergedEigenvalues();
 
43
  mode  = Prob.GetMode();
 
44
 
 
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;
 
47
  switch (mode) {
 
48
  case 1:
 
49
    *(GlobalVar.get_cout()) << "Regular mode" << std::endl;
 
50
    break;
 
51
  case 3:
 
52
    *(GlobalVar.get_cout()) << "Shift and invert mode" << std::endl;
 
53
  }
 
54
  *(GlobalVar.get_cout()) << std::endl;
 
55
 
 
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;
 
62
 
 
63
  if (Prob.EigenvaluesFound()) {
 
64
 
 
65
    // Printing eigenvalues.
 
66
 
 
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;
 
70
    }
 
71
    *(GlobalVar.get_cout()) << std::endl;
 
72
  }
 
73
 
 
74
  if (Prob.EigenvectorsFound()) {
 
75
 
 
76
    // Printing the residual norm || A*x - lambda*x ||
 
77
    // for the nconv accurately computed eigenvectors.
 
78
 
 
79
    Ax      = new ARFLOAT[n];
 
80
    ResNorm = new ARFLOAT[nconv+1];
 
81
 
 
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));
 
86
    }
 
87
 
 
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";
 
91
    }
 
92
    *(GlobalVar.get_cout()) << "\n";
 
93
 
 
94
    delete[] Ax;
 
95
    delete[] ResNorm;
 
96
 
 
97
  }
 
98
 
 
99
} // Solution
 
100
 
 
101
 
 
102
void MtxEigenValue_ARP(FlxMtx_baseS& Amtx, const int M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors)
 
103
{
 
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.
 
107
  const char uplo='U';
 
108
  
 
109
  ARdsSymMatrix<tdouble> matrix(n, A,'U');
 
110
  ARluSymStdEig<tdouble> dprob(M, matrix);
 
111
 
 
112
  // Finding eigenvalues and eigenvectors.
 
113
  dprob.FindEigenvectors();
 
114
  
 
115
  if (M != dprob.ConvergedEigenvalues() ) {
 
116
    std::ostringstream ssV;
 
117
    ssV << "Not all eigenvalues converged.";
 
118
    throw FlxException("MtxEigenValue_ARP_2", ssV.str() );
 
119
  }
 
120
  if (!dprob.EigenvaluesFound()) {
 
121
    std::ostringstream ssV;
 
122
    ssV << "Eigenvalues could not be obtained.";
 
123
    throw FlxException("MtxEigenValue_ARP_3", ssV.str() );
 
124
  }
 
125
  if (!dprob.EigenvectorsFound()) {
 
126
    std::ostringstream ssV;
 
127
    ssV << "Eigenvectors could not be obtained.";
 
128
    throw FlxException("MtxEigenValue_ARP_4", ssV.str() );
 
129
  }
 
130
 
 
131
  for (int i=0; i<M; i++) {
 
132
    EigenValues[M-1-i] = dprob.Eigenvalue(i);
 
133
  }
 
134
  
 
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];
 
139
    }
 
140
  }
 
141
  
 
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;
 
148
  
 
149
  // Printing solution.
 
150
//   Solution(matrix, dprob);
 
151
 
 
152
}
 
153
 
 
154
void MtxEigenValue_ARP(FlxMtx_baseS& Amtx, FlxMtx_baseS& Bmtx, const int M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors)
 
155
{
 
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();
 
160
  const char uplo='U';
 
161
  
 
162
  ARdsSymMatrix<tdouble> matrixA(n, A,uplo);
 
163
  ARdsSymMatrix<tdouble> matrixB(n, B,uplo);
 
164
  ARluSymGenEig<tdouble> dprob(M, matrixA, matrixB);
 
165
 
 
166
  // Finding eigenvalues and eigenvectors.
 
167
  dprob.FindEigenvectors();
 
168
  
 
169
  if (M != dprob.ConvergedEigenvalues() ) {
 
170
    std::ostringstream ssV;
 
171
    ssV << "Not all eigenvalues converged.";
 
172
    throw FlxException("MtxEigenValue_ARP_302", ssV.str() );
 
173
  }
 
174
  if (!dprob.EigenvaluesFound()) {
 
175
    std::ostringstream ssV;
 
176
    ssV << "Eigenvalues could not be obtained.";
 
177
    throw FlxException("MtxEigenValue_ARP_303", ssV.str() );
 
178
  }
 
179
  if (!dprob.EigenvectorsFound()) {
 
180
    std::ostringstream ssV;
 
181
    ssV << "Eigenvectors could not be obtained.";
 
182
    throw FlxException("MtxEigenValue_ARP_304", ssV.str() );
 
183
  }
 
184
 
 
185
  for (int i=0; i<M; i++) {
 
186
    EigenValues[M-1-i] = dprob.Eigenvalue(i);
 
187
  }
 
188
  
 
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];
 
193
    }
 
194
  }
 
195
}
 
196
 
 
197
#else
 
198
 
 
199
void MtxEigenValue_ARP(FlxMtx_baseS& Amtx, const int M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors)
 
200
{
 
201
  std::ostringstream ssV;
 
202
  ssV << "Fesslix has been compiled without ARPACK-support.";
 
203
  throw FlxException("MtxEigenValue_ARP_77", ssV.str() );
 
204
}
 
205
 
 
206
void MtxEigenValue_ARP(FlxMtx_baseS& Amtx, FlxMtx_baseS& Bmtx, const int M, flxVec& EigenValues, std::vector< flxVec >& Eigenvectors)
 
207
{
 
208
  std::ostringstream ssV;
 
209
  ssV << "Fesslix has been compiled without ARPACK-support.";
 
210
  throw FlxException("MtxEigenValue_ARP_78", ssV.str() );
 
211
}
 
212
 
 
213
#endif
 
214