diff --git a/org/qmckl_sherman_morrison_woodbury.org b/org/qmckl_sherman_morrison_woodbury.org index 000f06a..d562897 100644 --- a/org/qmckl_sherman_morrison_woodbury.org +++ b/org/qmckl_sherman_morrison_woodbury.org @@ -223,13 +223,13 @@ qmckl_exit_code qmckl_sherman_morrison_c(const qmckl_context context, double threshold = 0.0; qmckl_exit_code rc = qmckl_sherman_morrison_threshold_c(&threshold); - unsigned int l = 0; + uint64_t l = 0; // For each update while (l < N_updates) { // C = A^{-1} x U_l - for (unsigned int i = 0; i < Dim; i++) { + for (uint64_t i = 0; i < Dim; i++) { C[i] = 0; - for (unsigned int j = 0; j < Dim; j++) { + for (uint64_t j = 0; j < Dim; j++) { C[i] += Slater_inv[i * Dim + j] * Updates[l * Dim + j]; } } @@ -244,13 +244,13 @@ qmckl_exit_code qmckl_sherman_morrison_c(const qmckl_context context, double iden = 1 / den; // D = v^T x A^{-1} - for (unsigned int j = 0; j < Dim; j++) { + for (uint64_t j = 0; j < Dim; j++) { D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j]; } // A^{-1} = A^{-1} - C x D / den - for (unsigned int i = 0; i < Dim; i++) { - for (unsigned int j = 0; j < Dim; j++) { + for (uint64_t i = 0; i < Dim; i++) { + for (uint64_t j = 0; j < Dim; j++) { double update = C[i] * D[j] * iden; Slater_inv[i * Dim + j] -= update; } @@ -415,16 +415,16 @@ qmckl_exit_code qmckl_woodbury_2_c(const qmckl_context context, // std::cerr << "Called Woodbury 2x2 kernel" << std::endl; // #endif - const unsigned int row1 = (Updates_index[0] - 1); - const unsigned int row2 = (Updates_index[1] - 1); + const uint64_t row1 = (Updates_index[0] - 1); + const uint64_t row2 = (Updates_index[1] - 1); // Compute C = S_inv * U !! NON-STANDARD MATRIX MULTIPLICATION BECAUSE // OF LAYOUT OF 'Updates' !! double C[2 * Dim]; - for (unsigned int i = 0; i < Dim; i++) { - for (unsigned int j = 0; j < 2; j++) { + for (uint64_t i = 0; i < Dim; i++) { + for (uint64_t j = 0; j < 2; j++) { C[i * 2 + j] = 0; - for (unsigned int k = 0; k < Dim; k++) { + for (uint64_t k = 0; k < Dim; k++) { C[i * 2 + j] += Slater_inv[i * Dim + k] * Updates[Dim * j + k]; } } @@ -453,16 +453,16 @@ qmckl_exit_code qmckl_woodbury_2_c(const qmckl_context context, // Compute tmp = B^{-1} x (V.S^{-1}) double tmp[2 * Dim]; - for (unsigned int i = 0; i < 2; i++) { - for (unsigned int j = 0; j < Dim; j++) { + for (uint64_t i = 0; i < 2; i++) { + for (uint64_t j = 0; j < Dim; j++) { tmp[i * Dim + j] = Binv[i * 2] * Slater_inv[row1 * Dim + j]; tmp[i * Dim + j] += Binv[i * 2 + 1] * Slater_inv[row2 * Dim + j]; } } // Compute (S + U V)^{-1} = S^{-1} - C x tmp - for (unsigned int i = 0; i < Dim; i++) { - for (unsigned int j = 0; j < Dim; j++) { + for (uint64_t i = 0; i < Dim; i++) { + for (uint64_t j = 0; j < Dim; j++) { Slater_inv[i * Dim + j] -= C[i * 2] * tmp[j]; Slater_inv[i * Dim + j] -= C[i * 2 + 1] * tmp[Dim + j]; }