Department of Scientific Computing    Institute for Numerical Simulation    University of Bonn

Documentation

Download

Programming References

Bug Reports / Suggestions

FAQ

Authors

           

---

Matrix.hpp

//
//
//  Templates for  D-dimensional arrays  
//
//

#ifndef _DARRAY_
# define _DARRAY_

#include "Math.hpp"
#include "UDF.h"
#include "Wavelet.hpp"

#include "Function.hpp"
#include "SparseMatrixIO.hpp"
#include "UDFInformation.hpp"

template<class I,int DIM>
class Matrix {
public:
                 Matrix() ;
                 Matrix(int *begin,int *end) ;
        void     Init  (int *begin,int *end) ;
                ~Matrix() ;
        void     Clear () ;

struct iterator {
public:
        int      i[DIM] ;
                 iterator() ;
                 iterator(Matrix<I,DIM> *d) ;
                 iterator(Matrix<I,DIM> *d,int constcoordinate) ;
                 iterator(const int *begin,const int *end) ;
                 iterator(const int *begin,const int *end , int constcoordinate) ;
       iterator& operator++() ;
        bool     operator!=(const  Matrix<I,DIM>::iterator& x) const ;
        bool     operator<=(const  Matrix<I,DIM>::iterator& x) const ;
        bool     operator<=(const  int b[DIM]) const ;
        bool     operator>=(const  int b[DIM]) const ;
        void     Set      (const int *j) ;
//        void     Set      (const int *j , const int drop  ) ;
        void     SetBounds(const int *b,const int *e) ;
        void     Print() ;

friend class Matrix<I,DIM> ;
  // private:
   int b[DIM],e[DIM],a[DIM],dim ;
} ;

  inline I&       operator[](const iterator &it) ;
  inline I&       operator[](const int *in) ;
  inline I&       operator[](const unsigned int  *in) ;
  inline I&       operator[](const unsigned char *in) ;
  inline iterator begin()         const ;
  inline iterator begin(int constcoordinate) const ; 
  inline iterator end  ()         const ;
  inline int      size (const int i) const ;

  void            CalcLevel(int *J) ;
  void            SameSize (Matrix<I,DIM> *X) ;
  void            Print(int l=0) ;
  void            Set  (const I& i)  ;
  void            Set  (Function *F) ;
  void            Copy (Matrix<I,DIM> *X) ;
  void            Set  (Matrix<I,DIM> *X) ;

  void            Abs  (Matrix<I,DIM> *X) ;
  void            Pow  (Matrix<I,DIM> *X, double p) ;
  void            Sqr  (Matrix<I,DIM> *X) ;
  void            PPlus(Matrix<I,DIM> *X) ;
  void            PPlus(const double a , Matrix<I,DIM> *X) ;
  void            Add  (Matrix<I,DIM> *X , Matrix<I,DIM> *Y) ;
  void            Add  (const double a , Matrix<I,DIM> *X , const double b , Matrix<I,DIM> *Y) ;
  void            Add  (const double a , Matrix<I,DIM> *X , const double b , Matrix<I,DIM> *Y , const double c , Matrix<I,DIM> *Z) ;
  void            Add  (const double a , Matrix<I,DIM> *X , const double b , Matrix<I,DIM> *Y , const double c , Matrix<I,DIM> *Z , const double d , Matrix<I,DIM> *Q) ;
  void            Add  (const double a , Matrix<I,DIM> *X , const double b , Matrix<I,DIM> *Y , const double c , Matrix<I,DIM> *Z , const double d , Matrix<I,DIM> *Q , const double f , Matrix<I,DIM> *R) ;
  void            MMinus(Matrix<I,DIM> *X) ;
  void            MMinus(const double a, Matrix<I,DIM> *X) ;
  void            Sub(Matrix<I,DIM> *X , Matrix<I,DIM> *Y) ;
  double          InnerProd(Matrix<I,DIM> *X) ;
  void            Mul  (Matrix<I,DIM> *X , Matrix<I,DIM> *Y) ;
  void            Mul  (Matrix<I,DIM> *X , Matrix<I,DIM> *Y , const double f) ;
  void            Mul  (Matrix<I,DIM> *X , const double f) ;
  void            Mul  (const double f) ;

  void            Multiply(Matrix<I,DIM> *X , Matrix<I,DIM> *Y) ;

        void     XplusYmalZ (Matrix<I,DIM>*X , Matrix<I,DIM>*Y , Matrix<I,DIM>*Z) ;
        void     XminusYmalZ(Matrix<I,DIM>*X , Matrix<I,DIM>*Y , Matrix<I,DIM>*Z) ;

  double          Min() ;
  double          Max() ;
  double          MaxAbs() ;
  double          Sum()    ;
  double          SumAbs() ;
  double          SumSqr() ;
  void            Transpose(Matrix<I,DIM> *X) ;// DIM must be 2

  void            Smooth(Matrix<I,DIM> *X, int type=0) ;

  typedef void (*UniDirOperation)(double *d,int *BCd, int *dec, double *c,int *BCc, Wavelets *From,int J,int lev , double DX) ;
 
  void            WriteUDF(const char *filename, Extensions<DIM> *E, const char *ScalarOrVectorName,
                           Matrix<I,DIM> **M,int NM, bool castfloat=false, bool ascii=false) ;
  void            ReadUDF (const char *filename, Extensions<DIM> *E, const char *ScalarOrVectorComponentName) ;

  // writes a sparse matrix of all elements |Aij|>=eps
  // file can be read by ReadSparse.m
  //void           WriteMATLABSparse(const char *filename, double eps) ;
  void            ReadSparse (const char *filename, double *XA, double *XE , int Level0, int which) ;
  void            WriteSparse(const char *filename, double *XA, double *XE , int Level0, 
                              Matrix<I,DIM> **M=NULL, int N=0, bool CastToFloat=false) ;

  void           ReadMATLAB(const char *filename) ;

friend struct iterator ;

public:
 int b[DIM],e[DIM],o[DIM],o1 ;
 I*  s              ;
} ;


//
// iterator
//

//template<class I,int DIM>
//struct Matrix<I,DIM>::iterator 

/*****************************************/
/*                                       */
/*           Implementierungen           */
/*                                       */
/*****************************************/

//
// Matrix 
template<class I,int DIM>
Matrix<I,DIM>::Matrix()  
{ 
 for (int k=0;k<DIM;k++) { b[k]=0 , e[k]=-1 ;} 
 s=NULL ;
}

template<class I,int DIM>
Matrix<I,DIM>::Matrix(int *bb,int *ee) { s=NULL ; Init(bb,ee) ;}

template<class I,int DIM>
void Matrix<I,DIM>::Init(int *bb,int *ee) {
  int   c=1,k ; 

  if (s!=NULL) {
    bool resize=false ;
    for (k=0;k<DIM;k++) if ((ee[k]-bb[k])!=(e[k]-b[k])) resize=true ;
    if (resize) {
      std::cout<<"Matrix Resize!\n" ;
      Clear() ;
    }
  }

  for (k=0;k<DIM;k++) { b[k]=bb[k] , e[k]=ee[k] , c*=(e[k]-b[k]+1) ; assert (c>0) ;}
  c*=sizeof(I) ;

  if (s==0) s=(I *)malloc(c) ;
  if (s==NULL) {
    std::cout<<"Matrix<"<<DIM<<">::Init could not allocate "<<c<<" bytes \n" ;
    exit(-1) ;
  }

  o1=0 ;
  c=1  ;
  for (k=DIM-1;k>=0;k--) { o[k]=c , c*=(e[k]-b[k]+1) ; o1 +=o[k]*b[k] ;}  
}


template<class I,int DIM>
void Matrix<I,DIM>::Clear()  { free(s) ; for (int k=0;k<DIM;k++) b[k]=0 , e[k]=-1 ; s=NULL ; }

template<class I,int DIM>
Matrix<I,DIM>::~Matrix()  { 
  Clear() ;
}

template<class I,int DIM>
inline I& Matrix<I,DIM>::operator[](const Matrix<I,DIM>::iterator &it) { return (*this)[&it.i[0]] ;}

template<class I,int DIM>
inline I& Matrix<I,DIM>::operator[](const int *in) {
  int i=in[DIM-1]-o1 ;
  for (int k=DIM-2 ;k>=0;k--) i+= in[k]*o[k] ;
  return s[i] ;
}

template<class I,int DIM>
inline I& Matrix<I,DIM>::operator[](const unsigned int *in) {
  int i=in[DIM-1]-o1 ;
  for (int k=DIM-2 ;k>=0;k--) i+= in[k]*o[k] ;
  return s[i] ;
}

template<class I,int DIM>
inline I& Matrix<I,DIM>::operator[](const unsigned char *in) {
  int i=in[DIM-1]-o1 ;
  for (int k=DIM-2 ;k>=0;k--) i+= in[k]*o[k] ;
  return s[i] ;
}


template<class I,int DIM>
Matrix<I,DIM>::iterator Matrix<I,DIM>::begin() const { return begin(-1) ; } 

template<class I,int DIM>
Matrix<I,DIM>::iterator Matrix<I,DIM>::begin(int drop) const {
 iterator tmp ;
 tmp.dim=0 ;
 for (int k=0;k<DIM;k++) { tmp.b[k]=tmp.i[k]=b[k] ; tmp.e[k]=e[k] ; if(k!=drop) tmp.a[tmp.dim++]=k ;}
 return tmp ;
} 


template<class I,int DIM>
Matrix<I,DIM>::iterator Matrix<I,DIM>::end()   const {
 iterator tmp    ;
 for (int i=0;i<DIM;i++) tmp.i[i]=e[i] ;
 return tmp ;
} 

