27 use,
intrinsic :: iso_fortran_env
151 real(real64) function sfmin()
153 real(real64) function dlamch(cmach)
156 character(1),
intent(in) :: cmach
163 subroutine lalg_dgeev(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info)
164 character(1),
intent(in) :: jobvl, jobvr
165 integer,
intent(in) :: n, lda, ldvl, ldvr, lwork
166 real(real64),
intent(inout) :: a(:, :)
167 complex(real64),
intent(out) :: w(:)
168 real(real64),
intent(out) :: vl(:, :), vr(:, :)
169 real(real64),
intent(out) :: rwork(:)
170 real(real64),
intent(out) :: work(:)
171 integer,
intent(out) :: info
173 real(real64),
allocatable :: wr(:), wi(:)
177 safe_allocate(wr(1:n))
178 safe_allocate(wi(1:n))
180 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)
182 w(1:n) = cmplx(wr(1:n), wi(1:n), real64)
184 safe_deallocate_a(wr)
185 safe_deallocate_a(wi)
190 subroutine lalg_zgeev(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info)
191 character(1),
intent(in) :: jobvl, jobvr
192 integer,
intent(in) :: n, lda, ldvl, ldvr, lwork
193 complex(real64),
intent(inout) :: a(:, :)
194 complex(real64),
intent(out) :: w(:)
195 complex(real64),
intent(out) :: vl(:, :), vr(:, :)
196 real(real64),
intent(out) :: rwork(:)
197 complex(real64),
intent(out) :: work(:)
198 integer,
intent(out) :: info
202 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)
224 subroutine zlalg_exp(nn, pp, aa, ex, hermitian)
225 integer,
intent(in) :: nn
226 complex(real64),
intent(in) :: pp
227 complex(real64),
intent(in) :: aa(:, :)
228 complex(real64),
intent(inout) :: ex(:, :)
229 logical,
intent(in) :: hermitian
231 complex(real64),
allocatable :: evectors(:, :), zevalues(:)
232 real(real64),
allocatable :: evalues(:)
238 safe_allocate(evectors(1:nn, 1:nn))
241 safe_allocate(evalues(1:nn))
242 safe_allocate(zevalues(1:nn))
244 evectors(1:nn, 1:nn) = aa(1:nn, 1:nn)
248 zevalues(1:nn) =
exp(pp*evalues(1:nn))
251 ex(1:nn, ii) = zevalues(1:nn)*conjg(evectors(ii, 1:nn))
254 ex(:, :) = matmul(evectors(:, :), ex(:, :))
256 safe_deallocate_a(evalues)
257 safe_deallocate_a(zevalues)
259 safe_allocate(zevalues(1:nn))
261 evectors(1:nn, 1:nn) = aa(1:nn, 1:nn)
265 zevalues(1:nn) =
exp(pp*zevalues(1:nn))
267 ex(1:nn, 1:nn) = evectors(1:nn, 1:nn)
272 evectors(1:nn, ii) = zevalues(1:nn)*evectors(1:nn, ii)
275 ex(:, :) = matmul(ex(:, :), evectors(:, :))
277 safe_deallocate_a(zevalues)
280 safe_deallocate_a(evectors)
302 subroutine zlalg_phi(nn, pp, aa, ex, hermitian)
303 integer,
intent(in) :: nn
304 complex(real64),
intent(in) :: pp
305 complex(real64),
intent(in) :: aa(:, :)
306 complex(real64),
intent(inout) :: ex(:, :)
307 logical,
intent(in) :: hermitian
309 complex(real64),
allocatable :: evectors(:, :), zevalues(:)
310 real(real64),
allocatable :: evalues(:)
316 safe_allocate(evectors(1:nn, 1:nn))
319 safe_allocate(evalues(1:nn))
320 safe_allocate(zevalues(1:nn))
322 evectors(:, :) = aa(:, :)
327 zevalues(ii) = (
exp(pp*evalues(ii)) - m_z1) / (pp*evalues(ii))
331 ex(1:nn, ii) = zevalues(1:nn)*conjg(evectors(ii, 1:nn))
334 ex(:, :) = matmul(evectors(:, :), ex(:, :))
336 safe_deallocate_a(evalues)
337 safe_deallocate_a(zevalues)
339 safe_allocate(zevalues(1:nn))
341 evectors(:, :) = aa(:, :)
346 zevalues(ii) = (
exp(pp*zevalues(ii)) - m_z1) / (pp*zevalues(ii))
349 ex(:, :) = evectors(:, :)
354 evectors(1:nn, ii) = zevalues(1:nn)*evectors(1:nn, ii)
357 ex(:, :) = matmul(ex(:, :), evectors(:, :))
359 safe_deallocate_a(zevalues)
378 integer,
intent(in) :: n
379 complex(real64),
contiguous,
intent(in) :: mat(:, :)
380 complex(real64),
contiguous,
intent(out) :: zeigenvec(:, :)
381 complex(real64),
contiguous,
intent(out) :: zeigenval(:)
382 complex(real64),
contiguous,
intent(out) :: zmat(:, :, :)
384 integer :: i, alpha, beta
385 complex(real64),
allocatable :: knaught(:, :)
386 complex(real64),
allocatable :: lambdaminusdm(:, :)
387 complex(real64),
allocatable :: ilambdaminusdm(:, :)
388 complex(real64),
allocatable :: unit(:, :)
392 safe_allocate(unit(1:n, 1:n))
393 safe_allocate(knaught(1:n, 1:n))
394 safe_allocate(lambdaminusdm(1:n, 1:n))
395 safe_allocate(ilambdaminusdm(1:n, 1:n))
408 knaught(alpha, beta) = - zeigenvec(alpha, i)*conjg(zeigenvec(beta, i))
411 knaught = knaught + unit
412 lambdaminusdm = zeigenval(i)*unit - mat
415 zmat(:, :, i) = matmul(ilambdaminusdm, knaught)
418 safe_deallocate_a(unit)
419 safe_deallocate_a(knaught)
420 safe_deallocate_a(lambdaminusdm)
421 safe_deallocate_a(ilambdaminusdm)
430 integer,
intent(in) :: n
431 complex(real64),
contiguous,
intent(in) :: mat(:, :)
432 complex(real64),
contiguous,
intent(out) :: imat(:, :)
435 complex(real64),
allocatable :: u(:, :), vt(:, :), sigma(:, :)
436 real(real64),
allocatable :: sg_values(:)
440 safe_allocate(u(1:n, 1:n))
441 safe_allocate(vt(1:n, 1:n))
442 safe_allocate(sigma(1:n, 1:n))
443 safe_allocate(sg_values(1:n))
450 if (abs(sg_values(i)) <= 1.0e-12_real64 * maxval(abs(sg_values)))
then
453 sigma(i, i) = m_z1 / sg_values(i)
457 vt(:, :) = conjg(transpose(vt(:, :)))
458 u(:, :) = conjg(transpose(u(:, :)))
459 imat = matmul(vt, matmul(sigma, u))
462 vt(:, :) = matmul(mat(:, :), matmul(imat(:, :), mat(:, :))) - mat(:, :)
463 if (maxval(abs(vt)) > 1.0e-10_real64 * maxval(abs(mat)))
then
464 write(*, *) maxval(abs(vt))
467 write(*, *) 1.0e-10_real64 * maxval(abs(mat))
468 write(*, *) maxval(abs(vt)) > 1.0e-10_real64 * maxval(abs(mat))
470 write(message(1),
'(a)')
'Pseudoinverse failed.'
471 call messages_fatal(1)
475 safe_deallocate_a(vt)
476 safe_deallocate_a(sigma)
477 safe_deallocate_a(sg_values)
489 integer,
intent(in) :: n
490 complex(real64),
intent(in) :: mat(:, :)
492 integer :: alpha, beta, gamma, delta
493 complex(real64) :: deltah, zuder_direct, zder_direct, zuder_directplus, zuder_directminus
494 complex(real64),
allocatable :: zeigenvec(:, :), dm(:, :), zeigref_(:, :), zeigenval(:), mmatrix(:, :, :)
495 complex(real64),
allocatable :: zeigplus(:), zeigminus(:), zeig0(:), zeigplusminus(:), zeigminusplus(:)
499 safe_allocate(zeigenvec(1:n, 1:n))
500 safe_allocate(dm(1:n, 1:n))
501 safe_allocate(zeigref_(1:n, 1:n))
502 safe_allocate(zeigenval(1:n))
503 safe_allocate(mmatrix(1:n, 1:n, 1:n))
504 safe_allocate(zeigplus(1:n))
505 safe_allocate(zeigminus(1:n))
506 safe_allocate(zeig0(1:n))
507 safe_allocate(zeigplusminus(1:n))
508 safe_allocate(zeigminusplus(1:n))
516 deltah = (0.000001_real64, 0.000_real64)
520 zder_direct =
lalg_zdni(zeigref_(:, 2), alpha, beta)
521 zuder_direct =
lalg_zduialpha(zeigref_(:, 2), mmatrix(:, :, 2), 2, alpha, beta)
523 zeigenvec(:, :) = dm(:, :)
524 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
526 zeigenvec(:, 1) = zeigenvec(:, 1)/sum(conjg(zeigref_(1:2, 1))*zeigenvec(1:2, 1))
527 zeigenvec(:, 2) = zeigenvec(:, 2)/sum(conjg(zeigref_(1:2, 2))*zeigenvec(1:2, 2))
528 zuder_directplus = zeigenvec(2, 2)
530 zeigenvec(:, :) = dm(:, :)
531 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
533 zeigenvec(:, 1) = zeigenvec(:, 1)/sum(conjg(zeigref_(1:2, 1))*zeigenvec(1:2, 1))
534 zeigenvec(:, 2) = zeigenvec(:, 2)/sum(conjg(zeigref_(1:2, 2))*zeigenvec(1:2, 2))
535 zuder_directminus = zeigenvec(2, 2)
537 write(*,
'(2i1,4f24.12)') alpha, beta, zder_direct, (zeigplus(2) - zeigminus(2))/(m_two * deltah)
538 write(*,
'(2i1,4f24.12)') alpha, beta, &
539 zuder_direct, (zuder_directplus - zuder_directminus) / (m_two * deltah)
543 if (alpha == gamma .and. beta == delta)
then
544 zder_direct =
lalg_zd2ni(zeigref_(:, 1), mmatrix(:, :, 1), alpha, beta, gamma, delta)
546 zeigenvec(:, :) = dm(:, :)
549 zeigenvec(:, :) = dm(:, :)
550 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
553 zeigenvec(:, :) = dm(:, :)
554 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
557 write(*,
'(4i1,4f24.12)') alpha, beta, gamma, delta, &
558 zder_direct, (zeigplus(1) + zeigminus(1) - m_two*zeig0(1))/(deltah**2)
560 zder_direct =
lalg_zd2ni(zeigref_(:, 1), mmatrix(:, :, 1), alpha, beta, gamma, delta)
562 zeigenvec(:, :) = dm(:, :)
563 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
564 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) + deltah
567 zeigenvec(:, :) = dm(:, :)
568 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
569 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) - deltah
572 zeigenvec(:, :) = dm(:, :)
573 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
574 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) - deltah
577 zeigenvec(:, :) = dm(:, :)
578 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
579 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) + deltah
582 write(*,
'(4i1,4f24.12)') alpha, beta, gamma, delta, &
583 zder_direct, (zeigplus(1) + zeigminus(1) - zeigplusminus(1) - zeigminusplus(1))/(m_four*deltah**2)
593 safe_deallocate_a(zeigenval)
594 safe_deallocate_a(mmatrix)
595 safe_deallocate_a(zeigenvec)
596 safe_deallocate_a(dm)
597 safe_deallocate_a(zeigref_)
598 safe_deallocate_a(zeigplus)
599 safe_deallocate_a(zeigminus)
600 safe_deallocate_a(zeig0)
601 safe_deallocate_a(zeigplusminus)
602 safe_deallocate_a(zeigminusplus)
606 complex(real64) function lalg_zdni(eigenvec, alpha, beta)
607 integer,
intent(in) :: alpha, beta
608 complex(real64) :: eigenvec(2)
609 lalg_zdni = conjg(eigenvec(alpha)) * eigenvec(beta)
612 complex(real64) function lalg_zduialpha(eigenvec, mmatrix, alpha, gamma, delta)
613 integer,
intent(in) :: alpha, gamma, delta
614 complex(real64) :: eigenvec(2), mmatrix(2, 2)
618 complex(real64) function lalg_zd2ni(eigenvec, mmatrix, alpha, beta, gamma, delta)
619 integer,
intent(in) :: alpha, beta, gamma, delta
620 complex(real64) :: eigenvec(2), mmatrix(2, 2)
621 lalg_zd2ni = conjg(mmatrix(alpha, delta) * eigenvec(gamma)) * eigenvec(beta) + &
622 conjg(eigenvec(alpha)) * mmatrix(beta, gamma) * eigenvec(delta)
631 integer,
intent(in) :: m, n
632 real(real64),
intent(in) :: sg_values(:)
633 tol = m_epsilon * m * n * maxval(sg_values)
643 integer,
intent(in) :: n
644 real(real64),
intent(in) :: a(1:n, 1:n)
645 real(real64) :: p(1:n, 1:n)
647 real(real64) :: ata(1:n, 1:n)
651 ata = matmul(transpose(a), a)
658#include "complex.F90"
659#include "lalg_adv_lapack_inc.F90"
663#include "lalg_adv_lapack_inc.F90"
Note that lalg_determinant and lalg_inverse are just wrappers over the same routine.
double exp(double __x) __attribute__((__nothrow__
This module contains interfaces for BLACS routines Interfaces are from http:
This module provides the BLACS processor grid.
This module contains interfaces for BLAS routines You should not use these routines directly....
complex(real64) function, public lalg_zd2ni(eigenvec, mmatrix, alpha, beta, gamma, delta)
integer function zmatrix_rank_svd(a, preserve_mat, tol)
Compute the rank of the matrix A using SVD.
subroutine zlowest_geneigensolve(k, n, a, b, e, v, preserve_mat, bof, err_code)
Computes the k lowest eigenvalues and the eigenvectors of a real symmetric or complex Hermitian gener...
complex(real64) function zdeterminant(n, a, preserve_mat)
Invert a real symmetric or complex Hermitian square matrix a.
subroutine zeigensolve_nonh(n, a, e, err_code, side, sort_eigenvectors)
Computes all the eigenvalues and the right (left) eigenvectors of a real or complex (non-Hermitian) e...
subroutine zlinsyssolve(n, nrhs, a, b, x)
compute the solution to a complex system of linear equations A*X = B, where A is an N-by-N matrix and...
subroutine, public zlalg_exp(nn, pp, aa, ex, hermitian)
subroutine deigensolve_parallel(n, a, e, bof, err_code)
Computes all the eigenvalues and the eigenvectors of a real symmetric or complex Hermitian eigenprobl...
subroutine dgeneigensolve(n, a, b, e, preserve_mat, bof, err_code)
Computes all the eigenvalues and the eigenvectors of a real symmetric or complex Hermitian generalize...
subroutine lalg_zpseudoinverse(n, mat, imat)
Computes the Moore-Penrose pseudoinverse of a complex matrix.
subroutine zlalg_pseudo_inverse(a, threshold)
Invert a matrix with the Moore-Penrose pseudo-inverse.
subroutine zeigensolve(n, a, e, bof, err_code)
Computes all eigenvalues and eigenvectors of a real symmetric or hermitian square matrix A.
subroutine dsingular_value_decomp(m, n, a, u, vt, sg_values, preserve_mat)
Computes the singular value decomposition of a real M x N matrix a.
subroutine deigensolve_nonh(n, a, e, err_code, side, sort_eigenvectors)
Computes all the eigenvalues and the right (left) eigenvectors of a real or complex (non-Hermitian) e...
subroutine deigensolve_tridiagonal(n, a, e, bof, err_code)
Computes all eigenvalues and eigenvectors of a real symmetric tridiagonal matrix. For the Hermitian c...
subroutine zinverse(n, a, method, det, threshold, uplo)
An interface to different method to invert a matrix.
real(real64) function sfmin()
Auxiliary function.
subroutine dlinsyssolve(n, nrhs, a, b, x)
compute the solution to a real system of linear equations A*X = B, where A is an N-by-N matrix and X ...
subroutine zsingular_value_decomp(m, n, a, u, vt, sg_values, preserve_mat)
Computes the singular value decomposition of a complex MxN matrix a.
subroutine dcholesky(n, a, bof, err_code)
Compute the Cholesky decomposition of real symmetric or complex Hermitian positive definite matrix a,...
subroutine dlalg_matrix_function(n, factor, a, fun_a, fun, hermitian, tridiagonal)
This routine calculates a function of a matrix by using an eigenvalue decomposition.
subroutine deigensolve(n, a, e, bof, err_code)
Computes all eigenvalues and eigenvectors of a real symmetric or hermitian square matrix A.
subroutine, public lalg_zeigenderivatives(n, mat, zeigenvec, zeigenval, zmat)
subroutine zeigensolve_tridiagonal(n, a, e, bof, err_code)
Computes all eigenvalues and eigenvectors of a real symmetric tridiagonal matrix. For the Hermitian c...
subroutine zlowest_eigensolve(k, n, a, e, v, preserve_mat)
Computes the k lowest eigenvalues and the eigenvectors of a standard symmetric-definite eigenproblem,...
complex(real64) function, public lalg_zduialpha(eigenvec, mmatrix, alpha, gamma, delta)
subroutine zgeneigensolve(n, a, b, e, preserve_mat, bof, err_code)
Computes all the eigenvalues and the eigenvectors of a real symmetric or complex Hermitian generalize...
subroutine dleast_squares_vec(nn, aa, bb, xx, preserve_mat)
subroutine dlowest_geneigensolve(k, n, a, b, e, v, preserve_mat, bof, err_code)
Computes the k lowest eigenvalues and the eigenvectors of a real symmetric or complex Hermitian gener...
subroutine zcholesky(n, a, bof, err_code)
Compute the Cholesky decomposition of real symmetric or complex Hermitian positive definite matrix a,...
subroutine, public lalg_check_zeigenderivatives(n, mat)
complex(real64) function, public lalg_zdni(eigenvec, alpha, beta)
pure real(real64) function pseudoinverse_default_tolerance(m, n, sg_values)
Computes the default Moore-Penrose pseudoinverse tolerance for zeroing.
subroutine zleast_squares_vec(nn, aa, bb, xx, preserve_mat)
real(real64) function, dimension(1:n, 1:n), public lalg_remove_rotation(n, A)
Remove rotation from affine transformation A by computing the polar decomposition and discarding the ...
subroutine lalg_zgeev(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info)
subroutine, public zlalg_matrix_function(n, factor, a, fun_a, fun, hermitian, tridiagonal)
This routine calculates a function of a matrix by using an eigenvalue decomposition.
integer function dmatrix_rank_svd(a, preserve_mat, tol)
Compute the rank of the matrix A using SVD.
subroutine dinverse(n, a, method, det, threshold, uplo)
An interface to different method to invert a matrix.
subroutine zeigensolve_parallel(n, a, e, bof, err_code)
Computes all the eigenvalues and the eigenvectors of a real symmetric or complex Hermitian eigenprobl...
real(real64) function ddeterminant(n, a, preserve_mat)
Invert a real symmetric or complex Hermitian square matrix a.
subroutine lalg_dgeev(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info)
subroutine dlowest_eigensolve(k, n, a, e, v, preserve_mat)
Computes the k lowest eigenvalues and the eigenvectors of a standard symmetric-definite eigenproblem,...
subroutine dsvd_inverse(m, n, a, threshold)
Computes the inverse of a real M x N matrix, a, using the SVD decomposition.
subroutine, public zlalg_phi(nn, pp, aa, ex, hermitian)
subroutine dlalg_pseudo_inverse(a, threshold)
Invert a matrix with the Moore-Penrose pseudo-inverse.
subroutine zsvd_inverse(m, n, a, threshold)
Computes inverse of a complex MxN matrix, a, using the SVD decomposition.
This module contains interfaces for LAPACK routines.
This module is intended to contain "only mathematical" functions and procedures.
This module contains interfaces for ScaLAPACK routines Interfaces are from http:
This module is intended to contain "only mathematical" functions and procedures.
This module is intended to contain simple general-purpose utility functions and procedures.