MatrixVector.h 13.6 KB
Newer Older
1
2
3
4
5
6
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
7
// ==  TU Dresden                                                            ==
8
// ==                                                                        ==
9
10
11
// ==  Institut für Wissenschaftliches Rechnen                               ==
// ==  Zellescher Weg 12-14                                                  ==
// ==  01069 Dresden                                                         ==
12
13
14
15
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
16
// ==  https://gforge.zih.tu-dresden.de/projects/amdis/                      ==
17
18
19
20
21
22
23
24
25
26
// ==                                                                        ==
// ============================================================================

/** \file MatrixVector.h */

#ifndef AMDIS_MATRIXVECTOR_H
#define AMDIS_MATRIXVECTOR_H

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

namespace AMDiS {

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

48
49
    inline bool used() const 
    {
50
      return (valArray != NULL);
51
    }
52

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
72
    /// Destructor.
73
74
    virtual ~Vector() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
75
76
77
78
      if (valArray != NULL) {
	delete [] valArray; 
	valArray = NULL;
      }
79
    }
80

Thomas Witkowski's avatar
Thomas Witkowski committed
81
    /// Assignement operator
82
83
    inline const Vector<T>& operator=(const Vector<T>& rhs) 
    {
84
      TEST_EXIT_DBG(rhs.size == size)("invalid size\n");
85
      T *rhsIt, *thisIt;
86
87
      for (rhsIt = rhs.begin(), thisIt = this->begin();
	   rhsIt != rhs.end();
Thomas Witkowski's avatar
Thomas Witkowski committed
88
89
	   ++rhsIt, ++thisIt)
	*thisIt = *rhsIt;      
90

91
      return *this;
92
    }
93

Thomas Witkowski's avatar
Thomas Witkowski committed
94
    /// Assignement operator
95
96
    inline const Vector<T>& operator=(const T& scal) 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
97
      for (T *thisIt = this->begin(); thisIt != this->end(); ++thisIt)
98
99
	*thisIt = scal;

100
      return *this;
101
    }
102

Thomas Witkowski's avatar
Thomas Witkowski committed
103
    /// Assignement operator
104
105
    inline const Vector<T>& operator=(const T* vec) 
    {
106
107
      T *thisIt;
      const T *vecIt;
108
109
110
111
112
      for (thisIt = this->begin(), vecIt = &vec[0];
	   thisIt != this->end();
	   ++thisIt, ++vecIt)
	*thisIt = *vecIt;

113
      return *this;
114
    }
115

Thomas Witkowski's avatar
Thomas Witkowski committed
116
    /// Sets all entries to scal.
117
118
    inline const Vector<T>& set(const T& scal) 
    {
119
      return *this = scal;
120
    }
121

Thomas Witkowski's avatar
Thomas Witkowski committed
122
    /// Sets all entries.
123
124
    inline const Vector<T>& setValues(const T* values) 
    {
125
126
      T *thisIt;
      const T *valuesIt;
127
128
129
130
131
      for (thisIt = this->begin(), valuesIt = values; 
	   thisIt != this->end(); 
	   ++thisIt, ++valuesIt) 
	*thisIt = *valuesIt;
	
132
      return *this;
133
    }
134

Thomas Witkowski's avatar
Thomas Witkowski committed
135
    /// Sets all entries.
