KNMatrixOperation.cpp

Go to the documentation of this file.

#include "KNMatrixOperation.h"
#include "KNMPIManager.h"
#include "KNTimeMeasurement.h"
#include "KNIPCCUtility.h"

#include "CKNGlobal.h"

#include "XeonPhi_header.h"

#define LOOP_OPTIMIZE_COUNT     10

unsigned int CKNMatrixOperation::CKNCSR::MAX_INDEX = 0xffffffff;
CKNComplex* CKNMatrixOperation::pValueBuffer = NULL;
unsigned int* CKNMatrixOperation::pRow = NULL;
unsigned int* CKNMatrixOperation::pColumn = NULL;

CKNMatrixOperation::CKNVector::CKNVector()
{
    m_nValueCount = 0;
}

CKNMatrixOperation::CKNVector::~CKNVector()
{
}

void CKNMatrixOperation::CKNVector::SetSize(unsigned int nSize)
{
    if (nSize == m_nValueCount)
        return;

    m_nValueCount = nSize;
    m_vectValueRealBuffer.resize(nSize);
    m_vectValueImaginaryBuffer.resize(nSize);
}

void CKNMatrixOperation::CKNVector::SetAt(unsigned int nIndex, CKNComplex value)
{
    SetAt(nIndex, value.GetRealNumber(), value.GetImaginaryNumber());
}

void CKNMatrixOperation::CKNVector::SetAtEx(unsigned int nIndex, CKNComplex *pValue)
{
    SetAt(nIndex, pValue->GetRealNumber(), pValue->GetImaginaryNumber());
}

void CKNMatrixOperation::CKNVector::SetAt(unsigned int nIndex, double fReal, double fImaginary)
{
    if (nIndex > GetSize())
    {
        THROW_END_EXIT(ERROR_OUT_OF_RANGE);
    }

    m_vectValueRealBuffer[nIndex] = fReal;
    m_vectValueImaginaryBuffer[nIndex] = fImaginary;
}

CKNComplex* CKNMatrixOperation::CKNVector::GetAtPt(unsigned int nIndex)
{
    if (nIndex > GetSize())
    {
        throw ERROR_OUT_OF_RANGE;
        return NULL;
    }

    m_rtnTemp.SetComplexNumber(m_vectValueRealBuffer[nIndex], m_vectValueImaginaryBuffer[nIndex]);
    return &m_rtnTemp;
}

CKNComplex CKNMatrixOperation::CKNVector::GetAt(unsigned int nIndex)
{
    CKNComplex      rtnComplex;

    if (nIndex > GetSize())
    {
        throw ERROR_OUT_OF_RANGE;
        return rtnComplex;
    }

    rtnComplex.SetComplexNumber(m_vectValueRealBuffer[nIndex], m_vectValueImaginaryBuffer[nIndex]);
    return rtnComplex;
}

void CKNMatrixOperation::CKNVector::ResetValue()
{
    unsigned int        i;

    for (i = 0; i < GetSize(); i++)
    {
        m_vectValueRealBuffer[i] = 0.;
        m_vectValueImaginaryBuffer[i] = 0.;
    }
    
}

void CKNMatrixOperation::CKNVector::ScalarMultiple(CKNComplex Scalar)
{
    unsigned int            i;
    CKNComplex              tempComplex;

    for (i = 0; i < GetSize(); i++)
    {
        tempComplex.SetComplexNumber(m_vectValueRealBuffer[i], m_vectValueImaginaryBuffer[i]);
        tempComplex = tempComplex * Scalar;

        m_vectValueRealBuffer[i] = tempComplex.GetRealNumber();
        m_vectValueImaginaryBuffer[i] = tempComplex.GetImaginaryNumber();
    }

    return;
}

void CKNMatrixOperation::CKNVector::ScalarMultiple(double fScalar)
{
    //double                    fReal, fImaginary;
    unsigned int            i;

    for (i = 0; i < GetSize(); i++)
    {
        m_vectValueRealBuffer[i] *= fScalar;
        m_vectValueImaginaryBuffer[i] *= fScalar;
    }
}

void CKNMatrixOperation::CKNVector::ScalarDivision(CKNComplex Scalar)
{
    unsigned int            i;
    CKNComplex              tempComplex;

    for (i = 0; i < GetSize(); i++)
    {
        tempComplex.SetComplexNumber(m_vectValueRealBuffer[i], m_vectValueImaginaryBuffer[i]);
        tempComplex = tempComplex / Scalar;

        m_vectValueRealBuffer[i] = tempComplex.GetRealNumber();
        m_vectValueImaginaryBuffer[i] = tempComplex.GetImaginaryNumber();
    }
}

void CKNMatrixOperation::CKNVector::ScalarDivision(double fScalar)
{
    double                  *pReal = NULL, *pImaginary = NULL;
    unsigned int            i;
    unsigned int            nSize = GetSize();

    pReal = m_vectValueRealBuffer.data();
    pImaginary = m_vectValueImaginaryBuffer.data();

#pragma omp parallel for
    for (i = 0; i < nSize ; i++)
    {
        /*m_vectValueRealBuffer[i] /= fScalar;
        m_vectValueImaginaryBuffer[i] /= fScalar;*/
        pReal[i] /= fScalar;
        pImaginary[i] /= fScalar;
    }
}

void CKNMatrixOperation::CKNVector::BuildRandomVector()
{
    unsigned int            i;

    srand((unsigned int)time(NULL));
    for (i = 0; i < m_nValueCount; i++)
    {
        m_vectValueRealBuffer[i] = rand();
        m_vectValueImaginaryBuffer[i] = rand();
    }
}

void CKNMatrixOperation::CKNVector::Normalize(bool bMPI)
{
    double  fVectorSize = GetNorm(bMPI);
    unsigned int        i;

    for (i = 0; i < m_nValueCount; i++)
    {
        m_vectValueRealBuffer[i] /= fVectorSize;
        m_vectValueImaginaryBuffer[i] /= fVectorSize;
    }
}

double CKNMatrixOperation::CKNVector::GetNorm(bool bMPI)
{
    unsigned int    i;
    CKNComplex      PowerTotal;
    double          fTotal = 0.0;
    unsigned int    nLeft = m_nValueCount % LOOP_OPTIMIZE_COUNT;
    double          fNorm;
    double          *pReal = m_vectValueRealBuffer.data();
    double          *pImaginary = m_vectValueImaginaryBuffer.data();

#pragma omp parallel for reduction(+:fTotal)
    for (i = 0; i < m_nValueCount; i++)
    {
        const double        fReal = pReal[i];
        const double        fImaginary = pImaginary[i];

        //fNorm = sqrt(m_vectValueRealBuffer[i] * m_vectValueRealBuffer[i] + m_vectValueImaginaryBuffer[i] * m_vectValueImaginaryBuffer[i]);
        /*fNorm = sqrt(fReal * fReal + fImaginary * fImaginary);
        fTotal += (fNorm * fNorm);*/
        fNorm = fReal * fReal + fImaginary * fImaginary;
        fTotal += fNorm;
    }

    if (bMPI)
        fTotal = CKNMPIManager::AllReduceDouble(fTotal);
    return sqrt(fTotal);
}

void CKNMatrixOperation::CKNVector::ScalarMultiThanMinusVector(double fScalar, CKNVector *vector)
{
    unsigned int        i, nSize = GetSize();
    double              *pReal = NULL, *pImaginary = NULL;
    double              *pOperandReal = NULL, *pOperandImagianray = NULL;

    if (nSize != vector->GetSize())
    {
        throw ERROR_WRONG_ORDER_OPERATION;
        return;
    }

    pReal = m_vectValueRealBuffer.data();
    pImaginary = m_vectValueImaginaryBuffer.data();
    pOperandReal = vector->m_vectValueRealBuffer.data();
    pOperandImagianray = vector->m_vectValueImaginaryBuffer.data();

#pragma omp parallel for
    for (i = 0; i < nSize; i++)
    {
        pReal[i] = pReal[i] - fScalar * pOperandReal[i];
        pImaginary[i] = pImaginary[i] - fScalar * pOperandImagianray[i];
    }
}

void CKNMatrixOperation::CKNVector::ReorthogonalizationVector(CKNVector *pVector, CKNComplex complex)
{
    unsigned int        i, nSize = GetSize();
    double              *pReal = NULL, *pImaginary = NULL;
    double              *pOperandReal = NULL, *pOperandImagianray = NULL;
    double              fReal = complex.GetRealNumber(), fImaginary = complex.GetImaginaryNumber();

    if (nSize != pVector->GetSize())
    {
        throw ERROR_WRONG_ORDER_OPERATION;
        return;
    }

    pReal = m_vectValueRealBuffer.data();
    pImaginary = m_vectValueImaginaryBuffer.data();
    pOperandReal = pVector->m_vectValueRealBuffer.data();
    pOperandImagianray = pVector->m_vectValueImaginaryBuffer.data();

#pragma omp parallel for
    for (i = 0; i < nSize; i++)
    {
        pReal[i] = pReal[i] - fReal* pOperandReal[i] + fImaginary * pOperandImagianray[i];
        pImaginary[i] = pImaginary[i] - fReal * pOperandImagianray[i] - fImaginary * pOperandReal[i];
    }
}

void CKNMatrixOperation::CKNVector::MinusVector(CKNVector *vector)
{
    unsigned int        i, nSize = GetSize();
    CKNComplex          *pOperand1, *pOperand2;

    if (nSize != vector->GetSize())
    {
        throw ERROR_WRONG_ORDER_OPERATION;
        return;
    }

    for (i = 0; i < nSize; i++)
    {
        pOperand1 = GetAtPt(i);
        pOperand2 = vector->GetAtPt(i);

        SetAt(i, pOperand1->GetRealNumber() - pOperand2->GetRealNumber(), pOperand1->GetImaginaryNumber() - pOperand2->GetImaginaryNumber());
    }
}

void CKNMatrixOperation::CKNVector::PlusVector(CKNVector *vector)
{
    unsigned int        i, nSize = GetSize();
    CKNComplex          *pOperand1, *pOperand2;

    if (nSize != vector->GetSize())
    {
        throw ERROR_WRONG_ORDER_OPERATION;
        return;
    }

    for (i = 0; i < nSize; i++)
    {
        pOperand1 = GetAtPt(i);
        pOperand2 = vector->GetAtPt(i);

        SetAt(i, pOperand1->GetRealNumber() + pOperand2->GetRealNumber(), pOperand1->GetImaginaryNumber() + pOperand2->GetImaginaryNumber());
    }
}

void CKNMatrixOperation::CKNVector::Finalize()
{
    m_vectValueRealBuffer.clear();
    m_vectValueImaginaryBuffer.clear();

    m_nValueCount = 0;
}

bool CKNMatrixOperation::CKNVector::InsertVector(unsigned int nStartIndex, CKNMatrixOperation::CKNVector *pVector)
{
    bool                    bRtn = false;
    unsigned int            i;

    if (nStartIndex > m_nValueCount || nStartIndex + pVector->GetSize() > m_nValueCount)
        return bRtn;

    for (i = 0; i < pVector->GetSize(); i++)
    {
        m_vectValueRealBuffer[nStartIndex + i] = pVector->GetAt(i).GetRealNumber();;
        m_vectValueImaginaryBuffer[nStartIndex + i] = pVector->GetAt(i).GetImaginaryNumber();
    }

    bRtn = true;
    return bRtn;
}

bool CKNMatrixOperation::CKNVector::Serialize(double *pBuffer, bool bStore)
{
    double              *pReal = NULL, *pImaginariy = NULL;
    bool                bRtn = false;

    if( NULL == pBuffer )
        return bRtn;

    if( bStore)
    {
        pReal = m_vectValueRealBuffer.data();
        memcpy(pReal, pBuffer, m_nValueCount * sizeof(double));

        pImaginariy = m_vectValueImaginaryBuffer.data();
        memcpy(pImaginariy, pBuffer + m_nValueCount, m_nValueCount * sizeof(double));
    }
    else
    {
        pReal = m_vectValueRealBuffer.data();
        memcpy(pBuffer, pReal, m_nValueCount * sizeof(double));

        pImaginariy = m_vectValueImaginaryBuffer.data();
        memcpy(pBuffer + m_nValueCount, pImaginariy, m_nValueCount * sizeof(double));
    }

    bRtn = true;
    return bRtn;
}

CKNComplex CKNMatrixOperation::CKNVector::operator*(CKNVector * vector)
{
    return operator*((*vector));
}

CKNComplex CKNMatrixOperation::CKNVector::operator*(CKNVector & vector)
{
    CKNComplex Rtn;
    unsigned int        i, nSize = GetSize();

    if (nSize != vector.GetSize())
    {
        throw ERROR_WRONG_ORDER_OPERATION;
        return Rtn;
    }

    for (i = 0; i < nSize; i++)
    {
        Rtn += (GetAt(i) * vector.GetAt(i));
    }
    return Rtn;
}

