hfwt2ds.cpp

Go to the documentation of this file.

// hfwt2ds.cpp
//////////////////////////////////////////////////////////////////////
//          (c) Copyright 1994-2002 Brown Deer Technology, LLC.
//                        All rights reserved.
//////////////////////////////////////////////////////////////////////

#ifdef STORM_USEF
#include "hfwt2ds_.h"
#include "dhfwt2ds_.h"
#include "ihfwt2ds_.h"
#include "hfwt2dsws_.h"
#include "dhfwt2dsws_.h"
#include "ihfwt2dsws_.h"
#endif

#include "stormdef.h"
#define INVSQRT2 STORMDEF_MATHCONST_INVSQRT2

int hfwt2ds( 
          const int& Sizef2d,
          const int& ns,
          const int& size,
          const float* sdata,
                float* wdata
) { 

#ifdef STORM_USEF

   int iflag;
   hfwt2ds_( &Sizef2d, &ns, &size, sdata, wdata, &iflag );
   return iflag;

#else

   int ii,l,size2;
   int size1 = size;
   float *u = new float[Sizef2d];
   float *v = new float[Sizef2d];

#ifdef STORM_FWT_CCHECKUSAGE
   if ( size1 < ns ) return 1;
#endif

   for( int i = 0; i < size1; i++ ) {
      ii = i * size;
      for( int j = 0; j < size1; j++ ) {
         wdata[ii+j] = sdata[ii+j];
      }
   }

   if ( size1 == ns ) return 0;

   while ( size1 > ns ) {

      size2 = size1 / 2;

      for( int j = 0; j < size1; j++ ) {        // index 1
         for( int i = 0; i < size1; i++ ) {
            ii = i * size;
            u[i] = wdata[ii+j];
         }
         for( int i = 0; i < size2; i++ ) {
            l = 2 * i;
            v[i] = ( u[l] + u[l+1] ) * INVSQRT2;
            v[size2+i] = ( u[l] - u[l+1] ) * INVSQRT2;
         }
         for( int i = 0; i < size1; i++ ) {
            ii = i * size;
            wdata[ii+j] = v[i];
         }
      }

      for( int i = 0; i < size1; i++ ) {        // index 2
         ii = i * size;
         for( int j = 0; j < size1; j++ ) {
            u[j] = wdata[ii+j];
         }
         for( int j = 0; j < size2; j++ ) {
            l = 2 * j;
            v[j] = ( u[l] + u[l+1] ) * INVSQRT2;
            v[size2+j] = ( u[l] - u[l+1] ) * INVSQRT2;
         }
         for( int j = 0; j < size1; j++ ) {
            wdata[ii+j] = v[j];
         }
      }

      size1 = size2;

   }

   delete [] u;
   delete [] v;

#ifdef STORM_FWT_CCHECKUSAGE
   if (size1 != ns) return 2;
#endif

   return 0;

#endif
} 


int hfwt2ds( 
          const int& Sizef2d,
          const int& ns,
          const int& size,
          const double* sdata,
                double* wdata
) { 

#ifdef STORM_USEF

   int iflag;
   dhfwt2ds_( &Sizef2d, &ns, &size, sdata, wdata, &iflag );
   return iflag;

#else

   int ii,jj,l,size2;
   int size1 = size;
   double *u = new double[Sizef2d];
   double *v = new double[Sizef2d];

#ifdef STORM_FWT_CCHECKUSAGE
   if ( size1 < ns ) return 1; 
#endif

   for( int i = 0; i < size1; i++ ) {
      ii = i * size;
      for( int j = 0; j < size1; j++ ) {
         wdata[ii+j] = sdata[ii+j];
      }
   }

   if ( size1 == ns ) return 0;

   while ( size1 > ns ) {

      size2 = size1 / 2;

      for( int j = 0; j < size1; j++ ) {        // index 1
         for( int i = 0; i < size1; i++ ) {
            ii = i * size;
            u[i] = wdata[ii+j];
         }
         for( int i = 0; i < size2; i++ ) {
            l = 2 * i;
            v[i] = ( u[l] + u[l+1] ) * INVSQRT2;
            v[size2+i] = ( u[l] - u[l+1] ) * INVSQRT2;
         }
         for( int i = 0; i < size1; i++ ) {
            ii = i * size;
            wdata[ii+j] = v[i];
         }
      }

      for( int i = 0; i < size1; i++ ) {        // index 2
         ii = i * size;
         for( int j = 0; j < size1; j++ ) {
            u[j] = wdata[ii+j];
         }
         for( int j = 0; j < size2; j++ ) {
            l = 2 * j;
            v[j] = ( u[l] + u[l+1] ) * INVSQRT2;
            v[size2+j] = ( u[l] - u[l+1] ) * INVSQRT2;
         }
         for( int j = 0; j < size1; j++ ) {
            jj = j * size;
            wdata[ii+j] = v[j];
         }
      }

      size1 = size2;

   }

   delete [] u;
   delete [] v;

#ifdef STORM_FWT_CCHECKUSAGE
   if (size1 != ns) return 2;
#endif

   return 0;

#endif
} 


