C++ DistSparseMatrix类代码示例

void SymmetricRuizEquil
( DistSparseMatrix<F>& A, 
  DistMultiVec<Base<F>>& d, 
  Int maxIter, bool progress )
{
    DEBUG_CSE
    typedef Base<F> Real;
    const Int n = A.Height();
    mpi::Comm comm = A.Comm();
    d.SetComm( comm );
    Ones( d, n, 1 );

    DistMultiVec<Real> scales(comm);
    const Int indent = PushIndent();
    for( Int iter=0; iter<maxIter; ++iter )
    {
        // Rescale the columns (and rows)
        // ------------------------------
        ColumnMaxNorms( A, scales );
        EntrywiseMap( scales, function<Real(Real)>(DampScaling<Real>) );
        EntrywiseMap( scales, function<Real(Real)>(SquareRootScaling<Real>) );
        DiagonalScale( LEFT, NORMAL, scales, d );
        SymmetricDiagonalSolve( scales, A );
    }
    SetIndent( indent );
}

void SymmetricRuizEquil
( DistSparseMatrix<Field>& A,
  DistMultiVec<Base<Field>>& d,
  Int maxIter, bool progress )
{
    EL_DEBUG_CSE
    typedef Base<Field> Real;
    const Int n = A.Height();
    const Grid& grid = A.Grid();
    d.SetGrid( grid );
    Ones( d, n, 1 );

    DistMultiVec<Real> scales(grid);
    const Int indent = PushIndent();
    for( Int iter=0; iter<maxIter; ++iter )
    {
        // Rescale the columns (and rows)
        // ------------------------------
        ColumnMaxNorms( A, scales );
        EntrywiseMap( scales, MakeFunction(DampScaling<Real>) );
        EntrywiseMap( scales, MakeFunction(SquareRootScaling<Real>) );
        DiagonalScale( LEFT, NORMAL, scales, d );
        SymmetricDiagonalSolve( scales, A );
    }
    SetIndent( indent );
}

void EntrywiseMap( DistSparseMatrix<T>& A, function<T(T)> func )
{
    DEBUG_CSE
    T* vBuf = A.ValueBuffer();
    const Int numLocalEntries = A.NumLocalEntries();
    for( Int k=0; k<numLocalEntries; ++k )
        vBuf[k] = func(vBuf[k]);
}

void StackedRuizEquil
( DistSparseMatrix<Field>& A,
  DistSparseMatrix<Field>& B,
  DistMultiVec<Base<Field>>& dRowA,
  DistMultiVec<Base<Field>>& dRowB,
  DistMultiVec<Base<Field>>& dCol,
  bool progress )
{
    EL_DEBUG_CSE
    typedef Base<Field> Real;
    const Int mA = A.Height();
    const Int mB = B.Height();
    const Int n = A.Width();
    mpi::Comm comm = A.Comm();
    dRowA.SetComm( comm );
    dRowB.SetComm( comm );
    dCol.SetComm( comm );
    Ones( dRowA, mA, 1 );
    Ones( dRowB, mB, 1 );
    Ones( dCol, n, 1 );

    // TODO(poulson): Expose to control structure
    // For, simply hard-code a small number of iterations
    const Int maxIter = 4;

    DistMultiVec<Real> scales(comm), maxAbsValsB(comm);
    auto& scalesLoc = scales.Matrix();
    auto& maxAbsValsBLoc = maxAbsValsB.Matrix();
    const Int localHeight = scalesLoc.Height();
    const Int indent = PushIndent();
    for( Int iter=0; iter<maxIter; ++iter )
    {
        // Rescale the columns
        // -------------------
        ColumnMaxNorms( A, scales );
        ColumnMaxNorms( B, maxAbsValsB );
        for( Int jLoc=0; jLoc<localHeight; ++jLoc )
            scalesLoc(jLoc) = Max(scalesLoc(jLoc),maxAbsValsBLoc(jLoc));
        EntrywiseMap( scales, MakeFunction(DampScaling<Real>) );
        DiagonalScale( LEFT, NORMAL, scales, dCol );
        DiagonalSolve( RIGHT, NORMAL, scales, A );
        DiagonalSolve( RIGHT, NORMAL, scales, B );

        // Rescale the rows
        // ----------------
        RowMaxNorms( A, scales );
        EntrywiseMap( scales, MakeFunction(DampScaling<Real>) );
        DiagonalScale( LEFT, NORMAL, scales, dRowA );
        DiagonalSolve( LEFT, NORMAL, scales, A );

        RowMaxNorms( B, scales );
        EntrywiseMap( scales, MakeFunction(DampScaling<Real>) );
        DiagonalScale( LEFT, NORMAL, scales, dRowB );
        DiagonalSolve( LEFT, NORMAL, scales, B );
    }
    SetIndent( indent );
}

