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