From 3e212e0b3e2a7cd96ecbac731867db139e75b38d Mon Sep 17 00:00:00 2001
From: Maxence Gollier <134112149+MaxenceGollier@users.noreply.github.com>
Date: Tue, 26 Nov 2024 11:53:51 -0500
Subject: [PATCH 01/10] Add semi normal equations functions

---
 src/utils.jl | 25 ++++++++++++++++++++++++-
 1 file changed, 24 insertions(+), 1 deletion(-)

diff --git a/src/utils.jl b/src/utils.jl
index 7f37c95..3191ec4 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -22,4 +22,27 @@ function qrm_refine!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, x :: Abstract
   qrm_solve!(spfct, Δx, y, transp = transp)
   qrm_solve!(spfct, y, Δx, transp = 'n')
   @. x = x + Δx
-end
\ No newline at end of file
+end
+
+function qrm_min_norm_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, x :: AbstractVector{T}, b :: AbstractVector{T}, Δx :: AbstractVector{T}, y :: AbstractVector{T}) where T
+  #@assert length(x) == spfct.fct.n
+  transp = T <: Real ? 't' : 'c'
+  qrm_analyse!(spmat, spfct, transp = transp)
+  qrm_factorize!(spmat, spfct, transp = transp)
+  qrm_solve!(spfct, b, Δx, transp = transp)
+  qrm_solve!(spfct, Δx, y, transp = 'n')
+  #x = A^T y
+  qrm_spmat_mv!(spmat, T(1),  y, T(0), x, transp = transp)
+end
+
+function qrm_least_squares_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, x :: AbstractVector{T}, b :: AbstractVector{T}, z :: AbstractVector{T}, Δx :: AbstractVector{T}, y :: AbstractVector{T}) where T
+  transp = T <: Real ? 't' : 'c'
+  qrm_analyse!(spmat, spfct, transp  = 'n')
+  qrm_factorize!(spmat, spfct, transp = 'n')
+  # z = A^T b
+  qrm_spmat_mv!(spmat, T(1),  b, T(0), z, transp = transp)
+  qrm_solve!(spfct, z, y, transp = transp)
+  qrm_solve!(spfct, y, x, transp = 'n')
+
+  qrm_refine!(spmat, spfct, x, z, Δx, y)
+end

From 28f5a6afe8a863a6efefec285de1cade616bfef8 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <134112149+MaxenceGollier@users.noreply.github.com>
Date: Tue, 26 Nov 2024 11:54:26 -0500
Subject: [PATCH 02/10] link to qrm_semi_normal!

---
 src/QRMumps.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/QRMumps.jl b/src/QRMumps.jl
index f9001a1..03bf1aa 100644
--- a/src/QRMumps.jl
+++ b/src/QRMumps.jl
@@ -51,6 +51,7 @@ export qrm_spmat, qrm_spfct,
     qrm_spposv!, qrm_spposv,
     qrm_least_squares!, qrm_least_squares,
     qrm_min_norm!, qrm_min_norm,
+    qrm_min_norm_semi_normal!, qrm_least_squares_semi_normal!,
     qrm_residual_norm!, qrm_residual_norm,
     qrm_residual_orth!, qrm_residual_orth,
     qrm_refine!, qrm_refine, qrm_set, qrm_get

From e60062485b4e35088e86f510bf7cffc8f813151a Mon Sep 17 00:00:00 2001
From: Maxence Gollier <134112149+MaxenceGollier@users.noreply.github.com>
Date: Tue, 26 Nov 2024 11:55:25 -0500
Subject: [PATCH 03/10] Add unit tests

---
 test/test_qrm.jl | 32 ++++++++++++++++++++++++++++++++
 1 file changed, 32 insertions(+)

diff --git a/test/test_qrm.jl b/test/test_qrm.jl
index ac3fb6c..58dcfb8 100644
--- a/test/test_qrm.jl
+++ b/test/test_qrm.jl
@@ -544,5 +544,37 @@ end
     x = qrm_solve(spfct, x₁, transp = 'n')
     x_refined = qrm_refine(spmat, spfct, x, b)
     @test norm(b - A'*(A*x)) ≥ norm(b - A'*(A*x_refined))
+
+    tol = (real(T) == Float32) ? 1e-3 : 1e-12
+    qrm_init()
+    A = sprand(T, n, m, 0.3)
+    b = rand(T, n)
+
+    spmat = qrm_spmat_init(A)
+    spfct = qrm_spfct_init(spmat)
+    qrm_set(spfct, "qrm_keeph", 0)
+    x = zeros(T, m)
+    Δx = similar(x)
+    y = zeros(T, n)
+    qrm_min_norm_semi_normal!(spmat, spfct, x, b, Δx, y)
+
+    x2 = qrm_min_norm(spmat, b)
+
+    @test norm(A*x - b) ≤ tol
+    @test abs(norm(x) - norm(x2)) ≤ tol
+
+    A = sprand(T, m, n, 0.3)
+    spmat = qrm_spmat_init(A)
+    spfct = qrm_spfct_init(spmat)
+    qrm_set(spfct, "qrm_keeph", 0)
+    b = rand(T, m)
+    x = zeros(T, n)
+    Δx = similar(x)
+    z = similar(x)
+    y = zeros(T, m)
+    qrm_least_squares_semi_normal!(spmat, spfct, x, b, z, Δx, y)
+
+    x2 = qrm_least_squares(spmat, b)
+    @test abs(norm(x) - norm(x2)) ≤ tol
   end
 end

From 588b7ba83223ad1bb6dedf953da13cb9a6d2a86a Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Tue, 26 Nov 2024 12:29:05 -0500
Subject: [PATCH 04/10] add qrm_semi_normal

---
 src/QRMumps.jl   | 19 ++++++++++++++++++-
 src/utils.jl     | 39 ++++++++++++++++++++++++++++++++++++++-
 test/test_qrm.jl | 19 +++++--------------
 3 files changed, 61 insertions(+), 16 deletions(-)

diff --git a/src/QRMumps.jl b/src/QRMumps.jl
index 03bf1aa..13383c6 100644
--- a/src/QRMumps.jl
+++ b/src/QRMumps.jl
@@ -50,8 +50,9 @@ export qrm_spmat, qrm_spfct,
     qrm_spbackslash!, qrm_spbackslash, \,
     qrm_spposv!, qrm_spposv,
     qrm_least_squares!, qrm_least_squares,
+    qrm_least_squares_semi_normal!, qrm_least_squares_semi_normal,
     qrm_min_norm!, qrm_min_norm,
-    qrm_min_norm_semi_normal!, qrm_least_squares_semi_normal!,
+    qrm_min_norm_semi_normal!, qrm_min_norm_semi_normal,
     qrm_residual_norm!, qrm_residual_norm,
     qrm_residual_orth!, qrm_residual_orth,
     qrm_refine!, qrm_refine, qrm_set, qrm_get
@@ -349,6 +350,14 @@ function qrm_least_squares! end
 """
 function qrm_least_squares end
 
+@doc raw"""
+"""
+function qrm_least_squares_semi_normal! end
+
+@doc raw"""
+"""
+function qrm_least_squares_semi_normal end
+
 @doc raw"""
     qrm_min_norm!(spmat, b, x; transp='n')
 
@@ -380,6 +389,14 @@ function qrm_min_norm! end
 """
 function qrm_min_norm end
 
+@doc raw"""
+"""
+function qrm_min_norm_semi_normal! end
+
+@doc raw"""
+"""
+function qrm_min_norm_semi_normal end
+
 @doc raw"""
     qrm_residual_norm!(spmat, b, x, nrm; transp='n')
 
diff --git a/src/utils.jl b/src/utils.jl
index 3191ec4..2e0d164 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -24,9 +24,25 @@ function qrm_refine!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, x :: Abstract
   @. x = x + Δx
 end
 
+function qrm_min_norm_semi_normal(spmat :: qrm_spmat{T}, b :: AbstractVector{T}) where T
+  spfct = qrm_spfct_init(spmat)
+  qrm_set(spfct, "qrm_keeph", 0)
+
+  x = similar(b, spmat.mat.n)
+  Δx = similar(x)
+  y = similar(b)
+
+  qrm_min_norm_semi_normal!(spmat, spfct, x, b, Δx, y)
+  return x
+end
 function qrm_min_norm_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, x :: AbstractVector{T}, b :: AbstractVector{T}, Δx :: AbstractVector{T}, y :: AbstractVector{T}) where T
