27 use,
intrinsic :: iso_fortran_env
153 real(real64) function sfmin()
155 real(real64) function dlamch(cmach)
158 character(1),
intent(in) :: cmach
165 subroutine lalg_dgeev(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info)
166 character(1),
intent(in) :: jobvl, jobvr
167 integer,
intent(in) :: n, lda, ldvl, ldvr, lwork
168 real(real64),
intent(inout) :: a(:, :)
169 complex(real64),
intent(out) :: w(:)
170 real(real64),
intent(out) :: vl(:, :), vr(:, :)
171 real(real64),
intent(out) :: rwork(:)
172 real(real64),
intent(out) :: work(:)
173 integer,
intent(out) :: info
175 real(real64),
allocatable :: wr(:), wi(:)
179 safe_allocate(wr(1:n))
180 safe_allocate(wi(1:n))
182 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)
184 w(1:n) = cmplx(wr(1:n), wi(1:n), real64)
186 safe_deallocate_a(wr)
187 safe_deallocate_a(wi)
192 subroutine lalg_zgeev(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info)
193 character(1),
intent(in) :: jobvl, jobvr
194 integer,
intent(in) :: n, lda, ldvl, ldvr, lwork
195 complex(real64),
intent(inout) :: a(:, :)
196 complex(real64),
intent(out) :: w(:)
197 complex(real64),
intent(out) :: vl(:, :), vr(:, :)
198 real(real64),
intent(out) :: rwork(:)
199 complex(real64),
intent(out) :: work(:)
200 integer,
intent(out) :: info
204 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)
227 integer,
intent(in) :: nn
228 complex(real64),
intent(in) :: pp
229 complex(real64),
intent(in) :: aa(:, :)
230 complex(real64),
intent(inout) :: ex(:, :)
231 logical,
intent(in) :: hermitian
233 complex(real64),
allocatable :: evectors(:, :), zevalues(:)
234 real(real64),
allocatable :: evalues(:)
240 safe_allocate(evectors(1:nn, 1:nn))
243 safe_allocate(evalues(1:nn))
244 safe_allocate(zevalues(1:nn))
246 evectors(1:nn, 1:nn) = aa(1:nn, 1:nn)
250 zevalues(1:nn) =
exp(pp*evalues(1:nn))
253 ex(1:nn, ii) = zevalues(1:nn)*conjg(evectors(ii, 1:nn))
256 ex(:, :) = matmul(evectors(:, :), ex(:, :))
258 safe_deallocate_a(evalues)
259 safe_deallocate_a(zevalues)
261 safe_allocate(zevalues(1:nn))
263 evectors(1:nn, 1:nn) = aa(1:nn, 1:nn)
267 zevalues(1:nn) =
exp(pp*zevalues(1:nn))
269 ex(1:nn, 1:nn) = evectors(1:nn, 1:nn)
274 evectors(1:nn, ii) = zevalues(1:nn)*evectors(1:nn, ii)
277 ex(:, :) = matmul(ex(:, :), evectors(:, :))
279 safe_deallocate_a(zevalues)
282 safe_deallocate_a(evectors)
304 subroutine zlalg_phi(nn, pp, aa, ex, hermitian)
305 integer,
intent(in) :: nn
306 complex(real64),
intent(in) :: pp
307 complex(real64),
intent(in) :: aa(:, :)
308 complex(real64),
intent(inout) :: ex(:, :)
309 logical,
intent(in) :: hermitian
311 complex(real64),
allocatable :: evectors(:, :), zevalues(:)
312 real(real64),
allocatable :: evalues(:)
318 safe_allocate(evectors(1:nn, 1:nn))
321 safe_allocate(evalues(1:nn))
322 safe_allocate(zevalues(1:nn))
324 evectors(:, :) = aa(:, :)
329 zevalues(ii) = (
exp(pp*evalues(ii)) - m_z1) / (pp*evalues(ii))
333 ex(1:nn, ii) = zevalues(1:nn)*conjg(evectors(ii, 1:nn))
336 ex(:, :) = matmul(evectors(:, :), ex(:, :))
338 safe_deallocate_a(evalues)
339 safe_deallocate_a(zevalues)
341 safe_allocate(zevalues(1:nn))
343 evectors(:, :) = aa(:, :)
348 zevalues(ii) = (
exp(pp*zevalues(ii)) - m_z1) / (pp*zevalues(ii))
351 ex(:, :) = evectors(:, :)
356 evectors(1:nn, ii) = zevalues(1:nn)*evectors(1:nn, ii)
359 ex(:, :) = matmul(ex(:, :), evectors(:, :))
361 safe_deallocate_a(zevalues)
380 integer,
intent(in) :: n
381 complex(real64),
contiguous,
intent(in) :: mat(:, :)
382 complex(real64),
contiguous,
intent(out) :: zeigenvec(:, :)
383 complex(real64),
contiguous,
intent(out) :: zeigenval(:)
384 complex(real64),
contiguous,
intent(out) :: zmat(:, :, :)
386 integer :: i, alpha, beta
387 complex(real64),
allocatable :: knaught(:, :)
388 complex(real64),
allocatable :: lambdaminusdm(:, :)
389 complex(real64),
allocatable :: ilambdaminusdm(:, :)
390 complex(real64),
allocatable :: unit(:, :)
394 safe_allocate(unit(1:n, 1:n))
395 safe_allocate(knaught(1:n, 1:n))
396 safe_allocate(lambdaminusdm(1:n, 1:n))
397 safe_allocate(ilambdaminusdm(1:n, 1:n))
410 knaught(alpha, beta) = - zeigenvec(alpha, i)*conjg(zeigenvec(beta, i))
413 knaught = knaught + unit
414 lambdaminusdm = zeigenval(i)*unit - mat
417 zmat(:, :, i) = matmul(ilambdaminusdm, knaught)
420 safe_deallocate_a(unit)
421 safe_deallocate_a(knaught)
422 safe_deallocate_a(lambdaminusdm)
423 safe_deallocate_a(ilambdaminusdm)
432 integer,
intent(in) :: n
433 complex(real64),
contiguous,
intent(in) :: mat(:, :)
434 complex(real64),
contiguous,
intent(out) :: imat(:, :)
437 complex(real64),
allocatable :: u(:, :), vt(:, :), sigma(:, :)
438 real(real64),
allocatable :: sg_values(:)
442 safe_allocate(u(1:n, 1:n))
443 safe_allocate(vt(1:n, 1:n))
444 safe_allocate(sigma(1:n, 1:n))
445 safe_allocate(sg_values(1:n))
452 if (abs(sg_values(i)) <= 1.0e-12_real64 * maxval(abs(sg_values)))
then
455 sigma(i, i) = m_z1 / sg_values(i)
459 vt(:, :) = conjg(transpose(vt(:, :)))
460 u(:, :) = conjg(transpose(u(:, :)))
461 imat = matmul(vt, matmul(sigma, u))
464 vt(:, :) = matmul(mat(:, :), matmul(imat(:, :), mat(:, :))) - mat(:, :)
465 if (maxval(abs(vt)) > 1.0e-10_real64 * maxval(abs(mat)))
then
466 write(*, *) maxval(abs(vt))
469 write(*, *) 1.0e-10_real64 * maxval(abs(mat))
470 write(*, *) maxval(abs(vt)) > 1.0e-10_real64 * maxval(abs(mat))
472 write(message(1),
'(a)')
'Pseudoinverse failed.'
473 call messages_fatal(1)
477 safe_deallocate_a(vt)
478 safe_deallocate_a(sigma)
479 safe_deallocate_a(sg_values)
491 integer,
intent(in) :: n
492 complex(real64),
intent(in) :: mat(:, :)
494 integer :: alpha, beta, gamma, delta
495 complex(real64) :: deltah, zuder_direct, zder_direct, zuder_directplus, zuder_directminus
496 complex(real64),
allocatable :: zeigenvec(:, :), dm(:, :), zeigref_(:, :), zeigenval(:), mmatrix(:, :, :)
497 complex(real64),
allocatable :: zeigplus(:), zeigminus(:), zeig0(:), zeigplusminus(:), zeigminusplus(:)
501 safe_allocate(zeigenvec(1:n, 1:n))
502 safe_allocate(dm(1:n, 1:n))
503 safe_allocate(zeigref_(1:n, 1:n))
504 safe_allocate(zeigenval(1:n))
505 safe_allocate(mmatrix(1:n, 1:n, 1:n))
506 safe_allocate(zeigplus(1:n))
507 safe_allocate(zeigminus(1:n))
508 safe_allocate(zeig0(1:n))
509 safe_allocate(zeigplusminus(1:n))
510 safe_allocate(zeigminusplus(1:n))
518 deltah = (0.000001_real64, 0.000_real64)
522 zder_direct =
lalg_zdni(zeigref_(:, 2), alpha, beta)
523 zuder_direct =
lalg_zduialpha(zeigref_(:, 2), mmatrix(:, :, 2), 2, alpha, beta)
525 zeigenvec(:, :) = dm(:, :)
526 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
528 zeigenvec(:, 1) = zeigenvec(:, 1)/sum(conjg(zeigref_(1:2, 1))*zeigenvec(1:2, 1))
529 zeigenvec(:, 2) = zeigenvec(:, 2)/sum(conjg(zeigref_(1:2, 2))*zeigenvec(1:2, 2))
530 zuder_directplus = zeigenvec(2, 2)
532 zeigenvec(:, :) = dm(:, :)
533 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
535 zeigenvec(:, 1) = zeigenvec(:, 1)/sum(conjg(zeigref_(1:2, 1))*zeigenvec(1:2, 1))
536 zeigenvec(:, 2) = zeigenvec(:, 2)/sum(conjg(zeigref_(1:2, 2))*zeigenvec(1:2, 2))
537 zuder_directminus = zeigenvec(2, 2)
539 write(*,
'(2i1,4f24.12)') alpha, beta, zder_direct, (zeigplus(2) - zeigminus(2))/(m_two * deltah)
540 write(*,
'(2i1,4f24.12)') alpha, beta, &
541 zuder_direct, (zuder_directplus - zuder_directminus) / (m_two * deltah)
545 if (alpha == gamma .and. beta == delta)
then
546 zder_direct =
lalg_zd2ni(zeigref_(:, 1), mmatrix(:, :, 1), alpha, beta, gamma, delta)
548 zeigenvec(:, :) = dm(:, :)
551 zeigenvec(:, :) = dm(:, :)
552 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
555 zeigenvec(:, :) = dm(:, :)
556 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
559 write(*,
'(4i1,4f24.12)') alpha, beta, gamma, delta, &
560 zder_direct, (zeigplus(1) + zeigminus(1) - m_two*zeig0(1))/(deltah**2)
562 zder_direct =
lalg_zd2ni(zeigref_(:, 1), mmatrix(:, :, 1), alpha, beta, gamma, delta)
564 zeigenvec(:, :) = dm(:, :)
565 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
566 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) + deltah
569 zeigenvec(:, :) = dm(:, :)
570 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
571 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) - deltah
574 zeigenvec(:, :) = dm(:, :)
575 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
576 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) - deltah
579 zeigenvec(:, :) = dm(:, :)
580 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
581 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) + deltah
584 write(*,
'(4i1,4f24.12)') alpha, beta, gamma, delta, &
585 zder_direct, (zeigplus(1) + zeigminus(1) - zeigplusminus(1) - zeigminusplus(1))/(m_four*deltah**2)
595 safe_deallocate_a(zeigenval)
596 safe_deallocate_a(mmatrix)
597 safe_deallocate_a(zeigenvec)
598 safe_deallocate_a(dm)
599 safe_deallocate_a(zeigref_)
600 safe_deallocate_a(zeigplus)
601 safe_deallocate_a(zeigminus)
602 safe_deallocate_a(zeig0)
603 safe_deallocate_a(zeigplusminus)
604 safe_deallocate_a(zeigminusplus)
608 complex(real64) function lalg_zdni(eigenvec, alpha, beta)
609 integer,
intent(in) :: alpha, beta
610 complex(real64) :: eigenvec(2)
611 lalg_zdni = conjg(eigenvec(alpha)) * eigenvec(beta)
614 complex(real64) function lalg_zduialpha(eigenvec, mmatrix, alpha, gamma, delta)
615 integer,
intent(in) :: alpha, gamma, delta
616 complex(real64) :: eigenvec(2), mmatrix(2, 2)
620 complex(real64) function lalg_zd2ni(eigenvec, mmatrix, alpha, beta, gamma, delta)
621 integer,
intent(in) :: alpha, beta, gamma, delta
622 complex(real64) :: eigenvec(2), mmatrix(2, 2)
623 lalg_zd2ni = conjg(mmatrix(alpha, delta) * eigenvec(gamma)) * eigenvec(beta) + &
624 conjg(eigenvec(alpha)) * mmatrix(beta, gamma) * eigenvec(delta)
633 integer,
intent(in) :: m, n
634 real(real64),
intent(in) :: sg_values(:)
635 tol = m_epsilon * m * n * maxval(sg_values)
645 integer,
intent(in) :: n
646 real(real64),
intent(in) :: a(1:n, 1:n)
647 real(real64) :: p(1:n, 1:n)
649 real(real64) :: ata(1:n, 1:n)
653 ata = matmul(transpose(a), a)
660#include "complex.F90"
661#include "lalg_adv_lapack_inc.F90"
665#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....
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 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 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, 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.
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.