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