GCC Code Coverage Report

Directory: ./
File: openvdb/openvdb/math/ConjGradient.h
Date: 2022-07-25 17:40:05
Exec Total Coverage
Lines: 361 382 94.5%
Functions: 54 56 96.4%
Branches: 302 631 47.9%
Line Branch Exec Source
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3
4 /// @file ConjGradient.h
5 /// @authors D.J. Hill, Peter Cucka
6 /// @brief Preconditioned conjugate gradient solver (solves @e Ax = @e b using
7 /// the conjugate gradient method with one of a selection of preconditioners)
8
9 #ifndef OPENVDB_MATH_CONJGRADIENT_HAS_BEEN_INCLUDED
10 #define OPENVDB_MATH_CONJGRADIENT_HAS_BEEN_INCLUDED
11
12 #include <openvdb/Exceptions.h>
13 #include <openvdb/Types.h>
14 #include <openvdb/util/logging.h>
15 #include <openvdb/util/NullInterrupter.h>
16 #include "Math.h" // for Abs(), isZero(), Max(), Sqrt()
17 #include <tbb/parallel_for.h>
18 #include <tbb/parallel_reduce.h>
19 #include <algorithm> // for std::lower_bound()
20 #include <cassert>
21 #include <cmath> // for std::isfinite()
22 #include <limits>
23 #include <sstream>
24 #include <string>
25
26
27 namespace openvdb {
28 OPENVDB_USE_VERSION_NAMESPACE
29 namespace OPENVDB_VERSION_NAME {
30 namespace math {
31 namespace pcg {
32
33 using SizeType = Index32;
34
35 using SizeRange = tbb::blocked_range<SizeType>;
36
37 template<typename ValueType> class Vector;
38
39 template<typename ValueType, SizeType STENCIL_SIZE> class SparseStencilMatrix;
40
41 template<typename ValueType> class Preconditioner;
42 template<typename MatrixType> class JacobiPreconditioner;
43 template<typename MatrixType> class IncompleteCholeskyPreconditioner;
44
45 /// Information about the state of a conjugate gradient solution
46 struct State {
47 bool success;
48 int iterations;
49 double relativeError;
50 double absoluteError;
51 };
52
53
54 /// Return default termination conditions for a conjugate gradient solver.
55 template<typename ValueType>
56 inline State
57 terminationDefaults()
58 {
59 State s;
60
4/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
 9 s.success = false;
61 2 s.iterations = 50;
62
3/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
 5 s.relativeError = 1.0e-6;
63 2 s.absoluteError = std::numeric_limits<ValueType>::epsilon() * 100.0;
64 return s;
65 }
66
67
68 ////////////////////////////////////////
69
70
71 /// @brief Solve @e Ax = @e b via the preconditioned conjugate gradient method.
72 ///
73 /// @param A a symmetric, positive-definite, @e N x @e N matrix
74 /// @param b a vector of size @e N
75 /// @param x a vector of size @e N
76 /// @param preconditioner a Preconditioner matrix
77 /// @param termination termination conditions given as a State object with the following fields:
78 /// <dl>
79 /// <dt><i>success</i>
80 /// <dd>ignored
81 /// <dt><i>iterations</i>
82 /// <dd>the maximum number of iterations, with or without convergence
83 /// <dt><i>relativeError</i>
84 /// <dd>the relative error ||<i>b</i> &minus; <i>Ax</i>|| / ||<i>b</i>||
85 /// that denotes convergence
86 /// <dt><i>absoluteError</i>
87 /// <dd>the absolute error ||<i>b</i> &minus; <i>Ax</i>|| that denotes convergence
88 ///
89 /// @throw ArithmeticError if either @a x or @a b is not of the appropriate size.
90 template<typename PositiveDefMatrix>
91 inline State
92 solve(
93 const PositiveDefMatrix& A,
94 const Vector<typename PositiveDefMatrix::ValueType>& b,
95 Vector<typename PositiveDefMatrix::ValueType>& x,
96 Preconditioner<typename PositiveDefMatrix::ValueType>& preconditioner,
97 const State& termination = terminationDefaults<typename PositiveDefMatrix::ValueType>());
98
99
100 /// @brief Solve @e Ax = @e b via the preconditioned conjugate gradient method.
101 ///
102 /// @param A a symmetric, positive-definite, @e N x @e N matrix
103 /// @param b a vector of size @e N
104 /// @param x a vector of size @e N
105 /// @param preconditioner a Preconditioner matrix
106 /// @param termination termination conditions given as a State object with the following fields:
107 /// <dl>
108 /// <dt><i>success</i>
109 /// <dd>ignored
110 /// <dt><i>iterations</i>
111 /// <dd>the maximum number of iterations, with or without convergence
112 /// <dt><i>relativeError</i>
113 /// <dd>the relative error ||<i>b</i> &minus; <i>Ax</i>|| / ||<i>b</i>||
114 /// that denotes convergence
115 /// <dt><i>absoluteError</i>
116 /// <dd>the absolute error ||<i>b</i> &minus; <i>Ax</i>|| that denotes convergence
117 /// @param interrupter an object adhering to the util::NullInterrupter interface
118 /// with which computation can be interrupted
119 ///
120 /// @throw ArithmeticError if either @a x or @a b is not of the appropriate size.
121 /// @throw RuntimeError if the computation is interrupted.
122 template<typename PositiveDefMatrix, typename Interrupter>
123 inline State
124 solve(
125 const PositiveDefMatrix& A,
126 const Vector<typename PositiveDefMatrix::ValueType>& b,
127 Vector<typename PositiveDefMatrix::ValueType>& x,
128 Preconditioner<typename PositiveDefMatrix::ValueType>& preconditioner,
129 Interrupter& interrupter,
130 const State& termination = terminationDefaults<typename PositiveDefMatrix::ValueType>());
131
132
133 ////////////////////////////////////////
134
135
136 /// Lightweight, variable-length vector
137 template<typename T>
138 class Vector
139 {
140 public:
141 using ValueType = T;
142 using Ptr = SharedPtr<Vector>;
143
144 /// Construct an empty vector.
145 Vector(): mData(nullptr), mSize(0) {}
146 /// Construct a vector of @a n elements, with uninitialized values.
147
12/24
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 9 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 11 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 9 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✓ Branch 14 taken 1 times.
✓ Branch 15 taken 3 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 1 times.
✓ Branch 20 taken 3 times.
✗ Branch 21 not taken.
✓ Branch 23 taken 3 times.
✗ Branch 24 not taken.
✓ Branch 27 taken 2 times.
✗ Branch 28 not taken.
 34 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 15 Vector(SizeType n, const ValueType& val): mData(new T[n]), mSize(n) { this->fill(val); }
150
151
18/72
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 9 times.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 9 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 9 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 9 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✓ Branch 27 taken 1 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✓ Branch 30 taken 1 times.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✓ Branch 33 taken 1 times.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✓ Branch 36 taken 1 times.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 42 not taken.
✓ Branch 43 taken 2 times.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✓ Branch 49 taken 2 times.
✗ Branch 50 not taken.
✓ Branch 51 taken 1 times.
✓ Branch 52 taken 2 times.
✗ Branch 53 not taken.
✓ Branch 54 taken 1 times.
✗ Branch 55 not taken.
✓ Branch 57 taken 1 times.
✗ Branch 58 not taken.
✗ Branch 61 not taken.
✗ Branch 62 not taken.
✗ Branch 64 not taken.
✗ Branch 65 not taken.
✗ Branch 67 not taken.
✗ Branch 68 not taken.
✓ Branch 72 taken 1 times.
✗ Branch 73 not taken.
✓ Branch 75 taken 1 times.
✗ Branch 76 not taken.
✗ Branch 80 not taken.
✗ Branch 81 not taken.
✗ Branch 83 not taken.
✗ Branch 84 not taken.
✗ Branch 86 not taken.
✗ Branch 87 not taken.
 43 ~Vector() { mSize = 0; delete[] mData; mData = nullptr; }
152
153 /// Deep copy the given vector.
154 Vector(const Vector&);
155 /// Deep copy the given vector.
156 Vector& operator=(const Vector&);
157
158 /// Return the number of elements in this vector.
159
9/27
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 16 taken 2 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 2 times.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✓ Branch 28 taken 2 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 2 times.
✗ Branch 32 not taken.
✓ Branch 33 taken 7202 times.
✓ Branch 34 taken 3 times.
✗ Branch 35 not taken.
✓ Branch 36 taken 1 times.
✗ Branch 37 not taken.
 7216 SizeType size() const { return mSize; }
160 /// Return @c true if this vector has no elements.
161 bool empty() const { return (mSize == 0); }
162
163 /// @brief Reset this vector to have @a n elements, with uninitialized values.
164 /// @warning All of this vector's existing values will be lost.
165 void resize(SizeType n);
166
167 /// Swap internal storage with another vector, which need not be the same size.
168 void swap(Vector& other) { std::swap(mData, other.mData); std::swap(mSize, other.mSize); }
169
170 /// Set all elements of this vector to @a value.
171 void fill(const ValueType& value);
172
173 //@{
174 /// @brief Multiply each element of this vector by @a s.
175 template<typename Scalar> void scale(const Scalar& s);
176 template<typename Scalar> Vector& operator*=(const Scalar& s) { this->scale(s); return *this; }
177 //@}
178
179 /// Return the dot product of this vector with the given vector, which must be the same size.
180 ValueType dot(const Vector&) const;
181
182 /// Return the infinity norm of this vector.
183 ValueType infNorm() const;
184 /// Return the L2 norm of this vector.
185 471 ValueType l2Norm() const { return Sqrt(this->dot(*this)); }
186
187 /// Return @c true if every element of this vector has a finite value.
188 bool isFinite() const;
189
190 /// @brief Return @c true if this vector is equivalent to the given vector
191 /// to within the specified tolerance.
192 template<typename OtherValueType>
193 bool eq(const Vector<OtherValueType>& other,
194 ValueType eps = Tolerance<ValueType>::value()) const;
195
196 /// Return a string representation of this vector.
197 std::string str() const;
198
199 //@{
200 /// @brief Return the value of this vector's ith element.
201
20/50
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 1457 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 1457 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✓ Branch 14 taken 1922 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 1922 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✓ Branch 20 taken 1457 times.
✗ Branch 21 not taken.
✓ Branch 23 taken 1457 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✓ Branch 27 taken 22418 times.
✗ Branch 28 not taken.
✓ Branch 30 taken 29620 times.
✗ Branch 31 not taken.
✓ Branch 33 taken 22418 times.
✗ Branch 34 not taken.
✓ Branch 36 taken 22418 times.
✗ Branch 37 not taken.
✓ Branch 39 taken 22418 times.
✗ Branch 40 not taken.
✓ Branch 42 taken 22418 times.
✗ Branch 43 not taken.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✓ Branch 50 taken 12418 times.
✗ Branch 51 not taken.
✓ Branch 53 taken 12418 times.
✗ Branch 54 not taken.
✓ Branch 56 taken 12418 times.
✗ Branch 57 not taken.
✓ Branch 59 taken 12418 times.
✗ Branch 60 not taken.
✓ Branch 62 taken 12418 times.
✗ Branch 63 not taken.
✓ Branch 65 taken 12418 times.
✗ Branch 66 not taken.
 4078523 inline T& at(SizeType i) { return mData[i]; }
202
1/2
✓ Branch 4 taken 7202 times.
✗ Branch 5 not taken.
 3852433 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 2825 inline T* data() { return mData; }
210 708006416 inline const T* data() const { return mData; }
211 inline const T* constData() const { return mData; }
212 //@}
213
214 private:
215 // Functor for use with tbb::parallel_for()
216 template<typename Scalar> struct ScaleOp;
217 struct DeterministicDotProductOp;
218 // Functors for use with tbb::parallel_reduce()
219 template<typename OtherValueType> struct EqOp;
220 struct InfNormOp;
221 struct IsFiniteOp;
222
223 T* mData;
224 SizeType mSize;
225 };
226
227 using VectorS = Vector<float>;
228 using VectorD = Vector<double>;
229
230
231 ////////////////////////////////////////
232
233
234 /// @brief Sparse, square matrix representing a 3D stencil operator of size @a STENCIL_SIZE
235 /// @details The implementation is a variation on compressed row storage (CRS).
236 template<typename ValueType_, SizeType STENCIL_SIZE>
237 class SparseStencilMatrix
238 {
239 public:
240 using ValueType = ValueType_;
241 using VectorType = Vector<ValueType>;
242 using Ptr = SharedPtr<SparseStencilMatrix>;
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.
251 SparseStencilMatrix(SizeType n);
252
253 /// Deep copy the given matrix.
254 SparseStencilMatrix(const SparseStencilMatrix&);
255
256 //@{
257 /// Return the number of rows in this matrix.
258
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 SizeType numRows() const { return mNumRows; }
259
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 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>
304 bool eq(const SparseStencilMatrix<OtherValueType, STENCIL_SIZE>& other,
305 ValueType eps = Tolerance<ValueType>::value()) const;
306
307 /// Return @c true if every element of this matrix has a finite value.
308 bool isFinite() const;
309
310 /// Return a string representation of this matrix.
311 std::string str() const;
312
313 private:
314 struct RowData {
315 70502483 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 1157260973 ConstRowData(const ValueType* v, const SizeType* c, const SizeType& s):
321 1157260973 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 1227763456 RowBase(const DataType& data): mData(data) {}
335
336 929686646 bool empty() const { return (mData.mSize == 0); }
337 1889530866 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 6075874694 const ValueType& value(SizeType i) const { return mData.mVals[i]; }
368 6107290272 SizeType column(SizeType i) const { return mData.mCols[i]; }
369
370 /// @brief Return the array index of the first column index that is
371 /// equal to <i>or greater than</i> the given column index.
372 /// @note If @a columnIdx is larger than any existing column index,
373 /// the return value will point beyond the end of the array.
374 SizeType find(SizeType columnIdx) const;
375
376 DataType mData;
377 };
378
379 using ConstRowBase = RowBase<ConstRowData>;
380
381 public:
382 /// Iterator over the stored values in a row of this matrix
383 class ConstValueIter
384 {
385 public:
386 const ValueType& operator*() const
387 {
388
5/10
✓ Branch 0 taken 18 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 18 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 26 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 14169979 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 10592476 times.
✗ Branch 9 not taken.
 24762517 if (mData.mSize == 0) return SparseStencilMatrix::sZeroValue;
389 24762517 return mData.mVals[mCursor];
390 }
391
392 172613470 SizeType column() const { return mData.mCols[mCursor]; }
393
394 172613496 void increment() { mCursor++; }
395 172613496 ConstValueIter& operator++() { increment(); return *this; }
396 197516014 operator bool() const { return (mCursor < mData.mSize); }
397
398 3577503 void reset() { mCursor = 0; }
399
400 private:
401 friend class SparseStencilMatrix;
402 ConstValueIter(const RowData& d): mData(d.mVals, d.mCols, d.mSize), mCursor(0) {}
403 21325005 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
442 private:
443 // Functors for use with tbb::parallel_for()
444 struct MatrixCopyOp;
445 template<typename VecValueType> struct VecMultOp;
446 template<typename Scalar> struct RowScaleOp;
447
448 // Functors for use with tbb::parallel_reduce()
449 struct IsFiniteOp;
450 template<typename OtherValueType> struct EqOp;
451
452 const SizeType mNumRows;
453 std::unique_ptr<ValueType[]> mValueArray;
454 std::unique_ptr<SizeType[]> mColumnIdxArray;
455 std::unique_ptr<SizeType[]> mRowSizeArray;
456 }; // class SparseStencilMatrix
457
458
459 ////////////////////////////////////////
460
461
462 /// Base class for conjugate gradient preconditioners
463 template<typename T>
464 class Preconditioner
465 {
466 public:
467 using ValueType = T;
468 using Ptr = SharedPtr<Preconditioner>;
469
470 9 template<SizeType STENCIL_SIZE> Preconditioner(const SparseStencilMatrix<T, STENCIL_SIZE>&) {}
471 virtual ~Preconditioner() = default;
472
473 ✗ virtual bool isValid() const { return true; }
474
475 /// @brief Apply this preconditioner to a residue vector:
476 /// @e z = <i>M</i><sup><small>&minus;1</small></sup><i>r</i>
477 /// @param r residue vector
478 /// @param[out] z preconditioned residue vector
479 virtual void apply(const Vector<T>& r, Vector<T>& z) = 0;
480 };
481
482
483 ////////////////////////////////////////
484
485
486 namespace internal {
487
488 // Functor for use with tbb::parallel_for() to copy data from one array to another
489 template<typename T>
490 struct CopyOp
491 {
492 11 CopyOp(const T* from_, T* to_): from(from_), to(to_) {}
493
494 void operator()(const SizeRange& range) const {
495
6/8
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 208424 times.
✓ Branch 5 taken 90 times.
✓ Branch 6 taken 3369084 times.
✓ Branch 7 taken 922 times.
 3578540 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) to[n] = from[n];
496 }
497
498 const T* from;
499 T* to;
500 };
501
502
503 // Functor for use with tbb::parallel_for() to fill an array with a constant value
504 template<typename T>
505 struct FillOp
506 {
507 980 FillOp(T* data_, const T& val_): data(data_), val(val_) {}
508
509 void operator()(const SizeRange& range) const {
510
8/8
✓ Branch 0 taken 628072 times.
✓ Branch 1 taken 263 times.
✓ Branch 2 taken 10111654 times.
✓ Branch 3 taken 2899 times.
✓ Branch 4 taken 43822536 times.
✓ Branch 5 taken 9418 times.
✓ Branch 6 taken 687828130 times.
✓ Branch 7 taken 119409 times.
 742522381 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) data[n] = val;
511 }
512
513 T* data;
514 const T val;
515 };
516
517
518 // Functor for use with tbb::parallel_for() that computes a * x + y
519 template<typename T>
520 struct LinearOp
521 {
522 1404 LinearOp(const T& a_, const T* x_, const T* y_, T* out_): a(a_), x(x_), y(y_), out(out_) {}
523
524
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 261059 times.
 261059 void operator()(const SizeRange& range) const {
525
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 261059 times.
 261059 if (isExactlyEqual(a, T(1))) {
526 ✗ for (SizeType n = range.begin(), N = range.end(); n < N; ++n) out[n] = x[n] + y[n];
527
2/2
✓ Branch 0 taken 1318 times.
✓ Branch 1 taken 259741 times.
 261059 } else if (isExactlyEqual(a, T(-1))) {
528
2/2
✓ Branch 0 taken 3577508 times.
✓ Branch 1 taken 1318 times.
 3578826 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) out[n] = -x[n] + y[n];
529 } else {
530
2/2
✓ Branch 0 taken 1058425753 times.
✓ Branch 1 taken 259741 times.
 1058685494 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) out[n] = a * x[n] + y[n];
531 }
532 261059 }
533
534 const T a, *x, *y;
535 T* out;
536 };
537
538 } // namespace internal
539
540
541 ////////////////////////////////////////
542
543
544 inline std::ostream&
545 operator<<(std::ostream& os, const State& state)
546 {
547 os << (state.success ? "succeeded with " : "")
548 << "rel. err. " << state.relativeError << ", abs. err. " << state.absoluteError
549 << " after " << state.iterations << " iteration" << (state.iterations == 1 ? "" : "s");
550 return os;
551 }
552
553
554 ////////////////////////////////////////
555
556
557 template<typename T>
558 inline
559 Vector<T>::Vector(const Vector& other): mData(new T[other.mSize]), mSize(other.mSize)
560 {
561 tbb::parallel_for(SizeRange(0, mSize),
562 internal::CopyOp<T>(/*from=*/other.mData, /*to=*/mData));
563 }
564
565
566 template<typename T>
567 inline
568 9 Vector<T>& Vector<T>::operator=(const Vector<T>& other)
569 {
570 // Update the internal storage to the correct size
571
572
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
 9 if (mSize != other.mSize) {
573 ✗ mSize = other.mSize;
574 ✗ delete[] mData;
575 ✗ mData = new T[mSize];
576 }
577
578 // Deep copy the data
579 18 tbb::parallel_for(SizeRange(0, mSize),
580 9 internal::CopyOp<T>(/*from=*/other.mData, /*to=*/mData));
581
582 9 return *this;
583 }
584
585
586 template<typename T>
587 inline void
588 Vector<T>::resize(SizeType n)
589 {
590 if (n != mSize) {
591 if (mData) delete[] mData;
592 mData = new T[n];
593 mSize = n;
594 }
595 }
596
597
598 template<typename T>
599 inline void
600 Vector<T>::fill(const ValueType& value)
601 {
602 951 tbb::parallel_for(SizeRange(0, mSize), internal::FillOp<T>(mData, value));
603 }
604
605
606 template<typename T>
607 template<typename Scalar>
608 struct Vector<T>::ScaleOp
609 {
610 7 ScaleOp(T* data_, const Scalar& s_): data(data_), s(s_) {}
611
612 void operator()(const SizeRange& range) const {
613
4/4
✓ Branch 0 taken 223274 times.
✓ Branch 1 taken 57 times.
✓ Branch 2 taken 3354224 times.
✓ Branch 3 taken 848 times.
 3578403 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) data[n] *= s;
614 }
615
616 T* data;
617 const Scalar s;
618 };
619
620
621 template<typename T>
622 template<typename Scalar>
623 inline void
624 Vector<T>::scale(const Scalar& s)
625 {
626
8/32
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 2 times.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✓ Branch 43 taken 1 times.
✗ Branch 44 not taken.
✓ Branch 46 taken 1 times.
✗ Branch 47 not taken.
 7 tbb::parallel_for(SizeRange(0, mSize), ScaleOp<Scalar>(mData, s));
627 }
628
629
630 template<typename T>
631 struct Vector<T>::DeterministicDotProductOp
632 {
633 1258 DeterministicDotProductOp(const T* a_, const T* b_,
634 const SizeType binCount_, const SizeType arraySize_, T* reducetmp_):
635 1258 a(a_), b(b_), binCount(binCount_), arraySize(arraySize_), reducetmp(reducetmp_) {}
636
637 125800 void operator()(const SizeRange& range) const
638 {
639 125800 const SizeType binSize = arraySize / binCount;
640
641 // Iterate over bins (array segments)
642
2/2
✓ Branch 0 taken 125800 times.
✓ Branch 1 taken 125800 times.
 251600 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
643 125800 const SizeType begin = n * binSize;
644
2/2
✓ Branch 0 taken 1258 times.
✓ Branch 1 taken 124542 times.
 125800 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
2/2
✓ Branch 0 taken 1071899673 times.
✓ Branch 1 taken 125800 times.
 1072025473 for (SizeType i = begin; i < end; ++i) { sum += a[i] * b[i]; }
649 // Store the partial sum
650 125800 reducetmp[n] = sum;
651 }
652 125800 }
653
654
655 const T* a;
656 const T* b;
657 const SizeType binCount;
658 const SizeType arraySize;
659 T* reducetmp;
660 };
661
662 template<typename T>
663 inline T
664 1415 Vector<T>::dot(const Vector<T>& other) const
665 {
666
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1415 times.
 1415 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
2/2
✓ Branch 0 taken 157 times.
✓ Branch 1 taken 1258 times.
 1415 if (arraySize < 1024) {
676
677 // Compute the dot product in serial for small arrays
678
679
2/2
✓ Branch 0 taken 139090 times.
✓ Branch 1 taken 157 times.
 139247 for (SizeType n = 0; n < arraySize; ++n) {
680 139090 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 1258 tbb::parallel_for(SizeRange(0, binCount),
693 DeterministicDotProductOp(aData, bData, binCount, arraySize, partialSums));
694
695
2/2
✓ Branch 0 taken 125800 times.
✓ Branch 1 taken 1258 times.
 127058 for (SizeType n = 0; n < binCount; ++n) {
696 125800 result += partialSums[n];
697 }
698 }
699
700 1415 return result;
701 }
702
703
704 template<typename T>
705 struct Vector<T>::InfNormOp
706 {
707 489 InfNormOp(const T* data_): data(data_) {}
708
709 T operator()(const SizeRange& range, T maxValue) const
710 {
711
2/2
✓ Branch 0 taken 361156103 times.
✓ Branch 1 taken 92418 times.
 361248521 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
712
2/2
✓ Branch 0 taken 100691 times.
✓ Branch 1 taken 361055412 times.
 361256794 maxValue = Max(maxValue, Abs(data[n]));
713 }
714 92418 return maxValue;
715 }
716
717 const T* data;
718 };
719
720
721 template<typename T>
722 inline T
723 489 Vector<T>::infNorm() const
724 {
725 // Parallelize over the elements of this vector.
726 489 T result = tbb::parallel_reduce(SizeRange(0, this->size()), /*seed=*/zeroVal<T>(),
727 2845 InfNormOp(this->data()), /*join=*/[](T max1, T max2) { return Max(max1, max2); });
728 489 return result;
729 }
730
731
732 template<typename T>
733 struct Vector<T>::IsFiniteOp
734 {
735 9 IsFiniteOp(const T* data_): data(data_) {}
736
737 bool operator()(const SizeRange& range, bool finite) const
738 {
739
2/4
✓ Branch 0 taken 71 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 848 times.
✗ Branch 3 not taken.
 919 if (finite) {
740
4/4
✓ Branch 0 taken 223328 times.
✓ Branch 1 taken 71 times.
✓ Branch 2 taken 3354180 times.
✓ Branch 3 taken 848 times.
 3578427 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
741
2/4
✓ Branch 0 taken 223328 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 3354180 times.
✗ Branch 3 not taken.
 3577508 if (!std::isfinite(data[n])) return false;
742 }
743 }
744 return finite;
745 }
746
747 const T* data;
748 };
749
750
751 template<typename T>
752 inline bool
753 9 Vector<T>::isFinite() const
754 {
755 // Parallelize over the elements of this vector.
756 9 bool finite = tbb::parallel_reduce(SizeRange(0, this->size()), /*seed=*/true,
757 IsFiniteOp(this->data()),
758 9 /*join=*/[](bool finite1, bool finite2) { return (finite1 && finite2); });
759 9 return finite;
760 }
761
762
763 template<typename T>
764 template<typename OtherValueType>
765 struct Vector<T>::EqOp
766 {
767 2 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
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 10 if (equal) {
772
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 20 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
773
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 10 if (!isApproxEqual(a[n], b[n], eps)) return false;
774 }
775 }
776 return equal;
777 }
778
779 const T* a;
780 const OtherValueType* b;
781 const T eps;
782 };
783
784
785 template<typename T>
786 template<typename OtherValueType>
787 inline bool
788 2 Vector<T>::eq(const Vector<OtherValueType>& other, ValueType eps) const
789 {
790
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if (this->size() != other.size()) return false;
791 2 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 2 return equal;
795 }
796
797
798 template<typename T>
799 inline std::string
800 Vector<T>::str() const
801 {
802 std::ostringstream ostr;
803 ostr << "[";
804 std::string sep;
805 for (SizeType n = 0, N = this->size(); n < N; ++n) {
806 ostr << sep << (*this)[n];
807 sep = ", ";
808 }
809 ostr << "]";
810 return ostr.str();
811 }
812
813
814 ////////////////////////////////////////
815
816
817 template<typename ValueType, SizeType STENCIL_SIZE>
818 const ValueType SparseStencilMatrix<ValueType, STENCIL_SIZE>::sZeroValue = zeroVal<ValueType>();
819
820
821 template<typename ValueType, SizeType STENCIL_SIZE>
822 inline
823 58 SparseStencilMatrix<ValueType, STENCIL_SIZE>::SparseStencilMatrix(SizeType numRows)
824 : mNumRows(numRows)
825 58 , mValueArray(new ValueType[mNumRows * STENCIL_SIZE])
826
1/2
✓ Branch 1 taken 29 times.
✗ Branch 2 not taken.
 58 , mColumnIdxArray(new SizeType[mNumRows * STENCIL_SIZE])
