use num::{One, Zero};
#[cfg(feature = "abomonation-serialize")]
use std::io::{Result as IOResult, Write};
use approx::{AbsDiffEq, RelativeEq, UlpsEq};
use std::any::TypeId;
use std::cmp::Ordering;
use std::fmt;
use std::hash::{Hash, Hasher};
use std::marker::PhantomData;
use std::mem;
#[cfg(feature = "serde-serialize")]
use serde::{Deserialize, Deserializer, Serialize, Serializer};
#[cfg(feature = "abomonation-serialize")]
use abomonation::Abomonation;
use simba::scalar::{ClosedAdd, ClosedMul, ClosedSub, Field};
use simba::simd::SimdPartialOrd;
use crate::base::allocator::{Allocator, SameShapeAllocator, SameShapeC, SameShapeR};
use crate::base::constraint::{DimEq, SameNumberOfColumns, SameNumberOfRows, ShapeConstraint};
use crate::base::dimension::{Dim, DimAdd, DimSum, IsNotStaticOne, U1, U2, U3};
use crate::base::iter::{
ColumnIter, ColumnIterMut, MatrixIter, MatrixIterMut, RowIter, RowIterMut,
};
use crate::base::storage::{
ContiguousStorage, ContiguousStorageMut, Owned, SameShapeStorage, Storage, StorageMut,
};
use crate::base::{DefaultAllocator, MatrixMN, MatrixN, Scalar, Unit, VectorN};
use crate::SimdComplexField;
pub type SquareMatrix<N, D, S> = Matrix<N, D, D, S>;
pub type Vector<N, D, S> = Matrix<N, D, U1, S>;
pub type RowVector<N, D, S> = Matrix<N, U1, D, S>;
pub type MatrixSum<N, R1, C1, R2, C2> =
Matrix<N, SameShapeR<R1, R2>, SameShapeC<C1, C2>, SameShapeStorage<N, R1, C1, R2, C2>>;
pub type VectorSum<N, R1, R2> =
Matrix<N, SameShapeR<R1, R2>, U1, SameShapeStorage<N, R1, U1, R2, U1>>;
pub type MatrixCross<N, R1, C1, R2, C2> =
Matrix<N, SameShapeR<R1, R2>, SameShapeC<C1, C2>, SameShapeStorage<N, R1, C1, R2, C2>>;
#[repr(C)]
#[derive(Clone, Copy)]
pub struct Matrix<N: Scalar, R: Dim, C: Dim, S> {
pub data: S,
_phantoms: PhantomData<(N, R, C)>,
}
impl<N: Scalar, R: Dim, C: Dim, S: fmt::Debug> fmt::Debug for Matrix<N, R, C, S> {
fn fmt(&self, formatter: &mut fmt::Formatter) -> Result<(), fmt::Error> {
formatter
.debug_struct("Matrix")
.field("data", &self.data)
.finish()
}
}
impl<N, R, C, S> Default for Matrix<N, R, C, S>
where
N: Scalar,
R: Dim,
C: Dim,
S: Default,
{
fn default() -> Self {
Matrix {
data: Default::default(),
_phantoms: PhantomData,
}
}
}
#[cfg(feature = "serde-serialize")]
impl<N, R, C, S> Serialize for Matrix<N, R, C, S>
where
N: Scalar,
R: Dim,
C: Dim,
S: Serialize,
{
fn serialize<T>(&self, serializer: T) -> Result<T::Ok, T::Error>
where
T: Serializer,
{
self.data.serialize(serializer)
}
}
#[cfg(feature = "serde-serialize")]
impl<'de, N, R, C, S> Deserialize<'de> for Matrix<N, R, C, S>
where
N: Scalar,
R: Dim,
C: Dim,
S: Deserialize<'de>,
{
fn deserialize<D>(deserializer: D) -> Result<Self, D::Error>
where
D: Deserializer<'de>,
{
S::deserialize(deserializer).map(|x| Matrix {
data: x,
_phantoms: PhantomData,
})
}
}
#[cfg(feature = "abomonation-serialize")]
impl<N: Scalar, R: Dim, C: Dim, S: Abomonation> Abomonation for Matrix<N, R, C, S> {
unsafe fn entomb<W: Write>(&self, writer: &mut W) -> IOResult<()> {
self.data.entomb(writer)
}
unsafe fn exhume<'a, 'b>(&'a mut self, bytes: &'b mut [u8]) -> Option<&'b mut [u8]> {
self.data.exhume(bytes)
}
fn extent(&self) -> usize {
self.data.extent()
}
}
#[cfg(feature = "compare")]
impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> matrixcompare_core::Matrix<N>
for Matrix<N, R, C, S>
{
fn rows(&self) -> usize {
self.nrows()
}
fn cols(&self) -> usize {
self.ncols()
}
fn access(&self) -> matrixcompare_core::Access<N> {
matrixcompare_core::Access::Dense(self)
}
}
#[cfg(feature = "compare")]
impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> matrixcompare_core::DenseAccess<N>
for Matrix<N, R, C, S>
{
fn fetch_single(&self, row: usize, col: usize) -> N {
self.index((row, col)).clone()
}
}
impl<N: Scalar, R: Dim, C: Dim, S> Matrix<N, R, C, S> {
#[inline]
pub unsafe fn from_data_statically_unchecked(data: S) -> Matrix<N, R, C, S> {
Matrix {
data,
_phantoms: PhantomData,
}
}
}
impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
#[inline]
pub fn from_data(data: S) -> Self {
unsafe { Self::from_data_statically_unchecked(data) }
}
#[inline]
pub fn shape(&self) -> (usize, usize) {
let (nrows, ncols) = self.data.shape();
(nrows.value(), ncols.value())
}
#[inline]
pub fn nrows(&self) -> usize {
self.shape().0
}
#[inline]
pub fn ncols(&self) -> usize {
self.shape().1
}
#[inline]
pub fn strides(&self) -> (usize, usize) {
let (srows, scols) = self.data.strides();
(srows.value(), scols.value())
}
#[inline]
pub fn vector_to_matrix_index(&self, i: usize) -> (usize, usize) {
let (nrows, ncols) = self.shape();
if nrows == 1 {
(0, i)
} else if ncols == 1 {
(i, 0)
} else {
(i % nrows, i / nrows)
}
}
#[inline]
pub fn as_ptr(&self) -> *const N {
self.data.ptr()
}
#[inline]
pub fn relative_eq<R2, C2, SB>(
&self,
other: &Matrix<N, R2, C2, SB>,
eps: N::Epsilon,
max_relative: N::Epsilon,
) -> bool
where
N: RelativeEq,
R2: Dim,
C2: Dim,
SB: Storage<N, R2, C2>,
N::Epsilon: Copy,
ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>,
{
assert!(self.shape() == other.shape());
self.iter()
.zip(other.iter())
.all(|(a, b)| a.relative_eq(b, eps, max_relative))
}
#[inline]
pub fn eq<R2, C2, SB>(&self, other: &Matrix<N, R2, C2, SB>) -> bool
where
N: PartialEq,
R2: Dim,
C2: Dim,
SB: Storage<N, R2, C2>,
ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>,
{
assert!(self.shape() == other.shape());
self.iter().zip(other.iter()).all(|(a, b)| *a == *b)
}
#[inline]
pub fn into_owned(self) -> MatrixMN<N, R, C>
where
DefaultAllocator: Allocator<N, R, C>,
{
Matrix::from_data(self.data.into_owned())
}
#[inline]
pub fn into_owned_sum<R2, C2>(self) -> MatrixSum<N, R, C, R2, C2>
where
R2: Dim,
C2: Dim,
DefaultAllocator: SameShapeAllocator<N, R, C, R2, C2>,
ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>,
{
if TypeId::of::<SameShapeStorage<N, R, C, R2, C2>>() == TypeId::of::<Owned<N, R, C>>() {
unsafe {
let owned = self.into_owned();
let res = mem::transmute_copy(&owned);
mem::forget(owned);
res
}
} else {
self.clone_owned_sum()
}
}
#[inline]
pub fn clone_owned(&self) -> MatrixMN<N, R, C>
where
DefaultAllocator: Allocator<N, R, C>,
{
Matrix::from_data(self.data.clone_owned())
}
#[inline]
pub fn clone_owned_sum<R2, C2>(&self) -> MatrixSum<N, R, C, R2, C2>
where
R2: Dim,
C2: Dim,
DefaultAllocator: SameShapeAllocator<N, R, C, R2, C2>,
ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>,
{
let (nrows, ncols) = self.shape();
let nrows: SameShapeR<R, R2> = Dim::from_usize(nrows);
let ncols: SameShapeC<C, C2> = Dim::from_usize(ncols);
let mut res: MatrixSum<N, R, C, R2, C2> =
unsafe { Matrix::new_uninitialized_generic(nrows, ncols) };
for j in 0..res.ncols() {
for i in 0..res.nrows() {
unsafe {
*res.get_unchecked_mut((i, j)) = self.get_unchecked((i, j)).inlined_clone();
}
}
}
res
}
#[inline]
pub fn transpose_to<R2, C2, SB>(&self, out: &mut Matrix<N, R2, C2, SB>)
where
R2: Dim,
C2: Dim,
SB: StorageMut<N, R2, C2>,
ShapeConstraint: SameNumberOfRows<R, C2> + SameNumberOfColumns<C, R2>,
{
let (nrows, ncols) = self.shape();
assert!(
(ncols, nrows) == out.shape(),
"Incompatible shape for transpose-copy."
);
for i in 0..nrows {
for j in 0..ncols {
unsafe {
*out.get_unchecked_mut((j, i)) = self.get_unchecked((i, j)).inlined_clone();
}
}
}
}
#[inline]
#[must_use = "Did you mean to use transpose_mut()?"]
pub fn transpose(&self) -> MatrixMN<N, C, R>
where
DefaultAllocator: Allocator<N, C, R>,
{
let (nrows, ncols) = self.data.shape();
unsafe {
let mut res = Matrix::new_uninitialized_generic(ncols, nrows);
self.transpose_to(&mut res);
res
}
}
}
impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
#[inline]
pub fn map<N2: Scalar, F: FnMut(N) -> N2>(&self, mut f: F) -> MatrixMN<N2, R, C>
where
DefaultAllocator: Allocator<N2, R, C>,
{
let (nrows, ncols) = self.data.shape();
let mut res = unsafe { MatrixMN::new_uninitialized_generic(nrows, ncols) };
for j in 0..ncols.value() {
for i in 0..nrows.value() {
unsafe {
let a = self.data.get_unchecked(i, j).inlined_clone();
*res.data.get_unchecked_mut(i, j) = f(a)
}
}
}
res
}
#[inline]
pub fn fold_with<N2>(
&self,
init_f: impl FnOnce(Option<&N>) -> N2,
f: impl FnMut(N2, &N) -> N2,
) -> N2 {
let mut it = self.iter();
let init = init_f(it.next());
it.fold(init, f)
}
#[inline]
pub fn map_with_location<N2: Scalar, F: FnMut(usize, usize, N) -> N2>(
&self,
mut f: F,
) -> MatrixMN<N2, R, C>
where
DefaultAllocator: Allocator<N2, R, C>,
{
let (nrows, ncols) = self.data.shape();
let mut res = unsafe { MatrixMN::new_uninitialized_generic(nrows, ncols) };
for j in 0..ncols.value() {
for i in 0..nrows.value() {
unsafe {
let a = self.data.get_unchecked(i, j).inlined_clone();
*res.data.get_unchecked_mut(i, j) = f(i, j, a)
}
}
}
res
}
#[inline]
pub fn zip_map<N2, N3, S2, F>(&self, rhs: &Matrix<N2, R, C, S2>, mut f: F) -> MatrixMN<N3, R, C>
where
N2: Scalar,
N3: Scalar,
S2: Storage<N2, R, C>,
F: FnMut(N, N2) -> N3,
DefaultAllocator: Allocator<N3, R, C>,
{
let (nrows, ncols) = self.data.shape();
let mut res = unsafe { MatrixMN::new_uninitialized_generic(nrows, ncols) };
assert_eq!(
(nrows.value(), ncols.value()),
rhs.shape(),
"Matrix simultaneous traversal error: dimension mismatch."
);
for j in 0..ncols.value() {
for i in 0..nrows.value() {
unsafe {
let a = self.data.get_unchecked(i, j).inlined_clone();
let b = rhs.data.get_unchecked(i, j).inlined_clone();
*res.data.get_unchecked_mut(i, j) = f(a, b)
}
}
}
res
}
#[inline]
pub fn zip_zip_map<N2, N3, N4, S2, S3, F>(
&self,
b: &Matrix<N2, R, C, S2>,
c: &Matrix<N3, R, C, S3>,
mut f: F,
) -> MatrixMN<N4, R, C>
where
N2: Scalar,
N3: Scalar,
N4: Scalar,
S2: Storage<N2, R, C>,
S3: Storage<N3, R, C>,
F: FnMut(N, N2, N3) -> N4,
DefaultAllocator: Allocator<N4, R, C>,
{
let (nrows, ncols) = self.data.shape();
let mut res = unsafe { MatrixMN::new_uninitialized_generic(nrows, ncols) };
assert_eq!(
(nrows.value(), ncols.value()),
b.shape(),
"Matrix simultaneous traversal error: dimension mismatch."
);
assert_eq!(
(nrows.value(), ncols.value()),
c.shape(),
"Matrix simultaneous traversal error: dimension mismatch."
);
for j in 0..ncols.value() {
for i in 0..nrows.value() {
unsafe {
let a = self.data.get_unchecked(i, j).inlined_clone();
let b = b.data.get_unchecked(i, j).inlined_clone();
let c = c.data.get_unchecked(i, j).inlined_clone();
*res.data.get_unchecked_mut(i, j) = f(a, b, c)
}
}
}
res
}
#[inline]
pub fn fold<Acc>(&self, init: Acc, mut f: impl FnMut(Acc, N) -> Acc) -> Acc {
let (nrows, ncols) = self.data.shape();
let mut res = init;
for j in 0..ncols.value() {
for i in 0..nrows.value() {
unsafe {
let a = self.data.get_unchecked(i, j).inlined_clone();
res = f(res, a)
}
}
}
res
}
#[inline]
pub fn zip_fold<N2, R2, C2, S2, Acc>(
&self,
rhs: &Matrix<N2, R2, C2, S2>,
init: Acc,
mut f: impl FnMut(Acc, N, N2) -> Acc,
) -> Acc
where
N2: Scalar,
R2: Dim,
C2: Dim,
S2: Storage<N2, R2, C2>,
ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>,
{
let (nrows, ncols) = self.data.shape();
let mut res = init;
assert_eq!(
(nrows.value(), ncols.value()),
rhs.shape(),
"Matrix simultaneous traversal error: dimension mismatch."
);
for j in 0..ncols.value() {
for i in 0..nrows.value() {
unsafe {
let a = self.data.get_unchecked(i, j).inlined_clone();
let b = rhs.data.get_unchecked(i, j).inlined_clone();
res = f(res, a, b)
}
}
}
res
}
#[inline]
pub fn apply<F: FnMut(N) -> N>(&mut self, mut f: F)
where
S: StorageMut<N, R, C>,
{
let (nrows, ncols) = self.shape();
for j in 0..ncols {
for i in 0..nrows {
unsafe {
let e = self.data.get_unchecked_mut(i, j);
*e = f(e.inlined_clone())
}
}
}
}
#[inline]
pub fn zip_apply<N2, R2, C2, S2>(
&mut self,
rhs: &Matrix<N2, R2, C2, S2>,
mut f: impl FnMut(N, N2) -> N,
) where
S: StorageMut<N, R, C>,
N2: Scalar,
R2: Dim,
C2: Dim,
S2: Storage<N2, R2, C2>,
ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>,
{
let (nrows, ncols) = self.shape();
assert_eq!(
(nrows, ncols),
rhs.shape(),
"Matrix simultaneous traversal error: dimension mismatch."
);
for j in 0..ncols {
for i in 0..nrows {
unsafe {
let e = self.data.get_unchecked_mut(i, j);
let rhs = rhs.get_unchecked((i, j)).inlined_clone();
*e = f(e.inlined_clone(), rhs)
}
}
}
}
#[inline]
pub fn zip_zip_apply<N2, R2, C2, S2, N3, R3, C3, S3>(
&mut self,
b: &Matrix<N2, R2, C2, S2>,
c: &Matrix<N3, R3, C3, S3>,
mut f: impl FnMut(N, N2, N3) -> N,
) where
S: StorageMut<N, R, C>,
N2: Scalar,
R2: Dim,
C2: Dim,
S2: Storage<N2, R2, C2>,
N3: Scalar,
R3: Dim,
C3: Dim,
S3: Storage<N3, R3, C3>,
ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>,
ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>,
{
let (nrows, ncols) = self.shape();
assert_eq!(
(nrows, ncols),
b.shape(),
"Matrix simultaneous traversal error: dimension mismatch."
);
assert_eq!(
(nrows, ncols),
c.shape(),
"Matrix simultaneous traversal error: dimension mismatch."
);
for j in 0..ncols {
for i in 0..nrows {
unsafe {
let e = self.data.get_unchecked_mut(i, j);
let b = b.get_unchecked((i, j)).inlined_clone();
let c = c.get_unchecked((i, j)).inlined_clone();
*e = f(e.inlined_clone(), b, c)
}
}
}
}
}
impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
#[inline]
pub fn iter(&self) -> MatrixIter<N, R, C, S> {
MatrixIter::new(&self.data)
}
#[inline]
pub fn row_iter(&self) -> RowIter<N, R, C, S> {
RowIter::new(self)
}
#[inline]
pub fn column_iter(&self) -> ColumnIter<N, R, C, S> {
ColumnIter::new(self)
}
#[inline]
pub fn iter_mut(&mut self) -> MatrixIterMut<N, R, C, S>
where
S: StorageMut<N, R, C>,
{
MatrixIterMut::new(&mut self.data)
}
#[inline]
pub fn row_iter_mut(&mut self) -> RowIterMut<N, R, C, S>
where
S: StorageMut<N, R, C>,
{
RowIterMut::new(self)
}
#[inline]
pub fn column_iter_mut(&mut self) -> ColumnIterMut<N, R, C, S>
where
S: StorageMut<N, R, C>,
{
ColumnIterMut::new(self)
}
}
impl<N: Scalar, R: Dim, C: Dim, S: StorageMut<N, R, C>> Matrix<N, R, C, S> {
#[inline]
pub fn as_mut_ptr(&mut self) -> *mut N {
self.data.ptr_mut()
}
#[inline]
pub unsafe fn swap_unchecked(&mut self, row_cols1: (usize, usize), row_cols2: (usize, usize)) {
debug_assert!(row_cols1.0 < self.nrows() && row_cols1.1 < self.ncols());
debug_assert!(row_cols2.0 < self.nrows() && row_cols2.1 < self.ncols());
self.data.swap_unchecked(row_cols1, row_cols2)
}
#[inline]
pub fn swap(&mut self, row_cols1: (usize, usize), row_cols2: (usize, usize)) {
let (nrows, ncols) = self.shape();
assert!(
row_cols1.0 < nrows && row_cols1.1 < ncols,
"Matrix elements swap index out of bounds."
);
assert!(
row_cols2.0 < nrows && row_cols2.1 < ncols,
"Matrix elements swap index out of bounds."
);
unsafe { self.swap_unchecked(row_cols1, row_cols2) }
}
#[inline]
pub fn copy_from_slice(&mut self, slice: &[N]) {
let (nrows, ncols) = self.shape();
assert!(
nrows * ncols == slice.len(),
"The slice must contain the same number of elements as the matrix."
);
for j in 0..ncols {
for i in 0..nrows {
unsafe {
*self.get_unchecked_mut((i, j)) =
slice.get_unchecked(i + j * nrows).inlined_clone();
}
}
}
}
#[inline]
pub fn copy_from<R2, C2, SB>(&mut self, other: &Matrix<N, R2, C2, SB>)
where
R2: Dim,
C2: Dim,
SB: Storage<N, R2, C2>,
ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>,
{
assert!(
self.shape() == other.shape(),
"Unable to copy from a matrix with a different shape."
);
for j in 0..self.ncols() {
for i in 0..self.nrows() {
unsafe {
*self.get_unchecked_mut((i, j)) = other.get_unchecked((i, j)).inlined_clone();
}
}
}
}
#[inline]
pub fn tr_copy_from<R2, C2, SB>(&mut self, other: &Matrix<N, R2, C2, SB>)
where
R2: Dim,
C2: Dim,
SB: Storage<N, R2, C2>,
ShapeConstraint: DimEq<R, C2> + SameNumberOfColumns<C, R2>,
{
let (nrows, ncols) = self.shape();
assert!(
(ncols, nrows) == other.shape(),
"Unable to copy from a matrix with incompatible shape."
);
for j in 0..ncols {
for i in 0..nrows {
unsafe {
*self.get_unchecked_mut((i, j)) = other.get_unchecked((j, i)).inlined_clone();
}
}
}
}
#[inline]
pub fn apply_into<F: FnMut(N) -> N>(mut self, f: F) -> Self {
self.apply(f);
self
}
}
impl<N: Scalar, D: Dim, S: Storage<N, D>> Vector<N, D, S> {
#[inline]
pub unsafe fn vget_unchecked(&self, i: usize) -> &N {
debug_assert!(i < self.nrows(), "Vector index out of bounds.");
let i = i * self.strides().0;
self.data.get_unchecked_linear(i)
}
}
impl<N: Scalar, D: Dim, S: StorageMut<N, D>> Vector<N, D, S> {
#[inline]
pub unsafe fn vget_unchecked_mut(&mut self, i: usize) -> &mut N {
debug_assert!(i < self.nrows(), "Vector index out of bounds.");
let i = i * self.strides().0;
self.data.get_unchecked_linear_mut(i)
}
}
impl<N: Scalar, R: Dim, C: Dim, S: ContiguousStorage<N, R, C>> Matrix<N, R, C, S> {
#[inline]
pub fn as_slice(&self) -> &[N] {
self.data.as_slice()
}
}
impl<N: Scalar, R: Dim, C: Dim, S: ContiguousStorageMut<N, R, C>> Matrix<N, R, C, S> {
#[inline]
pub fn as_mut_slice(&mut self) -> &mut [N] {
self.data.as_mut_slice()
}
}
impl<N: Scalar, D: Dim, S: StorageMut<N, D, D>> Matrix<N, D, D, S> {
pub fn transpose_mut(&mut self) {
assert!(
self.is_square(),
"Unable to transpose a non-square matrix in-place."
);
let dim = self.shape().0;
for i in 1..dim {
for j in 0..i {
unsafe { self.swap_unchecked((i, j), (j, i)) }
}
}
}
}
impl<N: SimdComplexField, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
#[inline]
pub fn adjoint_to<R2, C2, SB>(&self, out: &mut Matrix<N, R2, C2, SB>)
where
R2: Dim,
C2: Dim,
SB: StorageMut<N, R2, C2>,
ShapeConstraint: SameNumberOfRows<R, C2> + SameNumberOfColumns<C, R2>,
{
let (nrows, ncols) = self.shape();
assert!(
(ncols, nrows) == out.shape(),
"Incompatible shape for transpose-copy."
);
for i in 0..nrows {
for j in 0..ncols {
unsafe {
*out.get_unchecked_mut((j, i)) = self.get_unchecked((i, j)).simd_conjugate();
}
}
}
}
#[inline]
#[must_use = "Did you mean to use adjoint_mut()?"]
pub fn adjoint(&self) -> MatrixMN<N, C, R>
where
DefaultAllocator: Allocator<N, C, R>,
{
let (nrows, ncols) = self.data.shape();
unsafe {
let mut res: MatrixMN<_, C, R> = Matrix::new_uninitialized_generic(ncols, nrows);
self.adjoint_to(&mut res);
res
}
}
#[deprecated(note = "Renamed `self.adjoint_to(out)`.")]
#[inline]
pub fn conjugate_transpose_to<R2, C2, SB>(&self, out: &mut Matrix<N, R2, C2, SB>)
where
R2: Dim,
C2: Dim,
SB: StorageMut<N, R2, C2>,
ShapeConstraint: SameNumberOfRows<R, C2> + SameNumberOfColumns<C, R2>,
{
self.adjoint_to(out)
}
#[deprecated(note = "Renamed `self.adjoint()`.")]
#[inline]
pub fn conjugate_transpose(&self) -> MatrixMN<N, C, R>
where
DefaultAllocator: Allocator<N, C, R>,
{
self.adjoint()
}
#[inline]
#[must_use = "Did you mean to use conjugate_mut()?"]
pub fn conjugate(&self) -> MatrixMN<N, R, C>
where
DefaultAllocator: Allocator<N, R, C>,
{
self.map(|e| e.simd_conjugate())
}
#[inline]
#[must_use = "Did you mean to use unscale_mut()?"]
pub fn unscale(&self, real: N::SimdRealField) -> MatrixMN<N, R, C>
where
DefaultAllocator: Allocator<N, R, C>,
{
self.map(|e| e.simd_unscale(real))
}
#[inline]
#[must_use = "Did you mean to use scale_mut()?"]
pub fn scale(&self, real: N::SimdRealField) -> MatrixMN<N, R, C>
where
DefaultAllocator: Allocator<N, R, C>,
{
self.map(|e| e.simd_scale(real))
}
}
impl<N: SimdComplexField, R: Dim, C: Dim, S: StorageMut<N, R, C>> Matrix<N, R, C, S> {
#[inline]
pub fn conjugate_mut(&mut self) {
self.apply(|e| e.simd_conjugate())
}
#[inline]
pub fn unscale_mut(&mut self, real: N::SimdRealField) {
self.apply(|e| e.simd_unscale(real))
}
#[inline]
pub fn scale_mut(&mut self, real: N::SimdRealField) {
self.apply(|e| e.simd_scale(real))
}
}
impl<N: SimdComplexField, D: Dim, S: StorageMut<N, D, D>> Matrix<N, D, D, S> {
#[deprecated(note = "Renamed to `self.adjoint_mut()`.")]
pub fn conjugate_transform_mut(&mut self) {
self.adjoint_mut()
}
pub fn adjoint_mut(&mut self) {
assert!(
self.is_square(),
"Unable to transpose a non-square matrix in-place."
);
let dim = self.shape().0;
for i in 0..dim {
for j in 0..i {
unsafe {
let ref_ij = self.get_unchecked_mut((i, j)) as *mut N;
let ref_ji = self.get_unchecked_mut((j, i)) as *mut N;
let conj_ij = (*ref_ij).simd_conjugate();
let conj_ji = (*ref_ji).simd_conjugate();
*ref_ij = conj_ji;
*ref_ji = conj_ij;
}
}
{
let diag = unsafe { self.get_unchecked_mut((i, i)) };
*diag = diag.simd_conjugate();
}
}
}
}
impl<N: Scalar, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S> {
#[inline]
pub fn diagonal(&self) -> VectorN<N, D>
where
DefaultAllocator: Allocator<N, D>,
{
self.map_diagonal(|e| e)
}
pub fn map_diagonal<N2: Scalar>(&self, mut f: impl FnMut(N) -> N2) -> VectorN<N2, D>
where
DefaultAllocator: Allocator<N2, D>,
{
assert!(
self.is_square(),
"Unable to get the diagonal of a non-square matrix."
);
let dim = self.data.shape().0;
let mut res = unsafe { VectorN::new_uninitialized_generic(dim, U1) };
for i in 0..dim.value() {
unsafe {
*res.vget_unchecked_mut(i) = f(self.get_unchecked((i, i)).inlined_clone());
}
}
res
}
#[inline]
pub fn trace(&self) -> N
where
N: Scalar + Zero + ClosedAdd,
{
assert!(
self.is_square(),
"Cannot compute the trace of non-square matrix."
);
let dim = self.data.shape().0;
let mut res = N::zero();
for i in 0..dim.value() {
res += unsafe { self.get_unchecked((i, i)).inlined_clone() };
}
res
}
}
impl<N: SimdComplexField, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S> {
#[inline]
pub fn symmetric_part(&self) -> MatrixMN<N, D, D>
where
DefaultAllocator: Allocator<N, D, D>,
{
assert!(
self.is_square(),
"Cannot compute the symmetric part of a non-square matrix."
);
let mut tr = self.transpose();
tr += self;
tr *= crate::convert::<_, N>(0.5);
tr
}
#[inline]
pub fn hermitian_part(&self) -> MatrixMN<N, D, D>
where
DefaultAllocator: Allocator<N, D, D>,
{
assert!(
self.is_square(),
"Cannot compute the hermitian part of a non-square matrix."
);
let mut tr = self.adjoint();
tr += self;
tr *= crate::convert::<_, N>(0.5);
tr
}
}
impl<N: Scalar + Zero + One, D: DimAdd<U1> + IsNotStaticOne, S: Storage<N, D, D>>
Matrix<N, D, D, S>
{
#[inline]
pub fn to_homogeneous(&self) -> MatrixN<N, DimSum<D, U1>>
where
DefaultAllocator: Allocator<N, DimSum<D, U1>, DimSum<D, U1>>,
{
assert!(
self.is_square(),
"Only square matrices can currently be transformed to homogeneous coordinates."
);
let dim = DimSum::<D, U1>::from_usize(self.nrows() + 1);
let mut res = MatrixN::identity_generic(dim, dim);
res.generic_slice_mut::<D, D>((0, 0), self.data.shape())
.copy_from(&self);
res
}
}
impl<N: Scalar + Zero, D: DimAdd<U1>, S: Storage<N, D>> Vector<N, D, S> {
#[inline]
pub fn to_homogeneous(&self) -> VectorN<N, DimSum<D, U1>>
where
DefaultAllocator: Allocator<N, DimSum<D, U1>>,
{
self.push(N::zero())
}
#[inline]
pub fn from_homogeneous<SB>(v: Vector<N, DimSum<D, U1>, SB>) -> Option<VectorN<N, D>>
where
SB: Storage<N, DimSum<D, U1>>,
DefaultAllocator: Allocator<N, D>,
{
if v[v.len() - 1].is_zero() {
let nrows = D::from_usize(v.len() - 1);
Some(v.generic_slice((0, 0), (nrows, U1)).into_owned())
} else {
None
}
}
}
impl<N: Scalar + Zero, D: DimAdd<U1>, S: Storage<N, D>> Vector<N, D, S> {
#[inline]
pub fn push(&self, element: N) -> VectorN<N, DimSum<D, U1>>
where
DefaultAllocator: Allocator<N, DimSum<D, U1>>,
{
let len = self.len();
let hnrows = DimSum::<D, U1>::from_usize(len + 1);
let mut res = unsafe { VectorN::<N, _>::new_uninitialized_generic(hnrows, U1) };
res.generic_slice_mut((0, 0), self.data.shape())
.copy_from(self);
res[(len, 0)] = element;
res
}
}
impl<N, R: Dim, C: Dim, S> AbsDiffEq for Matrix<N, R, C, S>
where
N: Scalar + AbsDiffEq,
S: Storage<N, R, C>,
N::Epsilon: Copy,
{
type Epsilon = N::Epsilon;
#[inline]
fn default_epsilon() -> Self::Epsilon {
N::default_epsilon()
}
#[inline]
fn abs_diff_eq(&self, other: &Self, epsilon: Self::Epsilon) -> bool {
self.iter()
.zip(other.iter())
.all(|(a, b)| a.abs_diff_eq(b, epsilon))
}
}
impl<N, R: Dim, C: Dim, S> RelativeEq for Matrix<N, R, C, S>
where
N: Scalar + RelativeEq,
S: Storage<N, R, C>,
N::Epsilon: Copy,
{
#[inline]
fn default_max_relative() -> Self::Epsilon {
N::default_max_relative()
}
#[inline]
fn relative_eq(
&self,
other: &Self,
epsilon: Self::Epsilon,
max_relative: Self::Epsilon,
) -> bool {
self.relative_eq(other, epsilon, max_relative)
}
}
impl<N, R: Dim, C: Dim, S> UlpsEq for Matrix<N, R, C, S>
where
N: Scalar + UlpsEq,
S: Storage<N, R, C>,
N::Epsilon: Copy,
{
#[inline]
fn default_max_ulps() -> u32 {
N::default_max_ulps()
}
#[inline]
fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool {
assert!(self.shape() == other.shape());
self.iter()
.zip(other.iter())
.all(|(a, b)| a.ulps_eq(b, epsilon, max_ulps))
}
}
impl<N, R: Dim, C: Dim, S> PartialOrd for Matrix<N, R, C, S>
where
N: Scalar + PartialOrd,
S: Storage<N, R, C>,
{
#[inline]
fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
if self.shape() != other.shape() {
return None;
}
if self.nrows() == 0 || self.ncols() == 0 {
return Some(Ordering::Equal);
}
let mut first_ord = unsafe {
self.data
.get_unchecked_linear(0)
.partial_cmp(other.data.get_unchecked_linear(0))
};
if let Some(first_ord) = first_ord.as_mut() {
let mut it = self.iter().zip(other.iter());
let _ = it.next();
for (left, right) in it {
if let Some(ord) = left.partial_cmp(right) {
match ord {
Ordering::Equal => { }
Ordering::Less => {
if *first_ord == Ordering::Greater {
return None;
}
*first_ord = ord
}
Ordering::Greater => {
if *first_ord == Ordering::Less {
return None;
}
*first_ord = ord
}
}
} else {
return None;
}
}
}
first_ord
}
#[inline]
fn lt(&self, right: &Self) -> bool {
assert_eq!(
self.shape(),
right.shape(),
"Matrix comparison error: dimensions mismatch."
);
self.iter().zip(right.iter()).all(|(a, b)| a.lt(b))
}
#[inline]
fn le(&self, right: &Self) -> bool {
assert_eq!(
self.shape(),
right.shape(),
"Matrix comparison error: dimensions mismatch."
);
self.iter().zip(right.iter()).all(|(a, b)| a.le(b))
}
#[inline]
fn gt(&self, right: &Self) -> bool {
assert_eq!(
self.shape(),
right.shape(),
"Matrix comparison error: dimensions mismatch."
);
self.iter().zip(right.iter()).all(|(a, b)| a.gt(b))
}
#[inline]
fn ge(&self, right: &Self) -> bool {
assert_eq!(
self.shape(),
right.shape(),
"Matrix comparison error: dimensions mismatch."
);
self.iter().zip(right.iter()).all(|(a, b)| a.ge(b))
}
}
impl<N, R: Dim, C: Dim, S> Eq for Matrix<N, R, C, S>
where
N: Scalar + Eq,
S: Storage<N, R, C>,
{
}
impl<N, R, R2, C, C2, S, S2> PartialEq<Matrix<N, R2, C2, S2>> for Matrix<N, R, C, S>
where
N: Scalar + PartialEq,
C: Dim,
C2: Dim,
R: Dim,
R2: Dim,
S: Storage<N, R, C>,
S2: Storage<N, R2, C2>,
{
#[inline]
fn eq(&self, right: &Matrix<N, R2, C2, S2>) -> bool {
self.shape() == right.shape() && self.iter().zip(right.iter()).all(|(l, r)| l == r)
}
}
macro_rules! impl_fmt {
($trait: path, $fmt_str_without_precision: expr, $fmt_str_with_precision: expr) => {
impl<N, R: Dim, C: Dim, S> $trait for Matrix<N, R, C, S>
where
N: Scalar + $trait,
S: Storage<N, R, C>,
DefaultAllocator: Allocator<usize, R, C>,
{
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
#[cfg(feature = "std")]
fn val_width<N: Scalar + $trait>(val: &N, f: &mut fmt::Formatter) -> usize {
match f.precision() {
Some(precision) => format!($fmt_str_with_precision, val, precision)
.chars()
.count(),
None => format!($fmt_str_without_precision, val).chars().count(),
}
}
#[cfg(not(feature = "std"))]
fn val_width<N: Scalar + $trait>(_: &N, _: &mut fmt::Formatter) -> usize {
4
}
let (nrows, ncols) = self.data.shape();
if nrows.value() == 0 || ncols.value() == 0 {
return write!(f, "[ ]");
}
let mut max_length = 0;
let mut lengths: MatrixMN<usize, R, C> = Matrix::zeros_generic(nrows, ncols);
let (nrows, ncols) = self.shape();
for i in 0..nrows {
for j in 0..ncols {
lengths[(i, j)] = val_width(&self[(i, j)], f);
max_length = crate::max(max_length, lengths[(i, j)]);
}
}
let max_length_with_space = max_length + 1;
writeln!(f)?;
writeln!(
f,
" ┌ {:>width$} ┐",
"",
width = max_length_with_space * ncols - 1
)?;
for i in 0..nrows {
write!(f, " │")?;
for j in 0..ncols {
let number_length = lengths[(i, j)] + 1;
let pad = max_length_with_space - number_length;
write!(f, " {:>thepad$}", "", thepad = pad)?;
match f.precision() {
Some(precision) => {
write!(f, $fmt_str_with_precision, (*self)[(i, j)], precision)?
}
None => write!(f, $fmt_str_without_precision, (*self)[(i, j)])?,
}
}
writeln!(f, " │")?;
}
writeln!(
f,
" └ {:>width$} ┘",
"",
width = max_length_with_space * ncols - 1
)?;
writeln!(f)
}
}
};
}
impl_fmt!(fmt::Display, "{}", "{:.1$}");
impl_fmt!(fmt::LowerExp, "{:e}", "{:.1$e}");
impl_fmt!(fmt::UpperExp, "{:E}", "{:.1$E}");
impl_fmt!(fmt::Octal, "{:o}", "{:1$o}");
impl_fmt!(fmt::LowerHex, "{:x}", "{:1$x}");
impl_fmt!(fmt::UpperHex, "{:X}", "{:1$X}");
impl_fmt!(fmt::Binary, "{:b}", "{:.1$b}");
impl_fmt!(fmt::Pointer, "{:p}", "{:.1$p}");
#[test]
fn lower_exp() {
let test = crate::Matrix2::new(1e6, 2e5, 2e-5, 1.);
assert_eq!(
format!("{:e}", test),
r"
┌ ┐
│ 1e6 2e5 │
│ 2e-5 1e0 │
└ ┘
"
)
}
impl<N: Scalar + ClosedAdd + ClosedSub + ClosedMul, R: Dim, C: Dim, S: Storage<N, R, C>>
Matrix<N, R, C, S>
{
#[inline]
pub fn perp<R2, C2, SB>(&self, b: &Matrix<N, R2, C2, SB>) -> N
where
R2: Dim,
C2: Dim,
SB: Storage<N, R2, C2>,
ShapeConstraint: SameNumberOfRows<R, U2>
+ SameNumberOfColumns<C, U1>
+ SameNumberOfRows<R2, U2>
+ SameNumberOfColumns<C2, U1>,
{
assert!(
self.shape() == (2, 1),
"2D perpendicular product requires (2, 1) vector but found {:?}",
self.shape()
);
unsafe {
self.get_unchecked((0, 0)).inlined_clone() * b.get_unchecked((1, 0)).inlined_clone()
- self.get_unchecked((1, 0)).inlined_clone()
* b.get_unchecked((0, 0)).inlined_clone()
}
}
#[inline]
pub fn cross<R2, C2, SB>(&self, b: &Matrix<N, R2, C2, SB>) -> MatrixCross<N, R, C, R2, C2>
where
R2: Dim,
C2: Dim,
SB: Storage<N, R2, C2>,
DefaultAllocator: SameShapeAllocator<N, R, C, R2, C2>,
ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>,
{
let shape = self.shape();
assert_eq!(shape, b.shape(), "Vector cross product dimension mismatch.");
assert!(
(shape.0 == 3 && shape.1 == 1) || (shape.0 == 1 && shape.1 == 3),
"Vector cross product dimension mismatch: must be (3, 1) or (1, 3) but found {:?}.",
shape
);
if shape.0 == 3 {
unsafe {
let nrows = SameShapeR::<R, R2>::from_usize(3);
let ncols = SameShapeC::<C, C2>::from_usize(1);
let mut res = Matrix::new_uninitialized_generic(nrows, ncols);
let ax = self.get_unchecked((0, 0));
let ay = self.get_unchecked((1, 0));
let az = self.get_unchecked((2, 0));
let bx = b.get_unchecked((0, 0));
let by = b.get_unchecked((1, 0));
let bz = b.get_unchecked((2, 0));
*res.get_unchecked_mut((0, 0)) = ay.inlined_clone() * bz.inlined_clone()
- az.inlined_clone() * by.inlined_clone();
*res.get_unchecked_mut((1, 0)) = az.inlined_clone() * bx.inlined_clone()
- ax.inlined_clone() * bz.inlined_clone();
*res.get_unchecked_mut((2, 0)) = ax.inlined_clone() * by.inlined_clone()
- ay.inlined_clone() * bx.inlined_clone();
res
}
} else {
unsafe {
let nrows = SameShapeR::<R, R2>::from_usize(1);
let ncols = SameShapeC::<C, C2>::from_usize(3);
let mut res = Matrix::new_uninitialized_generic(nrows, ncols);
let ax = self.get_unchecked((0, 0));
let ay = self.get_unchecked((0, 1));
let az = self.get_unchecked((0, 2));
let bx = b.get_unchecked((0, 0));
let by = b.get_unchecked((0, 1));
let bz = b.get_unchecked((0, 2));
*res.get_unchecked_mut((0, 0)) = ay.inlined_clone() * bz.inlined_clone()
- az.inlined_clone() * by.inlined_clone();
*res.get_unchecked_mut((0, 1)) = az.inlined_clone() * bx.inlined_clone()
- ax.inlined_clone() * bz.inlined_clone();
*res.get_unchecked_mut((0, 2)) = ax.inlined_clone() * by.inlined_clone()
- ay.inlined_clone() * bx.inlined_clone();
res
}
}
}
}
impl<N: Scalar + Field, S: Storage<N, U3>> Vector<N, U3, S>
where
DefaultAllocator: Allocator<N, U3>,
{
#[inline]
pub fn cross_matrix(&self) -> MatrixN<N, U3> {
MatrixN::<N, U3>::new(
N::zero(),
-self[2].inlined_clone(),
self[1].inlined_clone(),
self[2].inlined_clone(),
N::zero(),
-self[0].inlined_clone(),
-self[1].inlined_clone(),
self[0].inlined_clone(),
N::zero(),
)
}
}
impl<N: SimdComplexField, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
#[inline]
pub fn angle<R2: Dim, C2: Dim, SB>(&self, other: &Matrix<N, R2, C2, SB>) -> N::SimdRealField
where
SB: Storage<N, R2, C2>,
ShapeConstraint: DimEq<R, R2> + DimEq<C, C2>,
{
let prod = self.dotc(other);
let n1 = self.norm();
let n2 = other.norm();
if n1.is_zero() || n2.is_zero() {
N::SimdRealField::zero()
} else {
let cang = prod.simd_real() / (n1 * n2);
cang.simd_clamp(-N::SimdRealField::one(), N::SimdRealField::one())
.simd_acos()
}
}
}
impl<N, R: Dim, C: Dim, S> AbsDiffEq for Unit<Matrix<N, R, C, S>>
where
N: Scalar + AbsDiffEq,
S: Storage<N, R, C>,
N::Epsilon: Copy,
{
type Epsilon = N::Epsilon;
#[inline]
fn default_epsilon() -> Self::Epsilon {
N::default_epsilon()
}
#[inline]
fn abs_diff_eq(&self, other: &Self, epsilon: Self::Epsilon) -> bool {
self.as_ref().abs_diff_eq(other.as_ref(), epsilon)
}
}
impl<N, R: Dim, C: Dim, S> RelativeEq for Unit<Matrix<N, R, C, S>>
where
N: Scalar + RelativeEq,
S: Storage<N, R, C>,
N::Epsilon: Copy,
{
#[inline]
fn default_max_relative() -> Self::Epsilon {
N::default_max_relative()
}
#[inline]
fn relative_eq(
&self,
other: &Self,
epsilon: Self::Epsilon,
max_relative: Self::Epsilon,
) -> bool {
self.as_ref()
.relative_eq(other.as_ref(), epsilon, max_relative)
}
}
impl<N, R: Dim, C: Dim, S> UlpsEq for Unit<Matrix<N, R, C, S>>
where
N: Scalar + UlpsEq,
S: Storage<N, R, C>,
N::Epsilon: Copy,
{
#[inline]
fn default_max_ulps() -> u32 {
N::default_max_ulps()
}
#[inline]
fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool {
self.as_ref().ulps_eq(other.as_ref(), epsilon, max_ulps)
}
}
impl<N, R, C, S> Hash for Matrix<N, R, C, S>
where
N: Scalar + Hash,
R: Dim,
C: Dim,
S: Storage<N, R, C>,
{
fn hash<H: Hasher>(&self, state: &mut H) {
let (nrows, ncols) = self.shape();
(nrows, ncols).hash(state);
for j in 0..ncols {
for i in 0..nrows {
unsafe {
self.get_unchecked((i, j)).hash(state);
}
}
}
}
}