|
realroot_doc 0.1.1
|
|
realroot_doc 0.1.1
|
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 int i,p; 00090 min = src[(sz-1)*st]; 00091 max = src[(sz-1)*st]; 00092 minmaxu(min,max,src,sz,st); 00093 }; 00094 00095 template<typename real_t> 00096 real_t min( real_t const * src, int sz, int st = 1 ) 00097 { 00098 real_t r = src[0]; 00099 for ( int p = st; p != sz*st; p += st ) if ( r > src[p] ) r = src[p]; 00100 return r; 00101 }; 00102 00103 template<typename real_t> 00104 real_t max( real_t const * src, int sz, int st = 1 ) 00105 { 00106 real_t r = src[0]; 00107 for ( int p = st; p != sz*st; p += st ) if ( r < src[p] ) r = src[p]; 00108 return r; 00109 }; 00110 00111 00112 template<typename T> /* Uniform RANDom, store sz random values in [a,b] to data */ 00113 void urand( T * data, unsigned sz, const T& a, const T& b, int st = 1 ) 00114 { 00115 for ( int p = 0; p != sz*st; p += st ) 00116 { 00117 double t = ((double)rand()/(double)RAND_MAX); 00118 data[p] = T(((1.0-t)*a + t*b)); 00119 }; 00120 }; 00121 00122 template<typename real_t> /* MEAN value */ 00123 real_t mean( real_t const * const data, int sz, int st ) 00124 { real_t mv = 0.0; for ( int p = 0; p != sz*st; mv += data[p], p += st ) ; return mv/sz; }; 00125 00126 template<typename real_t> 00127 std::ostream& print( real_t const * const data, unsigned sz, int st = 1, std::ostream& out = std::cout ) 00128 { 00129 int p = 0; 00130 out << "[ "; 00131 for ( unsigned i = 0; i < sz-1; i++, p += st ) 00132 { 00133 out << data[p]; 00134 out << ", "; 00135 }; 00136 out << data[p] << " ]"; 00137 return out; 00138 }; 00139 00140 template<typename real_t> inline /* SiGN CouNT, count changes of sign */ 00141 unsigned sgncnt ( real_t const * b , unsigned sz, int st = 1 ) 00142 { 00143 real_t const * last = b + (sz-1)*st; 00144 unsigned c; 00145 for ( c = 0; b != last; b += st ) 00146 c += *b > 0 != *(b+st) > 0; 00147 return c; 00148 }; 00149 00150 template<typename real_t> inline /* SiGN CHanGes, changes of sign ? */ 00151 bool sgnchg( real_t const * const b, unsigned sz, int st = 1 ) 00152 { 00153 int p = st; 00154 if ( std::abs(b[0]) < numerics::epsilon<real_t>::result || 00155 std::abs(b[(sz-1)*st]) < numerics::epsilon<real_t>::result ) return true; 00156 bool pprv = b[0]>0;//numerics::epsilon<real_t>::result; 00157 // bool pcurr; 00158 for ( unsigned i = 1; i < sz; i++, p+= st ) 00159 if ( (b[p]>0) != pprv ) return true; 00160 return false; 00161 }; 00162 00163 template<typename real_t> inline /* squared euclidian distance */ 00164 real_t distance2( real_t const * const a, real_t const * const b, unsigned sz, int sta, int stb ) 00165 { 00166 real_t s; 00167 int pa,pb; 00168 for ( s = pb = pa = 0; pa != sz*sta; pa += sta, pb += stb ) s += (a[pa]-b[pb])*(a[pa]-b[pb]); 00169 return s; 00170 }; 00171 00172 template<typename real_t> inline 00173 real_t delta_max( real_t * a, real_t * b, unsigned sz, int st = 1 ) 00174 { 00175 int p = st; 00176 real_t tmp = a[0]-b[0]; 00177 for ( unsigned i = 1; i < sz; i++ , p += st ) 00178 if ( a[p]-b[p] > tmp ) tmp = a[p]-b[p]; 00179 return tmp; 00180 }; 00181 00182 template<typename real_t> inline 00183 int delta_max_index( real_t * a, real_t * b, unsigned sz, int st = 1 ) 00184 { 00185 int mx = 0; 00186 int p = st; 00187 real_t tmp = a[0]-b[0]; 00188 for ( unsigned i = 1; i < sz; i++ , p += st ) 00189 if ( a[p]-b[p] > tmp ) { mx = i; tmp = a[p]-b[p]; }; 00190 return mx; 00191 }; 00192 00193 template<typename real_t> inline 00194 real_t sum( real_t const * const src, unsigned sz, int st = 1 ) 00195 { real_t acc = (real_t)0.0; for ( int i = 0; i != sz*st; acc += src[i], i += st ) ; return acc;}; 00196 00197 template<typename real_t> inline 00198 unsigned set_conversion( real_t * src, unsigned sz, const real_t& epsilon = (real_t)0 ) 00199 { 00200 std::sort(src,src+sz); 00201 unsigned c = 0; 00202 for ( unsigned i = 1; i < sz; i++ ) 00203 if ( src[i]-src[c] > epsilon ) src[++c] = src[i]; 00204 return c+1; 00205 }; 00206 00207 template<typename real_t> inline 00208 void scale( real_t * src, unsigned sz, const real_t& sc = (real_t)(1.0), int st= 1 ) 00209 { real_t mn,mx; minmax(mn,mx,src,sz,st); scmul(sc/(mx-mn),src,sz,st); }; 00210 00211 template< typename U, typename Y, typename Z > /* convolution 1d */ 00212 void convolution(U*dst,Y const*a,Z const*b, unsigned sza, unsigned szb, int sta=1, int stb=1, int stout=1) 00213 { 00214 assert((dst!=a)&&(dst!=b)); 00215 int p,ia,ib,q; 00216 for ( p = 0; p != (sza+szb-1)*stout; dst[p] = (U)0, p += stout ) ; 00217 for ( p = ia =0; ia != sza*sta; ia += sta, p += stout ) 00218 for ( ib = 0, q = p; ib != szb*stb; dst[q] += a[ia]*b[ib], ib += stb, q += stout ) ; 00219 }; 00220 00221 template< typename U, typename Y, typename Z > 00222 void self_convolution( U * dst, Y const * a, unsigned sza, int sta = 1, int stout = 1 ) 00223 { 00224 int p,q,ia,ib; 00225 for ( p = 0; p != (2*sza-1); dst[p] = 0, p += stout ) ; 00226 for ( p = ia = 0; ia != sza*sta; dst[2*p] += a[ia]*a[ia], ia += sta, p += stout ) 00227 for ( ib = 0, q = p; ib != ia; dst[q] += 2*a[ia]*a[ib], ib += sta, q += stout ) ; 00228 }; 00229 }; 00230 //-------------------------------------------------------------------- 00231 } //namespace mmx 00232 /********************************************************************/ 00233 #endif //
