notebook.community

Edit and run 





In [75]:



    


function hilb(n::Integer)
    a=Array(typeof(1.0),n,n)
    for i=1:n
      for j=1:n
        a[j,i]=1.0/(i+j-1.0)
      end
    end
    return a
 end

m = 1000
n = 100

 # form an ill-conditioned matrix
a = randn(m,n) * hilb(n)
x0 = ones(n, 1)
b = a*x0;



    











In [76]:



    


#Compute errors and residuals using both DGELSY and DGELSD
residy = Float64[]
residd = Float64[]
errory = Float64[]
errord = Float64[]

for i=20:60
    xy, ry = Base.LinAlg.LAPACK.gelsy!(copy(a), copy(b), 2.0^-i)
    xd, rd = Base.LinAlg.LAPACK.gelsd!(copy(a), copy(b), 2.0^-i)
    push!(residy, norm(a*xy-b))
    push!(residd, norm(a*xd-b))
    push!(errory, norm(xy-x0))
    push!(errord, norm(xd-x0))
    println(i, "  ", norm(a*xy-b), "  ", norm(a*xd-b), "  ", norm(xy-x0), "  ", norm(xd-x0))
end



    


















    






20  3.9983388475875396e-7  3.5959545568796455e-7  0.010163407424271053  0.00944839953574666
21  3.9983388475875396e-7  1.9133671399592153e-8  0.010163407424271053  0.0036254909737808124
22  1.876925887647936e-8  1.9133671399592153e-8  0.0035617401257268637  0.0036254909737808124
23  1.876925887647936e-8  1.9133671399592153e-8  0.0035617401257268637  0.0036254909737808124
24  9.86782100063115e-10  9.38237140464082e-10  0.0014015087135101916  0.0013501348916897624
25  9.86782100063115e-10  9.38237140464082e-10  0.0014015087135101916  0.0013501348916897624
26  9.86782100063115e-10  9.38237140464082e-10  0.0014015087135101916  0.0013501348916897624
27  9.86782100063115e-10  4.501216309397148e-11  0.0014015087135101916  0.00047924661963947977
28  4.5380438445153544e-11  4.501216309397148e-11  0.0004826258019643997  0.00047924661963947977
29  4.5380438445153544e-11  4.501216309397148e-11  0.0004826258019643997  0.00047924661963947977
30  4.5380438445153544e-11  1.928385746965249e-12  0.0004826258019643997  0.00017043837146270087
31  2.094728504195077e-12  1.928385746965249e-12  0.00018205479424213356  0.00017043837146270087
32  2.094728504195077e-12  1.928385746965249e-12  0.00018205479424213356  0.00017043837146270087
33  3.031389561498526e-13  2.895589073262959e-13  5.964779758724967e-5  5.97189915062968e-5
34  3.031389561498526e-13  2.895589073262959e-13  5.964779758724967e-5  5.97189915062968e-5
35  3.031389561498526e-13  2.895589073262959e-13  5.964779758724967e-5  5.97189915062968e-5
36  2.7314827139603245e-13  2.988599326975199e-13  2.6196956059889193e-5  2.086033427559339e-5
37  2.7314827139603245e-13  2.988599326975199e-13  2.6196956059889193e-5  2.086033427559339e-5
38  2.7314827139603245e-13  2.988599326975199e-13  2.6196956059889193e-5  2.086033427559339e-5
39  2.7314827139603245e-13  3.673611658204464e-13  2.6196956059889193e-5  5.860160544963344e-5
40  2.67669791325659e-13  3.673611658204464e-13  1.9185107605319417e-5  5.860160544963344e-5
41  2.67669791325659e-13  3.673611658204464e-13  1.9185107605319417e-5  5.860160544963344e-5
42  2.67669791325659e-13  3.0028124508327845e-13  1.9185107605319417e-5  0.00043448069356473
43  2.77042069976064e-13  3.0028124508327845e-13  0.0019624885428488927  0.00043448069356473
44  2.77042069976064e-13  3.0028124508327845e-13  0.0019624885428488927  0.00043448069356473
45  2.631749703260356e-13  2.9705648594906537e-13  0.004086375683080754  0.0042076175429510455
46  2.631749703260356e-13  2.9705648594906537e-13  0.004086375683080754  0.0042076175429510455
47  2.631749703260356e-13  2.9705648594906537e-13  0.004086375683080754  0.0042076175429510455
48  2.631749703260356e-13  2.9705648594906537e-13  0.004086375683080754  0.0042076175429510455
49  2.631749703260356e-13  3.624092461081448e-13  0.004086375683080754  0.09187796069656362
50  3.0694417491259723e-13  3.624092461081448e-13  0.07709801273810697  0.09187796069656362
51  3.0694417491259723e-13  3.624092461081448e-13  0.07709801273810697  0.09187796069656362
52  3.0694417491259723e-13  2.861422398112231e-13  0.07709801273810697  0.2464050339073647
53  2.7117135693396676e-13  2.861422398112231e-13  2.940516677215117  0.2464050339073647
54  3.1456702835423373e-13  2.795849925388482e-13  64.00113443968891  29.59732691857252
55  3.0456647763262273e-13  2.795849925388482e-13  73.755046216312  29.59732691857252
56  3.375003860296157e-13  2.795849925388482e-13  93.91194090001832  29.59732691857252
57  3.0964726316673704e-13  2.795849925388482e-13  97.76908847864615  29.59732691857252
58  3.0964726316673704e-13  2.795849925388482e-13  97.76908847864615  29.59732691857252
59  3.0964726316673704e-13  2.795849925388482e-13  97.76908847864615  29.59732691857252
60  3.0964726316673704e-13  2.795849925388482e-13  97.76908847864615  29.59732691857252

















