name_modes Subroutine

subroutine name_modes(X, lambda, phi, names, match, tol)
  !! Replace the eigenvectors of each group of equal eigenvalues (within tol*max|lambda|,
  !! default 1e-8) by the projections of the best-matching catalog fields onto that
  !! eigenspace. match is the cosine between field and eigenspace.
  !
  ! The projections are exact eigenvectors, but not orthogonal within a group. Modes
  ! without a match keep an orthonormal remainder of the eigenspace, named "-".
  real(dp), intent(in) :: X(:,:), lambda(:)
  real(dp), intent(inout) :: phi(:,:)
  character(len=*), intent(out) :: names(:)
  real(dp), intent(out) :: match(:)
  real(dp), intent(in), optional :: tol
  real(dp), parameter :: min_match = 0.5_dp, min_independent = 0.1_dp
  real(dp), allocatable :: F(:,:), V(:,:), P(:,:), Q(:,:), frac(:), c(:)
  character(len=40), allocatable :: tnames(:)
  logical, allocatable :: avail(:)
  real(dp) :: t
  integer :: i0, i1, k, n, j

  t = 1e-8_dp
  if (present(tol)) t = tol
  call catalog(X, F, tnames)
  names = '-'
  match = 0
  i0 = 1
  do while (i0 <= size(lambda))
    i1 = i0
    do while (i1 < size(lambda))
      if (abs(lambda(i1+1) - lambda(i0)) > t*maxval(abs(lambda))) exit
      i1 = i1 + 1
    end do
    k = i1 - i0 + 1
    V = phi(:, i0:i1)
    P = matmul(transpose(V), F)   ! fields in eigenspace coordinates
    frac = norm2(P, 1)
    avail = frac > min_match
    allocate (Q(k, 0))            ! orthonormal basis of the accepted fields
    n = 0
    ! greedy: best match first, catalog order breaks ties
    do while (n < k .and. any(avail))
      j = findloc(avail .and. frac >= maxval(frac, mask=avail) - 1e-8_dp, .true., dim=1)
      avail(j) = .false.
      c = P(:,j)/frac(j)
      c = c - matmul(Q, matmul(transpose(Q), c))
      if (norm2(c) < min_independent) cycle   ! dependent on accepted fields
      Q = reshape([Q, c/norm2(c)], [k, n+1])
      n = n + 1
      phi(:, i0+n-1) = matmul(V, P(:,j))/frac(j)
      names(i0+n-1) = tnames(j)
      match(i0+n-1) = frac(j)
    end do
    ! fill the rest with the orthonormal complement of the accepted fields
    do j = 1, k
      if (n == k) exit
      c = spread(0.0_dp, 1, k)
      c(j) = 1
      c = c - matmul(Q, matmul(transpose(Q), c))
      if (norm2(c) < min_independent) cycle
      Q = reshape([Q, c/norm2(c)], [k, n+1])
      n = n + 1
      phi(:, i0+n-1) = matmul(V, Q(:,n))
    end do
    deallocate (Q)
    i0 = i1 + 1
  end do
end subroutine name_modes