Department of Scientific Computing    Institute for Numerical Simulation    University of Bonn

Documentation

Download

Programming References

Bug Reports / Suggestions

FAQ

Authors

           

---

Adaptive1DMatrix.cc

//: Adaptive1D.cc

#include "Adaptive1D.hpp"
#include "Buffer.hpp"

//
//  MatVec (double *T1,IndexSet *IT,int *BCT  ,  double *F,IndexSet *IF,int *BCF,int l)
//  multiplies *This(l)  by the adaptive vektor (F,IF,BCF). Only the the components from IT are calculated 
//  All other components are left unchanged
// 
//  In the current implementation, the complete matrix-vector--product is performed, and thereafter just the
//  components from IT are copied to the final result
//  This avoids the clumsy checks for 'j \in IT ?'
//

extern int debugRefine ;
void OperatorMatrix::MatVec (double *T1,IndexSet *IT,int *BCT  ,  double *F,IndexSet *IF,int *BCF,
                             int l,int clflag, double DX, doubleBuffer *buf)
{
 if (IT->nr==0) return ;

 Matrix2D *CL,*CLT,*CR,*CRT ;
 int       l2,M,N,l21       ;
 int       KL0,KL1,KR0,KR1,BL0,BL1,BR0,BR1 ;
 int       i,j,r,e,je,l24   ; 
 double    q=pow((1<<l)/DX , power),x,*b,*T=buf->lock(); 
 assert(T) ;

 CL=CLT=CR=CRT=NULL ;
 KL0=KL1=KR0=KR1=BL0=BL1=BR0=BR1=0 ;

 // get dimensions of (sub-)matrices
 l2 =M=N= 1<<l ; // {0,..,M}x{0,..,N}-matrix
 l21=l2-1 ;
 l24=4*l2 ;

 if (IsWT) M-- ;
 if (IsWF) N-- ;

 // clear result
 if (clflag) IT->VecClear(T)   ;
       else  IT->VecCopy(T,T1) ;
 
 // left boundary
 if (BCT[1]!=PERIODIC) {

   assert(BCF[1]!=PERIODIC) ;

   CL =&AL [BCT[0]+1][BCF[0]+1] ; // sub blocks
   CLT=&ALT[BCT[0]+1][BCF[0]+1] ;
   CR =&AR [BCT[1]+1][BCF[1]+1] ;
   CRT=&ART[BCT[1]+1][BCF[1]+1] ;

   KL0=CL->size[0] ; KL1=CLT->size[1] ;
   KR0=CR->size[0] ; KR1=CRT->size[1] ;

   BL1=CL->size[1]-1 ; BL0=KL0+CLT->size[0]-1 ;
   BR1=CR->size[1]-1 ; BR0=KR0+CRT->size[0]-1 ;

   assert((KL0+0+KR0)<=M) ;

   // very small matrices are build up explicitely 
   if ((BL0+1+KR0)>=M) {
     int in[2],im[2] ;
     Matrix2D A(M+1,N+1) ;
     // copy inner parts, CL and CR
     for (i=KL0; i<=M-KR0; i++) for (j=-O; j<=U  ; j++) { in[0]=i; in[1]=i+j; A[in]=a[U-j]    ; }
     for (i=0  ; i<KL0   ; i++) for (j= 0; j<=BL1; j++) { in[0]=i; in[1]=  j; A[in]=(*CL)[in] ; }
     for (i=0  ; i<KR0   ; i++) for (j= 0; j<=BR1; j++) { im[0]=i; im[1]=j; in[0]=M-KR0+1+i; in[1]=N-BR1+j; A[in]=(*CR)[im] ; }
     //A.Print() ; exit(-1)  ;
     for (r=0; r<IF->nr; r++) A.AdaptiveAPlusBMatVec (1.0,q,F,T, IF->a[r], IF->e[r] ) ;
     goto ende ;
   }
   
   KR0=M-KR0 ; KR1=N-KR1 ;
   BR0=M-BR0 ; BR1=N-BR1 ;

   // treat upper left and lower right block separately, since they may have some overlap
   for (r=0; r<IF->nr; r++) {
     CL->AdaptiveAPlusBMatVec (1.0,q,F,T, IF->a[r], min(IF->e[r],BL1) ) ;
     if (IF->e[r]>=BL1) break ;
   }    

   for (r=IF->nr-1; r>=0; r--) {
     CR->AdaptiveAPlusBMatVec (1.0,q,&F[BR1],&T[KR0+1],max(IF->a[r]-BR1,0)  , min(IF->e[r],N)-BR1) ;
     if (IF->a[r]<=BR1) break ;
   }

   r=0 ;
   while (r<IF->nr) {
     CLT->AdaptiveAPlusBMatVec(1.0,q,F,&T[KL0],IF->a[r],min(IF->e[r],KL1-1)) ;
     e=min(IF->e[r],BL1);
     for (i=max(KL1,IF->a[r]);i<=e;i++) {
       x=F[i]*q ;
       for (j=max(KL0,i-U);j<=min(KR0,i+O);j++) T[j] +=a[j-i+U]*x ;
     }
     if(e==BL1) goto innen ; // ATTENTION: i is either F->a[r], if e< BL1     OR     i==BL1+1 if e==BL1
     r++ ;
    }
 } else { // periodic BC
   r=0 ;
   while (r<IF->nr) {
     e=min(IF->e[r],BL1) ;
     for (i=IF->a[r];i<=e;i++) {
       x=F[i]*q ; je=i+O+l24; b=&a[-i+U-l24] ;
       for (j=i-U+l24;j<=je;j++) T[j&l21] +=b[j]*x ;
     }
     if(e==BL1) goto innen ;
     r++ ;
   }
 }

//
// Inneres 
 i=-MAXINT ; // should never be executed (if  r < IF->nr)
innen: ;

 // optimized inner loops
 if ((U==4) && (O==4)) {
   if (fabs(a[0]+a[8])<1e-15) {OperatorMatrixInnerLoop44Odd (&i,&r,F,T,IF,BR1-1,a,q);goto rechts;}
   if (fabs(a[0]-a[8])<1e-15) {OperatorMatrixInnerLoop44Even(&i,&r,F,T,IF,BR1-1,a,q);goto rechts;}
  }
 if ((U==1) && (O==2)) {
   if (fabs(a[0]+a[3])<1e-15) {OperatorMatrixInnerLoop12Odd (&i,&r,F,T,IF,BR1-1,a,q);goto rechts;}
   if (fabs(a[0]-a[3])<1e-15) {OperatorMatrixInnerLoop12Even(&i,&r,F,T,IF,BR1-1,a,q);goto rechts;}
  }

 if ((U==7) && (O==6)) {
   if (fabs(a[0]+a[13])<1e-15) {OperatorMatrixInnerLoop76Odd (&i,&r,F,T,IF,BR1-1,a,q);goto rechts;}
   if (fabs(a[0]-a[13])<1e-15) {OperatorMatrixInnerLoop76Even(&i,&r,F,T,IF,BR1-1,a,q);goto rechts;}
  }

 // ordinary inner loop
 while (r<IF->nr) {
   e=min(IF->e[r],BR1-1) ;
   while (i<=e) {
     x=F[i]*q ; b=&a[-i+U] ; je=i+O ;
     for (j=i-U;j<=je;j++) T[j] +=b[j]*x ; // Bem. leider sind die 'gemischten' OperatorMatrizen <dxPsi,Phi> nicht (anti-)symm. !!
     i++ ; 
   } 
   if (e==BR1-1) goto rechts ;
   r++; 
   i=IF->a[r] ;
 }

// right boundary
rechts: ;
 if (BCT[1]!=PERIODIC) {
   while (r < IF->nr) {
     CRT->AdaptiveAPlusBMatVec(1.0,q,&F[KR1+1],&T[BR0],max(IF->a[r]-KR1-1,0),min(IF->e[r],N)-KR1-1) ;
     for (i=max(i,IF->a[r]);i<=min(IF->e[r],KR1);i++) {
       x=F[i]*q ;
       for (j=max(KL0,i-U);j<=min(KR0,i+O);j++) T[j] +=a[j-i+U]*x ;
      }
     r++ ;
   }
 } else {  // periodic BC
   while (r < IF->nr) {
     for (i=max(i,IF->a[r]);i<=min(IF->e[r],l2-1);i++) {
       x=F[i]*q ; je=i+O+l24; b=&a[-i+U-l24] ;
       for (j=i-U+l24;j<=je;j++) T[j&l21] +=b[j]*x ;
      }
     r++ ;
   }
 }

ende:;
 // copy result 
 IT->VecCopy(T1,T) ;
 buf->unlock() ;
}

