From 884850bcbbdaf2f8f3fe67a4ad97bf8ec2dc4089 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 11 Jan 2017 17:23:10 +0000
Subject: [PATCH] ENH: GAMGSolver: use raw pointers in critical loops. Fixes
 #379.

---
 .../solvers/GAMG/GAMGSolverInterpolate.C      | 15 +++++++++-----
 .../lduMatrix/solvers/GAMG/GAMGSolverScale.C  | 20 +++++++++++++------
 2 files changed, 24 insertions(+), 11 deletions(-)

diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverInterpolate.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverInterpolate.C
index 022d91c74b..3896775f43 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverInterpolate.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverInterpolate.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -107,25 +107,30 @@ void Foam::GAMGSolver::interpolate
     const label nCells = m.diag().size();
     scalar* __restrict__ psiPtr = psi.begin();
     const scalar* const __restrict__ diagPtr = m.diag().begin();
+    const scalar* const __restrict__ psiCPtr = psiC.begin();
+
 
     const label nCCells = psiC.size();
     scalarField corrC(nCCells, 0);
+    scalar* __restrict__ corrCPtr = corrC.begin();
+
     scalarField diagC(nCCells, 0);
+    scalar* __restrict__ diagCPtr = diagC.begin();
 
     for (label celli=0; celli<nCells; celli++)
     {
-        corrC[restrictAddressing[celli]] += diagPtr[celli]*psiPtr[celli];
-        diagC[restrictAddressing[celli]] += diagPtr[celli];
+        corrCPtr[restrictAddressing[celli]] += diagPtr[celli]*psiPtr[celli];
+        diagCPtr[restrictAddressing[celli]] += diagPtr[celli];
     }
 
     for (label ccelli=0; ccelli<nCCells; ccelli++)
     {
-        corrC[ccelli] = psiC[ccelli] - corrC[ccelli]/diagC[ccelli];
+        corrCPtr[ccelli] = psiCPtr[ccelli] - corrCPtr[ccelli]/diagCPtr[ccelli];
     }
 
     for (label celli=0; celli<nCells; celli++)
     {
-        psiPtr[celli] += corrC[restrictAddressing[celli]];
+        psiPtr[celli] += corrCPtr[restrictAddressing[celli]];
     }
 }
 
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverScale.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverScale.C
index b81455ff4b..a587b667c7 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverScale.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverScale.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -48,13 +48,20 @@ void Foam::GAMGSolver::scale
         cmpt
     );
 
+
+    const label nCells = field.size();
+    scalar* __restrict__ fieldPtr = field.begin();
+    const scalar* const __restrict__ sourcePtr = source.begin();
+    const scalar* const __restrict__ AcfPtr = Acf.begin();
+
+
     scalar scalingFactorNum = 0.0;
     scalar scalingFactorDenom = 0.0;
 
-    forAll(field, i)
+    for (label i=0; i<nCells; i++)
     {
-        scalingFactorNum += source[i]*field[i];
-        scalingFactorDenom += Acf[i]*field[i];
+        scalingFactorNum += sourcePtr[i]*fieldPtr[i];
+        scalingFactorDenom += AcfPtr[i]*fieldPtr[i];
     }
 
     vector2D scalingVector(scalingFactorNum, scalingFactorDenom);
@@ -68,10 +75,11 @@ void Foam::GAMGSolver::scale
     }
 
     const scalarField& D = A.diag();
+    const scalar* const __restrict__ DPtr = D.begin();
 
-    forAll(field, i)
+    for (label i=0; i<nCells; i++)
     {
-        field[i] = sf*field[i] + (source[i] - sf*Acf[i])/D[i];
+        fieldPtr[i] = sf*fieldPtr[i] + (sourcePtr[i] - sf*AcfPtr[i])/DPtr[i];
     }
 }