Hey,
>> 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?
even though ILU0 is somewhat random and does not have a strong
mathematical theorey, the failure you report is certainly unexpected.
I'll look into some time next week.
> 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.
You might just need to run with a tighter tolerance. If the residual of
Ax=b is below a tolerance eps, then you can only expect the solution
vector x to be below eps * condition_number(A).
> 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.
Iterative solvers should better be run in double precision, single
precision is too often insufficient... So, if the test system has double
precision available, let me suggest to skip the tests in single
precision to save time and false positives. If double precision is not
available, use a smaller solver tolerance (or a larger test tolerance).
Best regards,
Karli
------------------------------------------------------------------------------
Want fast and easy access to all the code in your enterprise? Index and
search up to 200,000 lines of code with a free copy of Black Duck
Code Sight - the same software that powers the world's largest code
search on Ohloh, the Black Duck Open Hub! Try it now.
http://p.sf.net/sfu/bds
_______________________________________________
ViennaCL-devel mailing list
ViennaCL-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/viennacl-devel
