Hi Karli,
Karl Rupp <[email protected]> writes:
> The 'standard' criterion is indeed to monitor the residual norm
> ||b - Ax||
> and quit if this becomes several orders of magnitude smaller than the
> initial residual. Since the initial guess is typically the zero vector,
> the initial residual is easily computed as ||b||. Common values for the
> tolerance are in the range 10^{-6}, which should be okay even with
> single precision. To make your system matrix well-behaved, a good hint
> is to make your random matrix an M-matrix, i.e. you select all
> off-diagonal entries negative, while the diagonal entries are equal or
> larger than the negative sum of the off-diagonal entries. This way the
> condition number remains reasonable, so you also get useful
> approximations of the true solution. Anyway, if you just monitor the
> residual, you should be fine.
Yeah, constructing an M-matrix is a good idea, and helps with the
conditioning. But the preconditioners still seem to suffer induced
instability. Consider the following test output:
===================================== FAILURES
======================================
__________ test_gmres_with_ilu0_solve_compressed_matrix_A_vector_b_float32
__________
def _test():
solver_tag = solver_tag_type(tolerance=tol/10)
precond_tag = precond_tag_type()
vcl_system = sparse_type.generate_fdm_laplace(points_x_y, points_x_y,
dtype=dt)
#vcl_system = get_sparse_matrix(20, dtype=dt, sparse_type=sparse_type)
numpy_system = vcl_system.as_ndarray() # TODO: SciPy-ise
numpy_solution, vcl_solution = vector_getter(vcl_system.size1, dt,
vector=np.ones(vcl_system.size1).astype(dt))
numpy_rhs, vcl_rhs = vector_getter(vcl_system.size1, dt,
vector=vcl_system*vcl_solution)
# solve using pyviennacl
vcl_solution = p.solve(vcl_system, vcl_rhs, solver_tag, precond_tag)
# compare with numpy_solution
act_diff = math.fabs(diff(vcl_rhs, vcl_system.dot(vcl_solution)))
> assert act_diff <= tol, "diff was {} > tolerance {}".format(act_diff,
> tol)
E AssertionError: diff was 78631665664.0 > tolerance 0.001
Here, you can see that I'm using the generate_fdm_laplace routine
(points_x = points_y = points_x_y = 5), which I trust to give a
sufficiently well-conditioned symmetric positive definite matrix, though
I haven't checked the condition number. The LU decomposition exists, so
why should using ILU0 go so wrong?
If I use my matrix-generation routine, which now produces a 20x20
symmetric M-matrix with 52 nonzeros (eight nonzeros per row and 20 on
the diagonal, to give an approximate sparsity of 0.01) and diagonals
about twice the negative sum of the off-diagonals (which are all equal
to -1.0 where non-zero), or 1.0 (if there are no nonzero off-diagonals),
then the test usually passes; though it sometimes still fails with a
small error.
Note that in both cases I'm using single precision with test tolerance
0.001 and solver tolerance 0.0001. If I use a test tolerance of 0.01 and
solver tolerance of 0.001, then the test passes using the Laplace
matrix. I still get some failures on other iterative tests using my
random matrices, though.
Oh, and I'm now doing as you and Philippe suggest, and using ||Ax - b||
as my test criterion.
Thanks for your help; it's useful to talk with someone with more
expertise in the relevant theory!
Toby
>> I tried to check for what the ViennaCL test does, but I couldn't find
>> one for the iterative solvers!
>
> Yeah, these are currently tested in other Vienna* packages, but
> (unfortunately) there is no strict test in the Nightly test suite, only
> the solverbench.
>
> Best regards,
> Karli
>
> ------------------------------------------------------------------------------
> Infragistics Professional
> Build stunning WinForms apps today!
> Reboot your WinForms applications with our WinForms controls.
> Build a bridge from your legacy apps to the future.
> http://pubads.g.doubleclick.net/gampad/clk?id=153845071&iu=/4140/ostg.clktrk
--
Toby St Clere Smithe
http://tsmithe.net
------------------------------------------------------------------------------
Infragistics Professional
Build stunning WinForms apps today!
Reboot your WinForms applications with our WinForms controls.
Build a bridge from your legacy apps to the future.
http://pubads.g.doubleclick.net/gampad/clk?id=153845071&iu=/4140/ostg.clktrk
_______________________________________________
ViennaCL-devel mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/viennacl-devel
