hifwt3ds.cpp

Go to the documentation of this file.

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

#ifdef STORM_USEF
#include "hifwt3ds_.h"
#include "dhifwt3ds_.h"
#include "ihifwt3ds_.h"
#include "hifwt3dsws_.h"
#include "dhifwt3dsws_.h"
#include "ihifwt3dsws_.h"
#endif

#include "stormdef.h"
#define INVSQRT2 STORMDEF_MATHCONST_INVSQRT2

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

#ifdef STORM_USEF

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

#else

   int ii,jj,size2;
   int size1 = ns;
   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 < size; i++ ) {
      ii=i*size*size;
      for( int j = 0; j < size; j++ ) {
         jj=j*size;
         for( int k = 0; k < size; k++ ) {
            sdata[ii+jj+k] = wdata[ii+jj+k];
         }
      }
   }

   if ( size == ns ) return 0;

   while ( size1 < size ) {

      size2 = size1;
      size1 = 2 * size1;

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

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

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

   }
   
   delete [] u;
   delete [] v;

#ifdef STORM_FWT_CCHECKUSAGE
   if ( size1 > size ) return 2;
#endif

   return 0;

#endif
} 


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

#ifdef STORM_USEF

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

#else

   int ii,jj,size2;
   int size1 = ns;
   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 < size; i++ ) {
      ii=i*size*size;
      for( int j = 0; j < size; j++ ) {
         jj=j*size;
         for( int k = 0; k < size; k++ ) {
            sdata[ii+jj+k] = wdata[ii+jj+k];
         }
      }
   }

   if ( size == ns ) return 0;

   while ( size1 < size ) {

      size2 = size1;
      size1 = 2 * size1;

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

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

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

   }

   delete [] u;
   delete [] v;


#ifdef STORM_FWT_CCHECKUSAGE
   if ( size1 > size ) return 2;
#endif

   return 0;

#endif
} 


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

#ifdef STORM_USEF

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

#else

   int ii,jj,size2;
   int size1 = ns;
   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 < size; i++ ) {
      ii=i*size*size;
      for( int j = 0; j < size; j++ ) {
         jj=j*size;
         for( int k = 0; k < size; k++ ) {
            sdata[ii+jj+k] = wdata[ii+jj+k];
         }
      }
   }

   if ( size == ns ) return 0;

   while ( size1 < size ) {

      size2 = size1;
      size1 = 2 * size1;

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

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

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

   }
   
   delete [] u;
   delete [] v;

#ifdef STORM_FWT_CCHECKUSAGE
   if ( size1 > size ) return 2;
#endif

   return 0;

#endif
} 


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

#ifdef STORM_USEF

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

#else

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

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

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

   if ( size == ns ) return 0;

   while ( size1 < size ) {

      size2 = size1;
      size1 = 2 * size1;

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

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

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

   }
   
#ifdef STORM_FWT_CCHECKUSAGE
   if ( size1 > size ) return 2;
#endif

   return 0;

#endif
} 


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

#ifdef STORM_USEF

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

#else

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

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

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

   if ( size == ns ) return 0;

   while ( size1 < size ) {

      size2 = size1;
      size1 = 2 * size1;

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

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

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

   }

#ifdef STORM_FWT_CCHECKUSAGE
   if ( size1 > size ) return 2;
#endif

   return 0;

#endif
} 


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

#ifdef STORM_USEF

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

#else

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

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

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

   if ( size == ns ) return 0;

   while ( size1 < size ) {

      size2 = size1;
      size1 = 2 * size1;

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

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

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

   }
   
#ifdef STORM_FWT_CCHECKUSAGE
   if ( size1 > size ) return 2;
#endif

   return 0;

#endif
} 

#undef INVSQRT2

---