//                                                          //
//  Fast versions of inner loop for OperatorMatrix::MatVec  //
//                                                          //
//     for Interpolets(6) only up to now                    //
//                                                          //

void OperatorMatrixInnerLoop44Odd(int *ii  , int *rr   , double *F, double *T , IndexSet *IF , 
                                  int BR11 , double *a , double q) {
 int r=*rr,i=*ii,e ;
 double x,y/*,*Ti*/ ;

 while (r<IF->nr) {
   e=min(IF->e[r],BR11) ;
   while (i<=e) {
     x=F[i]*q ;
     y=a[0]*x ; T[i-4] +=y ; T[i+4] -=y ;
     y=a[1]*x ; T[i-3] +=y ; T[i+3] -=y ;
     y=a[2]*x ; T[i-2] +=y ; T[i+2] -=y ;
     y=a[3]*x ; T[i-1] +=y ; T[i+1] -=y ;
     /*
     Ti=&T[i] ;
     y=a[0]*x ; Ti[-4] +=y; Ti[4] -=y ;
     y=a[1]*x ; Ti[-3] +=y; Ti[3] -=y ;
     y=a[2]*x ; Ti[-2] +=y; Ti[2] -=y ;
     y=a[3]*x ; Ti[-1] +=y; Ti[1] -=y ;
     */
     i++ ;
   }
   if (e==BR11) { *ii=i , *rr=r ; return ; }
   r++; 
   i=IF->a[r] ;
 }

*ii=i , *rr=r ;
}


