MATH_MultiLevelsToeplitzMatrix.h

Go to the documentation of this file.

#ifndef MATH_MultiLevelsToeplitzMatrix_H
#define MATH_MultiLevelsToeplitzMatrix_H

#include "CORE_Object.h"

#include "CORE_Array.h"
#include "CORE_Run.h"

#include "FFTW_FFT.h"

#include "MATH_LeafBlockMatrixStorage.h"

DEFINE_SPTR(MATH_MultiLevelsToeplitzMatrix);

class MATH_MultiLevelsToeplitzMatrix : public CORE_Object  {
    SP_OBJECT(MATH_MultiLevelsToeplitzMatrix);
  
public:
    

    // ATTRIBUTES
public:
    static const tFlag CANONICAL_SPACE=0;
    static const tFlag SPECTRAL_SPACE=1;

private:

    //state of the matrix
    tFlag mSpace;
    

    //number of values in the canonical space
    tUIndex mCanonicalSpaceSize;
    
 
    //the number of levels
    CORE_UIndexArray mLevels;
    
    //the number of rows or column of the matrix
    tUIndex mDimension;
    
      
    
    
    // ASSOCIATIONS


    //the leaf matrix storage
    SP::MATH_LeafBlockMatrixStorage mStorage;

    //values of the toeplitz matrix
    SP::FFTW_FFT mFFT;
    
   
    
protected: 
    
    // METHOD
    
    // CONSTRUCTORS 
    
    MATH_MultiLevelsToeplitzMatrix(void);
    
    
    
    // DESTRUCTORS 
    
  
    virtual ~MATH_MultiLevelsToeplitzMatrix(void);   
    
    
public:
    virtual SP::MATH_MultiLevelsToeplitzMatrix NewInstance() const=0;
    
   
    
  
    //OPERATORS
protected:
    inline tFFTWComplex & operator[](const tUIndex& k) {
        return (*mFFT.get())[k];
    };
    inline const tFFTWComplex & operator[](const tUIndex& k) const {
        return (*mFFT.get())[k];
    };

public:
    inline const tFFTWComplex& operator()(const tUIndex& i,const tUIndex& j) const {
        if (mSpace!=CANONICAL_SPACE) throw CORE_Exception("emicrom/math",
                                                          "MATH_MultiLevelsToeplitzMatrix()",
                                                          "the space is not the CANONICAL SPACE");
        return (*mFFT.get())[getIndex(i,j)];
    };
    
    inline tFFTWComplex & operator()(const tUIndex& i,const tUIndex& j) {
        if (mSpace!=CANONICAL_SPACE) throw CORE_Exception("emicrom/math",
                                                          "MATH_MultiLevelsToeplitzMatrix()",
                                                          "the space is not the CANONICAL SPACE");
        return (*mFFT.get())[getIndex(i,j)];
    };
    
  
   
  
    
    virtual void copy(const MATH_MultiLevelsToeplitzMatrix& mat);
    

    // SET methods
    virtual void clear(); 
    
protected:
    inline void setSpace(const tFlag& space) {
        mSpace=space;
    };

    
    inline void setStorage(SP::MATH_LeafBlockMatrixStorage storage) {
        mStorage=storage;
    }
    
public:
   
    virtual void setLeafBlockDimension(const tUIndex& n) {
        mStorage->setDimension(n);
    }
    
    virtual void setLevels(const tUSInt& levelsNumber,const tUIndex levels[],const tUIndex& nD) {
        setLeafBlockDimension(nD);
        setLevels(levelsNumber,levels);
    }
    
    void setLevels(const tUSInt& levelsNumber,const tUIndex levels[]);


    
   
    
  
protected:

    inline void setDFT(SP::FFTW_FFT dft) {
        mFFT=dft;
        if (dft.get()==null) throw CORE_Exception("emicrom/math/toeplitz",
                                                  "MATH_MultiLevelsToeplitzMatrix::setDFT(null)",
                                                  "impossible to set the Discrete Fourier Transform: FFTW_ClassFactory not added in general class factory");
    }
    
    // GET method
    
public:
    virtual tULLInt getMemorySize() const {
        return (mFFT.get()==null)?0:mFFT->getMemorySize();
    }

protected:
    
