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

/*
 * Distributes A in such a way that
 *   Layer 0 <- A(:, 0:(n/h - 1))
 *   Layer 1 <- A(:, (n/h):(2n/h - 1))
 *     .
 *     .
 *     .
 *   Layer h-1 <- A(:, ((h-1)n/h):n)
 */
void DistributeCols
( const mpi::Comm& depthComm,
  const DistMatrix<double,MC,MR>& A, 
        DistMatrix<double,MC,MR>& B )
{
    const Grid& meshGrid = A.Grid();
    const int meshSize = meshGrid.Size();
    const int depthSize = mpi::CommSize( depthComm );
    const int depthRank = mpi::CommRank( depthComm );

    const int sendCount = A.LocalHeight()*A.LocalWidth();
    const int recvCount = sendCount / depthSize;

    // For now, we will make B as large as A...
    // TODO: NOT DO THIS
    if( A.LocalHeight() != A.LocalLDim() )
        throw std::logic_error("Local height did not match local ldim");
    B.Empty();
    B.AlignWith( A );
    Zeros( A.Height(), A.Width(), B );

    // Scatter
    const int localColOffset = (A.LocalWidth()/depthSize)*depthRank;
    mpi::Scatter
    ( A.LockedLocalBuffer(), recvCount, 
      B.LocalBuffer(0,localColOffset), recvCount, 0, depthComm );
}

inline void
BinaryFlat
( DistMatrix<T,U,V>& A, Int height, Int width, const std::string filename )
{
    DEBUG_ONLY(CallStackEntry cse("read::BinaryFlat"))
    std::ifstream file( filename.c_str(), std::ios::binary );
    if( !file.is_open() )
        RuntimeError("Could not open ",filename);

    const Int numBytes = FileSize( file );
    const Int numBytesExp = height*width*sizeof(T);
    if( numBytes != numBytesExp )
        RuntimeError
        ("Expected file to be ",numBytesExp," bytes but found ",numBytes);

    A.Resize( height, width );
    if( U == A.UGath && V == A.VGath )
    {
        if( A.CrossRank() == A.Root() )
        {
            if( A.Height() == A.LDim() )
                file.read( (char*)A.Buffer(), height*width*sizeof(T) );
            else
                for( Int j=0; j<width; ++j )
                    file.read( (char*)A.Buffer(0,j), height*sizeof(T) );
        }
    }
    else if( U == A.UGath )
    {
        const Int localWidth = A.LocalWidth();
        for( Int jLoc=0; jLoc<localWidth; ++jLoc )
        {
            const Int j = A.GlobalCol(jLoc);
            const Int localIndex = j*height;
            const std::streamoff pos = localIndex*sizeof(T);
            file.seekg( pos );
            file.read( (char*)A.Buffer(0,jLoc), height*sizeof(T) );
        }
    }
    else
    {
        const Int localHeight = A.LocalHeight();
        const Int localWidth = A.LocalWidth();
        for( Int jLoc=0; jLoc<localWidth; ++jLoc )
        {
            const Int j = A.GlobalCol(jLoc);
            for( Int iLoc=0; iLoc<localHeight; ++iLoc )
            {
                const Int i = A.GlobalRow(iLoc);
                const Int localIndex = i+j*height;
                const std::streamoff pos = localIndex*sizeof(T);
                file.seekg( pos );
                file.read( (char*)A.Buffer(iLoc,jLoc), sizeof(T) );
            }
        }
    }
}

inline void
MakeHilbert( DistMatrix<F,U,V>& A )
{
#ifndef RELEASE
    PushCallStack("MakeHilbert");
#endif
    const int m = A.Height();
    const int n = A.Width();
    if( m != n )
        throw std::logic_error("Cannot make a non-square matrix Hilbert");

    const F one = static_cast<F>(1);
    const int localHeight = A.LocalHeight();
    const int localWidth = A.LocalWidth();
    const int colShift = A.ColShift();
    const int rowShift = A.RowShift();
    const int colStride = A.ColStride();
    const int rowStride = A.RowStride();
    for( int jLocal=0; jLocal<localWidth; ++jLocal )
    {
        const int j = rowShift + jLocal*rowStride;
        for( int iLocal=0; iLocal<localHeight; ++iLocal )
        {
            const int i = colShift + iLocal*colStride;
            A.SetLocalEntry( iLocal, jLocal, one/(i+j+1) );
        }
    }
#ifndef RELEASE
    PopCallStack();
#endif
}