template<class I,int DIM>
int Matrix<I,DIM>::size(const int i) const {
  if ((0<=i) && (i<DIM)) return e[i]-b[i]+1 ;
  else                   return 0 ;
}

template<class I,int DIM>
void Matrix<I,DIM>::Print(int level) 
{int i ;
 if (level>=0) {
   std::cout << "Dimensionen :" ;
   for (i=0;i<DIM;i++) if (i<DIM-1) std::cout<<'['<<b[i]<<','<<e[i]<<"]x"; else std::cout<<'['<<b[i]<<','<<e[i]<<"]\n";
 }
 if ((level>=1) && (DIM==2)) {
   int i,j ; 
   for (i=b[0]; i<=e[0]; i++) {
   for (j=b[1]; j<=e[1]; j++) {
     int ij[2]={i,j} ;
     std::cout<<(*this)[ij]<<"  " ;    
    }
   std::cout<<'\n';
  }
 }
}

//                              //
//       Matrix::iterator       //
//                              // 


template<class I,int DIM>
Matrix<I,DIM>::iterator::iterator() { 
  dim=DIM ;
  for (int k=0;k<DIM;k++) {i[k]=b[k]=e[k]=0 , a[k]=k ;} 
  i[0]=-1 ;
}

template<class I,int DIM>
Matrix<I,DIM>::iterator::iterator(const int *bb,const int *ee) {
  SetBounds(bb,ee) ;
  dim=DIM ;
  for (int k=0;k<DIM;k++) a[k]=k ;
}

template<class I,int DIM>
Matrix<I,DIM>::iterator::iterator(Matrix<I,DIM> *d) {
  SetBounds(d->b,d->e) ;
  dim=DIM ;
  for (int k=0;k<DIM;k++) a[k]=k ;  
}

template<class I,int DIM>
Matrix<I,DIM>::iterator::iterator(Matrix<I,DIM> *d,int drop) { 
  SetBounds(d->b,d->e) ;
  dim=0 ;
  for (int k=0;k<DIM;k++) { if(k!=drop) a[dim++]=k ;} 
}

template<class I,int DIM>
Matrix<I,DIM>::iterator::iterator(const int *bb,const int *ee , int drop) { 
  SetBounds(bb,ee) ;
  dim=0 ;
  for (int k=0;k<DIM;k++) { if(k!=drop) a[dim++]=k ;} 
}


template<class I,int DIM>
void Matrix<I,DIM>::iterator::SetBounds(const int *bb,const int *ee) { 
  for (int k=0;k<DIM;k++) {i[k]=b[k]=bb[k] , e[k]=ee[k] ;} 
}

template<class I,int DIM>
void Matrix<I,DIM>::iterator::Set(const int *j) {
  for (int k=0;k<DIM;k++) i[k]=j[k] ; 
}

/*
template<class I,int DIM>
void Matrix<I,DIM>::iterator::Set(const int *j , const int drop) {
  // dim=0 ;
  for (int k=0;k<DIM;k++) { i[k]=j[k] }
}
*/


template<class I,int DIM>
Matrix<I,DIM>::iterator& Matrix<I,DIM>::iterator::operator++() 
{ int k=0,c ;
  if (dim==0) { i[0]=e[0]+1 ; return *this ;} 
  do { i[a[k]]++ ;
    if(c=(i[a[k]]>e[a[k]])) {i[a[k]]=b[a[k]] , k++;}
  } while (c && (k<dim)) ;
  if (k==dim) i[a[0]]=e[a[0]]+1 ; 
  return *this ; 
} 


template<class I,int DIM> 
bool Matrix<I,DIM>::iterator::operator!=(const  Matrix<I,DIM>::iterator& x) const {
  for (int k=0;k<dim;k++) if (i[a[k]] != x.i[a[k]]) return true ;
  return false ; 
}


template<class I,int DIM>
bool Matrix<I,DIM>::iterator::operator<=(const  Matrix<I,DIM>::iterator& x) const {
  return (*this)<=x.i ;
}

template<class I,int DIM>
bool Matrix<I,DIM>::iterator::operator<=(const  int b[DIM]) const {
  for (int k=0;k<dim;k++) if (i[a[k]]>b[a[k]]) return false ;
  if (dim==0) return i[0]<=b[0] ;
  return true ;  
}

template<class I,int DIM>
bool Matrix<I,DIM>::iterator::operator>=(const  int b[DIM]) const {
  for (int k=0;k<dim;k++) if (i[a[k]]<b[a[k]]) return false ;
  return true ;  
}

template<class I,int DIM>
void Matrix<I,DIM>::iterator::Print() {
  int j ;
  std::cout<<"Dim: "<<dim ;
  std::cout<<"\nb: "; for (j=0;j<D;j++) std::cout<<b[j]<<"  " ;
  std::cout<<"\ne: "; for (j=0;j<D;j++) std::cout<<e[j]<<"  " ;
  std::cout<<"\na: "; for (j=0;j<D;j++) std::cout<<a[j]<<"  " ;
  std::cout<<"\ni: "; for (j=0;j<D;j++) std::cout<<i[j]<<"  " ;
  std::cout<<'\n' ;
}



