HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ConjGradient.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 
4 /// @file ConjGradient.h
5 /// @authors D.J. Hill, Peter Cucka
6 /// @brief Preconditioned conjugate gradient solver (solves @e Ax = @e b using
7 /// the conjugate gradient method with one of a selection of preconditioners)
8 
9 #ifndef OPENVDB_MATH_CONJGRADIENT_HAS_BEEN_INCLUDED
10 #define OPENVDB_MATH_CONJGRADIENT_HAS_BEEN_INCLUDED
11 
12 #include <openvdb/Exceptions.h>
13 #include <openvdb/Types.h>
14 #include <openvdb/util/logging.h>
16 #include "Math.h" // for Abs(), isZero(), Max(), Sqrt()
17 #include <tbb/parallel_for.h>
18 #include <tbb/parallel_reduce.h>
19 #include <algorithm> // for std::lower_bound()
20 #include <cassert>
21 #include <cmath> // for std::isfinite()
22 #include <limits>
23 #include <sstream>
24 #include <string>
25 
26 
27 namespace openvdb {
29 namespace OPENVDB_VERSION_NAME {
30 namespace math {
31 namespace pcg {
32 
33 using SizeType = Index32;
34 
35 using SizeRange = tbb::blocked_range<SizeType>;
36 
37 template<typename ValueType> class Vector;
38 
39 template<typename ValueType, SizeType STENCIL_SIZE> class SparseStencilMatrix;
40 
41 template<typename ValueType> class Preconditioner;
42 template<typename MatrixType> class JacobiPreconditioner;
43 template<typename MatrixType> class IncompleteCholeskyPreconditioner;
44 
45 /// Information about the state of a conjugate gradient solution
46 struct State {
47  bool success;
49  double relativeError;
50  double absoluteError;
51 };
52 
53 
54 /// Return default termination conditions for a conjugate gradient solver.
55 template<typename ValueType>
56 inline State
58 {
59  State s;
60  s.success = false;
61  s.iterations = 50;
62  s.relativeError = 1.0e-6;
63  s.absoluteError = std::numeric_limits<ValueType>::epsilon() * 100.0;
64  return s;
65 }
66 
67 
68 ////////////////////////////////////////
69 
70 
71 /// @brief Solve @e Ax = @e b via the preconditioned conjugate gradient method.
72 ///
73 /// @param A a symmetric, positive-definite, @e N x @e N matrix
74 /// @param b a vector of size @e N
75 /// @param x a vector of size @e N
76 /// @param preconditioner a Preconditioner matrix
77 /// @param termination termination conditions given as a State object with the following fields:
78 /// <dl>
79 /// <dt><i>success</i>
80 /// <dd>ignored
81 /// <dt><i>iterations</i>
82 /// <dd>the maximum number of iterations, with or without convergence
83 /// <dt><i>relativeError</i>
84 /// <dd>the relative error ||<i>b</i> &minus; <i>Ax</i>|| / ||<i>b</i>||
85 /// that denotes convergence
86 /// <dt><i>absoluteError</i>
87 /// <dd>the absolute error ||<i>b</i> &minus; <i>Ax</i>|| that denotes convergence
88 ///
89 /// @throw ArithmeticError if either @a x or @a b is not of the appropriate size.
90 template<typename PositiveDefMatrix>
91 inline State
92 solve(
93  const PositiveDefMatrix& A,
94  const Vector<typename PositiveDefMatrix::ValueType>& b,
95  Vector<typename PositiveDefMatrix::ValueType>& x,
96  Preconditioner<typename PositiveDefMatrix::ValueType>& preconditioner,
97  const State& termination = terminationDefaults<typename PositiveDefMatrix::ValueType>());
98 
99 
100 /// @brief Solve @e Ax = @e b via the preconditioned conjugate gradient method.
101 ///
102 /// @param A a symmetric, positive-definite, @e N x @e N matrix
103 /// @param b a vector of size @e N
104 /// @param x a vector of size @e N
105 /// @param preconditioner a Preconditioner matrix
106 /// @param termination termination conditions given as a State object with the following fields:
107 /// <dl>
108 /// <dt><i>success</i>
109 /// <dd>ignored
110 /// <dt><i>iterations</i>
111 /// <dd>the maximum number of iterations, with or without convergence
112 /// <dt><i>relativeError</i>
113 /// <dd>the relative error ||<i>b</i> &minus; <i>Ax</i>|| / ||<i>b</i>||
114 /// that denotes convergence
115 /// <dt><i>absoluteError</i>
116 /// <dd>the absolute error ||<i>b</i> &minus; <i>Ax</i>|| that denotes convergence
117 /// @param interrupter an object adhering to the util::NullInterrupter interface
118 /// with which computation can be interrupted
119 ///
120 /// @throw ArithmeticError if either @a x or @a b is not of the appropriate size.
121 /// @throw RuntimeError if the computation is interrupted.
122 template<typename PositiveDefMatrix, typename Interrupter>
123 inline State
124 solve(
125  const PositiveDefMatrix& A,
126  const Vector<typename PositiveDefMatrix::ValueType>& b,
127  Vector<typename PositiveDefMatrix::ValueType>& x,
128  Preconditioner<typename PositiveDefMatrix::ValueType>& preconditioner,
129  Interrupter& interrupter,
130  const State& termination = terminationDefaults<typename PositiveDefMatrix::ValueType>());
131 
132 
133 ////////////////////////////////////////
134 
135 
136 /// Lightweight, variable-length vector
137 template<typename T>
138 class Vector
139 {
140 public:
141  using ValueType = T;
143 
144  /// Construct an empty vector.
145  Vector(): mData(nullptr), mSize(0) {}
146  /// Construct a vector of @a n elements, with uninitialized values.
147  Vector(SizeType n): mData(new T[n]), mSize(n) {}
148  /// Construct a vector of @a n elements and initialize each element to the given value.
149  Vector(SizeType n, const ValueType& val): mData(new T[n]), mSize(n) { this->fill(val); }
150 
151  ~Vector() { mSize = 0; delete[] mData; mData = nullptr; }
152 
153  /// Deep copy the given vector.
154  Vector(const Vector&);
155  /// Deep copy the given vector.
156  Vector& operator=(const Vector&);
157 
158  /// Return the number of elements in this vector.
159  SizeType size() const { return mSize; }
160  /// Return @c true if this vector has no elements.
161  bool empty() const { return (mSize == 0); }
162 
163  /// @brief Reset this vector to have @a n elements, with uninitialized values.
164  /// @warning All of this vector's existing values will be lost.
165  void resize(SizeType n);
166 
167  /// Swap internal storage with another vector, which need not be the same size.
168  void swap(Vector& other) { std::swap(mData, other.mData); std::swap(mSize, other.mSize); }
169 
170  /// Set all elements of this vector to @a value.
171  void fill(const ValueType& value);
172 
173  //@{
174  /// @brief Multiply each element of this vector by @a s.
175  template<typename Scalar> void scale(const Scalar& s);
176  template<typename Scalar> Vector& operator*=(const Scalar& s) { this->scale(s); return *this; }
177  //@}
178 
179  /// Return the dot product of this vector with the given vector, which must be the same size.
180  ValueType dot(const Vector&) const;
181 
182  /// Return the infinity norm of this vector.
183  ValueType infNorm() const;
184  /// Return the L2 norm of this vector.
185  ValueType l2Norm() const { return Sqrt(this->dot(*this)); }
186 
187  /// Return @c true if every element of this vector has a finite value.
188  bool isFinite() const;
189 
190  /// @brief Return @c true if this vector is equivalent to the given vector
191  /// to within the specified tolerance.
192  template<typename OtherValueType>
193  bool eq(const Vector<OtherValueType>& other,
194  ValueType eps = Tolerance<ValueType>::value()) const;
195 
196  /// Return a string representation of this vector.
197  std::string str() const;
198 
199  //@{
200  /// @brief Return the value of this vector's ith element.
201  inline T& at(SizeType i) { return mData[i]; }
202  inline const T& at(SizeType i) const { return mData[i]; }
203  inline T& operator[](SizeType i) { return this->at(i); }
204  inline const T& operator[](SizeType i) const { return this->at(i); }
205  //@}
206 
207  //@{
208  /// @brief Return a pointer to this vector's elements.
209  inline T* data() { return mData; }
210  inline const T* data() const { return mData; }
211  inline const T* constData() const { return mData; }
212  //@}
213 
214 private:
215  // Functor for use with tbb::parallel_for()
216  template<typename Scalar> struct ScaleOp;
217  struct DeterministicDotProductOp;
218  // Functors for use with tbb::parallel_reduce()
219  template<typename OtherValueType> struct EqOp;
220  struct InfNormOp;
221  struct IsFiniteOp;
222 
223  T* mData;
224  SizeType mSize;
225 };
226 
229 
230 
231 ////////////////////////////////////////
232 
233 
234 /// @brief Sparse, square matrix representing a 3D stencil operator of size @a STENCIL_SIZE
235 /// @details The implementation is a variation on compressed row storage (CRS).
236 template<typename ValueType_, SizeType STENCIL_SIZE>
238 {
239 public:
240  using ValueType = ValueType_;
243 
244  class ConstValueIter;
245  class ConstRow;
246  class RowEditor;
247 
248  static const ValueType sZeroValue;
249 
250  /// Construct an @a n x @a n matrix with at most @a STENCIL_SIZE nonzero elements per row.
252 
253  /// Deep copy the given matrix.
255 
256  //@{
257  /// Return the number of rows in this matrix.
258  SizeType numRows() const { return mNumRows; }
259  SizeType size() const { return mNumRows; }
260  //@}
261 
262  /// @brief Set the value at the given coordinates.
263  /// @warning It is not safe to set values in the same row simultaneously
264  /// from multiple threads.
265  void setValue(SizeType row, SizeType col, const ValueType&);
266 
267  //@{
268  /// @brief Return the value at the given coordinates.
269  /// @warning It is not safe to get values from a row while another thread
270  /// is setting values in that row.
271  const ValueType& getValue(SizeType row, SizeType col) const;
272  const ValueType& operator()(SizeType row, SizeType col) const;
273  //@}
274 
275  /// Return a read-only view onto the given row of this matrix.
276  ConstRow getConstRow(SizeType row) const;
277 
278  /// Return a read/write view onto the given row of this matrix.
279  RowEditor getRowEditor(SizeType row);
280 
281  //@{
282  /// @brief Multiply all elements in the matrix by @a s;
283  template<typename Scalar> void scale(const Scalar& s);
284  template<typename Scalar>
285  SparseStencilMatrix& operator*=(const Scalar& s) { this->scale(s); return *this; }
286  //@}
287 
288  /// @brief Multiply this matrix by @a inVec and return the result in @a resultVec.
289  /// @throw ArithmeticError if either @a inVec or @a resultVec is not of size @e N,
290  /// where @e N x @e N is the size of this matrix.
291  template<typename VecValueType>
292  void vectorMultiply(const Vector<VecValueType>& inVec, Vector<VecValueType>& resultVec) const;
293 
294  /// @brief Multiply this matrix by the vector represented by the array @a inVec
295  /// and return the result in @a resultVec.
296  /// @warning Both @a inVec and @a resultVec must have at least @e N elements,
297  /// where @e N x @e N is the size of this matrix.
298  template<typename VecValueType>
299  void vectorMultiply(const VecValueType* inVec, VecValueType* resultVec) const;
300 
301  /// @brief Return @c true if this matrix is equivalent to the given matrix
302  /// to within the specified tolerance.
303  template<typename OtherValueType>
305  ValueType eps = Tolerance<ValueType>::value()) const;
306 
307  /// Return @c true if every element of this matrix has a finite value.
308  bool isFinite() const;
309 
310  /// Return a string representation of this matrix.
311  std::string str() const;
312 
313 private:
314  struct RowData {
315  RowData(ValueType* v, SizeType* c, SizeType& s): mVals(v), mCols(c), mSize(s) {}
316  ValueType* mVals; SizeType* mCols; SizeType& mSize;
317  };
318 
319  struct ConstRowData {
320  ConstRowData(const ValueType* v, const SizeType* c, const SizeType& s):
321  mVals(v), mCols(c), mSize(s) {}
322  const ValueType* mVals; const SizeType* mCols; const SizeType& mSize;
323  };
324 
325  /// Base class for row accessors
326  template<typename DataType_ = RowData>
327  class RowBase
328  {
329  public:
330  using DataType = DataType_;
331 
332  static SizeType capacity() { return STENCIL_SIZE; }
333 
334  RowBase(const DataType& data): mData(data) {}
335 
336  bool empty() const { return (mData.mSize == 0); }
337  const SizeType& size() const { return mData.mSize; }
338 
339  const ValueType& getValue(SizeType columnIdx, bool& active) const;
340  const ValueType& getValue(SizeType columnIdx) const;
341 
342  /// Return an iterator over the stored values in this row.
343  ConstValueIter cbegin() const;
344 
345  /// @brief Return @c true if this row is equivalent to the given row
346  /// to within the specified tolerance.
347  template<typename OtherDataType>
348  bool eq(const RowBase<OtherDataType>& other,
349  ValueType eps = Tolerance<ValueType>::value()) const;
350 
351  /// @brief Return the dot product of this row with the first
352  /// @a vecSize elements of @a inVec.
353  /// @warning @a inVec must have at least @a vecSize elements.
354  template<typename VecValueType>
355  VecValueType dot(const VecValueType* inVec, SizeType vecSize) const;
356 
357  /// Return the dot product of this row with the given vector.
358  template<typename VecValueType>
359  VecValueType dot(const Vector<VecValueType>& inVec) const;
360 
361  /// Return a string representation of this row.
362  std::string str() const;
363 
364  protected:
365  friend class ConstValueIter;
366 
367  const ValueType& value(SizeType i) const { return mData.mVals[i]; }
368  SizeType column(SizeType i) const { return mData.mCols[i]; }
369 
370  /// @brief Return the array index of the first column index that is
371  /// equal to <i>or greater than</i> the given column index.
372  /// @note If @a columnIdx is larger than any existing column index,
373  /// the return value will point beyond the end of the array.
374  SizeType find(SizeType columnIdx) const;
375 
376  DataType mData;
377  };
378 
379  using ConstRowBase = RowBase<ConstRowData>;
380 
381 public:
382  /// Iterator over the stored values in a row of this matrix
384  {
385  public:
386  const ValueType& operator*() const
387  {
388  if (mData.mSize == 0) return SparseStencilMatrix::sZeroValue;
389  return mData.mVals[mCursor];
390  }
391 
392  SizeType column() const { return mData.mCols[mCursor]; }
393 
394  void increment() { mCursor++; }
395  ConstValueIter& operator++() { increment(); return *this; }
396  operator bool() const { return (mCursor < mData.mSize); }
397 
398  void reset() { mCursor = 0; }
399 
400  private:
401  friend class SparseStencilMatrix;
402  ConstValueIter(const RowData& d): mData(d.mVals, d.mCols, d.mSize), mCursor(0) {}
403  ConstValueIter(const ConstRowData& d): mData(d), mCursor(0) {}
404 
405  const ConstRowData mData;
406  SizeType mCursor;
407  };
408 
409 
410  /// Read-only accessor to a row of this matrix
411  class ConstRow: public ConstRowBase
412  {
413  public:
414  ConstRow(const ValueType* valueHead, const SizeType* columnHead, const SizeType& rowSize);
415  }; // class ConstRow
416 
417 
418  /// Read/write accessor to a row of this matrix
419  class RowEditor: public RowBase<>
420  {
421  public:
422  RowEditor(ValueType* valueHead, SizeType* columnHead, SizeType& rowSize, SizeType colSize);
423 
424  /// Set the number of entries in this row to zero.
425  void clear();
426 
427  /// @brief Set the value of the entry in the specified column.
428  /// @return the current number of entries stored in this row.
430 
431  //@{
432  /// @brief Scale all of the entries in this row.
433  template<typename Scalar> void scale(const Scalar&);
434  template<typename Scalar>
435  RowEditor& operator*=(const Scalar& s) { this->scale(s); return *this; }
436  //@}
437 
438  private:
439  const SizeType mNumColumns; // used only for bounds checking
440  }; // class RowEditor
441 
442 private:
443  // Functors for use with tbb::parallel_for()
444  struct MatrixCopyOp;
445  template<typename VecValueType> struct VecMultOp;
446  template<typename Scalar> struct RowScaleOp;
447 
448  // Functors for use with tbb::parallel_reduce()
449  struct IsFiniteOp;
450  template<typename OtherValueType> struct EqOp;
451 
452  const SizeType mNumRows;
453  std::unique_ptr<ValueType[]> mValueArray;
454  std::unique_ptr<SizeType[]> mColumnIdxArray;
455  std::unique_ptr<SizeType[]> mRowSizeArray;
456 }; // class SparseStencilMatrix
457 
458 
459 ////////////////////////////////////////
460 
461 
462 /// Base class for conjugate gradient preconditioners
463 template<typename T>
464 class Preconditioner
465 {
466 public:
467  using ValueType = T;
469 
470  template<SizeType STENCIL_SIZE> Preconditioner(const SparseStencilMatrix<T, STENCIL_SIZE>&) {}
471  virtual ~Preconditioner() = default;
472 
473  virtual bool isValid() const { return true; }
474 
475  /// @brief Apply this preconditioner to a residue vector:
476  /// @e z = <i>M</i><sup><small>&minus;1</small></sup><i>r</i>
477  /// @param r residue vector
478  /// @param[out] z preconditioned residue vector
479  virtual void apply(const Vector<T>& r, Vector<T>& z) = 0;
480 };
481 
482 
483 ////////////////////////////////////////
484 
485 
486 namespace internal {
487 
488 // Functor for use with tbb::parallel_for() to copy data from one array to another
489 template<typename T>
490 struct CopyOp
491 {
492  CopyOp(const T* from_, T* to_): from(from_), to(to_) {}
493 
494  void operator()(const SizeRange& range) const {
495  for (SizeType n = range.begin(), N = range.end(); n < N; ++n) to[n] = from[n];
496  }
497 
498  const T* from;
499  T* to;
500 };
501 
502 
503 // Functor for use with tbb::parallel_for() to fill an array with a constant value
504 template<typename T>
505 struct FillOp
506 {
507  FillOp(T* data_, const T& val_): data(data_), val(val_) {}
508 
509  void operator()(const SizeRange& range) const {
510  for (SizeType n = range.begin(), N = range.end(); n < N; ++n) data[n] = val;
511  }
512 
513  T* data;
514  const T val;
515 };
516 
517 
518 // Functor for use with tbb::parallel_for() that computes a * x + y
519 template<typename T>
520 struct LinearOp
521 {
522  LinearOp(const T& a_, const T* x_, const T* y_, T* out_): a(a_), x(x_), y(y_), out(out_) {}
523 
524  void operator()(const SizeRange& range) const {
525  if (isExactlyEqual(a, T(1))) {
526  for (SizeType n = range.begin(), N = range.end(); n < N; ++n) out[n] = x[n] + y[n];
527  } else if (isExactlyEqual(a, T(-1))) {
528  for (SizeType n = range.begin(), N = range.end(); n < N; ++n) out[n] = -x[n] + y[n];
529  } else {
530  for (SizeType n = range.begin(), N = range.end(); n < N; ++n) out[n] = a * x[n] + y[n];
531  }
532  }
533 
534  const T a, *x, *y;
535  T* out;
536 };
537 
538 } // namespace internal
539 
540 
541 ////////////////////////////////////////
542 
543 
544 inline std::ostream&
545 operator<<(std::ostream& os, const State& state)
546 {
547  os << (state.success ? "succeeded with " : "")
548  << "rel. err. " << state.relativeError << ", abs. err. " << state.absoluteError
549  << " after " << state.iterations << " iteration" << (state.iterations == 1 ? "" : "s");
550  return os;
551 }
552 
553 
554 ////////////////////////////////////////
555 
556 
557 template<typename T>
558 inline
559 Vector<T>::Vector(const Vector& other): mData(new T[other.mSize]), mSize(other.mSize)
560 {
561  tbb::parallel_for(SizeRange(0, mSize),
562  internal::CopyOp<T>(/*from=*/other.mData, /*to=*/mData));
563 }
564 
565 
566 template<typename T>
567 inline
569 {
570  // Update the internal storage to the correct size
571 
572  if (mSize != other.mSize) {
573  mSize = other.mSize;
574  delete[] mData;
575  mData = new T[mSize];
576  }
577 
578  // Deep copy the data
579  tbb::parallel_for(SizeRange(0, mSize),
580  internal::CopyOp<T>(/*from=*/other.mData, /*to=*/mData));
581 
582  return *this;
583 }
584 
585 
586 template<typename T>
587 inline void
589 {
590  if (n != mSize) {
591  if (mData) delete[] mData;
592  mData = new T[n];
593  mSize = n;
594  }
595 }
596 
597 
598 template<typename T>
599 inline void
601 {
602  tbb::parallel_for(SizeRange(0, mSize), internal::FillOp<T>(mData, value));
603 }
604 
605 
606 template<typename T>
607 template<typename Scalar>
608 struct Vector<T>::ScaleOp
609 {
610  ScaleOp(T* data_, const Scalar& s_): data(data_), s(s_) {}
611 
612  void operator()(const SizeRange& range) const {
613  for (SizeType n = range.begin(), N = range.end(); n < N; ++n) data[n] *= s;
614  }
615 
616  T* data;
617  const Scalar s;
618 };
619 
620 
621 template<typename T>
622 template<typename Scalar>
623 inline void
624 Vector<T>::scale(const Scalar& s)
625 {
626  tbb::parallel_for(SizeRange(0, mSize), ScaleOp<Scalar>(mData, s));
627 }
628 
629 
630 template<typename T>
632 {
633  DeterministicDotProductOp(const T* a_, const T* b_,
634  const SizeType binCount_, const SizeType arraySize_, T* reducetmp_):
635  a(a_), b(b_), binCount(binCount_), arraySize(arraySize_), reducetmp(reducetmp_) {}
636 
637  void operator()(const SizeRange& range) const
638  {
639  const SizeType binSize = arraySize / binCount;
640 
641  // Iterate over bins (array segments)
642  for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
643  const SizeType begin = n * binSize;
644  const SizeType end = (n == binCount-1) ? arraySize : begin + binSize;
645 
646  // Compute the partial sum for this array segment
647  T sum = zeroVal<T>();
648  for (SizeType i = begin; i < end; ++i) { sum += a[i] * b[i]; }
649  // Store the partial sum
650  reducetmp[n] = sum;
651  }
652  }
653 
654 
655  const T* a;
656  const T* b;
660 };
661 
662 template<typename T>
663 inline T
664 Vector<T>::dot(const Vector<T>& other) const
665 {
666  assert(this->size() == other.size());
667 
668  const T* aData = this->data();
669  const T* bData = other.data();
670 
671  SizeType arraySize = this->size();
672 
673  T result = zeroVal<T>();
674 
675  if (arraySize < 1024) {
676 
677  // Compute the dot product in serial for small arrays
678 
679  for (SizeType n = 0; n < arraySize; ++n) {
680  result += aData[n] * bData[n];
681  }
682 
683  } else {
684 
685  // Compute the dot product by segmenting the arrays into
686  // a predetermined number of sub arrays in parallel and
687  // accumulate the finial result in series.
688 
689  const SizeType binCount = 100;
690  T partialSums[100];
691 
692  tbb::parallel_for(SizeRange(0, binCount),
693  DeterministicDotProductOp(aData, bData, binCount, arraySize, partialSums));
694 
695  for (SizeType n = 0; n < binCount; ++n) {
696  result += partialSums[n];
697  }
698  }
699 
700  return result;
701 }
702 
703 
704 template<typename T>
706 {
707  InfNormOp(const T* data_): data(data_) {}
708 
709  T operator()(const SizeRange& range, T maxValue) const
710  {
711  for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
712  maxValue = Max(maxValue, Abs(data[n]));
713  }
714  return maxValue;
715  }
716 
717  const T* data;
718 };
719 
720 
721 template<typename T>
722 inline T
724 {
725  // Parallelize over the elements of this vector.
726  T result = tbb::parallel_reduce(SizeRange(0, this->size()), /*seed=*/zeroVal<T>(),
727  InfNormOp(this->data()), /*join=*/[](T max1, T max2) { return Max(max1, max2); });
728  return result;
729 }
730 
731 
732 template<typename T>
734 {
735  IsFiniteOp(const T* data_): data(data_) {}
736 
737  bool operator()(const SizeRange& range, bool finite) const
738  {
739  if (finite) {
740  for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
741  if (!std::isfinite(data[n])) return false;
742  }
743  }
744  return finite;
745  }
746 
747  const T* data;
748 };
749 
750 
751 template<typename T>
752 inline bool
754 {
755  // Parallelize over the elements of this vector.
756  bool finite = tbb::parallel_reduce(SizeRange(0, this->size()), /*seed=*/true,
757  IsFiniteOp(this->data()),
758  /*join=*/[](bool finite1, bool finite2) { return (finite1 && finite2); });
759  return finite;
760 }
761 
762 
763 template<typename T>
764 template<typename OtherValueType>
765 struct Vector<T>::EqOp
766 {
767  EqOp(const T* a_, const OtherValueType* b_, T e): a(a_), b(b_), eps(e) {}
768 
769  bool operator()(const SizeRange& range, bool equal) const
770  {
771  if (equal) {
772  for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
773  if (!isApproxEqual(a[n], b[n], eps)) return false;
774  }
775  }
776  return equal;
777  }
778 
779  const T* a;
780  const OtherValueType* b;
781  const T eps;
782 };
783 
784 
785 template<typename T>
786 template<typename OtherValueType>
787 inline bool
789 {
790  if (this->size() != other.size()) return false;
791  bool equal = tbb::parallel_reduce(SizeRange(0, this->size()), /*seed=*/true,
792  EqOp<OtherValueType>(this->data(), other.data(), eps),
793  /*join=*/[](bool eq1, bool eq2) { return (eq1 && eq2); });
794  return equal;
795 }
796 
797 
798 template<typename T>
799 inline std::string
801 {
802  std::ostringstream ostr;
803  ostr << "[";
804  std::string sep;
805  for (SizeType n = 0, N = this->size(); n < N; ++n) {
806  ostr << sep << (*this)[n];
807  sep = ", ";
808  }
809  ostr << "]";
810  return ostr.str();
811 }
812 
813 
814 ////////////////////////////////////////
815 
816 
817 template<typename ValueType, SizeType STENCIL_SIZE>
818 const ValueType SparseStencilMatrix<ValueType, STENCIL_SIZE>::sZeroValue = zeroVal<ValueType>();
819 
820 
821 template<typename ValueType, SizeType STENCIL_SIZE>
822 inline
824  : mNumRows(numRows)
825  , mValueArray(new ValueType[mNumRows * STENCIL_SIZE])
826  , mColumnIdxArray(new SizeType[mNumRows * STENCIL_SIZE])
827  , mRowSizeArray(new SizeType[mNumRows])
828 {
829  // Initialize the matrix to a null state by setting the size of each row to zero.
830  tbb::parallel_for(SizeRange(0, mNumRows),
831  internal::FillOp<SizeType>(mRowSizeArray.get(), /*value=*/0));
832 }
833 
834 
835 template<typename ValueType, SizeType STENCIL_SIZE>
837 {
839  from(&from_), to(&to_) {}
840 
841  void operator()(const SizeRange& range) const
842  {
843  const ValueType* fromVal = from->mValueArray.get();
844  const SizeType* fromCol = from->mColumnIdxArray.get();
845  ValueType* toVal = to->mValueArray.get();
846  SizeType* toCol = to->mColumnIdxArray.get();
847  for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
848  toVal[n] = fromVal[n];
849  toCol[n] = fromCol[n];
850  }
851  }
852 
854 };
855 
856 
857 template<typename ValueType, SizeType STENCIL_SIZE>
858 inline
860  : mNumRows(other.mNumRows)
861  , mValueArray(new ValueType[mNumRows * STENCIL_SIZE])
862  , mColumnIdxArray(new SizeType[mNumRows * STENCIL_SIZE])
863  , mRowSizeArray(new SizeType[mNumRows])
864 {
865  SizeType size = mNumRows * STENCIL_SIZE;
866 
867  // Copy the value and column index arrays from the other matrix to this matrix.
868  tbb::parallel_for(SizeRange(0, size), MatrixCopyOp(/*from=*/other, /*to=*/*this));
869 
870  // Copy the row size array from the other matrix to this matrix.
871  tbb::parallel_for(SizeRange(0, mNumRows),
872  internal::CopyOp<SizeType>(/*from=*/other.mRowSizeArray.get(), /*to=*/mRowSizeArray.get()));
873 }
874 
875 
876 template<typename ValueType, SizeType STENCIL_SIZE>
877 inline void
879  const ValueType& val)
880 {
881  assert(row < mNumRows);
882  this->getRowEditor(row).setValue(col, val);
883 }
884 
885 
886 template<typename ValueType, SizeType STENCIL_SIZE>
887 inline const ValueType&
889 {
890  assert(row < mNumRows);
891  return this->getConstRow(row).getValue(col);
892 }
893 
894 
895 template<typename ValueType, SizeType STENCIL_SIZE>
896 inline const ValueType&
898 {
899  return this->getValue(row,col);
900 }
901 
902 
903 template<typename ValueType, SizeType STENCIL_SIZE>
904 template<typename Scalar>
905 struct SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowScaleOp
906 {
907  RowScaleOp(SparseStencilMatrix& m, const Scalar& s_): mat(&m), s(s_) {}
908 
909  void operator()(const SizeRange& range) const
910  {
911  for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
912  RowEditor row = mat->getRowEditor(n);
913  row.scale(s);
914  }
915  }
916 
917  SparseStencilMatrix* mat;
918  const Scalar s;
919 };
920 
921 
922 template<typename ValueType, SizeType STENCIL_SIZE>
923 template<typename Scalar>
924 inline void
926 {
927  // Parallelize over the rows in the matrix.
928  tbb::parallel_for(SizeRange(0, mNumRows), RowScaleOp<Scalar>(*this, s));
929 }
930 
931 
932 template<typename ValueType, SizeType STENCIL_SIZE>
933 template<typename VecValueType>
934 struct SparseStencilMatrix<ValueType, STENCIL_SIZE>::VecMultOp
935 {
936  VecMultOp(const SparseStencilMatrix& m, const VecValueType* i, VecValueType* o):
937  mat(&m), in(i), out(o) {}
938 
939  void operator()(const SizeRange& range) const
940  {
941  for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
942  ConstRow row = mat->getConstRow(n);
943  out[n] = row.dot(in, mat->numRows());
944  }
945  }
946 
947  const SparseStencilMatrix* mat;
948  const VecValueType* in;
949  VecValueType* out;
950 };
951 
952 
953 template<typename ValueType, SizeType STENCIL_SIZE>
954 template<typename VecValueType>
955 inline void
957  const Vector<VecValueType>& inVec, Vector<VecValueType>& resultVec) const
958 {
959  if (inVec.size() != mNumRows) {
960  OPENVDB_THROW(ArithmeticError, "matrix and input vector have incompatible sizes ("
961  << mNumRows << "x" << mNumRows << " vs. " << inVec.size() << ")");
962  }
963  if (resultVec.size() != mNumRows) {
964  OPENVDB_THROW(ArithmeticError, "matrix and result vector have incompatible sizes ("
965  << mNumRows << "x" << mNumRows << " vs. " << resultVec.size() << ")");
966  }
967 
968  vectorMultiply(inVec.data(), resultVec.data());
969 }
970 
971 
972 template<typename ValueType, SizeType STENCIL_SIZE>
973 template<typename VecValueType>
974 inline void
976  const VecValueType* inVec, VecValueType* resultVec) const
977 {
978  // Parallelize over the rows in the matrix.
979  tbb::parallel_for(SizeRange(0, mNumRows),
980  VecMultOp<VecValueType>(*this, inVec, resultVec));
981 }
982 
983 
984 template<typename ValueType, SizeType STENCIL_SIZE>
985 template<typename OtherValueType>
986 struct SparseStencilMatrix<ValueType, STENCIL_SIZE>::EqOp
987 {
988  EqOp(const SparseStencilMatrix& a_,
990  a(&a_), b(&b_), eps(e) {}
991 
992  bool operator()(const SizeRange& range, bool equal) const
993  {
994  if (equal) {
995  for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
996  if (!a->getConstRow(n).eq(b->getConstRow(n), eps)) return false;
997  }
998  }
999  return equal;
1000  }
1001 
1002  const SparseStencilMatrix* a;
1003  const SparseStencilMatrix<OtherValueType, STENCIL_SIZE>* b;
1004  const ValueType eps;
1005 };
1006 
1007 
1008 template<typename ValueType, SizeType STENCIL_SIZE>
1009 template<typename OtherValueType>
1010 inline bool
1013 {
1014  if (this->numRows() != other.numRows()) return false;
1015  bool equal = tbb::parallel_reduce(SizeRange(0, this->numRows()), /*seed=*/true,
1016  EqOp<OtherValueType>(*this, other, eps),
1017  /*join=*/[](bool eq1, bool eq2) { return (eq1 && eq2); });
1018  return equal;
1019 }
1020 
1021 
1022 template<typename ValueType, SizeType STENCIL_SIZE>
1024 {
1025  IsFiniteOp(const SparseStencilMatrix& m): mat(&m) {}
1026 
1027  bool operator()(const SizeRange& range, bool finite) const
1028  {
1029  if (finite) {
1030  for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1031  const ConstRow row = mat->getConstRow(n);
1032  for (ConstValueIter it = row.cbegin(); it; ++it) {
1033  if (!std::isfinite(*it)) return false;
1034  }
1035  }
1036  }
1037  return finite;
1038  }
1039 
1041 };
1042 
1043 
1044 template<typename ValueType, SizeType STENCIL_SIZE>
1045 inline bool
1047 {
1048  // Parallelize over the rows of this matrix.
1049  bool finite = tbb::parallel_reduce(SizeRange(0, this->numRows()), /*seed=*/true,
1050  IsFiniteOp(*this), /*join=*/[](bool finite1, bool finite2) { return (finite1&&finite2); });
1051  return finite;
1052 }
1053 
1054 
1055 template<typename ValueType, SizeType STENCIL_SIZE>
1056 inline std::string
1058 {
1059  std::ostringstream ostr;
1060  for (SizeType n = 0, N = this->size(); n < N; ++n) {
1061  ostr << n << ": " << this->getConstRow(n).str() << "\n";
1062  }
1063  return ostr.str();
1064 }
1065 
1066 
1067 template<typename ValueType, SizeType STENCIL_SIZE>
1070 {
1071  assert(i < mNumRows);
1072  const SizeType head = i * STENCIL_SIZE;
1073  return RowEditor(&mValueArray[head], &mColumnIdxArray[head], mRowSizeArray[i], mNumRows);
1074 }
1075 
1076 
1077 template<typename ValueType, SizeType STENCIL_SIZE>
1080 {
1081  assert(i < mNumRows);
1082  const SizeType head = i * STENCIL_SIZE; // index for this row into main storage
1083  return ConstRow(&mValueArray[head], &mColumnIdxArray[head], mRowSizeArray[i]);
1084 }
1085 
1086 
1087 template<typename ValueType, SizeType STENCIL_SIZE>
1088 template<typename DataType>
1089 inline SizeType
1091 {
1092  if (this->empty()) return mData.mSize;
1093 
1094  // Get a pointer to the first column index that is equal to or greater than the given index.
1095  // (This assumes that the data is sorted by column.)
1096  const SizeType* colPtr = std::lower_bound(mData.mCols, mData.mCols + mData.mSize, columnIdx);
1097  // Return the offset of the pointer from the beginning of the array.
1098  return static_cast<SizeType>(colPtr - mData.mCols);
1099 }
1100 
1101 
1102 template<typename ValueType, SizeType STENCIL_SIZE>
1103 template<typename DataType>
1104 inline const ValueType&
1106  SizeType columnIdx, bool& active) const
1107 {
1108  active = false;
1109  SizeType idx = this->find(columnIdx);
1110  if (idx < this->size() && this->column(idx) == columnIdx) {
1111  active = true;
1112  return this->value(idx);
1113  }
1115 }
1116 
1117 template<typename ValueType, SizeType STENCIL_SIZE>
1118 template<typename DataType>
1119 inline const ValueType&
1121 {
1122  SizeType idx = this->find(columnIdx);
1123  if (idx < this->size() && this->column(idx) == columnIdx) {
1124  return this->value(idx);
1125  }
1127 }
1128 
1129 
1130 template<typename ValueType, SizeType STENCIL_SIZE>
1131 template<typename DataType>
1132 inline typename SparseStencilMatrix<ValueType, STENCIL_SIZE>::ConstValueIter
1133 SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::cbegin() const
1134 {
1135  return ConstValueIter(mData);
1136 }
1137 
1138 
1139 template<typename ValueType, SizeType STENCIL_SIZE>
1140 template<typename DataType>
1141 template<typename OtherDataType>
1142 inline bool
1144  const RowBase<OtherDataType>& other, ValueType eps) const
1145 {
1146  if (this->size() != other.size()) return false;
1147  for (ConstValueIter it = cbegin(), oit = other.cbegin(); it || oit; ++it, ++oit) {
1148  if (it.column() != oit.column()) return false;
1149  if (!isApproxEqual(*it, *oit, eps)) return false;
1150  }
1151  return true;
1152 }
1153 
1154 
1155 template<typename ValueType, SizeType STENCIL_SIZE>
1156 template<typename DataType>
1157 template<typename VecValueType>
1158 inline VecValueType
1160  const VecValueType* inVec, SizeType vecSize) const
1161 {
1162  VecValueType result = zeroVal<VecValueType>();
1163  for (SizeType idx = 0, N = std::min(vecSize, this->size()); idx < N; ++idx) {
1164  result += static_cast<VecValueType>(this->value(idx) * inVec[this->column(idx)]);
1165  }
1166  return result;
1167 }
1168 
1169 template<typename ValueType, SizeType STENCIL_SIZE>
1170 template<typename DataType>
1171 template<typename VecValueType>
1172 inline VecValueType
1174  const Vector<VecValueType>& inVec) const
1175 {
1176  return dot(inVec.data(), inVec.size());
1177 }
1178 
1179 
1180 template<typename ValueType, SizeType STENCIL_SIZE>
1181 template<typename DataType>
1182 inline std::string
1184 {
1185  std::ostringstream ostr;
1186  std::string sep;
1187  for (SizeType n = 0, N = this->size(); n < N; ++n) {
1188  ostr << sep << "(" << this->column(n) << ", " << this->value(n) << ")";
1189  sep = ", ";
1190  }
1191  return ostr.str();
1192 }
1193 
1194 
1195 template<typename ValueType, SizeType STENCIL_SIZE>
1196 inline
1198  const ValueType* valueHead, const SizeType* columnHead, const SizeType& rowSize):
1199  ConstRowBase(ConstRowData(const_cast<ValueType*>(valueHead),
1200  const_cast<SizeType*>(columnHead), const_cast<SizeType&>(rowSize)))
1201 {
1202 }
1203 
1204 
1205 template<typename ValueType, SizeType STENCIL_SIZE>
1206 inline
1208  ValueType* valueHead, SizeType* columnHead, SizeType& rowSize, SizeType colSize):
1209  RowBase<>(RowData(valueHead, columnHead, rowSize)), mNumColumns(colSize)
1210 {
1211 }
1212 
1213 
1214 template<typename ValueType, SizeType STENCIL_SIZE>
1215 inline void
1217 {
1218  // Note: since mSize is a reference, this modifies the underlying matrix.
1219  RowBase<>::mData.mSize = 0;
1220 }
1221 
1222 
1223 template<typename ValueType, SizeType STENCIL_SIZE>
1224 inline SizeType
1226  SizeType column, const ValueType& value)
1227 {
1228  assert(column < mNumColumns);
1229 
1230  RowData& data = RowBase<>::mData;
1231 
1232  // Get the offset of the first column index that is equal to or greater than
1233  // the column to be modified.
1234  SizeType offset = this->find(column);
1235 
1236  if (offset < data.mSize && data.mCols[offset] == column) {
1237  // If the column already exists, just update its value.
1238  data.mVals[offset] = value;
1239  return data.mSize;
1240  }
1241 
1242  // Check that it is safe to add a new column.
1243  assert(data.mSize < this->capacity());
1244 
1245  if (offset >= data.mSize) {
1246  // The new column's index is larger than any existing index. Append the new column.
1247  data.mVals[data.mSize] = value;
1248  data.mCols[data.mSize] = column;
1249  } else {
1250  // Insert the new column at the computed offset after shifting subsequent columns.
1251  for (SizeType i = data.mSize; i > offset; --i) {
1252  data.mVals[i] = data.mVals[i - 1];
1253  data.mCols[i] = data.mCols[i - 1];
1254  }
1255  data.mVals[offset] = value;
1256  data.mCols[offset] = column;
1257  }
1258  ++data.mSize;
1259 
1260  return data.mSize;
1261 }
1262 
1263 
1264 template<typename ValueType, SizeType STENCIL_SIZE>
1265 template<typename Scalar>
1266 inline void
1268 {
1269  for (int idx = 0, N = this->size(); idx < N; ++idx) {
1270  RowBase<>::mData.mVals[idx] *= s;
1271  }
1272 }
1273 
1274 
1275 ////////////////////////////////////////
1276 
1277 
1278 /// Diagonal preconditioner
1279 template<typename MatrixType>
1280 class JacobiPreconditioner: public Preconditioner<typename MatrixType::ValueType>
1281 {
1282 private:
1283  struct InitOp;
1284  struct ApplyOp;
1285 
1286 public:
1287  using ValueType = typename MatrixType::ValueType;
1291 
1292  JacobiPreconditioner(const MatrixType& A): BaseType(A), mDiag(A.numRows())
1293  {
1294  // Initialize vector mDiag with the values from the matrix diagonal.
1295  tbb::parallel_for(SizeRange(0, A.numRows()), InitOp(A, mDiag.data()));
1296  }
1297 
1298  ~JacobiPreconditioner() override = default;
1299 
1300  void apply(const Vector<ValueType>& r, Vector<ValueType>& z) override
1301  {
1302  const SizeType size = mDiag.size();
1303 
1304  assert(r.size() == z.size());
1305  assert(r.size() == size);
1306 
1307  tbb::parallel_for(SizeRange(0, size), ApplyOp(mDiag.data(), r.data(), z.data()));
1308  }
1309 
1310  /// Return @c true if all values along the diagonal are finite.
1311  bool isFinite() const { return mDiag.isFinite(); }
1312 
1313 private:
1314  // Functor for use with tbb::parallel_for()
1315  struct InitOp
1316  {
1317  InitOp(const MatrixType& m, ValueType* v): mat(&m), vec(v) {}
1318  void operator()(const SizeRange& range) const {
1319  for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1320  const ValueType val = mat->getValue(n, n);
1321  assert(!isApproxZero(val, ValueType(0.0001)));
1322  vec[n] = static_cast<ValueType>(1.0 / val);
1323  }
1324  }
1325  const MatrixType* mat; ValueType* vec;
1326  };
1327 
1328  // Functor for use with tbb::parallel_reduce()
1329  struct ApplyOp
1330  {
1331  ApplyOp(const ValueType* x_, const ValueType* y_, ValueType* out_):
1332  x(x_), y(y_), out(out_) {}
1333  void operator()(const SizeRange& range) const {
1334  for (SizeType n = range.begin(), N = range.end(); n < N; ++n) out[n] = x[n] * y[n];
1335  }
1336  const ValueType *x, *y; ValueType* out;
1337  };
1338 
1339  // The Jacobi preconditioner is a diagonal matrix
1340  VectorType mDiag;
1341 }; // class JacobiPreconditioner
1342 
1343 
1344 /// Preconditioner using incomplete Cholesky factorization
1345 template<typename MatrixType>
1346 class IncompleteCholeskyPreconditioner: public Preconditioner<typename MatrixType::ValueType>
1347 {
1348 private:
1349  struct CopyToLowerOp;
1350  struct TransposeOp;
1351 
1352 public:
1353  using ValueType = typename MatrixType::ValueType;
1358  using TriangleConstRow = typename TriangularMatrix::ConstRow;
1359  using TriangleRowEditor = typename TriangularMatrix::RowEditor;
1360 
1362  : BaseType(matrix)
1363  , mLowerTriangular(matrix.numRows())
1364  , mUpperTriangular(matrix.numRows())
1365  , mTempVec(matrix.numRows())
1366  {
1367  // Size of matrix
1368  const SizeType numRows = mLowerTriangular.numRows();
1369 
1370  // Copy the upper triangular part to the lower triangular part.
1371  tbb::parallel_for(SizeRange(0, numRows), CopyToLowerOp(matrix, mLowerTriangular));
1372 
1373  // Build the Incomplete Cholesky Matrix
1374  //
1375  // Algorithm:
1376  //
1377  // for (k = 0; k < size; ++k) {
1378  // A(k,k) = sqrt(A(k,k));
1379  // for (i = k +1, i < size; ++i) {
1380  // if (A(i,k) == 0) continue;
1381  // A(i,k) = A(i,k) / A(k,k);
1382  // }
1383  // for (j = k+1; j < size; ++j) {
1384  // for (i = j; i < size; ++i) {
1385  // if (A(i,j) == 0) continue;
1386  // A(i,j) -= A(i,k)*A(j,k);
1387  // }
1388  // }
1389  // }
1390 
1391  mPassedCompatibilityCondition = true;
1392 
1393  for (SizeType k = 0; k < numRows; ++k) {
1394 
1395  TriangleConstRow crow_k = mLowerTriangular.getConstRow(k);
1396  ValueType diagonalValue = crow_k.getValue(k);
1397 
1398  // Test if the matrix build has failed.
1399  if (diagonalValue < 1.e-5) {
1400  mPassedCompatibilityCondition = false;
1401  break;
1402  }
1403 
1404  diagonalValue = Sqrt(diagonalValue);
1405 
1406  TriangleRowEditor row_k = mLowerTriangular.getRowEditor(k);
1407  row_k.setValue(k, diagonalValue);
1408 
1409  // Exploit the fact that the matrix is symmetric.
1410  typename MatrixType::ConstRow srcRow = matrix.getConstRow(k);
1411  typename MatrixType::ConstValueIter citer = srcRow.cbegin();
1412  for ( ; citer; ++citer) {
1413  SizeType ii = citer.column();
1414  if (ii < k+1) continue; // look above diagonal
1415 
1416  TriangleRowEditor row_ii = mLowerTriangular.getRowEditor(ii);
1417 
1418  row_ii.setValue(k, *citer / diagonalValue);
1419  }
1420 
1421  // for (j = k+1; j < size; ++j) replaced by row iter below
1422  citer.reset(); // k,j entries
1423  for ( ; citer; ++citer) {
1424  SizeType j = citer.column();
1425  if (j < k+1) continue;
1426 
1427  TriangleConstRow row_j = mLowerTriangular.getConstRow(j);
1428  ValueType a_jk = row_j.getValue(k); // a_jk is non zero if a_kj is non zero
1429 
1430  // Entry (i,j) is non-zero if matrix(j,i) is nonzero
1431 
1432  typename MatrixType::ConstRow mask = matrix.getConstRow(j);
1433  typename MatrixType::ConstValueIter maskIter = mask.cbegin();
1434  for ( ; maskIter; ++maskIter) {
1435  SizeType i = maskIter.column();
1436  if (i < j) continue;
1437 
1438  TriangleConstRow crow_i = mLowerTriangular.getConstRow(i);
1439  ValueType a_ij = crow_i.getValue(j);
1440  ValueType a_ik = crow_i.getValue(k);
1441  TriangleRowEditor row_i = mLowerTriangular.getRowEditor(i);
1442  a_ij -= a_ik * a_jk;
1443 
1444  row_i.setValue(j, a_ij);
1445  }
1446  }
1447  }
1448 
1449  // Build the transpose of the IC matrix: mUpperTriangular
1450  tbb::parallel_for(SizeRange(0, numRows),
1451  TransposeOp(matrix, mLowerTriangular, mUpperTriangular));
1452  }
1453 
1454  ~IncompleteCholeskyPreconditioner() override = default;
1455 
1456  bool isValid() const override { return mPassedCompatibilityCondition; }
1457 
1458  void apply(const Vector<ValueType>& rVec, Vector<ValueType>& zVec) override
1459  {
1460  if (!mPassedCompatibilityCondition) {
1461  OPENVDB_THROW(ArithmeticError, "invalid Cholesky decomposition");
1462  }
1463 
1464  // Solve mUpperTriangular * mLowerTriangular * rVec = zVec;
1465 
1466  SizeType size = mLowerTriangular.numRows();
1467 
1468  zVec.fill(zeroVal<ValueType>());
1469  ValueType* zData = zVec.data();
1470 
1471  if (size == 0) return;
1472 
1473  assert(rVec.size() == size);
1474  assert(zVec.size() == size);
1475 
1476  // Allocate a temp vector
1477  mTempVec.fill(zeroVal<ValueType>());
1478  ValueType* tmpData = mTempVec.data();
1479  const ValueType* rData = rVec.data();
1480 
1481  // Solve mLowerTriangular * tmp = rVec;
1482  for (SizeType i = 0; i < size; ++i) {
1483  typename TriangularMatrix::ConstRow row = mLowerTriangular.getConstRow(i);
1484  ValueType diagonal = row.getValue(i);
1485  ValueType dot = row.dot(mTempVec);
1486  tmpData[i] = (rData[i] - dot) / diagonal;
1487  if (!std::isfinite(tmpData[i])) {
1488  OPENVDB_LOG_DEBUG_RUNTIME("1 diagonal was " << diagonal);
1489  OPENVDB_LOG_DEBUG_RUNTIME("1a diagonal " << row.getValue(i));
1490  }
1491  }
1492 
1493  // Solve mUpperTriangular * zVec = tmp;
1494  for (SizeType ii = 0; ii < size; ++ii) {
1495  SizeType i = size - 1 - ii;
1496  typename TriangularMatrix::ConstRow row = mUpperTriangular.getConstRow(i);
1497  ValueType diagonal = row.getValue(i);
1498  ValueType dot = row.dot(zVec);
1499  zData[i] = (tmpData[i] - dot) / diagonal;
1500  if (!std::isfinite(zData[i])) {
1501  OPENVDB_LOG_DEBUG_RUNTIME("2 diagonal was " << diagonal);
1502  }
1503  }
1504  }
1505 
1506  const TriangularMatrix& lowerMatrix() const { return mLowerTriangular; }
1507  const TriangularMatrix& upperMatrix() const { return mUpperTriangular; }
1508 
1509 private:
1510  // Functor for use with tbb::parallel_for()
1511  struct CopyToLowerOp
1512  {
1513  CopyToLowerOp(const MatrixType& m, TriangularMatrix& l): mat(&m), lower(&l) {}
1514  void operator()(const SizeRange& range) const {
1515  for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1516  typename TriangularMatrix::RowEditor outRow = lower->getRowEditor(n);
1517  outRow.clear();
1518  typename MatrixType::ConstRow inRow = mat->getConstRow(n);
1519  for (typename MatrixType::ConstValueIter it = inRow.cbegin(); it; ++it) {
1520  if (it.column() > n) continue; // skip above diagonal
1521  outRow.setValue(it.column(), *it);
1522  }
1523  }
1524  }
1525  const MatrixType* mat; TriangularMatrix* lower;
1526  };
1527 
1528  // Functor for use with tbb::parallel_for()
1529  struct TransposeOp
1530  {
1531  TransposeOp(const MatrixType& m, const TriangularMatrix& l, TriangularMatrix& u):
1532  mat(&m), lower(&l), upper(&u) {}
1533  void operator()(const SizeRange& range) const {
1534  for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1535  typename TriangularMatrix::RowEditor outRow = upper->getRowEditor(n);
1536  outRow.clear();
1537  // Use the fact that matrix is symmetric.
1538  typename MatrixType::ConstRow inRow = mat->getConstRow(n);
1539  for (typename MatrixType::ConstValueIter it = inRow.cbegin(); it; ++it) {
1540  const SizeType column = it.column();
1541  if (column < n) continue; // only set upper triangle
1542  outRow.setValue(column, lower->getValue(column, n));
1543  }
1544  }
1545  }
1546  const MatrixType* mat; const TriangularMatrix* lower; TriangularMatrix* upper;
1547  };
1548 
1549  TriangularMatrix mLowerTriangular;
1550  TriangularMatrix mUpperTriangular;
1551  Vector<ValueType> mTempVec;
1552  bool mPassedCompatibilityCondition;
1553 }; // class IncompleteCholeskyPreconditioner
1554 
1555 
1556 ////////////////////////////////////////
1557 
1558 
1559 namespace internal {
1560 
1561 /// Compute @e ax + @e y.
1562 template<typename T>
1563 inline void
1564 axpy(const T& a, const T* xVec, const T* yVec, T* resultVec, SizeType size)
1565 {
1566  tbb::parallel_for(SizeRange(0, size), LinearOp<T>(a, xVec, yVec, resultVec));
1567 }
1568 
1569 /// Compute @e ax + @e y.
1570 template<typename T>
1571 inline void
1572 axpy(const T& a, const Vector<T>& xVec, const Vector<T>& yVec, Vector<T>& result)
1573 {
1574  assert(xVec.size() == yVec.size());
1575  assert(xVec.size() == result.size());
1576  axpy(a, xVec.data(), yVec.data(), result.data(), xVec.size());
1577 }
1578 
1579 
1580 /// Compute @e r = @e b &minus; @e Ax.
1581 template<typename MatrixOperator, typename VecValueType>
1582 inline void
1583 computeResidual(const MatrixOperator& A, const VecValueType* x,
1584  const VecValueType* b, VecValueType* r)
1585 {
1586  // Compute r = A * x.
1587  A.vectorMultiply(x, r);
1588  // Compute r = b - r.
1589  tbb::parallel_for(SizeRange(0, A.numRows()), LinearOp<VecValueType>(-1.0, r, b, r));
1590 }
1591 
1592 /// Compute @e r = @e b &minus; @e Ax.
1593 template<typename MatrixOperator, typename T>
1594 inline void
1595 computeResidual(const MatrixOperator& A, const Vector<T>& x, const Vector<T>& b, Vector<T>& r)
1596 {
1597  assert(x.size() == b.size());
1598  assert(x.size() == r.size());
1599  assert(x.size() == A.numRows());
1600 
1601  computeResidual(A, x.data(), b.data(), r.data());
1602 }
1603 
1604 } // namespace internal
1605 
1606 
1607 ////////////////////////////////////////
1608 
1609 
1610 template<typename PositiveDefMatrix>
1611 inline State
1613  const PositiveDefMatrix& Amat,
1617  const State& termination)
1618 {
1619  util::NullInterrupter interrupter;
1620  return solve(Amat, bVec, xVec, precond, interrupter, termination);
1621 }
1622 
1623 
1624 template<typename PositiveDefMatrix, typename Interrupter>
1625 inline State
1627  const PositiveDefMatrix& Amat,
1631  Interrupter& interrupter,
1632  const State& termination)
1633 {
1634  using ValueType = typename PositiveDefMatrix::ValueType;
1635  using VectorType = Vector<ValueType>;
1636 
1637  State result;
1638  result.success = false;
1639  result.iterations = 0;
1640  result.relativeError = 0.0;
1641  result.absoluteError = 0.0;
1642 
1643  const SizeType size = Amat.numRows();
1644  if (size == 0) {
1645  OPENVDB_LOG_WARN("pcg::solve(): matrix has dimension zero");
1646  return result;
1647  }
1648  if (size != bVec.size()) {
1649  OPENVDB_THROW(ArithmeticError, "A and b have incompatible sizes"
1650  << size << "x" << size << " vs. " << bVec.size() << ")");
1651  }
1652  if (size != xVec.size()) {
1653  OPENVDB_THROW(ArithmeticError, "A and x have incompatible sizes"
1654  << size << "x" << size << " vs. " << xVec.size() << ")");
1655  }
1656 
1657  // Temp vectors
1658  VectorType zVec(size); // transformed residual (M^-1 r)
1659  VectorType pVec(size); // search direction
1660  VectorType qVec(size); // A * p
1661 
1662  // Compute norm of B (the source)
1663  const ValueType tmp = bVec.infNorm();
1664  const ValueType infNormOfB = isZero(tmp) ? ValueType(1) : tmp;
1665 
1666  // Compute rVec: residual = b - Ax.
1667  VectorType rVec(size); // vector of residuals
1668 
1669  internal::computeResidual(Amat, xVec, bVec, rVec);
1670 
1671  assert(rVec.isFinite());
1672 
1673  // Normalize the residual norm with the source norm and look for early out.
1674  result.absoluteError = static_cast<double>(rVec.infNorm());
1675  result.relativeError = static_cast<double>(result.absoluteError / infNormOfB);
1676  if (result.relativeError <= termination.relativeError) {
1677  result.success = true;
1678  return result;
1679  }
1680 
1681  // Iterations of the CG solve
1682 
1683  ValueType rDotZPrev(1); // inner product of <z,r>
1684 
1685  // Keep track of the minimum error to monitor convergence.
1687  ValueType l2Error;
1688 
1689  int iteration = 0;
1690  for ( ; iteration < termination.iterations; ++iteration) {
1691 
1692  if (interrupter.wasInterrupted()) {
1693  OPENVDB_THROW(RuntimeError, "conjugate gradient solver was interrupted");
1694  }
1695 
1696  OPENVDB_LOG_DEBUG_RUNTIME("pcg::solve() " << result);
1697 
1698  result.iterations = iteration + 1;
1699 
1700  // Apply preconditioner to residual
1701  // z_{k} = M^-1 r_{k}
1702  precond.apply(rVec, zVec);
1703 
1704  // <r,z>
1705  const ValueType rDotZ = rVec.dot(zVec);
1706  assert(std::isfinite(rDotZ));
1707 
1708  if (0 == iteration) {
1709  // Initialize
1710  pVec = zVec;
1711  } else {
1712  const ValueType beta = rDotZ / rDotZPrev;
1713  // p = beta * p + z
1714  internal::axpy(beta, pVec, zVec, /*result */pVec);
1715  }
1716 
1717  // q_{k} = A p_{k}
1718  Amat.vectorMultiply(pVec, qVec);
1719 
1720  // alpha = <r_{k-1}, z_{k-1}> / <p_{k},q_{k}>
1721  const ValueType pAp = pVec.dot(qVec);
1722  assert(std::isfinite(pAp));
1723 
1724  const ValueType alpha = rDotZ / pAp;
1725  rDotZPrev = rDotZ;
1726 
1727  // x_{k} = x_{k-1} + alpha * p_{k}
1728  internal::axpy(alpha, pVec, xVec, /*result=*/xVec);
1729 
1730  // r_{k} = r_{k-1} - alpha_{k-1} A p_{k}
1731  internal::axpy(-alpha, qVec, rVec, /*result=*/rVec);
1732 
1733  // update tolerances
1734  l2Error = rVec.l2Norm();
1735  minL2Error = Min(l2Error, minL2Error);
1736 
1737  result.absoluteError = static_cast<double>(rVec.infNorm());
1738  result.relativeError = static_cast<double>(result.absoluteError / infNormOfB);
1739 
1740  if (l2Error > 2 * minL2Error) {
1741  // The solution started to diverge.
1742  result.success = false;
1743  break;
1744  }
1745  if (!std::isfinite(result.absoluteError)) {
1746  // Total divergence of solution
1747  result.success = false;
1748  break;
1749  }
1750  if (result.absoluteError <= termination.absoluteError) {
1751  // Convergence
1752  result.success = true;
1753  break;
1754  }
1755  if (result.relativeError <= termination.relativeError) {
1756  // Convergence
1757  result.success = true;
1758  break;
1759  }
1760  }
1761  OPENVDB_LOG_DEBUG_RUNTIME("pcg::solve() " << result);
1762 
1763  return result;
1764 }
1765 
1766 } // namespace pcg
1767 } // namespace math
1768 } // namespace OPENVDB_VERSION_NAME
1769 } // namespace openvdb
1770 
1771 #endif // OPENVDB_MATH_CONJGRADIENT_HAS_BEEN_INCLUDED
bool isFinite() const
Return true if every element of this matrix has a finite value.
void apply(const Vector< ValueType > &rVec, Vector< ValueType > &zVec) override
void parallel_for(int64_t start, int64_t end, std::function< void(int64_t index)> &&task, parallel_options opt=parallel_options(0, Split_Y, 1))
Definition: parallel.h:127
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1221
std::string upper(string_view a)
Return an all-upper case version of a (locale-independent).
Definition: strutil.h:349
void scale(const Scalar &s)
Multiply each element of this vector by s.
Definition: ConjGradient.h:624
RowEditor getRowEditor(SizeType row)
Return a read/write view onto the given row of this matrix.
void scale(const Scalar &)
Scale all of the entries in this row.
Lightweight, variable-length vector.
Definition: ConjGradient.h:37
void resize(SizeType n)
Reset this vector to have n elements, with uninitialized values.
Definition: ConjGradient.h:588
T * data()
Return a pointer to this vector's elements.
Definition: ConjGradient.h:209
const ValueType & operator()(SizeType row, SizeType col) const
Return the value at the given coordinates.
Definition: ConjGradient.h:897
#define OPENVDB_LOG_WARN(mesg)
Definition: logging.h:277
void axpy(const T &a, const T *xVec, const T *yVec, T *resultVec, SizeType size)
Compute ax + y.
bool eq(const Vector< OtherValueType > &other, ValueType eps=Tolerance< ValueType >::value()) const
Return true if this vector is equivalent to the given vector to within the specified tolerance...
Definition: ConjGradient.h:788
const ValueType & getValue(SizeType row, SizeType col) const
Return the value at the given coordinates.
Definition: ConjGradient.h:888
bool eq(const SparseStencilMatrix< OtherValueType, STENCIL_SIZE > &other, ValueType eps=Tolerance< ValueType >::value()) const
Return true if this matrix is equivalent to the given matrix to within the specified tolerance...
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:444
void apply(const Vector< ValueType > &r, Vector< ValueType > &z) override
void computeResidual(const MatrixOperator &A, const Vector< T > &x, const Vector< T > &b, Vector< T > &r)
Compute r = b − Ax.
void swap(UT::ArraySet< Key, MULTI, MAX_LOAD_FACTOR_256, Clearer, Hash, KeyEqual > &a, UT::ArraySet< Key, MULTI, MAX_LOAD_FACTOR_256, Clearer, Hash, KeyEqual > &b)
Definition: UT_ArraySet.h:1629
void axpy(const T &a, const Vector< T > &xVec, const Vector< T > &yVec, Vector< T > &result)
Compute ax + y.
ConstRow(const ValueType *valueHead, const SizeType *columnHead, const SizeType &rowSize)
const GLfloat * c
Definition: glew.h:16631
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:178
**But if you need a or simply need to know when the task has note that the like this
Definition: thread.h:623
Iterator over the stored values in a row of this matrix.
Definition: ConjGradient.h:383
Information about the state of a conjugate gradient solution.
Definition: ConjGradient.h:46
SizeType size() const
Return the number of rows in this matrix.
Definition: ConjGradient.h:259
Dummy NOOP interrupter class defining interface.
std::string str() const
Return a string representation of this vector.
Definition: ConjGradient.h:800
Tolerance for floating-point comparison.
Definition: Math.h:147
SizeType numRows() const
Return the number of rows in this matrix.
Definition: ConjGradient.h:258
const T * data() const
Return a pointer to this vector's elements.
Definition: ConjGradient.h:210
const T & at(SizeType i) const
Return the value of this vector's ith element.
Definition: ConjGradient.h:202
GLuint in
Definition: glew.h:11552
ValueType l2Norm() const
Return the L2 norm of this vector.
Definition: ConjGradient.h:185
GLdouble l
Definition: glew.h:9164
DeterministicDotProductOp(const T *a_, const T *b_, const SizeType binCount_, const SizeType arraySize_, T *reducetmp_)
Definition: ConjGradient.h:633
Sparse, square matrix representing a 3D stencil operator of size STENCIL_SIZE.
Definition: ConjGradient.h:39
GLint GLenum GLint x
Definition: glcorearb.h:408
tbb::blocked_range< SizeType > SizeRange
Definition: ConjGradient.h:35
GLsizeiptr size
Definition: glcorearb.h:663
std::shared_ptr< T > SharedPtr
Definition: Types.h:110
T & operator[](SizeType i)
Return the value of this vector's ith element.
Definition: ConjGradient.h:203
SYS_FORCE_INLINE const_iterator end() const
bool operator()(const SizeRange &range, bool finite) const
T & at(SizeType i)
Return the value of this vector's ith element.
Definition: ConjGradient.h:201
GLuint64EXT * result
Definition: glew.h:14311
bool isApproxEqual(const Type &a, const Type &b, const Type &tolerance)
Return true if a is equal to b to within the given tolerance.
Definition: Math.h:407
ImageBuf OIIO_API min(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
std::ostream & operator<<(std::ostream &os, const State &state)
Definition: ConjGradient.h:545
Preconditioner using incomplete Cholesky factorization.
Definition: ConjGradient.h:43
Coord Abs(const Coord &xyz)
Definition: Coord.h:514
ValueType dot(const Vector &) const
Return the dot product of this vector with the given vector, which must be the same size...
Definition: ConjGradient.h:664
SparseStencilMatrix(SizeType n)
Construct an n x n matrix with at most STENCIL_SIZE nonzero elements per row.
Definition: ConjGradient.h:823
bool empty() const
Return true if this vector has no elements.
Definition: ConjGradient.h:161
State solve(const PositiveDefMatrix &A, const Vector< typename PositiveDefMatrix::ValueType > &b, Vector< typename PositiveDefMatrix::ValueType > &x, Preconditioner< typename PositiveDefMatrix::ValueType > &preconditioner, const State &termination=terminationDefaults< typename PositiveDefMatrix::ValueType >())
Solve Ax = b via the preconditioned conjugate gradient method.
GLint GLuint mask
Definition: glcorearb.h:123
FMT_CONSTEXPR bool find(Ptr first, Ptr last, T value, Ptr &out)
Definition: format.h:2929
MatrixCopyOp(const SparseStencilMatrix &from_, SparseStencilMatrix &to_)
Definition: ConjGradient.h:838
Vector(SizeType n)
Construct a vector of n elements, with uninitialized values.
Definition: ConjGradient.h:147
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
const GLdouble * v
Definition: glcorearb.h:836
void clear()
Set the number of entries in this row to zero.
fpreal64 dot(const CE_VectorT< T > &a, const CE_VectorT< T > &b)
Definition: CE_Vector.h:130
Base class for conjugate gradient preconditioners.
Definition: ConjGradient.h:41
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1221
GLuint GLuint end
Definition: glcorearb.h:474
GLsizei const GLchar *const * string
Definition: glcorearb.h:813
const T * constData() const
Return a pointer to this vector's elements.
Definition: ConjGradient.h:211
void computeResidual(const MatrixOperator &A, const VecValueType *x, const VecValueType *b, VecValueType *r)
Compute r = b − Ax.
ImageBuf OIIO_API max(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:764
GLdouble GLdouble GLdouble z
Definition: glcorearb.h:847
Vector & operator*=(const Scalar &s)
Multiply each element of this vector by s.
Definition: ConjGradient.h:176
void swap(Vector &other)
Swap internal storage with another vector, which need not be the same size.
Definition: ConjGradient.h:168
bool operator()(const SizeRange &range, bool finite) const
Definition: ConjGradient.h:737
ValueType infNorm() const
Return the infinity norm of this vector.
Definition: ConjGradient.h:723
const T & operator[](SizeType i) const
Return the value of this vector's ith element.
Definition: ConjGradient.h:204
RowEditor(ValueType *valueHead, SizeType *columnHead, SizeType &rowSize, SizeType colSize)
bool isApproxZero(const Type &x)
Return true if x is equal to zero to within the default floating-point comparison tolerance...
Definition: Math.h:350
void setValue(SizeType row, SizeType col, const ValueType &)
Set the value at the given coordinates.
Definition: ConjGradient.h:878
GLenum GLenum void void * column
Definition: glew.h:4995
Vector & operator=(const Vector &)
Deep copy the given vector.
Definition: ConjGradient.h:568
void fill(const ValueType &value)
Set all elements of this vector to value.
Definition: ConjGradient.h:600
SizeType size() const
Return the number of elements in this vector.
Definition: ConjGradient.h:159
GLuint GLenum matrix
Definition: glew.h:15055
GLdouble n
Definition: glcorearb.h:2007
void scale(const Scalar &s)
Multiply all elements in the matrix by s;.
Definition: ConjGradient.h:925
GLboolean * data
Definition: glcorearb.h:130
GLuint GLfloat * val
Definition: glcorearb.h:1607
#define OPENVDB_LOG_DEBUG_RUNTIME(mesg)
Definition: logging.h:281
void vectorMultiply(const Vector< VecValueType > &inVec, Vector< VecValueType > &resultVec) const
Multiply this matrix by inVec and return the result in resultVec.
Definition: ConjGradient.h:956
SparseStencilMatrix & operator*=(const Scalar &s)
Multiply all elements in the matrix by s;.
Definition: ConjGradient.h:285
ConstRow getConstRow(SizeType row) const
Return a read-only view onto the given row of this matrix.
SizeType setValue(SizeType column, const ValueType &value)
Set the value of the entry in the specified column.
std::string lower(string_view a)
Return an all-upper case version of a (locale-independent).
Definition: strutil.h:342
Preconditioner(const SparseStencilMatrix< T, STENCIL_SIZE > &)
Definition: ConjGradient.h:470
GLenum GLenum void * row
Definition: glew.h:4995
RowEditor & operator*=(const Scalar &s)
Scale all of the entries in this row.
Definition: ConjGradient.h:435
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:659
GA_API const UT_StringHolder N
const GLdouble * m
Definition: glew.h:9166
GLsizei const GLfloat * value
Definition: glcorearb.h:823
bool isFinite() const
Return true if all values along the diagonal are finite.
bool equal(T1 a, T2 b, T3 t)
Definition: ImathFun.h:143
GLenum GLint * range
Definition: glcorearb.h:1924
T operator()(const SizeRange &range, T maxValue) const
Definition: ConjGradient.h:709
bool isFinite() const
Return true if every element of this vector has a finite value.
Definition: ConjGradient.h:753
GLintptr offset
Definition: glcorearb.h:664
std::string str() const
Return a string representation of this matrix.
GLboolean r
Definition: glcorearb.h:1221
GLfloat GLfloat GLfloat alpha
Definition: glcorearb.h:111
GLdouble s
Definition: glew.h:1395
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:114
virtual void apply(const Vector< T > &r, Vector< T > &z)=0
Apply this preconditioner to a residue vector: z = M−1r
State terminationDefaults()
Return default termination conditions for a conjugate gradient solver.
Definition: ConjGradient.h:57
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:598
LinearOp(const T &a_, const T *x_, const T *y_, T *out_)
Definition: ConjGradient.h:522
GLint y
Definition: glcorearb.h:102
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:338
Vector(SizeType n, const ValueType &val)
Construct a vector of n elements and initialize each element to the given value.
Definition: ConjGradient.h:149
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74