OpenMEEG
vector.h
Go to the documentation of this file.
1 // Project Name: OpenMEEG (http://openmeeg.github.io)
2 // © INRIA and ENPC under the French open source license CeCILL-B.
3 // See full copyright notice in the file LICENSE.txt
4 // If you make a copy of this file, you must either:
5 // - provide also LICENSE.txt and modify this header to refer to it.
6 // - replace this header by the LICENSE.txt content.
7 
8 #pragma once
9 
10 #include <OMassert.H>
11 #include <cstdlib>
12 #include <string>
13 
14 #include <OpenMEEGMathsConfig.h>
15 #include <linop.h>
16 #include <MathsIO.H>
17 
18 namespace OpenMEEG {
19 
20  class Matrix;
21  class SymMatrix;
22 
23  class OPENMEEGMATHS_EXPORT Vector: public LinOp {
24 
25  LinOpValue value;
26 
27  public:
28 
29  Vector(): LinOp(0,1,FULL,1),value() { }
30 
31  Vector(const Dimension N): LinOp(N,1,FULL,1),value(N) { }
32  Vector(const Vector& A,const DeepCopy): LinOp(A.nlin(),1,FULL,1),value(A.size(),A.data()) { }
33 
34  explicit Vector(const Matrix& A);
35  explicit Vector(const SymMatrix& A);
36 
37  void alloc_data() { value = LinOpValue(size()); }
38  void reference_data(const double* array) { value = LinOpValue(size(),array); }
39 
40  size_t size() const { return nlin(); }
41 
42  bool empty() const { return value.empty(); }
43 
44  double* data() const { return value.get(); }
45 
46  double operator()(const Index i) const {
47  om_assert(i<nlin());
48  return value[i];
49  }
50 
51  double& operator()(const Index i) {
52  om_assert(i<nlin());
53  return value[i];
54  }
55 
56  Vector subvect(const Index istart,const Index isize) const;
57 
58  Vector operator+(const Vector& v) const;
59  Vector operator-(const Vector& v) const;
60 
61  Vector operator-() const {
62  // No Blas implementation ?
63  Vector res(nlin());
64  for (Index i=0; i<nlin(); i++ )
65  res.data()[i] = -data()[i];
66  return res;
67  }
68 
69  inline void operator+=(const Vector& v);
70  inline void operator-=(const Vector& v);
71  inline void operator*=(const double x);
72  void operator/=(const double x) { (*this) *= (1.0/x); }
73  Vector operator+(const double i) const;
74  Vector operator-(const double i) const;
75  inline Vector operator*(const double x) const;
76  Vector operator/(const double x) const { return (*this)*(1.0/x); }
77  inline double operator*(const Vector& v) const;
78  Vector operator*(const Matrix& m) const;
79 
80  Vector kmult(const Vector& x) const;
81  // Vector conv(const Vector& v) const;
82  // Vector conv_trunc(const Vector& v) const;
83  Matrix outer_product(const Vector& v) const;
84 
85  double norm() const;
86  double sum() const;
87  double mean() const { return sum()/size(); }
88 
89  void set(const double x);
90 
91  void save(const char *filename) const;
92  void load(const char *filename);
93 
94  void save(const std::string& s) const { save(s.c_str()); }
95  void load(const std::string& s) { load(s.c_str()); }
96 
97  void info() const;
98 
99  friend class SymMatrix;
100  friend class Matrix;
101  };
102 
103  OPENMEEGMATHS_EXPORT Vector operator*(const double d,const Vector& v);
104 
105  OPENMEEGMATHS_EXPORT std::ostream& operator<<(std::ostream& f,const Vector& M);
106  OPENMEEGMATHS_EXPORT std::istream& operator>>(std::istream& f,Vector& M);
107 
108  inline Vector Vector::subvect(const Index istart,const Index isize) const {
109  om_assert (istart+isize<=nlin());
110  Vector a(isize);
111  for (Index i=0; i<isize; ++i)
112  a(i) = (*this)(istart+i);
113  return a;
114  }
115 
116  inline Vector Vector::operator+(const Vector& v) const {
117  om_assert(nlin()==v.nlin());
118  Vector p(*this,DEEP_COPY);
119  #ifdef HAVE_BLAS
120  BLAS(daxpy,DAXPY)(sizet_to_int(nlin()),1,v.data(),1,p.data(),1);
121  #else
122  for (Index i=0; i<nlin(); ++i)
123  p.data()[i] = data()[i]+v.data()[i];
124  #endif
125  return p;
126  }
127 
128  inline Vector Vector::operator-(const Vector& v) const {
129  om_assert(nlin()==v.nlin());
130  Vector p(*this,DEEP_COPY);
131  #ifdef HAVE_BLAS
132  BLAS(daxpy,DAXPY)(sizet_to_int(nlin()),-1,v.data(),1,p.data(),1);
133  #else
134  for (Index i=0; i<nlin(); ++i)
135  p.data()[i] = data()[i]-v.data()[i];
136  #endif
137  return p;
138  }
139 
140  inline void Vector::operator+=(const Vector& v) {
141  om_assert(nlin()==v.nlin());
142  #ifdef HAVE_BLAS
143  BLAS(daxpy,DAXPY)(sizet_to_int(nlin()),1,v.data(),1,data(),1);
144  #else
145  for (Index i=0; i<nlin(); ++i)
146  data()[i] += v.data()[i];
147  #endif
148  }
149 
150  inline void Vector::operator-=(const Vector& v) {
151  om_assert(nlin()==v.nlin());
152  #ifdef HAVE_BLAS
153  BLAS(daxpy,DAXPY)(sizet_to_int(nlin()),-1,v.data(),1,data(),1);
154  #else
155  for (Index i=0; i<nlin(); ++i)
156  data()[i] -= v.data()[i];
157  #endif
158  }
159 
160  inline double Vector::operator*(const Vector& v) const {
161  om_assert(nlin()==v.nlin());
162  #ifdef HAVE_BLAS
163  return BLAS(ddot,DDOT)(sizet_to_int(nlin()),data(),1,v.data(),1);
164  #else
165  double s=0;
166  for (Index i=0; i<nlin(); ++i)
167  s += data()[i]*v.data()[i];
168  return s;
169  #endif
170  }
171 
172  inline Vector Vector::operator*(const double x) const {
173  #ifdef HAVE_BLAS
174  Vector p(*this,DEEP_COPY);
175  BLAS(dscal,DSCAL)(sizet_to_int(nlin()),x,p.data(),1);
176  #else
177  Vector p(nlin());
178  for (Index i=0; i<nlin(); ++i)
179  p.data()[i] = x*data()[i];
180  #endif
181  return p;
182  }
183 
184  inline void Vector::operator*=(const double x) {
185  #ifdef HAVE_BLAS
186  BLAS(dscal,DSCAL)(sizet_to_int(nlin()),x,data(),1);
187  #else
188  for (Index i=0; i<nlin(); ++i)
189  data()[i] *= x;
190  #endif
191  }
192 
193  inline double Vector::norm() const {
194  #ifdef HAVE_BLAS
195  return BLAS(dnrm2,DNRM2)(sizet_to_int(nlin()),data(),1);
196  #else
197  throw OpenMEEG::maths::LinearAlgebraError("'Vector::norm' not implemented, requires BLAS");
198  #endif
199  }
200 
201  // inline Vector Vector::conv(const Vector& v) const {
202  // if (v.nlin()<nlin())
203  // return v.conv(*this);
204  //
205  // Vector p(nlin()+v.nlin()-1);
206  // p.set(0);
207  // for (Index i=0; i<v.nlin(); ++i) {
208  // #ifdef HAVE_BLAS
209  // BLAS(daxpy,DAXPY)(sizet_to_int(nlin()),v(i),data(),1,p.data()+i,1);
210  // #else
211  // for (Index j=0; j<nlin(); ++j)
212  // p(i+j) += v(i)*data()[j];
213  // #endif
214  // }
215  // return p;
216  // }
217  //
218  // inline Vector Vector::conv_trunc(const Vector& v) const {
219  // Vector p(v.nlin());
220  // p.set(0);
221  // for (Index i=0; i<v.nlin(); ++i)
222  // {
223  // const Index m = std::min(nlin(),v.nlin()-i);
224  // #ifdef HAVE_BLAS
225  // BLAS(daxpy,DAXPY)((int)m,v(i),data(),1,p.data()+i,1);
226  // #else
227  // for (Index j=0; j<m; ++j)
228  // p(i+j) += v(i)*(*this)(j);
229  // #endif
230  // }
231  // return p;
232  // }
233 
234  // Operators.
235 
236  inline Vector operator*(const double d,const Vector& v) { return v*d; }
237 }
Dimension nlin() const
Definition: linop.h:48
Matrix class Matrix class.
Definition: matrix.h:28
Matrix outer_product(const Vector &v) const
void alloc_data()
Definition: vector.h:37
Vector operator/(const double x) const
Definition: vector.h:76
void operator*=(const double x)
Definition: vector.h:184
Vector operator*(const Matrix &m) const
void save(const std::string &s) const
Definition: vector.h:94
double operator()(const Index i) const
Definition: vector.h:46
void set(const double x)
Vector operator+(const Vector &v) const
Definition: vector.h:116
size_t size() const
Definition: vector.h:40
Vector kmult(const Vector &x) const
Vector operator-(const double i) const
void reference_data(const double *array)
Definition: vector.h:38
void operator+=(const Vector &v)
Definition: vector.h:140
void operator-=(const Vector &v)
Definition: vector.h:150
Vector operator+(const double i) const
Vector(const Vector &A, const DeepCopy)
Definition: vector.h:32
double norm() const
Definition: vector.h:193
Vector operator*(const double x) const
Definition: vector.h:172
double mean() const
Definition: vector.h:87
void load(const char *filename)
Vector subvect(const Index istart, const Index isize) const
Definition: vector.h:108
Vector(const Dimension N)
Definition: vector.h:31
Vector operator-() const
Definition: vector.h:61
void save(const char *filename) const
double & operator()(const Index i)
Definition: vector.h:51
void load(const std::string &s)
Definition: vector.h:95
void info() const
double sum() const
double * data() const
Definition: vector.h:44
Vector(const Matrix &A)
bool empty() const
Definition: vector.h:42
void operator/=(const double x)
Definition: vector.h:72
Vector(const SymMatrix &A)
Vect3 operator*(const double d, const Vect3 &V)
Definition: vect3.h:105
DeepCopy
Definition: linop.h:84
@ DEEP_COPY
Definition: linop.h:84
unsigned Dimension
Definition: linop.h:32
std::ostream & operator<<(std::ostream &os, const Conductivity< REP > &m)
std::istream & operator>>(std::istream &is, Conductivity< REP > &m)
unsigned Index
Definition: linop.h:33
BLAS_INT sizet_to_int(const unsigned &num)
Definition: linop.h:26