36 #ifndef OPENVDB_MATH_CONJGRADIENT_HAS_BEEN_INCLUDED 37 #define OPENVDB_MATH_CONJGRADIENT_HAS_BEEN_INCLUDED 39 #include <openvdb/Exceptions.h> 40 #include <openvdb/Types.h> 41 #include <openvdb/util/logging.h> 42 #include <openvdb/util/NullInterrupter.h> 44 #include <tbb/parallel_for.h> 45 #include <tbb/parallel_reduce.h> 64 template<
typename ValueType>
class Vector;
82 template<
typename ValueType>
90 s.
absoluteError = std::numeric_limits<ValueType>::epsilon() * 100.0;
117 template<
typename PositiveDefMatrix>
120 const PositiveDefMatrix& A,
124 const State& termination = terminationDefaults<typename PositiveDefMatrix::ValueType>());
149 template<
typename PositiveDefMatrix,
typename Interrupter>
152 const PositiveDefMatrix& A,
156 Interrupter& interrupter,
157 const State& termination = terminationDefaults<typename PositiveDefMatrix::ValueType>());
178 ~Vector() { mSize = 0;
delete[] mData; mData =
nullptr; }
188 bool empty()
const {
return (mSize == 0); }
195 void swap(
Vector& other) { std::swap(mData, other.mData); std::swap(mSize, other.mSize); }
201 template<
typename Scalar>
void scale(
const Scalar& s);
219 template<
typename OtherValueType>
224 std::string str()
const;
227 inline T& at(
SizeType i) {
return mData[i]; }
235 inline T* data() {
return mData; }
237 inline const T*
data()
const {
return mData; }
243 template<
typename Scalar>
struct ScaleOp;
244 struct DeterministicDotProductOp;
246 template<
typename OtherValueType>
struct EqOp;
263 template<
typename ValueType_, SizeType STENCIL_SIZE>
271 class ConstValueIter;
284 SizeType numRows()
const {
return mNumRows; }
303 ConstRow getConstRow(
SizeType row)
const;
306 RowEditor getRowEditor(
SizeType row);
309 template<
typename Scalar>
void scale(
const Scalar& s);
311 template<
typename Scalar>
318 template<
typename VecValueType>
325 template<
typename VecValueType>
326 void vectorMultiply(
const VecValueType* inVec, VecValueType* resultVec)
const;
330 template<
typename OtherValueType>
338 std::string str()
const;
346 struct ConstRowData {
348 mVals(v), mCols(c), mSize(s) {}
353 template<
typename DataType_ = RowData>
357 using DataType = DataType_;
359 static SizeType capacity() {
return STENCIL_SIZE; }
361 RowBase(
const DataType& data): mData(data) {}
363 bool empty()
const {
return (mData.mSize == 0); }
364 const SizeType& size()
const {
return mData.mSize; }
370 ConstValueIter cbegin()
const;
374 template<
typename OtherDataType>
375 bool eq(
const RowBase<OtherDataType>& other,
381 template<
typename VecValueType>
382 VecValueType dot(
const VecValueType* inVec,
SizeType vecSize)
const;
385 template<
typename VecValueType>
389 std::string str()
const;
392 friend class ConstValueIter;
406 using ConstRowBase = RowBase<ConstRowData>;
415 if (mData.mSize == 0)
return SparseStencilMatrix::sZeroValue;
416 return mData.mVals[mCursor];
423 operator bool()
const {
return (mCursor < mData.mSize); }
429 ConstValueIter(
const RowData& d): mData(d.mVals, d.mCols, d.mSize), mCursor(0) {}
432 const ConstRowData mData;
459 template<
typename Scalar>
void scale(
const Scalar&);
461 template<
typename Scalar>
472 template<
typename VecValueType>
struct VecMultOp;
473 template<
typename Scalar>
struct RowScaleOp;
477 template<
typename OtherValueType>
struct EqOp;
480 std::unique_ptr<ValueType[]> mValueArray;
481 std::unique_ptr<SizeType[]> mColumnIdxArray;
482 std::unique_ptr<SizeType[]> mRowSizeArray;
519 CopyOp(
const T* from_, T* to_): from(from_), to(to_) {}
522 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) to[n] = from[n];
534 FillOp(T* data_,
const T& val_): data(data_), val(val_) {}
537 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) data[n] = val;
549 LinearOp(
const T& a_,
const T* x_,
const T* y_, T* out_): a(a_), x(x_), y(y_), out(out_) {}
553 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) out[n] = x[n] + y[n];
555 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) out[n] = -x[n] + y[n];
557 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) out[n] = a * x[n] + y[n];
574 os << (state.
success ?
"succeeded with " :
"")
599 if (mSize != other.mSize) {
602 mData =
new T[mSize];
618 if (mData)
delete[] mData;
634 template<
typename Scalar>
637 ScaleOp(T* data_,
const Scalar& s_): data(data_), s(s_) {}
639 void operator()(
const SizeRange& range)
const {
640 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) data[n] *= s;
649 template<
typename Scalar>
653 tbb::parallel_for(
SizeRange(0, mSize), ScaleOp<Scalar>(mData, s));
662 a(a_), b(b_), binCount(binCount_), arraySize(arraySize_), reducetmp(reducetmp_) {}
666 const SizeType binSize = arraySize / binCount;
669 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
671 const SizeType end = (n == binCount-1) ? arraySize : begin + binSize;
674 T sum = zeroVal<T>();
675 for (
SizeType i = begin; i < end; ++i) { sum += a[i] * b[i]; }
693 assert(this->size() == other.
size());
695 const T* aData = this->data();
696 const T* bData = other.
data();
700 T result = zeroVal<T>();
702 if (arraySize < 1024) {
706 for (
SizeType n = 0; n < arraySize; ++n) {
707 result += aData[n] * bData[n];
719 tbb::parallel_for(
SizeRange(0, binCount),
720 DeterministicDotProductOp(aData, bData, binCount, arraySize, partialSums));
722 for (
SizeType n = 0; n < binCount; ++n) {
723 result += partialSums[n];
738 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
739 maxValue =
Max(maxValue,
Abs(data[n]));
753 T result = tbb::parallel_reduce(
SizeRange(0, this->size()), zeroVal<T>(),
754 InfNormOp(this->data()), [](T max1, T max2) {
return Max(max1, max2); });
767 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
768 if (!std::isfinite(data[n]))
return false;
783 bool finite = tbb::parallel_reduce(
SizeRange(0, this->size()),
true,
784 IsFiniteOp(this->data()),
785 [](
bool finite1,
bool finite2) {
return (finite1 && finite2); });
791 template<
typename OtherValueType>
794 EqOp(
const T* a_,
const OtherValueType* b_, T e): a(a_), b(b_), eps(e) {}
796 bool operator()(
const SizeRange& range,
bool equal)
const 799 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
807 const OtherValueType* b;
813 template<
typename OtherValueType>
817 if (this->size() != other.
size())
return false;
818 bool equal = tbb::parallel_reduce(
SizeRange(0, this->size()),
true,
819 EqOp<OtherValueType>(this->data(), other.
data(), eps),
820 [](
bool eq1,
bool eq2) {
return (eq1 && eq2); });
829 std::ostringstream ostr;
832 for (
SizeType n = 0, N = this->size(); n < N; ++n) {
833 ostr << sep << (*this)[n];
844 template<
typename ValueType, SizeType STENCIL_SIZE>
848 template<
typename ValueType, SizeType STENCIL_SIZE>
852 , mValueArray(new
ValueType[mNumRows * STENCIL_SIZE])
853 , mColumnIdxArray(new
SizeType[mNumRows * STENCIL_SIZE])
854 , mRowSizeArray(new
SizeType[mNumRows])
857 tbb::parallel_for(
SizeRange(0, mNumRows),
862 template<
typename ValueType, SizeType STENCIL_SIZE>
866 from(&from_), to(&to_) {}
870 const ValueType* fromVal = from->mValueArray.get();
871 const SizeType* fromCol = from->mColumnIdxArray.get();
872 ValueType* toVal = to->mValueArray.get();
873 SizeType* toCol = to->mColumnIdxArray.get();
874 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
875 toVal[n] = fromVal[n];
876 toCol[n] = fromCol[n];
884 template<
typename ValueType, SizeType STENCIL_SIZE>
887 : mNumRows(other.mNumRows)
888 , mValueArray(new
ValueType[mNumRows * STENCIL_SIZE])
889 , mColumnIdxArray(new
SizeType[mNumRows * STENCIL_SIZE])
890 , mRowSizeArray(new
SizeType[mNumRows])
892 SizeType size = mNumRows * STENCIL_SIZE;
895 tbb::parallel_for(
SizeRange(0, size), MatrixCopyOp(other, *
this));
898 tbb::parallel_for(
SizeRange(0, mNumRows),
903 template<
typename ValueType, SizeType STENCIL_SIZE>
908 assert(row < mNumRows);
909 this->getRowEditor(row).setValue(col, val);
913 template<
typename ValueType, SizeType STENCIL_SIZE>
917 assert(row < mNumRows);
918 return this->getConstRow(row).getValue(col);
922 template<
typename ValueType, SizeType STENCIL_SIZE>
926 return this->getValue(row,col);
930 template<
typename ValueType, SizeType STENCIL_SIZE>
931 template<
typename Scalar>
936 void operator()(
const SizeRange& range)
const 938 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
939 RowEditor row = mat->getRowEditor(n);
949 template<
typename ValueType, SizeType STENCIL_SIZE>
950 template<
typename Scalar>
955 tbb::parallel_for(
SizeRange(0, mNumRows), RowScaleOp<Scalar>(*
this, s));
959 template<
typename ValueType, SizeType STENCIL_SIZE>
960 template<
typename VecValueType>
964 mat(&m), in(i), out(o) {}
966 void operator()(
const SizeRange& range)
const 968 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
969 ConstRow row = mat->getConstRow(n);
970 out[n] = row.dot(in, mat->numRows());
975 const VecValueType* in;
980 template<
typename ValueType, SizeType STENCIL_SIZE>
981 template<
typename VecValueType>
986 if (inVec.
size() != mNumRows) {
988 << mNumRows <<
"x" << mNumRows <<
" vs. " << inVec.
size() <<
")");
990 if (resultVec.
size() != mNumRows) {
992 << mNumRows <<
"x" << mNumRows <<
" vs. " << resultVec.
size() <<
")");
995 vectorMultiply(inVec.
data(), resultVec.
data());
999 template<
typename ValueType, SizeType STENCIL_SIZE>
1000 template<
typename VecValueType>
1003 const VecValueType* inVec, VecValueType* resultVec)
const 1006 tbb::parallel_for(
SizeRange(0, mNumRows),
1007 VecMultOp<VecValueType>(*
this, inVec, resultVec));
1011 template<
typename ValueType, SizeType STENCIL_SIZE>
1012 template<
typename OtherValueType>
1017 a(&a_), b(&b_), eps(e) {}
1019 bool operator()(
const SizeRange& range,
bool equal)
const 1022 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1023 if (!a->getConstRow(n).eq(b->getConstRow(n), eps))
return false;
1031 const ValueType eps;
1035 template<
typename ValueType, SizeType STENCIL_SIZE>
1036 template<
typename OtherValueType>
1041 if (this->numRows() != other.
numRows())
return false;
1042 bool equal = tbb::parallel_reduce(
SizeRange(0, this->numRows()),
true,
1043 EqOp<OtherValueType>(*
this, other, eps),
1044 [](
bool eq1,
bool eq2) {
return (eq1 && eq2); });
1049 template<
typename ValueType, SizeType STENCIL_SIZE>
1057 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1058 const ConstRow row = mat->getConstRow(n);
1060 if (!std::isfinite(*it))
return false;
1071 template<
typename ValueType, SizeType STENCIL_SIZE>
1076 bool finite = tbb::parallel_reduce(
SizeRange(0, this->numRows()),
true,
1077 IsFiniteOp(*
this), [](
bool finite1,
bool finite2) {
return (finite1&&finite2); });
1082 template<
typename ValueType, SizeType STENCIL_SIZE>
1086 std::ostringstream ostr;
1087 for (
SizeType n = 0, N = this->size(); n < N; ++n) {
1088 ostr << n <<
": " << this->getConstRow(n).
str() <<
"\n";
1094 template<
typename ValueType, SizeType STENCIL_SIZE>
1098 assert(i < mNumRows);
1099 const SizeType head = i * STENCIL_SIZE;
1100 return RowEditor(&mValueArray[head], &mColumnIdxArray[head], mRowSizeArray[i], mNumRows);
1104 template<
typename ValueType, SizeType STENCIL_SIZE>
1108 assert(i < mNumRows);
1109 const SizeType head = i * STENCIL_SIZE;
1110 return ConstRow(&mValueArray[head], &mColumnIdxArray[head], mRowSizeArray[i]);
1114 template<
typename ValueType, SizeType STENCIL_SIZE>
1115 template<
typename DataType>
1119 if (this->empty())
return mData.mSize;
1123 const SizeType* colPtr = std::lower_bound(mData.mCols, mData.mCols + mData.mSize, columnIdx);
1125 return static_cast<SizeType>(colPtr - mData.mCols);
1129 template<
typename ValueType, SizeType STENCIL_SIZE>
1130 template<
typename DataType>
1131 inline const ValueType&
1133 SizeType columnIdx,
bool& active)
const 1136 SizeType idx = this->find(columnIdx);
1137 if (idx < this->size() && this->column(idx) == columnIdx) {
1139 return this->value(idx);
1141 return SparseStencilMatrix::sZeroValue;
1144 template<
typename ValueType, SizeType STENCIL_SIZE>
1145 template<
typename DataType>
1146 inline const ValueType&
1149 SizeType idx = this->find(columnIdx);
1150 if (idx < this->size() && this->column(idx) == columnIdx) {
1151 return this->value(idx);
1153 return SparseStencilMatrix::sZeroValue;
1157 template<
typename ValueType, SizeType STENCIL_SIZE>
1158 template<
typename DataType>
1162 return ConstValueIter(mData);
1166 template<
typename ValueType, SizeType STENCIL_SIZE>
1167 template<
typename DataType>
1168 template<
typename OtherDataType>
1171 const RowBase<OtherDataType>& other, ValueType eps)
const 1173 if (this->size() != other.size())
return false;
1174 for (ConstValueIter it = cbegin(), oit = other.cbegin(); it || oit; ++it, ++oit) {
1175 if (it.column() != oit.column())
return false;
1182 template<
typename ValueType, SizeType STENCIL_SIZE>
1183 template<
typename DataType>
1184 template<
typename VecValueType>
1187 const VecValueType* inVec,
SizeType vecSize)
const 1189 VecValueType result = zeroVal<VecValueType>();
1190 for (
SizeType idx = 0, N =
std::min(vecSize, this->size()); idx < N; ++idx) {
1191 result +=
static_cast<VecValueType
>(this->value(idx) * inVec[this->column(idx)]);
1196 template<
typename ValueType, SizeType STENCIL_SIZE>
1197 template<
typename DataType>
1198 template<
typename VecValueType>
1203 return dot(inVec.
data(), inVec.
size());
1207 template<
typename ValueType, SizeType STENCIL_SIZE>
1208 template<
typename DataType>
1212 std::ostringstream ostr;
1214 for (
SizeType n = 0, N = this->size(); n < N; ++n) {
1215 ostr << sep <<
"(" << this->column(n) <<
", " << this->value(n) <<
")";
1222 template<
typename ValueType, SizeType STENCIL_SIZE>
1225 const ValueType* valueHead,
const SizeType* columnHead,
const SizeType& rowSize):
1226 ConstRowBase(ConstRowData(const_cast<ValueType*>(valueHead),
1232 template<
typename ValueType, SizeType STENCIL_SIZE>
1236 RowBase<>(RowData(valueHead, columnHead, rowSize)), mNumColumns(colSize)
1241 template<
typename ValueType, SizeType STENCIL_SIZE>
1246 RowBase<>::mData.mSize = 0;
1250 template<
typename ValueType, SizeType STENCIL_SIZE>
1253 SizeType column,
const ValueType& value)
1255 assert(column < mNumColumns);
1257 RowData& data = RowBase<>::mData;
1261 SizeType offset = this->find(column);
1263 if (offset < data.mSize && data.mCols[offset] == column) {
1265 data.mVals[offset] = value;
1270 assert(data.mSize < this->capacity());
1272 if (offset >= data.mSize) {
1274 data.mVals[data.mSize] = value;
1275 data.mCols[data.mSize] = column;
1278 for (
SizeType i = data.mSize; i > offset; --i) {
1279 data.mVals[i] = data.mVals[i - 1];
1280 data.mCols[i] = data.mCols[i - 1];
1282 data.mVals[offset] = value;
1283 data.mCols[offset] = column;
1291 template<
typename ValueType, SizeType STENCIL_SIZE>
1292 template<
typename Scalar>
1296 for (
int idx = 0, N = this->size(); idx < N; ++idx) {
1297 RowBase<>::mData.mVals[idx] *= s;
1306 template<
typename MatrixType>
1314 using ValueType =
typename MatrixType::ValueType;
1322 tbb::parallel_for(
SizeRange(0, A.numRows()), InitOp(A, mDiag.data()));
1329 const SizeType size = mDiag.size();
1332 assert(r.
size() == size);
1334 tbb::parallel_for(
SizeRange(0, size), ApplyOp(mDiag.data(), r.
data(), z.
data()));
1344 InitOp(
const MatrixType& m, ValueType* v): mat(&m), vec(v) {}
1346 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1347 const ValueType val = mat->
getValue(n, n);
1349 vec[n] =
static_cast<ValueType
>(1.0 / val);
1352 const MatrixType* mat; ValueType* vec;
1358 ApplyOp(
const ValueType* x_,
const ValueType* y_, ValueType* out_):
1359 x(x_), y(y_), out(out_) {}
1361 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) out[n] = x[n] * y[n];
1363 const ValueType *x, *y; ValueType* out;
1372 template<
typename MatrixType>
1376 struct CopyToLowerOp;
1380 using ValueType =
typename MatrixType::ValueType;
1390 , mLowerTriangular(matrix.
numRows())
1391 , mUpperTriangular(matrix.
numRows())
1398 tbb::parallel_for(
SizeRange(0, numRows), CopyToLowerOp(matrix, mLowerTriangular));
1418 mPassedCompatibilityCondition =
true;
1423 ValueType diagonalValue = crow_k.
getValue(k);
1426 if (diagonalValue < 1.e-5) {
1427 mPassedCompatibilityCondition =
false;
1431 diagonalValue =
Sqrt(diagonalValue);
1434 row_k.setValue(k, diagonalValue);
1437 typename MatrixType::ConstRow srcRow = matrix.getConstRow(k);
1438 typename MatrixType::ConstValueIter citer = srcRow.cbegin();
1439 for ( ; citer; ++citer) {
1441 if (ii < k+1)
continue;
1445 row_ii.setValue(k, *citer / diagonalValue);
1450 for ( ; citer; ++citer) {
1452 if (j < k+1)
continue;
1455 ValueType a_jk = row_j.
getValue(k);
1459 typename MatrixType::ConstRow mask = matrix.getConstRow(j);
1460 typename MatrixType::ConstValueIter maskIter = mask.cbegin();
1461 for ( ; maskIter; ++maskIter) {
1463 if (i < j)
continue;
1466 ValueType a_ij = crow_i.
getValue(j);
1467 ValueType a_ik = crow_i.
getValue(k);
1469 a_ij -= a_ik * a_jk;
1471 row_i.setValue(j, a_ij);
1477 tbb::parallel_for(
SizeRange(0, numRows),
1478 TransposeOp(matrix, mLowerTriangular, mUpperTriangular));
1483 bool isValid()
const override {
return mPassedCompatibilityCondition; }
1487 if (!mPassedCompatibilityCondition) {
1493 SizeType size = mLowerTriangular.numRows();
1495 zVec.
fill(zeroVal<ValueType>());
1496 ValueType* zData = zVec.
data();
1498 if (size == 0)
return;
1500 assert(rVec.
size() == size);
1501 assert(zVec.
size() == size);
1504 mTempVec.fill(zeroVal<ValueType>());
1505 ValueType* tmpData = mTempVec.data();
1506 const ValueType* rData = rVec.
data();
1509 for (
SizeType i = 0; i < size; ++i) {
1510 typename TriangularMatrix::ConstRow row = mLowerTriangular.getConstRow(i);
1511 ValueType diagonal = row.
getValue(i);
1512 ValueType dot = row.dot(mTempVec);
1513 tmpData[i] = (rData[i] - dot) / diagonal;
1514 if (!std::isfinite(tmpData[i])) {
1521 for (
SizeType ii = 0; ii < size; ++ii) {
1523 typename TriangularMatrix::ConstRow row = mUpperTriangular.getConstRow(i);
1524 ValueType diagonal = row.
getValue(i);
1525 ValueType dot = row.dot(zVec);
1526 zData[i] = (tmpData[i] - dot) / diagonal;
1527 if (!std::isfinite(zData[i])) {
1538 struct CopyToLowerOp
1540 CopyToLowerOp(
const MatrixType& m,
TriangularMatrix& l): mat(&m), lower(&l) {}
1542 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1543 typename TriangularMatrix::RowEditor outRow = lower->getRowEditor(n);
1545 typename MatrixType::ConstRow inRow = mat->getConstRow(n);
1546 for (
typename MatrixType::ConstValueIter it = inRow.cbegin(); it; ++it) {
1547 if (it.column() > n)
continue;
1548 outRow.setValue(it.column(), *it);
1559 mat(&m), lower(&l), upper(&u) {}
1561 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1562 typename TriangularMatrix::RowEditor outRow = upper->getRowEditor(n);
1565 typename MatrixType::ConstRow inRow = mat->getConstRow(n);
1566 for (
typename MatrixType::ConstValueIter it = inRow.cbegin(); it; ++it) {
1567 const SizeType column = it.column();
1568 if (column < n)
continue;
1569 outRow.setValue(column, lower->getValue(column, n));
1579 bool mPassedCompatibilityCondition;
1586 namespace internal {
1589 template<
typename T>
1591 axpy(
const T& a,
const T* xVec,
const T* yVec, T* resultVec,
SizeType size)
1597 template<
typename T>
1601 assert(xVec.
size() == yVec.
size());
1602 assert(xVec.
size() == result.
size());
1608 template<
typename MatrixOperator,
typename VecValueType>
1611 const VecValueType* b, VecValueType* r)
1614 A.vectorMultiply(x, r);
1620 template<
typename MatrixOperator,
typename T>
1626 assert(x.
size() == A.numRows());
1637 template<
typename PositiveDefMatrix>
1640 const PositiveDefMatrix& Amat,
1644 const State& termination)
1647 return solve(Amat, bVec, xVec, precond, interrupter, termination);
1651 template<
typename PositiveDefMatrix,
typename Interrupter>
1654 const PositiveDefMatrix& Amat,
1658 Interrupter& interrupter,
1659 const State& termination)
1661 using ValueType =
typename PositiveDefMatrix::ValueType;
1670 const SizeType size = Amat.numRows();
1675 if (size != bVec.
size()) {
1677 << size <<
"x" << size <<
" vs. " << bVec.
size() <<
")");
1679 if (size != xVec.
size()) {
1681 << size <<
"x" << size <<
" vs. " << xVec.
size() <<
")");
1690 const ValueType tmp = bVec.
infNorm();
1698 assert(rVec.isFinite());
1710 ValueType rDotZPrev(1);
1717 for ( ; iteration < termination.
iterations; ++iteration) {
1719 if (interrupter.wasInterrupted()) {
1729 precond.
apply(rVec, zVec);
1732 const ValueType rDotZ = rVec.dot(zVec);
1733 assert(std::isfinite(rDotZ));
1735 if (0 == iteration) {
1739 const ValueType beta = rDotZ / rDotZPrev;
1745 Amat.vectorMultiply(pVec, qVec);
1748 const ValueType pAp = pVec.dot(qVec);
1749 assert(std::isfinite(pAp));
1751 const ValueType alpha = rDotZ / pAp;
1761 l2Error = rVec.l2Norm();
1762 minL2Error =
Min(l2Error, minL2Error);
1767 if (l2Error > 2 * minL2Error) {
1798 #endif // OPENVDB_MATH_CONJGRADIENT_HAS_BEEN_INCLUDED void operator()(const SizeRange &range) const
Definition: ConjGradient.h:521
bool isApproxZero(const Type &x)
Return true if x is equal to zero to within the default floating-point comparison tolerance...
Definition: Math.h:320
void increment()
Definition: ConjGradient.h:421
const SizeType arraySize
Definition: ConjGradient.h:685
std::string str() const
Return a string representation of this matrix.
Definition: ConjGradient.h:1084
T operator()(const SizeRange &range, T maxValue) const
Definition: ConjGradient.h:736
MatrixCopyOp(const SparseStencilMatrix &from_, SparseStencilMatrix &to_)
Definition: ConjGradient.h:865
OtherValueType ValueType
Definition: ConjGradient.h:267
void scale(const Scalar &)
Scale all of the entries in this row.
Definition: ConjGradient.h:1294
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:1639
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:610
void swap(Vector &other)
Swap internal storage with another vector, which need not be the same size.
Definition: ConjGradient.h:195
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:395
Coord Abs(const Coord &xyz)
Definition: Coord.h:513
Preconditioner(const SparseStencilMatrix< T, STENCIL_SIZE > &)
Definition: ConjGradient.h:497
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:549
const T val
Definition: ConjGradient.h:541
void computeResidual(const MatrixOperator &A, const Vector< T > &x, const Vector< T > &b, Vector< T > &r)
Compute r = b − Ax.
Definition: ConjGradient.h:1622
#define OPENVDB_LOG_DEBUG_RUNTIME(message)
Log a debugging message in both debug and optimized builds.
Definition: logging.h:294
void fill(const ValueType &value)
Set all elements of this vector to value.
Definition: ConjGradient.h:627
DeterministicDotProductOp(const T *a_, const T *b_, const SizeType binCount_, const SizeType arraySize_, T *reducetmp_)
Definition: ConjGradient.h:660
const T * constData() const
Return a pointer to this vector's elements.
Definition: ConjGradient.h:238
Definition: ConjGradient.h:760
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:109
IncompleteCholeskyPreconditioner(const MatrixType &matrix)
Definition: ConjGradient.h:1388
RowEditor & operator*=(const Scalar &s)
Scale all of the entries in this row.
Definition: ConjGradient.h:462
Vector & operator*=(const Scalar &s)
Multiply each element of this vector by s.
Definition: ConjGradient.h:203
Preconditioner using incomplete Cholesky factorization.
Definition: ConjGradient.h:70
Vector(SizeType n)
Construct a vector of n elements, with uninitialized values.
Definition: ConjGradient.h:174
void apply(const Vector< ValueType > &r, Vector< ValueType > &z) override
Apply this preconditioner to a residue vector: z = M−1r
Definition: ConjGradient.h:1327
std::shared_ptr< T > SharedPtr
Definition: Types.h:139
const SizeType binCount
Definition: ConjGradient.h:684
bool isFinite() const
Return true if all values along the diagonal are finite.
Definition: ConjGradient.h:1338
bool success
Definition: ConjGradient.h:74
ValueType l2Norm() const
Return the L2 norm of this vector.
Definition: ConjGradient.h:212
Read-only accessor to a row of this matrix.
Definition: ConjGradient.h:438
IsFiniteOp(const T *data_)
Definition: ConjGradient.h:762
const T * data
Definition: ConjGradient.h:774
SparseStencilMatrix * to
Definition: ConjGradient.h:880
bool isFinite(const float x)
Return true if x is finite.
Definition: Math.h:346
Information about the state of a conjugate gradient solution.
Definition: ConjGradient.h:73
tbb::blocked_range< SizeType > SizeRange
Definition: ConjGradient.h:62
Vector< ValueType > VectorType
Definition: ConjGradient.h:268
ValueType ValueType
Definition: ConjGradient.h:168
T * data()
Return a pointer to this vector's elements.
Definition: ConjGradient.h:236
Tolerance for floating-point comparison.
Definition: Math.h:117
T * reducetmp
Definition: ConjGradient.h:686
std::string str() const
Return a string representation of this vector.
Definition: ConjGradient.h:827
typename TriangularMatrix::ConstRow TriangleConstRow
Definition: ConjGradient.h:1385
SharedPtr< SparseStencilMatrix > Ptr
Definition: ConjGradient.h:269
double relativeError
Definition: ConjGradient.h:76
FillOp(T *data_, const T &val_)
Definition: ConjGradient.h:534
const T & at(SizeType i) const
Return the value of this vector's ith element.
Definition: ConjGradient.h:229
MatrixType::ValueType ValueType
Definition: ConjGradient.h:494
JacobiPreconditioner(const MatrixType &A)
Definition: ConjGradient.h:1319
Diagonal preconditioner.
Definition: ConjGradient.h:69
SharedPtr< Vector > Ptr
Definition: ConjGradient.h:169
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:136
T * out
Definition: ConjGradient.h:562
virtual bool isValid() const
Definition: ConjGradient.h:500
SizeType column() const
Definition: ConjGradient.h:419
InfNormOp(const T *data_)
Definition: ConjGradient.h:734
SharedPtr< IncompleteCholeskyPreconditioner > Ptr
Definition: ConjGradient.h:1383
const T * data
Definition: ConjGradient.h:744
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 >())
Solve Ax = b via the preconditioned conjugate gradient method.
Definition: ConjGradient.h:1653
Lightweight, variable-length vector.
Definition: ConjGradient.h:64
#define OPENVDB_LOG_WARN(message)
Log a warning message of the form 'someVar << "some text" << ...'.
Definition: logging.h:280
const ValueType & operator*() const
Definition: ConjGradient.h:413
void computeResidual(const MatrixOperator &A, const VecValueType *x, const VecValueType *b, VecValueType *r)
Compute r = b − Ax.
Definition: ConjGradient.h:1610
MatType scale(const Vec3< typename MatType::value_type > &s)
Return a matrix that scales by s.
Definition: Mat.h:647
const TriangularMatrix & upperMatrix() const
Definition: ConjGradient.h:1534
void apply(const Vector< ValueType > &rVec, Vector< ValueType > &zVec) override
Apply this preconditioner to a residue vector: z = M−1r
Definition: ConjGradient.h:1485
Definition: Exceptions.h:40
ValueType infNorm() const
Return the infinity norm of this vector.
Definition: ConjGradient.h:750
const T * b
Definition: ConjGradient.h:683
Definition: Exceptions.h:90
SizeType numRows() const
Return the number of rows in this matrix.
Definition: ConjGradient.h:285
Index32 SizeType
Definition: ConjGradient.h:60
Base class for conjugate gradient preconditioners.
Definition: ConjGradient.h:68
void operator()(const SizeRange &range) const
Definition: ConjGradient.h:536
static const ValueType sZeroValue
Definition: ConjGradient.h:273
SharedPtr< Preconditioner > Ptr
Definition: ConjGradient.h:495
bool isValid() const override
Definition: ConjGradient.h:1483
const ValueType & operator()(SizeType row, SizeType col) const
Return the value at the given coordinates.
Definition: ConjGradient.h:924
bool empty() const
Return true if this vector has no elements.
Definition: ConjGradient.h:188
SizeType size() const
Return the number of rows in this matrix.
Definition: ConjGradient.h:286
void operator()(const SizeRange &range) const
Definition: ConjGradient.h:551
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:715
bool operator()(const SizeRange &range, bool finite) const
Definition: ConjGradient.h:1054
int iterations
Definition: ConjGradient.h:75
State terminationDefaults()
Return default termination conditions for a conjugate gradient solver.
Definition: ConjGradient.h:84
Definition: ConjGradient.h:1050
Read/write accessor to a row of this matrix.
Definition: ConjGradient.h:446
CopyOp(const T *from_, T *to_)
Definition: ConjGradient.h:519
T * data
Definition: ConjGradient.h:540
const SparseStencilMatrix * mat
Definition: ConjGradient.h:1067
const TriangularMatrix & lowerMatrix() const
Definition: ConjGradient.h:1533
const ValueType & getValue(SizeType row, SizeType col) const
Return the value at the given coordinates.
Definition: ConjGradient.h:915
Vector()
Construct an empty vector.
Definition: ConjGradient.h:172
void reset()
Definition: ConjGradient.h:425
const T * data() const
Return a pointer to this vector's elements.
Definition: ConjGradient.h:237
const T & operator[](SizeType i) const
Return the value of this vector's ith element.
Definition: ConjGradient.h:231
void operator()(const SizeRange &range) const
Definition: ConjGradient.h:664
SizeType size() const
Return the number of elements in this vector.
Definition: ConjGradient.h:186
SizeType setValue(SizeType column, const ValueType &value)
Set the value of the entry in the specified column.
Definition: ConjGradient.h:1252
Definition: ConjGradient.h:658
T * to
Definition: ConjGradient.h:526
void clear()
Set the number of entries in this row to zero.
Definition: ConjGradient.h:1243
T & operator[](SizeType i)
Return the value of this vector's ith element.
Definition: ConjGradient.h:230
bool operator()(const SizeRange &range, bool finite) const
Definition: ConjGradient.h:764
const T * a
Definition: ConjGradient.h:682
LinearOp(const T &a_, const T *x_, const T *y_, T *out_)
Definition: ConjGradient.h:549
SharedPtr< JacobiPreconditioner > Ptr
Definition: ConjGradient.h:1317
Definition: ConjGradient.h:517
Vector(SizeType n, const ValueType &val)
Construct a vector of n elements and initialize each element to the given value.
Definition: ConjGradient.h:176
typename TriangularMatrix::RowEditor TriangleRowEditor
Definition: ConjGradient.h:1386
Definition: ConjGradient.h:863
bool isApproxEqual(const Type &a, const Type &b)
Return true if a is equal to b to within the default floating-point comparison tolerance.
Definition: Math.h:358
const T * y
Definition: ConjGradient.h:561
Iterator over the stored values in a row of this matrix.
Definition: ConjGradient.h:410
void axpy(const T &a, const T *xVec, const T *yVec, T *resultVec, SizeType size)
Compute ax + y.
Definition: ConjGradient.h:1591
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:188
SparseStencilMatrix & operator*=(const Scalar &s)
Multiply all elements in the matrix by s;.
Definition: ConjGradient.h:312
void operator()(const SizeRange &range) const
Definition: ConjGradient.h:868
~Vector()
Definition: ConjGradient.h:178
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:308
Definition: ConjGradient.h:732
Definition: ConjGradient.h:547
Sparse, square matrix representing a 3D stencil operator of size STENCIL_SIZE.
Definition: ConjGradient.h:66
double absoluteError
Definition: ConjGradient.h:77
void axpy(const T &a, const Vector< T > &xVec, const Vector< T > &yVec, Vector< T > &result)
Compute ax + y.
Definition: ConjGradient.h:1599
IsFiniteOp(const SparseStencilMatrix &m)
Definition: ConjGradient.h:1052
Definition: Exceptions.h:83
std::ostream & operator<<(std::ostream &os, const State &state)
Definition: ConjGradient.h:572
RowEditor(ValueType *valueHead, SizeType *columnHead, SizeType &rowSize, SizeType colSize)
Definition: ConjGradient.h:1234
Definition: ConjGradient.h:532
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:52
uint32_t Index32
Definition: Types.h:59
ConstValueIter & operator++()
Definition: ConjGradient.h:422
virtual void apply(const Vector< T > &r, Vector< T > &z)=0
Apply this preconditioner to a residue vector: z = M−1r
const T * from
Definition: ConjGradient.h:525