827
1/2
✓ Branch 2 taken 29 times.
✗ Branch 3 not taken.
 116 , mRowSizeArray(new SizeType[mNumRows])
828 {
829 // Initialize the matrix to a null state by setting the size of each row to zero.
830
1/4
✓ Branch 1 taken 29 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 58 tbb::parallel_for(SizeRange(0, mNumRows),
831 internal::FillOp<SizeType>(mRowSizeArray.get(), /*value=*/0));
832 58 }
833
834
835 template<typename ValueType, SizeType STENCIL_SIZE>
836 struct SparseStencilMatrix<ValueType, STENCIL_SIZE>::MatrixCopyOp
837 {
838 2 MatrixCopyOp(const SparseStencilMatrix& from_, SparseStencilMatrix& to_):
839 2 from(&from_), to(&to_) {}
840
841 40 void operator()(const SizeRange& range) const
842 {
843 40 const ValueType* fromVal = from->mValueArray.get();
844 const SizeType* fromCol = from->mColumnIdxArray.get();
845 40 ValueType* toVal = to->mValueArray.get();
846 SizeType* toCol = to->mColumnIdxArray.get();
847
2/2
✓ Branch 0 taken 40 times.
✓ Branch 1 taken 40 times.
 80 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
848 40 toVal[n] = fromVal[n];
849 40 toCol[n] = fromCol[n];
850 }
851 40 }
852
853 const SparseStencilMatrix* from; SparseStencilMatrix* to;
854 };
855
856
857 template<typename ValueType, SizeType STENCIL_SIZE>
858 inline
859 2 SparseStencilMatrix<ValueType, STENCIL_SIZE>::SparseStencilMatrix(const SparseStencilMatrix& other)
860 2 : mNumRows(other.mNumRows)
861 2 , mValueArray(new ValueType[mNumRows * STENCIL_SIZE])
862
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 , mColumnIdxArray(new SizeType[mNumRows * STENCIL_SIZE])
863
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 4 , mRowSizeArray(new SizeType[mNumRows])
864 {
865
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 SizeType size = mNumRows * STENCIL_SIZE;
866
867 // Copy the value and column index arrays from the other matrix to this matrix.
868
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 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
1/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 2 tbb::parallel_for(SizeRange(0, mNumRows),
872 internal::CopyOp<SizeType>(/*from=*/other.mRowSizeArray.get(), /*to=*/mRowSizeArray.get()));
873 2 }
874
875
876 template<typename ValueType, SizeType STENCIL_SIZE>
877 inline void
878 88 SparseStencilMatrix<ValueType, STENCIL_SIZE>::setValue(SizeType row, SizeType col,
879 const ValueType& val)
880 {
881
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 44 times.
 88 assert(row < mNumRows);
882 88 this->getRowEditor(row).setValue(col, val);
883 88 }
884
885
886 template<typename ValueType, SizeType STENCIL_SIZE>
887 inline const ValueType&
888 28339968 SparseStencilMatrix<ValueType, STENCIL_SIZE>::getValue(SizeType row, SizeType col) const
889 {
890
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 14169984 times.
 28339968 assert(row < mNumRows);
891 28339968 return this->getConstRow(row).getValue(col);
892 }
893
894
895 template<typename ValueType, SizeType STENCIL_SIZE>
896 inline const ValueType&
897 SparseStencilMatrix<ValueType, STENCIL_SIZE>::operator()(SizeType row, SizeType col) const
898 {
899 return this->getValue(row,col);
900 }
901
902
903 template<typename ValueType, SizeType STENCIL_SIZE>
904 template<typename Scalar>
905 struct SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowScaleOp
906 {
907 9 RowScaleOp(SparseStencilMatrix& m, const Scalar& s_): mat(&m), s(s_) {}
908
909 1428 void operator()(const SizeRange& range) const
910 {
911
2/2
✓ Branch 0 taken 3584700 times.
✓ Branch 1 taken 1428 times.
 3586128 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
912 3584700 RowEditor row = mat->getRowEditor(n);
913 row.scale(s);
914 }
915 1428 }
916
917 SparseStencilMatrix* mat;
918 const Scalar s;
919 };
920
921
922 template<typename ValueType, SizeType STENCIL_SIZE>
923 template<typename Scalar>
924 inline void
925 SparseStencilMatrix<ValueType, STENCIL_SIZE>::scale(const Scalar& s)
926 {
927 // Parallelize over the rows in the matrix.
928
10/32
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 2 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 2 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 2 times.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✓ Branch 43 taken 1 times.
✗ Branch 44 not taken.
✓ Branch 46 taken 1 times.
✗ Branch 47 not taken.
 9 tbb::parallel_for(SizeRange(0, mNumRows), RowScaleOp<Scalar>(*this, s));
929 }
930
931
932 template<typename ValueType, SizeType STENCIL_SIZE>
933 template<typename VecValueType>
934 struct SparseStencilMatrix<ValueType, STENCIL_SIZE>::VecMultOp
935 {
936 482 VecMultOp(const SparseStencilMatrix& m, const VecValueType* i, VecValueType* o):
937 482 mat(&m), in(i), out(o) {}
938
939 70612 void operator()(const SizeRange& range) const
940 {
941
2/2
✓ Branch 0 taken 357585797 times.
✓ Branch 1 taken 70356 times.
 357663611 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
942 357592999 ConstRow row = mat->getConstRow(n);
943 357592999 out[n] = row.dot(in, mat->numRows());
944 }
945 70612 }
946
947 const SparseStencilMatrix* mat;
948 const VecValueType* in;
949 VecValueType* out;
950 };
951
952
953 template<typename ValueType, SizeType STENCIL_SIZE>
954 template<typename VecValueType>
955 inline void
956 475 SparseStencilMatrix<ValueType, STENCIL_SIZE>::vectorMultiply(
957 const Vector<VecValueType>& inVec, Vector<VecValueType>& resultVec) const
958 {
959
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 473 times.
 475 if (inVec.size() != mNumRows) {
960 ✗ OPENVDB_THROW(ArithmeticError, "matrix and input vector have incompatible sizes ("
961 << mNumRows << "x" << mNumRows << " vs. " << inVec.size() << ")");
962 }
963
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 473 times.
 475 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 475 }
