MatrixVector.h 11.9 KB
Newer Older
1 2 3 4
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
5
// ==  http://www.amdis-fem.org                                              ==
6 7
// ==                                                                        ==
// ============================================================================
8 9 10 11 12 13 14 15 16 17 18 19
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


20 21 22 23 24 25 26 27

/** \file MatrixVector.h */

#ifndef AMDIS_MATRIXVECTOR_H
#define AMDIS_MATRIXVECTOR_H

#include <vector>
#include "Global.h"
28
#include "AMDiS_fwd.h"
29 30 31 32
#include "Serializable.h"

namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
33
  /// Class for efficient vector operations of fixed size numerical vectors.
34 35 36 37
  template<typename T>
  class Vector : public Serializable
  {
  public:
38 39 40
    
    typedef T value_type;
    
Thomas Witkowski's avatar
Thomas Witkowski committed
41
    /// Constructor.
42
    Vector(int i = 0) 
Thomas Witkowski's avatar
Thomas Witkowski committed
43 44
      : size(i),
	valArray(NULL)
45 46 47 48 49
    {
      if (size == 0) 
	valArray = NULL;
      else
	valArray = new T[size];
50
    }
51

52 53
    inline bool used() const 
    {
54
      return (valArray != NULL);
55
    }
56

Thomas Witkowski's avatar
Thomas Witkowski committed
57
    /// Change the size of the vector to newSize.
58 59
    inline void resize(int newSize) 
    {
60
      if (size != newSize) {
Thomas Witkowski's avatar
Thomas Witkowski committed
61
	if (valArray != NULL) 
62 63 64 65
	  delete [] valArray;
	valArray = new T[newSize];
	size = newSize;
      }
66
    }
67

Thomas Witkowski's avatar
Thomas Witkowski committed
68
    /// Copy constructor.
69
    Vector(const Vector<T>& rhs) 
Thomas Witkowski's avatar
Thomas Witkowski committed
70 71
      : Serializable(),
	size(rhs.size)
72 73 74
    {
      valArray = new T[rhs.size];
      *this = rhs; // uses operator=()
75
    }
76

Thomas Witkowski's avatar
Thomas Witkowski committed
77
    /// Destructor.
78 79
    virtual ~Vector() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
80 81 82 83
      if (valArray != NULL) {
	delete [] valArray; 
	valArray = NULL;
      }
84
    }
85

Thomas Witkowski's avatar
Thomas Witkowski committed
86
    /// Assignement operator
87 88
    inline const Vector<T>& operator=(const Vector<T>& rhs) 
    {
89 90 91
      TEST_EXIT_DBG(rhs.size == size)
	("Invalid sizes %d != %d!\n", rhs.size, size);

92
      T *rhsIt, *thisIt;
93 94
      for (rhsIt = rhs.begin(), thisIt = this->begin();
	   rhsIt != rhs.end();
Thomas Witkowski's avatar
Thomas Witkowski committed
95 96
	   ++rhsIt, ++thisIt)
	*thisIt = *rhsIt;      
97

98
      return *this;
99
    }
100

Thomas Witkowski's avatar
Thomas Witkowski committed
101
    /// Assignement operator
102 103
    inline const Vector<T>& operator=(const T& scal) 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
104
      for (T *thisIt = this->begin(); thisIt != this->end(); ++thisIt)
105 106
	*thisIt = scal;

107
      return *this;
108
    }
109

Thomas Witkowski's avatar
Thomas Witkowski committed
110
    /// Assignement operator
111 112
    inline const Vector<T>& operator=(const T* vec) 
    {
113 114
      T *thisIt;
      const T *vecIt;
115 116 117 118 119
      for (thisIt = this->begin(), vecIt = &vec[0];
	   thisIt != this->end();
	   ++thisIt, ++vecIt)
	*thisIt = *vecIt;

120
      return *this;
121
    }
122

Thomas Witkowski's avatar
Thomas Witkowski committed
123
    /// Sets all entries to scal.
124 125
    inline const Vector<T>& set(const T& scal) 
    {
126
      return *this = scal;
127
    }
128

Thomas Witkowski's avatar
Thomas Witkowski committed
129
    /// Sets all entries.
130 131
    inline const Vector<T>& setValues(const T* values) 
    {
132 133
      T *thisIt;
      const T *valuesIt;
134 135 136 137 138
      for (thisIt = this->begin(), valuesIt = values; 
	   thisIt != this->end(); 
	   ++thisIt, ++valuesIt) 
	*thisIt = *valuesIt;
	
139
      return *this;
140
    }
141

