Skip to content

Instantly share code, notes, and snippets.

@unixpickle

unixpickle/LU.swift

Created June 19, 2025 19:56
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save unixpickle/61b4a909ddaee4a2d0a745b53b5e267b to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save unixpickle/61b4a909ddaee4a2d0a745b53b5e267b to your computer and use it in GitHub Desktop.
Download ZIP

Revisions

  1. unixpickle created this gist Jun 19, 2025.
    66 changes: 66 additions & 0 deletions LU.swift
    View file
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,66 @@
    import Accelerate

    private let Epsilon: Double = 1e-8

    private enum LUFactorizationResult {
    case success(LUFactorization)
    case singular(Int)
    }

    private struct LUFactorization {
    let rows: Int
    let cols: Int
    var values: [Double]
    var pivots: [__CLPK_integer]

    mutating func apply(_ mat: inout ColMajorMatrix) {
    var nrhs = __CLPK_integer(mat.cols)
    var ldb = __CLPK_integer(mat.rows)
    var lda = __CLPK_integer(mat.cols)
    var n1 = __CLPK_integer(cols)
    var info: __CLPK_integer = 0
    var order = CChar("N")!
    dgetrs_(&order, &n1, &nrhs, &values, &lda, &pivots, &mat.values, &ldb, &info)
    assert(info == 0)
    }
    }

    private struct ColMajorMatrix {
    let rows: Int
    let cols: Int
    var values: [Double]

    init(rows: Int, cols: Int) {
    self.rows = rows
    self.cols = cols
    self.values = [Double](repeating: 0, count: rows * cols)
    }

    subscript(_ row: Int, _ col: Int) -> Double {
    get {
    assert(row >= 0 && row < rows && col >= 0 && col < cols)
    return values[row + col * rows]
    }
    set {
    assert(row >= 0 && row < rows && col >= 0 && col < cols)
    values[row + col * rows] = newValue
    }
    }

    func lu() -> LUFactorizationResult {
    var newValues = values
    // Compute LU factorization
    var n1 = __CLPK_integer(rows)
    var n2 = __CLPK_integer(cols)
    var lda = n1
    var pivots = [__CLPK_integer](repeating: 0, count: min(rows, cols))
    var info: __CLPK_integer = 0
    dgetrf_(&n1, &n2, &newValues, &lda, &pivots, &info)
    if info != 0 {
    assert(info > 0)
    return .singular(Int(info) - 1)
    } else {
    return .success(LUFactorization(rows: rows, cols: cols, values: newValues, pivots: pivots))
    }
    }
    }