MAST
lapack_dgeev_interface.h
Go to the documentation of this file.
1 /*
2  * MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
3  * Copyright (C) 2013-2019 Manav Bhatia
4  *
5  * This library is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU Lesser General Public
7  * License as published by the Free Software Foundation; either
8  * version 2.1 of the License, or (at your option) any later version.
9  *
10  * This library is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13  * Lesser General Public License for more details.
14  *
15  * You should have received a copy of the GNU Lesser General Public
16  * License along with this library; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18  */
19 
20 #ifndef __mast__lapack_dgeev_interface_h__
21 #define __mast__lapack_dgeev_interface_h__
22 
23 
24 // MAST includes
25 #include "base/mast_data_types.h"
26 
27 
28 extern "C" {
29 
30  /*
31  * =====================================================================
32  * Purpose
33  * =======
34  *
35  * DGEEV computes for an N-by-N real nonsymmetric matrix A, the
36  * eigenvalues and, optionally, the left and/or right eigenvectors.
37  *
38  * The right eigenvector v(j) of A satisfies
39  * A * v(j) = lambda(j) * v(j)
40  * where lambda(j) is its eigenvalue.
41  * The left eigenvector u(j) of A satisfies
42  * u(j)**H * A = lambda(j) * u(j)**H
43  * where u(j)**H denotes the conjugate transpose of u(j).
44  *
45  * The computed eigenvectors are normalized to have Euclidean norm
46  * equal to 1 and largest component real.
47  *
48  * Arguments
49  * =========
50  *
51  * JOBVL (input) CHARACTER*1
52  * = 'N': left eigenvectors of A are not computed;
53  * = 'V': left eigenvectors of A are computed.
54  *
55  * JOBVR (input) CHARACTER*1
56  * = 'N': right eigenvectors of A are not computed;
57  * = 'V': right eigenvectors of A are computed.
58  *
59  * N (input) INTEGER
60  * The order of the matrix A. N >= 0.
61  *
62  * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
63  * On entry, the N-by-N matrix A.
64  * On exit, A has been overwritten.
65  *
66  * LDA (input) INTEGER
67  * The leading dimension of the array A. LDA >= max(1,N).
68  *
69  * WR (output) DOUBLE PRECISION array, dimension (N)
70  * WI (output) DOUBLE PRECISION array, dimension (N)
71  * WR and WI contain the real and imaginary parts,
72  * respectively, of the computed eigenvalues. Complex
73  * conjugate pairs of eigenvalues appear consecutively
74  * with the eigenvalue having the positive imaginary part
75  * first.
76  *
77  * VL (output) DOUBLE PRECISION array, dimension (LDVL,N)
78  * If JOBVL = 'V', the left eigenvectors u(j) are stored one
79  * after another in the columns of VL, in the same order
80  * as their eigenvalues.
81  * If JOBVL = 'N', VL is not referenced.
82  * If the j-th eigenvalue is real, then u(j) = VL(:,j),
83  * the j-th column of VL.
84  * If the j-th and (j+1)-st eigenvalues form a complex
85  * conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
86  * u(j+1) = VL(:,j) - i*VL(:,j+1).
87  *
88  * LDVL (input) INTEGER
89  * The leading dimension of the array VL. LDVL >= 1; if
90  * JOBVL = 'V', LDVL >= N.
91  *
92  * VR (output) DOUBLE PRECISION array, dimension (LDVR,N)
93  * If JOBVR = 'V', the right eigenvectors v(j) are stored one
94  * after another in the columns of VR, in the same order
95  * as their eigenvalues.
96  * If JOBVR = 'N', VR is not referenced.
97  * If the j-th eigenvalue is real, then v(j) = VR(:,j),
98  * the j-th column of VR.
99  * If the j-th and (j+1)-st eigenvalues form a complex
100  * conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
101  * v(j+1) = VR(:,j) - i*VR(:,j+1).
102  *
103  * LDVR (input) INTEGER
104  * The leading dimension of the array VR. LDVR >= 1; if
105  * JOBVR = 'V', LDVR >= N.
106  *
107  * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
108  * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
109  *
110  * LWORK (input) INTEGER
111  * The dimension of the array WORK. LWORK >= max(1,3*N), and
112  * if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N. For good
113  * performance, LWORK must generally be larger.
114  *
115  * If LWORK = -1, then a workspace query is assumed; the routine
116  * only calculates the optimal size of the WORK array, returns
117  * this value as the first entry of the WORK array, and no error
118  * message related to LWORK is issued by XERBLA.
119  *
120  * INFO (output) INTEGER
121  * = 0: successful exit
122  * < 0: if INFO = -i, the i-th argument had an illegal value.
123  * > 0: if INFO = i, the QR algorithm failed to compute all the
124  * eigenvalues, and no eigenvectors have been computed;
125  * elements i+1:N of WR and WI contain eigenvalues which
126  * have converged.
127  * =====================================================================
128  */
129  extern void dgeev_(const char* jobvl,
130  const char* jobvr,
131  int* n,
132  double* a,
133  int* lda,
134  double* w_r,
135  double* w_i,
136  double* vl,
137  int* ldvl,
138  double* vr,
139  int* ldvr,
140  double* work,
141  int* lwork,
142  int* info);
143 
144 }
145 
146 
147 namespace MAST {
148 
149 
151 
152  public:
153 
155  info_val(-1)
156  { }
157 
162  void compute(const RealMatrixX& A,
163  bool computeEigenvectors = true);
164 
165  ComputationInfo info() const;
166 
167  const RealMatrixX& A() const {
168  libmesh_assert(info_val == 0);
169  return this->_A;
170  }
171 
172 
173  const ComplexVectorX& eig_vals() const {
174  libmesh_assert(info_val == 0);
175  return this->W;
176  }
177 
179  libmesh_assert(info_val == 0);
180  return this->VL;
181  }
182 
184  libmesh_assert(info_val == 0);
185  return this->VR;
186  }
187 
194  libmesh_assert(info_val == 0);
195 
196  // this product should be an identity matrix
197  ComplexMatrixX r = this->VL.conjugate().transpose() * this->VR;
198 
199  // scale the right eigenvectors by the inverse of the inner-product
200  // diagonal
201  Complex val;
202  for (unsigned int i=0; i<_A.cols(); i++) {
203  val = r(i,i);
204  if (std::abs(val) > 0.)
205  this->VR.col(i) *= (1./val);
206  }
207  }
208 
209  void print_inner_product(std::ostream& out) const {
210  libmesh_assert(info_val == 0);
211  ComplexMatrixX r;
212  r = this->VL.conjugate().transpose() * _A * this->VR;
213  out << "conj(VL)' * A * VR" << std::endl
214  << r << std::endl;
215 
216  r = this->VL.conjugate().transpose() * this->VR;
217  out << "conj(VL)' * B * VR" << std::endl
218  << r << std::endl;
219 
220  }
221 
222  protected:
223 
225 
227 
229 
231 
232  int info_val;
233  };
234 
235 }
236 
237 #endif // __mast__lapack_dgeev_interface_h__
const ComplexVectorX & eig_vals() const
ComputationInfo info() const
void scale_eigenvectors_to_identity_innerproduct()
Scales the right eigenvector so that the inner product with respect to the B matrix is equal to an Id...
const ComplexMatrixX & right_eigenvectors() const
Matrix< Complex, Dynamic, 1 > ComplexVectorX
void dgeev_(const char *jobvl, const char *jobvr, int *n, double *a, int *lda, double *w_r, double *w_i, double *vl, int *ldvl, double *vr, int *ldvr, double *work, int *lwork, int *info)
libMesh::Complex Complex
Matrix< Real, Dynamic, Dynamic > RealMatrixX
Matrix< Complex, Dynamic, Dynamic > ComplexMatrixX
void compute(const RealMatrixX &A, bool computeEigenvectors=true)
computes the eigensolution for .
const ComplexMatrixX & left_eigenvectors() const
const RealMatrixX & A() const
void print_inner_product(std::ostream &out) const