src/OpenFOAM/matrices/Matrix/Matrix.C

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