SR2k4 Assembler: g1ifwt3ds.cpp Source File

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