AdaptiveLevels.cc

00001 # include "AdaptiveLevels.hpp"
00002 
00003 extern int debugRefine ;
00004 
00006 //                  //
00007 //  AdaptiveLevels  //
00008 //                  //
00010 
00011 // functions for BlockAllocator
00012 void AdaptiveLevels::Init() {
00013   for (int k=0; k<=LMAX; k++) a[k] = new AdaptiveLevels::AdaptiveL() ;
00014 }
00015 
00016 void AdaptiveLevels::Clear() {
00017   for (int k=0; k<=LMAX; k++) (*a[k]).clear();
00018 }
00019 
00020 void AdaptiveLevels::Delete() {
00021   int k ;
00022   for (k=0; k<=LMAX; k++) { 
00023     delete a[k]; 
00024     a[k]=NULL ; 
00025   }
00026 }
00027 
00028 // other member functions
00029 size_t AdaptiveLevels::size() {
00030   size_t i,s=0 ;
00031   for (i=0; i<=LMAX; i++) s += (*a[i]).size() ;
00032   return s ;
00033 }
00034 
00035 // generate AdaptiveLevels:AdaptiveL with entries from *I
00036 size_t IndexSetWriteToAdaptiveL(AdaptiveLevels::AdaptiveL *al,size_t *counter,IndexSet *I) {
00037   int    *a=I->a,*e=I->e,nr=I->nr ;
00038   int    r,i=0 ;
00039   size_t c=(*counter),c0=(*counter) ;
00040   
00041   (*al).clear() ;
00042   for (r=0;r<nr;r++)
00043     for (i=a[r];i<=e[r]; i++) {
00044       pair<const int,long> p(i,(long)c) ;
00045       (*al).insert((*al).end(),p) ;
00046       c++ ;
00047      }
00048    
00049   *counter = c ;
00050   return c-c0 ;
00051 }
00052 
00054 //                       //
00055 //  AdaptiveLevelsIndex  //
00056 //                       //
00058 
00059 void AdaptiveLevelsIndex::Print() {
00060   int i;
00061   std::cout <<dim<<" (" ;
00062   if (dim>0) { 
00063     std::cout <<l[0] ;
00064     for (i=1; i<dim; i++) std::cout <<','<<l[i] ;
00065   }
00066   std::cout <<"),(" ;
00067   if (dim>0) { 
00068     std::cout <<t[0] ;
00069     for (i=1; i<dim; i++) std::cout <<','<<t[i] ;
00070   }
00071   std::cout << ")\n" ;
00072 }
00073 
00074 // comparison
00075 bool AdaptiveLevelsIndexCmp::operator()(const AdaptiveLevelsIndex *a,const AdaptiveLevelsIndex *b) const {
00076   int i ;
00077   const unsigned int  *at=a->t,*bt=b->t; 
00078   const unsigned char *al=a->l,*bl=b->l;
00079   assert (a->dim==b->dim) ;
00080   assert (a->dim>=0) ;
00081   for (i=a->dim-1; i>=0; i--) {
00082     if (at[i]<bt[i]) return true ;
00083     if (at[i]>bt[i]) return false ;
00084     //
00085     if (al[i]<bl[i]) return true ;
00086     if (al[i]>bl[i]) return false ;
00087   }
00088   return false ;
00089 }
00090 
00092 //                        //
00093 // AdaptiveLevelsPointer  //
00094 //                        //
00096 
00097 AdaptiveLevelsPointer::AdaptiveLevelsPointer() {
00098   s =0 ;
00099   a =NULL ;
00100   lt=NULL ;
00101 }
00102 
00103 void AdaptiveLevelsPointer::Copy(AdaptiveLevelsPointer *alp) {
00104   s =alp->s ;
00105   a =alp->a ;
00106   lt=alp->lt ;
00107 }
00108 
00109 // compare for AdaptiveLevelsPointer
00110 int AdaptiveLevelsPointercmp(const void *pa,const void *pb) {
00111   AdaptiveLevelsPointer *a=(AdaptiveLevelsPointer *)pa ,
00112                         *b=(AdaptiveLevelsPointer *)pb ;
00113   if (a->s <  b->s) return +1 ;
00114   if (a->s == b->s) return  0 ;
00115   return -1 ;
00116 }
00117 
00118 // return number of bits
00119 int nbits(const unsigned long x) {
00120  int n=0 ;
00121  unsigned long k=1 ;
00122  while (k) { if (x&k) n++ ; k=k*2 ; }
00123  return n ;
00124 }
00125 
00126 // special compare relation for unsigned longs
00127 int lcmp(const void *pa,const void *pb) {
00128   unsigned long a=*((unsigned long *)pa) , b=*((unsigned long *)pb) ; 
00129   int na=nbits(a) , nb=nbits(b) ;
00130   if (na<nb) return -1 ;
00131   if (na>nb) return  1 ;
00132   if ((a)<(b)) return -1 ;
00133   if ((a)>(b)) return  1 ;
00134   return 0 ;
00135 }
00136 
00138 //          //
00139 //  SMPjob  //
00140 //          //
00142 
00143 SMPjob::SMPjob() {
00144   nalp=0 ;
00145   alp=NULL ;
00146   AB1=AB2=NULL ;
00147   AGS=AG2=NULL ;
00148   op=dir=DIM=MAXINT ;
00149   WC=NULL ;
00150 
00151   in=out=wenoroe=NULL ;
00152 
00153   finished  =false ;
00154   terminate =false ;
00155   poll      =true  ;
00156   jobcount  =0     ; // id of current SMP job to be processed ;
00157   requests  =0     ;
00158 }
00159 
00160 void SMPjob::set(int operation, int direction, int dim,
00161                  AdaptiveLevelsPointer *alpstart, size_t num_alp,
00162                  vector<double> **input, 
00163                  vector<double> **output, 
00164                  vector<double> **roespeed) {
00165 
00166     op  =operation ;
00167     dir =direction ;
00168     DIM =dim       ;
00169     in  =input     ;
00170     out =output    ;
00171     alp =alpstart  ;
00172     nalp=num_alp   ;
00173   
00174     wenoroe  =roespeed ;
00175     finished =false    ;
00176     terminate=false    ;
00177 }
00178 
00179 bool SMPjob::isfinished() { return finished ;}
00180 
00181 // copy of template<class I,int DIM>
00182 // inline I& Matrix<I,DIM>::operator[](const int *in)
00183 // the underlying matrix is a [0...LMAX]^DIM matrix
00184 size_t SMPjob::leveloffset(int *l) {
00185   size_t i=l[DIM-1], c=(LMAX-0+1) ;
00186   for (int k=DIM-2; k>=0; k--) {
00187     i += l[k]*c ;
00188     c *= (LMAX-0+1) ;
00189   }
00190   return i ;
00191 }
00192 
00193 void SMPjob::setpoll() { poll=true   ;}
00194 void SMPjob::clrpoll() { poll=false  ;}
00195 bool SMPjob::getpoll() { requests++; return poll ;}
00196 
00197 void *ProcessSMPjob(void *jb) {
00198   SMPjob *job=(SMPjob *)jb ;
00199 
00200   job->jobcount++ ; // signal catched !
00201 
00202   size_t ac ;
00203   int   K0,K1,ld,i ;
00204   // local copies
00205   size_t num_alp=job->nalp;
00206   int dir       =job->dir ;
00207   int DIM       =job->DIM ;
00208   int op        =job->op  ;
00209   Wavelets *WC  =job->WC  ;
00210 
00211   vector<double> *ba ;
00212   AdaptiveB      *ABX, *AB1=job->AB1, *AB2=job->AB2 ;
00213   AdaptiveGS     *AGS      =job->AGS, *AG2=job->AG2 ; 
00214 
00215   for (ac=0; ac<num_alp; ac++) {
00216     AdaptiveLevels        *als = job->alp[ac].a ;
00217     AdaptiveLevelsIndex   *ali = job->alp[ac].lt ;
00218  
00219     // extract level
00220     int l[DMAX] ; // should be enough
00221     assert ((size_t)DIM < sizeof(l)/sizeof(l[0])) ;
00222     for (i=0  ; i<dir  ; i++) l[i  ]=(int)ali->l[i] ;
00223     for (i=dir; i<DIM-1; i++) l[i+1]=(int)ali->l[i] ;
00224 
00225     // get Data and apply operation
00226     if (op==PROJECTION) { // optimized code for PROJECTION
00227 
00228       // see Adaptive1DOperations.cc:
00229       bool hack=(WC->InterpoletFlag && (WC->N==2)) ;
00230         
00231       // load 1D adaptive basis, but only the boundary part 
00232       for (ld=(int)AB1->Level0; ld<=LMAX; ld++) {
00233         if (ld > AB1->Level0) { 
00234           K0=WC->KW0 ; 
00235           K1=(1<<(ld-1)) -1 - WC->KW1 ; 
00236         } else { 
00237           K0=WC->K0  ;
00238           K1=(1<<ld)        - WC->K1  ; 
00239          }
00240         l[dir]=ld ;
00241         // data from where ?
00242         ba = job->in[job->leveloffset(l)] ;
00243         if (hack) AB1->L[ld].IS->Load        ((*als)[ld] , ba , AB1->L[ld].d ) ;
00244         else      AB1->L[ld].IS->LoadBoundary((*als)[ld] , ba , AB1->L[ld].d , K0 , K1 ) ;
00245       }
00246          
00247       // projection 
00248       AB2->Projection(AB1 , WC , AGS , AG2) ;
00249         
00250       // store boundary part of basis
00251       for (ld=AB1->Level0; ld<=LMAX; ld++) {
00252         if (ld > AB1->Level0) { 
00253           K0= WC->KW0 ;
00254           K1=(1<<(ld-1)) -1 - WC->KW1 ; 
00255          } else { 
00256           K0= WC->K0  ; 
00257           K1=(1<<ld)      -  WC->K1  ; 
00258          }
00259         l[dir]=ld ;
00260         // data to ?
00261         ba = job->out[job->leveloffset(l)] ;
00262         if (hack) AB1->L[ld].IS->Store        ((*als)[ld] , ba , AB2->L[ld].d ) ;
00263         else      AB1->L[ld].IS->StoreBoundary((*als)[ld] , ba , AB2->L[ld].d , K0 , K1 ) ;
00264       }
00265       
00266     } else {  // other Operations
00267       
00268       // load 1D adaptive basis 
00269       for (ld=AB1->Level0; ld<=LMAX; ld++) {
00270         l[dir]=ld ;
00271         // data from where ?
00272         ba = job->in[job->leveloffset(l)] ;
00273         AB1->L[ld].IS->Load  ((*als)[ld] , ba , AB1->L[ld].d ) ; 
00274       }
00275 
00276       ApplyUnivariateOperators(AB1, AB2, &ABX, AGS, AG2, op, WC) ;
00277  
00278       // store result
00279       for (ld=AB1->Level0; ld<=LMAX; ld++) {
00280         l[dir]=ld ;
00281         ba = job->out[job->leveloffset(l)] ;
00282         ABX->L[ld].IS->Store((*als)[ld], ba , ABX->L[ld].d ) ;
00283       } // next ld  
00284     }   // other Operations 
00285   }
00286 
00287   // unlock
00288   job->finished=true ;
00289 
00290   return NULL ;
00291 }
00292 
00293 
00294 
00295 
00296 void *ProcessSMPjobWENO(void *jb) {
00297   SMPjob *job=(SMPjob *)jb ;
00298 
00299   job->jobcount++ ; // signal catched !
00300 
00301   size_t ac ;
00302   int    i,ld ;
00303 
00304   // local copies
00305   size_t num_alp=job->nalp;
00306   int dir       =job->dir ;
00307   int DIM       =job->DIM ;
00308   int op        =job->op  ;
00309   Wavelets *WC  =job->WC  ;     
00310 
00311   vector<double> *ba ;
00312   AdaptiveB      *AB1=job->AB1, *AB2=job->AB2 ;
00313   AdaptiveGS     *AGS=job->AGS, *AG2=job->AG2 ;
00314 
00315   for (ac=0; ac<num_alp; ac++) {
00316     AdaptiveLevels        *als = job->alp[ac].a ;
00317     AdaptiveLevelsIndex   *ali = job->alp[ac].lt ;
00318 
00319     // extract level
00320     int l[10] ; // should be enough
00321     assert ((size_t)DIM < sizeof(l)/sizeof(l[0])) ;
00322     for (i=0  ; i<dir  ; i++) l[i  ]=(int)ali->l[i] ;
00323     for (i=dir; i<DIM-1; i++) l[i+1]=(int)ali->l[i] ;
00324       
00325     // load 1D adaptive basis 
00326     for (ld=AB1->Level0; ld<=LMAX; ld++) {
00327       l[dir]=ld ;
00328       // data from where ?
00329       ba = job->in[job->leveloffset(l)] ;
00330       AB1->L[ld].IS->Load  ((*als)[ld] , ba , AB1->L[ld].d ) ; 
00331       
00332       ba = job->wenoroe[job->leveloffset(l)] ;
00333       AB2->L[ld].IS->Load  ((*als)[ld] , ba , AB2->L[ld].d ) ; 
00334     }
00335 
00336     // IndexSet for WENO-FD
00337     AGS->IndexSetForAdaptiveGS(AB1, WC) ;
00338 
00339     // prepare WENO
00340     switch (op) {
00341     case WENO3:
00342       AG2->IndexSetForAdaptiveGSFD(AB1, &WC->FDOperators[WC->FDOp(FD40)],  WC) ;
00343       break ;
00344     case WENO5:
00345       AG2->IndexSetForAdaptiveGSFD(AB1, &WC->FDOperators[WC->FDOp(FD60)],  WC) ;
00346       break ;
00347     default: 
00348       assert(0) ;
00349       break ;
00350     }
00351 
00352     // calc RoeSpeed GS
00353     AG2->RestoreFromInterpoletBasis(AB2 , AGS , WC) ;
00354     // Patch for non-validated RoeSpeed part ....
00355     //for (ld=G->Level0[dir]; ld<=LMAX;ld++) G->AG2->L[ld].I.VecFill(G->AG2->L[ld].d,1.0) ;
00356 
00357     // GS for FD Operators required for auxiliary GS 
00358     AGS->CopyIndexSets(AG2) ;
00359 
00360     // WENO
00361     AB2->ApplyWENO(AB1 , op , WC, AGS, AG2) ;
00362       
00363     // store result
00364     for (ld=AB1->Level0; ld<=LMAX; ld++) {
00365       l[dir]=ld ;
00366       ba = job->out[job->leveloffset(l)] ;
00367       AB2->L[ld].IS->Store((*als)[ld], ba , AB2->L[ld].d ) ;
00368     } // next ld
00369   } // next ac
00370 
00371   // unlock
00372   job->finished=true ;
00373 
00374   return NULL ;
00375 }
Generated at Mon Aug 19 10:02:31 2002 for AWFD by