Octopus
lalg_adv.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
21module lalg_adv_oct_m
22 use blas_oct_m
24 use blacs_oct_m
25 use debug_oct_m
26 use global_oct_m
27 use, intrinsic :: iso_fortran_env
28 use lapack_oct_m
30 use math_oct_m
32 use mpi_oct_m
35 use sort_oct_m
36 use utils_oct_m
37#ifdef HAVE_ELPA
38 use elpa
39#endif
40
41 implicit none
42
43 private
44 public :: &
59 zlalg_exp, &
60 zlalg_phi, &
65 lalg_zdni, &
67 lalg_zd2ni, &
71
72 interface lalg_cholesky
73 module procedure dcholesky, zcholesky
74 end interface lalg_cholesky
75
76 interface lalg_geneigensolve
77 module procedure dgeneigensolve, zgeneigensolve
78 end interface lalg_geneigensolve
79
80 interface lalg_eigensolve_nonh
81 module procedure zeigensolve_nonh, deigensolve_nonh
82 end interface lalg_eigensolve_nonh
83
84 interface lalg_eigensolve
85 module procedure deigensolve, zeigensolve
86 end interface lalg_eigensolve
87
90 end interface lalg_eigensolve_tridiagonal
91
94 end interface lalg_eigensolve_parallel
95
98
99 interface lalg_determinant
100 module procedure ddeterminant, zdeterminant
101 end interface lalg_determinant
102
103 interface lalg_inverse
104 module procedure dinverse, zinverse
105 end interface lalg_inverse
106
107 interface lalg_linsyssolve
108 module procedure dlinsyssolve, zlinsyssolve
109 end interface lalg_linsyssolve
110
113 end interface lalg_singular_value_decomp
114
115 interface lalg_pseudo_inverse
117 end interface lalg_pseudo_inverse
118
119 interface lalg_svd_inverse
120 module procedure dsvd_inverse, zsvd_inverse
121 end interface lalg_svd_inverse
122
125 end interface lalg_lowest_geneigensolve
126
127 interface lalg_lowest_eigensolve
128 module procedure dlowest_eigensolve, zlowest_eigensolve
129 end interface lalg_lowest_eigensolve
130
131 interface lapack_geev
132 module procedure lalg_dgeev, lalg_zgeev
133 end interface lapack_geev
134
135 interface lalg_least_squares
136 module procedure dleast_squares_vec, zleast_squares_vec
137 end interface lalg_least_squares
138
139 interface lalg_matrix_function
141 end interface lalg_matrix_function
142
143 interface lalg_matrix_rank_svd
144 module procedure dmatrix_rank_svd, zmatrix_rank_svd
145 end interface lalg_matrix_rank_svd
146
147contains
148
149 ! ---------------------------------------------------------
151 real(real64) function sfmin()
152 interface
153 real(real64) function dlamch(cmach)
154 import real64
155 implicit none
156 character(1), intent(in) :: cmach
157 end function dlamch
158 end interface
159
160 sfmin = dlamch('S')
161 end function sfmin
162
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
172
173 real(real64), allocatable :: wr(:), wi(:)
174
175 push_sub(lalg_dgeev)
176
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)
181
182 w(1:n) = cmplx(wr(1:n), wi(1:n), real64)
184 safe_deallocate_a(wr)
185 safe_deallocate_a(wi)
186
187 pop_sub(lalg_dgeev)
188 end subroutine lalg_dgeev
189
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
199
200 push_sub(lalg_zgeev)
201
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)
203
204 pop_sub(lalg_zgeev)
205 end subroutine lalg_zgeev
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(:)
233
234 integer :: ii
235
236 push_sub(zlalg_exp)
237
238 safe_allocate(evectors(1:nn, 1:nn))
239
240 if (hermitian) then
241 safe_allocate(evalues(1:nn))
242 safe_allocate(zevalues(1:nn))
243
244 evectors(1:nn, 1:nn) = aa(1:nn, 1:nn)
245
246 call lalg_eigensolve(nn, evectors, evalues)
247
248 zevalues(1:nn) = exp(pp*evalues(1:nn))
249
250 do ii = 1, nn
251 ex(1:nn, ii) = zevalues(1:nn)*conjg(evectors(ii, 1:nn))
252 end do
253
254 ex(:, :) = matmul(evectors(:, :), ex(:, :))
255
256 safe_deallocate_a(evalues)
257 safe_deallocate_a(zevalues)
258 else
259 safe_allocate(zevalues(1:nn))
260
261 evectors(1:nn, 1:nn) = aa(1:nn, 1:nn)
262
263 call lalg_eigensolve_nonh(nn, evectors, zevalues)
264
265 zevalues(1:nn) = exp(pp*zevalues(1:nn))
266
267 ex(1:nn, 1:nn) = evectors(1:nn, 1:nn)
268
269 call lalg_inverse(nn, evectors, 'dir')
270
271 do ii = 1, nn
272 evectors(1:nn, ii) = zevalues(1:nn)*evectors(1:nn, ii)
273 end do
274
275 ex(:, :) = matmul(ex(:, :), evectors(:, :))
276
277 safe_deallocate_a(zevalues)
278 end if
279
280 safe_deallocate_a(evectors)
281
282 pop_sub(zlalg_exp)
283 end subroutine zlalg_exp
284
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
308
309 complex(real64), allocatable :: evectors(:, :), zevalues(:)
310 real(real64), allocatable :: evalues(:)
311
312 integer :: ii
313
314 push_sub(zlalg_phi)
315
316 safe_allocate(evectors(1:nn, 1:nn))
317
318 if (hermitian) then
319 safe_allocate(evalues(1:nn))
320 safe_allocate(zevalues(1:nn))
321
322 evectors(:, :) = aa(:, :)
323
324 call lalg_eigensolve(nn, evectors, evalues)
325
326 do ii = 1, nn
327 zevalues(ii) = (exp(pp*evalues(ii)) - m_z1) / (pp*evalues(ii))
328 end do
329
330 do ii = 1, nn
331 ex(1:nn, ii) = zevalues(1:nn)*conjg(evectors(ii, 1:nn))
332 end do
333
334 ex(:, :) = matmul(evectors(:, :), ex(:, :))
335
336 safe_deallocate_a(evalues)
337 safe_deallocate_a(zevalues)
338 else
339 safe_allocate(zevalues(1:nn))
340
341 evectors(:, :) = aa(:, :)
342
343 call lalg_eigensolve_nonh(nn, evectors, zevalues)
344
345 do ii = 1, nn
346 zevalues(ii) = (exp(pp*zevalues(ii)) - m_z1) / (pp*zevalues(ii))
347 end do
348
349 ex(:, :) = evectors(:, :)
350
351 call lalg_inverse(nn, evectors, 'dir')
352
353 do ii = 1, nn
354 evectors(1:nn, ii) = zevalues(1:nn)*evectors(1:nn, ii)
355 end do
356
357 ex(:, :) = matmul(ex(:, :), evectors(:, :))
358
359 safe_deallocate_a(zevalues)
360 end if
361
362 pop_sub(zlalg_phi)
363 end subroutine zlalg_phi
364
365
377 subroutine lalg_zeigenderivatives(n, mat, zeigenvec, zeigenval, zmat)
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(:, :, :)
383
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(:, :)
389
390 push_sub(lalg_zeigenderivatives)
391
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))
396
397 zeigenvec = mat
398 call lalg_eigensolve_nonh(n, zeigenvec, zeigenval)
399
400 unit = m_z0
401 do i = 1, n
402 unit(i, i) = m_z1
403 end do
404
405 do i = 1, n
406 do alpha = 1, n
407 do beta = 1, n
408 knaught(alpha, beta) = - zeigenvec(alpha, i)*conjg(zeigenvec(beta, i))
409 end do
410 end do
411 knaught = knaught + unit
412 lambdaminusdm = zeigenval(i)*unit - mat
413 ! TODO(Alex) Swap this call
414 call lalg_zpseudoinverse(n, lambdaminusdm, ilambdaminusdm)
415 zmat(:, :, i) = matmul(ilambdaminusdm, knaught)
416 end do
417
418 safe_deallocate_a(unit)
419 safe_deallocate_a(knaught)
420 safe_deallocate_a(lambdaminusdm)
421 safe_deallocate_a(ilambdaminusdm)
423 end subroutine lalg_zeigenderivatives
424
425
426 ! TODO(Alex) Issue #1132.
427 ! * This routine is redunant. Its usage should be replaced with either zsvd_inverse or zlalg_pseudo_inverse.
429 subroutine lalg_zpseudoinverse(n, mat, imat)
430 integer, intent(in) :: n
431 complex(real64), contiguous, intent(in) :: mat(:, :)
432 complex(real64), contiguous, intent(out) :: imat(:, :)
433
434 integer :: i
435 complex(real64), allocatable :: u(:, :), vt(:, :), sigma(:, :)
436 real(real64), allocatable :: sg_values(:)
437
438 push_sub(lalg_zpseudoinverse)
439
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))
444
445 imat = mat
446 call lalg_singular_value_decomp(n, n, imat, u, vt, sg_values)
447
448 sigma = m_z0
449 do i = 1, n
450 if (abs(sg_values(i)) <= 1.0e-12_real64 * maxval(abs(sg_values))) then
451 sigma(i, i) = m_z0
452 else
453 sigma(i, i) = m_z1 / sg_values(i)
454 end if
455 end do
456
457 vt(:, :) = conjg(transpose(vt(:, :)))
458 u(:, :) = conjg(transpose(u(:, :)))
459 imat = matmul(vt, matmul(sigma, u))
460
461 ! Check if we truly have a pseudoinverse
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))
465 write(*, *) vt
466 write(*, *)
467 write(*, *) 1.0e-10_real64 * maxval(abs(mat))
468 write(*, *) maxval(abs(vt)) > 1.0e-10_real64 * maxval(abs(mat))
469 write(*, *) mat
470 write(message(1), '(a)') 'Pseudoinverse failed.'
471 call messages_fatal(1)
472 end if
473
474 safe_deallocate_a(u)
475 safe_deallocate_a(vt)
476 safe_deallocate_a(sigma)
477 safe_deallocate_a(sg_values)
478 pop_sub(lalg_zpseudoinverse)
479 end subroutine lalg_zpseudoinverse
480
481
488 subroutine lalg_check_zeigenderivatives(n, mat)
489 integer, intent(in) :: n
490 complex(real64), intent(in) :: mat(:, :)
491
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(:)
496
498
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))
509
510 assert(n == 2)
511
512 dm(:, :) = mat(:, :)
513 call lalg_zeigenderivatives(2, dm, zeigref_, zeigenval, mmatrix)
514
515
516 deltah = (0.000001_real64, 0.000_real64)
517 !deltah = M_z1 * maxval(abs(dm)) * 0.001_real64
518 do alpha = 1, n
519 do beta = 1, n
520 zder_direct = lalg_zdni(zeigref_(:, 2), alpha, beta)
521 zuder_direct = lalg_zduialpha(zeigref_(:, 2), mmatrix(:, :, 2), 2, alpha, beta)
522
523 zeigenvec(:, :) = dm(:, :)
524 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
525 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigplus)
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)
529
530 zeigenvec(:, :) = dm(:, :)
531 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
532 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigminus)
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)
536
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)
540
541 do gamma = 1, n
542 do delta = 1, n
543 if (alpha == gamma .and. beta == delta) then
544 zder_direct = lalg_zd2ni(zeigref_(:, 1), mmatrix(:, :, 1), alpha, beta, gamma, delta)
545
546 zeigenvec(:, :) = dm(:, :)
547 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeig0)
548
549 zeigenvec(:, :) = dm(:, :)
550 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
551 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigplus)
552
553 zeigenvec(:, :) = dm(:, :)
554 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
555 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigminus)
556
557 write(*, '(4i1,4f24.12)') alpha, beta, gamma, delta, &
558 zder_direct, (zeigplus(1) + zeigminus(1) - m_two*zeig0(1))/(deltah**2)
559 else
560 zder_direct = lalg_zd2ni(zeigref_(:, 1), mmatrix(:, :, 1), alpha, beta, gamma, delta)
561
562 zeigenvec(:, :) = dm(:, :)
563 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
564 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) + deltah
565 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigplus)
566
567 zeigenvec(:, :) = dm(:, :)
568 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
569 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) - deltah
570 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigminus)
571
572 zeigenvec(:, :) = dm(:, :)
573 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
574 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) - deltah
575 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigplusminus)
576
577 zeigenvec(:, :) = dm(:, :)
578 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
579 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) + deltah
580 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigminusplus)
581
582 write(*, '(4i1,4f24.12)') alpha, beta, gamma, delta, &
583 zder_direct, (zeigplus(1) + zeigminus(1) - zeigplusminus(1) - zeigminusplus(1))/(m_four*deltah**2)
584
585
586 end if
587 end do
588 end do
589
590 end do
591 end do
592
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)
604 end subroutine lalg_check_zeigenderivatives
605
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)
610 end function lalg_zdni
611
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)
615 lalg_zduialpha = mmatrix(alpha, gamma) * eigenvec(delta)
616 end function lalg_zduialpha
617
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)
623 end function lalg_zd2ni
624
625 ! ---------------------------------------------------------
630 pure real(real64) function pseudoinverse_default_tolerance(m, n, sg_values) result(tol)
631 integer, intent(in) :: m, n
632 real(real64), intent(in) :: sg_values(:)
633 tol = m_epsilon * m * n * maxval(sg_values)
635
636
642 function lalg_remove_rotation(n, A) result(P)
643 integer, intent(in) :: n
644 real(real64), intent(in) :: a(1:n, 1:n)
645 real(real64) :: p(1:n, 1:n)
646
647 real(real64) :: ata(1:n, 1:n)
648
649 push_sub(lalg_remove_rotation)
650
651 ata = matmul(transpose(a), a)
652 call lalg_matrix_function(n, m_one, ata, p, square_root, .true.)
653
654 pop_sub(lalg_remove_rotation)
655 end function lalg_remove_rotation
656
657#include "undef.F90"
658#include "complex.F90"
659#include "lalg_adv_lapack_inc.F90"
660
661#include "undef.F90"
662#include "real.F90"
663#include "lalg_adv_lapack_inc.F90"
664
665end module lalg_adv_oct_m
666
667!! Local Variables:
668!! mode: f90
669!! coding: utf-8
670!! End:
Note that lalg_determinant and lalg_inverse are just wrappers over the same routine.
Definition: lalg_adv.F90:194
double exp(double __x) __attribute__((__nothrow__
This module contains interfaces for BLACS routines Interfaces are from http:
Definition: blacs.F90:27
This module provides the BLACS processor grid.
This module contains interfaces for BLAS routines You should not use these routines directly....
Definition: blas.F90:120
complex(real64) function, public lalg_zd2ni(eigenvec, mmatrix, alpha, beta, gamma, delta)
Definition: lalg_adv.F90:714
integer function zmatrix_rank_svd(a, preserve_mat, tol)
Compute the rank of the matrix A using SVD.
Definition: lalg_adv.F90:1917
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...
Definition: lalg_adv.F90:1108
complex(real64) function zdeterminant(n, a, preserve_mat)
Invert a real symmetric or complex Hermitian square matrix a.
Definition: lalg_adv.F90:1467
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...
Definition: lalg_adv.F90:989
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...
Definition: lalg_adv.F90:1624
subroutine, public zlalg_exp(nn, pp, aa, ex, hermitian)
Definition: lalg_adv.F90:320
subroutine deigensolve_parallel(n, a, e, bof, err_code)
Computes all the eigenvalues and the eigenvectors of a real symmetric or complex Hermitian eigenprobl...
Definition: lalg_adv.F90:3579
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...
Definition: lalg_adv.F90:2389
subroutine lalg_zpseudoinverse(n, mat, imat)
Computes the Moore-Penrose pseudoinverse of a complex matrix.
Definition: lalg_adv.F90:525
subroutine zlalg_pseudo_inverse(a, threshold)
Invert a matrix with the Moore-Penrose pseudo-inverse.
Definition: lalg_adv.F90:1851
subroutine zeigensolve(n, a, e, bof, err_code)
Computes all eigenvalues and eigenvectors of a real symmetric or hermitian square matrix A.
Definition: lalg_adv.F90:1227
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.
Definition: lalg_adv.F90:3207
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...
Definition: lalg_adv.F90:2495
subroutine deigensolve_tridiagonal(n, a, e, bof, err_code)
Computes all eigenvalues and eigenvectors of a real symmetric tridiagonal matrix. For the Hermitian c...
Definition: lalg_adv.F90:2804
subroutine zinverse(n, a, method, det, threshold, uplo)
An interface to different method to invert a matrix.
Definition: lalg_adv.F90:2104
real(real64) function sfmin()
Auxiliary function.
Definition: lalg_adv.F90:247
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 ...
Definition: lalg_adv.F90:3129
subroutine zsingular_value_decomp(m, n, a, u, vt, sg_values, preserve_mat)
Computes the singular value decomposition of a complex MxN matrix a.
Definition: lalg_adv.F90:1702
subroutine dcholesky(n, a, bof, err_code)
Compute the Cholesky decomposition of real symmetric or complex Hermitian positive definite matrix a,...
Definition: lalg_adv.F90:2334
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.
Definition: lalg_adv.F90:3650
subroutine deigensolve(n, a, e, bof, err_code)
Computes all eigenvalues and eigenvectors of a real symmetric or hermitian square matrix A.
Definition: lalg_adv.F90:2734
subroutine, public lalg_zeigenderivatives(n, mat, zeigenvec, zeigenval, zmat)
Definition: lalg_adv.F90:473
subroutine zeigensolve_tridiagonal(n, a, e, bof, err_code)
Computes all eigenvalues and eigenvectors of a real symmetric tridiagonal matrix. For the Hermitian c...
Definition: lalg_adv.F90:1297
subroutine zlowest_eigensolve(k, n, a, e, v, preserve_mat)
Computes the k lowest eigenvalues and the eigenvectors of a standard symmetric-definite eigenproblem,...
Definition: lalg_adv.F90:1376
complex(real64) function, public lalg_zduialpha(eigenvec, mmatrix, alpha, gamma, delta)
Definition: lalg_adv.F90:708
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...
Definition: lalg_adv.F90:878
subroutine dleast_squares_vec(nn, aa, bb, xx, preserve_mat)
Definition: lalg_adv.F90:3507
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...
Definition: lalg_adv.F90:2614
subroutine zcholesky(n, a, bof, err_code)
Compute the Cholesky decomposition of real symmetric or complex Hermitian positive definite matrix a,...
Definition: lalg_adv.F90:823
subroutine, public lalg_check_zeigenderivatives(n, mat)
Definition: lalg_adv.F90:584
complex(real64) function, public lalg_zdni(eigenvec, alpha, beta)
Definition: lalg_adv.F90:702
pure real(real64) function pseudoinverse_default_tolerance(m, n, sg_values)
Computes the default Moore-Penrose pseudoinverse tolerance for zeroing.
Definition: lalg_adv.F90:726
subroutine zleast_squares_vec(nn, aa, bb, xx, preserve_mat)
Definition: lalg_adv.F90:2006
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 ...
Definition: lalg_adv.F90:738
subroutine lalg_zgeev(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info)
Definition: lalg_adv.F90:286
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.
Definition: lalg_adv.F90:2152
integer function dmatrix_rank_svd(a, preserve_mat, tol)
Compute the rank of the matrix A using SVD.
Definition: lalg_adv.F90:3418
subroutine dinverse(n, a, method, det, threshold, uplo)
An interface to different method to invert a matrix.
Definition: lalg_adv.F90:3602
subroutine zeigensolve_parallel(n, a, e, bof, err_code)
Computes all the eigenvalues and the eigenvectors of a real symmetric or complex Hermitian eigenprobl...
Definition: lalg_adv.F90:2081
real(real64) function ddeterminant(n, a, preserve_mat)
Invert a real symmetric or complex Hermitian square matrix a.
Definition: lalg_adv.F90:2972
subroutine lalg_dgeev(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info)
Definition: lalg_adv.F90:259
subroutine dlowest_eigensolve(k, n, a, e, v, preserve_mat)
Computes the k lowest eigenvalues and the eigenvectors of a standard symmetric-definite eigenproblem,...
Definition: lalg_adv.F90:2883
subroutine dsvd_inverse(m, n, a, threshold)
Computes the inverse of a real M x N matrix, a, using the SVD decomposition.
Definition: lalg_adv.F90:3286
subroutine, public zlalg_phi(nn, pp, aa, ex, hermitian)
Definition: lalg_adv.F90:398
subroutine dlalg_pseudo_inverse(a, threshold)
Invert a matrix with the Moore-Penrose pseudo-inverse.
Definition: lalg_adv.F90:3352
subroutine zsvd_inverse(m, n, a, threshold)
Computes inverse of a complex MxN matrix, a, using the SVD decomposition.
Definition: lalg_adv.F90:1785
This module contains interfaces for LAPACK routines.
Definition: lapack.F90:120
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:117
This module contains interfaces for ScaLAPACK routines Interfaces are from http:
Definition: scalapack.F90:133
This module is intended to contain "only mathematical" functions and procedures.
Definition: sort.F90:119
This module is intended to contain simple general-purpose utility functions and procedures.
Definition: utils.F90:120
int true(void)