33 template<
class MatrixType>
39 if (mode_ == modes::FULL)
43 else if (mode_ == modes::ECONOMY)
45 return min(
A.m(),
A.n());
52 template<
class MatrixType>
58 const label
n = AT.m();
59 const label m = AT.n();
61 List<cmptType> Rdiag(
n,
Zero);
63 for (label
k = 0;
k <
n; ++
k)
68 for (label i =
k; i < m; ++i)
83 for (label i =
k; i < m; ++i)
91 for (label j =
k + 1; j <
n; ++j)
95 for (label i =
k; i < m; ++i)
100 if (
mag(AT(
k,
k)) > SMALL)
105 for (label i =
k; i < m; ++i)
107 AT(j,i) +=
s*AT(
k,i);
121 template<
class MatrixType>
128 const label
n = AT.m();
129 const label m = AT.n();
130 const label sz =
min(m,
n);
136 List<scalar> norms(
n,
Zero);
137 for (label
k = 0;
k <
n; ++
k)
139 for (label i = 0; i < m; ++i)
145 List<cmptType> Rdiag(
n,
Zero);
147 for (label
k = 0;
k < sz; ++
k)
149 const auto it = std::next(norms.cbegin(),
k);
150 const label maxNormi =
154 std::max_element(it, norms.cend())
160 if (
mag(AT(
k + maxNormi,
k)) > SMALL && maxNormi != 0)
162 const RMatrix R1(AT.subRow(
k));
163 const RMatrix R2(AT.subRow(maxNormi +
k));
166 AT.subRow(maxNormi +
k) = R1;
168 std::swap(p_[
k], p_[maxNormi +
k]);
169 std::swap(norms[
k], norms[maxNormi +
k]);
176 for (label i =
k; i < m; ++i)
181 if (nrm != scalar(0))
184 if (
mag(AT(
k,
k)) < 0)
189 for (label i =
k; i < m; ++i)
194 AT(
k,
k) += scalar(1);
197 for (label j =
k + 1; j <
n; ++j)
201 for (label i =
k; i < m; ++i)
206 if (
mag(AT(
k,
k)) > SMALL)
211 for (label i =
k; i < m; ++i)
213 AT(j,i) +=
s*AT(
k,i);
225 for (
const auto& val : RMatrix(AT.subColumn(
k, q)))
239 template<
class MatrixType>
245 if (output_ == outputs::ONLY_R)
250 const label
n = AT.m();
251 const label m = AT.n();
255 MatrixType QT(Q_.transpose());
257 for (label
k = sz_ - 1;
k >= 0; --
k)
259 for (label i = 0; i < m; ++i)
266 for (label j =
k; j < sz_; ++j)
268 if (
n >
k &&
mag(AT(
k,
k)) > SMALL)
272 for (label i =
k; i < m; ++i)
274 s += AT(
k,i)*QT(j,i);
279 for (label i =
k; i < m; ++i)
291 template<
class MatrixType>
294 const MatrixType& AT,
295 const List<cmptType>&
diag 298 if (output_ == outputs::ONLY_Q)
303 const label
n = AT.m();
307 for (label i = 0; i < sz_; ++i)
309 for (label j = 0; j <
n; ++j)
324 template<
class MatrixType>
325 template<
template<
typename>
class ListContainer>
328 ListContainer<cmptType>&
x 331 const label
n = R_.n();
333 for (label i =
n - 1; 0 <= i; --i)
337 for (label j = i + 1; j <
n; ++j)
339 sum -=
x[j]*R_(i, j);
342 if (
mag(R_(i, i)) < SMALL)
345 <<
"Back-substitution failed due to small diagonal" 354 template<
class MatrixType>
355 template<
template<
typename>
class ListContainer>
359 const ListContainer<cmptType>& source
363 const label m = Q_.m();
366 for (label i = 0; i < m; ++i)
369 for (label j = 0; j < m; ++j)
371 x[i] += Q_(j, i)*source[j];
381 template<
class MatrixType>
392 sz_(calcShapeFactor(
A)),
398 MatrixType AT(
A.transpose());
414 template<
class MatrixType>
425 sz_(calcShapeFactor(
A)),
430 MatrixType AT(
A.transpose());
447 template<
class MatrixType>
454 solveImpl(
x, source);
458 template<
class MatrixType>
466 solveImpl(
x, source);
470 template<
class MatrixType>
477 auto tresult(Q_.Tmul(source));
479 solvex(tresult.ref());
485 template<
class MatrixType>
493 auto tresult(Q_.Tmul(source));
495 solvex(tresult.ref());
501 template<
class MatrixType>
508 const label mRows = R_.m();
509 const label nCols = R_.n();
510 const label pCols = rhs.n();
514 const label qRows = rhs.m();
518 <<
"Linear system is not solvable since the number of rows of" 519 <<
"A and rhs are not equal:" <<
tab << mRows <<
"vs." << qRows
524 for (
const auto& val :
diag)
526 if (
mag(val) < SMALL)
529 <<
"SMALL-valued diagonal elem in back-substitution." 536 RMatrix
b({nCols, pCols},
Zero);
538 for (label i = 0; i < pCols; ++i)
540 for (label j = mRows - 1; -1 < j; --j)
542 cmptType
alpha(rhs(j, i));
544 for (label
k = j + 1;
k < mRows; ++
k)
557 template<
class MatrixType>
561 const label m = Q_.
m();
566 for (label i = 0; i < m; ++i)
568 for (label j = 0; j < m; ++j)
573 inv.subColumn(i) =
x;
582 template<
class MatrixType>
589 typedef typename MatrixType::cmptType cmptType;
594 <<
"Empty matrix found." 600 if (
A(0,0) == cmptType(0))
602 return MatrixType({1, 1}, cmptType(0));
605 return MatrixType({1, 1}, cmptType(1)/(
A(0,0) + cmptType(VSMALL)));
608 QRMatrix<MatrixType> QRM
611 QRMatrix<MatrixType>::modes::FULL,
612 QRMatrix<MatrixType>::outputs::BOTH_QR,
613 QRMatrix<MatrixType>::pivoting::TRUE
615 const MatrixType&
R = QRM.R();
616 const MatrixType& Q = QRM.Q();
624 const List<cmptType>
diag(
R.diag());
626 auto lessThan = [=](
const cmptType&
x) {
return tol >
mag(
x); };
640 <<
"The largest diagonal element magnitude is nearly zero. " 641 <<
"Tightening the tolerance." 650 <<
"The largest diagonal element magnitude is nearly zero. " 651 <<
"Returning a zero matrix." 655 return MatrixType({
A.n(),
A.m()},
Zero);
666 const RectangularMatrix<cmptType> R1(
R.subMatrix(0, 0, elemi));
671 const SquareMatrix<cmptType>
C(R1&R1);
673 QRMatrix<SquareMatrix<cmptType>> QRSolve
676 QRMatrix<SquareMatrix<cmptType>>::modes::FULL,
677 QRMatrix<SquareMatrix<cmptType>>::outputs::BOTH_QR
680 RectangularMatrix<cmptType> R2
689 R2.resize(
R.m(),
R.n());
691 return MatrixType((QRM.P()^R2)^Q);
695 const SquareMatrix<cmptType>
C(R1^R1);
697 QRMatrix<SquareMatrix<cmptType>> QRSolve
700 QRMatrix<SquareMatrix<cmptType>>::modes::FULL,
701 QRMatrix<SquareMatrix<cmptType>>::outputs::BOTH_QR
704 RectangularMatrix<cmptType> R2
713 R2.resize(
R.m(),
R.n());
715 return MatrixType((QRM.P()^R2)^Q);
pivoting
Options for the column pivoting.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &f1)
std::enable_if< !std::is_same< complex, T >::value, const T &>::type conj(const T &val)
The 'conjugate' of non-complex returns itself (pass-through) it does not return a complex! ...
QRMatrix computes QR decomposition of a given scalar/complex matrix A into the following: ...
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
void solve(List< cmptType > &x, const UList< cmptType > &source) const
Solve the linear system with the given source and return the solution in the argument x...
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
static bool lessThan(const scalar &val, const scalar &upper)
scalar distance(const vector &p1, const vector &p2)
constexpr char tab
The tab '\t' character(0x09)
label k
Boltzmann constant.
Base for lists with indirect addressing, templated on the list contents type and the addressing type...
SMatrix inv() const
Return the inverse of (Q*R), solving x = (Q*R).inv()*source.
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
const dimensionedScalar e
Elementary charge.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i), works like std::iota() but returning a...
Generic templated field type.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
A templated (M x N) rectangular matrix of objects of <Type>, containing M*N elements, derived from Matrix.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
label m() const noexcept
The number of rows.
#define R(A, B, C, D, E, F, K, M)
#define WarningInFunction
Report a warning using Foam::Warning.
label find_if(const ListType &input, const UnaryPredicate &pred, const label start=0)
Find index of the first occurrence that satisfies the predicate.
static Ostream & output(Ostream &os, const IntRange< T > &range)
A class for managing temporary objects.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
dimensionedScalar hypot(const dimensionedScalar &x, const dimensionedScalar &y)
mode_t mode(const fileName &name, const bool followLink=true)
Return the file mode, normally following symbolic links.
QRMatrix()=delete
No default construct.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
static constexpr const zero Zero
Global zero (0)