//                              //
//       Matrix::E/A            //
//                              // 


template<class I,int DIM>
void Matrix<I,DIM>::Set(const I &x) {  
 Matrix<I,DIM>::iterator it,tend=(*this).end() ;
 for (it=(*this).begin();it<=tend;++it) (*this)[it]=x ;
}


template<class I,int DIM>
void Matrix<I,DIM>::Set(Function *F) {
 int i ;
 Matrix<I,DIM>::iterator it(this) ;

 for (it=(*this).begin(); it<=(*this).end(); ++it) {
   double x[DIM] ;
   for (i=0; i<DIM ;i++) x[i]=((double)it.i[i] - b[i])/(e[i]-b[i]) ;
   (*this)[it]=(I)F->Eval(x) ;
 }
}


template<class I,int DIM>
void Matrix<I,DIM>::WriteUDF(const char *filename, Extensions<DIM> *Ext, const char *ScalarOrVectorName,
                             Matrix<I,DIM> **M,int NM, bool castfloat, bool ascii)
{
 // scalars or vector ?
 Matrix<I,DIM> *MM[3] ;
 bool vec   =false ;
 UDFInformation<DIM> UDF ;

 int i,j ;
 if (M!=NULL) {
   vec=true ;
   assert(NM >0) ;
   assert(NM<=3) ;
   for (i=1; i<NM; i++) {
     for (j=0; j<DIM; j++) assert( (M[0]->e[j]==M[i]->e[j]) && (M[0]->b[j]==M[i]->b[j]) ) ;
     MM[i]=M[i] ;
   }
 } else {
   MM[0]=this ;
   NM=1 ;
 }

 // dimensions
 int eb[DIM+3]={0},count=1 ;
 for (i=0; i<DIM; i++) { eb[i]=MM[0]->e[i]-MM[0]->b[i]+1 ; count *=eb[i] ;}
 eb[i]=1; eb[i+1]=1 ;

 // write VTK header
 FILE *fid=fopen(filename,"w") ;
 if (fid==NULL) std::cout<<"could not open for write "<<filename<<'\n' ;
 assert(fid) ;

 const char *dataType = (castfloat) ? "float" : "double" ;

 fprintf(fid,"# vtk DataFile Version 2.0\n") ;
 fprintf(fid,"UDF-DataSetEXT00")             ;
 Ext->WriteUDF(fid) ;

 if (ascii) fprintf(fid,"ASCII\n")           ;
     else   fprintf(fid,"BINARY\n")          ;

 fprintf(fid,"DATASET RECTILINEAR_GRID\n")   ;
 fprintf(fid,"DIMENSIONS ") ;
 for (i=0; i<DIM; i++) fprintf(fid,"%d ",eb[i]) ;
 for (i=DIM; i<3; i++) { fprintf(fid,"%d ",1) ; eb[i]=1; }
 fprintf(fid,"\n") ;

 double x  ;
 float  fx ;
 for (i=0; i<max(DIM,3); i++) {
   switch (i) {
   case 0: fprintf(fid,"X_COORDINATES %d %s\n",eb[i],dataType) ;
   break ;
   case 1: fprintf(fid,"Y_COORDINATES %d %s\n",eb[i],dataType) ;
   break ;
   case 2: fprintf(fid,"Z_COORDINATES %d %s\n",eb[i],dataType) ;
   break ;
   default:fprintf(fid,"X%d_COORDINATES %d %s\n",i,eb[i],dataType) ;
   break ; 
  }
  
  double dx ;
  if (i<DIM) dx=(Ext->XE[i] - Ext->XA[i]) / (eb[i]-1) ;  else dx=1.0 ;
  for (j=0; j<eb[i]; j++) {
    if (i<DIM) x=Ext->XA[i] + j*dx ; else x=j*dx ;
    if (ascii) {
      if (castfloat) fprintf(fid,"%e\n",x) ;
          else       fprintf(fid,"%18.16e\n",x) ;
    } else {
      if (castfloat) { 
        fx=x ;
        if (UDF.do_swap) UDF_swap(&fx,sizeof(float),1) ;
        fwrite(&fx,sizeof(float),1,fid) ;
      } else {
        if (UDF.do_swap) UDF_swap(&x,sizeof(double),1) ;
        fwrite(&x,sizeof(double),1,fid) ;
      }
    }
  } // j
  if (!ascii) fprintf(fid,"\n");
 } // i 

 //
 // Write Data
 //
 fprintf (fid,"POINT_DATA %d\n",count) ;
 
 if (!vec) {
    fprintf (fid,"SCALARS %s %s\n",ScalarOrVectorName,dataType) ;
    fprintf (fid,"LOOKUP_TABLE default\n") ;
 } else {
    fprintf (fid,"VECTORS %s %s\n",ScalarOrVectorName,dataType) ;
 }

 Matrix<I,DIM>::iterator it(MM[0]), tend=(*MM[0]).end() ;

 float mi=1e+20,ma=-mi ;
 for (it=(*MM[0]).begin(); it<=tend; ++it) 
  {
    int je = (vec) ? max(NM,3) : NM;
    for (j=0; j<je; j++) {
      if (vec && (j>=NM)) x=0.0 ;
          else            x=(double)(*MM[j])[it] ;

      mi=x<mi ? x : mi ;
      ma=x>ma ? x : ma ;

      if (ascii) {
        if (castfloat) fprintf(fid,"%e\n",x) ;
            else       fprintf(fid,"%18.16e\n",x) ;
      } else {
        if (castfloat) { 
          fx=(float)x ;
          if (UDF.do_swap) UDF_swap(&fx,sizeof(float),1) ;
          fwrite(&fx,sizeof(float),1,fid) ;
        } else {
          if (UDF.do_swap) UDF_swap(&x,sizeof(double),1) ;
          fwrite(&x,sizeof(double),1,fid) ;
        }
      } // binary
    } // j
  } // it
 std::cout<< "WriteUDF Min= "<<scientific<<mi<<' '<<"Max= "<<scientific<<ma<<'\n' ;
 fclose(fid) ;
}


