From 4e5eb819488a23af86270bad0c2eb722aebe37ae Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Sun, 21 Jul 2024 12:56:16 +0200
Subject: [PATCH 01/19] add qless tutorial

---
 docs/src/tutorials/qless.md | 58 +++++++++++++++++++++++++++++++++++++
 1 file changed, 58 insertions(+)
 create mode 100644 docs/src/tutorials/qless.md

diff --git a/docs/src/tutorials/qless.md b/docs/src/tutorials/qless.md
new file mode 100644
index 0000000..08f6404
--- /dev/null
+++ b/docs/src/tutorials/qless.md
@@ -0,0 +1,58 @@
+```@example qless1
+using QRMumps, LinearAlgebra, SparseArrays, Printf
+
+irn = [1, 3, 5, 2, 3, 5, 1, 4, 4, 5, 2, 1, 3]
+jcn = [1, 1, 1, 2, 3, 3, 4, 4, 5, 5, 6, 7, 7]
+val = [1.0, 2.0, 3.0, 1.0, 1.0, 2.0, 4.0, 1.0, 5.0, 1.0, 3.0, 6.0, 1.0]
+
+A = sparse(irn, jcn, val, 5, 7)
+b = [40.0, 10.0, 44.0, 98.0, 87.0]
+y_star = [0.0, 1.0, 2.0, 3.0, 4.0]
+y₁ = similar(b,7)
+y₂ = similar(b)
+
+qrm_init()
+
+spmat = qrm_spmat_init(A)
+spfct = qrm_spfct_init(spmat)
+qrm_analyse!(spmat, spfct; transp='t')
+qrm_set(spfct, "qrm_keeph", 0)
+qrm_factorize!(spmat, spfct, transp='t')
+
+qrm_solve!(spfct, b, y₁; transp='t')
+qrm_solve!(spfct, y₁, y₂; transp='n')
+
+error_norm = norm(y₂ - y_star)
+
+@printf("Error norm ‖y* - y‖ = %10.5e\n", error_norm)
+```
+
+```@example qless2
+using QRMumps, LinearAlgebra, SparseArrays, Printf
+
+irn = [1, 1, 1, 2, 3, 3, 4, 4, 5, 5, 6, 7, 7]
+jcn = [1, 3, 5, 2, 3, 5, 1, 4, 4, 5, 2, 1, 3]
+val = [1.0, 2.0, 3.0, 1.0, 1.0, 2.0, 4.0, 1.0, 3.0, 1.0, 3.0, 2.0, 1.0]
+
+A = sparse(irn, jcn, val, 7, 5)
+b = [22.0, 2.0, 13.0, 8.0, 17.0, 6.0, 5.0]
+x_star = [1.0, 2.0, 3.0, 4.0, 5.0]
+
+qrm_init()
+
+spmat = qrm_spmat_init(A)
+spfct = qrm_spfct_init(spmat)
+
+qrm_set(spfct, "qrm_keeph", 0)
+
+qrm_analyse!(spmat, spfct)
+qrm_factorize!(spmat, spfct)
+
+z = A'*b
+x₁ = qrm_solve(spfct, z; transp = 't')
+x₂ = qrm_solve(spfct, x₁; transp = 'n')
+
+error_norm = norm(x₂ - x_star)
+
+@printf("Error norm ‖x* - x‖ = %10.5e\n", error_norm)
+```
\ No newline at end of file

From a98097c0c60e12dcbf6689bdee1d42448acfefa4 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Mon, 22 Jul 2024 10:22:07 +0200
Subject: [PATCH 02/19] Apply suggestions

---
 docs/src/tutorials/qless.md | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/docs/src/tutorials/qless.md b/docs/src/tutorials/qless.md
index 08f6404..7c9d5e6 100644
--- a/docs/src/tutorials/qless.md
+++ b/docs/src/tutorials/qless.md
@@ -1,22 +1,24 @@
 ```@example qless1
 using QRMumps, LinearAlgebra, SparseArrays, Printf
 
+m, n = 5, 7
 irn = [1, 3, 5, 2, 3, 5, 1, 4, 4, 5, 2, 1, 3]
 jcn = [1, 1, 1, 2, 3, 3, 4, 4, 5, 5, 6, 7, 7]
 val = [1.0, 2.0, 3.0, 1.0, 1.0, 2.0, 4.0, 1.0, 5.0, 1.0, 3.0, 6.0, 1.0]
 
-A = sparse(irn, jcn, val, 5, 7)
+A = sparse(irn, jcn, val, m, n)
 b = [40.0, 10.0, 44.0, 98.0, 87.0]
 y_star = [0.0, 1.0, 2.0, 3.0, 4.0]
-y₁ = similar(b,7)
-y₂ = similar(b)
+y₁ = zeros(n)
+y₂ = zeros(m)
 
 qrm_init()
 
 spmat = qrm_spmat_init(A)
 spfct = qrm_spfct_init(spmat)
-qrm_analyse!(spmat, spfct; transp='t')
+
 qrm_set(spfct, "qrm_keeph", 0)
+qrm_analyse!(spmat, spfct; transp='t')
 qrm_factorize!(spmat, spfct, transp='t')
 
 qrm_solve!(spfct, b, y₁; transp='t')
@@ -30,11 +32,12 @@ error_norm = norm(y₂ - y_star)
 ```@example qless2
 using QRMumps, LinearAlgebra, SparseArrays, Printf
 
+m, n = 7, 5
 irn = [1, 1, 1, 2, 3, 3, 4, 4, 5, 5, 6, 7, 7]
 jcn = [1, 3, 5, 2, 3, 5, 1, 4, 4, 5, 2, 1, 3]
 val = [1.0, 2.0, 3.0, 1.0, 1.0, 2.0, 4.0, 1.0, 3.0, 1.0, 3.0, 2.0, 1.0]
 
-A = sparse(irn, jcn, val, 7, 5)
+A = sparse(irn, jcn, val, m, n)
 b = [22.0, 2.0, 13.0, 8.0, 17.0, 6.0, 5.0]
 x_star = [1.0, 2.0, 3.0, 4.0, 5.0]
 

From 70e8e0cabdf92e44d4f120d7e544cbe4cb841f46 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Mon, 22 Jul 2024 10:37:27 +0200
Subject: [PATCH 03/19] Add comments to tutorial

---
 docs/src/tutorials/qless.md | 20 ++++++++++++++++++++
 1 file changed, 20 insertions(+)

