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
73
74 interface lalg_cholesky
75 module procedure dcholesky, zcholesky
76 end interface lalg_cholesky
77
78 interface lalg_geneigensolve
79 module procedure dgeneigensolve, zgeneigensolve
80 end interface lalg_geneigensolve
81
82 interface lalg_eigensolve_nonh
83 module procedure zeigensolve_nonh, deigensolve_nonh
84 end interface lalg_eigensolve_nonh
85
86 interface lalg_eigensolve
87 module procedure deigensolve, zeigensolve
88 end interface lalg_eigensolve
89
92 end interface lalg_eigensolve_tridiagonal
93
96 end interface lalg_eigensolve_parallel
97
100
101 interface lalg_determinant
102 module procedure ddeterminant, zdeterminant
103 end interface lalg_determinant
104
105 interface lalg_inverse
106 module procedure dinverse, zinverse
107 end interface lalg_inverse
108
109 interface lalg_linsyssolve
110 module procedure dlinsyssolve, zlinsyssolve
111 end interface lalg_linsyssolve
112
115 end interface lalg_singular_value_decomp
116
117 interface lalg_pseudo_inverse
119 end interface lalg_pseudo_inverse
120
121 interface lalg_svd_inverse
122 module procedure dsvd_inverse, zsvd_inverse
123 end interface lalg_svd_inverse
124
127 end interface lalg_lowest_geneigensolve
128
129 interface lalg_lowest_eigensolve
130 module procedure dlowest_eigensolve, zlowest_eigensolve
131 end interface lalg_lowest_eigensolve
132
133 interface lapack_geev
134 module procedure lalg_dgeev, lalg_zgeev
135 end interface lapack_geev
136
137 interface lalg_least_squares
138 module procedure dleast_squares_vec, zleast_squares_vec
139 end interface lalg_least_squares
140
141 interface lalg_matrix_norm2
142 module procedure dmatrix_norm2, zmatrix_norm2
143 end interface lalg_matrix_norm2
144
145 interface lalg_matrix_function
147 end interface lalg_matrix_function
148
149contains
150
151 ! ---------------------------------------------------------
153 real(real64) function sfmin()
154 interface
155 real(real64) function dlamch(cmach)
156 import real64
157 implicit none
158 character(1), intent(in) :: cmach
159 end function dlamch
160 end interface
161
162 sfmin = dlamch('S')
163 end function sfmin
164
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
174
175 real(real64), allocatable :: wr(:), wi(:)
176
177 push_sub(lalg_dgeev)
178
179 safe_allocate(wr(1:n))
180 safe_allocate(wi(1:n))
181
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)
185
186 safe_deallocate_a(wr)
187 safe_deallocate_a(wi)
188
189 pop_sub(lalg_dgeev)
190 end subroutine lalg_dgeev
191
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
201
202 push_sub(lalg_zgeev)
203
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)
205
206 pop_sub(lalg_zgeev)
207 end subroutine lalg_zgeev
208
226 subroutine zlalg_exp(nn, pp, aa, ex, hermitian)
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
232
233 complex(real64), allocatable :: evectors(:, :), zevalues(:)
234 real(real64), allocatable :: evalues(:)
235
236 integer :: ii
237
238 push_sub(zlalg_exp)
239
240 safe_allocate(evectors(1:nn, 1:nn))
241
242 if (hermitian) then
243 safe_allocate(evalues(1:nn))
244 safe_allocate(zevalues(1:nn))
245
246 evectors(1:nn, 1:nn) = aa(1:nn, 1:nn)
247
248 call lalg_eigensolve(nn, evectors, evalues)
249
250 zevalues(1:nn) = exp(pp*evalues(1:nn))
251
252 do ii = 1, nn
253 ex(1:nn, ii) = zevalues(1:nn)*conjg(evectors(ii, 1:nn))
254 end do
255
256 ex(:, :) = matmul(evectors(:, :), ex(:, :))
257
258 safe_deallocate_a(evalues)
259 safe_deallocate_a(zevalues)
260 else
261 safe_allocate(zevalues(1:nn))
262
263 evectors(1:nn, 1:nn) = aa(1:nn, 1:nn)
264
265 call lalg_eigensolve_nonh(nn, evectors, zevalues)
266
267 zevalues(1:nn) = exp(pp*zevalues(1:nn))
268
269 ex(1:nn, 1:nn) = evectors(1:nn, 1:nn)
270
271 call lalg_inverse(nn, evectors, 'dir')
272
273 do ii = 1, nn
274 evectors(1:nn, ii) = zevalues(1:nn)*evectors(1:nn, ii)
275 end do
276
277 ex(:, :) = matmul(ex(:, :), evectors(:, :))
278
279 safe_deallocate_a(zevalues)
280 end if
281
282 safe_deallocate_a(evectors)
283
284 pop_sub(zlalg_exp)
285 end subroutine zlalg_exp
286
287
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
310
311 complex(real64), allocatable :: evectors(:, :), zevalues(:)
312 real(real64), allocatable :: evalues(:)
313
314 integer :: ii
315
316 push_sub(zlalg_phi)
317
318 safe_allocate(evectors(1:nn, 1:nn))
320 if (hermitian) then
321 safe_allocate(evalues(1:nn))
322 safe_allocate(zevalues(1:nn))
323
324 evectors(:, :) = aa(:, :)
325
326 call lalg_eigensolve(nn, evectors, evalues)
327
328 do ii = 1, nn
329 zevalues(ii) = (exp(pp*evalues(ii)) - m_z1) / (pp*evalues(ii))
330 end do
331
332 do ii = 1, nn
333 ex(1:nn, ii) = zevalues(1:nn)*conjg(evectors(ii, 1:nn))
334 end do
335
336 ex(:, :) = matmul(evectors(:, :), ex(:, :))
337
338 safe_deallocate_a(evalues)
339 safe_deallocate_a(zevalues)
340 else
341 safe_allocate(zevalues(1:nn))
342
343 evectors(:, :) = aa(:, :)
344
345 call lalg_eigensolve_nonh(nn, evectors, zevalues)
346
347 do ii = 1, nn
348 zevalues(ii) = (exp(pp*zevalues(ii)) - m_z1) / (pp*zevalues(ii))
349 end do
350
351 ex(:, :) = evectors(:, :)
352
353 call lalg_inverse(nn, evectors, 'dir')
354
355 do ii = 1, nn
356 evectors(1:nn, ii) = zevalues(1:nn)*evectors(1:nn, ii)
357 end do
358
359 ex(:, :) = matmul(ex(:, :), evectors(:, :))
360
361 safe_deallocate_a(zevalues)
362 end if
363
364 pop_sub(zlalg_phi)
365 end subroutine zlalg_phi
366
367
379 subroutine lalg_zeigenderivatives(n, mat, zeigenvec, zeigenval, zmat)
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(:, :, :)
385
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(:, :)
391
392 push_sub(lalg_zeigenderivatives)
393
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))
398
399 zeigenvec = mat
400 call lalg_eigensolve_nonh(n, zeigenvec, zeigenval)
401
402 unit = m_z0
403 do i = 1, n
404 unit(i, i) = m_z1
405 end do
406
407 do i = 1, n
408 do alpha = 1, n
409 do beta = 1, n
410 knaught(alpha, beta) = - zeigenvec(alpha, i)*conjg(zeigenvec(beta, i))
411 end do
412 end do
413 knaught = knaught + unit
414 lambdaminusdm = zeigenval(i)*unit - mat
415 ! TODO(Alex) Swap this call
416 call lalg_zpseudoinverse(n, lambdaminusdm, ilambdaminusdm)
417 zmat(:, :, i) = matmul(ilambdaminusdm, knaught)
418 end do
419
420 safe_deallocate_a(unit)
421 safe_deallocate_a(knaught)
422 safe_deallocate_a(lambdaminusdm)
423 safe_deallocate_a(ilambdaminusdm)
425 end subroutine lalg_zeigenderivatives
426
427
428 ! TODO(Alex) Issue #1132.
429 ! * This routine is redunant. Its usage should be replaced with either zsvd_inverse or zlalg_pseudo_inverse.
431 subroutine lalg_zpseudoinverse(n, mat, imat)
432 integer, intent(in) :: n
433 complex(real64), contiguous, intent(in) :: mat(:, :)
434 complex(real64), contiguous, intent(out) :: imat(:, :)
435
436 integer :: i
437 complex(real64), allocatable :: u(:, :), vt(:, :), sigma(:, :)
438 real(real64), allocatable :: sg_values(:)
439
440 push_sub(lalg_zpseudoinverse)
441
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))
446
447 imat = mat
448 call lalg_singular_value_decomp(n, n, imat, u, vt, sg_values)
449
450 sigma = m_z0
451 do i = 1, n
452 if (abs(sg_values(i)) <= 1.0e-12_real64 * maxval(abs(sg_values))) then
453 sigma(i, i) = m_z0
454 else
455 sigma(i, i) = m_z1 / sg_values(i)
456 end if
457 end do
458
459 vt(:, :) = conjg(transpose(vt(:, :)))
460 u(:, :) = conjg(transpose(u(:, :)))
461 imat = matmul(vt, matmul(sigma, u))
462
463 ! Check if we truly have a pseudoinverse
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))
467 write(*, *) vt
468 write(*, *)
469 write(*, *) 1.0e-10_real64 * maxval(abs(mat))
470 write(*, *) maxval(abs(vt)) > 1.0e-10_real64 * maxval(abs(mat))
471 write(*, *) mat
472 write(message(1), '(a)') 'Pseudoinverse failed.'
473 call messages_fatal(1)
474 end if
475
476 safe_deallocate_a(u)
477 safe_deallocate_a(vt)
478 safe_deallocate_a(sigma)
479 safe_deallocate_a(sg_values)
480 pop_sub(lalg_zpseudoinverse)
481 end subroutine lalg_zpseudoinverse
482
483
490 subroutine lalg_check_zeigenderivatives(n, mat)
491 integer, intent(in) :: n
492 complex(real64), intent(in) :: mat(:, :)
493
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(:)
498
500
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))
511
512 assert(n == 2)
513
514 dm(:, :) = mat(:, :)
515 call lalg_zeigenderivatives(2, dm, zeigref_, zeigenval, mmatrix)
516
517
518 deltah = (0.000001_real64, 0.000_real64)
519 !deltah = M_z1 * maxval(abs(dm)) * 0.001_real64
520 do alpha = 1, n
521 do beta = 1, n
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
527 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigplus)
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)
531
532 zeigenvec(:, :) = dm(:, :)
533 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
534 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigminus)
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)
538
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)
542
543 do gamma = 1, n
544 do delta = 1, n
545 if (alpha == gamma .and. beta == delta) then
546 zder_direct = lalg_zd2ni(zeigref_(:, 1), mmatrix(:, :, 1), alpha, beta, gamma, delta)
547
548 zeigenvec(:, :) = dm(:, :)
549 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeig0)
550
551 zeigenvec(:, :) = dm(:, :)
552 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
553 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigplus)
554
555 zeigenvec(:, :) = dm(:, :)
556 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
557 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigminus)
558
559 write(*, '(4i1,4f24.12)') alpha, beta, gamma, delta, &
560 zder_direct, (zeigplus(1) + zeigminus(1) - m_two*zeig0(1))/(deltah**2)
561 else
562 zder_direct = lalg_zd2ni(zeigref_(:, 1), mmatrix(:, :, 1), alpha, beta, gamma, delta)
563
564 zeigenvec(:, :) = dm(:, :)
565 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
566 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) + deltah
567 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigplus)
568
569 zeigenvec(:, :) = dm(:, :)
570 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
571 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) - deltah
572 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigminus)
573
574 zeigenvec(:, :) = dm(:, :)
575 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
576 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) - deltah
577 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigplusminus)
578
579 zeigenvec(:, :) = dm(:, :)
580 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
581 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) + deltah
582 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigminusplus)
584 write(*, '(4i1,4f24.12)') alpha, beta, gamma, delta, &
585 zder_direct, (zeigplus(1) + zeigminus(1) - zeigplusminus(1) - zeigminusplus(1))/(m_four*deltah**2)
586
587
588 end if
589 end do
590 end do
591
592 end do
593 end do
594
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)
606 end subroutine lalg_check_zeigenderivatives
607
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)
612 end function lalg_zdni
613
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)
617 lalg_zduialpha = mmatrix(alpha, gamma) * eigenvec(delta)
618 end function lalg_zduialpha
619
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)
625 end function lalg_zd2ni
626
627 ! ---------------------------------------------------------
632 pure real(real64) function pseudoinverse_default_tolerance(m, n, sg_values) result(tol)
633 integer, intent(in) :: m, n
634 real(real64), intent(in) :: sg_values(:)
635 tol = m_epsilon * m * n * maxval(sg_values)
637
638
644 function lalg_remove_rotation(n, A) result(P)
645 integer, intent(in) :: n
646 real(real64), intent(in) :: a(1:n, 1:n)
647 real(real64) :: p(1:n, 1:n)
648
649 real(real64) :: ata(1:n, 1:n)
650
651 push_sub(lalg_remove_rotation)
652
653 ata = matmul(transpose(a), a)
654 call lalg_matrix_function(n, m_one, ata, p, square_root, .true.)
655
656 pop_sub(lalg_remove_rotation)
657 end function lalg_remove_rotation
658
659#include "undef.F90"
660#include "complex.F90"
661#include "lalg_adv_lapack_inc.F90"
662
663#include "undef.F90"
664#include "real.F90"
665#include "lalg_adv_lapack_inc.F90"
666
667end module lalg_adv_oct_m
668
669!! Local Variables:
670!! mode: f90
671!! coding: utf-8
672!! 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:118
subroutine zmatrix_norm2(m, n, a, norm_l2, preserve_mat)
Norm of a 2D matrix.
Definition: lalg_adv.F90:1582
complex(real64) function, public lalg_zd2ni(eigenvec, mmatrix, alpha, beta, gamma, delta)
Definition: lalg_adv.F90:714
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:1667
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:3689
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:2499
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:1894
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:3360
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:2605
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:2914
subroutine zinverse(n, a, method, det, threshold, uplo)
An interface to different method to invert a matrix.
Definition: lalg_adv.F90:2214
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:3282
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:1745
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:2444
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:3870
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:2844
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 dmatrix_norm2(m, n, a, norm_l2, preserve_mat)
Norm of a 2D matrix.
Definition: lalg_adv.F90:3197
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:3617
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:2724
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:2262
subroutine dinverse(n, a, method, det, threshold, uplo)
An interface to different method to invert a matrix.
Definition: lalg_adv.F90:3822
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:3082
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:2993
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:3439
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:3505
subroutine zsvd_inverse(m, n, a, threshold)
Computes inverse of a complex MxN matrix, a, using the SVD decomposition.
Definition: lalg_adv.F90:1828
This module contains interfaces for LAPACK routines.
Definition: lapack.F90:118
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:115
This module contains interfaces for ScaLAPACK routines Interfaces are from http:
Definition: scalapack.F90:131
This module is intended to contain "only mathematical" functions and procedures.
Definition: sort.F90:117
This module is intended to contain simple general-purpose utility functions and procedures.
Definition: utils.F90:118
int true(void)