C++ DistMatrix::LockedMatrix方法代碼示例

inline void LocalGer
( T alpha, const DistMatrix<T,xColDist,xRowDist>& x,
           const DistMatrix<T,yColDist,yRowDist>& y,
                 DistMatrix<T,AColDist,ARowDist>& A )
{
    DEBUG_ONLY(CallStackEntry cse("LocalGer"))
    // TODO: Add error checking here
    Ger( alpha, x.LockedMatrix(), y.LockedMatrix(), A.Matrix() );
}

inline void LocalGer
( T alpha, const DistMatrix<T,xColDist,xRowDist>& x,
  const DistMatrix<T,yColDist,yRowDist>& y,
  DistMatrix<T,AColDist,ARowDist>& A )
{
#ifndef RELEASE
    CallStackEntry entry("LocalGer");
    // TODO: Add error checking here
#endif
    Ger( alpha, x.LockedMatrix(), y.LockedMatrix(), A.Matrix() );
}

const DistMatrix<T,MD,STAR>&
DistMatrix<T,MD,STAR>::operator=( const DistMatrix<T,MD,STAR>& A )
{
#ifndef RELEASE
    CallStackEntry entry("[MD,* ] = [MD,* ]");
    this->AssertNotLocked();
    this->AssertSameGrid( A.Grid() );
#endif
    if( !this->Viewing() && !this->ConstrainedColAlignment() )
    {
        this->diagPath_ = A.diagPath_;
        this->colAlignment_ = A.colAlignment_;
        if( this->Participating() )
            this->colShift_ = A.ColShift();
    }
    this->ResizeTo( A.Height(), A.Width() );

    if( this->diagPath_ == A.diagPath_ && 
        this->colAlignment_ == A.colAlignment_ )
    {
        this->matrix_ = A.LockedMatrix();
    }
    else
    {
#ifdef UNALIGNED_WARNINGS
        if( this->Grid().Rank() == 0 )
            std::cerr << "Unaligned [MD,* ] <- [MD,* ]." << std::endl;
#endif
        LogicError("Unaligned [MD,* ] = [MD,* ] not yet implemented");
    }
    return *this;
}

void RowMaxNorms
( const DistMatrix<F,U,V>& A, DistMatrix<Base<F>,U,STAR>& norms )
{
    DEBUG_CSE
    norms.AlignWith( A );
    norms.Resize( A.Height(), 1 );
    RowMaxNorms( A.LockedMatrix(), norms.Matrix() );
    AllReduce( norms, A.RowComm(), mpi::MAX );
}

inline void
DiagonalScale
( LeftOrRight side, Orientation orientation,
  const DistMatrix<typename Base<T>::type,U,V>& d, DistMatrix<T,W,Z>& X )
{
#ifndef RELEASE
    PushCallStack("DiagonalScale");
#endif
    typedef typename Base<T>::type R;

    if( side == LEFT )
    {
        if( U == W && V == STAR && d.ColAlignment() == X.ColAlignment() )
        {
            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.ColAlignment() == X.RowAlignment() )
        {
            DiagonalScale( RIGHT, orientation, d.LockedMatrix(), X.Matrix() );
        }
        else
        {
            DistMatrix<R,Z,STAR> d_Z_STAR( X.Grid() );
            d_Z_STAR = d;
            DiagonalScale
            ( RIGHT, orientation, d_Z_STAR.LockedMatrix(), X.Matrix() );
        }
    }
#ifndef RELEASE
    PopCallStack();
#endif
}

inline void
Print
( const DistMatrix<T,CIRC,CIRC>& A, std::string title="",
  std::ostream& os=std::cout )
{
#ifndef RELEASE
    CallStackEntry entry("Print");
#endif
    if( A.Grid().VCRank() == A.Root() )
        Print( A.LockedMatrix(), title, os );
}

void AllGather
( const DistMatrix<T,        U,           V   >& A,
        DistMatrix<T,Collect<U>(),Collect<V>()>& B )
{
    EL_DEBUG_CSE
    AssertSameGrids( A, B );

    const Int height = A.Height();
    const Int width = A.Width();
    B.SetGrid( A.Grid() );
    B.Resize( height, width );

    if( A.Participating() )
    {
        if( A.DistSize() == 1 )
        {
            Copy( A.LockedMatrix(), B.Matrix() );
        }
        else
        {
            const Int colStride = A.ColStride();
            const Int rowStride = A.RowStride();
            const Int distStride = colStride*rowStride;
            const Int maxLocalHeight = MaxLength(height,colStride);
            const Int maxLocalWidth = MaxLength(width,rowStride);
            const Int portionSize = mpi::Pad( maxLocalHeight*maxLocalWidth );
            vector<T> buf;
            FastResize( buf, (distStride+1)*portionSize );
            T* sendBuf = &buf[0];
            T* recvBuf = &buf[portionSize];

            // Pack
            util::InterleaveMatrix
            ( A.LocalHeight(), A.LocalWidth(),
              A.LockedBuffer(), 1, A.LDim(),
              sendBuf,          1, A.LocalHeight() );

            // Communicate
            mpi::AllGather
            ( sendBuf, portionSize, recvBuf, portionSize, A.DistComm() );

            // Unpack
            util::StridedUnpack
            ( height, width,
              A.ColAlign(), colStride,
              A.RowAlign(), rowStride,
              recvBuf, portionSize,
              B.Buffer(), B.LDim() );
        }
    }
    if( A.Grid().InGrid() && A.CrossComm() != mpi::COMM_SELF )
        El::Broadcast( B, A.CrossComm(), A.Root() );
}

