KNTBMS_Solver.cpp

Go to the documentation of this file.

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <algorithm>
#include <mpi.h>
#include <omp.h>

#include "CKNGlobal.h"
#include "KNLanczosMethod.h"
#include "KNIPCCUtility.h"
#include "KNMPIManager.h"
#include "KNTimeMeasurement.h"
#include "KNLanczosResultAudit.h"
#include "KNCommandFileParser.h"
#include "KNGeometricAtom.h"
#include "KNGeometricAtomFactory.h"
#include "KNMPIManager.h"
#include "KNTBMS_Solver.h"
#include "KNHamiltonianBuilder.h"

#include "XeonPhi_header.h"
#ifndef _WIN32
#include <unistd.h>
#endif // _WIN32
int phi_tid = 0;


void CKNTBMS_Solver::getCurrentTime(char *pszBuffer)
{
    time_t          timer;
    struct tm       *t;

    timer = time(NULL);
    t = localtime(&timer);

    sprintf(pszBuffer, "%d:%d:%d", t->tm_hour, t->tm_min, t->tm_sec);
}

void CKNTBMS_Solver::GetKValues(CKNCommandFileParser::LPINPUT_CMD_PARAM lpParam, double *pKValue[3])
{
    double                  fInterval[3];
    unsigned int            i, j;

    for( i = 0; i < 3; ++i)
    {
        double              *pfKValueUnit = NULL;
        pKValue[i] = (double*)malloc(sizeof(double)*(unsigned int)lpParam->fKPoints);
        pfKValueUnit = pKValue[i];
        fInterval[i] = (lpParam->fKValueFinal[i] - lpParam->fKValueInit[i]) / (lpParam->fKPoints - 1);
        pfKValueUnit[0] = lpParam->fKValueInit[i];
        if (lpParam->fKPoints > 1 )
            pfKValueUnit[(int)lpParam->fKPoints-1] = lpParam->fKValueFinal[i];
        for( j = 1; j < (unsigned int)lpParam->fKPoints-1 ; ++j )
            pfKValueUnit[j] = pfKValueUnit[j-1] + fInterval[i];
    }
}

bool CKNTBMS_Solver::InitMPIEnv(bool &bMPI, CKNCommandFileParser::LPINPUT_CMD_PARAM lpParam)
{
    int                         rank;
    int                         world_size;
    bool                        bRtn = true;

    if (!bMPI)
        return bRtn;

    MPI_Init(NULL, NULL);
    MPI_Comm_size(MPI_COMM_WORLD, &world_size);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);

    CKNMPIManager::SetMPIEnviroment(rank, world_size);
    bRtn = CKNMPIManager::InitLevel(lpParam->nMPILevel, lpParam->nFindingDegeneratedEVCount);

    return bRtn;
}

void CKNTBMS_Solver::FinalEvn(CKNMatrixOperation::CKNCSR *lpResult, CKNCommandFileParser::LPINPUT_CMD_PARAM lpParam, double *pKValue[3], NEIGHBOR_MAP_INFO *lpMapInfo, CKNGeometricShape *pGeometricShape, bool bMPI)
{
    int                     i;

    if (NULL != lpResult)
        CKNMatrixOperation::FreeCSR(lpResult);

    CKNGeometricShape::FreeMapInfo(lpMapInfo);
    FREE_MEM(lpParam);
    for (i = 0; i < 3; ++i)
        FREE_MEM(pKValue[i]);

    if (bMPI)
    {
        CKNMPIManager::FinalizeManager();
        MPI_Finalize();
    }
}