diff --git a/docs/src/tutorials/qless.md b/docs/src/tutorials/qless.md
index 7c9d5e6..1d4ea3f 100644
--- a/docs/src/tutorials/qless.md
+++ b/docs/src/tutorials/qless.md
@@ -1,6 +1,7 @@
 ```@example qless1
 using QRMumps, LinearAlgebra, SparseArrays, Printf
 
+# Initialize data
 m, n = 5, 7
 irn = [1, 3, 5, 2, 3, 5, 1, 4, 4, 5, 2, 1, 3]
 jcn = [1, 1, 1, 2, 3, 3, 4, 4, 5, 5, 6, 7, 7]
@@ -12,18 +13,27 @@ y_star = [0.0, 1.0, 2.0, 3.0, 4.0]
 y₁ = zeros(n)
 y₂ = zeros(m)
 
+# Initialize QRMumps
 qrm_init()
 
+# Initialize data structures
 spmat = qrm_spmat_init(A)
 spfct = qrm_spfct_init(spmat)
 
+# Specify not storing Q for the QR factorization
 qrm_set(spfct, "qrm_keeph", 0)
+
+# Perform symbolic analysis of Aᵀ and factorize Aᵀ = QR
 qrm_analyse!(spmat, spfct; transp='t')
 qrm_factorize!(spmat, spfct, transp='t')
 
+# Solve Rᵀy₁ =  b  
 qrm_solve!(spfct, b, y₁; transp='t')
+
+# Solve Ry₂ = y₁
 qrm_solve!(spfct, y₁, y₂; transp='n')
 
+# Overall, RᵀRy₂ = b. Equivalently, RᵀQᵀQRy₂ = b or AAᵀy₂ = b
 error_norm = norm(y₂ - y_star)
 
 @printf("Error norm ‖y* - y‖ = %10.5e\n", error_norm)
@@ -32,6 +42,7 @@ error_norm = norm(y₂ - y_star)
 ```@example qless2
 using QRMumps, LinearAlgebra, SparseArrays, Printf
 
+# Initialize data
 m, n = 7, 5
 irn = [1, 1, 1, 2, 3, 3, 4, 4, 5, 5, 6, 7, 7]
 jcn = [1, 3, 5, 2, 3, 5, 1, 4, 4, 5, 2, 1, 3]
@@ -41,20 +52,29 @@ A = sparse(irn, jcn, val, m, n)
 b = [22.0, 2.0, 13.0, 8.0, 17.0, 6.0, 5.0]
 x_star = [1.0, 2.0, 3.0, 4.0, 5.0]
 
+# Initialize QRMumps
 qrm_init()
 
+# Initialize data structures
 spmat = qrm_spmat_init(A)
 spfct = qrm_spfct_init(spmat)
 
+# Specify not storing Q for the QR factorization
 qrm_set(spfct, "qrm_keeph", 0)
 
+# Perform symbolic analysis of A and factorize A = QR
 qrm_analyse!(spmat, spfct)
 qrm_factorize!(spmat, spfct)
 
 z = A'*b
+
+# Solve Rᵀx₁ = z = Aᵀb
 x₁ = qrm_solve(spfct, z; transp = 't')
+
+# Solve Rx₂ = x₁ 
 x₂ = qrm_solve(spfct, x₁; transp = 'n')
 
+# Overall, RᵀRx₂ = Aᵀb. Equivalently, RᵀQᵀQRx₂ = Aᵀb or AᵀAx₂ = Aᵀb
 error_norm = norm(x₂ - x_star)
 
 @printf("Error norm ‖x* - x‖ = %10.5e\n", error_norm)

From a67eb9ed086e283411cfd6b1732dcfe8906b5f71 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Thu, 1 Aug 2024 13:00:41 +0200
Subject: [PATCH 04/19] add least norm solution to first example

---
 docs/src/tutorials/qless.md | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/docs/src/tutorials/qless.md b/docs/src/tutorials/qless.md
index 1d4ea3f..a0c9d9a 100644
--- a/docs/src/tutorials/qless.md
+++ b/docs/src/tutorials/qless.md
@@ -34,9 +34,16 @@ qrm_solve!(spfct, b, y₁; transp='t')
 qrm_solve!(spfct, y₁, y₂; transp='n')
 
 # Overall, RᵀRy₂ = b. Equivalently, RᵀQᵀQRy₂ = b or AAᵀy₂ = b
+
+# Compute least norm solution of min ‖b - Ax‖
+x = A'*y₂
+
+# Compute error norm and residual norm
 error_norm = norm(y₂ - y_star)
+residual_norm = norm(b - A*x)
 
 @printf("Error norm ‖y* - y‖ = %10.5e\n", error_norm)
+@printf("Residual norm ‖b - Ax‖ = %10.5e\n", residual_norm)
 ```
 
 ```@example qless2

From abc8368263e77a2fe0a7b88cd4fc02044ed2efaa Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Tue, 8 Oct 2024 11:01:53 -0400
Subject: [PATCH 05/19] add suggestions

---
 docs/src/tutorials/qless.md | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/docs/src/tutorials/qless.md b/docs/src/tutorials/qless.md
index a0c9d9a..8b3d562 100644
--- a/docs/src/tutorials/qless.md
+++ b/docs/src/tutorials/qless.md
@@ -9,9 +9,10 @@ val = [1.0, 2.0, 3.0, 1.0, 1.0, 2.0, 4.0, 1.0, 5.0, 1.0, 3.0, 6.0, 1.0]
 
 A = sparse(irn, jcn, val, m, n)
 b = [40.0, 10.0, 44.0, 98.0, 87.0]
-y_star = [0.0, 1.0, 2.0, 3.0, 4.0]
+y_star = [16.0, 1.0, 10.0, 3.0, 19.0, 3.0, 2.0]
 y₁ = zeros(n)
 y₂ = zeros(m)
+y₃ = zeros(n)
 
 # Initialize QRMumps
 qrm_init()
@@ -36,11 +37,11 @@ qrm_solve!(spfct, y₁, y₂; transp='n')
 # Overall, RᵀRy₂ = b. Equivalently, RᵀQᵀQRy₂ = b or AAᵀy₂ = b
 
 # Compute least norm solution of min ‖b - Ax‖
-x = A'*y₂
+y₃ .= A'*y₂
 
 # Compute error norm and residual norm
