src/Imath/ImathMatrix.h

   1 //
   2 // SPDX-License-Identifier: BSD-3-Clause
   3 // Copyright Contributors to the OpenEXR Project.
   4 //
   5
   6 //
   7 // 2x2, 3x3, and 4x4 transformation matrix templates
   8 //
   9
  10 #ifndef INCLUDED_IMATHMATRIX_H
  11 #define INCLUDED_IMATHMATRIX_H
  12
  13 #include "ImathExport.h"
  14 #include "ImathNamespace.h"
  15
  16 #include "ImathFun.h"
  17 #include "ImathPlatform.h"
  18 #include "ImathShear.h"
  19 #include "ImathVec.h"
  20
  21 #include <cstring>
  22 #include <iomanip>
  23 #include <iostream>
  24 #include <limits>
  25 #include <string.h>
  26
  27 #if (defined _WIN32 || defined _WIN64) && defined _MSC_VER
  28 // suppress exception specification warnings
  29 #    pragma warning(disable : 4290)
  30 #endif
  31
  32 IMATH_INTERNAL_NAMESPACE_HEADER_ENTER
  33
  34 /// Enum used to indicate uninitialized construction of Matrix22,
  35 /// Matrix33, Matrix44
  36 enum IMATH_EXPORT_ENUM Uninitialized
  37 {
  38     UNINITIALIZED
  39 };
  40
  41 ///
  42 /// 2x2 transformation matrix
  43 ///
  44
  45 template <class T> class IMATH_EXPORT_TEMPLATE_TYPE Matrix22
  46 {
  47   public:
  48
  49     /// @{
  50     /// @name Direct access to elements
  51
  52     /// Matrix elements
  53     T x[2][2];
  54
  55     /// @}
  56
  57     /// Row access
  58     IMATH_HOSTDEVICE T* operator[] (int i) IMATH_NOEXCEPT;
  59
  60     /// Row access
  61     IMATH_HOSTDEVICE const T* operator[] (int i) const IMATH_NOEXCEPT;
  62
  63     /// @{
  64     /// @name Constructors and Assignment
  65
  66     /// Uninitialized
  67     IMATH_HOSTDEVICE Matrix22 (Uninitialized) IMATH_NOEXCEPT {}
  68
  69     /// Default constructor: initialize to identity
  70     ///
  71     ///     1 0
  72     ///     0 1
  73     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix22() IMATH_NOEXCEPT;
  74
  75     /// Initialize to scalar constant:
  76     ///
  77     ///     a a
  78     ///     a a
  79     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix22 (T a) IMATH_NOEXCEPT;
  80
  81     /// Construct from 2x2 array:
  82     ///
  83     ///     a[0][0] a[0][1]
  84     ///     a[1][0] a[1][1]
  85     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix22 (const T a[2][2]) IMATH_NOEXCEPT;
  86     /// Construct from given scalar values:
  87     ///
  88     ///     a b
  89     ///     c d
  90     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix22 (T a, T b, T c, T d) IMATH_NOEXCEPT;
  91
  92     /// Copy constructor
  93     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix22 (const Matrix22& v) IMATH_NOEXCEPT;
  94
  95     /// Construct from Matrix22 of another base type
  96     template <class S> IMATH_HOSTDEVICE IMATH_CONSTEXPR14 explicit Matrix22 (const Matrix22<S>& v) IMATH_NOEXCEPT;
  97
  98     /// Assignment
  99     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix22& operator= (const Matrix22& v) IMATH_NOEXCEPT;
 100
 101     /// Assignment from scalar
 102     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix22& operator= (T a) IMATH_NOEXCEPT;
 103
 104     /// Destructor
 105     ~Matrix22() IMATH_NOEXCEPT = default;
 106
 107     /// @}
 108
 109 #if IMATH_FOREIGN_VECTOR_INTEROP
 110     /// @{
 111     /// @name Interoperability with other matrix types
 112     ///
 113     /// Construction and assignment are allowed from other classes that
 114     /// appear to be equivalent matrix types, provided that they support
 115     /// double-subscript (i.e., `m[j][i]`) giving the same type as the
 116     /// elements of this matrix, and their total size appears to be the
 117     /// right number of matrix elements.
 118     ///
 119     /// This functionality is disabled for gcc 4.x, which seems to have a
 120     /// compiler bug that results in spurious errors. It can also be
 121     /// disabled by defining IMATH_FOREIGN_VECTOR_INTEROP to be 0 prior to
 122     /// including any Imath header files.
 123     ///
 124     template<typename M, IMATH_ENABLE_IF(has_double_subscript<M,T,2,2>::value)>
 125     IMATH_HOSTDEVICE explicit Matrix22 (const M& m)
 126         : Matrix22(T(m[0][0]), T(m[0][1]), T(m[1][0]), T(m[1][1]))
 127     { }
 128
 129     template<typename M, IMATH_ENABLE_IF(has_double_subscript<M,T,2,2>::value)>
 130     IMATH_HOSTDEVICE const Matrix22& operator= (const M& m)
 131     {
 132         *this = Matrix22(T(m[0][0]), T(m[0][1]), T(m[1][0]), T(m[1][1]));
 133         return *this;
 134     }
 135     /// @}
 136 #endif
 137
 138     /// @{
 139     /// @name Compatibility with Sb
 140
 141     /// Return a raw pointer to the array of values
 142     IMATH_HOSTDEVICE T* getValue() IMATH_NOEXCEPT;
 143
 144     /// Return a raw pointer to the array of values
 145     IMATH_HOSTDEVICE const T* getValue() const IMATH_NOEXCEPT;
 146
 147     /// Return the value in `v`
 148     template <class S> IMATH_HOSTDEVICE void getValue (Matrix22<S>& v) const IMATH_NOEXCEPT;
 149
 150     /// Set the value
 151     template <class S> IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix22& setValue (const Matrix22<S>& v) IMATH_NOEXCEPT;
 152
 153     /// Set the value
 154     template <class S>
 155     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix22& setTheMatrix (const Matrix22<S>& v) IMATH_NOEXCEPT;
 156
 157     /// @}
 158
 159     /// @{
 160     /// @name Arithmetic and Comparison
 161
 162     /// Equality
 163     IMATH_HOSTDEVICE constexpr bool operator== (const Matrix22& v) const IMATH_NOEXCEPT;
 164
 165     /// Inequality
 166     IMATH_HOSTDEVICE constexpr bool operator!= (const Matrix22& v) const IMATH_NOEXCEPT;
 167
 168     /// Compare two matrices and test if they are "approximately equal":
 169     /// @return True if the coefficients of this and `m` are the same
 170     /// with an absolute error of no more than e, i.e., for all i, j:
 171     ///
 172     ///     abs (this[i][j] - m[i][j]) <= e
 173     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 bool equalWithAbsError (const Matrix22<T>& v, T e) const IMATH_NOEXCEPT;
 174
 175     /// Compare two matrices and test if they are "approximately equal":
 176     /// @return True if the coefficients of this and m are the same with
 177     /// a relative error of no more than e, i.e., for all i, j:
 178     ///
 179     ///     abs (this[i] - v[i][j]) <= e * abs (this[i][j])
 180     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 bool equalWithRelError (const Matrix22<T>& v, T e) const IMATH_NOEXCEPT;
 181
 182     /// Component-wise addition
 183     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix22& operator+= (const Matrix22& v) IMATH_NOEXCEPT;
 184
 185     /// Component-wise addition
 186     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix22& operator+= (T a) IMATH_NOEXCEPT;
 187
 188     /// Component-wise addition
 189     IMATH_HOSTDEVICE constexpr Matrix22 operator+ (const Matrix22& v) const IMATH_NOEXCEPT;
 190
 191     /// Component-wise subtraction
 192     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix22& operator-= (const Matrix22& v) IMATH_NOEXCEPT;
 193
 194     /// Component-wise subtraction
 195     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix22& operator-= (T a) IMATH_NOEXCEPT;
 196
 197     /// Component-wise subtraction
 198     IMATH_HOSTDEVICE constexpr Matrix22 operator- (const Matrix22& v) const IMATH_NOEXCEPT;
 199
 200     /// Component-wise multiplication by -1
 201     IMATH_HOSTDEVICE constexpr Matrix22 operator-() const IMATH_NOEXCEPT;
 202
 203     /// Component-wise multiplication by -1
 204     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix22& negate() IMATH_NOEXCEPT;
 205
 206     /// Component-wise multiplication
 207     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix22& operator*= (T a) IMATH_NOEXCEPT;
 208
 209     /// Component-wise multiplication
 210     IMATH_HOSTDEVICE constexpr Matrix22 operator* (T a) const IMATH_NOEXCEPT;
 211
 212     /// Component-wise division
 213     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix22& operator/= (T a) IMATH_NOEXCEPT;
 214
 215     /// Component-wise division
 216     IMATH_HOSTDEVICE constexpr Matrix22 operator/ (T a) const IMATH_NOEXCEPT;
 217
 218     /// Matrix-matrix multiplication
 219     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix22& operator*= (const Matrix22& v) IMATH_NOEXCEPT;
 220
 221     /// Matrix-matrix multiplication
 222     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix22 operator* (const Matrix22& v) const IMATH_NOEXCEPT;
 223
 224     /// Vector * matrix multiplication
 225     /// @param[in] src Input vector
 226     /// @param[out] dst transformed vector
 227     template <class S> IMATH_HOSTDEVICE void multDirMatrix (const Vec2<S>& src, Vec2<S>& dst) const IMATH_NOEXCEPT;
 228
 229     /// @}
 230
 231     /// @{
 232     /// @name Maniplation
 233
 234     /// Set to the identity
 235     IMATH_HOSTDEVICE void makeIdentity() IMATH_NOEXCEPT;
 236
 237     /// Transpose
 238     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix22& transpose() IMATH_NOEXCEPT;
 239
 240     /// Return the transpose
 241     IMATH_HOSTDEVICE constexpr Matrix22 transposed() const IMATH_NOEXCEPT;
 242
 243     /// Invert in place
 244     /// @param singExc If true, throw an exception if the matrix cannot be inverted.
 245     /// @return const reference to this
 246     IMATH_CONSTEXPR14 const Matrix22& invert (bool singExc);
 247
 248     /// Invert in place
 249     /// @return const reference to this
 250     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix22& invert() IMATH_NOEXCEPT;
 251
 252     /// Return the inverse, leaving this unmodified.
 253     /// @param singExc If true, throw an exception if the matrix cannot be inverted.
 254     IMATH_CONSTEXPR14 Matrix22<T> inverse (bool singExc) const;
 255
 256     /// Return the inverse, leaving this unmodified.
 257     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix22<T> inverse() const IMATH_NOEXCEPT;
 258
 259     /// Determinant
 260     IMATH_HOSTDEVICE constexpr T determinant() const IMATH_NOEXCEPT;
 261
 262     /// Trace
 263     IMATH_HOSTDEVICE constexpr T trace() const IMATH_NOEXCEPT;
 264
 265     /// Set matrix to rotation by r (in radians)
 266     /// @return const referenced to this
 267     template <class S> IMATH_HOSTDEVICE const Matrix22& setRotation (S r) IMATH_NOEXCEPT;
 268
 269     /// Rotate the given matrix by r (in radians)
 270     /// @return const referenced to this
 271     template <class S> IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix22& rotate (S r) IMATH_NOEXCEPT;
 272
 273     /// Set matrix to scale by given uniform factor
 274     /// @return const referenced to this
 275     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix22& setScale (T s) IMATH_NOEXCEPT;
 276
 277     /// Set matrix to scale by given vector
 278     /// @return const referenced to this
 279     template <class S>
 280     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix22& setScale (const Vec2<S>& s) IMATH_NOEXCEPT;
 281
 282     // Scale the matrix by s
 283     /// @return const referenced to this
 284     template <class S> IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix22& scale (const Vec2<S>& s) IMATH_NOEXCEPT;
 285
 286     /// @}
 287
 288     /// @{
 289     /// @name Numeric Limits
 290
 291     /// Largest possible negative value
 292     IMATH_HOSTDEVICE constexpr static T baseTypeLowest() IMATH_NOEXCEPT { return std::numeric_limits<T>::lowest(); }
 293
 294     /// Largest possible positive value
 295     IMATH_HOSTDEVICE constexpr static T baseTypeMax() IMATH_NOEXCEPT { return std::numeric_limits<T>::max(); }
 296
 297     /// Smallest possible positive value
 298     IMATH_HOSTDEVICE constexpr static T baseTypeSmallest() IMATH_NOEXCEPT { return std::numeric_limits<T>::min(); }
 299
 300     /// Smallest possible e for which 1+e != 1
 301     IMATH_HOSTDEVICE constexpr static T baseTypeEpsilon() IMATH_NOEXCEPT { return std::numeric_limits<T>::epsilon(); }
 302
 303     /// @}
 304
 305     /// Return the number of the row and column dimensions, i.e. 2.
 306     IMATH_HOSTDEVICE constexpr static unsigned int dimensions() IMATH_NOEXCEPT { return 2; }
 307
 308     /// The base type: In templates that accept a parameter `V`, you
 309     /// can refer to `T` as `V::BaseType`
 310     typedef T BaseType;
 311
 312     /// The base vector type
 313     typedef Vec2<T> BaseVecType;
 314 };
 315
 316 ///
 317 /// 3x3 transformation matrix
 318 ///
 319
 320 template <class T> class IMATH_EXPORT_TEMPLATE_TYPE Matrix33
 321 {
 322   public:
 323
 324     /// @{
 325     /// @name Direct access to elements
 326
 327     /// Matrix elements
 328     T x[3][3];
 329
 330     /// @}
 331
 332     /// Row access
 333     IMATH_HOSTDEVICE T* operator[] (int i) IMATH_NOEXCEPT;
 334
 335     /// Row access
 336     IMATH_HOSTDEVICE const T* operator[] (int i) const IMATH_NOEXCEPT;
 337
 338     /// @{
 339     /// @name Constructors and Assignment
 340
 341     /// Uninitialized
 342     IMATH_HOSTDEVICE Matrix33 (Uninitialized) IMATH_NOEXCEPT {}
 343
 344     /// Default constructor: initialize to identity
 345     ///     1 0 0
 346     ///     0 1 0
 347     ///     0 0 1
 348     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix33() IMATH_NOEXCEPT;
 349
 350     /// Initialize to scalar constant
 351     ///     a a a
 352     ///     a a a
 353     ///     a a a
 354     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix33 (T a) IMATH_NOEXCEPT;
 355
 356     /// Construct from 3x3 array
 357     ///     a[0][0] a[0][1] a[0][2]
 358     ///     a[1][0] a[1][1] a[1][2]
 359     ///     a[2][0] a[2][1] a[2][2]
 360     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix33 (const T a[3][3]) IMATH_NOEXCEPT;
 361     /// Construct from given scalar values
 362     ///     a b c
 363     ///     d e f
 364     ///     g h i
 365     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix33 (T a, T b, T c, T d, T e, T f, T g, T h, T i) IMATH_NOEXCEPT;
 366
 367     /// Copy constructor
 368     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix33 (const Matrix33& v) IMATH_NOEXCEPT;
 369
 370     /// Construct from Matrix33 of another base type
 371     template <class S> IMATH_HOSTDEVICE IMATH_CONSTEXPR14 explicit Matrix33 (const Matrix33<S>& v) IMATH_NOEXCEPT;
 372
 373     /// Assignment operator
 374     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix33& operator= (const Matrix33& v) IMATH_NOEXCEPT;
 375
 376     /// Assignment from scalar
 377     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix33& operator= (T a) IMATH_NOEXCEPT;
 378
 379     /// Destructor
 380     ~Matrix33() IMATH_NOEXCEPT = default;
 381
 382     /// @}
 383
 384 #if IMATH_FOREIGN_VECTOR_INTEROP
 385     /// @{
 386     /// @name Interoperability with other matrix types
 387     ///
 388     /// Construction and assignment are allowed from other classes that
 389     /// appear to be equivalent matrix types, provided that they support
 390     /// double-subscript (i.e., `m[j][i]`) giving the same type as the
 391     /// elements of this matrix, and their total size appears to be the
 392     /// right number of matrix elements.
 393     ///
 394     /// This functionality is disabled for gcc 4.x, which seems to have a
 395     /// compiler bug that results in spurious errors. It can also be
 396     /// disabled by defining IMATH_FOREIGN_VECTOR_INTEROP to be 0 prior to
 397     /// including any Imath header files.
 398     ///
 399     template<typename M, IMATH_ENABLE_IF(has_double_subscript<M,T,3,3>::value)>
 400     IMATH_HOSTDEVICE explicit Matrix33 (const M& m)
 401         : Matrix33(T(m[0][0]), T(m[0][1]), T(m[0][2]),
 402                    T(m[1][0]), T(m[1][1]), T(m[1][2]),
 403                    T(m[2][0]), T(m[2][1]), T(m[2][2]))
 404     { }
 405
 406     /// Interoperability assignment from another type that behaves as if it
 407     /// were an equivalent matrix.
 408     template<typename M, IMATH_ENABLE_IF(has_double_subscript<M,T,3,3>::value)>
 409     IMATH_HOSTDEVICE const Matrix33& operator= (const M& m)
 410     {
 411         *this = Matrix33(T(m[0][0]), T(m[0][1]), T(m[0][2]),
 412                          T(m[1][0]), T(m[1][1]), T(m[1][2]),
 413                          T(m[2][0]), T(m[2][1]), T(m[2][2]));
 414         return *this;
 415     }
 416     /// @}
 417 #endif
 418
 419     /// @{
 420     /// @name Compatibility with Sb
 421
 422     /// Return a raw pointer to the array of values
 423     IMATH_HOSTDEVICE T* getValue() IMATH_NOEXCEPT;
 424
 425     /// Return a raw pointer to the array of values
 426     IMATH_HOSTDEVICE const T* getValue() const IMATH_NOEXCEPT;
 427
 428     /// Return the value in `v`
 429     template <class S> IMATH_HOSTDEVICE void getValue (Matrix33<S>& v) const IMATH_NOEXCEPT;
 430
 431     /// Set the value
 432     template <class S> IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix33& setValue (const Matrix33<S>& v) IMATH_NOEXCEPT;
 433
 434     /// Set the value
 435     template <class S>
 436     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix33& setTheMatrix (const Matrix33<S>& v) IMATH_NOEXCEPT;
 437
 438     /// @}
 439
 440     /// @{
 441     /// @name Arithmetic and Comparison
 442
 443     /// Equality
 444     IMATH_HOSTDEVICE constexpr bool operator== (const Matrix33& v) const IMATH_NOEXCEPT;
 445
 446     /// Inequality
 447     IMATH_HOSTDEVICE constexpr bool operator!= (const Matrix33& v) const IMATH_NOEXCEPT;
 448
 449     /// Compare two matrices and test if they are "approximately equal":
 450     /// @return True if the coefficients of this and `m` are the same
 451     /// with an absolute error of no more than e, i.e., for all i, j:
 452     ///
 453     ///     abs (this[i][j] - m[i][j]) <= e
 454     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 bool equalWithAbsError (const Matrix33<T>& v, T e) const IMATH_NOEXCEPT;
 455
 456     /// Compare two matrices and test if they are "approximately equal":
 457     /// @return True if the coefficients of this and m are the same with
 458     /// a relative error of no more than e, i.e., for all i, j:
 459     ///
 460     ///     abs (this[i] - v[i][j]) <= e * abs (this[i][j])
 461     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 bool equalWithRelError (const Matrix33<T>& v, T e) const IMATH_NOEXCEPT;
 462
 463     /// Component-wise addition
 464     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix33& operator+= (const Matrix33& v) IMATH_NOEXCEPT;
 465
 466     /// Component-wise addition
 467     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix33& operator+= (T a) IMATH_NOEXCEPT;
 468
 469     /// Component-wise addition
 470     IMATH_HOSTDEVICE constexpr Matrix33 operator+ (const Matrix33& v) const IMATH_NOEXCEPT;
 471
 472     /// Component-wise subtraction
 473     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix33& operator-= (const Matrix33& v) IMATH_NOEXCEPT;
 474
 475     /// Component-wise subtraction
 476     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix33& operator-= (T a) IMATH_NOEXCEPT;
 477
 478     /// Component-wise subtraction
 479     IMATH_HOSTDEVICE constexpr Matrix33 operator- (const Matrix33& v) const IMATH_NOEXCEPT;
 480
 481     /// Component-wise multiplication by -1
 482     IMATH_HOSTDEVICE constexpr Matrix33 operator-() const IMATH_NOEXCEPT;
 483
 484     /// Component-wise multiplication by -1
 485     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix33& negate() IMATH_NOEXCEPT;
 486
 487     /// Component-wise multiplication
 488     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix33& operator*= (T a) IMATH_NOEXCEPT;
 489
 490     /// Component-wise multiplication
 491     IMATH_HOSTDEVICE constexpr Matrix33 operator* (T a) const IMATH_NOEXCEPT;
 492
 493     /// Component-wise division
 494     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix33& operator/= (T a) IMATH_NOEXCEPT;
 495
 496     /// Component-wise division
 497     IMATH_HOSTDEVICE constexpr Matrix33 operator/ (T a) const IMATH_NOEXCEPT;
 498
 499     /// Matrix-matrix multiplication
 500     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix33& operator*= (const Matrix33& v) IMATH_NOEXCEPT;
 501
 502     /// Matrix-matrix multiplication
 503     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix33 operator* (const Matrix33& v) const IMATH_NOEXCEPT;
 504
 505     /// Vector-matrix multiplication: a homogeneous transformation
 506     /// by computing Vec3 (src.x, src.y, 1) * m and dividing by the
 507     /// result's third element.
 508     /// @param[in] src The input vector
 509     /// @param[out] dst The output vector
 510     template <class S> IMATH_HOSTDEVICE void multVecMatrix (const Vec2<S>& src, Vec2<S>& dst) const IMATH_NOEXCEPT;
 511
 512     /// Vector-matrix multiplication: multiply `src` by the upper left 2x2
 513     /// submatrix, ignoring the rest of matrix.
 514     /// @param[in] src The input vector
 515     /// @param[out] dst The output vector
 516     template <class S> IMATH_HOSTDEVICE void multDirMatrix (const Vec2<S>& src, Vec2<S>& dst) const IMATH_NOEXCEPT;
 517
 518     /// @}
 519
 520     /// @{
 521     /// @name Maniplation
 522
 523     /// Set to the identity matrix
 524     IMATH_HOSTDEVICE void makeIdentity() IMATH_NOEXCEPT;
 525
 526     /// Transpose
 527     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix33& transpose() IMATH_NOEXCEPT;
 528
 529     /// Return the transpose
 530     IMATH_HOSTDEVICE constexpr Matrix33 transposed() const IMATH_NOEXCEPT;
 531
 532     /// Invert in place using the determinant.
 533     /// @param singExc If true, throw an exception if the matrix cannot be inverted.
 534     /// @return const reference to this
 535     IMATH_CONSTEXPR14 const Matrix33& invert (bool singExc);
 536
 537     /// Invert in place using the determinant.
 538     /// @return const reference to this
 539     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix33& invert() IMATH_NOEXCEPT;
 540
 541     /// Return the inverse using the determinant, leaving this unmodified.
 542     /// @param singExc If true, throw an exception if the matrix cannot be inverted.
 543     IMATH_CONSTEXPR14 Matrix33<T> inverse (bool singExc) const;
 544
 545     /// Return the inverse using the determinant, leaving this unmodified.
 546     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix33<T> inverse() const IMATH_NOEXCEPT;
 547
 548     /// Invert in place using the Gauss-Jordan method. Significantly slower
 549     /// but more accurate than invert().
 550     /// @param singExc If true, throw an exception if the matrix cannot be inverted.
 551     /// @return const reference to this
 552     const Matrix33& gjInvert (bool singExc);
 553
 554     /// Invert in place using the Gauss-Jordan method. Significantly slower
 555     /// but more accurate than invert().
 556     /// @return const reference to this
 557     IMATH_HOSTDEVICE const Matrix33& gjInvert() IMATH_NOEXCEPT;
 558
 559     /// Return the inverse using the Gauss-Jordan method, leaving this
 560     /// unmodified. Significantly slower but more accurate than inverse().
 561     Matrix33<T> gjInverse (bool singExc) const;
 562
 563     /// Return the inverse using the Gauss-Jordan method. Significantly slower,
 564     /// leaving this unmodified. Slower but more accurate than inverse().
 565     IMATH_HOSTDEVICE Matrix33<T> gjInverse() const IMATH_NOEXCEPT;
 566
 567     /// Calculate the matrix minor of the (r,c) element
 568     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 T minorOf (const int r, const int c) const IMATH_NOEXCEPT;
 569
 570     /// Build a minor using the specified rows and columns
 571     IMATH_HOSTDEVICE
 572     constexpr T fastMinor (const int r0, const int r1, const int c0, const int c1) const IMATH_NOEXCEPT;
 573
 574     /// Determinant
 575     IMATH_HOSTDEVICE constexpr T determinant() const IMATH_NOEXCEPT;
 576
 577     /// Trace
 578     IMATH_HOSTDEVICE constexpr T trace() const IMATH_NOEXCEPT;
 579
 580     /// Set matrix to rotation by r (in radians)
 581     /// @return const referenced to this
 582     template <class S> IMATH_HOSTDEVICE const Matrix33& setRotation (S r) IMATH_NOEXCEPT;
 583
 584     // Rotate the given matrix by r (in radians)
 585     /// @return const referenced to this
 586     template <class S> IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix33& rotate (S r) IMATH_NOEXCEPT;
 587
 588     /// Set matrix to scale by given uniform factor
 589     /// @return const referenced to this
 590     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix33& setScale (T s) IMATH_NOEXCEPT;
 591
 592     /// Set matrix to scale by given vector
 593     /// @return const referenced to this
 594     template <class S>
 595     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix33& setScale (const Vec2<S>& s) IMATH_NOEXCEPT;
 596
 597     /// Scale the matrix by s
 598     /// @return const referenced to this
 599     template <class S> IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix33& scale (const Vec2<S>& s) IMATH_NOEXCEPT;
 600
 601     /// Set matrix to translation by given vector
 602     /// @return const referenced to this
 603     template <class S>
 604     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix33& setTranslation (const Vec2<S>& t) IMATH_NOEXCEPT;
 605
 606     /// Return the translation component
 607     IMATH_HOSTDEVICE constexpr Vec2<T> translation() const IMATH_NOEXCEPT;
 608
 609     /// Translate the matrix by t
 610     /// @return const referenced to this
 611     template <class S>
 612     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix33& translate (const Vec2<S>& t) IMATH_NOEXCEPT;
 613
 614     /// Set matrix to shear x for each y coord. by given factor xy
 615     /// @return const referenced to this
 616     template <class S> IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix33& setShear (const S& h) IMATH_NOEXCEPT;
 617
 618     /// Set matrix to shear x for each y coord. by given factor h.x
 619     /// and to shear y for each x coord. by given factor h.y
 620     /// @return const referenced to this
 621     template <class S>
 622     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix33& setShear (const Vec2<S>& h) IMATH_NOEXCEPT;
 623
 624     /// Shear the matrix in x for each y coord. by given factor xy
 625     /// @return const referenced to this
 626     template <class S> IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix33& shear (const S& xy) IMATH_NOEXCEPT;
 627
 628     /// Shear the matrix in x for each y coord. by given factor xy
 629     /// and shear y for each x coord. by given factor yx
 630     /// @return const referenced to this
 631     template <class S> IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix33& shear (const Vec2<S>& h) IMATH_NOEXCEPT;
 632
 633     /// @}
 634
 635     /// @{
 636     /// @name Numeric Limits
 637
 638     /// Largest possible negative value
 639     IMATH_HOSTDEVICE constexpr static T baseTypeLowest() IMATH_NOEXCEPT { return std::numeric_limits<T>::lowest(); }
 640
 641     /// Largest possible positive value
 642     IMATH_HOSTDEVICE constexpr static T baseTypeMax() IMATH_NOEXCEPT { return std::numeric_limits<T>::max(); }
 643
 644     /// Smallest possible positive value
 645     IMATH_HOSTDEVICE constexpr static T baseTypeSmallest() IMATH_NOEXCEPT { return std::numeric_limits<T>::min(); }
 646
 647     /// Smallest possible e for which 1+e != 1
 648     IMATH_HOSTDEVICE constexpr static T baseTypeEpsilon() IMATH_NOEXCEPT { return std::numeric_limits<T>::epsilon(); }
 649
 650     /// @}
 651
 652     /// Return the number of the row and column dimensions, i.e. 3.
 653     IMATH_HOSTDEVICE constexpr static unsigned int dimensions() IMATH_NOEXCEPT { return 3; }
 654
 655     /// The base type: In templates that accept a parameter `V` (could be a Color4), you can refer to `T` as `V::BaseType`
 656     typedef T BaseType;
 657
 658     /// The base vector type
 659     typedef Vec3<T> BaseVecType;
 660 };
 661
 662 ///
 663 /// 4x4 transformation matrix
 664 ///
 665
 666 template <class T> class IMATH_EXPORT_TEMPLATE_TYPE Matrix44
 667 {
 668   public:
 669
 670     /// @{
 671     /// @name Direct access to elements
 672
 673     /// Matrix elements
 674     T x[4][4];
 675
 676     /// @}
 677
 678     /// Row access
 679     IMATH_HOSTDEVICE T* operator[] (int i) IMATH_NOEXCEPT;
 680
 681     /// Row access
 682     IMATH_HOSTDEVICE const T* operator[] (int i) const IMATH_NOEXCEPT;
 683
 684     /// @{
 685     /// @name Constructors and Assignment
 686
 687     /// Uninitialized
 688     IMATH_HOSTDEVICE constexpr Matrix44 (Uninitialized) IMATH_NOEXCEPT {}
 689
 690     /// Default constructor: initialize to identity
 691     ///     1 0 0 0
 692     ///     0 1 0 0
 693     ///     0 0 1 0
 694     ///     0 0 0 1
 695     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix44() IMATH_NOEXCEPT;
 696
 697     /// Initialize to scalar constant
 698     ///     a a a a
 699     ///     a a a a
 700     ///     a a a a
 701     ///     a a a a
 702     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix44 (T a) IMATH_NOEXCEPT;
 703
 704     /// Construct from 4x4 array
 705     ///     a[0][0] a[0][1] a[0][2] a[0][3]
 706     ///     a[1][0] a[1][1] a[1][2] a[1][3]
 707     ///     a[2][0] a[2][1] a[2][2] a[2][3]
 708     ///     a[3][0] a[3][1] a[3][2] a[3][3]
 709     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix44 (const T a[4][4]) IMATH_NOEXCEPT;
 710     /// Construct from given scalar values
 711     ///     a b c d
 712     ///     e f g h
 713     ///     i j k l
 714     ///     m n o p
 715     IMATH_HOSTDEVICE IMATH_CONSTEXPR14
 716     Matrix44 (T a, T b, T c, T d, T e, T f, T g, T h, T i, T j, T k, T l, T m, T n, T o, T p) IMATH_NOEXCEPT;
 717
 718
 719     /// Construct from a 3x3 rotation matrix and a translation vector
 720     ///     r r r 0
 721     ///     r r r 0
 722     ///     r r r 0
 723     ///     t t t 1
 724     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix44 (Matrix33<T> r, Vec3<T> t) IMATH_NOEXCEPT;
 725
 726     /// Copy constructor
 727     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix44 (const Matrix44& v) IMATH_NOEXCEPT;
 728
 729     /// Construct from Matrix44 of another base type
 730     template <class S> IMATH_HOSTDEVICE IMATH_CONSTEXPR14 explicit Matrix44 (const Matrix44<S>& v) IMATH_NOEXCEPT;
 731
 732     /// Assignment operator
 733     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix44& operator= (const Matrix44& v) IMATH_NOEXCEPT;
 734
 735     /// Assignment from scalar
 736     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix44& operator= (T a) IMATH_NOEXCEPT;
 737
 738     /// Destructor
 739     ~Matrix44() IMATH_NOEXCEPT = default;
 740
 741     /// @}
 742
 743 #if IMATH_FOREIGN_VECTOR_INTEROP
 744     /// @{
 745     /// @name Interoperability with other matrix types
 746     ///
 747     /// Construction and assignment are allowed from other classes that
 748     /// appear to be equivalent matrix types, provided that they support
 749     /// double-subscript (i.e., `m[j][i]`) giving the same type as the
 750     /// elements of this matrix, and their total size appears to be the
 751     /// right number of matrix elements.
 752     ///
 753     /// This functionality is disabled for gcc 4.x, which seems to have a
 754     /// compiler bug that results in spurious errors. It can also be
 755     /// disabled by defining IMATH_FOREIGN_VECTOR_INTEROP to be 0 prior to
 756     /// including any Imath header files.
 757     ///
 758     template<typename M, IMATH_ENABLE_IF(has_double_subscript<M,T,4,4>::value)>
 759     IMATH_HOSTDEVICE explicit Matrix44 (const M& m)
 760         : Matrix44(T(m[0][0]), T(m[0][1]), T(m[0][2]), T(m[0][3]),
 761                    T(m[1][0]), T(m[1][1]), T(m[1][2]), T(m[1][3]),
 762                    T(m[2][0]), T(m[2][1]), T(m[2][2]), T(m[2][3]),
 763                    T(m[3][0]), T(m[3][1]), T(m[3][2]), T(m[3][3]))
 764     { }
 765
 766     /// Interoperability assignment from another type that behaves as if it
 767     /// were an equivalent matrix.
 768     template<typename M, IMATH_ENABLE_IF(has_double_subscript<M,T,4,4>::value)>
 769     IMATH_HOSTDEVICE const Matrix44& operator= (const M& m)
 770     {
 771         *this = Matrix44(T(m[0][0]), T(m[0][1]), T(m[0][2]), T(m[0][3]),
 772                          T(m[1][0]), T(m[1][1]), T(m[1][2]), T(m[1][3]),
 773                          T(m[2][0]), T(m[2][1]), T(m[2][2]), T(m[2][3]),
 774                          T(m[3][0]), T(m[3][1]), T(m[3][2]), T(m[3][3]));
 775         return *this;
 776     }
 777     /// @}
 778 #endif
 779
 780     /// @{
 781     /// @name Compatibility with Sb
 782
 783     /// Return a raw pointer to the array of values
 784     IMATH_HOSTDEVICE T* getValue() IMATH_NOEXCEPT;
 785
 786     /// Return a raw pointer to the array of values
 787     IMATH_HOSTDEVICE const T* getValue() const IMATH_NOEXCEPT;
 788
 789     /// Return the value in `v`
 790     template <class S> IMATH_HOSTDEVICE void getValue (Matrix44<S>& v) const IMATH_NOEXCEPT;
 791
 792     /// Set the value
 793     template <class S> IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix44& setValue (const Matrix44<S>& v) IMATH_NOEXCEPT;
 794
 795     /// Set the value
 796     template <class S>
 797     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix44& setTheMatrix (const Matrix44<S>& v) IMATH_NOEXCEPT;
 798
 799     /// @}
 800
 801     /// @{
 802     /// @name Arithmetic and Comparison
 803
 804     /// Equality
 805     IMATH_HOSTDEVICE constexpr bool operator== (const Matrix44& v) const IMATH_NOEXCEPT;
 806
 807     /// Inequality
 808     IMATH_HOSTDEVICE constexpr bool operator!= (const Matrix44& v) const IMATH_NOEXCEPT;
 809
 810     /// Compare two matrices and test if they are "approximately equal":
 811     /// @return True if the coefficients of this and `m` are the same
 812     /// with an absolute error of no more than e, i.e., for all i, j:
 813     ///
 814     ///     abs (this[i][j] - m[i][j]) <= e
 815     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 bool equalWithAbsError (const Matrix44<T>& v, T e) const IMATH_NOEXCEPT;
 816
 817     /// Compare two matrices and test if they are "approximately equal":
 818     /// @return True if the coefficients of this and m are the same with
 819     /// a relative error of no more than e, i.e., for all i, j:
 820     ///
 821     ///     abs (this[i] - v[i][j]) <= e * abs (this[i][j])
 822     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 bool equalWithRelError (const Matrix44<T>& v, T e) const IMATH_NOEXCEPT;
 823
 824     /// Component-wise addition
 825     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix44& operator+= (const Matrix44& v) IMATH_NOEXCEPT;
 826
 827     /// Component-wise addition
 828     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix44& operator+= (T a) IMATH_NOEXCEPT;
 829
 830     /// Component-wise addition
 831     IMATH_HOSTDEVICE constexpr Matrix44 operator+ (const Matrix44& v) const IMATH_NOEXCEPT;
 832
 833     /// Component-wise subtraction
 834     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix44& operator-= (const Matrix44& v) IMATH_NOEXCEPT;
 835
 836     /// Component-wise subtraction
 837     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix44& operator-= (T a) IMATH_NOEXCEPT;
 838
 839     /// Component-wise subtraction
 840     IMATH_HOSTDEVICE constexpr Matrix44 operator- (const Matrix44& v) const IMATH_NOEXCEPT;
 841
 842     /// Component-wise multiplication by -1
 843     IMATH_HOSTDEVICE constexpr Matrix44 operator-() const IMATH_NOEXCEPT;
 844
 845     /// Component-wise multiplication by -1
 846     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix44& negate() IMATH_NOEXCEPT;
 847
 848     /// Component-wise multiplication
 849     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix44& operator*= (T a) IMATH_NOEXCEPT;
 850
 851     /// Component-wise multiplication
 852     IMATH_HOSTDEVICE constexpr Matrix44 operator* (T a) const IMATH_NOEXCEPT;
 853
 854     /// Component-wise division
 855     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix44& operator/= (T a) IMATH_NOEXCEPT;
 856
 857     /// Component-wise division
 858     IMATH_HOSTDEVICE constexpr Matrix44 operator/ (T a) const IMATH_NOEXCEPT;
 859
 860     /// Matrix-matrix multiplication
 861     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix44& operator*= (const Matrix44& v) IMATH_NOEXCEPT;
 862
 863     /// Matrix-matrix multiplication
 864     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix44 operator* (const Matrix44& v) const IMATH_NOEXCEPT;
 865
 866     /// Matrix-matrix multiplication: compute c = a * b
 867     IMATH_HOSTDEVICE
 868     static void multiply (const Matrix44& a,     // assumes that
 869                           const Matrix44& b,     // &a != &c and
 870                           Matrix44& c) IMATH_NOEXCEPT; // &b != &c.
 871
 872     /// Matrix-matrix multiplication returning a result.
 873     IMATH_HOSTDEVICE
 874     static IMATH_CONSTEXPR14 Matrix44 multiply (const Matrix44& a, const Matrix44& b) IMATH_NOEXCEPT;
 875
 876     /// Vector-matrix multiplication: a homogeneous transformation
 877     /// by computing Vec3 (src.x, src.y, src.z, 1) * m and dividing by the
 878     /// result's third element.
 879     /// @param[in] src The input vector
 880     /// @param[out] dst The output vector
 881     template <class S> IMATH_HOSTDEVICE void multVecMatrix (const Vec3<S>& src, Vec3<S>& dst) const IMATH_NOEXCEPT;
 882
 883     /// Vector-matrix multiplication: multiply `src` by the upper left 2x2
 884     /// submatrix, ignoring the rest of matrix.
 885     /// @param[in] src The input vector
 886     /// @param[out] dst The output vector
 887     template <class S> IMATH_HOSTDEVICE void multDirMatrix (const Vec3<S>& src, Vec3<S>& dst) const IMATH_NOEXCEPT;
 888
 889     /// @}
 890
 891     /// @{
 892     /// @name Maniplation
 893
 894     /// Set to the identity matrix
 895     IMATH_HOSTDEVICE void makeIdentity() IMATH_NOEXCEPT;
 896
 897     /// Transpose
 898     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix44& transpose() IMATH_NOEXCEPT;
 899
 900     /// Return the transpose
 901     IMATH_HOSTDEVICE constexpr Matrix44 transposed() const IMATH_NOEXCEPT;
 902
 903     /// Invert in place using the determinant.
 904     /// @param singExc If true, throw an exception if the matrix cannot be inverted.
 905     /// @return const reference to this
 906     IMATH_CONSTEXPR14 const Matrix44& invert (bool singExc);
 907
 908     /// Invert in place using the determinant.
 909     /// @return const reference to this
 910     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix44& invert() IMATH_NOEXCEPT;
 911
 912     /// Return the inverse using the determinant, leaving this unmodified.
 913     /// @param singExc If true, throw an exception if the matrix cannot be inverted.
 914     IMATH_CONSTEXPR14 Matrix44<T> inverse (bool singExc) const;
 915
 916     /// Return the inverse using the determinant, leaving this unmodified.
 917     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Matrix44<T> inverse() const IMATH_NOEXCEPT;
 918
 919     /// Invert in place using the Gauss-Jordan method. Significantly slower
 920     /// but more accurate than invert().
 921     /// @param singExc If true, throw an exception if the matrix cannot be inverted.
 922     /// @return const reference to this
 923     IMATH_CONSTEXPR14 const Matrix44& gjInvert (bool singExc);
 924
 925     /// Invert in place using the Gauss-Jordan method. Significantly slower
 926     /// but more accurate than invert().
 927     /// @return const reference to this
 928     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix44& gjInvert() IMATH_NOEXCEPT;
 929
 930     /// Return the inverse using the Gauss-Jordan method, leaving this
 931     /// unmodified. Significantly slower but more accurate than inverse().
 932     Matrix44<T> gjInverse (bool singExc) const;
 933
 934     /// Return the inverse using the Gauss-Jordan method, leaving this
 935     /// unmodified Significantly slower but more accurate than inverse().
 936     IMATH_HOSTDEVICE Matrix44<T> gjInverse() const IMATH_NOEXCEPT;
 937
 938     /// Calculate the matrix minor of the (r,c) element
 939     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 T minorOf (const int r, const int c) const IMATH_NOEXCEPT;
 940
 941     /// Build a minor using the specified rows and columns
 942     IMATH_HOSTDEVICE
 943     constexpr T fastMinor (const int r0,
 944                            const int r1,
 945                            const int r2,
 946                            const int c0,
 947                            const int c1,
 948                            const int c2) const IMATH_NOEXCEPT;
 949
 950     /// Determinant
 951     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 T determinant() const IMATH_NOEXCEPT;
 952
 953     /// Trace
 954     IMATH_HOSTDEVICE constexpr T trace() const IMATH_NOEXCEPT;
 955
 956     /// Set matrix to rotation by XYZ euler angles (in radians)
 957     /// @return const referenced to this
 958     template <class S> IMATH_HOSTDEVICE const Matrix44& setEulerAngles (const Vec3<S>& r) IMATH_NOEXCEPT;
 959
 960     /// Set matrix to rotation around given axis by given angle (in radians)
 961     /// @return const referenced to this
 962     template <class S>
 963     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix44& setAxisAngle (const Vec3<S>& ax, S ang) IMATH_NOEXCEPT;
 964
 965     /// Rotate the matrix by XYZ euler angles in r (in radians)
 966     /// @return const referenced to this
 967     template <class S> IMATH_HOSTDEVICE const Matrix44& rotate (const Vec3<S>& r) IMATH_NOEXCEPT;
 968
 969     /// Set matrix to scale by given uniform factor
 970     /// @return const referenced to this
 971     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix44& setScale (T s) IMATH_NOEXCEPT;
 972
 973     /// Set matrix to scale by given vector
 974     /// @return const referenced to this
 975     template <class S>
 976     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix44& setScale (const Vec3<S>& s) IMATH_NOEXCEPT;
 977
 978     /// Scale the matrix by s
 979     /// @return const referenced to this
 980     template <class S> IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix44& scale (const Vec3<S>& s) IMATH_NOEXCEPT;
 981
 982     /// Set matrix to translation by given vector
 983     /// @return const referenced to this
 984     template <class S>
 985     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix44& setTranslation (const Vec3<S>& t) IMATH_NOEXCEPT;
 986
 987     /// Return translation component
 988     IMATH_HOSTDEVICE constexpr const Vec3<T> translation() const IMATH_NOEXCEPT;
 989
 990     /// Translate the matrix by t
 991     /// @return const referenced to this
 992     template <class S>
 993     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix44& translate (const Vec3<S>& t) IMATH_NOEXCEPT;
 994
 995     /// Set matrix to shear by given vector h.  The resulting matrix
 996     /// - will shear x for each y coord. by a factor of h[0] ;
 997     /// - will shear x for each z coord. by a factor of h[1] ;
 998     /// - will shear y for each z coord. by a factor of h[2] .
 999     /// @return const referenced to this
1000     template <class S>
1001     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix44& setShear (const Vec3<S>& h) IMATH_NOEXCEPT;
1002
1003     /// Set matrix to shear by given factors.  The resulting matrix
1004     /// - will shear x for each y coord. by a factor of h.xy ;
1005     /// - will shear x for each z coord. by a factor of h.xz ;
1006     /// - will shear y for each z coord. by a factor of h.yz ;
1007     /// - will shear y for each x coord. by a factor of h.yx ;
1008     /// - will shear z for each x coord. by a factor of h.zx ;
1009     /// - will shear z for each y coord. by a factor of h.zy .
1010     /// @return const referenced to this
1011     template <class S>
1012     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix44& setShear (const Shear6<S>& h) IMATH_NOEXCEPT;
1013
1014     /// Shear the matrix by given vector.  The composed matrix
1015     /// will be `shear` * `this`, where the shear matrix ...
1016     /// - will shear x for each y coord. by a factor of h[0] ;
1017     /// - will shear x for each z coord. by a factor of h[1] ;
1018     /// - will shear y for each z coord. by a factor of h[2] .
1019     /// @return const referenced to this
1020     template <class S> IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix44& shear (const Vec3<S>& h) IMATH_NOEXCEPT;
1021
1022     /// Shear the matrix by the given factors.  The composed matrix
1023     /// will be `shear` * `this`, where the shear matrix ...
1024     /// - will shear x for each y coord. by a factor of h.xy ;
1025     /// - will shear x for each z coord. by a factor of h.xz ;
1026     /// - will shear y for each z coord. by a factor of h.yz ;
1027     /// - will shear y for each x coord. by a factor of h.yx ;
1028     /// - will shear z for each x coord. by a factor of h.zx ;
1029     /// - will shear z for each y coord. by a factor of h.zy .
1030     /// @return const referenced to this
1031     template <class S>
1032     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Matrix44& shear (const Shear6<S>& h) IMATH_NOEXCEPT;
1033
1034     /// @}
1035
1036     /// @{
1037     /// @name Numeric Limits
1038
1039     /// Largest possible negative value
1040     IMATH_HOSTDEVICE constexpr static T baseTypeLowest() IMATH_NOEXCEPT { return std::numeric_limits<T>::lowest(); }
1041
1042     /// Largest possible positive value
1043     IMATH_HOSTDEVICE constexpr static T baseTypeMax() IMATH_NOEXCEPT { return std::numeric_limits<T>::max(); }
1044
1045     /// Smallest possible positive value
1046     IMATH_HOSTDEVICE constexpr static T baseTypeSmallest() IMATH_NOEXCEPT { return std::numeric_limits<T>::min(); }
1047
1048     /// Smallest possible e for which 1+e != 1
1049     IMATH_HOSTDEVICE constexpr static T baseTypeEpsilon() IMATH_NOEXCEPT { return std::numeric_limits<T>::epsilon(); }
1050
1051     /// @}
1052
1053     /// Return the number of the row and column dimensions, i.e. 4
1054     IMATH_HOSTDEVICE constexpr static unsigned int dimensions() IMATH_NOEXCEPT { return 4; }
1055
1056     /// The base type: In templates that accept a parameter `V` (could be a Color4), you can refer to `T` as `V::BaseType`
1057     typedef T BaseType;
1058
1059     /// The base vector type
1060     typedef Vec4<T> BaseVecType;
1061 };
1062
1063 /// Stream output, as:
1064 ///     (m00 m01
1065 ///      m10 m11)
1066 template <class T> std::ostream& operator<< (std::ostream& s, const Matrix22<T>& m);
1067
1068 /// Stream output, as:
1069 ///     (m00 m01 m02
1070 ///      m10 m11 m12
1071 ///      m20 m21 m22)
1072 template <class T> std::ostream& operator<< (std::ostream& s, const Matrix33<T>& m);
1073
1074 /// Stream output, as:
1075 ///
1076 ///     (m00 m01 m02 m03
1077 ///      m10 m11 m12 m13
1078 ///      m20 m21 m22 m23
1079 ///      m30 m31 m32 m33)
1080 template <class T> std::ostream& operator<< (std::ostream& s, const Matrix44<T>& m);
1081
1082 //---------------------------------------------
1083 // Vector-times-matrix multiplication operators
1084 //---------------------------------------------
1085
1086 /// Vector-matrix multiplication: v *= m
1087 template <class S, class T>
1088 IMATH_HOSTDEVICE inline const Vec2<S>& operator*= (Vec2<S>& v, const Matrix22<T>& m) IMATH_NOEXCEPT;
1089
1090 /// Vector-matrix multiplication: r = v * m
1091 template <class S, class T>
1092 IMATH_HOSTDEVICE inline Vec2<S> operator* (const Vec2<S>& v, const Matrix22<T>& m) IMATH_NOEXCEPT;
1093
1094 /// Vector-matrix multiplication: v *= m
1095 template <class S, class T>
1096 IMATH_HOSTDEVICE inline const Vec2<S>& operator*= (Vec2<S>& v, const Matrix33<T>& m) IMATH_NOEXCEPT;
1097
1098 /// Vector-matrix multiplication: r = v * m
1099 template <class S, class T>
1100 IMATH_HOSTDEVICE inline Vec2<S> operator* (const Vec2<S>& v, const Matrix33<T>& m) IMATH_NOEXCEPT;
1101
1102 /// Vector-matrix multiplication: v *= m
1103 template <class S, class T>
1104 IMATH_HOSTDEVICE inline const Vec3<S>& operator*= (Vec3<S>& v, const Matrix33<T>& m) IMATH_NOEXCEPT;
1105
1106 /// Vector-matrix multiplication: r = v * m
1107 template <class S, class T>
1108 IMATH_HOSTDEVICE inline Vec3<S> operator* (const Vec3<S>& v, const Matrix33<T>& m) IMATH_NOEXCEPT;
1109
1110 /// Vector-matrix multiplication: v *= m
1111 template <class S, class T>
1112 IMATH_HOSTDEVICE inline const Vec3<S>& operator*= (Vec3<S>& v, const Matrix44<T>& m) IMATH_NOEXCEPT;
1113
1114 /// Vector-matrix multiplication: r = v * m
1115 template <class S, class T>
1116 IMATH_HOSTDEVICE inline Vec3<S> operator* (const Vec3<S>& v, const Matrix44<T>& m) IMATH_NOEXCEPT;
1117
1118 /// Vector-matrix multiplication: v *= m
1119 template <class S, class T>
1120 IMATH_HOSTDEVICE inline const Vec4<S>& operator*= (Vec4<S>& v, const Matrix44<T>& m) IMATH_NOEXCEPT;
1121
1122 /// Vector-matrix multiplication: r = v * m
1123 template <class S, class T>
1124 IMATH_HOSTDEVICE inline Vec4<S> operator* (const Vec4<S>& v, const Matrix44<T>& m) IMATH_NOEXCEPT;
1125
1126 //-------------------------
1127 // Typedefs for convenience
1128 //-------------------------
1129
1130 /// 2x2 matrix of float
1131 typedef Matrix22<float> M22f;
1132
1133 /// 2x2 matrix of double
1134 typedef Matrix22<double> M22d;
1135
1136 /// 3x3 matrix of float
1137 typedef Matrix33<float> M33f;
1138
1139 /// 3x3 matrix of double
1140 typedef Matrix33<double> M33d;
1141
1142 /// 4x4 matrix of float
1143 typedef Matrix44<float> M44f;
1144
1145 /// 4x4 matrix of double
1146 typedef Matrix44<double> M44d;
1147
1148 //---------------------------
1149 // Implementation of Matrix22
1150 //---------------------------
1151
1152 template <class T>
1153 IMATH_HOSTDEVICE inline T*
1154 Matrix22<T>::operator[] (int i) IMATH_NOEXCEPT
1155 {
1156     return x[i];
1157 }
1158
1159 template <class T>
1160 IMATH_HOSTDEVICE inline const T*
1161 Matrix22<T>::operator[] (int i) const IMATH_NOEXCEPT
1162 {
1163     return x[i];
1164 }
1165
1166 template <class T> IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix22<T>::Matrix22() IMATH_NOEXCEPT
1167 {
1168     x[0][0] = 1;
1169     x[0][1] = 0;
1170     x[1][0] = 0;
1171     x[1][1] = 1;
1172 }
1173
1174 template <class T> IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix22<T>::Matrix22 (T a) IMATH_NOEXCEPT
1175 {
1176     x[0][0] = a;
1177     x[0][1] = a;
1178     x[1][0] = a;
1179     x[1][1] = a;
1180 }
1181
1182 template <class T>
1183 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix22<T>::Matrix22 (
1184     const T a[2][2]) IMATH_NOEXCEPT
1185 {
1186     // Function calls and aliasing issues can inhibit vectorization versus
1187     // straight assignment of data members, so instead of this:
1188     //     memcpy (x, a, sizeof (x));
1189     // we do this:
1190     x[0][0] = a[0][0];
1191     x[0][1] = a[0][1];
1192     x[1][0] = a[1][0];
1193     x[1][1] = a[1][1];
1194 }
1195
1196 template <class T>
1197 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix22<T>::Matrix22 (
1198     T a, T b, T c, T d) IMATH_NOEXCEPT
1199 {
1200     x[0][0] = a;
1201     x[0][1] = b;
1202     x[1][0] = c;
1203     x[1][1] = d;
1204 }
1205
1206 template <class T> IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix22<T>::Matrix22 (const Matrix22& v) IMATH_NOEXCEPT
1207 {
1208     // Function calls and aliasing issues can inhibit vectorization versus
1209     // straight assignment of data members, so we don't do this:
1210     //     memcpy (x, v.x, sizeof (x));
1211     // we do this:
1212     x[0][0] = v.x[0][0];
1213     x[0][1] = v.x[0][1];
1214     x[1][0] = v.x[1][0];
1215     x[1][1] = v.x[1][1];
1216 }
1217
1218 template <class T>
1219 template <class S>
1220 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix22<T>::Matrix22 (const Matrix22<S>& v) IMATH_NOEXCEPT
1221 {
1222     x[0][0] = T (v.x[0][0]);
1223     x[0][1] = T (v.x[0][1]);
1224     x[1][0] = T (v.x[1][0]);
1225     x[1][1] = T (v.x[1][1]);
1226 }
1227
1228 template <class T>
1229 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix22<T>&
1230 Matrix22<T>::operator= (const Matrix22& v) IMATH_NOEXCEPT
1231 {
1232     // Function calls and aliasing issues can inhibit vectorization versus
1233     // straight assignment of data members, so we don't do this:
1234     //     memcpy (x, v.x, sizeof (x));
1235     // we do this:
1236     x[0][0] = v.x[0][0];
1237     x[0][1] = v.x[0][1];
1238     x[1][0] = v.x[1][0];
1239     x[1][1] = v.x[1][1];
1240     return *this;
1241 }
1242
1243 template <class T>
1244 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix22<T>&
1245 Matrix22<T>::operator= (T a) IMATH_NOEXCEPT
1246 {
1247     x[0][0] = a;
1248     x[0][1] = a;
1249     x[1][0] = a;
1250     x[1][1] = a;
1251     return *this;
1252 }
1253
1254 template <class T>
1255 IMATH_HOSTDEVICE inline T*
1256 Matrix22<T>::getValue () IMATH_NOEXCEPT
1257 {
1258     return (T*) &x[0][0];
1259 }
1260
1261 template <class T>
1262 IMATH_HOSTDEVICE inline const T*
1263 Matrix22<T>::getValue() const IMATH_NOEXCEPT
1264 {
1265     return (const T*) &x[0][0];
1266 }
1267
1268 template <class T>
1269 template <class S>
1270 IMATH_HOSTDEVICE inline void
1271 Matrix22<T>::getValue (Matrix22<S>& v) const IMATH_NOEXCEPT
1272 {
1273     v.x[0][0] = x[0][0];
1274     v.x[0][1] = x[0][1];
1275     v.x[1][0] = x[1][0];
1276     v.x[1][1] = x[1][1];
1277 }
1278
1279 template <class T>
1280 template <class S>
1281 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix22<T>&
1282 Matrix22<T>::setValue (const Matrix22<S>& v) IMATH_NOEXCEPT
1283 {
1284     x[0][0] = v.x[0][0];
1285     x[0][1] = v.x[0][1];
1286     x[1][0] = v.x[1][0];
1287     x[1][1] = v.x[1][1];
1288     return *this;
1289 }
1290
1291 template <class T>
1292 template <class S>
1293 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix22<T>&
1294 Matrix22<T>::setTheMatrix (const Matrix22<S>& v) IMATH_NOEXCEPT
1295 {
1296     x[0][0] = v.x[0][0];
1297     x[0][1] = v.x[0][1];
1298     x[1][0] = v.x[1][0];
1299     x[1][1] = v.x[1][1];
1300     return *this;
1301 }
1302
1303 template <class T>
1304 IMATH_HOSTDEVICE inline void
1305 Matrix22<T>::makeIdentity() IMATH_NOEXCEPT
1306 {
1307     x[0][0] = 1;
1308     x[0][1] = 0;
1309     x[1][0] = 0;
1310     x[1][1] = 1;
1311 }
1312
1313 template <class T>
1314 IMATH_HOSTDEVICE constexpr inline bool
1315 Matrix22<T>::operator== (const Matrix22& v) const IMATH_NOEXCEPT
1316 {
1317     return x[0][0] == v.x[0][0] && x[0][1] == v.x[0][1] && x[1][0] == v.x[1][0] &&
1318            x[1][1] == v.x[1][1];
1319 }
1320
1321 template <class T>
1322 IMATH_HOSTDEVICE constexpr inline bool
1323 Matrix22<T>::operator!= (const Matrix22& v) const IMATH_NOEXCEPT
1324 {
1325     return x[0][0] != v.x[0][0] || x[0][1] != v.x[0][1] || x[1][0] != v.x[1][0] ||
1326            x[1][1] != v.x[1][1];
1327 }
1328
1329 template <class T>
1330 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline bool
1331 Matrix22<T>::equalWithAbsError (const Matrix22<T>& m, T e) const IMATH_NOEXCEPT
1332 {
1333     for (int i = 0; i < 2; i++)
1334         for (int j = 0; j < 2; j++)
1335             if (!IMATH_INTERNAL_NAMESPACE::equalWithAbsError ((*this).x[i][j], m.x[i][j], e))
1336                 return false;
1337
1338     return true;
1339 }
1340
1341 template <class T>
1342 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline bool
1343 Matrix22<T>::equalWithRelError (const Matrix22<T>& m, T e) const IMATH_NOEXCEPT
1344 {
1345     for (int i = 0; i < 2; i++)
1346         for (int j = 0; j < 2; j++)
1347             if (!IMATH_INTERNAL_NAMESPACE::equalWithRelError ((*this).x[i][j], m.x[i][j], e))
1348                 return false;
1349
1350     return true;
1351 }
1352
1353 template <class T>
1354 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix22<T>&
1355 Matrix22<T>::operator+= (const Matrix22<T>& v) IMATH_NOEXCEPT
1356 {
1357     x[0][0] += v.x[0][0];
1358     x[0][1] += v.x[0][1];
1359     x[1][0] += v.x[1][0];
1360     x[1][1] += v.x[1][1];
1361
1362     return *this;
1363 }
1364
1365 template <class T>
1366 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix22<T>&
1367 Matrix22<T>::operator+= (T a) IMATH_NOEXCEPT
1368 {
1369     x[0][0] += a;
1370     x[0][1] += a;
1371     x[1][0] += a;
1372     x[1][1] += a;
1373
1374     return *this;
1375 }
1376
1377 template <class T>
1378 IMATH_HOSTDEVICE constexpr inline Matrix22<T>
1379 Matrix22<T>::operator+ (const Matrix22<T>& v) const IMATH_NOEXCEPT
1380 {
1381     return Matrix22 (x[0][0] + v.x[0][0],
1382                      x[0][1] + v.x[0][1],
1383                      x[1][0] + v.x[1][0],
1384                      x[1][1] + v.x[1][1]);
1385 }
1386
1387 template <class T>
1388 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix22<T>&
1389 Matrix22<T>::operator-= (const Matrix22<T>& v) IMATH_NOEXCEPT
1390 {
1391     x[0][0] -= v.x[0][0];
1392     x[0][1] -= v.x[0][1];
1393     x[1][0] -= v.x[1][0];
1394     x[1][1] -= v.x[1][1];
1395
1396     return *this;
1397 }
1398
1399 template <class T>
1400 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix22<T>&
1401 Matrix22<T>::operator-= (T a) IMATH_NOEXCEPT
1402 {
1403     x[0][0] -= a;
1404     x[0][1] -= a;
1405     x[1][0] -= a;
1406     x[1][1] -= a;
1407
1408     return *this;
1409 }
1410
1411 template <class T>
1412 IMATH_HOSTDEVICE constexpr inline Matrix22<T>
1413 Matrix22<T>::operator- (const Matrix22<T>& v) const IMATH_NOEXCEPT
1414 {
1415     return Matrix22 (x[0][0] - v.x[0][0],
1416                      x[0][1] - v.x[0][1],
1417                      x[1][0] - v.x[1][0],
1418                      x[1][1] - v.x[1][1]);
1419 }
1420
1421 template <class T>
1422 IMATH_HOSTDEVICE constexpr inline Matrix22<T>
1423 Matrix22<T>::operator-() const IMATH_NOEXCEPT
1424 {
1425     return Matrix22 (-x[0][0], -x[0][1], -x[1][0], -x[1][1]);
1426 }
1427
1428 template <class T>
1429 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix22<T>&
1430 Matrix22<T>::negate() IMATH_NOEXCEPT
1431 {
1432     x[0][0] = -x[0][0];
1433     x[0][1] = -x[0][1];
1434     x[1][0] = -x[1][0];
1435     x[1][1] = -x[1][1];
1436
1437     return *this;
1438 }
1439
1440 template <class T>
1441 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix22<T>&
1442 Matrix22<T>::operator*= (T a) IMATH_NOEXCEPT
1443 {
1444     x[0][0] *= a;
1445     x[0][1] *= a;
1446     x[1][0] *= a;
1447     x[1][1] *= a;
1448
1449     return *this;
1450 }
1451
1452 template <class T>
1453 IMATH_HOSTDEVICE constexpr inline Matrix22<T>
1454 Matrix22<T>::operator* (T a) const IMATH_NOEXCEPT
1455 {
1456     return Matrix22 (x[0][0] * a, x[0][1] * a, x[1][0] * a, x[1][1] * a);
1457 }
1458
1459 /// Matrix-scalar multiplication
1460 template <class T>
1461 IMATH_HOSTDEVICE inline Matrix22<T>
1462 operator* (T a, const Matrix22<T>& v) IMATH_NOEXCEPT
1463 {
1464     return v * a;
1465 }
1466
1467 template <class T>
1468 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix22<T>&
1469 Matrix22<T>::operator*= (const Matrix22<T>& v) IMATH_NOEXCEPT
1470 {
1471     Matrix22 tmp (T (0));
1472
1473     for (int i = 0; i < 2; i++)
1474         for (int j = 0; j < 2; j++)
1475             for (int k = 0; k < 2; k++)
1476                 tmp.x[i][j] += x[i][k] * v.x[k][j];
1477
1478     *this = tmp;
1479     return *this;
1480 }
1481
1482 template <class T>
1483 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix22<T>
1484 Matrix22<T>::operator* (const Matrix22<T>& v) const IMATH_NOEXCEPT
1485 {
1486     Matrix22 tmp (T (0));
1487
1488     for (int i = 0; i < 2; i++)
1489         for (int j = 0; j < 2; j++)
1490             for (int k = 0; k < 2; k++)
1491                 tmp.x[i][j] += x[i][k] * v.x[k][j];
1492
1493     return tmp;
1494 }
1495
1496 template <class T>
1497 template <class S>
1498 IMATH_HOSTDEVICE inline void
1499 Matrix22<T>::multDirMatrix (const Vec2<S>& src, Vec2<S>& dst) const IMATH_NOEXCEPT
1500 {
1501     S a, b;
1502
1503     a = src.x * x[0][0] + src.y * x[1][0];
1504     b = src.x * x[0][1] + src.y * x[1][1];
1505
1506     dst.x = a;
1507     dst.y = b;
1508 }
1509
1510 template <class T>
1511 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix22<T>&
1512 Matrix22<T>::operator/= (T a) IMATH_NOEXCEPT
1513 {
1514     x[0][0] /= a;
1515     x[0][1] /= a;
1516     x[1][0] /= a;
1517     x[1][1] /= a;
1518
1519     return *this;
1520 }
1521
1522 template <class T>
1523 IMATH_HOSTDEVICE constexpr inline Matrix22<T>
1524 Matrix22<T>::operator/ (T a) const IMATH_NOEXCEPT
1525 {
1526     return Matrix22 (x[0][0] / a, x[0][1] / a, x[1][0] / a, x[1][1] / a);
1527 }
1528
1529 template <class T>
1530 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix22<T>&
1531 Matrix22<T>::transpose() IMATH_NOEXCEPT
1532 {
1533     Matrix22 tmp (x[0][0], x[1][0], x[0][1], x[1][1]);
1534     *this = tmp;
1535     return *this;
1536 }
1537
1538 template <class T>
1539 IMATH_HOSTDEVICE constexpr inline Matrix22<T>
1540 Matrix22<T>::transposed() const IMATH_NOEXCEPT
1541 {
1542     return Matrix22 (x[0][0], x[1][0], x[0][1], x[1][1]);
1543 }
1544
1545 template <class T>
1546 IMATH_CONSTEXPR14 inline const Matrix22<T>&
1547 Matrix22<T>::invert (bool singExc)
1548 {
1549     *this = inverse (singExc);
1550     return *this;
1551 }
1552
1553 template <class T>
1554 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix22<T>&
1555 Matrix22<T>::invert() IMATH_NOEXCEPT
1556 {
1557     *this = inverse();
1558     return *this;
1559 }
1560
1561 template <class T>
1562 IMATH_CONSTEXPR14 inline Matrix22<T>
1563 Matrix22<T>::inverse (bool singExc) const
1564 {
1565     Matrix22 s (x[1][1], -x[0][1], -x[1][0], x[0][0]);
1566
1567     T r = x[0][0] * x[1][1] - x[1][0] * x[0][1];
1568
1569     if (IMATH_INTERNAL_NAMESPACE::abs (r) >= 1)
1570     {
1571         for (int i = 0; i < 2; ++i)
1572         {
1573             for (int j = 0; j < 2; ++j)
1574             {
1575                 s[i][j] /= r;
1576             }
1577         }
1578     }
1579     else
1580     {
1581         T mr = IMATH_INTERNAL_NAMESPACE::abs (r) / std::numeric_limits<T>::min();
1582
1583         for (int i = 0; i < 2; ++i)
1584         {
1585             for (int j = 0; j < 2; ++j)
1586             {
1587                 if (mr > IMATH_INTERNAL_NAMESPACE::abs (s[i][j]))
1588                 {
1589                     s[i][j] /= r;
1590                 }
1591                 else
1592                 {
1593                     if (singExc)
1594                         throw std::invalid_argument ("Cannot invert "
1595                                                      "singular matrix.");
1596                     return Matrix22();
1597                 }
1598             }
1599         }
1600     }
1601     return s;
1602 }
1603
1604 template <class T>
1605 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix22<T>
1606 Matrix22<T>::inverse() const IMATH_NOEXCEPT
1607 {
1608     Matrix22 s (x[1][1], -x[0][1], -x[1][0], x[0][0]);
1609
1610     T r = x[0][0] * x[1][1] - x[1][0] * x[0][1];
1611
1612     if (IMATH_INTERNAL_NAMESPACE::abs (r) >= 1)
1613     {
1614         for (int i = 0; i < 2; ++i)
1615         {
1616             for (int j = 0; j < 2; ++j)
1617             {
1618                 s[i][j] /= r;
1619             }
1620         }
1621     }
1622     else
1623     {
1624         T mr = IMATH_INTERNAL_NAMESPACE::abs (r) / std::numeric_limits<T>::min();
1625
1626         for (int i = 0; i < 2; ++i)
1627         {
1628             for (int j = 0; j < 2; ++j)
1629             {
1630                 if (mr > IMATH_INTERNAL_NAMESPACE::abs (s[i][j]))
1631                 {
1632                     s[i][j] /= r;
1633                 }
1634                 else
1635                 {
1636                     return Matrix22();
1637                 }
1638             }
1639         }
1640     }
1641     return s;
1642 }
1643
1644 template <class T>
1645 IMATH_HOSTDEVICE constexpr inline T
1646 Matrix22<T>::determinant() const IMATH_NOEXCEPT
1647 {
1648     return x[0][0] * x[1][1] - x[1][0] * x[0][1];
1649 }
1650
1651 template <class T>
1652 IMATH_HOSTDEVICE constexpr inline T
1653 Matrix22<T>::trace () const IMATH_NOEXCEPT
1654 {
1655     return x[0][0] + x[1][1];
1656 }
1657
1658 template <class T>
1659 template <class S>
1660 IMATH_HOSTDEVICE inline const Matrix22<T>&
1661 Matrix22<T>::setRotation (S r) IMATH_NOEXCEPT
1662 {
1663     S cos_r, sin_r;
1664
1665     cos_r = cos ((T) r);
1666     sin_r = sin ((T) r);
1667
1668     x[0][0] = cos_r;
1669     x[0][1] = sin_r;
1670
1671     x[1][0] = -sin_r;
1672     x[1][1] = cos_r;
1673
1674     return *this;
1675 }
1676
1677 template <class T>
1678 template <class S>
1679 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix22<T>&
1680 Matrix22<T>::rotate (S r) IMATH_NOEXCEPT
1681 {
1682     *this *= Matrix22<T>().setRotation (r);
1683     return *this;
1684 }
1685
1686 template <class T>
1687 IMATH_CONSTEXPR14 inline const Matrix22<T>&
1688 Matrix22<T>::setScale (T s) IMATH_NOEXCEPT
1689 {
1690     //
1691     // Set the matrix to:
1692     //  | s 0 |
1693     //  | 0 s |
1694     //
1695
1696     x[0][0] = s;
1697     x[0][1] = static_cast<T> (0);
1698     x[1][0] = static_cast<T> (0);
1699     x[1][1] = s;
1700
1701     return *this;
1702 }
1703
1704 template <class T>
1705 template <class S>
1706 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix22<T>&
1707 Matrix22<T>::setScale (const Vec2<S>& s) IMATH_NOEXCEPT
1708 {
1709     //
1710     // Set the matrix to:
1711     //  | s.x  0  |
1712     //  |  0  s.y |
1713     //
1714
1715     x[0][0] = s.x;
1716     x[0][1] = static_cast<T> (0);
1717     x[1][0] = static_cast<T> (0);
1718     x[1][1] = s.y;
1719
1720     return *this;
1721 }
1722
1723 template <class T>
1724 template <class S>
1725 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix22<T>&
1726 Matrix22<T>::scale (const Vec2<S>& s) IMATH_NOEXCEPT
1727 {
1728     x[0][0] *= s.x;
1729     x[0][1] *= s.x;
1730
1731     x[1][0] *= s.y;
1732     x[1][1] *= s.y;
1733
1734     return *this;
1735 }
1736
1737 //---------------------------
1738 // Implementation of Matrix33
1739 //---------------------------
1740
1741 template <class T>
1742 IMATH_HOSTDEVICE inline T*
1743 Matrix33<T>::operator[] (int i) IMATH_NOEXCEPT
1744 {
1745     return x[i];
1746 }
1747
1748 template <class T>
1749 IMATH_HOSTDEVICE inline const T*
1750 Matrix33<T>::operator[] (int i) const IMATH_NOEXCEPT
1751 {
1752     return x[i];
1753 }
1754
1755 template <class T>
1756 IMATH_HOSTDEVICE inline IMATH_CONSTEXPR14
1757 Matrix33<T>::Matrix33() IMATH_NOEXCEPT
1758 {
1759     x[0][0] = 1;
1760     x[0][1] = 0;
1761     x[0][2] = 0;
1762     x[1][0] = 0;
1763     x[1][1] = 1;
1764     x[1][2] = 0;
1765     x[2][0] = 0;
1766     x[2][1] = 0;
1767     x[2][2] = 1;
1768 }
1769
1770 template <class T> IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix33<T>::Matrix33 (T a) IMATH_NOEXCEPT
1771 {
1772     x[0][0] = a;
1773     x[0][1] = a;
1774     x[0][2] = a;
1775     x[1][0] = a;
1776     x[1][1] = a;
1777     x[1][2] = a;
1778     x[2][0] = a;
1779     x[2][1] = a;
1780     x[2][2] = a;
1781 }
1782
1783 template <class T>
1784 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix33<T>::Matrix33 (
1785     const T a[3][3]) IMATH_NOEXCEPT
1786 {
1787     // Function calls and aliasing issues can inhibit vectorization versus
1788     // straight assignment of data members, so instead of this:
1789     //     memcpy (x, a, sizeof (x));
1790     // we do this:
1791     x[0][0] = a[0][0];
1792     x[0][1] = a[0][1];
1793     x[0][2] = a[0][2];
1794     x[1][0] = a[1][0];
1795     x[1][1] = a[1][1];
1796     x[1][2] = a[1][2];
1797     x[2][0] = a[2][0];
1798     x[2][1] = a[2][1];
1799     x[2][2] = a[2][2];
1800 }
1801
1802 template <class T>
1803 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix33<T>::Matrix33 (T a, T b, T c, T d, T e, T f, T g, T h, T i) IMATH_NOEXCEPT
1804 {
1805     x[0][0] = a;
1806     x[0][1] = b;
1807     x[0][2] = c;
1808     x[1][0] = d;
1809     x[1][1] = e;
1810     x[1][2] = f;
1811     x[2][0] = g;
1812     x[2][1] = h;
1813     x[2][2] = i;
1814 }
1815
1816 template <class T>
1817 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix33<T>::Matrix33 (const Matrix33& v) IMATH_NOEXCEPT
1818 {
1819     // Function calls and aliasing issues can inhibit vectorization versus
1820     // straight assignment of data members, so instead of this:
1821     //     memcpy (x, v.x, sizeof (x));
1822     // we do this:
1823     x[0][0] = v.x[0][0];
1824     x[0][1] = v.x[0][1];
1825     x[0][2] = v.x[0][2];
1826     x[1][0] = v.x[1][0];
1827     x[1][1] = v.x[1][1];
1828     x[1][2] = v.x[1][2];
1829     x[2][0] = v.x[2][0];
1830     x[2][1] = v.x[2][1];
1831     x[2][2] = v.x[2][2];
1832 }
1833
1834 template <class T>
1835 template <class S>
1836 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix33<T>::Matrix33 (const Matrix33<S>& v) IMATH_NOEXCEPT
1837 {
1838     x[0][0] = T (v.x[0][0]);
1839     x[0][1] = T (v.x[0][1]);
1840     x[0][2] = T (v.x[0][2]);
1841     x[1][0] = T (v.x[1][0]);
1842     x[1][1] = T (v.x[1][1]);
1843     x[1][2] = T (v.x[1][2]);
1844     x[2][0] = T (v.x[2][0]);
1845     x[2][1] = T (v.x[2][1]);
1846     x[2][2] = T (v.x[2][2]);
1847 }
1848
1849 template <class T>
1850 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix33<T>&
1851 Matrix33<T>::operator= (const Matrix33& v) IMATH_NOEXCEPT
1852 {
1853     // Function calls and aliasing issues can inhibit vectorization versus
1854     // straight assignment of data members, so instead of this:
1855     //     memcpy (x, v.x, sizeof (x));
1856     // we do this:
1857     x[0][0] = v.x[0][0];
1858     x[0][1] = v.x[0][1];
1859     x[0][2] = v.x[0][2];
1860     x[1][0] = v.x[1][0];
1861     x[1][1] = v.x[1][1];
1862     x[1][2] = v.x[1][2];
1863     x[2][0] = v.x[2][0];
1864     x[2][1] = v.x[2][1];
1865     x[2][2] = v.x[2][2];
1866     return *this;
1867 }
1868
1869 template <class T>
1870 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix33<T>&
1871 Matrix33<T>::operator= (T a) IMATH_NOEXCEPT
1872 {
1873     x[0][0] = a;
1874     x[0][1] = a;
1875     x[0][2] = a;
1876     x[1][0] = a;
1877     x[1][1] = a;
1878     x[1][2] = a;
1879     x[2][0] = a;
1880     x[2][1] = a;
1881     x[2][2] = a;
1882     return *this;
1883 }
1884
1885 template <class T>
1886 IMATH_HOSTDEVICE inline T*
1887 Matrix33<T>::getValue () IMATH_NOEXCEPT
1888 {
1889     return (T*) &x[0][0];
1890 }
1891
1892 template <class T>
1893 IMATH_HOSTDEVICE inline const T*
1894 Matrix33<T>::getValue() const IMATH_NOEXCEPT
1895 {
1896     return (const T*) &x[0][0];
1897 }
1898
1899 template <class T>
1900 template <class S>
1901 IMATH_HOSTDEVICE inline void
1902 Matrix33<T>::getValue (Matrix33<S>& v) const IMATH_NOEXCEPT
1903 {
1904     v.x[0][0] = x[0][0];
1905     v.x[0][1] = x[0][1];
1906     v.x[0][2] = x[0][2];
1907     v.x[1][0] = x[1][0];
1908     v.x[1][1] = x[1][1];
1909     v.x[1][2] = x[1][2];
1910     v.x[2][0] = x[2][0];
1911     v.x[2][1] = x[2][1];
1912     v.x[2][2] = x[2][2];
1913 }
1914
1915 template <class T>
1916 template <class S>
1917 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix33<T>&
1918 Matrix33<T>::setValue (const Matrix33<S>& v) IMATH_NOEXCEPT
1919 {
1920     x[0][0] = v.x[0][0];
1921     x[0][1] = v.x[0][1];
1922     x[0][2] = v.x[0][2];
1923     x[1][0] = v.x[1][0];
1924     x[1][1] = v.x[1][1];
1925     x[1][2] = v.x[1][2];
1926     x[2][0] = v.x[2][0];
1927     x[2][1] = v.x[2][1];
1928     x[2][2] = v.x[2][2];
1929     return *this;
1930 }
1931
1932 template <class T>
1933 template <class S>
1934 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix33<T>&
1935 Matrix33<T>::setTheMatrix (const Matrix33<S>& v) IMATH_NOEXCEPT
1936 {
1937     x[0][0] = v.x[0][0];
1938     x[0][1] = v.x[0][1];
1939     x[0][2] = v.x[0][2];
1940     x[1][0] = v.x[1][0];
1941     x[1][1] = v.x[1][1];
1942     x[1][2] = v.x[1][2];
1943     x[2][0] = v.x[2][0];
1944     x[2][1] = v.x[2][1];
1945     x[2][2] = v.x[2][2];
1946     return *this;
1947 }
1948
1949 template <class T>
1950 IMATH_HOSTDEVICE inline void
1951 Matrix33<T>::makeIdentity() IMATH_NOEXCEPT
1952 {
1953     x[0][0] = 1;
1954     x[0][1] = 0;
1955     x[0][2] = 0;
1956     x[1][0] = 0;
1957     x[1][1] = 1;
1958     x[1][2] = 0;
1959     x[2][0] = 0;
1960     x[2][1] = 0;
1961     x[2][2] = 1;
1962 }
1963
1964 template <class T>
1965 IMATH_HOSTDEVICE constexpr inline bool
1966 Matrix33<T>::operator== (const Matrix33& v) const IMATH_NOEXCEPT
1967 {
1968     return x[0][0] == v.x[0][0] && x[0][1] == v.x[0][1] && x[0][2] == v.x[0][2] &&
1969            x[1][0] == v.x[1][0] && x[1][1] == v.x[1][1] && x[1][2] == v.x[1][2] &&
1970            x[2][0] == v.x[2][0] && x[2][1] == v.x[2][1] && x[2][2] == v.x[2][2];
1971 }
1972
1973 template <class T>
1974 IMATH_HOSTDEVICE constexpr inline bool
1975 Matrix33<T>::operator!= (const Matrix33& v) const IMATH_NOEXCEPT
1976 {
1977     return x[0][0] != v.x[0][0] || x[0][1] != v.x[0][1] || x[0][2] != v.x[0][2] ||
1978            x[1][0] != v.x[1][0] || x[1][1] != v.x[1][1] || x[1][2] != v.x[1][2] ||
1979            x[2][0] != v.x[2][0] || x[2][1] != v.x[2][1] || x[2][2] != v.x[2][2];
1980 }
1981
1982 template <class T>
1983 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline bool
1984 Matrix33<T>::equalWithAbsError (const Matrix33<T>& m, T e) const IMATH_NOEXCEPT
1985 {
1986     for (int i = 0; i < 3; i++)
1987         for (int j = 0; j < 3; j++)
1988             if (!IMATH_INTERNAL_NAMESPACE::equalWithAbsError ((*this)[i][j], m[i][j], e))
1989                 return false;
1990
1991     return true;
1992 }
1993
1994 template <class T>
1995 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline bool
1996 Matrix33<T>::equalWithRelError (const Matrix33<T>& m, T e) const IMATH_NOEXCEPT
1997 {
1998     for (int i = 0; i < 3; i++)
1999         for (int j = 0; j < 3; j++)
2000             if (!IMATH_INTERNAL_NAMESPACE::equalWithRelError ((*this)[i][j], m[i][j], e))
2001                 return false;
2002
2003     return true;
2004 }
2005
2006 template <class T>
2007 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix33<T>&
2008 Matrix33<T>::operator+= (const Matrix33<T>& v) IMATH_NOEXCEPT
2009 {
2010     x[0][0] += v.x[0][0];
2011     x[0][1] += v.x[0][1];
2012     x[0][2] += v.x[0][2];
2013     x[1][0] += v.x[1][0];
2014     x[1][1] += v.x[1][1];
2015     x[1][2] += v.x[1][2];
2016     x[2][0] += v.x[2][0];
2017     x[2][1] += v.x[2][1];
2018     x[2][2] += v.x[2][2];
2019
2020     return *this;
2021 }
2022
2023 template <class T>
2024 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix33<T>&
2025 Matrix33<T>::operator+= (T a) IMATH_NOEXCEPT
2026 {
2027     x[0][0] += a;
2028     x[0][1] += a;
2029     x[0][2] += a;
2030     x[1][0] += a;
2031     x[1][1] += a;
2032     x[1][2] += a;
2033     x[2][0] += a;
2034     x[2][1] += a;
2035     x[2][2] += a;
2036
2037     return *this;
2038 }
2039
2040 template <class T>
2041 IMATH_HOSTDEVICE constexpr inline Matrix33<T>
2042 Matrix33<T>::operator+ (const Matrix33<T>& v) const IMATH_NOEXCEPT
2043 {
2044     return Matrix33 (x[0][0] + v.x[0][0],
2045                      x[0][1] + v.x[0][1],
2046                      x[0][2] + v.x[0][2],
2047                      x[1][0] + v.x[1][0],
2048                      x[1][1] + v.x[1][1],
2049                      x[1][2] + v.x[1][2],
2050                      x[2][0] + v.x[2][0],
2051                      x[2][1] + v.x[2][1],
2052                      x[2][2] + v.x[2][2]);
2053 }
2054
2055 template <class T>
2056 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix33<T>&
2057 Matrix33<T>::operator-= (const Matrix33<T>& v) IMATH_NOEXCEPT
2058 {
2059     x[0][0] -= v.x[0][0];
2060     x[0][1] -= v.x[0][1];
2061     x[0][2] -= v.x[0][2];
2062     x[1][0] -= v.x[1][0];
2063     x[1][1] -= v.x[1][1];
2064     x[1][2] -= v.x[1][2];
2065     x[2][0] -= v.x[2][0];
2066     x[2][1] -= v.x[2][1];
2067     x[2][2] -= v.x[2][2];
2068
2069     return *this;
2070 }
2071
2072 template <class T>
2073 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix33<T>&
2074 Matrix33<T>::operator-= (T a) IMATH_NOEXCEPT
2075 {
2076     x[0][0] -= a;
2077     x[0][1] -= a;
2078     x[0][2] -= a;
2079     x[1][0] -= a;
2080     x[1][1] -= a;
2081     x[1][2] -= a;
2082     x[2][0] -= a;
2083     x[2][1] -= a;
2084     x[2][2] -= a;
2085
2086     return *this;
2087 }
2088
2089 template <class T>
2090 IMATH_HOSTDEVICE constexpr inline Matrix33<T>
2091 Matrix33<T>::operator- (const Matrix33<T>& v) const IMATH_NOEXCEPT
2092 {
2093     return Matrix33 (x[0][0] - v.x[0][0],
2094                      x[0][1] - v.x[0][1],
2095                      x[0][2] - v.x[0][2],
2096                      x[1][0] - v.x[1][0],
2097                      x[1][1] - v.x[1][1],
2098                      x[1][2] - v.x[1][2],
2099                      x[2][0] - v.x[2][0],
2100                      x[2][1] - v.x[2][1],
2101                      x[2][2] - v.x[2][2]);
2102 }
2103
2104 template <class T>
2105 IMATH_HOSTDEVICE constexpr inline Matrix33<T>
2106 Matrix33<T>::operator-() const IMATH_NOEXCEPT
2107 {
2108     return Matrix33 (-x[0][0],
2109                      -x[0][1],
2110                      -x[0][2],
2111                      -x[1][0],
2112                      -x[1][1],
2113                      -x[1][2],
2114                      -x[2][0],
2115                      -x[2][1],
2116                      -x[2][2]);
2117 }
2118
2119 template <class T>
2120 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix33<T>&
2121 Matrix33<T>::negate() IMATH_NOEXCEPT
2122 {
2123     x[0][0] = -x[0][0];
2124     x[0][1] = -x[0][1];
2125     x[0][2] = -x[0][2];
2126     x[1][0] = -x[1][0];
2127     x[1][1] = -x[1][1];
2128     x[1][2] = -x[1][2];
2129     x[2][0] = -x[2][0];
2130     x[2][1] = -x[2][1];
2131     x[2][2] = -x[2][2];
2132
2133     return *this;
2134 }
2135
2136 template <class T>
2137 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix33<T>&
2138 Matrix33<T>::operator*= (T a) IMATH_NOEXCEPT
2139 {
2140     x[0][0] *= a;
2141     x[0][1] *= a;
2142     x[0][2] *= a;
2143     x[1][0] *= a;
2144     x[1][1] *= a;
2145     x[1][2] *= a;
2146     x[2][0] *= a;
2147     x[2][1] *= a;
2148     x[2][2] *= a;
2149
2150     return *this;
2151 }
2152
2153 template <class T>
2154 IMATH_HOSTDEVICE constexpr inline Matrix33<T>
2155 Matrix33<T>::operator* (T a) const IMATH_NOEXCEPT
2156 {
2157     return Matrix33 (x[0][0] * a,
2158                      x[0][1] * a,
2159                      x[0][2] * a,
2160                      x[1][0] * a,
2161                      x[1][1] * a,
2162                      x[1][2] * a,
2163                      x[2][0] * a,
2164                      x[2][1] * a,
2165                      x[2][2] * a);
2166 }
2167
2168 /// Matrix-scalar multiplication
2169 template <class T>
2170 IMATH_HOSTDEVICE inline Matrix33<T> constexpr
2171 operator* (T a, const Matrix33<T>& v) IMATH_NOEXCEPT
2172 {
2173     return v * a;
2174 }
2175
2176 template <class T>
2177 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix33<T>&
2178 Matrix33<T>::operator*= (const Matrix33<T>& v) IMATH_NOEXCEPT
2179 {
2180     // Avoid initializing with 0 values before immediately overwriting them,
2181     // and unroll all loops for the best autovectorization.
2182     Matrix33 tmp(IMATH_INTERNAL_NAMESPACE::UNINITIALIZED);
2183
2184     tmp.x[0][0] = x[0][0] * v.x[0][0] + x[0][1] * v.x[1][0] + x[0][2] * v.x[2][0];
2185     tmp.x[0][1] = x[0][0] * v.x[0][1] + x[0][1] * v.x[1][1] + x[0][2] * v.x[2][1];
2186     tmp.x[0][2] = x[0][0] * v.x[0][2] + x[0][1] * v.x[1][2] + x[0][2] * v.x[2][2];
2187
2188     tmp.x[1][0] = x[1][0] * v.x[0][0] + x[1][1] * v.x[1][0] + x[1][2] * v.x[2][0];
2189     tmp.x[1][1] = x[1][0] * v.x[0][1] + x[1][1] * v.x[1][1] + x[1][2] * v.x[2][1];
2190     tmp.x[1][2] = x[1][0] * v.x[0][2] + x[1][1] * v.x[1][2] + x[1][2] * v.x[2][2];
2191
2192     tmp.x[2][0] = x[2][0] * v.x[0][0] + x[2][1] * v.x[1][0] + x[2][2] * v.x[2][0];
2193     tmp.x[2][1] = x[2][0] * v.x[0][1] + x[2][1] * v.x[1][1] + x[2][2] * v.x[2][1];
2194     tmp.x[2][2] = x[2][0] * v.x[0][2] + x[2][1] * v.x[1][2] + x[2][2] * v.x[2][2];
2195
2196     *this = tmp;
2197     return *this;
2198 }
2199
2200 template <class T>
2201 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix33<T>
2202 Matrix33<T>::operator* (const Matrix33<T>& v) const IMATH_NOEXCEPT
2203 {
2204     // Avoid initializing with 0 values before immediately overwriting them,
2205     // and unroll all loops for the best autovectorization.
2206     Matrix33 tmp(IMATH_INTERNAL_NAMESPACE::UNINITIALIZED);
2207
2208     tmp.x[0][0] = x[0][0] * v.x[0][0] + x[0][1] * v.x[1][0] + x[0][2] * v.x[2][0];
2209     tmp.x[0][1] = x[0][0] * v.x[0][1] + x[0][1] * v.x[1][1] + x[0][2] * v.x[2][1];
2210     tmp.x[0][2] = x[0][0] * v.x[0][2] + x[0][1] * v.x[1][2] + x[0][2] * v.x[2][2];
2211
2212     tmp.x[1][0] = x[1][0] * v.x[0][0] + x[1][1] * v.x[1][0] + x[1][2] * v.x[2][0];
2213     tmp.x[1][1] = x[1][0] * v.x[0][1] + x[1][1] * v.x[1][1] + x[1][2] * v.x[2][1];
2214     tmp.x[1][2] = x[1][0] * v.x[0][2] + x[1][1] * v.x[1][2] + x[1][2] * v.x[2][2];
2215
2216     tmp.x[2][0] = x[2][0] * v.x[0][0] + x[2][1] * v.x[1][0] + x[2][2] * v.x[2][0];
2217     tmp.x[2][1] = x[2][0] * v.x[0][1] + x[2][1] * v.x[1][1] + x[2][2] * v.x[2][1];
2218     tmp.x[2][2] = x[2][0] * v.x[0][2] + x[2][1] * v.x[1][2] + x[2][2] * v.x[2][2];
2219
2220     return tmp;
2221 }
2222
2223 template <class T>
2224 template <class S>
2225 IMATH_HOSTDEVICE inline void
2226 Matrix33<T>::multVecMatrix (const Vec2<S>& src, Vec2<S>& dst) const IMATH_NOEXCEPT
2227 {
2228     S a, b, w;
2229
2230     a = src.x * x[0][0] + src.y * x[1][0] + x[2][0];
2231     b = src.x * x[0][1] + src.y * x[1][1] + x[2][1];
2232     w = src.x * x[0][2] + src.y * x[1][2] + x[2][2];
2233
2234     dst.x = a / w;
2235     dst.y = b / w;
2236 }
2237
2238 template <class T>
2239 template <class S>
2240 IMATH_HOSTDEVICE inline void
2241 Matrix33<T>::multDirMatrix (const Vec2<S>& src, Vec2<S>& dst) const IMATH_NOEXCEPT
2242 {
2243     S a, b;
2244
2245     a = src.x * x[0][0] + src.y * x[1][0];
2246     b = src.x * x[0][1] + src.y * x[1][1];
2247
2248     dst.x = a;
2249     dst.y = b;
2250 }
2251
2252 template <class T>
2253 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix33<T>&
2254 Matrix33<T>::operator/= (T a) IMATH_NOEXCEPT
2255 {
2256     x[0][0] /= a;
2257     x[0][1] /= a;
2258     x[0][2] /= a;
2259     x[1][0] /= a;
2260     x[1][1] /= a;
2261     x[1][2] /= a;
2262     x[2][0] /= a;
2263     x[2][1] /= a;
2264     x[2][2] /= a;
2265
2266     return *this;
2267 }
2268
2269 template <class T>
2270 IMATH_HOSTDEVICE constexpr inline Matrix33<T>
2271 Matrix33<T>::operator/ (T a) const IMATH_NOEXCEPT
2272 {
2273     return Matrix33 (x[0][0] / a,
2274                      x[0][1] / a,
2275                      x[0][2] / a,
2276                      x[1][0] / a,
2277                      x[1][1] / a,
2278                      x[1][2] / a,
2279                      x[2][0] / a,
2280                      x[2][1] / a,
2281                      x[2][2] / a);
2282 }
2283
2284 template <class T>
2285 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix33<T>&
2286 Matrix33<T>::transpose() IMATH_NOEXCEPT
2287 {
2288     Matrix33 tmp (x[0][0], x[1][0], x[2][0], x[0][1], x[1][1], x[2][1], x[0][2], x[1][2], x[2][2]);
2289     *this = tmp;
2290     return *this;
2291 }
2292
2293 template <class T>
2294 IMATH_HOSTDEVICE constexpr inline Matrix33<T>
2295 Matrix33<T>::transposed() const IMATH_NOEXCEPT
2296 {
2297     return Matrix33 (x[0][0],
2298                      x[1][0],
2299                      x[2][0],
2300                      x[0][1],
2301                      x[1][1],
2302                      x[2][1],
2303                      x[0][2],
2304                      x[1][2],
2305                      x[2][2]);
2306 }
2307
2308 template <class T>
2309 const inline Matrix33<T>&
2310 Matrix33<T>::gjInvert (bool singExc)
2311 {
2312     *this = gjInverse (singExc);
2313     return *this;
2314 }
2315
2316 template <class T>
2317 IMATH_HOSTDEVICE const inline Matrix33<T>&
2318 Matrix33<T>::gjInvert() IMATH_NOEXCEPT
2319 {
2320     *this = gjInverse();
2321     return *this;
2322 }
2323
2324 template <class T>
2325 inline Matrix33<T>
2326 Matrix33<T>::gjInverse (bool singExc) const
2327 {
2328     int i, j, k;
2329     Matrix33 s;
2330     Matrix33 t (*this);
2331
2332     // Forward elimination
2333
2334     for (i = 0; i < 2; i++)
2335     {
2336         int pivot = i;
2337
2338         T pivotsize = t.x[i][i];
2339
2340         if (pivotsize < 0)
2341             pivotsize = -pivotsize;
2342
2343         for (j = i + 1; j < 3; j++)
2344         {
2345             T tmp = t.x[j][i];
2346
2347             if (tmp < 0)
2348                 tmp = -tmp;
2349
2350             if (tmp > pivotsize)
2351             {
2352                 pivot     = j;
2353                 pivotsize = tmp;
2354             }
2355         }
2356
2357         if (pivotsize == 0)
2358         {
2359             if (singExc)
2360                 throw std::invalid_argument ("Cannot invert singular matrix.");
2361
2362             return Matrix33();
2363         }
2364
2365         if (pivot != i)
2366         {
2367             for (j = 0; j < 3; j++)
2368             {
2369                 T tmp;
2370
2371                 tmp           = t.x[i][j];
2372                 t.x[i][j]     = t.x[pivot][j];
2373                 t.x[pivot][j] = tmp;
2374
2375                 tmp           = s.x[i][j];
2376                 s.x[i][j]     = s.x[pivot][j];
2377                 s.x[pivot][j] = tmp;
2378             }
2379         }
2380
2381         for (j = i + 1; j < 3; j++)
2382         {
2383             T f = t.x[j][i] / t.x[i][i];
2384
2385             for (k = 0; k < 3; k++)
2386             {
2387                 t.x[j][k] -= f * t.x[i][k];
2388                 s.x[j][k] -= f * s.x[i][k];
2389             }
2390         }
2391     }
2392
2393     // Backward substitution
2394
2395     for (i = 2; i >= 0; --i)
2396     {
2397         T f;
2398
2399         if ((f = t[i][i]) == 0)
2400         {
2401             if (singExc)
2402                 throw std::invalid_argument ("Cannot invert singular matrix.");
2403
2404             return Matrix33();
2405         }
2406
2407         for (j = 0; j < 3; j++)
2408         {
2409             t.x[i][j] /= f;
2410             s.x[i][j] /= f;
2411         }
2412
2413         for (j = 0; j < i; j++)
2414         {
2415             f = t.x[j][i];
2416
2417             for (k = 0; k < 3; k++)
2418             {
2419                 t.x[j][k] -= f * t.x[i][k];
2420                 s.x[j][k] -= f * s.x[i][k];
2421             }
2422         }
2423     }
2424
2425     return s;
2426 }
2427
2428 template <class T>
2429 IMATH_HOSTDEVICE inline Matrix33<T>
2430 Matrix33<T>::gjInverse() const IMATH_NOEXCEPT
2431 {
2432     int i, j, k;
2433     Matrix33 s;
2434     Matrix33 t (*this);
2435
2436     // Forward elimination
2437
2438     for (i = 0; i < 2; i++)
2439     {
2440         int pivot = i;
2441
2442         T pivotsize = t.x[i][i];
2443
2444         if (pivotsize < 0)
2445             pivotsize = -pivotsize;
2446
2447         for (j = i + 1; j < 3; j++)
2448         {
2449             T tmp = t.x[j][i];
2450
2451             if (tmp < 0)
2452                 tmp = -tmp;
2453
2454             if (tmp > pivotsize)
2455             {
2456                 pivot     = j;
2457                 pivotsize = tmp;
2458             }
2459         }
2460
2461         if (pivotsize == 0)
2462         {
2463             return Matrix33();
2464         }
2465
2466         if (pivot != i)
2467         {
2468             for (j = 0; j < 3; j++)
2469             {
2470                 T tmp;
2471
2472                 tmp           = t.x[i][j];
2473                 t.x[i][j]     = t.x[pivot][j];
2474                 t.x[pivot][j] = tmp;
2475
2476                 tmp           = s.x[i][j];
2477                 s.x[i][j]     = s.x[pivot][j];
2478                 s.x[pivot][j] = tmp;
2479             }
2480         }
2481
2482         for (j = i + 1; j < 3; j++)
2483         {
2484             T f = t.x[j][i] / t.x[i][i];
2485
2486             for (k = 0; k < 3; k++)
2487             {
2488                 t.x[j][k] -= f * t.x[i][k];
2489                 s.x[j][k] -= f * s.x[i][k];
2490             }
2491         }
2492     }
2493
2494     // Backward substitution
2495
2496     for (i = 2; i >= 0; --i)
2497     {
2498         T f;
2499
2500         if ((f = t.x[i][i]) == 0)
2501         {
2502             return Matrix33();
2503         }
2504
2505         for (j = 0; j < 3; j++)
2506         {
2507             t.x[i][j] /= f;
2508             s.x[i][j] /= f;
2509         }
2510
2511         for (j = 0; j < i; j++)
2512         {
2513             f = t.x[j][i];
2514
2515             for (k = 0; k < 3; k++)
2516             {
2517                 t.x[j][k] -= f * t.x[i][k];
2518                 s.x[j][k] -= f * s.x[i][k];
2519             }
2520         }
2521     }
2522
2523     return s;
2524 }
2525
2526 template <class T>
2527 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix33<T>&
2528 Matrix33<T>::invert (bool singExc)
2529 {
2530     *this = inverse (singExc);
2531     return *this;
2532 }
2533
2534 template <class T>
2535 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix33<T>&
2536 Matrix33<T>::invert() IMATH_NOEXCEPT
2537 {
2538     *this = inverse();
2539     return *this;
2540 }
2541
2542 template <class T>
2543 IMATH_CONSTEXPR14 inline Matrix33<T>
2544 Matrix33<T>::inverse (bool singExc) const
2545 {
2546     if (x[0][2] != 0 || x[1][2] != 0 || x[2][2] != 1)
2547     {
2548         Matrix33 s (x[1][1] * x[2][2] - x[2][1] * x[1][2],
2549                     x[2][1] * x[0][2] - x[0][1] * x[2][2],
2550                     x[0][1] * x[1][2] - x[1][1] * x[0][2],
2551
2552                     x[2][0] * x[1][2] - x[1][0] * x[2][2],
2553                     x[0][0] * x[2][2] - x[2][0] * x[0][2],
2554                     x[1][0] * x[0][2] - x[0][0] * x[1][2],
2555
2556                     x[1][0] * x[2][1] - x[2][0] * x[1][1],
2557                     x[2][0] * x[0][1] - x[0][0] * x[2][1],
2558                     x[0][0] * x[1][1] - x[1][0] * x[0][1]);
2559
2560         T r = x[0][0] * s[0][0] + x[0][1] * s[1][0] + x[0][2] * s[2][0];
2561
2562         if (IMATH_INTERNAL_NAMESPACE::abs (r) >= 1)
2563         {
2564             for (int i = 0; i < 3; ++i)
2565             {
2566                 for (int j = 0; j < 3; ++j)
2567                 {
2568                     s.x[i][j] /= r;
2569                 }
2570             }
2571         }
2572         else
2573         {
2574             T mr = IMATH_INTERNAL_NAMESPACE::abs (r) / std::numeric_limits<T>::min();
2575
2576             for (int i = 0; i < 3; ++i)
2577             {
2578                 for (int j = 0; j < 3; ++j)
2579                 {
2580                     if (mr > IMATH_INTERNAL_NAMESPACE::abs (s.x[i][j]))
2581                     {
2582                         s.x[i][j] /= r;
2583                     }
2584                     else
2585                     {
2586                         if (singExc)
2587                             throw std::invalid_argument ("Cannot invert "
2588                                                          "singular matrix.");
2589                         return Matrix33();
2590                     }
2591                 }
2592             }
2593         }
2594
2595         return s;
2596     }
2597     else
2598     {
2599         Matrix33 s (x[1][1],
2600                     -x[0][1],
2601                     0,
2602
2603                     -x[1][0],
2604                     x[0][0],
2605                     0,
2606
2607                     0,
2608                     0,
2609                     1);
2610
2611         T r = x[0][0] * x[1][1] - x[1][0] * x[0][1];
2612
2613         if (IMATH_INTERNAL_NAMESPACE::abs (r) >= 1)
2614         {
2615             for (int i = 0; i < 2; ++i)
2616             {
2617                 for (int j = 0; j < 2; ++j)
2618                 {
2619                     s.x[i][j] /= r;
2620                 }
2621             }
2622         }
2623         else
2624         {
2625             T mr = IMATH_INTERNAL_NAMESPACE::abs (r) / std::numeric_limits<T>::min();
2626
2627             for (int i = 0; i < 2; ++i)
2628             {
2629                 for (int j = 0; j < 2; ++j)
2630                 {
2631                     if (mr > IMATH_INTERNAL_NAMESPACE::abs (s.x[i][j]))
2632                     {
2633                         s.x[i][j] /= r;
2634                     }
2635                     else
2636                     {
2637                         if (singExc)
2638                             throw std::invalid_argument ("Cannot invert "
2639                                                          "singular matrix.");
2640                         return Matrix33();
2641                     }
2642                 }
2643             }
2644         }
2645
2646         s.x[2][0] = -x[2][0] * s.x[0][0] - x[2][1] * s.x[1][0];
2647         s.x[2][1] = -x[2][0] * s.x[0][1] - x[2][1] * s.x[1][1];
2648
2649         return s;
2650     }
2651 }
2652
2653 template <class T>
2654 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix33<T>
2655 Matrix33<T>::inverse () const IMATH_NOEXCEPT
2656 {
2657     if (x[0][2] != 0 || x[1][2] != 0 || x[2][2] != 1)
2658     {
2659         Matrix33 s (x[1][1] * x[2][2] - x[2][1] * x[1][2],
2660                     x[2][1] * x[0][2] - x[0][1] * x[2][2],
2661                     x[0][1] * x[1][2] - x[1][1] * x[0][2],
2662
2663                     x[2][0] * x[1][2] - x[1][0] * x[2][2],
2664                     x[0][0] * x[2][2] - x[2][0] * x[0][2],
2665                     x[1][0] * x[0][2] - x[0][0] * x[1][2],
2666
2667                     x[1][0] * x[2][1] - x[2][0] * x[1][1],
2668                     x[2][0] * x[0][1] - x[0][0] * x[2][1],
2669                     x[0][0] * x[1][1] - x[1][0] * x[0][1]);
2670
2671         T r = x[0][0] * s.x[0][0] + x[0][1] * s.x[1][0] + x[0][2] * s.x[2][0];
2672
2673         if (IMATH_INTERNAL_NAMESPACE::abs (r) >= 1)
2674         {
2675             for (int i = 0; i < 3; ++i)
2676             {
2677                 for (int j = 0; j < 3; ++j)
2678                 {
2679                     s.x[i][j] /= r;
2680                 }
2681             }
2682         }
2683         else
2684         {
2685             T mr = IMATH_INTERNAL_NAMESPACE::abs (r) / std::numeric_limits<T>::min();
2686
2687             for (int i = 0; i < 3; ++i)
2688             {
2689                 for (int j = 0; j < 3; ++j)
2690                 {
2691                     if (mr > IMATH_INTERNAL_NAMESPACE::abs (s.x[i][j]))
2692                     {
2693                         s.x[i][j] /= r;
2694                     }
2695                     else
2696                     {
2697                         return Matrix33();
2698                     }
2699                 }
2700             }
2701         }
2702
2703         return s;
2704     }
2705     else
2706     {
2707         Matrix33 s (x[1][1],
2708                     -x[0][1],
2709                     0,
2710
2711                     -x[1][0],
2712                     x[0][0],
2713                     0,
2714
2715                     0,
2716                     0,
2717                     1);
2718
2719         T r = x[0][0] * x[1][1] - x[1][0] * x[0][1];
2720
2721         if (IMATH_INTERNAL_NAMESPACE::abs (r) >= 1)
2722         {
2723             for (int i = 0; i < 2; ++i)
2724             {
2725                 for (int j = 0; j < 2; ++j)
2726                 {
2727                     s.x[i][j] /= r;
2728                 }
2729             }
2730         }
2731         else
2732         {
2733             T mr = IMATH_INTERNAL_NAMESPACE::abs (r) / std::numeric_limits<T>::min();
2734
2735             for (int i = 0; i < 2; ++i)
2736             {
2737                 for (int j = 0; j < 2; ++j)
2738                 {
2739                     if (mr > IMATH_INTERNAL_NAMESPACE::abs (s.x[i][j]))
2740                     {
2741                         s.x[i][j] /= r;
2742                     }
2743                     else
2744                     {
2745                         return Matrix33();
2746                     }
2747                 }
2748             }
2749         }
2750
2751         s.x[2][0] = -x[2][0] * s.x[0][0] - x[2][1] * s.x[1][0];
2752         s.x[2][1] = -x[2][0] * s.x[0][1] - x[2][1] * s.x[1][1];
2753
2754         return s;
2755     }
2756 }
2757
2758 template <class T>
2759 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline T
2760 Matrix33<T>::minorOf (const int r, const int c) const IMATH_NOEXCEPT
2761 {
2762     int r0 = 0 + (r < 1 ? 1 : 0);
2763     int r1 = 1 + (r < 2 ? 1 : 0);
2764     int c0 = 0 + (c < 1 ? 1 : 0);
2765     int c1 = 1 + (c < 2 ? 1 : 0);
2766
2767     return x[r0][c0] * x[r1][c1] - x[r1][c0] * x[r0][c1];
2768 }
2769
2770 template <class T>
2771 IMATH_HOSTDEVICE constexpr inline T
2772 Matrix33<T>::fastMinor (const int r0, const int r1, const int c0, const int c1) const IMATH_NOEXCEPT
2773 {
2774     return x[r0][c0] * x[r1][c1] - x[r0][c1] * x[r1][c0];
2775 }
2776
2777 template <class T>
2778 IMATH_HOSTDEVICE constexpr inline T
2779 Matrix33<T>::determinant() const IMATH_NOEXCEPT
2780 {
2781     return x[0][0] * (x[1][1] * x[2][2] - x[1][2] * x[2][1]) +
2782            x[0][1] * (x[1][2] * x[2][0] - x[1][0] * x[2][2]) +
2783            x[0][2] * (x[1][0] * x[2][1] - x[1][1] * x[2][0]);
2784 }
2785
2786 template <class T>
2787 IMATH_HOSTDEVICE constexpr inline T
2788 Matrix33<T>::trace () const IMATH_NOEXCEPT
2789 {
2790     return x[0][0] + x[1][1] + x[2][2];
2791 }
2792
2793 template <class T>
2794 template <class S>
2795 IMATH_HOSTDEVICE inline const Matrix33<T>&
2796 Matrix33<T>::setRotation (S r) IMATH_NOEXCEPT
2797 {
2798     S cos_r, sin_r;
2799
2800     cos_r = cos ((T) r);
2801     sin_r = sin ((T) r);
2802
2803     x[0][0] = cos_r;
2804     x[0][1] = sin_r;
2805     x[0][2] = 0;
2806
2807     x[1][0] = -sin_r;
2808     x[1][1] = cos_r;
2809     x[1][2] = 0;
2810
2811     x[2][0] = 0;
2812     x[2][1] = 0;
2813     x[2][2] = 1;
2814
2815     return *this;
2816 }
2817
2818 template <class T>
2819 template <class S>
2820 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix33<T>&
2821 Matrix33<T>::rotate (S r) IMATH_NOEXCEPT
2822 {
2823     *this *= Matrix33<T>().setRotation (r);
2824     return *this;
2825 }
2826
2827 template <class T>
2828 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix33<T>&
2829 Matrix33<T>::setScale (T s) IMATH_NOEXCEPT
2830 {
2831     //
2832     // Set the matrix to a 2D homogeneous transform scale:
2833     //  | s 0 0 |
2834     //  | 0 s 0 |
2835     //  | 0 0 1 |
2836     //
2837
2838     x[0][0] = s;
2839     x[0][1] = 0;
2840     x[0][2] = 0;
2841     x[1][0] = 0;
2842     x[1][1] = s;
2843     x[1][2] = 0;
2844     x[2][0] = 0;
2845     x[2][1] = 0;
2846     x[2][2] = 1;
2847     return *this;
2848 }
2849
2850 template <class T>
2851 template <class S>
2852 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix33<T>&
2853 Matrix33<T>::setScale (const Vec2<S>& s) IMATH_NOEXCEPT
2854 {
2855     //
2856     // Set the matrix to a 2D homogeneous transform scale:
2857     //  | s.x  0   0 |
2858     //  |  0  s.y  0 |
2859     //  |  0   0   1 |
2860     //
2861
2862     x[0][0] = s.x;
2863     x[0][1] = 0;
2864     x[0][2] = 0;
2865     x[1][0] = 0;
2866     x[1][1] = s.y;
2867     x[1][2] = 0;
2868     x[2][0] = 0;
2869     x[2][1] = 0;
2870     x[2][2] = 1;
2871     return *this;
2872 }
2873
2874 template <class T>
2875 template <class S>
2876 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix33<T>&
2877 Matrix33<T>::scale (const Vec2<S>& s) IMATH_NOEXCEPT
2878 {
2879     x[0][0] *= s.x;
2880     x[0][1] *= s.x;
2881     x[0][2] *= s.x;
2882
2883     x[1][0] *= s.y;
2884     x[1][1] *= s.y;
2885     x[1][2] *= s.y;
2886
2887     return *this;
2888 }
2889
2890 template <class T>
2891 template <class S>
2892 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix33<T>&
2893 Matrix33<T>::setTranslation (const Vec2<S>& t) IMATH_NOEXCEPT
2894 {
2895     x[0][0] = 1;
2896     x[0][1] = 0;
2897     x[0][2] = 0;
2898
2899     x[1][0] = 0;
2900     x[1][1] = 1;
2901     x[1][2] = 0;
2902
2903     x[2][0] = t.x;
2904     x[2][1] = t.y;
2905     x[2][2] = 1;
2906
2907     return *this;
2908 }
2909
2910 template <class T>
2911 IMATH_HOSTDEVICE constexpr inline Vec2<T>
2912 Matrix33<T>::translation() const IMATH_NOEXCEPT
2913 {
2914     return Vec2<T> (x[2][0], x[2][1]);
2915 }
2916
2917 template <class T>
2918 template <class S>
2919 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix33<T>&
2920 Matrix33<T>::translate (const Vec2<S>& t) IMATH_NOEXCEPT
2921 {
2922     x[2][0] += t.x * x[0][0] + t.y * x[1][0];
2923     x[2][1] += t.x * x[0][1] + t.y * x[1][1];
2924     x[2][2] += t.x * x[0][2] + t.y * x[1][2];
2925
2926     return *this;
2927 }
2928
2929 template <class T>
2930 template <class S>
2931 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix33<T>&
2932 Matrix33<T>::setShear (const S& xy) IMATH_NOEXCEPT
2933 {
2934     x[0][0] = 1;
2935     x[0][1] = 0;
2936     x[0][2] = 0;
2937
2938     x[1][0] = xy;
2939     x[1][1] = 1;
2940     x[1][2] = 0;
2941
2942     x[2][0] = 0;
2943     x[2][1] = 0;
2944     x[2][2] = 1;
2945
2946     return *this;
2947 }
2948
2949 template <class T>
2950 template <class S>
2951 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix33<T>&
2952 Matrix33<T>::setShear (const Vec2<S>& h) IMATH_NOEXCEPT
2953 {
2954     x[0][0] = 1;
2955     x[0][1] = h.y;
2956     x[0][2] = 0;
2957
2958     x[1][0] = h.x;
2959     x[1][1] = 1;
2960     x[1][2] = 0;
2961
2962     x[2][0] = 0;
2963     x[2][1] = 0;
2964     x[2][2] = 1;
2965
2966     return *this;
2967 }
2968
2969 template <class T>
2970 template <class S>
2971 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix33<T>&
2972 Matrix33<T>::shear (const S& xy) IMATH_NOEXCEPT
2973 {
2974     //
2975     // In this case, we don't need a temp. copy of the matrix
2976     // because we never use a value on the RHS after we've
2977     // changed it on the LHS.
2978     //
2979
2980     x[1][0] += xy * x[0][0];
2981     x[1][1] += xy * x[0][1];
2982     x[1][2] += xy * x[0][2];
2983
2984     return *this;
2985 }
2986
2987 template <class T>
2988 template <class S>
2989 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix33<T>&
2990 Matrix33<T>::shear (const Vec2<S>& h) IMATH_NOEXCEPT
2991 {
2992     Matrix33<T> P (*this);
2993
2994     x[0][0] = P.x[0][0] + h.y * P.x[1][0];
2995     x[0][1] = P.x[0][1] + h.y * P.x[1][1];
2996     x[0][2] = P.x[0][2] + h.y * P.x[1][2];
2997
2998     x[1][0] = P.x[1][0] + h.x * P.x[0][0];
2999     x[1][1] = P.x[1][1] + h.x * P.x[0][1];
3000     x[1][2] = P.x[1][2] + h.x * P.x[0][2];
3001
3002     return *this;
3003 }
3004
3005 //---------------------------
3006 // Implementation of Matrix44
3007 //---------------------------
3008
3009 template <class T>
3010 IMATH_HOSTDEVICE inline T*
3011 Matrix44<T>::operator[] (int i) IMATH_NOEXCEPT
3012 {
3013     return x[i];
3014 }
3015
3016 template <class T>
3017 IMATH_HOSTDEVICE inline const T*
3018 Matrix44<T>::operator[] (int i) const IMATH_NOEXCEPT
3019 {
3020     return x[i];
3021 }
3022
3023 template <class T> IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix44<T>::Matrix44() IMATH_NOEXCEPT
3024 {
3025     x[0][0] = 1;
3026     x[0][1] = 0;
3027     x[0][2] = 0;
3028     x[0][3] = 0;
3029     x[1][0] = 0;
3030     x[1][1] = 1;
3031     x[1][2] = 0;
3032     x[1][3] = 0;
3033     x[2][0] = 0;
3034     x[2][1] = 0;
3035     x[2][2] = 1;
3036     x[2][3] = 0;
3037     x[3][0] = 0;
3038     x[3][1] = 0;
3039     x[3][2] = 0;
3040     x[3][3] = 1;
3041 }
3042
3043 template <class T> IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix44<T>::Matrix44 (T a) IMATH_NOEXCEPT
3044 {
3045     x[0][0] = a;
3046     x[0][1] = a;
3047     x[0][2] = a;
3048     x[0][3] = a;
3049     x[1][0] = a;
3050     x[1][1] = a;
3051     x[1][2] = a;
3052     x[1][3] = a;
3053     x[2][0] = a;
3054     x[2][1] = a;
3055     x[2][2] = a;
3056     x[2][3] = a;
3057     x[3][0] = a;
3058     x[3][1] = a;
3059     x[3][2] = a;
3060     x[3][3] = a;
3061 }
3062
3063 template <class T>
3064 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix44<T>::Matrix44 (
3065     const T a[4][4]) IMATH_NOEXCEPT
3066 {
3067     x[0][0] = a[0][0];
3068     x[0][1] = a[0][1];
3069     x[0][2] = a[0][2];
3070     x[0][3] = a[0][3];
3071     x[1][0] = a[1][0];
3072     x[1][1] = a[1][1];
3073     x[1][2] = a[1][2];
3074     x[1][3] = a[1][3];
3075     x[2][0] = a[2][0];
3076     x[2][1] = a[2][1];
3077     x[2][2] = a[2][2];
3078     x[2][3] = a[2][3];
3079     x[3][0] = a[3][0];
3080     x[3][1] = a[3][1];
3081     x[3][2] = a[3][2];
3082     x[3][3] = a[3][3];
3083 }
3084
3085 template <class T>
3086 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix44<
3087     T>::Matrix44 (T a, T b, T c, T d, T e, T f, T g, T h, T i, T j, T k, T l, T m, T n, T o, T p) IMATH_NOEXCEPT
3088 {
3089     x[0][0] = a;
3090     x[0][1] = b;
3091     x[0][2] = c;
3092     x[0][3] = d;
3093     x[1][0] = e;
3094     x[1][1] = f;
3095     x[1][2] = g;
3096     x[1][3] = h;
3097     x[2][0] = i;
3098     x[2][1] = j;
3099     x[2][2] = k;
3100     x[2][3] = l;
3101     x[3][0] = m;
3102     x[3][1] = n;
3103     x[3][2] = o;
3104     x[3][3] = p;
3105 }
3106
3107 template <class T> IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix44<T>::Matrix44 (Matrix33<T> r, Vec3<T> t) IMATH_NOEXCEPT
3108 {
3109     x[0][0] = r.x[0][0];
3110     x[0][1] = r.x[0][1];
3111     x[0][2] = r.x[0][2];
3112     x[0][3] = 0;
3113     x[1][0] = r.x[1][0];
3114     x[1][1] = r.x[1][1];
3115     x[1][2] = r.x[1][2];
3116     x[1][3] = 0;
3117     x[2][0] = r.x[2][0];
3118     x[2][1] = r.x[2][1];
3119     x[2][2] = r.x[2][2];
3120     x[2][3] = 0;
3121     x[3][0] = t.x;
3122     x[3][1] = t.y;
3123     x[3][2] = t.z;
3124     x[3][3] = 1;
3125 }
3126
3127 template <class T> IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix44<T>::Matrix44 (const Matrix44& v) IMATH_NOEXCEPT
3128 {
3129     x[0][0] = v.x[0][0];
3130     x[0][1] = v.x[0][1];
3131     x[0][2] = v.x[0][2];
3132     x[0][3] = v.x[0][3];
3133     x[1][0] = v.x[1][0];
3134     x[1][1] = v.x[1][1];
3135     x[1][2] = v.x[1][2];
3136     x[1][3] = v.x[1][3];
3137     x[2][0] = v.x[2][0];
3138     x[2][1] = v.x[2][1];
3139     x[2][2] = v.x[2][2];
3140     x[2][3] = v.x[2][3];
3141     x[3][0] = v.x[3][0];
3142     x[3][1] = v.x[3][1];
3143     x[3][2] = v.x[3][2];
3144     x[3][3] = v.x[3][3];
3145 }
3146
3147 template <class T>
3148 template <class S>
3149 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix44<T>::Matrix44 (const Matrix44<S>& v) IMATH_NOEXCEPT
3150 {
3151     x[0][0] = T (v.x[0][0]);
3152     x[0][1] = T (v.x[0][1]);
3153     x[0][2] = T (v.x[0][2]);
3154     x[0][3] = T (v.x[0][3]);
3155     x[1][0] = T (v.x[1][0]);
3156     x[1][1] = T (v.x[1][1]);
3157     x[1][2] = T (v.x[1][2]);
3158     x[1][3] = T (v.x[1][3]);
3159     x[2][0] = T (v.x[2][0]);
3160     x[2][1] = T (v.x[2][1]);
3161     x[2][2] = T (v.x[2][2]);
3162     x[2][3] = T (v.x[2][3]);
3163     x[3][0] = T (v.x[3][0]);
3164     x[3][1] = T (v.x[3][1]);
3165     x[3][2] = T (v.x[3][2]);
3166     x[3][3] = T (v.x[3][3]);
3167 }
3168
3169 template <class T>
3170 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix44<T>&
3171 Matrix44<T>::operator= (const Matrix44& v) IMATH_NOEXCEPT
3172 {
3173     x[0][0] = v.x[0][0];
3174     x[0][1] = v.x[0][1];
3175     x[0][2] = v.x[0][2];
3176     x[0][3] = v.x[0][3];
3177     x[1][0] = v.x[1][0];
3178     x[1][1] = v.x[1][1];
3179     x[1][2] = v.x[1][2];
3180     x[1][3] = v.x[1][3];
3181     x[2][0] = v.x[2][0];
3182     x[2][1] = v.x[2][1];
3183     x[2][2] = v.x[2][2];
3184     x[2][3] = v.x[2][3];
3185     x[3][0] = v.x[3][0];
3186     x[3][1] = v.x[3][1];
3187     x[3][2] = v.x[3][2];
3188     x[3][3] = v.x[3][3];
3189     return *this;
3190 }
3191
3192 template <class T>
3193 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix44<T>&
3194 Matrix44<T>::operator= (T a) IMATH_NOEXCEPT
3195 {
3196     x[0][0] = a;
3197     x[0][1] = a;
3198     x[0][2] = a;
3199     x[0][3] = a;
3200     x[1][0] = a;
3201     x[1][1] = a;
3202     x[1][2] = a;
3203     x[1][3] = a;
3204     x[2][0] = a;
3205     x[2][1] = a;
3206     x[2][2] = a;
3207     x[2][3] = a;
3208     x[3][0] = a;
3209     x[3][1] = a;
3210     x[3][2] = a;
3211     x[3][3] = a;
3212     return *this;
3213 }
3214
3215 template <class T>
3216 IMATH_HOSTDEVICE inline T*
3217 Matrix44<T>::getValue () IMATH_NOEXCEPT
3218 {
3219     return (T*) &x[0][0];
3220 }
3221
3222 template <class T>
3223 IMATH_HOSTDEVICE inline const T*
3224 Matrix44<T>::getValue() const IMATH_NOEXCEPT
3225 {
3226     return (const T*) &x[0][0];
3227 }
3228
3229 template <class T>
3230 template <class S>
3231 IMATH_HOSTDEVICE inline void
3232 Matrix44<T>::getValue (Matrix44<S>& v) const IMATH_NOEXCEPT
3233 {
3234     v.x[0][0] = x[0][0];
3235     v.x[0][1] = x[0][1];
3236     v.x[0][2] = x[0][2];
3237     v.x[0][3] = x[0][3];
3238     v.x[1][0] = x[1][0];
3239     v.x[1][1] = x[1][1];
3240     v.x[1][2] = x[1][2];
3241     v.x[1][3] = x[1][3];
3242     v.x[2][0] = x[2][0];
3243     v.x[2][1] = x[2][1];
3244     v.x[2][2] = x[2][2];
3245     v.x[2][3] = x[2][3];
3246     v.x[3][0] = x[3][0];
3247     v.x[3][1] = x[3][1];
3248     v.x[3][2] = x[3][2];
3249     v.x[3][3] = x[3][3];
3250 }
3251
3252 template <class T>
3253 template <class S>
3254 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix44<T>&
3255 Matrix44<T>::setValue (const Matrix44<S>& v) IMATH_NOEXCEPT
3256 {
3257     x[0][0] = T(v.x[0][0]);
3258     x[0][1] = T(v.x[0][1]);
3259     x[0][2] = T(v.x[0][2]);
3260     x[0][3] = T(v.x[0][3]);
3261     x[1][0] = T(v.x[1][0]);
3262     x[1][1] = T(v.x[1][1]);
3263     x[1][2] = T(v.x[1][2]);
3264     x[1][3] = T(v.x[1][3]);
3265     x[2][0] = T(v.x[2][0]);
3266     x[2][1] = T(v.x[2][1]);
3267     x[2][2] = T(v.x[2][2]);
3268     x[2][3] = T(v.x[2][3]);
3269     x[3][0] = T(v.x[3][0]);
3270     x[3][1] = T(v.x[3][1]);
3271     x[3][2] = T(v.x[3][2]);
3272     x[3][3] = T(v.x[3][3]);
3273     return *this;
3274 }
3275
3276 template <class T>
3277 template <class S>
3278 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix44<T>&
3279 Matrix44<T>::setTheMatrix (const Matrix44<S>& v) IMATH_NOEXCEPT
3280 {
3281     x[0][0] = v.x[0][0];
3282     x[0][1] = v.x[0][1];
3283     x[0][2] = v.x[0][2];
3284     x[0][3] = v.x[0][3];
3285     x[1][0] = v.x[1][0];
3286     x[1][1] = v.x[1][1];
3287     x[1][2] = v.x[1][2];
3288     x[1][3] = v.x[1][3];
3289     x[2][0] = v.x[2][0];
3290     x[2][1] = v.x[2][1];
3291     x[2][2] = v.x[2][2];
3292     x[2][3] = v.x[2][3];
3293     x[3][0] = v.x[3][0];
3294     x[3][1] = v.x[3][1];
3295     x[3][2] = v.x[3][2];
3296     x[3][3] = v.x[3][3];
3297     return *this;
3298 }
3299
3300 template <class T>
3301 IMATH_HOSTDEVICE inline void
3302 Matrix44<T>::makeIdentity() IMATH_NOEXCEPT
3303 {
3304     x[0][0] = 1;
3305     x[0][1] = 0;
3306     x[0][2] = 0;
3307     x[0][3] = 0;
3308     x[1][0] = 0;
3309     x[1][1] = 1;
3310     x[1][2] = 0;
3311     x[1][3] = 0;
3312     x[2][0] = 0;
3313     x[2][1] = 0;
3314     x[2][2] = 1;
3315     x[2][3] = 0;
3316     x[3][0] = 0;
3317     x[3][1] = 0;
3318     x[3][2] = 0;
3319     x[3][3] = 1;
3320 }
3321
3322 template <class T>
3323 IMATH_HOSTDEVICE constexpr inline bool
3324 Matrix44<T>::operator== (const Matrix44& v) const IMATH_NOEXCEPT
3325 {
3326     return x[0][0] == v.x[0][0] && x[0][1] == v.x[0][1] && x[0][2] == v.x[0][2] &&
3327            x[0][3] == v.x[0][3] && x[1][0] == v.x[1][0] && x[1][1] == v.x[1][1] &&
3328            x[1][2] == v.x[1][2] && x[1][3] == v.x[1][3] && x[2][0] == v.x[2][0] &&
3329            x[2][1] == v.x[2][1] && x[2][2] == v.x[2][2] && x[2][3] == v.x[2][3] &&
3330            x[3][0] == v.x[3][0] && x[3][1] == v.x[3][1] && x[3][2] == v.x[3][2] &&
3331            x[3][3] == v.x[3][3];
3332 }
3333
3334 template <class T>
3335 IMATH_HOSTDEVICE constexpr inline bool
3336 Matrix44<T>::operator!= (const Matrix44& v) const IMATH_NOEXCEPT
3337 {
3338     return x[0][0] != v.x[0][0] || x[0][1] != v.x[0][1] || x[0][2] != v.x[0][2] ||
3339            x[0][3] != v.x[0][3] || x[1][0] != v.x[1][0] || x[1][1] != v.x[1][1] ||
3340            x[1][2] != v.x[1][2] || x[1][3] != v.x[1][3] || x[2][0] != v.x[2][0] ||
3341            x[2][1] != v.x[2][1] || x[2][2] != v.x[2][2] || x[2][3] != v.x[2][3] ||
3342            x[3][0] != v.x[3][0] || x[3][1] != v.x[3][1] || x[3][2] != v.x[3][2] ||
3343            x[3][3] != v.x[3][3];
3344 }
3345
3346 template <class T>
3347 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline bool
3348 Matrix44<T>::equalWithAbsError (const Matrix44<T>& m, T e) const IMATH_NOEXCEPT
3349 {
3350     for (int i = 0; i < 4; i++)
3351         for (int j = 0; j < 4; j++)
3352             if (!IMATH_INTERNAL_NAMESPACE::equalWithAbsError ((*this).x[i][j], m.x[i][j], e))
3353                 return false;
3354
3355     return true;
3356 }
3357
3358 template <class T>
3359 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline bool
3360 Matrix44<T>::equalWithRelError (const Matrix44<T>& m, T e) const IMATH_NOEXCEPT
3361 {
3362     for (int i = 0; i < 4; i++)
3363         for (int j = 0; j < 4; j++)
3364             if (!IMATH_INTERNAL_NAMESPACE::equalWithRelError ((*this).x[i][j], m.x[i][j], e))
3365                 return false;
3366
3367     return true;
3368 }
3369
3370 template <class T>
3371 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix44<T>&
3372 Matrix44<T>::operator+= (const Matrix44<T>& v) IMATH_NOEXCEPT
3373 {
3374     x[0][0] += v.x[0][0];
3375     x[0][1] += v.x[0][1];
3376     x[0][2] += v.x[0][2];
3377     x[0][3] += v.x[0][3];
3378     x[1][0] += v.x[1][0];
3379     x[1][1] += v.x[1][1];
3380     x[1][2] += v.x[1][2];
3381     x[1][3] += v.x[1][3];
3382     x[2][0] += v.x[2][0];
3383     x[2][1] += v.x[2][1];
3384     x[2][2] += v.x[2][2];
3385     x[2][3] += v.x[2][3];
3386     x[3][0] += v.x[3][0];
3387     x[3][1] += v.x[3][1];
3388     x[3][2] += v.x[3][2];
3389     x[3][3] += v.x[3][3];
3390
3391     return *this;
3392 }
3393
3394 template <class T>
3395 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix44<T>&
3396 Matrix44<T>::operator+= (T a) IMATH_NOEXCEPT
3397 {
3398     x[0][0] += a;
3399     x[0][1] += a;
3400     x[0][2] += a;
3401     x[0][3] += a;
3402     x[1][0] += a;
3403     x[1][1] += a;
3404     x[1][2] += a;
3405     x[1][3] += a;
3406     x[2][0] += a;
3407     x[2][1] += a;
3408     x[2][2] += a;
3409     x[2][3] += a;
3410     x[3][0] += a;
3411     x[3][1] += a;
3412     x[3][2] += a;
3413     x[3][3] += a;
3414
3415     return *this;
3416 }
3417
3418 template <class T>
3419 IMATH_HOSTDEVICE constexpr inline Matrix44<T>
3420 Matrix44<T>::operator+ (const Matrix44<T>& v) const IMATH_NOEXCEPT
3421 {
3422     return Matrix44 (x[0][0] + v.x[0][0],
3423                      x[0][1] + v.x[0][1],
3424                      x[0][2] + v.x[0][2],
3425                      x[0][3] + v.x[0][3],
3426                      x[1][0] + v.x[1][0],
3427                      x[1][1] + v.x[1][1],
3428                      x[1][2] + v.x[1][2],
3429                      x[1][3] + v.x[1][3],
3430                      x[2][0] + v.x[2][0],
3431                      x[2][1] + v.x[2][1],
3432                      x[2][2] + v.x[2][2],
3433                      x[2][3] + v.x[2][3],
3434                      x[3][0] + v.x[3][0],
3435                      x[3][1] + v.x[3][1],
3436                      x[3][2] + v.x[3][2],
3437                      x[3][3] + v.x[3][3]);
3438 }
3439
3440 template <class T>
3441 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix44<T>&
3442 Matrix44<T>::operator-= (const Matrix44<T>& v) IMATH_NOEXCEPT
3443 {
3444     x[0][0] -= v.x[0][0];
3445     x[0][1] -= v.x[0][1];
3446     x[0][2] -= v.x[0][2];
3447     x[0][3] -= v.x[0][3];
3448     x[1][0] -= v.x[1][0];
3449     x[1][1] -= v.x[1][1];
3450     x[1][2] -= v.x[1][2];
3451     x[1][3] -= v.x[1][3];
3452     x[2][0] -= v.x[2][0];
3453     x[2][1] -= v.x[2][1];
3454     x[2][2] -= v.x[2][2];
3455     x[2][3] -= v.x[2][3];
3456     x[3][0] -= v.x[3][0];
3457     x[3][1] -= v.x[3][1];
3458     x[3][2] -= v.x[3][2];
3459     x[3][3] -= v.x[3][3];
3460
3461     return *this;
3462 }
3463
3464 template <class T>
3465 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix44<T>&
3466 Matrix44<T>::operator-= (T a) IMATH_NOEXCEPT
3467 {
3468     x[0][0] -= a;
3469     x[0][1] -= a;
3470     x[0][2] -= a;
3471     x[0][3] -= a;
3472     x[1][0] -= a;
3473     x[1][1] -= a;
3474     x[1][2] -= a;
3475     x[1][3] -= a;
3476     x[2][0] -= a;
3477     x[2][1] -= a;
3478     x[2][2] -= a;
3479     x[2][3] -= a;
3480     x[3][0] -= a;
3481     x[3][1] -= a;
3482     x[3][2] -= a;
3483     x[3][3] -= a;
3484
3485     return *this;
3486 }
3487
3488 template <class T>
3489 IMATH_HOSTDEVICE constexpr inline Matrix44<T>
3490 Matrix44<T>::operator- (const Matrix44<T>& v) const IMATH_NOEXCEPT
3491 {
3492     return Matrix44 (x[0][0] - v.x[0][0],
3493                      x[0][1] - v.x[0][1],
3494                      x[0][2] - v.x[0][2],
3495                      x[0][3] - v.x[0][3],
3496                      x[1][0] - v.x[1][0],
3497                      x[1][1] - v.x[1][1],
3498                      x[1][2] - v.x[1][2],
3499                      x[1][3] - v.x[1][3],
3500                      x[2][0] - v.x[2][0],
3501                      x[2][1] - v.x[2][1],
3502                      x[2][2] - v.x[2][2],
3503                      x[2][3] - v.x[2][3],
3504                      x[3][0] - v.x[3][0],
3505                      x[3][1] - v.x[3][1],
3506                      x[3][2] - v.x[3][2],
3507                      x[3][3] - v.x[3][3]);
3508 }
3509
3510 template <class T>
3511 IMATH_HOSTDEVICE constexpr inline Matrix44<T>
3512 Matrix44<T>::operator-() const IMATH_NOEXCEPT
3513 {
3514     return Matrix44 (-x[0][0],
3515                      -x[0][1],
3516                      -x[0][2],
3517                      -x[0][3],
3518                      -x[1][0],
3519                      -x[1][1],
3520                      -x[1][2],
3521                      -x[1][3],
3522                      -x[2][0],
3523                      -x[2][1],
3524                      -x[2][2],
3525                      -x[2][3],
3526                      -x[3][0],
3527                      -x[3][1],
3528                      -x[3][2],
3529                      -x[3][3]);
3530 }
3531
3532 template <class T>
3533 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix44<T>&
3534 Matrix44<T>::negate() IMATH_NOEXCEPT
3535 {
3536     x[0][0] = -x[0][0];
3537     x[0][1] = -x[0][1];
3538     x[0][2] = -x[0][2];
3539     x[0][3] = -x[0][3];
3540     x[1][0] = -x[1][0];
3541     x[1][1] = -x[1][1];
3542     x[1][2] = -x[1][2];
3543     x[1][3] = -x[1][3];
3544     x[2][0] = -x[2][0];
3545     x[2][1] = -x[2][1];
3546     x[2][2] = -x[2][2];
3547     x[2][3] = -x[2][3];
3548     x[3][0] = -x[3][0];
3549     x[3][1] = -x[3][1];
3550     x[3][2] = -x[3][2];
3551     x[3][3] = -x[3][3];
3552
3553     return *this;
3554 }
3555
3556 template <class T>
3557 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix44<T>&
3558 Matrix44<T>::operator*= (T a) IMATH_NOEXCEPT
3559 {
3560     x[0][0] *= a;
3561     x[0][1] *= a;
3562     x[0][2] *= a;
3563     x[0][3] *= a;
3564     x[1][0] *= a;
3565     x[1][1] *= a;
3566     x[1][2] *= a;
3567     x[1][3] *= a;
3568     x[2][0] *= a;
3569     x[2][1] *= a;
3570     x[2][2] *= a;
3571     x[2][3] *= a;
3572     x[3][0] *= a;
3573     x[3][1] *= a;
3574     x[3][2] *= a;
3575     x[3][3] *= a;
3576
3577     return *this;
3578 }
3579
3580 template <class T>
3581 IMATH_HOSTDEVICE constexpr inline Matrix44<T>
3582 Matrix44<T>::operator* (T a) const IMATH_NOEXCEPT
3583 {
3584     return Matrix44 (x[0][0] * a,
3585                      x[0][1] * a,
3586                      x[0][2] * a,
3587                      x[0][3] * a,
3588                      x[1][0] * a,
3589                      x[1][1] * a,
3590                      x[1][2] * a,
3591                      x[1][3] * a,
3592                      x[2][0] * a,
3593                      x[2][1] * a,
3594                      x[2][2] * a,
3595                      x[2][3] * a,
3596                      x[3][0] * a,
3597                      x[3][1] * a,
3598                      x[3][2] * a,
3599                      x[3][3] * a);
3600 }
3601
3602 /// Matrix-scalar multiplication
3603 template <class T>
3604 IMATH_HOSTDEVICE inline Matrix44<T>
3605 operator* (T a, const Matrix44<T>& v) IMATH_NOEXCEPT
3606 {
3607     return v * a;
3608 }
3609
3610
3611 template <class T>
3612 IMATH_HOSTDEVICE inline IMATH_CONSTEXPR14 Matrix44<T>
3613 Matrix44<T>::multiply (const Matrix44 &a, const Matrix44 &b) IMATH_NOEXCEPT
3614 {
3615     const auto a00 = a.x[0][0];
3616     const auto a01 = a.x[0][1];
3617     const auto a02 = a.x[0][2];
3618     const auto a03 = a.x[0][3];
3619
3620     const auto c00  = a00 * b.x[0][0] + a01 * b.x[1][0] + a02 * b.x[2][0] + a03 * b.x[3][0];
3621     const auto c01  = a00 * b.x[0][1] + a01 * b.x[1][1] + a02 * b.x[2][1] + a03 * b.x[3][1];
3622     const auto c02  = a00 * b.x[0][2] + a01 * b.x[1][2] + a02 * b.x[2][2] + a03 * b.x[3][2];
3623     const auto c03  = a00 * b.x[0][3] + a01 * b.x[1][3] + a02 * b.x[2][3] + a03 * b.x[3][3];
3624
3625     const auto a10 = a.x[1][0];
3626     const auto a11 = a.x[1][1];
3627     const auto a12 = a.x[1][2];
3628     const auto a13 = a.x[1][3];
3629
3630     const auto c10  = a10 * b.x[0][0] + a11 * b.x[1][0] + a12 * b.x[2][0] + a13 * b.x[3][0];
3631     const auto c11  = a10 * b.x[0][1] + a11 * b.x[1][1] + a12 * b.x[2][1] + a13 * b.x[3][1];
3632     const auto c12  = a10 * b.x[0][2] + a11 * b.x[1][2] + a12 * b.x[2][2] + a13 * b.x[3][2];
3633     const auto c13  = a10 * b.x[0][3] + a11 * b.x[1][3] + a12 * b.x[2][3] + a13 * b.x[3][3];
3634
3635     const auto a20 = a.x[2][0];
3636     const auto a21 = a.x[2][1];
3637     const auto a22 = a.x[2][2];
3638     const auto a23 = a.x[2][3];
3639
3640     const auto c20 = a20 * b.x[0][0] + a21 * b.x[1][0] + a22 * b.x[2][0] + a23 * b.x[3][0];
3641     const auto c21 = a20 * b.x[0][1] + a21 * b.x[1][1] + a22 * b.x[2][1] + a23 * b.x[3][1];
3642     const auto c22 = a20 * b.x[0][2] + a21 * b.x[1][2] + a22 * b.x[2][2] + a23 * b.x[3][2];
3643     const auto c23 = a20 * b.x[0][3] + a21 * b.x[1][3] + a22 * b.x[2][3] + a23 * b.x[3][3];
3644
3645     const auto a30 = a.x[3][0];
3646     const auto a31 = a.x[3][1];
3647     const auto a32 = a.x[3][2];
3648     const auto a33 = a.x[3][3];
3649
3650     const auto c30 = a30 * b.x[0][0] + a31 * b.x[1][0] + a32 * b.x[2][0] + a33 * b.x[3][0];
3651     const auto c31 = a30 * b.x[0][1] + a31 * b.x[1][1] + a32 * b.x[2][1] + a33 * b.x[3][1];
3652     const auto c32 = a30 * b.x[0][2] + a31 * b.x[1][2] + a32 * b.x[2][2] + a33 * b.x[3][2];
3653     const auto c33 = a30 * b.x[0][3] + a31 * b.x[1][3] + a32 * b.x[2][3] + a33 * b.x[3][3];
3654     return Matrix44(c00, c01, c02, c03,
3655                     c10, c11, c12, c13,
3656                     c20, c21, c22, c23,
3657                     c30, c31, c32, c33);
3658 }
3659
3660
3661 template <class T>
3662 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix44<T>&
3663 Matrix44<T>::operator*= (const Matrix44<T>& v) IMATH_NOEXCEPT
3664 {
3665     *this = multiply(*this, v);
3666     return *this;
3667 }
3668
3669 template <class T>
3670 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix44<T>
3671 Matrix44<T>::operator* (const Matrix44<T>& v) const IMATH_NOEXCEPT
3672 {
3673     return multiply(*this, v);
3674 }
3675
3676 template <class T>
3677 IMATH_HOSTDEVICE inline void
3678 Matrix44<T>::multiply (const Matrix44<T>& a, const Matrix44<T>& b, Matrix44<T>& c) IMATH_NOEXCEPT
3679 {
3680     c = multiply(a, b);
3681 }
3682
3683 template <class T>
3684 template <class S>
3685 IMATH_HOSTDEVICE inline void
3686 Matrix44<T>::multVecMatrix (const Vec3<S>& src, Vec3<S>& dst) const IMATH_NOEXCEPT
3687 {
3688     S a, b, c, w;
3689
3690     a = src.x * x[0][0] + src.y * x[1][0] + src.z * x[2][0] + x[3][0];
3691     b = src.x * x[0][1] + src.y * x[1][1] + src.z * x[2][1] + x[3][1];
3692     c = src.x * x[0][2] + src.y * x[1][2] + src.z * x[2][2] + x[3][2];
3693     w = src.x * x[0][3] + src.y * x[1][3] + src.z * x[2][3] + x[3][3];
3694
3695     dst.x = a / w;
3696     dst.y = b / w;
3697     dst.z = c / w;
3698 }
3699
3700 template <class T>
3701 template <class S>
3702 IMATH_HOSTDEVICE inline void
3703 Matrix44<T>::multDirMatrix (const Vec3<S>& src, Vec3<S>& dst) const IMATH_NOEXCEPT
3704 {
3705     S a, b, c;
3706
3707     a = src.x * x[0][0] + src.y * x[1][0] + src.z * x[2][0];
3708     b = src.x * x[0][1] + src.y * x[1][1] + src.z * x[2][1];
3709     c = src.x * x[0][2] + src.y * x[1][2] + src.z * x[2][2];
3710
3711     dst.x = a;
3712     dst.y = b;
3713     dst.z = c;
3714 }
3715
3716 template <class T>
3717 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix44<T>&
3718 Matrix44<T>::operator/= (T a) IMATH_NOEXCEPT
3719 {
3720     x[0][0] /= a;
3721     x[0][1] /= a;
3722     x[0][2] /= a;
3723     x[0][3] /= a;
3724     x[1][0] /= a;
3725     x[1][1] /= a;
3726     x[1][2] /= a;
3727     x[1][3] /= a;
3728     x[2][0] /= a;
3729     x[2][1] /= a;
3730     x[2][2] /= a;
3731     x[2][3] /= a;
3732     x[3][0] /= a;
3733     x[3][1] /= a;
3734     x[3][2] /= a;
3735     x[3][3] /= a;
3736
3737     return *this;
3738 }
3739
3740 template <class T>
3741 IMATH_HOSTDEVICE constexpr inline Matrix44<T>
3742 Matrix44<T>::operator/ (T a) const IMATH_NOEXCEPT
3743 {
3744     return Matrix44 (x[0][0] / a,
3745                      x[0][1] / a,
3746                      x[0][2] / a,
3747                      x[0][3] / a,
3748                      x[1][0] / a,
3749                      x[1][1] / a,
3750                      x[1][2] / a,
3751                      x[1][3] / a,
3752                      x[2][0] / a,
3753                      x[2][1] / a,
3754                      x[2][2] / a,
3755                      x[2][3] / a,
3756                      x[3][0] / a,
3757                      x[3][1] / a,
3758                      x[3][2] / a,
3759                      x[3][3] / a);
3760 }
3761
3762 template <class T>
3763 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix44<T>&
3764 Matrix44<T>::transpose() IMATH_NOEXCEPT
3765 {
3766     Matrix44 tmp (x[0][0],
3767                   x[1][0],
3768                   x[2][0],
3769                   x[3][0],
3770                   x[0][1],
3771                   x[1][1],
3772                   x[2][1],
3773                   x[3][1],
3774                   x[0][2],
3775                   x[1][2],
3776                   x[2][2],
3777                   x[3][2],
3778                   x[0][3],
3779                   x[1][3],
3780                   x[2][3],
3781                   x[3][3]);
3782     *this = tmp;
3783     return *this;
3784 }
3785
3786 template <class T>
3787 IMATH_HOSTDEVICE constexpr inline Matrix44<T>
3788 Matrix44<T>::transposed() const IMATH_NOEXCEPT
3789 {
3790     return Matrix44 (x[0][0],
3791                      x[1][0],
3792                      x[2][0],
3793                      x[3][0],
3794                      x[0][1],
3795                      x[1][1],
3796                      x[2][1],
3797                      x[3][1],
3798                      x[0][2],
3799                      x[1][2],
3800                      x[2][2],
3801                      x[3][2],
3802                      x[0][3],
3803                      x[1][3],
3804                      x[2][3],
3805                      x[3][3]);
3806 }
3807
3808 template <class T>
3809 IMATH_CONSTEXPR14 inline const Matrix44<T>&
3810 Matrix44<T>::gjInvert (bool singExc)
3811 {
3812     *this = gjInverse (singExc);
3813     return *this;
3814 }
3815
3816 template <class T>
3817 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix44<T>&
3818 Matrix44<T>::gjInvert() IMATH_NOEXCEPT
3819 {
3820     *this = gjInverse();
3821     return *this;
3822 }
3823
3824 template <class T>
3825 inline Matrix44<T>
3826 Matrix44<T>::gjInverse (bool singExc) const
3827 {
3828     int i, j, k;
3829     Matrix44 s;
3830     Matrix44 t (*this);
3831
3832     // Forward elimination
3833
3834     for (i = 0; i < 3; i++)
3835     {
3836         int pivot = i;
3837
3838         T pivotsize = t.x[i][i];
3839
3840         if (pivotsize < 0)
3841             pivotsize = -pivotsize;
3842
3843         for (j = i + 1; j < 4; j++)
3844         {
3845             T tmp = t.x[j][i];
3846
3847             if (tmp < 0)
3848                 tmp = -tmp;
3849
3850             if (tmp > pivotsize)
3851             {
3852                 pivot     = j;
3853                 pivotsize = tmp;
3854             }
3855         }
3856
3857         if (pivotsize == 0)
3858         {
3859             if (singExc)
3860                 throw std::invalid_argument ("Cannot invert singular matrix.");
3861
3862             return Matrix44();
3863         }
3864
3865         if (pivot != i)
3866         {
3867             for (j = 0; j < 4; j++)
3868             {
3869                 T tmp;
3870
3871                 tmp           = t.x[i][j];
3872                 t.x[i][j]     = t.x[pivot][j];
3873                 t.x[pivot][j] = tmp;
3874
3875                 tmp           = s.x[i][j];
3876                 s.x[i][j]     = s.x[pivot][j];
3877                 s.x[pivot][j] = tmp;
3878             }
3879         }
3880
3881         for (j = i + 1; j < 4; j++)
3882         {
3883             T f = t.x[j][i] / t.x[i][i];
3884
3885             for (k = 0; k < 4; k++)
3886             {
3887                 t.x[j][k] -= f * t.x[i][k];
3888                 s.x[j][k] -= f * s.x[i][k];
3889             }
3890         }
3891     }
3892
3893     // Backward substitution
3894
3895     for (i = 3; i >= 0; --i)
3896     {
3897         T f;
3898
3899         if ((f = t.x[i][i]) == 0)
3900         {
3901             if (singExc)
3902                 throw std::invalid_argument ("Cannot invert singular matrix.");
3903
3904             return Matrix44();
3905         }
3906
3907         for (j = 0; j < 4; j++)
3908         {
3909             t.x[i][j] /= f;
3910             s.x[i][j] /= f;
3911         }
3912
3913         for (j = 0; j < i; j++)
3914         {
3915             f = t.x[j][i];
3916
3917             for (k = 0; k < 4; k++)
3918             {
3919                 t.x[j][k] -= f * t.x[i][k];
3920                 s.x[j][k] -= f * s.x[i][k];
3921             }
3922         }
3923     }
3924
3925     return s;
3926 }
3927
3928 template <class T>
3929 IMATH_HOSTDEVICE inline Matrix44<T>
3930 Matrix44<T>::gjInverse() const IMATH_NOEXCEPT
3931 {
3932     int i, j, k;
3933     Matrix44 s;
3934     Matrix44 t (*this);
3935
3936     // Forward elimination
3937
3938     for (i = 0; i < 3; i++)
3939     {
3940         int pivot = i;
3941
3942         T pivotsize = t.x[i][i];
3943
3944         if (pivotsize < 0)
3945             pivotsize = -pivotsize;
3946
3947         for (j = i + 1; j < 4; j++)
3948         {
3949             T tmp = t.x[j][i];
3950
3951             if (tmp < 0)
3952                 tmp = -tmp;
3953
3954             if (tmp > pivotsize)
3955             {
3956                 pivot     = j;
3957                 pivotsize = tmp;
3958             }
3959         }
3960
3961         if (pivotsize == 0)
3962         {
3963             return Matrix44();
3964         }
3965
3966         if (pivot != i)
3967         {
3968             for (j = 0; j < 4; j++)
3969             {
3970                 T tmp;
3971
3972                 tmp           = t.x[i][j];
3973                 t.x[i][j]     = t.x[pivot][j];
3974                 t.x[pivot][j] = tmp;
3975
3976                 tmp           = s.x[i][j];
3977                 s.x[i][j]     = s.x[pivot][j];
3978                 s.x[pivot][j] = tmp;
3979             }
3980         }
3981
3982         for (j = i + 1; j < 4; j++)
3983         {
3984             T f = t.x[j][i] / t.x[i][i];
3985
3986             for (k = 0; k < 4; k++)
3987             {
3988                 t.x[j][k] -= f * t.x[i][k];
3989                 s.x[j][k] -= f * s.x[i][k];
3990             }
3991         }
3992     }
3993
3994     // Backward substitution
3995
3996     for (i = 3; i >= 0; --i)
3997     {
3998         T f;
3999
4000         if ((f = t.x[i][i]) == 0)
4001         {
4002             return Matrix44();
4003         }
4004
4005         for (j = 0; j < 4; j++)
4006         {
4007             t.x[i][j] /= f;
4008             s.x[i][j] /= f;
4009         }
4010
4011         for (j = 0; j < i; j++)
4012         {
4013             f = t.x[j][i];
4014
4015             for (k = 0; k < 4; k++)
4016             {
4017                 t.x[j][k] -= f * t.x[i][k];
4018                 s.x[j][k] -= f * s.x[i][k];
4019             }
4020         }
4021     }
4022
4023     return s;
4024 }
4025
4026 template <class T>
4027 IMATH_CONSTEXPR14 inline const Matrix44<T>&
4028 Matrix44<T>::invert (bool singExc)
4029 {
4030     *this = inverse (singExc);
4031     return *this;
4032 }
4033
4034 template <class T>
4035 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix44<T>&
4036 Matrix44<T>::invert() IMATH_NOEXCEPT
4037 {
4038     *this = inverse();
4039     return *this;
4040 }
4041
4042 template <class T>
4043 IMATH_CONSTEXPR14 inline Matrix44<T>
4044 Matrix44<T>::inverse (bool singExc) const
4045 {
4046     if (x[0][3] != 0 || x[1][3] != 0 || x[2][3] != 0 || x[3][3] != 1)
4047         return gjInverse (singExc);
4048
4049     Matrix44 s (x[1][1] * x[2][2] - x[2][1] * x[1][2],
4050                 x[2][1] * x[0][2] - x[0][1] * x[2][2],
4051                 x[0][1] * x[1][2] - x[1][1] * x[0][2],
4052                 0,
4053
4054                 x[2][0] * x[1][2] - x[1][0] * x[2][2],
4055                 x[0][0] * x[2][2] - x[2][0] * x[0][2],
4056                 x[1][0] * x[0][2] - x[0][0] * x[1][2],
4057                 0,
4058
4059                 x[1][0] * x[2][1] - x[2][0] * x[1][1],
4060                 x[2][0] * x[0][1] - x[0][0] * x[2][1],
4061                 x[0][0] * x[1][1] - x[1][0] * x[0][1],
4062                 0,
4063
4064                 0,
4065                 0,
4066                 0,
4067                 1);
4068
4069     T r = x[0][0] * s.x[0][0] + x[0][1] * s.x[1][0] + x[0][2] * s.x[2][0];
4070
4071     if (IMATH_INTERNAL_NAMESPACE::abs (r) >= 1)
4072     {
4073         for (int i = 0; i < 3; ++i)
4074         {
4075             for (int j = 0; j < 3; ++j)
4076             {
4077                 s.x[i][j] /= r;
4078             }
4079         }
4080     }
4081     else
4082     {
4083         T mr = IMATH_INTERNAL_NAMESPACE::abs (r) / std::numeric_limits<T>::min();
4084
4085         for (int i = 0; i < 3; ++i)
4086         {
4087             for (int j = 0; j < 3; ++j)
4088             {
4089                 if (mr > IMATH_INTERNAL_NAMESPACE::abs (s.x[i][j]))
4090                 {
4091                     s.x[i][j] /= r;
4092                 }
4093                 else
4094                 {
4095                     if (singExc)
4096                         throw std::invalid_argument ("Cannot invert singular matrix.");
4097
4098                     return Matrix44();
4099                 }
4100             }
4101         }
4102     }
4103
4104     s.x[3][0] = -x[3][0] * s.x[0][0] - x[3][1] * s.x[1][0] - x[3][2] * s.x[2][0];
4105     s.x[3][1] = -x[3][0] * s.x[0][1] - x[3][1] * s.x[1][1] - x[3][2] * s.x[2][1];
4106     s.x[3][2] = -x[3][0] * s.x[0][2] - x[3][1] * s.x[1][2] - x[3][2] * s.x[2][2];
4107
4108     return s;
4109 }
4110
4111 template <class T>
4112 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Matrix44<T>
4113 Matrix44<T>::inverse() const IMATH_NOEXCEPT
4114 {
4115     if (x[0][3] != 0 || x[1][3] != 0 || x[2][3] != 0 || x[3][3] != 1)
4116         return gjInverse();
4117
4118     Matrix44 s (x[1][1] * x[2][2] - x[2][1] * x[1][2],
4119                 x[2][1] * x[0][2] - x[0][1] * x[2][2],
4120                 x[0][1] * x[1][2] - x[1][1] * x[0][2],
4121                 0,
4122
4123                 x[2][0] * x[1][2] - x[1][0] * x[2][2],
4124                 x[0][0] * x[2][2] - x[2][0] * x[0][2],
4125                 x[1][0] * x[0][2] - x[0][0] * x[1][2],
4126                 0,
4127
4128                 x[1][0] * x[2][1] - x[2][0] * x[1][1],
4129                 x[2][0] * x[0][1] - x[0][0] * x[2][1],
4130                 x[0][0] * x[1][1] - x[1][0] * x[0][1],
4131                 0,
4132
4133                 0,
4134                 0,
4135                 0,
4136                 1);
4137
4138     T r = x[0][0] * s.x[0][0] + x[0][1] * s.x[1][0] + x[0][2] * s.x[2][0];
4139
4140     if (IMATH_INTERNAL_NAMESPACE::abs (r) >= 1)
4141     {
4142         for (int i = 0; i < 3; ++i)
4143         {
4144             for (int j = 0; j < 3; ++j)
4145             {
4146                 s.x[i][j] /= r;
4147             }
4148         }
4149     }
4150     else
4151     {
4152         T mr = IMATH_INTERNAL_NAMESPACE::abs (r) / std::numeric_limits<T>::min();
4153
4154         for (int i = 0; i < 3; ++i)
4155         {
4156             for (int j = 0; j < 3; ++j)
4157             {
4158                 if (mr > IMATH_INTERNAL_NAMESPACE::abs (s.x[i][j]))
4159                 {
4160                     s.x[i][j] /= r;
4161                 }
4162                 else
4163                 {
4164                     return Matrix44();
4165                 }
4166             }
4167         }
4168     }
4169
4170     s.x[3][0] = -x[3][0] * s.x[0][0] - x[3][1] * s.x[1][0] - x[3][2] * s.x[2][0];
4171     s.x[3][1] = -x[3][0] * s.x[0][1] - x[3][1] * s.x[1][1] - x[3][2] * s.x[2][1];
4172     s.x[3][2] = -x[3][0] * s.x[0][2] - x[3][1] * s.x[1][2] - x[3][2] * s.x[2][2];
4173
4174     return s;
4175 }
4176
4177 template <class T>
4178 IMATH_HOSTDEVICE constexpr inline T
4179 Matrix44<T>::fastMinor (const int r0,
4180                         const int r1,
4181                         const int r2,
4182                         const int c0,
4183                         const int c1,
4184                         const int c2) const IMATH_NOEXCEPT
4185 {
4186     return x[r0][c0] * (x[r1][c1] * x[r2][c2] - x[r1][c2] * x[r2][c1]) +
4187            x[r0][c1] * (x[r1][c2] * x[r2][c0] - x[r1][c0] * x[r2][c2]) +
4188            x[r0][c2] * (x[r1][c0] * x[r2][c1] - x[r1][c1] * x[r2][c0]);
4189 }
4190
4191 template <class T>
4192 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline T
4193 Matrix44<T>::minorOf (const int r, const int c) const IMATH_NOEXCEPT
4194 {
4195     int r0 = 0 + (r < 1 ? 1 : 0);
4196     int r1 = 1 + (r < 2 ? 1 : 0);
4197     int r2 = 2 + (r < 3 ? 1 : 0);
4198     int c0 = 0 + (c < 1 ? 1 : 0);
4199     int c1 = 1 + (c < 2 ? 1 : 0);
4200     int c2 = 2 + (c < 3 ? 1 : 0);
4201
4202     Matrix33<T> working (x[r0][c0],
4203                          x[r1][c0],
4204                          x[r2][c0],
4205                          x[r0][c1],
4206                          x[r1][c1],
4207                          x[r2][c1],
4208                          x[r0][c2],
4209                          x[r1][c2],
4210                          x[r2][c2]);
4211
4212     return working.determinant();
4213 }
4214
4215 template <class T>
4216 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline T
4217 Matrix44<T>::determinant() const IMATH_NOEXCEPT
4218 {
4219     T sum = (T) 0;
4220
4221     if (x[0][3] != 0.)
4222         sum -= x[0][3] * fastMinor (1, 2, 3, 0, 1, 2);
4223     if (x[1][3] != 0.)
4224         sum += x[1][3] * fastMinor (0, 2, 3, 0, 1, 2);
4225     if (x[2][3] != 0.)
4226         sum -= x[2][3] * fastMinor (0, 1, 3, 0, 1, 2);
4227     if (x[3][3] != 0.)
4228         sum += x[3][3] * fastMinor (0, 1, 2, 0, 1, 2);
4229
4230     return sum;
4231 }
4232
4233 template <class T>
4234 IMATH_HOSTDEVICE constexpr inline T
4235 Matrix44<T>::trace () const IMATH_NOEXCEPT
4236 {
4237     return x[0][0] + x[1][1] + x[2][2] + x[3][3];
4238 }
4239
4240 template <class T>
4241 template <class S>
4242 IMATH_HOSTDEVICE inline const Matrix44<T>&
4243 Matrix44<T>::setEulerAngles (const Vec3<S>& r) IMATH_NOEXCEPT
4244 {
4245     S cos_rz, sin_rz, cos_ry, sin_ry, cos_rx, sin_rx;
4246
4247     cos_rz = cos ((T) r.z);
4248     cos_ry = cos ((T) r.y);
4249     cos_rx = cos ((T) r.x);
4250
4251     sin_rz = sin ((T) r.z);
4252     sin_ry = sin ((T) r.y);
4253     sin_rx = sin ((T) r.x);
4254
4255     x[0][0] = cos_rz * cos_ry;
4256     x[0][1] = sin_rz * cos_ry;
4257     x[0][2] = -sin_ry;
4258     x[0][3] = 0;
4259
4260     x[1][0] = -sin_rz * cos_rx + cos_rz * sin_ry * sin_rx;
4261     x[1][1] = cos_rz * cos_rx + sin_rz * sin_ry * sin_rx;
4262     x[1][2] = cos_ry * sin_rx;
4263     x[1][3] = 0;
4264
4265     x[2][0] = sin_rz * sin_rx + cos_rz * sin_ry * cos_rx;
4266     x[2][1] = -cos_rz * sin_rx + sin_rz * sin_ry * cos_rx;
4267     x[2][2] = cos_ry * cos_rx;
4268     x[2][3] = 0;
4269
4270     x[3][0] = 0;
4271     x[3][1] = 0;
4272     x[3][2] = 0;
4273     x[3][3] = 1;
4274
4275     return *this;
4276 }
4277
4278 template <class T>
4279 template <class S>
4280 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix44<T>&
4281 Matrix44<T>::setAxisAngle (const Vec3<S>& axis, S angle) IMATH_NOEXCEPT
4282 {
4283     Vec3<S> unit (axis.normalized());
4284     S sine   = std::sin (angle);
4285     S cosine = std::cos (angle);
4286
4287     x[0][0] = unit.x * unit.x * (1 - cosine) + cosine;
4288     x[0][1] = unit.x * unit.y * (1 - cosine) + unit.z * sine;
4289     x[0][2] = unit.x * unit.z * (1 - cosine) - unit.y * sine;
4290     x[0][3] = 0;
4291
4292     x[1][0] = unit.x * unit.y * (1 - cosine) - unit.z * sine;
4293     x[1][1] = unit.y * unit.y * (1 - cosine) + cosine;
4294     x[1][2] = unit.y * unit.z * (1 - cosine) + unit.x * sine;
4295     x[1][3] = 0;
4296
4297     x[2][0] = unit.x * unit.z * (1 - cosine) + unit.y * sine;
4298     x[2][1] = unit.y * unit.z * (1 - cosine) - unit.x * sine;
4299     x[2][2] = unit.z * unit.z * (1 - cosine) + cosine;
4300     x[2][3] = 0;
4301
4302     x[3][0] = 0;
4303     x[3][1] = 0;
4304     x[3][2] = 0;
4305     x[3][3] = 1;
4306
4307     return *this;
4308 }
4309
4310 template <class T>
4311 template <class S>
4312 IMATH_HOSTDEVICE inline const Matrix44<T>&
4313 Matrix44<T>::rotate (const Vec3<S>& r) IMATH_NOEXCEPT
4314 {
4315     S cos_rz, sin_rz, cos_ry, sin_ry, cos_rx, sin_rx;
4316     S m00, m01, m02;
4317     S m10, m11, m12;
4318     S m20, m21, m22;
4319
4320     cos_rz = cos ((S) r.z);
4321     cos_ry = cos ((S) r.y);
4322     cos_rx = cos ((S) r.x);
4323
4324     sin_rz = sin ((S) r.z);
4325     sin_ry = sin ((S) r.y);
4326     sin_rx = sin ((S) r.x);
4327
4328     m00 = cos_rz * cos_ry;
4329     m01 = sin_rz * cos_ry;
4330     m02 = -sin_ry;
4331     m10 = -sin_rz * cos_rx + cos_rz * sin_ry * sin_rx;
4332     m11 = cos_rz * cos_rx + sin_rz * sin_ry * sin_rx;
4333     m12 = cos_ry * sin_rx;
4334     m20 = -sin_rz * -sin_rx + cos_rz * sin_ry * cos_rx;
4335     m21 = cos_rz * -sin_rx + sin_rz * sin_ry * cos_rx;
4336     m22 = cos_ry * cos_rx;
4337
4338     Matrix44<T> P (*this);
4339
4340     x[0][0] = P.x[0][0] * m00 + P.x[1][0] * m01 + P.x[2][0] * m02;
4341     x[0][1] = P.x[0][1] * m00 + P.x[1][1] * m01 + P.x[2][1] * m02;
4342     x[0][2] = P.x[0][2] * m00 + P.x[1][2] * m01 + P.x[2][2] * m02;
4343     x[0][3] = P.x[0][3] * m00 + P.x[1][3] * m01 + P.x[2][3] * m02;
4344
4345     x[1][0] = P.x[0][0] * m10 + P.x[1][0] * m11 + P.x[2][0] * m12;
4346     x[1][1] = P.x[0][1] * m10 + P.x[1][1] * m11 + P.x[2][1] * m12;
4347     x[1][2] = P.x[0][2] * m10 + P.x[1][2] * m11 + P.x[2][2] * m12;
4348     x[1][3] = P.x[0][3] * m10 + P.x[1][3] * m11 + P.x[2][3] * m12;
4349
4350     x[2][0] = P.x[0][0] * m20 + P.x[1][0] * m21 + P.x[2][0] * m22;
4351     x[2][1] = P.x[0][1] * m20 + P.x[1][1] * m21 + P.x[2][1] * m22;
4352     x[2][2] = P.x[0][2] * m20 + P.x[1][2] * m21 + P.x[2][2] * m22;
4353     x[2][3] = P.x[0][3] * m20 + P.x[1][3] * m21 + P.x[2][3] * m22;
4354
4355     return *this;
4356 }
4357
4358 template <class T>
4359 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix44<T>&
4360 Matrix44<T>::setScale (T s) IMATH_NOEXCEPT
4361 {
4362     //
4363     // Set the matrix to a 3D homogeneous transform scale:
4364     //  | s 0 0 0 |
4365     //  | 0 s 0 0 |
4366     //  | 0 0 s 0 |
4367     //  | 0 0 0 1 |
4368     //
4369
4370     x[0][0] = s;
4371     x[0][1] = 0;
4372     x[0][2] = 0;
4373     x[0][3] = 0;
4374     x[1][0] = 0;
4375     x[1][1] = s;
4376     x[1][2] = 0;
4377     x[1][3] = 0;
4378     x[2][0] = 0;
4379     x[2][1] = 0;
4380     x[2][2] = s;
4381     x[2][3] = 0;
4382     x[3][0] = 0;
4383     x[3][1] = 0;
4384     x[3][2] = 0;
4385     x[3][3] = 1;
4386     return *this;
4387 }
4388
4389 template <class T>
4390 template <class S>
4391 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix44<T>&
4392 Matrix44<T>::setScale (const Vec3<S>& s) IMATH_NOEXCEPT
4393 {
4394     //
4395     // Set the matrix to a 3D homogeneous transform scale:
4396     //  | s.x  0   0   0 |
4397     //  |  0  s.y  0   0 |
4398     //  |  0   0  s.z  0 |
4399     //  |  0   0   0   1 |
4400     //
4401
4402     x[0][0] = s.x;
4403     x[0][1] = 0;
4404     x[0][2] = 0;
4405     x[0][3] = 0;
4406     x[1][0] = 0;
4407     x[1][1] = s.y;
4408     x[1][2] = 0;
4409     x[1][3] = 0;
4410     x[2][0] = 0;
4411     x[2][1] = 0;
4412     x[2][2] = s.z;
4413     x[2][3] = 0;
4414     x[3][0] = 0;
4415     x[3][1] = 0;
4416     x[3][2] = 0;
4417     x[3][3] = 1;
4418     return *this;
4419 }
4420
4421 template <class T>
4422 template <class S>
4423 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix44<T>&
4424 Matrix44<T>::scale (const Vec3<S>& s) IMATH_NOEXCEPT
4425 {
4426     x[0][0] *= s.x;
4427     x[0][1] *= s.x;
4428     x[0][2] *= s.x;
4429     x[0][3] *= s.x;
4430
4431     x[1][0] *= s.y;
4432     x[1][1] *= s.y;
4433     x[1][2] *= s.y;
4434     x[1][3] *= s.y;
4435
4436     x[2][0] *= s.z;
4437     x[2][1] *= s.z;
4438     x[2][2] *= s.z;
4439     x[2][3] *= s.z;
4440
4441     return *this;
4442 }
4443
4444 template <class T>
4445 template <class S>
4446 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix44<T>&
4447 Matrix44<T>::setTranslation (const Vec3<S>& t) IMATH_NOEXCEPT
4448 {
4449     x[0][0] = 1;
4450     x[0][1] = 0;
4451     x[0][2] = 0;
4452     x[0][3] = 0;
4453
4454     x[1][0] = 0;
4455     x[1][1] = 1;
4456     x[1][2] = 0;
4457     x[1][3] = 0;
4458
4459     x[2][0] = 0;
4460     x[2][1] = 0;
4461     x[2][2] = 1;
4462     x[2][3] = 0;
4463
4464     x[3][0] = t.x;
4465     x[3][1] = t.y;
4466     x[3][2] = t.z;
4467     x[3][3] = 1;
4468
4469     return *this;
4470 }
4471
4472 template <class T>
4473 IMATH_HOSTDEVICE constexpr inline const Vec3<T>
4474 Matrix44<T>::translation() const IMATH_NOEXCEPT
4475 {
4476     return Vec3<T> (x[3][0], x[3][1], x[3][2]);
4477 }
4478
4479 template <class T>
4480 template <class S>
4481 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix44<T>&
4482 Matrix44<T>::translate (const Vec3<S>& t) IMATH_NOEXCEPT
4483 {
4484     x[3][0] += t.x * x[0][0] + t.y * x[1][0] + t.z * x[2][0];
4485     x[3][1] += t.x * x[0][1] + t.y * x[1][1] + t.z * x[2][1];
4486     x[3][2] += t.x * x[0][2] + t.y * x[1][2] + t.z * x[2][2];
4487     x[3][3] += t.x * x[0][3] + t.y * x[1][3] + t.z * x[2][3];
4488
4489     return *this;
4490 }
4491
4492 template <class T>
4493 template <class S>
4494 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix44<T>&
4495 Matrix44<T>::setShear (const Vec3<S>& h) IMATH_NOEXCEPT
4496 {
4497     x[0][0] = 1;
4498     x[0][1] = 0;
4499     x[0][2] = 0;
4500     x[0][3] = 0;
4501
4502     x[1][0] = h.x;
4503     x[1][1] = 1;
4504     x[1][2] = 0;
4505     x[1][3] = 0;
4506
4507     x[2][0] = h.y;
4508     x[2][1] = h.z;
4509     x[2][2] = 1;
4510     x[2][3] = 0;
4511
4512     x[3][0] = 0;
4513     x[3][1] = 0;
4514     x[3][2] = 0;
4515     x[3][3] = 1;
4516
4517     return *this;
4518 }
4519
4520 template <class T>
4521 template <class S>
4522 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix44<T>&
4523 Matrix44<T>::setShear (const Shear6<S>& h) IMATH_NOEXCEPT
4524 {
4525     x[0][0] = 1;
4526     x[0][1] = h.yx;
4527     x[0][2] = h.zx;
4528     x[0][3] = 0;
4529
4530     x[1][0] = h.xy;
4531     x[1][1] = 1;
4532     x[1][2] = h.zy;
4533     x[1][3] = 0;
4534
4535     x[2][0] = h.xz;
4536     x[2][1] = h.yz;
4537     x[2][2] = 1;
4538     x[2][3] = 0;
4539
4540     x[3][0] = 0;
4541     x[3][1] = 0;
4542     x[3][2] = 0;
4543     x[3][3] = 1;
4544
4545     return *this;
4546 }
4547
4548 template <class T>
4549 template <class S>
4550 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix44<T>&
4551 Matrix44<T>::shear (const Vec3<S>& h) IMATH_NOEXCEPT
4552 {
4553     //
4554     // In this case, we don't need a temp. copy of the matrix
4555     // because we never use a value on the RHS after we've
4556     // changed it on the LHS.
4557     //
4558
4559     for (int i = 0; i < 4; i++)
4560     {
4561         x[2][i] += h.y * x[0][i] + h.z * x[1][i];
4562         x[1][i] += h.x * x[0][i];
4563     }
4564
4565     return *this;
4566 }
4567
4568 template <class T>
4569 template <class S>
4570 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Matrix44<T>&
4571 Matrix44<T>::shear (const Shear6<S>& h) IMATH_NOEXCEPT
4572 {
4573     Matrix44<T> P (*this);
4574
4575     for (int i = 0; i < 4; i++)
4576     {
4577         x[0][i] = P.x[0][i] + h.yx * P.x[1][i] + h.zx * P.x[2][i];
4578         x[1][i] = h.xy * P.x[0][i] + P.x[1][i] + h.zy * P.x[2][i];
4579         x[2][i] = h.xz * P.x[0][i] + h.yz * P.x[1][i] + P.x[2][i];
4580     }
4581
4582     return *this;
4583 }
4584
4585 //--------------------------------
4586 // Implementation of stream output
4587 //--------------------------------
4588
4589 template <class T>
4590 std::ostream&
4591 operator<< (std::ostream& s, const Matrix22<T>& m)
4592 {
4593     std::ios_base::fmtflags oldFlags = s.flags();
4594     int width;
4595
4596     if (s.flags() & std::ios_base::fixed)
4597     {
4598         s.setf (std::ios_base::showpoint);
4599         width = static_cast<int> (s.precision()) + 5;
4600     }
4601     else
4602     {
4603         s.setf (std::ios_base::scientific);
4604         s.setf (std::ios_base::showpoint);
4605         width = static_cast<int> (s.precision()) + 8;
4606     }
4607
4608     s << "(" << std::setw (width) << m[0][0] << " " << std::setw (width) << m[0][1] << "\n"
4609       <<
4610
4611         " " << std::setw (width) << m[1][0] << " " << std::setw (width) << m[1][1] << ")\n";
4612
4613     s.flags (oldFlags);
4614     return s;
4615 }
4616
4617 template <class T>
4618 std::ostream&
4619 operator<< (std::ostream& s, const Matrix33<T>& m)
4620 {
4621     std::ios_base::fmtflags oldFlags = s.flags();
4622     int width;
4623
4624     if (s.flags() & std::ios_base::fixed)
4625     {
4626         s.setf (std::ios_base::showpoint);
4627         width = static_cast<int> (s.precision()) + 5;
4628     }
4629     else
4630     {
4631         s.setf (std::ios_base::scientific);
4632         s.setf (std::ios_base::showpoint);
4633         width = static_cast<int> (s.precision()) + 8;
4634     }
4635
4636     s << "(" << std::setw (width) << m[0][0] << " " << std::setw (width) << m[0][1] << " "
4637       << std::setw (width) << m[0][2] << "\n"
4638       <<
4639
4640         " " << std::setw (width) << m[1][0] << " " << std::setw (width) << m[1][1] << " "
4641       << std::setw (width) << m[1][2] << "\n"
4642       <<
4643
4644         " " << std::setw (width) << m[2][0] << " " << std::setw (width) << m[2][1] << " "
4645       << std::setw (width) << m[2][2] << ")\n";
4646
4647     s.flags (oldFlags);
4648     return s;
4649 }
4650
4651 template <class T>
4652 std::ostream&
4653 operator<< (std::ostream& s, const Matrix44<T>& m)
4654 {
4655     std::ios_base::fmtflags oldFlags = s.flags();
4656     int width;
4657
4658     if (s.flags() & std::ios_base::fixed)
4659     {
4660         s.setf (std::ios_base::showpoint);
4661         width = static_cast<int> (s.precision()) + 5;
4662     }
4663     else
4664     {
4665         s.setf (std::ios_base::scientific);
4666         s.setf (std::ios_base::showpoint);
4667         width = static_cast<int> (s.precision()) + 8;
4668     }
4669
4670     s << "(" << std::setw (width) << m[0][0] << " " << std::setw (width) << m[0][1] << " "
4671       << std::setw (width) << m[0][2] << " " << std::setw (width) << m[0][3] << "\n"
4672       <<
4673
4674         " " << std::setw (width) << m[1][0] << " " << std::setw (width) << m[1][1] << " "
4675       << std::setw (width) << m[1][2] << " " << std::setw (width) << m[1][3] << "\n"
4676       <<
4677
4678         " " << std::setw (width) << m[2][0] << " " << std::setw (width) << m[2][1] << " "
4679       << std::setw (width) << m[2][2] << " " << std::setw (width) << m[2][3] << "\n"
4680       <<
4681
4682         " " << std::setw (width) << m[3][0] << " " << std::setw (width) << m[3][1] << " "
4683       << std::setw (width) << m[3][2] << " " << std::setw (width) << m[3][3] << ")\n";
4684
4685     s.flags (oldFlags);
4686     return s;
4687 }
4688
4689 //---------------------------------------------------------------
4690 // Implementation of vector-times-matrix multiplication operators
4691 //---------------------------------------------------------------
4692
4693 template <class S, class T>
4694 IMATH_HOSTDEVICE inline const Vec2<S>&
4695 operator*= (Vec2<S>& v, const Matrix22<T>& m) IMATH_NOEXCEPT
4696 {
4697     S x = S (v.x * m.x[0][0] + v.y * m.x[1][0]);
4698     S y = S (v.x * m.x[0][1] + v.y * m.x[1][1]);
4699
4700     v.x = x;
4701     v.y = y;
4702
4703     return v;
4704 }
4705
4706 template <class S, class T>
4707 IMATH_HOSTDEVICE inline Vec2<S>
4708 operator* (const Vec2<S>& v, const Matrix22<T>& m) IMATH_NOEXCEPT
4709 {
4710     S x = S (v.x * m.x[0][0] + v.y * m.x[1][0]);
4711     S y = S (v.x * m.x[0][1] + v.y * m.x[1][1]);
4712
4713     return Vec2<S> (x, y);
4714 }
4715
4716 template <class S, class T>
4717 IMATH_HOSTDEVICE inline const Vec2<S>&
4718 operator*= (Vec2<S>& v, const Matrix33<T>& m) IMATH_NOEXCEPT
4719 {
4720     S x = S (v.x * m.x[0][0] + v.y * m.x[1][0] + m.x[2][0]);
4721     S y = S (v.x * m.x[0][1] + v.y * m.x[1][1] + m.x[2][1]);
4722     S w = S (v.x * m.x[0][2] + v.y * m.x[1][2] + m.x[2][2]);
4723
4724     v.x = x / w;
4725     v.y = y / w;
4726
4727     return v;
4728 }
4729
4730 template <class S, class T>
4731 IMATH_HOSTDEVICE inline Vec2<S>
4732 operator* (const Vec2<S>& v, const Matrix33<T>& m) IMATH_NOEXCEPT
4733 {
4734     S x = S (v.x * m.x[0][0] + v.y * m.x[1][0] + m.x[2][0]);
4735     S y = S (v.x * m.x[0][1] + v.y * m.x[1][1] + m.x[2][1]);
4736     S w = S (v.x * m.x[0][2] + v.y * m.x[1][2] + m.x[2][2]);
4737
4738     return Vec2<S> (x / w, y / w);
4739 }
4740
4741 template <class S, class T>
4742 IMATH_HOSTDEVICE inline const Vec3<S>&
4743 operator*= (Vec3<S>& v, const Matrix33<T>& m) IMATH_NOEXCEPT
4744 {
4745     S x = S (v.x * m.x[0][0] + v.y * m.x[1][0] + v.z * m.x[2][0]);
4746     S y = S (v.x * m.x[0][1] + v.y * m.x[1][1] + v.z * m.x[2][1]);
4747     S z = S (v.x * m.x[0][2] + v.y * m.x[1][2] + v.z * m.x[2][2]);
4748
4749     v.x = x;
4750     v.y = y;
4751     v.z = z;
4752
4753     return v;
4754 }
4755
4756 template <class S, class T>
4757 IMATH_HOSTDEVICE inline Vec3<S>
4758 operator* (const Vec3<S>& v, const Matrix33<T>& m) IMATH_NOEXCEPT
4759 {
4760     S x = S (v.x * m.x[0][0] + v.y * m.x[1][0] + v.z * m.x[2][0]);
4761     S y = S (v.x * m.x[0][1] + v.y * m.x[1][1] + v.z * m.x[2][1]);
4762     S z = S (v.x * m.x[0][2] + v.y * m.x[1][2] + v.z * m.x[2][2]);
4763
4764     return Vec3<S> (x, y, z);
4765 }
4766
4767 template <class S, class T>
4768 IMATH_HOSTDEVICE inline const Vec3<S>&
4769 operator*= (Vec3<S>& v, const Matrix44<T>& m) IMATH_NOEXCEPT
4770 {
4771     S x = S (v.x * m.x[0][0] + v.y * m.x[1][0] + v.z * m.x[2][0] + m.x[3][0]);
4772     S y = S (v.x * m.x[0][1] + v.y * m.x[1][1] + v.z * m.x[2][1] + m.x[3][1]);
4773     S z = S (v.x * m.x[0][2] + v.y * m.x[1][2] + v.z * m.x[2][2] + m.x[3][2]);
4774     S w = S (v.x * m.x[0][3] + v.y * m.x[1][3] + v.z * m.x[2][3] + m.x[3][3]);
4775
4776     v.x = x / w;
4777     v.y = y / w;
4778     v.z = z / w;
4779
4780     return v;
4781 }
4782
4783 template <class S, class T>
4784 IMATH_HOSTDEVICE inline Vec3<S>
4785 IMATH_HOSTDEVICE operator* (const Vec3<S>& v, const Matrix44<T>& m) IMATH_NOEXCEPT
4786 {
4787     S x = S (v.x * m.x[0][0] + v.y * m.x[1][0] + v.z * m.x[2][0] + m.x[3][0]);
4788     S y = S (v.x * m.x[0][1] + v.y * m.x[1][1] + v.z * m.x[2][1] + m.x[3][1]);
4789     S z = S (v.x * m.x[0][2] + v.y * m.x[1][2] + v.z * m.x[2][2] + m.x[3][2]);
4790     S w = S (v.x * m.x[0][3] + v.y * m.x[1][3] + v.z * m.x[2][3] + m.x[3][3]);
4791
4792     return Vec3<S> (x / w, y / w, z / w);
4793 }
4794
4795 template <class S, class T>
4796 IMATH_HOSTDEVICE inline const Vec4<S>&
4797 IMATH_HOSTDEVICE operator*= (Vec4<S>& v, const Matrix44<T>& m) IMATH_NOEXCEPT
4798 {
4799     S x = S (v.x * m.x[0][0] + v.y * m.x[1][0] + v.z * m.x[2][0] + v.w * m.x[3][0]);
4800     S y = S (v.x * m.x[0][1] + v.y * m.x[1][1] + v.z * m.x[2][1] + v.w * m.x[3][1]);
4801     S z = S (v.x * m.x[0][2] + v.y * m.x[1][2] + v.z * m.x[2][2] + v.w * m.x[3][2]);
4802     S w = S (v.x * m.x[0][3] + v.y * m.x[1][3] + v.z * m.x[2][3] + v.w * m.x[3][3]);
4803
4804     v.x = x;
4805     v.y = y;
4806     v.z = z;
4807     v.w = w;
4808
4809     return v;
4810 }
4811
4812 template <class S, class T>
4813 IMATH_HOSTDEVICE inline Vec4<S>
4814 IMATH_HOSTDEVICE operator* (const Vec4<S>& v, const Matrix44<T>& m) IMATH_NOEXCEPT
4815 {
4816     S x = S (v.x * m.x[0][0] + v.y * m.x[1][0] + v.z * m.x[2][0] + v.w * m.x[3][0]);
4817     S y = S (v.x * m.x[0][1] + v.y * m.x[1][1] + v.z * m.x[2][1] + v.w * m.x[3][1]);
4818     S z = S (v.x * m.x[0][2] + v.y * m.x[1][2] + v.z * m.x[2][2] + v.w * m.x[3][2]);
4819     S w = S (v.x * m.x[0][3] + v.y * m.x[1][3] + v.z * m.x[2][3] + v.w * m.x[3][3]);
4820
4821     return Vec4<S> (x, y, z, w);
4822 }
4823
4824 IMATH_INTERNAL_NAMESPACE_HEADER_EXIT
4825
4826 #endif // INCLUDED_IMATHMATRIX_H