MatrixVector.h 13.1 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
27
28
29
30
31
32
33
// ============================================================================
// ==                                                                        ==
// == 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"
#include "MemoryManager.h"
#include "Serializable.h"

namespace AMDiS {

  template<typename T> class DOFVector;

Thomas Witkowski's avatar
Thomas Witkowski committed
34
  /// Class for efficient vector operations of fixed size numerical vectors.
35
36
37
38
39
40
  template<typename T>
  class Vector : public Serializable
  {
  public:
    MEMORY_MANAGED(Vector<T>);

Thomas Witkowski's avatar
Thomas Witkowski committed
41
    /// Constructor.
42
43
44
45
46
47
48
    Vector(int i = 0) 
      : size(i) 
    {
      if (size == 0) 
	valArray = NULL;
      else
	valArray = new T[size];
49
    }
50

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
103
    /// Assignement operator
104
105
106
107
108
109
110
111
112
113
    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;
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) {
      return *this = scal;
119
    }
120

Thomas Witkowski's avatar
Thomas Witkowski committed
121
    /// Sets all entries.
122
123
124
125
126
127
128
129
130
131
    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;
132
    }
133

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

Thomas Witkowski's avatar
Thomas Witkowski committed
141
    /// Comparison operator.
142
143
144
145
146
147
148
149
150
151
    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;
152
    }
153

Thomas Witkowski's avatar
Thomas Witkowski committed
154
    /// Comparison operator.
