diff --git a/quest/src/core/utilities.cpp b/quest/src/core/utilities.cpp
index 999bd9a72..b49838d16 100644
--- a/quest/src/core/utilities.cpp
+++ b/quest/src/core/utilities.cpp
@@ -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"
 
@@ -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++) {
 
@@ -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
@@ -530,8 +527,15 @@ 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);
@@ -539,8 +543,15 @@ bool util_isUnitary(CompMatr m, qreal eps) {
 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);
@@ -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
@@ -609,8 +616,15 @@ 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);
@@ -618,8 +632,15 @@ bool util_isHermitian(CompMatr m, qreal eps) {
 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);
@@ -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
@@ -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++) {
 
@@ -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);
 }
 
diff --git a/quest/src/core/utilities.hpp b/quest/src/core/utilities.hpp
index 4b7fb5db6..1667a8265 100644
--- a/quest/src/core/utilities.hpp
+++ b/quest/src/core/utilities.hpp
@@ -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
  */
diff --git a/quest/src/gpu/gpu_subroutines.cpp b/quest/src/gpu/gpu_subroutines.cpp
index 5e18048f7..6c993ef45 100644
--- a/quest/src/gpu/gpu_subroutines.cpp
+++ b/quest/src/gpu/gpu_subroutines.cpp
@@ -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
+}
diff --git a/quest/src/gpu/gpu_subroutines.hpp b/quest/src/gpu/gpu_subroutines.hpp
index ff42c2239..8176a8d63 100644
--- a/quest/src/gpu/gpu_subroutines.hpp
+++ b/quest/src/gpu/gpu_subroutines.hpp
@@ -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
\ No newline at end of file
diff --git a/quest/src/gpu/gpu_thrust.cuh b/quest/src/gpu/gpu_thrust.cuh
index 9f8d8f1ab..a8404846b 100644
--- a/quest/src/gpu/gpu_thrust.cuh
+++ b/quest/src/gpu/gpu_thrust.cuh
@@ -59,6 +59,8 @@
 #include <thrust/iterator/transform_iterator.h>
 #include <thrust/system/system_error.h>
 
