dftfe/linear_algebra_operations_8h_source.html

// ---------------------------------------------------------------------

//

// Copyright (c) 2017-2025  The Regents of the University of Michigan and DFT-FE

// authors.

//

// This file is part of the DFT-FE code.

//

// The DFT-FE code is free software; you can use it, redistribute

// it, and/or modify it under the terms of the GNU Lesser General

// Public License as published by the Free Software Foundation; either

// version 2.1 of the License, or (at your option) any later version.

// The full text of the license can be found in the file LICENSE at

// the top level of the DFT-FE distribution.

//

// ---------------------------------------------------------------------


#ifndef linearAlgebraOperations_h

#define linearAlgebraOperations_h


#include <elpaScalaManager.h>

#include <headers.h>

#include <operator.h>

#include "process_grid.h"

#include "scalapackWrapper.h"

#include "dftParameters.h"

#include <BLASWrapper.h>

namespace dftfe

{

  //

  // extern declarations for blas-lapack routines

  //

#ifndef DOXYGEN_SHOULD_SKIP_THIS

  extern "C"

  {

    void

    dgemv_(const char         *TRANS,

           const unsigned int *M,

           const unsigned int *N,

           const double       *alpha,

           const double       *A,

           const unsigned int *LDA,

           const double       *X,

           const unsigned int *INCX,

           const double       *beta,

           double             *C,

           const unsigned int *INCY);


    void

    sgemv_(const char         *TRANS,

           const unsigned int *M,

           const unsigned int *N,

           const float        *alpha,

           const float        *A,

           const unsigned int *LDA,

           const float        *X,

           const unsigned int *INCX,

           const float        *beta,

           float              *C,

           const unsigned int *INCY);


    void

    zgemv_(const char                 *TRANS,

           const unsigned int         *M,

           const unsigned int         *N,

           const std::complex<double> *alpha,

           const std::complex<double> *A,

           const unsigned int         *LDA,

           const std::complex<double> *X,

           const unsigned int         *INCX,

           const std::complex<double> *beta,

           std::complex<double>       *C,

           const unsigned int         *INCY);


    void

    cgemv_(const char                *TRANS,

           const unsigned int        *M,

           const unsigned int        *N,

           const std::complex<float> *alpha,

           const std::complex<float> *A,

           const unsigned int        *LDA,

           const std::complex<float> *X,

           const unsigned int        *INCX,

           const std::complex<float> *beta,

           std::complex<float>       *C,

           const unsigned int        *INCY);

    void

    dsymv_(const char         *UPLO,

           const unsigned int *N,

           const double       *alpha,

           const double       *A,

           const unsigned int *LDA,

           const double       *X,

           const unsigned int *INCX,

           const double       *beta,

           double             *C,

           const unsigned int *INCY);

    void

    dgesv_(int    *n,

           int    *nrhs,

           double *a,

           int    *lda,

           int    *ipiv,

           double *b,

           int    *ldb,

           int    *info);

    void

    dsysv_(const char *UPLO,

           const int  *n,

           const int  *nrhs,

           double     *a,

           const int  *lda,

           int        *ipiv,

           double     *b,

           const int  *ldb,

           double     *work,

           const int  *lwork,

           int        *info);

    void

    dscal_(const unsigned int *n,

           const double       *alpha,

           double             *x,

           const unsigned int *inc);

    void

    sscal_(const unsigned int *n,

           const float        *alpha,

           float              *x,

           const unsigned int *inc);

    void

    zscal_(const unsigned int         *n,

           const std::complex<double> *alpha,

           std::complex<double>       *x,

           const unsigned int         *inc);

    void

    zdscal_(const unsigned int   *n,

            const double         *alpha,

            std::complex<double> *x,

            const unsigned int   *inc);

    void

    daxpy_(const unsigned int *n,

           const double       *alpha,

           const double       *x,

           const unsigned int *incx,

           double             *y,

           const unsigned int *incy);

    void

    saxpy_(const unsigned int *n,

           const float        *alpha,

           const float        *x,

           const unsigned int *incx,

           float              *y,

           const unsigned int *incy);

    void

    dgemm_(const char         *transA,

           const char         *transB,

           const unsigned int *m,

           const unsigned int *n,

           const unsigned int *k,

           const double       *alpha,

           const double       *A,

           const unsigned int *lda,

           const double       *B,

           const unsigned int *ldb,

           const double       *beta,

           double             *C,

           const unsigned int *ldc);

    void

    sgemm_(const char         *transA,

           const char         *transB,

           const unsigned int *m,

           const unsigned int *n,

           const unsigned int *k,

           const float        *alpha,

           const float        *A,

           const unsigned int *lda,

           const float        *B,

           const unsigned int *ldb,

           const float        *beta,

           float              *C,

           const unsigned int *ldc);

    void

    dsyevd_(const char         *jobz,

            const char         *uplo,

            const unsigned int *n,

            double             *A,

            const unsigned int *lda,

            double             *w,

            double             *work,

            const unsigned int *lwork,

            int                *iwork,

            const unsigned int *liwork,

            int                *info);

    void

    dsygvx_(const int    *itype,

            const char   *jobz,

            const char   *range,

            const char   *uplo,

            const int    *n,

            double       *a,

            const int    *lda,

            double       *b,

            const int    *ldb,

            const double *vl,

            const double *vu,

            const int    *il,

            const int    *iu,

            const double *abstol,

            int          *m,

            double       *w,

            double       *z,

            const int    *ldz,

            double       *work,

            const int    *lwork,

            int          *iwork,

            int          *ifail,

            int          *info);

    void

    dsyevx_(const char   *jobz,

            const char   *range,

            const char   *uplo,

            const int    *n,

            double       *a,

            const int    *lda,

            const double *vl,

            const double *vu,

            const int    *il,

            const int    *iu,

            const double *abstol,

            int          *m,

            double       *w,

            double       *z,

            const int    *ldz,

            double       *work,

            const int    *lwork,

            int          *iwork,

            int          *ifail,

            int          *info);

    double

    dlamch_(const char *cmach);

    void

    dsyevr_(const char         *jobz,

            const char         *range,

            const char         *uplo,

            const unsigned int *n,

            double             *A,

            const unsigned int *lda,

            const double       *vl,

            const double       *vu,

            const unsigned int *il,

            const unsigned int *iu,

            const double       *abstol,

            const unsigned int *m,

            double             *w,

            double             *Z,

            const unsigned int *ldz,

            unsigned int       *isuppz,

            double             *work,

            const int          *lwork,

            int                *iwork,

            const int          *liwork,

            int                *info);

    void

    dsyrk_(const char         *uplo,

           const char         *trans,

           const unsigned int *n,

           const unsigned int *k,

           const double       *alpha,

           const double       *A,

           const unsigned int *lda,

           const double       *beta,

           double             *C,

           const unsigned int *ldc);

    void

    dsyr_(const char         *uplo,

          const unsigned int *n,

          const double       *alpha,

          const double       *X,

          const unsigned int *incx,

          double             *A,

          const unsigned int *lda);

    void

    dsyr2_(const char         *uplo,

           const unsigned int *n,

           const double       *alpha,

           const double       *x,

           const unsigned int *incx,

           const double       *y,

           const unsigned int *incy,

           double             *a,

           const unsigned int *lda);

    void

    dcopy_(const unsigned int *n,

           const double       *x,

           const unsigned int *incx,

           double             *y,

           const unsigned int *incy);

    void

    scopy_(const unsigned int *n,

           const float        *x,

           const unsigned int *incx,

           float              *y,

           const unsigned int *incy);

    void

    zgemm_(const char                 *transA,

           const char                 *transB,

           const unsigned int         *m,

           const unsigned int         *n,

           const unsigned int         *k,

           const std::complex<double> *alpha,

           const std::complex<double> *A,

           const unsigned int         *lda,

           const std::complex<double> *B,

           const unsigned int         *ldb,

           const std::complex<double> *beta,

           std::complex<double>       *C,

           const unsigned int         *ldc);

    void

    cgemm_(const char                *transA,

           const char                *transB,

           const unsigned int        *m,

           const unsigned int        *n,

           const unsigned int        *k,

           const std::complex<float> *alpha,

           const std::complex<float> *A,

           const unsigned int        *lda,

           const std::complex<float> *B,

           const unsigned int        *ldb,

           const std::complex<float> *beta,

           std::complex<float>       *C,

           const unsigned int        *ldc);

    void

    zheevd_(const char           *jobz,

            const char           *uplo,

            const unsigned int   *n,

            std::complex<double> *A,

            const unsigned int   *lda,

            double               *w,

            std::complex<double> *work,

            const unsigned int   *lwork,

            double               *rwork,

            const unsigned int   *lrwork,

            int                  *iwork,

            const unsigned int   *liwork,

            int                  *info);

    void

    zheevr_(const char           *jobz,

            const char           *range,

            const char           *uplo,

            const unsigned int   *n,

            std::complex<double> *A,

            const unsigned int   *lda,

            const double         *vl,

            const double         *vu,

            const unsigned int   *il,

            const unsigned int   *iu,

            const double         *abstol,

            const unsigned int   *m,

            double               *w,

            std::complex<double> *Z,

            const unsigned int   *ldz,

            unsigned int         *isuppz,

            std::complex<double> *work,

            const int            *lwork,

            double               *rwork,

            const int            *lrwork,

            int                  *iwork,

            const int            *liwork,

            int                  *info);

    void

    zherk_(const char                 *uplo,

           const char                 *trans,

           const unsigned int         *n,

           const unsigned int         *k,

           const double               *alpha,

           const std::complex<double> *A,

           const unsigned int         *lda,

           const double               *beta,

           std::complex<double>       *C,

           const unsigned int         *ldc);

    void

    zcopy_(const unsigned int         *n,

           const std::complex<double> *x,

           const unsigned int         *incx,

           std::complex<double>       *y,

           const unsigned int         *incy);


    void

    ccopy_(const unsigned int        *n,

           const std::complex<float> *x,

           const unsigned int        *incx,

           std::complex<float>       *y,

           const unsigned int        *incy);


    std::complex<double>

    zdotc_(const unsigned int         *N,

           const std::complex<double> *X,

           const unsigned int         *INCX,

           const std::complex<double> *Y,

           const unsigned int         *INCY);

    std::complex<float>

    cdotc_(const unsigned int        *N,

           const std::complex<float> *X,

           const unsigned int        *INCX,

           const std::complex<float> *Y,

           const unsigned int        *INCY);

    double

    ddot_(const unsigned int *N,

          const double       *X,

          const unsigned int *INCX,

          const double       *Y,

          const unsigned int *INCY);

    float

    sdot_(const unsigned int *N,

          const float        *X,

          const unsigned int *INCX,

          const float        *Y,

          const unsigned int *INCY);

    double

    dnrm2_(const unsigned int *n, const double *x, const unsigned int *incx);


    double

    dznrm2_(const unsigned int         *n,

            const std::complex<double> *x,

            const unsigned int         *incx);

    void

    zaxpy_(const unsigned int         *n,

           const std::complex<double> *alpha,

           const std::complex<double> *x,

           const unsigned int         *incx,

           std::complex<double>       *y,

           const unsigned int         *incy);

    void

    caxpy_(const unsigned int        *n,

           const std::complex<float> *alpha,

           const std::complex<float> *x,

           const unsigned int        *incx,

           std::complex<float>       *y,

           const unsigned int        *incy);

    void

    dpotrf_(const char         *uplo,

            const unsigned int *n,

            double             *a,

            const unsigned int *lda,

            int                *info);

    void

    dpotri_(const char         *uplo,

            const unsigned int *n,

            double             *A,

            const unsigned int *lda,

            int                *info);


    void

    zpotrf_(const char           *uplo,

            const unsigned int   *n,

            std::complex<double> *a,

            const unsigned int   *lda,

            int                  *info);

    void

    dtrtri_(const char         *uplo,

            const char         *diag,

            const unsigned int *n,

            double             *a,

            const unsigned int *lda,

            int                *info);

    void

    ztrtri_(const char           *uplo,

            const char           *diag,

            const unsigned int   *n,

            std::complex<double> *a,

            const unsigned int   *lda,

            int                  *info);


    // LU decomoposition of a general matrix

    void

    dgetrf_(int *M, int *N, double *A, int *lda, int *IPIV, int *INFO);


    // generate inverse of a matrix given its LU decomposition

    void

    dgetri_(int    *N,

            double *A,

            int    *lda,

            int    *IPIV,

            double *WORK,

            int    *lwork,

            int    *INFO);

    // LU decomoposition of a general matrix

    void

    zgetrf_(int                  *M,

            int                  *N,

            std::complex<double> *A,

            int                  *lda,

            int                  *IPIV,

            int                  *INFO);


    // generate inverse of a matrix given its LU decomposition

    void

    zgetri_(int                  *N,

            std::complex<double> *A,

            int                  *lda,

            int                  *IPIV,

            std::complex<double> *WORK,

            int                  *lwork,

            int                  *INFO);

  }

#endif


