OpenMEEG
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
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 }
double * data() const
Definition: vector.h:77
size_t size() const
Definition: vector.h:73
utils::RCPtr< LinOpValue > value
Definition: vector.h:58
std::istream & operator>>(std::istream &is, Conductivity< REP > &m)
Vector operator/(double x) const
Definition: vector.h:99
Vector operator*(double x) const
Definition: vector.h:194
double & operator()(const size_t i)
Definition: vector.h:84
std::ostream & operator<<(std::ostream &os, const Conductivity< REP > &m)
Vect3 operator*(const double &d, const Vect3 &v)
Definition: vect3.h:151
void alloc_data()
Definition: vector.h:70
Vector(const Vector &A, const DeepCopy)
Definition: vector.h:65
void save(const std::string &s) const
Definition: vector.h:116
double norm() const
Definition: vector.h:215
void operator+=(const Vector &v)
Definition: vector.h:162
void operator*=(double x)
Definition: vector.h:206
double mean() const
Definition: vector.h:110
OPENMEEGMATHS_EXPORT BLAS_INT sizet_to_int(const size_t &num)
Definition: linop.h:56
Vector operator-(const Vector &v) const
Definition: vector.h:150
void reference_data(const double *array)
Definition: vector.h:71
bool empty() const
Definition: vector.h:75
Vector operator+(const Vector &v) const
Definition: vector.h:138
size_t nlin() const
Definition: linop.h:85
Vector subvect(size_t istart, size_t isize) const
Definition: vector.h:130
void operator-=(const Vector &v)
Definition: vector.h:172
Vector(const size_t N)
Definition: vector.h:64
double operator()(const size_t i) const
Definition: vector.h:79
DeepCopy
Definition: linop.h:126
Matrix class.
Definition: matrix.h:61
void operator/=(double x)
Definition: vector.h:95
void load(const std::string &s)
Definition: vector.h:117