From 745624c024fb6b83bef6017b32bfa35972420667 Mon Sep 17 00:00:00 2001 From: kuti Date: Thu, 23 May 2019 11:32:45 +0100 Subject: [PATCH] ENH: partial overhaul of Matrix type (#1220) - additional operators: + compound assignment + inner product: operator& + outer product: operator^ - additional functions: - MatrixBlock methods: subColumn, subRow, subMatrix - L2 norms for matrix or column - trace, diag, round, transpose - MatrixBlock methods: col(), block() are deprecated since their access patterns with (size, offset) are unnatural/unwieldy. - verifications by test/Matrix/Test-Matrix --- applications/test/Matrix/Test-Matrix.C | 1030 ++++++++++++++--- src/OpenFOAM/matrices/Matrix/Matrix.C | 445 +++++-- src/OpenFOAM/matrices/Matrix/Matrix.H | 386 ++++-- src/OpenFOAM/matrices/Matrix/MatrixI.H | 254 ++-- src/OpenFOAM/matrices/Matrix/MatrixTools.C | 136 +++ src/OpenFOAM/matrices/Matrix/MatrixTools.H | 90 ++ .../matrices/MatrixBlock/MatrixBlock.C | 72 +- .../matrices/MatrixBlock/MatrixBlock.H | 32 +- .../matrices/MatrixBlock/MatrixBlockI.H | 99 +- src/OpenFOAM/matrices/QRMatrix/QRMatrix.C | 2 +- .../RectangularMatrix/RectangularMatrix.H | 6 +- .../RectangularMatrix/RectangularMatrixI.H | 4 +- .../matrices/SquareMatrix/SquareMatrix.C | 5 +- .../matrices/SquareMatrix/SquareMatrix.H | 3 +- .../matrices/SquareMatrix/SquareMatrixI.H | 4 +- 15 files changed, 1974 insertions(+), 594 deletions(-) create mode 100644 src/OpenFOAM/matrices/Matrix/MatrixTools.C create mode 100644 src/OpenFOAM/matrices/Matrix/MatrixTools.H diff --git a/applications/test/Matrix/Test-Matrix.C b/applications/test/Matrix/Test-Matrix.C index 4e001d7c6a..dbccf225fa 100644 --- a/applications/test/Matrix/Test-Matrix.C +++ b/applications/test/Matrix/Test-Matrix.C @@ -25,17 +25,23 @@ License \*---------------------------------------------------------------------------*/ -#include "scalarMatrices.H" +#include "MatrixTools.H" #include "LUscalarMatrix.H" #include "LLTMatrix.H" -#include "QRMatrix.H" -#include "vector.H" -#include "tensor.H" -#include "IFstream.H" - +#include "Random.H" +#include "SortList.H" #include using namespace Foam; +using namespace Foam::MatrixTools; + +#define isEqual MatrixTools::equal + +void horizontalLine() +{ + Info<< "+---------+---------+---------+---------+---------+" << nl; +} + // Copy values into matrix template @@ -59,254 +65,906 @@ void assignMatrix std::copy(list.begin(), list.end(), mat.begin()); } -// Create matrix with values + template -MatrixType makeMatrix +MatrixType makeRandomMatrix ( const labelPair& dims, - std::initializer_list list + Random& rndGen ) { MatrixType mat(dims); - assignMatrix(mat, list); + std::generate + ( + mat.begin(), + mat.end(), + [&]{return rndGen.GaussNormal();} + ); return mat; } -// Create matrix with values -template -MatrixType makeMatrix -( - std::initializer_list list -) +template +Ostream& operator<<(Ostream& os, const ConstMatrixBlock& mat) { - MatrixType mat(labelPair(nRows, nCols)); - - assignMatrix(mat, list); - - return mat; + return MatrixTools::printMatrix(os, mat); } -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -// Main program: +template +Ostream& operator<<(Ostream& os, const MatrixBlock& mat) +{ + return MatrixTools::printMatrix(os, mat); +} + + +// * * * * * * * * * * * * * * * Main Program * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { - // SquareMatrix + typedef SquareMatrix SMatrix; + typedef SquareMatrix SCMatrix; + typedef RectangularMatrix RMatrix; + typedef RectangularMatrix RCMatrix; + + Random rndGen(1234); + + #if 1 { - Info<< nl << "Test SquareMatrix" << nl; + horizontalLine(); - SquareMatrix square1(3); + Info<< "## Matrix Constructors & Assignments:" << nl << nl; + Info<< "# SquareMatrix examples:" << nl; + { + const label mRows = 2; + const label nCols = 2; + Matrix MS0(mRows, nCols); + Matrix MS1(mRows, nCols, Zero); + Matrix MS2(mRows, nCols, 17.44); + Matrix MS3(MS2); + + SMatrix S0(mRows); + SMatrix S1(mRows, Zero); + SMatrix S2(mRows, I); + SMatrix S3(mRows, 17.44); + SMatrix S4(2, Zero); + SMatrix S5(labelPair{mRows, nCols}); + SMatrix S6(mRows, nCols, Zero); + SMatrix S7({mRows, nCols}, Zero); + SMatrix S8({mRows, nCols}, 17.44); + assignMatrix + ( + S8, + { + 1, 2, // 0th row + 3, 4 // 1st row and so on + } + ); + + S1 = Zero; + S2 = 13.13; + S3 = I; + } + + Info<< "# RectangularMatrix examples:" << nl; + { + const label mRows = 2; + const label nCols = 3; + Matrix MR0(mRows, nCols); + Matrix MR1(mRows, nCols, Zero); + Matrix MR2(mRows, nCols, 17.44); + + RMatrix R0(mRows, nCols); + RMatrix R1(labelPair{mRows, nCols}); + RMatrix R2(mRows, nCols, Zero); + RMatrix R3({mRows, nCols}, Zero); + RMatrix R4(mRows, nCols, -17.44); + RMatrix R5({mRows, nCols}, -17.44); + RMatrix R6(3, 4, Zero); + RMatrix R7({3, 4}, Zero); + assignMatrix + ( + R7, + { + 1, 2, 3, 4, + 5, 6, 7, 8.8, + 9, 10, 11, 12 + } + ); + R0 = Zero; + R2 = -13.13; + } + + Info<< "# SquareMatrix examples:" << nl; + { + const label mRows = 2; + const label nCols = 2; + Matrix MS0(mRows, nCols); + Matrix MS1(mRows, nCols, Zero); + Matrix MS2(mRows, nCols, complex(17.44, 44.17)); + Matrix MS3(MS2); + + SCMatrix S0(mRows); + SCMatrix S1(mRows, Zero); + SCMatrix S2(mRows, I); + SCMatrix S3(mRows, complex(17.44,0)); + SCMatrix S4(2, Zero); + SCMatrix S5(labelPair{mRows, nCols}); + SCMatrix S6(mRows, nCols, Zero); + SCMatrix S7({mRows, nCols}, Zero); + SCMatrix S8({mRows, nCols}, complex(17.44,0)); + assignMatrix + ( + S8, + { + complex(1,0), complex(2,2), + complex(4,2), complex(5,0) + } + ); + + S1 = Zero; + S2 = complex(13.13,0); + S3 = I; + } + + Info<< nl; + horizontalLine(); + } + #endif + + + #if 1 + { + horizontalLine(); + + Info<< "## Access Functions:" << nl; + + SMatrix S(3, Zero); assignMatrix ( - square1, + S, { - -3.0, 10.0, -4.0, - 2.0, 3.0, 10.0, - 2.0, 6.0, 1.0 + -3, 11, -4, + 2, 3, 10, + 2, 6, 1 + } + ); + S(0,1) = 10; + + Info<< "# SquareMatrix example:" << nl; + printMatrix(Info, S) << nl; + + const label mRows = S.m(); + const label nCols = S.n(); + const label size = S.size(); + const labelPair sizes = S.sizes(); + const bool isEmpty = S.empty(); + const long constPointer = long(S.cdata()); + long pointer = long(S.data()); + + Info + << "Number of rows =" << tab << mRows << nl + << "Number of columns = " << tab << nCols << nl + << "Number of elements = " << tab << size << nl + << "Number of rows/columns = " << tab << sizes << nl + << "Matrix is empty = " << tab << isEmpty << nl + << "Constant pointer = " << tab << constPointer << nl + << "Pointer = " << tab << pointer << nl + << nl; + + horizontalLine(); + } + #endif + + + #if 1 + { + horizontalLine(); + + Info<< "## Block Access Functions:" << nl; + + RMatrix R(4, 3, Zero); + assignMatrix + ( + R, + { + -3, 11, -4, + 2, 3, 10, + 2, 6, 1, + 1, 2, 3 } ); - Info<< "matrix: " << square1 - << " begin: " << uintptr_t(square1.cdata()) << nl; + Info<< "# RectangularMatrix example:" << nl; + printMatrix(Info, R) << nl; - // Test makeMatrix + // Indices begin at 0 + const label columnI = 1; + const label rowI = 2; + const label colStartI = 1; + const label rowStartI = 1; + const label szRows = 2; + const label szCols = 2; + + Info + << "column: " << R.subColumn(columnI) << nl + << "sub-column: " << R.subColumn(columnI, rowStartI) << nl + << "sub-sub-column: " << R.subColumn(columnI, rowStartI, szRows) << nl + << "row: " << R.subRow(rowI) << nl + << "sub-row: " << R.subRow(rowI, colStartI) << nl + << "sub-sub-row: " << R.subRow(rowI, colStartI, szCols) << nl + << "sub-block: " << R.subMatrix(rowStartI, colStartI) << nl + << "sub-block with row size: " << R.subMatrix(rowStartI, colStartI, szRows) << nl + << "sub-block with row|col size: " << R.subMatrix(rowStartI, colStartI, szRows, szCols) << nl; + + horizontalLine(); + } + #endif + + + #if 1 + { + horizontalLine(); + + Info<< "## Edit Functions:" << nl; + + RMatrix R(2, 2, Zero); + assignMatrix + ( + R, + { + 1, 2, + 3, 4 + } + ); + + Info<< "# Before Matrix.clear():" << nl << R << nl; + R.clear(); + Info<< "# After Matrix.clear():" << nl << R << nl; + + Info<< "# Before Matrix.resize(m, n):" << nl << R; + R.resize(3, 2); + Info<< "# After Matrix.resize(m, n):" << nl << R << nl << nl; + + // Steal matrix contents + Info<< "# Before Matrix.release():" << nl << R << nl; + List newOwner(R.release()); + Info<< "# After Matrix.release():" << nl << R + << "& corresponding List content:" << nl << newOwner << nl << nl; + + R.resize(3, 2); + assignMatrix + ( + R, + { + 1, 2, + 5, 6, + 9e-19, 1e-17 + } + ); + + Info<< "# Before Matrix.round():" << nl << R << nl; + R.round(); + Info<< "# After Matrix.round():" << nl << R << nl << nl; + + Info<< "# Before Matrix.T() (Transpose):" << nl << R << nl; + RMatrix RT = R.T(); + Info<< "# After Matrix.T():" << nl << RT << nl << nl; + + RCMatrix RC(3, 2, Zero); + assignMatrix + ( + RC, + { + complex(1, 0), complex(2,5), + complex(4, 3), complex(1,-1), + complex(1, 2), complex(2,9) + } + ); + + Info<< "# Before Matrix.T() (Hermitian transpose):" << nl << RC << nl; + RCMatrix RCH = RC.T(); + Info<< "# After Matrix.T():" << nl << RCH << nl << nl; + + Info<< "# Diagonal and trace:" << nl; + RMatrix R1(5, 7, 3.4); { - auto square1b + List diagR = R1.diag(); + scalar traceR = R1.trace(); + Info<< "RectangularMatrix:" + << nl << R1 << nl + << "Major diagonal of RectangularMatrix:" + << nl << diagR << nl + << "Trace of RectangularMatrix:" + << nl << traceR << nl << nl; + } + + Info<< "# List assignment to the main diagonal of Matrix:" << nl; + { + SMatrix S1(5, I); + Info<< "Before List assignment:" << nl; + printMatrix(Info, S1); + List lst(5, 2.0); + S1.diag(lst); + Info<< "After List assignment:" << nl; + printMatrix(Info, S1); + } + + Info<< "# Diagonal sort of Matrix:" << nl; + { + auto S1 ( - makeMatrix + makeRandomMatrix ( - {3, 3}, - { - -3.0, 10.0, -4.0, - 2.0, 3.0, 10.0, - 2.0, 6.0, 1.0 - } + {3, 3}, rndGen ) ); - Info<< "makeMatrix: " << square1b << nl; + List diag = S1.diag(); + SortList incrDiag(diag); + incrDiag.sort(); - auto square1c - ( - makeMatrix - ({ - -3.0, 10.0, -4.0, - 2.0, 3.0, 10.0, - 2.0, 6.0, 1.0 - }) - ); - - Info<< "makeMatrix: " << square1c << nl; + Info<< "Diagonal list:" << diag << nl + << "Ascending diagonal list: " << incrDiag << nl << nl; } - //Info<< square1 - 2.0*(-square1) << nl; - Info<< "min:" << min(square1) << " max:" << max(square1) << nl; - Info<< "min/max: " << minMax(square1) << nl; - - // Steal matrix contents - - List stole(square1.release()); - - Info<< "matrix: " << square1 - << " begin: " << uintptr_t(square1.cdata()) << nl; - - Info<< "List: " << stole << nl; - - Info<< "min/max: " << minMax(square1) << nl; - - square1 = 100; - - Info<< "matrix: " << square1 - << " begin: " << uintptr_t(square1.cdata()) << nl; - - - SquareMatrix square2(3, I); - - square1 = square2; - - Info<< nl << "After copy assign from Identity:" << nl << square1 << nl; - - square1 += 1.25*square2; - - Info<< nl << "After +=" << nl << square1 << nl; - - square1 -= square2.T(); - - Info<< nl << "After -=" << nl << square1 << nl; - - square1 *= 10; - - Info<< nl << "After *=" << nl << square1 << nl; - - square1 /= 8; - - Info<< nl << "After /=" << nl << square1 << nl; - - SquareMatrix square4; - square4 = square2; - - Info<< nl << square4 << endl; - - SquareMatrix square5; - square4 = square5; - Info<< nl << square5 << endl; + horizontalLine(); } + #endif - // RectangularMatrix + + #if 1 { - Info<< nl << "Test RectangularMatrix" << nl; + horizontalLine(); - RectangularMatrix rm1({5, 6}, 3.1); - rm1(0, 1) = 4.5; + Info<< "## Transpose Verifications:" << nl; - RectangularMatrix rm1b(rm1.block(2, 2, 0, 0)); + const label numberOfTests = 10; - Info // << "Full matrix " << rm1 << nl - << "block = " << rm1b << endl; + for (label i = 0; i < numberOfTests; ++i) + { + const label mRows = rndGen.position(1,10); + const label nCols = rndGen.position(1,10); + + Info << "# SquareMatrix example:" << nl; + { + Info<< "A(mxm) where m = " << mRows << nl; + auto A(makeRandomMatrix({mRows, mRows}, rndGen)); + auto B(makeRandomMatrix({mRows, mRows}, rndGen)); + + Info<< nl << "# (A.T()).T() = A:" << nl; + isEqual((A.T()).T(), A, true); + + Info<< nl << "# (A + B).T() = A.T() + B.T():" << nl; + isEqual((A + B).T(), A.T() + B.T(), true); + + Info<< nl << "# (A*B).T() = B.T()*A.T():" << nl; + isEqual((A*B).T(), B.T()*A.T(), true); + Info<< nl << nl; + } + + Info << "# RectangularMatrix example:" << nl; + { + Info<< "A(mxn) where" + << " m = " << mRows << "," + << " n = " << nCols << nl; + auto A(makeRandomMatrix({mRows, nCols}, rndGen)); + auto B(makeRandomMatrix({mRows, nCols}, rndGen)); + auto C(makeRandomMatrix({nCols, mRows}, rndGen)); + + Info<< nl << "# (A.T()).T() = A:" << nl; + isEqual((A.T()).T(), A, true); + + Info<< nl << "# (A + B).T() = A.T() + B.T():" << nl; + isEqual((A + B).T(), A.T() + B.T(), true); + + Info<< nl << "# (A*C).T() = C.T()*C.T():" << nl; + isEqual((A*C).T(), C.T()*A.T(), true); + Info<< nl << nl; + } + + Info << "# SquareMatrix example:" << nl; + { + Info<< "A(mxm) where m = " << mRows << nl; + SCMatrix A(mRows); + SCMatrix B(mRows); + + for (auto& val : A) + { + val.Re() = rndGen.GaussNormal(); + val.Im() = rndGen.GaussNormal(); + } + for (auto& val : B) + { + val.Re() = rndGen.GaussNormal(); + val.Im() = rndGen.GaussNormal(); + } + + Info<< nl << "# (A.T()).T() = A:" << nl; + isEqual((A.T()).T(), A, true); + + Info<< nl << "# (A + B).T() = A.T() + B.T():" << nl; + isEqual((A + B).T(), A.T() + B.T(), true); + + Info<< nl << "# (A*B).T() = B.T()*A.T():" << nl; + isEqual((A*B).T(), B.T()*A.T(), true); + Info<< nl << nl; + } + } + + horizontalLine(); } + #endif + + #if 1 { - scalarSymmetricSquareMatrix symmMatrix(3, Zero); + horizontalLine(); - symmMatrix(0, 0) = 4; - symmMatrix(1, 0) = 12; - symmMatrix(1, 1) = 37; - symmMatrix(2, 0) = -16; - symmMatrix(2, 1) = -43; - symmMatrix(2, 2) = 98; + Info<< "## Scalar Arithmetic Operations:" << nl; - Info<< "Symmetric Square Matrix = " << symmMatrix << endl; + if (true) + { + SMatrix S1(3, Zero); + assignMatrix + ( + S1, + { + 4.1, 12.5, -16.3, + -192.3, -9.1, -3.0, + 1.0, 5.02, -4.4 + } + ); + SMatrix S2 = S1; + SMatrix S3(3, Zero); - Info<< "Inverse = " << inv(symmMatrix) << endl; - Info<< "Determinant = " << det(symmMatrix) << endl; + Info<< "SquareMatrix S1 = " << S1 << nl; + Info<< "SquareMatrix S2 = " << S2 << nl; + Info<< "SquareMatrix S3 = " << S3 << nl; - scalarSymmetricSquareMatrix symmMatrix2(symmMatrix); - LUDecompose(symmMatrix2); + S3 = S1*S2; + Info<< "S1*S2 = " << S3 << nl; - Info<< "Inverse = " << invDecomposed(symmMatrix2) << endl; - Info<< "Determinant = " << detDecomposed(symmMatrix2) << endl; + S3 = S1 - S2; + Info<< "S1 - S2 = " << S3 << nl; - scalarDiagonalMatrix rhs(3, Zero); + S3 = S1 + S2; + Info<< "S1 + S2 = " << S3 << nl; + +// S3 *= S1; // Not Implemented + + S3 -= S1; + Info<< "S3 -= S1; S3 = " << S3 << nl; + + S3 += S1; + Info<< "S3 += S1; S3 = " << S3 << nl; + + + Info<< nl << "# Scalar broadcasting:" << nl; + + S3 *= 5.0; + Info<< "S3 *= 5.0; S3 = " << S3 << nl; + + S3 /= 5.0; + Info<< "S3 /= 5.0; S3 = " << S3 << nl; + + S3 -= 1.0; + Info<< "S3 -= 1.0; S3 = " << S3 << nl; + + S3 += 1.0; + Info<< "S3 += 1.0; S3 = " << S3 << nl; + + + Info<< nl << "# Operations between different matrix types:" << nl; + RMatrix R1 = S1; + Info<< "RectangularMatrix R1 = " << R1 << nl; + + // RectangularMatrix*SquareMatrix returns RectangularMatrix + // S3 = S3*R1; // Not implemented + R1 = S3*R1; + Info<< "R1 = S3*R1; R1 = " << R1 << nl; + + S3 = S3 - R1; + Info<< "S3 = S3 - R1; S3 = " << S3 << nl; + + S3 = S3 + R1; + Info<< "S3 = S3 + R1; S3 = " << S3 << nl; + + R1 = S3 + R1; + Info<< "R1 = S3 + R1; R1 = " << R1 << nl; + + + Info<< nl << "# Extrema operations:" << nl; + + Info<< "min(S3) = " << min(S3) << nl + << "max(S3) = " << max(S3) << nl + << "minMax(S3) = " << minMax(S3) << nl; + } + + if (true) + { + Info<< nl << "# Operations with ListTypes:" << nl; + + SMatrix A(3, Zero); + assignMatrix + ( + A, + { + 4.1, 12.5, -16.3, + -192.3, -9.1, -3.0, + 1.0, 5.02, -4.4 + } + ); + const List lst({1, 2, 3}); + RMatrix colVector(3, 1, Zero); + assignMatrix + ( + colVector, + {1, 2, 3} + ); + RMatrix rowVector(1, 3, Zero); + assignMatrix + ( + rowVector, + {1, 2, 3} + ); + + Info<< "A = " << A << nl; + Info<< "colVector = " << colVector << nl; + Info<< "rowVector = " << rowVector << nl; + + const Field field1(A.Amul(lst)); + const Field field2(A*lst); + const Field field3(A.Tmul(lst)); + const Field field4(lst*A); + + Info + << "Field1 = A*lst = A.Amul(lst):" << nl << field1 << nl + << "Field2 = A*lst:" << nl << field2 << nl + << "A*colVector:" << nl << A*colVector << nl + << "Field3 = lst*A = A.Tmul(lst):" << nl << field3 << nl + << "Field4 = lst*A:" << nl << field4 << nl + << "rowVector*A:" << nl << rowVector*A << nl + << nl; + } + + if (true) + { + Info<< nl << "# Implicit inner/outer products:" << nl; + + const scalar r = 2.0; + RMatrix A(3, 2, Zero); + assignMatrix + ( + A, + { + 1, 2, + 3, 4, + 5, 6 + } + ); + const RMatrix B(A); + const RMatrix C(B); + + Info<< nl << "# Inner product:" << nl; + { + Info<< "- Explicit vs Implicit => A.T()*B == A&B:" << nl; + isEqual(A.T()*B, A&B, true); + + Info<< "- Commutative => A&B == B&A:" << nl; + isEqual(A&B, B&A, true); + + Info<< "- Distributive => A&(B+C) == A&B + A&C:" << nl; + isEqual(A&(B + C), (A&B) + (A&C), true); + + Info<< "- Bilinear => A&(rB+C) == r(A&B) + (A&C):" << nl; + isEqual(A&(r*B + C), r*(A&B) + (A&C), true); + + Info<< "- Scalar multiplication => (rA)&(rB) == rr(A&B):" << nl; + isEqual((r*A)&(r*B), r*r*(A&B), true); + } + + Info<< nl << "# Outer product:" << nl; + { + Info<< "- Explicit vs Implicit => A*B.T() == A^B:" << nl; + isEqual(A*B.T(), A^B, true); + + Info<< "- Commutative => A^B == B^A:" << nl; + isEqual(A^B, B^A, true); + + Info<< "- Distributive => A^(B+C) == A^B + A^C:" << nl; + isEqual(A^(B + C), (A^B) + (A^C), true); + + Info<< "- Scalar multiplication => r*(A^B) == (r*A)^B:" << nl; + isEqual(r*(A^B), (r*A)^B, true); + } + } + + Info<< nl; + horizontalLine(); + } + #endif + + + #if 1 + { + horizontalLine(); + + Info<< "## Complex Arithmetic Operations:" << nl; + + if (true) + { + SCMatrix S1(3, Zero); + assignMatrix + ( + S1, + { + complex(4.1, 1.0), complex(12.5, -1), complex(-16.3,-3), + complex(-192.3, 5), complex(-9.1, 3), complex(-3.0, 1), + complex(1.0, 3), complex(5.02, 0.3), complex(-4.4, 1) + } + ); + SCMatrix S2 = S1; + SCMatrix S3(3, Zero); + + Info<< "SquareMatrix S1 = " << S1 << nl; + Info<< "SquareMatrix S2 = " << S2 << nl; + Info<< "SquareMatrix S3 = " << S3 << nl; + + S3 = S1*S2; + Info<< "S1*S2 = " << S3 << nl; + + S3 = S1 - S2; + Info<< "S1 - S2 = " << S3 << nl; + + S3 = S1 + S2; + Info<< "S1 + S2 = " << S3 << nl; + +// S3 *= S1; // Not Implemented + + S3 -= S1; + Info<< "S3 -= S1; S3 = " << S3 << nl; + + S3 += S1; + Info<< "S3 += S1; S3 = " << S3 << nl; + + + Info<< nl << "# Complex broadcasting:" << nl; + + S3 *= complex(5, 0); + Info<< "S3 *= complex(5, 0); S3 = " << S3 << nl; + + S3 /= complex(5, 0); + Info<< "S3 /= complex(5, 0); S3 = " << S3 << nl; + + S3 -= complex(1, 0); + Info<< "S3 -= complex(1, 0); S3 = " << S3 << nl; + + S3 += complex(1, 0); + Info<< "S3 += complex(1, 0); S3 = " << S3 << nl; + + + Info<< nl << "# Operations between different matrix types:" << nl; + RCMatrix R1 = S1; + Info<< "RectangularMatrix R1 = " << R1 << nl; + + R1 = S3*R1; + Info<< "R1 = S3*R1; R1 = " << R1 << nl; + + S3 = S3 - R1; + Info<< "S3 = S3 - R1; S3 = " << S3 << nl; + + S3 = S3 + R1; + Info<< "S3 = S3 + R1:" << S3 << nl; + + // Extrama operations // Not Implemented + } + + if (true) + { + Info<< nl << "# Operations with ListTypes:" << nl; + + SCMatrix A(3, Zero); + assignMatrix + ( + A, + { + complex(4.1, 1.0), complex(12.5, -1), complex(-16.3,-3), + complex(-192.3, 5), complex(-9.1, 3), complex(-3.0, 1), + complex(1.0, 3), complex(5.02, 0.3), complex(-4.4, 1) + } + ); + const List lst({complex(1,1), complex(2,2), complex(3,3)}); + RCMatrix colVector(3, 1, Zero); + assignMatrix + ( + colVector, + {complex(1,1), complex(2,2), complex(3,3)} + ); + RCMatrix rowVector(1, 3, Zero); + assignMatrix + ( + rowVector, + {complex(1,1), complex(2,2), complex(3,3)} + ); + + Info<< "A = " << A << nl; + Info<< "colVector = " << colVector << nl; + Info<< "rowVector = " << rowVector << nl; + + const Field field1(A.Amul(lst)); + const Field field2(A*lst); + const Field field3(A.Tmul(lst)); + const Field field4(lst*A); + + Info + << "Field1 = A*lst = A.Amul(lst):" << nl << field1 << nl + << "Field2 = A*lst:" << nl << field2 << nl + << "A*colVector:" << nl << A*colVector << nl + << "Field3 = lst*A = A.Tmul(lst):" << nl << field3 << nl + << "Field4 = lst*A:" << nl << field4 << nl + << "rowVector*A:" << nl << rowVector*A << nl + << nl; + } + + #if 1 + { + Info<< nl << "# Implicit inner/outer products:" << nl; + + const complex r(2.0,1.0); + RCMatrix A(3, 2, Zero); + assignMatrix + ( + A, + { + complex(1,0), complex(2,1), + complex(3,-0.3), complex(4,-1), + complex(0.5,-0.1), complex(6,0.4) + } + ); + const RCMatrix B(A); + const RCMatrix C(B); + + Info<< nl << "# Inner product:" << nl; + { + Info<< "- Explicit vs Implicit => A.T()*B == A&B:" << nl; + isEqual(A.T()*B, A&B, true); + + Info<< "- Commutative => A&B == B&A:" << nl; + isEqual(A&B, B&A, true); + + Info<< "- Distributive => A&(B+C) == A&B + A&C:" << nl; + isEqual(A&(B + C), (A&B) + (A&C), true); + + Info<< "- Bilinear => A&(rB+C) == r(A&B) + (A&C):" << nl; + isEqual(A&(r*B + C), r*(A&B) + (A&C), true); + } + + Info<< nl << "# Outer product:" << nl; + { + Info<< "- Explicit vs Implicit => A*B.T() == A^B:" << nl; + isEqual(A*B.T(), A^B, true); + + Info<< "- Commutative => A^B == B^A:" << nl; + isEqual(A^B, B^A, true); + + Info<< "- Distributive => A^(B+C) == A^B + A^C:" << nl; + isEqual(A^(B + C), (A^B) + (A^C), true); + + Info<< "- Scalar multiplication => r*(A^B) == (r*A)^B:" << nl; + isEqual(r*(A^B), (r*A)^B, true); + } + } + #endif + + Info<< nl; + horizontalLine(); + } + #endif + + + #if 1 + { + horizontalLine(); + + Info<< "## SymmetricSquareMatrix algorithms:" << nl; + + scalarSymmetricSquareMatrix symm(3, Zero); + + symm(0, 0) = 4; + symm(1, 0) = 12; + symm(1, 1) = 37; + symm(2, 0) = -16; + symm(2, 1) = -43; + symm(2, 2) = 98; + + Info<< "SymmetricSquareMatrix = " << nl << symm << nl + << "Inverse = " << nl << inv(symm) << nl + << "Determinant = " << det(symm) << nl; + + scalarSymmetricSquareMatrix symm2(symm); + LUDecompose(symm2); + Info<< "LU decomposition:" << nl + << "Inverse = " << nl << invDecomposed(symm2) << nl + << "Determinant = " << detDecomposed(symm2) << nl; + + scalarDiagonalMatrix rhs(3, 0); rhs[0] = 1; rhs[1] = 2; rhs[2] = 3; - LUsolve(symmMatrix, rhs); - - Info<< "Decomposition = " << symmMatrix << endl; - Info<< "Solution = " << rhs << endl; + LUsolve(symm, rhs); + Info<< "Solving linear system through LU decomposition:" << nl + << "Decomposition = " << nl << symm << nl + << "Solution = " << rhs << nl; } + #endif - scalarSquareMatrix squareMatrix(3, Zero); - - squareMatrix(0, 0) = 4; - squareMatrix(0, 1) = 12; - squareMatrix(0, 2) = -16; - squareMatrix(1, 0) = 12; - squareMatrix(1, 1) = 37; - squareMatrix(1, 2) = -43; - squareMatrix(2, 0) = -16; - squareMatrix(2, 1) = -43; - squareMatrix(2, 2) = 98; - - Info<< nl << "Square Matrix = " << squareMatrix << endl; - - const scalarField source(3, 1); - + #if 1 { + scalarSquareMatrix squareMatrix(3, Zero); + + scalarSquareMatrix S(3, Zero); + assignMatrix + ( + S, + { + 4, 12, -16, + 12, 37, -43, + -16, -43, 98 + } + ); + + const scalarField source(3, 1); + + Info<< nl << "Square Matrix = " << S << endl; + + if (true) { - scalarSquareMatrix sm(squareMatrix); - Info<< "det = " << det(sm) << endl; + { + scalarSquareMatrix S2(S); + Info<< "Determinant = " << det(S2) << nl; + } + + scalarSquareMatrix S2(S); + labelList rhs(3, 0); + label sign; + LUDecompose(S2, rhs, sign); + + Info<< "LU decomposition = " << S2 << nl + << "Pivots = " << rhs << nl + << "Sign = " << sign << nl + << "Determinant = " << detDecomposed(S2, sign) << nl; } - scalarSquareMatrix sm(squareMatrix); - labelList rhs(3, 0); - label sign; - LUDecompose(sm, rhs, sign); + if (true) + { + LUscalarMatrix LU(S); + scalarField x(LU.solve(source)); + scalarSquareMatrix inv(3); + LU.inv(inv); + Info<< "LU.solve(source) residual: " << (S*x - source) + << nl << "LU inv " << inv << nl + << "LU inv*S " << (inv*S) << nl; + } - Info<< "Decomposition = " << sm << endl; - Info<< "Pivots = " << rhs << endl; - Info<< "Sign = " << sign << endl; - Info<< "det = " << detDecomposed(sm, sign) << endl; + if (true) + { + LLTMatrix LLT(S); + scalarField x(LLT.solve(source)); + Info<< "LLT solve residual " << (S*x - source) << nl; + } + + horizontalLine(); } + #endif - { - LUscalarMatrix LU(squareMatrix); - scalarField x(LU.solve(source)); - Info<< "LU solve residual " << (squareMatrix*x - source) << endl; - scalarSquareMatrix inv(3); - LU.inv(inv); - Info<< "LU inv " << inv << endl; - Info<< "LU inv*squareMatrix " << (inv*squareMatrix) << endl; - } - - { - LLTMatrix LLT(squareMatrix); - scalarField x(LLT.solve(source)); - Info<< "LLT solve residual " << (squareMatrix*x - source) << endl; - } - - { - QRMatrix QR(squareMatrix); - scalarField x(QR.solve(source)); - - Info<< "QR solve residual " - << (squareMatrix*x - source) << endl; - - Info<< "QR inverse solve residual " - << (x - QR.inv()*source) << endl; - - Info<< "QR inv *squareMatrix " << (QR.inv()*squareMatrix) << endl; - } - - Info<< "\nEnd\n" << endl; + Info<< nl << "End" << nl; return 0; } diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.C b/src/OpenFOAM/matrices/Matrix/Matrix.C index 1b0908dc3c..072c78ffe7 100644 --- a/src/OpenFOAM/matrices/Matrix/Matrix.C +++ b/src/OpenFOAM/matrices/Matrix/Matrix.C @@ -334,16 +334,29 @@ void Foam::Matrix::resize(const label m, const label n) } +template +void Foam::Matrix::round(const scalar tol) +{ + for (Type& val : *this) + { + if (mag(val) < tol) + { + val = Zero; + } + } +} + + template Form Foam::Matrix::T() const { - Form At(n(), m()); + 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); + At(j, i) = Detail::conj((*this)(i, j)); } } @@ -351,6 +364,91 @@ Form Foam::Matrix::T() const } +template +Foam::List Foam::Matrix::diag() const +{ + const label len = Foam::min(mRows_, nCols_); + + List result(len); + + for (label i=0; i < len; ++i) + { + result[i] = (*this)(i, i); + } + + return result; +} + + +template +void Foam::Matrix::diag(const UList& list) +{ + const label len = Foam::min(mRows_, nCols_); + + #ifdef FULLDEBUG + if (list.size() != len) + { + FatalErrorInFunction + << "List size (" << list.size() + << ") incompatible with Matrix diagonal" << abort(FatalError); + } + #endif + + for (label i=0; i < len; ++i) + { + (*this)(i, i) = list[i]; + } +} + + +template +Type Foam::Matrix::trace() const +{ + const label len = Foam::min(mRows_, nCols_); + + Type val = Zero; + + for (label i=0; i < len; ++i) + { + val += (*this)(i, i); + } + + return val; +} + + +template +Foam::scalar Foam::Matrix::columnNorm +( + const label colIndex, + const bool noSqrt +) const +{ + scalar result = Zero; + + for (label i=0; i < mRows_; ++i) + { + result += magSqr((*this)(i, colIndex)); + } + + return noSqrt ? result : Foam::sqrt(result); +} + + +template +Foam::scalar Foam::Matrix::norm(const bool noSqrt) const +{ + scalar result = Zero; + + for (const Type& val : *this) + { + result += magSqr(val); + } + + return noSqrt ? result : Foam::sqrt(result); +} + + // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // template @@ -416,9 +514,9 @@ void Foam::Matrix::operator= const MatrixBlock& Mb ) { - for (label i=0; i < mRows_; ++i) + for (label i = 0; i < mRows_; ++i) { - for (label j=0; j < nCols_; ++j) + for (label j = 0; j < nCols_; ++j) { (*this)(i, j) = Mb(i, j); } @@ -443,6 +541,7 @@ void Foam::Matrix::operator=(const zero) template void Foam::Matrix::operator+=(const Matrix& other) { + #ifdef FULLDEBUG if (this == &other) { FatalErrorInFunction @@ -458,15 +557,13 @@ void Foam::Matrix::operator+=(const Matrix& other) << other.m() << ", " << other.n() << ')' << nl << abort(FatalError); } + #endif - Type* out = this->data(); - const Type* in = other.cdata(); - - const label len = this->size(); - - for (label idx = 0; idx < len; ++idx) + auto iter2 = other.cbegin(); + for (Type& val : *this) { - out[idx] += in[idx]; + val += *iter2; + ++iter2; } } @@ -474,6 +571,7 @@ void Foam::Matrix::operator+=(const Matrix& other) template void Foam::Matrix::operator-=(const Matrix& other) { + #ifdef FULLDEBUG if (this == &other) { FatalErrorInFunction @@ -489,21 +587,39 @@ void Foam::Matrix::operator-=(const Matrix& other) << other.m() << ", " << other.n() << ')' << nl << abort(FatalError); } + #endif - Type* out = this->data(); - const Type* in = other.cdata(); - - const label len = this->size(); - - for (label idx=0; idx < len; ++idx) + auto iter2 = other.cbegin(); + for (Type& val : *this) { - out[idx] -= in[idx]; + val -= *iter2; + ++iter2; } } template -void Foam::Matrix::operator*=(const scalar s) +void Foam::Matrix::operator+=(const Type& s) +{ + for (Type& val : *this) + { + val += s; + } +} + + +template +void Foam::Matrix::operator-=(const Type& s) +{ + for (Type& val : *this) + { + val -= s; + } +} + + +template +void Foam::Matrix::operator*=(const Type& s) { for (Type& val : *this) { @@ -513,7 +629,7 @@ void Foam::Matrix::operator*=(const scalar s) template -void Foam::Matrix::operator/=(const scalar s) +void Foam::Matrix::operator/=(const Type& s) { for (Type& val : *this) { @@ -522,7 +638,7 @@ void Foam::Matrix::operator/=(const scalar s) } -// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * // +// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * * // template const Type& Foam::max(const Matrix& mat) @@ -569,7 +685,7 @@ Foam::MinMax Foam::minMax(const Matrix& mat) template Form Foam::operator-(const Matrix& mat) { - Form result(mat.m(), mat.n()); + Form result(mat.sizes()); std::transform ( @@ -583,9 +699,14 @@ Form Foam::operator-(const Matrix& mat) } -template -Form Foam::operator+(const Matrix& A, const Matrix& B) +template +Form1 Foam::operator+ +( + const Matrix& A, + const Matrix& B +) { + #ifdef FULLDEBUG if (A.m() != B.m() || A.n() != B.n()) { FatalErrorInFunction @@ -594,27 +715,31 @@ Form Foam::operator+(const Matrix& A, const Matrix& B) << B.m() << ", " << B.n() << ')' << nl << abort(FatalError); } + #endif - Form AB(A.m(), A.n()); + Form1 result(A.sizes()); - Type* ABv = AB.data(); - const Type* Av = A.cdata(); - const Type* Bv = B.cdata(); + std::transform + ( + A.cbegin(), + A.cend(), + B.cbegin(), + result.begin(), + std::plus() + ); - const label len = A.size(); - - for (label idx = 0; idx < len; ++idx) - { - ABv[idx] = Av[idx] + Bv[idx]; - } - - return AB; + return result; } -template -Form Foam::operator-(const Matrix& A, const Matrix& B) +template +Form1 Foam::operator- +( + const Matrix& A, + const Matrix& B +) { + #ifdef FULLDEBUG if (A.m() != B.m() || A.n() != B.n()) { FatalErrorInFunction @@ -623,85 +748,117 @@ Form Foam::operator-(const Matrix& A, const Matrix& B) << B.m() << ", " << B.n() << ')' << nl << abort(FatalError); } + #endif - Form AB(A.m(), A.n()); + Form1 result(A.sizes()); - Type* ABv = AB.data(); - const Type* Av = A.cdata(); - const Type* Bv = B.cdata(); - - const label len = A.size(); - - for (label idx=0; idx < len; ++idx) - { - ABv[idx] = Av[idx] - Bv[idx]; - } - - return AB; -} - - -template -Form Foam::operator*(const scalar s, const Matrix& mat) -{ - Form result(mat.m(), mat.n()); - - const label len = mat.size(); - - if (len) - { - Type* out = result.data(); - const Type* in = mat.cdata(); - - for (label idx = 0; idx < len; ++idx) - { - out[idx] = s * in[idx]; - } - } + std::transform + ( + A.cbegin(), + A.cend(), + B.cbegin(), + result.begin(), + std::minus() + ); return result; } template -Form Foam::operator*(const Matrix& mat, const scalar s) +Form Foam::operator*(const Type& s, const Matrix& mat) { - Form result(mat.m(), mat.n()); + Form result(mat.sizes()); - const label len = mat.size(); - - if (len) - { - Type* out = result.data(); - const Type* in = mat.cdata(); - - for (label idx=0; idx < len; ++idx) - { - out[idx] = in[idx] * s; - } - } + std::transform + ( + mat.cbegin(), + mat.cend(), + result.begin(), + [&](const Type& val) { return s * val; } + ); return result; } template -Form Foam::operator/(const Matrix& mat, const scalar s) +Form Foam::operator+(const Type& s, const Matrix& mat) { - Form result(mat.m(), mat.n()); + Form result(mat.sizes()); - const label len = mat.size(); + std::transform + ( + mat.cbegin(), + mat.cend(), + result.begin(), + [&](const Type& val) { return s + val; } + ); - if (len) - { - Type* out = result.data(); - const Type* in = mat.cdata(); + return result; +} - for (label idx=0; idx < len; ++idx) - { - out[idx] = in[idx] / s; - } - } + +template +Form Foam::operator-(const Type& s, const Matrix& mat) +{ + Form result(mat.sizes()); + + std::transform + ( + mat.cbegin(), + mat.end(), + result.begin(), + [&](const Type& val) { return s - val; } + ); + + return result; +} + + +template +Form Foam::operator*(const Matrix& mat, const Type& s) +{ + return s*mat; +} + + +template +Form Foam::operator+(const Matrix& mat, const Type& s) +{ + return s + mat; +} + + +template +Form Foam::operator-(const Matrix& mat, const Type& s) +{ + Form result(mat.sizes()); + + std::transform + ( + mat.cbegin(), + mat.end(), + result.begin(), + [&](const Type& val) { return val - s; } + ); + + return result; +} + + +template +Form Foam::operator/(const Matrix& mat, const Type& s) +{ + Form result(mat.sizes()); + + std::transform + ( + mat.cbegin(), + mat.end(), + result.begin(), + [&](const Type& val) { return val / s; } + ); return result; } @@ -715,6 +872,7 @@ Foam::operator* const Matrix& B ) { + #ifdef FULLDEBUG if (A.n() != B.m()) { FatalErrorInFunction @@ -724,6 +882,7 @@ Foam::operator* << "The columns of A must equal rows of B" << abort(FatalError); } + #endif typename typeOfInnerProduct::type AB ( @@ -732,13 +891,97 @@ Foam::operator* Zero ); - for (label i=0; i < AB.m(); ++i) + for (label i = 0; i < AB.m(); ++i) { - for (label j=0; j < AB.n(); ++j) + for (label k = 0; k < B.m(); ++k) { - for (label k=0; k < B.m(); ++k) + for (label j = 0; j < AB.n(); ++j) { - AB(i, j) += A(i, k) * B(k, j); + AB(i, j) += A(i, k)*B(k, j); + } + } + } + + return AB; +} + + +template +typename Foam::typeOfInnerProduct::type +Foam::operator& +( + const Matrix& AT, + const Matrix& B +) +{ + #ifdef FULLDEBUG + if (AT.m() != B.m()) + { + FatalErrorInFunction + << "Attempt to multiply incompatible matrices:" << nl + << "Matrix A : (" << AT.m() << ", " << AT.n() << ')' << nl + << "Matrix B : (" << B.m() << ", " << B.n() << ')' << nl + << "The rows of A must equal rows of B" + << abort(FatalError); + } + #endif + + typename typeOfInnerProduct::type AB + ( + AT.n(), + B.n(), + Zero + ); + + for (label i = 0; i < AB.m(); ++i) + { + for (label k = 0; k < B.m(); ++k) + { + for (label j = 0; j < AB.n(); ++j) + { + AB(i, j) += Detail::conj(AT(k, i))*B(k, j); + } + } + } + + return AB; +} + + +template +typename Foam::typeOfInnerProduct::type +Foam::operator^ +( + const Matrix& A, + const Matrix& BT +) +{ + #ifdef FULLDEBUG + if (A.n() != BT.n()) + { + FatalErrorInFunction + << "Attempt to multiply incompatible matrices:" << nl + << "Matrix A : (" << A.m() << ", " << A.n() << ')' << nl + << "Matrix B : (" << BT.m() << ", " << BT.n() << ')' << nl + << "The columns of A must equal columns of B" + << abort(FatalError); + } + #endif + + typename typeOfInnerProduct::type AB + ( + A.m(), + BT.m(), + Zero + ); + + for (label i = 0; i < AB.m(); ++i) + { + for (label k = 0; k < BT.n(); ++k) + { + for (label j = 0; j < AB.n(); ++j) + { + 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 bfe1c4e51f..2aef390896 100644 --- a/src/OpenFOAM/matrices/Matrix/Matrix.H +++ b/src/OpenFOAM/matrices/Matrix/Matrix.H @@ -53,6 +53,7 @@ SourceFiles #include "Pair.H" #include "Field.H" #include "autoPtr.H" +#include "complex.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -74,7 +75,7 @@ class Matrix { // Private Data - //- Number of rows and columns in Matrix. + //- Number of rows and columns in Matrix label mRows_, nCols_; //- Vector of values of type Type @@ -86,11 +87,11 @@ class Matrix //- Allocate storage for the contents inline void doAlloc(); - //- Multiply matrix with vector (A * x) + //- Right-multiply Matrix by a column vector (A * x) template tmp> AmulImpl(const ListType& x) const; - //- Multiply matrix transpose with vector (AT * x, or x * A) + //- Left-multiply Matrix by a row vector (x * A) template tmp> TmulImpl(const ListType& x) const; @@ -138,7 +139,7 @@ public: //- Construct given number of rows/columns inline explicit Matrix(const labelPair& dims); - //- Construct given number of rows/columns. + //- Construct given number of rows/columns //- initializing all elements to zero inline Matrix(const labelPair& dims, const zero); @@ -152,15 +153,15 @@ public: //- Move construct Matrix(Matrix&& mat); - //- Copy constructor from matrix of a different form + //- Copy constructor from Matrix of a different form template explicit Matrix(const Matrix& mat); - //- Construct from a block of another matrix + //- Construct from a block of another Matrix template Matrix(const ConstMatrixBlock& Mb); - //- Construct from a block of another matrix + //- Construct from a block of another Matrix template Matrix(const MatrixBlock& Mb); @@ -179,34 +180,34 @@ public: // Access - //- Return the number of rows. + //- Return the number of rows inline label m() const noexcept; - //- Return the number of columns. + //- Return the number of columns inline label n() const noexcept; - //- Return the number of elements in matrix (m*n) + //- Return the number of elements in Matrix (m*n) inline label size() const; //- Return row/column sizes inline labelPair sizes() const; - //- Return true if the matrix is empty (ie, size() is zero) + //- Return true if Matrix is empty (i.e., size() is zero) inline bool empty() const noexcept; //- Return const pointer to the first data element, which can also - //- be used to address into the matrix contents + //- be used to address into Matrix contents inline const Type* cdata() const noexcept; //- Return pointer to the first data element, which can also - //- be used to address into the matrix contents + //- be used to address into Matrix contents inline Type* data() noexcept; - //- Return const pointer to data in the specified row. + //- Return const pointer to data in the specified row // Subscript checking only with FULLDEBUG inline const Type* rowData(const label irow) const; - //- Return pointer to data in the specified row. + //- Return pointer to data in the specified row // Subscript checking only with FULLDEBUG inline Type* rowData(const label irow); @@ -219,87 +220,108 @@ public: inline Type& at(const label idx); + // Block Access (const) - // Block Access + //- Return const column or column's subset of Matrix + // Return entire column by its index: M.subColumn(a); + // Return subset of a column starting from rowIndex: M.subColumn(a, b); + // Return subset of a column starting from rowIndex with szRows elems: + // M.subColumn(a, b, c); + inline ConstMatrixBlock subColumn + ( + const label colIndex, + const label rowIndex = 0, + label len = -1 + ) const; - inline ConstMatrixBlock block - ( - const label m, - const label n, - const label mStart, - const label nStart - ) const; + //- Return const row or const row's subset of Matrix + // Return entire row by its index: M.subRow(a); + // Return subset of a row starting from columnIndex: M.subRow(a,b); + // Return subset of a row starting from columnIndex with szCols elems: + // M.subRow(a, b, c); + inline ConstMatrixBlock subRow + ( + const label rowIndex, + const label colIndex = 0, + label len = -1 + ) const; - template - inline ConstMatrixBlock block - ( - const label mStart, - const label nStart - ) const; + //- Return const sub-block of Matrix + // Sub-block starting at columnIndex & rowIndex indices + inline ConstMatrixBlock subMatrix + ( + const label rowIndex, + const label colIndex, + label szRows = -1, + label szCols = -1 + ) const; - inline ConstMatrixBlock col - ( - const label m, - const label rowStart - ) const; - - inline ConstMatrixBlock col - ( - const label m, - const label mStart, - const label nStart - ) const; + //- Access Field as a ConstMatrixBlock + template + inline ConstMatrixBlock block + ( + const label rowIndex, + const label colIndex + ) const; - inline MatrixBlock block - ( - const label m, - const label n, - const label mStart, - const label nStart - ); + // Block Access (non-const) - template - inline MatrixBlock block - ( - const label mStart, - const label nStart - ); + //- Return column or column's subset of Matrix + inline MatrixBlock subColumn + ( + const label colIndex, + const label rowIndex = 0, + label len = -1 + ); - inline MatrixBlock col - ( - const label m, - const label rowStart - ); + //- Return row or row's subset of Matrix + inline MatrixBlock subRow + ( + const label rowIndex, + const label colIndex = 0, + label len = -1 + ); + + //- Return sub-block of Matrix + inline MatrixBlock subMatrix + ( + const label rowIndex, + const label colIndex, + label szRows = -1, + label szCols = -1 + ); + + //- Access Field as a MatrixBlock + template + inline MatrixBlock block + ( + const label rowIndex, + const label colIndex + ); - inline MatrixBlock col - ( - const label m, - const label mStart, - const label nStart - ); // Check - //- Check index i is within valid range (0 ... m-1). + //- Check index i is within valid range [0, m) inline void checki(const label irow) const; - //- Check index j is within valid range (0 ... n-1). + //- Check index j is within valid range [0, n) inline void checkj(const label jcol) const; //- Check that dimensions are positive, non-zero inline void checkSize() const; - //- True if all entries have identical values, and matrix is non-empty + //- True if all entries have identical values, and Matrix is non-empty inline bool uniform() const; // Edit - //- Clear the Matrix, i.e. set sizes to zero. + //- Clear Matrix, i.e. set sizes to zero void clear(); - //- Release storage management of the Matrix contents by transferring + //- Release storage management of Matrix contents by transferring //- management to a List List release(); @@ -307,44 +329,65 @@ public: void swap(Matrix& mat); //- Transfer the contents of the argument Matrix into this Matrix - //- and annul the argument Matrix. + //- and annul the argument Matrix void transfer(Matrix& mat); - //- Change the matrix dimensions, preserving the elements + //- Change Matrix dimensions, preserving the elements void resize(const label m, const label n); - //- Change the matrix dimensions, preserving the elements + //- Change Matrix dimensions, preserving the elements inline void setSize(const label m, const label n); - //- Resize the matrix without reallocating storage (unsafe) + //- Resize Matrix without reallocating storage (unsafe) inline void shallowResize(const label m, const label n); + //- Round to zero elements with magnitude smaller than tol (SMALL) + void round(const scalar tol = SMALL); + // Operations - //- Return the transpose of the matrix + //- Return (conjugate) transpose of Matrix Form T() const; - //- Multiply matrix with vector (A * x) + //- Right-multiply Matrix by a column vector (A * x) inline tmp> Amul(const UList& x) const; - //- Multiply matrix with vector (A * x) + //- Right-multiply Matrix by a column vector (A * x) template inline tmp> Amul ( const IndirectListBase& x ) const; - //- Multiply matrix transpose with vector (AT * x, or x * A) + //- Left-multiply Matrix by a row vector (x * A) inline tmp> Tmul(const UList& x) const; - //- Multiply matrix transpose with vector (AT * x, or x * A) + //- Left-multiply Matrix by a row vector (x * A) template inline tmp> Tmul ( const IndirectListBase& x ) const; + //- Extract the diagonal elements. Method may change in the future. + List diag() const; + + //- Assign diagonal of Matrix + void diag(const UList& list); + + //- Return the trace + Type trace() const; + + //- Return L2-Norm of chosen column + // Optional without sqrt for parallel usage. + scalar columnNorm(const label colIndex, const bool noSqrt=false) const; + + //- Return Frobenius norm of Matrix + // Optional without sqrt for parallel usage. + // https://en.wikipedia.org/wiki/Matrix_norm#Frobenius_norm + scalar norm(const bool noSqrt=false) const; + // Member Operators @@ -370,11 +413,11 @@ public: //- Move assignment void operator=(Matrix&& mat); - //- Assignment to a block of another matrix + //- Assignment to a block of another Matrix template void operator=(const ConstMatrixBlock& Mb); - //- Assignment to a block of another matrix + //- Assignment to a block of another Matrix template void operator=(const MatrixBlock& Mb); @@ -390,14 +433,20 @@ public: //- Matrix subtraction void operator-=(const Matrix& other); + //- Matrix scalar addition + void operator+=(const Type& s); + + //- Matrix scalar subtraction + void operator-=(const Type& s); + //- Matrix scalar multiplication - void operator*=(const scalar s); + void operator*=(const Type& s); //- Matrix scalar division - void operator/=(const scalar s); + void operator/=(const Type& s); - // Random access iterator (non-const) + // Random Access Iterator (non-const) //- Return an iterator to begin traversing a Matrix inline iterator begin(); @@ -406,7 +455,7 @@ public: inline iterator end(); - // Random access iterator (const) + // Random Access Iterator (const) //- Return const_iterator to begin traversing a constant Matrix inline const_iterator cbegin() const; @@ -447,10 +496,73 @@ public: { return v_; } + + //- Deprecated(2019-04) - use subMatrix() + // \deprecated(2019-04) - use subMatrix() + ConstMatrixBlock + FOAM_DEPRECATED_FOR(2019-04, "subMatrix() method") block + ( + const label m, + const label n, + const label mStart, + const label nStart + ) const + { + return this->subMatrix(mStart, nStart, m, n); + } + + //- Deprecated(2019-04) - use subMatrix() + // \deprecated(2019-04) - use subMatrix() + MatrixBlock + FOAM_DEPRECATED_FOR(2019-04, "subMatrix() method") block + ( + const label m, + const label n, + const label mStart, + const label nStart + ) + { + return this->subMatrix(mStart, nStart, m, n); + } + + + //- Deprecated(2019-04) - use subColumn() + // \deprecated(2019-04) - use subColumn() + ConstMatrixBlock + FOAM_DEPRECATED_FOR(2019-04, "subColumn() method") col + ( + const label m, + const label mStart, + const label nStart + ) const + { + return this->subColumn(nStart, mStart, m); + } + + //- Deprecated(2019-04) - use subColumn() + // \deprecated(2019-04) - use subColumn() + MatrixBlock + FOAM_DEPRECATED_FOR(2019-04, "subColumn() method") col + ( + const label m, + const label mStart, + const label nStart + ) + { + return this->subColumn(nStart, mStart, m); + } + + //- Deleted(2019-04) - use subColumn() + // \deprecated(2019-04) - use subColumn() + void col(const label m, const label rowStart) const = delete; + + //- Deleted(2019-04) - use subColumn() + // \deprecated(2019-04) - use subColumn() + void col(const label m, const label rowStart) = delete; }; -// IOstream Operators +// * * * * * * * * * * * * * * * IOstream Operator * * * * * * * * * * * * * // //- Read Matrix from Istream, discarding contents of existing Matrix. template @@ -462,70 +574,100 @@ template Ostream& operator<<(Ostream& os, const Matrix& mat); -// Global Functions, Operators +// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * * // -//- Find max value in the matrix +//- Find max value in Matrix template const Type& max(const Matrix& mat); -//- Find min value in the matrix +//- Find min value in Matrix template const Type& min(const Matrix& mat); -//- Find the min/max values of the matrix +//- Find the min/max values of Matrix template MinMax minMax(const Matrix& mat); + +// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * * // + //- Matrix negation template Form operator-(const Matrix& mat); -//- Matrix addition -template -Form operator+ +//- Matrix addition. Returns Matrix of the same form as the first parameter. +template +Form1 operator+ ( - const Matrix& A, - const Matrix& B + const Matrix& A, + const Matrix& B ); - -//- Matrix subtraction -template -Form operator- +//- Matrix subtraction. Returns Matrix of the same form as the first parameter. +template +Form1 operator- ( - const Matrix& A, - const Matrix& B + const Matrix& A, + const Matrix& B ); - -//- Scalar multiplication of a matrix +//- Scalar multiplication of Matrix template Form operator* ( - const scalar s, + const Type& s, const Matrix& mat ); +//- Scalar addition of Matrix +template +Form operator+ +( + const Type& s, + const Matrix& mat +); -//- Scalar multiplication of a matrix +//- Scalar subtraction of Matrix +template +Form operator- +( + const Type& s, + const Matrix& mat +); + +//- Scalar multiplication of Matrix template Form operator* ( const Matrix& mat, - const scalar s + const Type& s ); +//- Scalar addition of Matrix +template +Form operator+ +( + const Matrix& mat, + const Type& s +); -//- Scalar division of a matrix +//- Scalar subtraction of Matrix +template +Form operator- +( + const Matrix& mat, + const Type& s +); + +//- Scalar division of Matrix template Form operator/ ( const Matrix& mat, - const scalar s + const Type& s ); - -//- Matrix-matrix multiplication +//- Matrix-Matrix multiplication using ikj-algorithm template typename typeOfInnerProduct::type operator* @@ -534,7 +676,6 @@ operator* const Matrix& B ); - //- Matrix-vector multiplication (A * x), where x is a column vector template inline tmp> operator* @@ -551,8 +692,7 @@ inline tmp> operator* const IndirectListBase& x ); - -//- Vector-matrix multiplication (x * A), where x is a row vector +//- Vector-Matrix multiplication (x * A), where x is a row vector template inline tmp> operator* ( @@ -560,7 +700,7 @@ inline tmp> operator* const Matrix& mat ); -//- Vector-matrix multiplication (x * A), where x is a row vector +//- Vector-Matrix multiplication (x * A), where x is a row vector template inline tmp> operator* ( @@ -568,6 +708,24 @@ inline tmp> operator* const Matrix& mat ); +//- Implicit inner product of Matrix-Matrix, equivalent to A.T()*B +template +typename typeOfInnerProduct::type +operator& +( + const Matrix& ATranspose, + const Matrix& B +); + +//- Implicit outer product of Matrix-Matrix, equivalent to A*B.T() +template +typename typeOfInnerProduct::type +operator^ +( + const Matrix& A, + const Matrix& BTranspose +); + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/OpenFOAM/matrices/Matrix/MatrixI.H b/src/OpenFOAM/matrices/Matrix/MatrixI.H index d6e3d836c6..4c229fbd18 100644 --- a/src/OpenFOAM/matrices/Matrix/MatrixI.H +++ b/src/OpenFOAM/matrices/Matrix/MatrixI.H @@ -128,40 +128,36 @@ inline bool Foam::Matrix::empty() const noexcept template inline void Foam::Matrix::checki(const label i) const { - #ifdef FULLDEBUG if (!mRows_ || !nCols_) { FatalErrorInFunction << "Attempt to access element from empty matrix" << abort(FatalError); } - if (i < 0 || i >= mRows_) + if (i < 0 || mRows_ <= i) { FatalErrorInFunction << "Index " << i << " out of range 0 ... " << mRows_-1 << abort(FatalError); } - #endif } template inline void Foam::Matrix::checkj(const label j) const { - #ifdef FULLDEBUG if (!mRows_ || !nCols_) { FatalErrorInFunction << "Attempt to access element from empty matrix" << abort(FatalError); } - if (j < 0 || j >= nCols_) + if (j < 0 || nCols_ <= j) { FatalErrorInFunction << "index " << j << " out of range 0 ... " << nCols_-1 << abort(FatalError); } - #endif } @@ -238,7 +234,7 @@ template inline const Type& Foam::Matrix::at(const label idx) const { #ifdef FULLDEBUG - if (idx < 0 || idx >= this->size()) + if (idx < 0 || this->size() <= idx) { FatalErrorInFunction << "index " << idx << " out of range 0 ... " << this->size() @@ -253,7 +249,7 @@ template inline Type& Foam::Matrix::at(const label idx) { #ifdef FULLDEBUG - if (idx < 0 || idx >= this->size()) + if (idx < 0 || this->size() <= idx) { FatalErrorInFunction << "index " << idx << " out of range 0 ... " << this->size() @@ -266,21 +262,74 @@ inline Type& Foam::Matrix::at(const label idx) template inline Foam::ConstMatrixBlock> -Foam::Matrix::block +Foam::Matrix::subColumn ( - const label m, - const label n, - const label mStart, - const label nStart + const label colIndex, + const label rowIndex, + label len ) const { + if (len < 0) + { + len = mRows_ - rowIndex; + } + return ConstMatrixBlock ( *this, - m, - n, - mStart, - nStart + len, // rows + 1, + rowIndex, + colIndex + ); +} + + +template +inline Foam::ConstMatrixBlock> +Foam::Matrix::subRow +( + const label rowIndex, + const label colIndex, + label len +) const +{ + if (len < 0) + { + len = nCols_ - colIndex; + } + + return ConstMatrixBlock + ( + *this, + 1, + len, // columns + rowIndex, + colIndex + ); +} + + +template +inline Foam::ConstMatrixBlock> +Foam::Matrix::subMatrix +( + const label rowIndex, + const label colIndex, + label szRows, + label szCols +) const +{ + if (szRows < 0) szRows = mRows_ - rowIndex; + if (szCols < 0) szCols = nCols_ - colIndex; + + return ConstMatrixBlock + ( + *this, + szRows, + szCols, + rowIndex, + colIndex ); } @@ -290,8 +339,8 @@ template inline Foam::ConstMatrixBlock> Foam::Matrix::block ( - const label mStart, - const label nStart + const label rowIndex, + const label colIndex ) const { return ConstMatrixBlock @@ -299,68 +348,82 @@ Foam::Matrix::block *this, VectorSpace::mRows, VectorSpace::nCols, - mStart, - nStart - ); -} - - -template -inline Foam::ConstMatrixBlock> -Foam::Matrix::col -( - const label m, - const label mStart -) const -{ - return ConstMatrixBlock - ( - *this, - m, - 1, - mStart, - 0 - ); -} - - -template -inline Foam::ConstMatrixBlock> -Foam::Matrix::col -( - const label m, - const label mStart, - const label nStart -) const -{ - return ConstMatrixBlock - ( - *this, - m, - 1, - mStart, - nStart + rowIndex, + colIndex ); } template inline Foam::MatrixBlock> -Foam::Matrix::block +Foam::Matrix::subColumn ( - const label m, - const label n, - const label mStart, - const label nStart + const label colIndex, + const label rowIndex, + label len ) { + if (len < 0) + { + len = mRows_ - rowIndex; + } + return MatrixBlock ( *this, - m, - n, - mStart, - nStart + len, // rows + 1, + rowIndex, + colIndex + ); +} + + +template +inline Foam::MatrixBlock> +Foam::Matrix::subRow +( + const label rowIndex, + const label colIndex, + label len +) +{ + if (len < 0) + { + len = nCols_ - colIndex; + } + + return MatrixBlock + ( + *this, + 1, + len, // columns + rowIndex, + colIndex + ); +} + + +template +inline Foam::MatrixBlock> +Foam::Matrix::subMatrix +( + const label rowIndex, + const label colIndex, + label szRows, + label szCols +) +{ + if (szRows < 0) szRows = mRows_ - rowIndex; + if (szCols < 0) szCols = nCols_ - colIndex; + + return MatrixBlock + ( + *this, + szRows, + szCols, + rowIndex, + colIndex ); } @@ -368,50 +431,19 @@ Foam::Matrix::block template template inline Foam::MatrixBlock> -Foam::Matrix::block(const label mStart, const label nStart) +Foam::Matrix::block +( + const label rowIndex, + const label colIndex +) { return MatrixBlock ( *this, VectorSpace::mRows, VectorSpace::nCols, - mStart, - nStart - ); -} - - -template -inline Foam::MatrixBlock> -Foam::Matrix::col(const label m, const label mStart) -{ - return MatrixBlock - ( - *this, - m, - 1, - mStart, - 0 - ); -} - - -template -inline Foam::MatrixBlock> -Foam::Matrix::col -( - const label m, - const label mStart, - const label nStart -) -{ - return MatrixBlock - ( - *this, - m, - 1, - mStart, - nStart + rowIndex, + colIndex ); } @@ -532,8 +564,10 @@ inline const Type& Foam::Matrix::operator() const label jcol ) const { + #ifdef FULLDEBUG checki(irow); checkj(jcol); + #endif return v_[irow*nCols_ + jcol]; } @@ -545,8 +579,10 @@ inline Type& Foam::Matrix::operator() const label jcol ) { + #ifdef FULLDEBUG checki(irow); checkj(jcol); + #endif return v_[irow*nCols_ + jcol]; } @@ -554,7 +590,9 @@ inline Type& Foam::Matrix::operator() template inline const Type* Foam::Matrix::operator[](const label irow) const { + #ifdef FULLDEBUG checki(irow); + #endif return v_ + irow*nCols_; } @@ -562,7 +600,9 @@ inline const Type* Foam::Matrix::operator[](const label irow) const template inline Type* Foam::Matrix::operator[](const label irow) { + #ifdef FULLDEBUG checki(irow); + #endif return v_ + irow*nCols_; } diff --git a/src/OpenFOAM/matrices/Matrix/MatrixTools.C b/src/OpenFOAM/matrices/Matrix/MatrixTools.C new file mode 100644 index 0000000000..5ce736cdf3 --- /dev/null +++ b/src/OpenFOAM/matrices/Matrix/MatrixTools.C @@ -0,0 +1,136 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2019 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +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 . + +\*---------------------------------------------------------------------------*/ + +#include "MatrixTools.H" + +// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * * // + +template +bool Foam::MatrixTools::equal +( + const Matrix& A, + const Matrix& B, + const bool verbose, + const scalar relTol, + const scalar absTol +) +{ + const label len = A.size(); + + if (len != B.size()) + { + if (verbose) + { + Info<< "Matrices have different sizes: " + << len << " vs " << B.size() << nl; + } + return false; + } + + auto iter1 = A.cbegin(); + auto iter2 = B.cbegin(); + + for (label i = 0; i < len; ++i) + { + if ((absTol + relTol*mag(*iter2)) < Foam::mag(*iter1 - *iter2)) + { + if (verbose) + { + Info<< "Matrix element " << i + << " differs beyond tolerance: " + << *iter1 << " vs " << *iter2 << nl; + } + return false; + } + + ++iter1; + ++iter2; + } + + if (verbose) + { + Info<< "All elements equal within the tolerances" << nl; + } + + return true; +} + + +template +Foam::Ostream& Foam::MatrixTools::printMatrix +( + Ostream& os, + const Container& mat +) +{ + os << mat.m() << ' ' << mat.n(); + + if (mat.m() == 1) + { + // row-vector + os << " ("; + for (label j = 0; j < mat.n(); ++j) + { + if (j) os << ' '; + os << mat(0,j); + } + os << ')' << nl; + } + else if (mat.n() == 1) + { + // col-vector + + os << " ("; + for (label i = 0; i < mat.m(); ++i) + { + if (i) os << ' '; + os << mat(i,0); + } + os << ')' << nl; + } + else + { + // Regular + + os << nl << '(' << nl; + + for (label i = 0; i < mat.m(); ++i) + { + os << '('; + for (label j = 0; j < mat.n(); ++j) + { + if (j) os << ' '; + os << mat(i,j); + } + os << ')' << nl; + } + os << ')' << nl; + } + + return os; +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/Matrix/MatrixTools.H b/src/OpenFOAM/matrices/Matrix/MatrixTools.H new file mode 100644 index 0000000000..8c63efc0fe --- /dev/null +++ b/src/OpenFOAM/matrices/Matrix/MatrixTools.H @@ -0,0 +1,90 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2019 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +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 . + +Namespace + Foam::MatrixTools + +Description + Collection of static functions for matrix-related verifications. + +SourceFiles + MatrixTools.C + +\*---------------------------------------------------------------------------*/ + +#ifndef MatrixTools_H +#define MatrixTools_H + +#include "Matrix.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward declarations +class Ostream; + + +/*---------------------------------------------------------------------------*\ + Namespace MatrixTools Declaration +\*---------------------------------------------------------------------------*/ + +namespace MatrixTools +{ + +//- Compare matrix elements for absolute or relative equality +template +bool equal +( + const Matrix& A, + const Matrix& B, + const bool verbose = false, + const scalar relTol = 1e-5, + const scalar absTol = 1e-8 +); + + +//- Simple ASCII output of Matrix, MatrixBlock +template +Ostream& printMatrix(Ostream& os, const Container& mat); + + +} // End namespace MatrixTools + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository + #include "MatrixTools.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.C b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.C index 2784a342a3..7ced3ca8de 100644 --- a/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.C +++ b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | + \\ / A nd | Copyright (C) 2004-2010, 2019 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- | Copyright (C) 2016 OpenFOAM Foundation @@ -71,6 +71,76 @@ Foam::MatrixBlock::operator Field() const } +template +Foam::label Foam::ConstMatrixBlock::disallow +( + const char* what +) const +{ + FatalErrorInFunction + << "Block addresses " << what + << " outside matrix or invalid matrix components" + << abort(FatalError); + return 0; +} + + +template +Foam::label Foam::MatrixBlock::disallow +( + const char* what +) const +{ + FatalErrorInFunction + << "Block addresses " << what + << " outside matrix or invalid matrix components" + << abort(FatalError); + return 0; +} + + +template void Foam::ConstMatrixBlock::checkIndex +( + const label i, + const label j +) const +{ + if (i < 0 || i >= mRows_) + { + FatalErrorInFunction + << "Index " << i << " is out of range 0 ... " << mRows_ - 1 + << abort(FatalError); + } + else if (j < 0 || j >= nCols_) + { + FatalErrorInFunction + << "Index " << j << " is out of range 0 ... " << nCols_ - 1 + << abort(FatalError); + } +} + + +template void Foam::MatrixBlock::checkIndex +( + const label i, + const label j +) const +{ + if (i < 0 || i >= mRows_) + { + FatalErrorInFunction + << "Index " << i << " is out of range 0 ... " << mRows_ - 1 + << abort(FatalError); + } + else if (j < 0 || j >= nCols_) + { + FatalErrorInFunction + << "Index " << j << " is out of range 0 ... " << nCols_ - 1 + << abort(FatalError); + } +} + + // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // template diff --git a/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.H b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.H index 27e90c7a9f..bb3cf1d95e 100644 --- a/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.H +++ b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2019 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2010, 2019 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- | Copyright (C) 2016 OpenFOAM Foundation @@ -65,14 +65,21 @@ class ConstMatrixBlock //- Const reference to the parent matrix const MatrixType& matrix_; - // Block size + //- Block size const label mRows_; const label nCols_; - // Block location in parent matrix + //- Block location in parent matrix const label rowStart_; const label colStart_; + + // Private Member Functions + + //- Error message for failed sanity checks during matrix construction + label disallow(const char *what) const; + + public: typedef typename MatrixType::cmptType cmptType; @@ -110,6 +117,9 @@ public: //- Convert a column of a matrix to a Field operator Field() const; + + //- Check if (i, j) is within range of row-column limits + void checkIndex(const label i, const label j) const; }; @@ -125,14 +135,21 @@ class MatrixBlock //- Reference to the parent matrix MatrixType& matrix_; - // Block size + //- Block size const label mRows_; const label nCols_; - // Block location in parent matrix + //- Block location in parent matrix const label rowStart_; const label colStart_; + + // Private Member Functions + + //- Error message for failed sanity checks during matrix construction + label disallow(const char *what) const; + + public: typedef typename MatrixType::cmptType cmptType; @@ -174,8 +191,11 @@ public: //- Convert a column of a matrix to a Field operator Field() const; + //- Check if (i, j) is within range of row-column limits + void checkIndex(const label i, const label j) const; - // Member operators + + // Member Operators //- Assignment to a compatible matrix template diff --git a/src/OpenFOAM/matrices/MatrixBlock/MatrixBlockI.H b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlockI.H index 408f1393a1..e6a47e0ba4 100644 --- a/src/OpenFOAM/matrices/MatrixBlock/MatrixBlockI.H +++ b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlockI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2019 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2010, 2019 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- | Copyright (C) 2016 OpenFOAM Foundation @@ -38,24 +38,21 @@ Foam::ConstMatrixBlock::ConstMatrixBlock ) : matrix_(matrix), - mRows_(m), - nCols_(n), - rowStart_(mStart), - colStart_(nStart) -{ - #ifdef FULLDEBUG - if + mRows_(0 < m ? m : disallow("row dim")), + nCols_(0 < n ? n : disallow("col dim")), + rowStart_ ( - rowStart_ + mRows_ > matrix.m() - || colStart_ + nCols_ > matrix.n() + 0 <= mStart + || mStart + mRows_ <= matrix.m() + ? mStart : disallow("row start") + ), + colStart_ + ( + 0 <= nStart + || nStart + nCols_ <= matrix.n() + ? nStart : disallow("col start") ) - { - FatalErrorInFunction - << "Block addresses outside matrix" - << abort(FatalError); - } - #endif -} +{} template @@ -69,24 +66,21 @@ Foam::MatrixBlock::MatrixBlock ) : matrix_(matrix), - mRows_(m), - nCols_(n), - rowStart_(mStart), - colStart_(nStart) -{ - #ifdef FULLDEBUG - if + mRows_(0 < m ? m : disallow("row dim")), + nCols_(0 < n ? n : disallow("col dim")), + rowStart_ ( - rowStart_ + mRows_ > matrix.m() - || colStart_ + nCols_ > matrix.n() + 0 <= mStart + || mStart + mRows_ <= matrix.m() + ? mStart : disallow("row start") + ), + colStart_ + ( + 0 <= nStart + || nStart + nCols_ <= matrix.n() + ? nStart : disallow("col start") ) - { - FatalErrorInFunction - << "Block addresses outside matrix" - << abort(FatalError); - } - #endif -} +{} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // @@ -144,18 +138,7 @@ Foam::ConstMatrixBlock::operator() ) const { #ifdef FULLDEBUG - if (i<0 || i>=mRows_) - { - FatalErrorInFunction - << "Index " << i << " out of range 0 ... " << mRows_-1 - << abort(FatalError); - } - if (j<0 || j>=nCols_) - { - FatalErrorInFunction - << "Index " << j << " out of range 0 ... " << nCols_-1 - << abort(FatalError); - } + checkIndex(i, j); #endif return matrix_(i + rowStart_, j + colStart_); @@ -171,18 +154,7 @@ Foam::MatrixBlock::operator() ) const { #ifdef FULLDEBUG - if (i<0 || i>=mRows_) - { - FatalErrorInFunction - << "Index " << i << " out of range 0 ... " << mRows_-1 - << abort(FatalError); - } - if (j<0 || j>=nCols_) - { - FatalErrorInFunction - << "Index " << j << " out of range 0 ... " << nCols_-1 - << abort(FatalError); - } + checkIndex(i, j); #endif return matrix_(i + rowStart_, j + colStart_); @@ -198,18 +170,7 @@ Foam::MatrixBlock::operator() ) { #ifdef FULLDEBUG - if (i<0 || i>=mRows_) - { - FatalErrorInFunction - << "Index " << i << " out of range 0 ... " << mRows_-1 - << abort(FatalError); - } - if (j<0 || j>=nCols_) - { - FatalErrorInFunction - << "Index " << j << " out of range 0 ... " << nCols_-1 - << abort(FatalError); - } + checkIndex(i, j); #endif return matrix_(i + rowStart_, j + colStart_); diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C index ea3c437a86..624db7b7dc 100644 --- a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C +++ b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C @@ -282,7 +282,7 @@ Foam::QRMatrix::inv() const x[j] = Q_(i, j); } solvex(x); - inv.block(m, 1, 0, i) = x; + inv.subColumn(i) = x; } return inv; diff --git a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H index 68ecd56a59..668e0241b0 100644 --- a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H +++ b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H @@ -118,7 +118,7 @@ public: }; -// Global functions and operators +// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * // template class typeOfInnerProduct, RectangularMatrix> @@ -128,6 +128,7 @@ public: typedef RectangularMatrix type; }; + template class typeOfInnerProduct, SquareMatrix> { @@ -136,6 +137,7 @@ public: typedef RectangularMatrix type; }; + template class typeOfInnerProduct, RectangularMatrix> { @@ -145,6 +147,8 @@ public: }; +// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * // + template RectangularMatrix outer(const Field& f1, const Field& f2); diff --git a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H index da3e03b53f..1ff3147371 100644 --- a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H +++ b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H @@ -190,9 +190,9 @@ inline Foam::RectangularMatrix outer { RectangularMatrix f1f2T(f1.size(), f2.size()); - for (label i=0; i::one; - for (label i=0; i @@ -171,7 +171,6 @@ public: typedef SquareMatrix type; }; - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H index 4e1b3ca483..47897b61dc 100644 --- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H +++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H @@ -251,9 +251,9 @@ inline Foam::SquareMatrix symmOuter { SquareMatrix f1f2T(f1.size()); - for (label i=0; i