Skip to content

Changes in Gram-Schmidt functions #239

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
cval26 merged 7 commits into from
Dec 15, 2023
Merged

Changes in Gram-Schmidt functions #239

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
2ff3ee0
fix issues in Matrix::orthogonalize
cval26 Sep 27, 2023
698b157
add Matrix::orthogonalize_last method
cval26 Sep 27, 2023
d1c65a2
run astyle
cval26 Sep 27, 2023
9cbc970
add another test for Matrix::orthogonalize_last
cval26 Sep 27, 2023
6f1a113
rerun astyle properly
cval26 Sep 27, 2023
f4893d9
review n.1 changes
cval26 Sep 28, 2023
419ba38
review n.2
cval26 Nov 1, 2023
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
105 changes: 73 additions & 32 deletions lib/linalg/Matrix.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1788,47 +1788,88 @@ const
}

void
Matrix::orthogonalize()
Matrix::orthogonalize(double zero_tol)
{
for (int work = 1; work < d_num_cols; ++work) {
double tmp;
for (int col = 0; col < work; ++col) {
for (int work = 0; work < d_num_cols; ++work)
{
// Orthogonalize the column.
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,
MPI_COMM_WORLD);
}
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)
factor += item(i, col) * item(i, work);

if (d_distributed && d_num_procs > 1)
{
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);
}

// 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) {
MPI_Allreduce(&tmp, &norm, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);

for (int i = 0; i < d_num_rows; ++i)
norm += item(i, work) * item(i, work);

if (d_distributed && d_num_procs > 1)
{
CAROM_VERIFY( MPI_Allreduce(MPI_IN_PLACE, &norm, 1, MPI_DOUBLE,
MPI_SUM, MPI_COMM_WORLD) == MPI_SUCCESS );
}
else {
norm = tmp;
if (norm > zero_tol)
{
norm = 1.0 / sqrt(norm);
for (int i = 0; i < d_num_rows; ++i)
item(i, work) *= norm;
}
norm = sqrt(norm);
for (int i = 0; i < d_num_rows; ++i) {
item(i, work) /= norm;
}
}

void
Matrix::orthogonalize_last(int ncols, double zero_tol)
{
if (ncols == -1) ncols = d_num_cols;
CAROM_VERIFY((ncols > 0) && (ncols <= d_num_cols));

const int last_col = ncols - 1; // index of column to be orthonormalized

// Orthogonalize the column.
for (int col = 0; col < last_col; ++col)
{
double factor = 0.0;

for (int i = 0; i < d_num_rows; ++i)
factor += item(i, col) * item(i, last_col);

if (d_distributed && d_num_procs > 1)
{
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);
}

// Normalize the column.
double norm = 0.0;

for (int i = 0; i < d_num_rows; ++i)
norm += item(i, last_col) * item(i, last_col);

if (d_distributed && d_num_procs > 1)
{
CAROM_VERIFY( MPI_Allreduce(MPI_IN_PLACE, &norm, 1, MPI_DOUBLE,
MPI_SUM, MPI_COMM_WORLD) == MPI_SUCCESS );
}
if (norm > zero_tol)
{
norm = 1.0 / sqrt(norm);
for (int i = 0; i < d_num_rows; ++i)
item(i, last_col) *= norm;
}
}

Expand Down
27 changes: 25 additions & 2 deletions lib/linalg/Matrix.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -941,10 +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(double zero_tol = 1.0e-15);

/**
* @brief Orthonormalizes the matrix's last column, assuming the previous
* columns are already orthonormal.
*
* 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();
orthogonalize_last(int ncols = -1, double zero_tol = 1.0e-15);

/**
* @brief Rescale every matrix row by its maximum absolute value.
Expand Down
147 changes: 147 additions & 0 deletions unit_tests/test_Matrix.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -457,6 +457,153 @@ 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_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_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)
{
/**
Expand Down