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