void OperatorMatrixInnerLoop44Even(int *ii  , int *rr   , double *F, double *T , IndexSet *IF , 
                                   int BR11 , double *a , double q) {
 int r=*rr,i=*ii,e ;
 double x,y ;
 
 while (r<IF->nr) {
   e=min(IF->e[r],BR11) ;
   while (i<=e) {
     x=F[i]*q ;
     y=a[0]*x ; T[i-4] +=y ; T[i+4] +=y ;
     y=a[1]*x ; T[i-3] +=y ; T[i+3] +=y ;
     y=a[2]*x ; T[i-2] +=y ; T[i+2] +=y ;
     y=a[3]*x ; T[i-1] +=y ; T[i+1] +=y ;
     T[i] +=a[4]*x ;
     i++ ;
   }
   if (e==BR11) { *ii=i , *rr=r ; return ; }
   r++; 
   i=IF->a[r] ;
 }

*ii=i , *rr=r ;
}


void OperatorMatrixInnerLoop12Odd(int *ii  , int *rr   , double *F, double *T , IndexSet *IF , 
                                  int BR11 , double *a , double q) {
 int r=*rr,i=*ii,e ;
 double x,y ;
 
 while (r<IF->nr) {
   e=min(IF->e[r],BR11) ;
   while (i<=e) {
     x=F[i]*q ;
     y=a[0]*x ; T[i-1] +=y ; T[i+2] -=y ;
     y=a[1]*x ; T[i  ] +=y ; T[i+1] -=y ;
     i++ ;
   }
   if (e==BR11) { *ii=i , *rr=r ; return ; }
   r++; 
   i=IF->a[r] ;
 }

*ii=i , *rr=r ;
}


void OperatorMatrixInnerLoop12Even(int *ii  , int *rr   , double *F, double *T , IndexSet *IF , 
                                   int BR11 , double *a , double q) {
 int r=*rr,i=*ii,e ;
 double x,y ;

 while (r<IF->nr) {
   e=min(IF->e[r],BR11) ;
   while (i<=e) {
     x=F[i]*q ;
     y=a[0]*x ; T[i-1] +=y ; T[i+2] +=y ;
     y=a[1]*x ; T[i  ] +=y ; T[i+1] +=y ;
     i++ ;
   }
   if (e==BR11) { *ii=i , *rr=r ; return ; }
   r++; 
   i=IF->a[r] ;
 }

*ii=i , *rr=r ;
}


void OperatorMatrixInnerLoop76Odd(int *ii  , int *rr   , double *F, double *T , IndexSet *IF , 
                                  int BR11 , double *a , double q) {
 int r=*rr,i=*ii,e ;
 double x,y ;
 
 while (r<IF->nr) {
   e=min(IF->e[r],BR11) ;
   while (i<=e) {
     x=F[i]*q ;
     y=a[0]*x ; T[i-7] +=y ; T[i+6] -=y ;
     y=a[1]*x ; T[i-6] +=y ; T[i+5] -=y ;
     y=a[2]*x ; T[i-5] +=y ; T[i+4] -=y ;
     y=a[3]*x ; T[i-4] +=y ; T[i+3] -=y ;
     y=a[4]*x ; T[i-3] +=y ; T[i+2] -=y ;
     y=a[5]*x ; T[i-2] +=y ; T[i+1] -=y ;
     y=a[6]*x ; T[i-1] +=y ; T[i  ] -=y ;
     i++ ;
   }
   if (e==BR11) { *ii=i , *rr=r ; return ; }
   r++; 
   i=IF->a[r] ;
 }

*ii=i , *rr=r ;
}

