Loading...
Searching...
No Matches
ConjGradient.h
Go to the documentation of this file.
// Copyright Contributors to the OpenVDB Project
// SPDX-License-Identifier: MPL-2.0
/// @file ConjGradient.h
/// @authors D.J. Hill, Peter Cucka
/// @brief Preconditioned conjugate gradient solver (solves @e Ax = @e b using
/// the conjugate gradient method with one of a selection of preconditioners)
#ifndef OPENVDB_MATH_CONJGRADIENT_HAS_BEEN_INCLUDED
#define OPENVDB_MATH_CONJGRADIENT_HAS_BEEN_INCLUDED
#include <openvdb/Exceptions.h>
#include <openvdb/Types.h>
#include <openvdb/util/logging.h>
#include <openvdb/util/NullInterrupter.h>
#include "Math.h" // for Abs(), isZero(), Max(), Sqrt()
#include <tbb/parallel_for.h>
#include <tbb/parallel_reduce.h>
#include <algorithm> // for std::lower_bound()
#include <cassert>
#include <cmath> // for std::isfinite()
#include <limits>
#include <sstream>
#include <string>
namespace openvdb {
OPENVDB_USE_VERSION_NAMESPACE
namespace OPENVDB_VERSION_NAME {
namespace math {
namespace pcg {
using SizeType = Index32;
using SizeRange = tbb::blocked_range<SizeType>;
template<typename ValueType> class Vector;
template<typename ValueType, SizeType STENCIL_SIZE> class SparseStencilMatrix;
template<typename ValueType> class Preconditioner;
template<typename MatrixType> class JacobiPreconditioner;
template<typename MatrixType> class IncompleteCholeskyPreconditioner;
/// Information about the state of a conjugate gradient solution
struct State {
bool success;
int iterations;
double relativeError;
double absoluteError;
};
/// Return default termination conditions for a conjugate gradient solver.
template<typename ValueType>
inline State
terminationDefaults()
{
State s;
s.success = false;
s.iterations = 50;
s.relativeError = 1.0e-6;
s.absoluteError = std::numeric_limits<ValueType>::epsilon() * 100.0;
return s;
}
////////////////////////////////////////
/// @brief Solve @e Ax = @e b via the preconditioned conjugate gradient method.
///
/// @param A a symmetric, positive-definite, @e N x @e N matrix
/// @param b a vector of size @e N
/// @param x a vector of size @e N
/// @param preconditioner a Preconditioner matrix
/// @param termination termination conditions given as a State object with the following fields:
/// <dl>
/// <dt><i>success</i>
/// <dd>ignored
/// <dt><i>iterations</i>
/// <dd>the maximum number of iterations, with or without convergence
/// <dt><i>relativeError</i>
/// <dd>the relative error ||<i>b</i> &minus; <i>Ax</i>|| / ||<i>b</i>||
/// that denotes convergence
/// <dt><i>absoluteError</i>
/// <dd>the absolute error ||<i>b</i> &minus; <i>Ax</i>|| that denotes convergence
///
/// @throw ArithmeticError if either @a x or @a b is not of the appropriate size.
template<typename PositiveDefMatrix>
inline 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>());
/// @brief Solve @e Ax = @e b via the preconditioned conjugate gradient method.
///
/// @param A a symmetric, positive-definite, @e N x @e N matrix
/// @param b a vector of size @e N
/// @param x a vector of size @e N
/// @param preconditioner a Preconditioner matrix
/// @param termination termination conditions given as a State object with the following fields:
/// <dl>
/// <dt><i>success</i>
/// <dd>ignored
/// <dt><i>iterations</i>
/// <dd>the maximum number of iterations, with or without convergence
/// <dt><i>relativeError</i>
/// <dd>the relative error ||<i>b</i> &minus; <i>Ax</i>|| / ||<i>b</i>||
/// that denotes convergence
/// <dt><i>absoluteError</i>
/// <dd>the absolute error ||<i>b</i> &minus; <i>Ax</i>|| that denotes convergence
/// @param interrupter an object adhering to the util::NullInterrupter interface
/// with which computation can be interrupted
///
/// @throw ArithmeticError if either @a x or @a b is not of the appropriate size.
/// @throw RuntimeError if the computation is interrupted.
template<typename PositiveDefMatrix, typename Interrupter>
inline State
solve(
const PositiveDefMatrix& A,
const Vector<typename PositiveDefMatrix::ValueType>& b,
Vector<typename PositiveDefMatrix::ValueType>& x,
Preconditioner<typename PositiveDefMatrix::ValueType>& preconditioner,
Interrupter& interrupter,
const State& termination = terminationDefaults<typename PositiveDefMatrix::ValueType>());
////////////////////////////////////////
/// Lightweight, variable-length vector
template<typename T>
class Vector
{
public:
using ValueType = T;
using Ptr = SharedPtr<Vector>;
/// Construct an empty vector.
Vector(): mData(nullptr), mSize(0) {}
/// Construct a vector of @a n elements, with uninitialized values.
Vector(SizeType n): mData(new T[n]), mSize(n) {}
/// Construct a vector of @a n elements and initialize each element to the given value.
Vector(SizeType n, const ValueType& val): mData(new T[n]), mSize(n) { this->fill(val); }
~Vector() { mSize = 0; delete[] mData; mData = nullptr; }
/// Deep copy the given vector.
Vector(const Vector&);
/// Deep copy the given vector.
Vector& operator=(const Vector&);
/// Return the number of elements in this vector.
SizeType size() const { return mSize; }
/// Return @c true if this vector has no elements.
bool empty() const { return (mSize == 0); }
/// @brief Reset this vector to have @a n elements, with uninitialized values.
/// @warning All of this vector's existing values will be lost.
void resize(SizeType n);
/// Swap internal storage with another vector, which need not be the same size.
void swap(Vector& other) { std::swap(mData, other.mData); std::swap(mSize, other.mSize); }
/// Set all elements of this vector to @a value.
void fill(const ValueType& value);
//@{
/// @brief Multiply each element of this vector by @a s.
template<typename Scalar> void scale(const Scalar& s);
template<typename Scalar> Vector& operator*=(const Scalar& s) { this->scale(s); return *this; }
//@}
/// Return the dot product of this vector with the given vector, which must be the same size.
ValueType dot(const Vector&) const;
/// Return the infinity norm of this vector.
ValueType infNorm() const;
/// Return the L2 norm of this vector.
ValueType l2Norm() const { return Sqrt(this->dot(*this)); }
/// Return @c true if every element of this vector has a finite value.
bool isFinite() const;
/// @brief Return @c true if this vector is equivalent to the given vector
/// to within the specified tolerance.
template<typename OtherValueType>
bool eq(const Vector<OtherValueType>& other,
ValueType eps = Tolerance<ValueType>::value()) const;
/// Return a string representation of this vector.
std::string str() const;
//@{
/// @brief Return the value of this vector's ith element.
inline T& at(SizeType i) { return mData[i]; }
inline const T& at(SizeType i) const { return mData[i]; }
inline T& operator[](SizeType i) { return this->at(i); }
inline const T& operator[](SizeType i) const { return this->at(i); }
//@}
//@{
/// @brief Return a pointer to this vector's elements.
inline T* data() { return mData; }
inline const T* data() const { return mData; }
inline const T* constData() const { return mData; }
//@}
private:
// Functor for use with tbb::parallel_for()
template<typename Scalar> struct ScaleOp;
struct DeterministicDotProductOp;
// Functors for use with tbb::parallel_reduce()
template<typename OtherValueType> struct EqOp;
struct InfNormOp;
struct IsFiniteOp;
T* mData;
SizeType mSize;
};
using VectorS = Vector<float>;
using VectorD = Vector<double>;
////////////////////////////////////////
/// @brief Sparse, square matrix representing a 3D stencil operator of size @a STENCIL_SIZE
/// @details The implementation is a variation on compressed row storage (CRS).
template<typename ValueType_, SizeType STENCIL_SIZE>
class SparseStencilMatrix
{
public:
using ValueType = ValueType_;
using VectorType = Vector<ValueType>;
using Ptr = SharedPtr<SparseStencilMatrix>;
class ConstValueIter;
class ConstRow;
class RowEditor;
static const ValueType sZeroValue;
/// Construct an @a n x @a n matrix with at most @a STENCIL_SIZE nonzero elements per row.
SparseStencilMatrix(SizeType n);
/// Deep copy the given matrix.
SparseStencilMatrix(const SparseStencilMatrix&);
//@{
/// Return the number of rows in this matrix.
SizeType numRows() const { return mNumRows; }
SizeType size() const { return mNumRows; }
//@}
/// @brief Set the value at the given coordinates.
/// @warning It is not safe to set values in the same row simultaneously
/// from multiple threads.
void setValue(SizeType row, SizeType col, const ValueType&);
//@{
/// @brief Return the value at the given coordinates.
/// @warning It is not safe to get values from a row while another thread
/// is setting values in that row.
const ValueType& getValue(SizeType row, SizeType col) const;
const ValueType& operator()(SizeType row, SizeType col) const;
//@}
/// Return a read-only view onto the given row of this matrix.
ConstRow getConstRow(SizeType row) const;
/// Return a read/write view onto the given row of this matrix.
RowEditor getRowEditor(SizeType row);
//@{
/// @brief Multiply all elements in the matrix by @a s;
template<typename Scalar> void scale(const Scalar& s);
template<typename Scalar>
SparseStencilMatrix& operator*=(const Scalar& s) { this->scale(s); return *this; }
//@}
/// @brief Multiply this matrix by @a inVec and return the result in @a resultVec.
/// @throw ArithmeticError if either @a inVec or @a resultVec is not of size @e N,
/// where @e N x @e N is the size of this matrix.
template<typename VecValueType>
void vectorMultiply(const Vector<VecValueType>& inVec, Vector<VecValueType>& resultVec) const;
/// @brief Multiply this matrix by the vector represented by the array @a inVec
/// and return the result in @a resultVec.
/// @warning Both @a inVec and @a resultVec must have at least @e N elements,
/// where @e N x @e N is the size of this matrix.
template<typename VecValueType>
void vectorMultiply(const VecValueType* inVec, VecValueType* resultVec) const;
/// @brief Return @c true if this matrix is equivalent to the given matrix
/// to within the specified tolerance.
template<typename OtherValueType>
bool eq(const SparseStencilMatrix<OtherValueType, STENCIL_SIZE>& other,
ValueType eps = Tolerance<ValueType>::value()) const;
/// Return @c true if every element of this matrix has a finite value.
bool isFinite() const;
/// Return a string representation of this matrix.
std::string str() const;
private:
struct RowData {
RowData(ValueType* v, SizeType* c, SizeType& s): mVals(v), mCols(c), mSize(s) {}
ValueType* mVals; SizeType* mCols; SizeType& mSize;
};
struct ConstRowData {
ConstRowData(const ValueType* v, const SizeType* c, const SizeType& s):
mVals(v), mCols(c), mSize(s) {}
const ValueType* mVals; const SizeType* mCols; const SizeType& mSize;
};
/// Base class for row accessors
template<typename DataType_ = RowData>
class RowBase
{
public:
using DataType = DataType_;
static SizeType capacity() { return STENCIL_SIZE; }
RowBase(const DataType& data): mData(data) {}
bool empty() const { return (mData.mSize == 0); }
const SizeType& size() const { return mData.mSize; }
const ValueType& getValue(SizeType columnIdx, bool& active) const;
const ValueType& getValue(SizeType columnIdx) const;
/// Return an iterator over the stored values in this row.
ConstValueIter cbegin() const;
/// @brief Return @c true if this row is equivalent to the given row
/// to within the specified tolerance.
template<typename OtherDataType>
bool eq(const RowBase<OtherDataType>& other,
ValueType eps = Tolerance<ValueType>::value()) const;
/// @brief Return the dot product of this row with the first
/// @a vecSize elements of @a inVec.
/// @warning @a inVec must have at least @a vecSize elements.
template<typename VecValueType>
VecValueType dot(const VecValueType* inVec, SizeType vecSize) const;
/// Return the dot product of this row with the given vector.
template<typename VecValueType>
VecValueType dot(const Vector<VecValueType>& inVec) const;
/// Return a string representation of this row.
std::string str() const;
protected:
friend class ConstValueIter;
const ValueType& value(SizeType i) const { return mData.mVals[i]; }
SizeType column(SizeType i) const { return mData.mCols[i]; }
/// @brief Return the array index of the first column index that is
/// equal to <i>or greater than</i> the given column index.
/// @note If @a columnIdx is larger than any existing column index,
/// the return value will point beyond the end of the array.
SizeType find(SizeType columnIdx) const;
DataType mData;
};
using ConstRowBase = RowBase<ConstRowData>;
public:
/// Iterator over the stored values in a row of this matrix
class ConstValueIter
{
public:
const ValueType& operator*() const
{
if (mData.mSize == 0) return SparseStencilMatrix::sZeroValue;
return mData.mVals[mCursor];
}
SizeType column() const { return mData.mCols[mCursor]; }
void increment() { mCursor++; }
ConstValueIter& operator++() { increment(); return *this; }
operator bool() const { return (mCursor < mData.mSize); }
void reset() { mCursor = 0; }
private:
friend class SparseStencilMatrix;
ConstValueIter(const RowData& d): mData(d.mVals, d.mCols, d.mSize), mCursor(0) {}
ConstValueIter(const ConstRowData& d): mData(d), mCursor(0) {}
const ConstRowData mData;
SizeType mCursor;
};
/// Read-only accessor to a row of this matrix
class ConstRow: public ConstRowBase
{
public:
ConstRow(const ValueType* valueHead, const SizeType* columnHead, const SizeType& rowSize);
}; // class ConstRow
/// Read/write accessor to a row of this matrix
class RowEditor: public RowBase<>
{
public:
RowEditor(ValueType* valueHead, SizeType* columnHead, SizeType& rowSize, SizeType colSize);
/// Set the number of entries in this row to zero.
void clear();
/// @brief Set the value of the entry in the specified column.
/// @return the current number of entries stored in this row.
SizeType setValue(SizeType column, const ValueType& value);
//@{
/// @brief Scale all of the entries in this row.
template<typename Scalar> void scale(const Scalar&);
template<typename Scalar>
RowEditor& operator*=(const Scalar& s) { this->scale(s); return *this; }
//@}
private:
const SizeType mNumColumns; // used only for bounds checking
}; // class RowEditor
private:
// Functors for use with tbb::parallel_for()
struct MatrixCopyOp;
template<typename VecValueType> struct VecMultOp;
template<typename Scalar> struct RowScaleOp;
// Functors for use with tbb::parallel_reduce()
struct IsFiniteOp;
template<typename OtherValueType> struct EqOp;
const SizeType mNumRows;
std::unique_ptr<ValueType[]> mValueArray;
std::unique_ptr<SizeType[]> mColumnIdxArray;
std::unique_ptr<SizeType[]> mRowSizeArray;
}; // class SparseStencilMatrix
////////////////////////////////////////
/// Base class for conjugate gradient preconditioners
template<typename T>
class Preconditioner
{
public:
using ValueType = T;
using Ptr = SharedPtr<Preconditioner>;
template<SizeType STENCIL_SIZE> Preconditioner(const SparseStencilMatrix<T, STENCIL_SIZE>&) {}
virtual ~Preconditioner() = default;
virtual bool isValid() const { return true; }
/// @brief Apply this preconditioner to a residue vector:
/// @e z = <i>M</i><sup><small>&minus;1</small></sup><i>r</i>
/// @param r residue vector
/// @param[out] z preconditioned residue vector
virtual void apply(const Vector<T>& r, Vector<T>& z) = 0;
};
////////////////////////////////////////
namespace internal {
// Functor for use with tbb::parallel_for() to copy data from one array to another
template<typename T>
struct CopyOp
{
CopyOp(const T* from_, T* to_): from(from_), to(to_) {}
void operator()(const SizeRange& range) const {
for (SizeType n = range.begin(), N = range.end(); n < N; ++n) to[n] = from[n];
}
const T* from;
T* to;
};
// Functor for use with tbb::parallel_for() to fill an array with a constant value
template<typename T>
struct FillOp
{
FillOp(T* data_, const T& val_): data(data_), val(val_) {}
void operator()(const SizeRange& range) const {
for (SizeType n = range.begin(), N = range.end(); n < N; ++n) data[n] = val;
}
T* data;
const T val;
};
// Functor for use with tbb::parallel_for() that computes a * x + y
template<typename T>
struct LinearOp
{
LinearOp(const T& a_, const T* x_, const T* y_, T* out_): a(a_), x(x_), y(y_), out(out_) {}
void operator()(const SizeRange& range) const {
if (isExactlyEqual(a, T(1))) {
for (SizeType n = range.begin(), N = range.end(); n < N; ++n) out[n] = x[n] + y[n];
} else if (isExactlyEqual(a, T(-1))) {
for (SizeType n = range.begin(), N = range.end(); n < N; ++n) out[n] = -x[n] + y[n];
} else {
for (SizeType n = range.begin(), N = range.end(); n < N; ++n) out[n] = a * x[n] + y[n];
}
}
const T a, *x, *y;
T* out;
};
} // namespace internal
////////////////////////////////////////
inline std::ostream&
operator<<(std::ostream& os, const State& state)
{
os << (state.success ? "succeeded with " : "")
<< "rel. err. " << state.relativeError << ", abs. err. " << state.absoluteError
<< " after " << state.iterations << " iteration" << (state.iterations == 1 ? "" : "s");
return os;
}
////////////////////////////////////////
template<typename T>
inline
Vector<T>::Vector(const Vector& other): mData(new T[other.mSize]), mSize(other.mSize)
{
tbb::parallel_for(SizeRange(0, mSize),
internal::CopyOp<T>(/*from=*/other.mData, /*to=*/mData));
}
template<typename T>
inline
Vector<T>& Vector<T>::operator=(const Vector<T>& other)
{
// Update the internal storage to the correct size
if (mSize != other.mSize) {
mSize = other.mSize;
delete[] mData;
mData = new T[mSize];
}
// Deep copy the data
tbb::parallel_for(SizeRange(0, mSize),
internal::CopyOp<T>(/*from=*/other.mData, /*to=*/mData));
return *this;
}
template<typename T>
inline void
Vector<T>::resize(SizeType n)
{
if (n != mSize) {
delete[] mData;
mData = new T[n];
mSize = n;
}
}
template<typename T>
inline void
Vector<T>::fill(const ValueType& value)
{
tbb::parallel_for(SizeRange(0, mSize), internal::FillOp<T>(mData, value));
}
template<typename T>
template<typename Scalar>
struct Vector<T>::ScaleOp
{
ScaleOp(T* data_, const Scalar& s_): data(data_), s(s_) {}
void operator()(const SizeRange& range) const {
for (SizeType n = range.begin(), N = range.end(); n < N; ++n) data[n] *= s;
}
T* data;
const Scalar s;
};
template<typename T>
template<typename Scalar>
inline void
Vector<T>::scale(const Scalar& s)
{
tbb::parallel_for(SizeRange(0, mSize), ScaleOp<Scalar>(mData, s));
}
template<typename T>
struct Vector<T>::DeterministicDotProductOp
{
DeterministicDotProductOp(const T* a_, const T* b_,
const SizeType binCount_, const SizeType arraySize_, T* reducetmp_):
a(a_), b(b_), binCount(binCount_), arraySize(arraySize_), reducetmp(reducetmp_) {}
void operator()(const SizeRange& range) const
{
const SizeType binSize = arraySize / binCount;
// Iterate over bins (array segments)
for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
const SizeType begin = n * binSize;
const SizeType end = (n == binCount-1) ? arraySize : begin + binSize;
// Compute the partial sum for this array segment
T sum = zeroVal<T>();
for (SizeType i = begin; i < end; ++i) { sum += a[i] * b[i]; }
// Store the partial sum
reducetmp[n] = sum;
}
}
const T* a;
const T* b;
const SizeType binCount;
const SizeType arraySize;
T* reducetmp;
};
template<typename T>
inline T
Vector<T>::dot(const Vector<T>& other) const
{
assert(this->size() == other.size());
const T* aData = this->data();
const T* bData = other.data();
SizeType arraySize = this->size();
T result = zeroVal<T>();
if (arraySize < 1024) {
// Compute the dot product in serial for small arrays
for (SizeType n = 0; n < arraySize; ++n) {
result += aData[n] * bData[n];
}
} else {
// Compute the dot product by segmenting the arrays into
// a predetermined number of sub arrays in parallel and
// accumulate the finial result in series.
const SizeType binCount = 100;
T partialSums[100];
tbb::parallel_for(SizeRange(0, binCount),
DeterministicDotProductOp(aData, bData, binCount, arraySize, partialSums));
for (SizeType n = 0; n < binCount; ++n) {
result += partialSums[n];
}
}
return result;
}
template<typename T>
struct Vector<T>::InfNormOp
{
InfNormOp(const T* data_): data(data_) {}
T operator()(const SizeRange& range, T maxValue) const
{
for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
maxValue = Max(maxValue, Abs(data[n]));
}
return maxValue;
}
const T* data;
};
template<typename T>
inline T
Vector<T>::infNorm() const
{
// Parallelize over the elements of this vector.
T result = tbb::parallel_reduce(SizeRange(0, this->size()), /*seed=*/zeroVal<T>(),
InfNormOp(this->data()), /*join=*/[](T max1, T max2) { return Max(max1, max2); });
return result;
}
template<typename T>
struct Vector<T>::IsFiniteOp
{
IsFiniteOp(const T* data_): data(data_) {}
bool operator()(const SizeRange& range, bool finite) const
{
if (finite) {
for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
if (!std::isfinite(data[n])) return false;
}
}
return finite;
}
const T* data;
};
template<typename T>
inline bool
Vector<T>::isFinite() const
{
// Parallelize over the elements of this vector.
bool finite = tbb::parallel_reduce(SizeRange(0, this->size()), /*seed=*/true,
IsFiniteOp(this->data()),
/*join=*/[](bool finite1, bool finite2) { return (finite1 && finite2); });
return finite;
}
template<typename T>
template<typename OtherValueType>
struct Vector<T>::EqOp
{
EqOp(const T* a_, const OtherValueType* b_, T e): a(a_), b(b_), eps(e) {}
bool operator()(const SizeRange& range, bool equal) const
{
if (equal) {
for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
if (!isApproxEqual(a[n], b[n], eps)) return false;
}
}
return equal;
}
const T* a;
const OtherValueType* b;
const T eps;
};
template<typename T>
template<typename OtherValueType>
inline bool
Vector<T>::eq(const Vector<OtherValueType>& other, ValueType eps) const
{
if (this->size() != other.size()) return false;
bool equal = tbb::parallel_reduce(SizeRange(0, this->size()), /*seed=*/true,
EqOp<OtherValueType>(this->data(), other.data(), eps),
/*join=*/[](bool eq1, bool eq2) { return (eq1 && eq2); });
return equal;
}
template<typename T>
inline std::string
Vector<T>::str() const
{
std::ostringstream ostr;
ostr << "[";
std::string sep;
for (SizeType n = 0, N = this->size(); n < N; ++n) {
ostr << sep << (*this)[n];
sep = ", ";
}
ostr << "]";
return ostr.str();
}
////////////////////////////////////////
template<typename ValueType, SizeType STENCIL_SIZE>
const ValueType SparseStencilMatrix<ValueType, STENCIL_SIZE>::sZeroValue = zeroVal<ValueType>();
template<typename ValueType, SizeType STENCIL_SIZE>
inline
SparseStencilMatrix<ValueType, STENCIL_SIZE>::SparseStencilMatrix(SizeType numRows)
: mNumRows(numRows)
, mValueArray(new ValueType[mNumRows * STENCIL_SIZE])
, mColumnIdxArray(new SizeType[mNumRows * STENCIL_SIZE])
, mRowSizeArray(new SizeType[mNumRows])
{
// Initialize the matrix to a null state by setting the size of each row to zero.
tbb::parallel_for(SizeRange(0, mNumRows),
internal::FillOp<SizeType>(mRowSizeArray.get(), /*value=*/0));
}
template<typename ValueType, SizeType STENCIL_SIZE>
struct SparseStencilMatrix<ValueType, STENCIL_SIZE>::MatrixCopyOp
{
MatrixCopyOp(const SparseStencilMatrix& from_, SparseStencilMatrix& to_):
from(&from_), to(&to_) {}
void operator()(const SizeRange& range) const
{
const ValueType* fromVal = from->mValueArray.get();
const SizeType* fromCol = from->mColumnIdxArray.get();
ValueType* toVal = to->mValueArray.get();
SizeType* toCol = to->mColumnIdxArray.get();
for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
toVal[n] = fromVal[n];
toCol[n] = fromCol[n];
}
}
const SparseStencilMatrix* from; SparseStencilMatrix* to;
};
template<typename ValueType, SizeType STENCIL_SIZE>
inline
SparseStencilMatrix<ValueType, STENCIL_SIZE>::SparseStencilMatrix(const SparseStencilMatrix& other)
: mNumRows(other.mNumRows)
, mValueArray(new ValueType[mNumRows * STENCIL_SIZE])
, mColumnIdxArray(new SizeType[mNumRows * STENCIL_SIZE])
, mRowSizeArray(new SizeType[mNumRows])
{
SizeType size = mNumRows * STENCIL_SIZE;
// Copy the value and column index arrays from the other matrix to this matrix.
tbb::parallel_for(SizeRange(0, size), MatrixCopyOp(/*from=*/other, /*to=*/*this));
// Copy the row size array from the other matrix to this matrix.
tbb::parallel_for(SizeRange(0, mNumRows),
internal::CopyOp<SizeType>(/*from=*/other.mRowSizeArray.get(), /*to=*/mRowSizeArray.get()));
}
template<typename ValueType, SizeType STENCIL_SIZE>
inline void
SparseStencilMatrix<ValueType, STENCIL_SIZE>::setValue(SizeType row, SizeType col,
const ValueType& val)
{
assert(row < mNumRows);
this->getRowEditor(row).setValue(col, val);
}
template<typename ValueType, SizeType STENCIL_SIZE>
inline const ValueType&
SparseStencilMatrix<ValueType, STENCIL_SIZE>::getValue(SizeType row, SizeType col) const
{
assert(row < mNumRows);
return this->getConstRow(row).getValue(col);
}
template<typename ValueType, SizeType STENCIL_SIZE>
inline const ValueType&
SparseStencilMatrix<ValueType, STENCIL_SIZE>::operator()(SizeType row, SizeType col) const
{
return this->getValue(row,col);
}
template<typename ValueType, SizeType STENCIL_SIZE>
template<typename Scalar>
struct SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowScaleOp
{
RowScaleOp(SparseStencilMatrix& m, const Scalar& s_): mat(&m), s(s_) {}
void operator()(const SizeRange& range) const
{
for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
RowEditor row = mat->getRowEditor(n);
row.scale(s);
}
}
SparseStencilMatrix* mat;
const Scalar s;
};
template<typename ValueType, SizeType STENCIL_SIZE>
template<typename Scalar>
inline void
SparseStencilMatrix<ValueType, STENCIL_SIZE>::scale(const Scalar& s)
{
// Parallelize over the rows in the matrix.
tbb::parallel_for(SizeRange(0, mNumRows), RowScaleOp<Scalar>(*this, s));
}
template<typename ValueType, SizeType STENCIL_SIZE>
template<typename VecValueType>
struct SparseStencilMatrix<ValueType, STENCIL_SIZE>::VecMultOp
{
VecMultOp(const SparseStencilMatrix& m, const VecValueType* i, VecValueType* o):
mat(&m), in(i), out(o) {}
void operator()(const SizeRange& range) const
{
for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
ConstRow row = mat->getConstRow(n);
out[n] = row.dot(in, mat->numRows());
}
}
const SparseStencilMatrix* mat;
const VecValueType* in;
VecValueType* out;
};
template<typename ValueType, SizeType STENCIL_SIZE>
template<typename VecValueType>
inline void
SparseStencilMatrix<ValueType, STENCIL_SIZE>::vectorMultiply(
const Vector<VecValueType>& inVec, Vector<VecValueType>& resultVec) const
{
if (inVec.size() != mNumRows) {
OPENVDB_THROW(ArithmeticError, "matrix and input vector have incompatible sizes ("
<< mNumRows << "x" << mNumRows << " vs. " << inVec.size() << ")");
}
if (resultVec.size() != mNumRows) {
OPENVDB_THROW(ArithmeticError, "matrix and result vector have incompatible sizes ("
<< mNumRows << "x" << mNumRows << " vs. " << resultVec.size() << ")");
}
vectorMultiply(inVec.data(), resultVec.data());
}
template<typename ValueType, SizeType STENCIL_SIZE>
template<typename VecValueType>
inline void
SparseStencilMatrix<ValueType, STENCIL_SIZE>::vectorMultiply(
const VecValueType* inVec, VecValueType* resultVec) const
{
// Parallelize over the rows in the matrix.
tbb::parallel_for(SizeRange(0, mNumRows),
VecMultOp<VecValueType>(*this, inVec, resultVec));
}
template<typename ValueType, SizeType STENCIL_SIZE>
template<typename OtherValueType>
struct SparseStencilMatrix<ValueType, STENCIL_SIZE>::EqOp
{
EqOp(const SparseStencilMatrix& a_,
const SparseStencilMatrix<OtherValueType, STENCIL_SIZE>& b_, ValueType e):
a(&a_), b(&b_), eps(e) {}
bool operator()(const SizeRange& range, bool equal) const
{
if (equal) {
for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
if (!a->getConstRow(n).eq(b->getConstRow(n), eps)) return false;
}
}
return equal;
}
const SparseStencilMatrix* a;
const SparseStencilMatrix<OtherValueType, STENCIL_SIZE>* b;
const ValueType eps;
};
template<typename ValueType, SizeType STENCIL_SIZE>
template<typename OtherValueType>
inline bool
SparseStencilMatrix<ValueType, STENCIL_SIZE>::eq(
const SparseStencilMatrix<OtherValueType, STENCIL_SIZE>& other, ValueType eps) const
{
if (this->numRows() != other.numRows()) return false;
bool equal = tbb::parallel_reduce(SizeRange(0, this->numRows()), /*seed=*/true,
EqOp<OtherValueType>(*this, other, eps),
/*join=*/[](bool eq1, bool eq2) { return (eq1 && eq2); });
return equal;
}
template<typename ValueType, SizeType STENCIL_SIZE>
struct SparseStencilMatrix<ValueType, STENCIL_SIZE>::IsFiniteOp
{
IsFiniteOp(const SparseStencilMatrix& m): mat(&m) {}
bool operator()(const SizeRange& range, bool finite) const
{
if (finite) {
for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
const ConstRow row = mat->getConstRow(n);
for (ConstValueIter it = row.cbegin(); it; ++it) {
if (!std::isfinite(*it)) return false;
}
}
}
return finite;
}
const SparseStencilMatrix* mat;
};
template<typename ValueType, SizeType STENCIL_SIZE>
inline bool
SparseStencilMatrix<ValueType, STENCIL_SIZE>::isFinite() const
{
// Parallelize over the rows of this matrix.
bool finite = tbb::parallel_reduce(SizeRange(0, this->numRows()), /*seed=*/true,
IsFiniteOp(*this), /*join=*/[](bool finite1, bool finite2) { return (finite1&&finite2); });
return finite;
}
template<typename ValueType, SizeType STENCIL_SIZE>
inline std::string
SparseStencilMatrix<ValueType, STENCIL_SIZE>::str() const
{
std::ostringstream ostr;
for (SizeType n = 0, N = this->size(); n < N; ++n) {
ostr << n << ": " << this->getConstRow(n).str() << "\n";
}
return ostr.str();
}
template<typename ValueType, SizeType STENCIL_SIZE>
inline typename SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowEditor
SparseStencilMatrix<ValueType, STENCIL_SIZE>::getRowEditor(SizeType i)
{
assert(i < mNumRows);
const SizeType head = i * STENCIL_SIZE;
return RowEditor(&mValueArray[head], &mColumnIdxArray[head], mRowSizeArray[i], mNumRows);
}
template<typename ValueType, SizeType STENCIL_SIZE>
inline typename SparseStencilMatrix<ValueType, STENCIL_SIZE>::ConstRow
SparseStencilMatrix<ValueType, STENCIL_SIZE>::getConstRow(SizeType i) const
{
assert(i < mNumRows);
const SizeType head = i * STENCIL_SIZE; // index for this row into main storage
return ConstRow(&mValueArray[head], &mColumnIdxArray[head], mRowSizeArray[i]);
}
template<typename ValueType, SizeType STENCIL_SIZE>
template<typename DataType>
inline SizeType
SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::find(SizeType columnIdx) const
{
if (this->empty()) return mData.mSize;
// Get a pointer to the first column index that is equal to or greater than the given index.
// (This assumes that the data is sorted by column.)
const SizeType* colPtr = std::lower_bound(mData.mCols, mData.mCols + mData.mSize, columnIdx);
// Return the offset of the pointer from the beginning of the array.
return static_cast<SizeType>(colPtr - mData.mCols);
}
template<typename ValueType, SizeType STENCIL_SIZE>
template<typename DataType>
inline const ValueType&
SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::getValue(
SizeType columnIdx, bool& active) const
{
active = false;
SizeType idx = this->find(columnIdx);
if (idx < this->size() && this->column(idx) == columnIdx) {
active = true;
return this->value(idx);
}
return SparseStencilMatrix::sZeroValue;
}
template<typename ValueType, SizeType STENCIL_SIZE>
template<typename DataType>
inline const ValueType&
SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::getValue(SizeType columnIdx) const
{
SizeType idx = this->find(columnIdx);
if (idx < this->size() && this->column(idx) == columnIdx) {
return this->value(idx);
}
return SparseStencilMatrix::sZeroValue;
}
template<typename ValueType, SizeType STENCIL_SIZE>
template<typename DataType>
inline typename SparseStencilMatrix<ValueType, STENCIL_SIZE>::ConstValueIter
SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::cbegin() const
{
return ConstValueIter(mData);
}
template<typename ValueType, SizeType STENCIL_SIZE>
template<typename DataType>
template<typename OtherDataType>
inline bool
SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::eq(
const RowBase<OtherDataType>& other, ValueType eps) const
{
if (this->size() != other.size()) return false;
for (ConstValueIter it = cbegin(), oit = other.cbegin(); it || oit; ++it, ++oit) {
if (it.column() != oit.column()) return false;
if (!isApproxEqual(*it, *oit, eps)) return false;
}
return true;
}
template<typename ValueType, SizeType STENCIL_SIZE>
template<typename DataType>
template<typename VecValueType>
inline VecValueType
SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::dot(
const VecValueType* inVec, SizeType vecSize) const
{
VecValueType result = zeroVal<VecValueType>();
for (SizeType idx = 0, N = std::min(vecSize, this->size()); idx < N; ++idx) {
result += static_cast<VecValueType>(this->value(idx) * inVec[this->column(idx)]);
}
return result;
}
template<typename ValueType, SizeType STENCIL_SIZE>
template<typename DataType>
template<typename VecValueType>
inline VecValueType
SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::dot(
const Vector<VecValueType>& inVec) const
{
return dot(inVec.data(), inVec.size());
}
template<typename ValueType, SizeType STENCIL_SIZE>
template<typename DataType>
inline std::string
SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::str() const
{
std::ostringstream ostr;
std::string sep;
for (SizeType n = 0, N = this->size(); n < N; ++n) {
ostr << sep << "(" << this->column(n) << ", " << this->value(n) << ")";
sep = ", ";
}
return ostr.str();
}
template<typename ValueType, SizeType STENCIL_SIZE>
inline
SparseStencilMatrix<ValueType, STENCIL_SIZE>::ConstRow::ConstRow(
const ValueType* valueHead, const SizeType* columnHead, const SizeType& rowSize):
ConstRowBase(ConstRowData(const_cast<ValueType*>(valueHead),
const_cast<SizeType*>(columnHead), const_cast<SizeType&>(rowSize)))
{
}
template<typename ValueType, SizeType STENCIL_SIZE>
inline
SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowEditor::RowEditor(
ValueType* valueHead, SizeType* columnHead, SizeType& rowSize, SizeType colSize):
RowBase<>(RowData(valueHead, columnHead, rowSize)), mNumColumns(colSize)
{
}
template<typename ValueType, SizeType STENCIL_SIZE>
inline void
SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowEditor::clear()
{
// Note: since mSize is a reference, this modifies the underlying matrix.
RowBase<>::mData.mSize = 0;
}
template<typename ValueType, SizeType STENCIL_SIZE>
inline SizeType
SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowEditor::setValue(
SizeType column, const ValueType& value)
{
assert(column < mNumColumns);
RowData& data = RowBase<>::mData;
// Get the offset of the first column index that is equal to or greater than
// the column to be modified.
SizeType offset = this->find(column);
if (offset < data.mSize && data.mCols[offset] == column) {
// If the column already exists, just update its value.
data.mVals[offset] = value;
return data.mSize;
}
// Check that it is safe to add a new column.
assert(data.mSize < this->capacity());
if (offset >= data.mSize) {
// The new column's index is larger than any existing index. Append the new column.
data.mVals[data.mSize] = value;
data.mCols[data.mSize] = column;
} else {
// Insert the new column at the computed offset after shifting subsequent columns.
for (SizeType i = data.mSize; i > offset; --i) {
data.mVals[i] = data.mVals[i - 1];
data.mCols[i] = data.mCols[i - 1];
}
data.mVals[offset] = value;
data.mCols[offset] = column;
}
++data.mSize;
return data.mSize;
}
template<typename ValueType, SizeType STENCIL_SIZE>
template<typename Scalar>
inline void
SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowEditor::scale(const Scalar& s)
{
for (int idx = 0, N = this->size(); idx < N; ++idx) {
RowBase<>::mData.mVals[idx] *= s;
}
}
////////////////////////////////////////
/// Diagonal preconditioner
template<typename MatrixType>
class JacobiPreconditioner: public Preconditioner<typename MatrixType::ValueType>
{
private:
struct InitOp;
struct ApplyOp;
public:
using ValueType = typename MatrixType::ValueType;
using BaseType = Preconditioner<ValueType>;
using VectorType = Vector<ValueType>;
using Ptr = SharedPtr<JacobiPreconditioner>;
JacobiPreconditioner(const MatrixType& A): BaseType(A), mDiag(A.numRows())
{
// Initialize vector mDiag with the values from the matrix diagonal.
tbb::parallel_for(SizeRange(0, A.numRows()), InitOp(A, mDiag.data()));
}
~JacobiPreconditioner() override = default;
void apply(const Vector<ValueType>& r, Vector<ValueType>& z) override
{
const SizeType size = mDiag.size();
assert(r.size() == z.size());
assert(r.size() == size);
tbb::parallel_for(SizeRange(0, size), ApplyOp(mDiag.data(), r.data(), z.data()));
}
/// Return @c true if all values along the diagonal are finite.
bool isFinite() const { return mDiag.isFinite(); }
private:
// Functor for use with tbb::parallel_for()
struct InitOp
{
InitOp(const MatrixType& m, ValueType* v): mat(&m), vec(v) {}
void operator()(const SizeRange& range) const {
for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
const ValueType val = mat->getValue(n, n);
assert(!isApproxZero(val, ValueType(0.0001)));
vec[n] = static_cast<ValueType>(1.0 / val);
}
}
const MatrixType* mat; ValueType* vec;
};
// Functor for use with tbb::parallel_reduce()
struct ApplyOp
{
ApplyOp(const ValueType* x_, const ValueType* y_, ValueType* out_):
x(x_), y(y_), out(out_) {}
void operator()(const SizeRange& range) const {
for (SizeType n = range.begin(), N = range.end(); n < N; ++n) out[n] = x[n] * y[n];
}
const ValueType *x, *y; ValueType* out;
};
// The Jacobi preconditioner is a diagonal matrix
VectorType mDiag;
}; // class JacobiPreconditioner
/// Preconditioner using incomplete Cholesky factorization
template<typename MatrixType>
class IncompleteCholeskyPreconditioner: public Preconditioner<typename MatrixType::ValueType>
{
private:
struct CopyToLowerOp;
struct TransposeOp;
public:
using ValueType = typename MatrixType::ValueType;
using BaseType = Preconditioner<ValueType>;
using VectorType = Vector<ValueType>;
using Ptr = SharedPtr<IncompleteCholeskyPreconditioner>;
using TriangularMatrix = SparseStencilMatrix<ValueType, 4>;
using TriangleConstRow = typename TriangularMatrix::ConstRow;
using TriangleRowEditor = typename TriangularMatrix::RowEditor;
IncompleteCholeskyPreconditioner(const MatrixType& matrix)
: BaseType(matrix)
, mLowerTriangular(matrix.numRows())
, mUpperTriangular(matrix.numRows())
, mTempVec(matrix.numRows())
{
// Size of matrix
const SizeType numRows = mLowerTriangular.numRows();
// Copy the upper triangular part to the lower triangular part.
tbb::parallel_for(SizeRange(0, numRows), CopyToLowerOp(matrix, mLowerTriangular));
// Build the Incomplete Cholesky Matrix
//
// Algorithm:
//
// for (k = 0; k < size; ++k) {
// A(k,k) = sqrt(A(k,k));
// for (i = k +1, i < size; ++i) {
// if (A(i,k) == 0) continue;
// A(i,k) = A(i,k) / A(k,k);
// }
// for (j = k+1; j < size; ++j) {
// for (i = j; i < size; ++i) {
// if (A(i,j) == 0) continue;
// A(i,j) -= A(i,k)*A(j,k);
// }
// }
// }
mPassedCompatibilityCondition = true;
for (SizeType k = 0; k < numRows; ++k) {
TriangleConstRow crow_k = mLowerTriangular.getConstRow(k);
ValueType diagonalValue = crow_k.getValue(k);
// Test if the matrix build has failed.
if (diagonalValue < 1.e-5) {
mPassedCompatibilityCondition = false;
break;
}
diagonalValue = Sqrt(diagonalValue);
TriangleRowEditor row_k = mLowerTriangular.getRowEditor(k);
row_k.setValue(k, diagonalValue);
// Exploit the fact that the matrix is symmetric.
typename MatrixType::ConstRow srcRow = matrix.getConstRow(k);
typename MatrixType::ConstValueIter citer = srcRow.cbegin();
for ( ; citer; ++citer) {
SizeType ii = citer.column();
if (ii < k+1) continue; // look above diagonal
TriangleRowEditor row_ii = mLowerTriangular.getRowEditor(ii);
row_ii.setValue(k, *citer / diagonalValue);
}
// for (j = k+1; j < size; ++j) replaced by row iter below
citer.reset(); // k,j entries
for ( ; citer; ++citer) {
SizeType j = citer.column();
if (j < k+1) continue;
TriangleConstRow row_j = mLowerTriangular.getConstRow(j);
ValueType a_jk = row_j.getValue(k); // a_jk is non zero if a_kj is non zero
// Entry (i,j) is non-zero if matrix(j,i) is nonzero
typename MatrixType::ConstRow mask = matrix.getConstRow(j);
typename MatrixType::ConstValueIter maskIter = mask.cbegin();
for ( ; maskIter; ++maskIter) {
SizeType i = maskIter.column();
if (i < j) continue;
TriangleConstRow crow_i = mLowerTriangular.getConstRow(i);
ValueType a_ij = crow_i.getValue(j);
ValueType a_ik = crow_i.getValue(k);
TriangleRowEditor row_i = mLowerTriangular.getRowEditor(i);
a_ij -= a_ik * a_jk;
row_i.setValue(j, a_ij);
}
}
}
// Build the transpose of the IC matrix: mUpperTriangular
tbb::parallel_for(SizeRange(0, numRows),
TransposeOp(matrix, mLowerTriangular, mUpperTriangular));
}
~IncompleteCholeskyPreconditioner() override = default;
bool isValid() const override { return mPassedCompatibilityCondition; }
void apply(const Vector<ValueType>& rVec, Vector<ValueType>& zVec) override
{
if (!mPassedCompatibilityCondition) {
OPENVDB_THROW(ArithmeticError, "invalid Cholesky decomposition");
}
// Solve mUpperTriangular * mLowerTriangular * rVec = zVec;
SizeType size = mLowerTriangular.numRows();
zVec.fill(zeroVal<ValueType>());
ValueType* zData = zVec.data();
if (size == 0) return;
assert(rVec.size() == size);
assert(zVec.size() == size);
// Allocate a temp vector
mTempVec.fill(zeroVal<ValueType>());
ValueType* tmpData = mTempVec.data();
const ValueType* rData = rVec.data();
// Solve mLowerTriangular * tmp = rVec;
for (SizeType i = 0; i < size; ++i) {
typename TriangularMatrix::ConstRow row = mLowerTriangular.getConstRow(i);
ValueType diagonal = row.getValue(i);
ValueType dot = row.dot(mTempVec);
tmpData[i] = (rData[i] - dot) / diagonal;
if (!std::isfinite(tmpData[i])) {
OPENVDB_LOG_DEBUG_RUNTIME("1 diagonal was " << diagonal);
OPENVDB_LOG_DEBUG_RUNTIME("1a diagonal " << row.getValue(i));
}
}
// Solve mUpperTriangular * zVec = tmp;
for (SizeType ii = 0; ii < size; ++ii) {
SizeType i = size - 1 - ii;
typename TriangularMatrix::ConstRow row = mUpperTriangular.getConstRow(i);
ValueType diagonal = row.getValue(i);
ValueType dot = row.dot(zVec);
zData[i] = (tmpData[i] - dot) / diagonal;
if (!std::isfinite(zData[i])) {
OPENVDB_LOG_DEBUG_RUNTIME("2 diagonal was " << diagonal);
}
}
}
const TriangularMatrix& lowerMatrix() const { return mLowerTriangular; }
const TriangularMatrix& upperMatrix() const { return mUpperTriangular; }
private:
// Functor for use with tbb::parallel_for()
struct CopyToLowerOp
{
CopyToLowerOp(const MatrixType& m, TriangularMatrix& l): mat(&m), lower(&l) {}
void operator()(const SizeRange& range) const {
for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
typename TriangularMatrix::RowEditor outRow = lower->getRowEditor(n);
outRow.clear();
typename MatrixType::ConstRow inRow = mat->getConstRow(n);
for (typename MatrixType::ConstValueIter it = inRow.cbegin(); it; ++it) {
if (it.column() > n) continue; // skip above diagonal
outRow.setValue(it.column(), *it);
}
}
}
const MatrixType* mat; TriangularMatrix* lower;
};
// Functor for use with tbb::parallel_for()
struct TransposeOp
{
TransposeOp(const MatrixType& m, const TriangularMatrix& l, TriangularMatrix& u):
mat(&m), lower(&l), upper(&u) {}
void operator()(const SizeRange& range) const {
for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
typename TriangularMatrix::RowEditor outRow = upper->getRowEditor(n);
outRow.clear();
// Use the fact that matrix is symmetric.
typename MatrixType::ConstRow inRow = mat->getConstRow(n);
for (typename MatrixType::ConstValueIter it = inRow.cbegin(); it; ++it) {
const SizeType column = it.column();
if (column < n) continue; // only set upper triangle
outRow.setValue(column, lower->getValue(column, n));
}
}
}
const MatrixType* mat; const TriangularMatrix* lower; TriangularMatrix* upper;
};
TriangularMatrix mLowerTriangular;
TriangularMatrix mUpperTriangular;
Vector<ValueType> mTempVec;
bool mPassedCompatibilityCondition;
}; // class IncompleteCholeskyPreconditioner
////////////////////////////////////////
namespace internal {
/// Compute @e ax + @e y.
template<typename T>
inline void
axpy(const T& a, const T* xVec, const T* yVec, T* resultVec, SizeType size)
{
tbb::parallel_for(SizeRange(0, size), LinearOp<T>(a, xVec, yVec, resultVec));
}
/// Compute @e ax + @e y.
template<typename T>
inline void
axpy(const T& a, const Vector<T>& xVec, const Vector<T>& yVec, Vector<T>& result)
{
assert(xVec.size() == yVec.size());
assert(xVec.size() == result.size());
axpy(a, xVec.data(), yVec.data(), result.data(), xVec.size());
}
/// Compute @e r = @e b &minus; @e Ax.
template<typename MatrixOperator, typename VecValueType>
inline void
computeResidual(const MatrixOperator& A, const VecValueType* x,
const VecValueType* b, VecValueType* r)
{
// Compute r = A * x.
A.vectorMultiply(x, r);
// Compute r = b - r.
tbb::parallel_for(SizeRange(0, A.numRows()), LinearOp<VecValueType>(-1.0, r, b, r));
}
/// Compute @e r = @e b &minus; @e Ax.
template<typename MatrixOperator, typename T>
inline void
computeResidual(const MatrixOperator& A, const Vector<T>& x, const Vector<T>& b, Vector<T>& r)
{
assert(x.size() == b.size());
assert(x.size() == r.size());
assert(x.size() == A.numRows());
computeResidual(A, x.data(), b.data(), r.data());
}
} // namespace internal
////////////////////////////////////////
template<typename PositiveDefMatrix>
inline State
solve(
const PositiveDefMatrix& Amat,
const Vector<typename PositiveDefMatrix::ValueType>& bVec,
Vector<typename PositiveDefMatrix::ValueType>& xVec,
Preconditioner<typename PositiveDefMatrix::ValueType>& precond,
const State& termination)
{
util::NullInterrupter interrupter;
return solve(Amat, bVec, xVec, precond, interrupter, termination);
}
template<typename PositiveDefMatrix, typename Interrupter>
inline State
solve(
const PositiveDefMatrix& Amat,
const Vector<typename PositiveDefMatrix::ValueType>& bVec,
Vector<typename PositiveDefMatrix::ValueType>& xVec,
Preconditioner<typename PositiveDefMatrix::ValueType>& precond,
Interrupter& interrupter,
const State& termination)
{
using ValueType = typename PositiveDefMatrix::ValueType;
using VectorType = Vector<ValueType>;
State result;
result.success = false;
result.iterations = 0;
result.relativeError = 0.0;
result.absoluteError = 0.0;
const SizeType size = Amat.numRows();
if (size == 0) {
OPENVDB_LOG_WARN("pcg::solve(): matrix has dimension zero");
return result;
}
if (size != bVec.size()) {
OPENVDB_THROW(ArithmeticError, "A and b have incompatible sizes"
<< size << "x" << size << " vs. " << bVec.size() << ")");
}
if (size != xVec.size()) {
OPENVDB_THROW(ArithmeticError, "A and x have incompatible sizes"
<< size << "x" << size << " vs. " << xVec.size() << ")");
}
// Temp vectors
VectorType zVec(size); // transformed residual (M^-1 r)
VectorType pVec(size); // search direction
VectorType qVec(size); // A * p
// Compute norm of B (the source)
const ValueType tmp = bVec.infNorm();
const ValueType infNormOfB = isZero(tmp) ? ValueType(1) : tmp;
// Compute rVec: residual = b - Ax.
VectorType rVec(size); // vector of residuals
internal::computeResidual(Amat, xVec, bVec, rVec);
assert(rVec.isFinite());
// Normalize the residual norm with the source norm and look for early out.
result.absoluteError = static_cast<double>(rVec.infNorm());
result.relativeError = static_cast<double>(result.absoluteError / infNormOfB);
if (result.relativeError <= termination.relativeError) {
result.success = true;
return result;
}
// Iterations of the CG solve
ValueType rDotZPrev(1); // inner product of <z,r>
// Keep track of the minimum error to monitor convergence.
ValueType minL2Error = std::numeric_limits<ValueType>::max();
ValueType l2Error;
int iteration = 0;
for ( ; iteration < termination.iterations; ++iteration) {
if (interrupter.wasInterrupted()) {
OPENVDB_THROW(RuntimeError, "conjugate gradient solver was interrupted");
}
OPENVDB_LOG_DEBUG_RUNTIME("pcg::solve() " << result);
result.iterations = iteration + 1;
// Apply preconditioner to residual
// z_{k} = M^-1 r_{k}
precond.apply(rVec, zVec);
// <r,z>
const ValueType rDotZ = rVec.dot(zVec);
assert(std::isfinite(rDotZ));
if (0 == iteration) {
// Initialize
pVec = zVec;
} else {
const ValueType beta = rDotZ / rDotZPrev;
// p = beta * p + z
internal::axpy(beta, pVec, zVec, /*result */pVec);
}
// q_{k} = A p_{k}
Amat.vectorMultiply(pVec, qVec);
// alpha = <r_{k-1}, z_{k-1}> / <p_{k},q_{k}>
const ValueType pAp = pVec.dot(qVec);
assert(std::isfinite(pAp));
const ValueType alpha = rDotZ / pAp;
rDotZPrev = rDotZ;
// x_{k} = x_{k-1} + alpha * p_{k}
internal::axpy(alpha, pVec, xVec, /*result=*/xVec);
// r_{k} = r_{k-1} - alpha_{k-1} A p_{k}
internal::axpy(-alpha, qVec, rVec, /*result=*/rVec);
// update tolerances
l2Error = rVec.l2Norm();
minL2Error = Min(l2Error, minL2Error);
result.absoluteError = static_cast<double>(rVec.infNorm());
result.relativeError = static_cast<double>(result.absoluteError / infNormOfB);
if (l2Error > 2 * minL2Error) {
// The solution started to diverge.
result.success = false;
break;
}
if (!std::isfinite(result.absoluteError)) {
// Total divergence of solution
result.success = false;
break;
}
if (result.absoluteError <= termination.absoluteError) {
// Convergence
result.success = true;
break;
}
if (result.relativeError <= termination.relativeError) {
// Convergence
result.success = true;
break;
}
}
OPENVDB_LOG_DEBUG_RUNTIME("pcg::solve() " << result);
return result;
}
} // namespace pcg
} // namespace math
} // namespace OPENVDB_VERSION_NAME
} // namespace openvdb
#endif // OPENVDB_MATH_CONJGRADIENT_HAS_BEEN_INCLUDED
value
ValueT value
Definition GridBuilder.h:1290
Math.h
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
NullInterrupter.h
openvdb::v10_0::ArithmeticError
Definition Exceptions.h:56
openvdb::v10_0::RuntimeError
Definition Exceptions.h:63
openvdb::v10_0::math::pcg::IncompleteCholeskyPreconditioner
Preconditioner using incomplete Cholesky factorization.
Definition ConjGradient.h:1347
openvdb::v10_0::math::pcg::IncompleteCholeskyPreconditioner::TriangleRowEditor
typename TriangularMatrix::RowEditor TriangleRowEditor
Definition ConjGradient.h:1359
openvdb::v10_0::math::pcg::IncompleteCholeskyPreconditioner::Ptr
SharedPtr< IncompleteCholeskyPreconditioner > Ptr
Definition ConjGradient.h:1356
openvdb::v10_0::math::pcg::IncompleteCholeskyPreconditioner::IncompleteCholeskyPreconditioner
IncompleteCholeskyPreconditioner(const MatrixType &matrix)
Definition ConjGradient.h:1361
openvdb::v10_0::math::pcg::IncompleteCholeskyPreconditioner::TriangleConstRow
typename TriangularMatrix::ConstRow TriangleConstRow
Definition ConjGradient.h:1358
openvdb::v10_0::math::pcg::IncompleteCholeskyPreconditioner::BaseType
Preconditioner< ValueType > BaseType
Definition ConjGradient.h:1354
openvdb::v10_0::math::pcg::IncompleteCholeskyPreconditioner::TriangularMatrix
SparseStencilMatrix< ValueType, 4 > TriangularMatrix
Definition ConjGradient.h:1357
openvdb::v10_0::math::pcg::IncompleteCholeskyPreconditioner::apply
void apply(const Vector< ValueType > &rVec, Vector< ValueType > &zVec) override
Apply this preconditioner to a residue vector: z = M−1r.
Definition ConjGradient.h:1458
openvdb::v10_0::math::pcg::IncompleteCholeskyPreconditioner::~IncompleteCholeskyPreconditioner
~IncompleteCholeskyPreconditioner() override=default
openvdb::v10_0::math::pcg::IncompleteCholeskyPreconditioner::VectorType
Vector< ValueType > VectorType
Definition ConjGradient.h:1355
openvdb::v10_0::math::pcg::IncompleteCholeskyPreconditioner::ValueType
typename MatrixType::ValueType ValueType
Definition ConjGradient.h:1353
openvdb::v10_0::math::pcg::IncompleteCholeskyPreconditioner::upperMatrix
const TriangularMatrix & upperMatrix() const
Definition ConjGradient.h:1507
openvdb::v10_0::math::pcg::IncompleteCholeskyPreconditioner::isValid
bool isValid() const override
Definition ConjGradient.h:1456
openvdb::v10_0::math::pcg::IncompleteCholeskyPreconditioner::lowerMatrix
const TriangularMatrix & lowerMatrix() const
Definition ConjGradient.h:1506
openvdb::v10_0::math::pcg::JacobiPreconditioner
Diagonal preconditioner.
Definition ConjGradient.h:1281
openvdb::v10_0::math::pcg::JacobiPreconditioner::~JacobiPreconditioner
~JacobiPreconditioner() override=default
openvdb::v10_0::math::pcg::JacobiPreconditioner::Ptr
SharedPtr< JacobiPreconditioner > Ptr
Definition ConjGradient.h:1290
openvdb::v10_0::math::pcg::JacobiPreconditioner::apply
void apply(const Vector< ValueType > &r, Vector< ValueType > &z) override
Apply this preconditioner to a residue vector: z = M−1r.
Definition ConjGradient.h:1300
openvdb::v10_0::math::pcg::JacobiPreconditioner::BaseType
Preconditioner< ValueType > BaseType
Definition ConjGradient.h:1288
openvdb::v10_0::math::pcg::JacobiPreconditioner::VectorType
Vector< ValueType > VectorType
Definition ConjGradient.h:1289
openvdb::v10_0::math::pcg::JacobiPreconditioner::ValueType
typename MatrixType::ValueType ValueType
Definition ConjGradient.h:1287
openvdb::v10_0::math::pcg::JacobiPreconditioner::JacobiPreconditioner
JacobiPreconditioner(const MatrixType &A)
Definition ConjGradient.h:1292
openvdb::v10_0::math::pcg::JacobiPreconditioner::isFinite
bool isFinite() const
Return true if all values along the diagonal are finite.
Definition ConjGradient.h:1311
openvdb::v10_0::math::pcg::Preconditioner
Base class for conjugate gradient preconditioners.
Definition ConjGradient.h:465
openvdb::v10_0::math::pcg::Preconditioner::Preconditioner
Preconditioner(const SparseStencilMatrix< T, STENCIL_SIZE > &)
Definition ConjGradient.h:470
openvdb::v10_0::math::pcg::Preconditioner::isValid
virtual bool isValid() const
Definition ConjGradient.h:473
openvdb::v10_0::math::pcg::Preconditioner::apply
virtual void apply(const Vector< T > &r, Vector< T > &z)=0
Apply this preconditioner to a residue vector: z = M−1r.
openvdb::v10_0::math::pcg::Preconditioner::~Preconditioner
virtual ~Preconditioner()=default
openvdb::v10_0::math::pcg::Preconditioner::Ptr
SharedPtr< Preconditioner > Ptr
Definition ConjGradient.h:468
openvdb::v10_0::math::pcg::Preconditioner::ValueType
T ValueType
Definition ConjGradient.h:467
openvdb::v10_0::math::pcg::SparseStencilMatrix::ConstRow
Read-only accessor to a row of this matrix.
Definition ConjGradient.h:412
openvdb::v10_0::math::pcg::SparseStencilMatrix::ConstRow::ConstRow
ConstRow(const ValueType *valueHead, const SizeType *columnHead, const SizeType &rowSize)
Definition ConjGradient.h:1197
openvdb::v10_0::math::pcg::SparseStencilMatrix::ConstValueIter
Iterator over the stored values in a row of this matrix.
Definition ConjGradient.h:384
openvdb::v10_0::math::pcg::SparseStencilMatrix::ConstValueIter::operator*
const ValueType & operator*() const
Definition ConjGradient.h:386
openvdb::v10_0::math::pcg::SparseStencilMatrix::ConstValueIter::operator++
ConstValueIter & operator++()
Definition ConjGradient.h:395
openvdb::v10_0::math::pcg::SparseStencilMatrix::ConstValueIter::column
SizeType column() const
Definition ConjGradient.h:392
openvdb::v10_0::math::pcg::SparseStencilMatrix::ConstValueIter::reset
void reset()
Definition ConjGradient.h:398
openvdb::v10_0::math::pcg::SparseStencilMatrix::ConstValueIter::SparseStencilMatrix
friend class SparseStencilMatrix
Definition ConjGradient.h:401
openvdb::v10_0::math::pcg::SparseStencilMatrix::ConstValueIter::increment
void increment()
Definition ConjGradient.h:394
openvdb::v10_0::math::pcg::SparseStencilMatrix::RowEditor
Read/write accessor to a row of this matrix.
Definition ConjGradient.h:420
openvdb::v10_0::math::pcg::SparseStencilMatrix< ValueType, 4 >::RowEditor::setValue
SizeType setValue(SizeType column, const ValueType &value)
Definition ConjGradient.h:1225
openvdb::v10_0::math::pcg::SparseStencilMatrix::RowEditor::operator*=
RowEditor & operator*=(const Scalar &s)
Definition ConjGradient.h:435
openvdb::v10_0::math::pcg::SparseStencilMatrix< ValueType, 4 >::RowEditor::scale
void scale(const Scalar &)
Definition ConjGradient.h:1267
openvdb::v10_0::math::pcg::SparseStencilMatrix< ValueType, 4 >::RowEditor::clear
void clear()
Definition ConjGradient.h:1216
openvdb::v10_0::math::pcg::SparseStencilMatrix::RowEditor::RowEditor
RowEditor(ValueType *valueHead, SizeType *columnHead, SizeType &rowSize, SizeType colSize)
Definition ConjGradient.h:1207
openvdb::v10_0::math::pcg::SparseStencilMatrix
Sparse, square matrix representing a 3D stencil operator of size STENCIL_SIZE.
Definition ConjGradient.h:238
openvdb::v10_0::math::pcg::SparseStencilMatrix::size
SizeType size() const
Definition ConjGradient.h:259
openvdb::v10_0::math::pcg::SparseStencilMatrix::getValue
const ValueType & getValue(SizeType row, SizeType col) const
Return the value at the given coordinates.
Definition ConjGradient.h:888
openvdb::v10_0::math::pcg::SparseStencilMatrix::SparseStencilMatrix
SparseStencilMatrix(const SparseStencilMatrix &)
Deep copy the given matrix.
Definition ConjGradient.h:859
openvdb::v10_0::math::pcg::SparseStencilMatrix::setValue
void setValue(SizeType row, SizeType col, const ValueType &)
Set the value at the given coordinates.
Definition ConjGradient.h:878
openvdb::v10_0::math::pcg::SparseStencilMatrix::getConstRow
ConstRow getConstRow(SizeType row) const
Return a read-only view onto the given row of this matrix.
Definition ConjGradient.h:1079
openvdb::v10_0::math::pcg::SparseStencilMatrix::SparseStencilMatrix
SparseStencilMatrix(SizeType n)
Construct an n x n matrix with at most STENCIL_SIZE nonzero elements per row.
Definition ConjGradient.h:823
openvdb::v10_0::math::pcg::SparseStencilMatrix::operator()
const ValueType & operator()(SizeType row, SizeType col) const
Definition ConjGradient.h:897
openvdb::v10_0::math::pcg::SparseStencilMatrix::getRowEditor
RowEditor getRowEditor(SizeType row)
Return a read/write view onto the given row of this matrix.
Definition ConjGradient.h:1069
openvdb::v10_0::math::pcg::SparseStencilMatrix::ValueType
ValueType_ ValueType
Definition ConjGradient.h:240
openvdb::v10_0::math::pcg::SparseStencilMatrix::vectorMultiply
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
openvdb::v10_0::math::pcg::SparseStencilMatrix::operator*=
SparseStencilMatrix & operator*=(const Scalar &s)
Definition ConjGradient.h:285
openvdb::v10_0::math::pcg::SparseStencilMatrix< ValueType, 4 >::sZeroValue
static const ValueType sZeroValue
Definition ConjGradient.h:248
openvdb::v10_0::math::pcg::SparseStencilMatrix::scale
void scale(const Scalar &s)
Multiply all elements in the matrix by s;.
Definition ConjGradient.h:925
openvdb::v10_0::math::pcg::SparseStencilMatrix::VectorType
Vector< ValueType > VectorType
Definition ConjGradient.h:241
openvdb::v10_0::math::pcg::SparseStencilMatrix::vectorMultiply
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
openvdb::v10_0::math::pcg::SparseStencilMatrix::Ptr
SharedPtr< SparseStencilMatrix > Ptr
Definition ConjGradient.h:242
openvdb::v10_0::math::pcg::SparseStencilMatrix::eq
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
openvdb::v10_0::math::pcg::SparseStencilMatrix::str
std::string str() const
Return a string representation of this matrix.
Definition ConjGradient.h:1057
openvdb::v10_0::math::pcg::SparseStencilMatrix::numRows
SizeType numRows() const
Return the number of rows in this matrix.
Definition ConjGradient.h:258
openvdb::v10_0::math::pcg::SparseStencilMatrix::isFinite
bool isFinite() const
Return true if every element of this matrix has a finite value.
Definition ConjGradient.h:1046
openvdb::v10_0::math::pcg::Vector
Lightweight, variable-length vector.
Definition ConjGradient.h:139
openvdb::v10_0::math::pcg::Vector::size
SizeType size() const
Return the number of elements in this vector.
Definition ConjGradient.h:159
openvdb::v10_0::math::pcg::Vector::swap
void swap(Vector &other)
Swap internal storage with another vector, which need not be the same size.
Definition ConjGradient.h:168
openvdb::v10_0::math::pcg::Vector::l2Norm
ValueType l2Norm() const
Return the L2 norm of this vector.
Definition ConjGradient.h:185
openvdb::v10_0::math::pcg::Vector::dot
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
openvdb::v10_0::math::pcg::Vector::Ptr
SharedPtr< Vector > Ptr
Definition ConjGradient.h:142
openvdb::v10_0::math::pcg::Vector::data
const T * data() const
Definition ConjGradient.h:210
openvdb::v10_0::math::pcg::Vector::Vector
Vector()
Construct an empty vector.
Definition ConjGradient.h:145
openvdb::v10_0::math::pcg::Vector::eq
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
openvdb::v10_0::math::pcg::Vector::operator[]
const T & operator[](SizeType i) const
Definition ConjGradient.h:204
openvdb::v10_0::math::pcg::Vector::resize
void resize(SizeType n)
Reset this vector to have n elements, with uninitialized values.
Definition ConjGradient.h:588
openvdb::v10_0::math::pcg::Vector::at
const T & at(SizeType i) const
Definition ConjGradient.h:202
openvdb::v10_0::math::pcg::Vector::empty
bool empty() const
Return true if this vector has no elements.
Definition ConjGradient.h:161
openvdb::v10_0::math::pcg::Vector::at
T & at(SizeType i)
Return the value of this vector's ith element.
Definition ConjGradient.h:201
openvdb::v10_0::math::pcg::Vector::Vector
Vector(SizeType n, const ValueType &val)
Construct a vector of n elements and initialize each element to the given value.
Definition ConjGradient.h:149
openvdb::v10_0::math::pcg::Vector::Vector
Vector(const Vector &)
Deep copy the given vector.
Definition ConjGradient.h:559
openvdb::v10_0::math::pcg::Vector::scale
void scale(const Scalar &s)
Multiply each element of this vector by s.
Definition ConjGradient.h:624
openvdb::v10_0::math::pcg::Vector::constData
const T * constData() const
Definition ConjGradient.h:211
openvdb::v10_0::math::pcg::Vector::~Vector
~Vector()
Definition ConjGradient.h:151
openvdb::v10_0::math::pcg::Vector::operator*=
Vector & operator*=(const Scalar &s)
Definition ConjGradient.h:176
openvdb::v10_0::math::pcg::Vector::infNorm
ValueType infNorm() const
Return the infinity norm of this vector.
Definition ConjGradient.h:723
openvdb::v10_0::math::pcg::Vector< float >::fill
void fill(const ValueType &value)
openvdb::v10_0::math::pcg::Vector::operator=
Vector & operator=(const Vector &)
Deep copy the given vector.
Definition ConjGradient.h:568
openvdb::v10_0::math::pcg::Vector::Vector
Vector(SizeType n)
Construct a vector of n elements, with uninitialized values.
Definition ConjGradient.h:147
openvdb::v10_0::math::pcg::Vector::operator[]
T & operator[](SizeType i)
Definition ConjGradient.h:203
openvdb::v10_0::math::pcg::Vector::data
T * data()
Return a pointer to this vector's elements.
Definition ConjGradient.h:209
openvdb::v10_0::math::pcg::Vector::str
std::string str() const
Return a string representation of this vector.
Definition ConjGradient.h:800
openvdb::v10_0::math::pcg::Vector::ValueType
T ValueType
Definition ConjGradient.h:141
openvdb::v10_0::math::pcg::Vector::isFinite
bool isFinite() const
Return true if every element of this vector has a finite value.
Definition ConjGradient.h:753
logging.h
OPENVDB_LOG_WARN
#define OPENVDB_LOG_WARN(message)
Log a warning message of the form 'someVar << "some text" << ...'.
Definition logging.h:256
OPENVDB_LOG_DEBUG_RUNTIME
#define OPENVDB_LOG_DEBUG_RUNTIME(message)
Log a debugging message in both debug and optimized builds.
Definition logging.h:270
openvdb::v10_0::math::pcg::internal
Definition ConjGradient.h:486
openvdb::v10_0::math::pcg::internal::computeResidual
void computeResidual(const MatrixOperator &A, const VecValueType *x, const VecValueType *b, VecValueType *r)
Compute r = b − Ax.
Definition ConjGradient.h:1583
openvdb::v10_0::math::pcg::internal::axpy
void axpy(const T &a, const T *xVec, const T *yVec, T *resultVec, SizeType size)
Compute ax + y.
Definition ConjGradient.h:1564
openvdb::v10_0::math::pcg
Definition ConjGradient.h:31
openvdb::v10_0::math::pcg::VectorD
Vector< double > VectorD
Definition ConjGradient.h:228
openvdb::v10_0::math::pcg::SizeType
Index32 SizeType
Definition ConjGradient.h:33
openvdb::v10_0::math::pcg::operator<<
std::ostream & operator<<(std::ostream &os, const State &state)
Definition ConjGradient.h:545
openvdb::v10_0::math::pcg::VectorS
Vector< float > VectorS
Definition ConjGradient.h:227
openvdb::v10_0::math::pcg::solve
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
openvdb::v10_0::math::pcg::SizeRange
tbb::blocked_range< SizeType > SizeRange
Definition ConjGradient.h:35
openvdb::v10_0::math::pcg::terminationDefaults
State terminationDefaults()
Return default termination conditions for a conjugate gradient solver.
Definition ConjGradient.h:57
openvdb::v10_0::math::isApproxZero
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
openvdb::v10_0::math::Sqrt
float Sqrt(float x)
Return the square root of a floating-point value.
Definition Math.h:761
openvdb::v10_0::math::isApproxEqual
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
openvdb::v10_0::math::Max
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition Math.h:595
openvdb::v10_0::math::Min
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition Math.h:656
openvdb::v10_0::math::Abs
Coord Abs(const Coord &xyz)
Definition Coord.h:515
openvdb::v10_0::math::isExactlyEqual
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition Math.h:443
openvdb::v10_0::math::isZero
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition Math.h:337
openvdb::v10_0::Index32
uint32_t Index32
Definition Types.h:52
openvdb::v10_0::zeroVal
T zeroVal()
Return the value of type T that corresponds to zero.
Definition Math.h:70
openvdb::v10_0::SharedPtr
std::shared_ptr< T > SharedPtr
Definition Types.h:114
openvdb
Definition Exceptions.h:13
Exceptions.h
OPENVDB_THROW
#define OPENVDB_THROW(exception, message)
Definition Exceptions.h:74
Types.h
openvdb::v10_0::math::Tolerance::value
static T value()
Definition Math.h:148
openvdb::v10_0::math::pcg::SparseStencilMatrix::IsFiniteOp
Definition ConjGradient.h:1024
openvdb::v10_0::math::pcg::SparseStencilMatrix::IsFiniteOp::operator()
bool operator()(const SizeRange &range, bool finite) const
Definition ConjGradient.h:1027
openvdb::v10_0::math::pcg::SparseStencilMatrix::IsFiniteOp::mat
const SparseStencilMatrix * mat
Definition ConjGradient.h:1040
openvdb::v10_0::math::pcg::SparseStencilMatrix::IsFiniteOp::IsFiniteOp
IsFiniteOp(const SparseStencilMatrix &m)
Definition ConjGradient.h:1025
openvdb::v10_0::math::pcg::SparseStencilMatrix::MatrixCopyOp
Definition ConjGradient.h:837
openvdb::v10_0::math::pcg::SparseStencilMatrix::MatrixCopyOp::MatrixCopyOp
MatrixCopyOp(const SparseStencilMatrix &from_, SparseStencilMatrix &to_)
Definition ConjGradient.h:838
openvdb::v10_0::math::pcg::SparseStencilMatrix::MatrixCopyOp::to
SparseStencilMatrix * to
Definition ConjGradient.h:853
openvdb::v10_0::math::pcg::SparseStencilMatrix::MatrixCopyOp::operator()
void operator()(const SizeRange &range) const
Definition ConjGradient.h:841
openvdb::v10_0::math::pcg::SparseStencilMatrix::MatrixCopyOp::from
const SparseStencilMatrix * from
Definition ConjGradient.h:853
openvdb::v10_0::math::pcg::State
Information about the state of a conjugate gradient solution.
Definition ConjGradient.h:46
openvdb::v10_0::math::pcg::State::absoluteError
double absoluteError
Definition ConjGradient.h:50
openvdb::v10_0::math::pcg::State::iterations
int iterations
Definition ConjGradient.h:48
openvdb::v10_0::math::pcg::State::success
bool success
Definition ConjGradient.h:47
openvdb::v10_0::math::pcg::State::relativeError
double relativeError
Definition ConjGradient.h:49
openvdb::v10_0::math::pcg::Vector::DeterministicDotProductOp
Definition ConjGradient.h:632
openvdb::v10_0::math::pcg::Vector::DeterministicDotProductOp::reducetmp
T * reducetmp
Definition ConjGradient.h:659
openvdb::v10_0::math::pcg::Vector::DeterministicDotProductOp::operator()
void operator()(const SizeRange &range) const
Definition ConjGradient.h:637
openvdb::v10_0::math::pcg::Vector::DeterministicDotProductOp::binCount
const SizeType binCount
Definition ConjGradient.h:657
openvdb::v10_0::math::pcg::Vector::DeterministicDotProductOp::a
const T * a
Definition ConjGradient.h:655
openvdb::v10_0::math::pcg::Vector::DeterministicDotProductOp::arraySize
const SizeType arraySize
Definition ConjGradient.h:658
openvdb::v10_0::math::pcg::Vector::DeterministicDotProductOp::DeterministicDotProductOp
DeterministicDotProductOp(const T *a_, const T *b_, const SizeType binCount_, const SizeType arraySize_, T *reducetmp_)
Definition ConjGradient.h:633
openvdb::v10_0::math::pcg::Vector::DeterministicDotProductOp::b
const T * b
Definition ConjGradient.h:656
openvdb::v10_0::math::pcg::Vector::InfNormOp
Definition ConjGradient.h:706
openvdb::v10_0::math::pcg::Vector::InfNormOp::data
const T * data
Definition ConjGradient.h:717
openvdb::v10_0::math::pcg::Vector::InfNormOp::InfNormOp
InfNormOp(const T *data_)
Definition ConjGradient.h:707
openvdb::v10_0::math::pcg::Vector::InfNormOp::operator()
T operator()(const SizeRange &range, T maxValue) const
Definition ConjGradient.h:709
openvdb::v10_0::math::pcg::Vector::IsFiniteOp
Definition ConjGradient.h:734
openvdb::v10_0::math::pcg::Vector::IsFiniteOp::data
const T * data
Definition ConjGradient.h:747
openvdb::v10_0::math::pcg::Vector::IsFiniteOp::operator()
bool operator()(const SizeRange &range, bool finite) const
Definition ConjGradient.h:737
openvdb::v10_0::math::pcg::Vector::IsFiniteOp::IsFiniteOp
IsFiniteOp(const T *data_)
Definition ConjGradient.h:735
openvdb::v10_0::math::pcg::internal::CopyOp
Definition ConjGradient.h:491
openvdb::v10_0::math::pcg::internal::CopyOp::to
T * to
Definition ConjGradient.h:499
openvdb::v10_0::math::pcg::internal::CopyOp::CopyOp
CopyOp(const T *from_, T *to_)
Definition ConjGradient.h:492
openvdb::v10_0::math::pcg::internal::CopyOp::operator()
void operator()(const SizeRange &range) const
Definition ConjGradient.h:494
openvdb::v10_0::math::pcg::internal::CopyOp::from
const T * from
Definition ConjGradient.h:498
openvdb::v10_0::math::pcg::internal::FillOp
Definition ConjGradient.h:506
openvdb::v10_0::math::pcg::internal::FillOp::val
const T val
Definition ConjGradient.h:514
openvdb::v10_0::math::pcg::internal::FillOp::operator()
void operator()(const SizeRange &range) const
Definition ConjGradient.h:509
openvdb::v10_0::math::pcg::internal::FillOp::FillOp
FillOp(T *data_, const T &val_)
Definition ConjGradient.h:507
openvdb::v10_0::math::pcg::internal::FillOp::data
T * data
Definition ConjGradient.h:513
openvdb::v10_0::math::pcg::internal::LinearOp
Definition ConjGradient.h:521
openvdb::v10_0::math::pcg::internal::LinearOp::operator()
void operator()(const SizeRange &range) const
Definition ConjGradient.h:524
openvdb::v10_0::math::pcg::internal::LinearOp::LinearOp
LinearOp(const T &a_, const T *x_, const T *y_, T *out_)
Definition ConjGradient.h:522
openvdb::v10_0::math::pcg::internal::LinearOp::y
const T * y
Definition ConjGradient.h:534
openvdb::v10_0::math::pcg::internal::LinearOp::out
T * out
Definition ConjGradient.h:535
openvdb::v10_0::math::pcg::internal::LinearOp::x
const T * x
Definition ConjGradient.h:534
openvdb::v10_0::math::pcg::internal::LinearOp::a
const T a
Definition ConjGradient.h:534
openvdb::v10_0::util::NullInterrupter
Base class for interrupters.
Definition NullInterrupter.h:26
OPENVDB_VERSION_NAME
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition version.h.in:121
OPENVDB_USE_VERSION_NAMESPACE
#define OPENVDB_USE_VERSION_NAMESPACE
Definition version.h.in:212