Skip to content

[newchem-cpp] transcribe gaussj#393

Open
ChristopherBignamini wants to merge 20 commits into
grackle-project:newchem-cppfrom
ChristopherBignamini:newchem-cpp_transcribe_gaussj
Open

[newchem-cpp] transcribe gaussj#393
ChristopherBignamini wants to merge 20 commits into
grackle-project:newchem-cppfrom
ChristopherBignamini:newchem-cpp_transcribe_gaussj

Conversation

@ChristopherBignamini

@ChristopherBignamini ChristopherBignamini commented Aug 26, 2025

Copy link
Copy Markdown
Contributor

$\textcolor{red}{\textbf{MWA EDIT:}}$

This was originally proposed as brittonsmith#48


This should be reviewed after #386 is merged


This PR introduces a C++ function written "from-scratch" that is intended to replace the gaussj_g Fortran subroutine. The purpose was to avoid the licensing issues of the gaussj_g routine (since it had been ripped out of Numerical Recipes)

In the original PR, @ChristopherBignamini notes that "Implementation is quite basic and while it passes Grackle test it could have numerical issues in case of bad conditioned matrices."

@mabruzzo mabruzzo left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I took a quick look to investigate why CI is failing. (I haven't had a chance to look at everything in context yet, but let's get CI working first)

I made a few suggestions to address the issues:

  • you aren't including <cmath> (to take absolute values). Plus, I don't think you need to be including <numeric>
  • I suggested the use of std::fabs over std::abs (this isn't essential, but I generally think it is a little less confusing to people more familiar with C)

You are also going to need to adjust the loops. CI is complaining about mixing unsigned indices with signed loop bounds

After you do all that, you are going to need to format this file (this is why pre-commit.ci is failing). I would strongly recommend that you use pre-commit for this rather than directly invoking clang-format.1 The easiest thing to do is add a GitHub comment that just says:

pre-commit.ci autofix

The continuous integration will then push a commit to your branch to adjusti the formatting. (If you really want to run it locally, the documentation provides additional context)

Footnotes

  1. If you choose to manually invoke clang-format, you'll need to match the exact version of clang-format referenced in the pre-commit file.

Comment thread src/clib/gaussj_g.cpp Outdated
Comment thread src/clib/gaussj_g.cpp Outdated
Comment thread src/clib/gaussj_g.cpp Outdated
Comment thread src/clib/gaussj_g.cpp Outdated
@ChristopherBignamini

Copy link
Copy Markdown
Contributor Author

pre-commit.ci autofix

Comment thread src/clib/gaussj_g.hpp Outdated
Comment on lines +1 to +4
// See LICENSE file for license and copyright information

/// @file gaussj_g.hpp
/// @brief Implementation of Gauss-Jordan elimination for solving linear systems

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
// See LICENSE file for license and copyright information
/// @file gaussj_g.hpp
/// @brief Implementation of Gauss-Jordan elimination for solving linear systems
//===----------------------------------------------------------------------===//
//
// See the LICENSE file for license and copyright information
// SPDX-License-Identifier: NCSA AND BSD-3-Clause
//
//===----------------------------------------------------------------------===//
///
/// @file
/// Declares the function for Gauss-Jordan elimination
///
//===----------------------------------------------------------------------===//

Comment thread src/clib/gaussj_g.cpp
@@ -0,0 +1,77 @@
#include <cmath>

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
#include <cmath>
//===----------------------------------------------------------------------===//
//
// See the LICENSE file for license and copyright information
// SPDX-License-Identifier: NCSA AND BSD-3-Clause
//
//===----------------------------------------------------------------------===//
///
/// @file
/// Implementation of Gauss-Jordan elimination for solving linear systems
///
//===----------------------------------------------------------------------===//
#include <cmath>

@ChristopherBignamini

Copy link
Copy Markdown
Contributor Author

@mabruzzo @brittonsmith the 4x4 LinAlgCase was failing with an error of 4e-15, the same test now passes after increasing the tolerance to 5e-15. If you prefer, I can try to change the precision by using long double.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

[newchem-cpp] We need to replace the gaussj_g routine (for licensing reasons)

2 participants