C++-Programmierung/ Eine Matrix-Bibliothek – mitrax/ simple operator.hpp

Aus Wikibooks
Zur Navigation springen Zur Suche springen
Nuvola-inspired-terminal.svg
  1 #ifndef _mitrax_mitrax_simple_operator_hpp_INCLUDED_
  2 #define _mitrax_mitrax_simple_operator_hpp_INCLUDED_
  3 /// \file simple_operator.hpp
  4 ///
  5 /// \brief Operationen, die auf einer Matrix ausgeführt werden
  6 ///
  7 /// Diese Datei stellt einige grundlegenden Operationen bereit, welche auf einer Matrix ausgeführt
  8 /// werden können. Der Nutzer von mitrax kann in einem eigenen Namensraum beliebige weitere
  9 /// Operationen schreiben. Wenn diese per <code>using</code>-Deklaration oder
 10 /// <code>using</code>-Direktive verfügbar gemacht werden, lassen sie sich syntaktisch genau so
 11 /// verwenden wie die hier definierten Funktionen. Ist eine Funktion sowohl als Erweiterung, als
 12 /// auch als native (hier definierte Funktion) vorhanden, so wird der nativen Implementierung,
 13 /// aufgrund der parameterabhängigen Namenssuchen, der Vorzug gegeben.
 14 
 15 #include "exception.hpp"
 16 
 17 namespace mitrax{
 18 
 19 /// \brief Elementweise Negation einer Matrix
 20 template < typename T, typename Container >
 21 inline matrix< T, Container > const operator-(matrix< T, Container > op){
 22     for(typename matrix< T, Container >::iterator i = op.begin(); i != op.end(); ++i){
 23         *i = -*i;
 24     }
 25     return op;
 26 }
 27 
 28 /// \brief Vorzeichen Plus
 29 template < typename T, typename Container >
 30 inline matrix< T, Container > const operator+(matrix< T, Container > const& op){
 31     return op;
 32 }
 33 
 34 
 35 /// \brief Additionszuweisung einer Matrix an eine andere
 36 template < typename T, typename Container >
 37 inline matrix< T, Container >&
 38 operator+=(matrix< T, Container >& lhs, matrix< T, Container > const& rhs){
 39     if(lhs.dimension() != rhs.dimension()){
 40         throw error::dimension_unequal("mitrax::operator+=<>()", lhs.dimension(), rhs.dimension());
 41     }
 42     typename matrix< T, Container >::iterator       target = lhs.begin();
 43     typename matrix< T, Container >::const_iterator add    = rhs.cbegin();
 44     for(;target != lhs.end(); ++target, ++add){
 45         *target += *add;
 46     }
 47     return lhs;
 48 }
 49 
 50 /// \brief Subtraktionszuweisung einer Matrix an eine andere
 51 template < typename T, typename Container >
 52 inline matrix< T, Container >&
 53 operator-=(matrix< T, Container >& lhs, matrix< T, Container > const& rhs){
 54     if(lhs.dimension() != rhs.dimension()){
 55         throw error::dimension_unequal("mitrax::operator-=<>()", lhs.dimension(), rhs.dimension());
 56     }
 57     typename matrix< T, Container >::iterator       target = lhs.begin();
 58     typename matrix< T, Container >::const_iterator sub    = rhs.cbegin();
 59     for(;target != lhs.end(); ++target, ++sub){
 60         *target -= *sub;
 61     }
 62     return lhs;
 63 }
 64 
 65 
 66 /// \brief Elementweise Addition zweier Matrizen
 67 template < typename T, typename Container >
 68 inline matrix< T, Container > const
 69 operator+(matrix< T, Container > lhs, matrix< T, Container > const& rhs){
 70     return lhs += rhs;
 71 }
 72 
 73 /// \brief Elementweise Subtraktion zweier Matrizen
 74 template < typename T, typename Container >
 75 inline matrix< T, Container > const
 76 operator-(matrix< T, Container > lhs, matrix< T, Container > const& rhs){
 77     return lhs -= rhs;
 78 }
 79 
 80 
 81 namespace detail{
 82 
 83 /// \brief Multiplikation aller Elemente einer Matrix mit einem Skalar (Helferfunktion)
 84 template < typename T, typename Container >
 85 inline matrix< T, Container >& scalar_multiplication(matrix< T, Container >& lhs, T const& rhs){
 86     typedef typename matrix< T, Container >::iterator iterator;
 87     for(iterator target = lhs.begin(); target != lhs.end(); ++target){
 88         *target *= rhs;
 89     }
 90     return lhs;
 91 }
 92 
 93 /// \brief Division aller Elemente einer Matrix durch einen Skalar (Helferfunktion)
 94 template < typename T, typename Container >
 95 inline matrix< T, Container >& scalar_division(matrix< T, Container >& lhs, T const& rhs){
 96     typedef typename matrix< T, Container >::iterator iterator;
 97     for(iterator target = lhs.begin(); target != lhs.end(); ++target){
 98         *target /= rhs;
 99     }
100     return lhs;
101 }
102 
103 }
104 
105 /// \brief Multiplikation aller Elemente einer Matrix mit einem Skalar
106 template < typename T, typename Container, typename SkalarType >
107 inline matrix< T, Container >& operator*=(matrix< T, Container >& lhs, SkalarType const& rhs){
108     return detail::scalar_multiplication< T, Container >(lhs, rhs);
109 }
110 
111 /// \brief Division aller Elemente einer Matrix durch einen Skalar
112 template < typename T, typename Container, typename SkalarType >
113 inline matrix< T, Container >& operator/=(matrix< T, Container >& lhs, SkalarType const& rhs){
114     return detail::scalar_division< T, Container >(lhs, rhs);
115 }
116 
117 /// \brief Matrixmultiplikation mit Zuweisung
118 template < typename T, typename Container >
119 inline matrix< T, Container >&
120 operator*=(matrix< T, Container >& lhs, matrix< T, Container > const& rhs){
121     typedef typename matrix< T, Container >::size_type                    size_type;
122     typedef typename matrix< T, Container >::row_const_proxy::iterator    row_iterator;
123     typedef typename matrix< T, Container >::column_const_proxy::iterator column_iterator;
124 
125     // Dimensionskompatibilität prüfen
126     if(lhs.columns() != rhs.rows()){
127         throw error::dimension_incompatible(
128             "mitrax::operator*() - Matrix multiplication requires lhs.columns() == rhs.rows()",
129             lhs.dimension(), rhs.dimension()
130         );
131     }
132 
133     matrix< T, Container > result(lhs.rows(), rhs.columns());
134     for(size_type column = size_type(); column < rhs.columns(); ++column){
135         for(size_type row = size_type(); row < lhs.rows(); ++row){
136             result[row][column] = lhs[row][0] * rhs[0][column];
137             for(size_type line = size_type()+1; line < lhs.columns(); ++line){
138                 result[row][column] += lhs[row][line] * rhs[line][column];
139             }
140         }
141     }
142 
143     return lhs = result;
144 }
145 
146 /// \brief Skalarmultiplikation
147 template < typename T, typename Container, typename SkalarType >
148 inline matrix< T, Container > const operator*(matrix< T, Container > lhs, SkalarType const& rhs){
149     return lhs *= rhs;
150 }
151 
152 /// \brief Skalarmultiplikation
153 template < typename T, typename Container, typename SkalarType >
154 matrix< T, Container > const operator*(SkalarType const& lhs, matrix< T, Container > rhs){
155     return rhs *= lhs;
156 }
157 
158 /// \brief Skalardivision
159 template < typename T, typename Container, typename SkalarType >
160 inline matrix< T, Container > const operator/(matrix< T, Container > lhs, SkalarType const& rhs){
161     return lhs /= rhs;
162 }
163 
164 /// \brief Matrixmultiplikation
165 template < typename T, typename Container >
166 inline matrix< T, Container > const
167 operator*(matrix< T, Container > lhs, matrix< T, Container > const& rhs){
168     return lhs *= rhs;
169 }
170 
171 namespace detail{
172 
173 /// \brief Funktor zum Vertauschen zweier Matrixproxys
174 template < typename Proxy >
175 class proxy_swap_functor{
176 public:
177     /// \brief Initialisierung mit <code>false</code>
178     proxy_swap_functor():odd_(false){}
179 
180     /// \brief Vertauscht zwei Proxys und wechselt intern seinen Wert
181     void operator()(Proxy& lhs, Proxy& rhs){
182         using std::swap;
183         swap(lhs, rhs);
184         odd_ = !odd_;
185     }
186 
187     /// \brief <code>true</code> für eine ungerade Anzahl von erfolgten Vertauschungen
188     operator bool(){return odd_;}
189 
190 private:
191     /// \brief <code>true</code> für eine ungerade Anzahl von erfolgten Vertauschungen
192     bool odd_;
193 };
194 
195 }
196 
197 
198 /// \brief Determinantenberechnung (nur quadratische Matrizen)
199 template < typename T, typename Container >
200 inline T const det(matrix< T, Container > const& quad_matrix){
201     typedef typename matrix< T, Container >::size_type                    size_type;
202     typedef typename matrix< T, Container >::row_proxy                    row_proxy;
203     typedef typename matrix< T, Container >::column_const_proxy::iterator column_iterator;
204 
205     if(quad_matrix.rows() != quad_matrix.columns()){
206         throw error::not_quadratic("mitrax::det<>()", quad_matrix.dimension());
207     }
208 
209     size_type size = quad_matrix.rows();
210     matrix< T, Container > original(quad_matrix);
211 
212     // Liste von Zeilen-Proxys erstellen
213     std::vector< row_proxy > lines;
214     lines.reserve(size);
215     for(size_type line = size_type(); line < size; ++line){
216         lines.push_back(original.row(line));
217     }
218 
219     detail::proxy_swap_functor< row_proxy > proxy_swap;
220 
221     // Obere Dreiecksmatrix bilden
222     for(size_type line = size_type(); line < size; ++line){
223         // Zeile mit einer folgenden tauschen, falls nötig
224         if(lines[line][line] == identity< T >::additive()){
225             bool status = proxy_swap;
226             for(size_type row = line+1; row < size; ++row){
227                 if(lines[row][line] == identity< T >::additive()) continue;
228                 proxy_swap(lines[line], lines[row]);
229                 break;
230             }
231             // Kein Tausch durchgeführt? ja -> det == 0
232             if(status == proxy_swap){
233                 return identity< T >::additive();
234             }
235         }
236         // Alle folgenden Zeilen in dieser Spalte 0 machen
237         // Die Berechnung der Spalte selbst kann man sich sparen, deshalb line+1
238         for(size_type row = line+1; row < size; ++row){
239             T factor = (lines[row][line] / lines[line][line]);
240             for(size_type column = line+1; column < size; ++column){
241                 lines[row][column] -= factor * lines[line][column];
242             }
243         }
244     }
245 
246     // Diagonalelemente Multiplizieren
247     T result(lines[size_type()][size_type()]);
248     for(size_type line = size_type()+1; line < size; ++line){
249         result *= lines[line][line];
250     }
251 
252     // Negieren, falls Anzahl Vertauschungen ungerade
253     if(proxy_swap) result = -result;
254 
255     return result;
256 }
257 
258 /// \brief Inverse Matrix (nur reguläre quadratische Matrizen)
259 template < typename T, typename Container >
260 inline matrix< T, Container > const inverse(matrix< T, Container > const& quad_matrix){
261     typedef typename matrix< T, Container >::size_type                    size_type;
262     typedef typename matrix< T, Container >::row_proxy                    row_proxy;
263     typedef typename matrix< T, Container >::column_const_proxy::iterator column_iterator;
264 
265     if(quad_matrix.rows() != quad_matrix.columns()){
266         throw error::not_quadratic("mitrax::inverse<>()", quad_matrix.dimension());
267     }
268 
269     size_type size = quad_matrix.rows();
270     matrix< T, Container > original(quad_matrix);
271 
272     // Einheitsmatrix erzeugen
273     matrix< T, Container > result(quad_matrix.dimension());
274     {
275         typename matrix< T, Container >::iterator iter = result.begin();
276         for(size_type i = size_type(); i < size; ++i){
277             *iter = identity< T >::multiplicative();
278             iter += size + 1;
279         }
280     }
281 
282     // Liste von Zeilen-Proxys erstellen
283     std::vector< row_proxy > org_lines;
284     org_lines.reserve(size);
285     for(size_type line = size_type(); line < size; ++line){
286         org_lines.push_back(original.row(line));
287     }
288 
289     detail::proxy_swap_functor< row_proxy > proxy_swap;
290 
291     // Obere Dreiecksmatrix bilden
292     for(size_type line = size_type(); line < size; ++line){
293         // Zeile mit einer folgenden tauschen, falls nötig
294         if(org_lines[line][line] == identity< T >::additive()){
295             bool status = proxy_swap;
296             for(size_type row = line+1; row < size; ++row){
297                 if(org_lines[row][line] == identity< T >::additive()) continue;
298                 proxy_swap(org_lines[line], org_lines[row]);
299                 element_swap(result[line], result[row]);
300                 break;
301             }
302             // Kein Tausch durchgeführt? ja -> det == 0 -> nicht invertierbar
303             if(status == proxy_swap){
304                 throw error::singular_matrix("mitrax::inverse<>()");
305             }
306         }
307         // Alle folgenden Zeilen in dieser Spalte 0 machen
308         // (Elemente dieser Spalte müssen nicht berechnet werden)
309         for(size_type row = line+1; row < size; ++row){
310             T factor = (org_lines[row][line] / org_lines[line][line]);
311             for(size_type column = line+1; column < size; ++column){
312                 org_lines[row][column] -= factor * org_lines[line][column];
313             }
314             for(size_type column = size_type(); column < size; ++column){
315                 result[row][column] -= factor * result[line][column];
316             }
317         }
318         // Zeile normieren (Diagonalelemente müssen nicht berechnet werden)
319         for(size_type column = line+1; column < size; ++column){
320             org_lines[line][column] /= org_lines[line][line];
321         }
322         for(size_type column = size_type(); column < size; ++column){
323             result[line][column] /= org_lines[line][line];
324         }
325     }
326 
327     // Einheitsmatrix bilden (Elemente in org_lines müssen nicht berechnet werden)
328     for(size_type line = size-1; line > size_type(); --line){
329         for(size_type row = size_type(); row < line; ++row){
330             T factor(org_lines[row][line]);
331             for(size_type column = size_type(); column < size; ++column){
332                 result[row][column] -= factor * result[line][column];
333             }
334         }
335     }
336 
337     return result;
338 }
339 
340 
341 /// \brief <code>true</code> für gleiche Dimension und Elemente
342 template < typename T, typename Container >
343 inline bool operator==(matrix< T, Container > const& lhs, matrix< T, Container > const& rhs){
344     typedef typename matrix< T, Container >::const_iterator const_iterator;
345 
346     if(lhs.dimension() != rhs.dimension()){
347         return false;
348     }
349 
350     for(const_iterator lp = lhs.cbegin(), rp = rhs.cbegin(); lp != lhs.end(); ++lp, ++rp){
351         if(*lp != *rp) return false;
352     }
353 
354     return true;
355 }
356 
357 /// \brief <code>true</code> für ungleiche Dimension oder ungleiche Elemente
358 template < typename T, typename Container >
359 inline bool operator!=(dimension< T > const& lhs, dimension< T > const& rhs){
360     return !(lhs == rhs);
361 }
362 
363 
364 }
365 
366 #endif