Interface MC26

src/HSL.jl

@@ -46,6 +46,7 @@ include("hsl_mc64.jl")
 include("kb07.jl")
 include("mc19.jl")
 include("mc21.jl")
+include("mc26.jl")
 include("mc29.jl")
 include("mc30.jl")
 include("mc77.jl")

src/mc26.jl

@@ -0,0 +1,60 @@
+export mc26
+
+"""
+    jptr, iptr, ata = mc26(A::SparseMatrixCSC{Float32,Cint})
+    jptr, iptr, ata = mc26(A::SparseMatrixCSC{Float64,Cint})
+
+Compute the structure and optionally the values of `AᵀA`.
+
+Returns:
+- `jptr`: column pointer array for the compressed sparse column (CSC) format of `AᵀA`.
+- `iptr`: row indices corresponding to each nonzero entry in `AᵀA`.
+- `ata`: numerical values of `AᵀA` corresponding to the `(iptr, jptr)` structure.
+"""
+function mc26 end
+
+for (mc26ir, mc26ar, T) in ((:mc26i , :mc26a , :Float32),
+                            (:mc26id, :mc26ad, :Float64))
+  @eval begin
+    function mc26(A::SparseMatrixCSC{$T,Cint})
+      icntl = zeros(Cint, 5)
+      info = zeros(Cint, 5)
+      $mc26ir(icntl)  # Set default control parameters
+
+      m, n = size(A)
+      itype = Cint(1)  # CSC format
+      nz = nnz(A)
+      a = A.nzval
+      ind = A.rowval
+
+      lata = max(1, nz + n)  # initial guess
+      iata = Vector{Cint}(undef, lata)
+      ata  = Vector{$T}(undef, lata)
+      jptr = Vector{Cint}(undef, n + 1)
+      iptr = Vector{Cint}(undef, lata)
+      iw = zeros(Cint, max(1, 2 * n))
+      nzata = Ref{Cint}(0)
+      success = false
+
+      while !success
+        $mc26ar(m, n, itype, nz, a, ind, lata, iata, ata, jptr, iptr, nzata, iw, icntl, info)
+
+        if info[1] == -4
+          lata = max(lata + n, iata[1])
+          resize!(iata, lata)
+          resize!(ata,  lata)
+          resize!(iptr, lata)
+        elseif info[1] == 0
+          success = true
+        else
+          error("mc26a failed with info[1] = $(info[1])")
+        end
+      end
+
+      resize!(ata, nzata[])
+      resize!(iptr, nzata[])
+
+      return jptr, iptr, ata
+    end
+  end
+end

test/runtests.jl

@@ -13,6 +13,7 @@ if LIBHSL_isfunctional()
   include("test_kb07.jl")
   include("test_mc19.jl")
   include("test_mc21.jl")
+  # include("test_mc26.jl")
   include("test_mc29.jl")
   include("test_mc30.jl")
   include("test_mc77.jl")

test/test_mc26.jl

@@ -0,0 +1,19 @@
+@testset "mc26" begin
+    # Define the matrix A:
+    # [1.0 4.0]
+    # [0.0 3.0]
+    # [5.0 0.0]
+    rowval = Cint[1, 2, 3, 1]  # rows of nonzeros (Fortran-style indexing)
+    colptr = Cint[1, 3, 5]     # column pointer (CSC)
+    nzval  = [1.0, 3.0, 5.0, 4.0]  # nonzero values
+
+    @testset "Precision $T" for T in (Float32, Float64)
+        A = SparseMatrixCSC{T,Cint}(3, 2, colptr, rowval, T.(nzval))
+
+        jptr, iptr, ata = mc26(A)
+
+        @test jptr == Cint[1, 3, 4]   # two columns in AᵀA, jptr[3] = 4 means 3 nnz
+        @test iptr == Cint[1, 2, 2]   # row indices for columns
+        @test ata ≈ T[26.0, 4.0, 25.0] # AᵀA = [26 4; 4 25], stored by columns
+    end
+end