CKNMatrixOperation::CKNCSR* CKNTBMS_Solver::AllocateCSR(unsigned int &nRtn, NEIGHBOR_MAP_INFO *lpMapInfo, CKNCommandFileParser::LPINPUT_CMD_PARAM lpParam, CKNGeometricShape *pGeometricShape, bool bMPI)
{
    unsigned int                    nRowCount, nColumnCount;
    CKNMatrixOperation::CKNCSR      *pResult = NULL;

    pResult = new CKNMatrixOperation::CKNCSR();
    if (NULL == pResult)
        throw ERROR_MALLOC;

    nRowCount = (unsigned int)lpMapInfo->fItemCount * 10;
    nColumnCount = (unsigned int)lpMapInfo->fItemCount * 10;
    
    if (bMPI)
    {
        nColumnCount = (unsigned int)pGeometricShape->GetTotalAtomCount() * 10;
        CKNMPIManager::LoadBlancing(nRowCount);
        CKNMPIManager::InitCommunicationBufferMetric();
    }

    lpParam->nMatrixDemension = nColumnCount;
        
    pResult->SetRowCount(nRowCount);
    pResult->SetColumnCount(nColumnCount);
    pResult->BuildDataBuffer();
    pResult->SetFirstRowIndex(pGeometricShape->GetAtomStartID()*10);

    return pResult;
}

