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"

#ifdef WITH_SYMPACK
#include <sympack.hpp>
#include "pexsi/sympack_interf.hpp"
#endif

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_;           

  // Used for performing "CopyPattern"
  DistSparseMatrix<Real>     PatternMat_;

  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


  // Below specifically for the case when H and S are Hermitian
  DistSparseMatrix<Complex>  HComplexMat_;
  DistSparseMatrix<Complex>  SComplexMat_;

  DistSparseMatrix<Complex>     rhoComplexMat_;               // Density matrix 
  DistSparseMatrix<Complex>     rhoDrvMuComplexMat_;          // Derivative of the Fermi-Dirac with respect to mu
  DistSparseMatrix<Complex>     rhoDrvTComplexMat_;           // Derivative of the Fermi-Dirac with respect to T
  DistSparseMatrix<Complex>     freeEnergyDensityComplexMat_; // Helmholtz free energy density matrix
  DistSparseMatrix<Complex>     energyDensityComplexMat_;     // 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_;

#ifdef WITH_SYMPACK
  symPACK::symPACKMatrix<Real>*      symPACKRealMat_;
  symPACK::symPACKMatrix<Complex>*   symPACKComplexMat_;
  symPACK::symPACKOptions             symPACKOpt_;

  // Used for performing "CopyPattern"
  symPACK::DistSparseMatrix<Real>     symmPatternMat_;

  symPACK::DistSparseMatrix<Real>     symmHRealMat_;
  symPACK::DistSparseMatrix<Real>     symmSRealMat_;
  // Below specifically for the case when H and S are Hermitian
  symPACK::DistSparseMatrix<Complex>  symmHComplexMat_;
  symPACK::DistSparseMatrix<Complex>  symmSComplexMat_;

  symPACK::DistSparseMatrix<Real>     symmShiftRealMat_;
  symPACK::DistSparseMatrix<Complex>  symmShiftComplexMat_;

  symPACK::DistSparseMatrix<Real>     symmShiftInvRealMat_;
  symPACK::DistSparseMatrix<Complex>  symmShiftInvComplexMat_;
#endif

  SuperLUOptions             luOpt_;
  FactorizationOptions       factOpt_;
  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 LoadRealMatrix(
      Int           nrows,                        
      Int           nnz,                          
      Int           nnzLocal,                     
      Int           numColLocal,                  
      Int*          colptrLocal,                  
      Int*          rowindLocal,                  
      Real*         HnzvalLocal,                  
      Int           isSIdentity,                  
      Real*         SnzvalLocal,
    Int               solver,
      Int           verbosity );


  void LoadComplexMatrix(
      Int           nrows,                        
      Int           nnz,                          
      Int           nnzLocal,                     
      Int           numColLocal,                  
      Int*          colptrLocal,                  
      Int*          rowindLocal,                  
      Complex*      HnzvalLocal,                  
      Int           isSIdentity,                  
      Complex*      SnzvalLocal,
    Int               solver,
      Int           verbosity );


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

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



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

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



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

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


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

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



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

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



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


  void CalculateFermiOperatorComplex(
      Int   numPole, 
      Real  temperature,
      Real  gap,
      Real  deltaE,
      Real  mu,
      Real  numElectronExact, 
      Real  numElectronTolerance,
    Int               solver,
      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               solver,
      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               solver,
      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               solver,
      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<Complex>&   RhoComplexMat() const {return rhoComplexMat_;}

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

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

  Real   TotalEnergyH() const {return totalEnergyH_;}

  Real   TotalEnergyS() const {return totalEnergyS_;}

  Real   TotalFreeEnergy() const {return totalFreeEnergy_;}


}; // PPEXSIData


} // namespace PEXSI
#endif // _PPEXSI_HPP_