MatrixVector.h 11.8 KB
Newer Older
1
2
3
4
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
5
// ==  http://www.amdis-fem.org                                              ==
6
7
// ==                                                                        ==
// ============================================================================
8
9
10
11
12
13
14
15
16
17
18
19
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


20
21
22
23
24
25
26
27

/** \file MatrixVector.h */

#ifndef AMDIS_MATRIXVECTOR_H
#define AMDIS_MATRIXVECTOR_H

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

namespace AMDiS {

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
54
    /// Change the size of the vector to newSize.
55
56
    inline void resize(int newSize) 
    {
57
      if (size != newSize) {
Thomas Witkowski's avatar
Thomas Witkowski committed
58
	if (valArray != NULL) 
59
60
61
62
	  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() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
76
77
78
79
      if (valArray != NULL) {
	delete [] valArray; 
	valArray = NULL;
      }
80
    }
81

Thomas Witkowski's avatar
Thomas Witkowski committed
82
    /// Assignement operator
83
84
    inline const Vector<T>& operator=(const Vector<T>& rhs) 
    {
85
86
87
      TEST_EXIT_DBG(rhs.size == size)
	("Invalid sizes %d != %d!\n", rhs.size, size);

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

94
      return *this;
95
    }
96

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

103
      return *this;
104
    }
105

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

116
      return *this;
117
    }
118

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

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

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

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

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

158
      return true;
159
    }
160

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    T *resultIt, *mIt, *vIt;

    for (resultIt = result.begin(), mIt = m.begin(); 
	 resultIt != result.end(); 
346
347
	 ++resultIt) {
      *resultIt = 0;
348
      for (vIt = v.begin(); vIt != v.end(); ++vIt, ++mIt)
349
350
	*resultIt += *mIt * *vIt;
    }
351
352
353
354

    return result;
  }

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

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

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

    T *v1It, *v2It;
    for (v1It = v1.begin(), v2It = v2.begin();
	 v1It != v1.end();
379
	 ++v1It, ++v2It)
380
      result += *v1It * *v2It;
381
382

    return result;
383
  }
384

385
  /// Vector addition.
386
387
388
  template<typename T> 
  inline const Vector<T>& add(const Vector<T>& v1, const Vector<T>& v2, Vector<T>& result)
  {
389
390
    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");
391
392
393
394
    T *v1It, *v2It, *resultIt;
    for (v1It = v1.begin(), v2It = v2.begin(), resultIt = result.begin();
	 v1It != v1.end();
	 ++v1It, ++v2It, ++resultIt)
395
      *resultIt = *v1It + *v2It;
396
397

    return result;
398
  }
399

400
  /// scalar * vector
401
402
403
404
405
  template<typename T>
  inline const Vector<T>& mult(const T& scal, 
			       const Vector<T>& v, 
			       Vector<T>& result)
  {
406
    TEST_EXIT_DBG(v.getSize() == result.getSize())("invalid size\n");
407

408
409
410
    T *vIt, *resultIt;
    for (vIt = v.begin(), resultIt = result.begin();
	 vIt != v.end();
411
	 ++vIt, ++resultIt) 
412
      *resultIt = scal * *vIt;
413
414

    return result;
415
  }
416

417
  /// vector + scalar
418
419
420
  template<typename T>
  inline const Vector<T>& add(const Vector<T>& v, const T& scal, Vector<T>& result)
  {
421
    TEST_EXIT_DBG(v.getSize() == result.getSize())("invalid size\n");
422
423
424
425
    T *vIt, *resultIt;
    for (vIt = v.begin(), resultIt = result.begin();
	 vIt != v.end();
	 ++vIt, ++resultIt)
426
      *resultIt = *vIt + scal;
427
428

    return result;
429
  }
430

431
  /// y = a * x + y.
432
433
434
435
436
  template<typename T>
  inline const Vector<T>& axpy(const T& a,
			       const Vector<T> &x,
			       Vector<T> &y)
  {
437
    TEST_EXIT_DBG(x.getSize() == y.getSize())("invalid size\n");
438
439
440
    T *xIt, *yIt;
    for (xIt = x.begin(), yIt = y.begin();
	 xIt != x.end();
441
	 ++xIt, ++yIt) 
442
      *yIt += a * *xIt;
443
444

    return y;
445
  }
446
447

  template<typename T>
448
449
  inline const Vector<T>& operator*=(Vector<T>& v, const T& scal) 
  {
450
    return mult(scal, v, v);
451
  }
452
453

  template<typename T>
454
455
  inline Vector<T> operator*(const Vector<T>& v, const T& scal) 
  {
456
457
458
    Vector<T> result = v;
    result *= scal;
    return result;
459
  }
460
461

  template<typename T>
462
463
  inline const Vector<T>& operator+(const Vector<T>& v1, const T& scal) 
  {
464
465
    Vector<T> result(v1.getSize());
    return add(v1, scal, result);
466
  }
467
468

  template<typename T>
469
470
  inline const Vector<T>& operator+=(Vector<T>& v1, const Vector<T>& v2) 
  {
471
    return add(v1, v2, v1);
472
  }
473
474

  template<typename T>
475
476
  inline Vector<T> operator+(const Vector<T>& v1, const Vector<T>& v2) 
  {
477
478
479
    Vector<T> result = v1;
    result += v2;
    return result;
480
  }
481
482

  template<typename T>
483
484
  const Vector<T>& operator-=(Vector<T>& v1, const Vector<T>& v2)
  {
485
    return axpy(-1.0, v2, v1);
486
  }
487
488

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

  template<typename T>
  inline double norm(const Vector<T> *v)
  {
    T *vIt;
    double result = 0;
501
    for (vIt = v->begin(); vIt != v->end(); ++vIt)
502
503
      result += *vIt * *vIt;
    return sqrt(result);
504
  }
505
506
507
508
509
510

  template<typename T>
  void vectorProduct(const Vector<T>& x, 
		     const Vector<T>& y, 
		     Vector<T>& z)
  {
511
512
    FUNCNAME("vectorProduct()");
    TEST_EXIT_DBG(Global::getGeo(WORLD) == 3)("DIM_OF_WORLD != 3\n");
513
514
515
    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];
516
  }
517
518
519
520

}

#endif // AMDIS_MATRIXVECTOR_H