Department of Scientific Computing    Institute for Numerical Simulation    University of Bonn

Documentation

Download

Programming References

Bug Reports / Suggestions

FAQ

Authors

           

---

Adaptive1DIndexSet.cc

#include "Adaptive1D.hpp"

extern int debugRefine ;

/**********************/
/*                    */
/*     IndexSet::     */
/*                    */
/**********************/

IndexSet::IndexSet()          {a=e=NULL ;nr=0 ;}
IndexSet::IndexSet(int Size)  {a=e=NULL ;nr=0 ;Init(Size);}
IndexSet::~IndexSet()         {delete []a; delete []e;}

// Init. of IndexSet with number of intervals <= Size 
// to avoid UMR errors the arrays are allocated for Size+1 intervals and are explicitely int.
int IndexSet::Init(int Size) {
 Size=min(Size,_MAX_INDEXSET_) ;
 if (a==NULL) {
   a=new int[Size+1] ;
   e=new int[Size+1] ;
   for (int i=0; i<=Size; i++) a[i]=e[i]=0 ;
 }
 nr=0 ;
 return 1 ;
}

int IndexSet::Print(char *st) 
{ int c=0 ;
 if (st) std::cout <<st ;
 std::cout << "nr: "<<nr<<' '<<a<<','<<e<<'\n';
 for (int i=0; i<nr; i++) {  
   if (st) std::cout <<st ;
   std::cout << "r "<<i<<": "<<a[i]<<' '<<e[i]<<'\n';
   c+=(e[i]-a[i]+1) ; 
 }
 return c ;
}

// print three columns: *st, i*2^(-*l*) , l     //  *l* = (l==Level0) ? l : l-1 
void IndexSet::Print(char *st,int l,int Level0) 
{ int si = (l==Level0) ? 1<<l : 1<<(l-1) ;
 for (int r=0; r<nr; r++)
   for (int i=a[r];i<=e[r];i++)
     std::cout<<st<<' '<<((double)i)/si<<' '<<l<<'\n' ;
}

size_t IndexSet::Size() {
  size_t s=0 ;
  for (int r=0; r<nr; r++) s += (size_t)(e[r]-a[r]+1) ; 
  return s ;
}

void IndexSet::Copy(IndexSet *F)
{
 if (this==F) return ;
 nr=F->nr ; assert(nr<_MAX_INDEXSET_) ;
 for (int i=0;i<nr;i++) a[i]=F->a[i] , e[i]=F->e[i] ;
}

//
// essentially the given index set of scaling functions (!) is copied, 
// but if scaling functions on the interval are given, then the following condition
// is maintained: if any boundary scaling function is active then all boundary scaling
// functions are activated
void IndexSet::CopyWithBounds(IndexSet *B,int l, int *BC, Wavelets *W) {
  int *Ba=B->a , *Be=B->e ,Bn=B->nr , r, iend=(1<<l) ;
  int K0=-1,NKW1=iend+1 ;

  assert(Bn<_MAX_INDEXSET_) ;

  if (BC[1]!=PERIODIC) { K0=W->K0 , NKW1=iend-W->K1 ;}
  
  nr=0 ;
  if (Bn==0) return ;

  r=0 ;
  // left boundary
  if (Ba[0]<K0) {
    a[0]=0 ; e[0]=K0-1 ;
    while ((r<Bn)&&(Ba[r]<K0)) { e[0]=max(e[0],Be[r]) ;r++;}
    nr=1 ;
  }
  // center
  for (;r<Bn;r++) { a[nr]=Ba[r] ; e[nr++]=Be[r] ; }
  
  // right boundary
  r=nr-1 ;
  while ((r>=0) && (e[r]>NKW1)) r-- ;

  if (r+1<nr) {
    a[r+1]=min(a[r+1],NKW1) ;
    e[r+1]=iend;
    nr=r+2 ;
  }
  
 // Patch: if the last interval and it's predecessor touch merge them
 if (nr>1)
   if (a[nr-1]==e[nr-2]+1) { e[nr-2]=e[nr-1] ; nr--;}
}