void RowTwoNorms
( const DistMatrix<F,U,V>& A, DistMatrix<Base<F>,U,STAR>& norms )
{
    DEBUG_CSE
    norms.AlignWith( A );
    norms.Resize( A.Height(), 1 );
    if( A.Width() == 0 )
    {
        Zero( norms );
        return;
    }
    RowTwoNormsHelper( A.LockedMatrix(), norms.Matrix(), A.RowComm() );
}

inline void
LocalTrmm
( LeftOrRight side, UpperOrLower uplo,
  Orientation orientation, UnitOrNonUnit diag,
  T alpha, const DistMatrix<T,STAR,STAR>& A,
                 DistMatrix<T,BColDist,BRowDist>& B )
{
#ifndef RELEASE
    CallStackEntry entry("LocalTrmm");
    if( (side == LEFT && BColDist != STAR) ||
        (side == RIGHT && BRowDist != STAR) )
        LogicError
        ("Distribution of RHS must conform with that of triangle");
#endif
    Trmm
    ( side, uplo, orientation, diag, alpha, A.LockedMatrix(), B.Matrix() );
}

void GatherSubdiagonal
( const DistMatrix<F,MC,MR,BLOCK>& H,
  const IR& winInd,
        DistMatrix<Base<F>,STAR,STAR>& hSubWin )
{
    DEBUG_CSE
    const Int winSize = winInd.end - winInd.beg;
    const Int blockSize = H.BlockHeight();
    const Grid& grid = H.Grid();
    const auto& HLoc = H.LockedMatrix();
    DEBUG_ONLY(
      if( H.BlockHeight() != H.BlockWidth() )
          LogicError("Assumed square distribution blocks");
      if( H.ColCut() != H.RowCut() )
          LogicError("Assumed symmetric cuts");
      if( blockSize < 2 )
          LogicError("Assumed blocks of size at least two");
    )

void TransformRows
( const Matrix<F>& Z,
        DistMatrix<F,MC,MR,BLOCK>& H )
{
    DEBUG_CSE
    const Int height = H.Height();
    const Grid& grid = H.Grid();

    const Int blockHeight = H.BlockHeight();
    const Int firstBlockHeight = blockHeight - H.ColCut();
    if( height <= firstBlockHeight || grid.Height() == 1 )
    {
        if( grid.Row() == H.RowOwner(0) )
        {
            // This process row can locally update its portion of H
            Matrix<F> HLocCopy( H.Matrix() );
            Gemm( ADJOINT, NORMAL, F(1), Z, HLocCopy, H.Matrix() );
        }
    }
    else if( height <= firstBlockHeight + blockHeight )
    {
        const bool firstRow = H.RowOwner( 0 );
        const bool secondRow = H.RowOwner( firstBlockHeight );
        if( grid.Row() == firstRow )
        {
            // 
            // Replace H with 
            //
            //   | ZLeft, ZRight |' | HTop    |,
            //                      | HBottom |
            //
            // where HTop is owned by this process row and HBottom by the next.
            //
            auto ZLeft = Z( ALL, IR(0,firstBlockHeight) );

            // Partition space for the combined matrix
            Matrix<F> HCombine( height, H.LocalWidth() );
            auto HTop = HCombine( IR(0,firstBlockHeight), ALL );
            auto HBottom = HCombine( IR(firstBlockHeight,END), ALL );

            // Copy our portion into the combined matrix
            HTop = H.LockedMatrix();

            // Exchange the data
            El::SendRecv( HTop, HBottom, H.ColComm(), secondRow, secondRow );
            
            // Form our portion of the result
            Gemm( ADJOINT, NORMAL, F(1), ZLeft, HCombine, H.Matrix() );
        }
        else if( grid.Row() == secondRow )
        {
            // 
            // Replace H with 
            //
            //   | ZLeft, ZRight |' | HTop    |,
            //                      | HBottom |
            //
            // where HTop is owned by the previous process row and HBottom by
            // this one.
            //
            auto ZRight = Z( ALL, IR(firstBlockHeight,END) );

            // Partition space for the combined matrix
            Matrix<F> HCombine( height, H.LocalWidth() );
            auto HTop = HCombine( IR(0,firstBlockHeight), ALL );
            auto HBottom = HCombine( IR(firstBlockHeight,END), ALL );

            // Copy our portion into the combined matrix
            HBottom = H.LockedMatrix();

            // Exchange the data
            El::SendRecv( HBottom, HTop, H.ColComm(), firstRow, firstRow );
            
            // Form our portion of the result
            Gemm( ADJOINT, NORMAL, F(1), ZRight, HCombine, H.Matrix() );
        }
    }
    else
    {
        // Fall back to the entire process column interacting.
        // TODO(poulson): Only form the subset of the result that we need.
        DistMatrix<F,STAR,MR,BLOCK> H_STAR_MR( H );
        Matrix<F> HLocCopy( H_STAR_MR.Matrix() );
        Gemm( ADJOINT, NORMAL, F(1), Z, HLocCopy, H_STAR_MR.Matrix() );
        H = H_STAR_MR;
    }
}

const DistMatrix<T,STAR,STAR>&
DistMatrix<T,STAR,STAR>::operator=( const DistMatrix<T,STAR,STAR>& A )
{
#ifndef RELEASE
    CallStackEntry entry("[* ,* ] = [* ,* ]");
    this->AssertNotLocked();
#endif
    this->ResizeTo( A.Height(), A.Width() );

    if( this->Grid() == A.Grid() )
    {
        this->matrix_ = A.LockedMatrix();
    }
    else
    {
        // TODO: Remember why I wrote this...
        if( !mpi::CongruentComms( A.Grid().ViewingComm(),
                                  this->Grid().ViewingComm() ) )
            LogicError
            ("Redistributing between nonmatching grids currently requires"
             " the viewing communicators to match.");

        // Compute and allocate the amount of required memory
        Int requiredMemory = 0;
        if( A.Grid().VCRank() == 0 )
            requiredMemory += A.Height()*A.Width();
        if( this->Participating() )
            requiredMemory += A.Height()*A.Width();
        T* buffer = this->auxMemory_.Require( requiredMemory );
        Int offset = 0;
        T* sendBuf = &buffer[offset];
        if( A.Grid().VCRank() == 0 )
            offset += A.Height()*A.Width();
        T* bcastBuffer = &buffer[offset];

        // Send from the root of A to the root of this matrix's grid
        mpi::Request sendRequest;
        if( A.Grid().VCRank() == 0 )
        {
            for( Int j=0; j<A.Width(); ++j ) 
                for( Int i=0; i<A.Height(); ++i )
                    sendBuf[i+j*A.Height()] = A.GetLocal(i,j);
            const Int recvViewingRank = this->Grid().VCToViewingMap(0);
            mpi::ISend
            ( sendBuf, A.Height()*A.Width(), recvViewingRank,
              this->Grid().ViewingComm(), sendRequest );
        }

        // Receive on the root of this matrix's grid and then broadcast
        // over this matrix's owning communicator
        if( this->Participating() )
        {
            if( this->Grid().VCRank() == 0 )
            {
                const Int sendViewingRank = A.Grid().VCToViewingMap(0);
                mpi::Recv
                ( bcastBuffer, A.Height()*A.Width(), sendViewingRank,
                  this->Grid().ViewingComm() );
            }

            mpi::Broadcast
            ( bcastBuffer, A.Height()*A.Width(), 0, this->Grid().VCComm() );

            for( Int j=0; j<A.Width(); ++j )
                for( Int i=0; i<A.Height(); ++i )
                    this->SetLocal(i,j,bcastBuffer[i+j*A.Height()]);
        }

        if( A.Grid().VCRank() == 0 )
            mpi::Wait( sendRequest );
        this->auxMemory_.Release();
    }
    return *this;
}

inline void
ComposePivots
( const DistMatrix<Int,STAR,STAR>& p, 
  std::vector<Int>& image, std::vector<Int>& preimage )
{ ComposePivots( p.LockedMatrix(), image, preimage ); }

void TransformRows
( const Matrix<F>& V,
        DistMatrix<F,MC,MR,BLOCK>& A )
{
    DEBUG_CSE
    const Int height = A.Height();
    const Grid& grid = A.Grid();

    const Int blockHeight = A.BlockHeight();
    const Int firstBlockHeight = blockHeight - A.ColCut();
    if( height <= firstBlockHeight || grid.Height() == 1 )
    {
        if( grid.Row() == A.RowOwner(0) )
        {
            // This process row can locally update its portion of A
            TransformRows( V, A.Matrix() );
        }
    }
    else if( height <= firstBlockHeight + blockHeight )
    {
        const int firstRow = A.RowOwner( 0 );
        const int secondRow = A.RowOwner( firstBlockHeight );
        if( grid.Row() == firstRow )
        {
            //
            // Replace A with
            //
            //   | VLeft, VRight |' | ATop    |,
            //                      | ABottom |
            //
            // where ATop is owned by this process row and ABottom by the next.
            //
            auto VLeft = V( ALL, IR(0,firstBlockHeight) );

            // Partition space for the combined matrix
            Matrix<F> ACombine( height, A.LocalWidth() );
            auto ATop = ACombine( IR(0,firstBlockHeight), ALL );
            auto ABottom = ACombine( IR(firstBlockHeight,END), ALL );

            // Copy our portion into the combined matrix
            ATop = A.LockedMatrix();

            // Exchange the data
            El::SendRecv( ATop, ABottom, A.ColComm(), secondRow, secondRow );

            // Form our portion of the result
            Gemm( ADJOINT, NORMAL, F(1), VLeft, ACombine, A.Matrix() );
        }
        else if( grid.Row() == secondRow )
        {
            //
            // Replace A with
            //
            //   | VLeft, VRight |' | ATop    |,
            //                      | ABottom |
            //
            // where ATop is owned by the previous process row and ABottom by
            // this one.
            //
            auto VRight = V( ALL, IR(firstBlockHeight,END) );

            // Partition space for the combined matrix
            Matrix<F> ACombine( height, A.LocalWidth() );
            auto ATop = ACombine( IR(0,firstBlockHeight), ALL );
            auto ABottom = ACombine( IR(firstBlockHeight,END), ALL );

            // Copy our portion into the combined matrix
            ABottom = A.LockedMatrix();

            // Exchange the data
            El::SendRecv( ABottom, ATop, A.ColComm(), firstRow, firstRow );

            // Form our portion of the result
            Gemm( ADJOINT, NORMAL, F(1), VRight, ACombine, A.Matrix() );
        }
    }
    else
    {
        // Fall back to the entire process column interacting.
        // TODO(poulson): Only form the subset of the result that we need.
        DistMatrix<F,STAR,MR,BLOCK> A_STAR_MR( A );
        Matrix<F> ALocCopy( A_STAR_MR.Matrix() );
        Gemm( ADJOINT, NORMAL, F(1), V, ALocCopy, A_STAR_MR.Matrix() );
        A = A_STAR_MR;
    }
}

void Scatter
( const DistMatrix<T,CIRC,CIRC>& A,
        ElementalMatrix<T>& B )
{
    DEBUG_CSE
    AssertSameGrids( A, B );

    const Int m = A.Height();
    const Int n = A.Width();
    const Int colStride = B.ColStride();
    const Int rowStride = B.RowStride();
    B.Resize( m, n );
    if( B.CrossSize() != 1 || B.RedundantSize() != 1 )
    {
        // TODO:
        // Broadcast over the redundant communicator and use mpi::Translate
        // rank to determine whether a process is the root of the broadcast.
        GeneralPurpose( A, B ); 
        return;
    }

    const Int pkgSize = mpi::Pad(MaxLength(m,colStride)*MaxLength(n,rowStride));
    const Int recvSize = pkgSize;
    const Int sendSize = B.DistSize()*pkgSize;

    // Translate the root of A into the DistComm of B (if possible)
    const Int root = A.Root();
    const Int target = mpi::Translate( A.CrossComm(), root, B.DistComm() ); 
    if( target == mpi::UNDEFINED )
        return;

    if( B.DistSize() == 1 )
    {
        Copy( A.LockedMatrix(), B.Matrix() );
        return;
    }

    vector<T> buffer;
    T* recvBuf=0; // some compilers (falsely) warn otherwise
    if( A.CrossRank() == root )
    {
        FastResize( buffer, sendSize+recvSize );
        T* sendBuf = &buffer[0];
        recvBuf    = &buffer[sendSize];

        // Pack the send buffer
        copy::util::StridedPack
        ( m, n,
          B.ColAlign(), colStride,
          B.RowAlign(), rowStride,
          A.LockedBuffer(), A.LDim(),
          sendBuf,          pkgSize );

        // Scatter from the root
        mpi::Scatter
        ( sendBuf, pkgSize, recvBuf, pkgSize, target, B.DistComm() );
    }
    else
    {
        FastResize( buffer, recvSize );
        recvBuf = &buffer[0];

        // Perform the receiving portion of the scatter from the non-root
        mpi::Scatter
        ( static_cast<T*>(0), pkgSize,
          recvBuf,            pkgSize, target, B.DistComm() );
    }

    // Unpack
    copy::util::InterleaveMatrix
    ( B.LocalHeight(), B.LocalWidth(),
      recvBuf,    1, B.LocalHeight(),
      B.Buffer(), 1, B.LDim() );
}