lalg_adv.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module lalg_adv_oct_m
  use blacs_proc_grid_oct_m
  use blacs_oct_m
  use debug_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use lapack_oct_m
  use lalg_basic_oct_m, only: lalg_gemm_cc, lalg_scal
  use math_oct_m
  use messages_oct_m
  use mpi_oct_m
  use profiling_oct_m
  use scalapack_oct_m
  use sort_oct_m
  use utils_oct_m
#ifdef HAVE_ELPA
  use elpa
#endif
 
  implicit none
 
  private
  public ::                       &
    lalg_cholesky,                &
    lalg_geneigensolve,           &
    lalg_eigensolve,              &
    lalg_eigensolve_nonh,         &
    lalg_eigensolve_parallel,     &
    lalg_determinant,             &
    lalg_linsyssolve,             &
    lalg_singular_value_decomp,   &
    lalg_inverse,                 &
    lalg_svd_inverse,             &
    lalg_pseudo_inverse,          &
    lalg_lowest_geneigensolve,    &
    lalg_lowest_eigensolve,       &
    zlalg_exp,                    &
    zlalg_phi,                    &
    zlalg_matrix_function,        &
    lalg_matrix_function,         &
    lalg_zeigenderivatives,       &
    lalg_check_zeigenderivatives, &
    lalg_zdni,                    &
    lalg_zduialpha,               &
    lalg_zd2ni,                   &
    lalg_least_squares,           &
    lalg_matrix_norm2,            &
    lalg_remove_rotation
 
 
 
  interface lalg_cholesky
    module procedure dcholesky, zcholesky
  end interface lalg_cholesky
 
  interface lalg_geneigensolve
    module procedure dgeneigensolve, zgeneigensolve
  end interface lalg_geneigensolve
 
  interface lalg_eigensolve_nonh
    module procedure zeigensolve_nonh, deigensolve_nonh
  end interface lalg_eigensolve_nonh
 
  interface lalg_eigensolve
    module procedure deigensolve, zeigensolve
  end interface lalg_eigensolve
 
  interface lalg_eigensolve_parallel
    module procedure deigensolve_parallel, zeigensolve_parallel
  end interface lalg_eigensolve_parallel
 
 
  interface lalg_determinant
    module procedure ddeterminant, zdeterminant
  end interface lalg_determinant
 
  interface lalg_inverse
    module procedure dinverse, zinverse
  end interface lalg_inverse
 
  interface lalg_linsyssolve
    module procedure dlinsyssolve, zlinsyssolve
  end interface lalg_linsyssolve
 
  interface lalg_singular_value_decomp
    module procedure dsingular_value_decomp, zsingular_value_decomp
  end interface lalg_singular_value_decomp
 
  interface lalg_pseudo_inverse
    module procedure dlalg_pseudo_inverse, zlalg_pseudo_inverse
  end interface lalg_pseudo_inverse
 
  interface lalg_svd_inverse
    module procedure dsvd_inverse, zsvd_inverse
  end interface lalg_svd_inverse
 
  interface lalg_lowest_geneigensolve
    module procedure dlowest_geneigensolve, zlowest_geneigensolve
  end interface lalg_lowest_geneigensolve
 
  interface lalg_lowest_eigensolve
    module procedure dlowest_eigensolve, zlowest_eigensolve
  end interface lalg_lowest_eigensolve
 
  interface lapack_geev
    module procedure lalg_dgeev, lalg_zgeev
  end interface lapack_geev
 
  interface lalg_least_squares
    module procedure dleast_squares_vec, zleast_squares_vec
  end interface lalg_least_squares
 
  interface lalg_matrix_norm2
    module procedure dmatrix_norm2, zmatrix_norm2
  end interface lalg_matrix_norm2
 
  interface lalg_matrix_function
    module procedure dlalg_matrix_function, zlalg_matrix_function
  end interface lalg_matrix_function
 
