* LAPACKSupport::unusable after
* this operation.
*
+ * The optional arguments allow
+ * to compute left and right
+ * eigenvectors as well.
+ *
* @note Calls the LAPACK
* function Xgeev.
*/
- void compute_eigenvalues ();
+ void compute_eigenvalues (const bool right_eigenvectors = false,
+ const bool left_eigenvectors = false);
/**
* Retrieve eigenvalue after
/**
* Print the matrix and allow
* formatting of entries.
- *
+ *
* The parameters allow for a
* flexible setting of the output
* format:
* some LAPACK functions.
*/
mutable std::vector<number> work;
+
/**
* Real parts of
* eigenvalues. Filled by
*/
std::vector<number> wi;
+ /**
+ * Space where left eigenvectors
+ * can be stored.
+ */
+ std::vector<number> vl;
+
+ /**
+ * Space where right eigenvectors
+ * can be stored.
+ */
+ std::vector<number> vr;
};
/**@}*/
// $Id$
// Version: $Name$
//
-// Copyright (C) 2005 by the deal.II authors
+// Copyright (C) 2005, 2006 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
template <typename number>
void
-LAPACKFullMatrix<number>::compute_eigenvalues()
+LAPACKFullMatrix<number>::compute_eigenvalues(
+ const bool right,
+ const bool left)
{
const int nn = this->n_cols();
wr.resize(nn);
wi.resize(nn);
+ if (right) vr.resize(nn*nn);
+ if (left) vl.resize(nn*nn);
+
number* values = const_cast<number*> (this->data());
int info = 0;
int lwork = 1;
+ const char * const jobvr = (right) ? (&V) : (&N);
+ const char * const jobvl = (left) ? (&V) : (&N);
+
// Optimal workspace query:
// The LAPACK routine DGEEV requires
lwork = -1;
work.resize(1);
- geev(&N, &N, &nn, values, &nn,
+ geev(jobvl, jobvr, &nn, values, &nn,
&wr[0], &wi[0],
- 0, &one, 0, &one,
+ &vl[0], &nn, &vr[0], &nn,
&work[0], &lwork, &info);
// geev returns info=0 on
// success. Since we only queried
work.resize((unsigned int) lwork);
// Finally compute the eigenvalues.
- geev(&N, &N, &nn, values, &nn,
+ geev(jobvl, jobvr, &nn, values, &nn,
&wr[0], &wi[0],
- 0, &one, 0, &one,
+ &vl[0], &nn, &vr[0], &nn,
&work[0], &lwork, &info);
// Negative return value implies a
// wrong argument. This should be