-error_norm = norm(y₂ - y_star)
-residual_norm = norm(b - A*x)
+error_norm = norm(y₃ - y_star)
+residual_norm = norm(b - A*y₃)
 
 @printf("Error norm ‖y* - y‖ = %10.5e\n", error_norm)
 @printf("Residual norm ‖b - Ax‖ = %10.5e\n", residual_norm)

From c3661bf28086b3fbbcd17035867c12deecaf5bbf Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Thu, 31 Oct 2024 17:01:03 +0000
Subject: [PATCH 06/19] update qless1 for iterative refinement

---
 docs/src/tutorials/qless.md | 35 +++++++++++++++++++++++++++++++++--
 1 file changed, 33 insertions(+), 2 deletions(-)

diff --git a/docs/src/tutorials/qless.md b/docs/src/tutorials/qless.md
index 8b3d562..dc869b2 100644
--- a/docs/src/tutorials/qless.md
+++ b/docs/src/tutorials/qless.md
@@ -14,6 +14,11 @@ y₁ = zeros(n)
 y₂ = zeros(m)
 y₃ = zeros(n)
 
+Δy₁ = zeros(n)
+Δy₂ = zeros(m)
+Δy₃ = zeros(n)
+r = zeros(m)
+
 # Initialize QRMumps
 qrm_init()
 
@@ -36,7 +41,7 @@ qrm_solve!(spfct, y₁, y₂; transp='n')
 
 # Overall, RᵀRy₂ = b. Equivalently, RᵀQᵀQRy₂ = b or AAᵀy₂ = b
 
-# Compute least norm solution of min ‖b - Ax‖
+# Compute least norm solution of min ‖b - Ay‖
 y₃ .= A'*y₂
 
 # Compute error norm and residual norm
@@ -44,7 +49,33 @@ error_norm = norm(y₃ - y_star)
 residual_norm = norm(b - A*y₃)
 
 @printf("Error norm ‖y* - y‖ = %10.5e\n", error_norm)
-@printf("Residual norm ‖b - Ax‖ = %10.5e\n", residual_norm)
+@printf("Residual norm ‖b - Ay‖ = %10.5e\n", residual_norm)
+
+# We can improve this solution with iterative refinement: see Björck 1967 - Iterative refinement of linear least squares solutions I
+#                                                             in BIT Numerical Mathematics.
+# For this, we compute the least norm solution Δy₃ of min ‖r - AΔy‖, where r is the residual r = b - A*y₃.
+# We then update y₃ := y₃ + Δy₃
+r .= b - A*y₃
+
+# Solve RᵀΔy₁ =  r 
+qrm_solve!(spfct, r, Δy₁; transp='t')
+
+# Solve RΔy₂ = Δy₁
+qrm_solve!(spfct, Δy₁, Δy₂; transp='n')
+
+# Overall, RᵀRΔy₂ = r. Equivalently, RᵀQᵀQRΔy₂ = r or AAᵀΔy₂ = r
+
+# Compute least norm solution of min ‖r - Ay‖
+Δy₃ .= A'*Δy₂
+
+# Update the least norm solution
+y₃ .= y₃ .+ Δy₃
+
+error_norm = norm(y₃ - y_star)
+residual_norm = norm(b - A*y₃)
+
+@printf("Error norm - iterative refinement ‖y* - y‖ = %10.5e\n", error_norm)
+@printf("Residual norm - iterative refinement ‖b - Ay‖ = %10.5e\n", residual_norm)
 ```
 
 ```@example qless2

From 69d7760d03529b64669b057bd759d3b2ddd2df2f Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Thu, 31 Oct 2024 17:12:46 +0000
Subject: [PATCH 07/19] add iterative refinement on second example

---
 docs/src/tutorials/qless.md | 25 +++++++++++++++++++++++++
 1 file changed, 25 insertions(+)

diff --git a/docs/src/tutorials/qless.md b/docs/src/tutorials/qless.md
index dc869b2..f704d4d 100644
--- a/docs/src/tutorials/qless.md
+++ b/docs/src/tutorials/qless.md
@@ -115,6 +115,31 @@ x₂ = qrm_solve(spfct, x₁; transp = 'n')
 
 # Overall, RᵀRx₂ = Aᵀb. Equivalently, RᵀQᵀQRx₂ = Aᵀb or AᵀAx₂ = Aᵀb
 error_norm = norm(x₂ - x_star)
+residual_norm = norm(A*x₂ - b)
 
 @printf("Error norm ‖x* - x‖ = %10.5e\n", error_norm)
+@printf("Residual norm ‖Ax - b‖= %10.5e\n", residual_norm)
+
+# We can improve this solution with iterative refinement: see Björck 1967 - Iterative refinement of linear least squares solutions I
+#                                                             in BIT Numerical Mathematics.
+# For this, we compute the least norm solution Δx₂ of min ‖r - AΔx‖, where r is the residual r = Aᵀb - AᵀA*x₂.
+# We then update x₂ := x₂ + Δx₂
+r = A'*(b - A*x₂)
+
+# Solve RᵀΔx₁ =  r 
+Δx₁ = qrm_solve(spfct, r; transp='t')
+
+# Solve Rxy₂ = Δx₁
+Δx₂ = qrm_solve(spfct, Δx₁; transp='n')
+
+# Overall, RᵀRΔx₂ = r. Equivalently, RᵀQᵀQRx₂ = r or AᵀAx₂ = r
+
+# Update the least squares solution
+x₂ .= x₂ .+ Δx₂
+
+error_norm = norm(x₂ - x_star)
+residual_norm = norm(A*x₂ - b)
+
+@printf("Error norm ‖x* - x‖ = %10.5e\n", error_norm)
+@printf("Residual norm ‖Ax - b‖= %10.5e\n", residual_norm)
 ```
\ No newline at end of file

From d22443bafe1a68ffbef146e655af1836a8391ad2 Mon Sep 17 00:00:00 2001
From: MaxenceGollier <134112149+MaxenceGollier@users.noreply.github.com>
Date: Fri, 8 Nov 2024 15:09:41 -0500
Subject: [PATCH 08/19] Apply suggestions from code review

Co-authored-by: Dominique <dominique.orban@gmail.com>
---
 docs/src/tutorials/qless.md | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/docs/src/tutorials/qless.md b/docs/src/tutorials/qless.md
index f704d4d..317c04f 100644
--- a/docs/src/tutorials/qless.md
+++ b/docs/src/tutorials/qless.md
@@ -1,5 +1,6 @@
 ```@example qless1
