Forming a sparse approximation of a MatShell

[Barry Smith]

    Jed,

     The Mat coloring code can also be used directly, not through  
SNES. Once you have
the coloring for the matrix (you can get that with MatGetColoring(),  
of course, this assumes
you have already set a nonzero pattern for your matrix)). Call  
MatFDColoringCreate()
then MatFDColoringSetFunction(), MatFDColoringSetFromOptions() and  
then MatFDColoringApply().

    Good luck,

     Barry

On Apr 9, 2008, at 2:13 PM, Jed Brown wrote:
> I'm trying to improve the preconditioning of my spectral collocation  
> method for
> non-Newtonian incompressible Stokes flow.  My current algorithm uses  
> MatShell
> for the full Jacobian as well as each of its blocks [A B1'; B2 0]  
> and the Schur
> complement S = -B2*A*B1'.  I needed a preconditioner for A so I  
> thought I'd
> solve the same problem using finite differences on the Chebyshev  
> nodes.  In
> reality, the stencil is really ugly in 3D so I just used a simpler  
> elliptic
> operator.  This works okay, but it's performance decays  
> significantly as I
> increase the continuation parameter.  Also, dealing with general  
> boundary
> conditions is rather tricky and it seems to be a much weaker  
> preconditioner when
> I have mixed boundary conditions.  To rectify this, I tried a finite  
> element
> discretization on the Chebyshev nodes (using Q1 elements).  This  
> must be scaled
> by the inverse (lumped) mass matrix due to the collocation nature of  
> the
> spectral method.  Strangely, even though it captures all the terms  
> in the
> Jacobian, it is slightly weaker than the finite difference version.   
> At least it
> is less error-prone and boundary conditions are easier to get right.
> Regardless, forming the explicit matrix separately from the spectral  
> matrix
> causes a duplication of concepts that have to be kept in sync.  So I  
> started
> thinking, the spectral matrix is pretty cheap to apply a few times,  
> so perhaps I
> can use a coloring to compute a sparse approximation.  However, the
> documentation I found is using the function from the SNES context to  
> form the
> matrix.  In my case, the entire Jacobian doesn't help, I just want an
> approximation of A.  (A itself is full, but implemented via FFT.)   
> What is the
> correct way to do this?  Should I just stick with finite differences  
> or finite
> elements?
>
> Also, any ideas for preconditioning S?  It's condition number also  
> grows
> significantly with the continuation parameter.
>
> Thanks,
>
> Jed