Rtmra_Data_t_scalar.h

Go to the documentation of this file.

// Rtmra_Data_t_scalar.h
//////////////////////////////////////////////////////////////////////
//          (c) Copyright 2000-2002 Brown Deer Technology, LLC.
//                        All rights reserved.
//////////////////////////////////////////////////////////////////////

#ifndef _RTMRA_DATA_T_SCALAR_H
#define _RTMRA_DATA_T_SCALAR_H

#ifdef USE_STD
#include <cmath>
using namespace std;
#else
#include <math.h>
#endif

#include "stormdef.h"
#define INVSQRT2 STORMDEF_MATHCONST_INVSQRT2

// partial specialization Rtmra<Data_t,scalar>
template <typename Data_t>
class Rtmra<Data_t,scalar> {

   public:

      typedef Data_t DataType;
      static const int ArrayDim = 0;

      Rtmra( const int msize = 8);
      virtual ~Rtmra();
      virtual void remalloc( const int msize = 8);
      virtual void clear();
      void assign( const Data_t data ) { *_data = data; }
      Rtmra<Data_t,scalar>& operator=( const Data_t data);
      operator Data_t() { return *_data; }
      virtual void rtwt();
      int operator ++() { rtwt(); return 0; }
      int operator ++(int) { rtwt(); return 0; }
      int mmax() { return _mmax; }
      virtual void rtiwt();
      virtual void rtjiwt();
      int operator --() { rtiwt(); return 0; }
      int operator --(int) { rtiwt(); return 0; }
      int mrastep() { return _mrastep; }
      void mrastep( const int mrastep ) { _mrastep = mrastep; } //DEBUG
      virtual int mcycle();
      bool isvalid(const int m) { return ( m < mcycle() ) ; }
      virtual Data_t scal( const int ); // m
      virtual Data_t wav( const int );  // m
      virtual Data_t ave( const int );  // m
      virtual Data_t chg( const int );  // m
      virtual int imcycle();
      virtual void scal( const int m, Data_t scalcoef );
      virtual void wav( const int m, Data_t wavcoef );

   protected:
      int _Msize;
      int _mmax;
      int _mrastep;
      Data_t *_mracoef;
      Data_t *_data;

};

template <typename Data_t>
inline Rtmra<Data_t,scalar>::Rtmra( const int msize ) 
  : _Msize(msize), _mmax(msize-1), _mrastep(-1) 
{
   _mracoef = new Data_t[ 2*_Msize ];
   _data = new Data_t;
   for( int i = 0; i < 2*_Msize; i++) { _mracoef[i] = 0; }
   *_data = 0;
}

template <typename Data_t>
inline Rtmra<Data_t,scalar>::~Rtmra() 
{
   delete _data; _data = 0;
   delete [] _mracoef; _mracoef = 0;
}

template <typename Data_t>
inline void Rtmra<Data_t,scalar>::remalloc( const int msize ) 
{
   _Msize = msize; 
   _mmax = msize-1;
   _mrastep = -1;
   if (_mracoef) delete [] _mracoef; _mracoef = 0;
   if (_data) delete _data; _data = 0;
   _mracoef = new Data_t[ 2*_Msize ];
   _data = new Data_t;
   for( int i = 0; i< 2*_Msize; i++) { _mracoef[i] = 0; }
   *_data = 0;
}

template <typename Data_t>
inline void Rtmra<Data_t,scalar>::clear() 
{
   _mrastep = -1;
   for( int i = 0; i < 2*_Msize; i++) { _mracoef[i] = 0; }
   *_data = 0;
}

template <typename Data_t>
inline Rtmra<Data_t,scalar>&
  Rtmra<Data_t,scalar>::operator=( const Data_t data ) 
{
   assign(data);
   return *this;
}

template <typename Data_t>
inline void Rtmra<Data_t,scalar>::rtwt() 
{
#ifdef STORM_USEC
   rtmrartwt( _Msize, _mmax, _mrastep, _mracoef, *_data );
#else
   _mrastep++;
   int mcycle = 0;
   while ( ( _mrastep & (1<<mcycle) ) && ( mcycle < _mmax ) ) { mcycle++; }
   Data_t tmpold;
   Data_t tmp = *_data;
   int m = 1;
   while ( m <= mcycle ) {
      tmpold = _mracoef[ 2*(m-1) ];
      _mracoef[ 2*(m-1) ] = tmp;
      _mracoef[ 2*m + 1 ] = (tmpold-tmp) * INVSQRT2;
      tmp = (tmpold+tmp) * INVSQRT2;
      m++;
   }
   _mracoef[ 2*(m-1) ] = tmp;
#endif
}

