Interface qrm_spfct_unmqr

src/QRMumps.jl

@@ -42,7 +42,7 @@ export qrm_spmat, qrm_spfct,
     qrm_update!, qrm_factorize!,
     qrm_solve!, qrm_solve,
     qrm_apply!, qrm_apply,
-    qrm_spfct_get_rp, qrm_spfct_get_cp, qrm_spfct_get_r,
+    qrm_spfct_get_rp, qrm_spfct_get_cp, qrm_spfct_get_r, qrm_spfct_unmqr,
     qrm_spmat_mv!, mul!, qrm_spmat_nrm,
     qrm_vecnrm!, qrm_vecnrm,
     qrm_spbackslash!, qrm_spbackslash, \,
@@ -81,7 +81,7 @@ function qrm_finalize end
     qrm_spmat_init!(spmat, A; sym=false)
     spmat = qrm_spmat_init!(spmat, m, n, rows, cols, vals; sym=false)
 
-This routine initializes a **qrm_spmat** type data structure from a **sparseMatrixCSC**.
+This routine initializes a **qrm_spmat** type data structure from a **SparseMatrixCSC**.
 
 #### Input Arguments :
 
@@ -151,7 +151,7 @@ function qrm_spfct_destroy! end
     qrm_update!(spmat, A)
     qrm_update!(spmat, vals)
 
-This routine updates a **qrm_spmat** type data structure from a **sparseMatrixCSC**.
+This routine updates a **qrm_spmat** type data structure from a **SparseMatrixCSC**.
 `spmat` and `A` must have the same sparsity pattern.
 
 #### Input Arguments :
@@ -167,6 +167,20 @@ In the second form,
 """
 function qrm_update! end
 
+@doc raw"""
+    qrm_spfct_unmqr(spfct :: qrm_spfct, transp, b)
+
+#### Input Arguments :
+
+This routine computes `Qb`, `Qᵀb` or `Qᴴb` in place and overwrites b.
+It can only be executed once the factorization is done.
+
+* `spfct`: the sparse factorization object of type `qrm_spfct`.
+* `transp`: whether the input matrix should be transposed or not. Can be either `'t'`, `'c'` or `'n'`.
+* `b`: the vector or matrix on which we apply the orthogonal matrix.
+"""
+function qrm_spfct_unmqr end
+
 @doc raw"""
     qrm_analyse!(spmat, spfct; transp='n')
 

src/wrapper/qr_mumps_api.jl

@@ -1,3 +1,17 @@
+for (fname, elty) in ((:sqrm_spfct_unmqr_c, :Float32   ),
+                      (:dqrm_spfct_unmqr_c, :Float64   ),
+                      (:cqrm_spfct_unmqr_c, :ComplexF32),
+                      (:zqrm_spfct_unmqr_c, :ComplexF64))
+    @eval begin
+        function qrm_spfct_unmqr(spfct :: qrm_spfct{$elty}, trans::Char, b::VecOrMat{$elty})
+            nrhs = size(b, 2)
+            err = $fname(spfct, trans, b, nrhs)
+            qrm_check(err)
+            return b
+        end
+    end
+end
+
 for (fname, elty) in ((:sqrm_spfct_get_rp_c, :Float32   ),
                       (:dqrm_spfct_get_rp_c, :Float64   ),
                       (:cqrm_spfct_get_rp_c, :ComplexF32),

test/test_qrm.jl

@@ -1,17 +1,17 @@
 Random.seed!(1234)
 m = 200
 n = 100
-p = 5
+p = 2
 
 @testset "least-squares problems" begin
   for T in (Float32, Float64, ComplexF32, ComplexF64)
 
     tol = (real(T) == Float32) ? 1e-3 : 1e-12
     transp = (T <: Real) ? 't' : 'c'
 
-    for I in (Int32 , Int64)
+    for INT in (Int32 , Int64)
       A = sprand(T, m, n, 0.3)
-      A = convert(SparseMatrixCSC{T,I}, A)
+      A = convert(SparseMatrixCSC{T,INT}, A)
       b = rand(T, m)
       B = rand(T, m, p)
 
@@ -40,6 +40,15 @@ p = 5
       R = B - A * X
       @test norm(A' * R) ≤ tol
 
+      if VERSION ≥ v"1.10"
+        Q = Matrix{T}(I, m, m)
+        qrm_spfct_unmqr(spfct, 'n', Q)
+        @test Q * Q' ≈ I
+        Qᴴ = Matrix{T}(I, m, m)
+        qrm_spfct_unmqr(spfct, transp, Qᴴ)
+        @test Qᴴ * Qᴴ' ≈ I
+      end
+
       spmat = qrm_spmat_init(T)
       qrm_spmat_init!(spmat, A)
 
@@ -158,9 +167,9 @@ end
     tol = (real(T) == Float32) ? 1e-3 : 1e-12
     transp = (T <: Real) ? 't' : 'c'
 
-    for I in (Int32 , Int64)
+    for INT in (Int32 , Int64)
       A = sprand(T, n, m, 0.3)
-      A = convert(SparseMatrixCSC{T,I}, A)
+      A = convert(SparseMatrixCSC{T,INT}, A)
       b = rand(T, n)
       B = rand(T, n, p)
       spmat = qrm_spmat_init(A)
@@ -188,6 +197,15 @@ end
       R = B - A * X
       @test norm(R) ≤ tol
 
+      if VERSION ≥ v"1.10"
+        Q = Matrix{T}(I, m, m)
+        qrm_spfct_unmqr(spfct, 'n', Q)
+        @test Q * Q' ≈ I
+        Qᴴ = Matrix{T}(I, m, m)
+        qrm_spfct_unmqr(spfct, transp, Qᴴ)
+        @test Qᴴ * Qᴴ' ≈ I
+      end
+
       spmat = qrm_spmat_init(T)
       qrm_spmat_init!(spmat, A)
 
@@ -297,9 +315,9 @@ end
     transp = (T <: Real) ? 't' : 'c'
     Id = Diagonal(ones(T, n))
 
-    for I in (Int32 , Int64)
+    for INT in (Int32 , Int64)
       A = sprand(T, n, n, 0.01)
-      A = convert(SparseMatrixCSC{T,I}, A)
+      A = convert(SparseMatrixCSC{T,INT}, A)
       A = A * A' + Id
       A = (T <: Real) ? Symmetric(tril(A), :L) : Hermitian(tril(A), :L)
       b = rand(T, n)
@@ -331,7 +349,7 @@ end
       @test norm(R) ≤ tol
 
       A = sprand(T, n, n, 0.01)
-      A = convert(SparseMatrixCSC{T,I}, A)
+      A = convert(SparseMatrixCSC{T,INT}, A)
       A = A * A' + Id
       A = (T <: Real) ? Symmetric(triu(A), :U) : Hermitian(triu(A), :U)
       x = zeros(T, n)