MatrixVector.h 14.4 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
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
// ============================================================================
// ==                                                                        ==
// == 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;

  // ============================================================================
  // ===== class Vector =========================================================
  // ============================================================================

  /** \brief
   * Class for efficient vector operations of fixed size numerical vectors.
   */
  template<typename T>
  class Vector : public Serializable
  {
  public:
    MEMORY_MANAGED(Vector<T>);

    /** \brief
     * Constructor.
     */
    Vector(int i = 0) 
      : size(i) 
    {
      if (size == 0) 
	valArray = NULL;
      else
	valArray = new T[size];
57
    }
58

59
    inline bool used() const {
60
      return (valArray != NULL);
61
    }
62
63
64
65
66
67
68
69
70
71
72

    /** \brief
     * Change the size of the vector to newSize.
     */
    inline void resize(int newSize) {
      if (size != newSize) {
	if (valArray) 
	  delete [] valArray;
	valArray = new T[newSize];
	size = newSize;
      }
73
    }
74
75
76
77
78
79
80
81
82

    /** \brief
     * Copy constructor.
     */
    Vector(const Vector<T>& rhs) 
      : Serializable(),size(rhs.size)
    {
      valArray = new T[rhs.size];
      *this = rhs; // uses operator=()
83
    }
84
85
86
87
88
89

    /** \brief
     * Destructor.
     */
    virtual ~Vector() { 
      delete [] valArray; 
90
    }
91
92
93
94
95

    /** \brief
     * Assignement operator
     */
    inline const Vector<T>& operator=(const Vector<T>& rhs) {
96
      TEST_EXIT_DBG(rhs.size == size)("invalid size\n");
97
      T *rhsIt, *thisIt;
98
99
100
101
102
103
      for (rhsIt = rhs.begin(), thisIt = this->begin();
	   rhsIt != rhs.end();
	   ++rhsIt, ++thisIt) {
	*thisIt = *rhsIt;
      }

104
      return *this;
105
    }
106
107
108
109
110
111
112
113
114
115
116
117
118

    /** \brief
     * Assignement operator
     */
    inline const Vector<T>& operator=(const T& scal) {
      T *thisIt;
      for(thisIt = this->begin();
	  thisIt != this->end();
	  ++thisIt)
	{
	  *thisIt = scal;
	}
      return *this;
119
    }
120
121
122
123
124
125
126
127
128
129
130
131
132
133

    /** \brief
     * Assignement operator
     */
    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;
134
    }
135
136
137
138
139
140

    /** \brief
     * Sets all entries to scal.
     */
    inline const Vector<T>& set(const T& scal) {
      return *this = scal;
141
    }
142
143
144
145
146
147
148
149
150
151
152
153
154
155

    /** \brief
     * Sets all entries.
     */
    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;
156
    }
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179

    /** \brief
     * Sets all entries.
     */
    inline void fill(const T value) {
      for (T *thisIt = this->begin(); thisIt != this->end(); thisIt++) {
	*thisIt = value;
      }
    }

    /** \brief
     * Comparison operator.
     */
    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;
180
    }
181
182
183
184
185
186
187
188
189
190
191
192

    /** \brief
     * Comparison operator.
     */
    inline bool operator!=(const Vector<T>& rhs) const {
      return !(*this==rhs);
    }

    /** \brief
     * Access to the i-th vector element.
     */
    inline T& operator[](int i) {
193
      TEST_EXIT_DBG(i < size && i >= 0)("invalid index\n");
194
      return valArray[i];
195
    }
196
197
198
199
200

    /** \brief
     * Access to the i-th vector element for const vectors.
     */
    inline const T& operator[] (int i) const {
201
      TEST_EXIT_DBG(i < size && i >= 0)("invalid index\n");
202
      return valArray[i];
203
    }
204
205
206
207
208
209

    /** \brief
     * Returns pointer to the first vector element.
     */
    inline T *begin() const { 
      return valArray; 
210
    }
211
212
213
214
215
216

    /** \brief
     * Returns pointer after the last vector element.
     */
    inline T *end() const { 
      return valArray + size; 
217
    }
218
219
220
221
222
223

    /** \brief
     * Returns \ref size.
     */
    virtual int getSize() const { 
      return size; 
224
    }
225
226
227
228
229
230

    /** \brief
     * Returns \ref valArray as T-array
     */
    inline T *getValArray() { 
      return valArray; 
231
    }
