C++ DistMultiVec类代码示例

void SOCSquareRoot
( const DistMultiVec<Real>& x, 
        DistMultiVec<Real>& xRoot,
  const DistMultiVec<Int>& orders, 
  const DistMultiVec<Int>& firstInds, Int cutoff )
{
    DEBUG_ONLY(CSE cse("SOCSquareRoot"))

    DistMultiVec<Real> d(x.Comm());
    SOCDets( x, d, orders, firstInds );
    ConeBroadcast( d, orders, firstInds );

    auto roots = x;
    ConeBroadcast( roots, orders, firstInds );

    const Int localHeight = x.LocalHeight();
    xRoot.SetComm( x.Comm() );
    Zeros( xRoot, x.Height(), 1 );
    for( Int iLoc=0; iLoc<localHeight; ++iLoc )
    {
        const Int i = x.GlobalRow(iLoc);
        const Real x0 = roots.GetLocal(iLoc,0);
        const Real det = d.GetLocal(iLoc,0);
        const Real eta0 = Sqrt(x0+Sqrt(det))/Sqrt(Real(2));
        if( i == firstInds.GetLocal(iLoc,0) )
            xRoot.SetLocal( iLoc, 0, eta0 );
        else
            xRoot.SetLocal( iLoc, 0, x.GetLocal(iLoc,0)/(2*eta0) );
    }
}

Int Degree( const DistMultiVec<Int>& firstInds )
{
    DEBUG_CSE
    Int localDegree = 0;
    const Int localHeight = firstInds.LocalHeight();
    auto& firstIndsLoc = firstInds.LockedMatrix();
    for( Int iLoc=0; iLoc<localHeight; ++iLoc )
    {
        const Int i = firstInds.GlobalRow(iLoc);
        if( i == firstIndsLoc(iLoc) )
            ++localDegree;
    }
    return mpi::AllReduce( localDegree, firstInds.Comm() );
}

void PushPairInto
(       DistMultiVec<Real>& s, 
        DistMultiVec<Real>& z,
  const DistMultiVec<Real>& w,
  const DistMultiVec<Int>& orders, 
  const DistMultiVec<Int>& firstInds, 
  Real wMaxNormLimit, Int cutoff )
{
    DEBUG_ONLY(CSE cse("soc::PushPairInto"))

    DistMultiVec<Real> sLower(s.Comm()), zLower(z.Comm());
    soc::LowerNorms( s, sLower, orders, firstInds, cutoff );
    soc::LowerNorms( z, zLower, orders, firstInds, cutoff );

    const int localHeight = s.LocalHeight();
    for( Int iLoc=0; iLoc<localHeight; ++iLoc )
    {
        const Int i = s.GlobalRow(iLoc);
        const Real w0 = w.GetLocal(iLoc,0);
        if( i == firstInds.GetLocal(iLoc,0) && w0 > wMaxNormLimit )
        {
            // TODO: Switch to a non-adhoc modification     
            s.UpdateLocal( iLoc, 0, Real(1)/wMaxNormLimit );
            z.UpdateLocal( iLoc, 0, Real(1)/wMaxNormLimit );
        }
    }
}

void SymmetricRuizEquil
( DistSparseMatrix<Field>& A,
  DistMultiVec<Base<Field>>& d,
  Int maxIter, bool progress )
{
    EL_DEBUG_CSE
    typedef Base<Field> Real;
    const Int n = A.Height();
    const Grid& grid = A.Grid();
    d.SetGrid( grid );
    Ones( d, n, 1 );

    DistMultiVec<Real> scales(grid);
    const Int indent = PushIndent();
    for( Int iter=0; iter<maxIter; ++iter )
    {
        // Rescale the columns (and rows)
        // ------------------------------
        ColumnMaxNorms( A, scales );
        EntrywiseMap( scales, MakeFunction(DampScaling<Real>) );
        EntrywiseMap( scales, MakeFunction(SquareRootScaling<Real>) );
        DiagonalScale( LEFT, NORMAL, scales, d );
        SymmetricDiagonalSolve( scales, A );
    }
    SetIndent( indent );
}

