Fixed a stupid bug when exporting meshes with empty material slots.
[blender.git] / intern / iksolver / intern / TNT / vec.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 // Basic TNT  numerical vector (0-based [i] AND 1-based (i) indexing )
31 //
32
33 #ifndef VEC_H
34 #define VEC_H
35
36 #include "subscript.h"
37 #include <stdlib.h>
38 #include <assert.h>
39 #include <iostream>
40
41 namespace TNT
42 {
43
44 template <class T>
45 class Vector 
46 {
47
48
49   public:
50
51     typedef Subscript   size_type;
52     typedef         T   value_type;
53     typedef         T   element_type;
54     typedef         T*  pointer;
55     typedef         T*  iterator;
56     typedef         T&  reference;
57     typedef const   T*  const_iterator;
58     typedef const   T&  const_reference;
59
60     Subscript lbound() const { return 1;}
61  
62   protected:
63     T* v_;                  
64     T* vm1_;        // pointer adjustment for optimzied 1-offset indexing
65     Subscript n_;
66
67     // internal helper function to create the array
68     // of row pointers
69
70     void initialize(Subscript N)
71     {
72         // adjust pointers so that they are 1-offset:
73         // v_[] is the internal contiguous array, it is still 0-offset
74         //
75         assert(v_ == NULL);
76         v_ = new T[N];
77         assert(v_  != NULL);
78         vm1_ = v_-1;
79         n_ = N;
80     }
81    
82     void copy(const T*  v)
83     {
84         Subscript N = n_;
85         Subscript i;
86
87 #ifdef TNT_UNROLL_LOOPS
88         Subscript Nmod4 = N & 3;
89         Subscript N4 = N - Nmod4;
90
91         for (i=0; i<N4; i+=4)
92         {
93             v_[i] = v[i];
94             v_[i+1] = v[i+1];
95             v_[i+2] = v[i+2];
96             v_[i+3] = v[i+3];
97         }
98
99         for (i=N4; i< N; i++)
100             v_[i] = v[i];
101 #else
102
103         for (i=0; i< N; i++)
104             v_[i] = v[i];
105 #endif      
106     }
107
108     void set(const T& val)
109     {
110         Subscript N = n_;
111         Subscript i;
112
113 #ifdef TNT_UNROLL_LOOPS
114         Subscript Nmod4 = N & 3;
115         Subscript N4 = N - Nmod4;
116
117         for (i=0; i<N4; i+=4)
118         {
119             v_[i] = val;
120             v_[i+1] = val;
121             v_[i+2] = val;
122             v_[i+3] = val; 
123         }
124
125         for (i=N4; i< N; i++)
126             v_[i] = val;
127 #else
128
129         for (i=0; i< N; i++)
130             v_[i] = val;
131         
132 #endif      
133     }
134     
135
136
137     void destroy()
138     {     
139         /* do nothing, if no memory has been previously allocated */
140         if (v_ == NULL) return ;
141
142         /* if we are here, then matrix was previously allocated */
143         delete [] (v_);     
144
145         v_ = NULL;
146         vm1_ = NULL;
147     }
148
149
150   public:
151
152     // access
153
154     iterator begin() { return v_;}
155     iterator end()   { return v_ + n_; }
156     const iterator begin() const { return v_;}
157     const iterator end() const  { return v_ + n_; }
158
159     // destructor
160
161     ~Vector() 
162     {
163         destroy();
164     }
165
166     // constructors
167
168     Vector() : v_(0), vm1_(0), n_(0)  {};
169
170     Vector(const Vector<T> &A) : v_(0), vm1_(0), n_(0)
171     {
172         initialize(A.n_);
173         copy(A.v_);
174     }
175
176     Vector(Subscript N, const T& value = T()) :  v_(0), vm1_(0), n_(0)
177     {
178         initialize(N);
179         set(value);
180     }
181
182     Vector(Subscript N, const T* v) :  v_(0), vm1_(0), n_(0)
183     {
184         initialize(N);
185         copy(v);
186     }
187
188
189     // methods
190     // 
191     Vector<T>& newsize(Subscript N)
192     {
193         if (n_ == N) return *this;
194
195         destroy();
196         initialize(N);
197
198         return *this;
199     }
200
201
202     // assignments
203     //
204     Vector<T>& operator=(const Vector<T> &A)
205     {
206         if (v_ == A.v_)
207             return *this;
208
209         if (n_ == A.n_)         // no need to re-alloc
210             copy(A.v_);
211
212         else
213         {
214             destroy();
215             initialize(A.n_);
216             copy(A.v_);
217         }
218
219         return *this;
220     }
221         
222     Vector<T>& operator=(const T& scalar)
223     { 
224         set(scalar);  
225         return *this;
226     }
227
228     inline Subscript dim() const 
229     {
230         return  n_; 
231     }
232
233     inline Subscript size() const 
234     {
235         return  n_; 
236     }
237
238
239     inline reference operator()(Subscript i)
240     { 
241 #ifdef TNT_BOUNDS_CHECK
242         assert(1<=i);
243         assert(i <= n_) ;
244 #endif
245         return vm1_[i]; 
246     }
247
248     inline const_reference operator() (Subscript i) const
249     {
250 #ifdef TNT_BOUNDS_CHECK
251         assert(1<=i);
252         assert(i <= n_) ;
253 #endif
254         return vm1_[i]; 
255     }
256
257     inline reference operator[](Subscript i)
258     { 
259 #ifdef TNT_BOUNDS_CHECK
260         assert(0<=i);
261         assert(i < n_) ;
262 #endif
263         return v_[i]; 
264     }
265
266     inline const_reference operator[](Subscript i) const
267     {
268 #ifdef TNT_BOUNDS_CHECK
269         assert(0<=i) ;
270         assert(i < n_) ;
271 #endif
272         return v_[i]; 
273     }
274
275
276
277 };
278
279
280 /* ***************************  I/O  ********************************/
281
282 template <class T>
283 std::ostream& operator<<(std::ostream &s, const Vector<T> &A)
284 {
285     Subscript N=A.dim();
286
287     s <<  N << std::endl;
288
289     for (Subscript i=0; i<N; i++)
290         s   << A[i] << " " << std::endl;
291     s << std::endl;
292
293     return s;
294 }
295
296 template <class T>
297 std::istream & operator>>(std::istream &s, Vector<T> &A)
298 {
299
300     Subscript N;
301
302     s >> N;
303
304     if ( !(N == A.size() ))
305     {
306         A.newsize(N);
307     }
308
309
310     for (Subscript i=0; i<N; i++)
311             s >>  A[i];
312
313
314     return s;
315 }
316
317 // *******************[ basic matrix algorithms ]***************************
318
319
320 template <class T>
321 Vector<T> operator+(const Vector<T> &A, 
322     const Vector<T> &B)
323 {
324     Subscript N = A.dim();
325
326     assert(N==B.dim());
327
328     Vector<T> tmp(N);
329     Subscript i;
330
331     for (i=0; i<N; i++)
332             tmp[i] = A[i] + B[i];
333
334     return tmp;
335 }
336
337 template <class T>
338 Vector<T> operator-(const Vector<T> &A, 
339     const Vector<T> &B)
340 {
341     Subscript N = A.dim();
342
343     assert(N==B.dim());
344
345     Vector<T> tmp(N);
346     Subscript i;
347
348     for (i=0; i<N; i++)
349             tmp[i] = A[i] - B[i];
350
351     return tmp;
352 }
353
354 template <class T>
355 Vector<T> operator*(const Vector<T> &A, 
356     const Vector<T> &B)
357 {
358     Subscript N = A.dim();
359
360     assert(N==B.dim());
361
362     Vector<T> tmp(N);
363     Subscript i;
364
365     for (i=0; i<N; i++)
366             tmp[i] = A[i] * B[i];
367
368     return tmp;
369 }
370
371 template <class T>
372 Vector<T> operator*(const Vector<T> &A, 
373     const T &B)
374 {
375     Subscript N = A.dim();
376
377     Vector<T> tmp(N);
378     Subscript i;
379
380     for (i=0; i<N; i++)
381             tmp[i] = A[i] * B;
382
383     return tmp;
384 }
385
386
387 template <class T>
388 T dot_prod(const Vector<T> &A, const Vector<T> &B)
389 {
390     Subscript N = A.dim();
391     assert(N == B.dim());
392
393     Subscript i;
394     T sum = 0;
395
396     for (i=0; i<N; i++)
397         sum += A[i] * B[i];
398
399     return sum;
400 }
401
402 // inplace versions of the above template functions
403
404 // A = A + B
405
406 template <class T>
407 void vectoradd(
408         Vector<T> &A, 
409     const Vector<T> &B)
410 {
411     const Subscript N = A.dim();
412     assert(N==B.dim());
413     Subscript i;
414
415     for (i=0; i<N; i++)
416             A[i] += B[i];
417 }
418
419 // same with seperate output vector
420
421 template <class T>
422 void vectoradd(
423         Vector<T> &C,
424         const Vector<T> &A, 
425     const Vector<T> &B)
426 {
427     const Subscript N = A.dim();
428     assert(N==B.dim());
429     Subscript i;
430
431     for (i=0; i<N; i++)
432             C[i] = A[i] + B[i];
433 }
434
435 // A = A - B
436
437 template <class T>
438 void vectorsub(
439         Vector<T> &A, 
440     const Vector<T> &B)
441 {
442     const Subscript N = A.dim();
443     assert(N==B.dim());
444     Subscript i;
445
446     for (i=0; i<N; i++)
447             A[i] -= B[i];
448 }
449
450 template <class T>
451 void vectorsub(
452         Vector<T> &C,
453         const Vector<T> &A, 
454     const Vector<T> &B)
455 {
456     const Subscript N = A.dim();
457     assert(N==B.dim());
458     Subscript i;
459
460     for (i=0; i<N; i++)
461             C[i] = A[i] - B[i];
462 }
463
464 template <class T>
465 void vectorscale(
466         Vector<T> &C,
467         const Vector<T> &A, 
468     const T &B)
469 {
470     const Subscript N = A.dim();
471     Subscript i;
472
473     for (i=0; i<N; i++)
474             C[i] = A[i] * B;
475 }
476
477 template <class T>
478 void vectorscale(
479         Vector<T> &A, 
480     const T &B)
481 {
482     const Subscript N = A.dim();
483     Subscript i;
484
485     for (i=0; i<N; i++)
486             A[i] *= B;
487 }
488
489 }   /* namespace TNT */
490
491 #endif // VEC_H
492