core/vnl/vnl_c_vector.txx
Go to the documentation of this file.
00001 // This is core/vnl/vnl_c_vector.txx
00002 #ifndef vnl_c_vector_txx_
00003 #define vnl_c_vector_txx_
00004 //:
00005 // \file
00006 // \author Andrew W. Fitzgibbon, Oxford RRG
00007 // \date   12 Feb 1998
00008 //
00009 //-----------------------------------------------------------------------------
00010 
00011 #include "vnl_c_vector.h"
00012 #include <vcl_cmath.h>     // vcl_sqrt()
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 // non-conjugating "dot" product.
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 // conjugating "dot" product.
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 // conjugates one block of data into another block.
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 //: Returns max value of the vector.
00241 template<class T>
00242 T vnl_c_vector<T>::max_value(T const *src, unsigned n)
00243 {
00244   assert(n!=0); // max_value of an empty vector is undefined
00245   return vnl_sse<T>::max(src,n);
00246 }
00247 
00248 //: Returns min value of the vector.
00249 template<class T>
00250 T vnl_c_vector<T>::min_value(T const *src, unsigned n)
00251 {
00252   assert(n!=0); // min_value of an empty vector is undefined
00253   return vnl_sse<T>::min(src,n);
00254 }
00255 
00256 //: Returns location of max value of the vector.
00257 template<class T>
00258 unsigned vnl_c_vector<T>::arg_max(T const *src, unsigned n)
00259 {
00260   assert(n!=0); // max value of an empty vector is undefined
00261   return vnl_sse<T>::arg_max(src,n);
00262 }
00263 
00264 //: Returns location of min value of the vector.
00265 template<class T>
00266 unsigned vnl_c_vector<T>::arg_min(T const *src, unsigned n)
00267 {
00268   assert(n!=0); // min value of an empty vector is undefined
00269   return vnl_sse<T>::arg_min(src,n);
00270 }
00271 
00272 //: Sum of Differences squared.
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   // IMS: MSVC's optimiser does much better with *p++ than with p[i];
00300   // consistently about 30% better over vectors from 4 to 20000 dimensions.
00301   // PVr: with gcc 3.0 on alpha this is even a factor 3 faster!
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 // "T *" is POD, but "T" might not be.
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 // The compiler is confused
00396 # pragma option push -w-8057
00397 // Warning W8057 vnl/vnl_c_vector.txx 414:
00398 // Parameter 'p' is never used in function
00399 // vnl_c_vector_destruct<int>(int *,int)
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)   // For each index in vector
00441     s << ' ' << v[i];                   // Output data element
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_