C++ DistMatrix::GetLocal方法代码示例

void dotproduct(DistMatrix<R> &A,DistMatrix<R> &B){
  if(A.Height() != B.Height()){
    //ERROR!
  }
  if(A.Width() != B.Width()){
    //ERROR!
  }
  double temp,temp1;
 
  const int colShift = A.ColShift(); // first row we own
  const int rowShift = A.RowShift(); // first col we own
  const int colStride = A.ColStride();
  const int rowStride = A.RowStride();
  const int localHeight = A.LocalHeight();
  const int localWidth = A.LocalWidth();
  for( int iLocal=0; iLocal<localHeight; ++iLocal ){
    for( int jLocal=0; jLocal<localWidth; ++jLocal ){
      const int i = colShift + iLocal*colStride;
      const int j = rowShift + jLocal*rowStride;
      temp = A.GetLocal(iLocal, jLocal);
      temp1 = B.GetLocal(iLocal,jLocal);
      B.SetLocal(iLocal, jLocal, (temp*temp1));
    }
  }
}

void calc_x_else(const R alpha, const R beta,const R c, const R delta, DistMatrix<R> &x, DistMatrix<R> &phi,DistMatrix<R> &w){

R bphi;
R temp;
R temp2;
const R pi = 4*atan(1);

  const int colShift = x.ColShift(); // first row we own
  const int rowShift = x.RowShift(); // first col we own
  const int colStride = x.ColStride();
  const int rowStride = x.RowStride();
  const int localHeight = x.LocalHeight();
  const int localWidth = x.LocalWidth();
  for( int iLocal=0; iLocal<localHeight; ++iLocal ){
    for( int jLocal=0; jLocal<localWidth; ++jLocal ){
      const int i = colShift + iLocal*colStride;
      const int j = rowShift + jLocal*rowStride;
      temp = phi.GetLocal(0,jLocal);//phi
      temp2 = w.GetLocal(0,jLocal); //w
      bphi = (pi/2) + beta * temp;
      //cout<< (delta + c * (((2/pi) * (bphi * tan(temp) - beta * log((pi/2) * temp2 * cos(temp) / bphi))) + (beta * tan (pi *alpha/2))))<<endl;
      x.SetLocal(iLocal, jLocal, (-1* sqrt(abs(delta + c * (((2/pi) * (bphi * tan(temp) - beta * log((pi/2) * temp2 * cos(temp) / bphi))) + (beta * tan (pi *alpha/2)))))));
    }
  }
} 

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

void calc_x_if(const R alpha,const R beta,const R c, const R delta,DistMatrix<R> &x, DistMatrix<R> &phi,DistMatrix<R> &w){
 
  R temp; // phi
  R temp2; // w
  R aphi;
  R a1phi;
  const R pi = 4*atan(1);
  const R zeta = beta * tan(pi * alpha/2);

  const int colShift = x.ColShift(); // first row we own
  const int rowShift = x.RowShift(); // first col we own
  const int colStride = x.ColStride();
  const int rowStride = x.RowStride();
  const int localHeight = x.LocalHeight();
  const int localWidth = x.LocalWidth();
  for( int iLocal=0; iLocal<localHeight; ++iLocal ){
    for( int jLocal=0; jLocal<localWidth; ++jLocal ){
      const int i = colShift + iLocal*colStride;
      const int j = rowShift + jLocal*rowStride;
      temp = phi.GetLocal(iLocal,jLocal);//phi
      temp2 = w.GetLocal(iLocal,jLocal); //w
       aphi = alpha * temp;
       a1phi = (1-alpha)*temp;

       x.SetLocal(iLocal, jLocal, (-1 * sqrt(abs(delta + c *( ( (sin(aphi)+zeta * cos(aphi))/ cos(temp)) * -1 * pow( abs( ((cos(a1phi) + zeta * sin(a1phi))/ (temp2 * cos(temp))) ) ,((1-alpha)/alpha) ) + beta * tan(pi * alpha/2) ))))  );


    }
  }
}

