include/deal.II/lac/block_matrix_array.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: block_matrix_array.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2012 by the deal.II authors
00005 //
00006 //    This file is subject to QPL and may not be  distributed
00007 //    without copyright and license information. Please refer
00008 //    to the file deal.II/doc/license.html for the  text  and
00009 //    further information on this license.
00010 //
00011 //---------------------------------------------------------------------------
00012 #ifndef __deal2__block_matrix_array_h
00013 #define __deal2__block_matrix_array_h
00014 
00015 #include <deal.II/base/config.h>
00016 #include <deal.II/base/subscriptor.h>
00017 #include <deal.II/base/table.h>
00018 
00019 #include <deal.II/lac/pointer_matrix.h>
00020 #include <deal.II/lac/vector_memory.h>
00021 
00022 #include <vector>
00023 #include <map>
00024 #include <string>
00025 #include <memory>
00026 #include <sstream>
00027 
00028 DEAL_II_NAMESPACE_OPEN
00029 
00030 template <typename> class BlockVector;
00031 template <typename> class Vector;
00032 
00122 template <typename number = double>
00123 class BlockMatrixArray : public Subscriptor
00124 {
00125   public:
00132     BlockMatrixArray ();
00133 
00138     BlockMatrixArray (const unsigned int n_block_rows,
00139                       const unsigned int n_block_cols);
00140 
00151     void initialize (const unsigned int n_block_rows,
00152                      const unsigned int n_block_cols);
00153 
00161     BlockMatrixArray (const unsigned int n_block_rows,
00162                       const unsigned int n_block_cols,
00163                       VectorMemory<Vector<number> >& mem);
00164 
00175     void initialize (const unsigned int n_block_rows,
00176                      const unsigned int n_block_cols,
00177                      VectorMemory<Vector<number> >& mem);
00178 
00183     void reinit (const unsigned int n_block_rows,
00184                  const unsigned int n_block_cols);
00185 
00205     template <class MATRIX>
00206     void enter (const MATRIX      &matrix,
00207                 const unsigned int row,
00208                 const unsigned int col,
00209                 const double       prefix = 1.,
00210                 const bool         transpose = false);
00211 
00221     template <class MATRIX>
00222     void enter_aux (VectorMemory<Vector<number> >& mem,
00223                     const MATRIX      &matrix,
00224                     const unsigned int row,
00225                     const unsigned int col,
00226                     const double       prefix = 1.,
00227                     const bool         transpose = false);
00228 
00229 
00234     void clear();
00235 
00240     unsigned int n_block_rows () const;
00241 
00246     unsigned int n_block_cols () const;
00247 
00251     void vmult (BlockVector<number>& dst,
00252                 const BlockVector<number>& src) const;
00253 
00258     void vmult_add (BlockVector<number>& dst,
00259                     const BlockVector<number>& src) const;
00260 
00265     void Tvmult (BlockVector<number>& dst,
00266                  const BlockVector<number>& src) const;
00267 
00273     void Tvmult_add (BlockVector<number>& dst,
00274                      const BlockVector<number>& src) const;
00275 
00281     number matrix_scalar_product (const BlockVector<number>& u,
00282                                   const BlockVector<number>& v) const;
00283 
00290     number matrix_norm_square (const BlockVector<number>& u) const;
00291 
00309     template <class STREAM>
00310     void print_latex (STREAM& out) const;
00311 
00312   protected:
00325     class Entry
00326     {
00327       public:
00335         template<class MATRIX>
00336         Entry (const MATRIX& matrix,
00337                unsigned row, unsigned int col,
00338                double prefix, bool transpose);
00339 
00353         Entry(const Entry&);
00354 
00361         ~Entry();
00362 
00367         unsigned int row;
00368 
00373         unsigned int col;
00374 
00379         double prefix;
00380 
00386         bool transpose;
00387 
00391         PointerMatrixBase<Vector<number> >* matrix;
00392     };
00393 
00398     std::vector<Entry> entries;
00399 
00400   private:
00404     unsigned int block_rows;
00408     unsigned int block_cols;
00409 };
00410 
00471 template <typename number = double>
00472 class BlockTrianglePrecondition
00473   : private BlockMatrixArray<number>
00474 {
00475   public:
00482     BlockTrianglePrecondition ();
00483 
00490     BlockTrianglePrecondition (unsigned int n_blocks);
00491 
00501     BlockTrianglePrecondition (unsigned int n_block_rows,
00502                                VectorMemory<Vector<number> >& mem,
00503                                bool backward = false);
00504 
00515     void initialize (const unsigned int n_block_rows,
00516                      VectorMemory<Vector<number> >& mem,
00517                      bool backward = false);
00518 
00523     void reinit(const unsigned int n_block_rows);
00524 
00525 
00533     template <class MATRIX>
00534     void enter (const MATRIX      &matrix,
00535                 const unsigned int row,
00536                 const unsigned int col,
00537                 const double       prefix = 1.,
00538                 const bool         transpose = false);
00539 
00551     template <class MATRIX>
00552     void enter_aux (VectorMemory<Vector<double> >& mem,
00553                     const MATRIX      &matrix,
00554                     const unsigned int row,
00555                     const unsigned int col,
00556                     const double       prefix = 1.,
00557                     const bool         transpose = false);
00558 
00562     void vmult (BlockVector<number>& dst,
00563                 const BlockVector<number>& src) const;
00564 
00569     void vmult_add (BlockVector<number>& dst,
00570                     const BlockVector<number>& src) const;
00571 
00575     void Tvmult (BlockVector<number>& dst,
00576                  const BlockVector<number>& src) const;
00577 
00582     void Tvmult_add (BlockVector<number>& dst,
00583                      const BlockVector<number>& src) const;
00584 
00588     using BlockMatrixArray<number>::print_latex;
00589 
00593     using BlockMatrixArray<number>::n_block_rows;
00594 
00598     using BlockMatrixArray<number>::n_block_cols;
00599     using BlockMatrixArray<number>::clear;
00600     using BlockMatrixArray<number>::Subscriptor::subscribe;
00601     using BlockMatrixArray<number>::Subscriptor::unsubscribe;
00602 
00613     DeclException1(ExcNoDiagonal,
00614                    unsigned int,
00615                    << "No diagonal entry was added for block " << arg1);
00616 
00624     DeclException1(ExcMultipleDiagonal,
00625                    unsigned int,
00626                    << "Inverse diagonal entries may not be added in block "
00627                    << arg1);
00629   private:
00636     void do_row (BlockVector<number>& dst,
00637                  unsigned int row_num) const;
00638 
00642     bool backward;
00643 };
00644 
00645 
00646 #ifndef DOXYGEN
00647 //---------------------------------------------------------------------------
00648 
00649 template <typename number>
00650 template <class MATRIX>
00651 inline
00652 BlockMatrixArray<number>::Entry::Entry (
00653   const MATRIX& m,
00654   unsigned int row,
00655   unsigned int col,
00656   double prefix,
00657   bool transpose)
00658                 :
00659                 row (row),
00660                 col (col),
00661                 prefix (prefix),
00662                 transpose (transpose),
00663                 matrix (new_pointer_matrix_base(m, Vector<number>(), typeid(*this).name()))
00664 {}
00665 
00666 
00667 
00668 template <typename number>
00669 template <class MATRIX>
00670 inline
00671 void
00672 BlockMatrixArray<number>::enter (
00673   const MATRIX& matrix,
00674   unsigned int row,
00675   unsigned int col,
00676   double prefix,
00677   bool transpose)
00678 {
00679   Assert(row<n_block_rows(), ExcIndexRange(row, 0, n_block_rows()));
00680   Assert(col<n_block_cols(), ExcIndexRange(col, 0, n_block_cols()));
00681   entries.push_back(Entry(matrix, row, col, prefix, transpose));
00682 }
00683 
00684 
00685 template <typename number>
00686 template <class MATRIX>
00687 inline
00688 void
00689 BlockMatrixArray<number>::enter_aux (
00690   VectorMemory<Vector<number> >& mem,
00691   const MATRIX& matrix,
00692   unsigned int row,
00693   unsigned int col,
00694   double prefix,
00695   bool transpose)
00696 {
00697   Assert(row<n_block_rows(), ExcIndexRange(row, 0, n_block_rows()));
00698   Assert(col<n_block_cols(), ExcIndexRange(col, 0, n_block_cols()));
00699   entries.push_back(Entry(matrix, row, col, prefix, transpose, mem));
00700 }
00701 
00702 
00703 template<typename number>
00704 struct BlockMatrixArrayPointerMatrixLess
00705 {
00706     bool operator () (const PointerMatrixBase<Vector<number> >*a,
00707                       const PointerMatrixBase<Vector<number> >* b) const
00708       {
00709         return *a < *b;
00710       }
00711 };
00712 
00713 
00714 template <typename number>
00715 template <class STREAM>
00716 inline
00717 void
00718 BlockMatrixArray<number>::print_latex (STREAM& out) const
00719 {
00720   out << "\\begin{array}{"
00721       << std::string(n_block_cols(), 'c')
00722       << "}" << std::endl;
00723 
00724   Table<2,std::string> array(n_block_rows(), n_block_cols());
00725 
00726   typedef std::map<const PointerMatrixBase<Vector<number> >*,
00727     std::string, BlockMatrixArrayPointerMatrixLess<number> > NameMap;
00728   NameMap matrix_names;
00729 
00730   typename std::vector<Entry>::const_iterator m = entries.begin();
00731   typename std::vector<Entry>::const_iterator end = entries.end();
00732 
00733   unsigned int matrix_number = 0;
00734   for (; m != end ; ++m)
00735     {
00736       if (matrix_names.find(m->matrix) == matrix_names.end())
00737         {
00738           std::pair<typename NameMap::iterator, bool> x =
00739             matrix_names.insert(
00740               std::pair<const PointerMatrixBase<Vector<number> >*, std::string> (m->matrix,
00741                                                      std::string("M")));
00742           std::ostringstream stream;
00743           stream << matrix_number++;
00744 
00745           x.first->second += stream.str();
00746         }
00747 
00748       std::ostringstream stream;
00749 
00750       if (array(m->row, m->col) != "" && m->prefix >= 0)
00751         stream << "+";
00752       if (m->prefix != 1.)
00753         stream << m->prefix << 'x';
00754       stream << matrix_names.find(m->matrix)->second;
00755 //      stream << '(' << m->matrix << ')';
00756       if (m->transpose)
00757         stream << "^T";
00758 
00759       array(m->row, m->col) += stream.str();
00760     }
00761   for (unsigned int i=0;i<n_block_rows();++i)
00762     for (unsigned int j=0;j<n_block_cols();++j)
00763       {
00764         out << '\t' << array(i,j);
00765         if (j==n_block_cols()-1)
00766           out << "\\\\" << std::endl;
00767         else
00768           out << " &";
00769       }
00770   out << "\\end{array}" << std::endl;
00771 }
00772 
00773 template <typename number>
00774 template <class MATRIX>
00775 inline
00776 void
00777 BlockTrianglePrecondition<number>::enter (const MATRIX& matrix,
00778                                           unsigned row, unsigned int col,
00779                                           double prefix, bool transpose)
00780 {
00781   BlockMatrixArray<number>::enter(matrix, row, col, prefix, transpose);
00782 }
00783 
00784 
00785 
00786 template <typename number>
00787 template <class MATRIX>
00788 inline
00789 void
00790 BlockTrianglePrecondition<number>::enter_aux (
00791   VectorMemory<Vector<double> >& mem,
00792   const MATRIX& matrix,
00793   unsigned int row,
00794   unsigned int col,
00795   double prefix,
00796   bool transpose)
00797 {
00798   BlockMatrixArray<number>::enter_aux(mem, matrix, row, col, prefix, transpose);
00799 }
00800 
00801 
00802 
00803 #endif // DOXYGEN
00804 
00805 DEAL_II_NAMESPACE_CLOSE
00806 
00807 #endif
00808 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Thu May 17 2012 20:04:39 by doxygen 1.7.3