int hfwt2ds( 
          const int& Sizef2d,
          const int& ns,
          const int& size,
          const int* sdata,
                int* wdata
) { 

#ifdef STORM_USEF

   int iflag;
   ihfwt2ds_( &Sizef2d, &ns, &size, sdata, wdata, &iflag );
   return iflag;

#else

   int ii,l,size2;
   int size1 = size;
   int *u = new int[Sizef2d];
   int *v = new int[Sizef2d];

#ifdef STORM_FWT_CCHECKUSAGE
   if ( size1 < ns ) return 1;
#endif

   for( int i = 0; i < size1; i++ ) {
      ii = i * size;
      for( int j = 0; j < size1; j++ ) {
         wdata[ii+j] = sdata[ii+j];
      }
   }

   if ( size1 == ns ) return 0;

   while ( size1 > ns ) {

      size2 = size1 / 2;

      for( int j = 0; j < size1; j++ ) {        // index 1
         for( int i = 0; i < size1; i++ ) {
            ii = i * size;
            u[i] = wdata[ii+j];
         }
         for( int i = 0; i < size2; i++ ) {
            l = 2 * i;
            v[i] = ( u[l] + u[l+1] ) / 2;
            v[size2+i] = u[l] - u[l+1];
         }
         for( int i = 0; i < size1; i++ ) {
            ii = i * size;
            wdata[ii+j] = v[i];
         }
      }

      for( int i = 0; i < size1; i++ ) {        // index 2
         ii = i * size;
         for( int j = 0; j < size1; j++ ) {
            u[j] = wdata[ii+j];
         }
         for( int j = 0; j < size2; j++ ) {
            l = 2 * j;
            v[j] = ( u[l] + u[l+1] ) / 2;
            v[size2+j] = u[l] - u[l+1];
         }
         for( int j = 0; j < size1; j++ ) {
            wdata[ii+j] = v[j];
         }
      }

      size1 = size2;

   }

   delete [] u;
   delete [] v;

#ifdef STORM_FWT_CCHECKUSAGE
   if (size1 != ns) return 2;
#endif

   return 0;

#endif
} 


int hfwt2ds( 
          const int& Sizef2d,
          const int& ns,
          const int& size,
          const float* sdata,
                float* wdata,
                float* ws
) { 

#ifdef STORM_USEF

   int iflag;
   hfwt2dsws_( &Sizef2d, &ns, &size, sdata, wdata, ws, &iflag );
   return iflag;

#else

   int ii,l,size2;
   int size1 = size;
   float *u = ws;
   float *v = &ws[Sizef2d];

#ifdef STORM_FWT_CCHECKUSAGE
   if ( size1 < ns ) return 1;
#endif

   for( int i = 0; i < size1; i++ ) {
      ii = i * size;
      for( int j = 0; j < size1; j++ ) {
         wdata[ii+j] = sdata[ii+j];
      }
   }

   if ( size1 == ns ) return 0;

   while ( size1 > ns ) {

      size2 = size1 / 2;

      for( int j = 0; j < size1; j++ ) {        // index 1
         for( int i = 0; i < size1; i++ ) {
            ii = i * size;
            u[i] = wdata[ii+j];
         }
         for( int i = 0; i < size2; i++ ) {
               l = 2 * i;
               v[i] = ( u[l] + u[l+1] ) * INVSQRT2;
               v[size2+i] = ( u[l] - u[l+1] ) * INVSQRT2;
         }
         for( int i = 0; i < size1; i++ ) {
            ii = i * size;
            wdata[ii+j] = v[i];
         }  
      }

      for( int i = 0; i < size1; i++ ) {        // index 2
         ii = i * size;
         for( int j = 0; j < size1; j++ ) {
            u[j] = wdata[ii+j];
         }
         for( int j = 0; j < size2; j++ ) {
            l = 2 * j;
            v[j] = ( u[l] + u[l+1] ) * INVSQRT2;
            v[size2+j] = ( u[l] - u[l+1] ) * INVSQRT2;
         }
         for( int j = 0; j < size1; j++ ) {
            wdata[ii+j] = v[j];
         }
      }

      size1 = size2;

   }

#ifdef STORM_FWT_CCHECKUSAGE
   if (size1 != ns) return 2;
#endif

   return 0;

#endif
} 


