// $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
unusable = 0x8000
};
+ /**
+ * Function printing the name of a State.
+ */
+ inline const char* state_name(State s)
+ {
+ switch (s)
+ {
+ case matrix:
+ return "matrix";
+ case lu:
+ return "lu decomposition";
+ case eigenvalues:
+ return "eigenvalues";
+ case unusable:
+ return "unusable";
+ default:
+ return "unknown";
+ }
+ return "internal error";
+ }
/**
* A matrix can have certain features allowing for optimization, but
* Integer constant.
*/
static const int one = 1;
+
+ /**
+ * Exception thrown when a matrix
+ * is not in a suitable state for
+ * an operation. For instance, a
+ * LAPACK routine may have left the
+ * matrix in an unusable state,
+ * then vmult does not make sense
+ * anymore.
+ */
+ DeclException1(ExcState, State,
+ << "The function cannot be called while the matrix is in state "
+ << state_name(arg1));
+
}
#include <base/config.h>
#include <base/subscriptor.h>
#include <base/smartpointer.h>
+#include <lac/lapack_support.h>
#include <vector>
#include <iomanip>
///@name LAPACK operations
//@{
/**
- * Return the eigenvalues of the
- * matrix in the vector provided.
+ * Compute the eigenvalues of the
+ * symmetric tridiagonal matrix.
*
* @note This function requires
* configuration of deal.II with
* the matrix must use symmetric
* storage technique.
*/
- void eigenvalues(Vector<number>& eigenvalues) const;
+ void compute_eigenvalues();
+ /**
+ * After calling
+ * compute_eigenvalues(), you can
+ * access each eigenvalue here.
+ */
+ number eigenvalue(const unsigned i) const;
//@}
///@name Miscellanea
//@{
* matrix is assumed symmetric.
*/
bool is_symmetric;
+
+ /**
+ * The state of the
+ * matrix. Normally, the state is
+ * LAPACKSupport::matrix,
+ * indicating that the object can
+ * be used for regular matrix
+ * operations.
+ *
+ * See explanation of this data
+ * type for details.
+ */
+ LAPACKSupport::State state;
};
/**@}*/
#include <lac/tridiagonal_matrix.h>
#include <lac/vector.h>
+#include <lac/lapack_templates.h>
+
+using namespace LAPACKSupport;
template<typename number>
TridiagonalMatrix<number>::TridiagonalMatrix(
diagonal(size, 0.),
left((symmetric ? 0 : size), 0.),
right(size, 0.),
- is_symmetric(symmetric)
+ is_symmetric(symmetric),
+ state(matrix)
{}
diagonal.resize(size);
right.resize(size);
left.resize(symmetric ? 0 : size);
+ state = matrix;
}
bool
TridiagonalMatrix<number>::all_zero() const
{
+ Assert(state == matrix, ExcState(state));
+
typename std::vector<number>::const_iterator i;
typename std::vector<number>::const_iterator e;
const Vector<number> &v,
const bool adding) const
{
+ Assert(state == matrix, ExcState(state));
+
Assert(w.size() == n(), ExcDimensionMismatch(w.size(), n()));
Assert(v.size() == n(), ExcDimensionMismatch(v.size(), n()));
const Vector<number> &v,
const bool adding) const
{
+ Assert(state == matrix, ExcState(state));
+
Assert(w.size() == n(), ExcDimensionMismatch(w.size(), n()));
Assert(v.size() == n(), ExcDimensionMismatch(v.size(), n()));
const Vector<number>& w,
const Vector<number>& v) const
{
+ Assert(state == matrix, ExcState(state));
+
const unsigned int e=n()-1;
typename std::vector<number>::const_iterator d = diagonal.begin();
typename std::vector<number>::const_iterator r = right.begin();
template<typename number>
void
-TridiagonalMatrix<number>::eigenvalues(Vector<number>&) const
+TridiagonalMatrix<number>::compute_eigenvalues()
+{
+#ifdef HAVE_LIBLAPACK
+ Assert(state == matrix, ExcState(state));
+ Assert(is_symmetric, ExcNotImplemented());
+
+ const int nn = n();
+ int info;
+ stev (&N, &nn, &*diagonal.begin(), &*right.begin(), 0, &one, 0, &info);
+ Assert(info == 0, ExcInternalError());
+
+ state = eigenvalues;
+#else
+ Assert(false, ExcNeedsLAPACK());
+#endif
+}
+
+
+template<typename number>
+number
+TridiagonalMatrix<number>::eigenvalue(const unsigned int i) const
{
- Assert(false, ExcNotImplemented());
+ Assert(state == eigenvalues, ExcState(state));
+ Assert(i<n(), ExcIndexRange(i,0,n()));
+ return diagonal[i];
}
+
/*
template<typename number>
TridiagonalMatrix<number>::