src/OpenFOAM/matrices/Matrix/Matrix.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "Matrix.H"
  28
  29 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
  30
  31 template<class Form, class Type>
  32 void Foam::Matrix<Form, Type>::allocate()
  33 {
  34     if (n_ && m_)
  35     {
  36         v_ = new Type*[n_];
  37         v_[0] = new Type[n_*m_];
  38
  39         for (register label i=1; i<n_; i++)
  40         {
  41             v_[i] = v_[i-1] + m_;
  42         }
  43     }
  44 }
  45
  46
  47 // * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * * //
  48
  49 template<class Form, class Type>
  50 Foam::Matrix<Form, Type>::~Matrix()
  51 {
  52     if (v_)
  53     {
  54         delete[] (v_[0]);
  55         delete[] v_;
  56     }
  57 }
  58
  59
  60 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  61
  62 template<class Form, class Type>
  63 Foam::Matrix<Form, Type>::Matrix(const label n, const label m)
  64 :
  65     v_(NULL),
  66     n_(n),
  67     m_(m)
  68 {
  69     if (n_ < 0 || m_ < 0)
  70     {
  71         FatalErrorIn
  72         (
  73             "Matrix<Form, Type>::Matrix(const label n, const label m)"
  74         )   << "bad n, m " << n_ << ", " << m_
  75             << abort(FatalError);
  76     }
  77
  78     allocate();
  79 }
  80
  81
  82 template<class Form, class Type>
  83 Foam::Matrix<Form, Type>::Matrix(const label n, const label m, const Type& a)
  84 :
  85     v_(NULL),
  86     n_(n),
  87     m_(m)
  88 {
  89     if (n_ < 0 || m_ < 0)
  90     {
  91         FatalErrorIn
  92         (
  93             "Matrix<Form, Type>::Matrix"
  94             "(const label n, const label m, const T&)"
  95         )   << "bad n, m " << n_ << ", " << m_
  96             << abort(FatalError);
  97     }
  98
  99     allocate();
 100
 101     if (v_)
 102     {
 103         Type* v = v_[0];
 104
 105         label nm = n_*m_;
 106
 107         for (register label i=0; i<nm; i++)
 108         {
 109             v[i] = a;
 110         }
 111     }
 112 }
 113
 114
 115 template<class Form, class Type>
 116 Foam::Matrix<Form, Type>::Matrix(const Matrix<Form, Type>& a)
 117 :
 118     v_(NULL),
 119     n_(a.n_),
 120     m_(a.m_)
 121 {
 122     if (a.v_)
 123     {
 124         allocate();
 125         Type* v = v_[0];
 126         const Type* av = a.v_[0];
 127
 128         label nm = n_*m_;
 129         for (register label i=0; i<nm; i++)
 130         {
 131             v[i] = av[i];
 132         }
 133     }
 134 }
 135
 136
 137 template<class Form, class Type>
 138 void Foam::Matrix<Form, Type>::clear()
 139 {
 140     if (v_)
 141     {
 142         delete[] (v_[0]);
 143         delete[] v_;
 144     }
 145     n_ = 0;
 146     m_ = 0;
 147     v_ = NULL;
 148 }
 149
 150
 151 template<class Form, class Type>
 152 void Foam::Matrix<Form, Type>::transfer(Matrix<Form, Type>& a)
 153 {
 154     clear();
 155
 156     n_ = a.n_;
 157     a.n_ = 0;
 158
 159     m_ = a.m_;
 160     a.m_ = 0;
 161
 162     v_ = a.v_;
 163     a.v_ = NULL;
 164 }
 165
 166
 167 template<class Form, class Type>
 168 Form Foam::Matrix<Form, Type>::T() const
 169 {
 170     const Matrix<Form, Type>& A = *this;
 171     Form At(m(), n());
 172
 173     for (register label i=0; i<n(); i++)
 174     {
 175         for (register label j=0; j<m(); j++)
 176         {
 177             At[j][i] = A[i][j];
 178         }
 179     }
 180
 181     return At;
 182 }
 183
 184
 185 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 186
 187 template<class Form, class Type>
 188 void Foam::Matrix<Form, Type>::operator=(const Type& t)
 189 {
 190     if (v_)
 191     {
 192         Type* v = v_[0];
 193
 194         label nm = n_*m_;
 195         for (register label i=0; i<nm; i++)
 196         {
 197             v[i] = t;
 198         }
 199     }
 200 }
 201
 202
 203 // Assignment operator. Takes linear time.
 204 template<class Form, class Type>
 205 void Foam::Matrix<Form, Type>::operator=(const Matrix<Form, Type>& a)
 206 {
 207     if (this == &a)
 208     {
 209         FatalErrorIn
 210         (
 211             "Matrix<Form, Type>::operator=(const Matrix<Form, Type>&)"
 212         )   << "attempted assignment to self"
 213             << abort(FatalError);
 214     }
 215
 216     if (n_ != a.n_ || m_ != a.m_)
 217     {
 218         clear();
 219         n_ = a.n_;
 220         m_ = a.m_;
 221         allocate();
 222     }
 223
 224     if (v_)
 225     {
 226         Type* v = v_[0];
 227         const Type* av = a.v_[0];
 228
 229         label nm = n_*m_;
 230         for (register label i=0; i<nm; i++)
 231         {
 232             v[i] = av[i];
 233         }
 234     }
 235 }
 236
 237
 238 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 239
 240 template<class Form, class Type>
 241 const Type& Foam::max(const Matrix<Form, Type>& a)
 242 {
 243     label nm = a.n()*a.m();
 244
 245     if (nm)
 246     {
 247         label curMaxI = 0;
 248         const Type* v = a[0];
 249
 250         for (register label i=1; i<nm; i++)
 251         {
 252             if (v[i] > v[curMaxI])
 253             {
 254                 curMaxI = i;
 255             }
 256         }
 257
 258         return v[curMaxI];
 259     }
 260     else
 261     {
 262         FatalErrorIn("max(const Matrix<Form, Type>&)")
 263             << "matrix is empty"
 264             << abort(FatalError);
 265
 266         // Return in error to keep compiler happy
 267         return a[0][0];
 268     }
 269 }
 270
 271
 272 template<class Form, class Type>
 273 const Type& Foam::min(const Matrix<Form, Type>& a)
 274 {
 275     label nm = a.n()*a.m();
 276
 277     if (nm)
 278     {
 279         label curMinI = 0;
 280         const Type* v = a[0];
 281
 282         for (register label i=1; i<nm; i++)
 283         {
 284             if (v[i] < v[curMinI])
 285             {
 286                 curMinI = i;
 287             }
 288         }
 289
 290         return v[curMinI];
 291     }
 292     else
 293     {
 294         FatalErrorIn("min(const Matrix<Form, Type>&)")
 295             << "matrix is empty"
 296             << abort(FatalError);
 297
 298         // Return in error to keep compiler happy
 299         return a[0][0];
 300     }
 301 }
 302
 303
 304 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 305
 306 template<class Form, class Type>
 307 Form Foam::operator-(const Matrix<Form, Type>& a)
 308 {
 309     Form na(a.n(), a.m());
 310
 311     if (a.n() && a.m())
 312     {
 313         Type* nav = na[0];
 314         const Type* av = a[0];
 315
 316         label nm = a.n()*a.m();
 317         for (register label i=0; i<nm; i++)
 318         {
 319             nav[i] = -av[i];
 320         }
 321     }
 322
 323     return na;
 324 }
 325
 326
 327 template<class Form, class Type>
 328 Form Foam::operator+(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
 329 {
 330     if (a.n() != b.n())
 331     {
 332         FatalErrorIn
 333         (
 334             "Matrix<Form, Type>::operator+(const Matrix<Form, Type>&, "
 335             "const Matrix<Form, Type>&)"
 336         )   << "attempted add matrices with different number of rows: "
 337             << a.n() << ", " << b.n()
 338             << abort(FatalError);
 339     }
 340
 341     if (a.m() != b.m())
 342     {
 343         FatalErrorIn
 344         (
 345             "Matrix<Form, Type>::operator+(const Matrix<Form, Type>&, "
 346             "const Matrix<Form, Type>&)"
 347         )   << "attempted add matrices with different number of columns: "
 348             << a.m() << ", " << b.m()
 349             << abort(FatalError);
 350     }
 351
 352     Form ab(a.n(), a.m());
 353
 354     Type* abv = ab[0];
 355     const Type* av = a[0];
 356     const Type* bv = b[0];
 357
 358     label nm = a.n()*a.m();
 359     for (register label i=0; i<nm; i++)
 360     {
 361         abv[i] = av[i] + bv[i];
 362     }
 363
 364     return ab;
 365 }
 366
 367
 368 template<class Form, class Type>
 369 Form Foam::operator-(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
 370 {
 371     if (a.n() != b.n())
 372     {
 373         FatalErrorIn
 374         (
 375             "Matrix<Form, Type>::operator-(const Matrix<Form, Type>&, "
 376             "const Matrix<Form, Type>&)"
 377         )   << "attempted add matrices with different number of rows: "
 378             << a.n() << ", " << b.n()
 379             << abort(FatalError);
 380     }
 381
 382     if (a.m() != b.m())
 383     {
 384         FatalErrorIn
 385         (
 386             "Matrix<Form, Type>::operator-(const Matrix<Form, Type>&, "
 387             "const Matrix<Form, Type>&)"
 388         )   << "attempted add matrices with different number of columns: "
 389             << a.m() << ", " << b.m()
 390             << abort(FatalError);
 391     }
 392
 393     Form ab(a.n(), a.m());
 394
 395     Type* abv = ab[0];
 396     const Type* av = a[0];
 397     const Type* bv = b[0];
 398
 399     label nm = a.n()*a.m();
 400     for (register label i=0; i<nm; i++)
 401     {
 402         abv[i] = av[i] - bv[i];
 403     }
 404
 405     return ab;
 406 }
 407
 408
 409 template<class Form, class Type>
 410 Form Foam::operator*(const scalar s, const Matrix<Form, Type>& a)
 411 {
 412     Form sa(a.n(), a.m());
 413
 414     if (a.n() && a.m())
 415     {
 416         Type* sav = sa[0];
 417         const Type* av = a[0];
 418
 419         label nm = a.n()*a.m();
 420         for (register label i=0; i<nm; i++)
 421         {
 422             sav[i] = s*av[i];
 423         }
 424     }
 425
 426     return sa;
 427 }
 428
 429
 430 // * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //
 431
 432 #include "MatrixIO.C"
 433
 434 // ************************************************************************* //