  /**

   *  @brief Contains linear algebra functions used in the implementation of an eigen solver

   *

   *  @author Phani Motamarri, Sambit Das

   */


  namespace linearAlgebraOperations

  {

    /** @brief Compute inverse of serial matrix using LAPACK LU factorization

     */

    void

    inverse(double *A, dftfe::Int N);


    /** @brief Compute inverse of serial matrix using LAPACK LU factorization

     */

    void

    inverse(std::complex<double> *A, dftfe::Int N);


    /** @brief Calculates an estimate of lower and upper bounds of a matrix using

     *  k-step Generalised Lanczos method. Algo is present in PAW PRB paper

     *

     *  @param  operatorMatrix An object which has access to the given matrix

     *  @param  vect A dummy vector

     *  @return std::pair<double,double> An estimate of the lower and upper bound of the given matrix

     */

    template <typename T, dftfe::utils::MemorySpace memorySpace>

    std::pair<double, double>

    generalisedLanczosLowerUpperBoundEigenSpectrum(

      const std::shared_ptr<dftfe::linearAlgebra::BLASWrapper<memorySpace>>

                                                        &BLASWrapperPtr,

      operatorDFTClass<memorySpace>                     &operatorMatrix,

      dftfe::linearAlgebra::MultiVector<T, memorySpace> &X,

      dftfe::linearAlgebra::MultiVector<T, memorySpace> &Y,

      dftfe::linearAlgebra::MultiVector<T, memorySpace> &Z,

      dftfe::linearAlgebra::MultiVector<T, memorySpace> &tempVec,

      const dftParameters                               &dftParams);


