MAST
lapack_zggev_base.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_zggev_interface_base_h__
21 #define __mast__lapack_zggev_interface_base_h__
22 
23 
24 // MAST includes
25 #include "base/mast_data_types.h"
26 
27 
28 
29 
30 namespace MAST {
31 
33 
34  public:
35 
37  info_val(-1)
38  { }
39 
44  virtual void compute(const ComplexMatrixX& A,
45  const ComplexMatrixX& B,
46  bool computeEigenvectors = true) = 0;
47 
48  const ComplexMatrixX& A() const {
49  libmesh_assert(info_val == 0);
50  return this->_A;
51  }
52 
53 
54  const ComplexMatrixX& B() const {
55  libmesh_assert(info_val == 0);
56  return this->_B;
57  }
58 
59 
60  const ComplexVectorX& alphas() const {
61  libmesh_assert(info_val == 0);
62  return this->alpha;
63  }
64 
65  const ComplexVectorX& betas() const {
66  libmesh_assert(info_val == 0);
67  return this->beta;
68  }
69 
71  libmesh_assert(info_val == 0);
72  return this->VL;
73  }
74 
76  libmesh_assert(info_val == 0);
77  return this->VR;
78  }
79 
86  libmesh_assert(info_val == 0);
87 
88  // scale the right eigenvectors so that they all have unit norm
89  Real l2 = 0.;
90 
91  for (unsigned int i=0; i<this->VR.cols(); i++) {
92 
93  l2 = this->VR.col(i).norm();
94  libmesh_assert(l2 > 0.);
95 
96  // scale all right eigenvectors to unit length
97  this->VR.col(i) /= l2;
98  }
99 
100  // this product should be an identity matrix
101  ComplexMatrixX r = this->VL.conjugate().transpose() * _B * this->VR;
102 
103  // scale the right eigenvectors by the inverse of the inner-product
104  // diagonal
105  Complex val;
106  for (unsigned int i=0; i<_B.cols(); i++) {
107  val = r(i,i);
108  if (std::abs(val) > 0.)
109  this->VL.col(i) *= (1./val);
110  }
111  }
112 
113  void print_inner_product(std::ostream& out) const {
114  libmesh_assert(info_val == 0);
115  ComplexMatrixX r;
116  r = this->VL.conjugate().transpose() * _A * this->VR;
117  out << "conj(VL)' * A * VR" << std::endl
118  << r << std::endl;
119 
120  r = this->VL.conjugate().transpose() * _B * this->VR;
121  out << "conj(VL)' * B * VR" << std::endl
122  << r << std::endl;
123 
124  }
125 
126  protected:
127 
129 
131 
133 
135 
137 
139 
140  int info_val;
141  };
142 }
143 
144 
145 
146 #endif // __mast__lapack_zggev_interface_base_h__
const ComplexMatrixX & left_eigenvectors() const
Matrix< Complex, Dynamic, 1 > ComplexVectorX
virtual void compute(const ComplexMatrixX &A, const ComplexMatrixX &B, bool computeEigenvectors=true)=0
computes the eigensolution for .
void print_inner_product(std::ostream &out) const
const ComplexMatrixX & B() const
libMesh::Real Real
libMesh::Complex Complex
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, Dynamic > ComplexMatrixX
const ComplexVectorX & betas() const
const ComplexMatrixX & A() const
const ComplexVectorX & alphas() const