void OperatorMatrixInnerLoop76Even(int *ii  , int *rr   , double *F, double *T , IndexSet *IF , 
                                  int BR11 , double *a , double q) {
 int r=*rr,i=*ii,e ;
 double x,y ;
 
 while (r<IF->nr) {
   e=min(IF->e[r],BR11) ;
   while (i<=e) {
     x=F[i]*q ;
     y=a[0]*x ; T[i-7] +=y ; T[i+6] +=y ;
     y=a[1]*x ; T[i-6] +=y ; T[i+5] +=y ;
     y=a[2]*x ; T[i-5] +=y ; T[i+4] +=y ;
     y=a[3]*x ; T[i-4] +=y ; T[i+3] +=y ;
     y=a[4]*x ; T[i-3] +=y ; T[i+2] +=y ;
     y=a[5]*x ; T[i-2] +=y ; T[i+1] +=y ;
     y=a[6]*x ; T[i-1] +=y ; T[i  ] +=y ;
     i++ ;
   }
   if (e==BR11) { *ii=i , *rr=r ; return ; }
   r++; 
   i=IF->a[r] ;
 }

*ii=i , *rr=r ;
}



// alternative matrix-vector multiply, especially for 'completetly consistent' finite difference schemes
// There is no access on interpolated values.

void OperatorMatrix::MatVec2(double *T1,IndexSet *IT,int *BCT  ,  double *F,IndexSet *IF,int *BCF,
                             int l,int , double DX, doubleBuffer *)
{
 if (IT->nr==0) return ;

 Matrix2D *CL,*CR ;
 int       M,N ;
 int       K0,K1 ;
 int       i,j,r,er,ar ; 
 double    q=pow((1<<l)/DX , power),s ;

 CL=CR=NULL ;

 IT->VecClear(T1) ;

 M=N= 1<<l ;

 if (IsWT) M-- ;
 if (IsWF) N-- ;

 K0=AL[0][0].size[0] ;
 K1=AL[0][0].size[1] ;

 if (BCF[1]!=PERIODIC) {
   assert(BCT[1]!=PERIODIC) ;
   CL =&AL [0][0] ;
   CR =&AR [0][0] ;
 } 

 for (r=0; r<IF->nr; r++) {
   ar=IF->a[r]; er=IF->e[r]; 
   // check size of interval
   assert ((er-ar+1) >= K1) ;

   // left boundary
   CL->APlusBMatVec (0.0,q, &F[ar], &T1[ar] ) ;

   // right boundary
   CR->APlusBMatVec (0.0,q, &F[er-K1+1], &T1[er-K0+1] ) ;

   ar += K0 ;
   er -= K0 ;

   for (i=ar; i<=er ; i++) {
     s=0 ;
     for (j=-U; j<=O; j++) s += a[U+j]*F[i-j] ;
     T1[i]=s*q ;
   }
 }

}




//
//  TransformMatrix 
//

