HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
UT_Matrix.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: Matrix of arbitrary size (C++)
7  *
8  * COMMENTS: From Numerical Recipes
9  * See Matrix Computations, Gene H. Golub, Charles F. Van Loan
10  * Rewritten without row array by Jeff Lait
11  *
12  */
13 
14 #ifndef __UT_Matrix_H__
15 #define __UT_Matrix_H__
16 
17 #include "UT_API.h"
18 #include "UT_Assert.h"
19 #include "UT_Matrix2.h"
20 #include "UT_Vector.h"
21 #include "UT_VectorTypes.h"
22 #include <SYS/SYS_Types.h>
23 #include <stdlib.h>
24 #include <iosfwd>
25 
26 
27 template <typename T>
29 {
30 public:
31 
32  typedef T value_type;
33 
34  UT_MatrixT();
35 
36  // Input the row and column index range [nrl..nrh][ncl..nch].
37  // The "standard" way to set up a matrix of R rows and C columns
38  // would be:
39  // UT_MatrixF foo(1, R, 1, C);
40 
41  UT_MatrixT(int nrl, int nrh, int ncl, int nch);
42  UT_MatrixT(int nrl, int nrh, int ncl, int nch, T *array);
43  UT_MatrixT(int nrl, int ncl, UT_MatrixT &mat);
44  explicit UT_MatrixT(const UT_MatrixT<T> &mat); // deep copy constructor
45  ~UT_MatrixT();
46 
47  // Allocate space. Resulting memory is not initialized.
48  void init(int nrl, int nrh, int ncl, int nch);
49 
50  // Enlarge/reduce matrix to have a new number of rows and
51  // columns. NRL & NCL are unchanged. Afterwards, the
52  // matrix always owns its data!
53  void resize(int nrows, int ncols);
54 
55  // Enlarge by one row/column and copy given data in (if given, otherwise
56  // filled in with zeroes
57  void appendCol(T *new_col=0);
58  void appendRow(T *new_row=0);
59 
60  // insert new row of zeros (or nw_row, if given) at the specified row idx
61  //
62  // this has been tested and it works, we just don't use it anywhere.
63  // uncomment it if you need it.
64  // void insertRowAt(int at, T *new_row=0);
65 
66  int isInit() const { return myMatrix ? 1 : 0; }
67 
68  // Steal from another matrix or array
69  // Resulting matrix has origin 1, 1.
70  void submatrix(const UT_MatrixT<T> &A,
71  int nrl, int nrh, int ncl, int nch);
72  // &(*this)(nrl, ncl) == array. If array is wider than
73  // nch - ncl + 1, stride should be set.
74  void submatrix(T *array,
75  int nrl, int nrh, int ncl, int nch, int stride=-1);
76 
77  // Sets to zeros. The version with a range will set a sub-block
78  // to zero.
79  void zero(int nrl, int nrh, int ncl, int nch);
80  void zero()
81  { zero(myNRL, myNRH, myNCL, myNCH); }
82 
83  /// These methods allow one to read out and write into blocks of
84  /// the UT_Matrix using our other matrix classes. This is very useful
85  /// when working with block algorithms.
86  /// Keep in mind that the usual UT_MatrixF? methods are 0 based while
87  /// this class is usually 1 based.
88  /// Only double matrices are supported as they are the preferred
89  /// matrix format.
90  void getSubmatrix2(UT_Matrix2T<T> &mat, int row, int col) const;
91  void setSubmatrix2(int row, int col, const UT_Matrix2T<T> &mat);
92  void addSubmatrix2(int row, int col, const UT_Matrix2T<T> &mat);
93 
94  void getSubmatrix3(UT_Matrix3T<T> &mat, int row, int col) const;
95  void setSubmatrix3(int row, int col, const UT_Matrix3T<T> &mat);
96  void addSubmatrix3(int row, int col, const UT_Matrix3T<T> &mat);
97 
98  void getSubmatrix4(UT_Matrix4T<T> &mat, int row, int col) const;
99  void setSubmatrix4(int row, int col, const UT_Matrix4T<T> &mat);
100  void addSubmatrix4(int row, int col, const UT_Matrix4T<T> &mat);
101 
102  // Make identity:
103  void makeIdentity();
104 
105  // Negate the matrix.
106  void negate();
107 
108  // Get the row low index.
109  int getNRL() const { return myNRL; }
110 
111  // Get the row high index.
112  int getNRH() const { return myNRH; }
113 
114  // Get the column low index.
115  int getNCL() const { return myNCL; }
116 
117  // Get the column high index.
118  int getNCH() const { return myNCH; }
119 
120  // Get number of rows
121  int rows() const { return myNRH - myNRL + 1; }
122 
123  // Get number of columns
124  int columns() const { return myNCH - myNCL + 1; }
125 
126  // Change the row and column low indices,
127  // and high indices will adjust themselves.
128  void changeNRLAndNCL(int nrl, int ncl);
129 
130  // Change the indices in a very shallow manner, i.e. simply change
131  // the index without shifting the data. This is highly dangerous, yet
132  // very useful if used with care.
133  void setShallowNRL(int nrl) { myNRL = nrl; }
134  void setShallowNRH(int nrh) { myNRH = nrh; }
135  void setShallowNCL(int ncl) { myNCL = ncl; }
136  void setShallowNCH(int nch) { myNCH = nch; }
137 
138  // Retrieve data
139  T &operator()(int row, int col)
140  {
141  UT_ASSERT_P(row >= myNRL && row <= myNRH);
142  UT_ASSERT_P(col >= myNCL && col <= myNCH);
143  return myMatrix[row*myStride + col];
144  }
145  T operator()(int row, int col) const
146  {
147  UT_ASSERT_P(row >= myNRL && row <= myNRH);
148  UT_ASSERT_P(col >= myNCL && col <= myNCH);
149  return myMatrix[row*myStride + col];
150  }
151 
152  UT_MatrixT<T> &operator=(const UT_MatrixT<T> &m);
153 
154  // Component wise operators
155  UT_MatrixT<T> &operator+=(const UT_MatrixT<T> &m);
156  UT_MatrixT<T> &operator-=(const UT_MatrixT<T> &m);
157 
158  // left multiply this matrix[1..m][1..n] by the given matrix
159  // A[1..p][1..m] and puts the product in 'result[1..p][1..n]'
160  // NOTE: this, A, and result must all have the same lower index
161  // for rows and columns
162  void preMult(const UT_MatrixT<T> &A, UT_MatrixT<T> &result) const;
163 
164  // right multiply this matrix[1..m][1..n] by the given matrix
165  // A[1..n][1..p] and puts the product in 'result[1..m][1..p]'
166  // NOTE: this, A, and result must all have the same lower index
167  // for rows and columns
168  void postMult(const UT_MatrixT<T> &A, UT_MatrixT<T> &result) const;
169 
170  // Vector multiplication, before and after.
171  // result = xT * this
172  template <typename S>
173  void preMult(const UT_VectorT<S> &x, UT_VectorT<S> &result) const;
174  // result = this * x
175  template <typename S>
176  void postMult(const UT_VectorT<S> &x, UT_VectorT<S> &result) const
177  { multVec(x, result); }
178 
179  // Fills the upper-triangular portion of this matrix with A^T . A
180  void upperNormalUpdate(const UT_MatrixT<T> &A);
181 
182  // Scaled outer product update.
183  // this += b * x * yT
184  template <typename S>
185  void outerproductUpdate(T b,
186  const UT_VectorT<S> &x, const UT_VectorT<S> &y);
187 
188  // Add another matrix scaled appropriately:
189  void addScaledMatrix(const UT_MatrixT<T> &A, T scale);
190  // this = A*scale
191  void setAndScale(const UT_MatrixT<T> &A, T scale);
192 
193  // Givens functions: do a pre or post multiplication with
194  // the givens matrix.
195  void preMultGivensInPlace(T c, T s);
196  void postMultGivensInPlace(T c, T s);
197 
198  template <typename S>
199  void multVec(const UT_VectorT<S> &x, UT_VectorT<S> &result) const;
200 
201  // Calculate the L2-distance (Euclidean) between two rows for columns [cl..ch].
202  T rowsL2dist(int r1, int r2, int cl=-1, int ch=-1) const;
203 
204  // Norm calculations.
205  T normFrobenius() const;
206  T norm1() const;
207  T normInfinite() const;
208 
209  // place this matrix's transpose in the given result matrix
210  void transpose( UT_MatrixT<T> &result ) const;
211 
212  // Returns true if matrix is symmetric within given tolerance
213  bool isSymmetric(T tolerance = SYS_FTOLERANCE_D) const;
214 
215  // Get the pointer to a specific row.
216  // Note real data starts at row[myNCL]
217  T *row(int i) const
218  {
219  UT_ASSERT_P(i >= myNRL && i <= myNRH);
220  return &myMatrix[i * myStride];
221  }
222 
223  void clearAndDestroy();
224 
225  int save(std::ostream &os, int binary) const;
226 
227  // I/O friends
228  friend std::ostream &operator<<(std::ostream &os, const UT_MatrixT<T> &m )
229  {
230  m.outTo(os);
231  return os;
232  }
233 
234  int64 getMemoryUsage(bool inclusive) const
235  {
236  int64 mem = inclusive ? sizeof(*this) : 0;
237  if (myOwnData && myMatrix)
238  mem += (myNRH - myNRL + 1)*(myNCH - myNCL + 1)*sizeof(*myMatrix);
239  return mem;
240  }
241 
242 protected:
243  // I/O methods
244  void outTo( std::ostream &os ) const;
245 
246 private:
247  int myNRL, myNRH, myNCL, myNCH;
248  T *myMatrix;
249  int myStride, myOffset;
250  char myOwnData;
251 };
252 
256 typedef UT_MatrixT<fpreal64> UT_Matrix; // Match precision w/UT_Vector
257 
258 #endif
int getNRL() const
Definition: UT_Matrix.h:109
T & operator()(int row, int col)
Definition: UT_Matrix.h:139
png_bytep png_bytep new_row
Definition: png.h:2130
void setShallowNCL(int ncl)
Definition: UT_Matrix.h:135
const GLuint GLenum const void * binary
Definition: glcorearb.h:1923
void postMult(const UT_VectorT< S > &x, UT_VectorT< S > &result) const
Definition: UT_Matrix.h:176
void setShallowNRL(int nrl)
Definition: UT_Matrix.h:133
#define UT_API
Definition: UT_API.h:12
GLint y
Definition: glcorearb.h:102
#define SYS_FTOLERANCE_D
Definition: SYS_Types.h:197
UT_MatrixT< fpreal > UT_MatrixR
Definition: UT_Matrix.h:253
png_uint_32 i
Definition: png.h:2877
#define UT_ASSERT_P(ZZ)
Definition: UT_Assert.h:101
GLsizei GLboolean transpose
Definition: glcorearb.h:831
long long int64
Definition: SYS_Types.h:100
GA_API const UT_StringHolder scale
T * row(int i) const
Definition: UT_Matrix.h:217
int getNCH() const
Definition: UT_Matrix.h:118
GLint GLenum GLboolean GLsizei stride
Definition: glcorearb.h:871
int columns() const
Definition: UT_Matrix.h:124
void setShallowNCH(int nch)
Definition: UT_Matrix.h:136
png_bytepp row
Definition: png.h:1836
int nrows
Definition: png.h:1821
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1221
void zero()
Definition: UT_Matrix.h:80
int rows() const
Definition: UT_Matrix.h:121
UT_MatrixT< fpreal64 > UT_MatrixD
Definition: UT_Matrix.h:255
bool isSymmetric(const MatType &m)
Determine if a matrix is symmetric.
Definition: Mat.h:875
GLint GLenum GLint x
Definition: glcorearb.h:408
int getNCL() const
Definition: UT_Matrix.h:115
UT_MatrixT< fpreal32 > UT_MatrixF
Definition: UT_Matrix.h:254
UT_MatrixT< fpreal64 > UT_Matrix
Definition: UT_Matrix.h:256
T operator()(int row, int col) const
Definition: UT_Matrix.h:145
int isInit() const
Definition: UT_Matrix.h:66
int getNRH() const
Definition: UT_Matrix.h:112
void setShallowNRH(int nrh)
Definition: UT_Matrix.h:134
int64 getMemoryUsage(bool inclusive) const
Definition: UT_Matrix.h:234
T value_type
Definition: UT_Matrix.h:32