Lightning

Usage

Development

API

  1. Docs
  2. Program Listing for File UtilLinearAlg.hpp

Program Listing for File UtilLinearAlg.hpp

↰ Return to documentation for file (pennylane_lightning/core/utils/UtilLinearAlg.hpp)

// Copyright 2018-2024 Xanadu Quantum Technologies Inc.

// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at

//     http://www.apache.org/licenses/LICENSE-2.0

// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#pragma once

#include <algorithm>
#include <complex>
#include <string>
#include <vector>

#include "Error.hpp"

#ifdef SCIPY_OPENBLAS_ENABLED
#include "BLASLibLoaderManager.hpp"
#endif

namespace {
// Declare heev function pointers to access corresponding functions in
// OpenBLAS
using zheevPtr = void (*)(const char *, const char *, const int *,
                          std::complex<double> *, const int *, double *,
                          std::complex<double> *, const int *, double *, int *);
using cheevPtr = void (*)(const char *, const char *, const int *,
                          std::complex<float> *, const int *, float *,
                          std::complex<float> *, const int *, float *, int *);
} // namespace

namespace Pennylane::Util {
#ifdef SCIPY_OPENBLAS_ENABLED
template <typename T>
void compute_diagonalizing_gates(int n, int lda,
                                 const std::vector<std::complex<T>> &Ah,
                                 std::vector<T> &eigenVals,
                                 std::vector<std::complex<T>> &unitary) {
    auto &&blasLoader = BLASLibLoaderManager::getInstance();
    auto *blasLibLoader = blasLoader.getBLASLib();

    eigenVals.clear();
    eigenVals.resize(n);
    unitary = std::vector<std::complex<T>>(n * n, {0, 0});

    std::vector<std::complex<T>> ah(n * lda, {0.0, 0.0});

    // TODO optmize transpose
    for (std::size_t i = 0; i < static_cast<std::size_t>(n); i++) {
        for (std::size_t j = 0; j <= i; j++) {
            ah[j * n + i] = Ah[i * lda + j];
        }
    }

    constexpr char jobz =
        'V'; // Enable both eigenvalues and eigenvectors computation
    constexpr char uplo = 'L'; // Upper triangle of matrix is stored
    std::vector<std::complex<T>> work_query(1); // Vector for optimal size query
    int lwork = -1;                             // Optimal workspace size query
    std::vector<T> rwork(3 * n - 2);            // Real workspace array
    int info;

    if constexpr (std::is_same<T, float>::value) {
        auto cheev = blasLibLoader->getSymbol<cheevPtr>("scipy_cheev_");
        // Query optimal workspace size
        cheev(&jobz, &uplo, &n, ah.data(), &lda, eigenVals.data(),
              work_query.data(), &lwork, rwork.data(), &info);
        // Allocate workspace
        lwork = static_cast<int>(work_query[0].real());
        std::vector<std::complex<T>> work_optimal(lwork, {0, 0});
        // Perform eigenvalue and eigenvector computation
        cheev(&jobz, &uplo, &n, ah.data(), &lda, eigenVals.data(),
              work_optimal.data(), &lwork, rwork.data(), &info);
    } else {
        auto zheev = blasLibLoader->getSymbol<zheevPtr>("scipy_zheev_");
        // Query optimal workspace size
        zheev(&jobz, &uplo, &n, ah.data(), &lda, eigenVals.data(),
              work_query.data(), &lwork, rwork.data(), &info);
        // Allocate workspace
        lwork = static_cast<int>(work_query[0].real());
        std::vector<std::complex<T>> work_optimal(lwork, {0, 0});
        // Perform eigenvalue and eigenvector computation
        zheev(&jobz, &uplo, &n, ah.data(), &lda, eigenVals.data(),
              work_optimal.data(), &lwork, rwork.data(), &info);
    }

    std::transform(ah.begin(), ah.end(), unitary.begin(),
                   [](std::complex<T> value) {
                       return std::complex<T>{value.real(), -value.imag()};
                   });
}
#else
template <typename T>
void compute_diagonalizing_gates(
    [[maybe_unused]] int n, [[maybe_unused]] int lda,
    [[maybe_unused]] const std::vector<std::complex<T>> &Ah,
    [[maybe_unused]] std::vector<T> &eigenVals,
    [[maybe_unused]] std::vector<std::complex<T>> &unitary) {
    PL_ABORT("Decompose Hermitian matrix into diagonal matrix and unitaries is "
             "only available"
             "with scipy-openblas. Consider enabling this feature by setting "
             "ENABLE_SCIPY_OPENBLAS.");
}
#endif

} // namespace Pennylane::Util