Int ZeroNorm( const DistSparseMatrix<T>& A, Base<T> tol )
{
    DEBUG_ONLY(CSE cse("ZeroNorm"))
    Int numNonzeros = 0;
    const Int numLocalEntries = A.NumLocalEntries();
    for( Int k=0; k<numLocalEntries; ++k )
        if( Abs(A.Value(k)) > tol )
            ++numNonzeros;
    return mpi::AllReduce( numNonzeros, A.Comm() );
}

void ShiftDiagonal
( DistSparseMatrix<T>& A, S alphaPre, Int offset, bool existingDiag )
{
    EL_DEBUG_CSE
    const Int mLocal = A.LocalHeight();
    const Int n = A.Width();
    const T alpha = T(alphaPre);
    if( existingDiag ) 
    {
        T* valBuf = A.ValueBuffer();
        for( Int iLoc=0; iLoc<mLocal; ++iLoc )
        {
            const Int i = A.GlobalRow(iLoc);
            const Int e = A.Offset( iLoc, i+offset );
            valBuf[e] += alpha;
        } 
    }
    else
    {
        A.Reserve( mLocal );
        for( Int iLoc=0; iLoc<mLocal; ++iLoc )
        {
            const Int i = A.GlobalRow(iLoc);
            if( i+offset >= 0 && i+offset < n )
                A.QueueLocalUpdate( iLoc, i+offset, alpha );
        }
        A.ProcessLocalQueues();
    }
}

void Helmholtz( DistSparseMatrix<F>& H, Int nx, Int ny, F shift )
{
    DEBUG_CSE
    typedef Base<F> Real;
    const Int n = nx*ny;
    Zeros( H, n, n );

    const Real hxInv = nx+1;
    const Real hyInv = ny+1;
    const Real hxInvSquared = hxInv*hxInv;
    const Real hyInvSquared = hyInv*hyInv;
    const F mainTerm = 2*(hxInvSquared+hyInvSquared) - shift;

    const Int localHeight = H.LocalHeight();
    H.Reserve( 5*localHeight );
    for( Int iLoc=0; iLoc<localHeight; ++iLoc )
    {
        const Int i = H.GlobalRow(iLoc);
        const Int x = i % nx;
        const Int y = i/nx;

        H.QueueUpdate( i, i, mainTerm );
        if( x != 0 )
            H.QueueUpdate( i, i-1, -hxInvSquared );
        if( x != nx-1 )
            H.QueueUpdate( i, i+1, -hxInvSquared );
        if( y != 0 )
            H.QueueUpdate( i, i-nx, -hyInvSquared );
        if( y != ny-1 )
            H.QueueUpdate( i, i+nx, -hyInvSquared );
    }
    H.ProcessQueues();
}

Int NumOutside( const DistSparseMatrix<Real>& A )
{
    EL_DEBUG_CSE
    const Int numLocalEntries = A.NumLocalEntries();
    const Real* valBuf = A.LockedValueBuffer();

    Int numLocalNonPos = 0;
    for( Int k=0; k<numLocalEntries; ++k )
        if( valBuf[k] <= Real(0) )
            ++numLocalNonPos;

    return mpi::AllReduce( numLocalNonPos, A.Comm() );
}

void KKT
( const DistSparseMatrix<Real>& A,
  const DistSparseMatrix<Real>& G,
  const DistMultiVec<Real>& s,
  const DistMultiVec<Real>& z,
        DistSparseMatrix<Real>& J,
  bool onlyLower )
{
    EL_DEBUG_CSE
    const Int n = A.Width();
    DistSparseMatrix<Real> Q(A.Grid());
    Q.Resize( n, n );
    qp::affine::KKT( Q, A, G, s, z, J, onlyLower );
}

void AugmentedKKT
( const DistSparseMatrix<Real>& A,
        Real gamma,
        Real delta,
  const DistMultiVec<Real>& x,
  const DistMultiVec<Real>& z,
        DistSparseMatrix<Real>& J,
  bool onlyLower )
{
    EL_DEBUG_CSE
    const Int n = A.Width();
    DistSparseMatrix<Real> Q(A.Comm());
    Zeros( Q, n, n );
    qp::direct::AugmentedKKT( Q, A, gamma, delta, x, z, J, onlyLower );
}

Base<Field>
MinAbsNonzero( const DistSparseMatrix<Field>& A, Base<Field> upperBound )
{
    EL_DEBUG_CSE
    typedef Base<Field> Real;
    const Int numEntries = A.NumLocalEntries();
    Real minLocAbs = upperBound;
    for( Int e=0; e<numEntries; ++e )
    {
        const Real absVal = Abs(A.Value(e));
        if( absVal > Real(0) )
            minLocAbs = Min(minLocAbs,absVal);
    }
    return mpi::AllReduce( minLocAbs, mpi::MIN, A.Comm() );
}