inline void
Hankel( int m, int n, const std::vector<T>& a, DistMatrix<T,U,V>& A )
{
#ifndef RELEASE
    PushCallStack("Hankel");
#endif
    const int length = m+n-1;
    if( a.size() != (unsigned)length )
        throw std::logic_error("a was the wrong size");
    A.ResizeTo( m, n );

    const int localHeight = A.LocalHeight();
    const int localWidth = A.LocalWidth();
    const int colShift = A.ColShift();
    const int rowShift = A.RowShift();
    const int colStride = A.ColStride();
    const int rowStride = A.RowStride();
    for( int jLocal=0; jLocal<localWidth; ++jLocal )
    {
        const int j = rowShift + jLocal*rowStride;
        for( int iLocal=0; iLocal<localHeight; ++iLocal )
        {
            const int i = colShift + iLocal*colStride;
            A.SetLocal( iLocal, jLocal, a[i+j] );
        }
    }
#ifndef RELEASE
    PopCallStack();
#endif
}

inline void
MakeJordan( DistMatrix<T,U,V>& J, T lambda )
{
    DEBUG_ONLY(CallStackEntry cse("MakeJordan"))
    Zero( J.Matrix() );

    const Int localHeight = J.LocalHeight();
    const Int localWidth = J.LocalWidth();
    const Int colShift = J.ColShift();
    const Int rowShift = J.RowShift();
    const Int colStride = J.ColStride();
    const Int rowStride = J.RowStride();
    for( Int jLoc=0; jLoc<localWidth; ++jLoc )
    {
        const Int j = rowShift + jLoc*rowStride;
        for( Int iLoc=0; iLoc<localHeight; ++iLoc )
        {
            const Int i = colShift + iLoc*colStride;
            if( i == j )
                J.SetLocal( iLoc, jLoc, lambda );
            else if( i == j-1 )
                J.SetLocal( iLoc, jLoc, T(1) );
        }
    }
}

inline void
MakeLegendre( DistMatrix<F,U,V>& A )
{
#ifndef RELEASE
    CallStackEntry entry("MakeLegendre");
#endif
    if( A.Height() != A.Width() )
        LogicError("Cannot make a non-square matrix Legendre");
    MakeZeros( A );

    const Int localHeight = A.LocalHeight();
    const Int localWidth = A.LocalWidth();
    const Int colShift = A.ColShift();
    const Int rowShift = A.RowShift();
    const Int colStride = A.ColStride();
    const Int rowStride = A.RowStride();
    for( Int jLoc=0; jLoc<localWidth; ++jLoc )
    {
        const Int j = rowShift + jLoc*rowStride;
        for( Int iLoc=0; iLoc<localHeight; ++iLoc )
        {
            const Int i = colShift + iLoc*colStride;
            if( j == i+1 || j == i-1 )
            {
                const Int k = Max( i, j );
                const F gamma = F(1) / Pow( F(2)*k, F(2) );
                const F beta = F(1) / (2*Sqrt(F(1)-gamma));
                A.SetLocal( iLoc, jLoc, beta );
            }
        }
    }
}

inline typename Base<F>::type
HermitianEntrywiseOneNorm( UpperOrLower uplo, const DistMatrix<F>& A )
{
#ifndef RELEASE
    PushCallStack("HermitianEntrywiseOneNorm");
#endif
    if( A.Height() != A.Width() )
        throw std::logic_error("Hermitian matrices must be square.");

    const int r = A.Grid().Height();
    const int c = A.Grid().Width();
    const int colShift = A.ColShift();
    const int rowShift = A.RowShift();

    typedef typename Base<F>::type R;
    R localSum = 0;
    const int localWidth = A.LocalWidth();
    if( uplo == UPPER )
    {
        for( int jLocal=0; jLocal<localWidth; ++jLocal )
        {
            int j = rowShift + jLocal*c;
            int numUpperRows = Length(j+1,colShift,r);
            for( int iLocal=0; iLocal<numUpperRows; ++iLocal )
            {
                int i = colShift + iLocal*r;
                const R alpha = Abs(A.GetLocal(iLocal,jLocal));
                if( i ==j )
                    localSum += alpha;
                else
                    localSum += 2*alpha;
            }
        }
    }
    else
    {
        for( int jLocal=0; jLocal<localWidth; ++jLocal )
        {
            int j = rowShift + jLocal*c;
            int numStrictlyUpperRows = Length(j,colShift,r);
            for( int iLocal=numStrictlyUpperRows;
                 iLocal<A.LocalHeight(); ++iLocal )
            {
                int i = colShift + iLocal*r;
                const R alpha = Abs(A.GetLocal(iLocal,jLocal));
                if( i ==j )
                    localSum += alpha;
                else
                    localSum += 2*alpha;
            }
        }
    }

    R norm;
    mpi::AllReduce( &localSum, &norm, 1, mpi::SUM, A.Grid().VCComm() );
#ifndef RELEASE
    PopCallStack();
#endif
    return norm;
}