template <typename Data_t>
inline void Rtmra<Data_t,scalar>::rtiwt() 
{
#ifdef STORM_USEC
   rtmrartiwt( _Msize, _mmax, _mrastep, _mracoef, *_data );
#else
   _mrastep++;
   int itmp = (1<<_mmax) - 1;
   itmp = itmp - ( itmp & _mrastep );
   int imcycle = 0;
   while ( ( itmp & (1<<imcycle) ) && ( imcycle < _mmax ) ) { imcycle++; }
   int m = imcycle;
   if ( m < _mmax ) { 
      _mracoef[2*m] = ( _mracoef[ 2*(m+1) ] 
                        - _mracoef[ 2*(m+1) + 1 ] ) * INVSQRT2; 
   }
   m--;
   while ( m >= 0 ) {
      _mracoef[2*m] = ( _mracoef[ 2*(m+1) ] 
                        + _mracoef[ 2*(m+1) + 1 ]) * INVSQRT2;
      m--;
   }
   *_data = _mracoef[0];
#endif
}

template <typename Data_t>
inline void Rtmra<Data_t,scalar>::rtjiwt() 
{
//#ifdef STORM_USEC
   rtmrartjiwt( _Msize, _mmax, _mrastep, _mracoef, *_data );
//#else
//   _mrastep++;
//   int itmp = (1<<_mmax) - 1;
//   itmp = itmp - ( itmp & _mrastep );
//   int imcycle = 0;
//   while ( ( itmp & (1<<imcycle) ) && ( imcycle < _mmax ) ) { imcycle++; }
//   int m = imcycle;
//   if ( m < _mmax ) { 
//      _mracoef[2*m] = ( _mracoef[ 2*(m+1) ] 
//                        - _mracoef[ 2*(m+1) + 1 ] ) * INVSQRT2; 
//   }
//   m--;
//   while ( m >= 0 ) {
//      _mracoef[2*m] = ( _mracoef[ 2*(m+1) ] 
//                        + _mracoef[ 2*(m+1) + 1 ]) * INVSQRT2;
//      m--;
//   }
//   *_data = _mracoef[0];
//#endif
}

template <typename Data_t>
inline int Rtmra<Data_t,scalar>::mcycle() 
{
   int mcycle=0;
#ifdef STORM_USEC
   rtmramcycle( _Msize, _mmax, _mrastep, mcycle );
#else
   while ((_mrastep&(1<<mcycle))&&(mcycle<_mmax)) { mcycle++; }
#endif
   return mcycle;
}

template <typename Data_t>
inline Data_t Rtmra<Data_t,scalar>::scal( const int m ) 
{
   Data_t scalcoef;
#ifdef STORM_USEC
   rtmrascalcoef( _Msize, _mmax, _mrastep, _mracoef, m, scalcoef );
#else
//   if ( ( m >= 0 ) && ( m <= _mmax ) ) { 
      scalcoef = _mracoef[2*m]; 
//   }
//   else { 
//      scalcoef = 1.0/0; 
//   }
#endif
   return scalcoef;
}

template <typename Data_t>
inline Data_t Rtmra<Data_t,scalar>::wav( const int m ) 
{
   Data_t wavcoef;
#ifdef STORM_USEC
   rtmrawavcoef( _Msize, _mmax, _mrastep, _mracoef, m, wavcoef );
#else
//   if ( ( m >= 0 ) && ( m <= _mmax ) ) { 
      wavcoef = _mracoef[ 2*m + 1 ]; 
//   }
//   else { wavcoef = 1.0/0; }
#endif
   return wavcoef;
}

