dnrops
/
Surge
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Added `mulAdd(Matrix, Matrix)` (multiply-add) and turned existing addition/subtraction into shallow wrappers around it

This commit is contained in:
Vincent Esche 2019-09-23 14:12:28 +02:00
parent 1c274b3ee7
commit 46b1687a29
2 changed files with 288 additions and 70 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

120
Sources/Surge/Linear Algebra/Matrix.swift
View File

@ -267,23 +267,19 @@ public func ==<T> (lhs: Matrix<T>, rhs: Matrix<T>) -> Bool {
// MARK: - Addition
public func add(_ lhs: Matrix<Float>, _ rhs: Matrix<Float>) -> Matrix<Float> {
precondition(lhs.rows == rhs.rows && lhs.columns == rhs.columns, "Matrix dimensions not compatible with addition")
var result = lhs
result += rhs
return result
}
public func add(_ lhs: Matrix<Double>, _ rhs: Matrix<Double>) -> Matrix<Double> {
precondition(lhs.rows == rhs.rows && lhs.columns == rhs.columns, "Matrix dimensions not compatible with addition")
var result = lhs
var results = rhs
results.grid.withUnsafeMutableBufferPointer { pointer in
cblas_daxpy(Int32(lhs.grid.count), 1.0, lhs.grid, 1, pointer.baseAddress!, 1)
}
result += rhs
return results
return result
}
public func + (lhs: Matrix<Float>, rhs: Matrix<Float>) -> Matrix<Float> {
@ -297,31 +293,11 @@ public func + (lhs: Matrix<Double>, rhs: Matrix<Double>) -> Matrix<Double> {
// MARK: - Addition: In Place
func addInPlace(_ lhs: inout Matrix<Float>, _ rhs: Matrix<Float>) {
precondition(lhs.rows == rhs.rows && lhs.columns == rhs.columns, "Matrix dimensions not compatible with addition")
let gridSize = Int32(lhs.grid.count)
let alpha: Float = 1.0
let stride: Int32 = 1
lhs.grid.withUnsafeMutableBufferPointer { lhsPointer in
rhs.grid.withUnsafeBufferPointer { rhsPointer in
cblas_saxpy(gridSize, alpha, rhsPointer.baseAddress!, stride, lhsPointer.baseAddress!, stride)
}
}
mulAddInPlace(&lhs, rhs, 1.0)
}
func addInPlace(_ lhs: inout Matrix<Double>, _ rhs: Matrix<Double>) {
precondition(lhs.rows == rhs.rows && lhs.columns == rhs.columns, "Matrix dimensions not compatible with addition")
let gridSize = Int32(lhs.grid.count)
let alpha: Double = 1.0
let stride: Int32 = 1
lhs.grid.withUnsafeMutableBufferPointer { lhsPointer in
rhs.grid.withUnsafeBufferPointer { rhsPointer in
cblas_daxpy(gridSize, alpha, rhsPointer.baseAddress!, stride, lhsPointer.baseAddress!, stride)
}
}
mulAddInPlace(&lhs, rhs, 1.0)
}
public func += (lhs: inout Matrix<Float>, rhs: Matrix<Float>) {
@ -335,25 +311,17 @@ public func += (lhs: inout Matrix<Double>, rhs: Matrix<Double>) {
// MARK: - Subtraction
public func sub(_ lhs: Matrix<Float>, _ rhs: Matrix<Float>) -> Matrix<Float> {
precondition(lhs.rows == rhs.rows && lhs.columns == rhs.columns, "Matrix dimensions not compatible with subtraction")
var result = lhs
result -= rhs
var results = rhs
results.grid.withUnsafeMutableBufferPointer { pointer in
catlas_saxpby(Int32(lhs.grid.count), 1.0, lhs.grid, 1, -1, pointer.baseAddress!, 1)
}
return results
return result
}
public func sub(_ lhs: Matrix<Double>, _ rhs: Matrix<Double>) -> Matrix<Double> {
precondition(lhs.rows == rhs.rows && lhs.columns == rhs.columns, "Matrix dimensions not compatible with subtraction")
var result = lhs
result -= rhs
var results = rhs
results.grid.withUnsafeMutableBufferPointer { pointer in
catlas_daxpby(Int32(lhs.grid.count), 1.0, lhs.grid, 1, -1, pointer.baseAddress!, 1)
}
return results
return result
}
public func - (lhs: Matrix<Float>, rhs: Matrix<Float>) -> Matrix<Float> {
@ -364,6 +332,70 @@ public func - (lhs: Matrix<Double>, rhs: Matrix<Double>) -> Matrix<Double> {
return sub(lhs, rhs)
}
// MARK: - Subtraction: In Place
func subInPlace(_ lhs: inout Matrix<Float>, _ rhs: Matrix<Float>) {
mulAddInPlace(&lhs, rhs, -1.0)
}
func subInPlace(_ lhs: inout Matrix<Double>, _ rhs: Matrix<Double>) {
mulAddInPlace(&lhs, rhs, -1.0)
}
public func -= (lhs: inout Matrix<Float>, rhs: Matrix<Float>) {
return subInPlace(&lhs, rhs)
}
public func -= (lhs: inout Matrix<Double>, rhs: Matrix<Double>) {
return subInPlace(&lhs, rhs)
}
// MARK: - Multiply Addition
func mulAdd(_ lhs: Matrix<Float>, _ rhs: Matrix<Float>, _ alpha: Float) -> Matrix<Float> {
var result = lhs
mulAddInPlace(&result, rhs, alpha)
return result
}
func mulAdd(_ lhs: Matrix<Double>, _ rhs: Matrix<Double>, _ alpha: Double) -> Matrix<Double> {
var result = lhs
mulAddInPlace(&result, rhs, alpha)
return result
}
// MARK: - Multiply Addition: In Place
func mulAddInPlace(_ lhs: inout Matrix<Float>, _ rhs: Matrix<Float>, _ alpha: Float) {
precondition(lhs.rows == rhs.rows && lhs.columns == rhs.columns, "Matrix dimensions not compatible with addition")
let gridSize = Int32(lhs.grid.count)
let stride: Int32 = 1
lhs.grid.withUnsafeMutableBufferPointer { lhsPointer in
rhs.grid.withUnsafeBufferPointer { rhsPointer in
cblas_saxpy(gridSize, alpha, rhsPointer.baseAddress!, stride, lhsPointer.baseAddress!, stride)
}
}
}
func mulAddInPlace(_ lhs: inout Matrix<Double>, _ rhs: Matrix<Double>, _ alpha: Double) {
precondition(lhs.rows == rhs.rows && lhs.columns == rhs.columns, "Matrix dimensions not compatible with addition")
let gridSize = Int32(lhs.grid.count)
let stride: Int32 = 1
lhs.grid.withUnsafeMutableBufferPointer { lhsPointer in
rhs.grid.withUnsafeBufferPointer { rhsPointer in
cblas_daxpy(gridSize, alpha, rhsPointer.baseAddress!, stride, lhsPointer.baseAddress!, stride)
}
}
}
// MARK: - Multiplication
public func mul(_ lhs: Matrix<Float>, _ rhs: Matrix<Float>) -> Matrix<Float> {

238
Tests/SurgeTests/MatrixTests.swift
View File

@ -116,8 +116,10 @@ class MatrixTests: XCTestCase {
}
func test_init_identity() {
let actual: Matrix<Double> = Matrix.identity(size: 3)
let expected: Matrix<Double> = [
typealias Scalar = Double
let actual: Matrix<Scalar> = Matrix.identity(size: 3)
let expected: Matrix<Scalar> = [
[1.0, 0.0, 0.0],
[0.0, 1.0, 0.0],
[0.0, 0.0, 1.0],
@ -126,23 +128,25 @@ class MatrixTests: XCTestCase {
}
func test_init_eye() {
let actual_2x3: Matrix<Double> = Matrix.eye(rows: 2, columns: 3)
let expected_2x3: Matrix<Double> = [
typealias Scalar = Double
let actual_2x3: Matrix<Scalar> = Matrix.eye(rows: 2, columns: 3)
let expected_2x3: Matrix<Scalar> = [
[1.0, 0.0, 0.0],
[0.0, 1.0, 0.0],
]
XCTAssertEqual(actual_2x3, expected_2x3)
let actual_3x2: Matrix<Double> = Matrix.eye(rows: 3, columns: 2)
let expected_3x2: Matrix<Double> = [
let actual_3x2: Matrix<Scalar> = Matrix.eye(rows: 3, columns: 2)
let expected_3x2: Matrix<Scalar> = [
[1.0, 0.0],
[0.0, 1.0],
[0.0, 0.0],
]
XCTAssertEqual(actual_3x2, expected_3x2)
let actual_3x3: Matrix<Double> = Matrix.eye(rows: 3, columns: 3)
let expected_3x3: Matrix<Double> = [
let actual_3x3: Matrix<Scalar> = Matrix.eye(rows: 3, columns: 3)
let expected_3x3: Matrix<Scalar> = [
[1.0, 0.0, 0.0],
[0.0, 1.0, 0.0],
[0.0, 0.0, 1.0],
@ -151,23 +155,25 @@ class MatrixTests: XCTestCase {
}
func test_init_diagonal_repeatedValue() {
let actual_2x3: Matrix<Double> = Matrix.diagonal(rows: 2, columns: 3, repeatedValue: 42.0)
let expected_2x3: Matrix<Double> = [
typealias Scalar = Double
let actual_2x3: Matrix<Scalar> = Matrix.diagonal(rows: 2, columns: 3, repeatedValue: 42.0)
let expected_2x3: Matrix<Scalar> = [
[42.0, 0.0, 0.0],
[0.0, 42.0, 0.0],
]
XCTAssertEqual(actual_2x3, expected_2x3)
let actual_3x2: Matrix<Double> = Matrix.diagonal(rows: 3, columns: 2, repeatedValue: 42.0)
let expected_3x2: Matrix<Double> = [
let actual_3x2: Matrix<Scalar> = Matrix.diagonal(rows: 3, columns: 2, repeatedValue: 42.0)
let expected_3x2: Matrix<Scalar> = [
[42.0, 0.0],
[0.0, 42.0],
[0.0, 0.0],
]
XCTAssertEqual(actual_3x2, expected_3x2)
let actual_3x3: Matrix<Double> = Matrix.diagonal(rows: 3, columns: 3, repeatedValue: 42.0)
let expected_3x3: Matrix<Double> = [
let actual_3x3: Matrix<Scalar> = Matrix.diagonal(rows: 3, columns: 3, repeatedValue: 42.0)
let expected_3x3: Matrix<Scalar> = [
[42.0, 0.0, 0.0],
[0.0, 42.0, 0.0],
[0.0, 0.0, 42.0],
@ -176,23 +182,25 @@ class MatrixTests: XCTestCase {
}
func test_init_diagonal_scalars() {
let actual_2x3: Matrix<Double> = Matrix.diagonal(rows: 2, columns: 3, scalars: [1, 2])
let expected_2x3: Matrix<Double> = [
typealias Scalar = Double
let actual_2x3: Matrix<Scalar> = Matrix.diagonal(rows: 2, columns: 3, scalars: [1, 2])
let expected_2x3: Matrix<Scalar> = [
[1.0, 0.0, 0.0],
[0.0, 2.0, 0.0],
]
XCTAssertEqual(actual_2x3, expected_2x3)
let actual_3x2: Matrix<Double> = Matrix.diagonal(rows: 3, columns: 2, scalars: [1, 2])
let expected_3x2: Matrix<Double> = [
let actual_3x2: Matrix<Scalar> = Matrix.diagonal(rows: 3, columns: 2, scalars: [1, 2])
let expected_3x2: Matrix<Scalar> = [
[1.0, 0.0],
[0.0, 2.0],
[0.0, 0.0],
]
XCTAssertEqual(actual_3x2, expected_3x2)
let actual_3x3: Matrix<Double> = Matrix.diagonal(rows: 3, columns: 3, scalars: [1, 2, 3])
let expected_3x3: Matrix<Double> = [
let actual_3x3: Matrix<Scalar> = Matrix.diagonal(rows: 3, columns: 3, scalars: [1, 2, 3])
let expected_3x3: Matrix<Scalar> = [
[1.0, 0.0, 0.0],
[0.0, 2.0, 0.0],
[0.0, 0.0, 3.0],
@ -280,7 +288,7 @@ class MatrixTests: XCTestCase {
]
let actual = matrix + matrix
let expected: Matrix<Float> = [
let expected: Matrix<Scalar> = [
[2, 4, 6, 8],
[10, 12, 14, 16],
[18, 20, 22, 24],
@ -308,6 +316,50 @@ class MatrixTests: XCTestCase {
XCTAssertEqual(actual, expected)
}
// MARK: - Addition: In Place
func test_add_in_place_matrix_matrix_float() {
typealias Scalar = Float
let matrix: Matrix<Scalar> = [
[1, 2, 3, 4],
[5, 6, 7, 8],
[9, 10, 11, 12],
]
var actual = matrix
actual += matrix
let expected: Matrix<Scalar> = [
[2, 4, 6, 8],
[10, 12, 14, 16],
[18, 20, 22, 24],
]
XCTAssertEqual(actual, expected)
}
func test_add_in_place_matrix_matrix_double() {
typealias Scalar = Double
let matrix: Matrix<Scalar> = [
[1, 2, 3, 4],
[5, 6, 7, 8],
[9, 10, 11, 12],
]
var actual = matrix
actual += matrix
let expected: Matrix<Scalar> = [
[2, 4, 6, 8],
[10, 12, 14, 16],
[18, 20, 22, 24],
]
XCTAssertEqual(actual, expected)
}
// MARK: - Subtraction
func test_sub_matrix_matrix_float() {
@ -320,7 +372,7 @@ class MatrixTests: XCTestCase {
]
let actual = matrix - matrix
let expected: Matrix<Float> = [
let expected: Matrix<Scalar> = [
[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0],
@ -348,6 +400,132 @@ class MatrixTests: XCTestCase {
XCTAssertEqual(actual, expected)
}
// MARK: - Subtraction: In Place
func test_sub_in_place_matrix_matrix_float() {
typealias Scalar = Float
let matrix: Matrix<Scalar> = [
[1, 2, 3, 4],
[5, 6, 7, 8],
[9, 10, 11, 12],
]
var actual = matrix
actual -= matrix
let expected: Matrix<Scalar> = [
[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0],
]
XCTAssertEqual(actual, expected)
}
func test_sub_in_place_matrix_matrix_double() {
typealias Scalar = Double
let matrix: Matrix<Scalar> = [
[1, 2, 3, 4],
[5, 6, 7, 8],
[9, 10, 11, 12],
]
var actual = matrix
actual -= matrix
let expected: Matrix<Scalar> = [
[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0],
]
XCTAssertEqual(actual, expected)
}
// MARK: - Multiply Addition
func test_mul_add_matrix_matrix_scalar_float() {
typealias Scalar = Float
let matrix: Matrix<Scalar> = [
[1, 2, 3, 4],
[5, 6, 7, 8],
[9, 10, 11, 12],
]
let actual = matrix + matrix
let expected: Matrix<Scalar> = [
[2, 4, 6, 8],
[10, 12, 14, 16],
[18, 20, 22, 24],
]
XCTAssertEqual(actual, expected)
}
func test_mul_add_matrix_matrix_scalar_double() {
typealias Scalar = Double
let matrix: Matrix<Scalar> = [
[1, 2, 3, 4],
[5, 6, 7, 8],
[9, 10, 11, 12],
]
let actual = mulAdd(matrix, matrix, 2.0)
let expected: Matrix<Scalar> = [
[3, 6, 9, 12],
[15, 18, 21, 24],
[27, 30, 33, 36],
]
XCTAssertEqual(actual, expected)
}
// MARK: - Multiply Addition: In Place
func test_mul_add_in_place_matrix_matrix_scalar_float() {
typealias Scalar = Float
let matrix: Matrix<Scalar> = [
[1, 2, 3, 4],
[5, 6, 7, 8],
[9, 10, 11, 12],
]
let actual = mulAdd(matrix, matrix, 2.0)
let expected: Matrix<Scalar> = [
[3, 6, 9, 12],
[15, 18, 21, 24],
[27, 30, 33, 36],
]
XCTAssertEqual(actual, expected)
}
func test_mul_add_in_place_matrix_matrix_scalar_double() {
typealias Scalar = Double
let matrix: Matrix<Scalar> = [
[1, 2, 3, 4],
[5, 6, 7, 8],
[9, 10, 11, 12],
]
var actual = matrix
actual += matrix
let expected: Matrix<Scalar> = [
[2, 4, 6, 8],
[10, 12, 14, 16],
[18, 20, 22, 24],
]
XCTAssertEqual(actual, expected)
}
// MARK: - Multiplication
func test_mul_matrix_matrix_float() {
@ -829,7 +1007,9 @@ class MatrixTests: XCTestCase {
}
func test_eigen_decomposition_trivial() throws {
try eigen_decomposition_trivial_generic(defaultAccuracy: 1e-11) { (m: Matrix<Double>) throws -> MatrixEigenDecompositionResult<Double> in
typealias Scalar = Double
try eigen_decomposition_trivial_generic(defaultAccuracy: 1e-11) { (m: Matrix<Scalar>) throws -> MatrixEigenDecompositionResult<Double> in
try eigenDecompose(m)
}
try eigen_decomposition_trivial_generic(defaultAccuracy: 1e-8) { (m: Matrix<Float>) throws -> MatrixEigenDecompositionResult<Float> in
@ -856,7 +1036,9 @@ class MatrixTests: XCTestCase {
}
func test_eigen_decomposition_complex_results_numpy_example() throws {
try eigen_decomposition_complex_results_numpy_example_generic(defaultAccuracy: 1e-11) { (m: Matrix<Double>) throws -> MatrixEigenDecompositionResult<Double> in
typealias Scalar = Double
try eigen_decomposition_complex_results_numpy_example_generic(defaultAccuracy: 1e-11) { (m: Matrix<Scalar>) throws -> MatrixEigenDecompositionResult<Double> in
try eigenDecompose(m)
}
try eigen_decomposition_complex_results_numpy_example_generic(defaultAccuracy: 1e-8) { (m: Matrix<Float>) throws -> MatrixEigenDecompositionResult<Float> in
@ -905,7 +1087,9 @@ class MatrixTests: XCTestCase {
}
func test_eigen_decomposition_complex_results_dgeev_example() throws {
try eigen_decomposition_complex_results_dgeev_example_generic(defaultAccuracy: 1e-11) { (m: Matrix<Double>) throws -> MatrixEigenDecompositionResult<Double> in
typealias Scalar = Double
try eigen_decomposition_complex_results_dgeev_example_generic(defaultAccuracy: 1e-11) { (m: Matrix<Scalar>) throws -> MatrixEigenDecompositionResult<Double> in
try eigenDecompose(m)
}
try eigen_decomposition_complex_results_dgeev_example_generic(defaultAccuracy: 1e-8) { (m: Matrix<Float>) throws -> MatrixEigenDecompositionResult<Float> in
@ -914,7 +1098,9 @@ class MatrixTests: XCTestCase {
}
func test_eigen_decompose_throws_when_not_square() throws {
let matrix = Matrix<Double>([
typealias Scalar = Double
let matrix = Matrix<Scalar>([
[1, 0, 0],
[0, 2, 0],
])