CKNMatrixOperation::CKNVector CKNMatrixOperation::CKNVector::operator-(CKNVector * vector)
{
    return operator-((*vector));
}

CKNMatrixOperation::CKNVector CKNMatrixOperation::CKNVector::operator-(CKNVector & vector)
{
    CKNVector           rtnVector;
    unsigned int        i, nSize = GetSize();

    if (nSize != vector.GetSize())
    {
        throw ERROR_WRONG_ORDER_OPERATION;
        return rtnVector;
    }

    rtnVector.SetSize(nSize);

    for (i = 0; i < nSize; i++)
    {
        rtnVector.SetAt(i, GetAt(i) - vector.GetAt(i));
    }
    return rtnVector;
}

CKNMatrixOperation::CKNVector CKNMatrixOperation::CKNVector::operator+(CKNVector * vector)
{
    return operator+((*vector));
}

CKNMatrixOperation::CKNVector CKNMatrixOperation::CKNVector::operator+(CKNVector & vector)
{
    CKNVector       rtnVector;
    unsigned int    i, nSize = GetSize();

    if (nSize != vector.GetSize())
    {
        throw ERROR_WRONG_ORDER_OPERATION;
        return rtnVector;
    }

    rtnVector.SetSize(nSize);
    for (i = 0; i < nSize; i++)
    {
        rtnVector.SetAt(i, GetAt(i) + vector.GetAt(i));
    }

    return rtnVector;
}

void CKNMatrixOperation::CKNVector::operator=(CKNVector *vector)
{
    operator=((*vector));
}

void CKNMatrixOperation::CKNVector::operator=(CKNVector &vector)
{
    unsigned int        i, nSize = vector.GetSize();
    double              *pReal = NULL, *pImaginary = NULL;
    double              *pSourceReal = NULL, *pSourceImagianry = NULL;

    SetSize(nSize);

    pReal = m_vectValueRealBuffer.data();
    pImaginary = m_vectValueImaginaryBuffer.data();
    pSourceReal = vector.m_vectValueRealBuffer.data();
    pSourceImagianry = vector.m_vectValueImaginaryBuffer.data();

#pragma omp parallel for
    for (i = 0; i < nSize; i++)
    {
        pReal[i] = pSourceReal[i];
        pImaginary[i] = pSourceImagianry[i];
    }
}


CKNMatrixOperation::CKNDMatrix::CKNDMatrix()
{
}

CKNMatrixOperation::CKNDMatrix::~CKNDMatrix()
{
}

bool CKNMatrixOperation::CKNDMatrix::BuildMatrixFirst(unsigned int nRow, unsigned int nColumn)
{
    bool                bRtn = false;
    unsigned int        i, j;

    if (!m_vectValueBuffer.empty())
        return bRtn;

    bRtn = true;
    for (i = 0; i < nRow; i++)
    {
        for (j = 0; j < nColumn; j++)
        {
            CKNComplex          element;
            m_vectValueBuffer.push_back(element);
        }
    }

    m_nRowCount = nRow;
    m_nColumnCount = nColumn;

    return bRtn;
}

bool CKNMatrixOperation::CKNDMatrix::SetElement(unsigned int nRow, unsigned int nColumn, CKNComplex element)
{
    return SetElement(nRow, nColumn, element.GetRealNumber(), element.GetImaginaryNumber());
}

bool CKNMatrixOperation::CKNDMatrix::SetElement(unsigned int nRow, unsigned int nColumn, double fRealNumber, double fImageNumber)
{
    bool            bRtn = false;

    if (nRow > m_nRowCount || nColumn > m_nColumnCount)
        return bRtn;

    m_vectValueBuffer[m_nColumnCount*nRow + nColumn].SetComplexNumber(fRealNumber, fImageNumber);

    bRtn = true;
    return bRtn;
}

bool CKNMatrixOperation::CKNDMatrix::SetElement(unsigned int nRowStart, unsigned int nColumnStart, unsigned int nSrcRowStart, unsigned int nSrcColumnStart, unsigned int nSrcRowCount, unsigned int nSrcColumnCount, CKNDMatrix matrix)
{
    bool                    bRtn = false;
    unsigned int            i, j;

    if (nRowStart > m_nRowCount || nRowStart + nSrcRowCount > m_nRowCount
        || nColumnStart > m_nColumnCount || nColumnStart + nSrcColumnCount > m_nColumnCount)
        return bRtn;

    for (i = 0; i < nSrcRowCount; i++)
    {
        for (j = 0; j < nSrcColumnCount; j++)
        {
            CKNComplex          complexNumber = matrix.GetElement(nSrcRowStart + i, nSrcColumnStart + j);
            m_vectValueBuffer[m_nColumnCount*(nRowStart + i) + (nColumnStart + j)] = complexNumber;
        }
    }

    bRtn = true;
    return bRtn;
}

void CKNMatrixOperation::CKNDMatrix::AppendMatrix(APPEND_DRIECTION direction, unsigned int nCount)
{
    unsigned int                    nFormerRow = m_nRowCount;
    unsigned int                    nFormerColumn = m_nColumnCount;
    unsigned int                    i, j;
    std::vector<CKNComplex>         vectTemp;

    switch (direction)
    {
    case ROW_DIRECTION:
        m_nRowCount += nCount;
        break;
    case COLUMN_DIRECTION:
        m_nColumnCount += nCount;
        break;
    }

    for (i = 0; i < m_nRowCount * m_nColumnCount; ++i)
    {
        CKNComplex          element;
        vectTemp.push_back(element);
    }

    for (i = 0; i < nFormerRow; ++i)
    {
        for (j = 0; j < nFormerColumn; ++j)
        {
            CKNComplex          element = m_vectValueBuffer[nFormerColumn*i + j];
            vectTemp[m_nColumnCount*i + j].SetComplexNumber(element.GetRealNumber(), element.GetImaginaryNumber());
        }
    }
}

bool CKNMatrixOperation::CKNDMatrix::SetDiagonal(CKNVector vector)
{
    unsigned int                i;
    bool                        bRtn = false;

    if (vector.GetSize() != m_nRowCount || vector.GetSize() != m_nColumnCount)
        return bRtn;

    for (i = 0; i < m_nRowCount; ++i)
    {
        SetElement(i, i, vector.GetAt(i));
    }

    bRtn = true;
    return bRtn;
}

CKNComplex CKNMatrixOperation::CKNDMatrix::GetElement(unsigned int nRowIndex, unsigned int nColumnIndex)
{
    return m_vectValueBuffer[m_nColumnCount*nRowIndex + nColumnIndex];
}

bool CKNMatrixOperation::CKNDMatrix::SetColumnElement(CKNVector vector, unsigned int nColumnIndex)
{
    bool                bRtn = false;
    unsigned int        i;

    if (vector.GetSize() > m_nRowCount)
        return bRtn;

    for (i = 0; i < vector.GetSize(); ++i)
    {
        SetElement(i, nColumnIndex, vector.GetAt(i));
    }

    bRtn = true;
    return bRtn;
}

bool CKNMatrixOperation::CKNDMatrix::SetRowElement(CKNVector vector, unsigned int nRowIndex)
{
    bool                bRtn = false;
    unsigned int        i;

    if (vector.GetSize() > m_nColumnCount)
        return bRtn;

    for (i = 0; i < vector.GetSize(); ++i)
    {
        SetElement(nRowIndex, i, vector.GetAt(i));
    }

    bRtn = true;
    return bRtn;
}

void CKNMatrixOperation::CKNDMatrix::ScalarMultiple(CKNComplex Scalar)
{
    unsigned int            i;

    for (i = 0; i < m_vectValueBuffer.size(); i++)
    {
        m_vectValueBuffer[i] = m_vectValueBuffer[i] * Scalar;
    }
}

void CKNMatrixOperation::CKNDMatrix::ScalarMultiple(double fScalar)
{
    double                  fReal, fImaginary;
    unsigned int            i;

    for (i = 0; i < m_vectValueBuffer.size(); i++)
    {
        fReal = m_vectValueBuffer[i].GetRealNumber();
        fImaginary = m_vectValueBuffer[i].GetImaginaryNumber();
        m_vectValueBuffer[i].SetComplexNumber(fScalar*fReal, fScalar*fImaginary);
    }
}


bool CKNMatrixOperation::CKNDMatrix::GetColumnByVector(unsigned int nColumnIndex, CKNMatrixOperation::CKNVector *pVector)
{
    bool                        bRtn = false;
    unsigned int                i;

    if (nColumnIndex > m_nColumnCount)
        return bRtn;

    pVector->SetSize(m_nRowCount);
    for (i = 0; i < m_nRowCount; ++i)
        pVector->SetAt(i, m_vectValueBuffer[i*m_nColumnCount + nColumnIndex]);

    bRtn = true;
    return bRtn;
}

bool CKNMatrixOperation::CKNDMatrix::GetRowByVector(unsigned int nRowIndex, CKNMatrixOperation::CKNVector *pVector)
{
    bool                        bRtn = false;
    unsigned int                i;

    if (nRowIndex > m_nRowCount)
        return bRtn;

    pVector->SetSize(m_nColumnCount);
    for (i = 0; i < m_nColumnCount; ++i)
        pVector->SetAt(i, m_vectValueBuffer[nRowIndex*m_nColumnCount + i]);

    bRtn = true;
    return bRtn;
}

bool CKNMatrixOperation::CKNDMatrix::GetSmallMatrix(unsigned int nRowStartIndex, unsigned int nColumnStartIndex, unsigned int nRowCount, unsigned int nColumnCount, CKNMatrixOperation::CKNDMatrix *pMatrix)
{
    bool                    bRtn = false;
    unsigned int            i, j;

    if (nRowStartIndex + nRowCount > m_nRowCount)
        return bRtn;

    if (nColumnCount + nColumnCount > m_nColumnCount)
        return bRtn;

    if (pMatrix->GetColumnCount() < nColumnCount || pMatrix->GetRowCount() < nRowCount)
        return bRtn;


    for (i = nRowStartIndex; i < nRowStartIndex + nRowCount; i++)
    {
        for (j = nColumnStartIndex; j < nColumnStartIndex + nColumnCount; j++)
        {
            pMatrix->SetElement(i - nRowStartIndex, j - nColumnStartIndex, m_vectValueBuffer[i*m_nColumnCount + j]);
        }
    }

    bRtn = true;
    return bRtn;
}

void CKNMatrixOperation::CKNDMatrix::ScalarDivision(double fScalar)
{
    unsigned int            nCount = m_vectValueBuffer.size(), i;

    for (i = 0; i < nCount; ++i)
        m_vectValueBuffer[i].Division(2.0);
}

void CKNMatrixOperation::CKNDMatrix::operator+=(CKNDMatrix *matrix)
{
    operator+=((*matrix));
}

void CKNMatrixOperation::CKNDMatrix::operator+=(CKNDMatrix &matrix)
{
    if (matrix.GetColumnCount() != GetColumnCount() || matrix.GetRowCount() != GetRowCount())
        return;

    unsigned int        i, j;

    for (i = 0; i < m_nRowCount; ++i)
    {
        for (j = 0; j < m_nColumnCount; ++j)
        {
            m_vectValueBuffer[m_nColumnCount*i + j] = m_vectValueBuffer[m_nColumnCount*i + j] + matrix.GetElement(i, j);
        }
    }
}

void CKNMatrixOperation::CKNDMatrix::operator=(CKNDMatrix *matrix)
{
    operator=((*matrix));
}

void CKNMatrixOperation::CKNDMatrix::operator=(CKNDMatrix &matrix)
{
    unsigned int        i, j;

    BuildMatrixFirst(matrix.GetRowCount(), matrix.GetColumnCount());

    for (i = 0; i < m_nRowCount; i++)
    {
        for (j = 0; j < m_nColumnCount; j++)
        {
            SetElement(i, j, matrix.GetElement(i, j));
        }
    }
}

bool CKNMatrixOperation::CKNDMatrix::TrnasPos()
{
    bool                        bRtn = false;
    unsigned int                i, j, nTemp;
    CKNComplex                  tempNumber;


    if (m_nColumnCount == m_nRowCount)
    {
        for (i = 0; i < m_nRowCount; ++i)
        {
            for (j = 0; j < m_nColumnCount; ++j)
            {
                if (j <= i)
                    continue;

                tempNumber = m_vectValueBuffer[i*m_nColumnCount + j];
                m_vectValueBuffer[i*m_nColumnCount + j] = m_vectValueBuffer[j*m_nColumnCount + i];
                m_vectValueBuffer[j*m_nColumnCount + i] = tempNumber;
            }
        }
    }
    else
    {
        std::vector<CKNComplex>     tempVector;

        for (i = 0; i < m_nColumnCount; i++)
        {
            for (j = 0; j < m_nRowCount; j++)
            {
                CKNComplex          element;
                tempVector.push_back(element);
            }
        }


        for (i = 0; i < m_nRowCount; ++i)
            for (j = 0; j < m_nColumnCount; ++j)
                tempVector[j*m_nRowCount + i] = m_vectValueBuffer[i*m_nColumnCount + j];

        nTemp = m_vectValueBuffer.size();
        for (i = 0; i < nTemp; ++i)
            m_vectValueBuffer[i] = tempVector[i];

        nTemp = m_nRowCount;
        m_nRowCount = m_nColumnCount;
        m_nColumnCount = nTemp;
    }

    bRtn = true;
    return bRtn;
}