155
156
157
158
    inline bool operator!=(const Vector<T>& rhs) const {
      return !(*this==rhs);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
165
    /// Access to the i-th vector element for const vectors.
166
    inline const T& operator[] (int i) const {
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
    /// Returns pointer to the first vector element.
172
173
    inline T *begin() const { 
      return valArray; 
174
    }
175

Thomas Witkowski's avatar
Thomas Witkowski committed
176
    /// Returns pointer after the last vector element.
177
178
    inline T *end() const { 
      return valArray + size; 
179
    }
180

Thomas Witkowski's avatar
Thomas Witkowski committed
181
    /// Returns \ref size.
182
183
    virtual int getSize() const { 
      return size; 
184
    }
185

Thomas Witkowski's avatar
Thomas Witkowski committed
186
    /// Returns \ref valArray as T-array
187
188
    inline T *getValArray() { 
      return valArray; 
189
    }
190
191

    void print() const {
192
      std::cout << this->size << " vector: " << std::endl;
193
      for (int i = 0; i < size; i++) {
194
	std::cout << this->valArray[i] << " ";
195
      }
196
      std::cout << std::endl;
197
    }
198

199
    void serialize(std::ostream &out) {
200
201
      out.write(reinterpret_cast<const char*>(&size), sizeof(int));
      out.write(reinterpret_cast<const char*>(valArray), size * sizeof(T));
202
    }
203

204
    void deserialize(std::istream &in) {
205
206
      in.read(reinterpret_cast<char*>(&size), sizeof(int));
      in.read(reinterpret_cast<char*>(valArray), size * sizeof(T));
207
    }
208

209
    std::string getTypeName() const { 
210
      return "Vector"; 
211
    }
212
213

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
214
    /// Size of the vector.
215
216
    int size;

Thomas Witkowski's avatar
Thomas Witkowski committed
217
    /// Internal storage pointer.
218
219
220
221
    T *valArray;
  };


Thomas Witkowski's avatar
Thomas Witkowski committed
222
  /// Class for efficient matrix operations of fixed size numerical matrices.
223
224
225
226
227
228
  template<typename T>
  class Matrix : public Vector<T>
  {
  public:
    MEMORY_MANAGED(Matrix<T>);

Thomas Witkowski's avatar
Thomas Witkowski committed
229
    /// Constructor.
230
    Matrix(int r, int c)
Thomas Witkowski's avatar
Thomas Witkowski committed
231
      : Vector<T>(r * c), 
232
233
	rows(r),
	cols(c)
234
    {}
235

Thomas Witkowski's avatar
Thomas Witkowski committed
236
    /// Changes the size of the matrix to newRows x newCols.
237
238
239
240
241
242
    inline void resize(int newRows, int newCols) {
      if ((newRows != rows) || (newCols != cols)) {
	Vector<T>::resize(newRows * newCols);
	rows = newRows;
	cols = newCols;
      }
243
    }
244

Thomas Witkowski's avatar
Thomas Witkowski committed
245
    /// Assignement operator.
246
    inline const Matrix<T>& operator=(const T& scal) {
247
      return static_cast<const Matrix<T>&>(Vector<T>::operator=(scal));
248
    }
249

Thomas Witkowski's avatar
Thomas Witkowski committed
250
    ///
251
    inline bool operator==(const Matrix<T>& rhs) const {
252
253
254
      if (rows != rhs.getNumRows()) return false;
      if (cols != rhs.getNumCols()) return false;
      return Vector<T>::operator == (rhs);
255
    }
256

Thomas Witkowski's avatar
Thomas Witkowski committed
257
    /// Comparison operator.
258
259
    inline bool operator!=(const Matrix<T>& rhs) const {
      return !(*this == rhs);
260
261
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
267
    /// Acces to i-th matrix row for constant matrices.
268
    inline const T *operator[](int i) const {
269
      return this->valArray + cols * i;
270
    }
271

Thomas Witkowski's avatar
Thomas Witkowski committed
272
    /// Returns \ref rows.
273
    inline int getNumRows() const { 
274
      return rows; 
275
    }
276

Thomas Witkowski's avatar
Thomas Witkowski committed
277
    /// Return \ref cols.
278
    inline int getNumCols() const { 
279
      return cols; 
Thomas Witkowski's avatar
Thomas Witkowski committed
280
    }
281

Thomas Witkowski's avatar
Thomas Witkowski committed
282
    /// Returns \ref rows.
283
    inline int getSize() const {
284
      return rows; 
Thomas Witkowski's avatar
Thomas Witkowski committed
285
    }
286

Thomas Witkowski's avatar
Thomas Witkowski committed
287
    /// Returns pointer after the last vector element.
288
289
    inline T *end() const { 
      return this->valArray + (cols * rows); 
290
    }
291
292

    void print() const {
293
      std::cout << this->rows << " x " << this->cols << " matrix: " << std::endl;
294
295
      for (int i = 0; i < rows; i++) {
	for (int j = 0; j < cols; j++) {
296
	  std::cout << this->valArray[i * cols + j] << " ";
297
	}
298
	std::cout << std::endl;
299
      }
300
301
302
    }

    void invert2x2();
303
304

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
305
    /// Number of matrix rows.
306
307
    int rows;

Thomas Witkowski's avatar
Thomas Witkowski committed
308
    /// Number of matrix columns.
309
310
311
    int cols;
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
312
  /// Matrix vector multiplication.
313
314
315
  template<typename T>
  inline const Vector<T>& mv(const Matrix<T>& m, const Vector<T>& v, Vector<T>& result)
  {
316
317
    TEST_EXIT_DBG(m.getNumCols() == v.getSize())("m and v not compatible\n");
    TEST_EXIT_DBG(v.getSize() == result.getSize())("wrong result size\n");
318
319
320
321
322

    T *resultIt, *mIt, *vIt;

    for (resultIt = result.begin(), mIt = m.begin(); 
	 resultIt != result.end(); 
323
324
325
326
	 ++resultIt) {
      *resultIt = 0;
      for (vIt = v.begin(); vIt != v.end(); ++vIt, ++mIt) {
	*resultIt += *mIt * *vIt;
327
      }
328
    }
329
330
331
332
    return result;
  };

  /** \brief
333
   * Computes x(Ay)^T, with x and y vectors, and A a matrix.
334
335
   */
  template<typename T>
336
  inline T xAy(const Vector<T>& x, const Matrix<T>& a, const Vector<T>& y) 
337
  {
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
    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) {
361
362
363
364
365
366
367
    return mv(m, v, v);
  };

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


  /** \brief
   * Scalar product.
   */
  template<typename T> 
378
  inline double operator*(const Vector<T>& v1, const Vector<T>& v2) 
379
  {
380
    double result = 0.0;
381
382
383
384

    T *v1It, *v2It;
    for (v1It = v1.begin(), v2It = v2.begin();
	 v1It != v1.end();
385
386
387
	 ++v1It, ++v2It) {
      result += *v1It * *v2It;
    }
388
389
390
391
392
393
394
395
396
397

    return result;
  };

  /** \brief
   * Vector addition.
   */
  template<typename T> 
  inline const Vector<T>& add(const Vector<T>& v1, const Vector<T>& v2, Vector<T>& result)
  {
398
399
    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");
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
    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)
  {
419
    TEST_EXIT_DBG(v.getSize() == result.getSize())("invalid size\n");
420

421
422
423
    T *vIt, *resultIt;
    for (vIt = v.begin(), resultIt = result.begin();
	 vIt != v.end();
424
425
426
	 ++vIt, ++resultIt) {
      *resultIt = scal * *vIt;
    }
427
428
429
430
431
432
433
434
435
436

    return result;
  };

  /** \brief
   * vector + scalar
   */
  template<typename T>
  inline const Vector<T>& add(const Vector<T>& v, const T& scal, Vector<T>& result)
  {
437
    TEST_EXIT_DBG(v.getSize() == result.getSize())("invalid size\n");
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
    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)
  {
457
    TEST_EXIT_DBG(x.getSize() == y.getSize())("invalid size\n");
458
459
460
    T *xIt, *yIt;
    for (xIt = x.begin(), yIt = y.begin();
	 xIt != x.end();
461
462
463
	 ++xIt, ++yIt) {
      *yIt += a * *xIt;
    }
464
465
466
467
468

    return y;
  };

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

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

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

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

  template<typename T>
492
  inline Vector<T> operator+(const Vector<T>& v1, const Vector<T>& v2) {
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
518
519
520
521
522
523
524
525
    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)
  {
526
527
    FUNCNAME("vectorProduct()");
    TEST_EXIT_DBG(Global::getGeo(WORLD) == 3)("DIM_OF_WORLD != 3\n");
528
529
530
531
532
533
534
535
    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