int Corrupt( DistMatrix<F>& A, double probCorrupt )
{
#ifndef RELEASE
    CallStackEntry entry("Corrupt");
#endif
    typedef BASE(F) Real;

    Int numLocalCorrupt = 0;
    const Int localHeight = A.LocalHeight();
    const Int localWidth = A.LocalWidth();
    for( Int jLocal=0; jLocal<localWidth; ++jLocal )
    {
        for( Int iLocal=0; iLocal<localHeight; ++iLocal )
        {
            if( Uniform<Real>() <= probCorrupt )
            {
                ++numLocalCorrupt;
                const F perturb = SampleBall<F>();
                A.SetLocal( iLocal, jLocal, A.GetLocal(iLocal,jLocal)+perturb );
            }
        }
    }
    
    Int numCorrupt;
    mpi::AllReduce
    ( &numLocalCorrupt, &numCorrupt, 1, mpi::SUM, A.Grid().VCComm() );
    return numCorrupt;
}

inline void
Riemann( DistMatrix<T,U,V>& R, int n )
{
#ifndef RELEASE
    CallStackEntry entry("Riemann");
#endif
    R.ResizeTo( n, n );
    const int localHeight = R.LocalHeight();
    const int localWidth = R.LocalWidth();
    const int colShift = R.ColShift();
    const int rowShift = R.RowShift();
    const int colStride = R.ColStride();
    const int rowStride = R.RowStride();
    for( int jLocal=0; jLocal<localWidth; ++jLocal )
    {
        const int j = rowShift + jLocal*rowStride;
        for( int iLocal=0; iLocal<localHeight; ++iLocal )
        {
            const int i = colShift + iLocal*colStride;
            if( ((j+2)%(i+2))==0 )
                R.SetLocal( iLocal, jLocal, T(i+1) );
            else
                R.SetLocal( iLocal, jLocal, T(-1) );
        }
    }
}

    static void Func
    ( DistMatrix<T,STAR,MR>& A, T center, typename Base<T>::type radius )
    {
        const Grid& grid = A.Grid();
        const int m = A.Height();
        const int localWidth = A.LocalWidth();
        const int bufSize = m*localWidth;
        std::vector<T> buffer( bufSize );

        // Create random matrix on process row 0, then broadcast
        if( grid.Row() == 0 )
        {
            for( int j=0; j<localWidth; ++j )
                for( int i=0; i<m; ++i )
                    buffer[i+j*m] = center+radius*SampleUnitBall<T>();
        }
        mpi::Broadcast( &buffer[0], bufSize, 0, grid.ColComm() );

        // Unpack
        T* localBuffer = A.LocalBuffer();
        const int ldim = A.LocalLDim();
#ifdef HAVE_OPENMP
        #pragma omp parallel for
#endif
        for( int jLocal=0; jLocal<localWidth; ++jLocal )
        {
            const T* bufferCol = &buffer[jLocal*m];
            T* col = &localBuffer[jLocal*ldim];
            MemCopy( col, bufferCol, m );
        }
    }

inline void
MakeDiscreteFourier( DistMatrix<Complex<R>,U,V>& A )
{
#ifndef RELEASE
    CallStackEntry entry("MakeDiscreteFourier");
#endif
    typedef Complex<R> F;

    const int m = A.Height();
    const int n = A.Width();
    if( m != n )
        throw std::logic_error("Cannot make a non-square DFT matrix");

    const R pi = 4*Atan( R(1) );
    const F nSqrt = Sqrt( R(n) );
    const int localHeight = A.LocalHeight();
    const int localWidth = A.LocalWidth();
    const int colShift = A.ColShift();
    const int rowShift = A.RowShift();
    const int colStride = A.ColStride();
    const int rowStride = A.RowStride();
    for( int jLocal=0; jLocal<localWidth; ++jLocal )
    {
        const int j = rowShift + jLocal*rowStride;
        for( int iLocal=0; iLocal<localHeight; ++iLocal )
        {
            const int i = colShift + iLocal*colStride;
            A.SetLocal( iLocal, jLocal, Exp(-2*pi*i*j/n)/nSqrt );

            const R theta = -2*pi*i*j/n;
            const Complex<R> alpha( Cos(theta), Sin(theta) );
            A.SetLocal( iLocal, jLocal, alpha/nSqrt );
        }
    }
}

