hfwt3ds.cpp

Go to the documentation of this file.

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

#ifdef STORM_USEF
#include "hfwt3ds_.h"
#include "dhfwt3ds_.h"
#include "ihfwt3ds_.h"
#include "hfwt3dsws_.h"
#include "dhfwt3dsws_.h"
#include "ihfwt3dsws_.h"
#endif

#include "stormdef.h"
#define INVSQRT2 STORMDEF_MATHCONST_INVSQRT2

int hfwt3ds( 
          const int& Sizef3d,
          const int& ns,
          const int& size,
          const float* sdata,
                float* wdata
) { 

#ifdef STORM_USEF

   int iflag;
   hfwt3ds_( &Sizef3d, &ns, &size, sdata, wdata, &iflag );
   return iflag;

#else

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

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

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

   if ( size1 == ns ) return 0;

   while ( size1 > ns ) {

      size2 = size1 / 2;

      for( int j = 0; j < size1; j++ ) {        // index 1
         jj = j * size;
         for( int k = 0; k < size1; k++ ) {
            for( int i = 0; i < size1; i++ ) {
               ii = i * size * size;
               u[i] = wdata[ii+jj+k];
            }
            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 * size;
               wdata[ii+jj+k] = v[i];
            }
         }  
      }

      for( int i = 0; i < size1; i++ ) {        // index 2
         ii = i * size * size;
         for( int k = 0; k < size1; k++ ) {
            for( int j = 0; j < size1; j++ ) {
               jj = j * size;
               u[j] = wdata[ii+jj+k];
            }
            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+jj+k] = v[j];
            }
         }
      }

      for( int i = 0; i < size1; i++ ) {        // index 3
         ii = i * size * size;
         for( int j = 0; j < size1; j++ ) {
            jj = j * size;

            for( int k = 0; k < size1; k++ ) {
               u[k] = wdata[ii+jj+k];
            }

            for( int k = 0; k < size2; k++ ) {
               l = 2 * k;
               v[k] = ( u[l] + u[l+1] ) * INVSQRT2;
               v[size2+k] = ( u[j] - u[l+1] ) * INVSQRT2;
            }

            for( int k = 0; k < size1; k++ ) {
               wdata[ii+jj+k] = v[k];
            }
         }
      }

      size1 = size2;

   }

   delete [] u;
   delete [] v;

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

   return 0;

#endif
}
 

int hfwt3ds( 
          const int& Sizef3d,
          const int& ns,
          const int& size,
          const double* sdata,
                double* wdata
) { 

#ifdef STORM_USEF

   int iflag;
   dhfwt3ds_( &Sizef3d, &ns, &size, sdata, wdata, &iflag );
   return iflag;

#else

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

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

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

   if ( size1 == ns ) return 0;

   while ( size1 > ns ) {

      size2 = size1 / 2;

      for( int j = 0; j < size1; j++ ) {        // index 1
         jj = j * size;
         for( int k = 0; k < size1; k++ ) {
            for( int i = 0; i < size1; i++ ) {
               ii = i * size * size;
               u[i] = wdata[ii+jj+k];
            }
            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 * size;
               wdata[ii+jj+k] = v[i];
            }
         }  
      }

      for( int i = 0; i < size1; i++ ) {        // index 2
         ii = i * size * size;
         for( int k = 0; k < size1; k++ ) {
            for( int j = 0; j < size1; j++ ) {
               jj = j * size;
               u[j] = wdata[ii+jj+k];
            }
            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+jj+k] = v[j];
            }
         }
      }

      for( int i = 0; i < size1; i++ ) {        // index 3
         ii = i * size * size;
         for( int j = 0; j < size1; j++ ) {
            jj = j * size;

            for( int k = 0; k < size1; k++ ) {
               u[k] = wdata[ii+jj+k];
            }

            for( int k = 0; k < size2; k++ ) {
               l = 2 * k;
               v[k] = ( u[l] + u[l+1] ) * INVSQRT2;
               v[size2+k] = ( u[j] - u[l+1] ) * INVSQRT2;
            }

            for( int k = 0; k < size1; k++ ) {
               wdata[ii+jj+k] = v[k];
            }
         }
      }

      size1 = size2;

   }

   delete [] u;
   delete [] v;

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

   return 0;

