HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
UT_Vector.h
Go to the documentation of this file.
1 /*
2  * PROPRIETARY INFORMATION. This software is proprietary to
3  * Side Effects Software Inc., and is not to be reproduced,
4  * transmitted, or disclosed in any way without written permission.
5  *
6  * NAME: Vector of arbitrary size (C++)
7  *
8  * COMMENTS: From Numerical Recipes
9  * And contains a permutation vector class.
10  *
11  */
12 
13 #ifndef __UT_Vector_H__
14 #define __UT_Vector_H__
15 
16 #include "UT_API.h"
17 #include <stdlib.h>
18 #include <iosfwd>
19 #include <SYS/SYS_Types.h>
20 #include "UT_Assert.h"
21 #include "UT_ThreadedAlgorithm.h"
22 
23 #include "UT_VectorTypes.h"
24 
25 //////////////////////////////////////////////////////////////////////////////
26 //
27 // UT_VectorT
28 //
29 
30 template <typename T>
31 class UT_VectorT
32 {
33 public:
34 
35  typedef T value_type;
36 
37  /// Input the index range [nl..nh].
38  UT_VectorT() { myVector = 0; myOwnData = 0; }
39  UT_VectorT(exint nl, exint nh);
40  UT_VectorT(const UT_VectorT<T> &v);
41  ~UT_VectorT();
42 
43  /// Initialize nl, nh and allocate space.
44  void init(exint nl, exint nh);
45 
46  /// Steal from another vector, resulting vector has origin at 1.
47  void subvector(const UT_VectorT<T> &v, exint nl, exint nh);
48 
49  /// Assign from a float array into the designated subvector, performing a
50  /// deep copy.
51  void assign(const fpreal32 *data, exint nl, exint nh);
52  void assign(const fpreal64 *data, exint nl, exint nh);
53 
54  int isInit() const { return myVector ? 1 : 0; }
55 
56  /// Initialize to zeros.
57  void zero()
58  { zero(myNL, myNH); }
59  THREADED_METHOD2(UT_VectorT, nh - nl > 4999, zero,
60  exint, nl, exint, nh);
61  void zeroPartial(exint nl, exint nh, const UT_JobInfo &info);
62  /// Initialize to the given constant.
63  void constant(T c)
64  { constant(myNL, myNH, c); }
65  void constant(exint nl, exint nh, T c);
66 
67  /// These methods allow one to read out and write into subvectors
68  /// of UT_Vector using our other vector classes.
69  /// Keep in mind that the usual UT_Vector? methods are 0 based,
70  /// while this class is usually 1 based.
71  void getSubvector2(UT_Vector2 &v, exint idx) const;
72  void setSubvector2(exint idx, const UT_Vector2 &v);
73  void getSubvector3(UT_Vector3 &v, exint idx) const;
74  void setSubvector3(exint idx, const UT_Vector3 &v);
75  void getSubvector4(UT_Vector4 &v, exint idx) const;
76  void setSubvector4(exint idx, const UT_Vector4 &v);
77 
78  /// Get the low index.
79  exint getNL() const { return myNL; }
80 
81  /// Get the high index.
82  exint getNH() const { return myNH; }
83 
84  /// Get dimension of this vector
85  exint length() const { return myNH - myNL + 1; }
86 
87  /// Determines if we have a large enough vector to make micro
88  /// threading worthwhile.
89  /// Experimentation shows with 4 procs on NT, 5000 is a cut off for good
90  /// behaviour on the simplistic addScaledVec method.
91  bool shouldMultiThread() const
92  {
93  return (myNH - myNL) > 4999;
94  }
95 
96  /// Determines the [start,end) interval that should be worked on
97  void getPartialRange(exint &start, exint &end, const UT_JobInfo &info) const;
98 
99  /// Determines the [startblock,endblock) interval that should be worked on,
100  /// where block 0 starts at myNL and all blocks are blocksize long, except
101  /// possibly the last.
102  void getPartialBlockRange(exint &startblock, exint &endblock,
103  const exint blocksize, const UT_JobInfo &info) const;
104 
105  /// Change the low index, and the high index will adjust itself.
106  void changeNL(exint nl);
107 
108  /// Change the indices in a very shallow manner, i.e. simply change
109  /// the index without shifting the data. This is highly dangerous, yet
110  /// very useful if used with care.
111  void setShallowNL(exint nl) { myNL = nl; }
112  void setShallowNH(exint nh) { myNH = nh; }
113 
114  /// For retrieving data
116  {
117  UT_ASSERT_P(i >= myNL && i <= myNH);
118  return myVector[i];
119  }
120  T operator()(exint i) const
121  {
122  UT_ASSERT_P(i >= myNL && i <= myNH);
123  return myVector[i];
124  }
125 
126  /// Lp-norm
127  /// type: 0 L-infinity norm (ie. max abs)
128  /// 1 L1-norm (ie. sum of abs)
129  /// 2 L2-norm (ie. Euclidean distance)
130  T norm(int type=2) const;
131 
132  /// Square of L2-norm
133  T norm2() const;
134 
135  T distance2(const UT_VectorT<T> &v) const;
136 
137  /// Negate
141  const UT_VectorT<T> &, v)
142  void negPlusPartial(const UT_VectorT<T> &v, const UT_JobInfo &info);
143 
144  /// Add scaled vector
145  THREADED_METHOD2(UT_VectorT, shouldMultiThread(), addScaledVec,
146  T, s,
147  const UT_VectorT<T> &, v)
148  void addScaledVecPartial(T s, const UT_VectorT<T> &v,
149  const UT_JobInfo &info);
150 
151  /// Add scaled vector and compute the squared L2 norm
152  void addScaledVecNorm2(T s, const UT_VectorT<T> &v, fpreal64 *norm2);
153 
154  /// Add scaled vector and compute the squared L2 norm,
155  /// but only using components up to (but excluding) index normlimit
156  /// for the norm computation.
157  void addScaledVecNorm2UpTo(T s, const UT_VectorT<T> &v, fpreal64 *norm2, exint normlimit);
158 
159  /// Scale itself then add vector
160  THREADED_METHOD2(UT_VectorT<T>, shouldMultiThread(), scaleAddVec,
161  T, s,
162  const UT_VectorT<T> &, v)
163  void scaleAddVecPartial(T s, const UT_VectorT<T> &v,
164  const UT_JobInfo &info);
165 
166  /// Multiply two sources together, save to this. Requires
167  /// that three vectors already have the same size.
168  THREADED_METHOD2(UT_VectorT<T>, shouldMultiThread(), multAndSet,
169  const UT_VectorT<T> &, a,
170  const UT_VectorT<T> &, b)
171  void multAndSetPartial(const UT_VectorT<T> &a,
172  const UT_VectorT<T> &b,
173  const UT_JobInfo &info);
174 
175  /// Multiply two sources together, save to this. Requires
176  /// that three vectors already have the same size.
177  /// Computes dot(a*b, a)
178  /// This strange seeming pair is useful for jacobian preconditioners.
179  void multSetAndDotUpTo(const UT_VectorT<T> &a, const UT_VectorT<T> &b,
181 
182 
183  /// Divides two sources together, save to this. Requires
184  /// that three vectors already have the same size.
185  THREADED_METHOD2(UT_VectorT<T>, shouldMultiThread(), divAndSet,
186  const UT_VectorT<T> &, a,
187  const UT_VectorT<T> &, b)
188  void divAndSetPartial(const UT_VectorT<T> &a,
189  const UT_VectorT<T> &b,
190  const UT_JobInfo &info);
191 
192  /// Inverts the source and saves to this.
193  /// Requires that the vectors match.
194  THREADED_METHOD1(UT_VectorT<T>, shouldMultiThread(), safeInvertAndSet,
195  const UT_VectorT<T> &, a)
196  void safeInvertAndSetPartial(const UT_VectorT<T> &a,
197  const UT_JobInfo &info);
198 
199  THREADED_METHOD1(UT_VectorT<T>, shouldMultiThread(), invertAndSet,
200  const UT_VectorT<T> &, a)
201  void invertAndSetPartial(const UT_VectorT<T> &a,
202  const UT_JobInfo &info);
203 
204  T dot(const UT_VectorT<T> &v) const;
205 
206  /// Multithreaded copy. operator= uses this.
207  THREADED_METHOD1(UT_VectorT<T>, shouldMultiThread(), copyFrom,
208  const UT_VectorT<T> &, v)
209  void copyFromPartial(const UT_VectorT<T> &v,
210  const UT_JobInfo &info);
211 
212  /// Operators
213  /// NOTE: operator= requires the destination be a matching size and
214  /// layout!
215  UT_VectorT &operator= (const UT_VectorT<T> &v);
216  UT_VectorT &operator+= (const UT_VectorT<T> &v);
217  UT_VectorT &operator-= (const UT_VectorT<T> &v);
218 
219  /// Componentwise multiplication & division.
220  UT_VectorT &operator*= (const UT_VectorT<T> &v);
221  UT_VectorT &operator/= (const UT_VectorT<T> &v);
222 
223  /// Scalar multiplication and division.
224  UT_VectorT &operator*= (T scalar);
225  UT_VectorT &operator/= (T scalar);
226 
227  /// Equals
228  /// upls is the number of representable fpreal64's to accept between
229  /// each individual component. If you deal in fpreal32's and assign into
230  /// this fpreal64 class, then you will want a ulps that is scaled by
231  /// 2^32. eg. for 50upls at fpreal32 precision, you will need 50*2^32
232  /// or approximately 1e+11.
233  bool isEqual(const UT_VectorT<T> &v, int64 ulps);
234 
235  /// Output
236  std::ostream &save(std::ostream &os) const;
237  friend std::ostream &operator<<(std::ostream &os, const UT_VectorT<T> &v)
238  { v.save(os); return os; }
239 
240  T *getData() const { return &myVector[myNL]; }
241 
242  bool hasNan() const;
243  void testForNan() const;
244 
245 protected:
246  /// A multithreaded norm function that returns the non-normalized
247  /// norm (Ie, squared in the L2 case.) and is mulithreaded.
249  fpreal64 *, result,
250  int, type)
251  void normInternalPartial(fpreal64 *result, int type, const UT_JobInfo &info) const;
252 
254  fpreal64 *, result,
255  const UT_VectorT<T> &, v)
256  void distance2InternalPartial(fpreal64 *result, const UT_VectorT<T> &v, const UT_JobInfo &info) const;
257 
258  /// A multithreaded dot method
259  THREADED_METHOD2_CONST(UT_VectorT, shouldMultiThread(), dotInternal,
260  fpreal64 *, result,
261  const UT_VectorT<T> &, v)
262  void dotInternalPartial(fpreal64 *result,
263  const UT_VectorT<T> &v,
264  const UT_JobInfo &info) const;
265 
266  THREADED_METHOD3(UT_VectorT, shouldMultiThread(), addScaledVecNorm2Internal,
267  T, s,
268  const UT_VectorT<T> &, v,
269  fpreal64 *, norm2)
270  void addScaledVecNorm2InternalPartial(T s, const UT_VectorT<T> &v,
271  fpreal64 *norm2, const UT_JobInfo &info);
272 
273  THREADED_METHOD4(UT_VectorT, shouldMultiThread(), addScaledVecNorm2UpToInternal,
274  T, s,
275  const UT_VectorT<T> &, v,
276  fpreal64 *, norm2,
277  exint, normlimit)
278  void addScaledVecNorm2UpToInternalPartial(T s, const UT_VectorT<T> &v,
279  fpreal64 *norm2, exint normlimit, const UT_JobInfo &info);
280 
281  THREADED_METHOD4(UT_VectorT<T>, shouldMultiThread(), multSetAndDotUpToInternal,
282  const UT_VectorT<T> &, a,
283  const UT_VectorT<T> &, b,
284  fpreal64 *, dot_aba,
285  exint, dotlimit)
286  void multSetAndDotUpToInternalPartial(
287  const UT_VectorT<T> &a, const UT_VectorT<T> &b,
288  fpreal64 *dot_aba, exint dotlimit,
289  const UT_JobInfo &info);
290 
291  THREADED_METHOD3(UT_VectorT, nh - nl > 4999, constantInternal,
292  exint, nl, exint, nh, T, c);
293  void constantInternalPartial(exint nl, exint nh, T c, const UT_JobInfo &info);
294 
295  static const exint PARALLEL_BLOCK_SIZE = 1024;
296 
297 private:
298  exint myNL, myNH;
299  int myOwnData;
300  T *myVector;
301 };
302 
303 template <>
304 inline void UT_VectorT<fpreal32>::constant(exint nl, exint nh, fpreal32 c);
305 template <>
306 inline void UT_VectorT<int32>::constant(exint nl, exint nh, int32 c);
307 
308 
309 //////////////////////////////////////////////////////////////////////////////
310 //
311 // UT_Permutation
312 //
313 
314 
315 template <typename T>
317 {
318 public:
319  // Input the index range [nl..nh].
320  UT_PermutationT(exint nl, exint nh);
321  ~UT_PermutationT();
322 
323  // Assignment
324  UT_PermutationT(const UT_PermutationT<T> &p);
326 
327  // Initialize to zeros.
328  void zero();
329 
330  // Get the low index.
331  exint getNL() const { return myNL; }
332 
333  // Get the high index.
334  exint getNH() const { return myNH; }
335 
336  exint length() const { return myNH - myNL + 1; }
337 
338  // Change the low index, and the high index will adjust itself.
339  void changeNL(exint nl);
340 
341  // Change the indices in a very shallow manner, i.e. simply change
342  // the index without shifting the data. This is highly dangerous, yet
343  // very useful if used with care.
344  void setShallowNL(exint nl) { myNL = nl; }
345  void setShallowNH(exint nh) { myNH = nh; }
346 
347  // For retrieve data.
348  T &operator()(exint i)
349  {
350  UT_ASSERT_P(i >= myNL && i <= myNH);
351  return myVector[i];
352  }
353  T operator()(exint i) const
354  {
355  UT_ASSERT_P(i >= myNL && i <= myNH);
356  return myVector[i];
357  }
358 
359 private:
360  void clone(const UT_PermutationT<T> &p);
361  exint myNL, myNH;
362  T *myVector;
363 };
364 
365 #include "UT_Vector.C"
366 
367 
368 //////////////////////////////////////////////////////////////////////////////
369 //
370 // UT_Vector typedefs
371 //
372 
378 
379 //////////////////////////////////////////////////////////////////////////////
380 //
381 // UT_Vector inline binary operators
382 //
383 
384 // Free floating functions
385 inline fpreal64 dot(const UT_VectorD &v1, const UT_VectorD &v2);
386 inline fpreal64 distance2(const UT_VectorD &v1, const UT_VectorD &v2);
387 inline fpreal32 dot(const UT_VectorF &v1, const UT_VectorF &v2);
388 inline fpreal32 distance2(const UT_VectorF &v1, const UT_VectorF &v2);
389 
390 /// Dot product
391 inline fpreal64
392 dot(const UT_VectorD &v1, const UT_VectorD &v2)
393 {
394  return v1.dot(v2);
395 }
396 
397 /// Distance squared (L2) aka quadrance
398 inline fpreal64
400 {
401  return v1.distance2(v2);
402 }
403 
404 inline fpreal32
405 dot(const UT_VectorF &v1, const UT_VectorF &v2)
406 {
407  return v1.dot(v2);
408 }
409 
410 /// Distance squared (L2) aka quadrance
411 inline fpreal32
413 {
414  return v1.distance2(v2);
415 }
416 
417 // Overload for custom formatting of UT_VectorT<T> with UTformat.
418 template<typename T> inline size_t
419 format(char *buffer, size_t buffer_size, const UT_VectorT<T> &v);
420 
421 #endif
422 
GLdouble s
Definition: glew.h:1390
bool hasNan() const
Definition: UT_Vector.C:975
void zero()
Initialize to zeros.
Definition: UT_Vector.h:57
int const UT_JobInfo &info const
Definition: UT_Vector.h:251
void getPartialRange(exint &start, exint &end, const UT_JobInfo &info) const
Determines the [start,end) interval that should be worked on.
Definition: UT_Vector.C:136
exint getNH() const
Definition: UT_Vector.h:334
void negPartial(const UT_JobInfo &info)
Negate.
Definition: UT_Vector.C:436
int int32
Definition: SYS_Types.h:39
const UT_VectorT< T > fpreal64 * dot_aba
Definition: UT_Vector.h:287
void setShallowNH(exint nh)
Definition: UT_Vector.h:112
int isInit() const
Definition: UT_Vector.h:54
UT_VectorT< fpreal64 > UT_VectorD
Definition: UT_Vector.h:376
*get result *(waiting if necessary)*A common idiom is to fire a bunch of sub tasks at the and then *wait for them to all complete We provide a helper class
Definition: thread.h:643
void init(exint nl, exint nh)
Initialize nl, nh and allocate space.
Definition: UT_Vector.C:107
UT_VectorT< fpreal64 > UT_Vector
Definition: UT_Vector.h:377
T & operator()(exint i)
For retrieving data.
Definition: UT_Vector.h:115
GLboolean GLboolean GLboolean GLboolean a
Definition: glew.h:9477
exint getNH() const
Get the high index.
Definition: UT_Vector.h:82
void setShallowNL(exint nl)
Definition: UT_Vector.h:111
int64 exint
Definition: SYS_Types.h:125
void getPartialBlockRange(exint &startblock, exint &endblock, const exint blocksize, const UT_JobInfo &info) const
Definition: UT_Vector.C:146
void addScaledVecNorm2(T s, const UT_VectorT< T > &v, fpreal64 *norm2)
Add scaled vector and compute the squared L2 norm.
Definition: UT_Vector.C:503
THREADED_METHOD2_CONST(UT_VectorT, shouldMultiThread(), normInternal, fpreal64 *, result, int, type) void normInternalPartial(fpreal64 *result
T * getData() const
Definition: UT_Vector.h:240
void constantInternalPartial(exint nl, exint nh, T c, const UT_JobInfo &info)
Definition: UT_Vector.C:168
const GLdouble * v
Definition: glew.h:1391
UT_VectorT< fpreal32 > UT_VectorF
Definition: UT_Vector.h:375
UT_PermutationT< int > UT_Permutation
Definition: UT_Vector.h:373
exint length() const
Definition: UT_Vector.h:336
T dot(const UT_VectorT< T > &v) const
Definition: UT_Vector.C:906
void multSetAndDotUpTo(const UT_VectorT< T > &a, const UT_VectorT< T > &b, fpreal64 *dot_aba, exint dotlimit)
Definition: UT_Vector.C:674
float fpreal32
Definition: SYS_Types.h:200
void testForNan() const
Definition: UT_Vector.C:987
void addScaledVecNorm2UpTo(T s, const UT_VectorT< T > &v, fpreal64 *norm2, exint normlimit)
Definition: UT_Vector.C:562
THREADED_METHOD4(UT_VectorT, shouldMultiThread(), addScaledVecNorm2UpToInternal, T, s, const UT_VectorT< T > &, v, fpreal64 *, norm2, exint, normlimit) void addScaledVecNorm2UpToInternalPartial(T s
THREADED_METHOD3(UT_VectorT, shouldMultiThread(), addScaledVecNorm2Internal, T, s, const UT_VectorT< T > &, v, fpreal64 *, norm2) void addScaledVecNorm2InternalPartial(T s
double fpreal64
Definition: SYS_Types.h:201
T & operator()(exint i)
Definition: UT_Vector.h:348
T value_type
Definition: UT_Vector.h:35
~UT_VectorT()
Definition: UT_Vector.C:53
void changeNL(exint nl)
Change the low index, and the high index will adjust itself.
Definition: UT_Vector.C:267
THREADED_METHOD2(UT_VectorT, nh-nl > 4999, zero, exint, nl, exint, nh)
GLfloat GLfloat GLfloat v2
Definition: glew.h:1856
T norm(int type=2) const
Definition: UT_Vector.C:278
void zeroPartial(exint nl, exint nh, const UT_JobInfo &info)
Definition: UT_Vector.C:155
GLuint buffer
Definition: glew.h:1680
GLint GLenum GLsizei GLint GLsizei const void * data
Definition: glew.h:1379
#define UT_ASSERT_P(ZZ)
Definition: UT_Assert.h:134
UT_VectorT()
Input the index range [nl..nh].
Definition: UT_Vector.h:38
std::ostream & save(std::ostream &os) const
Output.
Definition: UT_Vector.C:949
void getSubvector4(UT_Vector4 &v, exint idx) const
Definition: UT_Vector.C:247
GLuint GLuint end
Definition: glew.h:1253
void
Definition: png.h:1083
T distance2(const UT_VectorT< T > &v) const
Definition: UT_Vector.C:390
const GLfloat * c
Definition: glew.h:16296
static const exint PARALLEL_BLOCK_SIZE
Definition: UT_Vector.h:295
T norm2() const
Square of L2-norm.
long long int64
Definition: SYS_Types.h:116
void setSubvector3(exint idx, const UT_Vector3 &v)
Definition: UT_Vector.C:238
GLuint GLuint GLsizei GLenum type
Definition: glew.h:1253
T operator()(exint i) const
Definition: UT_Vector.h:353
GLuint start
Definition: glew.h:1253
UT_VectorT< fpreal > UT_VectorR
Definition: UT_Vector.h:374
void getSubvector2(UT_Vector2 &v, exint idx) const
Definition: UT_Vector.C:213
size_t format(char *buffer, size_t buffer_size, const UT_VectorT< T > &v)
const UT_VectorT< T > fpreal64 exint dotlimit
Definition: UT_Vector.h:287
GLdouble GLdouble GLdouble b
Definition: glew.h:9122
GLfloat GLfloat p
Definition: glew.h:16321
void setShallowNL(exint nl)
Definition: UT_Vector.h:344
THREADED_METHOD1(UT_VectorT, shouldMultiThread(), negPlus, const UT_VectorT< T > &, v) void negPlusPartial(const UT_VectorT< T > &v
void setSubvector4(exint idx, const UT_Vector4 &v)
Definition: UT_Vector.C:257
void setSubvector2(exint idx, const UT_Vector2 &v)
Definition: UT_Vector.C:221
void getSubvector3(UT_Vector3 &v, exint idx) const
Definition: UT_Vector.C:229
const UT_VectorT< T > fpreal64 exint normlimit
Definition: UT_Vector.h:278
void subvector(const UT_VectorT< T > &v, exint nl, exint nh)
Steal from another vector, resulting vector has origin at 1.
Definition: UT_Vector.C:61
GLuint64EXT * result
Definition: glew.h:14007
exint getNL() const
Definition: UT_Vector.h:331
void assign(const fpreal32 *data, exint nl, exint nh)
Definition: UT_Vector.C:91
UT_VectorT & operator=(const UT_VectorT< T > &v)
Definition: UT_Vector.C:814
T operator()(exint i) const
Definition: UT_Vector.h:120
bool shouldMultiThread() const
Definition: UT_Vector.h:91
exint length() const
Get dimension of this vector.
Definition: UT_Vector.h:85
void constant(T c)
Initialize to the given constant.
Definition: UT_Vector.h:63
void setShallowNH(exint nh)
Definition: UT_Vector.h:345
#define THREADED_METHOD(CLASSNAME, DOMULTI, METHOD)
GLfloat GLfloat v1
Definition: glew.h:1852
const UT_JobInfo & info
Definition: UT_Vector.h:142
bool isEqual(const UT_VectorT< T > &v, int64 ulps)
Definition: UT_Vector.C:961
exint getNL() const
Get the low index.
Definition: UT_Vector.h:79