|| || || || || ||

---

ComputeKernelMatrix.h

Go to the documentation of this file.

#ifndef Impala_Core_Training_ComputeKernelMatrix_h
#define Impala_Core_Training_ComputeKernelMatrix_h

#ifdef CUDA
#include <cuda.h>
#include <cuda_runtime.h>
#include <cublas.h>
#endif
#include "Link/Cuda/PrecomputeMatrixConfig.h"
#include "Core/Feature/VirtualFeatureTableFactory.h"
#include "Core/Vector/Chi2Distance.h"
#include "Link/Cuda/Chi2Distance.h"
#include "Core/Vector/HistogramIntersection.h"
#include "Link/Cuda/HistogramIntersection.h"
#include "Util/TimeStats.h"

namespace Impala
{
namespace Core
{
namespace Training
{


#ifndef CUDA_SAFE_CALL
#define CUDA_SAFE_CALL(msg) {}
#endif

template <class FType, class PrecomputeTaskT>
void
ComputeKernelMatrix(PrecomputeTaskT* pt, int slabWidth, bool GPU)
{
    ILOG_VAR(Impala.Core.Training.ComputeKernelMatrix);
    ILOG_INFO("Computing matrix (fsize=" << sizeof(FType) << ")...");

    // ** compute matrix
    typedef Core::Matrix::VirtualMatrix VirtualMatrix;
    typedef Core::Feature::VirtualFeatureTable VirtualFeatureTable;
    typedef Core::Feature::VirtualFeatureTableFactory VFTFactory;
    std::vector<VirtualFeatureTable*> tablesA;
    std::vector<VirtualFeatureTable*> tablesB;

    // check input sizes
    Int64 vectorCountA = -1;
    Int64 vectorCountB = -1;
    int vectorSize = 0;
    ILOG_INFO("Opening 2 x " << pt->NrFeatures() << " feature tables");
    for (unsigned int feature = 0; feature < pt->NrFeatures(); feature++)
    {
        VirtualFeatureTable* tableA =
            VFTFactory::GetInstance().ConstructIOBufferReader
                                 (pt->GetFeatureLocatorA(feature), true);
        tablesA.push_back(tableA);
        VirtualFeatureTable* tableB =
            VFTFactory::GetInstance().ConstructIOBufferReader
                                 (pt->GetFeatureLocatorB(feature), false);
        tablesB.push_back(tableB);
        vectorSize = Impala::Max(vectorSize, tableA->GetVectorLength());
        if (vectorCountA == -1)
            vectorCountA = tableA->Size();
        if (vectorCountB == -1)
            vectorCountB = tableB->Size();
        if (vectorCountA != tableA->Size())
        {
            ILOG_ERROR("FATAL: Vector count mismatch");
            return;
        }
        if (vectorCountB != tableB->Size())
        {
            ILOG_ERROR("FATAL: Vector count mismatch");
            return;
        }
        ILOG_INFO("Did 2 x " << feature+1 << " feature(s): " <<
                  Process::MemoryInfo::GetUsageString());
    }
    ILOG_INFO("vectorCountA(devel)=" << vectorCountA <<
              ", vectorCountB(test)=" << vectorCountB);
    VirtualMatrix* matrix = pt->GetWritableMatrix(vectorCountB, vectorCountA);
    int alignFactor = (pt->IsChi2()) ? VECTORCHUNK_SIZE : 32;
    int VECTOR_SIZE = IntAlignUp(vectorSize, alignFactor);
    int vectorChunkCount = VECTOR_SIZE / VECTORCHUNK_SIZE;
    ILOG_INFO("VECTOR_SIZE=" << VECTOR_SIZE << "; vectorSize=" << vectorSize <<
              "; chunkCount=" << vectorChunkCount <<
              "; chunkSize=" << VECTORCHUNK_SIZE);
    if (VECTOR_SIZE > 16384)
    {
        slabWidth = Impala::Min(slabWidth, 512);
    }

    // ** allocs start

    // allocate host memory for matrices A and B
    unsigned int size_A = VECTOR_SIZE * slabWidth;
    unsigned int mem_size_A = sizeof(FType) * size_A;
    unsigned int size_B = VECTOR_SIZE * slabWidth;
    unsigned int mem_size_B = sizeof(FType) * size_B;
    FType* h_A = 0; //(FType*) malloc(mem_size_A);
    FType* h_B = 0; //(FType*) malloc(mem_size_B);
#ifdef CUDA
    cudaHostAlloc((void**)&h_A, mem_size_A, cudaHostAllocDefault);
    cudaHostAlloc((void**)&h_B, mem_size_B, cudaHostAllocDefault);
#else
    h_A = (FType*) malloc(mem_size_A);
    h_B = (FType*) malloc(mem_size_B);
#endif

    unsigned int size_C = slabWidth * slabWidth;
    unsigned int mem_size_C = sizeof(FType) * size_C;

    // allocate host memory for the intermediate result
    FType* h_C = 0;
    if (!GPU)
        h_C = (FType*) malloc(mem_size_C);
    // allocate host memory for the final result
    FType* h_D = (FType*) malloc(mem_size_C);

    // allocate device memory
    FType* d_A = 0;
    FType* d_B = 0;
    FType* d_Bextra = 0;
    FType* d_C = 0;
    FType* d_D = 0;
    if (GPU)
    {
        CUDA_SAFE_CALL(cudaMalloc((void**) &d_A, mem_size_A));
        CUDA_SAFE_CALL(cudaMalloc((void**) &d_B, mem_size_B));
        CUDA_SAFE_CALL(cudaMalloc((void**) &d_Bextra, mem_size_B));
        // allocate device memory for intermediate result
        CUDA_SAFE_CALL(cudaMalloc((void**) &d_C, mem_size_C));
        // allocate device memory for final result
        CUDA_SAFE_CALL(cudaMalloc((void**) &d_D, mem_size_C));
    }

    // ** allocs end

    Util::TimeStats statsSlab;
    statsSlab.AddGroup("read");
    statsSlab.AddGroup("process");
    Util::TimeStats statsOverall;
    statsOverall.AddGroupsFromSub(&statsSlab);
    statsOverall.AddGroup("copy");
    statsOverall.AddGroup("write");

    for (int slabA = 0; slabA < vectorCountA; slabA += slabWidth)
    {
        int slabB = 0;
        if (pt->IsSymmetric())
        {
            slabB = slabA;
        }
        for ( ; slabB < vectorCountB; slabB += slabWidth)
        {
            statsOverall.MeasureFirst();
            // initialize final result for slab to zero
            if (GPU)
            {
                CUDA_SAFE_CALL(cudaMemset(d_D, 0, mem_size_C));
            }
            else
            {
                for (int i = 0; i < size_C; i++)
                    h_D[i] = 0.0;
            }
            int slabSizeA = Impala::Min((Int64)slabWidth, vectorCountA - slabA);
            int slabSizeB = Impala::Min((Int64)slabWidth, vectorCountB - slabB);
            
            for (unsigned int feature = 0; feature < pt->NrFeatures(); feature++)
            {
                // read data
                statsSlab.MeasureFirst();
                tablesA[feature]->GetAlignedVectors(slabA, slabWidth, h_A, VECTOR_SIZE);
                tablesB[feature]->GetAlignedVectors(slabB, slabWidth, h_B, VECTOR_SIZE);
                statsSlab.MeasureNext();
                
                if (GPU)
                {
#ifdef CUDA
                    // copy host memory to device
                    CUDA_SAFE_CALL(cudaMemcpy(d_A, h_A, mem_size_A,
                                              cudaMemcpyHostToDevice) );
                    CUDA_SAFE_CALL(cudaMemcpy(d_B, h_B, mem_size_B,
                                              cudaMemcpyHostToDevice) );
                
                    // execute the kernel
                    if (sizeof(FType) == 8)
                    {
                        ILOG_ERROR("Cannot have doubles in GPU mode!");
                    }
                    if (pt->IsChi2())
                        //DoSlabChi2DistanceSlow(VECTORCHUNK_SIZE, slabSizeA, slabSizeB, d_C, d_A, d_B, slabSizeA, vectorChunkCount);
                        Link::Cuda::DoSlabChi2Distance(d_C, d_A, d_B, d_Bextra, VECTOR_SIZE, slabSizeA, slabSizeB, slabWidth);
                    else
                        Link::Cuda::DoSlabHistogramIntersection(d_C, d_A, d_B, d_Bextra, VECTOR_SIZE, slabSizeA, slabSizeB, slabWidth);

                    // task: D = (-weight/avg) * C + D
                    //       Y =  alpha        * X + Y
                    double alpha = pt->GetFeatureWeight(feature);
                    if (pt->IsChi2())
                        alpha = -alpha / pt->GetFeatureAverage(feature);
                    cublasSaxpy(slabWidth * slabWidth, float(alpha), d_C, 1, d_D, 1);

                    // check if kernel execution generated an error
                    CUT_CHECK_ERROR("Kernel execution failed");
#endif // CUDA
                }
                else
                {
                    if (pt->IsChi2())
                        Core::Vector::DoSlabChi2Distance(h_C, h_A, h_B, slabSizeA, slabSizeB, VECTOR_SIZE);
                    else
                        Core::Vector::DoSlabHistogramIntersection(h_C, h_A, h_B, slabSizeA, slabSizeB, VECTOR_SIZE);

                    // task: D = (-weight/avg) * C + D
                    //       Y =  alpha        * X + Y
                    double alpha = pt->GetFeatureWeight(feature);
                    if (pt->IsChi2())
                        alpha = -alpha / pt->GetFeatureAverage(feature);
                    for (int i = 0; i < slabWidth * slabWidth; i++)
                        h_D[i] += (FType)(alpha * h_C[i]);
                }
                statsSlab.MeasureLast();
            }
            ILOG_INFO("slab(" << slabA << "," << slabB << ") " <<
                      statsSlab.AsString());
            statsOverall.MeasureFromSub(&statsSlab);

            if (GPU)
            {
                // copy result from device to host
                CUDA_SAFE_CALL(cudaMemcpy(h_D, d_D, mem_size_C,
                                          cudaMemcpyDeviceToHost) );
            }
            statsOverall.MeasureNext();

            double divisor = pt->GetTotalFeatureWeight();
            for (Int64 b=0 ; b<slabSizeB ; b++)
            {
                FType* row = h_D + b * slabSizeA;
                if (pt->IsChi2())
                {
                    for (Int64 a=0 ; a<slabSizeA ; a++)
                        row[a] = exp(row[a] / divisor);
                }
                else
                {
                    for (Int64 a=0 ; a<slabSizeA ; a++)
                        row[a] = row[a] / divisor;
                }
                matrix->AddRowPart(slabB+b, slabA, row, slabSizeA);
            }

            if ((slabA != slabB) && pt->IsSymmetric())
            {
                FType* tmpBuf = new FType[slabSizeB];
                // do the symmmetric write
                for (Int64 a = 0; a < slabSizeA; a++)
                {
                    for (Int64 b = 0; b < slabSizeB; b++)
                        tmpBuf[b] = h_D[b * slabSizeA + a];
                    matrix->AddRowPart(slabA+a, slabB, tmpBuf, slabSizeB);
                }
                delete tmpBuf;
            }
            statsOverall.MeasureLast();
            ILOG_INFO("Overall " << statsOverall.AsString());
        }
    }
    pt->Finalize();
    // ** compute matrix ends
    
    // clean up memory
#ifdef CUDA
    cudaFreeHost(h_A);
    cudaFreeHost(h_B);
#else
    free(h_A);
    free(h_B);
#endif
    if (!GPU)
        free(h_C);
    free(h_D);
    if (GPU)
    {
        CUDA_SAFE_CALL(cudaFree(d_A));
        CUDA_SAFE_CALL(cudaFree(d_B));
        CUDA_SAFE_CALL(cudaFree(d_Bextra));
        CUDA_SAFE_CALL(cudaFree(d_C));
        CUDA_SAFE_CALL(cudaFree(d_D));
    }
    for (unsigned int feature=0 ; feature<pt->NrFeatures() ; feature++)
    {
        delete tablesA[feature];
        delete tablesB[feature];
    }
    tablesA.clear();
    tablesB.clear();
}

} // namespace Training
} // namespace Core
} // namespace Impala

#endif

---