-  #@assert length(x) == spfct.fct.n
+  @assert length(x) == spmat.mat.n
+  @assert length(b) == spmat.mat.m
+  @assert length(Δx) == spmat.mat.n
+  @assert length(y) == spmat.mat.m
+
   transp = T <: Real ? 't' : 'c'
+
   qrm_analyse!(spmat, spfct, transp = transp)
   qrm_factorize!(spmat, spfct, transp = transp)
   qrm_solve!(spfct, b, Δx, transp = transp)
@@ -35,8 +51,29 @@ function qrm_min_norm_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T},
   qrm_spmat_mv!(spmat, T(1),  y, T(0), x, transp = transp)
 end
 
+function qrm_least_squares_semi_normal(spmat :: qrm_spmat{T}, b :: AbstractVector{T}) where T
+
+  spfct = qrm_spfct_init(spmat)
+  qrm_set(spfct, "qrm_keeph", 0)
+
+  x = similar(b, spmat.mat.n)
+  z = similar(x)
+  Δx = similar(x)
+  y = similar(b)
+
+  qrm_least_squares_semi_normal!(spmat, spfct, x, b, z, Δx, y)
+  return x
+end
+
 function qrm_least_squares_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, x :: AbstractVector{T}, b :: AbstractVector{T}, z :: AbstractVector{T}, Δx :: AbstractVector{T}, y :: AbstractVector{T}) where T
+  @assert length(x) == spmat.mat.n
+  @assert length(b) == spmat.mat.m
+  @assert length(z) == spmat.mat.n
+  @assert length(Δx) == spmat.mat.n
+  @assert length(y) == spmat.mat.m
+
   transp = T <: Real ? 't' : 'c'
+  
   qrm_analyse!(spmat, spfct, transp  = 'n')
   qrm_factorize!(spmat, spfct, transp = 'n')
   # z = A^T b
diff --git a/test/test_qrm.jl b/test/test_qrm.jl
index 58dcfb8..633ab10 100644
--- a/test/test_qrm.jl
+++ b/test/test_qrm.jl
@@ -551,30 +551,21 @@ end
     b = rand(T, n)
 
     spmat = qrm_spmat_init(A)
-    spfct = qrm_spfct_init(spmat)
-    qrm_set(spfct, "qrm_keeph", 0)
-    x = zeros(T, m)
-    Δx = similar(x)
-    y = zeros(T, n)
-    qrm_min_norm_semi_normal!(spmat, spfct, x, b, Δx, y)
 
+    x = qrm_min_norm_semi_normal(spmat, b)
     x2 = qrm_min_norm(spmat, b)
 
     @test norm(A*x - b) ≤ tol
     @test abs(norm(x) - norm(x2)) ≤ tol
 
     A = sprand(T, m, n, 0.3)
-    spmat = qrm_spmat_init(A)
-    spfct = qrm_spfct_init(spmat)
-    qrm_set(spfct, "qrm_keeph", 0)
     b = rand(T, m)
-    x = zeros(T, n)
-    Δx = similar(x)
-    z = similar(x)
-    y = zeros(T, m)
-    qrm_least_squares_semi_normal!(spmat, spfct, x, b, z, Δx, y)
 
+    spmat = qrm_spmat_init(A)
+
+    x = qrm_least_squares_semi_normal(spmat, b)
     x2 = qrm_least_squares(spmat, b)
+
     @test abs(norm(x) - norm(x2)) ≤ tol
   end
 end

From 04dd809468a7fdb15e5bab4f47d8265e5910e676 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Tue, 26 Nov 2024 13:38:25 -0500
Subject: [PATCH 05/10] add documentation for qrm_semi_normal

---
 docs/src/api.md |  4 ++++
 src/QRMumps.jl  | 61 +++++++++++++++++++++++++++++++++++++++++++++++++
 src/utils.jl    |  8 +++----
 3 files changed, 69 insertions(+), 4 deletions(-)

diff --git a/docs/src/api.md b/docs/src/api.md
index a43d10e..8362e2d 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -46,8 +46,12 @@ qrm_spposv
 qrm_spposv!
 qrm_least_squares
 qrm_least_squares!
+qrm_least_squares_semi_normal
+qrm_least_squares_semi_normal!
 qrm_min_norm
 qrm_min_norm!
