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

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

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 RowMaxNorms( const SparseMatrix<F>& A, Matrix<Base<F>>& norms )
{
    DEBUG_CSE
    typedef Base<F> Real;
    const Int m = A.Height();
    const F* valBuf = A.LockedValueBuffer();
    const Int* offsetBuf = A.LockedOffsetBuffer();

    norms.Resize( m, 1 );
    for( Int i=0; i<m; ++i )
    {
        Real rowMax = 0;
        const Int offset = offsetBuf[i];
        const Int numConn = offsetBuf[i+1] - offset;
        for( Int e=offset; e<offset+numConn; ++e )
            rowMax = Max(rowMax,Abs(valBuf[e]));
        norms(i) = rowMax;
    }
}

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

Base<Field> LanczosDecomp
( const SparseMatrix<Field>& A,
        Matrix<Field>& V,
        Matrix<Base<Field>>& T,
        Matrix<Field>& v,
        Int basisSize )
{
    EL_DEBUG_CSE
    const Int n = A.Height();
    if( n != A.Width() )
        LogicError("A was not square");

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

void RowTwoNorms( const SparseMatrix<F>& A, Matrix<Base<F>>& norms )
{
    DEBUG_CSE
    typedef Base<F> Real;
    const Int m = A.Height();
    const F* valBuf = A.LockedValueBuffer();
    const Int* offsetBuf = A.LockedOffsetBuffer();

    norms.Resize( m, 1 );
    for( Int i=0; i<m; ++i )
    {
        Real scale = 0;
        Real scaledSquare = 1;
        const Int offset = offsetBuf[i];
        const Int numConn = offsetBuf[i+1] - offset;
        for( Int e=offset; e<offset+numConn; ++e )
            UpdateScaledSquare( valBuf[e], scale, scaledSquare );
        norms(i) = scale*Sqrt(scaledSquare);
    }
}

void GetMappedDiagonal
( const SparseMatrix<T>& A,
        Matrix<S>& d,
        function<S(const T&)> func,
        Int offset )
{
    EL_DEBUG_CSE
    const Int m = A.Height();
    const Int n = A.Width();
    const T* valBuf = A.LockedValueBuffer();
    const Int* colBuf = A.LockedTargetBuffer();

    const Int iStart = Max(-offset,0);
    const Int jStart = Max( offset,0);

    const Int diagLength = El::DiagonalLength(m,n,offset);
    d.Resize( diagLength, 1 );
    Zero( d );
    S* dBuf = d.Buffer();

    for( Int k=0; k<diagLength; ++k )
    {
        const Int i = iStart + k;
        const Int j = jStart + k;
        const Int thisOff = A.RowOffset(i);
        const Int nextOff = A.RowOffset(i+1);
        auto it = std::lower_bound( colBuf+thisOff, colBuf+nextOff, j );
        if( *it == j )
        {
            const Int e = it-colBuf;
            dBuf[Min(i,j)] = func(valBuf[e]);
        }
        else
            dBuf[Min(i,j)] = func(0);
    }
}

void SymmetricRuizEquil
( SparseMatrix<F>& A, 
  Matrix<Base<F>>& d,
  Int maxIter, bool progress )
{
    DEBUG_CSE
    typedef Base<F> Real;
    const Int n = A.Height();
    Ones( d, n, 1 );

    Matrix<Real> scales;
    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 );
}

//.........这里部分代码省略.........
   //    corresponding to fespace. Initialize x with initial guess of zero,
   //    which satisfies the boundary conditions.
   GridFunction x(fespace);
   x = 0.0;

   // 9. Set up the bilinear form a(.,.) on the finite element space
   //    corresponding to the linear elasticity integrator with piece-wise
   //    constants coefficient lambda and mu.
   Vector lambda(mesh->attributes.Max());
   lambda = 1.0;
   lambda(0) = lambda(1)*50;
   PWConstCoefficient lambda_func(lambda);
   Vector mu(mesh->attributes.Max());
   mu = 1.0;
   mu(0) = mu(1)*50;
   PWConstCoefficient mu_func(mu);

   BilinearForm *a = new BilinearForm(fespace);
   a->AddDomainIntegrator(new ElasticityIntegrator(lambda_func,mu_func));

   // 10. Assemble the bilinear form and the corresponding linear system,
   //     applying any necessary transformations such as: eliminating boundary
   //     conditions, applying conforming constraints for non-conforming AMR,
   //     static condensation, etc.
   cout << "matrix ... " << flush;
   if (static_cond) { a->EnableStaticCondensation(); }
   a->Assemble();

   SparseMatrix A;
   Vector B, X;
   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
   cout << "done." << endl;

   cout << "Size of linear system: " << A.Height() << endl;

#ifndef MFEM_USE_SUITESPARSE
   // 11. Define a simple symmetric Gauss-Seidel preconditioner and use it to
   //     solve the system Ax=b with PCG.
   GSSmoother M(A);
   PCG(A, M, B, X, 1, 500, 1e-8, 0.0);
#else
   // 11. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
   UMFPackSolver umf_solver;
   umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
   umf_solver.SetOperator(A);
   umf_solver.Mult(B, X);
#endif

   // 12. Recover the solution as a finite element grid function.
   a->RecoverFEMSolution(X, *b, x);

   // 13. For non-NURBS meshes, make the mesh curved based on the finite element
   //     space. This means that we define the mesh elements through a fespace
   //     based transformation of the reference element. This allows us to save
   //     the displaced mesh as a curved mesh when using high-order finite
   //     element displacement field. We assume that the initial mesh (read from
   //     the file) is not higher order curved mesh compared to the chosen FE
   //     space.
   if (!mesh->NURBSext)
   {
      mesh->SetNodalFESpace(fespace);
   }

   // 14. Save the displaced mesh and the inverted solution (which gives the
   //     backward displacements to the original grid). This output can be
   //     viewed later using GLVis: "glvis -m displaced.mesh -g sol.gf".