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