Matrix.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "Matrix.H"
27 
28 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
29 
30 template<class Form, class Type>
32 {
33  if (mRows_ && nCols_)
34  {
35  v_ = new Type[size()];
36  }
37 }
38 
39 
40 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
41 
42 template<class Form, class Type>
44 :
45  mRows_(m),
46  nCols_(n),
47  v_(nullptr)
48 {
49  if (mRows_ < 0 || nCols_ < 0)
50  {
52  << "Incorrect m, n " << mRows_ << ", " << nCols_
53  << abort(FatalError);
54  }
55 
56  allocate();
57 }
58 
59 
60 template<class Form, class Type>
62 :
63  mRows_(m),
64  nCols_(n),
65  v_(nullptr)
66 {
67  if (mRows_ < 0 || nCols_ < 0)
68  {
70  << "Incorrect m, n " << mRows_ << ", " << nCols_
71  << abort(FatalError);
72  }
73 
74  allocate();
75 
76  if (v_)
77  {
78  const label mn = size();
79  for (label i=0; i<mn; i++)
80  {
81  v_[i] = Zero;
82  }
83  }
84 }
85 
86 
87 template<class Form, class Type>
88 Foam::Matrix<Form, Type>::Matrix(const label m, const label n, const Type& s)
89 :
90  mRows_(m),
91  nCols_(n),
92  v_(nullptr)
93 {
94  if (mRows_ < 0 || nCols_ < 0)
95  {
97  << "Incorrect m, n " << mRows_ << ", " << nCols_
98  << abort(FatalError);
99  }
100 
101  allocate();
102 
103  if (v_)
104  {
105  const label mn = size();
106  for (label i=0; i<mn; i++)
107  {
108  v_[i] = s;
109  }
110  }
111 }
112 
113 
114 template<class Form, class Type>
116 :
117  mRows_(M.mRows_),
118  nCols_(M.nCols_),
119  v_(nullptr)
120 {
121  if (M.v_)
122  {
123  allocate();
124 
125  const label mn = size();
126  for (label i=0; i<mn; i++)
127  {
128  v_[i] = M.v_[i];
129  }
130  }
131 }
132 
133 
134 template<class Form, class Type>
135 template<class Form2>
137 :
138  mRows_(M.m()),
139  nCols_(M.n()),
140  v_(nullptr)
141 {
142  if (M.v())
143  {
144  allocate();
145 
146  const label mn = size();
147  for (label i=0; i<mn; i++)
148  {
149  v_[i] = M.v()[i];
150  }
151  }
152 }
153 
154 
155 template<class Form, class Type>
156 template<class MatrixType>
158 (
160 )
161 :
162  mRows_(Mb.m()),
163  nCols_(Mb.n())
164 {
165  allocate();
166 
167  for (label i=0; i<mRows_; i++)
168  {
169  for (label j=0; j<nCols_; j++)
170  {
171  (*this)(i,j) = Mb(i,j);
172  }
173  }
174 }
175 
176 
177 template<class Form, class Type>
178 template<class MatrixType>
180 (
181  const MatrixBlock<MatrixType>& Mb
182 )
183 :
184  mRows_(Mb.m()),
185  nCols_(Mb.n())
186 {
187  allocate();
188 
189  for (label i=0; i<mRows_; i++)
190  {
191  for (label j=0; j<nCols_; j++)
192  {
193  (*this)(i,j) = Mb(i,j);
194  }
195  }
196 }
197 
198 
199 // * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * * //
200 
201 template<class Form, class Type>
203 {
204  if (v_)
205  {
206  delete[] v_;
207  }
208 }
209 
210 
211 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
212 
213 template<class Form, class Type>
215 {
216  if (v_)
217  {
218  delete[] v_;
219  v_ = nullptr;
220  }
221 
222  mRows_ = 0;
223  nCols_ = 0;
224 }
225 
226 
227 template<class Form, class Type>
229 {
230  clear();
231 
232  mRows_ = M.mRows_;
233  M.mRows_ = 0;
234 
235  nCols_ = M.nCols_;
236  M.nCols_ = 0;
237 
238  v_ = M.v_;
239  M.v_ = nullptr;
240 }
241 
242 
243 template<class Form, class Type>
245 {
246  mType newMatrix(m, n, Zero);
247 
248  label minM = min(m, mRows_);
249  label minN = min(n, nCols_);
250 
251  for (label i=0; i<minM; i++)
252  {
253  for (label j=0; j<minN; j++)
254  {
255  newMatrix(i, j) = (*this)(i, j);
256  }
257  }
258 
259  transfer(newMatrix);
260 }
261 
262 
263 template<class Form, class Type>
265 {
266  const Matrix<Form, Type>& A = *this;
267  Form At(n(), m());
268 
269  for (label i=0; i<m(); i++)
270  {
271  for (label j=0; j<n(); j++)
272  {
273  At(j, i) = A(i, j);
274  }
275  }
276 
277  return At;
278 }
279 
280 
281 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
282 
283 template<class Form, class Type>
285 {
286  if (this == &M)
287  {
289  << "Attempted assignment to self"
290  << abort(FatalError);
291  }
292 
293  if (mRows_ != M.mRows_ || nCols_ != M.nCols_)
294  {
295  clear();
296  mRows_ = M.mRows_;
297  nCols_ = M.nCols_;
298  allocate();
299  }
300 
301  if (v_)
302  {
303  const label mn = size();
304  for (label i=0; i<mn; i++)
305  {
306  v_[i] = M.v_[i];
307  }
308  }
309 }
310 
311 
312 template<class Form, class Type>
313 template<class MatrixType>
314 void Foam::Matrix<Form, Type>::operator=
315 (
317 )
318 {
319  for (label i=0; i<mRows_; i++)
320  {
321  for (label j=0; j<nCols_; j++)
322  {
323  (*this)(i,j) = Mb(i,j);
324  }
325  }
326 }
327 
328 
329 template<class Form, class Type>
330 template<class MatrixType>
331 void Foam::Matrix<Form, Type>::operator=
332 (
333  const MatrixBlock<MatrixType>& Mb
334 )
335 {
336  for (label i=0; i<mRows_; i++)
337  {
338  for (label j=0; j<nCols_; j++)
339  {
340  (*this)(i,j) = Mb(i,j);
341  }
342  }
343 }
344 
345 
346 template<class Form, class Type>
348 {
349  if (v_)
350  {
351  const label mn = size();
352  for (label i=0; i<mn; i++)
353  {
354  v_[i] = s;
355  }
356  }
357 }
358 
359 
360 template<class Form, class Type>
362 {
363  if (v_)
364  {
365  const label mn = size();
366  for (label i=0; i<mn; i++)
367  {
368  v_[i] = Zero;
369  }
370  }
371 }
372 
373 
374 // * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
375 
376 template<class Form, class Type>
377 const Type& Foam::max(const Matrix<Form, Type>& M)
378 {
379  const label mn = M.size();
380 
381  if (mn)
382  {
383  label curMaxI = 0;
384  const Type* Mv = M.v();
385 
386  for (label i=1; i<mn; i++)
387  {
388  if (Mv[i] > Mv[curMaxI])
389  {
390  curMaxI = i;
391  }
392  }
393 
394  return Mv[curMaxI];
395  }
396  else
397  {
399  << "Matrix is empty"
400  << abort(FatalError);
401 
402  // Return in error to keep compiler happy
403  return M(0, 0);
404  }
405 }
406 
407 
408 template<class Form, class Type>
409 const Type& Foam::min(const Matrix<Form, Type>& M)
410 {
411  const label mn = M.size();
412 
413  if (mn)
414  {
415  label curMinI = 0;
416  const Type* Mv = M.v();
417 
418  for (label i=1; i<mn; i++)
419  {
420  if (Mv[i] < Mv[curMinI])
421  {
422  curMinI = i;
423  }
424  }
425 
426  return Mv[curMinI];
427  }
428  else
429  {
431  << "Matrix is empty"
432  << abort(FatalError);
433 
434  // Return in error to keep compiler happy
435  return M(0, 0);
436  }
437 }
438 
439 
440 // * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * //
441 
442 template<class Form, class Type>
444 {
445  Form nM(M.m(), M.n());
446 
447  if (M.m() && M.n())
448  {
449  Type* nMv = nM.v();
450  const Type* Mv = M.v();
451 
452  const label mn = M.size();
453  for (label i=0; i<mn; i++)
454  {
455  nMv[i] = -Mv[i];
456  }
457  }
458 
459  return nM;
460 }
461 
462 
463 template<class Form, class Type>
465 {
466  if (A.m() != B.m())
467  {
469  << "Attempt to add matrices with different numbers of rows: "
470  << A.m() << ", " << B.m()
471  << abort(FatalError);
472  }
473 
474  if (A.n() != B.n())
475  {
477  << "Attempt to add matrices with different numbers of columns: "
478  << A.n() << ", " << B.n()
479  << abort(FatalError);
480  }
481 
482  Form AB(A.m(), A.n());
483 
484  Type* ABv = AB.v();
485  const Type* Av = A.v();
486  const Type* Bv = B.v();
487 
488  const label mn = A.size();
489  for (label i=0; i<mn; i++)
490  {
491  ABv[i] = Av[i] + Bv[i];
492  }
493 
494  return AB;
495 }
496 
497 
498 template<class Form, class Type>
500 {
501  if (A.m() != B.m())
502  {
504  << "Attempt to add matrices with different numbers of rows: "
505  << A.m() << ", " << B.m()
506  << abort(FatalError);
507  }
508 
509  if (A.n() != B.n())
510  {
512  << "Attempt to add matrices with different numbers of columns: "
513  << A.n() << ", " << B.n()
514  << abort(FatalError);
515  }
516 
517  Form AB(A.m(), A.n());
518 
519  Type* ABv = AB.v();
520  const Type* Av = A.v();
521  const Type* Bv = B.v();
522 
523  const label mn = A.size();
524  for (label i=0; i<mn; i++)
525  {
526  ABv[i] = Av[i] - Bv[i];
527  }
528 
529  return AB;
530 }
531 
532 
533 template<class Form, class Type>
534 Form Foam::operator*(const scalar s, const Matrix<Form, Type>& M)
535 {
536  Form sM(M.m(), M.n());
537 
538  if (M.m() && M.n())
539  {
540  Type* sMv = sM.v();
541  const Type* Mv = M.v();
542 
543  const label mn = M.size();
544  for (label i=0; i<mn; i++)
545  {
546  sMv[i] = s*Mv[i];
547  }
548  }
549 
550  return sM;
551 }
552 
553 
554 template<class Form, class Type>
555 Form Foam::operator*(const Matrix<Form, Type>& M, const scalar s)
556 {
557  Form sM(M.m(), M.n());
558 
559  if (M.m() && M.n())
560  {
561  Type* sMv = sM.v();
562  const Type* Mv = M.v();
563 
564  const label mn = M.size();
565  for (label i=0; i<mn; i++)
566  {
567  sMv[i] = Mv[i]*s;
568  }
569  }
570 
571  return sM;
572 }
573 
574 
575 template<class Form, class Type>
576 Form Foam::operator/(const Matrix<Form, Type>& M, const scalar s)
577 {
578  Form sM(M.m(), M.n());
579 
580  if (M.m() && M.n())
581  {
582  Type* sMv = sM.v();
583  const Type* Mv = M.v();
584 
585  const label mn = M.size();
586  for (label i=0; i<mn; i++)
587  {
588  sMv[i] = Mv[i]/s;
589  }
590  }
591 
592  return sM;
593 }
594 
595 
596 template<class Form1, class Form2, class Type>
598 Foam::operator*
599 (
600  const Matrix<Form1, Type>& A,
601  const Matrix<Form2, Type>& B
602 )
603 {
604  if (A.n() != B.m())
605  {
607  << "Attempt to multiply incompatible matrices:" << nl
608  << "Matrix A : " << A.m() << " x " << A.n() << nl
609  << "Matrix B : " << B.m() << " x " << B.n() << nl
610  << "In order to multiply matrices, columns of A must equal "
611  << "rows of B"
612  << abort(FatalError);
613  }
614 
616  (
617  A.m(),
618  B.n(),
619  Zero
620  );
621 
622  for (label i=0; i<AB.m(); i++)
623  {
624  for (label j=0; j<AB.n(); j++)
625  {
626  for (label k=0; k<B.m(); k++)
627  {
628  AB(i, j) += A(i, k)*B(k, j);
629  }
630  }
631  }
632 
633  return AB;
634 }
635 
636 
637 template<class Form, class Type>
638 inline Foam::tmp<Foam::Field<Type>> Foam::operator*
639 (
640  const Matrix<Form, Type>& M,
641  const Field<Type>& f
642 )
643 {
644  if (M.n() != f.size())
645  {
647  << "Attempt to multiply incompatible matrix and field:" << nl
648  << "Matrix : " << M.m() << " x " << M.n() << nl
649  << "Field : " << f.size() << " rows" << nl
650  << "In order to multiply a matrix M and field f, "
651  "columns of M must equal rows of f"
652  << abort(FatalError);
653  }
654 
655  tmp<Field<Type>> tMf(new Field<Type>(M.m(), Zero));
656  Field<Type>& Mf = tMf.ref();
657 
658  for (label i=0; i<M.m(); i++)
659  {
660  for (label j=0; j<M.n(); j++)
661  {
662  Mf[i] += M(i, j)*f[j];
663  }
664  }
665 
666  return tMf;
667 }
668 
669 
670 // * * * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * //
671 
672 #include "MatrixIO.C"
673 
674 // ************************************************************************* //
label m() const
Return the number of rows in the block.
Definition: MatrixBlockI.H:93
label n() const
Return the number of columns.
Definition: MatrixI.H:64
Abstract template class to provide the form resulting from.
Definition: products.H:47
tmp< fvMatrix< Type > > operator*(const volScalarField::Internal &, const fvMatrix< Type > &)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
tUEqn clear()
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
Form T() const
Return the transpose of the matrix.
Definition: Matrix.C:264
label m() const
Return the number of rows in the block.
Definition: MatrixBlockI.H:107
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
A templated block of an (m x n) matrix of type <MatrixType>.
Definition: Matrix.H:71
void setSize(const label m, const label n)
Resize the matrix preserving the elements.
Definition: Matrix.C:244
label k
Boltzmann constant.
tmp< fvMatrix< Type > > operator/(const fvMatrix< Type > &, const volScalarField::Internal &)
Matrix()
Null constructor.
Definition: MatrixI.H:31
label n() const
Return the number of columns in the block.
Definition: MatrixBlockI.H:114
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
label size() const
Return the number of elements in matrix (m*n)
Definition: MatrixI.H:71
Pre-declare SubField and related Field type.
Definition: Field.H:56
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
tmp< fvMatrix< Type > > operator-(const fvMatrix< Type > &)
static const zero Zero
Definition: zero.H:97
tmp< fvMatrix< Type > > operator+(const fvMatrix< Type > &, const fvMatrix< Type > &)
void clear()
Clear the Matrix, i.e. set sizes to zero.
Definition: Matrix.C:214
errorManip< error > abort(error &err)
Definition: errorManip.H:131
A templated (m x n) matrix of objects of <T>.
static const char nl
Definition: Ostream.H:260
label n() const
Return the number of columns in the block.
Definition: MatrixBlockI.H:100
labelList f(nPoints)
void operator=(const mType &)
Assignment operator. Takes linear time.
Definition: Matrix.C:284
~Matrix()
Destructor.
Definition: Matrix.C:202
label m() const
Return the number of rows.
Definition: MatrixI.H:57
A class representing the concept of 0 used to avoid unnecessary manipulations for objects that are kn...
Definition: zero.H:49
label n
A class for managing temporary objects.
Definition: PtrList.H:53
const Type * v() const
Return element vector of the constant Matrix.
Definition: MatrixI.H:118
void transfer(mType &)
Transfer the contents of the argument Matrix into this Matrix.
Definition: Matrix.C:228
#define M(I)