VectorOperations.h 19.7 KB
Newer Older
1
2
3
4
5
6
/** \file VectorOperations.h */

#ifndef VECTOR_OPERATIONS_H
#define VECTOR_OPERATIONS_H

#include "AMDiS.h"
7
#include "ValueTypes.h"
8
9
10
11
12
13
14
15
16
17
18

#if HAVE_PARALLEL_DOMAIN_AMDIS
#include "parallel/StdMpi.h"
#endif

#include <ios>
#include <iostream>
#include <fstream>
#include <string>
#include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/itl/itl.hpp>
19
#include <boost/utility/enable_if.hpp>
20
#include <boost/type_traits.hpp>
21
22
23
24
25
26
27
28
29
30
31
32
33

using namespace AMDiS; 

template<typename T>
const WorldVector<T>& operator*=(WorldVector<T>& v1, const WorldVector<T>& v2)
{
  T* v1It, *v2It;
  for (v1It = v1.begin(), v2It = v2.begin();
	  v1It != v1.end(); 
	  v1It++, v2It++)
  *v1It *= *v2It;
  
  return v1;
34
}
35
36
37
38
39
40
41
42
43

// x ~= y
template<typename T>
class compareTol : public binary_function<T, T, bool> {
public:
	compareTol(double tol_) : tol(tol_) {}
	bool operator()(const T &x, const T &y) const { return sqrt((x-y)*(x-y))<tol; }
	double tol;
};
44

45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
template<typename T1, typename T2, int component>
struct comparePair : public binary_function<std::pair<T1, T2>, std::pair<T1, T2>, bool> 
{
  bool operator()(const std::pair<T1, T2> &x, const std::pair<T1, T2> &y) const 
  { 
    return (component==0?(x.first<y.first):(x.second<y.second));
  }
};

template<typename T1, typename T2>
struct compare0 : public binary_function<std::pair<T1, T2>, std::pair<T1, T2>, bool> 
{
  bool operator()(const std::pair<T1, T2> &x, const std::pair<T1, T2> &y) const 
  { 
    return x.first<y.first;
  }
};
62

63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
template<typename T1, typename T2>
struct compare1 : public binary_function<std::pair<T1, T2>, std::pair<T1, T2>, bool> 
{
  bool operator()(const std::pair<T1, T2> &x, const std::pair<T1, T2> &y) const 
  { 
    return x.second<y.second;
  }
};

class SortPair : public binary_function<std::pair<WorldVector<double>, double>, std::pair<WorldVector<double>, double>, bool> {
public:
	SortPair(int component_=1,int idx_=0) : component(component_),idx(idx_) {}
	bool operator()(const std::pair<WorldVector<double>, double> &x, 
		const std::pair<WorldVector<double>, double> &y) const { 
		return (component==1?(x.first)[idx]<(y.first)[idx]:x.second<y.second); }
	int component,idx;
};

namespace vector_operations {

83
  namespace traits {
84

85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
    /// General declaration, used to disable unsupported types
    template <typename Collection, class Enable = void>
    struct num_rows {};
    
    /// size implementation for STL vectors
    template <typename Value>
    struct num_rows< std::vector<Value> > 
    {
      typedef std::size_t   type;
      type operator()(const std::vector<Value>& v) { return v.size(); }
    };

    /// size implementation for (1D) arrays interpreted as vectors
    template <typename Value, unsigned Size>
    struct num_rows<Value[Size]>
    {
      typedef std::size_t   type;
      type operator()(const Value[Size]) { return Size; }
103
    };   
104
105
106
107
108
109
110

    /// size implementation for (2D and higher) arrays interpreted as matrices
    template <typename Value, unsigned Rows, unsigned Cols>
    struct num_rows<Value[Rows][Cols]>
    {
      typedef std::size_t   type;
      type operator()(const Value[Rows][Cols]) { return Rows; }
111
    }; 
112
    
113
    /// size implementation for AMDiS WorldVectors
114
115
116
117
118
    template <typename Value>
    struct num_rows< WorldVector<Value> > 
    {
      typedef std::size_t   type;
      type operator()(const WorldVector<Value>& v) { return static_cast<size_t>(v.getSize()); }
119
    };    
120
    
121
    /// size implementation for AMDiS WorldMatrices
122
123
124
125
126
    template <typename Value>
    struct num_rows< WorldMatrix<Value> > 
    {
      typedef std::size_t   type;
      type operator()(const WorldMatrix<Value>& v) { return static_cast<size_t>(v.getNumRows()); }
127
    };
128
    
129
    /// size implementation for arithmetic types (double, float, int, ...)
130
131
132
133
134
135
136
137
138
139
140
141
142
    template <typename Value>
    struct num_rows< Value, typename boost::enable_if< boost::is_arithmetic<Value> >::type > 
    {
      typedef std::size_t   type;
      type operator()(const Value& v) { return 1; }
    };
    
