| 825 | // T can be qcomp*** or vector<vector<vector<qcomp>>> |
| 826 | template <typename T> |
| 827 | void setSuperoperator(qcomp** superop, T matrices, int numMatrices, qindex logMatrixDim) { |
| 828 | |
| 829 | /// @todo |
| 830 | /// we initialise the superoperator completely serially, under the assumption that the |
| 831 | /// superoperator will be small in size and initialised infrequently. Still, it would |
| 832 | /// be better to provide backend initialisation functions (OpenMP and CUDA accelerated), |
| 833 | /// called when the superoperator size is above some threshold! |
| 834 | |
| 835 | qindex matrixDim = powerOf2(logMatrixDim); |
| 836 | qindex superopDim = matrixDim * matrixDim; |
| 837 | |
| 838 | // clear superoperator |
| 839 | for (qindex r=0; r<superopDim; r++) |
| 840 | for (qindex c=0; c<superopDim; c++) |
| 841 | superop[r][c] = 0; |
| 842 | |
| 843 | // add each matrix's contribution to the superoperator |
| 844 | for (int n=0; n<numMatrices; n++) { |
| 845 | auto matrix = matrices[n]; |
| 846 | |
| 847 | // superop += conj(matrix) (tensor) matrix |
| 848 | for (qindex i=0; i<matrixDim; i++) |
| 849 | for (qindex j=0; j<matrixDim; j++) |
| 850 | for (qindex k=0; k<matrixDim; k++) |
| 851 | for (qindex l=0; l<matrixDim; l++) { |
| 852 | qindex r = i*matrixDim + k; |
| 853 | qindex c = j*matrixDim + l; |
| 854 | superop[r][c] += std::conj(matrix[i][j]) * matrix[k][l]; |
| 855 | } |
| 856 | } |
| 857 | } |
| 858 | void util_setSuperoperator(qcomp** superop, vector<vector<vector<qcomp>>> matrices, int numMatrices, int numQubits) { |
| 859 | setSuperoperator(superop, matrices, numMatrices, numQubits); |
| 860 | } |
no test coverage detected