MatrixVector.h 13 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  crystal growth group                                                  ==
// ==                                                                        ==
// ==  Stiftung caesar                                                       ==
// ==  Ludwig-Erhard-Allee 2                                                 ==
// ==  53175 Bonn                                                            ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  http://www.caesar.de/cg/AMDiS                                         ==
// ==                                                                        ==
// ============================================================================

/** \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
    inline bool used() const {
48
      return (valArray != NULL);
49
    }
50

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

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

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

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

84
      return *this;
85
    }
86

Thomas Witkowski's avatar
Thomas Witkowski committed
87
    /// Assignement operator
88
89
90
91
92
93
94
95
96
    inline const Vector<T>& operator=(const T& scal) {
      T *thisIt;
      for(thisIt = this->begin();
	  thisIt != this->end();
	  ++thisIt)
	{
	  *thisIt = scal;
	}
      return *this;
97
    }
98

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

Thomas Witkowski's avatar
Thomas Witkowski committed
112
    /// Sets all entries to scal.
113
114
    inline const Vector<T>& set(const T& scal) {
      return *this = scal;
115
    }
116

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

Thomas Witkowski's avatar
Thomas Witkowski committed
130
    /// Sets all entries.
131
132
133
134
135
136
    inline void fill(const T value) {
      for (T *thisIt = this->begin(); thisIt != this->end(); thisIt++) {
	*thisIt = value;
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
137
    /// Comparison operator.
138
139
140
141
142
143
144
145
146
147
    inline bool operator==(const Vector<T>& rhs) const {
      if(size != rhs.size) return false;
      T *rhsIt, *thisIt;
      for(rhsIt = rhs.begin(), thisIt = this->begin();
	  rhsIt != rhs.end();
	  ++rhsIt, ++thisIt)
	{
	  if(*thisIt != *rhsIt) return false;
	}
      return true;
148
    }
149

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

Thomas Witkowski's avatar
Thomas Witkowski committed
155
    /// Access to the i-th vector element.
156
    inline T& operator[](int i) {
157
      TEST_EXIT_DBG(i < size && i >= 0)("invalid index\n");
158
      return valArray[i];
159
    }
160

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

Thomas Witkowski's avatar
Thomas Witkowski committed
167
    /// Returns pointer to the first vector element.
168
169
    inline T *begin() const { 
      return valArray; 
170
    }
171

Thomas Witkowski's avatar
Thomas Witkowski committed
172
    /// Returns pointer after the last vector element.
173
174
    inline T *end() const { 
      return valArray + size; 
175
    }
176

Thomas Witkowski's avatar
Thomas Witkowski committed
177
    /// Returns \ref size.
178
179
    virtual int getSize() const { 
      return size; 
180
    }
181

Thomas Witkowski's avatar
Thomas Witkowski committed
182
    /// Returns \ref valArray as T-array
183
184
    inline T *getValArray() { 
      return valArray; 
185
    }
186
187

    void print() const {
188
      std::cout << this->size << " vector: " << std::endl;
189
      for (int i = 0; i < size; i++) {
190
	std::cout << this->valArray[i] << " ";
191
      }
192
      std::cout << std::endl;
193
    }
194

195
    void serialize(std::ostream &out) {
196
197
      out.write(reinterpret_cast<const char*>(&size), sizeof(int));
      out.write(reinterpret_cast<const char*>(valArray), size * sizeof(T));
198
    }
199

200
    void deserialize(std::istream &in) {
201
202
      in.read(reinterpret_cast<char*>(&size), sizeof(int));
      in.read(reinterpret_cast<char*>(valArray), size * sizeof(T));
203
    }
204

205
    std::string getTypeName() const { 
206
      return "Vector"; 
207
    }
208
209

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
210
    /// Size of the vector.
211
212
    int size;

Thomas Witkowski's avatar
Thomas Witkowski committed
213
    /// Internal storage pointer.
214
215
216
217
    T *valArray;
  };


Thomas Witkowski's avatar
Thomas Witkowski committed
218
  /// Class for efficient matrix operations of fixed size numerical matrices.
219
220
221
222
  template<typename T>
  class Matrix : public Vector<T>
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
223
    /// Constructor.
224
    Matrix(int r, int c)
Thomas Witkowski's avatar
Thomas Witkowski committed
225
      : Vector<T>(r * c), 
226
227
	rows(r),
	cols(c)
228
    {}
229

Thomas Witkowski's avatar
Thomas Witkowski committed
230
    /// Changes the size of the matrix to newRows x newCols.
231
232
233
234
235
236
    inline void resize(int newRows, int newCols) {
      if ((newRows != rows) || (newCols != cols)) {
	Vector<T>::resize(newRows * newCols);
	rows = newRows;
	cols = newCols;
      }
237
    }
238

Thomas Witkowski's avatar
Thomas Witkowski committed
239
    /// Assignement operator.
240
    inline const Matrix<T>& operator=(const T& scal) {
241
      return static_cast<const Matrix<T>&>(Vector<T>::operator=(scal));
242
    }
243

Thomas Witkowski's avatar
Thomas Witkowski committed
244
    ///
245
    inline bool operator==(const Matrix<T>& rhs) const {
246
247
248
      if (rows != rhs.getNumRows()) return false;
      if (cols != rhs.getNumCols()) return false;
      return Vector<T>::operator == (rhs);
249
    }
250

Thomas Witkowski's avatar
Thomas Witkowski committed
251
    /// Comparison operator.
252
253
    inline bool operator!=(const Matrix<T>& rhs) const {
      return !(*this == rhs);
254
255
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
256
    /// Acces to i-th matrix row.
257
    inline T *operator[](int i) {
258
      return this->valArray + cols * i;
259
    }
260

Thomas Witkowski's avatar
Thomas Witkowski committed
261
    /// Acces to i-th matrix row for constant matrices.
262
    inline const T *operator[](int i) const {
263
      return this->valArray + cols * i;
264
    }
265

Thomas Witkowski's avatar
Thomas Witkowski committed
266
    /// Returns \ref rows.
267
    inline int getNumRows() const { 
268
      return rows; 
269
    }
270

Thomas Witkowski's avatar
Thomas Witkowski committed
271
    /// Return \ref cols.
272
    inline int getNumCols() const { 
273
      return cols; 
Thomas Witkowski's avatar
Thomas Witkowski committed
274
    }
275

Thomas Witkowski's avatar
Thomas Witkowski committed
276
    /// Returns \ref rows.
277
    inline int getSize() const {
278
      return rows; 
Thomas Witkowski's avatar
Thomas Witkowski committed
279
    }
280

Thomas Witkowski's avatar
Thomas Witkowski committed
281
    /// Returns pointer after the last vector element.
282
283
    inline T *end() const { 
      return this->valArray + (cols * rows); 
284
    }
285
286

    void print() const {
287
      std::cout << this->rows << " x " << this->cols << " matrix: " << std::endl;
288
289
      for (int i = 0; i < rows; i++) {
	for (int j = 0; j < cols; j++) {
290
	  std::cout << this->valArray[i * cols + j] << " ";
291
	}
292
	std::cout << std::endl;
293
      }
294
295
    }

296
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
297
    /// Number of matrix rows.
298
299
    int rows;

Thomas Witkowski's avatar
Thomas Witkowski committed
300
    /// Number of matrix columns.
301
302
303
    int cols;
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
304
  /// Matrix vector multiplication.
305
306
307
  template<typename T>
  inline const Vector<T>& mv(const Matrix<T>& m, const Vector<T>& v, Vector<T>& result)
  {
308
309
    TEST_EXIT_DBG(m.getNumCols() == v.getSize())("m and v not compatible\n");
    TEST_EXIT_DBG(v.getSize() == result.getSize())("wrong result size\n");
310
311
312
313
314

    T *resultIt, *mIt, *vIt;

    for (resultIt = result.begin(), mIt = m.begin(); 
	 resultIt != result.end(); 
315
316
317
318
	 ++resultIt) {
      *resultIt = 0;
      for (vIt = v.begin(); vIt != v.end(); ++vIt, ++mIt) {
	*resultIt += *mIt * *vIt;
319
      }
320
    }
321
322
323
324
    return result;
  };

  /** \brief
325
   * Computes x(Ay)^T, with x and y vectors, and A a matrix.
326
327
   */
  template<typename T>