    //________________________________________________________________________________    
    
    /// General declaration, used to disable unsupported types
    template <typename Collection, class Enable = void>
    struct num_cols {};
    
143
    /// size implementation for AMDiS WorldMatrices
144
145
146
147
148
    template <typename Value>
    struct num_cols< WorldMatrix<Value> > 
    {
      typedef std::size_t   type;
      type operator()(const WorldMatrix<Value>& v) { return static_cast<size_t>(v.getNumCols()); }
149
    };
150
151
152
153
154
155
156
    
    //________________________________________________________________________________    
    
    /// General declaration, used to disable unsupported types
    template <typename Collection, class Enable = void>
    struct resize {};
    
157
    /// change_dim implementation for AMDiS WorldVectors
158
159
160
161
162
163
164
165
    template <typename Value>
    struct resize< WorldVector<Value> > 
    {
      typedef void   vector_void_type;
      void operator()(WorldVector<Value>& v, size_t r) { 
	TEST_EXIT(Global::getGeo(WORLD) == r)
	  ("WorldVectors can not be resized!\n");
      }
166
    };  
167
    
168
    /// change_dim implementation for STL vectors
169
170
171
172
173
174
175
    template <typename Value>
    struct resize< std::vector<Value> > 
    {
      typedef void   vector_void_type;
      void operator()(std::vector<Value>& v, size_t r) { 
	v.resize(r);
      }
176
    };  
177
    
178
    /// change_dim implementation for MTL4 vectors
179
180
181
182
183
184
185
186
187
188
189
    template <typename Value>
    struct resize< mtl::dense_vector<Value> > 
    {
      typedef void   vector_void_type;
      void operator()(mtl::dense_vector<Value>& v, size_t r) { 
	v.change_dim(r);
      }
    };
        
    // _________________________________________________________________________________
        
190
    /// change_dim implementation for AMDiS WorldMatrices
191
192
193
194
195
    template <typename Value>
    struct resize< WorldMatrix<Value> > 
    {
      typedef void   matrix_void_type;
      void operator()(WorldMatrix<Value>& v, size_t r, size_t c) { 
196
197
	size_t dow = static_cast<size_t>(Global::getGeo(WORLD));
	TEST_EXIT(dow == r && dow == c)
198
199
200
201
	  ("WorldMatrices can not be resized!\n");
      }
    };   
    
202
    /// change_dim implementation for MTL4 matrix
203
204
205
206
207
208
209
210
211
    template <typename Value, typename Param>
    struct resize< mtl::matrix::base_matrix<Value, Param> > 
    {
      typedef void   matrix_void_type;
      void operator()(mtl::matrix::base_matrix<Value, Param>& v, size_t r, size_t c) { 
	v.change_dim(r,c);
      }
    };
    
212
    /// change_dim implementation for MTL4 matrix
213
214
215
216
217
218
219
220
221
    template <typename Value>
    struct resize< mtl::matrix::dense2D<Value> > 
    {
      typedef void   matrix_void_type;
      void operator()(mtl::matrix::dense2D<Value>& v, size_t r, size_t c) { 
	v.change_dim(r,c);
      }
    };
    
222
    /// change_dim implementation for MTL4 matrix
223
224
225
226
227
228
229
230
    template <typename Value>
    struct resize< mtl::matrix::compressed2D<Value> > 
    {
      typedef void   matrix_void_type;
      void operator()(mtl::matrix::compressed2D<Value>& v, size_t r, size_t c) { 
	v.change_dim(r,c);
      }
    };
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
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
    
    // _________________________________________________________________________________
    
    
    /// General declaration, used to disable unsupported types
    template <typename Collection, class Enable = void>
    struct at {};
    
