ppexsi.hpp

Go to the documentation of this file.

/*
   Copyright (c) 2012 The Regents of the University of California,
   through Lawrence Berkeley National Laboratory.  

   Author: Lin Lin

   This file is part of PEXSI. All rights reserved.

   Redistribution and use in source and binary forms, with or without
   modification, are permitted provided that the following conditions are met:

   (1) Redistributions of source code must retain the above copyright notice, this
   list of conditions and the following disclaimer.
   (2) Redistributions in binary form must reproduce the above copyright notice,
   this list of conditions and the following disclaimer in the documentation
   and/or other materials provided with the distribution.
   (3) Neither the name of the University of California, Lawrence Berkeley
   National Laboratory, U.S. Dept. of Energy nor the names of its contributors may
   be used to endorse or promote products derived from this software without
   specific prior written permission.

   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
   ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
   WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
   DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
   ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
   (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
   LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
   ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

   You are under no obligation whatsoever to provide any bug fixes, patches, or
   upgrades to the features, functionality or performance of the source code
   ("Enhancements") to anyone; however, if you choose to make your Enhancements
   available either publicly, or directly to Lawrence Berkeley National
   Laboratory, without imposing a separate written license agreement for such
   Enhancements, then you hereby grant the following license: a non-exclusive,
   royalty-free perpetual license to install, use, modify, prepare derivative
   works, incorporate into other computer software, distribute, and sublicense
   such enhancements or derivative works thereof, in binary and source code form.
*/
#ifndef _PPEXSI_HPP_
#define _PPEXSI_HPP_
#include "pexsi/environment.hpp"
#include "pexsi/sparse_matrix.hpp"
#include "pexsi/NumVec.hpp"
#include "pexsi/utility.hpp"
#include "pexsi/pole.hpp"
#include "pexsi/mpi_interf.hpp"
#include "pexsi/SuperLUGrid.hpp"
#include "pexsi/superlu_dist_interf.hpp"
#include "pexsi/pselinv.hpp"
#include        "pexsi/pselinv_unsym.hpp"
//#include "pexsi/ngchol_interf.hpp"
//#include "pexsi/c_pexsi_interface.h"

namespace PEXSI{

        class PPEXSIData{
        private:
                // *********************************************************************
                // Computational variables
                // *********************************************************************

                std::vector<Complex>  zshift_;      // Complex shift for the pole expansion
                std::vector<Complex>  zweightRho_;  // Complex weight for the pole expansion for density
                std::vector<Complex>  zweightRhoDrvMu_;  // Complex weight for the pole expansion for derivative of the Fermi-Dirac with respect to the chemical potential
                std::vector<Complex>  zweightRhoDrvT_;   // Complex weight for the pole expansion for derivative of the Fermi-Dirac with respect to the temperature T (1/beta, in au)
                std::vector<Complex>  zweightHelmholtz_;  // Complex shift for the pole expansion for Helmholtz free energy
                std::vector<Complex>  zweightForce_;  // Complex weight for the pole expansion for force

    // Outer layer communicator. Also used for distributing the
    // DistSparseMatrix.  Each DistSparseMatrix is replicated in the row
    // (numPoleGroup) direction of gridPole.
                const GridType*           gridPole_;          
                const GridType*           gridSelInv_;        // Inner layer communicator for SelInv
        
    // Inner layer communicator for SuperLU factorization
        const SuperLUGrid<Real>*       gridSuperLUReal_;           
        const SuperLUGrid<Complex>*    gridSuperLUComplex_;           


    DistSparseMatrix<Real>     HRealMat_;
    DistSparseMatrix<Real>     SRealMat_;

    DistSparseMatrix<Real>     shiftRealMat_;
    DistSparseMatrix<Complex>  shiftComplexMat_;

    DistSparseMatrix<Real>     shiftInvRealMat_;
    DistSparseMatrix<Complex>  shiftInvComplexMat_;

