C++ DistSparseMatrix::Comm方法代码示例

void StackedRuizEquil
( DistSparseMatrix<Field>& A,
  DistSparseMatrix<Field>& B,
  DistMultiVec<Base<Field>>& dRowA,
  DistMultiVec<Base<Field>>& dRowB,
  DistMultiVec<Base<Field>>& dCol,
  bool progress )
{
    EL_DEBUG_CSE
    typedef Base<Field> Real;
    const Int mA = A.Height();
    const Int mB = B.Height();
    const Int n = A.Width();
    mpi::Comm comm = A.Comm();
    dRowA.SetComm( comm );
    dRowB.SetComm( comm );
    dCol.SetComm( comm );
    Ones( dRowA, mA, 1 );
    Ones( dRowB, mB, 1 );
    Ones( dCol, n, 1 );

    // TODO(poulson): Expose to control structure
    // For, simply hard-code a small number of iterations
    const Int maxIter = 4;

    DistMultiVec<Real> scales(comm), maxAbsValsB(comm);
    auto& scalesLoc = scales.Matrix();
    auto& maxAbsValsBLoc = maxAbsValsB.Matrix();
    const Int localHeight = scalesLoc.Height();
    const Int indent = PushIndent();
    for( Int iter=0; iter<maxIter; ++iter )
    {
        // Rescale the columns
        // -------------------
        ColumnMaxNorms( A, scales );
        ColumnMaxNorms( B, maxAbsValsB );
        for( Int jLoc=0; jLoc<localHeight; ++jLoc )
            scalesLoc(jLoc) = Max(scalesLoc(jLoc),maxAbsValsBLoc(jLoc));
        EntrywiseMap( scales, MakeFunction(DampScaling<Real>) );
        DiagonalScale( LEFT, NORMAL, scales, dCol );
        DiagonalSolve( RIGHT, NORMAL, scales, A );
        DiagonalSolve( RIGHT, NORMAL, scales, B );

        // Rescale the rows
        // ----------------
        RowMaxNorms( A, scales );
        EntrywiseMap( scales, MakeFunction(DampScaling<Real>) );
        DiagonalScale( LEFT, NORMAL, scales, dRowA );
        DiagonalSolve( LEFT, NORMAL, scales, A );

        RowMaxNorms( B, scales );
        EntrywiseMap( scales, MakeFunction(DampScaling<Real>) );
        DiagonalScale( LEFT, NORMAL, scales, dRowB );
        DiagonalSolve( LEFT, NORMAL, scales, B );
    }
    SetIndent( indent );
}

Int ZeroNorm( const DistSparseMatrix<T>& A, Base<T> tol )
{
    DEBUG_ONLY(CSE cse("ZeroNorm"))
    Int numNonzeros = 0;
    const Int numLocalEntries = A.NumLocalEntries();
    for( Int k=0; k<numLocalEntries; ++k )
        if( Abs(A.Value(k)) > tol )
            ++numNonzeros;
    return mpi::AllReduce( numNonzeros, A.Comm() );
}

Int NumOutside( const DistSparseMatrix<Real>& A )
{
    EL_DEBUG_CSE
    const Int numLocalEntries = A.NumLocalEntries();
    const Real* valBuf = A.LockedValueBuffer();

    Int numLocalNonPos = 0;
    for( Int k=0; k<numLocalEntries; ++k )
        if( valBuf[k] <= Real(0) )
            ++numLocalNonPos;

    return mpi::AllReduce( numLocalNonPos, A.Comm() );
}

void AugmentedKKT
( const DistSparseMatrix<Real>& A,
        Real gamma,
        Real delta,
  const DistMultiVec<Real>& x,
  const DistMultiVec<Real>& z,
        DistSparseMatrix<Real>& J,
  bool onlyLower )
{
    EL_DEBUG_CSE
    const Int n = A.Width();
    DistSparseMatrix<Real> Q(A.Comm());
    Zeros( Q, n, n );
    qp::direct::AugmentedKKT( Q, A, gamma, delta, x, z, J, onlyLower );
}

Base<Field>
MinAbsNonzero( const DistSparseMatrix<Field>& A, Base<Field> upperBound )
{
    EL_DEBUG_CSE
    typedef Base<Field> Real;
    const Int numEntries = A.NumLocalEntries();
    Real minLocAbs = upperBound;
    for( Int e=0; e<numEntries; ++e )
    {
        const Real absVal = Abs(A.Value(e));
        if( absVal > Real(0) )
            minLocAbs = Min(minLocAbs,absVal);
    }
    return mpi::AllReduce( minLocAbs, mpi::MIN, A.Comm() );
}

Base<Field> ProductLanczosDecomp
( const DistSparseMatrix<Field>& A,
        DistMultiVec<Field>& V,
        AbstractDistMatrix<Base<Field>>& T,
        DistMultiVec<Field>& v,
        Int basisSize )
{
    EL_DEBUG_CSE
    const Int m = A.Height();
    const Int n = A.Width();
    mpi::Comm comm = A.Comm();

    DistMultiVec<Field> S(comm);

    // Cache the adjoint
    // -----------------
    DistSparseMatrix<Field> AAdj(comm);
    Adjoint( A, AAdj );

    if( m >= n )
    {
        auto applyA =
          [&]( const DistMultiVec<Field>& X, DistMultiVec<Field>& Y )
          {
              Zeros( S, m, X.Width() );
              Multiply( NORMAL, Field(1), A,    X, Field(0), S );
              Zeros( Y, n, X.Width() );
              Multiply( NORMAL, Field(1), AAdj, S, Field(0), Y );
          };
        return LanczosDecomp( n, applyA, V, T, v, basisSize );
    }
    else
    {
        auto applyA =
          [&]( const DistMultiVec<Field>& X, DistMultiVec<Field>& Y )
          {
              Zeros( S, n, X.Width() );
              Multiply( NORMAL, Field(1), AAdj, X, Field(0), S );
              Zeros( Y, m, X.Width() );
              Multiply( NORMAL, Field(1), A,    S, Field(0), Y );
          };
        return LanczosDecomp( m, applyA, V, T, v, basisSize );
    }
}

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;
    }
}

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;
    }
}

pair<Base<F>,Base<F>>
Helper( const DistSparseMatrix<F>& A, Int basisSize )
{
    typedef Base<F> Real;
    Grid grid( A.Comm() );
    DistMatrix<Real,STAR,STAR> T(grid);
    ProductLanczos( A, T, basisSize );
    const Int k = T.Height();
    if( k == 0 )
        return pair<Real,Real>(0,0);

    auto d = GetDiagonal( T.Matrix() );
    auto dSub = GetDiagonal( T.Matrix(), -1 );
    
    Matrix<Real> w;
    HermitianTridiagEig( d, dSub, w, ASCENDING );
    
    pair<Real,Real> extremal;
    extremal.first = Sqrt( Max(w.Get(0,0),Real(0)) );
    extremal.second = Sqrt( Max(w.Get(k-1,0),Real(0)) );
    return extremal;
}

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 IPM
( const DistSparseMatrix<Real>& A,
  const DistMultiVec<Real>& b, 
        Real lambda,
        DistMultiVec<Real>& x,
  const qp::affine::Ctrl<Real>& ctrl )
{
    DEBUG_CSE
    const Int m = A.Height();
    const Int n = A.Width();
    mpi::Comm comm = A.Comm();
    DistSparseMatrix<Real> Q(comm), AHat(comm), G(comm);
    DistMultiVec<Real> c(comm), h(comm);

    // Q := | 0 0 0 |
    //      | 0 0 0 |
    //      | 0 0 I |
    // ==============
    Zeros( Q, 2*n+m, 2*n+m );
    {
        Int numLocalUpdates = 0;
        for( Int iLoc=0; iLoc<Q.LocalHeight(); ++iLoc )
            if( Q.GlobalRow(iLoc) >= 2*n )
                ++numLocalUpdates;
        Q.Reserve( numLocalUpdates );
        for( Int iLoc=0; iLoc<Q.LocalHeight(); ++iLoc )
            if( Q.GlobalRow(iLoc) >= 2*n )
                Q.QueueLocalUpdate( iLoc, Q.GlobalRow(iLoc), Real(1) );
        Q.ProcessLocalQueues();
    }

    // c := lambda*[1;1;0]
    // ===================
    Zeros( c, 2*n+m, 1 );
    auto& cLoc = c.Matrix();
    for( Int iLoc=0; iLoc<c.LocalHeight(); ++iLoc )
        if( c.GlobalRow(iLoc) < 2*n )
            cLoc(iLoc) = lambda;

    // \hat A := [A, -A, I]
    // ====================
    // NOTE: Since A and \hat A are the same height and each distributed within
    //       columns, it is possible to form \hat A from A without communication
    const Int numLocalEntriesA = A.NumLocalEntries();
    Zeros( AHat, m, 2*n+m );
    AHat.Reserve( 2*numLocalEntriesA+AHat.LocalHeight() );
    for( Int e=0; e<numLocalEntriesA; ++e )
    {
        AHat.QueueUpdate( A.Row(e), A.Col(e),    A.Value(e) );
        AHat.QueueUpdate( A.Row(e), A.Col(e)+n, -A.Value(e) );
    }
    for( Int iLoc=0; iLoc<AHat.LocalHeight(); ++iLoc )
    {
        const Int i = AHat.GlobalRow(iLoc);
        AHat.QueueLocalUpdate( iLoc, i+2*n, Real(1) );
    }
    AHat.ProcessLocalQueues();

    // G := | -I  0 0 |
    //      |  0 -I 0 |
    // ================
    Zeros( G, 2*n, 2*n+m );
    G.Reserve( G.LocalHeight() );
    for( Int iLoc=0; iLoc<G.LocalHeight(); ++iLoc )
    {
        const Int i = G.GlobalRow(iLoc);
        G.QueueLocalUpdate( iLoc, i, Real(-1) );
    }
    G.ProcessLocalQueues();

    // h := 0
    // ======
    Zeros( h, 2*n, 1 );

    // Solve the affine QP
    // ===================
    DistMultiVec<Real> xHat(comm), y(comm), z(comm), s(comm);
    QP( Q, AHat, G, b, c, h, xHat, y, z, s, ctrl );

    // x := u - v
    // ==========
    Zeros( x, n, 1 );
    Int numRemoteUpdates = 0;
    for( Int iLoc=0; iLoc<xHat.LocalHeight(); ++iLoc )
        if( xHat.GlobalRow(iLoc) < 2*n )
            ++numRemoteUpdates; 
        else
            break;
    x.Reserve( numRemoteUpdates );
    auto& xHatLoc = xHat.LockedMatrix();
    for( Int iLoc=0; iLoc<xHat.LocalHeight(); ++iLoc )
    {
        const Int i = xHat.GlobalRow(iLoc);
        if( i < n )
            x.QueueUpdate( i, 0, xHatLoc(iLoc) );
        else if( i < 2*n )
            x.QueueUpdate( i-n, 0, -xHatLoc(iLoc) );
        else
            break;
    }
//.........这里部分代码省略.........

void RLS
( const DistSparseMatrix<Real>& A, 
  const DistMultiVec<Real>& b, 
        Real rho,
        DistMultiVec<Real>& x, 
  const socp::affine::Ctrl<Real>& ctrl )
{
    DEBUG_CSE
    const Int m = A.Height();
    const Int n = A.Width();
    mpi::Comm comm = A.Comm();

    DistMultiVec<Int> orders(comm), firstInds(comm);
    Zeros( orders, m+n+3, 1 );
    Zeros( firstInds, m+n+3, 1 );
    {
        const Int localHeight = orders.LocalHeight();
        for( Int iLoc=0; iLoc<localHeight; ++iLoc )
        {
            const Int i = orders.GlobalRow(iLoc);
            if( i < m+1 )
            {
                orders.SetLocal( iLoc, 0, m+1 );
                firstInds.SetLocal( iLoc, 0, 0 );
            }
            else
            {
                orders.SetLocal( iLoc, 0, n+2 ); 
                firstInds.SetLocal( iLoc, 0, m+1 );
            }
        }
    }

    // G := | -1  0  0 |
    //      |  0  0  A |
    //      |  0 -1  0 |
    //      |  0  0 -I |
    //      |  0  0  0 |
    DistSparseMatrix<Real> G(comm);
    {
        Zeros( G, m+n+3, n+2 );

        // Count the number of entries of G to reserve
        // -------------------------------------------
        Int numLocalUpdates = 0;
        const Int localHeight = G.LocalHeight();
        for( Int iLoc=0; iLoc<localHeight; ++iLoc )
        {
            const Int i = G.GlobalRow(iLoc);
            if( i == 0 || i == m+1 || (i>m+1 && i<m+n+2) )
                ++numLocalUpdates;
        }
        const Int numEntriesA = A.NumLocalEntries();
        G.Reserve( numLocalUpdates+numEntriesA, numEntriesA );

        // Queue the local updates
        // ----------------------- 
        for( Int iLoc=0; iLoc<localHeight; ++iLoc )
        {
            const Int i = G.GlobalRow(iLoc);
            if( i == 0 )
                G.QueueLocalUpdate( iLoc, 0, -1 );
            else if( i == m+1 )
                G.QueueLocalUpdate( iLoc, 1, -1 );
            else if( i > m+1 && i < m+n+2 )
                G.QueueLocalUpdate( iLoc, i-m, -1 );
        }

        // Queue the remote updates
        // ------------------------
        for( Int e=0; e<numEntriesA; ++e )
            G.QueueUpdate( A.Row(e)+1, A.Col(e)+2, A.Value(e) );

        G.ProcessQueues();
    }

    // h := | 0 |
    //      | b |
    //      | 0 |
    //      | 0 |
    //      | 1 |
    DistMultiVec<Real> h(comm);
    Zeros( h, m+n+3, 1 ); 
    auto& bLoc = b.LockedMatrix();
    {
        const Int bLocalHeight = b.LocalHeight();
        h.Reserve( bLocalHeight );
        for( Int iLoc=0; iLoc<bLocalHeight; ++iLoc )
            h.QueueUpdate( b.GlobalRow(iLoc)+1, 0, bLoc(iLoc) );
        h.ProcessQueues();
    }
    h.Set( END, 0, 1 );

    // c := [1; rho; 0]
    DistMultiVec<Real> c(comm);
    Zeros( c, n+2, 1 );
    c.Set( 0, 0, 1 );
    c.Set( 1, 0, rho );

    DistSparseMatrix<Real> AHat(comm);
//.........这里部分代码省略.........

void Tikhonov
( Orientation orientation,
  const DistSparseMatrix<F>& A,
  const DistMultiVec<F>& B,
  const DistSparseMatrix<F>& G,
        DistMultiVec<F>& X, 
  const LeastSquaresCtrl<Base<F>>& ctrl )
{
    DEBUG_CSE
    mpi::Comm comm = A.Comm();
    
    // Explicitly form W := op(A)
    // ==========================
    DistSparseMatrix<F> W(comm);
    if( orientation == NORMAL )
        W = A;
    else if( orientation == TRANSPOSE )
        Transpose( A, W );
    else
        Adjoint( A, W );

    const Int m = W.Height();
    const Int n = W.Width();
    const Int numRHS = B.Width();

    // Embed into a higher-dimensional problem via appending regularization
    // ====================================================================
    DistSparseMatrix<F> WEmb(comm);
    if( m >= n )
        VCat( W, G, WEmb ); 
    else
        HCat( W, G, WEmb );

    DistMultiVec<F> BEmb(comm);
    Zeros( BEmb, WEmb.Height(), numRHS );
    if( m >= n )
    {
        // BEmb := [B; 0]
        // --------------
        const Int mLocB = B.LocalHeight();
        BEmb.Reserve( mLocB*numRHS );
        for( Int iLoc=0; iLoc<mLocB; ++iLoc )
        {
            const Int i = B.GlobalRow(iLoc);
            for( Int j=0; j<numRHS; ++j )
                BEmb.QueueUpdate( i, j, B.GetLocal(iLoc,j) );
        }
        BEmb.ProcessQueues();
    }
    else
        BEmb = B;

    // Solve the higher-dimensional problem
    // ====================================
    DistMultiVec<F> XEmb(comm);
    LeastSquares( NORMAL, WEmb, BEmb, XEmb, ctrl );

    // Extract the solution
    // ====================
    if( m >= n )
        X = XEmb;
    else
        GetSubmatrix( XEmb, IR(0,n), IR(0,numRHS), X );
}