| 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 ]; ///< Typedef for a 4x4 matrix stored as a 2D array.
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 /// Provides access to the raw column-major matrix data as a one dimensional
126 /// array in a format suitable for Open GL matrix operations.
127 /// \return A pointer to the underlying matrix array.
128 inline Type* Data() { return _data; }
129
130 /// Provides access to the raw column-major matrix data as a one dimensional
131 /// array in a format suitable for Open GL matrix operations.
132 /// \return A const pointer to the underlying matrix array.
133 inline const Type* Data() const { return _data; }
134
135 //------------------------------------------------------------------------
136 /**
137 * Set this matrix to the identity matrix. The identity matrix is
138 * I = [Aij], where Aii = 1 and Aij = 0 for i != j. Thus, for the identity
139 * matrix, the terms off the main diagonal are all zero, and the terms on
140 * the main diagonal are all one.
141 */
142 //------------------------------------------------------------------------
143 void MakeIdentity();
144
145 //------------------------------------------------------------------------
146 /**
147 * Calculates the determinant of the matrix.
148 *
149 * \return the calculated value of the determinant.
150 */
151 //------------------------------------------------------------------------
152 Type Determinant() const;
153
154 //------------------------------------------------------------------------
155 /**
156 * Calculates the inverse of this matrix. The inverse of matrix A is B if
157 * AB = BA = I, where I is the identity matrix. The inverse of a matrix is
158 * the notion corresponding to the reciprocal of a nonzero real number.
159 *
160 * \return the inverse of this matrix.
161 */
162 //------------------------------------------------------------------------
163 GlsMatrix<Type, DIM> Inverse() const;
164
165 //------------------------------------------------------------------------
166 /**
167 * Sets this matrix equal to its inverse.
168 */
169 //------------------------------------------------------------------------
170 void Invert();
171
172 //------------------------------------------------------------------------
173 /**
174 * Calculates the transpose of this matrix. The transpose of a matrix
175 * is obtained by interchanging its rows and columns.
176 *
177 * \return the transpose of this matrix.
178 */
179 //------------------------------------------------------------------------
180 GlsMatrix<Type, DIM> Transposition() const;
181
182 //------------------------------------------------------------------------
183 /**
184 * Sets this matrix equal to its transpose.
185 */
186 //------------------------------------------------------------------------
187 void Transpose();
188
189 //------------------------------------------------------------------------
190 /**
191 * Determines if this is a diagonal matrix. An n x n matrix A = [Aij] is
192 * called a diagonal matrix if Aij = 0 for i != j. Thus, for a diagonal
193 * matrix, the terms off the main diagonal are all zero.
194 *
195 * \return true if the matrix is a diagonal matrix, false otherwise.
196 */
197 //------------------------------------------------------------------------
198 bool IsDiagonal() const;
199
200 //------------------------------------------------------------------------
201 /**
202 * Determines if this matrix is equal to the identity matrix. The
203 * identity matrix is I = [Aij], where Aii = 1 and Aij = 0 for i != j.
204 * Thus, for the identity matrix, the terms off the main diagonal are all
205 * zero, and the terms on the main diagonal are all one.
206 *
207 * \return true if the matrix is the identity matrix, false otherwise.
208 */
209 //------------------------------------------------------------------------
210 bool IsIdentity() const;
211
212 //------------------------------------------------------------------------
213 /**
214 * Determines if this is a scalar matrix. A scalar matrix is a diagonal
215 * matrix whose diagonal terms are equal.
216 *
217 * \return true if the matrix is a scalar matrix, false otherwise.
218 */
219 //------------------------------------------------------------------------
220 bool IsScalar() const;
221
222 //------------------------------------------------------------------------
223 /**
224 * Determines if this matrix is singular, or non-invertible. A singular
225 * matrix has no inverse.
226 *
227 * \return true if the matrix is singular, false otherwise.
228 */
229 //------------------------------------------------------------------------
230 bool IsSingular() const;
231
232 //------------------------------------------------------------------------
233 /**
234 * Determines if this matrix is symmetric. A matrix is symmetric if it is
235 * equal to its transpose, i.e. A = A'.
236 *
237 * \return true if the matrix is symmetric, false otherwise.
238 */
239 //------------------------------------------------------------------------
240 bool IsSymmetric() const;
241
242 //------------------------------------------------------------------------
243 /**
244 * Determines if this matrix is skew symmetric. A matrix is skew symmetric
245 * if its transpose is equal to its negative, i.e. -A = A'.
246 *
247 * \return true if the matrix is skew symmetric, false otherwise.
248 */
249 //------------------------------------------------------------------------
250 bool IsSkewSymmetric() const;
251
252 //------------------------------------------------------------------------
253 /**
254 * Determines if this matrix is lower triangular. An n x n matrix A = [Aij]
255 * is lower triangular if Aij = 0 for i < j.
256 *
257 * \return true if the matrix is lower triangular, false otherwise.
258 */
259 //------------------------------------------------------------------------
260 bool IsLowerTriangular() const;
261
262 //------------------------------------------------------------------------
263 /**
264 * Determines if this matrix is upper triangular. An n x n matrix A = [Aij]
265 * is upper triangular if Aij = 0 for i > j.
266 *
267 * \return true if the matrix is upper triangular, false otherwise.
268 */
269 //------------------------------------------------------------------------
270 bool IsUpperTriangular() const;
271
272 //------------------------------------------------------------------------
273 // Operators
274 //------------------------------------------------------------------------
275
276 /// Subscript operator (zero-based, no bounds checking).
277 /// \param row The row of the element to retrieve.
278 /// \param col The column of the element to retrieve.
279 /// \return A reference to the element at that position.
280 Type& operator()( unsigned row, unsigned col );
281
282 /// Subscript operator (zero-based, no bounds checking).
283 /// \param row The row of the element to retrieve.
284 /// \param col The column of the element to retrieve.
285 /// \return The value of the element at that position.
286 Type operator()( unsigned row, unsigned col ) const;
287
288 // Equality operator
289 /// \param rhs The matrix to compare.
290 /// \return Whether or not the matrices are equal.
291 bool operator==( const GlsMatrix<Type, DIM>& rhs ) const;
292
293 // Inequality operator
294 /// \param rhs The matrix to compare.
295 /// \return whether or not the matrices are not equal.
296 bool operator!=( const GlsMatrix<Type, DIM>& rhs ) const;
297
298 /// Assignment operator
299 /// \param m The matrix to copy from.
300 /// \return The resulting matrix (this).
301 GlsMatrix<Type, DIM>& operator=( const GlsMatrix<Type, DIM>& m );
302
303 /// Assignment operator
304 /// \param m The matrix to copy from.
305 /// \return The resulting matrix (this).
306 GlsMatrix<Type, DIM>& operator=( const CStyleMatrix m );
307
308 /// Unary + operator
309 /// \return This matrix.
310 GlsMatrix<Type, DIM> operator+() const { return *this; }
311
312 /// Unary - operator
313 /// \return Matrix with all elements multiplied by -1.
314 GlsMatrix<Type, DIM> operator-() const;
315
316 /// Combined assignment += operator
317 /// \param rhs Matrix to add.
318 /// \return The resulting matrix (this).
319 GlsMatrix<Type, DIM>& operator+=( const GlsMatrix<Type, DIM>& rhs );
320
321 /// Combined assignment -= operator
322 /// \param rhs Matrix to subtract.
323 /// \return The resulting matrix (this).
324 GlsMatrix<Type, DIM>& operator-=( const GlsMatrix<Type, DIM>& rhs );
325
326 /// Combined assignment *= operator
327 /// \param rhs Matrix to multiply by.
328 /// \return The resulting matrix (this).
329 GlsMatrix<Type, DIM>& operator*=( const GlsMatrix<Type, DIM>& rhs );
330
331 /// Combined assignment *= operator
332 /// \param scalar Value to multiply by.
333 /// \return The resulting matrix (this).
334 GlsMatrix<Type, DIM>& operator*=( Type scalar );
335
336 /// Combined assignment /= operator
337 /// \param scalar Value to divide by.
338 /// \return The resulting matrix (this).
339 GlsMatrix<Type, DIM>& operator/=( Type scalar );
340
341 // Mathematical + operator
342 /// \param rhs The matrix to add.
343 /// \return The resulting matrix.
344 GlsMatrix<Type, DIM> operator+( const GlsMatrix<Type, DIM>& rhs ) const;
345
346 // Mathematical - operator
347 /// \param rhs The matrix to subtract.
348 /// \return The resulting matrix.
349 GlsMatrix<Type, DIM> operator-( const GlsMatrix<Type, DIM>& rhs ) const;
350
351 // Mathematical * operator
352 /// \param rhs The matrix to multiply by.
353 /// \return The resulting matrix.
354 GlsMatrix<Type, DIM> operator*( const GlsMatrix<Type, DIM>& rhs ) const;
355
356 // Mathematical * operator
357 /// \param scalar The value to multiply by.
358 /// \return The resulting matrix.
359 GlsMatrix<Type, DIM> operator*( Type scalar ) const;
360
361 // Mathematical / operator
362 /// \param scalar The value to divide by.
363 /// \return The resulting matrix.
364 GlsMatrix<Type, DIM> operator/( Type scalar ) const;
365
366 protected:
367 /** A very tiny number, close to zero. Used to compare floats */
368 static const double PRECISION;
369
370 /**
371 * This is the actual matrix data stored in column-major ordering in
372 * a one-dimensional array of the appropriate size to store all matrix
373 * terms. The data is stored in this manner so that it can be passed
374 * to OpenGL matrix operations, which expect matrices to be ordered
375 * this way.
376 */
377 Type _data[ DIM * DIM ];
378
379 /**
380 * A column index array. Each element in this array "points"
381 * to a column of the matrix within the above data. By using this index,
382 * we can address matrix terms using two dimensions, column and row,
383 * without incurring a multiplication. Eg. matrix[0] points to a column
384 * within data, which is the address of the first term in the column.
385 * We can further address into the column and access a particular row
386 * by adding another dimension, eg. matrix[0][1]. The indices must
387 * be set up at construction time.
388 */
389 Type* _matrix[ DIM ];
390
391 //------------------------------------------------------------------------
392 /**
393 * This should be called from all constructors, as it initializes the
394 * matrix column index so that matrix elements can be addressed in two
395 * dimensions, but a multiplication is not incurred when doing so.
396 */
397 //------------------------------------------------------------------------
398 void Initialize();
399
400 /// Pivoting is used in conjunction with Gauss-Jordan reduction to find
401 /// the inverse of a matrix and in Gaussian Elimination to find the
402 /// determinant of a matrix. It involves the elementary row operation of
403 /// interchanging rows (or columns in our column-major matrix). Only the
404 /// row passed in and rows below it are considered. The row passed in
405 /// will be swapped with another row below it by selecting the desired
406 /// pivot element. The pivot is the largest term (in magnitude) in
407 /// the column equal to row for any of the considered rows. This will
408 /// always put the largest (in magnitude) term from the column on the
409 /// diagonal of the matrix.
410 /// \param row is the current row to pivot about. Rows above this are not
411 /// considered during pivoting, as these rows should already be in
412 /// their proper form.
413 /// \return the row number index that was interchanged, if the pivot is different to the incoming row
414 /// 0 if the pivot was successful on the incoming row
415 /// -1 if no pivot is possible
416 int Pivot( unsigned row );
417
418 /// Decides if the two values are close together determined by the constant PRECISION.
419 /// The formula is: lv - PRECISION < rv < lv + PRECISION
420 /// \param lv The left value.
421 /// \param rv The right value.
422 /// \return Whether or not the two values are close together.
423 bool closeValues( Type lv, Type rv ) const;
424};
425
426//============================================================================
427//
428// CLASS IMPLEMENTATION
429//
430//============================================================================
431
432//----------------------------------------------------------------------------
433 template<class Type, int DIM>
434 const double GlsMatrix<Type, DIM>::PRECISION = 1e-7;
435
436//----------------------------------------------------------------------------
437 template<class Type, int DIM>
438 GlsMatrix<Type, DIM>::GlsMatrix()
439{
440 Initialize();
441 MakeIdentity();
442}
443//----------------------------------------------------------------------------
444 template<class Type, int DIM>
445 GlsMatrix<Type, DIM>::GlsMatrix( const CStyleMatrix m )
446{
447 Initialize();
448 for( int row = 0; row < DIM; ++row )
449 {
450 for( int col = 0; col < DIM; ++col )
451 {
452 _matrix[ col ][ row ] = m[ row ][ col ];
453 }
454 }
455}
456//----------------------------------------------------------------------------
457 template<class Type, int DIM>
458 GlsMatrix<Type, DIM>::GlsMatrix( const Type* m )
459{
460 Initialize();
461 memcpy( _data, m, sizeof( Type ) * DIM * DIM );
462}
463//----------------------------------------------------------------------------
464 template<class Type, int DIM>
465 GlsMatrix<Type, DIM>::GlsMatrix( const GlsMatrix<Type, DIM>& m )
466{
467 Initialize();
468 memcpy( _data, m._data, sizeof( Type ) * DIM * DIM );
469}
470//----------------------------------------------------------------------------
471 template<class Type, int DIM>
472 GlsMatrix<Type, DIM>::~GlsMatrix()
473{
474}
475
476/// Mathematical * operator
477/// \param rhs The matrix to multiply by.
478/// \return The resulting matrix.
479 template<>
480 inline GlsMatrix<double, 4>
481 GlsMatrix<double, 4>::operator*( const GlsMatrix<double, 4>& rhs ) const
482{
483 GlsMatrix<double, 4> result;
484
485 // Same as regular, just unroll the loop
486 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 ];
487 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 ];
488 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 ];
489 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 ];
490 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 ];
491 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 ];
492 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 ];
493 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 ];
494 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 ];
495 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 ];
496 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 ];
497 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 ];
498 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 ];
499 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 ];
500 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 ];
501 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 ];
502
503 return result;
504}
505
506/// Mathematical * operator
507/// \param rhs The matrix to multiply by.
508/// \return The resulting matrix.
509 template<>
510 inline GlsMatrix<float, 4>
511 GlsMatrix<float, 4>::operator*( const GlsMatrix<float, 4>& rhs ) const
512{
513 GlsMatrix<float, 4> result;
514
515 // Same as regular, just unroll the loop
516 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 ];
517 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 ];
518 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 ];
519 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 ];
520 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 ];
521 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 ];
522 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 ];
523 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 ];
524 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 ];
525 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 ];
526 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 ];
527 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 ];
528 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 ];
529 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 ];
530 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 ];
531 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 ];
532
533 return result;
534}
535//----------------------------------------------------------------------------
536 template<class Type, int DIM>
537 inline GlsMatrix<Type, DIM>
538 GlsMatrix<Type, DIM>::operator*( const GlsMatrix<Type, DIM>& rhs ) const
539{
540 GlsMatrix<Type, DIM> result;
541 memset( result._data, 0, sizeof( Type ) * DIM * DIM );
542 for( int row = 0; row < DIM; ++row )
543 {
544 for( int col = 0; col < DIM; ++col )
545 {
546 for( int i = 0; i < DIM; ++i )
547 {
548 result._matrix[ col ][ row ] += _matrix[ i ][ row ] * rhs._matrix[ col ][ i ];
549 }
550 }
551 }
552 return result;
553}
554//----------------------------------------------------------------------------
555 template<class Type, int DIM>
556 inline GlsMatrix<Type, DIM> GlsMatrix<Type, DIM>::operator*( Type scalar ) const
557{
558 GlsMatrix<Type, DIM> result;
559 for( int row = 0; row < DIM; ++row )
560 {
561 for( int col = 0; col < DIM; ++col )
562 {
563 result._matrix[ col ][ row ] = _matrix[ col ][ row ] * scalar;
564 }
565 }
566 return result;
567}
568//----------------------------------------------------------------------------
569 template<class Type, int DIM>
570 void GlsMatrix<Type, DIM>::Dump() const
571{
572 for( int row = 0; row < DIM; ++row )
573 {
574 for( int col = 0; col < DIM; ++col )
575 {
576 std::cout.width( 11 );
577 std::cout << _matrix[ col ][ row ] << ' ';
578 }
579 std::cout << std::endl;
580 }
581 std::cout << std::endl;
582}
583//----------------------------------------------------------------------------
584 template<class Type, int DIM>
585 inline void GlsMatrix<Type, DIM>::MakeIdentity()
586{
587 Clear();
588 for( int i = 0; i < DIM; ++i )
589 {
590 _matrix[ i ][ i ] = Type( 1 );
591 }
592}
593//----------------------------------------------------------------------------
594 template<class Type, int DIM>
595Type GlsMatrix<Type, DIM>::Determinant() const
596{
597 // In this algorithm we use gaussian elimination with pivoting to reduce
598 // the matrix to an upper triangular form. The determinant of an
599 // upper triangular matrix is the product of the terms on the diagonal.
600 // To perform the elimination, we will use two types of elementary row
601 // operations on the matrix.
602 // 1) Interchanging rows of the matrix.
603 // 2) Add c times row i of the matrix to row j where i != j.
604 //
605 // The first operation of interchanging rows prevents numerical
606 // instability, and the second allows us to zero out terms below
607 // the diagonal by choosing c appropriately.
608
609 Type pivotFactor;
610 Type det = Type( 1 );
611
612 // Make a temporary copy of the matrix so the original is not destroyed
613 // during gaussian elimination and pivoting
614 GlsMatrix<Type, DIM> temp( *this );
615
616 // Loop through each row in the matrix
617 for( int k = 0; k < DIM; ++k )
618 {
619 // Perform pivoting
620 int index = temp.Pivot( k );
621
622 // If pivot returns a -1, then a zero exists on the diagonal
623 // so the determinant will be zero.
624 if( index == -1 )
625 return 0;
626 // If pivot returns a number other than zero, rows have been
627 // interchanged, so we need to change the sign of the determinant
628 // due to the theorem that states "If matrix B results from A by
629 // interchanging two rows (columns) of A, then det(B) = - det(A).
630 if( index != 0 )
631 det = -det;
632
633 // The determinant will be the product of the diagonal terms. The
634 // row that the [k][k] term is in already has the proper form
635 // for an upper triangular matrix, so go ahead and use it in the det.
636 det *= temp._matrix[ k ][ k ];
637
638 // Since we are using the diagonal element from row k, column k, we
639 // need all other terms below this element, in column k, to be
640 // zero. We do this by an elementary operation (2) above to each
641 // row below this diagonal element.
642 for( int i = k + 1; i < DIM; ++i )
643 {
644 // Choose the pivot factor that would zero out the [i][k]term
645 // when multiplied by the kth row and subtracted from the
646 // ith row.
647 pivotFactor = temp._matrix[ i ][ k ] / temp._matrix[ k ][ k ];
648
649 // Perform the operation on each term in the row to the right
650 // of the pivot. The terms to the left are assumed to be zero
651 // and we never reference them again so no need to waste time
652 // processing them.
653 for( int j = k + 1; j < DIM; ++j )
654 {
655 temp._matrix[ i ][ j ] -= pivotFactor * temp._matrix[ k ][ j ];
656 }
657 }
658 }
659 return det;
660}
661//----------------------------------------------------------------------------
662 template<class Type, int DIM>
663 GlsMatrix<Type, DIM> GlsMatrix<Type, DIM>::Inverse() const
664{
665 // In this algorithm we use Gauss-Jordan reduction with pivoting to
666 // transform the original matrix to the identity matrix, all the while
667 // performing the same operations on what started out as the identity
668 // matrix. What started out as the identity matrix will become the
669 // inverse of the original matrix. This is because A*Inv(A) = I and
670 // performing the same operations on either side of the equation to
671 // reduce A to the identity matrix, the equation becomes I*Inv(A) = Inv(A).
672 // To perform the elimination, we will use two types of elementary row
673 // operations on the matrix.
674 // 1) Interchanging rows of the matrix.
675 // 2) Multiply row i of the matrix by c, where c != 0.
676 // 2) Add c times row i of the matrix to row j where i != j.
677 //
678 // The first operation of interchanging rows prevents numerical
679 // instability. The second allows us to put 1's on the diagonal, by
680 // dividing the row by the value of the term on the diagonal.
681 // and the third to zero out terms that are not on the diagonal by
682 // choosing c appropriately.
683
684 Type a1_inv, a2;
685
686 // Make a temporary copy of the matrix so the original is not destroyed
687 // during Gauss-Jordan reduction and pivoting
688 GlsMatrix<Type, DIM> temp( *this );
689
690 // This starts out as the identity matrix. Every operation we perform
691 // on temp, we will also do to this matrix to maintain the equality in
692 // the equation we are starting out with: A*Inv(A) = I
693 GlsMatrix<Type, DIM> inverse;
694
695 // Loop through each row in the matrix
696 for( int k = 0; k < DIM; ++k )
697 {
698 // Perform pivoting
699 int index = temp.Pivot( k );
700
701 // If pivot returns a -1, then the matrix is singular or
702 // non-invertible. TBD: WHAT SHOULD WE DO HERE?
703 if( index == -1 )
704 return inverse;
705
706 // If pivoting resulted in the swapping of rows, then
707 // we need to do the same thing to the inverse matrix
708 if( index != 0 )
709 {
710 // Swap rows
711 Type rowptr[ DIM ];
712 size_t byteCount = sizeof( Type ) * DIM;
713 memcpy( rowptr, inverse._matrix[ k ], byteCount );
714 memcpy( inverse._matrix[ k ], inverse._matrix[ index ], byteCount );
715 memcpy( inverse._matrix[ index ], rowptr, byteCount );
716 }
717
718 // Divide the current row by the term on the diagonal so we
719 // get a 1 on the diagonal. Remember, do the same operation
720 // to the inverse matrix.
721 a1_inv = ( (Type)1.0 ) / temp._matrix[ k ][ k ];
722 for( int j = 0; j < DIM; ++j )
723 {
724 temp._matrix[ k ][ j ] *= a1_inv;
725 inverse._matrix[ k ][ j ] *= a1_inv;
726 }
727
728 // Now, loop through each row, skipping the current
729 // one (k).
730 // Choose the pivot factor (a2) that would zero out the [i][k]term
731 // when multiplied by the kth row and subtracted from the
732 // ith row. Since the kth row has a 1 in the diagonal term,
733 // the factor should be the [i][k] term itself. This loop results
734 // in zeros in the kth column, except for the diagonal term [k][k].
735 for( int i = 0; i < DIM; ++i )
736 {
737 if( i != k )
738 {
739 a2 = temp._matrix[ i ][ k ];
740 for( int j = 0; j < DIM; ++j )
741 {
742 temp._matrix[ i ][ j ] -= a2 * temp._matrix[ k ][ j ];
743 inverse._matrix[ i ][ j ] -= a2 * inverse._matrix[ k ][ j ];
744 }
745 }
746 }
747 }
748 return inverse;
749}
750//----------------------------------------------------------------------------
751 template<class Type, int DIM>
752 inline void GlsMatrix<Type, DIM>::Invert()
753{
754 *this = Inverse();
755}
756//----------------------------------------------------------------------------
757 template<class Type, int DIM>
758 inline GlsMatrix<Type, DIM> GlsMatrix<Type, DIM>::Transposition() const
759{
760 // Interchange rows with columns
761 GlsMatrix<Type, DIM> result;
762 for( int row = 0; row < DIM; ++row )
763 {
764 for( int col = 0; col < DIM; ++col )
765 {
766 result._matrix[ col ][ row ] = _matrix[ row ][ col ];
767 }
768 }
769 return result;
770}
771//----------------------------------------------------------------------------
772 template<class Type, int DIM>
773 inline void GlsMatrix<Type, DIM>::Transpose()
774{
775 GlsMatrix<Type, DIM> result;
776 result = Transposition();
777 *this = result;
778}
779//----------------------------------------------------------------------------
780 template<class Type, int DIM>
781 inline bool GlsMatrix<Type, DIM>::IsDiagonal() const
782{
783 // Look for the first term off the diagonal (where row != column) that
784 // is not equal to zero. If we find one, the matrix is not a diagonal.
785 for( int row = 0; row < DIM; ++row )
786 {
787 for( int col = 0; col < DIM; ++col )
788 {
789 if( col != row && !closeValues( _matrix[ col ][ row ], Type( 0 ) ) )
790 return false;
791 }
792 }
793 return true;
794}
795//----------------------------------------------------------------------------
796 template<class Type, int DIM>
797 inline bool GlsMatrix<Type, DIM>::IsIdentity() const
798{
799 // First do a quick check to see if it is an un-touched identity matrix
800 // This allows for zero tolerences, but should be much faster than
801 // doing all the tolerance checking every time.
802 static const GlsMatrix referenceIdentity; // Defaults to identity.
803 if( 0 == memcmp( _data, referenceIdentity._data, sizeof( Type ) * DIM * DIM ) )
804 return true;
805
806 // If the matrix is scalar, it is diagonal and all terms on the diagonal
807 // are equal. Make sure they are equal to 1 to be the identity matrix.
808 if( IsScalar() && closeValues( _matrix[ 0 ][ 0 ], Type( 1 ) ) )
809 return true;
810 else
811 return false;
812}
813//----------------------------------------------------------------------------
814 template<class Type, int DIM>
815 inline bool GlsMatrix<Type, DIM>::IsScalar() const
816{
817 // If the matrix is scalar, it is diagonal and all terms on the diagonal
818 // are equal.
819 if( !IsDiagonal() )
820 return false;
821
822 Type val = _matrix[ 0 ][ 0 ];
823 for( int i = 0; i < DIM; ++i )
824 {
825 if( !closeValues( _matrix[ i ][ i ], val ) )
826 return false;
827 }
828 return true;
829}
830//----------------------------------------------------------------------------
831 template<class Type, int DIM>
832 inline bool GlsMatrix<Type, DIM>::IsSingular() const
833{
834 return Determinant() == Type( 0 );
835}
836//----------------------------------------------------------------------------
837 template<class Type, int DIM>
838 inline bool GlsMatrix<Type, DIM>::IsSymmetric() const
839{
840 return Transposition() == *this;
841}
842//----------------------------------------------------------------------------
843 template<class Type, int DIM>
844 inline bool GlsMatrix<Type, DIM>::IsSkewSymmetric() const
845{
846 return Transposition() == -( *this );
847}
848//----------------------------------------------------------------------------
849 template<class Type, int DIM>
850 inline bool GlsMatrix<Type, DIM>::IsLowerTriangular() const
851{
852 // An n x n matrix A = [Aij] is lower triangular if Aij = 0 for i < j.
853 for( int row = 0; row < DIM - 1; ++row )
854 {
855 for( int col = row + 1; col < DIM; ++col )
856 {
857 if( !closeValues( _matrix[ col ][ row ], Type( 0 ) ) )
858 return false;
859 }
860 }
861 return true;
862}
863//----------------------------------------------------------------------------
864 template<class Type, int DIM>
865 inline bool GlsMatrix<Type, DIM>::IsUpperTriangular() const
866{
867 // An n x n matrix A = [Aij] is upper triangular if Aij = 0 for i > j.
868 for( int row = 1; row < DIM; ++row )
869 {
870 for( int col = 0; col < row; ++col )
871 {
872 if( !closeValues( _matrix[ col ][ row ], Type( 0 ) ) )
873 return false;
874 }
875 }
876 return true;
877}
878//----------------------------------------------------------------------------
879 template<class Type, int DIM>
880 inline Type& GlsMatrix<Type, DIM>::operator()( unsigned row, unsigned col )
881{
882 // transpose the row and column since the matrix is in column
883 // major order but we want users to be able to use row-major ordering.
884 return _matrix[ col ][ row ];
885}
886//----------------------------------------------------------------------------
887 template<class Type, int DIM>
888 inline Type GlsMatrix<Type, DIM>::operator()( unsigned row, unsigned col ) const
889{
890 // transpose the row and column since the matrix is in column
891 // major order but we want users to be able to use row-major ordering.
892 return _matrix[ col ][ row ];
893}
894//----------------------------------------------------------------------------
895 template<class Type, int DIM>
896 inline bool GlsMatrix<Type, DIM>::operator==( const GlsMatrix<Type, DIM>& rhs ) const
897{
898 static const int SIZE = DIM * DIM;
899 for( int i = 0; i < SIZE; ++i )
900 {
901 if( !closeValues( _data[ i ], rhs._data[ i ] ) )
902 return false;
903 }
904 return true;
905}
906//----------------------------------------------------------------------------
907 template<class Type, int DIM>
908 inline bool GlsMatrix<Type, DIM>::operator!=( const GlsMatrix<Type, DIM>& rhs ) const
909{
910 return !( *this == rhs );
911}
912//----------------------------------------------------------------------------
913 template<class Type, int DIM>
914 inline GlsMatrix<Type, DIM>&
915 GlsMatrix<Type, DIM>::operator=( const GlsMatrix<Type, DIM>& m )
916{
917 if( this == &m )
918 return *this;
919
920 memcpy( _data, m._data, sizeof( Type ) * DIM * DIM );
921
922 return *this;
923}
924//----------------------------------------------------------------------------
925 template<class Type, int DIM>
926 inline GlsMatrix<Type, DIM>&
927 GlsMatrix<Type, DIM>::operator=( const CStyleMatrix m )
928{
929 for( int row = 0; row < DIM; ++row )
930 {
931 for( int col = 0; col < DIM; ++col )
932 {
933 _matrix[ col ][ row ] = m[ row ][ col ];
934 }
935 }
936
937 return *this;
938}
939//----------------------------------------------------------------------------
940 template<class Type, int DIM>
941 inline GlsMatrix<Type, DIM> GlsMatrix<Type, DIM>::operator-() const
942{
943 GlsMatrix<Type, DIM> temp;
944
945 for( int row = 0; row < DIM; ++row )
946 {
947 for( int col = 0; col < DIM; ++col )
948 {
949 temp._matrix[ col ][ row ] = -_matrix[ col ][ row ];
950 }
951 }
952
953 return temp;
954}
955//----------------------------------------------------------------------------
956 template<class Type, int DIM>
957 inline GlsMatrix<Type, DIM>&
958 GlsMatrix<Type, DIM>::operator+=( const GlsMatrix<Type, DIM>& rhs )
959{
960 for( int row = 0; row < DIM; ++row )
961 {
962 for( int col = 0; col < DIM; ++col )
963 {
964 _matrix[ col ][ row ] += rhs._matrix[ col ][ row ];
965 }
966 }
967 return *this;
968}
969//----------------------------------------------------------------------------
970 template<class Type, int DIM>
971 inline GlsMatrix<Type, DIM>&
972 GlsMatrix<Type, DIM>::operator-=( const GlsMatrix<Type, DIM>& rhs )
973{
974 for( int row = 0; row < DIM; ++row )
975 {
976 for( int col = 0; col < DIM; ++col )
977 {
978 _matrix[ col ][ row ] -= rhs._matrix[ col ][ row ];
979 }
980 }
981 return *this;
982}
983
984/// Combined assignment *= operator
985/// \param rhs The matrix to multiply by.
986/// \return The resulting matrix (this).
987 template<>
988 inline GlsMatrix<double, 4>&
989 GlsMatrix<double, 4>::operator*=( const GlsMatrix<double, 4>& rhs )
990{
991 return *this = *this * rhs; // faster because it calls the specialized operator*
992}
993
994/// Combined assignment *= operator
995/// \param rhs The matrix to multiply by.
996/// \return The resulting matrix (this).
997 template<>
998 inline GlsMatrix<float, 4>&
999 GlsMatrix<float, 4>::operator*=( const GlsMatrix<float, 4>& rhs )
1000{
1001 return *this = *this * rhs; // faster because it calls the specialized operator*
1002}
1003//----------------------------------------------------------------------------
1004 template<class Type, int DIM>
1005 inline GlsMatrix<Type, DIM>&
1006 GlsMatrix<Type, DIM>::operator*=( const GlsMatrix<Type, DIM>& rhs )
1007{
1008 GlsMatrix<Type, DIM> result;
1009 memset( result._data, 0, sizeof( Type ) * DIM * DIM );
1010 for( int row = 0; row < DIM; ++row )
1011 {
1012 for( int col = 0; col < DIM; ++col )
1013 {
1014 for( int i = 0; i < DIM; ++i )
1015 {
1016 result._matrix[ col ][ row ] += _matrix[ i ][ row ] * rhs._matrix[ col ][ i ];
1017 }
1018 }
1019 }
1020 *this = result;
1021 return *this;
1022}
1023//----------------------------------------------------------------------------
1024 template<class Type, int DIM>
1025 inline GlsMatrix<Type, DIM>&
1026 GlsMatrix<Type, DIM>::operator*=( Type scalar )
1027{
1028 for( int row = 0; row < DIM; ++row )
1029 {
1030 for( int col = 0; col < DIM; ++col )
1031 {
1032 _matrix[ col ][ row ] *= scalar;
1033 }
1034 }
1035 return *this;
1036}
1037//----------------------------------------------------------------------------
1038 template<class Type, int DIM>
1039 inline GlsMatrix<Type, DIM>&
1040 GlsMatrix<Type, DIM>::operator/=( Type scalar )
1041{
1042 for( int row = 0; row < DIM; ++row )
1043 {
1044 for( int col = 0; col < DIM; ++col )
1045 {
1046 _matrix[ col ][ row ] /= scalar;
1047 }
1048 }
1049 return *this;
1050}
1051//----------------------------------------------------------------------------
1052 template<class Type, int DIM>
1053 inline GlsMatrix<Type, DIM>
1054 GlsMatrix<Type, DIM>::operator+( const GlsMatrix<Type, DIM>& rhs ) const
1055{
1056 GlsMatrix<Type, DIM> result;
1057 for( int row = 0; row < DIM; ++row )
1058 {
1059 for( int col = 0; col < DIM; ++col )
1060 {
1061 result._matrix[ col ][ row ] = _matrix[ col ][ row ] + rhs._matrix[ col ][ row ];
1062 }
1063 }
1064 return result;
1065}
1066//----------------------------------------------------------------------------
1067 template<class Type, int DIM>
1068 inline GlsMatrix<Type, DIM>
1069 GlsMatrix<Type, DIM>::operator-( const GlsMatrix<Type, DIM>& rhs ) const
1070{
1071 GlsMatrix<Type, DIM> result;
1072 for( int row = 0; row < DIM; ++row )
1073 {
1074 for( int col = 0; col < DIM; ++col )
1075 {
1076 result._matrix[ col ][ row ] = _matrix[ col ][ row ] - rhs._matrix[ col ][ row ];
1077 }
1078 }
1079 return result;
1080}
1081//----------------------------------------------------------------------------
1082 template<class Type, int DIM>
1083 inline GlsMatrix<Type, DIM> GlsMatrix<Type, DIM>::operator/( Type scalar ) const
1084{
1085 GlsMatrix<Type, DIM> result;
1086 for( int row = 0; row < DIM; ++row )
1087 {
1088 for( int col = 0; col < DIM; ++col )
1089 {
1090 result._matrix[ col ][ row ] = _matrix[ col ][ row ] / scalar;
1091 }
1092 }
1093 return result;
1094}
1095
1096/// Stream out operator
1097/// \param ostrm The stream to write to.
1098/// \param m The value to write to the stream.
1099/// \return The stream in its current state.
1100 template<class Type, int DIM>
1101 inline std::ostream& operator<<( std::ostream& ostrm, const GlsMatrix<Type, DIM>& m )
1102{
1103 for( int row = 0; row < DIM; ++row )
1104 {
1105 for( int col = 0; col < DIM; ++col )
1106 {
1107 ostrm.width( 11 );
1108 ostrm << m( row, col ) << ' ';
1109 }
1110 }
1111 return ostrm;
1112}
1113
1114/// Stream in operator
1115/// \param istrm The stream to read from.
1116/// \param m The returned value read from the stream.
1117/// \return The stream in its current state.
1118 template<class Type, int DIM>
1119 inline std::istream& operator>>( std::istream& istrm, GlsMatrix<Type, DIM>& m )
1120{
1121 Type x;
1122 for( int row = 0; row < DIM; ++row )
1123 {
1124 for( int col = 0; col < DIM; ++col )
1125 {
1126 istrm >> x;
1127 m( row, col ) = x;
1128 }
1129 }
1130 return istrm;
1131}
1132//----------------------------------------------------------------------------
1133 template<class Type, int DIM>
1134 inline void GlsMatrix<Type, DIM>::Initialize()
1135{
1136 // Set each element in the matrix index to point to a column of the
1137 // _matrix data.
1138 for( int i = 0, j = 0; i < DIM; ++i, j += DIM )
1139 {
1140 _matrix[ i ] = &( _data[ j ] );
1141 }
1142}
1143//----------------------------------------------------------------------------
1144 template<class Type, int DIM>
1145 int GlsMatrix<Type, DIM>::Pivot( unsigned row )
1146{
1147 int pivotRow;
1148 int k = int( row );
1149 double amax = -1;
1150 double temp;
1151
1152 // Go through each element in the column below the diagonal element
1153 // [row][row] to find the largest term (in magnitude). Save which row
1154 // has the largest term in the column.
1155 for( unsigned i = row; i < unsigned( DIM ); ++i )
1156 {
1157 if( ( temp = fabs( _matrix[ i ][ row ] ) ) > amax && temp != 0.0 )
1158 {
1159 amax = temp;
1160 k = i;
1161 }
1162 }
1163
1164 if( _matrix[ k ][ row ] == Type( 0 ) )
1165 {
1166 // If the largest term is zero, then return a -1 so the caller knows
1167 // we cannot pivot.
1168 pivotRow = -1;
1169 }
1170 else if( k != int( row ) )
1171 {
1172 // If k is not equal to the passed in row, it means another row
1173 // has the largest term in it, so lets interchange rows.
1174 // Return the row that was interchanged to let the caller know.
1175 Type rowptr[ DIM ];
1176 size_t byteCount = sizeof( Type ) * DIM;
1177 memcpy( rowptr, _matrix[ k ], byteCount );
1178 memcpy( _matrix[ k ], _matrix[ row ], byteCount );
1179 memcpy( _matrix[ row ], rowptr, byteCount );
1180 pivotRow = k;
1181 }
1182 else
1183 {
1184 // If k is equal to the passed in row, it means current row
1185 // has the largest term in it, so lets not interchange any rows.
1186 // return 0 to let the caller know.
1187 pivotRow = 0;
1188 }
1189
1190 return pivotRow;
1191}
1192//----------------------------------------------------------------------------
1193 template<class Type, int DIM>
1194 inline bool GlsMatrix<Type, DIM>::closeValues( Type lv, Type rv ) const
1195{
1196 return ( lv - PRECISION ) <= rv && rv <= ( lv + PRECISION );
1197}
1198
1199} // namespace disti
1200
1201#endif
disti::GlsMatrix
Definition: gls_matrix.h:62
disti::GlsMatrix::IsSkewSymmetric
bool IsSkewSymmetric() const
Definition: gls_matrix.h:844
disti::GlsMatrix::operator!=
bool operator!=(const GlsMatrix< Type, DIM > &rhs) const
Definition: gls_matrix.h:908
disti::GlsMatrix::CStyleMatrix
Type CStyleMatrix[DIM][DIM]
Typedef for a 4x4 matrix stored as a 2D array.
Definition: gls_matrix.h:64
disti::GlsMatrix::IsIdentity
bool IsIdentity() const
Definition: gls_matrix.h:797
disti::GlsMatrix::Transpose
void Transpose()
Definition: gls_matrix.h:773
disti::GlsMatrix::operator+
GlsMatrix< Type, DIM > operator+(const GlsMatrix< Type, DIM > &rhs) const
Definition: gls_matrix.h:1054
disti::GlsMatrix::operator+=
GlsMatrix< Type, DIM > & operator+=(const GlsMatrix< Type, DIM > &rhs)
Definition: gls_matrix.h:958
disti::GlsMatrix::operator()
Type & operator()(unsigned row, unsigned col)
Definition: gls_matrix.h:880
disti::GlsMatrix::operator-
GlsMatrix< Type, DIM > operator-() const
Definition: gls_matrix.h:941
disti::GlsMatrix::operator+
GlsMatrix< Type, DIM > operator+() const
Definition: gls_matrix.h:310
disti::GlsMatrix::Inverse
GlsMatrix< Type, DIM > Inverse() const
Definition: gls_matrix.h:663
disti::GlsMatrix::GlsMatrix
GlsMatrix()
Definition: gls_matrix.h:438
disti::GlsMatrix::~GlsMatrix
virtual ~GlsMatrix()
Definition: gls_matrix.h:472
disti::GlsMatrix::Determinant
Type Determinant() const
Definition: gls_matrix.h:595
disti::GlsMatrix::IsLowerTriangular
bool IsLowerTriangular() const
Definition: gls_matrix.h:850
disti::GlsMatrix::operator*=
GlsMatrix< Type, DIM > & operator*=(const GlsMatrix< Type, DIM > &rhs)
Definition: gls_matrix.h:1006
disti::GlsMatrix::operator==
bool operator==(const GlsMatrix< Type, DIM > &rhs) const
Definition: gls_matrix.h:896
disti::GlsMatrix::MakeIdentity
void MakeIdentity()
Definition: gls_matrix.h:585
disti::GlsMatrix::operator/
GlsMatrix< Type, DIM > operator/(Type scalar) const
Definition: gls_matrix.h:1083
disti::GlsMatrix::operator*=
GlsMatrix< Type, DIM > & operator*=(Type scalar)
Definition: gls_matrix.h:1026
disti::GlsMatrix::Data
Type * Data()
Definition: gls_matrix.h:128
disti::GlsMatrix::IsDiagonal
bool IsDiagonal() const
Definition: gls_matrix.h:781
disti::GlsMatrix::Pivot
int Pivot(unsigned row)
Definition: gls_matrix.h:1145
disti::GlsMatrix::operator=
GlsMatrix< Type, DIM > & operator=(const GlsMatrix< Type, DIM > &m)
Definition: gls_matrix.h:915
disti::GlsMatrix::Data
const Type * Data() const
Definition: gls_matrix.h:133
disti::GlsMatrix::operator=
GlsMatrix< Type, DIM > & operator=(const CStyleMatrix m)
Definition: gls_matrix.h:927
disti::GlsMatrix::Transposition
GlsMatrix< Type, DIM > Transposition() const
Definition: gls_matrix.h:758
disti::GlsMatrix::IsSingular
bool IsSingular() const
Definition: gls_matrix.h:832
disti::GlsMatrix::Dump
void Dump() const
Definition: gls_matrix.h:570
disti::GlsMatrix::IsUpperTriangular
bool IsUpperTriangular() const
Definition: gls_matrix.h:865
disti::GlsMatrix::PRECISION
static const double PRECISION
Definition: gls_matrix.h:368
disti::GlsMatrix::closeValues
bool closeValues(Type lv, Type rv) const
Definition: gls_matrix.h:1194
disti::GlsMatrix::Initialize
void Initialize()
Definition: gls_matrix.h:1134
disti::GlsMatrix::_matrix
Type * _matrix[DIM]
Definition: gls_matrix.h:389
disti::GlsMatrix::Clear
void Clear()
Definition: gls_matrix.h:123
disti::GlsMatrix::operator*
GlsMatrix< Type, DIM > operator*(const GlsMatrix< Type, DIM > &rhs) const
Definition: gls_matrix.h:538
disti::GlsMatrix::operator()
Type operator()(unsigned row, unsigned col) const
Definition: gls_matrix.h:888
disti::GlsMatrix::operator/=
GlsMatrix< Type, DIM > & operator/=(Type scalar)
Definition: gls_matrix.h:1040
disti::GlsMatrix::operator-
GlsMatrix< Type, DIM > operator-(const GlsMatrix< Type, DIM > &rhs) const
Definition: gls_matrix.h:1069
disti::GlsMatrix::GlsMatrix
GlsMatrix(const CStyleMatrix m)
Definition: gls_matrix.h:445
disti::GlsMatrix::Invert
void Invert()
Definition: gls_matrix.h:752
disti::GlsMatrix::GlsMatrix
GlsMatrix(const GlsMatrix< Type, DIM > &m)
Definition: gls_matrix.h:465
disti::GlsMatrix::_data
Type _data[DIM *DIM]
Definition: gls_matrix.h:377
disti::GlsMatrix::operator*
GlsMatrix< Type, DIM > operator*(Type scalar) const
Definition: gls_matrix.h:556
disti::GlsMatrix::operator-=
GlsMatrix< Type, DIM > & operator-=(const GlsMatrix< Type, DIM > &rhs)
Definition: gls_matrix.h:972
disti::GlsMatrix::GlsMatrix
GlsMatrix(const Type *m)
Definition: gls_matrix.h:458
disti::GlsMatrix::IsSymmetric
bool IsSymmetric() const
Definition: gls_matrix.h:838
disti::GlsMatrix::IsScalar
bool IsScalar() const
Definition: gls_matrix.h:815
disti
Definition: bmpimage.h:47
disti::operator!=
bool operator!=(const AttributeName &attr1, const AttributeName &attr2)
Definition: disti_metadata.h:165
disti::operator>>
std::istream & operator>>(std::istream &istrm, GlsMatrix< Type, DIM > &m)
Definition: gls_matrix.h:1119
disti::operator<<
std::ostream & operator<<(std::ostream &ostrm, const GlsMatrix< Type, DIM > &m)
Definition: gls_matrix.h:1101
disti::operator+
GlsPropString operator+(const GlsPropString &a, GlsPropString &b)
Definition: util.h:1792
disti::operator==
bool operator==(const AttributeName &attr1, const AttributeName &attr2)
Definition: disti_metadata.h:154
disti::operator+=
GlsPropString & operator+=(GlsPropString &a, const GlsPropString &b)
Definition: util.h:1846