    template <typename Value>
    struct at< WorldMatrix<Value> > 
    {
      typedef Value   type;
      type& operator()(WorldMatrix<Value>& v, size_t r, size_t c) { 
	return v[r][c];
      }
      type const& operator()(const WorldMatrix<Value>& v, size_t r, size_t c) { 
	return v[r][c];
      }
    };
    
    template <typename Value, typename Param>
    struct at< mtl::matrix::base_matrix<Value, Param> > 
    {
      typedef Value   type;
      type& operator()(mtl::matrix::base_matrix<Value, Param>& v, size_t r, size_t c) { 
	return v(r,c);
      }
      type const& operator()(const mtl::matrix::base_matrix<Value, Param>& v, size_t r, size_t c) { 
	return v(r,c);
      }
    };
    
    template <typename Value>
    struct at< mtl::matrix::dense2D<Value> > 
    {
      typedef Value   type;
      type& operator()(mtl::matrix::dense2D<Value>& v, size_t r, size_t c) { 
	return v(r,c);
      }
      type const& operator()(const mtl::matrix::dense2D<Value>& v, size_t r, size_t c) { 
	return v(r,c);
      }
    };
    
    template <typename Value>
    struct at< mtl::matrix::compressed2D<Value> > 
    {
      typedef Value   type;
      type& operator()(mtl::matrix::compressed2D<Value>& v, size_t r, size_t c) { 
	return v(r,c);
      }
      type const& operator()(const mtl::matrix::compressed2D<Value>& v, size_t r, size_t c) { 
	return v(r,c);
      }
    };
286
287
288
289
290
  }

  /// num_rows function for non-MTL types (uses implicit enable_if)
  template <typename Collection>
  typename traits::num_rows<Collection>::type 
291
  num_rows(const Collection& c)
292
293
294
  {  
    return traits::num_rows<Collection>()(c);
  }
295
  
296
297
298
  /// num_cols function for non-MTL types (uses implicit enable_if)
  template <typename Collection>
  typename traits::num_cols<Collection>::type 
299
  num_cols(const Collection& c)
300
301
302
  {  
    return traits::num_cols<Collection>()(c);
  }
303
  
304
305
  /// resize function for vectors
  template <typename Collection>
306
307
  typename traits::resize<Collection>::vector_void_type 
  resize(Collection& c, size_t rows)
308
309
310
  {  
    traits::resize<Collection>()(c, rows);
  }
311
  
312
313
  /// resize function for matrices
  template <typename Collection>
314
315
  typename traits::resize<Collection>::matrix_void_type 
  resize(Collection& c, size_t rows, size_t cols)
316
317
318
  {  
    traits::resize<Collection>()(c, rows, cols);
  }
319
  
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
  /// at function for matrices
  template <typename Collection>
  typename traits::at<Collection>::type const&
  at(const Collection& c, size_t rows, size_t cols)
  {  
    return traits::at<Collection>()(c, rows, cols);
  }
  
  /// at function for matrices
  template <typename Collection>
  typename traits::at<Collection>::type&
  at(Collection& c, size_t rows, size_t cols)
  {  
    return traits::at<Collection>()(c, rows, cols);
  }
  
  /// 2-norm for vectors
  template<typename Vector>
  typename boost::enable_if< is_vector<Vector>, typename ValueType<Vector>::type >::type
  norm(const Vector &b)
340
  {
341
342
343
    typename ValueType<Vector>::type value;
    nullify(value);
    for (size_t i = 0; i < num_rows(b); ++i)
344
345
      value += sqr(b[i]);
    return sqrt(value);
346
  }
347
  
348
349
  /// 2-norm for MTL4 vectors
  template<typename T>
350
  T norm(const mtl::dense_vector<T> &b) { return two_norm(b); }
351
352
  
