Skip to content

Commit 40e43f4

Browse files
eorbay-metueorbay-metu
committed
Apply clang-format style to CSysSolve nested solver
1 parent c26fe01 commit 40e43f4

File tree

2 files changed

+95
-102
lines changed

2 files changed

+95
-102
lines changed

‎Common/include/linear_algebra/CSysSolve.hpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -361,8 +361,8 @@ class CSysSolve {
361361
* \param[in] config - Definition of the particular problem.
362362
*/
363363
unsigned long FGMRESandBCGSTAB2_LinSolver(const VectorType& b, VectorType& x, const ProductType& mat_vec,
364-
const PrecondType& precond, ScalarType tol, unsigned long m, ScalarType& residual,
365-
bool monitoring, const CConfig* config) const;
364+
const PrecondType& precond, ScalarType tol, unsigned long m,
365+
ScalarType& residual, bool monitoring, const CConfig* config) const;
366366

367367
/*!
368368
* \brief Generic smoother (modified Richardson iteration with preconditioner)

‎Common/src/linear_algebra/CSysSolve.cpp‎

Lines changed: 93 additions & 100 deletions
Original file line numberDiff line numberDiff line change
@@ -739,119 +739,115 @@ unsigned long CSysSolve<ScalarType>::BCGSTAB_LinSolver(const CSysVector<ScalarTy
739739

740740
template <class ScalarType>
741741
void BCGSTABpre_parallel(const CSysVector<ScalarType>& a, CSysVector<ScalarType>& b_in,
742-
const CMatrixVectorProduct<ScalarType>& mat_vec, const CPreconditioner<ScalarType>& precond, const CConfig* config) {
743-
744-
// NOTE: norm0_in is currently unused. To avoid -Werror unused-variable warnings,
745-
// we mark it [[maybe_unused]] so that future changes can reuse it without
746-
// triggering a build failure.
747-
ScalarType norm_r_in = 0.0;
748-
[[maybe_unused]] ScalarType norm0_in = 0.0;
749-
unsigned long ii = 0;
750-
751-
CSysVector<ScalarType> A_z_i;
752-
CSysVector<ScalarType> r_0_in;
753-
CSysVector<ScalarType> r_in;
754-
CSysVector<ScalarType> p;
755-
CSysVector<ScalarType> v;
756-
CSysVector<ScalarType> z_i;
757-
758-
/*--- Allocate if not allocated yet ---*/
759-
BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
760-
auto nVar = a.GetNVar();
761-
auto nBlk = a.GetNBlk();
762-
auto nBlkDomain = a.GetNBlkDomain();
763-
764-
A_z_i.Initialize(nBlk, nBlkDomain, nVar, nullptr);
765-
r_0_in.Initialize(nBlk, nBlkDomain, nVar, nullptr);
766-
r_in.Initialize(nBlk, nBlkDomain, nVar, nullptr);
767-
p.Initialize(nBlk, nBlkDomain, nVar, nullptr);
768-
v.Initialize(nBlk, nBlkDomain, nVar, nullptr);
769-
z_i.Initialize(nBlk, nBlkDomain, nVar, nullptr);
770-
}
771-
END_SU2_OMP_SAFE_GLOBAL_ACCESS
742+
const CMatrixVectorProduct<ScalarType>& mat_vec, const CPreconditioner<ScalarType>& precond,
743+
const CConfig* config) {
744+
// NOTE: norm0_in is currently unused. To avoid -Werror unused-variable warnings,
745+
// we mark it [[maybe_unused]] so that future changes can reuse it without
746+
// triggering a build failure.
747+
ScalarType norm_r_in = 0.0;
748+
[[maybe_unused]] ScalarType norm0_in = 0.0;
749+
unsigned long ii = 0;
750+
751+
CSysVector<ScalarType> A_z_i;
752+
CSysVector<ScalarType> r_0_in;
753+
CSysVector<ScalarType> r_in;
754+
CSysVector<ScalarType> p;
755+
CSysVector<ScalarType> v;
756+
CSysVector<ScalarType> z_i;
772757

773-
/*--- Calculate the initial residual, compute norm, and check if system is already solved ---*/
774-
mat_vec(b_in, A_z_i);
775-
r_in = a - A_z_i;
758+
/*--- Allocate if not allocated yet ---*/
759+
BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
760+
auto nVar = a.GetNVar();
761+
auto nBlk = a.GetNBlk();
762+
auto nBlkDomain = a.GetNBlkDomain();
763+
764+
A_z_i.Initialize(nBlk, nBlkDomain, nVar, nullptr);
765+
r_0_in.Initialize(nBlk, nBlkDomain, nVar, nullptr);
766+
r_in.Initialize(nBlk, nBlkDomain, nVar, nullptr);
767+
p.Initialize(nBlk, nBlkDomain, nVar, nullptr);
768+
v.Initialize(nBlk, nBlkDomain, nVar, nullptr);
769+
z_i.Initialize(nBlk, nBlkDomain, nVar, nullptr);
770+
}
771+
END_SU2_OMP_SAFE_GLOBAL_ACCESS
776772

