MatrixI.H
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2016 OpenFOAM Foundation
9  Copyright (C) 2019-2021 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "MatrixBlock.H"
30 
31 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
32 
33 template<class Form, class Type>
35 {
36  const label len = size();
37 
38  if (len > 0)
39  {
40  // With sign-check to avoid spurious -Walloc-size-larger-than
41  v_ = new Type[len];
42  }
43 }
44 
45 
46 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
47 
48 template<class Form, class Type>
50 :
51  mRows_(0),
52  nCols_(0),
53  v_(nullptr)
54 {}
55 
56 
57 template<class Form, class Type>
59 :
60  Matrix<Form, Type>(dims.first(), dims.second())
61 {}
62 
63 
64 template<class Form, class Type>
66 :
67  Matrix<Form, Type>(dims.first(), dims.second(), Foam::zero{})
68 {}
69 
70 
71 template<class Form, class Type>
72 inline Foam::Matrix<Form, Type>::Matrix(const labelPair& dims, const Type& val)
73 :
74  Matrix<Form, Type>(dims.first(), dims.second(), val)
75 {}
76 
77 
78 template<class Form, class Type>
81 {
83 }
84 
85 
86 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
87 
88 template<class Form, class Type>
90 {
91  return NullObjectRef<Matrix<Form, Type>>();
92 }
93 
94 
95 template<class Form, class Type>
96 inline Foam::label Foam::Matrix<Form, Type>::size() const
97 {
98  return mRows_ * nCols_;
99 }
100 
101 
102 template<class Form, class Type>
104 {
105  return labelPair(mRows_, nCols_);
106 }
107 
108 
109 template<class Form, class Type>
111 {
112  return !mRows_ || !nCols_;
113 }
114 
115 
116 template<class Form, class Type>
117 inline void Foam::Matrix<Form, Type>::checki(const label i) const
118 {
119  if (!mRows_ || !nCols_)
120  {
122  << "Attempt to access element from empty matrix"
123  << abort(FatalError);
124  }
125  if (i < 0 || mRows_ <= i)
126  {
128  << "Index " << i << " out of range 0 ... " << mRows_-1
129  << abort(FatalError);
130  }
131 }
132 
133 
134 template<class Form, class Type>
135 inline void Foam::Matrix<Form, Type>::checkj(const label j) const
136 {
137  if (!mRows_ || !nCols_)
138  {
140  << "Attempt to access element from empty matrix"
141  << abort(FatalError);
142  }
143  if (j < 0 || nCols_ <= j)
144  {
146  << "index " << j << " out of range 0 ... " << nCols_-1
147  << abort(FatalError);
148  }
149 }
150 
151 
152 template<class Form, class Type>
153 inline void Foam::Matrix<Form, Type>::checkSize() const
154 {
155  if (mRows_ < 0 || nCols_ < 0)
156  {
158  << "Incorrect size (" << mRows_ << ", " << nCols_ << ')' << nl
160  }
161  // Could also check for odd sizes, like (0, N) and make (0,0)
162 }
163 
164 
165 template<class Form, class Type>
166 inline bool Foam::Matrix<Form, Type>::uniform() const
167 {
168  const label len = size();
169 
170  if (!len)
171  {
172  return false;
173  }
174 
175  // std::all_of()
176  for (label idx = 1; idx < len; ++idx)
177  {
178  if (v_[0] != v_[idx])
179  {
180  return false;
181  }
182  }
183 
184  return true;
185 }
186 
187 
188 template<class Form, class Type>
189 inline const Type* Foam::Matrix<Form, Type>::cdata() const noexcept
190 {
191  return v_;
192 }
193 
194 
195 template<class Form, class Type>
197 {
198  return v_;
199 }
200 
201 
202 template<class Form, class Type>
204 {
205  return reinterpret_cast<const char*>(v_);
206 }
207 
208 
209 template<class Form, class Type>
211 {
212  return reinterpret_cast<char*>(v_);
213 }
214 
215 
216 template<class Form, class Type>
217 inline std::streamsize Foam::Matrix<Form, Type>::size_bytes() const noexcept
218 {
219  return mRows_*nCols_*sizeof(Type);
220 }
221 
222 
223 template<class Form, class Type>
224 inline const Type* Foam::Matrix<Form, Type>::rowData(const label irow) const
225 {
226  #ifdef FULLDEBUG
227  checki(irow);
228  #endif
229  return (v_ + irow*nCols_);
230 }
231 
232 
233 template<class Form, class Type>
234 inline Type* Foam::Matrix<Form, Type>::rowData(const label irow)
235 {
236  #ifdef FULLDEBUG
237  checki(irow);
238  #endif
239  return (v_ + irow*nCols_);
240 }
241 
242 
243 template<class Form, class Type>
244 inline const Type& Foam::Matrix<Form, Type>::at(const label idx) const
245 {
246  #ifdef FULLDEBUG
247  if (idx < 0 || this->size() <= idx)
248  {
250  << "index " << idx << " out of range 0 ... " << this->size()
251  << abort(FatalError);
252  }
253  #endif
254  return *(v_ + idx);
255 }
256 
257 
258 template<class Form, class Type>
259 inline Type& Foam::Matrix<Form, Type>::at(const label idx)
260 {
261  #ifdef FULLDEBUG
262  if (idx < 0 || this->size() <= idx)
263  {
265  << "index " << idx << " out of range 0 ... " << this->size()
266  << abort(FatalError);
267  }
268  #endif
269  return *(v_ + idx);
270 }
271 
272 
273 template<class Form, class Type>
276 (
277  const label colIndex,
278  const label rowIndex,
279  label len
280 ) const
281 {
282  if (len < 0)
283  {
284  len = mRows_ - rowIndex;
285  }
286 
287  return ConstMatrixBlock<mType>
288  (
289  *this,
290  len, // rows
291  1,
292  rowIndex,
293  colIndex
294  );
295 }
296 
297 
298 template<class Form, class Type>
301 (
302  const label rowIndex,
303  const label colIndex,
304  label len
305 ) const
306 {
307  if (len < 0)
308  {
309  len = nCols_ - colIndex;
310  }
311 
312  return ConstMatrixBlock<mType>
313  (
314  *this,
315  1,
316  len, // columns
317  rowIndex,
318  colIndex
319  );
320 }
321 
322 
323 template<class Form, class Type>
326 (
327  const label rowIndex,
328  const label colIndex,
329  label szRows,
330  label szCols
331 ) const
332 {
333  if (szRows < 0) szRows = mRows_ - rowIndex;
334  if (szCols < 0) szCols = nCols_ - colIndex;
335 
337  (
338  *this,
339  szRows,
340  szCols,
341  rowIndex,
342  colIndex
343  );
344 }
345 
346 
347 template<class Form, class Type>
348 template<class VectorSpace>
351 (
352  const label rowIndex,
353  const label colIndex
354 ) const
355 {
357  (
358  *this,
359  VectorSpace::mRows,
360  VectorSpace::nCols,
361  rowIndex,
362  colIndex
363  );
364 }
365 
366 
367 template<class Form, class Type>
370 (
371  const label colIndex,
372  const label rowIndex,
373  label len
374 )
375 {
376  if (len < 0)
377  {
378  len = mRows_ - rowIndex;
379  }
380 
381  return MatrixBlock<mType>
382  (
383  *this,
384  len, // rows
385  1,
386  rowIndex,
387  colIndex
388  );
389 }
390 
391 
392 template<class Form, class Type>
395 (
396  const label rowIndex,
397  const label colIndex,
398  label len
399 )
400 {
401  if (len < 0)
402  {
403  len = nCols_ - colIndex;
404  }
405 
406  return MatrixBlock<mType>
407  (
408  *this,
409  1,
410  len, // columns
411  rowIndex,
412  colIndex
413  );
414 }
415 
416 
417 template<class Form, class Type>
420 (
421  const label rowIndex,
422  const label colIndex,
423  label szRows,
424  label szCols
425 )
426 {
427  if (szRows < 0) szRows = mRows_ - rowIndex;
428  if (szCols < 0) szCols = nCols_ - colIndex;
429 
430  return MatrixBlock<mType>
431  (
432  *this,
433  szRows,
434  szCols,
435  rowIndex,
436  colIndex
437  );
438 }
439 
440 
441 template<class Form, class Type>
442 template<class VectorSpace>
445 (
446  const label rowIndex,
447  const label colIndex
448 )
449 {
450  return MatrixBlock<mType>
451  (
452  *this,
453  VectorSpace::mRows,
454  VectorSpace::nCols,
455  rowIndex,
456  colIndex
457  );
458 }
459 
460 
461 template<class Form, class Type>
462 inline void Foam::Matrix<Form, Type>::setSize(const label m, const label n)
463 {
464  resize(m, n);
465 }
466 
467 
468 template<class Form, class Type>
469 void Foam::Matrix<Form, Type>::shallowResize(const label m, const label n)
470 {
471  mRows_ = m;
472  nCols_ = n;
473 }
474 
475 
476 template<class Form, class Type>
478 (
479  const UList<Type>& x
480 ) const
481 {
482  return this->AmulImpl(x);
483 }
484 
485 
486 template<class Form, class Type>
487 template<class Addr>
489 (
491 ) const
492 {
493  return this->AmulImpl(x);
494 }
495 
496 
497 template<class Form, class Type>
499 (
500  const UList<Type>& x
501 ) const
502 {
503  return this->TmulImpl(x);
504 }
505 
506 
507 template<class Form, class Type>
508 template<class Addr>
510 (
512 ) const
513 {
514  return this->TmulImpl(x);
515 }
516 
517 
518 // * * * * * * * * * * * * * * * * * Iterators * * * * * * * * * * * * * * * //
519 
520 template<class Form, class Type>
523 {
524  return v_;
525 }
526 
527 
528 template<class Form, class Type>
531 {
532  return v_ + (mRows_ * nCols_);
533 }
534 
535 
536 template<class Form, class Type>
539 {
540  return v_;
541 }
542 
543 
544 template<class Form, class Type>
547 {
548  return v_ + (mRows_ * nCols_);
549 }
550 
551 
552 template<class Form, class Type>
555 {
556  return v_;
557 }
558 
559 
560 template<class Form, class Type>
563 {
564  return v_ + (mRows_ * nCols_);
565 }
566 
567 
568 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
569 
570 template<class Form, class Type>
571 inline const Type& Foam::Matrix<Form, Type>::operator()
572 (
573  const label irow,
574  const label jcol
575 ) const
576 {
577  #ifdef FULLDEBUG
578  checki(irow);
579  checkj(jcol);
580  #endif
581  return v_[irow*nCols_ + jcol];
582 }
583 
584 
585 template<class Form, class Type>
587 (
588  const label irow,
589  const label jcol
590 )
591 {
592  #ifdef FULLDEBUG
593  checki(irow);
594  checkj(jcol);
595  #endif
596  return v_[irow*nCols_ + jcol];
597 }
598 
599 
600 template<class Form, class Type>
601 inline const Type* Foam::Matrix<Form, Type>::operator[](const label irow) const
602 {
603  #ifdef FULLDEBUG
604  checki(irow);
605  #endif
606  return v_ + irow*nCols_;
607 }
608 
609 
610 template<class Form, class Type>
611 inline Type* Foam::Matrix<Form, Type>::operator[](const label irow)
612 {
613  #ifdef FULLDEBUG
614  checki(irow);
615  #endif
616  return v_ + irow*nCols_;
617 }
618 
619 
620 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
621 
622 namespace Foam
623 {
624 
625 // * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * //
626 
627 //- Matrix-vector multiplication (A * x), where x is a column vector
628 template<class Form, class Type>
629 inline tmp<Field<Type>> operator*
630 (
631  const Matrix<Form, Type>& mat,
632  const UList<Type>& x
633 )
634 {
635  return mat.Amul(x);
636 }
637 
638 
639 //- Matrix-vector multiplication (A * x), where x is a column vector
640 template<class Form, class Type, class Addr>
641 inline tmp<Field<Type>> operator*
642 (
643  const Matrix<Form, Type>& mat,
645 )
646 {
647  return mat.Amul(x);
648 }
649 
650 
651 //- Vector-Matrix multiplication (x * A), where x is a row vector
652 template<class Form, class Type>
653 inline tmp<Field<Type>> operator*
654 (
655  const UList<Type>& x,
656  const Matrix<Form, Type>& mat
657 )
658 {
659  return mat.Tmul(x);
660 }
661 
662 
663 //- Vector-Matrix multiplication (x * A), where x is a row vector
664 template<class Form, class Type, class Addr>
665 inline tmp<Field<Type>> operator*
666 (
668  const Matrix<Form, Type>& mat
669 )
670 {
671  return mat.Tmul(x);
672 }
673 
674 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
675 
676 } // End namespace Foam
677 
678 // ************************************************************************* //
void shallowResize(const label m, const label n)
Resize Matrix without reallocating storage (unsafe)
Definition: MatrixI.H:462
ConstMatrixBlock< mType > subColumn(const label colIndex, const label rowIndex=0, label len=-1) const
Return const column or column&#39;s subset of Matrix.
Definition: MatrixI.H:269
ConstMatrixBlock< mType > subRow(const label rowIndex, const label colIndex=0, label len=-1) const
Return const row or const row&#39;s subset of Matrix.
Definition: MatrixI.H:294
autoPtr< mType > clone() const
Clone.
Definition: MatrixI.H:73
patchWriters resize(patchIds.size())
tmp< Field< Type > > Tmul(const UList< Type > &x) const
Left-multiply Matrix by a row vector (x * A)
Definition: MatrixI.H:492
error FatalError
Error stream (stdout output on all processes), with additional &#39;FOAM FATAL ERROR&#39; header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:598
constexpr char nl
The newline &#39;\n&#39; character (0x0a)
Definition: Ostream.H:50
tmp< Field< Type > > Amul(const UList< Type > &x) const
Right-multiply Matrix by a column vector (A * x)
Definition: MatrixI.H:471
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tf1, const word &name, const dimensionSet &dimensions, const bool initCopy=false)
Global function forwards to reuseTmpDimensionedField::New.
A templated block of an (m x n) matrix of type <MatrixType>.
Definition: Matrix.H:62
void setSize(const label m, const label n)
Change Matrix dimensions, preserving the elements.
Definition: MatrixI.H:455
const Type * operator[](const label irow) const
Return const pointer to data in the specified row - rowData().
Definition: MatrixI.H:594
iterator begin() noexcept
Return an iterator to begin traversing a Matrix.
Definition: MatrixI.H:515
Base for lists with indirect addressing, templated on the list contents type and the addressing type...
bool empty() const noexcept
Return true if Matrix is empty (i.e., size() is zero)
Definition: MatrixI.H:103
ConstMatrixBlock< mType > subMatrix(const label rowIndex, const label colIndex, label szRows=-1, label szCols=-1) const
Return const sub-block of Matrix.
Definition: MatrixI.H:319
char * data_bytes() noexcept
Return pointer to the underlying array serving as data storage, reinterpreted as byte data...
Definition: MatrixI.H:203
static const Matrix< Form, Type > & null()
Return a null Matrix.
Definition: MatrixI.H:82
label size() const
The number of elements in Matrix (m*n)
Definition: MatrixI.H:89
const Type & at(const label idx) const
Linear addressing const element access.
Definition: MatrixI.H:237
Type * iterator
Random access iterator for traversing a Matrix.
Definition: Matrix.H:137
labelPair sizes() const
Return row/column sizes.
Definition: MatrixI.H:96
bool uniform() const
True if all entries have identical values, and Matrix is non-empty.
Definition: MatrixI.H:159
iterator end() noexcept
Return an iterator to end traversing a Matrix.
Definition: MatrixI.H:523
errorManip< error > abort(error &err)
Definition: errorManip.H:139
const Type * cdata() const noexcept
Return const pointer to the first data element, which can also be used to address into Matrix content...
Definition: MatrixI.H:182
void checki(const label irow) const
Check index i is within valid range [0, m)
Definition: MatrixI.H:110
A templated (m x n) matrix of objects of <T>. The layout is (mRows x nCols) - row-major order: ...
const direction noexcept
Definition: Scalar.H:258
std::streamsize size_bytes() const noexcept
Number of contiguous bytes for the Matrix data, no runtime check that the type is actually contiguous...
Definition: MatrixI.H:210
Pair< label > labelPair
A pair of labels.
Definition: Pair.H:51
const char * cdata_bytes() const noexcept
Return pointer to the underlying array serving as data storage, reinterpreted as byte data...
Definition: MatrixI.H:196
const Type * const_iterator
Random access iterator for traversing a Matrix.
Definition: Matrix.H:142
const_iterator cend() const noexcept
Return const_iterator to end traversing a constant Matrix.
Definition: MatrixI.H:539
const Type * rowData(const label irow) const
Return const pointer to data in the specified row.
Definition: MatrixI.H:217
void checkj(const label jcol) const
Check index j is within valid range [0, n)
Definition: MatrixI.H:128
Matrix() noexcept
Default construct (empty matrix)
Definition: MatrixI.H:42
Type * data() noexcept
Return pointer to the first data element, which can also be used to address into Matrix contents...
Definition: MatrixI.H:189
A class representing the concept of 0 (zero) that can be used to avoid manipulating objects known to ...
Definition: zero.H:57
label n
Pointer management similar to std::unique_ptr, with some additional methods and type checking...
Definition: HashPtrTable.H:48
const_iterator cbegin() const noexcept
Return const_iterator to begin traversing a constant Matrix.
Definition: MatrixI.H:531
A class for managing temporary objects.
Definition: HashPtrTable.H:50
ConstMatrixBlock< mType > block(const label rowIndex, const label colIndex) const
Access Field as a ConstMatrixBlock.
Namespace for OpenFOAM.
void checkSize() const
Check that dimensions are positive, non-zero.
Definition: MatrixI.H:146