Suggestion: merge change_smoothers() and setup_smoothers()

[Ivan]

Hi guys!

I'm considering PyAMG for a university project and I'm currently
exploring its functionality.

I found that pyamg.relaxation.smoothing there are two functions with
quite similar functionality, which in my opinion can be merged.

1. As far as I understand, setup_smoothers() can be implemented with a
call to change_smoothers() (functionality of the former is a strict
subset of the latter)

2. I see no reason for these both functions to exist together. I think
that it would be better to delete setup_smoothers, (and probably
rename change_smoothers to the more appropriate name of
setup_smoothers).

What do you think?
Ivan

[Nathan Bell]

Hi Ivan,

Good suggestion. I've filed a report [1] to track this issue. We'll
try to resolve this in v1.1, which we hope to release soon.

[1] http://code.google.com/p/pyamg/issues/detail?id=108

--
Nathan Bell wnb...@gmail.com
http://www.wnbell.com/

[Ivan Tonchev]

On Sep 27, 9:51 am, Nathan Bell <wnb...@gmail.com> wrote:

OK, great!

Btw do you have developers list for pyamg? And do you accept user
contributions? I'd like to extend pyamg gallery with some non-stencil
matrices, write a note about how to use pyamg with custom aggregation
schemes/custom prolongation/restriction matrices, probably add some
aggregations/smoothers. Would you accept such contributions and how
would I go about it?

Btw I also noticed some things that could be improved:
I. In the CG routine:
1. stopping criterion should be in preconditioner energy norm (when a
preconditioner is available), or if you don't like this, at least rz
== 0 should be added to the stopping criterion (when a preconditioner
is available), since after rz == 0, PCG is undefined
2. there must be some sort of checking that rz > 0. PCG is great iff
the preconditioner is positive definite and you should warn the user
if you detect it isn't! Also, floating point roundoff sometimes
introduces negative numerical zeros (e.g. rz/rz0 = -1e-27), which
should be handled correctly.
3. It's not a good idea to always compute next residual as r -= alpha
* Ap. After several iterations, residual computation should be
restarted with r = b - A*x

II. condition number estimation routines (pyamg.util.linalg.condest/
_approximate_eigenvalues) could handle preconditioned matrices also.

If you have a developers list probably we can discuss dev matters
there?

Ivan

[Luke Olson]

> --
> You received this message because you are subscribed to the Google Groups "pyamg-user" group.
> To post to this group, send email to pyamg...@googlegroups.com.
> To unsubscribe from this group, send email to pyamg-user+...@googlegroups.com.
> For more options, visit this group at http://groups.google.com/group/pyamg-user?hl=en.
>

Ivan,

This would be great. I'll send you a separate message. Examples in particular are in high demand.

As for CG, you're correct, rz should be checked. I always found Matlab's version strange in that they check == instead of within some tolerance. I'll add a ticket to make sure the stopping criteria handles these cases. Also, the recomputation of r=b-Ax is a good idea. I'm guessing that all of our krylov methods could be a little more robust in this sense. I'm working on this now in fact as I'm ripping out a separate iterative solvers/preconditioners package.

For the condition number estimates, these were quick ones so there's certainly room for expansions. If we can get some more robust estimates and general routines in there, that would be great. This part is probably secondary though.

I think your suggestions here are very appropriate for this list. There is not a developers list specifically, but if there are additional details you could always email me separately and we'll go from there.

Luke

[Ivan Tonchev]

> > For more options, visit this group athttp://groups.google.com/group/pyamg-user?hl=en.

>
> Ivan,
>
> This would be great.  I'll send you a separate message.  Examples in particular are in high demand.
>
> As for CG, you're correct, rz should be checked.  I always found Matlab's version strange in that they check == instead of within some tolerance.  I'll add a ticket to make sure the stopping criteria handles these cases.  Also, the recomputation of r=b-Ax is a good idea.  I'm guessing that all of our krylov methods could be a little more robust in this sense.  I'm working on this now in fact as I'm ripping out a separate iterative solvers/preconditioners package.
>

Hi Luke,

If you are rewriting the iterative solvers package, I'd suggest to add
one more method: (preconditined) simple iteration. This is *the* test
for preconditioner optimality. Scientists/students (like me) who want
to test their preconditioner quality could potentially find great use
of it, since unlike Krylov methods, SI takes forever to converge even
if just one of the (preconditioned) matrix eigenvalues is close to
zero (or two).

What I mean by simple iteration is:
B*(xk+1 - xk) = b - Axk,
where B^-1 is the preconditioner (in PyAMG notation, B^-1 = M). Stop
criterion is the B^-1 norm of rk=b - Axk to be <= tol * B^-1 norm of
the right hand side b. This of course is assuming B is symmetric and
positive definite.

Actually, speaking about stop criteria, another useful addition would
be to add an 'exact' parameter to CG/Simple iteration (any solver)
routines. When the user knows what is the exact solution of the system
(and he is just testing the quality of his preconditioner), he could
supply exact != None, so the solver could use it to compute The True
error in whatever norm is appropriate (e.g. A-norm), and not some
approximation to it via a norm of the residual (e.g. the B^-1-norm).

Here is my code for CG and SI: http://pastie.org/1194773 . It is not
exactly clean, but works well for my tests and illustrates my ideas.
The part where I handle rs < 0 has to be rethought, but you could just
ignore it. Also, pardon my asserts (so far I've been the only user of
my code so I didn't care enough to write exceptions).

And another great thing to add to your solvers would be a verbosity
parameter. When the matrix is large (several millions of unknowns),
iterations could take seconds and solver routines minutes. It would be
nice if the user has some easy way to see that something *is actually
happening* (callbacks are not appropriate for such a simple usecase,
besides, it's not clear what parameters to pass to the callback). For
example, in my code, if verb==True, on each iteration I print current
residual, relative residual, true error norm (if it's available) and
reduction factor.

> For the condition number estimates, these were quick ones so there's certainly room for expansions.  If we can get some more robust estimates and general routines in there, that would be great.  This part is probably secondary though.
>

Yes, this is definitely secondary, but here is my Lanczos code:
http://pastie.org/1194784 . It is not numerically robust (I've
implemented just the straightforward Lanczos algorithm, not the
restarted one), but at least it handles the preconditioned case. Of
course it's only for symmetric, real matrices (for now my focus is on
the symmetric case).

Btw, here's how I check the symmetry and positive definiteness of my
matrices: http://pastie.org/1194791 . Similar routines could be useful
to PyAMG...

> I think your suggestions here are very appropriate for this list.  There is not a developers list specifically, but if there are additional details you could always email me separately and we'll go from there.

I'm glad I can help. Btw, I'd prefer to send patches for stuff similar
to what I've just pointed out, but the code in google SVN repository
looks old to me... What would be reasonable in the future -- should I
generate patches against it or you have a repository with more recent
code?

>
> Luke

Ivan