Go to the documentation of this file.
33 template<
class MatrixType>
39 const label nIter =
min(
A.m() - 1,
A.n());
42 for (label
k = 0;
k < nIter; ++
k)
44 const RMatrix u(householderReflector(
A.subColumn(
k,
k)));
46 applyHouseholder(
A, u,
k);
49 if (outputType_ == outputTypes::REDUCED_R)
51 A.resize(
A.n(),
A.n());
56 template<
class MatrixType>
62 const label nCols =
A.n();
63 const label nIter =
min(
A.m() - 1,
A.n());
69 List<scalar> colNorms(nCols);
70 for (label
k = 0;
k < nCols; ++
k)
72 colNorms[
k] =
A.columnNorm(
k,
true);
76 for (label
k = 0;
k < nIter; ++
k)
78 const labelRange colRange(
k, nCols);
79 const SubList<scalar> subColNorms(colNorms, colRange);
82 const label maxColNormi =
86 std::max_element(subColNorms.cbegin(), subColNorms.cend())
92 if (maxColNormi != 0 && SMALL <
mag(
A(
k,
k + maxColNormi)))
94 const RMatrix R1(
A.subColumn(
k));
95 const RMatrix R2(
A.subColumn(maxColNormi +
k));
97 A.subColumn(maxColNormi +
k) = R1;
99 Swap(P_[
k], P_[maxColNormi +
k]);
100 Swap(colNorms[
k], colNorms[maxColNormi +
k]);
104 const RMatrix u(householderReflector(
A.subColumn(
k,
k)));
106 applyHouseholder(
A, u,
k);
113 for (
const auto& val : RMatrix(
A.subRow(
k, q)))
115 colNorms[q] -=
magSqr(val);
121 if (outputType_ == outputTypes::REDUCED_R)
123 A.resize(
A.n(),
A.n());
128 template<
class MatrixType>
129 template<
template<
typename>
class ListContainer>
132 ListContainer<cmptType>&
x
135 const label
n = R_.n();
137 for (label i =
n - 1; 0 <= i; --i)
141 for (label j = i + 1; j <
n; ++j)
143 sum -=
x[j]*R_(i, j);
146 if (
mag(R_(i, i)) < SMALL)
149 <<
"Back-substitution failed due to small diagonal"
158 template<
class MatrixType>
159 template<
template<
typename>
class ListContainer>
163 const ListContainer<cmptType>& source
167 const label m = Q_.m();
170 for (label i = 0; i < m; ++i)
173 for (label j = 0; j < m; ++j)
175 x[i] += Q_(j, i)*source[j];
185 template<
class MatrixType>
194 template<
class MatrixType>
202 outputType_(outputType),
203 storeMethod_(storeMethod),
211 template<
class MatrixType>
220 outputType_(outputType),
221 storeMethod_(storeMethod),
231 template<
class MatrixType>
239 outputType_(outputType),
240 storeMethod_(storeMethods::OUT_OF_PLACE),
252 template<
class MatrixType>
259 if (
A.m() <=
A.n() && outputType_ == outputTypes::REDUCED_R)
261 outputType_ = outputTypes::FULL_R;
265 <<
"Reduced QR decomposition is by definition limited to matrices"
266 <<
" wherein rows > columns, thus computing FULL decomposition."
272 if (outputType_ == outputTypes::FULL_QR)
274 Q_ = MatrixType({
A.m(),
A.m()},
I);
278 switch (storeMethod_)
280 case storeMethods::IN_PLACE:
293 case storeMethods::OUT_OF_PLACE:
311 template<
class MatrixType>
317 if (storeMethod_ == storeMethods::IN_PLACE)
320 <<
"const type qualifier invalidates storeMethods::IN_PLACE." <<
nl;
324 if (
A.m() <=
A.n() && outputType_ == outputTypes::REDUCED_R)
326 outputType_ = outputTypes::FULL_R;
330 <<
"Reduced QR decomposition is by definition limited to matrices"
331 <<
" wherein rows > columns, thus computing FULL decomposition."
337 if (outputType_ == outputTypes::FULL_QR)
339 Q_ = MatrixType({
A.m(),
A.m()},
I);
356 template<
class MatrixType>
363 solveImpl(
x, source);
367 template<
class MatrixType>
375 solveImpl(
x, source);
379 template<
class MatrixType>
386 auto tresult(Q_.Tmul(source));
388 solvex(tresult.ref());
394 template<
class MatrixType>
402 auto tresult(Q_.Tmul(source));
404 solvex(tresult.ref());
410 template<
class MatrixType>
414 const label m = Q_.
m();
419 for (label i = 0; i < m; ++i)
421 for (label j = 0; j < m; ++j)
426 inv.subColumn(i) =
x;
433 template<
class MatrixType>
441 const label mRows =
A.m();
442 const label nCols =
A.n();
443 const label pCols = rhs.
n();
447 const label qRows = rhs.
m();
451 <<
"Linear system is not solvable since the number of rows of"
452 <<
"A and rhs are not equal:" <<
tab << mRows <<
"vs." << qRows
457 for (
const auto& val :
diag)
459 if (
mag(val) < SMALL)
462 <<
"SMALL-valued diagonal elem in back-substitution." <<
nl;
470 for (label i = 0; i < pCols; ++i)
472 for (label j = mRows - 1; -1 < j; --j)
476 for (label
k = j + 1;
k < mRows; ++
k)
491 template<
class MatrixType>
495 const scalar tolerance
498 scalar tol = tolerance;
499 typedef typename MatrixType::cmptType cmptType;
504 <<
"Empty matrix found."
510 if (
A(0,0) == cmptType(0))
512 return MatrixType({1, 1}, cmptType(0));
515 return MatrixType({1, 1}, cmptType(1)/(
A(0,0) + cmptType(VSMALL)));
524 const MatrixType&
R(QRM.
R());
525 const MatrixType& Q(QRM.
Q());
529 label firstZeroElemi = 0;
535 auto lessThan = [=](
const cmptType&
x) {
return tol >
mag(
x); };
544 if (firstZeroElemi == 0)
549 <<
"The largest diagonal element magnitude is nearly zero. "
550 <<
"Tightening the tolerance."
559 <<
"The largest diagonal element magnitude is nearly zero. "
560 <<
"Returning a zero matrix."
564 return MatrixType({
A.n(),
A.m()},
Zero);
593 QRSolve.
Q() & (R1.T())
600 return MatrixType((QRM.
P()^R2)^Q);
624 return MatrixType((QRM.
P()^R2)^Q);
void solve(List< cmptType > &x, const UList< cmptType > &source) const
void Swap(DynamicList< T, SizeMinA > &a, DynamicList< T, SizeMinB > &b)
A class for managing temporary objects.
label n() const noexcept
The number of columns.
static constexpr const zero Zero
Global zero (0)
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
MatrixType pinv(const MatrixType &A, const scalar tol=1e-5)
MatrixType::cmptType cmptType
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
Ostream & endl(Ostream &os)
Add newline and flush stream.
SMatrix inv() const
Return the inverse of (Q*R), so that solving x = (Q*R).inv()*source.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
const MatrixType & Q() const
Return the unitary similarity matrix.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
QRMatrix (i.e. QR decomposition, QR factorisation or orthogonal-triangular decomposition) decomposes ...
void resize(const label m, const label n)
Change Matrix dimensions, preserving the elements.
#define R(A, B, C, D, E, F, K, M)
RMatrix backSubstitution(const SMatrix &A, const RMatrix &b)
Solve a row-echelon-form linear system starting from the bottom.
Generic templated field type.
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
A templated (M x N) rectangular matrix of objects of <Type>, containing M*N elements,...
const MatrixType & R() const
Return the upper triangular matrix.
errorManip< error > abort(error &err)
scalar distance(const vector &p1, const vector &p2)
colPivoting
Options for the computation of column pivoting.
storeMethods
Options where to store R.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
label m() const noexcept
The number of rows.
QRMatrix()
Construct null.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
label k
Boltzmann constant.
const dimensionedScalar e
Elementary charge.
labelList identity(const label len, label start=0)
Create identity map of the given length with (map[i] == i)
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Base for lists with indirect addressing, templated on the list contents type and the addressing type....
outputTypes
Options for the output matrix forms of QRMatrix.
Graphite solid properties.
#define WarningInFunction
Report a warning using Foam::Warning.
static const Identity< scalar > I
void decompose(MatrixType &A)
Compute QR decomposition according to constructor settings.
static bool lessThan(const scalar &val, const scalar &upper)
SMatrix P() const
Create and return the permutation matrix.