// scalings(l) which span the given wavelets(l)
void IndexSet::InducedByWavelets(IndexSet *B,int l, int *BC, Wavelets *W)
{                                       
  if(l==W->Level0) { CopyWithBounds(B,W->Level0,BC,W) ;return ;} // Auf level0 ist nichts zu tun 
  //if(l==Level0) { this->Copy(B) ; return ;}

  struct MatrixShape H ;
  W->GtJ.GetShape(l,&H) ;
  H.Transpose()        ;
  CAM(B,BC,&H,TransformMatrixTransposeFirstNZEP , TransformMatrixTransposeLastNZEP ,
              TransformMatrixTransposeFirstNZE  , TransformMatrixTransposeLastNZE) ;

  assert(nr<_MAX_INDEXSET_) ;
} 

// scalings(l) influenced by the given scalings(l-1)
void IndexSet::InducedByCoarseScalings(IndexSet *B,int l, int *BC, Wavelets *W)
{  
  struct MatrixShape H ;
  W->HJ.GetShape(l,&H) ;
  H.Transpose()        ;
  CAM(B,BC,&H,TransformMatrixTransposeFirstNZEP , TransformMatrixTransposeLastNZEP ,
              TransformMatrixTransposeFirstNZE  , TransformMatrixTransposeLastNZE) ;

  assert(nr<_MAX_INDEXSET_) ;
} 

// scalings(l-1) required to span the given scalings(l)
void IndexSet::InducedByFineScalings(IndexSet *S,int l, int *BC, Wavelets *W) { 
  struct MatrixShape H ;
  W->HJ.GetShape(l,&H) ;
  CAM(S,BC,&H,TransformMatrixFirstNZEP , TransformMatrixLastNZEP ,
              TransformMatrixFirstNZE  , TransformMatrixLastNZE) ;
}

// dual scalings(l-1) required to span the given dual scalings(l)
void IndexSet::InducedByFineDualScalings(IndexSet *S,int l, int *BC, Wavelets *W) { 
  struct MatrixShape H ;
  W->HtJ.GetShape(l,&H) ;
  CAM(S,BC,&H,TransformMatrixFirstNZEP , TransformMatrixLastNZEP ,
              TransformMatrixFirstNZE  , TransformMatrixLastNZE) ;

  assert(nr<_MAX_INDEXSET_) ;
}

// wavelets(l) which coincide with the given scalings(l)
void IndexSet::WaveletsInducedByFineScalings(IndexSet *S,int l, int *BC, Wavelets *W) { 
  struct MatrixShape H ;
  //NEVER EVER:: if(l==Level0) { this->CopyWithBounds(S,Level0,BC,W) ;return ;} 
  if(l==W->Level0)  { Copy(S) ;return ;}
  W->GJ.GetShape(l,&H) ;
  CAM(S,BC,&H,TransformMatrixFirstNZEP , TransformMatrixLastNZEP ,
              TransformMatrixFirstNZE  , TransformMatrixLastNZE) ;

  assert(nr<_MAX_INDEXSET_) ;
}

// dual wavelets(l) influenced by the given dual scalings(l)
void IndexSet::WaveletsInducedByFineDualScalings(IndexSet *B,int l, int *BC, Wavelets *W) {
  struct MatrixShape H ;
  if(l==W->Level0)  { Copy(B) ;return ;}
  W->GtJ.GetShape(l,&H) ;
  //debugRefine=1 ;
  CAM(B,BC,&H,TransformMatrixFirstNZEP , TransformMatrixLastNZEP ,
              TransformMatrixFirstNZE  , TransformMatrixLastNZE) ;
  //debugRefine=0 ;
  //H.Print() ;
  //B->Print() ;
  //this->Print() ;

  assert(nr<_MAX_INDEXSET_) ;
}


void IndexSet::InducedByOperatorMatrix(IndexSet *S,int l, int *BC, OperatorMatrix *W) { 
  struct MatrixShape H ;
  W->GetShape(l,&H) ;
  CAM(S,BC,&H,OperatorMatrixFirstNZEP , OperatorMatrixLastNZEP ,
              OperatorMatrixFirstNZE  , OperatorMatrixLastNZE) ;

  assert(nr<_MAX_INDEXSET_) ;
}

