Implementation and test for matfree zeroRowsColumns #4
Conversation
| @@ -9,7 +9,6 @@ | |||
| from firedrake.petsc import PETSc | |||
| from firedrake.utils import cached_property | |||
|
|
|||
|
|
|||
There was a problem hiding this comment.
Don't change things if you don't have to.
| @@ -57,6 +56,88 @@ def find_sub_block(iset, ises): | |||
| return found | |||
|
|
|||
|
|
|||
| def find_element_of_which_sub_block(rows, ises): | |||
| # This function acts as lookup to find which block the indices belong in | |||
There was a problem hiding this comment.
Needs more explanation in the comment. Try to convey the bigger picture.
There was a problem hiding this comment.
Write what it takes in, what it returns, and why.
| for i in range(step): | ||
| rows_i = numpy.extract(belong_to_sub == i, rows) | ||
| rows_i = rows_i - i | ||
| rows_i = numpy.int32(rows_i/step) |
There was a problem hiding this comment.
Is this safe if you have more than 4 billion dofs? Strange to hardcode int32.
| @@ -106,6 +187,8 @@ def __init__(self, a, row_bcs=[], col_bcs=[], | |||
| self._y = function.Function(test_space) | |||
| self._x = function.Function(trial_space) | |||
|
|
|||
| # Temporary storage for holding the BC values during zeroRowsColumns | |||
| self._tmp_zeroRowsColumns = function.Function(test_space) | |||
There was a problem hiding this comment.
Functions are, on average, expensive. (Imagine you have a problem with a billion dofs.) Can we make this so that this is only allocated if needed?
| if active_rows.size == 0: | ||
| # This can happen when running in parallel. Every processor sharing the matrix | ||
| # must call zeroRowsColumns but the processor might not have any rows to zero. | ||
| # For now we just return without adding additional DirichletBC, is this the correct thing to do? |
There was a problem hiding this comment.
I think that bc application is collective -- you should make a DirichletBC object on every processor, even if it's empty.
| # must call zeroRowsColumns but the processor might not have any rows to zero. | ||
| # For now we just return without adding additional DirichletBC, is this the correct thing to do? | ||
| return | ||
| if not numpy.allclose(diag, 1.0): |
There was a problem hiding this comment.
What's wrong with if diag != 1?
| if test.function_space().shape == (): | ||
| rhs_form = inner(Constant(1), test)*dx | ||
| elif test.function_space().shape == (2, ): | ||
| rhs_form = inner(Constant((1, 1)), test)*dx |
There was a problem hiding this comment.
else: raise NotImplementedError
This is an implementation for a matfree
zeroRowsColumnsmethod. Since we are dealing with unassembled matrices, we cannot perform the algebraic operation of zeroing the row and column. So instead whenpetsc4pypasses along which row/column, needs to be zeroed, we find its associated dof and make aZeroRowsColumnsBCobject (which overloads theDirichletBCclass) to impose the same effect as zeroing the rows and column of that dof and adding a 1 onto the diagonal. Tests show that this behaves as expected.