+#include <algorithm>
+
 
 
 /*
@@ -641,6 +643,123 @@ struct functor_setRandomStateVecAmp : public thrust::unary_function<qindex,cu_qc
 };
 
 
+struct functor_compmatr_isUnitaryTerm : public thrust::unary_function<qindex,bool>{
+
+    cu_qcomp* matr;
+    qreal eps;
+    qindex dim;
+
+    functor_compmatr_isUnitaryTerm(cu_qcomp* matr, qreal eps, qindex dim) :
+        matr(matr), eps(eps), dim(dim)
+    {}
+
+    __host__ __device__ bool operator()(qindex i) {
+
+        qindex r = i / dim;
+        qindex c = i % dim;
+
+        cu_qcomp elem = getCuQcomp(0, 0);
+        for (qindex k=0; k<dim; k++)
+            elem = elem + (matr[r * dim + k] * getCompConj(matr[c * dim + k]));
+
+        cu_qcomp target = getCuQcomp(r == c, 0);
+        return getCompNorm(elem - target) <= eps;
+    }
+};
+
+struct functor_diagmatr_isUnitaryTerm : public thrust::unary_function<qindex,bool>{
+
+    cu_qcomp* diags;
+    qreal eps;
+
+    functor_diagmatr_isUnitaryTerm(cu_qcomp* diags, qreal eps) :
+        diags(diags), eps(eps)
+    {}
+
+    __host__ __device__ bool operator()(qindex i) {
+
+        // We want |sqrt(norm) - 1| <= eps;
+        // equivalent to (1-eps)^2 <= norm <= (1+eps)^2
+        // eps^2 is small, so approximate bound as
+        // [1-2eps, 1+2eps], or |norm - 1| <= 2*eps
+
+        qreal norm = getCompNorm(diags[i]);
+        return fabs(norm - 1) <= 2 * eps;
+    }
+};
+
+struct functor_compmatr_isHermitianTerm : public thrust::unary_function<qindex,bool>{
+
+    // check adjoint(elems) == elems
+
+    cu_qcomp* matr;
+    qreal eps;
+    qindex dim;
+
+    functor_compmatr_isHermitianTerm(cu_qcomp* matr, qreal eps, qindex dim) :
+        matr(matr), eps(eps), dim(dim)
+    {}
+
+    __host__ __device__ bool operator()(qindex i) {
+
+        qindex row = i / dim;
+        qindex col = i % dim;
+
+        
+        if (col >= row)
+            return true;
+
+        cu_qcomp elem = matr[row * dim + col];
+        cu_qcomp conjOfMirror = getCompConj(matr[col * dim + row]);
+
+        return getCompNorm(elem - conjOfMirror) <= eps;
+    }
+};
+
+struct functor_diagmatr_isHermitianTerm : public thrust::unary_function<qindex,bool>{
+
+    cu_qcomp* diags;
+    qreal eps;
+
+    functor_diagmatr_isHermitianTerm(cu_qcomp* diags, qreal eps) :
+       diags(diags), eps(eps)
+    {}
+
+    __host__ __device__ bool operator()(qindex i) {
+        qreal imag = diags[i].y;
+        return fabs(imag) <= eps;
+    }
+};
+
+struct functor_krausmap_isCPTPTerm : public thrust::unary_function<qindex,bool>{
+    
+    cu_qcomp* matr;
+    qreal eps;
+    qindex dim;
+    qindex numMatrices;
+
+    functor_krausmap_isCPTPTerm(cu_qcomp* matr, qreal eps, qindex dim, qindex numMatrices) :
+        matr(matr), eps(eps), dim(dim), numMatrices(numMatrices)
+    {}
+
+    __host__ __device__ bool operator()(qindex i) {
+
+        qindex row = i / dim;
+        qindex col = i % dim;
+
+        cu_qcomp elem = getCuQcomp(0, 0);
+        for (qindex n=0; n<numMatrices; n++)
+            for (qindex k=0; k<dim; k++) {
+                qindex base = n * dim * dim + k * dim;
+                elem = elem + (getCompConj(matr[base + row]) * matr[base + col]);
+            }
+
+        cu_qcomp target = getCuQcomp(row == col, 0);
+        return getCompNorm(elem - target) <= eps;
+    }
+};
+
+
 
 /*
  * MATRIX INITIALISATION
@@ -1083,4 +1202,104 @@ void thrust_statevec_initUnnormalisedUniformlyRandomPureStateAmps_sub(Qureg qure
 
 
 
+/*
+ * MATRIX PROPERTIES
+ */
+
+
+bool thrust_compmatr_isUnitary_sub(CompMatr matr, qreal eps) {
+
+    qindex dim = matr.numRows;
+
+    //functor accepts an index and returns a boolean
+    auto functor = functor_compmatr_isUnitaryTerm(toCuQcomps(matr.gpuElemsFlat), eps, dim);
+
+    auto indIter = thrust::make_counting_iterator(0);
+    qindex numIts = dim * dim;
+
+    return thrust::transform_reduce(
+        indIter, indIter + numIts,
+        functor, true, thrust::logical_and<bool>()
+    );
+}
+
+bool thrust_diagmatr_isUnitary_sub(DiagMatr matr, qreal eps){
+
+    qindex dim = matr.numElems;
+
+    //functor accepts an index and returns a boolean
+    auto functor = functor_diagmatr_isUnitaryTerm(toCuQcomps(matr.gpuElems), eps);
+
+    auto indIter = thrust::make_counting_iterator(0);
+    qindex numIts = dim;
+
+    return thrust::transform_reduce(
+        indIter, indIter + numIts,
+        functor, true, thrust::logical_and<bool>()
+    );
+}
+
+bool thrust_compmatr_isHermitian_sub(CompMatr matr, qreal eps){
+
+    qindex dim = matr.numRows;
+
+    //functor accepts an index and returns a boolean
+    auto functor = functor_compmatr_isHermitianTerm(toCuQcomps(matr.gpuElemsFlat), eps, dim);
+
+    auto indIter = thrust::make_counting_iterator(0);
+    qindex numIts = dim * dim;
+
+    return thrust::transform_reduce(
+        indIter, indIter + numIts,
+        functor, true, thrust::logical_and<bool>()
+    );
+}
+
+bool thrust_diagmatr_isHermitian_sub(DiagMatr matr, qreal eps){
+    
+    qindex dim = matr.numElems;
+
+    //functor accepts an index and returns a boolean
+    auto functor = functor_diagmatr_isHermitianTerm(toCuQcomps(matr.gpuElems), eps);
+
+    auto indIter = thrust::make_counting_iterator(0);
+    qindex numIts = dim;
+
+    return thrust::transform_reduce(
+        indIter, indIter + numIts,
+        functor, true, thrust::logical_and<bool>()
+    );
+}
+
+bool thrust_krausmap_isCPTP_sub(KrausMap map, qreal eps){
+
+    qindex dim = map.numRows;
+    int numMatrices = map.numMatrices;
+
+    //map.matrices is CPU-only
+    //flatten to copy over to GPU, but matrices and rows are non-contigious
+    vector<qcomp> hostMatrices(numMatrices * dim * dim);
+    for (int n=0; n<numMatrices; n++){
+        for (qindex k=0; k<dim; k++){
+            std::copy(map.matrices[n][k], map.matrices[n][k] + dim, hostMatrices.begin() + (n*dim*dim + k*dim));
+        }
+    }
+    
+    thrust::device_vector<qcomp> devMatrices(hostMatrices);
+    cu_qcomp* devMatricesPtr = toCuQcomps(thrust::raw_pointer_cast(devMatrices.data()));
+
+    //functor accepts an index and returns a boolean
+    auto functor = functor_krausmap_isCPTPTerm(devMatricesPtr, eps, dim, numMatrices);
+
+    auto indIter = thrust::make_counting_iterator(0);
+    qindex numIts = dim * dim;
+
+    return thrust::transform_reduce(
+        indIter, indIter + numIts,
+        functor, true, thrust::logical_and<bool>()
+    );
+}
+
+
+
 #endif // GPU_THRUST_HPP
\ No newline at end of file