// scalings(l-1) required to span the given lifting wavelets(l)
void IndexSet::ScalingsInducedByLifting(IndexSet *S,int l, int *BC, Wavelets *W) { 
  struct MatrixShape H  ;
  W->Q.GetShape(l-1,&H) ;
  CAM(S,BC,&H,LiftingMatrixTransposeFirstNZEP , LiftingMatrixTransposeLastNZEP ,
              LiftingMatrixTransposeFirstNZE  , LiftingMatrixTransposeLastNZE) ;

  assert(nr<_MAX_INDEXSET_) ;
}


void IndexSet::IndexSetForTypeIII(IndexSet *S,int l, int *BC, Wavelets *W) {
  struct MatrixShape H ;
  W->HtJ.GetShape(l,&H) ;
  EAM(S,BC,&H,TransformMatrixFirst2NZEP , TransformMatrixLast2NZEP ,
              TransformMatrixFirst2NZE  , TransformMatrixLast2NZE) ;

  assert(nr<_MAX_INDEXSET_) ;
}

//
// calculates the intersection of I with the set of boundary wavelets
void IndexSet::ForBoundaryWavelets(bool left, bool right, IndexSet *I, int level,Wavelets *W) {
  int *Ia=I->a , *Ie=I->e , n=I->nr ,r ,KW0=W->KW0 , NKW1=(1<<(level-1))-1-W->KW1 ;
  bool flag ;

  nr=0 ;
  if (left) 
    for (r=0; r<n ;r++) if (Ia[r]<KW0) { a[nr]=Ia[r] ; e[nr++]=min(Ie[r],KW0-1) ; }
                                  else r=n ;

  if (right) {
    r=n-1 ;
    flag=true ;
    while ((r>=0)&&flag) { 
      if   (Ie[r]>NKW1) r-- ;
                   else flag=false ;
     }

    for (r=r+1; r<n; r++) { a[nr]=max(Ia[r],NKW1+1) ; e[nr++]=Ie[r] ;}
  }

  assert(nr<_MAX_INDEXSET_) ;
}

//
//extern int debugRefine ;
void IndexSet::CAM(IndexSet *S, int *BC, 
                   MatrixShape *H,
                   int (*FirstNZEP)(int,MatrixShape *) , int (*LastNZEP )(int,MatrixShape *) ,
                   int (*FirstNZE )(int,MatrixShape *) , int (*LastNZE  )(int,MatrixShape *)) 
{                         
 int *Sa=S->a , *Se=S->e , Sn=S->nr ;
 int aa,ee,a2,e2,r ;
 
 nr=0 ;
 if (Sn==0) return ; 

 if (BC[1]==PERIODIC) { 
   int M=H->M  /* N=H->N , U=H->U , O=H->O*/  ;
   M=(M+1)&0xfffe ;                            // M shall be even
   int a3=M ;

   aa=FirstNZEP(Sa[0],H) ,  ee=min(LastNZEP(Se[0],H) , M-1) ; 

   if (aa<0) { a3=aa+M , aa=0 ; }              // Ueberlapp ans Ende
   
   e2=LastNZEP(Se[Sn-1],H) ;
   if (e2>=M) {                                // Ueberlapp an den Anfang
     e2-=M ;
     if (e2+1 < aa) { a[nr]=0 , e[nr++]=e2 ; } // aber kein Ueberlapp mit erstem Intervall
     else   aa=0 ;                             // Ueberlapp -> erstes Intervall faengt mit 0 an
    }
  
   for (r=1; r<Sn; r++) {
     a2=max(FirstNZEP(Sa[r],H),0)  ,  e2=min(LastNZEP(Se[r],H) , M-1) ; // abgeschnittenes Intervall holen
     if (ee+1 < a2) { a[nr]=aa , e[nr++]=ee , aa=a2 ;}                           // kein Ueberlapp mit angefangenem
     ee=e2 ;                                                              
   }

   if (ee+1<a3) {                                                                // angefangenes hat kein Ueberlapp mit 
     a[nr]=aa , e[nr++]=ee ;                                                     // moeglichem letztem a3 Intervall
     if (a3<M) { a[nr]=a3 , e[nr++]=M-1 ; }                                      // letztes Intervall eintragen
   } else {                                                                      
     a[nr]=aa , e[nr++]=M-1 ;                                                    // Ueberlapp mit a3 Intervall
   }
   return ;
 }

 // not periodic
 aa=FirstNZE(Sa[0],H) ; // first prospective interval
 ee=LastNZE (Se[0],H) ;

 for (r=1; r< Sn ;r++) {
   a2=FirstNZE(Sa[r],H) ;
   e2=LastNZE (Se[r],H) ;

   if (a2<=ee+1) ee=e2 ;
   else a[nr]=aa , e[nr]=ee , nr++ , aa=a2 , ee=e2 ;
 }
 if(aa<=ee) a[nr]=aa , e[nr]=ee , nr++ ;
}