  /// 1-norm for matrices
353
  template<typename Matrix>
354
355
  typename boost::enable_if< is_matrix<Matrix>, double >::type
  matrix_norm(const Matrix &M)
356
357
358
359
360
361
362
  {
    typename Matrix::value_type asum;
    asum = 0.0;
    for (size_t i = 0; i < num_rows(M); ++i) {
      for (size_t j = 0; j < num_cols(M); ++j)
	asum += abs(M[i][j]);
    }
363
    return static_cast<double>(asum);
364
  }
365
366
  
  
367
368
369
370
  /// 2-norm for vectors with eps-pertubation
  template<typename Vector>
  typename boost::enable_if< is_vector<Vector>, double >::type
  norm(const Vector &b, const double eps)
371
  {
372
    return std::max(static_cast<double>(norm(b)), eps);
373
  }
374

375
  /// trace of a matrix
376
  template<typename Matrix>
377
  typename boost::enable_if< is_matrix<Matrix>, double >::type
378
  trace(Matrix &M)
379
  {
380
    TEST_EXIT_DBG(num_rows(M) == num_cols(M))
381
      ("Matrix dimensions must be equal!\n");
382
    typename Matrix::value_type sum;
383
    nullify(sum);
384
    for (size_t i = 0; i < num_rows(M); ++i)
385
386
      sum += M[i][i];
    return static_cast<double>(sum);
387
  }
388
  
389
390
391
392
  /// equivalent to matlab .* for vectors
  template<typename Vector1, typename Vector2, typename ResultVector>
  typename boost::enable_if< boost::mpl::and_< is_vector<Vector1>, is_vector<Vector2> >, ResultVector >::type
  mult(const Vector1& v1, const Vector2& v2)
393
  {
394
395
396
397
398
    TEST_EXIT_DBG(num_rows(v1)==num_rows(v2))
      ("Vectors must have the same size!\n");
    ResultVector result(num_rows(v1));
    for (size_t i = 0; i < num_rows(v1); ++i)
      result[i] = v1[i] * v2[i];
399
    return result;
400
  }
401
  
402
403
404
  template<typename Matrix, typename Vector, typename ResultVector>
  typename boost::enable_if< boost::mpl::and_< is_matrix<Matrix>, is_vector<Vector> >, ResultVector >::type
  mv(const Matrix &A, const Vector &x)
405
  {
406
407
408
409
    TEST_EXIT_DBG(num_cols(A)==num_rows(x))
      ("Matrix and vector dimension must fit together!");
    ResultVector result(num_rows(A));
    for (size_t i = 0; i < num_rows(A); ++i) {
410
      result[i] = 0.0;
411
412
      for (size_t j = 0; j < num_cols(A); ++j)
	result[i] += at(A,i,j) * x[j];
413
414
    }
    return result;
415
  }
416
  
417
418
419
420
  template<typename Matrix, typename Vector1, typename Vector2, typename ResultVector>
  typename boost::enable_if< boost::mpl::and_< is_matrix<Matrix>, is_vector<Vector1>, is_vector<Vector2> >, 
			     ResultVector >::type
  Axpy(const Matrix &A, const Vector1 &x, const Vector2 &y)
421
  {
422
423
    TEST_EXIT_DBG(num_cols(A)==num_row(x) && num_rows(A)==num_rows(y))
      ("Matrix and vector dimensions must fit together!");
424
    ResultVector result;
425
426
    resize(result, num_rows(A));
    for (size_t i = 0; i < num_rows(A); ++i) {
427
      result[i] = y[i];
428
429
      for (size_t j = 0; j < num_cols(A); ++j)
	result[i] += at(A,i,j) * x[j];
430
431
    }
    return result;
432
  }
433

434
  // copy v -> w
435
436
437
  template<typename VectorIn, typename VectorOut>
  typename boost::enable_if< boost::mpl::and_< is_vector<VectorIn>, is_vector<VectorOut> >, void >::type
  copy(const VectorIn &v, VectorOut &w)
438
  {
439
440
    for (size_t i = 0; i < std::min(num_rows(v), num_rows(w)); i++)
      w[i] = v[i];
441
  }
442
443
444
445
446
447
448
449
450
451
452
  
  template<class MatrixIn, class MatrixOut>
  typename boost::enable_if< boost::mpl::and_< is_matrix<MatrixIn>, is_matrix<MatrixOut> >, void >::type
  fillMatrix(const MatrixIn &m_in, MatrixOut &m)
  {
    resize(m, num_rows(m_in), num_cols(m_in));
    for (size_t i = 0; i < num_rows(m_in); ++i) {
      for (size_t j = 0; j < num_cols(m_in); ++j) {
	at(m,i,j) = at(m_in,i,j);
      }
    }
453
  }
454
455
456
457
458
459
  
  template<class MatrixIn, class MatrixOut>
  typename boost::enable_if< boost::mpl::and_< is_matrix<MatrixIn>, is_matrix<MatrixOut> >, void >::type
  copy(const MatrixIn &m_in, MatrixOut &m)
  {
    fillMatrix(m_in, m);
460
  }
461
  
462
463
  // v3:=[v1, v2]
  template<typename Vector1, typename Vector2, typename VectorOut>
464
465
466
  typename boost::enable_if< boost::mpl::and_< is_vector<Vector1>, is_vector<Vector2>, is_vector<VectorOut> >, 
			     void >::type
  merge(Vector1 &v1, Vector2 &v2, VectorOut &v3)
467
  {
468
469
470
471
472
    resize(v3, num_rows(v1) + num_rows(v2));
    for (size_t i = 0; i < num_rows(v1); i++)
      v3[i] = v1[i];
    for (size_t j = 0; j < num_rows(v2); j++)
      v3[j + num_rows(v1)] = v2[j];
473
  }
474
  
475
  template<typename Vector>
476
  void getMin(const Vector &v, typename ValueType<Vector>::type &minVal, size_t &minIdx)
477
  {
478
479
480
481
482
483
484
485
486
    typedef typename ValueType<Vector>::type T;
    TEST_EXIT(num_rows(v) > 0)("getMin of empty vector!\n");
    
    minVal = v[0];
    minIdx = 0;
    for (size_t i = 1; i < num_rows(v); i++) {
      if (v[i] < minVal) {
	minVal = v[i];
	minIdx = i;
487
488
      }
    }
489
  }
490
  
491
  template<typename Vector>
492
  void getMax(const Vector &v, typename ValueType<Vector>::type &maxVal, size_t &maxIdx)
493
  {
494
495
496
497
498
499
500
501
502
    typedef typename ValueType<Vector>::type T;
    TEST_EXIT(num_rows(v) > 0)("getMax of empty vector!\n");
    
    maxVal = v[0];
    maxIdx = 0;
    for (size_t i = 1; i < num_rows(v); i++) {
      if (v[i] > maxVal) {
	maxVal = v[i];
	maxIdx = i;
503
504
      }
    }
505
  }
506
507

  template<typename Vector>
508
  typename ValueType<Vector>::type
509
  min(const Vector &vec)
510
  {
511
512
513
514
    typename ValueType<Vector>::type minVal;
    size_t minIdx;
    getMin(vec, minVal, minIdx);
    return minVal;
515
  }
516

517
  template<typename Vector>
518
  typename ValueType<Vector>::type
519
  max(const Vector &vec)
520
  {
521
522
523
    typename ValueType<Vector>::type maxVal;
    size_t maxIdx;
    getMax(vec, maxVal, maxIdx);
524
    return maxVal;
525
  }
526
527
  
  template<typename Vector>
528
  typename ValueType<Vector>::type
529
  absmin(const Vector &v)
530
  {
531
532
533
534
535
536
537
    typedef typename ValueType<Vector>::type T;
    TEST_EXIT(num_rows(v) > 0)("absmin of empty vector!\n");
    
    T minVal = abs(v[0]);
    for(size_t i = 1; i < num_rows(v); i++) {
      if (std::abs(v[i]) < minVal) {
	minVal = std::abs(v[i]);
538
539
540
      }
    }
    return minVal;
541
  }
542
  
543
  template<typename Vector>
544
  typename ValueType<Vector>::type
545
  absmax(const Vector &v)
546
  {
547
548
549
550
551
552
553
    typedef typename ValueType<Vector>::type T;
    TEST_EXIT(num_rows(v) > 0)("absmax of empty vector!\n");
    
    T maxVal = abs(v[0]);
    for(size_t i = 1; i < num_rows(v); i++) {
      if (std::abs(v[i]) > maxVal) {
	maxVal = std::abs(v[i]);
554
555
556
      }
    }
    return maxVal;
557
  }
558
  
559
  template<typename Vector>
560
  typename ValueType<Vector>::type
561
  sum(const Vector &v)
562
  {
563
564
565
566
567
568
569
    typedef typename ValueType<Vector>::type T;
    size_t N = num_rows(v);
    T value; nullify(value);
    if (N > 0) {
      value = v[0];
      for (size_t i = 1; i < N; ++i)
	value += v[i];
570
571
    }
    return value;
572
  }
573
  
574
575
576
  template<typename Vector>
  inline typename ValueType<Vector>::type
  mean(const Vector &v)
577
  {
578
579
    double N = static_cast<double>(num_rows(v));
    return sum(v) * (1.0 / N);
580
  }
581
  
582
583
584
585
  template<typename Vector>
  void getMean(const Vector &v, typename ValueType<Vector>::type &meanValue)
  {
    meanValue = mean(v);
586
  }
587
588
589
590
591
592
593
  
