From: guido Date: Wed, 23 Mar 2005 03:49:05 +0000 (+0000) Subject: starting LAPACK matrix X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1e5ffa92b6b6a5251c4b02911c370215ccc71996;p=dealii-svn.git starting LAPACK matrix git-svn-id: https://svn.dealii.org/trunk@10211 0785d39b-7218-0410-832d-ea1e28bc413d --- 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 index 0000000000..472d15aa3e --- /dev/null +++ b/deal.II/lac/include/lac/lapack_full_matrix.h @@ -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 +#include + +// forward declarations +template class Vector; +template class FullMatrix; + +/*! @addtogroup Matrix1 + *@{ + */ + + +/** + * A variant of FullMatrix using LAPACK functions whereever possible. + * + * @author Guido Kanschat, 2005 + */ +template +class LAPACKFullMatrix : public TransposeTable +{ + public: + /** + * Constructor. Initialize the + * matrix as a square matrix with + * dimension n. + * + * In order to avoid the implicit + * conversion of integers and + * other types to a matrix, this + * constructor is declared + * explicit. + * + * 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 explicit, + * since that prevents the use of + * this class in containers, such + * as std::vector. The + * responsibility to check + * performance of programs must + * therefore remain with the + * user of this class. + */ + LAPACKFullMatrix (const LAPACKFullMatrix&); + + /** + * Assignment operator. + */ + LAPACKFullMatrix & + operator = (const LAPACKFullMatrix&); + + /** + * Assignment operator for a + * regular FullMatrix. + */ + template + LAPACKFullMatrix & + operator = (const FullMatrix&); + + /** + * This operator assigns a scalar + * to a matrix. To avoid + * confusion with constructors, + * zero is the only value allowed + * for d + */ + LAPACKFullMatrix & + 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 + void copy_from (const MATRIX&); + + /** + * Fill rectangular block. + * + * A rectangular block of the + * matrix src is copied into + * this. The upper left + * corner of the block being + * copied is + * (src_offset_i,src_offset_j). + * The upper left corner of the + * copied block is + * (dst_offset_i,dst_offset_j). + * The size of the rectangular + * block being copied is the + * maximum size possible, + * determined either by the size + * of this or src. + */ + template + 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 + * adding determines, whether the + * result is stored in w or added + * to w. + * + * if (adding) + * w += A*v + * + * if (!adding) + * w = A*v + * + * Source and destination must + * not be the same vector. + */ + void vmult (Vector &w, + const Vector &v, + const bool adding=false) const; + /** + * Adding Matrix-vector-multiplication. + * w += A*v + * + * Source and destination must + * not be the same vector. + */ + void vmult_add (Vector &w, + const Vector &v) const; + + /** + * Transpose + * matrix-vector-multiplication. + * + * The optional parameter + * adding determines, whether the + * result is stored in w or added + * to w. + * + * if (adding) + * w += AT*v + * + * if (!adding) + * w = AT*v + * + * + * Source and destination must + * not be the same vector. + */ + void Tvmult (Vector &w, + const Vector &v, + const bool adding=false) const; + + /** + * Adding transpose + * matrix-vector-multiplication. + * w += AT*v + * + * Source and destination must + * not be the same vector. + */ + void Tvmult_add (Vector &w, + const Vector &v) const; + +}; + + +template +template +void +LAPACKFullMatrix::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 +template +void +LAPACKFullMatrix::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;icolumn() < 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 index 0000000000..45b7792c7c --- /dev/null +++ b/deal.II/lac/include/lac/lapack_support.h @@ -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 deal.II 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 index 0000000000..26c20cf0ef --- /dev/null +++ b/deal.II/lac/include/lac/lapack_templates.h @@ -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 diff --git a/deal.II/lac/include/lac/vector.h b/deal.II/lac/include/lac/vector.h index 9bd143ed07..98f7790064 100644 --- a/deal.II/lac/include/lac/vector.h +++ b/deal.II/lac/include/lac/vector.h @@ -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 class LAPACKFullMatrix; /*! @addtogroup Vectors *@{ @@ -748,11 +748,16 @@ class Vector */ Number *val; - /** + /* * Make all other vector types * friends. */ - template friend class Vector; + template friend class Vector; + /* + * LAPACK matrices need access to + * the data. + */ + friend class LAPACKFullMatrix; }; /*@}*/ diff --git a/deal.II/lac/source/lapack_full_matrix.cc b/deal.II/lac/source/lapack_full_matrix.cc new file mode 100644 index 0000000000..830fb889f2 --- /dev/null +++ b/deal.II/lac/source/lapack_full_matrix.cc @@ -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 +#include +#include +#include + +using namespace LAPACKSupport; + +template +LAPACKFullMatrix::LAPACKFullMatrix(const unsigned int n) + : + TransposeTable (n,n) +{} + + +template +LAPACKFullMatrix::LAPACKFullMatrix( + const unsigned int m, + const unsigned int n) + : + TransposeTable (m,n) +{} + + +template +LAPACKFullMatrix::LAPACKFullMatrix(const LAPACKFullMatrix &M) + : + TransposeTable (M) +{} + + +template +LAPACKFullMatrix & +LAPACKFullMatrix::operator = (const LAPACKFullMatrix& M) +{ + TransposeTable::operator=(M); + return *this; +} + + +template +template +LAPACKFullMatrix & +LAPACKFullMatrix::operator = (const FullMatrix& M) +{ + Assert (n() == M.m(), ExcDimensionMismatch(n(), M.m())); + Assert (n() == M.n(), ExcDimensionMismatch(n(), M.n())); + for (unsigned int i=0;i +LAPACKFullMatrix & +LAPACKFullMatrix::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 +void +LAPACKFullMatrix::vmult ( + Vector &w, + const Vector &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 +LAPACKFullMatrix::Tvmult ( + Vector &w, + const Vector &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 +void +LAPACKFullMatrix::vmult ( + Vector &, + const Vector &, + const bool ) const +{ + Assert(false, ExcNeedsBLAS()); +} + + +template +void +LAPACKFullMatrix::Tvmult ( + Vector &, + const Vector &, + const bool ) const +{ + Assert(false, ExcNeedsBLAS()); +} + + +#endif + +// template +// LAPACKFullMatrix::() +// {} + + +// template +// LAPACKFullMatrix::() +// {} + + +// template +// LAPACKFullMatrix::() +// {} + + +template +void +LAPACKFullMatrix::vmult_add ( + Vector &w, + const Vector &v) const +{ + vmult(w, v, true); +} + + +template +void +LAPACKFullMatrix::Tvmult_add ( + Vector &w, + const Vector &v) const +{ + Tvmult(w, v, true); +} + + +template LAPACKFullMatrix;