template<typename T>
void Herk
( UpperOrLower uplo, Orientation orientation,
  Base<T> alpha, const Matrix<T>& A, Base<T> beta, Matrix<T>& C )
{
    DEBUG_ONLY(CSE cse("Herk"))
    Syrk( uplo, orientation, T(alpha), A, T(beta), C, true );
}

template<typename T>
void Herk
( UpperOrLower uplo, Orientation orientation,
  Base<T> alpha, const Matrix<T>& A, Matrix<T>& C )
{
    DEBUG_ONLY(CSE cse("Herk"))
    const Int n = ( orientation==NORMAL ? A.Height() : A.Width() );
    C.Resize( n, n );
    Zero( C );
    Syrk( uplo, orientation, T(alpha), A, T(0), C, true );
}

template<typename T>
void Herk
( UpperOrLower uplo, Orientation orientation,
  Base<T> alpha, const ElementalMatrix<T>& A, 
  Base<T> beta,        ElementalMatrix<T>& C )
{
    DEBUG_ONLY(CSE cse("Herk"))
    Syrk( uplo, orientation, T(alpha), A, T(beta), C, true );
}

/*
   Copyright (c) 2009-2015, Jack Poulson
   All rights reserved.

   This file is part of Elemental and is under the BSD 2-Clause License, 
   which can be found in the LICENSE file in the root directory, or at 
   http://opensource.org/licenses/BSD-2-Clause
*/
#include "El.hpp"