//
//
void IndexSet::EAM(IndexSet *S, int *BC, 
                   MatrixShape *H,
                   int (*FirstNZEP)(int,MatrixShape *) , int (*LastNZEP )(int,MatrixShape *) ,
                   int (*FirstNZE )(int,MatrixShape *) , int (*LastNZE  )(int,MatrixShape *)) {

 int *Sa=S->a , *Se=S->e , Sn=S->nr ;
 int  a2,e2,a3=-MAXINT,e3=-MAXINT,r=0 ;
 bool flag=false ; 

 int M=H->M,N=H->N ;
 M=(M+1)&0xfffe ; // M shall be a power of 2
 N=(N+1)&0xfffe ;

 nr=0 ;
 if (Sn==0) return ; 

 if (BC[1]==PERIODIC) { 

   // first interval may be a periodic interval together with the last one
   // -> special treatment necessary
   if ((Sa[0]==0) && (Se[Sn-1]>=N-1)) {
     e2=FirstNZEP(Se[0   ],H)   ;
     a2=LastNZEP (Sa[Sn-1],H)-M ;

     if (a2<=e2) { // valid interval
       if ((a2<0) && (e2< 0)) { a3=a2+M, e3=e2+M, flag=true; } else 
       if ((a2<0) && (e2>=0)) { a3=a2+M, e3=M-1 , flag=true; a[0]=0, e[nr++]=e2; } else
       if  (a2>=0)            { a[0]=a2, e[nr++]=e2 ; }
     }
     r=1  ;
     Sn-- ;
   }
   
   // general case
   for (;r<Sn;r++) {
     a2=LastNZEP (Sa[r],H) ;
     e2=FirstNZEP(Se[r],H) ;
     if (a2<=e2) { a[nr]=a2, e[nr++]=e2 ;}
   }

   if (flag) { a[nr]=a3, e[nr++]=e3 ;}
   return ;
  }

 // not periodic:
 for (r=0; r<Sn; r++) {
   a2=LastNZE (Sa[r],H) ;
   e2=FirstNZE(Se[r],H) ;
   if (a2<=e2) { a[nr]=a2, e[nr++]=e2 ;}
 }

 assert(nr<_MAX_INDEXSET_) ;
}


//
//
void IndexSet::Union(IndexSet *A,IndexSet *B)
{ 
  int i,a1,e1 ;
  int *Aa=A->a , *Ae=A->e , Ai ; // Berechnet This= A vereinigt B
  int *Ba=B->a , *Be=B->e , Bi ;
 
  assert(this!=A) ;
  assert(this!=B) ;

  Ai=Bi=i=0 ;
  while ((Bi<B->nr) && (Ai<A->nr)) {

// start new run
    if (Ba[Bi] < Aa[Ai]) a1=Ba[Bi] , e1=Be[Bi] , Bi++ ;
                    else a1=Aa[Ai] , e1=Ae[Ai] , Ai++ ;

    int cond;
    while ( (cond=(Ai<A->nr) && (e1>=Aa[Ai]-1)) || ((Bi<B->nr)&&(e1>=Ba[Bi]-1)) )
     {                                             // Fortsetzung moeglich 
       if (cond) {                                 // aktueller Run ist  Ba,Be[Bi-1] , aber Ueberlappung mit A run
         if (e1<Ae[Ai]) e1=Ae[Ai] ;                // aha, nur teilweise Uberlappung des A runs
         Ai++ ;                                    // auf nachsten A setzen
       } else {                                    // e1 < Aa[Ai] aber e1>=Ba[Bi] , aktueller run ist Aa,Ae[Ai-1] , aber Uberlapp. mit B run
         if(e1<Be[Bi]) e1=Be[Bi] ;                 // ansonsten Ende des Ier nachruecken
         Bi++ ;                                    // auf naechsten B run setzen
        } 
     } // while 
   
   // aha entweder A/B zu Ende oder erstmal keine Ueberlappung mehr
    
    a[i]=a1 , e[i]=e1 , i++ ;                     // A U B schreiben
  } // while A und B noch nicht zu Ende 

  if(Bi==B->nr) // keine B runs mehr also alle weiteren A nach This , D  kopieren 
    while (Ai<A->nr) a[i]=Aa[Ai] , e[i++]=Ae[Ai++] ;
  else
    while (Bi<B->nr) a[i]=Ba[Bi] , e[i++]=Be[Bi++] ;

  nr  =i  ;
  assert(nr<_MAX_INDEXSET_) ;
}
 

void IndexSet::Difference(IndexSet *A,IndexSet *B) // Berechnet A\B
{
 int Ai,Bi ;
 int *Aa=A->a , *Ae=A->e , An=A->nr ;
 int *Ba=B->a , *Be=B->e , Bn=B->nr ;

 assert(this!=A) ;
 assert(this!=B) ;

 Ai=Bi=nr=0 ;
 for (Ai=0; Ai<An; Ai++) {

   while((Bi<Bn)&&(Be[Bi]<Aa[Ai])) Bi++ ; // B Intervalle uberspringen, die vor akt. A liegt

   if (Bi>=Bn) { a[nr]=Aa[Ai] , e[nr++]=Ae[Ai] ; } // kein B Intervall mehr da -> A kopieren
   else {                                          // :: Be[Bi]>=Aa[Ai]

     if (Aa[Ai]<Ba[Bi]) { a[nr]=Aa[Ai] , e[nr++]=min(Ae[Ai],Ba[Bi]-1) ;} // Anfang von A kopieren

     while((Be[Bi]<Ae[Ai])&&(Bi<Bn)) { 
       a[nr  ]=Be[Bi]+1 ;
       e[nr++]=(Bi+1<Bn) ? min(Ae[Ai],Ba[Bi+1]-1) : Ae[Ai] ;
       Bi++ ; 
     }
   } 
 } // next Ai

 assert(nr<_MAX_INDEXSET_) ;
}

bool IndexSet::IsSame(IndexSet *I) {
 if (this==I) return true;
 int *aa=I->a , *ee=I->e , n=I->nr ;
 if (n!=nr) return false ;
 
 for (int r=0;r<n ;r++) {
   if (aa[r]!=a[r]) return false ;
   if (ee[r]!=e[r]) return false ;
 }
 return true ;
}


bool IndexSet::Contains(int i,int *rr) {
  int r=(rr) ? *rr : 0 ;
  if (nr==0) return false ;
  if (r>=nr)  r=0 ; // invalid guess
  if (a[r]>i) r=0 ;

  if (a[0]>i) return false ;
  // here: a[r]<=i 
  do {
    if (i<=e[r]) { 
      if (rr) *rr=r; 
      return true ;
     }

    r++ ;
  } while ((r<nr) && (a[r]<=i)) ;
   
  return false ;  
}

double IndexSet::VecMaxAbs (double *t) {
  double mx=-1 ;
  int r,i ;
  for (r=0; r<nr; r++)  
    for (i=a[r]; i<=e[r]; i++) mx = (fabs(t[i]) > mx) ? fabs(t[i]) : mx ;
  return mx ;
}

double IndexSet::VecMax (double *t) {
  double mx=-HUGE ;
  int r,i ;
  for (r=0; r<nr; r++)  
    for (i=a[r]; i<=e[r]; i++) mx = (t[i] > mx) ? t[i] : mx ;
  return mx ;
}

double IndexSet::VecInnerProd(double *v1,double *v2) {
 int r,i ;
 double s=0.0 ; 
 for (r=0; r<nr; r++)  
 for (i=a[r]; i<=e[r]; i++) s+=v1[i]*v2[i] ;
 return s ;
}

void IndexSet::VecAdd (double *t,double *p,double *q) {
 int r,i ;
 for (r=0; r<nr; r++)  
 for (i=a[r]; i<=e[r]; i++) t[i] =q[i]+p[i] ;
}

void IndexSet::VecSub (double *t,double *p,double *q) {
 int r,i ;
 for (r=0; r<nr; r++)  
 for (i=a[r]; i<=e[r]; i++) t[i] =p[i]-q[i] ;
}

void IndexSet::VecPlus(double *t,double *p) {
 int r,i ;
 for (r=0; r<nr; r++)  
 for (i=a[r]; i<=e[r]; i++) t[i] +=p[i] ;
}

void IndexSet::VecMinus(double *t,double p) {
 int r,i ;
 for (r=0; r<nr; r++)  
 for (i=a[r]; i<=e[r]; i++) t[i] -=p ;
}

void IndexSet::VecMinus(double *t,double *p) {
 int r,i ;
 for (r=0; r<nr; r++)  
 for (i=a[r]; i<=e[r]; i++) t[i] -=p[i] ;
}

/*
void IndexSet::VecCopy(double *t,double *p) {
 int r,i,er ;
 for (r=0; r<nr; r++) {
   er=e[r] ;
   for (i=a[r]; i<=er; i++) t[i]=p[i] ;
 }
}
*/

void IndexSet::VecCopy(double *t,double *p) {
 int r,i,er,n,n4,n2 ;
 
 for (r=0; r<nr; r++) {
   n=e[r]-a[r]+1 ;
   n4=n/4 ; er=n4*4 ; n=n-er ; er +=a[r] ;
   n2=n/2 ; n=n-(n2*2) ;

   for (i=a[r]; i<er; i+=4 ) {
     t[i  ]=p[i  ] ;
     t[i+1]=p[i+1] ;
     t[i+2]=p[i+2] ;
     t[i+3]=p[i+3] ;
   }

   if (n2>0) {
     t[i  ]=p[i  ];
     t[i+1]=p[i+1];
     i+=2 ;
   }
   if (n>0) t[i]=p[i] ;
 }
}

void IndexSet::VecMul(double *t, const double f) {
 int r,i ;
 for (r=0; r<nr; r++)  
 for (i=a[r]; i<=e[r]; i++) t[i] *= f ;
}

void IndexSet::VecFill(double *t,const double f) {
 int r,i ;
 for (r=0; r<nr; r++)  
 for (i=a[r]; i<=e[r]; i++) t[i]=f ;
}

void IndexSet::VecClear(double *t) { VecFill(t,0) ;}

void IndexSet::VecPrint(const char *st,double *t) {
 int r,i ;
 for (r=0; r<nr; r++)  
   for (i=a[r]; i<=e[r]; i++) std::cout<<st<<' '<<i<<' '<<t[i]<<'\n' ;
}


//
// Load Operation including Treatment of Ghostcells
//

void IndexSet::Load(AdaptiveLE *al) {
 int    ee,i ;
 AdaptiveLE::iterator alit,alend=(*al).end() ;

 nr=0 ;
 ee=-10000000 ;
 for (alit=(*al).begin();alit!=alend ;alit++) {
   i=(*alit).first  ;
   if (ee<i-1) { a[nr]=i ; ee=e[nr]=i ; nr++ ;} 
         else    ee=e[nr-1]=i ;
  }
} 

void IndexSet::Load(AdaptiveLE *al, vector<double> *v, double *d) {
 int    ee,i ;
 long   j    ;
 AdaptiveLE::iterator alit,alend=(*al).end() ;

 nr=0 ;
 ee=-10000000 ;
 for (alit=(*al).begin();alit!=alend ;alit++) {
   i=(*alit).first  ;
   j=(*alit).second ;
   if (ee<i-1) { a[nr]=i ; ee=e[nr]=i ; nr++ ;} 
         else    ee=e[nr-1]=i ;
   d[i] = (j<0) ? (*v)[1-j] : (*v)[j] ;
  }
}