    /** @brief Apply Chebyshev filter to a given subspace

     *

     *  @param[in] operatorMatrix An object which has access to the given matrix

     *  @param[in,out]  X Given subspace as a dealii array representing multiple

     * fields as a flattened array. In-place update of the given subspace.

     *  @param[in]  numberComponents Number of multiple-fields

     *  @param[in]  m Chebyshev polynomial degree

     *  @param[in]  a lower bound of unwanted spectrum

     *  @param[in]  b upper bound of unwanted spectrum

     *  @param[in]  a0 lower bound of wanted spectrum

     */

    template <typename T, dftfe::utils::MemorySpace memorySpace>

    void

    chebyshevFilter(operatorDFTClass<memorySpace> &operatorMatrix,

                    dftfe::linearAlgebra::MultiVector<T, memorySpace> &X,

                    dftfe::linearAlgebra::MultiVector<T, memorySpace> &Y,

                    const dftfe::uInt                                  m,

                    const double                                       a,

                    const double                                       b,

                    const double                                       a0);


    /** @brief Apply Residual based Chebyshev filter to a given subspace

     *

     *  @param[in] operatorMatrix An object which has access to the given matrix

     *  @param[in,out]  X Given subspace as a dealii array representing multiple

     * fields as a flattened array. In-place update of the given subspace.

     *  @param[in]  eigenvalues estimate of eigenvalues, usually the eigenvalues

     * from previous pass

     *  @param[in]  m Chebyshev polynomial degree

     *  @param[in]  a lower bound of unwanted spectrum

     *  @param[in]  b upper bound of unwanted spectrum

     *  @param[in]  a0 lower bound of wanted spectrum

     *  @param[in]  approxOverlapMatrix to use approximate overlap matrix while

     * computing initial residual

     */


