C++ GenericMatrix::size方法代码示例

//---------------------------------------------------------------------------
void EigenMatrix::axpy(double a, const GenericMatrix& A,
                       bool same_nonzero_pattern)
{
  // Check for same size
  if (size(0) != A.size(0) or size(1) != A.size(1))
  {
    dolfin_error("EigenMatrix.cpp",
                 "perform axpy operation with Eigen matrix",
                 "Dimensions don't match");
  }

  _matA += (a)*(as_type<const EigenMatrix>(A).mat());
}

//-----------------------------------------------------------------------------
CoordinateMatrix::CoordinateMatrix(const GenericMatrix& A, bool symmetric,
                                   bool base_one)
  : _symmetric(symmetric), _base_one(base_one)
{
  _size[0] = A.size(0);
  _size[1] = A.size(1);

  // Iterate over local rows
  const std::pair<std::size_t, std::size_t> local_row_range = A.local_range(0);
  if (!_symmetric)
  {
    for (std::size_t i = local_row_range.first; i < local_row_range.second; ++i)
    {
      // Get column and value data for row
      std::vector<std::size_t> columns;
      std::vector<double> values;
      A.getrow(i, columns, values);

      // Insert data at end
      _rows.insert(_rows.end(), columns.size(), i);
      _cols.insert(_cols.end(), columns.begin(), columns.end());
      _vals.insert(_vals.end(), values.begin(), values.end());
    }
    assert(_rows.size() == _cols.size());
  }
  else
  {
    assert(_size[0] == _size[1]);
    for (std::size_t i = local_row_range.first; i < local_row_range.second; ++i)
    {
      // Get column and value data for row
      std::vector<std::size_t> columns;
      std::vector<double> values;
      A.getrow(i, columns, values);

      for (std::size_t j = 0; j < columns.size(); ++j)
      {
        if (columns[j] >= i)
        {
          _rows.push_back(i);
          _cols.push_back(columns[j]);
          _vals.push_back(values[j]);
        }
      }
    }
    assert(_rows.size() == _cols.size());
  }

  // Add 1 for Fortran-style indices
  if (base_one)
  {
    for (std::size_t i = 0; i < _cols.size(); ++i)
    {
      _rows[i]++;
      _cols[i]++;
    }
  }
}

//-----------------------------------------------------------------------------
void DirichletBC::zero_columns(GenericMatrix& A,
                               GenericVector& b,
                               double diag_val) const
{
  // Check arguments
  check_arguments(&A, &b, NULL, 1);

  // A map to hold the mapping from boundary dofs to boundary values
  Map bv_map;
  get_boundary_values(bv_map);

  // Create lookup table of dofs
  //const std::size_t nrows = A.size(0); // should be equal to b.size()
  const std::size_t ncols = A.size(1); // should be equal to max possible dof+1

  std::pair<std::size_t, std::size_t> rows = A.local_range(0);

  std::vector<char> is_bc_dof(ncols);
  std::vector<double> bc_dof_val(ncols);
  for (Map::const_iterator bv = bv_map.begin(); bv != bv_map.end(); ++bv)
  {
    is_bc_dof[bv->first] = 1;
    bc_dof_val[bv->first] = bv->second;
  }

  // Scan through all columns of all rows, setting to zero if
  // is_bc_dof[column]. At the same time, we collect corrections to
  // the RHS

  std::vector<std::size_t> cols;
  std::vector<double> vals;
  std::vector<double> b_vals;
  std::vector<dolfin::la_index> b_rows;

  for (std::size_t row = rows.first; row < rows.second; row++)
  {
    // If diag_val is nonzero, the matrix is a diagonal block
    // (nrows==ncols), and we can set the whole BC row
    if (diag_val != 0.0 && is_bc_dof[row])
    {
      A.getrow(row, cols, vals);
      for (std::size_t j = 0; j < cols.size(); j++)
        vals[j] = (cols[j] == row)*diag_val;
      A.setrow(row, cols, vals);
      A.apply("insert");
      b.setitem(row, bc_dof_val[row]*diag_val);
    }
    else // Otherwise, we scan the row for BC columns
    {
      A.getrow(row, cols, vals);
      bool row_changed = false;
      for (std::size_t j = 0; j < cols.size(); j++)
      {
        const std::size_t col = cols[j];

        // Skip columns that aren't BC, and entries that are zero
        if (!is_bc_dof[col] || vals[j] == 0.0)
          continue;

        // We're going to change the row, so make room for it
        if (!row_changed)
        {
          row_changed = true;
          b_rows.push_back(row);
          b_vals.push_back(0.0);
        }

        b_vals.back() -= bc_dof_val[col]*vals[j];
        vals[j] = 0.0;
      }
      if (row_changed)
      {
        A.setrow(row, cols, vals);
        A.apply("insert");
      }
    }
  }

  b.add_local(&b_vals.front(), b_rows.size(), &b_rows.front());
  b.apply("add");
}