OpenVDB 11.0.0
Loading...
Searching...
No Matches
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>
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
27namespace openvdb {
29namespace OPENVDB_VERSION_NAME {
30namespace math {
31namespace pcg {
32
34
35using SizeRange = tbb::blocked_range<SizeType>;
36
37template<typename ValueType> class Vector;
38
39template<typename ValueType, SizeType STENCIL_SIZE> class SparseStencilMatrix;
40
41template<typename ValueType> class Preconditioner;
42template<typename MatrixType> class JacobiPreconditioner;
43template<typename MatrixType> class IncompleteCholeskyPreconditioner;
44
45/// Information about the state of a conjugate gradient solution
46struct State {
47 bool success;
51};
52
53
54/// Return default termination conditions for a conjugate gradient solver.
55template<typename ValueType>
56inline 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.
90template<typename PositiveDefMatrix>
91inline State
92solve(
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.
122template<typename PositiveDefMatrix, typename Interrupter>
123inline State
124solve(
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
137template<typename T>
139{
140public:
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.
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.
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.
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,
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
214private:
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).
236template<typename ValueType_, SizeType STENCIL_SIZE>
238{
239public:
240 using ValueType = ValueType_;
243
244 class ConstValueIter;
245 class ConstRow;
246 class RowEditor;
247
248 static inline constexpr ValueType sZeroValue = zeroVal<ValueType>();
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>
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
313private:
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
381public:
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:
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.
429 SizeType setValue(SizeType column, const ValueType& value);
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
442private:
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
463template<typename T>
465{
466public:
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
486namespace internal {
487
488// Functor for use with tbb::parallel_for() to copy data from one array to another
489template<typename T>
490struct 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
504template<typename T>
505struct 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
514 const T val;
515};
516
517
518// Functor for use with tbb::parallel_for() that computes a * x + y
519template<typename T>
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
544inline std::ostream&
545operator<<(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
557template<typename T>
558inline
559Vector<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
566template<typename T>
567inline
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
586template<typename T>
587inline void
589{
590 if (n != mSize) {
591 delete[] mData;
592 mData = new T[n];
593 mSize = n;
594 }
595}
596
597
598template<typename T>
599inline void
601{
602 tbb::parallel_for(SizeRange(0, mSize), internal::FillOp<T>(mData, value));
603}
604
605
606template<typename T>
607template<typename Scalar>
608struct 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
621template<typename T>
622template<typename Scalar>
623inline void
624Vector<T>::scale(const Scalar& s)
625{
626 tbb::parallel_for(SizeRange(0, mSize), ScaleOp<Scalar>(mData, s));
627}
628
629
630template<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
662template<typename T>
663inline T
664Vector<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
704template<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
721template<typename T>
722inline 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
732template<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
751template<typename T>
752inline 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
763template<typename T>
764template<typename OtherValueType>
765struct 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
785template<typename T>
786template<typename OtherValueType>
787inline 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
798template<typename T>
799inline 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
817template<typename ValueType, SizeType STENCIL_SIZE>
818inline
820 : mNumRows(numRows)
821 , mValueArray(new ValueType[mNumRows * STENCIL_SIZE])
822 , mColumnIdxArray(new SizeType[mNumRows * STENCIL_SIZE])
823 , mRowSizeArray(new SizeType[mNumRows])
824{
825 // Initialize the matrix to a null state by setting the size of each row to zero.
826 tbb::parallel_for(SizeRange(0, mNumRows),
827 internal::FillOp<SizeType>(mRowSizeArray.get(), /*value=*/0));
828}
829
830
831template<typename ValueType, SizeType STENCIL_SIZE>
833{
835 from(&from_), to(&to_) {}
836
837 void operator()(const SizeRange& range) const
838 {
839 const ValueType* fromVal = from->mValueArray.get();
840 const SizeType* fromCol = from->mColumnIdxArray.get();
841 ValueType* toVal = to->mValueArray.get();
842 SizeType* toCol = to->mColumnIdxArray.get();
843 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
844 toVal[n] = fromVal[n];
845 toCol[n] = fromCol[n];
846 }
847 }
848
850};
851
852
853template<typename ValueType, SizeType STENCIL_SIZE>
854inline
856 : mNumRows(other.mNumRows)
857 , mValueArray(new ValueType[mNumRows * STENCIL_SIZE])
858 , mColumnIdxArray(new SizeType[mNumRows * STENCIL_SIZE])
859 , mRowSizeArray(new SizeType[mNumRows])
860{
861 SizeType size = mNumRows * STENCIL_SIZE;
862
863 // Copy the value and column index arrays from the other matrix to this matrix.
864 tbb::parallel_for(SizeRange(0, size), MatrixCopyOp(/*from=*/other, /*to=*/*this));
865
866 // Copy the row size array from the other matrix to this matrix.
867 tbb::parallel_for(SizeRange(0, mNumRows),
868 internal::CopyOp<SizeType>(/*from=*/other.mRowSizeArray.get(), /*to=*/mRowSizeArray.get()));
869}
870
871
872template<typename ValueType, SizeType STENCIL_SIZE>
873inline void
875 const ValueType& val)
876{
877 assert(row < mNumRows);
878 this->getRowEditor(row).setValue(col, val);
879}
880
881
882template<typename ValueType, SizeType STENCIL_SIZE>
883inline const ValueType&
885{
886 assert(row < mNumRows);
887 return this->getConstRow(row).getValue(col);
888}
889
890
891template<typename ValueType, SizeType STENCIL_SIZE>
892inline const ValueType&
894{
895 return this->getValue(row,col);
896}
897
898
899template<typename ValueType, SizeType STENCIL_SIZE>
900template<typename Scalar>
901struct SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowScaleOp
902{
903 RowScaleOp(SparseStencilMatrix& m, const Scalar& s_): mat(&m), s(s_) {}
904
905 void operator()(const SizeRange& range) const
906 {
907 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
908 RowEditor row = mat->getRowEditor(n);
909 row.scale(s);
910 }
911 }
912
914 const Scalar s;
915};
916
917
918template<typename ValueType, SizeType STENCIL_SIZE>
919template<typename Scalar>
920inline void
922{
923 // Parallelize over the rows in the matrix.
924 tbb::parallel_for(SizeRange(0, mNumRows), RowScaleOp<Scalar>(*this, s));
925}
926
927
928template<typename ValueType, SizeType STENCIL_SIZE>
929template<typename VecValueType>
930struct SparseStencilMatrix<ValueType, STENCIL_SIZE>::VecMultOp
931{
932 VecMultOp(const SparseStencilMatrix& m, const VecValueType* i, VecValueType* o):
933 mat(&m), in(i), out(o) {}
934
935 void operator()(const SizeRange& range) const
936 {
937 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
938 ConstRow row = mat->getConstRow(n);
939 out[n] = row.dot(in, mat->numRows());
940 }
941 }
942
943 const SparseStencilMatrix* mat;
944 const VecValueType* in;
945 VecValueType* out;
946};
947
948
949template<typename ValueType, SizeType STENCIL_SIZE>
950template<typename VecValueType>
951inline void
953 const Vector<VecValueType>& inVec, Vector<VecValueType>& resultVec) const
954{
955 if (inVec.size() != mNumRows) {
956 OPENVDB_THROW(ArithmeticError, "matrix and input vector have incompatible sizes ("
957 << mNumRows << "x" << mNumRows << " vs. " << inVec.size() << ")");
958 }
959 if (resultVec.size() != mNumRows) {
960 OPENVDB_THROW(ArithmeticError, "matrix and result vector have incompatible sizes ("
961 << mNumRows << "x" << mNumRows << " vs. " << resultVec.size() << ")");
962 }
963
964 vectorMultiply(inVec.data(), resultVec.data());
965}
966
967
968template<typename ValueType, SizeType STENCIL_SIZE>
969template<typename VecValueType>
970inline void
972 const VecValueType* inVec, VecValueType* resultVec) const
973{
974 // Parallelize over the rows in the matrix.
975 tbb::parallel_for(SizeRange(0, mNumRows),
976 VecMultOp<VecValueType>(*this, inVec, resultVec));
977}
978
979
980template<typename ValueType, SizeType STENCIL_SIZE>
981template<typename OtherValueType>
982struct SparseStencilMatrix<ValueType, STENCIL_SIZE>::EqOp
983{
984 EqOp(const SparseStencilMatrix& a_,
986 a(&a_), b(&b_), eps(e) {}
987
988 bool operator()(const SizeRange& range, bool equal) const
989 {
990 if (equal) {
991 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
992 if (!a->getConstRow(n).eq(b->getConstRow(n), eps)) return false;
993 }
994 }
995 return equal;
996 }
997
998 const SparseStencilMatrix* a;
999 const SparseStencilMatrix<OtherValueType, STENCIL_SIZE>* b;
1000 const ValueType eps;
1001};
1002
1003
1004template<typename ValueType, SizeType STENCIL_SIZE>
1005template<typename OtherValueType>
1006inline bool
1009{
1010 if (this->numRows() != other.numRows()) return false;
1011 bool equal = tbb::parallel_reduce(SizeRange(0, this->numRows()), /*seed=*/true,
1012 EqOp<OtherValueType>(*this, other, eps),
1013 /*join=*/[](bool eq1, bool eq2) { return (eq1 && eq2); });
1014 return equal;
1015}
1016
1017
1018template<typename ValueType, SizeType STENCIL_SIZE>
1020{
1021 IsFiniteOp(const SparseStencilMatrix& m): mat(&m) {}
1022
1023 bool operator()(const SizeRange& range, bool finite) const
1024 {
1025 if (finite) {
1026 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1027 const ConstRow row = mat->getConstRow(n);
1028 for (ConstValueIter it = row.cbegin(); it; ++it) {
1029 if (!std::isfinite(*it)) return false;
1030 }
1031 }
1032 }
1033 return finite;
1034 }
1035
1037};
1038
1039
1040template<typename ValueType, SizeType STENCIL_SIZE>
1041inline bool
1043{
1044 // Parallelize over the rows of this matrix.
1045 bool finite = tbb::parallel_reduce(SizeRange(0, this->numRows()), /*seed=*/true,
1046 IsFiniteOp(*this), /*join=*/[](bool finite1, bool finite2) { return (finite1&&finite2); });
1047 return finite;
1048}
1049
1050
1051template<typename ValueType, SizeType STENCIL_SIZE>
1052inline std::string
1054{
1055 std::ostringstream ostr;
1056 for (SizeType n = 0, N = this->size(); n < N; ++n) {
1057 ostr << n << ": " << this->getConstRow(n).str() << "\n";
1058 }
1059 return ostr.str();
1060}
1061
1062
1063template<typename ValueType, SizeType STENCIL_SIZE>
1066{
1067 assert(i < mNumRows);
1068 const SizeType head = i * STENCIL_SIZE;
1069 return RowEditor(&mValueArray[head], &mColumnIdxArray[head], mRowSizeArray[i], mNumRows);
1070}
1071
1072
1073template<typename ValueType, SizeType STENCIL_SIZE>
1076{
1077 assert(i < mNumRows);
1078 const SizeType head = i * STENCIL_SIZE; // index for this row into main storage
1079 return ConstRow(&mValueArray[head], &mColumnIdxArray[head], mRowSizeArray[i]);
1080}
1081
1082
1083template<typename ValueType, SizeType STENCIL_SIZE>
1084template<typename DataType>
1085inline SizeType
1087{
1088 if (this->empty()) return mData.mSize;
1089
1090 // Get a pointer to the first column index that is equal to or greater than the given index.
1091 // (This assumes that the data is sorted by column.)
1092 const SizeType* colPtr = std::lower_bound(mData.mCols, mData.mCols + mData.mSize, columnIdx);
1093 // Return the offset of the pointer from the beginning of the array.
1094 return static_cast<SizeType>(colPtr - mData.mCols);
1095}
1096
1097
1098template<typename ValueType, SizeType STENCIL_SIZE>
1099template<typename DataType>
1100inline const ValueType&
1102 SizeType columnIdx, bool& active) const
1103{
1104 active = false;
1105 SizeType idx = this->find(columnIdx);
1106 if (idx < this->size() && this->column(idx) == columnIdx) {
1107 active = true;
1108 return this->value(idx);
1109 }
1111}
1112
1113template<typename ValueType, SizeType STENCIL_SIZE>
1114template<typename DataType>
1115inline const ValueType&
1117{
1118 SizeType idx = this->find(columnIdx);
1119 if (idx < this->size() && this->column(idx) == columnIdx) {
1120 return this->value(idx);
1121 }
1123}
1124
1125
1126template<typename ValueType, SizeType STENCIL_SIZE>
1127template<typename DataType>
1128inline typename SparseStencilMatrix<ValueType, STENCIL_SIZE>::ConstValueIter
1129SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::cbegin() const
1130{
1131 return ConstValueIter(mData);
1132}
1133
1134
1135template<typename ValueType, SizeType STENCIL_SIZE>
1136template<typename DataType>
1137template<typename OtherDataType>
1138inline bool
1140 const RowBase<OtherDataType>& other, ValueType eps) const
1141{
1142 if (this->size() != other.size()) return false;
1143 for (ConstValueIter it = cbegin(), oit = other.cbegin(); it || oit; ++it, ++oit) {
1144 if (it.column() != oit.column()) return false;
1145 if (!isApproxEqual(*it, *oit, eps)) return false;
1146 }
1147 return true;
1148}
1149
1150
1151template<typename ValueType, SizeType STENCIL_SIZE>
1152template<typename DataType>
1153template<typename VecValueType>
1154inline VecValueType
1155SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::dot(
1156 const VecValueType* inVec, SizeType vecSize) const
1157{
1158 VecValueType result = zeroVal<VecValueType>();
1159 for (SizeType idx = 0, N = std::min(vecSize, this->size()); idx < N; ++idx) {
1160 result += static_cast<VecValueType>(this->value(idx) * inVec[this->column(idx)]);
1161 }
1162 return result;
1163}
1164
1165template<typename ValueType, SizeType STENCIL_SIZE>
1166template<typename DataType>
1167template<typename VecValueType>
1168inline VecValueType
1169SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::dot(
1170 const Vector<VecValueType>& inVec) const
1171{
1172 return dot(inVec.data(), inVec.size());
1173}
1174
1175
1176template<typename ValueType, SizeType STENCIL_SIZE>
1177template<typename DataType>
1178inline std::string
1180{
1181 std::ostringstream ostr;
1182 std::string sep;
1183 for (SizeType n = 0, N = this->size(); n < N; ++n) {
1184 ostr << sep << "(" << this->column(n) << ", " << this->value(n) << ")";
1185 sep = ", ";
1186 }
1187 return ostr.str();
1188}
1189
1190
1191template<typename ValueType, SizeType STENCIL_SIZE>
1192inline
1194 const ValueType* valueHead, const SizeType* columnHead, const SizeType& rowSize):
1195 ConstRowBase(ConstRowData(const_cast<ValueType*>(valueHead),
1196 const_cast<SizeType*>(columnHead), const_cast<SizeType&>(rowSize)))
1197{
1198}
1199
1200
1201template<typename ValueType, SizeType STENCIL_SIZE>
1202inline
1204 ValueType* valueHead, SizeType* columnHead, SizeType& rowSize, SizeType colSize):
1205 RowBase<>(RowData(valueHead, columnHead, rowSize)), mNumColumns(colSize)
1206{
1207}
1208
1209
1210template<typename ValueType, SizeType STENCIL_SIZE>
1211inline void
1213{
1214 // Note: since mSize is a reference, this modifies the underlying matrix.
1215 RowBase<>::mData.mSize = 0;
1216}
1217
1218
1219template<typename ValueType, SizeType STENCIL_SIZE>
1220inline SizeType
1222 SizeType column, const ValueType& value)
1223{
1224 assert(column < mNumColumns);
1225
1226 RowData& data = RowBase<>::mData;
1227
1228 // Get the offset of the first column index that is equal to or greater than
1229 // the column to be modified.
1230 SizeType offset = this->find(column);
1231
1232 if (offset < data.mSize && data.mCols[offset] == column) {
1233 // If the column already exists, just update its value.
1234 data.mVals[offset] = value;
1235 return data.mSize;
1236 }
1237
1238 // Check that it is safe to add a new column.
1239 assert(data.mSize < this->capacity());
1240
1241 if (offset >= data.mSize) {
1242 // The new column's index is larger than any existing index. Append the new column.
1243 data.mVals[data.mSize] = value;
1244 data.mCols[data.mSize] = column;
1245 } else {
1246 // Insert the new column at the computed offset after shifting subsequent columns.
1247 for (SizeType i = data.mSize; i > offset; --i) {
1248 data.mVals[i] = data.mVals[i - 1];
1249 data.mCols[i] = data.mCols[i - 1];
1250 }
1251 data.mVals[offset] = value;
1252 data.mCols[offset] = column;
1253 }
1254 ++data.mSize;
1255
1256 return data.mSize;
1257}
1258
1259
1260template<typename ValueType, SizeType STENCIL_SIZE>
1261template<typename Scalar>
1262inline void
1264{
1265 for (int idx = 0, N = this->size(); idx < N; ++idx) {
1266 RowBase<>::mData.mVals[idx] *= s;
1267 }
1268}
1269
1270
1271////////////////////////////////////////
1272
1273
1274/// Diagonal preconditioner
1275template<typename MatrixType>
1276class JacobiPreconditioner: public Preconditioner<typename MatrixType::ValueType>
1277{
1278private:
1279 struct InitOp;
1280 struct ApplyOp;
1281
1282public:
1283 using ValueType = typename MatrixType::ValueType;
1287
1288 JacobiPreconditioner(const MatrixType& A): BaseType(A), mDiag(A.numRows())
1289 {
1290 // Initialize vector mDiag with the values from the matrix diagonal.
1291 tbb::parallel_for(SizeRange(0, A.numRows()), InitOp(A, mDiag.data()));
1292 }
1293
1294 ~JacobiPreconditioner() override = default;
1295
1296 void apply(const Vector<ValueType>& r, Vector<ValueType>& z) override
1297 {
1298 const SizeType size = mDiag.size();
1299
1300 assert(r.size() == z.size());
1301 assert(r.size() == size);
1302
1303 tbb::parallel_for(SizeRange(0, size), ApplyOp(mDiag.data(), r.data(), z.data()));
1304 }
1305
1306 /// Return @c true if all values along the diagonal are finite.
1307 bool isFinite() const { return mDiag.isFinite(); }
1308
1309private:
1310 // Functor for use with tbb::parallel_for()
1311 struct InitOp
1312 {
1313 InitOp(const MatrixType& m, ValueType* v): mat(&m), vec(v) {}
1314 void operator()(const SizeRange& range) const {
1315 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1316 const ValueType val = mat->getValue(n, n);
1317 assert(!isApproxZero(val, ValueType(0.0001)));
1318 vec[n] = static_cast<ValueType>(1.0 / val);
1319 }
1320 }
1321 const MatrixType* mat; ValueType* vec;
1322 };
1323
1324 // Functor for use with tbb::parallel_reduce()
1325 struct ApplyOp
1326 {
1327 ApplyOp(const ValueType* x_, const ValueType* y_, ValueType* out_):
1328 x(x_), y(y_), out(out_) {}
1329 void operator()(const SizeRange& range) const {
1330 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) out[n] = x[n] * y[n];
1331 }
1332 const ValueType *x, *y; ValueType* out;
1333 };
1334
1335 // The Jacobi preconditioner is a diagonal matrix
1336 VectorType mDiag;
1337}; // class JacobiPreconditioner
1338
1339
1340/// Preconditioner using incomplete Cholesky factorization
1341template<typename MatrixType>
1342class IncompleteCholeskyPreconditioner: public Preconditioner<typename MatrixType::ValueType>
1343{
1344private:
1345 struct CopyToLowerOp;
1346 struct TransposeOp;
1347
1348public:
1349 using ValueType = typename MatrixType::ValueType;
1354 using TriangleConstRow = typename TriangularMatrix::ConstRow;
1355 using TriangleRowEditor = typename TriangularMatrix::RowEditor;
1356
1357 IncompleteCholeskyPreconditioner(const MatrixType& matrix)
1358 : BaseType(matrix)
1359 , mLowerTriangular(matrix.numRows())
1360 , mUpperTriangular(matrix.numRows())
1361 , mTempVec(matrix.numRows())
1362 {
1363 // Size of matrix
1364 const SizeType numRows = mLowerTriangular.numRows();
1365
1366 // Copy the upper triangular part to the lower triangular part.
1367 tbb::parallel_for(SizeRange(0, numRows), CopyToLowerOp(matrix, mLowerTriangular));
1368
1369 // Build the Incomplete Cholesky Matrix
1370 //
1371 // Algorithm:
1372 //
1373 // for (k = 0; k < size; ++k) {
1374 // A(k,k) = sqrt(A(k,k));
1375 // for (i = k +1, i < size; ++i) {
1376 // if (A(i,k) == 0) continue;
1377 // A(i,k) = A(i,k) / A(k,k);
1378 // }
1379 // for (j = k+1; j < size; ++j) {
1380 // for (i = j; i < size; ++i) {
1381 // if (A(i,j) == 0) continue;
1382 // A(i,j) -= A(i,k)*A(j,k);
1383 // }
1384 // }
1385 // }
1386
1387 mPassedCompatibilityCondition = true;
1388
1389 for (SizeType k = 0; k < numRows; ++k) {
1390
1391 TriangleConstRow crow_k = mLowerTriangular.getConstRow(k);
1392 ValueType diagonalValue = crow_k.getValue(k);
1393
1394 // Test if the matrix build has failed.
1395 if (diagonalValue < 1.e-5) {
1396 mPassedCompatibilityCondition = false;
1397 break;
1398 }
1399
1400 diagonalValue = Sqrt(diagonalValue);
1401
1402 TriangleRowEditor row_k = mLowerTriangular.getRowEditor(k);
1403 row_k.setValue(k, diagonalValue);
1404
1405 // Exploit the fact that the matrix is symmetric.
1406 typename MatrixType::ConstRow srcRow = matrix.getConstRow(k);
1407 typename MatrixType::ConstValueIter citer = srcRow.cbegin();
1408 for ( ; citer; ++citer) {
1409 SizeType ii = citer.column();
1410 if (ii < k+1) continue; // look above diagonal
1411
1412 TriangleRowEditor row_ii = mLowerTriangular.getRowEditor(ii);
1413
1414 row_ii.setValue(k, *citer / diagonalValue);
1415 }
1416
1417 // for (j = k+1; j < size; ++j) replaced by row iter below
1418 citer.reset(); // k,j entries
1419 for ( ; citer; ++citer) {
1420 SizeType j = citer.column();
1421 if (j < k+1) continue;
1422
1423 TriangleConstRow row_j = mLowerTriangular.getConstRow(j);
1424 ValueType a_jk = row_j.getValue(k); // a_jk is non zero if a_kj is non zero
1425
1426 // Entry (i,j) is non-zero if matrix(j,i) is nonzero
1427
1428 typename MatrixType::ConstRow mask = matrix.getConstRow(j);
1429 typename MatrixType::ConstValueIter maskIter = mask.cbegin();
1430 for ( ; maskIter; ++maskIter) {
1431 SizeType i = maskIter.column();
1432 if (i < j) continue;
1433
1434 TriangleConstRow crow_i = mLowerTriangular.getConstRow(i);
1435 ValueType a_ij = crow_i.getValue(j);
1436 ValueType a_ik = crow_i.getValue(k);
1437 TriangleRowEditor row_i = mLowerTriangular.getRowEditor(i);
1438 a_ij -= a_ik * a_jk;
1439
1440 row_i.setValue(j, a_ij);
1441 }
1442 }
1443 }
1444
1445 // Build the transpose of the IC matrix: mUpperTriangular
1446 tbb::parallel_for(SizeRange(0, numRows),
1447 TransposeOp(matrix, mLowerTriangular, mUpperTriangular));
1448 }
1449
1451
1452 bool isValid() const override { return mPassedCompatibilityCondition; }
1453
1454 void apply(const Vector<ValueType>& rVec, Vector<ValueType>& zVec) override
1455 {
1456 if (!mPassedCompatibilityCondition) {
1457 OPENVDB_THROW(ArithmeticError, "invalid Cholesky decomposition");
1458 }
1459
1460 // Solve mUpperTriangular * mLowerTriangular * rVec = zVec;
1461
1462 SizeType size = mLowerTriangular.numRows();
1463
1464 zVec.fill(zeroVal<ValueType>());
1465 ValueType* zData = zVec.data();
1466
1467 if (size == 0) return;
1468
1469 assert(rVec.size() == size);
1470 assert(zVec.size() == size);
1471
1472 // Allocate a temp vector
1473 mTempVec.fill(zeroVal<ValueType>());
1474 ValueType* tmpData = mTempVec.data();
1475 const ValueType* rData = rVec.data();
1476
1477 // Solve mLowerTriangular * tmp = rVec;
1478 for (SizeType i = 0; i < size; ++i) {
1479 typename TriangularMatrix::ConstRow row = mLowerTriangular.getConstRow(i);
1480 ValueType diagonal = row.getValue(i);
1481 ValueType dot = row.dot(mTempVec);
1482 tmpData[i] = (rData[i] - dot) / diagonal;
1483 if (!std::isfinite(tmpData[i])) {
1484 OPENVDB_LOG_DEBUG_RUNTIME("1 diagonal was " << diagonal);
1485 OPENVDB_LOG_DEBUG_RUNTIME("1a diagonal " << row.getValue(i));
1486 }
1487 }
1488
1489 // Solve mUpperTriangular * zVec = tmp;
1490 for (SizeType ii = 0; ii < size; ++ii) {
1491 SizeType i = size - 1 - ii;
1492 typename TriangularMatrix::ConstRow row = mUpperTriangular.getConstRow(i);
1493 ValueType diagonal = row.getValue(i);
1494 ValueType dot = row.dot(zVec);
1495 zData[i] = (tmpData[i] - dot) / diagonal;
1496 if (!std::isfinite(zData[i])) {
1497 OPENVDB_LOG_DEBUG_RUNTIME("2 diagonal was " << diagonal);
1498 }
1499 }
1500 }
1501
1502 const TriangularMatrix& lowerMatrix() const { return mLowerTriangular; }
1503 const TriangularMatrix& upperMatrix() const { return mUpperTriangular; }
1504
1505private:
1506 // Functor for use with tbb::parallel_for()
1507 struct CopyToLowerOp
1508 {
1509 CopyToLowerOp(const MatrixType& m, TriangularMatrix& l): mat(&m), lower(&l) {}
1510 void operator()(const SizeRange& range) const {
1511 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1512 typename TriangularMatrix::RowEditor outRow = lower->getRowEditor(n);
1513 outRow.clear();
1514 typename MatrixType::ConstRow inRow = mat->getConstRow(n);
1515 for (typename MatrixType::ConstValueIter it = inRow.cbegin(); it; ++it) {
1516 if (it.column() > n) continue; // skip above diagonal
1517 outRow.setValue(it.column(), *it);
1518 }
1519 }
1520 }
1521 const MatrixType* mat; TriangularMatrix* lower;
1522 };
1523
1524 // Functor for use with tbb::parallel_for()
1525 struct TransposeOp
1526 {
1527 TransposeOp(const MatrixType& m, const TriangularMatrix& l, TriangularMatrix& u):
1528 mat(&m), lower(&l), upper(&u) {}
1529 void operator()(const SizeRange& range) const {
1530 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1531 typename TriangularMatrix::RowEditor outRow = upper->getRowEditor(n);
1532 outRow.clear();
1533 // Use the fact that matrix is symmetric.
1534 typename MatrixType::ConstRow inRow = mat->getConstRow(n);
1535 for (typename MatrixType::ConstValueIter it = inRow.cbegin(); it; ++it) {
1536 const SizeType column = it.column();
1537 if (column < n) continue; // only set upper triangle
1538 outRow.setValue(column, lower->getValue(column, n));
1539 }
1540 }
1541 }
1542 const MatrixType* mat; const TriangularMatrix* lower; TriangularMatrix* upper;
1543 };
1544
1545 TriangularMatrix mLowerTriangular;
1546 TriangularMatrix mUpperTriangular;
1547 Vector<ValueType> mTempVec;
1548 bool mPassedCompatibilityCondition;
1549}; // class IncompleteCholeskyPreconditioner
1550
1551
1552////////////////////////////////////////
1553
1554
1555namespace internal {
1556
1557/// Compute @e ax + @e y.
1558template<typename T>
1559inline void
1560axpy(const T& a, const T* xVec, const T* yVec, T* resultVec, SizeType size)
1561{
1562 tbb::parallel_for(SizeRange(0, size), LinearOp<T>(a, xVec, yVec, resultVec));
1563}
1564
1565/// Compute @e ax + @e y.
1566template<typename T>
1567inline void
1568axpy(const T& a, const Vector<T>& xVec, const Vector<T>& yVec, Vector<T>& result)
1569{
1570 assert(xVec.size() == yVec.size());
1571 assert(xVec.size() == result.size());
1572 axpy(a, xVec.data(), yVec.data(), result.data(), xVec.size());
1573}
1574
1575
1576/// Compute @e r = @e b &minus; @e Ax.
1577template<typename MatrixOperator, typename VecValueType>
1578inline void
1579computeResidual(const MatrixOperator& A, const VecValueType* x,
1580 const VecValueType* b, VecValueType* r)
1581{
1582 // Compute r = A * x.
1583 A.vectorMultiply(x, r);
1584 // Compute r = b - r.
1585 tbb::parallel_for(SizeRange(0, A.numRows()), LinearOp<VecValueType>(-1.0, r, b, r));
1586}
1587
1588/// Compute @e r = @e b &minus; @e Ax.
1589template<typename MatrixOperator, typename T>
1590inline void
1591computeResidual(const MatrixOperator& A, const Vector<T>& x, const Vector<T>& b, Vector<T>& r)
1592{
1593 assert(x.size() == b.size());
1594 assert(x.size() == r.size());
1595 assert(x.size() == A.numRows());
1596
1597 computeResidual(A, x.data(), b.data(), r.data());
1598}
1599
1600} // namespace internal
1601
1602
1603////////////////////////////////////////
1604
1605
1606template<typename PositiveDefMatrix>
1607inline State
1609 const PositiveDefMatrix& Amat,
1613 const State& termination)
1614{
1615 util::NullInterrupter interrupter;
1616 return solve(Amat, bVec, xVec, precond, interrupter, termination);
1617}
1618
1619
1620template<typename PositiveDefMatrix, typename Interrupter>
1621inline State
1623 const PositiveDefMatrix& Amat,
1627 Interrupter& interrupter,
1628 const State& termination)
1629{
1630 using ValueType = typename PositiveDefMatrix::ValueType;
1632
1633 State result;
1634 result.success = false;
1635 result.iterations = 0;
1636 result.relativeError = 0.0;
1637 result.absoluteError = 0.0;
1638
1639 const SizeType size = Amat.numRows();
1640 if (size == 0) {
1641 OPENVDB_LOG_WARN("pcg::solve(): matrix has dimension zero");
1642 return result;
1643 }
1644 if (size != bVec.size()) {
1645 OPENVDB_THROW(ArithmeticError, "A and b have incompatible sizes"
1646 << size << "x" << size << " vs. " << bVec.size() << ")");
1647 }
1648 if (size != xVec.size()) {
1649 OPENVDB_THROW(ArithmeticError, "A and x have incompatible sizes"
1650 << size << "x" << size << " vs. " << xVec.size() << ")");
1651 }
1652
1653 // Temp vectors
1654 VectorType zVec(size); // transformed residual (M^-1 r)
1655 VectorType pVec(size); // search direction
1656 VectorType qVec(size); // A * p
1657
1658 // Compute norm of B (the source)
1659 const ValueType tmp = bVec.infNorm();
1660 const ValueType infNormOfB = isZero(tmp) ? ValueType(1) : tmp;
1661
1662 // Compute rVec: residual = b - Ax.
1663 VectorType rVec(size); // vector of residuals
1664
1665 internal::computeResidual(Amat, xVec, bVec, rVec);
1666
1667 assert(rVec.isFinite());
1668
1669 // Normalize the residual norm with the source norm and look for early out.
1670 result.absoluteError = static_cast<double>(rVec.infNorm());
1671 result.relativeError = static_cast<double>(result.absoluteError / infNormOfB);
1672 if (result.relativeError <= termination.relativeError) {
1673 result.success = true;
1674 return result;
1675 }
1676
1677 // Iterations of the CG solve
1678
1679 ValueType rDotZPrev(1); // inner product of <z,r>
1680
1681 // Keep track of the minimum error to monitor convergence.
1682 ValueType minL2Error = std::numeric_limits<ValueType>::max();
1683 ValueType l2Error;
1684
1685 int iteration = 0;
1686 for ( ; iteration < termination.iterations; ++iteration) {
1687
1688 if (interrupter.wasInterrupted()) {
1689 OPENVDB_THROW(RuntimeError, "conjugate gradient solver was interrupted");
1690 }
1691
1692 OPENVDB_LOG_DEBUG_RUNTIME("pcg::solve() " << result);
1693
1694 result.iterations = iteration + 1;
1695
1696 // Apply preconditioner to residual
1697 // z_{k} = M^-1 r_{k}
1698 precond.apply(rVec, zVec);
1699
1700 // <r,z>
1701 const ValueType rDotZ = rVec.dot(zVec);
1702 assert(std::isfinite(rDotZ));
1703
1704 if (0 == iteration) {
1705 // Initialize
1706 pVec = zVec;
1707 } else {
1708 const ValueType beta = rDotZ / rDotZPrev;
1709 // p = beta * p + z
1710 internal::axpy(beta, pVec, zVec, /*result */pVec);
1711 }
1712
1713 // q_{k} = A p_{k}
1714 Amat.vectorMultiply(pVec, qVec);
1715
1716 // alpha = <r_{k-1}, z_{k-1}> / <p_{k},q_{k}>
1717 const ValueType pAp = pVec.dot(qVec);
1718 assert(std::isfinite(pAp));
1719
1720 const ValueType alpha = rDotZ / pAp;
1721 rDotZPrev = rDotZ;
1722
1723 // x_{k} = x_{k-1} + alpha * p_{k}
1724 internal::axpy(alpha, pVec, xVec, /*result=*/xVec);
1725
1726 // r_{k} = r_{k-1} - alpha_{k-1} A p_{k}
1727 internal::axpy(-alpha, qVec, rVec, /*result=*/rVec);
1728
1729 // update tolerances
1730 l2Error = rVec.l2Norm();
1731 minL2Error = Min(l2Error, minL2Error);
1732
1733 result.absoluteError = static_cast<double>(rVec.infNorm());
1734 result.relativeError = static_cast<double>(result.absoluteError / infNormOfB);
1735
1736 if (l2Error > 2 * minL2Error) {
1737 // The solution started to diverge.
1738 result.success = false;
1739 break;
1740 }
1741 if (!std::isfinite(result.absoluteError)) {
1742 // Total divergence of solution
1743 result.success = false;
1744 break;
1745 }
1746 if (result.absoluteError <= termination.absoluteError) {
1747 // Convergence
1748 result.success = true;
1749 break;
1750 }
1751 if (result.relativeError <= termination.relativeError) {
1752 // Convergence
1753 result.success = true;
1754 break;
1755 }
1756 }
1757 OPENVDB_LOG_DEBUG_RUNTIME("pcg::solve() " << result);
1758
1759 return result;
1760}
1761
1762} // namespace pcg
1763} // namespace math
1764} // namespace OPENVDB_VERSION_NAME
1765} // namespace openvdb
1766
1767#endif // OPENVDB_MATH_CONJGRADIENT_HAS_BEEN_INCLUDED
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
Definition Exceptions.h:56
Definition Exceptions.h:63
Preconditioner using incomplete Cholesky factorization.
Definition ConjGradient.h:1343
typename TriangularMatrix::RowEditor TriangleRowEditor
Definition ConjGradient.h:1355
SharedPtr< IncompleteCholeskyPreconditioner > Ptr
Definition ConjGradient.h:1352
IncompleteCholeskyPreconditioner(const MatrixType &matrix)
Definition ConjGradient.h:1357
typename TriangularMatrix::ConstRow TriangleConstRow
Definition ConjGradient.h:1354
void apply(const Vector< ValueType > &rVec, Vector< ValueType > &zVec) override
Apply this preconditioner to a residue vector: z = M−1r
Definition ConjGradient.h:1454
typename MatrixType::ValueType ValueType
Definition ConjGradient.h:1349
const TriangularMatrix & upperMatrix() const
Definition ConjGradient.h:1503
bool isValid() const override
Definition ConjGradient.h:1452
const TriangularMatrix & lowerMatrix() const
Definition ConjGradient.h:1502
Diagonal preconditioner.
Definition ConjGradient.h:1277
SharedPtr< JacobiPreconditioner > Ptr
Definition ConjGradient.h:1286
void apply(const Vector< ValueType > &r, Vector< ValueType > &z) override
Apply this preconditioner to a residue vector: z = M−1r
Definition ConjGradient.h:1296
typename MatrixType::ValueType ValueType
Definition ConjGradient.h:1283
JacobiPreconditioner(const MatrixType &A)
Definition ConjGradient.h:1288
bool isFinite() const
Return true if all values along the diagonal are finite.
Definition ConjGradient.h:1307
Base class for conjugate gradient preconditioners.
Definition ConjGradient.h:465
Preconditioner(const SparseStencilMatrix< T, STENCIL_SIZE > &)
Definition ConjGradient.h:470
virtual bool isValid() const
Definition ConjGradient.h:473
virtual void apply(const Vector< T > &r, Vector< T > &z)=0
Apply this preconditioner to a residue vector: z = M−1r
SharedPtr< Preconditioner > Ptr
Definition ConjGradient.h:468
T ValueType
Definition ConjGradient.h:467
Read-only accessor to a row of this matrix.
Definition ConjGradient.h:412
ConstRow(const ValueType *valueHead, const SizeType *columnHead, const SizeType &rowSize)
Definition ConjGradient.h:1193
Iterator over the stored values in a row of this matrix.
Definition ConjGradient.h:384
const ValueType & operator*() const
Definition ConjGradient.h:386
ConstValueIter & operator++()
Definition ConjGradient.h:395
SizeType column() const
Definition ConjGradient.h:392
Read/write accessor to a row of this matrix.
Definition ConjGradient.h:420
SizeType setValue(SizeType column, const ValueType &value)
Set the value of the entry in the specified column.
Definition ConjGradient.h:1221
RowEditor & operator*=(const Scalar &s)
Definition ConjGradient.h:435
void scale(const Scalar &)
Scale all of the entries in this row.
Definition ConjGradient.h:1263
void clear()
Set the number of entries in this row to zero.
Definition ConjGradient.h:1212
RowEditor(ValueType *valueHead, SizeType *columnHead, SizeType &rowSize, SizeType colSize)
Definition ConjGradient.h:1203
Sparse, square matrix representing a 3D stencil operator of size STENCIL_SIZE.
Definition ConjGradient.h:238
SizeType size() const
Definition ConjGradient.h:259
const ValueType & getValue(SizeType row, SizeType col) const
Return the value at the given coordinates.
Definition ConjGradient.h:884
SparseStencilMatrix(const SparseStencilMatrix &)
Deep copy the given matrix.
Definition ConjGradient.h:855
void setValue(SizeType row, SizeType col, const ValueType &)
Set the value at the given coordinates.
Definition ConjGradient.h:874
ConstRow getConstRow(SizeType row) const
Return a read-only view onto the given row of this matrix.
Definition ConjGradient.h:1075
SparseStencilMatrix(SizeType n)
Construct an n x n matrix with at most STENCIL_SIZE nonzero elements per row.
Definition ConjGradient.h:819
const ValueType & operator()(SizeType row, SizeType col) const
Definition ConjGradient.h:893
RowEditor getRowEditor(SizeType row)
Return a read/write view onto the given row of this matrix.
Definition ConjGradient.h:1065
ValueType_ ValueType
Definition ConjGradient.h:240
void vectorMultiply(const Vector< VecValueType > &inVec, Vector< VecValueType > &resultVec) const
Multiply this matrix by inVec and return the result in resultVec.
Definition ConjGradient.h:952
SparseStencilMatrix & operator*=(const Scalar &s)
Definition ConjGradient.h:285
static constexpr ValueType sZeroValue
Definition ConjGradient.h:248
void scale(const Scalar &s)
Multiply all elements in the matrix by s;.
Definition ConjGradient.h:921
Vector< ValueType > VectorType
Definition ConjGradient.h:241
void vectorMultiply(const VecValueType *inVec, VecValueType *resultVec) const
Multiply this matrix by the vector represented by the array inVec and return the result in resultVec.
Definition ConjGradient.h:971
SharedPtr< SparseStencilMatrix > Ptr
Definition ConjGradient.h:242
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.
Definition ConjGradient.h:1007
std::string str() const
Return a string representation of this matrix.
Definition ConjGradient.h:1053
SizeType numRows() const
Return the number of rows in this matrix.
Definition ConjGradient.h:258
bool isFinite() const
Return true if every element of this matrix has a finite value.
Definition ConjGradient.h:1042
Lightweight, variable-length vector.
Definition ConjGradient.h:139
SizeType size() const
Return the number of elements in this vector.
Definition ConjGradient.h:159
void swap(Vector &other)
Swap internal storage with another vector, which need not be the same size.
Definition ConjGradient.h:168
ValueType l2Norm() const
Return the L2 norm of this vector.
Definition ConjGradient.h:185
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
SharedPtr< Vector > Ptr
Definition ConjGradient.h:142
const T * data() const
Definition ConjGradient.h:210
Vector()
Construct an empty vector.
Definition ConjGradient.h:145
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 T & operator[](SizeType i) const
Definition ConjGradient.h:204
void resize(SizeType n)
Reset this vector to have n elements, with uninitialized values.
Definition ConjGradient.h:588
const T & at(SizeType i) const
Definition ConjGradient.h:202
bool empty() const
Return true if this vector has no elements.
Definition ConjGradient.h:161
T & at(SizeType i)
Return the value of this vector's ith element.
Definition ConjGradient.h:201
Vector(SizeType n, const ValueType &val)
Construct a vector of n elements and initialize each element to the given value.
Definition ConjGradient.h:149
Vector(const Vector &)
Deep copy the given vector.
Definition ConjGradient.h:559
void scale(const Scalar &s)
Multiply each element of this vector by s.
Definition ConjGradient.h:624
const T * constData() const
Definition ConjGradient.h:211
~Vector()
Definition ConjGradient.h:151
Vector & operator*=(const Scalar &s)
Definition ConjGradient.h:176
ValueType infNorm() const
Return the infinity norm of this vector.
Definition ConjGradient.h:723
void fill(const ValueType &value)
Set all elements of this vector to value.
Definition ConjGradient.h:600
Vector & operator=(const Vector &)
Deep copy the given vector.
Definition ConjGradient.h:568
Vector(SizeType n)
Construct a vector of n elements, with uninitialized values.
Definition ConjGradient.h:147
T & operator[](SizeType i)
Definition ConjGradient.h:203
T * data()
Return a pointer to this vector's elements.
Definition ConjGradient.h:209
std::string str() const
Return a string representation of this vector.
Definition ConjGradient.h:800
T ValueType
Definition ConjGradient.h:141
bool isFinite() const
Return true if every element of this vector has a finite value.
Definition ConjGradient.h:753
#define OPENVDB_LOG_WARN(message)
Log a warning message of the form 'someVar << "some text" << ...'.
Definition logging.h:256
#define OPENVDB_LOG_DEBUG_RUNTIME(message)
Log a debugging message in both debug and optimized builds.
Definition logging.h:270
void computeResidual(const MatrixOperator &A, const VecValueType *x, const VecValueType *b, VecValueType *r)
Compute r = b − Ax.
Definition ConjGradient.h:1579
void axpy(const T &a, const T *xVec, const T *yVec, T *resultVec, SizeType size)
Compute ax + y.
Definition ConjGradient.h:1560
Index32 SizeType
Definition ConjGradient.h:33
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.
Definition ConjGradient.h:1608
tbb::blocked_range< SizeType > SizeRange
Definition ConjGradient.h:35
State terminationDefaults()
Return default termination conditions for a conjugate gradient solver.
Definition ConjGradient.h:57
bool isApproxZero(const Type &x)
Return true if x is equal to zero to within the default floating-point comparison tolerance.
Definition Math.h:349
float Sqrt(float x)
Return the square root of a floating-point value.
Definition Math.h:761
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:406
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition Math.h:595
MatType scale(const Vec3< typename MatType::value_type > &s)
Return a matrix that scales by s.
Definition Mat.h:615
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition Math.h:656
Coord Abs(const Coord &xyz)
Definition Coord.h:517
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition Math.h:443
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition Math.h:337
std::ostream & operator<<(std::ostream &ostr, const Metadata &metadata)
Write a Metadata to an output stream.
Definition Metadata.h:351
uint32_t Index32
Definition Types.h:52
std::shared_ptr< T > SharedPtr
Definition Types.h:114
Definition Exceptions.h:13
#define OPENVDB_THROW(exception, message)
Definition Exceptions.h:74
Tolerance for floating-point comparison.
Definition Math.h:148
bool operator()(const SizeRange &range, bool finite) const
Definition ConjGradient.h:1023
const SparseStencilMatrix * mat
Definition ConjGradient.h:1036
IsFiniteOp(const SparseStencilMatrix &m)
Definition ConjGradient.h:1021
MatrixCopyOp(const SparseStencilMatrix &from_, SparseStencilMatrix &to_)
Definition ConjGradient.h:834
void operator()(const SizeRange &range) const
Definition ConjGradient.h:837
const SparseStencilMatrix * from
Definition ConjGradient.h:849
Information about the state of a conjugate gradient solution.
Definition ConjGradient.h:46
double absoluteError
Definition ConjGradient.h:50
int iterations
Definition ConjGradient.h:48
bool success
Definition ConjGradient.h:47
double relativeError
Definition ConjGradient.h:49
void operator()(const SizeRange &range) const
Definition ConjGradient.h:637
const SizeType binCount
Definition ConjGradient.h:657
const SizeType arraySize
Definition ConjGradient.h:658
DeterministicDotProductOp(const T *a_, const T *b_, const SizeType binCount_, const SizeType arraySize_, T *reducetmp_)
Definition ConjGradient.h:633
const T * data
Definition ConjGradient.h:717
InfNormOp(const T *data_)
Definition ConjGradient.h:707
T operator()(const SizeRange &range, T maxValue) const
Definition ConjGradient.h:709
const T * data
Definition ConjGradient.h:747
bool operator()(const SizeRange &range, bool finite) const
Definition ConjGradient.h:737
IsFiniteOp(const T *data_)
Definition ConjGradient.h:735
Definition ConjGradient.h:491
T * to
Definition ConjGradient.h:499
CopyOp(const T *from_, T *to_)
Definition ConjGradient.h:492
void operator()(const SizeRange &range) const
Definition ConjGradient.h:494
const T * from
Definition ConjGradient.h:498
Definition ConjGradient.h:506
const T val
Definition ConjGradient.h:514
void operator()(const SizeRange &range) const
Definition ConjGradient.h:509
FillOp(T *data_, const T &val_)
Definition ConjGradient.h:507
T * data
Definition ConjGradient.h:513
void operator()(const SizeRange &range) const
Definition ConjGradient.h:524
LinearOp(const T &a_, const T *x_, const T *y_, T *out_)
Definition ConjGradient.h:522
T * out
Definition ConjGradient.h:535
const T a
Definition ConjGradient.h:534
Base class for interrupters.
Definition NullInterrupter.h:26
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition version.h.in:121
#define OPENVDB_USE_VERSION_NAMESPACE
Definition version.h.in:212