int hfwt2ds( 
          const int& Sizef2d,
          const int& ns,
          const int& size,
          const double* sdata,
                double* wdata,
                double* ws
) { 

#ifdef STORM_USEF

   int iflag;
   dhfwt2dsws_( &Sizef2d, &ns, &size, sdata, wdata, ws, &iflag );
   return iflag;

#else

   int ii,l,size2;
   int size1 = size;
   double* u = ws;
   double* v = &ws[Sizef2d];

#ifdef STORM_FWT_CCHECKUSAGE
   if ( size1 < ns ) return 1;
#endif

   for( int i = 0; i < size1; i++ ) {
      ii = i * size;
      for( int j = 0; j < size1; j++ ) {
         wdata[ii+j] = sdata[ii+j];
      }
   }

   if ( size1 == ns ) return 0;

   while ( size1 > ns ) {

      size2 = size1 / 2;

      for( int j = 0; j < size1; j++ ) {        // index 1
         for( int i = 0; i < size1; i++ ) {
            ii = i * size;
            u[i] = wdata[ii+j];
         }
         for( int i = 0; i < size2; i++ ) {
            l = 2 * i;
            v[i] = ( u[l] + u[l+1] ) * INVSQRT2;
            v[size2+i] = ( u[l] - u[l+1] ) * INVSQRT2;
         }
         for( int i = 0; i < size1; i++ ) {
            ii = i * size;
            wdata[ii+j] = v[i];
         }
      }

      for( int i = 0; i < size1; i++ ) {        // index 2
         ii = i * size;
         for( int j = 0; j < size1; j++ ) {
            u[j] = wdata[ii+j];
         }
         for( int j = 0; j < size2; j++ ) {
            l = 2 * j;
            v[j] = ( u[l] + u[l+1] ) * INVSQRT2;
            v[size2+j] = ( u[l] - u[l+1] ) * INVSQRT2;
         }
         for( int j = 0; j < size1; j++ ) {
            wdata[ii+j] = v[j];
         }
      }

      size1 = size2;

   }

#ifdef STORM_FWT_CCHECKUSAGE
   if (size1 != ns) return 2;
#endif

   return 0;

#endif
} 


int hfwt2ds( 
          const int& Sizef2d,
          const int& ns,
          const int& size,
          const int* sdata,
                int* wdata,
                int* ws
) { 

#ifdef STORM_USEF

   int iflag;
   ihfwt2dsws_( &Sizef2d, &ns, &size, sdata, wdata, ws, &iflag );
   return iflag;

#else

   int ii,l,size2;
   int size1 = size;
   int *u = ws;
   int *v = &ws[Sizef2d];

#ifdef STORM_FWT_CCHECKUSAGE
   if ( size1 < ns ) return 1;
#endif

   for( int i = 0; i < size1; i++ ) {
      ii = i * size;
      for( int j = 0; j < size1; j++ ) {
         wdata[ii+j] = sdata[ii+j];
      }
   }

   if ( size1 == ns ) return 0;

   while ( size1 > ns ) {

      size2 = size1 / 2;

      for( int j = 0; j < size1; j++ ) {        // index 1
         for( int i = 0; i < size1; i++ ) {
            ii = i * size;
            u[i] = wdata[ii+j];
         }
         for( int i = 0; i < size2; i++ ) {
               l = 2 * i;
               v[i] = ( u[l] + u[l+1] ) / 2;
               v[size2+i] = u[l] - u[l+1];
         }
         for( int i = 0; i < size1; i++ ) {
            ii = i * size;
            wdata[ii+j] = v[i];
         }  
      }

      for( int i = 0; i < size1; i++ ) {        // index 2
         ii = i * size;
         for( int j = 0; j < size1; j++ ) {
            u[j] = wdata[ii+j];
         }
         for( int j = 0; j < size2; j++ ) {
            l = 2 * j;
            v[j] = ( u[l] + u[l+1] ) / 2;
            v[size2+j] = u[l] - u[l+1];
         }
         for( int j = 0; j < size1; j++ ) {
            wdata[ii+j] = v[j];
         }
      }

      size1 = size2;

   }

#ifdef STORM_FWT_CCHECKUSAGE
   if (size1 != ns) return 2;
#endif

   return 0;

#endif
} 

#undef INVSQRT2

---