void IndexSet::LoadIS(AdaptiveLE *al) {
 int    ee,i ;
 long   j    ;
 AdaptiveLE::iterator alit,alend=(*al).end() ;

 nr=0 ;
 ee=-MAXINT ;
 for (alit=(*al).begin();alit!=alend ;alit++) {
   i=(*alit).first  ;
   j=(*alit).second ;
   if (j>=0) {
     if (ee<i-1) { a[nr]=i ; ee=e[nr]=i ; nr++ ;} 
           else    ee=e[nr-1]=i ;
    }
  }
}

/*
void IndexSet::Store(AdaptiveLE *al, vector<double> *v, double *d,int p) {
 int  ix ,sx=(*v).size() ;
 long jx ; 
 AdaptiveLE::iterator alit,alend=(*al).end() ;

 for (alit=(*al).begin();alit!=alend ;alit++) {
   ix=(*alit).first  ;
   jx=(*alit).second ;
   //assert(jx<sx) ;
   //if (jx>=sx) printf ("Store: %d,%d\n",jx,sx) ;
   //assert(jx>=0) ;
   //assert(ix>=0) ;
   //assert(ix<p)  ;
   if (jx>=0) (*v)[jx]=d[ix] ;
  }
}
*/

void IndexSet::Store(AdaptiveLE *al, vector<double> *v, double *d, bool add) {
 int  i ;
 long j ; 
 AdaptiveLE::iterator alit,alend=(*al).end() ;

 for (alit=(*al).begin();alit!=alend ;alit++) {
   i=(*alit).first  ;
   j=(*alit).second ;
   if (j>=0) 
     if (add) (*v)[j] +=d[i] ;
     else     (*v)[j]  =d[i] ;
  }
}

// load only the few coefficients near the boundaries, i.e. [0..l-1] , [r+1,..end()]
// it is assumed, that if any then the complete ranges [0..l-1] , [r+1,..end()] are active
void IndexSet::LoadBoundary(AdaptiveLE *al, vector<double> *v, double *d,int l,int r) {
 int  i ;
 long j ;
 AdaptiveLE::iterator          alit ;
 AdaptiveLE::reverse_iterator ralit ;

 nr=0 ;
 if ((*al).empty()) return;
 
 alit=(*al).begin(); 
 i=(*alit).first   ;

 if (i<l) { 
   assert(i==0) ;
   nr=1 ;
   a[0]=0;e[0]=l-1;
   j=(*alit).second ;
   d[0]=(*v)[j] ;
   alit++ ;
   for (; i<l-1; alit++) {
     i=(*alit).first  ;
     j=(*alit).second ;
     d[i] = (j<0) ? (*v)[1-j] : (*v)[j] ;
   }
 }

 ralit=(*al).rbegin(); 
 i=(*ralit).first   ;

 if (i>r) { 
   a[nr]=r+1;e[nr]=i;
   nr++ ;
   j=(*ralit).second ;
   d[i]=(*v)[j] ;
   ralit++ ;
   for (; i>r+1; ralit++) {
     i=(*ralit).first  ;
     j=(*ralit).second ;
     d[i] = (j<0) ? (*v)[1-j] : (*v)[j] ;
   }
 }
}


void IndexSet::LoadBoundaryIS(AdaptiveLE *al , int l,int r) {
 int  i ;
 AdaptiveLE::iterator          alit ;
 AdaptiveLE::reverse_iterator ralit ;

 nr=0 ;
 if ((*al).empty()) return;
 
 alit=(*al).begin(); 
 i=(*alit).first   ;
 if (i<l) { 
   assert(i==0) ;
   nr=1 ;
   a[0]=0;e[0]=l-1;
   alit++ ;
   for (; i<l-1; alit++) i=(*alit).first ;
 }

 ralit=(*al).rbegin(); 
 i=(*ralit).first   ;
 if (i>r) { 
   a[nr]=r+1;e[nr]=i;
   nr++ ;
   ralit++ ;
   for (; i>r+1; ralit++) i=(*ralit).first  ;
 }
}


