Dune Core Modules (2.8.0)

lfematrix.hh
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
11namespace 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
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
Dune namespace.
Definition: alignedallocator.hh:11
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition: field.hh:157
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 22, 23:33, 2024)