Skip to content

Commit

Permalink
Adding iterative refinement implementation into [Dist]SparseLDLFactor…
Browse files Browse the repository at this point in the history 
…ization
  • Loading branch information
@poulson
poulson committed Dec 28, 2016
1 parent f66fe79 commit 3ea8fc2
Show file tree
Hide file tree
Showing 17 changed files with 356 additions and 1,160 deletions.
93 changes: 8 additions & 85 deletions include/El/lapack_like/factor.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -252,19 +252,6 @@ void SolveAfter
( const DistNodeInfo& info,
const DistFront<Field>& front, DistMatrixNode<Field>& X );

template<typename Field>
Int SolveWithIterativeRefinement
( const SparseMatrix<Field>& A,
const vector<Int>& invMap, const NodeInfo& info,
const Front<Field>& front, Matrix<Field>& y,
Base<Field> relTolRefine, Int maxRefineIts );
template<typename Field>
Int SolveWithIterativeRefinement
( const DistSparseMatrix<Field>& A,
const DistMap& invMap, const DistNodeInfo& info,
const DistFront<Field>& front, DistMultiVec<Field>& y,
Base<Field> relTolRefine, Int maxRefineIts );

// Solve linear system with the implicit representations of L, D, and P
// --------------------------------------------------------------------
template<typename Field>
Expand Down Expand Up @@ -318,9 +305,7 @@ template<typename Field>
Int RegularizedSolveAfter
( const SparseMatrix<Field>& A,
const Matrix<Base<Field>>& reg,
const vector<Int>& invMap,
const ldl::NodeInfo& info,
const ldl::Front<Field>& front,
const SparseLDLFactorization<Field>& sparseLDLFact,
Matrix<Field>& B,
Base<Field> relTolRefine,
Int maxRefineIts,
Expand All @@ -330,36 +315,19 @@ template<typename Field>
Int RegularizedSolveAfter
( const DistSparseMatrix<Field>& A,
const DistMultiVec<Base<Field>>& reg,
const DistMap& invMap,
const ldl::DistNodeInfo& info,
const ldl::DistFront<Field>& front,
const DistSparseLDLFactorization<Field>& sparseLDLFact,
DistMultiVec<Field>& B,
Base<Field> relTolRefine,
Int maxRefineIts,
bool progress=false,
bool time=false );
template<typename Field>
Int RegularizedSolveAfter
( const DistSparseMatrix<Field>& A,
const DistMultiVec<Base<Field>>& reg,
const DistMap& invMap,
const ldl::DistNodeInfo& info,
const ldl::DistFront<Field>& front,
DistMultiVec<Field>& B,
ldl::DistMultiVecNodeMeta& meta,
Base<Field> relTolRefine,
Int maxRefineIts,
bool progress=false,
bool time=false );

template<typename Field>
Int RegularizedSolveAfter
( const SparseMatrix<Field>& A,
const Matrix<Base<Field>>& reg,
const Matrix<Base<Field>>& d,
const vector<Int>& invMap,
const ldl::NodeInfo& info,
const ldl::Front<Field>& front,
const SparseLDLFactorization<Field>& sparseLDLFact,
Matrix<Field>& B,
Base<Field> relTolRefine,
Int maxRefineIts,
Expand All @@ -370,88 +338,43 @@ Int RegularizedSolveAfter
( const DistSparseMatrix<Field>& A,
const DistMultiVec<Base<Field>>& reg,
const DistMultiVec<Base<Field>>& d,
const DistMap& invMap,
const ldl::DistNodeInfo& info,
const ldl::DistFront<Field>& front,
const DistSparseLDLFactorization<Field>& sparseLDLFact,
DistMultiVec<Field>& B,
Base<Field> relTolRefine,
Int maxRefineIts,
bool progress=false,
bool time=false );
template<typename Field>
Int RegularizedSolveAfter
( const DistSparseMatrix<Field>& A,
const DistMultiVec<Base<Field>>& reg,
const DistMultiVec<Base<Field>>& d,
const DistMap& invMap,
const ldl::DistNodeInfo& info,
const ldl::DistFront<Field>& front,
DistMultiVec<Field>& B,
ldl::DistMultiVecNodeMeta& meta,
Base<Field> relTolRefine,
Int maxRefineIts,
bool progress=false,
bool time=false );

template<typename Field>
Int SolveAfter
( const SparseMatrix<Field>& A,
const Matrix<Base<Field>>& reg,
const vector<Int>& invMap,
const ldl::NodeInfo& info,
const ldl::Front<Field>& front,
const SparseLDLFactorization<Field>& sparseLDLFact,
Matrix<Field>& B,
const RegSolveCtrl<Base<Field>>& ctrl );
template<typename Field>
Int SolveAfter
( const DistSparseMatrix<Field>& A,
const DistMultiVec<Base<Field>>& reg,
const DistMap& invMap,
const ldl::DistNodeInfo& info,
const ldl::DistFront<Field>& front,
DistMultiVec<Field>& B,
const RegSolveCtrl<Base<Field>>& ctrl );
template<typename Field>
Int SolveAfter
( const DistSparseMatrix<Field>& A,
const DistMultiVec<Base<Field>>& reg,
const DistMap& invMap,
const ldl::DistNodeInfo& info,
const ldl::DistFront<Field>& front,
const DistSparseLDLFactorization<Field>& sparseLDLFact,
DistMultiVec<Field>& B,
ldl::DistMultiVecNodeMeta& meta,
const RegSolveCtrl<Base<Field>>& ctrl );

template<typename Field>
Int SolveAfter
( const SparseMatrix<Field>& A,
const Matrix<Base<Field>>& reg,
const Matrix<Base<Field>>& d,
const vector<Int>& invMap,
const ldl::NodeInfo& info,
const ldl::Front<Field>& front,
const SparseLDLFactorization<Field>& sparseLDLFact,
Matrix<Field>& B,
const RegSolveCtrl<Base<Field>>& ctrl );
template<typename Field>
Int SolveAfter
( const DistSparseMatrix<Field>& A,
const DistMultiVec<Base<Field>>& reg,
const DistMultiVec<Base<Field>>& d,
const DistMap& invMap,
const ldl::DistNodeInfo& info,
const ldl::DistFront<Field>& front,
DistMultiVec<Field>& B,
const RegSolveCtrl<Base<Field>>& ctrl );
template<typename Field>
Int SolveAfter
( const DistSparseMatrix<Field>& A,
const DistMultiVec<Base<Field>>& reg,
const DistMultiVec<Base<Field>>& d,
const DistMap& invMap,
const ldl::DistNodeInfo& info,
const ldl::DistFront<Field>& front,
const DistSparseLDLFactorization<Field>& sparseLDLFact,
DistMultiVec<Field>& B,
ldl::DistMultiVecNodeMeta& meta,
const RegSolveCtrl<Base<Field>>& ctrl );

} // namespace reg_ldl
Expand Down
14 changes: 14 additions & 0 deletions include/El/lapack_like/factor/ldl/sparse/numeric.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -529,6 +529,13 @@ class SparseLDLFactorization
// Overwrite 'B' with the solution to 'A X = B'.
void Solve( Matrix<Field>& B ) const;

// Overwrite 'B' with the solution to 'A X = B' using Iterative Refinement.
void SolveWithIterativeRefinement
( const SparseMatrix<Field>& A,
Matrix<Field>& B,
const Base<Field>& relTolRefine,
Int maxRefineIts ) const;

ldl::Front<Field>& Front();
const ldl::Front<Field>& Front() const;

Expand Down Expand Up @@ -585,6 +592,13 @@ class DistSparseLDLFactorization
// Overwrite 'B' with the solution to 'A X = B'.
void Solve( DistMultiVec<Field>& B ) const;

// Overwrite 'B' with the solution to 'A X = B' using Iterative Refinement.
void SolveWithIterativeRefinement
( const DistSparseMatrix<Field>& A,
DistMultiVec<Field>& B,
const Base<Field>& relTolRefine,
Int maxRefineIts ) const;

ldl::DistFront<Field>& Front();
const ldl::DistFront<Field>& Front() const;

Expand Down
61 changes: 61 additions & 0 deletions src/lapack_like/factor/LDL/sparse/numeric/DistSparseLDLFactorization.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -108,6 +108,67 @@ void DistSparseLDLFactorization<Field>::Solve( DistMultiVec<Field>& B ) const
ldl::SolveAfter( inverseMap_, *info_, *front_, B );
}

