GL Studio C++ Runtime API
gls_matrix.h
Go to the documentation of this file.
1 /*! \file
2  \brief The GlsMatrix class.
3 
4  \par Copyright Information
5 
6  Copyright (c) 2017 by The DiSTI Corporation.<br>
7  11301 Corporate Blvd; Suite 100<br>
8  Orlando, Florida 32817<br>
9  USA<br>
10  <br>
11  All rights reserved.<br>
12 
13  This Software contains proprietary trade secrets of DiSTI and may not be
14 reproduced, in whole or part, in any form, or by any means of electronic,
15 mechanical, or otherwise, without the written permission of DiSTI. Said
16 permission may be derived through the purchase of applicable DiSTI product
17 licenses which detail the distribution rights of this content and any
18 Derivative Works based on this or other copyrighted DiSTI Software.
19 
20  NO WARRANTY. THE SOFTWARE IS PROVIDED "AS-IS," WITHOUT WARRANTY OF ANY KIND,
21 AND ANY USE OF THIS SOFTWARE PRODUCT IS AT YOUR OWN RISK. TO THE MAXIMUM EXTENT
22 PERMITTED BY APPLICABLE LAW, DISTI AND ITS SUPPLIERS DISCLAIM ALL WARRANTIES
23 AND CONDITIONS, EITHER EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
24 IMPLIED WARRANTIES AND CONDITIONS OF MERCHANTABILITY AND/OR FITNESS FOR A
25 PARTICULAR PURPOSE, TITLE, AND NON-INFRINGEMENT, WITH REGARD TO THE SOFTWARE.
26 
27  LIMITATION OF LIABILITY. TO THE MAXIMUM EXTENT PERMITTED BY APPLICABLE LAW,
28 IN NO EVENT SHALL DISTI OR ITS SUPPLIERS BE LIABLE FOR ANY SPECIAL, INCIDENTAL,
29 INDIRECT, OR CONSEQUENTIAL DAMAGES WHATSOEVER (INCLUDING, WITHOUT LIMITATION,
30 DAMAGES FOR LOSS OF BUSINESS PROFITS, BUSINESS INTERRUPTION, LOSS OF BUSINESS
31 INFORMATION, OR ANY OTHER PECUNIARY LOSS) ARISING OUT OF THE USE OF OR
32 INABILITY TO USE THE SOFTWARE, EVEN IF DISTI HAS BEEN ADVISED OF THE POSSIBILITY
33 OF SUCH DAMAGES. DISTI'S ENTIRE LIABILITY AND YOUR EXCLUSIVE REMEDY SHALL NOT
34 EXCEED FIVE DOLLARS (US$5.00).
35 
36  The aforementioned terms and restrictions are governed by the laws of the
37 State of Florida and the United States of America.
38 
39 */
40 #ifndef GLS_MATRIX_H
41 #define GLS_MATRIX_H
42 
43 #include <cstring>
44 #include <iostream>
45 #include <math.h>
46 #include <memory>
47 
48 namespace disti
49 {
50 //----------------------------------------------------------------------------
51 /**
52 * The GlsMatrix class provides support for an n x n square matrix of any
53 * numeric type. Most basic matrix operations are supported. Internally,
54 * the matrix terms are stored contiguously in column-major order so that the
55 * matrix data can be passed to OpenGL matrix operations. However, instances
56 * of this class can be treated as row-major matrices since the API performs
57 * the conversion to column-major ordering for the user.
58 */
59 //----------------------------------------------------------------------------
60 template<class Type, int DIM = 4>
61 class GlsMatrix
62 {
63 public:
64  typedef Type CStyleMatrix[ DIM ][ DIM ];
65 
66  //------------------------------------------------------------------------
67  /**
68  * Default constructor. The new matrix will be the identity matrix.
69  */
70  //------------------------------------------------------------------------
71  GlsMatrix();
72 
73  //------------------------------------------------------------------------
74  /**
75  * Constructor. This will transpose m so that the matrix is stored in
76  * column-major order, but a caller can pass one in that is in row-major
77  * order, which is more intuitive.
78  *
79  * \param m a row-major, two-dimensional array to initialize the matrix
80  * with.
81  */
82  //------------------------------------------------------------------------
83  GlsMatrix( const CStyleMatrix m );
84 
85  //------------------------------------------------------------------------
86  /**
87  * Constructor.
88  *
89  * \param m a column-major, one-dimensional array containing all values
90  * to initialize the matrix with.
91  */
92  //------------------------------------------------------------------------
93  GlsMatrix( const Type* m );
94 
95  //------------------------------------------------------------------------
96  /**
97  * Copy Constructor.
98  *
99  * \param m matrix to copy.
100  */
101  //------------------------------------------------------------------------
102  GlsMatrix( const GlsMatrix<Type, DIM>& m );
103 
104  //------------------------------------------------------------------------
105  /**
106  * Destructor.
107  */
108  //------------------------------------------------------------------------
109  virtual ~GlsMatrix();
110 
111  //------------------------------------------------------------------------
112  /**
113  * Dump the contents of the matrix to standard output.
114  */
115  //------------------------------------------------------------------------
116  void Dump() const;
117 
118  //------------------------------------------------------------------------
119  /**
120  * Zero out the matrix so that all terms are equal to zero.
121  */
122  //------------------------------------------------------------------------
123  inline void Clear() { memset( _data, 0, sizeof( Type ) * DIM * DIM ); }
124 
125  //------------------------------------------------------------------------
126  /**
127  * Provides access to the raw column-major matrix data as a one dimensional
128  * array in a format suitable for Open GL matrix operations.
129  */
130  //------------------------------------------------------------------------
131  inline Type* Data() { return _data; }
132 
133  //------------------------------------------------------------------------
134  /**
135  * Provides access to the raw column-major matrix data as a one dimensional
136  * array in a format suitable for Open GL matrix operations.
137  */
138  //------------------------------------------------------------------------
139  inline const Type* Data() const { return _data; }
140 
141  //------------------------------------------------------------------------
142  /**
143  * Set this matrix to the identity matrix. The identity matrix is
144  * I = [Aij], where Aii = 1 and Aij = 0 for i != j. Thus, for the identity
145  * matrix, the terms off the main diagonal are all zero, and the terms on
146  * the main diagonal are all one.
147  */
148  //------------------------------------------------------------------------
149  void MakeIdentity();
150 
151  //------------------------------------------------------------------------
152  /**
153  * Calculates the determinant of the matrix.
154  *
155  * \return the calculated value of the determinant.
156  */
157  //------------------------------------------------------------------------
158  Type Determinant() const;
159 
160  //------------------------------------------------------------------------
161  /**
162  * Calculates the inverse of this matrix. The inverse of matrix A is B if
163  * AB = BA = I, where I is the identity matrix. The inverse of a matrix is
164  * the notion corresponding to the reciprocal of a nonzero real number.
165  *
166  * \return the inverse of this matrix.
167  */
168  //------------------------------------------------------------------------
170 
171  //------------------------------------------------------------------------
172  /**
173  * Sets this matrix equal to its inverse.
174  */
175  //------------------------------------------------------------------------
176  void Invert();
177 
178  //------------------------------------------------------------------------
179  /**
180  * Calculates the transpose of this matrix. The transpose of a matrix
181  * is obtained by interchanging its rows and columns.
182  *
183  * \return the transpose of this matrix.
184  */
185  //------------------------------------------------------------------------
187 
188  //------------------------------------------------------------------------
189  /**
190  * Sets this matrix equal to its transpose.
191  */
192  //------------------------------------------------------------------------
193  void Transpose();
194 
195  //------------------------------------------------------------------------
196  /**
197  * Determines if this is a diagonal matrix. An n x n matrix A = [Aij] is
198  * called a diagonal matrix if Aij = 0 for i != j. Thus, for a diagonal
199  * matrix, the terms off the main diagonal are all zero.
200  *
201  * \return true if the matrix is a diagonal matrix, false otherwise.
202  */
203  //------------------------------------------------------------------------
204  bool IsDiagonal() const;
205 
206  //------------------------------------------------------------------------
207  /**
208  * Determines if this matrix is equal to the identity matrix. The
209  * identity matrix is I = [Aij], where Aii = 1 and Aij = 0 for i != j.
210  * Thus, for the identity matrix, the terms off the main diagonal are all
211  * zero, and the terms on the main diagonal are all one.
212  *
213  * \return true if the matrix is the identity matrix, false otherwise.
214  */
215  //------------------------------------------------------------------------
216  bool IsIdentity() const;
217 
218  //------------------------------------------------------------------------
219  /**
220  * Determines if this is a scalar matrix. A scalar matrix is a diagonal
221  * matrix whose diagonal terms are equal.
222  *
223  * \return true if the matrix is a scalar matrix, false otherwise.
224  */
225  //------------------------------------------------------------------------
226  bool IsScalar() const;
227 
228  //------------------------------------------------------------------------
229  /**
230  * Determines if this matrix is singular, or non-invertible. A singular
231  * matrix has no inverse.
232  *
233  * \return true if the matrix is singular, false otherwise.
234  */
235  //------------------------------------------------------------------------
236  bool IsSingular() const;
237 
238  //------------------------------------------------------------------------
239  /**
240  * Determines if this matrix is symmetric. A matrix is symmetric if it is
241  * equal to its transpose, i.e. A = A'.
242  *
243  * \return true if the matrix is symmetric, false otherwise.
244  */
245  //------------------------------------------------------------------------
246  bool IsSymmetric() const;
247 
248  //------------------------------------------------------------------------
249  /**
250  * Determines if this matrix is skew symmetric. A matrix is skew symmetric
251  * if its transpose is equal to its negative, i.e. -A = A'.
252  *
253  * \return true if the matrix is skew symmetric, false otherwise.
254  */
255  //------------------------------------------------------------------------
256  bool IsSkewSymmetric() const;
257 
258  //------------------------------------------------------------------------
259  /**
260  * Determines if this matrix is lower triangular. An n x n matrix A = [Aij]
261  * is lower triangular if Aij = 0 for i < j.
262  *
263  * \return true if the matrix is lower triangular, false otherwise.
264  */
265  //------------------------------------------------------------------------
266  bool IsLowerTriangular() const;
267 
268  //------------------------------------------------------------------------
269  /**
270  * Determines if this matrix is upper triangular. An n x n matrix A = [Aij]
271  * is upper triangular if Aij = 0 for i > j.
272  *
273  * \return true if the matrix is upper triangular, false otherwise.
274  */
275  //------------------------------------------------------------------------
276  bool IsUpperTriangular() const;
277 
278  //------------------------------------------------------------------------
279  // Operators
280  //------------------------------------------------------------------------
281  // Subscript operator (zero-based)
282  Type& operator()( unsigned row, unsigned col );
283  Type operator()( unsigned row, unsigned col ) const;
284 
285  // Equality operators
286  bool operator==( const GlsMatrix<Type, DIM>& rhs ) const;
287  bool operator!=( const GlsMatrix<Type, DIM>& rhs ) const;
288 
289  // Assignment operators
290  GlsMatrix<Type, DIM>& operator=( const GlsMatrix<Type, DIM>& m );
291  GlsMatrix<Type, DIM>& operator=( const CStyleMatrix m );
292 
293  // Unary operators
294  GlsMatrix<Type, DIM> operator+() const { return *this; }
295  GlsMatrix<Type, DIM> operator-() const;
296 
297  // Combined assignment - calculation operators
298  GlsMatrix<Type, DIM>& operator+=( const GlsMatrix<Type, DIM>& rhs );
299  GlsMatrix<Type, DIM>& operator-=( const GlsMatrix<Type, DIM>& rhs );
300  GlsMatrix<Type, DIM>& operator*=( const GlsMatrix<Type, DIM>& rhs );
301  GlsMatrix<Type, DIM>& operator*=( Type scalar );
302  GlsMatrix<Type, DIM>& operator/=( Type scalar );
303 
304  // Mathematical operators
305  GlsMatrix<Type, DIM> operator+( const GlsMatrix<Type, DIM>& rhs ) const;
306  GlsMatrix<Type, DIM> operator-( const GlsMatrix<Type, DIM>& rhs ) const;
307  GlsMatrix<Type, DIM> operator*( const GlsMatrix<Type, DIM>& rhs ) const;
308  GlsMatrix<Type, DIM> operator*( Type scalar ) const;
309  GlsMatrix<Type, DIM> operator/( Type scalar ) const;
310 
311 protected:
312  /** A very tiny number, close to zero. Used to compare floats */
313  static const double PRECISION;
314 
315  /**
316  * This is the actual matrix data stored in column-major ordering in
317  * a one-dimensional array of the appropriate size to store all matrix
318  * terms. The data is stored in this manner so that it can be passed
319  * to OpenGL matrix operations, which expect matrices to be ordered
320  * this way.
321  */
322  Type _data[ DIM * DIM ];
323 
324  /**
325  * A column index array. Each element in this array "points"
326  * to a column of the matrix within the above data. By using this index,
327  * we can address matrix terms using two dimensions, column and row,
328  * without incurring a multiplication. Eg. matrix[0] points to a column
329  * within data, which is the address of the first term in the column.
330  * We can further address into the column and access a particular row
331  * by adding another dimension, eg. matrix[0][1]. The indices must
332  * be set up at construction time.
333  */
334  Type* _matrix[ DIM ];
335 
336  //------------------------------------------------------------------------
337  /**
338  * This should be called from all constructors, as it initializes the
339  * matrix column index so that matrix elements can be addressed in two
340  * dimensions, but a multiplication is not incurred when doing so.
341  */
342  //------------------------------------------------------------------------
343  void Initialize();
344 
345  //------------------------------------------------------------------------
346  /**
347  * Pivoting is used in conjunction with Gauss-Jordan reduction to find
348  * the inverse of a matrix and in Gaussian Elimination to find the
349  * determinant of a matrix. It involves the elementary row operation of
350  * interchanging rows (or columns in our column-major matrix). Only the
351  * row passed in and rows below it are considered. The row passed in
352  * will be swapped with another row below it by selecting the desired
353  * pivot element. The pivot is the largest term (in magnitude) in
354  * the column equal to row for any of the considered rows. This will
355  * always put the largest (in magnitude) term from the column on the
356  * diagonal of the matrix.
357  *
358  * \param row is the current row to pivot about. Rows above this are not
359  * considered during pivoting, as these rows should already be in
360  * their proper form.
361  */
362  //------------------------------------------------------------------------
363  int Pivot( unsigned row );
364 
365  //------------------------------------------------------------------------
366  /**
367  * Decides if the two values are close together determined by the
368  * constant PRECISION. The formula is
369  * lv - PRECISION < rv < lv + PRECISION
370  *
371  * \param lv left value
372  * \param rv right value
373  */
374  //------------------------------------------------------------------------
375  bool closeValues( Type lv, Type rv ) const;
376 };
377 
378 //============================================================================
379 //
380 // CLASS IMPLEMENTATION
381 //
382 //============================================================================
383 
384 //----------------------------------------------------------------------------
385 template<class Type, int DIM>
386 const double GlsMatrix<Type, DIM>::PRECISION = 1e-7;
387 
388 //----------------------------------------------------------------------------
389 template<class Type, int DIM>
391 {
392  Initialize();
393  MakeIdentity();
394 }
395 //----------------------------------------------------------------------------
396 template<class Type, int DIM>
397 GlsMatrix<Type, DIM>::GlsMatrix( const CStyleMatrix m )
398 {
399  Initialize();
400  for( int row = 0; row < DIM; ++row )
401  {
402  for( int col = 0; col < DIM; ++col )
403  {
404  _matrix[ col ][ row ] = m[ row ][ col ];
405  }
406  }
407 }
408 //----------------------------------------------------------------------------
409 template<class Type, int DIM>
411 {
412  Initialize();
413  memcpy( _data, m, sizeof( Type ) * DIM * DIM );
414 }
415 //----------------------------------------------------------------------------
416 template<class Type, int DIM>
418 {
419  Initialize();
420  memcpy( _data, m._data, sizeof( Type ) * DIM * DIM );
421 }
422 //----------------------------------------------------------------------------
423 template<class Type, int DIM>
425 {
426 }
427 //----------------------------------------------------------------------------
428 template<>
431 {
432  GlsMatrix<double, 4> result;
433 
434  // Same as regular, just unroll the loop
435  result._matrix[ 0 ][ 0 ] = _matrix[ 0 ][ 0 ] * rhs._matrix[ 0 ][ 0 ] + _matrix[ 1 ][ 0 ] * rhs._matrix[ 0 ][ 1 ] + _matrix[ 2 ][ 0 ] * rhs._matrix[ 0 ][ 2 ] + _matrix[ 3 ][ 0 ] * rhs._matrix[ 0 ][ 3 ];
436  result._matrix[ 1 ][ 0 ] = _matrix[ 0 ][ 0 ] * rhs._matrix[ 1 ][ 0 ] + _matrix[ 1 ][ 0 ] * rhs._matrix[ 1 ][ 1 ] + _matrix[ 2 ][ 0 ] * rhs._matrix[ 1 ][ 2 ] + _matrix[ 3 ][ 0 ] * rhs._matrix[ 1 ][ 3 ];
437  result._matrix[ 2 ][ 0 ] = _matrix[ 0 ][ 0 ] * rhs._matrix[ 2 ][ 0 ] + _matrix[ 1 ][ 0 ] * rhs._matrix[ 2 ][ 1 ] + _matrix[ 2 ][ 0 ] * rhs._matrix[ 2 ][ 2 ] + _matrix[ 3 ][ 0 ] * rhs._matrix[ 2 ][ 3 ];
438  result._matrix[ 3 ][ 0 ] = _matrix[ 0 ][ 0 ] * rhs._matrix[ 3 ][ 0 ] + _matrix[ 1 ][ 0 ] * rhs._matrix[ 3 ][ 1 ] + _matrix[ 2 ][ 0 ] * rhs._matrix[ 3 ][ 2 ] + _matrix[ 3 ][ 0 ] * rhs._matrix[ 3 ][ 3 ];
439  result._matrix[ 0 ][ 1 ] = _matrix[ 0 ][ 1 ] * rhs._matrix[ 0 ][ 0 ] + _matrix[ 1 ][ 1 ] * rhs._matrix[ 0 ][ 1 ] + _matrix[ 2 ][ 1 ] * rhs._matrix[ 0 ][ 2 ] + _matrix[ 3 ][ 1 ] * rhs._matrix[ 0 ][ 3 ];
440  result._matrix[ 1 ][ 1 ] = _matrix[ 0 ][ 1 ] * rhs._matrix[ 1 ][ 0 ] + _matrix[ 1 ][ 1 ] * rhs._matrix[ 1 ][ 1 ] + _matrix[ 2 ][ 1 ] * rhs._matrix[ 1 ][ 2 ] + _matrix[ 3 ][ 1 ] * rhs._matrix[ 1 ][ 3 ];
441  result._matrix[ 2 ][ 1 ] = _matrix[ 0 ][ 1 ] * rhs._matrix[ 2 ][ 0 ] + _matrix[ 1 ][ 1 ] * rhs._matrix[ 2 ][ 1 ] + _matrix[ 2 ][ 1 ] * rhs._matrix[ 2 ][ 2 ] + _matrix[ 3 ][ 1 ] * rhs._matrix[ 2 ][ 3 ];
442  result._matrix[ 3 ][ 1 ] = _matrix[ 0 ][ 1 ] * rhs._matrix[ 3 ][ 0 ] + _matrix[ 1 ][ 1 ] * rhs._matrix[ 3 ][ 1 ] + _matrix[ 2 ][ 1 ] * rhs._matrix[ 3 ][ 2 ] + _matrix[ 3 ][ 1 ] * rhs._matrix[ 3 ][ 3 ];
443  result._matrix[ 0 ][ 2 ] = _matrix[ 0 ][ 2 ] * rhs._matrix[ 0 ][ 0 ] + _matrix[ 1 ][ 2 ] * rhs._matrix[ 0 ][ 1 ] + _matrix[ 2 ][ 2 ] * rhs._matrix[ 0 ][ 2 ] + _matrix[ 3 ][ 2 ] * rhs._matrix[ 0 ][ 3 ];
444  result._matrix[ 1 ][ 2 ] = _matrix[ 0 ][ 2 ] * rhs._matrix[ 1 ][ 0 ] + _matrix[ 1 ][ 2 ] * rhs._matrix[ 1 ][ 1 ] + _matrix[ 2 ][ 2 ] * rhs._matrix[ 1 ][ 2 ] + _matrix[ 3 ][ 2 ] * rhs._matrix[ 1 ][ 3 ];
445  result._matrix[ 2 ][ 2 ] = _matrix[ 0 ][ 2 ] * rhs._matrix[ 2 ][ 0 ] + _matrix[ 1 ][ 2 ] * rhs._matrix[ 2 ][ 1 ] + _matrix[ 2 ][ 2 ] * rhs._matrix[ 2 ][ 2 ] + _matrix[ 3 ][ 2 ] * rhs._matrix[ 2 ][ 3 ];
446  result._matrix[ 3 ][ 2 ] = _matrix[ 0 ][ 2 ] * rhs._matrix[ 3 ][ 0 ] + _matrix[ 1 ][ 2 ] * rhs._matrix[ 3 ][ 1 ] + _matrix[ 2 ][ 2 ] * rhs._matrix[ 3 ][ 2 ] + _matrix[ 3 ][ 2 ] * rhs._matrix[ 3 ][ 3 ];
447  result._matrix[ 0 ][ 3 ] = _matrix[ 0 ][ 3 ] * rhs._matrix[ 0 ][ 0 ] + _matrix[ 1 ][ 3 ] * rhs._matrix[ 0 ][ 1 ] + _matrix[ 2 ][ 3 ] * rhs._matrix[ 0 ][ 2 ] + _matrix[ 3 ][ 3 ] * rhs._matrix[ 0 ][ 3 ];
448  result._matrix[ 1 ][ 3 ] = _matrix[ 0 ][ 3 ] * rhs._matrix[ 1 ][ 0 ] + _matrix[ 1 ][ 3 ] * rhs._matrix[ 1 ][ 1 ] + _matrix[ 2 ][ 3 ] * rhs._matrix[ 1 ][ 2 ] + _matrix[ 3 ][ 3 ] * rhs._matrix[ 1 ][ 3 ];
449  result._matrix[ 2 ][ 3 ] = _matrix[ 0 ][ 3 ] * rhs._matrix[ 2 ][ 0 ] + _matrix[ 1 ][ 3 ] * rhs._matrix[ 2 ][ 1 ] + _matrix[ 2 ][ 3 ] * rhs._matrix[ 2 ][ 2 ] + _matrix[ 3 ][ 3 ] * rhs._matrix[ 2 ][ 3 ];
450  result._matrix[ 3 ][ 3 ] = _matrix[ 0 ][ 3 ] * rhs._matrix[ 3 ][ 0 ] + _matrix[ 1 ][ 3 ] * rhs._matrix[ 3 ][ 1 ] + _matrix[ 2 ][ 3 ] * rhs._matrix[ 3 ][ 2 ] + _matrix[ 3 ][ 3 ] * rhs._matrix[ 3 ][ 3 ];
451 
452  return result;
453 }
454 //----------------------------------------------------------------------------
455 template<>
456 inline GlsMatrix<float, 4>
457  GlsMatrix<float, 4>::operator*( const GlsMatrix<float, 4>& rhs ) const
458 {
459  GlsMatrix<float, 4> result;
460 
461  // Same as regular, just unroll the loop
462  result._matrix[ 0 ][ 0 ] = _matrix[ 0 ][ 0 ] * rhs._matrix[ 0 ][ 0 ] + _matrix[ 1 ][ 0 ] * rhs._matrix[ 0 ][ 1 ] + _matrix[ 2 ][ 0 ] * rhs._matrix[ 0 ][ 2 ] + _matrix[ 3 ][ 0 ] * rhs._matrix[ 0 ][ 3 ];
463  result._matrix[ 1 ][ 0 ] = _matrix[ 0 ][ 0 ] * rhs._matrix[ 1 ][ 0 ] + _matrix[ 1 ][ 0 ] * rhs._matrix[ 1 ][ 1 ] + _matrix[ 2 ][ 0 ] * rhs._matrix[ 1 ][ 2 ] + _matrix[ 3 ][ 0 ] * rhs._matrix[ 1 ][ 3 ];
464  result._matrix[ 2 ][ 0 ] = _matrix[ 0 ][ 0 ] * rhs._matrix[ 2 ][ 0 ] + _matrix[ 1 ][ 0 ] * rhs._matrix[ 2 ][ 1 ] + _matrix[ 2 ][ 0 ] * rhs._matrix[ 2 ][ 2 ] + _matrix[ 3 ][ 0 ] * rhs._matrix[ 2 ][ 3 ];
465  result._matrix[ 3 ][ 0 ] = _matrix[ 0 ][ 0 ] * rhs._matrix[ 3 ][ 0 ] + _matrix[ 1 ][ 0 ] * rhs._matrix[ 3 ][ 1 ] + _matrix[ 2 ][ 0 ] * rhs._matrix[ 3 ][ 2 ] + _matrix[ 3 ][ 0 ] * rhs._matrix[ 3 ][ 3 ];
466  result._matrix[ 0 ][ 1 ] = _matrix[ 0 ][ 1 ] * rhs._matrix[ 0 ][ 0 ] + _matrix[ 1 ][ 1 ] * rhs._matrix[ 0 ][ 1 ] + _matrix[ 2 ][ 1 ] * rhs._matrix[ 0 ][ 2 ] + _matrix[ 3 ][ 1 ] * rhs._matrix[ 0 ][ 3 ];
467  result._matrix[ 1 ][ 1 ] = _matrix[ 0 ][ 1 ] * rhs._matrix[ 1 ][ 0 ] + _matrix[ 1 ][ 1 ] * rhs._matrix[ 1 ][ 1 ] + _matrix[ 2 ][ 1 ] * rhs._matrix[ 1 ][ 2 ] + _matrix[ 3 ][ 1 ] * rhs._matrix[ 1 ][ 3 ];
468  result._matrix[ 2 ][ 1 ] = _matrix[ 0 ][ 1 ] * rhs._matrix[ 2 ][ 0 ] + _matrix[ 1 ][ 1 ] * rhs._matrix[ 2 ][ 1 ] + _matrix[ 2 ][ 1 ] * rhs._matrix[ 2 ][ 2 ] + _matrix[ 3 ][ 1 ] * rhs._matrix[ 2 ][ 3 ];
469  result._matrix[ 3 ][ 1 ] = _matrix[ 0 ][ 1 ] * rhs._matrix[ 3 ][ 0 ] + _matrix[ 1 ][ 1 ] * rhs._matrix[ 3 ][ 1 ] + _matrix[ 2 ][ 1 ] * rhs._matrix[ 3 ][ 2 ] + _matrix[ 3 ][ 1 ] * rhs._matrix[ 3 ][ 3 ];
470  result._matrix[ 0 ][ 2 ] = _matrix[ 0 ][ 2 ] * rhs._matrix[ 0 ][ 0 ] + _matrix[ 1 ][ 2 ] * rhs._matrix[ 0 ][ 1 ] + _matrix[ 2 ][ 2 ] * rhs._matrix[ 0 ][ 2 ] + _matrix[ 3 ][ 2 ] * rhs._matrix[ 0 ][ 3 ];
471  result._matrix[ 1 ][ 2 ] = _matrix[ 0 ][ 2 ] * rhs._matrix[ 1 ][ 0 ] + _matrix[ 1 ][ 2 ] * rhs._matrix[ 1 ][ 1 ] + _matrix[ 2 ][ 2 ] * rhs._matrix[ 1 ][ 2 ] + _matrix[ 3 ][ 2 ] * rhs._matrix[ 1 ][ 3 ];
472  result._matrix[ 2 ][ 2 ] = _matrix[ 0 ][ 2 ] * rhs._matrix[ 2 ][ 0 ] + _matrix[ 1 ][ 2 ] * rhs._matrix[ 2 ][ 1 ] + _matrix[ 2 ][ 2 ] * rhs._matrix[ 2 ][ 2 ] + _matrix[ 3 ][ 2 ] * rhs._matrix[ 2 ][ 3 ];
473  result._matrix[ 3 ][ 2 ] = _matrix[ 0 ][ 2 ] * rhs._matrix[ 3 ][ 0 ] + _matrix[ 1 ][ 2 ] * rhs._matrix[ 3 ][ 1 ] + _matrix[ 2 ][ 2 ] * rhs._matrix[ 3 ][ 2 ] + _matrix[ 3 ][ 2 ] * rhs._matrix[ 3 ][ 3 ];
474  result._matrix[ 0 ][ 3 ] = _matrix[ 0 ][ 3 ] * rhs._matrix[ 0 ][ 0 ] + _matrix[ 1 ][ 3 ] * rhs._matrix[ 0 ][ 1 ] + _matrix[ 2 ][ 3 ] * rhs._matrix[ 0 ][ 2 ] + _matrix[ 3 ][ 3 ] * rhs._matrix[ 0 ][ 3 ];
475  result._matrix[ 1 ][ 3 ] = _matrix[ 0 ][ 3 ] * rhs._matrix[ 1 ][ 0 ] + _matrix[ 1 ][ 3 ] * rhs._matrix[ 1 ][ 1 ] + _matrix[ 2 ][ 3 ] * rhs._matrix[ 1 ][ 2 ] + _matrix[ 3 ][ 3 ] * rhs._matrix[ 1 ][ 3 ];
476  result._matrix[ 2 ][ 3 ] = _matrix[ 0 ][ 3 ] * rhs._matrix[ 2 ][ 0 ] + _matrix[ 1 ][ 3 ] * rhs._matrix[ 2 ][ 1 ] + _matrix[ 2 ][ 3 ] * rhs._matrix[ 2 ][ 2 ] + _matrix[ 3 ][ 3 ] * rhs._matrix[ 2 ][ 3 ];
477  result._matrix[ 3 ][ 3 ] = _matrix[ 0 ][ 3 ] * rhs._matrix[ 3 ][ 0 ] + _matrix[ 1 ][ 3 ] * rhs._matrix[ 3 ][ 1 ] + _matrix[ 2 ][ 3 ] * rhs._matrix[ 3 ][ 2 ] + _matrix[ 3 ][ 3 ] * rhs._matrix[ 3 ][ 3 ];
478 
479  return result;
480 }
481 //----------------------------------------------------------------------------
482 template<class Type, int DIM>
483 inline GlsMatrix<Type, DIM>
484  GlsMatrix<Type, DIM>::operator*( const GlsMatrix<Type, DIM>& rhs ) const
485 {
486  GlsMatrix<Type, DIM> result;
487  memset( result._data, 0, sizeof( Type ) * DIM * DIM );
488  for( int row = 0; row < DIM; ++row )
489  {
490  for( int col = 0; col < DIM; ++col )
491  {
492  for( int i = 0; i < DIM; ++i )
493  {
494  result._matrix[ col ][ row ] += _matrix[ i ][ row ] * rhs._matrix[ col ][ i ];
495  }
496  }
497  }
498  return result;
499 }
500 //----------------------------------------------------------------------------
501 template<class Type, int DIM>
502 inline GlsMatrix<Type, DIM> GlsMatrix<Type, DIM>::operator*( Type scalar ) const
503 {
504  GlsMatrix<Type, DIM> result;
505  for( int row = 0; row < DIM; ++row )
506  {
507  for( int col = 0; col < DIM; ++col )
508  {
509  result._matrix[ col ][ row ] = _matrix[ col ][ row ] * scalar;
510  }
511  }
512  return result;
513 }
514 //----------------------------------------------------------------------------
515 template<class Type, int DIM>
517 {
518  for( int row = 0; row < DIM; ++row )
519  {
520  for( int col = 0; col < DIM; ++col )
521  {
522  std::cout.width( 11 );
523  std::cout << _matrix[ col ][ row ] << ' ';
524  }
525  std::cout << std::endl;
526  }
527  std::cout << std::endl;
528 }
529 //----------------------------------------------------------------------------
530 template<class Type, int DIM>
532 {
533  Clear();
534  for( int i = 0; i < DIM; ++i )
535  {
536  _matrix[ i ][ i ] = Type( 1 );
537  }
538 }
539 //----------------------------------------------------------------------------
540 template<class Type, int DIM>
542 {
543  // In this algorithm we use gaussian elimination with pivoting to reduce
544  // the matrix to an upper triangular form. The determinant of an
545  // upper triangular matrix is the product of the terms on the diagonal.
546  // To perform the elimination, we will use two types of elementary row
547  // operations on the matrix.
548  // 1) Interchanging rows of the matrix.
549  // 2) Add c times row i of the matrix to row j where i != j.
550  //
551  // The first operation of interchanging rows prevents numerical
552  // instability, and the second allows us to zero out terms below
553  // the diagonal by choosing c appropriately.
554 
555  Type pivotFactor;
556  Type det = Type( 1 );
557 
558  // Make a temporary copy of the matrix so the original is not destroyed
559  // during gaussian elimination and pivoting
560  GlsMatrix<Type, DIM> temp( *this );
561 
562  // Loop through each row in the matrix
563  for( int k = 0; k < DIM; ++k )
564  {
565  // Perform pivoting
566  int index = temp.Pivot( k );
567 
568  // If pivot returns a -1, then a zero exists on the diagonal
569  // so the determinant will be zero.
570  if( index == -1 )
571  return 0;
572  // If pivot returns a number other than zero, rows have been
573  // interchanged, so we need to change the sign of the determinant
574  // due to the theorem that states "If matrix B results from A by
575  // interchanging two rows (columns) of A, then det(B) = - det(A).
576  if( index != 0 )
577  det = -det;
578 
579  // The determinant will be the product of the diagonal terms. The
580  // row that the [k][k] term is in already has the proper form
581  // for an upper triangular matrix, so go ahead and use it in the det.
582  det *= temp._matrix[ k ][ k ];
583 
584  // Since we are using the diagonal element from row k, column k, we
585  // need all other terms below this element, in column k, to be
586  // zero. We do this by an elementary operation (2) above to each
587  // row below this diagonal element.
588  for( int i = k + 1; i < DIM; ++i )
589  {
590  // Choose the pivot factor that would zero out the [i][k]term
591  // when multiplied by the kth row and subtracted from the
592  // ith row.
593  pivotFactor = temp._matrix[ i ][ k ] / temp._matrix[ k ][ k ];
594 
595  // Perform the operation on each term in the row to the right
596  // of the pivot. The terms to the left are assumed to be zero
597  // and we never reference them again so no need to waste time
598  // processing them.
599  for( int j = k + 1; j < DIM; ++j )
600  {
601  temp._matrix[ i ][ j ] -= pivotFactor * temp._matrix[ k ][ j ];
602  }
603  }
604  }
605  return det;
606 }
607 //----------------------------------------------------------------------------
608 template<class Type, int DIM>
610 {
611  // In this algorithm we use Gauss-Jordan reduction with pivoting to
612  // transform the original matrix to the identity matrix, all the while
613  // performing the same operations on what started out as the identity
614  // matrix. What started out as the identity matrix will become the
615  // inverse of the original matrix. This is because A*Inv(A) = I and
616  // performing the same operations on either side of the equation to
617  // reduce A to the identity matrix, the equation becomes I*Inv(A) = Inv(A).
618  // To perform the elimination, we will use two types of elementary row
619  // operations on the matrix.
620  // 1) Interchanging rows of the matrix.
621  // 2) Multiply row i of the matrix by c, where c != 0.
622  // 2) Add c times row i of the matrix to row j where i != j.
623  //
624  // The first operation of interchanging rows prevents numerical
625  // instability. The second allows us to put 1's on the diagonal, by
626  // dividing the row by the value of the term on the diagonal.
627  // and the third to zero out terms that are not on the diagonal by
628  // choosing c appropriately.
629 
630  Type a1_inv, a2;
631 
632  // Make a temporary copy of the matrix so the original is not destroyed
633  // during Gauss-Jordan reduction and pivoting
634  GlsMatrix<Type, DIM> temp( *this );
635 
636  // This starts out as the identity matrix. Every operation we perform
637  // on temp, we will also do to this matrix to maintain the equality in
638  // the equation we are starting out with: A*Inv(A) = I
639  GlsMatrix<Type, DIM> inverse;
640 
641  // Loop through each row in the matrix
642  for( int k = 0; k < DIM; ++k )
643  {
644  // Perform pivoting
645  int index = temp.Pivot( k );
646 
647  // If pivot returns a -1, then the matrix is singular or
648  // non-invertible. TBD: WHAT SHOULD WE DO HERE?
649  if( index == -1 )
650  return inverse;
651 
652  // If pivoting resulted in the swapping of rows, then
653  // we need to do the same thing to the inverse matrix
654  if( index != 0 )
655  {
656  // Swap rows
657  Type rowptr[ DIM ];
658  size_t byteCount = sizeof( Type ) * DIM;
659  memcpy( rowptr, inverse._matrix[ k ], byteCount );
660  memcpy( inverse._matrix[ k ], inverse._matrix[ index ], byteCount );
661  memcpy( inverse._matrix[ index ], rowptr, byteCount );
662  }
663 
664  // Divide the current row by the term on the diagonal so we
665  // get a 1 on the diagonal. Remember, do the same operation
666  // to the inverse matrix.
667  a1_inv = ( (Type)1.0 ) / temp._matrix[ k ][ k ];
668  for( int j = 0; j < DIM; ++j )
669  {
670  temp._matrix[ k ][ j ] *= a1_inv;
671  inverse._matrix[ k ][ j ] *= a1_inv;
672  }
673 
674  // Now, loop through each row, skipping the current
675  // one (k).
676  // Choose the pivot factor (a2) that would zero out the [i][k]term
677  // when multiplied by the kth row and subtracted from the
678  // ith row. Since the kth row has a 1 in the diagonal term,
679  // the factor should be the [i][k] term itself. This loop results
680  // in zeros in the kth column, except for the diagonal term [k][k].
681  for( int i = 0; i < DIM; ++i )
682  {
683  if( i != k )
684  {
685  a2 = temp._matrix[ i ][ k ];
686  for( int j = 0; j < DIM; ++j )
687  {
688  temp._matrix[ i ][ j ] -= a2 * temp._matrix[ k ][ j ];
689  inverse._matrix[ i ][ j ] -= a2 * inverse._matrix[ k ][ j ];
690  }
691  }
692  }
693  }
694  return inverse;
695 }
696 //----------------------------------------------------------------------------
697 template<class Type, int DIM>
699 {
700  *this = Inverse();
701 }
702 //----------------------------------------------------------------------------
703 template<class Type, int DIM>
705 {
706  // Interchange rows with columns
707  GlsMatrix<Type, DIM> result;
708  for( int row = 0; row < DIM; ++row )
709  {
710  for( int col = 0; col < DIM; ++col )
711  {
712  result._matrix[ col ][ row ] = _matrix[ row ][ col ];
713  }
714  }
715  return result;
716 }
717 //----------------------------------------------------------------------------
718 template<class Type, int DIM>
720 {
721  GlsMatrix<Type, DIM> result;
722  result = Transposition();
723  *this = result;
724 }
725 //----------------------------------------------------------------------------
726 template<class Type, int DIM>
728 {
729  // Look for the first term off the diagonal (where row != column) that
730  // is not equal to zero. If we find one, the matrix is not a diagonal.
731  for( int row = 0; row < DIM; ++row )
732  {
733  for( int col = 0; col < DIM; ++col )
734  {
735  if( col != row && !closeValues( _matrix[ col ][ row ], Type( 0 ) ) )
736  return false;
737  }
738  }
739  return true;
740 }
741 //----------------------------------------------------------------------------
742 template<class Type, int DIM>
744 {
745  // First do a quick check to see if it is an un-touched identity matrix
746  // This allows for zero tolerences, but should be much faster than
747  // doing all the tolerance checking every time.
748  static const GlsMatrix referenceIdentity; // Defaults to identity.
749  if( 0 == memcmp( _data, referenceIdentity._data, sizeof( Type ) * DIM * DIM ) )
750  return true;
751 
752  // If the matrix is scalar, it is diagonal and all terms on the diagonal
753  // are equal. Make sure they are equal to 1 to be the identity matrix.
754  if( IsScalar() && closeValues( _matrix[ 0 ][ 0 ], Type( 1 ) ) )
755  return true;
756  else
757  return false;
758 }
759 //----------------------------------------------------------------------------
760 template<class Type, int DIM>
762 {
763  // If the matrix is scalar, it is diagonal and all terms on the diagonal
764  // are equal.
765  if( !IsDiagonal() )
766  return false;
767 
768  Type val = _matrix[ 0 ][ 0 ];
769  for( int i = 0; i < DIM; ++i )
770  {
771  if( !closeValues( _matrix[ i ][ i ], val ) )
772  return false;
773  }
774  return true;
775 }
776 //----------------------------------------------------------------------------
777 template<class Type, int DIM>
779 {
780  return Determinant() == Type( 0 );
781 }
782 //----------------------------------------------------------------------------
783 template<class Type, int DIM>
785 {
786  return Transposition() == *this;
787 }
788 //----------------------------------------------------------------------------
789 template<class Type, int DIM>
791 {
792  return Transposition() == -( *this );
793 }
794 //----------------------------------------------------------------------------
795 template<class Type, int DIM>
797 {
798  // An n x n matrix A = [Aij] is lower triangular if Aij = 0 for i < j.
799  for( int row = 0; row < DIM - 1; ++row )
800  {
801  for( int col = row + 1; col < DIM; ++col )
802  {
803  if( !closeValues( _matrix[ col ][ row ], Type( 0 ) ) )
804  return false;
805  }
806  }
807  return true;
808 }
809 //----------------------------------------------------------------------------
810 template<class Type, int DIM>
812 {
813  // An n x n matrix A = [Aij] is upper triangular if Aij = 0 for i > j.
814  for( int row = 1; row < DIM; ++row )
815  {
816  for( int col = 0; col < row; ++col )
817  {
818  if( !closeValues( _matrix[ col ][ row ], Type( 0 ) ) )
819  return false;
820  }
821  }
822  return true;
823 }
824 //----------------------------------------------------------------------------
825 template<class Type, int DIM>
826 inline Type& GlsMatrix<Type, DIM>::operator()( unsigned row, unsigned col )
827 {
828  // transpose the row and column since the matrix is in column
829  // major order but we want users to be able to use row-major ordering.
830  return _matrix[ col ][ row ];
831 }
832 //----------------------------------------------------------------------------
833 template<class Type, int DIM>
834 inline Type GlsMatrix<Type, DIM>::operator()( unsigned row, unsigned col ) const
835 {
836  // transpose the row and column since the matrix is in column
837  // major order but we want users to be able to use row-major ordering.
838  return _matrix[ col ][ row ];
839 }
840 //----------------------------------------------------------------------------
841 template<class Type, int DIM>
842 inline bool GlsMatrix<Type, DIM>::operator==( const GlsMatrix<Type, DIM>& rhs ) const
843 {
844  static const int SIZE = DIM * DIM;
845  for( int i = 0; i < SIZE; ++i )
846  {
847  if( !closeValues( _data[ i ], rhs._data[ i ] ) )
848  return false;
849  }
850  return true;
851 }
852 //----------------------------------------------------------------------------
853 template<class Type, int DIM>
854 inline bool GlsMatrix<Type, DIM>::operator!=( const GlsMatrix<Type, DIM>& rhs ) const
855 {
856  return !( *this == rhs );
857 }
858 //----------------------------------------------------------------------------
859 template<class Type, int DIM>
860 inline GlsMatrix<Type, DIM>&
861 GlsMatrix<Type, DIM>::operator=( const GlsMatrix<Type, DIM>& m )
862 {
863  if( this == &m )
864  return *this;
865 
866  memcpy( _data, m._data, sizeof( Type ) * DIM * DIM );
867 
868  return *this;
869 }
870 //----------------------------------------------------------------------------
871 template<class Type, int DIM>
872 inline GlsMatrix<Type, DIM>&
873 GlsMatrix<Type, DIM>::operator=( const CStyleMatrix m )
874 {
875  for( int row = 0; row < DIM; ++row )
876  {
877  for( int col = 0; col < DIM; ++col )
878  {
879  _matrix[ col ][ row ] = m[ row ][ col ];
880  }
881  }
882 
883  return *this;
884 }
885 //----------------------------------------------------------------------------
886 template<class Type, int DIM>
887 inline GlsMatrix<Type, DIM> GlsMatrix<Type, DIM>::operator-() const
888 {
889  GlsMatrix<Type, DIM> temp;
890 
891  for( int row = 0; row < DIM; ++row )
892  {
893  for( int col = 0; col < DIM; ++col )
894  {
895  temp._matrix[ col ][ row ] = -_matrix[ col ][ row ];
896  }
897  }
898 
899  return temp;
900 }
901 //----------------------------------------------------------------------------
902 template<class Type, int DIM>
903 inline GlsMatrix<Type, DIM>&
904 GlsMatrix<Type, DIM>::operator+=( const GlsMatrix<Type, DIM>& rhs )
905 {
906  for( int row = 0; row < DIM; ++row )
907  {
908  for( int col = 0; col < DIM; ++col )
909  {
910  _matrix[ col ][ row ] += rhs._matrix[ col ][ row ];
911  }
912  }
913  return *this;
914 }
915 //----------------------------------------------------------------------------
916 template<class Type, int DIM>
917 inline GlsMatrix<Type, DIM>&
918 GlsMatrix<Type, DIM>::operator-=( const GlsMatrix<Type, DIM>& rhs )
919 {
920  for( int row = 0; row < DIM; ++row )
921  {
922  for( int col = 0; col < DIM; ++col )
923  {
924  _matrix[ col ][ row ] -= rhs._matrix[ col ][ row ];
925  }
926  }
927  return *this;
928 }
929 //----------------------------------------------------------------------------
930 template<>
931 inline GlsMatrix<double, 4>&
932 GlsMatrix<double, 4>::operator*=( const GlsMatrix<double, 4>& rhs )
933 {
934  return *this = *this * rhs; // faster because it calls the specialized operator*
935 }
936 //----------------------------------------------------------------------------
937 template<>
938 inline GlsMatrix<float, 4>&
939 GlsMatrix<float, 4>::operator*=( const GlsMatrix<float, 4>& rhs )
940 {
941  return *this = *this * rhs; // faster because it calls the specialized operator*
942 }
943 //----------------------------------------------------------------------------
944 template<class Type, int DIM>
945 inline GlsMatrix<Type, DIM>&
946 GlsMatrix<Type, DIM>::operator*=( const GlsMatrix<Type, DIM>& rhs )
947 {
948  GlsMatrix<Type, DIM> result;
949  memset( result._data, 0, sizeof( Type ) * DIM * DIM );
950  for( int row = 0; row < DIM; ++row )
951  {
952  for( int col = 0; col < DIM; ++col )
953  {
954  for( int i = 0; i < DIM; ++i )
955  {
956  result._matrix[ col ][ row ] += _matrix[ i ][ row ] * rhs._matrix[ col ][ i ];
957  }
958  }
959  }
960  *this = result;
961  return *this;
962 }
963 //----------------------------------------------------------------------------
964 template<class Type, int DIM>
965 inline GlsMatrix<Type, DIM>&
966 GlsMatrix<Type, DIM>::operator*=( Type scalar )
967 {
968  for( int row = 0; row < DIM; ++row )
969  {
970  for( int col = 0; col < DIM; ++col )
971  {
972  _matrix[ col ][ row ] *= scalar;
973  }
974  }
975  return *this;
976 }
977 //----------------------------------------------------------------------------
978 template<class Type, int DIM>
979 inline GlsMatrix<Type, DIM>&
980 GlsMatrix<Type, DIM>::operator/=( Type scalar )
981 {
982  for( int row = 0; row < DIM; ++row )
983  {
984  for( int col = 0; col < DIM; ++col )
985  {
986  _matrix[ col ][ row ] /= scalar;
987  }
988  }
989  return *this;
990 }
991 //----------------------------------------------------------------------------
992 template<class Type, int DIM>
993 inline GlsMatrix<Type, DIM>
994 GlsMatrix<Type, DIM>::operator+( const GlsMatrix<Type, DIM>& rhs ) const
995 {
996  GlsMatrix<Type, DIM> result;
997  for( int row = 0; row < DIM; ++row )
998  {
999  for( int col = 0; col < DIM; ++col )
1000  {
1001  result._matrix[ col ][ row ] = _matrix[ col ][ row ] + rhs._matrix[ col ][ row ];
1002  }
1003  }
1004  return result;
1005 }
1006 //----------------------------------------------------------------------------
1007 template<class Type, int DIM>
1008 inline GlsMatrix<Type, DIM>
1009 GlsMatrix<Type, DIM>::operator-( const GlsMatrix<Type, DIM>& rhs ) const
1010 {
1011  GlsMatrix<Type, DIM> result;
1012  for( int row = 0; row < DIM; ++row )
1013  {
1014  for( int col = 0; col < DIM; ++col )
1015  {
1016  result._matrix[ col ][ row ] = _matrix[ col ][ row ] - rhs._matrix[ col ][ row ];
1017  }
1018  }
1019  return result;
1020 }
1021 //----------------------------------------------------------------------------
1022 template<class Type, int DIM>
1023 inline GlsMatrix<Type, DIM> GlsMatrix<Type, DIM>::operator/( Type scalar ) const
1024 {
1025  GlsMatrix<Type, DIM> result;
1026  for( int row = 0; row < DIM; ++row )
1027  {
1028  for( int col = 0; col < DIM; ++col )
1029  {
1030  result._matrix[ col ][ row ] = _matrix[ col ][ row ] / scalar;
1031  }
1032  }
1033  return result;
1034 }
1035 //----------------------------------------------------------------------------
1036 template<class Type, int DIM>
1037 inline std::ostream& operator<<( std::ostream& ostrm, const GlsMatrix<Type, DIM>& m )
1038 {
1039  for( int row = 0; row < DIM; ++row )
1040  {
1041  for( int col = 0; col < DIM; ++col )
1042  {
1043  ostrm.width( 11 );
1044  ostrm << m( row, col ) << ' ';
1045  }
1046  }
1047  return ostrm;
1048 }
1049 //----------------------------------------------------------------------------
1050 template<class Type, int DIM>
1051 inline std::istream& operator>>( std::istream& istrm, GlsMatrix<Type, DIM>& m )
1052 {
1053  Type x;
1054  for( int row = 0; row < DIM; ++row )
1055  {
1056  for( int col = 0; col < DIM; ++col )
1057  {
1058  istrm >> x;
1059  m( row, col ) = x;
1060  }
1061  }
1062  return istrm;
1063 }
1064 //----------------------------------------------------------------------------
1065 template<class Type, int DIM>
1067 {
1068  // Set each element in the matrix index to point to a column of the
1069  // _matrix data.
1070  for( int i = 0, j = 0; i < DIM; ++i, j += DIM )
1071  {
1072  _matrix[ i ] = &( _data[ j ] );
1073  }
1074 }
1075 //----------------------------------------------------------------------------
1076 template<class Type, int DIM>
1077 int GlsMatrix<Type, DIM>::Pivot( unsigned row )
1078 {
1079  int pivotRow;
1080  int k = int( row );
1081  double amax = -1;
1082  double temp;
1083 
1084  // Go through each element in the column below the diagonal element
1085  // [row][row] to find the largest term (in magnitude). Save which row
1086  // has the largest term in the column.
1087  for( unsigned i = row; i < unsigned( DIM ); ++i )
1088  {
1089  if( ( temp = fabs( _matrix[ i ][ row ] ) ) > amax && temp != 0.0 )
1090  {
1091  amax = temp;
1092  k = i;
1093  }
1094  }
1095 
1096  if( _matrix[ k ][ row ] == Type( 0 ) )
1097  {
1098  // If the largest term is zero, then return a -1 so the caller knows
1099  // we cannot pivot.
1100  pivotRow = -1;
1101  }
1102  else if( k != int( row ) )
1103  {
1104  // If k is not equal to the passed in row, it means another row
1105  // has the largest term in it, so lets interchange rows.
1106  // Return the row that was interchanged to let the caller know.
1107  Type rowptr[ DIM ];
1108  size_t byteCount = sizeof( Type ) * DIM;
1109  memcpy( rowptr, _matrix[ k ], byteCount );
1110  memcpy( _matrix[ k ], _matrix[ row ], byteCount );
1111  memcpy( _matrix[ row ], rowptr, byteCount );
1112  pivotRow = k;
1113  }
1114  else
1115  {
1116  // If k is equal to the passed in row, it means current row
1117  // has the largest term in it, so lets not interchange any rows.
1118  // return 0 to let the caller know.
1119  pivotRow = 0;
1120  }
1121 
1122  return pivotRow;
1123 }
1124 //----------------------------------------------------------------------------
1125 template<class Type, int DIM>
1126 inline bool GlsMatrix<Type, DIM>::closeValues( Type lv, Type rv ) const
1127 {
1128  return ( lv - PRECISION ) <= rv && rv <= ( lv + PRECISION );
1129 }
1130 
1131 } // namespace disti
1132 
1133 #endif
bool IsScalar() const
Definition: gls_matrix.h:761
bool IsDiagonal() const
Definition: gls_matrix.h:727
GlsMatrix< Type, DIM > Transposition() const
Definition: gls_matrix.h:704
Definition: gls_matrix.h:61
bool operator!=(const AttributeName &attr1, const AttributeName &attr2)
Definition: disti_metadata.h:149
bool IsLowerTriangular() const
Definition: gls_matrix.h:796
void Initialize()
Definition: gls_matrix.h:1066
static const double PRECISION
Definition: gls_matrix.h:313
GlsMatrix< Type, DIM > Inverse() const
Definition: gls_matrix.h:609
Type * _matrix[DIM]
Definition: gls_matrix.h:334
bool IsSkewSymmetric() const
Definition: gls_matrix.h:790
GlsMatrix()
Definition: gls_matrix.h:390
Type _data[DIM *DIM]
Definition: gls_matrix.h:322
bool operator==(const AttributeName &attr1, const AttributeName &attr2)
Definition: disti_metadata.h:141
Type Determinant() const
Definition: gls_matrix.h:541
bool IsIdentity() const
Definition: gls_matrix.h:743
void Clear()
Definition: gls_matrix.h:123
Type * Data()
Definition: gls_matrix.h:131
bool closeValues(Type lv, Type rv) const
Definition: gls_matrix.h:1126
const Type * Data() const
Definition: gls_matrix.h:139
void Dump() const
Definition: gls_matrix.h:516
bool IsSingular() const
Definition: gls_matrix.h:778
void Invert()
Definition: gls_matrix.h:698
int Pivot(unsigned row)
Definition: gls_matrix.h:1077
void MakeIdentity()
Definition: gls_matrix.h:531
void Transpose()
Definition: gls_matrix.h:719
bool IsUpperTriangular() const
Definition: gls_matrix.h:811
Definition: bmpimage.h:46
virtual ~GlsMatrix()
Definition: gls_matrix.h:424
bool IsSymmetric() const
Definition: gls_matrix.h:784