The OP wants extremely high precision and indicated that he was willing to
factor the matrix. I recommended iterative refinement which converges very
quickly, and exploits the state-of-the-art direct solvers. The solvers in
IterativeSolvers.jl are for a different domain, where the matrix is too
large or expensive to factor. To get high accuracy with them generally
requires tailored preconditioners, which are not "out of the box". In fact
one usually need a preconditioner to get any convergence with the
non-symmetric ones for interesting ranks. (I've been struggling for months
to find a good preconditioner for an application of GMRES in my work, so
this is a sore point.)
On Wednesday, August 10, 2016 at 10:56:22 PM UTC-4, Chris Rackauckas wrote:
>
> Yes textbook answer is, why do you want to use `\`? Iterative techniques
> are likely better suited for the problem. There's no need to roll you own,
> the package IterativeSolvers.jl has a good number of techniques implemented
> which are well-suited for the problem since A is a large sparse matrix.
> Their methods should work out of the box with Bigs, though you will likely
> want to adjust the tolerances.
>
> On Wednesday, August 10, 2016 at 7:37:35 PM UTC-7, Ralph Smith wrote:
>>
>> Here is a textbook answer. Appropriate choice of n depends on condition
>> of A.
>>
>> """
>>
>> iterimprove(A,b,n=1,verbose=true)
>>
>>
>>
>> Solve `A x = b` for `x` using iterative improvement
>>
>> """
>>
>> function iterimprove{T<:AbstractFloat}(A::SparseMatrixCSC{T},
>>> b::Vector{T},n=1,verbose=true)
>>>
>> eps(T) < eps(Float64) || throw(ArgumentError("wrong
>>> implementation"))
>>
>> A0 = SparseMatrixCSC{Float64}(A)
>>> F = factorize(A0)
>>> x = zeros(b)
>>> r = copy(b)
>>> for iter = 1:n+1
>>> y = F \ Vector{Float64}(r)
>>> for i in eachindex(x)
>>> x[i] += y[i]
>>> end
>>> r = b - A * x
>>> if verbose
>>> @printf "at iter %d resnorm = %.3g\n" iter norm(r)
>>> end
>>> end
>>> x
>>> end
>>
>>
>>
>> On Wednesday, August 10, 2016 at 3:47:10 PM UTC-4, Nicklas Andersen wrote:
>>>
>>> Hello
>>>
>>> I'm trying to solve a large, sparse and unsymmetrical linear system Ax =
>>> b.
>>> For this task I'm using Julias *SparseMatrixCSC *type for the
>>> definition of my matrices and Julias built in backslash ' \ ' operator for
>>> the solution of the system.
>>> I need *quadruple precision* and thus I've been trying to implement my
>>> routine with the *BigFloat *type together with the SparseMatrixCSC type.
>>>
>>> To illustrate this, I give a simple example here:
>>> set_bigfloat_precision(128);
>>> A = speye(BigFloat, 2, 2);
>>> b = ones(BigFloat, 2, 1);
>>> x = A\b;
>>>
>>> If I do this I either get a StackOverFlow error:
>>> ERROR: StackOverflowError:
>>> in copy at array.jl:100
>>> in float at sparse/sparsematrix.jl:234
>>> in call at essentials.jl:57 (repeats 254 times)
>>>
>>> or the solver seems to run forever and never terminates. As the second
>>> error indicates it seems like the sparse solver only accepts the normal
>>> *float* types.
>>> My question is then, is there a way to get quadruple precision with the
>>> standard solvers in Julia, in this case UMFpack I assume ? or should I look
>>> for something else (in this case any suggestions :) ) ?
>>>
>>> Regards Nicklas A.
>>>
>>>
