util/X86MAC64/cuda/samples/5_Simulations/smokeParticles/nvMatrix.h

   1 /*
   2  * Copyright 1993-2013 NVIDIA Corporation.  All rights reserved.
   3  *
   4  * Please refer to the NVIDIA end user license agreement (EULA) associated
   5  * with this source code for terms and conditions that govern your use of
   6  * this software. Any use, reproduction, disclosure, or distribution of
   7  * this software and related documentation outside the terms of the EULA
   8  * is strictly prohibited.
   9  *
  10  */
  11
  12 //
  13 // Template math library for common 3D functionality
  14 //
  15 // nvMatrix.h - template matrix code
  16 //
  17 // This code is in part deriver from glh, a cross platform glut helper library.
  18 // The copyright for glh follows this notice.
  19 //
  20 // Copyright (c) NVIDIA Corporation. All rights reserved.
  21 ////////////////////////////////////////////////////////////////////////////////
  22
  23 /*
  24     Copyright (c) 2000 Cass Everitt
  25     Copyright (c) 2000 NVIDIA Corporation
  26     All rights reserved.
  27
  28     Redistribution and use in source and binary forms, with or
  29     without modification, are permitted provided that the following
  30     conditions are met:
  31
  32      * Redistributions of source code must retain the above
  33        copyright notice, this list of conditions and the following
  34        disclaimer.
  35
  36      * Redistributions in binary form must reproduce the above
  37        copyright notice, this list of conditions and the following
  38        disclaimer in the documentation and/or other materials
  39        provided with the distribution.
  40
  41      * The names of contributors to this software may not be used
  42        to endorse or promote products derived from this software
  43        without specific prior written permission.
  44
  45        THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  46        ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  47        LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
  48        FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
  49        REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
  50        INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
  51        BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
  52        LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  53        CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  54        LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
  55        ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  56        POSSIBILITY OF SUCH DAMAGE.
  57
  58
  59     Cass Everitt - cass@r3.nu
  60 */
  61
  62 #ifndef NV_MATRIX_H
  63 #define NV_MATRIX_H
  64
  65 namespace nv
  66 {
  67
  68     template <class T> class vec2;
  69     template <class T> class vec3;
  70     template <class T> class vec4;
  71
  72     ////////////////////////////////////////////////////////////////////////////////
  73     //
  74     //  Matrix
  75     //
  76     ////////////////////////////////////////////////////////////////////////////////
  77     template<class T>
  78     class matrix4
  79     {
  80
  81         public:
  82
  83             matrix4()
  84             {
  85                 make_identity();
  86             }
  87
  88             matrix4(T t)
  89             {
  90                 set_value(t);
  91             }
  92
  93             matrix4(const T *m)
  94             {
  95                 set_value(m);
  96             }
  97
  98             matrix4(T a00, T a01, T a02, T a03,
  99                     T a10, T a11, T a12, T a13,
 100                     T a20, T a21, T a22, T a23,
 101                     T a30, T a31, T a32, T a33) :
 102                 _11(a00), _12(a01), _13(a02), _14(a03),
 103                 _21(a10), _22(a11), _23(a12), _24(a13),
 104                 _31(a20), _32(a21), _33(a22), _34(a23),
 105                 _41(a30), _42(a31), _43(a32), _44(a33)
 106             {}
 107
 108
 109             void get_value(T *mp) const
 110             {
 111                 int c = 0;
 112
 113                 for (int j=0; j < 4; j++)
 114                     for (int i=0; i < 4; i++)
 115                     {
 116                         mp[c++] = element(i,j);
 117                     }
 118             }
 119
 120             const T *get_value() const
 121             {
 122                 return _array;
 123             }
 124
 125             void set_value(T *mp)
 126             {
 127                 int c = 0;
 128
 129                 for (int j=0; j < 4; j++)
 130                     for (int i=0; i < 4; i++)
 131                     {
 132                         element(i,j) = mp[c++];
 133                     }
 134             }
 135
 136             void set_value(T r)
 137             {
 138                 for (int i=0; i < 4; i++)
 139                     for (int j=0; j < 4; j++)
 140                     {
 141                         element(i,j) = r;
 142                     }
 143             }
 144
 145             void make_identity()
 146             {
 147                 element(0,0) = 1.0;
 148                 element(0,1) = 0.0;
 149                 element(0,2) = 0.0;
 150                 element(0,3) = 0.0;
 151
 152                 element(1,0) = 0.0;
 153                 element(1,1) = 1.0;
 154                 element(1,2) = 0.0;
 155                 element(1,3) = 0.0;
 156
 157                 element(2,0) = 0.0;
 158                 element(2,1) = 0.0;
 159                 element(2,2) = 1.0;
 160                 element(2,3) = 0.0;
 161
 162                 element(3,0) = 0.0;
 163                 element(3,1) = 0.0;
 164                 element(3,2) = 0.0;
 165                 element(3,3) = 1.0;
 166             }
 167
 168             // set a uniform scale
 169             void set_scale(T s)
 170             {
 171                 element(0,0) = s;
 172                 element(1,1) = s;
 173                 element(2,2) = s;
 174             }
 175
 176             void set_scale(const vec3<T> &s)
 177             {
 178                 for (int i = 0; i < 3; i++)
 179                 {
 180                     element(i,i) = s[i];
 181                 }
 182             }
 183
 184
 185             void set_translate(const vec3<T> &t)
 186             {
 187                 for (int i = 0; i < 3; i++)
 188                 {
 189                     element(i,3) = t[i];
 190                 }
 191             }
 192
 193             void set_row(int r, const vec4<T> &t)
 194             {
 195                 for (int i = 0; i < 4; i++)
 196                 {
 197                     element(r,i) = t[i];
 198                 }
 199             }
 200
 201             void set_column(int c, const vec4<T> &t)
 202             {
 203                 for (int i = 0; i < 4; i++)
 204                 {
 205                     element(i,c) = t[i];
 206                 }
 207             }
 208
 209             vec4<T> get_row(int r) const
 210             {
 211                 vec4<T> v;
 212
 213                 for (int i = 0; i < 4; i++)
 214                 {
 215                     v[i] = element(r,i);
 216                 }
 217
 218                 return v;
 219             }
 220
 221             vec4<T> get_column(int c) const
 222             {
 223                 vec4<T> v;
 224
 225                 for (int i = 0; i < 4; i++)
 226                 {
 227                     v[i] = element(i,c);
 228                 }
 229
 230                 return v;
 231             }
 232
 233             friend matrix4 inverse(const matrix4 &m)
 234             {
 235                 matrix4 minv;
 236
 237                 T r1[8], r2[8], r3[8], r4[8];
 238                 T *s[4], *tmprow;
 239
 240                 s[0] = &r1[0];
 241                 s[1] = &r2[0];
 242                 s[2] = &r3[0];
 243                 s[3] = &r4[0];
 244
 245                 register int i,j,p,jj;
 246
 247                 for (i=0; i<4; i++)
 248                 {
 249                     for (j=0; j<4; j++)
 250                     {
 251                         s[i][j] = m.element(i,j);
 252
 253                         if (i==j)
 254                         {
 255                             s[i][j+4] = 1.0;
 256                         }
 257                         else
 258                         {
 259                             s[i][j+4] = 0.0;
 260                         }
 261                     }
 262                 }
 263
 264                 T scp[4];
 265
 266                 for (i=0; i<4; i++)
 267                 {
 268                     scp[i] = T(fabs(s[i][0]));
 269
 270                     for (j=1; j<4; j++)
 271                         if (T(fabs(s[i][j])) > scp[i])
 272                         {
 273                             scp[i] = T(fabs(s[i][j]));
 274                         }
 275
 276                     if (scp[i] == 0.0)
 277                     {
 278                         return minv;    // singular matrix!
 279                     }
 280                 }
 281
 282                 int pivot_to;
 283                 T scp_max;
 284
 285                 for (i=0; i<4; i++)
 286                 {
 287                     // select pivot row
 288                     pivot_to = i;
 289                     scp_max = T(fabs(s[i][i]/scp[i]));
 290
 291                     // find out which row should be on top
 292                     for (p=i+1; p<4; p++)
 293                         if (T(fabs(s[p][i]/scp[p])) > scp_max)
 294                         {
 295                             scp_max = T(fabs(s[p][i]/scp[p]));
 296                             pivot_to = p;
 297                         }
 298
 299                     // Pivot if necessary
 300                     if (pivot_to != i)
 301                     {
 302                         tmprow = s[i];
 303                         s[i] = s[pivot_to];
 304                         s[pivot_to] = tmprow;
 305                         T tmpscp;
 306                         tmpscp = scp[i];
 307                         scp[i] = scp[pivot_to];
 308                         scp[pivot_to] = tmpscp;
 309                     }
 310
 311                     T mji;
 312
 313                     // perform gaussian elimination
 314                     for (j=i+1; j<4; j++)
 315                     {
 316                         mji = s[j][i]/s[i][i];
 317                         s[j][i] = 0.0;
 318
 319                         for (jj=i+1; jj<8; jj++)
 320                         {
 321                             s[j][jj] -= mji*s[i][jj];
 322                         }
 323                     }
 324                 }
 325
 326                 if (s[3][3] == 0.0)
 327                 {
 328                     return minv;    // singular matrix!
 329                 }
 330
 331                 //
 332                 // Now we have an upper triangular matrix.
 333                 //
 334                 //  x x x x | y y y y
 335                 //  0 x x x | y y y y
 336                 //  0 0 x x | y y y y
 337                 //  0 0 0 x | y y y y
 338                 //
 339                 //  we'll back substitute to get the inverse
 340                 //
 341                 //  1 0 0 0 | z z z z
 342                 //  0 1 0 0 | z z z z
 343                 //  0 0 1 0 | z z z z
 344                 //  0 0 0 1 | z z z z
 345                 //
 346
 347                 T mij;
 348
 349                 for (i=3; i>0; i--)
 350                 {
 351                     for (j=i-1; j > -1; j--)
 352                     {
 353                         mij = s[j][i]/s[i][i];
 354
 355                         for (jj=j+1; jj<8; jj++)
 356                         {
 357                             s[j][jj] -= mij*s[i][jj];
 358                         }
 359                     }
 360                 }
 361
 362                 for (i=0; i<4; i++)
 363                     for (j=0; j<4; j++)
 364                     {
 365                         minv(i,j) = s[i][j+4] / s[i][i];
 366                     }
 367
 368                 return minv;
 369             }
 370
 371
 372             friend matrix4 transpose(const matrix4 &m)
 373             {
 374                 matrix4 mtrans;
 375
 376                 for (int i=0; i<4; i++)
 377                     for (int j=0; j<4; j++)
 378                     {
 379                         mtrans(i,j) = m.element(j,i);
 380                     }
 381
 382                 return mtrans;
 383             }
 384
 385             matrix4 &operator *= (const matrix4 &rhs)
 386             {
 387                 matrix4 mt(*this);
 388                 set_value(T(0));
 389
 390                 for (int i=0; i < 4; i++)
 391                     for (int j=0; j < 4; j++)
 392                         for (int c=0; c < 4; c++)
 393                         {
 394                             element(i,j) += mt(i,c) * rhs(c,j);
 395                         }
 396
 397                 return *this;
 398             }
 399
 400             friend matrix4 operator * (const matrix4 &lhs, const matrix4 &rhs)
 401             {
 402                 matrix4 r(T(0));
 403
 404                 for (int i=0; i < 4; i++)
 405                     for (int j=0; j < 4; j++)
 406                         for (int c=0; c < 4; c++)
 407                         {
 408                             r.element(i,j) += lhs(i,c) * rhs(c,j);
 409                         }
 410
 411                 return r;
 412             }
 413
 414             // dst = M * src
 415             vec4<T> operator *(const vec4<T> &src) const
 416             {
 417                 vec4<T> r;
 418
 419                 for (int i = 0; i < 4; i++)
 420                     r[i]  = (src[0] * element(i,0) + src[1] * element(i,1) +
 421                              src[2] * element(i,2) + src[3] * element(i,3));
 422
 423                 return r;
 424             }
 425
 426             // dst = src * M
 427             friend vec4<T> operator *(const vec4<T> &lhs, const matrix4 &rhs)
 428             {
 429                 vec4<T> r;
 430
 431                 for (int i = 0; i < 4; i++)
 432                     r[i]  = (lhs[0] * rhs.element(0,i) + lhs[1] * rhs.element(1,i) +
 433                              lhs[2] * rhs.element(2,i) + lhs[3] * rhs.element(3,i));
 434
 435                 return r;
 436             }
 437
 438             T &operator()(int row, int col)
 439             {
 440                 return element(row,col);
 441             }
 442
 443             const T &operator()(int row, int col) const
 444             {
 445                 return element(row,col);
 446             }
 447
 448             T &element(int row, int col)
 449             {
 450                 return _array[row | (col<<2)];
 451             }
 452
 453             const T &element(int row, int col) const
 454             {
 455                 return _array[row | (col<<2)];
 456             }
 457
 458             matrix4 &operator *= (const T &r)
 459             {
 460                 for (int i = 0; i < 4; ++i)
 461                 {
 462                     element(0,i) *= r;
 463                     element(1,i) *= r;
 464                     element(2,i) *= r;
 465                     element(3,i) *= r;
 466                 }
 467
 468                 return *this;
 469             }
 470
 471             matrix4 &operator += (const matrix4 &mat)
 472             {
 473                 for (int i = 0; i < 4; ++i)
 474                 {
 475                     element(0,i) += mat.element(0,i);
 476                     element(1,i) += mat.element(1,i);
 477                     element(2,i) += mat.element(2,i);
 478                     element(3,i) += mat.element(3,i);
 479                 }
 480
 481                 return *this;
 482             }
 483
 484
 485             friend bool operator == (const matrix4 &lhs, const matrix4 &rhs)
 486             {
 487                 bool r = true;
 488
 489                 for (int i = 0; i < 16; i++)
 490                 {
 491                     r &= lhs._array[i] == rhs._array[i];
 492                 }
 493
 494                 return r;
 495             }
 496
 497             friend bool operator != (const matrix4 &lhs, const matrix4 &rhs)
 498             {
 499                 bool r = true;
 500
 501                 for (int i = 0; i < 16; i++)
 502                 {
 503                     r &= lhs._array[i] != rhs._array[i];
 504                 }
 505
 506                 return r;
 507             }
 508
 509             union
 510             {
 511                 struct
 512                 {
 513                     T _11, _12, _13, _14;   // standard names for components
 514                     T _21, _22, _23, _24;   // standard names for components
 515                     T _31, _32, _33, _34;   // standard names for components
 516                     T _41, _42, _43, _44;   // standard names for components
 517                 };
 518                 T _array[16];     // array access
 519             };
 520     };
 521
 522 };
 523
 524 #endif