Skip to content

Omp/CUDA Utility Accelerations #813

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
jwilso39-uc wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Omp/CUDA Utility Accelerations #813

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
109 changes: 62 additions & 47 deletions quest/src/core/utilities.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@
#include "quest/src/core/utilities.hpp"
#include "quest/src/core/validation.hpp"
#include "quest/src/cpu/cpu_config.hpp"
#include "quest/src/gpu/gpu_subroutines.hpp"
#include "quest/src/comm/comm_config.hpp"
#include "quest/src/comm/comm_routines.hpp"

Expand Down Expand Up @@ -482,11 +483,9 @@ template <typename T>
bool getUnitarity(T elems, qindex dim, qreal eps) {
assert_utilsGivenNonZeroEpsilon(eps);

/// @todo
/// consider multithreading or GPU-accelerating this
/// when caller is big and e.g. has GPU memory

// check m * dagger(m) == identity
bool unitary = true;
#pragma omp parallel for reduction(&&:unitary) if(dim >= MIN_DIM_FOR_UTIL_DENSE_UNITARITY_MULTITHREADING && getQuESTEnv().isMultithreaded)
for (qindex r=0; r<dim; r++) {
for (qindex c=0; c<dim; c++) {

Expand All @@ -497,27 +496,25 @@ bool getUnitarity(T elems, qindex dim, qreal eps) {

// check if further than epsilon from identity[r,c]
if (!isApprox(elem, qcomp(r==c,0), eps))
return false;
unitary = false;
}
}

return true;
return unitary;
}

// diagonal version doesn't need templating because array decays to pointer, yay!
bool getUnitarity(qcomp* diags, qindex dim, qreal eps) {
assert_utilsGivenNonZeroEpsilon(eps);

/// @todo
/// consider multithreading or GPU-accelerating this
/// when caller is big and e.g. has GPU memory

// check every element has unit magnitude
bool unitary = true;
#pragma omp parallel for reduction(&&:unitary) if(dim >= MIN_DIM_FOR_UTIL_DIAG_UNITARITY_MULTITHREADING && getQuESTEnv().isMultithreaded)
for (qindex i=0; i<dim; i++)
if (!isApprox(std::abs(diags[i]), 1, eps))
return false;
unitary = false;

return true;
return unitary;
}

// unitarity of fixed-size matrices is always computed afresh
Expand All @@ -530,17 +527,31 @@ bool util_isUnitary(DiagMatr2 m, qreal eps) { return getUnitarity(m.elems, m.num
bool util_isUnitary(CompMatr m, qreal eps) {

// compute and record unitarity if not already known
if (*(m.isApproxUnitary) == validate_STRUCT_PROPERTY_UNKNOWN_FLAG)
*(m.isApproxUnitary) = getUnitarity(m.cpuElems, m.numRows, eps);
if (*(m.isApproxUnitary) == validate_STRUCT_PROPERTY_UNKNOWN_FLAG) {

if (util_isGpuAcceleratedMatrix(m) && m.numRows >= MIN_DIM_FOR_UTIL_DENSE_UNITARITY_GPU) {
*(m.isApproxUnitary) = gpu_compmatr_isUnitary_sub(m, eps);
}
else {
*(m.isApproxUnitary) = getUnitarity(m.cpuElems, m.numRows, eps);
}
}

// eps may have been ignored
return *(m.isApproxUnitary);
}
bool util_isUnitary(DiagMatr m, qreal eps) {

// compute and record unitarity if not already known
if (*(m.isApproxUnitary) == validate_STRUCT_PROPERTY_UNKNOWN_FLAG)
*(m.isApproxUnitary) = getUnitarity(m.cpuElems, m.numElems, eps);
if (*(m.isApproxUnitary) == validate_STRUCT_PROPERTY_UNKNOWN_FLAG) {

if (util_isGpuAcceleratedMatrix(m) && m.numElems >= MIN_DIM_FOR_UTIL_DIAG_UNITARITY_GPU) {
*(m.isApproxUnitary) = gpu_diagmatr_isUnitary_sub(m, eps);
}
else {
*(m.isApproxUnitary) = getUnitarity(m.cpuElems, m.numElems, eps);
}
}

// eps may have been ignored
return *(m.isApproxUnitary);
Expand Down Expand Up @@ -570,33 +581,29 @@ template <typename T>
bool getHermiticity(T elems, qindex dim, qreal eps) {
assert_utilsGivenNonZeroEpsilon(eps);

/// @todo
/// consider multithreading or GPU-accelerating this
/// when caller is big and e.g. has GPU memory

// check adjoint(elems) == elems
bool hermitian = true;
#pragma omp parallel for reduction(&&:hermitian) if(dim >= MIN_DIM_FOR_UTIL_DENSE_HERMITICITY_MULTITHREADING && getQuESTEnv().isMultithreaded)
for (qindex r=0; r<dim; r++)
for (qindex c=0; c<r; c++)
if (!isApprox(elems[r][c], std::conj(elems[c][r]), eps))
return false;
hermitian = false;

return true;
return hermitian;
}

// diagonal version doesn't need templating because array decays to pointer, yay!
bool getHermiticity(qcomp* diags, qindex dim, qreal eps) {
assert_utilsGivenNonZeroEpsilon(eps);

/// @todo
/// consider multithreading or GPU-accelerating this
/// when caller is big and e.g. has GPU memory

// check every element has a zero (or <eps) imaginary component
bool hermitian = true;
#pragma omp parallel for reduction(&&:hermitian) if(dim >= MIN_DIM_FOR_UTIL_DIAG_HERMITICITY_MULTITHREADING && getQuESTEnv().isMultithreaded)
for (qindex i=0; i<dim; i++)
if (!isApprox(std::imag(diags[i]), 0, eps))
return false;
hermitian = false;

return true;
return hermitian;
}

// hermiticity of fixed-size matrices is always computed afresh
Expand All @@ -609,17 +616,31 @@ bool util_isHermitian(DiagMatr2 m, qreal eps) { return getHermiticity(m.elems, m
bool util_isHermitian(CompMatr m, qreal eps) {

// compute and record hermiticity if not already known
if (*(m.isApproxHermitian) == validate_STRUCT_PROPERTY_UNKNOWN_FLAG)
*(m.isApproxHermitian) = getHermiticity(m.cpuElems, m.numRows, eps);
if (*(m.isApproxHermitian) == validate_STRUCT_PROPERTY_UNKNOWN_FLAG) {

if (util_isGpuAcceleratedMatrix(m) && m.numRows >= MIN_DIM_FOR_UTIL_DENSE_HERMITICITY_GPU) {
*(m.isApproxHermitian) = gpu_compmatr_isHermitian_sub(m, eps);
}
else {
*(m.isApproxHermitian) = getHermiticity(m.cpuElems, m.numRows, eps);
}
}

// eps may have been ignored
return *(m.isApproxHermitian);
}
bool util_isHermitian(DiagMatr m, qreal eps) {

// compute and record hermiticity if not already known
if (*(m.isApproxHermitian) == validate_STRUCT_PROPERTY_UNKNOWN_FLAG)
*(m.isApproxHermitian) = getHermiticity(m.cpuElems, m.numElems, eps);
if (*(m.isApproxHermitian) == validate_STRUCT_PROPERTY_UNKNOWN_FLAG) {

if (util_isGpuAcceleratedMatrix(m) && m.numElems >= MIN_DIM_FOR_UTIL_DIAG_HERMITICITY_GPU) {
*(m.isApproxHermitian) = gpu_diagmatr_isHermitian_sub(m, eps);
}
else {
*(m.isApproxHermitian) = getHermiticity(m.cpuElems, m.numElems, eps);
}
}

// eps may have been ignored
return *(m.isApproxHermitian);
Expand Down Expand Up @@ -647,10 +668,6 @@ bool util_isHermitian(FullStateDiagMatr m, qreal eps) {
bool getWhetherNonZero(qcomp* diags, qindex dim, qreal eps) {
assert_utilsGivenNonZeroEpsilon(eps);

/// @todo
/// consider multithreading or GPU-accelerating this
/// when caller is big and e.g. has GPU memory

for (qindex i=0; i<dim; i++) {

// check each complex element has non-zero abs
Expand Down Expand Up @@ -789,15 +806,14 @@ bool util_isCPTP(KrausMap map, qreal eps) {
if (*(map.isApproxCPTP) != validate_STRUCT_PROPERTY_UNKNOWN_FLAG)
return *(map.isApproxCPTP);

/// @todo
/// if KrausMap is GPU-accelerated, we should maybe
/// instead perform this calculation using the GPU.
/// otherwise, if matrix is large, we should potentially
/// use a multithreaded routine

*(map.isApproxCPTP) = 1;
if (util_isGpuAcceleratedMatrix(map) && map.numRows >= MIN_DIM_FOR_UTIL_CPTP_GPU) {
*(map.isApproxCPTP) = gpu_krausmap_isCPTP_sub(map, eps);
return *(map.isApproxCPTP);
}

// check whether each element satisfies Identity = sum dagger(m)*m
bool cptp = true;
#pragma omp parallel for reduction(&&:cptp) if(map.numRows >= MIN_DIM_FOR_UTIL_CPTP_MULTITHREADING && getQuESTEnv().isMultithreaded)
for (qindex r=0; r<map.numRows; r++) {
for (qindex c=0; c<map.numRows; c++) {

Expand All @@ -811,14 +827,13 @@ bool util_isCPTP(KrausMap map, qreal eps) {
qreal distSquared = std::norm(elem - (r==c));
if (distSquared > eps) {

// by recording the result and returning immediately
*(map.isApproxCPTP) = 0;
return *(map.isApproxCPTP);
// by recording the result (not returning immediately with many threads)
cptp = false;
}
}
}

// always true by this point
*(map.isApproxCPTP) = cptp;
return *(map.isApproxCPTP);
}

Expand Down
25 changes: 25 additions & 0 deletions quest/src/core/utilities.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -270,6 +270,31 @@ CompMatr2 util_getTranspose(CompMatr2 matrix);



/*
* OPENMP MULTITHREADING THRESHOLDS
*/

#define MIN_DIM_FOR_UTIL_DENSE_UNITARITY_MULTITHREADING 16 // In testing, only multithreading w/ >= 4 qubits gives performance improvement
#define MIN_DIM_FOR_UTIL_DIAG_UNITARITY_MULTITHREADING 1024 // >= 10 qubits

#define MIN_DIM_FOR_UTIL_DENSE_HERMITICITY_MULTITHREADING 128 // >= 7 qubits
#define MIN_DIM_FOR_UTIL_DIAG_HERMITICITY_MULTITHREADING 2048 // >= 11 qubits

#define MIN_DIM_FOR_UTIL_CPTP_MULTITHREADING 16 // >= 4 qubits

/*
* GPU MULTITHREADING THRESHOLDS
*/

#define MIN_DIM_FOR_UTIL_DENSE_UNITARITY_GPU 128 // In testing, only gpu acceleration w/ >= 7 qubits gives performance improvement
#define MIN_DIM_FOR_UTIL_DIAG_UNITARITY_GPU 4096 // >= 12 qubits

#define MIN_DIM_FOR_UTIL_DENSE_HERMITICITY_GPU 1024 // >= 10 qubits
#define MIN_DIM_FOR_UTIL_DIAG_HERMITICITY_GPU 16384 // >= 14 qubits

#define MIN_DIM_FOR_UTIL_CPTP_GPU 128 // >= 7 qubits


/*
* MATRIX PROPERTIES
*/
Expand Down
67 changes: 67 additions & 0 deletions quest/src/gpu/gpu_subroutines.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1895,3 +1895,70 @@ void gpu_statevec_initUnnormalisedUniformlyRandomPureStateAmps_sub(Qureg qureg)
error_gpuSimButGpuNotCompiled();
#endif
}



/*
* MATRIX PROPERTIES
*/


bool gpu_compmatr_isUnitary_sub(CompMatr matr, qreal eps) {

#if COMPILE_CUDA || COMPILE_CUQUANTUM

return thrust_compmatr_isUnitary_sub(matr, eps);

#else
error_gpuSimButGpuNotCompiled();
return false;
#endif
}

bool gpu_diagmatr_isUnitary_sub(DiagMatr matr, qreal eps){

#if COMPILE_CUDA || COMPILE_CUQUANTUM

return thrust_diagmatr_isUnitary_sub(matr, eps);

#else
error_gpuSimButGpuNotCompiled();
return false;
#endif
}

bool gpu_compmatr_isHermitian_sub(CompMatr matr, qreal eps){

#if COMPILE_CUDA || COMPILE_CUQUANTUM

return thrust_compmatr_isHermitian_sub(matr, eps);

#else
error_gpuSimButGpuNotCompiled();
return false;
#endif
}

bool gpu_diagmatr_isHermitian_sub(DiagMatr matr, qreal eps){

#if COMPILE_CUDA || COMPILE_CUQUANTUM

return thrust_diagmatr_isHermitian_sub(matr, eps);

#else
error_gpuSimButGpuNotCompiled();
return false;
#endif
}

bool gpu_krausmap_isCPTP_sub(KrausMap map, qreal eps){

#if COMPILE_CUDA || COMPILE_CUQUANTUM

return thrust_krausmap_isCPTP_sub(map, eps);

#else
error_gpuSimButGpuNotCompiled();
return false;
#endif
}
15 changes: 15 additions & 0 deletions quest/src/gpu/gpu_subroutines.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -196,4 +196,19 @@ void gpu_statevec_initDebugState_sub(Qureg qureg);
void gpu_statevec_initUnnormalisedUniformlyRandomPureStateAmps_sub(Qureg qureg);


/*
* MATRIX PROPERTIES
*/

bool gpu_compmatr_isUnitary_sub(CompMatr matr, qreal eps);

bool gpu_diagmatr_isUnitary_sub(DiagMatr matr, qreal eps);

bool gpu_compmatr_isHermitian_sub(CompMatr matr, qreal eps);

bool gpu_diagmatr_isHermitian_sub(DiagMatr matr, qreal eps);

bool gpu_krausmap_isCPTP_sub(KrausMap map, qreal eps);


#endif // GPU_SUBROUTINES_HPP
Loading