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.
>>
>>