void Mehrotra
( const DistSparseMatrix<Real>& A,
  const DistMultiVec<Real>& b,
  const DistMultiVec<Real>& c,
  const DistMultiVec<Int>& orders,
  const DistMultiVec<Int>& firstInds,
        DistMultiVec<Real>& x,
        DistMultiVec<Real>& y,
        DistMultiVec<Real>& z,
  const MehrotraCtrl<Real>& ctrl )
{
    EL_DEBUG_CSE
    const Int n = c.Height();
    const Grid& grid = A.Grid();

    DistSparseMatrix<Real> G(grid);
    Identity( G, n, n );
    G *= -1;

    DistMultiVec<Real> h(grid);
    Zeros( h, n, 1 );

    MehrotraCtrl<Real> affineCtrl = ctrl;
    affineCtrl.primalInit = false;
    affineCtrl.dualInit = false;

    DistMultiVec<Real> s(grid);
    socp::affine::Mehrotra(A,G,b,c,h,orders,firstInds,x,y,z,s,affineCtrl);
}

pair<Base<F>,Base<F>>
HermitianHelper( const DistSparseMatrix<F>& A, Int basisSize )
{
    typedef Base<F> Real;
    Grid grid( A.Comm() );

    DistMatrix<Real,STAR,STAR> T(grid);
    Lanczos( A, T, basisSize );
    const Int k = T.Height();
    if( k == 0 )
        return pair<Real,Real>(0,0);

    auto d = GetDiagonal( T.Matrix() );
    auto dSub = GetDiagonal( T.Matrix(), -1 );

    Matrix<Real> w;
    HermitianTridiagEig( d, dSub, w );

    pair<Real,Real> extremal;
    extremal.second = MaxNorm(w);
    extremal.first = extremal.second;
    for( Int i=0; i<k; ++i )
        extremal.first = Min(extremal.first,Abs(w.Get(i,0)));

    return extremal;
}

void Lanczos
( const DistSparseMatrix<F>& A,
        ElementalMatrix<Base<F>>& T,
        Int basisSize )
{
    DEBUG_CSE
    const Int n = A.Height();
    if( n != A.Width() )
        LogicError("A was not square");

    auto applyA =
      [&]( const DistMultiVec<F>& X, DistMultiVec<F>& Y )
      {
          Zeros( Y, n, X.Width() );
          Multiply( NORMAL, F(1), A, X, F(0), Y );
      };
    Lanczos<F>( n, applyA, T, basisSize );
}

void LAV
( const DistSparseMatrix<Real>& A,
  const DistMultiVec<Real>& b,
        DistMultiVec<Real>& x,
  const lp::affine::Ctrl<Real>& ctrl )
{
    EL_DEBUG_CSE
    const Int m = A.Height();
    const Int n = A.Width();
    const Grid& grid = A.Grid();

    DistSparseMatrix<Real> AHat(grid), G(grid);
    DistMultiVec<Real> c(grid), h(grid);

    // c := [0;1;1]
    // ============
    Zeros( c, n+2*m, 1 );
    for( Int iLoc=0; iLoc<c.LocalHeight(); ++iLoc )
        if( c.GlobalRow(iLoc) >= n )
            c.SetLocal( iLoc, 0, Real(1) );

    // \hat A := [A, I, -I]
    // ====================
    Zeros( AHat, m, n+2*m );
    const Int numLocalEntriesA = A.NumLocalEntries();
    AHat.Reserve( numLocalEntriesA + 2*AHat.LocalHeight() );
    for( Int e=0; e<numLocalEntriesA; ++e )
        AHat.QueueUpdate( A.Row(e), A.Col(e), A.Value(e) );
    for( Int iLoc=0; iLoc<AHat.LocalHeight(); ++iLoc )
    {
        const Int i = AHat.GlobalRow(iLoc);
        AHat.QueueLocalUpdate( iLoc, i+n,   Real( 1) );
        AHat.QueueLocalUpdate( iLoc, i+n+m, Real(-1) );
    }
    AHat.ProcessLocalQueues();

    // G := | 0 -I  0 |
    //      | 0  0 -I |
    // ================
    Zeros( G, 2*m, n+2*m );
    G.Reserve( G.LocalHeight() );
    for( Int iLoc=0; iLoc<G.LocalHeight(); ++iLoc )
        G.QueueLocalUpdate( iLoc, G.GlobalRow(iLoc)+n, Real(-1) );
    G.ProcessLocalQueues();

    // h := | 0 |
    //      | 0 |
    // ==========
    Zeros( h, 2*m, 1 );

    // Solve the affine QP
    // ===================
    DistMultiVec<Real> xHat(grid), y(grid), z(grid), s(grid);
    LP( AHat, G, b, c, h, xHat, y, z, s, ctrl );

    // Extract x
    // =========
    x = xHat( IR(0,n), ALL );
}