CKNMatrixOperation::CKNCSR::CKNCSR()
{
    m_nValueCount = 0; 
    m_nValueStackCount = 0; 
    m_nRowCount = 0; 
    m_nColumnCount = 0; 
    MAX_INDEX = -1; 
    m_fFirstRowIndex = 0;
    nComponentsFirstUnitCell = 0;
    nComponentsLastUnitCell = 0;
}

CKNMatrixOperation::CKNCSR::~CKNCSR()
{

}

void CKNMatrixOperation::CKNCSR::IncreaseNoneZeroCount()
{
    m_nValueCount++;
}

void CKNMatrixOperation::CKNCSR::BuildDataBuffer()
{
    for (unsigned int i = 0; i < m_nRowCount + 1; i++)
        m_vectRow.push_back(MAX_INDEX);
}

int compare(const void *pA, const void *pB)
{
    CKNMatrixOperation::LPFILL_MATRIX_DATA      lpA = (CKNMatrixOperation::LPFILL_MATRIX_DATA)pA;
    CKNMatrixOperation::LPFILL_MATRIX_DATA      lpB = (CKNMatrixOperation::LPFILL_MATRIX_DATA)pB;

    if (NULL == lpA->pMatrix)
        return 1;
    if (NULL == lpB->pMatrix)
        return -1;

    if (lpA->nColumnIndex > lpB->nColumnIndex)
        return 1;
    else
        return -1;

    return -1;
}

bool CKNMatrixOperation::CKNCSR::PushMatrixConcurrent(unsigned int nRow, CKNMatrixOperation::LPFILL_MATRIX_DATA lpData, bool bCopyZeroOnSite)
{
    bool                                bRtn = false;
    int                                 nValueCount = 0, i, j, k, nValidIndex[5];
    unsigned int                        nHonSiteIndex = lpData[0].nColumnIndex;
    CKNComplex                          tempNumber;

    for (i = 0; i < 4; ++i)
    {
        if (NULL != lpData[i + 1].pMatrix)
            nValidIndex[nValueCount++] = i + 1;
    }
    nValidIndex[nValueCount++] = 0; 
    qsort(lpData, 5, sizeof(CKNMatrixOperation::FILL_MATRIX_DATA), compare);

    for (i = 0; i < ORBITALS; ++i) 
    {
        for (j = 0; j < nValueCount; ++j) 
        {
            if( ATOM_DEFAULT_INDEX == lpData[j].nColumnIndex )
                continue;

            for (k = 0; k < ORBITALS; ++k) 
            {
                tempNumber = lpData[j].pMatrix->GetElement(i, k);
                if (nHonSiteIndex == lpData[j].nColumnIndex && false == bCopyZeroOnSite && 0 == tempNumber.GetRealNumber() && 0 == tempNumber.GetImaginaryNumber())
                    continue;

                PushNoneZeroValue(tempNumber.GetRealNumber(), tempNumber.GetImaginaryNumber(), nRow * ORBITALS + i, lpData[j].nColumnIndex *ORBITALS + k);
            }
        }
    }
    FinishPush();

    bRtn = true;
    return bRtn;
}

bool CKNMatrixOperation::CKNCSR::PushMatrix(unsigned int nRow, unsigned int nColumn, unsigned int nRowStart, unsigned int nColumnStart, unsigned int nRowCount, unsigned int nColumnCount, CKNMatrixOperation::CKNDMatrix *pMatrix, bool bCopyZero)
{
    bool                                bRtn = false;
    CKNComplex                          tempNumber;
    unsigned int                        i, j;

    if (m_nRowCount < nRow + nRowCount - nRowStart || m_nColumnCount < nColumn + nColumnCount - nColumnStart)
        return bRtn;

    for (i = nRowStart; i < nRowCount; i++)
    {
        for (j = nColumnStart; j < nColumnCount; j++)
        {
            tempNumber = pMatrix->GetElement(i, j);
            if (false == bCopyZero && 0 == tempNumber.GetRealNumber() && 0 == tempNumber.GetImaginaryNumber())
                continue;

            PushNoneZeroValue(tempNumber.GetRealNumber(), tempNumber.GetImaginaryNumber(), nRow + i - nRowStart, nColumn + j - nColumnStart);
        }
    }
    FinishPush();

    bRtn = true;
    return bRtn;
}

bool CKNMatrixOperation::CKNCSR::AreaScalarMultiple(unsigned int nRowStart, unsigned int nRowCount, unsigned int nColumnStart, unsigned int nColumnCount, CKNComplex Scalar)
{
    bool                    bRtn = false;
    unsigned int            i, j;

    for (i = nRowStart; i < nRowStart + nRowCount; ++i)
    {
        for (j = nColumnStart; j < nColumnStart + nColumnCount; ++j)
        {
            bRtn = ElementScalarMultiple(i, j, Scalar);
            if (!bRtn)
            {
                return bRtn;
            }
        }
    }

    return bRtn;
}

bool CKNMatrixOperation::CKNCSR::ElementScalarMultiple(unsigned int nRow, unsigned int nColumn, CKNComplex Scalar)
{
    CKNComplex              tempNumber;
    bool                    bResult, bRtn = false;;

    tempNumber = GetElement(nRow, nColumn, bResult);
    if (!bResult)
        return bRtn;

    tempNumber = tempNumber * Scalar;
    bRtn = SetAt(tempNumber, nRow, nColumn);

    return bRtn;
}

bool CKNMatrixOperation::CKNCSR::ElementScalarMultiple(unsigned int nRow, unsigned int nColumn, double fScalar)
{
    CKNComplex              tempNumber;
    bool                    bResult, bRtn = false;;

    tempNumber = GetElement(nRow, nColumn, bResult);
    if (!bResult)
        return bRtn;

    tempNumber = tempNumber * fScalar;
    bRtn = SetAt(tempNumber, nRow, nColumn);

    return bRtn;
}

bool CKNMatrixOperation::CKNCSR::InsertMatrix(unsigned int nRow, unsigned int nColumn, unsigned int nRowStart, unsigned int nColumnStart, unsigned int nRowCount, unsigned int nColumnCount, CKNMatrixOperation::CKNDMatrix *pMatrix, bool bCopyZero)
{
    unsigned int                        nMatrixRowCount = pMatrix->GetRowCount();
    unsigned int                        nMatrixColumnCount = pMatrix->GetColumnCount();
    bool                                bRtn = false;
    CKNComplex                          tempNumber;
    unsigned int                        i, j;

    if (m_nRowCount < nRow + nRowCount - nRowStart || m_nColumnCount < nColumn + nColumnCount - nColumnStart)
        return bRtn;

    for (i = nRowStart; i < nRowCount; i++)
    {
        for (j = nColumnStart; j < nColumnCount; j++)
        {
            tempNumber = pMatrix->GetElement(i, j);
            if (false == bCopyZero && 0 == tempNumber.GetRealNumber() && 0 == tempNumber.GetImaginaryNumber())
                continue;

            SetAt(tempNumber, nRow + i - nRowStart, nColumn + j - nColumnStart);
        }
    }

    bRtn = true;
    return bRtn;
}

unsigned int CKNMatrixOperation::CKNCSR::GetRowIndexNo(unsigned int nIndex)
{
    if (nIndex > GetRowCount() + 1)
    {
        throw ERROR_OUT_OF_RANGE;
        return MAX_INDEX;
    }

    return m_vectRow[nIndex];
}

unsigned int CKNMatrixOperation::CKNCSR::GetColIndexNo(unsigned int nIndex)
{
    if (nIndex > GetNoneZeroCount())
    {
        throw ERROR_OUT_OF_RANGE;
        return MAX_INDEX;
    }

    //return m_pnColum[nIndex];
    return m_vectColumn[nIndex];
}

CKNComplex* CKNMatrixOperation::CKNCSR::GetValue(unsigned int nIndex)
{
    CKNComplex      dumyValue;
    if (nIndex > GetNoneZeroCount())
    {
        throw ERROR_OUT_OF_RANGE;
        return NULL;
    }

    m_rtnTemp.SetComplexNumber(m_vectValueRealBuffer[nIndex], m_vectValueImaginaryBuffer[nIndex]);
    return &m_rtnTemp;
}

void CKNMatrixOperation::CKNCSR::PushNoneZeroValue(double fRealValue, double fImaginaryValue, unsigned int nRow, unsigned int nCol)
{
    m_vectValueRealBuffer.push_back(fRealValue);
    m_vectValueImaginaryBuffer.push_back(fImaginaryValue);
    m_vectColumn.push_back(nCol);

    if (MAX_INDEX == GetRowIndexNo(nRow))
        m_vectRow[nRow] = m_nValueStackCount;

    m_nValueStackCount++;
    m_nValueCount++;
}

void CKNMatrixOperation::CKNCSR::FinishPush()
{
    m_vectRow[GetRowCount()] = GetNoneZeroCount();
}

bool CKNMatrixOperation::CKNCSR::ConvertDoubleArray(unsigned int *pRowPtr, unsigned int *pColIndex, double *pNNZValueReal, double *pNNZValueImaginary, unsigned int nNNZSize, unsigned int nRowSize, unsigned int nColSize, unsigned int nFirstIndex, bool bZerobase)
{
    bool                bRtn = false;
    int                 nAdjustIndex = 0;
    int                 i, j;
    int                 nStartIndex, nEndIndex;
    double              fReal, fImaginary;

    SetRowCount(nRowSize);
    SetColumnCount(nColSize);
    BuildDataBuffer();

    if( !bZerobase )
        nAdjustIndex = -1;

    for( i = 0; i < nRowSize ; ++ i )
    {
        nStartIndex = pRowPtr[i] + nAdjustIndex;
        nEndIndex = pRowPtr[i+1] + nAdjustIndex;
        for( j = nStartIndex ; j < nEndIndex ; ++ j)
        {
            if( pNNZValueReal )
                fReal = pNNZValueReal[j];
            else
                fReal = 0.;
            
            if( pNNZValueImaginary )
                fImaginary = pNNZValueImaginary[j];
            else
                fImaginary = 0.;

            PushNoneZeroValue(fReal, fImaginary, i, pColIndex[j] + nAdjustIndex);
        }
    }
    
    FinishPush();
    SetFirstRowIndex(nFirstIndex + nAdjustIndex);

    bRtn = true;
    return bRtn;
}

void CKNMatrixOperation::CKNCSR::Finalize()
{
    if (NULL == this)
        return;

    if (m_vectColumn.size() > 0)
        m_vectColumn.clear();

    if (m_vectValueRealBuffer.size() > 0)
        m_vectValueRealBuffer.clear();

    if (m_vectValueImaginaryBuffer.size() > 0)
        m_vectValueImaginaryBuffer.clear();

    if (m_vectRow.size() > 0)
        m_vectRow.clear();
}

void CKNMatrixOperation::CKNCSR::ExpandMatrix(unsigned int nMulti, bool bRow, bool bColumn)
{
    unsigned int        nOriginSize, i;

    if (false == bRow && false == bColumn)
        return;

    if (bRow)
    {
        nOriginSize = GetRowCount();
        SetRowCount(nMulti*nOriginSize);

        uint_vector_t temp = m_vectRow;

        m_vectRow.clear();
        for (i = 0; i < GetRowCount(); i++)
            m_vectRow.push_back(MAX_INDEX);
        m_vectRow.push_back(GetNoneZeroCount());

        for (i = 0; i < temp.size() - 1; i++)
            m_vectRow[i*nMulti] = temp[i];
    }

    if (bColumn)
    {
        nOriginSize = GetColumnCount();
        SetColumnCount(nMulti*nOriginSize);
        for (i = 0; i < GetNoneZeroCount(); i++)
            m_vectColumn[i] *= nMulti;
    }
}

bool CKNMatrixOperation::CKNCSR::GetNextRowIndexValue(unsigned int nRowFrom, unsigned int &nValueIndex)
{
    bool            bRtn = false;
    unsigned int    nRowIndex = nRowFrom;

    if (nRowFrom > m_nRowCount)
        return bRtn;

    while (MAX_INDEX == m_vectRow[nRowIndex])
    {
        nRowIndex++;
        if (nRowIndex > m_nRowCount)
            return bRtn;
    }

    nValueIndex = m_vectRow[nRowIndex];

    bRtn = true;
    return bRtn;
}