inline typename Base<F>::type
HermitianMaxNorm( UpperOrLower uplo, const DistMatrix<F>& A )
{
#ifndef RELEASE
    PushCallStack("internal::HermitianMaxNorm");
#endif
    typedef typename Base<F>::type R;

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

    R localMaxAbs = 0;
    const int localWidth = A.LocalWidth();
    if( uplo == UPPER )
    {
        for( int jLocal=0; jLocal<localWidth; ++jLocal )
        {
            int j = rowShift + jLocal*c;
            int numUpperRows = LocalLength(j+1,colShift,r);
            for( int iLocal=0; iLocal<numUpperRows; ++iLocal )
            {
                const R thisAbs = Abs(A.GetLocal(iLocal,jLocal));
                localMaxAbs = std::max( localMaxAbs, thisAbs );
            }
        }
    }
    else
    {
        for( int jLocal=0; jLocal<localWidth; ++jLocal )
        {
            int j = rowShift + jLocal*c;
            int numStrictlyUpperRows = LocalLength(j,colShift,r);
            for( int iLocal=numStrictlyUpperRows; 
                 iLocal<A.LocalHeight(); ++iLocal )
            {
                const R thisAbs = Abs(A.GetLocal(iLocal,jLocal));
                localMaxAbs = std::max( localMaxAbs, thisAbs );
            }
        }
    }

    R maxAbs;
    mpi::AllReduce( &localMaxAbs, &maxAbs, 1, mpi::MAX, A.Grid().VCComm() );
#ifndef RELEASE
    PopCallStack();
#endif
    return maxAbs;
}

void GetMappedDiagonal
( const DistMatrix<T,U,V,BLOCK>& A,
        AbstractDistMatrix<S>& d,
        function<S(const T&)> func,
        Int offset )
{
    EL_DEBUG_CSE
    EL_DEBUG_ONLY(AssertSameGrids( A, d ))

    // TODO(poulson): Make this more efficient
    const Int diagLength = A.DiagonalLength(offset);
    d.Resize( diagLength, 1 );
    Zero( d );
    if( d.Participating() && A.RedundantRank() == 0 )
    {
        const Int iStart = Max(-offset,0);
        const Int jStart = Max( offset,0);
        for( Int k=0; k<diagLength; ++k )
        {
            if( A.IsLocal(iStart+k,jStart+k) )
            {
                const Int iLoc = A.LocalRow(iStart+k);
                const Int jLoc = A.LocalCol(jStart+k);
                d.QueueUpdate(k,0,func(A.GetLocal(iLoc,jLoc)));
            }
        }
    }
    d.ProcessQueues();
}

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

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 bool
PivotParity( const DistMatrix<int,VC,STAR>& p, int pivotOffset ) 
{
#ifndef RELEASE
    PushCallStack("PivotParity");
    if( p.Width() != 1 )
        throw std::logic_error("p must be a column vector");
    if( pivotOffset < 0 )
        throw std::logic_error("pivot offset cannot be negative");
#endif
    const int mLocal = p.LocalHeight();
    const Grid& g = p.Grid();

    bool isLocallyOdd = false;
    const int colShift = p.ColShift();
    const int gridSize = g.Size();
    for( int iLocal=0; iLocal<mLocal; ++iLocal )
    {
        const int i = colShift + iLocal*gridSize;
        if( p.GetLocal(iLocal,0) != i+pivotOffset )
            isLocallyOdd = !isLocallyOdd;
    }

    int localContribution = isLocallyOdd;
    int globalContribution;
    mpi::AllReduce
    ( &localContribution, &globalContribution, 1, MPI_SUM, g.VCComm() );
    globalContribution = globalContribution % 2;

    bool isOdd = globalContribution;
#ifndef RELEASE
    PopCallStack();
#endif
    return isOdd;
}

