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("Matrix<Form, Type>::operator=(const Matrix<Form, Type>&)")
 210             << "attempted assignment to self"
 211             << abort(FatalError);
 212     }
 213
 214     if (n_ != a.n_ || m_ != a.m_)
 215     {
 216         clear();
 217         n_ = a.n_;
 218         m_ = a.m_;
 219         allocate();
 220     }
 221
 222     if (v_)
 223     {
 224         Type* v = v_[0];
 225         const Type* av = a.v_[0];
 226
 227         label nm = n_*m_;
 228         for (register label i=0; i<nm; i++)
 229         {
 230             v[i] = av[i];
 231         }
 232     }
 233 }
 234
 235
 236 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 237
 238 template<class Form, class Type>
 239 const Type& Foam::max(const Matrix<Form, Type>& a)
 240 {
 241     label nm = a.n()*a.m();
 242
 243     if (nm)
 244     {
 245         label curMaxI = 0;
 246         const Type* v = a[0];
 247
 248         for (register label i=1; i<nm; i++)
 249         {
 250             if (v[i] > v[curMaxI])
 251             {
 252                 curMaxI = i;
 253             }
 254         }
 255
 256         return v[curMaxI];
 257     }
 258     else
 259     {
 260         FatalErrorIn("max(const Matrix<Form, Type>&)")
 261             << "matrix is empty"
 262             << abort(FatalError);
 263
 264         // Return in error to keep compiler happy
 265         return a[0][0];
 266     }
 267 }
 268
 269
 270 template<class Form, class Type>
 271 const Type& Foam::min(const Matrix<Form, Type>& a)
 272 {
 273     label nm = a.n()*a.m();
 274
 275     if (nm)
 276     {
 277         label curMinI = 0;
 278         const Type* v = a[0];
 279
 280         for (register label i=1; i<nm; i++)
 281         {
 282             if (v[i] < v[curMinI])
 283             {
 284                 curMinI = i;
 285             }
 286         }
 287
 288         return v[curMinI];
 289     }
 290     else
 291     {
 292         FatalErrorIn("min(const Matrix<Form, Type>&)")
 293             << "matrix is empty"
 294             << abort(FatalError);
 295
 296         // Return in error to keep compiler happy
 297         return a[0][0];
 298     }
 299 }
 300
 301
 302 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 303
 304 template<class Form, class Type>
 305 Form Foam::operator-(const Matrix<Form, Type>& a)
 306 {
 307     Form na(a.n(), a.m());
 308
 309     if (a.n() && a.m())
 310     {
 311         Type* nav = na[0];
 312         const Type* av = a[0];
 313
 314         label nm = a.n()*a.m();
 315         for (register label i=0; i<nm; i++)
 316         {
 317             nav[i] = -av[i];
 318         }
 319     }
 320
 321     return na;
 322 }
 323
 324
 325 template<class Form, class Type>
 326 Form Foam::operator+(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
 327 {
 328     if (a.n() != b.n())
 329     {
 330         FatalErrorIn
 331         (
 332             "Matrix<Form, Type>::operator+(const Matrix<Form, Type>&, const Matrix<Form, Type>&)"
 333         )   << "attempted add matrices with different number of rows: "
 334             << a.n() << ", " << b.n()
 335             << abort(FatalError);
 336     }
 337
 338     if (a.m() != b.m())
 339     {
 340         FatalErrorIn
 341         (
 342             "Matrix<Form, Type>::operator+(const Matrix<Form, Type>&, const Matrix<Form, Type>&)"
 343         )   << "attempted add matrices with different number of columns: "
 344             << a.m() << ", " << b.m()
 345             << abort(FatalError);
 346     }
 347
 348     Form ab(a.n(), a.m());
 349
 350     Type* abv = ab[0];
 351     const Type* av = a[0];
 352     const Type* bv = b[0];
 353
 354     label nm = a.n()*a.m();
 355     for (register label i=0; i<nm; i++)
 356     {
 357         abv[i] = av[i] + bv[i];
 358     }
 359
 360     return ab;
 361 }
 362
 363
 364 template<class Form, class Type>
 365 Form Foam::operator-(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
 366 {
 367     if (a.n() != b.n())
 368     {
 369         FatalErrorIn
 370         (
 371             "Matrix<Form, Type>::operator-(const Matrix<Form, Type>&, const Matrix<Form, Type>&)"
 372         )   << "attempted add matrices with different number of rows: "
 373             << a.n() << ", " << b.n()
 374             << abort(FatalError);
 375     }
 376
 377     if (a.m() != b.m())
 378     {
 379         FatalErrorIn
 380         (
 381             "Matrix<Form, Type>::operator-(const Matrix<Form, Type>&, const Matrix<Form, Type>&)"
 382         )   << "attempted add matrices with different number of columns: "
 383             << a.m() << ", " << b.m()
 384             << abort(FatalError);
 385     }
 386
 387     Form ab(a.n(), a.m());
 388
 389     Type* abv = ab[0];
 390     const Type* av = a[0];
 391     const Type* bv = b[0];
 392
 393     label nm = a.n()*a.m();
 394     for (register label i=0; i<nm; i++)
 395     {
 396         abv[i] = av[i] - bv[i];
 397     }
 398
 399     return ab;
 400 }
 401
 402
 403 template<class Form, class Type>
 404 Form Foam::operator*(const scalar s, const Matrix<Form, Type>& a)
 405 {
 406     Form sa(a.n(), a.m());
 407
 408     if (a.n() && a.m())
 409     {
 410         Type* sav = sa[0];
 411         const Type* av = a[0];
 412
 413         label nm = a.n()*a.m();
 414         for (register label i=0; i<nm; i++)
 415         {
 416             sav[i] = s*av[i];
 417         }
 418     }
 419
 420     return sa;
 421 }
 422
 423
 424 // * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //
 425
 426 #include "MatrixIO.C"
 427
 428 // ************************************************************************* //