[petsc-users] How to zeros all entries of a set of columns of a matrix

[Jed Brown]

On Thu, 25 Feb 2010 12:54:13 -0600, Barry Smith <bsmith at mcs.anl.gov> wrote:
> One can solve Dirichilet boundary condition problems  
> using symmetric methods two ways. (1) Manually call MatGetSubMatrix()  
> to extract your symmetric subproblem or (2) call MatZeroRows() then  
> use the PCREDISTRIBUTE preconditioner to solve the reduced system using
> symmetric methods.

Both of these use extra memory and make it more difficult to reuse data
structures if boundary conditions change to or from Dirichlet as the
simulation proceeds.  It's quite easy during assembly to just not put
values into columns that you don't want values in.  This is done by
applying the same selection process to the column indices as the row
indices.  I usually do something like (using Matlab-ish notation for
subvectors)

in F(X):

  Xd = X; Xd(Dirichlet) = known                 # put the known values in the Dirichlet slots
  F = Physics(Xd)                               # evaluate physics with the correct Dirichlet values
  F(Dirichlet) = scale*(X(Dirichlet) - known)   # residual in Dirichlet slots

which implies a symmetric Jacobian (since F(All-Dirichlet) is
independent of X(Dirichlet)), so in J(X) we have:

  Xd = X; Xd(Dirichlet) = known;
  J(All-Dirichlet,All-Dirichlet) = Physics'(Xd) # don't assemble into Dirichlet rows or columns
  J(Dirichlet,Dirichlet) = scale                # put this on the diagonal

You can choose scale=1, but it is usually better (mostly with geometric
multigrid) to scale by local viscosity or whatever your analogue is.

Jed