MatrixVector.h 12.9 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
39
40
41
42
43
44
    Vector(int i = 0) 
      : size(i) 
    {
      if (size == 0) 
	valArray = NULL;
      else
	valArray = new T[size];
45
    }
46

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
71
    /// Destructor.
72
73
    virtual ~Vector() 
    { 
74
      delete [] valArray; 
75
    }
76

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

88
      return *this;
89
    }
90

Thomas Witkowski's avatar
Thomas Witkowski committed
91
    /// Assignement operator
92
93
    inline const Vector<T>& operator=(const T& scal) 
    {
94
      T *thisIt;
95
96
97
98
99
      for (thisIt = this->begin();
	   thisIt != this->end();
	   ++thisIt)
	*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
    return result;
349
  }
350

351
  /// Computes x(Ay)^T, with x and y vectors, and A a matrix.
352
  template<typename T>
353
  inline T xAy(const Vector<T>& x, const Matrix<T>& a, const Vector<T>& y) 
354
  {
355
356
357
358
359
360
361
362
363
    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;
364
      for (yIt = y.begin(); yIt != y.end(); ++yIt, ++aIt)
365
366
367
368
369
370
371
	tmp += *aIt * *yIt;
      result += *xIt * tmp;
    }

    return result;
  }

372
  /// Matrix vector multiplication.
373
  template<typename T>
374
375
  inline const Vector<T>& operator*=(const Vector<T>& v, const Matrix<T>& m) 
  {
376
    return mv(m, v, v);
377
  }
378

379
  /// Matrix vector multiplication.
380
  template<typename T>
381
382
  inline Vector<T> operator*(const Matrix<T>& m, const Vector<T>& v) 
  {
383
384
    Vector<T> result(v.getSize());
    return mv(m, v, result);
385
  }
386

387
  /// Scalar product.
388
  template<typename T> 
389
  inline double operator*(const Vector<T>& v1, const Vector<T>& v2) 
390
  {
391
    double result = 0.0;
392
393
394
395

    T *v1It, *v2It;
    for (v1It = v1.begin(), v2It = v2.begin();
	 v1It != v1.end();
396
	 ++v1It, ++v2It)
397
      result += *v1It * *v2It;
398
399

    return result;
400
  }
401

402
  /// Vector addition.
403
404
405
  template<typename T> 
  inline const Vector<T>& add(const Vector<T>& v1, const Vector<T>& v2, Vector<T>& result)
  {
406
407
    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");
408
409
410
411
    T *v1It, *v2It, *resultIt;
    for (v1It = v1.begin(), v2It = v2.begin(), resultIt = result.begin();
	 v1It != v1.end();
	 ++v1It, ++v2It, ++resultIt)
412
      *resultIt = *v1It + *v2It;
413
414

    return result;
415
  }
416

417
  /// scalar * vector
418
419
420
421
422
  template<typename T>
  inline const Vector<T>& mult(const T& scal, 
			       const Vector<T>& v, 
			       Vector<T>& result)
  {
423
    TEST_EXIT_DBG(v.getSize() == result.getSize())("invalid size\n");
424

425
426
427
    T *vIt, *resultIt;
    for (vIt = v.begin(), resultIt = result.begin();
	 vIt != v.end();
428
	 ++vIt, ++resultIt) 
429
      *resultIt = scal * *vIt;
430
431

    return result;
432
  }
433

434
  /// vector + scalar
435
436
437
  template<typename T>
  inline const Vector<T>& add(const Vector<T>& v, const T& scal, Vector<T>& result)
  {
438
    TEST_EXIT_DBG(v.getSize() == result.getSize())("invalid size\n");
439
440
441
442
    T *vIt, *resultIt;
    for (vIt = v.begin(), resultIt = result.begin();
	 vIt != v.end();
	 ++vIt, ++resultIt)
443
      *resultIt = *vIt + scal;
444
445

    return result;
446
  }
447

448
  /// y = a * x + y.
449
450
451
452
453
  template<typename T>
  inline const Vector<T>& axpy(const T& a,
			       const Vector<T> &x,
			       Vector<T> &y)
  {
454
    TEST_EXIT_DBG(x.getSize() == y.getSize())("invalid size\n");
455
456
457
    T *xIt, *yIt;
    for (xIt = x.begin(), yIt = y.begin();
	 xIt != x.end();
458
	 ++xIt, ++yIt) 
459
      *yIt += a * *xIt;
460
461

    return y;
462
  }
463
464

  template<typename T>
465
466
  inline const Vector<T>& operator*=(Vector<T>& v, const T& scal) 
  {
467
    return mult(scal, v, v);
468
  }
469
470

  template<typename T>
471
472
  inline Vector<T> operator*(const Vector<T>& v, const T& scal) 
  {
473
474
475
    Vector<T> result = v;
    result *= scal;
    return result;
476
  }
477
478

  template<typename T>
479
480
  inline const Vector<T>& operator+(const Vector<T>& v1, const T& scal) 
  {
481
482
    Vector<T> result(v1.getSize());
    return add(v1, scal, result);
483
  }
484
485

  template<typename T>
486
487
  inline const Vector<T>& operator+=(Vector<T>& v1, const Vector<T>& v2) 
  {
488
    return add(v1, v2, v1);
489
  }
490
491

  template<typename T>
492
493
  inline Vector<T> operator+(const Vector<T>& v1, const Vector<T>& v2) 
  {
494
495
496
    Vector<T> result = v1;
    result += v2;
    return result;
497
  }
498
499

  template<typename T>
500
501
  const Vector<T>& operator-=(Vector<T>& v1, const Vector<T>& v2)
  {
502
    return axpy(-1.0, v2, v1);
503
  }
504
505

  template<typename T>
506
507
  Vector<T> operator-(const Vector<T>& v1, const Vector<T>& v2)
  {
508
509
510
    Vector<T> result = v1;
    result -= v2;
    return result;
511
  }
512
513
514
515
516
517

  template<typename T>
  inline double norm(const Vector<T> *v)
  {
    T *vIt;
    double result = 0;
518
    for (vIt = v->begin(); vIt != v->end(); ++vIt)
519
520
      result += *vIt * *vIt;
    return sqrt(result);
521
  }
522
523
524
525
526
527

  template<typename T>
  void vectorProduct(const Vector<T>& x, 
		     const Vector<T>& y, 
		     Vector<T>& z)
  {
528
529
    FUNCNAME("vectorProduct()");
    TEST_EXIT_DBG(Global::getGeo(WORLD) == 3)("DIM_OF_WORLD != 3\n");
530
531
532
    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];
533
  }
534
535
536
537

}

#endif // AMDIS_MATRIXVECTOR_H