C++ DistMultiVec::LocalHeight方法代码示例

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

void PushInto
(       DistMultiVec<Real>& x, 
  const DistMultiVec<Int>& orders, 
  const DistMultiVec<Int>& firstInds, 
  Real minDist, Int cutoff )
{
    DEBUG_ONLY(CSE cse("soc::PushInto"))

    DistMultiVec<Real> d(x.Comm());
    soc::LowerNorms( x, d, orders, firstInds, cutoff );

    const int localHeight = x.LocalHeight();
    for( Int iLoc=0; iLoc<localHeight; ++iLoc )
    {
        const Int i = x.GlobalRow(iLoc);
        const Real x0 = x.GetLocal(iLoc,0);
        const Real lowerNorm = d.GetLocal(iLoc,0);
        if( i == firstInds.GetLocal(iLoc,0) && x0-lowerNorm < minDist )
            x.UpdateLocal( iLoc, 0, minDist - (x0-lowerNorm) );
    }
}

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 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 SOCApply
( 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("SOCApply"))
    SOCDots( x, y, z, orders, firstInds );
    auto xRoots = x;
    auto yRoots = y;
    ConeBroadcast( xRoots, orders, firstInds );
    ConeBroadcast( yRoots, orders, firstInds );
    const Int localHeight = x.LocalHeight();
    for( Int iLoc=0; iLoc<localHeight; ++iLoc )
    {
        const Int i = x.GlobalRow(iLoc);
        const Int firstInd = firstInds.GetLocal(iLoc,0);
        if( i != firstInd )
            z.UpdateLocal
            ( iLoc, 0, xRoots.GetLocal(iLoc,0)*y.GetLocal(iLoc,0) +
                       yRoots.GetLocal(iLoc,0)*x.GetLocal(iLoc,0) );
    }
}

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