template<class I,int DIM>
void Matrix<I,DIM>::ReadUDF(const char *filename, Extensions<DIM> *Ext,
                            const char *) 
{
 int i ;
 char st[1000],fd[100],c,st1[1000] ;
 double x  ;
 float  fx ;

 UDFInformation<DIM> UDF ;

 FILE *fid=fopen(filename,"r") ;
 assert(fid) ;
 
 UDF.ReadHeader(fid, Ext) ;
 // Dimensions
 int A[DIM]={0},E[DIM] ;
 for (i=0; i<DIM; i++) E[i]=UDF.size[i]-1 ;
 Init(A,E) ;
 
 // Read Data
 fscanf(fid,"%s %s",st,st1)       ; // POINT_DATA xx
 fscanf(fid,"%s %s %s",st,st1,fd) ; // SCALARS bla
 fscanf(fid,"%s %s",st,st1)       ; // LOOKUP_TABLE default
 fscanf(fid,"%c",&c)              ; // EOL 
  
 Matrix<I,DIM>::iterator it(this) ;
 
 for (it=(*this).begin(); it<=(*this).end(); ++it) {

   if (UDF.ascii) {
     fscanf(fid,"%le",&x) ;
     (*this)[it]=(I)x ;
    } else {
     if (strcmp(fd,"float")==0) { 
       fread(&fx,sizeof(float),1,fid) ;
       if (UDF.do_swap) UDF_swap(&fx,sizeof(float),1) ; 
       (*this)[it]=(I)fx ;
     } else {
       fread(&x,sizeof(double),1,fid) ;
       if (UDF.do_swap) UDF_swap(&x,sizeof(double),1) ; 
       (*this)[it]=(I)x ;
     }
    } // binary

 } // it
 fclose(fid) ;
}


// SparseMatrix -> Matrix
template<class I,int DIM>
void Matrix<I,DIM>::ReadSparse(const char *name, double *XA, double *XE, int Level0, int which) {

 int           i,l[DIM],J[DIM],N ;
 unsigned long x [DIM],s[DIM]    ;
 long          j ;
 double        wd[100] ;
 float         wf[100] ;
 bool          insert  ;

 int header[200]={0},anz,jm1 ;

 bool CastToFloat, do_swap ;

 Set(0.0) ;
 CalcLevel(J) ;

 FILE *fid=fopen(name,"r") ;
 if (fid==NULL) {
   std::cout<<"Matrix<"<<DIM<<">::ReadSparse::Could not open "<<name<<'\n';
   exit(-1) ;
 }

 // read header
 ReadMATLABSparseHeader(fid,DIM, header, &anz, &jm1, &N,XA,XE,&CastToFloat,&do_swap) ;

 // read all entries
 for (i=0; i<anz; i++) {
   fread( x,sizeof(x[0]),DIM,fid) ;
   if (do_swap) UDF_swap( x,sizeof(x[0]),DIM) ;

   if (CastToFloat) {
     fread(wf,sizeof(float),N,fid) ;
     if (do_swap) UDF_swap(wf,sizeof(float),DIM) ;
     for (j=0; j<N; j++) wd[j]=(double)wf[j] ; 
   } else {
     fread(wd,sizeof(double),N,fid) ;
     if (do_swap) UDF_swap(wd,sizeof(double),DIM) ;
   }
   
   insert=true ;
   // level l ,index s decodieren
   for (j=0; j<DIM; j++) { 
     decode(&l[j],&s[j] , x[j] , jm1,Level0) ;
     if (l[j]>J[j]) insert=false ;
    }

   // insert (l,s)
   if (insert) {
     int t[DIM] ;
     for (j=0; j<DIM; j++) 
       t[j] = s[j] + ( (l[j]==Level0) ? 0 : 1<<(l[j]-1) ) ;
     
     (*this)[t]=(I)wd[which] ;
   }

 } // i 
 fclose(fid) ;
}