Thomas Witkowski's avatar
Thomas Witkowski committed
142
    /// Sets all entries.
143 144 145
    inline void fill(const T value) 
    {
      for (T *thisIt = this->begin(); thisIt != this->end(); thisIt++)
146 147 148
	*thisIt = value;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
149
    /// Comparison operator.
150 151 152 153 154
    inline bool operator==(const Vector<T>& rhs) const 
    {
      if (size != rhs.size) 
	return false;

155
      T *rhsIt, *thisIt;
156 157 158 159 160 161
      for (rhsIt = rhs.begin(), thisIt = this->begin();
	   rhsIt != rhs.end();
	   ++rhsIt, ++thisIt)
	if(*thisIt != *rhsIt) 
	  return false;

162
      return true;
163
    }
164

Thomas Witkowski's avatar
Thomas Witkowski committed
165
    /// Comparison operator.
166 167
    inline bool operator!=(const Vector<T>& rhs) const 
    {
168 169 170
      return !(*this==rhs);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
171
    /// Access to the i-th vector element.
172 173
    inline T& operator[](int i) 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
174
      TEST_EXIT_DBG(i < size && i >= 0)("Invalid index %d!\n", i);
175
      return valArray[i];
176
    }
177

Thomas Witkowski's avatar
Thomas Witkowski committed
178
    /// Access to the i-th vector element for const vectors.
179 180
    inline const T& operator[] (int i) const 
    {
181
      TEST_EXIT_DBG(i < size && i >= 0)("Invalid index %d!\n", i);
182
      return valArray[i];
183
    }
184

Thomas Witkowski's avatar
Thomas Witkowski committed
185
    /// Returns pointer to the first vector element.
186 187
    inline T *begin() const 
    { 
188
      return valArray; 
189
    }
190

Thomas Witkowski's avatar
Thomas Witkowski committed
191
    /// Returns pointer after the last vector element.
192 193
    inline T *end() const 
    { 
194
      return valArray + size; 
195
    }
196

Thomas Witkowski's avatar
Thomas Witkowski committed
197
    /// Returns \ref size.
198 199
    virtual int getSize() const 
    { 
200
      return size; 
201
    }
202

Thomas Witkowski's avatar
Thomas Witkowski committed
203
    /// Returns \ref valArray as T-array
204 205
    inline T *getValArray() 
    { 
206
      return valArray; 
207
    }
208

209 210
    void print() const 
    {
211
      std::cout << this->size << " vector: " << std::endl;
212
      for (int i = 0; i < size; i++) 
213 214
	std::cout << this->valArray[i] << " ";
      std::cout << std::endl;
215
    }
216

217 218
    void serialize(std::ostream &out) 
    {
219 220
      out.write(reinterpret_cast<const char*>(&size), sizeof(int));
      out.write(reinterpret_cast<const char*>(valArray), size * sizeof(T));
221
    }
222

223 224
    void deserialize(std::istream &in) 
    {
225 226
      in.read(reinterpret_cast<char*>(&size), sizeof(int));
      in.read(reinterpret_cast<char*>(valArray), size * sizeof(T));
227
    }
228

229 230
    std::string getTypeName() const 
    { 
231
      return "Vector"; 
232
    }
233 234

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
235
    /// Size of the vector.
236 237
    int size;

Thomas Witkowski's avatar
Thomas Witkowski committed
238
    /// Internal storage pointer.
239 240 241 242
    T *valArray;
  };


Thomas Witkowski's avatar
Thomas Witkowski committed
243
  /// Class for efficient matrix operations of fixed size numerical matrices.
244 245 246 247
  template<typename T>
  class Matrix : public Vector<T>
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
248
    /// Constructor.
249
    Matrix(int r, int c)
Thomas Witkowski's avatar
Thomas Witkowski committed
250
      : Vector<T>(r * c), 
251 252
	rows(r),
	cols(c)
253
    {}
254

Thomas Witkowski's avatar
Thomas Witkowski committed
255
    /// Changes the size of the matrix to newRows x newCols.
256 257
    inline void resize(int newRows, int newCols) 
    {
258 259 260 261 262
      if ((newRows != rows) || (newCols != cols)) {
	Vector<T>::resize(newRows * newCols);
	rows = newRows;
	cols = newCols;
      }
263
    }
264

Thomas Witkowski's avatar
Thomas Witkowski committed
265
    /// Assignement operator.
266 267
    inline const Matrix<T>& operator=(const T& scal) 
    {
268
      return static_cast<const Matrix<T>&>(Vector<T>::operator=(scal));
269
    }
270

