00001
00002 #ifndef vnl_c_vector_txx_
00003 #define vnl_c_vector_txx_
00004
00005
00006
00007
00008
00009
00010
00011 #include "vnl_c_vector.h"
00012 #include <vcl_cmath.h>
00013 #include <vcl_cassert.h>
00014 #include <vnl/vnl_math.h>
00015 #include <vnl/vnl_complex_traits.h>
00016 #include <vnl/vnl_numeric_traits.h>
00017
00018 #include <vnl/vnl_sse.h>
00019
00020 template <class T>
00021 T vnl_c_vector<T>::sum(T const* v, unsigned n)
00022 {
00023 return vnl_sse<T>::sum(v,n);
00024 }
00025
00026 template <class T>
00027 void vnl_c_vector<T>::normalize(T* v, unsigned n)
00028 {
00029 typedef typename vnl_numeric_traits<T>::abs_t abs_t;
00030 typedef typename vnl_numeric_traits<abs_t>::real_t real_t;
00031 abs_t tmp(0);
00032 for (unsigned i = 0; i < n; ++i)
00033 tmp += vnl_math_squared_magnitude(v[i]);
00034 if (tmp!=0)
00035 {
00036 tmp = abs_t(real_t(1) / vcl_sqrt(real_t(tmp)));
00037 for (unsigned i = 0; i < n; ++i)
00038 v[i] = T(tmp*v[i]);
00039 }
00040 }
00041
00042 template <class T>
00043 void vnl_c_vector<T>::apply(T const* v, unsigned n, T (*f)(T const&), T* v_out)
00044 {
00045 for (unsigned i = 0; i < n; ++i)
00046 v_out[i] = f(v[i]);
00047 }
00048
00049 template <class T>
00050 void vnl_c_vector<T>::apply(T const* v, unsigned n, T (*f)(T), T* v_out)
00051 {
00052 for (unsigned i = 0; i < n; ++i)
00053 v_out[i] = f(v[i]);
00054 }
00055
00056 template <class T>
00057 void vnl_c_vector<T>::copy(T const *src, T *dst, unsigned n)
00058 {
00059 for (unsigned i=0; i<n; ++i)
00060 dst[i] = src[i];
00061 }
00062
00063 template <class T>
00064 void vnl_c_vector<T>::scale(T const *x, T *y, unsigned n, T const &a_)
00065 {
00066 T a = a_;
00067 if (x == y)
00068 for (unsigned i=0; i<n; ++i)
00069 y[i] *= a;
00070 else
00071 for (unsigned i=0; i<n; ++i)
00072 y[i] = a*x[i];
00073 }
00074
00075
00076 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00077 #define impl_elmt_wise_commutative(op) \
00078 if (z == x) \
00079 for (unsigned i=0; i<n; ++i) \
00080 z[i] op##= y[i]; \
00081 else if (z == y) \
00082 for (unsigned i=0; i<n; ++i) \
00083 z[i] op##= x[i]; \
00084 else \
00085 for (unsigned i=0; i<n; ++i) \
00086 z[i] = x[i] op y[i];
00087
00088 #define impl_elmt_wise_non_commutative(op) \
00089 if (z == x) \
00090 for (unsigned i=0; i<n; ++i) \
00091 z[i] op##= y[i]; \
00092 else \
00093 for (unsigned i=0; i<n; ++i) \
00094 z[i] = x[i] op y[i];
00095
00096 #define impl_elmt_wise_commutative_a(op) \
00097 if (z == x) \
00098 for (unsigned i=0; i<n; ++i) \
00099 z[i] op##= y; \
00100 else \
00101 for (unsigned i=0; i<n; ++i) \
00102 z[i] = x[i] op y;
00103
00104 #define impl_elmt_wise_non_commutative_a(op) \
00105 if (z == x) \
00106 for (unsigned i=0; i<n; ++i) \
00107 z[i] op##= y; \
00108 else \
00109 for (unsigned i=0; i<n; ++i) \
00110 z[i] = x[i] op y;
00111 #endif // DOXYGEN_SHOULD_SKIP_THIS
00112
00113 template <class T>
00114 void vnl_c_vector<T>::add(T const *x, T const *y, T *z, unsigned n)
00115 {
00116 impl_elmt_wise_commutative(+);
00117 }
00118
00119 template <class T>
00120 void vnl_c_vector<T>::add(T const *x, T const& y, T *z, unsigned n)
00121 {
00122 impl_elmt_wise_commutative_a(+);
00123 }
00124
00125 template <class T>
00126 void vnl_c_vector<T>::subtract(T const *x, T const *y, T *z, unsigned n)
00127 {
00128 impl_elmt_wise_non_commutative(-);
00129 }
00130
00131 template <class T>
00132 void vnl_c_vector<T>::subtract(T const *x, T const& y, T *z, unsigned n)
00133 {
00134 impl_elmt_wise_commutative_a(-);
00135 }
00136
00137 template <class T>
00138 void vnl_c_vector<T>::multiply(T const *x, T const *y, T *z, unsigned n)
00139 {
00140 impl_elmt_wise_commutative(*);
00141 }
00142
00143 template <class T>
00144 void vnl_c_vector<T>::multiply(T const *x, T const& y, T *z, unsigned n)
00145 {
00146 impl_elmt_wise_commutative_a(*);
00147 }
00148
00149 template <class T>
00150 void vnl_c_vector<T>::divide(T const *x, T const *y, T *z, unsigned n)
00151 {
00152 impl_elmt_wise_non_commutative(/);
00153 }
00154
00155 template <class T>
00156 void vnl_c_vector<T>::divide(T const *x, T const& y, T *z, unsigned n)
00157 {
00158 impl_elmt_wise_commutative_a(/);
00159 }
00160
00161 #undef impl_elmt_wise_commutative
00162 #undef impl_elmt_wise_noncommutative
00163
00164
00165 template <class T>
00166 void vnl_c_vector<T>::negate(T const *x, T *y, unsigned n)
00167 {
00168 if (x == y)
00169 for (unsigned i=0; i<n; ++i)
00170 y[i] = -y[i];
00171 else
00172 for (unsigned i=0; i<n; ++i)
00173 y[i] = -x[i];
00174 }
00175
00176 template <class T>
00177 void vnl_c_vector<T>::invert(T const *x, T *y, unsigned n)
00178 {
00179 if (x == y)
00180 for (unsigned i=0; i<n; ++i)
00181 y[i] = T(1)/y[i];
00182 else
00183 for (unsigned i=0; i<n; ++i)
00184 y[i] = T(1)/x[i];
00185 }
00186
00187 template <class T>
00188 void vnl_c_vector<T>::saxpy(T const &a_, T const *x, T *y, unsigned n)
00189 {
00190 T a = a_;
00191 for (unsigned i=0; i<n; ++i)
00192 y[i] += a*x[i];
00193 }
00194
00195 template <class T>
00196 void vnl_c_vector<T>::fill(T *x, unsigned n, T const &v_)
00197 {
00198 T v = v_;
00199 for (unsigned i=0; i<n; ++i)
00200 x[i] = v;
00201 }
00202
00203 template <class T>
00204 void vnl_c_vector<T>::reverse(T *x, unsigned n)
00205 {
00206 for (unsigned i=0; 2*i+1<n; ++i) {
00207 T tmp = x[i];
00208 x[i] = x[n-1-i];
00209 x[n-1-i] = tmp;
00210 }
00211 }
00212
00213
00214 template<class T>
00215 T vnl_c_vector<T>::dot_product(T const *a, T const *b, unsigned n)
00216 {
00217 return vnl_sse<T>::dot_product(a,b,n);
00218 }
00219
00220
00221 template<class T>
00222 T vnl_c_vector<T>::inner_product(T const *a, T const *b, unsigned n)
00223 {
00224 T ip(0);
00225 for (unsigned i=0; i<n; ++i)
00226 ip += a[i] * vnl_complex_traits<T>::conjugate(b[i]);
00227 return ip;
00228 }
00229
00230
00231 template<class T>
00232 void vnl_c_vector<T>::conjugate(T const *src, T *dst, unsigned n)
00233 {
00234 for (unsigned i=0; i<n; ++i)
00235 dst[i] = vnl_complex_traits<T>::conjugate( src[i] );
00236 }
00237
00238
00239
00240
00241 template<class T>
00242 T vnl_c_vector<T>::max_value(T const *src, unsigned n)
00243 {
00244 assert(n!=0);
00245 return vnl_sse<T>::max(src,n);
00246 }
00247
00248
00249 template<class T>
00250 T vnl_c_vector<T>::min_value(T const *src, unsigned n)
00251 {
00252 assert(n!=0);
00253 return vnl_sse<T>::min(src,n);
00254 }
00255
00256
00257 template<class T>
00258 unsigned vnl_c_vector<T>::arg_max(T const *src, unsigned n)
00259 {
00260 assert(n!=0);
00261 return vnl_sse<T>::arg_max(src,n);
00262 }
00263
00264
00265 template<class T>
00266 unsigned vnl_c_vector<T>::arg_min(T const *src, unsigned n)
00267 {
00268 assert(n!=0);
00269 return vnl_sse<T>::arg_min(src,n);
00270 }
00271
00272
00273 template<class T>
00274 T vnl_c_vector<T>::euclid_dist_sq(T const *a, T const *b, unsigned n)
00275 {
00276 return vnl_sse<T>::euclid_dist_sq(a,b,n);
00277 }
00278
00279 template <class T>
00280 T vnl_c_vector<T>::sum_sq_diff_means(T const* v, unsigned n)
00281 {
00282 T sum(0);
00283 T sum_sq(0);
00284 for (unsigned i = 0; i < n; ++i, ++v)
00285 {
00286 sum += *v;
00287 sum_sq += *v * *v;
00288 }
00289 typedef typename vnl_numeric_traits<T>::abs_t abs_t;
00290 return sum_sq - sum*sum / abs_t(n);
00291 }
00292
00293
00294
00295 template <class T, class S>
00296 void vnl_c_vector_two_norm_squared(T const *p, unsigned n, S *out)
00297 {
00298 #if 1
00299
00300
00301
00302 S val = S(0);
00303 T const* end = p+n;
00304 while (p != end)
00305 val += S(vnl_math_squared_magnitude(*p++));
00306 *out = val;
00307 #else
00308 *out = 0;
00309 for (unsigned i=0; i<n; ++i)
00310 *out += vnl_math_squared_magnitude(p[i]);
00311 #endif
00312 }
00313
00314 template <class T, class S>
00315 void vnl_c_vector_rms_norm(T const *p, unsigned n, S *out)
00316 {
00317 vnl_c_vector_two_norm_squared(p, n, out);
00318 *out /= n;
00319 typedef typename vnl_numeric_traits<S>::real_t real_t;
00320 *out = S(vcl_sqrt(real_t(*out)));
00321 }
00322
00323 template <class T, class S>
00324 void vnl_c_vector_one_norm(T const *p, unsigned n, S *out)
00325 {
00326 *out = 0;
00327 T const* end = p+n;
00328 while (p != end)
00329 *out += vnl_math_abs(*p++);
00330 }
00331
00332 template <class T, class S>
00333 void vnl_c_vector_two_norm(T const *p, unsigned n, S *out)
00334 {
00335 vnl_c_vector_two_norm_squared(p, n, out);
00336 typedef typename vnl_numeric_traits<S>::real_t real_t;
00337 *out = S(vcl_sqrt(real_t(*out)));
00338 }
00339
00340 template <class T, class S>
00341 void vnl_c_vector_inf_norm(T const *p, unsigned n, S *out)
00342 {
00343 *out = 0;
00344 T const* end = p+n;
00345 while (p != end) {
00346 S v = vnl_math_abs(*p++);
00347 if (v > *out)
00348 *out = v;
00349 }
00350 }
00351
00352
00353
00354
00355
00356 inline void* vnl_c_vector_alloc(vcl_size_t n, unsigned size)
00357 {
00358 return vnl_sse_alloc(n,size);
00359 }
00360
00361
00362 inline void vnl_c_vector_dealloc(void* v, vcl_size_t n, unsigned size)
00363 {
00364 vnl_sse_dealloc(v,n,size);
00365 }
00366
00367 template<class T>
00368 T** vnl_c_vector<T>::allocate_Tptr(vcl_size_t n)
00369 {
00370 return (T**)vnl_c_vector_alloc(n, sizeof (T*));
00371 }
00372
00373 template<class T>
00374 void vnl_c_vector<T>::deallocate(T** v, vcl_size_t n)
00375 {
00376 vnl_c_vector_dealloc(v, n, sizeof (T*));
00377 }
00378
00379
00380 #include <vcl_new.h>
00381 template <class T> inline void vnl_c_vector_construct(T *p, vcl_size_t n)
00382 {
00383 for (vcl_size_t i=0; i<n; ++i)
00384 new (p+i) T();
00385 }
00386
00387 inline void vnl_c_vector_construct(float *, vcl_size_t) { }
00388 inline void vnl_c_vector_construct(double *, vcl_size_t) { }
00389 inline void vnl_c_vector_construct(long double *, vcl_size_t) { }
00390 inline void vnl_c_vector_construct(vcl_complex<float> *, vcl_size_t) { }
00391 inline void vnl_c_vector_construct(vcl_complex<double> *, vcl_size_t) { }
00392 inline void vnl_c_vector_construct(vcl_complex<long double> *, vcl_size_t) { }
00393
00394 #ifdef __BORLANDC__
00395
00396 # pragma option push -w-8057
00397
00398
00399
00400 #endif
00401
00402
00403 template <class T> inline void vnl_c_vector_destruct(T *p, vcl_size_t n)
00404 {
00405 for (vcl_size_t i=0; i<n; ++i)
00406 (p+i)->~T();
00407 }
00408
00409 #ifdef __BORLANDC__
00410 # pragma option pop
00411 #endif
00412
00413
00414 inline void vnl_c_vector_destruct(float *, vcl_size_t) { }
00415 inline void vnl_c_vector_destruct(double *, vcl_size_t) { }
00416 inline void vnl_c_vector_destruct(long double *, vcl_size_t) { }
00417 inline void vnl_c_vector_destruct(vcl_complex<float> *, vcl_size_t) { }
00418 inline void vnl_c_vector_destruct(vcl_complex<double> *, vcl_size_t) { }
00419 inline void vnl_c_vector_destruct(vcl_complex<long double> *, vcl_size_t) { }
00420
00421 template<class T>
00422 T* vnl_c_vector<T>::allocate_T(vcl_size_t n)
00423 {
00424 T *p = (T*)vnl_c_vector_alloc(n, sizeof (T));
00425 vnl_c_vector_construct(p, n);
00426 return p;
00427 }
00428
00429 template<class T>
00430 void vnl_c_vector<T>::deallocate(T* p, vcl_size_t n)
00431 {
00432 vnl_c_vector_destruct(p, n);
00433 vnl_c_vector_dealloc(p, n, sizeof (T));
00434 }
00435
00436 template<class T>
00437 vcl_ostream& print_vector(vcl_ostream& s, T const* v, unsigned size)
00438 {
00439 if (size != 0) s << v[0];
00440 for (unsigned i = 1; i < size; ++i)
00441 s << ' ' << v[i];
00442 return s;
00443 }
00444
00445
00446
00447 #define VNL_C_VECTOR_INSTANTIATE_norm(T, S) \
00448 template void vnl_c_vector_two_norm_squared(T const *, unsigned, S *); \
00449 template void vnl_c_vector_rms_norm(T const *, unsigned, S *); \
00450 template void vnl_c_vector_one_norm(T const *, unsigned, S *); \
00451 template void vnl_c_vector_two_norm(T const *, unsigned, S *); \
00452 template void vnl_c_vector_inf_norm(T const *, unsigned, S *)
00453
00454 #undef VNL_C_VECTOR_INSTANTIATE_ordered
00455 #define VNL_C_VECTOR_INSTANTIATE_ordered(T) \
00456 VNL_C_VECTOR_INSTANTIATE_norm(T, vnl_c_vector<T >::abs_t); \
00457 template class vnl_c_vector<T >; \
00458 template vcl_ostream& print_vector(vcl_ostream &,T const *,unsigned)
00459
00460 #undef VNL_C_VECTOR_INSTANTIATE_unordered
00461 #define VNL_C_VECTOR_INSTANTIATE_unordered(T) \
00462 VCL_DO_NOT_INSTANTIATE(T vnl_c_vector<T >::max_value(T const *, unsigned), T(0)); \
00463 VCL_DO_NOT_INSTANTIATE(T vnl_c_vector<T >::min_value(T const *, unsigned), T(0)); \
00464 VCL_DO_NOT_INSTANTIATE(unsigned vnl_c_vector<T >::arg_max(T const *, unsigned), 0U); \
00465 VCL_DO_NOT_INSTANTIATE(unsigned vnl_c_vector<T >::arg_min(T const *, unsigned), 0U); \
00466 VNL_C_VECTOR_INSTANTIATE_norm(T, vnl_c_vector<T >::abs_t); \
00467 template class vnl_c_vector<T >; \
00468 VCL_UNINSTANTIATE_SPECIALIZATION(T vnl_c_vector<T >::max_value(T const *, unsigned)); \
00469 VCL_UNINSTANTIATE_SPECIALIZATION(T vnl_c_vector<T >::min_value(T const *, unsigned)); \
00470 VCL_UNINSTANTIATE_SPECIALIZATION(unsigned vnl_c_vector<T >::arg_max(T const *, unsigned)); \
00471 VCL_UNINSTANTIATE_SPECIALIZATION(unsigned vnl_c_vector<T >::arg_min(T const *, unsigned))
00472
00473 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00474 #undef VNL_C_VECTOR_INSTANTIATE
00475 #define VNL_C_VECTOR_INSTANTIATE(T) extern "no such macro; use e.g. VNL_C_VECTOR_INSTANTIATE_ordered instead"
00476 #endif // DOXYGEN_SHOULD_SKIP_THIS
00477
00478 #endif // vnl_c_vector_txx_