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
85
86
      TEST_EXIT_DBG(rhs.size == size)
	("Invalid sizes %d != %d!\n", rhs.size, size);

87
      T *rhsIt, *thisIt;
88
89
      for (rhsIt = rhs.begin(), thisIt = this->begin();
	   rhsIt != rhs.end();
Thomas Witkowski's avatar
Thomas Witkowski committed
90
91
	   ++rhsIt, ++thisIt)
	*thisIt = *rhsIt;      
92

93
      return *this;
94
    }
95

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

102
      return *this;
103
    }
104

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

115
      return *this;
116
    }
117

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

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

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

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

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

157
      return true;
158
    }
159

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    T *resultIt, *mIt, *vIt;

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

    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;
    }

373
    return result;
374
  }
375

376
  /// Computes x(Ay)^T, with x and y vectors, and A a matrix.
377
  template<typename T>
378
  inline T xAy(const Vector<T>& x, const Matrix<T>& a, const Vector<T>& y) 
379
  {
380
381
382
383
384
385
386
387
388
    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;
389
      for (yIt = y.begin(); yIt != y.end(); ++yIt, ++aIt)
390
391
392
393
394
395
396
	tmp += *aIt * *yIt;
      result += *xIt * tmp;
    }

    return result;
  }

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

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

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

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

    return result;
425
  }
426

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

    return result;
440
  }
441

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

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

    return result;
457
  }
458

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

    return result;
471
  }
472

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

    return y;
487
  }
488
489

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

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

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

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

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

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

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

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

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

}

#endif // AMDIS_MATRIXVECTOR_H