void TransformMatrix::MatVec (double *T1,IndexSet *IT,double *F,IndexSet *IF,int *BC,int l,int clflag, doubleBuffer *buf)
{
 if (IT->nr==0) return ;

 Matrix2D *HGL,*HGR       ;
 int      BL0,BL1,BR0,BR1 ;
 int      l2,M,N,l21,l24  ;
 int      r,i,j,e,je      ; 
 double   x,*b,*T=buf->lock(); 
 assert(T);

 HGL=HGR=NULL ;

 // get dimensions H/GJ is a (0,..,M)x(0,..,N) matrix
 N=1<<l ;
 l2=M=1<<(l-1) ; if(IsWT) M-- ;
 l21=l2-1 ;
 l24=l2*4 ; 

 if (BC[1]!=PERIODIC) {i=BC[0]+1 ; j=BC[1]+1 ;}
                       else {i=j=0 ;}

 BL0=AL[i].size[0]-1 ;             // Groessen der Randbloecke 
 BL1=AL[i].size[1]-1 ;

 BR0=AR[j].size[0]-1 ; BR0=M-BR0 ;
 BR1=AR[j].size[1]-1 ; BR1=N-BR1 ;

 // dirty hack for hat functions
 if (BL1==-1) { BL1=0 ; BR1=N ; }

 // clear T
 if (clflag) IT->VecClear(T)   ;
       else  IT->VecCopy(T,T1) ;

 // left boundary
 r=0 ; 
 if (BC[1]!=PERIODIC) {
   HGL=&AL[BC[0]+1] ; 
   HGR=&AR[BC[1]+1] ;
   while (r<IF->nr) {
     HGL->AdaptiveAPlusBMatVec(1.0,1.0,F,T,IF->a[r],min(IF->e[r],BL1)) ;
     e=min(IF->e[r],BL1);
     for (i=IF->a[r];i<=e;i++) {
       x=F[i] ;
       for (j=max((i-O+1)/2,BL0+1);j<=min((i+U)/2,BR0-1);j++) T[j] +=a[i+U-2*j]*x ;
     }
     if(e==BL1) goto innen ; // Achtung: i ist entweder F->a[r], falls e< BL1     ODER     i==BL1+1 falls e==BL1
     r++ ;
    }
 } else { // pericodic BC
   while (r<IF->nr) {
     e=min(IF->e[r],BL1) ;
     for (i=IF->a[r];i<=e;i++) {
       x=F[i] ; je=(i+U)/2+l24 ; b=&a[i+U+2*l24] ;
       for (j=ceil2(i-O)+l24; j<=je; j++) T[j&l21] +=b[-2*j]*x ;
     }
     if(e==BL1) goto innen ;
     r++ ;
   }
 }

//
// Inneres 
 i=-MAXINT ; // duerfte nie durchlaufen werden, wenn r < IF->nr
innen: ;

 if ((U==0) && (O==0) && (fabs(a[0]-1.0)<1e-14)) { 
   TransformMatrixInnerLoop00(&i,&r, F,T , IF , BR1-1 , a) ; 
   goto rechts ;
  }

 while (r<IF->nr) {
   e=min(IF->e[r],BR1-1) ;
   while (i<=e) {
     x=F[i] ; b=&a[i+U] ; je=(i+U)/2 ;
     for (j=(i-O+1)/2;j<=je;j++) T[j] +=b[-2*j]*x ;
     i++ ; 
   } 
   if (e==BR1-1) goto rechts ;
   r++; 
   i=IF->a[r] ;
 }


// right boundary 

rechts: ;
 if (BC[1]!=PERIODIC) {
   while (r < IF->nr) {
     HGR->AdaptiveAPlusBMatVec (1.0,1.0,&F[BR1],&T[BR0],max(IF->a[r]-BR1,0),min(IF->e[r],N)-BR1) ;
     for (i=max(i,IF->a[r]);i<=IF->e[r];i++) {
       x=F[i] ;
       for (j=max((i-O+1)/2,BL0+1);j<=min((i+U)/2,BR0-1);j++) T[j] +=a[i+U-2*j]*x ;
      }
     r++ ;
   }
 } else {  // periodic BC
   while (r < IF->nr) {
     for (i=max(i,IF->a[r]);i<=min(IF->e[r],2*l2-1);i++) {
       x=F[i] ; je=(i+U)/2+l24 ; b=&a[i+U+2*l24] ;
       for (j=ceil2(i-O)+l24; j<=je; j++) T[j&l21] +=b[-2*j]*x ;
      }
     r++ ;
   }
 }

 //copy: ; 
 IT->VecCopy(T1,T) ;
 buf->unlock() ;
}

//
// Fast Inner Loops
//
void TransformMatrixInnerLoop00(int *ii  , int *rr   , double *F, double *T , IndexSet *IF , 
                                int BR11 , double *) {
 int r=*rr,i=*ii,e,ia ;
  
 while (r<IF->nr) {
   e=min(IF->e[r],BR11) ;
   ia=i ;
   i =(i+1) & 0xfffffe  ; // round (up) to even
   while (i<=e) {
     T[i>>1] += F[i] ;
     i +=2 ;
   }
   if (e==BR11) { *ii=max(e+1,ia) , *rr=r ; return ; }
   r++; 
   i=IF->a[r] ;
 }

*ii=i , *rr=r ;
}