-using QRMumps, LinearAlgebra, SparseArrays, Printf
+using LinearAlgebra, Printf, SparseArrays
+using QRMumps
 
 # Initialize data
 m, n = 5, 7
@@ -26,20 +27,20 @@ qrm_init()
 spmat = qrm_spmat_init(A)
 spfct = qrm_spfct_init(spmat)
 
-# Specify not storing Q for the QR factorization
+# Specify that we want the Q-less QR factorization
 qrm_set(spfct, "qrm_keeph", 0)
 
 # Perform symbolic analysis of Aᵀ and factorize Aᵀ = QR
 qrm_analyse!(spmat, spfct; transp='t')
 qrm_factorize!(spmat, spfct, transp='t')
 
-# Solve Rᵀy₁ =  b  
+# Solve RᵀR y = b in two steps:
+# 1. Solve Rᵀy₁ =  b  
 qrm_solve!(spfct, b, y₁; transp='t')
 
-# Solve Ry₂ = y₁
-qrm_solve!(spfct, y₁, y₂; transp='n')
+# 2. Solve Ry = y₁
+qrm_solve!(spfct, y₁, y; transp='n')
 
-# Overall, RᵀRy₂ = b. Equivalently, RᵀQᵀQRy₂ = b or AAᵀy₂ = b
 
 # Compute least norm solution of min ‖b - Ay‖
 y₃ .= A'*y₂

From 237e7462d4b080d543f838b9ff4b09cc58ee7edb Mon Sep 17 00:00:00 2001
From: MaxenceGollier <134112149+MaxenceGollier@users.noreply.github.com>
Date: Fri, 8 Nov 2024 15:14:30 -0500
Subject: [PATCH 09/19] Update docs/src/tutorials/qless.md

Co-authored-by: Dominique <dominique.orban@gmail.com>
---
 docs/src/tutorials/qless.md | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/docs/src/tutorials/qless.md b/docs/src/tutorials/qless.md
