HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
UT_SparseMatrix.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: Sparse Matrix Class (C++)
7  *
8  * COMMENTS:
9  *
10  * Storage format: Co-ordinate (COO). See
11  * Y. Saad. SPARSKIT: A basic toolkit for sparse matrix computations.
12  * Research Institute for Advanced Computer Science, RIACS-90-20,
13  * 1990.
14  * for details.
15  */
16 
17 #ifndef __UT_SparseMatrix_H__
18 #define __UT_SparseMatrix_H__
19 
20 #include "UT_API.h"
21 
22 #include "UT_IntArray.h"
23 #include "UT_ThreadedAlgorithm.h"
24 #include "UT_Vector.h"
25 #include "UT_Interrupt.h"
26 #include "UT_ParallelUtil.h"
27 
28 #include "UT_VectorTypes.h"
29 
30 #include <iterator>
31 #include <type_traits>
32 
33 class UT_IStream;
34 class UT_BitArray;
35 
36 template <typename T> class UT_API UT_SparseMatrixRowT;
37 
38 template <typename T, bool IsPaged>
40 {
41  class ut_MatrixCell
42  {
43  public:
44  int myRow;
45  int myCol;
46  T myValue;
47 
48  inline bool operator<(const ut_MatrixCell &o) const
49  {
50  if (myRow < o.myRow)
51  return true;
52  if (myRow == o.myRow && myCol < o.myCol)
53  return true;
54  return false;
55  }
56  };
57 
58  // For paged matrices, store the size of each page
59  // We always allocate full pages.
60  static const int CELL_PAGESIZE = 1024;
61  static const int CELL_PAGEMASK = 1023;
62  static const int CELL_PAGEBITS = 10;
63 
64  class ut_CellIterator : public std::iterator<std::random_access_iterator_tag, ut_MatrixCell>
65  {
66  public:
67  ut_CellIterator(const UT_SparseMatrixT<T, IsPaged> *matrix, int idx)
68  { myMatrix = matrix; myPos = idx; repage(); }
69 
70  ut_CellIterator operator+(ptrdiff_t n) { return ut_CellIterator(myMatrix, myPos+n); }
71  ut_CellIterator operator-(ptrdiff_t n) { return ut_CellIterator(myMatrix, myPos-n); }
72  ut_CellIterator &operator+=(ptrdiff_t n) { myPos += n; repage(); return *this; }
73  ut_CellIterator &operator-=(ptrdiff_t n) { myPos -= n; repage(); return *this; }
74  ut_CellIterator &operator++() { myPos++; myOffset++; myData++; if (myOffset >= CELL_PAGESIZE) repage(); return *this; }
75  ut_CellIterator &operator--() { myPos--; myOffset--; myData--; if (myOffset < 0) repage(); return *this; }
76  ut_CellIterator operator++(int) { ut_CellIterator result = *this; myPos++; myOffset++; myData++; if (myOffset >= CELL_PAGESIZE) repage(); return result; }
77  ut_CellIterator operator--(int) { ut_CellIterator result = *this; myPos--; myOffset--; myData--; if (myOffset < 0) repage(); return result; }
78  int operator-(ut_CellIterator b) const { return myPos - b.myPos; }
79  ut_MatrixCell &operator[](ptrdiff_t idx) { return myMatrix->getCell(idx + myPos); }
80 
81  bool operator<(ut_CellIterator b) const { return myPos < b.myPos; }
82  bool operator>(ut_CellIterator b) const { return myPos > b.myPos; }
83  bool operator>=(ut_CellIterator b) const { return myPos >= b.myPos; }
84  bool operator==(ut_CellIterator b) const { return myPos == b.myPos; }
85  bool operator!=(ut_CellIterator b) const { return myPos != b.myPos; }
86 
87  ut_MatrixCell &operator*() const { return *myData; }
88  ut_MatrixCell *operator->() const { return myData; }
89 
90  protected:
91 
92  void repage()
93  { myPage = myPos >> CELL_PAGEBITS; myOffset = myPos & CELL_PAGEMASK;
94  myData = 0;
95  if (myPage < myMatrix->myCellPages.entries())
96  myData = &myMatrix->myCellPages(myPage)[myOffset];
97  }
98 
99  const UT_SparseMatrixT<T, IsPaged> *myMatrix;
100  ut_MatrixCell *myData;
101  ptrdiff_t myPos;
102  int myPage, myOffset;
103  };
104 
105  // For compact compiled matrix cells, specify the number of values to
106  // pack into a single cell.
107  static const int CELLBITS = 2;
108  static const int CELLSIZE = 1 << CELLBITS;
109  static const int CELLMASK = CELLSIZE-1;
110 
111  class ut_4MatrixCell
112  {
113  public:
114  T myValue[CELLSIZE];
115  int myRow;
116  int myCol;
117  } SYS_ALIGN16;
118 
119 public:
121 
122  // creates a valid 'rows' x 'cols' zero matrix
123  UT_SparseMatrixT(int rows, int cols);
125  ~UT_SparseMatrixT();
126 
127  int getNumRows() const { return myNumRows; }
128  int getNumCols() const { return myNumCols; }
129 
130  /// Return the amount of memory used by this array.
131  int64 getMemoryUsage() const;
132  int64 getIdealMemoryUsage() const;
133 
134  /// Ensures at least this number of cells are available, useful
135  /// if you can predict your size.
136  void reserve(int numcells);
137 
138  /// Shrinks the array to exactly fit the number of cells present
139  void shrinkToFit();
140 
141  // sets equal to a 'rows' x 'cols' zero matrix
142  void init(int rows, int cols);
143 
144  // Sets all entries of the matrix to zero.
145  void zero();
146 
147  // Determines if we are large enough to justify mulithreading.
148  // The 5000 was determined from experiments with UT_Vector.
149  bool shouldMultiThread() const
150  {
151 #ifdef CELLBE
152  return false;
153 #else
154  return getNumRows() > 5000;
155 #endif
156  }
157 
158  // does M(row,col) += value
159  bool addToElement(int row, int col, T value);
160 
161  // Returns the index of the cell for which cell->row >= row is true
162  // first. This could be myCount if row is too high.
163  // myCells[result].myRow may be strictly greater than row if it isn't
164  // present.
165  int findCellFromRow(int row) const;
166 
167  // computes M * v
168  inline void multVec(const UT_VectorT<T> &v, UT_VectorT<T> &result) const
169  {
170  compile();
171  multVecInternal(v, result);
172  }
173 
174  /// Do result = result - (M * v)
175  inline void subtractMultVec(const UT_VectorT<T> &v, UT_VectorT<T> &result) const
176  {
177  compile();
178  subtractMultVecInternal(v, result);
179  }
180  inline void subtractMultVecNoThread(const UT_VectorT<T> &v, UT_VectorT<T> &result) const
181  {
182  compile();
183  subtractMultVecInternalNoThread(v, result);
184  }
185 
186  // computes M^T * v
187  void transposeMultVec(const UT_VectorT<T> &v, UT_VectorT<T> &result) const;
188 
189  /// Square of L2-norm of all columns
190  /// Use to form the diagonal part of A^tA of an overdetermined system A.
191  /// Needed for both Jacobi and Gauss-Seidel iterations.
192  void allColNorm2(UT_VectorT<T> &result) const;
193 
194  /// Initializes out to a matrix where the (i,j) elements are
195  /// (rowstart+i, colstart+j) elements of this. The rowend and
196  /// colend represent the exclusive end - [rowstart,rowend) will
197  /// be extracted for numrows = rowend - rowstart
198  /// This function will compile *this, if it wasn't before the call.
199  /// The resulting matrix will be compiled (by construction).
200  void extractSubMatrix(UT_SparseMatrixT<T, IsPaged> &out,
201  int rowstart, int rowend,
202  int colstart, int colend) const;
203  /// Extract submatrix without compiling *this.
204  /// The resulting submatrix won't be compiled either.
205  void extractSubMatrixUncompiled(UT_SparseMatrixT<T, IsPaged> &out,
206  int rowstart, int rowend,
207  int colstart, int colend) const;
208 
209  // Get fixed-size square submatrix
210  // This will trigger compilation
211  void getSmallSquareSubMatrix(T submatrix[],
212  const exint indices[], const int num_indices) const;
213 
214  // For each row, column, and value call visitor(row, column, value)
215  // visitor may be called with the same (row, column) pair multiple times
216  // the sum of all the associated values is the component in the matrix for (row, column)
217  template<typename Visitor>
218  void accept(Visitor& visitor) const
219  {
220  for(exint ci = 0; ci < myCount; ++ci)
221  {
222  ut_MatrixCell *cell = &getCell(ci);
223  visitor(cell->myRow, cell->myCol, cell->myValue);
224  }
225  }
226 
227  /// Clear the rows and columns specified by the given bit field to a row
228  /// from the identity matrix
229  void clearRowsAndColumns(const UT_BitArray &toclear);
230 
231  /// Extract the diagonal vector from a matrix. This does not
232  /// require the matrix to be compiled, nor does it try to compile it.
233  /// For entries (i, j), this method will extract the entries where
234  /// (i-j == idx).
235  void extractDiagonal(UT_VectorT<T> &out, int idx = 0) const;
236 
237  /// Extract everything but the diagonal:
238  /// out(i, j) = 0, if i == j,
239  // *this(i, j), otherwise
240  void extractNondiagonal(UT_SparseMatrixT<T, IsPaged> &out) const;
241 
242  /// Incomplete Cholesky Factorization.
243  /// Does a Cholesky Factorization but does not make any elements non-zero
244  /// The result of this operation is an UPPER triangular matrix G
245  /// such that A = Gt G.
246  /// The input must be a symmetric matrix.
247  /// Note that this factorization is not always stable.
248  /// Returns 0 on success, -1 if semi-definite (diagonal contained zeros
249  /// within the given tolerance) and -2 if not positive definite (diagonal
250  /// contained negative numbers)
251  int incompleteCholeskyFactorization(T tol=1e-5);
252 
253  /// Modified Incomplete Cholesky
254  /// Same as incomplte cholesky, except attempt to account for the
255  /// discarded entries by adjusting the diagonal.
256  /// tau is a tuning constant.
257  int modifiedIncompleteCholesky(T tau = 0.97,
258  T mindiagratio = 0.25,
259  T tol=1e-5);
260 
261  /// Assumes this is a lower triangular matrix. Solves the equation
262  /// A x = b
263  /// If the diagonal of A is zero within tolerance, the corresponding
264  /// x coordinate is zero.
265  /// Returned is the number of artifical zeros places into x. 0 means
266  /// the solution encountered no singularities, 10 would mean 10
267  /// singularities.
268  int solveLowerTriangular(UT_VectorT<T> &x, const UT_VectorT<T> &b,
269  T tol=1e-5) const;
270  int solveUpperTriangular(UT_VectorT<T> &x, const UT_VectorT<T> &b,
271  T tol=1e-5) const;
272  /// Given an upper triangular matrix, solves the lower triangular
273  /// transposed of it and negates the result.
274  int solveLowerTriangularTransposeNegate(UT_VectorT<T> &x,
275  const UT_VectorT<T> &b,
276  T tol=1e-5) const;
277 
278  // Use the conjugate gradient method to solve Ax=b.
279  // NB: 'x' should be initialized with the initial estimate
280  bool solveConjugateGradient(UT_VectorT<T> &x, const UT_VectorT<T> &b,
282  bool (*callback_func)(void *) = 0,
283  void *callback_data = 0, T tol2=1e-5,
284  int max_iters = -1) const;
285 
286  /// Transposes this matrix.
287  THREADED_METHOD(UT_SparseMatrixT, shouldMultiThread(), transpose)
288  void transposePartial(const UT_JobInfo &info);
289 
290  /// Makes this a transposed copy of source. By not working
291  /// in place we can avoid a sort if the source is compiled.
292  /// If source is not compiled, we compile it.
293  void transposeCompiled(const UT_SparseMatrixT<T, IsPaged> &src);
294 
295  /// *this = -*this.
296  THREADED_METHOD(UT_SparseMatrixT, shouldMultiThread(), negate)
297  void negatePartial(const UT_JobInfo &info);
298 
299  UT_SparseMatrixT<T, IsPaged> &operator=(const UT_SparseMatrixT<T, IsPaged> &m);
300  UT_SparseMatrixT<T, IsPaged> &operator*=(T scalar);
301  UT_SparseMatrixT<T, IsPaged> &operator+=(const UT_SparseMatrixT<T, IsPaged> &m);
302 
303  // Writes out the full matrix as an Row X Column array. Likely
304  // not too useful with big matrices...
305  void printFull(std::ostream &os) const;
306 
307  // Writes out all of the non-zero elements.
308  void printSparse(std::ostream &os) const;
309 
310  // Writes out all of the non-zero elements as a Matlab "sparse" expression.
311  void printSparseMatlab(std::ostream &os,
312  const UT_String &varname) const;
313 
314  // Binary save and load of a sparse matrix
315  void save(std::ostream &os) const;
316  void load(UT_IStream &is);
317 
318  /// Reorders the cell array to be sorted by row than column.
319  /// Consolidates any duplicate entries by adding them together.
320  /// Culls any zero entries.
321  /// compile will be a no-op if the matrix's structure has not changed
322  /// since the last invocation.
323  /// While compilation is expensive, it also performs the collapse/cull
324  /// operation, so you may want to invoke it explicitly if you are
325  /// repeatedly gathering similar matrices, say with +=.
326  void compile() const;
327  bool isCompiled() const { return myCompiledFlag; }
328  bool isStillSorted() const { return myStillSortedFlag; }
329 
330  bool hasNan() const;
331  void testForNan() const;
332 
333 private:
334  /// Always recompiles the matrix regardless of the compile flag setting.
335  void forceCompile() const;
336 
337  /// Extract a cell by its index.
338  inline ut_MatrixCell &getCell(int idx) const
339  {
340  if (IsPaged)
341  return myCellPages(idx >> CELL_PAGEBITS)[idx & CELL_PAGEMASK];
342  else
343  return myCells[idx];
344  }
345 
346  THREADED_METHOD2_CONST(UT_SparseMatrixT, shouldMultiThread(),
347  multVecInternal, const UT_VectorT<T> &, v,
348  UT_VectorT<T> &, result)
349  void multVecInternalPartial(const UT_VectorT<T> &v, UT_VectorT<T> &result,
350  const UT_JobInfo &info) const;
351 
352  THREADED_METHOD2_CONST(UT_SparseMatrixT, shouldMultiThread(),
353  subtractMultVecInternal,
354  const UT_VectorT<T>&, v,
355  UT_VectorT<T>&, result)
356  void subtractMultVecInternalPartial(const UT_VectorT<T>& v, UT_VectorT<T>& result,
357  const UT_JobInfo& info) const;
358 
359  /// Compile into 4-value cell data structure
360  void compile4() const;
361 
362  int myNumRows;
363  int myNumCols;
364  mutable UT_ValArray<ut_MatrixCell *> myCellPages;
365  mutable ut_MatrixCell *myCells;
366  mutable int myCount;
367  mutable int *my4RowOffsets;
368  mutable ut_4MatrixCell *my4Cells;
369  mutable int my4Count;
370  int myMaxSize;
371  mutable bool myCompiledFlag;
372  mutable bool myStillSortedFlag;
373  mutable bool my4CompiledFlag;
374 
375  friend class UT_SparseMatrixRowT<T>;
376 };
377 
378 typedef UT_SparseMatrixT<fpreal32, false> UT_SparseMatrixF;
379 typedef UT_SparseMatrixT<fpreal64, false> UT_SparseMatrixD;
380 typedef UT_SparseMatrixT<fpreal64, false> UT_SparseMatrix;
381 
382 
383 ///
384 /// This is a highly specialized varient of the SparseMatrix
385 /// which does not let you change the number or position of cells
386 /// It has, however, compiled fixed row offsets to make it fast
387 /// to look up specific rows in the matrix.
388 ///
389 template <typename T>
391 {
392  class ut_MatrixCell
393  {
394  public:
395  int myCol;
396  T myValue;
397  };
398 
399 public:
401  ~UT_SparseMatrixRowT();
402 
403  // Can't mix paged and non paged, and this is a non-paged algorithm
404  // Destroys the source matrix in the process as we reuse its
405  // myCells array.
406  void buildFrom(UT_SparseMatrixT<T, false> &m,
407  bool invertdiag = false,
408  T tol=1e-5f);
409 
410  int getNumRows() const { return myNumRows; }
411  int getNumCols() const { return myNumCols; }
412 
413  /// Return the amount of memory used by this array.
414  int64 getMemoryUsage() const;
415 
416  // The 5000 was determined from experiments with UT_Vector.
417  bool shouldMultiThread() const
418  {
419  return getNumRows() > 5000;
420  }
421 
422  // Returns the index of the cell for which cell->row >= row is true
423  // first. This could be myCount if row is too high.
424  // myCells[result].myRow may be strictly greater than row if it isn't
425  // present.
426  int findCellFromRow(int row) const
427  { return myRowOffsets(row); }
428 
429  // Get the inverted diagonal.
430  void getDiagonalInv(UT_VectorT<T> &out) const;
431 
432  // computes M * v
433  THREADED_METHOD2_CONST(UT_SparseMatrixRowT, shouldMultiThread(), multVec,
434  const UT_VectorT<T> &, v,
435  UT_VectorT<T> &, result)
436  void multVecPartial(const UT_VectorT<T> &v, UT_VectorT<T> &result,
437  const UT_JobInfo &info) const;
438 
439  // computes M * v and accumulates the dot product of v & result
440  void multVecAndDot(
441  const UT_VectorT<T> &v,
442  UT_VectorT<T> &result,
443  fpreal64 *dotpq) const;
444 
445  /// Assumes this is a lower triangular matrix. Solves the equation
446  /// A x = b
447  /// If the diagonal of A is zero within tolerance, the corresponding
448  /// x coordinate is zero.
449  /// Returned is the number of artifical zeros places into x. 0 means
450  /// the solution encountered no singularities, 10 would mean 10
451  /// singularities.
452  int solveUpperTriangular(UT_VectorT<T> &x, const UT_VectorT<T> &b,
453  T tol=1e-5) const;
454  /// Given an upper triangular matrix, solves the lower triangular
455  /// transposed of it and negates the result.
456  int solveLowerTriangularTransposeNegate(UT_VectorT<T> &x,
457  const UT_VectorT<T> &b,
458  T tol=1e-5) const;
459 
460  /// Solves conjugate gradient using our specialized functions.
461  /// These allow us to perform some normal and dot operations
462  /// while the cache is still hot.
463  /// This matrix is the matrix to solve. The provided GT matrix
464  /// is the upper triangular result of cholesky factoriziation.
465  /// Norm is hardcoded to 2.
466  /// If the GT matrix is null, no preconditioner will be used.
467  float solveConjugateGradient(UT_VectorT<T> &x, const UT_VectorT<T> &b,
468  const UT_SparseMatrixRowT<T> *GT, T tol2=1e-5,
469  int max_iters = -1, int *iterout = NULL) const;
470 private:
471  THREADED_METHOD3_CONST(UT_SparseMatrixRowT, shouldMultiThread(),
472  multVecAndDotInternal,
473  const UT_VectorT<T> &, v,
474  UT_VectorT<T> &, result,
475  fpreal64 *, dotpq)
476  void multVecAndDotInternalPartial(const UT_VectorT<T> &v, UT_VectorT<T> &result,
477  fpreal64 *dotpq,
478  const UT_JobInfo &info) const;
479 
480  static const exint PARALLEL_BLOCK_SIZE = 1024;
481 
482  /// Extract a cell by its index.
483  inline ut_MatrixCell &getCell(int idx) const
484  {
485  return myCells[idx];
486  }
487 
488  int myNumRows;
489  int myNumCols;
490  ut_MatrixCell *myCells;
491  int myCount;
492  UT_IntArray myRowOffsets;
493  UT_VectorT<T> myDiagonal;
494 };
495 
498 
499 
500 
501 /// This is another highly specialized varient of the SparseMatrix
502 /// which allows multithreaded addition of elements to individual rows and
503 /// supports a conjugate gradient solve. This can be useful for solvers which
504 /// need to create and solve large linear systems in parallel.
505 /// It uses the ELL format which is appropriate for matrices with a limited
506 /// but fairly consistent numer of non-zeros per row, for example the result
507 /// of discretizations across regular grids.
508 template <typename T, bool colmajor = false, bool useexint = false>
510 {
511  typedef typename std::conditional<useexint, exint, int>::type inttype;
512 
513 public:
514 
516 
517  UT_SparseMatrixELLT(inttype rows, int nzeros);
518 
519  void init(inttype rows, int nzeros);
520 
521  void zero();
522 
523  exint getNumRows() const { return myRows; }
524 
525  int getNonZerosPerRow() const { return myNonZeros; }
526 
527  const UT_ValArray<T> &getRowValues() const { return myRowVals; }
528  UT_ValArray<T> &getRowValues() { return myRowVals; }
529 
530  const UT_ValArray<inttype> &getColumns() const { return myColumns; }
531  UT_ValArray<inttype> &getColumns() { return myColumns; }
532 
533  exint index(inttype row, int nz) const
534  {
535  if (colmajor)
536  return (exint)nz * myRows + (exint)row;
537  return (exint)row * (exint)myNonZeros + (exint)nz;
538  }
539 
540  bool appendRowElement(inttype row, inttype col, T val, int &rowidx)
541  {
542  UT_ASSERT(rowidx < myNonZeros);
543  exint idx = index(row, rowidx);
544  myRowVals(idx) = val;
545  myColumns(idx) = col;
546  rowidx++;
547  return true;
548  }
549 
550  // Sort rows by column index.
551  THREADED_METHOD(UT_SparseMatrixELLT, shouldMultiThread(),
552  sortRows)
553  void sortRowsPartial(const UT_JobInfo &info);
554 
555  /// Every reference to given rows is removed from the diagonal.
556  /// Because searching the whole matrix is expensive, assumes it
557  /// is a standard fluid matrix, so for each specified row, it
558  /// looks at all columns that are not itself and removes those.
559  void removeEffectOfRows(const UT_ExintArray &rows);
560 
561  bool shouldMultiThread() const
562  {
563  return getNumRows() > 5000;
564  }
565 
566  /// Copy to the supplied SparseMatrix.
567  void copyToSparse(UT_SparseMatrixT<T, false> &A) const;
568 
569  // Extract the inverted diagonal
570  THREADED_METHOD1_CONST(UT_SparseMatrixELLT, shouldMultiThread(),
571  getDiagonalInv, UT_VectorT<T> &, out)
572 
573  /// Get the inverse of the diagonal. Required for the conjugate gradient
574  /// solve when using Jacobi pre-conditioning.
575  void getDiagonalInvPartial(UT_VectorT<T> &out,
576  const UT_JobInfo &info) const;
577 
578  // computes M * v
579  THREADED_METHOD2_CONST(UT_SparseMatrixELLT, shouldMultiThread(),
580  multVec, const UT_VectorT<T> &, v,
581  UT_VectorT<T> &, result)
582 
583  void multVecPartial(const UT_VectorT<T> &v, UT_VectorT<T> &result,
584  const UT_JobInfo &info) const;
585 
586  // computes M * v and accumulates the dot product of v & result
587  void multVecAndDot(
588  const UT_VectorT<T> &v,
589  UT_VectorT<T> &result,
590  fpreal64 *dotpq) const;
591 
592  // computes M * v and accumulates the dot product of v & result
593  void multVecAndDotUpTo(
594  const UT_VectorT<T> &v,
595  UT_VectorT<T> &result,
596  fpreal64 *dotpq,
597  exint solverbase) const;
598 
599  /// Solves conjugate gradient using our specialized functions.
600  /// These allow us to perform some normal and dot operations
601  /// while the cache is still hot.
602  /// This matrix is the matrix to solve. The provided GT matrix
603  /// is the upper triangular result of cholesky factoriziation.
604  /// Norm is hardcoded to 2.
605  /// If the GT matrix is null, no preconditioner will be used.
606  float solveConjugateGradient(UT_VectorT<T> &x, const UT_VectorT<T> &b,
607  const UT_SparseMatrixRowT<T> *GT, T tol2=1e-5,
608  int max_iters = -1, int *iterout = NULL) const;
609 
610  // Writes out all of the non-zero elements as a Matlab "sparse" expression.
611  void printSparseMatlab(std::ostream &os,
612  const UT_String &varname) const;
613 
614  // Writes out all of the non-zero elements in Matrix Market format.
615  void printSparseMatrixMarket(std::ostream &os) const;
616 
617 private:
618 
619  THREADED_METHOD3_CONST(UT_SparseMatrixELLT, shouldMultiThread(),
620  multVecAndDotInternal,
621  const UT_VectorT<T> &, v,
622  UT_VectorT<T> &, result,
623  fpreal64 *, dotpq)
624 
625  void multVecAndDotInternalPartial(const UT_VectorT<T> &v, UT_VectorT<T> &result,
626  fpreal64 *dotpq,
627  const UT_JobInfo &info) const;
628 
629  THREADED_METHOD4_CONST(UT_SparseMatrixELLT, shouldMultiThread(),
630  multVecAndDotUpToInternal,
631  const UT_VectorT<T> &, v,
632  UT_VectorT<T> &, result,
633  fpreal64 *, dotpq,
634  exint, solverbase)
635 
636  void multVecAndDotUpToInternalPartial(const UT_VectorT<T> &v, UT_VectorT<T> &result,
637  fpreal64 *dotpq, exint solverbase,
638  const UT_JobInfo &info) const;
639 
640  static const exint PARALLEL_BLOCK_SIZE = 1024;
641 
642 private:
643 
644  UT_ValArray<T> myRowVals;
645  UT_ValArray<inttype> myColumns;
646  exint myRows;
647  int myNonZeros;
648 };
649 
650 typedef UT_SparseMatrixELLT<fpreal32> UT_SparseMatrixELLF;
651 typedef UT_SparseMatrixELLT<fpreal64> UT_SparseMatrixELLD;
652 
653 template <typename T> class UT_API UT_SparseMatrixBuilderT;
654 
655 /// Sparse matrix class that uses the compressed sparse row storage scheme.
656 /// This is a general-purpose sparse matrix class that can contain an arbitrary
657 /// number of non-zero elements per row.
658 template <typename T>
660 {
661 public:
662  /// Triplet to construct matrix. Is used only for setValues and addValues
663  /// to construct the matrix.
664  struct Triplet {
665  exint myRow;
666  exint myCol;
668 
669  /// Construct a zero triplet.
670  Triplet() : myRow(0), myCol(0), myValue(0) {}
671 
672  /// Construct a triplet for the given row and column with value v.
673  Triplet(exint r, exint c, T v) : myRow(r), myCol(c), myValue(v) {}
674 
675  /// Returns true if the triplet's coordinate comes before the given
676  /// triplet's coordinate.
677  inline bool operator<(const Triplet &t) const {
678  return (myRow < t.myRow || (myRow == t.myRow && myCol < t.myCol));
679  }
680  };
681 
682  /// Defines one cell in the matrix. Can also be used to directly construct
683  /// the matrix with setValuesFromRows, which is easier to multithread.
684  struct ColumnValue {
685  exint myCol;
687 
688  /// Construct a zero column-value.
689  ColumnValue() : myCol(0), myValue(0) {}
690 
691  /// Construct a column value for the given column with value v.
692  ColumnValue(exint c, T v) : myCol(c), myValue(v) {}
693 
694  /// Returns true if the column-value's column comes before the given
695  /// column-value's column.
696  inline bool operator<(const ColumnValue &cv) const {
697  return myCol < cv.myCol;
698  }
699  };
700 
701  /// Constructs an uninitialized sparse matrix
702  UT_SparseMatrixCSRT();
703 
704  /// Constructs a sparse matrix initialized with the given rows and columns.
705  UT_SparseMatrixCSRT(exint rows, exint columns);
706 
707  /// Constructs a sparse matrix that is a deep copy of the given sparse
708  /// matrix.
709  explicit UT_SparseMatrixCSRT(const UT_SparseMatrixCSRT &matrix);
710 
711  /// Sets the size of the matrix with the given rows and columns and clears
712  /// the matrix.
713  void init(exint rows, exint columns);
714 
715  /// Clears all data in the matrix
716  void zero();
717 
718  /// Swaps size and data elements with the given sparse matrix.
719  void swap(UT_SparseMatrixCSRT &other);
720 
721  /// The number of rows in the matrix.
722  exint getNumRows() const { return myNumRows; }
723 
724  /// The number of columns in the matrix.
725  exint getNumCols() const { return myNumCols; }
726 
727  /// The total number of non-zero elements in the matrix.
728  exint getNumNonZeros() const { return myCells.size(); }
729 
730  /// The number of non-zero elements in the given row.
731  exint getNumCellsInRow(exint row) const { return rowEnd(row) - rowBegin(row); }
732 
733  /// Clears then fills the matrix with the given array of triplets.
734  /// This sorts the triplets array.
735  /// Triplets with the same indices will be summed together.
736  /// A cell whose value is smaller than the tolerance will be deleted.
737  void setValues(UT_Array<Triplet> &triplets, T tolerance = 0.0);
738 
739  /// Clears then fills the matrix with the given array of ColumnValue arrays.
740  /// rowValues must have the same size as the number of rows in this matrix.
741  /// This data layout makes it easier to multithread the construction of the
742  /// sparse matrix entries.
743  void setValuesFromRows(UT_Array<UT_Array<ColumnValue>> &rowValues,
744  T tolerance = 0.0);
745 
746  /// Adds the given cells to the matrix. Triplets with the same coordinates
747  /// as cells already present in the matrix will be added together.
748  /// The array of triplets will be sorted unless alreadySorted is true.
749  /// Cells whose magnitude is below the tolerance value after the addition
750  /// will be discarded
751  void addValues(UT_Array<Triplet> &triplets, bool alreadySorted = false,
752  T tolerance=1e-5);
753 
754  /// Writes the diagonal values of the matrix into the given vector
755  /// If invert is true, the values will be inverted (i.e. 1/x).
756  /// If a row diagonal's value <= tolerance, then defaultValue is placed
757  /// in its result row.
758  void extractDiagonal(UT_VectorT<T> &result, bool invert = false,
759  T tolerance = 0, T defaultValue = 0) const;
760 
761  /// Matrix-vector multiplication. Computes result = M * v.
762  void multVec(const UT_VectorT<T> &v, UT_VectorT<T> &result) const;
763 
764  /// Matrix-vector multiplication.Computes result = M * v and returns the dot
765  /// product of v and the result.
766  UT_PreciseT<T> multVecAndDot(const UT_VectorT<T> &v, UT_VectorT<T> &result) const;
767 
768  /// Matrix-vector multiplication. Computes result = y - Ax. result can refer
769  /// to the same vector as y.
770  void subtractMultVec(const UT_VectorT<T> &x, const UT_VectorT<T> &y,
771  UT_VectorT<T> &result) const;
772 
773  /// Matrix-matrix multiplication. Computes result = M * B. Resultant
774  /// elements whose values are smaller than the given tolerance are deleted.
775  void multMatrix(const UT_SparseMatrixCSRT<T> &B, UT_SparseMatrixCSRT &result,
776  T tolerance = 0.0) const;
777 
778  /// Computes the transpose of this matrix. result can refer to this matrix.
779  void transpose(UT_SparseMatrixCSRT<T> &result) const;
780 
781  /// Element-wise scalar multiplication. Every element in this matrix is
782  /// multiplied by the given scalar.
783  void scale(const T& scalar);
784 
785  /// Element-wise scalar multiplication. Every element in row i of the matrix
786  /// is multiplied by scalars[i]
787  void scaleRows(const UT_VectorT<T> &scalars);
788 
789  /// Multiples each row such that the row diagonal becomes 1. If the row
790  /// diagonal is below the tolerance, the row is unchanged.
791  void normalizeRows(T tolerance = 1e-5);
792 
793  /// Calculates the norms of each row. result must have the same size as the
794  /// number of rows in this matrix. p determines the exponent of each column.
795  /// i.e. result[i] = SUM( |Aij|^p ) for j in row(i)
796  /// If invert is true, then inverts each sum
797  /// i.e. result[i] = 1.0 / SUM( |Aij|^p ) for j in row(i)
798  /// If the row norm is smaller than the tolerance, defaultValue is
799  /// placed in its result entry instead.
800  void rowNorms(UT_VectorT<T> &result, int p = 2, bool invert = false,
801  T tolerance = 0, T defaultValue = 0) const;
802 
803  /// Analyzes the matrix and constructs a directed acyclic graph as a level
804  /// set, where each level contains rows that are independent of each other.
805  /// More specifically, level i of the levelSet contains rows that are only
806  /// dependent on rows present in levels [0, i-1].
807  /// The size of the levelSet is the number of levels.
808  /// If buildForUpper is true, this constructs a level set for the upper
809  /// triangular part of the matrix. Otherwise, this constructs a level set
810  /// for the lower triangular part of the matrix.
812  void buildDependencyLevelSet(LevelSet &levelSet, bool buildForUpper = false) const;
813 
814  /// Assumes this is a lower triangular matrix. Solves the equation Ax = b.
815  /// Values may be present in the upper triangular part of the matrix
816  /// but they will be ignored.
817  /// If unitDiagonal is true, then 1 will be used instead of the
818  /// current diagonal values.
819  /// If a diagonal of A is zero within the tolerance, 0 will be placed
820  /// in the corresponding row for x.
821  /// If negate is true, then x is negated
822  /// x can refer to the same vector as b to solve in-place.
823  void solveLowerTriangular(UT_VectorT<T> &x, const UT_VectorT<T> &b,
824  bool unitDiagonal = false, T tolerance = 1e-5,
825  bool negate = false) const;
826 
827  /// Assumes this is an upper triangular matrix. Solves the equation Ax = b.
828  /// Values may be present in the lower triangular part of the matrix
829  /// but they will be ignored.
830  /// If unitDiagonal is true, then 1 will be used instead of the
831  /// current diagonal values.
832  /// If a diagonal of A is zero within the tolerance, 0 will be placed
833  /// in the corresponding row for x.
834  /// If negate is true, then x is negated
835  /// x can refer to the same vector as b to solve in-place.
836  void solveUpperTriangular(UT_VectorT<T> &x, const UT_VectorT<T> &b,
837  bool unitDiagonal = false, T tolerance = 1e-5,
838  bool negate = false) const;
839 
840  /// Solves the lower triangular system A^T x = b.
841  /// A is assumed to be an upper triangular matrix.
842  /// Values may be present in the lower triangular part of the matrix
843  /// but they will be ignored.
844  /// If unitDiagonal is true, then 1 will be used instead of the
845  /// current diagonal values.
846  /// If a diagonal of A is zero within the tolerance, 0 will be placed
847  /// in the corresponding row for x.
848  /// If negate is true, then x is negated.
849  void solveLowerTriangularTranspose(UT_VectorT<T> &x, const UT_VectorT<T> &b,
850  bool unitDiagonal = false, T tolerance = 1e-5,
851  bool negate = false) const;
852 
853  /// Solves the upper triangular system A^T x = b.
854  /// A is assumed to be a lower triangular matrix.
855  /// Values may be present in the upper triangular part of the matrix
856  /// but they will be ignored.
857  /// If unitDiagonal is true, then 1 will be used instead of the
858  /// current diagonal values.
859  /// If a diagonal of A is zero within the tolerance, 0 will be placed
860  /// in the corresponding row for x.
861  /// If negate is true, then x is negated.
862  void solveUpperTriangularTranspose(UT_VectorT<T> &x, const UT_VectorT<T> &b,
863  bool unitDiagonal = false, T tolerance = 1e-5,
864  bool negate = false) const;
865 
866  /// Factorizes this matrix into a triangular matrix M such that M*M^T ~= A
867  /// This keeps the non-zero structure of the matrix. The matrix must be
868  /// symmetric.
869  /// This is unstable depending on the diagonal; Returns -1 if the diagonal
870  /// contains zeros within the given tolerance, and -2 if not positive
871  /// definite (diagonal has negative numbers).
872  /// Returns 0 on success.
873  /// If useModified = true, then this uses the Modified Incomplete Cholesky
874  /// Factorization, where discarded entries are accounted for by adjusting
875  /// the diagonal. tau is a tuning constant that scales the adjustment from
876  /// the missing entries. mindiagratio is the relative tolerance at which
877  /// the original diagonal will be used if the matrix becomes negative
878  /// definite as a result of the computation.
879  int incompleteCholeskyFactorization(T tol=1e-5, bool useModified = false,
880  T tau = 0.97, T mindiagratio = 0.25);
881 
882  /// Factorizes this matrix into L and U where LU approximates this matrix.
883  /// The non-zero structure of LU matches this matrix; that is, this produces
884  /// ILU(0) of this matrix. A lower triangular level set is needed for
885  /// multithreading.
886  /// e.g.
887  /// LevelSet levelSet;
888  /// A.buildDependencyLevelSet(levelSet);
889  /// A.incompleteLUFactorization(L, U, 1e-5, levelSet);
890  void incompleteLUFactorization(UT_SparseMatrixCSRT<T> &L,
891  UT_SparseMatrixCSRT<T> &U, T tolerance,
892  const LevelSet &workUnits, UT_Interrupt *boss = NULL) const;
893 
894  /// Initializes out to a matrix with size (rowend-rowstart, colend-colstart)
895  /// and fills it with the elements defined by the submatrix in rows
896  /// [rowstart, rowend), [colstart, colend).
897  void extractSubMatrix(UT_SparseMatrixCSRT<T> &out,
898  exint rowstart, exint rowend, exint colstart, exint colend) const;
899 
900  /// Calls inspect(row, col, value) on each cell. Each (row, col) is visited
901  /// exactly once. This is done in parallel over the rows by default, but can
902  /// be done serially by setting serial = true. inspect can modify the cell
903  /// values if it has the signature (exint, exint, &T).
904  template <typename CellFunction>
905  void visit(const CellFunction& inspect, bool serial = false)
906  {
907  parallelExecute(UT_BlockedRange<exint>(0, myNumRows), [this, &inspect]
908  (const UT_BlockedRange<exint> &r)
909  {
910  for (exint row = r.begin(); row < r.end(); ++row)
911  for (exint i = rowBegin(row), n = rowEnd(row); i < n; ++i)
912  inspect(row, cellColumn(i), cellValue(i));
913  }, serial, myLightRowOpSubRatio, myLightRowOpGrainSize);
914  }
915 
916  /// Condition is called for each cell exactly once. If the condition returns
917  /// true, then the cell is removed from the matrix. The function signature
918  /// of the condition should be (exint row, exint col, T value) -> bool.
919  /// condition can modify the cell vakyes if it has the signature
920  /// (exint, exint, &T)
921  template <typename CellFunction>
922  void removeIf(const CellFunction& condition)
923  {
924  exint idx = -1, i = rowBegin(0);
925  for (exint row = 0; row < myNumRows; ++row)
926  {
927  for (exint n = rowEnd(row); i < n; ++i)
928  if (!condition(row, cellColumn(i), cellValue(i)) && ++idx != i)
929  myCells[idx] = myCells[i];
930  myRowOffsets[row+1] = idx+1;
931  }
932  myCells.setSize(idx+1);
933  }
934 
935  /// Returns true if any non-zero element in this matrix has value NaN.
936  bool hasNan() const;
937 
938  /// Asserts that the matrix contains no element with value NaN.
939  void testForNan() const;
940 
941  /// Writes the full matrix, including zeros, to the given output stream.
942  void printFull(std::ostream &os) const;
943 
944  /// Serializes the matrix to the given output stream.
945  void save(std::ostream &os) const;
946 
947  /// Deserializes the matrix from the given input stream.
948  void load(UT_IStream &is);
949 
950  /// Deep copies all data from the given matrix.
951  UT_SparseMatrixCSRT<T> &operator=(const UT_SparseMatrixCSRT<T> &m);
952 
953  /// Equivalent to A.scale(scalar)
954  UT_SparseMatrixCSRT<T> &operator*=(const T& scalar) { scale(scalar); return *this; }
955 
956 private:
957 
958  // Helper functions to get value and column
959  inline T& cellValue(exint offset) { return myCells[offset].myValue; }
960  inline const T &cellValue(exint offset) const { return myCells[offset].myValue; }
961  inline const exint &cellColumn(exint offset) const { return myCells[offset].myCol; }
962 
963  // Helper functions to return the beginning and end indexes into myValues
964  // for the given row. This can be used to iterate over the non-zero values
965  // of a particular row.
966  // e.g. for (exint i = rowBegin(row); i < rowEnd(row); ++i)
967  // ... Do something with cellValue(i) and cellColumn(i)
968  inline const exint &rowBegin(exint row) const { return myRowOffsets[row]; }
969  inline const exint &rowEnd(exint row) const { return myRowOffsets[row + 1]; }
970 
971  // Multithreading constants
972  // If the number of rows is less than this number, all operations run
973  // serially
974  static const int myParallelThreshold = 1000;
975 
976  // The subratio and grainsize used for row operations with cheap work units
977  static const int myLightRowOpSubRatio = 0;
978  static const int myLightRowOpGrainSize = 64;
979 
980  // Executes the given body over the given range. Can force this method to
981  // run serially by setting serial = true
982  template <typename Body>
983  inline void parallelExecute(const UT_BlockedRange<exint> &range,
984  const Body &body, bool serial = false, int subratio = 0,
985  int grainsize = 32, exint parallelThreshold = -1) const {
986  // Use internal parallel threshold by default
987  if (parallelThreshold < 0)
988  parallelThreshold = myParallelThreshold;
989  // Perform serially if specified or if number of rows is below threshold
990  if (myNumRows < parallelThreshold || serial)
991  UTserialFor(range, body);
992  else
993  UTparallelFor(range, body, subratio, grainsize);
994  };
995 
996 private:
997  exint myNumRows;
998  exint myNumCols;
999 
1000  // Stores the non-zero cells of the matrix, sorted by row and column.
1001  // Has size Nz, where Nz is the number of non-zero values.
1002  UT_Array<ColumnValue> myCells;
1003 
1004  // Stores the offset of myValues at which each row starts.
1005  // Has size myNumRows + 1.
1006  // myRowOffsets[r] is the offset of myValues where row r begins.
1007  // myRowOffsets[r + 1] - myRowOffsets[r] is the number of
1008  // non-zero values in row r.
1009  // myRowOffsets[myNumRows] has the sentinel value Nz
1010  UT_ExintArray myRowOffsets;
1011 
1012  friend class UT_SparseMatrixBuilderT<T>;
1013 };
1014 
1017 
1018 /// Simple helper class to construct a SparseMatrixCSRT. This operates on the
1019 /// matrix data in-place and is generally much more performant than constructing
1020 /// via triplets, but the builder has the limitation of only being able to add
1021 /// to rows sequentially; i.e. each call to startRow must start the construction
1022 /// of a row greater than the last, while addToColumn can be called with any
1023 /// column ordering.
1024 /// 1. UT_SparseMatrixCSRD matrix(M, N);
1025 /// 2. UT_SparseMatrixBuilderT builder(matrix);
1026 /// 3. for (exint r = startRow; r < endRow; ++r)
1027 /// 4. builder.startRow(r);
1028 /// 5. for (any column)
1029 /// 6. builder.addToColumn(column, value);
1030 /// 7. builder.finish();
1031 template <typename T>
1032 class UT_API UT_SparseMatrixBuilderT
1033 {
1034 public:
1036 
1037  /// Construct an uninitialized sparse matrix builder.
1038  UT_SparseMatrixBuilderT() : myMatrix(NULL) {}
1039 
1040  /// Construct a sparse matrix builder for the given matrix.
1041  UT_SparseMatrixBuilderT(SparseMatrixType &matrix) { init(matrix); }
1042 
1043  /// Initializes the builder to start building the given matrix
1044  void init(SparseMatrixType &matrix);
1045 
1046  /// Sets the minimum capacity of non-zero elements for the given matrix.
1047  /// This is used to preallocate memory for the matrix to avoid excessive
1048  /// reallocation as elements are added.
1049  void setCapacity(exint capacity);
1050 
1051  /// Advances the current row to the given row. The row must be greater than
1052  /// the current row.
1053  void startRow(exint row);
1054 
1055  /// Adds the given value at the given column to the current row. If a cell
1056  /// is already present for the given column, the given value is added to it.
1057  void addToColumn(exint col, T value);
1058 
1059  /// Completes the remaining rows and finalizes the matrix.
1060  void finish();
1061 
1062 private:
1063  exint myCurrentRow;
1064  SparseMatrixType *myMatrix;
1065 };
1066 
1069 
1070 #endif
int getNumRows() const
int getNumRows() const
#define THREADED_METHOD4_CONST(CLASSNAME, DOMULTI, METHOD, PARMTYPE1, PARMNAME1, PARMTYPE2, PARMNAME2, PARMTYPE3, PARMNAME3, PARMTYPE4, PARMNAME4)
Color4< T > operator*(S a, const Color4< T > &v)
Definition: ImathColor.h:728
exint getNumCols() const
The number of columns in the matrix.
ColumnValue()
Construct a zero column-value.
GLsizei GLenum const void * indices
Definition: glcorearb.h:405
GLenum GLint * range
Definition: glcorearb.h:1924
bool operator<(const ColumnValue &cv) const
UT_SparseMatrixRowT< fpreal64 > UT_SparseMatrixRowD
GLboolean invert
Definition: glcorearb.h:548
exint getNumNonZeros() const
The total number of non-zero elements in the matrix.
exint getNumRows() const
const GLdouble * v
Definition: glcorearb.h:836
bool shouldMultiThread() const
typedef void(APIENTRYP PFNGLCULLFACEPROC)(GLenum mode)
const UT_ValArray< inttype > & getColumns() const
UT_SparseMatrixCSRT< T > SparseMatrixType
Plane3< T > operator-(const Plane3< T > &plane)
Definition: ImathPlane.h:249
ColumnValue(exint c, T v)
Construct a column value for the given column with value v.
class UT_API UT_SparseMatrixRowT
#define UT_API
Definition: UT_API.h:13
GLint y
Definition: glcorearb.h:102
bool isStillSorted() const
png_uint_32 i
Definition: png.h:2877
void subtractMultVecNoThread(const UT_VectorT< T > &v, UT_VectorT< T > &result) const
int getNumCols() const
UT_SparseMatrixRowT< fpreal32 > UT_SparseMatrixRowF
bool operator<(const Triplet &t) const
UT_SparseMatrixCSRT< T > & operator*=(const T &scalar)
Equivalent to A.scale(scalar)
GLsizei GLboolean transpose
Definition: glcorearb.h:831
void multVec(const UT_VectorT< T > &v, UT_VectorT< T > &result) const
bool operator==(const BaseDimensions< T > &a, const BaseDimensions< Y > &b)
Definition: Dimensions.h:137
long long int64
Definition: SYS_Types.h:107
GA_API const UT_StringHolder scale
GLdouble n
Definition: glcorearb.h:2007
GLfloat f
Definition: glcorearb.h:1925
exint getNumRows() const
The number of rows in the matrix.
void removeIf(const CellFunction &condition)
int64 exint
Definition: SYS_Types.h:116
UT_SparseMatrixBuilderT(SparseMatrixType &matrix)
Construct a sparse matrix builder for the given matrix.
const UT_ValArray< T > & getRowValues() const
double fpreal64
Definition: SYS_Types.h:192
GLintptr offset
Definition: glcorearb.h:664
typename UT_TypePromoteT< T >::PreciseType UT_PreciseT
bool operator<(const GU_TetrahedronFacet &a, const GU_TetrahedronFacet &b)
#define SYS_ALIGN16
Definition: SYS_Align.h:100
UT_SparseMatrixBuilderT< fpreal32 > UT_SparseMatrixBuilderF
void visit(const CellFunction &inspect, bool serial=false)
bool appendRowElement(inttype row, inttype col, T val, int &rowidx)
png_bytepp row
Definition: png.h:1836
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1221
Triplet()
Construct a zero triplet.
bool operator>(const Tuple< SIZE, T0 > &t0, const Tuple< SIZE, T1 > &t1)
Definition: Tuple.h:219
void subtractMultVec(const UT_VectorT< T > &v, UT_VectorT< T > &result) const
Do result = result - (M * v)
#define THREADED_METHOD2_CONST(CLASSNAME, DOMULTI, METHOD, PARMTYPE1, PARMNAME1, PARMTYPE2, PARMNAME2)
#define THREADED_METHOD1_CONST(CLASSNAME, DOMULTI, METHOD, PARMTYPE1, PARMNAME1)
GLsizei const GLfloat * value
Definition: glcorearb.h:823
Triplet(exint r, exint c, T v)
Construct a triplet for the given row and column with value v.
exint getNumCellsInRow(exint row) const
The number of non-zero elements in the given row.
void accept(Visitor &visitor) const
UT_ValArray< inttype > & getColumns()
UT_SparseMatrixBuilderT< fpreal64 > UT_SparseMatrixBuilderD
int findCellFromRow(int row) const
GLuint index
Definition: glcorearb.h:785
GLint GLenum GLint x
Definition: glcorearb.h:408
GLuint GLfloat * val
Definition: glcorearb.h:1607
GLint GLint GLsizei GLint GLenum GLenum type
Definition: glcorearb.h:107
UT_ValArray< T > & getRowValues()
void swap(TfErrorTransport &l, TfErrorTransport &r)
#define UT_ASSERT(ZZ)
Definition: UT_Assert.h:126
bool operator>=(const TfRefPtr< T > &p, std::nullptr_t)
Definition: refPtr.h:1245
GLboolean r
Definition: glcorearb.h:1221
bool operator!=(const BaseDimensions< T > &a, const BaseDimensions< Y > &b)
Definition: Dimensions.h:165
#define const
Definition: zconf.h:214
UT_SparseMatrixCSRT< fpreal32 > UT_SparseMatrixCSRF
int getNonZerosPerRow() const
#define THREADED_METHOD3_CONST(CLASSNAME, DOMULTI, METHOD, PARMTYPE1, PARMNAME1, PARMTYPE2, PARMNAME2, PARMTYPE3, PARMNAME3)
UT_SparseMatrixBuilderT()
Construct an uninitialized sparse matrix builder.
Declare prior to use.
exint index(inttype row, int nz) const
float fpreal32
Definition: SYS_Types.h:191
void UTparallelFor(const Range &range, const Body &body, const int subscribe_ratio=2, const int min_grain_size=1)
#define THREADED_METHOD(CLASSNAME, DOMULTI, METHOD)
Quat< T > operator+(const Quat< T > &q1, const Quat< T > &q2)
Definition: ImathQuat.h:919
void UTserialFor(const Range &range, const Body &body)
bool shouldMultiThread() const
int getNumCols() const
UT_Array< UT_ExintArray > LevelSet
GLenum src
Definition: glcorearb.h:1792
UT_SparseMatrixCSRT< fpreal64 > UT_SparseMatrixCSRD