26 use,
intrinsic :: iso_fortran_env
147 real(real64) function sfmin()
149 real(real64) function dlamch(cmach)
152 character(1),
intent(in) :: cmach
159 subroutine lalg_dgeev(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info)
160 character(1),
intent(in) :: jobvl, jobvr
161 integer,
intent(in) :: n, lda, ldvl, ldvr, lwork
162 real(real64),
intent(inout) :: a(:,:)
163 complex(real64),
intent(out) :: w(:)
164 real(real64),
intent(out) :: vl(:,:), vr(:,:)
165 real(real64),
intent(out) :: rwork(:)
166 real(real64),
intent(out) :: work(:)
167 integer,
intent(out) :: info
169 real(real64),
allocatable :: wr(:), wi(:)
173 safe_allocate(wr(1:n))
174 safe_allocate(wi(1:n))
176 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)
178 w(1:n) = cmplx(wr(1:n), wi(1:n), real64)
180 safe_deallocate_a(wr)
181 safe_deallocate_a(wi)
186 subroutine lalg_zgeev(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info)
187 character(1),
intent(in) :: jobvl, jobvr
188 integer,
intent(in) :: n, lda, ldvl, ldvr, lwork
189 complex(real64),
intent(inout) :: a(:,:)
190 complex(real64),
intent(out) :: w(:)
191 complex(real64),
intent(out) :: vl(:,:), vr(:,:)
192 real(real64),
intent(out) :: rwork(:)
193 complex(real64),
intent(out) :: work(:)
194 integer,
intent(out) :: info
198 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)
221 integer,
intent(in) :: nn
222 complex(real64),
intent(in) :: pp
223 complex(real64),
intent(in) :: aa(:, :)
224 complex(real64),
intent(inout) :: ex(:, :)
225 logical,
intent(in) :: hermitian
227 complex(real64),
allocatable :: evectors(:, :), zevalues(:)
228 real(real64),
allocatable :: evalues(:)
234 safe_allocate(evectors(1:nn, 1:nn))
237 safe_allocate(evalues(1:nn))
238 safe_allocate(zevalues(1:nn))
240 evectors(1:nn, 1:nn) = aa(1:nn, 1:nn)
244 zevalues(1:nn) =
exp(pp*evalues(1:nn))
247 ex(1:nn, ii) = zevalues(1:nn)*conjg(evectors(ii, 1:nn))
250 ex(1:nn, 1:nn) = matmul(evectors(1:nn, 1:nn), ex(1:nn, 1:nn))
252 safe_deallocate_a(evalues)
253 safe_deallocate_a(zevalues)
255 safe_allocate(zevalues(1:nn))
257 evectors(1:nn, 1:nn) = aa(1:nn, 1:nn)
261 zevalues(1:nn) =
exp(pp*zevalues(1:nn))
263 ex(1:nn, 1:nn) = evectors(1:nn, 1:nn)
268 evectors(1:nn, ii) = zevalues(1:nn)*evectors(1:nn, ii)
271 ex(1:nn, 1:nn) = matmul(ex(1:nn, 1:nn), evectors(1:nn, 1:nn))
273 safe_deallocate_a(zevalues)
276 safe_deallocate_a(evectors)
298 subroutine zlalg_phi(nn, pp, aa, ex, hermitian)
299 integer,
intent(in) :: nn
300 complex(real64),
intent(in) :: pp
301 complex(real64),
intent(in) :: aa(:, :)
302 complex(real64),
intent(inout) :: ex(:, :)
303 logical,
intent(in) :: hermitian
305 complex(real64),
allocatable :: evectors(:, :), zevalues(:)
306 real(real64),
allocatable :: evalues(:)
312 safe_allocate(evectors(1:nn, 1:nn))
315 safe_allocate(evalues(1:nn))
316 safe_allocate(zevalues(1:nn))
323 zevalues(ii) = (
exp(pp*evalues(ii)) - m_z1) / (pp*evalues(ii))
327 ex(1:nn, ii) = zevalues(1:nn)*conjg(evectors(ii, 1:nn))
330 ex(1:nn, 1:nn) = matmul(evectors(1:nn, 1:nn), ex(1:nn, 1:nn))
332 safe_deallocate_a(evalues)
333 safe_deallocate_a(zevalues)
335 safe_allocate(zevalues(1:nn))
337 evectors(1:nn, 1:nn) = aa(1:nn, 1:nn)
342 zevalues(ii) = (
exp(pp*zevalues(ii)) - m_z1) / (pp*zevalues(ii))
345 ex(1:nn, 1:nn) = evectors(1:nn, 1:nn)
350 evectors(1:nn, ii) = zevalues(1:nn)*evectors(1:nn, ii)
353 ex(1:nn, 1:nn) = matmul(ex(1:nn, 1:nn), evectors(1:nn, 1:nn))
355 safe_deallocate_a(zevalues)
374 integer,
intent(in) :: n
375 complex(real64),
intent(in) :: mat(:, :)
376 complex(real64),
intent(out) :: zeigenvec(:, :)
377 complex(real64),
intent(out) :: zeigenval(:)
378 complex(real64),
intent(out) :: zmat(:, :, :)
380 integer :: i, alpha, beta
381 complex(real64),
allocatable :: knaught(:, :)
382 complex(real64),
allocatable :: lambdaminusdm(:, :)
383 complex(real64),
allocatable :: ilambdaminusdm(:, :)
384 complex(real64),
allocatable :: unit(:, :)
388 safe_allocate(unit(1:n, 1:n))
389 safe_allocate(knaught(1:n, 1:n))
390 safe_allocate(lambdaminusdm(1:n, 1:n))
391 safe_allocate(ilambdaminusdm(1:n, 1:n))
404 knaught(alpha, beta) = - zeigenvec(alpha, i)*conjg(zeigenvec(beta, i))
407 knaught = knaught + unit
408 lambdaminusdm = zeigenval(i)*unit - mat
411 zmat(:, :, i) = matmul(ilambdaminusdm, knaught)
414 safe_deallocate_a(unit)
415 safe_deallocate_a(knaught)
416 safe_deallocate_a(lambdaminusdm)
417 safe_deallocate_a(ilambdaminusdm)
426 integer,
intent(in) :: n
427 complex(real64),
intent(in) :: mat(:, :)
428 complex(real64),
intent(out) :: imat(:, :)
431 complex(real64),
allocatable :: u(:, :), vt(:, :), sigma(:, :)
432 real(real64),
allocatable :: sg_values(:)
436 safe_allocate(u(1:n, 1:n))
437 safe_allocate(vt(1:n, 1:n))
438 safe_allocate(sigma(1:n, 1:n))
439 safe_allocate(sg_values(1:n))
446 if (abs(sg_values(i)) <= 1.0e-12_real64 * maxval(abs(sg_values)))
then
449 sigma(i, i) = m_z1 / sg_values(i)
453 vt = conjg(transpose(vt))
454 u = conjg(transpose(u))
455 imat = matmul(vt, matmul(sigma, u))
458 vt = matmul(mat, matmul(imat, mat)) - mat
459 if (maxval(abs(vt)) > 1.0e-10_real64 * maxval(abs(mat)))
then
460 write(*, *) maxval(abs(vt))
463 write(*, *) 1.0e-10_real64 * maxval(abs(mat))
464 write(*, *) maxval(abs(vt)) > 1.0e-10_real64 * maxval(abs(mat))
466 write(message(1),
'(a)')
'Pseudoinverse failed.'
467 call messages_fatal(1)
471 safe_deallocate_a(vt)
472 safe_deallocate_a(sigma)
473 safe_deallocate_a(sg_values)
485 integer,
intent(in) :: n
486 complex(real64),
intent(in) :: mat(:, :)
488 integer :: alpha, beta, gamma, delta
489 complex(real64) :: deltah, zuder_direct, zder_direct, zuder_directplus, zuder_directminus
490 complex(real64),
allocatable :: zeigenvec(:, :), dm(:, :), zeigref_(:, :), zeigenval(:), mmatrix(:, :, :)
491 complex(real64),
allocatable :: zeigplus(:), zeigminus(:), zeig0(:), zeigplusminus(:), zeigminusplus(:)
495 safe_allocate(zeigenvec(1:n, 1:n))
496 safe_allocate(dm(1:n, 1:n))
497 safe_allocate(zeigref_(1:n, 1:n))
498 safe_allocate(zeigenval(1:n))
499 safe_allocate(mmatrix(1:n, 1:n, 1:n))
500 safe_allocate(zeigplus(1:n))
501 safe_allocate(zeigminus(1:n))
502 safe_allocate(zeig0(1:n))
503 safe_allocate(zeigplusminus(1:n))
504 safe_allocate(zeigminusplus(1:n))
512 deltah = (0.000001_real64, 0.000_real64)
516 zder_direct =
lalg_zdni(zeigref_(:, 2), alpha, beta)
517 zuder_direct =
lalg_zduialpha(zeigref_(:, 2), mmatrix(:, :, 2), 2, alpha, beta)
520 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
522 zeigenvec(:, 1) = zeigenvec(:, 1)/sum(conjg(zeigref_(1:2, 1))*zeigenvec(1:2, 1))
523 zeigenvec(:, 2) = zeigenvec(:, 2)/sum(conjg(zeigref_(1:2, 2))*zeigenvec(1:2, 2))
524 zuder_directplus = zeigenvec(2, 2)
527 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
529 zeigenvec(:, 1) = zeigenvec(:, 1)/sum(conjg(zeigref_(1:2, 1))*zeigenvec(1:2, 1))
530 zeigenvec(:, 2) = zeigenvec(:, 2)/sum(conjg(zeigref_(1:2, 2))*zeigenvec(1:2, 2))
531 zuder_directminus = zeigenvec(2, 2)
533 write(*,
'(2i1,4f24.12)') alpha, beta, zder_direct, (zeigplus(2) - zeigminus(2))/(m_two * deltah)
534 write(*,
'(2i1,4f24.12)') alpha, beta, &
535 zuder_direct, (zuder_directplus - zuder_directminus) / (m_two * deltah)
539 if (alpha == gamma .and. beta == delta)
then
540 zder_direct =
lalg_zd2ni(zeigref_(:, 1), mmatrix(:, :, 1), alpha, beta, gamma, delta)
546 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
550 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
553 write(*,
'(4i1,4f24.12)') alpha, beta, gamma, delta, &
554 zder_direct, (zeigplus(1) + zeigminus(1) - m_two*zeig0(1))/(deltah**2)
556 zder_direct =
lalg_zd2ni(zeigref_(:, 1), mmatrix(:, :, 1), alpha, beta, gamma, delta)
559 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
560 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) + deltah
564 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
565 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) - deltah
569 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
570 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) - deltah
574 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
575 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) + deltah
578 write(*,
'(4i1,4f24.12)') alpha, beta, gamma, delta, &
579 zder_direct, (zeigplus(1) + zeigminus(1) - zeigplusminus(1) - zeigminusplus(1))/(m_four*deltah**2)
589 safe_deallocate_a(zeigenval)
590 safe_deallocate_a(mmatrix)
591 safe_deallocate_a(zeigenvec)
592 safe_deallocate_a(dm)
593 safe_deallocate_a(zeigref_)
594 safe_deallocate_a(zeigplus)
595 safe_deallocate_a(zeigminus)
596 safe_deallocate_a(zeig0)
597 safe_deallocate_a(zeigplusminus)
598 safe_deallocate_a(zeigminusplus)
602 complex(real64) function lalg_zdni(eigenvec, alpha, beta)
603 integer,
intent(in) :: alpha, beta
604 complex(real64) :: eigenvec(2)
605 lalg_zdni = conjg(eigenvec(alpha)) * eigenvec(beta)
608 complex(real64) function lalg_zduialpha(eigenvec, mmatrix, alpha, gamma, delta)
609 integer,
intent(in) :: alpha, gamma, delta
610 complex(real64) :: eigenvec(2), mmatrix(2, 2)
614 complex(real64) function lalg_zd2ni(eigenvec, mmatrix, alpha, beta, gamma, delta)
615 integer,
intent(in) :: alpha, beta, gamma, delta
616 complex(real64) :: eigenvec(2), mmatrix(2, 2)
617 lalg_zd2ni = conjg(mmatrix(alpha, delta) * eigenvec(gamma)) * eigenvec(beta) + &
618 conjg(eigenvec(alpha)) * mmatrix(beta, gamma) * eigenvec(delta)
627 integer,
intent(in) :: m, n
628 real(real64),
intent(in) :: sg_values(:)
629 tol = m_epsilon * m * n * maxval(sg_values)
639 integer,
intent(in) :: n
640 real(real64),
intent(in) :: a(1:n, 1:n)
641 real(real64) :: p(1:n, 1:n)
643 real(real64) :: ata(1:n, 1:n)
647 ata = matmul(transpose(a), a)
654#include "complex.F90"
655#include "lalg_adv_lapack_inc.F90"
659#include "lalg_adv_lapack_inc.F90"
Note that lalg_determinant and lalg_inverse are just wrappers over the same routine.
lalg_gemm with both the (Hermitian) transpose of A and B.
scales a vector by a constant
double exp(double __x) __attribute__((__nothrow__
This module contains interfaces for BLACS routines Interfaces are from http:
This module provides the BLACS processor grid.
subroutine zmatrix_norm2(m, n, a, norm_l2, preserve_mat)
Norm of a 2D matrix.
complex(real64) function, public lalg_zd2ni(eigenvec, mmatrix, alpha, beta, gamma, delta)
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 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 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 dmatrix_norm2(m, n, a, norm_l2, preserve_mat)
Norm of a 2D matrix.
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 dlalg_matrix_function(n, factor, a, fun_a, fun, hermitian)
This routine calculates a function of a matrix by using an eigenvalue decomposition.
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, public zlalg_matrix_function(n, factor, a, fun_a, fun, hermitian)
This routine calculates a function of a matrix by using an eigenvalue decomposition.
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.