diff --git a/applications/test/TestTools/TestTools.H b/applications/test/TestTools/TestTools.H
index 3fed3b7d5dfe40db2f912528de5a20d6ed8885fc..82dfe5e38b11e903675821854043468352700acf 100644
--- a/applications/test/TestTools/TestTools.H
+++ b/applications/test/TestTools/TestTools.H
@@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
- Copyright (C) 2020-2021 OpenCFD Ltd.
+ Copyright (C) 2020-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -35,6 +35,7 @@ using namespace Foam;
#include "complex.H"
#include "Matrix.H"
#include "Random.H"
+#include <chrono>
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -46,6 +47,27 @@ unsigned nTest_ = 0;
unsigned nFail_ = 0;
+// Return the absolute tolerance value for bitwise comparisons of floatScalars
+floatScalar getTol(floatScalar)
+{
+ return 1e-2;
+}
+
+
+// Return the absolute tolerance value for bitwise comparisons of doubleScalars
+doubleScalar getTol(doubleScalar)
+{
+ return 1e-10;
+}
+
+
+// Return the absolute tolerance value for bitwise comparisons of doubleScalars
+doubleScalar getTol(complex)
+{
+ return 1e-10;
+}
+
+
// Create a non-complex random Matrix.
template<class MatrixType>
typename std::enable_if
@@ -174,10 +196,14 @@ typename std::enable_if
const Type& x,
const Type& y,
const scalar absTol = 0, //<! useful for cmps near zero
- const scalar relTol = 1e-8 //<! are values the same within 8 decimals
+ const scalar relTol = 1e-8, //<! are values the same within 8 decimals
+ const bool verbose = false
)
{
- Info<< msg << x << "?=" << y << endl;
+ if (verbose)
+ {
+ Info<< msg << x << "?=" << y << endl;
+ }
unsigned nFail = 0;
@@ -212,10 +238,14 @@ typename std::enable_if
const Type& x,
const Type& y,
const scalar absTol = 0,
- const scalar relTol = 1e-8
+ const scalar relTol = 1e-8,
+ const bool verbose = false
)
{
- Info<< msg << x << "?=" << y << endl;
+ if (verbose)
+ {
+ Info<< msg << x << "?=" << y << endl;
+ }
unsigned nFail = 0;
@@ -253,10 +283,14 @@ typename std::enable_if
const Type1& x,
const Type2& y,
const scalar absTol = 0,
- const scalar relTol = 1e-8
+ const scalar relTol = 1e-8,
+ const bool verbose = false
)
{
- Info<< msg << x << "?=" << y << endl;
+ if (verbose)
+ {
+ Info<< msg << x << "?=" << y << endl;
+ }
unsigned nFail = 0;
@@ -284,10 +318,14 @@ void cmp
(
const word& msg,
const bool x,
- const bool y
+ const bool y,
+ const bool verbose = false
)
{
- Info<< msg << x << "?=" << y << endl;
+ if (verbose)
+ {
+ Info<< msg << x << "?=" << y << endl;
+ }
unsigned nFail = 0;
@@ -306,4 +344,39 @@ void cmp
}
+// 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;
+ }
+ }
+}
+
+
// ************************************************************************* //
diff --git a/applications/test/matrices/QRMatrix/Make/options b/applications/test/matrices/QRMatrix/Make/options
index 4e772fdf9d7bc94221d127458f9d2ca32850fe69..b23590707c8d302fe75ce411d5bf921c029a7ff3 100644
--- a/applications/test/matrices/QRMatrix/Make/options
+++ b/applications/test/matrices/QRMatrix/Make/options
@@ -1,2 +1,2 @@
-/* EXE_INC = -I$(LIB_SRC)/finiteVolume/lnInclude */
+EXE_INC = -I../../TestTools
/* EXE_LIBS = -lfiniteVolume */
diff --git a/applications/test/matrices/QRMatrix/Test-QRMatrix.C b/applications/test/matrices/QRMatrix/Test-QRMatrix.C
index 679072715257733728458b9f2810891f0fa04053..9e51d67eb023484d000252082c79bade54d44cf1 100644
--- a/applications/test/matrices/QRMatrix/Test-QRMatrix.C
+++ b/applications/test/matrices/QRMatrix/Test-QRMatrix.C
@@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
- Copyright (C) 2019-2020 OpenCFD Ltd.
+ Copyright (C) 2020-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -23,946 +23,338 @@ License
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+Application
+ Test-QRMatrix
+
+Description
+ Tests for \c QRMatrix constructors, and member functions
+ using \c doubleScalar and \c complex base types.
+
\*---------------------------------------------------------------------------*/
#include "MatrixTools.H"
#include "QRMatrix.H"
-#include "Random.H"
+#include "complex.H"
#include "IOmanip.H"
+#include "TestTools.H"
-using namespace Foam;
-using namespace Foam::MatrixTools;
+#define PRINTALL false
-#define equal MatrixTools::equal
-#define RUNALL true
-const bool verbose = true;
+using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-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
-)
+// Verify if Q is a unitary matrix for rectangular matrices
+template<class Type>
+void verify_Q(const RectangularMatrix<Type>& Q)
{
- Info<< "Empty " << objName << " = ";
- if (A.empty())
- {
- Info<< "true" << nl;
- }
- else
- {
- Info<< "false" << nl;
- }
+ // non-square unitary matrices are not possible - do nothing
}
-// Checks if Q matrix is a unitary matrix
-template<class MatrixType>
-void cross_check_QRMatrix
+// Verify if Q is a unitary matrix for square matrices
+template<class Type>
+void verify_Q
(
- const MatrixType& Aorig,
- const MatrixType& Q
+ const SquareMatrix<Type>& 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);
+ const SquareMatrix<Type> QTQ(Q&Q);
+ const SquareMatrix<Type> QQT(Q^Q);
+ SquareMatrix<Type> IMatrix(Q.sizes(), Zero);
for (label i = 0; i < IMatrix.n(); ++i)
{
- IMatrix(i,i) = pTraits<typename MatrixType::cmptType>::one;
+ IMatrix(i,i) = pTraits<Type>::one;
}
- equal(QTQ, IMatrix, verbose);
- equal(QQT, IMatrix, verbose);
- equal(QTQ, QQT, verbose);
+ cmp("# Q.T()*Q = I: ", flt(QTQ), flt(IMatrix), getTol(Type(0)));
+ cmp("# Q*Q.T() = I: ", flt(QQT), flt(IMatrix), getTol(Type(0)));
+ cmp("# QTQ = QQT: ", flt(QTQ), flt(QQT), getTol(Type(0)));
}
-// Checks if R matrix is an upper-triangular matrix
+// Verify if R is an upper-triangular matrix
template<class MatrixType>
-void cross_check_QRMatrix
+void verify_R
(
const MatrixType& R
)
{
- InfoNumPyFormat(R, "R");
-
+ DynamicList<scalar> ltriR;
+ label k = 0;
// 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);
+ ltriR.append(mag(R(i, j)));
+ ++k;
}
}
}
+ const List<scalar> ltri(k, Zero);
+ cmp("# R(i, j) = 0 for i > j: ", ltri, ltriR);
}
-// Checks if the given matrix can be reconstructed by Q and R matrices
+// Verify if R is an upper-triangular matrix with diag elems non-increasing
template<class MatrixType>
-void cross_check_QRMatrix
+void verify_R_pivoting
(
- 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);
+ // w.wiki/54m (Retrieved:20-06-19)
+ const List<scalar> diag0(mag(R.diag()));
+ List<scalar> diag1(diag0);
+ Foam::sort(diag1, std::greater<scalar>());
+ cmp
+ (
+ "# Diag elems of R non-increasing: |R11| >= |R22| >= ... >= |Rnn|: ",
+ diag0,
+ diag1
+ );
}
-// Checks if R matrix is an upper-triangular matrix, and obeys column pivoting
+// Verify if the given matrix can be reconstructed by Q and R matrices
template<class MatrixType>
-void cross_check_QRMatrix
+void verify_QR
(
- const MatrixType& Aorig,
- const labelList& orderP,
- const SquareMatrix<typename MatrixType::cmptType>& P,
+ const MatrixType& A,
+ const MatrixType& Q,
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);
- }
+ // mathworld.wolfram.com/QRDecomposition.html (Retrieved:18-06-19)
+ const MatrixType AReconstruct(Q*R);
+ cmp("# Q*R = A: ", flt(A), flt(AReconstruct));
}
-// Checks if the given matrix can be reconstructed by column-pivoted Q and R
+// Verify if the given matrix can be reconstructed by Q, R and P matrices
template<class MatrixType>
-void cross_check_QRMatrix
+void verify_QR_pivoting
(
- const MatrixType& Aorig,
+ const MatrixType& A,
const MatrixType& Q,
- const labelList& orderP,
- const SquareMatrix<typename MatrixType::cmptType>& P,
- const MatrixType& R
+ const MatrixType& R,
+ const SquareMatrix<typename MatrixType::cmptType>& P
)
{
- 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);
+ const MatrixType AP(A*P);
+ cmp("# Q*R = A*P: ", flt(AReconstruct), flt(AP));
}
-// Checks each constructor of QRMatrix type
+// Verify the given QR decomposition
template<class MatrixType>
-void verification_QRMatrix
+void verify
(
- MatrixType& A
+ const QRMatrix<MatrixType>& QRM,
+ const 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
+ const MatrixType& Q = QRM.Q();
+ const MatrixType& R = QRM.R();
- #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);
- }
+ #if PRINTALL
+ InfoNumPyFormat(A, "A");
+ InfoNumPyFormat(Q, "Q");
+ InfoNumPyFormat(R, "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
+ verify_Q(Q);
+ verify_R(R);
+ verify_QR(A, Q, R);
+ Info<< endl;
+}
- #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
+// Verify the given QR decomposition with column pivoting
+template<class MatrixType>
+void verify_pivoting
+(
+ const QRMatrix<MatrixType>& QRM,
+ const MatrixType& A
+)
+{
+ const MatrixType& Q = QRM.Q();
+ const MatrixType& R = QRM.R();
+ const SquareMatrix<typename MatrixType::cmptType>& P = QRM.P();
- #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);
- }
+ #if PRINTALL
+ InfoNumPyFormat(A, "A");
+ InfoNumPyFormat(Q, "Q");
+ InfoNumPyFormat(R, "R");
+ InfoNumPyFormat(P, "P");
#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
+ verify_Q(Q);
+ verify_R(R);
+ verify_R_pivoting(R);
+ verify_QR_pivoting(A, Q, R, P);
+ Info<< endl;
+}
- #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)
+// Verify various QR decompositions
+template<class MatrixType>
+void verify_decomposition
+(
+ const MatrixType& A
+)
+{
+ Info<< "## Verify various QR decompositions" << nl << endl;
{
- Info<< "# A | FULL_QR | IN_PLACE | colPivot = on" << nl;
- MatrixType A0(A);
+ Info<< "# modes::FULL" << nl;
QRMatrix<MatrixType> QRM
(
- A0,
- QRMatrix<MatrixType>::outputTypes::FULL_QR,
- QRMatrix<MatrixType>::storeMethods::IN_PLACE,
- QRMatrix<MatrixType>::colPivoting::TRUE
+ A,
+ QRMatrix<MatrixType>::modes::FULL
);
- 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);
+ verify(QRM, A);
}
- #endif
-
- // constructors with the const type specifier
- #if (0 | RUNALL)
{
- Info<< "# const A | FULL_R | colPivot = off" << nl;
- const MatrixType A0(A);
+ Info<< "# modes::ECONOMY" << nl;
QRMatrix<MatrixType> QRM
(
- A0,
- QRMatrix<MatrixType>::outputTypes::FULL_R,
- QRMatrix<MatrixType>::colPivoting::FALSE
+ A,
+ QRMatrix<MatrixType>::modes::ECONOMY
);
- const MatrixType& Q(QRM.Q());
- const MatrixType& R(QRM.R());
- isEmpty(Q, "Q");
- cross_check_QRMatrix(R);
+ verify(QRM, A);
}
- #endif
-
- #if (0 | RUNALL)
{
- Info<< "# const A | FULL_QR | colPivot = off" << nl;
- const MatrixType A0(A);
+ Info<< "# modes::FULL, outputs::BOTH_QR, column-pivoting" << nl;
QRMatrix<MatrixType> QRM
(
- A0,
- QRMatrix<MatrixType>::outputTypes::FULL_QR,
- QRMatrix<MatrixType>::colPivoting::FALSE
+ A,
+ QRMatrix<MatrixType>::modes::FULL,
+ QRMatrix<MatrixType>::outputs::BOTH_QR,
+ QRMatrix<MatrixType>::pivoting::TRUE
);
- const MatrixType& Q(QRM.Q());
- const MatrixType& R(QRM.R());
- cross_check_QRMatrix(Aorig, Q, R);
+ verify_pivoting(QRM, A);
}
- #endif
-
- #if (0 | RUNALL)
{
- Info<< "# const A | FULL_R | colPivot = on" << nl;
- const MatrixType A0(A);
+ Info<< "# modes::ECONOMY, outputs::BOTH_QR, column-pivoting" << nl;
QRMatrix<MatrixType> QRM
(
- A0,
- QRMatrix<MatrixType>::outputTypes::FULL_R,
- QRMatrix<MatrixType>::colPivoting::TRUE
+ A,
+ QRMatrix<MatrixType>::modes::ECONOMY,
+ QRMatrix<MatrixType>::outputs::BOTH_QR,
+ QRMatrix<MatrixType>::pivoting::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);
+ verify_pivoting(QRM, A);
}
- #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
+// Verify the Moore-Penrose pseudo-inverse function
+template<class MatrixType>
+void verify_pinv
+(
+ const MatrixType& A
+)
+{
+ const scalar absTol = 1e-6;
- #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
+ Info<< "## Verify Moore-Penrose pseudo-inverse: pinv()" << nl << endl;
+ const MatrixType APlus(MatrixTools::pinv(A));
- #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);
- }
+ #if PRINTALL
+ InfoNumPyFormat(A, "A");
+ InfoNumPyFormat(APlus, "A^+");
#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);
+ const MatrixType APlusPlus(MatrixTools::pinv(APlus));
+ cmp("# (A^+)^+ = A: ", flt(APlusPlus), flt(A), absTol);
}
- #endif
- #if (0 | RUNALL)
+ // The Moore-Penrose conditions
+ // w.wiki/5RL (Retrieved: 29-06-19)
{
- 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);
+ const MatrixType AAPlusA((A*APlus)*A);
+ cmp("# A*(A^+)*A = A: ", flt(AAPlusA), flt(A), absTol);
}
- #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);
+ const MatrixType APlusAAPlus((APlus*A)*APlus);
+ cmp("# (A^+)*A*(A^+) = A: ", flt(APlusAAPlus), flt(APlus), absTol);
}
- #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);
+ const MatrixType AAPlus(A*APlus);
+ const MatrixType AAPlusT(AAPlus.T());
+ cmp("# (A*(A^+)).T() = A*(A^+): ", flt(AAPlusT), flt(AAPlus), absTol);
}
- #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);
+ const MatrixType APlusA(APlus*A);
+ const MatrixType APlusAT(APlusA.T());
+ cmp("# ((A^+)*A = ((A^+)*A).T(): ", flt(APlusA), flt(APlusAT), absTol);
}
- #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
+// Verify the direct solution of a linear system by QRMatrix
+void verify_solve
+(
+ const SquareMatrix<complex>& A
+)
+{
+ // do nothing
}
-// Checks the direct solution of a linear system by QRMatrix
-template<class MatrixType>
-void verification_QRMatrixSolve
+// Verify the direct solution of a linear system by QRMatrix
+void verify_solve
(
- MatrixType& A
+ const SquareMatrix<scalar>& A
)
{
- typedef SquareMatrix<typename MatrixType::cmptType> SMatrix;
+ Info<< "## Verify direct solution of A*x = b with A = Q*R:" << nl << endl;
- // Create working copies of matrix A
- const MatrixType Aorig(A);
+ typedef SquareMatrix<scalar> SMatrix;
+ const scalar absTol = 1e-6;
// Direct solution of a linear system by QR decomposition
- #if (0 | RUNALL)
- {
- MatrixType A0(A);
- QRMatrix<MatrixType> QRM
- (
- A0,
- QRMatrix<MatrixType>::outputTypes::FULL_QR
- );
+ QRMatrix<SMatrix> QRM
+ (
+ A,
+ QRMatrix<SMatrix>::modes::FULL,
+ QRMatrix<SMatrix>::outputs::BOTH_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 residual(A.m(), 0);
+ const scalarField source(A.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);
+ const scalarField x(QRM.solve(source));
+ const scalarField AxMinusSource(A*x - source);
+ const scalarField xMinusQRinvSource(x - QRM.inv()*source);
+ const SMatrix QRinvA(QRM.inv()*A);
+ const SMatrix identity(A.m(), I);
- forAll(residual, i)
- {
- isEqual(A0xMinusSource[i], residual[i], verbose);
- isEqual(xMinusQRinvSource[i], residual[i], verbose);
- }
- equal(QRinvA0, identity, verbose);
- }
- #endif
+ cmp("# A*x - b = residual: ", AxMinusSource, residual, absTol);
+ cmp("# x - inv(QR)*b = residual: ", xMinusQRinvSource, residual, absTol);
+ cmp("# inv(QR)*A = I: ", flt(QRinvA), flt(identity), absTol);
}
-// Checks the parallel tall-skinny QR decomposition
+// Verify the parallel tall-skinny QR decomposition
template<class Type>
-void verification_tsqr
+void verify_tsqr
(
Random& rndGen
)
{
+ Info<< "## Verify parallel direct tall-skinny QR decomposition:"
+ << nl << endl;
+
typedef RectangularMatrix<Type> RMatrix;
// Size of the full matrix and its partitions
@@ -978,7 +370,7 @@ void verification_tsqr
mRowsSum += mRows;
}
- # if 0
+ #if PRINTALL
Info<< "mRowsSum = " << mRowsSum << nl;
Info<< "nColsSum = " << nColsSum << nl;
Info<< "mRowsList = " << mRowsList << nl;
@@ -999,8 +391,9 @@ void verification_tsqr
QRMatrix<RMatrix> QRM
(
A,
- QRMatrix<RMatrix>::outputTypes::REDUCED_R,
- QRMatrix<RMatrix>::colPivoting::FALSE
+ QRMatrix<RMatrix>::modes::ECONOMY,
+ QRMatrix<RMatrix>::outputs::ONLY_R,
+ QRMatrix<RMatrix>::pivoting::FALSE
);
const RMatrix& R = QRM.R();
@@ -1023,8 +416,9 @@ void verification_tsqr
QRMatrix<RMatrix> subQRM
(
subA,
- QRMatrix<RMatrix>::outputTypes::REDUCED_R,
- QRMatrix<RMatrix>::colPivoting::FALSE
+ QRMatrix<RMatrix>::modes::ECONOMY,
+ QRMatrix<RMatrix>::outputs::ONLY_R,
+ QRMatrix<RMatrix>::pivoting::FALSE
);
const RMatrix& subR = subQRM.R();
@@ -1036,8 +430,9 @@ void verification_tsqr
QRMatrix<RMatrix> TSQR
(
subRs,
- QRMatrix<RMatrix>::outputTypes::REDUCED_R,
- QRMatrix<RMatrix>::colPivoting::FALSE
+ QRMatrix<RMatrix>::modes::ECONOMY,
+ QRMatrix<RMatrix>::outputs::ONLY_R,
+ QRMatrix<RMatrix>::pivoting::FALSE
);
const RMatrix& TSQRR = TSQR.R();
@@ -1048,433 +443,161 @@ void verification_tsqr
}
// Compare magnitude of R, since R is not unique up to the sign
- equal(RMag, TSQRRMag, verbose);
+ cmp("# TSQR: ", flt(RMag), flt(TSQRRMag));
- #if 0
+ #if PRINTALL
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, 10, 1e-3, 1e-3);
- }
- #endif
-}
-
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-// Checks the Moore-Penrose pseudo-inverse function
template<class MatrixType>
-void verification_pinv
-(
- const MatrixType& A
-)
+void test_constructors(MatrixType)
{
- 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);
+ Info<< "# Construct from a Matrix and modes" << nl;
+ MatrixType A({5,5}, Zero);
+ const QRMatrix<MatrixType> QRM
+ (
+ A,
+ QRMatrix<MatrixType>::modes::FULL
+ );
}
-
- 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);
+ Info<< "# Construct from a const Matrix and modes" << nl;
+ const MatrixType A({5,5}, Zero);
+ const QRMatrix<MatrixType> QRM
+ (
+ A,
+ QRMatrix<MatrixType>::modes::FULL
+ );
}
}
-// * * * * * * * * * * * * * * * Main Program * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
+template<class MatrixType>
+void test_decomposition(MatrixType)
{
- typedef SquareMatrix<scalar> SMatrix;
- typedef RectangularMatrix<scalar> RMatrix;
- typedef SquareMatrix<complex> SCMatrix;
- typedef RectangularMatrix<complex> RCMatrix;
+ typedef typename MatrixType::cmptType cmptType;
+ typedef SquareMatrix<cmptType> SMatrix;
+ typedef RectangularMatrix<cmptType> RMatrix;
- Info<< setprecision(15);
Random rndGen(1234);
- label numberOfTests = 10;
+ const label szTests = 20;
+ const label min = 10;
+ const label max = 50;
- Info<< "### QR Decomposition:" << nl << nl;
-
- // SquareMatrix<scalar>
- #if (0 | RUNALL)
{
- horizontalLine();
-
- Info<< "## " << numberOfTests << " SquareMatrix<scalar> A:" << nl;
-
- for (label i = 0; i < numberOfTests; ++i)
+ for (label i = 0; i < szTests; ++i)
{
- const label mRows = rndGen.position<label>(1, 50);
- Info<< nl << nl << "# Random A with random mRows = " << mRows << nl;
+ const label m = rndGen.position(min, max);
+ const label n = rndGen.position(min, max);
- SMatrix A(makeRandomMatrix<SMatrix>({mRows, mRows}, rndGen));
- #if (0 | RUNALL)
- {
- const SMatrix B(A);
- verification_pinv(B);
- }
- #endif
- verification_QRMatrix(A);
- verification_QRMatrixSolve(A);
+ const RMatrix A
+ (
+ makeRandomMatrix<RMatrix>({m,n}, rndGen)
+ );
+ verify_decomposition(A);
+ verify_pinv(A);
}
-
- horizontalLine();
}
- #endif
-
- // RectangularMatrix<scalar>
- #if (0 | RUNALL)
{
- horizontalLine();
-
- Info<< "## " << numberOfTests << " RectangularMatrix<scalar> A:" << nl;
-
- for (label i = 0; i < numberOfTests; ++i)
+ for (label i = 0; i < szTests; ++i)
{
- const label mRows = rndGen.position<label>(1, 50);
- const label nCols = rndGen.position<label>(1, 50);
- Info<< nl << nl << "# Random matrix A with"
- << " random mRows = " << mRows
- << " random nCols = " << nCols << nl;
-
- RMatrix A(makeRandomMatrix<RMatrix>({mRows, nCols}, rndGen));
+ const label m = rndGen.position(min, max);
- #if (0 | RUNALL)
- {
- const RMatrix B(A);
- verification_pinv(B);
- }
- #endif
+ const SMatrix A
+ (
+ makeRandomMatrix<SMatrix>({m,m}, rndGen)
+ );
- verification_QRMatrix(A);
+ verify_decomposition(A);
+ verify_pinv(A);
+ verify_solve(A);
}
-
- horizontalLine();
}
- #endif
-
- // SquareMatrix<complex(scalar, 0.0)>
- #if (0 | RUNALL)
+ for (label i = 0; i < szTests; ++i)
{
- horizontalLine();
-
- Info<< "## " << numberOfTests << " SquareMatrix<complex, 0> A:" << nl;
-
- for (label i = 0; i < numberOfTests; ++i)
- {
- const label mRows = rndGen.position<label>(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();
+ verify_tsqr<cmptType>(rndGen);
}
- #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<label>(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
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Do compile-time recursion over the given types
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I == sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID){}
- // RectangularMatrix<complex(scalar, scalar)>
- #if (0 | RUNALL)
- {
- horizontalLine();
- Info<< "## " << numberOfTests << " RectangularMatrix<complex> A:" << nl;
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I < sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID)
+{
+ Info<< nl << " ## Test constructors: "<< typeID[I] <<" ##" << nl;
+ test_constructors(std::get<I>(types));
- for (label i = 0; i < numberOfTests; ++i)
- {
- const label mRows = rndGen.position<label>(1, 50);
- const label nCols = rndGen.position<label>(1, 50);
- Info<< nl << nl << "# Random matrix A with"
- << " random mRows = " << mRows
- << " random nCols = " << nCols << nl;
+ Info<< nl << " ## Test decomposition: "<< typeID[I] <<" ##" << nl;
+ test_decomposition(std::get<I>(types));
- RCMatrix A(makeRandomMatrix<RCMatrix>({mRows, nCols}, rndGen));
+ run_tests<I + 1, Tp...>(types, typeID);
+}
- #if (0 | RUNALL)
- {
- const RCMatrix B(A);
- verification_pinv(B);
- }
- #endif
- verification_QRMatrix(A);
- }
+// * * * * * * * * * * * * * * * Main Program * * * * * * * * * * * * * * * //
- horizontalLine();
- }
- #endif
+int main()
+{
+ Info<< setprecision(15);
+ const std::tuple
+ <
+ RectangularMatrix<doubleScalar>,
+ RectangularMatrix<complex>,
+ SquareMatrix<doubleScalar>,
+ SquareMatrix<complex>
+ > types
+ (
+ std::make_tuple(Zero, Zero, Zero, Zero)
+ );
- #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);
- }
+ const List<word> typeID
+ ({
+ "RectangularMatrix<doubleScalar>",
+ "RectangularMatrix<complex>",
+ "SquareMatrix<doubleScalar>",
+ "SquareMatrix<complex>"
+ });
- Info<< nl << "## " << numberOfTests << " <complex> tests:" << nl;
- for (label i = 0; i < numberOfTests; ++i)
- {
- verification_tsqr<complex>(rndGen);
- }
+ std::chrono::steady_clock::time_point begin =
+ std::chrono::steady_clock::now();
- horizontalLine();
- }
- #endif
+ run_tests(types, typeID);
+ std::chrono::steady_clock::time_point end =
+ std::chrono::steady_clock::now();
+ Info<< "Time elapsed = "
+ << std::chrono::duration_cast<std::chrono::microseconds>(end - begin)
+ .count()
+ << "[µs]" << endl;
- #if (0 | RUNALL)
+ if (nFail_)
{
- 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<label>(20, 50);
- const label pCols = rndGen.position<label>(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<label>(20, 50);
- const label pCols = rndGen.position<label>(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
+ Info<< nl << " #### "
+ << "Failed in " << nFail_ << " tests "
+ << "out of total " << nTest_ << " tests "
+ << "####\n" << endl;
+ return 1;
}
- #endif
-
- Info<< nl << "End" << nl;
+ Info<< nl << " #### Passed all " << nTest_ <<" tests ####\n" << endl;
return 0;
}
diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.C b/src/OpenFOAM/matrices/Matrix/Matrix.C
index d27013e8aeca8d6d7d9113112bfb00630a804fc0..0ba6168888428ee34a199c22e2374a50c2b8ebd2 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.C
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.C
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
- Copyright (C) 2019-2021 OpenCFD Ltd.
+ Copyright (C) 2019-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -405,6 +405,23 @@ Form Foam::Matrix<Form, Type>::T() const
}
+template<class Form, class Type>
+Form Foam::Matrix<Form, Type>::transpose() const
+{
+ Form At(labelPair{n(), m()});
+
+ for (label i = 0; i < m(); ++i)
+ {
+ for (label j = 0; j < n(); ++j)
+ {
+ At(j, i) = (*this)(i, j);
+ }
+ }
+
+ return At;
+}
+
+
template<class Form, class Type>
Foam::List<Type> Foam::Matrix<Form, Type>::diag() const
{
@@ -1003,9 +1020,9 @@ operator&
Zero
);
- for (label i = 0; i < AB.m(); ++i)
+ for (label k = 0; k < B.m(); ++k)
{
- for (label k = 0; k < B.m(); ++k)
+ for (label i = 0; i < AB.m(); ++i)
{
for (label j = 0; j < AB.n(); ++j)
{
@@ -1048,9 +1065,9 @@ operator^
for (label i = 0; i < AB.m(); ++i)
{
- for (label k = 0; k < BT.n(); ++k)
+ for (label j = 0; j < AB.n(); ++j)
{
- for (label j = 0; j < AB.n(); ++j)
+ for (label k = 0; k < BT.n(); ++k)
{
AB(i, j) += A(i, k)*Detail::conj(BT(j, k));
}
diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.H b/src/OpenFOAM/matrices/Matrix/Matrix.H
index c11d7a27ec0f544a4abe26526f2ea608dab4f53c..eb34e064d31dd36cace4fd1fc451f6fa35407e24 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.H
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.H
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
- Copyright (C) 2019-2021 OpenCFD Ltd.
+ Copyright (C) 2019-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -372,9 +372,12 @@ public:
// Operations
- //- Return (conjugate) transpose of Matrix
+ //- Return conjugate transpose of Matrix
Form T() const;
+ //- Return non-conjugate transpose of Matrix
+ Form transpose() const;
+
//- Right-multiply Matrix by a column vector (A * x)
inline tmp<Field<Type>> Amul(const UList<Type>& x) const;
diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
index 0b63f2017b03ed289be06edb2e6a381fd93e664d..f7bf7db3ebfa98f9b7b7360af7375661aba7eeae 100644
--- a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
+++ b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016 OpenFOAM Foundation
- Copyright (C) 2019-2021 OpenCFD Ltd.
+ Copyright (C) 2019-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -31,96 +31,292 @@ License
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class MatrixType>
-void Foam::QRMatrix<MatrixType>::qr
+Foam::label Foam::QRMatrix<MatrixType>::calcShapeFactor
(
- MatrixType& A
+ const MatrixType& A
+) const
+{
+ if (mode_ == modes::FULL)
+ {
+ return A.m();
+ }
+ else if (mode_ == modes::ECONOMY)
+ {
+ return min(A.m(), A.n());
+ }
+
+ return 0;
+}
+
+
+template<class MatrixType>
+void Foam::QRMatrix<MatrixType>::decompose
+(
+ MatrixType& AT
)
{
- const label nIter = min(A.m() - 1, A.n());
+ const label n = AT.m();
+ const label m = AT.n();
+
+ List<cmptType> Rdiag(n, Zero);
- // Loop through all subcolumns of which the diagonal elem is the first elem
- for (label k = 0; k < nIter; ++k)
+ for (label k = 0; k < n; ++k)
{
- const RMatrix u(householderReflector(A.subColumn(k, k)));
+ // Compute 2-norm of k-th column without under/overflow
+ scalar nrm = 0;
- applyHouseholder(A, u, k);
- }
+ for (label i = k; i < m; ++i)
+ {
+ nrm = Foam::hypot(nrm, mag(AT(k,i)));
+ }
- if (outputType_ == outputTypes::REDUCED_R)
- {
- A.resize(A.n(), A.n());
+ // Avoid divide-by-zero. Use compare with == 0 not mag or VSMALL etc,
+ // otherwise wrong results (may need more investigation)
+ if (nrm != scalar(0))
+ {
+ // Form k-th Householder vector
+ if (mag(AT(k,k)) < 0)
+ {
+ nrm *= -1;
+ }
+
+ for (label i = k; i < m; ++i)
+ {
+ AT(k,i) /= nrm;
+ }
+
+ AT(k,k) += scalar(1);
+
+ // Apply transformation to remaining columns
+ for (label j = k + 1; j < n; ++j)
+ {
+ cmptType s(Zero);
+
+ for (label i = k; i < m; ++i)
+ {
+ s += Detail::conj(AT(k,i))*AT(j,i);
+ }
+
+ if (mag(AT(k,k)) > SMALL)
+ {
+ s /= -AT(k,k);
+ }
+
+ for (label i = k; i < m; ++i)
+ {
+ AT(j,i) += s*AT(k,i);
+ }
+ }
+ }
+
+ Rdiag[k] = -nrm;
}
+
+ calcQ(AT);
+
+ calcR(AT, Rdiag);
}
template<class MatrixType>
-void Foam::QRMatrix<MatrixType>::qrPivot
+void Foam::QRMatrix<MatrixType>::decompose
(
- MatrixType& A
+ MatrixType& AT,
+ const bool pivoting
)
{
- const label nCols = A.n();
- const label nIter = min(A.m() - 1, A.n());
+ const label n = AT.m();
+ const label m = AT.n();
+ const label sz = min(m,n);
// Initialise permutation vector, and column norm vector
- P_ = identity(nCols);
+ p_ = identity(n);
// Initialise vector norms of each column of A
- List<scalar> colNorms(nCols);
- for (label k = 0; k < nCols; ++k)
+ List<scalar> norms(n, Zero);
+ for (label k = 0; k < n; ++k)
{
- colNorms[k] = A.columnNorm(k, true);
+ for (label i = 0; i < m; ++i)
+ {
+ norms[k] += magSqr(AT(k, i));
+ }
}
- // Loop through all subcolumns of which the diagonal elem is the first elem
- for (label k = 0; k < nIter; ++k)
- {
- const labelRange colRange(k, nCols);
- const SubList<scalar> subColNorms(colNorms, colRange);
+ List<cmptType> Rdiag(n, Zero);
- // Column pivoting
- const label maxColNormi =
+ for (label k = 0; k < sz; ++k)
+ {
+ const auto it = std::next(norms.cbegin(), k);
+ const label maxNormi =
std::distance
(
- subColNorms.cbegin(),
- std::max_element(subColNorms.cbegin(), subColNorms.cend())
+ it,
+ std::max_element(it, norms.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)))
+ // Swap A, P and norms according to pivot column if the current
+ // column is not the max norm column. Also avoid any column swaps
+ // where the leading element is very small
+ if (mag(AT(k + maxNormi,k)) > SMALL && maxNormi != 0)
{
- const RMatrix R1(A.subColumn(k));
- const RMatrix R2(A.subColumn(maxColNormi + k));
- A.subColumn(k) = R2;
- A.subColumn(maxColNormi + k) = R1;
+ const RMatrix R1(AT.subRow(k));
+ const RMatrix R2(AT.subRow(maxNormi + k));
+
+ AT.subRow(k) = R2;
+ AT.subRow(maxNormi + k) = R1;
- Swap(P_[k], P_[maxColNormi + k]);
- Swap(colNorms[k], colNorms[maxColNormi + k]);
+ Swap(p_[k], p_[maxNormi + k]);
+ Swap(norms[k], norms[maxNormi + k]);
}
{
- const RMatrix u(householderReflector(A.subColumn(k, k)));
+ // Compute 2-norm of k-th column without under/overflow
+ scalar nrm = 0;
- applyHouseholder(A, u, k);
+ for (label i = k; i < m; ++i)
+ {
+ nrm = Foam::hypot(nrm, mag(AT(k,i)));
+ }
+
+ if (nrm != scalar(0))
+ {
+ // Form k-th Householder vector
+ if (mag(AT(k,k)) < 0)
+ {
+ nrm *= -1;
+ }
+
+ for (label i = k; i < m; ++i)
+ {
+ AT(k,i) /= nrm;
+ }
+
+ AT(k,k) += scalar(1);
+
+ // Apply transformation to remaining columns
+ for (label j = k + 1; j < n; ++j)
+ {
+ cmptType s(Zero);
+
+ for (label i = k; i < m; ++i)
+ {
+ s += Detail::conj(AT(k,i))*AT(j,i);
+ }
+
+ if (mag(AT(k,k)) > SMALL)
+ {
+ s /= -AT(k,k);
+ }
+
+ for (label i = k; i < m; ++i)
+ {
+ AT(j,i) += s*AT(k,i);
+ }
+ }
+ }
+
+ Rdiag[k] = -nrm;
}
// Update norms
- if (k < nIter - 1)
+ if (k < sz - 1)
{
label q = k + 1;
- for (const auto& val : RMatrix(A.subRow(k, q)))
+ for (const auto& val : RMatrix(AT.subColumn(k, q)))
{
- colNorms[q] -= magSqr(val);
+ norms[q] -= magSqr(val);
++q;
}
}
}
- if (outputType_ == outputTypes::REDUCED_R)
+ calcQ(AT);
+
+ calcR(AT, Rdiag);
+}
+
+
+template<class MatrixType>
+void Foam::QRMatrix<MatrixType>::calcQ
+(
+ const MatrixType& AT
+)
+{
+ if (output_ == outputs::ONLY_R)
{
- A.resize(A.n(), A.n());
+ return;
+ }
+
+ const label n = AT.m();
+ const label m = AT.n();
+
+ Q_.resize(m, sz_);
+
+ MatrixType QT(Q_.transpose());
+
+ for (label k = sz_ - 1; k >= 0; --k)
+ {
+ for (label i = 0; i < m; ++i)
+ {
+ QT(k,i) = Zero;
+ }
+
+ QT(k,k) = scalar(1);
+
+ for (label j = k; j < sz_; ++j)
+ {
+ if (n > k && mag(AT(k,k)) > SMALL)
+ {
+ cmptType s(Zero);
+
+ for (label i = k; i < m; ++i)
+ {
+ s += AT(k,i)*QT(j,i);
+ }
+
+ s /= -AT(k,k);
+
+ for (label i = k; i < m; ++i)
+ {
+ QT(j,i) += s*Detail::conj(AT(k,i));
+ }
+ }
+ }
+ }
+
+ Q_ = QT.T();
+}
+
+
+template<class MatrixType>
+void Foam::QRMatrix<MatrixType>::calcR
+(
+ const MatrixType& AT,
+ const List<cmptType>& diag
+)
+{
+ if (output_ == outputs::ONLY_Q)
+ {
+ return;
+ }
+
+ const label n = AT.m();
+
+ R_.resize(sz_, n);
+
+ for (label i = 0; i < sz_; ++i)
+ {
+ for (label j = 0; j < n; ++j)
+ {
+ if (i < j)
+ {
+ R_(i,j) = AT(j,i);
+ }
+ else if (i == j)
+ {
+ R_(i,j) = diag[i];
+ }
+ }
}
}
@@ -182,49 +378,36 @@ 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
+ const modes mode,
+ const outputs output,
+ const bool pivoting,
+ MatrixType& A
)
:
- outputType_(outputType),
- storeMethod_(storeMethod),
- colPivot_(colPivot),
+ mode_(mode),
+ output_(output),
+ sz_(calcShapeFactor(A)),
Q_(),
R_(),
- P_()
-{}
+ p_()
+{
+ // Non-conjugate transpose of input matrix
+ MatrixType AT(A.transpose());
+ A.clear();
+ if (pivoting)
+ {
+ decompose(AT, pivoting);
+ }
+ else
+ {
+ decompose(AT);
+ }
-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);
+ AT.clear();
}
@@ -232,126 +415,34 @@ template<class MatrixType>
Foam::QRMatrix<MatrixType>::QRMatrix
(
const MatrixType& A,
- const outputTypes outputType,
- const colPivoting colPivot
+ const modes mode,
+ const outputs output,
+ const bool pivoting
)
:
- outputType_(outputType),
- storeMethod_(storeMethods::OUT_OF_PLACE),
- colPivot_(colPivot),
+ mode_(mode),
+ output_(output),
+ sz_(calcShapeFactor(A)),
Q_(),
R_(),
- P_()
+ p_()
{
- decompose(A);
-}
-
+ MatrixType AT(A.transpose());
-// * * * * * * * * * * * * * * 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)
+ if (pivoting)
{
- 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
+ decompose(AT, pivoting);
}
-
- // Allocate resources for Q_, if need be
- if (outputType_ == outputTypes::FULL_QR)
+ else
{
- Q_ = MatrixType({A.m(), A.m()}, I);
+ decompose(AT);
}
- // 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;
- }
- }
+ AT.clear();
}
-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_);
- }
-}
-
+// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
template<class MatrixType>
void Foam::QRMatrix<MatrixType>::solve
@@ -407,39 +498,15 @@ Foam::QRMatrix<MatrixType>::solve
}
-template<class MatrixType>
-typename Foam::QRMatrix<MatrixType>::SMatrix
-Foam::QRMatrix<MatrixType>::inv() const
-{
- const label m = Q_.m();
-
- Field<cmptType> x(m);
- SMatrix inv(m);
-
- for (label i = 0; i < m; ++i)
- {
- for (label j = 0; j < m; ++j)
- {
- x[j] = Q_(i, j);
- }
- solvex(x);
- inv.subColumn(i) = x;
- }
-
- return inv;
-}
-
-
template<class MatrixType>
Foam::RectangularMatrix<typename MatrixType::cmptType>
-Foam::QRMatrix<MatrixType>::backSubstitution
+Foam::QRMatrix<MatrixType>::solve
(
- const SMatrix& A,
const RMatrix& rhs
)
{
- const label mRows = A.m();
- const label nCols = A.n();
+ const label mRows = R_.m();
+ const label nCols = R_.n();
const label pCols = rhs.n();
#ifdef FULLDEBUG
@@ -453,19 +520,20 @@ Foam::QRMatrix<MatrixType>::backSubstitution
<< abort(FatalError);
}
- const List<cmptType> diag(A.diag());
+ const List<cmptType> diag(R_.diag());
for (const auto& val : diag)
{
if (mag(val) < SMALL)
{
WarningInFunction
- << "SMALL-valued diagonal elem in back-substitution." << nl;
+ << "SMALL-valued diagonal elem in back-substitution."
+ << endl;
}
}
}
#endif
- RMatrix X({nCols, pCols}, Zero);
+ RMatrix b({nCols, pCols}, Zero);
for (label i = 0; i < pCols; ++i)
{
@@ -475,27 +543,49 @@ Foam::QRMatrix<MatrixType>::backSubstitution
for (label k = j + 1; k < mRows; ++k)
{
- alpha -= X(k, i)*A(j, k);
+ alpha -= b(k, i)*R_(j, k);
}
- X(j, i) = alpha/A(j, j);
+ b(j, i) = alpha/R_(j, j);
+ }
+ }
+
+ return b;
+}
+
+
+template<class MatrixType>
+typename Foam::QRMatrix<MatrixType>::SMatrix
+Foam::QRMatrix<MatrixType>::inv() const
+{
+ const label m = Q_.m();
+
+ Field<cmptType> x(m);
+ SMatrix inv(m);
+
+ for (label i = 0; i < m; ++i)
+ {
+ for (label j = 0; j < m; ++j)
+ {
+ x[j] = Q_(i, j);
}
+ solvex(x);
+ inv.subColumn(i) = x;
}
- return X;
+ return inv;
}
// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
template<class MatrixType>
-MatrixType Foam::pinv
+MatrixType Foam::MatrixTools::pinv
(
const MatrixType& A,
- const scalar tolerance
+ scalar tol
)
{
- scalar tol = tolerance;
typedef typename MatrixType::cmptType cmptType;
if (A.empty())
@@ -518,15 +608,16 @@ MatrixType Foam::pinv
QRMatrix<MatrixType> QRM
(
A,
- QRMatrix<MatrixType>::outputTypes::FULL_QR,
- QRMatrix<MatrixType>::colPivoting::TRUE
+ QRMatrix<MatrixType>::modes::FULL,
+ QRMatrix<MatrixType>::outputs::BOTH_QR,
+ QRMatrix<MatrixType>::pivoting::TRUE
);
- const MatrixType& R(QRM.R());
- const MatrixType& Q(QRM.Q());
+ 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;
+ label elemi = 0;
label i = 0;
while (i < 2)
{
@@ -534,14 +625,14 @@ MatrixType Foam::pinv
auto lessThan = [=](const cmptType& x) { return tol > mag(x); };
- firstZeroElemi =
+ elemi =
std::distance
(
diag.cbegin(),
std::find_if(diag.cbegin(), diag.cend(), lessThan)
);
- if (firstZeroElemi == 0)
+ if (elemi == 0)
{
if (i == 0)
{
@@ -572,7 +663,7 @@ MatrixType Foam::pinv
// 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));
+ const RectangularMatrix<cmptType> R1(R.subMatrix(0, 0, elemi));
// R2 (KP:p. 648)
if (R1.n() < R1.m())
@@ -582,15 +673,15 @@ MatrixType Foam::pinv
QRMatrix<SquareMatrix<cmptType>> QRSolve
(
C,
- QRMatrix<SquareMatrix<cmptType>>::outputTypes::FULL_QR
+ QRMatrix<SquareMatrix<cmptType>>::modes::FULL,
+ QRMatrix<SquareMatrix<cmptType>>::outputs::BOTH_QR
);
RectangularMatrix<cmptType> R2
(
- QRSolve.backSubstitution
+ QRSolve.solve
(
- QRSolve.R(),
- QRSolve.Q() & (R1.T())
+ QRSolve.Q() & R1.T()
)
);
@@ -606,14 +697,14 @@ MatrixType Foam::pinv
QRMatrix<SquareMatrix<cmptType>> QRSolve
(
C,
- QRMatrix<SquareMatrix<cmptType>>::outputTypes::FULL_QR
+ QRMatrix<SquareMatrix<cmptType>>::modes::FULL,
+ QRMatrix<SquareMatrix<cmptType>>::outputs::BOTH_QR
);
RectangularMatrix<cmptType> R2
(
- QRSolve.backSubstitution
+ QRSolve.solve
(
- QRSolve.R(),
QRSolve.Q() & R1
)
);
diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H
index ee9bd62edf2ae4328c618aeb6083543a52050fde..dcb32cc6db14207853b1559c6c0bbca26b4ed71c 100644
--- a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H
+++ b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016 OpenFOAM Foundation
- Copyright (C) 2019-2021 OpenCFD Ltd.
+ Copyright (C) 2019-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -28,29 +28,36 @@ Class
Foam::QRMatrix
Description
- QRMatrix (i.e. QR decomposition, QR factorisation or orthogonal-triangular
- decomposition) decomposes a scalar/complex matrix \c A into the following
- matrix product:
+ \c QRMatrix computes QR decomposition of a given
+ scalar/complex matrix \c A into the following:
\verbatim
- A = Q*R,
+ A = Q R
\endverbatim
- where
- \c Q is a unitary similarity matrix,
- \c R is an upper triangular matrix.
+ or in case of QR decomposition with column pivoting:
- Reference:
\verbatim
- QR decomposition:
- mathworld.wolfram.com/QRDecomposition.html (Retrieved:17-06-19)
- w.wiki/52X (Retrieved: 17-06-19)
+ A P = Q R
+ \endverbatim
- QR decomposition with Householder reflector (tag:M):
- Monahan, J. F. (2011).
- Numerical methods of statistics.
- Cambridge: Cambridge University Press.
- DOI:10.1017/CBO9780511977176
+ where
+ \vartable
+ Q | Unitary/orthogonal matrix
+ R | Upper triangular matrix
+ P | Permutation matrix
+ \endvartable
+
+ References:
+ \verbatim
+ TNT implementation:
+ Pozo, R. (1997).
+ Template Numerical Toolkit for linear algebra:
+ High performance programming with C++
+ and the Standard Template Library.
+ The International Journal of Supercomputer Applications
+ and High Performance Computing, 11(3), 251-263.
+ DOI:10.1177/109434209701100307
QR decomposition with column pivoting (tag:QSB):
Quintana-OrtÃ, G., Sun, X., & Bischof, C. H. (1998).
@@ -75,54 +82,47 @@ Description
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
+ \heading Input:
+ \vartable
+ A | RectangularMatrix<Type> or SquareMatrix<Type>
+ \endvartable
+
+ \heading Options for the decomposition mode:
+ \vartable
+ modes::FULL | compute full-size decomposition
+ modes::ECONOMY | compute economy-size decomposition
+ \endvartable
+
+ \heading Options for the output types:
+ \vartable
+ outputs::ONLY\_Q | compute only Q
+ outputs::ONLY\_R | compute only R
+ outputs::BOTH\_QR | compute both Q and R
+ \endvartable
+
+ \heading Options for the column pivoting:
+ \vartable
+ pivoting::FALSE | switch off column pivoting
+ pivoting::TRUE | switch on column pivoting
+ \endvartable
+
+ \heading Output:
+ \vartable
+ Q | m-by-m (FULL) or m-by-k (ECONOMY) with k = min(m,n)
+ R | m-by-n (FULL) or k-by-n (ECONOMY) with k = min(m,n)
+ p | n-element label list
+ P | n-by-n permutation matrix
+ \endvartable
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.
+ - \c QRMatrix involves modified implementations of the public-domain
+ library \c TNT, which is available at https://math.nist.gov/tnt/index.html.
+ - \c Q and \c R are always of the same type of the matrix \c A.
+ - \c Type can be \c scalar or \c complex.
See also
Test-QRMatrix.C
@@ -133,12 +133,11 @@ SourceFiles
\*---------------------------------------------------------------------------*/
-#ifndef QRMatrix_H
-#define QRMatrix_H
+#ifndef Foam_QRMatrix_H
+#define Foam_QRMatrix_H
#include "RectangularMatrix.H"
#include "SquareMatrix.H"
-#include "complex.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -146,39 +145,38 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
- Class QRMatrix Declaration
+ Class QRMatrix Declaration
\*---------------------------------------------------------------------------*/
template<class MatrixType>
class QRMatrix
{
-
public:
typedef typename MatrixType::cmptType cmptType;
typedef SquareMatrix<cmptType> SMatrix;
typedef RectangularMatrix<cmptType> RMatrix;
- //- Options for the output matrix forms of QRMatrix
- enum outputTypes : uint8_t
+ //- Options for the decomposition mode
+ enum modes : 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
+ FULL = 1, //!< compute full-size decomposition
+ ECONOMY = 2, //!< compute economy-size decomposition
};
- //- Options where to store R
- enum storeMethods : uint8_t
+ //- Options for the output types
+ enum outputs : uint8_t
{
- IN_PLACE = 1, //!< replaces input matrix content with \c R
- OUT_OF_PLACE = 2 //!< creates new object of \c R
+ ONLY_Q = 1, //!< compute only Q
+ ONLY_R = 2, //!< compute only R
+ BOTH_QR = 3 //!< compute both Q and R
};
- //- Options for the computation of column pivoting
- enum colPivoting : bool
+ //- Options for the column pivoting
+ enum pivoting : bool
{
- FALSE = false, //!< switches off column pivoting
- TRUE = true //!< switches on column pivoting
+ FALSE = false, //!< switch off column pivoting
+ TRUE = true //!< switch on column pivoting
};
@@ -186,53 +184,50 @@ private:
// Private Data
- //- Selected option for the output matrix forms of QRMatrix
- outputTypes outputType_;
+ //- Decomposition mode
+ const modes mode_;
- //- Selected option where to store R
- const storeMethods storeMethod_;
+ //- Output type
+ const outputs output_;
- //- Selected option for the computation of column pivoting
- const colPivoting colPivot_;
+ //- Output shape factor based on the decomposition mode
+ const label sz_;
- //- Unitary similarity matrix
+ //- Unitary/orthogonal matrix
MatrixType Q_;
//- Upper triangular matrix
MatrixType R_;
//- Permutation list (if column-pivoting is on)
- labelList P_;
+ labelList p_;
// Private Member Functions
- //- Compute the QR decomposition without the column pivoting
- void qr
- (
- MatrixType& A
- );
+ // Evaluation
- //- Compute the QR decomposition with the column pivoting
- // (QSB:Fig. 1)
- void qrPivot
- (
- MatrixType& A
- );
+ //- Calculate output shape factor based on the decomposition mode
+ label calcShapeFactor(const MatrixType& A) const;
- //- Compute output matrices in selected forms
- //- using Householder reflectors
- // (M:Section 4)
- inline void applyHouseholder
- (
- MatrixType& A,
- const RMatrix& reflector,
- const label k = 0
- );
+ //- Compute QR decomposition
+ void decompose(MatrixType& AT);
- //- 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
+ //- Compute QR decomposition with column pivoting
+ void decompose(MatrixType& AT, const bool pivot);
+
+ //- Calculate Q
+ void calcQ(const MatrixType& AT);
+
+ //- Calculate R
+ void calcR(const MatrixType& AT, const List<cmptType>& diag);
+
+
+ // Linear system solution
+
+ //- 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
template<template<typename> class ListContainer>
void solvex(ListContainer<cmptType>& x) const;
@@ -246,33 +241,11 @@ private:
) const;
-protected:
+ //- No copy construct
+ QRMatrix(const QRMatrix&) = delete;
- // 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
- );
+ //- No copy assignment
+ QRMatrix& operator=(const QRMatrix&) = delete;
public:
@@ -280,65 +253,66 @@ public:
// Constructors
//- Construct null
- QRMatrix();
+ QRMatrix() = delete;
- //- Construct QRMatrix without performing the decomposition
- QRMatrix
+ //- Construct with a matrix and perform QR decomposition
+ explicit 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
+ const modes mode,
+ const outputs output,
+ const bool pivoting,
+ MatrixType& A
);
- //- Construct QRMatrix and perform the QR decomposition
- QRMatrix
+ //- Construct with a const matrix and perform QR decomposition
+ explicit QRMatrix
(
const MatrixType& A,
- const outputTypes outputType,
- const colPivoting colPivot = colPivoting::FALSE
+ const modes mode,
+ const outputs output = outputs::BOTH_QR,
+ const bool pivoting = false
);
// Member Functions
- // 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;
+ // Access
+
+ //- Return const reference to Q
+ const MatrixType& Q() const noexcept
+ {
+ return Q_;
+ }
+
+ //- Return reference to Q
+ MatrixType& Q() noexcept
+ {
+ return Q_;
+ }
+
+ //- Return const reference to R
+ const MatrixType& R() const noexcept
+ {
+ return R_;
+ }
+
+ //- Return reference to R
+ MatrixType& R() noexcept
+ {
+ return R_;
+ }
+
+ //- Return const reference to p
+ const labelList& p() const noexcept
+ {
+ return p_;
+ }
//- Create and return the permutation matrix
inline SMatrix P() const;
- // Algorithm
-
- //- Compute QR decomposition according to constructor settings
- void decompose
- (
- MatrixType& A
- );
-
- //- Compute QR decomposition according to constructor settings
- void decompose
- (
- const MatrixType& A
- );
+ // Linear system solution
//- Solve the linear system with the given source
//- and return the solution in the argument x
@@ -372,41 +346,39 @@ public:
const IndirectListBase<cmptType, Addr>& source
) const;
- //- Return the inverse of (Q*R), so that solving x = (Q*R).inv()*source
- 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
+ //- Solve a row-echelon-form linear system (Ax = b)
+ //- starting from the bottom by back substitution
+ // A = R: Non-singular upper-triangular square matrix (m-by-m)
+ // b: Source (m-by-p)
+ // x: Solution (m-by-p)
+ RMatrix solve
(
- const SMatrix& A,
const RMatrix& b
);
+
+ //- Return the inverse of (Q*R), solving x = (Q*R).inv()*source
+ SMatrix inv() const;
};
+namespace MatrixTools
+{
// * * * * * * * * * * * * * * 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.
+//- Moore-Penrose inverse of singular/non-singular square/rectangular
+//- scalar/complex matrices (KPP:p. 9834; KP:p. 648)
+// The tolerance to ensure the R1 matrix full-rank is set to 1e-5
+// by (TA; mentioned in (KPP:p. 9832)) in contrast to 1e-13 (KPP:p. 9834).
template<class MatrixType>
MatrixType pinv
(
const MatrixType& A,
- const scalar tol = 1e-5
+ scalar tol = 1e-5
);
-
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+} // End namespace MatrixTools
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H b/src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H
index 929f8b6fc89f986388923bf6f51d52718343de9f..bf61a6b4e908908076f35a3eee039f40ed28c917 100644
--- a/src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H
+++ b/src/OpenFOAM/matrices/QRMatrix/QRMatrixI.H
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016 OpenFOAM Foundation
- Copyright (C) 2019 OpenCFD Ltd.
+ Copyright (C) 2019-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -26,166 +26,20 @@ License
\*---------------------------------------------------------------------------*/
-// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
-
-template<class MatrixType>
-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 MatrixType& Foam::QRMatrix<MatrixType>::Q() const
-{
- const_cast<MatrixType&>(Q_).round();
-
- return Q_;
-}
-
-
-template<class MatrixType>
-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>
+inline Foam::SquareMatrix<typename MatrixType::cmptType>
Foam::QRMatrix<MatrixType>::P() const
{
- SquareMatrix<cmptType> permMat(P_.size(), Zero);
+ SquareMatrix<cmptType> P(p_.size(), Zero);
- forAll(P_, jcol)
+ forAll(p_, jcol)
{
- permMat(P_[jcol], jcol) = pTraits<cmptType>::one;
+ P(p_[jcol], jcol) = pTraits<cmptType>::one;
}
- return permMat;
+ return P;
}
diff --git a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H
index 367f6c5d10a95ea095465804be0453c505ef2476..3e90bff571523f5f49e1537df847091296a2a851 100644
--- a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H
+++ b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
- Copyright (C) 2019-2020 OpenCFD Ltd.
+ Copyright (C) 2019-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -135,6 +135,9 @@ public:
// Member Operators
+ //- Move assignment
+ inline void operator=(RectangularMatrix<Type>&& mat);
+
//- Assign all elements to zero
inline void operator=(const Foam::zero);
diff --git a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H
index e3164483107c1b4ed9ed76a79ce1afb504f1b233..e4f67dbbb74dbe74dddcfa7d02b1c3d608632f8d 100644
--- a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H
+++ b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
- Copyright (C) 2019-2020 OpenCFD Ltd.
+ Copyright (C) 2019-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -171,6 +171,16 @@ Foam::RectangularMatrix<Type>::clone() const
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
+template<class Type>
+inline void Foam::RectangularMatrix<Type>::operator=
+(
+ RectangularMatrix<Type>&& mat
+)
+{
+ this->transfer(mat);
+}
+
+
template<class Type>
inline void Foam::RectangularMatrix<Type>::operator=(const Foam::zero)
{
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
index 8937cfc775b10d667c9c11df6a64d0f83d5a93cd..fc0a7313e4b3dfbf2813c29693404429cf1f40a8 100644
--- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
+++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
- Copyright (C) 2019-2020 OpenCFD Ltd.
+ Copyright (C) 2019-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -181,6 +181,9 @@ public:
// Member Operators
+ //- Move assignment
+ inline void operator=(SquareMatrix<Type>&& mat);
+
//- Assign all elements to zero
inline void operator=(const Foam::zero);
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
index f72d26fe775c9f4ea5ceb36862c77b863bd59561..42cd6666eac9a627b9223ef67a24784268fdd9e0 100644
--- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
+++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
- Copyright (C) 2019-2020 OpenCFD Ltd.
+ Copyright (C) 2019-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -290,6 +290,13 @@ inline bool Foam::SquareMatrix<Type>::tridiagonal() const
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
+template<class Type>
+inline void Foam::SquareMatrix<Type>::operator=(SquareMatrix<Type>&& mat)
+{
+ this->transfer(mat);
+}
+
+
template<class Type>
inline void Foam::SquareMatrix<Type>::operator=(const Foam::zero)
{
diff --git a/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMD.C b/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMD.C
index bbbdf4e4da8361c99ea70bef64db197e4fa22ad3..7f26f6927f04a2807ccffe384e18597d5992c0f1 100644
--- a/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMD.C
+++ b/src/functionObjects/field/DMD/DMDModels/derived/STDMD/STDMD.C
@@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
- Copyright (C) 2020-2021 OpenCFD Ltd.
+ Copyright (C) 2020-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -159,10 +159,11 @@ reducedKoopmanOperator()
QRMatrix<RMatrix> QRM
(
Qupper_,
- QRMatrix<RMatrix>::outputTypes::REDUCED_R,
- QRMatrix<RMatrix>::colPivoting::FALSE
+ QRMatrix<RMatrix>::modes::ECONOMY,
+ QRMatrix<RMatrix>::outputs::ONLY_R
);
- Rx = QRM.R();
+ RMatrix& R = QRM.R();
+ Rx.transfer(R);
}
Rx.round();
@@ -228,14 +229,18 @@ reducedKoopmanOperator()
}
}
- // Apply interim QR decomposition on Rx of receiver subset masters
- QRMatrix<RMatrix> QRM
- (
- Rx,
- QRMatrix<RMatrix>::outputTypes::REDUCED_R,
- QRMatrix<RMatrix>::storeMethods::IN_PLACE,
- QRMatrix<RMatrix>::colPivoting::FALSE
- );
+ {
+ // Apply interim QR decomposition
+ // on Rx of receiver subset masters
+ QRMatrix<RMatrix> QRM
+ (
+ Rx,
+ QRMatrix<RMatrix>::modes::ECONOMY,
+ QRMatrix<RMatrix>::outputs::ONLY_R
+ );
+ RMatrix& R = QRM.R();
+ Rx.transfer(R);
+ }
Rx.round();
}
@@ -281,26 +286,23 @@ reducedKoopmanOperator()
}
Info<< tab << "Computing TSQR" << endl;
- QRMatrix<RMatrix> QRM
- (
- Rx,
- QRMatrix<RMatrix>::outputTypes::REDUCED_R,
- QRMatrix<RMatrix>::storeMethods::IN_PLACE,
- QRMatrix<RMatrix>::colPivoting::FALSE
- );
- Rx.round();
-
- // Rx produced by TSQR is unique up to the sign, hence the revert
- for (scalar& x : Rx)
{
- x *= -1;
+ QRMatrix<RMatrix> QRM
+ (
+ Rx,
+ QRMatrix<RMatrix>::modes::ECONOMY,
+ QRMatrix<RMatrix>::outputs::ONLY_R
+ );
+ RMatrix& R = QRM.R();
+ Rx.transfer(R);
}
+ Rx.round();
Info<< tab << "Computing RxInv" << endl;
- RxInv_ = pinv(Rx);
+ RxInv_ = MatrixTools::pinv(Rx);
Info<< tab << "Computing A1" << endl;
- A1 = RxInv_*RMatrix(pinv(Rx*(G_^Rx)));
+ A1 = RxInv_*RMatrix(MatrixTools::pinv(Rx*(G_^Rx)));
Rx.clear();
}
Pstream::scatter(RxInv_);
@@ -323,10 +325,10 @@ reducedKoopmanOperator()
else
{
Info<< tab << "Computing RxInv" << endl;
- RxInv_ = pinv(Rx);
+ RxInv_ = MatrixTools::pinv(Rx);
Info<< tab << "Computing A1" << endl;
- A1 = RxInv_*RMatrix(pinv(Rx*(G_^Rx)));
+ A1 = RxInv_*RMatrix(MatrixTools::pinv(Rx*(G_^Rx)));
Rx.clear();
Info<< tab << "Computing A2" << endl;
@@ -486,7 +488,7 @@ void Foam::DMDModels::STDMD::amplitudes()
Info<< tab << "Computing amplitudes" << endl;
amps_.resize(R.m());
- const RCMatrix pEvecs(pinv(evecs_));
+ const RCMatrix pEvecs(MatrixTools::pinv(evecs_));
// amps_ = pEvecs*R;
for (label i = 0; i < amps_.size(); ++i)