+qrm_min_norm_semi_normal
+qrm_min_norm_semi_normal!
 qrm_refine
 qrm_refine!
 ```
diff --git a/src/QRMumps.jl b/src/QRMumps.jl
index 13383c6..7b9bbb1 100644
--- a/src/QRMumps.jl
+++ b/src/QRMumps.jl
@@ -351,10 +351,43 @@ function qrm_least_squares! end
 function qrm_least_squares end
 
 @doc raw"""
+qrm_least_squares!(spmat, b, x, z, Δx, y; transp='n')
+
+This function can be used to solve a linear least squares problem
+
+```math
+\min \|Ax − b\|_2
+```
+
+in the case where the input matrix is square or overdetermined.
+Contrary to `qrm_least_squares!`, this function allows to solve the problem without storing the Q-factor in the QR factorization of A.
+
+It is a shortcut for the sequence
+
+    qrm_analyse!(spmat, spfct, transp  = 'n')
+    qrm_factorize!(spmat, spfct, transp = 'n')
+    qrm_spmat_mv!(spmat, T(1),  b, T(0), z, transp = 't')
+    qrm_solve!(spfct, z, y, transp = 't')
+    qrm_solve!(spfct, y, x, transp = 'n')
+
+    qrm_refine!(spmat, spfct, x, z, Δx, y)
+
+Remark that the Q-factor is not used in this sequence but rather A and R. 
+
+#### Input Arguments :
+
+* `spmat`: the input matrix.
+* `spfct`: a sparse factorization object of type `qrm_spfct`.
+* `b`: the right-hand side.
+* `x`: the solution vector.
+* `Δx`: an auxiliary vector used to compute the solution, the size of this vector is the same as x.
+* `z`: an auxiliary vector used to store the value z = Aᵀb, the size of this vector is the same as x and can be unitialized when the function is called.
+* `y`: an auxiliary vector used to compute the solution, the size of this vector is the same as b.
 """
 function qrm_least_squares_semi_normal! end
 
 @doc raw"""
+    x = qrm_least_squares_semi_normal(spmat, b)
 """
 function qrm_least_squares_semi_normal end
 
@@ -390,10 +423,38 @@ function qrm_min_norm! end
 function qrm_min_norm end
 
 @doc raw"""
+    qrm_min_norm_semi_normal!(spmat, spfct, b, x, Δx, y; transp='n')
+
+This function can be used to solve a linear minimum norm problem
+
+```math
+\min \|x\|_2 \quad s.t. \quad Ax = b
+```
+in the case where the input matrix is square or underdetermined.
+Contrary to `qrm_min_norm!`, this function allows to solve the problem without storing the Q-factor in the QR factorization of Aᵀ.
+It is a shortcut for the sequence
+
+    qrm_analyse!(spmat, spfct, transp = 't')
+    qrm_factorize!(spmat, spfct, transp = 't')
+    qrm_solve!(spfct, b, Δx, transp = 't')
+    qrm_solve!(spfct, Δx, y, transp = 'n')
+    qrm_spmat_mv!(spmat, T(1),  y, T(0), x, transp = 't')
+
+Remark that the Q-factor is not used in this sequence but rather A and R. 
+
+#### Input Arguments :
+
+* `spmat`: the input matrix.
+* `spfct`: a sparse factorization object of type `qrm_spfct`.
+* `b`: the right-hand side.
+* `x`: the solution vector.
+* `Δx`: an auxiliary vector used to compute the solution, the size of this vector is the same as x.
+* `y`: an auxiliary vector used to compute the solution, the size of this vector is the same as b.
 """
 function qrm_min_norm_semi_normal! end
 
 @doc raw"""
+    x = qrm_min_norm_semi_normal(spmat, b)
 """
 function qrm_min_norm_semi_normal end
 
diff --git a/src/utils.jl b/src/utils.jl
index 2e0d164..4107a2c 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -32,10 +32,10 @@ function qrm_min_norm_semi_normal(spmat :: qrm_spmat{T}, b :: AbstractVector{T})
   Δx = similar(x)
   y = similar(b)
 
-  qrm_min_norm_semi_normal!(spmat, spfct, x, b, Δx, y)
+  qrm_min_norm_semi_normal!(spmat, spfct, b, x, Δx, y)
   return x
 end
-function qrm_min_norm_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, x :: AbstractVector{T}, b :: AbstractVector{T}, Δx :: AbstractVector{T}, y :: AbstractVector{T}) where T
+function qrm_min_norm_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, b :: AbstractVector{T}, x :: AbstractVector{T}, Δx :: AbstractVector{T}, y :: AbstractVector{T}) where T
   @assert length(x) == spmat.mat.n
   @assert length(b) == spmat.mat.m
   @assert length(Δx) == spmat.mat.n
@@ -61,11 +61,11 @@ function qrm_least_squares_semi_normal(spmat :: qrm_spmat{T}, b :: AbstractVecto
   Δx = similar(x)
   y = similar(b)
 
-  qrm_least_squares_semi_normal!(spmat, spfct, x, b, z, Δx, y)
+  qrm_least_squares_semi_normal!(spmat, spfct, b, x, z, Δx, y)
   return x
 end
 
-function qrm_least_squares_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, x :: AbstractVector{T}, b :: AbstractVector{T}, z :: AbstractVector{T}, Δx :: AbstractVector{T}, y :: AbstractVector{T}) where T
+function qrm_least_squares_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, b :: AbstractVector{T}, x :: AbstractVector{T}, z :: AbstractVector{T}, Δx :: AbstractVector{T}, y :: AbstractVector{T}) where T
   @assert length(x) == spmat.mat.n
   @assert length(b) == spmat.mat.m
   @assert length(z) == spmat.mat.n

From d8df44769c914ab98aac401136347ee681f8b18e Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Tue, 26 Nov 2024 14:26:57 -0500
Subject: [PATCH 06/10] add qrm_semi_normal and qrm_refine for matrix entries

---
 src/utils.jl     | 89 +++++++++++++++++++++++++++++++++++++++++++++---
 test/test_qrm.jl | 64 ++++++++++++++++++++++------------
 2 files changed, 126 insertions(+), 27 deletions(-)

diff --git a/src/utils.jl b/src/utils.jl
index 4107a2c..371b10c 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -7,6 +7,15 @@ function qrm_refine(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, x :: AbstractV
   return x_refined
 end
 
+function qrm_refine(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, x :: AbstractMatrix{T}, z :: AbstractMatrix{T}) where T
+  Δx = similar(x)
+  y = similar(x, (spfct.fct.m, size(x,2)))
+  x_refined = copy(x)
+
+  qrm_refine!(spmat, spfct, x_refined, z, Δx, y)
+  return x_refined
+end
+
 function qrm_refine!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, x :: AbstractVector{T}, z :: AbstractVector{T}, Δx :: AbstractVector{T}, y :: AbstractVector{T}) where T
   @assert length(x) == spfct.fct.n
   @assert length(z) == spfct.fct.n
@@ -24,17 +33,42 @@ function qrm_refine!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, x :: Abstract
   @. x = x + Δx
 end
 
-function qrm_min_norm_semi_normal(spmat :: qrm_spmat{T}, b :: AbstractVector{T}) where T
+function qrm_refine!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, x :: AbstractMatrix{T}, z :: AbstractMatrix{T}, Δx :: AbstractMatrix{T}, y :: AbstractMatrix{T}) where T
+  @assert size(x, 1) == spfct.fct.n
+  @assert size(z, 1) == spfct.fct.n
+  @assert size(Δx, 1) == spfct.fct.n
+  @assert size(y, 1) == spfct.fct.m
+  @assert size(x, 2) == size(z, 2)
+  @assert size(Δx, 2) == size(z, 2)
+  @assert size(y, 2) == size(z, 2)
+
+  transp = T <: Real ? 't' : 'c'
+
+  qrm_spmat_mv!(spmat, T(1), x, T(0), y, transp = 'n')
+  qrm_spmat_mv!(spmat, T(1), y, T(0), Δx, transp = transp)
+  @. Δx = z - Δx 
+
+  qrm_solve!(spfct, Δx, y, transp = transp)
+  qrm_solve!(spfct, y, Δx, transp = 'n')
+  @. x = x + Δx
+end
+
+function qrm_min_norm_semi_normal(spmat :: qrm_spmat{T}, b :: AbstractVecOrMat{T}) where T
   spfct = qrm_spfct_init(spmat)
   qrm_set(spfct, "qrm_keeph", 0)
 
-  x = similar(b, spmat.mat.n)
+  if typeof(b) <: AbstractVector{T}
+    x = similar(b, spmat.mat.n)
+  else 
+    x = similar(b, (spmat.mat.n, size(b, 2)))
+  end
   Δx = similar(x)
   y = similar(b)
 
   qrm_min_norm_semi_normal!(spmat, spfct, b, x, Δx, y)
   return x
 end
+
 function qrm_min_norm_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, b :: AbstractVector{T}, x :: AbstractVector{T}, Δx :: AbstractVector{T}, y :: AbstractVector{T}) where T
   @assert length(x) == spmat.mat.n
   @assert length(b) == spmat.mat.m
@@ -51,12 +85,34 @@ function qrm_min_norm_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T},
   qrm_spmat_mv!(spmat, T(1),  y, T(0), x, transp = transp)
 end
 
-function qrm_least_squares_semi_normal(spmat :: qrm_spmat{T}, b :: AbstractVector{T}) where T
+function qrm_min_norm_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, b :: AbstractMatrix{T}, x :: AbstractMatrix{T}, Δx :: AbstractMatrix{T}, y :: AbstractMatrix{T}) where T
+  @assert size(x, 1) == spmat.mat.n
+  @assert size(b, 1) == spmat.mat.m
+  @assert size(Δx, 1) == spmat.mat.n
+  @assert size(y, 1) == spmat.mat.m
+  @assert size(x, 2) == size(b, 2)
+  @assert size(Δx, 2) == size(b, 2)
+  @assert size(y, 2) == size(b, 2)
+
+  transp = T <: Real ? 't' : 'c'
+
+  qrm_analyse!(spmat, spfct, transp = transp)
+  qrm_factorize!(spmat, spfct, transp = transp)
+  qrm_solve!(spfct, b, Δx, transp = transp)
+  qrm_solve!(spfct, Δx, y, transp = 'n')
+  #x = A^T y
+  qrm_spmat_mv!(spmat, T(1),  y, T(0), x, transp = transp)
+end
+
+function qrm_least_squares_semi_normal(spmat :: qrm_spmat{T}, b :: AbstractVecOrMat{T}) where T
 
   spfct = qrm_spfct_init(spmat)
   qrm_set(spfct, "qrm_keeph", 0)
-
-  x = similar(b, spmat.mat.n)
+  if typeof(b) <: AbstractVector{T}
+    x = similar(b, spmat.mat.n)
+  else 
+    x = similar(b, (spmat.mat.n, size(b, 2)))
+  end
   z = similar(x)
   Δx = similar(x)
   y = similar(b)
@@ -83,3 +139,26 @@ function qrm_least_squares_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfc
 
   qrm_refine!(spmat, spfct, x, z, Δx, y)
 end
+
+function qrm_least_squares_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, b :: AbstractMatrix{T}, x :: AbstractMatrix{T}, z :: AbstractMatrix{T}, Δx :: AbstractMatrix{T}, y :: AbstractMatrix{T}) where T
+  @assert size(x, 1) == spmat.mat.n
+  @assert size(b, 1) == spmat.mat.m
+  @assert size(z, 1) == spmat.mat.n
+  @assert size(Δx, 1) == spmat.mat.n
+  @assert size(y, 1) == spmat.mat.m
+  @assert size(x, 2) == size(b, 2)
+  @assert size(z, 2) == size(b, 2)
+  @assert size(Δx, 2) == size(b, 2)
+  @assert size(y, 2) == size(b, 2)
+
+  transp = T <: Real ? 't' : 'c'
+  
+  qrm_analyse!(spmat, spfct, transp  = 'n')
+  qrm_factorize!(spmat, spfct, transp = 'n')
+  # z = A^T b
+  qrm_spmat_mv!(spmat, T(1),  b, T(0), z, transp = transp)
+  qrm_solve!(spfct, z, y, transp = transp)
+  qrm_solve!(spfct, y, x, transp = 'n')
+
+  qrm_refine!(spmat, spfct, x, z, Δx, y)
+end
diff --git a/test/test_qrm.jl b/test/test_qrm.jl
index 633ab10..fc800c5 100644
--- a/test/test_qrm.jl
+++ b/test/test_qrm.jl
@@ -83,20 +83,42 @@ p = 5
       r = b - A * x
       @test norm(A' * r) ≤ tol
 
+      x = qrm_least_squares_semi_normal(spmat, b)
+      r = b - A * x
+      @test norm(A' * r) ≤ tol
+
       X = qrm_least_squares(spmat, B)
       R = B - A * X
       @test norm(A' * R) ≤ tol
 
+      X = qrm_least_squares_semi_normal(spmat, B)
+      R = B - A * X
+      @test norm(A' * R) ≤ tol
+
       bc = copy(b)
       qrm_least_squares!(spmat, bc, x)
       r = b - A * x
       @test norm(A' * r) ≤ tol
 
+      z = similar(x)
+      Δx = similar(x)
+      y = similar(b)
+      x = qrm_least_squares_semi_normal!(spmat, spfct, b, x, z, Δx, y)
+      r = b - A * x
+      @test norm(A' * r) ≤ tol
+
       Bc = copy(B)
       qrm_least_squares!(spmat, Bc, X)
       R = B - A * X
       @test norm(A' * R) ≤ tol
 
+      Z = similar(X)
+      ΔX = similar(X)
+      Y = similar(B)
+      qrm_least_squares_semi_normal!(spmat, spfct, B, X, Z, ΔX, Y)
+      R = B - A * X
+      @test norm(A' * R) ≤ tol
+
       x = qrm_spbackslash(spmat, b)
       r = b - A * x
       @test norm(A' * r) ≤ tol
@@ -227,18 +249,38 @@ end
       r = b - A * x
       @test norm(r) ≤ tol
 
+      x = qrm_min_norm_semi_normal(spmat, b)
+      r = b - A * x
+      @test norm(r) ≤ tol
+
       X = qrm_min_norm(spmat, B)
       R = B - A * X
       @test norm(R) ≤ tol
+      
+      X = qrm_min_norm_semi_normal(spmat, B)
+      R = B - A * X
+      @test norm(R) ≤ tol
 
       qrm_min_norm!(spmat, b, x)
       r = b - A * x
       @test norm(r) ≤ tol
+      
+      Δx = similar(x)
+      y = similar(b)
+      qrm_min_norm_semi_normal!(spmat, spfct, b, x, Δx, y)
+      r = b - A * x
+      @test norm(r) ≤ tol
 
       qrm_min_norm!(spmat, B, X)
       R = B - A * X
       @test norm(R) ≤ tol
 
+      ΔX = similar(X)
+      Y = similar(B)
+      qrm_min_norm_semi_normal!(spmat, spfct, B, X, ΔX, Y)
+      R = B - A * X
+      @test norm(R) ≤ tol
+
       x = qrm_spbackslash(spmat, b)
       r = b - A * x
       @test norm(r) ≤ tol
@@ -545,27 +587,5 @@ end
     x_refined = qrm_refine(spmat, spfct, x, b)
     @test norm(b - A'*(A*x)) ≥ norm(b - A'*(A*x_refined))
 
-    tol = (real(T) == Float32) ? 1e-3 : 1e-12
-    qrm_init()
-    A = sprand(T, n, m, 0.3)
-    b = rand(T, n)
-
-    spmat = qrm_spmat_init(A)
-
-    x = qrm_min_norm_semi_normal(spmat, b)
-    x2 = qrm_min_norm(spmat, b)
-
-    @test norm(A*x - b) ≤ tol
-    @test abs(norm(x) - norm(x2)) ≤ tol
-
-    A = sprand(T, m, n, 0.3)
-    b = rand(T, m)
-
-    spmat = qrm_spmat_init(A)
-
-    x = qrm_least_squares_semi_normal(spmat, b)
-    x2 = qrm_least_squares(spmat, b)
-
-    @test abs(norm(x) - norm(x2)) ≤ tol
   end
 end

From 8d43d6926d4a2521a1a20b40fef7730db0f52ed7 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Tue, 26 Nov 2024 14:37:01 -0500
Subject: [PATCH 07/10] update utils documentation for matrices input

---
 src/QRMumps.jl | 28 ++++++++++++++--------------
 1 file changed, 14 insertions(+), 14 deletions(-)

diff --git a/src/QRMumps.jl b/src/QRMumps.jl
index 7b9bbb1..4e0dcd0 100644
--- a/src/QRMumps.jl
+++ b/src/QRMumps.jl
@@ -351,7 +351,7 @@ function qrm_least_squares! end
 function qrm_least_squares end
 
 @doc raw"""
-qrm_least_squares!(spmat, b, x, z, Δx, y; transp='n')
+qrm_least_squares_semi_normal!(spmat, b, x, z, Δx, y; transp='n')
 
 This function can be used to solve a linear least squares problem
 
@@ -378,11 +378,11 @@ Remark that the Q-factor is not used in this sequence but rather A and R.
 
 * `spmat`: the input matrix.
 * `spfct`: a sparse factorization object of type `qrm_spfct`.
-* `b`: the right-hand side.
-* `x`: the solution vector.
-* `Δx`: an auxiliary vector used to compute the solution, the size of this vector is the same as x.
-* `z`: an auxiliary vector used to store the value z = Aᵀb, the size of this vector is the same as x and can be unitialized when the function is called.
-* `y`: an auxiliary vector used to compute the solution, the size of this vector is the same as b.
+* `b`: the right-hand side(s).
+* `x`: the solution vector(s).
+* `Δx`: an auxiliary vector (or matrix if b and x are matrices) used to compute the solution, the size of this vector (resp. matrix) is the same as x.
+* `z`: an auxiliary vector (or matrix if b and x are matrices) used to store the value z = Aᵀb, the size of this vector (resp. matrix) is the same as x and can be unitialized when the function is called.
+* `y`: an auxiliary vector (or matrix if b and x are matrices) used to compute the solution, the size of this vector (resp. matrix) is the same as b.
 """
 function qrm_least_squares_semi_normal! end
 
@@ -446,10 +446,10 @@ Remark that the Q-factor is not used in this sequence but rather A and R.
 
 * `spmat`: the input matrix.
 * `spfct`: a sparse factorization object of type `qrm_spfct`.
-* `b`: the right-hand side.
-* `x`: the solution vector.
-* `Δx`: an auxiliary vector used to compute the solution, the size of this vector is the same as x.
-* `y`: an auxiliary vector used to compute the solution, the size of this vector is the same as b.
+* `b`: the right-hand side(s).
+* `x`: the solution vector(s).
+* `Δx`: an auxiliary vector (or matrix if x and b are matrices) used to compute the solution, the size of this vector (resp. matrix) is the same as x.
+* `y`: an auxiliary vector (or matrix if x and b are matrices) used to compute the solution, the size of this vector (resp. matrix) is the same as b.
 """
 function qrm_min_norm_semi_normal! end
 
@@ -506,10 +506,10 @@ Given an approximate solution x of the linear system RᵀRx ≈ z where R is the
 
 * `spmat`: the input matrix.
 * `spfct`: a sparse factorization object of type `qrm_spfct`.
-* `x`: the approximate solution vector, the size of this vector is n.
-* `z`: the RHS vector of the linear system, the size of this vector is n.
-* `Δx`: an auxiliary vector used to compute the refinement, the size of this vector is n.
-* `y`: an auxiliary vector used to compute the refinement, the size of this vector is m.
+* `x`: the approximate solution vector(s), the size of this vector is n (or n×k if there are multiple solutions).
+* `z`: the RHS vector(s) of the linear system, the size of this vector is n (or n×k if there are multiple RHS).
+* `Δx`: an auxiliary vector (or matrix if x and z are matrices) used to compute the refinement, the size of this vector (resp. matrix) is n (resp. n×k).
+* `y`: an auxiliary vector (or matrix if x and z are matrices) used to compute the refinement, the size of this vector (resp. matrix) is m (resp. m×k).
 """
 function qrm_refine! end
 

From 186f07cac247492f2c6408c3ebd252e31898b2d3 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Thu, 19 Dec 2024 16:31:54 +0100
Subject: [PATCH 08/10] add transpose keyword argument for
 qrm_least_squares_semi_normal

---
 src/QRMumps.jl   |  1 +
 src/utils.jl     | 98 ++++++++++++++++++++++++++++--------------------
 test/test_qrm.jl | 15 ++++++++
 3 files changed, 73 insertions(+), 41 deletions(-)

diff --git a/src/QRMumps.jl b/src/QRMumps.jl
index 4e0dcd0..942fae6 100644
--- a/src/QRMumps.jl
+++ b/src/QRMumps.jl
@@ -383,6 +383,7 @@ Remark that the Q-factor is not used in this sequence but rather A and R.
 * `Δx`: an auxiliary vector (or matrix if b and x are matrices) used to compute the solution, the size of this vector (resp. matrix) is the same as x.
 * `z`: an auxiliary vector (or matrix if b and x are matrices) used to store the value z = Aᵀb, the size of this vector (resp. matrix) is the same as x and can be unitialized when the function is called.
 * `y`: an auxiliary vector (or matrix if b and x are matrices) used to compute the solution, the size of this vector (resp. matrix) is the same as b.
+* `transp`: whether to use A, Aᵀ or Aᴴ. Can be either `'t'`, `'c'` or `'n'`.
 """
 function qrm_least_squares_semi_normal! end
 
diff --git a/src/utils.jl b/src/utils.jl
index 371b10c..fe48682 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -22,13 +22,15 @@ function qrm_refine!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, x :: Abstract
   @assert length(Δx) == spfct.fct.n
   @assert length(y) == spfct.fct.m
 
-  transp = T <: Real ? 't' : 'c'
+  t = T <: Real ? 't' : 'c'
+  transp = spfct.fct.n == spmat.mat.n ? 'n' : t # Check whether it is A^T = QR or A = QR.
+  ntransp = transp == 't' || transp == 'c' ? 'n' : t
 
-  qrm_spmat_mv!(spmat, T(1), x, T(0), y, transp = 'n')
-  qrm_spmat_mv!(spmat, T(1), y, T(0), Δx, transp = transp)
+  qrm_spmat_mv!(spmat, T(1), x, T(0), y, transp = transp)
+  qrm_spmat_mv!(spmat, T(1), y, T(0), Δx, transp = ntransp)
   @. Δx = z - Δx 
 
-  qrm_solve!(spfct, Δx, y, transp = transp)
+  qrm_solve!(spfct, Δx, y, transp = t)
   qrm_solve!(spfct, y, Δx, transp = 'n')
   @. x = x + Δx
 end
@@ -42,18 +44,20 @@ function qrm_refine!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, x :: Abstract
   @assert size(Δx, 2) == size(z, 2)
   @assert size(y, 2) == size(z, 2)
 
-  transp = T <: Real ? 't' : 'c'
+  t = T <: Real ? 't' : 'c'
+  transp = spfct.fct.n == spmat.mat.n ? 'n' : t # Check whether it is A^T = QR or A = QR.
+  ntransp = transp == 't' || transp == 'c' ? 'n' : t
 
-  qrm_spmat_mv!(spmat, T(1), x, T(0), y, transp = 'n')
-  qrm_spmat_mv!(spmat, T(1), y, T(0), Δx, transp = transp)
+  qrm_spmat_mv!(spmat, T(1), x, T(0), y, transp = transp)
+  qrm_spmat_mv!(spmat, T(1), y, T(0), Δx, transp = ntransp)
   @. Δx = z - Δx 
 
-  qrm_solve!(spfct, Δx, y, transp = transp)
+  qrm_solve!(spfct, Δx, y, transp = t)
   qrm_solve!(spfct, y, Δx, transp = 'n')
   @. x = x + Δx
 end
 
-function qrm_min_norm_semi_normal(spmat :: qrm_spmat{T}, b :: AbstractVecOrMat{T}) where T
+function qrm_min_norm_semi_normal(spmat :: qrm_spmat{T}, b :: AbstractVecOrMat{T}; transp :: Char = 'n') where T
   spfct = qrm_spfct_init(spmat)
   qrm_set(spfct, "qrm_keeph", 0)
 
@@ -69,7 +73,7 @@ function qrm_min_norm_semi_normal(spmat :: qrm_spmat{T}, b :: AbstractVecOrMat{T
   return x
 end
 
-function qrm_min_norm_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, b :: AbstractVector{T}, x :: AbstractVector{T}, Δx :: AbstractVector{T}, y :: AbstractVector{T}) where T
+function qrm_min_norm_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, b :: AbstractVector{T}, x :: AbstractVector{T}, Δx :: AbstractVector{T}, y :: AbstractVector{T}; transp :: Char = 'n') where T
   @assert length(x) == spmat.mat.n
   @assert length(b) == spmat.mat.m
   @assert length(Δx) == spmat.mat.n
@@ -85,7 +89,7 @@ function qrm_min_norm_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T},
   qrm_spmat_mv!(spmat, T(1),  y, T(0), x, transp = transp)
 end
 
-function qrm_min_norm_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, b :: AbstractMatrix{T}, x :: AbstractMatrix{T}, Δx :: AbstractMatrix{T}, y :: AbstractMatrix{T}) where T
+function qrm_min_norm_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, b :: AbstractMatrix{T}, x :: AbstractMatrix{T}, Δx :: AbstractMatrix{T}, y :: AbstractMatrix{T}; transp :: Char = 'n') where T
   @assert size(x, 1) == spmat.mat.n
   @assert size(b, 1) == spmat.mat.m
   @assert size(Δx, 1) == spmat.mat.n
@@ -104,61 +108,73 @@ function qrm_min_norm_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T},
   qrm_spmat_mv!(spmat, T(1),  y, T(0), x, transp = transp)
 end
 
-function qrm_least_squares_semi_normal(spmat :: qrm_spmat{T}, b :: AbstractVecOrMat{T}) where T
+function qrm_least_squares_semi_normal(spmat :: qrm_spmat{T}, b :: AbstractVecOrMat{T}; transp :: Char = 'n') where T
+
+  n = transp == 'n' ? spmat.mat.n : spmat.mat.m
 
   spfct = qrm_spfct_init(spmat)
   qrm_set(spfct, "qrm_keeph", 0)
   if typeof(b) <: AbstractVector{T}
-    x = similar(b, spmat.mat.n)
+    x = similar(b, n)
   else 
-    x = similar(b, (spmat.mat.n, size(b, 2)))
+    x = similar(b, (n, size(b, 2)))
   end
   z = similar(x)
   Δx = similar(x)
   y = similar(b)
 
-  qrm_least_squares_semi_normal!(spmat, spfct, b, x, z, Δx, y)
+  qrm_least_squares_semi_normal!(spmat, spfct, b, x, z, Δx, y, transp = transp)
   return x
 end
 
-function qrm_least_squares_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, b :: AbstractVector{T}, x :: AbstractVector{T}, z :: AbstractVector{T}, Δx :: AbstractVector{T}, y :: AbstractVector{T}) where T
-  @assert length(x) == spmat.mat.n
-  @assert length(b) == spmat.mat.m
-  @assert length(z) == spmat.mat.n
-  @assert length(Δx) == spmat.mat.n
-  @assert length(y) == spmat.mat.m
-
-  transp = T <: Real ? 't' : 'c'
+function qrm_least_squares_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, b :: AbstractVector{T}, x :: AbstractVector{T}, z :: AbstractVector{T}, Δx :: AbstractVector{T}, y :: AbstractVector{T}; transp :: Char = 'n') where T
   
-  qrm_analyse!(spmat, spfct, transp  = 'n')
-  qrm_factorize!(spmat, spfct, transp = 'n')
+  n = transp == 'n' ? spmat.mat.n : spmat.mat.m
+  m = transp == 'n' ? spmat.mat.m : spmat.mat.n
+  t = T <: Real ? 't' : 'c'
+  ntransp = transp == 't' || transp == 'c' ? 'n' : t
+
+  @assert length(x) == n
+  @assert length(b) == m
+  @assert length(z) == n
+  @assert length(Δx) == n
+  @assert length(y) == m
+
+  qrm_analyse!(spmat, spfct, transp  = transp)
+  qrm_factorize!(spmat, spfct, transp = transp)
   # z = A^T b
-  qrm_spmat_mv!(spmat, T(1),  b, T(0), z, transp = transp)
-  qrm_solve!(spfct, z, y, transp = transp)
+  qrm_spmat_mv!(spmat, T(1),  b, T(0), z, transp = ntransp)
+  #R^T y = z
+  qrm_solve!(spfct, z, y, transp = t)
   qrm_solve!(spfct, y, x, transp = 'n')
 
-  qrm_refine!(spmat, spfct, x, z, Δx, y)
+  qrm_refine!(spmat, spfct, x, z, Δx, y, transp = transp)
 end
 
-function qrm_least_squares_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, b :: AbstractMatrix{T}, x :: AbstractMatrix{T}, z :: AbstractMatrix{T}, Δx :: AbstractMatrix{T}, y :: AbstractMatrix{T}) where T
-  @assert size(x, 1) == spmat.mat.n
-  @assert size(b, 1) == spmat.mat.m
-  @assert size(z, 1) == spmat.mat.n
-  @assert size(Δx, 1) == spmat.mat.n
-  @assert size(y, 1) == spmat.mat.m
+function qrm_least_squares_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, b :: AbstractMatrix{T}, x :: AbstractMatrix{T}, z :: AbstractMatrix{T}, Δx :: AbstractMatrix{T}, y :: AbstractMatrix{T}; transp :: Char = 'n') where T
+  
+  n = transp == 'n' ? spmat.mat.n : spmat.mat.m
+  m = transp == 'n' ? spmat.mat.m : spmat.mat.n
+  t = T <: Real ? 't' : 'c'
+  ntransp = transp == 't' || transp == 'c' ? 'n' : t
+
+  @assert size(x, 1) == n
+  @assert size(b, 1) == m
+  @assert size(z, 1) == n
+  @assert size(Δx, 1) == n
+  @assert size(y, 1) == m
+
   @assert size(x, 2) == size(b, 2)
   @assert size(z, 2) == size(b, 2)
   @assert size(Δx, 2) == size(b, 2)
   @assert size(y, 2) == size(b, 2)
 
-  transp = T <: Real ? 't' : 'c'
-  
-  qrm_analyse!(spmat, spfct, transp  = 'n')
-  qrm_factorize!(spmat, spfct, transp = 'n')
+  qrm_analyse!(spmat, spfct, transp  = transp)
+  qrm_factorize!(spmat, spfct, transp = transp)
   # z = A^T b
-  qrm_spmat_mv!(spmat, T(1),  b, T(0), z, transp = transp)
-  qrm_solve!(spfct, z, y, transp = transp)
+  qrm_spmat_mv!(spmat, T(1),  b, T(0), z, transp = ntransp)
+  qrm_solve!(spfct, z, y, transp = t)
   qrm_solve!(spfct, y, x, transp = 'n')
 
-  qrm_refine!(spmat, spfct, x, z, Δx, y)
+  qrm_refine!(spmat, spfct, x, z, Δx, y, transp = transp)
 end
diff --git a/test/test_qrm.jl b/test/test_qrm.jl
index fc800c5..a1dc620 100644
--- a/test/test_qrm.jl
+++ b/test/test_qrm.jl
@@ -12,6 +12,10 @@ p = 5
     for I in (Int32 , Int64)
       A = sprand(T, m, n, 0.3)
       A = convert(SparseMatrixCSC{T,I}, A)
+
+      A_transp = sprand(T, n, m, 0.3)
+      A_transp = convert(SparseMatrixCSC{T,I}, A_transp)
+
       b = rand(T, m)
       B = rand(T, m, p)
 
@@ -43,6 +47,9 @@ p = 5
       spmat = qrm_spmat_init(T)
       qrm_spmat_init!(spmat, A)
 
+      spmat_transp = qrm_spmat_init(T)
+      qrm_spmat_init!(spmat_transp, A_transp)
+
       spfct = qrm_spfct_init(spmat)
       qrm_analyse!(spmat, spfct)
       qrm_factorize!(spmat, spfct)
@@ -87,6 +94,10 @@ p = 5
       r = b - A * x
       @test norm(A' * r) ≤ tol
 
+      x = qrm_least_squares_semi_normal(spmat_transp, b, transp = transp)
+      r = b - A_transp' * x
+      @test norm(A_transp * r) ≤ tol
+
       X = qrm_least_squares(spmat, B)
       R = B - A * X
       @test norm(A' * R) ≤ tol
@@ -95,6 +106,10 @@ p = 5
       R = B - A * X
       @test norm(A' * R) ≤ tol
 
+      X = qrm_least_squares_semi_normal(spmat_transp, B, transp = transp)
+      R = B - A_transp' * X
+      @test norm(A_transp * r) ≤ tol
+
       bc = copy(b)
       qrm_least_squares!(spmat, bc, x)
       r = b - A * x

From ab9b7a5af3f01e63b59434e4bf5f7306c1f53732 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Thu, 19 Dec 2024 16:48:50 +0100
Subject: [PATCH 09/10] add transpose keyword argument for
 qrm_min_norm_semi_normal

---
 src/QRMumps.jl   |  1 +
 src/utils.jl     | 61 ++++++++++++++++++++++++++++--------------------
 test/test_qrm.jl | 17 ++++++++++++++
 3 files changed, 54 insertions(+), 25 deletions(-)

diff --git a/src/QRMumps.jl b/src/QRMumps.jl
index 942fae6..18e4905 100644
--- a/src/QRMumps.jl
+++ b/src/QRMumps.jl
@@ -451,6 +451,7 @@ Remark that the Q-factor is not used in this sequence but rather A and R.
 * `x`: the solution vector(s).
 * `Δx`: an auxiliary vector (or matrix if x and b are matrices) used to compute the solution, the size of this vector (resp. matrix) is the same as x.
 * `y`: an auxiliary vector (or matrix if x and b are matrices) used to compute the solution, the size of this vector (resp. matrix) is the same as b.
+* `transp`: whether to use A, Aᵀ or Aᴴ. Can be either `'t'`, `'c'` or `'n'`.
 """
 function qrm_min_norm_semi_normal! end
 
diff --git a/src/utils.jl b/src/utils.jl
index fe48682..40b9d86 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -58,54 +58,65 @@ function qrm_refine!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, x :: Abstract
 end
 
 function qrm_min_norm_semi_normal(spmat :: qrm_spmat{T}, b :: AbstractVecOrMat{T}; transp :: Char = 'n') where T
+
+  n = transp == 'n' ? spmat.mat.n : spmat.mat.m
+
   spfct = qrm_spfct_init(spmat)
   qrm_set(spfct, "qrm_keeph", 0)
 
   if typeof(b) <: AbstractVector{T}
-    x = similar(b, spmat.mat.n)
+    x = similar(b, n)
   else 
-    x = similar(b, (spmat.mat.n, size(b, 2)))
+    x = similar(b, (n, size(b, 2)))
   end
   Δx = similar(x)
   y = similar(b)
 
-  qrm_min_norm_semi_normal!(spmat, spfct, b, x, Δx, y)
+  qrm_min_norm_semi_normal!(spmat, spfct, b, x, Δx, y, transp = transp)
   return x
 end
 
 function qrm_min_norm_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, b :: AbstractVector{T}, x :: AbstractVector{T}, Δx :: AbstractVector{T}, y :: AbstractVector{T}; transp :: Char = 'n') where T
-  @assert length(x) == spmat.mat.n
-  @assert length(b) == spmat.mat.m
-  @assert length(Δx) == spmat.mat.n
-  @assert length(y) == spmat.mat.m
-
-  transp = T <: Real ? 't' : 'c'
+  
+  n = transp == 'n' ? spmat.mat.n : spmat.mat.m
+  m = transp == 'n' ? spmat.mat.m : spmat.mat.n
+  t = T <: Real ? 't' : 'c'
+  ntransp = transp == 't' || transp == 'c' ? 'n' : t
+  
+  @assert length(x) == n
+  @assert length(b) == m
+  @assert length(Δx) == n
+  @assert length(y) == m
 
-  qrm_analyse!(spmat, spfct, transp = transp)
-  qrm_factorize!(spmat, spfct, transp = transp)
-  qrm_solve!(spfct, b, Δx, transp = transp)
+  qrm_analyse!(spmat, spfct, transp = ntransp)
+  qrm_factorize!(spmat, spfct, transp = ntransp)
+  qrm_solve!(spfct, b, Δx, transp = t)
   qrm_solve!(spfct, Δx, y, transp = 'n')
   #x = A^T y
-  qrm_spmat_mv!(spmat, T(1),  y, T(0), x, transp = transp)
+  qrm_spmat_mv!(spmat, T(1),  y, T(0), x, transp = ntransp)
 end
 
 function qrm_min_norm_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, b :: AbstractMatrix{T}, x :: AbstractMatrix{T}, Δx :: AbstractMatrix{T}, y :: AbstractMatrix{T}; transp :: Char = 'n') where T
-  @assert size(x, 1) == spmat.mat.n
-  @assert size(b, 1) == spmat.mat.m
-  @assert size(Δx, 1) == spmat.mat.n
-  @assert size(y, 1) == spmat.mat.m
+  
+  n = transp == 'n' ? spmat.mat.n : spmat.mat.m
+  m = transp == 'n' ? spmat.mat.m : spmat.mat.n
+  t = T <: Real ? 't' : 'c'
+  ntransp = transp == 't' || transp == 'c' ? 'n' : t
+  
+  @assert size(x, 1) == n
+  @assert size(b, 1) == m
+  @assert size(Δx, 1) == n
+  @assert size(y, 1) == m
   @assert size(x, 2) == size(b, 2)
   @assert size(Δx, 2) == size(b, 2)
   @assert size(y, 2) == size(b, 2)
 
-  transp = T <: Real ? 't' : 'c'
-
-  qrm_analyse!(spmat, spfct, transp = transp)
-  qrm_factorize!(spmat, spfct, transp = transp)
-  qrm_solve!(spfct, b, Δx, transp = transp)
+  qrm_analyse!(spmat, spfct, transp = ntransp)
+  qrm_factorize!(spmat, spfct, transp = ntransp)
+  qrm_solve!(spfct, b, Δx, transp = t)
   qrm_solve!(spfct, Δx, y, transp = 'n')
   #x = A^T y
-  qrm_spmat_mv!(spmat, T(1),  y, T(0), x, transp = transp)
+  qrm_spmat_mv!(spmat, T(1),  y, T(0), x, transp = ntransp)
 end
 
 function qrm_least_squares_semi_normal(spmat :: qrm_spmat{T}, b :: AbstractVecOrMat{T}; transp :: Char = 'n') where T
@@ -148,7 +159,7 @@ function qrm_least_squares_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfc
   qrm_solve!(spfct, z, y, transp = t)
   qrm_solve!(spfct, y, x, transp = 'n')
 
-  qrm_refine!(spmat, spfct, x, z, Δx, y, transp = transp)
+  qrm_refine!(spmat, spfct, x, z, Δx, y)
 end
 
 function qrm_least_squares_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, b :: AbstractMatrix{T}, x :: AbstractMatrix{T}, z :: AbstractMatrix{T}, Δx :: AbstractMatrix{T}, y :: AbstractMatrix{T}; transp :: Char = 'n') where T
@@ -176,5 +187,5 @@ function qrm_least_squares_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfc
   qrm_solve!(spfct, z, y, transp = t)
   qrm_solve!(spfct, y, x, transp = 'n')
 
-  qrm_refine!(spmat, spfct, x, z, Δx, y, transp = transp)
+  qrm_refine!(spmat, spfct, x, z, Δx, y)
 end
diff --git a/test/test_qrm.jl b/test/test_qrm.jl
index a1dc620..e9f74e8 100644
--- a/test/test_qrm.jl
+++ b/test/test_qrm.jl
@@ -198,10 +198,19 @@ end
     for I in (Int32 , Int64)
       A = sprand(T, n, m, 0.3)
       A = convert(SparseMatrixCSC{T,I}, A)
+
+      A_transp = sprand(T, m, n, 0.3)
+      A_transp = convert(SparseMatrixCSC{T,I}, A_transp)
+
       b = rand(T, n)
       B = rand(T, n, p)
+
       spmat = qrm_spmat_init(A)
       spfct = qrm_analyse(spmat, transp=transp)
+
+      spmat_transp = qrm_spmat_init(T)
+      qrm_spmat_init!(spmat_transp, A_transp)
+
       qrm_factorize!(spmat, spfct, transp=transp)
       spfct2 = (T <: Real) ? Transpose(spfct) : Adjoint(spfct)
 
@@ -268,6 +277,10 @@ end
       r = b - A * x
       @test norm(r) ≤ tol
 
+      x = qrm_min_norm_semi_normal(spmat_transp, b, transp = transp)
+      r = b - A_transp' * x
+      @test norm(r) ≤ tol
+
       X = qrm_min_norm(spmat, B)
       R = B - A * X
       @test norm(R) ≤ tol
@@ -276,6 +289,10 @@ end
       R = B - A * X
       @test norm(R) ≤ tol
 
+      X = qrm_min_norm_semi_normal(spmat_transp, B, transp = transp)
+      R = B - A_transp' * X
+      @test norm(R) ≤ tol
+
       qrm_min_norm!(spmat, b, x)
       r = b - A * x
       @test norm(r) ≤ tol

From aa934fc040120213d759102b489010f323613135 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Thu, 19 Dec 2024 16:59:09 +0100
Subject: [PATCH 10/10] add seminormal keyword argument for qrm_min_norm and
 qrm_least_square

---
 src/wrapper/qr_mumps_api.jl | 32 ++++++++++++++++++++------------
 test/test_qrm.jl            | 16 ++++++++++++++++
 2 files changed, 36 insertions(+), 12 deletions(-)

diff --git a/src/wrapper/qr_mumps_api.jl b/src/wrapper/qr_mumps_api.jl
index 2e1f12d..b114917 100644
--- a/src/wrapper/qr_mumps_api.jl
+++ b/src/wrapper/qr_mumps_api.jl
@@ -622,7 +622,9 @@ for (fname, elty) in ((:sqrm_least_squares_c, :Float32   ),
             return nothing
         end
 
-        function qrm_least_squares(spmat :: qrm_spmat{$elty}, b :: Vector{$elty}; transp :: Char='n')
+        function qrm_least_squares(spmat :: qrm_spmat{$elty}, b :: Vector{$elty}; transp :: Char='n', seminormal :: Bool = false)
+            seminormal && return qrm_least_squares_semi_normal(spmat, b, transp = transp)
+
             nrhs = 1
             if transp == 'n'
                 x = zeros($elty, spmat.mat.n)
@@ -635,7 +637,9 @@ for (fname, elty) in ((:sqrm_least_squares_c, :Float32   ),
             return x
         end
 
-        function qrm_least_squares(spmat :: qrm_spmat{$elty}, b :: Matrix{$elty}; transp :: Char='n')
+        function qrm_least_squares(spmat :: qrm_spmat{$elty}, b :: Matrix{$elty}; transp :: Char='n', seminormal :: Bool = false)
+            seminormal && return qrm_least_squares_semi_normal(spmat, b, transp = transp)
+
             nrhs = size(b, 2)
             if transp == 'n'
                 x = zeros($elty, spmat.mat.n, nrhs)
@@ -651,14 +655,14 @@ for (fname, elty) in ((:sqrm_least_squares_c, :Float32   ),
         @inline qrm_least_squares!(spmat :: Transpose{$elty,qrm_spmat{$elty}}, b :: Vector{$elty}, x :: Vector{$elty}) = qrm_least_squares!(spmat.parent, b, x, transp='t')
         @inline qrm_least_squares!(spmat :: Transpose{$elty,qrm_spmat{$elty}}, b :: Matrix{$elty}, x :: Matrix{$elty}) = qrm_least_squares!(spmat.parent, b, x, transp='t')
 
-        @inline qrm_least_squares(spmat  :: Transpose{$elty,qrm_spmat{$elty}}, b :: Vector{$elty}) = qrm_least_squares(spmat.parent, b, transp='t')
-        @inline qrm_least_squares(spmat  :: Transpose{$elty,qrm_spmat{$elty}}, b :: Matrix{$elty}) = qrm_least_squares(spmat.parent, b, transp='t')
+        @inline qrm_least_squares(spmat  :: Transpose{$elty,qrm_spmat{$elty}}, b :: Vector{$elty}; seminormal :: Bool = false) = qrm_least_squares(spmat.parent, b, transp='t', seminormal = false)
+        @inline qrm_least_squares(spmat  :: Transpose{$elty,qrm_spmat{$elty}}, b :: Matrix{$elty}; seminormal :: Bool = false) = qrm_least_squares(spmat.parent, b, transp='t', seminormal = false)
 
         @inline qrm_least_squares!(spmat :: Adjoint{$elty,qrm_spmat{$elty}}, b :: Vector{$elty}, x :: Vector{$elty}) = qrm_least_squares!(spmat.parent, b, x, transp='c')
         @inline qrm_least_squares!(spmat :: Adjoint{$elty,qrm_spmat{$elty}}, b :: Matrix{$elty}, x :: Matrix{$elty}) = qrm_least_squares!(spmat.parent, b, x, transp='c')
 
-        @inline qrm_least_squares(spmat  :: Adjoint{$elty,qrm_spmat{$elty}}, b :: Vector{$elty}) = qrm_least_squares(spmat.parent, b, transp='c')
-        @inline qrm_least_squares(spmat  :: Adjoint{$elty,qrm_spmat{$elty}}, b :: Matrix{$elty}) = qrm_least_squares(spmat.parent, b, transp='c')
+        @inline qrm_least_squares(spmat  :: Adjoint{$elty,qrm_spmat{$elty}}, b :: Vector{$elty}; seminormal :: Bool = false) = qrm_least_squares(spmat.parent, b, transp='c', seminormal = false)
+        @inline qrm_least_squares(spmat  :: Adjoint{$elty,qrm_spmat{$elty}}, b :: Matrix{$elty}; seminormal :: Bool = false) = qrm_least_squares(spmat.parent, b, transp='c', seminormal = false)
 
     end
 end
@@ -692,7 +696,9 @@ for (fname, elty) in ((:sqrm_min_norm_c, :Float32   ),
             return nothing
         end
 
-        function qrm_min_norm(spmat :: qrm_spmat{$elty}, b :: Vector{$elty}; transp :: Char='n')
+        function qrm_min_norm(spmat :: qrm_spmat{$elty}, b :: Vector{$elty}; transp :: Char='n', seminormal :: Bool = false)
+            seminormal && return qrm_min_norm_semi_normal(spmat, b, transp = transp)
+
             nrhs = 1
             if transp == 'n'
                 x = zeros($elty, spmat.mat.n)
@@ -704,7 +710,9 @@ for (fname, elty) in ((:sqrm_min_norm_c, :Float32   ),
             return x
         end
 
-        function qrm_min_norm(spmat :: qrm_spmat{$elty}, b :: Matrix{$elty}; transp :: Char='n')
+        function qrm_min_norm(spmat :: qrm_spmat{$elty}, b :: Matrix{$elty}; transp :: Char='n', seminormal :: Bool = false)
+            seminormal && return qrm_min_norm_semi_normal(spmat, b, transp = transp)
+
             nrhs = size(b, 2)
             if transp == 'n'
                 x = zeros($elty, spmat.mat.n, nrhs)
@@ -719,14 +727,14 @@ for (fname, elty) in ((:sqrm_min_norm_c, :Float32   ),
         @inline qrm_min_norm!(spmat :: Transpose{$elty,qrm_spmat{$elty}}, b :: Vector{$elty}, x :: Vector{$elty}) = qrm_min_norm!(spmat.parent, b, x, transp='t')
         @inline qrm_min_norm!(spmat :: Transpose{$elty,qrm_spmat{$elty}}, b :: Matrix{$elty}, x :: Matrix{$elty}) = qrm_min_norm!(spmat.parent, b, x, transp='t')
 
-        @inline qrm_min_norm(spmat  :: Transpose{$elty,qrm_spmat{$elty}}, b :: Vector{$elty}) = qrm_min_norm(spmat.parent, b, transp='t')
-        @inline qrm_min_norm(spmat  :: Transpose{$elty,qrm_spmat{$elty}}, b :: Matrix{$elty}) = qrm_min_norm(spmat.parent, b, transp='t')
+        @inline qrm_min_norm(spmat  :: Transpose{$elty,qrm_spmat{$elty}}, b :: Vector{$elty}; seminormal :: Bool = false) = qrm_min_norm(spmat.parent, b, transp='t', seminormal = seminormal)
+        @inline qrm_min_norm(spmat  :: Transpose{$elty,qrm_spmat{$elty}}, b :: Matrix{$elty}; seminormal :: Bool = false) = qrm_min_norm(spmat.parent, b, transp='t', seminormal = seminormal)
 
         @inline qrm_min_norm!(spmat :: Adjoint{$elty,qrm_spmat{$elty}}, b :: Vector{$elty}, x :: Vector{$elty}) = qrm_min_norm!(spmat.parent, b, x, transp='c')
         @inline qrm_min_norm!(spmat :: Adjoint{$elty,qrm_spmat{$elty}}, b :: Matrix{$elty}, x :: Matrix{$elty}) = qrm_min_norm!(spmat.parent, b, x, transp='c')
 
-        @inline qrm_min_norm(spmat  :: Adjoint{$elty,qrm_spmat{$elty}}, b :: Vector{$elty}) = qrm_min_norm(spmat.parent, b, transp='c')
-        @inline qrm_min_norm(spmat  :: Adjoint{$elty,qrm_spmat{$elty}}, b :: Matrix{$elty}) = qrm_min_norm(spmat.parent, b, transp='c')
+        @inline qrm_min_norm(spmat  :: Adjoint{$elty,qrm_spmat{$elty}}, b :: Vector{$elty}; seminormal :: Bool = false) = qrm_min_norm(spmat.parent, b, transp='c', seminormal = seminormal)
+        @inline qrm_min_norm(spmat  :: Adjoint{$elty,qrm_spmat{$elty}}, b :: Matrix{$elty}; seminormal :: Bool = false) = qrm_min_norm(spmat.parent, b, transp='c', seminormal = seminormal)
 
     end
 end
diff --git a/test/test_qrm.jl b/test/test_qrm.jl
index e9f74e8..f57ee58 100644
--- a/test/test_qrm.jl
+++ b/test/test_qrm.jl
@@ -90,6 +90,10 @@ p = 5
       r = b - A * x
       @test norm(A' * r) ≤ tol
 
+      x = qrm_least_squares(spmat, b, seminormal = true)
+      r = b - A * x
+      @test norm(A' * r) ≤ tol
+
       x = qrm_least_squares_semi_normal(spmat, b)
       r = b - A * x
       @test norm(A' * r) ≤ tol
@@ -102,6 +106,10 @@ p = 5
       R = B - A * X
       @test norm(A' * R) ≤ tol
 
+      X = qrm_least_squares(spmat, B, seminormal = true)
+      R = B - A * X
+      @test norm(A' * R) ≤ tol
+
       X = qrm_least_squares_semi_normal(spmat, B)
       R = B - A * X
       @test norm(A' * R) ≤ tol
@@ -273,6 +281,10 @@ end
       r = b - A * x
       @test norm(r) ≤ tol
 
+      x = qrm_min_norm(spmat, b, seminormal = true)
+      r = b - A * x
+      @test norm(r) ≤ tol
+
       x = qrm_min_norm_semi_normal(spmat, b)
       r = b - A * x
       @test norm(r) ≤ tol
@@ -284,6 +296,10 @@ end
       X = qrm_min_norm(spmat, B)
       R = B - A * X
       @test norm(R) ≤ tol
+
+      X = qrm_min_norm(spmat, B, seminormal = true)
+      R = B - A * X
+      @test norm(R) ≤ tol
       
       X = qrm_min_norm_semi_normal(spmat, B)
       R = B - A * X