Add semi-normal solvers

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!
 ```

src/QRMumps.jl

@@ -50,7 +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_min_norm_semi_normal,
     qrm_residual_norm!, qrm_residual_norm,
     qrm_residual_orth!, qrm_residual_orth,
     qrm_refine!, qrm_refine, qrm_set, qrm_get
@@ -348,6 +350,48 @@ function qrm_least_squares! end
 """
 function qrm_least_squares end
 
+@doc raw"""
+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
+
+```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(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.
+* `transp`: whether to use A, Aᵀ or Aᴴ. Can be either `'t'`, `'c'` or `'n'`.
+"""
+function qrm_least_squares_semi_normal! end
+
+@doc raw"""
+    x = qrm_least_squares_semi_normal(spmat, b)
+"""
+function qrm_least_squares_semi_normal end
+
 @doc raw"""
     qrm_min_norm!(spmat, b, x; transp='n')
 
@@ -379,6 +423,43 @@ 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(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.
+* `transp`: whether to use A, Aᵀ or Aᴴ. Can be either `'t'`, `'c'` or `'n'`.
+"""
+function qrm_min_norm_semi_normal! end
+
+@doc raw"""
+    x = qrm_min_norm_semi_normal(spmat, b)
+"""
+function qrm_min_norm_semi_normal end
+
 @doc raw"""
     qrm_residual_norm!(spmat, b, x, nrm; transp='n')
 
@@ -427,10 +508,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
 

src/utils.jl

@@ -7,19 +7,185 @@ 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
   @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
+end
+
+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)
+
+  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 = transp)
+  qrm_spmat_mv!(spmat, T(1), y, T(0), Δx, transp = ntransp)
+  @. Δx = z - Δx 
+
+  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}; 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, n)
+  else 
+    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, 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
+
+  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 = 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 = 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
+
+  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)
+
+  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 = ntransp)
+end
+
+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, n)
+  else 
+    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, 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}; 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 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 = 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)
+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
+
+  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)
+
+  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 = ntransp)
+  qrm_solve!(spfct, z, y, transp = t)
+  qrm_solve!(spfct, y, x, transp = 'n')
+
+  qrm_refine!(spmat, spfct, x, z, Δx, y)
+end