NonAdaptive.cc

00001 /*
00002 
00003   Wavelet-Routines 
00004 
00005 */
00006 
00007 #include "NonAdaptive.hpp"
00008 #include "Buffer.hpp"
00009 
00010 /*********************************************/
00011 /*                                           */
00012 /*    1D        Transformationen             */
00013 /*                                           */
00014 /*********************************************/
00015 
00016 
00017 //
00018 // orthogonale Wavelets
00019 //
00020 void bwt(double *d,int *BCd , int *dec , double *c,int *BCc , Wavelets *WC,int J,int lev , double)
00021 {
00022  int i,j,nJ,nJ2,a,e ;
00023  int bc0=BCc[0] , bc1=BCc[1] ;
00024  int L  =WC->L  , R  =WC->R  ;
00025  int K0 =WC->K0 , K1 =WC->K1 , KW0=WC->KW0 , KW1=WC->KW1 ; 
00026 
00027  BCd[0]=bc0 ; BCd[1]=bc1 ;
00028 
00029  double *d0,s ;
00030 
00031  assert(lev>=0) ;
00032  assert(bc0>=-1) ;
00033  assert(bc1>=-1) ;
00034 
00035  if (strcmp(WC->type,"Interpolet")==0) { hbt(d,BCd, dec,c,BCc, WC,J,lev ,0 ) ;return ;}
00036   
00037  nJ=(1<<J)+1 ;
00038 
00039  /******************************************** 
00040     if no further decomposition alowed, stop */
00041  if ( ((nJ/2+1) < WC->MX)  || (lev==0)) { for (i=0 ;i<nJ ;i++) d[i]=c[i] ;*dec=0 ;return ;}
00042 
00043 
00044  // allocate Memory
00045  double *c0=new double[(1<<J)+1] ;
00046 
00047  /********************************************
00048     perform one further decomposition        */  
00049  nJ2=(1<<(J-1))+1 ;
00050  d0 =&d[nJ2] ;
00051   
00052  /******************************************** 
00053     WaveletTrafo                             */
00054  if (bc1 !=PERIODIC) {
00055  
00056   /* c0=Hh0*c */
00057    WC->HJ.AL[bc0+1].MatVec(c,c0) ;
00058  
00059   /* c0 += Hh1*c */
00060    WC->HJ.AR[bc1+1].MatVec(&c[nJ-WC->HJ.AR[bc1+1].size[1]],&c0[nJ2-WC->HJ.AR[bc1+1].size[0]]) ;
00061    
00062   /* c0 += h*c */
00063    for (i=K0 ;i<=(1<<(J-1))-K1 ;i++) {
00064      a= (       K0 >= 2*i-L) ?        K0 : 2*i-L ;
00065      e= ((1<<J)-K1 <= 2*i+R) ? (1<<J)-K1 : 2*i+R ;
00066      s=0 ;for (j=a ;j<=e ;j++) s+=WC->HJ.a[j-2*i+L]*c[j] ;
00067      c0[i]=s ;
00068     }
00069 
00070   /* d=Gg0*c   */
00071    WC->GJ.AL[bc0+1].MatVec(c,d0) ;
00072 
00073   /* d +=Gg1*c */
00074    WC->GJ.AR[bc1+1].MatVec(&c [nJ-WC->GJ.AR[bc1+1].size[1]],&d0[(1<<(J-1))-WC->GJ.AR[bc1+1].size[0]]) ;
00075 
00076   /* d0 += g*c */
00077    for (i=KW0 ;i<(1<<(J-1))-KW1 ;i++) {
00078      a= (K0        >= 2*i-R+1) ?        K0 : 2*i-R+1 ;
00079      e= ((1<<J)-K1 <= 2*i+L+1) ? (1<<J)-K1 : 2*i+L+1 ;
00080      s=0 ; for (j=a ; j <=e ;j++) s+=WC->GJ.a[j-2*i+R-1]*c[j] ;
00081      d0[i]=s ; 
00082     }
00083    
00084  } else { /* PERIODIC */
00085    
00086    double *h,*g ;
00087    int    a,e,nJ1=1<<(J-1)   ;
00088    int    MM=max((L+1)/2,R/2) , MP=MM+1 ;
00089    
00090    nJ=(1<<J) ;
00091         
00092    for (i=0 ;i < MM ;i++) {
00093      s=0 ;a=2*i-L; e=2*i+R  ;h=&WC->HJ.a[L-2*i];
00094      for (j=a  ; j<=e  ; j++) s+=h[j]*c[(j+128*nJ)%nJ]; c0[i]=s ;
00095  
00096      s=0 ;a=2*i-R+1; e=2*i+L+1;g=&WC->GJ.a[R-1-2*i] ; 
00097      for (j=a; j<=e; j++) s+=g[j]*c[(j+128*nJ)%nJ]; d0[i]=s ; 
00098     }
00099 
00100    for (i=MM ;i <= nJ1-MP ;i++) {
00101      s=0 ;a=2*i-L; e=2*i+R  ;h=&WC->HJ.a[L-2*i];
00102      for (j=a  ; j<=e  ; j++) s+=h[j]*c[j]; c0[i]=s ;
00103  
00104      s=0 ;a=2*i-R+1; e=2*i+L+1;g=&WC->GJ.a[R-1-2*i] ; 
00105      for (j=a; j<=e; j++) s+=g[j]*c[j]; d0[i]=s ; 
00106     }
00107    
00108 
00109    for (i=nJ1-MP+1 ;i < nJ1 ;i++) {
00110      s=0 ;a=2*i-L; e=2*i+R  ;h=&WC->HJ.a[L-2*i];
00111      for (j=a  ; j<=e  ; j++) s+=h[j]*c[(j+128*nJ)%nJ]; c0[i]=s ;
00112  
00113      s=0 ;a=2*i-R+1; e=2*i+L+1;g=&WC->GJ.a[R-1-2*i] ; 
00114      for (j=a; j<=e; j++) s+=g[j]*c[(j+128*nJ)%nJ]; d0[i]=s ; 
00115     }
00116    c0[nJ1]=0.0 ;
00117  }
00118 
00119  bwt(d,BCd,dec,c0,BCc,WC,J-1,lev-1 , 0) ;
00120 
00121  (*dec)++ ;
00122  delete []c0 ;
00123 }
00124 
00125 
00126 void ibwt(double *c,int *BCc , int *, double *d,int *BCd , Wavelets *WC,int J,int lev, double)
00127 {
00128  int i,j,nJ,nJ2,a,e ;
00129 
00130  int bc0=BCc[0] , bc1=BCc[1] ;
00131  int L  =WC->L  , R  =WC->R  ;
00132  int K0 =WC->K0 , K1 =WC->K1 , KW0=WC->KW0 , KW1=WC->KW1 ;
00133 
00134  double *d0,s ;
00135 
00136  assert((lev>=0) && (bc0>=-1) && (bc1>=-1)) ;
00137 
00138  BCc[0]=bc0 ; BCc[1]=bc1 ;
00139 
00140  if (strcmp(WC->type,"Interpolet")==0) { ihbt(c,BCc,NULL,d,BCd,WC,J,lev ,0) ;return ;}
00141 
00142  nJ=(1<<J) ;
00143 
00144  if (lev==0) { for (i=0 ;i<=nJ ;i++) c[i]=d[i] ;return ;}
00145 
00146  // allocate memory and regular case
00147  double *c0=new double[(1<<J)+1];
00148 
00149  nJ2=(1<<(J-1))+1 ;
00150  d0 =&d[nJ2] ;
00151  
00152  ibwt(c0,BCc,NULL,d,BCd,WC,J-1,lev-1 ,0) ; 
00153 
00154  /*********************************** 
00155     inverse WaveletTrafo            */
00156  
00157  if (bc1 !=PERIODIC) { 
00158    for (i=0 ;i<=nJ ;i++) c[i]=0 ;
00159 
00160    WC->HJ.AL[bc0+1].APlusBMatTVec(1.0,1.0,c0,c) ;
00161    WC->HJ.AR[bc1+1].APlusBMatTVec(1.0,1.0,&c0[nJ2-WC->HJ.AR[bc1+1].size[0]],
00162                                   &c[nJ+1-WC->HJ.AR[bc1+1].size[1]]) ;
00163 
00164    WC->GJ.AL[bc0+1].APlusBMatTVec(1.0,1.0,d0,c) ;
00165    WC->GJ.AR[bc1+1].APlusBMatTVec(1.0,1.0,&d0[(1<<(J-1))-WC->GJ.AR[bc1+1].size[0]],
00166                                   &c[nJ+1-WC->GJ.AR[bc1+1].size[1]]) ;
00167 
00168    for (i=K0 ;i<=(1<<J)-K1 ;i++) {
00169        s=c[i] ;
00170        a=(int)((i-R+1)/2) ; a = (a<=           K0) ?            K0 : a ; 
00171        e=(int)((i+L  )/2) ; e = (e>=(1<<(J-1))-K1) ? (1<<(J-1))-K1 : e ; 
00172        for (j=a ;j<=e ;j++) s += WC->HJ.a[i-2*j+L]*c0[j] ;
00173 
00174        a=(int)((i-L  )/2) ; a = (a<=             KW0) ?              KW0 : a ; 
00175        e=(int)((i+R-1)/2) ; e = (e>=(1<<(J-1))-1-KW1) ? (1<<(J-1))-1-KW1 : e ; 
00176        for (j=a ;j<=e ;j++) s += WC->GJ.a[i-2*j+R-1]*d0[j] ; 
00177        c[i]=s ;
00178      }
00179 
00180  } else { /* PERIODIC */
00181  
00182    double *h,*g  ;
00183    int    MM=max(L,R-1) ; 
00184    
00185    nJ2=(1<<(J-1)) ;
00186 
00187    for (i=0 ;i<MM ;i++) {
00188      s=0.0 ;h=&WC->HJ.a[i+L] ; g=&WC->GJ.a[i+R-1] ;
00189 
00190      for (j=floor2(i-R+1); j<=(i+L  )/2; j++) 
00191          s+=h[-2*j]*c0[(j+128*nJ2)%nJ2] ;
00192      for (j=floor2(i-L  );j<=(i+R-1)/2; j++) 
00193          s+=g[-2*j]*d0[(j+128*nJ2)%nJ2] ;
00194      c[i]=s ;
00195     }  
00196 
00197    for (i=MM ;i<=nJ-2-MM ;i++) {
00198      s=0.0 ; h=&WC->HJ.a[i+L] ; g=&WC->GJ.a[i+R-1] ;
00199 
00200      for (j=(i-R+1)/2; j<=(i+L)/2   ;j++) 
00201          s+=h[-2*j]*c0[j] ;
00202      for (j=(i-L  )/2; j<=(i+R-1)/2 ;j++) 
00203          s+=g[-2*j]*d0[j] ;
00204      c[i]=s ;
00205     }  
00206 
00207    for (i=nJ-1-MM ;i<nJ ;i++) {
00208      s=0.0 ; h=&WC->HJ.a[i+L]  ; g=&WC->GJ.a[i+R-1] ;
00209 
00210      for (j=floor2(i-R+1); j<=(i+L  )/2; j++) 
00211          s+=h[-2*j]*c0[(j+128*nJ2)%nJ2] ;
00212      for (j=floor2(i-L  ); j<=(i+R-1)/2;j++) 
00213          s+=g[-2*j]*d0[(j+128*nJ2)%nJ2] ;
00214      c[i]=s ;
00215     }
00216    c[nJ]=0.0 ;
00217  }
00218 
00219  delete []c0 ;
00220 }
00221 
00222 
00223 //
00224 // Interpolets 
00225 //
00226 void hbt(double *d,int *BCd , int *dec , double *c,int *BCc , Wavelets *WC,int J,int lev, double)
00227 {
00228  int i,j,nJ,nJ2 ;
00229  int bc0=BCc[0] , bc1=BCc[1] ;
00230  int L  =WC->GtJ.U , R=WC->GtJ.O , KW=WC->KW0 ; 
00231 
00232  double *d0,s,*g=WC->GtJ.a ;
00233  
00234  BCd[0]=bc0 ; BCd[1]=bc1 ;
00235 
00236  nJ=(1<<J) ;
00237 
00238  /******************************************** 
00239     if no further decomposition alowed, stop */
00240  if ( ((nJ/2) < WC->MX)  || (lev==0)) {
00241    double q=(WC->LiftingFlag) ? sqrt((double)(1<<J)) : 1 ;
00242    for (i=0 ;i<=nJ ;i++) d[i]=c[i]/q ;
00243    *dec=0 ;
00244    return ;
00245   }
00246 
00247 
00248  /********************************************
00249     perform one further decomposition        */  
00250  double *c0=new double[(1<<J)+1] ;
00251 
00252  nJ2=(1<<(J-1)) ;
00253  d0 =&d[nJ2+1]  ;
00254   
00255  /******************************************** 
00256     WaveletTrafo                             */
00257 
00258  // Tiefpassfilter
00259  for (i=0; i<=nJ2; i++) c0[i]=c[2*i] ;
00260  if (bc0==1) c0[1    ]=0 ;
00261  if (bc1==1) c0[nJ2-1]=0 ;
00262  
00263  // Hochpassfilter 
00264  if (bc1==PERIODIC) {
00265    for (i=0; i<nJ2; i++) {
00266      s=c[2*i+1] ;
00267      for (j=-L;j<=R ;j+=2) s+=g[j+L]*c[(2*i+j+nJ)%nJ] ;
00268      d0[i]=s ;
00269    }
00270  } else { // nicht periodisch
00271    for (i=KW; i<nJ2-KW; i++) {
00272      s=c[2*i+1] ;
00273      for (j=-L; j<=R; j+=2) s+=g[j+L]*c[2*i+j] ;
00274      d0[i]=s ;
00275    }
00276 
00277    WC->GtJ.AL[bc0+1].MatVec(c,d0) ;
00278    WC->GtJ.AR[bc1+1].MatVec(&c[nJ-WC->GtJ.AR[bc1+1].size[1]+1],&d0[nJ2-WC->GtJ.AR[bc1+1].size[0]]) ;
00279   
00280  }
00281 
00282  // Lifting
00283  if (WC->LiftingFlag) {
00284    WC->Q.MatTVec(c0 , d0 , BCc , J-1 , -1) ;
00285    double q=sqrt((double)(1<<J)) ;
00286    for (i=0; i<nJ2; i++) d0[i]/=q ;
00287  } 
00288 
00289  hbt(d,BCd,dec,c0,BCc,WC,J-1,lev-1 , 0) ;
00290 
00291  (*dec)++ ;
00292  delete []c0 ;
00293 }
00294 
00295 
00296 void ihbt(double *c,int *BCc , int *, double *d,int *BCd , Wavelets *WC,int J,int lev, double)
00297 {
00298  int i,j,nJ,nJ2 ;
00299 
00300  int bc0=BCd[0] , bc1=BCd[1] ;
00301  int L  =WC->HJ.U  , R=WC->HJ.O , K=WC->K0 ;
00302 
00303  double *d0,s,*g=WC->HJ.a ;
00304  BCc[0]=bc0 ; BCc[1]=bc1 ;
00305 
00306  nJ=(1<<J) ;
00307 
00308  if (lev==0) {
00309    double q=(WC->LiftingFlag) ? sqrt((double)(1<<J)) : 1 ;
00310    for (i=0 ;i<=nJ ;i++) c[i]=d[i]*q ;
00311    return ;
00312   }
00313 
00314 
00315  double *c0=new double[(1<<J)+1];
00316  nJ2=(1<<(J-1)) ;
00317  d0 =&d[nJ2+1]  ;
00318  
00319  ihbt(c0,BCc,NULL,d,BCd,WC,J-1,lev-1 , 0) ; 
00320 
00321  double *dd=d0 ;
00322 
00323  // Undo Lifting
00324  if (WC->LiftingFlag) {
00325    double q= sqrt((double)(1<<J)) ;
00326    dd=new double[nJ2] ;
00327    for (i=0; i<nJ2; i++) dd[i]=d0[i]*q ;
00328    WC->Q.MatTVec(c0 , dd , BCc , J-1 , +1) ;
00329  }
00330 
00331  /*********************************** 
00332     inverse WaveletTrafo            */
00333  
00334  for (i=0;i<=nJ;i++) c[i]=0 ;
00335 
00336  // Rekonstruktion der low-level Skalierungsfunktionen
00337  //
00338 
00339  if (bc1==PERIODIC) {
00340    for (i=0; i<nJ2; i++) {
00341      c[2*i]=s=c0[i] ;
00342      for (j=-L;j<=R ;j+=2) c[(2*i+j+nJ)%nJ] +=g[j+L]*s ;
00343    }
00344    c[nJ]=0 ;
00345  } else { // nicht periodisch
00346    for (i=K; i<=nJ2-K; i++) {
00347      c[2*i]=s=c0[i] ;
00348      for (j=-L;j<=R ;j+=2) c[2*i+j] +=g[j+L]*s ;
00349    }
00350 
00351    WC->HJ.AL[bc0+1].APlusBMatTVec(1.0,1.0,c0,c) ;
00352    WC->HJ.AR[bc1+1].APlusBMatTVec(1.0,1.0,
00353                                   &c0[nJ2-WC->HJ.AR[bc1+1].size[0]+1],
00354                                   &c [nJ -WC->HJ.AR[bc1+1].size[1]+1]) ;
00355 
00356   }
00357 
00358  // Rekonstruktion der Wavelets
00359  //
00360  for (i=1;i<nJ2-1;i++) c[2*i+1]+=dd[i] ;
00361  if (bc0==1) c[2]+=dd[0] ;
00362     else     c[1]+=dd[0] ;
00363  
00364  if (bc1==1) c[nJ-2]+=dd[nJ2-1] ;
00365     else     c[nJ-1]+=dd[nJ2-1] ;
00366 
00367  delete []c0 ;
00368  if (WC->LiftingFlag) delete []dd ;
00369 }
00370 
00371 void ihbtT(double *d,int *BCd , int *dec , double *c,int *BCc , Wavelets *WC,int J,int lev, double)
00372 {
00373  int i,j,nJ,nJ2 ;
00374  int bc0=BCc[0] , bc1=BCc[1] ;
00375  int L  =WC->GtJ.U , R=WC->GtJ.O , KW=WC->KW0 ; 
00376 
00377  double *d0,s,*g=WC->GtJ.a ;
00378 
00379  assert((BCd[1]==PERIODIC) && (BCc[1]==PERIODIC)) ;
00380 
00381  BCd[0]=bc0 ; BCd[1]=bc1 ;
00382 
00383  nJ=(1<<J) ;
00384 
00385  /******************************************** 
00386     if no further decomposition alowed, stop */
00387  if ( ((nJ/2) < WC->MX)  || (lev==0)) { for (i=0 ;i<=nJ ;i++) d[i]=c[i] ;*dec=0 ;return ;}
00388 
00389  /********************************************
00390     perform one further decomposition        */
00391  double *c0=new double[(1<<J)+1];
00392  nJ2=(1<<(J-1)) ;
00393  d0 =&d[nJ2+1]  ;
00394   
00395  /******************************************** 
00396     WaveletTrafo                             */
00397 
00398  // Tiefpassfilter
00399  for (i=0; i<nJ2; i++) d0[i]=c[2*i+1] ;
00400  
00401  // Hochpassfilter 
00402  if (bc1==PERIODIC) {
00403    for (i=0; i<nJ2; i++) {
00404      s=c[2*i] ;
00405      for (j=-L;j<=R ;j+=2) s-=g[j+L]*c[(2*i+1-j+nJ)%nJ] ;
00406      c0[i]=s ;
00407    }
00408  } else { // nicht periodisch
00409    for (i=KW; i<nJ2-KW; i++) {
00410      s=c[2*i+1] ;
00411      for (j=-L; j<=R; j+=2) s+=g[j+L]*c[2*i+j] ;
00412      d0[i]=s ;
00413    }
00414 
00415    WC->GtJ.AL[bc0+1].MatVec(c,d0) ;
00416    WC->GtJ.AR[bc1+1].MatVec(&c[nJ-WC->GtJ.AR[bc1+1].size[1]+1],&d0[nJ2-WC->GtJ.AR[bc1+1].size[0]]) ;
00417  }
00418 
00419  // Lifting
00420  if (WC->LiftingFlag) {
00421    switch (WC->N) {
00422     case 2:
00423       for (i=0 ;i<nJ2 ;i++)  d0[i] -=0.25*(c0[(i+1)%nJ2]+c0[i]) ;
00424       break ;
00425     case 4:
00426       break ;
00427     case 6:
00428       for (i=0 ;i<nJ2 ;i++) d0[i] -=( +2.929687499999999e-01*(c0[(i+1)%nJ2]+c0[i  ])
00429                                       -4.882812500000000e-02*(c0[(i+2)%nJ2]+c0[(i+nJ2-1)%nJ2])
00430                                       +5.859375000000001e-03*(c0[(i+3)%nJ2]+c0[(i+nJ2-2)%nJ2])) ;
00431       break ;
00432     }
00433  }
00434 
00435  ihbtT(d,BCd,dec,c0,BCc,WC,J-1,lev-1 , 0) ;
00436 
00437  (*dec)++ ;
00438 
00439  delete []c0 ;
00440 }
00441 
00442 
00443 //
00444 // Lifting prewavelets
00445 //
00446 void prwt  (double *d,int *BCd , int *dec, double *c,int *BCc , Wavelets * ,int J,int lev, double)
00447 {
00448  int    i,nJ=(1<<J), n2=(1<<(J-1)) ;
00449  double *c0,*d0=&d[n2+1] ;
00450  
00451  assert((BCd[1]==PERIODIC) && (BCc[1]==PERIODIC)) ;
00452 
00453  if (J<3) { for (i=0 ;i<=nJ ;i++) d[i]=c[i] ;*dec=0 ;return ;}
00454 
00455  c0=new double[n2] ;
00456 
00457 /* HB-Trafo   */
00458  for (i=0 ;i<nJ ;i +=2) c0[i/2]=c[i] ; 
00459  for (i=0 ;i<n2 ;i++  ) d0[i]  =c[2*i+1]-0.5*(c[(2*i+2)%nJ]+c[2*i]) ;
00460 
00461 /* Lifting */
00462  for (i=0 ;i<n2 ;i++) c0[i] +=0.25*(d0[(i+n2-1)%n2]+d0[i]) ;
00463  c0[n2]=0 ;
00464 
00465 /* Rekursiver Aufruf */
00466  prwt(d,BCd,dec,c0,BCc,NULL,J-1,lev-1 ,0) ;
00467 
00468  delete []c0 ;
00469  (*dec)++ ;
00470 }
00471 
00472 
00473 void iprwt (double *c,int *BCc , int *dec, double *d,int *BCd , Wavelets * ,int J,int lev, double)
00474 {
00475  int    i,nJ=(1<<J), n2=(1<<(J-1)) ;
00476  double *c0, *d0=&d[n2+1] ;
00477  
00478  assert((BCd[1]==PERIODIC) && (BCc[1]==PERIODIC)) ;
00479 
00480  if (lev==0) { for (i=0 ;i<=nJ ;i++) c[i]=d[i] ;return ;}
00481 
00482 /* Rekursiver Aufruf */
00483  c0=new double[n2] ;
00484 
00485  iprwt(c0,BCc,dec,d,BCd,NULL,J-1,lev-1 , 0) ;
00486 
00487 /* Undo Lifting */
00488 
00489  for (i=0 ;i<n2 ;i++) c0[i] -=0.25*(d0[(i+n2-1)%n2]+d0[i]); 
00490 
00491  /* Undo HB Filter  */
00492 
00493  for (i=0 ;i< n2 ;i++) { c[2*i]=c0[i]; c[2*i+1]=d0[i]+0.5*(c0[i]+c0[(i+1)%n2]) ; }
00494  
00495  delete []c0 ;
00496 }
00497 
00498 
00499 void iprwtT(double *d,int *BCd , int *dec, double *c,int *BCc , Wavelets * ,int J,int lev, double)
00500 {
00501  int    i,nJ=(1<<J), n2=(1<<(J-1)) ;
00502  double *c0,*d0=&d[n2+1] ;
00503 
00504  assert((BCd[1]==PERIODIC) && (BCc[1]==PERIODIC)) ;
00505 
00506  if (lev==0) { for (i=0 ;i<=nJ ;i++) d[i]=c[i] ;*dec=0 ;return ;}
00507 
00508  c0=new double[n2] ;
00509 
00510 /* HB-Trafo   */
00511 
00512  for (i=0 ;i<n2 ;i++) d0[i]=c[2*i+1] ;
00513  for (i=0 ;i<n2 ;i++) c0[i]=c[2*i]+0.5*(d0[i]+d0[(i+n2-1)%n2]) ; 
00514 
00515 /* Lifting */
00516  for (i=0 ;i<n2 ;i++) d0[i] -=0.25*(c0[(i+1)%n2]+c0[i]) ;
00517 
00518 /* Rekursiver Aufruf */
00519 
00520  iprwtT(d,BCd,dec,c0,BCc,NULL,J-1,lev-1 , 0) ;
00521 
00522  delete c0 ;
00523 
00524  (*dec)++ ;
00525 }
00526 
00527 //
00528 //
00529 void massHat(int *BC , double *d , double *c , int J , Wavelets *WC, double ) {
00530  assert(BC[1]==PERIODIC) ;
00531  WC->MassMatrix.MatVec(d,BC, c,BC , J , 1.0) ;
00532 }
00533 
00534 void stiffHat(int *BC , double *d , double *c , int J , Wavelets *WC, double) {
00535  assert(BC[1]==PERIODIC) ;
00536  WC->StiffnessMatrix.MatVec(d,BC, c,BC , J , 1.0) ;
00537 }
00538 
00539 
00540 /****************************************************/
00541 /* The following quadratures are supported:         */
00542 /*                                                  */
00543 /*  op== 1 : polynomial exact Quadrature            */
00544 /*  op==-1 : exact inverse of polynomial Q.         */
00545 /*                                                  */
00546 /*  op== 2 : exact interpolation using nodal values */ 
00547 /*  op==-2 : evaluation of nodal values             */
00548 /*                                                  */
00549 /*  op== 3 : short quadrature on the boundary only, no scaling */
00550 /*  op==-3 : short inverse quadrature                          */
00551 /*                                                  */
00552 /*  op== 4 : onepoint quadrature                    */
00553 /*  op==-4 : inverse onepoint quadrature            */
00554 /*                                                  */
00555 /****************************************************/
00556 
00557 void Quadrature(double *d,int *BCd , int * , double *c,int *BCc , Wavelets *WC,int J,int op, double DX)
00558 {
00559  int i,j,nJ ,iop=-MAXINT ;
00560 
00561  int bc0=BCc[0] , bc1=BCc[1] ;
00562  
00563  TransformMatrix *W=NULL ;
00564  double   *w,s,wg,*AT0=NULL,*AT1=NULL ;
00565  int      la,le,K0=WC->K0,K1=WC->K1 ;
00566 
00567  BCd[0]=bc0 ; BCd[1]=bc1 ;
00568 
00569  /*********************************
00570             check op              */
00571 
00572   assert((op==    POLYQUAD)||(op==INVERSEPOLYQUAD)||
00573          (op==   NODALQUAD)||(op==INVERSENODALQUAD)||
00574          (op==   SHORTQUAD)||(op==INVERSESHORTQUAD)||
00575          (op==ONEPOINTQUAD)||(op==INVERSEONEPOINTQUAD)) ;
00576 
00577  /*********************************
00578           set values              */
00579 
00580  nJ=(1<<J) ;
00581  wg=(1<<J) ;
00582  wg=sqrt(wg) ;
00583 
00584  /**************************************************
00585    implicitely given quadratures as inverse of ... * 
00586    solution by preconditioned BiCGStab             */ 
00587   
00588  if ((op==NODALQUAD)||(op==INVERSEPOLYQUAD)) { 
00589 
00590    double *r = new double [nJ+1] ;
00591    double *s = new double [nJ+1] ;
00592    double *t = new double [nJ+1] ;
00593    double *x = new double [nJ+1] ;
00594    double *y = new double [nJ+1] ;
00595    double *z = new double [nJ+1] ;
00596    double *v = new double [nJ+1] ;
00597    double *p = new double [nJ+1] ;
00598    double *kt= new double [nJ+1] ;
00599    double *rd= new double [nJ+1] ;
00600    double *bs= new double [nJ+1] ;
00601 
00602    double rr=0.0,alpha,beta,omega,rv,rvo,rho,rho1,best=1e+100 ;
00603    int    it=0 ;
00604    
00605    int    PrecQuad=-MAXINT ;
00606 
00607  /*****************************
00608      select preconditioner    */
00609    
00610    if (op==      NODALQUAD) { PrecQuad=POLYQUAD         ; iop=INVERSENODALQUAD ; }
00611    if (op==INVERSEPOLYQUAD) { PrecQuad=INVERSENODALQUAD ; iop=POLYQUAD ; }
00612 
00613  /*****************************
00614    first guess                */
00615 
00616    Quadrature(x,BCd,NULL,c,BCc,WC,J,PrecQuad,DX) ;  
00617    
00618    Quadrature(r,BCd,NULL,x,BCc,WC,J,iop,DX) ;
00619    for (i=0; i<=nJ; i++) 
00620      r [i] =c[i]-r[i] , /* r=b-Ax */
00621      rr   +=r[i]*r[i] , 
00622      rd[i] =r[i]      ,
00623      v [i] =0.0       ,
00624      p [i] =0.0 ;              
00625 
00626    rho=alpha=omega=1 ;
00627       
00628    while ((rr>1e-28) && (it<30)) {
00629      it++ ;
00630      rho1=rho ;
00631      rho=0.0; for (i=0; i<=nJ; i++) rho += rd[i]*r[i] ;      /* rho=(rd,r) */
00632      beta=(rho*alpha)/(rho1*omega) ;
00633 
00634      for (i=0; i<=nJ; i++) p[i]=r[i]+beta*(p[i]-omega*v[i]) ; /* p=r+beta*(p-omega*v) */
00635            
00636      Quadrature(y,BCd,NULL,p,BCc,WC,J,PrecQuad,DX)  ;                   /* y=K^(-1)p */
00637      Quadrature(v,BCd,NULL,y,BCc,WC,J,iop,DX)  ;                   /* v=Ay      */
00638      
00639      rv   =0.0; 
00640      for (i=0; i<=nJ; i++) rv += rd[i]*v[i] ;                   /* alpha=rho/(rd,v) */
00641      alpha=rho/rv ;
00642 
00643      for (i=0; i<=nJ; i++) s[i]=r[i]-alpha*v[i] ;               /* s=r-alpha*v */
00644 
00645      Quadrature(z ,BCd,NULL,s,BCc,WC,J,PrecQuad,DX)  ;                  /* z=K^(-1)s  */
00646      Quadrature(t ,BCd,NULL,z,BCc,WC,J,iop,DX)       ;                  /* t=Az       */
00647      Quadrature(kt,BCd,NULL,t,BCc,WC,J,PrecQuad,DX)  ;                  /* kt=K^(-1)t */
00648 
00649      omega=rvo=0.0; 
00650      for (i=0; i<=nJ; i++) { 
00651        omega += kt[i]*z [i] ;
00652        rvo   += kt[i]*kt[i] ;
00653      }
00654      
00655      if ((fabs(omega)<1e-20) && (fabs(rvo)<1e-20))
00656        { 
00657          omega=rvo=1.0 ; }
00658 
00659      omega=omega/rvo ;
00660 
00661      rr=0.0; 
00662      for (i=0; i<=nJ; i++) { 
00663        x[i]=x[i]+alpha*y[i]+omega*z[i] ;
00664        r[i]=s[i]-omega*t[i] ; 
00665        rr+=r[i]*r[i] ;
00666      }
00667 
00668      if (rr<best) { 
00669        best=rr ;
00670        for (i=0; i<=nJ; i++) bs[i]=x[i] ;
00671      }
00672    }  
00673    
00674    if (best<=rr) for (i=0; i<=nJ; i++) d[i]=bs[i] ; 
00675            else  for (i=0; i<=nJ; i++) d[i]=x [i] ;
00676 
00677    delete []r  ;
00678    delete []s  ;
00679    delete []t  ;
00680    delete []x  ;
00681    delete []y  ;
00682    delete []z  ;
00683    delete []v  ;
00684    delete []p  ;
00685    delete []rd ;
00686    delete []kt ;
00687    delete []bs ;
00688 
00689    return ; 
00690  }
00691 
00692 
00693 /****************************************************
00694   Quadratures, which are explicitely given: op==1,-2 , 3,-3,4,-4 */
00695   
00696  if (op==POLYQUAD) {   /* polynomial Quadrature      */
00697    W=&WC->W ;
00698    wg=1/wg ;
00699  } 
00700  if (op==INVERSENODALQUAD) {   /* evaluation of nodal values */
00701    W=&WC->Phi ;
00702  }
00703 
00704  if (op==SHORTQUAD) {   /* polynomial Quadrature only at boundary     */
00705    W=&WC->W ;
00706    wg=1.0  ;
00707  }
00708 
00709  if (op==INVERSESHORTQUAD) {   /* polynomial inverse Quadrature only at boundary */
00710    W=&WC->IW ;
00711    wg=1.0 ;
00712  }
00713 
00714  if ((op==ONEPOINTQUAD)||(op==INVERSEONEPOINTQUAD)) {
00715    AT0=&WC->ATau0[bc0+1].a[0][0] ;
00716    AT1=&WC->ATau1[bc1+1].a[0][0] ;
00717    if (op==ONEPOINTQUAD) wg=1/wg ;
00718 
00719    W=&WC->IW ; // Dummy
00720  }
00721 
00722  w  =W->a ;
00723  la =W->U ;
00724  le =W->O ;
00725 
00726  // short quadrature just for the boundaries and no scaling 
00727  if ( (op==   SHORTQUAD)||(op==INVERSESHORTQUAD)||
00728       (op==ONEPOINTQUAD)||(op==INVERSEONEPOINTQUAD) ) {
00729    if (d!=c) for (i=0 ;i<=nJ; i++) d[i]=c[i]*wg ;
00730 
00731    if (bc1==PERIODIC) { d[nJ]=0.0 ; return ; }
00732 
00733    if ((op==SHORTQUAD)||(op==INVERSESHORTQUAD)) {
00734      W->AL[bc0+1].MatVec(c,d) ;    
00735      W->AR[bc1+1].MatVec(&c[nJ+1-W->AR[bc1+1].size[1]],&d[nJ+1-W->AR[bc1+1].size[0]]) ;
00736     }
00737    
00738    if (op== ONEPOINTQUAD) {
00739      for (i=0;i<K0;i++) d[i        ] /=AT0[i] ;
00740      for (i=0;i<K1;i++) d[nJ-K1+1+i] /=AT1[i] ;
00741    }
00742 
00743    if (op==INVERSEONEPOINTQUAD) {
00744      for (i=0;i<K0;i++) d[i        ] *=AT0[i] ;
00745      for (i=0;i<K1;i++) d[nJ-K1+1+i] *=AT1[i] ;
00746    }
00747 
00748    return ;
00749  }
00750 
00751  if (d==c) {
00752    std::cout << "1D nonadaptive quadrature: argument and result must not coincide\n";
00753    exit(-1) ;
00754  }
00755 
00756  if (bc1==PERIODIC) {
00757    if (op==POLYQUAD) 
00758      for (i=0 ;i <nJ; i++)
00759        { s=0 ;for (j=-la;j<=le;j++) s+=c[(j+i+128*nJ)%nJ]*w[j+la] ;
00760          d[i]=s*wg ;}
00761    else 
00762      for (i=0 ;i <nJ; i++)
00763        { s=0 ;for (j=-la;j<=le;j++) s+=c[(-j+i+128*nJ)%nJ]*w[j+la] ;   
00764          d[i]=s*wg ;}
00765    d[nJ]=0.0 ;
00766    return ;
00767   }
00768  
00769  for (i=0 ;i<=nJ; i++) d[i]=0.0 ;
00770 
00771  /* d=W0*c */
00772 
00773  W->AL[bc0+1].APlusBMatVec(0.0,wg,c,d) ;
00774  
00775  /* d = W1*c */
00776  W->AR[bc1+1].APlusBMatVec(0.0,wg,&c[nJ+1-W->AR[bc1+1].size[1]],&d[nJ+1-W->AR[bc1+1].size[0]]) ;
00777  
00778  if (op==POLYQUAD)
00779    for (i=K0 ;i <=(1<<J)-K1; i++) 
00780      { s=0; for (j=-la;j<= le;j++) s+=c[j+i]*w[j+la] ;      
00781        d[i]=s*wg ; }
00782  else 
00783    for (i=K0 ;i <=(1<<J)-K1; i++)
00784      { s=c[i]*wg; for (j=i-la ; j<=i+le ;j++)  d[j] +=s*w[j-i+la] ;}
00785 }
00786 
00787 
00788 /*********************************************/
00789 /*                                           */
00790 /*    1D    Operatoren                       */
00791 /*                                           */
00792 /* Folgende Operatoren werden unterstuetzt:  */
00793 /*                                           */
00794 /* op==1 : GalerkinProj. der ersten Ableitung*/
00795 /* op==2 : GalerkinSteifigkeitsmatrix        */
00796 /*                                           */
00797 /*********************************************/
00798 
00799 void applyOp1W(double *d,int *BCd , int * ,double *c,int *BCc , Wavelets *WC,int J,int op, double DX)
00800 {
00801  double *t=WC->Buffer->lock(); assert(t);
00802 
00803  ibwt   (d,BCd,NULL,c,BCc,WC,J,op, 0) ;
00804  applyOp(t,BCd,NULL,d,BCd,WC,J,1 ,DX) ;
00805  bwt    (d,BCc,NULL,t,BCd,WC,J,op, 0) ;
00806 
00807  WC->Buffer->unlock();
00808 }
00809 
00810 void applyOp(double *d,int *BCd , int * , double *c,int *BCc , Wavelets *WC,int J,int op, double DX)
00811 {
00812  assert(d!=c) ;
00813 
00814  // smooth operators
00815  if (op==SMOOTH0) {
00816    assert(BCc[1]==PERIODIC) ;
00817    int nJ=1<<J ;
00818    for (int i=1; i<nJ-1; i++) d[i]=0.25*(c[i-1]+2*c[i]+c[i+1]) ;
00819    d[0   ]=0.25*(c[nJ-1]+2*c[0   ]+c[1]) ;
00820    d[nJ-1]=0.25*(c[nJ-2]+2*c[nJ-1]+c[0]) ;
00821    return ;
00822  }
00823 
00824  //special Operators: stabilized div / grad
00825  if ((op==GRAD) || (op==DIV)) {
00826    double *f=WC->Buffer->lock(); assert(f); 
00827 
00828    if (BCd[1]>BCc[1]) {
00829      applyOp   (d,BCc,NULL,c,BCc, WC,J, (WC->N==6) ? FD80 : FD60 , DX) ;
00830      Projection(f,BCd,NULL,d,BCc, WC,J, 1 , 0)   ;
00831    } else { 
00832      Projection(d,BCd,NULL,c,BCc, WC,J, 1 , 0)   ;
00833      applyOp   (f,BCd,NULL,d,BCd, WC,J, (WC->N==6) ? FD80 : FD60 , DX) ;
00834    } 
00835 
00836    applyOp   (d,BCd,NULL,c,BCc, WC,J, FIRSTDERIVATIVE, DX) ;
00837  
00838    int i,nJ=1<<J ;
00839    double alpha=300.0 ;
00840    if (BCc[1]!=PERIODIC) alpha=300.0 ;
00841    if (op==GRAD) for (i=0; i<=nJ; i++) d[i] += alpha*f[i] ;
00842             else for (i=0; i<=nJ; i++) d[i] -= alpha*f[i] ; 
00843 
00844    WC->Buffer->unlock() ;
00845    return ;
00846  }
00847 
00848  // if FD-Op, Patch (homogeneous Dirichlet) BC
00849  int BCF[2]={-1,-1} , BCT[2]={-1,-1} ;
00850  if (BCc[1]==PERIODIC) BCF[1]=BCT[1]=PERIODIC ;
00851  if (BCc[0]==0) c[0   ]=0 ;
00852  if (BCc[1]==0) c[1<<J]=0 ;
00853 
00854  //
00855  // all other operators
00856  // 
00857 
00858  switch (op) {
00859  case FD12:    case FD14: case FD16: case FD22: case FD24: case FD40: case FD60: case FD80: case GRADFD4: case DIVFD4:
00860    WC->FDOperators[WC->FDOp(op)].MatVec(d,BCT  ,  c,BCF,  J , DX) ;
00861    break ;
00862  case FIRSTDERIVATIVE:     
00863    if ( WC->InterpoletFlag && (WC->N==2) ) {
00864      applyOp(d,BCd , NULL , c,BCc , WC,J, FD12 , DX) ;
00865      return ;
00866    } else {
00867      WC->Operators[0][0][0].MatVec(d,BCd  ,  c,BCc,  J , DX) ;
00868    }
00869    break ;
00870  case STIFFNESS: 
00871    if ( WC->InterpoletFlag && (WC->N<=4) ) {
00872      if (WC->N==2) applyOp(d,BCd , NULL , c,BCc , WC,J, FD22 , DX) ;
00873      if (WC->N==4) applyOp(d,BCd , NULL , c,BCc , WC,J, FD24 , DX) ;
00874      return ;
00875    } else {
00876      WC->Operators[1][0][0].MatVec(d,BCd  ,  c,BCc,  J , DX) ;
00877    }
00878    break ;
00879  case PROJECTION:
00880    Projection(d,BCd,NULL,c,BCc, WC,J, 0 , DX) ;
00881    break ;
00882  default: assert(0) ;
00883  }
00884  
00885  // patch honogeneous Dirichlet BC. This has only some effect for FD Operators !!
00886  if  (BCd[0]==0)                            d[0   ]=0 ;
00887  if ((BCd[1]==0)||(BCd[1]==PERIODIC)) d[1<<J]=0 ;
00888 }
00889 
00890 
00891 /*******************************************************/
00892 /*  Multiplikation im Raum der SkalierungsFunktionen   */
00893 /*******************************************************/
00894 
00895 void Multiply(double *uv,double *u,double *v,int *BC, struct Wavelets *W,int J,int which)
00896 {
00897  int nJ=1<<J ,i ,K0=W->K0 , K1=W->K1 ;
00898  double q=sqrt((double)(nJ)),*A ;
00899 
00900  // perform one pint quadrature for inner scaling f. and scale properly 
00901  for (i=0;i<=nJ;i++) uv[i]=u[i]*v[i]*q ;
00902  
00903  // periodic -> finish
00904  if (BC[1]==PERIODIC) return ;
00905 
00906  //
00907  // boundary treatment
00908  assert( (which==ONE_POINT_MULTIPLY) || (which==BOUND_MULTIPLY) ) ;
00909 
00910  //
00911  // modified one point rule 
00912  if (which==ONE_POINT_MULTIPLY) {
00913    A=&W->ATau0[BC[0]+1].a[0][0] ;
00914    for (i=0;i<K0;i++) uv[i]*= A[i] ;
00915 
00916    A=&W->ATau1[BC[1]+1].a[0][0] ; 
00917    for (i=0;i<K1;i++) uv[nJ-K1+1+i] *= A[i] ;
00918  } 
00919  
00920  //
00921  // polynomial exact rule
00922  if (which==BOUND_MULTIPLY) {
00923 
00924    double DX=1.0 ;
00925    Quadrature(u,BC,NULL,u,BC,W,J,-3 , DX) ;
00926    if (u!=v) Quadrature(v,BC,NULL,v,BC,W,J,-3, DX) ;
00927   
00928    int d=nJ-K1+1 ;
00929    for (i=0;i<K0;i++) uv[i  ]=u[i]*v[i]*q ;
00930    for (i=0;i<K1;i++) uv[d+i]=u[d+i]*v[d+i]*q ;
00931 
00932    Quadrature(u,BC, NULL, u,BC,W,J, 3 , DX) ;
00933    if (u!=v) Quadrature(v,BC, NULL,v,BC, W,J, 3 , DX) ;
00934    Quadrature(uv,BC, NULL,uv,BC ,W,J, 3 , DX) ;
00935  }
00936 
00937 }
00938 
00939 /*********************************************/
00940 /*                                           */
00941 /*            Projection                     */
00942 /*                                           */
00943 /*********************************************/
00944 
00945 
00946 void Projection(double *d,int *BCd , int * , double *c,int *BCc , Wavelets *WC,int J,int , double)
00947 {
00948  int i,nJF ;
00949  int K1 =WC->K1 ;
00950 
00951  int bcF0=BCc[0] , bcF1=BCc[1] ; // from
00952  int bcT0=BCd[0] , bcT1=BCd[1] ; // to
00953   
00954  nJF=(1<<J)+1 ;
00955 //**********************************
00956 // for periodic BC just simply copy 
00957  
00958  if ((bcF1==PERIODIC)||(bcT1==PERIODIC)) {
00959    assert(bcF1==bcT1) ;
00960    if (d!=c) for (i=0;i<nJF;i++) d[i]=c[i] ;
00961    d[1<<J]=0 ;
00962   }
00963  
00964 //**********************************
00965 //  projections for wavelets on the interval
00966  
00967  if (d!=c) for (i=0; i<=(1<<J); i++) d[i]=c[i] ;
00968 
00969  if (!WC->InterpoletFlag) {
00970    if ((bcF0>= 0)&&(bcT0==-1)) WC->OP0[bcF0].MatTVec(&c[0],&d[0]) ; 
00971    if ((bcF1>= 0)&&(bcT1==-1)) WC->OP1[bcF1].MatTVec(&c[(1<<J)+1-K1],&d[(1<<J)+1-K1]) ;
00972 
00973    if ((bcF0==-1)&&(bcT0>=0 )) WC->OP0[bcT0].MatVec(&c[0],&d[0]) ;
00974    if ((bcF1==-1)&&(bcT1>=0 )) WC->OP1[bcT1].MatVec(&c[(1<<J)+1-K1],&d[(1<<J)+1-K1]) ;
00975  } else {
00976    if ((bcF0>= 0)&&(bcT0==-1)) {WC->OP0[bcF0].MatVec(&c[0],&d[0])                     ; } 
00977    if ((bcF1>= 0)&&(bcT1==-1)) {WC->OP1[bcF1].MatVec(&c[(1<<J)+1-K1],&d[(1<<J)+1-K1]) ; }    
00978 
00979    if ((bcF0==-1)&&(bcT0>=0)) d[bcT0]=0 ;
00980    if ((bcF1==-1)&&(bcT1>=0)) d[(1<<J)-bcT1]=0 ;
00981  }
00982 }
00983 
00984 
00985 /*************************************/
00986 /*                                   */
00987 /* Sortier-Funktionen fuer HB Trafo  */
00988 /*                                   */
00989 /*************************************/
00990 
00991 void SortHB(double *d,int *, int * , double *c,int *, Wavelets *,int J,int lev0 , double)
00992 {
00993  int i,l,nL,n,o,f ;
00994  assert (d!=c) ;
00995 
00996  l=lev0 ;
00997  f=1<<(J-l) ;
00998  nL=1<<l    ;
00999  for (i=0; i<=nL; i++) d[i*f]=c[i] ;
01000 
01001  for (l=lev0+1; l<=J; l++) {
01002    f =1<<(J-l) ;
01003    nL=1<<l    ;
01004    n =1<<(l-1) ;
01005    o =n+1      ;
01006 
01007    for (i=0; i<n; i++) d[(2*i+1)*f]=c[i+o] ;
01008  }
01009 }
Generated at Mon Aug 19 10:02:32 2002 for AWFD by