自学教程：C++ AbstractDistMatrix类代码示例

void LUNMedium( const AbstractDistMatrix<F>& UPre, AbstractDistMatrix<F>& XPre,   bool checkIfSingular ){    DEBUG_CSE    const Int m = XPre.Height();    const Int bsize = Blocksize();    const Grid& g = UPre.Grid();    DistMatrixReadProxy<F,F,MC,MR> UProx( UPre );    DistMatrixReadWriteProxy<F,F,MC,MR> XProx( XPre );    auto& U = UProx.GetLocked();    auto& X = XProx.Get();    DistMatrix<F,MC,  STAR> U01_MC_STAR(g);    DistMatrix<F,STAR,STAR> U11_STAR_STAR(g);    DistMatrix<F,MR,  STAR> X1Trans_MR_STAR(g);    const Int kLast = LastOffset( m, bsize );    Int k=kLast, kOld=m;    while( true )    {        const bool in2x2 = ( k>0 && U.Get(k,k-1) != F(0) );        if( in2x2 )            --k;        const Int nb = kOld-k;        const Range<Int> ind0( 0, k    ),                         ind1( k, k+nb );        auto U01 = U( ind0, ind1 );        auto U11 = U( ind1, ind1 );        auto X0 = X( ind0, ALL );        auto X1 = X( ind1, ALL );        U11_STAR_STAR = U11; // U11[* ,* ] <- U11[MC,MR]        X1Trans_MR_STAR.AlignWith( X0 );        Transpose( X1, X1Trans_MR_STAR );                // X1^T[MR,* ] := X1^T[MR,* ] U11^-T[* ,* ]        //              = (U11^-1[* ,* ] X1[* ,MR])^T        LocalQuasiTrsm        ( RIGHT, UPPER, TRANSPOSE,          F(1), U11_STAR_STAR, X1Trans_MR_STAR, checkIfSingular );        Transpose( X1Trans_MR_STAR, X1 );        U01_MC_STAR.AlignWith( X0 );        U01_MC_STAR = U01;  // U01[MC,* ] <- U01[MC,MR]        // X0[MC,MR] -= U01[MC,* ] X1[* ,MR]        LocalGemm        ( NORMAL, TRANSPOSE, F(-1), U01_MC_STAR, X1Trans_MR_STAR, F(1), X0 );        if( k == 0 )            break;        kOld = k;        k -= Min(bsize,k);    }}

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

void UN_C( T alpha,  const AbstractDistMatrix<T>& APre,         AbstractDistMatrix<T>& CPre,  bool conjugate=false ){    EL_DEBUG_CSE    const Int r = APre.Width();    const Int bsize = Blocksize();    const Grid& g = APre.Grid();    DistMatrixReadProxy<T,T,MC,MR> AProx( APre );    DistMatrixReadWriteProxy<T,T,MC,MR> CProx( CPre );    auto& A = AProx.GetLocked();    auto& C = CProx.Get();    // Temporary distributions    DistMatrix<T,MC,  STAR> A1_MC_STAR(g);    DistMatrix<T,VR,  STAR> A1_VR_STAR(g);    DistMatrix<T,STAR,MR  > A1Trans_STAR_MR(g);    A1_MC_STAR.AlignWith( C );    A1_VR_STAR.AlignWith( C );    A1Trans_STAR_MR.AlignWith( C );    for( Int k=0; k<r; k+=bsize )    {        const Int nb = Min(bsize,r-k);        auto A1 = A( ALL, IR(k,k+nb) );        A1_VR_STAR = A1_MC_STAR = A1;        Transpose( A1_VR_STAR, A1Trans_STAR_MR, conjugate );        LocalTrrk( UPPER, alpha, A1_MC_STAR, A1Trans_STAR_MR, T(1), C );     }}

void MinEig( const AbstractDistMatrix<Real>& xPre,  AbstractDistMatrix<Real>& minEigsPre,  const AbstractDistMatrix<Int>& orders,  const AbstractDistMatrix<Int>& firstIndsPre,  Int cutoff ){    EL_DEBUG_CSE    AssertSameGrids( xPre, minEigsPre, orders, firstIndsPre );    ElementalProxyCtrl ctrl;    ctrl.colConstrain = true;    ctrl.colAlign = 0;    DistMatrixReadProxy<Real,Real,VC,STAR>    xProx( xPre, ctrl );    DistMatrixWriteProxy<Real,Real,VC,STAR>    minEigsProx( minEigsPre, ctrl );    DistMatrixReadProxy<Int,Int,VC,STAR>    firstIndsProx( firstIndsPre, ctrl );    auto& x = xProx.GetLocked();    auto& minEigs = minEigsProx.Get();    auto& firstInds = firstIndsProx.GetLocked();    const Int height = x.Height();    const Int localHeight = x.LocalHeight();    EL_DEBUG_ONLY(        if( x.Width() != 1 || orders.Width() != 1 || firstInds.Width() != 1 )        LogicError("x, orders, and firstInds should be column vectors");        if( orders.Height() != height || firstInds.Height() != height )            LogicError("orders and firstInds should be of the same height as x");        )

void Mehrotra( const AbstractDistMatrix<Real>& A,  const AbstractDistMatrix<Real>& b,  const AbstractDistMatrix<Real>& c,  const AbstractDistMatrix<Int>& orders,  const AbstractDistMatrix<Int>& firstInds,        AbstractDistMatrix<Real>& x,        AbstractDistMatrix<Real>& y,        AbstractDistMatrix<Real>& z,  const MehrotraCtrl<Real>& ctrl ){    EL_DEBUG_CSE    const Int n = c.Height();    const Grid& grid = c.Grid();    DistMatrix<Real> G(grid);    Identity( G, n, n );    G *= -1;    DistMatrix<Real> h(grid);    Zeros( h, n, 1 );    MehrotraCtrl<Real> affineCtrl = ctrl;    affineCtrl.primalInit = false;    affineCtrl.dualInit = false;    DistMatrix<Real> s(grid);    socp::affine::Mehrotra(A,G,b,c,h,orders,firstInds,x,y,z,s,affineCtrl);}

void ExplicitTriang( AbstractDistMatrix<F>& A ){    DEBUG_ONLY(CallStackEntry cse("rq::ExplicitTriang"))    DistMatrix<F,MD,STAR> t(A.Grid());    DistMatrix<Base<F>,MD,STAR> d(A.Grid());    Householder( A, t, d );    MakeTrapezoidal( UPPER, A, A.Width()-A.Height() );}

inline void AssertSameDists( const AbstractDistMatrix<T>& A1, const AbstractDistMatrix<T>& A2,  Args&... args ) {    if( A1.ColDist() != A2.ColDist() || A1.RowDist() != A2.RowDist() )        LogicError("Distributions did not match");    AssertSameDists( A2, args... );}

inline voidAbstractDistMatrix<T,Int>::AssertSameSizeAsTranspose( const AbstractDistMatrix<U,Int>& A ) const{    if( Height() != A.Width() || Width() != A.Height() )        throw std::logic_error        ("Assertion that matrices be the same size (after trans.) failed");}

void ApplyQ( LeftOrRight side,  Orientation orientation,  const AbstractDistMatrix<F>& APre,  const AbstractDistMatrix<F>& householderScalars,  const AbstractDistMatrix<Base<F>>& signature,        AbstractDistMatrix<F>& BPre ){    EL_DEBUG_CSE    const bool normal = (orientation==NORMAL);    const bool onLeft = (side==LEFT);    const bool applyDFirst = normal==onLeft;    const Int minDim = Min(APre.Height(),APre.Width());    const ForwardOrBackward direction = ( normal==onLeft ? BACKWARD : FORWARD );    const Conjugation conjugation = ( normal ? CONJUGATED : UNCONJUGATED );    DistMatrixReadProxy<F,F,MC,MR> AProx( APre );    DistMatrixReadWriteProxy<F,F,MC,MR> BProx( BPre );    auto& A = AProx.GetLocked();    auto& B = BProx.Get();    const Int m = B.Height();    const Int n = B.Width();    if( applyDFirst )    {        if( onLeft )        {            auto BTop = B( IR(0,minDim), IR(0,n) );            DiagonalScale( side, orientation, signature, BTop );        }        else        {            auto BLeft = B( IR(0,m), IR(0,minDim) );            DiagonalScale( side, orientation, signature, BLeft );        }    }    ApplyPackedReflectors    ( side, LOWER, VERTICAL, direction, conjugation, 0,      A, householderScalars, B );    if( !applyDFirst )    {        if( onLeft )        {            auto BTop = B( IR(0,minDim), IR(0,n) );            DiagonalScale( side, orientation, signature, BTop );        }        else        {            auto BLeft = B( IR(0,m), IR(0,minDim) );            DiagonalScale( side, orientation, signature, BLeft );        }    }}

void LAV( const AbstractDistMatrix<Real>& A,  const AbstractDistMatrix<Real>& b,        AbstractDistMatrix<Real>& xPre,  const lp::affine::Ctrl<Real>& ctrl ){    EL_DEBUG_CSE    DistMatrixWriteProxy<Real,Real,MC,MR> xProx( xPre );    auto& x = xProx.Get();    const Int m = A.Height();    const Int n = A.Width();    const Grid& g = A.Grid();    const Range<Int> xInd(0,n), uInd(n,n+m), vInd(n+m,n+2*m);    DistMatrix<Real> c(g), AHat(g), G(g), h(g);    // c := [0;1;1]    // ============    Zeros( c, n+2*m, 1 );    auto cuv = c( IR(n,n+2*m), ALL );    Fill( cuv, Real(1) );    // /hat A := [A, I, -I]    // ====================    Zeros( AHat, m, n+2*m );    auto AHatx = AHat( IR(0,m), xInd );    auto AHatu = AHat( IR(0,m), uInd );    auto AHatv = AHat( IR(0,m), vInd );    AHatx = A;    FillDiagonal( AHatu, Real( 1) );    FillDiagonal( AHatv, Real(-1) );    // G := | 0 -I  0 |    //      | 0  0 -I |    // ================    Zeros( G, 2*m, n+2*m );    auto Guv = G( IR(0,2*m), IR(n,n+2*m) );    FillDiagonal( Guv, Real(-1) );    // h := | 0 |    //      | 0 |    // ==========    Zeros( h, 2*m, 1 );    // Solve the affine linear program    // ===============================    DistMatrix<Real> xHat(g), y(g), z(g), s(g);    LP( AHat, G, b, c, h, xHat, y, z, s, ctrl );    // Extract x    // =========    x = xHat( xInd, ALL );}

void Her2k( UpperOrLower uplo, Orientation orientation,  T alpha, const AbstractDistMatrix<T>& A, const AbstractDistMatrix<T>& B,                 AbstractDistMatrix<T>& C ){    EL_DEBUG_CSE    const Int n = ( orientation==NORMAL ? A.Height() : A.Width() );    C.Resize( n, n );    Zero( C );    Syr2k( uplo, orientation, alpha, A, B, T(0), C, true );}

void LLNMedium( const AbstractDistMatrix<F>& LPre,        AbstractDistMatrix<F>& XPre,   bool checkIfSingular ){    DEBUG_CSE    const Int m = XPre.Height();    const Int bsize = Blocksize();    const Grid& g = LPre.Grid();    DistMatrixReadProxy<F,F,MC,MR> LProx( LPre );    DistMatrixReadWriteProxy<F,F,MC,MR> XProx( XPre );    auto& L = LProx.GetLocked();    auto& X = XProx.Get();    DistMatrix<F,STAR,STAR> L11_STAR_STAR(g);    DistMatrix<F,MC,  STAR> L21_MC_STAR(g);    DistMatrix<F,MR,  STAR> X1Trans_MR_STAR(g);    for( Int k=0; k<m; k+=bsize )    {        const Int nbProp = Min(bsize,m-k);        const bool in2x2 = ( k+nbProp<m && L.Get(k+nbProp-1,k+nbProp) != F(0) );        const Int nb = ( in2x2 ? nbProp+1 : nbProp );        const Range<Int> ind1( k,    k+nb ),                         ind2( k+nb, m    );        auto L11 = L( ind1, ind1 );        auto L21 = L( ind2, ind1 );        auto X1 = X( ind1, ALL );        auto X2 = X( ind2, ALL );        L11_STAR_STAR = L11; // L11[* ,* ] <- L11[MC,MR]        X1Trans_MR_STAR.AlignWith( X2 );        Transpose( X1, X1Trans_MR_STAR );        // X1^T[MR,* ] := X1^T[MR,* ] L11^-T[* ,* ]        //              = (L11^-1[* ,* ] X1[* ,MR])^T        LocalQuasiTrsm        ( RIGHT, LOWER, TRANSPOSE,          F(1), L11_STAR_STAR, X1Trans_MR_STAR, checkIfSingular );        Transpose( X1Trans_MR_STAR, X1 );        L21_MC_STAR.AlignWith( X2 );        L21_MC_STAR = L21;                   // L21[MC,* ] <- L21[MC,MR]                // X2[MC,MR] -= L21[MC,* ] X1[* ,MR]        LocalGemm        ( NORMAL, TRANSPOSE, F(-1), L21_MC_STAR, X1Trans_MR_STAR, F(1), X2 );    }}

void LUNMedium( UnitOrNonUnit diag,   const AbstractDistMatrix<F>& UPre,        AbstractDistMatrix<F>& XPre,  bool checkIfSingular ){    EL_DEBUG_CSE    const Int m = XPre.Height();    const Int bsize = Blocksize();    const Grid& g = UPre.Grid();    DistMatrixReadProxy<F,F,MC,MR> UProx( UPre );    DistMatrixReadWriteProxy<F,F,MC,MR> XProx( XPre );    auto& U = UProx.GetLocked();    auto& X = XProx.Get();    DistMatrix<F,MC,  STAR> U01_MC_STAR(g);    DistMatrix<F,STAR,STAR> U11_STAR_STAR(g);    DistMatrix<F,MR,  STAR> X1Trans_MR_STAR(g);    const Int kLast = LastOffset( m, bsize );    for( Int k=kLast; k>=0; k-=bsize )    {        const Int nb = Min(bsize,m-k);        const Range<Int> ind0( 0, k    ),                         ind1( k, k+nb );        auto U01 = U( ind0, ind1 );        auto U11 = U( ind1, ind1 );        auto X0 = X( ind0, ALL );        auto X1 = X( ind1, ALL );        U11_STAR_STAR = U11; // U11[* ,* ] <- U11[MC,MR]        X1Trans_MR_STAR.AlignWith( X0 );        Transpose( X1, X1Trans_MR_STAR );                // X1^T[MR,* ] := X1^T[MR,* ] U11^-T[* ,* ]        //              = (U11^-1[* ,* ] X1[* ,MR])^T        LocalTrsm        ( RIGHT, UPPER, TRANSPOSE, diag,           F(1), U11_STAR_STAR, X1Trans_MR_STAR, checkIfSingular );        Transpose( X1Trans_MR_STAR, X1 );        U01_MC_STAR.AlignWith( X0 );        U01_MC_STAR = U01;  // U01[MC,* ] <- U01[MC,MR]        // X0[MC,MR] -= U01[MC,* ] X1[* ,MR]        LocalGemm        ( NORMAL, TRANSPOSE, F(-1), U01_MC_STAR, X1Trans_MR_STAR, F(1), X0 );    }}

Int ZeroNorm( const AbstractDistMatrix<T>& A, Base<T> tol ){    DEBUG_ONLY(CSE cse("ZeroNorm"))    Int numNonzeros;    if( A.Participating() )    {        const Int numLocalNonzeros = ZeroNorm( A.LockedMatrix(), tol );        numNonzeros = mpi::AllReduce( numLocalNonzeros, A.DistComm() );    }    mpi::Broadcast( numNonzeros, A.Root(), A.CrossComm() );    return numNonzeros;}

void Lotkin( AbstractDistMatrix<F>& A, Int n ){    DEBUG_ONLY(CallStackEntry cse("Lotkin"))    Hilbert( A, n );    // Set first row to all ones    if( A.ColShift() == 0 )    {        const Int localWidth = A.LocalWidth();        for( Int jLoc=0; jLoc<localWidth; ++jLoc )            A.SetLocal( 0, jLoc, F(1) );    } }

void Lauchli( AbstractDistMatrix<T>& A, Int n, T mu ){    DEBUG_ONLY(CallStackEntry cse("Lauchli"))    Zeros( A, n+1, n );    // Set the first row to all ones    unique_ptr<AbstractDistMatrix<T>> a0( A.Construct(A.Grid(),A.Root()) );    View( *a0, A, IR(0,1), IR(0,n) );    Fill( *a0, T(1) );    // Set the subdiagonal to mu    FillDiagonal( A, mu, -1 );}

void ExplicitTriang( AbstractDistMatrix<F>& A ){    EL_DEBUG_CSE    const Grid& g = A.Grid();    DistMatrix<F,MD,STAR> householderScalars(g);    DistMatrix<Base<F>,MD,STAR> signature(g);    LQ( A, householderScalars, signature );    const Int m = A.Height();    const Int n = A.Width();    const Int minDim = Min(m,n);    A.Resize( m, minDim );    MakeTrapezoidal( LOWER, A );}

Base<F> Coherence( const AbstractDistMatrix<F>& A ){    DEBUG_ONLY(CallStackEntry cse("Coherence"))    DistMatrix<F> B( A );    DistMatrix<Base<F>,MR,STAR> norms(B.Grid());    ColumnNorms( B, norms );    DiagonalSolve( RIGHT, NORMAL, norms, B, true );    DistMatrix<F> C(B.Grid());    Identity( C, A.Width(), A.Width() );    Herk( UPPER, ADJOINT, Base<F>(-1), B, Base<F>(1), C );    return HermitianMaxNorm( UPPER, C );}

void RestoreOrdering( const AbstractDistMatrix<Int>& preimage,        AbstractDistMatrix<T>& x ){    EL_DEBUG_CSE    DistMatrix<Int,STAR,STAR> preimageCopy( preimage );    DistMatrix<T,STAR,STAR> xCopy( x );    const Int numShifts = preimage.Height();    // TODO(poulson): Significantly lower the latency    for( Int j=0; j<numShifts; ++j )    {        const Int dest = preimageCopy.Get(j,0);        x.Set( dest, 0, xCopy.Get(j,0) );    }}

void Explicit( AbstractDistMatrix<F>& L, AbstractDistMatrix<F>& APre ){    EL_DEBUG_CSE    const Grid& g = APre.Grid();    DistMatrixReadWriteProxy<F,F,MC,MR> AProx( APre );    auto& A = AProx.Get();    DistMatrix<F,MD,STAR> householderScalars(g);    DistMatrix<Base<F>,MD,STAR> signature(g);    LQ( A, householderScalars, signature );    const Int m = A.Height();    const Int n = A.Width();    const Int minDim = Min(m,n);    auto AL = A( IR(0,m), IR(0,minDim) );    Copy( AL, L );    MakeTrapezoidal( LOWER, L );    // TODO: Replace this with an in-place expansion of Q    DistMatrix<F> Q(g);    Identity( Q, A.Height(), A.Width() );    lq::ApplyQ( RIGHT, NORMAL, A, householderScalars, signature, Q );    Copy( Q, APre );}

void Hilbert( AbstractDistMatrix<F>& A, Int n ){    DEBUG_ONLY(CSE cse("Hilbert"))    A.Resize( n, n );    auto hilbertFill = []( Int i, Int j ) { return F(1)/F(i+j+1); };    IndexDependentFill( A, function<F(Int,Int)>(hilbertFill) );}

void Forsythe( AbstractDistMatrix<T>& J, Int n, T alpha, T lambda ){    DEBUG_ONLY(CSE cse("Forsythe"))    Jordan( J, n, lambda );    if( n > 0 )        J.Set( n-1, 0, alpha );}

void MinIJ( AbstractDistMatrix<T>& M, Int n ){    DEBUG_ONLY(CSE cse("MinIJ"))    M.Resize( n, n );    auto minIJFill = []( Int i, Int j ) { return T(Min(i+1,j+1)); };    IndexDependentFill( M, function<T(Int,Int)>(minIJFill) );}