OpenMEEG
vector.h
Go to the documentation of this file.
1 /*
2 Project Name : OpenMEEG
3 
4 © INRIA and ENPC (contributors: Geoffray ADDE, Maureen CLERC, Alexandre
5 GRAMFORT, Renaud KERIVEN, Jan KYBIC, Perrine LANDREAU, Théodore PAPADOPOULO,
6 Emmanuel OLIVI
7 Maureen.Clerc.AT.inria.fr, keriven.AT.certis.enpc.fr,
8 kybic.AT.fel.cvut.cz, papadop.AT.inria.fr)
9 
10 The OpenMEEG software is a C++ package for solving the forward/inverse
11 problems of electroencephalography and magnetoencephalography.
12 
13 This software is governed by the CeCILL-B license under French law and
14 abiding by the rules of distribution of free software. You can use,
15 modify and/ or redistribute the software under the terms of the CeCILL-B
16 license as circulated by CEA, CNRS and INRIA at the following URL
17 "http://www.cecill.info".
18 
19 As a counterpart to the access to the source code and rights to copy,
20 modify and redistribute granted by the license, users are provided only
21 with a limited warranty and the software's authors, the holders of the
22 economic rights, and the successive licensors have only limited
23 liability.
24 
25 In this respect, the user's attention is drawn to the risks associated
26 with loading, using, modifying and/or developing or reproducing the
27 software by the user in light of its specific status of free software,
28 that may mean that it is complicated to manipulate, and that also
29 therefore means that it is reserved for developers and experienced
30 professionals having in-depth computer knowledge. Users are therefore
31 encouraged to load and test the software's suitability as regards their
32 requirements in conditions enabling the security of their systems and/or
33 data to be ensured and, more generally, to use and operate it in the
34 same conditions as regards security.
35 
36 The fact that you are presently reading this means that you have had
37 knowledge of the CeCILL-B license and that you accept its terms.
38 */
39 
40 #pragma once
41 
42 #include <OMassert.H>
43 #include <cstdlib>
44 #include <string>
45 
46 #include <OpenMEEGMathsConfig.h>
47 #include <linop.h>
48 #include <RC.H>
49 #include <MathsIO.H>
50 
51 namespace OpenMEEG {
52 
53  class Matrix;
54  class SymMatrix;
55 
56  class OPENMEEGMATHS_EXPORT Vector: public LinOp {
57 
58  utils::RCPtr<LinOpValue> value;
59 
60  public:
61 
62  Vector(): LinOp(0,1,FULL,1),value() { }
63 
64  Vector(const size_t N): LinOp(N,1,FULL,1),value(new LinOpValue(size())) { }
65  Vector(const Vector& A,const DeepCopy): LinOp(A.nlin(),1,FULL,1),value(new LinOpValue(A.size(),A.data())) { }
66 
67  explicit Vector(Matrix& A);
68  explicit Vector(SymMatrix& A);
69 
70  void alloc_data() { value = new LinOpValue(size()); }
71  void reference_data(const double* array) { value = new LinOpValue(size(),array); }
72 
73  size_t size() const { return nlin(); }
74 
75  bool empty() const { return value->empty(); }
76 
77  double* data() const { return value->data; }
78 
79  inline double operator()(const size_t i) const {
80  om_assert(i<nlin());
81  return value->data[i];
82  }
83 
84  inline double& operator()(const size_t i) {
85  om_assert(i<nlin());
86  return value->data[i];
87  }
88 
89  Vector subvect(size_t istart, size_t isize) const;
90  Vector operator+(const Vector& v) const;
91  Vector operator-(const Vector& v) const;
92  void operator+=(const Vector& v);
93  void operator-=(const Vector& v);
94  void operator*=(double x);
95  void operator/=(double x) { (*this) *= (1.0/x); }
96  Vector operator+(double i) const;
97  Vector operator-(double i) const;
98  Vector operator*(double x) const;
99  Vector operator/(double x) const { return (*this)*(1.0/x); }
100  double operator*(const Vector& v) const;
101  Vector operator*(const Matrix& m) const;
102 
103  Vector kmult(const Vector& x) const;
104  // Vector conv(const Vector& v) const;
105  // Vector conv_trunc(const Vector& v) const;
106  Matrix outer_product(const Vector& v) const;
107 
108  double norm() const;
109  double sum() const;
110  double mean() const { return sum()/size(); }
111 
112  void set(double x);
113  void save(const char *filename) const;
114  void load(const char *filename);
115 
116  void save(const std::string& s) const { save(s.c_str()); }
117  void load(const std::string& s) { load(s.c_str()); }
118 
119  void info() const;
120 
121  friend class SymMatrix;
122  friend class Matrix;
123  };
124 
125  OPENMEEGMATHS_EXPORT Vector operator*(const double &d, const Vector &v);
126 
127  OPENMEEGMATHS_EXPORT std::ostream& operator<<(std::ostream& f,const Vector &M);
128  OPENMEEGMATHS_EXPORT std::istream& operator>>(std::istream& f,Vector &M);
129 
130  inline Vector Vector::subvect(size_t istart, size_t isize) const {
131  om_assert (istart+isize<=nlin());
132  Vector a(isize);
133  for (size_t i=0; i<isize; i++)
134  a(i) = (*this)(istart+i);
135  return a;
136  }
137 
138  inline Vector Vector::operator+(const Vector& v) const {
139  om_assert(nlin()==v.nlin());
140  Vector p(*this,DEEP_COPY);
141  #ifdef HAVE_BLAS
142  BLAS(daxpy,DAXPY)(sizet_to_int(nlin()),1,v.data(),1,p.data(),1);
143  #else
144  for( size_t i=0; i<nlin(); i++ )
145  p.data()[i]=data()[i]+v.data()[i];
146  #endif
147  return p;
148  }
149 
150  inline Vector Vector::operator-(const Vector& v) const {
151  om_assert(nlin()==v.nlin());
152  Vector p(*this,DEEP_COPY);
153  #ifdef HAVE_BLAS
154  BLAS(daxpy,DAXPY)(sizet_to_int(nlin()),-1,v.data(),1,p.data(),1);
155  #else
156  for( size_t i=0; i<nlin(); i++ )
157  p.data()[i]=data()[i]-v.data()[i];
158  #endif
159  return p;
160  }
161 
162  inline void Vector::operator+=(const Vector& v) {
163  om_assert(nlin()==v.nlin());
164  #ifdef HAVE_BLAS
165  BLAS(daxpy,DAXPY)(sizet_to_int(nlin()),1,v.data(),1,data(),1);
166  #else
167  for( size_t i=0; i<nlin(); i++ )
168  data()[i]+=v.data()[i];
169  #endif
170  }
171 
172  inline void Vector::operator-=(const Vector& v) {
173  om_assert(nlin()==v.nlin());
174  #ifdef HAVE_BLAS
175  BLAS(daxpy,DAXPY)(sizet_to_int(nlin()),-1,v.data(),1,data(),1);
176  #else
177  for( size_t i=0; i<nlin(); i++ )
178  data()[i]-=v.data()[i];
179  #endif
180  }
181 
182  inline double Vector::operator*(const Vector& v) const {
183  om_assert(nlin()==v.nlin());
184  #ifdef HAVE_BLAS
185  return BLAS(ddot,DDOT)(sizet_to_int(nlin()),data(),1,v.data(),1);
186  #else
187  double s=0;
188  for( size_t i=0; i<nlin(); i++ )
189  s+=data()[i]*v.data()[i];
190  return s;
191  #endif
192  }
193 
194  inline Vector Vector::operator*(double x) const {
195  #ifdef HAVE_BLAS
196  Vector p(*this,DEEP_COPY);
197  BLAS(dscal,DSCAL)(sizet_to_int(nlin()),x,p.data(),1);
198  #else
199  Vector p(nlin());
200  for( size_t i=0; i<nlin(); i++ )
201  p.data()[i]=x*data()[i];
202  #endif
203  return p;
204  }
205 
206  inline void Vector::operator*=(double x) {
207  #ifdef HAVE_BLAS
208  BLAS(dscal,DSCAL)(sizet_to_int(nlin()),x,data(),1);
209  #else
210  for( size_t i=0; i<nlin(); i++ )
211  data()[i]*=x;
212  #endif
213  }
214 
215  inline double Vector::norm() const
216  {
217  #ifdef HAVE_BLAS
218  return BLAS(dnrm2,DNRM2)(sizet_to_int(nlin()),data(),1);
219  #else
220  std::cout << "'Vector::norm' not implemented" << std::endl;
221  exit(1);
222  return 0;
223  #endif
224  }
225 
226  // inline Vector Vector::conv(const Vector& v) const {
227  // if (v.nlin()<nlin()) return v.conv(*this);
228  //
229  // Vector p(nlin()+v.nlin()-1);
230  // p.set(0);
231  // for (size_t i=0; i<v.nlin(); i++) {
232  // #ifdef HAVE_BLAS
233  // BLAS(daxpy,DAXPY)(sizet_to_int(nlin()), v(i), data(), 1, p.data()+i, 1);
234  // #else
235  // for (size_t j=0;j<nlin();j++)
236  // p(i+j)+=v(i)*data()[j];
237  // #endif
238  // }
239  // return p;
240  // }
241  //
242  // inline Vector Vector::conv_trunc(const Vector& v) const {
243  // Vector p(v.nlin());
244  // p.set(0);
245  // for (size_t i=0; i<v.nlin(); i++)
246  // {
247  // size_t m = std::min(nlin(),v.nlin()-i);
248  // #ifdef HAVE_BLAS
249  // BLAS(daxpy,DAXPY)((int)m, v(i), data(), 1, p.data()+i, 1);
250  // #else
251  // for (size_t j=0;j<m;j++)
252  // p(i+j)+=v(i)*data()[j];
253  // #endif
254  // }
255  // return p;
256  // }
257 
258  // Operators.
259 
260  OPENMEEGMATHS_EXPORT inline Vector operator*(const double &d, const Vector &v) { return v*d; }
261 }
OpenMEEG::Vector::operator()
double operator()(const size_t i) const
Definition: vector.h:79
OpenMEEG::Vector::save
void save(const std::string &s) const
Definition: vector.h:116
OpenMEEG::DEEP_COPY
@ DEEP_COPY
Definition: linop.h:113
OpenMEEG::Vector::operator-
Vector operator-(const Vector &v) const
Definition: vector.h:150
OpenMEEG::Vector::subvect
Vector subvect(size_t istart, size_t isize) const
Definition: vector.h:130
linop.h
OpenMEEG::Vector::operator+=
void operator+=(const Vector &v)
Definition: vector.h:162
OpenMEEG::Vector::operator()
double & operator()(const size_t i)
Definition: vector.h:84
OpenMEEG::Vector::reference_data
void reference_data(const double *array)
Definition: vector.h:71
OpenMEEG::LinOp
Definition: linop.h:100
OpenMEEG::LinOpInfo::nlin
size_t nlin() const
Definition: linop.h:77
OpenMEEG
Definition: analytics.h:45
OpenMEEG::operator>>
std::istream & operator>>(std::istream &is, Conductivity< REP > &m)
Definition: PropertiesSpecialized.h:60
OpenMEEG::Vector::norm
double norm() const
Definition: vector.h:215
OpenMEEG::Vector::load
void load(const std::string &s)
Definition: vector.h:117
OpenMEEG::Vector::Vector
Vector(const Vector &A, const DeepCopy)
Definition: vector.h:65
OpenMEEG::Vector::operator*=
void operator*=(double x)
Definition: vector.h:206
OpenMEEG::Vector::value
utils::RCPtr< LinOpValue > value
Definition: vector.h:58
OpenMEEG::Vector::Vector
Vector()
Definition: vector.h:62
OpenMEEG::DeepCopy
DeepCopy
Definition: linop.h:113
OpenMEEG::Vector::size
size_t size() const
Definition: vector.h:73
OpenMEEG::LinOpValue
Definition: linop.h:115
OpenMEEG::Vector
Definition: vector.h:56
OpenMEEG::Matrix
Matrix class.
Definition: matrix.h:61
OpenMEEG::Vector::operator*
Vector operator*(double x) const
Definition: vector.h:194
OpenMEEGMathsConfig.h
OpenMEEG::Vector::Vector
Vector(const size_t N)
Definition: vector.h:64
OpenMEEG::Vector::empty
bool empty() const
Definition: vector.h:75
OpenMEEG::operator<<
std::ostream & operator<<(std::ostream &os, const Conductivity< REP > &m)
Definition: PropertiesSpecialized.h:63
OpenMEEG::sizet_to_int
OPENMEEGMATHS_EXPORT BLAS_INT sizet_to_int(const size_t &num)
Definition: linop.h:56
OpenMEEG::Vector::operator/=
void operator/=(double x)
Definition: vector.h:95
OpenMEEG::Vector::mean
double mean() const
Definition: vector.h:110
OpenMEEG::Vector::operator+
Vector operator+(const Vector &v) const
Definition: vector.h:138
OpenMEEG::operator*
Vect3 operator*(const double &d, const Vect3 &v)
Definition: vect3.h:151
OpenMEEG::Vector::data
double * data() const
Definition: vector.h:77
OpenMEEG::Vector::operator-=
void operator-=(const Vector &v)
Definition: vector.h:172
OpenMEEG::Vector::operator/
Vector operator/(double x) const
Definition: vector.h:99
OpenMEEG::SymMatrix
Definition: symmatrix.h:53
OpenMEEG::Vector::alloc_data
void alloc_data()
Definition: vector.h:70