232
233

    void print() const {
234
      std::cout << this->size << " vector: " << std::endl;
235
      for (int i = 0; i < size; i++) {
236
	std::cout << this->valArray[i] << " ";
237
      }
238
      std::cout << std::endl;
239
    }
240
241
242

    // ===== Serializable implementation =====
  
243
    void serialize(std::ostream &out) {
244
245
      out.write(reinterpret_cast<const char*>(&size), sizeof(int));
      out.write(reinterpret_cast<const char*>(valArray), size * sizeof(T));
246
    }
247

248
    void deserialize(std::istream &in) {
249
250
      in.read(reinterpret_cast<char*>(&size), sizeof(int));
      in.read(reinterpret_cast<char*>(valArray), size * sizeof(T));
251
    }
252

253
    std::string getTypeName() const { 
254
      return "Vector"; 
255
    }
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290

  protected:
    /** \brief
     * Size of the vector.
     */
    int size;

    /** \brief
     * Internal storage pointer.
     */
    T *valArray;
  };



  // ============================================================================
  // ===== class Matrix =========================================================
  // ============================================================================

  /** \brief
   * Class for efficient matrix operations of fixed size numerical matrices.
   */
  template<typename T>
  class Matrix : public Vector<T>
  {
  public:
    MEMORY_MANAGED(Matrix<T>);

    /** \brief
     * Constructor.
     */
    Matrix(int r, int c)
      : Vector<T>(r*c), 
	rows(r),
	cols(c)
291
    {}
292
293
294
295
296
297
298
299
300
301

    /** \brief
     * Changes the size of the matrix to newRows x newCols.
     */
    inline void resize(int newRows, int newCols) {
      if ((newRows != rows) || (newCols != cols)) {
	Vector<T>::resize(newRows * newCols);
	rows = newRows;
	cols = newCols;
      }
302
    }
303
304
305
306

    /** \brief
     * Assignement operator.
     */
307
    inline const Matrix<T>& operator=(const T& scal) {
308
      return static_cast<const Matrix<T>&>(Vector<T>::operator=(scal));
309
    }
310

311
    inline bool operator==(const Matrix<T>& rhs) const {
312
313
314
      if (rows != rhs.getNumRows()) return false;
      if (cols != rhs.getNumCols()) return false;
      return Vector<T>::operator == (rhs);
315
    }
316
317
318
319

    /** \brief
     * Comparison operator.
     */
320
321
    inline bool operator!=(const Matrix<T>& rhs) const {
      return !(*this == rhs);
322
323
324
325
326
    }

    /** \brief
     * Acces to i-th matrix row.
     */
327
    inline T *operator[](int i) {
328
      return this->valArray + cols * i;
329
    }
330
331
332
333

    /** \brief
     * Acces to i-th matrix row for constant matrices.
     */
334
    inline const T *operator[](int i) const {
335
      return this->valArray + cols * i;
336
    }
337
338
339
340

    /** \brief
     * Returns \ref rows.
     */
341
    inline int getNumRows() const { 
342
      return rows; 
343
    }
344
345
346
347

    /** \brief
     * Return \ref cols.
     */
348
    inline int getNumCols() const { 
349
350
351
352
353
354
      return cols; 
    };

    /** \brief
     * Returns \ref rows.
     */
355
    inline int getSize() const {
356
357
358
359
360
361
362
363
      return rows; 
    }; 

    /** \brief
     * Returns pointer after the last vector element.
     */
    inline T *end() const { 
      return this->valArray + (cols * rows); 
364
    }
365
366

    void print() const {
367
      std::cout << this->rows << " x " << this->cols << " matrix: " << std::endl;
368
369
      for (int i = 0; i < rows; i++) {
	for (int j = 0; j < cols; j++) {
370
	  std::cout << this->valArray[i * cols + j] << " ";
371
	}
372
	std::cout << std::endl;
373
      }
374
375
376
    }

    void invert2x2();
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395

  protected:
    /** \brief
     * Number of matrix rows.
     */
    int rows;

    /** \brief
     * Number of matrix columns.
     */
    int cols;
  };

  /** \brief
   * Matrix vector multiplication.
   */
  template<typename T>
  inline const Vector<T>& mv(const Matrix<T>& m, const Vector<T>& v, Vector<T>& result)
  {
396
397
    TEST_EXIT_DBG(m.getNumCols() == v.getSize())("m and v not compatible\n");
    TEST_EXIT_DBG(v.getSize() == result.getSize())("wrong result size\n");
398
399
400
401
402

    T *resultIt, *mIt, *vIt;

    for (resultIt = result.begin(), mIt = m.begin(); 
	 resultIt != result.end(); 
403
404
405
406
	 ++resultIt) {
      *resultIt = 0;
      for (vIt = v.begin(); vIt != v.end(); ++vIt, ++mIt) {
	*resultIt += *mIt * *vIt;
407
      }
408
    }
409
410
411
412
    return result;
  };

  /** \brief
413
   * Computes x(Ay)^T, with x and y vectors, and A a matrix.
414
415
   */
  template<typename T>