void CKNMatrixOperation::CKNCSR::IncreaseRowIndex(unsigned int nRowFrom)
{
    unsigned int                i;

    for (i = nRowFrom; i < GetRowCount(); i++)
    {
        if (MAX_INDEX != m_vectRow[i])
            m_vectRow[i]++;
    }
}

bool CKNMatrixOperation::CKNCSR::SetAt(CKNComplex number, unsigned int nRow, unsigned int nColumn)
{
    bool                                    bRtn = false;
    unsigned int                            nIndex;
    unsigned int                            nRowStart;
    CKNComplex                              bufferElement = number, tempElement;
    unsigned int                            bufferColIndex = nColumn;
    double_vector_t::iterator   valueIndex = m_vectValueRealBuffer.begin();
    uint_vector_t::iterator     columnIndex = m_vectColumn.begin();


    if (nRow >= GetRowCount() || nColumn >= GetColumnCount())
        return bRtn;

    if (IsNonzeroElement(nRow, nColumn, nIndex))
    {
        m_vectValueRealBuffer[nIndex] = number.GetRealNumber();
        m_vectValueImaginaryBuffer[nIndex] = number.GetImaginaryNumber();
    }
    else
    {
        nRowStart = GetRowIndexNo(nRow);
        if (MAX_INDEX != nRowStart)
        {
            unsigned int            nInsertPos;

            if (GetColIndexNo(nRowStart) > nColumn)
            {
                nInsertPos = nRowStart;
            }
            else
            {
                nInsertPos = GetNextNonzeroValueIndex(nRow, nColumn);
            }

            if (MAX_INDEX == nInsertPos)
            {
                PushNoneZeroValue(number.GetRealNumber(), number.GetImaginaryNumber(), nRow, nColumn);
                FinishPush();
            }
            else
            {
                m_vectValueRealBuffer.insert(valueIndex + nInsertPos, bufferElement.GetRealNumber());
                m_vectValueImaginaryBuffer.insert(valueIndex + nInsertPos, bufferElement.GetImaginaryNumber());
                m_vectColumn.insert(columnIndex + nInsertPos, nColumn);
                IncreaseRowIndex(nRowStart + 1);
                IncreaseNoneZeroCount();
                FinishPush();
            }
        }
        else
        {
            if (GetNextRowIndexValue(nRow + 1, nRowStart))
            {
                m_vectValueRealBuffer.insert(valueIndex + nRowStart, bufferElement.GetRealNumber());
                m_vectValueImaginaryBuffer.insert(valueIndex + nRowStart, bufferElement.GetImaginaryNumber());
                m_vectColumn.insert(columnIndex + nRowStart, nColumn);
                m_vectRow[nRow] = nRowStart;
                IncreaseRowIndex(nRowStart + 1);
                IncreaseNoneZeroCount();
                FinishPush();
            }
            else
            {
                PushNoneZeroValue(number.GetRealNumber(), number.GetImaginaryNumber(), nRow, nColumn);
                FinishPush();
            }
        }
    }

    bRtn = true;
    return bRtn;
}

CKNComplex CKNMatrixOperation::CKNCSR::GetElement(unsigned int nRow, unsigned int nColumn, bool &bResult)
{
    CKNComplex              elementValue;
    unsigned int            nSubStart = GetRowIndexNo(nRow), nSubEnd = GetRowIndexNo(nRow + 1), i;

    bResult = false;
    if (MAX_INDEX == nSubStart)
        return elementValue;

    if (MAX_INDEX == nSubEnd)
    {
        unsigned int nIndex = nRow + 2;

        while (MAX_INDEX == nSubEnd && nIndex <= GetRowCount())
            nSubEnd = GetRowIndexNo(nIndex++);
    }

    for (i = nSubStart; i < nSubEnd; i++)
    {
        if (nColumn == m_vectColumn[i])
        {
            bResult = true;
            elementValue.SetComplexNumber(m_vectValueRealBuffer[i], m_vectValueImaginaryBuffer[i]);
            return elementValue;
        }
    }

    return elementValue;
}


bool CKNMatrixOperation::CKNCSR::IsNonzeroElement(unsigned int nRow, unsigned int nColumn, unsigned int &nIndex)
{
    bool            bRtn = false;
    unsigned int    i;

    unsigned int    nSubStart = GetRowIndexNo(nRow);
    unsigned int    nSubEnd = GetRowIndexNo(nRow + 1);

    if (0 == GetNoneZeroCount())
        return bRtn;

    if (MAX_INDEX == nSubStart)
        return bRtn;

    if (MAX_INDEX == nSubEnd)
    {
        if (!GetNextRowIndexValue(nRow + 2, nSubEnd))
        {
            return bRtn;
        }
    }

    for (i = nSubStart; i < nSubEnd; i++)
    {
        unsigned int nColIndex = GetColIndexNo(i);
        if (nColIndex == nColumn)
        {
            nIndex = i;
            bRtn = true;
            return bRtn;
        }
    }

    return bRtn;
}

unsigned int CKNMatrixOperation::CKNCSR::GetNextNonzeroValueIndex(unsigned int nRow, unsigned int nColumn)
{
    unsigned int    i;

    unsigned int    nSubStart = GetRowIndexNo(nRow), nSubEnd = GetRowIndexNo(nRow + 1);
    unsigned int    nTemp = nRow;

    if (MAX_INDEX == nSubEnd)
    {
        if (!GetNextRowIndexValue(nRow + 2, nSubEnd))
        {
            return MAX_INDEX;
        }
    }

    for (i = nSubStart; i < nSubEnd; i++)
    {
        unsigned int nColIndex = GetColIndexNo(i);
        if (nColIndex > nColumn)
        {
            return i;
        }
    }
    return nSubEnd;
}

bool CKNMatrixOperation::CKNCSR::InsertRowBefore(unsigned int nRow)
{
    bool            bRtn = false;
    unsigned int    i;
    unsigned int    tempRowData, bufferRowData;

    if (nRow >= GetRowCount())
        return bRtn;

    bufferRowData = -1;
    for (i = nRow; i < GetRowCount(); i++)
    {
        tempRowData = GetRowIndexNo(i);
        m_vectRow[i] = bufferRowData;
        bufferRowData = tempRowData;
    }
    m_vectRow.push_back(bufferRowData);
    SetRowCount(GetRowCount() + 1);

    bRtn = true;
    return bRtn;
}

bool CKNMatrixOperation::CKNCSR::InsertRowAtEnd()
{
    m_vectRow[GetRowCount()] = -1;
    SetRowCount(GetRowCount() + 1);
    m_vectRow.push_back(GetNoneZeroCount());

    return true;
}

bool CKNMatrixOperation::CKNCSR::InsertColumnBefore(unsigned int nColumn)
{
    bool            bRtn = false;
    unsigned int    i, nNonZeroValueCount;

    if (nColumn >= GetColumnCount())
        return bRtn;

    nNonZeroValueCount = GetNoneZeroCount();

    for (i = 0; i < nNonZeroValueCount; i++)
    {
        if (m_vectColumn[i] >= nColumn)
            m_vectColumn[i]++;
    }

    SetColumnCount(GetColumnCount() + 1);

    bRtn = true;
    return bRtn;
}

bool CKNMatrixOperation::CKNCSR::InsertColumnAtEnd()
{
    SetColumnCount(GetColumnCount() + 1);
    return true;
}

CKNMatrixOperation::CKNVector CKNMatrixOperation::CKNCSR::operator*(CKNVector  & vector)
{
    CKNMatrixOperation::CKNVector   rtnVector;
    unsigned int                    i, j, nSize = GetRowCount();

    if (nSize != vector.GetSize())
    {
        throw ERROR_WRONG_ORDER_OPERATION;
        return rtnVector;
    }

    rtnVector.SetSize(GetRowCount());
    for (i = 0; i < nSize; i++)
    {
        CKNComplex          fSubTotal;
        unsigned int        nSubStart = GetRowIndexNo(i), nSubEnd = GetRowIndexNo(i + 1);

        if (MAX_INDEX == nSubStart)
            continue;

        if (MAX_INDEX == nSubEnd)
        {
            unsigned int nIndex = i + 2;

            while (MAX_INDEX == nSubEnd && nIndex <= GetRowCount())
                nSubEnd = GetRowIndexNo(nIndex++);
        }

        for (j = nSubStart; j < nSubEnd; j++)
        {
            unsigned int nColIndex = GetColIndexNo(j);
            fSubTotal += (*(GetValue(j))*vector.GetAt(nColIndex));
        }
        rtnVector.SetAt(i, fSubTotal);
    }

    return rtnVector;
}

CKNMatrixOperation::CKNVector CKNMatrixOperation::CKNCSR::operator*(CKNVector  * vector)
{
    return operator*((*vector));
}


CKNMatrixOperation::CKNMatrixOperation()
{
}

CKNMatrixOperation::~CKNMatrixOperation()
{
}

#define         REPEAT_COUNT        1000

CKNMatrixOperation::CKNCSR * CKNMatrixOperation::BuildCSRFromFileTemp(FILE *fDataFile, unsigned int nRowOrder, unsigned int nColumnOrder, int nDataCount)
{
    CKNMatrixOperation::CKNCSR              *pCSR = new CKNMatrixOperation::CKNCSR();
    size_t                                  readSize;
    unsigned int                            i;
    CKNMatrixOperation::CSR_ELEMENT_DATA    Data[REPEAT_COUNT];

    if (NULL == pCSR)
        return pCSR;

    pCSR->SetRowCount(nRowOrder);
    pCSR->SetColumnCount(nColumnOrder);
    pCSR->BuildDataBuffer();

    while (0 != (readSize = fread(Data, sizeof(CSR_ELEMENT_DATA), REPEAT_COUNT, fDataFile)))
    {
        for (i = 0; i < readSize; i++)
            pCSR->PushNoneZeroValue(Data[i].fReal, Data[i].fImaginary, (unsigned int)Data[i].nRow - 1, (unsigned int)Data[i].nColumn - 1);
    }

    pCSR->FinishPush();

    return pCSR;
}

CKNMatrixOperation::CKNCSR * CKNMatrixOperation::BuildCSRFromFile_(FILE *fDataFile, unsigned int nRowOrder, unsigned int nColumnOrder, int nDataCount)
{
    CKNMatrixOperation::CKNCSR              *pCSR = new CKNMatrixOperation::CKNCSR();
    size_t                                  readSize;
    unsigned int                            i;
    CKNMatrixOperation::CSR_ELEMENT_DATA    Data[REPEAT_COUNT];
    memset(Data, NULL, sizeof(CKNMatrixOperation::CSR_ELEMENT_DATA)*REPEAT_COUNT);

    if (NULL == pCSR)
        return pCSR;

    pCSR->SetRowCount(nRowOrder);
    pCSR->SetColumnCount(nColumnOrder);
    pCSR->BuildDataBuffer();

    while (0 != (readSize = fread(Data, sizeof(CSR_ELEMENT_DATA), REPEAT_COUNT, fDataFile)))
    {
        for (i = 0; i < readSize; i++)
        {
            if (0 == Data[i].nRow && 0 == Data[i].nColumn)
                break;

            pCSR->PushNoneZeroValue(Data[i].fReal, Data[i].fImaginary, (unsigned int)Data[i].nRow - 1, (unsigned int)Data[i].nColumn - 1);
        }
    }

    pCSR->FinishPush();

    return pCSR;
}

int CKNMatrixOperation::Compare(const void *pA, const void *pB)
{
    CKNMatrixOperation::LPCSR_ELEMENT_DATA              lpA = (CKNMatrixOperation::LPCSR_ELEMENT_DATA)pA;
    CKNMatrixOperation::LPCSR_ELEMENT_DATA              lpB = (CKNMatrixOperation::LPCSR_ELEMENT_DATA)pB;

    if (lpA->nRow > lpB->nRow)
        return 1;
    else if (lpA->nRow < lpB->nRow)
        return -1;

    if (lpA->nRow == lpB->nRow)
    {
        if (lpA->nColumn > lpB->nColumn)
            return 1;
        else if (lpA->nColumn < lpB->nColumn)
            return -1;
        else if (lpB->nColumn == lpB->nColumn)
            return 0;
    }

    return -1;
}

