Cleanup: remove redundant doxygen \file argument
[blender.git] / intern / smoke / intern / tnt / tnt_cmat.h
1 /** \file \ingroup smoke
2  */
3 /*
4 *
5 * Template Numerical Toolkit (TNT)
6 *
7 * Mathematical and Computational Sciences Division
8 * National Institute of Technology,
9 * Gaithersburg, MD USA
10 *
11 *
12 * This software was developed at the National Institute of Standards and
13 * Technology (NIST) by employees of the Federal Government in the course
14 * of their official duties. Pursuant to title 17 Section 105 of the
15 * United States Code, this software is not subject to copyright protection
16 * and is in the public domain. NIST assumes no responsibility whatsoever for
17 * its use by other parties, and makes no guarantees, expressed or implied,
18 * about its quality, reliability, or any other characteristic.
19 *
20 */
21
22
23 // C compatible matrix: row-oriented, 0-based [i][j] and 1-based (i,j) indexing
24 //
25
26 #ifndef TNT_CMAT_H
27 #define TNT_CMAT_H
28
29 #include "tnt_subscript.h"
30 #include "tnt_vec.h"
31 #include <cstdlib>
32 #include <cassert>
33 #include <iostream>
34 #include <sstream>
35
36 namespace TNT
37 {
38
39
40 template <class T>
41 class Matrix 
42 {
43
44
45   public:
46
47     typedef Subscript   size_type;
48     typedef         T   value_type;
49     typedef         T   element_type;
50     typedef         T*  pointer;
51     typedef         T*  iterator;
52     typedef         T&  reference;
53     typedef const   T*  const_iterator;
54     typedef const   T&  const_reference;
55
56     Subscript lbound() const { return 1;}
57  
58   protected:
59     Subscript m_;
60     Subscript n_;
61     Subscript mn_;      // total size
62     T* v_;                  
63     T** row_;           
64     T* vm1_ ;       // these point to the same data, but are 1-based 
65     T** rowm1_;
66
67     // internal helper function to create the array
68     // of row pointers
69
70     void initialize(Subscript M, Subscript N)
71     {
72         mn_ = M*N;
73         m_ = M;
74         n_ = N;
75
76         v_ = new T[mn_]; 
77         row_ = new T*[M];
78         rowm1_ = new T*[M];
79
80         assert(v_  != NULL);
81         assert(row_  != NULL);
82         assert(rowm1_ != NULL);
83
84         T* p = v_;              
85         vm1_ = v_ - 1;
86         for (Subscript i=0; i<M; i++)
87         {
88             row_[i] = p;
89             rowm1_[i] = p-1;
90             p += N ;
91             
92         }
93
94         rowm1_ -- ;     // compensate for 1-based offset
95     }
96    
97     void copy(const T*  v)
98     {
99         Subscript N = m_ * n_;
100         Subscript i;
101
102 #ifdef TNT_UNROLL_LOOPS
103         Subscript Nmod4 = N & 3;
104         Subscript N4 = N - Nmod4;
105
106         for (i=0; i<N4; i+=4)
107         {
108             v_[i] = v[i];
109             v_[i+1] = v[i+1];
110             v_[i+2] = v[i+2];
111             v_[i+3] = v[i+3];
112         }
113
114         for (i=N4; i< N; i++)
115             v_[i] = v[i];
116 #else
117
118         for (i=0; i< N; i++)
119             v_[i] = v[i];
120 #endif      
121     }
122
123     void set(const T& val)
124     {
125         Subscript N = m_ * n_;
126         Subscript i;
127
128 #ifdef TNT_UNROLL_LOOPS
129         Subscript Nmod4 = N & 3;
130         Subscript N4 = N - Nmod4;
131
132         for (i=0; i<N4; i+=4)
133         {
134             v_[i] = val;
135             v_[i+1] = val;
136             v_[i+2] = val;
137             v_[i+3] = val; 
138         }
139
140         for (i=N4; i< N; i++)
141             v_[i] = val;
142 #else
143
144         for (i=0; i< N; i++)
145             v_[i] = val;
146         
147 #endif      
148     }
149     
150
151     
152     void destroy()
153     {     
154         /* do nothing, if no memory has been previously allocated */
155         if (v_ == NULL) return ;
156
157         /* if we are here, then matrix was previously allocated */
158         if (v_ != NULL) delete [] (v_);     
159         if (row_ != NULL) delete [] (row_);
160
161         /* return rowm1_ back to original value */
162         rowm1_ ++;
163         if (rowm1_ != NULL ) delete [] (rowm1_);
164     }
165
166
167   public:
168
169     operator T**(){ return  row_; }
170     operator T**() const { return row_; }
171
172
173     Subscript size() const { return mn_; }
174
175     // constructors
176
177     Matrix() : m_(0), n_(0), mn_(0), v_(0), row_(0), vm1_(0), rowm1_(0) {};
178
179     Matrix(const Matrix<T> &A)
180     {
181         initialize(A.m_, A.n_);
182         copy(A.v_);
183     }
184
185     Matrix(Subscript M, Subscript N, const T& value = T())
186     {
187         initialize(M,N);
188         set(value);
189     }
190
191     Matrix(Subscript M, Subscript N, const T* v)
192     {
193         initialize(M,N);
194         copy(v);
195     }
196
197     Matrix(Subscript M, Subscript N, const char *s)
198     {
199         initialize(M,N);
200         //std::istrstream ins(s);
201         std::istringstream ins(s);
202
203         Subscript i, j;
204
205         for (i=0; i<M; i++)
206             for (j=0; j<N; j++)
207                 ins >> row_[i][j];
208     }
209
210     // destructor
211     //
212     ~Matrix()
213     {
214         destroy();
215     }
216
217
218     // reallocating
219     //
220     Matrix<T>& newsize(Subscript M, Subscript N)
221     {
222         if (num_rows() == M && num_cols() == N)
223             return *this;
224
225         destroy();
226         initialize(M,N);
227         
228         return *this;
229     }
230
231
232
233
234     // assignments
235     //
236     Matrix<T>& operator=(const Matrix<T> &A)
237     {
238         if (v_ == A.v_)
239             return *this;
240
241         if (m_ == A.m_  && n_ == A.n_)      // no need to re-alloc
242             copy(A.v_);
243
244         else
245         {
246             destroy();
247             initialize(A.m_, A.n_);
248             copy(A.v_);
249         }
250
251         return *this;
252     }
253         
254     Matrix<T>& operator=(const T& scalar)
255     { 
256         set(scalar); 
257         return *this;
258     }
259
260
261     Subscript dim(Subscript d) const 
262     {
263 #ifdef TNT_BOUNDS_CHECK
264         assert( d >= 1);
265         assert( d <= 2);
266 #endif
267         return (d==1) ? m_ : ((d==2) ? n_ : 0); 
268     }
269
270     Subscript num_rows() const { return m_; }
271     Subscript num_cols() const { return n_; }
272
273
274
275
276     inline T* operator[](Subscript i)
277     {
278 #ifdef TNT_BOUNDS_CHECK
279         assert(0<=i);
280         assert(i < m_) ;
281 #endif
282         return row_[i];
283     }
284
285     inline const T* operator[](Subscript i) const
286     {
287 #ifdef TNT_BOUNDS_CHECK
288         assert(0<=i);
289         assert(i < m_) ;
290 #endif
291         return row_[i];
292     }
293
294     inline reference operator()(Subscript i)
295     { 
296 #ifdef TNT_BOUNDS_CHECK
297         assert(1<=i);
298         assert(i <= mn_) ;
299 #endif
300         return vm1_[i]; 
301     }
302
303     inline const_reference operator()(Subscript i) const
304     { 
305 #ifdef TNT_BOUNDS_CHECK
306         assert(1<=i);
307         assert(i <= mn_) ;
308 #endif
309         return vm1_[i]; 
310     }
311
312
313
314     inline reference operator()(Subscript i, Subscript j)
315     { 
316 #ifdef TNT_BOUNDS_CHECK
317         assert(1<=i);
318         assert(i <= m_) ;
319         assert(1<=j);
320         assert(j <= n_);
321 #endif
322         return  rowm1_[i][j]; 
323     }
324
325
326     
327     inline const_reference operator() (Subscript i, Subscript j) const
328     {
329 #ifdef TNT_BOUNDS_CHECK
330         assert(1<=i);
331         assert(i <= m_) ;
332         assert(1<=j);
333         assert(j <= n_);
334 #endif
335         return rowm1_[i][j]; 
336     }
337
338
339
340
341 };
342
343
344 /* ***************************  I/O  ********************************/
345
346 template <class T>
347 std::ostream& operator<<(std::ostream &s, const Matrix<T> &A)
348 {
349     Subscript M=A.num_rows();
350     Subscript N=A.num_cols();
351
352     s << M << " " << N << "\n";
353
354     for (Subscript i=0; i<M; i++)
355     {
356         for (Subscript j=0; j<N; j++)
357         {
358             s << A[i][j] << " ";
359         }
360         s << "\n";
361     }
362
363
364     return s;
365 }
366
367 template <class T>
368 std::istream& operator>>(std::istream &s, Matrix<T> &A)
369 {
370
371     Subscript M, N;
372
373     s >> M >> N;
374
375     if ( !(M == A.num_rows() && N == A.num_cols() ))
376     {
377         A.newsize(M,N);
378     }
379
380
381     for (Subscript i=0; i<M; i++)
382         for (Subscript j=0; j<N; j++)
383         {
384             s >>  A[i][j];
385         }
386
387
388     return s;
389 }
390
391 // *******************[ basic matrix algorithms ]***************************
392
393
394 template <class T>
395 Matrix<T> operator+(const Matrix<T> &A, 
396     const Matrix<T> &B)
397 {
398     Subscript M = A.num_rows();
399     Subscript N = A.num_cols();
400
401     assert(M==B.num_rows());
402     assert(N==B.num_cols());
403
404     Matrix<T> tmp(M,N);
405     Subscript i,j;
406
407     for (i=0; i<M; i++)
408         for (j=0; j<N; j++)
409             tmp[i][j] = A[i][j] + B[i][j];
410
411     return tmp;
412 }
413
414 template <class T>
415 Matrix<T> operator-(const Matrix<T> &A, 
416     const Matrix<T> &B)
417 {
418     Subscript M = A.num_rows();
419     Subscript N = A.num_cols();
420
421     assert(M==B.num_rows());
422     assert(N==B.num_cols());
423
424     Matrix<T> tmp(M,N);
425     Subscript i,j;
426
427     for (i=0; i<M; i++)
428         for (j=0; j<N; j++)
429             tmp[i][j] = A[i][j] - B[i][j];
430
431     return tmp;
432 }
433
434 template <class T>
435 Matrix<T> mult_element(const Matrix<T> &A, 
436     const Matrix<T> &B)
437 {
438     Subscript M = A.num_rows();
439     Subscript N = A.num_cols();
440
441     assert(M==B.num_rows());
442     assert(N==B.num_cols());
443
444     Matrix<T> tmp(M,N);
445     Subscript i,j;
446
447     for (i=0; i<M; i++)
448         for (j=0; j<N; j++)
449             tmp[i][j] = A[i][j] * B[i][j];
450
451     return tmp;
452 }
453
454
455 template <class T>
456 Matrix<T> transpose(const Matrix<T> &A)
457 {
458     Subscript M = A.num_rows();
459     Subscript N = A.num_cols();
460
461     Matrix<T> S(N,M);
462     Subscript i, j;
463
464     for (i=0; i<M; i++)
465         for (j=0; j<N; j++)
466             S[j][i] = A[i][j];
467
468     return S;
469 }
470
471
472     
473 template <class T>
474 inline Matrix<T> matmult(const Matrix<T>  &A, 
475     const Matrix<T> &B)
476 {
477
478 #ifdef TNT_BOUNDS_CHECK
479     assert(A.num_cols() == B.num_rows());
480 #endif
481
482     Subscript M = A.num_rows();
483     Subscript N = A.num_cols();
484     Subscript K = B.num_cols();
485
486     Matrix<T> tmp(M,K);
487     T sum;
488
489     for (Subscript i=0; i<M; i++)
490     for (Subscript k=0; k<K; k++)
491     {
492         sum = 0;
493         for (Subscript j=0; j<N; j++)
494             sum = sum +  A[i][j] * B[j][k];
495
496         tmp[i][k] = sum; 
497     }
498
499     return tmp;
500 }
501
502 template <class T>
503 inline Matrix<T> operator*(const Matrix<T>  &A, 
504     const Matrix<T> &B)
505 {
506     return matmult(A,B);
507 }
508
509 template <class T>
510 inline int matmult(Matrix<T>& C, const Matrix<T>  &A, 
511     const Matrix<T> &B)
512 {
513
514     assert(A.num_cols() == B.num_rows());
515
516     Subscript M = A.num_rows();
517     Subscript N = A.num_cols();
518     Subscript K = B.num_cols();
519
520     C.newsize(M,K);
521
522     T sum;
523
524     const T* row_i;
525     const T* col_k;
526
527     for (Subscript i=0; i<M; i++)
528     for (Subscript k=0; k<K; k++)
529     {
530         row_i  = &(A[i][0]);
531         col_k  = &(B[0][k]);
532         sum = 0;
533         for (Subscript j=0; j<N; j++)
534         {
535             sum  += *row_i * *col_k;
536             row_i++;
537             col_k += K;
538         }
539         C[i][k] = sum; 
540     }
541
542     return 0;
543 }
544
545
546 template <class T>
547 Vector<T> matmult(const Matrix<T>  &A, const Vector<T> &x)
548 {
549
550 #ifdef TNT_BOUNDS_CHECK
551     assert(A.num_cols() == x.dim());
552 #endif
553
554     Subscript M = A.num_rows();
555     Subscript N = A.num_cols();
556
557     Vector<T> tmp(M);
558     T sum;
559
560     for (Subscript i=0; i<M; i++)
561     {
562         sum = 0;
563         const T* rowi = A[i];
564         for (Subscript j=0; j<N; j++)
565             sum = sum +  rowi[j] * x[j];
566
567         tmp[i] = sum; 
568     }
569
570     return tmp;
571 }
572
573 template <class T>
574 inline Vector<T> operator*(const Matrix<T>  &A, const Vector<T> &x)
575 {
576     return matmult(A,x);
577 }
578
579 } // namespace TNT
580
581 #endif
582 // CMAT_H