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