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

void QP
( const SparseMatrix<Real>& A, 
  const Matrix<Real>& B, 
        Matrix<Real>& X, 
  const qp::direct::Ctrl<Real>& ctrl )
{
    DEBUG_CSE

    const Int n = A.Width();
    const Int k = B.Width();
    SparseMatrix<Real> Q, AHat;
    Matrix<Real> bHat, c;

    Herk( LOWER, ADJOINT, Real(1), A, Q );
    MakeHermitian( LOWER, Q );
    Zeros( AHat, 0, n );
    Zeros( bHat, 0, 1 );
    Zeros( X,    n, k );

    Matrix<Real> y, z;
    for( Int j=0; j<k; ++j )
    {
        auto x = X( ALL, IR(j) );
        auto b = B( ALL, IR(j) );

        Zeros( c, n, 1 );
        Multiply( ADJOINT, Real(-1), A, b, Real(0), c );

        El::QP( Q, AHat, bHat, c, x, y, z, ctrl );
    }
}

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

T OneNorm(const SparseMatrix<T>& A)
{
    // compute the max absolute column sum
    const unsigned int* cols_a = A.LockedColBuffer();
    const T*            data_a = A.LockedDataBuffer();
    const unsigned int width_a = A.Width();

    T max_col_sum = T(0), col_sum;
    for (unsigned int c=0; c != width_a; ++c)
    {
        unsigned int start = cols_a[c];
        unsigned int   end = cols_a[c+1];
        col_sum = T(0);
        for (unsigned int offset=start; offset != end; ++offset)
        {
            T val = fabs(data_a[offset]);
            col_sum += val;
        }

        if (col_sum > max_col_sum)
            max_col_sum = col_sum;
    }

    return max_col_sum;
}

bool WriteMatrixMarketFile(const std::string& file_path,
                           const SparseMatrix<T>& A,
                           const unsigned int precision)
{
    // Write a MatrixMarket file with no comments.  Note that the
    // MatrixMarket format uses 1-based indexing for rows and columns.

    std::ofstream outfile(file_path);
    if (!outfile)
        return false;

    unsigned int height = A.Height();
    unsigned int width  = A.Width();
    unsigned int nnz    = A.Size();
    
    // write the 'banner'
    outfile << MM_BANNER << " matrix coordinate real general" << std::endl;

    // write matrix dimensions and number of nonzeros
    outfile << height << " " << width << " " << nnz << std::endl;

    outfile << std::fixed;
    outfile.precision(precision);
    
    const unsigned int* cols_a = A.LockedColBuffer();
    const unsigned int* rows_a = A.LockedRowBuffer();
    const T*            data_a = A.LockedDataBuffer();
    unsigned int width_a = A.Width();

    for (unsigned int c=0; c != width_a; ++c)
    {
        unsigned int start = cols_a[c];
        unsigned int end   = cols_a[c+1];
        for (unsigned int offset=start; offset != end; ++offset)
        {
            unsigned int r = rows_a[offset];
            T val = data_a[offset];
            outfile << r+1 << " " << c+1 << " " << val << std::endl;
        }
    }

    outfile.close();
    return true;
}

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 Print(const SparseMatrix<T>& M)
{
    // Print a SparseMatrix to the screen.

    const unsigned int* col_buf = M.LockedColBuffer();
    const unsigned int* row_buf = M.LockedRowBuffer();
    const T*                buf = M.LockedDataBuffer();

    if (0 == M.Size())
    {
        std::cout << "Matrix is empty." << std::endl;
        return;
    }

    for (unsigned int c=0; c != M.Width(); ++c)
    {
        unsigned int start = col_buf[c];
        unsigned int end   = col_buf[c+1];
        for (unsigned int offset=start; offset != end; ++offset)
        {
            assert(offset >= 0);
            assert(offset < M.Size());
            unsigned int row_index = row_buf[offset];
            T                 data = buf[offset];
            std::cout << "(" << row_index << ", " << c << "): " << data << std::endl;
        }
    }

    std::cout << "Col indices: "; std::cout.flush();
    for (unsigned int i=0; i != M.Width(); ++i)
        std::cout << col_buf[i] << ", ";
    std::cout << col_buf[M.Width()] << std::endl;

    std::cout << "Row indices: "; std::cout.flush();
    for (unsigned int i=0; i != M.Size(); ++i)
        std::cout << row_buf[i] << ", ";
    std::cout << std::endl;

    std::cout << "Data:        "; std::cout.flush();
    for (unsigned int i=0; i != M.Size(); ++i)
        std::cout << buf[i] << ", ";
    std::cout << std::endl;
}