// multiply by transpose matrix
void TransformMatrix::MatTVec (double *T1,IndexSet *IT,double *F,IndexSet *IF,int *BC,int l,int clflag, doubleBuffer *buf)
{
 if (IT->nr==0) return ;

 Matrix2D *HGL=NULL,*HGR=NULL ;
 int       l2,M,N,l21,l24 ;
 int       BL1,BR0,BR1 ;
 int       i,j,e,r,je ;
 double    x,*b, *T=buf->lock(); assert(T);

 // get dimensions
 M =1<<(l+1) ;
 l2=N=1<< l  ; if(IsWT) N-- ;
 l21=M-1; l24=M*4 ;

 if (BC[1]!=PERIODIC) {i=BC[0]+1; j=BC[1]+1;}
                       else {i=j=0;}

 //BL0=AL[0].size[1]-1 ;             // Groessen der Randbloecke 
 BL1=AL[i].size[0]-1 ;
 BR0=AR[j].size[1]-1 ; BR0=M-BR0 ;
 BR1=AR[j].size[0]-1 ; BR1=N-BR1 ;

 // clear T
 if (clflag) IT->VecClear(T)   ;
       else  IT->VecCopy(T,T1) ;

 //left 
 r=0 ; 
 if (BC[1]!=PERIODIC) {
   HGL=&AL[BC[0]+1] ; 
   HGR=&AR[BC[1]+1] ;
   while (r<IF->nr) {
     HGL->AdaptiveAPlusBMatTVec(1.0,1.0,F,T,IF->a[r],min(IF->e[r],BL1)) ;
     e=min(IF->e[r],BL1);
     for (i=max(IF->a[r],BL1+1);i<=e;i++) {
       x=F[i] ;
       for (j=2*i-U;j<=2*i+O;j++) T[j] +=a[j-2*i+U]*x ;
     }
     if(e==BL1) goto innen ; // Achtung: i ist entweder F->a[r], falls e< BL1     ODER     i==BL1+1 falls e==BL1
     r++ ;
    }
 } else { // periodic BC
   while (r<IF->nr) {
     e=min(IF->e[r],BL1) ;
     for (i=IF->a[r];i<=e;i++) {
       x=F[i] ; je=2*i+O+l24 ; b=&a[-2*i+U -l24] ;
       for (j=2*i-U+l24; j<=je; j++) T[j&l21] +=b[j]*x ; 
     }
     if(e==BL1) goto innen ;
     r++ ;
   }
 }


//
// Inneres 
 i=-MAXINT ; // duerfte nie durchlaufen werden, wenn r < IF->nr
innen: ;
 if ((U==5) && (O==5) && (fabs(a[0]-a[10])<1e-14)) { 
   TransformMatrixTransposeInnerLoop55(&i,&r, F,T , IF , BR1-1 , a) ;
   goto rechts ;
  }
 if ((U==-1) && (O==1) && (fabs(a[0]-1.0)<1e-14)) { 
   TransformMatrixTransposeInnerLoopM11(&i,&r, F,T , IF , BR1-1 , a) ;
   goto rechts ;
  }


 while (r<IF->nr) {
   e=min(IF->e[r],BR1-1) ;
   while (i<=e) {
     x=F[i] ; b=&a[-2*i+U] ; je=2*i+O ;
     for (j=2*i-U;j<=je;j++) T[j] +=b[j]*x ;
     i++ ; 
   } 
   if (e==BR1-1) goto rechts ;
   r++; 
   i=IF->a[r] ;
 }

// right
rechts: ;
 if (BC[1]!=PERIODIC) {
   while (r < IF->nr) {
     HGR->AdaptiveAPlusBMatTVec (1.0,1.0,&F[BR1],&T[BR0],max(IF->a[r]-BR1,0),min(IF->e[r],N)-BR1) ;
     for (i=max(i,IF->a[r]);i<=min(IF->e[r],BR1-1);i++) {
       x=F[i] ;
       for (j=2*i-U;j<=2*i+O;j++) T[j] +=a[j-2*i+U]*x ;
      }
     r++ ;
   }
 } else {  // periodic BC
   while (r < IF->nr) {
     for (i=max(i,IF->a[r]);i<=min(IF->e[r],l2-1);i++) {
       x=F[i] ; je=2*i+O+l24 ; b=&a[-2*i+U -l24] ;
       for (j=2*i-U+l24; j<=je; j++) T[j&l21] +=b[j]*x ;
      }
     r++ ;
   }
 }

 IT->VecCopy(T1,T) ;
 buf->unlock() ;
}