#endif
} 


int hfwt3ds( 
          const int& Sizef3d,
          const int& ns,
          const int& size,
          const int* sdata,
                int* wdata
) { 

#ifdef STORM_USEF

   int iflag;
   ihfwt3ds_( &Sizef3d, &ns, &size, sdata, wdata, &iflag );
   return iflag;

#else

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

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

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

   if ( size1 == ns ) return 0;

   while ( size1 > ns ) {

      size2 = size1 / 2;

      for( int j = 0; j < size1; j++ ) {        // index 1
         jj = j * size;
         for( int k = 0; k < size1; k++ ) {
            for( int i = 0; i < size1; i++ ) {
               ii = i * size * size;
               u[i] = wdata[ii+jj+k];
            }
            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 * size;
               wdata[ii+jj+k] = v[i];
            }
         }  
      }

      for( int i = 0; i < size1; i++ ) {        // index 2
         ii = i * size * size;
         for( int k = 0; k < size1; k++ ) {
            for( int j = 0; j < size1; j++ ) {
               jj = j * size;
               u[j] = wdata[ii+jj+k];
            }
            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++ ) {
               jj = j * size;
               wdata[ii+jj+k] = v[j];
            }
         }
      }

      for( int i = 0; i < size1; i++ ) {        // index 3
         ii = i * size * size;
         for( int j = 0; j < size1; j++ ) {
            jj = j * size;

            for( int k = 0; k < size1; k++ ) {
               u[k] = wdata[ii+jj+k];
            }

            for( int k = 0; k < size2; k++ ) {
               l = 2 * k;
               v[k] = ( u[l] + u[l+1] ) / 2;
               v[size2+k] = u[j] - u[l+1];
            }

            for( int k = 0; k < size1; k++ ) {
               wdata[ii+jj+k] = v[k];
            }
         }
      }

      size1 = size2;

   }

   delete [] u;
   delete [] v;

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

   return 0;

#endif
} 


int hfwt3ds( 
          const int& Sizef3d,
          const int& ns,
          const int& size,
          const float* sdata,
                float* wdata,
                float* ws
) { 

#ifdef STORM_USEF

   int iflag;
   hfwt3dsws_( &Sizef3d, &ns, &size, sdata, wdata, ws, &iflag );
   return iflag;

#else

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

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

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

   if ( size1 == ns ) return 0;

   while ( size1 > ns ) {

      size2 = size1 / 2;

      for( int j = 0; j < size1; j++ ) {        // index 1
         jj = j * size;
         for( int k = 0; k < size1; k++ ) {
            for( int i = 0; i < size1; i++ ) {
               ii = i * size * size;
               u[i] = wdata[ii+jj+k];
            }
            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 * size;
               wdata[ii+jj+k] = v[i];
            }
         }  
      }

      for( int i = 0; i < size1; i++ ) {        // index 2
         ii = i * size * size;
         for( int k = 0; k < size1; k++ ) {
            for( int j = 0; j < size1; j++ ) {
               jj = j * size;
               u[j] = wdata[ii+jj+k];
            }
            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+jj+k] = v[j];
            }
         }
      }

      for( int i = 0; i < size1; i++ ) {        // index 3
         ii = i * size * size;
         for( int j = 0; j < size1; j++ ) {
            jj = j * size;

            for( int k = 0; k < size1; k++ ) {
               u[k] = wdata[ii+jj+k];
            }

            for( int k = 0; k < size2; k++ ) {
               l = 2 * k;
               v[k] = ( u[l] + u[l+1] ) * INVSQRT2;
               v[size2+k] = ( u[j] - u[l+1] ) * INVSQRT2;
            }

            for( int k = 0; k < size1; k++ ) {
               wdata[ii+jj+k] = v[k];
            }
         }
      }

      size1 = size2;

   }

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

   return 0;