unsigned int CKNTBMS_Solver::Launching_TBMS_Solver(char *pszInputCommnadFileName, bool bMPI, bool bShowMsg, bool bRemovePrevResult)
{

    CKNGeometricShape                           GeometricShape;
    CKNMatrixOperation::CKNCSR                  *pHamiltonian = NULL;
    CKNMatrixOperation::CKNCSR                  *plocalH = NULL;
    CKNMatrixOperation::CKNCSR                  *plocalHm1 = NULL;
    CKNMatrixOperation::CKNCSR                  *plocalHp1 = NULL;
    CKNCommandFileParser::LPINPUT_CMD_PARAM     lpParam = NULL;
    NEIGHBOR_MAP_INFO                           mapInfo;
    unsigned int                                nRtnCode = CALCULATION_SUCCESS;
    unsigned int                                i, nRepeatCount = 1;
    CKNLanczosMethod::LPEIGENVALUE_RESULT       lpResult = NULL;
    CKNLanczosMethod                            lanczos;
    char                                        szMsg[1024];
    double                                      *pKValue[3] = {NULL, NULL, NULL};

    if (bShowMsg)
        CKNIPCCUtility::SetShow(true);

    lpParam = CKNCommandFileParser::ParsingInputCommand(pszInputCommnadFileName); 
    if (NULL == lpParam)
    {
        if (CKNMPIManager::IsRootRank())
        {
            sprintf(szMsg, "Terminated program. '%s' file doesn't existed\n", pszInputCommnadFileName);
            CKNIPCCUtility::ShowMsg(szMsg);
        }
        return CAN_NOT_FIND_COMMAND_FILE;
    }
    if(!CKNLanczosMethod::CheckingCalculationCondition(lpParam->bCalculateEigenVectors, lpParam->bCalculateWaveFunction, lpParam->nFindingDegeneratedEVCount))
    {
        SHOW_SIMPLE_MSG("Adjusting input parameter first!\n");
        return NEED_TO_ADJUST_MPI_CONFIG;
    }
    
    if( !InitMPIEnv(bMPI, lpParam)) 
    {
        SHOW_SIMPLE_MSG("Adjusting MPI level and node count first!\n");
        return NEED_TO_ADJUST_MPI_CONFIG;
    }

    if( CKNMPIManager::IsRootRank() )
    {
        if (bRemovePrevResult)
#ifdef _WIN32 
            system("rmdir /s /q result");
#else //WIN32
            system("rm -rf result");
#endif //WIN32

#pragma omp parallel
{
        if(omp_get_num_threads() > 1)
        {
            if(omp_get_thread_num() == 0)
            {
                sprintf(szMsg, "-Run with MPI (%d cores) and OpenMP (%d threads)\n", CKNMPIManager::GetTotalNodeCount(), omp_get_num_threads());
                CKNIPCCUtility::ShowMsg(szMsg);
            }
        }
        else
        {
            sprintf(szMsg, "-Run with MPI only (%d cores)\n", CKNMPIManager::GetTotalNodeCount());
            CKNIPCCUtility::ShowMsg(szMsg);
        }
}

        sprintf(szMsg, "-Command file: %s\n-Direction: [%03d], Periodic condition [%d %d %d]\n", pszInputCommnadFileName, lpParam->nDirectionSingle, lpParam->bConsiderBoundaryCondition[0] ? 1 : 0, lpParam->bConsiderBoundaryCondition[1] ? 1 : 0, lpParam->bConsiderBoundaryCondition[2] ? 1 : 0);
        CKNIPCCUtility::ShowMsg(szMsg);
        CKNIPCCUtility::ShowMsg("-Start Building Geometric construction\n");
    }

    CKNGeometricShape::SetShapeInformation(GeometricShape, lpParam); 
    if (!CKNGeometricShape::SetAtomAndNeighborInformation(lpParam)) 
    {
        SHOW_SIMPLE_MSG("Terminated program. Direction information doesn't existed.\n");
        nRtnCode = CAN_NOT_ALLOC_RANK;
        FREE_MEM(lpParam);
        return nRtnCode;
    }
    if(!GeometricShape.FillUnitcell(lpParam)) 
    {
        SHOW_SIMPLE_MSG("Terminated program. Node counts should be smaller than assigend unitcell count with X axis.\n");
        nRtnCode = CAN_NOT_ALLOC_RANK;
        FREE_MEM(lpParam);
        return nRtnCode;
    }

    GeometricShape.PeriodicUnitCellNumbering(); 
    if (bMPI)
    {
        GeometricShape.ExchangeAtomInfoBetweenNode(); 
        GeometricShape.PeriodicUnitCellNumbering(true); 
    }
    GeometricShape.BuildNeighborInformation(); 
    GeometricShape.ConstructMapInfo(&mapInfo, lpParam); 

    double      fSurfaceCount = (double)GeometricShape.GetSurfaceAtomList()->size();
    double      fTotalSurfaceCount = CKNMPIManager::AllReduceDouble(fSurfaceCount);

    if(CKNMPIManager::IsRootRank())
    {
        sprintf(szMsg, "-Surface atom counts: %d\n", (int)fTotalSurfaceCount);
        CKNIPCCUtility::ShowMsg(szMsg);
    }
    
    if (!strcmp(lpParam->szStructureType, CUBIC))
    {
        SHOW_SIMPLE_MSG("\n-Geometry construction completed!\n\n");

        GeometricShape.FinalShape();
        FinalEvn(pHamiltonian, lpParam, pKValue, &mapInfo, &GeometricShape, bMPI);

        return nRtnCode;
    }

    pHamiltonian = AllocateCSR(nRtnCode, &mapInfo, lpParam, &GeometricShape, bMPI); 
    
        
    SHOW_SIMPLE_MSG("-Building Hamiltonian Matrix\n");
    if (!CKNHamiltonianBuilder::BuildHamiltonian(lpParam, &mapInfo, pHamiltonian, GeometricShape.GetAtomStartID())) 
    {
        SHOW_SIMPLE_MSG("Terminated program. Can't load meterial param information.\n");
        nRtnCode = CAN_NOT_LOAD_METERIAL_PARAM;
        FinalEvn(pHamiltonian, lpParam, pKValue, &mapInfo, &GeometricShape, bMPI);
        return nRtnCode;
    }

    pHamiltonian->nComponentsFirstUnitCell = GeometricShape.m_nAtomFirstLayer*ORBITALS;
    pHamiltonian->nComponentsLastUnitCell = GeometricShape.m_nAtomLastLayer*ORBITALS;


    if( lpParam->bSaveHamiltonian )
    {
        SHOW_SIMPLE_MSG("-Save Hamiltonian matrix to file\n");
        CKNIPCCUtility::DumpCSR(pHamiltonian, "Hamiltonian.dat", GeometricShape.GetAtomStartID());
    }

#ifndef DISABLE_MPI_ROUTINE
    phi_tid = CKNMPIManager::GetCurrentRank() % 2;
// FIXME test code
#ifndef _WIN32
    char host_name[256];
    gethostname(host_name, 256);
    printf("-[Host:%s] (MPI:%03d) mapped to (MIC:%03d)\n", host_name, CKNMPIManager::GetCurrentRank(), phi_tid);
#endif //_WIN32
#else
    phi_tid = 0;
#endif

    double *local_real = NULL;
    double *local_imaginary = NULL;
    unsigned int *local_row = NULL;
    unsigned int *local_col = NULL;
    unsigned int local_size = 0;
    unsigned int local_row_size = 0;
    unsigned int local_col_size = 0;

    local_real = pHamiltonian->m_vectValueRealBuffer.data();
    local_imaginary = pHamiltonian->m_vectValueImaginaryBuffer.data();
    local_row = pHamiltonian->m_vectRow.data();
    local_col = pHamiltonian->m_vectColumn.data();
    local_size = pHamiltonian->GetNoneZeroCount();
    local_row_size = pHamiltonian->m_vectRow.size();
    local_col_size = pHamiltonian->m_vectColumn.size();
#pragma offload_transfer target(mic:phi_tid) \
    nocopy(local_real[0:local_size]      : ALLOC) \
    nocopy(local_imaginary[0:local_size] : ALLOC) \
    nocopy(local_row[0:local_row_size]   : ALLOC) \
    nocopy(local_col[0:local_col_size]   : ALLOC)

    double impurity_position[3];
    int tempnCell;
    tempnCell = floor(lpParam->fShapeLength[0][0] / lpParam->fUnitcellLength[0]);
    impurity_position[0] = (double)tempnCell / 2.0 * lpParam->fUnitcellLength[0];
    tempnCell = floor(lpParam->fShapeLength[0][1] / lpParam->fUnitcellLength[1]);
    impurity_position[1] = (double)tempnCell / 2.0 * lpParam->fUnitcellLength[1];
    tempnCell = floor(lpParam->fShapeLength[0][2] / lpParam->fUnitcellLength[2]);
    impurity_position[2] = (double)tempnCell / 2.0 * lpParam->fUnitcellLength[2];
    //ApplyPhPotential(impurity_position, pHamiltonian, mapInfo, true);
    
    
    GetKValues(lpParam, pKValue); 
    for (i = 0; i < (unsigned int)lpParam->fKPoints ; i++)
    {
        double              fKValue[3] = { 0., 0., 0. };
    
        fKValue[0] = pKValue[0][i]; fKValue[1] = pKValue[1][i]; fKValue[2] = pKValue[2][i];
        if(CKNMPIManager::IsRootRank())
        {
            sprintf(szMsg, "[%d] -Start Lanczos eigenvalue solver for Hamiltonian. (k = %f, %f, %f)\n", i, fKValue[0], fKValue[1], fKValue[2]);
            CKNIPCCUtility::ShowMsg(szMsg);
        }

        if (false == GeometricShape.RefillPeriodicBinding(pHamiltonian, lpParam, &mapInfo, i, fKValue)) 
            continue;

#pragma offload_transfer target(mic:phi_tid) \
    in(local_real[0:local_size]      : REUSE) \
    in(local_imaginary[0:local_size] : REUSE) \
    in(local_row[0:local_row_size]   : REUSE) \
    in(local_col[0:local_col_size]   : REUSE)

        CKNTimeMeasurement::InitTimer();
        CKNTimeMeasurement::TotalMeasurementStart();
        lpResult = lanczos.DoLanczosMethod(pHamiltonian,
                                            lpParam->nLanczosIterationCount,
                                            lpParam->nCheckEigenvalueInterval,
                                            lpParam->nFindingEigenValueCount,
                                            lpParam->fevMin,
                                            lpParam->fevMax,
                                            lpParam->fConvergeceCriteria,
                                            lpParam->bDoSelectiveReorthogonalization,
                                            lpParam->bCalculateEigenVectors,
                                            lpParam->bCalculateWaveFunction,
                                            lpParam->load_in_MIC,
                                            plocalH,   // ADDED
                                            plocalHm1, // ADDED
                                            plocalHp1);// ADDED

        if( CKNMPIManager::IsMultiLevelMPI() )
        {
            if(CKNMPIManager::IsRootRank())
                CKNIPCCUtility::ShowMsg("-Calculating degenerated eigenvalues\n");
            CKNLanczosMethod::MergeDegeneratedEigenvalues(lpResult, lpParam->nFindingDegeneratedEVCount, pHamiltonian, plocalH, plocalHm1, plocalHp1);
            if( lpParam->bCalculateWaveFunction )
                CKNLanczosMethod::RecalcuWaveFunction(lpResult);
        }

        if( CKNMPIManager::IsDeflationRoot() )
            CKNLanczosMethod::SortSolution(lpResult);
        
        if(CKNMPIManager::IsRootRank())
            CKNIPCCUtility::ShowMsg("\n-Save result into file\n");

        CKNLanczosMethod::SaveLanczosResult(lpResult, lpParam->bCalculateEigenVectors, lpParam->bCalculateWaveFunction, fKValue, i);
        CKNTimeMeasurement::TotalMeasurementEnd();

        CKNLanczosMethod::ShowLanczosResult(lpResult, lpParam->bCalculateEigenVectors, lpParam->bCalculateWaveFunction, fKValue, i);


        if(CKNMPIManager::IsRootRank())
        {
            sprintf(szMsg, "\n[%d] -Finished Lanczos eigenvalue solver for Hamiltonian. (k = %f, %f, %f)\n\n", i, fKValue[0], fKValue[1], fKValue[2]);
            CKNIPCCUtility::ShowMsg(szMsg);
        }
                
        CKNLanczosMethod::ReleaseResult(lpResult, true); 
        CKNHamiltonianBuilder::ResetPeriodicBinding(pHamiltonian, &mapInfo); 
    }
    CKNHamiltonianBuilder::FinalizeMatrixBuffer();

    SHOW_SIMPLE_MSG("- TBMS Solver completed!\n\n");


#pragma offload_transfer target(mic:phi_tid) \
    nocopy(local_real      : FREE) \
    nocopy(local_imaginary : FREE) \
    nocopy(local_row       : FREE) \
    nocopy(local_col       : FREE)

    FinalEvn(pHamiltonian, lpParam, pKValue, &mapInfo, &GeometricShape, bMPI);

    return nRtnCode;
}

