! ============================================================================= !> Submodule containing the implementation of the QR-invert algorithm. !! The original problem is written as a standard eigenvalue problem through !! $$ \mathcal{B}^{-1}\mathcal{A}\textbf{X} = \omega\textbf{X}\ $$. This is !! done using a LU decomposition via LAPACKS's <tt>zgbsv</tt>. !! Eventually a call to LAPACK's <tt>zgeev</tt> routine is done to obtain !! all eigenvalues and eigenvectors. submodule (mod_solvers) smod_qr_invert use mod_banded_matrix, only: banded_matrix_t, new_banded_matrix use mod_transform_matrix, only: matrix_to_banded implicit none contains !> Solves the eigenvalue problem by rewriting it to a standard form !! through inversion of the B-matrix. !! @warning Throws an error if <tt>matrix_A</tt> or <tt>matrix_B</tt> !! is not a square matrix. @endwarning module procedure qr_invert !> full array containing the \(B^{-1}A\)-matrix complex(dp), allocatable :: array_B_invA(:, :) !> pivoting array integer :: ipiv(matrix_B%matrix_dim) !> banded B-matrix type(banded_matrix_t) :: B_band !> order of matrices integer :: N !> calculate right eigenvectors if "V", omit if "N" character :: jobvr !> dimension of work array integer :: lwork !> work array complex(dp), allocatable :: work(:) !> info parameter, 0 on successful exit integer :: info !> second work array real(dp), allocatable :: rwork(:) !> dummy for left eigenvectors, jobvl = "N" so this is never referenced complex(dp) :: vl(2, 2) !> number of superdiagonals integer :: ku !> number of subdiagonals integer :: kl integer :: irow, icol call matrix_B%get_nb_diagonals(ku=ku, kl=kl) call logger%debug( & "B has " // str(ku) // " superdiagonals and " // str(kl) // " subdiagonals" & ) N = matrix_B%matrix_dim call logger%debug("converting B to banded form") ! for zgbsv later on we need double the amount of subdiagonals B_band = new_banded_matrix(rows=N, cols=N, subdiags=2*kl, superdiags=ku) ! fill banded matrix, see documentation: "The j-th column of B is stored in the ! j-th column of the array AB as follows: ! AB(KL+KU+1+i-j,j) = B(i,j) for max(1,j-KU)<=i<=min(N,j+KL)" do icol = 1, N do irow = max(1, icol - ku), min(N, icol + kl) B_band%AB(kl+ku+1+irow-icol, icol) = matrix_B%get_complex_element(irow, icol) end do end do allocate(array_B_invA(matrix_A%matrix_dim, matrix_A%matrix_dim)) call logger%debug("converting A to dense form") call matrix_to_array(matrix=matrix_A, array=array_B_invA) ! zgbsv calculates X from BX = A, where B is a banded matrix. ! solving this system yields X = B^{-1}A without explicitly inverting call logger%debug("calling LAPACK's zgbsv") call zgbsv( & B_band%m, & kl, & ku, & matrix_A%matrix_dim, & B_band%AB, & size(B_band%AB, dim=1), & ipiv, & array_B_invA, & size(array_B_invA, dim=1), & info & ) if (info /= 0) then call logger%error("LAPACK's zgbsv failed with info = " // str(info)) return end if call B_band%destroy() ! calculate eigenvectors, we don't use the left ones jobvr = "N" if (settings%io%should_compute_eigenvectors()) jobvr = "V" ! allocate rwork array allocate(rwork(2 * N)) ! get lwork call logger%debug("calculating optimal length of work array") allocate(work(1)) call zgeev( & "N", & jobvr, & N, & array_B_invA, & N, & omega, & vl, & size(vl, dim=1), & vr, & size(vr, dim=1), & work, & -1, & rwork, & info & ) lwork = int(work(1)) deallocate(work) ! allocate work array call logger%debug("allocating work array with N = " // str(lwork)) allocate(work(lwork)) ! solve eigenvalue problem call logger%debug("solving evp using QR algorithm zgeev (LAPACK)") call zgeev( & "N", & jobvr, & N, & array_B_invA, & N, & omega, & vl, & size(vl, dim=1), & vr, & size(vr, dim=1), & work, & lwork, & rwork, & info & ) if (info /= 0) then ! LCOV_EXCL_START call logger%warning("LAPACK routine zgeev failed!") call logger%warning("value for the info parameter: " // str(info)) end if ! LCOV_EXCL_STOP ! tear down work arrays deallocate(array_B_invA) deallocate(work) deallocate(rwork) end procedure qr_invert end submodule smod_qr_invert