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