Calculate dot product using Real type from Numerics package

Related Projects Swift Numerics
1

I have a Vector struct that I use for linear algebra calculations where the generic values conform to the Real type provided by the Numerics package. The example shown below is a simplified version of the struct. I can calculate the dot product of two vectors using the cblas_sdot() and cblas_ddot() functions from Accelerate which are only for Float and Double types respectively.

import Accelerate
import Numerics

struct Vector<T: Real> {

    let size: Int
    var values: [T]

    init(_ values: [T]) {
        self.size = values.count
        self.values = values
    }
}

func dotProduct(_ x: Vector<Float>, _ y: Vector<Float>) -> Float {
    precondition(x.size == y.size, "Vectors must be same size")
    var dot: Float = 0

    x.values.withUnsafeBufferPointer { xPointer in
        y.values.withUnsafeBufferPointer { yPointer in
            dot = cblas_sdot(
                Int32(x.size), .init(xPointer.baseAddress), 1, .init(yPointer.baseAddress), 1
            )
        }
    }

    return dot
}

I would like to have one generic dot product function instead of defining separate functions for Float and Double. I attempted this with the function shown below but I don't know how to convey the type information to the vector structs. I also can't pass a Real value type to Accelerate so I have to switch on the type to use the appropriate function.

func dotProduct<T: Real>(_ x: Vector<T>, _ y: Vector<T>) -> T {
    precondition(x.size == y.size, "Vectors must be same size")

    switch T.self {
    case is Float.Type:
        var dot: Float = 0

        x.values.withUnsafeBufferPointer { xPointer in
            y.values.withUnsafeBufferPointer { yPointer in
                dot = cblas_sdot(
                    Int32(x.size), .init(xPointer.baseAddress), 1, .init(yPointer.baseAddress), 1
                )
            }
        }

        return dot
    case is Double.Type:
        var dot: Double = 0

        x.values.withUnsafeBufferPointer { xPointer in
            y.values.withUnsafeBufferPointer { yPointer in
                dot = cblas_ddot(
                    Int32(x.size), .init(xPointer.baseAddress), 1, .init(yPointer.baseAddress), 1
                )
            }
        }

        return dot
    default:
        fatalError("Only Float and Double are supported")
    }
}

Does the Numerics package have any features that would help me write generic functions that use the Accelerate framework? I know I have asked similar questions in the forum on this topic but those questions didn't involve the use of the Numerics package.

2

You could declare a new protocol and conform Float and Double to this protocol. Then the switch over the type could call a protocol requirement instead:

protocol Scalar: Real {
  static func dot(_ x: Vector<Self>, _ y: Vector<Self>) -> Self
}

(However Real is maybe too strict here since the dot product is defined on vectors of complex numbers as well. Also Vector is a type I made up here and not the Vector being pitched right now.)

You could also add a default implementation in a protocol extension to compute the dot product in an unaccelerated way, and conform other types to your protocol as needed.

3

Based on your suggestion, the following code uses a Scalar protocol to extend Int, Float, and Double.

import Accelerate
import Numerics

struct Vector<T> {

    let size: Int
    var values: [T]

    init(_ values: [T]) {
        self.size = values.count
        self.values = values
    }
}

protocol Scalar {
    static func dot(_ x: Vector<Self>, _ y: Vector<Self>) -> Self
}

extension Int: Scalar {

    static func dot(_ x: Vector<Int>, _ y: Vector<Int>) -> Int {
        precondition(x.size == y.size, "Vectors must be same size")
        var result = 0

        for (a, b) in zip(x.values, y.values) {
            result += a * b
        }

        return result
    }
}

extension Float: Scalar {

    static func dot(_ x: Vector<Float>, _ y: Vector<Float>) -> Float {
        precondition(x.size == y.size, "Vectors must be same size")
        var result: Float = 0

        x.values.withUnsafeBufferPointer { xPointer in
            y.values.withUnsafeBufferPointer { yPointer in
                result = cblas_sdot(
                    Int32(x.size), .init(xPointer.baseAddress), 1, .init(yPointer.baseAddress), 1
                )
            }
        }

        return result
    }
}


extension Double: Scalar {

    static func dot(_ x: Vector<Double>, _ y: Vector<Double>) -> Double {
        precondition(x.size == y.size, "Vectors must be same size")
        var result = 0.0

        x.values.withUnsafeBufferPointer { xPointer in
            y.values.withUnsafeBufferPointer { yPointer in
                result = cblas_ddot(
                    Int32(x.size), .init(xPointer.baseAddress), 1, .init(yPointer.baseAddress), 1
                )
            }
        }

        return result
    }
}

