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 global_oct_m 23 use lapack_oct_m 24 use messages_oct_m 25 use profiling_oct_m 26 use sort_oct_m 27 use utils_oct_m 28 29 implicit none 30 31 private 32 public :: & 33 lalg_cholesky, & 34 lalg_geneigensolve, & 35 lalg_eigensolve, & 36 lalg_eigensolve_nonh, & 37 lalg_determinant, & 38 lalg_inverter, & 39 lalg_sym_inverter, & 40 lalg_linsyssolve, & 41 lalg_singular_value_decomp, & 42 lalg_svd_inverse, & 43 lalg_invert_upper_triangular, & 44 lalg_lowest_geneigensolve, & 45 lalg_lowest_eigensolve, & 46 zlalg_exp, & 47 zlalg_phi, & 48 lalg_zpseudoinverse, & 49 lalg_zeigenderivatives, & 50 lalg_check_zeigenderivatives, & 51 lalg_zdni, & 52 lalg_zduialpha, & 53 lalg_zd2ni, & 54 lalg_least_squares 55 56 type(profile_t), save :: cholesky_prof, eigensolver_prof 57 58 interface lalg_cholesky 59 module procedure dcholesky, zcholesky 60 end interface lalg_cholesky 61 62 interface lalg_geneigensolve 63 module procedure dgeneigensolve, zgeneigensolve 64 end interface lalg_geneigensolve 65 66 interface lalg_eigensolve_nonh 67 module procedure zeigensolve_nonh, deigensolve_nonh 68 end interface lalg_eigensolve_nonh 69 70 interface lalg_eigensolve 71 module procedure deigensolve, zeigensolve 72 end interface lalg_eigensolve 73 74 interface lalg_determinant 75 module procedure ddeterminant, zdeterminant 76 end interface lalg_determinant 77 78 interface lalg_inverter 79 module procedure dinverter, zinverter 80 end interface lalg_inverter 81 82 interface lalg_sym_inverter 83 module procedure dsym_inverter, zsym_inverter 84 end interface lalg_sym_inverter 85 86 interface lalg_linsyssolve 87 module procedure dlinsyssolve, zlinsyssolve 88 end interface lalg_linsyssolve 89 90 interface lalg_singular_value_decomp 91 module procedure dsingular_value_decomp, zsingular_value_decomp 92 end interface lalg_singular_value_decomp 93 94 interface lalg_svd_inverse 95 module procedure dsvd_inverse, zsvd_inverse 96 end interface lalg_svd_inverse 97 98 interface lalg_invert_upper_triangular 99 module procedure dinvert_upper_triangular, zinvert_upper_triangular 100 end interface lalg_invert_upper_triangular 101 102 interface lalg_lowest_geneigensolve 103 module procedure dlowest_geneigensolve, zlowest_geneigensolve 104 end interface lalg_lowest_geneigensolve 105 106 interface lalg_lowest_eigensolve 107 module procedure dlowest_eigensolve, zlowest_eigensolve 108 end interface lalg_lowest_eigensolve 109 110 interface lapack_geev 111 module procedure lalg_dgeev, lalg_zgeev 112 end interface lapack_geev 113 114 interface lalg_least_squares 115 module procedure dleast_squares_vec, zleast_squares_vec 116 end interface lalg_least_squares 117 118contains 119 120 ! --------------------------------------------------------- 121 !> Auxiliary function 122 FLOAT function sfmin() 123 interface 124 FLOAT function dlamch(cmach) 125 implicit none 126 character(1), intent(in) :: cmach 127 end function dlamch 128 end interface 129 130 sfmin = dlamch('S') 131 end function sfmin 132 133 subroutine lalg_dgeev(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info) 134 character(1), intent(in) :: jobvl, jobvr 135 integer, intent(in) :: n, lda, ldvl, ldvr, lwork 136 real(8), intent(inout) :: a(:,:) !< a(lda,n) 137 complex(8), intent(out) :: w(:) !< w(n) 138 real(8), intent(out) :: vl(:,:), vr(:,:) !< vl(ldvl,n), vl(ldvr,n) 139 real(8), intent(out) :: rwork(:) !< rwork(max(1,2n)) 140 real(8), intent(out) :: work(:) !< work(lwork) 141 integer, intent(out) :: info 142 143 FLOAT, allocatable :: wr(:), wi(:) 144 145 PUSH_SUB(lalg_dgeev) 146 147 SAFE_ALLOCATE(wr(1:n)) 148 SAFE_ALLOCATE(wi(1:n)) 149 150 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) 151 152 w(1:n) = TOCMPLX(wr(1:n), wi(1:n)) 153 154 SAFE_DEALLOCATE_A(wr) 155 SAFE_DEALLOCATE_A(wi) 156 157 POP_SUB(lalg_dgeev) 158 end subroutine lalg_dgeev 159 160 subroutine lalg_zgeev(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info) 161 character(1), intent(in) :: jobvl, jobvr 162 integer, intent(in) :: n, lda, ldvl, ldvr, lwork 163 complex(8), intent(inout) :: a(:,:) !< a(lda,n) 164 complex(8), intent(out) :: w(:) !< w(n) 165 complex(8), intent(out) :: vl(:,:), vr(:,:) !< vl(ldvl,n), vl(ldvr,n) 166 real(8), intent(out) :: rwork(:) !< rwork(max(1,2n)) 167 complex(8), intent(out) :: work(:) !< work(lwork) 168 integer, intent(out) :: info 169 170 PUSH_SUB(lalg_zgeev) 171 172 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) 173 174 POP_SUB(lalg_zgeev) 175 end subroutine lalg_zgeev 176 177!>------------------------------------------------- 178 !! 179 !! This routine calculates the exponential of a matrix by using an 180 !! eigenvalue decomposition. 181 !! 182 !! For the hermitian case: 183 !! 184 !! A = V D V^T => exp(A) = V exp(D) V^T 185 !! 186 !! and in general 187 !! 188 !! A = V D V^-1 => exp(A) = V exp(D) V^-1 189 !! 190 !! This is slow but it is simple to implement, and for the moment it 191 !! does not affect performance. 192 !! 193 !!--------------------------------------------- 194 subroutine zlalg_exp(nn, pp, aa, ex, hermitian) 195 integer, intent(in) :: nn 196 CMPLX, intent(in) :: pp 197 CMPLX, intent(in) :: aa(:, :) 198 CMPLX, intent(inout) :: ex(:, :) 199 logical, intent(in) :: hermitian 200 201 CMPLX, allocatable :: evectors(:, :), zevalues(:) 202 FLOAT, allocatable :: evalues(:) 203 204 integer :: ii 205 206 PUSH_SUB(zlalg_exp) 207 208 SAFE_ALLOCATE(evectors(1:nn, 1:nn)) 209 210 if(hermitian) then 211 SAFE_ALLOCATE(evalues(1:nn)) 212 SAFE_ALLOCATE(zevalues(1:nn)) 213 214 evectors(1:nn, 1:nn) = aa(1:nn, 1:nn) 215 216 call lalg_eigensolve(nn, evectors, evalues) 217 218 zevalues(1:nn) = exp(pp*evalues(1:nn)) 219 220 do ii = 1, nn 221 ex(1:nn, ii) = zevalues(1:nn)*conjg(evectors(ii, 1:nn)) 222 end do 223 224 ex(1:nn, 1:nn) = matmul(evectors(1:nn, 1:nn), ex(1:nn, 1:nn)) 225 226 SAFE_DEALLOCATE_A(evalues) 227 SAFE_DEALLOCATE_A(zevalues) 228 else 229 SAFE_ALLOCATE(zevalues(1:nn)) 230 231 evectors(1:nn, 1:nn) = aa(1:nn, 1:nn) 232 233 call lalg_eigensolve_nonh(nn, evectors, zevalues) 234 235 zevalues(1:nn) = exp(pp*zevalues(1:nn)) 236 237 ex(1:nn, 1:nn) = evectors(1:nn, 1:nn) 238 239 call lalg_inverter(nn, evectors) 240 241 do ii = 1, nn 242 evectors(1:nn, ii) = zevalues(1:nn)*evectors(1:nn, ii) 243 end do 244 245 ex(1:nn, 1:nn) = matmul(ex(1:nn, 1:nn), evectors(1:nn, 1:nn)) 246 247 SAFE_DEALLOCATE_A(zevalues) 248 end if 249 250 SAFE_DEALLOCATE_A(evectors) 251 252 POP_SUB(zlalg_exp) 253 end subroutine zlalg_exp 254 255 256 !>------------------------------------------------- 257 !! 258 !! This routine calculates phi(pp*A), where A is a matrix, 259 !! pp is any complex number, and phi is the function: 260 !! 261 !! phi(x) = (e^x - 1)/x 262 !! 263 !! For the Hermitian case, for any function f: 264 !! 265 !! A = V D V^T => f(A) = V f(D) V^T 266 !! 267 !! and in general 268 !! 269 !! A = V D V^-1 => f(A) = V f(D) V^-1 270 !! 271 !!--------------------------------------------- 272 subroutine zlalg_phi(nn, pp, aa, ex, hermitian) 273 integer, intent(in) :: nn 274 CMPLX, intent(in) :: pp 275 CMPLX, intent(in) :: aa(:, :) 276 CMPLX, intent(inout) :: ex(:, :) 277 logical, intent(in) :: hermitian 278 279 CMPLX, allocatable :: evectors(:, :), zevalues(:) 280 FLOAT, allocatable :: evalues(:) 281 282 integer :: ii 283 284 PUSH_SUB(zlalg_phi) 285 286 SAFE_ALLOCATE(evectors(1:nn, 1:nn)) 287 288 if(hermitian) then 289 SAFE_ALLOCATE(evalues(1:nn)) 290 SAFE_ALLOCATE(zevalues(1:nn)) 291 292 evectors = aa 293 294 call lalg_eigensolve(nn, evectors, evalues) 295 296 do ii = 1, nn 297 zevalues(ii) = (exp(pp*evalues(ii)) - M_z1) / (pp*evalues(ii)) 298 end do 299 300 do ii = 1, nn 301 ex(1:nn, ii) = zevalues(1:nn)*conjg(evectors(ii, 1:nn)) 302 end do 303 304 ex(1:nn, 1:nn) = matmul(evectors(1:nn, 1:nn), ex(1:nn, 1:nn)) 305 306 SAFE_DEALLOCATE_A(evalues) 307 SAFE_DEALLOCATE_A(zevalues) 308 else 309 SAFE_ALLOCATE(zevalues(1:nn)) 310 311 evectors(1:nn, 1:nn) = aa(1:nn, 1:nn) 312 313 call lalg_eigensolve_nonh(nn, evectors, zevalues) 314 315 do ii = 1, nn 316 zevalues(ii) = (exp(pp*zevalues(ii)) - M_z1) / (pp*zevalues(ii)) 317 end do 318 319 ex(1:nn, 1:nn) = evectors(1:nn, 1:nn) 320 321 call lalg_inverter(nn, evectors) 322 323 do ii = 1, nn 324 evectors(1:nn, ii) = zevalues(1:nn)*evectors(1:nn, ii) 325 end do 326 327 ex(1:nn, 1:nn) = matmul(ex(1:nn, 1:nn), evectors(1:nn, 1:nn)) 328 329 SAFE_DEALLOCATE_A(zevalues) 330 end if 331 332 POP_SUB(zlalg_phi) 333 end subroutine zlalg_phi 334 335 336 !>------------------------------------------------- 337 !! Computes the necessary ingredients to obtain, 338 !! later, the first and second derivatives of the 339 !! eigenvalues of a Hermitean complex matrix zmat, 340 !! and the first derivatives of the eigenvectors. 341 !! 342 !! This follows the scheme of J. R. Magnus, 343 !! Econometric Theory 1, 179 (1985), restricted to 344 !! Hermitean matrices, although probably this can be 345 !! found in other sources. 346 !!--------------------------------------------- 347 subroutine lalg_zeigenderivatives(n, mat, zeigenvec, zeigenval, zmat) 348 integer, intent(in) :: n 349 CMPLX, intent(in) :: mat(:, :) 350 CMPLX, intent(out) :: zeigenvec(:, :) 351 CMPLX, intent(out) :: zeigenval(:) 352 CMPLX, intent(out) :: zmat(:, :, :) 353 354 integer :: i, alpha, beta 355 CMPLX, allocatable :: knaught(:, :) 356 CMPLX, allocatable :: lambdaminusdm(:, :) 357 CMPLX, allocatable :: ilambdaminusdm(:, :) 358 CMPLX, allocatable :: unit(:, :) 359 360 PUSH_SUB(lalg_zeigenderivatives) 361 362 SAFE_ALLOCATE(unit(1:n, 1:n)) 363 SAFE_ALLOCATE(knaught(1:n, 1:n)) 364 SAFE_ALLOCATE(lambdaminusdm(1:n, 1:n)) 365 SAFE_ALLOCATE(ilambdaminusdm(1:n, 1:n)) 366 367 zeigenvec = mat 368 call lalg_eigensolve_nonh(n, zeigenvec, zeigenval) 369 370 unit = M_z0 371 do i = 1, n 372 unit(i, i) = M_z1 373 end do 374 375 do i = 1, n 376 do alpha = 1, n 377 do beta = 1, n 378 knaught(alpha, beta) = - zeigenvec(alpha, i)*conjg(zeigenvec(beta, i)) 379 end do 380 end do 381 knaught = knaught + unit 382 lambdaminusdm = zeigenval(i)*unit - mat 383 call lalg_zpseudoinverse(n, lambdaminusdm, ilambdaminusdm) 384 zmat(:, :, i) = matmul(ilambdaminusdm, knaught) 385 end do 386 387 SAFE_DEALLOCATE_A(unit) 388 SAFE_DEALLOCATE_A(knaught) 389 SAFE_DEALLOCATE_A(lambdaminusdm) 390 SAFE_DEALLOCATE_A(ilambdaminusdm) 391 POP_SUB(lalg_zeigenderivatives) 392 end subroutine lalg_zeigenderivatives 393 394 395 !>------------------------------------------------- 396 !! Computes the Moore-Penrose pseudoinverse of a 397 !! complex matrix. 398 !!--------------------------------------------- 399 subroutine lalg_zpseudoinverse(n, mat, imat) 400 integer, intent(in) :: n 401 CMPLX, intent(in) :: mat(:, :) 402 CMPLX, intent(out) :: imat(:, :) 403 404 integer :: i 405 CMPLX, allocatable :: u(:, :), vt(:, :), sigma(:, :) 406 FLOAT, allocatable :: sg_values(:) 407 408 PUSH_SUB(lalg_zpseudoinverse) 409 410 SAFE_ALLOCATE(u(1:n, 1:n)) 411 SAFE_ALLOCATE(vt(1:n, 1:n)) 412 SAFE_ALLOCATE(sigma(1:n, 1:n)) 413 SAFE_ALLOCATE(sg_values(1:n)) 414 415 imat = mat 416 call lalg_singular_value_decomp(n, n, imat, u, vt, sg_values) 417 418 sigma = M_z0 419 do i = 1, n 420 if(abs(sg_values(i)) <= CNST(1.0e-12) * maxval(abs(sg_values))) then 421 sigma(i, i) = M_z0 422 else 423 sigma(i, i) = M_z1 / sg_values(i) 424 end if 425 end do 426 427 vt = conjg(transpose(vt)) 428 u = conjg(transpose(u)) 429 imat = matmul(vt, matmul(sigma, u)) 430 431 ! Check if we truly have a pseudoinverse 432 vt = matmul(mat, matmul(imat, mat)) - mat 433 if( maxval(abs(vt)) > CNST(1.0e-10) * maxval(abs(mat))) then 434 write(*, *) maxval(abs(vt)) 435 write(*, *) vt 436 write(*, *) 437 write(*, *) CNST(1.0e-10) * maxval(abs(mat)) 438 write(*, *) maxval(abs(vt)) > CNST(1.0e-10) * maxval(abs(mat)) 439 write(*, *) mat 440 write(message(1), '(a)') 'Pseudoinverse failed.' 441 call messages_fatal(1) 442 end if 443 444 SAFE_DEALLOCATE_A(u) 445 SAFE_DEALLOCATE_A(vt) 446 SAFE_DEALLOCATE_A(sigma) 447 SAFE_DEALLOCATE_A(sg_values) 448 POP_SUB(lalg_zpseudoinverse) 449 end subroutine lalg_zpseudoinverse 450 451 452 !>------------------------------------------------- 453 !! The purpose of this routine is to check that "lalg_zeigenderivatives" 454 !! is working properly, and therefore, it is not really called anywhere 455 !! in the code. It is here only for debugging purposes (perhaps it will 456 !! disappear in the future...) 457 !!------------------------------------------------- 458 subroutine lalg_check_zeigenderivatives(n, mat) 459 integer, intent(in) :: n 460 CMPLX, intent(in) :: mat(:, :) 461 462 integer :: alpha, beta, gamma, delta 463 CMPLX :: deltah, zuder_direct, zder_direct, zuder_directplus, zuder_directminus 464 CMPLX, allocatable :: zeigenvec(:, :), dm(:, :), zeigref_(:, :), zeigenval(:), mmatrix(:, :, :) 465 CMPLX, allocatable :: zeigplus(:), zeigminus(:), zeig0(:), zeigplusminus(:), zeigminusplus(:) 466 467 PUSH_SUB(lalg_check_zeigenderivatives) 468 469 SAFE_ALLOCATE(zeigenvec(1:n, 1:n)) 470 SAFE_ALLOCATE(dm(1:n, 1:n)) 471 SAFE_ALLOCATE(zeigref_(1:n, 1:n)) 472 SAFE_ALLOCATE(zeigenval(1:n)) 473 SAFE_ALLOCATE(mmatrix(1:n, 1:n, 1:n)) 474 SAFE_ALLOCATE(zeigplus(1:n)) 475 SAFE_ALLOCATE(zeigminus(1:n)) 476 SAFE_ALLOCATE(zeig0(1:n)) 477 SAFE_ALLOCATE(zeigplusminus(1:n)) 478 SAFE_ALLOCATE(zeigminusplus(1:n)) 479 480 ASSERT(n == 2) 481 482 dm = mat 483 call lalg_zeigenderivatives(2, dm, zeigref_, zeigenval, mmatrix) 484 485 486 deltah = (CNST(0.000001), CNST(0.000)) 487 !deltah = M_z1 * maxval(abs(dm)) * CNST(0.001) 488 do alpha = 1, n 489 do beta = 1, n 490 zder_direct = lalg_zdni(zeigref_(:, 2), alpha, beta) 491 zuder_direct = lalg_zduialpha(zeigref_(:, 2), mmatrix(:, :, 2), 2, alpha, beta) 492 493 zeigenvec = dm 494 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah 495 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigplus) 496 zeigenvec(:, 1) = zeigenvec(:, 1)/sum(conjg(zeigref_(1:2, 1))*zeigenvec(1:2, 1)) 497 zeigenvec(:, 2) = zeigenvec(:, 2)/sum(conjg(zeigref_(1:2, 2))*zeigenvec(1:2, 2)) 498 zuder_directplus = zeigenvec(2, 2) 499 500 zeigenvec = dm 501 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah 502 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigminus) 503 zeigenvec(:, 1) = zeigenvec(:, 1)/sum(conjg(zeigref_(1:2, 1))*zeigenvec(1:2, 1)) 504 zeigenvec(:, 2) = zeigenvec(:, 2)/sum(conjg(zeigref_(1:2, 2))*zeigenvec(1:2, 2)) 505 zuder_directminus = zeigenvec(2, 2) 506 507 write(*, '(2i1,4f24.12)') alpha, beta, zder_direct, (zeigplus(2) - zeigminus(2))/(M_TWO * deltah) 508 write(*, '(2i1,4f24.12)') alpha, beta, & 509 zuder_direct, (zuder_directplus - zuder_directminus) / (M_TWO * deltah) 510 511 do gamma = 1, n 512 do delta = 1, n 513 if(alpha == gamma .and. beta == delta) then 514 zder_direct = lalg_zd2ni(zeigref_(:, 1), mmatrix(:, :, 1), alpha, beta, gamma, delta) 515 516 zeigenvec = dm 517 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeig0) 518 519 zeigenvec = dm 520 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah 521 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigplus) 522 523 zeigenvec = dm 524 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah 525 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigminus) 526 527 write(*, '(4i1,4f24.12)') alpha, beta, gamma, delta, & 528 zder_direct, (zeigplus(1) + zeigminus(1) - M_TWO*zeig0(1))/(deltah**2) 529 else 530 zder_direct = lalg_zd2ni(zeigref_(:, 1), mmatrix(:, :, 1), alpha, beta, gamma, delta) 531 532 zeigenvec = dm 533 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah 534 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) + deltah 535 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigplus) 536 537 zeigenvec = dm 538 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah 539 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) - deltah 540 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigminus) 541 542 zeigenvec = dm 543 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) + deltah 544 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) - deltah 545 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigplusminus) 546 547 zeigenvec = dm 548 zeigenvec(alpha, beta) = zeigenvec(alpha, beta) - deltah 549 zeigenvec(gamma, delta) = zeigenvec(gamma, delta) + deltah 550 call lalg_eigensolve_nonh(2, zeigenvec(:, :), zeigminusplus) 551 552 write(*, '(4i1,4f24.12)') alpha, beta, gamma, delta, & 553 zder_direct, (zeigplus(1) + zeigminus(1) - zeigplusminus(1) - zeigminusplus(1))/(M_FOUR*deltah**2) 554 555 556 end if 557 end do 558 end do 559 560 end do 561 end do 562 563 SAFE_DEALLOCATE_A(zeigenval) 564 SAFE_DEALLOCATE_A(mmatrix) 565 SAFE_DEALLOCATE_A(zeigenvec) 566 SAFE_DEALLOCATE_A(dm) 567 SAFE_DEALLOCATE_A(zeigref_) 568 SAFE_DEALLOCATE_A(zeigplus) 569 SAFE_DEALLOCATE_A(zeigminus) 570 SAFE_DEALLOCATE_A(zeig0) 571 SAFE_DEALLOCATE_A(zeigplusminus) 572 SAFE_DEALLOCATE_A(zeigminusplus) 573 POP_SUB(lalg_check_zeigenderivatives) 574 end subroutine lalg_check_zeigenderivatives 575 576 CMPLX function lalg_zdni(eigenvec, alpha, beta) 577 integer, intent(in) :: alpha, beta 578 CMPLX :: eigenvec(2) 579 lalg_zdni = conjg(eigenvec(alpha)) * eigenvec(beta) 580 end function lalg_zdni 581 582 CMPLX function lalg_zduialpha(eigenvec, mmatrix, alpha, gamma, delta) 583 integer, intent(in) :: alpha, gamma, delta 584 CMPLX :: eigenvec(2), mmatrix(2, 2) 585 lalg_zduialpha = mmatrix(alpha, gamma) * eigenvec(delta) 586 end function lalg_zduialpha 587 588 CMPLX function lalg_zd2ni(eigenvec, mmatrix, alpha, beta, gamma, delta) 589 integer, intent(in) :: alpha, beta, gamma, delta 590 CMPLX :: eigenvec(2), mmatrix(2, 2) 591 lalg_zd2ni = conjg(mmatrix(alpha, delta) * eigenvec(gamma)) * eigenvec(beta) + & 592 conjg(eigenvec(alpha)) * mmatrix(beta, gamma) * eigenvec(delta) 593 end function lalg_zd2ni 594 595#include "undef.F90" 596#include "complex.F90" 597#include "lalg_adv_lapack_inc.F90" 598 599#include "undef.F90" 600#include "real.F90" 601#include "lalg_adv_lapack_inc.F90" 602 603end module lalg_adv_oct_m 604 605!! Local Variables: 606!! mode: f90 607!! coding: utf-8 608!! End: 609