loops_vctops.hpp

Go to the documentation of this file.

/********************************************************************
 *   This file is part of the source code of the realroot library.
 *   Author(s): J.P. Pavone, GALAAD, INRIA
 *   $Id: vctops.hpp,v 1.1 2005/07/11 10:03:56 jppavone Exp $
 ********************************************************************/
#ifndef realroot_SOLVE_SBDSLV_LOOPS_VCTOPS_H
#define realroot_SOLVE_SBDSLV_LOOPS_VCTOPS_H
//--------------------------------------------------------------------
#include <iostream>
#include <stdlib.h>
#include <realroot/system_epsilon.hpp>
//--------------------------------------------------------------------
namespace mmx {
//--------------------------------------------------------------------
namespace vctops
  {
    /* SCalarSUBstract (in place) */
    template<typename real_t> inline
    void scsub( const real_t& sc, real_t * data, int sz, int st = 1)
    { for ( int p = 0; p !=sz*st; p += st ) data[p] -= sc;  };

    template<typename real_t> inline
    void scsub( const real_t& sc, real_t * data, const real_t * src, int sz, int sta = 1, int stb = 1 )
    { int pd,p; for ( pd = p = 0; p !=sz*sta; p += sta, pd += stb ) data[pd] = src[p]-sc; };
  
    template<typename real_t> inline
    void accscmul( real_t * dst, const real_t& sc, real_t const * _src_, int sz, int stdst = 1, int stsrc = 1 )
    { 
      real_t const * src = _src_;
      while ( src != _src_ + sz*stsrc ) { *dst += *src*sc; src += stsrc; dst += stdst; };
    };
    
    /* SCalarMULtiply (in place) */
    template<typename real_t> inline
    void scmul( const real_t& sc, real_t * data, int sz, int st = 1)
    { for ( int p = 0; p != sz*st; p += st ) data[p] *= sc; };
    /* SCalarDIVide   (in place) */
    template<typename real_t> inline
    void scdiv( const real_t& sc, real_t * data, int sz, int st = 1)
    { for ( int p = 0; p != sz*st; p += st ) data[p] /= sc; };
    
    /* DOTPRODuct */
    template<typename real_t> inline 
    real_t dotprod( const real_t * a, const real_t * b, int sz, int sta = 1, int stb = 1)
    { 
      real_t acc = 0.0; 
      for ( int p = 0; p != sz*sta; p += sta, b += stb ) acc += a[p]*(*b); 
      return acc; 
    };
  
    template<typename real_t> inline /* ParallelADD */
    void padd( real_t * a, const real_t * b, unsigned sz, int sta = 1, int stb = 1  )
    { for ( int p = 0; p != sz*sta; p += sta, b += stb ) a[p] += *b; };
    
    template<typename real_t> inline /* ParallelSUB */
    void psub( real_t * a, const real_t * b, unsigned sz, int sta = 1, int stb = 1  )
    { for ( int p = 0; p != sz*sta; p += sta, b += stb ) a[p] -= *b; };
    
    template<typename real_t> inline /* ParallelMUL */
    void pmul( real_t * a, const real_t * b, unsigned sz, int sta = 1, int stb = 1  )
    { for ( int p = 0; p != sz*sta; p += sta, b += stb ) a[p] *= *b; };
    
    template<typename real_t> inline /* ParallelDIV */
    void pdiv( real_t * a, const real_t * b, unsigned sz, int sta = 1, int stb = 1  )
    { int p; for ( p  = 0; p != sz*sta; a[p] /= *b, p += sta, b += stb ) ; };
    
    template<typename real0, typename real_t> /* MIN and MAX Update */
    void minmaxu( real0& min, real0& max, real_t * src, int sz, int st=1 )
    {
      int i,p;
      for ( p = 0, i = 0; i < sz/2; i ++, p += 2*st )
    {
      if ( src[p] < src[p+st] ) 
        {
          if ( src[p]    < min ) min = src[p];
          if ( src[p+st] > max ) max = src[p+st];
        }
      else 
        {
          if ( src[p]    > max ) max = src[p];
          if ( src[p+st] < min ) min = src[p+st];
        };
    };
    };
  
    template<typename real0, typename real_t> /* MIN and MAX */
    void minmax( real0& min, real0& max, real_t * src, int sz, int st=1 )
    {
      min = src[(sz-1)*st];
      max = src[(sz-1)*st];
      minmaxu(min,max,src,sz,st);
    };
    
    template<typename real_t>
    real_t min( real_t const * src, int sz, int st = 1 )
    {
      real_t r = src[0];
      for ( int p = st; p != sz*st; p += st ) if ( r > src[p] ) r = src[p];    
      return r;
    };
    
    template<typename real_t>
    real_t max( real_t const * src, int sz, int st = 1 )
    {
      real_t r = src[0];
      for ( int p = st; p != sz*st; p += st ) if ( r < src[p] ) r = src[p];    
      return r;
    };
    
    
    template<typename T> /* Uniform RANDom, store sz random values in [a,b] to data  */
    void urand( T * data, unsigned sz, const T& a, const T& b, int st = 1 )
    {
      for ( int p = 0; p != sz*st; p += st )
    {
      double t = ((double)rand()/(double)RAND_MAX);
      data[p] = T(((1.0-t)*a + t*b));
    };
    };
    
