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