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
23 use blacs_oct_m
24 use debug_oct_m
25 use global_oct_m
26 use, intrinsic :: iso_fortran_env
27 use lapack_oct_m
29 use math_oct_m
31 use mpi_oct_m
34 use sort_oct_m
35 use utils_oct_m
36#ifdef HAVE_ELPA
37 use elpa
38#endif
39
40 implicit none
41
42 private
43 public :: &
57 zlalg_exp, &
58 zlalg_phi, &
63 lalg_zdni, &
65 lalg_zd2ni, &
69
70
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_parallel
91
94
95 interface lalg_determinant
96 module procedure ddeterminant, zdeterminant
97 end interface lalg_determinant
98
99 interface lalg_inverse
100 module procedure dinverse, zinverse
101 end interface lalg_inverse
102
103 interface lalg_linsyssolve
104 module procedure dlinsyssolve, zlinsyssolve
105 end interface lalg_linsyssolve
106
109 end interface lalg_singular_value_decomp
110
111 interface lalg_pseudo_inverse
113 end interface lalg_pseudo_inverse
115 interface lalg_svd_inverse
116 module procedure dsvd_inverse, zsvd_inverse
117 end interface lalg_svd_inverse
118
121 end interface lalg_lowest_geneigensolve
122
123 interface lalg_lowest_eigensolve
124 module procedure dlowest_eigensolve, zlowest_eigensolve
125 end interface lalg_lowest_eigensolve
126
127 interface lapack_geev
128 module procedure lalg_dgeev, lalg_zgeev
129 end interface lapack_geev
130
131 interface lalg_least_squares
132 module procedure dleast_squares_vec, zleast_squares_vec
133 end interface lalg_least_squares
134
135 interface lalg_matrix_norm2
136 module procedure dmatrix_norm2, zmatrix_norm2
137 end interface lalg_matrix_norm2
138
139 interface lalg_matrix_function
141 end interface lalg_matrix_function
142
143contains
144
145 ! ---------------------------------------------------------
147 real(real64) function sfmin()
148 interface
149 real(real64) function dlamch(cmach)
150 import real64
151 implicit none
152 character(1), intent(in) :: cmach
153 end function dlamch
154 end interface
155
156 sfmin = dlamch('S')
157 end function sfmin
158
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
168
169 real(real64), allocatable :: wr(:), wi(:)
170
171 push_sub(lalg_dgeev)
172
173 safe_allocate(wr(1:n))
174 safe_allocate(wi(1:n))
175
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)
179
180 safe_deallocate_a(wr)
181 safe_deallocate_a(wi)
182
183 pop_sub(lalg_dgeev)
184 end subroutine lalg_dgeev
185
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
195
196 push_sub(lalg_zgeev)
197
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)
199
200 pop_sub(lalg_zgeev)
201 end subroutine lalg_zgeev
202
220 subroutine zlalg_exp(nn, pp, aa, ex, hermitian)
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
226
227 complex(real64), allocatable :: evectors(:, :), zevalues(:)
228 real(real64), allocatable :: evalues(:)
229
230 integer :: ii
231
232 push_sub(zlalg_exp)
233
234 safe_allocate(evectors(1:nn, 1:nn))
235
236 if (hermitian) then
237 safe_allocate(evalues(1:nn))
238 safe_allocate(zevalues(1:nn))
239
240 evectors(1:nn, 1:nn) = aa(1:nn, 1:nn)
241
242 call lalg_eigensolve(nn, evectors, evalues)
243
244 zevalues(1:nn) = exp(pp*evalues(1:nn))
245
246 do ii = 1, nn
247 ex(1:nn, ii) = zevalues(1:nn)*conjg(evectors(ii, 1:nn))
248 end do
249
250 ex(1:nn, 1:nn) = matmul(evectors(1:nn, 1:nn), ex(1:nn, 1:nn))
251
252 safe_deallocate_a(evalues)
253 safe_deallocate_a(zevalues)
254 else
255 safe_allocate(zevalues(1:nn))
256
257 evectors(1:nn, 1:nn) = aa(1:nn, 1:nn)
258
259 call lalg_eigensolve_nonh(nn, evectors, zevalues)
260
261 zevalues(1:nn) = exp(pp*zevalues(1:nn))
262
263 ex(1:nn, 1:nn) = evectors(1:nn, 1:nn)
264
265 call lalg_inverse(nn, evectors, 'dir')
266
267 do ii = 1, nn
268 evectors(1:nn, ii) = zevalues(1:nn)*evectors(1:nn, ii)
269 end do
270
271 ex(1:nn, 1:nn) = matmul(ex(1:nn, 1:nn), evectors(1:nn, 1:nn))
272
273 safe_deallocate_a(zevalues)
274 end if
275
276 safe_deallocate_a(evectors)
277
278 pop_sub(zlalg_exp)
279 end subroutine zlalg_exp
280
281
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
304
305 complex(real64), allocatable :: evectors(:, :), zevalues(:)
306 real(real64), allocatable :: evalues(:)
307
308 integer :: ii
309
310 push_sub(zlalg_phi)
311
312 safe_allocate(evectors(1:nn, 1:nn))
314 if (hermitian) then
315 safe_allocate(evalues(1:nn))
316 safe_allocate(zevalues(1:nn))
317
318 evectors = aa
319
320 call lalg_eigensolve(nn, evectors, evalues)
321
322 do ii = 1, nn
323 zevalues(ii) = (exp(pp*evalues(ii)) - m_z1) / (pp*evalues(ii))
324 end do
325
326 do ii = 1, nn
327 ex(1:nn, ii) = zevalues(1:nn)*conjg(evectors(ii, 1:nn))
328 end do
329
330 ex(1:nn, 1:nn) = matmul(evectors(1:nn, 1:nn), ex(1:nn, 1:nn))
331
332 safe_deallocate_a(evalues)
333 safe_deallocate_a(zevalues)
334 else
335 safe_allocate(zevalues(1:nn))
336
337 evectors(1:nn, 1:nn) = aa(1:nn, 1:nn)
338
339 call lalg_eigensolve_nonh(nn, evectors, zevalues)
340
341 do ii = 1, nn
342 zevalues(ii) = (exp(pp*zevalues(ii)) - m_z1) / (pp*zevalues(ii))
343 end do
344
345 ex(1:nn, 1:nn) = evectors(1:nn, 1:nn)
346
347 call lalg_inverse(nn, evectors, 'dir')
348
349 do ii = 1, nn
350 evectors(1:nn, ii) = zevalues(1:nn)*evectors(1:nn, ii)
351 end do
352
353 ex(1:nn, 1:nn) = matmul(ex(1:nn, 1:nn), evectors(1:nn, 1:nn))
354
355 safe_deallocate_a(zevalues)
356 end if
357
358 pop_sub(zlalg_phi)
359 end subroutine zlalg_phi
360
361
373 subroutine lalg_zeigenderivatives(n, mat, zeigenvec, zeigenval, zmat)
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(:, :, :)
379
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(:, :)
385
386 push_sub(lalg_zeigenderivatives)
387
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))
392
393 zeigenvec = mat
394 call lalg_eigensolve_nonh(n, zeigenvec, zeigenval)
395
396 unit = m_z0
397 do i = 1, n
398 unit(i, i) = m_z1
399 end do
400
401 do i = 1, n
402 do alpha = 1, n
403 do beta = 1, n
404 knaught(alpha, beta) = - zeigenvec(alpha, i)*conjg(zeigenvec(beta, i))
405 end do
406 end do
407 knaught = knaught + unit
408 lambdaminusdm = zeigenval(i)*unit - mat
409 ! TODO(Alex) Swap this call
410 call lalg_zpseudoinverse(n, lambdaminusdm, ilambdaminusdm)
411 zmat(:, :, i) = matmul(ilambdaminusdm, knaught)
412 end do
413
414 safe_deallocate_a(unit)
415 safe_deallocate_a(knaught)
416 safe_deallocate_a(lambdaminusdm)
417 safe_deallocate_a(ilambdaminusdm)
419 end subroutine lalg_zeigenderivatives
420
421
422 ! TODO(Alex) Issue #1132.
423 ! * This routine is redunant. Its usage should be replaced with either zsvd_inverse or zlalg_pseudo_inverse.
425 subroutine lalg_zpseudoinverse(n, mat, imat)
426 integer, intent(in) :: n
427 complex(real64), intent(in) :: mat(:, :)
428 complex(real64), intent(out) :: imat(:, :)
429
430 integer :: i
431 complex(real64), allocatable :: u(:, :), vt(:, :), sigma(:, :)
432 real(real64), allocatable :: sg_values(:)
433
434 push_sub(lalg_zpseudoinverse)
435
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))
440
441 imat = mat
442 call lalg_singular_value_decomp(n, n, imat, u, vt, sg_values)
443
444 sigma = m_z0
445 do i = 1, n
446 if (abs(sg_values(i)) <= 1.0e-12_real64 * maxval(abs(sg_values))) then
447 sigma(i, i) = m_z0
448 else
449 sigma(i, i) = m_z1 / sg_values(i)
450 end if
451 end do
452
453 vt = conjg(transpose(vt))
454 u = conjg(transpose(u))
455 imat = matmul(vt, matmul(sigma, u))
456
457 ! Check if we truly have a pseudoinverse
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))
461 write(*, *) vt
462 write(*, *)
463 write(*, *) 1.0e-10_real64 * maxval(abs(mat))
464 write(*, *) maxval(abs(vt)) > 1.0e-10_real64 * maxval(abs(mat))
465 write(*, *) mat
466 write(message(1), '(a)') 'Pseudoinverse failed.'
467 call messages_fatal(1)
468 end if
469
470 safe_deallocate_a(u)
471 safe_deallocate_a(vt)
472 safe_deallocate_a(sigma)
473 safe_deallocate_a(sg_values)
474 pop_sub(lalg_zpseudoinverse)
475 end subroutine lalg_zpseudoinverse
476
477
484 subroutine lalg_check_zeigenderivatives(n, mat)
485 integer, intent(in) :: n
486 complex(real64), intent(in) :: mat(:, :)
487
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(:)
492
494
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))
505
506 assert(n == 2)
507
508 dm = mat
509 call lalg_zeigenderivatives(2, dm, zeigref_, zeigenval, mmatrix)
510
511
512 deltah = (0.000001_real64, 0.000_real64)
513 !deltah = M_z1 * maxval(abs(dm)) * 0.001_real64
514 do alpha = 1, n
515 do beta = 1, n
516 zder_direct = lalg_zdni(zeigref_(:, 2), alpha, beta)
517 zuder_direct = lalg_zduialpha(zeigref_(:, 2), mmatrix(:, :, 2), 2, alpha, beta)
519 zeigenvec = dm
520 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
521 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigplus)
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)
525
526 zeigenvec = dm
527 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
528 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigminus)
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)
532
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)
536
537 do gamma = 1, n
538 do delta = 1, n
539 if (alpha == gamma .and. beta == delta) then
540 zder_direct = lalg_zd2ni(zeigref_(:, 1), mmatrix(:, :, 1), alpha, beta, gamma, delta)
541
542 zeigenvec = dm
543 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeig0)
544
545 zeigenvec = dm
546 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
547 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigplus)
548
549 zeigenvec = dm
550 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
551 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigminus)
552
553 write(*, '(4i1,4f24.12)') alpha, beta, gamma, delta, &
554 zder_direct, (zeigplus(1) + zeigminus(1) - m_two*zeig0(1))/(deltah**2)
555 else
556 zder_direct = lalg_zd2ni(zeigref_(:, 1), mmatrix(:, :, 1), alpha, beta, gamma, delta)
557
558 zeigenvec = dm
559 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
560 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) + deltah
561 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigplus)
562
563 zeigenvec = dm
564 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
565 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) - deltah
566 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigminus)
567
568 zeigenvec = dm
569 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah
570 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) - deltah
571 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigplusminus)
572
573 zeigenvec = dm
574 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah
575 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) + deltah
576 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigminusplus)
578 write(*, '(4i1,4f24.12)') alpha, beta, gamma, delta, &
579 zder_direct, (zeigplus(1) + zeigminus(1) - zeigplusminus(1) - zeigminusplus(1))/(m_four*deltah**2)
580
581
582 end if
583 end do
584 end do
585
586 end do
587 end do
588
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)
600 end subroutine lalg_check_zeigenderivatives
601
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)
606 end function lalg_zdni
607
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)
611 lalg_zduialpha = mmatrix(alpha, gamma) * eigenvec(delta)
612 end function lalg_zduialpha
613
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)
619 end function lalg_zd2ni
620
621 ! ---------------------------------------------------------
626 pure real(real64) function pseudoinverse_default_tolerance(m, n, sg_values) result(tol)
627 integer, intent(in) :: m, n
628 real(real64), intent(in) :: sg_values(:)
629 tol = m_epsilon * m * n * maxval(sg_values)
631
632
638 function lalg_remove_rotation(n, A) result(P)
639 integer, intent(in) :: n
640 real(real64), intent(in) :: a(1:n, 1:n)
641 real(real64) :: p(1:n, 1:n)
642
643 real(real64) :: ata(1:n, 1:n)
644
645 push_sub(lalg_remove_rotation)
646
647 ata = matmul(transpose(a), a)
648 call lalg_matrix_function(n, m_one, ata, p, square_root, .true.)
649
650 pop_sub(lalg_remove_rotation)
651 end function lalg_remove_rotation
652
653#include "undef.F90"
654#include "complex.F90"
655#include "lalg_adv_lapack_inc.F90"
656
657#include "undef.F90"
658#include "real.F90"
659#include "lalg_adv_lapack_inc.F90"
660
661end module lalg_adv_oct_m
662
663!! Local Variables:
664!! mode: f90
665!! coding: utf-8
666!! End:
Note that lalg_determinant and lalg_inverse are just wrappers over the same routine.
Definition: lalg_adv.F90:188
lalg_gemm with both the (Hermitian) transpose of A and B.
Definition: lalg_basic.F90:252
scales a vector by a constant
Definition: lalg_basic.F90:157
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.
subroutine zmatrix_norm2(m, n, a, norm_l2, preserve_mat)
Norm of a 2D matrix.
Definition: lalg_adv.F90:1488
complex(real64) function, public lalg_zd2ni(eigenvec, mmatrix, alpha, beta, gamma, delta)
Definition: lalg_adv.F90:708
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:1091
complex(real64) function zdeterminant(n, a, preserve_mat)
Invert a real symmetric or complex Hermitian square matrix a.
Definition: lalg_adv.F90:1373
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:972
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:1573
subroutine, public zlalg_exp(nn, pp, aa, ex, hermitian)
Definition: lalg_adv.F90:314
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:3490
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:2381
subroutine lalg_zpseudoinverse(n, mat, imat)
Computes the Moore-Penrose pseudoinverse of a complex matrix.
Definition: lalg_adv.F90:519
subroutine zlalg_pseudo_inverse(a, threshold)
Invert a matrix with the Moore-Penrose pseudo-inverse.
Definition: lalg_adv.F90:1800
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:1210
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:3161
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:2483
subroutine zinverse(n, a, method, det, threshold, uplo)
An interface to different method to invert a matrix.
Definition: lalg_adv.F90:2120
real(real64) function sfmin()
Auxiliary function.
Definition: lalg_adv.F90:241
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:3083
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:1651
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:2328
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:2722
subroutine, public lalg_zeigenderivatives(n, mat, zeigenvec, zeigenval, zmat)
Definition: lalg_adv.F90:467
subroutine dmatrix_norm2(m, n, a, norm_l2, preserve_mat)
Norm of a 2D matrix.
Definition: lalg_adv.F90:2998
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:1282
complex(real64) function, public lalg_zduialpha(eigenvec, mmatrix, alpha, gamma, delta)
Definition: lalg_adv.F90:702
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:870
subroutine dleast_squares_vec(nn, aa, bb, xx, preserve_mat)
Definition: lalg_adv.F90:3418
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:2602
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:817
subroutine, public lalg_check_zeigenderivatives(n, mat)
Definition: lalg_adv.F90:578
complex(real64) function, public lalg_zdni(eigenvec, alpha, beta)
Definition: lalg_adv.F90:696
pure real(real64) function pseudoinverse_default_tolerance(m, n, sg_values)
Computes the default Moore-Penrose pseudoinverse tolerance for zeroing.
Definition: lalg_adv.F90:720
subroutine zleast_squares_vec(nn, aa, bb, xx, preserve_mat)
Definition: lalg_adv.F90:1912
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:732
subroutine lalg_zgeev(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info)
Definition: lalg_adv.F90:280
subroutine dlalg_matrix_function(n, factor, a, fun_a, fun, hermitian)
This routine calculates a function of a matrix by using an eigenvalue decomposition.
Definition: lalg_adv.F90:3671
subroutine dinverse(n, a, method, det, threshold, uplo)
An interface to different method to invert a matrix.
Definition: lalg_adv.F90:3623
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:1987
real(real64) function ddeterminant(n, a, preserve_mat)
Invert a real symmetric or complex Hermitian square matrix a.
Definition: lalg_adv.F90:2883
subroutine lalg_dgeev(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info)
Definition: lalg_adv.F90:253
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.
Definition: lalg_adv.F90:2168
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:2794
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:3240
subroutine, public zlalg_phi(nn, pp, aa, ex, hermitian)
Definition: lalg_adv.F90:392
subroutine dlalg_pseudo_inverse(a, threshold)
Invert a matrix with the Moore-Penrose pseudo-inverse.
Definition: lalg_adv.F90:3306
subroutine zsvd_inverse(m, n, a, threshold)
Computes inverse of a complex MxN matrix, a, using the SVD decomposition.
Definition: lalg_adv.F90:1734
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)