void Tikhonov
( Orientation orientation,
  const SparseMatrix<F>& A,
  const Matrix<F>& B,
  const SparseMatrix<F>& G,
        Matrix<F>& X, 
  const LeastSquaresCtrl<Base<F>>& ctrl )
{
    DEBUG_CSE
    
    // Explicitly form W := op(A)
    // ==========================
    SparseMatrix<F> W;
    if( orientation == NORMAL )
        W = A;
    else if( orientation == TRANSPOSE )
        Transpose( A, W );
    else
        Adjoint( A, W );

    const Int m = W.Height();
    const Int n = W.Width();
    const Int numRHS = B.Width();

    // Embed into a higher-dimensional problem via appending regularization
    // ====================================================================
    SparseMatrix<F> WEmb;
    if( m >= n )
        VCat( W, G, WEmb ); 
    else
        HCat( W, G, WEmb );
    Matrix<F> BEmb;
    Zeros( BEmb, WEmb.Height(), numRHS );
    if( m >= n )
    {
        auto BEmbT = BEmb( IR(0,m), IR(0,numRHS) );
        BEmbT = B;
    }
    else
        BEmb = B;

    // Solve the higher-dimensional problem
    // ====================================
    Matrix<F> XEmb;
    LeastSquares( NORMAL, WEmb, BEmb, XEmb, ctrl );

    // Extract the solution
    // ====================
    if( m >= n )
        X = XEmb;
    else
        X = XEmb( IR(0,n), IR(0,numRHS) ); 
}

void StackedRuizEquil
( SparseMatrix<Field>& A,
  SparseMatrix<Field>& B,
  Matrix<Base<Field>>& dRowA,
  Matrix<Base<Field>>& dRowB,
  Matrix<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();
    Ones( dRowA, mA, 1 );
    Ones( dRowB, mB, 1 );
    Ones( dCol, n, 1 );

    // TODO(poulson): Expose these as control parameters
    // For now, simply hard-code the number of iterations
    const Int maxIter = 4;

    Matrix<Real> scales, maxAbsValsB;
    const Int indent = PushIndent();
    for( Int iter=0; iter<maxIter; ++iter )
    {
        // Rescale the columns
        // -------------------
        ColumnMaxNorms( A, scales );
        ColumnMaxNorms( B, maxAbsValsB );
        for( Int j=0; j<n; ++j )
            scales(j) = Max(scales(j),maxAbsValsB(j));
        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 );
}

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

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

void KKT
( const SparseMatrix<Real>& A,
        Real gamma,
        Real delta,
        Real beta,
  const Matrix<Real>& x,
  const Matrix<Real>& z,
        SparseMatrix<Real>& J, bool onlyLower )
{
    EL_DEBUG_CSE
    const Int n = A.Width();
    SparseMatrix<Real> Q;
    Q.Resize( n, n );
    qp::direct::KKT( Q, A, gamma, delta, beta, x, z, J, onlyLower );
}

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 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 Lanczos
( const SparseMatrix<F>& A,
        Matrix<Base<F>>& T,
        Int basisSize )
{
    DEBUG_CSE
    const Int n = A.Height();
    if( n != A.Width() )
        LogicError("A was not square");
    
    auto applyA =
      [&]( const Matrix<F>& X, Matrix<F>& Y )
      {
          Zeros( Y, n, X.Width() );
          Multiply( NORMAL, F(1), A, X, F(0), Y );
      };
    Lanczos<F>( n, applyA, T, basisSize );
}

Base<Field> ProductLanczosDecomp
( const SparseMatrix<Field>& A,
        Matrix<Field>& V,
        Matrix<Base<Field>>& T,
        Matrix<Field>& v,
        Int basisSize )
{
    EL_DEBUG_CSE
    const Int m = A.Height();
    const Int n = A.Width();

    Matrix<Field> S;

    // Cache the adjoint
    // -----------------
    SparseMatrix<Field> AAdj;
    Adjoint( A, AAdj );

    if( m >= n )
    {
        auto applyA =
          [&]( const Matrix<Field>& X, Matrix<Field>& Y )
          {
              Zeros( S, m, X.Width() );
              Multiply( NORMAL, Field(1), A,    X, Field(0), S );
              Zeros( Y, n, X.Width() );
              Multiply( NORMAL, Field(1), AAdj, S, Field(0), Y );
          };
        return LanczosDecomp( n, applyA, V, T, v, basisSize );
    }
    else
    {
        auto applyA =
          [&]( const Matrix<Field>& X, Matrix<Field>& Y )
          {
              Zeros( S, n, X.Width() );
              Multiply( NORMAL, Field(1), AAdj, X, Field(0), S );
              Zeros( Y, m, X.Width() );
              Multiply( NORMAL, Field(1), A,    S, Field(0), Y );
          };
        return LanczosDecomp( m, applyA, V, T, v, basisSize );
    }
}