pd4ifwt3ds.cpp

Go to the documentation of this file.

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

#ifdef STORM_USEF
#include "pd4ifwt3ds_.h"
#include "dpd4ifwt3ds_.h"
#include "pd4ifwt3dsws_.h"
#include "dpd4ifwt3dsws_.h"
#endif

#include "stormdef.h"
#define H0 STORMDEF_MATHCONST_D4H0
#define H1 STORMDEF_MATHCONST_D4H1
#define H2 STORMDEF_MATHCONST_D4H2
#define H3 STORMDEF_MATHCONST_D4H3

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

#ifdef STORM_USEF

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

#else

   int ii,jj;
   int i1,j1,k1;
   int 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++ ) {
               i1 = (i-1+size2)%size2;
               v[2*i] = H0*u[i] + H2*u[i1]
                 + H3*u[size2+i] + H1*u[size2+i1];
               v[2*i+1] = H1*u[i] + H3*u[i1]
                 - H2*u[size2+i] - H0*u[size2+i1];
            }
            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++ ) {
               j1 = (j-1+size2)%size2;
               v[2*j] = H0*u[j] + H2*u[j1]
                 + H3*u[size2+j] + H1*u[size2+j1];
               v[2*j+1] = H1*u[j] + H3*u[j1]
                 - H2*u[size2+j] - H0*u[size2+j1];
            }
            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++ ) {
               k1 = (k-1+size2)%size2;
               v[2*k] = H0*u[k] + H2*u[k1]
                 + H3*u[size2+k] + H1*u[size2+k1];
               v[2*k+1] = H1*u[k] + H3*u[k1]
                 - H2*u[size2+k] - H0*u[size2+k1];
            }

            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 pd4ifwt3ds( 
          const int& Sizef3d,
          const int& ns,
          const int& ks,
          const int& size,
          const double* wdata,
                double* sdata
) { 

#ifdef STORM_USEF

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

#else

   int ii,jj;
   int i1,j1,k1;
   int 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++ ) {
               i1 = (i-1+size2)%size2;
               v[2*i] = H0*u[i] + H2*u[i1]
                 + H3*u[size2+i] + H1*u[size2+i1];
               v[2*i+1] = H1*u[i] + H3*u[i1]
                 - H2*u[size2+i] - H0*u[size2+i1];
            }
            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++ ) {
               j1 = (j-1+size2)%size2;
               v[2*j] = H0*u[j] + H2*u[j1]
                 + H3*u[size2+j] + H1*u[size2+j1];
               v[2*j+1] = H1*u[j] + H3*u[j1]
                 - H2*u[size2+j] - H0*u[size2+j1];
            }
            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++ ) {
               k1 = (k-1+size2)%size2;
               v[2*k] = H0*u[k] + H2*u[k1]
                 + H3*u[size2+k] + H1*u[size2+k1];
               v[2*k+1] = H1*u[k] + H3*u[k1]
                 - H2*u[size2+k] - H0*u[size2+k1];
            }

            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 pd4ifwt3ds( 
          const int& Sizef3d,
          const int& ns,
          const int& ks,
          const int& size,
          const float* wdata,
                float* sdata,
                float* ws
) { 

#ifdef STORM_USEF

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

#else

   int ii,jj;
   int i1,j1,k1;
   int 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++ ) {
               i1 = (i-1+size2)%size2;
               v[2*i] = H0*u[i] + H2*u[i1]
                 + H3*u[size2+i] + H1*u[size2+i1];
               v[2*i+1] = H1*u[i] + H3*u[i1]
                 - H2*u[size2+i] - H0*u[size2+i1];
            }
            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++ ) {
               j1 = (j-1+size2)%size2;
               v[2*j] = H0*u[j] + H2*u[j1]
                 + H3*u[size2+j] + H1*u[size2+j1];
               v[2*j+1] = H1*u[j] + H3*u[j1]
                 - H2*u[size2+j] - H0*u[size2+j1];
            }
            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++ ) {
               k1 = (k-1+size2)%size2;
               v[2*k] = H0*u[k] + H2*u[k1]
                 + H3*u[size2+k] + H1*u[size2+k1];
               v[2*k+1] = H1*u[k] + H3*u[k1]
                 - H2*u[size2+k] - H0*u[size2+k1];
            }

            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 pd4ifwt3ds( 
          const int& Sizef3d,
          const int& ns,
          const int& ks,
          const int& size,
          const double* wdata,
                double* sdata,
                double* ws
) { 

#ifdef STORM_USEF

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

#else

   int ii,jj;
   int i1,j1,k1;
   int 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++ ) {
               i1 = (i-1+size2)%size2;
               v[2*i] = H0*u[i] + H2*u[i1]
                 + H3*u[size2+i] + H1*u[size2+i1];
               v[2*i+1] = H1*u[i] + H3*u[i1]
                 - H2*u[size2+i] - H0*u[size2+i1];
            }
            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++ ) {
               j1 = (j-1+size2)%size2;
               v[2*j] = H0*u[j] + H2*u[j1]
                 + H3*u[size2+j] + H1*u[size2+j1];
               v[2*j+1] = H1*u[j] + H3*u[j1]
                 - H2*u[size2+j] - H0*u[size2+j1];
            }
            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++ ) {
               k1 = (k-1+size2)%size2;
               v[2*k] = H0*u[k] + H2*u[k1]
                 + H3*u[size2+k] + H1*u[size2+k1];
               v[2*k+1] = H1*u[k] + H3*u[k1]
                 - H2*u[size2+k] - H0*u[size2+k1];
            }

            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 H0
#undef H1
#undef H2
#undef H3

---