970
971
972 template<typename ValueType, SizeType STENCIL_SIZE>
973 template<typename VecValueType>
974 inline void
975 SparseStencilMatrix<ValueType, STENCIL_SIZE>::vectorMultiply(
976 const VecValueType* inVec, VecValueType* resultVec) const
977 {
978 // Parallelize over the rows in the matrix.
979 482 tbb::parallel_for(SizeRange(0, mNumRows),
980 VecMultOp<VecValueType>(*this, inVec, resultVec));
981 }
982
983
984 template<typename ValueType, SizeType STENCIL_SIZE>
985 template<typename OtherValueType>
986 struct SparseStencilMatrix<ValueType, STENCIL_SIZE>::EqOp
987 {
988 2 EqOp(const SparseStencilMatrix& a_,
989 const SparseStencilMatrix<OtherValueType, STENCIL_SIZE>& b_, ValueType e):
990 2 a(&a_), b(&b_), eps(e) {}
991
992 10 bool operator()(const SizeRange& range, bool equal) const
993 {
994
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 10 if (equal) {
995
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 20 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
996
1/2
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
 10 if (!a->getConstRow(n).eq(b->getConstRow(n), eps)) return false;
997 }
998 }
999 return equal;
1000 }
1001
1002 const SparseStencilMatrix* a;
1003 const SparseStencilMatrix<OtherValueType, STENCIL_SIZE>* b;
1004 const ValueType eps;
1005 };
1006
1007
1008 template<typename ValueType, SizeType STENCIL_SIZE>
1009 template<typename OtherValueType>
1010 inline bool
1011 2 SparseStencilMatrix<ValueType, STENCIL_SIZE>::eq(
1012 const SparseStencilMatrix<OtherValueType, STENCIL_SIZE>& other, ValueType eps) const
1013 {
1014
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if (this->numRows() != other.numRows()) return false;
1015 2 bool equal = tbb::parallel_reduce(SizeRange(0, this->numRows()), /*seed=*/true,
1016 EqOp<OtherValueType>(*this, other, eps),
1017 1 /*join=*/[](bool eq1, bool eq2) { return (eq1 && eq2); });
1018 2 return equal;
1019 }
1020
1021
1022 template<typename ValueType, SizeType STENCIL_SIZE>
1023 struct SparseStencilMatrix<ValueType, STENCIL_SIZE>::IsFiniteOp
1024 {
1025 2 IsFiniteOp(const SparseStencilMatrix& m): mat(&m) {}
1026
1027 10 bool operator()(const SizeRange& range, bool finite) const
1028 {
1029
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 10 if (finite) {
1030
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 20 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1031 10 const ConstRow row = mat->getConstRow(n);
1032
2/2
✓ Branch 0 taken 26 times.
✓ Branch 1 taken 10 times.
 36 for (ConstValueIter it = row.cbegin(); it; ++it) {
1033
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 26 times.
 26 if (!std::isfinite(*it)) return false;
1034 }
1035 }
1036 }
1037 return finite;
1038 }
1039
1040 const SparseStencilMatrix* mat;
1041 };
1042
1043
1044 template<typename ValueType, SizeType STENCIL_SIZE>
1045 inline bool
1046 2 SparseStencilMatrix<ValueType, STENCIL_SIZE>::isFinite() const
1047 {
1048 // Parallelize over the rows of this matrix.
1049 2 bool finite = tbb::parallel_reduce(SizeRange(0, this->numRows()), /*seed=*/true,
1050 ✗ IsFiniteOp(*this), /*join=*/[](bool finite1, bool finite2) { return (finite1&&finite2); });
1051 2 return finite;
1052 }
1053
1054
1055 template<typename ValueType, SizeType STENCIL_SIZE>
1056 inline std::string
1057 SparseStencilMatrix<ValueType, STENCIL_SIZE>::str() const
1058 {
1059 std::ostringstream ostr;
1060 for (SizeType n = 0, N = this->size(); n < N; ++n) {
1061 ostr << n << ": " << this->getConstRow(n).str() << "\n";
1062 }
1063 return ostr.str();
1064 }
1065
1066
1067 template<typename ValueType, SizeType STENCIL_SIZE>
1068 inline typename SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowEditor
1069 141004966 SparseStencilMatrix<ValueType, STENCIL_SIZE>::getRowEditor(SizeType i)
1070 {
1071
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 70502483 times.
 141004966 assert(i < mNumRows);
1072 14338852 const SizeType head = i * STENCIL_SIZE;
1073 141004966 return RowEditor(&mValueArray[head], &mColumnIdxArray[head], mRowSizeArray[i], mNumRows);
1074 }
1075
1076
1077 template<typename ValueType, SizeType STENCIL_SIZE>
1078 inline typename SparseStencilMatrix<ValueType, STENCIL_SIZE>::ConstRow
1079 2314521946 SparseStencilMatrix<ValueType, STENCIL_SIZE>::getConstRow(SizeType i) const
1080 {
1081
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1157260973 times.
 2314521946 assert(i < mNumRows);
1082 757821594 const SizeType head = i * STENCIL_SIZE; // index for this row into main storage
1083 2314521946 return ConstRow(&mValueArray[head], &mColumnIdxArray[head], mRowSizeArray[i]);
1084 }
1085
1086
1087 template<typename ValueType, SizeType STENCIL_SIZE>
1088 template<typename DataType>
1089 inline SizeType
1090 1859373292 SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::find(SizeType columnIdx) const
1091 {
1092
2/2
✓ Branch 0 taken 918946920 times.
✓ Branch 1 taken 10739726 times.
 1859373292 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 1837893840 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 1837893840 return static_cast<SizeType>(colPtr - mData.mCols);
1099 }
1100
1101
1102 template<typename ValueType, SizeType STENCIL_SIZE>
1103 template<typename DataType>
1104 inline const ValueType&
1105 SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::getValue(
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 }
1114 return SparseStencilMatrix::sZeroValue;
1115 }
1116
1117 template<typename ValueType, SizeType STENCIL_SIZE>
1118 template<typename DataType>
1119 inline const ValueType&
1120 1640716430 SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::getValue(SizeType columnIdx) const
1121 {
1122 1640716430 SizeType idx = this->find(columnIdx);
1123
3/4
✓ Branch 0 taken 820358215 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 788942637 times.
✓ Branch 3 taken 31415578 times.
 1640716430 if (idx < this->size() && this->column(idx) == columnIdx) {
1124 1577885274 return this->value(idx);
1125 }
1126 return SparseStencilMatrix::sZeroValue;
1127 }
1128
1129
1130 template<typename ValueType, SizeType STENCIL_SIZE>
1131 template<typename DataType>
1132 inline typename SparseStencilMatrix<ValueType, STENCIL_SIZE>::ConstValueIter
1133 SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::cbegin() const
1134 {
1135 21325005 return ConstValueIter(mData);
1136 }
1137
1138
1139 template<typename ValueType, SizeType STENCIL_SIZE>
1140 template<typename DataType>
1141 template<typename OtherDataType>
1142 inline bool
1143 10 SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::eq(
1144 const RowBase<OtherDataType>& other, ValueType eps) const
1145 {
1146
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 10 if (this->size() != other.size()) return false;
1147
3/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10 times.
 28 for (ConstValueIter it = cbegin(), oit = other.cbegin(); it || oit; ++it, ++oit) {
1148
1/2
✓ Branch 0 taken 18 times.
✗ Branch 1 not taken.
 18 if (it.column() != oit.column()) return false;
1149
1/2
✓ Branch 0 taken 18 times.
✗ Branch 1 not taken.
 18 if (!isApproxEqual(*it, *oit, eps)) return false;
1150 }
1151 10 return true;
1152 }
1153
1154
1155 template<typename ValueType, SizeType STENCIL_SIZE>
1156 template<typename DataType>
1157 template<typename VecValueType>
1158 inline VecValueType
1159
2/2
✓ Branch 0 taken 1065587933 times.
✓ Branch 1 taken 8 times.
 2131175882 SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::dot(
1160 const VecValueType* inVec, SizeType vecSize) const
1161 {
1162 VecValueType result = zeroVal<VecValueType>();
1163
2/2
✓ Branch 0 taken 5286932057 times.
✓ Branch 1 taken 1065587941 times.
 12705039996 for (SizeType idx = 0, N = std::min(vecSize, this->size()); idx < N; ++idx) {
1164 10573864114 result += static_cast<VecValueType>(this->value(idx) * inVec[this->column(idx)]);
1165 }
1166 2131175882 return result;
1167 }
1168
1169 template<typename ValueType, SizeType STENCIL_SIZE>
1170 template<typename DataType>
1171 template<typename VecValueType>
1172 inline VecValueType
1173 SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::dot(
1174 const Vector<VecValueType>& inVec) const
1175 {
1176 708002144 return dot(inVec.data(), inVec.size());
1177 }
1178
1179
1180 template<typename ValueType, SizeType STENCIL_SIZE>
1181 template<typename DataType>
1182 inline std::string
1183 SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::str() const
1184 {
1185 std::ostringstream ostr;
1186 std::string sep;
1187 for (SizeType n = 0, N = this->size(); n < N; ++n) {
1188 ostr << sep << "(" << this->column(n) << ", " << this->value(n) << ")";
1189 sep = ", ";
1190 }
1191 return ostr.str();
1192 }
1193
1194
1195 template<typename ValueType, SizeType STENCIL_SIZE>
1196 inline
1197 SparseStencilMatrix<ValueType, STENCIL_SIZE>::ConstRow::ConstRow(
1198 const ValueType* valueHead, const SizeType* columnHead, const SizeType& rowSize):
1199 ConstRowBase(ConstRowData(const_cast<ValueType*>(valueHead),
1200 1157260973 const_cast<SizeType*>(columnHead), const_cast<SizeType&>(rowSize)))
1201 {
1202 }
1203
1204
1205 template<typename ValueType, SizeType STENCIL_SIZE>
1206 inline
1207 70502483 SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowEditor::RowEditor(
1208 ValueType* valueHead, SizeType* columnHead, SizeType& rowSize, SizeType colSize):
1209 70502483 RowBase<>(RowData(valueHead, columnHead, rowSize)), mNumColumns(colSize)
1210 {
1211 }
1212
1213
1214 template<typename ValueType, SizeType STENCIL_SIZE>
1215 inline void
1216 SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowEditor::clear()
1217 {
1218 // Note: since mSize is a reference, this modifies the underlying matrix.
1219 7155006 RowBase<>::mData.mSize = 0;
1220 }
1221
1222
1223 template<typename ValueType, SizeType STENCIL_SIZE>
1224 inline SizeType
1225 218656862 SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowEditor::setValue(
1226 SizeType column, const ValueType& value)
1227 {
1228
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 109328431 times.
 218656862 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 218656862 SizeType offset = this->find(column);
1235
1236
4/4
✓ Branch 0 taken 59318229 times.
✓ Branch 1 taken 50010202 times.
✓ Branch 2 taken 56178033 times.
✓ Branch 3 taken 3140196 times.
 218656862 if (offset < data.mSize && data.mCols[offset] == column) {
1237 // If the column already exists, just update its value.
1238 112356066 data.mVals[offset] = value;
1239 112356066 return data.mSize;
1240 }
1241
1242 // Check that it is safe to add a new column.
1243
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 53150398 times.
 106300796 assert(data.mSize < this->capacity());
1244
1245
2/2
✓ Branch 0 taken 50010202 times.
✓ Branch 1 taken 3140196 times.
 106300796 if (offset >= data.mSize) {
1246 // The new column's index is larger than any existing index. Append the new column.
1247 100020404 data.mVals[data.mSize] = value;
1248 100020404 data.mCols[data.mSize] = column;
1249 } else {
1250 // Insert the new column at the computed offset after shifting subsequent columns.
1251
2/2
✓ Branch 0 taken 4811545 times.
✓ Branch 1 taken 3140196 times.
 15903482 for (SizeType i = data.mSize; i > offset; --i) {
1252 9623090 data.mVals[i] = data.mVals[i - 1];
1253 9623090 data.mCols[i] = data.mCols[i - 1];
1254 }
1255 6280392 data.mVals[offset] = value;
1256 6280392 data.mCols[offset] = column;
1257 }
1258 106300796 ++data.mSize;
1259
1260 106300796 return data.mSize;
1261 }
1262
1263
1264 template<typename ValueType, SizeType STENCIL_SIZE>
1265 template<typename Scalar>
1266 inline void
1267 SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowEditor::scale(const Scalar& s)
1268 {
1269
2/2
✓ Branch 0 taken 24810396 times.
✓ Branch 1 taken 3584700 times.
 28395096 for (int idx = 0, N = this->size(); idx < N; ++idx) {
1270 24810396 RowBase<>::mData.mVals[idx] *= s;
1271 }
1272 }
1273
1274
1275 ////////////////////////////////////////
1276
1277
1278 /// Diagonal preconditioner
1279 template<typename MatrixType>
1280 class JacobiPreconditioner: public Preconditioner<typename MatrixType::ValueType>
1281 {
1282 private:
1283 struct InitOp;
1284 struct ApplyOp;
1285
1286 public:
1287 using ValueType = typename MatrixType::ValueType;
1288 using BaseType = Preconditioner<ValueType>;
1289 using VectorType = Vector<ValueType>;
1290 using Ptr = SharedPtr<JacobiPreconditioner>;
1291
1292 1 JacobiPreconditioner(const MatrixType& A): BaseType(A), mDiag(A.numRows())
1293 {
1294 // Initialize vector mDiag with the values from the matrix diagonal.
1295
1/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 1 tbb::parallel_for(SizeRange(0, A.numRows()), InitOp(A, mDiag.data()));
1296 1 }
1297
1298
2/6
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 2 ~JacobiPreconditioner() override = default;
1299
1300 3 void apply(const Vector<ValueType>& r, Vector<ValueType>& z) override
1301 {
1302 const SizeType size = mDiag.size();
1303
1304
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 3 assert(r.size() == z.size());
1305
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 3 assert(r.size() == size);
1306
1307 3 tbb::parallel_for(SizeRange(0, size), ApplyOp(mDiag.data(), r.data(), z.data()));
1308 3 }
1309
1310 /// Return @c true if all values along the diagonal are finite.
1311 bool isFinite() const { return mDiag.isFinite(); }
1312
1313 private:
1314 // Functor for use with tbb::parallel_for()
1315 struct InitOp
1316 {
1317
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 InitOp(const MatrixType& m, ValueType* v): mat(&m), vec(v) {}
1318 5 void operator()(const SizeRange& range) const {
1319
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 10 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1320
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
 5 const ValueType val = mat->getValue(n, n);
1321 ✗ assert(!isApproxZero(val, ValueType(0.0001)));
1322 5 vec[n] = static_cast<ValueType>(1.0 / val);
1323 }
1324 5 }
1325 const MatrixType* mat; ValueType* vec;
1326 };
1327
1328 // Functor for use with tbb::parallel_reduce()
1329 struct ApplyOp
1330 {
1331 3 ApplyOp(const ValueType* x_, const ValueType* y_, ValueType* out_):
1332 3 x(x_), y(y_), out(out_) {}
1333 void operator()(const SizeRange& range) const {
1334
2/4
✓ Branch 0 taken 15 times.
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 30 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) out[n] = x[n] * y[n];
1335 }
1336 const ValueType *x, *y; ValueType* out;
1337 };
1338
1339 // The Jacobi preconditioner is a diagonal matrix
1340 VectorType mDiag;
1341 }; // class JacobiPreconditioner
1342
1343
1344 /// Preconditioner using incomplete Cholesky factorization
1345 template<typename MatrixType>
1346 class IncompleteCholeskyPreconditioner: public Preconditioner<typename MatrixType::ValueType>
1347 {
1348 private:
1349 struct CopyToLowerOp;
1350 struct TransposeOp;
1351
1352 public:
1353 using ValueType = typename MatrixType::ValueType;
1354 using BaseType = Preconditioner<ValueType>;
1355 using VectorType = Vector<ValueType>;
1356 using Ptr = SharedPtr<IncompleteCholeskyPreconditioner>;
1357 using TriangularMatrix = SparseStencilMatrix<ValueType, 4>;
1358 using TriangleConstRow = typename TriangularMatrix::ConstRow;
1359 using TriangleRowEditor = typename TriangularMatrix::RowEditor;
1360
1361 8 IncompleteCholeskyPreconditioner(const MatrixType& matrix)
1362 : BaseType(matrix)
1363 , mLowerTriangular(matrix.numRows())
1364 , mUpperTriangular(matrix.numRows())
1365
1/2
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
 8 , 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
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 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 8 mPassedCompatibilityCondition = true;
1392
1393
2/2
✓ Branch 0 taken 3577503 times.
✓ Branch 1 taken 8 times.
 3577511 for (SizeType k = 0; k < numRows; ++k) {
1394
1395 3577503 TriangleConstRow crow_k = mLowerTriangular.getConstRow(k);
1396 3577503 ValueType diagonalValue = crow_k.getValue(k);
1397
1398 // Test if the matrix build has failed.
1399
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3577503 times.
 3577503 if (diagonalValue < 1.e-5) {
1400 ✗ mPassedCompatibilityCondition = false;
1401 ✗ break;
1402 }
1403
1404 3577503 diagonalValue = Sqrt(diagonalValue);
1405
1406 3577503 TriangleRowEditor row_k = mLowerTriangular.getRowEditor(k);
1407 3577503 row_k.setValue(k, diagonalValue);
1408
1409 // Exploit the fact that the matrix is symmetric.
1410 3577503 typename MatrixType::ConstRow srcRow = matrix.getConstRow(k);
1411 typename MatrixType::ConstValueIter citer = srcRow.cbegin();
1412
2/2
✓ Branch 0 taken 24762455 times.
✓ Branch 1 taken 3577503 times.
 28339958 for ( ; citer; ++citer) {
1413 SizeType ii = citer.column();
1414
2/2
✓ Branch 0 taken 14169979 times.
✓ Branch 1 taken 10592476 times.
 24762455 if (ii < k+1) continue; // look above diagonal
1415
1416 10592476 TriangleRowEditor row_ii = mLowerTriangular.getRowEditor(ii);
1417
1418 10592476 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
2/2
✓ Branch 0 taken 24762455 times.
✓ Branch 1 taken 3577503 times.
 28339958 for ( ; citer; ++citer) {
1424 SizeType j = citer.column();
1425
2/2
✓ Branch 0 taken 14169979 times.
✓ Branch 1 taken 10592476 times.
 24762455 if (j < k+1) continue;
1426
1427 10592476 TriangleConstRow row_j = mLowerTriangular.getConstRow(j);
1428 10592476 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 10592476 typename MatrixType::ConstRow mask = matrix.getConstRow(j);
1433 typename MatrixType::ConstValueIter maskIter = mask.cbegin();
1434
2/2
✓ Branch 0 taken 73563632 times.
✓ Branch 1 taken 10592476 times.
 84156108 for ( ; maskIter; ++maskIter) {
1435 SizeType i = maskIter.column();
1436
2/2
✓ Branch 0 taken 31555578 times.
✓ Branch 1 taken 42008054 times.
 73563632 if (i < j) continue;
1437
1438 42008054 TriangleConstRow crow_i = mLowerTriangular.getConstRow(i);
1439 42008054 ValueType a_ij = crow_i.getValue(j);
1440 42008054 ValueType a_ik = crow_i.getValue(k);
1441 42008054 TriangleRowEditor row_i = mLowerTriangular.getRowEditor(i);
1442 42008054 a_ij -= a_ik * a_jk;
1443
1444 42008054 row_i.setValue(j, a_ij);
1445 }
1446 }
1447 }
1448
1449 // Build the transpose of the IC matrix: mUpperTriangular
1450
1/4
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 8 tbb::parallel_for(SizeRange(0, numRows),
1451 TransposeOp(matrix, mLowerTriangular, mUpperTriangular));
1452 8 }
1453
1454
1/2
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
 46 ~IncompleteCholeskyPreconditioner() override = default;
1455
1456 7 bool isValid() const override { return mPassedCompatibilityCondition; }
1457
1458 468 void apply(const Vector<ValueType>& rVec, Vector<ValueType>& zVec) override
1459 {
1460
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 468 times.
 468 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
1/2
✓ Branch 0 taken 468 times.
✗ Branch 1 not taken.
 468 if (size == 0) return;
1472
1473
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 468 times.
 468 assert(rVec.size() == size);
1474
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 468 times.
 468 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
2/2
✓ Branch 0 taken 354001072 times.
✓ Branch 1 taken 468 times.
 354001540 for (SizeType i = 0; i < size; ++i) {
1483 354001072 typename TriangularMatrix::ConstRow row = mLowerTriangular.getConstRow(i);
1484 354001072 ValueType diagonal = row.getValue(i);
1485 ValueType dot = row.dot(mTempVec);
1486 354001072 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
2/2
✓ Branch 0 taken 354001072 times.
✓ Branch 1 taken 468 times.
 354001540 for (SizeType ii = 0; ii < size; ++ii) {
1495 354001072 SizeType i = size - 1 - ii;
1496 354001072 typename TriangularMatrix::ConstRow row = mUpperTriangular.getConstRow(i);
1497 354001072 ValueType diagonal = row.getValue(i);
1498 ValueType dot = row.dot(zVec);
1499 354001072 zData[i] = (tmpData[i] - dot) / diagonal;
1500 if (!std::isfinite(zData[i])) {
1501 OPENVDB_LOG_DEBUG_RUNTIME("2 diagonal was " << diagonal);
1502 }
1503 }
1504 }
1505
1506 const TriangularMatrix& lowerMatrix() const { return mLowerTriangular; }
1507 const TriangularMatrix& upperMatrix() const { return mUpperTriangular; }
1508
1509 private:
1510 // Functor for use with tbb::parallel_for()
1511 struct CopyToLowerOp
1512 {
1513
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 CopyToLowerOp(const MatrixType& m, TriangularMatrix& l): mat(&m), lower(&l) {}
1514 946 void operator()(const SizeRange& range) const {
1515
2/2
✓ Branch 0 taken 3577503 times.
✓ Branch 1 taken 946 times.
 3578449 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1516 3577503 typename TriangularMatrix::RowEditor outRow = lower->getRowEditor(n);
1517 outRow.clear();
1518 3577503 typename MatrixType::ConstRow inRow = mat->getConstRow(n);
1519
2/2
✓ Branch 0 taken 24762455 times.
✓ Branch 1 taken 3577503 times.
 28339958 for (typename MatrixType::ConstValueIter it = inRow.cbegin(); it; ++it) {
1520
2/2
✓ Branch 0 taken 10592476 times.
✓ Branch 1 taken 14169979 times.
 24762455 if (it.column() > n) continue; // skip above diagonal
1521 14169979 outRow.setValue(it.column(), *it);
1522 }
1523 }
1524 946 }
1525 const MatrixType* mat; TriangularMatrix* lower;
1526 };
1527
1528 // Functor for use with tbb::parallel_for()
1529 struct TransposeOp
1530 {
1531 8 TransposeOp(const MatrixType& m, const TriangularMatrix& l, TriangularMatrix& u):
1532
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 mat(&m), lower(&l), upper(&u) {}
1533 971 void operator()(const SizeRange& range) const {
1534
2/2
✓ Branch 0 taken 3577503 times.
✓ Branch 1 taken 971 times.
 3578474 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1535 3577503 typename TriangularMatrix::RowEditor outRow = upper->getRowEditor(n);
1536 outRow.clear();
1537 // Use the fact that matrix is symmetric.
1538 3577503 typename MatrixType::ConstRow inRow = mat->getConstRow(n);
1539
2/2
✓ Branch 0 taken 24762455 times.
✓ Branch 1 taken 3577503 times.
 28339958 for (typename MatrixType::ConstValueIter it = inRow.cbegin(); it; ++it) {
1540 const SizeType column = it.column();
1541
2/2
✓ Branch 0 taken 10592476 times.
✓ Branch 1 taken 14169979 times.
 24762455 if (column < n) continue; // only set upper triangle
1542 14169979 outRow.setValue(column, lower->getValue(column, n));
1543 }
1544 }
1545 971 }
1546 const MatrixType* mat; const TriangularMatrix* lower; TriangularMatrix* upper;
1547 };
1548
1549 TriangularMatrix mLowerTriangular;
1550 TriangularMatrix mUpperTriangular;
1551 Vector<ValueType> mTempVec;
1552 bool mPassedCompatibilityCondition;
1553 }; // class IncompleteCholeskyPreconditioner
1554
1555
1556 ////////////////////////////////////////
1557
1558
1559 namespace internal {
1560
1561 /// Compute @e ax + @e y.
1562 template<typename T>
1563 inline void
1564 axpy(const T& a, const T* xVec, const T* yVec, T* resultVec, SizeType size)
1565 {
1566 1404 tbb::parallel_for(SizeRange(0, size), LinearOp<T>(a, xVec, yVec, resultVec));
1567 }
1568
1569 /// Compute @e ax + @e y.
1570 template<typename T>
1571 inline void
1572 1404 axpy(const T& a, const Vector<T>& xVec, const Vector<T>& yVec, Vector<T>& result)
1573 {
1574
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1404 times.
 1404 assert(xVec.size() == yVec.size());
1575
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1404 times.
 1404 assert(xVec.size() == result.size());
1576 axpy(a, xVec.data(), yVec.data(), result.data(), xVec.size());
1577 1404 }
1578
1579
1580 /// Compute @e r = @e b &minus; @e Ax.
1581 template<typename MatrixOperator, typename VecValueType>
1582 inline void
1583 9 computeResidual(const MatrixOperator& A, const VecValueType* x,
1584 const VecValueType* b, VecValueType* r)
1585 {
1586 // Compute r = A * x.
1587 A.vectorMultiply(x, r);
1588 // Compute r = b - r.
1589 9 tbb::parallel_for(SizeRange(0, A.numRows()), LinearOp<VecValueType>(-1.0, r, b, r));
1590 9 }
1591
1592 /// Compute @e r = @e b &minus; @e Ax.
1593 template<typename MatrixOperator, typename T>
1594 inline void
1595 9 computeResidual(const MatrixOperator& A, const Vector<T>& x, const Vector<T>& b, Vector<T>& r)
1596 {
1597
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
 9 assert(x.size() == b.size());
1598
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
 9 assert(x.size() == r.size());
1599
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
 9 assert(x.size() == A.numRows());
1600
1601 9 computeResidual(A, x.data(), b.data(), r.data());
1602 9 }
1603
1604 } // namespace internal
1605
1606
1607 ////////////////////////////////////////
1608
1609
1610 template<typename PositiveDefMatrix>
1611 inline State
1612 solve(
1613 const PositiveDefMatrix& Amat,
1614 const Vector<typename PositiveDefMatrix::ValueType>& bVec,
1615 Vector<typename PositiveDefMatrix::ValueType>& xVec,
1616 Preconditioner<typename PositiveDefMatrix::ValueType>& precond,
1617 const State& termination)
1618 {
1619 2 util::NullInterrupter interrupter;
1620
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 return solve(Amat, bVec, xVec, precond, interrupter, termination);
1621 }
1622
1623
1624 template<typename PositiveDefMatrix, typename Interrupter>
1625 inline State
1626 9 solve(
1627 const PositiveDefMatrix& Amat,
1628 const Vector<typename PositiveDefMatrix::ValueType>& bVec,
1629 Vector<typename PositiveDefMatrix::ValueType>& xVec,
1630 Preconditioner<typename PositiveDefMatrix::ValueType>& precond,
1631 Interrupter& interrupter,
1632 const State& termination)
1633 {
1634 using ValueType = typename PositiveDefMatrix::ValueType;
1635 using VectorType = Vector<ValueType>;
1636
1637 State result;
1638 9 result.success = false;
1639 9 result.iterations = 0;
1640 9 result.relativeError = 0.0;
1641 9 result.absoluteError = 0.0;
1642
1643 const SizeType size = Amat.numRows();
1644
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
 9 if (size == 0) {
1645 OPENVDB_LOG_WARN("pcg::solve(): matrix has dimension zero");
1646 return result;
1647 }
1648
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
 9 if (size != bVec.size()) {
1649 ✗ OPENVDB_THROW(ArithmeticError, "A and b have incompatible sizes"
1650 << size << "x" << size << " vs. " << bVec.size() << ")");
1651 }
1652
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
 9 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
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 const ValueType tmp = bVec.infNorm();
1664
1/2
✓ Branch 0 taken 9 times.
✗ Branch 1 not taken.
 9 const ValueType infNormOfB = isZero(tmp) ? ValueType(1) : tmp;
1665
1666 // Compute rVec: residual = b - Ax.
1667 VectorType rVec(size); // vector of residuals
1668
1669
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 internal::computeResidual(Amat, xVec, bVec, rVec);
1670
1671
2/4
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 9 times.
 9 assert(rVec.isFinite());
1672
1673 // Normalize the residual norm with the source norm and look for early out.
1674
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 result.absoluteError = static_cast<double>(rVec.infNorm());
1675 9 result.relativeError = static_cast<double>(result.absoluteError / infNormOfB);
1676
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
 9 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 9 ValueType minL2Error = std::numeric_limits<ValueType>::max();
1687 ValueType l2Error;
1688
1689 int iteration = 0;
1690
1/2
✓ Branch 0 taken 471 times.
✗ Branch 1 not taken.
 471 for ( ; iteration < termination.iterations; ++iteration) {
1691
1692
2/4
✓ Branch 1 taken 471 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 471 times.
 471 if (interrupter.wasInterrupted()) {
1693 ✗ OPENVDB_THROW(RuntimeError, "conjugate gradient solver was interrupted");
1694 }
1695
1696 OPENVDB_LOG_DEBUG_RUNTIME("pcg::solve() " << result);
1697
1698 471 result.iterations = iteration + 1;
1699
1700 // Apply preconditioner to residual
1701 // z_{k} = M^-1 r_{k}
1702
1/2
✓ Branch 1 taken 471 times.
✗ Branch 2 not taken.
 471 precond.apply(rVec, zVec);
1703
1704 // <r,z>
1705
1/2
✓ Branch 1 taken 471 times.
✗ Branch 2 not taken.
 471 const ValueType rDotZ = rVec.dot(zVec);
1706
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 471 times.
 471 assert(std::isfinite(rDotZ));
1707
1708
2/2
✓ Branch 0 taken 9 times.
✓ Branch 1 taken 462 times.
 471 if (0 == iteration) {
1709 // Initialize
1710
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 pVec = zVec;
1711 } else {
1712 462 const ValueType beta = rDotZ / rDotZPrev;
1713 // p = beta * p + z
1714
1/2
✓ Branch 1 taken 462 times.
✗ Branch 2 not taken.
 462 internal::axpy(beta, pVec, zVec, /*result */pVec);
1715 }
1716
1717 // q_{k} = A p_{k}
1718
1/2
✓ Branch 1 taken 471 times.
✗ Branch 2 not taken.
 471 Amat.vectorMultiply(pVec, qVec);
1719
1720 // alpha = <r_{k-1}, z_{k-1}> / <p_{k},q_{k}>
1721
1/2
✓ Branch 1 taken 471 times.
✗ Branch 2 not taken.
 471 const ValueType pAp = pVec.dot(qVec);
1722
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 471 times.
 471 assert(std::isfinite(pAp));
1723
1724 471 const ValueType alpha = rDotZ / pAp;
1725 rDotZPrev = rDotZ;
1726
1727 // x_{k} = x_{k-1} + alpha * p_{k}
1728
1/2
✓ Branch 1 taken 471 times.
✗ Branch 2 not taken.
 471 internal::axpy(alpha, pVec, xVec, /*result=*/xVec);
1729
1730 // r_{k} = r_{k-1} - alpha_{k-1} A p_{k}
1731
2/6
✓ Branch 1 taken 471 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 471 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 471 internal::axpy(-alpha, qVec, rVec, /*result=*/rVec);
1732
1733 // update tolerances
1734
2/2
✓ Branch 0 taken 99 times.
✓ Branch 1 taken 372 times.
 471 l2Error = rVec.l2Norm();
1735 471 minL2Error = Min(l2Error, minL2Error);
1736
1737
1/2
✓ Branch 1 taken 471 times.
✗ Branch 2 not taken.
 471 result.absoluteError = static_cast<double>(rVec.infNorm());
1738 471 result.relativeError = static_cast<double>(result.absoluteError / infNormOfB);
1739
1740
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 471 times.
 471 if (l2Error > 2 * minL2Error) {
1741 // The solution started to diverge.
1742 ✗ result.success = false;
1743 9 break;
1744 }
1745
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 471 times.
 471 if (!std::isfinite(result.absoluteError)) {
1746 // Total divergence of solution
1747 ✗ result.success = false;
1748 ✗ break;
1749 }
1750
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 465 times.
 471 if (result.absoluteError <= termination.absoluteError) {
1751 // Convergence
1752 6 result.success = true;
1753 6 break;
1754 }
1755
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 462 times.
 465 if (result.relativeError <= termination.relativeError) {
1756 // Convergence
1757 3 result.success = true;
1758 3 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
1772