Skip to content

Bugfix hamiltonianreader #74

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 4 commits into from
Jul 22, 2025
Merged

Bugfix hamiltonianreader #74

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
64e318b
Bugfix: Hamiltonian reader on multiple processors.
steffi7574 Jul 22, 2025
3c77c42
Extend regtest for hamiltonian reader for multiple parallel initconds.
steffi7574 Jul 22, 2025
bc602be
Broadcast matrix rows/cols with MPIU_INT
steffi7574 Jul 22, 2025
950af22
Broadcast matrix rows/cols with MPIU_INT
steffi7574 Jul 22, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
152 changes: 110 additions & 42 deletions src/hamiltonianfilereader.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,11 @@ void HamiltonianFileReader::receiveHsys(Mat& Ad, Mat& Bd){
MatSetOption(Bd, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);


// Only rank 0 reads the file and sets all values
// Read and broadcast matrix data
std::vector<PetscInt> rows, cols;
std::vector<PetscScalar> real_vals, imag_vals;

// Only rank 0 reads the file
if (mpirank_world == 0) {
std::ifstream infile(hamiltonian_file_Hsys);
if (!infile.is_open()) {
Expand All @@ -36,28 +40,58 @@ void HamiltonianFileReader::receiveHsys(Mat& Ad, Mat& Bd){
PetscInt row, col;
double real, imag;
if (!(iss >> row >> col >> real >> imag)) continue;
// Assemble: Ad = Real(-i*Hsys) = Imag(Hsys)
// Assemble: Bd = Imag(-i*Hsys) = -Real(Hsys)
if (lindbladtype == LindbladType::NONE) {
// Schroedinger
if (fabs(imag) > 1e-15) MatSetValue(Ad, row, col, imag, INSERT_VALUES);
if (fabs(real) > 1e-15) MatSetValue(Bd, row, col, -real, INSERT_VALUES);
} else {
// Lindblad: Vectorize I_N \kron X - X^T \kron I_N
for (PetscInt k = 0; k < dim_rho; k++) {
PetscInt rowk = row + dim_rho * k;
PetscInt colk = col + dim_rho * k;
if (fabs(imag) > 1e-15) MatSetValue(Ad, rowk, colk, imag, INSERT_VALUES);
if (fabs(real) > 1e-15) MatSetValue(Bd, rowk, colk, -real, INSERT_VALUES);
rowk = col * dim_rho + k;
colk = row * dim_rho + k;
if (fabs(imag) > 1e-15) MatSetValue(Ad, rowk, colk, -imag, INSERT_VALUES);
if (fabs(real) > 1e-15) MatSetValue(Bd, rowk, colk, real, INSERT_VALUES);
}
}
rows.push_back(row);
cols.push_back(col);
real_vals.push_back(real);
imag_vals.push_back(imag);
}
infile.close();
}

// Broadcast the size and data to all ranks
int num_entries = rows.size();
MPI_Bcast(&num_entries, 1, MPI_INT, 0, MPI_COMM_WORLD);

if (mpirank_world != 0) {
rows.resize(num_entries);
cols.resize(num_entries);
real_vals.resize(num_entries);
imag_vals.resize(num_entries);
}

MPI_Bcast(rows.data(), num_entries, MPIU_INT, 0, MPI_COMM_WORLD);
MPI_Bcast(cols.data(), num_entries, MPIU_INT, 0, MPI_COMM_WORLD);
MPI_Bcast(real_vals.data(), num_entries, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast(imag_vals.data(), num_entries, MPI_DOUBLE, 0, MPI_COMM_WORLD);

// Now all ranks set the values
for (int i = 0; i < num_entries; i++) {
PetscInt row = rows[i];
PetscInt col = cols[i];
double real = real_vals[i];
double imag = imag_vals[i];

// Assemble: Ad = Real(-i*Hsys) = Imag(Hsys)
// Assemble: Bd = Imag(-i*Hsys) = -Real(Hsys)
if (lindbladtype == LindbladType::NONE) {
// Schroedinger
if (fabs(imag) > 1e-15) MatSetValue(Ad, row, col, imag, INSERT_VALUES);
if (fabs(real) > 1e-15) MatSetValue(Bd, row, col, -real, INSERT_VALUES);
} else {
// Lindblad: Vectorize I_N \kron X - X^T \kron I_N
for (PetscInt k = 0; k < dim_rho; k++) {
PetscInt rowk = row + dim_rho * k;
PetscInt colk = col + dim_rho * k;
if (fabs(imag) > 1e-15) MatSetValue(Ad, rowk, colk, imag, INSERT_VALUES);
if (fabs(real) > 1e-15) MatSetValue(Bd, rowk, colk, -real, INSERT_VALUES);
rowk = col * dim_rho + k;
colk = row * dim_rho + k;
if (fabs(imag) > 1e-15) MatSetValue(Ad, rowk, colk, -imag, INSERT_VALUES);
if (fabs(real) > 1e-15) MatSetValue(Bd, rowk, colk, real, INSERT_VALUES);
}
}
}

MatAssemblyBegin(Ad, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(Ad, MAT_FINAL_ASSEMBLY);
MatAssemblyBegin(Bd, MAT_FINAL_ASSEMBLY);
Expand All @@ -73,7 +107,11 @@ void HamiltonianFileReader::receiveHc(std::vector<Mat>& Ac_vec, std::vector<Mat>
MatSetOption(Bc_vec[i], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
}

// Only rank 0 reads the file and sets all values
// Read and broadcast matrix data
std::vector<PetscInt> rows, cols, oscs;
std::vector<PetscScalar> real_vals, imag_vals;

// Only rank 0 reads the file
if (mpirank_world == 0) {
std::ifstream infile(hamiltonian_file_Hc);
if (!infile.is_open()) {
Expand All @@ -88,31 +126,61 @@ void HamiltonianFileReader::receiveHc(std::vector<Mat>& Ac_vec, std::vector<Mat>
PetscInt row, col;
double real, imag;
if (!(iss >> osc >> row >> col >> real >> imag)) continue;
// printf("osc %d row %d col %d real %f imag %f\n", osc, row, col, real, imag);
// Assemble: Ac = Real(-i*Hc) = Imag(Hc)
// Assemble: Bc = Imag(-i*Hc) = -Real(Hc)
if (lindbladtype == LindbladType::NONE) {
// Schroedinger
if (fabs(imag) > 1e-15) MatSetValue(Ac_vec[osc], row, col, imag, INSERT_VALUES);
if (fabs(real) > 1e-15) MatSetValue(Bc_vec[osc], row, col, -real, INSERT_VALUES);
} else {
// Lindblad: Vectorize I_N \kron X - X^T \kron I_N
for (PetscInt k = 0; k < dim_rho; k++) {
PetscInt rowk = row + dim_rho * k;
PetscInt colk = col + dim_rho * k;
if (fabs(imag) > 1e-15) MatSetValue(Ac_vec[osc], rowk, colk, imag, INSERT_VALUES);
if (fabs(real) > 1e-15) MatSetValue(Bc_vec[osc], rowk, colk, -real, INSERT_VALUES);
rowk = col * dim_rho + k;
colk = row * dim_rho + k;
if (fabs(imag) > 1e-15) MatSetValue(Ac_vec[osc], rowk, colk, -imag, INSERT_VALUES);
if (fabs(real) > 1e-15) MatSetValue(Bc_vec[osc], rowk, colk, real, INSERT_VALUES);
}
}
oscs.push_back(osc);
rows.push_back(row);
cols.push_back(col);
real_vals.push_back(real);
imag_vals.push_back(imag);
}
infile.close();
}

// All ranks participate in assembly for all oscillators
// Broadcast the size and data to all ranks
int num_entries = rows.size();
MPI_Bcast(&num_entries, 1, MPI_INT, 0, MPI_COMM_WORLD);

if (mpirank_world != 0) {
rows.resize(num_entries);
cols.resize(num_entries);
oscs.resize(num_entries);
real_vals.resize(num_entries);
imag_vals.resize(num_entries);
}

MPI_Bcast(rows.data(), num_entries, MPIU_INT, 0, MPI_COMM_WORLD);
MPI_Bcast(cols.data(), num_entries, MPIU_INT, 0, MPI_COMM_WORLD);
MPI_Bcast(oscs.data(), num_entries, MPIU_INT, 0, MPI_COMM_WORLD);
MPI_Bcast(real_vals.data(), num_entries, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast(imag_vals.data(), num_entries, MPI_DOUBLE, 0, MPI_COMM_WORLD);

// Now all ranks set the values
for (int i = 0; i < num_entries; i++) {
PetscInt row = rows[i];
PetscInt col = cols[i];
PetscInt osc = oscs[i];
double real = real_vals[i];
double imag = imag_vals[i];
// Assemble: Ac = Real(-i*Hc) = Imag(Hc)
// Assemble: Bc = Imag(-i*Hc) = -Real(Hc)
if (lindbladtype == LindbladType::NONE) {
// Schroedinger
if (fabs(imag) > 1e-15) MatSetValue(Ac_vec[osc], row, col, imag, INSERT_VALUES);
if (fabs(real) > 1e-15) MatSetValue(Bc_vec[osc], row, col, -real, INSERT_VALUES);
} else {
// Lindblad: Vectorize I_N \kron X - X^T \kron I_N
for (PetscInt k = 0; k < dim_rho; k++) {
PetscInt rowk = row + dim_rho * k;
PetscInt colk = col + dim_rho * k;
if (fabs(imag) > 1e-15) MatSetValue(Ac_vec[osc], rowk, colk, imag, INSERT_VALUES);
if (fabs(real) > 1e-15) MatSetValue(Bc_vec[osc], rowk, colk, -real, INSERT_VALUES);
rowk = col * dim_rho + k;
colk = row * dim_rho + k;
if (fabs(imag) > 1e-15) MatSetValue(Ac_vec[osc], rowk, colk, -imag, INSERT_VALUES);
if (fabs(real) > 1e-15) MatSetValue(Bc_vec[osc], rowk, colk, real, INSERT_VALUES);
}
}
}

for (size_t i = 0; i < Ac_vec.size(); ++i) {
MatAssemblyBegin(Ac_vec[i], MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(Ac_vec[i], MAT_FINAL_ASSEMBLY);
Expand Down
Loading
Loading