777-
/*--- Only compute the residuals in full communication mode. ---*/
778-
if (config->GetComm_Level() == COMM_FULL) {
779-
norm_r_in = r_in.norm();
780-
norm0_in = a.norm();
781-
/*--- Set the norm to the initial residual value ---*/
782-
/*--- if (tol_type == LinearToleranceType::RELATIVE) norm0_in = norm_r_in; ---*/
783-
}
773+
/*--- Calculate the initial residual, compute norm, and check if system is already solved ---*/
774+
mat_vec(b_in, A_z_i);
775+
r_in = a - A_z_i;
784776

785-
/*--- Initialization ---*/
786-
ScalarType alpha = 1.0, omega = 1.0, rho = 1.0, rho_prime = 1.0;
787-
p = ScalarType(0.0);
788-
v = ScalarType(0.0);
789-
r_0_in = r_in;
777+
/*--- Only compute the residuals in full communication mode. ---*/
778+
if (config->GetComm_Level() == COMM_FULL) {
779+
norm_r_in = r_in.norm();
780+
norm0_in = a.norm();
781+
/*--- Set the norm to the initial residual value ---*/
782+
/*--- if (tol_type == LinearToleranceType::RELATIVE) norm0_in = norm_r_in; ---*/
783+
}
790784

791-
ScalarType tolerance = 1e-5; // Tolerance for the residual norm
785+
/*--- Initialization ---*/
786+
ScalarType alpha = 1.0, omega = 1.0, rho = 1.0, rho_prime = 1.0;
787+
p = ScalarType(0.0);
788+
v = ScalarType(0.0);
789+
r_0_in = r_in;
792790

793-
/*--- Loop over all search directions ---*/
794-
while (ii < 1000) { // Arbitrary high iteration limit for safety
795-
/*--- Compute rho_prime ---*/
796-
rho_prime = rho;
791+
ScalarType tolerance = 1e-5; // Tolerance for the residual norm
797792

798-
/*--- Compute rho_i ---*/
799-
rho = r_in.dot(r_0_in);
793+
/*--- Loop over all search directions ---*/
794+
while (ii < 1000) { // Arbitrary high iteration limit for safety
795+
/*--- Compute rho_prime ---*/
796+
rho_prime = rho;
800797

801-
/*--- Compute beta ---*/
802-
ScalarType beta_in = (rho / rho_prime) * (alpha / omega);
798+
/*--- Compute rho_i ---*/
799+
rho = r_in.dot(r_0_in);
803800

804-
/*--- Update p ---*/
805-
p = beta_in * (p - omega * v) + r_in;
801+
/*--- Compute beta ---*/
802+
ScalarType beta_in = (rho / rho_prime) * (alpha / omega);
806803

807-
/*--- Preconditioning step ---*/
808-
precond(p, z_i);
809-
mat_vec(z_i, v);
804+
/*--- Update p ---*/
805+
p = beta_in * (p - omega * v) + r_in;
810806

811-
/*--- Calculate step-length alpha ---*/
812-
ScalarType r_0_v = r_0_in.dot(v);
813-
alpha = rho / r_0_v;
807+
/*--- Preconditioning step ---*/
808+
precond(p, z_i);
809+
mat_vec(z_i, v);
814810

815-
/*--- Update solution and residual ---*/
816-
b_in += alpha * z_i;
817-
r_in -= alpha * v;
811+
/*--- Calculate step-length alpha ---*/
812+
ScalarType r_0_v = r_0_in.dot(v);
813+
alpha = rho / r_0_v;
818814

819-
/*--- Preconditioning step ---*/
820-
precond(r_in, z_i);
821-
mat_vec(z_i, A_z_i);
815+
/*--- Update solution and residual ---*/
816+
b_in += alpha * z_i;
817+
r_in -= alpha * v;
822818

823-
/*--- Calculate step-length omega, avoid division by 0. ---*/
824-
omega = A_z_i.squaredNorm();
825-
if (omega == ScalarType(0)) break;
826-
omega = A_z_i.dot(r_in) / omega;
819+
/*--- Preconditioning step ---*/
820+
precond(r_in, z_i);
821+
mat_vec(z_i, A_z_i);
827822

828-
/*--- Update solution and residual ---*/
829-
b_in += omega * z_i;
830-
r_in -= omega * A_z_i;
823+
/*--- Calculate step-length omega, avoid division by 0. ---*/
824+
omega = A_z_i.squaredNorm();
825+
if (omega == ScalarType(0)) break;
826+
omega = A_z_i.dot(r_in) / omega;
831827

832-
/*--- Update the residual norm ---*/
833-
norm_r_in = r_in.norm();
828+
/*--- Update solution and residual ---*/
829+
b_in += omega * z_i;
830+
r_in -= omega * A_z_i;
834831

835-
/*--- Check if residual norm is below tolerance ---*/
836-
if (norm_r_in < tolerance) {
837-
break; // Stop if the residual norm is below the desired tolerance
838-
}
832+
/*--- Update the residual norm ---*/
833+
norm_r_in = r_in.norm();
839834

840-
ii++; // Increment iteration counter
835+
/*--- Check if residual norm is below tolerance ---*/
836+
if (norm_r_in < tolerance) {
837+
break; // Stop if the residual norm is below the desired tolerance
841838
}
842-
}
843-
844-
845839

840+
ii++; // Increment iteration counter
841+
}
842+
}
846843