index 317c04f..b3342a8 100644
--- a/docs/src/tutorials/qless.md
+++ b/docs/src/tutorials/qless.md
@@ -1,4 +1,18 @@
 ```@example qless1
+# The Q-less QR factorization may be used to solve the least-norm problem
+#
+#   minimize ‖x‖  subject to  Ax=b
+#
+# while saving storage because Q is not formed.
+# Thus it is appropriate for large problems where storage is at a premium.
+# The normal equations of the second kind AAᵀy = b are the optimality conditions of the least-norm problems, where x = Aᵀy.
+# If Aᵀ = QR, they can be equivalently written RᵀRy = b.
+#
+# This procedure is backward stable if we perform one step of iterative refinement---see
+#
+# Å. Björck, Stability analysis of the method of seminormal equations for linear least squares problems,
+# Linear Algebra and its Applications, 88–89, pp. 31-48, 1987, DOI 10.1016/0024-3795(87)90101-7.
+
 using LinearAlgebra, Printf, SparseArrays
 using QRMumps
 

From bed0e6ee6f3c6b39ef2b824fbcb4b9286940de60 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Mon, 11 Nov 2024 11:32:07 -0500
Subject: [PATCH 10/19] update first example, add citation to paige

---
 docs/src/tutorials/qless.md | 55 ++++++++-----------------------------
 1 file changed, 12 insertions(+), 43 deletions(-)

diff --git a/docs/src/tutorials/qless.md b/docs/src/tutorials/qless.md
index b3342a8..c905c57 100644
--- a/docs/src/tutorials/qless.md
+++ b/docs/src/tutorials/qless.md
@@ -8,10 +8,10 @@
 # The normal equations of the second kind AAᵀy = b are the optimality conditions of the least-norm problems, where x = Aᵀy.
 # If Aᵀ = QR, they can be equivalently written RᵀRy = b.
 #
-# This procedure is backward stable if we perform one step of iterative refinement---see
+# The stability of this procedure is comparable to the method that uses Q---see
 #
-# Å. Björck, Stability analysis of the method of seminormal equations for linear least squares problems,
-# Linear Algebra and its Applications, 88–89, pp. 31-48, 1987, DOI 10.1016/0024-3795(87)90101-7.
+# C. C. Paige, An error analysis of a method for solving matrix equations,
+# Mathematics of Computations, 27, pp. 355-359, 1973, DOI 10.2307/2005623.
 
 using LinearAlgebra, Printf, SparseArrays
 using QRMumps
@@ -24,15 +24,10 @@ val = [1.0, 2.0, 3.0, 1.0, 1.0, 2.0, 4.0, 1.0, 5.0, 1.0, 3.0, 6.0, 1.0]
 
 A = sparse(irn, jcn, val, m, n)
 b = [40.0, 10.0, 44.0, 98.0, 87.0]
-y_star = [16.0, 1.0, 10.0, 3.0, 19.0, 3.0, 2.0]
+x_star = [16.0, 1.0, 10.0, 3.0, 19.0, 3.0, 2.0]
 y₁ = zeros(n)
-y₂ = zeros(m)
-y₃ = zeros(n)
-
-Δy₁ = zeros(n)
-Δy₂ = zeros(m)
-Δy₃ = zeros(n)
-r = zeros(m)
+y = zeros(m)
+x = zeros(n)
 
 # Initialize QRMumps
 qrm_init()
@@ -56,41 +51,15 @@ qrm_solve!(spfct, b, y₁; transp='t')
 qrm_solve!(spfct, y₁, y; transp='n')
 
 
-# Compute least norm solution of min ‖b - Ay‖
-y₃ .= A'*y₂
+# Compute least norm solution of Ax = b
+x .= A'*y
 
 # Compute error norm and residual norm
-error_norm = norm(y₃ - y_star)
-residual_norm = norm(b - A*y₃)
-
-@printf("Error norm ‖y* - y‖ = %10.5e\n", error_norm)
-@printf("Residual norm ‖b - Ay‖ = %10.5e\n", residual_norm)
-
-# We can improve this solution with iterative refinement: see Björck 1967 - Iterative refinement of linear least squares solutions I
-#                                                             in BIT Numerical Mathematics.
-# For this, we compute the least norm solution Δy₃ of min ‖r - AΔy‖, where r is the residual r = b - A*y₃.
-# We then update y₃ := y₃ + Δy₃
-r .= b - A*y₃
-
-# Solve RᵀΔy₁ =  r 
-qrm_solve!(spfct, r, Δy₁; transp='t')
+error_norm = norm(x - x_star)
+residual_norm = norm(b - A*x)
 
-# Solve RΔy₂ = Δy₁
-qrm_solve!(spfct, Δy₁, Δy₂; transp='n')
-
-# Overall, RᵀRΔy₂ = r. Equivalently, RᵀQᵀQRΔy₂ = r or AAᵀΔy₂ = r
-
-# Compute least norm solution of min ‖r - Ay‖
-Δy₃ .= A'*Δy₂
-
-# Update the least norm solution
-y₃ .= y₃ .+ Δy₃
-
-error_norm = norm(y₃ - y_star)
-residual_norm = norm(b - A*y₃)
-
-@printf("Error norm - iterative refinement ‖y* - y‖ = %10.5e\n", error_norm)
-@printf("Residual norm - iterative refinement ‖b - Ay‖ = %10.5e\n", residual_norm)
+@printf("Error norm ‖x* - x‖ = %10.5e\n", error_norm)
+@printf("Residual norm ‖b - Ax‖ = %10.5e\n", residual_norm)
 ```
 
 ```@example qless2

From 717a0ccb42537a7d7a69cdc4090a3fc5baf7789f Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Mon, 11 Nov 2024 12:54:32 -0500
Subject: [PATCH 11/19] update second example

---
 docs/src/tutorials/qless.md | 66 ++++++++++++++++++++++---------------
 1 file changed, 40 insertions(+), 26 deletions(-)

diff --git a/docs/src/tutorials/qless.md b/docs/src/tutorials/qless.md
index c905c57..764aba8 100644
--- a/docs/src/tutorials/qless.md
+++ b/docs/src/tutorials/qless.md
@@ -44,14 +44,14 @@ qrm_analyse!(spmat, spfct; transp='t')
 qrm_factorize!(spmat, spfct, transp='t')
 
 # Solve RᵀR y = b in two steps:
-# 1. Solve Rᵀy₁ =  b  
+# 1. Solve Rᵀy₁ = b  
 qrm_solve!(spfct, b, y₁; transp='t')
 
 # 2. Solve Ry = y₁
 qrm_solve!(spfct, y₁, y; transp='n')
 
 
-# Compute least norm solution of Ax = b
+# Compute the least norm solution of Ax = b
 x .= A'*y
 
 # Compute error norm and residual norm
@@ -63,7 +63,22 @@ residual_norm = norm(b - A*x)
 ```
 
 ```@example qless2
-using QRMumps, LinearAlgebra, SparseArrays, Printf
+# The Q-less QR factorization may be used to solve the least-square problem
+#
+#   minimize ‖Ax - b‖
+#
+# while saving storage because Q is not formed.
+# Thus it is appropriate for large problems where storage is at a premium.
+# The normal equations AᵀAx = Aᵀb are the optimality conditions of the least-square problems.
+# If A = QR, they can be equivalently written RᵀRx = Aᵀb.
+#
+# This procedure is backward stable if we perform one step of iterative refinement---see
+#
+# Å. Björck, Stability analysis of the method of seminormal equations for linear least squares problems,
+# Linear Algebra and its Applications, 88–89, pp. 31-48, 1987, DOI 10.1016/0024-3795(87)90101-7.
+
+using LinearAlgebra, Printf, SparseArrays 
+using QRMumps
 
 # Initialize data
 m, n = 7, 5
@@ -82,48 +97,47 @@ qrm_init()
 spmat = qrm_spmat_init(A)
 spfct = qrm_spfct_init(spmat)
 
-# Specify not storing Q for the QR factorization
+# Specify that we want the Q-less QR factorization
 qrm_set(spfct, "qrm_keeph", 0)
 
 # Perform symbolic analysis of A and factorize A = QR
 qrm_analyse!(spmat, spfct)
 qrm_factorize!(spmat, spfct)
 
+# Compute the RHS of the semi-normal equations
 z = A'*b
 
-# Solve Rᵀx₁ = z = Aᵀb
+# Solve RᵀR x = z = Aᵀb in two steps:
+# 1. Solve Rᵀx₁ = z  
 x₁ = qrm_solve(spfct, z; transp = 't')
 
-# Solve Rx₂ = x₁ 
-x₂ = qrm_solve(spfct, x₁; transp = 'n')
+# 2. Solve Rx = x₁
+x = qrm_solve(spfct, x₁; transp = 'n')
 
-# Overall, RᵀRx₂ = Aᵀb. Equivalently, RᵀQᵀQRx₂ = Aᵀb or AᵀAx₂ = Aᵀb
-error_norm = norm(x₂ - x_star)
-residual_norm = norm(A*x₂ - b)
+
+error_norm = norm(x - x_star)
+residual_norm = norm(A*x - b)
 
 @printf("Error norm ‖x* - x‖ = %10.5e\n", error_norm)
 @printf("Residual norm ‖Ax - b‖= %10.5e\n", residual_norm)
 
-# We can improve this solution with iterative refinement: see Björck 1967 - Iterative refinement of linear least squares solutions I
-#                                                             in BIT Numerical Mathematics.
-# For this, we compute the least norm solution Δx₂ of min ‖r - AΔx‖, where r is the residual r = Aᵀb - AᵀA*x₂.
-# We then update x₂ := x₂ + Δx₂
-r = A'*(b - A*x₂)
-
-# Solve RᵀΔx₁ =  r 
-Δx₁ = qrm_solve(spfct, r; transp='t')
+# As such, this method is not backward stable and we need to add an iterative refinement step:                                                          
+# For this, we compute the least-square solution Δx of min ‖r - AΔx‖, where r is the residual r = Aᵀb - AᵀA*x.
+# We then update x := x + Δx
 
-# Solve Rxy₂ = Δx₁
-Δx₂ = qrm_solve(spfct, Δx₁; transp='n')
+# Compute the residual
+r = A'*(b - A*x)
 
-# Overall, RᵀRΔx₂ = r. Equivalently, RᵀQᵀQRx₂ = r or AᵀAx₂ = r
+# Solve the semi-normal equations as before
+Δx₁ = qrm_solve(spfct, r; transp='t')
+Δx = qrm_solve(spfct, Δx₁; transp='n')
 
 # Update the least squares solution
-x₂ .= x₂ .+ Δx₂
+x .= x .+ Δx
 
-error_norm = norm(x₂ - x_star)
-residual_norm = norm(A*x₂ - b)
+error_norm = norm(x - x_star)
+residual_norm = norm(A*x - b)
 
-@printf("Error norm ‖x* - x‖ = %10.5e\n", error_norm)
-@printf("Residual norm ‖Ax - b‖= %10.5e\n", residual_norm)
+@printf("Error norm (iterative refinement step) ‖x* - x‖ = %10.5e\n", error_norm)
+@printf("Residual norm (iterative refinement step) ‖Ax - b‖= %10.5e\n", residual_norm)
 ```
