dune-localfunctions  2.5.1
lfematrix.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_LOCALFUNCTIONS_UTILITY_LFEMATRIX_HH
4 #define DUNE_LOCALFUNCTIONS_UTILITY_LFEMATRIX_HH
5 
6 #include <cassert>
7 #include <vector>
8 
9 #include "field.hh"
10 
11 namespace Dune
12 {
13 
14  template< class F >
15  class LFEMatrix
16  {
17  typedef LFEMatrix< F > This;
18  typedef std::vector< F > Row;
19  typedef std::vector<Row> RealMatrix;
20 
21  public:
22  typedef F Field;
23 
24  operator const RealMatrix & () const
25  {
26  return matrix_;
27  }
28 
29  operator RealMatrix & ()
30  {
31  return matrix_;
32  }
33 
34  template <class Vector>
35  void row( const unsigned int row, Vector &vec ) const
36  {
37  assert(row<rows());
38  for (int i=0; i<cols(); ++i)
39  field_cast(matrix_[row][i], vec[i]);
40  }
41 
42  const Field &operator() ( const unsigned int row, const unsigned int col ) const
43  {
44  assert(row<rows());
45  assert(col<cols());
46  return matrix_[ row ][ col ];
47  }
48 
49  Field &operator() ( const unsigned int row, const unsigned int col )
50  {
51  assert(row<rows());
52  assert(col<cols());
53  return matrix_[ row ][ col ];
54  }
55 
56  unsigned int rows () const
57  {
58  return rows_;
59  }
60 
61  unsigned int cols () const
62  {
63  return cols_;
64  }
65 
66  const Field *rowPtr ( const unsigned int row ) const
67  {
68  assert(row<rows());
69  return &(matrix_[row][0]);
70  }
71 
72  Field *rowPtr ( const unsigned int row )
73  {
74  assert(row<rows());
75  return &(matrix_[row][0]);
76  }
77 
78  void resize ( const unsigned int rows, const unsigned int cols )
79  {
80  matrix_.resize(rows);
81  for (unsigned int i=0; i<rows; ++i)
82  matrix_[i].resize(cols);
83  rows_ = rows;
84  cols_ = cols;
85  }
86 
87  bool invert ()
88  {
89  assert( rows() == cols() );
90  std::vector<unsigned int> p(rows());
91  for (unsigned int j=0; j<rows(); ++j)
92  p[j] = j;
93  for (unsigned int j=0; j<rows(); ++j)
94  {
95  // pivot search
96  unsigned int r = j;
97  Field max = std::abs( (*this)(j,j) );
98  for (unsigned int i=j+1; i<rows(); ++i)
99  {
100  if ( std::abs( (*this)(i,j) ) > max )
101  {
102  max = std::abs( (*this)(i,j) );
103  r = i;
104  }
105  }
106  if (max == Zero<Field>())
107  return false;
108  // row swap
109  if (r > j)
110  {
111  for (unsigned int k=0; k<cols(); ++k)
112  std::swap( (*this)(j,k), (*this)(r,k) );
113  std::swap( p[j], p[r] );
114  }
115  // transformation
116  Field hr = Unity<Field>()/(*this)(j,j);
117  for (unsigned int i=0; i<rows(); ++i)
118  (*this)(i,j) *= hr;
119  (*this)(j,j) = hr;
120  for (unsigned int k=0; k<cols(); ++k)
121  {
122  if (k==j) continue;
123  for (unsigned int i=0; i<rows(); ++i)
124  {
125  if (i==j) continue;
126  (*this)(i,k) -= (*this)(i,j)*(*this)(j,k);
127  }
128  (*this)(j,k) *= -hr;
129  }
130  }
131  // column exchange
132  Row hv(rows());
133  for (unsigned int i=0; i<rows(); ++i)
134  {
135  for (unsigned int k=0; k<rows(); ++k)
136  hv[ p[k] ] = (*this)(i,k);
137  for (unsigned int k=0; k<rows(); ++k)
138  (*this)(i,k) = hv[k];
139  }
140  return true;
141  }
142 
143  private:
144  RealMatrix matrix_;
145  unsigned int cols_,rows_;
146  };
147 
148  template< class Field >
149  inline std::ostream &operator<<(std::ostream &out, const LFEMatrix<Field> &mat)
150  {
151  for (unsigned int r=0; r<mat.rows(); ++r)
152  {
153  out << field_cast<double>(mat(r,0));
154  for (unsigned int c=1; c<mat.cols(); ++c)
155  {
156  out << " , " << field_cast<double>(mat(r,c));
157  }
158  out << std::endl;
159  }
160  return out;
161  }
162 }
163 
164 #endif // #ifndef DUNE_LOCALFUNCTIONS_UTILITY_LFEMATRIX_HH
A class representing the zero of a given Field.
Definition: field.hh:76
const Field * rowPtr(const unsigned int row) const
Definition: lfematrix.hh:66
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:15
unsigned int rows() const
Definition: lfematrix.hh:56
Field * rowPtr(const unsigned int row)
Definition: lfematrix.hh:72
Definition: lfematrix.hh:15
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition: field.hh:157
const Field & operator()(const unsigned int row, const unsigned int col) const
Definition: lfematrix.hh:42
A class representing the unit of a given Field.
Definition: field.hh:27
void resize(const unsigned int rows, const unsigned int cols)
Definition: lfematrix.hh:78
F Field
Definition: lfematrix.hh:22
bool invert()
Definition: lfematrix.hh:87
void row(const unsigned int row, Vector &vec) const
Definition: lfematrix.hh:35
unsigned int cols() const
Definition: lfematrix.hh:61