847844
template <class ScalarType>
848-
unsigned long CSysSolve<ScalarType>::FGMRESandBCGSTAB2_LinSolver(const CSysVector<ScalarType>& b, CSysVector<ScalarType>& x,
849-
const CMatrixVectorProduct<ScalarType>& mat_vec,
850-
const CPreconditioner<ScalarType>& precond, ScalarType tol,
851-
unsigned long m, ScalarType& residual, bool monitoring,
852-
const CConfig* config) const {
853-
854-
845+
unsigned long CSysSolve<ScalarType>::FGMRESandBCGSTAB2_LinSolver(const CSysVector<ScalarType>& b,
846+
CSysVector<ScalarType>& x,
847+
const CMatrixVectorProduct<ScalarType>& mat_vec,
848+
const CPreconditioner<ScalarType>& precond,
849+
ScalarType tol, unsigned long m, ScalarType& residual,
850+
bool monitoring, const CConfig* config) const {
855851
const bool masterRank = (SU2_MPI::GetRank() == MASTER_NODE);
856852
const bool flexible = !precond.IsIdentity();
857853
/*--- If we call the solver outside of a parallel region, but the number of threads allows,
@@ -939,7 +935,7 @@ unsigned long CSysSolve<ScalarType>::FGMRESandBCGSTAB2_LinSolver(const CSysVecto
939935

940936
if (flexible) {
941937
/*--- Use BCGSTAB as inner iteration ---*/
942-
BCGSTABpre_parallel(W[i], Z[i], mat_vec, precond, config);
938+
BCGSTABpre_parallel(W[i], Z[i], mat_vec, precond, config);
943939

944940
mat_vec(Z[i], W[i + 1]);
945941
} else {
@@ -1005,9 +1001,6 @@ unsigned long CSysSolve<ScalarType>::FGMRESandBCGSTAB2_LinSolver(const CSysVecto
10051001
return i;
10061002
}
10071003

1008-
1009-
1010-
10111004
template <class ScalarType>
10121005
unsigned long CSysSolve<ScalarType>::Smoother_LinSolver(const CSysVector<ScalarType>& b, CSysVector<ScalarType>& x,
10131006
const CMatrixVectorProduct<ScalarType>& mat_vec,
@@ -1245,8 +1238,8 @@ unsigned long CSysSolve<ScalarType>::Solve(CSysMatrix<ScalarType>& Jacobian, con
12451238
ScreenOutput, config);
12461239
break;
12471240
case FGMRESandBCGSTAB2:
1248-
IterLinSol = FGMRESandBCGSTAB2_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, residual,
1249-
ScreenOutput, config);
1241+
IterLinSol = FGMRESandBCGSTAB2_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter,
1242+
residual, ScreenOutput, config);
12501243
break;
12511244
case CONJUGATE_GRADIENT:
12521245
IterLinSol = CG_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, residual,
@@ -1408,8 +1401,8 @@ unsigned long CSysSolve<ScalarType>::Solve_b(CSysMatrix<ScalarType>& Jacobian, c
14081401
ScreenOutput, config);
14091402
break;
14101403
case FGMRESandBCGSTAB2:
1411-
IterLinSol = FGMRESandBCGSTAB2_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, residual,
1412-
ScreenOutput, config);
1404+
IterLinSol = FGMRESandBCGSTAB2_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter,
1405+
residual, ScreenOutput, config);
14131406
break;
14141407
case CONJUGATE_GRADIENT:
14151408
IterLinSol = CG_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, residual,

0 commit comments

Comments
 (0)