\ No newline at end of file

From a75766e009928eff02b7744e20858b56effb4532 Mon Sep 17 00:00:00 2001
From: MaxenceGollier <134112149+MaxenceGollier@users.noreply.github.com>
Date: Tue, 12 Nov 2024 10:37:27 -0500
Subject: [PATCH 12/19] Apply suggestions from code review

Co-authored-by: Dominique <dominique.orban@gmail.com>
---
 docs/src/tutorials/qless.md | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/docs/src/tutorials/qless.md b/docs/src/tutorials/qless.md
index 764aba8..3640ad6 100644
--- a/docs/src/tutorials/qless.md
+++ b/docs/src/tutorials/qless.md
@@ -114,7 +114,6 @@ x₁ = qrm_solve(spfct, z; transp = 't')
 # 2. Solve Rx = x₁
 x = qrm_solve(spfct, x₁; transp = 'n')
 
-
 error_norm = norm(x - x_star)
 residual_norm = norm(A*x - b)
 
@@ -122,7 +121,7 @@ residual_norm = norm(A*x - b)
 @printf("Residual norm ‖Ax - b‖= %10.5e\n", residual_norm)
 
 # As such, this method is not backward stable and we need to add an iterative refinement step:                                                          
-# For this, we compute the least-square solution Δx of min ‖r - AΔx‖, where r is the residual r = Aᵀb - AᵀA*x.
+# For this, we compute the least-squares solution Δx of min ‖r - AΔx‖, where r is the residual r = Aᵀb - AᵀA*x.
 # We then update x := x + Δx
 
 # Compute the residual

From 92f7afcf631da8a2d56328ed8c8b130277da3098 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Tue, 12 Nov 2024 11:01:27 -0500
Subject: [PATCH 13/19] make qless2 example non allocating

---
 docs/src/tutorials/qless.md | 32 +++++++++++++++++++++++---------
 1 file changed, 23 insertions(+), 9 deletions(-)

diff --git a/docs/src/tutorials/qless.md b/docs/src/tutorials/qless.md
index 3640ad6..9f2bf53 100644
--- a/docs/src/tutorials/qless.md
+++ b/docs/src/tutorials/qless.md
@@ -69,7 +69,7 @@ residual_norm = norm(b - A*x)
 #
 # while saving storage because Q is not formed.
 # Thus it is appropriate for large problems where storage is at a premium.
-# The normal equations AᵀAx = Aᵀb are the optimality conditions of the least-square problems.
+# The normal equations AᵀAx = Aᵀb are the optimality conditions of the least squares problems.
 # If A = QR, they can be equivalently written RᵀRx = Aᵀb.
 #
 # This procedure is backward stable if we perform one step of iterative refinement---see
@@ -90,6 +90,15 @@ A = sparse(irn, jcn, val, m, n)
 b = [22.0, 2.0, 13.0, 8.0, 17.0, 6.0, 5.0]
 x_star = [1.0, 2.0, 3.0, 4.0, 5.0]
 
+z = zeros(n)
+x₁ = zeros(m)
+x = zeros(n)
+
+y = zeros(m)
+r = zeros(n)
+Δx₁ = zeros(m)
+Δx = zeros(n)
+
 # Initialize QRMumps
 qrm_init()
 
@@ -105,14 +114,14 @@ qrm_analyse!(spmat, spfct)
 qrm_factorize!(spmat, spfct)
 
 # Compute the RHS of the semi-normal equations