    template <typename T1, typename T2, dftfe::utils::MemorySpace memorySpace>

    void

    reformulatedChebyshevFilter(

      const std::shared_ptr<dftfe::linearAlgebra::BLASWrapper<memorySpace>>

                                                         &BLASWrapperPtr,

      operatorDFTClass<memorySpace>                      &operatorMatrix,

      dftfe::linearAlgebra::MultiVector<T1, memorySpace> &X,

      dftfe::linearAlgebra::MultiVector<T1, memorySpace> &Y,

      dftfe::linearAlgebra::MultiVector<T2, memorySpace> &Residual,

      dftfe::linearAlgebra::MultiVector<T2, memorySpace> &ResidualNew,

      std::vector<double>                                 eigenvalues,

      const dftfe::uInt                                   m,

      const double                                        a,

      const double                                        b,

      const double                                        a0,

      const bool                                          approxOverlapMatrix);


  } // namespace linearAlgebraOperations


} // namespace dftfe

#endif

BLASWrapper.h

dftfe::dftParameters
Namespace which declares the input parameters and the functions to parse them from the input paramete...
Definition dftParameters.h:36

dftfe::linearAlgebra::BLASWrapper
Definition BLASWrapper.h:35

dftfe::linearAlgebra::MultiVector
An class template to encapsulate a MultiVector. A MultiVector is a collection of  vectors belonging t...
Definition MultiVector.h:127