template<class I,int DIM>
void Matrix<I,DIM>::WriteSparse(const char *filename, double *XA, double *XE , int Level0, 
                                      Matrix<I,DIM> **MM, int N, bool CastToFloat) {

  int  i,header[100]={0},anz, hs,J[DIM],JM ;
  bool do_swap = UDF_CPUisBigEndian() > 0 ;

  Matrix<I,DIM> *M[100] ;
  if (MM==NULL) { M[0]=this; N=1 ;}
  else {
    for (i=0; i<N; i++) {
      assert(MM[i]) ;
      M[i]=MM[i] ;
    }
  }
  
  FILE *fid=fopen(filename,"w") ;
  if (fid==NULL) {
    std::cout<<"WriteSparse::Could not open "<<filename<<'\n' ;
    exit(-1) ;
  }

  // number of elements
  CalcLevel(J) ;
  JM =0 ;
  anz=1 ;
  for (i=0; i<DIM; i++) {
    anz *=size(i) ;
    JM=max(JM,J[i]) ;
  }

  hs=header[0]=sizeof(header)/sizeof(header[0]) ;
  header[1]=DIM ;
  header[2]=anz ;
  header[3]=JM  ;
 
  float *p=(float *)(&(header[4])) ;
  assert(sizeof(int)==sizeof(float)) ;

  for (i=0; i<DIM; i++) {
    *p=(float)XA[i]; p++ ;
    *p=(float)XE[i]; p++ ;
  }
  
  header[4+2*DIM] = N ; // number of components
  header[5+2*DIM] = (CastToFloat ? 1 : 0) ;
  
  if (do_swap) UDF_swap(header,sizeof(int),hs) ;
  fwrite(header,sizeof(int) , hs , fid) ;


  //
  // write actual data
  // 

  int      s=0,l=0,j ;
  unsigned long x[DIM] ;
    
  Matrix<I,DIM>::iterator t(this) ;
  for (t=begin(); t<=end(); ++t) {
    // calc level & index  and code
    for (i=0; i<DIM; i++) {
      if (t.i[i]<=(1<<Level0)) { l=Level0; s=t.i[i]; }
      else 
        for (j=Level0+1; j<=JM; j++) 
          if (t.i[i]<=(1<<j))  { l=j     ; s=t.i[i]-(1<<(j-1)) -1 ; break ; }

      code(l,s,&x[i],JM,Level0) ;
    }

    if (do_swap) UDF_swap( x,sizeof(x[0]),DIM) ;
    fwrite( x,sizeof(x[0]),DIM,fid) ;

    if (CastToFloat) {
      float  w[100] ;
      for (i=0; i<N; i++) w[i]=(float)(*M[i])[t] ; 

      if (do_swap) UDF_swap(w,sizeof(float),DIM) ;
      fwrite(w,sizeof(float),N,fid) ;

    } else {

      double  w[100] ;
      for (i=0; i<N; i++) w[i]=(double)(*M[i])[t] ; 

      if (do_swap) UDF_swap(w,sizeof(double),DIM) ;
      fwrite(w,sizeof(double),N,fid) ;
    }
    
  } // i 
  fclose(fid) ;
  
}


// read file in outdated MATLAB format, maintained for compatibility reasons only
template<class I,int DIM>
void Matrix<I,DIM>::ReadMATLAB(const char *filename)
{
 FILE *fid=fopen(filename,"r") ;
 if (fid==NULL) { std::cout<<"Matrix<"<<DIM<<">::ReadMATLAB: could not not open file: "<<filename<<'\n'; exit(-1) ; }
  
 int dim,begin[DIM],end[DIM] ;

 // read dimensions
 fread(&dim,sizeof(int),1,fid) ; 
 assert(dim==DIM)              ;

 fread(begin,sizeof(int),dim,fid) ; 
 fread(end  ,sizeof(int),dim,fid) ;
  
 Init(begin,end) ;

 iterator it(this),tend=(*this).end() ;

 for (it=(*this).begin();it<=tend;++it) 
  { double x ; 
    fread(&x,sizeof(double),1,fid) ;
    (*this)[it]=(I)x ;
  } 
 fclose(fid) ;
}



//
//
//

template<class I,int DIM>
void Matrix<I,DIM>::CalcLevel(int *J) {
 int h,d ;
 for (d=0; d<DIM; d++) {
    h=e[d]-b[d]+2 ;
    J[d]=0 ; while(h=(h>>1)) J[d]++ ; 
  }
}

//
// Linear Algebra
//
template<class I,int DIM>
void Matrix<I,DIM>::SameSize(Matrix<I,DIM> *X) {
  for (int i=0;i<DIM;i++) assert((b[i]==X->b[i])&&(e[i]==X->e[i])) ;
}

