MAST
lapack_zggevx_interface.cpp
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 
21 // MAST includes
23 
24 
25 void
27  const ComplexMatrixX &B,
28  bool computeEigenvectors) {
29 
30  libmesh_assert(A.cols() == A.rows() &&
31  B.cols() == A.rows() &&
32  B.cols() == B.rows());
33 
34  _A = A;
35  _B = B;
36 
38  Amat = _A,
39  Bmat = _B;
40 
41  int n = (int)A.cols();
42 
43  char BAL='B', L='N',R='N', S='E';
44 
45  if (computeEigenvectors) {
46 
47  L = 'V'; R = 'V'; S='B';
48  VL.setZero(n, n);
49  VR.setZero(n, n);
50  }
51 
52  int
53  lwork =4*(n*n+n),
54  l_rwork =8*n,
55  ilo = 0,
56  ihi = 0;
57  info_val =-1,
58 
59  alpha.setZero(n);
60  beta.setZero(n);
62  work = ComplexVectorX::Zero(lwork);
64  rwork = RealVectorX::Zero(l_rwork),
65  lscale = RealVectorX::Zero(n),
66  rscale = RealVectorX::Zero(n),
67  rconde = RealVectorX::Zero(n),
68  rcondv = RealVectorX::Zero(n);
69 
70  Complex
71  *a_vals = Amat.data(),
72  *b_vals = Bmat.data(),
73  *alpha_v = alpha.data(),
74  *beta_v = beta.data(),
75  *VL_v = VL.data(),
76  *VR_v = VR.data(),
77  *work_v = work.data();
78 
79  Real
80  *rwork_v = rwork.data(),
81  *lscale_v = lscale.data(),
82  *rscale_v = rscale.data(),
83  *rconde_v = rconde.data(),
84  *rcondv_v = rcondv.data(),
85  abnrm = 0.,
86  bbnrm = 0.;
87 
88  std::vector<int>
89  iwork(n+2, 0);
90 
91  bool
92  bwork[n];
93 
94  zggevx_(&BAL, &L, &R, &S, &n,
95  &(a_vals[0]), &n,
96  &(b_vals[0]), &n,
97  &(alpha_v[0]), &(beta_v[0]),
98  &(VL_v[0]), &n,
99  &(VR_v[0]), &n,
100  &ilo, &ihi,
101  &(lscale_v[0]), &(rscale_v[0]),
102  &abnrm, &bbnrm,
103  &(rconde_v[0]), &(rcondv_v[0]),
104  &(work_v[0]), &lwork,
105  &(rwork_v[0]),
106  &(iwork[0]),
107  &(bwork[0]),
108  &info_val);
109 
110  if (info_val != 0)
111  libMesh::out
112  << "Warning!! ZGGEVX returned with nonzero info = "
113  << info_val << std::endl;
114 }
115 
116 
int zggevx_(char *balanc, char *jobvl, char *jobvr, char *sens, int *n, std::complex< double > *a, int *lda, std::complex< double > *b, int *ldb, std::complex< double > *alpha, std::complex< double > *beta, std::complex< double > *vl, int *ldvl, std::complex< double > *vr, int *ldvr, int *ilo, int *ihi, double *lscale, double *rscale, double *abnrm, double *bbnrm, double *rconde, double *rcondv, std::complex< double > *work, int *lwork, double *rwork, int *iwork, bool *bwork, int *info)
virtual void compute(const ComplexMatrixX &A, const ComplexMatrixX &B, bool computeEigenvectors=true)
computes the eigensolution for .
Matrix< Complex, Dynamic, 1 > ComplexVectorX
const ComplexMatrixX & B() const
libMesh::Real Real
libMesh::Complex Complex
Matrix< Complex, Dynamic, Dynamic > ComplexMatrixX
Matrix< Real, Dynamic, 1 > RealVectorX
const ComplexMatrixX & A() const