inline void
BinaryFlat
( DistMatrix<T,STAR,V>& A, Int height, Int width, const std::string filename )
{
    DEBUG_ONLY(CallStackEntry cse("read::BinaryFlat"))
    std::ifstream file( filename.c_str(), std::ios::binary );
    if( !file.is_open() )
        RuntimeError("Could not open ",filename);

    const Int numBytes = FileSize( file );
    const Int numBytesExp = height*width*sizeof(T);
    if( numBytes != numBytesExp )
        RuntimeError
        ("Expected file to be ",numBytesExp," bytes but found ",numBytes);

    A.Resize( height, width );
    const Int localWidth = A.LocalWidth();
    const Int rowShift = A.RowShift();
    const Int rowStride = A.RowStride();
    for( Int jLoc=0; jLoc<localWidth; ++jLoc )
    {
        const Int j = rowShift + jLoc*rowStride;
        const Int localIndex = j*height;
        const std::streamoff pos = localIndex*sizeof(T);
        file.seekg( pos );
        file.read( (char*)A.Buffer(0,jLoc), height*sizeof(T) );
    }
}

inline void
Hankel( DistMatrix<T,U,V>& A, Int m, Int n, const std::vector<T>& a )
{
#ifndef RELEASE
    CallStackEntry entry("Hankel");
#endif
    const Int length = m+n-1;
    if( a.size() != (Unsigned)length )
        LogicError("a was the wrong size");
    A.ResizeTo( m, n );

    const Int localHeight = A.LocalHeight();
    const Int localWidth = A.LocalWidth();
    const Int colShift = A.ColShift();
    const Int rowShift = A.RowShift();
    const Int colStride = A.ColStride();
    const Int rowStride = A.RowStride();
    for( Int jLoc=0; jLoc<localWidth; ++jLoc )
    {
        const Int j = rowShift + jLoc*rowStride;
        for( Int iLoc=0; iLoc<localHeight; ++iLoc )
        {
            const Int i = colShift + iLoc*colStride;
            A.SetLocal( iLoc, jLoc, a[i+j] );
        }
    }
}

inline void
Wilkinson( DistMatrix<T,U,V>& A, int k )
{
#ifndef RELEASE
    CallStackEntry entry("Wilkinson");
#endif
    const int n = 2*k+1;
    A.ResizeTo( n, n );
    MakeZeros( A );

    const int localHeight = A.LocalHeight();
    const int localWidth = A.LocalWidth();
    const int colShift = A.ColShift();
    const int rowShift = A.RowShift();
    const int colStride = A.ColStride();
    const int rowStride = A.RowStride();
    for( int jLocal=0; jLocal<localWidth; ++jLocal )
    {
        const int j = rowShift + jLocal*rowStride;
        for( int iLocal=0; iLocal<localHeight; ++iLocal )
        {
            const int i = colShift + iLocal*colStride;
            if( i == j )
            {
                if( j <= k )
                    A.SetLocal( iLocal, jLocal, T(k-j) );
                else
                    A.SetLocal( iLocal, jLocal, T(j-k) );
            }
            else if( i == j-1 || i == j+1 )
                A.SetLocal( iLocal, jLocal, T(1) );
        }
    }
}

InfinityNorm( const DistMatrix<F,U,V>& A )
{
#ifndef RELEASE
    CallStackEntry entry("InfinityNorm");
#endif
    // Compute the partial row sums defined by our local matrix, A[U,V]
    typedef BASE(F) R;
    const Int localHeight = A.LocalHeight();
    const Int localWidth = A.LocalWidth();
    std::vector<R> myPartialRowSums( localHeight );
    for( Int iLoc=0; iLoc<localHeight; ++iLoc )
    {
        myPartialRowSums[iLoc] = 0;
        for( Int jLoc=0; jLoc<localWidth; ++jLoc )
            myPartialRowSums[iLoc] += Abs(A.GetLocal(iLoc,jLoc));
    }

    // Sum our partial row sums to get the row sums over A[U,* ]
    std::vector<R> myRowSums( localHeight );
    mpi::Comm rowComm = ReduceRowComm<U,V>( A.Grid() );
    mpi::AllReduce( &myPartialRowSums[0], &myRowSums[0], localHeight, rowComm );

    // Find the maximum out of the row sums
    R myMaxRowSum = 0;
    for( Int iLoc=0; iLoc<localHeight; ++iLoc )
        myMaxRowSum = std::max( myMaxRowSum, myRowSums[iLoc] );

    // Find the global maximum row sum by searching over the U team
    mpi::Comm colComm = ReduceColComm<U,V>( A.Grid() );
    return mpi::AllReduce( myMaxRowSum, mpi::MAX, colComm );
}