//
// Fast Inner Loops
//
void TransformMatrixTransposeInnerLoop55(int *ii  , int *rr   , double *F, double *T , IndexSet *IF , 
                                         int BR11 , double *a) {
 int r=*rr,i=*ii,e,i2 ;
 double x,y ;
 
 while (r<IF->nr) {
   e=min(IF->e[r],BR11) ;
   while (i<=e) {
     x=F[i] ;
     i2=2*i ;
     y=a[0]*x ; T[i2-5] +=y ; T[i2+5] +=y ;
     y=a[2]*x ; T[i2-3] +=y ; T[i2+3] +=y ;
     y=a[4]*x ; T[i2-1] +=y ; T[i2+1] +=y ;
     T[i2] +=x ; 
     i++ ;
   }
   if (e==BR11) { *ii=i , *rr=r ; return ; }
   r++; 
   i=IF->a[r] ;
 }

*ii=i , *rr=r ;
}

void TransformMatrixTransposeInnerLoopM11(int *ii  , int *rr   , double *F, double *T , IndexSet *IF , 
                                          int BR11 , double *) {
 int r=*rr,i=*ii,e ;
 double *T1=T+1 ;
  
 while (r<IF->nr) {
   e=min(IF->e[r],BR11) ;
   while (i<=e) {
     T1[2*i] +=F[i] ;
     i++ ;
   }
   if (e==BR11) { *ii=i , *rr=r ; return ; }
   r++; 
   i=IF->a[r] ;
 }

*ii=i , *rr=r ;
}


//
// non-adaptive Matrix-Vector multiplies
//

void OperatorMatrix::MatVec(double *T,int *BCT  ,  double *F,int *BCF,  int J, double DX) {
 int i,j, nJ=1<<J ,b,e ;

 double *au=&a[U] , q=pow( nJ/DX , power) , s ;

 if (BCF[1]==PERIODIC) {
 
   assert(BCT[1]==PERIODIC) ;

   int mask=nJ-1 ;
   for (i=nJ ;i<2*nJ ;i++) {
     s = 0 ; 
     for (j=-U; j<=O; j++) s+=au[j]*F[(i-j)&mask] ;
     T[i-nJ]=s*q ;
    }
   T[nJ]=0.0 ;

 } else {

   int KL0,KL1,KR0,KR1 , BL0,BL1,BR0,BR1 ;
   KL0=AL[0][0].size[0] ; KL1=ALT[0][0].size[1] ;
   KR0=AR[0][0].size[0] ; KR1=ART[0][0].size[1] ;

   BL1=AL[0][0].size[1]-1 ; BL0=KL0+ALT[0][0].size[0]-1 ;
   BR1=AR[0][0].size[1]-1 ; BR0=KR0+ART[0][0].size[0]-1 ;

   assert((BL0+0+KR0)<=nJ) ;
   assert((BL1+0+KR1)<=nJ) ;

   KR0=nJ-KR0 ; KR1=nJ-KR1 ;
   BR0=nJ-BR0 ; BR1=nJ-BR1 ;
 
   Matrix2D *ALB,*ALBT,*ARB,*ARBT ;
   ALB =&AL [BCT[0]+1][BCF[0]+1] ;
   ALBT=&ALT[BCT[0]+1][BCF[0]+1] ;
   ARB =&AR [BCT[1]+1][BCF[1]+1] ;
   ARBT=&ART[BCT[1]+1][BCF[1]+1] ;

   // Inner scaling functions
 
   for (i=KL0; i<=KR1; i++) {
     b=max(i-U , KL1) ;
     e=min(i+O , KR1) ; 
     s=0 ; 
     for (j=b; j<=e; j++) s+=au[i-j]*F[j] ;
     T[i]=s*q ;
    }

   // 
   // Multiply for boundaries

    ALB->APlusBMatVec(0.0,q,F,T) ;
   ALBT->APlusBMatVec(1.0,q,F,&T[KL0]) ;

    ARB->APlusBMatVec(0.0,q,&F[BR1  ],&T[KR0+1]) ;
   ARBT->APlusBMatVec(1.0,q,&F[KR1+1],&T[BR0  ]) ;

 }

}

