From 8bfab7aa80b818b0188b6e4d7d039db0921578a9 Mon Sep 17 00:00:00 2001
From: Kutalmis Bercin <kutalmis.bercin@esi-group.com>
Date: Mon, 12 Apr 2021 11:59:45 +0100
Subject: [PATCH] ENH: QRMatrix: improve stability in pinv
---
src/OpenFOAM/matrices/QRMatrix/QRMatrix.C | 53 +++++++++++++++++------
src/OpenFOAM/matrices/QRMatrix/QRMatrix.H | 2 +-
2 files changed, 41 insertions(+), 14 deletions(-)
diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
index fa9f6833afa..0b63f2017b0 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-2020 OpenCFD Ltd.
+ Copyright (C) 2019-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -492,9 +492,10 @@ template<class MatrixType>
MatrixType Foam::pinv
(
const MatrixType& A,
- const scalar tol
+ const scalar tolerance
)
{
+ scalar tol = tolerance;
typedef typename MatrixType::cmptType cmptType;
if (A.empty())
@@ -506,7 +507,12 @@ MatrixType Foam::pinv
if (A.size() == 1)
{
- return MatrixType({1, 1}, cmptType(1.0)/(A(0,0) + cmptType(VSMALL)));
+ if (A(0,0) == cmptType(0))
+ {
+ return MatrixType({1, 1}, cmptType(0));
+ }
+
+ return MatrixType({1, 1}, cmptType(1)/(A(0,0) + cmptType(VSMALL)));
}
QRMatrix<MatrixType> QRM
@@ -521,27 +527,47 @@ MatrixType Foam::pinv
// R1 (KP:p. 648)
// Find the first diagonal element index with (almost) zero value in R
label firstZeroElemi = 0;
+ label i = 0;
+ while (i < 2)
{
const List<cmptType> diag(R.diag());
- auto lessT = [&](const cmptType& x) { return mag(x) < tol; };
+ auto lessThan = [=](const cmptType& x) { return tol > mag(x); };
firstZeroElemi =
std::distance
(
diag.cbegin(),
- std::find_if(diag.cbegin(), diag.cend(), lessT)
+ std::find_if(diag.cbegin(), diag.cend(), lessThan)
);
- }
- if (firstZeroElemi == 0)
- {
- WarningInFunction
- << "The largest (magnitude) diagonal element is (almost) zero."
- << nl << "Returning a zero matrix."
- << endl;
+ if (firstZeroElemi == 0)
+ {
+ if (i == 0)
+ {
+ WarningInFunction
+ << "The largest diagonal element magnitude is nearly zero. "
+ << "Tightening the tolerance."
+ << endl;
- return MatrixType(A.sizes(), Zero);
+ tol = 1e-13;
+ ++i;
+ }
+ else
+ {
+ WarningInFunction
+ << "The largest diagonal element magnitude is nearly zero. "
+ << "Returning a zero matrix."
+ << endl;
+ ++i;
+
+ return MatrixType({A.n(), A.m()}, Zero);
+ }
+ }
+ else
+ {
+ i += 2;
+ }
}
// Exclude the first (almost) zero diagonal row and the rows below
@@ -599,4 +625,5 @@ MatrixType Foam::pinv
}
}
+
// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H
index 415fbba853c..ee9bd62edf2 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 OpenCFD Ltd.
+ Copyright (C) 2019-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
--
GitLab