SM_PackedBlockGeneralMatrix.hpp

#ifndef SM_PackedBlockGeneralMatrix_HPP
#define SM_PackedBlockGeneralMatrix_HPP
 
template<typename T,tUCInt P> 
void SM_PackedBlockGeneralMatrix<T,P>::vectorProduct(const T* vX,T* vY) const {
 
    //P: dimension of a particle
    //number of row-blocks
    const tInteger &N=this->getRowBlocksNumber();
    
    //number of column blocks
    const tInteger &M=this->getColumnBlocksNumber();
 
    //Y=0
    memset(vY,0,sizeof(T)*N*P);
    
    //iterator on X  
    const T* Xj=null;
    
    
    //indices of block
    tIndex ib=0,jb=0;
    
    
    //iterator on block
    const T* vB=&this->getBlocksValues()[0];
    
    //iterator on block
    const T* Bij=vB;
 
    //iterator on Y
    T *Yi=vY;
    T* eYi=Yi;eYi+=P;
 
    //work pointers
    T* w1=null,*w2=null,*w3=null,*w4=null;
    
    for (ib=0;ib<N;ib++) {//row index of the block
 
        Xj=vX;
        for (jb=0;jb<M;jb++) {//column index of the block
            SM_PackedBlockSymmetricMatrix<T,P>::BlockVectorProduct(Bij,Xj,Yi,eYi,w1,w2,w3,w4);
            
            //std::cout<<"M("<<functions_array::toString(SM_PackedBlockSymmetricMatrix<T,P>::BLOCK_SIZE,Bij)<<")."<<functions_array::toString(3,Xj)<<"="<<functions_array::toString(3,Yi)<<"\n";
            //iterators at block at next column
            Bij+=SM_PackedBlockSymmetricMatrix<T,P>::BLOCK_SIZE;
            Xj+=P;
        }
 
        //iterators at block at next row
        Yi=eYi;
        eYi+=P;
    }
   
}
    
 
 
#endif