SM_PackedBlockSymmetricMatrix.hpp

#ifndef SM_PackedBlockSymmetricMatrix_HPP
#define SM_PackedBlockSymmetricMatrix_HPP
 
template<typename T,tUCInt P> 
void SM_PackedBlockSymmetricMatrix<T,P>::vectorProduct(const T* vX,T* vY) const {
 
    const tIndex &N=this->getRowBlocksNumber();
    
    //vY=0
    memset(vY,0,sizeof(T)*P*N);
  
  
    //diagonla values
    const T* vD=&this->getDiagonalValues()[0];
    
    //iterator on block
    const T* vB=&this->getSupBlocksValues()[0];
    
    
    
    
    //work indices of element in block of size PxP
    //tUCInt i=0,j=0;
    //tUCInt incB=0;//increment on block
    //const T* iXj=null;//iterator of Xj inside block
    
    //iterator on X  
    const T* Xj=0;
    
    
    //indices of block
    tIndex ib=0,jb=0;
    
    //iterator on diagonal values
    const T* Di=vD;
    
    //iterator on block
    const T* Bij=vB;
    
    //block at row 0 and column ib
    const T* B0i=Bij;
    
    //increment of block for next column
    tIndex incBij=0;
    
    
    //Y values at block start
    T* Yi_d,*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
        
        
        //Yi+=Bij.Xj
        Xj=vX;
        
        //Yi_x=&Yi[0];
        //Yi_y=&Yi[1];
        //Yi_z=&Yi[2];
        
        Bij=B0i;
        for (jb=0;jb<ib;jb++) {//column index of the block
            
            //Block product Yi + =B_ij . X_j
            //BlockVectorProduct(Bij,Xj,Yi,
            //                   i,j,incB,iXj);
            //Xj_x=Xj[0];
            //Xj_y=Xj[1];
            //Xj_z=Xj[2];
            //BlockVectorProduct(Bij,Xj_x,Xj_y,Xj_z,*Yi_x,*Yi_y,*Yi_z);
            
            SM_PackedBlockSymmetricMatrix<T,P>::BlockVectorProduct(Bij,Xj,Yi,eYi,w1,w2,w3,w4);
            
            //iterators at block at next row
            Bij+=SM_PackedBlockSymmetricMatrix<T,P>::BLOCK_SIZE;
            
            Xj+=P;
            
        }//end loop on column index < i index
        
        
        //diagonal term
        //jb=ib Yi+=Di.Xi
        Yi_d=Yi;
        while (Yi_d!=eYi) {
            (*Yi_d)+=(*Di)*(*Xj);
            Yi_d++;
            Xj++;
        }
        Di++;
        
        //jb>ib : jb=ib+1;
        jb=ib;jb++;
        //block for next column
        B0i=Bij;
        //jb-th block of next colum 
        incBij=ib;
        incBij*=SM_PackedBlockSymmetricMatrix<T,P>::BLOCK_SIZE;
        Bij+=incBij;
        while (jb<N) {//column index of the block
            
            //Block product
            
            SM_PackedBlockSymmetricMatrix<T,P>::BlockVectorProduct(Bij,Xj,Yi,eYi,w1,w2,w3,w4);
            
            
            //next block for X
            Xj+=P;
            
            //iterators to block at next column +=jb*BLOCK_SIZE
            incBij+=SM_PackedBlockSymmetricMatrix<T,P>::BLOCK_SIZE;
            Bij+=incBij;
            
            //next column
            jb++;
        }//end loop on column index >= i index
        
        //Yi at next row block
        Yi=eYi;
        eYi+=P;
        
    }//end loop on row index
 
}
    
 
 
#endif