double* work, const int* lwork,
int* info);
+
+// Symmetric tridiagonal matrix
+void dstev_ (const char* jobz, const int* n,
+ double* d, double* e, double* z,
+ const int* ldz, double* work,
+ int* info);
+
// This file was automatically generated from blas.h.in
// See blastemplates in the deal.II contrib directory
//
-// Copyright (C) 2005 by the deal.II authors
+// Copyright (C) 2005 by the deal authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
int* ilo, int* ihi,
double* scale, double* abnrm,
double* rconde, double* rcondv,
- double* work, int* lwork,
+ double* work, const int* lwork,
int* iwork, int* info);
void sgeevx_ (const char* balanc, const char* jobvl, const char* jobvr,
const char* sense,
int* ilo, int* ihi,
float* scale, float* abnrm,
float* rconde, float* rcondv,
- float* work, int* lwork,
+ float* work, const int* lwork,
int* iwork, int* info);
// Compute singular value decomposition
void dgesvd_ (int* jobu, int* jobvt,
float* vt, const int* ldvt,
float* work, const int* lwork,
int* info);
+// Symmetric tridiagonal matrix
+void dstev_ (const char* jobz, const int* n,
+ double* d, double* e, double* z,
+ const int* ldz, double* work,
+ int* info);
+void sstev_ (const char* jobz, const int* n,
+ float* d, float* e, float* z,
+ const int* ldz, float* work,
+ int* info);
}
inline void
-geevx (const char* balanc, const char* jobvl, const char* jobvr, const char* sense, const int* n, double* A, const int* lda, double* lambda_re, double* lambda_im, double* vl, const int* ldvl, double* vr, const int* ldvr, int* ilo, int* ihi, double* scale, double* abnrm, double* rconde, double* rcondv, double* work, int* lwork, int* iwork, int* info)
+geevx (const char* balanc, const char* jobvl, const char* jobvr, const char* sense, const int* n, double* A, const int* lda, double* lambda_re, double* lambda_im, double* vl, const int* ldvl, double* vr, const int* ldvr, int* ilo, int* ihi, double* scale, double* abnrm, double* rconde, double* rcondv, double* work, const int* lwork, int* iwork, int* info)
{
dgeevx_ (balanc,jobvl,jobvr,sense,n,A,lda,lambda_re,lambda_im,vl,ldvl,vr,ldvr,ilo,ihi,scale,abnrm,rconde,rcondv,work,lwork,iwork,info);
}
inline void
-geevx (const char* balanc, const char* jobvl, const char* jobvr, const char* sense, const int* n, float* A, const int* lda, float* lambda_re, float* lambda_im, float* vl, const int* ldvl, float* vr, const int* ldvr, int* ilo, int* ihi, float* scale, float* abnrm, float* rconde, float* rcondv, float* work, int* lwork, int* iwork, int* info)
+geevx (const char* balanc, const char* jobvl, const char* jobvr, const char* sense, const int* n, float* A, const int* lda, float* lambda_re, float* lambda_im, float* vl, const int* ldvl, float* vr, const int* ldvr, int* ilo, int* ihi, float* scale, float* abnrm, float* rconde, float* rcondv, float* work, const int* lwork, int* iwork, int* info)
{
sgeevx_ (balanc,jobvl,jobvr,sense,n,A,lda,lambda_re,lambda_im,vl,ldvl,vr,ldvr,ilo,ihi,scale,abnrm,rconde,rcondv,work,lwork,iwork,info);
}
}
+inline void
+stev (const char* jobz, const int* n, double* d, double* e, double* z, const int* ldz, double* work, int* info)
+{
+ dstev_ (jobz,n,d,e,z,ldz,work,info);
+}
+
+
+inline void
+stev (const char* jobz, const int* n, float* d, float* e, float* z, const int* ldz, float* work, int* info)
+{
+ sstev_ (jobz,n,d,e,z,ldz,work,info);
+}
+
+
#endif
* A quadratic tridiagonal matrix. That is, a matrix where all entries
* are zero, except the diagonal and the entries left and right of it.
*
- * @note Only data management and entry functions are implemented
- * directly. All more complex functions require LAPACK support.
- */
+ * The matrix has an additional symmetric mode, in which case only the
+ * upper triangle of the matrix is stored and mirrored to the lower
+ * one for matrix vector operations.
+ *
+ * @ingroup Matrix1
+ * @author Guido Kanschat, 2005, 2006
+*/
template<typename number>
class TridiagonalMatrix
{
* empty matrix of dimension
* <tt>n</tt>.
*/
- TridiagonalMatrix(unsigned int n = 0);
+ TridiagonalMatrix(unsigned int n = 0,
+ bool symmetric = false);
+
+ /**
+ * Reinitialize the matrix to a
+ * new size and reset all entries
+ * to zero. The symmetry
+ * properties may be set as well.
+ */
+ void reinit(unsigned int n,
+ bool symmetric = false);
+
//@}
///@name Non-modifying operators
* value. This is restricted to
* the case where <i>|i-j| <=
* 1</i>.
+ *
+ * @note In case of symmetric
+ * storage technique, the entries
+ * <i>(i,j)</i> and <i>(j,i)</i>
+ * are identified and <b>both</b>
+ * exist. This must be taken into
+ * account if adding up is used
+ * for matrix assembling in order
+ * not to obtain doubled entries.
*/
number& operator()(unsigned int i, unsigned int j);
*/
number relative_symmetry_norm2 () const;
//@}
+///@name LAPACK operations
+//@{
+ /**
+ * Return the eigenvalues of the
+ * matrix in the vector provided.
+ *
+ * @note This function requires
+ * configuration of deal.II with
+ * LAPACK support. Additionally,
+ * the matrix must use symmetric
+ * storage technique.
+ */
+ void eigenvalues(Vector<number>& eigenvalues) const;
+//@}
///@name Miscellanea
-//@{
+//@{
/**
* Output of the matrix in
* user-defined format.
* the diagonal. Therefore, the
* length of this vector is the
* same as that of #diagonal.
+ *
+ * The length of this vector is
+ * zero for symmetric storage. In
+ * this case, the second element
+ * of #left is identified with
+ * the first element of #right.
*/
std::vector<number> left;
/**
* same as that of #diagonal.
*/
std::vector<number> right;
+
+ /**
+ * If this flag is true, only the
+ * entries to the right of the
+ * diagonal are stored and the
+ * matrix is assumed symmetric.
+ */
+ bool is_symmetric;
};
/**@}*/
if (j==i)
return diagonal[i];
if (j==i-1)
- return left[i];
+ if (is_symmetric)
+ return right[i-1];
+ else
+ return left[i];
if (j==i+1)
return right[i];
AssertThrow(false, ExcInternalError());
if (j==i)
return diagonal[i];
if (j==i-1)
- return left[i];
+ if (is_symmetric)
+ return right[i-1];
+ else
+ return left[i];
if (j==i+1)
return right[i];
AssertThrow(false, ExcInternalError());
// $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
#include <lac/vector.h>
template<typename number>
-TridiagonalMatrix<number>::TridiagonalMatrix(unsigned int size)
+TridiagonalMatrix<number>::TridiagonalMatrix(
+ unsigned int size,
+ bool symmetric)
:
diagonal(size, 0.),
- left(size, 0.),
- right(size, 0.)
+ left((symmetric ? 0 : size), 0.),
+ right(size, 0.),
+ is_symmetric(symmetric)
{}
+template<typename number>
+void
+TridiagonalMatrix<number>::reinit(
+ unsigned int size,
+ bool symmetric)
+{
+ is_symmetric = symmetric;
+ diagonal.resize(size);
+ right.resize(size);
+ left.resize(symmetric ? 0 : size);
+}
+
+
template<typename number>
bool
TridiagonalMatrix<number>::all_zero() const
Assert(v.size() == n(), ExcDimensionMismatch(v.size(), n()));
if (n()==0) return;
+
+ // The actual loop skips the first
+ // and last row
+ const unsigned int e=n()-1;
+ // Let iterators point to the first
+ // entry of each diagonal
+ typename std::vector<number>::const_iterator d = diagonal.begin();
+ typename std::vector<number>::const_iterator r = right.begin();
+ // The left diagonal starts one
+ // later or is equal to the right
+ // one for symmetric storage
+ typename std::vector<number>::const_iterator l = left.begin();
+ if (is_symmetric)
+ l = r;
+ else
+ ++l;
if (adding)
{
- w(0) += diagonal[0]*v(0) + right[0]*v(1);
- const unsigned int e=n()-1;
- for (unsigned int i=1;i<e;++i)
- w(i) += left[i]*v(i-1)+diagonal[i]*v(i)+right[i]*v(i+1);
- w(e) += left[e]*v(e-1)+diagonal[e]*v(e);
+ // Treat first row separately
+ w(0) += (*d) * v(0) + (*r) * v(1);
+ ++d; ++r;
+ // All rows with three entries
+ for (unsigned int i=1;i<e;++i,++d,++r,++l)
+ w(i) += (*l) * v(i-1) + (*d) * v(i) + (*r) * v(i+1);
+ // Last row is special again
+ w(e) += (*l) * v(e-1) + (*d) * v(e);
}
else
{
- w(0) = diagonal[0]*v(0) + right[0]*v(1);
- const unsigned int e=n()-1;
- for (unsigned int i=1;i<e;++i)
- w(i) = left[i]*v(i-1)+diagonal[i]*v(i)+right[i]*v(i+1);
- w(e) = left[e]*v(e-1)+diagonal[e]*v(e);
+ w(0) = (*d) * v(0) + (*r) * v(1);
+ ++d; ++r;
+ for (unsigned int i=1;i<e;++i,++d,++r,++l)
+ w(i) = (*l) * v(i-1) + (*d) * v(i) + (*r) * v(i+1);
+ w(e) = (*l) * v(e-1) + (*d) * v(e);
}
}
Assert(v.size() == n(), ExcDimensionMismatch(v.size(), n()));
if (n()==0) return;
-//TODO:[GK] Check this!!!
+
+ const unsigned int e=n()-1;
+ typename std::vector<number>::const_iterator d = diagonal.begin();
+ typename std::vector<number>::const_iterator r = right.begin();
+ typename std::vector<number>::const_iterator l = left.begin();
+ if (is_symmetric)
+ l = r;
+ else
+ ++l;
+
if (adding)
{
- w(0) += diagonal[0]*v(0) + left[1]*v(1);
- const unsigned int e=n()-1;
- for (unsigned int i=1;i<e;++i)
- w(i) += left[i+1]*v(i+1)+diagonal[i]*v(i)+right[i-1]*v(i-1);
- w(e) += right[e-1]*v(e-1)+diagonal[e]*v(e);
+ w(0) += (*d) * v(0) + (*l) * v(1);
+ ++d; ++l;
+ for (unsigned int i=1;i<e;++i,++d,++r,++l)
+ w(i) += (*l) * v(i+1) + (*d) * v(i) + (*r) * v(i-1);
+ w(e) += (*d) * v(e) + (*r) * v(e-1);
}
else
{
- w(0) = diagonal[0]*v(0) + left[1]*v(1);
- const unsigned int e=n()-1;
- for (unsigned int i=1;i<e;++i)
- w(i) = left[i+1]*v(i+1)+diagonal[i]*v(i)+right[i-1]*v(i-1);
- w(e) = right[e-1]*v(e-1)+diagonal[e]*v(e);
+ w(0) = (*d) * v(0) + (*l) * v(1);
+ ++d; ++l;
+ for (unsigned int i=1;i<e;++i,++d,++r,++l)
+ w(i) = (*l) * v(i+1) + (*d) * v(i) + (*r) * v(i-1);
+ w(e) = (*d) * v(e) + (*r) * v(e-1);
}
}
Tvmult(w, v, true);
}
-/*
+
template<typename number>
-TridiagonalMatrix<number>::
+number
+TridiagonalMatrix<number>::matrix_scalar_product(
+ const Vector<number>& w,
+ const Vector<number>& v) const
{
+ const unsigned int e=n()-1;
+ typename std::vector<number>::const_iterator d = diagonal.begin();
+ typename std::vector<number>::const_iterator r = right.begin();
+ typename std::vector<number>::const_iterator l = left.begin();
+ if (is_symmetric)
+ l = r;
+ else
+ ++l;
+
+ number result = w(0) * ((*d) * v(0) + (*r) * v(1));
+ ++d; ++r;
+ for (unsigned int i=1;i<e;++i,++d,++r,++l)
+ result += w(i) * ((*l) * v(i-1)+ (*d) * v(i)+ (*r) * v(i+1));
+ result += w(e) * ((*l) * v(e-1) + (*d) * v(e));
+ return result;
+}
+
+
+template<typename number>
+number
+TridiagonalMatrix<number>::matrix_norm_square(
+ const Vector<number>& v) const
+{
+ return matrix_scalar_product(v,v);
}
+template<typename number>
+void
+TridiagonalMatrix<number>::eigenvalues(Vector<double>& lambda) const
+{
+ Assert(false, ExcNotImplemented());
+}
+
+/*
template<typename number>
TridiagonalMatrix<number>::
{