inline typename Base<F>::type
FrobeniusNorm( const DistMatrix<F,U,V>& A )
{
#ifndef RELEASE
    PushCallStack("internal::FrobeniusNorm");
#endif
    typedef typename Base<F>::type R;
    mpi::Comm comm = NormComm( A );

    R localScale = 0;
    R localScaledSquare = 1;
    const int localHeight = A.LocalHeight();
    const int localWidth = A.LocalWidth();
    for( int jLocal=0; jLocal<localWidth; ++jLocal )
    {
        for( int iLocal=0; iLocal<localHeight; ++iLocal )
        {
            const R alphaAbs = Abs(A.GetLocal(iLocal,jLocal));
            if( alphaAbs != 0 )
            {
                if( alphaAbs <= localScale )
                {
                    const R relScale = alphaAbs/localScale;
                    localScaledSquare += relScale*relScale;
                }
                else
                {
                    const R relScale = localScale/alphaAbs;
                    localScaledSquare = localScaledSquare*relScale*relScale + 1;
                    localScale = alphaAbs; 
                }
            }
        }
    }

    // Find the maximum relative scale
    R scale;
    mpi::AllReduce( &localScale, &scale, 1, mpi::MAX, comm );

    R norm = 0;
    if( scale != 0 )
    {
        // Equilibrate our local scaled sum to the maximum scale
        R relScale = localScale/scale;
        localScaledSquare *= relScale*relScale;

        // The scaled square is now simply the sum of the local contributions
        R scaledSquare;
        mpi::AllReduce( &localScaledSquare, &scaledSquare, 1, mpi::SUM, comm );

        norm = scale*Sqrt(scaledSquare);
    }
#ifndef RELEASE
    PopCallStack();
#endif
    return norm;
}

void NormalizeEntries( DistMatrix<F,U,V>& A )
{
    const Int localHeight = A.LocalHeight();
    const Int localWidth = A.LocalWidth();
    for( Int jLocal=0; jLocal<localWidth; ++jLocal )
    {
        for( Int iLocal=0; iLocal<localHeight; ++iLocal )
        {
            const F alpha = A.GetLocal( iLocal, jLocal );
            A.SetLocal( iLocal, jLocal, alpha/Abs(alpha) );
        }
    }
}

inline void HermitianSVD
( UpperOrLower uplo, DistMatrix<F>& A, 
  DistMatrix<BASE(F),VR,STAR>& s, DistMatrix<F>& U, DistMatrix<F>& V )
{
#ifndef RELEASE
    CallStackEntry entry("HermitianSVD");
#endif
#ifdef HAVE_PMRRR
    typedef BASE(F) R;

    // Grab an eigenvalue decomposition of A
    HermitianEig( uplo, A, s, V );

    // Redistribute the singular values into an [MR,* ] distribution
    const Grid& grid = A.Grid();
    DistMatrix<R,MR,STAR> s_MR_STAR( grid );
    s_MR_STAR.AlignWith( V.DistData() );
    s_MR_STAR = s;

    // Set the singular values to the absolute value of the eigenvalues
    const Int numLocalVals = s.LocalHeight();
    for( Int iLoc=0; iLoc<numLocalVals; ++iLoc )
    {
        const R sigma = s.GetLocal(iLoc,0);
        s.SetLocal(iLoc,0,Abs(sigma));
    }

    // Copy V into U (flipping the sign as necessary)
    U.AlignWith( V );
    U.ResizeTo( V.Height(), V.Width() );
    const Int localHeight = V.LocalHeight();
    const Int localWidth = V.LocalWidth();
    for( Int jLoc=0; jLoc<localWidth; ++jLoc )
    {
        const R sigma = s_MR_STAR.GetLocal( jLoc, 0 );
        F* UCol = U.Buffer( 0, jLoc );
        const F* VCol = V.LockedBuffer( 0, jLoc );
        if( sigma >= 0 )
            for( Int iLoc=0; iLoc<localHeight; ++iLoc )
                UCol[iLoc] = VCol[iLoc];
        else
            for( Int iLoc=0; iLoc<localHeight; ++iLoc )
                UCol[iLoc] = -VCol[iLoc];
    }
#else
    U = A;
    MakeHermitian( uplo, U );
    SVD( U, s, V );
#endif // ifdef HAVE_PMRRR
}

