SR2k4 Assembler: g1fwt3ds.cpp Source File

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