diff --git a/applications/test/Matrix/Test-Matrix.C b/applications/test/Matrix/Test-Matrix.C
index 4e001d7c6ac60a12901104e63a38f753417cf5e2..dbccf225fabf27df0105deb71fa58e9740d6b051 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 <algorithm>
using namespace Foam;
+using namespace Foam::MatrixTools;
+
+#define isEqual MatrixTools::equal
+
+void horizontalLine()
+{
+ Info<< "+---------+---------+---------+---------+---------+" << nl;
+}
+
// Copy values into matrix
template<class Form, class Type>
@@ -59,254 +65,906 @@ void assignMatrix
std::copy(list.begin(), list.end(), mat.begin());
}
-// Create matrix with values
+
template<class MatrixType>
-MatrixType makeMatrix
+MatrixType makeRandomMatrix
(
const labelPair& dims,
- std::initializer_list<typename MatrixType::cmptType> list
+ Random& rndGen
)
{
MatrixType mat(dims);
- assignMatrix(mat, list);
+ std::generate
+ (
+ mat.begin(),
+ mat.end(),
+ [&]{return rndGen.GaussNormal<typename MatrixType::cmptType>();}
+ );
return mat;
}
-// Create matrix with values
-template<class MatrixType, Foam::label nRows, Foam::label nCols>
-MatrixType makeMatrix
-(
- std::initializer_list<typename MatrixType::cmptType> list
-)
+template<class MatrixType>
+Ostream& operator<<(Ostream& os, const ConstMatrixBlock<MatrixType>& mat)
{
- MatrixType mat(labelPair(nRows, nCols));
+ return MatrixTools::printMatrix(os, mat);
+}
- assignMatrix(mat, list);
- return mat;
+template<class MatrixType>
+Ostream& operator<<(Ostream& os, const MatrixBlock<MatrixType>& mat)
+{
+ return MatrixTools::printMatrix(os, mat);
}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-// Main program:
+// * * * * * * * * * * * * * * * Main Program * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
- // SquareMatrix
+ typedef SquareMatrix<scalar> SMatrix;
+ typedef SquareMatrix<complex> SCMatrix;
+ typedef RectangularMatrix<scalar> RMatrix;
+ typedef RectangularMatrix<complex> RCMatrix;
+
+ Random rndGen(1234);
+
+ #if 1
+ {
+ horizontalLine();
+
+ Info<< "## Matrix Constructors & Assignments:" << nl << nl;
+
+ Info<< "# SquareMatrix<scalar> examples:" << nl;
+ {
+ const label mRows = 2;
+ const label nCols = 2;
+ Matrix<SMatrix, scalar> MS0(mRows, nCols);
+ Matrix<SMatrix, scalar> MS1(mRows, nCols, Zero);
+ Matrix<SMatrix, scalar> MS2(mRows, nCols, 17.44);
+ Matrix<SMatrix, scalar> 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<scalar> examples:" << nl;
+ {
+ const label mRows = 2;
+ const label nCols = 3;
+ Matrix<RMatrix, scalar> MR0(mRows, nCols);
+ Matrix<RMatrix, scalar> MR1(mRows, nCols, Zero);
+ Matrix<RMatrix, scalar> 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<complex> examples:" << nl;
+ {
+ const label mRows = 2;
+ const label nCols = 2;
+ Matrix<SCMatrix, complex> MS0(mRows, nCols);
+ Matrix<SCMatrix, complex> MS1(mRows, nCols, Zero);
+ Matrix<SCMatrix, complex> MS2(mRows, nCols, complex(17.44, 44.17));
+ Matrix<SCMatrix, complex> 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
+ (
+ S,
+ {
+ -3, 11, -4,
+ 2, 3, 10,
+ 2, 6, 1
+ }
+ );
+ S(0,1) = 10;
+
+ Info<< "# SquareMatrix<scalar> 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<< "# RectangularMatrix<scalar> example:" << nl;
+ printMatrix(Info, R) << nl;
+
+ // 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
{
- Info<< nl << "Test SquareMatrix" << nl;
+ 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;
- SquareMatrix<scalar> square1(3);
+ 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<scalar> 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
(
- square1,
+ RC,
{
- -3.0, 10.0, -4.0,
- 2.0, 3.0, 10.0,
- 2.0, 6.0, 1.0
+ complex(1, 0), complex(2,5),
+ complex(4, 3), complex(1,-1),
+ complex(1, 2), complex(2,9)
}
);
- Info<< "matrix: " << square1
- << " begin: " << uintptr_t(square1.cdata()) << nl;
+ 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);
+ {
+ List<scalar> diagR = R1.diag();
+ scalar traceR = R1.trace();
+ Info<< "RectangularMatrix<scalar>:"
+ << 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<scalar> lst(5, 2.0);
+ S1.diag(lst);
+ Info<< "After List assignment:" << nl;
+ printMatrix(Info, S1);
+ }
- // Test makeMatrix
+ Info<< "# Diagonal sort of Matrix:" << nl;
{
- auto square1b
+ auto S1
(
- makeMatrix<scalarSquareMatrix>
+ makeRandomMatrix<SMatrix>
(
- {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<scalar> diag = S1.diag();
+ SortList<scalar> incrDiag(diag);
+ incrDiag.sort();
+
+ Info<< "Diagonal list:" << diag << nl
+ << "Ascending diagonal list: " << incrDiag << nl << nl;
+ }
+
+ horizontalLine();
+ }
+ #endif
+
+
+ #if 1
+ {
+ horizontalLine();
+
+ Info<< "## Transpose Verifications:" << nl;
+
+ const label numberOfTests = 10;
+
+ for (label i = 0; i < numberOfTests; ++i)
+ {
+ const label mRows = rndGen.position(1,10);
+ const label nCols = rndGen.position(1,10);
+
+ Info << "# SquareMatrix<scalar> example:" << nl;
+ {
+ Info<< "A(mxm) where m = " << mRows << nl;
+ auto A(makeRandomMatrix<SMatrix>({mRows, mRows}, rndGen));
+ auto B(makeRandomMatrix<SMatrix>({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<scalar> example:" << nl;
+ {
+ Info<< "A(mxn) where"
+ << " m = " << mRows << ","
+ << " n = " << nCols << nl;
+ auto A(makeRandomMatrix<RMatrix>({mRows, nCols}, rndGen));
+ auto B(makeRandomMatrix<RMatrix>({mRows, nCols}, rndGen));
+ auto C(makeRandomMatrix<RMatrix>({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<complex> example:" << nl;
+ {
+ Info<< "A(mxm) where m = " << mRows << nl;
+ SCMatrix A(mRows);
+ SCMatrix B(mRows);
+
+ for (auto& val : A)
+ {
+ val.Re() = rndGen.GaussNormal<scalar>();
+ val.Im() = rndGen.GaussNormal<scalar>();
+ }
+ for (auto& val : B)
+ {
+ val.Re() = rndGen.GaussNormal<scalar>();
+ val.Im() = rndGen.GaussNormal<scalar>();
+ }
+
+ 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
- auto square1c
+
+ #if 1
+ {
+ horizontalLine();
+
+ Info<< "## Scalar Arithmetic Operations:" << nl;
+
+ if (true)
+ {
+ SMatrix S1(3, Zero);
+ assignMatrix
(
- makeMatrix<scalarSquareMatrix, 3, 3>
- ({
- -3.0, 10.0, -4.0,
- 2.0, 3.0, 10.0,
- 2.0, 6.0, 1.0
- })
+ 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<< "makeMatrix: " << square1c << nl;
- }
+ Info<< "SquareMatrix<scalar> S1 = " << S1 << nl;
+ Info<< "SquareMatrix<scalar> S2 = " << S2 << nl;
+ Info<< "SquareMatrix<scalar> S3 = " << S3 << nl;
- //Info<< square1 - 2.0*(-square1) << nl;
- Info<< "min:" << min(square1) << " max:" << max(square1) << nl;
- Info<< "min/max: " << minMax(square1) << nl;
+ S3 = S1*S2;
+ Info<< "S1*S2 = " << S3 << nl;
- // Steal matrix contents
+ S3 = S1 - S2;
+ Info<< "S1 - S2 = " << S3 << nl;
- List<scalar> stole(square1.release());
+ S3 = S1 + S2;
+ Info<< "S1 + S2 = " << S3 << nl;
- Info<< "matrix: " << square1
- << " begin: " << uintptr_t(square1.cdata()) << nl;
+// S3 *= S1; // Not Implemented
- Info<< "List: " << stole << nl;
+ S3 -= S1;
+ Info<< "S3 -= S1; S3 = " << S3 << nl;
- Info<< "min/max: " << minMax(square1) << nl;
+ S3 += S1;
+ Info<< "S3 += S1; S3 = " << S3 << nl;
- square1 = 100;
- Info<< "matrix: " << square1
- << " begin: " << uintptr_t(square1.cdata()) << nl;
+ Info<< nl << "# Scalar broadcasting:" << nl;
+ S3 *= 5.0;
+ Info<< "S3 *= 5.0; S3 = " << S3 << nl;
- SquareMatrix<scalar> square2(3, I);
+ S3 /= 5.0;
+ Info<< "S3 /= 5.0; S3 = " << S3 << nl;
- square1 = square2;
+ S3 -= 1.0;
+ Info<< "S3 -= 1.0; S3 = " << S3 << nl;
- Info<< nl << "After copy assign from Identity:" << nl << square1 << nl;
+ S3 += 1.0;
+ Info<< "S3 += 1.0; S3 = " << S3 << nl;
- square1 += 1.25*square2;
- Info<< nl << "After +=" << nl << square1 << nl;
+ Info<< nl << "# Operations between different matrix types:" << nl;
+ RMatrix R1 = S1;
+ Info<< "RectangularMatrix<scalar> R1 = " << R1 << nl;
- square1 -= square2.T();
+ // RectangularMatrix*SquareMatrix returns RectangularMatrix
+ // S3 = S3*R1; // Not implemented
+ R1 = S3*R1;
+ Info<< "R1 = S3*R1; R1 = " << R1 << nl;
- Info<< nl << "After -=" << nl << square1 << nl;
+ S3 = S3 - R1;
+ Info<< "S3 = S3 - R1; S3 = " << S3 << nl;
- square1 *= 10;
+ S3 = S3 + R1;
+ Info<< "S3 = S3 + R1; S3 = " << S3 << nl;
- Info<< nl << "After *=" << nl << square1 << nl;
+ R1 = S3 + R1;
+ Info<< "R1 = S3 + R1; R1 = " << R1 << nl;
- square1 /= 8;
- Info<< nl << "After /=" << nl << square1 << nl;
+ Info<< nl << "# Extrema operations:" << nl;
- SquareMatrix<scalar> square4;
- square4 = square2;
+ Info<< "min(S3) = " << min(S3) << nl
+ << "max(S3) = " << max(S3) << nl
+ << "minMax(S3) = " << minMax(S3) << nl;
+ }
- Info<< nl << square4 << endl;
+ if (true)
+ {
+ Info<< nl << "# Operations with ListTypes<scalar>:" << nl;
- SquareMatrix<scalar> square5;
- square4 = square5;
- Info<< nl << square5 << endl;
- }
+ 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<scalar> lst({1, 2, 3});
+ RMatrix colVector(3, 1, Zero);
+ assignMatrix
+ (
+ colVector,
+ {1, 2, 3}
+ );
+ RMatrix rowVector(1, 3, Zero);
+ assignMatrix
+ (
+ rowVector,
+ {1, 2, 3}
+ );
- // RectangularMatrix
- {
- Info<< nl << "Test RectangularMatrix" << nl;
+ Info<< "A = " << A << nl;
+ Info<< "colVector = " << colVector << nl;
+ Info<< "rowVector = " << rowVector << nl;
+
+ const Field<scalar> field1(A.Amul(lst));
+ const Field<scalar> field2(A*lst);
+ const Field<scalar> field3(A.Tmul(lst));
+ const Field<scalar> 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);
- RectangularMatrix<scalar> rm1({5, 6}, 3.1);
- rm1(0, 1) = 4.5;
+ 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);
- RectangularMatrix<scalar> rm1b(rm1.block(2, 2, 0, 0));
+ 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 // << "Full matrix " << rm1 << nl
- << "block = " << rm1b << endl;
+ Info<< 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<< "## Complex Arithmetic Operations:" << nl;
- Info<< "Symmetric Square Matrix = " << symmMatrix << endl;
+ 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<< "Inverse = " << inv(symmMatrix) << endl;
- Info<< "Determinant = " << det(symmMatrix) << endl;
+ Info<< "SquareMatrix<complex> S1 = " << S1 << nl;
+ Info<< "SquareMatrix<complex> S2 = " << S2 << nl;
+ Info<< "SquareMatrix<complex> 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);
- rhs[0] = 1;
- rhs[1] = 2;
- rhs[2] = 3;
+ S3 = S1 + S2;
+ Info<< "S1 + S2 = " << S3 << nl;
- LUsolve(symmMatrix, rhs);
+// S3 *= S1; // Not Implemented
- Info<< "Decomposition = " << symmMatrix << endl;
- Info<< "Solution = " << rhs << endl;
- }
+ S3 -= S1;
+ Info<< "S3 -= S1; S3 = " << S3 << nl;
+ S3 += S1;
+ Info<< "S3 += S1; S3 = " << S3 << nl;
- 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 << "# Complex broadcasting:" << nl;
- Info<< nl << "Square Matrix = " << squareMatrix << endl;
+ S3 *= complex(5, 0);
+ Info<< "S3 *= complex(5, 0); S3 = " << S3 << nl;
- const scalarField source(3, 1);
+ 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<complex> 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)
{
- scalarSquareMatrix sm(squareMatrix);
- Info<< "det = " << det(sm) << endl;
+ Info<< nl << "# Operations with ListTypes<complex>:" << 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<complex> 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<complex> field1(A.Amul(lst));
+ const Field<complex> field2(A*lst);
+ const Field<complex> field3(A.Tmul(lst));
+ const Field<complex> 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;
}
- scalarSquareMatrix sm(squareMatrix);
- labelList rhs(3, 0);
- label sign;
- LUDecompose(sm, rhs, sign);
+ #if 1
+ {
+ Info<< nl << "# Implicit inner/outer products:" << nl;
- Info<< "Decomposition = " << sm << endl;
- Info<< "Pivots = " << rhs << endl;
- Info<< "Sign = " << sign << endl;
- Info<< "det = " << detDecomposed(sm, sign) << endl;
- }
+ 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);
- {
- 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;
+ 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
{
- LLTMatrix<scalar> LLT(squareMatrix);
- scalarField x(LLT.solve(source));
- Info<< "LLT solve residual " << (squareMatrix*x - source) << endl;
+ horizontalLine();
+
+ Info<< "## SymmetricSquareMatrix<scalar> 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(symm, rhs);
+ Info<< "Solving linear system through LU decomposition:" << nl
+ << "Decomposition = " << nl << symm << nl
+ << "Solution = " << rhs << nl;
}
+ #endif
+
+ #if 1
{
- QRMatrix<scalarSquareMatrix> QR(squareMatrix);
- scalarField x(QR.solve(source));
+ scalarSquareMatrix squareMatrix(3, Zero);
- Info<< "QR solve residual "
- << (squareMatrix*x - source) << endl;
+ scalarSquareMatrix S(3, Zero);
+ assignMatrix
+ (
+ S,
+ {
+ 4, 12, -16,
+ 12, 37, -43,
+ -16, -43, 98
+ }
+ );
+
+ const scalarField source(3, 1);
- Info<< "QR inverse solve residual "
- << (x - QR.inv()*source) << endl;
+ Info<< nl << "Square Matrix = " << S << endl;
- Info<< "QR inv *squareMatrix " << (QR.inv()*squareMatrix) << endl;
+ if (true)
+ {
+ {
+ 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;
+ }
+
+ 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;
+ }
+
+ if (true)
+ {
+ LLTMatrix<scalar> LLT(S);
+ scalarField x(LLT.solve(source));
+ Info<< "LLT solve residual " << (S*x - source) << nl;
+ }
+
+ horizontalLine();
}
+ #endif
+
- 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 1b0908dc3c665d46d3c83b7ecd14343e1271f4fe..072c78ffe7e51fa0e97cf1583c8d7f4978271096 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.C
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.C
@@ -334,16 +334,29 @@ void Foam::Matrix<Form, Type>::resize(const label m, const label n)
}
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::round(const scalar tol)
+{
+ for (Type& val : *this)
+ {
+ if (mag(val) < tol)
+ {
+ val = Zero;
+ }
+ }
+}
+
+
template<class Form, class Type>
Form Foam::Matrix<Form, Type>::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<Form, Type>::T() const
}
+template<class Form, class Type>
+Foam::List<Type> Foam::Matrix<Form, Type>::diag() const
+{
+ const label len = Foam::min(mRows_, nCols_);
+
+ List<Type> result(len);
+
+ for (label i=0; i < len; ++i)
+ {
+ result[i] = (*this)(i, i);
+ }
+
+ return result;
+}
+
+
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::diag(const UList<Type>& 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<class Form, class Type>
+Type Foam::Matrix<Form, Type>::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<class Form, class Type>
+Foam::scalar Foam::Matrix<Form, Type>::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<class Form, class Type>
+Foam::scalar Foam::Matrix<Form, Type>::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<class Form, class Type>
@@ -416,9 +514,9 @@ void Foam::Matrix<Form, Type>::operator=
const MatrixBlock<MatrixType>& 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<Form, Type>::operator=(const zero)
template<class Form, class Type>
void Foam::Matrix<Form, Type>::operator+=(const Matrix<Form, Type>& other)
{
+ #ifdef FULLDEBUG
if (this == &other)
{
FatalErrorInFunction
@@ -458,15 +557,13 @@ void Foam::Matrix<Form, Type>::operator+=(const Matrix<Form, Type>& 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<Form, Type>::operator+=(const Matrix<Form, Type>& other)
template<class Form, class Type>
void Foam::Matrix<Form, Type>::operator-=(const Matrix<Form, Type>& other)
{
+ #ifdef FULLDEBUG
if (this == &other)
{
FatalErrorInFunction
@@ -489,21 +587,39 @@ void Foam::Matrix<Form, Type>::operator-=(const Matrix<Form, Type>& other)
<< other.m() << ", " << other.n() << ')' << nl
<< abort(FatalError);
}
+ #endif
+
+ auto iter2 = other.cbegin();
+ for (Type& val : *this)
+ {
+ val -= *iter2;
+ ++iter2;
+ }
+}
- Type* out = this->data();
- const Type* in = other.cdata();
- const label len = this->size();
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::operator+=(const Type& s)
+{
+ for (Type& val : *this)
+ {
+ val += s;
+ }
+}
+
- for (label idx=0; idx < len; ++idx)
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::operator-=(const Type& s)
+{
+ for (Type& val : *this)
{
- out[idx] -= in[idx];
+ val -= s;
}
}
template<class Form, class Type>
-void Foam::Matrix<Form, Type>::operator*=(const scalar s)
+void Foam::Matrix<Form, Type>::operator*=(const Type& s)
{
for (Type& val : *this)
{
@@ -513,7 +629,7 @@ void Foam::Matrix<Form, Type>::operator*=(const scalar s)
template<class Form, class Type>
-void Foam::Matrix<Form, Type>::operator/=(const scalar s)
+void Foam::Matrix<Form, Type>::operator/=(const Type& s)
{
for (Type& val : *this)
{
@@ -522,7 +638,7 @@ void Foam::Matrix<Form, Type>::operator/=(const scalar s)
}
-// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * * //
template<class Form, class Type>
const Type& Foam::max(const Matrix<Form, Type>& mat)
@@ -569,7 +685,7 @@ Foam::MinMax<Type> Foam::minMax(const Matrix<Form, Type>& mat)
template<class Form, class Type>
Form Foam::operator-(const Matrix<Form, Type>& mat)
{
- Form result(mat.m(), mat.n());
+ Form result(mat.sizes());
std::transform
(
@@ -583,9 +699,14 @@ Form Foam::operator-(const Matrix<Form, Type>& mat)
}
-template<class Form, class Type>
-Form Foam::operator+(const Matrix<Form, Type>& A, const Matrix<Form, Type>& B)
+template<class Form1, class Form2, class Type>
+Form1 Foam::operator+
+(
+ const Matrix<Form1, Type>& A,
+ const Matrix<Form2, Type>& B
+)
{
+ #ifdef FULLDEBUG
if (A.m() != B.m() || A.n() != B.n())
{
FatalErrorInFunction
@@ -594,27 +715,31 @@ Form Foam::operator+(const Matrix<Form, Type>& A, const Matrix<Form, Type>& B)
<< B.m() << ", " << B.n() << ')' << nl
<< abort(FatalError);
}
+ #endif
- Form AB(A.m(), A.n());
-
- Type* ABv = AB.data();
- const Type* Av = A.cdata();
- const Type* Bv = B.cdata();
-
- const label len = A.size();
+ Form1 result(A.sizes());
- for (label idx = 0; idx < len; ++idx)
- {
- ABv[idx] = Av[idx] + Bv[idx];
- }
+ std::transform
+ (
+ A.cbegin(),
+ A.cend(),
+ B.cbegin(),
+ result.begin(),
+ std::plus<Type>()
+ );
- return AB;
+ return result;
}
-template<class Form, class Type>
-Form Foam::operator-(const Matrix<Form, Type>& A, const Matrix<Form, Type>& B)
+template<class Form1, class Form2, class Type>
+Form1 Foam::operator-
+(
+ const Matrix<Form1, Type>& A,
+ const Matrix<Form2, Type>& B
+)
{
+ #ifdef FULLDEBUG
if (A.m() != B.m() || A.n() != B.n())
{
FatalErrorInFunction
@@ -623,85 +748,117 @@ Form Foam::operator-(const Matrix<Form, Type>& A, const Matrix<Form, Type>& 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::minus<Type>()
+ );
- const label len = A.size();
+ return result;
+}
- for (label idx=0; idx < len; ++idx)
- {
- ABv[idx] = Av[idx] - Bv[idx];
- }
- return AB;
+template<class Form, class Type>
+Form Foam::operator*(const Type& s, const Matrix<Form, Type>& mat)
+{
+ Form result(mat.sizes());
+
+ std::transform
+ (
+ mat.cbegin(),
+ mat.cend(),
+ result.begin(),
+ [&](const Type& val) { return s * val; }
+ );
+
+ return result;
}
template<class Form, class Type>
-Form Foam::operator*(const scalar s, const Matrix<Form, Type>& mat)
+Form Foam::operator+(const Type& s, const Matrix<Form, Type>& mat)
{
- Form result(mat.m(), mat.n());
+ Form result(mat.sizes());
+
+ std::transform
+ (
+ mat.cbegin(),
+ mat.cend(),
+ result.begin(),
+ [&](const Type& val) { return s + val; }
+ );
- const label len = mat.size();
+ return result;
+}
- if (len)
- {
- Type* out = result.data();
- const Type* in = mat.cdata();
- for (label idx = 0; idx < len; ++idx)
- {
- out[idx] = s * in[idx];
- }
- }
+template<class Form, class Type>
+Form Foam::operator-(const Type& s, const Matrix<Form, Type>& mat)
+{
+ Form result(mat.sizes());
+
+ std::transform
+ (
+ mat.cbegin(),
+ mat.end(),
+ result.begin(),
+ [&](const Type& val) { return s - val; }
+ );
return result;
}
template<class Form, class Type>
-Form Foam::operator*(const Matrix<Form, Type>& mat, const scalar s)
+Form Foam::operator*(const Matrix<Form, Type>& mat, const Type& s)
{
- Form result(mat.m(), mat.n());
+ return s*mat;
+}
- const label len = mat.size();
- if (len)
- {
- Type* out = result.data();
- const Type* in = mat.cdata();
+template<class Form, class Type>
+Form Foam::operator+(const Matrix<Form, Type>& mat, const Type& s)
+{
+ return s + mat;
+}
- for (label idx=0; idx < len; ++idx)
- {
- out[idx] = in[idx] * s;
- }
- }
+
+template<class Form, class Type>
+Form Foam::operator-(const Matrix<Form, Type>& 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<class Form, class Type>
-Form Foam::operator/(const Matrix<Form, Type>& mat, const scalar s)
+Form Foam::operator/(const Matrix<Form, Type>& mat, const Type& s)
{
- Form result(mat.m(), mat.n());
-
- const label len = mat.size();
-
- if (len)
- {
- Type* out = result.data();
- const Type* in = mat.cdata();
+ Form result(mat.sizes());
- for (label idx=0; idx < len; ++idx)
- {
- out[idx] = in[idx] / s;
- }
- }
+ std::transform
+ (
+ mat.cbegin(),
+ mat.end(),
+ result.begin(),
+ [&](const Type& val) { return val / s; }
+ );
return result;
}
@@ -715,6 +872,7 @@ Foam::operator*
const Matrix<Form2, Type>& 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, Form1, Form2>::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 k = 0; k < B.m(); ++k)
+ {
+ for (label j = 0; j < AB.n(); ++j)
+ {
+ AB(i, j) += A(i, k)*B(k, j);
+ }
+ }
+ }
+
+ return AB;
+}
+
+
+template<class Form1, class Form2, class Type>
+typename Foam::typeOfInnerProduct<Type, Form1, Form2>::type
+Foam::operator&
+(
+ const Matrix<Form1, Type>& AT,
+ const Matrix<Form2, Type>& 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, Form1, Form2>::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<class Form1, class Form2, class Type>
+typename Foam::typeOfInnerProduct<Type, Form1, Form2>::type
+Foam::operator^
+(
+ const Matrix<Form1, Type>& A,
+ const Matrix<Form2, Type>& 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, Form1, Form2>::type AB
+ (
+ A.m(),
+ BT.m(),
+ Zero
+ );
+
+ for (label i = 0; i < AB.m(); ++i)
{
- for (label j=0; j < AB.n(); ++j)
+ for (label k = 0; k < BT.n(); ++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)*Detail::conj(BT(j, k));
}
}
}
diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.H b/src/OpenFOAM/matrices/Matrix/Matrix.H
index bfe1c4e51fecf2d56cad216f2ae09f79e0ce199d..2aef3908967f739eee8ca85a2a04e9d2e016b3de 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<class ListType>
tmp<Field<Type>> 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<class ListType>
tmp<Field<Type>> 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<Form, Type>&& mat);
- //- Copy constructor from matrix of a different form
+ //- Copy constructor from Matrix of a different form
template<class Form2>
explicit Matrix(const Matrix<Form2, Type>& mat);
- //- Construct from a block of another matrix
+ //- Construct from a block of another Matrix
template<class MatrixType>
Matrix(const ConstMatrixBlock<MatrixType>& Mb);
- //- Construct from a block of another matrix
+ //- Construct from a block of another Matrix
template<class MatrixType>
Matrix(const MatrixBlock<MatrixType>& 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)
+
+ //- 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<mType> subColumn
+ (
+ const label colIndex,
+ const label rowIndex = 0,
+ label len = -1
+ ) 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<mType> subRow
+ (
+ const label rowIndex,
+ const label colIndex = 0,
+ label len = -1
+ ) const;
+
+ //- Return const sub-block of Matrix
+ // Sub-block starting at columnIndex & rowIndex indices
+ inline ConstMatrixBlock<mType> subMatrix
+ (
+ const label rowIndex,
+ const label colIndex,
+ label szRows = -1,
+ label szCols = -1
+ ) const;
+
+ //- Access Field as a ConstMatrixBlock
+ template<class VectorSpace>
+ inline ConstMatrixBlock<mType> block
+ (
+ const label rowIndex,
+ const label colIndex
+ ) const;
+
+
+ // Block Access (non-const)
+
+ //- Return column or column's subset of Matrix
+ inline MatrixBlock<mType> subColumn
+ (
+ const label colIndex,
+ const label rowIndex = 0,
+ label len = -1
+ );
+
+ //- Return row or row's subset of Matrix
+ inline MatrixBlock<mType> subRow
+ (
+ const label rowIndex,
+ const label colIndex = 0,
+ label len = -1
+ );
+
+ //- Return sub-block of Matrix
+ inline MatrixBlock<mType> subMatrix
+ (
+ const label rowIndex,
+ const label colIndex,
+ label szRows = -1,
+ label szCols = -1
+ );
+
+ //- Access Field as a MatrixBlock
+ template<class VectorSpace>
+ inline MatrixBlock<mType> block
+ (
+ const label rowIndex,
+ const label colIndex
+ );
- // Block Access
-
- inline ConstMatrixBlock<mType> block
- (
- const label m,
- const label n,
- const label mStart,
- const label nStart
- ) const;
-
- template<class VectorSpace>
- inline ConstMatrixBlock<mType> block
- (
- const label mStart,
- const label nStart
- ) const;
-
- inline ConstMatrixBlock<mType> col
- (
- const label m,
- const label rowStart
- ) const;
-
- inline ConstMatrixBlock<mType> col
- (
- const label m,
- const label mStart,
- const label nStart
- ) const;
-
-
- inline MatrixBlock<mType> block
- (
- const label m,
- const label n,
- const label mStart,
- const label nStart
- );
-
- template<class VectorSpace>
- inline MatrixBlock<mType> block
- (
- const label mStart,
- const label nStart
- );
-
- inline MatrixBlock<mType> col
- (
- const label m,
- const label rowStart
- );
-
- inline MatrixBlock<mType> 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<Type> release();
@@ -307,44 +329,65 @@ public:
void swap(Matrix<Form, Type>& mat);
//- Transfer the contents of the argument Matrix into this Matrix
- //- and annul the argument Matrix.
+ //- and annul the argument Matrix
void transfer(Matrix<Form, Type>& 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<Field<Type>> Amul(const UList<Type>& x) const;
- //- Multiply matrix with vector (A * x)
+ //- Right-multiply Matrix by a column vector (A * x)
template<class Addr>
inline tmp<Field<Type>> Amul
(
const IndirectListBase<Type, Addr>& x
) const;
- //- Multiply matrix transpose with vector (AT * x, or x * A)
+ //- Left-multiply Matrix by a row vector (x * A)
inline tmp<Field<Type>> Tmul(const UList<Type>& x) const;
- //- Multiply matrix transpose with vector (AT * x, or x * A)
+ //- Left-multiply Matrix by a row vector (x * A)
template<class Addr>
inline tmp<Field<Type>> Tmul
(
const IndirectListBase<Type, Addr>& x
) const;
+ //- Extract the diagonal elements. Method may change in the future.
+ List<Type> diag() const;
+
+ //- Assign diagonal of Matrix
+ void diag(const UList<Type>& 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<Form, Type>&& mat);
- //- Assignment to a block of another matrix
+ //- Assignment to a block of another Matrix
template<class MatrixType>
void operator=(const ConstMatrixBlock<MatrixType>& Mb);
- //- Assignment to a block of another matrix
+ //- Assignment to a block of another Matrix
template<class MatrixType>
void operator=(const MatrixBlock<MatrixType>& Mb);
@@ -390,14 +433,20 @@ public:
//- Matrix subtraction
void operator-=(const Matrix<Form, Type>& 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<mType>
+ 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<mType>
+ 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<mType>
+ 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<mType>
+ 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<class Form, class Type>
@@ -462,70 +574,100 @@ template<class Form, class Type>
Ostream& operator<<(Ostream& os, const Matrix<Form, Type>& mat);
-// Global Functions, Operators
+// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * * //
-//- Find max value in the matrix
+//- Find max value in Matrix
template<class Form, class Type>
const Type& max(const Matrix<Form, Type>& mat);
-//- Find min value in the matrix
+//- Find min value in Matrix
template<class Form, class Type>
const Type& min(const Matrix<Form, Type>& mat);
-//- Find the min/max values of the matrix
+//- Find the min/max values of Matrix
template<class Form, class Type>
MinMax<Type> minMax(const Matrix<Form, Type>& mat);
+
+// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * * //
+
//- Matrix negation
template<class Form, class Type>
Form operator-(const Matrix<Form, Type>& mat);
-//- Matrix addition
-template<class Form, class Type>
-Form operator+
+//- Matrix addition. Returns Matrix of the same form as the first parameter.
+template<class Form1, class Form2, class Type>
+Form1 operator+
(
- const Matrix<Form, Type>& A,
- const Matrix<Form, Type>& B
+ const Matrix<Form1, Type>& A,
+ const Matrix<Form2, Type>& B
);
+//- Matrix subtraction. Returns Matrix of the same form as the first parameter.
+template<class Form1, class Form2, class Type>
+Form1 operator-
+(
+ const Matrix<Form1, Type>& A,
+ const Matrix<Form2, Type>& B
+);
-//- Matrix subtraction
+//- Scalar multiplication of Matrix
template<class Form, class Type>
-Form operator-
+Form operator*
(
- const Matrix<Form, Type>& A,
- const Matrix<Form, Type>& B
+ const Type& s,
+ const Matrix<Form, Type>& mat
);
-
-//- Scalar multiplication of a matrix
+//- Scalar addition of Matrix
template<class Form, class Type>
-Form operator*
+Form operator+
(
- const scalar s,
+ const Type& s,
const Matrix<Form, Type>& mat
);
+//- Scalar subtraction of Matrix
+template<class Form, class Type>
+Form operator-
+(
+ const Type& s,
+ const Matrix<Form, Type>& mat
+);
-//- Scalar multiplication of a matrix
+//- Scalar multiplication of Matrix
template<class Form, class Type>
Form operator*
(
const Matrix<Form, Type>& mat,
- const scalar s
+ const Type& s
);
+//- Scalar addition of Matrix
+template<class Form, class Type>
+Form operator+
+(
+ const Matrix<Form, Type>& mat,
+ const Type& s
+);
-//- Scalar division of a matrix
+//- Scalar subtraction of Matrix
template<class Form, class Type>
-Form operator/
+Form operator-
(
const Matrix<Form, Type>& mat,
- const scalar s
+ const Type& s
);
+//- Scalar division of Matrix
+template<class Form, class Type>
+Form operator/
+(
+ const Matrix<Form, Type>& mat,
+ const Type& s
+);
-//- Matrix-matrix multiplication
+//- Matrix-Matrix multiplication using ikj-algorithm
template<class Form1, class Form2, class Type>
typename typeOfInnerProduct<Type, Form1, Form2>::type
operator*
@@ -534,7 +676,6 @@ operator*
const Matrix<Form2, Type>& B
);
-
//- Matrix-vector multiplication (A * x), where x is a column vector
template<class Form, class Type>
inline tmp<Field<Type>> operator*
@@ -551,8 +692,7 @@ inline tmp<Field<Type>> operator*
const IndirectListBase<Type, Addr>& x
);
-
-//- Vector-matrix multiplication (x * A), where x is a row vector
+//- Vector-Matrix multiplication (x * A), where x is a row vector
template<class Form, class Type>
inline tmp<Field<Type>> operator*
(
@@ -560,7 +700,7 @@ inline tmp<Field<Type>> operator*
const Matrix<Form, Type>& mat
);
-//- Vector-matrix multiplication (x * A), where x is a row vector
+//- Vector-Matrix multiplication (x * A), where x is a row vector
template<class Form, class Type, class Addr>
inline tmp<Field<Type>> operator*
(
@@ -568,6 +708,24 @@ inline tmp<Field<Type>> operator*
const Matrix<Form, Type>& mat
);
+//- Implicit inner product of Matrix-Matrix, equivalent to A.T()*B
+template<class Form1, class Form2, class Type>
+typename typeOfInnerProduct<Type, Form1, Form2>::type
+operator&
+(
+ const Matrix<Form1, Type>& ATranspose,
+ const Matrix<Form2, Type>& B
+);
+
+//- Implicit outer product of Matrix-Matrix, equivalent to A*B.T()
+template<class Form1, class Form2, class Type>
+typename typeOfInnerProduct<Type, Form1, Form2>::type
+operator^
+(
+ const Matrix<Form1, Type>& A,
+ const Matrix<Form2, Type>& BTranspose
+);
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/OpenFOAM/matrices/Matrix/MatrixI.H b/src/OpenFOAM/matrices/Matrix/MatrixI.H
index d6e3d836c674c1f80a0c6c2ec94cf9f382033e32..4c229fbd1879636c12dbf90a3f3238c7f406021e 100644
--- a/src/OpenFOAM/matrices/Matrix/MatrixI.H
+++ b/src/OpenFOAM/matrices/Matrix/MatrixI.H
@@ -128,40 +128,36 @@ inline bool Foam::Matrix<Form, Type>::empty() const noexcept
template<class Form, class Type>
inline void Foam::Matrix<Form, Type>::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<class Form, class Type>
inline void Foam::Matrix<Form, Type>::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<class Form, class Type>
inline const Type& Foam::Matrix<Form, Type>::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<class Form, class Type>
inline Type& Foam::Matrix<Form, Type>::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,152 +262,188 @@ inline Type& Foam::Matrix<Form, Type>::at(const label idx)
template<class Form, class Type>
inline Foam::ConstMatrixBlock<Foam::Matrix<Form, Type>>
-Foam::Matrix<Form, Type>::block
+Foam::Matrix<Form, Type>::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<mType>
(
*this,
- m,
- n,
- mStart,
- nStart
+ len, // rows
+ 1,
+ rowIndex,
+ colIndex
);
}
template<class Form, class Type>
-template<class VectorSpace>
inline Foam::ConstMatrixBlock<Foam::Matrix<Form, Type>>
-Foam::Matrix<Form, Type>::block
+Foam::Matrix<Form, Type>::subRow
(
- const label mStart,
- const label nStart
+ const label rowIndex,
+ const label colIndex,
+ label len
) const
{
+ if (len < 0)
+ {
+ len = nCols_ - colIndex;
+ }
+
return ConstMatrixBlock<mType>
(
*this,
- VectorSpace::mRows,
- VectorSpace::nCols,
- mStart,
- nStart
+ 1,
+ len, // columns
+ rowIndex,
+ colIndex
);
}
template<class Form, class Type>
inline Foam::ConstMatrixBlock<Foam::Matrix<Form, Type>>
-Foam::Matrix<Form, Type>::col
+Foam::Matrix<Form, Type>::subMatrix
(
- const label m,
- const label mStart
+ 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<mType>
(
*this,
- m,
- 1,
- mStart,
- 0
+ szRows,
+ szCols,
+ rowIndex,
+ colIndex
);
}
template<class Form, class Type>
+template<class VectorSpace>
inline Foam::ConstMatrixBlock<Foam::Matrix<Form, Type>>
-Foam::Matrix<Form, Type>::col
+Foam::Matrix<Form, Type>::block
(
- const label m,
- const label mStart,
- const label nStart
+ const label rowIndex,
+ const label colIndex
) const
{
return ConstMatrixBlock<mType>
(
*this,
- m,
- 1,
- mStart,
- nStart
+ VectorSpace::mRows,
+ VectorSpace::nCols,
+ rowIndex,
+ colIndex
);
}
template<class Form, class Type>
inline Foam::MatrixBlock<Foam::Matrix<Form, Type>>
-Foam::Matrix<Form, Type>::block
+Foam::Matrix<Form, Type>::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<mType>
(
*this,
- m,
- n,
- mStart,
- nStart
+ len, // rows
+ 1,
+ rowIndex,
+ colIndex
);
}
template<class Form, class Type>
-template<class VectorSpace>
inline Foam::MatrixBlock<Foam::Matrix<Form, Type>>
-Foam::Matrix<Form, Type>::block(const label mStart, const label nStart)
+Foam::Matrix<Form, Type>::subRow
+(
+ const label rowIndex,
+ const label colIndex,
+ label len
+)
{
+ if (len < 0)
+ {
+ len = nCols_ - colIndex;
+ }
+
return MatrixBlock<mType>
(
*this,
- VectorSpace::mRows,
- VectorSpace::nCols,
- mStart,
- nStart
+ 1,
+ len, // columns
+ rowIndex,
+ colIndex
);
}
template<class Form, class Type>
inline Foam::MatrixBlock<Foam::Matrix<Form, Type>>
-Foam::Matrix<Form, Type>::col(const label m, const label mStart)
+Foam::Matrix<Form, Type>::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<mType>
(
*this,
- m,
- 1,
- mStart,
- 0
+ szRows,
+ szCols,
+ rowIndex,
+ colIndex
);
}
template<class Form, class Type>
+template<class VectorSpace>
inline Foam::MatrixBlock<Foam::Matrix<Form, Type>>
-Foam::Matrix<Form, Type>::col
+Foam::Matrix<Form, Type>::block
(
- const label m,
- const label mStart,
- const label nStart
+ const label rowIndex,
+ const label colIndex
)
{
return MatrixBlock<mType>
(
*this,
- m,
- 1,
- mStart,
- nStart
+ VectorSpace::mRows,
+ VectorSpace::nCols,
+ rowIndex,
+ colIndex
);
}
@@ -532,8 +564,10 @@ inline const Type& Foam::Matrix<Form, Type>::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<Form, Type>::operator()
const label jcol
)
{
+ #ifdef FULLDEBUG
checki(irow);
checkj(jcol);
+ #endif
return v_[irow*nCols_ + jcol];
}
@@ -554,7 +590,9 @@ inline Type& Foam::Matrix<Form, Type>::operator()
template<class Form, class Type>
inline const Type* Foam::Matrix<Form, Type>::operator[](const label irow) const
{
+ #ifdef FULLDEBUG
checki(irow);
+ #endif
return v_ + irow*nCols_;
}
@@ -562,7 +600,9 @@ inline const Type* Foam::Matrix<Form, Type>::operator[](const label irow) const
template<class Form, class Type>
inline Type* Foam::Matrix<Form, Type>::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 0000000000000000000000000000000000000000..5ce736cdf3980a003a81c8bf5635474d125953d2
--- /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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "MatrixTools.H"
+
+// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * * //
+
+template<class Form1, class Form2, class Type>
+bool Foam::MatrixTools::equal
+(
+ const Matrix<Form1, Type>& A,
+ const Matrix<Form2, Type>& 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<class Container>
+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 0000000000000000000000000000000000000000..8c63efc0fe0d2ba349caefb4e73d26fb350ac97e
--- /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 <http://www.gnu.org/licenses/>.
+
+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<class Form1, class Form2, class Type>
+bool equal
+(
+ const Matrix<Form1, Type>& A,
+ const Matrix<Form2, Type>& B,
+ const bool verbose = false,
+ const scalar relTol = 1e-5,
+ const scalar absTol = 1e-8
+);
+
+
+//- Simple ASCII output of Matrix, MatrixBlock
+template<class Container>
+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 2784a342a369fec7379c0691da68aa302285eb3c..7ced3ca8de125530acbffb654e374810de26d053 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<MatrixType>::operator Field<cmptType>() const
}
+template<class MatrixType>
+Foam::label Foam::ConstMatrixBlock<MatrixType>::disallow
+(
+ const char* what
+) const
+{
+ FatalErrorInFunction
+ << "Block addresses " << what
+ << " outside matrix or invalid matrix components"
+ << abort(FatalError);
+ return 0;
+}
+
+
+template<class MatrixType>
+Foam::label Foam::MatrixBlock<MatrixType>::disallow
+(
+ const char* what
+) const
+{
+ FatalErrorInFunction
+ << "Block addresses " << what
+ << " outside matrix or invalid matrix components"
+ << abort(FatalError);
+ return 0;
+}
+
+
+template<class MatrixType> void Foam::ConstMatrixBlock<MatrixType>::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<class MatrixType> void Foam::MatrixBlock<MatrixType>::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<class MatrixType>
diff --git a/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.H b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.H
index 27e90c7a9f4d77811d82e8fcd720a17201734e86..bb3cf1d95e9194888cbda316528628445c7e7867 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<cmptType>() 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<cmptType>() 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<class Form>
diff --git a/src/OpenFOAM/matrices/MatrixBlock/MatrixBlockI.H b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlockI.H
index 408f1393a10213472dbbc78186f4e129f6b59829..e6a47e0ba4c63e6c7f363750a7bb1ce689d8b959 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<MatrixType>::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_
+ (
+ 0 <= mStart
+ || mStart + mRows_ <= matrix.m()
+ ? mStart : disallow("row start")
+ ),
+ colStart_
(
- rowStart_ + mRows_ > matrix.m()
- || colStart_ + nCols_ > matrix.n()
+ 0 <= nStart
+ || nStart + nCols_ <= matrix.n()
+ ? nStart : disallow("col start")
)
- {
- FatalErrorInFunction
- << "Block addresses outside matrix"
- << abort(FatalError);
- }
- #endif
-}
+{}
template<class MatrixType>
@@ -69,24 +66,21 @@ Foam::MatrixBlock<MatrixType>::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_
+ (
+ 0 <= mStart
+ || mStart + mRows_ <= matrix.m()
+ ? mStart : disallow("row start")
+ ),
+ colStart_
(
- rowStart_ + mRows_ > matrix.m()
- || colStart_ + nCols_ > matrix.n()
+ 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<MatrixType>::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<MatrixType>::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<MatrixType>::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 ea3c437a860cbedd2f96289777dee65ed188d8e3..624db7b7dcd0dc6e87378bb16c940bee028e8f18 100644
--- a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
+++ b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
@@ -282,7 +282,7 @@ Foam::QRMatrix<MatrixType>::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 68ecd56a598cfadac6e473c1fe3fc2949ff033d3..668e0241b013dbd326acb5444144bbbe184b60b1 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 Type>
class typeOfInnerProduct<Type, RectangularMatrix<Type>, RectangularMatrix<Type>>
@@ -128,6 +128,7 @@ public:
typedef RectangularMatrix<Type> type;
};
+
template<class Type>
class typeOfInnerProduct<Type, RectangularMatrix<Type>, SquareMatrix<Type>>
{
@@ -136,6 +137,7 @@ public:
typedef RectangularMatrix<Type> type;
};
+
template<class Type>
class typeOfInnerProduct<Type, SquareMatrix<Type>, RectangularMatrix<Type>>
{
@@ -145,6 +147,8 @@ public:
};
+// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
+
template<class Type>
RectangularMatrix<Type> outer(const Field<Type>& f1, const Field<Type>& f2);
diff --git a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H
index da3e03b53f3025f61cad5ab1b2313baf84e4c00e..1ff31473717866cc5d894e7a972226fbe59e6751 100644
--- a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H
+++ b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H
@@ -190,9 +190,9 @@ inline Foam::RectangularMatrix<Type> outer
{
RectangularMatrix<Type> f1f2T(f1.size(), f2.size());
- for (label i=0; i<f1f2T.m(); i++)
+ for (label i = 0; i < f1f2T.m(); ++i)
{
- for (label j=0; j<f1f2T.n(); j++)
+ for (label j = 0; j < f1f2T.n(); ++j)
{
f1f2T(i, j) = f1[i]*f2[j];
}
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C
index b48d01545c376a421fde9fb38c1d4d9f4afa925f..848b4901b32ad2d6b51b17b56bdfd7abb6755b1b 100644
--- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C
+++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.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) 2013-2016 OpenFOAM Foundation
@@ -26,6 +26,7 @@ License
\*---------------------------------------------------------------------------*/
#include "SquareMatrix.H"
+#include "RectangularMatrix.H"
#include "labelList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -39,7 +40,7 @@ Foam::scalar Foam::detDecomposed
{
Type diagProduct = pTraits<Type>::one;
- for (label i=0; i<matrix.m(); i++)
+ for (label i = 0; i < matrix.m(); ++i)
{
diagProduct *= matrix(i, i);
}
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
index dec80275bd28fdb80584b8ccb1e4b66c7d404cb6..1bd82812157b4d15bd8ef9edd38932e1436a9505 100644
--- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
+++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
@@ -149,7 +149,7 @@ public:
};
-// Global functions and operators
+// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
//- Return the LU decomposed SquareMatrix det
template<class Type>
@@ -171,7 +171,6 @@ public:
typedef SquareMatrix<Type> type;
};
-
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
index 4e1b3ca483b02b06388b1a4bbbf97cf69db86098..47897b61dc033b4df7751e52797874e06719cb97 100644
--- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
+++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
@@ -251,9 +251,9 @@ inline Foam::SquareMatrix<Type> symmOuter
{
SquareMatrix<Type> f1f2T(f1.size());
- for (label i=0; i<f1f2T.m(); i++)
+ for (label i = 0; i < f1f2T.m(); ++i)
{
- for (label j=0; j<f1f2T.n(); j++)
+ for (label j = 0; j < f1f2T.n(); ++j)
{
f1f2T(i, j) = f1[i]*f2[j];
}