FrobeniusNorm( const DistMatrix<F,U,V>& A )
{
#ifndef RELEASE
    CallStackEntry entry("FrobeniusNorm");
#endif
    typedef BASE(F) R;
    R localScale = 0;
    R localScaledSquare = 1;
    const Int localHeight = A.LocalHeight();
    const Int localWidth = A.LocalWidth();
    for( Int jLoc=0; jLoc<localWidth; ++jLoc )
    {
        for( Int iLoc=0; iLoc<localHeight; ++iLoc )
        {
            const R alphaAbs = Abs(A.GetLocal(iLoc,jLoc));
            if( alphaAbs != 0 )
            {
                if( alphaAbs <= localScale )
                {
                    const R relScale = alphaAbs/localScale;
                    localScaledSquare += relScale*relScale;
                }
                else
                {
                    const R relScale = localScale/alphaAbs;
                    localScaledSquare = localScaledSquare*relScale*relScale + 1;
                    localScale = alphaAbs; 
                }
            }
        }
    }

    // Find the maximum relative scale
    mpi::Comm comm = ReduceComm<U,V>( A.Grid() );
    const R scale = mpi::AllReduce( localScale, mpi::MAX, comm );

    R norm = 0;
    if( scale != 0 )
    {
        // Equilibrate our local scaled sum to the maximum scale
        R relScale = localScale/scale;
        localScaledSquare *= relScale*relScale;

        // The scaled square is now simply the sum of the local contributions
        const R scaledSquare = mpi::AllReduce( localScaledSquare, comm );
        norm = scale*Sqrt(scaledSquare);
    }
    return norm;
}

inline void HermitianSVD
( UpperOrLower uplo, 
  DistMatrix<F>& A, DistMatrix<typename Base<F>::type,VR,STAR>& s, 
  DistMatrix<F>& U, DistMatrix<F>& V )
{
#ifndef RELEASE
    PushCallStack("HermitianSVD");
#endif
    typedef typename Base<F>::type R;

    // Grab an eigenvalue decomposition of A
    HermitianEig( uplo, A, s, V ); 
    
    // Redistribute the singular values into an [MR,* ] distribution
    const Grid& grid = A.Grid();
    DistMatrix<R,MR,STAR> s_MR_STAR( grid );
    s_MR_STAR.AlignWith( V );
    s_MR_STAR = s;

    // Set the singular values to the absolute value of the eigenvalues
    const int numLocalVals = s.LocalHeight();
    for( int iLocal=0; iLocal<numLocalVals; ++iLocal )
    {
        const R sigma = s.GetLocal(iLocal,0);
        s.SetLocal(iLocal,0,Abs(sigma));
    }

    // Copy V into U (flipping the sign as necessary)
    U.AlignWith( V );
    U.ResizeTo( V.Height(), V.Width() );
    const int localHeight = V.LocalHeight();
    const int localWidth = V.LocalWidth();
    for( int jLocal=0; jLocal<localWidth; ++jLocal )
    {
        const R sigma = s_MR_STAR.GetLocal( jLocal, 0 );
        F* UCol = U.LocalBuffer( 0, jLocal );
        const F* VCol = V.LockedLocalBuffer( 0, jLocal );
        if( sigma >= 0 )
            for( int iLocal=0; iLocal<localHeight; ++iLocal )
                UCol[iLocal] = VCol[iLocal];
        else
            for( int iLocal=0; iLocal<localHeight; ++iLocal )
                UCol[iLocal] = -VCol[iLocal];
    }
#ifndef RELEASE
    PopCallStack();
#endif
}

void repmat_as2(DistMatrix<R> &vector, DistMatrix<R> &repmatrix){
  //old 1 x d
  //new nSim x d
  //Must set dimensions of New before being passed to function
  double temp;
  const int colShift = repmatrix.ColShift(); // first row we own
  const int rowShift = repmatrix.RowShift(); // first col we own
  const int colStride = repmatrix.ColStride();
  const int rowStride = repmatrix.RowStride();
  const int localHeight = repmatrix.LocalHeight();
  const int localWidth = repmatrix.LocalWidth();
  for( int iLocal=0; iLocal<localHeight; ++iLocal ){
    for( int jLocal=0; jLocal<localWidth; ++jLocal ){
      const int i = colShift + iLocal*colStride;
      const int j = rowShift + jLocal*rowStride;
      temp = vector.GetLocal(jLocal, 0);
      repmatrix.SetLocal(iLocal, jLocal, temp);
    }
  }
}