00001
00002
00003
00004
00005
00006
00007
00008 #ifndef synaps_linalg_vector_H
00009 #define synaps_linalg_vector_H
00010
00011 #include <cassert>
00012 #include <synaps/init.h>
00013 #include <synaps/arithm/let.h>
00014
00015 __BEGIN_NAMESPACE_SYNAPS
00016
00018
00039 namespace VECTOR {
00040
00042 template <class OS, class R>
00043 inline OS & print(OS & os, const R & v) {
00044 os <<"[";
00045 typename R::const_iterator i=v.begin();
00046 if (v.size()){
00047 os<<*i; ++i;
00048 for(; i != v.end(); ++i) os <<","<<*i;
00049 }
00050 os <<"]";
00051 return os;
00052 }
00053
00055 template <class OS, class R>
00056 inline OS & print(OS & os, const R & v, unsigned sz) {
00057 os <<"[";
00058 unsigned i=0;
00059 if (sz){
00060 os<<v; ++i;
00061 for(; i <sz; ++i) os <<","<<v[i];
00062 }
00063 os <<"]";
00064 return os;
00065 }
00066
00071 template <class IS, class R>
00072 inline IS & read(IS & is, R & V) {
00073 typedef typename R::size_type size_type;
00074 size_type s;
00075 is >> s;
00076 V.rep().resize(s);
00077 for(size_type i=0; i< s; ++i) is >> V;
00078 return(is);
00079 }
00080
00081
00082 template <class V1, class V2> inline
00083 bool equal_n( const V1& v1, const V2 & v2, int n )
00084 {
00085 int i=0;
00086 while(i < n && v1[i] == v2[i]) i++;
00087 return i == n;
00088 };
00089
00091 template <class V1, class V2> inline
00092 bool equal(const V1 & v1, const V2 & v2)
00093 { return (v1.size() == v2.size()) && equal_n(v1,v2,v1.size()); }
00094
00095
00096
00097 template <class V, class C> inline
00098 void init_n(V & a, const C & v, int n ) {
00099 for ( int i = 0; i < n; a = v, i ++ ); };
00101 template <class V, class C>
00102 void init(V & a, const C & v ) { init_n(a,v,(int)a.size()); };
00103
00104
00106 template <class V>
00107 void reverse_n( V & v, int a, int n )
00108 {
00109 int k = n - a;
00110 for ( int i = 0; i < k/2; std::swap(v[a+i],v[k-1-i]), i ++ );
00111 };
00112
00113 template <class V>
00114 void reverse( V & v, int n ) { reverse_n( v, 0, n ); };
00115
00117 template <class V,class I>
00118 void reverse(V & v, I d) { for(I i=0;i<d/2; std::swap(v[i],v[d-1-i]),++i); }
00119
00120
00121 template<class V, class W> inline
00122 void neg_range( V & r, const W & w, int a, int b )
00123 { for ( int i = a; i < b; r = -w, i++ ); };
00124 template<class V, class W> inline
00125 void neg_range( V & r, int a, int b )
00126 { for ( int i = a; i < b; r = -r, i++ ); };
00127
00128 template<class VI> inline
00129 void neg_i( VI a, VI eda ) { for ( ; a != eda; *a = -(*a), a++ ); };
00130 template<class VI> inline
00131 void neg_i( VI a, VI eda, VI b ) { for ( ; a != eda; *a = -(*b), a++, b++ ); };
00132 template<class V> inline
00133 void neg( V & v ) { neg_i(v.begin(),v.end()); };
00134 template<class V> inline
00135 void neg( V & a, const V & b )
00136 {
00137 a.resize( b.size() );
00138 neg_i( a.begin(), a.end(), b.begin() );
00139 };
00140
00141 template<class V, class W> inline
00142 void copy_range( V & v, const W & w, int a, int b ) { for ( int i = a; i < b; v = w, i ++ ); }
00143 template< class V, class W >
00144 void add_n(V & a, const W & b, int n ) {for ( int i = 0; i < n; a += b, i ++ ); }
00145 template<class VI, class WI>
00146 void add_i( VI a, VI ea, WI b ) { for ( ; a != ea; *a += *b, ++a, ++b ); };
00147
00149 template <class V, class W>
00150 void add(V & a, const W & b)
00151 {
00152 int asz, bsz;
00153 asz = a.size(), bsz = b.size();
00154 if ( asz < bsz )
00155 {
00156 a.resize(bsz);
00157 add_i(a.begin(),a.end(),b.begin());
00158 copy_range(a,b,asz,bsz);
00159 }
00160 else
00161 add_i(a.begin(),a.begin()+bsz,b.begin());
00162 }
00163
00164 template< class V, class W, class X > inline
00165 void add_n(V & a, const W & b, const X& c, int n ) { for ( int i = 0; i < n; a = b+c, i ++ ); };
00166
00167
00170 template <class V, class W,class X>
00171 void add_g(V & r, const W & a, const X & b)
00172 {
00173 int asz, bsz, rsz;
00174 asz = a.size(), bsz = b.size(), rsz = r.size();
00175
00176 if ( asz > bsz )
00177 {
00178 if ( asz > rsz ) r.resize(asz);
00179 add_n(r,a,b,bsz);
00180 copy_range(r,a,bsz,asz);
00181 }
00182 else
00183 {
00184 if ( bsz > rsz ) r.resize(bsz);
00185 add_n(r,a,b,asz);
00186 copy_range(r,b,asz,bsz);
00187 };
00188 }
00189
00192 template <class V, class W,class X> inline
00193 void add(V & r, const W & a, const X & b)
00194 {
00195 add_g(r,a,b);
00196 }
00197
00200 template <class V> inline
00201 void add(V & r, const V & a, const V & b)
00202 {
00203 if(&r==&a) { add(r,b); return; };
00204 if(&r==&b) { add(r,a); return; };
00205 add_g(r,a,b);
00206 }
00207
00208 template< class V, class W >
00209 void sub_n(V & a, const W & b, int n ) { for ( int i = 0; i < n; a -= b, i ++ ); };
00210
00211
00213 template <class V, class W>
00214 void sub(V & a, const W & b)
00215 {
00216 int asz, bsz;
00217 asz = a.size(), bsz = b.size();
00218 if ( asz < bsz )
00219 {
00220 a.resize(bsz);
00221 sub_n(a,b,asz);
00222 neg_range(a,b,asz,bsz);
00223 }
00224 else
00225 sub_n(a,b,bsz);
00226 }
00227
00228 template< class V, class W, class X > inline
00229 void sub_n(V & a, const W & b, const X& c, int n ) { for ( int i = 0; i < n; a = b-c, i ++ ); };
00230
00231
00232 template <class V, class W,class X>
00233 void sub_g(V & r, const W & a, const X & b)
00234 {
00235 int asz, bsz, rsz;
00236 asz = a.size(), bsz = b.size(), rsz = r.size();
00237
00238 if ( asz > bsz )
00239 {
00240 if ( asz > rsz ) r.resize(asz);
00241 sub_n(r,a,b,bsz);
00242 copy_range(r,a,bsz,asz);
00243 }
00244 else
00245 {
00246 if ( bsz > rsz ) r.resize(bsz);
00247 sub_n(r,a,b,asz);
00248 neg_range(r,b,asz,bsz);
00249 };
00250 }
00253 template <class V, class W,class X>
00254 void sub(V & r, const W & a, const X & b)
00255 {
00256 sub_g(r,a,b);
00257 }
00258
00261 template <class V>
00262 void sub(V & r, const V & a, const V & b)
00263 {
00264 if(&r==&a) { sub(r,b); return; };
00265
00266 sub_g(r,a,b);
00267 }
00268
00269 template<class V, class W> inline
00270 void mul_ext_n( V & a, const W & c, int n )
00271 { for ( int i = 0; i < n; a = a*c, i ++ ); };
00272 template<class V, class W, class S> inline
00273 void mul_ext_n( V & a, const W & b, const S& c, int n )
00274 { for ( int i = 0; i < n; a = b*c, i++ ); };
00275 template<class VI, class W> inline
00276 void mul_ext_i( VI bga, VI eda, const W & c )
00277 { for ( VI i = bga; i != eda; *i *= c, ++i ); };
00278 template<class VI, class VII, class W> inline
00279 void mul_ext_i( VI bga, VI eda, VII bgb, const W & c )
00280 { for ( VI a = bga; a != eda; *a = *bgb * c, ++bgb, ++a ); };
00281
00283 template <class V,class W>
00284 void mul_ext(V & a, const V & b, const W & c)
00285 {
00286 if(&a != &b)
00287 {
00288 a.resize( b.size() );
00289 mul_ext_i(a.begin(),a.end(),b.begin(),c);
00290 }
00291 else mul_ext_i(a.begin(),a.end(),c);
00292 }
00293
00294
00295 template<class V, class SC> inline
00296 void div_ext_n( V & a, const SC & sc, int n )
00297 { for ( int i = 0; i < n; a /= sc, i++); }
00298 template<class V1, class V2, class SC>
00299 void div_ext_n( V1 & a, const V2 & b, const SC& sc, int n )
00300 { for ( int i = 0; i < n; a = b/sc, i ++ ); };
00301 template<class VI, class SC> inline
00302 void div_ext_i( VI bga, VI eda, const SC & sc )
00303 { for ( ; bga < eda; *bga /= sc, bga++); }
00304 template<class VI, class SC> inline
00305 void div_ext_i( VI a, VI eda, VI b, const SC & sc )
00306 { for ( ; a < eda; *a = *b / sc, a++, b ++ ); };
00308 template <class V, class SC>
00309 void div_ext( V & a, const V & b, const SC & sc)
00310 {
00311 if(&a != &b) {
00312 a.resize( b.size() );
00313 div_ext_i(a.begin(),a.end(),b.begin(),sc);
00314 }
00315 else div_ext_i(a.begin(),a.end(),sc);
00316 }
00318 template <class V, class W>
00319 void div_ext(V & a, const W & c)
00320 { div_ext_i(a.begin(),a.end(),c); }
00321
00323 template <class V, class W>
00324 void div(V & a, const W & c)
00325 {
00326 div_ext(a,c);
00327
00328 }
00329
00330
00331 template <class R, class S, class C> inline
00332 void apply_mult(R & result, const S& a, const C & m)
00333 {
00334 typedef typename S::const_iterator const_iterator;
00335 typedef typename R::iterator iterator;
00336 result.resize(a.size());
00337 const_iterator b=a.begin(); iterator i=result.begin();
00338 for(; b!=a.end(); ++i, ++b) *i = (*b) * m;
00339 }
00340
00342 template <class C, class R>
00343 C norm(const R & v)
00344 {
00345 C r=0;
00346 for(typename R::const_iterator it=v.begin(); it !=v.end(); ++it) r+=(*it)*(*it);
00347 return sqrt(r);
00348 }
00349
00351 template <class C,class R>
00352 C norm(const R & v, int p)
00353 {
00354 C r=0;
00355 for(typename R::const_iterator it=v.begin(); it !=v.end(); ++it) r+=pow(*it,p);
00356 return pow(r,1/p);
00357 }
00358
00360 template <class R, class S>
00361 typename R::value_type innerprod_n( const S & v, const R & w, unsigned n )
00362 {
00363 std::cout << "innerprod \n";
00364 typename R::value_type r = 0;
00365 for ( unsigned i = 0; i < n; r += v*w, i++ );
00366 return r;
00367 };
00368
00369 template <class R, class S>
00370 typename R::value_type innerprod(const S & v, const R & w)
00371 {
00372 typename R::value_type r =0;
00373 for(unsigned int i=0; i <v.size(); i++){r += (v[i]*w[i]);}
00374 return r;
00375 }
00376
00377 template<class R>
00378 typename R::value_type max_abs(const R& v)
00379 {
00380 typename R::value_type r=0;
00381 for(unsigned int i=0; i <v.size(); i++)
00382 {
00383 if(r>v[i]) r=v;
00384 else if(r>-v[i]) r=-v;
00385 }
00386 return r;
00387 }
00388
00389
00390 template<class Ia, class Ib> inline
00391 void assign_i( Ia a, Ia eda, Ib b )
00392 {
00393 for ( ; a != eda; ++a, ++b ) let::assign(*a,*b);
00394 };
00395
00396 template<class Ra, class Rb> inline
00397 void assign( Ra & a, const Rb & b )
00398 {
00399 a.resize( b.size() );
00400 assign_i( a.begin(), a.end(), b.begin() );
00401 };
00402
00403 template<class R, class C> inline
00404 void vaddsc( R & r, const C & x )
00405 { for ( typename R::iterator i = r.begin(); i != r.end(); *i += x, i ++ ); };
00406
00407 }
00408
00409 __END_NAMESPACE_SYNAPS
00410
00411 #endif // SYNAPS_LINALG_VECTOR_H
00412
