*/
void compute_eigenvalues (const bool right_eigenvectors = false,
const bool left_eigenvectors = false);
+
+ /**
+ * Compute eigenvalues and
+ * eigenvectors of a real symmetric
+ * matrix. Only eigenvalues in the
+ * interval (lower_bound, upper_bound]
+ * are computed with the absolute
+ * tolerance abs_accuracy. An approximate
+ * eigenvalue is accepted as converged
+ * when it is determined to lie in an
+ * interval [a,b] of width less than or
+ * equal to abs_accuracy + eps * max( |a|,|b| ),
+ * where eps is the machine precision.
+ * If abs_accuracy is less than
+ * or equal to zero, then eps*|t| will
+ * be used in its place, where |t| is the
+ * 1-norm of the tridiagonal matrix obtained
+ * by reducing A to tridiagonal form.
+ * Eigenvalues will be computed most accurately
+ * when abs_accuracy is set to twice the
+ * underflow threshold, not zero.
+ * After this routine has
+ * been called, all eigenvalues in
+ * (lower_bound, upper_bound] will be
+ * stored in eigenvalues and the
+ * corresponding eigenvectors will be stored
+ * in the columns of eigenvectors, whose
+ * dimension is set accordingly.
+ *
+ * @note Calls the LAPACK function
+ * Xsyevx. For this to work, ./configure
+ * has to be told to use LAPACK.
+ */
+ void compute_eigenvalues_symmetric(
+ const number lower_bound,
+ const number upper_bound,
+ const number abs_accuracy,
+ Vector<number> & eigenvalues,
+ FullMatrix<number> & eigenvectors);
/**
* Compute generalized eigenvalues
void dsyev_ (const char *jobz, const char *uplo, const int *n,
double *A, const int *lda, double *w,
double *work, const int *lwork, int *info);
+// Same functionality as dsyev_ but with more options: E.g.
+// Compute only eigenvalues in a specific interval,
+// Compute only eigenvalues with a specific index,
+// Set tolerance for eigenvalue computation
+void dsyevx_ (const char* jobz, const char* range,
+ const char* uplo, const int* n, double* A, const int* lda,
+ const double* vl, const double* vu,
+ const int* il, const int* iu, const double* abstol,
+ int* m, double* w, double* z,
+ const int* ldz, double* work, const int* lwork, int* iwork,
+ int* ifail, int* info);
// Generalized eigenvalues and eigenvectors of
// 1: A*x = lambda*B*x; 2: A*B*x = lambda*x; 3: B*A*x = lambda*x
// A and B are symmetric and B is definite
}
+template <typename number>
+void
+LAPACKFullMatrix<number>::compute_eigenvalues_symmetric(
+ const number lower_bound,
+ const number upper_bound,
+ const number abs_accuracy,
+ Vector<number> & eigenvalues,
+ FullMatrix<number> & eigenvectors)
+{
+ Assert(state == matrix, ExcState(state));
+ const int nn = (this->n_cols() > 0 ? this->n_cols() : 1);
+ Assert(static_cast<unsigned int>(nn) == this->n_rows(), ExcNotQuadratic());
+
+ wr.resize(nn);
+ LAPACKFullMatrix<number> matrix_eigenvectors(nn, nn);
+
+ number* values_A = const_cast<number*> (this->data());
+ number* values_eigenvectors = const_cast<number*> (matrix_eigenvectors.data());
+
+ int info(0),
+ lwork(1),
+ n_eigenpairs(0);
+ const char * const jobz(&V);
+ const char * const uplo(&U);
+ const char * const range(&V);
+ const int * const dummy(&one);
+ std::vector<int> iwork(static_cast<unsigned int> (5*nn));
+ std::vector<int> ifail(static_cast<unsigned int> (nn));
+
+
+ // Optimal workspace query:
+
+ // The LAPACK routine ?SYEVX requires
+ // a sufficient large workspace variable,
+ // minimum requirement is
+ // work.size>=3*nn-1.
+ // However, to improve performance, a
+ // somewhat larger workspace may be needed.
+
+ // SOME implementations of the LAPACK routine
+ // provide a workspace query call,
+ // info:=0, lwork:=-1
+ // which returns an optimal value for the
+ // size of the workspace array
+ // (the PETSc 2.3.0 implementation does NOT
+ // provide this functionality!).
+
+ // define the DEAL_II_LIBLAPACK_NOQUERYMODE flag to
+ // disable the workspace query.
+#ifndef DEAL_II_LIBLAPACK_NOQUERYMODE
+ lwork = -1;
+ work.resize(1);
+
+ syevx (jobz, range,
+ uplo, &nn, values_A, &nn,
+ &lower_bound, &upper_bound,
+ dummy, dummy, &abs_accuracy,
+ &n_eigenpairs, &wr[0], values_eigenvectors,
+ &nn, &work[0], &lwork, &iwork[0],
+ &ifail[0], &info);
+ // syevx returns info=0 on
+ // success. Since we only queried
+ // the optimal size for work,
+ // everything else would not be
+ // acceptable.
+ Assert (info == 0, ExcInternalError());
+ // Allocate working array according
+ // to suggestion.
+ lwork = (int) (work[0]+.1);
+#else
+ lwork = 8*nn > 1 ? 8*nn : 1; // no query mode
+#endif
+ // resize workspace arrays
+ work.resize(static_cast<unsigned int> (lwork));
+
+ // Finally compute the eigenvalues.
+ syevx (jobz, range,
+ uplo, &nn, values_A, &nn,
+ &lower_bound, &upper_bound,
+ dummy, dummy, &abs_accuracy,
+ &n_eigenpairs, &wr[0], values_eigenvectors,
+ &nn, &work[0], &lwork, &iwork[0],
+ &ifail[0], &info);
+
+ // Negative return value implies a
+ // wrong argument. This should be
+ // internal.
+ Assert (info >=0, ExcInternalError());
+ if (info != 0)
+ std::cerr << "LAPACK error in syevx" << std::endl;
+
+ eigenvalues.reinit(n_eigenpairs);
+ eigenvectors.reinit(nn, n_eigenpairs, true);
+ for(unsigned int i=0; i < static_cast<unsigned int> (n_eigenpairs); ++i)
+ {
+ eigenvalues(i) = wr[i];
+ unsigned int col_begin(i*nn);
+ for (unsigned int j=0; j < static_cast<unsigned int> (nn); ++j)
+ {
+ eigenvectors(j,i) = values_eigenvectors[col_begin+j];
+ }
+ }
+
+ state = LAPACKSupport::State(unusable);
+}
+
+
template <typename number>
void
LAPACKFullMatrix<number>::compute_generalized_eigenvalues_symmetric(
// internal.
Assert (info >=0, ExcInternalError());
if (info != 0)
- std::cerr << "LAPACK error in sygv" << std::endl;
+ std::cerr << "LAPACK error in sygvx" << std::endl;
eigenvalues.reinit(n_eigenpairs);
eigenvectors.resize(n_eigenpairs);