write an ITensor to file
#include "itensor/all.h" | |
#include <H5Cpp.h> | |
using namespace H5; | |
using namespace itensor; | |
const std::complex<double> im(0,1); //1i: should be better way | |
/* Export the data tensor A to filename filename; data is reshaped to | |
match the lengths of the indices. | |
Suppose A has indices | |
i1, i2, ..., ir | |
with sizes | |
n1, n2, ..., nr. | |
If A is real, this writes A as an n1 x n2 x ... x nr matrix; if A | |
is complex, this writes A as an n1 x n2 x ... x nr x 2 matrix (real | |
and complex parts). | |
My use case is exporting an I(Q)Tensor of results for analysis in | |
Julia; if you want to do something else you're probably going to | |
have to go learn some HDF5 and write some more code. | |
Based on | |
https://support.hdfgroup.org/ftp/HDF5/current/src/unpacked/c++/examples/h5tutr_rdwt.cpp | |
https://support.hdfgroup.org/ftp/HDF5/current/src/unpacked/c++/examples/h5tutr_crtdat.cpp ; | |
flags are passed straight to H5File initialization function, | |
documented at | |
https://support.hdfgroup.org/HDF5/doc/cpplus_RM/class_h5_1_1_h5_file.html#af25054898de738072217e274217a278c . | |
You probably want H5F_ACC_TRUNC (silently overwrites on existing | |
data) or H5F_ACC_EXCL (errors if file exists). I define versions | |
with default H5F_ACC_TRUNC below. | |
*/ | |
int export_hdf5(ITensor A, std::string filename, int flags) | |
{ | |
int rank = A.r(); | |
std::vector<hsize_t> dims; | |
for(auto i: A.inds()) { dims.push_back(i.m()); } | |
try { | |
Exception::dontPrint(); | |
H5File file(filename, flags); | |
/* This is basically a closure: I have to pass in the indices | |
and file somehow, and I don't know the applyFunction API well | |
enough to know the Right Way to do that. | |
Makes sense for the write function to depend on whether the | |
tensor is real or complex---have to do different things on | |
both the hdf5 side and the itensor side. | |
These guys rely on two nontrivial facts about C++ and one | |
about the HDF5 C++ API: | |
1. a std::complex<double> is in memory a (contiguous) array | |
of two doubles, | |
2. a std::vector<T> stores its data as a contiguous array | |
of Ts, and | |
3. The HDF5 API just takes a pointer to a contiguous array | |
of data and does the requisite reshaping. | |
These three facts together imply that the HDF5 library will | |
take the memory store for a std::vector<std::complex<double>> | |
and just reshape it into an n1 x n2 x ... x nr x 2 array. | |
All of this results from the general C++ culture/design goal | |
of being backwards-compatible with C. Nonetheless, it gives | |
me the *howling fantods*: this sort of magic makes it so, so | |
easy to make mistakes that the compiler can't catch. It's | |
*not OK*.*/ | |
auto writeReal = [&rank,&dims,&file](Dense<Real> const& d) { | |
DataSpace dataspace(rank, dims.data()); | |
DataSet dataset = file.createDataSet("tensor", PredType::NATIVE_DOUBLE, dataspace); | |
dataset.write(d.store.data(), PredType::NATIVE_DOUBLE); | |
return 1; | |
}; | |
auto writeCplx = [&rank,&dims,&file](Dense<Cplx> const& d) { | |
rank += 1; | |
dims.push_back(2); | |
DataSpace dataspace(rank, dims.data()); | |
DataSet dataset = file.createDataSet("tensor", PredType::NATIVE_DOUBLE, dataspace); | |
dataset.write(d.store.data(), PredType::NATIVE_DOUBLE); | |
return 1; | |
}; | |
if (isReal(A)) { applyFunc( writeReal, A.store());} | |
else { applyFunc( writeCplx, A.store());} | |
} | |
catch(FileIException error) { | |
error.printError(); | |
return -1; | |
} | |
catch(DataSetIException error) { | |
error.printError(); | |
return -1; | |
} | |
catch(DataSpaceIException error) { | |
error.printError(); | |
return -1; | |
} | |
return 0; | |
} | |
int export_hdf5(ITensor A, std::string filename) | |
{ return export_hdf5(A, filename, H5F_ACC_TRUNC); } | |
/* export an IQTensor by just turning it into an ITensor */ | |
int export_hdf5(IQTensor A, std::string filename) | |
{ return export_hdf5(toITensor(A), filename, H5F_ACC_TRUNC); } | |
int export_hdf5(IQTensor A, std::string filename, int flags) | |
{ return export_hdf5(toITensor(A), filename, flags); } | |
int main(void) { | |
int N = 3; | |
/* create a tensor to write */ | |
Index ini = Index("i", N, Link, 0); | |
Index inj = Index("j", N, Link, 0); | |
Index ink = Index("k", N, Link, 0); | |
auto tensor = ITensor(ini, inj, ink); | |
/* populate this tensor in such a way that elements are unique and | |
predictable */ | |
for(int i = 1; i <= N; i++){ | |
for(int j = 1; j <= N; j++){ | |
for(int k = 1; k <= N; k++) { | |
tensor.set(ini(i), inj(j), ink(k), (i-1)*pow(N, 2) + (j-1)*N + (k-1)); | |
} | |
} | |
} | |
/* write the tensor */ | |
export_hdf5(tensor, "test-tensor.h5"); | |
tensor *= im; | |
export_hdf5(tensor, "test-tensor.cplx.h5"); | |
return 0; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment