Initial wrapper for blasfeo.

.github/dependabot.yml

@@ -0,0 +1,7 @@
+# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates
+version: 2
+updates:
+  - package-ecosystem: "github-actions"
+    directory: "/" # Location of package manifests
+    schedule:
+      interval: "weekly"

.github/workflows/CI.yml

@@ -0,0 +1,44 @@
+name: CI
+on:
+  push:
+    branches:
+      - master
+    tags: ['*']
+  pull_request:
+  workflow_dispatch:
+concurrency:
+  # Skip intermediate builds: always.
+  # Cancel intermediate builds: only if it is a pull request build.
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }}
+jobs:
+  test:
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }}
+    runs-on: ${{ matrix.os }}
+    timeout-minutes: 60
+    permissions: # needed to allow julia-actions/cache to proactively delete old caches that it has created
+      actions: write
+      contents: read
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - '1.10'
+          - '1.12'
+          - 'pre'
+        os:
+          - ubuntu-latest
+          - macos-latest # apple_m1
+          - macos-15-intel # intel based macos
+          - windows-latest # windows
+        arch:
+          - x64
+    steps:
+      - uses: actions/checkout@v6
+      - uses: julia-actions/setup-julia@v2
+        with:
+          version: ${{ matrix.version }}
+          arch: ${{ matrix.arch }}
+      - uses: julia-actions/cache@v2
+      - uses: julia-actions/julia-buildpkg@v1
+      - uses: julia-actions/julia-runtest@v1

.github/workflows/CompatHelper.yml

@@ -0,0 +1,18 @@
+name: CompatHelper
+on:
+  schedule:
+    - cron: 0 0 * * *
+  workflow_dispatch:
+jobs:
+  CompatHelper:
+    runs-on: ubuntu-latest
+    steps:
+      - name: Install CompatHelper
+        run: using Pkg; Pkg.add("CompatHelper")
+        shell: julia --color=yes {0}
+      - name: Run CompatHelper
+        run: using CompatHelper; CompatHelper.main()
+        shell: julia --color=yes {0}
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          COMPATHELPER_PRIV: ${{ secrets.DOCUMENTER_KEY }}

.github/workflows/TagBot.yml

@@ -0,0 +1,19 @@
+name: TagBot
+on:
+  issue_comment:
+    types:
+      - created
+  workflow_dispatch:
+jobs:
+  TagBot:
+    if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot'
+    runs-on: ubuntu-latest
+    steps:
+      - uses: JuliaRegistries/TagBot@v1
+        with:
+          token: ${{ secrets.GITHUB_TOKEN }}
+          # For commits that modify workflow files: SSH key enables tagging, but
+          # releases require manual creation. For full automation of such commits,
+          # use a PAT with `workflow` scope instead of GITHUB_TOKEN.
+          # See: https://github.com/JuliaRegistries/TagBot#commits-that-modify-workflow-files
+          ssh: ${{ secrets.DOCUMENTER_KEY }}

Project.toml

@@ -1,10 +1,16 @@
 name = "BLASFEO"
 uuid = "2d913c56-f67c-11e9-3afb-9119491b9f4c"
-authors = ["Ian McInerney", "Imperial College London"]
 version = "0.1.0"
+authors = ["Anton Pozharskiy <apozharski@gmail.com>", "Ian McInerney <ian.s.mcinerney@ieee.org>"]
+
+[deps]
+CEnum = "fa961155-64e5-5f13-b03f-caf6b980ea82"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+blasfeo_jll = "6b574d4a-bb57-5a4b-b7e6-1c794c903646"
 
 [compat]
-julia = "1"
+blasfeo_jll = "0.1.4"
+julia = "1.10"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

README.md

@@ -1,7 +1,32 @@
 # BLASFEO
 
