]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
add symmetry option for TridiagonalMatrix
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 6 Aug 2006 22:01:07 +0000 (22:01 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 6 Aug 2006 22:01:07 +0000 (22:01 +0000)
git-svn-id: https://svn.dealii.org/trunk@13610 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/contrib/blastemplates/lapack_templates.h.in
deal.II/lac/include/lac/lapack_templates.h
deal.II/lac/include/lac/tridiagonal_matrix.h
deal.II/lac/source/tridiagonal_matrix.cc

index af3b9a98141ea93bfb26c3aa9cd720bff69368d3..1b23d0d64f443da55ecdd9b4d88e87be6f822dc9 100644 (file)
@@ -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);
+
index 9fad4a7386f9b81be714ed58a8c1a66c92655f32..4266d4eb75bca316f04bdd2aa17ab67ed143c44c 100644 (file)
@@ -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
index db6609dc71687733f38725cdb3b1201932b25109..4fdcffecfa7d3b9c18479dcd80bf9cdd0b6b793c 100644 (file)
@@ -33,9 +33,13 @@ template<typename number> 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<typename number>
 class TridiagonalMatrix
 {
@@ -48,7 +52,18 @@ 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
@@ -98,6 +113,15 @@ class TridiagonalMatrix
                                      * 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);
     
@@ -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<number>& 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<number> left;
                                     /**
@@ -295,6 +339,14 @@ class TridiagonalMatrix
                                      * 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;
 };
 
 /**@}*/
@@ -332,7 +384,10 @@ TridiagonalMatrix<number>::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<number>::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());
index 915c873c29f8fa13e71bbd15f5f296c386824c1d..4a3c588824f73dd045db9f1aa9cbaeb7e58a16d1 100644 (file)
@@ -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
 #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
@@ -56,22 +72,41 @@ TridiagonalMatrix<number>::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<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);
     }
 }
 
@@ -97,22 +132,31 @@ TridiagonalMatrix<number>::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<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);
     }
 }
 
@@ -126,13 +170,48 @@ TridiagonalMatrix<number>::Tvmult_add (
   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>::
 {

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.