! ============================================================================= !> Module containing the implementation for the ARPACK general-type solver, that is, !! given the general eigenvalue problem $$ AX = \omega BX, $$ find \(k\) eigenvalues !! that satisfy a given criterion. submodule (mod_solvers:smod_arpack_main) smod_arpack_general use mod_banded_matrix, only: banded_matrix_t use mod_transform_matrix, only: matrix_to_banded use mod_banded_operations, only: multiply implicit none contains !> Implementation of the ARPACK general solver. module procedure solve_arpack_general !> contains the basis vectors complex(dp) :: basis_vectors(arpack_cfg%get_evpdim(), arpack_cfg%get_ncv()) !> work array of length 3N complex(dp) :: workd(3 * arpack_cfg%get_evpdim()) !> work array complex(dp) :: workl(arpack_cfg%get_lworkl()) !> work array of length ncv real(dp) :: rwork(arpack_cfg%get_ncv()) !> integer array with pointers to mark work array locations integer :: ipntr(14) !> logical array of dimension ncv, sets which Ritz vectors to compute logical :: select_vectors(arpack_cfg%get_ncv()) !> work array of length 2*ncv complex(dp) :: workev(2 * arpack_cfg%get_ncv()) integer :: diags logical :: converged type(banded_matrix_t) :: amat_banded, bmat_banded integer :: xstart, xend, ystart, yend complex(dp) :: axvector(arpack_cfg%get_evpdim()) ! we fill at most the full quadblock, so -1 diagonals ! IDEA: maybe do a pass over the B-matrix first to figure out how many diagonals ! we need? This is problem-dependent so in some cases there might be quite some ! room here for optimisations diags = settings%dims%get_dim_quadblock() - 1 call logger%debug("converting A-matrix into banded structure") call matrix_to_banded( & matrix=matrix_A, subdiags=diags, superdiags=diags, banded=amat_banded & ) call logger%debug("converting B-matrix into banded structure") call matrix_to_banded( & matrix=matrix_B, subdiags=diags, superdiags=diags, banded=bmat_banded & ) call logger%debug("doing Arnoldi iteration") converged = .false. ! we keep iterating for as long as the eigenvalues are not converged. ! If convergence is achieved or the maximum number of iterations is reached, ! ARPACK sets ido=99 so the while-loop will always break. do while (.not. converged) call znaupd( & arpack_cfg%ido, & arpack_cfg%get_bmat(), & arpack_cfg%get_evpdim(), & arpack_cfg%get_which(), & arpack_cfg%get_nev(), & arpack_cfg%get_tolerance(), & arpack_cfg%residual, & arpack_cfg%get_ncv(), & basis_vectors, & size(basis_vectors, dim=1), & arpack_cfg%iparam, & ipntr, & workd, & workl, & arpack_cfg%get_lworkl(), & rwork, & arpack_cfg%info & ) ! in the following y is given by workd(ipntr(2)), x by workd(ipntr(1)). ! we need to explicitly set start and end to avoid incompatible ranks xstart = ipntr(1) xend = xstart + arpack_cfg%get_evpdim() - 1 ystart = ipntr(2) yend = ystart + arpack_cfg%get_evpdim() - 1 select case(arpack_cfg%ido) case(-1, 1) ! get y <--- OP*x, we do not calculate OP*x explicitly. ! We need R = OP*x = inv[B]*A*x, so do the following: ! 1. calculate u = A*x ! 2. solve linear system B * R = u for R axvector = multiply(amat_banded, workd(xstart:xend)) workd(ystart:yend) = solve_linear_system_complex_banded( & bandmatrix=bmat_banded, vector=axvector & ) case default ! when convergence is achieved or maxiter is reached exit end select end do call arpack_cfg%parse_znaupd_info(converged) ! if we have a normal exit, extract the eigenvalues through zneupd call zneupd( & .true., & ! always calculate eigenvectors, negligible additional cost in ARPACK "A", & ! calculate Ritz vectors select_vectors, & omega(1:arpack_cfg%get_nev()), & vr(:, 1:arpack_cfg%get_nev()), & size(vr, dim=1), & (0.0d0, 0.0d0), & ! sigma value, not needed here workev, & arpack_cfg%get_bmat(), & arpack_cfg%get_evpdim(), & arpack_cfg%get_which(), & arpack_cfg%get_nev(), & arpack_cfg%get_tolerance(), & arpack_cfg%residual, & arpack_cfg%get_ncv(), & basis_vectors, & size(basis_vectors, dim=1), & arpack_cfg%iparam, & ipntr, & workd, & workl, & arpack_cfg%get_lworkl(), & rwork, & arpack_cfg%info & ) call arpack_cfg%parse_zneupd_info() call arpack_cfg%parse_finished_stats() end procedure solve_arpack_general end submodule smod_arpack_general