    inline const FFTW_FFT& getDFT() const {
        return *mFFT.get();
    }
    inline FFTW_FFT& getDFT()  {
        return *mFFT.get();
    }
  
  
    
public:
    inline tReal getRealValue(const tUIndex& i,const tUIndex& j) const {
        return (*this)(i,j)[0];
    }
    
  
    inline const tUIndex&  getDimension() const {
        return mDimension;
    };
    inline const tUIndex&  getRowsNumber() const {
        return mDimension;
    };
    inline const tUIndex&  getColumnsNumber() const {
        return mDimension;
    };
   

public:

    
    inline FFTW_ComplexArray& getValues()  {
        return mFFT->getArray();
    };
    inline const FFTW_ComplexArray& getValues() const {
        return mFFT->getArray();
    };
    inline const tFFTWComplex* getLeafBlockValues(const tIndex indices[]) const {
        return &getValues()[getToeplitzBlockIndex(mLevels.getSize(),&mLevels[0],indices)*mStorage->getSize()];
        
    }
    inline tFFTWComplex* getLeafBlockValues(const tIndex indices[]) {
        return &getValues()[getToeplitzBlockIndex(mLevels.getSize(),&mLevels[0],indices)*mStorage->getSize()];
        
    }
    
    inline const CORE_UIndexArray& getLevels() const {
        return mLevels;
    };

    inline tUSInt getLevelsNumber() const {
        return mLevels.getSize();
    }

  

    inline const tUIndex& getCanonicalSpaceSize() const {
        return mCanonicalSpaceSize;
    }
   

    inline tUIndex getLeafBlockSize() const {
        return mStorage->getSize();
    }
    inline const tUIndex& getLeafBlockDimension() const {
        return mStorage->getDimension();
    }
    inline tUIndex getBlocksNumber() const {
        return mCanonicalSpaceSize/getLeafBlockSize();
    }
    inline const tUIndex& getBlocksNumberPerRow(const tUSInt& l) const {
        return mLevels[l-1];
    }
    
    inline tUIndex getBlocksNumber(const tUIndex& l) const {
        return 2*mLevels[l-1]-1;
    }
    

    
public:
    inline const tFlag& getSpace() const {
        return mSpace;
    }
    
    inline const MATH_LeafBlockMatrixStorage& getStorage() const {
        return *mStorage.get();
    }
    inline MATH_LeafBlockMatrixStorage& getStorage()  {
        return *mStorage.get();
    }
    
    inline tUIndex getIndex(const tUIndex& i,const tUIndex& j) const {
        tUIndex nValues=0;
        return getIndex(i,j,mLevels.getSize(),getRowsNumber(),nValues);
        
    };
    
    tUIndex getToeplitzBlockIndex(const tUSInt& nLevels,const tUIndex N[],const tIndex I[]) const;
    inline tUIndex getToeplitzBlockIndex(const tIndex I[]) const {
        if (mLevels.getSize()==0) return 0;
        return getToeplitzBlockIndex(mLevels.getSize(),&mLevels[0],I);
    }
    
private:
    tUIndex getIndex(const tUIndex& i,
                     const tUIndex& j,
                     const tUSInt& level,
                     const tUIndex& nRows,
                     tUIndex& nValues) const;
protected:
    
    inline tUIndex getBlockIndex(const tUIndex& ib,const tUIndex& jb,const tUIndex& N) const {
        //\f$ {T_{-N+1},...,T_{-1},T_0,T_1,.....T_{N-1} \f$
        return ib-jb+N-1;
      
    }
  

    
    
    
public:
     
    //product in spectral space methods
    //==================================
  
    void projection(const tFlag& space) {
        if (space==getSpace()) return;
    
        switch(space) {
        case SPECTRAL_SPACE:
            diagonalize();
            break;
        case CANONICAL_SPACE:
            recover();
            break;
        }
        setSpace(space);
    }
    
protected:
    virtual void diagonalize()=0;
    
    
    virtual void recover()=0;
     
public:
    virtual void vectorProduct(const CORE_RealArray& X,
                               CORE_RealArray& Y) const {
        Y.setSize(X.getSize());
        vectorProduct(X.getSize(),&X[0],&Y[0]);
    }
    virtual void vectorProduct(const tUIndex& N,
                               const tReal* X,
                               tReal* Y) const=0;
    
    virtual tString toString() const; 

    
};
#endif