From: kanschat Date: Sun, 6 Aug 2006 22:01:07 +0000 (+0000) Subject: add symmetry option for TridiagonalMatrix X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=aadea6075b7f750db885b0f1313378530f654205;p=dealii-svn.git add symmetry option for TridiagonalMatrix git-svn-id: https://svn.dealii.org/trunk@13610 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/contrib/blastemplates/lapack_templates.h.in b/deal.II/contrib/blastemplates/lapack_templates.h.in index af3b9a9814..1b23d0d64f 100644 --- a/deal.II/contrib/blastemplates/lapack_templates.h.in +++ b/deal.II/contrib/blastemplates/lapack_templates.h.in @@ -33,3 +33,10 @@ void dgesvd_ (int* jobu, int* jobvt, 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); + diff --git a/deal.II/lac/include/lac/lapack_templates.h b/deal.II/lac/include/lac/lapack_templates.h index 9fad4a7386..4266d4eb75 100644 --- a/deal.II/lac/include/lac/lapack_templates.h +++ b/deal.II/lac/include/lac/lapack_templates.h @@ -5,7 +5,7 @@ // 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 @@ -54,7 +54,7 @@ void dgeevx_ (const char* balanc, const char* jobvl, const char* jobvr, 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, @@ -65,7 +65,7 @@ void sgeevx_ (const char* balanc, const char* jobvl, const char* jobvr, 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, @@ -82,6 +82,15 @@ void sgesvd_ (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); } @@ -116,14 +125,14 @@ geev (const char* jobvl, const char* jobvr, const int* n, float* A, const int* l 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); } @@ -143,4 +152,18 @@ gesvd (int* jobu, int* jobvt, const int* n, const int* m, float* A, const int* l } +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 diff --git a/deal.II/lac/include/lac/tridiagonal_matrix.h b/deal.II/lac/include/lac/tridiagonal_matrix.h index db6609dc71..4fdcffecfa 100644 --- a/deal.II/lac/include/lac/tridiagonal_matrix.h +++ b/deal.II/lac/include/lac/tridiagonal_matrix.h @@ -33,9 +33,13 @@ template class Vector; * 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 class TridiagonalMatrix { @@ -48,7 +52,18 @@ class TridiagonalMatrix * empty matrix of dimension * n. */ - 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 @@ -98,6 +113,15 @@ class TridiagonalMatrix * value. This is restricted to * the case where |i-j| <= * 1. + * + * @note In case of symmetric + * storage technique, the entries + * (i,j) and (j,i) + * are identified and both + * 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); @@ -251,8 +275,22 @@ class TridiagonalMatrix */ 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& eigenvalues) const; +//@} ///@name Miscellanea -//@{ +//@{ /** * Output of the matrix in * user-defined format. @@ -283,6 +321,12 @@ class TridiagonalMatrix * 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 left; /** @@ -295,6 +339,14 @@ class TridiagonalMatrix * same as that of #diagonal. */ std::vector 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; }; /**@}*/ @@ -332,7 +384,10 @@ TridiagonalMatrix::operator()(unsigned int i, unsigned int j) const 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()); @@ -353,7 +408,10 @@ TridiagonalMatrix::operator()(unsigned int i, unsigned int j) 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()); diff --git a/deal.II/lac/source/tridiagonal_matrix.cc b/deal.II/lac/source/tridiagonal_matrix.cc index 915c873c29..4a3c588824 100644 --- a/deal.II/lac/source/tridiagonal_matrix.cc +++ b/deal.II/lac/source/tridiagonal_matrix.cc @@ -2,7 +2,7 @@ // $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 @@ -15,14 +15,30 @@ #include template -TridiagonalMatrix::TridiagonalMatrix(unsigned int size) +TridiagonalMatrix::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 +void +TridiagonalMatrix::reinit( + unsigned int size, + bool symmetric) +{ + is_symmetric = symmetric; + diagonal.resize(size); + right.resize(size); + left.resize(symmetric ? 0 : size); +} + + template bool TridiagonalMatrix::all_zero() const @@ -56,22 +72,41 @@ TridiagonalMatrix::vmult ( 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::const_iterator d = diagonal.begin(); + typename std::vector::const_iterator r = right.begin(); + // The left diagonal starts one + // later or is equal to the right + // one for symmetric storage + typename std::vector::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::Tvmult ( 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::const_iterator d = diagonal.begin(); + typename std::vector::const_iterator r = right.begin(); + typename std::vector::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::Tvmult_add ( Tvmult(w, v, true); } -/* + template -TridiagonalMatrix:: +number +TridiagonalMatrix::matrix_scalar_product( + const Vector& w, + const Vector& v) const { + const unsigned int e=n()-1; + typename std::vector::const_iterator d = diagonal.begin(); + typename std::vector::const_iterator r = right.begin(); + typename std::vector::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 +number +TridiagonalMatrix::matrix_norm_square( + const Vector& v) const +{ + return matrix_scalar_product(v,v); } +template +void +TridiagonalMatrix::eigenvalues(Vector& lambda) const +{ + Assert(false, ExcNotImplemented()); +} + +/* template TridiagonalMatrix:: {