realroot: include/realroot/loops_vctops.hpp Source File

00001 /********************************************************************
00002  *   This file is part of the source code of the realroot library.
00003  *   Author(s): J.P. Pavone, GALAAD, INRIA
00004  *   $Id: vctops.hpp,v 1.1 2005/07/11 10:03:56 jppavone Exp $
00005  ********************************************************************/
00006 #ifndef realroot_SOLVE_SBDSLV_LOOPS_VCTOPS_H
00007 #define realroot_SOLVE_SBDSLV_LOOPS_VCTOPS_H
00008 //--------------------------------------------------------------------
00009 #include <iostream>
00010 #include <stdlib.h>
00011 #include <realroot/system_epsilon.hpp>
00012 //--------------------------------------------------------------------
00013 namespace mmx {
00014 //--------------------------------------------------------------------
00015 namespace vctops
00016   {
00017     /* SCalarSUBstract (in place) */
00018     template<typename real_t> inline
00019     void scsub( const real_t& sc, real_t * data, int sz, int st = 1)
00020     { for ( int p = 0; p !=sz*st; p += st ) data[p] -= sc;  };
00021 
00022     template<typename real_t> inline
00023     void scsub( const real_t& sc, real_t * data, const real_t * src, int sz, int sta = 1, int stb = 1 )
00024     { int pd,p; for ( pd = p = 0; p !=sz*sta; p += sta, pd += stb ) data[pd] = src[p]-sc; };
00025   
00026     template<typename real_t> inline
00027     void accscmul( real_t * dst, const real_t& sc, real_t const * _src_, int sz, int stdst = 1, int stsrc = 1 )
00028     { 
00029       real_t const * src = _src_;
00030       while ( src != _src_ + sz*stsrc ) { *dst += *src*sc; src += stsrc; dst += stdst; };
00031     };
00032     
00033     /* SCalarMULtiply (in place) */
00034     template<typename real_t> inline
00035     void scmul( const real_t& sc, real_t * data, int sz, int st = 1)
00036     { for ( int p = 0; p != sz*st; p += st ) data[p] *= sc; };
00037     /* SCalarDIVide   (in place) */
00038     template<typename real_t> inline
00039     void scdiv( const real_t& sc, real_t * data, int sz, int st = 1)
00040     { for ( int p = 0; p != sz*st; p += st ) data[p] /= sc; };
00041     
00042     /* DOTPRODuct */
00043     template<typename real_t> inline 
00044     real_t dotprod( const real_t * a, const real_t * b, int sz, int sta = 1, int stb = 1)
00045     { 
00046       real_t acc = 0.0; 
00047       for ( int p = 0; p != sz*sta; p += sta, b += stb ) acc += a[p]*(*b); 
00048       return acc; 
00049     };
00050   
00051     template<typename real_t> inline /* ParallelADD */
00052     void padd( real_t * a, const real_t * b, unsigned sz, int sta = 1, int stb = 1  )
00053     { for ( int p = 0; p != sz*sta; p += sta, b += stb ) a[p] += *b; };
00054     
00055     template<typename real_t> inline /* ParallelSUB */
00056     void psub( real_t * a, const real_t * b, unsigned sz, int sta = 1, int stb = 1  )
00057     { for ( int p = 0; p != sz*sta; p += sta, b += stb ) a[p] -= *b; };
00058     
00059     template<typename real_t> inline /* ParallelMUL */
00060     void pmul( real_t * a, const real_t * b, unsigned sz, int sta = 1, int stb = 1  )
00061     { for ( int p = 0; p != sz*sta; p += sta, b += stb ) a[p] *= *b; };
00062     
00063     template<typename real_t> inline /* ParallelDIV */
00064     void pdiv( real_t * a, const real_t * b, unsigned sz, int sta = 1, int stb = 1  )
00065     { int p; for ( p  = 0; p != sz*sta; a[p] /= *b, p += sta, b += stb ) ; };
00066     
00067     template<typename real0, typename real_t> /* MIN and MAX Update */
00068     void minmaxu( real0& min, real0& max, real_t * src, int sz, int st=1 )
00069     {
00070       int i,p;
00071       for ( p = 0, i = 0; i < sz/2; i ++, p += 2*st )
00072         {
00073           if ( src[p] < src[p+st] ) 
00074             {
00075               if ( src[p]    < min ) min = src[p];
00076               if ( src[p+st] > max ) max = src[p+st];
00077             }
00078           else 
00079             {
00080               if ( src[p]    > max ) max = src[p];
00081               if ( src[p+st] < min ) min = src[p+st];
00082             };
00083         };
00084     };
00085   
00086     template<typename real0, typename real_t> /* MIN and MAX */
00087     void minmax( real0& min, real0& max, real_t * src, int sz, int st=1 )
00088     {
00089       min = src[(sz-1)*st];
00090       max = src[(sz-1)*st];
00091       minmaxu(min,max,src,sz,st);
00092     };
00093     
00094     template<typename real_t>
00095     real_t min( real_t const * src, int sz, int st = 1 )
00096     {
00097       real_t r = src[0];
00098       for ( int p = st; p != sz*st; p += st ) if ( r > src[p] ) r = src[p];    
00099       return r;
00100     };
00101     
00102     template<typename real_t>
00103     real_t max( real_t const * src, int sz, int st = 1 )
00104     {
00105       real_t r = src[0];
00106       for ( int p = st; p != sz*st; p += st ) if ( r < src[p] ) r = src[p];    
00107       return r;
00108     };
00109     
00110     
00111     template<typename T> /* Uniform RANDom, store sz random values in [a,b] to data  */
00112     void urand( T * data, unsigned sz, const T& a, const T& b, int st = 1 )
00113     {
00114       for ( int p = 0; p != sz*st; p += st )
00115         {
00116           double t = ((double)rand()/(double)RAND_MAX);
00117           data[p] = T(((1.0-t)*a + t*b));
00118         };
00119     };
00120     
00121     template<typename real_t> /* MEAN value */
00122     real_t mean( real_t const * const data, int sz, int st )
00123     { real_t mv   = 0.0; for ( int p = 0; p != sz*st; mv += data[p], p += st ) ; return mv/sz; };
00124     
00125     template<typename real_t>
00126     std::ostream& print( real_t const * const data, unsigned sz, int st = 1, std::ostream& out = std::cout )
00127     {
00128       int p = 0;
00129       out << "[ ";
00130       for ( unsigned i = 0; i < sz-1; i++, p += st )
00131         {
00132           out << data[p];
00133           out << ", ";
00134         };
00135       out << data[p] << " ]";
00136       return out;
00137     };
00138     
00139     template<typename real_t>  inline /* SiGN CouNT, count changes of sign */
00140     unsigned sgncnt ( real_t const *  b , unsigned sz, int st = 1 )
00141     {
00142       real_t const * last = b + (sz-1)*st;
00143       unsigned c;
00144       for ( c = 0; b != last; b += st )
00145         c += *b > 0 != *(b+st) > 0;
00146       return c;
00147     };
00148     
00149     template<typename real_t> inline /* SiGN CHanGes, changes of sign ? */
00150     bool sgnchg( real_t const * const  b, unsigned sz, int st = 1 )
00151     {
00152       int p = st;
00153       if ( std::abs(b[0]) < numerics::epsilon<real_t>::result || 
00154            std::abs(b[(sz-1)*st]) < numerics::epsilon<real_t>::result ) return true;
00155       bool pprv = b[0]>0;//numerics::epsilon<real_t>::result;
00156       //      bool pcurr;
00157       for ( unsigned i = 1; i < sz; i++, p+= st )  
00158         if ( (b[p]>0) != pprv )  return true;
00159       return false;
00160     };
00161     
00162     template<typename real_t> inline /* squared euclidian distance */
00163     real_t distance2( real_t const * const a, real_t const * const  b, unsigned sz, int sta, int stb )
00164     { 
00165       real_t s; 
00166       int pa,pb;
00167       for ( s = pb = pa = 0; pa != sz*sta; pa += sta, pb += stb ) s += (a[pa]-b[pb])*(a[pa]-b[pb]); 
00168       return s; 
00169     };
00170     
00171     template<typename real_t> inline
00172     real_t delta_max( real_t * a, real_t * b, unsigned sz, int st = 1 )
00173     {
00174       int p = st;
00175       real_t tmp = a[0]-b[0];
00176       for ( unsigned i = 1; i < sz; i++ , p += st )
00177         if ( a[p]-b[p] > tmp ) tmp = a[p]-b[p];
00178       return tmp;
00179     };
00180     
00181     template<typename real_t> inline
00182     int delta_max_index( real_t * a, real_t * b, unsigned sz, int st = 1  )
00183     {
00184       int mx = 0;
00185       int p = st;
00186       real_t tmp = a[0]-b[0];
00187       for ( unsigned i = 1; i < sz; i++ , p += st )
00188         if ( a[p]-b[p] > tmp ) { mx = i; tmp = a[p]-b[p]; };
00189       return mx;
00190     };
00191     
00192     template<typename real_t> inline
00193     real_t sum( real_t const * const src, unsigned sz, int st = 1 )
00194     {  real_t acc = (real_t)0.0; for ( int i = 0; i != sz*st; acc += src[i], i += st ) ; return acc;}; 
00195     
00196     template<typename real_t> inline
00197     unsigned set_conversion( real_t * src, unsigned sz, const real_t& epsilon = (real_t)0 )
00198     {
00199       std::sort(src,src+sz);
00200       unsigned c = 0;
00201       for ( unsigned i = 1; i < sz; i++ )
00202         if ( src[i]-src[c] > epsilon ) src[++c] = src[i];
00203       return c+1;
00204     };
00205     
00206     template<typename real_t> inline 
00207     void scale( real_t * src, unsigned sz, const real_t& sc = (real_t)(1.0), int st= 1 )
00208     { real_t mn,mx; minmax(mn,mx,src,sz,st); scmul(sc/(mx-mn),src,sz,st); };
00209   
00210     template< typename U, typename Y, typename Z  > /* convolution 1d */
00211     void convolution(U*dst,Y const*a,Z const*b, unsigned sza, unsigned szb, int sta=1, int stb=1, int stout=1)
00212     {
00213       assert((dst!=a)&&(dst!=b));
00214       int p,ia,ib,q;
00215       for ( p = 0; p != (sza+szb-1)*stout; dst[p] = (U)0, p += stout ) ;
00216       for ( p = ia =0; ia != sza*sta; ia += sta, p += stout ) 
00217         for ( ib = 0, q = p; ib != szb*stb; dst[q] += a[ia]*b[ib], ib += stb, q += stout ) ;
00218     };
00219     
00220     template< typename U, typename Y, typename Z  >
00221     void self_convolution( U * dst, Y const * a, unsigned sza,  int sta = 1,  int stout = 1  )
00222     {
00223       int p,q,ia,ib;
00224       for ( p = 0; p != (2*sza-1); dst[p] = 0, p += stout ) ;
00225       for ( p = ia = 0; ia != sza*sta; dst[2*p] += a[ia]*a[ia], ia += sta, p += stout )
00226         for ( ib = 0, q = p; ib != ia; dst[q] += 2*a[ia]*a[ib], ib += sta, q += stout ) ;
00227     };
00228   };
00229 //--------------------------------------------------------------------
00230 } //namespace mmx
00231 /********************************************************************/
00232 #endif //