void IndexSet::StoreBoundary(AdaptiveLE *al, vector<double> *v, double *d,int l,int r,bool add) {
 int  i ;
 long j ;
 AdaptiveLE::iterator          alit ;
 AdaptiveLE::reverse_iterator ralit ;

 if ((*al).empty()) {assert(nr==0); return;}
 
 alit=(*al).begin();
 i=(*alit).first   ;
 if (i<l) {
   j=(*alit).second ;
   if (add) (*v)[j]+=d[i]; else (*v)[j]=d[i] ;
   alit++ ;
   for (; i<l-1; alit++) {
     i=(*alit).first  ;
     j=(*alit).second ;
     if (add) (*v)[j]+=d[i]; else (*v)[j]=d[i] ;
   }
 }

 ralit=(*al).rbegin(); 
 i=(*ralit).first    ;
 if (i>r) {
   j=(*ralit).second ;
   if (add) (*v)[j]+=d[i]; else (*v)[j]=d[i] ;
   ralit++ ;
   for (; i>r+1; ralit++) {
     i=(*ralit).first  ;
     j=(*ralit).second ;
     if (add) (*v)[j]+=d[i]; else (*v)[j]=d[i] ;
   }
 }
}





/*
//
// fast insert of IndexSet into the AdaptiveLevels:AdaptiveL datastructure
void IndexSet::WriteToAdaptiveLE(AdaptiveLE *al,int *counter) {
  int    r,i=0 ;
  size_t c=(size_t)(*counter) ;
  //AdaptiveLE::iterator hint,ins  ;
 
  
  //(*al).clear() ;
  hint=(*al).end() ;
  
  for (r=nr-1;r>=0;r--)
    for (i=e[r];i>=a[r]; i--) {
      pair<const int,size_t> p(i,c) ;
      ins=(*al).insert(hint,p) ;
      hint=ins ;
      //(*al)[i]=c ;
      c++ ;
     }
  
  *counter=c ;
}
*/


void IndexSet::GenerateRandom(bool boundflag , int KW , int J , bool isW , bool dorandom) {
  int i,r,m ;
  int nJ=1<<J , M= isW ? nJ-1 : nJ ;

  // generate full indexset
  if (!dorandom) {
    M = (!boundflag) ? nJ-1 : nJ ;
    if (isW) M=nJ-1 ;
    nr=1 ;
    a[0]=0; 
    e[0]=M ;
    return ;
  }


  m=(int)(mydrand()*8) ;
  nr=0 ;
  a[0]=(int)(mydrand()*8) ;
  for (i=0; i<m ;i++) {
    e[nr]=min(a[nr] + (int)(mydrand()*10) , M) ;
    if (e[nr] >= a[nr]) { 
      nr++ ;
      a[nr] = e[nr-1] + (int)(mydrand()*10) + 1  ;
    }
  }
 
  if (boundflag && (nr>0)) {
    r=0 ;
    while ((r<nr) && (a[r]<KW)) r++ ;
    if (r>0) {
      r-- ;
      for (i=r; i<nr ;i++) a[i-r]=a[i] , e[i-r]=e[i] ;
      a[0]=0 ;
      nr=nr-r  ;
    }
  
    if (e[0]<KW-1) {
      e[0]=KW-1 ;
      if ( (nr>1) && (a[1]==KW) ) {
        for (i=1; i<nr ;i++) a[i-1]=a[i] , e[i-1]=e[i] ;
        a[0]=0  ;
        nr--    ;
      }
    }

    r=nr-1 ;
    while ((r>=0) && (e[r]>M-KW)) r-- ;

    if (r<nr-1) {
      r++ ;
      e[r]=M ;
      nr=r+1  ;
   
      if (a[nr-1]>M-KW+1) {
        a[nr-1]=M-KW+1 ;
         if ((nr>1)&&(e[nr-2]==M-KW)) { 
           e[nr-2]=M ;
           nr-- ;
          } 
       }
    }
  }
}

---