All Classes Namespaces Functions Typedefs Enumerations Pages
linear_algebra.h
1 /* Copyright (C) 2003-2005 Martin Förg, Rober Huber
2 
3  * This library is free software; you can redistribute it and/or
4  * modify it under the terms of the GNU Lesser General Public
5  * License as published by the Free Software Foundation; either
6  * version 2.1 of the License, or (at your option) any later version.
7  *
8  * This library is distributed in the hope that it will be useful,
9  * but WITHOUT ANY WARRANTY; without even the implied warranty of
10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11  * Lesser General Public License for more details.
12  *
13  * You should have received a copy of the GNU Lesser General Public
14  * License along with this library; if not, write to the Free Software
15  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
16  *
17  * Contact:
18  * martin.o.foerg@googlemail.com
19  *
20  */
21 
22 #ifndef linear_algebra_h
23 #define linear_algebra_h
24 
25 #include <math.h>
26 
27 #include "square_matrix.h"
28 #include "vector.h"
29 #include "row_vector.h"
30 #include "fixed_square_matrix.h"
31 #include "fixed_vector.h"
32 #include "fixed_row_vector.h"
33 #include "fixed_symmetric_matrix.h"
34 #include "var_square_matrix.h"
35 #include "var_symmetric_matrix.h"
36 #include "var_vector.h"
37 #include "var_row_vector.h"
38 #include "var_fixed_general_matrix.h"
39 #include "fixed_var_general_matrix.h"
40 #include <cmath>
41 
42 namespace fmatvec {
43 
45 
51  template <class AT>
53 
54 #ifndef FMATVEC_NO_SIZE_CHECK
55  assert(y.size() == x.size());
56 #endif
57 
58  Vector<Ref, AT> z(x.size(), NONINIT);
59 
60  if (y.transposed()) {
61  if (x.transposed())
62  for (int i = 0; i < x.size(); i++)
63  z.er(i) = x.et(i) + y.et(i);
64  else
65  for (int i = 0; i < x.size(); i++)
66  z.er(i) = x.er(i) + y.et(i);
67  }
68  else {
69  if (x.transposed())
70  for (int i = 0; i < x.size(); i++)
71  z.er(i) = x.et(i) + y.er(i);
72  else
73  for (int i = 0; i < x.size(); i++)
74  z.er(i) = x.er(i) + y.er(i);
75  }
76 
77  return z;
78  }
79 
80  // Vector-Vector
81  template <class Row1, class Row2, class Row3, class AT>
82  inline void add(const Vector<Row1, AT> &a1, const Vector<Row2, AT> &a2, Vector<Row3, AT> &a3) {
83 #ifndef FMATVEC_NO_SIZE_CHECK
84  assert(a1.size() == a2.size());
85 #endif
86  for (int i = 0; i < a1.size(); i++)
87  a3.e(i) = a1.e(i) + a2.e(i);
88  }
89 
90  template <class Row1, class Row2, class AT>
91  inline void add(Vector<Row1, AT> &a1, const Vector<Row2, AT> &a2) {
92 #ifndef FMATVEC_NO_SIZE_CHECK
93  assert(a1.size() == a2.size());
94 #endif
95  for (int i = 0; i < a1.size(); i++)
96  a1.e(i) += a2.e(i);
97  }
98 
99  template <class Row1, class Row2, class Row3, class AT>
100  inline void sub(const Vector<Row1, AT> &a1, const Vector<Row2, AT> &a2, Vector<Row3, AT> &a3) {
101 #ifndef FMATVEC_NO_SIZE_CHECK
102  assert(a1.size() == a2.size());
103 #endif
104  for (int i = 0; i < a1.size(); i++)
105  a3.e(i) = a1.e(i) - a2.e(i);
106  }
107 
108  template <class Row1, class Row2, class AT>
109  inline void sub(Vector<Row1, AT> &a1, const Vector<Row2, AT> &a2) {
110 #ifndef FMATVEC_NO_SIZE_CHECK
111  assert(a1.size() == a2.size());
112 #endif
113  for (int i = 0; i < a1.size(); i++)
114  a1.e(i) -= a2.e(i);
115  }
116 
117  // Addition
118  template <class AT, class Row1, class Row2>
119  inline Vector<Var, AT> operator+(const Vector<Row1, AT> &a, const Vector<Row2, AT> &b) {
120  Vector<Var, AT> c(a.size(), NONINIT);
121  add(a, b, c);
122  return c;
123  }
124  template <class AT, class Row>
125  inline Vector<Row, AT> operator+(const Vector<Row, AT> &a, const Vector<Row, AT> &b) {
126  Vector<Row, AT> c(a.size(), NONINIT);
127  add(a, b, c);
128  return c;
129  }
130  template <class AT, int M, class Row>
131  inline Vector<Fixed<M>, AT> operator+(const Vector<Fixed<M>, AT> &a, const Vector<Row, AT> &b) {
132  Vector<Fixed<M>, AT> c(a.size(), NONINIT);
133  add(a, b, c);
134  return c;
135  }
136  template <class AT, int M>
137  inline Vector<Fixed<M>, AT> operator+(const Vector<Fixed<M>, AT> &a, const Vector<Fixed<M>, AT> &b) {
138  Vector<Fixed<M>, AT> c(a.size(), NONINIT);
139  add(a, b, c);
140  return c;
141  }
142  template <class AT, class Row1, class Row2>
143  inline Vector<Row1, AT> operator+=(const Vector<Row1, AT> &a_, const Vector<Row2, AT> &b) {
144  Vector<Row1, AT> &a = const_cast<Vector<Row1, AT> &>(a_);
145  add(a, b);
146  return a;
147  }
148 
150  template <class AT, class Row1, class Row2>
152  Vector<Var, AT> c(a.size(), NONINIT);
153  sub(a, b, c);
154  return c;
155  }
156  template <class AT, class Row>
157  inline Vector<Row, AT> operator-(const Vector<Row, AT> &a, const Vector<Row, AT> &b) {
158  Vector<Row, AT> c(a.size(), NONINIT);
159  sub(a, b, c);
160  return c;
161  }
162  template <class AT, int M, class Row2>
163  inline Vector<Fixed<M>, AT> operator-(const Vector<Fixed<M>, AT> &a, const Vector<Row2, AT> &b) {
164  Vector<Fixed<M>, AT> c(a.size(), NONINIT);
165  sub(a, b, c);
166  return c;
167  }
168  template <class AT, int M>
169  inline Vector<Fixed<M>, AT> operator-(const Vector<Fixed<M>, AT> &a, const Vector<Fixed<M>, AT> &b) {
170  Vector<Fixed<M>, AT> c(a.size(), NONINIT);
171  sub(a, b, c);
172  return c;
173  }
174  template <class AT, class Row1, class Row2>
175  inline Vector<Row1, AT> operator-=(const Vector<Row1, AT> &a_, const Vector<Row2, AT> &b) {
176  Vector<Row1, AT> &a = const_cast<Vector<Row1, AT> &>(a_);
177  sub(a, b);
178  return a;
179  }
180 
181  // RowVector-Vector
182  template <class Col1, class Col2, class Col3, class AT>
183  inline void add(const RowVector<Col1, AT> &a1, const RowVector<Col2, AT> &a2, RowVector<Col3, AT> &a3) {
184 #ifndef FMATVEC_NO_SIZE_CHECK
185  assert(a1.size() == a2.size());
186 #endif
187  for (int i = 0; i < a1.size(); i++)
188  a3.e(i) = a1.e(i) + a2.e(i);
189  }
190 
191  template <class Col1, class Col2, class AT>
192  inline void add(RowVector<Col1, AT> &a1, const RowVector<Col2, AT> &a2) {
193 #ifndef FMATVEC_NO_SIZE_CHECK
194  assert(a1.size() == a2.size());
195 #endif
196  for (int i = 0; i < a1.size(); i++)
197  a1.e(i) += a2.e(i);
198  }
199 
200  template <class Col1, class Col2, class Col3, class AT>
201  inline void sub(const RowVector<Col1, AT> &a1, const RowVector<Col2, AT> &a2, RowVector<Col3, AT> &a3) {
202 #ifndef FMATVEC_NO_SIZE_CHECK
203  assert(a1.size() == a2.size());
204 #endif
205  for (int i = 0; i < a1.size(); i++)
206  a3.e(i) = a1.e(i) - a2.e(i);
207  }
208 
209  template <class Col1, class Col2, class AT>
210  inline void sub(RowVector<Col1, AT> &a1, const RowVector<Col2, AT> &a2) {
211 #ifndef FMATVEC_NO_SIZE_CHECK
212  assert(a1.size() == a2.size());
213 #endif
214  for (int i = 0; i < a1.size(); i++)
215  a1.e(i) -= a2.e(i);
216  }
217 
218  // Addition
219  template <class AT, class Col1, class Col2>
220  inline RowVector<Var, AT> operator+(const RowVector<Col1, AT> &a, const RowVector<Col2, AT> &b) {
221  RowVector<Var, AT> c(a.size(), NONINIT);
222  add(a, b, c);
223  return c;
224  }
225  template <class AT, class Col>
226  inline RowVector<Col, AT> operator+(const RowVector<Col, AT> &a, const RowVector<Col, AT> &b) {
227  RowVector<Col, AT> c(a.size(), NONINIT);
228  add(a, b, c);
229  return c;
230  }
231  template <class AT, int N, class Col2>
232  inline RowVector<Fixed<N>, AT> operator+(const RowVector<Fixed<N>, AT> &a, const RowVector<Col2, AT> &b) {
233  RowVector<Fixed<N>, AT> c(a.size(), NONINIT);
234  add(a, b, c);
235  return c;
236  }
237  template <class AT, int N>
238  inline RowVector<Fixed<N>, AT> operator+(const RowVector<Fixed<N>, AT> &a, const RowVector<Fixed<N>, AT> &b) {
239  RowVector<Fixed<N>, AT> c(a.size(), NONINIT);
240  add(a, b, c);
241  return c;
242  }
243  template <class AT, class Col1, class Col2>
244  inline RowVector<Col1, AT> operator+=(const RowVector<Col1, AT> &a_, const RowVector<Col2, AT> &b) {
245  RowVector<Col1, AT> &a = const_cast<RowVector<Col1, AT> &>(a_);
246  add(a, b);
247  return a;
248  }
249 
251  template <class AT, class Col1, class Col2>
253  RowVector<Var, AT> c(a.size(), NONINIT);
254  sub(a, b, c);
255  return c;
256  }
257  template <class AT, class Col>
258  inline RowVector<Col, AT> operator-(const RowVector<Col, AT> &a, const RowVector<Col, AT> &b) {
259  RowVector<Col, AT> c(a.size(), NONINIT);
260  sub(a, b, c);
261  return c;
262  }
263  template <class AT, int N, class Col2>
264  inline RowVector<Fixed<N>, AT> operator-(const RowVector<Fixed<N>, AT> &a, const RowVector<Col2, AT> &b) {
265  RowVector<Fixed<N>, AT> c(a.size(), NONINIT);
266  sub(a, b, c);
267  return c;
268  }
269  template <class AT, int N>
270  inline RowVector<Fixed<N>, AT> operator-(const RowVector<Fixed<N>, AT> &a, const RowVector<Fixed<N>, AT> &b) {
271  RowVector<Fixed<N>, AT> c(a.size(), NONINIT);
272  sub(a, b, c);
273  return c;
274  }
275  template <class AT, class Col1, class Col2>
276  inline RowVector<Col1, AT> operator-=(const RowVector<Col1, AT> &a_, const RowVector<Col2, AT> &b) {
277  RowVector<Col1, AT> &a = const_cast<RowVector<Col1, AT> &>(a_);
278  sub(a, b);
279  return a;
280  }
281 
282  // Matrix-Matrix
283  template <class Type1, class Row1, class Col1, class Type2, class Row2, class Col2, class Type3, class Row3, class Col3, class AT>
284  inline void add(const Matrix<Type1, Row1, Col1, AT> &A1, const Matrix<Type2, Row2, Col2, AT> &A2, Matrix<Type3, Row3, Col3, AT> &A3) {
285 #ifndef FMATVEC_NO_SIZE_CHECK
286  assert(A1.rows() == A2.rows());
287  assert(A1.cols() == A2.cols());
288 #endif
289  for (int i = 0; i < A1.rows(); i++)
290  for (int j = 0; j < A2.cols(); j++)
291  A3.e(i, j) = A1.e(i, j) + A2.e(i, j);
292  }
293 
294  template <class Row1, class Row2, class Row3, class AT>
295  inline void add(const Matrix<Symmetric, Row1, Row1, AT> &A1, const Matrix<Symmetric, Row2, Row2, AT> &A2, Matrix<Symmetric, Row3, Row3, AT> &A3) {
296 #ifndef FMATVEC_NO_SIZE_CHECK
297  assert(A1.size() == A2.size());
298 #endif
299  for (int i = 0; i < A1.size(); i++)
300  for (int j = i; j < A2.size(); j++)
301  A3.ej(i, j) = A1.ej(i, j) + A2.ej(i, j);
302  }
303 
304  template <class Row1, class Row2, class Row3, class AT>
305  inline void add(const Matrix<Diagonal, Row1, Row1, AT> &A1, const Matrix<Diagonal, Row2, Row2, AT> &A2, Matrix<Diagonal, Row3, Row3, AT> &A3) {
306 #ifndef FMATVEC_NO_SIZE_CHECK
307  assert(A1.size() == A2.size());
308 #endif
309  for (int i = 0; i < A3.size(); i++)
310  A3.e(i) = A1.e(i) + A2.e(i);
311  }
312 
313  template <class Type1, class Row1, class Col1, class Type2, class Row2, class Col2, class AT>
314  inline void add(Matrix<Type1, Row1, Col1, AT> &A1, const Matrix<Type2, Row2, Col2, AT> &A2) {
315 #ifndef FMATVEC_NO_SIZE_CHECK
316  assert(A1.rows() == A2.rows());
317  assert(A1.cols() == A2.cols());
318 #endif
319  for (int i = 0; i < A1.rows(); i++)
320  for (int j = 0; j < A2.cols(); j++)
321  A1.e(i, j) += A2.e(i, j);
322  }
323 
324  template <class Row1, class Row2, class AT>
325  inline void add(Matrix<Symmetric, Row1, Row1, AT> &A1, const Matrix<Symmetric, Row2, Row2, AT> &A2) {
326 #ifndef FMATVEC_NO_SIZE_CHECK
327  assert(A1.size() == A2.size());
328 #endif
329  for (int i = 0; i < A1.size(); i++)
330  for (int j = i; j < A2.size(); j++)
331  A1.ej(i, j) += A2.ej(i, j);
332  }
333 
334  template <class Type1, class Row1, class Col1, class Type2, class Row2, class Col2, class Type3, class Row3, class Col3, class AT>
335  inline void sub(const Matrix<Type1, Row1, Col1, AT> &A1, const Matrix<Type2, Row2, Col2, AT> &A2, Matrix<Type3, Row3, Col3, AT> &A3) {
336 #ifndef FMATVEC_NO_SIZE_CHECK
337  assert(A1.rows() == A2.rows());
338  assert(A1.cols() == A2.cols());
339 #endif
340  for (int i = 0; i < A1.rows(); i++)
341  for (int j = 0; j < A2.cols(); j++)
342  A3.e(i, j) = A1.e(i, j) - A2.e(i, j);
343  }
344 
345  template <class Row1, class Row2, class Row3, class AT>
346  inline void sub(const Matrix<Symmetric, Row1, Row1, AT> &A1, const Matrix<Symmetric, Row2, Row2, AT> &A2, Matrix<Symmetric, Row3, Row3, AT> &A3) {
347 #ifndef FMATVEC_NO_SIZE_CHECK
348  assert(A1.size() == A2.size());
349 #endif
350  for (int i = 0; i < A1.size(); i++)
351  for (int j = i; j < A2.size(); j++)
352  A3.ej(i, j) = A1.ej(i, j) - A2.ej(i, j);
353  }
354 
355  template <class Row1, class Row2, class Row3, class AT>
356  inline void sub(const Matrix<Diagonal, Row1, Row1, AT> &A1, const Matrix<Diagonal, Row2, Row2, AT> &A2, Matrix<Diagonal, Row3, Row3, AT> &A3) {
357 #ifndef FMATVEC_NO_SIZE_CHECK
358  assert(A1.size() == A2.size());
359 #endif
360  for (int i = 0; i < A3.size(); i++)
361  A3.e(i) = A1.e(i) - A2.e(i);
362  }
363 
364  template <class Type1, class Row1, class Col1, class Type2, class Row2, class Col2, class AT>
365  inline void sub(Matrix<Type1, Row1, Col1, AT> &A1, const Matrix<Type2, Row2, Col2, AT> &A2) {
366 #ifndef FMATVEC_NO_SIZE_CHECK
367  assert(A1.rows() == A2.rows());
368  assert(A1.cols() == A2.cols());
369 #endif
370  for (int i = 0; i < A1.rows(); i++)
371  for (int j = 0; j < A2.cols(); j++)
372  A1.e(i, j) -= A2.e(i, j);
373  }
374 
375  template <class Row1, class Row2, class AT>
376  inline void sub(Matrix<Symmetric, Row1, Row1, AT> &A1, const Matrix<Symmetric, Row2, Row2, AT> &A2) {
377 #ifndef FMATVEC_NO_SIZE_CHECK
378  assert(A1.size() == A2.size());
379 #endif
380  for (int i = 0; i < A1.size(); i++)
381  for (int j = i; j < A2.size(); j++)
382  A1.ej(i, j) -= A2.ej(i, j);
383  }
384 
385  // Addition
386  // Type1 Type2
387  template <class AT, class Type1, class Type2, class Row1, class Col1, class Row2, class Col2>
388  inline Matrix<General, Var, Var, AT> operator+(const Matrix<Type1, Row1, Col1, AT> &A, const Matrix<Type2, Row2, Col2, AT> &B) {
389  Matrix<General, Var, Var, AT> C(A.rows(), A.cols(), NONINIT);
390  add(A, B, C);
391  return C;
392  }
393  template <class AT, class Type1, class Type2, class Row, class Col>
394  inline Matrix<General, Row, Col, AT> operator+(const Matrix<Type1, Row, Col, AT> &A, const Matrix<Type2, Row, Col, AT> &B) {
395  Matrix<General, Row, Col, AT> C(A.rows(), A.cols(), NONINIT);
396  add(A, B, C);
397  return C;
398  }
399  template <class AT, int M, class Type1, class Type2, class Row, class Col>
400  inline Matrix<General, Fixed<M>, Var, AT> operator+(const Matrix<Type1, Fixed<M>, Var, AT> &A, const Matrix<Type2, Row, Col, AT> &B) {
401  Matrix<General, Fixed<M>, Var, AT> C(A.rows(), A.cols(), NONINIT);
402  add(A, B, C);
403  return C;
404  }
405  template <class AT, int M, class Type1, class Type2>
406  inline Matrix<General, Fixed<M>, Var, AT> operator+(const Matrix<Type1, Fixed<M>, Var, AT> &A, const Matrix<Type2, Fixed<M>, Var, AT> &B) {
407  Matrix<General, Fixed<M>, Var, AT> C(A.rows(), A.cols(), NONINIT);
408  add(A, B, C);
409  return C;
410  }
411  template <class AT, int N, class Type1, class Type2, class Row, class Col>
412  inline Matrix<General, Var, Fixed<N>, AT> operator+(const Matrix<Type1, Var, Fixed<N>, AT> &A, const Matrix<Type2, Row, Col, AT> &B) {
413  Matrix<General, Var, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
414  add(A, B, C);
415  return C;
416  }
417  template <class AT, int N, class Type1, class Type2>
418  inline Matrix<General, Var, Fixed<N>, AT> operator+(const Matrix<Type1, Var, Fixed<N>, AT> &A, const Matrix<Type2, Var, Fixed<N>, AT> &B) {
419  Matrix<General, Var, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
420  add(A, B, C);
421  return C;
422  }
423  template <class AT, int M, int N, class Type1, class Type2, class Row, class Col>
424  inline Matrix<General, Fixed<M>, Fixed<N>, AT> operator+(const Matrix<Type1, Fixed<M>, Fixed<N>, AT> &A, const Matrix<Type2, Row, Col, AT> &B) {
425  Matrix<General, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
426  add(A, B, C);
427  return C;
428  }
429  template <class AT, int M, int N, class Type1, class Type2>
430  inline Matrix<General, Fixed<M>, Fixed<N>, AT> operator+(const Matrix<Type1, Fixed<M>, Fixed<N>, AT> &A, const Matrix<Type2, Fixed<M>, Var, AT> &B) {
431  Matrix<General, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
432  add(A, B, C);
433  return C;
434  }
435  template <class AT, int M, int N, class Type1, class Type2>
436  inline Matrix<General, Fixed<M>, Fixed<N>, AT> operator+(const Matrix<Type1, Fixed<M>, Fixed<N>, AT> &A, const Matrix<Type2, Var, Fixed<N>, AT> &B) {
437  Matrix<General, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
438  add(A, B, C);
439  return C;
440  }
441  template <class AT, int M, int N, class Type1, class Type2>
442  inline Matrix<General, Fixed<M>, Fixed<N>, AT> operator+(const Matrix<Type1, Fixed<M>, Fixed<N>, AT> &A, const Matrix<Type2, Fixed<M>, Fixed<N>, AT> &B) {
443  Matrix<General, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
444  add(A, B, C);
445  return C;
446  }
447  template <class AT, int M, int N, class Type1, class Type2, class Row, class Col>
448  inline Matrix<General, Fixed<M>, Fixed<N>, AT> operator+(const Matrix<Type2, Row, Col, AT> &A, const Matrix<Type1, Fixed<M>, Fixed<N>, AT> &B) {
449  Matrix<General, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
450  add(A, B, C);
451  return C;
452  }
453  template <class AT, int M, int N, class Type1, class Type2>
454  inline Matrix<General, Fixed<M>, Fixed<N>, AT> operator+(const Matrix<Type2, Fixed<M>, Var, AT> &A, const Matrix<Type1, Fixed<M>, Fixed<N>, AT> &B) {
455  Matrix<General, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
456  add(A, B, C);
457  return C;
458  }
459  template <class AT, int M, int N, class Type1, class Type2>
460  inline Matrix<General, Fixed<M>, Fixed<N>, AT> operator+(const Matrix<Type2, Var, Fixed<N>, AT> &A, const Matrix<Type1, Fixed<M>, Fixed<N>, AT> &B) {
461  Matrix<General, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
462  add(A, B, C);
463  return C;
464  }
465  // Type
466  template <class AT, class Type, class Row1, class Col1, class Row2, class Col2>
467  inline Matrix<Type, Var, Var, AT> operator+(const Matrix<Type, Row1, Col1, AT> &A, const Matrix<Type, Row2, Col2, AT> &B) {
468  Matrix<Type, Var, Var, AT> C(A.rows(), A.cols(), NONINIT);
469  add(A, B, C);
470  return C;
471  }
472  template <class AT, class Type, class Row, class Col>
473  inline Matrix<Type, Row, Col, AT> operator+(const Matrix<Type, Row, Col, AT> &A, const Matrix<Type, Row, Col, AT> &B) {
474  Matrix<Type, Row, Col, AT> C(A.rows(), A.cols(), NONINIT);
475  add(A, B, C);
476  return C;
477  }
478  template <class AT, int M, class Type, class Row, class Col>
479  inline Matrix<Type, Fixed<M>, Var, AT> operator+(const Matrix<Type, Fixed<M>, Var, AT> &A, const Matrix<Type, Row, Col, AT> &B) {
480  Matrix<Type, Fixed<M>, Var, AT> C(A.rows(), A.cols(), NONINIT);
481  add(A, B, C);
482  return C;
483  }
484  template <class AT, int M, class Type>
485  inline Matrix<Type, Fixed<M>, Var, AT> operator+(const Matrix<Type, Fixed<M>, Var, AT> &A, const Matrix<Type, Fixed<M>, Var, AT> &B) {
486  Matrix<Type, Fixed<M>, Var, AT> C(A.rows(), A.cols(), NONINIT);
487  add(A, B, C);
488  return C;
489  }
490  template <class AT, int N, class Type, class Row, class Col>
491  inline Matrix<Type, Var, Fixed<N>, AT> operator+(const Matrix<Type, Var, Fixed<N>, AT> &A, const Matrix<Type, Row, Col, AT> &B) {
492  Matrix<Type, Var, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
493  add(A, B, C);
494  return C;
495  }
496  template <class AT, int N, class Type>
497  inline Matrix<Type, Var, Fixed<N>, AT> operator+(const Matrix<Type, Var, Fixed<N>, AT> &A, const Matrix<Type, Var, Fixed<N>, AT> &B) {
498  Matrix<Type, Var, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
499  add(A, B, C);
500  return C;
501  }
502  template <class AT, int M, int N, class Type, class Row, class Col>
503  inline Matrix<Type, Fixed<M>, Fixed<N>, AT> operator+(const Matrix<Type, Fixed<M>, Fixed<N>, AT> &A, const Matrix<Type, Row, Col, AT> &B) {
504  Matrix<Type, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
505  add(A, B, C);
506  return C;
507  }
508  template <class AT, int M, int N, class Type>
509  inline Matrix<Type, Fixed<M>, Fixed<N>, AT> operator+(const Matrix<Type, Fixed<M>, Fixed<N>, AT> &A, const Matrix<Type, Fixed<M>, Var, AT> &B) {
510  Matrix<Type, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
511  add(A, B, C);
512  return C;
513  }
514  template <class AT, int M, int N, class Type>
515  inline Matrix<Type, Fixed<M>, Fixed<N>, AT> operator+(const Matrix<Type, Fixed<M>, Fixed<N>, AT> &A, const Matrix<Type, Var, Fixed<N>, AT> &B) {
516  Matrix<Type, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
517  add(A, B, C);
518  return C;
519  }
520  template <class AT, int M, int N, class Type>
521  inline Matrix<Type, Fixed<M>, Fixed<N>, AT> operator+(const Matrix<Type, Fixed<M>, Fixed<N>, AT> &A, const Matrix<Type, Fixed<M>, Fixed<N>, AT> &B) {
522  Matrix<Type, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
523  add(A, B, C);
524  return C;
525  }
526  template <class AT, int M, int N, class Type, class Row, class Col>
527  inline Matrix<Type, Fixed<M>, Fixed<N>, AT> operator+(const Matrix<Type, Row, Col, AT> &A, const Matrix<Type, Fixed<M>, Fixed<N>, AT> &B) {
528  Matrix<Type, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
529  add(A, B, C);
530  return C;
531  }
532  template <class AT, int M, int N, class Type>
533  inline Matrix<Type, Fixed<M>, Fixed<N>, AT> operator+(const Matrix<Type, Fixed<M>, Var, AT> &A, const Matrix<Type, Fixed<M>, Fixed<N>, AT> &B) {
534  Matrix<Type, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
535  add(A, B, C);
536  return C;
537  }
538  template <class AT, int M, int N, class Type>
539  inline Matrix<Type, Fixed<M>, Fixed<N>, AT> operator+(const Matrix<Type, Var, Fixed<N>, AT> &A, const Matrix<Type, Fixed<M>, Fixed<N>, AT> &B) {
540  Matrix<Type, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
541  add(A, B, C);
542  return C;
543  }
544  template <class AT, class Type1, class Row1, class Col1, class Type2, class Row2, class Col2>
545  inline Matrix<Type1, Row1, Col1, AT>& operator+=(const Matrix<Type1, Row1, Col1, AT> &A_, const Matrix<Type2, Row2, Col2, AT> &B) {
546  Matrix<Type1, Row1, Col1, AT> &A = const_cast<Matrix<Type1, Row1, Col1, AT> &>(A_);
547  add(A, B);
548  return A;
549  }
550 
551  // Subtraction
552  // Type1 Type2
553  template <class AT, class Type1, class Type2, class Row1, class Col1, class Row2, class Col2>
554  inline Matrix<General, Var, Var, AT> operator-(const Matrix<Type1, Row1, Col1, AT> &A, const Matrix<Type2, Row2, Col2, AT> &B) {
555  Matrix<General, Var, Var, AT> C(A.rows(), A.cols(), NONINIT);
556  sub(A, B, C);
557  return C;
558  }
559  template <class AT, class Type1, class Type2, class Row, class Col>
560  inline Matrix<General, Row, Col, AT> operator-(const Matrix<Type1, Row, Col, AT> &A, const Matrix<Type2, Row, Col, AT> &B) {
561  Matrix<General, Row, Col, AT> C(A.rows(), A.cols(), NONINIT);
562  sub(A, B, C);
563  return C;
564  }
565  template <class AT, int M, class Type1, class Type2, class Row, class Col>
566  inline Matrix<General, Fixed<M>, Var, AT> operator-(const Matrix<Type1, Fixed<M>, Var, AT> &A, const Matrix<Type2, Row, Col, AT> &B) {
567  Matrix<General, Fixed<M>, Var, AT> C(A.rows(), A.cols(), NONINIT);
568  sub(A, B, C);
569  return C;
570  }
571  template <class AT, int M, class Type1, class Type2>
572  inline Matrix<General, Fixed<M>, Var, AT> operator-(const Matrix<Type1, Fixed<M>, Var, AT> &A, const Matrix<Type2, Fixed<M>, Var, AT> &B) {
573  Matrix<General, Fixed<M>, Var, AT> C(A.rows(), A.cols(), NONINIT);
574  sub(A, B, C);
575  return C;
576  }
577  template <class AT, int N, class Type1, class Type2, class Row, class Col>
578  inline Matrix<General, Var, Fixed<N>, AT> operator-(const Matrix<Type1, Var, Fixed<N>, AT> &A, const Matrix<Type2, Row, Col, AT> &B) {
579  Matrix<General, Var, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
580  sub(A, B, C);
581  return C;
582  }
583  template <class AT, int N, class Type1, class Type2>
584  inline Matrix<General, Var, Fixed<N>, AT> operator-(const Matrix<Type1, Var, Fixed<N>, AT> &A, const Matrix<Type2, Var, Fixed<N>, AT> &B) {
585  Matrix<General, Var, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
586  sub(A, B, C);
587  return C;
588  }
589  template <class AT, int M, int N, class Type1, class Type2, class Row, class Col>
590  inline Matrix<General, Fixed<M>, Fixed<N>, AT> operator-(const Matrix<Type1, Fixed<M>, Fixed<N>, AT> &A, const Matrix<Type2, Row, Col, AT> &B) {
591  Matrix<General, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
592  sub(A, B, C);
593  return C;
594  }
595  template <class AT, int M, int N, class Type1, class Type2>
596  inline Matrix<General, Fixed<M>, Fixed<N>, AT> operator-(const Matrix<Type1, Fixed<M>, Fixed<N>, AT> &A, const Matrix<Type2, Fixed<M>, Var, AT> &B) {
597  Matrix<General, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
598  sub(A, B, C);
599  return C;
600  }
601  template <class AT, int M, int N, class Type1, class Type2>
602  inline Matrix<General, Fixed<M>, Fixed<N>, AT> operator-(const Matrix<Type1, Fixed<M>, Fixed<N>, AT> &A, const Matrix<Type2, Var, Fixed<N>, AT> &B) {
603  Matrix<General, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
604  sub(A, B, C);
605  return C;
606  }
607  template <class AT, int M, int N, class Type1, class Type2>
608  inline Matrix<General, Fixed<M>, Fixed<N>, AT> operator-(const Matrix<Type1, Fixed<M>, Fixed<N>, AT> &A, const Matrix<Type2, Fixed<M>, Fixed<N>, AT> &B) {
609  Matrix<General, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
610  sub(A, B, C);
611  return C;
612  }
613  template <class AT, int M, int N, class Type1, class Type2, class Row, class Col>
614  inline Matrix<General, Fixed<M>, Fixed<N>, AT> operator-(const Matrix<Type2, Row, Col, AT> &A, const Matrix<Type1, Fixed<M>, Fixed<N>, AT> &B) {
615  Matrix<General, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
616  sub(A, B, C);
617  return C;
618  }
619  template <class AT, int M, int N, class Type1, class Type2>
620  inline Matrix<General, Fixed<M>, Fixed<N>, AT> operator-(const Matrix<Type2, Fixed<M>, Var, AT> &A, const Matrix<Type1, Fixed<M>, Fixed<N>, AT> &B) {
621  Matrix<General, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
622  sub(A, B, C);
623  return C;
624  }
625  template <class AT, int M, int N, class Type1, class Type2>
626  inline Matrix<General, Fixed<M>, Fixed<N>, AT> operator-(const Matrix<Type2, Var, Fixed<N>, AT> &A, const Matrix<Type1, Fixed<M>, Fixed<N>, AT> &B) {
627  Matrix<General, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
628  sub(A, B, C);
629  return C;
630  }
631  // Type
632  template <class AT, class Type, class Row1, class Col1, class Row2, class Col2>
633  inline Matrix<Type, Var, Var, AT> operator-(const Matrix<Type, Row1, Col1, AT> &A, const Matrix<Type, Row2, Col2, AT> &B) {
634  Matrix<Type, Var, Var, AT> C(A.rows(), A.cols(), NONINIT);
635  sub(A, B, C);
636  return C;
637  }
638  template <class AT, class Type, class Row, class Col>
639  inline Matrix<Type, Row, Col, AT> operator-(const Matrix<Type, Row, Col, AT> &A, const Matrix<Type, Row, Col, AT> &B) {
640  Matrix<Type, Row, Col, AT> C(A.rows(), A.cols(), NONINIT);
641  sub(A, B, C);
642  return C;
643  }
644  template <class AT, int M, class Type, class Row, class Col>
645  inline Matrix<Type, Fixed<M>, Var, AT> operator-(const Matrix<Type, Fixed<M>, Var, AT> &A, const Matrix<Type, Row, Col, AT> &B) {
646  Matrix<Type, Fixed<M>, Var, AT> C(A.rows(), A.cols(), NONINIT);
647  sub(A, B, C);
648  return C;
649  }
650  template <class AT, int M, class Type>
651  inline Matrix<Type, Fixed<M>, Var, AT> operator-(const Matrix<Type, Fixed<M>, Var, AT> &A, const Matrix<Type, Fixed<M>, Var, AT> &B) {
652  Matrix<Type, Fixed<M>, Var, AT> C(A.rows(), A.cols(), NONINIT);
653  sub(A, B, C);
654  return C;
655  }
656  template <class AT, int N, class Type, class Row, class Col>
657  inline Matrix<Type, Var, Fixed<N>, AT> operator-(const Matrix<Type, Var, Fixed<N>, AT> &A, const Matrix<Type, Row, Col, AT> &B) {
658  Matrix<Type, Var, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
659  sub(A, B, C);
660  return C;
661  }
662  template <class AT, int N, class Type>
663  inline Matrix<Type, Var, Fixed<N>, AT> operator-(const Matrix<Type, Var, Fixed<N>, AT> &A, const Matrix<Type, Var, Fixed<N>, AT> &B) {
664  Matrix<Type, Var, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
665  sub(A, B, C);
666  return C;
667  }
668  template <class AT, int M, int N, class Type, class Row, class Col>
669  inline Matrix<Type, Fixed<M>, Fixed<N>, AT> operator-(const Matrix<Type, Fixed<M>, Fixed<N>, AT> &A, const Matrix<Type, Row, Col, AT> &B) {
670  Matrix<Type, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
671  sub(A, B, C);
672  return C;
673  }
674  template <class AT, int M, int N, class Type>
675  inline Matrix<Type, Fixed<M>, Fixed<N>, AT> operator-(const Matrix<Type, Fixed<M>, Fixed<N>, AT> &A, const Matrix<Type, Fixed<M>, Var, AT> &B) {
676  Matrix<Type, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
677  sub(A, B, C);
678  return C;
679  }
680  template <class AT, int M, int N, class Type>
681  inline Matrix<Type, Fixed<M>, Fixed<N>, AT> operator-(const Matrix<Type, Fixed<M>, Fixed<N>, AT> &A, const Matrix<Type, Var, Fixed<N>, AT> &B) {
682  Matrix<Type, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
683  sub(A, B, C);
684  return C;
685  }
686  template <class AT, int M, int N, class Type>
687  inline Matrix<Type, Fixed<M>, Fixed<N>, AT> operator-(const Matrix<Type, Fixed<M>, Fixed<N>, AT> &A, const Matrix<Type, Fixed<M>, Fixed<N>, AT> &B) {
688  Matrix<Type, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
689  sub(A, B, C);
690  return C;
691  }
692  template <class AT, int M, int N, class Type, class Row, class Col>
693  inline Matrix<Type, Fixed<M>, Fixed<N>, AT> operator-(const Matrix<Type, Row, Col, AT> &A, const Matrix<Type, Fixed<M>, Fixed<N>, AT> &B) {
694  Matrix<Type, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
695  sub(A, B, C);
696  return C;
697  }
698  template <class AT, int M, int N, class Type>
699  inline Matrix<Type, Fixed<M>, Fixed<N>, AT> operator-(const Matrix<Type, Fixed<M>, Var, AT> &A, const Matrix<Type, Fixed<M>, Fixed<N>, AT> &B) {
700  Matrix<Type, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
701  sub(A, B, C);
702  return C;
703  }
704  template <class AT, int M, int N, class Type>
705  inline Matrix<Type, Fixed<M>, Fixed<N>, AT> operator-(const Matrix<Type, Var, Fixed<N>, AT> &A, const Matrix<Type, Fixed<M>, Fixed<N>, AT> &B) {
706  Matrix<Type, Fixed<M>, Fixed<N>, AT> C(A.rows(), A.cols(), NONINIT);
707  sub(A, B, C);
708  return C;
709  }
710  template <class AT, class Type1, class Row1, class Col1, class Type2, class Row2, class Col2>
711  inline Matrix<Type1, Row1, Col1, AT>& operator-=(const Matrix<Type1, Row1, Col1, AT> &A_, const Matrix<Type2, Row2, Col2, AT> &B) {
712  Matrix<Type1, Row1, Col1, AT> &A = const_cast<Matrix<Type1, Row1, Col1, AT> &>(A_);
713  sub(A, B);
714  return A;
715  }
716 
717 // SquareMatrix-SquareMatrix
718  template <class AT, class Row1, class Row2>
719  inline SquareMatrix<Var, AT> operator+(const SquareMatrix<Row1, AT> &A1, const SquareMatrix<Row2, AT> &A2) {
720  SquareMatrix<Var, AT> A3(A1.size(), NONINIT);
721  add(A1, A2, A3);
722  return A3;
723  }
724  template <class AT, class Row>
725  inline SquareMatrix<Row, AT> operator+(const SquareMatrix<Row, AT> &A1, const SquareMatrix<Row, AT> &A2) {
726  SquareMatrix<Row, AT> A3(A1.size(), NONINIT);
727  add(A1, A2, A3);
728  return A3;
729  }
730 
731 // SquareMatrix-SquareMatrix
732  template <class AT, class Row1, class Row2>
733  inline SquareMatrix<Var, AT> operator-(const SquareMatrix<Row1, AT> &A1, const SquareMatrix<Row2, AT> &A2) {
734  SquareMatrix<Var, AT> A3(A1.size(), NONINIT);
735  sub(A1, A2, A3);
736  return A3;
737  }
738  template <class AT, class Row>
739  inline SquareMatrix<Row, AT> operator-(const SquareMatrix<Row, AT> &A1, const SquareMatrix<Row, AT> &A2) {
740  SquareMatrix<Row, AT> A3(A1.size(), NONINIT);
741  sub(A1, A2, A3);
742  return A3;
743  }
744 
750  template <class Type1, class Row1, class Col1, class Row2, class Row3, class AT>
751  inline void mult(const Matrix<Type1, Row1, Col1, AT> &A, const Vector<Row2, AT> &x, Vector<Row3, AT> &y) {
752 #ifndef FMATVEC_NO_SIZE_CHECK
753  assert(A.cols() == x.size());
754 #endif
755  for (int i = 0; i < y.size(); i++) {
756  y.e(i) = 0;
757  for (int j = 0; j < x.size(); j++)
758  y.e(i) += A.e(i, j) * x.e(j);
759  }
760  }
761 
762  template <class Row1, class Row2, class Row3, class AT>
763  inline void mult(const Matrix<Symmetric, Row1, Row1, AT> &A, const Vector<Row2, AT> &x, Vector<Row3, AT> &y) {
764 #ifndef FMATVEC_NO_SIZE_CHECK
765  assert(A.cols() == x.size());
766 #endif
767  for (int i = 0; i < y.size(); i++) {
768  y.e(i) = 0;
769  for (int j = 0; j < i; j++)
770  y.e(i) += A.ei(i, j) * x.e(j);
771  for (int j = i; j < x.size(); j++)
772  y.e(i) += A.ej(i, j) * x.e(j);
773  }
774  }
775 
776  // Matrix-Vector
777  template <class AT, class Type1, class Row1, class Col1, class Row2>
778  inline Vector<Var, AT> operator*(const Matrix<Type1, Row1, Col1, AT> &A, const Vector<Row2, AT> &x) {
779  Vector<Var, AT> y(A.rows(), NONINIT);
780  mult(A, x, y);
781  return y;
782  }
783  template <class AT, class Type1, class Row, class Col>
784  inline Vector<Row, AT> operator*(const Matrix<Type1, Row, Col, AT> &A, const Vector<Row, AT> &x) {
785  Vector<Row, AT> y(A.rows(), NONINIT);
786  mult(A, x, y);
787  return y;
788  }
789  template <class AT, int M, class Type1, class Col1, class Row2>
790  inline Vector<Fixed<M>, AT> operator*(const Matrix<Type1, Fixed<M>, Col1, AT> &A, const Vector<Row2, AT> &x) {
791  Vector<Fixed<M>, AT> y(A.rows(), NONINIT);
792  mult(A, x, y);
793  return y;
794  }
795  template <class AT, int M, class Type1>
796  inline Vector<Fixed<M>, AT> operator*(const Matrix<Type1, Fixed<M>, Fixed<M>, AT> &A, const Vector<Fixed<M>, AT> &x) {
797  Vector<Fixed<M>, AT> y(A.rows(), NONINIT);
798  mult(A, x, y);
799  return y;
800  }
801  template <class AT, int M, int N, class Type, class Row2>
802  inline Vector<Fixed<M>, AT> operator*(const Matrix<Type, Fixed<M>, Fixed<N>, AT> &A, const Vector<Row2, AT> &x) {
803  Vector<Fixed<M>, AT> y(A.rows(), NONINIT);
804  mult(A, x, y);
805  return y;
806  }
807  template <class AT, int M, int N, class Type>
808  inline Vector<Fixed<M>, AT> operator*(const Matrix<Type, Fixed<M>, Fixed<N>, AT> &A, const Vector<Fixed<N>, AT> &x) {
809  Vector<Fixed<M>, AT> y(A.rows(), NONINIT);
810  mult(A, x, y);
811  return y;
812  }
813 
814  // RowVector-Matrix
815  template <class Col1, class Type2, class Row2, class Col2, class Type3, class Col3, class AT>
816  inline void mult(const RowVector<Col1, AT> &x, const Matrix<Type2, Row2, Col2, AT> &A, RowVector<Col3, AT> &y) {
817 #ifndef FMATVEC_NO_SIZE_CHECK
818  assert(x.size() == A.rows());
819 #endif
820  for (int i = 0; i < y.size(); i++) {
821  y.e(i) = 0;
822  for (int j = 0; j < x.size(); j++)
823  y.e(i) += x.e(j) * A.e(j, i);
824  }
825  }
826 
827  template <class Col1, class Row2, class Col3, class AT>
828  inline void mult(const RowVector<Col1, AT> &x, const Matrix<Symmetric, Row2, Row2, AT> &A, RowVector<Col3, AT> &y) {
829 #ifndef FMATVEC_NO_SIZE_CHECK
830  assert(x.size() == A.rows());
831 #endif
832  for (int i = 0; i < y.size(); i++) {
833  y.e(i) = 0;
834  for (int j = 0; j < i; j++)
835  y.e(i) += x.e(j) * A.ej(j, i);
836  for (int j = i; j < x.size(); j++)
837  y.e(i) += x.e(j) * A.ei(j, i);
838  }
839  }
840 
841  template <class AT, class Type2, class Col2, class Row1, class Col1>
842  inline RowVector<Var, AT> operator*(const RowVector<Col2, AT> &x, const Matrix<Type2, Row1, Col1, AT> &A) {
843  RowVector<Var, AT> y(A.cols(), NONINIT);
844  mult(x, A, y);
845  return y;
846  }
847  template <class AT, class Type2, class Row, class Col>
848  inline RowVector<Col, AT> operator*(const RowVector<Col, AT> &x, const Matrix<Type2, Row, Col, AT> &A) {
849  RowVector<Col, AT> y(A.cols(), NONINIT);
850  mult(x, A, y);
851  return y;
852  }
853  template <class AT, int N, class Type2, class Col1, class Row2>
854  inline RowVector<Fixed<N>, AT> operator*(const RowVector<Col1, AT> &x, const Matrix<Type2, Row2, Fixed<N>, AT> &A) {
855  RowVector<Fixed<N>, AT> y(A.cols(), NONINIT);
856  mult(x, A, y);
857  return y;
858  }
859  template <class AT, int N, class Type2>
860  inline RowVector<Fixed<N>, AT> operator*(const RowVector<Fixed<N>, AT> &x, const Matrix<Type2, Fixed<N>, Fixed<N>, AT> &A) {
861  RowVector<Fixed<N>, AT> y(A.cols(), NONINIT);
862  mult(x, A, y);
863  return y;
864  }
865 
866  // Matrix-Matrix
867  template <class Type1, class Row1, class Col1, class Type2, class Row2, class Col2, class Type3, class Row3, class Col3, class AT>
868  inline void mult(const Matrix<Type1, Row1, Col1, AT> &A1, const Matrix<Type2, Row2, Col2, AT> &A2, Matrix<Type3, Row3, Col3, AT> &A3) {
869 #ifndef FMATVEC_NO_SIZE_CHECK
870  assert(A1.cols() == A2.rows());
871 #endif
872  for (int i = 0; i < A3.rows(); i++) {
873  for (int k = 0; k < A3.cols(); k++) {
874  A3.e(i, k) = 0;
875  for (int j = 0; j < A1.cols(); j++)
876  A3.e(i, k) += A1.e(i, j) * A2.e(j, k);
877  }
878  }
879  }
880  template <class Row1, class Type2, class Row2, class Col2, class Type3, class Row3, class Col3, class AT>
881  inline void mult(const Matrix<Symmetric, Row1, Row1, AT> &A1, const Matrix<Type2, Row2, Col2, AT> &A2, Matrix<Type3, Row3, Col3, AT> &A3) {
882 #ifndef FMATVEC_NO_SIZE_CHECK
883  assert(A1.size() == A2.rows());
884 #endif
885  for (int i = 0; i < A3.rows(); i++) {
886  for (int k = 0; k < A3.cols(); k++) {
887  A3.e(i, k) = 0;
888  for (int j = 0; j < i; j++)
889  A3.e(i, k) += A1.ei(i, j) * A2.e(j, k);
890  for (int j = i; j < A1.cols(); j++)
891  A3.e(i, k) += A1.ej(i, j) * A2.e(j, k);
892  }
893  }
894  }
895  template <class Row1, class Row2, class Row3, class AT>
896  inline void mult(const Matrix<Diagonal, Row1, Row1, AT> &A1, const Matrix<Diagonal, Row2, Row2, AT> &A2, Matrix<Diagonal, Row3, Row3, AT> &A3) {
897 #ifndef FMATVEC_NO_SIZE_CHECK
898  assert(A1.size() == A2.size());
899 #endif
900  for (int i = 0; i < A3.size(); i++)
901  A3.e(i) = A1.e(i) * A2.e(i);
902  }
903 
904  template <class AT, class Type1, class Type2, class Row1, class Col1, class Row2, class Col2>
905  inline Matrix<General, Var, Var, AT> operator*(const Matrix<Type1, Row1, Col1, AT> &A1, const Matrix<Type2, Row2, Col2, AT> &A2) {
906  Matrix<General, Var, Var, AT> A3(A1.rows(), A2.cols(), NONINIT);
907  mult(A1, A2, A3);
908  return A3;
909  }
910  template <class AT, class Type1, class Type2, class Row, class Col>
911  inline Matrix<General, Row, Col, AT> operator*(const Matrix<Type1, Row, Col, AT> &A1, const Matrix<Type2, Row, Col, AT> &A2) {
912  Matrix<General, Row, Col, AT> A3(A1.rows(), A2.cols(), NONINIT);
913  mult(A1, A2, A3);
914  return A3;
915  }
916  template <class AT, int M, class Type1, class Type2, class Col1, class Row2, class Col2>
917  inline Matrix<General, Fixed<M>, Var, AT> operator*(const Matrix<Type1, Fixed<M>, Col1, AT> &A1, const Matrix<Type2, Row2, Col2, AT> &A2) {
918  Matrix<General, Fixed<M>, Var, AT> A3(A1.rows(), A2.cols(), NONINIT);
919  mult(A1, A2, A3);
920  return A3;
921  }
922  template <class AT, int N, class Type1, class Row1, class Col1, class Type2, class Row2>
923  inline Matrix<General, Var, Fixed<N>, AT> operator*(const Matrix<Type1, Row1, Col1, AT> &A1, const Matrix<Type2, Row2, Fixed<N>, AT> &A2) {
924  Matrix<General, Var, Fixed<N>, AT> A3(A1.rows(), A2.cols(), NONINIT);
925  mult(A1, A2, A3);
926  return A3;
927  }
928  template <class AT, int M, int N, class Type1, class Type2, class Col1, class Row2>
929  inline Matrix<General, Fixed<M>, Fixed<N>, AT> operator*(const Matrix<Type1, Fixed<M>, Col1, AT> &A1, const Matrix<Type2, Row2, Fixed<N>, AT> &A2) {
930  Matrix<General, Fixed<M>, Fixed<N>, AT> A3(A1.rows(), A2.cols(), NONINIT);
931  mult(A1, A2, A3);
932  return A3;
933  }
934  template <class AT, int M, class Type1, class Type2>
935  inline Matrix<General, Fixed<M>, Fixed<M>, AT> operator*(const Matrix<Type1, Fixed<M>, Fixed<M>, AT> &A1, const Matrix<Type2, Fixed<M>, Fixed<M>, AT> &A2) {
936  Matrix<General, Fixed<M>, Fixed<M>, AT> A3(A1.rows(), A2.cols(), NONINIT);
937  mult(A1, A2, A3);
938  return A3;
939  }
940  template <class AT, class Type, class Row1, class Col1, class Row2, class Col2>
941  inline Matrix<Type, Var, Var, AT> operator*(const Matrix<Type, Row1, Col1, AT> &A1, const Matrix<Type, Row2, Col2, AT> &A2) {
942  Matrix<Type, Var, Var, AT> A3(A1.rows(), A2.cols(), NONINIT);
943  mult(A1, A2, A3);
944  return A3;
945  }
946  template <class AT, class Type, class Row, class Col>
947  inline Matrix<Type, Row, Col, AT> operator*(const Matrix<Type, Row, Col, AT> &A1, const Matrix<Type, Row, Col, AT> &A2) {
948  Matrix<Type, Row, Col, AT> A3(A1.rows(), A2.cols(), NONINIT);
949  mult(A1, A2, A3);
950  return A3;
951  }
952  template <class AT, int M, class Type, class Col1, class Row2, class Col2>
953  inline Matrix<Type, Fixed<M>, Var, AT> operator*(const Matrix<Type, Fixed<M>, Col1, AT> &A1, const Matrix<Type, Row2, Col2, AT> &A2) {
954  Matrix<Type, Fixed<M>, Var, AT> A3(A1.rows(), A2.cols(), NONINIT);
955  mult(A1, A2, A3);
956  return A3;
957  }
958  template <class AT, int N, class Type, class Row1, class Col1, class Row2>
959  inline Matrix<Type, Var, Fixed<N>, AT> operator*(const Matrix<Type, Row1, Col1, AT> &A1, const Matrix<Type, Row2, Fixed<N>, AT> &A2) {
960  Matrix<Type, Var, Fixed<N>, AT> A3(A1.rows(), A2.cols(), NONINIT);
961  mult(A1, A2, A3);
962  return A3;
963  }
964  template <class AT, int M, int N, class Type, class Col1, class Row2>
965  inline Matrix<Type, Fixed<M>, Fixed<N>, AT> operator*(const Matrix<Type, Fixed<M>, Col1, AT> &A1, const Matrix<Type, Row2, Fixed<N>, AT> &A2) {
966  Matrix<Type, Fixed<M>, Fixed<N>, AT> A3(A1.rows(), A2.cols(), NONINIT);
967  mult(A1, A2, A3);
968  return A3;
969  }
970  template <class AT, int M, class Type>
971  inline Matrix<Type, Fixed<M>, Fixed<M>, AT> operator*(const Matrix<Type, Fixed<M>, Fixed<M>, AT> &A1, const Matrix<Type, Fixed<M>, Fixed<M>, AT> &A2) {
972  Matrix<Type, Fixed<M>, Fixed<M>, AT> A3(A1.rows(), A2.cols(), NONINIT);
973  mult(A1, A2, A3);
974  return A3;
975  }
976 
977 // SquareMatrix-SquareMatrix
978  template <class AT, class Row1, class Row2>
979  inline SquareMatrix<Var, AT> operator*(const SquareMatrix<Row1, AT> &A1, const SquareMatrix<Row2, AT> &A2) {
980  SquareMatrix<Var, AT> A3(A1.size(), NONINIT);
981  mult(A1, A2, A3);
982  return A3;
983  }
984  template <class AT, class Row>
985  inline SquareMatrix<Row, AT> operator*(const SquareMatrix<Row, AT> &A1, const SquareMatrix<Row, AT> &A2) {
986  SquareMatrix<Row, AT> A3(A1.size(), NONINIT);
987  mult(A1, A2, A3);
988  return A3;
989  }
990 
991  // SymmetricMatrix-SymmetricMatrix
992  template <class AT, class Row1, class Row2>
993  inline Matrix<General, Var, Var, AT> operator*(const Matrix<Symmetric, Row1, Row1, AT> &A1, const Matrix<Symmetric, Row2, Row2, AT> &A2) {
994  Matrix<General, Var, Var, AT> A3(A1.rows(), A2.cols(), NONINIT);
995  mult(A1, A2, A3);
996  return A3;
997  }
998  template <class AT, class Row>
999  inline Matrix<General, Row, Row, AT> operator*(const Matrix<Symmetric, Row, Row, AT> &A1, const Matrix<Symmetric, Row, Row, AT> &A2) {
1000  Matrix<General, Row, Row, AT> A3(A1.rows(), A2.cols(), NONINIT);
1001  mult(A1, A2, A3);
1002  return A3;
1003  }
1004 
1006 
1013  template <class Row, class AT>
1014  Vector<Row, AT> operator*(const Vector<Row, AT> &x, const AT& alpha) {
1015 
1016  Vector<Row, AT> y(x.size(), NONINIT);
1017 
1018  for (int i = 0; i < x.size(); i++)
1019  y.e(i) = x.e(i) * alpha;
1020 
1021  return y;
1022  }
1023 
1028  template <class Row, class AT>
1029  Vector<Row, AT> operator*(const AT& alpha, const Vector<Row, AT> &x) {
1030 
1031  Vector<Row, AT> y(x.size(), NONINIT);
1032 
1033  for (int i = 0; i < x.size(); i++)
1034  y.e(i) = x.e(i) * alpha;
1035 
1036  return y;
1037  }
1038 
1039  template <class Row, class AT>
1040  inline Vector<Row, AT> operator*=(const Vector<Row, AT> &x_, const AT &alpha) {
1041  Vector<Row, AT> &x = const_cast<Vector<Row, AT> &>(x_);
1042  for (int i = 0; i < x.size(); i++)
1043  x.e(i) *= alpha;
1044  return x;
1045  }
1046 
1047  template <class Row, class AT>
1048  inline Vector<Row, AT> operator/(const Vector<Row, AT> &x, const AT &alpha) {
1049  Vector<Row, AT> y(x.size(), NONINIT);
1050  for (int i = 0; i < x.size(); i++)
1051  y.e(i) = x.e(i) / alpha;
1052  return y;
1053  }
1054 
1055  template <class Row, class AT>
1056  inline Vector<Row, AT> operator/=(const Vector<Row, AT> &x_, const AT &alpha) {
1057  Vector<Row, AT> &x = const_cast<Vector<Row, AT> &>(x_);
1058  for (int i = 0; i < x.size(); i++)
1059  x.e(i) /= alpha;
1060  return x;
1061  }
1068  template <class Col, class AT>
1069  RowVector<Col, AT> operator*(const RowVector<Col, AT> &x, const AT& alpha) {
1070 
1071  RowVector<Col, AT> y(x.size(), NONINIT);
1072 
1073  for (int i = 0; i < x.size(); i++)
1074  y.e(i) = x.e(i) * alpha;
1075 
1076  return y;
1077  }
1078 
1083  template <class Col, class AT>
1084  RowVector<Col, AT> operator*(const AT &alpha, const RowVector<Col, AT> &x) {
1085 
1086  RowVector<Col, AT> y(x.size(), NONINIT);
1087 
1088  for (int i = 0; i < x.size(); i++)
1089  y.e(i) = x.e(i) * alpha;
1090 
1091  return y;
1092  }
1093 
1094  template <class Col, class AT>
1095  inline RowVector<Col, AT> operator*=(const RowVector<Col, AT> &x_, const AT &alpha) {
1096  RowVector<Col, AT> &x = const_cast<RowVector<Col, AT> &>(x_);
1097  for (int i = 0; i < x.size(); i++)
1098  x.e(i) *= alpha;
1099  return x;
1100  }
1101 
1102  template <class Col, class AT>
1103  inline RowVector<Col, AT> operator/(const RowVector<Col, AT> &x, const AT &a) {
1104  RowVector<Col, AT> y(x.size(), NONINIT);
1105  for (int i = 0; i < x.size(); i++)
1106  y.e(i) = x.e(i) / a;
1107  return y;
1108  }
1109 
1110  template <class Col, class AT>
1111  inline RowVector<Col, AT> operator/=(const RowVector<Col, AT> &x_, const AT &a) {
1112  RowVector<Col, AT> &x = const_cast<RowVector<Col, AT> &>(x_);
1113  for (int i = 0; i < x.size(); i++)
1114  x.e(i) /= a;
1115  return x;
1116  }
1117 
1119 
1121 
1128  template <class Col, class Row, class AT>
1129  AT operator*(const RowVector<Col, AT> &x, const Vector<Row, AT> &y) {
1130 
1131 #ifndef FMATVEC_NO_SIZE_CHECK
1132  assert(x.size() == y.size());
1133 #endif
1134 
1135  AT res = 0;
1136 
1137  for (int i = 0; i < x.size(); i++)
1138  res += x.e(i) * y.e(i);
1139 
1140  return res;
1141  }
1142 
1148  template <class Row, class AT>
1150 
1151 #ifndef FMATVEC_NO_SIZE_CHECK
1152  assert(x.size() == y.size());
1153 #endif
1154 
1155  AT res = 0;
1156 
1157  for (int i = 0; i < x.size(); i++)
1158  res += x.e(i) * y.e(i);
1159 
1160  return res;
1161  }
1162 
1168  template <class Row, class AT>
1170 
1171 #ifndef FMATVEC_NO_SIZE_CHECK
1172  assert(x.size() == 3);
1173  assert(y.size() == 3);
1174 #endif
1175 
1176  Vector<Row, AT> z(3, NONINIT);
1177 
1178  z.e(0) = x.e(1) * y.e(2) - x.e(2) * y.e(1);
1179  z.e(1) = x.e(2) * y.e(0) - x.e(0) * y.e(2);
1180  z.e(2) = x.e(0) * y.e(1) - x.e(1) * y.e(0);
1181 
1182  return z;
1183  }
1184 
1190  template <class Row, class AT>
1191  double tripleProduct(const Vector<Row, AT> &a, const Vector<Row, AT> &x, const Vector<Row, AT> &y) {
1192 
1193 #ifndef FMATVEC_NO_SIZE_CHECK
1194  assert(a.size() == 3);
1195  assert(x.size() == 3);
1196  assert(y.size() == 3);
1197 #endif
1198 
1199  return a.e(0) * (x.e(1) * y.e(2) - x.e(2) * y.e(1)) + a.e(1) * (x.e(2) * y.e(0) - x.e(0) * y.e(2)) + a.e(2) * (x.e(0) * y.e(1) - x.e(1) * y.e(0));
1200  }
1201 
1203  //
1205 
1214  template <class Type, class Row, class Col, class AT>
1215  Matrix<Type, Row, Col, AT> operator*(const Matrix<Type, Row, Col, AT> &A, const AT &alpha) {
1216 
1217  Matrix<Type, Row, Col, AT> B(A.rows(), A.cols(), NONINIT);
1218 
1219  for (int i = 0; i < A.rows(); i++)
1220  for (int j = 0; j < A.cols(); j++)
1221  B.e(i, j) = A.e(i, j) * alpha;
1222 
1223  return B;
1224  }
1225 
1226  template <class Type, class Row, class Col, class AT>
1227  Matrix<Type, Row, Col, AT> operator*(const AT &alpha, const Matrix<Type, Row, Col, AT> &A) {
1228 
1229  Matrix<Type, Row, Col, AT> B(A.rows(), A.cols(), NONINIT);
1230 
1231  for (int i = 0; i < A.rows(); i++)
1232  for (int j = 0; j < A.cols(); j++)
1233  B.e(i, j) = A.e(i, j) * alpha;
1234 
1235  return B;
1236  }
1237 
1238  template <class Row, class AT>
1239  inline Matrix<Symmetric, Row, Row, AT> operator*(const AT &alpha, const Matrix<Symmetric, Row, Row, AT> &A) {
1240  Matrix<Symmetric, Row, Row, AT> B(A.size(), A.size(), NONINIT);
1241  for (int i = 0; i < A.size(); i++)
1242  for (int j = i; j < A.size(); j++)
1243  B.ej(i, j) = A.ej(i, j) * alpha;
1244  return B;
1245  }
1246 
1247  template <class Row, class AT>
1248  inline Matrix<Symmetric, Row, Row, AT> operator*(const Matrix<Symmetric, Row, Row, AT> &A, const AT &alpha) {
1249  Matrix<Symmetric, Row, Row, AT> B(A.size(), A.size(), NONINIT);
1250  for (int i = 0; i < A.size(); i++)
1251  for (int j = i; j < A.size(); j++)
1252  B.ej(i, j) = A.ej(i, j) * alpha;
1253  return B;
1254  }
1255 
1256  template <class Row, class AT>
1257  inline Matrix<Diagonal, Row, Row, AT> operator*(const AT &alpha, const Matrix<Diagonal, Row, Row, AT> &A) {
1258  Matrix<Diagonal, Row, Row, AT> B(A.size(), A.size(), NONINIT);
1259  for (int i = 0; i < A.size(); i++)
1260  B.e(i) = A.e(i) * alpha;
1261  return B;
1262  }
1263 
1264  template <class Row, class AT>
1265  inline Matrix<Diagonal, Row, Row, AT> operator*(const Matrix<Diagonal, Row, Row, AT> &A, const AT &alpha) {
1266  Matrix<Diagonal, Row, Row, AT> B(A.size(), A.size(), NONINIT);
1267  for (int i = 0; i < A.size(); i++)
1268  B.e(i) = A.e(i) * alpha;
1269  return B;
1270  }
1271 
1272  template <class Type, class Row, class Col, class AT>
1273  Matrix<Type, Row, Col, AT> operator/(const Matrix<Type, Row, Col, AT> &A, const AT &alpha) {
1274 
1275  Matrix<Type, Row, Col, AT> B(A.rows(), A.cols(), NONINIT);
1276 
1277  for (int i = 0; i < A.rows(); i++)
1278  for (int j = 0; j < A.cols(); j++)
1279  B.e(i, j) = A.e(i, j) / alpha;
1280 
1281  return B;
1282  }
1283 
1284  template <class Row, class AT>
1285  inline Matrix<Symmetric, Row, Row, AT> operator/(const Matrix<Symmetric, Row, Row, AT> &A, const AT &alpha) {
1286  Matrix<Symmetric, Row, Row, AT> B(A.size(), A.size(), NONINIT);
1287  for (int i = 0; i < A.size(); i++)
1288  for (int j = i; j < A.size(); j++)
1289  B.ej(i, j) = A.ej(i, j) / alpha;
1290  return B;
1291  }
1292 
1293  template <class Row, class AT>
1294  inline Matrix<Diagonal, Row, Row, AT> operator/(const Matrix<Diagonal, Row, Row, AT> &A, const AT &alpha) {
1295  Matrix<Diagonal, Row, Row, AT> B(A.size(), A.size(), NONINIT);
1296  for (int i = 0; i < A.size(); i++)
1297  B.e(i) = A.e(i) / alpha;
1298  return B;
1299  }
1300 
1301  template <class Row, class AT>
1302  SquareMatrix<Row, AT> operator*(const SquareMatrix<Row, AT> &A, const AT &alpha) {
1303 
1304  SquareMatrix<Row, AT> B(A.size(), NONINIT);
1305 
1306  for (int i = 0; i < A.size(); i++)
1307  for (int j = 0; j < A.size(); j++)
1308  B.e(i, j) = A.e(i, j) * alpha;
1309 
1310  return B;
1311  }
1312 
1313  template <class Row, class AT>
1314  SquareMatrix<Row, AT> operator*(const AT &alpha, const SquareMatrix<Row, AT> &A) {
1315 
1316  SquareMatrix<Row, AT> B(A.size(), NONINIT);
1317 
1318  for (int i = 0; i < A.size(); i++)
1319  for (int j = 0; j < A.size(); j++)
1320  B.e(i, j) = A.e(i, j) * alpha;
1321 
1322  return B;
1323  }
1324 
1325  template <class Row, class AT>
1326  SquareMatrix<Row, AT> operator/(const SquareMatrix<Row, AT> &A, const AT &alpha) {
1327 
1328  SquareMatrix<Row, AT> B(A.size(), NONINIT);
1329 
1330  for (int i = 0; i < A.size(); i++)
1331  for (int j = 0; j < A.size(); j++)
1332  B.e(i, j) = A.e(i, j) / alpha;
1333 
1334  return B;
1335  }
1336 
1337  template <class Type, class Row, class Col, class AT>
1338  Matrix<Type, Row, Col, AT> operator*=(const Matrix<Type, Row, Col, AT> &A_, const AT &alpha) {
1339  Matrix<Type, Row, Col, AT> &A = const_cast<Matrix<Type, Row, Col, AT> &>(A_);
1340  for (int i = 0; i < A.rows(); i++)
1341  for (int j = 0; j < A.cols(); j++)
1342  A.e(i, j) *= alpha;
1343 
1344  return A;
1345  }
1346 
1347  template <class Row, class AT>
1348  inline Matrix<Symmetric, Row, Row, AT> operator*=(const Matrix<Symmetric, Row, Row, AT> &A_, const AT &alpha) {
1349  Matrix<Symmetric, Row, Row, AT> &A = const_cast<Matrix<Symmetric, Row, Row, AT> &>(A_);
1350  for (int i = 0; i < A.size(); i++)
1351  for (int j = i; j < A.size(); j++)
1352  A.ej(i, j) *= alpha;
1353  return A;
1354  }
1355 
1356  template <class Row, class AT>
1357  inline Matrix<Diagonal, Row, Row, AT> operator*=(const Matrix<Diagonal, Row, Row, AT> &A_, const AT &alpha) {
1358  Matrix<Diagonal, Row, Row, AT> &A = const_cast<Matrix<Diagonal, Row, Row, AT> &>(A_);
1359  for (int i = 0; i < A.size(); i++)
1360  A.e(i) *= alpha;
1361  return A;
1362  }
1363 
1364  template <class Type, class Row, class Col, class AT>
1365  Matrix<Type, Row, Col, AT> operator/=(const Matrix<Type, Row, Col, AT> &A_, const AT &alpha) {
1366  Matrix<Type, Row, Col, AT> A = const_cast<Matrix<Type, Row, Col, AT> &>(A_);
1367  for (int i = 0; i < A.rows(); i++)
1368  for (int j = 0; j < A.cols(); j++)
1369  A.e(i, j) /= alpha;
1370 
1371  return A;
1372  }
1373 
1374  template <class Row, class AT>
1375  inline Matrix<Symmetric, Row, Row, AT> operator/=(const Matrix<Symmetric, Row, Row, AT> &A_, const AT &alpha) {
1376  Matrix<Symmetric, Row, Row, AT> A = const_cast<Matrix<Symmetric, Row, Row, AT> &>(A_);
1377  for (int i = 0; i < A.size(); i++)
1378  for (int j = i; j < A.size(); j++)
1379  A.ej(i, j) /= alpha;
1380  return A;
1381  }
1382 
1383  template <class Row, class AT>
1384  inline Matrix<Diagonal, Row, Row, AT> operator/=(const Matrix<Diagonal, Row, Row, AT> &A_, const AT &alpha) {
1385  Matrix<Diagonal, Row, Row, AT> A = const_cast<Matrix<Diagonal, Row, Row, AT> &>(A_);
1386  for (int i = 0; i < A.size(); i++)
1387  A.e(i) /= alpha;
1388  return A;
1389  }
1390 
1392 
1394 
1400  template <class Row, class AT>
1402 
1403  Vector<Row, AT> y(x.size(), NONINIT);
1404 
1405  for (int i = 0; i < x.size(); i++)
1406  y.e(i) = -x.e(i);
1407 
1408  return y;
1409  }
1410 
1416  template <class Col, class AT>
1418 
1419  RowVector<Col, AT> c(a.size(), NONINIT);
1420 
1421  for (int i = 0; i < a.size(); i++)
1422  c.e(i) = -a.e(i);
1423 
1424  return c;
1425  }
1426 
1432  template <class Row, class AT>
1434 
1435  SquareMatrix<Row, AT> B(A.size(), NONINIT);
1436 
1437  for (int i = 0; i < A.size(); i++)
1438  for (int j = 0; j < A.size(); j++)
1439  B.e(i, j) = -A.e(i, j);
1440 
1441  return B;
1442  }
1443 
1449  template <class Type, class Row, class Col, class AT>
1451 
1452  Matrix<Type, Row, Col, AT> B(A.rows(), A.cols(), NONINIT);
1453 
1454  for (int i = 0; i < A.rows(); i++)
1455  for (int j = 0; j < A.cols(); j++)
1456  B.e(i, j) = -A.e(i, j);
1457 
1458  return B;
1459  }
1460 
1463 
1469  template <class AT>
1471 
1472  return RowVector<Ref, AT>(x.m, x.lda, x.tp ? false : true, x.memory, x.ele).copy();
1473  }
1474 
1482  template <class AT>
1484 
1485  return Vector<Ref, AT>(x.n, x.lda, x.tp ? false : true, x.memory, x.ele).copy();
1486  }
1487 
1496  template <class AT>
1498 
1499  return Matrix<General, Ref, Ref, AT>(A.n, A.m, A.lda, A.tp ? false : true, A.memory, A.ele).copy();
1500  }
1501 
1510  template <class AT>
1512 
1513  return SquareMatrix<Ref, AT>(A.n, A.lda, A.tp ? false : true, A.memory, A.ele).copy();
1514  }
1515 
1516  template <class Type, class Row, class Col, class AT>
1517  inline Matrix<Type, Col, Row, AT> trans(const Matrix<Type, Row, Col, AT> &A) {
1518  Matrix<Type, Col, Row, AT> B(A.cols(), A.rows(), NONINIT);
1519  for (int i = 0; i < B.rows(); i++)
1520  for (int j = 0; j < B.cols(); j++)
1521  B.e(i, j) = A.e(j, i);
1522  return B;
1523  }
1524 
1525  template <class Row, class AT>
1526  inline SquareMatrix<Row, AT> trans(const SquareMatrix<Row, AT> &A) {
1527  SquareMatrix<Row, AT> B(A.size(), NONINIT);
1528  for (int i = 0; i < B.size(); i++)
1529  for (int j = 0; j < B.size(); j++)
1530  B.e(i, j) = A.e(j, i);
1531  return B;
1532  }
1533 
1534  template <class Row, class AT>
1535  inline RowVector<Row, AT> trans(const Vector<Row, AT> &x) {
1536  RowVector<Row, AT> y(x.size(), NONINIT);
1537  for (int i = 0; i < y.size(); i++)
1538  y.e(i) = x.e(i);
1539  return y;
1540  }
1541 
1542  template <class Row, class AT>
1543  inline Vector<Row, AT> trans(const RowVector<Row, AT> &x) {
1544  Vector<Row, AT> y(x.size(), NONINIT);
1545  for (int i = 0; i < y.size(); i++)
1546  y.e(i) = x.e(i);
1547  return y;
1548  }
1549 
1551 
1553  template <class Row, class AT>
1554  inline AT nrmInf(const Vector<Row, AT> &x) {
1555  AT c = 0;
1556  for (int i = 0; i < x.size(); i++)
1557  if (c < fabs(x.e(i)))
1558  c = fabs(x.e(i));
1559  return c;
1560  }
1561 
1562  template <class Row, class AT>
1563  inline AT nrmInf(const RowVector<Row, AT> &x) {
1564  AT c = 0;
1565  for (int i = 0; i < x.size(); i++)
1566  if (c < fabs(x.e(i)))
1567  c = fabs(x.e(i));
1568  return c;
1569  }
1570 
1571  template <class Row, class AT>
1572  inline AT nrm2(const Vector<Row, AT> &x) {
1573  AT c = 0;
1574  for (int i = 0; i < x.size(); i++)
1575  c += pow(x.e(i), 2);
1576  return sqrt(c);
1577  }
1578 
1579  template <class Col, class AT>
1580  inline AT nrm2(const RowVector<Col, AT> &x) {
1581  AT c = 0;
1582  for (int i = 0; i < x.size(); i++)
1583  c += pow(x.e(i), 2);
1584  return sqrt(c);
1585  }
1586 
1588 
1600  template <class Row, class AT>
1602 
1603  SquareMatrix<Row, AT> B(x.size(), NONINIT);
1604 
1605  B.e(0, 0) = 0;
1606  B.e(1, 1) = 0;
1607  B.e(2, 2) = 0;
1608  B.e(0, 1) = -x.e(2);
1609  B.e(0, 2) = x.e(1);
1610  B.e(1, 0) = x.e(2);
1611  B.e(1, 2) = -x.e(0);
1612  B.e(2, 0) = -x.e(1);
1613  B.e(2, 1) = x.e(0);
1614 
1615  return B;
1616  }
1617 
1623  template <class AT, class Type1, class Row1, class Col1, class Type2, class Row2, class Col2>
1625 
1626  if (&A == &B)
1627  return true;
1628 
1629  if (A.rows() != B.rows() || A.cols() != B.cols())
1630  return false;
1631 
1632  for (int i = 0; i < A.rows(); i++)
1633  for (int j = 0; j < A.cols(); j++)
1634  if (A(i, j) != B(i, j))
1635  return false;
1636 
1637  return true;
1638  }
1639 
1645  template <class AT, class Type1, class Row1, class Col1, class Type2, class Row2, class Col2>
1647 
1648  if (A.rows() != B.rows() || A.cols() != B.cols())
1649  return true;
1650 
1651  for (int i = 0; i < A.rows(); i++)
1652  for (int j = 0; j < A.cols(); j++)
1653  if (A(i, j) != B(i, j))
1654  return true;
1655 
1656  return false;
1657  }
1658 
1664  template <class Row, class AT>
1665  AT max(const Vector<Row, AT> &x) {
1666 
1667 #ifndef FMATVEC_NO_SIZE_CHECK
1668  assert(x.size() > 0);
1669 #endif
1670  AT maximum = x.e(0);
1671  for (int i = 1; i < x.size(); i++) {
1672  if (x.e(i) > maximum)
1673  maximum = x.e(i);
1674  }
1675  return maximum;
1676  }
1677 
1683  template <class Row, class AT>
1684  int maxIndex(const Vector<Row, AT> &x) {
1685 
1686 #ifndef FMATVEC_NO_SIZE_CHECK
1687  assert(x.size() > 0);
1688 #endif
1689  AT maximum = x.e(0);
1690  int index = 0;
1691  for (int i = 1; i < x.size(); i++) {
1692  if (x.e(i) > maximum) {
1693  maximum = x.e(i);
1694  index = i;
1695  }
1696  }
1697  return index;
1698  }
1699 
1705  template <class Row, class AT>
1706  AT min(const Vector<Row, AT> &x) {
1707 
1708 #ifndef FMATVEC_NO_SIZE_CHECK
1709  assert(x.size() > 0);
1710 #endif
1711  AT minimum = x.e(0);
1712  for (int i = 1; i < x.size(); i++) {
1713  if (x.e(i) < minimum)
1714  minimum = x.e(i);
1715  }
1716  return minimum;
1717  }
1718 
1724  template <class Row, class AT>
1725  int minIndex(const Vector<Row, AT> &x) {
1726 
1727 #ifndef FMATVEC_NO_SIZE_CHECK
1728  assert(x.size() > 0);
1729 #endif
1730  AT minimum = x.e(0);
1731  int index = 0;
1732  for (int i = 1; i < x.size(); i++) {
1733  if (x.e(i) < minimum) {
1734  minimum = x.e(i);
1735  index = i;
1736  }
1737  }
1738  return index;
1739  }
1740 
1741  // HR 28.09.2006
1748  template <class AT>
1750 
1751 #ifndef FMATVEC_NO_SIZE_CHECK
1752  assert(A_.rows() > 0);
1753  assert(A_.cols() > 0);
1754  assert(A_.cols() > PivotCol);
1755 #endif
1757  int i, j, N;
1758  RowVector<Ref, AT> tmp(A.cols(), NONINIT);
1759  N = A.rows();
1760  for (i = 1; i <= N - 1; i++) {
1761  for (j = 0; j < N - 1; j++) {
1762  if (A(j, PivotCol) > A(j + 1, PivotCol)) {
1763  tmp = A.row(j);
1764  A.row(j) = A.row(j + 1);
1765  A.row(j + 1) = tmp;
1766  }
1767  }
1768  }
1769  return A;
1770  }
1771 
1772  // HR 28.09.2006
1774  template <class AT>
1775  void quicksortmedian_intern(Matrix<General, Ref, Ref, AT> &A, int PivotCol, RowVector<Ref, AT> &tmp, int l, int r) {
1776 
1777  if (r > l) {
1778  int i = l - 1, j = r;
1779  //RowVec tmp;
1780  if (r - l > 3) { //Median of three
1781  int m = l + (r - l) / 2;
1782  if (A(l, PivotCol) > A(m, PivotCol)) {
1783  tmp = A.row(l);
1784  A.row(l) = A.row(m);
1785  A.row(m) = tmp;
1786  }
1787  if (A(l, PivotCol) > A(r, PivotCol)) {
1788  tmp = A.row(l);
1789  A.row(l) = A.row(r);
1790  A.row(r) = tmp;
1791  }
1792  else if (A(r, PivotCol) > A(m, PivotCol)) {
1793  tmp = A.row(r);
1794  A.row(r) = A.row(m);
1795  A.row(m) = tmp;
1796  }
1797  }
1798 
1799  for (;;) {
1800  while(A(++i,PivotCol) < A(r,PivotCol));
1801  while(A(--j,PivotCol) > A(r,PivotCol) && j>i);
1802  if(i>=j) break;
1803  tmp = A.row(i);
1804  A.row(i) = A.row(j);
1805  A.row(j) = tmp;
1806  }
1807  tmp = A.row(i);
1808  A.row(i) = A.row(r);
1809  A.row(r) = tmp;
1810  quicksortmedian_intern(A, PivotCol, tmp, l, i - 1);
1811  quicksortmedian_intern(A, PivotCol, tmp, i + 1, r);
1812  }
1813  }
1814 
1815  // HR 28.09.2006
1824  template <class AT>
1827  int N = A.rows();
1828  RowVector<Ref, AT> tmp(A.cols(), NONINIT);
1829  quicksortmedian_intern(A, PivotCol, tmp, 0, N - 1);
1830  return A;
1831  }
1832 
1839  template <class AT> int countElements(const SquareMatrix<Ref, AT> &A) {
1840  int k = 0;
1841  for (int i = 0; i < A.size(); i++) {
1842  for (int j = 0; j < A.size(); j++) {
1843  if (fabs(A(i, j)) > 1e-16 || i == j) {
1844  k++;
1845  }
1846  }
1847  }
1848  return k;
1849  }
1850 
1857  template <class AT> int countElementsLT(const Matrix<Symmetric,Ref,Ref,AT> &A) {
1858  int k = 0;
1859  for (int j = 0; j < A.cols(); j++) {
1860  for (int i = j; i < A.rows(); i++) {
1861  if (fabs(A(i, j)) > 1e-17 || i == j) {
1862  k++;
1863  }
1864  }
1865  }
1866  return k;
1867  }
1868 
1875  template <class AT> int countElementsLT(const Matrix<Symmetric,Var,Var,AT> &A) {
1876  int k = 0;
1877  for (int j = 0; j < A.cols(); j++) {
1878  for (int i = j; i < A.rows(); i++) {
1879  if (fabs(A(i, j)) > 1e-17 || i == j) {
1880  k++;
1881  }
1882  }
1883  }
1884  return k;
1885  }
1886 
1887  template <class Row, class Col, class AT>
1888  inline Matrix<Symmetric, Col, Col, AT> JTJ(const Matrix<General, Row, Col, AT> &A) {
1889  Matrix<Symmetric, Col, Col, AT> S(A.cols(), A.cols(), NONINIT);
1890  for (int i = 0; i < A.cols(); i++) {
1891  for (int k = i; k < A.cols(); k++) {
1892  S.ej(i, k) = 0;
1893  for (int j = 0; j < A.rows(); j++)
1894  S.ej(i, k) += A.e(j, i) * A.e(j, k);
1895  }
1896  }
1897  return S;
1898  }
1899 
1900  template <class Row, class Col, class AT>
1901  inline Matrix<Symmetric, Col, Col, AT> JTMJ(const Matrix<Symmetric, Row, Row, AT> &B, const Matrix<General, Row, Col, AT> &A) {
1902 
1903  Matrix<Symmetric, Col, Col, AT> S(A.cols(), A.cols(), NONINIT);
1904  Matrix<General, Row, Col, AT> C = B * A;
1905 
1906  for (int i = 0; i < A.cols(); i++) {
1907  for (int k = i; k < A.cols(); k++) {
1908  S.ej(i, k) = 0;
1909  for (int j = 0; j < A.rows(); j++)
1910  S.ej(i, k) += A.e(j, i) * C.e(j, k);
1911  }
1912  }
1913  return S;
1914  }
1915 
1916  template <class Row, class Col, class AT>
1917  inline Matrix<Symmetric, Row, Row, AT> JMJT(const Matrix<General, Row, Col, AT> &A, const Matrix<Symmetric, Col, Col, AT> &B) {
1918 
1919  Matrix<Symmetric, Row, Row, AT> S(A.rows(), A.rows(), NONINIT);
1920  Matrix<General, Row, Col, AT> C = A * B;
1921 
1922  for (int i = 0; i < S.size(); i++) {
1923  for (int k = i; k < S.size(); k++) {
1924  S.ej(i, k) = 0;
1925  for (int j = 0; j < A.cols(); j++)
1926  S.ej(i, k) += C.e(k, j) * A.e(i, j);
1927  }
1928  }
1929  return S;
1930  }
1931 
1932 }
1933 
1934 #endif
This is a matrix class for symmetric matrices.
Definition: var_symmetric_matrix.h:37
Matrix< General, Ref, Ref, AT > quickSortMedian(const Matrix< General, Ref, Ref, AT > &A_, int PivotCol)
Modified QuickSort Algorithm (unstable sorting Algorithm )
Definition: linear_algebra.h:1825
This is the basic matrix class for arbitrary matrices.
Definition: matrix.h:56
int rows() const
Number of rows.
Definition: var_symmetric_matrix.h:231
This is a vector class of general shape in dense storage format.
Definition: var_vector.h:39
int rows() const
Number of rows.
Definition: general_matrix.h:270
void quicksortmedian_intern(Matrix< General, Ref, Ref, AT > &A, int PivotCol, RowVector< Ref, AT > &tmp, int l, int r)
Definition: linear_algebra.h:1775
int minIndex(const Vector< Row, AT > &x)
Index of minimum value.
Definition: linear_algebra.h:1725
int rows() const
Number of rows.
Definition: symmetric_matrix.h:264
Definition: matrix.h:218
Vector< Var, AT > operator-(const Vector< Row1, AT > &a, const Vector< Row2, AT > &b)
// Subtraction
Definition: linear_algebra.h:151
RowVector< Ref, AT > row(int i)
Row operator.
Definition: general_matrix.h:680
bool operator==(const Index &I, const Index &J)
Equality operator for indices.
Definition: index.cc:26
int rows() const
Number of rows.
int countElements(const SquareMatrix< Ref, AT > &A)
Count nonzero elements.
Definition: linear_algebra.h:1839
int cols() const
Number of columns.
Vector< Ref, AT > operator+(const Vector< Ref, AT > &x, const Vector< Ref, AT > &y)
Vector-vector addition.
Definition: linear_algebra.h:52
Matrix< General, Ref, Ref, AT > copy() const
Matrix duplicating.
Definition: general_matrix.h:702
int size() const
Size.
Definition: square_matrix.h:183
This is a rowvector class of general shape in dense storage format.
Definition: row_vector.h:36
Definition: matrix.h:212
AT min(const Vector< Row, AT > &x)
Minimum value.
Definition: linear_algebra.h:1706
Definition: matrix.h:215
Matrix< General, Ref, Ref, AT > bubbleSort(const Matrix< General, Ref, Ref, AT > &A_, int PivotCol)
Bubble Sort Algorithm (stable sorting Algorithm )
Definition: linear_algebra.h:1749
int maxIndex(const Vector< Row, AT > &x)
Index of maximum value.
Definition: linear_algebra.h:1684
int cols() const
Number of columns.
Definition: symmetric_matrix.h:270
This is a matrix class for general matrices.
Definition: general_matrix.h:40
bool operator!=(const Matrix< Type1, Row1, Col1, AT > &A, const Matrix< Type2, Row2, Col2, AT > &B)
Matrix-matrix comparison.
Definition: linear_algebra.h:1646
int size() const
Size.
Definition: vector.h:260
Vector< Row, AT > crossProduct(const Vector< Row, AT > &x, const Vector< Row, AT > &y)
Cross product.
Definition: linear_algebra.h:1169
RowVector< Ref, AT > trans(const Vector< Ref, AT > &x)
Transpose of a vector.
Definition: linear_algebra.h:1470
AT max(const Vector< Row, AT > &x)
Maximum value.
Definition: linear_algebra.h:1665
This is a matrix class for symmetric matrices.
Definition: symmetric_matrix.h:39
bool transposed() const
Transposed status.
Definition: general_matrix.h:289
This is a vector class of general shape in dense storage format.
Definition: vector.h:39
double tripleProduct(const Vector< Row, AT > &a, const Vector< Row, AT > &x, const Vector< Row, AT > &y)
Triple product.
Definition: linear_algebra.h:1191
This is a matrix class of general quadratic matrices.
Definition: square_matrix.h:38
SquareMatrix< Row, AT > tilde(const Vector< Row, AT > &x)
Tilde operator.
Definition: linear_algebra.h:1601
int cols() const
Number of columns.
Definition: var_symmetric_matrix.h:237
int countElementsLT(const Matrix< Symmetric, Ref, Ref, AT > &A)
Count nonzero elements of the low triangular part of a symmetric matrix.
Definition: linear_algebra.h:1857
AT scalarProduct(const Vector< Row, AT > &x, const Vector< Row, AT > &y)
Scalar product.
Definition: linear_algebra.h:1149
This is a vector class of general shape in dense storage format.
Definition: var_row_vector.h:39
int cols() const
Number of columns.
Definition: general_matrix.h:276

Impressum / Disclaimer / Datenschutz Generated by doxygen 1.8.5 Valid HTML