#endif
} 


int hfwt3ds( 
          const int& Sizef3d,
          const int& ns,
          const int& size,
          const double* sdata,
                double* wdata,
                double* ws
) { 

#ifdef STORM_USEF

   int iflag;
   dhfwt3dsws_( &Sizef3d, &ns, &size, sdata, wdata, ws, &iflag );
   return iflag;

#else

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

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

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

   if ( size1 == ns ) return 0;

   while ( size1 > ns ) {

      size2 = size1 / 2;

      for( int j = 0; j < size1; j++ ) {        // index 1
         jj = j * size;
         for( int k = 0; k < size1; k++ ) {
            for( int i = 0; i < size1; i++ ) {
               ii = i * size * size;
               u[i] = wdata[ii+jj+k];
            }
            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 * size;
               wdata[ii+jj+k] = v[i];
            }
         }  
      }

      for( int i = 0; i < size1; i++ ) {        // index 2
         ii = i * size * size;
         for( int k = 0; k < size1; k++ ) {
            for( int j = 0; j < size1; j++ ) {
               jj = j * size;
               u[j] = wdata[ii+jj+k];
            }
            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+jj+k] = v[j];
            }
         }
      }

      for( int i = 0; i < size1; i++ ) {        // index 3
         ii = i * size * size;
         for( int j = 0; j < size1; j++ ) {
            jj = j * size;

            for( int k = 0; k < size1; k++ ) {
               u[k] = wdata[ii+jj+k];
            }

            for( int k = 0; k < size2; k++ ) {
               l = 2 * k;
               v[k] = ( u[l] + u[l+1] ) * INVSQRT2;
               v[size2+k] = ( u[j] - u[l+1] ) * INVSQRT2;
            }

            for( int k = 0; k < size1; k++ ) {
               wdata[ii+jj+k] = v[k];
            }
         }
      }

      size1 = size2;

   }

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

   return 0;

#endif
} 


int hfwt3ds( 
          const int& Sizef3d,
          const int& ns,
          const int& size,
          const int* sdata,
                int* wdata,
                int* ws
) { 

#ifdef STORM_USEF

   int iflag;
   ihfwt3dsws_( &Sizef3d, &ns, &size, sdata, wdata, ws, &iflag );
   return iflag;

#else

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

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

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

   if ( size1 == ns ) return 0;

   while ( size1 > ns ) {

      size2 = size1 / 2;

      for( int j = 0; j < size1; j++ ) {        // index 1
         jj = j * size;
         for( int k = 0; k < size1; k++ ) {
            for( int i = 0; i < size1; i++ ) {
               ii = i * size * size;
               u[i] = wdata[ii+jj+k];
            }
            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 * size;
               wdata[ii+jj+k] = v[i];
            }
         }  
      }

      for( int i = 0; i < size1; i++ ) {        // index 2
         ii = i * size * size;
         for( int k = 0; k < size1; k++ ) {
            for( int j = 0; j < size1; j++ ) {
               jj = j * size;
               u[j] = wdata[ii+jj+k];
            }
            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++ ) {
               jj = j * size;
               wdata[ii+jj+k] = v[j];
            }
         }
      }

      for( int i = 0; i < size1; i++ ) {        // index 3
         ii = i * size * size;
         for( int j = 0; j < size1; j++ ) {
            jj = j * size;

            for( int k = 0; k < size1; k++ ) {
               u[k] = wdata[ii+jj+k];
            }

            for( int k = 0; k < size2; k++ ) {
               l = 2 * k;
               v[k] = ( u[l] + u[l+1] ) / 2;
               v[size2+k] = u[j] - u[l+1];
            }

            for( int k = 0; k < size1; k++ ) {
               wdata[ii+jj+k] = v[k];
            }
         }
      }

      size1 = size2;

   }

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

   return 0;

#endif
} 

#undef INVSQRT2

---