328
  inline T xAy(const Vector<T>& x, const Matrix<T>& a, const Vector<T>& y) 
329
  {
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
    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;
      for (yIt = y.begin(); yIt != y.end(); ++yIt, ++aIt) {
	tmp += *aIt * *yIt;
      }
      result += *xIt * tmp;
    }

    return result;
  }

  /** \brief
   * Matrix vector multiplication.
   */
  template<typename T>
  inline const Vector<T>& operator*=(const Vector<T>& v, const Matrix<T>& m) {
353
354
355
356
357
358
359
    return mv(m, v, v);
  };

  /** \brief
   * Matrix vector multiplication.
   */
  template<typename T>
360
  inline Vector<T> operator*(const Matrix<T>& m, const Vector<T>& v) {
361
362
363
364
365
366
367
368
369
    Vector<T> result(v.getSize());
    return mv(m, v, result);
  };


  /** \brief
   * Scalar product.
   */
  template<typename T> 
370
  inline double operator*(const Vector<T>& v1, const Vector<T>& v2) 
371
  {
372
    double result = 0.0;
373
374
375
376

    T *v1It, *v2It;
    for (v1It = v1.begin(), v2It = v2.begin();
	 v1It != v1.end();
377
378
379
	 ++v1It, ++v2It) {
      result += *v1It * *v2It;
    }
380
381
382
383
384
385
386
387
388
389

    return result;
  };

  /** \brief
   * Vector addition.
   */
  template<typename T> 
  inline const Vector<T>& add(const Vector<T>& v1, const Vector<T>& v2, Vector<T>& result)
  {
390
391
    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");
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
    T *v1It, *v2It, *resultIt;
    for (v1It = v1.begin(), v2It = v2.begin(), resultIt = result.begin();
	 v1It != v1.end();
	 ++v1It, ++v2It, ++resultIt)
      {
	*resultIt = *v1It + *v2It;
      }

    return result;
  };

  /** \brief
   * scalar * vector
   */
  template<typename T>
  inline const Vector<T>& mult(const T& scal, 
			       const Vector<T>& v, 
			       Vector<T>& result)
  {
411
    TEST_EXIT_DBG(v.getSize() == result.getSize())("invalid size\n");
412

413
414
415
    T *vIt, *resultIt;
    for (vIt = v.begin(), resultIt = result.begin();
	 vIt != v.end();
416
417
418
	 ++vIt, ++resultIt) {
      *resultIt = scal * *vIt;
    }
419
420
421
422
423
424
425
426
427
428

    return result;
  };

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

    return result;
  };

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

    return y;
  };

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

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

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

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

  template<typename T>
484
  inline Vector<T> operator+(const Vector<T>& v1, const Vector<T>& v2) {
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
    Vector<T> result = v1;
    result += v2;
    return result;
  };

  template<typename T>
  const Vector<T>& operator-=(Vector<T>& v1, const Vector<T>& v2){
    return axpy(-1.0, v2, v1);
  };

  template<typename T>
  Vector<T> operator-(const Vector<T>& v1, const Vector<T>& v2){
    Vector<T> result = v1;
    result -= v2;
    return result;
  };

  template<typename T>
  inline double norm(const Vector<T> *v)
  {
    T *vIt;
    double result = 0;
    for (vIt = v->begin(); vIt != v->end(); ++vIt) {
      result += *vIt * *vIt;
    }
    return sqrt(result);
  };

  template<typename T>
  void vectorProduct(const Vector<T>& x, 
		     const Vector<T>& y, 
		     Vector<T>& z)
  {
518
519
    FUNCNAME("vectorProduct()");
    TEST_EXIT_DBG(Global::getGeo(WORLD) == 3)("DIM_OF_WORLD != 3\n");
520
521
522
523
524
525
526
527
    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];
  };

}

#endif // AMDIS_MATRIXVECTOR_H