Skip to content

Commit

Permalink
Modifying Bunch-Kaufman to return the subdiagonal entries of D so tha…
Browse files Browse the repository at this point in the history 
…t only L's main diagonal need be implicitly modified, and returning pivot vector in standard form which does not use negative entries to signal 2x2 pivots. QuasiDiagonalSolve was then added, as well as ldl::SolveAfter for BunchKaufman.
  • Loading branch information
@poulson
poulson committed Sep 30, 2013
1 parent f288c2b commit d2c79c6
Show file tree
Hide file tree
Showing 17 changed files with 1,102 additions and 684 deletions.
2 changes: 2 additions & 0 deletions TODO
View file
@@ -1,5 +1,6 @@
Functionality
=============
- QuasiDiagonalScale (and finish all options for QuasiDiagonalSolve)
- Distribute between different grids for any distribution
- Consistently implement 'GetDiagonal' functionality for DistMatrix
- Ability to change default grid with command-line options
Expand Down Expand Up @@ -38,6 +39,7 @@ Functionality

Maintenance
===========
- Documentation for QuasiDiagonalSolve and Symmetric2x2Solve
- Make transpose-options of LocalTrr(2)k more consistent with Trr(2)k
- README's in each directory explaining the contents of the folder
- Support for OpenBLAS [-D MATH_LIBS="-lopenblas;-lpthread;-lm;-lgfortran"]
Expand Down
38 changes: 16 additions & 22 deletions doc/source/lapack-like/factorizations.rst
View file
Expand Up @@ -41,11 +41,13 @@ be thrown.
Overwrite the `uplo` triangle of the potentially singular matrix `A` with
its Cholesky factor.

:math:`LDL^H` factorization
---------------------------
:math:`LDL` factorization
-------------------------

.. cpp:function:: void LDLH( Matrix<F>& A, Matrix<Int>& p )
.. cpp:function:: void LDLH( DistMatrix<F>& A, DistMatrix<Int,VC,STAR>& p )
.. cpp:function:: void LDLH( Matrix<F>& A, Matrix<F>& dSub, Matrix<int>& p )
.. cpp:function:: void LDLT( Matrix<F>& A, Matrix<F>& dSub, Matrix<int>& p )
.. cpp:function:: void LDLH( DistMatrix<F>& A, DistMatrix<F,MD,STAR>& dSub, DistMatrix<int,VC,STAR>& p )
.. cpp:function:: void LDLT( DistMatrix<F>& A, DistMatrix<F,MD,STAR>& dSub, DistMatrix<int,VC,STAR>& p )

Blocked Bunch-Kaufman. More details to come.

Expand All @@ -60,34 +62,26 @@ attempted, then a :cpp:type:`SingularMatrixException` will be thrown.
The following routines do not pivot, so please use with caution.

.. cpp:function:: void LDLH( Matrix<F>& A )
.. cpp:function:: void LDLT( Matrix<F>& A )
.. cpp:function:: void LDLH( DistMatrix<F>& A )
.. cpp:function:: void LDLT( DistMatrix<F>& A )

Overwrite the strictly lower triangle of :math:`A` with the strictly lower
portion of :math:`L` (:math:`L` implicitly has ones on its diagonal) and
the diagonal with :math:`D`.

:math:`LDL^T` factorization
---------------------------

.. cpp:function:: void LDLT( Matrix<F>& A, Matrix<Int>& p )
.. cpp:function:: void LDLT( DistMatrix<F>& A, DistMatrix<Int,VC,STAR>& p )

Blocked Bunch-Kaufman. More details to come.

While the :math:`LDL^H` factorization targets Hermitian matrices, the
:math:`LDL^T` factorization targets symmetric matrices. If a zero pivot is
attempted, then a :cpp:type:`SingularMatrixException` will be thrown.
Detailed interface
^^^^^^^^^^^^^^^^^^

.. warning::
.. cpp:function:: ldl::SolveAfter( const Matrix<F>& A, Matrix<F>& B, bool conjugated=false )
.. cpp:function:: ldl::SolveAfter( const DistMatrix<F>& A, DistMatrix<F>& B, bool conjugated=false )

