From 2ff3ee061e55b55f93626c1920da990f237bec59 Mon Sep 17 00:00:00 2001 From: Chris Vales Date: Wed, 27 Sep 2023 08:32:39 -0700 Subject: [PATCH 1/7] fix issues in Matrix::orthogonalize --- lib/linalg/Matrix.cpp | 55 +++++++++++++++++++++----------------- unit_tests/test_Matrix.cpp | 54 +++++++++++++++++++++++++++++++++++++ 2 files changed, 85 insertions(+), 24 deletions(-) diff --git a/lib/linalg/Matrix.cpp b/lib/linalg/Matrix.cpp index d0b6514d3..ac2d0a810 100644 --- a/lib/linalg/Matrix.cpp +++ b/lib/linalg/Matrix.cpp @@ -1790,44 +1790,51 @@ const void Matrix::orthogonalize() { - for (int work = 1; work < d_num_cols; ++work) { + for (int work = 0; work < d_num_cols; ++work) + { + // Orthogonalize the column. double tmp; - for (int col = 0; col < work; ++col) { + for (int col = 0; col < work; ++col) + { double factor = 0.0; tmp = 0.0; - for (int i = 0; i < d_num_rows; ++i) { - tmp += item(i, col)*item(i, work); - } - if (d_num_procs > 1) { - MPI_Allreduce(&tmp, - &factor, - 1, - MPI_DOUBLE, - MPI_SUM, + + for (int i = 0; i < d_num_rows; ++i) + tmp += item(i, col) * item(i, work); + + if (d_num_procs > 1) + { + MPI_Allreduce(&tmp, &factor, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); } - else { + else + { factor = tmp; } - - for (int i = 0; i < d_num_rows; ++i) { - item(i, work) -= factor*item(i, col); - } + for (int i = 0; i < d_num_rows; ++i) + item(i, work) -= factor * item(i, col); } + + // Normalize the column. double norm = 0.0; tmp = 0.0; - for (int i = 0; i < d_num_rows; ++i) { - tmp += item(i, work)*item(i, work); - } - if (d_num_procs > 1) { + + for (int i = 0; i < d_num_rows; ++i) + tmp += item(i, work) * item(i, work); + + if (d_num_procs > 1) + { MPI_Allreduce(&tmp, &norm, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); } - else { + else + { norm = tmp; } - norm = sqrt(norm); - for (int i = 0; i < d_num_rows; ++i) { - item(i, work) /= norm; + if (norm > 1.0e-15) + { + norm = 1.0 / sqrt(norm); + for (int i = 0; i < d_num_rows; ++i) + item(i, work) *= norm; } } } diff --git a/unit_tests/test_Matrix.cpp b/unit_tests/test_Matrix.cpp index 7534d5b19..7f86486e6 100644 --- a/unit_tests/test_Matrix.cpp +++ b/unit_tests/test_Matrix.cpp @@ -457,6 +457,60 @@ TEST(MatrixSerialTest, Test_Matrix_rescale_cols_max2) "(i, j) = (" << i << ", " << j << ")"; } +TEST(MatrixSerialTest, Test_Matrix_orthogonalize) +{ + // Matrix data to orthonormalize. + double d_mat[16] = {3.5, 7.1, 0.0, 0.0, + 0.0, 1.9, 8.3, 0.0, + 0.0, 0.0, 5.7, 4.6, + 0.0, 0.0, 0.0, 3.2}; + + // Target matrix data. + double d_mat2[16] = {1.0, 0.0, 0.0, 0.0, + 0.0, 1.0, 0.0, 0.0, + 0.0, 0.0, 1.0, 0.0, + 0.0, 0.0, 0.0, 1.0}; + + CAROM::Matrix matrix(d_mat, 4, 4, false); + CAROM::Matrix target(d_mat2, 4, 4, false); + + double abs_error = 1.0e-15; // absolute error threshold + + matrix.orthogonalize(); + + for (int i = 0; i < 4; i++) + for (int j = 0; j < 4; j++) + EXPECT_NEAR(matrix.item(i, j), target.item(i, j), abs_error) << + "(i, j) = (" << i << ", " << j << ")"; +} + +TEST(MatrixSerialTest, Test_Matrix_orthogonalize2) +{ + // Matrix data to orthonormalize. + double d_mat[16] = {3.5, 7.1, 0.0, 0.0, + 0.0, 1.9, 8.3, 1e-14, + 0.0, 0.0, 5.7, 1.0+1.0e-14, + 0.0, 0.0, 0.0, 0.0}; + + // Target matrix data. + double d_mat2[16] = {1.0, 0.0, 0.0, 0.0, + 0.0, 1.0, 0.0, 0.0, + 0.0, 0.0, 1.0, 0.0, + 0.0, 0.0, 0.0, 0.0}; + + CAROM::Matrix matrix(d_mat, 4, 4, false); + CAROM::Matrix target(d_mat2, 4, 4, false); + + double abs_error = 1.0e-15; // absolute error threshold + + matrix.orthogonalize(); + + for (int i = 0; i < 4; i++) + for (int j = 0; j < 4; j++) + EXPECT_NEAR(matrix.item(i, j), target.item(i, j), abs_error) << + "(i, j) = (" << i << ", " << j << ")"; +} + TEST(MatrixSerialTest, Test_pMatrix_mult_reference) { /** From 698b15792c37f0693dea90f07cc06111a138c227 Mon Sep 17 00:00:00 2001 From: Chris Vales Date: Wed, 27 Sep 2023 08:34:40 -0700 Subject: [PATCH 2/7] add Matrix::orthogonalize_last method --- lib/linalg/Matrix.cpp | 54 +++++++++++++++++++++++++++++++++++++ lib/linalg/Matrix.h | 13 +++++++++ unit_tests/test_Matrix.cpp | 55 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 122 insertions(+) diff --git a/lib/linalg/Matrix.cpp b/lib/linalg/Matrix.cpp index ac2d0a810..8274fce69 100644 --- a/lib/linalg/Matrix.cpp +++ b/lib/linalg/Matrix.cpp @@ -1839,6 +1839,60 @@ Matrix::orthogonalize() } } +void +Matrix::orthogonalize_last(int ncols) +{ + if (ncols == -1) ncols = d_num_cols; + CAROM_ASSERT((ncols > 0) && (ncols <= d_num_cols)); + + const int last_col = ncols - 1; // index of column to be orthonormalized + double tmp; + + // Orthogonalize the column. + for (int col = 0; col < last_col; ++col) + { + double factor = 0.0; + tmp = 0.0; + + for (int i = 0; i < d_num_rows; ++i) + tmp += item(i, col) * item(i, last_col); + + if (d_num_procs > 1) + { + MPI_Allreduce(&tmp, &factor, 1, MPI_DOUBLE, MPI_SUM, + MPI_COMM_WORLD); + } + else + { + factor = tmp; + } + for (int i = 0; i < d_num_rows; ++i) + item(i, last_col) -= factor * item(i, col); + } + + // Normalize the column. + double norm = 0.0; + tmp = 0.0; + + for (int i = 0; i < d_num_rows; ++i) + tmp += item(i, last_col) * item(i, last_col); + + if (d_num_procs > 1) + { + MPI_Allreduce(&tmp, &norm, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + } + else + { + norm = tmp; + } + if (norm > 1.0e-15) + { + norm = 1.0 / sqrt(norm); + for (int i = 0; i < d_num_rows; ++i) + item(i, last_col) *= norm; + } +} + void Matrix::rescale_rows_max() { diff --git a/lib/linalg/Matrix.h b/lib/linalg/Matrix.h index 4f7f74aba..3cd6925b8 100644 --- a/lib/linalg/Matrix.h +++ b/lib/linalg/Matrix.h @@ -946,6 +946,19 @@ class Matrix void orthogonalize(); + /** + * @brief Orthonormalizes the matrix's last column, assuming the previous + * columns are already orthonormal. + * + * By default, the function considers the whole matrix. + * If input parameter ncols < d_num_cols, then a subset of the matrix + * is considered. + * This allows one to reorthonormalize the matrix every time a new column + * is added, assuming the previous columns have remained unchanged. + */ + void + orthogonalize_last(int ncols = -1); + /** * @brief Rescale every matrix row by its maximum absolute value. */ diff --git a/unit_tests/test_Matrix.cpp b/unit_tests/test_Matrix.cpp index 7f86486e6..0f111c812 100644 --- a/unit_tests/test_Matrix.cpp +++ b/unit_tests/test_Matrix.cpp @@ -511,6 +511,61 @@ TEST(MatrixSerialTest, Test_Matrix_orthogonalize2) "(i, j) = (" << i << ", " << j << ")"; } +TEST(MatrixSerialTest, Test_Matrix_orthogonalize_last) +{ + // Matrix data to orthonormalize. + double d_mat[16] = {1.0, 0.0, 0.0, 1.3, + 0.0, 1.0, 0.0, 4.7, + 0.0, 0.0, 1.0, 2.5, + 0.0, 0.0, 0.0, 7.3}; + + // Target matrix data. + double d_mat2[16] = {1.0, 0.0, 0.0, 0.0, + 0.0, 1.0, 0.0, 0.0, + 0.0, 0.0, 1.0, 0.0, + 0.0, 0.0, 0.0, 1.0}; + + CAROM::Matrix matrix(d_mat, 4, 4, false); + CAROM::Matrix target(d_mat2, 4, 4, false); + + double abs_error = 1.0e-15; // absolute error threshold + + matrix.orthogonalize_last(); + + for (int i = 0; i < 4; i++) + for (int j = 0; j < 4; j++) + EXPECT_NEAR(matrix.item(i, j), target.item(i, j), abs_error) << + "(i, j) = (" << i << ", " << j << ")"; +} + +TEST(MatrixSerialTest, Test_Matrix_orthogonalize_last2) +{ + // Matrix data to orthonormalize. + double d_mat[16] = {1.0, 0.0, 3.8, 1.3, + 0.0, 1.0, 5.6, 4.7, + 0.0, 0.0, 9.8, 2.5, + 0.0, 0.0, 0.0, 7.3}; + + // Target matrix data. + double d_mat2[16] = {1.0, 0.0, 0.0, 0.0, + 0.0, 1.0, 0.0, 0.0, + 0.0, 0.0, 1.0, 0.0, + 0.0, 0.0, 0.0, 1.0}; + + CAROM::Matrix matrix(d_mat, 4, 4, false); + CAROM::Matrix target(d_mat2, 4, 4, false); + + double abs_error = 1.0e-15; // absolute error threshold + + matrix.orthogonalize_last(3); + matrix.orthogonalize_last(4); + + for (int i = 0; i < 4; i++) + for (int j = 0; j < 4; j++) + EXPECT_NEAR(matrix.item(i, j), target.item(i, j), abs_error) << + "(i, j) = (" << i << ", " << j << ")"; +} + TEST(MatrixSerialTest, Test_pMatrix_mult_reference) { /** From d1c65a22871067f9b154556e5023a163ac73f7f2 Mon Sep 17 00:00:00 2001 From: Chris Vales Date: Wed, 27 Sep 2023 11:40:08 -0400 Subject: [PATCH 3/7] run astyle --- lib/linalg/Matrix.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/linalg/Matrix.cpp b/lib/linalg/Matrix.cpp index 8274fce69..a8877d75a 100644 --- a/lib/linalg/Matrix.cpp +++ b/lib/linalg/Matrix.cpp @@ -1818,7 +1818,7 @@ Matrix::orthogonalize() // Normalize the column. double norm = 0.0; tmp = 0.0; - + for (int i = 0; i < d_num_rows; ++i) tmp += item(i, work) * item(i, work); @@ -1844,7 +1844,7 @@ Matrix::orthogonalize_last(int ncols) { if (ncols == -1) ncols = d_num_cols; CAROM_ASSERT((ncols > 0) && (ncols <= d_num_cols)); - + const int last_col = ncols - 1; // index of column to be orthonormalized double tmp; From 9cbc97021f89cfe4fbbaece817c96ffa3080eaa8 Mon Sep 17 00:00:00 2001 From: Chris Vales Date: Wed, 27 Sep 2023 09:52:57 -0700 Subject: [PATCH 4/7] add another test for Matrix::orthogonalize_last --- unit_tests/test_Matrix.cpp | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/unit_tests/test_Matrix.cpp b/unit_tests/test_Matrix.cpp index 0f111c812..557cc243c 100644 --- a/unit_tests/test_Matrix.cpp +++ b/unit_tests/test_Matrix.cpp @@ -566,6 +566,34 @@ TEST(MatrixSerialTest, Test_Matrix_orthogonalize_last2) "(i, j) = (" << i << ", " << j << ")"; } +TEST(MatrixSerialTest, Test_Matrix_orthogonalize_last3) +{ + // Matrix data to orthonormalize. + double d_mat[16] = {3.5, 7.1, 0.0, 0.0, + 0.0, 1.9, 8.3, 0.0, + 0.0, 0.0, 5.7, 4.6, + 0.0, 0.0, 0.0, 3.2}; + + // Target matrix data. + double d_mat2[16] = {1.0, 0.0, 0.0, 0.0, + 0.0, 1.0, 0.0, 0.0, + 0.0, 0.0, 1.0, 0.0, + 0.0, 0.0, 0.0, 1.0}; + + CAROM::Matrix matrix(d_mat, 4, 4, false); + CAROM::Matrix target(d_mat2, 4, 4, false); + + for (int i = 0; i < 4; i++) + matrix.orthogonalize_last(i+1); + + double abs_error = 1.0e-15; // absolute error threshold + + for (int i = 0; i < 4; i++) + for (int j = 0; j < 4; j++) + EXPECT_NEAR(matrix.item(i, j), target.item(i, j), abs_error) << + "(i, j) = (" << i << ", " << j << ")"; +} + TEST(MatrixSerialTest, Test_pMatrix_mult_reference) { /** From 6f1a113712490472a11194c8f93cb74cc228671d Mon Sep 17 00:00:00 2001 From: Chris Vales Date: Wed, 27 Sep 2023 10:10:13 -0700 Subject: [PATCH 5/7] rerun astyle properly --- lib/linalg/Matrix.h | 20 +++++++-------- unit_tests/test_Matrix.cpp | 50 +++++++++++++++++++++++--------------- 2 files changed, 40 insertions(+), 30 deletions(-) diff --git a/lib/linalg/Matrix.h b/lib/linalg/Matrix.h index 3cd6925b8..23a86d795 100644 --- a/lib/linalg/Matrix.h +++ b/lib/linalg/Matrix.h @@ -948,16 +948,16 @@ class Matrix /** * @brief Orthonormalizes the matrix's last column, assuming the previous - * columns are already orthonormal. - * - * By default, the function considers the whole matrix. - * If input parameter ncols < d_num_cols, then a subset of the matrix - * is considered. - * This allows one to reorthonormalize the matrix every time a new column - * is added, assuming the previous columns have remained unchanged. - */ - void - orthogonalize_last(int ncols = -1); + * columns are already orthonormal. + * + * By default, the function considers the whole matrix. + * If input parameter ncols < d_num_cols, then a subset of the matrix + * is considered. + * This allows one to reorthonormalize the matrix every time a new column + * is added, assuming the previous columns have remained unchanged. + */ + void + orthogonalize_last(int ncols = -1); /** * @brief Rescale every matrix row by its maximum absolute value. diff --git a/unit_tests/test_Matrix.cpp b/unit_tests/test_Matrix.cpp index 557cc243c..8f3d77503 100644 --- a/unit_tests/test_Matrix.cpp +++ b/unit_tests/test_Matrix.cpp @@ -463,13 +463,15 @@ TEST(MatrixSerialTest, Test_Matrix_orthogonalize) double d_mat[16] = {3.5, 7.1, 0.0, 0.0, 0.0, 1.9, 8.3, 0.0, 0.0, 0.0, 5.7, 4.6, - 0.0, 0.0, 0.0, 3.2}; + 0.0, 0.0, 0.0, 3.2 + }; // Target matrix data. double d_mat2[16] = {1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, - 0.0, 0.0, 0.0, 1.0}; + 0.0, 0.0, 0.0, 1.0 + }; CAROM::Matrix matrix(d_mat, 4, 4, false); CAROM::Matrix target(d_mat2, 4, 4, false); @@ -480,8 +482,8 @@ TEST(MatrixSerialTest, Test_Matrix_orthogonalize) for (int i = 0; i < 4; i++) for (int j = 0; j < 4; j++) - EXPECT_NEAR(matrix.item(i, j), target.item(i, j), abs_error) << - "(i, j) = (" << i << ", " << j << ")"; + EXPECT_NEAR(matrix.item(i, j), target.item(i, j), abs_error) << + "(i, j) = (" << i << ", " << j << ")"; } TEST(MatrixSerialTest, Test_Matrix_orthogonalize2) @@ -490,13 +492,15 @@ TEST(MatrixSerialTest, Test_Matrix_orthogonalize2) double d_mat[16] = {3.5, 7.1, 0.0, 0.0, 0.0, 1.9, 8.3, 1e-14, 0.0, 0.0, 5.7, 1.0+1.0e-14, - 0.0, 0.0, 0.0, 0.0}; + 0.0, 0.0, 0.0, 0.0 + }; // Target matrix data. double d_mat2[16] = {1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, - 0.0, 0.0, 0.0, 0.0}; + 0.0, 0.0, 0.0, 0.0 + }; CAROM::Matrix matrix(d_mat, 4, 4, false); CAROM::Matrix target(d_mat2, 4, 4, false); @@ -507,8 +511,8 @@ TEST(MatrixSerialTest, Test_Matrix_orthogonalize2) for (int i = 0; i < 4; i++) for (int j = 0; j < 4; j++) - EXPECT_NEAR(matrix.item(i, j), target.item(i, j), abs_error) << - "(i, j) = (" << i << ", " << j << ")"; + EXPECT_NEAR(matrix.item(i, j), target.item(i, j), abs_error) << + "(i, j) = (" << i << ", " << j << ")"; } TEST(MatrixSerialTest, Test_Matrix_orthogonalize_last) @@ -517,13 +521,15 @@ TEST(MatrixSerialTest, Test_Matrix_orthogonalize_last) double d_mat[16] = {1.0, 0.0, 0.0, 1.3, 0.0, 1.0, 0.0, 4.7, 0.0, 0.0, 1.0, 2.5, - 0.0, 0.0, 0.0, 7.3}; + 0.0, 0.0, 0.0, 7.3 + }; // Target matrix data. double d_mat2[16] = {1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, - 0.0, 0.0, 0.0, 1.0}; + 0.0, 0.0, 0.0, 1.0 + }; CAROM::Matrix matrix(d_mat, 4, 4, false); CAROM::Matrix target(d_mat2, 4, 4, false); @@ -534,8 +540,8 @@ TEST(MatrixSerialTest, Test_Matrix_orthogonalize_last) for (int i = 0; i < 4; i++) for (int j = 0; j < 4; j++) - EXPECT_NEAR(matrix.item(i, j), target.item(i, j), abs_error) << - "(i, j) = (" << i << ", " << j << ")"; + EXPECT_NEAR(matrix.item(i, j), target.item(i, j), abs_error) << + "(i, j) = (" << i << ", " << j << ")"; } TEST(MatrixSerialTest, Test_Matrix_orthogonalize_last2) @@ -544,13 +550,15 @@ TEST(MatrixSerialTest, Test_Matrix_orthogonalize_last2) double d_mat[16] = {1.0, 0.0, 3.8, 1.3, 0.0, 1.0, 5.6, 4.7, 0.0, 0.0, 9.8, 2.5, - 0.0, 0.0, 0.0, 7.3}; + 0.0, 0.0, 0.0, 7.3 + }; // Target matrix data. double d_mat2[16] = {1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, - 0.0, 0.0, 0.0, 1.0}; + 0.0, 0.0, 0.0, 1.0 + }; CAROM::Matrix matrix(d_mat, 4, 4, false); CAROM::Matrix target(d_mat2, 4, 4, false); @@ -562,8 +570,8 @@ TEST(MatrixSerialTest, Test_Matrix_orthogonalize_last2) for (int i = 0; i < 4; i++) for (int j = 0; j < 4; j++) - EXPECT_NEAR(matrix.item(i, j), target.item(i, j), abs_error) << - "(i, j) = (" << i << ", " << j << ")"; + EXPECT_NEAR(matrix.item(i, j), target.item(i, j), abs_error) << + "(i, j) = (" << i << ", " << j << ")"; } TEST(MatrixSerialTest, Test_Matrix_orthogonalize_last3) @@ -572,13 +580,15 @@ TEST(MatrixSerialTest, Test_Matrix_orthogonalize_last3) double d_mat[16] = {3.5, 7.1, 0.0, 0.0, 0.0, 1.9, 8.3, 0.0, 0.0, 0.0, 5.7, 4.6, - 0.0, 0.0, 0.0, 3.2}; + 0.0, 0.0, 0.0, 3.2 + }; // Target matrix data. double d_mat2[16] = {1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, - 0.0, 0.0, 0.0, 1.0}; + 0.0, 0.0, 0.0, 1.0 + }; CAROM::Matrix matrix(d_mat, 4, 4, false); CAROM::Matrix target(d_mat2, 4, 4, false); @@ -590,8 +600,8 @@ TEST(MatrixSerialTest, Test_Matrix_orthogonalize_last3) for (int i = 0; i < 4; i++) for (int j = 0; j < 4; j++) - EXPECT_NEAR(matrix.item(i, j), target.item(i, j), abs_error) << - "(i, j) = (" << i << ", " << j << ")"; + EXPECT_NEAR(matrix.item(i, j), target.item(i, j), abs_error) << + "(i, j) = (" << i << ", " << j << ")"; } TEST(MatrixSerialTest, Test_pMatrix_mult_reference) From f4893d9aa7e6989b8ad24f8b64bcab69d07b5e15 Mon Sep 17 00:00:00 2001 From: Chris Vales Date: Thu, 28 Sep 2023 09:20:57 -0700 Subject: [PATCH 6/7] review n.1 changes --- lib/linalg/Matrix.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/lib/linalg/Matrix.cpp b/lib/linalg/Matrix.cpp index a8877d75a..5842ff9e8 100644 --- a/lib/linalg/Matrix.cpp +++ b/lib/linalg/Matrix.cpp @@ -1802,7 +1802,7 @@ Matrix::orthogonalize() for (int i = 0; i < d_num_rows; ++i) tmp += item(i, col) * item(i, work); - if (d_num_procs > 1) + if (d_distributed && d_num_procs > 1) { MPI_Allreduce(&tmp, &factor, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); @@ -1822,7 +1822,7 @@ Matrix::orthogonalize() for (int i = 0; i < d_num_rows; ++i) tmp += item(i, work) * item(i, work); - if (d_num_procs > 1) + if (d_distributed && d_num_procs > 1) { MPI_Allreduce(&tmp, &norm, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); } @@ -1857,7 +1857,7 @@ Matrix::orthogonalize_last(int ncols) for (int i = 0; i < d_num_rows; ++i) tmp += item(i, col) * item(i, last_col); - if (d_num_procs > 1) + if (d_distributed && d_num_procs > 1) { MPI_Allreduce(&tmp, &factor, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); @@ -1877,7 +1877,7 @@ Matrix::orthogonalize_last(int ncols) for (int i = 0; i < d_num_rows; ++i) tmp += item(i, last_col) * item(i, last_col); - if (d_num_procs > 1) + if (d_distributed && d_num_procs > 1) { MPI_Allreduce(&tmp, &norm, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); } From 419ba382e4924135c5788b84a4e0a57afd1e8390 Mon Sep 17 00:00:00 2001 From: Chris Vales Date: Wed, 1 Nov 2023 08:27:41 -0700 Subject: [PATCH 7/7] review n.2 --- lib/linalg/Matrix.cpp | 54 ++++++++++++++----------------------------- lib/linalg/Matrix.h | 20 ++++++++++++---- 2 files changed, 32 insertions(+), 42 deletions(-) diff --git a/lib/linalg/Matrix.cpp b/lib/linalg/Matrix.cpp index 5842ff9e8..4918b1ba1 100644 --- a/lib/linalg/Matrix.cpp +++ b/lib/linalg/Matrix.cpp @@ -1788,28 +1788,22 @@ const } void -Matrix::orthogonalize() +Matrix::orthogonalize(double zero_tol) { for (int work = 0; work < d_num_cols; ++work) { // Orthogonalize the column. - double tmp; for (int col = 0; col < work; ++col) { double factor = 0.0; - tmp = 0.0; for (int i = 0; i < d_num_rows; ++i) - tmp += item(i, col) * item(i, work); + factor += item(i, col) * item(i, work); if (d_distributed && d_num_procs > 1) { - MPI_Allreduce(&tmp, &factor, 1, MPI_DOUBLE, MPI_SUM, - MPI_COMM_WORLD); - } - else - { - factor = tmp; + CAROM_VERIFY( MPI_Allreduce(MPI_IN_PLACE, &factor, 1, + MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD) == MPI_SUCCESS ); } for (int i = 0; i < d_num_rows; ++i) item(i, work) -= factor * item(i, col); @@ -1817,20 +1811,16 @@ Matrix::orthogonalize() // Normalize the column. double norm = 0.0; - tmp = 0.0; for (int i = 0; i < d_num_rows; ++i) - tmp += item(i, work) * item(i, work); + norm += item(i, work) * item(i, work); if (d_distributed && d_num_procs > 1) { - MPI_Allreduce(&tmp, &norm, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + CAROM_VERIFY( MPI_Allreduce(MPI_IN_PLACE, &norm, 1, MPI_DOUBLE, + MPI_SUM, MPI_COMM_WORLD) == MPI_SUCCESS ); } - else - { - norm = tmp; - } - if (norm > 1.0e-15) + if (norm > zero_tol) { norm = 1.0 / sqrt(norm); for (int i = 0; i < d_num_rows; ++i) @@ -1840,31 +1830,25 @@ Matrix::orthogonalize() } void -Matrix::orthogonalize_last(int ncols) +Matrix::orthogonalize_last(int ncols, double zero_tol) { if (ncols == -1) ncols = d_num_cols; - CAROM_ASSERT((ncols > 0) && (ncols <= d_num_cols)); + CAROM_VERIFY((ncols > 0) && (ncols <= d_num_cols)); const int last_col = ncols - 1; // index of column to be orthonormalized - double tmp; // Orthogonalize the column. for (int col = 0; col < last_col; ++col) { double factor = 0.0; - tmp = 0.0; for (int i = 0; i < d_num_rows; ++i) - tmp += item(i, col) * item(i, last_col); + factor += item(i, col) * item(i, last_col); if (d_distributed && d_num_procs > 1) { - MPI_Allreduce(&tmp, &factor, 1, MPI_DOUBLE, MPI_SUM, - MPI_COMM_WORLD); - } - else - { - factor = tmp; + CAROM_VERIFY( MPI_Allreduce(MPI_IN_PLACE, &factor, 1, MPI_DOUBLE, + MPI_SUM, MPI_COMM_WORLD) == MPI_SUCCESS ); } for (int i = 0; i < d_num_rows; ++i) item(i, last_col) -= factor * item(i, col); @@ -1872,20 +1856,16 @@ Matrix::orthogonalize_last(int ncols) // Normalize the column. double norm = 0.0; - tmp = 0.0; for (int i = 0; i < d_num_rows; ++i) - tmp += item(i, last_col) * item(i, last_col); + norm += item(i, last_col) * item(i, last_col); if (d_distributed && d_num_procs > 1) { - MPI_Allreduce(&tmp, &norm, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - } - else - { - norm = tmp; + CAROM_VERIFY( MPI_Allreduce(MPI_IN_PLACE, &norm, 1, MPI_DOUBLE, + MPI_SUM, MPI_COMM_WORLD) == MPI_SUCCESS ); } - if (norm > 1.0e-15) + if (norm > zero_tol) { norm = 1.0 / sqrt(norm); for (int i = 0; i < d_num_rows; ++i) diff --git a/lib/linalg/Matrix.h b/lib/linalg/Matrix.h index 23a86d795..e5f5bc988 100644 --- a/lib/linalg/Matrix.h +++ b/lib/linalg/Matrix.h @@ -941,23 +941,33 @@ class Matrix int pivots_requested) const; /** - * @brief Orthogonalizes the matrix. + * @brief Orthonormalizes the matrix. + * + * If the norm of a matrix column is below the value of zero_tol then it + * is considered to be zero, and we do not divide by it. + * Therefore, that column is considered to be zero and is not normalized. + * By default, zero_tol = 1.0e-15. */ void - orthogonalize(); + orthogonalize(double zero_tol = 1.0e-15); /** * @brief Orthonormalizes the matrix's last column, assuming the previous * columns are already orthonormal. * - * By default, the function considers the whole matrix. - * If input parameter ncols < d_num_cols, then a subset of the matrix + * By default, ncols == -1, and the function considers the whole matrix. + * If ncols != -1 and ncols < d_num_cols, then a subset of the matrix * is considered. * This allows one to reorthonormalize the matrix every time a new column * is added, assuming the previous columns have remained unchanged. + * + * If the norm of a matrix column is below the value of zero_tol then it + * is considered to be zero, and we do not divide by it. + * Therefore, that column is considered to be zero and is not normalized. + * By default, zero_tol = 1.0e-15. */ void - orthogonalize_last(int ncols = -1); + orthogonalize_last(int ncols = -1, double zero_tol = 1.0e-15); /** * @brief Rescale every matrix row by its maximum absolute value.