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