CKNMatrixOperation::CKNCSR * CKNMatrixOperation::BuildCSRFromFileUnsortdata(FILE *fDataFile, unsigned int nRowOrder, unsigned int nColumnOrder, int nDataCount)
{
    LPCSR_ELEMENT_DATA                      lpData = (LPCSR_ELEMENT_DATA)malloc(sizeof(CSR_ELEMENT_DATA)* nDataCount);
    CKNMatrixOperation::CKNCSR              *pCSR = new CKNMatrixOperation::CKNCSR();
    int                                     i;

    if (lpData == NULL)
        return pCSR;

    for (i = 0; i < nDataCount; i++)
    {
        fread(&lpData[i], sizeof(CSR_ELEMENT_DATA), 1, fDataFile);
    }

    pCSR->SetRowCount(nRowOrder);
    pCSR->SetColumnCount(nColumnOrder);
    pCSR->BuildDataBuffer();

    qsort(lpData, nDataCount, sizeof(CSR_ELEMENT_DATA), CKNMatrixOperation::Compare);

    for (i = 0; i < nDataCount; i++)
    {
        pCSR->PushNoneZeroValue(lpData[i].fReal, lpData[i].fImaginary, (unsigned int)lpData[i].nRow - 1, (unsigned int)lpData[i].nColumn - 1);
    }
    pCSR->FinishPush();

    return pCSR;

}
CKNMatrixOperation::CKNCSR * CKNMatrixOperation::BuildCSRFromOneDimArray(double *pReal, double *pImaginary, unsigned int nRowOrder, unsigned int nColumnOrder)
{
    unsigned int                row, col;
    CKNMatrixOperation::CKNCSR  *pCSR = new CKNMatrixOperation::CKNCSR();

    if (NULL == pCSR)
        return pCSR;

    pCSR->SetRowCount(nRowOrder);
    pCSR->SetColumnCount(nColumnOrder);
    pCSR->BuildDataBuffer();

    unsigned int nIndex = 0;
    bool bPushedInRow = false;
    for (row = 0; row < nRowOrder; row++)
    {
        for (col = 0; col < nColumnOrder; col++)
        {
            if (0 != pReal[row*nColumnOrder + col] || 0 != pImaginary[row*nColumnOrder + col])
            {
                pCSR->PushNoneZeroValue(pReal[nIndex], pImaginary[nIndex], row, col);
                bPushedInRow = true;
            }
            nIndex++;
        }
    }
    pCSR->FinishPush();

    return pCSR;
}

CKNMatrixOperation::CKNCSR* CKNMatrixOperation::CKNCSR::SplitCSR(int nStart, int nEnd)
{
    unsigned int                        i, j;
    unsigned int                        nSubStart, nSubEnd;
    CKNMatrixOperation::CKNCSR          *pCSR = new CKNMatrixOperation::CKNCSR();

    pCSR->SetColumnCount(GetColumnCount());
    pCSR->SetRowCount(nEnd - nStart);

    pCSR->BuildDataBuffer();

    int nRowIndex = 0;
    for (i = nStart; i < (unsigned int)nEnd; i++)
    {
        nSubStart = GetRowIndexNo(i);
        nSubEnd = GetRowIndexNo(i + 1);
        for (j = nSubStart; j < nSubEnd; j++)
        {
            //CKNComplex        element = m_vectValueBuffer[j];
            //pCSR->PushNoneZeroValue(element.GetRealNumber(), element.GetImaginaryNumber(), nRowIndex, m_vectColumn[j]);

            pCSR->PushNoneZeroValue(m_vectValueRealBuffer[i], m_vectValueImaginaryBuffer[i], nRowIndex, m_vectColumn[j]);
        }
        nRowIndex++;
    }
    pCSR->FinishPush();

    return pCSR;

}

bool CKNMatrixOperation::CKNCSR::DiagonalOperation(CKNVector *pOperand, OPERATION_TYPE type, bool bUseSplitVector)
{
    bool                    bRtn = false, bResult;
    unsigned int            i, j, nRowIndex;
    unsigned int            nVectorIndex = 0, nVectorEnd = pOperand->GetSize();
    
    if( bUseSplitVector )
    {
        if( pOperand->GetSize() != m_nRowCount / 10 )
            return bRtn;
    }
    else
    {
        if( pOperand->GetSize() != m_nColumnCount / 10 )
            return bRtn;
    
        nVectorIndex = (unsigned int)m_fFirstRowIndex/10;
        nVectorEnd = nVectorIndex + m_nRowCount / 10;
    }

    nRowIndex = 0;
    for ( i = nVectorIndex; i < nVectorEnd ; ++ i )
    {
        CKNComplex          number = pOperand->GetAt(i);
        for( j = 0 ; j < ORBITALS ; ++ j )
        {
            CKNComplex      csrElement = GetElement(nRowIndex, nRowIndex + (unsigned int)m_fFirstRowIndex, bResult);
            if( !bResult )
                return bRtn;
            
            switch(type)
            {
            case PLUS:
                csrElement = csrElement + number;
                break;
            case MINUS:
                csrElement = csrElement - number;
                break;
            case MULTIPLE:
                csrElement = csrElement * number;
                break;
            case DIVISION:
                csrElement = csrElement / number;
                break;
            }

            SetAt(csrElement, nRowIndex, nRowIndex + (unsigned int)m_fFirstRowIndex);
            nRowIndex++;

        }
    }

    bRtn = true;
    return bRtn;
}

void CKNMatrixOperation::FreeCSR(CKNMatrixOperation::CKNCSR *pCSR)
{
    if (NULL == pCSR)
        return;

    pCSR->Finalize();

    delete pCSR;
    pCSR = NULL;
}

void CKNMatrixOperation::CKNCSR::DumpCSR(const char *pstrFileName)
{
    FILE        *out;
    unsigned int            i, nCount;

    out = fopen(pstrFileName, "wt");
    if (NULL != out)
    {
        fputs("None Zero Value\r\n", out);
        fputs("------------------------------------------\r\n", out);
        nCount = GetNoneZeroCount();
        for (i = 0; i < nCount; i++)
            fprintf(out, "%f + %fi\r\n", m_vectValueRealBuffer[i], m_vectValueImaginaryBuffer[i]);

        fputs("Column\r\n", out);
        fputs("------------------------------------------\r\n", out);
        nCount = GetNoneZeroCount();
        for (i = 0; i < nCount; i++)
            fprintf(out, "%u\r\n", m_vectColumn[i]);

        fputs("Row\r\n", out);
        fputs("------------------------------------------\r\n", out);
        nCount = GetRowCount() + 1;
        for (i = 0; i < nCount; i++)
            fprintf(out, "%u\r\n", GetRowIndexNo(i));

        fclose(out);
    }
}

void CKNMatrixOperation::MVMul(CKNCSR *pAMatrix, CKNVector *pVector, CKNVector *pResult)
{
    unsigned int                    i, j, nSize = pAMatrix->GetColumnCount();
    CKNVector                       *pOperandVector = NULL, VOperand;
    double                          *pMatrixReal = NULL, *pMatrixImaginary = NULL;
    double                          *pVectorReal = NULL, *pVectorImaginary = NULL;
    double                          *pResultReal = NULL, *pResultImaginary = NULL;
    unsigned int                    *pMatrixRow = NULL, *pMatrixColumn = NULL;

    VOperand = *pVector;
    pOperandVector = &VOperand;

#ifndef DISABLE_MPI_ROUTINE
    
    VOperand.SetSize(pAMatrix->GetColumnCount());
    CKNMPIManager::MergeVector(pVector, &VOperand, pAMatrix->GetColumnCount());
    pOperandVector = &VOperand;
#else
    pOperandVector = pVector;
#endif

    if (nSize != pOperandVector->GetSize())
    {
        throw ERROR_WRONG_ORDER_OPERATION;
        return;
    }

    nSize = pAMatrix->GetRowCount();
    pMatrixReal = pAMatrix->m_vectValueRealBuffer.data();
    pMatrixImaginary = pAMatrix->m_vectValueImaginaryBuffer.data();
    pMatrixRow = pAMatrix->m_vectRow.data();
    pMatrixColumn = pAMatrix->m_vectColumn.data();
    pVectorReal = pOperandVector->m_vectValueRealBuffer.data();
    pVectorImaginary = pOperandVector->m_vectValueImaginaryBuffer.data();
    pResultReal = pResult->m_vectValueRealBuffer.data();
    pResultImaginary = pResult->m_vectValueImaginaryBuffer.data();

    unsigned int input_real_size       = pOperandVector->m_vectValueRealBuffer.size();
    unsigned int input_imaginary_size  = pOperandVector->m_vectValueImaginaryBuffer.size();
    unsigned int output_real_size      = pResult->m_vectValueRealBuffer.size();
    unsigned int output_imaginary_size = pResult->m_vectValueImaginaryBuffer.size();

    #pragma omp parallel for
    for ( i = 0; i < nSize; i++)
    {
        double real_sum      = 0.0;
        double imaginary_sum = 0.0;
        const unsigned int nSubStart = pMatrixRow[i];
        const unsigned int nSubEnd   = pMatrixRow[i + 1];
        

        for ( j = nSubStart; j < nSubEnd; j++)
        {
            const unsigned int nColIndex = pMatrixColumn[j];
            const double m_real      = pMatrixReal[j];
            const double m_imaginary = pMatrixImaginary[j];
            const double v_real      = pVectorReal[nColIndex];
            const double v_imaginary = pVectorImaginary[nColIndex];

            real_sum      += m_real * v_real      - m_imaginary * v_imaginary;
            imaginary_sum += m_real * v_imaginary + m_imaginary * v_real;
        }

        pResultReal[i] = real_sum;
        pResultImaginary[i] = imaginary_sum;
    }
#ifndef DISABLE_MPI_ROUTINE
    VOperand.Finalize();
#endif //DISABLE_MPI_ROUTINE
}

void CKNMatrixOperation::MVMulEx_AsyncCommWithLocalBlocks(CKNMatrixOperation::CKNCSR *mylocalblock, CKNMatrixOperation::CKNCSR *leftlocalblock, CKNMatrixOperation::CKNCSR *rightlocalblock, CKNVector *pVector, CKNVector *pResult, double *X, double *Xrt, double *Xlt)
{
    unsigned int                            nSize, B, Brt, Blt;
    double                                  *pMatrixValueReal = NULL, *pOperandVectorReal = NULL, *pResultVectorReal = NULL;
    double                                  *pMatrixValueImaginary = NULL, *pOperandVectorImaginary = NULL, *pResultVectorImaginary = NULL;
    int                                     tag = 1002;
    int                                     myrank = CKNMPIManager::GetCurrentRank();
    int                                     ncpus = CKNMPIManager::GetTotalNodeCount();
    int                                     left_neighbor = (myrank - 1 + ncpus) % ncpus;
    int                                     right_neighbor = (myrank + 1) % ncpus;
    MPI_Request                             req_sr[2];
    MPI_Status                              stat_sr[2];
    // XXX jinpil: pRow, pColumn should be local in offload directives
    unsigned int                            *pRow = mylocalblock->m_vectRow.data();
    unsigned int                            *pColumn = mylocalblock->m_vectColumn.data();

#ifndef _WIN32
    __assume_aligned(X,   64);
    __assume_aligned(X,   64);
    __assume_aligned(Xrt, 64);

    __assume_aligned(pMatrixValueReal,        64);
    __assume_aligned(pMatrixValueImaginary,   64);
    __assume_aligned(pOperandVectorReal,      64);
    __assume_aligned(pOperandVectorImaginary, 64);
    __assume_aligned(pResultVectorReal,       64);
    __assume_aligned(pResultVectorImaginary,  64);
    __assume_aligned(pRow,                    64);
    __assume_aligned(pColumn,                 64);
#endif //_WIN32

    pMatrixValueReal = mylocalblock->m_vectValueRealBuffer.data();
    pOperandVectorReal = pVector->m_vectValueRealBuffer.data();
    pResultVectorReal = pResult->m_vectValueRealBuffer.data();

    pMatrixValueImaginary = mylocalblock->m_vectValueImaginaryBuffer.data();
    pOperandVectorImaginary = pVector->m_vectValueImaginaryBuffer.data();
    pResultVectorImaginary = pResult->m_vectValueImaginaryBuffer.data();;

    nSize = mylocalblock->GetRowCount();

    if (nSize != pVector->GetSize())
    {
        throw ERROR_WRONG_ORDER_OPERATION;
        return;
    }

    B = nSize;
    Brt = 0;
    Blt = 0;

    for (int ii = 0; ii< nSize; ii++)
    {
        X[2 * ii] = pOperandVectorReal[ii];
        X[2 * ii + 1] = pOperandVectorImaginary[ii];
    }

    CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::MV_COMM);

    MPI_Irecv(&Brt, 1, MPI_INT, right_neighbor, tag, CKNMPIManager::GetMPIComm(), &req_sr[0]);
    MPI_Isend(&B, 1, MPI_INT, left_neighbor, tag, CKNMPIManager::GetMPIComm(), &req_sr[1]);
    MPI_Waitall(2, req_sr, stat_sr); // now Brt has B of right neighbor.

    MPI_Irecv(Xrt, 2 * Brt, MPI_DOUBLE, right_neighbor, tag, CKNMPIManager::GetMPIComm(), &req_sr[0]); // receive from bottom neighbor
    MPI_Isend(X, 2 * B, MPI_DOUBLE, left_neighbor, tag, CKNMPIManager::GetMPIComm(), &req_sr[1]);   // send to top neighbor

    CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::MV_COMM);

    unsigned int input_size = X_largest * 2;
    unsigned int output_real_size = pResult->m_vectValueRealBuffer.size();
    unsigned int output_imaginary_size = pResult->m_vectValueImaginaryBuffer.size();