template<class I,int DIM>
double Matrix<I,DIM>::Min() {
  double mx=1e+200,x ;
  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) { 
    x =(double)(*this)[it] ;
    mx= (mx<x) ? mx : x ;
  }
 return mx ;
}

template<class I,int DIM>
double Matrix<I,DIM>::Max() {
  double mx=-1e+200,x ;
  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) { 
    x =(double)(*this)[it] ;
    mx= (mx>x) ? mx : x ;  
  }
 return mx ;
}

template<class I,int DIM>
double Matrix<I,DIM>::MaxAbs() {
  double mx=0.0 ,x  ;
  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) { 
    x =fabs((double)(*this)[it]) ;
    mx= (mx>x) ? mx : x ;  
  }
 return mx ;
}

template<class I,int DIM>
double Matrix<I,DIM>::Sum() {
  double mx=0.0 ;
  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) mx +=(double)(*this)[it] ;
  return mx ;
}

template<class I,int DIM>
double Matrix<I,DIM>::SumAbs() {
  double mx=0.0 ;
  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) mx +=fabs((double)(*this)[it]) ;
  return mx ;
}

template<class I,int DIM>
double Matrix<I,DIM>::SumSqr() {
  double mx=0.0 ;
  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) mx +=((double)(*this)[it])*((double)(*this)[it]) ;
  return mx ;
}

template<class I,int DIM>
void Matrix<I,DIM>::Copy(Matrix<I,DIM> *X) {
  SameSize(X) ;
  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) (*this)[it]=(*X)[it] ;
}


template<class I,int DIM>
void Matrix<I,DIM>::Set(Matrix<I,DIM> *X) {
  assert(X!=this) ;
  Set((I)0) ;
  int aa[DIM],ee[DIM],i ;
  for (i=0; i<DIM; i++) {
    aa[i]=max(X->b[i], b[i]) ;
    ee[i]=min(X->e[i], e[i]) ;
  }

  Matrix<I,DIM>::iterator it(aa,ee) ;
  for (it.Set(aa); it<=ee; ++it) (*this)[it.i]=(*X)[it.i] ;
}

template<class I,int DIM>
void Matrix<I,DIM>::Abs(Matrix<I,DIM> *X) {
  SameSize(X) ;
  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) (*this)[it]=fabs((*X)[it]) ;
}

template<class I,int DIM>
void Matrix<I,DIM>::Pow(Matrix<I,DIM> *X, double p) {
  SameSize(X) ;
  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) (*this)[it]=pow((*X)[it],p) ;
}

template<class I,int DIM>
void Matrix<I,DIM>::Sqr(Matrix<I,DIM> *X) {
  SameSize(X) ;
  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) (*this)[it]=(*X)[it]*(*X)[it] ;
}

template<class I,int DIM>
void Matrix<I,DIM>::PPlus(Matrix<I,DIM> *X) {
  SameSize(X) ;
  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) (*this)[it] +=(*X)[it] ;
}

template<class I,int DIM>
void Matrix<I,DIM>::PPlus(const double a , Matrix<I,DIM> *X) {
  SameSize(X) ;
  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) (*this)[it] +=a*(*X)[it] ;
}

template<class I,int DIM>
void Matrix<I,DIM>::Add(Matrix<I,DIM> *X , Matrix<I,DIM> *Y) {
  SameSize(X) ;
  SameSize(Y) ;
  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) (*this)[it] =(*X)[it] + (*Y)[it] ;
}


template<class I,int DIM>
void Matrix<I,DIM>::Add(const double a, Matrix<I,DIM> *X ,const double b, Matrix<I,DIM> *Y) {
  SameSize(X) ;
  SameSize(Y) ;
  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) (*this)[it] =a*(*X)[it] + b*(*Y)[it] ;
}

template<class I,int DIM>
void Matrix<I,DIM>::Add(const double a, Matrix<I,DIM> *X ,const double b, Matrix<I,DIM> *Y , const double c, Matrix<I,DIM> *Z) {
  SameSize(X) ;
  SameSize(Y) ;
  SameSize(Z) ;

  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) (*this)[it] =a*(*X)[it] + b*(*Y)[it] +c*(*Z)[it] ;
}

template<class I,int DIM>
void Matrix<I,DIM>::Add(const double a, Matrix<I,DIM> *X ,const double b, Matrix<I,DIM> *Y , const double c, Matrix<I,DIM> *Z , const double d, Matrix<I,DIM> *Q ) {
  SameSize(X) ;
  SameSize(Y) ;
  SameSize(Z) ;
  SameSize(Q) ;

  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) (*this)[it] =a*(*X)[it] + b*(*Y)[it] +c*(*Z)[it] +d*(*Q)[it];
}