contains
 
  ! ---------------------------------------------------------
  real(real64) function sfmin()
    interface
      real(real64) function dlamch(cmach)
        import real64
        implicit none
        character(1), intent(in) :: cmach
      end function dlamch
    end interface
 
    sfmin = dlamch('S')
  end function sfmin
 
  subroutine lalg_dgeev(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info)
    character(1),      intent(in)    :: jobvl, jobvr
    integer,           intent(in)    :: n, lda, ldvl, ldvr, lwork
    real(real64),      intent(inout) :: a(:,:)
    complex(real64),   intent(out)   :: w(:)
    real(real64),      intent(out)   :: vl(:,:), vr(:,:)
    real(real64),      intent(out)   :: rwork(:)
    real(real64),      intent(out)   :: work(:)
    integer,           intent(out)   :: info
 
    real(real64), allocatable :: wr(:), wi(:)
 
    push_sub(lalg_dgeev)
 
    safe_allocate(wr(1:n))
    safe_allocate(wi(1:n))
 
    call dgeev(jobvl, jobvr, n, a(1, 1), lda, wr(1), wi(1), vl(1, 1), ldvl, vr(1, 1), ldvr, work(1), lwork, rwork(1), info)
 
    w(1:n) = cmplx(wr(1:n), wi(1:n), real64)
 
    safe_deallocate_a(wr)
    safe_deallocate_a(wi)
 
    pop_sub(lalg_dgeev)
  end subroutine lalg_dgeev
 
  subroutine lalg_zgeev(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info)
    character(1),      intent(in)    :: jobvl, jobvr
    integer,           intent(in)    :: n, lda, ldvl, ldvr, lwork
    complex(real64),   intent(inout) :: a(:,:)
    complex(real64),   intent(out)   :: w(:)
    complex(real64),   intent(out)   :: vl(:,:), vr(:,:)
    real(real64),      intent(out)   :: rwork(:)
    complex(real64),   intent(out)   :: work(:)
    integer,           intent(out)   :: info
 
    push_sub(lalg_zgeev)
 
    call zgeev(jobvl, jobvr, n, a(1, 1), lda, w(1), vl(1, 1), ldvl, vr(1, 1), ldvr, work(1), lwork, rwork(1), info)
 
    pop_sub(lalg_zgeev)
  end subroutine lalg_zgeev
 
  subroutine zlalg_exp(nn, pp, aa, ex, hermitian)
    integer,           intent(in)      :: nn
    complex(real64),   intent(in)      :: pp
    complex(real64),   intent(in)      :: aa(:, :)
    complex(real64),   intent(inout)   :: ex(:, :)
    logical,           intent(in)      :: hermitian
 
    complex(real64), allocatable :: evectors(:, :), zevalues(:)
    real(real64), allocatable :: evalues(:)
 
    integer :: ii
 
    push_sub(zlalg_exp)
 
    safe_allocate(evectors(1:nn, 1:nn))
 
    if (hermitian) then
      safe_allocate(evalues(1:nn))
      safe_allocate(zevalues(1:nn))
 
      evectors(1:nn, 1:nn) = aa(1:nn, 1:nn)
 
      call lalg_eigensolve(nn, evectors, evalues)
 
      zevalues(1:nn) = exp(pp*evalues(1:nn))
 
      do ii = 1, nn
        ex(1:nn, ii) = zevalues(1:nn)*conjg(evectors(ii, 1:nn))
      end do
 
      ex(1:nn, 1:nn) = matmul(evectors(1:nn, 1:nn), ex(1:nn, 1:nn))
 
      safe_deallocate_a(evalues)
      safe_deallocate_a(zevalues)
    else
      safe_allocate(zevalues(1:nn))
 
      evectors(1:nn, 1:nn) = aa(1:nn, 1:nn)
 
      call lalg_eigensolve_nonh(nn, evectors, zevalues)
 
      zevalues(1:nn) = exp(pp*zevalues(1:nn))
 
      ex(1:nn, 1:nn) = evectors(1:nn, 1:nn)
 
      call lalg_inverse(nn, evectors, 'dir')
 
      do ii = 1, nn
        evectors(1:nn, ii) = zevalues(1:nn)*evectors(1:nn, ii)
      end do
 
      ex(1:nn, 1:nn) = matmul(ex(1:nn, 1:nn), evectors(1:nn, 1:nn))
 
      safe_deallocate_a(zevalues)
    end if
 
    safe_deallocate_a(evectors)
 
    pop_sub(zlalg_exp)
  end subroutine zlalg_exp
 
 
  subroutine zlalg_phi(nn, pp, aa, ex, hermitian)
    integer,           intent(in)      :: nn
    complex(real64),   intent(in)      :: pp
    complex(real64),   intent(in)      :: aa(:, :)
    complex(real64),   intent(inout)   :: ex(:, :)
    logical,           intent(in)      :: hermitian
 
    complex(real64), allocatable :: evectors(:, :), zevalues(:)
    real(real64), allocatable :: evalues(:)
 
    integer :: ii
 
    push_sub(zlalg_phi)
 
    safe_allocate(evectors(1:nn, 1:nn))
 
    if (hermitian) then
      safe_allocate(evalues(1:nn))
      safe_allocate(zevalues(1:nn))
 
      evectors = aa
 
      call lalg_eigensolve(nn, evectors, evalues)
 
      do ii = 1, nn
        zevalues(ii) = (exp(pp*evalues(ii)) - m_z1) / (pp*evalues(ii))
      end do
 
      do ii = 1, nn
        ex(1:nn, ii) = zevalues(1:nn)*conjg(evectors(ii, 1:nn))
      end do
 
      ex(1:nn, 1:nn) = matmul(evectors(1:nn, 1:nn), ex(1:nn, 1:nn))
 
      safe_deallocate_a(evalues)
      safe_deallocate_a(zevalues)
    else
      safe_allocate(zevalues(1:nn))
 
      evectors(1:nn, 1:nn) = aa(1:nn, 1:nn)
 
      call lalg_eigensolve_nonh(nn, evectors, zevalues)
 
      do ii = 1, nn
        zevalues(ii) = (exp(pp*zevalues(ii)) - m_z1) / (pp*zevalues(ii))
      end do
 
      ex(1:nn, 1:nn) = evectors(1:nn, 1:nn)
 
      call lalg_inverse(nn, evectors, 'dir')
 
      do ii = 1, nn
        evectors(1:nn, ii) = zevalues(1:nn)*evectors(1:nn, ii)
      end do
 
      ex(1:nn, 1:nn) = matmul(ex(1:nn, 1:nn), evectors(1:nn, 1:nn))
 
      safe_deallocate_a(zevalues)
    end if
 
    pop_sub(zlalg_phi)
  end subroutine zlalg_phi
 
 
  subroutine lalg_zeigenderivatives(n, mat, zeigenvec, zeigenval, zmat)
    integer, intent(in) :: n
    complex(real64), intent(in) :: mat(:, :)
    complex(real64), intent(out) :: zeigenvec(:, :)
    complex(real64), intent(out) :: zeigenval(:)
    complex(real64), intent(out) :: zmat(:, :, :)
 
    integer :: i, alpha, beta
    complex(real64), allocatable :: knaught(:, :)
    complex(real64), allocatable :: lambdaminusdm(:, :)
    complex(real64), allocatable :: ilambdaminusdm(:, :)
    complex(real64), allocatable :: unit(:, :)
 
    push_sub(lalg_zeigenderivatives)
 
    safe_allocate(unit(1:n, 1:n))
    safe_allocate(knaught(1:n, 1:n))
    safe_allocate(lambdaminusdm(1:n, 1:n))
    safe_allocate(ilambdaminusdm(1:n, 1:n))
 
    zeigenvec = mat
    call lalg_eigensolve_nonh(n, zeigenvec, zeigenval)
 
    unit = m_z0
    do i = 1, n
      unit(i, i) = m_z1
    end do
 
    do i = 1, n
      do alpha = 1, n
        do beta = 1, n
          knaught(alpha, beta) = - zeigenvec(alpha, i)*conjg(zeigenvec(beta, i))
        end do
      end do
      knaught = knaught + unit
      lambdaminusdm = zeigenval(i)*unit - mat
      ! TODO(Alex) Swap this call
      call lalg_zpseudoinverse(n, lambdaminusdm, ilambdaminusdm)
      zmat(:, :, i) = matmul(ilambdaminusdm, knaught)
    end do
 
    safe_deallocate_a(unit)
    safe_deallocate_a(knaught)
    safe_deallocate_a(lambdaminusdm)
    safe_deallocate_a(ilambdaminusdm)
    pop_sub(lalg_zeigenderivatives)
  end subroutine lalg_zeigenderivatives
 
 
  ! TODO(Alex) Issue #1132.
  ! * This routine is redunant. Its usage should be replaced with either zsvd_inverse or zlalg_pseudo_inverse.
  subroutine lalg_zpseudoinverse(n, mat, imat)
    integer, intent(in) :: n
    complex(real64), intent(in) :: mat(:, :)
    complex(real64), intent(out) :: imat(:, :)
 
    integer :: i
    complex(real64), allocatable :: u(:, :), vt(:, :), sigma(:, :)
    real(real64), allocatable :: sg_values(:)
 
    push_sub(lalg_zpseudoinverse)
 
    safe_allocate(u(1:n, 1:n))
    safe_allocate(vt(1:n, 1:n))
    safe_allocate(sigma(1:n, 1:n))
    safe_allocate(sg_values(1:n))
 
    imat = mat
    call lalg_singular_value_decomp(n, n, imat, u, vt, sg_values)
 
    sigma = m_z0
    do i = 1, n
      if (abs(sg_values(i)) <= 1.0e-12_real64 * maxval(abs(sg_values))) then
        sigma(i, i) = m_z0
      else
        sigma(i, i) = m_z1 / sg_values(i)
      end if
    end do
 
    vt = conjg(transpose(vt))
    u = conjg(transpose(u))
    imat = matmul(vt, matmul(sigma, u))
 
    ! Check if we truly have a pseudoinverse
    vt = matmul(mat, matmul(imat, mat)) - mat
    if (maxval(abs(vt)) > 1.0e-10_real64 * maxval(abs(mat))) then
      write(*, *) maxval(abs(vt))
      write(*, *) vt
      write(*, *)
      write(*, *) 1.0e-10_real64 * maxval(abs(mat))
      write(*, *) maxval(abs(vt)) > 1.0e-10_real64 * maxval(abs(mat))
      write(*, *) mat
      write(message(1), '(a)') 'Pseudoinverse failed.'
      call messages_fatal(1)
    end if
 
    safe_deallocate_a(u)
    safe_deallocate_a(vt)
    safe_deallocate_a(sigma)
    safe_deallocate_a(sg_values)
    pop_sub(lalg_zpseudoinverse)
  end subroutine lalg_zpseudoinverse
 
 
  subroutine lalg_check_zeigenderivatives(n, mat)
    integer, intent(in) :: n
    complex(real64), intent(in) :: mat(:, :)
 
    integer :: alpha, beta, gamma, delta
    complex(real64) :: deltah, zuder_direct, zder_direct, zuder_directplus, zuder_directminus
    complex(real64), allocatable :: zeigenvec(:, :), dm(:, :), zeigref_(:, :), zeigenval(:), mmatrix(:, :, :)
    complex(real64), allocatable :: zeigplus(:), zeigminus(:), zeig0(:), zeigplusminus(:), zeigminusplus(:)
 
    push_sub(lalg_check_zeigenderivatives)
 
    safe_allocate(zeigenvec(1:n, 1:n))
    safe_allocate(dm(1:n, 1:n))
    safe_allocate(zeigref_(1:n, 1:n))
    safe_allocate(zeigenval(1:n))
    safe_allocate(mmatrix(1:n, 1:n, 1:n))
    safe_allocate(zeigplus(1:n))
    safe_allocate(zeigminus(1:n))
    safe_allocate(zeig0(1:n))
    safe_allocate(zeigplusminus(1:n))
    safe_allocate(zeigminusplus(1:n))
 
    assert(n == 2)
 
    dm = mat
    call lalg_zeigenderivatives(2, dm, zeigref_, zeigenval, mmatrix)
 
 
    deltah = (0.000001_real64, 0.000_real64)
    !deltah = M_z1 * maxval(abs(dm)) * 0.001_real64
    do alpha = 1, n
      do beta = 1, n
        zder_direct = lalg_zdni(zeigref_(:, 2), alpha, beta)
        zuder_direct = lalg_zduialpha(zeigref_(:, 2), mmatrix(:, :, 2), 2, alpha, beta)
 
        zeigenvec = dm
        zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
        call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigplus)
        zeigenvec(:, 1) = zeigenvec(:, 1)/sum(conjg(zeigref_(1:2, 1))*zeigenvec(1:2, 1))
        zeigenvec(:, 2) = zeigenvec(:, 2)/sum(conjg(zeigref_(1:2, 2))*zeigenvec(1:2, 2))
        zuder_directplus = zeigenvec(2, 2)
 
        zeigenvec = dm
        zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
        call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigminus)
        zeigenvec(:, 1) = zeigenvec(:, 1)/sum(conjg(zeigref_(1:2, 1))*zeigenvec(1:2, 1))
        zeigenvec(:, 2) = zeigenvec(:, 2)/sum(conjg(zeigref_(1:2, 2))*zeigenvec(1:2, 2))
        zuder_directminus = zeigenvec(2, 2)
 
        write(*, '(2i1,4f24.12)') alpha, beta, zder_direct, (zeigplus(2) - zeigminus(2))/(m_two * deltah)
        write(*, '(2i1,4f24.12)') alpha, beta, &
          zuder_direct, (zuder_directplus - zuder_directminus) / (m_two * deltah)
 
        do gamma = 1, n
          do delta = 1, n
            if (alpha == gamma .and. beta == delta) then
              zder_direct = lalg_zd2ni(zeigref_(:, 1), mmatrix(:, :, 1), alpha, beta, gamma, delta)
 
              zeigenvec = dm
              call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeig0)
 
              zeigenvec = dm
              zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
              call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigplus)
 
              zeigenvec = dm
              zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
              call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigminus)
 
              write(*, '(4i1,4f24.12)') alpha, beta, gamma, delta, &
                zder_direct, (zeigplus(1) + zeigminus(1) - m_two*zeig0(1))/(deltah**2)
            else
              zder_direct = lalg_zd2ni(zeigref_(:, 1), mmatrix(:, :, 1), alpha, beta, gamma, delta)
 
              zeigenvec = dm
              zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
              zeigenvec(gamma, delta) = zeigenvec(gamma, delta) + deltah
              call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigplus)
 
              zeigenvec = dm
              zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
              zeigenvec(gamma, delta) = zeigenvec(gamma, delta) - deltah
              call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigminus)
 
              zeigenvec = dm
              zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
              zeigenvec(gamma, delta) = zeigenvec(gamma, delta) - deltah
              call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigplusminus)
 
              zeigenvec = dm
              zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
              zeigenvec(gamma, delta) = zeigenvec(gamma, delta) + deltah
              call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigminusplus)
 
              write(*, '(4i1,4f24.12)') alpha, beta, gamma, delta, &
                zder_direct, (zeigplus(1) + zeigminus(1) - zeigplusminus(1) - zeigminusplus(1))/(m_four*deltah**2)
 
 
            end if
          end do
        end do
 
      end do
    end do
 
    safe_deallocate_a(zeigenval)
    safe_deallocate_a(mmatrix)
    safe_deallocate_a(zeigenvec)
    safe_deallocate_a(dm)
    safe_deallocate_a(zeigref_)
    safe_deallocate_a(zeigplus)
    safe_deallocate_a(zeigminus)
    safe_deallocate_a(zeig0)
    safe_deallocate_a(zeigplusminus)
    safe_deallocate_a(zeigminusplus)
    pop_sub(lalg_check_zeigenderivatives)
  end subroutine lalg_check_zeigenderivatives
 
  complex(real64) function lalg_zdni(eigenvec, alpha, beta)
    integer, intent(in) :: alpha, beta
    complex(real64) :: eigenvec(2)
    lalg_zdni = conjg(eigenvec(alpha)) * eigenvec(beta)
  end function lalg_zdni
 
  complex(real64) function lalg_zduialpha(eigenvec, mmatrix, alpha, gamma, delta)
    integer, intent(in) :: alpha, gamma, delta
    complex(real64) :: eigenvec(2), mmatrix(2, 2)
    lalg_zduialpha = mmatrix(alpha, gamma) * eigenvec(delta)
  end function lalg_zduialpha
 
  complex(real64) function lalg_zd2ni(eigenvec, mmatrix, alpha, beta, gamma, delta)
    integer, intent(in) :: alpha, beta, gamma, delta
    complex(real64) :: eigenvec(2), mmatrix(2, 2)
    lalg_zd2ni  = conjg(mmatrix(alpha, delta) * eigenvec(gamma)) * eigenvec(beta) + &
      conjg(eigenvec(alpha)) * mmatrix(beta, gamma) * eigenvec(delta)
  end function lalg_zd2ni
 
  ! ---------------------------------------------------------
  pure real(real64) function pseudoinverse_default_tolerance(m, n, sg_values) result(tol)
    integer, intent(in) :: m, n
    real(real64),   intent(in) :: sg_values(:)
    tol = m_epsilon * m * n * maxval(sg_values)
  end function pseudoinverse_default_tolerance
 
 
  function lalg_remove_rotation(n, A) result(P)
    integer,      intent(in) :: n
    real(real64), intent(in) :: a(1:n, 1:n)
    real(real64) :: p(1:n, 1:n)
 
    real(real64) :: ata(1:n, 1:n)
 
    push_sub(lalg_remove_rotation)
 
    ata = matmul(transpose(a), a)
    call lalg_matrix_function(n, m_one, ata, p, square_root, .true.)
 
    pop_sub(lalg_remove_rotation)
  end function lalg_remove_rotation
 
#include "undef.F90"
#include "complex.F90"
#include "lalg_adv_lapack_inc.F90"
 
#include "undef.F90"
#include "real.F90"
#include "lalg_adv_lapack_inc.F90"
 
end module lalg_adv_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: