1 /*=========================================================================
4 Module: $RCSfile: vector.cxx,v $
6 Date: $Date: 2008/10/31 16:32:56 $
7 Version: $Revision: 1.1 $
9 Copyright: (c) 2002, 2003
12 This software is distributed WITHOUT ANY WARRANTY; without even
13 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
14 PURPOSE. See the above copyright notice for more information.
16 =========================================================================*/
20 #include <gsl/gsl_blas.h>
22 // ---------------------------------------------------------------------
23 std::ostream& operator<<( std::ostream& os, const gslobj_vector& v )
25 for( int i = 0; i < v.size( ); i++ ) os << " " << v( i );
29 // ---------------------------------------------------------------------
30 gslobj_vector::gslobj_vector( size_t s )
33 _gsl = gsl_vector_alloc( s );
36 // ---------------------------------------------------------------------
37 gslobj_vector::gslobj_vector( const gslobj_vector& v )
40 _gsl = gsl_vector_alloc( v._gsl->size );
41 gsl_vector_memcpy( _gsl, v._gsl );
44 // ---------------------------------------------------------------------
45 gslobj_vector::gslobj_vector( double* v, size_t s )
48 _view = gsl_vector_view_array( v, s );
49 _gsl = &( _view.vector );
52 // ---------------------------------------------------------------------
53 gslobj_vector::gslobj_vector( gsl_vector* v )
59 // ---------------------------------------------------------------------
60 /*==PS========================================
61 void gslobj_vector::resize( size_t s )
66 gsl_vector_free( _gsl );
67 _gsl = gsl_vector_alloc( s );
71 ==PS========================================*/
72 // ---------------------------------------------------------------------
73 gslobj_vector& gslobj_vector::operator=( const gslobj_vector& o )
75 gsl_vector_memcpy( _gsl, o._gsl );
79 // ---------------------------------------------------------------------
80 gslobj_vector& gslobj_vector::operator=( double o )
82 gsl_vector_set_all( _gsl, o );
86 // ---------------------------------------------------------------------
87 gslobj_vector& gslobj_vector::operator=( double* o )
89 memcpy( _gsl->data, o, _gsl->size * sizeof( double ) );
93 // ---------------------------------------------------------------------
94 gslobj_vector& gslobj_vector::operator=( gsl_vector* o )
96 gsl_vector_memcpy( _gsl, o );
100 // ---------------------------------------------------------------------
101 bool gslobj_vector::operator==( const gslobj_vector& o ) const
105 if( _gsl->size != o._gsl->size )
108 for( int i = 0; i < _gsl->size && ret; i++ )
109 ret = ret && ( _gsl->data[ i ] != o._gsl->data[ i ] );
113 // ---------------------------------------------------------------------
114 gslobj_vector gslobj_vector::operator+( const gslobj_vector& o )
116 gslobj_vector r( *this );
117 gsl_vector_add( r._gsl, o._gsl );
121 // ---------------------------------------------------------------------
122 gslobj_vector gslobj_vector::operator+( double o )
124 gslobj_vector r( *this );
125 gsl_vector_add_constant( r._gsl, o );
129 // ---------------------------------------------------------------------
130 gslobj_vector gslobj_vector::operator+( double* o )
132 gslobj_vector r( *this );
133 gsl_vector_view vw = gsl_vector_view_array( o, _gsl->size );
134 gsl_vector_add( r._gsl, &( vw.vector ) );
138 // ---------------------------------------------------------------------
139 gslobj_vector gslobj_vector::operator+( gsl_vector* o )
141 gslobj_vector r( *this );
142 gsl_vector_add( r._gsl, o );
146 // ---------------------------------------------------------------------
147 gslobj_vector& gslobj_vector::operator+=( const gslobj_vector& o )
149 gsl_vector_add( _gsl, o._gsl );
153 // ---------------------------------------------------------------------
154 gslobj_vector& gslobj_vector::operator+=( double o )
156 gsl_vector_add_constant( _gsl, o );
160 // ---------------------------------------------------------------------
161 gslobj_vector& gslobj_vector::operator+=( double* o )
163 gsl_vector_view vw = gsl_vector_view_array( o, _gsl->size );
164 gsl_vector_add( _gsl, &( vw.vector ) );
168 // ---------------------------------------------------------------------
169 gslobj_vector& gslobj_vector::operator+=( gsl_vector* o )
171 gsl_vector_add( _gsl, o );
175 // ---------------------------------------------------------------------
176 gslobj_vector gslobj_vector::operator-( const gslobj_vector& o )
178 gslobj_vector r( *this );
179 gsl_vector_sub( r._gsl, o._gsl );
183 // ---------------------------------------------------------------------
184 gslobj_vector gslobj_vector::operator-( double o )
186 gslobj_vector r( *this );
187 gsl_vector_add_constant( r._gsl, -o );
191 // ---------------------------------------------------------------------
192 gslobj_vector gslobj_vector::operator-( double* o )
194 gslobj_vector r( *this );
195 gsl_vector_view vw = gsl_vector_view_array( o, _gsl->size );
196 gsl_vector_sub( r._gsl, &( vw.vector ) );
200 // ---------------------------------------------------------------------
201 gslobj_vector gslobj_vector::operator-( gsl_vector* o )
203 gslobj_vector r( *this );
204 gsl_vector_sub( r._gsl, o );
208 // ---------------------------------------------------------------------
209 gslobj_vector& gslobj_vector::operator-=( const gslobj_vector& o )
211 gsl_vector_sub( _gsl, o._gsl );
215 // ---------------------------------------------------------------------
216 gslobj_vector& gslobj_vector::operator-=( double o )
218 gsl_vector_add_constant( _gsl, -o );
222 // ---------------------------------------------------------------------
223 gslobj_vector& gslobj_vector::operator-=( double* o )
225 gsl_vector_view vw = gsl_vector_view_array( o, _gsl->size );
226 gsl_vector_sub( _gsl, &( vw.vector ) );
230 // ---------------------------------------------------------------------
231 gslobj_vector& gslobj_vector::operator-=( gsl_vector* o )
233 gsl_vector_sub( _gsl, o );
237 // ---------------------------------------------------------------------
238 gslobj_vector gslobj_vector::operator*( double o )
240 gslobj_vector r( *this );
241 gsl_vector_scale( r._gsl, o );
245 // ---------------------------------------------------------------------
246 gslobj_vector& gslobj_vector::operator*=( double o )
248 gsl_vector_scale( _gsl, o );
252 // ---------------------------------------------------------------------
253 double gslobj_vector::dot( const gslobj_vector& o )
256 gsl_blas_ddot( _gsl, o._gsl, &r );
260 // ---------------------------------------------------------------------
261 double gslobj_vector::dot( double* o )
264 gsl_vector_view vw = gsl_vector_view_array( o, _gsl->size );
265 gsl_blas_ddot( _gsl, &( vw.vector ), &r );
269 // ---------------------------------------------------------------------
270 double gslobj_vector::dot( gsl_vector* o )
273 gsl_blas_ddot( _gsl, o, &r );
277 // ---------------------------------------------------------------------
278 gslobj_vector gslobj_vector::cross( const gslobj_vector& o )
280 return( cross( o._gsl->data ) );
283 // ---------------------------------------------------------------------
284 gslobj_vector gslobj_vector::cross( double* o )
286 gslobj_vector r( *this );
290 S = _gsl->data[ 1 ] * o[ 2 ];
291 s = _gsl->data[ 2 ] * o[ 1 ];
292 r._gsl->data[ 0 ] = S - s;
295 S = _gsl->data[ 2 ] * o[ 0 ];
296 s = _gsl->data[ 0 ] * o[ 2 ];
297 r._gsl->data[ 1 ] = S - s;
300 S = _gsl->data[ 0 ] * o[ 1 ];
301 s = _gsl->data[ 1 ] * o[ 0 ];
302 r._gsl->data[ 2 ] = S - s;
306 // ---------------------------------------------------------------------
307 gslobj_vector gslobj_vector::cross( gsl_vector* o )
309 return( cross( o->data ) );
312 // ---------------------------------------------------------------------
313 int gslobj_vector::scross( const gslobj_vector& o )
319 // ---------------------------------------------------------------------
320 int gslobj_vector::scross( double* o )
326 // ---------------------------------------------------------------------
327 int gslobj_vector::scross( gsl_vector* o )
333 // ---------------------------------------------------------------------
334 gslobj_vector gslobj_vector::operator/( double o )
336 gslobj_vector r( *this );
337 gsl_vector_scale( r._gsl, 1.0 / o );
341 // ---------------------------------------------------------------------
342 gslobj_vector& gslobj_vector::operator/=( double o )
344 gsl_vector_scale( _gsl, 1.0 / o );
348 // ---------------------------------------------------------------------
349 /*==PS========================================
350 double gslobj_vector::norm1( )
352 return( gsl_blas_dasum( _gsl ) );
354 ==PS========================================*/
355 // ---------------------------------------------------------------------
356 double gslobj_vector::norm2( )
358 return( gsl_blas_dnrm2( _gsl ) );
361 // ---------------------------------------------------------------------
362 gslobj_vector gslobj_vector::normalize( )
364 gslobj_vector r( *this );
365 gsl_vector_scale( r._gsl, 1.0 / gsl_blas_dnrm2( _gsl ) );
369 // ---------------------------------------------------------------------
370 int gslobj_vector::snormalize( )
372 return( gsl_vector_scale( _gsl, 1.0 / gsl_blas_dnrm2( _gsl ) ) );
375 // eof - gslobj_vector.cxx