dftfe::operatorDFTClass
Base class for building the DFT operator and the action of operator on a vector.
Definition operator.h:43

dftParameters.h

elpaScalaManager.h

headers.h

dftfe::dftUtils::dgesv_
void dgesv_(int *N, int *NRHS, double *A, int *LDA, int *IPIV, double *B, int *LDB, int *INFO)

dftfe::linearAlgebraOperations
Contains linear algebra functions used in the implementation of an eigen solver.
Definition linearAlgebraOperations.h:513

dftfe::linearAlgebraOperations::inverse
void inverse(double *A, dftfe::Int N)
Compute inverse of serial matrix using LAPACK LU factorization.

dftfe::linearAlgebraOperations::chebyshevFilter
void chebyshevFilter(operatorDFTClass< memorySpace > &operatorMatrix, dftfe::linearAlgebra::MultiVector< T, memorySpace > &X, dftfe::linearAlgebra::MultiVector< T, memorySpace > &Y, const dftfe::uInt m, const double a, const double b, const double a0)
Apply Chebyshev filter to a given subspace.

dftfe::linearAlgebraOperations::reformulatedChebyshevFilter
void reformulatedChebyshevFilter(const std::shared_ptr< dftfe::linearAlgebra::BLASWrapper< memorySpace > > &BLASWrapperPtr, operatorDFTClass< memorySpace > &operatorMatrix, dftfe::linearAlgebra::MultiVector< T1, memorySpace > &X, dftfe::linearAlgebra::MultiVector< T1, memorySpace > &Y, dftfe::linearAlgebra::MultiVector< T2, memorySpace > &Residual, dftfe::linearAlgebra::MultiVector< T2, memorySpace > &ResidualNew, std::vector< double > eigenvalues, const dftfe::uInt m, const double a, const double b, const double a0, const bool approxOverlapMatrix)
Apply Residual based Chebyshev filter to a given subspace.

dftfe::linearAlgebraOperations::generalisedLanczosLowerUpperBoundEigenSpectrum
std::pair< double, double > generalisedLanczosLowerUpperBoundEigenSpectrum(const std::shared_ptr< dftfe::linearAlgebra::BLASWrapper< memorySpace > > &BLASWrapperPtr, operatorDFTClass< memorySpace > &operatorMatrix, dftfe::linearAlgebra::MultiVector< T, memorySpace > &X, dftfe::linearAlgebra::MultiVector< T, memorySpace > &Y, dftfe::linearAlgebra::MultiVector< T, memorySpace > &Z, dftfe::linearAlgebra::MultiVector< T, memorySpace > &tempVec, const dftParameters &dftParams)
Calculates an estimate of lower and upper bounds of a matrix using k-step Generalised Lanczos method....

dftfe
Definition pseudoPotentialToDftfeConverter.cc:34

dftfe::HubbardOccFieldType::Residual
@ Residual
Definition hubbardClass.h:59

dftfe::uInt
std::uint32_t uInt
Definition TypeConfig.h:10

dftfe::ExcFamilyType::LDA
@ LDA
Definition ExcSSDFunctionalBaseClass.h:35

dftfe::Int
std::int32_t Int
Definition TypeConfig.h:11

operator.h

process_grid.h

scalapackWrapper.h