-[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://juliaembedded.github.io/BLASFEO.jl/stable)
-[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://juliaembedded.github.io/BLASFEO.jl/dev)
-[![Build Status](https://travis-ci.org/JuliaEmbedded/BLASFEO.jl.svg?branch=master)](https://travis-ci.org/JuliaEmbedded/BLASFEO.jl)
-[![Build Status](https://ci.appveyor.com/api/projects/status/github/juliaembedded/BLASFEO.jl?svg=true)](https://ci.appveyor.com/project/imciner2/BLASFEO-jl)
-[![Coveralls](https://coveralls.io/repos/github/JuliaEmbedded/BLASFEO.jl/badge.svg?branch=master)](https://coveralls.io/github/JuliaEmbedded/BLASFEO.jl?branch=master)
+| **License** | **Build Status** | **Coverage** |
+|:-----------:|:-----------------:|:----------------:|
+| [![License: MIT][license-img]][license-url] | [![build-gh][build-gh-img]][build-gh-url] | [![codecov][codecov-img]][codecov-url] |
+
+[license-img]: https://img.shields.io/badge/License-MIT-yellow.svg
+[license-url]: https://github.com/JuliaEmbedded/BLASFEO.jl/blob/master/LICENSE
+[build-gh-img]: https://github.com/JuliaEmbedded/BLASFEO.jl/actions/workflows/CI.yml/badge.svg
+[build-gh-url]: https://github.com/JuliaEmbedded/BLASFEO.jl/actions/workflows/CI.yml
+[codecov-img]: https://codecov.io/gh/JuliaEmbedded/BLASFEO.jl/branch/master/graph/badge.svg?token=MBxH2AAu8Z
+[codecov-url]: https://codecov.io/gh/JuliaEmbedded/BLASFEO.jl
+
+A Julia package wrapping the [`blasfeo`](https://github.com/giaf/blasfeo) linear algebra library. It is a wrapper around the binary package `blasfeo_jll` which builds architecture specific versions of `blasfeo` for the following hosts:
+- `x86_64-[linux|macos|windows]-*` with the following microarchitectures
+  - `x86_64`: A fallback option, which does not have platform specific code.
+  - `avx`: For old processors supporting only the original avx extension.
+  - `avx2`: Modern processors supporting the avx2 extension.
+  - `avx512`: For platforms supporting `avx512`. **Warning**: this only provides meaningful speedups for >= workstation class AMD ZEN5 CPUs. For others it may even degrade performance.
+- `apple_m1-macos-*` for `>=m1` apple architectures.
+
+Currently other `arm` architectures are unsupported but may be in the future if there are usecases.
+
+## Installation
+
+To install `BLASFEO.jl` using
+```julia
+pkg> add https://github.com/JuliaEmbedded/BLASFEO.jl
+```
+
+## Usage
+TODO: add examples

docs/make.jl

@@ -1,14 +1,20 @@
 using Documenter, BLASFEO
 
-makedocs(;
+makedocs(
+    ;
     modules=[BLASFEO],
     format=Documenter.HTML(),
     pages=[
         "Home" => "index.md",
     ],
     repo="https://github.com/JuliaEmbedded/BLASFEO.jl/blob/{commit}{path}#L{line}",
     sitename="BLASFEO.jl",
-    authors="Ian McInerney, Imperial College London",
+    authors="Anton Pozharskiy <apozharski@gmail.com>, Ian McInerney <ian.s.mcinerney@ieee.org>",
+    format=Documenter.HTML(
+        ;
+        edit_link="master",
+        assets=String[],
+    ),
     assets=String[],
 )
 

examples/riccati.jl

@@ -0,0 +1,159 @@
+using BLASFEO
+using BLASFEO.LibBlasfeo
+
+struct OcpStructuredLinearSystem
+    N::Int
+    x0::BlasfeoDvec
+    nx::Vector{Int}
+    nu::Vector{Int}
+    BAbt::Vector{BlasfeoDmat}
+    RSQrq::Vector{BlasfeoDmat}
+    L::Vector{BlasfeoDmat}
+    ux::Vector{BlasfeoDvec}
+    pi::Vector{BlasfeoDvec}
+    work_mat::BlasfeoDmat
+    work_vec::BlasfeoDvec
+end
+
+function OcpStructuredLinearSystem(
+    As::Vector{<:AbstractArray{Float64}},
+    Bs::Vector{<:AbstractArray{Float64}},
+    bs::Vector{<:AbstractArray{Float64}},
+    Qs::Vector{<:AbstractArray{Float64}},
+    qs::Vector{<:AbstractArray{Float64}},
+    Rs::Vector{<:AbstractArray{Float64}},
+    rs::Vector{<:AbstractArray{Float64}},
+    Ss::Vector{<:AbstractArray{Float64}},
+    x0::Vector{Float64}
+    )
+    # Do no size checks :P
+    N = length(Bs)
+    nx = [0; [size(A,1) for A in As]]
+    nu = [[size(B,2) for B in Bs]; 0]
+    ux = [BlasfeoDvec(nx[ii]+nu[ii]+1) for ii in 1:N+1]
+    pi = [BlasfeoDvec(nx[ii]) for ii in 1:N+1]
+
+    # Adjust bs[1]
+    bs[1] += As[1]*x0
+
+    BAbt = Vector{BlasfeoDmat}(undef, N)
+    RSQrq = Vector{BlasfeoDmat}(undef, N+1)
+    L = [BlasfeoDmat(nx[ii]+nu[ii]+1, nx[ii]+nu[ii]) for ii in 1:N+1]
+    BAbt[1] = BlasfeoDmat(permutedims(hcat(Bs[1],bs[1])))
+    RSQrq[1] = BlasfeoDmat([Rs[1]; rs[1]])
+    for ii in 2:N
+        BAbt[ii] = BlasfeoDmat(permutedims(hcat(Bs[ii],As[ii],bs[ii])))
+        RSQ = vcat(
+            hcat(Rs[ii],zeros(nu[ii],nx[ii])),
+            hcat(Ss[ii], Qs[ii]),
+            permutedims(vcat(rs[ii], qs[ii])),
+        )
+        RSQrq[ii] = BlasfeoDmat(RSQ)
+    end
+    RSQrq[N+1] = BlasfeoDmat(vcat(Qs[N+1],transpose(qs[N+1])))
+
+    return OcpStructuredLinearSystem(
+        N,
+        BlasfeoDvec(x0),
+        nx,
+        nu,
+        BAbt,
+        RSQrq,
+        L,
+        ux,
+        pi,
+        BlasfeoDmat(maximum(nu) + maximum(nx) + 1, maximum(nx)),
+        BlasfeoDvec(maximum(nx)),
+    )
+end
+
+
+function simple_model(N;x0=[-1;-1.])
+    As = [[1. 1.; 0. 1.] for i in 1:N]
+    Bs = [[0.; 1;;] for i in 1:N]
+    bs = [[0.1; 0.1] for i in 1:N]
+    Qs = push!([[1. 0.; 0. 1.] for i in 1:N],[10. 0.;0. 10.])
+    qs = push!([[0.1; 0.1] for i in 1:N],[10. ; 10.])
+    Ss = [[0.; 0.] for i in 1:N]
+    Rs = [[1.;;] for i in 1:N]
+    rs = [[0.1] for i in 1:N]
+
+    return As, Bs, bs, Qs, qs, Rs, rs, Ss, x0
+end
+
+"""
+  An example riccati recursion solver via the thin-wrapper interface
+"""
+function riccati(kkt::OcpStructuredLinearSystem)
+    N = kkt.N
+    nx = kkt.nx
+    nu = kkt.nu
+    RSQrq = kkt.RSQrq
+    BAbt = kkt.BAbt
+    ux = kkt.ux
+    pi = kkt.pi
+    work_mat = kkt.work_mat
+    work_vec = kkt.work_vec
+    L = kkt.L
+    nx_end = kkt.nx[N+1]
+    # Factorize last stage Q
+    blasfeo_dpotrf_l_mn(nx_end+1,nx_end, RSQrq[N+1], 0, 0, L[N+1], 0, 0)
+
+    for ii in N:-1:1
+        blasfeo_dtrmm_rlnn(
+            nu[ii]+nx[ii]+1,
+            nx[ii+1],
+            1.0, L[ii+1], nu[ii+1], nu[ii+1],
+            BAbt[ii], 0, 0,
+            work_mat, 0, 0,
+        )
+        blasfeo_dgead(
+			      1, nx[ii+1],
+			      1.0,
+			      L[ii+1], nu[ii+1]+nx[ii+1], nu[ii+1],
+			      work_mat, nu[ii]+nx[ii], 0,
+			  )
+
+        blasfeo_dsyrk_dpotrf_ln_mn(
+			      nu[ii]+nx[ii]+1,
+			      nu[ii]+nx[ii],
+			      nx[ii+1],
+			      work_mat, 0, 0,
+			      work_mat, 0, 0,
+			      RSQrq[ii], 0, 0,
+			      L[ii], 0, 0
+        )
+    end
+
+    blasfeo_drowex(nu[1]+nx[1], -1.0, L[1], nu[1]+nx[1], 0, ux[1], 0)
+	  blasfeo_dtrsv_ltn(nu[1]+nx[1], L[1], 0, 0, ux[1], 0, ux[1], 0)
+	  blasfeo_drowex(nx[2], 1.0, BAbt[1], nu[1]+nx[1], 0, ux[2], nu[2])
+	  blasfeo_dgemv_t(nu[1]+nx[1], nx[2], 1.0, BAbt[1], 0, 0, ux[1], 0, 1.0, ux[2], nu[2], ux[2], nu[2])
+	  blasfeo_dveccp(nx[2], ux[2], nu[2], pi[1], 0)
+	  blasfeo_drowex(nx[2], 1.0, L[2], nu[2]+nx[2], nu[2], work_vec, 0)
+	  blasfeo_dtrmv_ltn(nx[2], L[2], nu[2], nu[2], pi[1], 0, pi[1], 0)
+	  blasfeo_daxpy(nx[2], 1.0, work_vec, 0, pi[1], 0, pi[1], 0)
+	  blasfeo_dtrmv_lnn(nx[2], L[2], nu[2], nu[2], pi[1], 0, pi[1], 0)
+
+    for ii in 2:N
+        blasfeo_drowex(nu[ii], -1.0, L[ii], nu[ii]+nx[ii], 0, ux[ii], 0)
+		    blasfeo_dtrsv_ltn_mn(nu[ii]+nx[ii], nu[ii], L[ii], 0, 0, ux[ii], 0, ux[ii], 0)
+		    blasfeo_drowex(nx[ii+1], 1.0, BAbt[ii], nu[ii]+nx[ii], 0, ux[ii+1], nu[ii+1])
+		    blasfeo_dgemv_t(nu[ii]+nx[ii], nx[ii+1], 1.0, BAbt[ii], 0, 0, ux[ii], 0, 1.0, ux[ii+1], nu[ii+1], ux[ii+1], nu[ii+1])
+		    blasfeo_dveccp(nx[ii+1], ux[ii+1], nu[ii+1], pi[ii], 0)
+		    blasfeo_drowex(nx[ii+1], 1.0, L[ii+1], nu[ii+1]+nx[ii+1], nu[ii+1], work_vec, 0)
+		    blasfeo_dtrmv_ltn(nx[ii+1], L[ii+1], nu[ii+1], nu[ii+1], pi[ii], 0, pi[ii], 0)
+		    blasfeo_daxpy(nx[ii+1], 1.0, work_vec, 0, pi[ii], 0, pi[ii], 0)
+		    blasfeo_dtrmv_lnn(nx[ii+1], L[ii+1], nu[ii+1], nu[ii+1], pi[ii], 0, pi[ii], 0)
+    end
+end
+
+function get_u(kkt::OcpStructuredLinearSystem)
+    return [kkt.ux[ii][1:kkt.nu[ii]] for ii in 1:kkt.N]
+end
+
+function get_x(kkt::OcpStructuredLinearSystem)
+    x = [kkt.ux[ii][kkt.nu[ii]+1:kkt.nu[ii]+kkt.nx[ii]] for ii in 1:kkt.N+1]
+    x[1] = kkt.x0
+    return x
+end