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 (r0: 1, r1: 2, r2: 3, c0: 0, c1: 1, c2: 2);
4223	if (x[1][3] != 0.)
4224	sum += x[1][3] * fastMinor (r0: 0, r1: 2, r2: 3, c0: 0, c1: 1, c2: 2);
4225	if (x[2][3] != 0.)
4226	sum -= x[2][3] * fastMinor (r0: 0, r1: 1, r2: 3, c0: 0, c1: 1, c2: 2);
4227	if (x[3][3] != 0.)
4228	sum += x[3][3] * fastMinor (r0: 0, r1: 1, r2: 2, c0: 0, c1: 1, c2: 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 (fmtfl: 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 (fmtfl: 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 (fmtfl: 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
4827