OpenMEEG
Loading...
Searching...
No Matches
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
18namespace 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 & operator()(const Index i)
Definition: vector.h:51
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
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::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
std::ostream & operator<<(std::ostream &os, const Conductivity< REP > &m)