template <typename Data_t>
inline Data_t Rtmra<Data_t,scalar>::ave( const int m ) 
{
   Data_t scala;
#ifdef STORM_USEC
   rtmrascala( _Msize, _mmax, _mrastep, _mracoef, m, scala );
#else
//   if ( ( m >= 0 ) && ( m <= _mmax ) ) { 
      scala = _mracoef[2*m] * pow(INVSQRT2,m); 
//   }
//   else { 
//      scala = 1.0/0; 
//   }
#endif
   return scala;
}

template <typename Data_t>
inline Data_t Rtmra<Data_t,scalar>::chg( const int m ) 
{
   Data_t wavd;
#ifdef STORM_USEC
   rtmrawavd( _Msize, _mmax, _mrastep, _mracoef, m, wavd );
#else
//   if ( ( m >= 0 ) && ( m <= _mmax ) ) { 
      wavd = _mracoef[ 2*m + 1 ] * pow(INVSQRT2,m); 
//   }
//   else { 
//      wavd = 1.0/0; 
//   }
#endif
   return wavd;
}

template <typename Data_t>
inline int Rtmra<Data_t,scalar>::imcycle() 
{
   int imcycle = 0;
#ifdef STORM_USEC
   rtmraimcycle( _Msize, _mmax, _mrastep, imcycle );
#else
   int itmp = (1<<_mmax) - 1;
   itmp = itmp - ( itmp & ( _mrastep + 1) );
   while ( ( itmp & (1<<imcycle) ) && ( imcycle < _mmax ) ) { imcycle++; }
#endif
   return imcycle;
}

template <typename Data_t>
inline void Rtmra<Data_t,scalar>::scal( const int m, Data_t scalcoef ) 
{
#ifdef STORM_USEC
   rtmraputscalcoef( _Msize, _mmax, _mrastep, _mracoef, m, scalcoef );
#else
//   if ( ( m >= 0 ) && ( m <= _mmax ) ) { 
   _mracoef[2*m] = scalcoef; 
//    }
#endif
}

template <typename Data_t>
inline void Rtmra<Data_t,scalar>::wav( const int m, Data_t wavcoef ) 
{
#ifdef STORM_USEC
   rtmraputwavcoef( _Msize, _mmax, _mrastep, _mracoef, m, wavcoef );
#else
//   if ( ( m >= 0 ) && ( m <= _mmax ) ) { 
   _mracoef[ 2*m + 1 ] = wavcoef; 
//   }
#endif
}

template <>
inline void Rtmra<int,scalar>::rtwt() 
{
#ifdef STORM_USEC
   rtmrartwt( _Msize, _mmax, _mrastep, _mracoef, *_data );
#else
   _mrastep++;
   int mcycle = 0;
   while ( ( _mrastep & (1<<mcycle) ) && ( mcycle < _mmax ) ) { mcycle++; }
   int tmpold;
   int tmp = *_data;
   int m = 1;
   while ( m <= mcycle ) {
      tmpold = _mracoef[ 2*(m-1) ];
      _mracoef[ 2*(m-1) ] = tmp;
      _mracoef[ 2*m + 1 ] = (tmpold-tmp);
      tmp = (tmpold+tmp) >> 1;
      m++;
   }
   _mracoef[ 2*(m-1) ] = tmp;
#endif
}

template <>
inline void Rtmra<int,scalar>::rtiwt() 
{
#ifdef STORM_USEC
   rtmrartiwt( _Msize, _mmax, _mrastep, _mracoef, *_data );
//cout << "SPECIAL" << endl;
#else
   _mrastep++;
   int itmp = (1<<_mmax) - 1;
   itmp = itmp - ( itmp & _mrastep );
   int imcycle = 0;
   while ( ( itmp & (1<<imcycle) ) && ( imcycle < _mmax ) ) { imcycle++; }
   int m = imcycle;
   if ( m < _mmax ) { 
      _mracoef[2*m] = _mracoef[ 2*(m+1) ] 
                  - ( _mracoef[ 2*(m+1) + 1 ] >> 1 );
   }
   m--;
   while ( m >= 0 ) {
      _mracoef[2*m] = _mracoef[ 2*(m+1) ] 
                        + ( ( _mracoef[ 2*(m+1) + 1 ] + 1 ) >> 1 );
      m--;
   }
   *_data = _mracoef[0];
#endif
}

#endif // _RTMRA_DATA_T_SCALAR_H

---