The following routines do not pivot, so please use with caution.
Solve linear systems using an unpivoted LDL factorization.

.. cpp:function:: void LDLT( Matrix<F>& A )
.. cpp:function:: void LDLT( DistMatrix<F>& A )
.. cpp:function:: ldl::SolveAfter( const Matrix<F>& A, const Matrix<F>& dSub, const Matrix<int>& p, Matrix<F>& B, bool conjugated=false )
.. cpp:function:: ldl::SolveAfter( const DistMatrix<F>& A, const DistMatrix<F,MD,STAR>& dSub, const DistMatrix<int,VC,STAR>& p, DistMatrix<F>& B, bool conjugated=false )

Overwrite the strictly lower triangle of :math:`A` with the strictly lower
portion of :math:`L` (:math:`L` implicitly has ones on its diagonal) and
the diagonal with :math:`D`.
Solve linear systems using a pivoted LDL factorization.

:math:`LU` factorization
------------------------
Expand Down
48 changes: 43 additions & 5 deletions examples/lapack-like/BunchKaufman.cpp
View file
Expand Up @@ -14,10 +14,13 @@
#include "elemental/blas-like/level1/MakeTriangular.hpp"
#include "elemental/blas-like/level1/SetDiagonal.hpp"
#include "elemental/blas-like/level1/Transpose.hpp"
#include "elemental/blas-like/level3/Symm.hpp"
#include "elemental/lapack-like/ApplyPackedReflectors/Util.hpp"
#include "elemental/lapack-like/LDL.hpp"
#include "elemental/lapack-like/Norm/Frobenius.hpp"
#include "elemental/matrices/HermitianUniformSpectrum.hpp"
#include "elemental/matrices/Uniform.hpp"
#include "elemental/matrices/Zeros.hpp"
using namespace std;
using namespace elem;

Expand All @@ -34,10 +37,12 @@ main( int argc, char* argv[] )
{
const Int n = Input("--size","size of matrix to factor",100);
const Int nb = Input("--nb","algorithmic blocksize",96);
const Int numRhs = Input("--numRhs","number of random r.h.s.",100);
const double realMean = Input("--realMean","real mean",0.);
const double imagMean = Input("--imagMean","imag mean",0.);
const double stddev = Input("--stddev","standard dev.",1.);
const bool conjugate = Input("--conjugate","LDL^H?",false);
const bool print = Input("--print","print matrices?",false);
ProcessInput();
PrintInputReport();

Expand All @@ -59,15 +64,48 @@ main( int argc, char* argv[] )

// Make a copy of A and then overwrite it with its LDL factorization
DistMatrix<Int,VC,STAR> p;
DistMatrix<C,MD,STAR> dSub;
DistMatrix<C> factA( A );
MakeTriangular( LOWER, factA );
if( conjugate )
LDLH( factA, p );
LDLH( factA, dSub, p );
else
LDLT( factA, p );
Print( A, "A" );
Print( factA, "factA" );
Print( p, "p" );
LDLT( factA, dSub, p );
if( print )
{
Print( A, "A" );
Print( factA, "factA" );
Print( dSub, "dSub" );
Print( p, "p" );
}

// Generate a random set of vectors
DistMatrix<C> X;
Uniform( X, n, numRhs );
DistMatrix<C> B;
Zeros( B, n, numRhs );
Symm( LEFT, LOWER, C(1), A, X, C(0), B, conjugate );
if( print )
{
Print( X, "X" );
Print( B, "B" );
}
ldl::SolveAfter( factA, dSub, p, B, conjugate );
const Real AFrob = HermitianFrobeniusNorm( LOWER, A );
const Real XFrob = FrobeniusNorm( X );
Axpy( C(-1), B, X );
const Real errFrob = FrobeniusNorm( X );
if( print )
{
Print( B, "XComp" );
Print( X, "E" );
}
if( mpi::WorldRank() == 0 )
{
std::cout << "|| A ||_F = " << AFrob << "\n"
<< "|| X ||_F = " << XFrob << "\n"
<< "|| X - inv(A) A X ||_F = " << errFrob << std::endl;
}
}
catch( exception& e ) { ReportException(e); }

