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