#! /bin/sh
-# From configure.in Revision: 24994 .
+# From configure.in Revision: 24995 .
# Guess values for system-dependent variables and create Makefiles.
# Generated by GNU Autoconf 2.68 for deal.II 7.2.pre.
#
fi
done
-for ac_func in dgetrs_ sgetrs_ dstev_ sstev_ dsygv_ ssygv_
+for ac_func in dgetrs_ sgetrs_ dstev_ sstev_ dsygv_ ssygv_ dsygvx_ ssygvx_
do :
as_ac_var=`$as_echo "ac_cv_func_$ac_func" | $as_tr_sh`
ac_fn_cxx_check_func "$LINENO" "$ac_func" "$as_ac_var"
AC_CHECK_FUNCS([daxpy_ saxpy_ dgemv_ sgemv_ dgeev_ sgeev_ dgeevx_ sgeevx_])
AC_CHECK_FUNCS([dgemm_ sgemm_ dgetrf_ sgetrf_ dgetri_ sgetri_])
AC_CHECK_FUNCS([dgeqrf_ sgeqrf_ dormqr_ sormqr_ dorgqr_ sorgqr_ dtrtrs_ strtrs_])
-AC_CHECK_FUNCS([dgetrs_ sgetrs_ dstev_ sstev_ dsygv_ ssygv_])
+AC_CHECK_FUNCS([dgetrs_ sgetrs_ dstev_ sstev_ dsygv_ ssygv_ dsygvx_ ssygvx_])
dnl Singular value decomposition
AC_CHECK_FUNCS([dgesvd_ sgesvd_ dgesdd_ sgesdd_ dgelsd_ sgelsd_])
dnl Check Bessel functions in GNU libc
/* Define to 1 if you have the `dstev_' function. */
#undef HAVE_DSTEV_
+/* Define to 1 if you have the `dsygvx_' function. */
+#undef HAVE_DSYGVX_
+
/* Define to 1 if you have the `dsygv_' function. */
#undef HAVE_DSYGV_
/**
* Compute generalized eigenvalues
- * and (optionally) eigenvectors of
+ * and eigenvectors of
* a real generalized symmetric
* eigenproblem of the form
* itype = 1: $Ax=\lambda B x$
* where A is this matrix.
* A and B are assumed to be symmetric,
* and B has to be positive definite.
+ * 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, eigenvalues can
- * be retrieved using the
- * eigenvalue() function. The
- * matrix itself will be
- * LAPACKSupport::unusable after
- * this operation. The number of
- * computed eigenvectors is equal
- * to eigenvectors.size()
+ * been called, all eigenvalues in
+ * (lower_bound, upper_bound] will be
+ * stored in eigenvalues and the
+ * corresponding eigenvectors will be stored
+ * in eigenvectors, whose dimension is set
+ * accordingly.
*
- * Note that the function does
+ * @note Calls the LAPACK
+ * function Xsygvx. For this to
+ * work, ./configure has to
+ * be told to use LAPACK.
+ */
+ void compute_generalized_eigenvalues_symmetric(
+ LAPACKFullMatrix<number> & B,
+ const number lower_bound,
+ const number upper_bound,
+ const number abs_accuracy,
+ Vector<number> & eigenvalues,
+ std::vector<Vector<number> > & eigenvectors,
+ const int itype = 1);
+
+ /**
+ * Same as the other
+ * compute_generalized_eigenvalues_symmetric
+ * function except that all
+ * eigenvalues are computed
+ * and the tolerance is set
+ * automatically.
+ * Note that this function does
* not return the computed
* eigenvalues right away since
* that involves copying data
* eigenvalues and have a
* separate function that returns
* whatever eigenvalue is
- * requested.
+ * requested. Eigenvalues can
+ * be retrieved using the
+ * eigenvalue() function.
+ * The number of computed
+ * eigenvectors is equal
+ * to eigenvectors.size()
*
* @note Calls the LAPACK
* function Xsygv. For this to
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);
-// General eigenvalues and eigenvectors of
+// 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
void dsygv_ (const int* itype, const char* jobz, const char* uplo,
const int* n, double* A, const int* lda, double* B,
const int* ldb, double* w, double* work,
const int* lwork, int* info);
+// Same functionality as dsygv_ 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 dsygvx_ (const int* itype, const char* jobz, const char* range,
+ const char* uplo, const int* n, double* A, const int* lda,
+ double* B, const int* ldb, 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);
// Compute singular value decomposition using divide and conquer
void dgesdd_ (const char* jobz,
}
+template <typename number>
+void
+LAPACKFullMatrix<number>::compute_generalized_eigenvalues_symmetric(
+ LAPACKFullMatrix<number> & B,
+ const number lower_bound,
+ const number upper_bound,
+ const number abs_accuracy,
+ Vector<number> & eigenvalues,
+ std::vector<Vector<number> > & eigenvectors,
+ const int itype)
+{
+ 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());
+ Assert(B.n_rows() == B.n_cols(), ExcNotQuadratic());
+ Assert(static_cast<unsigned int>(nn) == B.n_cols(),
+ ExcDimensionMismatch (nn, B.n_cols()));
+
+ wr.resize(nn);
+ LAPACKFullMatrix<number> matrix_eigenvectors(nn, nn);
+
+ number* values_A = const_cast<number*> (this->data());
+ number* values_B = const_cast<number*> (B.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 ?SYGVX 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);
+
+ sygvx (&itype, jobz, range, uplo, &nn, values_A, &nn,
+ values_B, &nn, &lower_bound, &upper_bound,
+ dummy, dummy, &abs_accuracy, &n_eigenpairs,
+ &wr[0], values_eigenvectors, &nn, &work[0],
+ &lwork, &iwork[0], &ifail[0], &info);
+ // sygvx 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 generalized
+ // eigenvalues.
+ sygvx (&itype, jobz, range, uplo, &nn, values_A, &nn,
+ values_B, &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 sygv" << std::endl;
+
+ eigenvalues.reinit(n_eigenpairs);
+ eigenvectors.resize(n_eigenpairs);
+ for(unsigned int i=0; i < static_cast<unsigned int> (n_eigenpairs); ++i)
+ {
+ eigenvalues(i) = wr[i];
+ unsigned int col_begin(i*nn);
+ eigenvectors[i].reinit(nn, true);
+ for (unsigned int j=0; j < static_cast<unsigned int> (nn); ++j)
+ {
+ eigenvectors[i](j) = values_eigenvectors[col_begin+j];
+ }
+ }
+
+ state = LAPACKSupport::State(unusable);
+}
+
+
template <typename number>
void
LAPACKFullMatrix<number>::compute_generalized_eigenvalues_symmetric (