Expand Down
48 changes: 34 additions & 14 deletions examples/lapack-like/SequentialBunchKaufman.cpp
View file
Expand Up @@ -14,10 +14,13 @@
#include "elemental/blas-like/level1/MakeTriangular.hpp"
#include "elemental/blas-like/level1/SetDiagonal.hpp"
#include "elemental/blas-like/level1/Transpose.hpp"
#include "elemental/blas-like/level3/Symm.hpp"
#include "elemental/lapack-like/ApplyPackedReflectors/Util.hpp"
#include "elemental/lapack-like/LDL.hpp"
#include "elemental/lapack-like/Norm/Frobenius.hpp"
#include "elemental/matrices/HermitianUniformSpectrum.hpp"
#include "elemental/matrices/Uniform.hpp"
#include "elemental/matrices/Zeros.hpp"
using namespace std;
using namespace elem;

Expand All @@ -34,10 +37,12 @@ main( int argc, char* argv[] )
{
const Int n = Input("--size","size of matrix to factor",100);
const Int nb = Input("--nb","algorithmic blocksize",96);
const Int numRhs = Input("--numRhs","number of random r.h.s.",100);
const double realMean = Input("--realMean","real mean",0.);
const double imagMean = Input("--imagMean","imag mean",0.);
const double stddev = Input("--stddev","standard dev.",1.);
const bool conjugate = Input("--conjugate","LDL^H?",false);
const bool print = Input("--print","print matrices?",false);
ProcessInput();
PrintInputReport();

Expand All @@ -59,31 +64,46 @@ main( int argc, char* argv[] )

// Make a copy of A and then overwrite it with its LDL factorization
Matrix<Int> p;
Matrix<C> factA( A );
Matrix<C> dSub, factA( A );
MakeTriangular( LOWER, factA );
if( conjugate )
LDLH( factA, p );
LDLH( factA, dSub, p );
else
LDLT( factA, p );
if( mpi::WorldRank() == 0 )
LDLT( factA, dSub, p );
if( print && mpi::WorldRank()==0 )
{
Print( A, "A" );
Print( factA, "factA" );
Print( dSub, "dSub" );
Print( p, "p" );
}

// Do the same thing with the unblocked algorithm
Matrix<Int> pUnb;
Matrix<C> factAUnb( A );
MakeTriangular( LOWER, factAUnb );
if( conjugate )
ldl::UnblockedPivoted( ADJOINT, factAUnb, pUnb );
else
ldl::UnblockedPivoted( TRANSPOSE, factAUnb, pUnb );
// Generate a random set of vectors
Matrix<C> X;
Uniform( X, n, numRhs );
Matrix<C> B;
Zeros( B, n, numRhs );
Symm( LEFT, LOWER, C(1), A, X, C(0), B, conjugate );
if( print && mpi::WorldRank()==0 )
{
Print( X, "X" );
Print( B, "B" );
}
ldl::SolveAfter( factA, dSub, p, B, conjugate );
const Real AFrob = HermitianFrobeniusNorm( LOWER, A );
const Real XFrob = FrobeniusNorm( X );
Axpy( C(-1), B, X );
const Real errFrob = FrobeniusNorm( X );
if( print && mpi::WorldRank() == 0 )
{
Print( B, "XComp" );
Print( X, "E" );
}
if( mpi::WorldRank() == 0 )
{
Print( factAUnb, "factAUnb" );
Print( pUnb, "pUnb" );
std::cout << "|| A ||_F = " << AFrob << "\n"
<< "|| X ||_F = " << XFrob << "\n"
<< "|| X - inv(A) A X ||_F = " << errFrob << std::endl;
}
}
catch( exception& e ) { ReportException(e); }
Expand Down
1 change: 1 addition & 0 deletions include/elemental/blas-like/level1.hpp
View file
Expand Up @@ -26,6 +26,7 @@
#include "./level1/MakeTriangular.hpp"
#include "./level1/Max.hpp"
#include "./level1/Nrm2.hpp"
#include "./level1/QuasiDiagonalSolve.hpp"
#include "./level1/Scale.hpp"
#include "./level1/ScaleTrapezoid.hpp"
#include "./level1/SetDiagonal.hpp"
Expand Down
93 changes: 8 additions & 85 deletions include/elemental/blas-like/level1/DiagonalScale.hpp
View file
Expand Up @@ -12,11 +12,11 @@