func dotProduct<T: Scalar>(_ x: Vector<T>, _ y: Vector<T>) -> T {
    T.dot(x, y)
}

This allows me to use the Vector struct as follows:

let a = Vector<Float>([1, 2, 3, 4, 5])
let b = Vector<Float>([3, 4, 5, 6, 7])
let c = dotProduct(a, b)
print(c)

let s = Vector([1, 2, 3, 4, 5.0])
let t = Vector([3, 4, 5, 6, 17.0])
let u = dotProduct(s, t)
print(u)

let x = Vector([1, 2, 3, 4, 5])
let y = Vector([3, 4, 5, 6, 7])
let z = dotProduct(x, y)
print(z)

But I'm not using Real from the Numerics package in any of this code. I thought the purpose of the Numerics package was to help write more generic numerical code but in this case it doesn't seem necessary. Is there something I'm missing to properly use Numerics package for this type of problem?

You also mentioned using a default implementation in a protocol extension but how is that done with a protocol? I tried the following but I don't see how to access the function. I can't call the function with Scalar.dot(x, y).

extension Scalar {

    static func dot(_ x: Vector<Int>, _ y: Vector<Int>) -> Int {
        precondition(x.size == y.size, "Vectors must be same size")
        var result = 0

        for (a, b) in zip(x.values, y.values) {
            result += a * b
        }

        return result
    }
}
4

You don't need to use Real to call the Accelerate functions because there you're operating on concrete Float and Double types.

However,

