Commit 30dcd8e2 authored by Pavel Burovskiy's avatar Pavel Burovskiy
Browse files

Clean up to the code

parent 9371a8ef
///////////////////////////////////////////////////////////////////////////////
//
// File Vmath.hpp
// File: Vmath.hpp
//
// For more information, please see: http://www.nektar.info
//
......@@ -492,9 +492,9 @@ namespace Vmath
}
}
template LIB_UTILITIES_EXPORT void Svtvp(int n, const double alpha, const double *x,
const int incx, const double *y, const int incy,
double *z, const int incz);
template LIB_UTILITIES_EXPORT void Svtvp(int n, const Nektar::NekDouble alpha, const Nektar::NekDouble *x,
const int incx, const Nektar::NekDouble *y, const int incy,
Nektar::NekDouble *z, const int incz);
/// \brief svtvp (scalar times vector plus vector): z = alpha*x - y
......
......@@ -34,7 +34,6 @@
///////////////////////////////////////////////////////////////////////////////
#include <MultiRegions/GlobalLinSysIterative.h>
#include <LibUtilities/BasicUtils/Timer.h>
namespace Nektar
{
......@@ -362,7 +361,19 @@ namespace Nektar
/**  
* Solve a global linear system using the conjugate gradient method.  
* We solve only for the non-Dirichlet modes. The operator is evaluated  
* using an auxiliary function v_DoMatrixMultiply defined by the  
* specific solver. Distributed math routines are used to support  
* parallel execution of the solver.  
*  
* The implemented algorithm uses a reduced-communication reordering of  
* the standard PCG method (Demmel, Heath and Vorst, 1993)  
*  
* @param pInput Input residual of all DOFs.  
* @param pOutput Solution vector of all DOFs.  
*/  
void GlobalLinSysIterative::DoConjugateGradient(
const int nGlobal,
const Array<OneD,const NekDouble> &pInput,
......@@ -395,7 +406,6 @@ namespace Nektar
Array<OneD, NekDouble> r_A (nNonDir, 0.0);
Array<OneD, NekDouble> q_A (nNonDir, 0.0);
Array<OneD, NekDouble> tmp;
Array<OneD, NekDouble> tmp_A (nNonDir, 0.0);
// Create NekVector wrappers for linear algebra operations
NekVector<NekDouble> in (nNonDir,pInput + nDir, eWrapper);
......@@ -479,17 +489,13 @@ namespace Nektar
"Exceeded maximum number of iterations (5000)");
// Compute new search direction p_k, q_k
//p = w + beta * p;
//q = s + beta * q;
Vmath::Svtvp(nNonDir, beta, &p_A[0], 1, &w_A[nDir], 1, &p_A[0], 1);
Vmath::Svtvp(nNonDir, beta, &q_A[0], 1, &s_A[nDir], 1, &q_A[0], 1);
// Update solution x_{k+1}
//out = out + alpha * p;
Vmath::Svtvp(nNonDir, alpha, &p_A[0], 1, &pOutput[nDir], 1, &pOutput[nDir], 1);
// Update residual vector r_{k+1}
//r = r - alpha * q;
Vmath::Svtvp(nNonDir, -alpha, &q_A[0], 1, &r_A[0], 1, &r_A[0], 1);
// Apply preconditioner
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment