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