]> https://gitweb.dealii.org/ - dealii.git/commitdiff
starting LAPACK matrix
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Wed, 23 Mar 2005 03:49:05 +0000 (03:49 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Wed, 23 Mar 2005 03:49:05 +0000 (03:49 +0000)
git-svn-id: https://svn.dealii.org/trunk@10211 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/lapack_full_matrix.h [new file with mode: 0644]
deal.II/lac/include/lac/lapack_support.h [new file with mode: 0644]
deal.II/lac/include/lac/lapack_templates.h [new file with mode: 0644]
deal.II/lac/include/lac/vector.h
deal.II/lac/source/lapack_full_matrix.cc [new file with mode: 0644]

diff --git a/deal.II/lac/include/lac/lapack_full_matrix.h b/deal.II/lac/include/lac/lapack_full_matrix.h
new file mode 100644 (file)
index 0000000..472d15a
--- /dev/null
@@ -0,0 +1,251 @@
+//-------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998 - 2005 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-------------------------------------------------------------------
+#ifndef __deal2__lapack_full_matrix_h
+#define __deal2__lapack_full_matrix_h
+
+#include <base/table.h>
+#include <lac/lapack_support.h>
+
+// forward declarations
+template<typename number> class Vector;
+template<typename number> class FullMatrix;
+
+/*! @addtogroup Matrix1
+ *@{
+ */
+
+
+/**
+ * A variant of FullMatrix using LAPACK functions whereever possible.
+ *
+ * @author Guido Kanschat, 2005
+ */
+template <typename number>
+class LAPACKFullMatrix : public TransposeTable<number>
+{
+  public:
+                                    /**
+                                     * Constructor. Initialize the
+                                     * matrix as a square matrix with
+                                     * dimension <tt>n</tt>.
+                                     *
+                                     * In order to avoid the implicit
+                                     * conversion of integers and
+                                     * other types to a matrix, this
+                                     * constructor is declared
+                                     * <tt>explicit</tt>.
+                                     *
+                                     * By default, no memory is
+                                     * allocated.
+                                     */
+    explicit LAPACKFullMatrix (const unsigned int n = 0);
+    
+                                    /**
+                                     * Constructor. Initialize the
+                                     * matrix as a rectangular
+                                     * matrix.
+                                     */
+    LAPACKFullMatrix (const unsigned int rows,
+                     const unsigned int cols);
+    
+                                    /** 
+                                     * Copy constructor. This
+                                     * constructor does a deep copy
+                                     * of the matrix. Therefore, it
+                                     * poses a possible efficiency
+                                     * problem, if for example,
+                                     * function arguments are passed
+                                     * by value rather than by
+                                     * reference. Unfortunately, we
+                                     * can't mark this copy
+                                     * constructor <tt>explicit</tt>,
+                                     * since that prevents the use of
+                                     * this class in containers, such
+                                     * as <tt>std::vector</tt>. The
+                                     * responsibility to check
+                                     * performance of programs must
+                                     * therefore remain with the
+                                     * user of this class.
+                                     */
+    LAPACKFullMatrix (const LAPACKFullMatrix&);
+
+                                    /**
+                                     * Assignment operator.
+                                     */
+    LAPACKFullMatrix<number> &
+    operator = (const LAPACKFullMatrix<number>&);
+    
+                                    /**
+                                     * Assignment operator for a
+                                     * regular FullMatrix.
+                                     */
+    template <typename number2>
+    LAPACKFullMatrix<number> &
+    operator = (const FullMatrix<number2>&);
+    
+                                    /**
+                                     * This operator assigns a scalar
+                                     * to a matrix. To avoid
+                                     * confusion with constructors,
+                                     * zero is the only value allowed
+                                     * for <tt>d</tt>
+                                     */
+    LAPACKFullMatrix<number> &
+    operator = (const double d);
+
+                                     /**
+                                     * Assignment from different
+                                     * matrix classes. This
+                                     * assignment operator uses
+                                     * iterators of the class
+                                     * MATRIX. Therefore, sparse
+                                     * matrices are possible sources.
+                                     */
+    template <class MATRIX>
+    void copy_from (const MATRIX&);
+    
+                                    /**
+                                     * Fill rectangular block.
+                                     *
+                                     * A rectangular block of the
+                                     * matrix <tt>src</tt> is copied into
+                                     * <tt>this</tt>. The upper left
+                                     * corner of the block being
+                                     * copied is
+                                     * <tt>(src_offset_i,src_offset_j)</tt>.
+                                     * The upper left corner of the
+                                     * copied block is
+                                     * <tt>(dst_offset_i,dst_offset_j)</tt>.
+                                     * The size of the rectangular
+                                     * block being copied is the
+                                     * maximum size possible,
+                                     * determined either by the size
+                                     * of <tt>this</tt> or <tt>src</tt>.
+                                     */
+    template<class MATRIX>
+    void fill (const MATRIX &src,
+              const unsigned int dst_offset_i = 0,
+              const unsigned int dst_offset_j = 0,
+              const unsigned int src_offset_i = 0,
+              const unsigned int src_offset_j = 0);
+    
+                                    /**
+                                     * Matrix-vector-multiplication.
+                                     *
+                                     * The optional parameter
+                                     * <tt>adding</tt> determines, whether the
+                                     * result is stored in <tt>w</tt> or added
+                                     * to <tt>w</tt>.
+                                     *
+                                     * if (adding)
+                                     *  <i>w += A*v</i>
+                                     *
+                                     * if (!adding)
+                                     *  <i>w = A*v</i>
+                                      *
+                                      * Source and destination must
+                                      * not be the same vector.
+                                     */
+    void vmult (Vector<number>   &w,
+               const Vector<number> &v,
+               const bool            adding=false) const;
+                                    /**
+                                     * Adding Matrix-vector-multiplication.
+                                     *  <i>w += A*v</i>
+                                      *
+                                      * Source and destination must
+                                      * not be the same vector.
+                                     */
+    void vmult_add (Vector<number>       &w,
+                   const Vector<number> &v) const;
+    
+                                    /**
+                                     * Transpose
+                                     * matrix-vector-multiplication.
+                                     *
+                                     * The optional parameter
+                                     * <tt>adding</tt> determines, whether the
+                                     * result is stored in <tt>w</tt> or added
+                                     * to <tt>w</tt>.
+                                     *
+                                     * if (adding)
+                                     *  <i>w += A<sup>T</sup>*v</i>
+                                     *
+                                     * if (!adding)
+                                     *  <i>w = A<sup>T</sup>*v</i>
+                                      *
+                                      *
+                                      * Source and destination must
+                                      * not be the same vector.
+                                     */
+    void Tvmult (Vector<number>       &w,
+                const Vector<number> &v,
+                const bool            adding=false) const;
+
+                                    /**
+                                     * Adding transpose
+                                     * matrix-vector-multiplication.
+                                     *  <i>w += A<sup>T</sup>*v</i>
+                                      *
+                                      * Source and destination must
+                                      * not be the same vector.
+                                     */
+    void Tvmult_add (Vector<number>       &w,
+                    const Vector<number> &v) const;
+
+};
+
+
+template <typename number>
+template <class MATRIX>
+void
+LAPACKFullMatrix<number>::copy_from (const MATRIX& M)
+{
+  this->reinit (M.m(), M.n());
+  const typename MATRIX::const_iterator end = M.end();
+  for (typename MATRIX::const_iterator entry = M.begin();
+       entry != end; ++entry)
+    this->el(entry->row(), entry->column()) = entry->value();
+}
+
+
+
+template <typename number>
+template <class MATRIX>
+void
+LAPACKFullMatrix<number>::fill (
+  const MATRIX& M,
+  const unsigned int dst_offset_i,
+  const unsigned int dst_offset_j,
+  const unsigned int src_offset_i,
+  const unsigned int src_offset_j)
+{
+  const unsigned int endrow = src_offset_i + m();
+  const unsigned int endcol = src_offset_j + n();
+
+  for (unsigned int i=0;i<m();++i)
+    {
+      const typename MATRIX::const_iterator end = M.end(src_offset_i+i);
+      for (typename MATRIX::const_iterator entry = M.begin(src_offset_i+i);
+          entry != end; ++entry)
+       if (entry->column() < endcol)
+         this->el(dst_offset_i+i,
+                  dst_offset_j+entry->column()-src_offset_j)
+           = entry->value();
+    }
+}
+
+
+
+
+#endif
diff --git a/deal.II/lac/include/lac/lapack_support.h b/deal.II/lac/include/lac/lapack_support.h
new file mode 100644 (file)
index 0000000..45b7792
--- /dev/null
@@ -0,0 +1,76 @@
+//-------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2005 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-------------------------------------------------------------------
+#ifndef __deal2__lapack_support_h
+#define __deal2__lapack_support_h
+
+namespace LAPACKSupport
+{
+/**
+ * Most LAPACK functions change the contents of the matrix applied to
+ * to something which is not a matrix anymore. Therefore, LAPACK
+ * matrix classes in <tt>deal.II</tt> have a state flag indicating
+ * what happened to them.
+ *
+ * @author Guido Kanschat, 2005
+ */
+  enum State
+  {
+                                        /// Contents is something useless.
+       unusable,
+                                        /// Contents is actually a matrix.
+       matrix,
+                                        /// Contents is an LU decomposition.
+       lu
+  };
+  
+  
+/**
+ * A matrix can have certain features allowing for optimization, but
+ * hard to test. These are listed here.
+ */
+  enum Properties
+  {
+                                        /// No special properties
+       general = 0,
+                                        /// Matrix is symmetric
+       symmetric = 1,
+                                        /// Matrix is upper triangular
+       upper_triangle = 2,
+                                        /// Matrix is lower triangular
+       lower_triangle = 4,
+                                        /// Matrix is diagonal
+       diagonal = 6,
+                                        /// Matrix is in upper Hessenberg form
+       hessenberg = 8
+  };
+  
+                                  /**
+                                   * Character constant.
+                                   */
+  extern const char N = 'N';
+                                  /**
+                                   * Character constant.
+                                   */
+  extern const char T = 'T';
+                                  /**
+                                   * Interger constant.
+                                   */
+  extern const int zero = 0;
+                                  /**
+                                   * Interger constant.
+                                   */
+  extern const int one = 1;
+};
+
+
+#endif
diff --git a/deal.II/lac/include/lac/lapack_templates.h b/deal.II/lac/include/lac/lapack_templates.h
new file mode 100644 (file)
index 0000000..26c20cf
--- /dev/null
@@ -0,0 +1,144 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    This file was automatically generated from blas.h.in
+//    See blastemplates in the deal.II contrib directory
+//
+//    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
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#ifndef __LAPACK_TEMPLATES_H
+#define __LAPACK_TEMPLATES_H
+
+extern "C"
+{
+// General Matrix
+// Matrix vector product
+void dgemv_ (const char* trans, const int* m, const int* n,
+            const double* alpha, const double* A, const int* lda,
+            const double* x, const int* incx,
+            const double* b, double* y, const int* incy);
+void sgemv_ (const char* trans, const int* m, const int* n,
+            const float* alpha, const float* A, const int* lda,
+            const float* x, const int* incx,
+            const float* b, float* y, const int* incy);
+// Compute eigenvalues and vectors
+void dgeev_ (char* jobvl, char* jobvr,
+            int* n, double* A, int* lda,
+            double* lambda_re, double* lambda_im,
+            double* vl, int* ldvl,
+            double* vr, int* ldva,
+            double* work, int* lwork,
+            int* info);
+void sgeev_ (char* jobvl, char* jobvr,
+            int* n, float* A, int* lda,
+            float* lambda_re, float* lambda_im,
+            float* vl, int* ldvl,
+            float* vr, int* ldva,
+            float* work, int* lwork,
+            int* info);
+// Compute eigenvalues and vectors (expert)
+void dgeevx_ (char* balanc, char* jobvl, char* jobvr, char* sense,
+             int* n, double* A, int* lda,
+             double* lambda_re, double* lambda_im,
+             double* vl, int* ldvl,
+             double* vr, int* ldvr,
+             int* ilo, int* ihi,
+             double* scale, double* abnrm,
+             double* rconde, double* rcondv,
+             double* work, int* lwork,
+             int* iwork, int* info);
+void sgeevx_ (char* balanc, char* jobvl, char* jobvr, char* sense,
+             int* n, float* A, int* lda,
+             float* lambda_re, float* lambda_im,
+             float* vl, int* ldvl,
+             float* vr, int* ldvr,
+             int* ilo, int* ihi,
+             float* scale, float* abnrm,
+             float* rconde, float* rcondv,
+             float* work, int* lwork,
+             int* iwork, int* info);
+// Compute singular value decomposition
+void dgesvd_ (int* jobu, int* jobvt,
+             int* n, int* m, double* A, int* lda,
+             double* s,
+             double* u, int* ldu,
+             double* vt, int* ldvt,
+             double* work, int* lwork,
+             int* info);
+void sgesvd_ (int* jobu, int* jobvt,
+             int* n, int* m, float* A, int* lda,
+             float* s,
+             float* u, int* ldu,
+             float* vt, int* ldvt,
+             float* work, int* lwork,
+             int* info);
+
+}
+
+
+
+inline void
+gemv (const char* trans, const int* m, const int* n, const double* alpha, const double* A, const int* lda, const double* x, const int* incx, const double* b, double* y, const int* incy)
+{
+  dgemv_ (trans,m,n,alpha,A,lda,x,incx,b,y,incy);
+}
+
+
+inline void
+gemv (const char* trans, const int* m, const int* n, const float* alpha, const float* A, const int* lda, const float* x, const int* incx, const float* b, float* y, const int* incy)
+{
+  sgemv_ (trans,m,n,alpha,A,lda,x,incx,b,y,incy);
+}
+
+
+inline void
+geev (char* jobvl, char* jobvr, int* n, double* A, int* lda, double* lambda_re, double* lambda_im, double* vl, int* ldvl, double* vr, int* ldva, double* work, int* lwork, int* info)
+{
+  dgeev_ (jobvl,jobvr,n,A,lda,lambda_re,lambda_im,vl,ldvl,vr,ldva,work,lwork,info);
+}
+
+
+inline void
+geev (char* jobvl, char* jobvr, int* n, float* A, int* lda, float* lambda_re, float* lambda_im, float* vl, int* ldvl, float* vr, int* ldva, float* work, int* lwork, int* info)
+{
+  sgeev_ (jobvl,jobvr,n,A,lda,lambda_re,lambda_im,vl,ldvl,vr,ldva,work,lwork,info);
+}
+
+
+inline void
+geevx (char* balanc, char* jobvl, char* jobvr, char* sense, int* n, double* A, int* lda, double* lambda_re, double* lambda_im, double* vl, int* ldvl, double* vr, int* ldvr, int* ilo, int* ihi, double* scale, double* abnrm, double* rconde, double* rcondv, double* work, 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 (char* balanc, char* jobvl, char* jobvr, char* sense, int* n, float* A, int* lda, float* lambda_re, float* lambda_im, float* vl, int* ldvl, float* vr, int* ldvr, int* ilo, int* ihi, float* scale, float* abnrm, float* rconde, float* rcondv, float* work, 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
+gesvd (int* jobu, int* jobvt, int* n, int* m, double* A, int* lda, double* s, double* u, int* ldu, double* vt, int* ldvt, double* work, int* lwork, int* info)
+{
+  dgesvd_ (jobu,jobvt,n,m,A,lda,s,u,ldu,vt,ldvt,work,lwork,info);
+}
+
+
+inline void
+gesvd (int* jobu, int* jobvt, int* n, int* m, float* A, int* lda, float* s, float* u, int* ldu, float* vt, int* ldvt, float* work, int* lwork, int* info)
+{
+  sgesvd_ (jobu,jobvt,n,m,A,lda,s,u,ldu,vt,ldvt,work,lwork,info);
+}
+
+
+#endif
index 9bd143ed07f733ea9178c4b709e951b9fb99d4ed..98f779006475c652656328c0f517b100115c46ce 100644 (file)
@@ -1,15 +1,15 @@
-//----------------------------  vector.h  ---------------------------
+//---------------------------------------------------------------------------
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 by the deal.II authors
+//    Copyright (C) 1998 - 2005 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
 //    to the file deal.II/doc/license.html for the  text  and
 //    further information on this license.
 //
-//----------------------------  vector.h  ---------------------------
+//---------------------------------------------------------------------------
 #ifndef __deal2__vector_h
 #define __deal2__vector_h
 
@@ -31,7 +31,7 @@ namespace PETScWrappers
 }
 #endif
 
-
+template<typename number> class LAPACKFullMatrix;
 
 /*! @addtogroup Vectors
  *@{
@@ -748,11 +748,16 @@ class Vector
                                      */
     Number *val;
 
-                                     /**
+                                     /*
                                       * Make all other vector types
                                       * friends.
                                       */
-    template <typename Number2> friend class Vector;    
+    template <typename Number2> friend class Vector;
+                                    /*
+                                     * LAPACK matrices need access to
+                                     * the data.
+                                     */
+    friend class LAPACKFullMatrix<Number>;
 };
 
 /*@}*/
diff --git a/deal.II/lac/source/lapack_full_matrix.cc b/deal.II/lac/source/lapack_full_matrix.cc
new file mode 100644 (file)
index 0000000..830fb88
--- /dev/null
@@ -0,0 +1,174 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2005 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#include <lac/lapack_full_matrix.h>
+#include <lac/lapack_templates.h>
+#include <lac/lapack_support.h>
+#include <lac/vector.h>
+
+using namespace LAPACKSupport;
+
+template <typename number>
+LAPACKFullMatrix<number>::LAPACKFullMatrix(const unsigned int n)
+               :
+               TransposeTable<number> (n,n)
+{}
+
+
+template <typename number>
+LAPACKFullMatrix<number>::LAPACKFullMatrix(
+  const unsigned int m,
+  const unsigned int n)
+               :
+               TransposeTable<number> (m,n)
+{}
+
+
+template <typename number>
+LAPACKFullMatrix<number>::LAPACKFullMatrix(const LAPACKFullMatrix &M)
+               :
+               TransposeTable<number> (M)
+{}
+
+
+template <typename number>
+LAPACKFullMatrix<number> &
+LAPACKFullMatrix<number>::operator = (const LAPACKFullMatrix<number>& M)
+{
+  TransposeTable<number>::operator=(M);
+  return *this;
+}
+
+
+template <typename number>
+template <typename number2>
+LAPACKFullMatrix<number> &
+LAPACKFullMatrix<number>::operator = (const FullMatrix<number2>& M)
+{
+  Assert (n() == M.m(), ExcDimensionMismatch(n(), M.m()));
+  Assert (n() == M.n(), ExcDimensionMismatch(n(), M.n()));
+  for (unsigned int i=0;i<m;++i)
+    for (unsigned int j=0;j<n;++j)
+      (*this)(i,j) = M(i,j);
+  return *this;
+}
+
+
+template <typename number>
+LAPACKFullMatrix<number> &
+LAPACKFullMatrix<number>::operator = (const double d)
+{
+  Assert (d==0, ExcScalarAssignmentOnlyForZeroValue());
+  
+  if (this->n_elements() != 0)
+    std::fill_n (this->val, this->n_elements(), number());
+  
+  return *this;
+}
+
+#ifdef HAVE_LIBLAPACK
+
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::vmult (
+  Vector<number>       &w,
+  const Vector<number> &v,
+  const bool            adding) const
+{
+  const int mm = this->n_rows();
+  const int nn = this->n_cols();
+  const number alpha = 1.;
+  const number beta = (adding ? 1. : 0.);
+  
+  gemv(&N, &mm, &nn, &alpha, data(), &mm, v.val, &one, &beta, w.val, &one);
+}
+
+
+template <typename number>
+LAPACKFullMatrix<number>::Tvmult (
+  Vector<number>       &w,
+  const Vector<number> &v,
+  const bool            adding) const
+{
+  const int mm = this->n_rows();
+  const int nn = this->n_cols();
+  const number alpha = 1.;
+  const number beta = (adding ? 1. : 0.);
+  
+  gemv(&T, &mm, &nn, &alpha, data(), &mm, v.val, &one, &beta, w.val, &one);
+}
+
+#else
+
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::vmult (
+  Vector<number>       &,
+  const Vector<number> &,
+  const bool            ) const
+{
+  Assert(false, ExcNeedsBLAS());
+}
+
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::Tvmult (
+  Vector<number>       &,
+  const Vector<number> &,
+  const bool            ) const
+{
+  Assert(false, ExcNeedsBLAS());
+}
+
+
+#endif
+
+// template <typename number>
+// LAPACKFullMatrix<number>::()
+// {}
+
+
+// template <typename number>
+// LAPACKFullMatrix<number>::()
+// {}
+
+
+// template <typename number>
+// LAPACKFullMatrix<number>::()
+// {}
+
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::vmult_add (
+  Vector<number>       &w,
+  const Vector<number> &v) const
+{
+  vmult(w, v, true);
+}
+
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::Tvmult_add (
+  Vector<number>       &w,
+  const Vector<number> &v) const
+{
+  Tvmult(w, v, true);
+}
+
+
+template LAPACKFullMatrix<double>;

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.