This is a nasty bug. The code:
import mpmath as mpm
N=64
A=mpm.matrix(mpm.linspace(-N//2,N//2,N))+1J*mpm.ones(N,1)
b=mpm.ones(N,1)
x,res=mpm.qr_solve(A,b)
x_correct = (A.H*b)/(A.H*A)[0]
print(x[0], '\t', mpm.norm(mpm.residual(A,x,b)))
print(x_correct[0], '\t', mpm.norm(mpm.residual(A,x_correct,b)))
should print two identical lines, but prints this:
(-3.50817994558216e-5 + 0.00319170755575594j) 8.03982720027231
(-5.0297598862402e-18 - 0.00283150309367931j) 7.98866595884473
Also,
print(res,'\n',mpm.norm(mpm.residual(A,x,b)))
should print two equal values, but prints
8.02818594205587
8.03982720027231
lu_solve works correctly:
print(mpm.lu_solve(A,b), '\n', x_correct)
prints
[(-5.0297598862402e-18 - 0.00283150309367931j)]
[(-5.0297598862402e-18 - 0.00283150309367931j)]
This is a nasty bug. The code:
should print two identical lines, but prints this:
Also,
should print two equal values, but prints
lu_solve works correctly:
prints