namespace elem {

template<typename T>
template<typename T,typename TDiag>
inline void
DiagonalScale
( LeftOrRight side, Orientation orientation,
const Matrix<T>& d, Matrix<T>& X )
const Matrix<TDiag>& d, Matrix<T>& X )
{
#ifndef RELEASE
CallStackEntry entry("DiagonalScale");
Expand Down Expand Up @@ -54,48 +54,13 @@ DiagonalScale
}
}

template<typename T>
template<typename T,typename TDiag,
Distribution U,Distribution V,
Distribution W,Distribution Z>
inline void
DiagonalScale
( LeftOrRight side, Orientation orientation,
const Matrix<BASE(T)>& d, Matrix<T>& X )
{
#ifndef RELEASE
CallStackEntry entry("DiagonalScale");
#endif
typedef BASE(T) R;

const Int m = X.Height();
const Int n = X.Width();
const Int ldim = X.LDim();
if( side == LEFT )
{
for( Int i=0; i<m; ++i )
{
const R delta = d.Get(i,0);
T* XBuffer = X.Buffer(i,0);
for( Int j=0; j<n; ++j )
XBuffer[j*ldim] *= delta;
}
}
else
{
for( Int j=0; j<n; ++j )
{
const R delta = d.Get(j,0);
T* XBuffer = X.Buffer(0,j);
for( Int i=0; i<m; ++i )
XBuffer[i] *= delta;
}
}
}

template<typename T,Distribution U,Distribution V,
Distribution W,Distribution Z>
inline void
DiagonalScale
( LeftOrRight side, Orientation orientation,
const DistMatrix<T,U,V>& d, DistMatrix<T,W,Z>& X )
const DistMatrix<TDiag,U,V>& d, DistMatrix<T,W,Z>& X )
{
#ifndef RELEASE
CallStackEntry entry("DiagonalScale");
Expand All @@ -108,7 +73,7 @@ DiagonalScale
}
else
{
DistMatrix<T,W,STAR> d_W_STAR( X.Grid() );
DistMatrix<TDiag,W,STAR> d_W_STAR( X.Grid() );
d_W_STAR = d;
DiagonalScale
( LEFT, orientation,
Expand All @@ -124,49 +89,7 @@ DiagonalScale
}
else
{
DistMatrix<T,Z,STAR> d_Z_STAR( X.Grid() );
d_Z_STAR = d;
DiagonalScale
( RIGHT, orientation, d_Z_STAR.LockedMatrix(), X.Matrix() );
}
}
}

template<typename T,Distribution U,Distribution V,
Distribution W,Distribution Z>
inline void
DiagonalScale
( LeftOrRight side, Orientation orientation,
const DistMatrix<BASE(T),U,V>& d, DistMatrix<T,W,Z>& X )
{
#ifndef RELEASE
CallStackEntry entry("DiagonalScale");
#endif
typedef BASE(T) R;

if( side == LEFT )
{
if( U == W && V == STAR && d.ColAlign() == X.ColAlign() )
{
DiagonalScale( LEFT, orientation, d.LockedMatrix(), X.Matrix() );
}
else
{
DistMatrix<R,W,STAR> d_W_STAR( X.Grid() );
d_W_STAR = d;
DiagonalScale
( LEFT, orientation, d_W_STAR.LockedMatrix(), X.Matrix() );
}
}
else
{
if( U == Z && V == STAR && d.ColAlign() == X.RowAlign() )
{
DiagonalScale( RIGHT, orientation, d.LockedMatrix(), X.Matrix() );
}
else
{
DistMatrix<R,Z,STAR> d_Z_STAR( X.Grid() );
DistMatrix<TDiag,Z,STAR> d_Z_STAR( X.Grid() );
d_Z_STAR = d;
DiagonalScale
( RIGHT, orientation, d_Z_STAR.LockedMatrix(), X.Matrix() );
Expand Down

0 comments on commit d2c79c6

Please sign in to comment.