Commit c12b0358 authored by Dave Moxey's avatar Dave Moxey
Browse files

Merge branch 'fix/blas-calls-in-SL-staticcond' of localhost:nektar

parents aaecbf1c 58010c0f
......@@ -110,6 +110,7 @@ ADD_NEKTAR_TEST(Helmholtz2D_HDG_P7_Modes_AllBCs)
ADD_NEKTAR_TEST(Helmholtz3D_CG_Hex)
ADD_NEKTAR_TEST(Helmholtz3D_CG_Hex_AllBCs)
ADD_NEKTAR_TEST(Helmholtz3D_CG_Hex_AllBCs_iter_ml)
ADD_NEKTAR_TEST(Helmholtz3D_CG_Hex_AllBCs_iter_sc_cont)
ADD_NEKTAR_TEST(Helmholtz3D_CG_Tet)
ADD_NEKTAR_TEST(Helmholtz3D_CG_Prism)
ADD_NEKTAR_TEST(Helmholtz3D_CG_Prism_iter_ml)
......
<?xml version="1.0" encoding="utf-8"?>
<test>
<description>Helmholtz 3D CG, hexes, mixed BCs, iterative SC, contiguous storage strategy</description>
<executable>Helmholtz3D</executable>
<parameters>-I GlobalSysSoln=IterativeStaticCond -I LocalMatrixStorageStrategy=Contiguous Helmholtz3D_Hex_AllBCs_P6.xml</parameters>
<files>
<file description="Session File">Helmholtz3D_Hex_AllBCs_P6.xml</file>
</files>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-8">0.000416575</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-8">0.000871589</value>
</metric>
</metrics>
</test>
///////////////////////////////////////////////////////////////////////////////
//
// File GlobalLinSysIterativeStaticCond.cpp
// File: GlobalLinSysIterativeStaticCond.cpp
//
// For more information, please see: http://www.nektar.info
//
......@@ -62,6 +62,22 @@ namespace Nektar
GlobalLinSysIterativeStaticCond::create,
"Iterative multi-level static condensation.");
std::string GlobalLinSysIterativeStaticCond::storagedef =
LibUtilities::SessionReader::RegisterDefaultSolverInfo(
"LocalMatrixStorageStrategy",
"Non-contiguous");
std::string GlobalLinSysIterativeStaticCond::storagelookupIds[2] = {
LibUtilities::SessionReader::RegisterEnumValue(
"LocalMatrixStorageStrategy",
"Contiguous",
MultiRegions::eContiguous),
LibUtilities::SessionReader::RegisterEnumValue(
"LocalMatrixStorageStrategy",
"Non-contiguous",
MultiRegions::eNonContiguous),
};
/**
* For a matrix system of the form @f[
* \left[ \begin{array}{cc}
......@@ -167,8 +183,8 @@ namespace Nektar
int nLocBndDofs = pLocToGloMap->GetNumLocalBndCoeffs();
int nIntDofs = pLocToGloMap->GetNumGlobalCoeffs()
- nGlobBndDofs;
Array<OneD, NekDouble> F = m_wsp + nLocBndDofs;
Array<OneD, NekDouble> F = m_wsp + 2*nLocBndDofs;
Array<OneD, NekDouble> tmp;
if(nDirBndDofs && dirForcCalculated)
{
......@@ -367,8 +383,8 @@ namespace Nektar
{
int nLocalBnd = m_locToGloMap->GetNumLocalBndCoeffs();
int nGlobal = m_locToGloMap->GetNumGlobalCoeffs();
m_wsp = Array<OneD, NekDouble>(nLocalBnd + nGlobal);
m_wsp = Array<OneD, NekDouble>(2*nLocalBnd + nGlobal);
if(pLocToGloMap->AtLastLevel())
{
// decide whether to assemble schur complement globally
......@@ -381,22 +397,61 @@ namespace Nektar
{
AssembleSchurComplement(pLocToGloMap);
}
int nbdry, nblks;
unsigned int esize[1];
int nBlk = m_schurCompl->GetNumberOfBlockRows();
m_schurComplBlock = Array<OneD, DNekScalBlkMatSharedPtr>(nBlk);
for (int i = 0; i < nBlk; ++i)
else
{
nbdry = m_schurCompl->GetBlock(i,i)->GetRows();
nblks = 1;
esize[0] = nbdry;
m_schurComplBlock[i] = MemoryManager<DNekScalBlkMat>
::AllocateSharedPtr(nblks, nblks, esize, esize);
m_schurComplBlock[i]->SetBlock(
0, 0, m_schurCompl->GetBlock(i,i));
}
LocalMatrixStorageStrategy storageStrategy = m_expList.lock()->GetSession()->
GetSolverInfoAsEnum<LocalMatrixStorageStrategy>("LocalMatrixStorageStrategy");
size_t storageSize = 0;
int nBlk = m_schurCompl->GetNumberOfBlockRows();
m_scale = Array<OneD, NekDouble> (nBlk, 1.0);
m_rows = Array<OneD, unsigned int> (nBlk, 0U);
// Determine storage requirements for dense blocks.
for (int i = 0; i < nBlk; ++i)
{
m_rows[i] = m_schurCompl->GetBlock(i,i)->GetRows();
m_scale[i] = m_schurCompl->GetBlock(i,i)->Scale();
storageSize += m_rows[i] * m_rows[i];
}
// Assemble dense storage blocks.
DNekScalMatSharedPtr loc_mat;
m_denseBlocks.resize(nBlk);
double *ptr = 0;
if (MultiRegions::eContiguous == storageStrategy)
{
m_storage.resize (storageSize);
ptr = &m_storage[0];
}
for (unsigned int n = 0; n < nBlk; ++n)
{
loc_mat = m_schurCompl->GetBlock(n,n);
if (MultiRegions::eContiguous == storageStrategy)
{
int loc_lda = loc_mat->GetRows();
int blockSize = loc_lda * loc_lda;
m_denseBlocks[n] = ptr;
for(int i = 0; i < loc_lda; ++i)
{
for(int j = 0; j < loc_lda; ++j)
{
ptr[j*loc_lda+i] = (*loc_mat)(i,j);
}
}
ptr += blockSize;
}
else
{
m_denseBlocks[n] = loc_mat->GetRawPtr();
}
}
} // if (doGlobalOp)
}
else
{
......@@ -404,16 +459,26 @@ namespace Nektar
pLocToGloMap->GetNextLevelLocalToGlobalMap());
}
}
int GlobalLinSysIterativeStaticCond::v_GetNumBlocks()
{
return m_schurCompl->GetNumberOfBlockRows();
}
DNekScalBlkMatSharedPtr GlobalLinSysIterativeStaticCond::
v_GetStaticCondBlock(unsigned int n)
{
return m_schurComplBlock[n];
DNekScalBlkMatSharedPtr schurComplBlock;
DNekScalMatSharedPtr localMat = m_schurCompl->GetBlock(n,n);
int nbdry = localMat->GetRows();
int nblks = 1;
unsigned int esize[1] = {nbdry};
schurComplBlock = MemoryManager<DNekScalBlkMat>
::AllocateSharedPtr(nblks, nblks, esize, esize);
schurComplBlock->SetBlock(0, 0, localMat);
return schurComplBlock;
}
/**
......@@ -538,7 +603,7 @@ namespace Nektar
m_S1Blk->SetBlock(n,n, tmp_mat = loc_mat->GetBlock(0,0));
m_RBlk->SetBlock(n,n, tmp_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,m_R));
m_RTBlk->SetBlock(n,n, tmp_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,m_RT));
}
}
}
/**
......@@ -864,13 +929,19 @@ namespace Nektar
}
else
{
// Do matrix multiply locally
// Do matrix multiply locally, using direct BLAS calls
m_locToGloMap->GlobalToLocalBnd(pInput, m_wsp);
NekVector<NekDouble> loc(nLocal, m_wsp, eWrapper);
loc = (*m_schurCompl)*loc;
m_locToGloMap->AssembleBnd(m_wsp, pOutput);
int i, cnt;
Array<OneD, NekDouble> tmpout = m_wsp + nLocal;
for (i = cnt = 0; i < m_denseBlocks.size(); cnt += m_rows[i], ++i)
{
const int rows = m_rows[i];
Blas::Dgemv('N', rows, rows,
m_scale[i], m_denseBlocks[i], rows,
m_wsp.get()+cnt, 1,
0.0, tmpout.get()+cnt, 1);
}
m_locToGloMap->AssembleBnd(tmpout, pOutput);
}
}
......
///////////////////////////////////////////////////////////////////////////////
//
// File GlobalLinSysIterativeStaticCond.h
// File: GlobalLinSysIterativeStaticCond.h
//
// For more information, please see: http://www.nektar.info
//
......@@ -49,6 +49,19 @@ namespace Nektar
typedef boost::shared_ptr<GlobalLinSysIterativeStaticCond>
GlobalLinSysIterativeStaticCondSharedPtr;
enum LocalMatrixStorageStrategy
{
eContiguous,
eNonContiguous,
};
const char* const LocalMatrixStorageStrategyMap[] =
{
"Contiguous",
"Non-contiguous",
};
/// A global linear system.
class GlobalLinSysIterativeStaticCond : public GlobalLinSysIterative
{
......@@ -97,13 +110,21 @@ namespace Nektar
GlobalLinSysIterativeStaticCondSharedPtr m_recursiveSchurCompl;
/// Block Schur complement matrix.
DNekScalBlkMatSharedPtr m_schurCompl;
/// Dense storage for block Schur complement matrix.
std::vector<double> m_storage;
/// Vector of pointers to local matrix data
std::vector<const double*> m_denseBlocks;
/// Ranks of local matrices
Array<OneD, unsigned int> m_rows;
/// Scaling factors for local matrices
Array<OneD, NekDouble> m_scale;
/// Block \f$ BD^{-1} \f$ matrix.
DNekScalBlkMatSharedPtr m_BinvD;
/// Block \f$ C \f$ matrix.
DNekScalBlkMatSharedPtr m_C;
/// Block \f$ D^{-1} \f$ matrix.
DNekScalBlkMatSharedPtr m_invD;
// Block matrices for low energy
/// Block matrices for low energy
DNekScalBlkMatSharedPtr m_RBlk;
DNekScalBlkMatSharedPtr m_RTBlk;
DNekScalBlkMatSharedPtr m_S1Blk;
......@@ -115,8 +136,11 @@ namespace Nektar
Array<OneD, NekDouble> m_wsp;
/// Preconditioner object.
PreconditionerSharedPtr m_precon;
/// Wrapper for block matrices.
Array<OneD, DNekScalBlkMatSharedPtr> m_schurComplBlock;
/// Utility strings
static std::string storagedef;
static std::string storagelookupIds[];
/// Solve the linear system for given input and output vectors
/// using a specified local to global map.
......
Supports Markdown
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