C++ SparseMatrix::Col方法代码示例

void LAV
( const SparseMatrix<Real>& A,
  const Matrix<Real>& b,
        Matrix<Real>& x,
  const lp::affine::Ctrl<Real>& ctrl )
{
    EL_DEBUG_CSE
    const Int m = A.Height();
    const Int n = A.Width();
    const Range<Int> xInd(0,n), uInd(n,n+m), vInd(n+m,n+2*m);
    SparseMatrix<Real> AHat, G;
    Matrix<Real> c, h;

    // 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 );
    const Int numEntriesA = A.NumEntries();
    AHat.Reserve( numEntriesA + 2*m );
    for( Int e=0; e<numEntriesA; ++e )
        AHat.QueueUpdate( A.Row(e), A.Col(e), A.Value(e) );
    for( Int i=0; i<m; ++i )
    {
        AHat.QueueUpdate( i, i+n,   Real( 1) );
        AHat.QueueUpdate( i, i+n+m, Real(-1) );
    }
    AHat.ProcessQueues();

    // G := | 0 -I  0 |
    //      | 0  0 -I |
    // ================
    Zeros( G, 2*m, n+2*m );
    G.Reserve( G.Height() );
    for( Int i=0; i<2*m; ++i )
        G.QueueUpdate( i, i+n, Real(-1) );
    G.ProcessQueues();

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

    // Solve the affine linear program
    // ===============================
    Matrix<Real> xHat, y, z, s;
    LP( AHat, G, b, c, h, xHat, y, z, s, ctrl );

    // Extract x
    // =========
    x = xHat( xInd, ALL );
}

void MatrixMarket( const SparseMatrix<T>& A, string basename="matrix" )
{
    EL_DEBUG_CSE
    
    string filename = basename + "." + FileExtension(MATRIX_MARKET);
    ofstream file( filename.c_str(), std::ios::binary );
    if( !file.is_open() )
        RuntimeError("Could not open ",filename);

    // Write the header
    // ================
    {
        ostringstream os;
        os << "%%MatrixMarket matrix coordinate ";
        if( IsComplex<T>::value )
            os << "complex "; 
        else
            os << "real ";
        os << "general\n";
        file << os.str();
    }
    
    // Write the size line
    // ===================
    const Int m = A.Height();
    const Int n = A.Width();
    const Int numNonzeros = A.NumEntries();
    file << BuildString(m," ",n," ",numNonzeros,"\n");
    
    // Write the entries
    // =================
    for( Int e=0; e<numNonzeros; ++e )
    {
        const Int i = A.Row(e);
        const Int j = A.Col(e);
        const T value = A.Value(e);
        if( IsComplex<T>::value )
        {
            file << 
              BuildString(i," ",j," ",RealPart(value)," ",ImagPart(value),"\n");
        }
        else
        {
            file << BuildString(i," ",j," ",RealPart(value),"\n");
        }
    }
}

void IPM
( const SparseMatrix<Real>& A,
  const Matrix<Real>& b, 
        Real lambda,
        Matrix<Real>& x,
  const qp::affine::Ctrl<Real>& ctrl )
{
    DEBUG_CSE
    const Int m = A.Height();
    const Int n = A.Width();
    const Range<Int> uInd(0,n), vInd(n,2*n), rInd(2*n,2*n+m);
    SparseMatrix<Real> Q, AHat, G;
    Matrix<Real> c, h;

    // Q := | 0 0 0 |
    //      | 0 0 0 |
    //      | 0 0 I |
    // ==============
    Zeros( Q, 2*n+m, 2*n+m );
    Q.Reserve( m );
    for( Int e=0; e<m; ++e )
        Q.QueueUpdate( 2*n+e, 2*n+e, Real(1) );
    Q.ProcessQueues();

    // c := lambda*[1;1;0]
    // ===================
    Zeros( c, 2*n+m, 1 );
    auto cuv = c( IR(0,2*n), ALL );
    Fill( cuv, lambda );

    // \hat A := [A, -A, I]
    // ====================
    const Int numEntriesA = A.NumEntries();
    Zeros( AHat, m, 2*n+m );
    AHat.Reserve( 2*numEntriesA+m );
    for( Int e=0; e<numEntriesA; ++e )
    {
        AHat.QueueUpdate( A.Row(e), A.Col(e),    A.Value(e) );
        AHat.QueueUpdate( A.Row(e), A.Col(e)+n, -A.Value(e) );
    }
    for( Int e=0; e<m; ++e )
        AHat.QueueUpdate( e, e+2*n, Real(1) );
    AHat.ProcessQueues();

    // G := | -I  0 0 |
    //      |  0 -I 0 |
    // ================
    Zeros( G, 2*n, 2*n+m );
    G.Reserve( 2*m );
    for( Int e=0; e<2*m; ++e )
        G.QueueUpdate( e, e, Real(-1) );
    G.ProcessQueues();

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

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

    // x := u - v
    // ==========
    x = xHat( uInd, ALL );
    x -= xHat( vInd, ALL );
}

