--- /dev/null
+//-------------------------------------------------------------------
+// $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
--- /dev/null
+//-------------------------------------------------------------------
+// $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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
-//---------------------------- 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
}
#endif
-
+template<typename number> class LAPACKFullMatrix;
/*! @addtogroup Vectors
*@{
*/
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>;
};
/*@}*/
--- /dev/null
+//---------------------------------------------------------------------------
+// $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>;