eigh.f90 Source File

submodule(specialmatrices_poisson2D) poisson2D_eigh
   use stdlib_constants, only: pi => pi_dp
   use stdlib_sorting, only: sort, sort_index
   use stdlib_linalg, only: outer_product
   implicit none(type, external)
contains
   module procedure eigvalsh_rdp
   integer(ilp) :: i, j
   real(dp), pointer :: lambda_mat(:, :)

   !> Eigenvalues for Dx and Dy.
   real(dp), allocatable :: lambda_x(:), lambda_y(:)
   lambda_x = -eigvalsh(Strang(A%nx))/A%dx**2
   lambda_y = -eigvalsh(Strang(A%ny))/A%dy**2

   !> Eigenvalues of the 2D Poisson operator.
   allocate (lambda(A%nx*A%ny)); lambda_mat(1:A%nx, 1:A%ny) => lambda
   do concurrent(i=1:A%nx, j=1:A%ny)
      lambda_mat(i, j) = lambda_x(i) + lambda_y(j)
   end do

   !> Sort eigenvalues.
   call sort(lambda)
   end procedure eigvalsh_rdp

   module procedure eigh_rdp
   integer(ilp) :: i, j, index(A%nx*A%ny)
   integer(ilp) :: nx, ny, n, counter
   real(dp), pointer :: lambda_mat(:, :)
   real(dp), allocatable :: lambda_x(:), lambda_y(:)
   real(dp), allocatable :: vecs_x(:, :), vecs_y(:, :)

   nx = A%nx; ny = A%ny; n = nx*ny

   !> Eigenvalues and eigevectors for Dx and Dy.
   call eigh(Strang(A%nx), lambda_x, vecs_x)
   call eigh(Strang(A%ny), lambda_y, vecs_y)
   !> Scale eigenvalues.
   lambda_x = -lambda_x/A%dx**2
   lambda_y = -lambda_y/A%dy**2

   !> Eigenvalues of the 2D Poisson operator.
   allocate (lambda(n)); lambda_mat(1:nx, 1:ny) => lambda
   do concurrent(i=1:nx, j=1:ny)
      lambda_mat(i, j) = lambda_x(i) + lambda_y(j)
   end do

   !> Eigenvectors of the 2D Poisson operator.
   allocate (vectors(n, n)); vectors = 0.0_dp; counter = 1
   do j = 1, ny
      do i = 1, nx
         vectors(:, counter) = pack(outer_product(vecs_x(:, i), vecs_y(:, j)), .true.)
         counter = counter + 1
      end do
   end do

   !> Sort eigenvalues and eigenvectors.
   call sort_index(lambda, index); vectors = vectors(:, index)
   end procedure eigh_rdp
end submodule poisson2D_eigh