Thomas Witkowski's avatar
Thomas Witkowski committed
271
    ///
272 273
    inline bool operator==(const Matrix<T>& rhs) const 
    {
274 275 276
      if (rows != rhs.getNumRows()) return false;
      if (cols != rhs.getNumCols()) return false;
      return Vector<T>::operator == (rhs);
277
    }
278

Thomas Witkowski's avatar
Thomas Witkowski committed
279
    /// Comparison operator.
280 281
    inline bool operator!=(const Matrix<T>& rhs) const 
    {
282
      return !(*this == rhs);
283 284
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
285
    /// Access to i-th matrix row.
286 287
    inline T *operator[](int i) 
    {
288
      return this->valArray + cols * i;
289
    }
290

Thomas Witkowski's avatar
Thomas Witkowski committed
291
    /// Access to i-th matrix row for constant matrices.
292 293
    inline const T *operator[](int i) const 
    {
294
      return this->valArray + cols * i;
295
    }
296

Thomas Witkowski's avatar
Thomas Witkowski committed
297
    /// Returns \ref rows.
298 299
    inline int getNumRows() const 
    { 
300
      return rows; 
301
    }
302

Thomas Witkowski's avatar
Thomas Witkowski committed
303
    /// Return \ref cols.
304 305
    inline int getNumCols() const 
    { 
306
      return cols; 
Thomas Witkowski's avatar
Thomas Witkowski committed
307
    }
308

Thomas Witkowski's avatar
Thomas Witkowski committed
309
    /// Returns \ref rows.
310 311
    inline int getSize() const 
    {
312
      return rows; 
Thomas Witkowski's avatar
Thomas Witkowski committed
313
    }
314

Thomas Witkowski's avatar
Thomas Witkowski committed
315
    /// Returns pointer after the last vector element.
316 317
    inline T *end() const 
    { 
318
      return this->valArray + (cols * rows); 
319
    }
320

321 322
    void print() const 
    {
323
      std::cout << this->rows << " x " << this->cols << " matrix: " << std::endl;
324
      for (int i = 0; i < rows; i++) {
325
	for (int j = 0; j < cols; j++)
326 327
	  std::cout << this->valArray[i * cols + j] << " ";
	std::cout << std::endl;
328
      }
329 330
    }

331
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
332
    /// Number of matrix rows.
333 334
    int rows;

Thomas Witkowski's avatar
Thomas Witkowski committed
335
    /// Number of matrix columns.
336 337 338
    int cols;
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
339
  /// Matrix vector multiplication.
340 341 342
  template<typename T>
  inline const Vector<T>& mv(const Matrix<T>& m, const Vector<T>& v, Vector<T>& result)
  {
343 344
    TEST_EXIT_DBG(m.getNumCols() == v.getSize())("m and v not compatible\n");
    TEST_EXIT_DBG(v.getSize() == result.getSize())("wrong result size\n");
345 346 347 348 349

    T *resultIt, *mIt, *vIt;

    for (resultIt = result.begin(), mIt = m.begin(); 
	 resultIt != result.end(); 
350 351
	 ++resultIt) {
      *resultIt = 0;
352
      for (vIt = v.begin(); vIt != v.end(); ++vIt, ++mIt)
353 354
	*resultIt += *mIt * *vIt;
    }
355 356 357 358

    return result;
  }

359
  /// Matrix vector multiplication.
360
  template<typename T>
361 362
  inline const Vector<T>& operator*=(const Vector<T>& v, const Matrix<T>& m) 
  {
363
    return mv(m, v, v);
364
  }
365

366
  /// Matrix vector multiplication.
367
  template<typename T>
368 369
  inline Vector<T> operator*(const Matrix<T>& m, const Vector<T>& v) 
  {
370 371
    Vector<T> result(v.getSize());
    return mv(m, v, result);
372
  }
373

374
  /// Scalar product.
375
  template<typename T> 
376
  inline double operator*(const Vector<T>& v1, const Vector<T>& v2) 
377
  {
378
    double result = 0.0;
379 380 381 382

    T *v1It, *v2It;
    for (v1It = v1.begin(), v2It = v2.begin();
	 v1It != v1.end();
383
	 ++v1It, ++v2It)
384
      result += *v1It * *v2It;
385 386

    return result;
387
  }
388

389
  /// Vector addition.
390 391 392
  template<typename T> 
  inline const Vector<T>& add(const Vector<T>& v1, const Vector<T>& v2, Vector<T>& result)
  {
393 394
    TEST_EXIT_DBG(v1.getSize() == v2.getSize())("invalid size in test v1 == v2\n");
    TEST_EXIT_DBG(v2.getSize() == result.getSize())("invalid size in test v2 == result\n");
395 396 397 398
    T *v1It, *v2It, *resultIt;
    for (v1It = v1.begin(), v2It = v2.begin(), resultIt = result.begin();
	 v1It != v1.end();
	 ++v1It, ++v2It, ++resultIt)
399
      *resultIt = *v1It + *v2It;
400 401

    return result;
402
  }
403

404
  /// scalar * vector
405 406
  template<typename T, typename S>
  inline const Vector<T>& mult(const S& scal,
407 408 409
			       const Vector<T>& v, 
			       Vector<T>& result)
  {
410
    TEST_EXIT_DBG(v.getSize() == result.getSize())("invalid size\n");
411

412 413 414
    T *vIt, *resultIt;
    for (vIt = v.begin(), resultIt = result.begin();
	 vIt != v.end();
415
	 ++vIt, ++resultIt) 
416
      *resultIt = scal * *vIt;
417 418

    return result;
419
  }
420

421
  /// vector + scalar
422 423 424
  template<typename T>
  inline const Vector<T>& add(const Vector<T>& v, const T& scal, Vector<T>& result)
  {
425
    TEST_EXIT_DBG(v.getSize() == result.getSize())("invalid size\n");
426 427 428 429
    T *vIt, *resultIt;
    for (vIt = v.begin(), resultIt = result.begin();
	 vIt != v.end();
	 ++vIt, ++resultIt)
430
      *resultIt = *vIt + scal;
431 432

    return result;
433
  }
434

435
  /// y = a * x + y.
436 437
  template<typename T, typename S>
  inline const Vector<T>& axpy(const S& a,
438 439 440
			       const Vector<T> &x,
			       Vector<T> &y)
  {
441
    TEST_EXIT_DBG(x.getSize() == y.getSize())("invalid size\n");
442 443 444
    T *xIt, *yIt;
    for (xIt = x.begin(), yIt = y.begin();
	 xIt != x.end();
445
	 ++xIt, ++yIt) 
446
      *yIt += a * *xIt;
447 448

    return y;
449
  }
450

451 452
  template<typename T, typename S>
  inline const Vector<T>& operator*=(Vector<T>& v, const S& scal)
453
  {
454
    return mult(scal, v, v);
455
  }
456 457

  template<typename T>
458 459
  inline Vector<T> operator*(const Vector<T>& v, const T& scal) 
  {
460 461 462
    Vector<T> result = v;
    result *= scal;
    return result;
463
  }
464 465

  template<typename T>
466 467
  inline const Vector<T>& operator+(const Vector<T>& v1, const T& scal) 
  {
468 469
    Vector<T> result(v1.getSize());
    return add(v1, scal, result);
470
  }
471 472

  template<typename T>
473 474
  inline const Vector<T>& operator+=(Vector<T>& v1, const Vector<T>& v2) 
  {
475
    return add(v1, v2, v1);
476
  }
477 478

  template<typename T>
479 480
  inline Vector<T> operator+(const Vector<T>& v1, const Vector<T>& v2) 
  {
481 482 483
    Vector<T> result = v1;
    result += v2;
    return result;
484
  }
485 486

  template<typename T>
487 488
  const Vector<T>& operator-=(Vector<T>& v1, const Vector<T>& v2)
  {
489
    return axpy(-1.0, v2, v1);
490
  }
491 492

  template<typename T>
493 494
  Vector<T> operator-(const Vector<T>& v1, const Vector<T>& v2)
  {
495 496 497
    Vector<T> result = v1;
    result -= v2;
    return result;
498
  }
499 500 501 502 503 504

  template<typename T>
  inline double norm(const Vector<T> *v)
  {
    T *vIt;
    double result = 0;
505
    for (vIt = v->begin(); vIt != v->end(); ++vIt)
506 507
      result += *vIt * *vIt;
    return sqrt(result);
508
  }
509 510 511 512 513 514

  template<typename T>
  void vectorProduct(const Vector<T>& x, 
		     const Vector<T>& y, 
		     Vector<T>& z)
  {
515 516
    FUNCNAME("vectorProduct()");
    TEST_EXIT_DBG(Global::getGeo(WORLD) == 3)("DIM_OF_WORLD != 3\n");
517 518 519
    z[0] = x[1] * y[2] - x[2] * y[1];
    z[1] = x[2] * y[0] - x[0] * y[2];
    z[2] = x[0] * y[1] - x[1] * y[0];
520
  }
521 522 523
}

#endif // AMDIS_MATRIXVECTOR_H