  template<typename T1, typename T2, typename Compare>
  void sort(std::vector<T1> &vec1, std::vector<T2> &vec2, Compare &comp)
  {
    std::vector<pair<T1,T2> > order(vec1.size());
    
    // create vector of pairs
594
    for (size_t i=0; i<vec1.size(); ++i)
595
596
597
598
599
600
      order[i] = make_pair(vec1[i], vec2[i]);
    
    // sort vector of pairs
    sort(order.begin(), order.end(), comp);
    
    // copy sorted pairs back to vectors
601
    for (size_t i=0; i<vec1.size(); ++i) {
602
603
604
      vec1[i] = order[i].first;
      vec2[i] = order[i].second;
    }
605
  }
606
607
608
609
610
  
  template<typename T1, typename T2, typename Compare>
  void sort(std::vector<std::pair<T1,T2> > &vec, Compare &comp)
  {
    std::sort(vec.begin(), vec.end(), comp);
611
  }
612
613
614
615
616
617
618
619
620
621
622
623
624
625
  
  inline void unique(std::vector<WorldVector<double> > &vec, double tol, std::vector<unsigned> &ind)
  {
    compareTol<WorldVector<double> > comp(tol);
    unsigned newVec=0;
    for(unsigned i=0; i<vec.size();++i) {
      bool inNew=false;
      for(unsigned j=0;j<newVec;++j) {
	inNew=inNew||comp(vec[i],vec[j]);
	if(inNew) break;
      }
      if(!inNew) {vec[newVec]=vec[i]; newVec++; ind.push_back(i);}
    }
    vec.erase(vec.begin()+newVec,vec.end());
626
  }
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
  
  template<typename T1, typename T2>
  void unique(std::vector<std::pair<T1,T2> > &vec, double tol, int uniqueBy=1)
  {
//     comparePair comp = comparePair(tol,uniqueBy);
//     unsigned newVec=0;
//     for(unsigned i=0; i<vec.size();++i) {
//       bool inNew=false;
//       for(unsigned j=0;j<newVec;++j) {
// 	inNew=inNew||comp(vec[i],vec[j]);
// 	if(inNew) break;
//       }
//       if(!inNew) {vec[newVec]=vec[i]; newVec++;}
//     }
//     vec.erase(vec.begin()+newVec,vec.end());
642
  }
643
  
644
645
  template<typename T>
  inline void filter(std::vector<T> &x, const std::vector<bool> &filter)
646
  {
647
648
649
650
651
652
    TEST_EXIT_DBG(num_rows(x) == num_rows(filter))
      ("number of rows of x and filter must be equal!\n");
    size_t j = 0;
    for (size_t i = 0; i < num_rows(x); ++i) {
      if (filter[i]) {
	x[j] = x[i];
653
654
655
656
	j++;
      }
    }
    x.erase(x.begin()+j,x.end());
657
  }
658
  
659
660
  template<typename T>
  void filter(std::vector<T> &x, const std::vector<int> &filter)
661
  {
662
663
664
665
    std::vector<T> x_(num_rows(filter));
    for (size_t i = 0; i < num_rows(filter); ++i)
      x_[i] = x[filter[i]];
    swap(x, x_);
666
667
  }
 
668
669
670
  template<typename Vector1, typename Vector2>
  typename ProductType<typename ValueType<Vector1>::type, typename ValueType<Vector2>::type>::type
  dot(const Vector1 &vec1, const Vector2 &vec2)
671
  {
672
673
674
675
676
677
678
679
    typedef typename ProductType<typename ValueType<Vector1>::type, typename ValueType<Vector2>::type>::type value_type;
    value_type value; nullify(value);
    
    size_t N = num_rows(vec1);
    if (N > 0) {
      value = vec1[0] * vec2[0];
      for (size_t i = 1; i < N; ++i)
	value += vec1[i] * vec2[i];
680
    }
681
    return value;
682
  }
683
684
685
}

#endif // VECTOR_OPERATIONS_H