void LiftingMatrix::MatTVec(double *T1,IndexSet *IT  , double *F,IndexSet *IF, int *BC, int J , int add, doubleBuffer *buf) {
 int i,j,r,e,p, N=1<<J , M=N-1 , b0=BC[0] , b1=BC[1] ;
 
 // nothing to do ?
 if (IT->nr==0) return ;
 if (IF->nr==0) {
   if (add==0) IT->VecClear(T1) ;
   return ;
 }

 double *T=buf->lock(); assert(T);
 double s, *b, q = (add>=0) ? 1 : -1 ;

 // clear T
 if (add==0) IT->VecClear(T)   ;
       else  IT->VecCopy(T,T1) ;
 
 if (b1==PERIODIC) {
   for (r=0; r<IF->nr; r++) {
     e=min(IF->e[r],M) ;
     for (i=IF->a[r]; i<=e; i++) {
       s = q*F[i] ;
       b = &a[U-N-i] ;
       p = i+O+N ;
       for (j=i-U+N; j<=p; j++) T[j&M] +=b[j]*s ;
     }
    }
   IT->VecCopy(T1,T) ;
   buf->unlock() ;
   return ;
 }

 //
 // not periodic

 int BWL=QL[b0+1].size[0]-1 ,
   //  BL =QL[b0+1].size[1]-1 ,
     BVL=BWL+PL[b0+1].size[0] ,
     BWR=M-QR[b1+1].size[0]+1 ,
     BR =N-QR[b1+1].size[1]+1 ,
     BVR=BWR-PR[b1+1].size[0] ;

 int n=IF->nr ;

 //
 // Boundary Multiplication, at most the first run can touch the boundary and if it does, it covers completely
 // the boundary block
 r=0 ;
 i=IF->a[0] ;

 // left boundary
 if (i<=BWL) {
   e=min(IF->e[0] , BWL) ;
   QL[b0+1].AdaptiveAPlusBMatTVec(1.0 , q ,  F , T , i,e) ;
   if (e==IF->e[0]) { r++ , i=IF->a[r] ; }
               else   i=e+1 ;
  }

 // inner
 while (r<n) {
   e=min(IF->e[r],BWR-1) ;
   while (i<=e) {
     b=&a[U-i] ;
     if (i<=BVL) b=&PL[b0+1].a[i-BWL-1][U-i] ;
     if (i>=BVR) b=&PR[b1+1].a[i-BVR  ][U-i] ;
     s =q*F[i] ;
     p =i+O ;
     for (j=i-U; j<=p; j++) T[j] +=b[j]*s ;
     i++ ;
   }
   r++ ;
   i=IF->a[r] ;
 }

 // right boundary
 i=IF->a[n-1];
 e=min(M,IF->e[n-1]) ;
 if (e>=BWR) QR[b1+1].AdaptiveAPlusBMatTVec(1.0 , q ,  &F[BWR] , &T[BR] , max(i-BWR,0) , e-BWR) ;
 
 IT->VecCopy(T1,T) ;
 buf->unlock() ;
}

void LiftingMatrix::MatTVec(double *T, double *F , int *BC , int J , int add) {
 int i,j,e, N=1<<J , M=N-1 , b0=BC[0] , b1=BC[1] ;

 double s,*b,q = (add>=0) ? 1 : -1 ;
   
 if (add==0) for (i=0; i<=N; i++) T[i]=0 ;

 if (b1==PERIODIC) {
   for (i=0; i<N; i++) {
     s = q*F[i] ;
     b = &a[U-N-i] ;
     e = i+O+N ;
     for (j=i-U+N; j<=e; j++) T[j&M] +=b[j]*s ;
   }
   T[N]=0 ;
   return ;
 }

 int BWL=QL[b0+1].size[0]-1 ,
   //     BL =QL[b0+1].size[1]-1 ,
     BVL=BWL+PL[b0+1].size[0] ,
     BWR=M-QR[b1+1].size[0]+1 ,
     BR =N-QR[b1+1].size[1]+1 ,
     BVR=BWR-PR[b1+1].size[0] ;

 //
 // Rand-Multiplikation
 QL[b0+1].APlusBMatTVec(1.0 , q ,  F    , T) ;
 QR[b1+1].APlusBMatTVec(1.0 , q , &F[BWR],&T[BR]) ;
 
 // innen
 for (i=BWL+1; i<=BWR-1; i++) {
   b=&a[U] ;
   if (i<=BVL) b=&PL[b0+1].a[i-BWL-1][U] ;
   if (i>=BVR) b=&PR[b1+1].a[i-BVR  ][U] ;
   s =q*F[i] ;
   e =i+O ;
   for (j=i-U; j<=e; j++) T[j] +=b[j-i]*s ;
 }
}

---