template<typename Field>
void DistSparseLDLFactorization<Field>::SolveWithIterativeRefinement
( const DistSparseMatrix<Field>& A,
DistMultiVec<Field>& B,
const Base<Field>& minReductionFactor,
Int maxRefineIts ) const
{
EL_DEBUG_CSE
const Grid& grid = B.Grid();
// TODO(poulson): Generalize this implementation
if( B.Width() > 1 )
LogicError("Iterative Refinement currently only supported for one RHS");

DistMultiVec<Field> BOrig(grid);
BOrig = B;

// Compute the initial guess
// =========================
DistMultiVec<Field> X(B);
Solve( X );

Int refineIt = 0;
if( maxRefineIts > 0 )
{
DistMultiVec<Field> dX(grid), XCand(grid);
Multiply( NORMAL, Field(-1), A, X, Field(1), B );
Base<Field> errorNorm = FrobeniusNorm( B );
for( ; refineIt<maxRefineIts; ++refineIt )
{
// Compute the proposed update to the solution
// -------------------------------------------
dX = B;
Solve( dX );
XCand = X;
XCand += dX;

// If the proposed update lowers the residual, accept it
// -----------------------------------------------------
B = BOrig;
Multiply( NORMAL, Field(-1), A, XCand, Field(1), B );
Base<Field> newErrorNorm = FrobeniusNorm( B );
if( minReductionFactor*newErrorNorm < errorNorm )
{
X = XCand;
errorNorm = newErrorNorm;
}
else if( newErrorNorm < errorNorm )
{
X = XCand;
errorNorm = newErrorNorm;
break;
}
else
break;
}
}
// Store the final result
// ======================
B = X;
}

template<typename Field>
ldl::DistFront<Field>& DistSparseLDLFactorization<Field>::Front()
{
Expand Down
Loading

0 comments on commit 3ea8fc2

Please sign in to comment.