#ifdef _ENABLE_PAPI
    long long papi_values[4];
    PAPI_start(papi_event_set);
#endif

#pragma offload target(mic:phi_tid) \
    nocopy(*pMatrixValueReal       : REUSE) \
    nocopy(*pMatrixValueImaginary  : REUSE) \
    nocopy(*pRow                   : REUSE) \
    nocopy(*pColumn                : REUSE) \
    in(X[0:input_size]             : REUSE) \
    nocopy(*pResultVectorReal      : REUSE) \
    nocopy(*pResultVectorImaginary : REUSE)
    //*/
#pragma omp parallel for
    for (unsigned int i = 0; i < nSize; i++)
    {
        double real_sum      = 0.0;
        double imaginary_sum = 0.0;
        const unsigned int nSubStart = pRow[i    ];
        const unsigned int nSubEnd = pRow[i + 1];
        for (unsigned int j = nSubStart; j < nSubEnd; j++)
        {
            const unsigned int nColIndex = pColumn[j];
            const double m_real = pMatrixValueReal[j];
            const double m_imaginary = pMatrixValueImaginary[j];
            const double v_real = X[2 * nColIndex];
            const double v_imaginary = X[2 * nColIndex + 1];

            real_sum += m_real * v_real - m_imaginary * v_imaginary;
            imaginary_sum += m_real * v_imaginary + m_imaginary * v_real;
        }

        pResultVectorReal[i] = real_sum;
        pResultVectorImaginary[i] = imaginary_sum;
    }

#ifdef _ENABLE_PAPI
    PAPI_stop(papi_event_set, papi_values);
    printf("[LOCAL] L2 access = %lld | L2 miss = %lld | L2 miss rate = %g %\n", papi_values[0], papi_values[1], (papi_values[1] * 100.0) / papi_values[0]);
    printf("[LOCAL] L3 access = %lld | L3 miss = %lld | L3 miss rate = %g %\n", papi_values[2], papi_values[3], (papi_values[3] * 100.0) / papi_values[2]);
#endif

    CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::MV_COMM);

    MPI_Waitall(2, req_sr, stat_sr); // now Xrt has X of right neighbor.

    MPI_Irecv(&Blt, 1, MPI_INT, left_neighbor, tag, CKNMPIManager::GetMPIComm(), &req_sr[0]);
    MPI_Isend(&B, 1, MPI_INT, right_neighbor, tag, CKNMPIManager::GetMPIComm(), &req_sr[1]);
    MPI_Waitall(2, req_sr, stat_sr); // now Blt has B of left neighbor.


    MPI_Irecv(Xlt, 2 * Blt, MPI_DOUBLE, left_neighbor, tag, CKNMPIManager::GetMPIComm(), &req_sr[0]);
    MPI_Isend(X, 2 * B, MPI_DOUBLE, right_neighbor, tag, CKNMPIManager::GetMPIComm(), &req_sr[1]);

    CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::MV_COMM);

    pMatrixValueReal = rightlocalblock->m_vectValueRealBuffer.data();
    pMatrixValueImaginary = rightlocalblock->m_vectValueImaginaryBuffer.data();
    pRow = rightlocalblock->m_vectRow.data();
    pColumn = rightlocalblock->m_vectColumn.data();

#ifdef _ENABLE_PAPI
    PAPI_start(papi_event_set);
#endif

#pragma offload target(mic:phi_tid) \
    nocopy(*pMatrixValueReal       : REUSE) \
    nocopy(*pMatrixValueImaginary  : REUSE) \
    nocopy(*pRow                   : REUSE) \
    nocopy(*pColumn                : REUSE) \
    in(Xrt[0:input_size]           : REUSE) \
    nocopy(*pResultVectorReal      : REUSE) \
    nocopy(*pResultVectorImaginary : REUSE)
    //*/
#pragma omp parallel for
    for (unsigned int i = 0; i < nSize; i++)
    {
        double real_sum = 0.0;
        double imaginary_sum = 0.0;
        const unsigned int nSubStart = pRow[i    ];
        const unsigned int nSubEnd = pRow[i + 1];
        for (unsigned int j = nSubStart; j < nSubEnd; j++)
        {
            const unsigned int nColIndex = pColumn[j];
            const double m_real = pMatrixValueReal[j];
            const double m_imaginary = pMatrixValueImaginary[j];
            const double v_real = Xrt[2 * nColIndex];
            const double v_imaginary = Xrt[2 * nColIndex + 1];

            real_sum += m_real * v_real - m_imaginary * v_imaginary;
            imaginary_sum += m_real * v_imaginary + m_imaginary * v_real;
        }

        pResultVectorReal[i] += real_sum;
        pResultVectorImaginary[i] += imaginary_sum;
    }

#ifdef _ENABLE_PAPI
    PAPI_stop(papi_event_set, papi_values);
    printf("[RIGHT] L2 access = %lld | L2 miss = %lld | L2 miss rate = %g %\n", papi_values[0], papi_values[1], (papi_values[1] * 100.0) / papi_values[0]);
    printf("[RIGHT] L3 access = %lld | L3 miss = %lld | L3 miss rate = %g %\n", papi_values[2], papi_values[3], (papi_values[3] * 100.0) / papi_values[2]);
#endif

    CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::MV_COMM);
    MPI_Waitall(2, req_sr, stat_sr); // now Xlt has X of left neighbor.
    CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::MV_COMM);

    pMatrixValueReal = leftlocalblock->m_vectValueRealBuffer.data();
    pMatrixValueImaginary = leftlocalblock->m_vectValueImaginaryBuffer.data();
    pRow = leftlocalblock->m_vectRow.data();
    pColumn = leftlocalblock->m_vectColumn.data();

#ifdef _ENABLE_PAPI
    PAPI_start(papi_event_set);
#endif

#pragma offload target(mic:phi_tid) \
    nocopy(*pMatrixValueReal      : REUSE) \
    nocopy(*pMatrixValueImaginary : REUSE) \
    nocopy(*pRow                  : REUSE) \
    nocopy(*pColumn               : REUSE) \
    in(Xlt[0:input_size]                                : REUSE) \
    out(pResultVectorReal[0:output_real_size]           : REUSE) \
    out(pResultVectorImaginary[0:output_imaginary_size] : REUSE)
    //*/
#pragma omp parallel for
    for (unsigned int i = 0; i < nSize; i++)
    {
        double real_sum = 0.0;
        double imaginary_sum = 0.0;
        const unsigned int nSubStart = pRow[i    ];
        const unsigned int nSubEnd = pRow[i + 1];
        for (unsigned int j = nSubStart; j < nSubEnd; j++)
        {
            const unsigned int nColIndex = pColumn[j];
            const double m_real = pMatrixValueReal[j];
            const double m_imaginary = pMatrixValueImaginary[j];
            const double v_real = Xlt[2 * nColIndex];
            const double v_imaginary = Xlt[2 * nColIndex + 1];

            real_sum += m_real * v_real - m_imaginary * v_imaginary;
            imaginary_sum += m_real * v_imaginary + m_imaginary * v_real;
        }

        pResultVectorReal[i] += real_sum;
        pResultVectorImaginary[i] += imaginary_sum;
    }

#ifdef _ENABLE_PAPI
    PAPI_stop(papi_event_set, papi_values);
    printf("[LEFT] L2 access = %lld | L2 miss = %lld | L2 miss rate = %g %\n", papi_values[0], papi_values[1], (papi_values[1] * 100.0) / papi_values[0]);
    printf("[LEFT] L3 access = %lld | L3 miss = %lld | L3 miss rate = %g %\n", papi_values[2], papi_values[3], (papi_values[3] * 100.0) / papi_values[2]);
#endif
}

void CKNMatrixOperation::MVMulEx_Optimal(CKNCSR *pAMatrix, CKNVector *pVector, CKNVector *pResult, unsigned int nSizeFromPrevRank, unsigned int nSizeFromNextRank, CKNVector *VTemp, int nSizePHI)
{
    unsigned int                    nSize = pAMatrix->GetColumnCount();
    CKNVector                       *pOperandVector = NULL;
    double                          *pMatrixReal = NULL, *pMatrixImaginary = NULL;
    double                          *pVectorReal = NULL, *pVectorImaginary = NULL;
    double                          *pResultReal = NULL, *pResultImaginary = NULL;
    unsigned int                    *pMatrixRow = NULL, *pMatrixColumn = NULL;
    char                            sigval;
    unsigned int                    input_size1, input_size2, input_size3;
    unsigned int                    input_offset1, input_offset2, input_offset3, offsettmp[3];
    unsigned int                    output_size, output_offset;

#ifndef DISABLE_MPI_ROUTINE
    pOperandVector = VTemp;
    pVectorReal = pOperandVector->m_vectValueRealBuffer.data();
    pVectorImaginary = pOperandVector->m_vectValueImaginaryBuffer.data();

    //  memset(pVectorReal, 0, sizeof(double)*pOperandVector->GetSize());
    //  memset(pVectorImaginary, 0, sizeof(double)*pOperandVector->GetSize());  

    CKNMPIManager::MergeVectorEx_Optimal(pVector, pOperandVector, pAMatrix->GetColumnCount(), pAMatrix->m_fFirstRowIndex, nSizeFromPrevRank, nSizeFromNextRank, pAMatrix->nComponentsFirstUnitCell, pAMatrix->nComponentsLastUnitCell, offsettmp);

#else
    pOperandVector = pVector;
#endif

    if (nSize != pOperandVector->GetSize())
    {
        throw ERROR_WRONG_ORDER_OPERATION;
        return;
    }

    pMatrixReal = pAMatrix->m_vectValueRealBuffer.data();
    pMatrixImaginary = pAMatrix->m_vectValueImaginaryBuffer.data();
    pMatrixRow = pAMatrix->m_vectRow.data();
    pMatrixColumn = pAMatrix->m_vectColumn.data();
    pResultReal = pResult->m_vectValueRealBuffer.data();
    pResultImaginary = pResult->m_vectValueImaginaryBuffer.data();
    nSize = pAMatrix->GetRowCount();

    if (CKNMPIManager::GetTotalNodeCount() <= 3)
    {
        input_size1 = pOperandVector->m_vectValueRealBuffer.size();
        input_offset1 = 0;

#pragma offload_transfer target(mic:phi_tid) in(pVectorReal[input_offset1:input_size1]      : REUSE)
#pragma offload_transfer target(mic:phi_tid) in(pVectorImaginary[input_offset1:input_size1] : REUSE)
    }
    else
    {
        input_size1 = nSizeFromPrevRank;
        input_size2 = nSize;
        input_size3 = nSizeFromNextRank;
        input_offset1 = offsettmp[0];
        input_offset2 = offsettmp[1];
        input_offset3 = offsettmp[2];

#pragma offload_transfer target(mic:phi_tid) in(pVectorReal[input_offset1:input_size1]      : REUSE)
#pragma offload_transfer target(mic:phi_tid) in(pVectorImaginary[input_offset1:input_size1] : REUSE)
#pragma offload_transfer target(mic:phi_tid) in(pVectorReal[input_offset2:input_size2]      : REUSE)
#pragma offload_transfer target(mic:phi_tid) in(pVectorImaginary[input_offset2:input_size2] : REUSE)
#pragma offload_transfer target(mic:phi_tid) in(pVectorReal[input_offset3:input_size3]      : REUSE)
#pragma offload_transfer target(mic:phi_tid) in(pVectorImaginary[input_offset3:input_size3] : REUSE)
    }

    output_size = nSizePHI;
    output_offset = 0;

    // FIXME jinpil:
    // nocopy(pMatrixReal : REUSE)
    // correct directive syntax, but segmentation fault without *
    // Xeon Phi device cannot find the correct pointer value
    // compiler bug?
    // FIXME allocate pVectorReal, pVectorImaginary outside the Lanczos loop

#pragma offload target(mic:phi_tid) \
    nocopy(*pMatrixReal                 : REUSE) \
    nocopy(*pMatrixImaginary            : REUSE) \
    nocopy(*pMatrixRow                  : REUSE) \
    nocopy(*pMatrixColumn               : REUSE) \
    nocopy(*pVectorReal             : REUSE) \
    nocopy(*pVectorImaginary            : REUSE) \
    out(pResultReal[output_offset:output_size]      : REUSE) \
    out(pResultImaginary[output_offset:output_size] : REUSE) \
    signal(&sigval)
    //*/

#pragma omp parallel for
    for (unsigned int i = 0; i < nSizePHI; i++)
    {
        double real_sum = 0.0;
        double imaginary_sum = 0.0;
        const unsigned int nSubStart = pMatrixRow[i];
        const unsigned int nSubEnd = pMatrixRow[i + 1];


        for (unsigned int j = nSubStart; j < nSubEnd; j++)
        {
            const unsigned int nColIndex = pMatrixColumn[j];
            const double m_real = pMatrixReal[j];
            const double m_imaginary = pMatrixImaginary[j];
            const double v_real = pVectorReal[nColIndex];
            const double v_imaginary = pVectorImaginary[nColIndex];

            real_sum += m_real * v_real - m_imaginary * v_imaginary;
            imaginary_sum += m_real * v_imaginary + m_imaginary * v_real;
        }

        pResultReal[i] = real_sum;
        pResultImaginary[i] = imaginary_sum;
    }

#pragma omp parallel for
    for (unsigned int i = nSizePHI; i < nSize; i++)
    {
        double real_sum = 0.0;
        double imaginary_sum = 0.0;
        const unsigned int nSubStart = pMatrixRow[i];
        const unsigned int nSubEnd = pMatrixRow[i + 1];


        for (unsigned int j = nSubStart; j < nSubEnd; j++)
        {
            const unsigned int nColIndex = pMatrixColumn[j];
            const double m_real = pMatrixReal[j];
            const double m_imaginary = pMatrixImaginary[j];
            const double v_real = pVectorReal[nColIndex];
            const double v_imaginary = pVectorImaginary[nColIndex];

            real_sum += m_real * v_real - m_imaginary * v_imaginary;
            imaginary_sum += m_real * v_imaginary + m_imaginary * v_real;
        }

        pResultReal[i] = real_sum;
        pResultImaginary[i] = imaginary_sum;
    }

#pragma offload_wait target(mic:phi_tid) wait(&sigval)

}

void CKNMatrixOperation::MVMulOptimal(CKNCSR *pAMatrix, CKNVector *pVector, CKNVector *pResult)
{
    unsigned int                    i, j, nSize = pAMatrix->GetColumnCount();
    CKNVector                       *pOperandVector = NULL;
    double                          *pMatrixReal = NULL, *pMatrixImaginary = NULL;
    double                          *pVectorReal = NULL, *pVectorImaginary = NULL;
    double                          *pResultReal = NULL, *pResultImaginary = NULL;
    unsigned int                    *pMatrixRow = NULL, *pMatrixColumn = NULL;
#ifndef DISABLE_MPI_ROUTINE
    CKNVector                       VOperand;

    VOperand.SetSize(pAMatrix->GetColumnCount());
    CKNMPIManager::MergeVectorOptimal(pVector, &VOperand, pAMatrix->GetColumnCount(), pAMatrix->m_fFirstRowIndex);
    pOperandVector = &VOperand;
#else
    pOperandVector = pVector;
#endif

    if (nSize != pOperandVector->GetSize())
    {
        throw ERROR_WRONG_ORDER_OPERATION;
        return;
    }

    nSize = pAMatrix->GetRowCount();
    pMatrixReal = pAMatrix->m_vectValueRealBuffer.data();
    pMatrixImaginary = pAMatrix->m_vectValueImaginaryBuffer.data();
    pMatrixRow = pAMatrix->m_vectRow.data();
    pMatrixColumn = pAMatrix->m_vectColumn.data();
    pVectorReal = pOperandVector->m_vectValueRealBuffer.data();
    pVectorImaginary = pOperandVector->m_vectValueImaginaryBuffer.data();
    pResultReal = pResult->m_vectValueRealBuffer.data();
    pResultImaginary = pResult->m_vectValueImaginaryBuffer.data();

    unsigned int input_real_size       = pOperandVector->m_vectValueRealBuffer.size();
    unsigned int input_imaginary_size  = pOperandVector->m_vectValueImaginaryBuffer.size();
    unsigned int output_real_size      = pResult->m_vectValueRealBuffer.size();
    unsigned int output_imaginary_size = pResult->m_vectValueImaginaryBuffer.size();

// FIXME jinpil:
// nocopy(pMatrixReal : REUSE)
// correct directive syntax, but segmentation fault without *
// Xeon Phi device cannot find the correct pointer value
// compiler bug?
// FIXME allocate pVectorReal, pVectorImaginary outside the Lanczos loop
#pragma offload target(mic:phi_tid) \
    nocopy(*pMatrixReal      : REUSE) \
    nocopy(*pMatrixImaginary : REUSE) \
    nocopy(*pMatrixRow       : REUSE) \
    nocopy(*pMatrixColumn    : REUSE) \
    in(pVectorReal[0:input_real_size]             : LOCAL) \
    in(pVectorImaginary[0:input_imaginary_size]   : LOCAL) \
    out(pResultReal[0:output_real_size]           : REUSE) \
    out(pResultImaginary[0:output_imaginary_size] : REUSE)
//*/
#pragma omp parallel for
    for ( i = 0; i < nSize; i++)
    {
        double real_sum      = 0.0;
        double imaginary_sum = 0.0;
        const unsigned int nSubStart = pMatrixRow[i];
        const unsigned int nSubEnd   = pMatrixRow[i + 1];
        

        for ( j = nSubStart; j < nSubEnd; j++)
        {
            const unsigned int nColIndex = pMatrixColumn[j];
            const double m_real      = pMatrixReal[j];
            const double m_imaginary = pMatrixImaginary[j];
            const double v_real      = pVectorReal[nColIndex];
            const double v_imaginary = pVectorImaginary[nColIndex];

            real_sum      += m_real * v_real      - m_imaginary * v_imaginary;
            imaginary_sum += m_real * v_imaginary + m_imaginary * v_real;
        }

        pResultReal[i] = real_sum;
        pResultImaginary[i] = imaginary_sum;
    }
#ifndef DISABLE_MPI_ROUTINE
    VOperand.Finalize();
#endif //DISABLE_MPI_ROUTINE
}

bool CKNMatrixOperation::VVDot(CKNVector *pVector1, CKNVector *pVector2, CKNComplex *pResult)
{
    double              *pOp1Real = NULL, *pOp1Imaginary = NULL;
    double              *pOp2Real = NULL, *pOp2Imaginary = NULL;
    unsigned int        i, nSize = pVector1->GetSize();

    if (nSize != pVector2->GetSize())
    {
        throw ERROR_WRONG_ORDER_OPERATION;
        return false;
    }

    pOp1Real = pVector1->m_vectValueRealBuffer.data();
    pOp1Imaginary = pVector1->m_vectValueImaginaryBuffer.data();

    pOp2Real = pVector2->m_vectValueRealBuffer.data();
    pOp2Imaginary = pVector2->m_vectValueImaginaryBuffer.data();

    double          fReal = 0., fImaginary = 0.;
#pragma omp parallel for reduction(+:fReal, fImaginary)
    for (i = 0; i < nSize; i++)
    {
        /*fReal += pOp1Real[i] * pOp2Real[i] - pOp1Imaginary[i] * (-1 * pOp2Imaginary[i]);
        fImaginary += pOp1Real[i] * (-1 * pOp2Imaginary[i]) + pOp1Imaginary[i] * pOp2Real[i];*/
        fReal += pOp1Real[i] * pOp2Real[i] - (-1*pOp1Imaginary[i]) * pOp2Imaginary[i];
        fImaginary += pOp1Real[i] * pOp2Imaginary[i] + (-1*pOp1Imaginary[i]) * pOp2Real[i];
    }

    pResult->SetComplexNumber(fReal, fImaginary);
    CKNMPIManager::AllReduceComlex(pResult, CKNTimeMeasurement::VV_COMM);

    return true;

}

void CKNMatrixOperation::MVMul(CKNDMatrix *pMatrix, CKNVector *pVector, CKNVector *pResult)
{
    int                 i, j, nRow, nColumn;

    if (pMatrix->GetColumnCount() != pVector->GetSize())
        return;

    pResult->SetSize(pVector->GetSize());

    nRow = pMatrix->GetRowCount();
    nColumn = pMatrix->GetColumnCount();
    for (i = 0; i < nRow; ++i)
    {
        CKNComplex          result;
        for (j = 0; j < nColumn; ++j)
        {
            result.Add(CKNComplex::MulltiplyComplex(pMatrix->GetElement(i, j), pVector->GetAt(j)));
        }
        pResult->SetAt(i, result);
    }
}

void CKNMatrixOperation::MMMul(CKNDMatrix *pMatrix, CKNDMatrix *pMatrixOperand, CKNDMatrix *pResult)
{
    int                 i, j, k;
    int                 nRow, nColumn, nL;
    if (pMatrix->GetColumnCount() != pMatrixOperand->GetRowCount())
        return;

    pResult->BuildMatrixFirst(pMatrix->GetRowCount(), pMatrixOperand->GetColumnCount());

    nL = pMatrixOperand->GetColumnCount();
    nRow = pMatrix->GetRowCount();
    nColumn = pMatrix->GetColumnCount();
    for (k = 0; k < nL; ++k)
    {
        for (i = 0; i < nRow; ++i)
        {
            CKNComplex          result;
            for (j = 0; j < nColumn; ++j)
            {
                result.Add(CKNComplex::MulltiplyComplex(pMatrix->GetElement(i, j), pMatrixOperand->GetElement(j, k)));
            }
            pResult->SetElement(i, k, result);
        }
    }
}

bool CKNMatrixOperation::IsSame(double operand1, double operand2, double tol)
{
    if (fabs(operand1 - operand2) > tol)
        return false;
    else
        return true;
}

bool CKNMatrixOperation::IsSameA(double operand1, double operand2, double tol)
{
    if (fabs(fabs(operand1) - fabs(operand2)) > tol)
        return false;
    else
        return true;
}

int CKNMatrixOperation::Gram_schmidt(CKNVector *pVect1, CKNVector *pVect2)
{
    CKNComplex              result;
    CKNVector               tempVector;

    CKNMatrixOperation::VVDot(pVect1, pVect2, &result);
    if( IsSameA(result.GetRealNumber(), 0, GENERAL_TOLERANCE ))
        return 1;
        
    tempVector = *pVect1;
    tempVector.Normalize(true);
    pVect2->Normalize(true);
    
    CKNMatrixOperation::VVDot(&tempVector, pVect2, &result);
    tempVector.ScalarMultiple(result);

    pVect2->MinusVector(&tempVector);
    pVect2->Normalize(true);

    CKNMatrixOperation::VVDot(pVect1, pVect2, &result);
    if( IsSameA(result.GetRealNumber(), 0, GENERAL_TOLERANCE ))
        return 1;
    else
        return 0;


}

CKNMatrixOperation::CKNVector operator*(const CKNComplex fScalar, CKNMatrixOperation::CKNVector & vector)
{
    vector.ScalarMultiple(fScalar);
    return vector;
}

CKNMatrixOperation::CKNVector operator/(CKNMatrixOperation::CKNVector & vector, const CKNComplex fScalar)
{
    vector.ScalarDivision(fScalar);
    return vector;
}

void CKNMatrixOperation::AllocateLocalCSR(CKNMatrixOperation::CKNCSR **mine, CKNMatrixOperation::CKNCSR **left, CKNMatrixOperation::CKNCSR **right)
{
    unsigned int                    nRowCount, nColumnCount;
    int                             myrank = CKNMPIManager::GetCurrentRank();
    int                             ncpus = CKNMPIManager::GetTotalNodeCount();
    int                             left_neighbor = (myrank-1+ncpus)%ncpus; // top neighbor
    int                             right_neighbor = (myrank+1)%ncpus; // bottom neighbor
    unsigned int                    temp;

    nRowCount = CKNMPIManager::GetLoadBalanceCount(myrank);
    *mine = NULL; *left = NULL; *right = NULL;
    
    // Allocate mine
    *mine = new CKNMatrixOperation::CKNCSR();
    if (*mine == NULL)
        throw ERROR_MALLOC;
    nColumnCount = CKNMPIManager::GetLoadBalanceCount(myrank); 
    (*mine)->SetRowCount(nRowCount);
    (*mine)->SetColumnCount(nColumnCount);
    (*mine)->BuildDataBuffer(); temp = 0;
    for (int jj=0; jj<myrank; jj++)
        temp += CKNMPIManager::GetLoadBalanceCount(jj);
    (*mine)->SetFirstRowIndex((double)temp); // FirstRowIndex will be used in a bit different way: Starting "column" index.

//  printf("%d %d %d %d %d %d\n", myrank, left_neighbor, right_neighbor, (int)(*mine)->GetFirstRowIndex(), (*mine)->GetColumnCount(), (*mine)->GetNoneZeroCount());

    // Allocate left
    *left = new CKNMatrixOperation::CKNCSR();
    if (*left == NULL)
        throw ERROR_MALLOC;
    nColumnCount = CKNMPIManager::GetLoadBalanceCount(left_neighbor);
    (*left)->SetRowCount(nRowCount);
    (*left)->SetColumnCount(nColumnCount);
    (*left)->BuildDataBuffer(); temp = 0;
    for (int jj=0; jj<left_neighbor; jj++)
        temp += CKNMPIManager::GetLoadBalanceCount(jj);
    (*left)->SetFirstRowIndex((double)temp); // FirstRowIndex will be used in a bit different way: Starting "column" index.

    //Allocate right
    *right = new CKNMatrixOperation::CKNCSR();
    if (*right == NULL)
        throw ERROR_MALLOC;
    nColumnCount = CKNMPIManager::GetLoadBalanceCount(right_neighbor);
    (*right)->SetRowCount(nRowCount);
    (*right)->SetColumnCount(nColumnCount);
    (*right)->BuildDataBuffer(); temp = 0;
    for(int jj=0; jj<right_neighbor; jj++)
        temp += CKNMPIManager::GetLoadBalanceCount(jj);
    (*right)->SetFirstRowIndex((double)temp); // FirstRowIndex will be used in a bit different way: Starting "column" index.

    MPI_Barrier(CKNMPIManager::GetMPIComm());   
}

