Merge lp:~jobh/dolfin/blockmatrix into lp:dolfin/1.0.x

[Joachim Haga]

New functionality that I use in my current block preconditioning / block LA work

[Garth Wells (garth-wells)]

I'm not keen on the use of GenericMatrix::getrow (I want to remove this function) because it's ambiguous in parallel.

[Joachim Haga (jobh)]

Hi Garth, I'll look into it later this week. I can easily get rid of use
of GM::setrow, but I'll need functionality similar to GM::getrow
(whether or not this actual interface is used). Would it be helpful if
GM::setrow is excised even if GM::getrow is still used? It seems that
the latter is better supported from the backends in parallel.

I imagine at a later stage this can be changed to work on the element
tensor level, but that is a larger change.

-j.

On 01/31/2011 07:12 PM, Garth Wells wrote:
> I'm not keen on the use of GenericMatrix::getrow (I want to remove this function) because it's ambiguous in parallel.

=== modified file 'dolfin/fem/DirichletBC.cpp'
--- dolfin/fem/DirichletBC.cpp	2011-01-01 22:26:49 +0000
+++ dolfin/fem/DirichletBC.cpp	2011-01-28 22:59:27 +0000
@@ -209,6 +209,82 @@
   A.apply("insert");
 }
 //-----------------------------------------------------------------------------
+void DirichletBC::zero_columns(GenericMatrix& A, GenericVector& b, double diag_val) const
+{
+  Map bv_map;
+  get_boundary_values(bv_map, _method);
+
+  // Create lookup table of dofs
+  //const uint nrows = A.size(0); // should be equal to b.size()
+  const uint ncols = A.size(1); // should be equal to max possible dof+1
+
+  std::pair<uint,uint> rows = A.local_range(0);
+  //std::pair<uint,uint> cols = A.local_range(1);
+
+  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<uint>   cols;
+  std::vector<double> vals;
+  std::vector<double> b_vals;
+  std::vector<uint>   b_rows;
+
+  for (uint 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 (uint 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 (uint j=0; j<cols.size(); j++)
+      {
+        const uint 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(&b_vals.front(), b_rows.size(), &b_rows.front());
+  b.apply("add");
+}
+//-----------------------------------------------------------------------------
 const std::vector<std::pair<dolfin::uint, dolfin::uint> >& DirichletBC::markers()
 {
   return facets;

=== modified file 'dolfin/fem/DirichletBC.h'
--- dolfin/fem/DirichletBC.h	2010-12-31 21:41:41 +0000
+++ dolfin/fem/DirichletBC.h	2011-01-28 22:59:27 +0000
@@ -169,6 +169,11 @@
     /// non-diagonal matrices in a block matrix.
     void zero(GenericMatrix& A) const;
 
+    /// Make columns associated with boundary conditions zero, and
+    /// update the RHS to reflect the changes. Useful for non-diagonals.
+    /// The diag_val parameter would normally be -1, 0 or 1.
+    void zero_columns(GenericMatrix& A, GenericVector& b, double diag_val=0) const;
+
     /// Return boundary markers (facets stored as pairs of cells and local
     /// facet numbers)
     const std::vector<std::pair<uint, uint> >& markers();

=== modified file 'dolfin/la/BlockMatrix.cpp'
--- dolfin/la/BlockMatrix.cpp	2011-01-22 23:04:29 +0000
+++ dolfin/la/BlockMatrix.cpp	2011-01-28 22:59:27 +0000
@@ -132,3 +132,34 @@
   }
 }
 //-----------------------------------------------------------------------------
+boost::shared_ptr<GenericMatrix> BlockMatrix::schur_approximation(double symmetry) const
+{
+  // Currently returns [diag(C * diag(A)^-1 * B) - D]
+  if (symmetry==0)
+  {
+    error("only implemented for symmetry != 0");
+  }
+  assert(matrices.size()==2 && matrices[0].size()==2 && matrices[1].size()==2);
+
+  GenericMatrix &A = *matrices[0][0];
+  GenericMatrix &C = *matrices[1][0];
+  GenericMatrix &D = *matrices[1][1];
+
+  boost::shared_ptr<GenericMatrix> S(D.copy());
+
+  std::vector<uint> cols_i;
+  std::vector<double> vals_i;
+  for (uint i=0; i<D.size(0); i++)
+  {
+    C.getrow(i, cols_i, vals_i);
+    double diag_ii=0;
+    for (uint k=0; k<cols_i.size(); k++)
+    {
+      const uint j=cols_i[k];
+      const double val=vals_i[k];
+      diag_ii -= val*val/A(j,j);
+    }
+    S->add(&diag_ii, 1, &i, 1, &i);
+  }
+  return S;
+}

=== modified file 'dolfin/la/BlockMatrix.h'
--- dolfin/la/BlockMatrix.h	2011-01-22 23:04:29 +0000
+++ dolfin/la/BlockMatrix.h	2011-01-28 22:59:27 +0000
@@ -53,6 +53,10 @@
     /// Matrix-vector product, y = Ax
     void mult(const BlockVector& x, BlockVector& y, bool transposed=false) const;
 
+    /// Create a crude explicit Schur approximation of S = D - C A^-1 B of (A B; C D)
+    /// If symmetry != 0, then the caller promises that B = symmetry * transpose(C).
+    boost::shared_ptr<GenericMatrix> schur_approximation(double symmetry=1) const;
+
   private:
 
     boost::multi_array<boost::shared_ptr<GenericMatrix>, 2> matrices;

=== modified file 'dolfin/la/EpetraMatrix.cpp'
--- dolfin/la/EpetraMatrix.cpp	2011-01-24 14:28:43 +0000
+++ dolfin/la/EpetraMatrix.cpp	2011-01-28 22:59:27 +0000
@@ -469,7 +469,15 @@
 void EpetraMatrix::setrow(uint row, const std::vector<uint>& columns,
                           const std::vector<double>& values)
 {
-  dolfin_not_implemented();
+  static bool print_msg_once=true;
+  if (print_msg_once)
+  {
+    info("EpetraMatrix::setrow is implemented inefficiently");
+    print_msg_once = false;
+  }
+
+  for (uint i=0; i<columns.size(); i++)
+    set(&values[i], 1, &row, 1, &columns[i]);
 }
 //-----------------------------------------------------------------------------
 LinearAlgebraFactory& EpetraMatrix::factory() const