src/foam/matrices/Matrix/Matrix.C

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