template<class I,int DIM>
void Matrix<I,DIM>::Add(const double a, Matrix<I,DIM> *X ,const double b, Matrix<I,DIM> *Y , const double c, Matrix<I,DIM> *Z , const double d, Matrix<I,DIM> *Q , const double e, Matrix<I,DIM> *R ) {
  SameSize(X) ;
  SameSize(Y) ;
  SameSize(Z) ;
  SameSize(Q) ;
  SameSize(R) ;

  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) (*this)[it] =a*(*X)[it] + b*(*Y)[it] +c*(*Z)[it] +d*(*Q)[it] +e*(*R)[it]  ;
}

template<class I,int DIM>
void Matrix<I,DIM>::MMinus(Matrix<I,DIM> *X) {
  SameSize(X) ;
  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) (*this)[it] -=(*X)[it] ;
}

template<class I,int DIM>
void Matrix<I,DIM>::MMinus(const double a,Matrix<I,DIM> *X) {
  SameSize(X) ;
  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) (*this)[it] -=a*(*X)[it] ;
}

template<class I,int DIM>
void Matrix<I,DIM>::Sub(Matrix<I,DIM> *X , Matrix<I,DIM> *Y) {
  SameSize(X) ;
  SameSize(Y) ;
  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) (*this)[it] =(*X)[it] - (*Y)[it] ;
}

template<class I,int DIM>
void Matrix<I,DIM>::Mul(Matrix<I,DIM> *X , Matrix<I,DIM> *Y) {
  SameSize(X) ;
  SameSize(Y) ;
  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) (*this)[it] =(*X)[it] * (*Y)[it] ;
}

template<class I,int DIM>
void Matrix<I,DIM>::Mul(Matrix<I,DIM> *X , Matrix<I,DIM> *Y , const double f) {
  SameSize(X) ;
  SameSize(Y) ;
  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) (*this)[it] =(*X)[it] * (*Y)[it] *f ;
}

template<class I,int DIM>
void Matrix<I,DIM>::Mul(Matrix<I,DIM> *X , const double f) {
  SameSize(X) ;
  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) (*this)[it] =(*X)[it] *f ;
}

template<class I,int DIM>
void Matrix<I,DIM>::Mul(const double f) {
  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) (*this)[it] =(*this)[it] *f ;
}


template<class I,int DIM>
double Matrix<I,DIM>::InnerProd(Matrix<I,DIM> *X) {
  SameSize(X) ;
  double s=0.0 ;
  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) s +=(*this)[it]*(*X)[it] ;
  return s;
}

template<class I,int DIM>
void Matrix<I,DIM>::XplusYmalZ(Matrix<I,DIM> *X , Matrix<I,DIM> *Y , Matrix<I,DIM> *Z) {
  SameSize(X) ;
  SameSize(Y) ;
  SameSize(Z) ;

  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) (*this)[it] =(*X)[it] + (*Y)[it]*(*Z)[it] ;
}

template<class I,int DIM>
void Matrix<I,DIM>::XminusYmalZ(Matrix<I,DIM> *X , Matrix<I,DIM> *Y , Matrix<I,DIM> *Z) {
  SameSize(X) ;
  SameSize(Y) ;
  SameSize(Z) ;

  Matrix<I,DIM>::iterator it(this) ;
  for (it=begin(); it<=end(); ++it) (*this)[it] =(*X)[it] + (*Y)[it]*(*Z)[it] ;
}

// usual product of two 2D compatible matrices 
template<class I,int DIM>
void Matrix<I,DIM>::Multiply(Matrix<I,DIM> *X , Matrix<I,DIM> *Y) {
  assert(DIM==2) ;
  assert(b[0]==X->b[0]) ; assert(e[0]==X->e[0]) ; 
  assert(b[1]==Y->b[1]) ; assert(e[1]==Y->e[1]) ; 
  assert(X->b[1]==Y->b[0]) ; assert(X->e[1]==Y->e[0]) ; 

  double s;
  int i,j,k,ij[2],ik[2],kj[2] ;
  for (i=b[0]; i<=e[0]; i++)
    for (j=b[1]; j<=e[1]; j++) {
      ij[0]=i,ij[1]=j ;
      s=0 ;
      for (k=X->b[1]; k<=X->e[1]; k++) {
        ik[0]=i,ik[1]=k ;
        kj[0]=k,kj[1]=j ;
        s += (*X)[ik]*(*Y)[kj] ;
       }
      (*this)[ij]=s ;
    }
}

template<class I,int DIM>
void Matrix<I,DIM>::Transpose(Matrix<I,DIM> *X) {
  assert(DIM==2) ;
  assert(b[0]==X->b[0]) ; 
  assert(e[1]==X->e[1]) ; 

  Matrix<I,DIM> *Y=X ;
  if (this==X) {
    Y=new Matrix<I,DIM>(X->b, X->e) ;
    Y.Copy(X) ;
  }
 
  int i,j ;
  for (i=b[0]; i<=e[0]; i++)
    for (j=b[1]; j<=e[1]; j++) {
      int ij[2]={i,j}, ji[2]={j,i} ; 
      (*this)[ij]=(*Y)[ji] ;
    }

  if (this==C) delete Y ;
}
#endif

---