9#ifndef OPENVDB_MATH_CONJGRADIENT_HAS_BEEN_INCLUDED
10#define OPENVDB_MATH_CONJGRADIENT_HAS_BEEN_INCLUDED
17#include <tbb/parallel_for.h>
18#include <tbb/parallel_reduce.h>
37template<
typename ValueType>
class Vector;
55template<
typename ValueType>
63 s.
absoluteError = std::numeric_limits<ValueType>::epsilon() * 100.0;
90template<
typename PositiveDefMatrix>
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>());
122template<
typename PositiveDefMatrix,
typename Interrupter>
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>());
151 ~Vector() { mSize = 0;
delete[] mData; mData =
nullptr; }
161 bool empty()
const {
return (mSize == 0); }
168 void swap(
Vector& other) { std::swap(mData, other.mData); std::swap(mSize, other.mSize); }
175 template<
typename Scalar>
void scale(
const Scalar& s);
192 template<
typename OtherValueType>
209 inline T*
data() {
return mData; }
210 inline const T*
data()
const {
return mData; }
216 template<
typename Scalar>
struct ScaleOp;
217 struct DeterministicDotProductOp;
219 template<
typename OtherValueType>
struct EqOp;
236template<
typename ValueType_, SizeType STENCIL_SIZE>
248 static inline constexpr ValueType sZeroValue = zeroVal<ValueType>();
283 template<
typename Scalar>
void scale(
const Scalar& s);
284 template<
typename Scalar>
291 template<
typename VecValueType>
298 template<
typename VecValueType>
303 template<
typename OtherValueType>
319 struct ConstRowData {
320 ConstRowData(
const ValueType* v,
const SizeType* c,
const SizeType& s):
321 mVals(v), mCols(c), mSize(s) {}
322 const ValueType* mVals;
const SizeType* mCols;
const SizeType& mSize;
326 template<
typename DataType_ = RowData>
330 using DataType = DataType_;
332 static SizeType capacity() {
return STENCIL_SIZE; }
334 RowBase(
const DataType& data): mData(data) {}
336 bool empty()
const {
return (mData.mSize == 0); }
337 const SizeType& size()
const {
return mData.mSize; }
339 const ValueType& getValue(SizeType columnIdx,
bool& active)
const;
340 const ValueType& getValue(SizeType columnIdx)
const;
343 ConstValueIter cbegin()
const;
347 template<
typename OtherDataType>
348 bool eq(
const RowBase<OtherDataType>& other,
349 ValueType eps = Tolerance<ValueType>::value())
const;
354 template<
typename VecValueType>
355 VecValueType dot(
const VecValueType* inVec, SizeType vecSize)
const;
358 template<
typename VecValueType>
359 VecValueType dot(
const Vector<VecValueType>& inVec)
const;
362 std::string str()
const;
365 friend class ConstValueIter;
367 const ValueType& value(SizeType i)
const {
return mData.mVals[i]; }
368 SizeType column(SizeType i)
const {
return mData.mCols[i]; }
374 SizeType find(SizeType columnIdx)
const;
379 using ConstRowBase = RowBase<ConstRowData>;
388 if (mData.mSize == 0)
return SparseStencilMatrix::sZeroValue;
389 return mData.mVals[mCursor];
396 operator bool()
const {
return (mCursor < mData.mSize); }
402 ConstValueIter(
const RowData& d): mData(d.mVals, d.mCols, d.mSize), mCursor(0) {}
405 const ConstRowData mData;
433 template<
typename Scalar>
void scale(
const Scalar&);
434 template<
typename Scalar>
445 template<
typename VecValueType>
struct VecMultOp;
446 template<
typename Scalar>
struct RowScaleOp;
450 template<
typename OtherValueType>
struct EqOp;
453 std::unique_ptr<ValueType[]> mValueArray;
454 std::unique_ptr<SizeType[]> mColumnIdxArray;
455 std::unique_ptr<SizeType[]> mRowSizeArray;
492 CopyOp(
const T* from_, T* to_): from(from_), to(to_) {}
495 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) to[n] = from[n];
507 FillOp(T* data_,
const T& val_): data(data_), val(val_) {}
510 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) data[n] = val;
522 LinearOp(
const T& a_,
const T* x_,
const T* y_, T* out_): a(a_), x(x_), y(y_), out(out_) {}
526 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) out[n] = x[n] + y[n];
528 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) out[n] = -x[n] + y[n];
530 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) out[n] = a * x[n] + y[n];
547 os << (state.
success ?
"succeeded with " :
"")
572 if (mSize != other.mSize) {
575 mData =
new T[mSize];
607template<
typename Scalar>
613 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n)
data[n] *= s;
622template<
typename Scalar>
635 a(a_), b(b_), binCount(binCount_), arraySize(arraySize_), reducetmp(reducetmp_) {}
639 const SizeType binSize = arraySize / binCount;
642 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
644 const SizeType end = (n == binCount-1) ? arraySize : begin + binSize;
648 for (
SizeType i = begin; i < end; ++i) { sum += a[i] * b[i]; }
666 assert(this->size() == other.
size());
668 const T* aData = this->data();
669 const T* bData = other.
data();
675 if (arraySize < 1024) {
679 for (
SizeType n = 0; n < arraySize; ++n) {
680 result += aData[n] * bData[n];
692 tbb::parallel_for(
SizeRange(0, binCount),
695 for (
SizeType n = 0; n < binCount; ++n) {
696 result += partialSums[n];
711 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
727 InfNormOp(this->data()), [](T max1, T max2) {
return Max(max1, max2); });
740 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
741 if (!std::isfinite(
data[n]))
return false;
756 bool finite = tbb::parallel_reduce(
SizeRange(0, this->size()),
true,
758 [](
bool finite1,
bool finite2) {
return (finite1 && finite2); });
764template<
typename OtherValueType>
767 EqOp(
const T* a_,
const OtherValueType* b_, T e): a(a_), b(b_), eps(e) {}
772 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
780 const OtherValueType*
b;
786template<
typename OtherValueType>
790 if (this->size() != other.
size())
return false;
791 bool equal = tbb::parallel_reduce(
SizeRange(0, this->size()),
true,
793 [](
bool eq1,
bool eq2) { return (eq1 && eq2); });
802 std::ostringstream ostr;
805 for (
SizeType n = 0, N = this->size(); n < N; ++n) {
806 ostr << sep << (*this)[n];
817template<
typename ValueType, SizeType STENCIL_SIZE>
821 , mValueArray(new
ValueType[mNumRows * STENCIL_SIZE])
822 , mColumnIdxArray(new
SizeType[mNumRows * STENCIL_SIZE])
823 , mRowSizeArray(new
SizeType[mNumRows])
826 tbb::parallel_for(
SizeRange(0, mNumRows),
831template<
typename ValueType, SizeType STENCIL_SIZE>
835 from(&from_), to(&to_) {}
839 const ValueType* fromVal = from->mValueArray.get();
840 const SizeType* fromCol = from->mColumnIdxArray.get();
841 ValueType* toVal = to->mValueArray.get();
842 SizeType* toCol = to->mColumnIdxArray.get();
843 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
844 toVal[n] = fromVal[n];
845 toCol[n] = fromCol[n];
853template<
typename ValueType, SizeType STENCIL_SIZE>
856 : mNumRows(other.mNumRows)
857 , mValueArray(new
ValueType[mNumRows * STENCIL_SIZE])
858 , mColumnIdxArray(new
SizeType[mNumRows * STENCIL_SIZE])
859 , mRowSizeArray(new
SizeType[mNumRows])
867 tbb::parallel_for(
SizeRange(0, mNumRows),
872template<
typename ValueType, SizeType STENCIL_SIZE>
877 assert(row < mNumRows);
878 this->getRowEditor(row).setValue(col, val);
882template<
typename ValueType, SizeType STENCIL_SIZE>
883inline const ValueType&
886 assert(row < mNumRows);
887 return this->getConstRow(row).getValue(col);
891template<
typename ValueType, SizeType STENCIL_SIZE>
892inline const ValueType&
895 return this->getValue(row,col);
899template<
typename ValueType, SizeType STENCIL_SIZE>
900template<
typename Scalar>
907 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
918template<
typename ValueType, SizeType STENCIL_SIZE>
919template<
typename Scalar>
928template<
typename ValueType, SizeType STENCIL_SIZE>
929template<
typename VecValueType>
933 mat(&m), in(i), out(o) {}
937 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
939 out[n] = row.dot(in, mat->numRows());
944 const VecValueType*
in;
949template<
typename ValueType, SizeType STENCIL_SIZE>
950template<
typename VecValueType>
955 if (inVec.
size() != mNumRows) {
957 << mNumRows <<
"x" << mNumRows <<
" vs. " << inVec.
size() <<
")");
959 if (resultVec.
size() != mNumRows) {
961 << mNumRows <<
"x" << mNumRows <<
" vs. " << resultVec.
size() <<
")");
964 vectorMultiply(inVec.
data(), resultVec.
data());
968template<
typename ValueType, SizeType STENCIL_SIZE>
969template<
typename VecValueType>
972 const VecValueType* inVec, VecValueType* resultVec)
const
975 tbb::parallel_for(
SizeRange(0, mNumRows),
980template<
typename ValueType, SizeType STENCIL_SIZE>
981template<
typename OtherValueType>
986 a(&a_), b(&b_), eps(e) {}
991 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
992 if (!a->getConstRow(n).eq(b->getConstRow(n), eps))
return false;
1004template<
typename ValueType, SizeType STENCIL_SIZE>
1005template<
typename OtherValueType>
1010 if (this->numRows() != other.
numRows())
return false;
1011 bool equal = tbb::parallel_reduce(
SizeRange(0, this->numRows()),
true,
1013 [](
bool eq1,
bool eq2) {
return (eq1 && eq2); });
1018template<
typename ValueType, SizeType STENCIL_SIZE>
1026 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1027 const ConstRow row = mat->getConstRow(n);
1029 if (!std::isfinite(*it))
return false;
1040template<
typename ValueType, SizeType STENCIL_SIZE>
1045 bool finite = tbb::parallel_reduce(
SizeRange(0, this->numRows()),
true,
1046 IsFiniteOp(*
this), [](
bool finite1,
bool finite2) {
return (finite1&&finite2); });
1051template<
typename ValueType, SizeType STENCIL_SIZE>
1055 std::ostringstream ostr;
1056 for (
SizeType n = 0, N = this->size(); n < N; ++n) {
1057 ostr << n <<
": " << this->getConstRow(n).str() <<
"\n";
1063template<
typename ValueType, SizeType STENCIL_SIZE>
1067 assert(i < mNumRows);
1068 const SizeType head = i * STENCIL_SIZE;
1069 return RowEditor(&mValueArray[head], &mColumnIdxArray[head], mRowSizeArray[i], mNumRows);
1073template<
typename ValueType, SizeType STENCIL_SIZE>
1077 assert(i < mNumRows);
1078 const SizeType head = i * STENCIL_SIZE;
1079 return ConstRow(&mValueArray[head], &mColumnIdxArray[head], mRowSizeArray[i]);
1083template<
typename ValueType, SizeType STENCIL_SIZE>
1084template<
typename DataType>
1088 if (this->empty())
return mData.mSize;
1092 const SizeType* colPtr = std::lower_bound(mData.mCols, mData.mCols + mData.mSize, columnIdx);
1094 return static_cast<SizeType>(colPtr - mData.mCols);
1098template<
typename ValueType, SizeType STENCIL_SIZE>
1099template<
typename DataType>
1100inline const ValueType&
1102 SizeType columnIdx,
bool& active)
const
1105 SizeType idx = this->find(columnIdx);
1106 if (idx < this->size() && this->column(idx) == columnIdx) {
1108 return this->value(idx);
1113template<
typename ValueType, SizeType STENCIL_SIZE>
1114template<
typename DataType>
1115inline const ValueType&
1118 SizeType idx = this->find(columnIdx);
1119 if (idx < this->size() && this->column(idx) == columnIdx) {
1120 return this->value(idx);
1126template<
typename ValueType, SizeType STENCIL_SIZE>
1127template<
typename DataType>
1128inline typename SparseStencilMatrix<ValueType, STENCIL_SIZE>::ConstValueIter
1129SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::cbegin()
const
1131 return ConstValueIter(mData);
1135template<
typename ValueType, SizeType STENCIL_SIZE>
1136template<
typename DataType>
1137template<
typename OtherDataType>
1140 const RowBase<OtherDataType>& other, ValueType eps)
const
1142 if (this->size() != other.size())
return false;
1143 for (ConstValueIter it = cbegin(), oit = other.cbegin(); it || oit; ++it, ++oit) {
1144 if (it.column() != oit.column())
return false;
1151template<
typename ValueType, SizeType STENCIL_SIZE>
1152template<
typename DataType>
1153template<
typename VecValueType>
1155SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::dot(
1156 const VecValueType* inVec,
SizeType vecSize)
const
1159 for (
SizeType idx = 0, N = std::min(vecSize, this->size()); idx < N; ++idx) {
1160 result +=
static_cast<VecValueType
>(this->value(idx) * inVec[this->column(idx)]);
1165template<
typename ValueType, SizeType STENCIL_SIZE>
1166template<
typename DataType>
1167template<
typename VecValueType>
1169SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::dot(
1170 const Vector<VecValueType>& inVec)
const
1172 return dot(inVec.data(), inVec.size());
1176template<
typename ValueType, SizeType STENCIL_SIZE>
1177template<
typename DataType>
1181 std::ostringstream ostr;
1183 for (
SizeType n = 0, N = this->size(); n < N; ++n) {
1184 ostr << sep <<
"(" << this->column(n) <<
", " << this->value(n) <<
")";
1191template<
typename ValueType, SizeType STENCIL_SIZE>
1195 ConstRowBase(ConstRowData(const_cast<
ValueType*>(valueHead),
1201template<
typename ValueType, SizeType STENCIL_SIZE>
1205 RowBase<>(RowData(valueHead, columnHead, rowSize)), mNumColumns(colSize)
1210template<
typename ValueType, SizeType STENCIL_SIZE>
1215 RowBase<>::mData.mSize = 0;
1219template<
typename ValueType, SizeType STENCIL_SIZE>
1224 assert(column < mNumColumns);
1226 RowData& data = RowBase<>::mData;
1230 SizeType offset = this->find(column);
1232 if (offset < data.mSize && data.mCols[offset] == column) {
1234 data.mVals[offset] = value;
1239 assert(data.mSize < this->capacity());
1241 if (offset >= data.mSize) {
1243 data.mVals[data.mSize] = value;
1244 data.mCols[data.mSize] = column;
1247 for (
SizeType i = data.mSize; i > offset; --i) {
1248 data.mVals[i] = data.mVals[i - 1];
1249 data.mCols[i] = data.mCols[i - 1];
1251 data.mVals[offset] = value;
1252 data.mCols[offset] = column;
1260template<
typename ValueType, SizeType STENCIL_SIZE>
1261template<
typename Scalar>
1265 for (
int idx = 0, N = this->
size(); idx < N; ++idx) {
1266 RowBase<>::mData.mVals[idx] *= s;
1275template<
typename MatrixType>
1291 tbb::parallel_for(
SizeRange(0, A.numRows()), InitOp(A, mDiag.data()));
1298 const SizeType size = mDiag.size();
1301 assert(r.
size() == size);
1303 tbb::parallel_for(
SizeRange(0, size), ApplyOp(mDiag.data(), r.
data(), z.
data()));
1313 InitOp(
const MatrixType& m, ValueType* v): mat(&m), vec(v) {}
1314 void operator()(
const SizeRange& range)
const {
1315 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1316 const ValueType val = mat->getValue(n, n);
1318 vec[n] =
static_cast<ValueType
>(1.0 / val);
1321 const MatrixType* mat; ValueType* vec;
1327 ApplyOp(
const ValueType* x_,
const ValueType* y_, ValueType* out_):
1328 x(x_), y(y_), out(out_) {}
1329 void operator()(
const SizeRange& range)
const {
1330 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) out[n] = x[n] * y[n];
1332 const ValueType *x, *y; ValueType* out;
1341template<
typename MatrixType>
1345 struct CopyToLowerOp;
1359 , mLowerTriangular(matrix.numRows())
1360 , mUpperTriangular(matrix.numRows())
1361 , mTempVec(matrix.numRows())
1364 const SizeType numRows = mLowerTriangular.numRows();
1367 tbb::parallel_for(
SizeRange(0, numRows), CopyToLowerOp(matrix, mLowerTriangular));
1387 mPassedCompatibilityCondition =
true;
1389 for (
SizeType k = 0; k < numRows; ++k) {
1392 ValueType diagonalValue = crow_k.getValue(k);
1395 if (diagonalValue < 1.e-5) {
1396 mPassedCompatibilityCondition =
false;
1400 diagonalValue =
Sqrt(diagonalValue);
1403 row_k.setValue(k, diagonalValue);
1406 typename MatrixType::ConstRow srcRow = matrix.getConstRow(k);
1407 typename MatrixType::ConstValueIter citer = srcRow.cbegin();
1408 for ( ; citer; ++citer) {
1410 if (ii < k+1)
continue;
1414 row_ii.setValue(k, *citer / diagonalValue);
1419 for ( ; citer; ++citer) {
1421 if (j < k+1)
continue;
1428 typename MatrixType::ConstRow mask = matrix.getConstRow(j);
1429 typename MatrixType::ConstValueIter maskIter = mask.cbegin();
1430 for ( ; maskIter; ++maskIter) {
1432 if (i < j)
continue;
1438 a_ij -= a_ik * a_jk;
1440 row_i.setValue(j, a_ij);
1446 tbb::parallel_for(
SizeRange(0, numRows),
1447 TransposeOp(matrix, mLowerTriangular, mUpperTriangular));
1452 bool isValid()
const override {
return mPassedCompatibilityCondition; }
1456 if (!mPassedCompatibilityCondition) {
1462 SizeType size = mLowerTriangular.numRows();
1464 zVec.
fill(zeroVal<ValueType>());
1467 if (size == 0)
return;
1469 assert(rVec.
size() == size);
1470 assert(zVec.
size() == size);
1473 mTempVec.fill(zeroVal<ValueType>());
1478 for (
SizeType i = 0; i < size; ++i) {
1479 typename TriangularMatrix::ConstRow row = mLowerTriangular.getConstRow(i);
1482 tmpData[i] = (rData[i] - dot) / diagonal;
1483 if (!std::isfinite(tmpData[i])) {
1490 for (
SizeType ii = 0; ii < size; ++ii) {
1492 typename TriangularMatrix::ConstRow row = mUpperTriangular.getConstRow(i);
1495 zData[i] = (tmpData[i] - dot) / diagonal;
1496 if (!std::isfinite(zData[i])) {
1507 struct CopyToLowerOp
1509 CopyToLowerOp(
const MatrixType& m, TriangularMatrix& l): mat(&m), lower(&l) {}
1510 void operator()(
const SizeRange& range)
const {
1511 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1512 typename TriangularMatrix::RowEditor outRow = lower->getRowEditor(n);
1514 typename MatrixType::ConstRow inRow = mat->getConstRow(n);
1515 for (
typename MatrixType::ConstValueIter it = inRow.cbegin(); it; ++it) {
1516 if (it.column() > n)
continue;
1517 outRow.setValue(it.column(), *it);
1521 const MatrixType* mat; TriangularMatrix* lower;
1527 TransposeOp(
const MatrixType& m,
const TriangularMatrix& l, TriangularMatrix& u):
1528 mat(&m), lower(&l), upper(&u) {}
1529 void operator()(
const SizeRange& range)
const {
1530 for (SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1531 typename TriangularMatrix::RowEditor outRow = upper->getRowEditor(n);
1534 typename MatrixType::ConstRow inRow = mat->getConstRow(n);
1535 for (
typename MatrixType::ConstValueIter it = inRow.cbegin(); it; ++it) {
1536 const SizeType column = it.column();
1537 if (column < n)
continue;
1538 outRow.setValue(column, lower->getValue(column, n));
1542 const MatrixType* mat;
const TriangularMatrix* lower; TriangularMatrix* upper;
1545 TriangularMatrix mLowerTriangular;
1546 TriangularMatrix mUpperTriangular;
1547 Vector<ValueType> mTempVec;
1548 bool mPassedCompatibilityCondition;
1570 assert(xVec.
size() == yVec.
size());
1571 assert(xVec.
size() == result.
size());
1577template<
typename MatrixOperator,
typename VecValueType>
1580 const VecValueType* b, VecValueType* r)
1583 A.vectorMultiply(x, r);
1589template<
typename MatrixOperator,
typename T>
1595 assert(x.
size() == A.numRows());
1606template<
typename PositiveDefMatrix>
1609 const PositiveDefMatrix& Amat,
1613 const State& termination)
1616 return solve(Amat, bVec, xVec, precond, interrupter, termination);
1620template<
typename PositiveDefMatrix,
typename Interrupter>
1623 const PositiveDefMatrix& Amat,
1627 Interrupter& interrupter,
1628 const State& termination)
1630 using ValueType =
typename PositiveDefMatrix::ValueType;
1646 <<
size <<
"x" <<
size <<
" vs. " << bVec.
size() <<
")");
1650 <<
size <<
"x" <<
size <<
" vs. " << xVec.
size() <<
")");
1667 assert(rVec.isFinite());
1682 ValueType minL2Error = std::numeric_limits<ValueType>::max();
1686 for ( ; iteration < termination.
iterations; ++iteration) {
1688 if (interrupter.wasInterrupted()) {
1698 precond.
apply(rVec, zVec);
1702 assert(std::isfinite(rDotZ));
1704 if (0 == iteration) {
1708 const ValueType beta = rDotZ / rDotZPrev;
1714 Amat.vectorMultiply(pVec, qVec);
1718 assert(std::isfinite(pAp));
1730 l2Error = rVec.l2Norm();
1731 minL2Error =
Min(l2Error, minL2Error);
1736 if (l2Error > 2 * minL2Error) {
OPENVDB_API std::ostream & operator<<(std::ostream &os, half h)
Output h to os, formatted as a float.
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
Definition Exceptions.h:56
Definition Exceptions.h:63
Preconditioner using incomplete Cholesky factorization.
Definition ConjGradient.h:1343
typename TriangularMatrix::RowEditor TriangleRowEditor
Definition ConjGradient.h:1355
SharedPtr< IncompleteCholeskyPreconditioner > Ptr
Definition ConjGradient.h:1352
IncompleteCholeskyPreconditioner(const MatrixType &matrix)
Definition ConjGradient.h:1357
typename TriangularMatrix::ConstRow TriangleConstRow
Definition ConjGradient.h:1354
void apply(const Vector< ValueType > &rVec, Vector< ValueType > &zVec) override
Apply this preconditioner to a residue vector: z = M−1r
Definition ConjGradient.h:1454
~IncompleteCholeskyPreconditioner() override=default
typename MatrixType::ValueType ValueType
Definition ConjGradient.h:1349
const TriangularMatrix & upperMatrix() const
Definition ConjGradient.h:1503
bool isValid() const override
Definition ConjGradient.h:1452
const TriangularMatrix & lowerMatrix() const
Definition ConjGradient.h:1502
Diagonal preconditioner.
Definition ConjGradient.h:1277
~JacobiPreconditioner() override=default
SharedPtr< JacobiPreconditioner > Ptr
Definition ConjGradient.h:1286
void apply(const Vector< ValueType > &r, Vector< ValueType > &z) override
Apply this preconditioner to a residue vector: z = M−1r
Definition ConjGradient.h:1296
typename MatrixType::ValueType ValueType
Definition ConjGradient.h:1283
JacobiPreconditioner(const MatrixType &A)
Definition ConjGradient.h:1288
bool isFinite() const
Return true if all values along the diagonal are finite.
Definition ConjGradient.h:1307
Base class for conjugate gradient preconditioners.
Definition ConjGradient.h:465
Preconditioner(const SparseStencilMatrix< T, STENCIL_SIZE > &)
Definition ConjGradient.h:470
virtual bool isValid() const
Definition ConjGradient.h:473
virtual void apply(const Vector< T > &r, Vector< T > &z)=0
Apply this preconditioner to a residue vector: z = M−1r
virtual ~Preconditioner()=default
SharedPtr< Preconditioner > Ptr
Definition ConjGradient.h:468
T ValueType
Definition ConjGradient.h:467
Read-only accessor to a row of this matrix.
Definition ConjGradient.h:412
ConstRow(const ValueType *valueHead, const SizeType *columnHead, const SizeType &rowSize)
Definition ConjGradient.h:1193
Iterator over the stored values in a row of this matrix.
Definition ConjGradient.h:384
const ValueType & operator*() const
Definition ConjGradient.h:386
ConstValueIter & operator++()
Definition ConjGradient.h:395
SizeType column() const
Definition ConjGradient.h:392
void reset()
Definition ConjGradient.h:398
void increment()
Definition ConjGradient.h:394
Read/write accessor to a row of this matrix.
Definition ConjGradient.h:420
SizeType setValue(SizeType column, const ValueType &value)
Set the value of the entry in the specified column.
Definition ConjGradient.h:1221
RowEditor & operator*=(const Scalar &s)
Definition ConjGradient.h:435
void scale(const Scalar &)
Scale all of the entries in this row.
Definition ConjGradient.h:1263
void clear()
Set the number of entries in this row to zero.
Definition ConjGradient.h:1212
RowEditor(ValueType *valueHead, SizeType *columnHead, SizeType &rowSize, SizeType colSize)
Definition ConjGradient.h:1203
Sparse, square matrix representing a 3D stencil operator of size STENCIL_SIZE.
Definition ConjGradient.h:238
SizeType size() const
Definition ConjGradient.h:259
const ValueType & getValue(SizeType row, SizeType col) const
Return the value at the given coordinates.
Definition ConjGradient.h:884
SparseStencilMatrix(const SparseStencilMatrix &)
Deep copy the given matrix.
Definition ConjGradient.h:855
void setValue(SizeType row, SizeType col, const ValueType &)
Set the value at the given coordinates.
Definition ConjGradient.h:874
ConstRow getConstRow(SizeType row) const
Return a read-only view onto the given row of this matrix.
Definition ConjGradient.h:1075
SparseStencilMatrix(SizeType n)
Construct an n x n matrix with at most STENCIL_SIZE nonzero elements per row.
Definition ConjGradient.h:819
const ValueType & operator()(SizeType row, SizeType col) const
Definition ConjGradient.h:893
RowEditor getRowEditor(SizeType row)
Return a read/write view onto the given row of this matrix.
Definition ConjGradient.h:1065
ValueType_ ValueType
Definition ConjGradient.h:240
void vectorMultiply(const Vector< VecValueType > &inVec, Vector< VecValueType > &resultVec) const
Multiply this matrix by inVec and return the result in resultVec.
Definition ConjGradient.h:952
SparseStencilMatrix & operator*=(const Scalar &s)
Definition ConjGradient.h:285
static constexpr ValueType sZeroValue
Definition ConjGradient.h:248
void scale(const Scalar &s)
Multiply all elements in the matrix by s;.
Definition ConjGradient.h:921
Vector< ValueType > VectorType
Definition ConjGradient.h:241
void vectorMultiply(const VecValueType *inVec, VecValueType *resultVec) const
Multiply this matrix by the vector represented by the array inVec and return the result in resultVec.
Definition ConjGradient.h:971
SharedPtr< SparseStencilMatrix > Ptr
Definition ConjGradient.h:242
bool eq(const SparseStencilMatrix< OtherValueType, STENCIL_SIZE > &other, ValueType eps=Tolerance< ValueType >::value()) const
Return true if this matrix is equivalent to the given matrix to within the specified tolerance.
Definition ConjGradient.h:1007
std::string str() const
Return a string representation of this matrix.
Definition ConjGradient.h:1053
SizeType numRows() const
Return the number of rows in this matrix.
Definition ConjGradient.h:258
bool isFinite() const
Return true if every element of this matrix has a finite value.
Definition ConjGradient.h:1042
Lightweight, variable-length vector.
Definition ConjGradient.h:139
SizeType size() const
Return the number of elements in this vector.
Definition ConjGradient.h:159
void swap(Vector &other)
Swap internal storage with another vector, which need not be the same size.
Definition ConjGradient.h:168
ValueType l2Norm() const
Return the L2 norm of this vector.
Definition ConjGradient.h:185
ValueType dot(const Vector &) const
Return the dot product of this vector with the given vector, which must be the same size.
Definition ConjGradient.h:664
SharedPtr< Vector > Ptr
Definition ConjGradient.h:142
const T * data() const
Definition ConjGradient.h:210
Vector()
Construct an empty vector.
Definition ConjGradient.h:145
bool eq(const Vector< OtherValueType > &other, ValueType eps=Tolerance< ValueType >::value()) const
Return true if this vector is equivalent to the given vector to within the specified tolerance.
Definition ConjGradient.h:788
const T & operator[](SizeType i) const
Definition ConjGradient.h:204
void resize(SizeType n)
Reset this vector to have n elements, with uninitialized values.
Definition ConjGradient.h:588
const T & at(SizeType i) const
Definition ConjGradient.h:202
bool empty() const
Return true if this vector has no elements.
Definition ConjGradient.h:161
T & at(SizeType i)
Return the value of this vector's ith element.
Definition ConjGradient.h:201
Vector(SizeType n, const ValueType &val)
Construct a vector of n elements and initialize each element to the given value.
Definition ConjGradient.h:149
Vector(const Vector &)
Deep copy the given vector.
Definition ConjGradient.h:559
void scale(const Scalar &s)
Multiply each element of this vector by s.
Definition ConjGradient.h:624
const T * constData() const
Definition ConjGradient.h:211
~Vector()
Definition ConjGradient.h:151
Vector & operator*=(const Scalar &s)
Definition ConjGradient.h:176
ValueType infNorm() const
Return the infinity norm of this vector.
Definition ConjGradient.h:723
void fill(const ValueType &value)
Set all elements of this vector to value.
Definition ConjGradient.h:600
Vector & operator=(const Vector &)
Deep copy the given vector.
Definition ConjGradient.h:568
Vector(SizeType n)
Construct a vector of n elements, with uninitialized values.
Definition ConjGradient.h:147
T & operator[](SizeType i)
Definition ConjGradient.h:203
T * data()
Return a pointer to this vector's elements.
Definition ConjGradient.h:209
std::string str() const
Return a string representation of this vector.
Definition ConjGradient.h:800
T ValueType
Definition ConjGradient.h:141
bool isFinite() const
Return true if every element of this vector has a finite value.
Definition ConjGradient.h:753
#define OPENVDB_LOG_WARN(message)
Log a warning message of the form 'someVar << "some text" << ...'.
Definition logging.h:256
#define OPENVDB_LOG_DEBUG_RUNTIME(message)
Log a debugging message in both debug and optimized builds.
Definition logging.h:270
void computeResidual(const MatrixOperator &A, const VecValueType *x, const VecValueType *b, VecValueType *r)
Compute r = b − Ax.
Definition ConjGradient.h:1579
void axpy(const T &a, const T *xVec, const T *yVec, T *resultVec, SizeType size)
Compute ax + y.
Definition ConjGradient.h:1560
Index32 SizeType
Definition ConjGradient.h:33
State solve(const PositiveDefMatrix &A, const Vector< typename PositiveDefMatrix::ValueType > &b, Vector< typename PositiveDefMatrix::ValueType > &x, Preconditioner< typename PositiveDefMatrix::ValueType > &preconditioner, const State &termination=terminationDefaults< typename PositiveDefMatrix::ValueType >())
Solve Ax = b via the preconditioned conjugate gradient method.
Definition ConjGradient.h:1608
tbb::blocked_range< SizeType > SizeRange
Definition ConjGradient.h:35
State terminationDefaults()
Return default termination conditions for a conjugate gradient solver.
Definition ConjGradient.h:57
bool isApproxZero(const Type &x)
Return true if x is equal to zero to within the default floating-point comparison tolerance.
Definition Math.h:349
float Sqrt(float x)
Return the square root of a floating-point value.
Definition Math.h:761
bool isApproxEqual(const Type &a, const Type &b, const Type &tolerance)
Return true if a is equal to b to within the given tolerance.
Definition Math.h:406
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition Math.h:595
MatType scale(const Vec3< typename MatType::value_type > &s)
Return a matrix that scales by s.
Definition Mat.h:615
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition Math.h:656
Coord Abs(const Coord &xyz)
Definition Coord.h:517
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition Math.h:443
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition Math.h:337
constexpr T zeroVal()
Return the value of type T that corresponds to zero.
Definition Math.h:70
uint32_t Index32
Definition Types.h:52
std::shared_ptr< T > SharedPtr
Definition Types.h:114
Definition Exceptions.h:13
#define OPENVDB_THROW(exception, message)
Definition Exceptions.h:74
Tolerance for floating-point comparison.
Definition Math.h:148
Definition ConjGradient.h:983
bool operator()(const SizeRange &range, bool equal) const
Definition ConjGradient.h:988
const ValueType eps
Definition ConjGradient.h:1000
const SparseStencilMatrix< OtherValueType, STENCIL_SIZE > * b
Definition ConjGradient.h:999
EqOp(const SparseStencilMatrix &a_, const SparseStencilMatrix< OtherValueType, STENCIL_SIZE > &b_, ValueType e)
Definition ConjGradient.h:984
const SparseStencilMatrix * a
Definition ConjGradient.h:998
Definition ConjGradient.h:1020
bool operator()(const SizeRange &range, bool finite) const
Definition ConjGradient.h:1023
const SparseStencilMatrix * mat
Definition ConjGradient.h:1036
IsFiniteOp(const SparseStencilMatrix &m)
Definition ConjGradient.h:1021
Definition ConjGradient.h:833
MatrixCopyOp(const SparseStencilMatrix &from_, SparseStencilMatrix &to_)
Definition ConjGradient.h:834
void operator()(const SizeRange &range) const
Definition ConjGradient.h:837
const SparseStencilMatrix * from
Definition ConjGradient.h:849
Definition ConjGradient.h:902
const Scalar s
Definition ConjGradient.h:914
void operator()(const SizeRange &range) const
Definition ConjGradient.h:905
RowScaleOp(SparseStencilMatrix &m, const Scalar &s_)
Definition ConjGradient.h:903
SparseStencilMatrix * mat
Definition ConjGradient.h:913
Definition ConjGradient.h:931
const VecValueType * in
Definition ConjGradient.h:944
const SparseStencilMatrix * mat
Definition ConjGradient.h:943
void operator()(const SizeRange &range) const
Definition ConjGradient.h:935
VecMultOp(const SparseStencilMatrix &m, const VecValueType *i, VecValueType *o)
Definition ConjGradient.h:932
VecValueType * out
Definition ConjGradient.h:945
Information about the state of a conjugate gradient solution.
Definition ConjGradient.h:46
double absoluteError
Definition ConjGradient.h:50
int iterations
Definition ConjGradient.h:48
bool success
Definition ConjGradient.h:47
double relativeError
Definition ConjGradient.h:49
Definition ConjGradient.h:632
T * reducetmp
Definition ConjGradient.h:659
void operator()(const SizeRange &range) const
Definition ConjGradient.h:637
const SizeType binCount
Definition ConjGradient.h:657
const T * a
Definition ConjGradient.h:655
const SizeType arraySize
Definition ConjGradient.h:658
DeterministicDotProductOp(const T *a_, const T *b_, const SizeType binCount_, const SizeType arraySize_, T *reducetmp_)
Definition ConjGradient.h:633
const T * b
Definition ConjGradient.h:656
Definition ConjGradient.h:766
const OtherValueType * b
Definition ConjGradient.h:780
bool operator()(const SizeRange &range, bool equal) const
Definition ConjGradient.h:769
EqOp(const T *a_, const OtherValueType *b_, T e)
Definition ConjGradient.h:767
const T * a
Definition ConjGradient.h:779
const T eps
Definition ConjGradient.h:781
Definition ConjGradient.h:706
const T * data
Definition ConjGradient.h:717
InfNormOp(const T *data_)
Definition ConjGradient.h:707
T operator()(const SizeRange &range, T maxValue) const
Definition ConjGradient.h:709
Definition ConjGradient.h:734
const T * data
Definition ConjGradient.h:747
bool operator()(const SizeRange &range, bool finite) const
Definition ConjGradient.h:737
IsFiniteOp(const T *data_)
Definition ConjGradient.h:735
Definition ConjGradient.h:609
const Scalar s
Definition ConjGradient.h:617
void operator()(const SizeRange &range) const
Definition ConjGradient.h:612
ScaleOp(T *data_, const Scalar &s_)
Definition ConjGradient.h:610
T * data
Definition ConjGradient.h:616
Definition ConjGradient.h:491
T * to
Definition ConjGradient.h:499
CopyOp(const T *from_, T *to_)
Definition ConjGradient.h:492
void operator()(const SizeRange &range) const
Definition ConjGradient.h:494
const T * from
Definition ConjGradient.h:498
Definition ConjGradient.h:506
const T val
Definition ConjGradient.h:514
void operator()(const SizeRange &range) const
Definition ConjGradient.h:509
FillOp(T *data_, const T &val_)
Definition ConjGradient.h:507
T * data
Definition ConjGradient.h:513
Definition ConjGradient.h:521
void operator()(const SizeRange &range) const
Definition ConjGradient.h:524
LinearOp(const T &a_, const T *x_, const T *y_, T *out_)
Definition ConjGradient.h:522
T * out
Definition ConjGradient.h:535
const T a
Definition ConjGradient.h:534
Base class for interrupters.
Definition NullInterrupter.h:26
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition version.h.in:121
#define OPENVDB_USE_VERSION_NAMESPACE
Definition version.h.in:212