OpenVDB 10.0.1
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 const ValueType sZeroValue;
249
250 /// Construct an @a n x @a n matrix with at most @a STENCIL_SIZE nonzero elements per row.
252
253 /// Deep copy the given matrix.
255
256 //@{
257 /// Return the number of rows in this matrix.
258 SizeType numRows() const { return mNumRows; }
259 SizeType size() const { return mNumRows; }
260 //@}
261
262 /// @brief Set the value at the given coordinates.
263 /// @warning It is not safe to set values in the same row simultaneously
264 /// from multiple threads.
265 void setValue(SizeType row, SizeType col, const ValueType&);
266
267 //@{
268 /// @brief Return the value at the given coordinates.
269 /// @warning It is not safe to get values from a row while another thread
270 /// is setting values in that row.
271 const ValueType& getValue(SizeType row, SizeType col) const;
272 const ValueType& operator()(SizeType row, SizeType col) const;
273 //@}
274
275 /// Return a read-only view onto the given row of this matrix.
276 ConstRow getConstRow(SizeType row) const;
277
278 /// Return a read/write view onto the given row of this matrix.
279 RowEditor getRowEditor(SizeType row);
280
281 //@{
282 /// @brief Multiply all elements in the matrix by @a s;
283 template<typename Scalar> void scale(const Scalar& s);
284 template<typename Scalar>
285 SparseStencilMatrix& operator*=(const Scalar& s) { this->scale(s); return *this; }
286 //@}
287
288 /// @brief Multiply this matrix by @a inVec and return the result in @a resultVec.
289 /// @throw ArithmeticError if either @a inVec or @a resultVec is not of size @e N,
290 /// where @e N x @e N is the size of this matrix.
291 template<typename VecValueType>
292 void vectorMultiply(const Vector<VecValueType>& inVec, Vector<VecValueType>& resultVec) const;
293
294 /// @brief Multiply this matrix by the vector represented by the array @a inVec
295 /// and return the result in @a resultVec.
296 /// @warning Both @a inVec and @a resultVec must have at least @e N elements,
297 /// where @e N x @e N is the size of this matrix.
298 template<typename VecValueType>
299 void vectorMultiply(const VecValueType* inVec, VecValueType* resultVec) const;
300
301 /// @brief Return @c true if this matrix is equivalent to the given matrix
302 /// to within the specified tolerance.
303 template<typename OtherValueType>
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>
818const ValueType SparseStencilMatrix<ValueType, STENCIL_SIZE>::sZeroValue = zeroVal<ValueType>();
819
820
821template<typename ValueType, SizeType STENCIL_SIZE>
822inline
824 : mNumRows(numRows)
825 , mValueArray(new ValueType[mNumRows * STENCIL_SIZE])
826 , mColumnIdxArray(new SizeType[mNumRows * STENCIL_SIZE])
827 , mRowSizeArray(new SizeType[mNumRows])
828{
829 // Initialize the matrix to a null state by setting the size of each row to zero.
830 tbb::parallel_for(SizeRange(0, mNumRows),
831 internal::FillOp<SizeType>(mRowSizeArray.get(), /*value=*/0));
832}
833
834
835template<typename ValueType, SizeType STENCIL_SIZE>
837{
839 from(&from_), to(&to_) {}
840
841 void operator()(const SizeRange& range) const
842 {
843 const ValueType* fromVal = from->mValueArray.get();
844 const SizeType* fromCol = from->mColumnIdxArray.get();
845 ValueType* toVal = to->mValueArray.get();
846 SizeType* toCol = to->mColumnIdxArray.get();
847 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
848 toVal[n] = fromVal[n];
849 toCol[n] = fromCol[n];
850 }
851 }
852
854};
855
856
857template<typename ValueType, SizeType STENCIL_SIZE>
858inline
860 : mNumRows(other.mNumRows)
861 , mValueArray(new ValueType[mNumRows * STENCIL_SIZE])
862 , mColumnIdxArray(new SizeType[mNumRows * STENCIL_SIZE])
863 , mRowSizeArray(new SizeType[mNumRows])
864{
865 SizeType size = mNumRows * STENCIL_SIZE;
866
867 // Copy the value and column index arrays from the other matrix to this matrix.
868 tbb::parallel_for(SizeRange(0, size), MatrixCopyOp(/*from=*/other, /*to=*/*this));
869
870 // Copy the row size array from the other matrix to this matrix.
871 tbb::parallel_for(SizeRange(0, mNumRows),
872 internal::CopyOp<SizeType>(/*from=*/other.mRowSizeArray.get(), /*to=*/mRowSizeArray.get()));
873}
874
875
876template<typename ValueType, SizeType STENCIL_SIZE>
877inline void
879 const ValueType& val)
880{
881 assert(row < mNumRows);
882 this->getRowEditor(row).setValue(col, val);
883}
884
885
886template<typename ValueType, SizeType STENCIL_SIZE>
887inline const ValueType&
889{
890 assert(row < mNumRows);
891 return this->getConstRow(row).getValue(col);
892}
893
894
895template<typename ValueType, SizeType STENCIL_SIZE>
896inline const ValueType&
898{
899 return this->getValue(row,col);
900}
901
902
903template<typename ValueType, SizeType STENCIL_SIZE>
904template<typename Scalar>
905struct SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowScaleOp
906{
907 RowScaleOp(SparseStencilMatrix& m, const Scalar& s_): mat(&m), s(s_) {}
908
909 void operator()(const SizeRange& range) const
910 {
911 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
912 RowEditor row = mat->getRowEditor(n);
913 row.scale(s);
914 }
915 }
916
917 SparseStencilMatrix* mat;
918 const Scalar s;
919};
920
921
922template<typename ValueType, SizeType STENCIL_SIZE>
923template<typename Scalar>
924inline void
926{
927 // Parallelize over the rows in the matrix.
928 tbb::parallel_for(SizeRange(0, mNumRows), RowScaleOp<Scalar>(*this, s));
929}
930
931
932template<typename ValueType, SizeType STENCIL_SIZE>
933template<typename VecValueType>
934struct SparseStencilMatrix<ValueType, STENCIL_SIZE>::VecMultOp
935{
936 VecMultOp(const SparseStencilMatrix& m, const VecValueType* i, VecValueType* o):
937 mat(&m), in(i), out(o) {}
938
939 void operator()(const SizeRange& range) const
940 {
941 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
942 ConstRow row = mat->getConstRow(n);
943 out[n] = row.dot(in, mat->numRows());
944 }
945 }
946
947 const SparseStencilMatrix* mat;
948 const VecValueType* in;
949 VecValueType* out;
950};
951
952
953template<typename ValueType, SizeType STENCIL_SIZE>
954template<typename VecValueType>
955inline void
957 const Vector<VecValueType>& inVec, Vector<VecValueType>& resultVec) const
958{
959 if (inVec.size() != mNumRows) {
960 OPENVDB_THROW(ArithmeticError, "matrix and input vector have incompatible sizes ("
961 << mNumRows << "x" << mNumRows << " vs. " << inVec.size() << ")");
962 }
963 if (resultVec.size() != mNumRows) {
964 OPENVDB_THROW(ArithmeticError, "matrix and result vector have incompatible sizes ("
965 << mNumRows << "x" << mNumRows << " vs. " << resultVec.size() << ")");
966 }
967
968 vectorMultiply(inVec.data(), resultVec.data());
969}
970
971
972template<typename ValueType, SizeType STENCIL_SIZE>
973template<typename VecValueType>
974inline void
976 const VecValueType* inVec, VecValueType* resultVec) const
977{
978 // Parallelize over the rows in the matrix.
979 tbb::parallel_for(SizeRange(0, mNumRows),
980 VecMultOp<VecValueType>(*this, inVec, resultVec));
981}
982
983
984template<typename ValueType, SizeType STENCIL_SIZE>
985template<typename OtherValueType>
986struct SparseStencilMatrix<ValueType, STENCIL_SIZE>::EqOp
987{
988 EqOp(const SparseStencilMatrix& a_,
990 a(&a_), b(&b_), eps(e) {}
991
992 bool operator()(const SizeRange& range, bool equal) const
993 {
994 if (equal) {
995 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
996 if (!a->getConstRow(n).eq(b->getConstRow(n), eps)) return false;
997 }
998 }
999 return equal;
1000 }
1001
1002 const SparseStencilMatrix* a;
1003 const SparseStencilMatrix<OtherValueType, STENCIL_SIZE>* b;
1004 const ValueType eps;
1005};
1006
1007
1008template<typename ValueType, SizeType STENCIL_SIZE>
1009template<typename OtherValueType>
1010inline bool
1013{
1014 if (this->numRows() != other.numRows()) return false;
1015 bool equal = tbb::parallel_reduce(SizeRange(0, this->numRows()), /*seed=*/true,
1016 EqOp<OtherValueType>(*this, other, eps),
1017 /*join=*/[](bool eq1, bool eq2) { return (eq1 && eq2); });
1018 return equal;
1019}
1020
1021
1022template<typename ValueType, SizeType STENCIL_SIZE>
1024{
1025 IsFiniteOp(const SparseStencilMatrix& m): mat(&m) {}
1026
1027 bool operator()(const SizeRange& range, bool finite) const
1028 {
1029 if (finite) {
1030 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1031 const ConstRow row = mat->getConstRow(n);
1032 for (ConstValueIter it = row.cbegin(); it; ++it) {
1033 if (!std::isfinite(*it)) return false;
1034 }
1035 }
1036 }
1037 return finite;
1038 }
1039
1041};
1042
1043
1044template<typename ValueType, SizeType STENCIL_SIZE>
1045inline bool
1047{
1048 // Parallelize over the rows of this matrix.
1049 bool finite = tbb::parallel_reduce(SizeRange(0, this->numRows()), /*seed=*/true,
1050 IsFiniteOp(*this), /*join=*/[](bool finite1, bool finite2) { return (finite1&&finite2); });
1051 return finite;
1052}
1053
1054
1055template<typename ValueType, SizeType STENCIL_SIZE>
1056inline std::string
1058{
1059 std::ostringstream ostr;
1060 for (SizeType n = 0, N = this->size(); n < N; ++n) {
1061 ostr << n << ": " << this->getConstRow(n).str() << "\n";
1062 }
1063 return ostr.str();
1064}
1065
1066
1067template<typename ValueType, SizeType STENCIL_SIZE>
1070{
1071 assert(i < mNumRows);
1072 const SizeType head = i * STENCIL_SIZE;
1073 return RowEditor(&mValueArray[head], &mColumnIdxArray[head], mRowSizeArray[i], mNumRows);
1074}
1075
1076
1077template<typename ValueType, SizeType STENCIL_SIZE>
1080{
1081 assert(i < mNumRows);
1082 const SizeType head = i * STENCIL_SIZE; // index for this row into main storage
1083 return ConstRow(&mValueArray[head], &mColumnIdxArray[head], mRowSizeArray[i]);
1084}
1085
1086
1087template<typename ValueType, SizeType STENCIL_SIZE>
1088template<typename DataType>
1089inline SizeType
1091{
1092 if (this->empty()) return mData.mSize;
1093
1094 // Get a pointer to the first column index that is equal to or greater than the given index.
1095 // (This assumes that the data is sorted by column.)
1096 const SizeType* colPtr = std::lower_bound(mData.mCols, mData.mCols + mData.mSize, columnIdx);
1097 // Return the offset of the pointer from the beginning of the array.
1098 return static_cast<SizeType>(colPtr - mData.mCols);
1099}
1100
1101
1102template<typename ValueType, SizeType STENCIL_SIZE>
1103template<typename DataType>
1104inline const ValueType&
1106 SizeType columnIdx, bool& active) const
1107{
1108 active = false;
1109 SizeType idx = this->find(columnIdx);
1110 if (idx < this->size() && this->column(idx) == columnIdx) {
1111 active = true;
1112 return this->value(idx);
1113 }
1115}
1116
1117template<typename ValueType, SizeType STENCIL_SIZE>
1118template<typename DataType>
1119inline const ValueType&
1121{
1122 SizeType idx = this->find(columnIdx);
1123 if (idx < this->size() && this->column(idx) == columnIdx) {
1124 return this->value(idx);
1125 }
1127}
1128
1129
1130template<typename ValueType, SizeType STENCIL_SIZE>
1131template<typename DataType>
1132inline typename SparseStencilMatrix<ValueType, STENCIL_SIZE>::ConstValueIter
1133SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::cbegin() const
1134{
1135 return ConstValueIter(mData);
1136}
1137
1138
1139template<typename ValueType, SizeType STENCIL_SIZE>
1140template<typename DataType>
1141template<typename OtherDataType>
1142inline bool
1144 const RowBase<OtherDataType>& other, ValueType eps) const
1145{
1146 if (this->size() != other.size()) return false;
1147 for (ConstValueIter it = cbegin(), oit = other.cbegin(); it || oit; ++it, ++oit) {
1148 if (it.column() != oit.column()) return false;
1149 if (!isApproxEqual(*it, *oit, eps)) return false;
1150 }
1151 return true;
1152}
1153
1154
1155template<typename ValueType, SizeType STENCIL_SIZE>
1156template<typename DataType>
1157template<typename VecValueType>
1158inline VecValueType
1159SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::dot(
1160 const VecValueType* inVec, SizeType vecSize) const
1161{
1162 VecValueType result = zeroVal<VecValueType>();
1163 for (SizeType idx = 0, N = std::min(vecSize, this->size()); idx < N; ++idx) {
1164 result += static_cast<VecValueType>(this->value(idx) * inVec[this->column(idx)]);
1165 }
1166 return result;
1167}
1168
1169template<typename ValueType, SizeType STENCIL_SIZE>
1170template<typename DataType>
1171template<typename VecValueType>
1172inline VecValueType
1173SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::dot(
1174 const Vector<VecValueType>& inVec) const
1175{
1176 return dot(inVec.data(), inVec.size());
1177}
1178
1179
1180template<typename ValueType, SizeType STENCIL_SIZE>
1181template<typename DataType>
1182inline std::string
1184{
1185 std::ostringstream ostr;
1186 std::string sep;
1187 for (SizeType n = 0, N = this->size(); n < N; ++n) {
1188 ostr << sep << "(" << this->column(n) << ", " << this->value(n) << ")";
1189 sep = ", ";
1190 }
1191 return ostr.str();
1192}
1193
1194
1195template<typename ValueType, SizeType STENCIL_SIZE>
1196inline
1198 const ValueType* valueHead, const SizeType* columnHead, const SizeType& rowSize):
1199 ConstRowBase(ConstRowData(const_cast<ValueType*>(valueHead),
1200 const_cast<SizeType*>(columnHead), const_cast<SizeType&>(rowSize)))
1201{
1202}
1203
1204
1205template<typename ValueType, SizeType STENCIL_SIZE>
1206inline
1208 ValueType* valueHead, SizeType* columnHead, SizeType& rowSize, SizeType colSize):
1209 RowBase<>(RowData(valueHead, columnHead, rowSize)), mNumColumns(colSize)
1210{
1211}
1212
1213
1214template<typename ValueType, SizeType STENCIL_SIZE>
1215inline void
1217{
1218 // Note: since mSize is a reference, this modifies the underlying matrix.
1219 RowBase<>::mData.mSize = 0;
1220}
1221
1222
1223template<typename ValueType, SizeType STENCIL_SIZE>
1224inline SizeType
1226 SizeType column, const ValueType& value)
1227{
1228 assert(column < mNumColumns);
1229
1230 RowData& data = RowBase<>::mData;
1231
1232 // Get the offset of the first column index that is equal to or greater than
1233 // the column to be modified.
1234 SizeType offset = this->find(column);
1235
1236 if (offset < data.mSize && data.mCols[offset] == column) {
1237 // If the column already exists, just update its value.
1238 data.mVals[offset] = value;
1239 return data.mSize;
1240 }
1241
1242 // Check that it is safe to add a new column.
1243 assert(data.mSize < this->capacity());
1244
1245 if (offset >= data.mSize) {
1246 // The new column's index is larger than any existing index. Append the new column.
1247 data.mVals[data.mSize] = value;
1248 data.mCols[data.mSize] = column;
1249 } else {
1250 // Insert the new column at the computed offset after shifting subsequent columns.
1251 for (SizeType i = data.mSize; i > offset; --i) {
1252 data.mVals[i] = data.mVals[i - 1];
1253 data.mCols[i] = data.mCols[i - 1];
1254 }
1255 data.mVals[offset] = value;
1256 data.mCols[offset] = column;
1257 }
1258 ++data.mSize;
1259
1260 return data.mSize;
1261}
1262
1263
1264template<typename ValueType, SizeType STENCIL_SIZE>
1265template<typename Scalar>
1266inline void
1268{
1269 for (int idx = 0, N = this->size(); idx < N; ++idx) {
1270 RowBase<>::mData.mVals[idx] *= s;
1271 }
1272}
1273
1274
1275////////////////////////////////////////
1276
1277
1278/// Diagonal preconditioner
1279template<typename MatrixType>
1280class JacobiPreconditioner: public Preconditioner<typename MatrixType::ValueType>
1281{
1282private:
1283 struct InitOp;
1284 struct ApplyOp;
1285
1286public:
1287 using ValueType = typename MatrixType::ValueType;
1291
1292 JacobiPreconditioner(const MatrixType& A): BaseType(A), mDiag(A.numRows())
1293 {
1294 // Initialize vector mDiag with the values from the matrix diagonal.
1295 tbb::parallel_for(SizeRange(0, A.numRows()), InitOp(A, mDiag.data()));
1296 }
1297
1298 ~JacobiPreconditioner() override = default;
1299
1300 void apply(const Vector<ValueType>& r, Vector<ValueType>& z) override
1301 {
1302 const SizeType size = mDiag.size();
1303
1304 assert(r.size() == z.size());
1305 assert(r.size() == size);
1306
1307 tbb::parallel_for(SizeRange(0, size), ApplyOp(mDiag.data(), r.data(), z.data()));
1308 }
1309
1310 /// Return @c true if all values along the diagonal are finite.
1311 bool isFinite() const { return mDiag.isFinite(); }
1312
1313private:
1314 // Functor for use with tbb::parallel_for()
1315 struct InitOp
1316 {
1317 InitOp(const MatrixType& m, ValueType* v): mat(&m), vec(v) {}
1318 void operator()(const SizeRange& range) const {
1319 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1320 const ValueType val = mat->getValue(n, n);
1321 assert(!isApproxZero(val, ValueType(0.0001)));
1322 vec[n] = static_cast<ValueType>(1.0 / val);
1323 }
1324 }
1325 const MatrixType* mat; ValueType* vec;
1326 };
1327
1328 // Functor for use with tbb::parallel_reduce()
1329 struct ApplyOp
1330 {
1331 ApplyOp(const ValueType* x_, const ValueType* y_, ValueType* out_):
1332 x(x_), y(y_), out(out_) {}
1333 void operator()(const SizeRange& range) const {
1334 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) out[n] = x[n] * y[n];
1335 }
1336 const ValueType *x, *y; ValueType* out;
1337 };
1338
1339 // The Jacobi preconditioner is a diagonal matrix
1340 VectorType mDiag;
1341}; // class JacobiPreconditioner
1342
1343
1344/// Preconditioner using incomplete Cholesky factorization
1345template<typename MatrixType>
1346class IncompleteCholeskyPreconditioner: public Preconditioner<typename MatrixType::ValueType>
1347{
1348private:
1349 struct CopyToLowerOp;
1350 struct TransposeOp;
1351
1352public:
1353 using ValueType = typename MatrixType::ValueType;
1358 using TriangleConstRow = typename TriangularMatrix::ConstRow;
1359 using TriangleRowEditor = typename TriangularMatrix::RowEditor;
1360
1361 IncompleteCholeskyPreconditioner(const MatrixType& matrix)
1362 : BaseType(matrix)
1363 , mLowerTriangular(matrix.numRows())
1364 , mUpperTriangular(matrix.numRows())
1365 , mTempVec(matrix.numRows())
1366 {
1367 // Size of matrix
1368 const SizeType numRows = mLowerTriangular.numRows();
1369
1370 // Copy the upper triangular part to the lower triangular part.
1371 tbb::parallel_for(SizeRange(0, numRows), CopyToLowerOp(matrix, mLowerTriangular));
1372
1373 // Build the Incomplete Cholesky Matrix
1374 //
1375 // Algorithm:
1376 //
1377 // for (k = 0; k < size; ++k) {
1378 // A(k,k) = sqrt(A(k,k));
1379 // for (i = k +1, i < size; ++i) {
1380 // if (A(i,k) == 0) continue;
1381 // A(i,k) = A(i,k) / A(k,k);
1382 // }
1383 // for (j = k+1; j < size; ++j) {
1384 // for (i = j; i < size; ++i) {
1385 // if (A(i,j) == 0) continue;
1386 // A(i,j) -= A(i,k)*A(j,k);
1387 // }
1388 // }
1389 // }
1390
1391 mPassedCompatibilityCondition = true;
1392
1393 for (SizeType k = 0; k < numRows; ++k) {
1394
1395 TriangleConstRow crow_k = mLowerTriangular.getConstRow(k);
1396 ValueType diagonalValue = crow_k.getValue(k);
1397
1398 // Test if the matrix build has failed.
1399 if (diagonalValue < 1.e-5) {
1400 mPassedCompatibilityCondition = false;
1401 break;
1402 }
1403
1404 diagonalValue = Sqrt(diagonalValue);
1405
1406 TriangleRowEditor row_k = mLowerTriangular.getRowEditor(k);
1407 row_k.setValue(k, diagonalValue);
1408
1409 // Exploit the fact that the matrix is symmetric.
1410 typename MatrixType::ConstRow srcRow = matrix.getConstRow(k);
1411 typename MatrixType::ConstValueIter citer = srcRow.cbegin();
1412 for ( ; citer; ++citer) {
1413 SizeType ii = citer.column();
1414 if (ii < k+1) continue; // look above diagonal
1415
1416 TriangleRowEditor row_ii = mLowerTriangular.getRowEditor(ii);
1417
1418 row_ii.setValue(k, *citer / diagonalValue);
1419 }
1420
1421 // for (j = k+1; j < size; ++j) replaced by row iter below
1422 citer.reset(); // k,j entries
1423 for ( ; citer; ++citer) {
1424 SizeType j = citer.column();
1425 if (j < k+1) continue;
1426
1427 TriangleConstRow row_j = mLowerTriangular.getConstRow(j);
1428 ValueType a_jk = row_j.getValue(k); // a_jk is non zero if a_kj is non zero
1429
1430 // Entry (i,j) is non-zero if matrix(j,i) is nonzero
1431
1432 typename MatrixType::ConstRow mask = matrix.getConstRow(j);
1433 typename MatrixType::ConstValueIter maskIter = mask.cbegin();
1434 for ( ; maskIter; ++maskIter) {
1435 SizeType i = maskIter.column();
1436 if (i < j) continue;
1437
1438 TriangleConstRow crow_i = mLowerTriangular.getConstRow(i);
1439 ValueType a_ij = crow_i.getValue(j);
1440 ValueType a_ik = crow_i.getValue(k);
1441 TriangleRowEditor row_i = mLowerTriangular.getRowEditor(i);
1442 a_ij -= a_ik * a_jk;
1443
1444 row_i.setValue(j, a_ij);
1445 }
1446 }
1447 }
1448
1449 // Build the transpose of the IC matrix: mUpperTriangular
1450 tbb::parallel_for(SizeRange(0, numRows),
1451 TransposeOp(matrix, mLowerTriangular, mUpperTriangular));
1452 }
1453
1455
1456 bool isValid() const override { return mPassedCompatibilityCondition; }
1457
1458 void apply(const Vector<ValueType>& rVec, Vector<ValueType>& zVec) override
1459 {
1460 if (!mPassedCompatibilityCondition) {
1461 OPENVDB_THROW(ArithmeticError, "invalid Cholesky decomposition");
1462 }
1463
1464 // Solve mUpperTriangular * mLowerTriangular * rVec = zVec;
1465
1466 SizeType size = mLowerTriangular.numRows();
1467
1468 zVec.fill(zeroVal<ValueType>());
1469 ValueType* zData = zVec.data();
1470
1471 if (size == 0) return;
1472
1473 assert(rVec.size() == size);
1474 assert(zVec.size() == size);
1475
1476 // Allocate a temp vector
1477 mTempVec.fill(zeroVal<ValueType>());
1478 ValueType* tmpData = mTempVec.data();
1479 const ValueType* rData = rVec.data();
1480
1481 // Solve mLowerTriangular * tmp = rVec;
1482 for (SizeType i = 0; i < size; ++i) {
1483 typename TriangularMatrix::ConstRow row = mLowerTriangular.getConstRow(i);
1484 ValueType diagonal = row.getValue(i);
1485 ValueType dot = row.dot(mTempVec);
1486 tmpData[i] = (rData[i] - dot) / diagonal;
1487 if (!std::isfinite(tmpData[i])) {
1488 OPENVDB_LOG_DEBUG_RUNTIME("1 diagonal was " << diagonal);
1489 OPENVDB_LOG_DEBUG_RUNTIME("1a diagonal " << row.getValue(i));
1490 }
1491 }
1492
1493 // Solve mUpperTriangular * zVec = tmp;
1494 for (SizeType ii = 0; ii < size; ++ii) {
1495 SizeType i = size - 1 - ii;
1496 typename TriangularMatrix::ConstRow row = mUpperTriangular.getConstRow(i);
1497 ValueType diagonal = row.getValue(i);
1498 ValueType dot = row.dot(zVec);
1499 zData[i] = (tmpData[i] - dot) / diagonal;
1500 if (!std::isfinite(zData[i])) {
1501 OPENVDB_LOG_DEBUG_RUNTIME("2 diagonal was " << diagonal);
1502 }
1503 }
1504 }
1505
1506 const TriangularMatrix& lowerMatrix() const { return mLowerTriangular; }
1507 const TriangularMatrix& upperMatrix() const { return mUpperTriangular; }
1508
1509private:
1510 // Functor for use with tbb::parallel_for()
1511 struct CopyToLowerOp
1512 {
1513 CopyToLowerOp(const MatrixType& m, TriangularMatrix& l): mat(&m), lower(&l) {}
1514 void operator()(const SizeRange& range) const {
1515 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1516 typename TriangularMatrix::RowEditor outRow = lower->getRowEditor(n);
1517 outRow.clear();
1518 typename MatrixType::ConstRow inRow = mat->getConstRow(n);
1519 for (typename MatrixType::ConstValueIter it = inRow.cbegin(); it; ++it) {
1520 if (it.column() > n) continue; // skip above diagonal
1521 outRow.setValue(it.column(), *it);
1522 }
1523 }
1524 }
1525 const MatrixType* mat; TriangularMatrix* lower;
1526 };
1527
1528 // Functor for use with tbb::parallel_for()
1529 struct TransposeOp
1530 {
1531 TransposeOp(const MatrixType& m, const TriangularMatrix& l, TriangularMatrix& u):
1532 mat(&m), lower(&l), upper(&u) {}
1533 void operator()(const SizeRange& range) const {
1534 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1535 typename TriangularMatrix::RowEditor outRow = upper->getRowEditor(n);
1536 outRow.clear();
1537 // Use the fact that matrix is symmetric.
1538 typename MatrixType::ConstRow inRow = mat->getConstRow(n);
1539 for (typename MatrixType::ConstValueIter it = inRow.cbegin(); it; ++it) {
1540 const SizeType column = it.column();
1541 if (column < n) continue; // only set upper triangle
1542 outRow.setValue(column, lower->getValue(column, n));
1543 }
1544 }
1545 }
1546 const MatrixType* mat; const TriangularMatrix* lower; TriangularMatrix* upper;
1547 };
1548
1549 TriangularMatrix mLowerTriangular;
1550 TriangularMatrix mUpperTriangular;
1551 Vector<ValueType> mTempVec;
1552 bool mPassedCompatibilityCondition;
1553}; // class IncompleteCholeskyPreconditioner
1554
1555
1556////////////////////////////////////////
1557
1558
1559namespace internal {
1560
1561/// Compute @e ax + @e y.
1562template<typename T>
1563inline void
1564axpy(const T& a, const T* xVec, const T* yVec, T* resultVec, SizeType size)
1565{
1566 tbb::parallel_for(SizeRange(0, size), LinearOp<T>(a, xVec, yVec, resultVec));
1567}
1568
1569/// Compute @e ax + @e y.
1570template<typename T>
1571inline void
1572axpy(const T& a, const Vector<T>& xVec, const Vector<T>& yVec, Vector<T>& result)
1573{
1574 assert(xVec.size() == yVec.size());
1575 assert(xVec.size() == result.size());
1576 axpy(a, xVec.data(), yVec.data(), result.data(), xVec.size());
1577}
1578
1579
1580/// Compute @e r = @e b &minus; @e Ax.
1581template<typename MatrixOperator, typename VecValueType>
1582inline void
1583computeResidual(const MatrixOperator& A, const VecValueType* x,
1584 const VecValueType* b, VecValueType* r)
1585{
1586 // Compute r = A * x.
1587 A.vectorMultiply(x, r);
1588 // Compute r = b - r.
1589 tbb::parallel_for(SizeRange(0, A.numRows()), LinearOp<VecValueType>(-1.0, r, b, r));
1590}
1591
1592/// Compute @e r = @e b &minus; @e Ax.
1593template<typename MatrixOperator, typename T>
1594inline void
1595computeResidual(const MatrixOperator& A, const Vector<T>& x, const Vector<T>& b, Vector<T>& r)
1596{
1597 assert(x.size() == b.size());
1598 assert(x.size() == r.size());
1599 assert(x.size() == A.numRows());
1600
1601 computeResidual(A, x.data(), b.data(), r.data());
1602}
1603
1604} // namespace internal
1605
1606
1607////////////////////////////////////////
1608
1609
1610template<typename PositiveDefMatrix>
1611inline State
1613 const PositiveDefMatrix& Amat,
1617 const State& termination)
1618{
1619 util::NullInterrupter interrupter;
1620 return solve(Amat, bVec, xVec, precond, interrupter, termination);
1621}
1622
1623
1624template<typename PositiveDefMatrix, typename Interrupter>
1625inline State
1627 const PositiveDefMatrix& Amat,
1631 Interrupter& interrupter,
1632 const State& termination)
1633{
1634 using ValueType = typename PositiveDefMatrix::ValueType;
1636
1637 State result;
1638 result.success = false;
1639 result.iterations = 0;
1640 result.relativeError = 0.0;
1641 result.absoluteError = 0.0;
1642
1643 const SizeType size = Amat.numRows();
1644 if (size == 0) {
1645 OPENVDB_LOG_WARN("pcg::solve(): matrix has dimension zero");
1646 return result;
1647 }
1648 if (size != bVec.size()) {
1649 OPENVDB_THROW(ArithmeticError, "A and b have incompatible sizes"
1650 << size << "x" << size << " vs. " << bVec.size() << ")");
1651 }
1652 if (size != xVec.size()) {
1653 OPENVDB_THROW(ArithmeticError, "A and x have incompatible sizes"
1654 << size << "x" << size << " vs. " << xVec.size() << ")");
1655 }
1656
1657 // Temp vectors
1658 VectorType zVec(size); // transformed residual (M^-1 r)
1659 VectorType pVec(size); // search direction
1660 VectorType qVec(size); // A * p
1661
1662 // Compute norm of B (the source)
1663 const ValueType tmp = bVec.infNorm();
1664 const ValueType infNormOfB = isZero(tmp) ? ValueType(1) : tmp;
1665
1666 // Compute rVec: residual = b - Ax.
1667 VectorType rVec(size); // vector of residuals
1668
1669 internal::computeResidual(Amat, xVec, bVec, rVec);
1670
1671 assert(rVec.isFinite());
1672
1673 // Normalize the residual norm with the source norm and look for early out.
1674 result.absoluteError = static_cast<double>(rVec.infNorm());
1675 result.relativeError = static_cast<double>(result.absoluteError / infNormOfB);
1676 if (result.relativeError <= termination.relativeError) {
1677 result.success = true;
1678 return result;
1679 }
1680
1681 // Iterations of the CG solve
1682
1683 ValueType rDotZPrev(1); // inner product of <z,r>
1684
1685 // Keep track of the minimum error to monitor convergence.
1686 ValueType minL2Error = std::numeric_limits<ValueType>::max();
1687 ValueType l2Error;
1688
1689 int iteration = 0;
1690 for ( ; iteration < termination.iterations; ++iteration) {
1691
1692 if (interrupter.wasInterrupted()) {
1693 OPENVDB_THROW(RuntimeError, "conjugate gradient solver was interrupted");
1694 }
1695
1696 OPENVDB_LOG_DEBUG_RUNTIME("pcg::solve() " << result);
1697
1698 result.iterations = iteration + 1;
1699
1700 // Apply preconditioner to residual
1701 // z_{k} = M^-1 r_{k}
1702 precond.apply(rVec, zVec);
1703
1704 // <r,z>
1705 const ValueType rDotZ = rVec.dot(zVec);
1706 assert(std::isfinite(rDotZ));
1707
1708 if (0 == iteration) {
1709 // Initialize
1710 pVec = zVec;
1711 } else {
1712 const ValueType beta = rDotZ / rDotZPrev;
1713 // p = beta * p + z
1714 internal::axpy(beta, pVec, zVec, /*result */pVec);
1715 }
1716
1717 // q_{k} = A p_{k}
1718 Amat.vectorMultiply(pVec, qVec);
1719
1720 // alpha = <r_{k-1}, z_{k-1}> / <p_{k},q_{k}>
1721 const ValueType pAp = pVec.dot(qVec);
1722 assert(std::isfinite(pAp));
1723
1724 const ValueType alpha = rDotZ / pAp;
1725 rDotZPrev = rDotZ;
1726
1727 // x_{k} = x_{k-1} + alpha * p_{k}
1728 internal::axpy(alpha, pVec, xVec, /*result=*/xVec);
1729
1730 // r_{k} = r_{k-1} - alpha_{k-1} A p_{k}
1731 internal::axpy(-alpha, qVec, rVec, /*result=*/rVec);
1732
1733 // update tolerances
1734 l2Error = rVec.l2Norm();
1735 minL2Error = Min(l2Error, minL2Error);
1736
1737 result.absoluteError = static_cast<double>(rVec.infNorm());
1738 result.relativeError = static_cast<double>(result.absoluteError / infNormOfB);
1739
1740 if (l2Error > 2 * minL2Error) {
1741 // The solution started to diverge.
1742 result.success = false;
1743 break;
1744 }
1745 if (!std::isfinite(result.absoluteError)) {
1746 // Total divergence of solution
1747 result.success = false;
1748 break;
1749 }
1750 if (result.absoluteError <= termination.absoluteError) {
1751 // Convergence
1752 result.success = true;
1753 break;
1754 }
1755 if (result.relativeError <= termination.relativeError) {
1756 // Convergence
1757 result.success = true;
1758 break;
1759 }
1760 }
1761 OPENVDB_LOG_DEBUG_RUNTIME("pcg::solve() " << result);
1762
1763 return result;
1764}
1765
1766} // namespace pcg
1767} // namespace math
1768} // namespace OPENVDB_VERSION_NAME
1769} // namespace openvdb
1770
1771#endif // OPENVDB_MATH_CONJGRADIENT_HAS_BEEN_INCLUDED
ValueT value
Definition GridBuilder.h:1290
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:1347
typename TriangularMatrix::RowEditor TriangleRowEditor
Definition ConjGradient.h:1359
SharedPtr< IncompleteCholeskyPreconditioner > Ptr
Definition ConjGradient.h:1356
IncompleteCholeskyPreconditioner(const MatrixType &matrix)
Definition ConjGradient.h:1361
typename TriangularMatrix::ConstRow TriangleConstRow
Definition ConjGradient.h:1358
void apply(const Vector< ValueType > &rVec, Vector< ValueType > &zVec) override
Apply this preconditioner to a residue vector: z = M−1r
Definition ConjGradient.h:1458
typename MatrixType::ValueType ValueType
Definition ConjGradient.h:1353
const TriangularMatrix & upperMatrix() const
Definition ConjGradient.h:1507
bool isValid() const override
Definition ConjGradient.h:1456
const TriangularMatrix & lowerMatrix() const
Definition ConjGradient.h:1506
Diagonal preconditioner.
Definition ConjGradient.h:1281
SharedPtr< JacobiPreconditioner > Ptr
Definition ConjGradient.h:1290
void apply(const Vector< ValueType > &r, Vector< ValueType > &z) override
Apply this preconditioner to a residue vector: z = M−1r
Definition ConjGradient.h:1300
typename MatrixType::ValueType ValueType
Definition ConjGradient.h:1287
JacobiPreconditioner(const MatrixType &A)
Definition ConjGradient.h:1292
bool isFinite() const
Return true if all values along the diagonal are finite.
Definition ConjGradient.h:1311
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:1197
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:1225
RowEditor & operator*=(const Scalar &s)
Definition ConjGradient.h:435
void scale(const Scalar &)
Scale all of the entries in this row.
Definition ConjGradient.h:1267
void clear()
Set the number of entries in this row to zero.
Definition ConjGradient.h:1216
RowEditor(ValueType *valueHead, SizeType *columnHead, SizeType &rowSize, SizeType colSize)
Definition ConjGradient.h:1207
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:888
SparseStencilMatrix(const SparseStencilMatrix &)
Deep copy the given matrix.
Definition ConjGradient.h:859
void setValue(SizeType row, SizeType col, const ValueType &)
Set the value at the given coordinates.
Definition ConjGradient.h:878
ConstRow getConstRow(SizeType row) const
Return a read-only view onto the given row of this matrix.
Definition ConjGradient.h:1079
SparseStencilMatrix(SizeType n)
Construct an n x n matrix with at most STENCIL_SIZE nonzero elements per row.
Definition ConjGradient.h:823
const ValueType & operator()(SizeType row, SizeType col) const
Definition ConjGradient.h:897
RowEditor getRowEditor(SizeType row)
Return a read/write view onto the given row of this matrix.
Definition ConjGradient.h:1069
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:956
SparseStencilMatrix & operator*=(const Scalar &s)
Definition ConjGradient.h:285
static const ValueType sZeroValue
Definition ConjGradient.h:248
void scale(const Scalar &s)
Multiply all elements in the matrix by s;.
Definition ConjGradient.h:925
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:975
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:1011
std::string str() const
Return a string representation of this matrix.
Definition ConjGradient.h:1057
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:1046
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:1583
void axpy(const T &a, const T *xVec, const T *yVec, T *resultVec, SizeType size)
Compute ax + y.
Definition ConjGradient.h:1564
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:1612
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:515
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:1027
const SparseStencilMatrix * mat
Definition ConjGradient.h:1040
IsFiniteOp(const SparseStencilMatrix &m)
Definition ConjGradient.h:1025
MatrixCopyOp(const SparseStencilMatrix &from_, SparseStencilMatrix &to_)
Definition ConjGradient.h:838
void operator()(const SizeRange &range) const
Definition ConjGradient.h:841
const SparseStencilMatrix * from
Definition ConjGradient.h:853
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