diff --git a/applications/test/Matrix/Test-Matrix.C b/applications/test/Matrix/Test-Matrix.C
index c372cf9b9a7a03ffc137bb2419b9bcb7d86d57ce..1cbc60254538df916b33a17d12bb5a579795fe18 100644
--- a/applications/test/Matrix/Test-Matrix.C
+++ b/applications/test/Matrix/Test-Matrix.C
@@ -37,6 +37,7 @@ License
using namespace Foam;
using namespace Foam::MatrixTools;
+#define RUNALL true
#define isEqual MatrixTools::equal
void horizontalLine()
@@ -113,7 +114,7 @@ int main(int argc, char *argv[])
Random rndGen(1234);
- #if 1
+ #if (0 | RUNALL)
{
horizontalLine();
@@ -218,7 +219,7 @@ int main(int argc, char *argv[])
#endif
- #if 1
+ #if (0 | RUNALL)
{
horizontalLine();
@@ -254,7 +255,7 @@ int main(int argc, char *argv[])
#endif
- #if 1
+ #if (0 | RUNALL)
{
horizontalLine();
@@ -299,7 +300,7 @@ int main(int argc, char *argv[])
#endif
- #if 1
+ #if (0 | RUNALL)
{
horizontalLine();
@@ -410,7 +411,7 @@ int main(int argc, char *argv[])
#endif
- #if 1
+ #if (0 | RUNALL)
{
horizontalLine();
@@ -494,7 +495,7 @@ int main(int argc, char *argv[])
#endif
- #if 1
+ #if (0 | RUNALL)
{
horizontalLine();
@@ -683,7 +684,7 @@ int main(int argc, char *argv[])
#endif
- #if 1
+ #if (0 | RUNALL)
{
horizontalLine();
@@ -804,7 +805,7 @@ int main(int argc, char *argv[])
<< nl;
}
- #if 1
+ #if (0 | RUNALL)
{
Info<< nl << "# Implicit inner/outer products:" << nl;
@@ -860,7 +861,65 @@ int main(int argc, char *argv[])
#endif
- #if 1
+ #if (0 | RUNALL)
+ {
+ horizontalLine();
+
+ Info<< "## Query to determine if a square matrix is symmetric:" << nl;
+
+ {
+ const label mRows = 100;
+ SMatrix A(makeRandomMatrix<SMatrix>({mRows, mRows}, rndGen));
+
+ Info<< "# SquareMatrix<scalar>.symmetric():" << nl
+ << A.symmetric() << nl;
+
+ // Symmetrise
+ for (label n = 0; n < A.n() - 1; ++n)
+ {
+ for (label m = A.m() - 1; n < m; --m)
+ {
+ A(n, m) = A(m, n);
+ }
+ }
+
+ Info<< "# SquareMatrix<scalar>.symmetric():" << nl
+ << A.symmetric() << nl;
+ }
+
+ {
+ const label mRows = 100;
+
+ SCMatrix A(mRows);
+
+ for (auto& val : A)
+ {
+ val.Re() = rndGen.GaussNormal<scalar>();
+ val.Im() = rndGen.GaussNormal<scalar>();
+ }
+
+ Info<< "# SquareMatrix<complex>.symmetric():" << nl
+ << A.symmetric() << nl;
+
+ // Symmetrise
+ for (label n = 0; n < A.n() - 1; ++n)
+ {
+ for (label m = A.m() - 1; n < m; --m)
+ {
+ A(n, m) = A(m, n);
+ }
+ }
+
+ Info<< "# SquareMatrix<complex>.symmetric():" << nl
+ << A.symmetric() << nl;
+ }
+
+ horizontalLine();
+ }
+ #endif
+
+
+ #if (0 | RUNALL)
{
horizontalLine();
@@ -898,7 +957,7 @@ int main(int argc, char *argv[])
#endif
- #if 1
+ #if (0 | RUNALL)
{
scalarSquareMatrix squareMatrix(3, Zero);
diff --git a/applications/test/QRMatrix/Make/files b/applications/test/QRMatrix/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..9825bbfcf351043bd03d02a5d50dd924e73a1d41
--- /dev/null
+++ b/applications/test/QRMatrix/Make/files
@@ -0,0 +1,3 @@
+Test-QRMatrix.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-QRMatrix
diff --git a/applications/test/QRMatrix/Make/options b/applications/test/QRMatrix/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..4e772fdf9d7bc94221d127458f9d2ca32850fe69
--- /dev/null
+++ b/applications/test/QRMatrix/Make/options
@@ -0,0 +1,2 @@
+/* EXE_INC = -I$(LIB_SRC)/finiteVolume/lnInclude */
+/* EXE_LIBS = -lfiniteVolume */
diff --git a/applications/test/QRMatrix/Test-QRMatrix.C b/applications/test/QRMatrix/Test-QRMatrix.C
new file mode 100644
index 0000000000000000000000000000000000000000..fb22d78bde3f9c8316a2a51485fdd9d61a5daaf7
--- /dev/null
+++ b/applications/test/QRMatrix/Test-QRMatrix.C
@@ -0,0 +1,1484 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "MatrixTools.H"
+#include "QRMatrix.H"
+#include "Random.H"
+#include "IOmanip.H"
+
+using namespace Foam;
+using namespace Foam::MatrixTools;
+
+#define equal MatrixTools::equal
+#define RUNALL true
+const bool verbose = true;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void horizontalLine()
+{
+ Info<< "+---------+---------+---------+---------+---------+" << nl;
+}
+
+
+// Create random scalar-type matrix
+template<class MatrixType>
+typename std::enable_if
+<
+ !std::is_same<complex, typename MatrixType::cmptType>::value,
+ MatrixType
+>::type makeRandomMatrix
+(
+ const labelPair& dims,
+ Random& rndGen
+)
+{
+ MatrixType mat(dims);
+
+ std::generate
+ (
+ mat.begin(),
+ mat.end(),
+ [&]{ return rndGen.GaussNormal<typename MatrixType::cmptType>(); }
+ );
+
+ return mat;
+}
+
+
+// Create random complex-type matrix
+template<class MatrixType>
+typename std::enable_if
+<
+ std::is_same<complex, typename MatrixType::cmptType>::value,
+ MatrixType
+>::type makeRandomMatrix
+(
+ const labelPair& dims,
+ Random& rndGen
+)
+{
+ MatrixType mat(dims);
+
+ std::generate
+ (
+ mat.begin(),
+ mat.end(),
+ [&]
+ {
+ return complex
+ (
+ rndGen.GaussNormal<scalar>(),
+ rndGen.GaussNormal<scalar>()
+ );
+ }
+ );
+
+ return mat;
+}
+
+
+// Print OpenFOAM matrix in NumPy format
+template<class MatrixType>
+void InfoNumPyFormat
+(
+ const MatrixType& mat,
+ const word objName
+)
+{
+ Info<< objName << ": " << mat.m() << "x" << mat.n() << nl;
+ for (label m = 0; m < mat.m(); ++m)
+ {
+ Info<< "[";
+ for (label n = 0; n < mat.n(); ++n)
+ {
+ if (n == mat.n() - 1)
+ {
+ Info<< mat(m,n);
+ }
+ else
+ {
+ Info<< mat(m,n) << ",";
+ }
+ }
+ if (m == mat.m() - 1)
+ {
+ Info << "]" << nl;
+ }
+ else
+ {
+ Info << "]," << nl;
+ }
+ }
+}
+
+
+// Returns true if two scalars are equal within a given tolerance, and v.v.
+bool isEqual
+(
+ const scalar a,
+ const scalar b,
+ const bool verbose = false,
+ const scalar relTol = 1e-5,
+ const scalar absTol = 1e-8
+)
+{
+ if ((absTol + relTol*mag(b)) < mag(a - b))
+ {
+ if (verbose)
+ {
+ Info<< "Elements are not close in terms of tolerances:"
+ << nl << a << tab << b << nl;
+ }
+ return false;
+ }
+
+ if (verbose)
+ {
+ Info<< "All elems are the same within the tolerances" << nl;
+ }
+
+ return true;
+}
+
+
+// Prints true if a given matrix is empty, and v.v.
+template<class MatrixType>
+void isEmpty
+(
+ const MatrixType& A,
+ const word objName
+)
+{
+ Info<< "Empty " << objName << " = ";
+ if (A.empty())
+ {
+ Info<< "true" << nl;
+ }
+ else
+ {
+ Info<< "false" << nl;
+ }
+}
+
+
+// Checks if Q matrix is a unitary matrix
+template<class MatrixType>
+void cross_check_QRMatrix
+(
+ const MatrixType& Aorig,
+ const MatrixType& Q
+)
+{
+ InfoNumPyFormat(Q, "Q");
+
+ // mathworld.wolfram.com/UnitaryMatrix.html (Retrieved:18-06-19)
+ Info<< nl << "# Q.T()*Q = Q*Q.T() = I:" << nl;
+ const MatrixType QTQ(Q&Q);
+ const MatrixType QQT(Q^Q);
+ MatrixType IMatrix({Q.m(), Q.n()}, Zero);
+ for (label i = 0; i < IMatrix.n(); ++i)
+ {
+ IMatrix(i,i) = pTraits<typename MatrixType::cmptType>::one;
+ }
+ equal(QTQ, IMatrix, verbose);
+ equal(QQT, IMatrix, verbose);
+ equal(QTQ, QQT, verbose);
+}
+
+
+// Checks if R matrix is an upper-triangular matrix
+template<class MatrixType>
+void cross_check_QRMatrix
+(
+ const MatrixType& R
+)
+{
+ InfoNumPyFormat(R, "R");
+
+ // mathworld.wolfram.com/UpperTriangularMatrix.html (Retrieved:18-06-19)
+ Info<< nl << "# R(i, j) = 0 for i > j:" << nl;
+ for (label i = 0; i < R.m(); ++i)
+ {
+ for (label j = 0; j < R.n(); ++j)
+ {
+ if (j < i)
+ {
+ isEqual(mag(R(i, j)), 0.0, verbose);
+ }
+ }
+ }
+}
+
+
+// Checks if the given matrix can be reconstructed by Q and R matrices
+template<class MatrixType>
+void cross_check_QRMatrix
+(
+ const MatrixType& Aorig,
+ const MatrixType& Q,
+ const MatrixType& R
+)
+{
+ InfoNumPyFormat(Aorig, "Aorig");
+
+ cross_check_QRMatrix(Aorig, Q);
+
+ cross_check_QRMatrix(R);
+
+ // mathworld.wolfram.com/QRDecomposition.html (Retrieved:18-06-19)
+ Info<< nl << "# Q*R = A:" << nl;
+ const MatrixType AReconstruct(Q*R);
+ equal(Aorig, AReconstruct, verbose);
+}
+
+
+// Checks if R matrix is an upper-triangular matrix, and obeys column pivoting
+template<class MatrixType>
+void cross_check_QRMatrix
+(
+ const MatrixType& Aorig,
+ const labelList& orderP,
+ const SquareMatrix<typename MatrixType::cmptType>& P,
+ const MatrixType& R
+)
+{
+ InfoNumPyFormat(Aorig, "Aorig");
+
+ InfoNumPyFormat(P, "P");
+
+ cross_check_QRMatrix(R);
+
+ // w.wiki/54m (Retrieved:20-06-19)
+ Info<< nl << "# Diag elems of R non-increasing:"
+ << "|R11| >= |R22| >= ... >= |Rnn|:" << nl;
+ const List<scalar> diag0(mag(R.diag()));
+ List<scalar> diag1(diag0);
+ Foam::sort(diag1, std::greater<scalar>());
+ forAll(diag0, i)
+ {
+ isEqual(diag0[i], diag1[i], verbose);
+ }
+}
+
+
+// Checks if the given matrix can be reconstructed by column-pivoted Q and R
+template<class MatrixType>
+void cross_check_QRMatrix
+(
+ const MatrixType& Aorig,
+ const MatrixType& Q,
+ const labelList& orderP,
+ const SquareMatrix<typename MatrixType::cmptType>& P,
+ const MatrixType& R
+)
+{
+ cross_check_QRMatrix(Aorig, orderP, P, R);
+ cross_check_QRMatrix(Aorig, Q);
+
+ // w.wiki/54m (Retrieved:20-06-19)
+ Info<< nl << "# Q*R = A*P:" << nl;
+ const MatrixType AReconstruct(Q*R);
+ const MatrixType AP(Aorig*P);
+ equal(AP, AReconstruct, verbose);
+}
+
+
+// Checks each constructor of QRMatrix type
+template<class MatrixType>
+void verification_QRMatrix
+(
+ MatrixType& A
+)
+{
+ typedef SquareMatrix<typename MatrixType::cmptType> SMatrix;
+
+ // Create working copies of matrix A
+ const MatrixType Aorig(A);
+
+ // QRMatrix Constructors
+ #if (0 | RUNALL)
+ {
+ QRMatrix<MatrixType> QRNull();
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# FULL_R | OUT_OF_PLACE" << nl;
+ QRMatrix<MatrixType> QRM
+ (
+ QRMatrix<MatrixType>::outputTypes::FULL_R,
+ QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE
+ );
+ QRM.decompose(A);
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ isEmpty(Q, "Q");
+ cross_check_QRMatrix(R);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# FULL_R | IN_PLACE" << nl;
+ QRMatrix<MatrixType> QRM
+ (
+ QRMatrix<MatrixType>::outputTypes::FULL_R,
+ QRMatrix<MatrixType>::storeMethods::IN_PLACE
+ );
+ MatrixType A0(A);
+ QRM.decompose(A0);
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ isEmpty(Q, "Q");
+ isEmpty(R, "R");
+ cross_check_QRMatrix(A0);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# FULL_QR | OUT_OF_PLACE" << nl;
+ QRMatrix<MatrixType> QRM
+ (
+ QRMatrix<MatrixType>::outputTypes::FULL_QR,
+ QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE
+ );
+ QRM.decompose(A);
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ cross_check_QRMatrix(Aorig, Q, R);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# FULL_QR | IN_PLACE" << nl;
+ QRMatrix<MatrixType> QRM
+ (
+ QRMatrix<MatrixType>::outputTypes::FULL_QR,
+ QRMatrix<MatrixType>::storeMethods::IN_PLACE
+ );
+ MatrixType A0(A);
+ QRM.decompose(A0);
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ isEmpty(R, "R");
+ cross_check_QRMatrix(Aorig, Q, A0);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# FULL_R | OUT_OF_PLACE | colPivot = on" << nl;
+ QRMatrix<MatrixType> QRM
+ (
+ QRMatrix<MatrixType>::outputTypes::FULL_R,
+ QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE,
+ QRMatrix<MatrixType>::colPivoting::TRUE
+ );
+ QRM.decompose(A);
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ const labelList& PList(QRM.orderP());
+ const SMatrix P(QRM.P());
+ isEmpty(Q, "Q");
+ cross_check_QRMatrix(Aorig, PList, P, R);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# FULL_R | IN_PLACE | colPivot = on" << nl;
+ QRMatrix<MatrixType> QRM
+ (
+ QRMatrix<MatrixType>::outputTypes::FULL_R,
+ QRMatrix<MatrixType>::storeMethods::IN_PLACE,
+ QRMatrix<MatrixType>::colPivoting::TRUE
+ );
+ MatrixType A0(A);
+ QRM.decompose(A0);
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ const labelList& PList(QRM.orderP());
+ const SMatrix P(QRM.P());
+ isEmpty(Q, "Q");
+ isEmpty(R, "R");
+ cross_check_QRMatrix(Aorig, PList, P, A0);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# FULL_QR | OUT_OF_PLACE | colPivot = on" << nl;
+ QRMatrix<MatrixType> QRM
+ (
+ QRMatrix<MatrixType>::outputTypes::FULL_QR,
+ QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE,
+ QRMatrix<MatrixType>::colPivoting::TRUE
+ );
+ QRM.decompose(A);
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ const labelList& PList(QRM.orderP());
+ const SMatrix P(QRM.P());
+ cross_check_QRMatrix(Aorig, Q, PList, P, R);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# FULL_QR | IN_PLACE | colPivot = on" << nl;
+ QRMatrix<MatrixType> QRM
+ (
+ QRMatrix<MatrixType>::outputTypes::FULL_QR,
+ QRMatrix<MatrixType>::storeMethods::IN_PLACE,
+ QRMatrix<MatrixType>::colPivoting::TRUE
+ );
+ MatrixType A0(A);
+ QRM.decompose(A0);
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ const labelList& PList(QRM.orderP());
+ const SMatrix P(QRM.P());
+ isEmpty(R, "R");
+ cross_check_QRMatrix(Aorig, Q, PList, P, A0);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# A | FULL_R | OUT_OF_PLACE | colPivot = off" << nl;
+ MatrixType A0(A);
+ QRMatrix<MatrixType> QRM
+ (
+ A0,
+ QRMatrix<MatrixType>::outputTypes::FULL_R,
+ QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE,
+ QRMatrix<MatrixType>::colPivoting::FALSE
+ );
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ isEmpty(Q, "Q");
+ cross_check_QRMatrix(R);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# A | FULL_R | IN_PLACE | colPivot = off" << nl;
+ MatrixType A0(A);
+ QRMatrix<MatrixType> QRM
+ (
+ A0,
+ QRMatrix<MatrixType>::outputTypes::FULL_R,
+ QRMatrix<MatrixType>::storeMethods::IN_PLACE,
+ QRMatrix<MatrixType>::colPivoting::FALSE
+ );
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ isEmpty(Q, "Q");
+ isEmpty(R, "R");
+ cross_check_QRMatrix(A0);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# A | FULL_QR | OUT_OF_PLACE | colPivot = off" << nl;
+ MatrixType A0(A);
+ QRMatrix<MatrixType> QRM
+ (
+ A0,
+ QRMatrix<MatrixType>::outputTypes::FULL_QR,
+ QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE,
+ QRMatrix<MatrixType>::colPivoting::FALSE
+ );
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ cross_check_QRMatrix(Aorig, Q, R);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# A | FULL_QR | IN_PLACE | colPivot = off" << nl;
+ MatrixType A0(A);
+ QRMatrix<MatrixType> QRM
+ (
+ A0,
+ QRMatrix<MatrixType>::outputTypes::FULL_QR,
+ QRMatrix<MatrixType>::storeMethods::IN_PLACE,
+ QRMatrix<MatrixType>::colPivoting::FALSE
+ );
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ isEmpty(R, "R");
+ cross_check_QRMatrix(Aorig, Q, A0);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# A | FULL_R | OUT_OF_PLACE | colPivot = on" << nl;
+ MatrixType A0(A);
+ QRMatrix<MatrixType> QRM
+ (
+ A0,
+ QRMatrix<MatrixType>::outputTypes::FULL_R,
+ QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE,
+ QRMatrix<MatrixType>::colPivoting::TRUE
+ );
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ const labelList& PList(QRM.orderP());
+ const SMatrix P(QRM.P());
+ isEmpty(Q, "Q");
+ cross_check_QRMatrix(Aorig, PList, P, R);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# A | FULL_R | IN_PLACE | colPivot = on" << nl;
+ MatrixType A0(A);
+ QRMatrix<MatrixType> QRM
+ (
+ A0,
+ QRMatrix<MatrixType>::outputTypes::FULL_R,
+ QRMatrix<MatrixType>::storeMethods::IN_PLACE,
+ QRMatrix<MatrixType>::colPivoting::TRUE
+ );
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ const labelList& PList(QRM.orderP());
+ const SMatrix P(QRM.P());
+ isEmpty(Q, "Q");
+ isEmpty(R, "R");
+ cross_check_QRMatrix(Aorig, PList, P, A0);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# A | FULL_QR | OUT_OF_PLACE | colPivot = on" << nl;
+ MatrixType A0(A);
+ QRMatrix<MatrixType> QRM
+ (
+ A0,
+ QRMatrix<MatrixType>::outputTypes::FULL_QR,
+ QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE,
+ QRMatrix<MatrixType>::colPivoting::TRUE
+ );
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ const labelList& PList(QRM.orderP());
+ const SMatrix P(QRM.P());
+ cross_check_QRMatrix(Aorig, Q, PList, P, R);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# A | FULL_QR | IN_PLACE | colPivot = on" << nl;
+ MatrixType A0(A);
+ QRMatrix<MatrixType> QRM
+ (
+ A0,
+ QRMatrix<MatrixType>::outputTypes::FULL_QR,
+ QRMatrix<MatrixType>::storeMethods::IN_PLACE,
+ QRMatrix<MatrixType>::colPivoting::TRUE
+ );
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ const labelList& PList(QRM.orderP());
+ const SMatrix P(QRM.P());
+ isEmpty(R, "R");
+ cross_check_QRMatrix(Aorig, Q, PList, P, A0);
+ }
+ #endif
+
+ // constructors with the const type specifier
+ #if (0 | RUNALL)
+ {
+ Info<< "# const A | FULL_R | colPivot = off" << nl;
+ const MatrixType A0(A);
+ QRMatrix<MatrixType> QRM
+ (
+ A0,
+ QRMatrix<MatrixType>::outputTypes::FULL_R,
+ QRMatrix<MatrixType>::colPivoting::FALSE
+ );
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ isEmpty(Q, "Q");
+ cross_check_QRMatrix(R);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# const A | FULL_QR | colPivot = off" << nl;
+ const MatrixType A0(A);
+ QRMatrix<MatrixType> QRM
+ (
+ A0,
+ QRMatrix<MatrixType>::outputTypes::FULL_QR,
+ QRMatrix<MatrixType>::colPivoting::FALSE
+ );
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ cross_check_QRMatrix(Aorig, Q, R);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# const A | FULL_R | colPivot = on" << nl;
+ const MatrixType A0(A);
+ QRMatrix<MatrixType> QRM
+ (
+ A0,
+ QRMatrix<MatrixType>::outputTypes::FULL_R,
+ QRMatrix<MatrixType>::colPivoting::TRUE
+ );
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ const labelList& PList(QRM.orderP());
+ const SMatrix P(QRM.P());
+ isEmpty(Q, "Q");
+ cross_check_QRMatrix(Aorig, PList, P, R);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# const A | FULL_QR | colPivot = on" << nl;
+ const MatrixType A0(A);
+ QRMatrix<MatrixType> QRM
+ (
+ A0,
+ QRMatrix<MatrixType>::outputTypes::FULL_QR,
+ QRMatrix<MatrixType>::colPivoting::TRUE
+ );
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ const labelList& PList(QRM.orderP());
+ const SMatrix P(QRM.P());
+ cross_check_QRMatrix(Aorig, Q, PList, P, R);
+ }
+ #endif
+
+ //
+ #if (0 | RUNALL)
+ {
+ Info<< "# REDUCED_R | OUT_OF_PLACE" << nl;
+ QRMatrix<MatrixType> QRM
+ (
+ QRMatrix<MatrixType>::outputTypes::REDUCED_R,
+ QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE
+ );
+ QRM.decompose(A);
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ isEmpty(Q, "Q");
+ cross_check_QRMatrix(R);
+
+ // check if full R and reduced R give the same answer
+ if (A.n() < A.m())
+ {
+ QRMatrix<MatrixType> fullQRM
+ (
+ QRMatrix<MatrixType>::outputTypes::FULL_QR,
+ QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE
+ );
+ fullQRM.decompose(A);
+ const MatrixType& fullR(fullQRM.R());
+ MatrixType fullR2(fullR.subMatrix(0,0,R.m()));
+ equal(fullR2, R, verbose);
+ }
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# REDUCED_R | IN_PLACE" << nl;
+ QRMatrix<MatrixType> QRM
+ (
+ QRMatrix<MatrixType>::outputTypes::REDUCED_R,
+ QRMatrix<MatrixType>::storeMethods::IN_PLACE
+ );
+ MatrixType A0(A);
+ QRM.decompose(A0);
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ isEmpty(Q, "Q");
+ isEmpty(R, "R");
+ cross_check_QRMatrix(A0);
+
+ // check if full R and reduced R give the same answer
+ if (A.n() < A.m())
+ {
+ QRMatrix<MatrixType> fullQRM
+ (
+ QRMatrix<MatrixType>::outputTypes::FULL_QR,
+ QRMatrix<MatrixType>::storeMethods::IN_PLACE
+ );
+ MatrixType A1(A);
+ fullQRM.decompose(A1);
+ MatrixType fullR(A1.subMatrix(0,0,A0.m()));
+ equal(fullR, A0, verbose);
+ }
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# REDUCED_R | OUT_OF_PLACE | colPivot = on" << nl;
+ QRMatrix<MatrixType> QRM
+ (
+ QRMatrix<MatrixType>::outputTypes::REDUCED_R,
+ QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE,
+ QRMatrix<MatrixType>::colPivoting::TRUE
+ );
+ QRM.decompose(A);
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ const labelList& PList(QRM.orderP());
+ const SMatrix P(QRM.P());
+ isEmpty(Q, "Q");
+ cross_check_QRMatrix(Aorig, PList, P, R);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# REDUCED_R | IN_PLACE | colPivot = on" << nl;
+ QRMatrix<MatrixType> QRM
+ (
+ QRMatrix<MatrixType>::outputTypes::REDUCED_R,
+ QRMatrix<MatrixType>::storeMethods::IN_PLACE,
+ QRMatrix<MatrixType>::colPivoting::TRUE
+ );
+ MatrixType A0(A);
+ QRM.decompose(A0);
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ const labelList& PList(QRM.orderP());
+ const SMatrix P(QRM.P());
+ isEmpty(Q, "Q");
+ isEmpty(R, "R");
+ cross_check_QRMatrix(Aorig, PList, P, A0);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# A | REDUCED_R | OUT_OF_PLACE | colPivot = off" << nl;
+ MatrixType A0(A);
+ QRMatrix<MatrixType> QRM
+ (
+ A0,
+ QRMatrix<MatrixType>::outputTypes::REDUCED_R,
+ QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE,
+ QRMatrix<MatrixType>::colPivoting::FALSE
+ );
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ isEmpty(Q, "Q");
+ cross_check_QRMatrix(R);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# A | REDUCED_R | IN_PLACE | colPivot = off" << nl;
+ MatrixType A0(A);
+ QRMatrix<MatrixType> QRM
+ (
+ A0,
+ QRMatrix<MatrixType>::outputTypes::REDUCED_R,
+ QRMatrix<MatrixType>::storeMethods::IN_PLACE,
+ QRMatrix<MatrixType>::colPivoting::FALSE
+ );
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ isEmpty(Q, "Q");
+ isEmpty(R, "R");
+ cross_check_QRMatrix(A0);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# A | REDUCED_R | OUT_OF_PLACE | colPivot = on" << nl;
+ MatrixType A0(A);
+ QRMatrix<MatrixType> QRM
+ (
+ A0,
+ QRMatrix<MatrixType>::outputTypes::REDUCED_R,
+ QRMatrix<MatrixType>::storeMethods::OUT_OF_PLACE,
+ QRMatrix<MatrixType>::colPivoting::TRUE
+ );
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ const labelList& PList(QRM.orderP());
+ const SMatrix P(QRM.P());
+ isEmpty(Q, "Q");
+ cross_check_QRMatrix(Aorig, PList, P, R);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# A | REDUCED_R | IN_PLACE | colPivot = on" << nl;
+ MatrixType A0(A);
+ QRMatrix<MatrixType> QRM
+ (
+ A0,
+ QRMatrix<MatrixType>::outputTypes::REDUCED_R,
+ QRMatrix<MatrixType>::storeMethods::IN_PLACE,
+ QRMatrix<MatrixType>::colPivoting::TRUE
+ );
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ const labelList& PList(QRM.orderP());
+ const SMatrix P(QRM.P());
+ isEmpty(Q, "Q");
+ isEmpty(R, "R");
+ cross_check_QRMatrix(Aorig, PList, P, A0);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# const A | REDUCED_R | colPivot = off" << nl;
+ const MatrixType A0(A);
+ QRMatrix<MatrixType> QRM
+ (
+ A0,
+ QRMatrix<MatrixType>::outputTypes::REDUCED_R,
+ QRMatrix<MatrixType>::colPivoting::FALSE
+ );
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ isEmpty(Q, "Q");
+ cross_check_QRMatrix(R);
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< "# const A | REDUCED_R | colPivot = on" << nl;
+ const MatrixType A0(A);
+ QRMatrix<MatrixType> QRM
+ (
+ A0,
+ QRMatrix<MatrixType>::outputTypes::REDUCED_R,
+ QRMatrix<MatrixType>::colPivoting::TRUE
+ );
+ const MatrixType& Q(QRM.Q());
+ const MatrixType& R(QRM.R());
+ const labelList& PList(QRM.orderP());
+ const SMatrix P(QRM.P());
+ isEmpty(Q, "Q");
+ cross_check_QRMatrix(Aorig, PList, P, R);
+ }
+ #endif
+}
+
+
+// Checks the direct solution of a linear system by QRMatrix
+template<class MatrixType>
+void verification_QRMatrixSolve
+(
+ MatrixType& A
+)
+{
+ typedef SquareMatrix<typename MatrixType::cmptType> SMatrix;
+
+ // Create working copies of matrix A
+ const MatrixType Aorig(A);
+
+ // Direct solution of a linear system by QR decomposition
+ #if (0 | RUNALL)
+ {
+ MatrixType A0(A);
+ QRMatrix<MatrixType> QRM
+ (
+ A0,
+ QRMatrix<MatrixType>::outputTypes::FULL_QR
+ );
+
+ Info<< nl << "# Solution of A*x = b with A = Q*R:" << nl;
+ const scalarField residual(A0.m(), 0);
+ const scalarField source(A0.m(), 1);
+
+ const scalarField x(QRM.solve(source));
+ const scalarField A0xMinusSource(A0*x - source);
+ const scalarField xMinusQRinvSource(x - QRM.inv()*source);
+ const SMatrix QRinvA0(QRM.inv()*A0);
+ const SMatrix identity(A0.m(), I);
+
+ forAll(residual, i)
+ {
+ isEqual(A0xMinusSource[i], residual[i], verbose);
+ isEqual(xMinusQRinvSource[i], residual[i], verbose);
+ }
+ equal(QRinvA0, identity, verbose);
+ }
+ #endif
+}
+
+
+// Checks the parallel tall-skinny QR decomposition
+template<class Type>
+void verification_tsqr
+(
+ Random& rndGen
+)
+{
+ typedef RectangularMatrix<Type> RMatrix;
+
+ // Size of the full matrix and its partitions
+ const label nColsSum = rndGen.position(1, 100);
+ const label qParts = rndGen.position(10, 30);
+ List<label> mRowsList(qParts, Zero);
+
+ label mRowsSum = 0;
+ for (label i = 0; i < qParts; ++i)
+ {
+ const label mRows = rndGen.position(nColsSum, 10*nColsSum);
+ mRowsList[i] = mRows;
+ mRowsSum += mRows;
+ }
+
+ # if 0
+ Info<< "mRowsSum = " << mRowsSum << nl;
+ Info<< "nColsSum = " << nColsSum << nl;
+ Info<< "mRowsList = " << mRowsList << nl;
+ #endif
+
+ if (mRowsSum < nColsSum)
+ {
+ FatalErrorInFunction
+ << "Tall skinny QR decomposition cannot be executed for matrices "
+ << "with number of columns (nCols) = " << nColsSum << " > "
+ << "than number of rows (mRows) = " << mRowsSum << "."
+ << "Yet nCols > mRows is allowed for a submatrix."
+ << exit(FatalError);
+ }
+
+ // Perform the QR decomposition for the full matrix
+ const RMatrix A(makeRandomMatrix<RMatrix>({mRowsSum, nColsSum}, rndGen));
+ QRMatrix<RMatrix> QRM
+ (
+ A,
+ QRMatrix<RMatrix>::outputTypes::REDUCED_R,
+ QRMatrix<RMatrix>::colPivoting::FALSE
+ );
+ const RMatrix& R = QRM.R();
+
+ RMatrix RMag(R);
+ for (auto& val : RMag)
+ {
+ val = mag(val);
+ }
+
+ // Partition the full matrix,
+ // and perform the QR decomposition on each partition
+ RMatrix subRs({qParts*nColsSum, nColsSum}, Zero);
+ label mRowsIndex = 0;
+ for (label i = 0; i < qParts; ++i)
+ {
+ RMatrix subA({mRowsList[i], nColsSum}, Zero);
+ subA = A.subMatrix(mRowsIndex, 0, mRowsList[i]);
+ mRowsIndex += mRowsList[i];
+
+ QRMatrix<RMatrix> subQRM
+ (
+ subA,
+ QRMatrix<RMatrix>::outputTypes::REDUCED_R,
+ QRMatrix<RMatrix>::colPivoting::FALSE
+ );
+ const RMatrix& subR = subQRM.R();
+
+ // Append subR
+ subRs.subMatrix(i*nColsSum, 0, nColsSum) = subR;
+ }
+
+ // Perform the tall-skinny QR decomposition
+ QRMatrix<RMatrix> TSQR
+ (
+ subRs,
+ QRMatrix<RMatrix>::outputTypes::REDUCED_R,
+ QRMatrix<RMatrix>::colPivoting::FALSE
+ );
+ const RMatrix& TSQRR = TSQR.R();
+
+ RMatrix TSQRRMag(TSQRR);
+ for (auto& val : TSQRRMag)
+ {
+ val = mag(val);
+ }
+
+ // Compare magnitude of R, since R is not unique up to the sign
+ equal(RMag, TSQRRMag, verbose);
+
+ #if 0
+ InfoNumPyFormat(R, "R");
+ InfoNumPyFormat(TSQRR, "TSQRR");
+ #endif
+}
+
+
+// Checks the back substitution function
+template<class Type>
+void verification_backSubstitution
+(
+ const SquareMatrix<Type>& A,
+ const RectangularMatrix<Type>& b
+)
+{
+ // Ax = b
+ typedef RectangularMatrix<Type> RMatrix;
+
+ QRMatrix<RMatrix> QRNull;
+ const RMatrix X(QRNull.backSubstitution(A, b));
+
+ #if 1
+ {
+ InfoNumPyFormat(A, "A");
+ InfoNumPyFormat(b, "b");
+ InfoNumPyFormat(X, "x");
+ }
+ #endif
+
+ #if (0 | RUNALL)
+ {
+ Info<< nl << "# A*X = b:" << nl;
+ const RMatrix AX(A*X);
+ equal(AX, b, verbose, 1e-3, 1e-3);
+ }
+ #endif
+}
+
+
+// Checks the Moore-Penrose pseudo-inverse function
+template<class MatrixType>
+void verification_pinv
+(
+ const MatrixType& A
+)
+{
+ Info<< "### Moore-Penrose pseudo-inverse: pinv()" << nl << nl;
+ const MatrixType APlus(pinv(A));
+
+ #if 0
+ InfoNumPyFormat(A, "A");
+ InfoNumPyFormat(APlus, "A^+");
+ #endif
+
+ if (0 | RUNALL)
+ {
+ Info<< nl << "# (A^+)^+ = A:" << nl;
+ const MatrixType APlusPlus(pinv(APlus));
+ equal(APlusPlus, A, verbose);
+ }
+
+ // The four Moore-Penrose conditions
+ // w.wiki/5RL (Retrieved: 29-06-19)
+ if (0 | RUNALL)
+ {
+ Info<< nl << "# A*(A^+)*A = A:" << nl;
+ const MatrixType AAPlusA((A*APlus)*A);
+ equal(AAPlusA, A, verbose);
+ }
+
+ if (0 | RUNALL)
+ {
+ Info<< nl << "# (A^+)*A*(A^+) = A:" << nl;
+ const MatrixType APlusAAPlus((APlus*A)*APlus);
+ equal(APlusAAPlus, APlus, verbose);
+ }
+
+ if (0 | RUNALL)
+ {
+ Info<< nl << "# (A*(A^+)).T() = A*(A^+):" << nl;
+ const MatrixType AAPlus(A*APlus);
+ const MatrixType AAPlusT(AAPlus.T());
+ equal(AAPlusT, AAPlus, verbose);
+ }
+
+ if (0 | RUNALL)
+ {
+ Info<< nl << "# ((A^+)*A = ((A^+)*A).T():" << nl;
+ const MatrixType APlusA(APlus*A);
+ const MatrixType APlusAT(APlusA.T());
+ equal(APlusA, APlusAT, verbose);
+ }
+}
+
+
+// * * * * * * * * * * * * * * * Main Program * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+ typedef SquareMatrix<scalar> SMatrix;
+ typedef RectangularMatrix<scalar> RMatrix;
+ typedef SquareMatrix<complex> SCMatrix;
+ typedef RectangularMatrix<complex> RCMatrix;
+
+ Info<< setprecision(15);
+ Random rndGen(1234);
+ label numberOfTests = 10;
+
+ Info<< "### QR Decomposition:" << nl << nl;
+
+ // SquareMatrix<scalar>
+ #if (0 | RUNALL)
+ {
+ horizontalLine();
+
+ Info<< "## " << numberOfTests << " SquareMatrix<scalar> A:" << nl;
+
+ for (label i = 0; i < numberOfTests; ++i)
+ {
+ const label mRows = rndGen.position(1, 50);
+ Info<< nl << nl << "# Random A with random mRows = " << mRows << nl;
+
+ SMatrix A(makeRandomMatrix<SMatrix>({mRows, mRows}, rndGen));
+ #if (0 | RUNALL)
+ {
+ const SMatrix B(A);
+ verification_pinv(B);
+ }
+ #endif
+ verification_QRMatrix(A);
+ verification_QRMatrixSolve(A);
+
+ }
+
+ horizontalLine();
+ }
+ #endif
+
+
+ // RectangularMatrix<scalar>
+ #if (0 | RUNALL)
+ {
+ horizontalLine();
+
+ Info<< "## " << numberOfTests << " RectangularMatrix<scalar> A:" << nl;
+
+ for (label i = 0; i < numberOfTests; ++i)
+ {
+ const label mRows = rndGen.position(1, 50);
+ const label nCols = rndGen.position(1, 50);
+ Info<< nl << nl << "# Random matrix A with"
+ << " random mRows = " << mRows
+ << " random nCols = " << nCols << nl;
+
+ RMatrix A(makeRandomMatrix<RMatrix>({mRows, nCols}, rndGen));
+
+ #if (0 | RUNALL)
+ {
+ const RMatrix B(A);
+ verification_pinv(B);
+ }
+ #endif
+
+ verification_QRMatrix(A);
+ }
+
+ horizontalLine();
+ }
+ #endif
+
+
+ // SquareMatrix<complex(scalar, 0.0)>
+ #if (0 | RUNALL)
+ {
+ horizontalLine();
+
+ Info<< "## " << numberOfTests << " SquareMatrix<complex, 0> A:" << nl;
+
+ for (label i = 0; i < numberOfTests; ++i)
+ {
+ const label mRows = rndGen.position(1, 50);
+ Info<< nl << nl << "# Random A with random mRows = " << mRows << nl;
+
+ SCMatrix A({mRows, mRows}, complex(0, 0));
+
+ for (auto& val : A)
+ {
+ val.Re() = rndGen.GaussNormal<scalar>();
+ }
+
+ verification_QRMatrix(A);
+ }
+
+ horizontalLine();
+ }
+ #endif
+
+
+ // SquareMatrix<complex(scalar, scalar)>
+ #if (0 | RUNALL)
+ {
+ horizontalLine();
+
+ Info<< "## " << numberOfTests << " SquareMatrix<complex> A:" << nl;
+
+ for (label i = 0; i < numberOfTests; ++i)
+ {
+ const label mRows = rndGen.position(1, 50);
+ Info<< nl << nl << "# Random A with random mRows = " << mRows << nl;
+
+ SCMatrix A(makeRandomMatrix<SCMatrix>({mRows, mRows}, rndGen));
+
+ #if (0 | RUNALL)
+ {
+ const SCMatrix B(A);
+ verification_pinv(B);
+ }
+ #endif
+
+ verification_QRMatrix(A);
+ }
+
+ horizontalLine();
+ }
+ #endif
+
+
+ // RectangularMatrix<complex(scalar, scalar)>
+ #if (0 | RUNALL)
+ {
+ horizontalLine();
+
+ Info<< "## " << numberOfTests << " RectangularMatrix<complex> A:" << nl;
+
+ for (label i = 0; i < numberOfTests; ++i)
+ {
+ const label mRows = rndGen.position(1, 50);
+ const label nCols = rndGen.position(1, 50);
+ Info<< nl << nl << "# Random matrix A with"
+ << " random mRows = " << mRows
+ << " random nCols = " << nCols << nl;
+
+ RCMatrix A(makeRandomMatrix<RCMatrix>({mRows, nCols}, rndGen));
+
+ #if (0 | RUNALL)
+ {
+ const RCMatrix B(A);
+ verification_pinv(B);
+ }
+ #endif
+
+ verification_QRMatrix(A);
+ }
+
+ horizontalLine();
+ }
+ #endif
+
+
+ #if (0 | RUNALL)
+ {
+ horizontalLine();
+
+ Info<< "### Parallel direct tall-skinny QR decomposition:" << nl;
+ /*
+ Parallel direct tall-skinny QR decomposition:
+ Benson, A. R., Gleich, D. F., & Demmel, J. (2013).
+ Direct QR factorizations for tall-and-skinny matrices in MapReduce
+ architectures.
+ 2013 IEEE International Conference on Big Data.
+ DOI:10.1109/bigdata.2013.6691583
+
+ Demmel, J., Grigori, L., Hoemmen, M., & Langou, J. (2012).
+ Communication-optimal parallel and sequential QR and LU
+ factorizations.
+ SIAM Journal on Scientific Computing, 34(1), A206-A239.
+ DOI:10.1137/080731992
+ */
+ // (Benson et al. 2013, Fig. 5) & (Demmel et al. 2012, Fig. 2)
+
+ Info<< "## " << numberOfTests << " <scalar> tests:" << nl;
+ for (label i = 0; i < numberOfTests; ++i)
+ {
+ verification_tsqr<scalar>(rndGen);
+ }
+
+ Info<< nl << "## " << numberOfTests << " <complex> tests:" << nl;
+ for (label i = 0; i < numberOfTests; ++i)
+ {
+ verification_tsqr<complex>(rndGen);
+ }
+
+ horizontalLine();
+ }
+ #endif
+
+
+ #if (0 | RUNALL)
+ {
+ Info<< "### Back substitution:" << nl << nl;
+
+ numberOfTests = 100;
+
+ // <scalar>
+ #if (0 | RUNALL)
+ {
+ horizontalLine();
+
+ Info<< "## " << numberOfTests << " <scalar>:" << nl;
+
+ for (label i = 0; i < numberOfTests; ++i)
+ {
+ const label mRows = rndGen.position(20, 50);
+ const label pCols = rndGen.position(20, 50);
+ Info<< nl << nl << "# Random square matrix A with"
+ << " random mRows = " << mRows
+ << " random nCols = " << mRows << nl;
+
+ SMatrix A(makeRandomMatrix<SMatrix>({mRows, mRows}, rndGen));
+
+ // Zeroise the lower diagonal,
+ // so that A becomes an upper-triangular matrix
+ for (label i = 0; i < mRows; ++i)
+ {
+ for (label j = 0; j < i; ++j)
+ {
+ A(i, j) = 0.0;
+ }
+ }
+
+ // det(triangularA) = prod(diag(triangularA))
+ scalar det = 1;
+ const List<scalar> diag0(A.diag());
+ for (const auto& val : diag0)
+ {
+ det *= val;
+ }
+
+ if (1e-8 < det)
+ {
+ Info<< "# Random matrix b with"
+ << " random mRows = " << mRows
+ << " random nCols = " << pCols << nl;
+
+ const RMatrix b
+ (
+ makeRandomMatrix<RMatrix>({mRows, pCols}, rndGen)
+ );
+ Info<< "det(A) = " << det << nl;
+ verification_backSubstitution(A, b);
+ }
+ else
+ {
+ Info<< "Random matrix A is a nearly-singular matrix."
+ << " Skipping back-substitution" << nl;
+ }
+ }
+
+ horizontalLine();
+ }
+ #endif
+
+ // <complex>
+ #if (0 | RUNALL)
+ {
+ horizontalLine();
+
+ Info<< "## " << numberOfTests << " <complex>:" << nl;
+
+ for (label i = 0; i < numberOfTests; ++i)
+ {
+ const label mRows = rndGen.position(20, 50);
+ const label pCols = rndGen.position(20, 50);
+ Info<< nl << nl << "# Random square matrix A with"
+ << " random mRows = " << mRows
+ << " random nCols = " << mRows << nl;
+
+ SCMatrix A(makeRandomMatrix<SCMatrix>({mRows, mRows}, rndGen));
+
+ // Zeroise the lower diagonal,
+ // so that A becomes an upper-triangular matrix
+ for (label i = 0; i < mRows; ++i)
+ {
+ for (label j = 0; j < i; ++j)
+ {
+ A(i, j) = pTraits<complex>::zero;
+ }
+ }
+
+ // det(triangularA) = prod(diag(triangularA))
+ complex det = pTraits<complex>::one;
+ const List<complex> diag0(A.diag());
+ for (const auto& val : diag0)
+ {
+ det *= val;
+ }
+
+ if (1e-8 < mag(det))
+ {
+ Info<< "# Random matrix b with"
+ << " random mRows = " << mRows
+ << " random nCols = " << pCols << nl;
+
+ const RCMatrix b
+ (
+ makeRandomMatrix<RCMatrix>({mRows, pCols}, rndGen)
+ );
+ Info<< "det(A) = " << det << nl;
+ verification_backSubstitution(A, b);
+ }
+ else
+ {
+ Info<< "Random matrix A is a nearly-singular matrix."
+ << " Skipping back-substitution" << nl;
+ }
+ }
+
+ horizontalLine();
+ }
+ #endif
+ }
+ #endif
+
+
+ Info<< nl << "End" << nl;
+
+ return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/Matrix/MatrixTools.C b/src/OpenFOAM/matrices/Matrix/MatrixTools.C
index cf9cb1813072ad487e31c564e4729cac557b107f..30dbce057147b3a16bb106bf0c25f44162e35e78 100644
--- a/src/OpenFOAM/matrices/Matrix/MatrixTools.C
+++ b/src/OpenFOAM/matrices/Matrix/MatrixTools.C
@@ -35,6 +35,7 @@ bool Foam::MatrixTools::equal
const Matrix<Form1, Type>& A,
const Matrix<Form2, Type>& B,
const bool verbose,
+ const label lenDiffs,
const scalar relTol,
const scalar absTol
)
@@ -54,6 +55,8 @@ bool Foam::MatrixTools::equal
auto iter1 = A.cbegin();
auto iter2 = B.cbegin();
+ label j = 0;
+
for (label i = 0; i < len; ++i)
{
if ((absTol + relTol*mag(*iter2)) < Foam::mag(*iter1 - *iter2))
@@ -63,8 +66,15 @@ bool Foam::MatrixTools::equal
Info<< "Matrix element " << i
<< " differs beyond tolerance: "
<< *iter1 << " vs " << *iter2 << nl;
+ ++j;
+ }
+ if (lenDiffs < j)
+ {
+ Info<< "Number of different elements exceeds = " << lenDiffs
+ << " Ceasing comparisons for the remaining of elements."
+ << nl;
+ return false;
}
- return false;
}
++iter1;
diff --git a/src/OpenFOAM/matrices/Matrix/MatrixTools.H b/src/OpenFOAM/matrices/Matrix/MatrixTools.H
index 89e74388570d30e89d54ac8a2c4bb596ac397190..58bdac3ef482662637ae340201dfd1568d9cab50 100644
--- a/src/OpenFOAM/matrices/Matrix/MatrixTools.H
+++ b/src/OpenFOAM/matrices/Matrix/MatrixTools.H
@@ -62,6 +62,7 @@ bool equal
const Matrix<Form1, Type>& A,
const Matrix<Form2, Type>& B,
const bool verbose = false,
+ const label lenDiffs = 10,
const scalar relTol = 1e-5,
const scalar absTol = 1e-8
);
diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
index f4e5fbb996b90590c0d9b694cfa9f58c75147f2d..4394700e0e282b7b13d302ef78f08bbeff47ea8b 100644
--- a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
+++ b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
@@ -28,132 +28,99 @@ License
#include "QRMatrix.H"
-// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class MatrixType>
-Foam::QRMatrix<MatrixType>::QRMatrix(const MatrixType& mat)
+void Foam::QRMatrix<MatrixType>::qr
+(
+ MatrixType& A
+)
{
- decompose(mat);
-}
+ const label nIter = min(A.m() - 1, A.n());
+ // Loop through all subcolumns of which the diagonal elem is the first elem
+ for (label k = 0; k < nIter; ++k)
+ {
+ const RMatrix u(householderReflector(A.subColumn(k, k)));
+
+ applyHouseholder(A, u, k);
+ }
+
+ if (outputType_ == outputTypes::REDUCED_R)
+ {
+ A.resize(A.n(), A.n());
+ }
+}
-// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
template<class MatrixType>
-void Foam::QRMatrix<MatrixType>::decompose(const MatrixType& mat)
+void Foam::QRMatrix<MatrixType>::qrPivot
+(
+ MatrixType& A
+)
{
- const label m = mat.m();
- const label n = mat.n();
-
- // Initialize the R-matrix to M
- R_ = mat;
+ const label nCols = A.n();
+ const label nIter = min(A.m() - 1, A.n());
- // Initialize the Q-matrix to I
- Q_.setSize(m);
- Q_ = I;
+ // Initialise permutation vector, and column norm vector
+ P_ = identity(nCols);
- // Pre-allocate temporary storage for the Householder steps
- QMatrixType Qk(m);
- QMatrixType Rk(m);
- Field<cmptType> uk(m);
+ // Initialise vector norms of each column of A
+ List<scalar> colNorms(nCols);
+ for (label k = 0; k < nCols; ++k)
+ {
+ colNorms[k] = A.columnNorm(k, true);
+ }
- for (label k=0; k<n; k++)
+ // Loop through all subcolumns of which the diagonal elem is the first elem
+ for (label k = 0; k < nIter; ++k)
{
- // alpha = -|column k of Rk|
- cmptType alpha = Zero;
- for (label bi=k; bi<m; bi++)
+ const labelRange colRange(k, nCols);
+ const SubList<scalar> subColNorms(colNorms, colRange);
+
+ // Column pivoting
+ const label maxColNormi =
+ std::distance
+ (
+ subColNorms.cbegin(),
+ std::max_element(subColNorms.cbegin(), subColNorms.cend())
+ );
+
+ // Swap R_, P_ and colNorms_ according to pivot column if the current
+ // column is not the max norm column by avoiding any column swap where
+ // the leading elem is very small
+ if (maxColNormi != 0 && SMALL < mag(A(k, k + maxColNormi)))
{
- alpha += sqr(R_(bi, k));
- }
- alpha = sqrt(alpha);
+ const RMatrix R1(A.subColumn(k));
+ const RMatrix R2(A.subColumn(maxColNormi + k));
+ A.subColumn(k) = R2;
+ A.subColumn(maxColNormi + k) = R1;
- if (R_(k, k) > 0)
- {
- alpha = -alpha;
+ Swap(P_[k], P_[maxColNormi + k]);
+ Swap(colNorms[k], colNorms[maxColNormi + k]);
}
- // uk = column k of Rk - alpha*ek
- // rSumSqrUk = 2/sum(sqr(uk))
- uk[k] = R_(k, k) - alpha;
- cmptType rSumSqrUk = sqr(uk[k]);
- for (label bi=k+1; bi<m; bi++)
{
- uk[bi] = R_(bi, k);
- rSumSqrUk += sqr(uk[bi]);
- }
- rSumSqrUk = 2/rSumSqrUk;
+ const RMatrix u(householderReflector(A.subColumn(k, k)));
- // Qk = I - 2*u*uT/sum(sqr(uk))
- for (label bi=k; bi<m; bi++)
- {
- for (label bj=k; bj<m; bj++)
- {
- Qk(bi, bj) = -rSumSqrUk*uk[bi]*uk[bj];
- }
- Qk(bi, bi) += 1;
+ applyHouseholder(A, u, k);
}
- // Rk = Qk*R
- for (label bi=k; bi<m; bi++)
+ // Update norms
+ if (k < nIter - 1)
{
- for (label bk=k; bk<n; bk++)
+ label q = k + 1;
+ for (const auto& val : RMatrix(A.subRow(k, q)))
{
- Rk(bi, bk) = Zero;
- for (label bj=k; bj<n; bj++)
- {
- Rk(bi, bk) += Qk(bi, bj)*R_(bj, bk);
- }
- }
- }
-
- // Ensure diagonal is positive
- if (R_(k, k) < 0)
- {
- // R = -Rk
- // Q = -Q
- for (label bi=k; bi<m; bi++)
- {
- for (label bj=k; bj<n; bj++)
- {
- R_(bi, bj) = -Rk(bi, bj);
- }
- for (label bj=k; bj<m; bj++)
- {
- Q_(bi, bj) = -Q_(bi, bj);
- }
- }
- }
- else
- {
- // R = Rk
- for (label bi=k; bi<m; bi++)
- {
- for (label bj=k; bj<n; bj++)
- {
- R_(bi, bj) = Rk(bi, bj);
- }
+ colNorms[q] -= magSqr(val);
+ ++q;
}
}
+ }
- // Q = Q*Qk (using Rk as temporary storage)
- for (label bi=0; bi<m; bi++)
- {
- for (label bk=k; bk<m; bk++)
- {
- Rk(bi, bk) = Zero;
- for (label bj=k; bj<m; bj++)
- {
- Rk(bi, bk) += Q_(bi, bj)*Qk(bj, bk);
- }
- }
- }
- for (label bi=0; bi<m; bi++)
- {
- for (label bj=k; bj<n; bj++)
- {
- Q_(bi, bj) = Rk(bi, bj);
- }
- }
+ if (outputType_ == outputTypes::REDUCED_R)
+ {
+ A.resize(A.n(), A.n());
}
}
@@ -167,7 +134,7 @@ void Foam::QRMatrix<MatrixType>::solvex
{
const label n = R_.n();
- for (label i = n - 1; i >= 0; --i)
+ for (label i = n - 1; 0 <= i; --i)
{
cmptType sum = x[i];
@@ -199,7 +166,7 @@ void Foam::QRMatrix<MatrixType>::solveImpl
// Assert (&x != &source) ?
const label m = Q_.m();
- // x = Q_.T()*source; ie, Q_.Tmul(source)
+ // x = Q_.T()*source; i.e., Q_.Tmul(source)
for (label i = 0; i < m; ++i)
{
x[i] = 0;
@@ -213,6 +180,179 @@ void Foam::QRMatrix<MatrixType>::solveImpl
}
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+template<class MatrixType>
+Foam::QRMatrix<MatrixType>::QRMatrix()
+:
+ outputType_(),
+ storeMethod_(),
+ colPivot_()
+{}
+
+
+template<class MatrixType>
+Foam::QRMatrix<MatrixType>::QRMatrix
+(
+ const outputTypes outputType,
+ const storeMethods storeMethod,
+ const colPivoting colPivot
+)
+:
+ outputType_(outputType),
+ storeMethod_(storeMethod),
+ colPivot_(colPivot),
+ Q_(),
+ R_(),
+ P_()
+{}
+
+
+template<class MatrixType>
+Foam::QRMatrix<MatrixType>::QRMatrix
+(
+ MatrixType& A,
+ const outputTypes outputType,
+ const storeMethods storeMethod,
+ const colPivoting colPivot
+)
+:
+ outputType_(outputType),
+ storeMethod_(storeMethod),
+ colPivot_(colPivot),
+ Q_(),
+ R_(),
+ P_()
+{
+ decompose(A);
+}
+
+
+template<class MatrixType>
+Foam::QRMatrix<MatrixType>::QRMatrix
+(
+ const MatrixType& A,
+ const outputTypes outputType,
+ const colPivoting colPivot
+)
+:
+ outputType_(outputType),
+ storeMethod_(storeMethods::OUT_OF_PLACE),
+ colPivot_(colPivot),
+ Q_(),
+ R_(),
+ P_()
+{
+ decompose(A);
+}
+
+
+// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
+
+template<class MatrixType>
+void Foam::QRMatrix<MatrixType>::decompose
+(
+ MatrixType& A
+)
+{
+ // Check whether settings and input are consistent for reduced QRMatrix
+ if (A.m() <= A.n() && outputType_ == outputTypes::REDUCED_R)
+ {
+ outputType_ = outputTypes::FULL_R;
+
+ #if FULLDEBUG
+ WarningInFunction
+ << "Reduced QR decomposition is by definition limited to matrices"
+ << " wherein rows > columns, thus computing FULL decomposition."
+ << nl;
+ #endif
+ }
+
+ // Allocate resources for Q_, if need be
+ if (outputType_ == outputTypes::FULL_QR)
+ {
+ Q_ = MatrixType({A.m(), A.m()}, I);
+ }
+
+ // Allocate resources for R_ and execute decomposition
+ switch (storeMethod_)
+ {
+ case storeMethods::IN_PLACE:
+ {
+ if (colPivot_)
+ {
+ qrPivot(A);
+ }
+ else
+ {
+ qr(A);
+ }
+ break;
+ }
+
+ case storeMethods::OUT_OF_PLACE:
+ {
+ R_ = A;
+
+ if (colPivot_)
+ {
+ qrPivot(R_);
+ }
+ else
+ {
+ qr(R_);
+ }
+ break;
+ }
+ }
+}
+
+
+template<class MatrixType>
+void Foam::QRMatrix<MatrixType>::decompose
+(
+ const MatrixType& A
+)
+{
+ if (storeMethod_ == storeMethods::IN_PLACE)
+ {
+ WarningInFunction
+ << "const type qualifier invalidates storeMethods::IN_PLACE." << nl;
+ }
+
+ // Check whether settings and input are consistent for reduced QRMatrix
+ if (A.m() <= A.n() && outputType_ == outputTypes::REDUCED_R)
+ {
+ outputType_ = outputTypes::FULL_R;
+
+ #if FULLDEBUG
+ WarningInFunction
+ << "Reduced QR decomposition is by definition limited to matrices"
+ << " wherein rows > columns, thus computing FULL decomposition."
+ << nl;
+ #endif
+ }
+
+ // Allocate resources for Q_, if need be
+ if (outputType_ == outputTypes::FULL_QR)
+ {
+ Q_ = MatrixType({A.m(), A.m()}, I);
+ }
+
+ // Allocate resources for R_ and execute decomposition
+ R_ = A;
+
+ if (colPivot_)
+ {
+ qrPivot(R_);
+ }
+ else
+ {
+ qr(R_);
+ }
+}
+
+
template<class MatrixType>
void Foam::QRMatrix<MatrixType>::solve
(
@@ -268,13 +408,13 @@ Foam::QRMatrix<MatrixType>::solve
template<class MatrixType>
-typename Foam::QRMatrix<MatrixType>::QMatrixType
+typename Foam::QRMatrix<MatrixType>::SMatrix
Foam::QRMatrix<MatrixType>::inv() const
{
const label m = Q_.m();
Field<cmptType> x(m);
- QMatrixType inv(m);
+ SMatrix inv(m);
for (label i = 0; i < m; ++i)
{
@@ -290,4 +430,168 @@ Foam::QRMatrix<MatrixType>::inv() const
}
+template<class MatrixType>
+Foam::RectangularMatrix<typename MatrixType::cmptType>
+Foam::QRMatrix<MatrixType>::backSubstitution
+(
+ const SMatrix& A,
+ const RMatrix& rhs
+)
+{
+ const label mRows = A.m();
+ const label nCols = A.n();
+ const label pCols = rhs.n();
+
+ #ifdef FULLDEBUG
+ {
+ const label qRows = rhs.m();
+ if (mRows != qRows)
+ {
+ FatalErrorInFunction
+ << "Linear system is not solvable since the number of rows of"
+ << "A and rhs are not equal:" << tab << mRows << "vs." << qRows
+ << abort(FatalError);
+ }
+
+ const List<cmptType> diag(A.diag());
+ for (const auto& val : diag)
+ {
+ if (mag(val) < SMALL)
+ {
+ WarningInFunction
+ << "SMALL-valued diagonal elem in back-substitution." << nl;
+ }
+ }
+ }
+ #endif
+
+ RMatrix X({nCols, pCols}, Zero);
+
+ for (label i = 0; i < pCols; ++i)
+ {
+ for (label j = mRows - 1; -1 < j; --j)
+ {
+ cmptType alpha(rhs(j, i));
+
+ for (label k = j + 1; k < mRows; ++k)
+ {
+ alpha -= X(k, i)*A(j, k);
+ }
+
+ X(j, i) = alpha/A(j, j);
+ }
+ }
+
+ return X;
+}
+
+
+// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
+
+template<class MatrixType>
+MatrixType Foam::pinv
+(
+ const MatrixType& A,
+ const scalar tol
+)
+{
+ typedef typename MatrixType::cmptType cmptType;
+
+ #if FULLDEBUG
+ if (A.empty())
+ {
+ FatalErrorInFunction
+ << "Empty matrix." << abort(FatalError);
+ }
+
+ // Check if A is full-rank
+ #endif
+
+ QRMatrix<MatrixType> QRM
+ (
+ A,
+ QRMatrix<MatrixType>::outputTypes::FULL_QR,
+ QRMatrix<MatrixType>::colPivoting::TRUE
+ );
+ const MatrixType& R(QRM.R());
+ const MatrixType& Q(QRM.Q());
+
+ // R1 (KP:p. 648)
+ // Find the first diagonal element index with (almost) zero value in R
+ label firstZeroElemi = 0;
+ {
+ const List<cmptType> diag(R.diag());
+
+ auto lessT = [&](const cmptType& x) { return mag(x) < tol; };
+
+ firstZeroElemi =
+ std::distance
+ (
+ diag.cbegin(),
+ std::find_if(diag.cbegin(), diag.cend(), lessT)
+ );
+ }
+
+ if (firstZeroElemi == 0)
+ {
+ FatalErrorInFunction
+ << "The largest (magnitude) diagonal element is (almost) zero."
+ << abort(FatalError);
+ }
+
+ // Exclude the first (almost) zero diagonal row and the rows below
+ // since R diagonal is already descending due to the QR column pivoting
+ const RectangularMatrix<cmptType> R1(R.subMatrix(0, 0, firstZeroElemi));
+
+ // R2 (KP:p. 648)
+ if (R1.n() < R1.m())
+ {
+ const SquareMatrix<cmptType> C(R1&R1);
+
+ QRMatrix<SquareMatrix<cmptType>> QRSolve
+ (
+ C,
+ QRMatrix<SquareMatrix<cmptType>>::outputTypes::FULL_QR
+ );
+
+ RectangularMatrix<cmptType> R2
+ (
+ QRSolve.backSubstitution
+ (
+ QRSolve.R(),
+ QRSolve.Q() & (R1.T())
+ )
+ );
+
+ // R3 (KP:p. 648)
+ R2.resize(R.m(), R.n());
+
+ return MatrixType((QRM.P()^R2)^Q);
+ }
+ else
+ {
+ const SquareMatrix<cmptType> C(R1^R1);
+
+ QRMatrix<SquareMatrix<cmptType>> QRSolve
+ (
+ C,
+ QRMatrix<SquareMatrix<cmptType>>::outputTypes::FULL_QR
+ );
+
+ RectangularMatrix<cmptType> R2
+ (
+ QRSolve.backSubstitution
+ (
+ QRSolve.R(),
+ QRSolve.Q() & R1
+ )
+ );
+
+ // R3
+ R2.resize(R.m(), R.n());
+
+ return MatrixType((QRM.P()^R2)^Q);
+ }
+}
+
// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H
index 9058cd027ea6c234936344d49efd3f77ddbd584b..415fbba853c3cc73527a4974df120212d2989fe4 100644
--- a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H
+++ b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H
@@ -28,19 +28,117 @@ Class
Foam::QRMatrix
Description
- Class templated on matrix type to perform the QR decomposition using
- Householder reflections on a square or rectangular matrix.
+ QRMatrix (i.e. QR decomposition, QR factorisation or orthogonal-triangular
+ decomposition) decomposes a scalar/complex matrix \c A into the following
+ matrix product:
+
+ \verbatim
+ A = Q*R,
+ \endverbatim
+
+ where
+ \c Q is a unitary similarity matrix,
+ \c R is an upper triangular matrix.
+
+ Reference:
+ \verbatim
+ QR decomposition:
+ mathworld.wolfram.com/QRDecomposition.html (Retrieved:17-06-19)
+ w.wiki/52X (Retrieved: 17-06-19)
+
+ QR decomposition with Householder reflector (tag:M):
+ Monahan, J. F. (2011).
+ Numerical methods of statistics.
+ Cambridge: Cambridge University Press.
+ DOI:10.1017/CBO9780511977176
+
+ QR decomposition with column pivoting (tag:QSB):
+ Quintana-OrtÃ, G., Sun, X., & Bischof, C. H. (1998).
+ A BLAS-3 version of the QR factorization with column pivoting.
+ SIAM Journal on Scientific Computing, 19(5), 1486-1494.
+ DOI:10.1137/S1064827595296732
+
+ Moore-Penrose inverse algorithm (tags:KP; KPP):
+ Katsikis, V. N., & Pappas, D. (2008).
+ Fast computing of the Moore-Penrose inverse matrix.
+ Electronic Journal of Linear Algebra, 17(1), 637-650.
+ DOI:10.13001/1081-3810.1287
+
+ Katsikis, V. N., Pappas, D., & Petralias, A. (2011).
+ An improved method for the computation of
+ the Moore–Penrose inverse matrix.
+ Applied Mathematics and Computation, 217(23), 9828-9834.
+ DOI:10.1016/j.amc.2011.04.080
+
+ Tolerance for the Moore-Penrose inverse algorithm (tag:TA):
+ Toutounian, F., & Ataei, A. (2009).
+ A new method for computing Moore–Penrose inverse matrices.
+ Journal of Computational and applied Mathematics, 228(1), 412-417.
+ DOI:10.1016/j.cam.2008.10.008
+
+ Operands of QR decomposition:
+ mathworld.wolfram.com/UnitaryMatrix.html (Retrieved:13-06-19)
+ mathworld.wolfram.com/UpperTriangularMatrix.html (Retrieved:16-06-19)
+ mathworld.wolfram.com/PermutationMatrix.html (Retrieved:20-06-19)
+
+ Back substitution:
+ mathworld.wolfram.com/GaussianElimination.html (Retrieved:15-06-19)
+ \endverbatim
+
+Usage
+ Input types:
+ - \c A can be a \c SquareMatrix<Type> or \c RectangularMatrix<Type>
+
+ Output types:
+ - \c Q is always of the type of the matrix \c A
+ - \c R is always of the type of the matrix \c A
+
+ Options for the output forms of \c QRMatrix (for an (m-by-n) input matrix
+ \c A with k = min(m, n)):
+ - outputTypes::FULL_R: computes only \c R (m-by-n)
+ - outputTypes::FULL_QR: computes both \c R and \c Q (m-by-m)
+ - outputTypes::REDUCED_R: computes only reduced \c R (k-by-n)
+
+ Options where to store \c R:
+ - storeMethods::IN_PLACE: replaces input matrix content with \c R
+ - storeMethods::OUT_OF_PLACE: creates new object of \c R
+
+ Options for the computation of column pivoting:
+ - colPivoting::FALSE: switches off column pivoting
+ - colPivoting::TRUE: switches on column pivoting
+
+ Direct solution of linear systems A x = b is possible by solve() alongside
+ the following limitations:
+ - \c A = a scalar square matrix
+ - output type = outputTypes::FULL_QR
+ - store method = storeMethods::IN_PLACE
+
+Notes
+ - QR decomposition is not unique if \c R is not positive diagonal \c R.
+ - The option combination:
+ - outputTypes::REDUCED_R
+ - storeMethods::IN_PLACE
+ will not modify the rows of input matrix \c A after its nth row.
+ - Both FULL_R and REDUCED_R QR decompositions execute the same number of
+ operations. Yet REDUCED_R QR decomposition returns only the first n rows
+ of \c R if m > n for an input m-by-n matrix \c A.
+ - For m <= n, FULL_R and REDUCED_R will produce the same matrices.
+
+See also
+ Test-QRMatrix.C
SourceFiles
- QRMatrixI.H
QRMatrix.C
+ QRMatrixI.H
\*---------------------------------------------------------------------------*/
#ifndef QRMatrix_H
#define QRMatrix_H
+#include "RectangularMatrix.H"
#include "SquareMatrix.H"
+#include "complex.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -48,32 +146,90 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
- Class QRMatrix Declaration
+ Class QRMatrix Declaration
\*---------------------------------------------------------------------------*/
template<class MatrixType>
class QRMatrix
{
+
public:
typedef typename MatrixType::cmptType cmptType;
- typedef SquareMatrix<cmptType> QMatrixType;
- typedef MatrixType RMatrixType;
+ typedef SquareMatrix<cmptType> SMatrix;
+ typedef RectangularMatrix<cmptType> RMatrix;
+
+ //- Options for the output matrix forms of QRMatrix
+ enum outputTypes : uint8_t
+ {
+ FULL_R = 1, //!< computes only \c R
+ FULL_QR = 2, //!< computes both \c R and \c Q
+ REDUCED_R = 3 //!< computes only reduced \c R
+ };
+
+ //- Options where to store R
+ enum storeMethods : uint8_t
+ {
+ IN_PLACE = 1, //!< replaces input matrix content with \c R
+ OUT_OF_PLACE = 2 //!< creates new object of \c R
+ };
+
+ //- Options for the computation of column pivoting
+ enum colPivoting : bool
+ {
+ FALSE = false, //!< switches off column pivoting
+ TRUE = true //!< switches on column pivoting
+ };
private:
// Private Data
- //- The Q-matrix
- QMatrixType Q_;
+ //- Selected option for the output matrix forms of QRMatrix
+ outputTypes outputType_;
+
+ //- Selected option where to store R
+ const storeMethods storeMethod_;
+
+ //- Selected option for the computation of column pivoting
+ const colPivoting colPivot_;
- //- The R-matrix
- RMatrixType R_;
+ //- Unitary similarity matrix
+ MatrixType Q_;
+
+ //- Upper triangular matrix
+ MatrixType R_;
+
+ //- Permutation list (if column-pivoting is on)
+ labelList P_;
// Private Member Functions
+ //- Compute the QR decomposition without the column pivoting
+ void qr
+ (
+ MatrixType& A
+ );
+
+ //- Compute the QR decomposition with the column pivoting
+ // (QSB:Fig. 1)
+ void qrPivot
+ (
+ MatrixType& A
+ );
+
+ //- Compute output matrices in selected forms
+ //- using Householder reflectors
+ // (M:Section 4)
+ inline void applyHouseholder
+ (
+ MatrixType& A,
+ const RMatrix& reflector,
+ const label k = 0
+ );
+
//- Solve the linear system with the Field argument x initialized to
//- the appropriate transformed source (e.g. Q.T()*source)
//- and return the solution in x
@@ -90,27 +246,99 @@ private:
) const;
+protected:
+
+ // Protected Member Functions
+
+ //- Compute Householder reflector on a given matrix column, u
+ // (M:Eq.6.814)
+ inline RMatrix householderReflector
+ (
+ RMatrix u
+ );
+
+ //- Apply (in-place) Householder reflectors from the left side: u*A
+ inline void applyLeftReflector
+ (
+ MatrixType& A,
+ const RMatrix& reflector,
+ const label k = 0,
+ const label k1 = 0
+ );
+
+ //- Apply (in-place) Householder reflectors from the right side: (u*A)*u
+ inline void applyRightReflector
+ (
+ MatrixType& A,
+ const RMatrix& reflector,
+ const label k = 0
+ );
+
+
public:
// Constructors
//- Construct null
- inline QRMatrix();
+ QRMatrix();
- //- Construct decomposing given matrix
- QRMatrix(const MatrixType& M);
+ //- Construct QRMatrix without performing the decomposition
+ QRMatrix
+ (
+ const outputTypes outputType,
+ const storeMethods storeMethod = storeMethods::OUT_OF_PLACE,
+ const colPivoting colPivot = colPivoting::FALSE
+ );
+
+ //- Construct QRMatrix and perform the QR decomposition
+ QRMatrix
+ (
+ MatrixType& A,
+ const outputTypes outputType,
+ const storeMethods storeMethod = storeMethods::OUT_OF_PLACE,
+ const colPivoting colPivot = colPivoting::FALSE
+ );
+
+ //- Construct QRMatrix and perform the QR decomposition
+ QRMatrix
+ (
+ const MatrixType& A,
+ const outputTypes outputType,
+ const colPivoting colPivot = colPivoting::FALSE
+ );
// Member Functions
- //- Return Q-matrix
- inline const QMatrixType& Q() const;
+ // Information
+
+ //- Return the unitary similarity matrix
+ // Includes implicit round-to-zero as mutable operation
+ inline const MatrixType& Q() const;
+
+ //- Return the upper triangular matrix
+ inline const MatrixType& R() const;
+
+ //- Return the permutation order (P) as a list
+ inline const labelList& orderP() const;
+
+ //- Create and return the permutation matrix
+ inline SMatrix P() const;
- //- Return R-matrix
- inline const RMatrixType& R() const;
- //- Decompose given matrix
- void decompose(const MatrixType& M);
+ // Algorithm
+
+ //- Compute QR decomposition according to constructor settings
+ void decompose
+ (
+ MatrixType& A
+ );
+
+ //- Compute QR decomposition according to constructor settings
+ void decompose
+ (
+ const MatrixType& A
+ );
//- Solve the linear system with the given source
//- and return the solution in the argument x
@@ -145,10 +373,38 @@ public:
) const;
//- Return the inverse of (Q*R), so that solving x = (Q*R).inv()*source
- QMatrixType inv() const;
+ SMatrix inv() const;
+
+ //- Solve a row-echelon-form linear system starting from the bottom
+ // Back substitution: Ax = b where:
+ // A = Non-singular upper-triangular square matrix (m-by-m)
+ // x = Solution (m-by-p)
+ // b = Source (m-by-p)
+ RMatrix backSubstitution
+ (
+ const SMatrix& A,
+ const RMatrix& b
+ );
};
+// * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * * //
+
+//- (Out-of-place) Moore-Penrose inverse of singular/non-singular
+//- square/rectangular scalar/complex matrices
+// (KPP:p. 9834; KP:p. 648)
+// The tolerance (i.e. tol) to ensure the R1 matrix full-rank is given as 1e-13
+// in the original paper (KPP:p. 9834). Nonetheless, the tolerance = 1e-5
+// given by (TA; mentioned in (KPP:p. 9832)) was observed to be more robust
+// for the tested scenarios.
+template<class MatrixType>
+MatrixType pinv
+(
+ const MatrixType& A,
+ const scalar tol = 1e-5
+);
+
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H b/src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H
index 69102fa722e9459ad63c09d6f50573da292f777a..929f8b6fc89f986388923bf6f51d52718343de9f 100644
--- a/src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H
+++ b/src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H
@@ -26,29 +26,167 @@ License
\*---------------------------------------------------------------------------*/
-// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class MatrixType>
-inline Foam::QRMatrix<MatrixType>::QRMatrix()
-{}
+inline void Foam::QRMatrix<MatrixType>::applyHouseholder
+(
+ MatrixType& A,
+ const RMatrix& reflector,
+ const label k
+)
+{
+ applyLeftReflector(A, reflector, k, k);
+
+ if (outputType_ == outputTypes::FULL_QR)
+ {
+ applyRightReflector(Q_, reflector, k);
+ }
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
+template<class MatrixType>
+inline Foam::RectangularMatrix<typename MatrixType::cmptType>
+Foam::QRMatrix<MatrixType>::householderReflector
+(
+ RMatrix u
+)
+{
+ #ifdef FULLDEBUG
+ // Check if the given RectangularMatrix is effectively a column vector
+ if (u.n() != 1)
+ {
+ FatalErrorInFunction
+ << "Input matrix is not a column vector." << exit(FatalError);
+ }
+ #endif
+
+ scalar magnitude(mag(u(0,0)));
+ if (magnitude < VSMALL)
+ {
+ magnitude = SMALL;
+ #if FULLDEBUG
+ FatalErrorInFunction
+ << "Almost zero leading elem in Householder matrix."
+ << abort(FatalError);
+ #endif
+ }
+
+ u(0,0) += u(0,0)/magnitude*u.columnNorm(0);
+
+ scalar colNorm(u.columnNorm(0));
+ if (colNorm < VSMALL)
+ {
+ colNorm = SMALL;
+ #if FULLDEBUG
+ FatalErrorInFunction
+ << "Almost zero norm in the Householder matrix."
+ << abort(FatalError);
+ #endif
+ }
+
+ u /= cmptType(colNorm);
+
+ return u;
+}
+
+
+template<class MatrixType>
+inline void Foam::QRMatrix<MatrixType>::applyLeftReflector
+(
+ MatrixType& A,
+ const RMatrix& reflector,
+ const label k,
+ const label k1
+)
+{
+ // const RMatrix& A0(A.subMatrix(k1, k));
+ // A0 -= (cmptType(2)*reflector)*(reflector & A0);
+
+ for (label j = k; j < A.n(); ++j)
+ {
+ cmptType sum = Zero;
+ for (label i = 0; i < reflector.m(); ++i)
+ {
+ sum += Detail::conj(reflector(i, 0))*A(i + k1, j);
+ }
+
+ sum *= cmptType(2);
+ for (label i = 0; i < reflector.m(); ++i)
+ {
+ A(i + k1, j) -= reflector(i, 0)*sum;
+ }
+ }
+}
+
+
+template<class MatrixType>
+inline void Foam::QRMatrix<MatrixType>::applyRightReflector
+(
+ MatrixType& A,
+ const RMatrix& reflector,
+ const label k
+)
+{
+ // const RMatrix A0(A.subMatrix(0, k));
+ // A0 -= ((A0*reflector)^(cmptType(2)*reflector))
+
+ for (label i = 0; i < A.m(); ++i)
+ {
+ cmptType sum = Zero;
+ for (label j = 0; j < reflector.m(); ++j)
+ {
+ sum += A(i, j + k)*reflector(j, 0);
+ }
+
+ sum *= cmptType(2);
+ for (label j = 0; j < reflector.m(); ++j)
+ {
+ A(i, j + k) -= Detail::conj(reflector(j, 0))*sum;
+ }
+ }
+}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class MatrixType>
-inline const typename Foam::QRMatrix<MatrixType>::QMatrixType&
-Foam::QRMatrix<MatrixType>::Q() const
+inline const MatrixType& Foam::QRMatrix<MatrixType>::Q() const
{
+ const_cast<MatrixType&>(Q_).round();
+
return Q_;
}
template<class MatrixType>
-inline const typename Foam::QRMatrix<MatrixType>::RMatrixType&
-Foam::QRMatrix<MatrixType>::R() const
+inline const MatrixType& Foam::QRMatrix<MatrixType>::R() const
{
return R_;
}
+template<class MatrixType>
+inline const Foam::labelList& Foam::QRMatrix<MatrixType>::orderP() const
+{
+ return P_;
+}
+
+
+template<class MatrixType>
+inline Foam::SquareMatrix<typename Foam::QRMatrix<MatrixType>::cmptType>
+Foam::QRMatrix<MatrixType>::P() const
+{
+ SquareMatrix<cmptType> permMat(P_.size(), Zero);
+
+ forAll(P_, jcol)
+ {
+ permMat(P_[jcol], jcol) = pTraits<cmptType>::one;
+ }
+
+ return permMat;
+}
+
+
// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H
index bb1877e078b585dab9565e6d0eab2a42bd1724be..5b09ea7568b0e153706f494c4b40122a2a4b2843 100644
--- a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H
+++ b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H
@@ -80,6 +80,11 @@ public:
//- initializing all elements to the given value
inline RectangularMatrix(const label m, const label n, const Type& val);
+ //- Construct for given number of rows/columns
+ //- initializing all elements to zero, and diagonal to one
+ template<class AnyType>
+ inline RectangularMatrix(const labelPair& dims, const Identity<AnyType>);
+
//- Construct given number of rows/columns
inline explicit RectangularMatrix(const labelPair& dims);
diff --git a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H
index b014e290699e85f54a09937c4ca13d57d56f0db7..7cef3593e77c223e6e24c17022a14b6797ab9374 100644
--- a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H
+++ b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H
@@ -80,6 +80,23 @@ inline Foam::RectangularMatrix<Type>::RectangularMatrix
{}
+template<class Type>
+template<class AnyType>
+inline Foam::RectangularMatrix<Type>::RectangularMatrix
+(
+ const labelPair& dims,
+ const Identity<AnyType>
+)
+:
+ Matrix<RectangularMatrix<Type>, Type>(dims.first(), dims.second(), Zero)
+{
+ for (label i = 0; i < min(dims.first(), dims.second()); ++i)
+ {
+ this->operator()(i, i) = pTraits<Type>::one;
+ }
+}
+
+
template<class Type>
inline Foam::RectangularMatrix<Type>::RectangularMatrix
(
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
index b3ac5eb9f3bafa3378396e04b07ac70abad96563..e404910e069a8e4fad048a10d2056e064fcd4ed8 100644
--- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
+++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
@@ -85,6 +85,9 @@ public:
template<class AnyType>
inline SquareMatrix(const label n, const Identity<AnyType>);
+ template<class AnyType>
+ inline SquareMatrix(const labelPair& dims, const Identity<AnyType>);
+
//- Construct given number of rows/columns (checked to be equal)
// For constructor consistency in rectangular matrices
inline explicit SquareMatrix(const labelPair& dims);
@@ -129,12 +132,23 @@ public:
//- Resize the matrix preserving the elements
inline void resize(const label m);
+ //- Resize the matrix preserving the elements (compatibility)
+ inline void resize(const label m, const label n);
+
//- Resize the matrix preserving the elements
inline void setSize(const label m);
//- Resize the matrix without reallocating storage (unsafe)
inline void shallowResize(const label m);
+ // Check
+
+ //- Return true if the square matrix is effectively symmetric/Hermitian
+ inline bool symmetric() const;
+
+ //- Return true if the square matrix is reduced tridiagonal
+ inline bool tridiagonal() const;
+
// Member Operators
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
index b41b5e97bea1814744e284cb82e6676f6a6aacd5..75e893e068f8c9d91e32d83aaec652364b967455 100644
--- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
+++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
@@ -90,6 +90,25 @@ inline Foam::SquareMatrix<Type>::SquareMatrix
}
+template<class Type>
+template<class AnyType>
+inline Foam::SquareMatrix<Type>::SquareMatrix
+(
+ const labelPair& dims,
+ const Identity<AnyType>
+)
+:
+ Matrix<SquareMatrix<Type>, Type>(dims, Zero)
+{
+ CHECK_MATRIX_IS_SQUARE(dims.first(), dims.second());
+
+ for (label i = 0; i < dims.first(); ++i)
+ {
+ this->operator()(i, i) = pTraits<Type>::one;
+ }
+}
+
+
template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
@@ -206,6 +225,19 @@ inline void Foam::SquareMatrix<Type>::resize(const label m)
}
+template<class Type>
+inline void Foam::SquareMatrix<Type>::resize(const label m, const label n)
+{
+ if (m != n)
+ {
+ FatalErrorInFunction<< "Disallowed use of resize() for SquareMatrix"
+ << abort(FatalError);
+ }
+
+ Matrix<SquareMatrix<Type>, Type>::resize(m, m);
+}
+
+
template<class Type>
inline void Foam::SquareMatrix<Type>::setSize(const label m)
{
@@ -220,6 +252,49 @@ inline void Foam::SquareMatrix<Type>::shallowResize(const label m)
}
+template<class Type>
+inline bool Foam::SquareMatrix<Type>::symmetric() const
+{
+ for (label n = 0; n < this->n() - 1; ++n)
+ {
+ for (label m = this->m() - 1; n < m; --m)
+ {
+ if (SMALL < mag((*this)(n, m) - (*this)(m, n)))
+ {
+ return false;
+ }
+ }
+ }
+ return true;
+}
+
+
+template<class Type>
+inline bool Foam::SquareMatrix<Type>::tridiagonal() const
+{
+ for (label i = 0; i < this->m(); ++i)
+ {
+ for (label j = 0; j < this->n(); ++j)
+ {
+ const Type& val = (*this)(i, j);
+
+ if ((i == j) || (i - 1 == j) || (i + 1 == j))
+ {
+ if (mag(val) < SMALL)
+ {
+ return false;
+ }
+ }
+ else if (SMALL < mag(val))
+ {
+ return false;
+ }
+ }
+ }
+ return true;
+}
+
+
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class Type>
diff --git a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/kkLOmega/kkLOmega.C b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/kkLOmega/kkLOmega.C
index 34462139a55c0637b3b01bffc1bab06286e63291..31498a6be4415c91e7ca642997be32532db413c0 100644
--- a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/kkLOmega/kkLOmega.C
+++ b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/kkLOmega/kkLOmega.C
@@ -636,7 +636,7 @@ void kkLOmega::correct()
*fINT()
*Cmu(sqrt(S2))*sqrt(ktS)*lambdaEff
);
- const volScalarField Pkt(nuts*S2);
+ volScalarField Pkt(this->GName(), nuts*S2);
const volScalarField ktL(kt_ - ktS);
const volScalarField ReOmega(sqr(y_)*Omega/nu());
diff --git a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/kkLOmega/kkLOmega.H b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/kkLOmega/kkLOmega.H
index 934911e6ea5315cc56ce7e1a4c686222246fb8bd..b7fe7c295c6b056d6ab304382b75b94c8f2a5163 100644
--- a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/kkLOmega/kkLOmega.H
+++ b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/kkLOmega/kkLOmega.H
@@ -34,25 +34,22 @@ Description
Low Reynolds-number k-kl-omega turbulence model for
incompressible flows.
- This turbulence model is described in:
+ Reference:
\verbatim
- Walters, D. K., & Cokljat, D. (2008).
- A three-equation eddy-viscosity model for Reynolds-averaged
- Navier–Stokes simulations of transitional flow.
- Journal of Fluids Engineering, 130(12), 121401.
+ Standard model:
+ Walters, D. K., & Cokljat, D. (2008).
+ A three-equation eddy-viscosity model for Reynolds-averaged
+ Navier–Stokes simulations of transitional flow.
+ Journal of Fluids Engineering, 130(12), 121401.
+ DOI:10.1115/1.2979230
+
+ Typo corrections:
+ Furst, J. (2013).
+ Numerical simulation of transitional flows
+ with laminar kinetic energy.
+ Engineering MECHANICS, 20(5), 379-388.
\endverbatim
- however the paper contains several errors which must be corrected for the
- model to operation correctly as explained in
-
- \verbatim
- Furst, J. (2013).
- Numerical simulation of transitional flows with laminar kinetic energy.
- Engineering MECHANICS, 20(5), 379-388.
- \endverbatim
-
- All these corrections and updates are included in this implementation.
-
The default model coefficients are
\verbatim
kkLOmegaCoeffs
diff --git a/src/TurbulenceModels/turbulenceModels/laminar/Maxwell/Maxwell.H b/src/TurbulenceModels/turbulenceModels/laminar/Maxwell/Maxwell.H
index 98fefe5ad0fb27c1fb7caa570f86506878110d02..1f3ce3c6863ed322761e1bcc4b7a2d2a5306698a 100644
--- a/src/TurbulenceModels/turbulenceModels/laminar/Maxwell/Maxwell.H
+++ b/src/TurbulenceModels/turbulenceModels/laminar/Maxwell/Maxwell.H
@@ -47,6 +47,7 @@ Description
Comparison of different formulations for the numerical calculation
of unsteady incompressible viscoelastic fluid flow.
Adv. Appl. Math. Mech, 4, 483-502.
+ DOI:10.4208/aamm.10-m1010
\endverbatim
SourceFiles
@@ -78,7 +79,7 @@ class Maxwell
protected:
- // Protected data
+ // Protected Data
// Model coefficients
@@ -153,7 +154,7 @@ public:
) const;
//- Solve the turbulence equations and correct eddy-Viscosity and
- // related properties
+ //- related properties
virtual void correct();
};
diff --git a/src/TurbulenceModels/turbulenceModels/laminar/Stokes/Stokes.H b/src/TurbulenceModels/turbulenceModels/laminar/Stokes/Stokes.H
index 3ef35711287e25f6123cc676a0fe5bbcd6e9b91c..a2da356cd3a7a9f1414f1aa3fa9e20a3810c1e89 100644
--- a/src/TurbulenceModels/turbulenceModels/laminar/Stokes/Stokes.H
+++ b/src/TurbulenceModels/turbulenceModels/laminar/Stokes/Stokes.H
@@ -25,7 +25,10 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
- Foam::Stokes
+ Foam::laminarModels::Stokes
+
+Group
+ grpLaminar
Description
Turbulence model for Stokes flow.
@@ -127,7 +130,7 @@ public:
virtual tmp<volScalarField> k() const;
//- Return the turbulence kinetic energy dissipation rate,
- // i.e. 0 for Stokes flow
+ //- i.e. 0 for Stokes flow
virtual tmp<volScalarField> epsilon() const;
//- Return the Reynolds stress tensor, i.e. 0 for Stokes flow
diff --git a/tutorials/combustion/XiDyMFoam/annularCombustorTurbine/Allrun b/tutorials/combustion/XiDyMFoam/annularCombustorTurbine/Allrun
index 20ddac1d57e38d7142c839a8399a7398ba48a3c9..3e78f257360ab88b5b01d28577b875ad4180c81a 100755
--- a/tutorials/combustion/XiDyMFoam/annularCombustorTurbine/Allrun
+++ b/tutorials/combustion/XiDyMFoam/annularCombustorTurbine/Allrun
@@ -11,8 +11,6 @@ runApplication decomposePar -force
runParallel potentialFoam -pName pPotential -initialiseUBCs
-rm -f processor*/0/phi
-
runParallel XiDyMFoam
#------------------------------------------------------------------------------
diff --git a/tutorials/incompressible/adjointOptimisationFoam/motorBike/Allrun b/tutorials/incompressible/adjointOptimisationFoam/motorBike/Allrun
index 3b20001a5058dccc502892ff555d5ae05bbc714a..46c60b78fc43150f4cdd38f64eb9a489adff7fae 100755
--- a/tutorials/incompressible/adjointOptimisationFoam/motorBike/Allrun
+++ b/tutorials/incompressible/adjointOptimisationFoam/motorBike/Allrun
@@ -33,7 +33,7 @@ runParallel $decompDict snappyHexMesh -overwrite
restore0Dir -processor
runParallel $decompDict patchSummary
-runParallel $decompDict potentialFoam
+runParallel $decompDict potentialFoam -writephi
runParallel $decompDict checkMesh -writeFields '(nonOrthoAngle)' -constant
runParallel $decompDict $(getApplication)
diff --git a/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/Allclean b/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/Allclean
index d670ed7d367cbac27a704cb51c2587d991e01cea..cc7e6894a0f7174c1708520e12a307fa97022eb7 100755
--- a/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/Allclean
+++ b/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/Allclean
@@ -4,8 +4,7 @@ cd "${0%/*}" || exit # Run from this directory
#------------------------------------------------------------------------------
cleanCase
-rm -rf *.dat validation/*.eps
-
+rm -rf *.{dat,png,txt} system/controlDict constant/turbulenceProperties results
wclean validation/WatersKing
#------------------------------------------------------------------------------
diff --git a/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/Allrun b/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/Allrun
index 893cbf8c2a16f8ca0c2205557cc42fd7ffab2a96..9b2ab05be8ec105b1d3a174faa103ab68e4bb86e 100755
--- a/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/Allrun
+++ b/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/Allrun
@@ -1,14 +1,97 @@
#!/bin/sh
cd "${0%/*}" || exit # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions # Tutorial run functions
+. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions # Tutorial clean functions
#------------------------------------------------------------------------------
-runApplication blockMesh
-runApplication $(getApplication)
+models="
+Maxwell
+Stokes
+"
-wmake validation/WatersKing
-runApplication WatersKing
+endTime=2
+if notTest $@
+then
+ endTime=25
+fi
-( cd validation && ./createGraph )
+#------------------------------------------------------------------------------
+
+# Compute the case in 'serial' mode, and collect the data
+#
+# $1 = endTime
+# $* = models
+# ----
+serialRun() {
+ endTime=$1
+ shift 1
+ models=$*
+
+ sed "s|END_TIME|$endTime|g" system/controlDict.template \
+ > system/controlDict
+ resultDir="results"
+
+ runApplication blockMesh
+ wmake validation/WatersKing
+
+ for model in $models
+ do
+ sed "s|LAMINAR_MODEL|$model|g" constant/turbulenceProperties.template \
+ > constant/turbulenceProperties
+
+ # Compute numerical data
+ runApplication $(getApplication)
+ tail -n +4 postProcessing/probes/0/U | \
+ tr -s " " | tr -d '(' | cut -d " " -f2-3 > "${model}.txt"
+
+ # Collect numerical data
+ modelDir="$resultDir/$model"
+ [ -d "$modelDir" ] || mkdir -p "$modelDir"
+ mv -f postProcessing log.* "$modelDir"
+ cleanTimeDirectories
+ done
+
+ runApplication WatersKing
+}
+
+
+# Plot streamwise flow speed at y=1.0 [m] as a function of time
+#
+# $* = models
+# ----
+plot() {
+ # Require gnuplot
+ command -v gnuplot >/dev/null || {
+ echo "gnuplot not found - skipping graph creation" 1>&2
+ exit 1
+ }
+
+ models=$*
+ endTime=$(foamDictionary -entry endTime -value system/controlDict)
+
+ gnuplot<<PLT
+ set terminal pngcairo font "helvetica,16" size 800,600
+ set output "planarPoiseuille.png"
+ set grid
+ set key right top
+ set xrange [0:"$endTime"]
+ set yrange [0:8]
+ set xlabel "t [s]"
+ set ylabel "U_x [m/s]" rotate by 0 offset 3,0,0
+
+ results=system("ls *.txt")
+ names="${models[*]}"
+ plot \
+ "WatersKing.dat" w lines t "Analytical" lt -1, \
+ for [i=1:words(results)] word(results, i) t word(names, i) \
+ w linespoints pointinterval 100 lt i pt 6 ps 1.5
+
+PLT
+}
+
+#------------------------------------------------------------------------------
+
+serialRun $endTime $models
+plot $models
#------------------------------------------------------------------------------
diff --git a/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/constant/turbulenceProperties b/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/constant/turbulenceProperties.template
similarity index 96%
rename from tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/constant/turbulenceProperties
rename to tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/constant/turbulenceProperties.template
index c782dcde8b1f160fab70d0526e4c9011a5f661d8..aefd35ba5d94d9dce492e05f9608191c09c22c07 100644
--- a/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/constant/turbulenceProperties
+++ b/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/constant/turbulenceProperties.template
@@ -18,7 +18,7 @@ simulationType laminar;
laminar
{
- laminarModel Maxwell;
+ laminarModel LAMINAR_MODEL;
MaxwellCoeffs
{
diff --git a/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/system/controlDict b/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/system/controlDict.template
similarity index 97%
rename from tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/system/controlDict
rename to tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/system/controlDict.template
index f26ec0685b2cb6d7660a061c37480546e5c0b880..6e8a4bada10ba9ab0bbf1d9f872d7f4778ba3a8f 100644
--- a/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/system/controlDict
+++ b/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/system/controlDict.template
@@ -22,7 +22,7 @@ startTime 0;
stopAt endTime;
-endTime 25;
+endTime END_TIME;
deltaT 5e-3;
diff --git a/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/system/fvSchemes b/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/system/fvSchemes
index 8b06c3af277151032b975cbf3962f7f4fd70c318..15a53aa38a4b12b2b2908e8a611f6db49e00c537 100644
--- a/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/system/fvSchemes
+++ b/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/system/fvSchemes
@@ -27,13 +27,12 @@ gradSchemes
divSchemes
{
default none;
-
div(phi,U) Gauss linearUpwind grad(U);
div(phi,sigma) Gauss vanAlbada;
-
div(sigma) Gauss linear;
div((nu*dev2(T(grad(U))))) Gauss linear;
div((nuM*grad(U))) Gauss linear;
+ div((nuEff*dev2(T(grad(U))))) Gauss linear;
}
laplacianSchemes
diff --git a/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/validation/WatersKing/WatersKing.C b/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/validation/WatersKing/WatersKing.C
index b215bc8d1b5742f7f134229f814f80430886c8e3..8e78c17f15892b0dc43b21cb7a9e5f4c2d827f85 100644
--- a/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/validation/WatersKing/WatersKing.C
+++ b/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/validation/WatersKing/WatersKing.C
@@ -37,9 +37,10 @@ Description
Rheologica Acta, 9, 345-355.
Amoreira, L. J., & Oliveira, P. J. (2010).
- Comparison of different formulations for the numerical
- calculation of unsteady incompressible viscoelastic fluid
- flow. Adv. Appl. Math. Mech, 4, 483-502.
+ Comparison of different formulations for the numerical calculation of
+ unsteady incompressible viscoelastic fluid flow.
+ Adv. Appl. Math. Mech, 4, 483-502.
+ DOI:10.4208/aamm.10-m1010
\endverbatim
\*---------------------------------------------------------------------------*/
diff --git a/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/validation/createGraph b/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/validation/createGraph
deleted file mode 100755
index 77a23fedf430171167919a073cbe4700fc7e4fed..0000000000000000000000000000000000000000
--- a/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/validation/createGraph
+++ /dev/null
@@ -1,27 +0,0 @@
-#!/bin/sh
-
-tail -n +4 ../postProcessing/probes/0/U | \
- tr -s " " | tr -d '(' | cut -d " " -f2-3 > ../Numerical.dat
-
-# Require gnuplot
-command -v gnuplot >/dev/null || {
- echo "gnuplot not found - skipping graph creation" 1>&2
- exit 1
-}
-
-gnuplot<<EOF
- set terminal postscript eps color enhanced "Helvetica,20"
- set output "planarPoiseuille.eps"
- set xlabel "Time / [s]" font "Helvetica,24"
- set ylabel "Velocity / [m/s]" font "Helvetica,24"
- set grid
- set key right top
- set xrange [0:25]
- set yrange [0:8]
- plot \
- "../Numerical.dat" t "OpenFOAM (every 100 pts)" \
- with linespoints pointinterval 100 lt 1 pt 6 ps 1.5, \
- "../WatersKing.dat" with lines t "Analytical" lt -1
-EOF
-
-#------------------------------------------------------------------------------
diff --git a/tutorials/incompressible/pisoFoam/LES/motorBike/motorBike/Allrun b/tutorials/incompressible/pisoFoam/LES/motorBike/motorBike/Allrun
index 78c665597a5ab2611c7c04aa995e05ec36b33b96..22cc513515f7491f02fc135d62d47fb6882decd5 100755
--- a/tutorials/incompressible/pisoFoam/LES/motorBike/motorBike/Allrun
+++ b/tutorials/incompressible/pisoFoam/LES/motorBike/motorBike/Allrun
@@ -22,7 +22,7 @@ find . -iname '*level*' -type f -delete
restore0Dir -processor
runParallel renumberMesh -overwrite
-runParallel potentialFoam -initialiseUBCs
+runParallel potentialFoam -initialiseUBCs -writephi
runParallel checkMesh -writeFields '(nonOrthoAngle)' -constant
runParallel $(getApplication)
diff --git a/tutorials/incompressible/simpleFoam/motorBike/Allrun b/tutorials/incompressible/simpleFoam/motorBike/Allrun
index 5646a4fbb34f762071744f10769e67ea952c5dd9..1dcc3f559be43e31dcee5f778a9c8496a1e52e7e 100755
--- a/tutorials/incompressible/simpleFoam/motorBike/Allrun
+++ b/tutorials/incompressible/simpleFoam/motorBike/Allrun
@@ -33,7 +33,7 @@ runParallel $decompDict snappyHexMesh -overwrite
restore0Dir -processor
runParallel $decompDict patchSummary
-runParallel $decompDict potentialFoam
+runParallel $decompDict potentialFoam -writephi
runParallel $decompDict checkMesh -writeFields '(nonOrthoAngle)' -constant
runParallel $decompDict $(getApplication)
diff --git a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/Allrun b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/Allrun
index 815dc6f3d0a6ce5f5ff4032f110b294f18ccfcfa..9ec0596bab62acca9301776a93baa8fb237ea0cb 100755
--- a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/Allrun
+++ b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/Allrun
@@ -11,8 +11,6 @@ restore0Dir
# initialise with potentialFoam solution
runApplication potentialFoam
-rm -f 0/phi
-
# run the solver
runApplication $(getApplication)
diff --git a/tutorials/lagrangian/reactingParcelFoam/verticalChannelLTS/Allrun b/tutorials/lagrangian/reactingParcelFoam/verticalChannelLTS/Allrun
index 815dc6f3d0a6ce5f5ff4032f110b294f18ccfcfa..9ec0596bab62acca9301776a93baa8fb237ea0cb 100755
--- a/tutorials/lagrangian/reactingParcelFoam/verticalChannelLTS/Allrun
+++ b/tutorials/lagrangian/reactingParcelFoam/verticalChannelLTS/Allrun
@@ -11,8 +11,6 @@ restore0Dir
# initialise with potentialFoam solution
runApplication potentialFoam
-rm -f 0/phi
-
# run the solver
runApplication $(getApplication)
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/Allrun b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/Allrun
index eaf8f1bdb01339e76aad4962d03efb534c97e523..a5d9a2cb057a3a418c92f0fb6eac32e2eef1c21f 100755
--- a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/Allrun
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/Allrun
@@ -9,9 +9,6 @@ runApplication blockMesh
runApplication potentialFoam
-# Remove incompatible (volumetric) flux field
-rm -f 0/phi
-
runApplication $(getApplication)
#------------------------------------------------------------------------------
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/Allrun-parallel b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/Allrun-parallel
index 2392ee3aec70c2a5a8951f2ac33fe27582e3b718..8bb64db8194a2d9f23ba608c072195988d995a34 100755
--- a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/Allrun-parallel
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/Allrun-parallel
@@ -9,9 +9,6 @@ runApplication blockMesh
runApplication potentialFoam
-# Remove incompatible (volumetric) flux field
-rm -f 0/phi
-
runApplication decomposePar
runParallel $(getApplication)
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/Allrun b/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/Allrun
index b2e91a27a0ea2ed5f85acf1bf125f42448055156..41ff036dccd17c69303d84fa7bf09a835c05b957 100755
--- a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/Allrun
+++ b/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/Allrun
@@ -15,7 +15,7 @@ runApplication snappyHexMesh -overwrite
restore0Dir
# Initialise with potentialFoam solution
-runApplication potentialFoam -pName p_rgh
+runApplication potentialFoam -pName p_rgh -writephi
# Run the solver
runApplication $(getApplication)
diff --git a/tutorials/preProcessing/createZeroDirectory/motorBike/Allrun b/tutorials/preProcessing/createZeroDirectory/motorBike/Allrun
index 77d273e307d5a0e4bdc32ce7d840f7e54a6f0ca1..b99d37cc556fcb3dc651b5781677e9d23f5bfba6 100755
--- a/tutorials/preProcessing/createZeroDirectory/motorBike/Allrun
+++ b/tutorials/preProcessing/createZeroDirectory/motorBike/Allrun
@@ -15,7 +15,7 @@ runParallel snappyHexMesh -overwrite
runParallel createZeroDirectory
runParallel patchSummary
-runParallel potentialFoam -noFunctionObjects
+runParallel potentialFoam -noFunctionObjects -writephi
runParallel $(getApplication)
runApplication reconstructParMesh -constant