MatrixVector.h 12.1 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
70
71
72
73
    Vector(const Vector<T>& rhs) 
      : Serializable(),size(rhs.size)
    {
      valArray = new T[rhs.size];
      *this = rhs; // uses operator=()
74
    }
75

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

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

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

97
      return *this;
98
    }
99

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

106
      return *this;
107
    }
108

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

119
      return *this;
120
    }
121

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

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

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

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

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

161
      return true;
162
    }
163

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    T *resultIt, *mIt, *vIt;

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

    return result;
  }

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

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

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

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

    return result;
386
  }
387

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

    return result;
401
  }
402

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

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

    return result;
418
  }
419

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

    return result;
432
  }
433

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

    return y;
448
  }
449

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

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

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

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

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

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

  template<typename T>
492
493
  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
500
501
502
503

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
521
522
523
524

  /// Returns a vector with the FE spaces of each row in the matrix. Thus, the
  /// resulting vector may contain the same FE space multiple times.
  vector<const FiniteElemSpace*> getFeSpaces(Matrix<DOFMatrix*> &mat);
525
526
527
}

#endif // AMDIS_MATRIXVECTOR_H