                DistSparseMatrix<Real>     rhoRealMat_;                   // Density matrix 
          DistSparseMatrix<Real>     rhoDrvMuRealMat_;              // Derivative of the Fermi-Dirac with respect to mu
                DistSparseMatrix<Real>     rhoDrvTRealMat_;               // Derivative of the Fermi-Dirac with respect to T
                DistSparseMatrix<Real>     freeEnergyDensityRealMat_;     // Helmholtz free energy density matrix
                DistSparseMatrix<Real>     energyDensityRealMat_;         // Energy density matrix for computing the Pulay force

    // SuperLUMatrix and PMatrix structures These structures are saved
    // to avoid repetitive symbolic factorization process, and saved in
    // pointer form because of the constructors.
    SuperLUMatrix<Real>*       luRealMat_;
    SuperLUMatrix<Complex>*    luComplexMat_;

    SuperLUOptions             luOpt_;
    PSelInvOptions             selinvOpt_;

    PMatrix<Real>*             PMRealMat_;
    PMatrix<Complex>*          PMComplexMat_;
    PMatrixUnsym<Real>*        PMRealUnsymMat_;
    PMatrixUnsym<Complex>*     PMComplexUnsymMat_;

    // Whether the matrices have been loaded into HRealMat_ and
    // SRealMat_
    bool                       isMatrixLoaded_;
    // Whether the matrices (luMat and PMat) have obtained symbolic
    // information
    bool                       isRealSymmetricSymbolicFactorized_;
    bool                       isComplexSymmetricSymbolicFactorized_;
    bool                       isRealUnsymmetricSymbolicFactorized_;
    bool                       isComplexUnsymmetricSymbolicFactorized_;
    // Supernode partition for the real matrix
    SuperNodeType              superReal_;             
    // Supernode partition for the complex matrix
    SuperNodeType              superComplex_;             

                // Saves all the indices of diagonal elements in H, so that
                // H.nzvalLocal(diagIdxLocal_[j]) are diagonal elements for all j.
                // This is manly used when S is implicitly given as an identity matrix.
                std::vector<Int>           diagIdxLocal_;    

    // Energy computed from Tr[H*DM]
    Real                       totalEnergyH_;
    // Energy computed from Tr[S*EDM]
    Real                       totalEnergyS_;
    // Free energy 
    Real                       totalFreeEnergy_;


                // *********************************************************************
                // Saved variables for nonlinear iterations
                // *********************************************************************

        public:
    PPEXSIData(
        MPI_Comm   comm,
        Int        numProcRow, 
        Int        numProcCol, 
        Int        outputFileIndex );

                ~PPEXSIData();

    void LoadRealSymmetricMatrix(
        Int           nrows,                        
        Int           nnz,                          
        Int           nnzLocal,                     
        Int           numColLocal,                  
        Int*          colptrLocal,                  
        Int*          rowindLocal,                  
        Real*         HnzvalLocal,                  
        Int           isSIdentity,                  
        Real*         SnzvalLocal,
        Int           verbosity );

    void LoadRealUnsymmetricMatrix(
        Int           nrows,                        
        Int           nnz,                          
        Int           nnzLocal,                     
        Int           numColLocal,                  
        Int*          colptrLocal,                  
        Int*          rowindLocal,                  
        Real*         HnzvalLocal,                  
        Int           isSIdentity,                  
        Real*         SnzvalLocal,
        Int           verbosity );


    
    void SymbolicFactorizeRealSymmetricMatrix(
                                std::string                    ColPerm,
                                Int                            numProcSymbFact,
        Int                            verbosity );
 
    void SymbolicFactorizeRealUnsymmetricMatrix(
                                std::string                    ColPerm,
        std::string                    RowPerm,
                                Int                            numProcSymbFact,
        Int                            Transpose,
        double*                        AnzvalLocal,                  
        Int                            verbosity );

   
    void SymbolicFactorizeComplexSymmetricMatrix(
                                std::string                    ColPerm,
                                Int                            numProcSymbFact,
        Int                            verbosity );