void IPM
( const SparseMatrix<Real>& A,
  const Matrix<Real>& d,
        Real lambda,
        Matrix<Real>& x,
  const qp::affine::Ctrl<Real>& ctrl )
{
    EL_DEBUG_CSE
    const Int m = A.Height();
    const Int n = A.Width();
    const Range<Int> wInd(0,n), betaInd(n,n+1), zInd(n+1,n+m+1);

    SparseMatrix<Real> Q, AHat, G;
    Matrix<Real> c, b, h;

    // Q := | I 0 0 |
    //      | 0 0 0 |
    //      | 0 0 0 |
    // ==============
    Zeros( Q, n+m+1, n+m+1 );
    Q.Reserve( n );
    for( Int e=0; e<n; ++e )
        Q.QueueUpdate( e, e, Real(1) );
    Q.ProcessQueues();

    // c := [0;0;lambda]
    // =================
    Zeros( c, n+m+1, 1 );
    auto cz = c( zInd, ALL );
    Fill( cz, lambda );

    // AHat = []
    // =========
    Zeros( AHat, 0, n+m+1 );

    // b = []
    // ======
    Zeros( b, 0, 1 );

    // G := |-diag(d) A, -d, -I|
    //      |      0,     0, -I|
    // =========================
    Zeros( G, 2*m, n+m+1 );
    const Int numEntriesA = A.NumEntries();
    G.Reserve( numEntriesA+3*m );
    for( Int e=0; e<numEntriesA; ++e )
        G.QueueUpdate( A.Row(e), A.Col(e), -d(A.Row(e))*A.Value(e) );
    for( Int e=0; e<m; ++e )
        G.QueueUpdate( e, n, -d(e) );
    for( Int e=0; e<m; ++e )
    {
        G.QueueUpdate( e,   e+n+1, Real(-1) );
        G.QueueUpdate( e+m, e+n+1, Real(-1) );
    }
    G.ProcessQueues();

    // h := [-ones(m,1); zeros(m,1)]
    // =============================
    Zeros( h, 2*m, 1 );
    auto h0 = h( IR(0,m), ALL );
    Fill( h0, Real(-1) );

    // Solve the affine QP
    // ===================
    Matrix<Real> y, z, s;
    QP( Q, AHat, G, b, c, h, x, y, z, s, ctrl );
}

void RLS
( const SparseMatrix<Real>& A, 
  const Matrix<Real>& b, 
        Real rho,
        Matrix<Real>& x, 
  const socp::affine::Ctrl<Real>& ctrl )
{
    DEBUG_CSE
    const Int m = A.Height();
    const Int n = A.Width();

    Matrix<Int> orders, firstInds;
    Zeros( orders, m+n+3, 1 );
    Zeros( firstInds, m+n+3, 1 );
    for( Int i=0; i<m+1; ++i )
    {
        orders.Set( i, 0, m+1 );
        firstInds.Set( i, 0, 0 );
    }
    for( Int i=0; i<n+2; ++i )
    {
        orders.Set( i+m+1, 0, n+2 ); 
        firstInds.Set( i+m+1, 0, m+1 );
    }

    // G := | -1  0  0 |
    //      |  0  0  A |
    //      |  0 -1  0 |
    //      |  0  0 -I |
    //      |  0  0  0 |
    SparseMatrix<Real> G;
    {
        const Int numEntriesA = A.NumEntries();
        Zeros( G, m+n+3, n+2 );
        G.Reserve( numEntriesA+n+2 );
        G.QueueUpdate( 0, 0, -1 );
        for( Int e=0; e<numEntriesA; ++e )
            G.QueueUpdate( A.Row(e)+1, A.Col(e)+2, A.Value(e) );
        G.QueueUpdate( m+1, 1, -1 );
        for( Int j=0; j<n; ++j )
            G.QueueUpdate( j+m+2, j+2, -1 );
        G.ProcessQueues();
    }

    // h := | 0 |
    //      | b |
    //      | 0 |
    //      | 0 |
    //      | 1 |
    Matrix<Real> h;
    Zeros( h, m+n+3, 1 ); 
    auto hb = h( IR(1,m+1), ALL );
    hb = b;
    h.Set( END, 0, 1 );

    // c := [1; rho; 0]
    Matrix<Real> c;
    Zeros( c, n+2, 1 );
    c.Set( 0, 0, 1 );
    c.Set( 1, 0, rho );

    SparseMatrix<Real> AHat;
    Zeros( AHat, 0, n+2 );
    Matrix<Real> bHat;
    Zeros( bHat, 0, 1 );

    Matrix<Real> xHat, y, z, s;
    SOCP( AHat, G, bHat, c, h, orders, firstInds, xHat, y, z, s, ctrl );
    x = xHat( IR(2,END), ALL );
}