    template<typename real_t> /* MEAN value */
    real_t mean( real_t const * const data, int sz, int st )
    {
        real_t mv   = 0.0;
        for ( int p = 0; p != sz*st; mv += data[p], p += st ) ;
        return mv/sz;
    }
    
    template<typename real_t>
    std::ostream& print( real_t const * const data, unsigned sz, int st = 1, std::ostream& out = std::cout )
    {
      int p = 0;
      out << "[ ";
      for ( unsigned i = 0; i < sz-1; i++, p += st )
    {
      out << data[p];
      out << ", ";
    };
      out << data[p] << " ]";
      return out;
    };
    
    template<typename real_t>  inline /* SiGN CouNT, count changes of sign */
    unsigned sgncnt ( real_t const *  b , unsigned sz, int st = 1 )
    {
      real_t const * last = b + (sz-1)*st;
      unsigned c;
      for ( c = 0; b != last; b += st )
    c += *b > 0 != *(b+st) > 0;
      return c;
    };
    
    template<typename real_t> inline /* SiGN CHanGes, changes of sign ? */
    bool sgnchg( real_t const * const  b, unsigned sz, int st = 1 )
    {
      int p = st;
      if ( std::abs(b[0]) < numerics::epsilon<real_t>::result || 
       std::abs(b[(sz-1)*st]) < numerics::epsilon<real_t>::result ) return true;
      bool pprv = b[0]>0;//numerics::epsilon<real_t>::result;
      //      bool pcurr;
      for ( unsigned i = 1; i < sz; i++, p+= st )  
    if ( (b[p]>0) != pprv )  return true;
      return false;
    };
    
    template<typename real_t> inline /* squared euclidian distance */
    real_t distance2( real_t const * const a, real_t const * const  b, unsigned sz, int sta, int stb )
    { 
      real_t s; 
      int pa,pb;
      for ( s = pb = pa = 0; pa != sz*sta; pa += sta, pb += stb ) s += (a[pa]-b[pb])*(a[pa]-b[pb]); 
      return s; 
    };
    
    template<typename real_t> inline
    real_t delta_max( real_t * a, real_t * b, unsigned sz, int st = 1 )
    {
      int p = st;
      real_t tmp = a[0]-b[0];
      for ( unsigned i = 1; i < sz; i++ , p += st )
    if ( a[p]-b[p] > tmp ) tmp = a[p]-b[p];
      return tmp;
    };
    
    template<typename real_t> inline
    int delta_max_index( real_t * a, real_t * b, unsigned sz, int st = 1  )
    {
      int mx = 0;
      int p = st;
      real_t tmp = a[0]-b[0];
      for ( unsigned i = 1; i < sz; i++ , p += st )
    if ( a[p]-b[p] > tmp ) { mx = i; tmp = a[p]-b[p]; };
      return mx;
    };
    
    template<typename real_t> inline
    real_t sum( real_t const * const src, unsigned sz, int st = 1 )
    {  
      real_t acc = (real_t)0.0; 
      for ( int i = 0; i != sz*st; acc += src[i], i += st ) 
        ; 
      return acc;
    } 
    
    template<typename real_t> inline
    unsigned set_conversion( real_t * src, unsigned sz, const real_t& epsilon = (real_t)0 )
    {
      std::sort(src,src+sz);
      unsigned c = 0;
      for ( unsigned i = 1; i < sz; i++ )
    if ( src[i]-src[c] > epsilon ) src[++c] = src[i];
      return c+1;
    };
    
    template<typename real_t> inline 
    void scale( real_t * src, unsigned sz, const real_t& sc = (real_t)(1.0), int st= 1 )
    { real_t mn,mx; minmax(mn,mx,src,sz,st); scmul(sc/(mx-mn),src,sz,st); };
  
    template< typename U, typename Y, typename Z  > /* convolution 1d */
    void convolution(U*dst,Y const*a,Z const*b, unsigned sza, unsigned szb, int sta=1, int stb=1, int stout=1)
    {
      assert((dst!=a)&&(dst!=b));
      int p,ia,ib,q;
      for ( p = 0; p != (sza+szb-1)*stout; dst[p] = (U)0, p += stout ) ;
      for ( p = ia =0; ia != sza*sta; ia += sta, p += stout ) 
    for ( ib = 0, q = p; ib != szb*stb; dst[q] += a[ia]*b[ib], ib += stb, q += stout ) ;
    };
    
    template< typename U, typename Y, typename Z  >
    void self_convolution( U * dst, Y const * a, unsigned sza,  int sta = 1,  int stout = 1  )
    {
      int p,q,ia,ib;
      for ( p = 0; p != (2*sza-1); dst[p] = 0, p += stout ) ;
      for ( p = ia = 0; ia != sza*sta; dst[2*p] += a[ia]*a[ia], ia += sta, p += stout )
    for ( ib = 0, q = p; ib != ia; dst[q] += 2*a[ia]*a[ib], ib += sta, q += stout ) ;
    };
  };
//--------------------------------------------------------------------
} //namespace mmx
/********************************************************************/
#endif //