void SymmetricRuizEquil
( DistSparseMatrix<F>& A, 
  DistMultiVec<Base<F>>& d, 
  Int maxIter, bool progress )
{
    DEBUG_CSE
    typedef Base<F> Real;
    const Int n = A.Height();
    mpi::Comm comm = A.Comm();
    d.SetComm( comm );
    Ones( d, n, 1 );

    DistMultiVec<Real> scales(comm);
    const Int indent = PushIndent();
    for( Int iter=0; iter<maxIter; ++iter )
    {
        // Rescale the columns (and rows)
        // ------------------------------
        ColumnMaxNorms( A, scales );
        EntrywiseMap( scales, function<Real(Real)>(DampScaling<Real>) );
        EntrywiseMap( scales, function<Real(Real)>(SquareRootScaling<Real>) );
        DiagonalScale( LEFT, NORMAL, scales, d );
        SymmetricDiagonalSolve( scales, A );
    }
    SetIndent( indent );
}

inline void
DistNodalMultiVec<F>::Pull
( const DistMap& inverseMap, const DistSymmInfo& info,
  const DistMultiVec<F>& X )
{
    DEBUG_ONLY(CallStackEntry cse("DistNodalMultiVec::Pull"))
    height_ = X.Height();
    width_ = X.Width();

    // Traverse our part of the elimination tree to see how many indices we need
    int numRecvInds=0;
    const int numLocal = info.localNodes.size();
    for( int s=0; s<numLocal; ++s )
        numRecvInds += info.localNodes[s].size;
    const int numDist = info.distNodes.size();
    for( int s=1; s<numDist; ++s )
        numRecvInds += info.distNodes[s].multiVecMeta.localSize;
    
    // Fill the set of indices that we need to map to the original ordering
    int off=0;
    std::vector<int> mappedInds( numRecvInds );
    for( int s=0; s<numLocal; ++s )
    {
        const SymmNodeInfo& nodeInfo = info.localNodes[s];
        for( int t=0; t<nodeInfo.size; ++t )
            mappedInds[off++] = nodeInfo.off+t;
    }
    for( int s=1; s<numDist; ++s )
    {
        const DistSymmNodeInfo& nodeInfo = info.distNodes[s];
        const Grid& grid = *nodeInfo.grid;
        const int gridSize = grid.Size();
        const int gridRank = grid.VCRank();
        const int alignment = 0;
        const int shift = Shift( gridRank, alignment, gridSize );
        for( int t=shift; t<nodeInfo.size; t+=gridSize )
            mappedInds[off++] = nodeInfo.off+t;
    }
    DEBUG_ONLY(
        if( off != numRecvInds )
            LogicError("mappedInds was filled incorrectly");
    )

void QP
( const DistSparseMatrix<Real>& A, 
  const DistMultiVec<Real>& B, 
        DistMultiVec<Real>& X, 
  const qp::direct::Ctrl<Real>& ctrl )
{
    DEBUG_CSE

    const Int m = A.Height();
    const Int n = A.Width();
    const Int k = B.Width();
    mpi::Comm comm = A.Comm();
    DistSparseMatrix<Real> Q(comm), AHat(comm);
    DistMultiVec<Real> bHat(comm), c(comm);

    Herk( LOWER, ADJOINT, Real(1), A, Q );
    MakeHermitian( LOWER, Q );
    Zeros( AHat, 0, n );
    Zeros( bHat, 0, 1 );
    Zeros( X,    n, k );

    DistMultiVec<Real> q(comm), y(comm), z(comm);
    auto& qLoc = q.Matrix();
    auto& XLoc = X.Matrix();
    auto& BLoc = B.LockedMatrix();
    for( Int j=0; j<k; ++j )
    {
        auto xLoc = XLoc( ALL, IR(j) );
        auto bLoc = BLoc( ALL, IR(j) );

        Zeros( c, n, 1 );
        Zeros( q, m, 1 );
        qLoc = bLoc;
        Multiply( ADJOINT, Real(-1), A, q, Real(0), c );

        Zeros( q, n, 1 );
        qLoc = xLoc;
        El::QP( Q, AHat, bHat, c, q, y, z, ctrl );
        xLoc = qLoc;
    }
}

Real ComplementRatio
( const DistMultiVec<Real>& s,
  const DistMultiVec<Real>& z )
{
    DEBUG_CSE
    const Int localHeight = s.LocalHeight();
    const Real* sBuf = s.LockedMatrix().LockedBuffer();
    const Real* zBuf = z.LockedMatrix().LockedBuffer();

    Real maxLocProd = 0;
    for( Int iLoc=0; iLoc<localHeight; ++iLoc )
        maxLocProd = Max( sBuf[iLoc]*zBuf[iLoc], maxLocProd );
    const Real maxProd = mpi::AllReduce( maxLocProd, mpi::MAX, s.Comm() );

    Real minLocProd = maxProd;
    for( Int iLoc=0; iLoc<localHeight; ++iLoc )
        minLocProd = Min( sBuf[iLoc]*zBuf[iLoc], minLocProd );
    const Real minProd = mpi::AllReduce( minLocProd, mpi::MIN, s.Comm() );

    return maxProd/minProd;
}

void EntrywiseMap
( const DistMultiVec<S>& A,
        DistMultiVec<T>& B,
        function<T(S)> func )
{
    DEBUG_CSE
    B.SetComm( A.Comm() );
    B.Resize( A.Height(), A.Width() );
    EntrywiseMap( A.LockedMatrix(), B.Matrix(), func );
}

void SOCApply
( const DistMultiVec<Real>& x, 
        DistMultiVec<Real>& y,
  const DistMultiVec<Int>& orders, 
  const DistMultiVec<Int>& firstInds, Int cutoff )
{
    DEBUG_ONLY(CSE cse("SOCApply"))
    // TODO?: Optimize
    DistMultiVec<Real> z(x.Comm());
    SOCApply( x, y, z, orders, firstInds, cutoff );
    y = z;
}

void RowMaxNorms( const DistSparseMatrix<F>& A, DistMultiVec<Base<F>>& norms )
{
    DEBUG_CSE
    typedef Base<F> Real;
    const Int localHeight = A.LocalHeight();
    const F* valBuf = A.LockedValueBuffer();
    const Int* offsetBuf = A.LockedOffsetBuffer();

    norms.SetComm( A.Comm() );
    norms.Resize( A.Height(), 1 );
    auto& normsLoc = norms.Matrix();
    for( Int iLoc=0; iLoc<localHeight; ++iLoc )
    {
        Real rowMax = 0;
        const Int offset = offsetBuf[iLoc];
        const Int numConn = offsetBuf[iLoc+1] - offset;
        for( Int e=offset; e<offset+numConn; ++e )
            rowMax = Max(rowMax,Abs(valBuf[e]));
        normsLoc(iLoc) = rowMax;
    }
}

void Apply
( const DistMultiVec<Real>& x, 
  const DistMultiVec<Real>& y,
        DistMultiVec<Real>& z,
  const DistMultiVec<Int>& orders, 
  const DistMultiVec<Int>& firstInds,
  Int cutoff )
{
    DEBUG_ONLY(CSE cse("soc::Apply"))
    soc::Dots( x, y, z, orders, firstInds );
    auto xRoots = x;
    auto yRoots = y;
    cone::Broadcast( xRoots, orders, firstInds );
    cone::Broadcast( yRoots, orders, firstInds );

    const Int firstLocalRow = x.FirstLocalRow();
    const Int localHeight = x.LocalHeight();
    const Real* xBuf     = x.LockedMatrix().LockedBuffer();
    const Real* xRootBuf = xRoots.LockedMatrix().LockedBuffer();
    const Real* yBuf     = y.LockedMatrix().LockedBuffer();
    const Real* yRootBuf = yRoots.LockedMatrix().LockedBuffer();
          Real* zBuf     = z.Matrix().Buffer();
    const Int* firstIndBuf = firstInds.LockedMatrix().LockedBuffer();

    for( Int iLoc=0; iLoc<localHeight; ++iLoc )
    {
        const Int i = iLoc + firstLocalRow;
        const Int firstInd = firstIndBuf[iLoc];
        if( i != firstInd )
            zBuf[iLoc] += xRootBuf[iLoc]*yBuf[iLoc] + yRootBuf[iLoc]*xBuf[iLoc];
    }
}

void RowTwoNorms( const DistSparseMatrix<F>& A, DistMultiVec<Base<F>>& norms )
{
    DEBUG_CSE
    typedef Base<F> Real;
    const Int localHeight = A.LocalHeight();
    const F* valBuf = A.LockedValueBuffer();
    const Int* offsetBuf = A.LockedOffsetBuffer();

    norms.SetComm( A.Comm() );
    norms.Resize( A.Height(), 1 );
    auto& normLoc = norms.Matrix();
    for( Int iLoc=0; iLoc<localHeight; ++iLoc )
    {
        Real scale = 0;
        Real scaledSquare = 1;
        const Int offset = offsetBuf[iLoc];
        const Int numConn = offsetBuf[iLoc+1] - offset;
        for( Int e=offset; e<offset+numConn; ++e )
            UpdateScaledSquare( valBuf[e], scale, scaledSquare );
        normLoc(iLoc) = scale*Sqrt(scaledSquare);
    }
}

void GetMappedDiagonal
( const DistSparseMatrix<T>& A,
        DistMultiVec<S>& d,
        function<S(const T&)> func,
        Int offset )
{
    EL_DEBUG_CSE
    const Int m = A.Height();
    const Int n = A.Width();
    const T* valBuf = A.LockedValueBuffer();
    const Int* colBuf = A.LockedTargetBuffer();

    if( m != n )
        LogicError("DistSparseMatrix GetMappedDiagonal assumes square matrix");
    if( offset != 0 )
        LogicError("DistSparseMatrix GetMappedDiagonal assumes offset=0");

    d.SetGrid( A.Grid() );
    d.Resize( El::DiagonalLength(m,n,offset), 1 );
    Fill( d, S(1) );

    S* dBuf = d.Matrix().Buffer();
    const Int dLocalHeight = d.LocalHeight();
    for( Int iLoc=0; iLoc<dLocalHeight; ++iLoc )
    {
        const Int i = d.GlobalRow(iLoc);
        const Int thisOff = A.RowOffset(iLoc);
        const Int nextOff = A.RowOffset(iLoc+1);
        auto it = std::lower_bound( colBuf+thisOff, colBuf+nextOff, i );
        if( *it == i )
        {
            const Int e = it-colBuf;
            dBuf[iLoc] = func(valBuf[e]);
        }
        else
            dBuf[iLoc] = func(0);
    }
}

void NesterovTodd
( const DistMultiVec<Real>& s, 
  const DistMultiVec<Real>& z,
        DistMultiVec<Real>& w )
{
    DEBUG_CSE
    w.SetComm( s.Comm() );
    w.Resize( s.Height(), 1 );
    const Real* sBuf = s.LockedMatrix().LockedBuffer();
    const Real* zBuf = z.LockedMatrix().LockedBuffer();
          Real* wBuf = w.Matrix().Buffer();

    const Int localHeight = w.LocalHeight();
    for( Int iLoc=0; iLoc<localHeight; ++iLoc )
        wBuf[iLoc] = Sqrt(sBuf[iLoc]/zBuf[iLoc]);
}