SMOMPI_PackedBlockSymmetricMatrix.hpp

#ifndef SMOMPI_PackedBlockSymmetricMatrix_HPP
#define SMOMPI_PackedBlockSymmetricMatrix_HPP
 
#include "MPI_Run.h"
#include "MPI_WorldToWorldMessage.h"
 
template<typename T,tUCInt P> 
void SMOMPI_PackedBlockSymmetricMatrix<T,P>::vectorProductBbB(const T* vX,T* vY ) const {
 
    const tIndex &N=this->getRowBlocksNumber();
    
    //vY=0
    memset(vY,0,sizeof(T)*P*N);
    
    //MPI environment
    const MPI_Environment& mpiEnv=MPI_Run::GetEnvironment();
    
    //current id of the MPI core
    const tMPICoreId& coreId=mpiEnv.getCoreId();
    
    //number of MPI cores
    const tMPICoreId& coresNumber=mpiEnv.getCoresNumber();
 
    std::valarray<int> valuesLengthPerCore;
    valuesLengthPerCore.resize(coresNumber);
 
    std::valarray<int> valuesOffset;
    valuesOffset.resize(coresNumber);
        
    //start index for diagonal values
    tIndex start=coreId*N/coresNumber;
    
    //end index for diagonal values
    tIndex end=(coreId+1)*N/coresNumber;
    
    //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[start];
    
    //iterator on block
    const T* Bij=&vB[(start==0)?0:SM_PackedBlockSymmetricMatrix<T,P>::GetBlockIndex(0,start)];
    
    //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[start*P];
    T* eYi=Yi;eYi+=P;
    
    //work pointers
    T* w1,*w2,*w3,*w4;
    
    for (ib=start;ib<end;ib++) {//row index of the block in [start,end[
        
        
        //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
 
    //dispatch all values
    MPI_WorldToWorldMessage::AllGather(mpiEnv,(int)(end-start)*P,valuesLengthPerCore);
    MPI_WorldToWorldMessage::BuildOffset(valuesLengthPerCore,valuesOffset);
    
    MPI_WorldToWorldMessage::AllGather(mpiEnv,&vY[start*P],valuesLengthPerCore[coreId],
                                       &vY[0],valuesLengthPerCore,valuesOffset);
    
}
 
 
 
template<typename T,tUCInt P>
void SMOMPI_PackedBlockSymmetricMatrix<T,P>::vectorProduct3D(const T* vX,T* vY ) const {
 
    const tIndex &N=this->getRowBlocksNumber();
    
    //vY=0
    memset(vY,0,sizeof(T)*P*N);
 
   
    //MPI environment
    const MPI_Environment& mpiEnv=MPI_Run::GetEnvironment();
    
    //current id of the MPI core
    const tMPICoreId& coreId=mpiEnv.getCoreId();
    
    //number of MPI cores
    const tMPICoreId& coresNumber=mpiEnv.getCoresNumber();
 
    std::valarray<int> valuesLengthPerCore;
    valuesLengthPerCore.resize(coresNumber);
    std::valarray<int> valuesOffset;
    valuesOffset.resize(coresNumber);
        
    //start index for diagonal values
    tIndex start=coreId*N/coresNumber;
    
    //end index for diagonal values
    tIndex end=(coreId+1)*N/coresNumber;
    
    //diagonal 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[start];
    
    //iterator on block
    const T* Bij=&vB[(start==0)?0:SM_PackedBlockSymmetricMatrix<T,P>::GetBlockIndex(0,start)];
    
    //block at row 0 and column ib
    const T* B0i=Bij;
    
    //increment of block for next column
    tIndex incBij=0;
    
    //iterator on Yi at coordinate d
    //T* Yi_d=null;
    
    //Y values at block start
    T* Yi=&vY[start*P];
    
    tReal Xj_x,Xj_y,Xj_z;
    tReal *Yi_x,*Yi_y,*Yi_z;
    for (ib=start;ib<end;ib++) {//row index of the block in [start,end[
        
        
        //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];
            
            //Bij[0]
            (*Yi_x)+=(*Bij)*Xj_x;
            Bij++;
            
            //Bij[1]
            (*Yi_x)+=(*Bij)*Xj_y;
            (*Yi_y)+=(*Bij)*Xj_x;
            Bij++;
            
            //Bij[2]
            (*Yi_y)+=(*Bij)*Xj_y;
            Bij++;
            
            //Bij[3]
            (*Yi_x)+=(*Bij)*Xj_z;
            (*Yi_z)+=(*Bij)*Xj_x;
            Bij++;
            
            //Bij[4]
            (*Yi_y)+=(*Bij)*Xj_z;
            (*Yi_z)+=(*Bij)*Xj_y;
            Bij++;
            
            //Bij[5]
            (*Yi_z)+=(*Bij)*Xj_z;
            Bij++;
            
            //iterators at block at next row
            //Bij+=BLOCK_SIZE;
            
            Xj+=P;
            
        }//end loop on column index < i index
        
        //jb=ib Yi+=Di.Xi
        // Yi_d=Yi;
        // for(jb=0;jb<P;jb++) {
        //     (*Yi_d)+=(*Di)*(*Xj);
        //     Xj++;
        //     Yi_d++;
        // }
        Xj_x=Xj[0];
        Xj_y=Xj[1];
        Xj_z=Xj[2];
        (*Yi_x)+=(*Di)*Xj_x;
        (*Yi_y)+=(*Di)*Xj_y;
        (*Yi_z)+=(*Di)*Xj_z;
        Xj+=P;
        // if (ib==2) {
        //     std::cout<<"<"<<ib<<","<<jb<<">["<<(*Di)<<"]";
        //     std::cout<<"["<<Xj[-3]<<",";
        //     std::cout<<Xj[-2]<<",";
        //     std::cout<<Xj[-1]<<"]";
        //     std::cout<<"=["<<Yi[0]<<","<<Yi[1]<<","<<Yi[2]<<"]\n";
                    
        // }
        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
            //BlockVectorProduct(Bij,Xj,Yi,
            //                   i,j,incB,iXj);
            Xj_x=Xj[0];
            Xj_y=Xj[1];
            Xj_z=Xj[2];
            
            //Bij[0]
            (*Yi_x)+=(*Bij)*Xj_x;
            Bij++;
            
            //Bij[1]
            (*Yi_x)+=(*Bij)*Xj_y;
            (*Yi_y)+=(*Bij)*Xj_x;
            Bij++;
            
            //Bij[2]
            (*Yi_y)+=(*Bij)*Xj_y;
            Bij++;
            
            //Bij[3]
            (*Yi_x)+=(*Bij)*Xj_z;
            (*Yi_z)+=(*Bij)*Xj_x;
            Bij++;
            
            //Bij[4]
            (*Yi_y)+=(*Bij)*Xj_z;
            (*Yi_z)+=(*Bij)*Xj_y;
            Bij++;
            
            //Bij[5]
            (*Yi_z)+=(*Bij)*Xj_z;
            Bij++;
            
            
            // if (ib==2) {
            //     std::cout<<"<"<<ib<<","<<jb<<">["<<Bij[0]<<",";
            //     std::cout<<Bij[1]<<",";
            //     std::cout<<Bij[2]<<",";
            //     std::cout<<Bij[3]<<",";
            //     std::cout<<Bij[4]<<",";
            //     std::cout<<Bij[5]<<"].";
            //     std::cout<<"["<<Xj[0]<<",";
            //     std::cout<<Xj[1]<<",";
            //     std::cout<<Xj[2]<<"]";
            //     std::cout<<"=["<<Yi[0]<<","<<Yi[1]<<","<<Yi[2]<<"]\n";
                    
            // }
            //next block for X
            Xj+=P;
            //iterators to block at next column +=jb*BLOCK_SIZE
            Bij+=incBij;
            incBij+=SM_PackedBlockSymmetricMatrix<T,P>::BLOCK_SIZE;
            //next column
            jb++;
        }//end loop on column index >= i index
            
        //Yi at next row block
        Yi+=P;
            
    }//end loop on row index
    
    //dispatch all values
    MPI_WorldToWorldMessage::AllGather(mpiEnv,(int) (end-start)*P,valuesLengthPerCore);
    MPI_WorldToWorldMessage::BuildOffset(valuesLengthPerCore,valuesOffset);
    
    MPI_WorldToWorldMessage::AllGather(mpiEnv,&vY[start*P],valuesLengthPerCore[coreId],
                                       &vY[0],valuesLengthPerCore,valuesOffset);
}
 
 
 
 
#endif