namespace El {

template<typename F> 
Base<F> OneCondition( const Matrix<F>& A )
{
    DEBUG_ONLY(CallStackEntry cse("OneCondition"))
    typedef Base<F> Real;
    Matrix<F> B( A );
    const Real oneNorm = OneNorm( B );
    try { Inverse( B ); }
    catch( SingularMatrixException& e ) 
    { return std::numeric_limits<Real>::infinity(); }
    const Real oneNormInv = OneNorm( B );
    return oneNorm*oneNormInv;
}

template<typename F> 
Base<F> OneCondition( const AbstractDistMatrix<F>& A )
{
    DEBUG_ONLY(CallStackEntry cse("OneCondition"))
    typedef Base<F> Real;

//   | -I 0   (-z <> x) | | dz | = | z <> rmu |
//
// where 
//
//   rc = A^T y - z + c,
//   rb = A x - b,
//   rmu = x o z - tau e

template<typename Real>
void KKT
( const Matrix<Real>& A, 
  const Matrix<Real>& x,
  const Matrix<Real>& z,
        Matrix<Real>& J, bool onlyLower )
{
    DEBUG_ONLY(CSE cse("lp::direct::KKT"))
    const Int m = A.Height();
    const Int n = A.Width();

    Zeros( J, 2*n+m, 2*n+m );
    const IR xInd(0,n), yInd(n,n+m), zInd(n+m,2*n+m);
    auto Jxx = J(xInd,xInd); auto Jxy = J(xInd,yInd); auto Jxz = J(xInd,zInd); 
    auto Jyx = J(yInd,xInd); auto Jyy = J(yInd,yInd); auto Jyz = J(yInd,zInd); 
    auto Jzx = J(zInd,xInd); auto Jzy = J(zInd,yInd); auto Jzz = J(zInd,zInd); 

    // Jyx := A
    // ========
    Jyx = A;

    // Jzx := -I
    // =========

//   s.t. |  A  -1 | | x | <= |  b |
//        | -A  -1 | | t |    | -b |
//
// NOTE: There is likely an appropriate citation, but the derivation is 
//       trivial. If one is found, it will be added.

namespace El {

template<typename Real>
void CP
( const Matrix<Real>& A,
  const Matrix<Real>& b, 
        Matrix<Real>& x,
  const lp::affine::Ctrl<Real>& ctrl )
{
    DEBUG_ONLY(CSE cse("CP"))
    const Int m = A.Height();
    const Int n = A.Width();
    Matrix<Real> c, AHat, bHat, G, h;

    // c := [zeros(n,1);1]
    // ===================
    Zeros( c, n+1, 1 );
    c.Set( n, 0, Real(1) );

    // \hat A := zeros(0,n+1)
    // ====================== 
    Zeros( AHat, 0, n+1 );

    // \hat b := zeros(0,1)
    // ====================

#pragma once
#ifndef ELEM_BIDIAG_APPLY_HPP
#define ELEM_BIDIAG_APPLY_HPP

#include ELEM_APPLYPACKEDREFLECTORS_INC

namespace elem {
namespace bidiag {

template<typename F>
inline void
ApplyQ
( LeftOrRight side, Orientation orientation, 
  const Matrix<F>& A, const Matrix<F>& t, Matrix<F>& B )
{
    DEBUG_ONLY(CallStackEntry cse("bidiag::ApplyQ"))
    const bool normal = (orientation==NORMAL);
    const bool onLeft = (side==LEFT);
    const ForwardOrBackward direction = ( normal==onLeft ? BACKWARD : FORWARD );
    const Conjugation conjugation = ( normal ? CONJUGATED : UNCONJUGATED );
    const Int offset = ( A.Height()>=A.Width() ? 0 : -1 );
    ApplyPackedReflectors
    ( side, LOWER, VERTICAL, direction, conjugation, offset, A, t, B );
}

template<typename F>
inline void
ApplyP
( LeftOrRight side, Orientation orientation, 
  const Matrix<F>& A, const Matrix<F>& t, Matrix<F>& B )
{

   All rights reserved.
   
   This file is part of Elemental and is under the BSD 2-Clause License, 
   which can be found in the LICENSE file in the root directory, or at 
   http://opensource.org/licenses/BSD-2-Clause
*/
#include "El.hpp"

namespace El {
namespace ldl {

template<typename F>
void DiagonalScale
( const NodeInfo& info, const Front<F>& front, MatrixNode<F>& X )
{
    DEBUG_ONLY(CSE cse("ldl::DiagonalScale"))

    const Int numChildren = info.children.size();
    for( Int c=0; c<numChildren; ++c )
        DiagonalScale( *info.children[c], *front.children[c], *X.children[c] );

    if( PivotedFactorization(front.type) )
        QuasiDiagonalScale
        ( LEFT, LOWER, front.diag, front.subdiag, 
          X.matrix, front.isHermitian );
    else
        DiagonalScale( LEFT, NORMAL, front.diag, X.matrix );
}

template<typename F>
void DiagonalScale

   This file is part of Elemental and is under the BSD 2-Clause License, 
   which can be found in the LICENSE file in the root directory, or at 
   http://opensource.org/licenses/BSD-2-Clause
*/
#include <El-lite.hpp>
#include <El/blas_like/level1.hpp>
#include <El/matrices.hpp>

namespace El {

// Draw each entry from a normal PDF
template<typename F>
void MakeGaussian( Matrix<F>& A, F mean, Base<F> stddev )
{
    DEBUG_ONLY(CSE cse("MakeGaussian"))
    auto sampleNormal = [=]() { return SampleNormal(mean,stddev); };
    EntrywiseFill( A, function<F()>(sampleNormal) );
}

template<typename F>
void MakeGaussian( AbstractDistMatrix<F>& A, F mean, Base<F> stddev )
{
    DEBUG_ONLY(CSE cse("MakeGaussian"))
    if( A.RedundantRank() == 0 )
        MakeGaussian( A.Matrix(), mean, stddev );
    Broadcast( A, A.RedundantComm(), 0 );
}

template<typename F>
void MakeGaussian( DistMultiVec<F>& A, F mean, Base<F> stddev )

inline
DistNodalMultiVec<F>::DistNodalMultiVec( const DistNodalMatrix<F>& X )
{
    DEBUG_ONLY(CallStackEntry cse("DistNodalMultiVec::DistNodalMultiVec"))
    *this = X;
}

   This file is part of Elemental and is under the BSD 2-Clause License, 
   which can be found in the LICENSE file in the root directory, or at 
   http://opensource.org/licenses/BSD-2-Clause
*/
#pragma once
#ifndef EL_TRSTRM_LLN_HPP
#define EL_TRSTRM_LLN_HPP

namespace El {
namespace trstrm {

template<typename F>
inline void
LLNUnb( UnitOrNonUnit diag, F alpha, const Matrix<F>& L, Matrix<F>& X )
{
    DEBUG_ONLY(CSE cse("trstrm::LLNUnb"))
    const bool isUnit = ( diag==UNIT );
    const Int n = L.Height();
    const Int LLDim = L.LDim();
    const Int XLDim = X.LDim();
    const F* LBuffer = L.LockedBuffer();
    F* XBuffer = X.Buffer();

    // X := alpha X
    if( alpha != F(1) )
        for( Int j=0; j<n; ++j ) 
            for( Int i=j; i<n; ++i )
                XBuffer[i+j*XLDim] *= alpha;

    for( Int i=0; i<n; ++i )
    {

BDM& BDM::operator=( const DistMatrix<T,STAR,MR,BLOCK_CYCLIC>& A )
{ 
    DEBUG_ONLY(CSE cse("[STAR,VR] = [STAR,MR]"))
    copy::PartialRowFilter( A, *this );
    return *this;
}

#pragma once
#ifndef EL_BIDIAG_L_HPP
#define EL_BIDIAG_L_HPP

#include "./LUnb.hpp"
#include "./LPan.hpp"

namespace El {
namespace bidiag {

// NOTE: Very little is changed versus the upper case. Perhaps they should be
//       combined.
template<typename F>
inline void L( Matrix<F>& A, Matrix<F>& tP, Matrix<F>& tQ )
{
    DEBUG_ONLY(CSE cse("bidiag::L"))
    const Int m = A.Height();
    const Int n = A.Width();
    DEBUG_ONLY(
      if( m > n )
          LogicError("A must be at least as wide as it is tall");
      // Are these requirements necessary?!?
      if( tP.Viewing() || tQ.Viewing() )
          LogicError("tP and tQ must not be views");
    )
    const Int tPHeight = m;
    const Int tQHeight = Max(m-1,0);
    tP.Resize( tPHeight, 1 );
    tQ.Resize( tQHeight, 1 );

    Matrix<F> X, Y;

void ColAllToAllDemote
( const DistMatrix<T,Partial<U>(),PartialUnionRow<U,V>()>& A,
        DistMatrix<T,        U,                     V   >& B )
{
    DEBUG_ONLY(CSE cse("copy::ColAllToAllDemote"))
    AssertSameGrids( A, B );

    const Int height = A.Height();
    const Int width = A.Width();
    B.AlignColsAndResize( A.ColAlign(), height, width, false, false );
    if( !B.Participating() )
        return;

    const Int colAlign = B.ColAlign();
    const Int rowAlignA = A.RowAlign();

    const Int colStride = B.ColStride();
    const Int colStridePart = B.PartialColStride();
    const Int colStrideUnion = B.PartialUnionColStride();
    const Int colRankPart = B.PartialColRank();
    const Int colDiff = (colAlign%colStridePart) - A.ColAlign();

    const Int colShiftA = A.ColShift();

    const Int localHeightB = B.LocalHeight();
    const Int localWidthA = A.LocalWidth();
    const Int maxLocalHeight = MaxLength(height,colStride);
    const Int maxLocalWidth = MaxLength(width,colStrideUnion);
    const Int portionSize = mpi::Pad( maxLocalHeight*maxLocalWidth );

    vector<T> buffer( 2*colStrideUnion*portionSize );
    T* firstBuf  = &buffer[0];
    T* secondBuf = &buffer[colStrideUnion*portionSize];

    if( colDiff == 0 )
    {
        // Pack            
        util::PartialColStridedPack
        ( height, localWidthA,
          colAlign, colStride, 
          colStrideUnion, colStridePart, colRankPart,
          colShiftA,
          A.LockedBuffer(), A.LDim(),
          firstBuf,         portionSize );

        // Simultaneously Scatter in columns and Gather in rows
        mpi::AllToAll
        ( firstBuf,  portionSize,
          secondBuf, portionSize, B.PartialUnionColComm() );

        // Unpack
        util::RowStridedUnpack
        ( localHeightB, width,
          rowAlignA, colStrideUnion,
          secondBuf, portionSize,
          B.Buffer(), B.LDim() );
    }
    else
    {
#ifdef EL_UNALIGNED_WARNINGS
        if( B.Grid().Rank() == 0 )
            cerr << "Unaligned ColAllToAllDemote" << endl;
#endif
        const Int sendColRankPart = Mod( colRankPart+colDiff, colStridePart );
        const Int recvColRankPart = Mod( colRankPart-colDiff, colStridePart );

        // Pack
        util::PartialColStridedPack
        ( height, localWidthA,
          colAlign, colStride, 
          colStrideUnion, colStridePart, sendColRankPart,
          colShiftA,
          A.LockedBuffer(), A.LDim(),
          secondBuf,        portionSize );

        // Simultaneously Scatter in columns and Gather in rows
        mpi::AllToAll
        ( secondBuf, portionSize,
          firstBuf,  portionSize, B.PartialUnionColComm() );

        // Realign the result
        mpi::SendRecv
        ( firstBuf,  colStrideUnion*portionSize, sendColRankPart,
          secondBuf, colStrideUnion*portionSize, recvColRankPart,
          B.PartialColComm() );

        // Unpack
        util::RowStridedUnpack
        ( localHeightB, width,
          rowAlignA, colStrideUnion,
          secondBuf, portionSize,
          B.Buffer(), B.LDim() );
    }
}

#include "./Write/Ascii.hpp"
#include "./Write/AsciiMatlab.hpp"
#include "./Write/Binary.hpp"
#include "./Write/BinaryFlat.hpp"
#include "./Write/Image.hpp"
#include "./Write/MatrixMarket.hpp"

namespace El {

template<typename T>
void Write
( const Matrix<T>& A, 
  std::string basename, FileFormat format, std::string title )
{
    DEBUG_ONLY(CallStackEntry cse("Write"))
    switch( format )
    {
    case ASCII:         write::Ascii( A, basename, title );       break;
    case ASCII_MATLAB:  write::AsciiMatlab( A, basename, title ); break;
    case BINARY:        write::Binary( A, basename );             break;
    case BINARY_FLAT:   write::BinaryFlat( A, basename );         break;
    case MATRIX_MARKET: write::MatrixMarket( A, basename );       break;
    case BMP:
    case JPG:
    case JPEG:
    case PNG:
    case PPM:
    case XBM:
    case XPM:
        write::Image( A, basename, format ); break;

/*
   Copyright (c) 2009-2015, Jack Poulson
   All rights reserved.

   This file is part of Elemental and is under the BSD 2-Clause License, 
   which can be found in the LICENSE file in the root directory, or at 
   http://opensource.org/licenses/BSD-2-Clause
*/
#include "El.hpp"

namespace El {

template<typename T> 
void Wilkinson( Matrix<T>& A, Int k )
{
    DEBUG_ONLY(CSE cse("Wilkinson"))
    const Int n = 2*k+1;
    Zeros( A, n, n );
    FillDiagonal( A, T(1), -1 );
    FillDiagonal( A, T(1),  1 );
    
    for( Int j=0; j<=k; ++j )
        A.Set( j, j, T(k-j) );
    for( Int j=k+1; j<n; ++j )
        A.Set( j, j, T(j-k) );
}

template<typename T>
void Wilkinson( AbstractDistMatrix<T>& A, Int k )
{
    DEBUG_ONLY(CSE cse("Wilkinson"))

    // Set the singular values to the absolute value of the eigenvalues
    auto absLambda = []( Real sigma ) { return Abs(sigma); };
    EntrywiseMap( s, function<Real(Real)>(absLambda) );

    // TODO: Descending sort of triplets
}

// Return the singular values
// ==========================
// NOTE: A is ovewritten with its packed reduction to tridiagonal form

template<typename F>
void HermitianSVD( UpperOrLower uplo, Matrix<F>& A, Matrix<Base<F>>& s )
{
    DEBUG_ONLY(CSE cse("HermitianSVD"))
#if 1
    typedef Base<F> Real;
    // Grab the eigenvalues of A
    HermitianEig( uplo, A, s );

    // Set the singular values to the absolute value of the eigenvalues
    auto absLambda = []( Real sigma ) { return Abs(sigma); };
    EntrywiseMap( s, function<Real(Real)>(absLambda) );

    Sort( s, DESCENDING );
#else
    MakeHermitian( uplo, A );
    SVD( A, s );
#endif 
}

/*
   Copyright (c) 2009-2015, Jack Poulson
   All rights reserved.

   This file is part of Elemental and is under the BSD 2-Clause License, 
   which can be found in the LICENSE file in the root directory, or at 
   http://opensource.org/licenses/BSD-2-Clause
*/
#include "El.hpp"

namespace El {

template<typename T>
void Copy( const Matrix<T>& A, Matrix<T>& B )
{
    DEBUG_ONLY(CallStackEntry cse("Copy"))
    const Int height = A.Height();
    const Int width = A.Width();
    B.Resize( height, width ); 

    const Int ALDim = A.LDim();
    const Int BLDim = B.LDim();
    const T* ABuf = A.LockedBuffer();
    T* BBuf = B.Buffer();
    EL_PARALLEL_FOR
    for( Int j=0; j<width; ++j )
        MemCopy( &BBuf[j*BLDim], &ABuf[j*ALDim], height );
}

template<typename S,typename T>
void Copy( const Matrix<S>& A, Matrix<T>& B )

   This file is part of Elemental and is under the BSD 2-Clause License, 
   which can be found in the LICENSE file in the root directory, or at 
   http://opensource.org/licenses/BSD-2-Clause
*/
#include "El.hpp"

namespace El {

template<typename FDiag,typename F>
void DiagonalSolve
( LeftOrRight side, Orientation orientation,
  const Matrix<FDiag>& d, 
        Matrix<F>& A, 
  bool checkIfSingular )
{
    DEBUG_ONLY(CSE cse("DiagonalSolve"))
    const Int m = A.Height();
    const Int n = A.Width();
    const bool conj = ( orientation == ADJOINT );
    F* ABuf = A.Buffer();
    const Int ALDim = A.LDim();
    const FDiag* dBuf = d.LockedBuffer();
    if( side == LEFT )
    {
        for( Int i=0; i<m; ++i )
        {
            const F delta = ( conj ? Conj(dBuf[i]) : dBuf[i] );
            if( checkIfSingular && delta == F(0) )
                throw SingularMatrixException();
            const F deltaInv = F(1)/delta;
            for( Int j=0; j<n; ++j )

   Copyright (c) 2009-2015, Jack Poulson
   All rights reserved.

   This file is part of Elemental and is under the BSD 2-Clause License, 
   which can be found in the LICENSE file in the root directory, or at 
   http://opensource.org/licenses/BSD-2-Clause
*/
#include "El.hpp"

namespace El {

template<typename F>
pair<Base<F>,Base<F>>
ExtremalSingValEst( const SparseMatrix<F>& A, Int basisSize )
{
    DEBUG_ONLY(CSE cse("ExtremalSingValEst"))
    typedef Base<F> Real;
    Matrix<Real> T;
    ProductLanczos( A, T, basisSize );
    const Int k = T.Height();
    if( k == 0 )
        return pair<Real,Real>(0,0);

    Matrix<Real> d, dSub;
    d = GetDiagonal( T );
    dSub = GetDiagonal( T, -1 );
    
    Matrix<Real> w;
    HermitianTridiagEig( d, dSub, w, ASCENDING );
    
    pair<Real,Real> extremal;

void Zeros( AbstractBlockDistMatrix<T>& A, Int m, Int n )
{
    DEBUG_ONLY(CallStackEntry cse("Zeros"))
    A.Resize( m, n );
    Zero( A );
}

/*
   Copyright (c) 2009-2016, Jack Poulson
   All rights reserved.

   This file is part of Elemental and is under the BSD 2-Clause License, 
   which can be found in the LICENSE file in the root directory, or at 
   http://opensource.org/licenses/BSD-2-Clause
*/
#include "El.hpp"

namespace El {

template<typename F>
Base<F> LogDetDiv( UpperOrLower uplo, const Matrix<F>& A, const Matrix<F>& B )
{
    DEBUG_ONLY(CSE cse("LogDetDiv"))
    if( A.Height() != A.Width() || B.Height() != B.Width() ||
        A.Height() != B.Height() )
        LogicError("A and B must be square matrices of the same size");

    typedef Base<F> Real;
    const Int n = A.Height();

    Matrix<F> ACopy( A ), BCopy( B );
    Cholesky( uplo, ACopy );
    Cholesky( uplo, BCopy );

    if( uplo == LOWER )
    {
        Trstrm( LEFT, uplo, NORMAL, NON_UNIT, F(1), BCopy, ACopy );
    }