MatrixVector.h 11.9 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:
38
39
40
    
    typedef T value_type;
    
Thomas Witkowski's avatar
Thomas Witkowski committed
41
    /// Constructor.
42
    Vector(int i = 0) 
Thomas Witkowski's avatar
Thomas Witkowski committed
43
44
      : size(i),
	valArray(NULL)
45
46
47
48
49
    {
      if (size == 0) 
	valArray = NULL;
      else
	valArray = new T[size];
50
    }
51

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
77
    /// Destructor.
78
79
    virtual ~Vector() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
80
81
82
83
      if (valArray != NULL) {
	delete [] valArray; 
	valArray = NULL;
      }
84
    }
85

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

92
      T *rhsIt, *thisIt;
93
94
      for (rhsIt = rhs.begin(), thisIt = this->begin();
	   rhsIt != rhs.end();
Thomas Witkowski's avatar
Thomas Witkowski committed
95
96
	   ++rhsIt, ++thisIt)
	*thisIt = *rhsIt;      
97

98
      return *this;
99
    }
100

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

107
      return *this;
108
    }
109

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

120
      return *this;
121
    }
122

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

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

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

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

155
      T *rhsIt, *thisIt;
156
157
158
159
160
161
      for (rhsIt = rhs.begin(), thisIt = this->begin();
	   rhsIt != rhs.end();
	   ++rhsIt, ++thisIt)
	if(*thisIt != *rhsIt) 
	  return false;

162
      return true;
163
    }
164

Thomas Witkowski's avatar
Thomas Witkowski committed
165
    /// Comparison operator.
166
167
    inline bool operator!=(const Vector<T>& rhs) const 
    {
168
169
170
      return !(*this==rhs);
    }

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
185
    /// Returns pointer to the first vector element.
186
187
    inline T *begin() const 
    { 
188
      return valArray; 
189
    }
190

Thomas Witkowski's avatar
Thomas Witkowski committed
191
    /// Returns pointer after the last vector element.
192
193
    inline T *end() const 
    { 
194
      return valArray + size; 
195
    }
196

Thomas Witkowski's avatar
Thomas Witkowski committed
197
    /// Returns \ref size.
198
199
    virtual int getSize() const 
    { 
200
      return size; 
201
    }
202

Thomas Witkowski's avatar
Thomas Witkowski committed
203
    /// Returns \ref valArray as T-array
204
205
    inline T *getValArray() 
    { 
206
      return valArray; 
207
    }
208

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

217
218
    void serialize(std::ostream &out) 
    {
219
220
      out.write(reinterpret_cast<const char*>(&size), sizeof(int));
      out.write(reinterpret_cast<const char*>(valArray), size * sizeof(T));
221
    }
222

223
224
    void deserialize(std::istream &in) 
    {
225
226
      in.read(reinterpret_cast<char*>(&size), sizeof(int));
      in.read(reinterpret_cast<char*>(valArray), size * sizeof(T));
227
    }
228

229
230
    std::string getTypeName() const 
    { 
231
      return "Vector"; 
232
    }
233
234

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
235
    /// Size of the vector.
236
237
    int size;

Thomas Witkowski's avatar
Thomas Witkowski committed
238
    /// Internal storage pointer.
239
240
241
242
    T *valArray;
  };


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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
265
    /// Assignement operator.
266
267
    inline const Matrix<T>& operator=(const T& scal) 
    {
268
      return static_cast<const Matrix<T>&>(Vector<T>::operator=(scal));
269
    }
270

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

Thomas Witkowski's avatar
Thomas Witkowski committed
279
    /// Comparison operator.
280
281
    inline bool operator!=(const Matrix<T>& rhs) const 
    {
282
      return !(*this == rhs);
283
284
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
291
    /// Access to i-th matrix row for constant matrices.
292
293
    inline const T *operator[](int i) const 
    {
294
      return this->valArray + cols * i;
295
    }
296

Thomas Witkowski's avatar
Thomas Witkowski committed
297
    /// Returns \ref rows.
298
299
    inline int getNumRows() const 
    { 
300
      return rows; 
301
    }
302

Thomas Witkowski's avatar
Thomas Witkowski committed
303
    /// Return \ref cols.
304
305
    inline int getNumCols() const 
    { 
306
      return cols; 
Thomas Witkowski's avatar
Thomas Witkowski committed
307
    }
308

Thomas Witkowski's avatar
Thomas Witkowski committed
309
    /// Returns \ref rows.
310
311
    inline int getSize() const 
    {
312
      return rows; 
Thomas Witkowski's avatar
Thomas Witkowski committed
313
    }
314

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

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

331
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
332
    /// Number of matrix rows.
333
334
    int rows;

Thomas Witkowski's avatar
Thomas Witkowski committed
335
    /// Number of matrix columns.
336
337
338
    int cols;
  };

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

    T *resultIt, *mIt, *vIt;

    for (resultIt = result.begin(), mIt = m.begin(); 
	 resultIt != result.end(); 
350
351
	 ++resultIt) {
      *resultIt = 0;
352
      for (vIt = v.begin(); vIt != v.end(); ++vIt, ++mIt)
353
354
	*resultIt += *mIt * *vIt;
    }
355
356
357
358

    return result;
  }

359
  /// Matrix vector multiplication.
360
  template<typename T>
361
362
  inline const Vector<T>& operator*=(const Vector<T>& v, const Matrix<T>& m) 
  {
363
    return mv(m, v, v);
364
  }
365

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

374
  /// Scalar product.
375
  template<typename T> 
376
  inline double operator*(const Vector<T>& v1, const Vector<T>& v2) 
377
  {
378
    double result = 0.0;
379
380
381
382

    T *v1It, *v2It;
    for (v1It = v1.begin(), v2It = v2.begin();
	 v1It != v1.end();
383
	 ++v1It, ++v2It)
384
      result += *v1It * *v2It;
385
386

    return result;
387
  }
388

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

    return result;
402
  }
403

404
  /// scalar * vector
405
406
  template<typename T, typename S>
  inline const Vector<T>& mult(const S& scal,
407
408
409
			       const Vector<T>& v, 
			       Vector<T>& result)
  {
410
    TEST_EXIT_DBG(v.getSize() == result.getSize())("invalid size\n");
411

412
413
414
    T *vIt, *resultIt;
    for (vIt = v.begin(), resultIt = result.begin();
	 vIt != v.end();
415
	 ++vIt, ++resultIt) 
416
      *resultIt = scal * *vIt;
417
418

    return result;
419
  }
420

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

    return result;
433
  }
434

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

    return y;
449
  }
450

451
452
  template<typename T, typename S>
  inline const Vector<T>& operator*=(Vector<T>& v, const S& scal)
453
  {
454
    return mult(scal, v, v);
455
  }
456
457

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

  template<typename T>
466
467
  inline const Vector<T>& operator+(const Vector<T>& v1, const T& scal) 
  {
468
469
    Vector<T> result(v1.getSize());
    return add(v1, scal, result);
470
  }
471
472

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

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

  template<typename T>
487
488
  const Vector<T>& operator-=(Vector<T>& v1, const Vector<T>& v2)
  {
489
    return axpy(-1.0, v2, v1);
490
  }
491
492

  template<typename T>
493
494
  Vector<T> operator-(const Vector<T>& v1, const Vector<T>& v2)
  {
495
496
497
    Vector<T> result = v1;
    result -= v2;
    return result;
498
  }
499
500
501
502
503
504

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

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

#endif // AMDIS_MATRIXVECTOR_H