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>
37 template<
typename ValueType>
class Vector;
55 template<
typename ValueType>
63 s.
absoluteError = std::numeric_limits<ValueType>::epsilon() * 100.0;
90 template<
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>());
122 template<
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); }
171 void fill(
const ValueType& value);
174 template<
typename Scalar>
void scale(
const Scalar& s);
180 ValueType dot(
const Vector&)
const;
183 ValueType infNorm()
const;
192 template<
typename OtherValueType>
197 std::string str()
const;
200 inline T& at(
SizeType i) {
return mData[i]; }
208 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;
236 template<
typename ValueType_, SizeType STENCIL_SIZE>
244 class ConstValueIter;
257 SizeType numRows()
const {
return mNumRows; }
276 ConstRow getConstRow(
SizeType row)
const;
279 RowEditor getRowEditor(
SizeType row);
282 template<
typename Scalar>
void scale(
const Scalar& s);
284 template<
typename Scalar>
291 template<
typename VecValueType>
298 template<
typename VecValueType>
299 void vectorMultiply(
const VecValueType* inVec, VecValueType* resultVec)
const;
303 template<
typename OtherValueType>
311 std::string str()
const;
315 RowData(ValueType* v,
SizeType* c,
SizeType& s): mVals(v), mCols(c), mSize(s) {}
319 struct ConstRowData {
321 mVals(v), mCols(c), mSize(s) {}
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]; }
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;
432 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];
591 if (mData)
delete[] mData;
607 template<
typename Scalar>
610 ScaleOp(T* data_,
const Scalar& s_): data(data_), s(s_) {}
612 void operator()(
const SizeRange& range)
const {
613 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) data[n] *= s;
622 template<
typename Scalar>
626 tbb::parallel_for(
SizeRange(0, mSize), ScaleOp<Scalar>(mData, s));
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;
647 T sum = zeroVal<T>();
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();
673 T result = zeroVal<T>();
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),
693 DeterministicDotProductOp(aData, bData, binCount, arraySize, partialSums));
695 for (
SizeType n = 0; n < binCount; ++n) {
696 result += partialSums[n];
711 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
712 maxValue =
Max(maxValue,
Abs(data[n]));
726 T result = tbb::parallel_reduce(
SizeRange(0, this->size()), zeroVal<T>(),
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,
757 IsFiniteOp(this->data()),
758 [](
bool finite1,
bool finite2) {
return (finite1 && finite2); });
764 template<
typename OtherValueType>
767 EqOp(
const T* a_,
const OtherValueType* b_, T e): a(a_), b(b_), eps(e) {}
769 bool operator()(
const SizeRange& range,
bool equal)
const
772 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
780 const OtherValueType* b;
786 template<
typename OtherValueType>
790 if (this->size() != other.size())
return false;
791 bool equal = tbb::parallel_reduce(
SizeRange(0, this->size()),
true,
792 EqOp<OtherValueType>(this->data(), other.data(), eps),
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];
817 template<
typename ValueType, SizeType STENCIL_SIZE>
821 template<
typename ValueType, SizeType STENCIL_SIZE>
825 , mValueArray(new
ValueType[mNumRows * STENCIL_SIZE])
826 , mColumnIdxArray(new
SizeType[mNumRows * STENCIL_SIZE])
827 , mRowSizeArray(new
SizeType[mNumRows])
830 tbb::parallel_for(
SizeRange(0, mNumRows),
835 template<
typename ValueType, SizeType STENCIL_SIZE>
839 from(&from_), to(&to_) {}
843 const ValueType* fromVal = from->mValueArray.get();
844 const SizeType* fromCol = from->mColumnIdxArray.get();
845 ValueType* toVal = to->mValueArray.get();
846 SizeType* toCol = to->mColumnIdxArray.get();
847 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
848 toVal[n] = fromVal[n];
849 toCol[n] = fromCol[n];
857 template<
typename ValueType, SizeType STENCIL_SIZE>
860 : mNumRows(other.mNumRows)
861 , mValueArray(new
ValueType[mNumRows * STENCIL_SIZE])
862 , mColumnIdxArray(new
SizeType[mNumRows * STENCIL_SIZE])
863 , mRowSizeArray(new
SizeType[mNumRows])
865 SizeType size = mNumRows * STENCIL_SIZE;
868 tbb::parallel_for(
SizeRange(0, size), MatrixCopyOp(other, *
this));
871 tbb::parallel_for(
SizeRange(0, mNumRows),
876 template<
typename ValueType, SizeType STENCIL_SIZE>
881 assert(row < mNumRows);
882 this->getRowEditor(row).setValue(col, val);
886 template<
typename ValueType, SizeType STENCIL_SIZE>
887 inline const ValueType&
890 assert(row < mNumRows);
891 return this->getConstRow(row).getValue(col);
895 template<
typename ValueType, SizeType STENCIL_SIZE>
896 inline const ValueType&
899 return this->getValue(row,col);
903 template<
typename ValueType, SizeType STENCIL_SIZE>
904 template<
typename Scalar>
909 void operator()(
const SizeRange& range)
const
911 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
912 RowEditor row = mat->getRowEditor(n);
917 SparseStencilMatrix* mat;
922 template<
typename ValueType, SizeType STENCIL_SIZE>
923 template<
typename Scalar>
928 tbb::parallel_for(
SizeRange(0, mNumRows), RowScaleOp<Scalar>(*
this, s));
932 template<
typename ValueType, SizeType STENCIL_SIZE>
933 template<
typename VecValueType>
937 mat(&m), in(i), out(o) {}
939 void operator()(
const SizeRange& range)
const
941 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
942 ConstRow row = mat->getConstRow(n);
943 out[n] = row.dot(in, mat->numRows());
947 const SparseStencilMatrix* mat;
948 const VecValueType* in;
953 template<
typename ValueType, SizeType STENCIL_SIZE>
954 template<
typename VecValueType>
959 if (inVec.size() != mNumRows) {
961 << mNumRows <<
"x" << mNumRows <<
" vs. " << inVec.size() <<
")");
963 if (resultVec.size() != mNumRows) {
965 << mNumRows <<
"x" << mNumRows <<
" vs. " << resultVec.size() <<
")");
968 vectorMultiply(inVec.data(), resultVec.data());
972 template<
typename ValueType, SizeType STENCIL_SIZE>
973 template<
typename VecValueType>
976 const VecValueType* inVec, VecValueType* resultVec)
const
979 tbb::parallel_for(
SizeRange(0, mNumRows),
980 VecMultOp<VecValueType>(*
this, inVec, resultVec));
984 template<
typename ValueType, SizeType STENCIL_SIZE>
985 template<
typename OtherValueType>
990 a(&a_), b(&b_), eps(e) {}
992 bool operator()(
const SizeRange& range,
bool equal)
const
995 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
996 if (!a->getConstRow(n).eq(b->getConstRow(n), eps))
return false;
1002 const SparseStencilMatrix* a;
1003 const SparseStencilMatrix<OtherValueType, STENCIL_SIZE>* b;
1004 const ValueType eps;
1008 template<
typename ValueType, SizeType STENCIL_SIZE>
1009 template<
typename OtherValueType>
1014 if (this->numRows() != other.
numRows())
return false;
1015 bool equal = tbb::parallel_reduce(
SizeRange(0, this->numRows()),
true,
1016 EqOp<OtherValueType>(*
this, other, eps),
1017 [](
bool eq1,
bool eq2) {
return (eq1 && eq2); });
1022 template<
typename ValueType, SizeType STENCIL_SIZE>
1030 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1031 const ConstRow row = mat->getConstRow(n);
1033 if (!std::isfinite(*it))
return false;
1044 template<
typename ValueType, SizeType STENCIL_SIZE>
1049 bool finite = tbb::parallel_reduce(
SizeRange(0, this->numRows()),
true,
1050 IsFiniteOp(*
this), [](
bool finite1,
bool finite2) {
return (finite1&&finite2); });
1055 template<
typename ValueType, SizeType STENCIL_SIZE>
1059 std::ostringstream ostr;
1060 for (
SizeType n = 0, N = this->size(); n < N; ++n) {
1061 ostr << n <<
": " << this->getConstRow(n).
str() <<
"\n";
1067 template<
typename ValueType, SizeType STENCIL_SIZE>
1071 assert(i < mNumRows);
1072 const SizeType head = i * STENCIL_SIZE;
1073 return RowEditor(&mValueArray[head], &mColumnIdxArray[head], mRowSizeArray[i], mNumRows);
1077 template<
typename ValueType, SizeType STENCIL_SIZE>
1081 assert(i < mNumRows);
1082 const SizeType head = i * STENCIL_SIZE;
1083 return ConstRow(&mValueArray[head], &mColumnIdxArray[head], mRowSizeArray[i]);
1087 template<
typename ValueType, SizeType STENCIL_SIZE>
1088 template<
typename DataType>
1092 if (this->empty())
return mData.mSize;
1096 const SizeType* colPtr = std::lower_bound(mData.mCols, mData.mCols + mData.mSize, columnIdx);
1098 return static_cast<SizeType>(colPtr - mData.mCols);
1102 template<
typename ValueType, SizeType STENCIL_SIZE>
1103 template<
typename DataType>
1104 inline const ValueType&
1105 SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::getValue(
1106 SizeType columnIdx,
bool& active)
const
1109 SizeType idx = this->find(columnIdx);
1110 if (idx < this->size() && this->column(idx) == columnIdx) {
1112 return this->value(idx);
1114 return SparseStencilMatrix::sZeroValue;
1117 template<
typename ValueType, SizeType STENCIL_SIZE>
1118 template<
typename DataType>
1119 inline const ValueType&
1120 SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::getValue(
SizeType columnIdx)
const
1122 SizeType idx = this->find(columnIdx);
1123 if (idx < this->size() && this->column(idx) == columnIdx) {
1124 return this->value(idx);
1126 return SparseStencilMatrix::sZeroValue;
1130 template<
typename ValueType, SizeType STENCIL_SIZE>
1131 template<
typename DataType>
1132 inline typename SparseStencilMatrix<ValueType, STENCIL_SIZE>::ConstValueIter
1133 SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::cbegin()
const
1135 return ConstValueIter(mData);
1139 template<
typename ValueType, SizeType STENCIL_SIZE>
1140 template<
typename DataType>
1141 template<
typename OtherDataType>
1143 SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::eq(
1144 const RowBase<OtherDataType>& other, ValueType eps)
const
1146 if (this->size() != other.size())
return false;
1147 for (ConstValueIter it = cbegin(), oit = other.cbegin(); it || oit; ++it, ++oit) {
1148 if (it.column() != oit.column())
return false;
1155 template<
typename ValueType, SizeType STENCIL_SIZE>
1156 template<
typename DataType>
1157 template<
typename VecValueType>
1159 SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::dot(
1160 const VecValueType* inVec,
SizeType vecSize)
const
1162 VecValueType result = zeroVal<VecValueType>();
1163 for (
SizeType idx = 0, N =
std::min(vecSize, this->size()); idx < N; ++idx) {
1164 result +=
static_cast<VecValueType
>(this->value(idx) * inVec[this->column(idx)]);
1169 template<
typename ValueType, SizeType STENCIL_SIZE>
1170 template<
typename DataType>
1171 template<
typename VecValueType>
1173 SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::dot(
1174 const Vector<VecValueType>& inVec)
const
1176 return dot(inVec.data(), inVec.size());
1180 template<
typename ValueType, SizeType STENCIL_SIZE>
1181 template<
typename DataType>
1183 SparseStencilMatrix<ValueType, STENCIL_SIZE>::RowBase<DataType>::str()
const
1185 std::ostringstream ostr;
1187 for (
SizeType n = 0, N = this->size(); n < N; ++n) {
1188 ostr << sep <<
"(" << this->column(n) <<
", " << this->value(n) <<
")";
1195 template<
typename ValueType, SizeType STENCIL_SIZE>
1199 ConstRowBase(ConstRowData(const_cast<
ValueType*>(valueHead),
1205 template<
typename ValueType, SizeType STENCIL_SIZE>
1209 RowBase<>(RowData(valueHead, columnHead, rowSize)), mNumColumns(colSize)
1214 template<
typename ValueType, SizeType STENCIL_SIZE>
1219 RowBase<>::mData.mSize = 0;
1223 template<
typename ValueType, SizeType STENCIL_SIZE>
1228 assert(column < mNumColumns);
1230 RowData& data = RowBase<>::mData;
1234 SizeType offset = this->find(column);
1236 if (offset < data.mSize && data.mCols[offset] == column) {
1238 data.mVals[offset] = value;
1243 assert(data.mSize < this->capacity());
1245 if (offset >= data.mSize) {
1247 data.mVals[data.mSize] = value;
1248 data.mCols[data.mSize] = column;
1251 for (
SizeType i = data.mSize; i > offset; --i) {
1252 data.mVals[i] = data.mVals[i - 1];
1253 data.mCols[i] = data.mCols[i - 1];
1255 data.mVals[offset] = value;
1256 data.mCols[offset] = column;
1264 template<
typename ValueType, SizeType STENCIL_SIZE>
1265 template<
typename Scalar>
1269 for (
int idx = 0, N = this->
size(); idx < N; ++idx) {
1270 RowBase<>::mData.mVals[idx] *= s;
1279 template<
typename MatrixType>
1295 tbb::parallel_for(
SizeRange(0, A.numRows()), InitOp(A, mDiag.data()));
1317 InitOp(
const MatrixType& m,
ValueType* v): mat(&m), vec(v) {}
1319 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1320 const ValueType val = mat->getValue(n, n);
1322 vec[n] =
static_cast<ValueType>(1.0 / val);
1332 x(x_), y(y_), out(out_) {}
1334 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) out[n] = x[n] * y[n];
1345 template<
typename MatrixType>
1346 class IncompleteCholeskyPreconditioner:
public Preconditioner<typename MatrixType::ValueType>
1349 struct CopyToLowerOp;
1363 , mLowerTriangular(matrix.
numRows())
1364 , mUpperTriangular(matrix.
numRows())
1371 tbb::parallel_for(
SizeRange(0,
numRows), CopyToLowerOp(matrix, mLowerTriangular));
1391 mPassedCompatibilityCondition =
true;
1396 ValueType diagonalValue = crow_k.getValue(k);
1399 if (diagonalValue < 1.e-5) {
1400 mPassedCompatibilityCondition =
false;
1404 diagonalValue =
Sqrt(diagonalValue);
1407 row_k.setValue(k, diagonalValue);
1410 typename MatrixType::ConstRow srcRow = matrix.getConstRow(k);
1411 typename MatrixType::ConstValueIter citer = srcRow.cbegin();
1412 for ( ; citer; ++citer) {
1414 if (ii < k+1)
continue;
1418 row_ii.setValue(k, *citer / diagonalValue);
1423 for ( ; citer; ++citer) {
1425 if (j < k+1)
continue;
1432 typename MatrixType::ConstRow mask = matrix.getConstRow(j);
1433 typename MatrixType::ConstValueIter maskIter = mask.cbegin();
1434 for ( ; maskIter; ++maskIter) {
1436 if (i < j)
continue;
1442 a_ij -= a_ik * a_jk;
1444 row_i.setValue(j, a_ij);
1451 TransposeOp(matrix, mLowerTriangular, mUpperTriangular));
1456 bool isValid()
const override {
return mPassedCompatibilityCondition; }
1460 if (!mPassedCompatibilityCondition) {
1468 zVec.
fill(zeroVal<ValueType>());
1471 if (
size == 0)
return;
1477 mTempVec.fill(zeroVal<ValueType>());
1483 typename TriangularMatrix::ConstRow row = mLowerTriangular.getConstRow(i);
1486 tmpData[i] = (rData[i] - dot) / diagonal;
1487 if (!std::isfinite(tmpData[i])) {
1496 typename TriangularMatrix::ConstRow row = mUpperTriangular.getConstRow(i);
1499 zData[i] = (tmpData[i] - dot) / diagonal;
1500 if (!std::isfinite(zData[i])) {
1511 struct CopyToLowerOp
1513 CopyToLowerOp(
const MatrixType& m, TriangularMatrix& l): mat(&m), lower(&l) {}
1515 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1516 typename TriangularMatrix::RowEditor outRow = lower->getRowEditor(n);
1518 typename MatrixType::ConstRow inRow = mat->getConstRow(n);
1519 for (
typename MatrixType::ConstValueIter it = inRow.cbegin(); it; ++it) {
1520 if (it.column() > n)
continue;
1521 outRow.setValue(it.column(), *it);
1525 const MatrixType* mat; TriangularMatrix* lower;
1531 TransposeOp(
const MatrixType& m,
const TriangularMatrix& l, TriangularMatrix& u):
1532 mat(&m), lower(&l), upper(&u) {}
1534 for (
SizeType n = range.begin(), N = range.end(); n < N; ++n) {
1535 typename TriangularMatrix::RowEditor outRow = upper->getRowEditor(n);
1538 typename MatrixType::ConstRow inRow = mat->getConstRow(n);
1539 for (
typename MatrixType::ConstValueIter it = inRow.cbegin(); it; ++it) {
1540 const SizeType column = it.column();
1541 if (column < n)
continue;
1542 outRow.setValue(column, lower->getValue(column, n));
1546 const MatrixType* mat;
const TriangularMatrix* lower; TriangularMatrix* upper;
1549 TriangularMatrix mLowerTriangular;
1550 TriangularMatrix mUpperTriangular;
1551 Vector<ValueType> mTempVec;
1552 bool mPassedCompatibilityCondition;
1559 namespace internal {
1562 template<
typename T>
1570 template<
typename T>
1574 assert(xVec.size() == yVec.size());
1575 assert(xVec.size() == result.size());
1576 axpy(a, xVec.data(), yVec.data(), result.data(), xVec.size());
1581 template<
typename MatrixOperator,
typename VecValueType>
1584 const VecValueType* b, VecValueType* r)
1587 A.vectorMultiply(x, r);
1593 template<
typename MatrixOperator,
typename T>
1597 assert(x.size() == b.size());
1598 assert(x.size() == r.size());
1599 assert(x.size() == A.numRows());
1610 template<
typename PositiveDefMatrix>
1613 const PositiveDefMatrix& Amat,
1617 const State& termination)
1620 return solve(Amat, bVec, xVec, precond, interrupter, termination);
1624 template<
typename PositiveDefMatrix,
typename Interrupter>
1627 const PositiveDefMatrix& Amat,
1631 Interrupter& interrupter,
1632 const State& termination)
1634 using ValueType =
typename PositiveDefMatrix::ValueType;
1648 if (
size != bVec.size()) {
1650 <<
size <<
"x" <<
size <<
" vs. " << bVec.size() <<
")");
1652 if (
size != xVec.size()) {
1654 <<
size <<
"x" <<
size <<
" vs. " << xVec.size() <<
")");
1671 assert(rVec.isFinite());
1690 for ( ; iteration < termination.
iterations; ++iteration) {
1692 if (interrupter.wasInterrupted()) {
1702 precond.apply(rVec, zVec);
1706 assert(std::isfinite(rDotZ));
1708 if (0 == iteration) {
1712 const ValueType beta = rDotZ / rDotZPrev;
1718 Amat.vectorMultiply(pVec, qVec);
1722 assert(std::isfinite(pAp));
1734 l2Error = rVec.l2Norm();
1735 minL2Error =
Min(l2Error, minL2Error);
1740 if (l2Error > 2 * minL2Error) {
1771 #endif // OPENVDB_MATH_CONJGRADIENT_HAS_BEEN_INCLUDED