Merge pull request #78 from LLNL/bugfix-hamiltonianreader

Bugfix: Hamiltonian reader for Petsc-parallel Lindblad solver

Author: steffi7574

src/hamiltonianfilereader.cpp

@@ -64,7 +64,9 @@ void HamiltonianFileReader::receiveHsys(Mat& Ad, Mat& Bd){
   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
+  // Now all ranks set their local values
+  PetscInt ilow, iupp;
+  MatGetOwnershipRange(Ad, &ilow, &iupp); 
   for (int i = 0; i < num_entries; i++) {
     PetscInt row = rows[i]; 
     PetscInt col = cols[i];
@@ -75,19 +77,21 @@ void HamiltonianFileReader::receiveHsys(Mat& Ad, Mat& Bd){
     // 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);
+      if (fabs(imag) > 1e-15 && ilow <= row && row < iupp) MatSetValue(Ad, row, col, imag, INSERT_VALUES);
+      if (fabs(real) > 1e-15 && ilow <= row && row < iupp) MatSetValue(Bd, row, col, -real, INSERT_VALUES);
     } else {
-      // Lindblad: Vectorize I_N \kron X - X^T \kron I_N
+      // Lindblad: 
+      // Vectorized Ad = Real(I_N \kron (-iH) - (-iH)^T \kron I_N) = I_N \kron Im(Hsys) - Im(Hsys)^T \kron I_N
+      // Vectorized Bd = Imag(I_N \kron (-iH) - (-iH)^T \kron I_N) = - I_n \kron Real(Hsys) + Real(Hsys)^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);
+        if (fabs(imag) > 1e-15 && ilow <= row && row < iupp) MatSetValue(Ad, rowk, colk, imag, ADD_VALUES);
+        if (fabs(real) > 1e-15 && ilow <= row && row < iupp) MatSetValue(Bd, rowk, colk, -real, ADD_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);
+        if (fabs(imag) > 1e-15 && ilow <= row && row < iupp) MatSetValue(Ad, rowk, colk, -imag, ADD_VALUES);
+        if (fabs(real) > 1e-15 && ilow <= row && row < iupp) MatSetValue(Bd, rowk, colk, real, ADD_VALUES);
       }
     }
   }
@@ -153,7 +157,9 @@ void HamiltonianFileReader::receiveHc(std::vector<Mat>& Ac_vec, std::vector<Mat>
   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
+  // Now all ranks set their local values
+  PetscInt ilow, iupp;
+  MatGetOwnershipRange(Ac_vec[0], &ilow, &iupp); 
   for (int i = 0; i < num_entries; i++) {
     PetscInt row = rows[i]; 
     PetscInt col = cols[i];
@@ -164,19 +170,21 @@ void HamiltonianFileReader::receiveHc(std::vector<Mat>& Ac_vec, std::vector<Mat>
     // 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);
+      if (fabs(imag) > 1e-15 && ilow <= row && row < iupp) MatSetValue(Ac_vec[osc], row, col, imag, ADD_VALUES);
+      if (fabs(real) > 1e-15 && ilow <= row && row < iupp) MatSetValue(Bc_vec[osc], row, col, -real, ADD_VALUES);
     } else {
-      // Lindblad: Vectorize I_N \kron X - X^T \kron I_N
+      // Lindblad: 
+      // Vectorized Ac = Real(I_N \kron (-iH) - (-iH)^T \kron I_N) = I_N \kron Im(Hsys) - Im(Hsys)^T \kron I_N
+      // Vectorized Bc = Imag(I_N \kron (-iH) - (-iH)^T \kron I_N) = - I_n \kron Real(Hsys) + Real(Hsys)^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);
+        if (fabs(imag) > 1e-15 && ilow <= row && row < iupp) MatSetValue(Ac_vec[osc], rowk, colk, imag, ADD_VALUES);
+        if (fabs(real) > 1e-15 && ilow <= row && row < iupp) MatSetValue(Bc_vec[osc], rowk, colk, -real, ADD_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);
+        if (fabs(imag) > 1e-15 && ilow <= row && row < iupp) MatSetValue(Ac_vec[osc], rowk, colk, -imag, ADD_VALUES);
+        if (fabs(real) > 1e-15 && ilow <= row && row < iupp) MatSetValue(Bc_vec[osc], rowk, colk, real, ADD_VALUES);
       }
     }
   }