136
137
138
    inline void fill(const T value) 
    {
      for (T *thisIt = this->begin(); thisIt != this->end(); thisIt++)
139
140
141
	*thisIt = value;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
142
    /// Comparison operator.
143
144
145
146
147
    inline bool operator==(const Vector<T>& rhs) const 
    {
      if (size != rhs.size) 
	return false;

148
      T *rhsIt, *thisIt;
149
150
151
152
153
154
      for (rhsIt = rhs.begin(), thisIt = this->begin();
	   rhsIt != rhs.end();
	   ++rhsIt, ++thisIt)
	if(*thisIt != *rhsIt) 
	  return false;

155
      return true;
156
    }
157

Thomas Witkowski's avatar
Thomas Witkowski committed
158
    /// Comparison operator.
159
160
    inline bool operator!=(const Vector<T>& rhs) const 
    {
161
162
163
      return !(*this==rhs);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
164
    /// Access to the i-th vector element.
165
166
    inline T& operator[](int i) 
    {
167
      TEST_EXIT_DBG(i < size && i >= 0)("invalid index\n");
168
      return valArray[i];
169
    }
170

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

Thomas Witkowski's avatar
Thomas Witkowski committed
178
    /// Returns pointer to the first vector element.
179
180
    inline T *begin() const 
    { 
181
      return valArray; 
182
    }
183

Thomas Witkowski's avatar
Thomas Witkowski committed
184
    /// Returns pointer after the last vector element.
185
186
    inline T *end() const 
    { 
187
      return valArray + size; 
188
    }
189

Thomas Witkowski's avatar
Thomas Witkowski committed
190
    /// Returns \ref size.
191
192
    virtual int getSize() const 
    { 
193
      return size; 
194
    }
195

Thomas Witkowski's avatar
Thomas Witkowski committed
196
    /// Returns \ref valArray as T-array
197
198
    inline T *getValArray() 
    { 
199
      return valArray; 
200
    }
201

202
203
    void print() const 
    {
204
      std::cout << this->size << " vector: " << std::endl;
205
      for (int i = 0; i < size; i++) 
206
207
	std::cout << this->valArray[i] << " ";
      std::cout << std::endl;
208
    }
209

210
211
    void serialize(std::ostream &out) 
    {
212
213
      out.write(reinterpret_cast<const char*>(&size), sizeof(int));
      out.write(reinterpret_cast<const char*>(valArray), size * sizeof(T));
214
    }
215

216
217
    void deserialize(std::istream &in) 
    {
218
219
      in.read(reinterpret_cast<char*>(&size), sizeof(int));
      in.read(reinterpret_cast<char*>(valArray), size * sizeof(T));
220
    }
221

222
223
    std::string getTypeName() const 
    { 
224
      return "Vector"; 
225
    }
226
227

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
228
    /// Size of the vector.
229
230
    int size;

Thomas Witkowski's avatar
Thomas Witkowski committed
231
    /// Internal storage pointer.
232
233
234
235
    T *valArray;
  };


Thomas Witkowski's avatar
Thomas Witkowski committed
236
  /// Class for efficient matrix operations of fixed size numerical matrices.
237
238
239
240
  template<typename T>
  class Matrix : public Vector<T>
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
241
    /// Constructor.
242
    Matrix(int r, int c)
Thomas Witkowski's avatar
Thomas Witkowski committed
243
      : Vector<T>(r * c), 
244
245
	rows(r),
	cols(c)
246
    {}
247

Thomas Witkowski's avatar
Thomas Witkowski committed
248
    /// Changes the size of the matrix to newRows x newCols.
249
250
    inline void resize(int newRows, int newCols) 
    {
251
252
253
254
255
      if ((newRows != rows) || (newCols != cols)) {
	Vector<T>::resize(newRows * newCols);
	rows = newRows;
	cols = newCols;
      }
256
    }
257

Thomas Witkowski's avatar
Thomas Witkowski committed
258
    /// Assignement operator.
259
260
    inline const Matrix<T>& operator=(const T& scal) 
    {
261
      return static_cast<const Matrix<T>&>(Vector<T>::operator=(scal));
262
    }
263

Thomas Witkowski's avatar
Thomas Witkowski committed
264
    ///
265
266
    inline bool operator==(const Matrix<T>& rhs) const 
    {
267
268
269
      if (rows != rhs.getNumRows()) return false;
      if (cols != rhs.getNumCols()) return false;
      return Vector<T>::operator == (rhs);
270
    }
271

Thomas Witkowski's avatar
Thomas Witkowski committed
272
    /// Comparison operator.
273
274
    inline bool operator!=(const Matrix<T>& rhs) const 
    {
275
      return !(*this == rhs);
276
277
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
278
    /// Acces to i-th matrix row.
279
280
    inline T *operator[](int i) 
    {
281
      return this->valArray + cols * i;
282
    }
283

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

Thomas Witkowski's avatar
Thomas Witkowski committed
290
    /// Returns \ref rows.
291
292
    inline int getNumRows() const 
    { 
293
      return rows; 
294
    }
295

Thomas Witkowski's avatar
Thomas Witkowski committed
296
    /// Return \ref cols.
297
298
    inline int getNumCols() const 
    { 
299
      return cols; 
Thomas Witkowski's avatar
Thomas Witkowski committed
300
    }
301

Thomas Witkowski's avatar
Thomas Witkowski committed
302
    /// Returns \ref rows.
303
304
    inline int getSize() const 
    {
305
      return rows; 
Thomas Witkowski's avatar
Thomas Witkowski committed
306
    }
307

Thomas Witkowski's avatar
Thomas Witkowski committed
308
    /// Returns pointer after the last vector element.
309
310
    inline T *end() const 
    { 
311
      return this->valArray + (cols * rows); 
312
    }
313

314
315
    void print() const 
    {
316
      std::cout << this->rows << " x " << this->cols << " matrix: " << std::endl;
317
      for (int i = 0; i < rows; i++) {
318
	for (int j = 0; j < cols; j++)
319
320
	  std::cout << this->valArray[i * cols + j] << " ";
	std::cout << std::endl;
321
      }
322
323
    }

324
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
325
    /// Number of matrix rows.
326
327
    int rows;

Thomas Witkowski's avatar
Thomas Witkowski committed
328
    /// Number of matrix columns.
329
330
331
    int cols;
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
332
  /// Matrix vector multiplication.
333
334
335
  template<typename T>
  inline const Vector<T>& mv(const Matrix<T>& m, const Vector<T>& v, Vector<T>& result)
  {
336
337
    TEST_EXIT_DBG(m.getNumCols() == v.getSize())("m and v not compatible\n");
    TEST_EXIT_DBG(v.getSize() == result.getSize())("wrong result size\n");
338
339
340
341
342

    T *resultIt, *mIt, *vIt;

    for (resultIt = result.begin(), mIt = m.begin(); 
	 resultIt != result.end(); 
343
344
	 ++resultIt) {
      *resultIt = 0;
345
      for (vIt = v.begin(); vIt != v.end(); ++vIt, ++mIt)
346
347
	*resultIt += *mIt * *vIt;
    }
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370

    return result;
  }

  /// Scaled matrix vector multiplication: result = a * (m * v)
  template<typename T>
  inline const Vector<T>& amv(double a, const Matrix<T>& m, const Vector<T>& v, 
			      Vector<T>& result)
  {
    TEST_EXIT_DBG(m.getNumCols() == v.getSize())("m and v not compatible\n");
    TEST_EXIT_DBG(v.getSize() == result.getSize())("wrong result size\n");

    T *resultIt, *mIt, *vIt;

    for (resultIt = result.begin(), mIt = m.begin(); 
	 resultIt != result.end(); 
	 ++resultIt) {
      *resultIt = 0;
      for (vIt = v.begin(); vIt != v.end(); ++vIt, ++mIt)
	*resultIt += *mIt * *vIt;
      *resultIt *= a;
    }

371
    return result;
372
  }
373

374
  /// Computes x(Ay)^T, with x and y vectors, and A a matrix.
375
  template<typename T>
376
  inline T xAy(const Vector<T>& x, const Matrix<T>& a, const Vector<T>& y) 
377
  {
378
379
380
381
382
383
384
385
386
    TEST_EXIT_DBG(a.getNumCols() == x.getSize())("A and x not compatible\n");
    TEST_EXIT_DBG(a.getNumCols() == y.getSize())("A and y not compatible\n");

    T result = 0;
    T tmp = 0;
    T *aIt, *xIt, *yIt;

    for (aIt = a.begin(), xIt = x.begin(); aIt != a.end(); ++xIt) {
      tmp = 0;
387
      for (yIt = y.begin(); yIt != y.end(); ++yIt, ++aIt)
388
389
390
391
392
393
394
	tmp += *aIt * *yIt;
      result += *xIt * tmp;
    }

    return result;
  }

395
  /// Matrix vector multiplication.
396
  template<typename T>
397
398
  inline const Vector<T>& operator*=(const Vector<T>& v, const Matrix<T>& m) 
  {
399
    return mv(m, v, v);
400
  }
401

402
  /// Matrix vector multiplication.
403
  template<typename T>
404
405
  inline Vector<T> operator*(const Matrix<T>& m, const Vector<T>& v) 
  {
406
407
    Vector<T> result(v.getSize());
    return mv(m, v, result);
408
  }
409

410
  /// Scalar product.
411
  template<typename T> 
412
  inline double operator*(const Vector<T>& v1, const Vector<T>& v2) 
413
  {
414
    double result = 0.0;
415
416
417
418

    T *v1It, *v2It;
    for (v1It = v1.begin(), v2It = v2.begin();
	 v1It != v1.end();
419
	 ++v1It, ++v2It)
420
      result += *v1It * *v2It;
421
422

    return result;
423
  }
424

425
  /// Vector addition.
426
427
428
  template<typename T> 
  inline const Vector<T>& add(const Vector<T>& v1, const Vector<T>& v2, Vector<T>& result)
  {
429
430
    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");
431
432
433
434
    T *v1It, *v2It, *resultIt;
    for (v1It = v1.begin(), v2It = v2.begin(), resultIt = result.begin();
	 v1It != v1.end();
	 ++v1It, ++v2It, ++resultIt)
435
      *resultIt = *v1It + *v2It;
436
437

    return result;
438
  }
439

440
  /// scalar * vector
441
442
443
444
445
  template<typename T>
  inline const Vector<T>& mult(const T& scal, 
			       const Vector<T>& v, 
			       Vector<T>& result)
  {
446
    TEST_EXIT_DBG(v.getSize() == result.getSize())("invalid size\n");
447

448
449
450
    T *vIt, *resultIt;
    for (vIt = v.begin(), resultIt = result.begin();
	 vIt != v.end();
451
	 ++vIt, ++resultIt) 
452
      *resultIt = scal * *vIt;
453
454

    return result;
455
  }
456

457
  /// vector + scalar
458
459
460
  template<typename T>
  inline const Vector<T>& add(const Vector<T>& v, const T& scal, Vector<T>& result)
  {
461
    TEST_EXIT_DBG(v.getSize() == result.getSize())("invalid size\n");
462
463
464
465
    T *vIt, *resultIt;
    for (vIt = v.begin(), resultIt = result.begin();
	 vIt != v.end();
	 ++vIt, ++resultIt)
466
      *resultIt = *vIt + scal;
467
468

    return result;
469
  }
470

471
  /// y = a * x + y.
472
473
474
475
476
  template<typename T>
  inline const Vector<T>& axpy(const T& a,
			       const Vector<T> &x,
			       Vector<T> &y)
  {
477
    TEST_EXIT_DBG(x.getSize() == y.getSize())("invalid size\n");
478
479
480
    T *xIt, *yIt;
    for (xIt = x.begin(), yIt = y.begin();
	 xIt != x.end();
481
	 ++xIt, ++yIt) 
482
      *yIt += a * *xIt;
483
484

    return y;
485
  }
486
487

  template<typename T>
488
489
  inline const Vector<T>& operator*=(Vector<T>& v, const T& scal) 
  {
490
    return mult(scal, v, v);
491
  }
492
493

  template<typename T>
494
495
  inline Vector<T> operator*(const Vector<T>& v, const T& scal) 
  {
496
497
498
    Vector<T> result = v;
    result *= scal;
    return result;
499
  }
500
501

  template<typename T>
502
503
  inline const Vector<T>& operator+(const Vector<T>& v1, const T& scal) 
  {
504
505
    Vector<T> result(v1.getSize());
    return add(v1, scal, result);
506
  }
507
508

  template<typename T>
509
510
  inline const Vector<T>& operator+=(Vector<T>& v1, const Vector<T>& v2) 
  {
511
    return add(v1, v2, v1);
512
  }
513
514

  template<typename T>
515
516
  inline Vector<T> operator+(const Vector<T>& v1, const Vector<T>& v2) 
  {
517
518
519
    Vector<T> result = v1;
    result += v2;
    return result;
520
  }
521
522

  template<typename T>
523
524
  const Vector<T>& operator-=(Vector<T>& v1, const Vector<T>& v2)
  {
525
    return axpy(-1.0, v2, v1);
526
  }
527
528

  template<typename T>
529
530
  Vector<T> operator-(const Vector<T>& v1, const Vector<T>& v2)
  {
531
532
533
    Vector<T> result = v1;
    result -= v2;
    return result;
534
  }
535
536
537
538
539
540

  template<typename T>
  inline double norm(const Vector<T> *v)
  {
    T *vIt;
    double result = 0;
541
    for (vIt = v->begin(); vIt != v->end(); ++vIt)
542
543
      result += *vIt * *vIt;
    return sqrt(result);
544
  }
545
546
547
548
549
550

  template<typename T>
  void vectorProduct(const Vector<T>& x, 
		     const Vector<T>& y, 
		     Vector<T>& z)
  {
551
552
    FUNCNAME("vectorProduct()");
    TEST_EXIT_DBG(Global::getGeo(WORLD) == 3)("DIM_OF_WORLD != 3\n");
553
554
555
    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];
556
  }
557
558
559
560

}

#endif // AMDIS_MATRIXVECTOR_H