void CKNTBMS_Solver::ApplyPhPotential(double *impurity_position, CKNMatrixOperation::CKNCSR *matrix, NEIGHBOR_MAP_INFO mapInfo, bool bUseSplitVector)
{
    CKNMatrixOperation::CKNVector potential;
    CKNComplex curval;
    potential.SetSize((unsigned int)mapInfo.fItemCount);
    double inveps = 1.0 / 11.9;
    double esqr = 1.4399766;
    double myx, myy, myz, distance;

    for (unsigned int i = 0; i < (unsigned int)mapInfo.fItemCount; i++)
    {
            myx = (mapInfo.pfX_Coordination[i] - impurity_position[0]);
            myy = (mapInfo.pfY_Coordination[i] - impurity_position[1]);
            myz = (mapInfo.pfZ_Coordination[i] - impurity_position[2]);
            
            distance = sqrt(myx*myx + myy*myy + myz*myz);
            
            if (distance < 1.0e-5)
                    curval.SetRealNumber(-3.782);
            else
                    curval.SetRealNumber(-inveps*esqr / distance);
            potential.SetAt(i, curval);
    }

    matrix->DiagonalOperation(&potential, CKNMatrixOperation::CKNCSR::PLUS, bUseSplitVector);

    if( CKNMPIManager::IsRootRank() )
        CKNIPCCUtility::ShowMsg("-Phospher potential plungged\n");
}