fixes

10-12/lesson.ipynb

@@ -1282,7 +1282,7 @@
     "include(\"DIS.jl\")\n",
     "include(\"UNM.jl\")\n",
     "\n",
-    "plt = plot(yscale = :log,\n",
+    "plt = plot(yscale = :log10,\n",
     "                   xlims=(0, 35),\n",
     "                   ylims=(1e-15, Inf),\n",
     "                   guidefontsize=16)\n",
@@ -1299,7 +1299,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Julia 1.9.3",
+   "display_name": "Julia 1.9.4",
    "language": "julia",
    "name": "julia-1.9"
   },
@@ -1307,7 +1307,7 @@
    "file_extension": ".jl",
    "mimetype": "application/julia",
    "name": "julia",
-   "version": "1.9.3"
+   "version": "1.9.4"
   }
  },
  "nbformat": 4,

11-09/NWTN.jl

@@ -12,7 +12,7 @@ function NWTN(f;
         MInf::Real=-Inf,
         mina::Real=1e-16,
         plt::Union{Plots.Plot, Nothing}=nothing,
-        plotatend::Bool=true,
+        plotatend::Bool=false,
         Plotf::Integer=0,
         printing::Bool=true)::Tuple{AbstractArray, String}
 
@@ -99,7 +99,7 @@ function NWTN(f;
         am = 0
         a = as
         phipm = phip0
-        while (feval ≤ MaxFeval ) && ((as - am)) > mina && (phips > 1e-12)
+        while (feval ≤ MaxFeval ) && (as - am) > mina && (phips > 1e-12)
         
             # compute the new value by safeguarded quadratic interpolation
             a = (am * phips - as * phipm) / (phips - phipm)
@@ -107,7 +107,7 @@ function NWTN(f;
         
             # compute phi(a)
             phia, phip = f2phi(a, true)
-            
+
             if (phia ≤ phi0 + m1 * a * phip0) && (abs(phip) ≤ -m2 * phip0)
                 break  # Armijo + strong Wolfe satisfied, we are done
             end
@@ -372,12 +372,12 @@ function NWTN(f;
         end
     end
 
+    if Plotf ≥ 2
+        plot!(plt, gap)
+    elseif Plotf == 1 && n == 2
+        plot!(plt, PXY[1, :], PXY[2, :])
+    end
     if plotatend
-        if Plotf ≥ 2
-            plot!(plt, gap)
-        elseif Plotf == 1 && n == 2
-            plot!(plt, PXY[1, :], PXY[2, :])
-        end
         display(plt)
     end
 

project/thinQR.jl

@@ -1,7 +1,7 @@
 module thinQR
 
 import Base: size, show, getproperty, getfield, propertynames, \, *
-using LinearAlgebra: norm, I, triu, diagm, ldiv!
+using LinearAlgebra: norm, I, triu, diagm, ldiv!, dot
 
 export QRthin, qrhous!, qyhous, qyhoust
 
@@ -18,16 +18,22 @@ function qrhous!(A::AbstractMatrix{T})::QRthin{T} where T
     m, n = size(A)
 
     d = zeros(n)
-    for j ∈ 1:n
+    @inbounds begin
-        s = norm(A[j:m, j])
+        for j ∈ 1:n
-        # iszero(s) && throw(ArgumentError("The matrix A is singular"))
+            s = norm(A[j:m, j])
-        d[j]= copysign(s, -A[j,j])
+            # iszero(s) && throw(ArgumentError("The matrix A is singular"))
-        fak = sqrt(s * (s + abs(A[j,j])))
+            d[j]= copysign(s, -A[j,j])
-        A[j,j] -= d[j]
+            fak = sqrt(s * (s + abs(A[j,j])))
-        A[j:m, j] ./= fak
+            A[j,j] -= d[j]
+            
+            A[j:m, j] ./= fak
 
-        if j < n
+            if j < n
-            A[j:m, j+1:n] -= A[j:m, j] * (A[j:m, j]' * A[j:m, j+1:n])
+                A[j:m, j+1:n] -= A[j:m, j] * (A[j:m, j]' * A[j:m, j+1:n])
+                #for k ∈ j+1:n
+                #    A[j:m, k] -= A[j:m, j] .* dot(A[j:m, j], A[j:m, k])
+                #end
+            end
         end
     end
 
@@ -47,7 +53,8 @@ function qyhoust(A::QRthin{T}, y::AbstractArray{T}) where T
     m, n = size(A.A)
     z = deepcopy(y)
     for j ∈ 1:n
-        z[j:m] = z[j:m] - A.A[j:m, j] * (A.A[j:m, j]' * z[j:m])
+        # z[j:m] = z[j:m] - A.A[j:m, j] * (A.A[j:m, j]' * z[j:m])
+        z[j:m] -= A.A[j:m, j] .* dot(A.A[j:m, j], z[j:m])
     end
     return z
 end
@@ -101,7 +108,7 @@ function (\)(A::QRthin{T}, b::AbstractVector{T}) where T
     v = qyhoust(A, b)
     x = zeros(m)
     for j ∈ m:-1:1
-        x[j] = (v[j] - x[j+1:m]' * A.A[j, j+1:m]) * A.d[j]^-1
+        x[j] = (v[j] - dot(x[j+1:m], A.A[j, j+1:m])) * A.d[j]^-1
     end
     return x
 end