Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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