If your protocol inherits from Real, then you can use math operators in the default implementation of dot(), to implement it for Real-conforming types other than Float and Double. (Perhaps there aren't any though?)

Don't forget to replace Int with Self, since this is a protocol extension.

5

If I replace Int with Self as done in the following code

extension Scalar {

    static func dot(_ x: Vector<Self>, _ y: Vector<Self>) -> Self {
        precondition(x.size == y.size, "Vectors must be same size")
        var result = 0

        for (a, b) in zip(x.values, y.values) {
            result += a * b
        }

        return result
    }
}

I get warnings that the binary operator += can't be applied and the return type can't be converted to Self when I use integer values such as

let x = Vector([1, 2, 3, 4, 5])
let y = Vector([3, 4, 5, 6, 7])
let z = dotProduct(x, y)

The only way I can get this to work is if I conform Int to Scalar like I did for Float and Double. So I don't see how a default implementation can be written in the protocol extension to handle other types.

6

The inferred type of result is always Int here which is not what you want. Instead try this:

let result = Self.zero

Int cannot conform to Real because Real inherits from FloatingPoint. However Numeric suffices for your purposes, to which Int conforms. You don't need any of the additional requirements of Real here.

The idea is that you'll still need to actually conform Int or whatever to your protocol, but you'll be able to use the default unless a specialized implementation is needed:

extension Int: Scalar {}
7

Using Self.zero does not work. How do I get the value from Self.zero? Here is my latest code:

import Accelerate
import Numerics

struct Vector<T> {
    let size: Int
    var values: [T]

    init(_ values: [T]) {
        self.size = values.count
        self.values = values
    }
}

protocol Scalar {
    static var zero: Self { get }
    static func dot(_ x: Vector<Self>, _ y: Vector<Self>) -> Self
}

extension Scalar {

    static func dot(_ x: Vector<Self>, _ y: Vector<Self>) -> Self {
        precondition(x.size == y.size, "Vectors must be same size")

        var result = Self.zero    // This does not get a zero value

        for (a, b) in zip(x.values, y.values) {
            result += a * b       // This fails because binary operator += cannot be applied
        }

        return result
    }
}

extension Int: Scalar { }

extension Float: Scalar {

    static func dot(_ x: Vector<Float>, _ y: Vector<Float>) -> Float {
        precondition(x.size == y.size, "Vectors must be same size")
        var result: Float = 0

        x.values.withUnsafeBufferPointer { xPointer in
            y.values.withUnsafeBufferPointer { yPointer in
                result = cblas_sdot(
                    Int32(x.size), .init(xPointer.baseAddress), 1, .init(yPointer.baseAddress), 1
                )
            }
        }

        return result
    }
}


extension Double: Scalar {

    static func dot(_ x: Vector<Double>, _ y: Vector<Double>) -> Double {
        precondition(x.size == y.size, "Vectors must be same size")
        var result = 0.0

        x.values.withUnsafeBufferPointer { xPointer in
            y.values.withUnsafeBufferPointer { yPointer in
                result = cblas_ddot(
                    Int32(x.size), .init(xPointer.baseAddress), 1, .init(yPointer.baseAddress), 1
                )
            }
        }

        return result
    }
}

func dotProduct<T: Scalar>(_ x: Vector<T>, _ y: Vector<T>) -> T {
    T.dot(x, y)
}
8

zero is a requirement of the AdditiveArithmetic protocol:

static var zero: Self { get }

You also need * to implement a generic dot product, so you want your protocol to inherit from Numeric, which inherits AdditiveArithmetic:

protocol Scalar: Numeric {...}
1 Like
9

Ok, I gave up on using Real and created a Scalar protocol that conforms to Numeric as you suggested. Everything is working fine now. Thank you very much for the help. See below for the working code.

import Accelerate

// --- Protocol

protocol Scalar: Numeric {
    static func add(_ a: Vector<Self>, _ b: Vector<Self>) -> Vector<Self>
    static func subtract(_ a: Vector<Self>, _ b: Vector<Self>) -> Vector<Self>
    static func dot(_ a: Vector<Self>, _ b: Vector<Self>) -> Self
}

// --- Protocol extension

extension Scalar {
    static func add(_ a: Vector<Self>, _ b: Vector<Self>) -> Vector<Self> {
        var vec = Vector(like: a)
        for i in 0..<a.size {
            vec[i] = a[i] + b[i]
        }
        return vec
    }

    static func subtract(_ a: Vector<Self>, _ b: Vector<Self>) -> Vector<Self> {
        var vec = Vector(like: a)
        for i in 0..<a.size {
            vec[i] = a[i] - b[i]
        }
        return vec
    }

    static func dot(_ x: Vector<Self>, _ y: Vector<Self>) -> Self {
        var result = Self.zero
        for (a, b) in zip(x.values, y.values) {
            result += a * b
        }
        return result
    }
}

// --- Type extensions

extension Int: Scalar { }

extension Int16: Scalar { }

extension Float: Scalar {
    static func add(_ a: Vector<Float>, _ b: Vector<Float>) -> Vector<Float> {
        var vec = Vector(like: a)
        vDSP.add(a.values, b.values, result: &vec.values)
        return vec
    }

    static func subtract(_ a: Vector<Float>, _ b: Vector<Float>) -> Vector<Float> {
        var vec = Vector(like: a)
        vDSP.subtract(a.values, b.values, result: &vec.values)
        return vec
    }

    static func dot(_ x: Vector<Float>, _ y: Vector<Float>) -> Float {
        var result: Float = 0
        x.values.withUnsafeBufferPointer { xPointer in
            y.values.withUnsafeBufferPointer { yPointer in
                result = cblas_sdot(
                    Int32(x.size), .init(xPointer.baseAddress), 1, .init(yPointer.baseAddress), 1
                )
            }
        }
        return result
    }
}

extension Double: Scalar {
    static func add(_ a: Vector<Double>, _ b: Vector<Double>) -> Vector<Double> {
        var vec = Vector(like: a)
        vDSP.add(a.values, b.values, result: &vec.values)
        return vec
    }

    static func subtract(_ a: Vector<Double>, _ b: Vector<Double>) -> Vector<Double> {
        var vec = Vector(like: a)
        vDSP.subtract(a.values, b.values, result: &vec.values)
        return vec
    }

    static func dot(_ x: Vector<Double>, _ y: Vector<Double>) -> Double {
        var result: Double = 0
        x.values.withUnsafeBufferPointer { xPointer in
            y.values.withUnsafeBufferPointer { yPointer in
                result = cblas_ddot(
                    Int32(x.size), .init(xPointer.baseAddress), 1, .init(yPointer.baseAddress), 1
                )
            }
        }
        return result
    }
}

// --- Vector numeric type

struct Vector<T: Scalar> {

    let size: Int
    var values: [T]

    init(_ values: [T]) {
        self.size = values.count
        self.values = values
    }

    init(like vector: Self) {
        self.size = vector.size
        self.values = vector.values
    }

    subscript(item: Int) -> T {
        get { return self.values[item] }
        set { self.values[item] = newValue }
    }
}

// --- Vector addition and subtraction

extension Vector: AdditiveArithmetic {

    static var zero: Vector { Vector([0, 0, 0, 0, 0]) }

    static func + (lhs: Vector, rhs: Vector) -> Vector {
        T.add(lhs, rhs)
    }

    static func - (lhs: Vector, rhs: Vector) -> Vector {
        T.subtract(lhs, rhs)
    }
}

// --- Vector dot product

func dotProduct<T: Scalar>(_ x: Vector<T>, _ y: Vector<T>) -> T {
    T.dot(x, y)
}

Where should I put the precondition statements? Should they go in the top-level functions defined in the vector extension, in the functions defined in the value type extensions, or both?

In the Swift Numerics package, I noticed a lot of the code using the @inlinable attribute. If I put the above code into a package, where would I apply this attribute?

1 Like