-z = A'*b
+mul!(z, A', b)
 
 # Solve RᵀR x = z = Aᵀb in two steps:
 # 1. Solve Rᵀx₁ = z  
-x₁ = qrm_solve(spfct, z; transp = 't')
+qrm_solve!(spfct, z, x₁; transp = 't')
 
 # 2. Solve Rx = x₁
-x = qrm_solve(spfct, x₁; transp = 'n')
+qrm_solve!(spfct, x₁, x; transp = 'n')
 
 error_norm = norm(x - x_star)
 residual_norm = norm(A*x - b)
@@ -124,15 +133,20 @@ residual_norm = norm(A*x - b)
 # For this, we compute the least-squares solution Δx of min ‖r - AΔx‖, where r is the residual r = Aᵀb - AᵀA*x.
 # We then update x := x + Δx
 
-# Compute the residual
-r = A'*(b - A*x)
+# Compute the residual in two steps to prevent allocating memory:
+# 1. Compute y = b - Ax
+mul!(y, A, x)
+@. y = b - y
+
+# 2. Compute r = Aᵀy = Aᵀ(b - A*x)
+mul!(r, A', y)
 
 # Solve the semi-normal equations as before
-Δx₁ = qrm_solve(spfct, r; transp='t')
-Δx = qrm_solve(spfct, Δx₁; transp='n')
+qrm_solve!(spfct, r, Δx₁; transp='t')
+qrm_solve!(spfct, Δx₁, Δx; transp='n')
 
 # Update the least squares solution
-x .= x .+ Δx
+@. x = x + Δx
 
 error_norm = norm(x - x_star)
 residual_norm = norm(A*x - b)

From ed555a23e798d68fd941f8655098ad538bd097c1 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Tue, 12 Nov 2024 11:04:42 -0500
Subject: [PATCH 14/19] modify residual in qless2 example

---
 docs/src/tutorials/qless.md | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/docs/src/tutorials/qless.md b/docs/src/tutorials/qless.md
index 9f2bf53..1a83e58 100644
--- a/docs/src/tutorials/qless.md
+++ b/docs/src/tutorials/qless.md
@@ -124,10 +124,10 @@ qrm_solve!(spfct, z, x₁; transp = 't')
 qrm_solve!(spfct, x₁, x; transp = 'n')
 
 error_norm = norm(x - x_star)
-residual_norm = norm(A*x - b)
+residual_norm = norm(A'*(A*x - b))
 
 @printf("Error norm ‖x* - x‖ = %10.5e\n", error_norm)
-@printf("Residual norm ‖Ax - b‖= %10.5e\n", residual_norm)
+@printf("Residual norm ‖Aᵀ(Ax - b)‖= %10.5e\n", residual_norm)
 
 # As such, this method is not backward stable and we need to add an iterative refinement step:                                                          
 # For this, we compute the least-squares solution Δx of min ‖r - AΔx‖, where r is the residual r = Aᵀb - AᵀA*x.
@@ -149,8 +149,8 @@ qrm_solve!(spfct, Δx₁, Δx; transp='n')
 @. x = x + Δx
 
 error_norm = norm(x - x_star)
-residual_norm = norm(A*x - b)
+residual_norm = norm(A'*(A*x - b))
 
 @printf("Error norm (iterative refinement step) ‖x* - x‖ = %10.5e\n", error_norm)
-@printf("Residual norm (iterative refinement step) ‖Ax - b‖= %10.5e\n", residual_norm)
+@printf("Residual norm (iterative refinement step) ‖Aᵀ(Ax - b)‖= %10.5e\n", residual_norm)
 ```
\ No newline at end of file

From 27c26faa924b029e622c4debd18ad852fce39994 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Tue, 12 Nov 2024 11:29:14 -0500
Subject: [PATCH 15/19] change to 'normal equations residual norm' in qless2
 example

---
 docs/src/tutorials/qless.md | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/docs/src/tutorials/qless.md b/docs/src/tutorials/qless.md
index 1a83e58..6187b6e 100644
--- a/docs/src/tutorials/qless.md
+++ b/docs/src/tutorials/qless.md
@@ -124,10 +124,10 @@ qrm_solve!(spfct, z, x₁; transp = 't')
 qrm_solve!(spfct, x₁, x; transp = 'n')
 
 error_norm = norm(x - x_star)
-residual_norm = norm(A'*(A*x - b))
+normal_residual_norm = norm(A'*(A*x - b))
 
 @printf("Error norm ‖x* - x‖ = %10.5e\n", error_norm)
-@printf("Residual norm ‖Aᵀ(Ax - b)‖= %10.5e\n", residual_norm)
+@printf("Normal equations residual norm ‖Aᵀ(Ax - b)‖= %10.5e\n", normal_residual_norm)
 
 # As such, this method is not backward stable and we need to add an iterative refinement step:                                                          
 # For this, we compute the least-squares solution Δx of min ‖r - AΔx‖, where r is the residual r = Aᵀb - AᵀA*x.
@@ -149,8 +149,8 @@ qrm_solve!(spfct, Δx₁, Δx; transp='n')
 @. x = x + Δx
 
 error_norm = norm(x - x_star)
-residual_norm = norm(A'*(A*x - b))
+normal_residual_norm = norm(A'*(A*x - b))
 
 @printf("Error norm (iterative refinement step) ‖x* - x‖ = %10.5e\n", error_norm)
-@printf("Residual norm (iterative refinement step) ‖Aᵀ(Ax - b)‖= %10.5e\n", residual_norm)
+@printf("Normal equations residual norm (iterative refinement step) ‖Aᵀ(Ax - b)‖= %10.5e\n", normal_residual_norm)
 ```
\ No newline at end of file

From 807ef27b3593ef248f8cd0b220fee895aa0c31ee Mon Sep 17 00:00:00 2001
From: MaxenceGollier <134112149+MaxenceGollier@users.noreply.github.com>
Date: Tue, 12 Nov 2024 11:41:09 -0500
Subject: [PATCH 16/19] rename variable normal_residual_norm in qless2 example

---
 docs/src/tutorials/qless.md | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/docs/src/tutorials/qless.md b/docs/src/tutorials/qless.md
index 6187b6e..e2b22ef 100644
--- a/docs/src/tutorials/qless.md
+++ b/docs/src/tutorials/qless.md
@@ -124,10 +124,10 @@ qrm_solve!(spfct, z, x₁; transp = 't')
 qrm_solve!(spfct, x₁, x; transp = 'n')
 
 error_norm = norm(x - x_star)
-normal_residual_norm = norm(A'*(A*x - b))
+Aresidual_norm = norm(A'*(A*x - b))
 
 @printf("Error norm ‖x* - x‖ = %10.5e\n", error_norm)
-@printf("Normal equations residual norm ‖Aᵀ(Ax - b)‖= %10.5e\n", normal_residual_norm)
+@printf("Normal equations residual norm ‖Aᵀ(Ax - b)‖= %10.5e\n", Aresidual_norm)
 
 # As such, this method is not backward stable and we need to add an iterative refinement step:                                                          
 # For this, we compute the least-squares solution Δx of min ‖r - AΔx‖, where r is the residual r = Aᵀb - AᵀA*x.
@@ -149,8 +149,8 @@ qrm_solve!(spfct, Δx₁, Δx; transp='n')
 @. x = x + Δx
 
 error_norm = norm(x - x_star)
-normal_residual_norm = norm(A'*(A*x - b))
+Aresidual_norm = norm(A'*(A*x - b))
 
 @printf("Error norm (iterative refinement step) ‖x* - x‖ = %10.5e\n", error_norm)
-@printf("Normal equations residual norm (iterative refinement step) ‖Aᵀ(Ax - b)‖= %10.5e\n", normal_residual_norm)
-```
\ No newline at end of file
+@printf("Normal equations residual norm (iterative refinement step) ‖Aᵀ(Ax - b)‖= %10.5e\n", Aresidual_norm)
+```

From 9404cc8dafedb5295794c820fd0857885601701e Mon Sep 17 00:00:00 2001
From: MaxenceGollier <134112149+MaxenceGollier@users.noreply.github.com>
Date: Wed, 13 Nov 2024 15:27:36 -0500
Subject: [PATCH 17/19] Update docs/src/tutorials/qless.md

Co-authored-by: Dominique <dominique.orban@gmail.com>
---
 docs/src/tutorials/qless.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/tutorials/qless.md b/docs/src/tutorials/qless.md
index e2b22ef..931e617 100644
--- a/docs/src/tutorials/qless.md
+++ b/docs/src/tutorials/qless.md
@@ -69,7 +69,7 @@ residual_norm = norm(b - A*x)
 #
 # while saving storage because Q is not formed.
 # Thus it is appropriate for large problems where storage is at a premium.
-# The normal equations AᵀAx = Aᵀb are the optimality conditions of the least squares problems.
+# The normal equations AᵀAx = Aᵀb are the optimality conditions of the least-squares problem.
 # If A = QR, they can be equivalently written RᵀRx = Aᵀb.
 #
 # This procedure is backward stable if we perform one step of iterative refinement---see

From 67565fec012c151ebdb0d1e04d887762dc31c31b Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Wed, 13 Nov 2024 15:32:06 -0500
Subject: [PATCH 18/19] rename variables in qless example 2

---
 docs/src/tutorials/qless.md | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/docs/src/tutorials/qless.md b/docs/src/tutorials/qless.md
index 931e617..b0bce52 100644
--- a/docs/src/tutorials/qless.md
+++ b/docs/src/tutorials/qless.md
@@ -94,8 +94,8 @@ z = zeros(n)
 x₁ = zeros(m)
 x = zeros(n)
 
-y = zeros(m)
-r = zeros(n)
+r = zeros(m)
+Ar = zeros(n)
 Δx₁ = zeros(m)
 Δx = zeros(n)
 
@@ -133,16 +133,16 @@ Aresidual_norm = norm(A'*(A*x - b))
 # For this, we compute the least-squares solution Δx of min ‖r - AΔx‖, where r is the residual r = Aᵀb - AᵀA*x.
 # We then update x := x + Δx
 
-# Compute the residual in two steps to prevent allocating memory:
-# 1. Compute y = b - Ax
-mul!(y, A, x)
-@. y = b - y
+# Compute the normal equations residual in two steps to prevent allocating memory:
+# 1. Compute r = b - Ax
+mul!(r, A, x)
+@. r = b - r
 
-# 2. Compute r = Aᵀy = Aᵀ(b - A*x)
-mul!(r, A', y)
+# 2. Compute Aᵀr = Aᵀ(b - A*x)
+mul!(Ar, A', r)
 
 # Solve the semi-normal equations as before
-qrm_solve!(spfct, r, Δx₁; transp='t')
+qrm_solve!(spfct, Ar, Δx₁; transp='t')
 qrm_solve!(spfct, Δx₁, Δx; transp='n')
 
 # Update the least squares solution

From 963ba874532001e13744e5e6706e5109847fb6ce Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Wed, 13 Nov 2024 15:43:10 -0500
Subject: [PATCH 19/19] make examples non allocating when computing errors in
 qless

---
 docs/src/tutorials/qless.md | 42 +++++++++++++++++++++++++------------
 1 file changed, 29 insertions(+), 13 deletions(-)

diff --git a/docs/src/tutorials/qless.md b/docs/src/tutorials/qless.md
index b0bce52..7f54980 100644
--- a/docs/src/tutorials/qless.md
+++ b/docs/src/tutorials/qless.md
@@ -29,6 +29,9 @@ y₁ = zeros(n)
 y = zeros(m)
 x = zeros(n)
 
+e = zeros(n)
+r = zeros(m)
+
 # Initialize QRMumps
 qrm_init()
 
@@ -55,8 +58,12 @@ qrm_solve!(spfct, y₁, y; transp='n')
 x .= A'*y
 
 # Compute error norm and residual norm
-error_norm = norm(x - x_star)
-residual_norm = norm(b - A*x)
+@. e = x - x_star
+error_norm = norm(e)
+
+mul!(r, A, x)
+@. r = b - r
+residual_norm = norm(r)
 
 @printf("Error norm ‖x* - x‖ = %10.5e\n", error_norm)
 @printf("Residual norm ‖b - Ax‖ = %10.5e\n", residual_norm)
@@ -96,6 +103,7 @@ x = zeros(n)
 
 r = zeros(m)
 Ar = zeros(n)
+e = zeros(n) 
 Δx₁ = zeros(m)
 Δx = zeros(n)
 
@@ -123,15 +131,9 @@ qrm_solve!(spfct, z, x₁; transp = 't')
 # 2. Solve Rx = x₁
 qrm_solve!(spfct, x₁, x; transp = 'n')
 
-error_norm = norm(x - x_star)
-Aresidual_norm = norm(A'*(A*x - b))
-
-@printf("Error norm ‖x* - x‖ = %10.5e\n", error_norm)
-@printf("Normal equations residual norm ‖Aᵀ(Ax - b)‖= %10.5e\n", Aresidual_norm)
-
-# As such, this method is not backward stable and we need to add an iterative refinement step:                                                          
-# For this, we compute the least-squares solution Δx of min ‖r - AΔx‖, where r is the residual r = Aᵀb - AᵀA*x.
-# We then update x := x + Δx
+# Compute errors
+@. e = x - x_star
+error_norm = norm(e)
 
 # Compute the normal equations residual in two steps to prevent allocating memory:
 # 1. Compute r = b - Ax
@@ -140,6 +142,14 @@ mul!(r, A, x)
 
 # 2. Compute Aᵀr = Aᵀ(b - A*x)
 mul!(Ar, A', r)
+Aresidual_norm = norm(Ar)
+
+@printf("Error norm ‖x* - x‖ = %10.5e\n", error_norm)
+@printf("Normal equations residual norm ‖Aᵀ(Ax - b)‖= %10.5e\n", Aresidual_norm)
+
+# As such, this method is not backward stable and we need to add an iterative refinement step:                                                          
+# For this, we compute the least-squares solution Δx of min ‖Aᵀr - AΔx‖, where r is the residual r = b - A*x.
+# We then update x := x + Δx
 
 # Solve the semi-normal equations as before
 qrm_solve!(spfct, Ar, Δx₁; transp='t')
@@ -148,8 +158,14 @@ qrm_solve!(spfct, Δx₁, Δx; transp='n')
 # Update the least squares solution
 @. x = x + Δx
 
-error_norm = norm(x - x_star)
-Aresidual_norm = norm(A'*(A*x - b))
+# Compute errors as before
+@. e = x - x_star
+error_norm = norm(e)
+
+mul!(r, A, x)
+@. r = b - r
+mul!(Ar, A', r)
+Aresidual_norm = norm(Ar)
 
 @printf("Error norm (iterative refinement step) ‖x* - x‖ = %10.5e\n", error_norm)
 @printf("Normal equations residual norm (iterative refinement step) ‖Aᵀ(Ax - b)‖= %10.5e\n", Aresidual_norm)