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 blas_oct_m
  use blacs_proc_grid_oct_m
  use blacs_oct_m
  use debug_oct_m
  use global_oct_m
#ifndef NDEBUG
  use, intrinsic :: ieee_exceptions
#endif
  use, intrinsic :: iso_fortran_env
  use lapack_oct_m
  use lalg_basic_oct_m
  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_cholesky_pivoted,             &
    lalg_geneigensolve,           &
    lalg_eigensolve,              &
    lalg_eigensolve_tridiagonal,  &
    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_zdni,                    &
    lalg_zduialpha,               &
    lalg_zd2ni,                   &
    lalg_least_squares,           &
    lalg_matrix_rank_svd,         &
    lalg_remove_rotation,         &
    pseudoinverse_default_tolerance
 
 
 
  interface lalg_cholesky
    module procedure dcholesky, zcholesky
  end interface lalg_cholesky
 
  interface lalg_cholesky_pivoted
    module procedure dcholesky_pivoted, zcholesky_pivoted
  end interface lalg_cholesky_pivoted
 
  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_tridiagonal
    module procedure deigensolve_tridiagonal, zeigensolve_tridiagonal
  end interface lalg_eigensolve_tridiagonal
 
  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_function
    module procedure dlalg_matrix_function, zlalg_matrix_function
  end interface lalg_matrix_function
 
  interface lalg_matrix_rank_svd
    module procedure dmatrix_rank_svd, zmatrix_rank_svd
  end interface lalg_matrix_rank_svd
 
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)   :: work(:)
    real(real64),      intent(out)   :: rwork(:)
    integer,           intent(out)   :: info
 
    real(real64), allocatable :: wr(:), wi(:)
#ifndef NDEBUG
    logical :: halting_mode(size(ieee_all))
#endif
 
    push_sub(lalg_dgeev)
 
    safe_allocate(wr(1:n))
    safe_allocate(wi(1:n))
 
#ifndef NDEBUG
    ! Optimized LAPACK backends (e.g. MKL) may raise spurious floating-point
    ! exceptions while still returning correct results, so suspend trapping
    ! for builds that enable it
    call ieee_get_halting_mode(ieee_all, halting_mode)
    call ieee_set_halting_mode(ieee_all, .false.)
#endif
    call dgeev(jobvl, jobvr, n, a(1, 1), lda, wr(1), wi(1), vl(1, 1), ldvl, vr(1, 1), ldvr, work(1), lwork, info)
#ifndef NDEBUG
    call ieee_set_halting_mode(ieee_all, halting_mode)
#endif
    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(:, :) = matmul(evectors(:, :), ex(:, :))
 
      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(:, :) = matmul(ex(:, :), evectors(:, :))
 
      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(:, :) = matmul(evectors(:, :), ex(:, :))
 
      safe_deallocate_a(evalues)
      safe_deallocate_a(zevalues)
    else
      safe_allocate(zevalues(1:nn))
 
      evectors(:, :) = aa(:, :)
 
      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(:, :) = evectors(:, :)
 
      call lalg_inverse(nn, evectors, 'dir')
 
      do ii = 1, nn
        evectors(1:nn, ii) = zevalues(1:nn)*evectors(1:nn, ii)
      end do
 
      ex(:, :) = matmul(ex(:, :), evectors(:, :))
 
      safe_deallocate_a(zevalues)
    end if
 
    pop_sub(zlalg_phi)
  end subroutine zlalg_phi
 
  complex(real64) function lalg_zdni(eigenvec, alpha, beta)
    integer,         intent(in) :: alpha, beta
    complex(real64), intent(in) :: 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),intent(in) :: 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), intent(in) :: 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: