SR2k4 Assembler: pd4fwt3ds.cpp Source File

00001 // pd4fwt3ds.cpp
00002 //////////////////////////////////////////////////////////////////////
00003 //          (c) Copyright 1994-2002 Brown Deer Technology, LLC.
00004 //                        All rights reserved.
00005 //////////////////////////////////////////////////////////////////////
00006 
00007 #ifdef STORM_USEF
00008 #include "pd4fwt3ds_.h"
00009 #include "dpd4fwt3ds_.h"
00010 #include "pd4fwt3dsws_.h"
00011 #include "dpd4fwt3dsws_.h"
00012 #endif
00013 
00014 #include "stormdef.h"
00015 #define H0 STORMDEF_MATHCONST_D4H0
00016 #define H1 STORMDEF_MATHCONST_D4H1
00017 #define H2 STORMDEF_MATHCONST_D4H2
00018 #define H3 STORMDEF_MATHCONST_D4H3
00019 
00020 int pd4fwt3ds( 
00021           const int& Sizef3d,
00022           const int& ns,
00023           const int& ks,
00024           const int& size,
00025           const float* sdata,
00026                 float* wdata
00027 ) { 
00028 
00029 #ifdef STORM_USEF
00030 
00031    int iflag;
00032    pd4fwt3ds_( &Sizef3d, &ns, &ks, &size, sdata, wdata, &iflag );
00033    return iflag;
00034 
00035 #else
00036 
00037    int ii,jj,l;
00038    int l1,l2,l3;
00039    int size2;
00040    int size1 = size;
00041    float *u = new float[Sizef3d];
00042    float *v = new float[Sizef3d];
00043 
00044 #ifdef STORM_FWT_CCHECKUSAGE
00045    if ( size1 < ns ) return 1;
00046 #endif
00047 
00048    for( int i = 0; i < size1; i++ ) {
00049       ii=i*size*size;
00050       for( int j = 0; j < size1; j++ ) {
00051          jj=j*size;
00052          for( int k = 0; k < size1; k++ ) {
00053             wdata[ii+jj+k] = sdata[ii+jj+k];
00054          }
00055       }
00056    }
00057 
00058    if ( size1 == ns ) return 0;
00059 
00060    while ( size1 > ns ) {
00061 
00062       size2 = size1 / 2;
00063 
00064       for( int j = 0; j < size1; j++ ) {        // index 1
00065          jj = j * size;
00066          for( int k = 0; k < size1; k++ ) {
00067             for( int i = 0; i < size1; i++ ) {
00068                ii = i * size * size;
00069                u[i] = wdata[ii+jj+k];
00070             }
00071             for( int i = 0; i < size2; i++ ) {
00072                l = 2 * i;
00073                l1 = (l+1)%size1;
00074                l2 = (l+2)%size1;
00075                l3 = (l+3)%size1;
00076                v[i] = H0*u[l] + H1*u[l1] + H2*u[l2] + H3*u[l3];
00077                v[size2+i] = H3*u[l] - H2*u[l1] + H1*u[l2] - H0*u[l3]; 
00078             }
00079             for( int i = 0; i < size1; i++ ) {
00080                ii = i * size * size;
00081                wdata[ii+jj+k] = v[i];
00082             }
00083          }  
00084       }
00085 
00086       for( int i = 0; i < size1; i++ ) {        // index 2
00087          ii = i * size * size;
00088          for( int k = 0; k < size1; k++ ) {
00089             for( int j = 0; j < size1; j++ ) {
00090                jj = j * size;
00091                u[j] = wdata[ii+jj+k];
00092             }
00093             for( int j = 0; j < size2; j++ ) {
00094                l = 2 * j;
00095                l1 = (l+1)%size1;
00096                l2 = (l+2)%size1;
00097                l3 = (l+3)%size1;
00098                v[j] = H0*u[l] + H1*u[l1] + H2*u[l2] + H3*u[l3];
00099                v[size2+j] = H3*u[l] - H2*u[l1] + H1*u[l2] - H0*u[l3]; 
00100             }
00101             for( int j = 0; j < size1; j++ ) {
00102                jj = j * size;
00103                wdata[ii+jj+k] = v[j];
00104             }
00105          }
00106       }
00107 
00108       for( int i = 0; i < size1; i++ ) {        // index 3
00109          ii = i * size * size;
00110          for( int j = 0; j < size1; j++ ) {
00111             jj = j * size;
00112 
00113             for( int k = 0; k < size1; k++ ) {
00114                u[k] = wdata[ii+jj+k];
00115             }
00116 
00117             for( int k = 0; k < size2; k++ ) {
00118                l = 2 * k;
00119                l1 = (l+1)%size1;
00120                l2 = (l+2)%size1;
00121                l3 = (l+3)%size1;
00122                v[k] = H0*u[l] + H1*u[l1] + H2*u[l2] + H3*u[l3];
00123                v[size2+k] = H3*u[l] - H2*u[l1] + H1*u[l2] - H0*u[l3]; 
00124             }
00125 
00126             for( int k = 0; k < size1; k++ ) {
00127                wdata[ii+jj+k] = v[k];
00128             }
00129          }
00130       }
00131 
00132       size1 = size2;
00133 
00134    }
00135 
00136    delete [] u;
00137    delete [] v;
00138 
00139 #ifdef STORM_FWT_CCHECKUSAGE
00140    if (size1 != ns) return 2;
00141 #endif
00142 
00143    return 0;
00144 
00145 #endif
00146 }
00147  
00148 
00149 int pd4fwt3ds( 
00150           const int& Sizef3d,
00151           const int& ns,
00152           const int& ks,
00153           const int& size,
00154           const double* sdata,
00155                 double* wdata
00156 ) { 
00157 
00158 #ifdef STORM_USEF
00159 
00160    int iflag;
00161    dpd4fwt3ds_( &Sizef3d, &ns, &ks, &size, sdata, wdata, &iflag );
00162    return iflag;
00163 
00164 #else
00165 
00166    int ii,jj,l;
00167    int l1,l2,l3;
00168    int size2;
00169    int size1 = size;
00170    double *u = new double[Sizef3d];
00171    double *v = new double[Sizef3d];
00172 
00173 #ifdef STORM_FWT_CCHECKUSAGE
00174    if ( size1 < ns ) return 1;
00175 #endif
00176 
00177    for( int i = 0; i < size1; i++ ) {
00178       ii=i*size*size;
00179       for( int j = 0; j < size1; j++ ) {
00180          jj=j*size;
00181          for( int k = 0; k < size1; k++ ) {
00182             wdata[ii+jj+k] = sdata[ii+jj+k];
00183          }
00184       }
00185    }
00186 
00187    if ( size1 == ns ) return 0;
00188 
00189    while ( size1 > ns ) {
00190 
00191       size2 = size1 / 2;
00192 
00193       for( int j = 0; j < size1; j++ ) {        // index 1
00194          jj = j * size;
00195          for( int k = 0; k < size1; k++ ) {
00196             for( int i = 0; i < size1; i++ ) {
00197                ii = i * size * size;
00198                u[i] = wdata[ii+jj+k];
00199             }
00200             for( int i = 0; i < size2; i++ ) {
00201                l = 2 * i;
00202                l1 = (l+1)%size1;
00203                l2 = (l+2)%size1;
00204                l3 = (l+3)%size1;
00205                v[i] = H0*u[l] + H1*u[l1] + H2*u[l2] + H3*u[l3];
00206                v[size2+i] = H3*u[l] - H2*u[l1] + H1*u[l2] - H0*u[l3]; 
00207             }
00208             for( int i = 0; i < size1; i++ ) {
00209                ii = i * size * size;
00210                wdata[ii+jj+k] = v[i];
00211             }
00212          }  
00213       }
00214 
00215       for( int i = 0; i < size1; i++ ) {        // index 2
00216          ii = i * size * size;
00217          for( int k = 0; k < size1; k++ ) {
00218             for( int j = 0; j < size1; j++ ) {
00219                jj = j * size;
00220                u[j] = wdata[ii+jj+k];
00221             }
00222             for( int j = 0; j < size2; j++ ) {
00223                l = 2 * j;
00224                l1 = (l+1)%size1;
00225                l2 = (l+2)%size1;
00226                l3 = (l+3)%size1;
00227                v[j] = H0*u[l] + H1*u[l1] + H2*u[l2] + H3*u[l3];
00228                v[size2+j] = H3*u[l] - H2*u[l1] + H1*u[l2] - H0*u[l3]; 
00229             }
00230             for( int j = 0; j < size1; j++ ) {
00231                jj = j * size;
00232                wdata[ii+jj+k] = v[j];
00233             }
00234          }
00235       }
00236 
00237       for( int i = 0; i < size1; i++ ) {        // index 3
00238          ii = i * size * size;
00239          for( int j = 0; j < size1; j++ ) {
00240             jj = j * size;
00241 
00242             for( int k = 0; k < size1; k++ ) {
00243                u[k] = wdata[ii+jj+k];
00244             }
00245 
00246             for( int k = 0; k < size2; k++ ) {
00247                l = 2 * k;
00248                l1 = (l+1)%size1;
00249                l2 = (l+2)%size1;
00250                l3 = (l+3)%size1;
00251                v[k] = H0*u[l] + H1*u[l1] + H2*u[l2] + H3*u[l3];
00252                v[size2+k] = H3*u[l] - H2*u[l1] + H1*u[l2] - H0*u[l3]; 
00253             }
00254 
00255             for( int k = 0; k < size1; k++ ) {
00256                wdata[ii+jj+k] = v[k];
00257             }
00258          }
00259       }
00260 
00261       size1 = size2;
00262 
00263    }
00264 
00265    delete [] u;
00266    delete [] v;
00267 
00268 #ifdef STORM_FWT_CCHECKUSAGE
00269    if (size1 != ns) return 2;
00270 #endif
00271 
00272    return 0;
00273 
00274 #endif
00275 } 
00276 
00277 
00278 
00279 int pd4fwt3ds( 
00280           const int& Sizef3d,
00281           const int& ns,
00282           const int& ks,
00283           const int& size,
00284           const float* sdata,
00285                 float* wdata,
00286                 float* ws
00287 ) { 
00288 
00289 #ifdef STORM_USEF
00290 
00291    int iflag;
00292    pd4fwt3dsws_( &Sizef3d, &ns, &ks, &size, sdata, wdata, ws, &iflag );
00293    return iflag;
00294 
00295 #else
00296 
00297    int ii,jj,l;
00298    int l1,l2,l3;
00299    int size2;
00300    int size1 = size;
00301    float *u = ws;
00302    float *v = &ws[Sizef3d];
00303 
00304 #ifdef STORM_FWT_CCHECKUSAGE
00305    if ( size1 < ns ) return 1;
00306 #endif
00307 
00308    for( int i = 0; i < size1; i++ ) {
00309       ii=i*size*size;
00310       for( int j = 0; j < size1; j++ ) {
00311          jj=j*size;
00312          for( int k = 0; k < size1; k++ ) {
00313             wdata[ii+jj+k] = sdata[ii+jj+k];
00314          }
00315       }
00316    }
00317 
00318    if ( size1 == ns ) return 0;
00319 
00320    while ( size1 > ns ) {
00321 
00322       size2 = size1 / 2;
00323 
00324       for( int j = 0; j < size1; j++ ) {        // index 1
00325          jj = j * size;
00326          for( int k = 0; k < size1; k++ ) {
00327             for( int i = 0; i < size1; i++ ) {
00328                ii = i * size * size;
00329                u[i] = wdata[ii+jj+k];
00330             }
00331             for( int i = 0; i < size2; i++ ) {
00332                l = 2 * i;
00333                l1 = (l+1)%size1;
00334                l2 = (l+2)%size1;
00335                l3 = (l+3)%size1;
00336                v[i] = H0*u[l] + H1*u[l1] + H2*u[l2] + H3*u[l3];
00337                v[size2+i] = H3*u[l] - H2*u[l1] + H1*u[l2] - H0*u[l3]; 
00338             }
00339             for( int i = 0; i < size1; i++ ) {
00340                ii = i * size * size;
00341                wdata[ii+jj+k] = v[i];
00342             }
00343          }  
00344       }
00345 
00346       for( int i = 0; i < size1; i++ ) {        // index 2
00347          ii = i * size * size;
00348          for( int k = 0; k < size1; k++ ) {
00349             for( int j = 0; j < size1; j++ ) {
00350                jj = j * size;
00351                u[j] = wdata[ii+jj+k];
00352             }
00353             for( int j = 0; j < size2; j++ ) {
00354                l = 2 * j;
00355                l1 = (l+1)%size1;
00356                l2 = (l+2)%size1;
00357                l3 = (l+3)%size1;
00358                v[j] = H0*u[l] + H1*u[l1] + H2*u[l2] + H3*u[l3];
00359                v[size2+j] = H3*u[l] - H2*u[l1] + H1*u[l2] - H0*u[l3]; 
00360             }
00361             for( int j = 0; j < size1; j++ ) {
00362                jj = j * size;
00363                wdata[ii+jj+k] = v[j];
00364             }
00365          }
00366       }
00367 
00368       for( int i = 0; i < size1; i++ ) {        // index 3
00369          ii = i * size * size;
00370          for( int j = 0; j < size1; j++ ) {
00371             jj = j * size;
00372 
00373             for( int k = 0; k < size1; k++ ) {
00374                u[k] = wdata[ii+jj+k];
00375             }
00376 
00377             for( int k = 0; k < size2; k++ ) {
00378                l = 2 * k;
00379                l1 = (l+1)%size1;
00380                l2 = (l+2)%size1;
00381                l3 = (l+3)%size1;
00382                v[k] = H0*u[l] + H1*u[l1] + H2*u[l2] + H3*u[l3];
00383                v[size2+k] = H3*u[l] - H2*u[l1] + H1*u[l2] - H0*u[l3]; 
00384             }
00385 
00386             for( int k = 0; k < size1; k++ ) {
00387                wdata[ii+jj+k] = v[k];
00388             }
00389          }
00390       }
00391 
00392       size1 = size2;
00393 
00394    }
00395 
00396 #ifdef STORM_FWT_CCHECKUSAGE
00397    if (size1 != ns) return 2;
00398 #endif
00399 
00400    return 0;
00401 
00402 #endif
00403 } 
00404 
00405 
00406 int pd4fwt3ds( 
00407           const int& Sizef3d,
00408           const int& ns,
00409           const int& ks,
00410           const int& size,
00411           const double* sdata,
00412                 double* wdata,
00413                 double* ws
00414 ) { 
00415 
00416 #ifdef STORM_USEF
00417 
00418    int iflag;
00419    dpd4fwt3dsws_( &Sizef3d, &ns, &ks, &size, sdata, wdata, ws, &iflag );
00420    return iflag;
00421 
00422 #else
00423 
00424    int ii,jj,l;
00425    int l1,l2,l3;
00426    int size2;
00427    int size1 = size;
00428    double* u = ws;
00429    double* v = &ws[Sizef3d];
00430 
00431 #ifdef STORM_FWT_CCHECKUSAGE
00432    if ( size1 < ns ) return 1;
00433 #endif
00434 
00435    for( int i = 0; i < size1; i++ ) {
00436       ii=i*size*size;
00437       for( int j = 0; j < size1; j++ ) {
00438          jj=j*size;
00439          for( int k = 0; k < size1; k++ ) {
00440             wdata[ii+jj+k] = sdata[ii+jj+k];
00441          }
00442       }
00443    }
00444 
00445    if ( size1 == ns ) return 0;
00446 
00447    while ( size1 > ns ) {
00448 
00449       size2 = size1 / 2;
00450 
00451       for( int j = 0; j < size1; j++ ) {        // index 1
00452          jj = j * size;
00453          for( int k = 0; k < size1; k++ ) {
00454             for( int i = 0; i < size1; i++ ) {
00455                ii = i * size * size;
00456                u[i] = wdata[ii+jj+k];
00457             }
00458             for( int i = 0; i < size2; i++ ) {
00459                l = 2 * i;
00460                l1 = (l+1)%size1;
00461                l2 = (l+2)%size1;
00462                l3 = (l+3)%size1;
00463                v[i] = H0*u[l] + H1*u[l1] + H2*u[l2] + H3*u[l3];
00464                v[size2+i] = H3*u[l] - H2*u[l1] + H1*u[l2] - H0*u[l3]; 
00465             }
00466             for( int i = 0; i < size1; i++ ) {
00467                ii = i * size * size;
00468                wdata[ii+jj+k] = v[i];
00469             }
00470          }  
00471       }
00472 
00473       for( int i = 0; i < size1; i++ ) {        // index 2
00474          ii = i * size * size;
00475          for( int k = 0; k < size1; k++ ) {
00476             for( int j = 0; j < size1; j++ ) {
00477                jj = j * size;
00478                u[j] = wdata[ii+jj+k];
00479             }
00480             for( int j = 0; j < size2; j++ ) {
00481                l = 2 * j;
00482                l1 = (l+1)%size1;
00483                l2 = (l+2)%size1;
00484                l3 = (l+3)%size1;
00485                v[j] = H0*u[l] + H1*u[l1] + H2*u[l2] + H3*u[l3];
00486                v[size2+j] = H3*u[l] - H2*u[l1] + H1*u[l2] - H0*u[l3]; 
00487             }
00488             for( int j = 0; j < size1; j++ ) {
00489                jj = j * size;
00490                wdata[ii+jj+k] = v[j];
00491             }
00492          }
00493       }
00494 
00495       for( int i = 0; i < size1; i++ ) {        // index 3
00496          ii = i * size * size;
00497          for( int j = 0; j < size1; j++ ) {
00498             jj = j * size;
00499 
00500             for( int k = 0; k < size1; k++ ) {
00501                u[k] = wdata[ii+jj+k];
00502             }
00503 
00504             for( int k = 0; k < size2; k++ ) {
00505                l = 2 * k;
00506                l1 = (l+1)%size1;
00507                l2 = (l+2)%size1;
00508                l3 = (l+3)%size1;
00509                v[k] = H0*u[l] + H1*u[l1] + H2*u[l2] + H3*u[l3];
00510                v[size2+k] = H3*u[l] - H2*u[l1] + H1*u[l2] - H0*u[l3]; 
00511             }
00512 
00513             for( int k = 0; k < size1; k++ ) {
00514                wdata[ii+jj+k] = v[k];
00515             }
00516          }
00517       }
00518 
00519       size1 = size2;
00520 
00521    }
00522 
00523 #ifdef STORM_FWT_CCHECKUSAGE
00524    if (size1 != ns) return 2;
00525 #endif
00526 
00527    return 0;
00528 
00529 #endif
00530 } 
00531 
00532 #undef H0
00533 #undef H1
00534 #undef H2
00535 #undef H3
00536