/*
* Copyright (c) Huawei Technologies Co., Ltd. 2024-2024. All rights reserved.
*/
package matrix4cj
import std.math.*
import std.convert.*
import std.random.Random as RandomGenerator
/**
<P>
The CJ Matrix Class provides the fundamental operations of numerical
linear algebra. Various constructors create Matrices from two dimensional
arrays of double precision floating point numbers. Various "gets" and
"sets" provide access to submatrices and matrix elements. Several methods
implement basic matrix arithmetic, including matrix addition and
multiplication, matrix norms, and element-by-element array operations.
Methods for reading and printing matrices are also included. All the
operations in this version of the Matrix Class involve real matrices.
Complex matrices may be handled in a future version.
<P>
Five fundamental matrix decompositions, which consist of pairs or triples
of matrices, permutation vectors, and the like, produce results in five
decomposition classes. These decompositions are accessed by the Matrix
class to compute solutions of simultaneous linear equations, determinants,
inverses and other matrix functions. The five decompositions are:
<P><UL>
<LI>Cholesky Decomposition of symmetric, positive definite matrices.
<LI>LU Decomposition of rectangular matrices.
<LI>QR Decomposition of rectangular matrices.
<LI>Singular Value Decomposition of rectangular matrices.
<LI>Eigenvalue Decomposition of both symmetric and nonsymmetric square matrices.
</UL>
<DL>
<DT><B>Example of use:</B></DT>
<P>
<DD>Solve a linear system A x = b and compute the residual norm, ||b - A x||.
<P><PRE>
let vals: Array<Array<Float64>> = [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 10.0]]
let A: Matrix = Matrix(vals)
let b: Matrix = Matrix.random(3, 1)
let x: Matrix = A.solve(b)
let r: Matrix = A.times(x) - b
let rnorm: Float64 = r.normInf()
</PRE></DD>
</DL>
*/
public interface IMatrixArithmetic<Base> where Base <: Number<Float64> {
operator func *(other: Matrix): Matrix
operator func +(other: Matrix): Matrix
operator func -(other: Matrix): Matrix
operator func /(other: Matrix): Matrix
}
public enum StorageOrder {
/** Column-major order */
ColumnMajor |
/** Row-major order */
RowMajor
public operator func == (other: StorageOrder): Bool {
match ((this, other)) {
case (ColumnMajor, ColumnMajor) |
(RowMajor, RowMajor) => true
case _ => false
}
}
public operator func != (other: StorageOrder): Bool {
return !(this == other)
}
}
public enum Axis {
/** Row axis */
Row |
/** Column axis */
Column
public operator func == (other: Axis): Bool {
match ((this, other)) {
case (Row, Row) |
(Column, Column) => true
case _ => false
}
}
public operator func != (other: Axis): Bool {
return !(this == other)
}
}
public class Matrix <: IElementWiseOperation<Matrix> & Equatable<Matrix> {
/** Array for internal storage of elements.
* @serial internal array storage.
*/
private let _data: Array<Float64>
/** Row dimensions*/
private let _rows: Int64
/**column dimensions*/
private let _cols: Int64
/** Storage order of the matrix */
private let _storageOrder: StorageOrder
public prop rowNum: Int64 {
get() { _rows }
}
public prop colNum: Int64 {
get() { _cols }
}
public prop data: Array<Float64> {
get() { _data }
}
public prop storageOrder: StorageOrder {
get() { _storageOrder }
}
/**
* Construct an m-by-n constant matrix.
* @param m Number of rows.
* @param n Number of colums.
* @param value Fill the matrix with this scalar value.
*/
public init(m: Int64, n: Int64, value!: Float64 = 0.0) {
if (m <= 0 || n <= 0) {
throw IllegalArgumentException("Matrix dimensions must be greater than zero.")
}
_data = Array<Float64>(m * n, repeat: value)
_rows = m
_cols = n
_storageOrder = StorageOrder.RowMajor
}
/**
* Construct a matrix from a 2-D array of doubles.
* @param data Two-dimensional array of doubles.
* @exception IllegalArgumentException All rows must have the same length
*/
public init(data: Array<Array<Float64>>, order!: StorageOrder = StorageOrder.RowMajor) {
_storageOrder = order
if (data.isEmpty() || data[0].isEmpty()) {
throw IllegalArgumentException("cannot create Matrix from empty array")
}
match (order) {
case StorageOrder.RowMajor =>
_rows = data.size
_cols = data[0].size
case StorageOrder.ColumnMajor =>
_cols = data.size
_rows = data[0].size
}
this._data = data.flatten()
// Check if all rows have the same length
if ( this._data.size != _rows * _cols) {
throw IllegalArgumentException("Expected a consistent row lengths 2D array with ${_rows} rows and ${_cols} columns, but got a varying lengths in some rows.")
}
}
/** Construct a matrix from a one-dimensional array
* @param vals One-dimensional array of floats
* @param rowNum Number of rows.
* @param colNum Number of columns.
* @param order Storage order of the matrix (default is RowMajor).
* @exception IllegalArgumentException Array length must be a multiple of m.
*/
public init(data: Array<Float64>, rowNum!: Int64, colNum!: Int64 = -1, order!: StorageOrder = StorageOrder.RowMajor) {
if (data.isEmpty()) {
throw IllegalArgumentException("cannot create Matrix from empty array")
}
this._rows = rowNum
this._cols = if (colNum == -1 ) { data.size / rowNum } else { colNum }
if (this._rows <= 0 || this._cols <= 0) {
throw IllegalArgumentException("Matrix dimensions must be greater than zero.")
}
if (data.size != this._rows * this._cols) {
throw IllegalArgumentException("Expected a 1D array with size=${this._rows * this._cols}, but got size=${data.size}.")
}
this._data = data.clone()
_storageOrder = order
}
/** Construct a matrix with specified dimensions and storage order.
* <p> Note: this func will not clone the data array, it will use the reference directly. </p>
* @param data One-dimensional array of doubles.
* @param rowNum Number of rows.
* @param colNum Number of columns.
* @param order Storage order of the matrix (default is RowMajor).
* @exception IllegalArgumentException Array length must be a multiple of m.
*/
internal init(data: Array<Float64>, rowNum: Int64, colNum: Int64, order: StorageOrder, withoutClone!: Bool) {
if (data.isEmpty()) {
throw IllegalArgumentException("cannot create Matrix from empty array")
}
this._data = data
this._rows = rowNum
this._cols = colNum
this._storageOrder = order
if (this._rows <= 0 || this._cols <= 0) {
throw IllegalArgumentException("Matrix dimensions must be greater than zero.")
}
if (data.size % _cols != 0) {
throw IllegalArgumentException("Array length must be a multiple of m.")
}
}
/** Construct a matrix from another matrix.
* @param other Another matrix to copy from.
* @exception IllegalArgumentException Array length must be a multiple of m.
*/
public init(other: Matrix) {
this._data = other._data.clone()
this._rows = other._rows
this._cols = other._cols
this._storageOrder = other._storageOrder
}
/**
* Clone the Matrix object.
*/
public func clone(): Matrix {
return Matrix(this)
}
/**
* Get a one-dimensional array representation of the matrix.
* @param order Storage order of the matrix (default is RowMajor).
* @return A one-dimensional array containing the elements of the matrix.
*/
public func getOneDimensionalArray(order!: StorageOrder = StorageOrder.RowMajor): Array<Float64> {
if (this._storageOrder == order) {
return this._data.clone()
}
return this.transpose(fast: false).data
}
public operator func == (other: Matrix): Bool {
if (this._rows != other._rows || this._cols != other._cols || this._storageOrder != other._storageOrder) {
return false
}
if (this._data.size != other._data.size) {
return false
}
for (i in 0..this._rows) {
for (j in 0..this._cols) {
if (this[i, j] != other[i, j]) {
return false
}
}
}
return true
}
public operator func != (other: Matrix): Bool {
return !(this == other)
}
/** Get a single element.
* @param row Row index of the element.
* @param col Column index of the element.
* @return A(i,j)
* @exception IndexOutOfBoundsException
*/
public func get(row: Int64, col: Int64): Float64 {
this[row, col]
}
/** Set a single element.
* @param row Row index of the element.
* @param col Column index of the element.
* @param value Value to set at the specified position.
* @exception IndexOutOfBoundsException
*/
public func set(row: Int64, col: Int64, value: Float64): Unit {
this[row, col] = value
}
public operator func [] (row: Int64, col: Int64): Float64 {
if (row < 0 || row >= this._rows || col < 0 || col >= this._cols) {
throw IndexOutOfBoundsException("Index out of bounds: (${row}, ${col})")
}
match (this._storageOrder) {
case StorageOrder.RowMajor => return this._data[row * _cols + col]
case StorageOrder.ColumnMajor => return this._data[col * _rows + row]
}
}
public operator func [] (row: Int64, col: Int64, value!: Float64): Unit {
if (row < 0 || row >= this._rows || col < 0 || col >= this._cols) {
throw IndexOutOfBoundsException("Index out of bounds: (${row}, ${col})")
}
match (this._storageOrder) {
case StorageOrder.RowMajor => this._data[row * _cols + col] = value
case StorageOrder.ColumnMajor => this._data[col * _rows + row] = value
}
}
public operator func [] (row: Int64, colRange: Range<Int64>): Array<Float64> {
if (colRange.step != 1) {
throw IllegalArgumentException("Column range must have a step of 1.")
}
let cols = colRange.iterator().count()
if (row < 0 || row >= this._rows) {
throw IndexOutOfBoundsException("Index out of bounds with row: ${row}.")
}
if (colRange.start < 0 || colRange.start + cols > this._cols) {
throw IndexOutOfBoundsException("Column range out of bounds: ${this._cols} to ${colRange.start + cols}.")
}
if (cols <= 0) {
throw IllegalArgumentException("Invalid column range for extraction.")
}
match (this._storageOrder) {
case StorageOrder.RowMajor => return this.data.slice(row * this._cols + colRange.start, cols).clone()
case StorageOrder.ColumnMajor => return Array<Float64>(cols) { i =>
this.data[(colRange.start + i) * this._rows + row]
}
}
}
public operator func [] (row: Int64, colRange: Range<Int64>, value!: Array<Float64>): Unit {
if (colRange.step != 1) {
throw IllegalArgumentException("Column range must have a step of 1.")
}
let cols = colRange.iterator().count()
if (cols <= 0) {
throw IllegalArgumentException("Invalid column range for assignment.")
}
if (cols != value.size) {
throw IllegalArgumentException("Value array size must match the column range size.")
}
if (row < 0 || row >= this._rows) {
throw IndexOutOfBoundsException("Index out of bounds: (${row}, ${cols})")
}
match (this._storageOrder) {
case StorageOrder.RowMajor => value.copyTo(this.data, 0, row * this._cols + colRange.start, cols)
case StorageOrder.ColumnMajor =>
for (i in 0..cols) {
this.data[(colRange.start + i) * this._rows + row] = value[i]
}
}
}
public operator func [] (rowRange: Range<Int64>, col: Int64): Array<Float64> {
if (rowRange.step != 1) {
throw IllegalArgumentException("Row range must have a step of 1.")
}
let rows = rowRange.iterator().count()
if (col < 0 || col >= this._cols) {
throw IndexOutOfBoundsException("Index out of bounds: (${rows}, ${col})")
}
match (this._storageOrder) {
case StorageOrder.RowMajor => return Array<Float64>(rows) { i =>
this.data[(rowRange.start + i) * this._cols + col]
}
case StorageOrder.ColumnMajor => return this.data.slice(col * this._rows + rowRange.start, rows).clone()
}
}
public operator func [] (rowRange: Range<Int64>, col: Int64, value!: Array<Float64>): Unit {
if (rowRange.step != 1) {
throw IllegalArgumentException("Row range must have a step of 1.")
}
let rows = rowRange.iterator().count()
if (rows <= 0) {
throw IllegalArgumentException("Invalid row range for assignment.")
}
if (rows != value.size) {
throw IllegalArgumentException("Value array size must match the row range size.")
}
if (col < 0 || col >= this._cols) {
throw IndexOutOfBoundsException("Index out of bounds: (${rows}, ${col})")
}
match (this._storageOrder) {
case StorageOrder.RowMajor =>
for (i in 0..rows) {
this.data[(rowRange.start + i) * this._cols + col] = value[i]
}
case StorageOrder.ColumnMajor => value.copyTo(this.data, 0, col * this._rows + rowRange.start, rows)
}
}
public operator func [] (rowRange: Range<Int64>, colRange: Range<Int64>): Matrix {
if (colRange.step != 1 || rowRange.step != 1) {
throw IllegalArgumentException("Range must have a step of 1.")
}
let rows = rowRange.iterator().count()
let cols = colRange.iterator().count()
if (rows <= 0 || cols <= 0) {
throw IllegalArgumentException("Invalid range for block extraction.")
}
let subBlock = if (this._storageOrder == StorageOrder.RowMajor) {
Array<Array<Float64>>(rows) { i =>
this.data.slice((rowRange.start + i) * this._cols + colRange.start, cols).clone()
}
} else {
Array<Array<Float64>>(cols) { j =>
Array<Float64>(rows) { i =>
this.data[(colRange.start + j) * this._rows + (rowRange.start + i)]
}
}
}
return Matrix(subBlock, order: this._storageOrder)
}
public operator func [] (rowRange: Range<Int64>, colRange: Range<Int64>, value!: Matrix): Unit {
if (colRange.step != 1 || rowRange.step != 1) {
throw IllegalArgumentException("Range must have a step of 1.")
}
let rows = rowRange.iterator().count()
let cols = colRange.iterator().count()
if (rows <= 0 || cols <= 0) {
throw IllegalArgumentException("Invalid range for block extraction.")
}
if (value.rowNum != rows || value.colNum != cols) {
throw IllegalArgumentException("Value matrix dimensions must match the specified row and column ranges.")
}
if (this._storageOrder == value.storageOrder) {
if (this._storageOrder == StorageOrder.RowMajor) {
for (i in rowRange) {
value.data.copyTo(this.data, (i - rowRange.start) * cols, i * this._cols + colRange.start, cols)
}
} else {
for (i in colRange) {
value.data.copyTo(this.data, (i - colRange.start) * rows, i * this._rows + rowRange.start, rows)
}
}
} else {
for (i in 0..rows) {
for (j in 0..cols) {
this[rowRange.start + i, colRange.start + j] = value[i, j]
}
}
}
}
public operator func [] (rowIndices: Array<Int64>, colRange: Range<Int64>): Matrix {
if (colRange.step != 1) {
throw IllegalArgumentException("Column range must have a step of 1.")
}
let rowArray = Array<Array<Float64>>(rowIndices.size) { i =>
this[rowIndices[i], colRange]
}
return Matrix(rowArray)
}
public operator func [] (rowIndices: Array<Int64>, colRange: Range<Int64>, value!: Matrix): Unit {
if (colRange.step != 1) {
throw IllegalArgumentException("Column range must have a step of 1.")
}
if (rowIndices.size != value.rowNum || (colRange.iterator().count()) != value.colNum) {
throw IllegalArgumentException("Value matrix dimensions must match the specified row indices and column range.")
}
for (i in 0..rowIndices.size) {
this[rowIndices[i], colRange] = value[i, 0..value.colNum]
}
}
public operator func [] (rowRange: Range<Int64>, colIndices: Array<Int64>): Matrix {
let rows = rowRange.iterator().count()
if (rowRange.step != 1) {
throw IllegalArgumentException("Row range must have a step of 1.")
}
let colArray = Array<Array<Float64>>(rows) { i =>
Array<Float64>(colIndices.size) { j =>
this[rowRange.start + i, colIndices[j]]
}
}
return Matrix(colArray)
}
public operator func [] (rowRange: Range<Int64>, colIndices: Array<Int64>, value!: Matrix): Unit {
if (rowRange.step != 1) {
throw IllegalArgumentException("Row range must have a step of 1.")
}
if (colIndices.size != value.colNum || (rowRange.iterator().count()) != value.rowNum) {
throw IllegalArgumentException("Value matrix dimensions must match the specified row range and column indices.")
}
for (j in 0..colIndices.size) {
this[rowRange, colIndices[j]] = value[0..value.rowNum, j]
}
}
public operator func [] (rowIndices: Array<Int64>, colIndices: Array<Int64>): Matrix {
let newData = if (this._storageOrder == StorageOrder.RowMajor) {
Array<Array<Float64>>(rowIndices.size) { i =>
Array<Float64>(colIndices.size) { j =>
this[rowIndices[i], colIndices[j]]
}
}
} else {
Array<Array<Float64>>(colIndices.size) { j =>
Array<Float64>(rowIndices.size) { i =>
this[rowIndices[i], colIndices[j]]
}
}
}
return Matrix(newData, order: this._storageOrder)
}
public operator func [] (rowIndices: Array<Int64>, colIndices: Array<Int64>, value!: Matrix): Unit {
if (rowIndices.size != value.rowNum || colIndices.size != value.colNum) {
throw IllegalArgumentException("Value matrix dimensions must match the specified row and column indices.")
}
if (this._storageOrder == StorageOrder.RowMajor) {
for (i in 0..rowIndices.size) {
for (j in 0..colIndices.size) {
this[rowIndices[i], colIndices[j]] = value[i, j]
}
}
} else {
for (j in 0..colIndices.size) {
for (i in 0..rowIndices.size) {
this[rowIndices[i], colIndices[j]] = value[i, j]
}
}
}
}
/**
* Assert that the other matrix has the same size as this matrix.
* @param other Another matrix to compare with.
* @param axis Axis along which to check the size (0 for rows, 1 for columns).
* @exception IllegalArgumentException Matrices must have the same dimensions.
*/
public func concat(other: Matrix, axis!: Axis): Matrix {
let newData = if (axis == Axis.Row) { // Concatenate along rows
if (this._cols != other.colNum) {
throw IllegalArgumentException("Row sizes must match for column concatenation.")
}
match ((this._storageOrder, other.storageOrder)) {
case (StorageOrder.RowMajor, StorageOrder.RowMajor) =>
this.data.concat(other.data)
case (StorageOrder.RowMajor, StorageOrder.ColumnMajor) =>
Array<Array<Float64>>(this._rows + other.rowNum) { i =>
let rowData = Array<Float64>(this._cols, repeat: 0.0)
if (i < this._rows) {
this.data.slice(i * this._cols, this._cols).copyTo(rowData)
} else {
for (j in 0..other.colNum) {
rowData[j] = other[i - this._rows, j]
}
}
return rowData
}
case (StorageOrder.ColumnMajor, StorageOrder.ColumnMajor) =>
Array<Array<Float64>>(this._cols) { i =>
let rowData = Array<Float64>(this._rows + other.rowNum, repeat: 0.0)
this._data.slice(i * this._rows, this._rows).copyTo(rowData)
other._data.slice(i * other._rows, other._rows).copyTo(rowData, 0, this._rows, other.rowNum)
return rowData
}
case (StorageOrder.ColumnMajor, StorageOrder.RowMajor) =>
Array<Array<Float64>>(this.colNum) { i =>
let rowData = Array<Float64>(this._rows + other.rowNum, repeat: 0.0)
this._data.slice(i * this._rows, this._rows).copyTo(rowData)
for (j in 0..other.rowNum) {
rowData[this._rows + j] = other[j, i]
}
return rowData
}
}
} else if (axis == Axis.Column) { // Concatenate along columns
if (this._rows != other.rowNum) {
throw IllegalArgumentException("Column sizes must match for row concatenation.")
}
match ((this._storageOrder, other.storageOrder)) {
case (StorageOrder.RowMajor, StorageOrder.RowMajor) =>
Array<Array<Float64>>(this.rowNum) { i =>
let colData = Array<Float64>(this._cols + other.colNum, repeat: 0.0)
this._data.slice(i * this._cols, this._cols).copyTo(colData)
other.data.slice(i * other._cols, other._cols).copyTo(colData, 0, this._cols, other.colNum)
return colData
}.flatten()
case (StorageOrder.RowMajor, StorageOrder.ColumnMajor) =>
Array<Array<Float64>>(this.rowNum) { i =>
let colData = Array<Float64>(this._cols + other.colNum, repeat: 0.0)
this._data.slice(i * this._cols, this._cols).copyTo(colData)
for (j in 0..other.colNum) {
colData[this._cols + j] = other[i, j]
}
return colData
}.flatten()
case (StorageOrder.ColumnMajor, StorageOrder.ColumnMajor) =>
this._data.concat(other._data)
case (StorageOrder.ColumnMajor, StorageOrder.RowMajor) =>
Array<Array<Float64>>(this.colNum + other.colNum) { i =>
let colData = Array<Float64>(this._rows, repeat: 0.0)
if (i < this._cols) {
this._data.slice(i * this._rows, this._rows).copyTo(colData)
} else {
for (j in 0..other.rowNum) {
colData[j] = other[j, i - this._cols]
}
}
return colData
}.flatten()
}
} else {
throw IllegalArgumentException("Axis must be 0 or 1.")
}
if (newData is Array<Array<Float64>>) {
return Matrix((newData as Array<Array<Float64>>).getOrThrow(), order: this._storageOrder)
} else {
return Matrix((newData as Array<Float64>).getOrThrow(), this.rowNum, this.colNum + other.colNum, this._storageOrder, withoutClone: true)
}
}
/**
* Matrix transpose.
* @param fast If true, will only change the storage order without transposing the data. (default is false)
* @return Transposed matrix.
* @exception IllegalArgumentException If the matrix is empty.
*/
public func transpose(fast!: Bool = false): Matrix {
if (fast) {
let newStorageOrder = match (this._storageOrder) {
case StorageOrder.RowMajor => StorageOrder.ColumnMajor
case StorageOrder.ColumnMajor => StorageOrder.RowMajor
}
return Matrix(this._data.clone(), this._cols, this._rows, newStorageOrder, withoutClone: true)
}
if (this._rows == 1 || this._cols == 1) {
return Matrix(this._data.clone(), this._cols, this._rows, this._storageOrder, withoutClone: true)
}
let transposedData = Array<Float64>(_rows * _cols, repeat: 0.0)
match (this._storageOrder) {
case StorageOrder.RowMajor =>
// Convert to column-major order
for (i in 0..this._rows) {
for (j in 0..this._cols) {
transposedData[j * this._rows + i] = this._data[i * this._cols + j]
}
}
case StorageOrder.ColumnMajor =>
// Convert to row-major order
for (i in 0..this._rows) {
for (j in 0..this._cols) {
transposedData[i * this._cols + j] = this._data[j * this._rows + i]
}
}
}
return Matrix(transposedData, this._cols, this._rows, this._storageOrder, withoutClone: true)
}
public prop T: Matrix {
get() { this.transpose(fast: true) }
}
public prop I: Matrix {
get() { inverse() }
}
/**
* Unary minus
* @return -A
*/
public operator func - (): Matrix {
return this * -1.0
}
/**
* C = A + B
* @param B another matrix
* @return A + B
*/
public operator func + (other: Matrix): Matrix {
assertSameSize(other)
let resData = Array<Array<Float64>>(_rows) {
currentRow => Array<Float64>(_cols) {
currentCol =>
this[currentRow, currentCol] + other[currentRow, currentCol]
}
}
return Matrix(resData.flatten(), _rows, _cols, this._storageOrder, withoutClone: true)
}
/**
* C = A + constant
* @param value constant to add
* @return A + value
*/
public operator func + (value: Float64): Matrix {
return Matrix(this._data.map({item: Float64 => item + value}), _rows, _cols, this._storageOrder, withoutClone: true)
}
/**
* C = A - B
* @param B another matrix
* @return A - B
*/
public operator func - (other: Matrix): Matrix {
assertSameSize(other)
let resData = Array<Array<Float64>>(_rows) {
currentRow => Array<Float64>(_cols) {
currentCol =>
this[currentRow, currentCol] - other[currentRow, currentCol]
}
}
return Matrix(resData.flatten(), _rows, _cols, this._storageOrder, withoutClone: true)
}
/**
* C = A - constant
* @param value constant to subtract
* @return A - value
*/
public operator func - (value: Float64): Matrix {
return Matrix(this._data.map({item: Float64 => item - value}), _rows, _cols, this._storageOrder, withoutClone: true)
}
/**
* Multiply a matrix by a scalar, C = s*A
* @param s scalar
* @return s*A
*/
public operator func * (s: Float64): Matrix {
return Matrix(this._data.map({item: Float64 => item * s}), _rows, _cols, this._storageOrder, withoutClone: true)
}
/**
* Linear algebraic matrix multiplication, A * B
* @param B another matrix
* @return Matrix product, A * B
* @exception IllegalArgumentException Matrix inner dimensions must agree.
*/
public operator func * (other: Matrix): Matrix {
return this.matMul(other)
}
/**
* Linear algebraic matrix multiplication, A * B
* @param B another matrix
* @return Matrix product, A * B
* @exception IllegalArgumentException Matrix inner dimensions must agree.
*/
public func matMul(other: Matrix): Matrix {
let M = this._rows
let K = this._cols
let N = other._cols
if (K != other.rowNum) {
throw IllegalArgumentException("Matrix inner dimensions must agree.")
}
let result = Matrix(_rows, other.colNum)
match ((this._storageOrder, other._storageOrder)) {
case (StorageOrder.RowMajor, StorageOrder.RowMajor) =>
// Case 1: A行主序, B行主序
let A = this
let BT = other.transpose(fast: false)
for (i in 0..M) {
for (j in 0..N) {
var sum: Float64 = 0.0
for (k in 0..K) {
sum += A[i, k] * BT[j, k]
}
result[i, j] = sum
}
}
case (StorageOrder.ColumnMajor, StorageOrder.ColumnMajor) =>
// Case 2: A列主序, B列主序
let AT = this.transpose(fast: false)
let B = other
for (i in 0..M) {
for (j in 0..N) {
var sum: Float64 = 0.0
for (k in 0..K) {
sum += AT[k, i] * B[k, j]
}
result[i, j] = sum
}
}
case (StorageOrder.ColumnMajor, StorageOrder.RowMajor) =>
// Case 3: A列主序, B行主序
// 使用 ikj 顺序优化缓存访问
let A = this
let B = other
for (i in 0..M) {
for (k in 0..K) {
let a_ik = A[i, k]
for (j in 0..N) {
result[i, j] += a_ik * B[k, j]
}
}
}
case (StorageOrder.RowMajor, StorageOrder.ColumnMajor) =>
// Case 4: A行主序, B列主序
// 最优情况：A按行访问，B按列访问，都是连续的
let A = this
let B = other
for (i in 0..M) {
for (j in 0..N) {
var sum: Float64 = 0.0
for (k in 0..K) {
sum += A[i, k] * B[k, j]
}
result[i, j] = sum
}
}
}
return result
}
/**
* Multiply an Identity matrix
* @param rowNum Number of rows.
* @param colNum Number of columns.
* @return Identity matrix of size rowNum x colNum.
* @exception IllegalArgumentException If rowNum or colNum is less than or equal
* to zero.
*/
public static func Identity(rowNum: Int64, colNum: Int64): Matrix {
if (rowNum <= 0 || colNum <= 0) {
throw IllegalArgumentException("Matrix dimensions must be greater than zero.")
}
let data = Array<Float64>(rowNum * colNum, repeat: 0.0)
for (i in 0..min(rowNum, colNum)) {
data[i * colNum + i] = 1.0
}
return Matrix(data, rowNum, colNum, StorageOrder.RowMajor, withoutClone: true)
}
public func identity(): Matrix {
return Matrix.Identity(this._rows, this._cols)
}
/**
* Create a random matrix with specified dimensions.
* @param rowNum Number of rows.
* @param colNum Number of columns.
* @return A matrix filled with random values in the range [-1.0, 1.0].
* @exception IllegalArgumentException If rowNum or colNum is less than or equal to zero.
*/
public static func Random(rowNum: Int64, colNum: Int64): Matrix {
if (rowNum <= 0 || colNum <= 0) {
throw IllegalArgumentException("Matrix dimensions must be greater than zero.")
}
let data = Array<Float64>(rowNum * colNum) {
_ => RandomGenerator().nextFloat64()
}
return Matrix(data, rowNum, colNum, StorageOrder.RowMajor, withoutClone: true)
}
public func random(): Matrix {
return Matrix.Random(rowNum, colNum)
}
/**
* L1 norm
* @return maximum column sum.
*/
public func normL1(): Float64 {
var f: Float64 = 0.0
for (j in 0..this._cols) {
var s: Float64 = 0.0
for (i in 0..this._rows) {
s += abs(this[i, j])
}
f = max(f, s)
}
return f
}
/**
* L2 norm
* @return maximum singular value.
*/
public func normL2(): Float64 {
return (SingularValueDecomposition(this).norm2())
}
/**
* Infinity norm
* @return maximum row sum.
*/
public func normInf(): Float64 {
var f = 0.0
for (i in 0..this._rows) {
var s = 0.0
for (j in 0..this._cols) {
s += abs(this[i, j])
}
f = max(f, s)
}
return f
}
/**
* Frobenius norm
* @return sqrt of sum of squares of all elements.
*/
public func normFrobenius(): Float64 {
var f = 0.0
for (i in 0..this._rows) {
for (j in 0..this._cols) {
f = Maths.hypot(f, this[i, j])
}
}
return f
}
/**
* LU Decomposition
* @return LUDecomposition
* @see LUDecomposition
*/
public func lu(): LUDecomposition {
return LUDecomposition(this)
}
/**
* QR Decomposition
* @return QRDecomposition
* @see QRDecomposition
*/
public func qr(): QRDecomposition {
return QRDecomposition(this)
}
/**
* Cholesky Decomposition
* @return CholeskyDecomposition
* @see CholeskyDecomposition
*/
public func chol(): CholeskyDecomposition {
return CholeskyDecomposition(this)
}
/**
* Singular Value Decomposition
* @return SingularValueDecomposition
* @see SingularValueDecomposition
*/
public func svd(): SingularValueDecomposition {
return SingularValueDecomposition(this)
}
/**
* Eigenvalue Decomposition
* @return EigenvalueDecomposition
* @see EigenvalueDecomposition
*/
public func eig(): EigenvalueDecomposition {
return EigenvalueDecomposition(this)
}
/**
* Solve A*X = B
* @param B right hand side
* @return solution if A is square, least squares solution otherwise
*/
public func solve(B: Matrix): Matrix {
return (if (_rows == _cols) {
(LUDecomposition(this)).solve(B)
} else {
(QRDecomposition(this)).solve(B)
})
}
/**
* Solve X*A = B, which is also A'*X' = B'
* @param B right hand side
* @return solution if A is square, least squares solution otherwise.
*/
public func solveTranspose(B: Matrix): Matrix {
return transpose().solve(B.transpose())
}
/**
* Matrix inverse or pseudoinverse
* @return inverse(A) if A is square, pseudoinverse otherwise.
*/
public func inverse(): Matrix {
return solve(Matrix.Identity(_rows, _rows))
}
/**
* Matrix determinant
* @return determinant
*/
public func det(): Float64 {
return LUDecomposition(this).det()
}
/**
* Matrix rank
* @return effective numerical rank, obtained from SVD.
*/
public func rank(): Int64 {
return SingularValueDecomposition(this).rank()
}
/**
* Matrix condition (2 norm)
* @return ratio of largest to smallest singular value.
*/
public func cond(): Float64 {
return SingularValueDecomposition(this).cond()
}
/**
* Matrix trace.
* @return sum of the diagonal elements.
*/
public func trace(): Float64 {
var t: Float64 = 0.0
for (i in 0..min(_rows, _cols)) {
t += this[i, i]
}
return t
}
/**
* Element-wise operations wrapper.
* @return ElementWiseWapper for this matrix.
*/
public func elementWise(): ElementWiseWapper<Matrix> {
return ElementWiseWapper(this)
}
public func elementWiseMul(other: Matrix): Matrix {
assertSameSize(other)
let resData = Array<Array<Float64>>(_rows) {
currentRow => Array<Float64>(_cols) {
currentCol =>
this[currentRow, currentCol] * other[currentRow, currentCol]
}
}
return Matrix(resData.flatten(), _rows, _cols, this._storageOrder, withoutClone: true)
}
public func elementWiseDiv(other: Matrix): Matrix {
assertSameSize(other)
let resData = Array<Array<Float64>>(_rows) {
currentRow => Array<Float64>(_cols) {
currentCol =>
this[currentRow, currentCol] / other[currentRow, currentCol]
}
}
return Matrix(resData.flatten(), _rows, _cols, this._storageOrder, withoutClone: true)
}
public operator func *(other: ElementWiseWapper<Matrix>): Matrix {
return this.elementWiseMul(other.ref)
}
public operator func /(other: ElementWiseWapper<Matrix>): Matrix {
return this.elementWiseDiv(other.ref)
}
private func assertSameSize(other: Matrix): Unit {
if (this.rowNum != other.rowNum || this.colNum != other.colNum) {
throw IllegalArgumentException("Matrix dimensions must agree.")
}
}
/** Construct a matrix quickly without checking arguments.
* @param A Two-dimensional array of doubles.
* @param m Number of rows.
* @param n Number of colums.
*/
@Deprecated["Use Matrix(data: Array<Array<Float64>>) instead.", since: "2.0.0", strict: true]
public init(A: Array<Array<Float64>>, m: Int64, n: Int64) {
this(A)
}
/** Construct a matrix from a one-dimensional packed array
* @param vals One-dimensional array of doubles, packed by columns (ala Fortran).
* @param m Number of rows.
* @exception IllegalArgumentException Array length must be a multiple of m.
*/
@Deprecated["Use Matrix(data: Array<Float64>, rowNum!: ?Int64, colNum!: ?Int64, order!: StorageOrder ) instead.", since: "2.0.0", strict: true]
public init(vals: Array<Float64>, m: Int64) {
this(vals, rowNum: m)
}
/** Construct a matrix from a copy of a 2-D array.
* @param A Two-dimensional array of doubles.
* @exception IllegalArgumentException All rows must have the same length
*/
@Deprecated["All Constructors will copy the data", since: "2.0.0", strict: true]
public static func constructWithCopy(A: Array<Array<Float64>>): Matrix {
return Matrix(A)
}
@Deprecated["Use Matrix::clone() instead.", since: "2.0.0", strict: true]
public func copy(): Matrix {
return this.clone()
}
/** Access the internal two-dimensional array.
* @return Pointer to the two-dimensional array of matrix elements.
*/
@Deprecated["Get the internal array are not allowed.", since: "2.0.0", strict: true]
public func getArray(): Array<Array<Float64>> {
return Array<Array<Float64>>()
}
/** Copy the internal two-dimensional array.
* @return Two-dimensional array copy of matrix elements.
*/
@Deprecated["Get the internal array are not allowed.", since: "2.0.0", strict: true]
public func getArrayCopy(): Array<Array<Float64>> {
return Array<Array<Float64>>()
}
/** Make a one-dimensional column packed copy of the internal array.
* @return Matrix elements packed in a one-dimensional array by columns.
*/
@Deprecated["Use Matrix::getOneDimensionalArray(order: StorageOrder) instead.", since: "2.0.0", strict: true]
public func getColumnPackedCopy(): Array<Float64> {
return getOneDimensionalArray(order: StorageOrder.ColumnMajor)
}
/** Make a one-dimensional row packed copy of the internal array.
* @return Matrix elements packed in a one-dimensional array by rows.
*/
@Deprecated["Use Matrix::getOneDimensionalArray(order: StorageOrder) instead.", since: "2.0.0", strict: true]
public func getRowPackedCopy(): Array<Float64> {
return getOneDimensionalArray(order: StorageOrder.RowMajor)
}
/** Get row dimension.
* @return m, the number of rows.
*/
@Deprecated["Use Matrix::rowNum instead.", since: "2.0.0", strict: true]
public func getRowDimension(): Int64 {
return rowNum
}
/** Get column dimension.
* @return n, the number of columns.
*/
@Deprecated["Use Matrix::colNum instead.", since: "2.0.0", strict: true]
public func getColumnDimension(): Int64 {
return colNum
}
/** Get a submatrix.
* @param i0 Initial row index
* @param i1 Final row index
* @param j0 Initial column index
* @param j1 Final column index
* @return A(i0:i1,j0:j1)
* @exception IndexOutOfBoundsException Submatrix indices
*/
@Deprecated["Use operator [rowRange: Range<Int64>, colRange: Range<Int64>] instead.", since: "2.0.0", strict: true]
public func getMatrix(i0: Int64, i1: Int64, j0: Int64, j1: Int64): Matrix {
return Matrix(this[Range<Int64>(i0, i1, 1, true, true, true), Range<Int64>(j0, j1, 1, true, true, true)])
}
/** Get a submatrix.
* @param r Array of row indices.
* @param c Array of column indices.
* @return A(r(:),c(:))
* @exception IndexOutOfBoundsException Submatrix indices
*/
@Deprecated["Use operator [rowRange: Range<Int64>, colRange: Range<Int64>] instead.", since: "2.0.0", strict: true]
public func getMatrix(r: Array<Int64>, c: Array<Int64>): Matrix {
throw UnsupportedException("Use operator [rowRange: Range<Int64>, colRange: Range<Int64>] instead.")
return Matrix(1,1)
}
/** Get a submatrix.
* @param i0 Initial row index
* @param i1 Final row index
* @param c Array of column indices.
* @return A(i0:i1,c(:))
* @exception IndexOutOfBoundsException Submatrix indices
*/
@Deprecated["Use operator [rowRange: Range<Int64>, colRange: Range<Int64>] instead.", since: "2.0.0", strict: true]
public func getMatrix(i0: Int64, i1: Int64, c: Array<Int64>): Matrix {
throw UnsupportedException("Use operator [rowRange: Range<Int64>, colRange: Range<Int64>] instead.")
return Matrix(1,1)
}
/** Get a submatrix.
* @param r Array of row indices.
* @param j0 Initial column index
* @param j1 Final column index
* @return A(r(:),j0:j1)
* @exception IndexOutOfBoundsException Submatrix indices
*/
@Deprecated["Use operator [rowRange: Range<Int64>, colRange: Range<Int64>] instead.", since: "2.0.0", strict: true]
public func getMatrix(r: Array<Int64>, j0: Int64, j1: Int64): Matrix {
throw UnsupportedException("Use operator [rowRange: Range<Int64>, colRange: Range<Int64>] instead.")
return Matrix(1,1)
}
/** Set a submatrix.
* @param i0 Initial row index
* @param i1 Final row index
* @param j0 Initial column index
* @param j1 Final column index
* @param X A(i0:i1,j0:j1)
* @exception IndexOutOfBoundsException Submatrix indices
*/
@Deprecated["Use operator [rowRange: Range<Int64>, colRange: Range<Int64>, value!: Matrix] instead.", since: "2.0.0", strict: true]
public func setMatrix(i0: Int64, i1: Int64, j0: Int64, j1: Int64, X: Matrix): Unit {
throw UnsupportedException("Use operator [rowRange: Range<Int64>, colRange: Range<Int64>, value!: Matrix] instead.")
}
/** Set a submatrix.
* @param r Array of row indices.
* @param c Array of column indices.
* @param X A(r(:),c(:))
* @exception ArrayIndexOutOfBoundsException Submatrix indices
*/
@Deprecated["Use operator [rowRange: Range<Int64>, colRange: Range<Int64>, value!: Matrix] instead.", since: "2.0.0", strict: true]
public func setMatrix(r: Array<Int64>, c: Array<Int64>, X: Matrix): Unit {
throw UnsupportedException("Use operator [rowRange: Range<Int64>, colRange: Range<Int64>, value!: Matrix] instead.")
}
/** Set a submatrix.
* @param r Array of row indices.
* @param j0 Initial column index
* @param j1 Final column index
* @param X A(r(:),j0:j1)
* @exception IndexOutOfBoundsException Submatrix indices
*/
@Deprecated["Use operator [rowRange: Range<Int64>, colRange: Range<Int64>, value!: Matrix] instead.", since: "2.0.0", strict: true]
public func setMatrix(r: Array<Int64>, j0: Int64, j1: Int64, X: Matrix): Unit {
throw UnsupportedException("Use operator [rowRange: Range<Int64>, colRange: Range<Int64>, value!: Matrix] instead.")
}
/** Set a submatrix.
* @param i0 Initial row index
* @param i1 Final row index
* @param c Array of column indices.
* @param X A(i0:i1,c(:))
* @exception IndexOutOfBoundsException Submatrix indices
*/
@Deprecated["Use operator [rowRange: Range<Int64>, colRange: Range<Int64>, value!: Matrix] instead.", since: "2.0.0", strict: true]
public func setMatrix(i0: Int64, i1: Int64, c: Array<Int64>, X: Matrix): Unit {
throw UnsupportedException("Use operator [rowRange: Range<Int64>, colRange: Range<Int64>, value!: Matrix] instead.")
}
/**
* Unary minus
* @return -A
*/
@Deprecated["Use operator -() instead.", since: "2.0.0", strict: true]
public func uminus(): Matrix {
return -this
}
/**
* C = A + B
* @param B another matrix
* @return A + B
*/
@Deprecated["Use operator +() instead.", since: "2.0.0", strict: true]
public func plus(B: Matrix): Matrix {
return this + B
}
/**
* A = A - B
* @param B another matrix
* @return A - B
*/
@Deprecated["Use A -= B instead.", since: "2.0.0", strict: true]
public func minusEquals(B: Matrix): Matrix {
throw UnsupportedException("Use A -= B instead.")
return this
}
/**
* Element-by-element multiplication, C = A.*B
* @param B another matrix
* @return A.*B
*/
@Deprecated["Use A.elementWise() * B instead.", since: "2.0.0", strict: true]
public func arrayTimes(B: Matrix): Matrix {
throw UnsupportedException("Use Matrix::elementWise() * B instead.")
return this
}
/**
* Element-by-element multiplication in place, A = A.*B
* @param B another matrix
* @return A.*B
*/
@Deprecated["Use A *= B.elementWise() instead.", since: "2.0.0", strict: true]
public func arrayTimesEquals(B: Matrix): Matrix {
throw UnsupportedException("Use A *= B.elementWise() instead.")
return this
}
/**
* Element-by-element right division, C = A./B
* @param B another matrix
* @return A./B
*/
@Deprecated["Use B.elementWise() / A instead.", since: "2.0.0", strict: true]
public func arrayRightDivide(B: Matrix): Matrix {
throw UnsupportedException("Use B.elementWise() / A instead.")
return this
}
/**
* Element-by-element right division in place, A = A./B
* @param B another matrix
* @return A./B
*/
@Deprecated["Use B /= A.elementWise() instead.", since: "2.0.0", strict: true]
public func arrayRightDivideEquals(B: Matrix): Matrix {
throw UnsupportedException("Use B /= A.elementWise() instead.")
return this
}
/**
* Element-by-element left division, C = A.\B
* @param B another matrix
* @return A.\B
*/
@Deprecated["Use A.elementWise() / B instead.", since: "2.0.0", strict: true]
public func arrayLeftDivide(B: Matrix): Matrix {
throw UnsupportedException("Use A.elementWise() / B instead.")
return this
}
/**
* Element-by-element left division in place, A = A.\B
* @param B another matrix
* @return A.\B
*/
@Deprecated["Use A /= B.elementWise() instead.", since: "2.0.0", strict: true]
public func arrayLeftDivideEquals(B: Matrix): Matrix {
throw UnsupportedException("Use A /= B.elementWise() instead.")
return this
}
/**
* Multiply a matrix by a scalar, C = s*A
* @param s scalar
* @return s*A
*/
@Deprecated["Use operator *() instead.", since: "2.0.0", strict: true]
public func times(s: Float64): Matrix {
return this * s
}
/**
* Multiply a matrix by a scalar in place, A = s*A
* @param s scalar
* @return replace A by s*A
*/
@Deprecated["Use A *= s instead.", since: "2.0.0", strict: true]
public func timesEquals(s: Float64): Matrix {
throw UnsupportedException("Use A *= s instead.")
return this
}
/**
* Linear algebraic matrix multiplication, A * B
* @param B another matrix
* @return Matrix product, A * B
* @exception IllegalArgumentException Matrix inner dimensions must agree.
*/
@Deprecated["Use operator *() instead.", since: "2.0.0", strict: true]
public func times(B: Matrix): Matrix {
return this * B
}
@Deprecated["Deprecated", since: "2.0.0", strict: true]
private func checkMatrixDimensions(B: Matrix): Unit {
""
}
}
extend Matrix <: ToString {
public func toString(): String {
toString("7.2")
}
public func toString(format: String): String {
let sb = StringBuilder()
for (i in 0..rowNum) {
sb.append(this[i, 0..colNum].map({item: Float64 => item.format(format) + ", "}))
if (i < rowNum - 1) {
sb.append("\n")
}
}
sb.toString()
}
}
