dune-localfunctions  2.8.0
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  using std::abs;
90  assert( rows() == cols() );
91  std::vector<unsigned int> p(rows());
92  for (unsigned int j=0; j<rows(); ++j)
93  p[j] = j;
94  for (unsigned int j=0; j<rows(); ++j)
95  {
96  // pivot search
97  unsigned int r = j;
98  Field max = abs( (*this)(j,j) );
99  for (unsigned int i=j+1; i<rows(); ++i)
100  {
101  if ( abs( (*this)(i,j) ) > max )
102  {
103  max = abs( (*this)(i,j) );
104  r = i;
105  }
106  }
107  if (max == Zero<Field>())
108  return false;
109  // row swap
110  if (r > j)
111  {
112  for (unsigned int k=0; k<cols(); ++k)
113  std::swap( (*this)(j,k), (*this)(r,k) );
114  std::swap( p[j], p[r] );
115  }
116  // transformation
117  Field hr = Unity<Field>()/(*this)(j,j);
118  for (unsigned int i=0; i<rows(); ++i)
119  (*this)(i,j) *= hr;
120  (*this)(j,j) = hr;
121  for (unsigned int k=0; k<cols(); ++k)
122  {
123  if (k==j) continue;
124  for (unsigned int i=0; i<rows(); ++i)
125  {
126  if (i==j) continue;
127  (*this)(i,k) -= (*this)(i,j)*(*this)(j,k);
128  }
129  (*this)(j,k) *= -hr;
130  }
131  }
132  // column exchange
133  Row hv(rows());
134  for (unsigned int i=0; i<rows(); ++i)
135  {
136  for (unsigned int k=0; k<rows(); ++k)
137  hv[ p[k] ] = (*this)(i,k);
138  for (unsigned int k=0; k<rows(); ++k)
139  (*this)(i,k) = hv[k];
140  }
141  return true;
142  }
143 
144  private:
145  RealMatrix matrix_;
146  unsigned int cols_,rows_;
147  };
148 
149  template< class Field >
150  inline std::ostream &operator<<(std::ostream &out, const LFEMatrix<Field> &mat)
151  {
152  for (unsigned int r=0; r<mat.rows(); ++r)
153  {
154  out << field_cast<double>(mat(r,0));
155  for (unsigned int c=1; c<mat.cols(); ++c)
156  {
157  out << " , " << field_cast<double>(mat(r,c));
158  }
159  out << std::endl;
160  }
161  return out;
162  }
163 }
164 
165 #endif // #ifndef DUNE_LOCALFUNCTIONS_UTILITY_LFEMATRIX_HH
Definition: bdfmcube.hh:16
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition: field.hh:157
std::ostream & operator<<(std::ostream &out, const LFEMatrix< Field > &mat)
Definition: lfematrix.hh:150
A class representing the unit of a given Field.
Definition: field.hh:28
A class representing the zero of a given Field.
Definition: field.hh:77
Definition: lfematrix.hh:16
const Field & operator()(const unsigned int row, const unsigned int col) const
Definition: lfematrix.hh:42
unsigned int cols() const
Definition: lfematrix.hh:61
void resize(const unsigned int rows, const unsigned int cols)
Definition: lfematrix.hh:78
void row(const unsigned int row, Vector &vec) const
Definition: lfematrix.hh:35
unsigned int rows() const
Definition: lfematrix.hh:56
const Field * rowPtr(const unsigned int row) const
Definition: lfematrix.hh:66
bool invert()
Definition: lfematrix.hh:87
Field * rowPtr(const unsigned int row)
Definition: lfematrix.hh:72
F Field
Definition: lfematrix.hh:22