Commit ee6dc08d authored by Pavel Burovskiy's avatar Pavel Burovskiy
Browse files

Adding code to test different Svtvp implementations

parent 7047e194
......@@ -28,7 +28,32 @@ SET(NodalTetElecSources
SET(TimeIntegrationDemoSources
TimeIntegrationDemo.cpp)
SET(SvtvpHeaderInlineSources
Svtvp-header-inline.hpp
Svtvp-header-inline.cpp)
SET(SvtvpHeaderSources
Svtvp-header.hpp
Svtvp-header.cpp)
SET(SvtvpExplicitInstLibSources
Svtvp-explicit-inst-inc.cpp
Svtvp-explicit-inst.cpp)
SET(SvtvpExplicitInstLibHeaders
Svtvp-explicit-inst-inc.hpp
Svtvp-explicit-inst.hpp)
SET(SvtvpExplicitInstIncSources
Svtvp-explicit-inst-inc-main.cpp)
SET(SvtvpExplicitInstSources
Svtvp-explicit-inst-main.cpp)
ADD_NEKTAR_LIBRARY(SvtvpExplicitInstLib lib ${NEKTAR_LIBRARY_TYPE} ${SvtvpExplicitInstLibSources} ${SvtvpExplicitInstLibHeaders})
#ADD_NEKTAR_EXECUTABLE(Graph demos GraphSources )
......@@ -58,3 +83,13 @@ SET_LAPACK_LINK_LIBRARIES(NodalTetElecDemo)
ADD_NEKTAR_EXECUTABLE(TimeIntegrationDemo demos TimeIntegrationDemoSources)
SET_LAPACK_LINK_LIBRARIES(TimeIntegrationDemo)
ADD_NEKTAR_EXECUTABLE(SvtvpHeaderInlineDemo demos SvtvpHeaderInlineSources)
ADD_NEKTAR_EXECUTABLE(SvtvpHeaderDemo demos SvtvpHeaderSources)
ADD_NEKTAR_EXECUTABLE(SvtvpExplicitInstIncDemo demos SvtvpExplicitInstIncSources)
ADD_NEKTAR_EXECUTABLE(SvtvpExplicitInstDemo demos SvtvpExplicitInstSources)
TARGET_LINK_LIBRARIES(SvtvpExplicitInstIncDemo SvtvpExplicitInstLib)
TARGET_LINK_LIBRARIES(SvtvpExplicitInstDemo SvtvpExplicitInstLib)
#include"Svtvp-explicit-inst-inc.hpp"
#include<vector>
#include<iostream>
#include<cstdlib>
using namespace std;
double experiment(int size)
{
vector<double> a(size, 1.0/rand());
vector<double> b(size, 2.0/rand());
vector<double> c(size, 3.0/rand());
SvtvpExplInstInc::Svtvp<double>(size, 5.0, &a[0], 1, &b[0], 1, &c[0], 1);
return c[size-1];
}
int main()
{
int n = 1000000;
for (int i = 0; i < n; i++)
{
experiment(10000);
}
}
#include"Svtvp-explicit-inst-inc.hpp"
namespace SvtvpExplInstInc
{
template<class T> void Svtvp(int n, const T alpha, const T *x,
const int incx, const T *y, const int incy,
T *z, const int incz)
{
++n;
if (incx == 1 && incy == 1 && incz == 1)
{
while( --n )
{
*z = alpha * (*x) + (*y);
++x;
++y;
++z;
}
}
else
{
while( --n )
{
*z = alpha * (*x) + (*y);
x += incx;
y += incy;
z += incz;
}
}
}
// explicit instantiation
template void Svtvp(int n, const double alpha, const double *x,
const int incx, const double *y, const int incy,
double *z, const int incz);
}
\ No newline at end of file
namespace SvtvpExplInstInc
{
/// \brief svtvp (scalar times vector plus vector): z = alpha*x + y
template<class T> void Svtvp(int n, const T alpha, const T *x,
const int incx, const T *y, const int incy,
T *z, const int incz);
}
\ No newline at end of file
#include"Svtvp-explicit-inst.hpp"
#include<vector>
#include<iostream>
#include<cstdlib>
using namespace std;
double experiment(int size)
{
vector<double> a(size, 1.0/rand());
vector<double> b(size, 2.0/rand());
vector<double> c(size, 3.0/rand());
SvtvpExplInst::Svtvp<double>(size, 5.0, &a[0], &b[0], &c[0]);
return c[size-1];
}
int main()
{
int n = 1000000;
for (int i = 0; i < n; i++)
{
experiment(10000);
}
}
#include"Svtvp-explicit-inst.hpp"
namespace SvtvpExplInst
{
template<class T> void Svtvp(int n, const T alpha, const T *x, const T *y, T *z)
{
++n;
while( --n )
{
*z = alpha * (*x) + (*y);
++x;
++y;
++z;
}
}
// explicit instantiation
template void Svtvp(int n, const double alpha,
const double *x,
const double *y,
double *z);
}
\ No newline at end of file
namespace SvtvpExplInst
{
/// \brief svtvp (scalar times vector plus vector): z = alpha*x + y
template<class T> void Svtvp(int n, const T alpha, const T *x, const T *y, T *z);
}
\ No newline at end of file
#include"Svtvp-header-inline.hpp"
#include<vector>
#include<iostream>
#include<cstdlib>
using namespace std;
double experiment(int size)
{
vector<double> a(size, 1.0/rand());
vector<double> b(size, 2.0/rand());
vector<double> c(size, 3.0/rand());
SvtvpHeaderInline::Svtvp<double>(size, 5.0, &a[0], 1, &b[0], 1, &c[0], 1);
return c[size-1];
}
int main()
{
int n = 1000000;
for (int i = 0; i < n; i++)
{
experiment(10000);
}
}
namespace SvtvpHeaderInline
{
/// \brief svtvp (scalar times vector plus vector): z = alpha*x + y
template<class T> inline void Svtvp(int n, const T alpha, const T *x,
const int incx, const T *y, const int incy,
T *z, const int incz)
{
++n;
if (incx == 1 && incy == 1 && incz == 1)
{
while( --n )
{
*z = alpha * (*x) + (*y);
++x;
++y;
++z;
}
}
else
{
while( --n )
{
*z = alpha * (*x) + (*y);
x += incx;
y += incy;
z += incz;
}
}
}
}
\ No newline at end of file
#include"Svtvp-header.hpp"
#include<LibUtilities/BasicUtils/Vmath.hpp>
#include<vector>
#include<iostream>
#include<cstdlib>
using namespace std;
double experiment(int size)
{
vector<double> a(size, 1.0/rand());
vector<double> b(size, 2.0/rand());
vector<double> c(size, 3.0/rand());
SvtvpHeader::Svtvp<double>(size, 5.0, &a[0], 1, &b[0], 1, &c[0], 1);
Vmath::Svtvp(size, 5.0, &a[0], 1, &b[0], 1, &c[0], 1);
return c[size-1];
}
double experimentVmath(int size)
{
vector<double> a(size, 1.0/rand());
vector<double> b(size, 2.0/rand());
vector<double> c(size, 3.0/rand());
Vmath::Svtvp(size, 5.0, &a[0], 1, &b[0], 1, &c[0], 1);
return c[size-1];
}
int main()
{
int n = 1000000;
for (int i = 0; i < n; i++)
{
experiment(10000);
}
}
namespace SvtvpHeader
{
/// \brief svtvp (scalar times vector plus vector): z = alpha*x + y
template<class T> void Svtvp(int n, const T alpha, const T *x,
const int incx, const T *y, const int incy,
T *z, const int incz)
{
++n;
if (incx == 1 && incy == 1 && incz == 1)
{
while( --n )
{
*z = alpha * (*x) + (*y);
++x;
++y;
++z;
}
}
else
{
while( --n )
{
*z = alpha * (*x) + (*y);
x += incx;
y += incy;
z += incz;
}
}
}
}
\ No newline at end of file
......@@ -463,8 +463,9 @@ namespace Vmath
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
template<class T> void Svtvp(int n, const T alpha, const T *x,
template<class T> LIB_UTILITIES_EXPORT void Svtvp(int n, const T alpha, const T *x,
const int incx, const T *y, const int incy,
T *z, const int incz)
{
......@@ -491,9 +492,26 @@ namespace Vmath
}
}
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
template<class T> LIB_UTILITIES_EXPORT void Svtvp(int n, const T alpha, const T *x, const T *y, T *z)
{
++n;
while( --n )
{
*z = alpha * (*x) + (*y);
++x;
++y;
++z;
}
}
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 double alpha, const double *x,
const double *y, double *z);
/// \brief svtvp (scalar times vector plus vector): z = alpha*x - y
template<class T> void Svtvm(int n, const T alpha, const T *x,
......
///////////////////////////////////////////////////////////////////////////////
//
// File Vmath.hpp
// File: Vmath.hpp
//
// For more information, please see: http://www.nektar.info
//
......@@ -153,6 +153,35 @@ namespace Vmath
template<class T> LIB_UTILITIES_EXPORT void Svtvp(int n, const T alpha, const T *x,
const int incx, const T *y, const int incy,
T *z, const int incz);
template<class T> LIB_UTILITIES_EXPORT void Svtvp(int n, const T alpha, const T *x, const T *y, T *z);
/// \brief svtvp (scalar times vector plus vector): z = alpha*x + y
template<class T> LIB_UTILITIES_EXPORT void Svtvph(int n, const T alpha, const T *x,
const int incx, const T *y, const int incy,
T *z, const int incz)
{
++n;
if (incx == 1 && incy == 1 && incz == 1)
{
while( --n )
{
*z = alpha * (*x) + (*y);
++x;
++y;
++z;
}
}
else
{
while( --n )
{
*z = alpha * (*x) + (*y);
x += incx;
y += incy;
z += incz;
}
}
}
/// \brief svtvp (scalar times vector plus vector): z = alpha*x - y
template<class T> LIB_UTILITIES_EXPORT void Svtvm(int n, const T alpha, const T *x,
......
......@@ -34,6 +34,7 @@
///////////////////////////////////////////////////////////////////////////////
#include <MultiRegions/GlobalLinSysIterative.h>
#include <LibUtilities/BasicUtils/Timer.h>
namespace Nektar
{
......@@ -361,19 +362,7 @@ 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,
......@@ -406,6 +395,7 @@ 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);
......@@ -463,6 +453,36 @@ namespace Nektar
v_DoMatrixMultiply(w_A, s_A);
k = 0;
Timer t;
int rank = vComm->GetRank();
cout << "MPI " << rank << ": calling Svtvp 1000 times on array with " << nNonDir << " elements" << endl;
t.Start();
for (int i = 0 ; i < 1000; i++)
{
Vmath::Svtvp(nNonDir, beta, &p_A[0], 1, &w_A[nDir], 1, &tmp_A[0], 1);
}
t.Stop();
cout << "MPI " << rank << ": One call to cpp Svtvp with increment takes " << t.TimePerTest(1000) << endl;
t.Start();
for (int i = 0 ; i < 1000; i++)
{
Vmath::Svtvph(nNonDir, beta, &p_A[0], 1, &w_A[nDir], 1, &tmp_A[0], 1);
}
t.Stop();
cout << "MPI " << rank << ": One call to hpp Svtvp with increment takes " << t.TimePerTest(1000) << endl;
t.Start();
for (int i = 0 ; i < 1000; i++)
{
Vmath::Svtvp(nNonDir, beta, &p_A[0], &w_A[nDir], &tmp_A[0]);
}
t.Stop();
cout << "MPI " << rank << ": One call to cpp Svtvp without increment takes " << t.TimePerTest(1000) << endl;
vExchange[0] = Vmath::Dot2(nNonDir,
r_A,
w_A + nDir,
......@@ -493,14 +513,18 @@ namespace Nektar
//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);
//Vmath::Svtvp(nNonDir, beta, &p_A[0], &w_A[nDir], &p_A[0]);
//Vmath::Svtvp(nNonDir, beta, &q_A[0], &s_A[nDir], &q_A[0]);
// Update solution x_{k+1}
//out = out + alpha * p;
Vmath::Svtvp(nNonDir, alpha, &p_A[0], 1, &pOutput[nDir], 1, &pOutput[nDir], 1);
//Vmath::Svtvp(nNonDir, alpha, &p_A[0], &pOutput[nDir], &pOutput[nDir]);
// 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);
//Vmath::Svtvp(nNonDir, -alpha, &q_A[0], &r_A[0], &r_A[0]);
// Apply preconditioner
m_precon->DoPreconditioner(r_A, tmp = w_A + nDir);
......@@ -558,7 +582,6 @@ namespace Nektar
}
}
void GlobalLinSysIterative::Set_Rhs_Magnitude(const NekVector<NekDouble> &pIn)
{
......
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