SM_PackedBlockSymmetricMatrix.h

#ifndef SM_PackedBlockSymmetricMatrix_H
#define SM_PackedBlockSymmetricMatrix_H
 
//base classes
#include "SM_PackedBlockMatrix.h"
 
 
template<typename T,tUCInt P>
class SM_PackedBlockSymmetricMatrix : public SM_PackedBlockMatrix<T,P> {
 
public:
 
 
private :
 
    //type class
    typedef SM_PackedBlockMatrix<T,P> SuperClass;
    typedef SM_PackedBlockSymmetricMatrix<T,P> SelfClass;
    
    //single attributes
    //-----------------
    
 
 
    
    //multiple attributes
    //------------------
 
    //values of superior blocks of the matrix
    std::valarray<T> mSupBlocksValues;
    //diagonal values : Diagonal[i].I_P
    std::valarray<T> mDiagonalValues;
 
 
    tIndex mWIndex;
    
    //single associations
    //-------------------
 
 
    //multiple associations
    //---------------------
 
    
   
protected:
    // CONSTRUCTORS 
    SM_PackedBlockSymmetricMatrix(void) {
      
    }
    
    // DESTRUCTORS 
    virtual ~SM_PackedBlockSymmetricMatrix(void) {
    }
   
    
public : 
 
    //Instance building
    //=================
 
   
    virtual tMemSize getMemorySize() const {
        return sizeof(*this)+this->getContentsMemorySize();
    }
 
    virtual tMemSize getContentsMemorySize() const {
        tMemSize mem=SuperClass::getContentsMemorySize();
        mem+=mSupBlocksValues.size()*sizeof(T);
        mem+=mDiagonalValues.size()*sizeof(T);
        return mem;
    }
 
    inline static CORE_UniquePointer<SelfClass> New() {
        return CORE_UniquePointer<SelfClass>(new SelfClass(),
                                             SM_Object::Delete());
    }
 
public:
   
    //SET & GET
    //=========
 
    virtual void reset() final {
        memset(&mSupBlocksValues[0],0,mSupBlocksValues.size()*sizeof(T));
        memset(&mDiagonalValues[0],0,mDiagonalValues.size()*sizeof(T));
        
    }
    inline void setSize(const tIndex& nRowBlocks)  {
        setSize(nRowBlocks,nRowBlocks);
    }
    
    virtual void setSize(const tIndex& nRowBlocks,const tIndex& nColBlocks) final {
        ASSERT_IN(nRowBlocks==nColBlocks);
        
        SuperClass::setSize(nRowBlocks,nColBlocks);
 
        //number of elements of the sup block values
        tIndex s=GetSupBlocksValuesNumber(nRowBlocks);
        if (s!=mSupBlocksValues.size()) {
            try {
                mSupBlocksValues.resize(s);
                memset(&mSupBlocksValues[0],0,sizeof(T)*s);
            } catch (std::exception& e) {
                throw CORE_Exception("stochmagnet/operators/dipolar",
                                     "SM_PackedBlockSymmetricMatrix::setSize("+std::to_string(nRowBlocks)+")",
                                     "not enough size to allocate "+std::to_string(s)+" elements for sup matrix : memory size needed:  "+std::to_string((sizeof(T)*s)/(1024*1024*1024))+" Go");
            }
        }
        if (nRowBlocks!=mDiagonalValues.size()) {
            try {
                mDiagonalValues.resize(nRowBlocks);
                memset(&mDiagonalValues[0],0,sizeof(T)*nRowBlocks);
            } catch (std::exception& e) {
                throw CORE_Exception("stochmagnet/operators/dipolar",
                                     "SM_PackedBlockSymmetricMatrix::setSize("+std::to_string(nRowBlocks)+")",
                                     "not enough size to allocate "+std::to_string(s)+" elements for diagonla matrix : memory size needed:  "+std::to_string((sizeof(T)*s)/(1024*1024*1024))+" Go");
                
            }
        }
    }
 
    virtual void setValue(const tIndex& i,const tIndex& j,const T& v) final{
        (*this)(i,j)=v;
    }
   
    virtual  const T& getValue(const tIndex& i,const tIndex& j) const final {
        return (*this)(i,j);
    }
    virtual  T& getValue(const tIndex& i,const tIndex& j)  final {
        return (*this)(i,j);
    }
    
public:
    
    inline std::valarray<T>& getSupBlocksValues() {
        return mSupBlocksValues;
    }
    inline std::valarray<T>& getDiagonalValues(){
        return mDiagonalValues;
    }
 
    inline const std::valarray<T>& getSupBlocksValues() const {
        return mSupBlocksValues;
    }
    inline const std::valarray<T>& getDiagonalValues() const {
        return mDiagonalValues;
    }
    
    
    //OPERATORS
    //==========
public:
    inline const T& operator()(tIndex i,tIndex j) const {
        tIndex  index=this->GetPackedIndex(i%P,j%P);
        //index of the block
        i/=P;
        j/=P;
        return (i==j)?mDiagonalValues[i]:mSupBlocksValues[GetSupBlockIndex<P>(i,j)+index];
    }
    inline T& operator()(tIndex i,tIndex j) {
        mWIndex=this->GetPackedIndex(i%P,j%P);
        i/=P;
        j/=P;
        return (i==j)?mDiagonalValues[i]:mSupBlocksValues[GetSupBlockIndex<P>(i,j)+mWIndex];
    }
  
    
private:
    inline  static tIndex GetSupBlocksValuesNumber(const tIndex& n) {
        //number of blocks
        tIndex p=n;
        p*=(n-1);
        p/=2;
        //number of elements
        p*= SM_PackedBlockMatrix<T,P>::BLOCK_SIZE;
        return p;
    }
   
   
public:
    inline static tIndex GetBlockIndex(const tIndex& ib ,const tIndex& jb) {
        return GetSupBlockIndex<SM_PackedBlockMatrix<T,P>::BLOCK_SIZE>(ib,jb);
    }
private:
    template<tUCInt s>
    inline static tIndex GetBlockIndex(const tIndex& i ,const tIndex& j) {
        ASSERT_IN(i<=j);
        tIndex index=j;
        index++;
        index*=j;
        index/=2;
        index+=i;
        index*=s;
        return index;
    }
    template<tUCInt s>
    inline static tIndex GetSupBlockIndex(const tIndex& ib ,const tIndex& jb) {
        ASSERT_IN(ib<jb);
        tIndex index=jb;
        index--;
        index*=jb;
        index/=2;
        index+=ib;
        index*=s;
        return index;
    }
    
 
public:
  
    virtual void vectorProduct(const T* vX,T* vY) const;
    
public:
    
    virtual tString toString() const override {
        std::stringstream ret;
        ret<<SuperClass::toString()<<"\n";
        //print matrix 6 first blocks
        ret<<"sup blocks values:"<<functions_array::toString(mSupBlocksValues)<<"\n";
        ret<<"diagonal values:"<<functions_array::toString(mDiagonalValues)<<"\n";
        return ret.str();
    }
  
};
 
#include "SM_PackedBlockSymmetricMatrix.hpp"
 
#endif