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

void Helmholtz( SparseMatrix<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;

    H.Reserve( 5*n );
    for( Int i=0; i<n; ++i )
    {
        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();
}

void ShiftDiagonal
( SparseMatrix<T>& A, S alphaPre, Int offset, bool existingDiag )
{
    EL_DEBUG_CSE
    const Int m = A.Height();
    const Int n = A.Width();
    const T alpha = T(alphaPre);
    if( existingDiag )
    {
        T* valBuf = A.ValueBuffer();
        for( Int i=Max(0,-offset); i<Min(m,n-offset); ++i )
        {
            const Int e = A.Offset( i, i+offset );
            valBuf[e] += alpha;
        } 
    }
    else
    {
        const Int diagLength = Min(m,n-offset) - Max(0,-offset);
        A.Reserve( diagLength );
        for( Int i=Max(0,-offset); i<Min(m,n-offset); ++i )
            A.QueueUpdate( i, i+offset, alpha );
        A.ProcessQueues();
    }
}

void Fill( SparseMatrix<T>& A, T alpha )
{
    EL_DEBUG_CSE
    const Int m = A.Height();
    const Int n = A.Width();
    A.Resize( m, n );
    Zero( A );
    if( alpha != T(0) )
    {
        A.Reserve( m*n ); 
        for( Int i=0; i<m; ++i )
            for( Int j=0; j<n; ++j )
                A.QueueUpdate( i, j, alpha );
        A.ProcessQueues();
    }
}

void JordanCholesky( SparseMatrix<T>& A, Int n )
{
    DEBUG_ONLY(CSE cse("JordanCholesky"))
    Zeros( A, n, n );
    A.Reserve( 3*n );
    
    for( Int e=0; e<n; ++e )
    {
        if( e == 0 )
            A.QueueUpdate( e, e, T(1) );
        else
            A.QueueUpdate( e, e, T(5) );
        if( e > 0 )
            A.QueueUpdate( e, e-1, T(2) );
        if( e < n-1 )
            A.QueueUpdate( e, e+1, T(2) );
    }
    A.ProcessQueues();
}

bool LoadMatrixMarketFile(const std::string& file_path, 
                          SparseMatrix<T>& A,
                          unsigned int& height,
                          unsigned int& width,
                          unsigned int& nnz)
{
    std::ifstream infile(file_path);
    if (!infile)
        return false;

    char mm_typecode[4];

    // read the matrix market banner (header)
    if (0 != mm_read_banner(infile, mm_typecode))
        return false;

    if (!mm_is_valid(mm_typecode))
        return false;

    // this reader supports these matrix types:
    //
    //  sparse, real/integer/pattern, general/symm/skew
    //

    if (!mm_is_sparse(mm_typecode))
    {
        std::cerr << "Only sparse MatrixMarket files are supported." << std::endl;
        return false;
    }

    if (!mm_is_real(mm_typecode) && !mm_is_integer(mm_typecode) && !mm_is_pattern(mm_typecode))
    {
        std::cerr << "Only real, integer, and pattern MatrixMarket formats are supported." << std::endl;
        return false;
    }

    if (!mm_is_general(mm_typecode) && !mm_is_symmetric(mm_typecode) && !mm_is_skew(mm_typecode))
    {
        std::cerr << "Only general, symmetric, and skew-symmetric MatrixMarket formats are supported." 
                  << std::endl;
        return false;
    }

    // read the number of rows, cols, nonzeros
    if (0 != mm_read_mtx_crd_size(infile, height, width, nnz))
    {
        std::cerr << "could not read matrix coordinate information" << std::endl;
        height = width = nnz = 0;
        return false;
    }

    // read the data according to the type 

    bool is_real      = mm_is_real(mm_typecode);
    bool is_int       = mm_is_integer(mm_typecode);
    bool is_symmetric = mm_is_symmetric(mm_typecode);
    bool is_skew      = mm_is_skew(mm_typecode);

    std::string line;
    unsigned int reserve_size = nnz;
    if (is_symmetric || is_skew)
        reserve_size *= 2;

    A.Clear();
    A.Reserve(height, width, reserve_size);

    // load num random entries of A
    A.BeginLoad();
    
    unsigned int row, col, count;

    if (is_real)
    {
        double val;
        for (count=0; count != nnz; ++count)
        {
            infile >> row; assert(row >= 1);
            infile >> col; assert(col >= 1);
            infile >> val;
            
            // convert to 0-based indexing
            row -= 1;
            col -= 1;
            A.Load(row, col, val);

            if (row != col)
            {
                if (is_symmetric)
                    A.Load(col, row, val);
                else if (is_skew)
                    A.Load(col, row, -val);
            }
        }
    }
    else if (is_int)

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

void SparseMatrix<T>::SubMatrixColsCompact(SparseMatrix<T>& result,
                                           const std::vector<unsigned int>& col_indices,
                                           std::vector<unsigned int>& old_to_new_rows,
                                           std::vector<unsigned int>& new_to_old_rows) const
{
    const unsigned int UNUSED_ROW = 0xFFFFFFFF;

    // extract entire columns from the source matrix to form the dest matrix

    unsigned int new_width = col_indices.size();
    if (0u == new_width)
        throw std::logic_error("SparseMatrix::SubMatrixColsCompact: empty column set");
    
    // need one entry per original row in the row_map
    if (old_to_new_rows.size() < height_)
        old_to_new_rows.resize(height_);

    std::fill(old_to_new_rows.begin(), old_to_new_rows.begin() + height_, UNUSED_ROW);
    // no need to fill 'new_to_old_rows'

    // check the column indices for validty and count nonzeros

    unsigned int new_size = 0;
    for (auto it=col_indices.begin(); it != col_indices.end(); ++it)
    {
        // index of next source column
        unsigned int c = *it;
        if (c >= width_)
            throw std::logic_error("SparseMatrix::SubMatrixColsCompact: column index out of range");
        
        // number of nonzeros in source column c
        new_size += (col_offsets_[c+1] - col_offsets_[c]);
    }

    if (0u == new_size)
        throw std::logic_error("SparseMatrix::SubMatrixColsCompact: submatrix is the zero matrix");

    // allocate memory in the result; won't allocate if sufficient memory available
    result.Reserve(height_, new_width, new_size);

    unsigned int* cols_b = result.ColBuffer();
    unsigned int* rows_b = result.RowBuffer();
    T* data_b            = result.DataBuffer();

    unsigned int c_dest = 0;
    unsigned int elt_count = 0;
    for (auto it=col_indices.begin(); it != col_indices.end(); ++it)
    {
        // index of the next source column
        unsigned int c = *it;

        // set element offset for the next dest column
        cols_b[c_dest] = elt_count;

        unsigned int start = col_offsets_[c];
        unsigned int end   = col_offsets_[c+1];
        for (unsigned int offset = start; offset != end; ++offset)
        {
            unsigned int row = row_indices_[offset];
            old_to_new_rows[row] = row;
            rows_b[elt_count] = row;
            data_b[elt_count] = data_[offset];
            ++elt_count;
        }

        // have now completed another column
        ++c_dest;
    }

    cols_b[new_width] = elt_count;
    assert(new_width == c_dest);

    // set the size of the new matrix explicitly, since Load() has been bypassed
    result.SetSize(elt_count);

    // determine the new height of the submatrix
    unsigned int new_height = 0;
    for (unsigned int r=0; r != height_; ++r)
    {
        if (UNUSED_ROW != old_to_new_rows[r])
            ++new_height;
    }
    
    new_to_old_rows.resize(new_height);

    // renumber the rows in the submatrix

    unsigned int new_r = 0;
    for (unsigned int r=0; r<height_; ++r)
    {
        if (UNUSED_ROW != old_to_new_rows[r])
        {
            old_to_new_rows[r] = new_r;
            new_to_old_rows[new_r] = r;
            ++new_r;
        }
    }

    // set the height of the new matrix explicitly
    assert(new_r == new_height);
//.........这里部分代码省略.........