    void SymbolicFactorizeComplexUnsymmetricMatrix(
                                std::string                    ColPerm,
        std::string                    RowPerm,
                                Int                            numProcSymbFact,
        Int                            Transpose,
        double*                        AnzvalLocal,                  
        Int                            verbosity );



    void SelInvRealSymmetricMatrix(
          double*           AnzvalLocal,                  
          Int               verbosity,
          double*           AinvnzvalLocal );

    void SelInvRealUnsymmetricMatrix(
          double*           AnzvalLocal,                  
          Int               verbosity,
          double*           AinvnzvalLocal );


    void SelInvComplexSymmetricMatrix(
          double*           AnzvalLocal,                  
          Int               verbosity,
          double*           AinvnzvalLocal );

    void SelInvComplexUnsymmetricMatrix(
          double*           AnzvalLocal,                  
          Int               verbosity,
          double*           AinvnzvalLocal );



                void CalculateNegativeInertiaReal(
                                const std::vector<Real>&       shiftVec, 
                                std::vector<Real>&             inertiaVec,
        Int                            verbosity );


                void CalculateFermiOperatorReal(
                                Int   numPole, 
                                Real  temperature,
                                Real  gap,
                                Real  deltaE,
                                Real  mu,
        Real  numElectronExact, 
        Real  numElectronTolerance,
        Int   verbosity,
        Real& numElectron,
        Real& numElectronDrvMu );


    void DFTDriver(
        Real       numElectronExact,
        Real       temperature,
        Real       gap,
        Real       deltaE,
        Int        numPole, 
        Int        isInertiaCount,
        Int        maxPEXSIIter,
        Real       muMin0,
        Real       muMax0,
        Real       mu0,
        Real       muInertiaTolerance,
        Real       muInertiaExpansion,
        Real       muPEXSISafeGuard,
        Real       numElectronPEXSITolerance,
        Int        matrixType,
        Int        isSymbolicFactorize,
        Int        ordering,
        Int        numProcSymbFact,
        Int        verbosity,
        Real&      muPEXSI,                   
        Real&      numElectronPEXSI,         
        Real&      muMinInertia,              
        Real&      muMaxInertia,             
        Int&       numTotalInertiaIter,   
        Int&       numTotalPEXSIIter );



                void CalculateFermiOperatorReal2(
                                Int   numPole, 
                                Real  temperature,
                                Real  gap,
                                Real  deltaE,
        Real  numElectronExact, 
        Real  numElectronTolerance,
        Real  muMinPEXSI,
        Real  muMaxPEXSI,
        Int   verbosity,
                                Real& mu,
        Real& numElectron, 
        bool& isPEXSIConverged );

    void DFTDriver2(
        Real       numElectronExact,
        Real       temperature,
        Real       gap,
        Real       deltaE,
        Int        numPole, 
        Int        isInertiaCount,
        Real       muMin0,
        Real       muMax0,
        Real       mu0,
        Real       muInertiaTolerance,
        Real       muInertiaExpansion,
        Real       numElectronPEXSITolerance,
        Int        matrixType,
        Int        isSymbolicFactorize,
        Int        ordering,
        Int        numProcSymbFact,
        Int        verbosity,
        Real&      muPEXSI,                   
        Real&      numElectronPEXSI,         
        Real&      muMinInertia,              
        Real&      muMaxInertia,             
        Int&       numTotalInertiaIter );




    // *********************************************************************
    // Access data
    // *********************************************************************

    const GridType*  GridPole() const {return gridPole_;}

    const DistSparseMatrix<Real>&   RhoRealMat() const {return rhoRealMat_;}

    const DistSparseMatrix<Real>&   EnergyDensityRealMat() const {return energyDensityRealMat_;}

    const DistSparseMatrix<Real>&   FreeEnergyDensityRealMat() const {return freeEnergyDensityRealMat_;}

    Real   TotalEnergyH() const {return totalEnergyH_;}
    
    Real   TotalEnergyS() const {return totalEnergyS_;}

    Real   TotalFreeEnergy() const {return totalFreeEnergy_;}


        }; // PPEXSIData


} // namespace PEXSI
#endif // _PPEXSI_HPP_