From c04bb7bfacbbcb3a312004d3a69ec22f44e04a66 Mon Sep 17 00:00:00 2001 From: Thomas Loke Date: Mon, 20 Aug 2018 22:22:57 +0800 Subject: [PATCH] Basic C++ implementation of C99 hafnian source code Only implemented for complex (lhafnian.c) support -- real (rhafnian.c) support not added --- .gitignore | 9 +- src/hafnian-cpp/hafnian-cpp/lhafnian.cpp | 352 +++++++++++++++++++++++ src/hafnian-cpp/hafnian-cpp/lhafnian.h | 12 + src/hafnian-cpp/hafnian-cpp/stdafx.h | 15 + src/hafnian-cpp/hafnian-cpp/timing.cpp | 33 +++ 5 files changed, 419 insertions(+), 2 deletions(-) create mode 100644 src/hafnian-cpp/hafnian-cpp/lhafnian.cpp create mode 100644 src/hafnian-cpp/hafnian-cpp/lhafnian.h create mode 100644 src/hafnian-cpp/hafnian-cpp/stdafx.h create mode 100644 src/hafnian-cpp/hafnian-cpp/timing.cpp diff --git a/.gitignore b/.gitignore index 7be6ba5af..36e75be6f 100644 --- a/.gitignore +++ b/.gitignore @@ -12,5 +12,10 @@ tests/__pycache__/* .coverage coverage_html_report hafnian/lib - - +src/hafnian-cpp/Debug +src/hafnian-cpp/Release +src/hafnian-cpp/hafnian-cpp.* +src/hafnian-cpp/hafnian-cpp/Debug +src/hafnian-cpp/hafnian-cpp/Release +src/hafnian-cpp/hafnian-cpp/*.vcxproj* +src/hafnian-cpp/.vs \ No newline at end of file diff --git a/src/hafnian-cpp/hafnian-cpp/lhafnian.cpp b/src/hafnian-cpp/hafnian-cpp/lhafnian.cpp new file mode 100644 index 000000000..155408acf --- /dev/null +++ b/src/hafnian-cpp/hafnian-cpp/lhafnian.cpp @@ -0,0 +1,352 @@ +// Copyright 2018 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. + +#include "lhafnian.h" +using namespace std; + +/* + Given a complex matrix z of dimensions n x n it calculates traces[k] = tr(z^j) for 1 <= j <= l +*/ +void powtrace(CplxType* z, int n, CplxType* traces, int l) { + CplxType* vals = new CplxType[n]; + CplxType* pvals = new CplxType[n]; + + // work arrays + int lwork = 2 * n; + CplxType* work = new CplxType[lwork]; + double* rwork = new double[n]; + + if (n != 0) { + //evals(z, vals, n, work, lwork, rwork); // TODO: Uncomment and link to actual method? + } + + for (int j = 0; j < n; j++) { + pvals[j] = vals[j]; + } + for (int i = 0; i < l; i++) { + CplxType sum = 0.0; + for (int j = 0; j < n; j++) { + sum += pvals[j]; + } + traces[i] = sum; + for (int j = 0; j < n; j++) { + pvals[j] = pvals[j] * vals[j]; + } + } + + delete[] vals, pvals, work, rwork; +} + +/* + Convert decimal number in x to character vector dst of length len representing binary number +*/ +void dec2bin(char* dst, unsigned long long int x, Byte len) { + char i; // this variable cannot be unsigned + for (i = len - 1; i >= 0; i--) { + *dst++ = x >> i & 1; + } +} + +/* + Given a string of length len it finds in which positions it has a one and stores its position i, + as 2*i and 2*i+1 in consecutive slots of the array pos. + It also returns (twice) the number of ones in array dst +*/ +Byte find2(char* dst, Byte len, Byte* pos) { + Byte j = 0; + for (Byte i = 0; i < len; i++) { + if (1 != dst[i]) + break; + pos[2 * j] = 2 * i; + pos[2 * j + 1] = 2 * i + 1; + j++; + } + return 2 * j; +} + +// Functions for computing Hafnians without considering loops (diagonal elements). + +/* + This function calculates adds parts X to X+chunksize of Cygan and Pilipczuk formula for the Hafnian of matrix mat +*/ +CplxType do_chunk(CplxType* mat, int n, unsigned long long int X, unsigned long long int chunksize) { + CplxType res = 0.0; + Byte m = n / 2; + + #pragma omp parallel for + for (unsigned long long int x = X; x < X + chunksize; x++) { + CplxType summand = 0.0; + CplxType factor; + CplxType powfactor; + + char* dst = new char[m]; + Byte* pos = new Byte[n]; + dec2bin(dst, x, m); + + char sum = find2(dst, m, pos); + CplxType* B = new CplxType[sum * sum]; + for (int i = 0; i < sum; i++) { + for (int j = 0; j < sum; j++) { + B[i * sum + j] = (CplxType) mat[pos[i] * n + ((pos[j]) ^ 1)]; + } + } + + CplxType* traces = new CplxType[m]; + powtrace(B, sum, traces, m); + + char cnt; + Byte cntindex; + cnt = 1; + cntindex = 0; + + CplxType** comb = new CplxType*[2]; + for (int i = 0; i < 2; i++) { + comb[i] = new CplxType[m + 1]; + for (int j = 0; j < m + 1; j++) { + comb[i][j] = 0; + } + } + comb[0][0] = 1.0; + + for (int i = 1; i <= n / 2; i++) { + factor = traces[i - 1] / ((double) 2 * i); + powfactor = 1.0; + //the following line is significantly different than in Octave + //It basically reassings the first column of comb to the second and viceversa + cnt = -cnt; + cntindex = (1 + cnt) / 2; + for (int j = 0; j < n / 2 + 1; j++) { + comb[1 - cntindex][j] = comb[cntindex][j]; + } + for (int j = 1; j <= (n / (2 * i)); j++) { + powfactor = powfactor * factor / (double) j; + for (int k = i * j + 1; k <= n / 2 + 1; k++) { + comb[1 - cntindex][k - 1] += comb[cntindex][k - i * j - 1] * powfactor; + } + } + } + if (((sum / 2) % 2) == (n / 2 % 2)) { + summand = comb[1 - cntindex][n / 2]; + } else { + summand = -comb[1 - cntindex][n / 2]; + } + + #pragma omp critical // TODO: There's better OpenMP constructs for handling this... + res += summand; + + // TODO: Can probably optimise array allocation to be outside the loop -- dimensions are independent of x + delete[] comb[0], comb[1]; + delete[] dst, pos, B, traces, comb; + } + + return res; +} + +/* + This is a wrapper for the function hafnian. Instead of returning the Hafnian of the matrix mat by value + it does by reference with two doubles for the real an imaginary respectively. +*/ +void haf(CplxType* mat, int n, double* res) { + CplxType result = hafnian(mat, n); + res[0] = real(result); + res[1] = imag(result); +} + +void dhaf(double* mat, int n, double *res) { + haf((CplxType*) mat, n, res); // XXX: This looks pretty sketchy -- it doesn't IOOBE in C99?? Presumably does in C++... +} + +/* + Add all the terms necessary to calculate the Hafnian of the matrix mat of size n + It additionally uses OpenMP to parallelize the summation +*/ +CplxType hafnian(CplxType* mat, int n) { + assert(n % 2 == 0); + Byte m = n / 2; + Byte mm = m / 2; + + unsigned long long int pow1 = ((unsigned long long int)pow(2.0, (double)m)); //-1; + unsigned long long int workers = ((unsigned long long int)pow(2.0, (double)mm)); + workers = min(workers, pow1); + + unsigned long long int chunksize = pow1; + unsigned long long int rank = 0; + + return do_chunk(mat, n, rank, chunksize); + +} + +CplxType dhafnian(double *mat, int n) { + return hafnian((CplxType*)mat, n); // XXX: Same sketchiness... +} + +// Functions for computing Hafnians WITH considering loops (diagonal elements). These functions are named with a suffix "_loops" to differentiate them from previous ones. + +/* + This function calculates adds parts X to X+chunksize of Cygan and Pilipczuk formula for the Hafnian of matrix mat +*/ +CplxType do_chunk_loops(CplxType* mat, CplxType* C, CplxType* D, int n, unsigned long long int X, unsigned long long int chunksize) { + CplxType res = 0.0; + + #pragma omp parallel for + for (unsigned long long int x = X * chunksize; x < (X + 1) * chunksize; x++) { + CplxType summand = 0.0; + Byte m = n / 2; + CplxType factor, powfactor; + + char* dst = new char[m]; + Byte* pos = new Byte[n]; + dec2bin(dst, x, m); + Byte sum = find2(dst, m, pos); + + CplxType* B = new CplxType[sum * sum]; + CplxType* B_powtrace = new CplxType[sum * sum]; + CplxType* C1 = new CplxType[sum]; + CplxType* D1 = new CplxType[sum]; + + for (int i = 0; i < sum; i++) { + for (int j = 0; j < sum; j++) { + B[i * sum + j] = mat[pos[i] * n + ((pos[j]) ^ 1)]; + B_powtrace[i * sum + j] = mat[pos[i] * n + ((pos[j]) ^ 1)]; + } + C1[i] = C[pos[i]]; + D1[i] = D[pos[i]]; + } + + CplxType* traces = new CplxType[m]; + powtrace(B_powtrace, sum, traces, m); + char cnt = 1; + Byte cntindex = 0; + CplxType** comb = new CplxType*[2]; + for (int i = 0; i < 2; i++) { + comb[i] = new CplxType[m + 1]; + for (int j = 0; j < m + 1; j++) { + comb[i][j] = 0; + } + } + comb[0][0] = 1.0; + + for (int i = 1; i <= n / 2; i++) { + factor = traces[i - 1] / (2. * i); + + CplxType tmpn = 0.0; + for (int k = 0; k < sum; k++) { + tmpn += C1[k] * D1[k]; + } + factor += 0.5 * tmpn; + + CplxType* tmp_c1 = new CplxType[sum]; + CplxType tmp = 0.0; + for (int k = 0; k < sum; k++) { + tmp = 0.0; + for (int j = 0; j < sum; j++) { + tmp += C1[j] * B[j * sum + k]; + } + tmp_c1[k] = tmp; + C1[k] = tmp_c1[k]; + } + + powfactor = 1.0; + //the following line is significantly different than in Octave + //It basically reassigns the first column of comb to the second and viceversa + cnt = -cnt; + cntindex = (1 + cnt) / 2; + for (int j = 0; j < n / 2 + 1; j++) { + comb[1 - cntindex][j] = comb[cntindex][j]; + } + + for (int j = 1; j <= (n / (2 * i)); j++) { + powfactor = powfactor * factor / (CplxType) j; + for (int k = i * j + 1; k <= n / 2 + 1; k++) { + comb[1 - cntindex][k - 1] = comb[1 - cntindex][k - 1] + comb[cntindex][k - i * j - 1] * powfactor; + } + } + + delete[] tmp_c1; + } + + if (((sum / 2) % 2) == (n / 2 % 2)) { + summand = comb[1 - cntindex][n / 2]; + } else { + summand = -comb[1 - cntindex][n / 2]; + } + #pragma omp critical + res += summand; + + // TODO: More potential array allocation optimisations... + delete[] comb[0], comb[1]; + delete[] dst, pos, B, B_powtrace, C1, D1, traces, comb; + } + + return res; +} + +/* + Add all the terms necessary to calculate the Hafnian of the matrix mat of size n + It additionally uses open MP to parallelize the summation +*/ +CplxType hafnian_loops(CplxType* mat, int n) { + assert(n % 2 == 0); + Byte m = n / 2; + Byte mm = m / 2; + unsigned long long int pow1 = ((unsigned long long int)pow(2.0, (double)m)); //-1; + unsigned long long int workers = ((unsigned long long int)pow(2.0, (double)mm)); + workers = min(workers, pow1); + + unsigned long long int chunksize = pow1; + unsigned long long int rank = 0; + + // CplxType D_mat[n][n] ; + CplxType* C = new CplxType[n]; + CplxType* D = new CplxType[n]; + for (int i = 0; i < n; i++) { + C[i] = 0.0; + D[i] = 0.0; + } + + for (int i = 0; i < n; i++) + { + // D_mat[i][i] = mat[i * n + i]; + D[i] = mat[i * n + i]; + } + + for (int i = 0; i < n; i += 2) { + C[i] = D[i + 1]; + C[i + 1] = D[i]; + } + + CplxType summand = do_chunk_loops(mat, C, D, n, rank, chunksize); + + delete[] C, D; + + return summand; +} + +CplxType dhafnian_loops(double *mat, int n) { + return hafnian_loops((CplxType*) mat, n); // XXX: Yet more sketchiness... +} + +/* + This is a wrapper for the function hafnian. Instead of returning the Hafnian of the matrix + mat by value it does by reference with two doubles for the real an imaginary respectively. +*/ +void haf_loops(CplxType* mat, int n, double* res) { + CplxType result = hafnian_loops(mat, n); + res[0] = result.real(); + res[1] = result.imag(); +} + +void dhaf_loops(double *mat, int n, double *res) { + haf_loops((CplxType*) mat, n, res); // XXX: Viva la sketchiness... +} \ No newline at end of file diff --git a/src/hafnian-cpp/hafnian-cpp/lhafnian.h b/src/hafnian-cpp/hafnian-cpp/lhafnian.h new file mode 100644 index 000000000..650dc107c --- /dev/null +++ b/src/hafnian-cpp/hafnian-cpp/lhafnian.h @@ -0,0 +1,12 @@ +#pragma once + +#include "stdafx.h" +using namespace std; + +typedef complex CplxType; +typedef unsigned char Byte; + +void haf(CplxType* mat, int n, double* res); +CplxType hafnian(CplxType* mat, int n); +//void evals (CplxType* z, CplxType* vals, int n); +CplxType hafnian_loops(CplxType* mat, int n); \ No newline at end of file diff --git a/src/hafnian-cpp/hafnian-cpp/stdafx.h b/src/hafnian-cpp/hafnian-cpp/stdafx.h new file mode 100644 index 000000000..a80b1dd3a --- /dev/null +++ b/src/hafnian-cpp/hafnian-cpp/stdafx.h @@ -0,0 +1,15 @@ +// stdafx.h : include file for standard system include files, +// or project specific include files that are used frequently, but +// are changed infrequently +// + +#pragma once + +#include +#include + +#include +#include +#include +#include +#include \ No newline at end of file diff --git a/src/hafnian-cpp/hafnian-cpp/timing.cpp b/src/hafnian-cpp/hafnian-cpp/timing.cpp new file mode 100644 index 000000000..638e1980d --- /dev/null +++ b/src/hafnian-cpp/hafnian-cpp/timing.cpp @@ -0,0 +1,33 @@ +// Copyright 2018 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. + +#include "lhafnian.h" +using namespace std; + +int main() { + int nmax = 10; + for (int m = 1; m <= nmax; m++) { + int n = 2 * m; + CplxType* mat = new CplxType[n * n]; + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + mat[n * i + j] = 1.0; + } + } + + CplxType hafval = hafnian_loops(mat, n); + printf("%d %lf %lf\n", n, hafval.real(), hafval.imag()); + } + return 0; +} \ No newline at end of file