In [77]:



    


using Gadfly
plot(layer(x=map(x->2.0^-x, 20:60), y=residy, Geom.line, Theme(default_color=color("red"))),
     layer(x=map(x->2.0^-x, 20:60), y=errory, Geom.line, Theme(default_color=color("red"))),
     layer(x=map(x->2.0^-x, 20:60), y=residd, Geom.line),
     layer(x=map(x->2.0^-x, 20:60), y=errord, Geom.line),
     layer(xintercept=[eps(), eps()*(1+1.8*√n*0+1)*1.9*√n*(n+1), sqrt(eps())], Geom.vline(color=color("orange"))),
     Guide.xlabel("Tolerance"), Guide.ylabel("Error or residual"), Scale.x_log2, Scale.y_log2,
)



    


















    
Out[77]:









  
    Tolerance
  
  
    2-110
    2-100
    2-90
    2-80
    2-70
    2-60
    2-50
    2-40
    2-30
    2-20
    2-10
    20
    210
    220
    230
    2-100
    2-98
    2-96
    2-94
    2-92
    2-90
    2-88
    2-86
    2-84
    2-82
    2-80
    2-78
    2-76
    2-74
    2-72
    2-70
    2-68
    2-66
    2-64
    2-62
    2-60
    2-58
    2-56
    2-54
    2-52
    2-50
    2-48
    2-46
    2-44
    2-42
    2-40
    2-38
    2-36
    2-34
    2-32
    2-30
    2-28
    2-26
    2-24
    2-22
    2-20
    2-18
    2-16
    2-14
    2-12
    2-10
    2-8
    2-6
    2-4
    2-2
    20
    22
    24
    26
    28
    210
    212
    214
    216
    218
    220
    2-100
    2-50
    20
    250
    2-100
    2-95
    2-90
    2-85
    2-80
    2-75
    2-70
    2-65
    2-60
    2-55
    2-50
    2-45
    2-40
    2-35
    2-30
    2-25
    2-20
    2-15
    2-10
    2-5
    20
    25
    210
    215
    220
  
  
    Error or residual
  
  
    
      
    
    
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
    
    
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
    
    
      
        
        
        
      
      
        
      
      
        
      
      
        
      
      
        
      
    
    
      
        
        
      
      
        
      
      
        
      
      
        
        
      
    
  
  
    2-120
    2-110
    2-100
    2-90
    2-80
    2-70
    2-60
    2-50
    2-40
    2-30
    2-20
    2-10
    20
    210
    220
    230
    240
    250
    260
    270
    280
    2-110
    2-108
    2-106
    2-104
    2-102
    2-100
    2-98
    2-96
    2-94
    2-92
    2-90
    2-88
    2-86
    2-84
    2-82
    2-80
    2-78
    2-76
    2-74
    2-72
    2-70
    2-68
    2-66
    2-64
    2-62
    2-60
    2-58
    2-56
    2-54
    2-52
    2-50
    2-48
    2-46
    2-44
    2-42
    2-40
    2-38
    2-36
    2-34
    2-32
    2-30
    2-28
    2-26
    2-24
    2-22
    2-20
    2-18
    2-16
    2-14
    2-12
    2-10
    2-8
    2-6
    2-4
    2-2
    20
    22
    24
    26
    28
    210
    212
    214
    216
    218
    220
    222
    224
    226
    228
    230
    232
    234
    236
    238
    240
    242
    244
    246
    248
    250
    252
    254
    256
    258
    260
    262
    264
    266
    268
    270
    2-200
    2-100
    20
    2100
    2-110
    2-105
    2-100
    2-95
    2-90
    2-85
    2-80
    2-75
    2-70
    2-65
    2-60
    2-55
    2-50
    2-45
    2-40
    2-35
    2-30
    2-25
    2-20
    2-15
    2-10
    2-5
    20
    25
    210
    215
    220
    225
    230
    235
    240
    245
    250
    255
    260
    265
    270
  



  





















In [78]:

Content source: jiahao/jin

Similar notebooks:

notebook.community | gallery | about