void CKNMatrixOperation::BuildLocalCSR(CKNMatrixOperation::CKNCSR *source, CKNMatrixOperation::CKNCSR *mine, CKNMatrixOperation::CKNCSR *left, CKNMatrixOperation::CKNCSR *right)
{
    CKNMatrixOperation::pRow                = source->m_vectRow.data();
    CKNMatrixOperation::pColumn             = source->m_vectColumn.data();
    //CKNMemoryManager::LPVECTOR_ELEMENTS       lpMatrixValueElement = NULL;
    CKNComplex                              *pData = NULL;
    unsigned int my_nnz, left_nnz, right_nnz;
    int isthisrowfilled;

    //lpMatrixValueElement = source->GetValueElement();

    // 1. Build left local block.
    
    left_nnz = 0;   

    for (unsigned int ii = 0; ii < source->GetRowCount(); ii++)
    {
        isthisrowfilled = -1;
        unsigned int nSubStart = CKNMatrixOperation::pRow[ii], nSubEnd = CKNMatrixOperation::pRow[ii+1];
        unsigned int startColIndex = (int)left->GetFirstRowIndex(), endColIndex = startColIndex + left->GetColumnCount() - 1;

        for (unsigned int jj = nSubStart; jj < nSubEnd; jj++)
        {
            unsigned int nColIndex = CKNMatrixOperation::pColumn[jj];

            if(startColIndex <= nColIndex && nColIndex <= endColIndex)
            {   
                isthisrowfilled = 0;
                left_nnz++;
                pData = source->GetValue(jj);
                //left->PushNoneZeroValue(lpMatrixValueElement->pfReal[jj], lpMatrixValueElement->pfImaginary[jj], ii, nColIndex-startColIndex);
                left->PushNoneZeroValue(pData->GetRealNumber(), pData->GetImaginaryNumber(), ii, nColIndex-startColIndex);
            }
        }

        if(isthisrowfilled == -1)
        {
            left_nnz++;
            left->PushNoneZeroValue(0.0, 0.0, ii, 0);
        }
    }

    left->FinishPush();

//  if(CKNMPIManager::IsRootRank())
//      printf("Left block conversion completed: left_nnz = %d (computed), %d (CSR-allocated)\n", left_nnz, left->GetNoneZeroCount());  
                        
    // 2. Build right block
    
    right_nnz = 0;

    for (unsigned int ii = 0; ii < source->GetRowCount(); ii++)
    {
        isthisrowfilled = -1;
        unsigned int nSubStart = CKNMatrixOperation::pRow[ii], nSubEnd = CKNMatrixOperation::pRow[ii+1];
        unsigned int startColIndex = (int)right->GetFirstRowIndex(), endColIndex = startColIndex + right->GetColumnCount() - 1;
        
        for (unsigned int jj = nSubStart; jj < nSubEnd; jj++)
        {
            unsigned int nColIndex = CKNMatrixOperation::pColumn[jj];
            if(startColIndex <= nColIndex && nColIndex <= endColIndex)
            {   
                isthisrowfilled = 0;
                right_nnz++;
                pData = source->GetValue(jj);
                //right->PushNoneZeroValue(lpMatrixValueElement->pfReal[jj], lpMatrixValueElement->pfImaginary[jj], ii, nColIndex-startColIndex);
                right->PushNoneZeroValue(pData->GetRealNumber(), pData->GetImaginaryNumber(), ii, nColIndex-startColIndex);
            }
        }

        if(isthisrowfilled == -1)
        {
            right_nnz++;
            right->PushNoneZeroValue(0.0, 0.0, ii, 0);
        }
    }
                    
    right->FinishPush();

//  if(CKNMPIManager::IsRootRank())
//      printf("Right block conversion completed: right_nnz = %d (computed), %d (CSR-allocated)\n", right_nnz, right->GetNoneZeroCount());

    // 3. Build my block
    
    my_nnz = 0;

    for (unsigned int ii = 0; ii < source->GetRowCount(); ii++)
    {
        unsigned int nSubStart = CKNMatrixOperation::pRow[ii], nSubEnd = CKNMatrixOperation::pRow[ii+1];
        unsigned int startColIndex = (int)mine->GetFirstRowIndex(), endColIndex = startColIndex + mine->GetColumnCount() - 1;

        for(unsigned int jj = nSubStart; jj < nSubEnd; jj++)
        {
            unsigned int nColIndex = CKNMatrixOperation::pColumn[jj];
            if(startColIndex <= nColIndex && nColIndex <= endColIndex)
            {
                my_nnz++;
                pData = source->GetValue(jj);
                //mine->PushNoneZeroValue(lpMatrixValueElement->pfReal[jj], lpMatrixValueElement->pfImaginary[jj], ii, nColIndex-startColIndex);
                mine->PushNoneZeroValue(pData->GetRealNumber(), pData->GetImaginaryNumber(), ii, nColIndex-startColIndex);
            }               
        }
    }

    mine->FinishPush();

//  if(CKNMPIManager::IsRootRank())
//      printf("My block conversion completed: my_nnz = %d (computed), %d (CSR-allocated)\n", my_nnz, mine->GetNoneZeroCount());

}

void CKNMatrixOperation::UpdateLocalCSR(CKNMatrixOperation::CKNCSR *source, CKNMatrixOperation::CKNCSR *mine, CKNMatrixOperation::CKNCSR *left, CKNMatrixOperation::CKNCSR *right)
{   
    CKNMatrixOperation::pRow                = source->m_vectRow.data();
    CKNMatrixOperation::pColumn             = source->m_vectColumn.data();
    //CKNMemoryManager::LPVECTOR_ELEMENTS       lpMatrixValueElement = NULL;
    CKNComplex                              *pData = NULL;
    unsigned int my_nnz, left_nnz, right_nnz;
    CKNComplex curval; 
    int isthisrowfilled;

    //lpMatrixValueElement = source->GetValueElement();

    // 1. Update left block
    
    left_nnz = 0;   

    for (unsigned int ii = 0; ii < source->GetRowCount(); ii++)
    {
        isthisrowfilled = -1;
        unsigned int nSubStart = CKNMatrixOperation::pRow[ii], nSubEnd = CKNMatrixOperation::pRow[ii+1];
        unsigned int startColIndex = (int)left->GetFirstRowIndex(), endColIndex = startColIndex + left->GetColumnCount() - 1;

        for (unsigned int jj = nSubStart; jj < nSubEnd; jj++)
        {
            unsigned int nColIndex = CKNMatrixOperation::pColumn[jj];

            if(startColIndex <= nColIndex && nColIndex <= endColIndex)
            {   
                isthisrowfilled = 0;
                left_nnz++;
                //curval.SetRealNumber(lpMatrixValueElement->pfReal[jj]); curval.SetImaginaryNumber(lpMatrixValueElement->pfImaginary[jj]);
                pData = source->GetValue(jj);
                curval.SetRealNumber(pData->GetRealNumber()); curval.SetImaginaryNumber(pData->GetImaginaryNumber());
                left->SetAt(curval, ii, nColIndex-startColIndex);
            }
        }

        if(isthisrowfilled == -1)
        {
            left_nnz++;
            curval.SetRealNumber(0.0); curval.SetImaginaryNumber(0.0);
            left->SetAt(curval, ii, 0);
        }
    }

//  if(CKNMPIManager::IsRootRank())
//      printf("Left block update completed: left_nnz = %d (computed), %d (CSR-allocated)\n", left_nnz, left->GetNoneZeroCount());  
                        
    // 2. Figure out nnz: right block
    
    right_nnz = 0;

    for (unsigned int ii = 0; ii < source->GetRowCount(); ii++)
    {
        isthisrowfilled = -1;
        unsigned int nSubStart = CKNMatrixOperation::pRow[ii], nSubEnd = CKNMatrixOperation::pRow[ii+1];
        unsigned int startColIndex = (int)right->GetFirstRowIndex(), endColIndex = startColIndex + right->GetColumnCount() - 1;
        
        for (unsigned int jj = nSubStart; jj < nSubEnd; jj++)
        {
            unsigned int nColIndex = CKNMatrixOperation::pColumn[jj];
            if(startColIndex <= nColIndex && nColIndex <= endColIndex)
            {   
                isthisrowfilled = 0;
                right_nnz++;
                //curval.SetRealNumber(lpMatrixValueElement->pfReal[jj]); curval.SetImaginaryNumber(lpMatrixValueElement->pfImaginary[jj]);
                pData = source->GetValue(jj);
                curval.SetRealNumber(pData->GetRealNumber()); curval.SetImaginaryNumber(pData->GetImaginaryNumber());
                right->SetAt(curval, ii, nColIndex-startColIndex);
            }

        }

        if(isthisrowfilled == -1)
        {
            right_nnz++;
            curval.SetRealNumber(0.0); curval.SetImaginaryNumber(0.0);
            right->SetAt(curval, ii, 0);
        }
    }

//  if(CKNMPIManager::IsRootRank())
//      printf("Right block update completed: right_nnz = %d (computed), %d (CSR-allocated)\n", right_nnz, right->GetNoneZeroCount());

    // 3. Figure out nnz: my block
    
    my_nnz = 0;

    for (unsigned int ii = 0; ii < source->GetRowCount(); ii++)
    {
        unsigned int nSubStart = CKNMatrixOperation::pRow[ii], nSubEnd = CKNMatrixOperation::pRow[ii+1];
        unsigned int startColIndex = (int)mine->GetFirstRowIndex(), endColIndex = startColIndex + mine->GetColumnCount() - 1;

        for(unsigned int jj = nSubStart; jj < nSubEnd; jj++)
        {
            unsigned int nColIndex = CKNMatrixOperation::pColumn[jj];
            if(startColIndex <= nColIndex && nColIndex <= endColIndex)
            {
                my_nnz++;
                //curval.SetRealNumber(lpMatrixValueElement->pfReal[jj]); curval.SetImaginaryNumber(lpMatrixValueElement->pfImaginary[jj]);
                pData = source->GetValue(jj);
                curval.SetRealNumber(pData->GetRealNumber()); curval.SetImaginaryNumber(pData->GetImaginaryNumber());
                mine->SetAt(curval, ii, nColIndex-startColIndex);
            }               
        }
    }

//  if(CKNMPIManager::IsRootRank())
//      printf("My block update completed: my_nnz = %d (computed), %d (CSR-allocated)\n", my_nnz, mine->GetNoneZeroCount());
};

void CKNMatrixOperation::FreeLocalCSR(CKNMatrixOperation::CKNCSR *mine, CKNMatrixOperation::CKNCSR *left, CKNMatrixOperation::CKNCSR *right)
{
    if (mine != NULL)
    {
        CKNMatrixOperation::FreeCSR(mine);
        mine = NULL;
    }
    if (left != NULL)
    {
        CKNMatrixOperation::FreeCSR(left);
        left = NULL;
    }
    if (right != NULL)
    {
        CKNMatrixOperation::FreeCSR(right);
        right = NULL;
    }
}

bool CKNMatrixOperation::IsSame(CKNComplex operand1, CKNComplex operand2, double tol)
{
    if( fabs(fabs(operand1.GetRealNumber()) - fabs(operand2.GetRealNumber())) > tol )
        return false;

    if( fabs(fabs(operand1.GetImaginaryNumber()) - fabs(operand2.GetImaginaryNumber())) > tol )
        return false;
    else
        return true;
}

bool CKNMatrixOperation::IsSame(CKNVector *pVector1, CKNVector *pVector2)
{
    CKNVector       vectorTemp = *pVector1;
    double          fNorm = 1.0;

    vectorTemp.MinusVector(pVector2);
#ifndef DISABLE_MPI_ROUTINE
    fNorm = vectorTemp.GetNorm(true);
#else //DISABLE_MPI_ROUTINE
    fNorm = vectorTemp.GetNorm();
#endif //DISABLE_MPI_ROUTINE
    
    if (IsSame(fNorm, 0.0, GENERAL_TOLERANCE))
        return true;
    else
        return false;
}