416
  inline T xAy(const Vector<T>& x, const Matrix<T>& a, const Vector<T>& y) 
417
  {
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
    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) {
441
442
443
444
445
446
447
    return mv(m, v, v);
  };

  /** \brief
   * Matrix vector multiplication.
   */
  template<typename T>
448
  inline Vector<T> operator*(const Matrix<T>& m, const Vector<T>& v) {
449
450
451
452
453
454
455
456
457
    Vector<T> result(v.getSize());
    return mv(m, v, result);
  };


  /** \brief
   * Scalar product.
   */
  template<typename T> 
458
  inline double operator*(const Vector<T>& v1, const Vector<T>& v2) 
459
  {
460
    double result = 0.0;
461
462
463
464

    T *v1It, *v2It;
    for (v1It = v1.begin(), v2It = v2.begin();
	 v1It != v1.end();
465
466
467
	 ++v1It, ++v2It) {
      result += *v1It * *v2It;
    }
468
469
470
471
472
473
474
475
476
477

    return result;
  };

  /** \brief
   * Vector addition.
   */
  template<typename T> 
  inline const Vector<T>& add(const Vector<T>& v1, const Vector<T>& v2, Vector<T>& result)
  {
478
479
    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");
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
    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)
  {
499
    TEST_EXIT_DBG(v.getSize() == result.getSize())("invalid size\n");
500

501
502
503
    T *vIt, *resultIt;
    for (vIt = v.begin(), resultIt = result.begin();
	 vIt != v.end();
504
505
506
	 ++vIt, ++resultIt) {
      *resultIt = scal * *vIt;
    }
507
508
509
510
511
512
513
514
515
516

    return result;
  };

  /** \brief
   * vector + scalar
   */
  template<typename T>
  inline const Vector<T>& add(const Vector<T>& v, const T& scal, Vector<T>& result)
  {
517
    TEST_EXIT_DBG(v.getSize() == result.getSize())("invalid size\n");
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
    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)
  {
537
    TEST_EXIT_DBG(x.getSize() == y.getSize())("invalid size\n");
538
539
540
    T *xIt, *yIt;
    for (xIt = x.begin(), yIt = y.begin();
	 xIt != x.end();
541
542
543
	 ++xIt, ++yIt) {
      *yIt += a * *xIt;
    }
544
545
546
547
548

    return y;
  };

  template<typename T>
549
  inline const Vector<T>& operator*=(Vector<T>& v, const T& scal) {
550
551
552
553
    return mult(scal, v, v);
  };

  template<typename T>
554
  inline Vector<T> operator*(const Vector<T>& v, const T& scal) {
555
556
557
558
559
560
    Vector<T> result = v;
    result *= scal;
    return result;
  };

  template<typename T>
561
  inline const Vector<T>& operator+(const Vector<T>& v1, const T& scal) {
562
563
564
565
566
    Vector<T> result(v1.getSize());
    return add(v1, scal, result);
  };

  template<typename T>
567
  inline const Vector<T>& operator+=(Vector<T>& v1, const Vector<T>& v2) {
568
569
570
571
    return add(v1, v2, v1);
  };

  template<typename T>
572
  inline Vector<T> operator+(const Vector<T>& v1, const Vector<T>& v2) {
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
    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)
  {
606
607
    FUNCNAME("vectorProduct()");
    TEST_EXIT_DBG(Global::getGeo(WORLD) == 3)("DIM_OF_WORLD != 3\n");
608
609
610
611
612
613
614
615
    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